Probability Prediction of Crack Driving Force of Crack Tips of TP304 Stainless Steel Based on Kriging Surrogate Model

Abstract

To explore the influence of randomness of materials and loads on the crack driving force of TP304 stainless steel, a probability prediction for crack driving force through the elastic-plastic finite element method (EPFEM) coupled with the Kriging surrogate model was proposed. To improve the efficiency of finite element analysis, MATLAB was used to further develop the pre-processing and post-processing procedures of ABAQUS software to realize the automatic change of random specimens, batch calculation, and automatic analysis of probability prediction results. The statistical distribution law of the crack driving force of TP304 stainless steel material under the action of random factors was obtained, as well as other probability characteristics, including failure probability, failure probability density function, cumulative probability density function, etc. The sensitivity of each random factor was analyzed. Finally, the effectiveness and efficiency of the proposed method were analyzed, compared with those of the Monte Carlo method. Results show that the randomness of load and material parameters can significantly influence the driving force of crack tips of TP304 stainless steel, thereby affecting the failure probability of TP304 stainless steel. The load and strain hardening exponent present the most obvious effect on the dispersion of crack driving force of austenitic TP304 stainless steel.

Keywords

probability prediction; crack driving force; J-integral; Kriging surrogate model

Nuclear energy as a low-carbon energy source plays an essential role in the energy conservation^[1]. However, there are some longstanding concerns about the long-term operation of nuclear power plants (NPPs), such as structural reliability and integrity^[2–3]. It is well-known that the corrosive environment, sensitivity of material, and high tensile stress near the crack tip are critical factors influencing the stress corrosion cracking (SCC) in high-temperature aqueous environment^[4]. Since the austenitic stainless steel is widely employed in the vessels and pipes of NPPs, SCC sensitivity (residual stress or mechanical heterogeneity) of the materials has been researched^[5–8]. The driving force of crack tips of TP304stainless steel is one of the critical factors influencing the stress corrosion cracking behavior. However, the crack driving force near the tip region of TP304 stainless steel is still obscure.

Many mechanisms of SCC crack propagation are proposed. Among them, the crack driving forces including the plastic strain and J-integral near the crack tip are optimal fracture parameters to describe the mechanical state of the crack tip. Currently, the deterministic models of linear-elastic fracture mechanics (LEFM) and elastic-plastic fracture mechanics (EPFM) are widely used to investigate the crack driving forces^[9–10]. Probability analysis is also conducted to estimate the statistics and reliability of crack driving forces^[11–12]. Concerning the complexity and time-consuming of the micro-mechanical state calculations, the probability analysis based on complex EPFM models is rarely employed.

Various probabilistic methods are proposed for the fracture mechanics, such as first- and second-order moment methods, Latin hypercube sampling method (LHSM), and Monte Carlo method (MCM)^[13–15]. Among these methods, the estimation formulas of the crack driving force are a function related to the material properties, crack size, and load. However, these methods cannot be applied to complex crack driving force analysis. Traditionally, the evaluation of crack driving force can be performed by the elastic-plastic finite element method (EPFEM) under any crack shape and loading conditions^[16]. Since Kriging surrogate model only requires simple evalua-tions in the deterministic model, it has high computational efficiency in dealing with the long-running and complex engineering models^[17–19]. However, the crack driving force analyses involving EPFEM and Kriging surrogate models are rarely reported.

In this research, a computational method for randomness prediction of the crack driving force through EPFEM coupled with Kriging surrogate models was proposed, which could improve the prediction accuracy of driving force of crack tips of TP304 stainless steel in the essential structures for NPPs.

1 Calculation Model

1.1 Geometric model

One-inch compact tension (1T-CT) specimen suffering a constant load is standardized by ASTM in the SCC experiments under high-temperature aqueous environments^[20]. The schematic diagram of 1T-CT specimen is shown in Fig.1. Therefore, to investigate the probability distribution characteristics of the driving force of crack tips of TP304 stainless steel, the simulated numerical tests with 1T-CT specimens were conducted according to ASTM standards^[21]. The geometry and size of 1T-CT specimen are shown in Fig.2.

Fig.1  Schematic diagram of 1T-CT specimen

Fig.2  Geometry and size of 1T-CT specimen (W=50 mm, a=0.5W, c=1.5 mm)

1.2 Material model

To explore the probability distribution characteristics of crack driving force through EPFEM, the material model should be defined. The constitutive law of the stress (σ)- strain (ε) relationship beyond the yielding stage can be expressed by the Ramberg-Osgood relationship, as follows:

$$\frac{\epsilon }{{\epsilon }_{\mathrm{0}}}=\frac{\sigma }{{\sigma }_{\mathrm{0}}}+\alpha {\left(\frac{\sigma }{{\sigma }_{\mathrm{0}}}\right)}^n$$

(1)

where ${\sigma }_{\mathrm{0}}$ is the yield strength; ${\epsilon }_{\mathrm{0}}$ is the reference strain; α is the dimensional material constant; n is the material strain hardening exponent. If the modulus of elasticity is expressed as E, ${\sigma }_{\mathrm{0}}$=E${\epsilon }_{\mathrm{0}}$. The specimen material was set as TP304 stainless steel and the high temperature was set as 288 °C. In this research, the reference stress ${\sigma }_{\mathrm{0}}$=154.78 MPa and Poisson's ratio $\nu$=0.3 were set as the deterministic parameters. Thus, there are only two independent variables in Eq.(1). Therefore, the random parameters involved E, α, n, and load P. To explore the influence of load on the crack driving force, the load was varied from 0 N to 1000 N. The mean coefficient of variation (COV) and probability distributions of these random variables are shown in Table 1^[11].

Table 1  Statistical properties of random variables^[11]

Random parameter	Mean value	COV	Probability distribution
Elastic modulus, E	206.8 GPa	0.05	Gaussian
Constant, α	8.073	0.439	Lognormal
Strain hardening exponent, n	3.8	0.146	Lognormal
Load, P	Variable	0.10	Gaussian

1.3 Finite element model

1T-CT specimen was adopted in the experiments, and the loading process was simulated by the commercial finite element code ABAQUS (Version 6.14). Since the front-end of crack along the thickness direction of specimen was mainly controlled by the plane strain condition in the pipeline, a two-dimensional plane-strain model was analyzed. A seam with length of 1.5 mm was defined as the initial crack. The finite element mesh was constructed, as shown in Fig.3. 5176 second-order elements from the ABAQUS element library were employed, and the focused elements were adopted in

Fig.3  Schematic diagrams of finite element mesh of 1T-CT specimen: (a) global mesh model and (b) refined mesh of crack tip

the vicinity of the crack tip. The element type CPE8RH was adopted, and a mixed formulation element was typically employed to address the incompressibility constraint. The deformation theory was employed in the material model.

2 Research Methods

2.1 LHSM theory

LHSM involves the multi-dimensional stratified-random sampling method based on its variance reduction technique. Compared with simple sampling method, full coverage of the range of variables can be satisfied^[22–23]. If the variables are independent of each other, the procedure for this method can be defined as follows.

(1) Divide the cumulative distribution of each variable into N equal probability intervals.

(2) A value is selected from each interval randomly. Based on Ref.[24], the cumulative probability for the j-th interval can be defined as $P_j=\left(\frac{\mathrm{1}}{N}\right)r_i+\frac{\left(j-\mathrm{1}\right)}{N}$, where $r_i$ is a random number of 0–1.

(3) Using the inverse value of distribution function F^-1, the probability values are transformed into x, namely, $x=F^{-\mathrm{1}}\left(P\right)$.

(4) Each variable x has N values. Different variables have different N values.

2.2 Kriging theory

Kriging theory is usually used to predict the response values of discrete input design points^[25]. Based on the observation points, the response of unobserved points can be predicted by the Kriging model. For the function $\widehat{Y}\left(X\right)$, the Kriging model can be represented as follows:

$$\widehat{Y}\left(X\right)=f_{\mathrm{r}}{\left(X\right)}^T{\beta }_{\mathrm{K}}+\epsilon \left(X\right)$$

(2)

where $f_{\mathrm{r}}\left(X\right)$ is a vector of regression functions; ${\beta }_{\mathrm{K}}$ is a vector of regression coefficients; $\epsilon \left(X\right)$ is a stationary Gauss process. Thus, the covariance can be expressed as follows:

$$\mathrm{C}\mathrm{o}\mathrm{v}\left\{\epsilon \left[X^i,\epsilon \left(X^j\right)\right]\right\}={\sigma }_{\mathrm{\epsilon }}^{\mathrm{2}}R\left(X^i,X^j\right)$$

(3)

where ${\sigma }_{\mathrm{\epsilon }}^{\mathrm{2}}$ represents the variance of Gauss process; $R\left(X^i,X^j\right)$ denotes the correlation function between $X^i$ and $X^j$. $R\left(X^i,X^j\right)$ is a relationship function which satisfies specific conditions (symmetric; positive semi-definite). In this research, the anisotropic Gaussian model was adopted, as follows:

$$R\left(X^i,X^j\right)={\displaystyle \prod _{m=\mathrm{1}}^n}\mathrm{e}\mathrm{x}\mathrm{p}\left[-{\theta }_m{\left(x_m^i-x_m^j\right)}^{\mathrm{2}}\right]$$

(4)

where $\theta$ is the length; n is the number of random variables.

2.3 Proposed procedure

The framework of the proposed procedure is shown in Fig.4. The main procedure involves EPFEM coupled with Kriging surrogate model. The proposed procedure has the advantage of high efficiency which is attributed to Kriging surrogate model with acceptable computational efforts. Kriging model was constructed by DACE tool in MATLAB^[26]. Multiple sampling calculations were required, and MATLAB was used to further develop the pre-processing and post-processing of ABAQUS software, which could improve the efficiency of EPFEM simulation. Table 2 shows the detailed step of this proposed method.

Fig.4  Framework of proposed method

Table 2  Step description of proposed method

Step	Description
1	Create a deterministic model (noted as Md) of 1T-CT specimen by EPFEM
2	Further development of pre-processing and post-processing of ABAQUS with MATLAB
3	Generate a large number of input specimens by LHSM and get the corresponding model responses adopting Md
4	Construct a Kriging model with current specimens and responses
5	Predict the model response of specimens according to the distribution defined in Table 1 using the current Kriging model
6	Probability prediction of crack driving force through the obtained Kriging model

3 Results and Discussion

3.1 Plastic zone and plastic strain in crack tip

Since the crack driving forces play an essential role in SCC behavior, the local plastic strain and plastic zone around the crack tip are usually adopted as the main affecting factors to predict the SCC growth rate^[27–28]. Therefore, the influence of random parameters on the local mechanical parameters was investigated in this study.

A circular area with diameter of 3 mm around the crack tip was concerned. The load was set as 300 N, and other parameters are presented in Table 1. Firstly, the input uncertainties were represented by LHSM theory, and then EPFEM model was constructed. Secondly, the further developed programs were prepared. Then, a Kriging surrogate model was created. In this case, 800 deterministic evaluations are required to establish a reasonable Kriging model. Finally, a large number of specimens (10⁴) were set by LHSM, and the probability analysis was conducted by the surrogate model. Equivalent plastic strain of 0.2% was used as the yield strain^[29], and the plastic zone in the crack tip of TP304 stainless steel is shown in Fig.5, where the contours 1, 2, and 3 indicate the boundaries of plastic zone. Contour 1 indicates the minimum plastic zone when the random parameters are E=199.7 GPa, α=2.773, n=3.087, and P=233.0 N; contour 2 indicates the plastic zone when all the parameters are their mean values; contour 3 indicates the maximum plastic zone when the random parameters are E=202.4 GPa, α=13.226, n=2.892, and P=352.7 N. The results suggest that the impact of random factors on the plastic zone of crack tip is very significant.

Fig.5  Schematic diagram of plastic zones of crack tip (equivalent plastic strain=0.2%)

The Ford-Andresen model^[30] is widely accepted as a rational description of SCC and usually used to obtain the SCC growth rate of structure material in high-temperature oxygenated aqueous water. In this model, the crack growth rate is directly proportional to the crack tip strain rate. Xue et al^[31] found that the variation in tensile plastic strain (${\epsilon }_{\mathrm{p}}^{\mathrm{22}}$) at the characteristic distance r₀ can replace the crack tip strain rate. Therefore, it is necessary to explore the probability distribution of ${\epsilon }_{\mathrm{p}}^{\mathrm{22}}$ at a certain distance from the front-end of crack tip. Fig.6 shows ${\epsilon }_{\mathrm{p}}^{\mathrm{22}}$ values of the front-end of crack tip. It can be observed that the ${\epsilon }_{\mathrm{p}}^{\mathrm{22}}$ is rapidly decreased with increasing the distance from the crack tip. Since the characteristic distance is a critical parameter in the prediction of SCC crack growth rate, collecting the probability results of ${\epsilon }_{\mathrm{p}}^{\mathrm{22}}$ as much as possible is beneficial for engineers to improve the precision prediction of crack growth rate.

Fig.6  Relationship between tensile plastic strains of crack tip and the distance from crack tip

3.2 J-integral in the crack tip

The J-integral indicates a vital crack driving force to describe the mechanical state of the crack tip. Traditionally, the J-integral can be calculated by EPFEM or simplified estimations. For more precise calculation, the deterministic J-integral was calculated with the mean values of parameters in Table 1. A nouveau contour was used because of the numerical accuracy in the mechanical state at crack tip. Fig.7 shows the comparison of J-integral values obtained by EPFEM and the method in Ref.[32] under the load of 100–2000 N. It is observed that the J-integral results of the proposed method are in good agreement with those obtained by the method in Ref.[32]: the J-integral value is increased rapidly with increasing the load.

Fig.7  Comparison of J-integral values of TP304 stainless steel obtained by different methods

Afterwards, the probability prediction was conducted by the proposed method to calculate the maximum value, mean value, and minimum value when the load varies from 100 N to 2000 N. Fig.8 shows the variation of J-integral values under

Fig.8  Maximum, mean, and minimum J-integral values under different loads

different loads. It is observed that the J-integral value is random and rapidly increased with increasing the load. Therefore, the safe margin can be deduced based on the current input parameters.

Since a large number of J-integral results are available by the proposed method, it is possible to obtain the statistical characteristics, including moment, density function, etc. When the load is 300 N, the probability density function (PDF) and cumulative distribution function (CDF) of the J-integral values are shown in Fig.9. It can be seen that the J-integral distribution is not symmetric. The most possible J-integral value varies between 0.8 and 1.4. CDF curve presents the probability that J-integral value is lower than a given threshold. For example, when the J-integral value is smaller than 1 kJ/m², the probability is around 35%.

Fig.9  PDF and CDF curves of J-integral values

For the small uncracked ligament, the crack initiation at flaws can be characterized by the J-integral value which exceeds the material fracture toughness (J_Ic). The number of specimens used in LHSM is 10⁴. The mean J_Ic is set as 1242.6 kJ/m², and COV of J_Ic is set as 0.47^[11]. Failure probability (P_F) was also calculated by the proposed method. The variation of P_F with load is shown in Fig.10. It can be seen that the failure probability is increased with increasing the load. The maximum value of P_F is 0.76 when the load is 2000 N.

Fig.10  Relationship between failure probability and load of 1T-CT specimen

3.3 Global sensitivity analysis of input variables

The purpose of global sensitivity analysis (GSA) is to evaluate the influence of an input parameter on the output variance. Since the parameter importance is mainly investigated, GSA is of high interest in the probability analysis. Based on the surrogate model obtained in Ref.[18], the global sensitivity index of each input parameter is evaluated by the Sobol indices, and Fig.11 shows the sensitivity index results. It can be seen that the load P has a great influence on the J-integral value, presenting the Sobol index of 0.65, which indicates that the dispersion of load P should be reduced firstly to decrease the dispersion of J-integral value. The variables n and α are also important parameters. The most negligible contribution originates from the elastic modulus E, whose Sobol index is lower than 0.1.

Fig.11  Sensitivity index of different parameters

3.4 Efficiency and accuracy of proposed method

Probability prediction for the crack driving force has been extensively researched. However, its application is still restricted due to the complexity of nuclear engineering structures and the non-existence of probability methods in specialized simulation software. Based on EPFEM coupled with Kriging surrogate model, the proposed method was further developed by ABAQUS with MATLAB software. EPFM for complex structures can be analyzed by numerical simulation software, and the Kriging surrogate model can be constructed by many packages, such as the DACE package in MATLAB. Then, the probability prediction can be achieved by a simple operation.

In order to obtain the efficiency and accuracy of the proposed method, the direct MCM, direct LHSM coupled with the computational model, and the proposed method were compared, and Table 3 shows the comparison results. It can be seen that all results obtained by the proposed method are in good agreement with those obtained by the direct MCM and direct LHSM, which indicates that this method can accurately predict the crack driving force and conduct the probability analysis. It is also found that the proposed method is more efficient than the direct MCM and LHSM. Additionally, the Kriging model establishment occupies most operation time rather than the calculation when the proposed method is conducted.

Table 3  Comparison of the proposed procedure with different direct sampling methods

Method	Number of specimens	Time	Mean J-integral value, μJ	Stand deviation of J-integral value, σJ
Direct MCM	105	12.5 d	1.2312	0.3985
Direct LHSM	103	3 h	1.2472	0.4004
Proposed method	103	10 s	1.2358	0.3963

The comparison of J-integral values obtained by the proposed method and MCM is shown in Fig.12. Firstly, 1000 input specimens were obtained according to the distribution of the random parameters in Table 1. The mean load is 300 N. Then, the J-integral value of each input specimen was calculated by the proposed method and MCM. It is found that the J-integral values vary mainly between 0.5 and 2 kJ/m². The slope of the linear fitting line is 0.9834; the root mean square error (R_MSE) and the goodness (R²) of fitting line are 0.00198 and 0.9874, respectively. The results indicate that the proposed method can present the predictions of crack driving force which are similar to those of the computational model by MCM. In conclusion, EPFEM coupled with the Kriging surrogate model is capable of probability prediction of the crack driving force on crack tips of TP304 stainless steel.

Fig.12  Comparison of J-integral values obtained by proposed method and MCM

4 Conclusions

1) A probability prediction for crack driving force through the elastic-plastic finite element method (EPFEM) coupled with the Kriging surrogate model was proposed, which presents good accuracy and high efficiency. Based on the further development of numerical analysis software, many probability parameters, such as failure probability, probability density function, cumulative distribution function, and global sensitivity analysis of the crack driving force can be obtained by the proposed method.

2) The response of deterministic model cannot fully characterize the variation of the crack driving force. Obtaining the probability prediction results as much as possible is beneficial to better understand and accurately predict the mechanical states.

3) The crack driving force of TP304 stainless steel can be influenced by the elastic modulus (E), constant (α), strain hardening exponent (n), and load (P). The load P and strain hardening exponent n have obvious effects on the crack driving force, whereas E has neglectable contribution.

4) The proposed method can predict the crack driving force with good accuracy and high efficiency, compared with the direct sampling methods.