Y Jin

O A González-Estrada

O Pierard

S P A Bordas

Error-controlled adaptive extended finite element method for 3D linear elastic crack propagation

We present a simple error estimation and mesh adaptation approach for 3D linear elastic crack propagation simulations using the eXtended Finite Element Method (X-FEM). A global extended recovery technique [1] is used to quantify the interpolation error. Based on this error distribution, four strategies relying on two different mesh optimality criteria are compared. The first aims at homogenizing the error distribution. The second minimizes the total number of elements given a target global error level. We study the behaviour of these criteria in the context of cracks treated by an X-FE approach. In particular, we investigate the convergence rates at the element-level depending its enrichment type. We conclude on the most suitable refinement criterion and propose and verify a strategy for mesh adaptation on 3D damage tolerance assessment problems.

Introduction 1.Motivation

Accurate fracture fatigue predictions are crucial for the assessment of the in-service life of numerous industrial components. In industrial practice, linear elastic assumptions are often made. Based on this stress analysis, a crack is inserted within the component where the principal stresses are maximum, and this crack is grown under the fatigue loading this component is expected to experience. This raises a number of difficulties:

1. The mesh used for stress analysis is much coarser than a suitable mesh for fracture mechanics simulations, because of the large stress gradients in the vicinity of the crack front, due to the singularity of the solution along the front which requires the use of special elements [2].

2. This (fine) mesh must be regenerated as the crack evolves and must conform to the crack faces and crack front.

Given these two difficulties, a major challenge for damage tolerance simulations is to provide reliable solutions at the lowest computational cost, whilst requiring minimum user intervention. The eXtended Finite Element Method (X-FEM) [3] is an approach able to lift the second difficulty mentioned above as it allows modelling cracks without a conforming mesh and with much coarser meshes than otherwise necessary, through the use of specially tailored functions, able to reproduce the solution behaviour in the vicinity of the crack.

However, whilst the X-FEM provides certain simplifications in the treatment of propagating cracks, a number of hurdles remain for such approaches to be usable in practice, as noted in [4]:

• The crack path is in general unknown a priori. Therefore, a "fine enough" mesh is necessary within a possibly large 3D subregion of the component under consideration, where the analyst expects the crack to propagate.

• There exists very few methods to estimate the discretisation error committed when solving the 3D crack problem using the X-FEM, and none have been developed to optimise meshes used for industrial damage tolerance assessment. Although very recent progress has been made on this front [5,6], the stress intensity factors provided by the standard X-FEM are generally oscillatory on the crack front [4,7] The aim of this paper is to provide a simple approach to error-controlled adaptive simulations for the extended finite element method, so that the same mesh as that used for stress analysis can be used as a starting point for damage tolerance assessment, followed by an automatic adaptive mesh refinement to minimise the computational cost given a target error level. Our aim is, given a global measure of the error, to ensure the optimality of the mesh for a given refinement criterion. Our experience is that, in spite of the early belief that X-FEM would provide a tool for crack propagation without remeshing adaptation is a "marriage of convenience" because the adapted mesh need not conform to the geometry of the cracks as they evolve.

Aside of the X-FEM approach to solve fracture mechanics problems, it is worth noticing that mesh adaptation has recently been used to control crack nucleation and propagation during damage-driven simulations [8,9]

Previous work

We now provide a brief review of the literature and previous work done in the field of adaptive methods for fracture mechanics using the X-FEM.

The X-FEM uses the Partition of Unity Method [10] to enrich the classical basis of shape functions with so-called enrichment functions to describe the discontinuity of the displacement field across the crack surface as well as the singular stress field along the crack front. The X-FEM was shown to provide optimal global convergence rates [11][12][13] and, more recently, [5,6].

Although the X-FEM is known for its capacity to describe singularities within coarse meshes, the method shows significant mesh size dependence. For fracture mechanics problems in particular, as discussed in [14] for academic 3D fracture problems and in [4,7,15], [5,6] for industrial applications, the accuracy of the stress intensity factors computed along the crack front is strongly sensitive to the mesh size.

This work on three-dimensional crack propagation using the X-FEM motivated the development of error estimates aiming at the minimisation of the error in energy or on a quantity of interest, so that the mesh can be optimised for a given target error level. As for the standard finite element method, the error estimation tools for X-FEM are classified mainly into two categories: residual type and recovery type. In residual-based approaches, the goal is to compute the error in the satisfaction of the governing equations in each element by measuring the residual. In recovery approaches, an approximation of the "true" solution is built through a smoothing operation, under the assumption that the true solution is smoother than its (extended) finite element counterpart, as discussed in the seminal papers of Zienkiewicz [16,17].

The first attempts at recovery based error estimation for fracture mechanics simulations using X-FEM are due to [18,19] and [20][21][22], using a local recovery technique and [1], using a global recovery technique, also adopted in the present paper. The basic idea of these papers is to develop enriched recovery techniques, which are able to capture the special features introduced through partition of unity enrichment within the approximation field. The authors show that through such an enrichment, an optimal convergence of the approximate error to the exact error can be obtained, which is not possible using the original, unenriched recovery techniques of [16,17]. A similar idea was introduced in [23], where the analytical expressions of the pure mode I and II stress components are used in the construction of the enhanced stress field. González et al. [24] compared the generalized patch recovery technique of [20][21][22] with the extended moving least squares recovery approach of [18,19] and showed the importance of enforcing statical, kinematical and boundary condition constraints in the recovery procedure. The recovery approaches described above were extended to the framework of the Stable Generalized FEM (SGFEM) in [25].

Residual based error estimation aims at deriving an upper bound of the global error. Thanks to a solid mathematical basis, the resulting estimated error is guaranteed to converge at the same rate as the exact error [26,27]. Recovery based error estimators are thought to be preferred by practitioners due to the simplicity of their implementation and robustness, and their ability to provide excellent effectivities, i.e. approximate errors which are very close to the exact errors. We therefore chose such a recovery approach, which we implemented in the commercial software Morfeo Crack. It has to be acknowledged, however, that residual based error estimates are a perfectly adequate choice, which we are currently investigating.

Once an estimate of the approximation error is known, the mesh should be accordingly optimized to concentrate the computational cost where it is most needed. It is shown in [23,27] that by coupling h-adaptation strategies with error estimators, the number of degrees of freedom used to reach a certain global error level can be significantly reduced. However, in these two studies, which are based respectively on residual [27] and global recovery [23] error estimation, the mesh is refined when the local estimated error is larger than a pre-defined threshold. This basic approach involving only a local error parameter neither allows to capture the convergence rate of the global error norm nor to coarsen the mesh in the low-error areas. Therefore, the optimality of the mesh cannot be guaranteed.

Contributions and organisation of the paper

The present study aims at filling this gap and provide a simple, robust and efficient approach to adapt the mesh on the fly during 3D crack propagation by the extended finite element method and to investigate a number of adaptation criteria to define mesh optimality. In the literature, two such criteria are introduced to define the optimality of a mesh. For the first criterion [28], the optimality criterion is based on the uniformity of the error distribution. In [29], a mesh is considered as optimal with respect to a measure of the global error if a minimal number of elements is used. In the present study, both criteria are analysed for the first time in the framework of X-FEM.

A large number of studies describe numerical tools for 3D crack propagation simulations [30][31][32], but adaptive error estimation and mesh optimization during crack propagation have rarely been addressed, particularly in the context of X-FE approximations. In this work a new algorithm is proposed to perform error-controlled adaptive crack propagation simulations with the X-FEM. In particular, we investigate an issue which has never been tackled in this context, i.e convergence rates at the element level in various regions of the mesh, and as a function of the type of enrichment applied onto the element considered. For simplicity, as the method is implemented within our commercial code Morfeo Crack, the global recovery technique proposed in [1] (XGR) is used for error estimation. Any other error estimation technique suited to X-FEM could be used.

The paper is organized as follows: in Section 2, the theoretical bases of X-FEM and crack propagation techniques are briefly recalled. Section 3 presents XGR as well as the different mesh adaptation strategies. In Section 4, the mesh adaptation strategies are analysed on a two-dimensional fracture problem. The strategies that are found to be best-suited to X-FE approximations are applied in Section 5 to perform a realistic three-dimensional crack propagation simulation. Section 6 closes with some conclusions and directions of future work.

Crack propagation simulation with the X-FEM

This section presents the simulation framework valid for all mesh adaptation strategies and simulations performed in subsequent sections. First, the enrichment strategy within the X-FEM context is recalled. In addition, the description of the stress intensity factor computation and crack propagation technique are presented.

X-FEM for linear elastic fracture mechanics

The eXtended Finite Element Method (X-FEM) was initiated by Moës et al. in 1999 [3]. Even if dedicated initially to fracture mechanics simulations, the method has been used since then for a wide variety of applications in mechanics [7,33,34], nanomechanics [35] and biology [36]. Its first applications for industrial fracture problems were presented in Bordas and Moran [4] and Wyart et al. [37]. By respecting the partition of unity method [10], in addition to the classical finite element approximation, some nodes are enriched by special functions to represent particular features of the solution. In the context of fracture mechanics, two classical enrichment types are used:

• To represent the displacement jump across the crack, the nodes whose support is completely cut by the crack are enriched by modified Heaviside functions;

• To improve convergence rates, the nodes located within a given distance to the crack front are enriched with dedicated functions, which represent the near-front displacement analytical field. This deteriorates conditioning significantly unless an adequate remedy technique is used, e.g. a more global approach to enrichment as proposed by Agathos et al. [5,6], the vectorial enrichment functions [38] and/or a pre-conditioner [11].

Modified Heaviside functions h(x) takes a positive value +1 if point x lies above the crack and a negative one -1 otherwise.

Several expressions for the near-front functions have been proposed along the years with different impacts on the numerical error and/or the conditioning of the system [38][39][40]. Within this paper, the following functions are used:

F(r, θ) = [ √ r cos(θ/2), √ r sin(θ/2), √ r sin(θ/2) sin θ, √ r cos(θ/2) sin(θ)]. (1)
The final expression of the displacement field evaluated at a material point x thus reads:

û(x) = n I=1 φ I (x)a I + J∈J h φ J (x)h(x)b J + K∈J T φ K (x) 4 =1 F (x)c K , (2)
where φ I are classical finite element shape functions associated to node I and F is one of the four near-front enrichment function. Associated coefficients are a I , b J and c K , respectively. J h is the set of nodes enriched by Heaviside functions and J T is the set of nodes enriched by near-front ones. A standard Galerkin approach is used to determine all the coefficients. Stress Intensity Factors (SIFs) K i (where i corresponds to one of the three fracture modes) computed in this work are obtained from the numerical solution with the interaction integral method [41] written in a domain form as detailed in [30]. This approach guarantees pathindependence of the integral, as introduced by Eshelby [42] and Rice [43].

Crack propagation technique

All crack propagation simulations performed in this work are based on the simple Paris law [44], which writes:

da dN = C(∆K) m , (3)
where da/dN is the crack growth rate, da is the crack length increment, dN is the increment in number of cycles, C and m are empirical material parameters1 and ∆K is the equivalent stress intensity factor (such that mode-I energy release rate is equivalent to the real one) [45] between maximum and minimum load which writes, in absence of crack closure effect, as:

∆K = (K max I -K min I ) cos 3 (θ/2) -3(K max II -K min II ) cos 2 (θ/2) sin(θ/2) 2 + E * 2µ (K max III -K min III ) 2 1/2 , (4)
where

K max i (resp. K min i
) is the SIF i corresponding to the maximum load (resp. minimum), θ is the angle computed at each point of the crack front with the maximum hoop stress criteria, E * is the equivalent Young modulus (e.g. E * = E in plane stress conditions) and µ is the shear modulus. At a given point along the front of the crack, the propagation angle θ is assumed normal to the local maximum hoop stress.

Once the crack growth rates and propagation angles are computed all over the crack front, the implicit representation of the crack with the two level sets is updated according to [30]. We refer the interested reader to [4] for details on 3D crack propagation algorithms and implementation aspects as well as [46] for a review on 3D crack growth simulations using partition of unity methods.

Error estimation and automatic mesh quality control

In this section, the basic principles of the extended global derivative recovery introduced in [1] will be briefly recalled. Then different element refinement laws based on the estimated error will be introduced. At the end, we will describe the general algorithm used in the present work for adaptive error control in 3D crack propagation simulation with X-FEM.

Extended global derivative recovery

The a posteriori error estimator introduced in [1] is applied in the present study for simplicity, as it is implemented within the commercial code Morfeo Crack. The approximate (raw) strain is calculated by the symmetric gradient of the approximate displacement field (û):

ε (û(x)) = 1 2 ∇ + ∇ T ⊗ û(x), (5)
where ∇ denotes the gradient operator and ⊗ is the tensor product. In this estimator, the standard finite element approximation for the strain field is enriched with the following special functions [1] spanning the near-tip strain field:

Then, the smoothed strain field ε is expressed as

ε(x) = n I=1 φ I (x)d I + J∈J h φ J (x)h(x)e J + K∈J K φ K (x) 4 =1 G (x)f K , (7)
where d I , e J and f K are scalar coefficients that are determined by minimising the L 2 norm of the difference between the X-FEM strain field (ε( û)) and the smoothed strain field ( ε), over the whole simulation domain Ω. The corresponding functional is defined as:

Ω ε (û) -ε 2 dΩ. (8)
As in the standard recovery based error estimation [47], the smoothed strain field is compared with the original (raw) X-FEM strains (Eq.5). The difference between the original and smoothed strain fields is used as an error indicator to drive adaptive strategies.

Element refinement criterion

The mesh adaptation procedure can be defined according to two criteria: Equal-Distribution Criterion and Min-Number Criterion. The first criterion aims at reaching a certain level of global accuracy and equally-distributed local errors [47]. As for the Min-Number Criterion, its objective is to obtain a given level of accuracy with the minimal number of elements [28,29]. In this subsection, the theoretical bases of these two criteria as well as different mesh adaptation strategies based on these criteria will be presented.

Notations

Energy norm

The energy norm is defined as:

ε(û) Ω = Ω σ(û) : ε(û)dV , ( 9
)
where σ is the Cauchy stress tensor.

In the context of linear elasticity, this norm can be expressed as

ε(û) Ω = Ω (C : ε(û)) : ε(û)dV , ( 10
)
where C is the fourth-order Hooke's tensor.

Energy norm error

The smoothed strain field obtained by the extended global recovery approach XGR is denoted by ε. The error in the energy norm e Ω is estimated by εε(û) Ω :

e Ω = ε -ε(û) Ω = Ω (C : (ε -ε(û))) : (ε -ε(û))dV . (11)
Different rules are available in the literature to perform refinement in order to either distribute equally the error over the domain (see section 3.2.2) or to minimise the total number of elements to reach a given error level (see section 3.2.3).

Mesh adaptation process

Mesh adaptation is an iterative process. At each iteration, each element in the mesh is tagged for refinement or coarsening, depending on the value of the error indicator and the mesh adaptation criterion. Coarsening is possible only if using a local error parameter -see section 3.2.2. To start this mesh adaptation process, an updated mesh size h new i must be defined in each element i of the domain Ω. The element refinement rule for element i can be written in a general form as:

h new i = ξ i h old i , (12)
where ξ i stands, for element i, as the relative mesh size reduction. For the first order triangular or tetrahedral elements considered in this study, the mesh size h i is calculated from the element area s i (h i = √ 2s i ) or similarly from the element volume. Rules presented in subsequent sections aim at defining this parameter ξ i . MAdLib open source library is used to perform the mesh adaptation by mixing different basic operations as detailed in section 3.2.4.

Equal-Distribution Criterion

This criterion aims at reducing the global error and balancing the local error throughout the domain. These two objectives are described by two parameters, as introduced in Onãte et al. [48], namely the global error parameter ξ g i and the local error parameter ξ i . These two parameters are defined respectively by global and local absolute errors ( e Ω and e i ).

Global error parameter

ξ g i = e r e Ω . ( 13
)
e r is the required global error, which is defined as follows in the present study:

e r = max (η 1 e Ω , η 2 ε(û) Ω ) . (14)
The target global error after one iteration of the error-controlled mesh adaptation process is thus function of:

• the current global error e Ω , weighted by factor η 1 . The objective of this factor η 1 is to prevent a sharp increase of the number of elements between two iterations of the error-controlled mesh adaptation process.

• the target global error η 2 ε(û) Ω at the end of the mesh adaptation process; The global energy norm is recomputed after each iteration.

The role of the two above factors in mesh adaptation process will be further explained in Section 4 on the 2D Westergaard problem.

Local error parameter

ξ i = e r i e i , (15)
where e i is the current local error while e r i is the requested local error, which can be calculated by the following two strategies.

• Element based local error

The objective of this approach is to achieve an equal-distribution of the energy error norm among all elements. Thus, the required error norm is defined as:

e e r i = e Ω √ N , ( 16
)
where N is the total number of elements in Ω. Then the local error parameter calculated from e e r i is denoted by ξ le i .

• Volume based local error

The objective of this approach is to obtain the same error norm per unit area (volume). The required error is expressed as:

e v r i = e Ω Ω i Ω 1 2 , ( 17
)
where Ω i is the area (volume) of element i. Then the local error parameter calculated by using e v r i is denoted by ξ v i .

Combination

ξ i = (ξ g i ) 1 n ξ e i 2 2n+d . (18)
• Volume

ξ i = ξ g i ξ v i 1 n . ( 19
)
n is the degree of the interpolation function for u, and d is the dimension of the problem. In addition, Zienkiewicz et al [47] proposed a simpler expression of ξ i by combining ξ g i and ξ e i :

ξ i = ξ g i ξ e i 1 n . ( 20
)

Min-Number Criterion

An element refinement law based on the Min-Number Criterion is proposed by Coorevits et al [29]. It is shown in [29,49] that this criterion aims at minimizing the total number of elements N while reducing the global relative error θ Ω to a given level (θ 0 ):

θ Ω = ε -ε (û) Ω ε + ε (û) Ω = e Ω ε + ε (û) Ω . (21)
This constrained problem of finding all ξ i can be written as:

Minimize N = i 1 (ξ i ) 2 , with i ξ 2q i θ 2 i = θ 2 0 . ( 22
)
where q is the global convergence rate in the energy norm (e.g. q equals 1 for 3-noded triangle elements using first order Lagrange interpolants) and θ i is the local relative error computed over element i The constrained problem ( 22) can be solved by introducing a Lagrange multiplier λ:

Minimize i 1 (ξ i ) 2 + λ i ξ 2q i θ 2 i -θ 2 0 . ( 23
)
The solution is explicit and given by [29]:

ξ i = θ 1/q 0 θ 1/(q+1) i i θ 2/(q+1) i 1/2q . (24)

Mesh adaptation strategies

Five mesh refinement strategies based on the different laws presented previously are used to define the element size factor ξ i .

Strategy 1: uniform-distribution Combination of the global error and the local, element based, error (Eq. 18). Strategy 4: uniform mesh refinement The element size factor ξ i is only determined from the global error parameter, as obtained by setting the local error ξ le i to 1 in ( 18) or ( 19). This means that the size of all the elements are reduced by a uniform ratio.

Strategy 5: minimum number of elements (Eq. 24).

The overall algorithm of the error controlled adaptive crack propagation simulation is summarized in Fig. 1. Main iterative process is the crack propagation loop. For a given crack configuration within a propagation step, if the mesh has been adapted to satisfy the error criteria, the two level-sets representing implicitly the crack must be transferred. As level-sets are defined by nodal values on the old mesh, an interpolation field is constructed by using the same shape functions as those used for the displacement field. The level-set values at the new nodes can be obtained by a simple evaluation of this field. This means that implicit crack representation might slightly deteriorate in regions where coarsening occurs in presence of curved crack surfaces. Once this is done, crack propagation can be performed as described in section 2.2 by computing crack growth rate and direction before updating the implicit representation.

Within each propagation step, an error control loop is performed to adapt the mesh in order to satisfy the error criteria. In addition to solve the classical boundary value problem, the XGR one is also solved to obtained the smoothed strain field. From these two resolutions, the error in energy norm is computed as well as the required mesh size field. This mesh size field is constructed either globally on the whole computational domain, either locally on a predefined region on which the mesh adaptation is taking place. An optional smoothing step can be performed in order to avoid large discrepancy in terms of element sizes. If current mesh does not satisfies this required size field, the open source mesh adaptation library MAdLib is used [50] which can adapt meshes made with first order triangles or tetrahedra. This consists of a third iterative process (mesh adaptation loop) in order to reach the targeted size field, detect and eliminate badly shaped elements and keep the boundaries of the structure intact. Each iteration is made of the following operations: "edge collapse" to, if needed, decrease the number of nodes, "edge swap" to improve element quality, "remove slivers elements" by using various basic operators and finally "edge split" to reduce the element sizes. In [50], a complete description of the algorithm is given.

4 Error controlled mesh adaptation on the 2D Westergaard problem (Mode I)

In this section, the classical Westergaard problem (an edge crack subjected to uniform tension applied at infinity) is used to illustrate the different mesh adaptation procedures. Traction boundary conditions representing the exact stress state around the crack tip are applied on the boundary of a [2 × 2] domain, whereas the bottom left corner is restrained in both directions and the bottom right corner in the vertical (y) direction. The x direction is denoted by 1, the y direction by 2. The expression of the tractions imposed on the four (4) edges of the domain is given as follows:

T i = σ ij • n j , i, j ∈ {1, 2}, (25a)
σ 11 (r, θ) = K I √ 2πr cos θ 2 1 -sin θ 2 sin 3θ 2 (25b) σ 22 (r, θ) = K I √ 2πr cos θ 2 1 + sin θ 2 sin 3θ 2 (25c) σ 12 (r, θ) = K I √ 2πr sin θ 2 cos θ 2 cos 3θ 2 (25d)
where σ is the analytical stress field of the Westergaard problem and n j is the normal vector to the edge j. The expressions of the components in the x (1) and y (2) directions of σ are given by Eq. 25b-25d. The stress intensity factor of the first mode K I is set to be 0.001 in the present study.

For the simulations presented in this section, the interpolation degree m is 1 and a geometrical enrichment strategy is applied with an enrichment radius set to 0.1. All five strategies described above are applied to the first test case. The corresponding numerical strain field as well as the absolute energy norm difference between the numerical strain and the analytical solution obtained from the X-FEM on a mesh with uniform mesh size are given in Fig. 2. This (a) sketch of the problem on the initial mesh, the enrichment area is highlighted by the red circle; (b) the numerical strain (yy) using the mesh obtained from the first remesh increment; (c) the difference between the numerical strain and the analytical strain (yy).

section will be organized as follows. Firstly, the strategies based on the Equal-Distribution Criterion are compared with uniform refinement. Then the mesh adaptation sequences using Strategy 5, which is based on the Min-Number Criterion, are presented. These two approaches are compared in the third subsection. To conclude the section, the possible influence of the enrichment radius on mesh adaptation is discussed.

Equal-Distribution Criterion

Mesh adaptation procedure

The mesh adaptation sequences of the first four strategies are described in Fig. 3 and Fig. 4.

For the global required error e r , η 1 = 0.6 and η 2 = 0.02. The initial mesh is the same for all considered strategies and has 1134 first-order triangles. The numbers of elements of the meshes shown in Fig. 3 iterations are performed.

For the three strategies based on the Equal-Distribution Criterion (S1, S2 and S3 ), the finite element meshes are refined continuously for the first three increments. From the fourth increment, the total number of elements stops increasing. The error in energy norm weighted by the energy norm is defined as:

Err = e Ω ε(û) Ω . ( 26
)
Its evolution is illustrated in Fig. 5. The relative error in the energy norm (energy norm error weighted by the energy norm) on the initial mesh is about 0.1. As illustrated in Fig. 5(a), the relative error is reduced by 40% in each of the first three increments as defined by the mesh size evolution ratio (η 1 in Eq.14). Once the error reaches 2% of the total strain energy, which corresponds to the target global error level (η 2 in Eq.14), η 2 ε Ω is necessarily bigger than η 1 e Ω . Consequently the global error parameter ξ g i equals 1. Then, the mesh adaptation law is determined by the local element error parameter ξ le i , which is used for homogenizing the error distribution. Accordingly, the mesh starts to adapt in the regions where the error is highest e.g. the crack tip and the boundary of the tip enrichment area. This evolution is consistent with the results in the literature [48,51]. Meanwhile, using Strategy 3, the total number of elements increases even faster than uniform refinement during the first iterations, which is not consistent with the results in literature obtained with the classical finite element method. A possible explanation comes from the lower effectivity index we obtain, as described in the next paragraph. In addition, for the second strategy, due to the difference between global and local convergence rates, re-and derefinement oscillations are observed in the mesh adaptation process, as reported by [48] (see Fig. 6).

The global effectivity index, which equals the ratio between the estimated error ε-ε(û) Ω and the real error e Ω [1], is expressed in Fig. 5(b) as a function of the mesh adaptation increment. Generally speaking, the estimated error provides an overall accurate description of the real error. For S1, S2, at the beginning of the mesh adaptation process, the effectivity converges towards 1 as the mesh is refined, which is consistent with the results in [1]. After the target error level is reached, the mesh is adapted locally to homogenize the elementary error level and the average mesh size stays at a constant level, which leads to a constant overall effectivity. As for S3, despite a sharp increase in number, the effectivity index for the first iterations is less accurate compared to S1 and S2. Moreover, the non-optimized mesh adaptation strategy apparently jeopardizes the convergence of the estimated error. The global effectivity is thus not solely determined by the average mesh size. The distribution of the 6) of Figure 3). Re-and de-refinement oscillation is highlighted by the white framework.

norm error has the same magnitude as the degree of interpolation p, since

e Ω = Ω CO (h p ) O (h p ) dΩ 1 2 O (h p ) , ( 27
)
where h is the average mesh size. Meanwhile since the volume (surface) of element i depends on h i , the elemental error norm converges faster,

e i = Ω i CO (h p i ) O (h p i ) dΩ 1 2 O (h p i ) Ω 1 2 i O h p+ d 2 i . ( 28
)
This difference in terms of convergence rate should be considered for finite element simulations.

X-FEM was shown to be able to result in an optimal global convergence rate for problems involving singularities [11,12] using a fixed (independent of the mesh size) enrichment radius. This comes at the price of ill-conditioning of the system matrix, which can be addressed using more "global" enrichment based on cut-off functions or integral matching (see [13,52,53] and [5,6] for 3D crack problems).

The mathematical proof of this optimality can be found in [13,52,53]. In [13], the authors prove that with the so-called "cut-off" function, the discretisation error associated with the regular part of the displacement approximation converges at the optimal rate both globally and elementarily. However, the behaviour of the elementary error of the fully discretized displacement field including both the regular and the singular part is not addressed in [13].

In an enriched element containing the crack tip, the theoretical strain tends to infinity and the element error norm may not converge at the optimal rate due to integration errors of the singular function. The convergence rate and error level may also depend on the position of the crack tip within the tip element. In Fig. 8(a), the element error norm ( e i ) in the element containing the crack tip using the uniform mesh adaptation strategy is expressed as a function of the average mesh size h. As mentioned above, the element partitioning technique is used for integration in the crack tip elements. However, as reported in [38], this technique leads to oscillations in the simulation accuracy with respect to the relative position of the crack tip in the tip element. To avoid this oscillatory behaviour, a structured array of background cells can be used for integration. Other integration approaches were proposed by Ventura et al. [54] and Natarajan et al. using boundary integration [55,56] and stabilized conforming integration in Bordas et al. [57,58]. In addition, only one configuration where the crack tip is located on an inclined edge is considered. Since the crack tip is shared by several sub-elements, the energy norm of the error in all these elements is taken into account. The configuration used is illustrated in Fig. 7.

Local versus global convergence rates Contrary to the elemental convergence rate in energy norm in the finite element framework (n + d 2 = 2; n is the interpolation degree and d is the problem dimension), the convergence rate in this element is suboptimal and approaches unity (1), i.e. the global convergence rate in energy norm. However, only the tip element has a different elementary convergence behaviour. The convergence rate in the energy norm for an element within the enriched zone equals 2 see Fig. 8(a).

Figure 7: Crack tip (the red point) on an inclined edge in a structured mesh Besides, the estimated errors converge identically as the exact error (Fig. 8(a)). The effectivity indices are shown as a function of the mesh size h in Fig. 8(b). Mesh refinement allows for the estimation of both the global and the elementary errors. However, in the crack tip element, the effectivity converges slower than in other enriched elements, which shows that the estimated error is less accurate in this element. The lack of accuracy is likely induced by the singularity of the strain field, although this remains to be proven.

For Strategy 3, Eq.17 is obtained under the condition that [48]: We observe that optimal convergence rates are retrieved. Best effectivity index is obtained for elements outside the enrichment area.

e i Ω 1/d i O (h n i ) . (29)
However, as explained above, in the element containing the crack tip,

e i = O h n+ d 2 i
. Thus Eq.29 is no longer valid,

e i Ω 1/d i O (1) , (30)
which means that the error per unit volume (surface) in the crack tip element does not converge quadratically as the mesh is refined. This inconsistency results in the degraded performance of Strategy 3 in the present study. It is worth noticing that even though the real convergence rate of the crack tip element in X-FEM is different from the theoretical elementary convergence rate in FEM, the performance of Strategy 1 does not seem to be affected. This point will be further discussed in the third part of this section.

Min-Number Criterion

The mesh adaptation sequences of the fifth strategy are given in Fig. 9. Using Strategy 5, the mesh adaptation process is similar to the ones based on the Equal-distribution criterion. The mesh is refined for the first two increments, then it is optimized automatically by minimizing the total number of elements. In this strategy, a fixed relative error (θ 0 ) value is targeted.

In the present study, θ 0 is set to be 0.01, which means the required global error should be 1% of ε + ε (û) . Therefore, unlike the error evolution shown in Fig. 5, using Min-number criterion, the error drops sharply from the first increment to reach the aimed level as illustrated in Fig. 10(a). Once the target error level is reached, the mesh adapts the high-error areas by decreasing the total number of elements, thanks to the minimization approach adopted by the criterion. Moreover, the effectivity index remains nearly constant after a sharp increase in the first adaptation increment as shown in Fig. 10(b).

It must be underlined that since only the global convergence rate is used in this strategy, the sub-optimal convergence observed in section 4.1.2 for the tip element is avoided.

(1) 10232

(2) 14607

(3) 14672 ( 4) 14244 ( 5) 13960 (6) 13852 The target global error level is nearly reached after the first mesh adaptation increment and followed by mesh optimization.

Comparison of Equal-distribution criterion and Min-number criterion

The energy norm relative error, which is defined by Eq.26, is plotted as a function of the element size in Fig. 11. Except for the uniform refinement (Strategy 4 ), only the results of the last increment are shown for the four other strategies, since the average mesh size does not decrease continuously during the mesh adaptation process. Five meshes leading to a comparable error level (highlighted by the blue rectangle in Fig. 11) are chosen to further compare the performance of these five strategies. The error in the stress intensity factor of the first mode as well as the effectivity index of the error estimation on these meshes are given in Tab.2. Relative error in K I is obtained by the ratio of the difference (in absolute value) of numerical and exact K I over the exact value. Err), the relative error in K I and the effectivity on the meshes leading to a similar error level. For strategies 1, 2, 3 and 5, the results are obtained on mesh (6) as illustrated in Figs. 3 and4. As for S 4 , mesh (4) of Fig. 4 is used. We observe that, with a limited number of elements, S1, S2 and S5 have the lowest relative error in K I .

With strategies 1, 2 and 5, only one third of the elements are required compared to a uniform mesh refinement to reach the same level of accuracy. In an extended finite element framework, using geometrical enrichment, the convergence rate error in energy norm within the element containing the crack tip equals 1 instead of 2 in the standard finite element method. Meanwhile, Strategy 1 and Strategy 2 only consider a uniform elementary convergence rate. The inadequacy of Strategy 1 is thus localized in the crack tip element, which does not have an important effect on the global performance. As for Strategy 2, the re-and de-refinement oscillations do not seem to interfere with its efficiency in decreasing the global error. However, Strategy 3, which aims at homogenizing the error over the whole domain, is even less efficient than uniform refinement. Theoretically, the volume based error defined by Eq.16 in the crack tip element does not converge as explained previously. Consequently, the crack tip is "overadapted" with this approach, which leads to much more degrees of freedom than the other strategies.

Figure 12: The absolute energy norm error distribution of the simulation performed on the mesh obtained via Strategy 1 -6th iteration (see Fig. 3). Uniform error distribution is observed.

The errors in the stress intensity factor (K I ) on these five meshes are at a very low level (less than 0.1%). However, it is worth underlying that they are not correlated to the energy norm error as illustrated in Tab.2. The mesh adaptation strategies are defined to decrease the energy error norm.2 . In the present study, the J-integral is calculated within a circle with the radius equal to 0.15. Since the elements in the integration domain do not necessarily carry the highest error values in the simulation domain, they are not refined as much as other elements during the remeshing procedure, e.g. the elements within the enrichment zone distant from the crack tip. The distribution of the energy norm error throughout the domain3 S1 (6) is illustrated in Fig. 12. A rather uniform error distribution is observed. At the same time, Strategy 1 and 5, which are most adapted to the X-FEM, provide higher accuracy in the calculation of K I than the other strategies. Using these two strategies, the mesh is better adapted within the high-error zones, which results in a more homogeneous error distribution for a given global error level (about 2% in the present case). This uniform distribution can indeed improve the consistency of the calculation of K I , in particular its integration domain independence.

In addition to the simulation accuracy, the effectivity is also calculated for the different strategies (Tab.2). The convergence of the effectivity with different strategies are presented in Fig. 5 and 10. We fix the error level and compare uniform mesh refinement to our three adaptive strategies (strategies 1, 2 and 5). We note that for all three adaptive strategies, the error level is reached for fewer elements and that better effectivity indices are obtained compared to uniform refinement. The corresponding local effectivities of these three strategies as well as uniform refinement (Strategy 4 ) are shown in Fig. 13. It is clear from this Figure that all four approaches provide accurate error estimation on most elements. Note that the uniform mesh refinement approach (S4) leads to the most uniform effectivity distribution. For the latter, the error estimation is slightly less reliable in the zones where the error is highest, i.e. the crack tip and the boundary of the enrichment area, which is consistent with the analysis performed in section 4.1.2. With all three mesh adaptation approaches, we notice that the mesh is adapted within these high-error zones which leads to a better estimation of the error level, thus a better effectivity (closer to the optimal value of unity 1.0). For these adaptive meshes, which are also non-uniform, it is clear that the effectivity of the estimators is poorer in the larger elements in the mesh. elements locate is also relatively low, as shown in Fig. 2. Partial conclusions From the above discussion, we can conclude that the three mesh adaptation approaches we developed increase the reliability of the error estimator in the elements with high error level at the expense of a worse estimation in the elements where the error is already low in the initial, unrefined mesh. On the whole, the error is better estimated with a sequence of adaptive steps than with uniform refinement. At the same time, for Strategy 2, the oscillations occurring during mesh adaptation (successive refinement-coarsening operations) jeopardise the mesh quality, which could also have an effect on error estimation4 . For the purpose of quantifying this effect, the area where the relative estimated error is lower than 10% is calculated for each strategy. This calculation is performed on a circular area around the crack tip with a radius equal to 0.5. For strategies 1 and 5, this area (weighted by the area of the domain) is about 58%. But it drops to 50% for Strategy 2. This can explain why this strategy leads to lower global error estimation accuracy (see Tab.2).

Conclusion on the suitability of the proposed strategies From the above analysis, we see that none of Strategy 1 or 2 is adapted to the X-FEM. As mentioned previously, in the X-FEM context, the element convergence rate in the crack tip element (1) is lower than the rate in other elements (2). This difference leads to the inadequacy of the first two strategies in element size calculation. More precisely, Strategy 1, which uses an uniform element convergence rate to calculate the new element size, does not count for the difference between the crack tip element and other elements in terms of convergence rate. In Strategy 2, only the global convergence rate is applied to calculate the element size.

To summarize, for X-FEM simulations, the global convergence rate equals 1. In the crack tip element, the element convergence rate is same as the global one, while in other element the element energy norm error converges more rapidly (rate is equal to 2). Based on this observation, we propose a hybrid strategy, which addresses the speciality of X-FEM. Strategy 2, which does not distinguish element and global convergence rate, is used in the element containing the crack tip, and Strategy 1 is used elsewhere. The mesh obtained after six adaptation iterations by this hybrid strategy (Fig. 14) is almost the same as the one obtained by Strategy 1 (see Fig. 3). Only a slight increase in element number is observed. Since the exact element convergence rate is used (see Fig. 8 (b)), the mesh is adapted better around the crack tip with the combined strategy than with any of Strategies 1 and 2. As illustrated by Fig. 15, the size of the crack tip element is about a half of its value in Strategy 1 (the smallest element size is 4.61•10 -4 for the hybrid approach compared to 9.62•10 -4 in Strategy 1 ). Accordingly, the highest energy norm element error, which corresponds to this element, is reduced by about 50% (the maximum strain difference indicated in Fig. 14 (b) and Fig. 12 i.e., 5.97 • 10 -7 vs 1.29 • 10 -6 ). However, this local improvement does contribute to the global calculation accuracy. The overall energy norm error on this mesh equals 1.93%, which is very close to the value on the sixth mesh of Strategy 1.

Influence of the enrichment radius

The enrichment radius is increased from 0.1 to 0.3 to test its influence on mesh adaptation. In the previous section, evidence is shown that Strategy 4 using Min-number criterion is a suitable approach when used in the context of X-FEM. At the same time, despite the local incompatibility in the crack tip element, Strategy 1 also leads to acceptable optimised meshes after approximately six adaptation iterations without compromising the simulation accuracy. In this section, these two strategies are further compared with uniform refinement in terms of effectivity and efficiency. In particular, we consider the influence of the enrichment radius. The Westergaard problem is firstly solved on the uniform mesh (S4 (4)), since the simulation performed on this mesh has a similar accuracy compared to those performed on the mesh obtained at the sixth adaptation iteration of Strategy 1 and 5 for the first test case. The relative error in energy norm is 0.0118, which is about 60% of the error with the smaller enrichment radius (0.1). Accordingly, η 2 and θ 0 , which represent the target error levels for Strategy 1 and 5, are reduced to 0.012 and 0.006, respectively.5 The meshes obtained with these two strategies after six adaptation increments are shown in Fig. 16. It is noticed that some local refinement still occurs outside the transition region between the enriched and standard part of the computational domain6 . The corresponding energy norm relative error, the error in K I as well as the effectivity are given in Tab.3. Strategy 5 leads to meshes which are better refined within the high-error zones i.e. the crack tip and the boundary between the enriched and non-enriched zones. This approach leads to fewer elements. However, just as in the first test case, Strategy 5 is not able to reduce the error in the energy norm to the same level as Strategy 1 does. Nevertheless, it is worth noticing that the difference between these two approaches in terms of total number of elements as well as other indicators indicated in Tab.3 is less than 5%, which can be considered as negligible given the computational gain provided by both these approaches compared to uniform refinement. The enrichment zone is highlighted by the blue dashed line. We notice that the local refinement still occurs along the transition region between the enriched and standard areas. This is believed to be due to the lack of blending correction used in our computational approach.

To reach the same error level, both adaptive strategies require indeed much fewer elements than uniform refinement. Furthermore, the increase in the number of degrees of freedom driven by the increase of the enrichment radius in these two strategies is smaller than the increase in the number of degrees of freedom for a uniform mesh thanks to an adequate adaptation of the high-error areas. Therefore, the adaptive strategies are even more advantageous with respect to the uniform refinement when one wants to extend the enrichment zone to reach higher precision 7 8 .

Meanwhile unlike the energy norm relative error, the error in the stress intensity factor (K I ) does not decrease compared to the values in the first test case. In the present study, the J-integral is calculated within the disk of radius 0.15. Since the enrichment radius is increased from 0.1 to 0.3, the integration domain in the second test case is not fully refined by the adaptive strategies. Therefore the precision of the stress intensity factor calculation degrades9 . 7 It should be noted that enlarging the enrichment radius also leads to poorer conditioning of the system matrix. This can be alleviated using cut-off or integral matching enrichment schemes, see for example [11, 12] in 2D and [5, 6] in 3D, or by using "stable" versions of the X/G-FEM [39].

8 The cost of computing the error and element sizes to use in the adaptive scheme should obviously be taken into account. The global approach we use, for practical implementation reasons within a commercial code Morfeo Crack, is expensive. This is a drawback of the approach outlined in this paper, but we reiterate that the adaptive approach can be used with any (local or global) error indicator able to give an indication of the error distribution within the domain.

Discussion

None of the four strategies based on the Equal-Distribution Criterion is fully adapted to the X-FEM. This inadequacy results from the fact that the energy error norm in the crack tip element converges at a different rate with respect to the other elements. Despite this, Strategy 1, which aims at homogenizing the elementary energy norm error, requires fewer elements for a given error level compared to uniform refinement because of its ability to refine only in the high-error areas. In addition, we show evidence that by combining Strategy 1 and 2, using similar number of elements, the local energy norm error in the crack tip element can be effectively reduced. However, this hybrid strategy results in an error level that is very close to the one obtained by Strategy 1. This further demonstrates that the inadequacy of the first approach is limited to the crack tip element, which does not have a large effect on the global performance of this technique.

On the other hand, the Min-number criterion is fit for the X-FEM, since only the global convergence rate is applied. The two strategies (Strategy 1 and 5 ), which are well adapted to the X-FEM, result in a more precise calculation of K I in the two test cases compared to other strategies. Since the high-error areas are adequately refined using these two strategies, the error is more homogeneously distributed. Accordingly, for a given global error level, a uniform error distribution leads to more consistent values of K I . It is worth mentioning that using the Equal-Distribution Criterion, for the first few increments, the mesh is progressively refined to gradually decrease the global error; this being controlled by the mesh size evolution ratio η 1 . However, applying the Min-number criterion, the mesh is aggressively refined to reach the targeted global error level in the first adaptation increment. Then mesh refinement adjustments are made within the high-error areas. Consequently, the average number of elements affected by refinement at each mesh adaptation increment is higher, which increases the numerical cost. As illustrated by the two test cases, by increasing the enrichment radius, the global energy norm error can be effectively decreased. However, since the calculation path of the domain integral used to compute the stress intensity factor is defined independently of the enrichment domain, the elements within this domain are not necessarily refined by the adaptive approaches, which aim at reducing the energy norm error, not the error on the quantity of interest K I . This is evidence that the accuracy of the stress intensity factors cannot be correlated with the global error level. Therefore, a new adaptive technique based on a goal oriented error estimation is being considered in the future to address this issue [60,61].

5 Error controlled mesh adaptation for 3D Crack propagation

As illustrated by the 2D test case, Strategy 1 and 5 are more adapted to the X-FEM. In this section, these two strategies are used to perform a 3D crack propagation simulation of a four-point bending test. The initial configuration is given in Fig. 17. A pre-cracked rectangular right parallelepiped-shaped sample with a capillary. The size of the sample is [110mm, 20mm, 12mm]. Its center is located at (0mm, 0mm, 0mm). The crack, which is a quarter of an ellipse, is illustrated in Fig. 17(b). The center of the ellipse is at (1.5mm, -10mm, 6mm). The major and minor radii equal are 3.4mm and 2.1mm respectively. The major axis is oriented along the z direction and the plane of the initial crack is assumed orthogonal to the (x-z) plane.

The total force imposed on each side of the top surface is 6.5kN. Typical material properties corresponding to Titanium alloys are used (E = 122.5GPa and ν = 0.342). The mesh is made up of linear tetrahedral finite elements. An enrichment radius of 0.15mm is used around the crack front. For the crack propagation simulations, Paris law is used: C = 2.8977•10 -18 GPa -n •mm 1-0.5n and n = 4.76 (experimental validation on Titanium alloy; unpublished work). The load ratio equals 0.1 and the maximal crack propagation increment is set to be 0.2mm. The initial mesh used for crack propagation simulation with mesh adaptation is illustrated in Fig. 18. The mesh is refined around the crack and the capillary initially. The element size in the refined area equals 0.5mm while the coarse mesh size is 5mm. In this study, we focus on the performance of the mesh adaptation tool around the crack. Accordingly, the mesh adaptation is only performed in a parallelepipedic domain around the crack and crossing the sample as highlighted by the yellow domain in Fig. 18. The crack propagation path as well as the mesh adaptation process of the two strategies are described in Fig. 19 and Fig. 20 respectively. For Strategy 1, the error control parameters are η 1 = 0.9 and η 2 = 0.17. For Strategy 2, θ 0 = 0.07. To limit the numerical cost, for each crack propagation step, three mesh adaptation increments are performed. With these two strategies, the crack paths are nearly identical. After 32 steps, the crack intersects the capillary. The mesh is adapted around the crack front during the crack propagation. As the crack becomes larger, more elements are used around the crack front. The number of elements are given below the corresponding mesh in Fig. 19 and 20. As in the 2D simulations, the total number of elements required by the two strategies are comparable. Cuts of the mesh at cross section x = 1.4mm, after 32 crack propagation steps, for each of the two strategies are compared in Fig. 21. Overall, these two approaches lead to very similar meshes. The mesh size around the crack front approaches 0.05mm.

To further evaluate these two strategies, the results are compared with a simulation performed on a fixed mesh. The fixed mesh is refined around the capillary and the crack. The refined mesh size is 0.1mm. This mesh size is smaller than the enrichment radius and the maximal crack propagation increment to guarantee the simulation accuracy. The coarse mesh size is 5mm. The simulation results on this fixed mesh after 32 steps is given in Fig. 22. The crack fronts after 30 steps (before reaching the capillary) in these three simulations are compared in Fig. 23. The crack front shapes obtained by the two mesh adaptation approaches are very close to the one on the fixed mesh. The two fronts on adapted meshes are slightly closer to the capillary. This similarity is further confirmed by quantitative analyzes. The crack length is quantified by two parameters: the height (y direction) of the extremity of the crack on the front surface and the horizontal distance (z direction) on the bottom surface. As shown in Fig. 24 (a), these two characteristic lengths evolve identically as a function of the

Step 1 (23749)

Step 10 (51864)

Step 20 (125031)

Step 32 (296055)

EqualDist

Figure 19: Crack propagation after 1, 10, 20 and 32 propagation steps with mesh adaptation applying the Strategy 1 ; the total element number is indicated below each corresponding mesh.

The finite element mesh is adapted around the crack front as the crack propagates. As the crack front becomes longer, a larger number of elements are required to follow the crack and minimise the error.

Step 1 (21112)

Step 10 (47866)

Step 20 (114705)

Step 32 (276692) The cumulative number of load cycles is compared in Fig. 24(c) and the crack length evolution is expressed by the number of cycles in Fig. 24(d). The two mesh adaptation strategies lead to same number of cycles after 32 propagation steps. The difference from the result on the fixed mesh is less than 5%. The two mesh adaptation strategies result in similar crack paths and fatigue behaviour compared to the fixed mesh. Applying these two approaches, the mesh after 32 propagation steps has less than 300, 000 elements in total, while there are about 1.7 • 10 6 elements in the fixed mesh. However, the numerical cost of mesh adaptation and error estimation cannot be ignored. In 3D, to calculate the smooth strain field through XGR, the system to solve is three times as big as the finite element problem in terms of number of degrees of freedom. As mentioned, within each propagation step, three mesh adaptation iterations are performed to adapt the mesh, which means that the XGR system needs to be solved three times as well. At the beginning of crack propagation, this process is numerically less expensive compared to the cost induced by using more elements in the fixed mesh. However, as the crack propagates, the crack front becomes longer and more elements are required to represent it. Consequently, the cost of the mesh adaptation process becomes larger mostly due to the calculation of the smooth strain field. It is worth mentioning that the mesh adaptation strategies analysed in the present study can be used with any other approach to error estimation, in particular local approaches such as [18, 20] and [27, 59]. The application of a local recovery technique [22] or a residual error estimator [27] can be a possible route to reduce the numerical cost of the error controlled adaptation process, which we are currently investigating. It should also be noted that beyond driving mesh refinement, error indication can be used to provide the engineer with confidence on the quality of the results, which is an important aspect of error control.

MinNum

Furthermore, as explained in the previous section, to assess the error of the J-integral calculation more accurately, a goal-oriented error estimator is necessary. An error estimator based on the error in the J-integral can be more efficient and effective to improve the accuracy of crack propagation simulations compared to the one based on the energy norm error. Using such dual weighted residual approaches, it would also be possible to compute meshes which are optimal for the computational of the crack driving force and crack propagation direction, both within discrete crack models and phase field models. Development of an appropriate goaloriented tool is an important point that we are addressing within commercial code Morfeo Crack based on the work of [59,61]

Conclusions and future work

In this study, we propose an automatic error estimation and mesh adaptation procedure for 3D crack propagation simulation with the extended finite element method (X-FEM). The extended global recovery error estimator (XGR) is applied on the present work to assess the interpolation error with the X-FEM. Different mesh adaptation strategies based on the Equal Distribution Criterion and Minimum Number Criterion are analysed on the classic Westergaard problem.

The two strategies that are most adapted to X-FEM are applied on a 3D crack propagation simulation. We draw the following conclusions:

1. Thanks to the X-FEM, the global convergence rate in energy used in the Min-Num criterion for a problem involving discontinuities and singularities is optimal and equal to unity (1) for linear interpolation of the displacement field. This constant (independent of the crack length) global convergence rate provided by singular enrichment makes the mesh adaptation strategy based on the Min-Num criterion is well-suited to solutions provided through enriched approximations such as X-FEM. In particular, enrichment of the finite element approximation precludes the need for the procedure proposed in [29], which is required when the global convergence rate is a function of the crack size.

considering both global and elementary convergence rates to unify the elementary energynorm error, it is more adapted to X-FEM than the two others. We show evidence that this inadequacy is limited to the element containing the crack tip.

3. The application of the adapted mesh adaptation strategies can considerably reduce the total number of degrees of freedom required to reach a target global error level, typically by one order of magnitude in 3D. In addition, it can lead to higher effectivity of the error estimation approach.

4. The two strategies that are fit for X-FEM are applied on a 3D crack propagation simulation. They provide nearly identical simulation results in terms of crack path, cumulative cycle number and element number. The two simulations with mesh adaptation are compared with a simulation performed on a fixed mesh. The crack front morphologies after 32 propagation steps are very close. The difference in cumulative number of cycles is less than 5% which is a remarkable result, given that the error in stress intensity factors is compounded by large values of the Paris exponent. This similarity illustrates the applicability of the present mesh adaptation procedure.

A drawback of our approach, which needs further scrutiny is the added cost of computing the error and adapting the mesh. In particular due to our use of a global error estimator, this cost is clearly not negligible, especially for 3D cases and when the cracks have propagated significantly, i.e. when the adapted mesh is relatively fine. Besides, it is shown in the present study that the energy norm error of the strain field is not necessarily consistent with the error in a specific quantity of interest, e.g. the stress intensity factors. In other words, an optimal mesh with respect to the energy error, as computed within the current work, will not be optimal when it comes to minimising the error on the crack driving force, the crack propagation length or the crack propagation direction.

To address these two problems, combining the analysed mesh adaptation strategies with local goal-oriented error estimation techniques is a natural perspective of the present study which is being carried out within commercial code Morfeo Crack.

Strategy 2 :

 2 uniform-distribution-simplified Simplified version of the combination of the global error and the element based local one (Eq. 20). Strategy 3: uniform-distribution-global-local Combination of the global error and the volume-based local error (Eq. 19).

Figure 1 :

 1 Figure 1: Mesh adaptation algorithm.

Figure 2 :

 2 Figure 2: The numerical strain field compared with the analytical solution using Strategy 4.(a) sketch of the problem on the initial mesh, the enrichment area is highlighted by the red circle; (b) the numerical strain (yy) using the mesh obtained from the first remesh increment; (c) the difference between the numerical strain and the analytical strain (yy).

Figure 5 :Figure 6 :

 56 Figure 5: The energy norm error evolution weighted by the energy norm (a) and the effectivity index of the error estimation (b) expressed as a function of the mesh adaptation increments.With the first three adaptation strategies, the target error is reached after three steps.

Figure 8 :

 8 Figure8: Convergence analysis for the crack tip element, an enriched element, an element beyond the enrichment area and the overall simulation domain; (a) Estimated and exact energy norm relative errors versus mesh refinement (red and blue triangles describe respectively the convergence rates equal to 1 and 2); (b) Effectivity indices versus mesh refinement. We observe that optimal convergence rates are retrieved. Best effectivity index is obtained for elements outside the enrichment area.

Figure 9 :

 9 Figure9: The mesh adaptation sequences of Strategy 5 ; the total number of elements is indicated below each mesh. Mesh refinement if rapidly followed by mesh optimization iterations.

Figure 10 :

 10 Figure10: Evolution of the ernergy norm error (a) and the effectivity index of the error estimation (b) as a function of the mesh adaptation increments. The target global error level is nearly reached after the first mesh adaptation increment and followed by mesh optimization.

Figure 13 :

 13 Figure 13: The local effectivity of (a) S1 (5); (b) S2 (5); (c) S4 (3); (d) S5

Figure 14 :

 14 Figure 14: (a) The adapted mesh obtained by the proposed hybrid strategy (the total number of elements is indicated below the mesh); (b) the corresponding absolute energy norm error by element. Similar results as for Strategy 1 are observed.

Figure 15 :

 15 Figure 15: Elements for which the size is smaller than 0.01 (larger ones are omitted) in (a) the sixth mesh adaptation iteration of proposed hybrid strategy (b) the sixth mesh adaptation iteration of Strategy 1. Hybrid strategy generates smaller elements at crack tip.

Figure 16 :

 16 Figure 16: The meshes with an enrichment radius of 0.3 obtained after six adaptation increments by (a) Strategy 1 and (b) Strategy 5.The enrichment zone is highlighted by the blue dashed line. We notice that the local refinement still occurs along the transition region between the enriched and standard areas. This is believed to be due to the lack of blending correction used in our computational approach.

Figure 17 :

 17 Figure 17: Initial configuration of the 3D crack propagation case. Geometric description (a) and norm of the displacement field (b); the initial crack shape is a quarter of an ellipse (b)

Figure 18 :

 18 Figure 18: Initial mesh for the 3D mesh adaptation test case.

Figure 20 :Figure 21 : 1793540 Figure 22 :

 2021179354022 Figure 20: Crack propagation after 1, 10, 20, and 32 propagation steps with mesh adaptation applying the Strategy 5 ; the total number of elements is indicated below each corresponding mesh.

  .

Figure 23 :Figure 24 :

 2324 Figure 23: Comparison of Crack tip front morphologies obtained via different mesh adaptation strategies; (a)the crack fronts in the 3D geometry; (b)zoom of the crack fronts. Very similar crack fronts are observed.

  of global and local errors Based on the global and local error parameters, different mesh size adaptation rules were proposed in the literature. By taking into account the different convergence rates of the local and global error, Onãte et al. [48] proposed the following two strategies with element and volume based local errors

	• Element

No No Yes Yes No Crack propagation loop Error control loop Mesh adaptation loop

		Start
		Initial Crack
		Max. propagation steps reached ?	Stop
		Yes
		Current mesh meets
	Interpolate level-sets on new mesh	stopping criteria (error, iterations,...) ?
	Compute crack growth
	rate and direction
		Solve FE BVP problem
	Update level-sets	Get displacement, strain, stress fields
	to propagate crack	(see sec.2.1.)
		Solve XGR problem
		Get smoothed strain field
		(see sec. 3.1.)
		Compute error in energy norm
		(see sec. 3.2.)
		Compute required element size map
		(see sec. )
		Mesh respects
		prescribed
		size map ?
		Edge collapse
		Edge swap
		Remove slivers
		Edge split

Table 1 :

 1 Evaluation of number of elements at the iterations of the error control loop using different adaptation strategies. Contrary to uniform refinement, S1-S3 result in a stagnation of the increase in the number of elements after three adaptation steps.

	trategy	Increment	(1)	(2)	(3)	(4)	(5)	(6)
	S1 -Equal-distribution criterion	2090	5256	13793 15965 15620 15402
	Global/element based local errors
	S2 -Equal-distribution criterion	2778	7206	17893 19254 18640 17837
	Simplified global/element based local errors
	S3 -Equal-distribution criterion	2839	15667 45445 53721 53637 53920
	Global/volume based local errors
	S4 -Uniform mesh refinement	2836	7344	19758 48980 124336

are indicated in Tab.1. For uniform refinement, only five adaptation `````````````S

Table 2 :

 2 Comparison of the five mesh adaptation strategies in terms of energy norm error over energy norm (

	0.01

Figure

11

: The energy norm relative error expressed as a function of the average mesh size; the red triangle describes the theoretical convergence rate, which equals 1 and the blue rectangle highlights the fact that a comparable error level is achieved by the five strategies.

Table 3 :

 3 Comparison of the five used mesh adaptation strategies in terms of energy norm error (weighted by the total energy norm), the error in K I and the effectivity on the meshes leading to a similar error level.

	Strategy Number of elements Number of Dofs Err	Error in K I Effectivity
	S1	20728	24816	0.0116 1.5E-4	0.961
	S4	48980	64934	0.0118 2.6E-4	0.927
	S5	20052	23876	0.0119 1.5E-4	0.96

note that m is typically larger than unity and can reach values above 5. In such cases, any error on the stress intensity factor is proportionally compounded, leading to significant errors on the crack growth increment

Note that an energy-error-optimal mesh may not be an SIF-optimal mesh which explains the significant work which has been expended in devising goal-oriented (SIF-oriented) error estimates, both of residual type[59, 60] and recovery-based [61]

the mesh used is that obtained via Strategy 1 and given for reference in Fig.3

Methods based on strain smoothing were developed, also for enriched approximations, which allow to somewhat control the deleterious effects of mesh distortion[57]

η 1 , which controls the mesh size evolution ratio between two iterations, still equals 0.6.

Note that we do not perform any blending correction in those simulations, which is the likely cause of this refinement at the interface between the tip-enriched and standard regions.

Note that in order to further decrease the error on the stress intensity factors, goal-oriented error estimators would be a sensible approach. Those have been developed for fracture mechanics treated by an X-FEM approach and are being implemented within Morfeo Crack to improve the results obtained in this paper.

Three different mesh adaptation strategies based on the Equal-Dist criterion are analysed on the Westergaard problem. Despite a local inadequacy of the adaptation strategy

Acknowledgements

This work was partially funded by the EPSRC grant EP/G042705/1 "Increased Reliability for Industrially Relevant Automatic Crack Growth Simulation with the eXtended Finite Element Method" and an SBO Project grant 110070 "Effective Structural Health Monitoring with Additive Manufacturing" of the Agency for Innovation by Science and Technology (IWT).

Stéphane Bordas also thanks partial funding for his time provided by the European Research Council Starting Independent Research Grant (ERC Stg grant agreement No. 279578) "RealCut Towards real time multiscale simulation of cut-ting in non-linear materials with applications to surgical simulation and computer guided surgery".

The authors would like to thank A. François, N. Poletz, E. Wyart and K. Agathos for many insightful discussions about the implementation of the numerical framework.

Convergence rate analysis

For the classical finite element method, theoretically, the convergence rate of the global energy