A modified theta projection model for creep behavior of RPV steel 16MND5

Journal of Materials Science & Technology 47 (2020) 231–242

Contents lists available at ScienceDirect

Journal of Materials Science & Technology journal homepage: www.jmst.org

Research Article

A modified theta projection model for creep behavior of RPV steel 16MND5 Peng Yu ∗ , Weimin Ma Division of Nuclear Power Safety, Royal Institute of Technology (KTH), Stockholm 106 91, Sweden

a r t i c l e

i n f o

Article history: Received 8 August 2019 Received in revised form 22 November 2019 Accepted 2 December 2019 Available online 21 February 2020 Keywords: 16MND5 steel Creep modelling Tertiary stage Reactor pressure vessel Theta projection model

a b s t r a c t During a hypothetical severe accident of light water reactors, the reactor pressure vessel (RPV) could fail due to its creep under the influence of high-temperature corium. Hence, modelling of creep behavior of the RPV is paramount to reactor safety analysis since it predicts the transition point of accident progression from in-vessel to ex-vessel phase. In the present study we proposed a new creep model for the classical French RPV steel 16MND5, which is adapted from the “theta-projection model” and contains all three stages of a creep process. Creep curves are expressed as a function of time with five model parameters i (i = 1 − 4 and m). A model parameter dataset was constructed by fitting experimental creep curves into this function. To correlate the creep curves for different temperatures and stress loads, we directly interpolate the model’s parameters i (i = 1 − 4 and m) from this dataset, in contrast to the conventional “theta-projection model” which employs an extra single correlation for each i (i = 1 − 4 and m), to better accommodate all experimental curves over the wide ranges of temperature and stress loads. We also put a constraint on the trend of the creep strain that it would monotonically increase with temperature and stress load. A good agreement was achieved between each experimental creep curve and corresponding model’s prediction. The widely used time-hardening and strain-hardening models were performing reasonably well in the new method. © 2020 Published by Elsevier Ltd on behalf of The editorial office of Journal of Materials Science & Technology.

1. Introduction Creep is an irreversible time-dependent inelastic deformation process which happens as long as there is a stress. The molecular mechanism of creep can be generally classified as diffusion and dislocation, the former of which normally occurs due to the presence of vacancies within the crystal lattice (or the grain) while the latter occurs if there is a slip by the passage of crystal dislocations [1]. There are two types of diffusion creep [2], called Coble creep (favored at lower temperatures) if the diffusion paths are predominantly through the grain boundaries, and Nabarro-Herring creep (favored at higher temperatures) if through the grains themselves. This mechanism of creep tends to dominate at high stresses and relatively low temperatures. Dislocations can move by gliding in a slip plane, a process requiring little thermal activation. However, the rate-determining step for their motion is often a climb process, which requires diffusion and is thus time-dependent and favored by higher temperatures. The creep behavior is highly dependent

∗ Corresponding author. E-mail address: [email protected] (P. Yu).

on temperature and stress: creep can be negligibly small under low temperature and stress; it occurs faster at higher temperature and stress. Generally, it is expected that creep rate would have an exponential dependence on temperature for steady-state creep as shown in Eq. (1) and a power dependence on stress as expressed in Eq. (2): ε˙ ∝ exp(−Q/RT ε˙ ∝ 

n

(1) (2)

where ε˙ is the creep strain rate, Q is the activation energy associated with the rate, R is the universal gas constant (= 8.31 J/mol/K), T is the temperature in K, and  is the stress. A creep curve describes the creep strain development with time under constant stress and temperature as shown in Fig. 1. Typical creep curves can be characterized by three stages: primary stage where the creep strain rate decreases with time, secondary stage where the creep strain rate almost keeps constant and tertiary stage where the creep strain rate drastically increases with time. In the engineering modelling of kinematics for deformation and motion, the creep effect is considered by introducing extra terms into the conservation equations. To close the equation sets,

https://doi.org/10.1016/j.jmst.2020.02.016 1005-0302/© 2020 Published by Elsevier Ltd on behalf of The editorial office of Journal of Materials Science & Technology.

232

P. Yu, W. Ma / Journal of Materials Science & Technology 47 (2020) 231–242

Fig. 1. Creep curve under constant temperature and stress.

a creep model (a constitutive equation describing the creep strain or creep strain rate as functions of stress, temperature, time and creep strain.) is required. There are many creep models available for creep simulations, ranging from simple-phenomenological to complex-constitutive ones. Commonly used models include the typical time hardening and strain hardening models for primary stage of creep; Norton’ Law [3] and Exponential form for the secondary stage of creep; Norton-Bailey, RCC-MR [4] and Bartsch [5] for both the primary and secondary stages; Theta Projection [1], Rabotnov-Kachanov [6], Baker-Cane [7] and Dyson-McLean [8] for all three stages. The commercial software ANSYS [9] and ABAQUS [10] have default standard creep models that only include the primary and/or the secondary stage(s). The two codes offer the interfaces to implement user-defined creep laws through either ANSYS User Program Functions or ABAQUSs user subroutines. One reason why the tertiary stage is not modeled by default might be that for most industrial applications the failure analysis is beyond interest, and therefore the built-in creep laws in the commercial software is supposed to be sufficient [11]. The extra user-defined creep law should be implemented when a user is interested in the tertiary creep stage and/or wants more accurate creep modelling. The steel 16MND5 is a typical material widely used to construct the reactor pressure vessel (RPV). 16MND5 is a low-alloyed ferritic steel whose chemical composition is given in Table 1 [12,13]. It has undergone several heat treatments during its elaboration: two austenizations followed by water quench, a tempering and a stress-relief treatment which then lead to a ‘tempered bainite’ structure. Since the RPV also serves as the second physical barrier to contain the core material and prevent radioactivity release under operations and accidents, its integrity is of great importance to nuclear power safety. A severe accident with core meltdown may threaten the RPV integrity, since the relocated reactor core melt (corium) in the RPV lower head will exert thermal and mechanical loads on the vessel wall, which may ultimately lead to vessel failure due to creep under the high temperature and significant stress. The “vessel failure due to creep” means that the lower head wall does not only go through the primary and secondary stages, but also reaches the tertiary stage of creep. Hence, accurate modelling of creep in three stages is important to RPV integrity assessment during a hypothetical severe accident. For the creep behavior of the 16MND5 steel, the REVISA Experiment [14] at CEA performed the tensile and creep tests using 16MND5 specimens at a wide range of temperature from 600 ◦ C to 1300 ◦ C. Four tests with different stress loads were carried out for each temperature. It provided the most important experimental data of the creep behavior of 16MND5. Several treatments have been explored to correlate the experimental data into numerical analysis of structural integrity. Iknoen [15] fitted the data of the REVISA experiment into a three-stage creep model: three creep stages are distinguished by separate correlations, and creep at

each stage is described by a time-hardening model with different temperature-dependent coefficients. Altstadt and Moessner [16] transferred the experimental curves to a discretized database: under each given temperature T and stress ␴, the creep curve was ˙ indicating under such therdiscretized into finite point pairs (␧, ε) mal and mechanical loads the creep strain rate ε˙ is subjected to creep strain ␧. The creep strain rate for any other case can be interpolated over known temperature, stress, and strain from the database. Due to the limited information, we do not know the performance of such approach. Willschutz et al. [17] adopted this database idea and used in their creep failure simulation. A different treatment [18] was later proposed to predict the creep strain acceleration in the tertiary stage. A correction factor 1/(1 − D) was introduced to the creep strain increment, where D is a damage factor (0 < D < 1). (The usage of such damage variables was initially proposed by Kachanov [19] and Rabotnov [20]. A related review can also be found in [21].) In this research, D was defined by both creep and plasticity information. As damage parameter grows close to 1, the creep strain rate would be enlarged more significantly as is multiplied by this factor 1/(1 − D). The comparison [22] between model predictions and experimental data shows that, the model can predict three-stage creep behavior and generally agrees with the experimental data. However, large discrepancies exist in a good number of creep curves: regarding the time reaching the same strain, 8 out of 32 curves show a discrepancy of ∼100% between prediction and experiment; 5 more curves show a discrepancy of ∼25%. Moreover, it is noted this damage factor is dependent on plastic performance, and therefore a change in plastic model might numerically have influence on the creep prediction. As the creep model plays a pivotal role in modelling vessel failure during a severe accident, there is a clear need to improve its fidelity for a more credible RPV integrity assessment. The present study is intended to fill in such a gap by developing a new creep model for the 16MND5 steel which will not only contain all three creep stages, but also fits into all experimental data with minimum discrepancy. The new model was adapted from the idea of the theta projection method. Model modifications were introduced to accommodate the wide range of temperature and stress in experimental data sets. We also considered monotonic increase of creep strain with temperature and stress as a consider for developing the model. The results indicated that a good agreement was achieved between the model’s prediction and experimental curve for each pair of temperature and stress load. The widely used time-hardening and strain-hardening models were performing reasonably well. 2. Modelling 2.1. Classical theta projection method and several modified theta projection methods Theta projection method was first proposed by Evans et al. [23] to accurately predict long-term creep and creep-fracture properties. The general idea is to represent the creep strain curve under constant temperature and stress as follows ε = 1 (1 − e−2 t ) + 3 (e4 t − 1

(3)

where ε is the creep strain, t is time and 1 -4 are non-negative model parameters. As illustrated in Fig. 2, the first term on right hand side increases slowly with time, and the slope decreases gradually to zero, which is suitable to represent the creep strain change in primary stage. The second term is negligibly small at the early stage and increases exponentially with time at the later stage, which would be suitable to represent the tertiary stage. In the transition of dominancy from the first term to the second term,

P. Yu, W. Ma / Journal of Materials Science & Technology 47 (2020) 231–242

233

Table 1 Chemical composition of 16MND5 steel (weight-%). C

S

P

Mn

Si

Ni

Cr

Mo

V

Cu

Co

Al

N

O

Sn

As

0.159

0.008

0.005

1.37

0.24

0.7

0.17

0.5

<0.01

0.06

<0.01

0.023

0.07

35–36 ppm

50 ppm

180 ppm

m = A n e−Q/RT

(6)

where A and n are coefficients to be determined, Q is the activation energy associated with the rate parameter, R is the universal gas constant (8.31 Jmol-1 K-1). Due to very limited descriptions of the model, its performance can not be assessed. This model was suggested to be used for low/high alloy ferritic, austenitic steels, Al-alloys, Al-matrix composites. Day and Gordon [25] used Eq. (7) to describe the creep of the Nickel-based Super-alloy.



ε = 3 (1 − e + 7

−t/1

) 8 + 4

2 t/2 t+ − 1) (e 5e



8

t 6 (t/5 −1) e 8 5

(7)

8 is called ‘overall scaling’ which is expressed as, Fig. 2. Illustration of theta-projection method (e.g. ␪1 = 0.052, ␪2 = 0.0006, ␪3 = 0.000015, ␪4 = 0.0018).

8 = (1 − e−/0 )

an approximate stable creep strain rate is obtained, which can represent the secondary stage. The influence of temperature T and stress  are reflected in the model parameters i (i = 1, 2, 3, 4). The model assumes the relation between T ,  and  can be expressed as

The creep strains were scaled and compared with scaled experimental data. Overall, it agreed well with experiment data except for some extreme cases showing discrepancy by several times. No further suggestion was given for application to other materials. Kumar et al. [26] included the yield strength and make the model stress ratio dependent to model the metal and alloys. The model expression is shown in Eq. (9). It was applied for two aluminum alloys Al 2124 and Al 7075.

lni = ai + bi  + ci T + di T

(4)

where ai , bi , ci and di (i = 1, 2, 3, 4) are constants. As is explained in [1], this fitting Eq. (4) is suitable when the temperature range of experiment data sets are very narrow such that −1/T is almost proportional to T . The development of this model can be summarized as following two fitting procedures: (1) Fit the experimental creep curves to Eq. (3), after which we can obtain several groups of 1 − 4 for each creep curve under certain temperature and stress. We denote (Tj , j , 1, j , 2, j , 3, j , 4, j ) as parameters for the jth group of experimental data. (2) Fit to Eq. (4) as a function of temperature and stress providing the obtained data (Tj , j , 1, j ). Repeat this fitting procedure for 2 , 3 and 4 . (3) After the two steps we would get 4 groups of (ai , bi , ci , di ) (i = 1, 2, 3, 4) for 1 − 4 which are the parameters finally implemented in model. When using this model, given (Tnew , new ), corresponding 1, new − 4, new can be obtained by plugging Tnew and new to Eq. (4) and then corresponding creep strain curve can be obtained by plugging 1, new − 4, new to Eq. (3). Above method have been implemented to describe creep of pure metals like polycrystalline copper and particle-hardened alloys like 1 Cr 1 Mo 1 V ferritic steel. 2 2

4



εc = 1 1 − e−2 t

(8)

1−  y



+ 3 (e4 t − 1)

(9)

The method gave more accurate prediction for both alloys than the original theta projection model. It was suggested for the application to pure metals and class M alloys (similar creep behavior to pure metals) with intermediate range of stress (40%–70% of the yield stress) and temperature (40%–60% of the melting temperature). Harrison et al. [27] correlated the parameters as power function of the normalized stress /TS . 1 and 3 employ the form of Eq. (10) while 2 and 4 employ the form of Eq. (11). This model was applied for the material ␥-Titanium Aluminide Ti-45Al-2Mn-2Nb. i = Ai i = Ai

  ni

(10)

TS

  ni TS

 exp

−Qi∗

 (11)

RT

where Ai is the scale factors, Qi* the activation energy for i . The predicted strain curves fitted test data well providing well-represented minimum creep rate values. The creep ductility was also successfully predicted. No further suggestion was given for application to other materials. Song et al. [28] expressed the creep strain as







ε = 1 ln 2 t + 1 − 3 ln 1 − 4 t



(12)

Variations of this theta-projection model have been developed by other researchers to better accommodate the experimental data for other different materials, which are all then called a ‘modified theta projection model’ independently. Modified models differ from original model in various ways. The method in [24] introduced an extra linear term in the model which is then read,

The model parameters satisfy Eq. (4). This model was applied to the steel 12Cr1MoV. It showed good flexibility in fitting the entire test data of the creep-resistant steel 12Cr1MoV over the three creep stages. No further suggestion was given for application to other materials. Fu et al. [29] expressed the creep strain as

ε = 1 (1 − e−2 t ) + m t + 3 (e4 t − 1

ε = 1 (1 − e−2 t ) + 3 e4 e

(5)



5 εt



−1

(13)

234

P. Yu, W. Ma / Journal of Materials Science & Technology 47 (2020) 231–242

It was validated with the K465 and DZ125 superalloys. The predicted curves were in good agreement with the experimental data and was much better over the original model as it reduces the prediction uncertainty from 10%–20% to below 5%. This method was suggested for the application to other high temperature structural materials. 2.2. A new modified theta projection method for 16MND5 steel Though we already have the theta-projection method and several modified projection methods, problems still exist when generating a creep model for 16MND5. The experimental data for 16MND5 covers a wide temperature range varying from 600 ◦ C to 1300 ◦ C and a wide stress load range varying from 0.8 MPa to 248 MPa while Eq. (4) in the existing methods is based on a narrow temperature span. Poor parameter fitting results would be expected if we directly fit the parameters 1 − 4 with Eq. (4) as used in the original methods. Therefore, we developed a new modified theta projection model for 16MND5 in which we employed Eq. (5) to express the creep strain under certain strain and temperature. The main difference compared to the original [23] and [24] is that, numerical interpolation will be used for model parameters instead of fitting, i.e., instead of fitting the model parameters as functions of stress and temperature with Eq. (4), the model parameters’ distribution is assumed to be piecewise in the temperature-stress space and they will be interpolated directly from known parameters obtained from experiments. The reason why we did not use Eq. (3) is that the 16MND5 steel shows fairly long secondary stages in many experimental creep curves, which can be better fitted by Eq. (5) than Eq. (3). Researchers (e.g. in [28,30]) sometimes mention the absence of secondary creep stage in some creep curves, i.e. the acceleration can occur just after end of the creep strain rate decrement without the observation of creep stages with constant creep strain rate. Eq. (5) is also able to describe such curves e.g. by employing larger 3 and 4 , such that the tertiary stage would become significant just after or even before the end that primary term affects. Following steps were carried out to develop this model: (1) Experiment creep curves were fitted to Eq. (5) to generate an initial dataset of model parameters. (2) Model parameters of a creep curve for any other temperatures and stresses are interpolated from this dataset. (3) The dataset was iterated several times with minor adjustment to satisfy the monotonicity requirement that creep curves monotonically increase with temperature and stress. (4) Lastly, both the time hardening model and strain hardening model were implemented and verified. Details are described in following subsections. For simplicity, we denote (T, ) as an operating condition with temperature T and stress load  for a creep process. 2.2.1. Construction of model parameter dataset from experimental creep curve fittings The REVISA experiment provides creep data operated under 600 ◦ C, 700 ◦ C, 800 ◦ C, 900 ◦ C, 1000 ◦ C, 1100 ◦ C, 1200 ◦ C and 1300 ◦ C. Under each temperature, 4 tensile experiments with different stresses were done, providing 32 (4 stress loads per temperature×8 temperatures) groups creep curves in total with different temperatures and stresses. These creep curves were fitted to Eq. (5). The fitting of each curve generated a set of model parameters 1 , 2 , m , 3 , 4 which can be thought as a representation of this creep curve in the form of Eq. (5). Fitting all the experimental curves then generated 32 sets of model parameters which form a model





dataset and can be denoted as Tj , j , 1,j , 2,j , m,j , 3,j , 4,j for the jth group where j = 1, . . ., 32. The fitting process is done with use of MATLAB. During the fitting, we applied the constrains that all i (i = 1, 2, 3, 4, m) are positive and for the 4 creep tests under the same temperature, tests with higher stress loads gives higher fitted model parameter i (i = 1, 2, 3, 4, m). The former constrain generally guarantees that three terms represent the primary, secondary and tertiary stages properly. The latter constrain guarantees that under same temperature, creep strain monotonically increases with stress load.

2.2.2. Interpolation The above fitting process provided 32 groups of model parameters with different temperature and stress from experimental curves. However, for most of the cases, the stress and temperature can not be the same with the experimental test conditions. Creep curves for those new conditions would be then interpolated or extrapolated. In original model, model parameters fitted from experimental curves were used to fit Eq. (4) and then model parameters for a new temperature and stress case would be got from these new fitted functions. In our modified model, model parameters are directly interpolated/extrapolated for any creep curve (T, ) from the known model parameter dataset Tj , j , 1,j , 2,j , m,j , 3,j , 4,j . This can be understood as, instead of approximating the distribution of i (i = 1, 2, 3, 4, m) with only one single equation (e.g. Eq. (4)), piecewise functions in the same or similar form are used here. We have several considerations on adopting this new treatment. Firstly, the fitting Eq. (4) in original theta projection method is suitable when the temperature range of experiment data sets are very narrow while the creep tests for 16MND5 steel were conducted with a large temperature interval ranging from 600 ◦ C to 1300 ◦ C. Besides, as structural performance can be highly nonlinear under high temperature (which is possible for the 16MND5 steel), one single equation like Eq. (4) may fail to reflect the nonlinearity or capture the performance over the whole range. Lastly, these engineering modelling are more or less empirical and different interpretations can exist as also seen in other modified theta projection methods. As an example to show the difference between parameter fitting and parameter interpolation, we fix the temperature to 1100 ◦ C for simplicity (under which 4 groups of experimental data are available), and demonstrate parameter fitting and interpolations over different stresses. The model parameters of experimental data are shown in Fig. 3 as blue circles at stress 4 MPa, 5.5 MPa, 9 MPa and 16 MPa. In the idea of the original theta project model i.e. Eq. (4), an exponential curve (because when T is fixed, Eq. (4) degrades to an exponential function of ) should be fitted from these known points as the red curve shows in Fig. 3. New  for any new  are then got from this fitted curve. This curve clearly captures the trend of , but the function values at given known points are not necessarily equal to the given known values. This discrepancy can be larger if the fitting is done with a 32-group data set. Using an interpolation method instead can well-avoid this problem as it keeps values of known points unchanged. The blue dash curve, green dash curve and black solid curve in figure show the interpolation with the classic Piecewise Cubic Hermite Interpolating Polynomial (PCHIP), piecewise linear function and piecewise exponential function respectively. The PCHIP agrees well with the experimental data and is smooth at the known points (mathematically continuous in first derivative at given known points [31]). But for the range  < 4 MPa,  decreases as the stress increases which might lead to a lower creep under a higher stress and is thought to be unphysical. Compared to the linear interpolation the piecewise exponential function approximates

P. Yu, W. Ma / Journal of Materials Science & Technology 47 (2020) 231–242

235

The interpolation expression for m = A3  n with two    over stress  known points 1 , m1 and 2 , m2 would be m = exp

 ln − ln

1

ln2 − ln1

lnm2 +

ln2 − ln lnm1 ln2 − ln1

 (18)

for The interpolation expression  over temperature    m = A4 e−b4 /T with two known points T1 , m1 and T2 , m2 would be m = m2

Fig. 3. Comparison of coefficient fitting and coefficient interpolation.

the fitted curve best. With this regard, the piecewise exponential function is adopted. For general cases in this model, if the function f (t, , T ) is used in the fitting method, the same function is then used in the interpolation method but in a piecewise form. Since nonlinear interpolation process over temperature and stress is complex (which may also involve matrix operations), we split the interpolation for any (T ’ , ’ ) into two sub-interpolation steps: interpolate over stress with fixed temperature, and then interpolate over temperature with fixed stress. In the first interpolation, we find out the two temperatures Tknown,small and Tknown,large most close to T ’ from the experimental temperature set {600, 700, 800,900, 1000, 1100, 1200, 1300}, then the first interpolation is done over stress load with the model parameter dataset obtained in Subsection 2.2.1 under   each  temperature to output the parameters for Tknown,small , ’ and Tknown,large , ’ . In the sec themodel outputsparameters for (T ’ , ’ ) from ond interpolation, Tknown,small , ’ and Tknown,large , ’ . The choice of interpolation functions depends on which i (i = 1, 2, 3, 4, m) we are interpolating and whether we are interpolating over temperature or stress. For 1 , 2 , 3 , and 4 , we assume a logarithm relation between i , temperature T and stress  which is similar to Eq. (4),

 

ln i = (a + b) (c + dT )

(14)

For the first interpolation over stress with fixed temperature, Eqn. (14) reduces to i = A1 eb1  where A1 and b1 are simplified coefficients after ignoring influence of temperature. For the second interpolation over temperature with fixed stress, Eqn. (14) reduces to i = A2 eb2 T where A2 and b2 are simplified coefficients after ignoring influence of stress. b  Theinterpolation  expression  for i = A1 e 1 with two known points 1 , i,1 and 2 , i,2 is −2

1 −

i = i,1 1 −2 · i,2 1 −2 Theinterpolation  expression  for i = A2 points T1 , i, 1 and T2 , i,2 is T −T2

(15) eb2 T

with two known

T1 −T

i = i,1 T1 −T2 ·  i,2 T1 −T2

(16)

Similarly, we assume following relation between m , temperature T and stress , n −b/T

m = A e

(17) n

The interpolation is done over the stress with m = A3 by fixing temperature, and then over temperature with m = A4 e−b4 /T by fixing stress.

1

T1



− T1 / T1 − T1 1 2



· m1

1 T

− T1

2

1

/ T − T1 1 2



(19)

2.2.3. Monotonicity assumption Intuitively, we would think materials deform more under higher stress and temperature load than those under lower stress and temperature. This is most frequently observed in the tensile experiments. Though, exceptions exist, one example of which is the effect of phase change: some materials under certain temperature around the phase change temperature may show much weaker structural performance. In this model, we just simply assumed that the creep strain will monotonically increase with temperature and stress loads. Or specifically, with the same creep time, the creep strain under higher temperature or stress should be larger than the creep strain with lower temperature or stress (the creep curve under higher temperature or stress should be above the curve with lower temperature or stress). When fitting the experimental data as described in Section 2.2.1, we already kept the trend that cases with higher stress load would have higher creep model parameters and give higher creep under same temperature. Obviously, this trend would also hold for the interpolated temperatures. Therefore, creep strain generally monotonically increases with stress load under same temperature. Then we only need to ensure the assumption that the creep strain monotonically increases with temperature given same stress load. However, as interpolation method is used (Section 2.2.2), this assumption could not always be satisfied, especially in the extrapolation range of the stress load. One example is illustrated in Fig. 4 (a). The red dots and blue dots are known parameters gained from experimental creep curves for corresponding high temperature TH and low temperature TL respectively. The red dash line and blue dash line (if y coordinate axis is in logarithm, then exponential functions can be plotted as lines) are the corresponding interpolation curves. These interpolation curves intersect at Point C when they extrapolate to a smaller stress range. Recall the interpolation steps, to get the model parameters for (TM , M ) by interpolation (assuming TL < TM < TH ), we fix temperature and interpolate parameters for (TH , M ) and (TL , M ) (shown as green diamonds) by the red and blue dash curves first, then fix stress and interpolate parameters for (TM , M ) from these two points (TH , M ) and (TL , M ). Then following situations are expected: 1) If M > C , the higher TM is, the larger the model parameters will be; 2) If M = C , the model parameters are fixed and independent of TM ; 3) If M < C , the higher TM is, the smaller the model parameter will be. It means under same stress, the model would predict larger creep strain for lower temperatures which violates our monotonicity assumption. To fix this issue, we introduced a parameter adjustment process as Fig. 4 (b), the idea of which is to minorly adjust the values of the four red and blue points (the model parameters from experimental curves) such that the intersection point C can be sufficiently small and the last situation would barely occur in the extrapolation range. Therefore, it ensures that under same stress load, model param-

236

P. Yu, W. Ma / Journal of Materials Science & Technology 47 (2020) 231–242

Fig. 4. Illustration of model parameter interpolation before (a) and after (b) parameter adjustment.

described by the time t. Therefore, it is a time-hardening model in default. The creep strain increment at time t during time interval dt is given by, dε = (1 2 e−2 t + m + 3 4 e4 t )dt

(20)

A strain-hardening model can be derived from Eq. (5) in an implicit form. Recall that the goal of strain-hardening model is to predict the creep strain rate given the current creep strain under certain temperature and stress. Specifically, when temperature and stress load change to (T, ) at the time t and creep strain ε0 , the model will get the creep curve of constant (T, ), and output how this curve will react under the current strain ε0 . Since the creep curve shows creep strain monotonically increases with time, the time and creep strain can be determined by each other. Therefore, in creep curve of (T, ), given the creep strain ε0 , we can uniquely find this reference time tref satisfying Eq. (21). One should notice that the reference time tref can be different from the current physical time. ε0 = 1 (1 − e−2 tref ) + m tref + 3 (e4 tref − 1 Fig. 5. Parameter adjustment process.

eters with higher temperature will be higher than that of lower temperature for most of the cases. We designed a process to verify that the model predicts higher creep strain under higher temperature. Fig. 5 illustrates the general verification process: for any (T, ) in the interested load range, the creep curve of (T + ıT, ) would be larger than the creep curve of (T, ). If not, parameter dataset would be adjusted and the verification would be repeated. This process would be iterated several times until this requirement is reached. 2.2.4. Hardening models Instead of being constant loads all the time, temperature and stress are also dynamically changing with time in real cases. Hardening models describe how creep strain rate reacts to the change of temperature and stress during a creep process. Time-hardening model and strain-hardening model are two mostly used hardening models in industry. Time-hardening assumes that creep strain rate depends only upon the time from the beginning of the creep process (and also temperature and stress if used in creep model) and ignores how large the previous creep strain is. Strain-hardening assumes that the creep strain rate depends only upon the current creep strain (and also temperature and stress if used in creep model). [32] In our creep model, both time-hardening model and strainhardening model were implemented. Either model can be activated for a simulation by user’s option. In Eq. (5), Creep strain is directly

(21)

where all s are model parameters of creep curve (T, ). Eq. (21) is a non-linear equation of tref . It can be solved with methods e.g. Newton-Ramphson [31]. Then the creep strain increment at this exact tref during time interval dt can be given as dε = (1 2 e−2 tref + m + 3 4 e4 tref )dt

(22)

Therefore, combining Eq. (21) and (22) gives the creep strain rate at creep strain ε0 . 3. Results and discussions Using the approach described above, we implemented our model in ANSYS Mechanical. Predicted creep curves were compared with experimental data. Monotonicity was tested over wide range of temperature and stress loads. Performances of interpolation with different load and hardening were also tested and demonstrated by two examples. 3.1. Model parameters Figs. 6–10 present the final model parameters i (i = 1, 2, 3, 4, m) under different temperatures and stresses that satisfy the monotonicity assumption discussed in Section 2.2.3. All the solid points correspond to conditions where experimental curves are available. The piecewise solid curves between each two points represent the interpolated model parameters’ distributions with stress under the same temperature. The curves are extrapolated to cover the stresses and temperatures outside the

P. Yu, W. Ma / Journal of Materials Science & Technology 47 (2020) 231–242

237

Table 2 Model parameters of 16MND5 steel. T, ◦ C

, MPa

1

2

m

3

4

600 600 600 600 700 700 700 700 800 800 800 800 900 900 900 900 1000 1000 1000 1000 1100 1100 1100 1100 1200 1200 1200 1200 1300 1300 1300

115 150 190 248 25 40 70 90 18 25 43 65 13 20 28 38 9.5 15.1 17.1 26 4 5.5 9 16 2.4 3 4.5 9 0.8 1.5 2.4

4.1834E-03 5.9140E-03 8.4217E-03 1.6079E-02 2.9780E-03 3.9922E-03 6.7311E-03 9.7005E-03 1.0348E-02 1.2902E-02 2.4751E-02 5.2346E-02 1.0083E-02 1.2275E-02 1.4513E-02 2.2112E-02 1.0875E-02 1.3650E-02 1.4596E-02 2.2740E-02 1.7639E-02 2.4542E-02 3.2077E-02 5.8488E-02 2.2000E-02 2.8215E-02 3.0899E-02 6.0000E-02 1.0491E-02 1.1000E-02 1.2522E-02

1.7232E-04 4.8234E-04 1.6858E-03 9.5421E-03 1.8967E-04 5.7390E-04 5.7364E-03 2.3568E-02 2.1541E-04 5.2565E-04 1.3645E-03 5.9998E-03 1.9601E-04 6.7889E-04 2.1769E-03 4.5980E-03 2.6436E-04 1.2748E-03 1.6571E-03 9.4482E-03 2.5526E-04 6.8978E-04 2.8725E-03 1.6143E-02 1.6208E-04 2.9000E-04 5.8918E-04 1.0000E-02 1.7527E-04 3.1231E-04 8.6418E-04

2.7065E-07 7.5762E-07 3.9798E-06 2.8533E-05 3.4986E-07 1.5199E-06 1.2006E-05 3.7239E-05 7.1048E-07 3.2562E-06 8.5905E-06 6.8685E-05 1.4445E-06 5.0377E-06 1.8192E-05 3.9623E-05 1.9002E-06 9.5097E-06 1.4492E-05 7.0730E-05 6.3946E-07 2.0702E-06 1.0056E-05 6.2898E-05 8.4246E-07 1.1371E-06 6.0454E-06 4.7687E-05 1.0585E-06 4.4331E-06 9.5580E-06

2.8597E-06 2.2255E-05 1.2075E-04 2.6533E-03 4.3127E-05 1.0535E-04 2.5094E-04 4.9625E-04 3.7357E-06 6.2670E-06 1.0143E-05 3.1538E-05 4.1699E-09 3.6222E-08 8.4617E-08 3.0168E-07 2.8196E-12 7.8363E-12 1.7249E-11 5.4360E-11 5.4089E-05 1.9320E-04 2.6409E-04 1.2582E-03 7.0622E-05 1.2289E-04 2.7503E-04 4.3082E-03 2.4026E-05 1.9025E-04 3.7761E-04

2.8597E-05 8.1706E-05 2.5015E-04 1.8175E-03 1.3250E-05 5.8358E-05 4.2144E-04 1.3444E-03 3.4132E-05 1.1508E-04 2.2648E-04 1.7876E-03 5.2644E-05 1.7578E-04 7.1897E-04 1.1517E-03 1.1092E-04 5.4466E-04 7.1105E-04 4.3141E-03 2.9073E-05 5.0273E-05 3.7349E-04 3.3856E-03 2.3447E-05 2.6945E-05 2.0327E-04 2.2080E-03 3.9530E-05 9.6177E-05 2.0316E-04

Fig. 6. ␪1 distribution with stress under different temperature.

Fig. 7. ␪2 distribution with stress under different temperature.

Fig. 8. ␪m distribution with stress under different temperature.

Fig. 9. ␪3 distribution with stress under different temperature.

238

P. Yu, W. Ma / Journal of Materials Science & Technology 47 (2020) 231–242

Fig. 10. ␪4 distribution with stress under different temperature.

experimental range. For 1 and 2 (Figs. 6 and 7), the parameters show exponential dependence on the stress for same temperature (linear curves in the linear-logarithm coordinate.). The dependence on stress also becomes steeper under higher temperature. The parameters also generally increase with temperature. It indicates generally the primary term in the model would be more drastic and decay faster at higher temperature and stress. For m (Fig. 8), each curve at different temperature shows linear trend in the log-log coordinate system which corresponds to a power dependence on stress. Recall that the creep strain rate is a power function of stress for the stable creep stage in Eq. (2). The secondary stage term reflects such pattern to some degree. As temperature increases, power trend is kept and the dependence on stress generally becomes more drastic. 3 and 4 determine the tertiary term that accelerates with time. 4 (Fig.10) generally increases with temperature and stress. For higher temperature, it increases more drastic with stress. 3 (Fig. 9) shows monotonic increase as the stress load, but not strictly monotonic increase as the temperature. It is also noticed that the parameters are not strictly monotonic in e.g. [26–28]. Essentially, monotonicity of all model parameters over temperature and stress is also not strictly required to satisfy the monotonicity of creep strain over temperature and stress. In our model, we also only require the overall monotonicity of creep strain other than the model parameters. Since 4 is positive and in the power of the tertiary term, the monotonic 4 would be more significant than 3 which may compensate the effect of 3 and result in the monotonicity of the overall creep strain. The detailed values of the model parameters are also listed in Table 2.

Fig. 11. Creep simulation predictions compared with experimental creep tests at T = 600 ◦ C (873 K).

Fig. 12. Creep simulation predictions compared with experimental creep tests at T = 700 ◦ C (973 K).

3.2. Simulation prediction vs. Experimental data Figs. 11–18 illustrate the 32 groups of creep curves predicted in ANSYS Mechanical with model parameters shown above compared with experimental data. The predicted curves proved to show the three stages creep process: the creep strain rate decreases at the early stage, then keeps stable and finally increase drastically. Good agreements have been achieved between the model predictions and experiment data for all the curves. 3.3. Interpolation performance Interpolation performance of the model was demonstrated by a following example. In this example, the model was to predict the creep curve of (950 ◦ C,14.0 MPa). Model parameters for cases (900 ◦ C, 13.0 MPa), (900 ◦ C, 20.0 MPa), (1000 ◦ C, 9.5 MPa) and (1000 ◦ C, 15.1 MPa) are available in the dataset for interpolation. The model firstly interpolated the parameters of (900 ◦ C, 14.0 MPa)

Fig. 13. Creep simulation predictions compared with experimental creep tests at T = 800◦ C (1073 K).

P. Yu, W. Ma / Journal of Materials Science & Technology 47 (2020) 231–242

239

Fig. 14. Creep simulation predictions compared with experimental creep tests at T = 900 ◦ C (1173 K).

Fig. 16. Creep simulation predictions compared with experimental creep tests at T = 1100◦ C (1373 K).

Fig. 15. Creep simulation predictions compared with experimental creep tests at T = 1000 ◦ C (1273 K).

Fig. 17. Creep simulation predictions compared with experimental creep tests at T = 1200 ◦ C (1473 K).

from parameters of (900 ◦ C, 13.0 MPa) and (900 ◦ C, 20.0 MPa), as well as the parameters of (1000 ◦ C, 14.0 MPa) from parameters of (1000 ◦ C, 9.5 MPa) and (1000 ◦ C, 15.1 MPa), then interpolated the parameters of (950 ◦ C, 14.0 MPa) from parameters of (900 ◦ C, 14.0 MPa) and (1000 ◦ C, 14.0 MPa). The predicted case was plotted as the light blue dash curve in Fig. 19, together with which all other intermediate cases were plotted as references. The red solid creep curve of (900 ◦ C, 14.0 MPa) lies between the blue solid creep curve of (900 ◦ C, 13.0 MPa) and the blue dash curve of (900 ◦ C, 20.0 MPa). Similarly, the red dash creep curve of (1000 ◦ C, 14.0 MPa) lies between the green solid creep curve of (1000 ◦ C, 9.5 MPa) and green dash curve of (1000 ◦ C, 15.1 MPa). The creep curve of (950 ◦ C, 14.0 MPa) lies between the red solid curve of (900 ◦ C, 14.0 MPa) and red dash curve of (1000 ◦ C, 14.0 MPa). These interpolated curves contain three stage creep processes, meanwhile they also obey the monotonicity assumption that higher stress and temperature result in larger creep. This indicates a reasonable interpolation prediction as we expected.

3.4. Monotonicity As mentioned in Section 2.2.3, the monotonic increase of creep strain with stress under same temperature is already satisfied and a detailed verification process was done to test the monotonic increase of creep strain with temperature under same stress. The achieved model performance is as follows. 1) For any (T, ) in the temperature range 273–1073 K with an interval of 1 K and the stress range 1–260 MPa with an interval of 1 MPa, its creep strain is higher than that of (T − 1K, ) for the first 12,000,000 s. 2) For any (T, ) in the temperature range 1073–1673 K with an interval of 1 K and the stress range 1–80 MPa with an interval of 1 MPa, its creep strain is higher than that of (T − 1K, ) for the first 12,000,000 s. 3) For any (T, ) in the stress range 0.1–26 MPa with a smaller interval of 0.1 MPa, for the full temperature range 273–1673 K

240

P. Yu, W. Ma / Journal of Materials Science & Technology 47 (2020) 231–242

Fig. 19. Interpolated creep curve predictions in ANSYS Mechanical. Fig. 18. Creep simulation predictions compared with experimental creep tests at T = 1300 ◦ C (1573 K).

with an interval of 1 K, its creep strain is higher than that of (T − 1, ) for the first 12,000,000 s. The verification process covered a wide range of temperature and stress, especially the lower ranges where the stress and temperature were small and extrapolation is used. Here we did not do the verification process one step further to cover even higher stress range or temperature range, because creep increases drastically under such high stress and temperature, and the results would become meaningless. With these regards, we think the monotonicity requirement is reached. 3.5. Hardening model performance The performance of time-hardening model and strainhardening model was also demonstrated through an example. The operating stress of the tensile creep simulation was 15.1 MPa. The operating temperature was 900 ◦ C at the beginning and jumped suddenly from 900 ◦ C to 1000 ◦ C at 36,000 s and that temperature kept till the end. Simulation results are shown in Fig. 20. The orange dash curve shows the change of temperature with time. The blue solid curve (the former part was identical to and covered by the green solid curve) and blue dash curve are the creep curves

for constant operating temperature T = 900 ◦ C and T = 1000 ◦ C situations, respectively. The green solid curve and green dash curve show the model response for this temperature shifting with timehardening model and strain-hardening model, respectively. Before the temperature shifting, the creep behavior with either hardening model was identical to the creep behavior of the T = 900 ◦ C till Point A. In time hardening model, the creep strain rate then shifted from Point A to Point B (which has the same time with A) and behaved the same as the case T = 1000 ◦ C as shown in the following green solid curve. In strain hardening model, the creep strain rate then shifted from Point A to Point C (which has the same creep strain with A) and behaved the same as the case T = 1000 ◦ C as shown in the following green dash curve. It should be noted that significant difference of both models was observed, partially because the temperature jump was drastic. In real cases, the temperature changes more mildly and the difference might be smaller. Overall, this prediction is exactly what we expect from the time-hardening model and strain-hardening model behaviors. 4. Conclusions Creep is an important phenomenon to be considered for structural components which are exposed to high temperature and stress. During a hypothetical nuclear power plant severe accident, the RPV would be experiencing creep deformation which may even

Fig. 20. Time hardening model and strain hardening in ANSYS Mechanical.

P. Yu, W. Ma / Journal of Materials Science & Technology 47 (2020) 231–242

lead to a vessel failure. A good creep model for the RPV steels would benefit the nuclear safety analysis and reactor design. In the present study we developed a new creep model for RPV steel 16MND5. It is adapted from the theta projection method. In this new model, creep curves are described as functions of time with 5 model parameters i (i = 1 − 4 and m). A dataset of model parameters was initially constructed by fitting all experimental curves into this function. To determine creep curves for new temperature and stress cases, a direct interpolation method was introduced to obtain corresponding model parameters from the dataset instead of using the original parameters’ fitting process which was only suitable for narrow temperature range. In the model development, we also assumed that the model would predict higher creep strain under higher temperature and stress loads. The model parameters in the dataset were minorly adjusted in such a way that the assumption was satisfied. The model has been successfully implemented into the ANSYS Mechanical. Based on the predictions of the new model, the following conclusions on model performance can be summarized:

241

This expression describes that ln i is a linear functionof  with  two unknown coefficients lnA1 and b1 . Two known points 1 , i,1





and 2 , i,2 would satisfy Eq. (A.1), which gives ln i,1 = lnA1 + b1 1

(A.2)

ln i,2 = lnA1 + b1 2

(A.3)

Eq. (A.2) and Eq. (A.3) uniquely determine lnA1 and b1 as follows, lnA1 = b1 =

2 lni,1 − 1 lni,2 2 − 1

(A.4)

ln i,1 − ln i,2 1 − 2

(A.5)

Plugging Eq. (A.4) and (A.5) back to Eq. (A.1) gives ln i =

2 ln i,1 − 1 ln i,2 ln i,1 − ln i,2 +  2 − 1 1 − 2

(A.6)

Alternatively, we can use the two-point linear interpolation correlation [31] which directly gives Eq.(A.7) from Eq. (A.2) and (A.3).  − 2 1 −  + ln i,2 1 − 2 1 − 2

(1) The model successfully captures all three creep stages. Good agreements between the prediction and experiment have been achieved for the total 32 groups of creep curves. (2) The performance of the new model in interpolating creep curves for new temperature and stress loads was demonstrated through a test example, and reasonable predictions were observed. This property ensures the model to predict reasonable creep curves under any interested temperature and stress loads. (3) Both time-hardening model and strain hardening model were realized through the creep model, with an expected performance. This property ensures that the model reacts reasonably to the transient change of temperature and stress loads during a creep process.

ln i = ln i,1

Above properties fulfill the requirements of a creep model for structural analysis. Over other previous creep models of 16MND5, this model shows the advantages of modelling three-stage creep and better agreement with experimental data. Further study could apply this model into analysis of RPV failure during a severe accident.

i = i,1 T1 −T2 ·  i,2 T1 −T2

Acknowledgments

It shows that ln m is a linear function of ln. The two-point linear interpolation correlation gives

This study was made possible by the support from the research programs of APRI, ENSI and NKS. The authors thank the staff at the Nuclear Power Safety Division of KTH for their support. The authors are also grateful to the support of the scholarship awarded by the China Scholarship Council (CSC). Appendix A This section gives detailed derivations of the interpolation expressions Eq. (15),(16),(17) and (18). For exponential interpolation function y = Aebx and power interpolation function y = Axb with two known points, the two known points can uniquely determine the two unknown coefficients A and b. The key steps of the derivation are to, (1) take the log of both sides of the exponential/power equation which then results in a linear equation in the logarithm form, (2) then get the expression of linear interpolation with the two known points, and (3) finally take the exponential of both sides of the interpolation expression. For the interpolation function i = A1 eb1  , the unknowns are A1 and b1 . The logarithm form reads ln i = lnA1 + b1 

(A.1)

(A.7)

Eq. (A.6) and (A.7) are essentially the same. Taking the exponential of both sides of Eq. (A.7) gives exp(ln i ) = exp(ln i,1

 − 2 1 −  + ln i,2 ) 1 − 2 1 − 2

(A.8)

Rearranging Eq. (A.8) gives 1 −

−2

i = i,1 1 −2 · i,2 1 −2

(A.9)

Eq. (A.9) is identical to Eq. (15).     For i = A2 eb2 T with two known points T1 , 1 and T2 , 2 , the derivation procedure is similar except that the variable is temperature T instead of stress . We directly replace all  with, then we get Eq. (A.10) which is identical to Eq. (16). T −T2

T1 −T

(A.10)

function m = A3  n with two known points  power    For the

1 , m1 and 2 , m2 , we take a log of both sides of the equation and will get

ln m = lnA3 + n ln

ln m = ln m1

(A.11)

ln2 − ln ln − ln1 + ln m2 ln2 − ln1 ln2 − ln1

(A.12)

Taking the exponential of above equations gives m = exp(

ln − ln1 ln2 − ln lnm2 + lnm1 ) ln2 − ln1 ln2 − ln1

which is identical to Eq. (18). For m = A4 e−b4 /T with two known points





(A.13)



T1 , m1



and

T2 , m2 , the derivation procedure is similar to that of Eq. (15)

except that the variable here is − T1 instead of stress  and i should be replaced by m . Then we will get m = m2

1

T1



− T1 / T1 − T1 1 2



· m1

1 T

− T1

2

1

/ T − T1 1 2



(A.14)

which is identical to Eq. (19). References [1] R.W. Evans, B. Wilshire, Introduction to Creep, Institute of Materials, London, 1993. [2] DoITPoMS, Creep Deformation of Metals, University of Cambridge, 2006, https://www.doitpoms.ac.uk/tlplib/creep/index.php. [3] F.N. Norton, The Creep of Steel at High Temperature, McGraw-Hill, 1929.

242

P. Yu, W. Ma / Journal of Materials Science & Technology 47 (2020) 231–242

[4] RCC-MR, Design and Construction Rules for Mechanical Components of FBR Nuclear Islands, AFCEN, Paris, 1985. [5] H. Bartsch, Steel Res. 66 (1995) 384–388. [6] L.M. Kachanov, Introduction to Continuum Damage Mechanics, Martinus Nijhoff Publisher, 1986. [7] A.J. Baker, M.P. O’Donnell, The Proc. 2nd Intern. Conf. on Integrity of High Temperature Welds, London, Oct., 2003, pp. 10–12. [8] B.F. Dyson, M. McLean, et al., in: A. Strang (Ed.), Microstructural Stability of Creep Resistant Alloys for High Temperature Applications, 1998, pp. 371–393. [9] ANSYS, ANSYS Mechanical APDL Material Reference Release 16.2, 2015, pp. 67–69. [10] SIMULIA, Abaqus 6.12 Analysis User’s Manual, Vol 3, 2012, pp. 264–275, Materials, Chapter 23.2.4. [11] S. Imaoka, User Creep Subroutine STI0704A. ANSYS Release:11.0, 2007. [12] J.-P. Mathieu, O. Diard, K. Inal, S. Berveiller, Proceedings of PVP2008, Chicago, USA, July, 2008, pp. 27–31. [13] R. Pesci, K. Inal, R. Masson, Mater. Sci. Eng. A 527 (2009) 376–386. [14] C. Sainte Catherine, Tensile and Creep Tests Material Characterization of Pressure Vessel Steel(16MND5) at High Temperatures (20 up to 1300 ◦ C), CEA, France, 1998. [15] K. Ikonen, Creep Model Fitting Derived From REVISA Creep, Tensile and Relaxation Measurements, VTT-Energy, Espoo, Finland, 1999. [16] E. Altstadt, T. Moessner, Extension of the ANSYS Creep and Damage Simulation Capabilities, Forschungszentrum Rossendorf, 2000, pp. 6–11. [17] H.-G. Willschütz, E. Altstadt, B.R. Sehgal, F.-P. Weiss, Nucl. Eng. Des. 208 (2001) 265–282. [18] H.-G. Willschütz, E. Altstadt, B.R. Sehgal, F.-P. Weiss, Ann. Nucl. Energy 30 (2003) 1033–1063.

[19] L.M. Kachanov, Izv. Akad. Nauk. SSR, Otd Tekh. Nauk. 8 (1958) 26–31. [20] Y.N. Rabotnov, Creep Problems in Structural Members, North-Holland, Amsterdam, 1969. [21] J.L. Chaboche, Nucl. Eng. Des. 105 (1987) 19–33. [22] E. Altstadt, H.-G. Willschuetz, B.R. Sehgal, F.-P. Weiss, Modelling of In-vessel Retention After Relocation of Corium Into the Lower Plenum, Forschungszentrum Rossendorf, 2005, pp. 96–100. [23] R.W. Evans, J.D. Parker, B. Wilshire, Int. J. Press. Vessel. Pip. 50 (1992) 147–160. [24] ECCC, Recommendations and Guidance for the Assessment of Creep Strain and Creep Strength Data, Vol. 5, ECCC Recommendations, 2013, pp. 38, Part 1 b [Issue 2]. [25] W.D. Day, A.P. Gordon, Proceedings of ASME Turbo Expo 2014, in: Turbine Technical Conference and Exposition, Düsseldorf, Germany, June, 2014, pp. 16–20. [26] M. Kumar, I.V. Singh, B.K. Mishra, S. Ahmad, A. Venugopal Rao, V. Kumar, J. of Materi Eng and Perform 25 (2016) 3985–3992. [27] W. Harrison, Z. Abdallah, M. Whittaker, Materials 7 (2014) 2194–2209. [28] M. Song, T. Xu, Q. Wang, W. Wang, Y. Zhou, M. Gong, C. Sun, Int. J. Press. Vessel. Pip. 165 (2018) 224–228. [29] C. Fu, Y. Chen, X. Yuan, S. Tin, S. Antonov, K. Yagi, Q. Feng, J. Mater. Sci. Technol. 35 (2019) 223–230. [30] C. Fu, Y. Chen, X. Yuan, S. Tin, S. Antonov, K. Yagi, Q. Feng, J. Mater. Sci. Technol. 35 (2019) 687–694. [31] Q.Y. Li, N.C. Wang, D.Y. Yi, Numerical Analysis, 5th ed., Tsinghua University Press, 2008, pp. 24, 222-226. [32] ANSYS, ANSYS Mechanical Advanced Nonlinear Materials V16.0, Lecture 2: Rate Dependent Creep, ANSYS Training Materials, 2015.