Three‐dimensional simulation of crack with curved front with direct estimation of stress intensity factors

This paper consists of an extension of simulation with direct estimation of stress intensity factors to the three‐dimensional case. Here, it combines X‐FEM with localized multigrids and direct estimation of quantities of interest along the crack front (SIF, T‐stress, etc.) based on crack tip asymptotic series expansion. In practice, a three‐dimensional patch is introduced locally with a truncated basis of Williams series expansion and is linked in a weak sense with the X‐FEM localized multigrids. Some examples (with available analytical solutions) illustrate the efficiency and the robustness of the method. These examples consider planar cracks with curved front, but the proposed method aims to apply to any continuously curved crack. Copyright © 2014 John Wiley & Sons, Ltd.


INTRODUCTION
Finite element simulations of three-dimensional cracks in complex bodies have a convergence rate generally governed by the singularity [1]. The accuracy of such methods is enhanced by introducing an analytical basis able to capture the singularity; for instance, Tong et al. [2] built a super-element to locally introduce the singularity in FEM framework. Barsoum, in [3], proposed an even simpler method using quadratic finite element with nodes on the quarter of an edge. However, these methods need to have a mesh compatible with the crack surface. Therefore, heavy re-meshing [4] is needed to simulate propagation. This drawback has been tackled in the eXtended Finite Element Method (X-FEM) [5][6][7][8] by adding the crack discontinuity and singularity in the mechanical fields, thanks to FEM partition of unity property [9]. This method, very close to the Generalized Finite Element Method [10], allows to consider cracks non-conforming with the mesh that simplifies propagation simulations. Considering minor modifications [11], the X-FEM convergence rate for crack simulation sticks to the FEM convergence rate for regular problems.
Quantifying the crack tip singularity is physically useful in order to simulate propagation or failure in brittle materials or under small-scale yielding assumption [12]. This singularity can be decomposed on three modes (i.e., I the opening one, II the in plane shear one, and III the tearing out-of-plane one) that are meaningful for propagation velocity and orientation. A Stress Intensity Factor (SIF) is associated for each of these modes. These three SIFs are the first coefficients of a crack tip asymptotic series expansion introduced by Muskhelishvili and Williams [13,14]. Some other terms of this series expansion have be proven meaningful, such as the T-stress [15,16] or some super-singular terms [17,18]. In FEM and X-FEM, the stress intensity factors are usually recovered after the simulation by a post-processing method. The most accurate methods are based on the projection of the mechanical field on a surface or domain around the crack front. Arising from energetic approaches, path independent contour integrals, expressed as a domain/volume integral, have been proven efficient in mixed modes three-dimensional problems in an X-FEM context [6,[19][20][21]. Such a projection can also be applied to recover the T-stress as proposed by [22]. Besides, a method based on direct projection (in the least squares sense) of the displacement on the asymptotic series expansion was developed in experimental context [18].
Although X-FEM increases the mechanical field accuracy and enhances the convergence rate, the enriched degrees of freedom associated with the singularity are neither proportional to the energy release rate nor to the stress intensity factors. In References [23,24], Xiao and Karihaloo proposed a partition of unity (PU) based method that provides accurate SIFs for two-dimensional problems, where higher order terms are introduced (i.e., the enrichment is a truncation of the asymptotic series expansion) and constrained to have the same value in the crack tip vicinity. However, this method does not provide accurate estimation of the higher order terms. These terms (and the SIFs) can be obtained with the introduction of Hybrid Crack Element (HCE) [2,25,26]. In this approach, the crack tip discretization is a truncation of asymptotic series expansion alone.
Réthoré et al. [27] proposed a domain decomposition method where the crack tip vicinity is treated with the asymptotic series expansion only and provides two-dimensional stress intensity factors more accurately than the PU based approach. This domain treated with the asymptotic series expansion is coupled to an X-FEM domain with an Arlequin matching [28], but the coupling can be simplified with an integral matching and made compatible with a multigrid X-FEM approach [29,30]. Such an approach provides accurate evaluation of the SIFs (hence its name DEK-FEM for Direct Estimation of generalized stress intensity factors .K/ Finite Element Method) but also of higher order terms.
In this paper, this last method is extended to three-dimensional problems. In three-dimensions, no general series expansion exists for a crack problem; thus, the plane strain definition is generally considered in the bulk. For a planar crack with a straight front, such an approximation can be legitimate when far from the free surfaces. However, for general crack geometry, a common approximation in interaction integral methods [6,20,31] is to use the plane strain definition in local coordinates that follow the crack. In the present, such an approximation was not satisfying; the approximation investigated here is to derive the displacement series expansion in the local coordinate system. Moreover, the simulation of a complex three-dimensional crack requires local refinement that can lead to a high degrees of freedom number. To deal with this problem, the proposed method is coupled with the X-FEM localized multigrid approach proposed by Rannou et al. [32].

Two-dimensions development
The two-dimensional analytical study of the mechanical fields in a crack tip vicinity provides [13,14] first order (i.e., the singularity only) or higher series expansion. These asymptotic series expansions are presented in this section and then used in numerical simulations. [13,14]. Let us consider a straight crack in a linear elastic two-dimensional domain (with plane strain/stress assumption with respect to the crack plane .e n ; e t / defined in Figure 1). Under small perturbations hypothesis, the resolution of the equilibrium equation provides an expression of the stress field. If the crack faces are free . n D 0/, the stress can be developed [14] in r n=2 1 in the crack tip vicinity (i.e., r close to 0), in a polar coordinate system defined in Figure 1. From this stress series expansion, the displacement field can be computed "

The Williams series expansion
in the crack coordinate system defined Figure 1. is the shear modulus and Ä is the Kolossov's constant (function of Poisson's ratio ), which takes the value Ä D 3 4 in a plane strain state and Ä D .3 /=.1 C / for plane stress conditions. The asymptotic coefficients b n i define the stress and displacement fields and mainly depend on non-local effects (boundary conditions, etc.). b 0 I and b 0 II are in-plane rigid body translations, and b 2 II is the in-plane rigid body rotation. b 1 I and b 1 II are proportional to the SIF, and b 2 I is related to the T -stress, with In the crack tip vicinity, all the terms f n i are vanishing except the SIF (singular stress), the Tstress (uniform stress), and the rigid body translations. These terms (and especially the SIF) are of great interest to characterize the crack singularity and its influence on crack propagation (Paris, K IC ). On the contrary, the higher order terms have a far-field influence ./ r n=2 /.

Out-of-plane series expansion: mode III.
A tearing mode of rupture is added; it results from the same geometry with out-of-plane displacements u D u.r; Â/e 3 , where e 3 is the out-ofplane direction. The asymptotic study of such displacement under linear elastic fracture mechanics hypothesis provides series expansions of the displacement and stress fields [33]. They can be expressed as in Equation (2) Two out-of-plane rigid body motions are added-b 0 III , the out-of-plane rigid body translation, and b 2 III , the out-of-plane rigid body rotation. The coefficient b 1 III is proportional to the mode III stress intensity factor. 638 C. ROUX LANGLOIS ET AL.

Three-dimensions extension
The three-dimensional asymptotic study is much more complicated. Leblond and Torlai [34] proposed to account for the crack curvature and stress intensity factor evolution along the front as perturbations and computed some first order corrections to the classical Williams series expansion. Chaudhuri and Xie [35] also proposed an expression of the asymptotic displacement and stress field near the surface corner point. However, there is no three-dimensional series expansion available; thus, the two-dimensional plane strain modes of fracture are commonly used in three-dimensional crack simulations. The three modes of fracture are considered, and the asymptotic coefficients are expected to evolve continuously. They are generally written in local coordinates [6,20,31] without accounting for the curvature and evolution along the front. However, some recent studies show that accounting for the curvature in the definition of the auxiliary field can increase the SIF accuracy [21,36]. This approximation has been proven efficient for post processing methods such as the interaction integral. However, in the proposed approach, the asymptotic fields are introduced locally in the equilibrium problem; therefore, the displacement field is made compatible with the considered stiffness. A primal approach is considered, from the classical displacements in three-dimensions where s is the position along the front. Furthermore, the missing third rigid body rotation is added. From this statement, the evolution of h n i .Â; s/ along the curvilinear abscissa has to be chosen; a proposition of evolution is developed in Section 3.1.3. The asymptotic stress functions f n i cannot be used directly since they do not account for s variation; the corresponding displacement gradient is computed to obtain the stress field. This computation is complex since the local basis .e n ; e Â ; e s / is following the crack and the front, which are curved in the general case. This basis is illustrated in Figure 1 The displacement in Equation (5) is a function of local parameters .n; Â; s/; the computation of D (its derivative), with respect to these parameters, is therefore straightforward. D can be linked to the displacement gradient, thanks to the Jacobian matrix J . This gradient is needed in the global basis .e 1 ; e 2 ; e 3 /, in which the displacement is decomposed as u D V 1 e 1 C V 2 e 2 C V 3 e 3 r u D 2 6 6 6 4 The primal definition proposed here provides stress fields that are compatible with the displacement and verify the constitutive law but does not exactly ensure the local equilibrium. To the best of our knowledge, there is no three-dimensional fields satisfying all the local equations. The auxiliary fields of the interaction integral [6,20] suffer the same difficulty, and for curved cracks, the auxiliary stress field is not compatible to a displacement field. Nevertheless, this approximation is assumed and its accuracy studied on some classical benchmarks.

Mono-scale formulation
The considered problem is a material body in , a domain of the three-dimensional space, with well-posed boundary conditions: prescribed displacements u d on @ 1 and forces F d on @ 2 , with @ D @ 1 [ @ 2 and @ 1 \ @ 2 D ¿. A micro-structurally long crack is considered in this body under the assumption of small perturbations. The material is homogeneous, and its behavior is linear elastic and isotropic. Therefore, the state of the structure is defined by the displacement u alone. This displacement field belongs to the admissible displacement space The virtual work equation without body forces, equivalent to the governing equations for the cracked body as shown in [37], is given by where r s u is the symmetric part of the gradient of u, corresponding to the strain. H is the Hooke tensor. v is a test function in U 0 , the vector space corresponding to U with u d D 0.

Domain decomposition.
The crack tip vicinity is discretized with the asymptotic displacements of Equation (5). This approach allows for the direct estimation of the stress intensity factors (and the other asymptotic coefficients). Indeed, the associated degrees of freedom are directly proportional to the asymptotic coefficients. These functions are specific to the crack tip with the singularity, the discontinuity, and higher order terms. However, this development relies on strong hypotheses, a semi infinite rectilinear crack in three-dimensions (i.e., planar crack with straight front) in an infinite body. Therefore, the asymptotic fields can only be used in a small area, around the crack tip W , where boundary conditions are far and the crack almost plane. The subscript W will be used for this domain. The complementary domain X , such as D X [ W and ¿ D X \ W , is treated with the X-FEM. The subscript X will be used for this domain.
The Equation (7) can be generalized to such an arbitrary decomposition. In order to close the problem, the displacement continuity at the interface W , between X and W , is ensured by mortar approach as proposed in [29]. Such a matching is necessary because the considered displacements in X and W are not compatible. This condition is enforced in a weak sense with Lagrange The virtual work equation without body forces for the coupled two domains and the displacement matching can be written as follows, seeking .u Howewer, the asymptotic series expansion would not be able to handle directly most of the Dirichlet boundary condition; thus, situations where @ 1 W D ¿ are considered.
This equation is then numerically solved on the two sub-domains. The mechanical fields are discretized, using the X-FEM method on the domain X and the Williams series expansion on the domain W .

X-FEM domain X .
A classical three-dimensional extended finite element method, as introduced in [38], is used in the domain X . The domain is meshed, and a finite element interpolation is around the crack. The elements cut by the crack have their degrees of freedom N c enriched with the discontinuity generated by the crack The function H.x/ is a the generalized Heaviside function, representing the crack discontinuity. In the mono-scale approach with the analytical patch, the crack tip vicinity does not belong to X . Therefore, no singular enrichment is needed.
Let U X be the degrees of freedom vector containing u i and u 0 i . Equation (10) can be rewritten as a matrix u X D N T X .x/U X . The stiffness matrix K X is then 3.1.3. Analytical patch. The resolution of the considered problem on W is done, thanks to a Galerkin method where the mechanical basis is a truncation of the asymptotic series expansion, Equation (5). The upper-boundary of this truncation is referred to as n max . Since W is in the crack tip vicinity, the first order terms represent well the mechanical fields. The continuous evolution along the curvilinear abscissa s of the coefficients b n i of Equation (5) is discretized with one-dimensional finite elements where ' k .s/ are one-dimensional shape functions of order one, as represented Figure 2(b). b ink are the b n i coefficients at each of the k 2 N s nodes. The curvilinear abscissa s on the front is computed in the whole patch W . The border shape functions are extended to keep the partition of unity property. Thanks to this finite element decomposition, a continuous evolution of the asymptotic coefficients along the front is evaluated as degrees of freedom. However, a new parameter is introduced for the patch: the discretization increment along the front ds D s kC1 s k ; k 2 N s .
From this displacement field, the symmetric part of the gradient field is computed in order to compute the stiffness matrix K W such as Since the three variables r, Â, and s are uncoupled, the computation of the matrix D is easily performed.

Interface bounding.
The numerical treatment of the interface consists of bounding the two displacement fields U X and U W at their interface with Lagrange multipliers. This dual space is written as the linear span of some (here l) chosen shape functions L l The choice of the discretization of these multipliers influences the accuracy of the SIF and the convergence rate of the multigrid algorithm [29]. In two-dimensions, these considerations lead to the choice of a space composed of the non-zero normal Williams stress field on the interface and the Williams rigid body displacement modes. In three-dimensions, the Lagrange multipliers approximation space is chosen to be the patch displacement on the interface W W N L .x/ D N W .x/; 8.x/ 2 W . The discrete version of Equation (8) can be rewritten in the matrix form The construction of the two sub-domains is represented in Figure 2(a). It is realized from a FEM discretization of the whole domain. Then, the elements partially inside a radius r W of the patch corresponds to the analytical domain W . And the remaining elements define the X-FEM domain X and are directly used as X-FEM elements. The interface W is composed of the element faces between the two sub-domains. The standard Gauss quadrature is not well fitted to integrate discontinuous functions of the X-FEM; this difficulty is known from the beginning of the method [5] and can be overcome by subdividing the cut elements. The mechanical fields are also discontinuous and not linear in the analytical patch. The integration is performed with the finite element Gauss points used to build the two sub-domains. The discontinuity can be handled in the same way, and a higher order Gauss quadrature is performed to reduce the error introduced by the integration of asymptotic functions. The computation of the projectors C X and C W requires to integrate products of shape functions on the interface W . This interface is composed of three-dimensional element faces with high order Gauss quadrature.
Once discretized, the three-dimensional problem (9) is similar to the two-dimensional one [29]: find U X , U W , and ƒ such that 2 where F is the generalized force vector. However, in three-dimensions, K W is bigger than in twodimensions because of the discretization along the front and of mode III. It has been shown in two-dimensions that a patch radius of r W D 2 finite elements, and a truncation of the Williams series expansion to n max D 7, was sufficient for an accurate estimate of the SIFs [29]. The same assumption will be used in the three-dimensional case.

Multi-grid formulation
The simulation of a cracked body requires local refinement especially in three-dimensions and for complex crack geometries. Indeed, the stress is singular near the crack tip; local refinement is needed to increase the accuracy of the computed fracture parameters (i.e., SIF). If the crack is curved, the analytical patch has to be small enough to make the quasilinear crack hypothesis acceptable. Since the minimum patch size is linked to the surrounding X-FEM elements size, and as illustrated in Section 4, local refinement can be needed to match the required accuracy or where the crack curvature is important.
The mesh is locally refined with a X-FEM localized multigrid method, as introduced in [32], and made compatible with the analytical patch in a very similar manner than it was proposed in  two-dimensions in [29,30]. The whole crack front has to be in the finest grid with a surrounding layer of elements superior to r W . The analytical patch is introduced in the finest grid and interacts with this grid only.
Within the localized multigrid framework, the crack tip is included in the X-FEM grids. Indeed, the grids are overlapping as illustrated in Figure 3. Therefore, singular enrichments are also considered. The elements containing the front have their degrees of freedom N e enriched with the singularity in p r. The approximation space is then The functions j .x/ are enrichment functions designed to introduce the singular behavior (term n D 1 in Equation (5)) in the approximation space U X [5] j .x/ 2

Localized multigrid algorithm with X-FEM grids.
The principle of the multigrid approach is to solve the considered problem (introduced in Section 3.1) on a fine discretization through an iterative process on coarser grids. The idea is that the global displacement field, requiring many preconditioned conjugate gradient iterations on a fine grid, can be estimated in less operations on a coarse grid. A ratio of two between two successive grids has been shown to be quasi-optimal. The construction of the grids is done from an initial coarse mesh of the whole structure; each coarse element is subdivided in two and recursively (n g 1 times), until the needed refinement in the crack tip vicinity is reached. As represented in Figure 3, the coarser mesh is named M 1 and then two consecutive grids are named M g and M gC1 (respectively, the coarse and the fine). Since the refinement is local, a boundary g exists between the refined domain and the non-refined one (the boundary of the finer grid). The quantities Q on a grid g are referred with a g subscript Q g . The jB subscript of the quantity Q refers to the restriction to a zone where a finer grid exists; Q jB is padded with zeros such that it has the same size than Q.
The interpolation of displacements from a coarse grid to a fine one is performed with the prolongation operator P, and the return from fine to coarse one is performed with a restriction operator P T . This P operator is computed by a direct collocation method presented in [32]. Once initiated (step 0), the iterative algorithm (exponent i refers to the considered step) is described underneath in the general case for n g grids. It is composed of two steps: (1) resolutions from fine grids bring a corrective second member R to the coarser problem on which 1 conjugate gradient (CG) iterations called relaxations are performed and (2) 2 conjugate gradient iterations are done with the corrective residual term and the displacement is projected onto finer grid. The number of multigrid sub-cycles for each stage is taken to g D 1 8g 2 ¹n g 1; : : : ; 2º. 0. Initiation: the displacements are set to 0, g D n g ; c g D 1 8g 2 ¹n g 1; : : : ; 2º and R i n g D 0.
(1) Corrective residual term from fine grids, g D g 1 with boundary condition U i g j g 1 D P g U i g 1 j g 1 : (a) If g D 1 : go to 2. (b) If c g D g W c g D 1 and return to step 1 Else: c g D c g C 1 and g D g C 1.
(2) Relaxation of fine grid correction on coarse grids Else: (i) Interpolate displacements from coarse to fine mesh U i C1=2 (b) If g D n g W Return to step 1 Else: Return to 2 with g D g C 1.
This method allows to solve an equivalent problem on a refined mesh, represented in Figure 3, while only handling and solving smaller problems, either coarser or localized.

DEK-FEM coupled with the localized X-FEM multigrid method.
The analytical patch is inserted in the finest grid g D n g . Equation (22) is modified for the finer grid that includes the analytical patch. As introduced in [29,30], the X-FEM model of the finer grid exists in the patch volume but is inactivated at convergence by K n g jW U i n g , where K n g jW is the finer grid X-FEM stiffness matrix of the area overlapping the patch. This inactivation allows to decompose Equation (17) into an X-FEM one (coupled with the analytical patch) and a resolution on the analytical patch Patch resolution: The analytical patch influences the X-FEM model through ƒ i . In [30], the authors noticed that ƒ i C1 can be expressed as a function of the coarse problem unknowns U n g only where H W is known as the Steklov-Poincaré operator [39]. From this expression, the resolution of the finest grid including the patch is reduced to 2 This approach has been shown efficient in [30]. The patch resolution, Equation (23), is only treated at the end of the multigrid iterations to recover the Williams' coefficients. In practice, neither K W nor H W is inversed; the computation of H XW D C T X H 1 W C X is performed in two steps:

NUMERICAL SIMULATIONS
Some classical cracked configurations are studied here, and the identified stress intensity factors K DEK i ; i D I; II; III and T -stress are compared to former analytical and numerical solutions The crack geometries considered here are planar crack with either straight or circular fronts. For straight fronts, the Jacobian matrix J of Equation (6) is the identity. When the front is circular, the displacement gradient is easily constructed in cylindrical coordinates. This approach can be generalized, as a first order approximation to planar crack with continuously curved front using the local curvature of the front.
The circular cracks are considered first since analytical studies of elliptical cracks provide analytical SIF evolution. Most of these studies are performed in an infinite medium with remote boundary conditions, where the considered cracks are small by comparison with the considered domain (i.e., the crack radius a is about a tenth of the length-scale l). As it was done in [6,38], these two different length-scales are taken into account, thanks to a graded mesh. A penny-shaped crack under remote tension or shear loading is considered, and the evaluated SIF and T -stress are compared to analytical values. The crack is then inclined at 45°to activate all the fracture mode and to have variations along the front. Eventually, a corner edged crack under remote tension is considered. If not mentioned otherwise, the patch parameters used for the simulations are n max D 7, r W D 2, and ds ' 2 finite elements. In a second time, the localized multigrid approach developed in Section 3.2 is applied to a straight crack in order to illustrate the efficiency of the DEK-FEM with multigrids in three-dimensions (see [29] for the two-dimensional case)

Penny shaped crack under remote tension
As reminded in [40], under the Linear Elastic Fracture Mechanics (LEFM) hypothesis, the symmetric problem (hence mode I ) of a penny shaped crack of radius a under remote tension . 0 / in an infinite domain is subject to where is the Poisson's ratio. For the simulation, a cube of size 2l submitted to uniform tensile stress 0 is considered with an embedded penny-shaped crack of radius a D l 10 . Such a ratio is considered to match the infinite domain hypothesis. Two graded meshes with local refinement around the crack are considered (see Figure 4)-(1) A coarse one, as fine as in [38] in the crack tip area. For this mesh, there is 3.5 elements in the crack radius; therefore, the mesh size is H D a 3:5 .
(2) And a fine one, twice as fine. For this mesh, the crack radius is about seven finite elements W h D a 7 . As illustrated Figure 5, the fine mesh provides a good approximation .< 3%/ of the T -stress and a uniform mode I SIF. The mode I SIF estimation is uniform but presents an offset .' 12%/ with the analytical value. The mesh refinement increases the accuracy of the stress intensity factor (from ' 18% to ' 12%). On the other hand, the T -stress identification is not very sensitive to mesh refinement; indeed, its area of influence is larger than the singular mode.
The influence of the patch dimension r w is evaluated on the fine mesh with three different radius. The identified and projected SIF and T-stress are represented Figure 6. As it has been shown in twodimensions in [29], a patch of less than two elements of radius is too small for a proper identification of the asymptotic parameters. Indeed, the displacement on the border of the analytical patch W      is too poor (not enough degrees of freedoms) to make the two sub-domains communicate. On the other hand, when the patch is too large, it includes an area where the Williams expansion hypotheses are no longer valid and the quality of the identification decreases. Therefore, the dimension used is r w 2 finite elements.

Penny shaped crack under remote shear
A penny shaped crack under shear loading is considered. As illustrated in Figure 7 . This example is interesting since the non-symmetric modes II and III are activated and evolve along the front. An analytical value, also reminded in [41], exists for such a crack in an infinite body As displayed on Figure 7(b), the non-symmetric coefficients evolutions and amplitudes are in good agreement .< 5%/ with the analytical values even with the considered raw mesh.

Inclined penny shaped crack under remote tension
Let us consider the geometry and loading of Section 4.1 with an inclined crack (of angle D 45 ı with respect to the tension direction), as represented in Figure 8. This problem also has an analytical solution in the general case for an elliptical crack, reminded in [42], which can be reduced to the case of a circle This configuration is particularly interesting since the three modes are activated and both the shear mode and the anti-plane mode evolve along the front. The same two meshes as used in Section 4.1 are considered.  The evolution of the identified stress intensity factors are represented in Figure 9. The evolution of K DEK II and K DEK III is correctly identified, but the mode II amplitude is overestimated. For mode I SIF, the results are even worst since it oscillates strongly around the analytical value. The mode I oscillations and the K DEK II overestimation do not decrease with the mesh refinement and seem coupled with the other modes or other orders.
The mode III stress field computed from the displacement given Equation (12) is coupled with the Mode I . Indeed, as it was noticed in [34], the crack curvature and the evolution along the front (i.e., ' s .s/ in the present) bring out-of-plane stress for mode I and in-plane stress for mode III. A simple solution to reduce this coupling is to use a different discretization along the front for mode III For the sake of simple numerical implementation, ' 0 k .s/ is chosen to be a finite element decomposition of order 0 with the same nodes. The asymptotic coefficients identified thanks to this method are represented in Figure 10 for the two discretizations. For the fine discretization, the identified evolutions are satisfactory but an overestimation remains, ' 10% for K DEK II and ' 15% for K DEK III . On the other hand, K I is accurate < 4%, and the T-stress is oscillating as cos.2˛/. When the mesh is coarser, the 0 order discretization of mode III leads to uneven K DEK III , which follows the analytical tendency. For the three modes, the discretization enhances the accuracy. This uncoupling method 648 C. ROUX LANGLOIS ET AL.  (based on Equation (30)) works well in this case, and for all the previous test-cases. In the next parts, this approach will be considered.

Corner edge cracked bar
A square bar containing a corner crack loaded with remote tension 0 is now studied. The crack is perpendicular to the tensile loading as represented Figure 11(a). The bar has a square section of edge w, and its height is 2h D 2w. It includes a quarter penny shaped corner crack of dimension a D w 10 . This ratio is chosen in agreement with available closed-form solutions, reminded in [42] for such a corner edge crack This expression is reported to be accurate at 3% between =18 and 4 =9 and will be considered as reference. A graded mesh represented in Figure 11  crack radius. The T -stress intensity is compared to tabulated values proposed in [43] for elliptical cracks in plate. This reference is for a pretty different geometry, but the bulk value and the tendency near free surfaces of the T -stress have the same numerical behavior. The evolution of mode I SIF from the simulation, represented Figure 12(a), compares well with the reference. As noticed for the penny shaped crack under tension (Section 4.1), K DEK I presents an offset (an underestimation of about 3%). The T -stress on the other hand has the same negative intensity than the approximate reference (see Figure 12(b)).
At the vertices, it has been shown in [44] that the order of singularity is no longer 1 2 when the crack is perpendicular to the surface. Therefore, the considered asymptotic evolution is no longer accurate to describe the mechanical fields. However, a diminution of the stress intensity factor and and the increase of the T -stress are expected. The decrease of K DEK I is observed, and the T -stress tends to increase, expressing an increasing of the stress triaxiality.
The proposed description of the mechanical fields in the crack tip vicinity is an approximation which shows promising results. Yet, as illustrated in Sections 4.1 and 4.3, this approximation may require some local refinement near the crack tip to provide accurate asymptotic stress intensity factors. That is why the DEK-FEM strategy has been associated with an X-FEM localized multigrid approach. At this stage of development, this strategy is efficient for straight cracks. It is illustrated on the following example.

Single edge-crack tension specimen
A simple geometry, the single edge crack tension specimen is considered to evaluate the stress intensity factor identification method. It consists of a rectangular bar half cut by a single plane crack with a straight front. The considered loading is a unit normal stress, in order to activate the opening mode .I /. The considered geometry is defined in Figure 13(a); the width is W D a, the thickness t D 3a, and the height h D 3:5a. This example is also studied by Li et al. [41] with a boundary element method and by Sukumar et al. [38] with an extended finite element method analysis and a domain integral. These two previous study provides numerical references for K I .
As shown in Figure 13(b), the mode I stress intensity factor evaluated with DEK-FEM is close to References [38,41] in the bulk, with the same mesh in [38]. At the vertices, the stress intensity factors are less accurate but decrease as expected.
The mutligrid approach introduced in Section 3.2 has been developed for straight cracks and is tested here. Three simulations are presented, a monogrid one with cubic mesh of 24 36 42 elements 3 and two multigrid ones with two and three grids from a coarse grid of 12 8 21 elements 3 . As represented in Figure 14(a), the fine grids are localized around the crack tip in order to include at least one layer of elements around the finer grid and two layers around the analytical patch for the finest grid.   The multigrid simulation with two grids has a discretization as fine as the monogrid one around the crack tip and provides very similar K I in the bulk. The refinement provided by the third grid is more effective where the stress intensity factor evolves and increases the accuracy of the identification as expected. Another lever to improve the accuracy is the choice of the basis of the analytical domain.

CONCLUSION
In the proposed method, the simulation of three-dimensional cracked bodies with curved cracks and the identification of meaningful asymptotic coefficients are performed jointly. The asymptotic series expansion is used in the crack tip vicinity since it is shaped to represent well the singularity. The considered bases evolve continuously along the front, providing smooth asymptotic coefficients.

651
In three-dimensions, an underlying difficulty resides in the stress intensity factors estimate. Indeed, there is no general expression of the singular field for a three-dimensional crack with complex geometries. The two-dimensional asymptotic fields are used, but some local equations have to be approximated. Here, a method to construct these approximate fields is proposed. This method is shown to be efficient when the analytical patch is small enough with regards to crack size parameter. In order to have a patch as small as needed, a numerically efficient local refinement method is proposed.
This strategy is applied to numerous classical three-dimension benchmarks with known stress intensity factors, in which planar cracks with circular front are considered. This first study is promising; the evolution of the asymptotic coefficients (i.e., stress intensity factors and T -stress) is evaluated correctly with rather coarse meshes. However, the method aims to apply to any continuously curved cracks in further developments. Moreover, the higher order coefficients are also estimated as it was done in [29,45], but their accuracy is still to be considered.
The proposed method can easily be adapted to some other crack tip discretizations. With the one proposed, the convergence rate, with respect to the underlying X-FEM mesh, is still to be evaluated in detail. This method, based on the asymptotic Williams series expansion, is close to the one proposed by Lachambre et al. [46] to identify three-dimensional experimental asymptotic coefficients. Then, work is in progress to build a strong link between three-dimensional imaging and three-dimensional X-FEM in order to capture the behavior of three-dimensional cracks.