Stress Intensity Factors Evaluation at Tips of Multi-Site Cracks in Unstiffened 2024-T 3 Aluminium Panel Using XFEM

The aim of this paper was to establish and demonstrate significant capacity and performances of extended finite element method (XFEM) to calculate stress intensity factors (SIFs) histories versus crack length for problems involving multiple, interacting cracks resulting from multiple site damage (MSD). A typical aero structural configuration was analysed: unstiffened flat panel made of aluminium Al 2024-T3, containing 11 holes, each of which is a site for crack growth. Analysed model makes a unique 3D configuration with 22 cracks, propagating at the same time. Short theoretical background information is provided on the XFEM, as well as representation of cracks and SIFs computation. The computations were carried out in Morfeo/Crack for Abaqus software which relies on the implementation of XFEM. The accuracy of these computations was verified through FRANC2D/L software and superposition based approximate method. The conducted analysis showed that XFEM is efficient tool for the simulation of crack propagation even in the case of 3D configurations with MSD, and that the obtained solutions can be used for the prediction of SIFs in analysed MSD configuration with acceptable accuracy.


INTRODUCTION
Aircrafts should sustain damage up to the limit load for the whole service life, avoiding catastrophic failure at any cost.This is a design philosophy, which sometimes fails due to real life conditions.Namely, competition sometimes forces airlines to use aircrafts even after their original design service.By doing so, they are raising possibility of fatigue crack growth and widespread fatigue damage (WFD), including multiple site damage (MSD), as its common form.The Aloha accident in 1988 initially focussed the attention of fracture mechanics research community to this phenomenon, [1].Multiple site damage is the interaction of a major crack with several short cracks located at various sites of the same component.It presents a typical problem for ageing aircraft, often occurring in longitudinal and circumferential riveted lap joints.When cracks initiate and grow under fuselage pressure cycling fatigue loads.When critical conditions are met, MSD cracks interact and sudden crack link-up may occur reducing significantly the overall structural integrity, [2].So, prediction of the crack growth rate and residual strength of the structure with MSD require an accurate calculation of stress intensity factors (SIFs).
Many studies and analyses of this phenomenon have been carried out in recent decades.Overviews on analytical methods for MSD are provided in many scientific papers, such as [3].As technology and computer sciences were developing and became widely available, numerical analysis started to be regularly used for SIFs assessments.The use of these analyses came as a demand, since it was required to conduct various tests during exploitation to detect fatigue damage resulting from repeated working load and to evaluate residual strength and structure integrity of the structural elements.There are many examples of different numerical techniques and methods used for this kind of analysis, and we will mention here just a few.D. Partl and J. Shijve carried out a simple computational procedure using the principle of fracture mechanics and compounding K-solution method to predicting fatigue crack growth on 2024-T3 aluminium alloy sheet specimens with a collinear row of three central holes with different geometrical configurations [4].R.
Elangovan et al. studied the effect of interaction between cracks using numerical estimation of stress intensity factors at a crack tip in unstiffened curved panel with secondary crack in the vicinity of a primary crack.The results are presented in the form of design charts [5].G. Kastratovic et al. proposed approximation method for the determination of stress intensity factors (SIF) in case of MSD based on existing solution for SIF in case of two unequal cracks in infinite plate, subjected to remote uniform stress, [6].A research on a finite element calculation of SIF in structures with MSD using approximate procedure was published in [7], where results were compared against values obtained by FRANC2D software with acceptable accuracy.
It should be mentioned that the finite element method (FEM), the most popular numerical method nowadays, was also often used for crack modelling and propagation investigations.But, over time the restrictions in FEM crack propagation simulations discouraged researchers who then used numerous variations of FEM, [8,9], until the extended finite element method (XFEM) was developed for modelling discontinuities in 1D, 2D and 3D domains.In [10] authors used XFEM for analysing models containing several cracks with voids, cracks with multiple branches, and cracks emanating from the holes.J. H. Kim et al investigate the effect of compression stresses, stress level and stress order on fatigue crack growth of multiple site damage.In this investigation, XFEM was applied to predict lifetime under constant amplitude cyclic loading on several MSD specimens made of Al 2024-T3.Based on this, the multiple crack growths under service stress spectra are calculated to assess effects of compressive stress, stress order and sequence cyclic loading on stress level by using Forman and NASGROW equations [11].
In many studies, [12,13,14], XFEM was used to analyse various effects of multiple micro-cracks on a macrocrack.So-called homogenized XFEM has been proposed by S. Kumar et al. [15], for the evaluation of fatigue life of an edge crack plate in the presence of multiple discontinuities.To improve the results, a homogenization scheme based on strain energy density approach is used.All mentioned studies and research activities, regarding XFEM, were concentrated on its implementation on relatively simple configurations, with several cracks in 2D models of panels and sheets, where genuine MSD was not really an issue.Also XFEM itself in these studies was carried out through specifically made "custom code" numerical calculations, which make their application hard and reduced to a very small number of potential users.It should be mention that there are also several studies where XFEM was used for simulating 3D fatigue crack propagations, but only with one crack [16,17].
In this study, the XFEM is used to calculate the SIF versus crack length histories in a typical aero structural configuration: unstiffened flat panel made of aluminium Al 2024-T3, containing 11 holes, i.e. sites for crack growth.The analysed model makes a unique 3D configuration with 22 cracks, propagating at the same time, whereas stress intensity factors were computed along the crack front for each of 22 cracks.This kind of 3D multisite damage model with this number of cracks emanating from the rivet holes, and propagating at the same time until their link-up, was never analysed in the literature that was available to the authors.The computations were carried out in Morfeo/Crack for Abaqus software which relies on the implementation of XFEM.The procedure accuracy has been verified by comparison with solutions obtained by using FRANC2D software and approximate method, [6].

THE ADVANTAGES OF EXTENDED FINITE ELEMENT METHOD
The finite element method (FEM) has been used for decades to provide approximate solutions of partial differential equations.Anyhow, standard FEM is based on the approximation of polynomials and becomes expensive if the solution contains discontinuities in the displacement field.The non-smooth displacement near the crack tip is basically captured by refining the mesh locally which may lead to the significant increase in the number of degrees of freedom in 3D applications.Furthermore, in crack propagation simulations the finite element mesh needs to be updated after each propagation step in order to track the crack path.
Therefore, development of a partition of unity (PU) based enrichment method for discontinuous fields, using special functions, referred to as the XFEM, was a significant improvement in crack modelling, [18].For crack modelling, a discontinuous function such as the Heaviside step function and the two-dimensional linear elastic asymptotic crack tip displacement fields are used.No mesh updates are needed and domain can be modelled by finite elements without explicitly meshing the crack surfaces.The position of crack discontinuity can be anywhere, while the finite element framework and its properties are retained.This makes the XFEM a powerful tool for cracks and dislocations where near-field solutions can be embedded by the PU method to increase significantly the accuracy of relatively coarse meshes.The technique enables efficient analysis of phenomena such as grain boundaries, phase interfaces, surface effects in nanomechanics and voids.Recently, XFEM and its coupling with level set method (LSM) were intensively studied.LSM is used to define the location of non-smooth features and it complements the XFEM extremely well providing the information where and how to enrich.
The XFEM method has greatly enhanced the power of the FEM, so many attempts have been made to incorporate the modelling of discontinuities independent of the FE mesh by either a plug-in or native support.Cenaero [19] has developed a crack growth sumilation add-in Morfeo/Crack, based on the XFEM method in Abaqus software.Problems with static cracks, evolving cracks and cracks emanating from voids were numerically studied in order to demonstrate the robustness of the XFEM and precision of Morfeo/Crack for Abaqus [20].

THE ANALYSIS OF MSD CONFIGURATION USING XFEM
In  The FE model was based on the configuration of panels tested by Luzar in 1997 [21].The 2024-T3 clad MSD panels were made and tested to determine fatigue crack growth values of airframe structure with MSD.The panel thickness, fastener hole diameter and pitch were selected as for a real aircraft lap joint configuration.The test panels were held by the grip end fixtures resulting in all of the test loads being transmitted through joint friction.The loading was arranged to produce uniform stresses and displacements throughout the test section.Following properties for linear elastic behaviour of 2024-T3 aluminum alloy are used: Young's modulus 73.000 MPa, Poisson's ratio 0.33.Uniform tensile stresses were applied at the top and bottom edges of the model (Fig. 2).Each hole in the panel had two radial cracks, numbered from 1 to 22 in Fig. 1 and positioned as shown in Fig. 3. Fig. 3 also shows a portion of the mesh around the holes, and as it can be seen the mesh was refined around the cracks at the edges of each hole and a uniform template of elements around each crack front was used.Final mesh consisted of 163228 linear hexahedral elements of type C3D8R and 206800 nodes.This mesh was used to calculate SIFs at crack nodes as a function of crack length.The crack growth simulation capability of Morfeo/Crack for Abaqus was used for this purpose.Morfeo/Crack calls Abaqus at each step and also between the steps and reads the Abaqus solution, recovers a richer, improved XFEM solution in a small area surrounding the crack and computes the SIFs at crack front nodes [22].Finally, it defines the appropriate crack growth increment, extends each crack and then performs the next solution step.Results are presented in Fig. 4-6 (step 1, 14 and 40 resp.).The initial crack length was 0.5 mm.
As mentioned before, SIFs were calculated for each crack front and different crack sizes.Maximum values of SIFs calculated along the crack fronts were used as a reference.The SIF results shown in Fig. 7 represent the solutions for cracks 1, 2, 11, 12, 21 and 22 only.These six cracks were selected because of their unique positions.That is, at these positions the influences of adjacent cracks interaction are either minimal (at the first and the eleventh hole; cracks 1, 2, 21 and 22 respectively) or maximal (the fifth hole; cracks 11 and 12).
As can be seen in Fig. 7 the highest values of SIFs are for the cracks 11 and 12, that is, for the cracks emanating from the fifth hole.This becomes more obvious with cracks' growth.Also, this occurs for all three load cases.The values of SIFs for cracks at the first and eleventh hole are smaller and similar.This can be explained by the fact that the crack interaction effect is more dominant for the cracks at middle holes, than for the cracks at "edge" holes.That is precisely why the cracks 1, 2, 11, 12, 21 and 22 are selected for results presentation.

APPROXIMATE METHOD FOR SIFs CALCULATION IN A CASE OF MULTIPLE CRACKS
The XFEM results obtained in these simulations were compared against the results calculated by approximate method, based on principle of superposition, presented in [6].According to this procedure, the SIF for opening mode in a configuration with n cracks (Fig. 8) can be estimated as: where K IjA, B is the SIF for tip A, or B of analysed crack in presence of all other cracks in configuration; K Ii is individual SIF for all cracks in configuration, i.e., stress intensity factors of auxiliary configurations; c ib,d is the coefficient representing effect of the i th crack on SIF of analysed crack, (the influential coefficient of the analysed crack on itself is c jb, d = 1).

Figure 8 Configuration with multiple cracks
Therefore, analysed complex 2D or 3D configuration can be represented as a combination of simpler (auxiliary) configurations, usually containing only one crack.In that way, the calculation of SIF at analysed crack tip is based on effect of every crack in the initial configuration on the analysed one.The influential coefficients were estimated for vast number of crack lengths and distances between them, i.e. for their combinations, in [23].Here, eleven auxiliary configurations were used, which were all the same: thin plate with central circular hole with two radial cracks subjected to uniform uniaxial tensile stress, with well-known Bowie's solution for stress intensity factor [24].One should notice that the SIFs were calculated for models with different crack sizes for all the cracks in the configuration, but with same crack increment for all cracks, since all cracks in MSD are approximately of the same length (so called "catch-up" phenomenon), [25].
The SIFs for analysed configuration were also calculated in FRANC2D/L software, [26], which has the crack growth simulation capability.Contrary to Abaqus, FRANC2D/L calculates the SIF values at the tip of 2D crack using the J-integral, determines the appropriate crack growth increment, extends the crack, re-meshes around new crack tip, and then performs next solution step with new FE mesh.This procedure was performed 15 times in order to simulate incremental crack growth, using the 1/4point singular elements.

ANALYSIS OF RESULTS
In order to verify the results obtained by means of XFEM, SIFs calculated in Morfeo/Crack for Abaqus are compared with results obtained using classical FEM (FRANC2D/L software) and previously described approximate method.The results are presented through normalized SIF (geometry factor π Comparing the results obtained by XFEM (Morfeo/Crack for Abaqus) and FEM (FRANC2D/L) it can be seen that SIFs histories for selected cracks differ up to 14.7% (in the case of crack 21).The maximum differences occur for initial crack size, due to the proximity of the hole boundary.For the cracks on the first two holes the difference between XFEM and FEM results decreases with cracks growth to 6.1% (crack 1), i.e. 4.1% (crack 2).For cracks at fifth hole this difference starts to decrease with cracks growth to 0.5% (crack 11), i.e. 1% (crack 12), but with further cracks growth it increases up to 6.27% (crack 11), i.e. 8.1% (crack 12).The cracks 21 and 22 at eleventh hole exhibit the similar trend.In all of these cases XFEM gives slightly higher results than FEM.This is due to the initial crack size being less than the thickness of the plate in the case of FRANC2D/L, whilst in the case of Morfeo/Crack 3D solid model was used.The trend of increasing differences with cracks growth may be explained by the different manner in which these two methods implement mutual cracks interaction.Anyway, the difference can probably be reduced by increasing the number of nodes in the regions among the cracks on 3D model; nevertheless, the obtained solutions for SIFs are acceptable from an engineering point of view.
Comparing the XFEM results with the results calculated by means of approximate method better agreement can be seen (Fig. 9).SIFs histories for selected cracks differ up to 10.75% (crack 22), also for initial crack size.But, with the cracks' growth, these differences decrease, and for the cracks' length of 2 mm they are well under 4%.Since used approximation method defines the influential coefficients that take into consideration cracks interaction effect for every crack in the analysed structure, this leads to conclusion that XFEM also takes cracks interaction effect into consideration in a very good manner.

CONCLUSION
The prediction of crack growth rate and residual strength of cracked structure demands accurate calculation of stress intensity factors.There are continuous researches that deal with the problem of SIFs determination.Some solutions are available, but only for simple geometry configurations with few cracks.As technology and computer sciences develop and become more available, the researchers are trying to introduce and apply new computational methods and techniques in order to find solutions for complex geometries with multiple cracks.This study represents an effort in that direction.
Here, the SIFs calculations based on implementation of XFEM in Abaqus, are conducted for a typical aero structural configuration with MSD.Analysed model is a unique 3D configuration with 22 cracks that propagate at the same time, whereas stress intensity factors are computed along the crack fronts for all 22 cracks.
As it can be seen in Figs. 4 to 6 cracks spread all the time through unchanged mesh and can reach long lengths.With XFEM the need for new mesh creation after each step is eliminated.
In order to check the significance of obtained values, SIFs were also obtained by means of FRANC2D/L software, and by superposition based approximate method.The analysis of the results has shown that the values obtained by means of XFEM can be used for the evaluation of fatigue life of analysed configuration with acceptable accuracy.
It is interesting to note that -when compared with FRANC2D/L values -there is a trend of increasing differences between obtained SIFs as cracks grow, which may be related to the differences in methods for simulating crack interaction.Also, FRANC2D/L analysis was carried out on the two-dimensional structure.This kind of analysis does not take into account the effect of three-dimensional stress state in unstiffened panel which is actual 3D structure.
Although, XFEM is time consuming and requires substantial computer resources (especially when multiple cracks are involved), this paper showed that with a well-defined mesh, and well-set boundary conditions 3D simulations of multiple cracks' growth give good results.It also showed that for three-dimensional cracked structures, this kind of analysis can be fully acknowledged.This is of great importance for fatigue life estimation of aircraft structures, such as fuselage or wing with riveted joints, which are very often exposed to multiple site damage.

Figure 1 Figure 2
Figure 1 Analysed configuration with multiple cracks (not scaled)

Figure 3
Figure 3 Cracks' positions and numbering

Figure 4 Figure 5 Figure 6 Figure 7
Figure 4 XFEM model of MSD panel after 1st step of cracks propagation, σ = 200 MPa) shown in Fig.9.Factor β is a function of the SIF for mode I (K I ), stress σ and crack length a.

Figure 9
Figure 9 Comparison of SIF values obtained by three different methods (σ=200 MPa)