# An Effective Bayesian Method for Probability Fatigue Crack Propagation Modeling through Test Data

Authors: Xu Du; Yuting He; Tianyu Zhang; Chao Gao; Teng Zhang; Sheng Zhang
Journal: Mathematical Problems in Engineering
 (2019)

---

## Abstract

Fatigue crack growth test for 2A12-T4 aluminum alloy was conducted under constant amplitude loading, and the scatter of fatigue crack growth was analyzed by using experimental data based on mathematical statistics. A probabilistic modeling method was introduced to describe the crack growth behavior of 2A12-T4 aluminum alloy. The posterior distribution of model parameter is obtained based on diffuse prior distribution and fatigue crack test data, which is through Bayesian updating. Based on posterior samples of model parameter, the simulation steps and approach give us the crack length exceedance probability, the cumulative distribution function of loading cycle number, and scatter of crack length and loading cycle number, of which simulation results were used to verify the veracity and superiority of the proposed model versus the experimental results. In the present study, it can be used for the reliability assessment of aircraft cracked structures.


---

## Body

## 1. Introduction

Based on aircraft structural life management viewpoint, the mechanical performance of metal materials is often considered homogeneous. However, a considerable amount of scatter has been observed in fatigue life (includes crack initiation life and crack propagation life) even under the same service (loading and environment) conditions.Fatigue crack propagation life prediction is a main work for the damage tolerance design of key structures [1, 2]. Due to the considerable scatter and stochastic of fatigue crack propagation behavior, many researchers have carried out a lot of works and some achievements have been obtained [3–5]. And it has been indicated that the scatter is mainly affected by the dispersity of material’s performance for each single structure [6, 7].As we all know, a lot of fatigue crack growth data are needed to conduct the statistical analysis, and it is a time-consuming work to get these meaningful fatigue crack growth data [8]. To capture the statistical nature of fatigue crack growth, many stochastic probabilistic models have been proposed in literature [9–11]. For the proposed stochastic probabilistic models, the scatter of fatigue crack growth behavior must be taken into account.To continue with the investigation, in this paper, the fatigue crack growth tests are carried out under constant amplitude loading condition. And one meaningful fatigue crack growth dataset is obtained. Based on proposed method and Bayesian updating, the crack length exceedance probability, loading cycle number probability, and scatter of crack growth behavior determination approach are introduced briefly. Some comments on results of the analyses are made at the end of this article.


## 2. Probability Fatigue Crack Growth Model

### 2.1. Fatigue Crack Growth Model

Paris’ law relates the stress intensity factor range to subcritical crack growth under a fatigue stress regime [12, 13]. The basic formula reads(1)dadN=CΔKmwhere a is the crack length (in mm), da/dN is the crack growth rate per fatigue loading cycle,N is the number of fatigue loading cycles, and ΔK is the stress intensity factor range (in MPam) that is related to the applied load, crack length, and specimen/structures geometry. C and m are constants that depend on material, environment, and stress ratio, which can be obtained by testing standard specimens.With the difficulty for the determination of stress intensity factor rangeΔK, the Paris’ law is sometimes modified or simplified [10, 14]. Yang’s randomized Paris model after investigation of crack propagation in fastener holes of aircraft under loading spectra is introduced [15]. Though a lot of works about crack damage prognosis approach were done, the accurate da/dN is difficult to obtain under realistic service conditions.


### 2.2. Probability Modeling for Crack Growth

The crack length afterN number of loading cycles is written as a(N). Based on the Paris’ law, a(N) can be expressed by [13](2)aN=a01-a0θ2θ1θ2N1/θ2where a0 is initial crack length and (θ1, θ2) are the model parameters of single specimen crack growth under constant amplitude loading.The feasibility of using (2) to describe the fatigue crack growth behavior about the metal material has been validated in Michale S. Hamada’s work [13, 16], which will be validated and analyzed in the following work of this study. This equation about a(N) is a semiempirical crack growth fitting model, which has no physical meaning.Equation (2) is a nonlinear function about the loading cycle number N. The real crack propagation curve for every single specimen can be expressed as(3)aiN=a0i1-a0iθ2iθ1iθ2iN1/θ2iwhere subscript index i represents the sequence number of specimen. So a0i is the initial crack length of No.i specimen; θ⃑i=(θ1i,θ2i) is the crack model parameter of No.i specimen.(4)aiN=DN,θ⃑iwhereD(…) represents a function about the loading cycle numberN and random model parameter matrix θ⃑i. So, the scatter of simplified fatigue crack growth model (as (4)) can be described by the distribution of the matrix θ⃑i. And the inverse function of D(N,θ⃑i) expresses as D-1(a,θ⃑i).Based on the simplified model, for a given crack lengthac, the probability function of loading cycles numberN can be expressed as(5)RN=PθDN,θ⃑≥acFor a given loading cycle numberNc, the probability function of crack length a can be expressed as(6)Ra=PθD-1a,θ⃑≤NcBased on (5) and (6), we can know that the reliability function of N and the reliability function of a can be determined by experimental crack propagation curves, which do not take the measure error and model parameter error into account. In other words, the probability distribution can be determined by the probability distribution of model matrix parameter θ⃑=(θ1,θ2). It is difficult to obtain the analytical solution form of (5) and (6), so the simulation approach will be given in the following work of this study.


### 2.3. Prior Distribution

Taking natural e logarithms in (2), we can obtain the transformed observed crack length as(7)ln⁡aNa0=-1θ2ln⁡1-a0θ2θ1θ2NFor the transformed observed crack lengthY, it includes the measure error ε. Yij of every single specimen can be expressed as(8)Yij=-1θ2iln⁡1-a0iθ2iθ1iθ2iNij+εijwhere the subscript index i and j still represent the sequence number of specimen and the sequence number of experimental observation point, respectively. Yij and εij are the observed transformed crack length and observation error for jth observation point of No.i specimen, respectively. The εij subjects to normal distribution (simply as N(0,σ2)) and error εij for different specimen and different observation point are i.i.d. (independent identically distribution).Likelihood of transformed observed crack length includes the probability density function (PDF) ofYij and the PDF of vector matrix parameter θ⃑i=(θ1i,θ2i).(9)Yij~normal-1θ2iln⁡1-a0iθ2iθ1iθ2iNij,σ2θ⃑i=θ1i,θ2i~multi-normalμ⃑,Σ⃑In case of a complete lack of prior information, the diffuse distribution can be used to describe the prior information of parameterσ, μ⃑, and Σ⃑. Assumption matrix parameter μ⃑ subjects to multinormal distribution, expressed as(10)μ⃑~multi-normalμ⃑0,Σ⃑μ0,μ⃑0=00,Σ⃑μ0=1000001000Assumption parameterΣ⃑ is subject to inverse Wishart distribution, expressed as(11)Σ⃑~inverse-WishartΣ⃑0,2,Σ⃑0=100010And assumption parameterσ2 subjects to inverse gamma distribution, expressed as(12)σ2~inverse-gamma3,0.001It is worth noting that the inverse Wishart distribution is a probability distribution defined on real-valued positive-definite matrices [17]. The inverse gamma distribution is a two-parameter family of continuous probability distributions on the positive real line, which is the distribution of the reciprocal of a variable distributed according to the gamma distribution [18]. In Bayesian statistics, inverse Wishart distribution and inverse gamma distribution are commonly used as the conjugate prior for the covariance matrix of a multivariate normal distribution.


## 2.1. Fatigue Crack Growth Model

Paris’ law relates the stress intensity factor range to subcritical crack growth under a fatigue stress regime [12, 13]. The basic formula reads(1)dadN=CΔKmwhere a is the crack length (in mm), da/dN is the crack growth rate per fatigue loading cycle,N is the number of fatigue loading cycles, and ΔK is the stress intensity factor range (in MPam) that is related to the applied load, crack length, and specimen/structures geometry. C and m are constants that depend on material, environment, and stress ratio, which can be obtained by testing standard specimens.With the difficulty for the determination of stress intensity factor rangeΔK, the Paris’ law is sometimes modified or simplified [10, 14]. Yang’s randomized Paris model after investigation of crack propagation in fastener holes of aircraft under loading spectra is introduced [15]. Though a lot of works about crack damage prognosis approach were done, the accurate da/dN is difficult to obtain under realistic service conditions.


## 2.2. Probability Modeling for Crack Growth

The crack length afterN number of loading cycles is written as a(N). Based on the Paris’ law, a(N) can be expressed by [13](2)aN=a01-a0θ2θ1θ2N1/θ2where a0 is initial crack length and (θ1, θ2) are the model parameters of single specimen crack growth under constant amplitude loading.The feasibility of using (2) to describe the fatigue crack growth behavior about the metal material has been validated in Michale S. Hamada’s work [13, 16], which will be validated and analyzed in the following work of this study. This equation about a(N) is a semiempirical crack growth fitting model, which has no physical meaning.Equation (2) is a nonlinear function about the loading cycle number N. The real crack propagation curve for every single specimen can be expressed as(3)aiN=a0i1-a0iθ2iθ1iθ2iN1/θ2iwhere subscript index i represents the sequence number of specimen. So a0i is the initial crack length of No.i specimen; θ⃑i=(θ1i,θ2i) is the crack model parameter of No.i specimen.(4)aiN=DN,θ⃑iwhereD(…) represents a function about the loading cycle numberN and random model parameter matrix θ⃑i. So, the scatter of simplified fatigue crack growth model (as (4)) can be described by the distribution of the matrix θ⃑i. And the inverse function of D(N,θ⃑i) expresses as D-1(a,θ⃑i).Based on the simplified model, for a given crack lengthac, the probability function of loading cycles numberN can be expressed as(5)RN=PθDN,θ⃑≥acFor a given loading cycle numberNc, the probability function of crack length a can be expressed as(6)Ra=PθD-1a,θ⃑≤NcBased on (5) and (6), we can know that the reliability function of N and the reliability function of a can be determined by experimental crack propagation curves, which do not take the measure error and model parameter error into account. In other words, the probability distribution can be determined by the probability distribution of model matrix parameter θ⃑=(θ1,θ2). It is difficult to obtain the analytical solution form of (5) and (6), so the simulation approach will be given in the following work of this study.


## 2.3. Prior Distribution

Taking natural e logarithms in (2), we can obtain the transformed observed crack length as(7)ln⁡aNa0=-1θ2ln⁡1-a0θ2θ1θ2NFor the transformed observed crack lengthY, it includes the measure error ε. Yij of every single specimen can be expressed as(8)Yij=-1θ2iln⁡1-a0iθ2iθ1iθ2iNij+εijwhere the subscript index i and j still represent the sequence number of specimen and the sequence number of experimental observation point, respectively. Yij and εij are the observed transformed crack length and observation error for jth observation point of No.i specimen, respectively. The εij subjects to normal distribution (simply as N(0,σ2)) and error εij for different specimen and different observation point are i.i.d. (independent identically distribution).Likelihood of transformed observed crack length includes the probability density function (PDF) ofYij and the PDF of vector matrix parameter θ⃑i=(θ1i,θ2i).(9)Yij~normal-1θ2iln⁡1-a0iθ2iθ1iθ2iNij,σ2θ⃑i=θ1i,θ2i~multi-normalμ⃑,Σ⃑In case of a complete lack of prior information, the diffuse distribution can be used to describe the prior information of parameterσ, μ⃑, and Σ⃑. Assumption matrix parameter μ⃑ subjects to multinormal distribution, expressed as(10)μ⃑~multi-normalμ⃑0,Σ⃑μ0,μ⃑0=00,Σ⃑μ0=1000001000Assumption parameterΣ⃑ is subject to inverse Wishart distribution, expressed as(11)Σ⃑~inverse-WishartΣ⃑0,2,Σ⃑0=100010And assumption parameterσ2 subjects to inverse gamma distribution, expressed as(12)σ2~inverse-gamma3,0.001It is worth noting that the inverse Wishart distribution is a probability distribution defined on real-valued positive-definite matrices [17]. The inverse gamma distribution is a two-parameter family of continuous probability distributions on the positive real line, which is the distribution of the reciprocal of a variable distributed according to the gamma distribution [18]. In Bayesian statistics, inverse Wishart distribution and inverse gamma distribution are commonly used as the conjugate prior for the covariance matrix of a multivariate normal distribution.


## 3. Fatigue Crack Growth Test

### 3.1. Materials and Specimens

2A12-T4 aluminum alloy is used widely for aircraft structures in China. The mechanical properties and chemical composition of 2A12-T4 are listed in Tables1 and 2, respectively.Table 1
Chemical compositions for 2A12-T4 aluminum alloy (%).
Si Fe Cu Mn Mg Cr Zn Al 0.50 0.50 3.8-4.9 0.30-1.0 1.2-1.8 0.10 0.25 ReTable 2
Mechanical performance for 2A12-T4 aluminum alloy.
Yield strength Tensile strength Young’s modulus Elongation (%) 364 MPa 457 MPa 72 MPa 19.9Single edge cracked specimens used for the fatigue crack growth experiment was shown in Figure1. The width of specimen is 50 mm; the thickness is 2 mm. There was a set of 10 specimens. All specimens were made in L-T plane, with loading aligned in the longitudinal direction. In order to remove the effect of initial crack life, all specimens were made with 8 mm precrack.Figure 1
Shape and dimensions of specimen.


### 3.2. Experimental Procedure

The fatigue crack growth experiments were carried out at room temperature on a MTS material test machine equipped with a PC (personal computer). All experiments were conducted under the same constant amplitude loading cycle (under 20Hz). These experiments were performed under zero-tension condition with stress ratioR=0.06. During the testing, fatigue crack growth data were automatically recorded in the PC, as shown in the Figure 2.Figure 2
Fatigue experiment and crack length measure.


### 3.3. Experimental Results

The crack growth tests are under constant amplitude stress. The maximum stress is 65 MPa and the stress rationR is 0.06. The experimental fatigue crack length is shown in Figure 3.Figure 3
The experimentala-N curves of specimens.Crack growth propagation region can be identified into three regions: near-threshold (∆Kth) region, stable crack propagation region, and unstable (fast) crack propagation region. In order to get more stable data to investigate the crack propagation feature, the data of stable crack propagation region are chosen (crack length: 1.5mm~18mm). The partial fatigue crack growth a-N data of stable propagation region was shown in Figure 4. It is worth noting that the crack length in the Figures 3 and 4 does not include the precrack length (8 mm).Figure 4
The partiala-N curves of specimens.It has been known that, even under the same loading conditions, the number of loading fatigue cycles (at a given crack length) and the crack length (under same loading cycle number) all have a wide range of scatter [19, 20]. A considerable scatter of fatigue crack growth can be seen from the data of Figure 4.It is very necessary to investigate the scatter of fatigue crack growth behavior. The mean valuex-, standard deviation (SD), and coefficients of variation (CV) are used to describe the variation of statistics. The mean value x- can be expressed as (13)x-=1nx1+x2+···+xnwhere n is the total number of experimental data and xi(i=1,…,n) is the interesting data of No.i specimen. In this study, the interesting data is crack length or fatigue cycle number.The SD and CV are expressed as (14) and (15), respectively. It is worth noting that (14) has better applicability for small samples.(14)SD=∑i=1nxi2-nx-2n-1(15)CV=SDx¯The SD and CV of random loading cycle number (to reach a given crack length) for the set specimens were calculated; the calculating results were shown in Figures5 and 6, respectively.Figure 5
The SD of loading cycle under same crack length.Figure 6
The CV of loading cycle under same crack length.As shown in the Figure5, to reach a specific crack length, the SD of loading cycle number for different specimen gradual is increasing with the crack growth, besides the third data point (crack size is 9 mm). Figure 6 shows that the CV of loading cycle number is gradually decreasing with the crack growth, which means that the scatter degree of fatigue loading cycle number decreases as the crack growth.Under specific loading cycle number, the SD and CV of crack length of different specimen were calculated, and the results were shown in Figures7 and 8, respectively.Figure 7
The SD of crack length under same loading cycle.Figure 8
The CV of crack length under same loading cycle.Figures7 and 8 show that, under specific loading cycle number, the SD and CV of crack length increase with the increase of loading cycle number N.Based on the above discussion, the scatter of fatigue crack growth data must be taken into the work of damage tolerance structural design and aircraft structural service life management.


### 3.4. Validating of Probability Model

It is necessary to validate the accuracy of using the probability model to describe crack growth behavior of 2A12-T4 aluminum alloy. The value of matrix parameter (θ1, θ2) can be obtained based on crack growth experimental a-N data through maximum likelihood estimation (MLE) method or other methods. The estimating mean value of parameters (θ1, θ2) for this matrix specimens is listed in Table 3.Table 3
Estimation value of (θ1, θ2) for every single specimen.
Parameter No.1 No.2 No.3 No.4 No.5 θ 1 0.72 0.76 0.87 0.42 0.84 θ 2 0.23 0.31 0.12 0.68 0.11  Parameter No.6 No.7 No.8 No.9 No.10  θ 1 0.50 1.19 1.05 1.10 1.05 θ 2 0.48 -0.04 0.04 0.08 0.14For every single specimen, the fatigue crack growth curve can be obtained by substituting the corresponding estimation value of (θ1, θ2) to (2). In this study, the obtained curve is called crack growth fitting curve. The crack growth fitting curves of every single specimen are shown in Figures 9(a) and 9(b). The sampling rate for plotting Figures 9 and 10 is 1000 fatigue cycles.Figure 9
Crack growth fitting curves of specimens.
(a) (b)Figure 10
The error of crack growth fitting curve.
(a) (b)The error between fitting curve result and experimental observed result is defined as(16)Error=pridictioncurve-experimentcurveexperimentcurve×100%The error of crack growth fitting curve is shown in Figures10(a) and 10(b).As is shown in Figure10, the error between the fitting value and the experimental results was less than 8%. It means that the probabilistic model and the parameter estimating method are satisfactorily in using engineering. As a result, the probabilistic model is a feasible and accurate analytical solution to describe the crack growth behavior of 2A12-T4 aluminum alloy.


## 3.1. Materials and Specimens

2A12-T4 aluminum alloy is used widely for aircraft structures in China. The mechanical properties and chemical composition of 2A12-T4 are listed in Tables1 and 2, respectively.Table 1
Chemical compositions for 2A12-T4 aluminum alloy (%).
Si Fe Cu Mn Mg Cr Zn Al 0.50 0.50 3.8-4.9 0.30-1.0 1.2-1.8 0.10 0.25 ReTable 2
Mechanical performance for 2A12-T4 aluminum alloy.
Yield strength Tensile strength Young’s modulus Elongation (%) 364 MPa 457 MPa 72 MPa 19.9Single edge cracked specimens used for the fatigue crack growth experiment was shown in Figure1. The width of specimen is 50 mm; the thickness is 2 mm. There was a set of 10 specimens. All specimens were made in L-T plane, with loading aligned in the longitudinal direction. In order to remove the effect of initial crack life, all specimens were made with 8 mm precrack.Figure 1
Shape and dimensions of specimen.


## 3.2. Experimental Procedure

The fatigue crack growth experiments were carried out at room temperature on a MTS material test machine equipped with a PC (personal computer). All experiments were conducted under the same constant amplitude loading cycle (under 20Hz). These experiments were performed under zero-tension condition with stress ratioR=0.06. During the testing, fatigue crack growth data were automatically recorded in the PC, as shown in the Figure 2.Figure 2
Fatigue experiment and crack length measure.


## 3.3. Experimental Results

The crack growth tests are under constant amplitude stress. The maximum stress is 65 MPa and the stress rationR is 0.06. The experimental fatigue crack length is shown in Figure 3.Figure 3
The experimentala-N curves of specimens.Crack growth propagation region can be identified into three regions: near-threshold (∆Kth) region, stable crack propagation region, and unstable (fast) crack propagation region. In order to get more stable data to investigate the crack propagation feature, the data of stable crack propagation region are chosen (crack length: 1.5mm~18mm). The partial fatigue crack growth a-N data of stable propagation region was shown in Figure 4. It is worth noting that the crack length in the Figures 3 and 4 does not include the precrack length (8 mm).Figure 4
The partiala-N curves of specimens.It has been known that, even under the same loading conditions, the number of loading fatigue cycles (at a given crack length) and the crack length (under same loading cycle number) all have a wide range of scatter [19, 20]. A considerable scatter of fatigue crack growth can be seen from the data of Figure 4.It is very necessary to investigate the scatter of fatigue crack growth behavior. The mean valuex-, standard deviation (SD), and coefficients of variation (CV) are used to describe the variation of statistics. The mean value x- can be expressed as (13)x-=1nx1+x2+···+xnwhere n is the total number of experimental data and xi(i=1,…,n) is the interesting data of No.i specimen. In this study, the interesting data is crack length or fatigue cycle number.The SD and CV are expressed as (14) and (15), respectively. It is worth noting that (14) has better applicability for small samples.(14)SD=∑i=1nxi2-nx-2n-1(15)CV=SDx¯The SD and CV of random loading cycle number (to reach a given crack length) for the set specimens were calculated; the calculating results were shown in Figures5 and 6, respectively.Figure 5
The SD of loading cycle under same crack length.Figure 6
The CV of loading cycle under same crack length.As shown in the Figure5, to reach a specific crack length, the SD of loading cycle number for different specimen gradual is increasing with the crack growth, besides the third data point (crack size is 9 mm). Figure 6 shows that the CV of loading cycle number is gradually decreasing with the crack growth, which means that the scatter degree of fatigue loading cycle number decreases as the crack growth.Under specific loading cycle number, the SD and CV of crack length of different specimen were calculated, and the results were shown in Figures7 and 8, respectively.Figure 7
The SD of crack length under same loading cycle.Figure 8
The CV of crack length under same loading cycle.Figures7 and 8 show that, under specific loading cycle number, the SD and CV of crack length increase with the increase of loading cycle number N.Based on the above discussion, the scatter of fatigue crack growth data must be taken into the work of damage tolerance structural design and aircraft structural service life management.


## 3.4. Validating of Probability Model

It is necessary to validate the accuracy of using the probability model to describe crack growth behavior of 2A12-T4 aluminum alloy. The value of matrix parameter (θ1, θ2) can be obtained based on crack growth experimental a-N data through maximum likelihood estimation (MLE) method or other methods. The estimating mean value of parameters (θ1, θ2) for this matrix specimens is listed in Table 3.Table 3
Estimation value of (θ1, θ2) for every single specimen.
Parameter No.1 No.2 No.3 No.4 No.5 θ 1 0.72 0.76 0.87 0.42 0.84 θ 2 0.23 0.31 0.12 0.68 0.11  Parameter No.6 No.7 No.8 No.9 No.10  θ 1 0.50 1.19 1.05 1.10 1.05 θ 2 0.48 -0.04 0.04 0.08 0.14For every single specimen, the fatigue crack growth curve can be obtained by substituting the corresponding estimation value of (θ1, θ2) to (2). In this study, the obtained curve is called crack growth fitting curve. The crack growth fitting curves of every single specimen are shown in Figures 9(a) and 9(b). The sampling rate for plotting Figures 9 and 10 is 1000 fatigue cycles.Figure 9
Crack growth fitting curves of specimens.
(a) (b)Figure 10
The error of crack growth fitting curve.
(a) (b)The error between fitting curve result and experimental observed result is defined as(16)Error=pridictioncurve-experimentcurveexperimentcurve×100%The error of crack growth fitting curve is shown in Figures10(a) and 10(b).As is shown in Figure10, the error between the fitting value and the experimental results was less than 8%. It means that the probabilistic model and the parameter estimating method are satisfactorily in using engineering. As a result, the probabilistic model is a feasible and accurate analytical solution to describe the crack growth behavior of 2A12-T4 aluminum alloy.


## 4. Posterior Distribution and Numerical Simulation

### 4.1. Posterior Distribution

Theμ1, μ2, Σ11, Σ12, Σ21, Σ22, and σ are the distribution parameters of the parameters θ1, θ2, and ε. The posterior distribution of those parameters can be obtained by using Bayesian updating. Briefly, OpenBUGS is a software package for Bayesian analysis of complex multiparameter statistical models with the Gibbs sampling based, which is a MCMC method. For more information about background, application, and introduction to OpenBUGS refer to Amiri [21] and Modarres [22] and Cowles. The posterior distribution for those distribution parameters of fatigue crack growth model is listed in Table 4, which are obtained by updating experimental crack growth length data.Table 4
Posterior distribution of parametersμ⃑, Σ⃑, and σ.
Parameter Mean SD MC-error (%) Different percentile 0.025 0.5 0.975 μ 1 0.846 0.086 0.061 0.675 0.846 1.017 μ 2 0.214 0.071 0.058 0.074 0.214 0.359 Σ 11 0.070 0.043 0.034 0.026 0.06 0.177 Σ 12 -0.053 0.034 0.031 -0.137 -0.045 -0.017 Σ 21 -0.053 0.034 0.031 -0.137 -0.045 -0.017 Σ 22 0.048 0.03 0.032 0.016 0.040 0.124 σ 0.033 0.004 0.005 0.027 0.032 0.041Calculating the Monte Carlo error (MC-error) of posterior distribution for each interested parameter is an effective way to assess the accuracy of estimation processing. As a rule of thumb, the simulation should be run until the Monte Carlo error for each interested parameter is less than about 5% of the sample standard deviation. As listed in Table4, the MC-error for all interested parameters is less than 5%.Convergence diagnosis is an important work for the using of MCMC. Contrasting deviation method is a common method to diagnose convergences. Based on the contrasting deviation method proposed by Gelman and Rubin [23], a more briefness and convenience one was established by Brooks and Gelman [24]. The basic idea and approach can reference to literature [25]. In this study, the convergent of interested parameters is good, and the trend of model parameter with iteration increase is convergent.


### 4.2. Random Cycle Number Distribution

Shapiro-Wilk test gives a good approach to assess the fitting degree of normal distribution. [26]. Shapiro-Wilk test tests the null hypothesis that a sample x=(x1,x2,…,xn) came from a normally distributed population. The test statistic W is(17)W=∑k=1nskxk2∑k=1nxk-x-2where x(k) (with parentheses enclosing the subscript index k) is the kth order statistic, i.e., the kth-smallest number in the sample x. x-=(x1+x2+⋯+xn)/n is the sample mean and constant vector matrix s given by (18)s=s1,s2,⋯,sn=qTV-1qTV-1V-1q1/2where q=(q1, q2, …, qn)T and q1, q2, …, qn are the expected values of the order statistics of independent and identically distributed random normal distribution and V is the covariance matrix of those order statistics.The variableW and variable p-value are the key variables in Shapiro-Wilk test. And the R software provides a common method for the calculating of W and p-value. In this study, the loading cycle number data for specific crack length can be obtained by tests. Before the calculation of variable W and variable p-value, the loading cycle number N was taken in logarithmic base 10. The calculating value of variable W and variable p-value is listed in Table 5.Table 5
W and p-value of loading cycle lg-number.
Crack length (mm) 3 6 9 12 15 18 W 0.93 0.95 0.89 0.87 0.85 0.86  p-value 0.394 0.643 0.161 0.109 0.052 0.073As a rule of thumb, if the calculationp-value is less than the chosen alpha level, then the null hypothesis is rejected and there is an evidence that the data tested are not from a normally distributed population; in other words, the data are not normal. On the contrary, the data came from a normally distributed population cannot be rejected.In the airplane design, in order to insure aircraft structural service safety, the requirement significance levelγ is 0.05. As listed in Table 5, the calculation p-value is all greater than 0.05. Based on the above analysis, for the observed crack length, the loading cycles number N subjects to lg-normal distribution.For a given crack length, the lg-normal distribution parameter (mean and standard deviation) of loading cycles numberN can be obtained by MLE or other methods, and the cumulative distribution function (CDF) can be obtained. For example, when the given crack length is 3mm, the cumulative probability distribution of experimental loading cycle number is shown in Figure 11(a).Figure 11
Cycle number distribution to reach specified crack length: (a) 3mm, (b) 6mm, (c) 9mm, (d) 12mm, (e) 15mm, and (f) 18mm.
(a) (b) (c) (d) (e) (f)When the given crack length is 3mm, the cumulative distribution of random loading cycle number can be obtained by the following simulation steps:(1) Draw a sample of matrix distribution parameter(μ⃑,Σ⃑) form the posterior samples (generated by MCMC) as the distribution parameter of θ⃑=(θ1,θ2), written as (μ⃑,Σ⃑)l.(2) Based on the multinormal distribution(μ⃑,Σ⃑)l, generate 1000 random samples as the value of parameter (θ1, θ2).(3) Based on the above 1000 simulation random (θ1, θ2) samples, estimate the probability value of R(N), which satisfies a-1(ac,θ⃑)≤N for different loading cycle number. For this case, ac is the given crack length, in this condition ac=3mm.(4) The approximate estimation value ofR(N) is obtained by cycling the steps (1)-(3) 10000 numbers.The approximateR(N) is the CDF of loading cycle number. When crack grows to 3mm, the analytical CDF of loading cycle number can be obtained on the above frameworks, which is shown in Figure 11(a). So, the comparison diagram (Figure 11(a)) intuitionisticly reflects the goodness of fit between experimental CDF and analytical CDF of loading cycle number N.The experimental and analytical comparison for the observed crack length (6mm, 9mm, 12mm, 15mm, and 18mm) can be obtained through the same analysis approach. Figures11(b), 11(c), 11(d), 11(e), and 11(f) corresponding to the given crack length is 6mm, 9mm, 12mm, 15mm, and 18mm, respectively.As is shown in Figure11, for the observed crack length, the experimental CDF of loading cycle number fits the analytical CDF very well. It means that the modified probabilistic fatigue crack growth model and above simulation approach can describe the crack growth behavior pretty well.It is worth noting that the distribution of random cycle number for the other crack lengths can be obtained by the above simulation approach.


### 4.3. Probability Exceedance of Crack Length

Under observed loading cycle number, the crack length data for the specimen can be obtained by measurement.Inspection distribution type is the firstly work to determination probability of crack exceedance under specific loading cycle number. Shapiro-Wilk test is still used to the distribution type inspection of crack length. The calculation value of statistic variablesW and p-value is listed in Table 6.Table 6
W and p-value of crack length.
Cycle number (104×N) 3 6 9 12 15 18 W 0.94 0.95 0.95 0.93 0.95 0.92  P-value 0.535 0.662 0.626 0.470 0.667 0.371As is listed in Table6, when the requirement significance level γ is 0.05, the estimation p-value is all greater than γ=0.05. Thus, under the observed loading cycle number, the random crack length data subjects to normal distribution. The distribution parameter of experimental crack length can be obtained by MLE method or other methods. So the experimental probability of crack exceedance for the observed loading cycle numbers can be obtained.When the given loading cycle number is 30000, the probability exceedance of crack length as well as probability function of loading cycle number can be obtained by the following simulation steps:(1) Draw a random sample of matrix parameter(μ⃑,Σ⃑) form the posterior samples as the distribution parameter value of θ⃑=(θ1,θ2), writing as (μ⃑,Σ⃑)l.(2) Based on the multinormal distribution(μ⃑,Σ⃑)l, generate 1000 random samples as the parameter (θ1, θ2) value.(3) Based on the 1000 random model parameter (θ1, θ2) value, estimate the value of R(a) which satisfies a(Nc,θ⃑)≥a for different crack length a, where the Nc represents the given crack length.(4) The approximate estimation value ofR(a) can be obtained by cycling steps (1)–(3) 10000 numbers.The approximateR(a) is a probability of crack exceedance. Thus, the analytical probability of crack exceedance for loading cycle 30000 number can be obtained, which is shown in Figure 12(a). The comparison diagram (Figure 12(a)) intuitionisticly reflects the goodness of fit between experimental probability exceedance and analytical probability exceedance.Figure 12
Probability of crack l exceedance for loading cycles (104N) is (a) 3, (b) 6, (c) 9, (d) 12, (e) 15, and (f) 18.
(a) (b) (c) (d) (e) (f)Figures12(b), 12(c), 12(d), 12(e), and 12(f) corresponding to the loading cycle number are 60000, 90000, 120000, 150000, and 180000, respectively.As is shown in Figure12, the experimental probability of crack exceedance fits the analytical probability exceedance very well. It means that the modified fatigue crack growth model and above simulation approach can predicate the actual crack growth behavior pretty well.It is worth noting that the probability exceedance of crack length for the other loading cycle numberN can obtain through the above simulation approach by giving the loading cycle number.


### 4.4. Analytical Scatter Analysis

The validity of the introduced analytical solution was substantiated by convergence diagnose of the interested parameter. Accuracy of the probabilistic crack growth model was substantiated by designing simulation steps and approach, whose simulation results were compared with the experimental results. However, the scatter analysis is an important part for the validation of this study.The scatter of loading cycle numberN (to reach a given crack length) and the scatter of crack length (under same loading cycle number) can be obtained by the following simulation steps:(1) Draw a random sample of matrix distribution parameter(μ⃑,Σ⃑) form the posterior samples (generated by MCMC) as the distribution parameter value of θ⃑=(θ1,θ2), written as (μ⃑,Σ⃑)l.(2) Based on the multinormal distribution(μ⃑,Σ⃑)l, generate 1000 random samples as the value of model parameter (θ1, θ2).(3) Based on the every (θ1, θ2) simulation random samples, calculate a(N) for different loading cycle number N (0<N≤180000, and simulation step is 1000). Under the different loading cycle number N, calculate the mean of a(N) (1000 samples). And the mean a(N) is corresponding to (μ⃑,Σ⃑)l.(4) Based on the every (θ1,θ2) simulation random samples, calculate N(a) for different loading cycle number a (0<a≤18mm, and simulation step is 0.1mm). Under the different crack length a, calculate the mean of N(a) (1000 samples). And the mean N(a) is corresponding to (μ⃑,Σ⃑)l.(5) Cycling steps (1)-(4) 10000 numbers, we can obtain 10000a(N) and N(a).(6) Calculatethe mean, SD, and CV of above 10000a(N) and N(a).The analytical scatter of loading cycle number and crack length can be obtained. Combining with experimental scatter calculation for the two cases, the analytical and experimental SD and CV of crack length under same loading cycle number are shown in Figure13. The analytical and experimental SD and CV of loading cycle number to reach a specified crack length are shown in Figure 14.Figure 13
Analytical statistic variable value of crack length: (a) SD and (b) CV.
(a) (b)Figure 14
Analytical statistic variable value of loading cycle number: (a) SD and (b) CV.
(a) (b)As is shown in Figures13 and 14, the analytical scatter value is in good agreement with the experimental value. And the validity and accuracy of proposed crack growth model were further substantiated.It is founded that the introduced probabilistic fatigue crack growth modeling method and simulation approach can predict the crack growth behavior pretty well. The proposed simulation framework provides a good approach to obtain the CDF of loading cycle number, the probability of crack exceedance, and the scatter of crack growth behavior.


## 4.1. Posterior Distribution

Theμ1, μ2, Σ11, Σ12, Σ21, Σ22, and σ are the distribution parameters of the parameters θ1, θ2, and ε. The posterior distribution of those parameters can be obtained by using Bayesian updating. Briefly, OpenBUGS is a software package for Bayesian analysis of complex multiparameter statistical models with the Gibbs sampling based, which is a MCMC method. For more information about background, application, and introduction to OpenBUGS refer to Amiri [21] and Modarres [22] and Cowles. The posterior distribution for those distribution parameters of fatigue crack growth model is listed in Table 4, which are obtained by updating experimental crack growth length data.Table 4
Posterior distribution of parametersμ⃑, Σ⃑, and σ.
Parameter Mean SD MC-error (%) Different percentile 0.025 0.5 0.975 μ 1 0.846 0.086 0.061 0.675 0.846 1.017 μ 2 0.214 0.071 0.058 0.074 0.214 0.359 Σ 11 0.070 0.043 0.034 0.026 0.06 0.177 Σ 12 -0.053 0.034 0.031 -0.137 -0.045 -0.017 Σ 21 -0.053 0.034 0.031 -0.137 -0.045 -0.017 Σ 22 0.048 0.03 0.032 0.016 0.040 0.124 σ 0.033 0.004 0.005 0.027 0.032 0.041Calculating the Monte Carlo error (MC-error) of posterior distribution for each interested parameter is an effective way to assess the accuracy of estimation processing. As a rule of thumb, the simulation should be run until the Monte Carlo error for each interested parameter is less than about 5% of the sample standard deviation. As listed in Table4, the MC-error for all interested parameters is less than 5%.Convergence diagnosis is an important work for the using of MCMC. Contrasting deviation method is a common method to diagnose convergences. Based on the contrasting deviation method proposed by Gelman and Rubin [23], a more briefness and convenience one was established by Brooks and Gelman [24]. The basic idea and approach can reference to literature [25]. In this study, the convergent of interested parameters is good, and the trend of model parameter with iteration increase is convergent.


## 4.2. Random Cycle Number Distribution

Shapiro-Wilk test gives a good approach to assess the fitting degree of normal distribution. [26]. Shapiro-Wilk test tests the null hypothesis that a sample x=(x1,x2,…,xn) came from a normally distributed population. The test statistic W is(17)W=∑k=1nskxk2∑k=1nxk-x-2where x(k) (with parentheses enclosing the subscript index k) is the kth order statistic, i.e., the kth-smallest number in the sample x. x-=(x1+x2+⋯+xn)/n is the sample mean and constant vector matrix s given by (18)s=s1,s2,⋯,sn=qTV-1qTV-1V-1q1/2where q=(q1, q2, …, qn)T and q1, q2, …, qn are the expected values of the order statistics of independent and identically distributed random normal distribution and V is the covariance matrix of those order statistics.The variableW and variable p-value are the key variables in Shapiro-Wilk test. And the R software provides a common method for the calculating of W and p-value. In this study, the loading cycle number data for specific crack length can be obtained by tests. Before the calculation of variable W and variable p-value, the loading cycle number N was taken in logarithmic base 10. The calculating value of variable W and variable p-value is listed in Table 5.Table 5
W and p-value of loading cycle lg-number.
Crack length (mm) 3 6 9 12 15 18 W 0.93 0.95 0.89 0.87 0.85 0.86  p-value 0.394 0.643 0.161 0.109 0.052 0.073As a rule of thumb, if the calculationp-value is less than the chosen alpha level, then the null hypothesis is rejected and there is an evidence that the data tested are not from a normally distributed population; in other words, the data are not normal. On the contrary, the data came from a normally distributed population cannot be rejected.In the airplane design, in order to insure aircraft structural service safety, the requirement significance levelγ is 0.05. As listed in Table 5, the calculation p-value is all greater than 0.05. Based on the above analysis, for the observed crack length, the loading cycles number N subjects to lg-normal distribution.For a given crack length, the lg-normal distribution parameter (mean and standard deviation) of loading cycles numberN can be obtained by MLE or other methods, and the cumulative distribution function (CDF) can be obtained. For example, when the given crack length is 3mm, the cumulative probability distribution of experimental loading cycle number is shown in Figure 11(a).Figure 11
Cycle number distribution to reach specified crack length: (a) 3mm, (b) 6mm, (c) 9mm, (d) 12mm, (e) 15mm, and (f) 18mm.
(a) (b) (c) (d) (e) (f)When the given crack length is 3mm, the cumulative distribution of random loading cycle number can be obtained by the following simulation steps:(1) Draw a sample of matrix distribution parameter(μ⃑,Σ⃑) form the posterior samples (generated by MCMC) as the distribution parameter of θ⃑=(θ1,θ2), written as (μ⃑,Σ⃑)l.(2) Based on the multinormal distribution(μ⃑,Σ⃑)l, generate 1000 random samples as the value of parameter (θ1, θ2).(3) Based on the above 1000 simulation random (θ1, θ2) samples, estimate the probability value of R(N), which satisfies a-1(ac,θ⃑)≤N for different loading cycle number. For this case, ac is the given crack length, in this condition ac=3mm.(4) The approximate estimation value ofR(N) is obtained by cycling the steps (1)-(3) 10000 numbers.The approximateR(N) is the CDF of loading cycle number. When crack grows to 3mm, the analytical CDF of loading cycle number can be obtained on the above frameworks, which is shown in Figure 11(a). So, the comparison diagram (Figure 11(a)) intuitionisticly reflects the goodness of fit between experimental CDF and analytical CDF of loading cycle number N.The experimental and analytical comparison for the observed crack length (6mm, 9mm, 12mm, 15mm, and 18mm) can be obtained through the same analysis approach. Figures11(b), 11(c), 11(d), 11(e), and 11(f) corresponding to the given crack length is 6mm, 9mm, 12mm, 15mm, and 18mm, respectively.As is shown in Figure11, for the observed crack length, the experimental CDF of loading cycle number fits the analytical CDF very well. It means that the modified probabilistic fatigue crack growth model and above simulation approach can describe the crack growth behavior pretty well.It is worth noting that the distribution of random cycle number for the other crack lengths can be obtained by the above simulation approach.


## 4.3. Probability Exceedance of Crack Length

Under observed loading cycle number, the crack length data for the specimen can be obtained by measurement.Inspection distribution type is the firstly work to determination probability of crack exceedance under specific loading cycle number. Shapiro-Wilk test is still used to the distribution type inspection of crack length. The calculation value of statistic variablesW and p-value is listed in Table 6.Table 6
W and p-value of crack length.
Cycle number (104×N) 3 6 9 12 15 18 W 0.94 0.95 0.95 0.93 0.95 0.92  P-value 0.535 0.662 0.626 0.470 0.667 0.371As is listed in Table6, when the requirement significance level γ is 0.05, the estimation p-value is all greater than γ=0.05. Thus, under the observed loading cycle number, the random crack length data subjects to normal distribution. The distribution parameter of experimental crack length can be obtained by MLE method or other methods. So the experimental probability of crack exceedance for the observed loading cycle numbers can be obtained.When the given loading cycle number is 30000, the probability exceedance of crack length as well as probability function of loading cycle number can be obtained by the following simulation steps:(1) Draw a random sample of matrix parameter(μ⃑,Σ⃑) form the posterior samples as the distribution parameter value of θ⃑=(θ1,θ2), writing as (μ⃑,Σ⃑)l.(2) Based on the multinormal distribution(μ⃑,Σ⃑)l, generate 1000 random samples as the parameter (θ1, θ2) value.(3) Based on the 1000 random model parameter (θ1, θ2) value, estimate the value of R(a) which satisfies a(Nc,θ⃑)≥a for different crack length a, where the Nc represents the given crack length.(4) The approximate estimation value ofR(a) can be obtained by cycling steps (1)–(3) 10000 numbers.The approximateR(a) is a probability of crack exceedance. Thus, the analytical probability of crack exceedance for loading cycle 30000 number can be obtained, which is shown in Figure 12(a). The comparison diagram (Figure 12(a)) intuitionisticly reflects the goodness of fit between experimental probability exceedance and analytical probability exceedance.Figure 12
Probability of crack l exceedance for loading cycles (104N) is (a) 3, (b) 6, (c) 9, (d) 12, (e) 15, and (f) 18.
(a) (b) (c) (d) (e) (f)Figures12(b), 12(c), 12(d), 12(e), and 12(f) corresponding to the loading cycle number are 60000, 90000, 120000, 150000, and 180000, respectively.As is shown in Figure12, the experimental probability of crack exceedance fits the analytical probability exceedance very well. It means that the modified fatigue crack growth model and above simulation approach can predicate the actual crack growth behavior pretty well.It is worth noting that the probability exceedance of crack length for the other loading cycle numberN can obtain through the above simulation approach by giving the loading cycle number.


## 4.4. Analytical Scatter Analysis

The validity of the introduced analytical solution was substantiated by convergence diagnose of the interested parameter. Accuracy of the probabilistic crack growth model was substantiated by designing simulation steps and approach, whose simulation results were compared with the experimental results. However, the scatter analysis is an important part for the validation of this study.The scatter of loading cycle numberN (to reach a given crack length) and the scatter of crack length (under same loading cycle number) can be obtained by the following simulation steps:(1) Draw a random sample of matrix distribution parameter(μ⃑,Σ⃑) form the posterior samples (generated by MCMC) as the distribution parameter value of θ⃑=(θ1,θ2), written as (μ⃑,Σ⃑)l.(2) Based on the multinormal distribution(μ⃑,Σ⃑)l, generate 1000 random samples as the value of model parameter (θ1, θ2).(3) Based on the every (θ1, θ2) simulation random samples, calculate a(N) for different loading cycle number N (0<N≤180000, and simulation step is 1000). Under the different loading cycle number N, calculate the mean of a(N) (1000 samples). And the mean a(N) is corresponding to (μ⃑,Σ⃑)l.(4) Based on the every (θ1,θ2) simulation random samples, calculate N(a) for different loading cycle number a (0<a≤18mm, and simulation step is 0.1mm). Under the different crack length a, calculate the mean of N(a) (1000 samples). And the mean N(a) is corresponding to (μ⃑,Σ⃑)l.(5) Cycling steps (1)-(4) 10000 numbers, we can obtain 10000a(N) and N(a).(6) Calculatethe mean, SD, and CV of above 10000a(N) and N(a).The analytical scatter of loading cycle number and crack length can be obtained. Combining with experimental scatter calculation for the two cases, the analytical and experimental SD and CV of crack length under same loading cycle number are shown in Figure13. The analytical and experimental SD and CV of loading cycle number to reach a specified crack length are shown in Figure 14.Figure 13
Analytical statistic variable value of crack length: (a) SD and (b) CV.
(a) (b)Figure 14
Analytical statistic variable value of loading cycle number: (a) SD and (b) CV.
(a) (b)As is shown in Figures13 and 14, the analytical scatter value is in good agreement with the experimental value. And the validity and accuracy of proposed crack growth model were further substantiated.It is founded that the introduced probabilistic fatigue crack growth modeling method and simulation approach can predict the crack growth behavior pretty well. The proposed simulation framework provides a good approach to obtain the CDF of loading cycle number, the probability of crack exceedance, and the scatter of crack growth behavior.


## 5. Conclusion

The fatigue crack growth behavior of 2A12-T4 aluminum alloy is modeling with an introduced probabilistic model with Bayesian updating. We can get the following meaningful conclusions:(1) The introduced probabilistic model can describe the fatigue crack growth behavior of 2A12-T4 aluminum alloy accurately.(2) The scatter of crack growth behavior is described by introducing multinormal distribution of parameter (θ1, θ2), which can take into full consideration the scatter feature of 2A12-T4 aluminum alloy (with thickness 2mm) through crack growth behavior.(3) This study proposed a probabilistic crack growth modeling method to determine the probability exceedance and the scatter of random crack length, as well as the CDF and scatter of random loading cycle numberN which can be determined.



---

*Source: 7843627-2019-03-03.xml*