Application of X-FEM to 3D Real Cracks and Elastic-Plastic Fatigue Crack Growth

Summary. In a general point of view, X-FEM plus level set representation of the interfaces reveals to be of great interest in the aim to couple experimental data with numerical simulations. This can be highly illustrated in the case of 3D fatigue crack growth simulations where the initial 3D “real crack” is extracted from tomographic images. The experimentally observed fatigue crack propagation is compared to numerical simulations. Good agreement is found when a linear variation of closure stress along the crack front is taken into account in the “3D crack propagation law” used for the simulation. Furthermore, in order to take into account plasticity during fatigue crack growth, one develops an augmented Lagrangian formulation in the X-FEM framework that is able to deal with elastic-plastic crack growth with treatment of contact and friction. Numerical issues such as contact treatment and numerical integration are addressed, and ﬁnally numerical examples are shown to validate the method.


Introduction
It is now well established that the eXtended Finite Element Method is of great interest for evolving discontinuities, in particular for industrial applications [13,14].Indeed, no initial mesh and remeshing techniques are needed during the evolution of the interfaces [2,3,12].Furthermore, in a more general point of view, X-FEM plus level set representation of the interfaces reveals to be of great interest in the aim to couple experimental data with numerical simulations [1,2,12].This can be highly illustrated in the case of 3D fatigue crack growth simulations where the initial 3D "real crack" is extracted from tomographic images with a spatial resolution of the order of 1 µm [4].In such a case, X-FEM is coupled with a level set representation of the crack and a specific algorithm is developed in the aim to define the initial level sets from tomographic images.One can also notice that it is not necessary to impose the compatibility between the meshes of the structure and the support of the level sets.Furthermore, a robust and accurate technique is proposed to compute the stress intensity factors [1].Based on experimental observations and quantitative analysis of crack propagation in the bulk, a "3D" crack growth law can be established to predict the crack front shape evolution in the bulk of the material.From a more general point of view, this study shows that coupling X-ray microtomography to X-FEM provides a promising tool to assess the 3D behaviour of arbitrary shaped cracks and to perform direct comparisons of "experimental" and "simulated" crack shapes during propagation.
In a second example, it is shown that even with elastic-plastic behaviour coupled with contact and friction, non-remeshing property can be preserved for instance for 2D mixed-mode plastic fatigue crack growth [4,5].For that purpose, an augmented Lagrangian formulation in the X-FEM framework for the treatment of contact and friction with elastic-plastic behaviour is proposed.On the one hand, the numerical integration is adapted in order to properly integrate the high order terms in the enrichment basis, and to have a fine knowledge of the stress state around the tip to precisely model plasticity.On the other hand, the strategy of enrichment used for linear elastic X-FEM fatigue simulation has to be modified [5].In this respect, the integration strategy is designed in order to avoid the projection of stresses and internal variables as the crack evolves to ensure the reliability of the method.Numerical issues are addressed.In particular, one shows the ability of the method to model the phenomenon of crack closure under cyclic tension [6,9].Indeed, this can have a great influence on the propagation of the crack with fatigue loading.
2 Application of X-FEM to 3D Real Cracks

Experimental Observations
The propagation of a semi-elliptical crack in the bulk of an ultrafine-grained Al-Li alloy has been investigated using synchrotron radiation X-ray microtomography.In this material, the studied crack, despite its small dimension, can be considered as a "microstructurally long" and described in the frame of the linear elastic fracture mechanics.The main advantage of using this alloy for the present study is that it exhibits an exceptionally linear crack path compared to ingot metallurgy Al alloys.The ultrafine grains promote homogeneous deformation and prevent crystallographic cracking so that the crack shape is not disturbed by microstuctural features, at least not at the level of the spatial resolution employed in this study.In situ fatigue tests monitored by X-ray microtomography were carried out on beamline ID 191 of the European synchrotron radiation facility (ESRF) in Grenoble, France.This experimental station is dedicated to high-resolution X-ray imaging and features a highly coherent X-ray beam, a precision mechanics sample stage and a high-resolution detector system.During a 3D tomographic scan, the sample is rotated over 180 in order to acquire 1500 two-dimensional (2D) projection radiographs.From this set of projections, a quantitative 3D map of the attenuation coefficient distribution within the sample is produced by means of a standard filtered backprojection tomographic reconstruction algorithm.The spatial resolution obtained in the reconstructed images is of the order of 1 µm, a value comparable to the resolution of an optical microscope.One can also notice that the 3D scans are done with the maximal load: it ensures that the crack is open.However, even with a locally closed crack, diffraction effect allows to measure a crack opening until 100 nm.The sample was imaged in a dedicated fatigue machine designed at INSA Lyon to perform in situ fatigue test at the ESRF.The mechanical design of the cyclic tension loading mechanism enables operation at cycling frequencies of up to 80 Hz, minimising thereby the cycling time required for a fatigue test.
The geometry of the sample used is represented in Figure 1.A thin (2 µm) rectangular notch, 100 µmw i d ean d20µm deep, was produced in the sample using focused ion beam machining.This notch is located at the centre of one of the specimen faces and acts as a crack initiation site.
The in situ fatigue test was performed in air, at constant stress amplitude with σ max = 220 MPa, a stress ratio R =0 .1, a frequency of 40 Hz and at room temperature.In the fine-grained alloy studied here, the crack size (which is in the range 100-500 µm) is at least 100 times the grain size, and can hence be considered as microstructurally long.

Implicit Representation of the Crack Using the Level Set Method
From the previous experiments, 3D images are obtained and exploited in order to proceed to a simulation of the crack propagation.In this respect, one represents the crack by the use of two level sets.
Indeed, an implicit representation of the crack using the level sets method is very well suited to the eXtended Finite Element Method [1,2].
The crack front is defined as the intersection between a surface C crack that defines the crack surface and a surface C front that defines the crack front.These two surfaces, defined for a semi-elliptical crack, are displayed in Figure 4.I n the general case, these surfaces can have arbitrary 3D shape, and are defined by means of level sets functions.In order to couple the level set representation of the crack with the 3D images which come from X-ray microtomography, a specific strategy is proposed and is summarized in Figures 2 and 3.
The first step consists in a segmentation of the 3D image in order to extract the crack.The result is a set of voxels which defines a volume linked to the crack.The second step consists in the definition of a Boolean field calculated from the segmentation: a voxel outside the crack is initialized to 0 and a voxel inside the crack is initialized to 1.The main difficulty is now to be able to define a surface from this crack volume.In this respect, for the third step, it is needed first to "model" the crack as a surface with a predefined length scale.As a consequence, a length parameter is defined, and an intermediate triangulation step is done according to this length scale (in practice greater than the thickness of the crack).
Indeed, it seems essential in practice to model the crack as an explicit non-planar surface before to define the level sets.From the set of triangles, the 1D elements of the front are extracted with the local basis along the crack front.The fourth and the fifth steps consist then in defining the crack and the front level sets from the set of triangles and the set of 1D elements with their respective local basis.
In this respect, reinitialization equations (to the sign distance function) and orthogonalization equations proposed in [2] are used according to Figure 3 in order to compute the two initial 3D scalar fields linked to the crack and the front (see Figure 4).In particular, this approach can be done in a narrow band close to the crack.

Strategy of Enrichment and SIFS Calculation
In FE methods, the presence of cracks in a structure must be taken into account in mesh generation: the mesh must adequately define the crack geometry.Special elements and considerable mesh refinement near the crack tip are necessary to capture accurately the asymptotic displacement fields.The X-FEM alleviates the shortcomings associated with meshing the crack surface by representing the crack geometry using additional functions called enrichment functions.The enrichment method, as described in [1], can be summarised as follows.The displacement field calculated in the structure "without crack" is locally modified by adding, at specific nodes, the nodal values of the enrichment functions.These functions define the crack geometry by "modelling" the discontinuity introduced by the crack in the displacement field.As the discontinuities at the crack front and along the crack surface are different, two enrichment functions are necessary to model the entire crack [1].
The function H, defined as a generalised Heaviside function, models the displacement field along the crack surface.It is used to enrich the nodes of the element that are cut by the crack surface.
The enrichment F J is used to enrich the nodes of the elements that contain the crack front.F J consists of a span of functions which incorporate the  radial and angular behaviour of the asymptotic crack-tip displacement fields, where the local coordinate system can be extracted from the level sets [2].Furthermore, the nodal values of the crack level set and the front level set give the precise location of the crack and control whether a node has to be enriched or not (see Figure 5).
Concerning the calculation of the local stress intensity factor values (K I , K II and K III ) along the crack front, the interaction integral is used on a Jdomain defined in the local basis of the crack front [1,2].The J-domain is completely independent of the mesh of the structure and its size depends on the curvature of the crack front and the size of the finite elements cut by the front.On the boundary, due to plane stress effect, the interaction integral is not used and stress intensity factors in the bulk are extrapolated from the bulk to the boundary [4].In Figure 6, we illustrate the non-deformed mesh and the deformed mesh of the specimen submitted to a fatigue load after 14000 cycles.

Development of a 3D Local Paris Law
After 14,000 fatigue cycles, a crack initiated at the notch was detected on a 2D radiograph.The in situ fatigue test was performed at constant stress amplitude with σ max = 220 MPa, a stress ratio R =0.1, and a frequency of 40 Hz.A tomographic scan of the part of the specimen containing the crack was acquired and corresponds to the "first scan".During crack propagation, nine scans were recorded, reconstructed and segmented to obtain 3D renditions of the crack at different stages of its evolution.Projections in the (X-Y )planeof seven 3D renditions (chosen among the nine) are displayed in Figure 7.T h e dotted lines on the 3D crack renditions represent the location of the crack front at the previous step.The position along the crack front is defined by the angle x as shown in Figure 7(a).Qualitative and quantitative information on the evolution of the 3D crack geometry can be obtained.
The crack size at the surface (2c), for ω =0 • , and in the bulk (a), for ω =9 0 • , are measured on the 3D renditions of the crack.The fatigue crack growth rates at the surface, dc/dN , and in the bulk, da/dN , as a function of the stress intensity factor range ∆K =(K max − K min )aresh o w ni nFigure 9.The values of K max and K min are calculated for ω =0 • and for ω =90 • using domain integral techniques [1].
Thus for the same ∆K, the crack propagates faster in the bulk than at the surface implying that, at least for the specimen geometry used here, the crack growth behaviour is anisotropic.Thus, a fatigue crack growth law determined at the surface does not account for the crack behaviour in the bulk.Furthermore, if used for fatigue life calculation, the surface experimental law would lead to a non-conservative prediction.
One possible reason for the observed crack growth anisotropy between the surface and the bulk is variation of the closure stress along the crack front: the crack will then propagate faster in the bulk because the closure stress in this region of the crack front is smaller and hence the effective driving force is higher than at the surface.
Elber [6] proposed a modified Paris equation to account for the effect of closure on the crack growth rate: da/dN = C(∆K eff ) m where ∆K eff =  As a first approach, the variation of K op (ω) is taken to be linear with the angle ω and is therefore given by K op (ω)=aω + b.The value of the constants a and b are determined from the values of K op at the surface K op (0 • )a n da t the deepest point K op (90 • ).Data from the literature are used to determine K op (0 • ): the closure response of an Al-Li powder metallurgy alloy, very similar to the alloy studied here, was investigated in [7].The value of K op measured at the surface of CT specimens, for R = K max /K min =0 .1, is found to be equal to 0.4K max .We assume here that the ratio K op (ω =0 • )/K max remains constant during the fatigue test and that there is no closure at the deepest point for ω =9 0 • ,w h i c hg i v e sK op (90 • )=0 .1Kmax (90 • ).The evolution of the closure stress along the crack front is given by the linear relation: The 3D crack growth law taking into account the variation of the closure stress along the crack front is given by with C =1 0 −9.2 mm.cycle −1 .(MPa.m 0.5 ) −m and m =3 .51 for the previous aluminium alloy.In Equations ( 3) and ( 4), the angle ω, expressed in degrees, is defined for 0 • <ω<90 • and K max (ω) is obtained from a polynomial interpolation of the values calculated by the X-FEM as described in the previous section.This means that, for the sample geometry investigated here, a single Paris law can be used to predict the observed crack growth anisotropy provided the variation of the closure stress along the crack front is taken into account.
From a more general point of view, this study has shown that coupling X-ray microtomography to X-FEM provides a promising tool to assess the 3D behaviour of arbitrary shaped cracks and to perform direct comparisons of "experimental" and "simulated" crack shapes during propagation.Such data, to the best of the authors' knowledge, are currently lacking in the literature.
In this study, closure effect is taken into account empirically in the Paris law.An extension of this work consists in using explicitly a elastic-plastic behaviour with contact and friction in the simulation.This is the aim of the next section.

Strategy of Enrichment
The main purpose of this paragraph is to treat the case of bulk and interface non-linearities in the framework of the eXtended Finite Element Method.The presented method will focus on the case of plasticity combined with frictional contact and is applied to fatigue crack growth analysis.The aim of the method is to use the X-FEM to model the well known fatigue crack closure under cyclic tension phenomenon first observed by Elber [8].
In the case of material non-linearities, several issues have to be addressed.In a previous paper [5] an appropriate elastic-plastic enrichment basis was developed to allow the X-FEM to deal with plasticity.In the present paper this approach is coupled with the treatment of frictional contact, and the problems associated with plastic crack propagation are explored.The aim of the study is to develop a general strategy for a wide family of elastic-plastic behaviour (kinematics and/or isotropic hardening, etc.).In this approach, neither remeshing nor thermodynamic field interpolation during propagation of the crack are needed.Furthermore, the crack front is enriched with specific basis in order to use sufficiently coarse mesh compared to FEM.In fact, the aim is to be able to take into account closure effect and confined plasticity with a good accuracy during propagation of the crack: indeed, this can have a great influence on the fatigue crack growth.In order to accurately model the confined plasticity, a geometrical subdividing zone is defined with an estimation of the plastic zone which can be smaller than the spatial discretization assumed to be fixed (see Figure 11).In the same way, during the propagation of the crack, one defines new sub-elements between the two configurations according to the estimation of the new plastic zone (see Figure 12): only the Gauss points of the new sub-elements are needed (which allows a good accuracy of the plastic behaviour during the crack growth with no remeshing and no interpolation field) [3].Concerning the strategy of enrichment of the element containing the crack tip, an appropriate enrichment basis for fatigue crack growth simulation with confined plasticity is used [5]: Fig. 12. Evolution of the sub-elements between the two configurations submitted to loading and unloading even with coarse mesh for a wide range of elastic-plastic hardening [5].
In order to take into account contact and friction between the crack faces, interface elements are implemented along the crack with their own displacem e n ta n dt r a c t i o nfi e l d s( t and w).In this respect, a three fields Augmented Lagrangian formulation is introduced which allows to couple the interface elements with the displacement and the Cauchy stress fields (u and σ) [ 3]: where σ and ε are respectively the Cauchy stress and the strain fields; t and w respectively the load and displacement fields along the crack faces; F d the prescribed load; u * , w * and Λ * respectively the virtual displacement field in the bulk, the virtual displacement field on the crack faces and the virtual Lagrange multiplier field; α the penalty term (for the normal and the tangential problem) and Λ the Lagrange multiplier of the Augmented Lagrangian formulation.Finally, concerning the update of the local enrichment during propagation the following strategy is retained: in the aim to preserve the history of the thermodynamics variables, all enrichments are retained during crack growth [15].As a consequence, the change of discretization is performed by imposing the new crack segment to be closed (new enrichments are initialized to zero; the new t and w must be compatible with u and σ).Furthermore, the new stress and internal variables are initialized with elastic conditions thanks to the moving subdividing zone (see Figure 13).One considers as a first example a Compact Tension Shear (CTS) specimen, submitted to mode I fatigue crack growth with interspersed overload.A numerical and experimental investigation of that specimen is fully presented in [9].In this example, one wants to compare the results obtained by a classical FE calculation presented in [9] with the proposed method.In Figure 14,t h e geometry and boundary conditions are presented, and one also illustrates the mesh of the specimen with a zoom around the crack tip for the FE simulations and the X-FEM simulations.
The material is chosen similar to the one in [9]: 2.1011 Pa for the Young's modulus, 0.3 for the Poisson's ratio, 200 MPa for the yield strength, 534 MPa for the hardening coefficient, 0.27 for the hardening exponent n,and 90 MPa.m 1/2 for the critical mode I SIF.In this comparison, the X-FEM mesh as been chosen 10 times coarser than the FEM mesh.Furthermore, a standard von Mises plasticity with isotropic hardening (Newton+radial return scheme) has been considered here [3].The specimen is submitted to a total crack growth of ∆a =0 .5mm( 2 0s t e p si n FEM, 2 steps in X-FEM), then submitted to an overload with a ratio of 2.5 and then a total growth of delta = 0.5 mm.One can conclude on this example that, even with a 10 time coarser element length and a 10 time coarser time discretization, very good agreements can be obtained between FEM and X-FEM on the vertical displacement close to the crack tip before and after the propagation of the crack.One can also notice the closure of the crack due to the overload and the confined plasticity (after the propagation).
As a second example, one considers a mode I Compact Tension (CT) specimen with a loading ratio of 0.1.The geometry and dimensions of the specimen are given in Figure 17.The material is chosen to be similar to the previous example.The specimen is submitted to a tension cyclic loading in order to have a stabilized stress state around the tip, then a growth of ∆a =0.05 mm is imposed at maximum load, and so on.In Figure 18, one illustrates the ability of the method to simulate closure effect close to the crack tip.In Figure 19, one illustrates stress-strain curves obtained on a Gauss point close to the crack tip during the propagation and the cycles of stabilization.One can notice that, due to the fact that no remeshing is done during the propagation of the crack, it is very easy to follow the stress-strain evolution of the material in a fixed point, which can not be possible with the same accuracy if you proceed to remeshing and interpolations.
Concerning the time dependence of the problem, this case can be considered as a two scale approach: a fixed "coarse" scale linked to the geometry (structure and crack) and a "fine" scale linked to the non-linear properties of the material (confined plasticity, contact and friction between the crack faces) which follows the crack tip during the propagation.In this respect, no remeshing is needed; however, new integration points are needed (it consists to initialize the new internal variables to zero according to the linear elastic behaviour outside the crack tip plastic zone), and new interface elements are needed.
Numerical simulations show that the proposed basis of enrichment can give accurate results for specimens submitted to loading and unloading even with coarse meshes for a wide range of hardening.However, because of conditioning issues, the crack extension can not be small compared to the element size (discrete successive plastic zone).For cyclic loadings with elastic-plastic crack growth, the question of the optimal enrichment basis is still open.One can also notice that the effects of non-singular terms of the asymptotic fields like T-stress are not taken into account.Furthermore, numerical issues suck as locking, plastic incompressibility and convergence rate need specific studies.

Conclusions
The propagation of a semi-elliptical crack in the bulk of an ultrafine-grained Al-Li alloy has been investigated using synchrotron radiation X-ray microtomography.In this material, the studied crack, despite its small dimension, can be considered as a "microstructurally long" and described in the frame of the linear elastic fracture mechanics.The extended finite element method was used to calculate the stress intensity factors along the crack front taking into account the three-dimensional geometry extracted from the tomographic images.For the same nominal value of the stress intensity factor range, crack propagation was faster in the bulk than at the surface.The observed anisotropy is attributed to the variation of the closure stress along the crack front between surface and bulk.The experimentally observed fatigue crack propagation is compared to numerical simulations.Good agreement is found when a linear variation of closure stress along the crack front is taken into account in the "3D crack propagation law" used for the simulation.In a second example, it is shown that even with elastic-plastic behaviour coupled with contact and friction, non-remeshing property can be preserved for instance for 2D mixedmode plastic fatigue crack growth.The main idea consists in using a specific nonlinear enrichment basis which allows to take into account the asymptotic behaviour around the crack tip.Because of the treatment of multiple nonlinearities (plasticity, contact), the numerical integration is adapted in order to properly integrate the enrichment basis and to have a fine knowledge of the stress state around the crack tip.Second, the propagation strategy consists in keeping the history of the variables such as stresses and plastic strains.Finally, an original strategy was designed in order to avoid projection of stresses and internal variables as the crack evolves to ensure the reliability of the method.In the presented examples, the method shows its ability to model the physical phenomena that are present in fatigue crack growth.

Fig. 1 .
Fig. 1.Geometry of the sample used for the in situ fatigue test monitored by X-ray microtomography

Fig. 3 .
Fig. 3. Initialization of crack and front level sets

Fig. 4 .
Fig. 4. Representation of the 3D scalar fields linked to the crack and front level sets

Fig. 5 .
Fig. 5. Set of elements completely cut by the crack; set of elements cut by the front; 2D representation of the crack extracted from the crack level set

Fig. 6 .
Fig. 6.Representation of the non-deformed and deformed mesh of the specimen submitted to fatigue loading

Fig. 11 .
Fig. 11.Geometrical subdividing zone with an estimation of the plastic zone

Fig. 15 .
Fig. 15.Mesh of the CTS specimen with zoom around the crack tip

Fig. 17 .
Fig. 17.Geometry and X-FEM mesh for the mode I CT specimen