Robust and direct evaluation of J2 in linear elastic fracture mechanics with the X‐FEM

The aim of the present paper is to study the accuracy and the robustness of the evaluation of Jk‐integrals in linear elastic fracture mechanics using the extended finite element method (X‐FEM) approach. X‐FEM is a numerical method based on the partition of unity framework that allows the representation of discontinuity surfaces such as cracks, material inclusions or holes without meshing them explicitly. The main focus in this contribution is to compare various approaches for the numerical evaluation of the J2‐integral. These approaches have been proposed in the context of both classical and enriched finite elements. However, their convergence and the robustness have not yet been studied, which are the goals of this contribution. It is shown that the approaches that were used previously within the enriched finite element context do not converge numerically and that this convergence can be recovered with an improved strategy that is proposed in this paper. Copyright © 2008 John Wiley & Sons, Ltd.


INTRODUCTION
The numerical simulation of fracture mechanics is of major importance in the design of structures such as aircrafts or dams. For complex structures, robustness and accuracy of predictions are essential pre-requisites to use numerical methods for new parts design. The finite element method (FEM) is classically considered for the simulation of such problems. However, cracks must be explicitly meshed in order to incorporate the discontinuity of the displacement field across their faces. In addition, this process has to be handled at each step of propagation of the cracks. The main drawback of this approach is the computational and human cost of remeshing. Nevertheless, numerical tools improved much over the past years. Meshless methods (such as EFGM [1] or hp-clouds [2]) were proposed to avoid meshing of the domain studied and to improve the quality of the results by enriching the approximation. FEMs were also thoroughly improved using the concept of partition of unity [3], which was first employed in the context of meshless methods by Oden and Duarte [2,4]. Among the class of partition of unity FEMs, the generalized finite element method (GFEM) and the extended finite element method (X-FEM) are the most advanced. The GFEM was introduced by Strouboulis et al. [5,6] and was applied to the simulation of problems with complex micro-structures (see also the work of Babuška et al. [7] and Oden et al. [8]). The method was further extended to employ the idea of mesh-based numerically constructed handbook functions by Strouboulis et al. [9,10].
Concerning the X-FEM, it was first proposed as a solution to the remeshing issue for crack propagation in linear elastic fracture mechanics [11,12]. It uses the partition of unity in two ways: first to take into account the displacement jump across the crack faces far away from the crack tip and second to enrich the approximation with analytical asymptotic fields close to the tip. Many fracture mechanics problems have been successfully solved with the X-FEM approach: 2D [12][13][14][15][16], 3D [17][18][19][20], plates and shells [21][22][23], cohesive zone modelling [24,25], dynamic fracture [26], non-linear fracture mechanics [27][28][29][30][31] among others. For a wider review on partition of unity finite elements applied to fracture mechanics, the reader could also refer to [32].
A major issue in fracture mechanics is the prediction of crack kinking. In order to predict crack growth direction, several criteria have been proposed: one can cite, for example, the maximum hoop stress criterion [33], the maximum tangential stress criterion or the symmetry principle [34]. The maximum energy release rate criterion [35], also called as the 'vectorial J -integral' criterion [36], is an alternative, which was initially proposed by Hellen and Blackburn [36]. The main difference in evaluating this criterion with respect to the classical J -integral relies on the difficulties in the evaluation of the second component J 2 of the energy release rate vector J. It is well known that J 1 (which is basically the J -integral) is path-independent. Concerning J 2 , Herrmann and Herrmann [37] have shown that it is path-independent only in a modified manner, i.e. by adding an extra contribution on the crack faces. However, this new contribution requires an accurate determination of the mechanical fields near the crack tip. It is well known that the quality of these fields near a singularity is not ensured in the finite element context. This is why some authors proposed specific strategies for an accurate numerical evaluation of J 2 . The pioneering contribution is by Eischen [38] who proposed to decompose the integration into two parts: The first one far from the tip, evaluated using finite elements, and the second one close from the crack tip, calculated using near-tip asymptotic fields. Evaluating two contour integrals, it is then possible to calculate J 2 . One drawback of this approach is that it involves contour integrals, which are not well-adapted to finite elements. This contrasts with the evaluation of the J -integral that can be performed considering a domain integral by means of a virtual velocity field [39]. This is why Chang and Yeh [40] proposed an approach similar to that in [38], but with the use of domain integrals.
For enriched FEMs, the only study on this subject was recently proposed by Heintz [41] in the context of configurational forces. His approach is directly related to the one proposed in [40], except that no decomposition of the integral is considered to avoid the near-tip integration. The major advantages of this direct approach are that it does not necessitate near-tip fields assumptions and that it requires only one computation of the domain integral.
The accurate prediction of the crack kinking angle highly depends on the accuracy and robustness of the J-integral calculation. Nevertheless, no study has been proposed to evaluate the accuracy of the above-mentioned methods. This contrasts with the numerical evaluation of stress intensity factors within X-FEM, which has already been investigated in 2D and 3D [19,20,23,42].
The aim of the present paper is to compare various numerical approaches for the evaluation of the J-integral. It will be demonstrated that some of these approaches do not converge. Improvements will be proposed in order to overcome the lack of convergence. This paper is organized as follow. First, the governing equations are presented. Then, X-FEM basics are recalled. The third section presents the classical approaches used to compute the J-integral; finally, a new method is proposed in order to ensure convergence of the calculation.

GOVERNING EQUATIONS
Consider the static response of a 2D cracked elastic body that occupies a bounded domain with a sufficiently smooth boundary * split into two disjoint parts: * u where displacements are prescribed (Dirichlet boundary conditions) and * t where tractions are prescribed (Neumann boundary conditions) (see Figure 1); * t includes the crack C . The body is initially in an undeformed and stress-free configuration. The governing equations are as follows: where is the Cauchy stress tensor, b is the body force, u d is the prescribed displacement field, T d are the prescribed tractions, n is the outward unit normal to the boundary * , is the linearized strain tensor and C is the fourth-order elasticity tensor. For a linear isotropic elastic material, the constitutive equation depends only on two scalar parameters and (the Lamé coefficients) and can be expressed as The above strong form of the governing equation (1) can be expressed into the following weak form:

THE X-FEM
The above weak formulation (3) will be numerically solved considering the approximation of displacement and test fields. With classical finite elements, the approximation of a vector field u(x) on an element e is given as where n is the number of coefficients describing the approximation over the element, u is the th nodal value of this approximation and N a is the classical vectorial shape function associated with the degree of freedom (dof) u . Within the partition of unity, the approximation can be enriched as where n e is the number of enrichment modes, a is the additional dof associated with dof and stands for the th scalar enrichment function. Two different enrichments are considered in fracture mechanics [11,12]: • The first enrichment is a discontinuous enrichment for nodes whose support is bisected by the crack. In this case, the interpolation of the displacement field is discontinuous across the crack faces. A Heaviside jump function is considered to model this discontinuity: it equals +1 on one side of the crack and −1 on the other side. The associated dof manages the magnitude of the displacement discontinuity. The Heaviside function is computed using the level set representation of the crack [17,43]. • The second enrichment is a near-tip enrichment for nodes whose support contains the crack tip. In linear elastic fracture mechanics, these enrichment functions are determined using the asymptotic displacement field near the tip of the crack.
Finally, it is important to note that the modified Gauss quadrature scheme described in [12] is used to integrate both discontinuous and non-polynomial functions over the elements. For a more detailed presentation of fracture mechanic applications within X-FEM, the reader can refer to [12,17,18,[43][44][45].

NUMERICAL EVALUATION OF THE VECTORIAL J -INTEGRAL
In this section, the definition of the J-integral is first recalled. Then, strategies are presented in order to numerically evaluate this quantity. Next, attention is paid to the accuracy of these approaches in the context of enriched FEMs.

J and (J k ) k=1,2 integrals
Consider a 2D cracked body and the crack local coordinate system (e 1 , e 2 ) as depicted in Figure 2. Consider a contour that encloses the crack tip and n denotes its unit outward normal. The J -integral is defined as where n 1 is the projection of n on e 1 and W is the strain energy density of the material. This integral is a contour integral that was seen to be path-independent [46]. A second integral called J 2 can be introduced such that both J k (k = 1, 2) integrals are defined as where M = W I −∇u T is the Eshelby stress tensor [47]. Using this expression, the J -integral is the first component of a vector J whose components on (e 1 , e 2 ) are J 1 and J 2 , respectively. J 2 is not path-independent. Nevertheless, it is said to be path-independent in a modified sense [38]: If does not tend to the tip anymore, an extra boundary term has to be taken into account on the crack faces, and J can be expressed as where 'W ( stands for the jump in the strain energy density across the crack faces. Note that this extra term influences only J 2 because n·e 1 = 0 on C (for a straight crack ‡ ). ‡ Otherwise, e 1 has to be defined as a variable along the crack faces [17].

Numerical evaluation of J
As mentioned above, the extra boundary term in Equation (10) has not to be evaluated for J 1 . However, for J 2 the integration of the jump in the strain energy across the crack faces should be considered: this contribution does not vanish in general. In fact, the stress field near the crack tip can be expressed ass where x0 is classically referred to as the T-stress [48]. Eischen [38] showed that, in this case, the jump in the strain energy is given by The jump will be zero only for pure mode I loading conditions (K II = 0) or for x0 = 0. This also shows that the integration of the boundary term will be quite difficult because of its 1/ √ r behaviour. This is why various approaches have been proposed in order to accurately evaluate J 2 . The first proposal was given by Eischen and is developed in the following paragraph. [38] stated that the difficulties in the evaluation of J 2 is only due to the integration close to the crack tip. As the asymptotic behaviour of the mechanical fields is known in this region, he proposed to split the term into two integrals:

Accurate evaluation of J 2 : Eischen's proposal. Eischen
where represents the length of the crack faces where the singular behaviour dominates: it defines two subsets of C ( * C and C ) as depicted in Figure 3. The author proposed to evaluate the first term of the right-hand side of Equation (13) by using the finite element approximation and the second term by considering the asymptotic energy jump (12). Then, Equation (13) where the multiplier depends on the geometry of the problem, the loading conditions and the material moduli (i.e. is invariant for a given problem). Eischen proposed to evaluate J 1 and J 2 considering two values of , i.e. 1 and 2 . Then, with three computations (J 1 , J 2 ( 1 ) and J 2 ( 2 )), it is possible to evaluate the stress intensity factors, the T-Stress and finally J. This approach seems to lead to satisfactory results even for rather large values of . However, no convergence study was proposed in order to assess the accuracy of the method. Moreover, this approach is based on contour integrals that are not as accurate as domain integrals (see hereunder) for FEMs. A similar method with domain integrals was proposed by Chang and Yeh [40]; it is presented in the following section.

Accurate evaluation of J 2 : domain integral approach.
The domain integral method is a classical way to compute energy release rates within FEM and X-FEM. It was first introduced by Destuynder et al. [39] (see also Li et al. [49] and Shih et al. [50]). This method is very robust and has a convergence rate in O(h), even with classical finite elements. The method introduces a virtual velocity field Q = Qe 1 , e 1 being the unit vector tangent to the crack faces, and Q a sufficiently smooth scalar function, which equals 1 within , and 0 on the contour e of the domain (see Figure 4 for the notations). Then, the J -integral around can be expressed in an equivalent domain form as where V denotes the area between e and , which is equal to the area inside e as tends to zero. Following [40], this approach can be further extended to compute J rather than J . In this case, two virtual velocity fields Q 1 and Q 2 are introduced Q 1 = e 1 and Q 2 = e 2 with = 1 on and = 0 on e Vector Q 1 (resp. Q 2 ) is tangent (resp. normal) to the crack faces. The kth component (k = 1, 2) in the (e 1 , e 2 ) basis of the domain J-integral is In this expression, the term J B k (resp. J V k ) will be called a boundary term (resp. domain term) in the following. Note that the weight function is chosen so that it is equal to 1 inside and is equal to zero on e (plateau function). Finally, this approach can be easily extended to 3D problems [17]. The boundary term that appears in Equation (17) § is equivalent to that in Equation (10) and vanishes for the J 1 -integral but not for J 2 . The accuracy for the computation of the domain term J V 2 has been demonstrated both in the X-FEM and in the classical finite element context. Thus, the key issue for the calculation of the vectorial J -integral is the accuracy of the computation of the boundary term J B 2 . Chang and Yeh [40] proposed to proceed like Eischen, i.e. by splitting J B k into two parts. They considered the cases of both linear and non-linear elasticity. In the last case, a system of non-linear equations that results from the data set (geometry, loading conditions) with a number of different values of should be solved.
Here, we study the accuracy of the domain integral approach. However, we will not consider the previous works [38,40], because we would like to evaluate J without too many computations: two domain integrals, then the stress intensity factors and finally the J-integral should be computed. If one wants to compute the J-integral by means of the stress intensity factors, a simpler approach should consist in evaluating them directly using a domain interaction integral, which is known to be robust and accurate. Hence, we will directly evaluate Equation (15) similar to that by Heintz [41] who applied this approach in the context of enriched finite element approaches. In this work, the author derived a similar expression in the context of configurational mechanics. However, no convergence study was proposed so that the accuracy of the approach is not demonstrated. Moreover, no asymptotic enrichment was considered in the near-tip region, which means that the quality of the finite element fields may not be sufficient to evaluate J B k in a robust way (i.e. with accuracy, irrespective of the mesh topology).

Imprecisions with the domain integral approach
It has been shown that the evaluation of the domain J-integral was leading to the sum of a domain term J V and a boundary term J B (Equation (17)). The latter term involves the integration of the energy jump across the crack faces. However, the energy evolves in r −1 near the crack tip (even if the jump follows an r −1/2 evolution): if the r −1 terms do not vanish in the numerical evaluation of 'W (, the boundary term cannot be integrated. ¶ In this section, we first compare the errors due ¶ As the integral between and 1 of 1/r does not converge when goes to zero. to J V and J B using analytical fields and second exhibit the presence of spurious r −1 terms in J B when enriched finite element fields are considered. Consider a 2×2 square-structured mesh of 882 triangular elements as depicted in Figure 8(a). This mesh is fully enriched with the near-tip basis Equation (6). The displacement field corresponding to a cracked plate in a mixed mode (K I = K II = 1.0, E = 1.0 MPa, = 0 and x0 = 2.0 in Equation (11)) is then projected onto the enriched finite element displacement field by means of an L 2 projection: As the exact solution field is contained in the enriched finite element approximation, the resulting numerical field exactly represents the analytical field. The J-integral is then evaluated using a domain of radius 0.2. The errors due to the different terms of J are summarized in Table I. The two components of J can be accurately computed, especially J 1 . The presence of spurious r −1 terms is tested by changing the number of Gauss points for the integration of J B . If r −1 terms are still present in 'W (, then the value of J B 2 should degrade as the number of Gauss points increases. Figure 5 shows that J B 2 converges with a O(n) rate (n being the number of Gauss points). This demonstrates that the computation of J using domain integral is accurate when using mechanical fields of good quality. Next, the case where the mechanical field is obtained from a finite element solution is considered: the analytical stress field (Equation (11)) corresponding to the displacement Table I. Integration error when J is evaluated after L 2 projection of the analytical displacement field on the finite element displacement field (analytical value J 1 = J 2 = 2.0).  Figure 5. Convergence of J B 2 versus the number of Gauss points per element for a given mesh size (displacement field obtained using an L 2 projection).  field used in the L 2 projection is applied on the boundary of the mesh and the problem is solved. The radius for the near-tip enrichment is set as 0.4 and J is evaluated using the same domain radius (equal to 0.2) as in the projection. Table II shows that the degradation in the accuracy of J 1 is small, whereas it is huge for J 2 , and that this degradation is mainly due to the boundary term J B 2 . Finally, the convergence of J B 2 with respect to the number of Gauss points shows that J B 2 does not converge to its exact value, either because of the quality of the mechanical fields or because of the remaining spurious r −1 terms (see Figure 6). We can then conclude that the imprecision in the computation of J B 2 is a consequence of a global lack of quality in the finite element solution. However, it is still not clear whether convergence can be obtained using finite element fields.

Improvement in accuracy and robustness for the computation of J
An alternative approach is proposed for the evaluation of the J-integral at the crack tip. Equation (17) shows that the computation of J B 2 leads to the integration of the energy jump between the crack faces. As shown in Equation (12), this term vanishes if the (regular) T-stress x0 is zero. Moreover, if the T-stress can be computed, then the boundary term of the domain integral vanishes if the regular stress field x0 e 1 ⊗e 1 is subtracted from the singular stress field (see Equations (11)). This operation does not change the value of the J-integral [48], but allows one to remove the boundary term from Equation (17). The T-stress can be evaluated using a path-independent interaction integral I A+B [51]: where the superscript A denotes the finite element field and the superscript B denotes an auxiliary field as originally proposed by Michell [52]. This auxiliary field corresponds to the analytical solution to an infinite cracked plate loaded at its tip by a force f = f e 1 (see Figure 7). Sladek and Sladek [51] have shown that the interaction integral Equation (19) is related to the T-stress by where E * = E/(1− 2 ) and E * = E for plane strain and plane stress problems, respectively. After having determined the T-stress, the J-integral can be evaluated using only domain integrals and without considering the boundary term, which is the major source of error.

Convergence study
The convergence of the two domain approaches presented above is now studied by considering : with E * = E/(1− 2 ) and E * = E for plane strain and plane stress problems, respectively.
The T-stress x0 is set to as 2.0 (the boundary term is then non-zero in Equation (17)), and the convergence study is performed using a domain of radius 0.2 for the computation of domain integrals. Four cases are considered: 1. domain integral method (radius 0.2) using X-FEM without near-tip enrichment (the crack tip ends on the edge of an element); 2. domain integral method (radius 0.2) using X-FEM and geometrical near-tip enrichment (radius 0.1, basis given by (6)); 3. T-stress method (radius 0.2) using X-FEM without near-tip enrichment; 4. T-stress method (radius 0.2) using X-FEM and geometrical (radius 0.1) near-tip enrichment.
The first and second cases allow one to study the influence of the quality of the finite element field on the convergence properties of the domain integral method. The results are presented in Figure 9.
As shown in Figure 9 Integral X-FEM (H only)'). In fact, this last computational setting corresponds to the one used in [41]. Finally, even with near-tip enrichment, convergence rates are different for the two components of the J-integral: it is optimal for J 1 but not for J 2 .
Consider now the convergence of the method proposed in Section 5.1 as presented in Figure 10. It is seen now that both J 1 and J 2 converge at the same rate (O(h 2 )) with a very low error level Table III. Summary of the convergence rates of J k -integrals.

Geometrical domain
Geometrical domain Geometrical T-stress integral X-FEM integral +bnd term X-FEM domain integral X-FEM even for coarse meshes. The figure also shows that convergence is obtained even without near-tip enrichment, which was not the case with the classical domain integrals. Results obtained with this benchmark are summarized in Table III, and they clearly exhibit the accuracy of the proposed approach.

CONCLUSION
In this paper, the accuracy of the numerical computation of path-independent J k -integrals using domain integrals has been investigated. It has been recalled that a boundary term (related to the near-tip stress field) had to be considered to evaluate these quantities. However, the lack of accuracy of the numerical scheme adopted for the integration of this boundary term can lead to a loss of convergence. In particular, J k -domain integrals do not converge without near-tip enrichment (contrary to what was assumed in [41]). Moreover, even if it converges, optimal rate cannot be obtained for J 2 . To overcome these limitations, an alternative approach based on the evaluation of the T-Stress has been proposed to avoid the integration of the boundary term. In this case, optimal rates of convergence are recovered for J 2 . However, this paper considered only the linear-elastic bidimensional case. Investigations are still needed to extend these results to non-linear elastic-plastic fracture, and to the 3D setting. In the latter case, the T-stress is no longer scalar, but becomes a vectorial parameter, which leads to some extra difficulties. In addition, it is important to note that the results presented here imply consequences in the context of configurational mechanics [55] (for generalities on this theory, see Gurtin [56], Maugin [57,58] and Steinmann [59, 60] among others). As a matter of fact, it can be shown (see Appendix A) that nodal configurational forces, introduced by Steinmann [59,60], are similar to the domain J-integral on the support of the node that contains the crack tip, but without the boundary term, which is mandatory for convergence. Thus, nodal configurational forces cannot converge to J.

APPENDIX A: RELATIONSHIP BETWEEN DOMAIN J-INTEGRAL AND NODAL CONFIGURATIONAL FORCES
Recall the expression of the kth component of the nodal material force acting on node I as depicted in Figure A1 [59]: where N Elem is the number of elements connected to node I. Next, we consider the kth component of the domain J-integral (without the boundary term) expressed on the support (denoted as V I ) Crack Figure A1. Configurational force acting on node I and support of the node.
of node I where the plateau test function is interpolated using the finite element shape functions Simply comparing Equations (A1) and (A5), we have J V k = −F s k . This demonstrates that the nodal configurational force is not the J-integral, because the boundary term J B k is missing (see Equation (17)).
This missing term corresponds to a configurational Neumann boundary condition on C : this boundary condition should not vanish as a traction-free surface in the physical space does not lead to a configurational-free surface.