Three-DimensionalNumerical Simulation of CreepCrackGrowth Behavior for 316H Steel Using a Stress-Dependent Model

Based on detailed three-dimensional numerical simulation, creep crack growth behavior of C(T) specimen with dierent thicknesses of 316H steel was predicted using a stress-dependent creep ductility and strain rate model. ree regions were observed in the relation of creep crack growth rate versus fracture parameterC∗. e C(T) specimen with higher thickness exhibits higher CCG rate. e turning point 1 location from low C∗ region to transition C∗ region increases with increasing thickness, while that of turning point 2 seems to be independent of specimen thickness. Based on the nite element results, constraintdependent turning point 1 location and creep crack growth rate equations were tted. More accurate and realistic life assessment may be made when the stress-dependent model and the constraint eect were considered for creep life assessments of hightemperature components subjecting to a low applied load.

1. Introduction e in uence of constraint level on creep crack initiation (CCI) and creep crack growth (CCG) plays an important role and has been considered a vital issue in high-temperature component life assessments. Because the constraint can dramatically alter fracture behavior of materials and structures, it is indispensable to quantify constraint accurately and to incorporate constraint e ects in structural integrity assessments of high-temperature components. For this purpose, constraint parameters, such as Q [1][2][3], T z [4], R [5], R * [6], and Ac [7], have been proposed leading to twoparameter approach C * -Q, C * -T z , C * -R, C * -R * , and C * -A c and three-parameter approach C(t)-T z -Q to describe the creep-tip stress and strain rate elds.
Many experimental and theoretical assessments have shown that the constraint can a ect creep crack growth (CCG) rate. e CCG rates increase with the increase of crack depth [8] and specimen thickness [9][10][11]. e investigation of specimen type e ect on CCG rates demonstrated that C(T) specimen with high constraint has the fastest CCG rate and center cracked panels (CCP) with low constraint have the lowest CCG rate [3,12]. It also has been indicated that the CCG rate examined in C(T) specimens is signicantly faster than that of M(T) specimens at a given C * value for various steels [12,13].
Signi cant work has been done in recent decades to test and predict the creep crack growth behavior of type 316H stainless steel. Creep ductility exhaustion model has been extensively employed to predict CCG rate [14][15][16][17], and theuniaxial creep ductility is usually assumed to be a constant for a given temperature. However, a lot of experimental results and analyses have shown that the creep fracture mechanism depends on stress levels, which leads to the stress dependence of the creep ductility of materials [18][19][20]. To obtain an accurate prediction of CCG rate under a wide range of stress level, the stress-dependent ductility and strain rate model need to be employed in nite element (FE) simulation.
In recent works [20,21], the CCG behavior of 316H steel was predicted using the stress-dependent creep ductility and strain rate model. Two-dimensional nite element analyses were conducted for C(T) [20] and CS(T), SEN(B), SEN(T), DEN(T), and M(T) specimens [21], and good agreements have been found between the predicted CCG rates and available experimental data. Knowing that the CCG rates in steels generally increase with the increasing specimen thickness [9][10][11], investigation of the out-of-plane constraint effect induced by specimen thickness on CCG behavior of 316H steel using stress-dependent creep ductility and strain rate model has not been examined. erefore, in this work, three-dimensional finite element simulations based on stress-dependent model have been conducted for C(T) specimen with different thicknesses to account for the effect of out-of-plane constraint on CCG behavior of 316H steel.

Material
Properties. e elastic-plastic behavior of 316H at 550°C is shown by the true stress-strain curve in Figure 1 [22]. Young's modulus E and yield stress σ y (which is often taken as 0.2% proof stress σ 0.2 ) at 550°C are 140 GPa and 170 MPa, respectively. Poisson's ratio ] is assumed as 0.3, which is the same as that in literature [16][17][18]23]. In general, Poisson's ratio changes very little with temperature for hightemperature steels [24,25]. A stress-dependent average creep strain rate over a wide range of stresses has been observed for many engineering steels and alloys, which is determined by different creep deformation mechanisms [20,[26][27][28][29][30][31]. Two-regime Norton creep model has been proposed to characterize the different creep deformation behaviors, where the average creep strain rate is correlated by normalized stress σ/σ 0.2 in recent studies [19-21, 26, 32-35]. For given stress and temperature, the average creep strain rate is defined as the ratio of creep ductility ε f to rupture time t r . e stress dependency of the average creep strain rate _ ε A can be described using a tworegime Norton power law relationship: where _ ε A is in h −1 and A A is in MPa -n h −1 . e stress-dependent average creep strain rate for the 316H material at 550°C is illustrated in Figure 2 [20], and power law coefficients and exponents for low stress (σ/σ 0.2 < 1.035) and high stress (σ/σ 0.2 ≥ 1.035) are listed in Table 1 [20].

Finite Element Models.
ree-dimensional (3D) FE analyses have been conducted on plain sided C(T) specimens with different thicknesses. e geometry of the C(T) specimen is shown in Figure 3. Here, a 0 is the initial crack depth, B is the specimen thickness, W is the specimen width, and H is the distance between the clamps. e specimen width W is fixed at 50 mm, and H equals 0.71 W. e initial crack length a 0 /W is set as 0.5, and four thicknesses B/W � 1/8, 1/4, 1/2, and 1 are adopted to investigate the effect of specimen thickness on CCG behavior using stress-dependent model.
Due to symmetry in geometry, only a quarter of C(T) specimen was built. e symmetry boundary condition is      [20]. applied on the un-cracked ligament and middle plane (z/ B � 0) of the specimen. e load is applied to the loading hole as a distributed load. e typical 3D finite element model constructed for C(T) specimen is illustrated in Figure 4(a), and the local mesh distribution around the crack tip is shown in Figure 4(b). e FE simulations were carried out using ABAQUS code [36] with eight-node linear reduced integration element (C3D8R). e element size of 100 µm was chosen in the crack growth zone, which is similar to the average grain sizes of 316H stainless steel [17].

Creep Damage Model and Creep Crack Growth
Simulation. e creep damage accumulation ahead of the crack tip based on creep ductility exhaustion approach has been widely used in CCG simulations [14][15][16][17]. e accumulated creep damage parameter ω is defined as follows: where _ ω is the creep damage rate, _ ε c is equivalent creep strain rate, and ε * f is the multiaxial creep ductility. e damage parameter ω equals 0 at time t � 0 and 1 when a failure occurs. When the value of ω at a given element reaches 1, it is considered fully damaged and its load-carrying capacity is reduced to a value of close to zero. In this paper, the elastic modulus is reduced to 1 MPa to simulate the loss of loadcarrying capacity which has been used in recent works [14,17,20,21,[37][38][39].
Many models have been proposed to estimate the multiaxial creep ductility ε * f from uniaxial creep ductility ε f . e one which has been used in this study is the Cocks and Ashby void growth model [40] which is shown as where n is the creep exponent for power law creep and σ m /σ e is the ratio between mean (hydrostatic) stress and equivalent (von Mises) stress which is known as stress triaxiality.
A lot of experimental results and analyses have shown that the creep fracture mechanism depends on stress level, which leads to the stress-dependent creep failure strain (creep ductility) [20,33,34,41,42]. At high stress, the creep fracture progress is dominated by plasticity controlled void growth [43,44]. e cavity growth rate is linear with creep strain rate and inversely proportional to creep rupture time, resulting in a constant higher creep ductility [41]. At low stress, the fracture mechanism is constrained diffusion cavity growth. e local deformation due to the growth of inter-granular voids exceeds the deformation rate of the surrounding material. By virtue of redistributed stresses, the cavitating boundaries transverse to the maximum principal stress unload until the local strain rate equals the remote strain rate caused by the applied stress [45]. e cavity growth rate is similar to that at high stress, resulting in a constant lower creep ductility [41]. At a narrow transition stress region, the fracture mechanism is diffusion controlled cavity growth and the cavity growth is linear with stress, giving the creep ductility as a function of stress [41]. is stressdependent creep ductility has been identified for 316H steel [20], Cr-Mo-V steel [35], and T92, P92, T122, and P122 steels [46][47][48].
e stress-dependent creep ductility normalized at the lower shelf (LS), the upper shelf (US), and stress-dependent transition region in between is shown in Figure 5 [20]. A large ductility test data scatter collected from NIMS [49] and EDF Energy [50] has been found for 316H due to various material casts and batch differences. An average value of ductility has been calculated at high normalized stress shown in Figure 5. However, the ductility test data at low stress are not available. An estimated creep ductility at low stress has been obtained by assuming the same ductility difference between US ductility and LS ductility at 650°C and 550°C [20]. en, the stress-dependent region between high stress and low stress was generated and illustrated by dashed lines in Figure 5. e stress-dependent creep ductility has been used to predict creep crack growth and creep crack initiation behaviors. A good agreement has been obtained for both short-and long-term creep crack growth test data [20,21] and creep crack initiation times [32] for different specimen  geometries. e creep ductility values for a wide stress range have been fitted and expressed in e CCG simulations were performed using ABAQUS [36] code with a user subroutine USDFLD. e stress-dependent average creep strain rate and creep ductility illustrated in Figures 2 and 5 and given in Table 1 and equation (4) are implemented in the subroutine to calculate the creep damage at each time interval. When an element is fully creep damaged, the crack advances by reducing the load-carrying capability of that element until numerical convergence occurs. For each specimen, the range of initial stress intensity factor K in for CCG simulations is in the region between 3.46 MPa m 1/2 to 21.19 MPa m 1/2 . e creep crack growth length was estimated by numbers of fully damaged elements. e C * value was calculated by using the equation in ASTM E1457 [51]: where F is the applied load, _ Δ is the creep load line displacement rate, B n is the specimen net thickness, W is the specimen width, and a is the current crack length. In equation (5), H 1 and η are nondimensional geometry-dependent parameters and the solutions of which can be found in ASTM E 1457 [51]. e criteria in ASTM E 1457 [51] were applied to ensure that C * is a valid parameter to describe the CCG rate.

NSW Creep Crack Growth Model.
e well-known NSW model proposed by Nikbin, Smith, and Webster [52,53] has been widely used to predict the steady-state CCG rate from uniaxial creep data. e NSW model is based on a creep ductility exhaustion approach, and the CCG is predicted to occur when a critical level of damage is attained at a characteristic distance, r c , ahead of the crack tip. e NSW CCG model provides different prediction lines for plane stress (PS) and plain strain (PE) conditions due to its stress state dependence. Under steady-state creep condition, the CCG rate can be predicted by the NSW model denoted as _ a NSW and expressed as follows: where ε * f is the multiaxial creep ductility and A, n are power law creep coefficient and exponent. It is suggested that ε * f in (6) is taken as the uniaxial creep ductility ε f for plane stress condition and ε f /30 for plane strain condition [54]. e creep process zone distance r c in NSW model is usually taken as the average material grain size, thus is taken as 0.1 mm in this study [17]. Note that the _ a NSW is not sensitive to this value [3,52,53,55]. e dimensionless constant I n can be estimated by equations in [56] for PS and PE conditions.

Validation of the Simulation Method
To investigate CCG behavior of C(T) specimen with different thicknesses in 316H steel, the simulation method based on a stress-dependent model should be validated firstly. Comparisons between available CCG test data and FE results for short-term CCG test specimen (1TC(T) CT1 [17,57]) and long-term CCG test specimen (CT1B and CT20 [57]) in literature were made. e short-term [12,17,58,59] and long-term creep test band [60,61] in Figure 6 is gathered from extensive C(T) specimens creep crack growth test data of 316H at 550°C. Note that the scatter band is so wide up to a factor of around 10 for short-term test data [17]. e creep constraint effect plays a key role of leading such a scatter band. erefore, it is quite necessary to consider creep constraint effect in structural integrity assessments and give more accurate remaining life predictions of structures or components. e comparison of CCG rate between FE results and creep test data is shown in Figure 6. e FE results were calculated based on a stress-dependent model. It can be seen in Figure 6(a) that both FE results and short-term creep test data are within the short-term test data band. At the beginning of CCG, the FE result corresponds well with test data. As the crack advances, the CCG rate of FE result is slightly lower than that of test data. Due to high applied load in short-term CCG test, large plasticity occurs at the crack tip. In FE analysis, the crack growth may be dominated by high stress constitutive equation and upper shelf creep ductility (ε f �13.6%), which is higher than the test data of creep ductility (ε f �6%-9.3% [17]). Higher creep ductility leads to lower CCG rate. Similarly, for long-term CCG, the crack growth simulation was dominated by the low stress constitutive equation and lower shelf creep ductility (ε f �0.9%), which shows higher CCG rate than test data. Note that the CCG rate of FE result lies above the long-term test band. In general, the CCG results of FE simulations based on the stress-dependent model are consistent with that of test data. erefore, this simulation method can be used to investigate the effect of specimen thickness on CCG behavior of 316H. Figure 7 shows the typical creep damage distribution (CCG profile) of plane sided C(T) specimen along the thickness direction. e parameter z is the distance from middle plane to surface plane, which means the middle plane and surface plane are denoted as z/B � 0 and z/B � 0.5, respectively. It can be seen in Figure 7 that the largest CCG length is observed at the middle plane z/B � 0 while the smallest one is clearly seen at surface plane z/B � 0.5. is is caused by different C * and constraint levels through the specimen thickness. e specimen region with higher crack tip fracture parameter C * and constraint level has higher CCG rate. e distribution of crack tip fracture parameter C * along thickness direction has been given with the highest C * value at mid-thickness plane and the lowest C * value at the surface plane [62]. e constraint level decreases as z/B increases [23], leading to the highest constraint level at the middle plane (z/B � 0) and the lowest constraint level at the surface plane (z/B � 0.5). Figure 8 shows the creep crack growth rate versus C * from about 1 × 10 −8 MPa m/h to 1 × 10 −2 MPa m/h. It can be seen that similar creep crack growth trends are obtained for C(T) specimen with B � 12.5 and B � 25. e whole C * region can be divided into low C * region, transition C * region (between point 1 and point 2), and high C * region. is is consistent with the experimental observation of 316H steel  Mathematical Problems in Engineering 5 [60] and Cr-Mo-V steel [10]. e simulation results of Cr-Mo-V steel in previous works also show such a trend [35,37,38,63]. e CCG rate in low C * region is higher than that of extrapolation from high C * region, which means nonconservative prediction will be made if the extrapolated CCG rate from the lab test data is used for life prediction of practical components operated at low load level. us, accurate CCG prediction should be made by using actual CCG rate data at low C * region. e CCG behavior of C(T) specimen with B � 6.25 mm and B � 50 mm is similar to that in Figure 8.

Effect of Specimen
ickness on CCG Behavior. Figure 9 shows the creep crack growth behavior of C(T) specimen with different specimen thicknesses for a wide range of C * level. e NSW model prediction lines are also included in Figure 9. e prediction lines calculated using the constant lower shelf (LS) creep ductility and average creep strain rate exponent at low stress level for PE and PS conditions are denoted as NSW-PE-LS and NSW-PS-LS, respectively. In addition, the NSW prediction lines corresponding to the upper shelf (US) creep ductility and average creep strain rate exponent at high stress level for PE and PS conditions are denoted as NSW-PE-US and NSW-PS-US, respectively.
It can be seen that different CCG rates were obtained at different load levels C * . At both low C * and transition C * regions, the CCG rate increases with specimen thickness, which is the same as test data in the literature [9,10]. e location of turning point 1 from low C * to transition C * increases with specimen thickness. For turning point 2 between transition C * region and high C * region, the effect of thickness on location variation is less conspicuous than that of point 1. At high C * region, with increasing load level C * , the CCG rate initially increases with thickness and then slightly decreases after point 3, which was observed in creep test of Cr-Mo-V steel [10].
is is caused by a loss of   constraint at high C * levels where large plastic deformation occurs [61,64]. Note that the NSW CCG prediction is strongly dependent on creep ductility and thus higher CCG rates are predicted using NSW PE prediction lines as much lower creep ductility is employed in NSW models for the PE condition than that for PS condition. For short-term CCG data at high C * region, the NSW-PS-US prediction line falls close to the FE results, while more conservative results are predicted using NSW-PE-LS line for long-term CCG data at low C * region.

Constraint-Dependent CCG Rate Formula
Recently, a new load-independent constraint parameter R * has been proposed based on crack tip stress field by authors [6,65], which can fully characterize out-of-plane constraint when B/W ≤ 1 [23]. In this paper, the constraint parameter R * is used to establish a constraint-dependent CCG rate formula. e method for establishing the correlation of CCG rate with creep constraint parameters has been given in detail in the recent work of authors [6,7]. e load-independent creep constraint parameter R * at steady-state creep has been investigated and expressed as follows [6,65]: where σ 22 and C * 1 are the opening stress and C * value in the under-evaluated cracked specimen or component, respectively, and σ 22 , C(T) and C * 2 are the opening stress and C * value of the standard C(T) specimen in PE condition for obtaining reference stress field, respectively.
For 3D practical analyses, an average value of constraint parameter along crack front is usually used to characterize constraint in addition to that at mid-plane [65][66][67].
erefore, an average value of the parameter R * along specimen thickness, denoted as R * avg , was used to correlate with CCG rate. In this paper, 10 elements have been meshed along 3D specimen thickness. e R * value at every thickness layer from the mid-plane (z/B � 0) to the surface plane (z/B � 0.5) was obtained by (7). en, the average value R * avg can be calculated as follows [6]: where z is the distance from middle plane along specimen thickness. Figure 10 shows the relation between the turning point locations in Figure 9 and the constraint parameter R * avg . As actual structures or components are subject to low load in practical use, only the locations of turning points 1 and 2 are shown in Figure 10. A linear relation exists between the v value and constraint parameter R * avg for turning point 1, while the C * values for turning point 2 are constant, which is constraint independent. An approximate C * value of 5 × 10 −5 MPa m/h was observed for point 2. If the formula between turning point 1 location and constraint parameter R * avg was fitted, it can be predicted for specimens with other thicknesses and then the size of the transition region can be obtained. e fitted equation between turning point 1 and constraint parameter R * avg is expressed as follows: Figure 11 shows the relation between CCG rate ratio _ a/ _ a 0 and constraint parameter R * avg . It can be seen that a linear relation exists between _ a/ _ a 0 and (1 − R * avg ) on a log-log scale. e fitting formula is shown as follows: Transition C * region C *=2E-5 MPam/h (b) Figure 11: e relation between CCG rate ratio and constraint parameter for C(T) specimen at (a) low C * region and (b) transition C * region.
Mathematical Problems in Engineering 7 where _ a 0 is the CCG rate of standard reference C(T) specimen with W � 50 mm and a/W � 0.5 in PE condition. It has been validated that the relation is load-independent in a specific C * region in previous work [7]. is means the CCG rate can be predicted under different C * levels within a certain C * region. e formula above may be used to predict CCG rate data for various specimens with different thicknesses and creep life of high-temperature components accounting for constraint effect. As long as the constraint parameter R * avg is obtained at a certain load level C * using FE analysis, then the CCG behavior can be predicted by (10). To validate the formula above, a new C(T) specimen with a/W = 0.5 and B = 20 mm is modeled. A comparison between the CCG rates calculated based on FE results obtained by stress-dependent model and that predicted by (10) with constraint parameter R * avg was made and is shown in Figure 12. e NSW model prediction lines are also included. It shows a good agreement of CCG rate data between those two methods at low C * region and transition region, while the NSW PE prediction lines provide a more conservative CCG rate for a wide range of load level. In addition, the prediction of point 1 location by equation (8) (C * =3.55 × 10 −6 MPa m/h) is very close to that obtained by FE simulation based on stress-dependent model. It means the CCG rate and location of the turning point at low C * region can be successfully predicted by constraint-dependent equation.
erefore, more accurate CCG data and reasonable life assessment can be made for high-temperature components subjecting to a low applied load.

Conclusion
In this paper, three-dimensional finite element analyses have been conducted on C(T) specimen with different thicknesses using stress-dependent creep ductility and strain rate model. Based on the result, the main conclusions are as follows. e whole C * region of CCG rate can be divided into low C * region, transition C * region, and high C * region. e CCG rate in low C * region is higher than that of extrapolation from high C * region, which means nonconservative prediction will be made if the extrapolation CCG rate from the laboratory test data is used for life prediction of practical components operated at low load level.
At both low C * and transition C * regions, the CCG rate increases with specimen thickness, while a shift from increase to decrease with increasing thickness appears at high C * region. e C * values at the location of turning point 1 increase with specimen thickness. e variation of C * values at turning point 2 is less conspicuous than that of point 1, leading to the smaller size of transition region with increasing thickness.
Based on constraint parameter R * avg , the constraint-dependent equations of turning point 1 location and CCG rate at low C * region and transition C * region were fitted. By using the equations, the C * value at turning point 1 from low C * region to transition region may be predicted and the transition region size can be approximately obtained. e CCG rate can be easily predicted by using the constraintdependent CCG rate formula at both low C * region and transition C * region.

Data Availability
e data used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest
e authors declare that they have no conflicts of interest.