Monte-Carlo simulations of life expectancy using the dual boundary element method

Engineering Fracture Mechanics 68 (2001) 1371±1384

www.elsevier.com/locate/engfracmech

Monte-Carlo simulations of life expectancy using the dual boundary element method Hocine Kebir a,*, Jean Marc Roelandt a, Jocelyn Gaudin b a b

LG2mS UPRESA CNRS 6066, Universite de Technologie de Compiegne, BP 20529, F-60205 Compiegne, France AEROSPATIALE MATRA, Centre Commun de Recherches 12, rue Pasteur BP 76, 92152 Suresnes Cedex, France Received 27 March 2000; received in revised form 26 February 2001; accepted 13 March 2001

Abstract Ageing aircraft structures are susceptible to a very typical and dangerous form of fatigue damage named ``multiple site damage''. An overall assessment method has been developed in order to improve the understanding of this phenomenon and to provide inspection programs. A full deterministic approach has been completed ®rst. An automatic numerical process, using the dual boundary element method and fracture mechanics laws, computes the crack evolution until failure of the structure. Then, in order to take materials scatter into account, a probabilistic approach has been developed based on Monte-Carlo simulations. Initial cracks scenarios are randomly de®ned and cracks evolution is computed for each of them. Thus statistical results are provided. The results obtained are in good agreement with experimental tests carried out at the Aerospatiale Matra research center. Ó 2001 Elsevier Science Ltd. All rights reserved. Keywords: Dual boundary element method; Fracture mechanics; Crack propagation; Monte-Carlo simulation

1. Introduction Multiple site damage (MSD) is a type of cracking that may be found in ageing airplanes which can adversely a€ect the damage tolerance of an airframe structure. It is de®ned as the simultaneous occurrence of fatigue cracks in the same structural element. This phenomenon generally occurs in areas where a lot of critical sites (fastener holes for example) are subjected to the same stress level, as in the riveted joints of a fuselage. Multi-site crack growth can greatly weaken any part of the structure. The occurrence of MSD in older aircraft was highlighted by the in-¯ight failure of a portion of the fuselage of an Aloha Airlines Boeing 737 in April 1988 (Fig. 1). At the time of the incident, the plane had accumulated nearly 90 000 ¯ights. The failure was precipitated by the linkup of small fatigue cracks emanating from adjacent rivet holes in the lap joint of the fuselage [1]. Currently, approximately 50% of the jet airplanes in the ¯eet are over 15 years old. The increasing number of in-service ageing aircraft as well as new forthcoming airworthiness requirements have incited *

Corresponding author. Tel.: +33-44-23-44-23; fax: +33-44-23-46-89. E-mail address: [email protected] (H. Kebir).

0013-7944/01/$ - see front matter Ó 2001 Elsevier Science Ltd. All rights reserved. PII: S 0 0 1 3 - 7 9 4 4 ( 0 1 ) 0 0 0 4 2 - X

1372

H. Kebir et al. / Engineering Fracture Mechanics 68 (2001) 1371±1384

Fig. 1. In-¯ight failure of a portion of the fuselage of an Aloha Airlines Boeing 737 [1].

manufacturers to develop methods to better understand the behavior of multiple cracked structures. The main objective is to provide engineering models to assess MSD regarding the problems of maintenance programs of ageing aircraft as well as the design of new ones. As in the case of single crack, the MSD problem may be divided into di€erent phenomena: crack initiation, crack propagation and residual strength. The main particularity in the case of MSD is that these three stages are strongly linked. Behavior of one crack will depend on the behavior of adjacent cracks. Residual strength of the structure will be deeply in¯uenced by the initial cracks con®guration and the development of cracks patterns [7]. In this paper, an overall assessment method has been developed in order to improve the understanding of this phenomenon. A full deterministic approach has been completed ®rst. An automatic numerical process, using the dual boundary element method (DBEM) and fracture mechanics laws, computes the crack evolution until failure of the structure. Then, in order to take materials scatter into account, a probabilistic approach has been developed based on Monte-Carlo simulations. Initial cracks scenarios are randomly de®ned and cracks evolution is computed for each of them. Thus statistical results are provided.

2. Dual boundary integral equations 2.1. Formulation The dual equations, on which the DBEM is based, are the displacement and the traction boundary integral equations (2) and (3). In the absence of body forces, the boundary integral representation of the displacement components ui , at a boundary point P, not on crack surfaces, is given by Z Z Tij …P ; Q†uj …Q† ds…Q† ˆ Uij …P ; Q†tj …Q† ds…Q† …1† Cij …P †uj …P † ‡ oX

oX

where oX indicates the problem boundary (including the crack faces), i and j denote Cartesian components, Tij …P ; Q† and Uij …P ; Q† represent the Kelvin traction and displacement fundamental solutions, respectively, at a boundary point Q [2,3,10].

H. Kebir et al. / Engineering Fracture Mechanics 68 (2001) 1371±1384

1373

The stress components rij are obtained by di€erentiation of Eq. (1), followed by the application of Hooke's law; they are given by Z Z 1 r …P † ‡ Sijk …P ; Q†uk …Q† ds …Q† ˆ Dijk …P ; Q†tk …Q† ds …Q† …2† 2 ij oX

oX

and the traction components, tj , are given by Z Z 1 t …P † ‡ n …P † S …P ; Q†u …Q† ds …Q† ˆ n …P † i ijk k i 2 j oX

oX

Dijk …P ; Q†tk …Q† ds …Q†

…3†

where ni denotes the ith component of the unit outward normal to the boundary, at point P. If P is on the crack surfaces (P ˆ P ‡ or P ˆ P ), the displacement equation becomes: Z Z ‡ ‡ 1 1 u …P † ‡ u …P † T …P ; Q†u …Q† ds …Q† ˆ Uij …P ‡ ; Q†tj …Q† ds …Q† …4† ij j 2 i 2 i oX

oX

and the traction equation becomes: Z Z ‡ ‡ ‡ ‡ 1 1 t …P † ‡ t …P † ‡ n …P † S …P ; Q†u …Q† ds …Q† ˆ n …P † i i i ijk k i 2 2 oX

oX

Dijk …P ‡ ; Q†tk …Q† ds

…5†

p The crack-tip is modeled by singular elements that exactly represent the strain ®eld singularity 1= r (Fig. 2). The local coordinates of the nodes are n ˆ 2=3, 0 and 2=3. The shape functions of this element are: p p 2 …3 15†n ‡ 2 1 ‡ n 2 p p N1 …n† ˆ 3 15 ‡ 3 6 p p p p p 3… 15 3†n 12 1 ‡ n ‡ … 15 ‡ 3† p p N2 …e† ˆ …6† 2… 15 ‡ 3 6† p p 2 … 3 3†n ‡ 2 1 ‡ n 2 p p N3 …e† ˆ 3 15 ‡ 3 6 The improper integrals, that arise in the dual integral equations, are handled analytically [5].

Fig. 2. Crack-tip elements.

1374

H. Kebir et al. / Engineering Fracture Mechanics 68 (2001) 1371±1384

2.2. The stress intensity factor evaluation There are several approaches to compute the stress intensity factors (SIF) using a dual boundary element formulation. In the present work, they are computed from the crack opening displacements at nodes which are extremely close to the tip. # p r " l p 3 5 N3 M3 N2 M2 …u2 KI ˆ u2 † …7† 5…u2 u2 † k‡1 l 5 and l KII ˆ k‡1

r " p 5…uN1 2 l

2 uM 1 †

p 3 5 N3 …u1 5

where l is the shear modulus and k ˆ 3 where m is the Poisson's ratio.

# 3 uM 1 †

…8†

4g; for plane strain g ˆ m and for plane stress g ˆ m=…m ‡ 1†,

Fig. 3. W ohler curve (S±N) for crack initiation.

Fig. 4. Touching plastic zones as a criterion for ligament failure.

H. Kebir et al. / Engineering Fracture Mechanics 68 (2001) 1371±1384

1375

We have shown in Ref. [5] that this approach gave very good results and it is not necessary to use the Rice integral method.

Fig. 5. Experimental and numerical life expectancy of a plane plate with 14 free holes.

Fig. 6. Simulation process.

1376

H. Kebir et al. / Engineering Fracture Mechanics 68 (2001) 1371±1384

2.3. Crack modeling strategy The general modeling strategy can be summarized as follows [4]: · all boundaries are modeled with discontinuous quadratic elements, except the crack-tip which is modeled by a singular element; · the displacement equation (4) is applied to collocation nodes on one of the crack surfaces; · the traction equation (5) is applied to collocation nodes on the other crack surface; · the displacement equation (1) is applied to collocation nodes on all non-crack boundaries; · the SIF are computed from the crack opening displacement at collocation points of the singular element; · the new crack-extension increment is modeled with new singular boundary elements. They will generate new equations and update the ones already existing with new unknowns. 3. Simulation model There are di€erent approaches to take into account the probabilistic character of the MSD problem. However, in the present paper, a mixed model is used, i.e. a probabilistic approach is used for the determination of a starting fatigue life scenario; while the subsequent steps (e.g. the damage accumulation, crack growth and crack coalescence) are calculated by a deterministic model. 3.1. Crack initiation In this model, it is assumed that the mean value of the fatigue life for crack initiation and the scatter factor (the standard deviation) for a single point are known (W ohler curve, Fig. 3). For the aluminum 2024 T3 alloy the mean value of the stress load initiation cycles is  N ˆ 10

5

S Slim FQI Slim

P …9†

Fig. 7. First scenario of Monte-Carlo simulations.

H. Kebir et al. / Engineering Fracture Mechanics 68 (2001) 1371±1384

1377

where FQI ˆ 176 MPa (fatigue quality index: N …S ˆ FQI† ˆ 105 cycles) [9]. Slim ˆ 59 MPa (fatigue limit) p ˆ 2:28 and the standard deviation is rˆ

a b ‡ ‡c S2 S

…10†

where: a ˆ 2169:9, b ˆ 1:2987, c ˆ 0:0279. The fatigue life is distributed according to a log±normal distribution, this may be done by means of an ordinary random number in the interval [0;1], and by inversion of the log±normal distribution by an approximate method. Once the fatigue life Ni of all fatigue critical locations ``i'' (i ˆ 1; M) has been determined, we compute the damage accumulation ``Di '' in each location ``i ''as: Di ˆ

Nk Ni

…11†

where: Nk ˆ Min…Ni † …i ˆ 1; M†.

Fig. 8. Second scenario of Monte-Carlo simulations.

Fig. 9. Rectangular plate with one hole.

1378

H. Kebir et al. / Engineering Fracture Mechanics 68 (2001) 1371±1384

Fig. 10. Statistical results of life expectancy.

It is assumed that the fatigue damage at the location ``k'' is high enough to start a crack growth while all other locations ``i '' (i ˆ 1; M and i 6ˆ k) have to accumulate more fatigue damage until the crack starts (Di ˆ 1). The origins of the new cracks are located at boundaries where the damage accumulation Di ˆ 1, the new cracks are given a small ®nite initial length a0 and a direction that is perpendicular to the boundary. They are then allowed to grow during subsequent loading steps. 3.2. Propagation Generally, the crack propagation is calculated by means of the SIF (KI and KII ) and the Paris equation: da ˆ C…DKeff †n ; …12† dN where da is the increment of crack growth in dN stress cycles. The parameters C and n are material constants (C 1:20110 6 ; n ˆ 3:211) and DKeff is the range ……Keff †max …Keff †min † of the SIF during these cycles qˆ  …Keff ˆ KI2 ‡ 2K2II † [6]. 3.3. Coalescence In this model, mainly the detection of the link-up of relatively small cracks is essential. This may be treated in di€erent ways. One of the very simple models for this problem seems to be the model proposed by Swift. It actually checks a criterion which is very similar to a ``net section yielding'' criterion [8].

H. Kebir et al. / Engineering Fracture Mechanics 68 (2001) 1371±1384

1379

The main idea is to use the radius of the plastic zone in front of each crack-tip   1 Keff Rp ˆ 2P r0:2

…13†

and to apply the connection of both plastic zones as a criterion for the link-up of the cracks (Fig. 4): Rp …a† ‡ Rp …b† ˆ L:

…14†

4. Algorithm The overall problem of the assessment of MSD in aircraft structures are approached by Monte-Carlo simulation. Initial crack scenarios are de®ned using a random process (Fig. 5). For each of them, a deterministic computation of crack growth until fracture is completed. Detection probability aspect is also taken into account in the whole process. Fig. 6 shows a diagram of the general approach implementation. We present here the overall structure of the model, where ``a'' indicates the crack length and D indicates the accumulated damage. · Determination of the M fatigue critical locations and the values of the stresses Si …i ˆ 1; M† · Loop over ``n'' scenarios ± Random run for the determination of the initiation number of load cycles Ni for each site i by using the W ohler curve and the value of Si . ± Compute the initial damage: Di ˆ Nk =Ni , where Nk ˆ Min…Ni † …i ˆ 1; M† ± Introduce a small crack a0 in the site k (in this simulation a0 ˆ 1 mm)  loop over the load cycles ``l''  loop over all sites ``i '', if the location is still in the damage accumulation phase (Di < 1) then: ± compute the new stress Si0 using the DBEM ± Di ˆ Di ‡ DN Ni0 (where Ni0 is computed by using the W ohler curve and the value of Si0 ) if Di ˆ 1 then: introduce a new small crack a0 in the site ``i '' endif else if the location is an active crack then: ± calculation of the SIF Keff n ± Dai ˆ C…DKeff † DN ± ai ˆ ai ‡ Dai ± calculation of plastic zone size ± check link-up criterion endif  check overall residual strength of the structure  calculation of the life expectancy for a given con®guration · calculate statical values for the problem: mean value, standard deviation, etc. To illustrate our approach, we consider the example of a rectangular plate with two holes. We present the two ®rst scenarios of the Monte-Carlo simulation in Figs. 7 and 8.

1380

H. Kebir et al. / Engineering Fracture Mechanics 68 (2001) 1371±1384

5. Numerical results 5.1. Rectangular plate with one hole Consider a rectangular plate with one hole as shown in Fig. 9. Initiation and crack growth tests were carried out at Aerospatiale Matra research center under R ˆ 0:1 constant amplitude loading at a frequency of 15 Hz (Smax ˆ 100 MPa). The numerical results were collected for 200 di€erent con®gurations, the mean value and the standard deviation of life expectancy were compared with the experimental one. We note that the results given by this model are in good agreement with the experimental values (Fig. 10). 5.2. Rectangular plate with 14 free holes In order to validate the global probabilistic approach, fatigue tests on a plane plate with 14 free holes were achieved at the Aerospatiale-Matra laboratory in Suresnes (France) (Fig. 11). The used material is an

Fig. 11. Rectangular plate with 14 free holes.

H. Kebir et al. / Engineering Fracture Mechanics 68 (2001) 1371±1384

1381

Table 1 Mechanical properties of the material E (GPa) 72.7

R0:2 (MPa) 318

Rm (MPa) 478

A (%) 22.2

aluminum alloy 2024T3 sheet with a thickness of 1.6 mm. The load was applied on the transversal direction. The properties of the material are presented in Table 1. The tests are done under imposed load with a constant amplitude of a ratio R ˆ 0:1, and with a maximum stress level of 100 MPa. Six specimens were tested. The experimental results are expressed in Fig. 12. Fig. 13 illustrates the chronology of the major steps of the ®rst specimen evolution. The given life expectancies of initiation correspond to the detection of the ®rst crack of a length between 0.5 and 1 mm. The minimum detection length of cracks during the tests is around 0.3 mm. The evolution of the crack lengths with the number of loading cycles was measured using an optical microscope. 5.3. Numerical simulation The initial structure is discretized with 176 elements. The mesh and the boundary conditions are depicted in Fig. 14. The numerical progress step of the major crack is 2 mm, …Keff …majorcrack† ˆ max…Keff …ith crack††. i This value is a good compromise between the results accuracy and acceptable time of computation. For the other cracks, the progress step is proportional to the corresponding crack growth rate da=dN . To take into account the expansion process applied to external holes, we attributed to concerned sites a large potential life expectancy. A ®rst computation is carried out with the deterministic approach by using the W ohler's curve at 50%. The probabilistic approach is achieved with 200 con®gurations. The obtained results are presented in the graph (Fig. 15).

Fig. 12. Synthesis of the test results.

1382

H. Kebir et al. / Engineering Fracture Mechanics 68 (2001) 1371±1384

Fig. 13. Chronology of failure.

Fig. 14. Initial and deformed boundary element meshes.

With the deterministic approach, we note that all the cracks begin at the same time, since all the sites are undergoing the same stress level. This explains why the propagation phase is so short. It is clear from this results that the probabilistic approach gives a better estimate of the life expectancy in comparison with the deterministic approach. If we observe the results of the stochastic analysis in Fig. 15, we can see a slight correlation among the parameters Nbeginning (number of cycles corresponding to the initiation of the ®rst crack) and Npropagation (number of cycles from the initiation of the ®rst crack to the failure). We can explain this fact physically: if a ®rst initiation of a crack takes place so early, which is a priori a scare event, there is a high probability so that the crack is still insulated a certain time before that others appear. The propagation phase is long and

H. Kebir et al. / Engineering Fracture Mechanics 68 (2001) 1371±1384

1383

Fig. 15.

we are in the neighboring of a situation called ``mono-crack''. In the other side, more the ®rst initiation comes later, more the probability that the other cracks begin just after is high. The phase of propagation is shorter since we have a multi-crack situation.

6. Conclusion An overall assessment method has been developed in order to improve the understanding of MSD phenomenon. A full deterministic approach has been completed ®rst. An automatic numerical process, using the DBEM and fracture mechanics laws, computes the crack evolution until failure of the structure. Then, in order to take materials scatter into account, a probabilistic approach has been developed based on Monte-Carlo simulations. The results obtained are in good agreement with experimental and show that the probabilistic approach gives a better estimate of the life expectancy in comparison with the deterministic

1384

H. Kebir et al. / Engineering Fracture Mechanics 68 (2001) 1371±1384

approach. The work on this method will be continued in order to re®ne and extend the model to the bolted joints. References [1] Aircraft accident report: Aloha Airlines, Flight 243, Boeing 737-200, N73711, near Maui, Hawaii, 28 April, 1988. Washington: US national transportation safety board, 1989. [2] Brebbia CA, Dominguez J. Boundary elements: an introductory course. New York: Mc Graw-Hill; 1989. [3] Portela A, Aliabadi MH. The dual boundary element method: e€ective implementation for crack problems. Int J Numer Meth Engng 1992;33:1269±87. [4] Portela A, Aliabadi MH, Rooke DP. Dual boundary element incremental analysis of crack propagation. Comput Struct 1993;46:237±47. [5] Kebir H, Roelandt JM, Foulquier J. A new singular boundary element for crack problems: application to bolted joints. Engng Fract Mech 1999;V62:497±510. [6] Kebir H. Roelandt JM, Foulquier J. Dual boundary element incremental analysis of crack growth in bolted joints. MAT-TEC 96 Improvement of materials, 1996. p. 230±33. [7] Schijve J. Multiple-site damage in aircraft fuselage structures. Fatigue Fract Engng Mater Struct 1995:329±44. [8] Jeong DY, Brewer JC. On the linkup of multiple cracks. Engng Fract Mech 1995;51(2):233±8. [9] Foulquier J. Etude du comportement de multi®ssures sous sollicitations cycliques et statiques. Report, Centre Commun de Recherches Aerospatiale, 1996. [10] Aliabadi MH. A new generation of boundary element methods in fracture. Int J Fract 1997;86:91±125.