Probabilistic fatigue damage prognosis using ... - Yongming Liu
Probabilistic fatigue damage prognosis using maximum entropy approach Xuefei Guan, Ratneshwar Jha Department of Mechanical & Aerospace Engineering, Clarkson University, Potsdam, NY 136995725 USA Yongming Liu* Department of Civil & Environmental Engineering, Clarkson University, Potsdam, NY 136995710 USA Abstract: A general framework for probabilistic fatigue damage prognosis using maximum entropy concept is proposed in this paper. The fatigue damage is calculated using a physics-based crack growth model. Due to the stochastic nature of crack growth process, uncertainties arising from the underlying physical model, parameters of the model, and the measurement noise are considered and integrated into this framework. A maximum relative entropy (MRE) approach is proposed to update the prognosis result and confidence bounds incorporating various types of uncertainties from measuring, modeling, and maintenance. Markov Chain Monte Carlo (MCMC) sampling is employed to generate the posterior probability distribution of model parameters and provide statistical information for the maximum relative entropy updating procedure. Numerical examples are used to demonstrate the proposed MRE prognosis methodology. Experimental data for aluminum alloys are used to validate model predictions under uncertainty. Following this, a detailed comparison between the proposed MRE approach and classical Bayesian updating method is performed to illustrate advantages of the proposed MRE approach.
Keywords: Maximum relative entropy, Fatigue crack growth, Confidence intervals, Markov Chain Monte Carlo, Life prediction, Model updating.
1. Introduction Fatigue damage is critical for the reliability evaluation and health management of many mechanical/structural components. Experimental data indicate that the change in the crack growth rate is not a smooth, stable, well ordered process (Virkler et al., 1979). Therefore, fatigue damage accumulation is stochastic in nature. Several physical models have been proposed to describe this process and are used to predict the crack growth curve, among which the Paris type of crack growth laws (Paris and Erdogan, 1963; Foreman et al., 1967; Walker, 1970) are commonly used. The parameters in those models are generally determined using experimental 1 * Corresponding author, Tel.: 315-268-2341; Fax: 315-268-7985; Email:[email protected]
data by regression analysis or empirical judgments. However, various uncertainties of the prediction from the underlying physical model, parameter estimation and future inspection and maintenance need to be carefully included for risk assessment and decision-making. In order to describe the stochastic process of fatigue crack growth, measurements during this process and other empirical information may be used to calibrate the model under uncertainty. Several studies have been reported to combine the information obtained from inspection with the physical model to update the crack growth and fatigue life prediction. Bayes’ theorem is widely used for probability updating in engineering. Madsen developed the idea of updating probability using non-destructive inspection within the Bayesian framework (Madsen, 1997). Zhang and Mahadevan (2000) extended this idea to incorporate the uncertainties of mechanical model, probabilistic model, and distribution statistics into a Bayesian updating procedure. Perrin et al. (2007) applied the Markov Chain Monte Carlo approach together with Bayesian updating for fracture analysis Entropy methods, such as Maximum Entropy (MaxEnt) (Jaynes, 1957; Jaynes, 1979; Skilling, 1988) principle and relative entropy methods, are alternative approaches for assigning and updating probabilities when new information is given in the form of a constraint. Caticha and Giffin reported that the maximum relative entropy (MRE) approach reproduces every aspect of Bayesian updating rules and Jeffrey’s updating rule is also a special case of MRE method (Caticha and Giffin, 2006). The maximum relative entropy approach can address problems that cannot be handled by either the MaxEnt or classical Bayesian methods individually (Giffin and Caticha, 2007). For example, updating can be performed given the data in form of expected values constraints. MRE approach has been used in physics for statistical mechanics (Caticha and Preuss, 2004; Tseng and Caticha, 2008). The objective of this paper is to develop a general maximum relative entropy approach for model updating and prediction of fatigue damage. The underlying physics-based crack growth model is used to describe the deterministic damage accumulation. Some model parameters are assumed to follow a statistical distribution based on measurement data. Maximum relative entropy theory is used to update the model parameters and prognosis results using additional information from inspection or health monitoring system. A cascade Markov Chain Monte Carlo (MCMC) sampling technique is used to generate the posterior probability distribution. The statistics of the samples is then integrated into the posterior as a moment 2
constraint in the maximum relative entropy framework. A real application case shows that maximum relative entropy approach yields a more accurate confidence interval of crack size prediction compared to the classical Bayesian approach.
2. Maximum Relative Entropy for model updating The basic principle of the maximum relative entropy is introduced in this section and simple numerical example is presented for probabilistic residual life prediction using MRE approach. 2.1. MRE principles Given two normalized probability density functions (PDF) f1 (θ ) and f 2 (θ ) , the relative information entropy I ( f1 : f 2 ) of f1 with respect to f 2 , also referred to as the Kullback-Leibler
divergence (Kullback and Leibler, 1951), is defined by
∫
I ( f 1 : f 2 ) = − dθ ⋅ f 1 (θ ) log Θ
f 1 (θ ) f 2 (θ )
(1)
It is shown that the relative information entropy is zero only if f1 (θ ) = f 2 (θ ) . The “axioms” of maximum entropy indicates that the form of Eq. (1) is the unique entropy to be used in inductive inference. The three axioms (Skilling, 1988) are: Axiom 1: Locality - Local information has local effects. Axiom 2: Coordinative invariance - The ranking of the two probability densities should not depend on the system coordinates. This indicates that the coordinates carry no information. Axiom 3: Consistency for independent subsystem - For a system composed of subsystems that are believed to be independent, it should not make a difference whether the inference treats them separately or jointly. For an inverse problem, the posterior of a quantity θ ∈ Θ is inferred on the basis of three pieces of information: prior information about θ (the prior PDF), the known relationship between x and θ (the mathematical/physical model), and the observed values of the data x ∈ Χ . The search space for posterior of θ is Χ × Θ . Let µ (x, θ ) be the prior joint distribution and p( x, θ ) be the posterior joint distribution. 3
According to the maximum relative entropy axioms, the selected joint posterior is the one that maximizes the relative entropy I ( p : µ ) of p( x, θ ) with respect to the prior µ ( x, θ ) , subject to the constraints. I ( p : µ ) = − ∫ dxdθ ⋅ p(x, θ ) log
p ( x, θ ) µ ( x, θ )
(2)
In Eq. (2), µ (x, θ ) = µ (θ )µ (x | θ ) contains all prior information. µ ( x | θ ) is the likelihood function and µ (θ ) is the prior distribution of θ . When new information is available in the form of constraints, the updating procedure will search in the space of Χ × Θ for a posterior which maximizes I ( p : µ ) . Any information can be represented as a constraint can be used to update the posterior. For example, the observation of response variable is one type of constraint that can be used in Bayesian and the maximum relative entropy updating. However, other types of constraints, usually in the form of mean values or statistical moments of an interested parameter cannot be easily used in Bayesian updating, but they can be used in maximum relative entropy updating. Following the derivation of the maximum relative entropy posterior (Caticha and Giffin, 2006), if a new observation data point x′ is obtained, the maximum relative entropy updating will treat this information as a constraint on the posterior. The posteriors that reflect the fact x is now known to be x′ is a constraint such that c1 : p(x ) = ∫ dθ ⋅ p( x, θ ) = δ ( x − x′)
(3)
Other information in the form of a moment constraint, such as the expected value of some function g (θ ) , can be represented as c2 : ∫ dxdθ ⋅ p(x, θ ) ⋅ g (θ ) = g (θ ) = G
(4)
The normalization constraint can be expressed as
∫
c 3 : dxdθ ⋅ p(x, θ ) = 1
(5)
Maximizing Eq. (2) subject to the above constraints Eqs. (3-5),
[
] [
]
I + α dxdθ ⋅ p(x, θ ) − 1 + β dxdθ ⋅ p(x, θ ) ⋅ g (θ ) − G ∫ ∫ =0 L + ∫ dx ⋅ λ (x ) ∫ dθ ⋅ p(x, θ ) − δ (x − x′)
δ
[
]
(6)
the posterior can be obtained as p (θ ) ∝ µ (θ )µ (x′ | θ )e β ⋅ g (θ )
4
(7)
The normalization of Eq. (7) yields p(θ ) =
µ (θ )µ (x′ | θ )e β ⋅ g (θ ) Z (x′, β )
(8)
where Z (x′, β ) = ∫ dθ ⋅ µ (θ )µ (x′ | θ )e β ⋅ g (θ ) is the normalization constant. The Lagrange multiplier β is determined by ∂ ln Z ( x′, β ) =G ∂β
(9)
The right side of Eq. (7) consists of three terms. µ (θ ) is the parameter prior probability distribution, µ ( x′ | θ ) is the likelihood function, and e β ⋅ g (θ ) is the exponential term introduced by moment constraints. Eq. (7) is similar to Bayesian posterior except for the additional exponential term for moment constraints. This equation indicates that, if no moment constraint is available (i.e., β is zero), MRE updating will be identical to Bayesian updating procedure. In other words, Bayesian updating is a special case of MRE updating given no information on the moment constraints.
2.2. A numerical example of MRE approach The above updating procedure and the relationship between MRE and Bayesian updating can be illustrated using a numerical example. Given that a system component residual life follows an exponential distribution, µ (t | θ ) = θe −θ t for (t > 0 ) . The probability distribution of parameter θ is a random variable following a gamma distribution µ (θ ) = θe −θ for (θ > 0 ) . The prior joint distribution is
µ (t , θ ) = µ (t | θ )µ (θ ) = θ 2e −θ (t +1) The plot of prior distribution of θ is shown in Fig. 1.
5
(10) Formatted: Font: Not Bold Deleted: Fig. 1
Prior probability distribution of parameter θ 1
PDF(θ)
0.8
0.6
0.4
0.2
0
0
1
2
3
4
5
6
θ
Fig. 1 The prior probability distribution Different updating scenarios are discussed here and the results are compared with each other.
Case 1: Updating with observed data Suppose one measurement of the component life, t ′ = 1 , is available, MRE updating procedure for θ can be performed according to Eq. (8) with β = 0 . The updated posterior is p(θ ) = µ (θ )
µ (t ′ | θ ) = 4θ 2e − 2θ µ (t ′)
(11)
This is the same result with Bayesian rule because no moment constraints are given.
Case 2: Updating with moment information It is assumed that the expected value of θ ( θ = G = 0.8 ) is obtained by other experimental measurements and no other information is given. The maximum relative entropy updating procedure can be performed using Eq. (8), omitting the likelihood term µ ( x′ | θ ) . In this case, g (θ ) = θ . β is determined to be -1.5 by Eq. (9). The normalization yields a posterior given by Eq. (12).
p(θ ) = 6.25θ e −2.5θ
(12)
Case 3: Updating with observed data and moment information simultaneously Combining one piece of observed data t ′ = 1 and the moment information θ = G = 0.8 ,
6
the normalized posterior is updated to be
p (θ ) = 26.3671875θ 2e −3.75θ
(13)
with β = −1.75 . The prior and posteriors updated with different constraints are plotted in Fig. 2. The figure shows that moment information adds additional constraints on posteriors. Compared with Bayesian rules characterized by posterior updated with observed data, MRE updating procedure generates different posteriors in this example. The different posteriors give different predictions in practice. Prior and posteriors of parameter θ Prior Posterior (updating with t'=1)
1
Posterior (updating with 〈θ〉 =0.8) Posterior (updating with t'=1 and 〈θ〉 =0.8)
PDF(θ)
0.8
0.6
0.4
0.2
0
0
1
2
3
4
5
6
θ
Fig. 2 The prior and posteriors updated with different constraints
3. MRE updating for fatigue prognosis The fatigue prognosis is important for structure health management and residual life prediction. From the previous sections, a general form of the entropy posterior for is known to be in the form of Eq. (7). In order to use MRE updating for fatigue prognosis, a physical model is required to deterministically describe the fatigue crack growth. The prior statistical distributions of model parameters are required to incorporate the parameter estimation uncertainties and a likelihood function needs to be constructed using the physical model to characterize the model and measurement uncertainties. Detailed derivation and implementation of MRE updating to fatigue prognosis are shown below.
7
Deleted: Fig. 2
3.1. Deterministic crack growth model The deterministic model is used to describe the underlying physics of the fatigue damage accumulation. It is noted that the proposed MRE updating can be applied to different physical models depending on the specific problem. For demonstration and validation purpose, Paris law for the fatigue crack growth is chosen as the deterministic model (Paris and Erdogan, 1963). It is widely used to describe crack growth rate curve under constant amplitude cyclic loading:
da m = c(∆K ) dN
(14)
In Eq. (14), a is the crack length, ∆K is the stress intensity factor range, and (c, m ) are material parameters. For a plate specimen with a center through crack, ∆K is defined as:
a ∆K = πa ∆σF w
(15)
where w is the specimen width and F (⋅) is geometry correction factor for the specimen. ∆σ is the applied stress range. Given c , m and cycle number N , the crack size a can be calculated by solving the ODE equation numerically.
da m = c(∆K ) dN
πa ∆σ =c a F w
m
(16)
3.2. Building the likelihood equation Consider a general model prediction equation, d = y +τ
(17)
where d is the observed data, y is the model prediction, and τ is the general error term, including the measurement error ε and the modeling error e . The model prediction is given by Eq. (18) as y = M (x | θ ) + e
(18)
where M ( x | θ ) is the deterministic model and e is the model prediction error term, x is the prediction value and θ is the deterministic model parameter vector. Equation (17) can be rewritten as 8
d = M (x | θ ) + e + ε
(19)
In general, the probability distributions of the two uncorrelated error terms e and ε can be described using two Gaussian distributions with standard deviations σ e and σ ε , respectively. Let d be the measure of the response variable. The probability distributions of the two error terms
are
p(e | θ ) =
[d − M ( x | θ )]2 1 exp− 2σ e2 2π σ e
(20)
ε2 exp − 2 2π σ ε 2σ ε
(21)
and 1
p(ε | θ ) =
The likelihood function can be constructed as
L (D | θ ) = p (d | θ ) =
[d − M (x | θ )]2 exp− 2 σ e2 + σ ε2 2π σ e2 + σ ε2 1
(
)
(22)
Using general error term τ ~ Gaussian (0, σ τ ) , Eq. (22) can be rewritten as [d − M (x | θ )]2 1 exp− 2σ τ2 2π σ τ
L (D | θ ) = p (d | θ ) =
(23)
For multiple observed data, the likelihood is the product of the form L ( DN | θ ) =
n
∏
p(d i | θ ) =
i =1
1
(
2π σ τ
)
n
n [d − M (xi | θ )]2 exp − i 2σ τ2 i =1
∑
(24)
Substituting Eq. (24) in Eq. (7) and omitting the constant terms, the posterior is derived to be
p (θ ) ∝ µ (θ ) ⋅
n [d i − M (xi | θ )]2 β ⋅g (θ ) exp − ⋅e σ τn 2σ τ2 i =1 1
∑
(25)
3.3. MCMC sampling The analytical posterior probability distribution of c and m in fatigue prognosis is difficult to evaluate because the normalization of Eq. (25) involves multi-dimensional integration. In order to circumvent the difficulty of direct evaluation of the posterior probability distribution described by Eq. (25), a cascade MCMC sampling technique is employed in this paper. MCMC simulation was first introduced by Metropolis et al. (1953) as a method to 9
simulate a discrete-time homogeneous Markov chain. The random walk algorithm of MCMC ensures the samples converge to the posterior probability distribution from an arbitrary initial value. The widely used algorithm of Metropolis-Hastings (Hastings, 1970) is summarized below. Detailed derivation can be found in the referenced article. The transition between two successive samples xt and xt +1 is defined as x ~ q ( X | xt ) with probabilit y α (xt , ~ x) ~ xt +1 = x else t
(26)
where q( X | xt ) is the transition distribution. The acceptance probability α ( xt , ~ x ) is
α (xt , ~x ) = min (1, r )
(27)
The Metropolis ratio is given by r=
p (~ x ) q( xt | ~ x) p ( xt ) q(~ x | xt )
(28)
where p (⋅) is the posterior probability representation. For a symmetric transition distribution
q(⋅) , for example the normal distribution, the property of
q( xt | ~ x ) = q (~ x | xt ) simplifies the
Metropolis ratio to be
r=
p (~ x) p( xt )
(29)
The merit of MCMC is that it overcomes the normalization in Bayesian and the proposed MRE updating procedure. In this paper, a cascade MCMC (Tarantola, 2005) is employed to generate the posterior samples of Eq. (25).
4. Demonstration example and model validation A demonstration example of fatigue crack growth prognosis is illustrated in this section. Fatigue crack growth measurement data are used to demonstrate the application of the proposed maximum relative entropy updating procedure. The accuracy and precision of the prediction results are compared with experimental observations.
10
4.1. Crack growth data A large fatigue crack growth database of 2024-T3 aluminum alloy was reported by Virkler et al. (Virkler et al., 1979). The data set consists of 68 sample trajectories, each containing 164 measurement points. The entire specimen has the same geometry, i.e., an initial crack size ai = 9mm , length L = 558.8mm , width w = 152.4mm and thickness d = 2.54mm . The stress range during each experiment is constant ∆σ = 48.28MPa , and the stress ratio is R = 0.2 . The failure criterion is that the crack size equals 49.8mm . Kotulski reported the statistics of the Paris model parameters (c, m ) . The mean and coefficient of variance (COV) of m are 2.874 and 0.057, respectively. The mean and COV of log(c) are -26.155 and 0.037, respectively (Kotulski, 1998). The geometry correction factor for Virkler’s experiments is
a F = w
1
a cos π w
, for
a < 0.7 w
(30)
Based on Kotulski’s report, the prior PDF of c can be expressed as,
(log(c ) − ζ c )2 exp− 2σ c2 2π σ cc 1
µ (c ) =
(31)
and the prior of the parameter m is
µ (m ) =
(m − ζ m )2 exp− 2σ m2 2π σ m 1
(32)
with ζ c = −26.155 , σ c = COVc ⋅ ζ c = 0.037 ⋅ 26.155 = 0.968 and ζ m = 2.874 , σ m = COVm ⋅ ζ m = 0.057 ⋅ 2.874 = 0.164
4.2. MRE updating with MCMC Substituting the prior µ (θ ) in Eq. (25) with Eq. (31) and Eq. (32), the final form of maximum entropy posterior for the fatigue crack growth problem can be expressed as
11
1 log(c ) − ζ 2 1 1 m − ζ 2 1 c m + β c g c (c ) ⋅ + β m g m (m ) p(c, m ) ∝ exp− exp− σ cc σc 2 σ m 2 σ m 1 n d − M (a | c, m, N ) 2 i i ⋅ n exp− ∑ i στ 2 σ i =1 τ
(33)
1
where M (ai | c, m, N i ) is the Paris deterministic model, d i is the measure of crack size corresponding to cycles N i , and ai is the model prediction of crack size with parameters (c, m ) at N i cycles. The total error term σ τ is set to be 0.1mm . The exponential terms are g c (c ) = log(c ) and g m (m ) = m . The maximum relative entropy updating procedure using the cascade MCMC is described in the following steps. When one observation of d i , N i is available, the proposed MRE updating procedure begins with initial values of β c = β m = 0 . Step 1: Generate statistically sufficient number of (c, m ) samples by the cascade MCMC. Step 2: Calculate crack size a at any interested cycle count N using MCMC samples generated in step 1. This step is for crack size prediction purpose. Step 3: Calculate the coefficient β c , β m in the exponential term with the statistics of MCMC samples for the next updating. If additional observed data become available, the above mentioned steps are repeated till the final failure. One crack growth trajectory in Virkler’s dataset is selected for crack growth prediction
Deleted: arbitrarily
updating for demonstration purpose from previous publication (Ostergaard and Hillberry, 1983). Six points are randomly chosen to simulate the measured values of crack length a obtained by Deleted: Table 1
monitoring systems. The selected observations are shown in Table 1. Table 1 Data used for updating (Virkler’s dataset) Number
1
2
3
4
5
6
Crack size (a )
9.8273
10.5830
11.4121
12.1658
14.5744
20.1360
Cycle (N )
23988
44227
59220
73088
105697
153673
One hundred thousand samples with 5% burn-in period are generated by MCMC sampling in 12
step 1. 99.9% confidence intervals of predicted crack size a are calculated. The prognosis results using different measurements are shown in Fig. 3. It is observed that the confidence intervals of crack length prediction become narrower as more and more measurement data become available for MRE updating. It is also shown that the early stage measurement data can be used to improve model prediction significantly which is essential for condition-based maintenance and unit replacement. (a) Updating with data 1
(b) Updating with data 1-2
50
50 99.9% interval prediction Median prediction Experimental data Updating points
45
45 40
35
35
Crack length (mm)
Crack length (mm)
40
30 25 20
30 25 20
15
15
10
10
5
0
0.5
99.9% interval prediction Median prediction Experimental data Updating points
1 1.5 Number of cycles
5
2
0
0.5
5
x 10
(c) Updating with data 1-3
5
x 10
50 99.9% interval prediction Median prediction Experimental data Updating points
45 40
40
35 30 25 20
35 30 25 20
15
15
10
10
0
0.5
1 1.5 Number of cycles
99.9% interval prediction Median prediction Experimental data Updating points
45
Crack length (mm)
Crack length (mm)
2
(d) Updating with data 1-4
50
5
1 1.5 Number of cycles
5
2 5
x 10
13
0
0.5
1 1.5 Number of cycles
2 5
x 10
(e) Updating with data 1-5
(f) Updating with data 1-6
50
50 99.9% interval prediction Median prediction Experimental data Updating points
45
45 40
35
Crack length (mm)
Crack length (mm)
40
30 25 20
35 30 25 20
15
15
10
10
5
0
0.5
99.9% interval prediction Median prediction Experimental data Updating points
1 1.5 Number of cycles
5
2
0
0.5
5
x 10
1 1.5 Number of cycles
2 5
x 10
Fig. 3 Prognosis results using different measurements
5. Comparison with Bayesian approach Bayesian approach is widely used for reliability updating of engineering systems. A detailed comparison between the proposed MRE updating and the classical Bayesian updating is given in this section to show the advantages of the proposed MRE approach. Considering a Bayesian posterior pB (c, m ) for the same application example: pB (c, m ) ∝
1 log(c ) − ζ 2 1 1 m − ζ 2 1 c m ⋅ exp− exp− σ cc σ σ σ 2 2 c m m 1 n d − M (a | c, m, N ) 2 i i ⋅ n exp− ∑ i στ στ 2 i =1
(34)
1
Compared with MRE posterior in Eq. (33), the Bayesian posterior has no exponential terms introduced by moment constraints. In order to give a distinct comparison between Bayesian and MRE updating, all the parameters for prior and MCMC sampling in Bayesian updating procedure are the same with the ones used in the proposed MRE updating. Besides the dataset collected by Virkler’s experiments, another dataset acquired by McMaster on aluminum alloy T351 is also employed for extensive comparison between these two methodologies.
Formatted Formatted: heading2, Left, Level 1, Widow/Orphan control
5.1 Virkler’s dataset
All the physical and experimental description about the samples are referred to the previous section. The predictions of crack size using Bayesian approach and MRE approach are plotted in 14
Deleted: Fig. 4
Fig. 4. (a) Updating with data 1
(b) Updating with data 1-2
50
50 99.9% interval prediction (Bayasian) Median prediction (Bayesian) 99.9% interval prediction (MrE) Median prediction (MrE) Experimental data Updating points
Crack length (mm)
40 35
45 40
Crack length (mm)
45
30 25 20
35 30 25 20
15
15
10
10
5
0
0.5
1 1.5 Number of cycles
99.9% interval prediction (Bayasian) Median prediction (Bayesian) 99.9% interval prediction (MrE) Median prediction (MrE) Experimental data Updating points
5
2
0
0.5
5
x 10
(c) Updating with data 1-3
35
40
30 25 20
35 30 25 20
15
15
10
10
0
0.5
1 1.5 Number of cycles
99.9% interval prediction (Bayasian) Median prediction (Bayesian) 99.9% interval prediction (MrE) Median prediction (MrE) Experimental data Updating points
45
Crack length (mm)
Crack length (mm)
40
5
2
0
0.5
5
x 10
(e) Updating with data 1-5
1 1.5 Number of cycles
2 5
x 10
(f) Updating with data 1-6
50
50 99.9% interval prediction (Bayasian) Median prediction (Bayesian) 99.9% interval prediction (MrE) Median prediction (MrE) Experimental data Updating points
40 35
40
30 25 20
35 30 25 20
15
15
10
10
0
0.5
1 1.5 Number of cycles
99.9% interval prediction (Bayasian) Median prediction (Bayesian) 99.9% interval prediction (MrE) Median prediction (MrE) Experimental data Updating points
45
Crack length (mm)
45
Crack length (mm)
5
x 10
50 99.9% interval prediction (Bayasian) Median prediction (Bayesian) 99.9% interval prediction (MrE) Median prediction (MrE) Experimental data Updating points
45
5
2
(d) Updating with data 1-4
50
5
1 1.5 Number of cycles
5
2 5
x 10
0
0.5
1 1.5 Number of cycles
2 5
x 10
Fig. 4 Comparison of prognosis results between Bayesian updating and MRE updating (Virkler’s dataset) The comparison between MRE updating and the classical Bayesian updating in the application example shows that at the very beginning stage of crack growth, in particular the 15
prediction generated by first two measurements of response shown in Fig. 4 (a) and (b), both
Deleted: Fig. 4
MRE and Bayesian updating have similar predictions. However, as the crack grows to the early stage, characterized by response points 3 - 5 shown in Fig. 4 (a - c), MRE gives more satisfactory
Deleted: Fig. 4
prediction of confidence intervals of crack size, which is significant for residual life prediction. The narrower confidence intervals are due to the additional moment constraints introduced by MCMC samples. As shown in Fig. 4 (f), MRE still gives a narrower confidence interval
Deleted: Fig. 4
prediction compared with Bayesian but the difference is not as prominent as is in the previous updatings in Fig. 4 (c - e). One possible explanation is that the last measurement used in the
Deleted: Fig. 4
updating procedure is near the middle of the crack growth curve, and all the data points are sufficient enough to dominate the posteriors. It is clear that if more and more data become available, this crack curve can be statistically determined using MRE updating procedure.
Formatted Formatted: heading2, Left, Level 1, Indent: First line: 0", Widow/Orphan control
5.2 McMaster’s dataset
McMaster and Simth (1999) collected a set of 2024-T351 aluminum alloy experimental data under constant and variable loading conditions. The experimental data of center-cracked specimens with length L = 250mm , width w = 100mm and thickness t = 6mm under constant loading with ∆σ = 65.7 MPa and stress ratio R = 0.1 are used here for comparison using Paris’ law. The details of experimental setup and conditions can be found in referenced articles. Six points are picked up to simulate the measures of response in Table 2. One hundred thousand samples
Formatted: Lowered by 3 pt Formatted: Lowered by 3 pt Formatted: Lowered by 3 pt Formatted: Lowered by 3 pt Deleted: Table 2
are generated for both Bayesian and MRE updatings. Prediction curves at sixth updating together with a prior prediction are given in Fig. 5 in order to concentrate the results on the comparison. Table 2 Data used for updating (McMaster’s dataset) Number
1
2
3
4
5
6
Crack size (a )
11.3611
11.9282
12.3254
13.8563
14.8771
19.5841
Cycle (N )
4875
8474
11550
17775
21375
34500
Deleted: Fig. 5 Formatted: Check spelling and grammar
Formatted: Justified
Despite of an inaccurate prior curve presented in Fig. 5, MRE updating exhibits narrower intervals than that of Bayesian because of the moment constraints added by MCMC samples in MRE framework. In statistical point of view, as more data are available for such a probabilistic inference. Posterior PDFs given by these two methods will finally converge to the actual target 16
Deleted: Fig. 5 Formatted: Check spelling and grammar Formatted: Justified, Indent: First line: 0.5"
distribution. However, in practice, particularly for large and complex system, the measurements of response variables are expensive and hard to acquire and statistically sufficient amount of data are virtually impossible. At the same time, time-based preventive maintenance might cost too much for a complex system. In this case, it is beneficial to use MRE updating procedure to generate narrower confidence intervals such that early stage warnings can be generated based on those predictions to facilitate condition-based maintenance of the system. Formatted: Keep with next 99.9% interval prediction (Bayesian) Median prediction (Bayesian) 99.9% interval prediction (MRE) Median prediction (MRE) Experimental data Prior estimation Updating points
35
Crack length (mm)
30
25
20
15
10 0
0.5
1
1.5
2 2.5 3 Number of cycles
3.5
4
4.5
5 4
x 10
Fig. 5 Comparion of preciton between Bayesian and MRE updating (McMaster’s dataset)
Formatted: Centered, Level 1, Indent: First line: 0" Deleted: 5
6. Conclusion
Field Code Changed
A general maximum relative entropy approach is developed in this paper to include all available information and uncertainties for updating and prediction. For the unstable and stochastic fatigue crack growth process, the uncertainties include model prediction uncertainty, measurement uncertainty, and model parameter uncertainty. The maximum relative entropy updating procedure based on the measures of response quantities is presented in detail. The posterior samples are generated using a cascade MCMC random walk algorithm. Beyond the classical Bayesian updating procedure, the statistics of MCMC samples is employed to add additional constraints on the posterior. Real application examples reveal that the posterior model
Deleted: A r Deleted: s
obtained by maximum relative entropy approach significantly improves the prediction. Furthermore, comparisons of crack growth prediction between maximum relative entropy and Bayesian updating have been performed to show the benefits of maximum relative entropy 17
Deleted: a c Deleted: s
approach. The comparison shows that MRE updating in general gives more satisfactory prediction. It further indicates the potential that the maximum relative entropy approach can be applied to engineering systems as an improved alternative to the classical Bayesian approach.
Acknowledgement The research reported in this paper was supported by funds from National Aeronautics and Space Administration (NASA) (Contract No. NNX09AY54A). The support is gratefully acknowledged.
References Kullback, S., Leibler, R.A. (1951). On Information and Sufficiency. The Annals of Mathematical Statistics, 22(1):7986. Metropolis, N., Rosenbluth, A.W., Rosenbluth, M.N., Teller, A.H., Teller, E. (1953). Equations of state calculations by fast computing machine. Journal of Chemical Physics, 21(6):1087-1092. Jaynes, E.T. (1957). Information theory and statistical mechanics. Physical Review, 106: 620- 630 and 108: 171190. Paris, P.C., Erdogan, F. (1963). A critical analysis of crack propagation laws. Journal of Basic Engineering, Transactions of American Society of Mechanical Engineers, 85:528-534.
Foreman, R.G., Kearney, V.E., Engle, R.M. (1967). Numerical analysis of crack propagation in a cyclic-loaded structure. Journal of Basic Engineering, Transactions of American Society of Mechanical Engineers, 89D: 459-464. Walker, K. (1970). The effect of stress ratio during crack propagation and fatigue for 2024-T3 and 7075-T6 aluminum. In Effects of environment and complex load history on fatigue life. ASTM STP 462, (pp. 1-14). American Society for Testing and Materials. Hastings, W.K. (1970). Monte Carlo Sampling Methods Using Markov Chains and Their Applications, Biometrika, 57(1):97-109. Jaynes, E.T. (1979). Where do we stand on Maximum Entropy. In R.D. Levine and M. Tribus (Eds.), The Maximum Entropy Formalism (pp. 15-119). M.I.T. Press, Cambridge, MA.
Virkler, D.A., Hillberry, B.M., Goel, P.K. (1979). The statistical nature of fatigue crack propagation, ASME Journal of Engineering Materials and Technology, 101:148-153.
Ostergaard, D.F., Hillberry, B.M. (1983). Characterization of the variability in fatigue crack propagation data, In J.M. Bloom and J.C. Ekvall (Eds.), Probabilistic Fracture Mechanics and Fatigue Methods: Application for Structural Design and Maintenance, ASTM STP 798 (pp. 97-115). American Society for Testing and Materials.
Skilling, J. (1988). The axioms of maximum entropy. In G.J. Erickson and C.R. Smith (Eds.), Maximum Entropy and Bayesian Methods, in Science and Engineering, volume 1, (pp. 173-187). Kluwer Academic Publishers, Netherlands.
18
Madsen, H. (1997). Stochastic modeling of fatigue crack growth and inspection. In C.G. Soares (Ed.), Probabilistic methods for structure design, collection Solid Mechanics and its applications, (pp. 59-83). Kluwer Academic
Publishers, Netherlands. Kotulski, Z.A. (1998). On Efficiency of Identification of a Stochastic Crack Propagation Model based on Virkler Experimental Data. Archives of Mechanics, 50(5):829-847. McMaster, F.J., Smith, D.J. (1999). Effect of load excursions and specimen thickness on crack closure measurements, In R.C. McClung and J.C. Newman, Jr., (Eds) Advances in Fatigue Crack Closure Measurements and Analysis: Second Volume, ASTM STP 1343, (pp. 246-264). West Conshohocken, Pennsylvania: American
Society for Testing and Materials. Zhang, R., Mahadevan, S. (2000). Model uncertainty and Bayesian updating in reliability-based inspection. Structure Safety, 22(2):145-160.
Caticha, A., Preuss, R. (2004). Maximum entropy and Bayesian data analysis: Entropic prior distributions. Physical Review E. 70(4), 046127.
Tarantola, A. (2005). Inverse Problem Theory and Methods for Model Parameter Estimation. Society for Industrial and Applied Mathematics, 3600 University City Science Center, Philadelphia, USA. Caticha, A. & Giffin, A. (2006). Updating probabilities. In A.M. Djafari (Ed.), Bayesian Inference and Maximum Entropy Methods in Science and Engineering, AIP Conference Proceedings, 31: 872.
Perrin, F., Sudret, B., Pendola, M. (2007). Bayesian updating of mechanical models – Application in fracture mechanics. In 18ème Congrès Français de Mécanique, CFM 2007, Grenoble. Giffin, A. & Caticha, A. (2007) Updating probabilities with data and moments. In K. Knuth (Ed.), Bayesian Inference and Maximum Entropy Methods in Science and Engineering, AIP Conference Proceedings, 954:74
Tseng, C., Caticha, A. (2008). Using relative entropy to find optimal approximations: An application to simple fluids. Physica A, 387:6759-6770.
19
Keywords: Maximum relative entropy, Fatigue crack growth, Confidence intervals, Markov Chain Monte Carlo, Life prediction, Model updating.
1. Introduction Fatigue damage is critical for the reliability evaluation and health management of many mechanical/structural components. Experimental data indicate that the change in the crack growth rate is not a smooth, stable, well ordered process (Virkler et al., 1979). Therefore, fatigue damage accumulation is stochastic in nature. Several physical models have been proposed to describe this process and are used to predict the crack growth curve, among which the Paris type of crack growth laws (Paris and Erdogan, 1963; Foreman et al., 1967; Walker, 1970) are commonly used. The parameters in those models are generally determined using experimental 1 * Corresponding author, Tel.: 315-268-2341; Fax: 315-268-7985; Email:[email protected]
data by regression analysis or empirical judgments. However, various uncertainties of the prediction from the underlying physical model, parameter estimation and future inspection and maintenance need to be carefully included for risk assessment and decision-making. In order to describe the stochastic process of fatigue crack growth, measurements during this process and other empirical information may be used to calibrate the model under uncertainty. Several studies have been reported to combine the information obtained from inspection with the physical model to update the crack growth and fatigue life prediction. Bayes’ theorem is widely used for probability updating in engineering. Madsen developed the idea of updating probability using non-destructive inspection within the Bayesian framework (Madsen, 1997). Zhang and Mahadevan (2000) extended this idea to incorporate the uncertainties of mechanical model, probabilistic model, and distribution statistics into a Bayesian updating procedure. Perrin et al. (2007) applied the Markov Chain Monte Carlo approach together with Bayesian updating for fracture analysis Entropy methods, such as Maximum Entropy (MaxEnt) (Jaynes, 1957; Jaynes, 1979; Skilling, 1988) principle and relative entropy methods, are alternative approaches for assigning and updating probabilities when new information is given in the form of a constraint. Caticha and Giffin reported that the maximum relative entropy (MRE) approach reproduces every aspect of Bayesian updating rules and Jeffrey’s updating rule is also a special case of MRE method (Caticha and Giffin, 2006). The maximum relative entropy approach can address problems that cannot be handled by either the MaxEnt or classical Bayesian methods individually (Giffin and Caticha, 2007). For example, updating can be performed given the data in form of expected values constraints. MRE approach has been used in physics for statistical mechanics (Caticha and Preuss, 2004; Tseng and Caticha, 2008). The objective of this paper is to develop a general maximum relative entropy approach for model updating and prediction of fatigue damage. The underlying physics-based crack growth model is used to describe the deterministic damage accumulation. Some model parameters are assumed to follow a statistical distribution based on measurement data. Maximum relative entropy theory is used to update the model parameters and prognosis results using additional information from inspection or health monitoring system. A cascade Markov Chain Monte Carlo (MCMC) sampling technique is used to generate the posterior probability distribution. The statistics of the samples is then integrated into the posterior as a moment 2
constraint in the maximum relative entropy framework. A real application case shows that maximum relative entropy approach yields a more accurate confidence interval of crack size prediction compared to the classical Bayesian approach.
2. Maximum Relative Entropy for model updating The basic principle of the maximum relative entropy is introduced in this section and simple numerical example is presented for probabilistic residual life prediction using MRE approach. 2.1. MRE principles Given two normalized probability density functions (PDF) f1 (θ ) and f 2 (θ ) , the relative information entropy I ( f1 : f 2 ) of f1 with respect to f 2 , also referred to as the Kullback-Leibler
divergence (Kullback and Leibler, 1951), is defined by
∫
I ( f 1 : f 2 ) = − dθ ⋅ f 1 (θ ) log Θ
f 1 (θ ) f 2 (θ )
(1)
It is shown that the relative information entropy is zero only if f1 (θ ) = f 2 (θ ) . The “axioms” of maximum entropy indicates that the form of Eq. (1) is the unique entropy to be used in inductive inference. The three axioms (Skilling, 1988) are: Axiom 1: Locality - Local information has local effects. Axiom 2: Coordinative invariance - The ranking of the two probability densities should not depend on the system coordinates. This indicates that the coordinates carry no information. Axiom 3: Consistency for independent subsystem - For a system composed of subsystems that are believed to be independent, it should not make a difference whether the inference treats them separately or jointly. For an inverse problem, the posterior of a quantity θ ∈ Θ is inferred on the basis of three pieces of information: prior information about θ (the prior PDF), the known relationship between x and θ (the mathematical/physical model), and the observed values of the data x ∈ Χ . The search space for posterior of θ is Χ × Θ . Let µ (x, θ ) be the prior joint distribution and p( x, θ ) be the posterior joint distribution. 3
According to the maximum relative entropy axioms, the selected joint posterior is the one that maximizes the relative entropy I ( p : µ ) of p( x, θ ) with respect to the prior µ ( x, θ ) , subject to the constraints. I ( p : µ ) = − ∫ dxdθ ⋅ p(x, θ ) log
p ( x, θ ) µ ( x, θ )
(2)
In Eq. (2), µ (x, θ ) = µ (θ )µ (x | θ ) contains all prior information. µ ( x | θ ) is the likelihood function and µ (θ ) is the prior distribution of θ . When new information is available in the form of constraints, the updating procedure will search in the space of Χ × Θ for a posterior which maximizes I ( p : µ ) . Any information can be represented as a constraint can be used to update the posterior. For example, the observation of response variable is one type of constraint that can be used in Bayesian and the maximum relative entropy updating. However, other types of constraints, usually in the form of mean values or statistical moments of an interested parameter cannot be easily used in Bayesian updating, but they can be used in maximum relative entropy updating. Following the derivation of the maximum relative entropy posterior (Caticha and Giffin, 2006), if a new observation data point x′ is obtained, the maximum relative entropy updating will treat this information as a constraint on the posterior. The posteriors that reflect the fact x is now known to be x′ is a constraint such that c1 : p(x ) = ∫ dθ ⋅ p( x, θ ) = δ ( x − x′)
(3)
Other information in the form of a moment constraint, such as the expected value of some function g (θ ) , can be represented as c2 : ∫ dxdθ ⋅ p(x, θ ) ⋅ g (θ ) = g (θ ) = G
(4)
The normalization constraint can be expressed as
∫
c 3 : dxdθ ⋅ p(x, θ ) = 1
(5)
Maximizing Eq. (2) subject to the above constraints Eqs. (3-5),
[
] [
]
I + α dxdθ ⋅ p(x, θ ) − 1 + β dxdθ ⋅ p(x, θ ) ⋅ g (θ ) − G ∫ ∫ =0 L + ∫ dx ⋅ λ (x ) ∫ dθ ⋅ p(x, θ ) − δ (x − x′)
δ
[
]
(6)
the posterior can be obtained as p (θ ) ∝ µ (θ )µ (x′ | θ )e β ⋅ g (θ )
4
(7)
The normalization of Eq. (7) yields p(θ ) =
µ (θ )µ (x′ | θ )e β ⋅ g (θ ) Z (x′, β )
(8)
where Z (x′, β ) = ∫ dθ ⋅ µ (θ )µ (x′ | θ )e β ⋅ g (θ ) is the normalization constant. The Lagrange multiplier β is determined by ∂ ln Z ( x′, β ) =G ∂β
(9)
The right side of Eq. (7) consists of three terms. µ (θ ) is the parameter prior probability distribution, µ ( x′ | θ ) is the likelihood function, and e β ⋅ g (θ ) is the exponential term introduced by moment constraints. Eq. (7) is similar to Bayesian posterior except for the additional exponential term for moment constraints. This equation indicates that, if no moment constraint is available (i.e., β is zero), MRE updating will be identical to Bayesian updating procedure. In other words, Bayesian updating is a special case of MRE updating given no information on the moment constraints.
2.2. A numerical example of MRE approach The above updating procedure and the relationship between MRE and Bayesian updating can be illustrated using a numerical example. Given that a system component residual life follows an exponential distribution, µ (t | θ ) = θe −θ t for (t > 0 ) . The probability distribution of parameter θ is a random variable following a gamma distribution µ (θ ) = θe −θ for (θ > 0 ) . The prior joint distribution is
µ (t , θ ) = µ (t | θ )µ (θ ) = θ 2e −θ (t +1) The plot of prior distribution of θ is shown in Fig. 1.
5
(10) Formatted: Font: Not Bold Deleted: Fig. 1
Prior probability distribution of parameter θ 1
PDF(θ)
0.8
0.6
0.4
0.2
0
0
1
2
3
4
5
6
θ
Fig. 1 The prior probability distribution Different updating scenarios are discussed here and the results are compared with each other.
Case 1: Updating with observed data Suppose one measurement of the component life, t ′ = 1 , is available, MRE updating procedure for θ can be performed according to Eq. (8) with β = 0 . The updated posterior is p(θ ) = µ (θ )
µ (t ′ | θ ) = 4θ 2e − 2θ µ (t ′)
(11)
This is the same result with Bayesian rule because no moment constraints are given.
Case 2: Updating with moment information It is assumed that the expected value of θ ( θ = G = 0.8 ) is obtained by other experimental measurements and no other information is given. The maximum relative entropy updating procedure can be performed using Eq. (8), omitting the likelihood term µ ( x′ | θ ) . In this case, g (θ ) = θ . β is determined to be -1.5 by Eq. (9). The normalization yields a posterior given by Eq. (12).
p(θ ) = 6.25θ e −2.5θ
(12)
Case 3: Updating with observed data and moment information simultaneously Combining one piece of observed data t ′ = 1 and the moment information θ = G = 0.8 ,
6
the normalized posterior is updated to be
p (θ ) = 26.3671875θ 2e −3.75θ
(13)
with β = −1.75 . The prior and posteriors updated with different constraints are plotted in Fig. 2. The figure shows that moment information adds additional constraints on posteriors. Compared with Bayesian rules characterized by posterior updated with observed data, MRE updating procedure generates different posteriors in this example. The different posteriors give different predictions in practice. Prior and posteriors of parameter θ Prior Posterior (updating with t'=1)
1
Posterior (updating with 〈θ〉 =0.8) Posterior (updating with t'=1 and 〈θ〉 =0.8)
PDF(θ)
0.8
0.6
0.4
0.2
0
0
1
2
3
4
5
6
θ
Fig. 2 The prior and posteriors updated with different constraints
3. MRE updating for fatigue prognosis The fatigue prognosis is important for structure health management and residual life prediction. From the previous sections, a general form of the entropy posterior for is known to be in the form of Eq. (7). In order to use MRE updating for fatigue prognosis, a physical model is required to deterministically describe the fatigue crack growth. The prior statistical distributions of model parameters are required to incorporate the parameter estimation uncertainties and a likelihood function needs to be constructed using the physical model to characterize the model and measurement uncertainties. Detailed derivation and implementation of MRE updating to fatigue prognosis are shown below.
7
Deleted: Fig. 2
3.1. Deterministic crack growth model The deterministic model is used to describe the underlying physics of the fatigue damage accumulation. It is noted that the proposed MRE updating can be applied to different physical models depending on the specific problem. For demonstration and validation purpose, Paris law for the fatigue crack growth is chosen as the deterministic model (Paris and Erdogan, 1963). It is widely used to describe crack growth rate curve under constant amplitude cyclic loading:
da m = c(∆K ) dN
(14)
In Eq. (14), a is the crack length, ∆K is the stress intensity factor range, and (c, m ) are material parameters. For a plate specimen with a center through crack, ∆K is defined as:
a ∆K = πa ∆σF w
(15)
where w is the specimen width and F (⋅) is geometry correction factor for the specimen. ∆σ is the applied stress range. Given c , m and cycle number N , the crack size a can be calculated by solving the ODE equation numerically.
da m = c(∆K ) dN
πa ∆σ =c a F w
m
(16)
3.2. Building the likelihood equation Consider a general model prediction equation, d = y +τ
(17)
where d is the observed data, y is the model prediction, and τ is the general error term, including the measurement error ε and the modeling error e . The model prediction is given by Eq. (18) as y = M (x | θ ) + e
(18)
where M ( x | θ ) is the deterministic model and e is the model prediction error term, x is the prediction value and θ is the deterministic model parameter vector. Equation (17) can be rewritten as 8
d = M (x | θ ) + e + ε
(19)
In general, the probability distributions of the two uncorrelated error terms e and ε can be described using two Gaussian distributions with standard deviations σ e and σ ε , respectively. Let d be the measure of the response variable. The probability distributions of the two error terms
are
p(e | θ ) =
[d − M ( x | θ )]2 1 exp− 2σ e2 2π σ e
(20)
ε2 exp − 2 2π σ ε 2σ ε
(21)
and 1
p(ε | θ ) =
The likelihood function can be constructed as
L (D | θ ) = p (d | θ ) =
[d − M (x | θ )]2 exp− 2 σ e2 + σ ε2 2π σ e2 + σ ε2 1
(
)
(22)
Using general error term τ ~ Gaussian (0, σ τ ) , Eq. (22) can be rewritten as [d − M (x | θ )]2 1 exp− 2σ τ2 2π σ τ
L (D | θ ) = p (d | θ ) =
(23)
For multiple observed data, the likelihood is the product of the form L ( DN | θ ) =
n
∏
p(d i | θ ) =
i =1
1
(
2π σ τ
)
n
n [d − M (xi | θ )]2 exp − i 2σ τ2 i =1
∑
(24)
Substituting Eq. (24) in Eq. (7) and omitting the constant terms, the posterior is derived to be
p (θ ) ∝ µ (θ ) ⋅
n [d i − M (xi | θ )]2 β ⋅g (θ ) exp − ⋅e σ τn 2σ τ2 i =1 1
∑
(25)
3.3. MCMC sampling The analytical posterior probability distribution of c and m in fatigue prognosis is difficult to evaluate because the normalization of Eq. (25) involves multi-dimensional integration. In order to circumvent the difficulty of direct evaluation of the posterior probability distribution described by Eq. (25), a cascade MCMC sampling technique is employed in this paper. MCMC simulation was first introduced by Metropolis et al. (1953) as a method to 9
simulate a discrete-time homogeneous Markov chain. The random walk algorithm of MCMC ensures the samples converge to the posterior probability distribution from an arbitrary initial value. The widely used algorithm of Metropolis-Hastings (Hastings, 1970) is summarized below. Detailed derivation can be found in the referenced article. The transition between two successive samples xt and xt +1 is defined as x ~ q ( X | xt ) with probabilit y α (xt , ~ x) ~ xt +1 = x else t
(26)
where q( X | xt ) is the transition distribution. The acceptance probability α ( xt , ~ x ) is
α (xt , ~x ) = min (1, r )
(27)
The Metropolis ratio is given by r=
p (~ x ) q( xt | ~ x) p ( xt ) q(~ x | xt )
(28)
where p (⋅) is the posterior probability representation. For a symmetric transition distribution
q(⋅) , for example the normal distribution, the property of
q( xt | ~ x ) = q (~ x | xt ) simplifies the
Metropolis ratio to be
r=
p (~ x) p( xt )
(29)
The merit of MCMC is that it overcomes the normalization in Bayesian and the proposed MRE updating procedure. In this paper, a cascade MCMC (Tarantola, 2005) is employed to generate the posterior samples of Eq. (25).
4. Demonstration example and model validation A demonstration example of fatigue crack growth prognosis is illustrated in this section. Fatigue crack growth measurement data are used to demonstrate the application of the proposed maximum relative entropy updating procedure. The accuracy and precision of the prediction results are compared with experimental observations.
10
4.1. Crack growth data A large fatigue crack growth database of 2024-T3 aluminum alloy was reported by Virkler et al. (Virkler et al., 1979). The data set consists of 68 sample trajectories, each containing 164 measurement points. The entire specimen has the same geometry, i.e., an initial crack size ai = 9mm , length L = 558.8mm , width w = 152.4mm and thickness d = 2.54mm . The stress range during each experiment is constant ∆σ = 48.28MPa , and the stress ratio is R = 0.2 . The failure criterion is that the crack size equals 49.8mm . Kotulski reported the statistics of the Paris model parameters (c, m ) . The mean and coefficient of variance (COV) of m are 2.874 and 0.057, respectively. The mean and COV of log(c) are -26.155 and 0.037, respectively (Kotulski, 1998). The geometry correction factor for Virkler’s experiments is
a F = w
1
a cos π w
, for
a < 0.7 w
(30)
Based on Kotulski’s report, the prior PDF of c can be expressed as,
(log(c ) − ζ c )2 exp− 2σ c2 2π σ cc 1
µ (c ) =
(31)
and the prior of the parameter m is
µ (m ) =
(m − ζ m )2 exp− 2σ m2 2π σ m 1
(32)
with ζ c = −26.155 , σ c = COVc ⋅ ζ c = 0.037 ⋅ 26.155 = 0.968 and ζ m = 2.874 , σ m = COVm ⋅ ζ m = 0.057 ⋅ 2.874 = 0.164
4.2. MRE updating with MCMC Substituting the prior µ (θ ) in Eq. (25) with Eq. (31) and Eq. (32), the final form of maximum entropy posterior for the fatigue crack growth problem can be expressed as
11
1 log(c ) − ζ 2 1 1 m − ζ 2 1 c m + β c g c (c ) ⋅ + β m g m (m ) p(c, m ) ∝ exp− exp− σ cc σc 2 σ m 2 σ m 1 n d − M (a | c, m, N ) 2 i i ⋅ n exp− ∑ i στ 2 σ i =1 τ
(33)
1
where M (ai | c, m, N i ) is the Paris deterministic model, d i is the measure of crack size corresponding to cycles N i , and ai is the model prediction of crack size with parameters (c, m ) at N i cycles. The total error term σ τ is set to be 0.1mm . The exponential terms are g c (c ) = log(c ) and g m (m ) = m . The maximum relative entropy updating procedure using the cascade MCMC is described in the following steps. When one observation of d i , N i is available, the proposed MRE updating procedure begins with initial values of β c = β m = 0 . Step 1: Generate statistically sufficient number of (c, m ) samples by the cascade MCMC. Step 2: Calculate crack size a at any interested cycle count N using MCMC samples generated in step 1. This step is for crack size prediction purpose. Step 3: Calculate the coefficient β c , β m in the exponential term with the statistics of MCMC samples for the next updating. If additional observed data become available, the above mentioned steps are repeated till the final failure. One crack growth trajectory in Virkler’s dataset is selected for crack growth prediction
Deleted: arbitrarily
updating for demonstration purpose from previous publication (Ostergaard and Hillberry, 1983). Six points are randomly chosen to simulate the measured values of crack length a obtained by Deleted: Table 1
monitoring systems. The selected observations are shown in Table 1. Table 1 Data used for updating (Virkler’s dataset) Number
1
2
3
4
5
6
Crack size (a )
9.8273
10.5830
11.4121
12.1658
14.5744
20.1360
Cycle (N )
23988
44227
59220
73088
105697
153673
One hundred thousand samples with 5% burn-in period are generated by MCMC sampling in 12
step 1. 99.9% confidence intervals of predicted crack size a are calculated. The prognosis results using different measurements are shown in Fig. 3. It is observed that the confidence intervals of crack length prediction become narrower as more and more measurement data become available for MRE updating. It is also shown that the early stage measurement data can be used to improve model prediction significantly which is essential for condition-based maintenance and unit replacement. (a) Updating with data 1
(b) Updating with data 1-2
50
50 99.9% interval prediction Median prediction Experimental data Updating points
45
45 40
35
35
Crack length (mm)
Crack length (mm)
40
30 25 20
30 25 20
15
15
10
10
5
0
0.5
99.9% interval prediction Median prediction Experimental data Updating points
1 1.5 Number of cycles
5
2
0
0.5
5
x 10
(c) Updating with data 1-3
5
x 10
50 99.9% interval prediction Median prediction Experimental data Updating points
45 40
40
35 30 25 20
35 30 25 20
15
15
10
10
0
0.5
1 1.5 Number of cycles
99.9% interval prediction Median prediction Experimental data Updating points
45
Crack length (mm)
Crack length (mm)
2
(d) Updating with data 1-4
50
5
1 1.5 Number of cycles
5
2 5
x 10
13
0
0.5
1 1.5 Number of cycles
2 5
x 10
(e) Updating with data 1-5
(f) Updating with data 1-6
50
50 99.9% interval prediction Median prediction Experimental data Updating points
45
45 40
35
Crack length (mm)
Crack length (mm)
40
30 25 20
35 30 25 20
15
15
10
10
5
0
0.5
99.9% interval prediction Median prediction Experimental data Updating points
1 1.5 Number of cycles
5
2
0
0.5
5
x 10
1 1.5 Number of cycles
2 5
x 10
Fig. 3 Prognosis results using different measurements
5. Comparison with Bayesian approach Bayesian approach is widely used for reliability updating of engineering systems. A detailed comparison between the proposed MRE updating and the classical Bayesian updating is given in this section to show the advantages of the proposed MRE approach. Considering a Bayesian posterior pB (c, m ) for the same application example: pB (c, m ) ∝
1 log(c ) − ζ 2 1 1 m − ζ 2 1 c m ⋅ exp− exp− σ cc σ σ σ 2 2 c m m 1 n d − M (a | c, m, N ) 2 i i ⋅ n exp− ∑ i στ στ 2 i =1
(34)
1
Compared with MRE posterior in Eq. (33), the Bayesian posterior has no exponential terms introduced by moment constraints. In order to give a distinct comparison between Bayesian and MRE updating, all the parameters for prior and MCMC sampling in Bayesian updating procedure are the same with the ones used in the proposed MRE updating. Besides the dataset collected by Virkler’s experiments, another dataset acquired by McMaster on aluminum alloy T351 is also employed for extensive comparison between these two methodologies.
Formatted Formatted: heading2, Left, Level 1, Widow/Orphan control
5.1 Virkler’s dataset
All the physical and experimental description about the samples are referred to the previous section. The predictions of crack size using Bayesian approach and MRE approach are plotted in 14
Deleted: Fig. 4
Fig. 4. (a) Updating with data 1
(b) Updating with data 1-2
50
50 99.9% interval prediction (Bayasian) Median prediction (Bayesian) 99.9% interval prediction (MrE) Median prediction (MrE) Experimental data Updating points
Crack length (mm)
40 35
45 40
Crack length (mm)
45
30 25 20
35 30 25 20
15
15
10
10
5
0
0.5
1 1.5 Number of cycles
99.9% interval prediction (Bayasian) Median prediction (Bayesian) 99.9% interval prediction (MrE) Median prediction (MrE) Experimental data Updating points
5
2
0
0.5
5
x 10
(c) Updating with data 1-3
35
40
30 25 20
35 30 25 20
15
15
10
10
0
0.5
1 1.5 Number of cycles
99.9% interval prediction (Bayasian) Median prediction (Bayesian) 99.9% interval prediction (MrE) Median prediction (MrE) Experimental data Updating points
45
Crack length (mm)
Crack length (mm)
40
5
2
0
0.5
5
x 10
(e) Updating with data 1-5
1 1.5 Number of cycles
2 5
x 10
(f) Updating with data 1-6
50
50 99.9% interval prediction (Bayasian) Median prediction (Bayesian) 99.9% interval prediction (MrE) Median prediction (MrE) Experimental data Updating points
40 35
40
30 25 20
35 30 25 20
15
15
10
10
0
0.5
1 1.5 Number of cycles
99.9% interval prediction (Bayasian) Median prediction (Bayesian) 99.9% interval prediction (MrE) Median prediction (MrE) Experimental data Updating points
45
Crack length (mm)
45
Crack length (mm)
5
x 10
50 99.9% interval prediction (Bayasian) Median prediction (Bayesian) 99.9% interval prediction (MrE) Median prediction (MrE) Experimental data Updating points
45
5
2
(d) Updating with data 1-4
50
5
1 1.5 Number of cycles
5
2 5
x 10
0
0.5
1 1.5 Number of cycles
2 5
x 10
Fig. 4 Comparison of prognosis results between Bayesian updating and MRE updating (Virkler’s dataset) The comparison between MRE updating and the classical Bayesian updating in the application example shows that at the very beginning stage of crack growth, in particular the 15
prediction generated by first two measurements of response shown in Fig. 4 (a) and (b), both
Deleted: Fig. 4
MRE and Bayesian updating have similar predictions. However, as the crack grows to the early stage, characterized by response points 3 - 5 shown in Fig. 4 (a - c), MRE gives more satisfactory
Deleted: Fig. 4
prediction of confidence intervals of crack size, which is significant for residual life prediction. The narrower confidence intervals are due to the additional moment constraints introduced by MCMC samples. As shown in Fig. 4 (f), MRE still gives a narrower confidence interval
Deleted: Fig. 4
prediction compared with Bayesian but the difference is not as prominent as is in the previous updatings in Fig. 4 (c - e). One possible explanation is that the last measurement used in the
Deleted: Fig. 4
updating procedure is near the middle of the crack growth curve, and all the data points are sufficient enough to dominate the posteriors. It is clear that if more and more data become available, this crack curve can be statistically determined using MRE updating procedure.
Formatted Formatted: heading2, Left, Level 1, Indent: First line: 0", Widow/Orphan control
5.2 McMaster’s dataset
McMaster and Simth (1999) collected a set of 2024-T351 aluminum alloy experimental data under constant and variable loading conditions. The experimental data of center-cracked specimens with length L = 250mm , width w = 100mm and thickness t = 6mm under constant loading with ∆σ = 65.7 MPa and stress ratio R = 0.1 are used here for comparison using Paris’ law. The details of experimental setup and conditions can be found in referenced articles. Six points are picked up to simulate the measures of response in Table 2. One hundred thousand samples
Formatted: Lowered by 3 pt Formatted: Lowered by 3 pt Formatted: Lowered by 3 pt Formatted: Lowered by 3 pt Deleted: Table 2
are generated for both Bayesian and MRE updatings. Prediction curves at sixth updating together with a prior prediction are given in Fig. 5 in order to concentrate the results on the comparison. Table 2 Data used for updating (McMaster’s dataset) Number
1
2
3
4
5
6
Crack size (a )
11.3611
11.9282
12.3254
13.8563
14.8771
19.5841
Cycle (N )
4875
8474
11550
17775
21375
34500
Deleted: Fig. 5 Formatted: Check spelling and grammar
Formatted: Justified
Despite of an inaccurate prior curve presented in Fig. 5, MRE updating exhibits narrower intervals than that of Bayesian because of the moment constraints added by MCMC samples in MRE framework. In statistical point of view, as more data are available for such a probabilistic inference. Posterior PDFs given by these two methods will finally converge to the actual target 16
Deleted: Fig. 5 Formatted: Check spelling and grammar Formatted: Justified, Indent: First line: 0.5"
distribution. However, in practice, particularly for large and complex system, the measurements of response variables are expensive and hard to acquire and statistically sufficient amount of data are virtually impossible. At the same time, time-based preventive maintenance might cost too much for a complex system. In this case, it is beneficial to use MRE updating procedure to generate narrower confidence intervals such that early stage warnings can be generated based on those predictions to facilitate condition-based maintenance of the system. Formatted: Keep with next 99.9% interval prediction (Bayesian) Median prediction (Bayesian) 99.9% interval prediction (MRE) Median prediction (MRE) Experimental data Prior estimation Updating points
35
Crack length (mm)
30
25
20
15
10 0
0.5
1
1.5
2 2.5 3 Number of cycles
3.5
4
4.5
5 4
x 10
Fig. 5 Comparion of preciton between Bayesian and MRE updating (McMaster’s dataset)
Formatted: Centered, Level 1, Indent: First line: 0" Deleted: 5
6. Conclusion
Field Code Changed
A general maximum relative entropy approach is developed in this paper to include all available information and uncertainties for updating and prediction. For the unstable and stochastic fatigue crack growth process, the uncertainties include model prediction uncertainty, measurement uncertainty, and model parameter uncertainty. The maximum relative entropy updating procedure based on the measures of response quantities is presented in detail. The posterior samples are generated using a cascade MCMC random walk algorithm. Beyond the classical Bayesian updating procedure, the statistics of MCMC samples is employed to add additional constraints on the posterior. Real application examples reveal that the posterior model
Deleted: A r Deleted: s
obtained by maximum relative entropy approach significantly improves the prediction. Furthermore, comparisons of crack growth prediction between maximum relative entropy and Bayesian updating have been performed to show the benefits of maximum relative entropy 17
Deleted: a c Deleted: s
approach. The comparison shows that MRE updating in general gives more satisfactory prediction. It further indicates the potential that the maximum relative entropy approach can be applied to engineering systems as an improved alternative to the classical Bayesian approach.
Acknowledgement The research reported in this paper was supported by funds from National Aeronautics and Space Administration (NASA) (Contract No. NNX09AY54A). The support is gratefully acknowledged.
References Kullback, S., Leibler, R.A. (1951). On Information and Sufficiency. The Annals of Mathematical Statistics, 22(1):7986. Metropolis, N., Rosenbluth, A.W., Rosenbluth, M.N., Teller, A.H., Teller, E. (1953). Equations of state calculations by fast computing machine. Journal of Chemical Physics, 21(6):1087-1092. Jaynes, E.T. (1957). Information theory and statistical mechanics. Physical Review, 106: 620- 630 and 108: 171190. Paris, P.C., Erdogan, F. (1963). A critical analysis of crack propagation laws. Journal of Basic Engineering, Transactions of American Society of Mechanical Engineers, 85:528-534.
Foreman, R.G., Kearney, V.E., Engle, R.M. (1967). Numerical analysis of crack propagation in a cyclic-loaded structure. Journal of Basic Engineering, Transactions of American Society of Mechanical Engineers, 89D: 459-464. Walker, K. (1970). The effect of stress ratio during crack propagation and fatigue for 2024-T3 and 7075-T6 aluminum. In Effects of environment and complex load history on fatigue life. ASTM STP 462, (pp. 1-14). American Society for Testing and Materials. Hastings, W.K. (1970). Monte Carlo Sampling Methods Using Markov Chains and Their Applications, Biometrika, 57(1):97-109. Jaynes, E.T. (1979). Where do we stand on Maximum Entropy. In R.D. Levine and M. Tribus (Eds.), The Maximum Entropy Formalism (pp. 15-119). M.I.T. Press, Cambridge, MA.
Virkler, D.A., Hillberry, B.M., Goel, P.K. (1979). The statistical nature of fatigue crack propagation, ASME Journal of Engineering Materials and Technology, 101:148-153.
Ostergaard, D.F., Hillberry, B.M. (1983). Characterization of the variability in fatigue crack propagation data, In J.M. Bloom and J.C. Ekvall (Eds.), Probabilistic Fracture Mechanics and Fatigue Methods: Application for Structural Design and Maintenance, ASTM STP 798 (pp. 97-115). American Society for Testing and Materials.
Skilling, J. (1988). The axioms of maximum entropy. In G.J. Erickson and C.R. Smith (Eds.), Maximum Entropy and Bayesian Methods, in Science and Engineering, volume 1, (pp. 173-187). Kluwer Academic Publishers, Netherlands.
18
Madsen, H. (1997). Stochastic modeling of fatigue crack growth and inspection. In C.G. Soares (Ed.), Probabilistic methods for structure design, collection Solid Mechanics and its applications, (pp. 59-83). Kluwer Academic
Publishers, Netherlands. Kotulski, Z.A. (1998). On Efficiency of Identification of a Stochastic Crack Propagation Model based on Virkler Experimental Data. Archives of Mechanics, 50(5):829-847. McMaster, F.J., Smith, D.J. (1999). Effect of load excursions and specimen thickness on crack closure measurements, In R.C. McClung and J.C. Newman, Jr., (Eds) Advances in Fatigue Crack Closure Measurements and Analysis: Second Volume, ASTM STP 1343, (pp. 246-264). West Conshohocken, Pennsylvania: American
Society for Testing and Materials. Zhang, R., Mahadevan, S. (2000). Model uncertainty and Bayesian updating in reliability-based inspection. Structure Safety, 22(2):145-160.
Caticha, A., Preuss, R. (2004). Maximum entropy and Bayesian data analysis: Entropic prior distributions. Physical Review E. 70(4), 046127.
Tarantola, A. (2005). Inverse Problem Theory and Methods for Model Parameter Estimation. Society for Industrial and Applied Mathematics, 3600 University City Science Center, Philadelphia, USA. Caticha, A. & Giffin, A. (2006). Updating probabilities. In A.M. Djafari (Ed.), Bayesian Inference and Maximum Entropy Methods in Science and Engineering, AIP Conference Proceedings, 31: 872.
Perrin, F., Sudret, B., Pendola, M. (2007). Bayesian updating of mechanical models – Application in fracture mechanics. In 18ème Congrès Français de Mécanique, CFM 2007, Grenoble. Giffin, A. & Caticha, A. (2007) Updating probabilities with data and moments. In K. Knuth (Ed.), Bayesian Inference and Maximum Entropy Methods in Science and Engineering, AIP Conference Proceedings, 954:74
Tseng, C., Caticha, A. (2008). Using relative entropy to find optimal approximations: An application to simple fluids. Physica A, 387:6759-6770.
19
