Development of a finite element enriched method adapted for the case of multiple cracked structure

The present method is based on the eXtended Finite Element Method and the partition of unity properties which allows the addition of enriched functions in the displacement field discretization. Those correspond to the analytical displacement field near a straight crack. The additionnal degrees of freedom correspond directly to the Stress Intensity Factors. This method is validated in mode 1 and mixed mode. A parametric study has been performed to analyse the effect of different parameters on the accuracy of the method. The case of two shading cracks has also been investigated.


Introduction
The eXtended Finite Element Method (X-FEM) allows crack modeling independently of the mesh ( (Moes et al., 1999) and ). Thanks to the partition of unity properties of finite element shape functions (Melenk et al., 1996), enrichment functions can be added in the displacement formulation. Thus discontinuities and crack tip singularities can be captured more accurately. Heavyside enrichment function has been used to represent the discontinuity (HX-FEM (Moes et al., 1999)). The crack tip singularity is given by the basis of the analytical displacement field of Westergaard near the crack tip. This formulation leads to additional degrees of freedom. The enrichment is applied at the nodes of the elements which are cut by the crack. To obtain a more accurate estimation of the stress intensity factors a wider enriched zone should be used around the crack tip (Laborde et al., 2005). In the case of a structure with many cracks, the numbers of additional DOF can be high. In this kind of problem, the calculation of the stress intensity factors is very interesting to point out the critical cracks in the structure. The calculation of the stress intensity factors is usually performed using the interaction integral method (Yau et al., 1980) or the least square identification (McNeill et al., 1987). The use of these methods could be very difficult in the case of crack interacting closely or even intersecting. Specific discontinuous basis functions have been used to make calculation of intersecting cracks ,  and (Daux et al., 2000). However if the number of crack is high and if the crack tips become too closer the calculation could become very complicated. Analytical studies had been made for interacting crack (Kachanov, 1987) and (Li et al., 2003) or intersecting cracks (Cheung et al., 1992) using superposition principle to estimate stress intensity factors or equivalent elastic properties of a multiple cracked structure (Helsing, 2000). These examples are made for semi-infinite crack in infinite medium. To combine the advantages of the X-FEM in making calculation of finite cracked structure and those of the analytical method in the accuracy of stress intensity factors evaluation, the enrichment functions can be taken as the analytical displacement fields (Liu et al., 2004) such as the formalism of Williams' series (Williams, 1957). The method presented in this paper is based on a specific analytical enrichment for small straight cracks which provides directly the stress intensity factors with only two additional degrees of freedom per crack.
In the next section, the displacement formulation is detailed. In Section 3, the capabilities of the method are tested in mode 1 and mixed mode. The Section 4 discusses about the influence of different calculation parameters on the accuracy. Finally in Section 5, the interaction of two shading cracks is investigated.

Enrichment formulation
The enrichment developed herein is based on the analytical displacement field of Griffith-Inglis (Weertman, 1996). In this model cracks are accounted for such as a continuous sum of dislocations. The discontinuous displacement field (U d(x, y)) in mode 1 and 2 linearly depends on the stress intensity factors at the crack tips (K i Equation [1]). In this formulation cracks are symmetric. Thus the stress intensity factors are the same at both tips of the crack and then only one SIF is considered for each mode.
The total displacement field is decomposed into a discontinuous and a continuous part [2].
The continuous part is described using standard finite elements (U c(x, y)) and the discontinuous one (U d(x, y)) with the analytical enrichment functions (F j (x, y)) by using the partition of the unity properties (Melenk et al., 1996) N i (x, y) are the standard finite element shape functions. N is the total number of nodes in the structure. A domain around the crack tip is chosen and N H defines the number of nodes inside this domain. Two additional degrees of freedom are considered (K i 1 and K i 2 ) at each enriched node i. The following condition ∀i ∈ N H , K i j = K j is enforced to decrease the number of enriched degrees of freedom. Eventually there are only two additional degrees of freedom per crack. The corresponding eXtended Finite Element approximation is given by Equation [4]. Φ(x, y) is a cut-off function (Chahine et al., 2008). It is equal to 1 in all elements containing N H nodes and decreases linearly to 0 in adjacent elements. The displacement field is finally written as, The formulation leads to the calculation of the global elastic stiffness matrix K (Equation [5]). It is shared into a 2N × 2N square matrix which is the classical finite element stiffness matrix K c , a 2 × 2 square matrix which components represent the enriched terms (K e ) and a 2N × 2 which is the coupling matrix (C) (Equation [5]). The integration of enriched functions has ever been studied by several authors. The first developed idea was to introduce integration subcells into enriched elements. The construction of the subcells depends on the position of the crack into the element and it is often complicated to elaborate. Integration of discontinuous function using the contour integral transformation as proposed in (Ventura, 2006) and in (Ventura et al., 2009) has been used to circumvent this drawback. In our case the continuous components of the matrix K c are integrated with four gauss points. The enriched function are not given in a simple way and the integration cannot be evaluated exactly with usual quadrature rules. To integrate the coupled and the enriched components (K e and C) the elements are subdivided into a given number of square subcells in the reference frame. Inside each of these subcells, the integration is performed with four gauss points (Figure 1 and Equation [6]). This integration is an approximation but it gives accurate results with at least 10 × 10 subcells (n cell ). In this formulation the additional degrees of freedom correspond directly to stress intensity factors K j which are directly calculated during the resolution of the linear system without any postprocessing.

Study of the analytical extended finite element method in mode 1
The validity of the method is verified with a simple traction test on a sample containing a crack inclined by the angle β (Figure 2). In mode 1, β is equal to 0 and the relative crack length λ = a W increases from 0.02 to 0.25. a and W are respectively, the half of the crack length, and the width of the structure. The finite size correction on K 1 depends upon λ, it is given for example in (Lemaitre J., 1985) (Equation [7]). In mixed mode, stress intensity factors depends on β. They are given by Equation [8]. In this case, the crack length is fixed, λ = 0.1 and β varies between 0 and π. The material is chosen with a Young's Modulus of E = 210GP a and a Poisson's ratio of ν = 0.3. A unit applied stress is used σ A = 1M P a. A uniformed mesh of 45 × 45 quadrangles is used. Cut-off function support dimensions are defined by a radius equal to r = 0.4W . The evolution of the stress intensity factors with λ is compared with those given by Equation [7] in Figure 3. A good agreement is observed between the two curves. However a bigger difference is observed when λ is high. This is due to the fact that the crack length becomes similar to the cut-off function diameter. The accuracy depends on the relative size of the crack, the structure and the cut-off function radius. Figure 4 compares the stress intensity factors evaluated by the proposed method and by Equation [8]. No fluctuation is observed in mixed mode. K 1 is closed to 1 and K 2 closed to 0 when β = 0. With β = π 4 , K 2 is maximal and equal to K 1 . This study shows the validity of the method. However dimensional parameters seem to influence the accuracy of the stress intensity factors calculation.

Parametric study of the enrichment
All the studied parameters are listed in the table 1. The convergence, as a function of the number of element along one direction h/W, where h is the mesh width and W is the width of the structure, is studied with a fixed crack length and a fixed cut-off function size (λ = a/w = 0.1 and γ = r/w = 0.2). The mesh is uniform and it consists in quadrangle elements. The convergence is studied with respect to the analytical displacement applied on the boundary nodes of the structure. This prescribed displacement is given by the same function as the one we use as an enrichment. The boundary conditions thus depends on prescribed stress intensity factors. A relative error is calculated between the prescribed stress intensity factors and the computed ones. The relative error (ǫ r ) is plotted in Figure 5. The error is less than 10% even with a coarse mesh of 10 × 10 elements. The convergence rate is about 0.62 which corresponds to the order of magnitude for usual X-FEM approach with direct extraction of stress intensity factors (Nicaise et al., 2009).  The study in traction mode emphasises boundary effects as the crack length and the cut-off function radius become similar to the structure width. The accuracy of the calculation is studied with the simultaneous variation of relative crack length (λ) and the relative cut-off function radius (γ). The values of interest are given in Table 1. When γ increases and λ decreases, The relative error (ε r ) on the evaluation of the stress intensity factors decreases (Figure 6). With a fixed value of λ, the precision increases with the size of the enriched zone. The cut-off function has to cover a sufficiently large area around the crack to ensure a correct accuracy. The ratio γ/λ has to be controlled in order to ensure a good precision. In Figure 6 a plane corresponding to ε r = 0.05 indicates the values of λ and γ to provide a relative error lower than 5%.

Multiple cracks
A multiple cracked structure is considered. The material is supposed to be linear elastic and isotropic. The displacement field falls into the continuous part given by the standard finite element basis and the discontinuous part given by the analytical enrichment functions. Each of the discontinuous displacement is given by the corresponding enrichment function F k j . j corresponds to the mode and k is the number of the crack. Thus, the displacement field is searched based on the following decomposition : The number of degrees of freedom is given by the standard finite elements part and two additional degrees of freedom per crack. Those correspond directly to the stress intensity factor in mode j (1 or 2) for the crack k (K kj ) as in the previous sections of the paper. Each crack has a cut-off function Φ k (x, y). These functions can be separated or not. Two shading interacting cracks are considered (Figure 7). δ = ∆/W is the parameter which controls the relative distance between the two cracks. Figure 7 also illustrates how two intersecting cut-off functions are built. Figure 8 gives the evolution of K 1 with δ. The parameter λ is chosen equal to 0.1 thus the crack length is five time smaller than the structure width. To evaluate the interaction between the cracks, two calculations have been performed. The first one considers one crack with a distance of ∆ from the median line. A second one is carried out with two cracks such as described in Figure 7. In the first case K 1 is not constant with δ due to the boundary effects observed in section 3. When the crack is far from the boundaries K 1 remain almost constant. The calculation with two cracks shows that the second crack has a strong effect on the values of stress intensity factors. The closer the cracks, the lower the stress intensity factors are. Moreover the values tend to 0.5 which means that the two cracks tend a lonely crack with a K 1 equal to 1. However when δ is high the two curves are close because the interaction between the two cracks is small. This is confirmed by the amplified deformed shapes (pictures 9 (a) and (b)). When δ is small the two cracks are very close and it seems that the crack opening is lower. This observation is in agreement with the evolution of the stress intensity factor with δ.
The cut-off functions defines the numerically influenced area of each crack. Despite the fact that the two cut-off functions intersect or not, the evolution of the stress intensity factor remains continuous. This observation shows that the cut-off function enables a good mechanical connection between the continuous and the discontinuous parts.

Conclusion and discussion
This work shows that the X-FEM method is also a useful tool for multiple crack problems when analytical enrichment functions are adopted. The analytical enrichment used in the present paper corresponds to the displacement around a crack of finite size under symmetric loading (the values of the stress intensity factors are the same at both tips of the crack). The procedure has first been validated in a structure with one crack. Results show a good accuracy in the estimation of stress intensity factors. However calculation parameters such as crack length, cut-off function width, have to be controlled in order to ensure the quality of the numerical results. Cracks have to stay relatively small compared to the structure length. However, the convergence rate in term of stress intensity factors is the same as for other X-FEM based methods with direct estimation of stress intensity factors. It is shown that crack interactions can be efficiently investigated using the proposed method. As only two additional degrees of freedom are required per crack, one can imagine to study structures with a large number of cracks for a very limited computational cost.