CC BY-NC-ND
0Available online at www.sciencedirect.com
Procedía
ELSEVIER
Procedía IUTAM 6 (2013) 123 - 131
www.elsevier.com/locate/procedia
IUTAM Symposium on Multiscale Problems in Stochastic Mechanics 2012
Damage tolerance and reliability assessment under random Markovian loads
Jean-Marc Bourinetab*, Cecile Mattrandab
aClermont Université, IFMA, Institut Pascal, BP 10448, 63171 Aubiere, France bCNRS, UMR 6602, Institut Pascal, 63171 Aubiere, France
Abstract
Accounting for load uncertainties plays an important role in the design of safe structural components of aircrafts under damage tolerance requirements. The purpose of this paper is to develop a reliability assessment technique for cracked structures submitted to non-stationary random fatigue loads modeled by first-order Markov chains with discrete state space and identified from in-flight measurements. The strategy based on a multi-level version of the cross-entropy method consists in progressively updating the transition probability matrix in order to generate load sequences of increasing severity which are likely to cause failure. The proposed method is applied to a cracked M(T) specimen under the defined random fatigue loads. Load cycle interactions and retardation effects are accounted for by means of the PREFFAS crack closure model. The efficiency of the proposed approach in terms of computational cost is clearly observed for rare failure events in comparison with direct Monte Carlo simulations. In addition to the failure probability estimate, the multi-level cross-entropy method provides the analyst with information on the most probable load sequences at failure.
© 2013 The Authors. Published by Elsevier B.V.
Selection and(or peer review under responsibility of Karlsruhe Institute of Technology (KIT) Institute of the Engineering Mechanics. Keywords: damage tolerance ; reliability assessment; importance sampling ; cross-entropy method ; first-order Markov chains
1. Introduction
The damage tolerance design of aircraft structural components is mainly based on the assumption of a deterministic approach of the crack growth process. Real environmental conditions however exhibit various and significant sources of uncertainty which need to be accounted for. The main sources are known to come from material properties, length/location/orientation of initial cracks and applied loads. The work presented here focuses on the latter source of uncertainty of real importance for structures subjected to significantly scattered loads and when interaction effects between cycles along the random loading sequence need to be assessed in order to avoid overconservative designs. Non-stationary variable amplitude fatigue load sequences are here modeled by means of first-order Markov chains
* Corresponding author. Tel.: +33-4-73-28-81-16 ; fax: +33-4-73-28-80-27 . E-mail address: bourinet@ifma.fr
2230-9838 © 2033 The Authors. Published by Elsevier B.V.
Selection and/or peer review under responsibility of Karlsruhe Institute of Technology (KIT) Institute of the Engineering Mechanics. doi:30.3036/j.piutam.2033.03.034
with discrete state space. Their parameters are identified from in-flight measurements recorded on a fleet of fighter aircrafts. The constructed random fatigue load model is then used for the reliability assessment of a simplified damage tolerance problem. For this purpose, an importance sampling strategy based on the Cross-Entropy (CE) method [1] is developed for an efficient and accurate estimation of low failure probabilities. The algorithm consists in updating the transition matrix of the Markov chain in a multi-level approach so that the fatigue load sequences generated according to the new transition matrix progressively lead to an increasing number of failures. The proposed method is applied to the reliability assessment of a cracked M(T) specimen under random fatigue loads. Results obtained by the proposed approach are compared both in terms of accuracy and efficiency with those obtained with a crude Monte Carlo simulation (MCS) considered as reference results.
The paper is organized as follows. Section 2 briefly recalls the probabilistic model identified from the recorded data, based on previous works of the authors [2,3]. The reliability problem is formulated in Section 3 and the solving strategy based on the CE method is detailed in Section 4. The application example is presented in Section 5, where the strengths of the CE method are highlighted. A conclusion is drawn in Section 6.
2. Probabilistic model of the fatigue loading
A recourse to continuous time processes for modeling random load sequences has been a common practice in probabilistic fracture mechanics. Several works were carried out in the fields of fatigue initiation [4, 5, 6] and fatigue crack propagation [7, 8, 9, 10, 11, 12, 13, 14, 15]. A stationary Gaussian random process is most often conveniently assumed and its parameters (e.g. the shape of the power spectral density function) are selected or varied regardless of any real load data. Conversely, when load data are available, Markov processes have been considered as interesting alternatives [16, 17, 6]. This type of models has been investigated by the present authors in previous works [2, 3] for the following reasons:
• the identification of a suitable continuous time process was not straighforward due to the form of the recorded data (variable-length time steps selected according to the aircraft activity, variable flight durations),
• the load time series recorded on a fleet of fighter aircraft were obviously non-stationary,
• Markov processes are convenient for a direct modeling of sequences of max-min pairs expected for the purpose of crack propagation (with continuous-time processes, an additional filtering stage is required for obtaining such sequence from time trajectories),
• the recorded data were provided in a sufficient amount for the identification of a Markov process.
It is worth mentioning that the main objective here is the identification of a suitable process or a few processes from the data, without any intention to model the physical and complex loading phenomena involved during the flights. In this paper and for the sake of simplicity, we will focus on the randomness of fatigue loads observed in a single flight for a selected type of mission of the aircraft.
In the sequel we consider a N-cycle load sequence X ,...,Xn,... ,Xn ) where Xn is taken here as a cycle Xn = (Mn, mn) and Mn (resp. mn) stands for the peak (resp. trough) stress value of the nth cycle. We assume a first-order dependence in terms of load cycles, which implies a greater dependence in terms of stress levels. The First-order Markov Chain (FMC) with discrete state space E therefore satisfies the following property, for all n £ N*:
where Si, Sj represent two given stress levels taken from a set of Kc levels. These Kc levels are selected according to the shape of the empirical distribution of the measured stress levels, see Fig. 1 (a). The state space E encompasses K = Kc (Kc — 1) /2 cycle states such that a valley mn is systematically followed by a peak Mn. The present work is based on group "B" data of reference [3]: Si, Sj £ {0.039, 0.113, 0.248, 0.507, 0.840} and we therefore have Kc = 5 and K = 10.
P ( Xn+1 = Xn+1 I Xn = x^ ... X1 = X1 ) = P ( Xn+1 = Xn+1 I Xn = Xn )
The state space E is composed of a finite number K of load cycle states:
E = { ek = (si, Sj), Si > Sj and i, j e{1,... } }
0.09 0.08 0.07 0.06 , 0.05 0.04 0.03 0.02 0.01 0
0.014 0.012 0.01 0.008 " 0.006 0.004 0.002 0
25 50 75 100 125 150 Number of cycles per flight
Fig. 1. (a) Distribution of normalized peak/trough stresses (all stresses divided by the maximum value of the recorded data);(b) Transition matrix P (load cycle states ordered as follows: {ei = (s2, si); e2 = (S3, si); «3 = (S3, S2); •••; e9 = (S5, S3); eio = (S5, S4)}); (c) Distribution of Xi; (d) Distribution of N.
The probabilities of moving from a given load cycle to the next one are gathered in a K x K transition probability matrix P supposed constant over time n (time-homogeneous Markov chain), represented in Fig. 1 and which satisfies the following conditions:
0 <pi,j < 1 for i, j e{1,...,K]
XPi,j = 1 for ie{i,...,K>
where pt,j = P (Xn+1 = ej | Xn = e\) = P (X2 = ej | X1 = e\), for all n e N* and any (e\, ej) e E x E.
The FMC is therefore completely defined by its transition matrix P, its initial distribution X1 and its length N (see Fig. 1).
3. Formulation of the reliability problem
0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
G / G..
We assume that failure is caused by a crack growth instability, when the stress intensity factor K in opening mode (mode I) exceeds the fracture toughness Kc at a reference stress level Oi as defined in Eq. 4. This stress level Oi could be seen here as representative of the limit load of the aircraft as defined by aviation authorities.
K (a, ol) > Kc (4)
The crack size a in Eq. 4 represents the damage accumulation from an initial crack size a0 under a given trajectory of the previuously defined FMC, which consists of a variable amplitude load sequence. The crack extension is evaluated by means of a crack growth prediction model. Two models are used in the present work: the Paris-Erdogan law and the PREFFAS crack closure model [18]. This latter model presents the major advantage to account for potential retardation or acceleration effects on the crack growth process due to overloads/underloads in the variable amplitude load sequences. Such effects are important in the damage tolerance design of structural aircraft components and they need to be assessed in order to avoid overconservative designs.
For a reliability problem of a reduced complexity, we will assume that all load sequences (or flights) are composed of a given deterministic number of cycles N and that all FMC start from the same initial cycle x3 = e2 = (s3, s3) = (0.248, 0.039). The limit-state function g is expressed as follows, for x € EN:
g(x) = g(xi<n<N) = ac - a (xi<n<N) (5)
where ac is the critical crack length solution of K(ac, al) = Kc and a (x3<n<N) is the crack length resulting from the propagation under the N-cycle load sequence x3<n<N.
The failure probability is given by the following expression:
Pf = \ px (dx) = f px (x) dx (6)
IN jDfN
where PX denotes the probability measure over the observable space EN, pX refers to the probability density function (pdf) of the Markov chain X and Dfn = { x € En : g(x) < 0 } represents the failure domain.
According to the first-order Markov property, the probability of the homogeneous FMC X reads, by successive conditioning:
P (X) = P (X3 = xi) n P (Xn+3 = xn+31 Xn = xn) (7)
The pdf of the FMC pX can therefore be written as:
px (x)= px (x,P)= n P>,Jnij(x) (8)
i, J=3
where niy j (X) denotes the number of transitions of the Markov chain X = (X3,... ,XN) from cycle ei to cycle ej. It is assumed that P(X3 = x3) = 3, i.e. all trajectories of the first-order Markov chain X start from a unique cycle state x3.
4. Solution based on the CE method
The probability of failure Pf solution of Eq. 6 is now expressed as the following expectation:
Pf = jD px(x) dx = JEN tDfN(x) px(x) dx = Epx [lDfN(X)] (9)
Pf = jENlDfN(x)^qx(x)6x = ^qx [lD/w(X)W(X)] (10)
where lDf is the indicator function of the failure domain DfN such that lDf (x) = 3 if x € DfN and lDf (x) = 0 otherwise.
A direct estimation of the failure probability by a crude Monte Carlo Simulation (MCS) from Eq. 9 is known to be too costly in the case of rare failure events. Importance Sampling (IS) constitutes an efficient alternative method for such a purpose. With IS, the reliability problem is rewritten in the following form:
Px(x) qx(x)
where qX is an instrumental pdf which must dominate 1 Df pX and W (x) = pX (x) /qX (x) is called the likelihood ratio. The corresponding IS statistical estimator of Pf from a set of Ns samples is given by:
where X(3),...,X(Ns) are Ns i.i.d. copies of X with pdf qX.
The pdf which leads to a zero-variance of Pf is denoted by qX. This optimal solution has the following expression: Pf
The CE method introduced by Rubinstein [1] and used in the present work consists in selecting the importance sampling density qX which is the closest to the optimal density q*X within the family of the nominal density {pX (•, Q)} indexed by the parameter transition matrix Q. The CE problem therefore consists in choosing the optimal parameter matrix Q such that the Kullback-Leibler distance between the densities q*X and pX (•, Q) is minimal:
Q* = argminD (q*x, px (•,Q)) Q
= argmin Eqx
logf-^-
px (x, Q)
After some straightforward calculations which make use of the expression of qX defined at Eq. (12), the optimal solution Q* takes the following form:
Q* = argmax Epx (^ [tDfN (X) log px (X, Q) 0 < qij < 1 for i, j e{1,...,K}
^ !>, j = 1 for i e {1,-.-,K} (14)
where the constraints are added so that Q is a transition matrix.
In the case of rare failure events, it is expected that most of the realizations of the indicator function 1Df are zeros when X is sampled from the original probability distribution pX (•, P), which leads to poor estimations of Q*. A multi-level CE procedure has been developped to circumvent this issue [19, 20]. The key idea consists in building a sequence of reference parameters {Qt, t > 0} and a sequence of thresholds {Y, t > 1} which are adaptively updated. At each level t, a CE solution is derived by IS. Assuming some arbitrary proposal density pX (•,R) for IS, Eq. (14) rewrites as follows:
Q* = argmax Epx (^ \tDfN (x) W (x, P, R) log px (x, Q) 0 < qij < 1 for i, j e {1,2,...,K}
Iqu = 1 for i e{1,2,...,K} (15)
',7=1 ' * '"' '
The statistical estimator QQ* is easily derived from the solution of the maximization problem given in Eq. 15:
where the likelihood ratio takes the form W (x, P, R) = px {x, P) /px {x, R) = nf ' —
1 Ns ivI^/J ivs k=1 x(k) (A < v / ni, j(x(k))
1 Ns — y lDf N DfN Ns k=1 (x (k)) (A i, j=1 ni(x(k))
for i, j e{1,...,K} (16)
where x(1),...,x(Ns) are sampled from px (•, R) distribution and ni(x) = Xj=i ni,j(x) represents the number of transitions of the Markov chain x starting from state ei. The multi-level CE algorithm is described below:
1. Set Qo = P and t = 1.
2. Generate Ns samples x(1),... ,x(k),... ,x(Ns) from px (•,and evaluate the corresponding limit-state values g(x(k)). Define % as the sample p-quantile of g(x(k)) where p is a not very small chosen parameter (p = 0.1 in the present work). 0
3. Use the same samples x(3),... ,x(Ns) to estimate q*j from Eq. (36) with ritj = qt-3. for i, j € {3,...,K}. The indicator function 1DfN (x(k)^ is replaced by 1D^ (t) (x(k)^ in Eq. (36), where DfN(t) = { x € EN : g(x) < Y }. Denote qti J the solution of Eq. (36).
4. If Y < 0, stop the algorithm and proceed with step 5, set t = t + 3 otherwise and reiterate from step 2.
5. Estimate the failure probability by IS:
where t denotes here the final number of iterations (i.e. number of levels used).
For an insufficient number of samples Ns, some transitions ei ^ ej may appear unobserved at some iteration of the algorithm despite non-zero (but small) transition probabilities. The corresponding estimated transition probabilities qtij then remain equal to zero until the end of the multi-level CE algorithm, which results in a transition matrix Qt obtained at the last level potentially far from the true optimal solution. Weighting the solution obtained at each level with either the nominal probabilities P or the solution of the previous level Qt-3 has been proposed to address such an issue [20]. This latter solution is used in the subsequent application treated in Section 5, with a smoothing parameter a set to 0.6:
, = (1 - a)c¡t-i¡j + a
5. Application example
f 1 w (x(k)) w(x(k), P ,Qt) m,j (X(k)) k=1 _
f idfnt) (x(k)) w(x(k),P,Qt) n,(X(k))
The reliability analysis is applied to the crack propagation of a M(T) specimen 350 mm width, 2 mm thickness, made of a 2024-T353 aluminum alloy and submitted to the random fatigue loading modeled by FMC as described in Section 2. Two crack growth models are investigated: the Paris-Erdogan law (case 3 and 2) and the PREFFAS crack closure model (case 3), see Table 3. The following parameters are used: C = 2.437 • 30-33 and m = 3.42 in the Paris law (SIF range AK consistent with MPaVmm and crack growth rate da/dN with mm/cycle), A = 0.45 and B = 0.55 for Elber constants in the PREFFAS model. An initial crack size a0 = 5 mm is assumed for all three cases. The critical crack length ac is arbitrary fixed for sufficiently low failure probabilities. For case 3, the whole load sequence is composed of a randomly generated 525-cycle load subsequence repeated 300 times in order to meet the stationary spectra assumption of the PREFFAS model [38]. Works are underway to improve the PREFFAS model and relax such an assumption in order to assess the propagation of crack under fully random load sequences.
Table 1. Application examples.
Crack growth model
Number of applied cycles N
Initial crack length ac (mm)
Critical crack length ac (mm)
Case 1 Case 2 Case 3
Paris law Paris law PREFFAS
500 102000 52500
Results obtained with the multi-level CE algorithm and crude MCS are listed in Table 2. For the CE method, the coefficient of variation (c.o.v.) is obtained empirically from 30 independent runs of the algorithm. The CE method clearly outperforms crude MCS. For an equivalent accuracy given in terms of c.o.v., the CE method requires a number of calls to the limit-state function which is several orders or magnitude less than the one of a crude MCS. The gain increases with low failure probabilities. The averaged value of the 30 estimates of the failure probability obtained with the CE method is close to the single MCS estimate, which allows us to conclude that the CE method has no significant bias.
Table 2. Reliability results.
Method Total number of samples Failure probability Pf c.o.v. of Pf (in %)
Case 1 MCS 107 8.39 10-6 10.9
CE(Ns = 1000;4 levels) 4 x 103 7.69 x 10-6 W 26.2
CE(JVS = 5000; 4 levels) 2 ■ 104 7.71 ■ 10~61» 5.82
Case 2 MCS 105 1.10 10"3 9.53
CE(Ns = 1000;3 levels) 3 x 103 1.11 x 10-3 W 10
CE(JVS = 5000; 3 levels) 1.5 104 1.11 ■ 10-31» All
Case 3 MCS 105 5.31 10-4 13.7
CE(Ns = 1000;3 levels) 3 x 103 5.06 x 10-4« 17.1
CE(JVS = 5000; 3 levels) 1.5 104 4.84 ■ 1(T4<*> 4.31
(,): averaged value over 30 independent runs of the multi-level CE algorithm
Fig. 2. Evolution of (Qt - P) with iteration t (PREFFAS model). (a) t = 1, Y = 2.08;(b) t = 2, Y = 1.12;(c) t = 3, Y = -0.20.
The biased transition matrix Qt at the last iteration of the multi-level CE algorithm brings also some additional information. The load sequences simulated with this matrix induce severe crack growths which are likely to cause failure. The evolution of Qt with iterations is represented in Fig. 2 for case 3 which involves the PREFFAS model (plotting of the difference between Qt and P). With the PREFFAS model, the crack growth is especially sensitive to the maximum stress of the loading sequence, which is repeatedly applied as pointed out in [3]. In this model, the maximum stress is responsible for the overall amount of retardation but it also contributes more than the other maxima to the crack growth when it is applied. Transition probabilities which increase at failure are therefore those involving the maximum stress, here s5, as naturally expected.
6. Conclusion
This paper presents a method for assessing failure probabilities of cracked structures submitted to random loads modeled by first-order Markov chains with discrete state space. This problem is of real interest for a safe design of structural components of aircrafts under damage tolerance requirements, for which it is of importance to account for fatigue load uncertainties and interactions effects between cycles. The reliability problem is solved by means of the CE method and its multi-level algorithm. It is worth pointing that the problem tackled here differs from those addressed in existing works, most of them in the field of communication networks. Failure is obtained here from a complex mechanical model whose input is a Markov chain whereas failure is directly given by the state of a Markov chain in other works of the literature. The results obtained clearly demonstrate the accuracy and efficiency of the CE method. Beside the failure probability estimate, this method also provides us with the most probable load sequences at failure which brings further details about failure. The application of this method to Markov chains with continuous state space will be the subject of a forthcoming paper.
Acknowledgements
The work presented in this paper was financially supported by DGA which is gratefully acknowledged. The authors
wish also to thank DGA/Techniques Aeronautiques in Toulouse for providing them with the load data.
References
[1] Rubinstein RY. Optimization of computer simulation models with rare events. European Journal of Operational Research. 1997;99(1):89-112.
[2] Mattrand C, Bourinet JM. Random load sequences and stochastic crack growth based on measured load data. Engineering Fracture Mechanics. 2011;78(17):3030-3048.
[3] Mattrand C, Bourinet JM, Lemaire M, Bernard P, Fogli M. Modeling and simulation of stochastic fatigue load sequences derived from in-flight load data. In: Proc. 13th AIAA Non-Deterministic Approaches Conference, Denver, Colorado, USA; 2011..
[4] Pitoiset X. Methodes spectrales pour une analyse en fatigue des structures metalliques sous chargements aleatoires multiaxiaux [Ph.D. thesis]. Universite Libre de Bruxelles; 2001.
[5] Benasciutti D, Tovo R. Spectral methods for lifetime prediction under wide-band stationary random processes. International Journal of Fatigue. 2005;27(8):867-877.
[6] Genet G. A statistical approach to multi-input equivalent fatigue loads for the durability of automotive structures [Ph.D. thesis]. Göteborg University; 2006.
[7] Domínguez J, Zapatero J. Effect of the loading spectrum and history length on fatigue life distribution under random loading. Engineering Fracture Mechanics. 1992 aug;42:925-933.
[8] Domínguez J, Zapatero J, Pascual J. Effect of load histories on scatter of fatigue crack growth in aluminum alloy 2024-T351. Engineering Fracture Mechanics. 1997;56(1):65-76.
[9] Zapatero J, Moreno B, Domínguez J. On the use of the strip-yield model to predict fatigue crack growth under irregular loading. Fatigue & Fracture of Engineering Materials & Structures. 1997;20(5):759-770.
[10] Zheng R, Ellingwood BR. Stochastic fatigue crack growth in steel structures subject to random loading. Structural Safety. 1998;20(4):303-323.
[11] Ustilovsky S, Arone R. Random fatigue crack growth in aluminum alloys. International Journal of Fatigue. 1999;21, Supplement 1:S275-S282.
[12] Moreno B, Zapatero J, Domínguez J. An experimental analysis of fatigue crack growth under random loading. International Journal of Fatigue. 2003;25(7):597-608.
[13] Beck AT, Melchers RE. Overload failure of structural components under random crack propagation and loading - a random process approach. Structural Safety. 2004;26(4):471-488.
[14] Zapatero J, Moreno B, Gonzalez-Herrera A, Domínguez J. Numerical and experimental analysis of fatigue crack growth under random loading. International Journal of Fatigue. 2005;27(8):878-890.
[15] Wu WF, Ni CC. Statistical aspects of some fatigue crack growth data. Engineering Fracture Mechanics. 2007;74(18):2952-2963.
[16] Rychlik I. Simulation of load sequences from rainflow matrices: Markov method. International Journal of Fatigue. 1996;18:429-438.
[37] Johannesson P. Rainflow analysis of switching Markov loads [Ph.D. thesis]. Mathematical Statistics, Centre for Mathematical Sciences, Lund Institute of Technology; 3999.
[38] Aliaga D, Davy A, Schaff H. A simple crack closure model for predicting fatigue crack growth under flight simulation loading. ASTM STP. 3988;982:493-504.
[39] de Boer PT, Kroese D, Mannor S, Rubinstein RY. A tutorial on the cross-entropy method. Annals of Operations Research. 2005;334:39-67.
[20] Ridder AD. Importance sampling simulations of Markovian reliability systems using cross-entropy. Annals of Operations Research. 2005;334:3 39-336.
