Multiple crack propagation by DBEM in a riveted butt-joint : a simplified bidimensional approach

A Multi-Site Damage (MSD) crack growth simulation is presented, carried out by means of Dual Boundary Element Method (DBEM), in a two-dimensional analysis of a cracked butt-joint made of aluminium 2024 T3. An equivalent crack length is proposed for an approximated 2D analysis of a 3D problem where the crack front assumes a part elliptical shape due to secondary bending effects. The assumptions made to perform such simplified bidimensional analyses are validated by comparing numerical results with experimental data, the latter obtained from a fatigue tested riveted butt-joint.


INTRODUCTION
omplex engineering structures will often present microscopic flaws that can be caused by the manufacturing process or by fatigue, impact or environmental effects such as corrosion.It is therefore crucial to be able to predict how these cracks will affect the integrity of the structure, for example to determine whether a particular crack will grow at all under service loading and, in case, what the life of the component before failure.It is evident from this considerations that the field of Fracture Mechanics will be an integral part of the design process, particularly in the aircraft industry due to the high strength but low crack resistance materials used in weight critical applications.Experimental studies have shown that, if the hypothesis of Linear Elastic Fracture Mechanics (LEFM) holds true, failure occurs when the Stress Intensity Factor (SIF) reaches a critical value.Hence the predictions of crack behaviour and the integrity of a structure are dependent on the accurate calculation of SIFs.

C
In this work a Multi Site Damage (MSD) crack growth simulation is presented, carried out by means of Dual Boundary Element Method (DBEM) [1][2][3][4] as implemented in the commercial code BEASY [5], in a two-dimensional analysis of a butt-joint, made of aluminium alloy 2024 T3, undergoing a traction fatigue load.The bidimensional model represents an approximation of the real phenomena because of the secondary bending effects, illustrated in Fig. 1 with reference to a lap-joint (but the phenomenon is completely equivalent for a butt-joint) and responsible for a non-straight crack propagating front.Such phenomena would prevent a bidimensional modelling approach, but this drawback can be circumvented, in a preliminary analysis, by adopting an equivalent straight initial crack front as explained in the following.

J INTEGRAL
he J integral was developed by Rice and Cherepanov to characterize fracture for two dimensional configuration.The material is elastic or non-linear elastic.The J integral is defined as: where W is the strain energy density, ti are components of the traction vector and ui are components of the displacement vector.The vector n is the unit normal to the contour  .The strain energy is defined as 0 ( )  is the infinitesimal strain tensor component.The J integral is integrated along a contour  surrounding the crack tip (the ends of the  contour are on opposite faces of the crack and the contour encloses the crack tip).The crack is parallel to the x 1 axis.On addition or subtraction of quantities above and below the crack axis, symmetric and anti-symmetric stress and strain derivatives are found from which symmetric stress and strain products are calculated respectively.These products can be decomposed into two parts: one consisting of symmetric stress and strain derivatives (symmetric integrands), the other consisting of anti-symmetric stress and strain derivatives (anti-symmetric integrands).
The symmetric integrands can be integrated over the contour area to obtain the mode I J-integral; the anti-symmetric integrands, when integrated, will produce the mode II J-integral [6].

FATIGUE CRACK GROWTH
he well-known Paris law is adopted for the crack growth rate assessment: with the calibration parameter C and n for Al 2024 T3 given by n=2.144 and C=4.72*10 -10 [7], (Tanaka formula [1]) expressed in MPamm and da in mm.
It has been checked that, during the propagation, the SIF's range is included in the interval of validity of Paris law, skipping the very initial crack growth phase because of the non-modelled threshold phenomenon and the final crack growth phase, which is dominated by very high growth rates, with Kmax approaching Kc (Kc is the material fracture toughness).

PROBLEM DESCRIPTION AND PRELIMINARY RESULTS ON UNDAMAGED PANEL
n the following a Multi-Site damage (MSD) crack growth simulation for a butt-joint is presented, carried out by DBEM in a two-dimensional analysis.Numerical results are then compared with corresponding experimental outcomes, obtained from a fatigue tested riveted butt-joint specimen of aluminium 2024 T3 [7].

T T I
The modelled specimen is representative of the riveted joint existing between the upper shell and the lower shell on the front fuselage section, manufactured by AerMacchi for the Dornier Do328 aircraft (Fig. 2a).
The corresponding DBEM model is illustrated in Fig. 2b, where internal points are placed in correspondence of the strain gauge locations (Fig. 2a) in order to compare the numerical and experimental strains and assess the model accuracy.
The distance from the panel edge to the first hole was chosen from the experimenter equal to 10 mm, giving a total panel width W=280 mm (Fig. 2a).The four holes close to the plate border were cold worked in order to avoid crack initiation during the fatigue test (in Fig. 2a, holes with crack tips N. 1, 28, 29, 56) or at least strongly retard crack propagation [8].
The remote applied load is t=100 MPa (Fig. 2b) and the material properties are: Young modulus E=72500 MPa and Poisson ratio .33   .The butt-joint is designed and loaded in the experimental tests, in such a way to reproduce, as close as possible, the real in-service stress state, like, for example, a uniform distribution of longitudinal stresses  yy along a transversal section of the undamaged joint sufficiently far from the riveted area.Along such section a row of strain gauges was placed and correspondingly a row of internal points was introduced in the DBEM model (N.6:17 in Fig. 2b) in order to monitor the level of aforementioned stresses.The same was done to cross check the longitudinal stress distribution in a longitudinal section (internal points N. 1-5 in Fig. 2b).In particular, for the undamaged panel, the DBEM model provides the longitudinal stress distribution at strain gauge locations illustrated in Figs.3a-b: the numerical stresses are in good agreement with the experimental stress state [7].It is evident that, whenever experimental strains were available from both sides of the panel they should had been averaged in order to be comparable with numerical results (strain values are slightly different on the two panel sides because of the secondary bending).This check turned out very useful in order to calibrate the model and in particular to decide how many rivets to be explicitly modelled by inserting a pin in the corresponding hole: such pin would be constrained against y translation and disconnected from the hole upper surface by means of internal spring of negligible stiffness (Fig. 2b).Alternatively, the holes not directly involved by crack propagation can be modelled by just introducing longitudinal constraints on 180 degrees of the hole boundary and skipping the explicit pin modelling (Fig. 2b).Gap elements have also been introduced, to better tackle contact conditions but the solution improvement has been judged quite negligible (less than 2%), except in case of very short cracks initiated from the holes, more sensitive to pinhole contact conditions (but not present in this work).For this reason, and due to the computational effort of a non-linear analysis, they have not been used anymore.
The J-integral technique is adopted for SIF's evaluation, being more stable than Crack Opening Displacement method, against a variable refinement of crack mesh.On the J-integral path, 33 integration points are used (the increment of accuracy with 66 points turn out to be negligible).
The mesh used for the butt-joint is based on about 327 quadratic elements: a p-convergence study has been realized showing that cubic elements provide an accuracy improvement of less than 2% and that 2 quadratic elements per 90 degrees are sufficient on the cracked hole, except for very short cracks, where 3 elements are recommendable (possibly with a scaling ratio).All the undamaged holes are modelled with 6 elements, constrained in y-direction.

NUMERICAL CRACK PROPAGATION
he crack growth has been simulated considering the interval ranging between 73500 and 113885 fatigue cycles.
The initial crack length in the model, at left or right hole sides, is not taken equal to that (size b in Fig. 4a) visually observed on the free surfaces at 73500 cycles [7] because of the following reasons.The crack length monitored in the experimental phase is related to the external surface ("penetrated crack"), but the crack first appear on the faying surface because here undergoes the primary tensile stress and, in addition, the secondary bending stress (Fig. 1).As a matter of fact, due to the superimposed secondary bending the crack assumes a part elliptical shape (Fig. 4), in such a way that it is already abundantly extended on the faying surface when appearing on the external surface.It is possible to make a forecast of the hidden part crack length, right before the appearance of the external crack front (Fig. 4b), observing that the ratios between the two ellipse semi-axis, crack depth a and crack length c and between a and t (specimen thickness), even if variable at initiation, assume always the same values for a given bending factor k=max tension/max bending when the crack is on the verge to become through the thickness (e.g. for k=1 we obtain a/c=0.575and a/t=0.88)[9].In our particular case, for a plate thickness t=1.2 mm, when each crack appears on the visible plate side with a length b (as reported by the experimenter), by using the aforementioned relationships, it is possible to know the size of the part elliptical crack on the faying side, obtaining c=a/0.575=0.88*t/0.575=1.84mm (Figs.4a-b).Keeping into account these data and the shape of the real hole (Fig. 4a) against the simplified modelled geometry (a cylinder of radius 2.4 mm), when modelling an equivalent straight initial crack we have prolonged of 1 mm the initial crack length value b measured by the experimenter.
During the propagation some cracks will link-up and the criteria chosen to assess such condition is the overlapping of the crack tip plastic zones: L=rp1+rp2 where L is the residual ligament and rp=(Keq 2 /y 2 )/; y is the yield stress.When two cracks link-up a new crack initiation is postulated from the hole side opposite to the crack.In Tab. 1 and Tab. 2 the crack length vs. cycles and the related equivalent SIFs are respectively shown.With the aforementioned bidimensional approach a good agreement between numerical and experimental [7] crack scenarios is obtained and it is possible to capture the essential features of the response even if the secondary bending is not modelled.

CONCLUSIONS
t is important to point out the need, in a bidimensional crack propagation model for a butt-joint, to model the part elliptical crack front at initiation with an "equivalent" prolonged straight crack.As a matter of fact, in such case the penetrated front experiences a lower stress, compared with the hidden front surface and consequently there is not the "catch up" behaviour that is typical of the pure tensile case, where the "penetrated" crack front reaches immediately the more prominent front surface because of the higher SIFs [9].Even with the illustrated simplified bidimensional approach it is possible, for such kind of problems, to obtain a satisfactory agreement with experimental crack growth rates.Moreover, very short run times are needed to run the whole propagation and an easy preprocessing phase is enabled by the DBEM approach: an automatic remeshing is possible as the crack grows and the manual intervention is just necessary to initiate new cracks.A more accurate DBEM bidimensional approach (with respect to the simplified approach adopted in this work) to model the butt joint assembly is also possible when needed: each layer can be considered as an individual two-dimensional structure; individual layers can be explicitly modelled and connected with rivets; by gap elements can be used at the interface pin-hole… [10].The very basic approach presented here, with related very short run times, becomes mandatory in case of a probabilistic approach to crack propagation simulation where hundred thousands of such simulations are to be performed (e.g. when resorting to Monte Carlo method…) [11].

Figure 1 :
Deformed lap-joint undergoing a traction load with highlight of sites undergoing maximum bending stresses (a) and stress distribution through the thickness of the sheet (b).

Figure 2 :
Figure 2: Butt joint sketch (a) with highlight of strain gauge positions and crack tip numbering (N.1:28 in the lower hole row and N. 29:56 in the upper cracked hole row); DBEM uncracked model (b) with highlight of internal spring (yellow symbols), internal points and hole constraints (green symbols).

Figure 3 :
Figure 3: Membrane longitudinal stresses  yy on transversal section at y=54.25 mm (a); longitudinal stresses  yy on longitudinal section at x=110 mm (b).

T
Intermediate deformed plots, during the multi side crack propagation, are presented in Figs.5a-g, where each new crack introduced in the model comes from the experimental recordings or is postulated to appear when two cracks link-up (as previously said). in the hole to model the action of the overlapped plate (no pin modelling)

Figure 5 :
Figure 5: Deformed shape (scale factor 20) of multiple crack scenario after a number of cycles equal to: 73500 (a), with initial cracks (tip N. 41 and 42) as recorded by the experimenter; 84000 (b) with a new recorded initiation (tip N. 43); 99500 (c) with propagation of the previously indicated cracks and a new recorded initiation (tip N. 47); 105770 (d) with a link-up between tips N. 42 and 43 and corresponding new postulated initiation (tip N. 44); 108000 (e) with new recorded initiations (tip N. 39, 40 and 45); 110500 (f) with a link-up between tips N. 40 and 41 and a new recorded initiations (Tip N. 51); 112000 (g) with new recorded initiations (Tip N. 38).