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

C. Roux Langlois, A. Gravouil, Marie-Christine Baietto, Julien Réthoré

► To cite this version:

C. Roux Langlois, A. Gravouil, Marie-Christine Baietto, Julien Réthoré. Three-dimensional simulation of crack with curved front with direct estimation of stress intensity factors. International Journal for Numerical Methods in Engineering, 2015, 101 (8), pp.635-652. 10.1002/nme.4811. hal-01951564

HAL Id: hal-01951564

https://hal.science/hal-01951564

Submitted on 10 May 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Three-dimensional simulation of crack with curved front with
direct estimation of stress intensity factors

C. Roux Langlois¹, A. Gravouil¹,²,*,†, M.-C. Baietto¹ and J. Réthoré¹

¹Université de Lyon, INSA-Lyon, LaMCoS, CNRS UMR 5259, F-69621 Villeurbanne, (France)

²Institut Universitaire de France, Paris, France

SUMMARY

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.

KEY WORDS: 3D fracture mechanics; mixed-mode; Williams series expansion; X-FEM; multigrid

1. 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 sim-
pler 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–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 prop-
agation 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 Muskelishvili and Williams [13, 14]. Some
other terms of this series expansion have been proven meaningful, such as the T-stress [15, 16] or some
super-singular terms [17, 18].

*Correspondence to: A. Gravouil, LaMCoS (Laboratoire de Mécanique des Contacts et des Structures), Bat. J. d'Alembert, 18 rue des Sciences, F-69621 Villeurbanne Cedex, France

†E-mail: anthony.gravouil@insa-lyon.fr
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–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].

## 2. ASYMPTOTIC FIELDS IN THREE-DIMENSIONS

### 2.1. 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.

#### 2.1.1. The Williams series expansion [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 ($\sigma \cdot n = 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

$$ \sigma(r, \theta) = \sum_{i=1, II} \sum_{n=0}^{\infty} b_i^n r^{n/2-1} f_i^n(\theta) \quad \text{and} \quad \mathbf{u}(r, \theta) = \sum_{i=1, II} \sum_{n=0}^{\infty} b_i^n r^{n/2} \mathbf{g}_i^n(\theta), \quad (1) $$
<figure><img src="image_2.png" /><figcaption>Figure 1. Local coordinate system around the crack tip. (a) Two-dimensions. (b) Three-dimensions.</figcaption></figure>

$$
\begin{align}
f_I^n(\theta) &= \begin{bmatrix} (f_I^n)_{tt} \\ (f_I^n)_{nn} \\ (f_I^n)_{tn} \end{bmatrix} = \frac{n}{2} \left[ \begin{bmatrix} (2 + \frac{n}{2} + (-1)^n) \cos\left[\left(\frac{n}{2} - 1\right)\theta\right] - \left(\frac{n}{2} - 1\right) \cos\left[\left(\frac{n}{2} - 3\right)\theta\right] \\ (2 - \frac{n}{2} - (-1)^n) \cos\left[\left(\frac{n}{2} - 1\right)\theta\right] + \left(\frac{n}{2} - 1\right) \cos\left[\left(\frac{n}{2} - 3\right)\theta\right] \\ \left(\frac{n}{2} - 1\right) \sin\left[\left(\frac{n}{2} - 3\right)\theta\right] - \left(\frac{n}{2} + (-1)^n\right) \sin\left[\left(\frac{n}{2} - 1\right)\theta\right] \end{bmatrix}_{(e_t, e_n)}, \\
f_{II}^n(\theta) &= \begin{bmatrix} (f_{II}^n)_{tt} \\ (f_{II}^n)_{nn} \\ (f_{II}^n)_{tn} \end{bmatrix} = \frac{n}{2} \left[ \begin{bmatrix} -(2 + \frac{n}{2} - (-1)^n) \sin\left[\left(\frac{n}{2} - 1\right)\theta\right] + \left(\frac{n}{2} - 1\right) \sin\left[\left(\frac{n}{2} - 3\right)\theta\right] \\ -(2 - \frac{n}{2} + (-1)^n) \sin\left[\left(\frac{n}{2} - 1\right)\theta\right] - \left(\frac{n}{2} - 1\right) \sin\left[\left(\frac{n}{2} - 3\right)\theta\right] \\ \left(\frac{n}{2} - 1\right) \cos\left[\left(\frac{n}{2} - 3\right)\theta\right] - \left(\frac{n}{2} - (-1)^n\right) \cos\left[\left(\frac{n}{2} - 1\right)\theta\right] \end{bmatrix}_{(e_t, e_n)}, \tag{2} \\
g_I^n(\theta) &= \frac{1}{2\mu} \left[ \begin{bmatrix} (\kappa + n/2 + (-1)^n) \cos\left[(n/2)\theta\right] - (n/2) \cos\left[(n/2 - 2)\theta\right] \\ (\kappa - n/2 - (-1)^n) \sin\left[(n/2)\theta\right] + (n/2) \sin\left[(n/2 - 2)\theta\right] \end{bmatrix}_{(e_t, e_n)}, \right. \\
g_{II}^n(\theta) &= \frac{1}{2\mu} \left[ \begin{bmatrix} -(\kappa + n/2 - (-1)^n) \sin\left[(n/2)\theta\right] + (n/2) \sin\left[(n/2 - 2)\theta\right] \\ (\kappa - n/2 - (-1)^n) \cos\left[(n/2)\theta\right] + (n/2) \cos\left[(n/2 - 2)\theta\right] \end{bmatrix}_{(e_t, e_n)}, \right]
\end{align}
$$

in the crack coordinate system defined Figure 1. $\mu$ is the shear modulus and $\kappa$ is the Kolossov's constant (function of Poisson's ratio $\nu$), which takes the value $\kappa = 3 - 4\nu$ in a plane strain state and $\kappa = (3 - \nu)/(1 + \nu)$ for plane stress conditions. The asymptotic coefficients $b_i^n$ define the stress and displacement fields and mainly depend on non-local effects (boundary conditions, etc.). $b_I^0$ and $b_{II}^0$ are in-plane rigid body translations, and $b_{II}^2$ is the in-plane rigid body rotation. $b_I^1$ and $b_{II}^1$ are proportional to the SIF, and $b_I^2$ is related to the T-stress, with

$$b_I^1 = \frac{K_I}{\sqrt{2\pi}} \quad \text{and} \quad b_{II}^1 = \frac{K_{II}}{\sqrt{2\pi}} \quad \text{and} \quad b_I^2 = \frac{T}{4}. \qquad (3)$$

In the crack tip vicinity, all the terms $f_I^n$ are vanishing except the SIF (singular stress), the T-stress (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 ($\propto r^{n/2}$).

2.1.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 = u(r, \theta)e_3$, where $e_3$ is the out-of-plane 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) by extending the summation index to $i = I, II, III$

$$
\begin{align}
f_{III}^{n}(\theta) &= \begin{bmatrix} (f_{III}^{n})_{3n} \\ (f_{III}^{n})_{3t} \end{bmatrix} = \begin{bmatrix} -n \sin\left(\frac{n}{2}(\theta - \pi) - \theta\right) \\ n \cos\left(\frac{n}{2}(\theta - \pi) - \theta\right) \end{bmatrix}_{(e_n, e_t, e_3)}, \tag{4} \\
g_{III}^{n}(r, \theta) &= \frac{2}{\mu} \cos\left(\frac{n}{2}(\theta - \pi)\right) e_3.
\end{align}
$$

Two out-of-plane rigid body motions are added—$b_{III}^0$, the out-of-plane rigid body translation, and $b_{III}^2$, the out-of-plane rigid body rotation. The coefficient $b_{III}^1$ is proportional to the mode III stress intensity factor.
## 2.2. 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

$$ \mathbf{u}(r, \theta, s) = \sum_{i=I, II, III} \sum_{n=0}^{\infty} b_i^n r^{n/2} \mathbf{h}_i^n(\theta, s), \quad (5) $$

where s is the position along the front.

Furthermore, the missing third rigid body rotation is added. From this statement, the evolution of $\mathbf{h}_i^n(\theta, 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_i^n$ 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_\theta, e_s)$ is following the crack and the front, which are curved in the general case. This basis is illustrated in Figure 1(b).

The displacement in Equation (5) is a function of local parameters $(n, \theta, s)$; the computation of $\mathbf{D}$ (its derivative), with respect to these parameters, is therefore straightforward. $\mathbf{D}$ can be linked to the displacement gradient, thanks to the Jacobian matrix $\mathbf{J}$. This gradient is needed in the global basis $(e_1, e_2, e_3)$, in which the displacement is decomposed as $\mathbf{u} = V_1e_1 + V_2e_2 + V_3e_3$

$$ \nabla \mathbf{u} = \begin{bmatrix} \frac{\partial V_1}{\partial x_1} & \frac{\partial V_1}{\partial x_2} & \frac{\partial V_1}{\partial x_3} \\ \frac{\partial V_2}{\partial x_1} & \frac{\partial V_2}{\partial x_2} & \frac{\partial V_2}{\partial x_3} \\ \frac{\partial V_3}{\partial x_1} & \frac{\partial V_3}{\partial x_2} & \frac{\partial V_3}{\partial x_3} \end{bmatrix}_{(e_1, e_2, e_3)} $$

$$ \nabla \mathbf{u} = \mathbf{D} \cdot \mathbf{J} \quad \text{with} \quad \mathbf{D} = \begin{bmatrix} \frac{\partial V_1}{\partial r} & \frac{\partial V_1}{r \partial \theta} & \frac{\partial V_1}{r \partial s} \\ \frac{\partial V_2}{\partial r} & \frac{\partial V_2}{r \partial \theta} & \frac{\partial V_2}{r \partial s} \\ \frac{\partial V_3}{\partial r} & \frac{\partial V_3}{r \partial \theta} & \frac{\partial V_3}{r \partial s} \end{bmatrix}_{(e_1, e_2, e_3)} \quad \text{and} \quad \mathbf{J} = \begin{bmatrix} \frac{\partial r}{\partial x_1} & \frac{\partial r}{\partial x_2} & \frac{\partial r}{\partial x_3} \\ \frac{r \partial \theta}{\partial x_1} & \frac{r \partial \theta}{\partial x_2} & \frac{r \partial \theta}{\partial x_3} \\ \frac{\partial s}{\partial x_1} & \frac{\partial s}{\partial x_2} & \frac{\partial s}{\partial x_3} \end{bmatrix} \quad (6) $$

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.

# 3. FORMULATION

## 3.1. Mono-scale formulation

The considered problem is a material body in $\Omega$, a domain of the three-dimensional space, with well-posed boundary conditions: prescribed displacements $\mathbf{u}_d$ on $\partial_1\Omega$ and forces $\mathbf{F}_d$ on $\partial_2\Omega$,
with $\partial\Omega = \partial_1\Omega \cup \partial_2\Omega$ and $\partial_1\Omega \cap \partial_2\Omega = \emptyset$. 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 $\mathbf{u}$ alone. This displacement field belongs to the admissible displacement space $\mathcal{U} = \{\mathbf{u} \in H^1(\Omega) | \mathbf{u} = \mathbf{u}_d \text{ on } \partial_1\Omega\}$. The virtual work equation without body forces, equivalent to the governing equations for the cracked body as shown in [37], is given by

$$ \int_{\Omega} \mathbf{H} \nabla^s \mathbf{u} : \nabla^s \mathbf{v} d\Omega = \int_{\partial_2\Omega} \mathbf{v} \cdot \mathbf{F}_d dS \quad \forall \mathbf{v} \in \mathcal{U}^0 \quad (7) $$

where $\nabla^s \mathbf{u}$ is the symmetric part of the gradient of $\mathbf{u}$, corresponding to the strain. $\mathbf{H}$ is the Hooke tensor. $\mathbf{v}$ is a test function in $\mathcal{U}^0$, the vector space corresponding to $\mathcal{U}$ with $\mathbf{u}_d = 0$.

### 3.1.1. 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 $\Omega_W$, where boundary conditions are far and the crack almost plane. The subscript $W$ will be used for this domain. The complementary domain $\Omega_X$, such as $\Omega = \Omega_X \cup \Omega_W$ and $\emptyset = \Omega_X \cap \Omega_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 $\Gamma_W$, between $\Omega_X$ and $\Omega_W$, is ensured by mortar approach as proposed in [29]. Such a matching is necessary because the considered displacements in $\Omega_X$ and $\Omega_W$ are not compatible. This condition is enforced in a weak sense with Lagrange multipliers $\lambda \in \mathcal{L}$:

$$ \int_{\Gamma_W} \lambda (u_W - u_X) = 0 \quad \forall \lambda \in \mathcal{L} \quad (8) $$

The virtual work equation without body forces for the coupled two domains and the displacement matching can be written as follows, seeking $(u_W, u_X, \lambda) \in \mathcal{U}_W \times \mathcal{U}_X \times \mathcal{L}$:

$$ \begin{aligned} \int_{\Omega_W} \mathbf{H} \nabla^s \mathbf{u}_W : \nabla^s \mathbf{v}_W d\Omega &+ \int_{\Omega_X} \mathbf{H} \nabla^s \mathbf{u}_X : \nabla^s \mathbf{v}_X d\Omega + \int_{\Gamma_W} \mu (u_W - u_X) dS + \int_{\Gamma_W} \lambda (v_W - v_X) dS \\ = & \int_{\partial_2 \Omega_W} \mathbf{v}_W \cdot \mathbf{F}_d dS + \int_{\partial_2 \Omega_X} \mathbf{v}_X \cdot \mathbf{F}_d dS \quad \forall (v_W, v_X, \mu) \in \mathcal{U}_W^0 \times \mathcal{U}_X^0 \times \mathcal{L}. \end{aligned} \quad (9) $$

$\mathcal{U}_W^0$ refers to the admissible displacement space on $\Omega_W$ where imposed displacement is zero $\mathcal{U}_W^0 = \{\mathbf{u} \in H^1(\Omega_W) | \mathbf{u} = \mathbf{0} \text{ on } \partial_1\Omega_W\}$. However, the asymptotic series expansion would not be able to handle directly most of the Dirichlet boundary condition; thus, situations where $\partial_1\Omega_W = \emptyset$ 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 $\Omega_X$ and the Williams series expansion on the domain $\Omega_W$.

### 3.1.2. X-FEM domain $\Omega_X$.
A classical three-dimensional extended finite element method, as introduced in [38], is used in the domain $\Omega_X$. The domain is meshed, and a finite element interpolation is set $\mathbf{u}_X(\mathbf{x}) = \sum_{i \in N} u_i \varphi_i(\mathbf{x})\varphi_i(\mathbf{x})$ are the finite element shape functions, and $u_i$ the degrees of freedom. The partition of unity property of this method is used to enrich the mechanical fields
around the crack. The elements cut by the crack have their degrees of freedom $N_c$ enriched with the discontinuity generated by the crack

$$ \mathbf{u}_X(\mathbf{x}) = \sum_{i \in N} u_i \varphi_i(\mathbf{x}) + \sum_{i \in N_c} u'_i \varphi_i(\mathbf{x}) H(\mathbf{x}). \quad (10) $$

The function $H(\mathbf{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 $\Omega_X$. Therefore, no singular enrichment is needed.

Let $U_X$ be the degrees of freedom vector containing $u_i$ and $u'_i$. Equation (10) can be rewritten as a matrix $\mathbf{u}_X = N_X^T(\mathbf{x})U_X$. The stiffness matrix $\mathbf{K}_X$ is then

$$ \mathbf{K}_X = \int_{\Omega_X} H \nabla^s N_X^T(\mathbf{x}) : \nabla^s N_X^T(\mathbf{x}) d\Omega. \quad (11) $$

### 3.1.3. Analytical patch.
The resolution of the considered problem on $\Omega_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 $\Omega_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_i^n$ of Equation (5) is discretized with one-dimensional finite elements

$$ \mathbf{u}_W(\mathbf{x}) = \mathbf{u}_W(r, \theta, s) = \sum_{k \in N_s} \left[ \sum_{n=0}^{n_{max}} \sum_{i=1,II,III} b_{ink} \mathbf{g}_i^n(\theta) r^{\frac{n}{2}} \right] \varphi_k(s) \quad (12) $$

$$ = \sum_{j \in \{1, \dots, N_s \cdot (n_{max}+1) \cdot 3\}} b_j \psi_j = N_W^T(\mathbf{x}) U_W, \quad (13) $$

where $\varphi_k(s)$ are one-dimensional shape functions of order one, as represented Figure 2(b). $b_{ink}$ are the $b_i^n$ coefficients at each of the $k \in N_s$ nodes. The curvilinear abscissa $s$ on the front is computed in the whole patch $\Omega_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 = s_{k+1} - s_k$, $k \in N_s$.

From this displacement field, the symmetric part of the gradient field is computed in order to compute the stiffness matrix $\mathbf{K}_W$ such as

$$ \mathbf{K}_W = \int_{\Omega_W} \mathrm{Tr} \left[ \mathbf{H} \nabla^s N_W^T : \nabla^s N_W^T \right] d\Omega. \quad (14) $$

Since the three variables $r$, $\theta$, and $s$ are uncoupled, the computation of the matrix $\mathbf{D}$ is easily performed.

<figure><img src="image_3.png" /><figcaption>Figure 2. Example of the numerical implementation of the patch, in the three-dimensional case. (a) Mesh of the patch, in a plan containing the front. (b) One-dimensional finite element shape functions along $s$.</figcaption></figure>
### 3.1.4. 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$

$$ \lambda = \sum_l \lambda_l L_l = N_L^T \Lambda \qquad (15) $$

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 $\Gamma_W: N_L(x) = N_W(x), \forall(x) \in \Gamma_W$. The discrete version of Equation (8) can be rewritten in the matrix form

$$ C_X U_X + C_W U_W = 0 \quad \text{with} \quad C_X = \int_{\Gamma_W} N_X^T N_L dS \quad \text{and} \quad C_W = - \int_{\Gamma_W} N_W^T N_L dS. \quad (16) $$

### 3.1.5. Practical implementation.
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 $\Omega_W$. And the remaining elements define the X-FEM domain $\Omega_X$ and are directly used as X-FEM elements. The interface $\Gamma_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 $\Gamma_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 $\Lambda$ such that

$$ \begin{bmatrix} K_X & 0 & -C_X^T \\ 0 & K_W & -C_W^T \\ -C_X & -C_W & 0 \end{bmatrix} \begin{pmatrix} U_X \\ U_W \\ \Lambda \end{pmatrix} = \begin{pmatrix} F \\ 0 \\ 0 \end{pmatrix}, \qquad (17) $$

where $F$ is the generalized force vector. However, in three-dimensions, $K_W$ is bigger than in two-dimensions 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 = 2$ finite elements, and a truncation of the Williams series expansion to $n_{max} = 7$, was sufficient for an accurate estimate of the SIFs [29]. The same assumption will be used in the three-dimensional case.

## 3.2. 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
<figure><img src="image_4.png" /><figcaption>Figure 3. Recursive grids example for a square cracked problem with four grids ($n_g = 4$).</figcaption></figure>

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 $\sqrt{r}$. The approximation space is then

$$ \mathbf{u}_X(\mathbf{x}) = \sum_{i \in N} u_i \varphi_i(\mathbf{x}) + \sum_{i \in N_c} u'_i \varphi_i(\mathbf{x}) H(\mathbf{x}) + \sum_{i \in N_e} \sum_{j} u''_{ij} \varphi_i(\mathbf{x}) \varphi_j(\mathbf{x}). \quad (18) $$

The functions $\varphi_j(\mathbf{x})$ are enrichment functions designed to introduce the singular behavior (term $n = 1$ in Equation (5)) in the approximation space $\mathcal{U}_X$ [5]

$$ \varphi_j(\mathbf{x}) \in \left\{ \sqrt{r} \sin \frac{\theta}{2}, \sqrt{r} \cos \frac{\theta}{2}, \sqrt{r} \sin \frac{\theta}{2} \sin \theta, \sqrt{r} \cos \frac{\theta}{2} \sin \theta \right\}. \quad (19) $$

### 3.2.1. 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_{g+1}$ (respectively, the coarse and the fine). Since the refinement is local, a boundary $\Gamma_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 $|B$ subscript of the quantity $Q$ refers to the restriction to a zone where a finer grid exists; $Q|_B$ 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 $v_1$ conjugate gradient (CG) iterations called relaxations are performed and (2) $v_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 $\gamma_g = 1 \quad \forall g \in \{n_g - 1, \dots, 2\}$.

0. Initiation: the displacements are set to 0, $g = n_g, c_g = 1 \quad \forall g \in \{n_g - 1, \dots, 2\}$ and $R_{n_g}^i = 0$.

(1) Corrective residual term from fine grids, $g = g - 1$

$$ R_g^i = P_g^T [(K_{g+1} - K_{g+1|B}) U_{g+1}^i + R_{g+1}^i], \quad (20) $$
$$ v_1 \text{ relaxations of } K_g U_g^{i+1} = F_g - R_g^i + K_{g|B} U_g^i \tag{21} $$
with boundary condition $U_g^i|_{\Gamma_{g-1}} = P_g U_{g-1}^i |_{\Gamma_{g-1}}.$

(a) If $g = 1$: go to 2.

(b) If $c_g = \gamma_g$: $c_g = 1$ and return to step 1

Else: $c_g = c_g + 1$ and $g = g + 1$.

(2) Relaxation of fine grid correction on coarse grids

(a) If $g = 1$:

(i) $U_1^{i+1}$ from exact resolution of $K_1 U_1^{i+1} = F_1 - R_1^i + K_{1|B} U_1^i,$

(ii) **Convergence test:** if $\max_g \frac{\|U_g^{i+1} - U_g^i\|_E}{\|U_1^{i+1}\|_E} < \epsilon$ stop, $\|U_g\|_E = \frac{1}{2} U_g^T K_g U_g$.

Else:

(i) Interpolate displacements from coarse to fine mesh $U_g^{i+1/2} = U_g^i + P_g (U_{g-1}^{i+1} - U_{g-1}^i),$

(ii) $U_g^{i+1}$ obtained from $v_2$ CG iterations on $U_g^{i+1/2}$

$$ K_g U_g^{i+1} = F_g - R_g^i + K_{g|B} U_g^i \text{ with boundary } U_g^i|_{\Gamma_{g-1}} = P_g U_{g-1}^i |_{\Gamma_{g-1}}. \tag{22} $$

(b) If $g = n_g$: Return to step 1

Else: Return to 2 with $g = g + 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.

### 3.2.2. DEK-FEM coupled with the localized X-FEM multigrid method.
The analytical patch is inserted in the finest grid $g = 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|W} U_{n_g}^i$, where $K_{n_g|W}$ 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

$$ K_{n_g} U_{n_g}^{i+1} = F_{n_g} \rightarrow \begin{cases} \text{X-FEM model: } K_{n_g} U_{n_g}^{i+1} = F_{n_g} + C_X^T \Lambda^i + K_{n_g|W} U_{n_g}^i, \\ \text{Patch resolution: } \begin{bmatrix} K_W & -C_W^T \\ -C_W & \mathbf{0} \end{bmatrix} \begin{pmatrix} U_W^{i+1} \\ \Lambda^{i+1} \end{pmatrix} = \begin{bmatrix} \mathbf{0} \\ C_X U_{n_g}^{i+1} \end{bmatrix}. \end{cases} \tag{23} $$

The analytical patch influences the X-FEM model through $\Lambda^i$. In [30], the authors noticed that $\Lambda^{i+1}$ can be expressed as a function of the coarse problem unknowns $U_{n_g}$ only

$$ \underbrace{C_W K_W^{-1} C_W^T}_{H_W} \Lambda^{i+1} = -C_X U_{n_g}^{i+1} \rightarrow \Lambda^{i+1} = -H_W^{-1} C_X U_{n_g}^{i+1}, \tag{24} $$

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

$$ \left[ K_{n_g} - K_{n_g|W} + \underbrace{C_X^T H_W^{-1} C_X}_{H_{XW}} \right] U_{n_g}^{i+1} = F_{n_g}. \tag{25} $$
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} = C_X^T H_W^{-1} C_X$ is performed in two steps:

(1) Find $Y$ such that $K_W Y = C_W^T$; thus, $H_W = C_W Y$.

(2) Find $Z$ such that $H_W Z = C_X^T$; thus, $H_{XW} = C_X^T Z$.

## 4. NUMERICAL SIMULATIONS

Some classical cracked configurations are studied here, and the identified stress intensity factors ($K_i^{DEK}$, $i = I, II, III$) and T-stress are compared to former analytical and numerical solutions ($K_i^{ref}$). 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^\circ$ 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} = 7$, $r_W = 2$, and $ds \simeq 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).

### 4.1. 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 ($\sigma_0$) in an infinite domain is subject to

$$K_I = 2\sigma_0 \sqrt{\frac{a}{\pi}} \quad \text{and} \quad T = -\frac{1+2\nu}{2} \sigma_0, \qquad (26)$$

where $\nu$ is the Poisson's ratio.

For the simulation, a cube of size $2l$ submitted to uniform tensile stress $\sigma_0$ is considered with an embedded penny-shaped crack of radius $a = 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 = a/3.5$. (2) And a fine one, twice as fine. For this mesh, the crack radius is about seven finite elements: $h = 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 two-dimensions 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 $\Gamma_W$
<figure><img src="image_5.png" /><figcaption>Figure 4. Faces of the three-dimensional graded meshes used for penny shaped crack under remote tension. The considered domain is cubic, and the mesh is the same in the three directions. (a) A face of the coarse mesh. (b) A face of the fine mesh.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Figure 5. Penny shaped crack under remote tension. Numerical results of the comparison between a coarse mesh and a fine mesh, $ds \simeq 2$ finite elements of the coarse mesh, and four finite elements of the fine mesh. (a) Evaluated mode I SIF. (b) Evaluated T-stress.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Figure 6. Penny shaped crack under remote tension. Numerical results of the study of the influence of $r_w$ for a fine graded mesh. (a) Evaluated mode I SIF. (b) Evaluated T-stress.</figcaption></figure>

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 \approx 2$ finite elements.
<figure><img src="image_8.png" /><figcaption>Figure 7. Penny shaped crack under remote shear. (a) Geometry and boundary conditions. (b) Numerical results: evaluated SIF.</figcaption></figure>

## 4.2. Penny shaped crack under remote shear

A penny shaped crack under shear loading is considered. As illustrated in Figure 7(a), tangential displacement is enforced on the borders parallel to the crack plane ($\tau_0 = \frac{E}{2+(1+\nu)} \arctan \frac{2u_d}{l}$). 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

$$K_{II}^{\text{ref}} = \tau_0 \frac{4(1-\nu)\sqrt{\pi a}}{\pi(2-\nu)} \cos \alpha \quad \text{and} \quad K_{III}^{\text{ref}} = \tau_0 \frac{4\sqrt{\pi a}}{\pi(2-\nu)} \sin \alpha. \qquad (27)$$

Where $\alpha$ is the angle along the front. Here, considering a prism of height $l$ and square section of edge $2l$ with $l = a$ is sufficient in order to compare the simulation to the analytical solution. A structured mesh of $40 \times 40 \times 20$ finite elements is used.

The identified symmetric coefficients (mode I) are almost zeros,

$$\left| \frac{K_I}{\tau_0 \sqrt{a}} \right| < 0.3\% \quad \text{and} \quad \left| \frac{T}{\tau_0} \right| < 1.3\%. \qquad (28)$$

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.

## 4.3. Inclined penny shaped crack under remote tension

Let us consider the geometry and loading of Section 4.1 with an inclined crack (of angle $\gamma = 45^\circ$ 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

$$\begin{align}
K_I^{\text{ref}} &= 2\sigma_0 \sin^2 \gamma \sqrt{\frac{a}{\pi}}, \notag \\
K_{II}^{\text{ref}} &= -2\sigma_0 \sin \gamma \cos \gamma \sqrt{\frac{a}{\pi}} \cos \alpha, \tag{29} \\
K_{III}^{\text{ref}} &= 2\sigma_0 \sin \gamma \cos \gamma (1-\nu) \sqrt{\frac{a}{\pi}} \sin \alpha. \notag
\end{align}$$

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.
<figure><img src="image_9.png" /><figcaption>Figure 8. Geometry and boundary conditions of the inclined (45°) penny shaped crack under remote tension.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Figure 9. Inclined (45°) penny shaped crack under remote tension. Numerical results for a graded mesh (T-stress is normalized by σ). An order one finite elements discretization is used along the front for the three modes. (a) Evaluated asymptotic coefficients for coarse mesh. (b) Evaluated asymptotic coefficients for fine mesh.</figcaption></figure>

The evolution of the identified stress intensity factors are represented in Figure 9. The evolution of $K_{II}^{DEK}$ and $K_{III}^{DEK}$ 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_{II}^{DEK}$ 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., $\varphi_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

$$u_W(r, \theta, s) = \sum_{k \in N_s} \left[ \sum_{n=0}^{n_{max}} \sum_{i=I,II} b_{ink} g_i^n(\theta) r^{\frac{n}{2}} \right] \varphi_k(s) + \sum_{k \in N_s} \sum_{n=0}^{n_{max}} b_{IIInk} g_{III}^n(\theta) r^{\frac{n}{2}} \varphi'_k(s). \quad (30)$$

For the sake of simple numerical implementation, $\varphi'_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_{II}^{DEK}$ and ≈ 15% for $K_{III}^{DEK}$. On the other hand, $K_I$ is accurate < 4%, and the T-stress is oscillating as $\cos(2\alpha)$. When the mesh is coarser, the 0 order discretization of mode III leads to uneven $K_{III}^{DEK}$, which follows the analytical tendency. For the three modes, the discretization enhances the accuracy. This uncoupling method
<figure><img src="image_11.png" /><figcaption>Figure 10. Inclined (45°) penny shaped crack under remote tension. Numerical results for a graded mesh (T-stress is normalized by σ). The direct extraction is based on Equation (30) where φs(s) are order one finite elements shape functions and φs'(s) are order zero. (a) Evaluated asymptotic coefficients for coarse mesh. (b) Evaluated asymptotic coefficients for fine mesh.</figcaption></figure>

(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.

### 4.4. Corner edge cracked bar

A square bar containing a corner crack loaded with remote tension σ₀ 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 = 2w. It includes a quarter penny shaped corner crack of dimension a = w/10. This ratio is chosen in agreement with available closed-form solutions, reminded in [42] for such a corner edge crack

$$K_I^{\text{ref}} = 2\sigma_0 \sqrt{\frac{a}{\pi}} (1.211 - 0.186\sqrt{\sin\alpha}) (1.211 - 0.186\sqrt{\cos\alpha}). \quad (31)$$

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(b) is considered, with 10 finite elements in the

<figure><img src="image_12.png" /><figcaption>Figure 11. Corner edge cracked bar. (a) Geometry and boundary condition. (b) A face (cut by the crack) of the considered mesh.</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>Figure 12. Corner edge cracked bar. Numerical result of the DEK-FEM simulation. (a) Evaluated mode I SIF. (b) T-stress, reference ≈ geometry.</figcaption></figure>

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_I^{DEK}$ 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 $\frac{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 the increase of the *T*-stress are expected. The decrease of $K_I^{DEK}$ 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.

### 4.5. 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 = a$, the thickness $t = 3a$, and the height $h = 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³ and two multigrid ones with two and three grids from a coarse grid of 12 × 8 × 21 elements³. 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.
<figure><img src="image_14.png" /><figcaption>Figure 13. Single edge-crack tension specimen, with the monogrid mesh used in Sukumar *et al.* [38]. (a) Single edge-crack tension specimen. (b) Evaluated SIF for a structured mesh 20³.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>Figure 14. Single edge-crack tension specimen, with the localized multigrid method. (a) Multigrid mesh. (b) Comparison of the $K_I$.</figcaption></figure>

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.

## 5. 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.
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.

## ACKNOWLEDGEMENTS

This work is supported by the French National Research Agency (grant ANR-09-BLAN-0009-01-RUPXCUBE).

## REFERENCES

1. Tong P, Pian THH. On the convergence of the finite element method for problems with singularity. *International Journal of Solids and Structures* 1973; **9**(3):313-321.

2. Tong P, Pian THH, Lasry SJ. A hybrid-element approach to crack problems in plane elasticity. *International Journal for Numerical Methods in Engineering* 1973; **7**(3):297-308.

3. Barsoum RS. On the use of isoparametric finite elements in linear fracture mechanics. *International Journal for Numerical Methods in Engineering* 1976; **10**(1):25-37.

4. Chiaruttini V, Riolo V, Feyel F. Advanced remeshing techniques for complex 3D crack propagation. *ICF13*, Beijing, China, 2013; 547-555.

5. Belytschko T, Black T. Elastic crack growth in finite elements with minimal remeshing. *International Journal for Numerical Methods in Engineering* 1999; **45**(5):601-620.

6. Moës N, Gravouil A, Belytschko T. Non-planar 3D crack growth by the extended finite element and level sets—Part I: Mechanical model. *International Journal for Numerical Methods in Engineering* 2002; **53**(11):2549-2568.

7. Gravouil A, Moës N, Belytschko T. Non-planar 3D crack growth by the extended finite element and level sets—part II: level set update. *International Journal for Numerical Methods in Engineering* 2002; **53**(11):2569-2586.

8. Fries TP, Belytschko T. The extended/generalized finite element method: an overview of the method and its applications. *International Journal for Numerical Methods in Engineering* 2010; **84**(3):253-304.

9. Babuška I, Melenk J. The partition of unity method. *International Journal for Numerical Methods in Engineering* 1997; **40**:727-758.

10. Strouboulis T, Babuška I, Copps K. The design and analysis of the generalized finite element method. *Computer Methods in Applied Mechanics and Engineering* 2000; **181**(1-3):43-69.

11. Laborde P, Pommier J, Renard Y, Salaün M. High-order extended finite element method for cracked domains. *International Journal for Numerical Methods in Engineering* July 2005; **64**:354-381.

12. Suresh S. *Fatigue of Materials*, Cambridge Solid State Science Series. Cambridge University Press: Cambridge, 1998.

13. Muskhelishvili NI. *Some Basic Problems of the Mathematical Theory of Elasticity*. P. Noordhoff: Leyden, 1953.

14. Williams M. On the stress distribution at the base of a stationary crack. *ASME Journal of Applied Mechanics* 1957; **24**:109-114.

15. Larsson SG, Carlsson AJ. Influence of non-singular stress terms and specimen geometry on small-scale yielding at crack tips in elastic-plastic materials. *Journal of the Mechanics and Physics of Solids* 1973; **21**(4):263-277.

16. Rice JR. Limitations to small scale yieldg approximation for crack tip plasticity. *Journal of the Mechanics and Physics of Solids* 1974; **22**:17-26.

17. Hui CY, Ruina A. Why K? High order singularities and small scale yielding. *International Journal of Fracture* 1995; **72**(2):97-120.
18. Roux S, Réthoré J, Hild F. Digital image correlation and fracture: an advanced technique for estimating stress intensity factors of 2D and 3D cracks. *Journal of Physics D: Applied Physics* 2009; **42**(214004):21 pages.

19. Moran B, Shih C. A general treatment of crack tip contour integrals. *International Journal of Fracture* 1987; **35**(4): 295–310.

20. Gosz M, Moran B. An interaction energy integral method for computation of mixed-mode stress intensity factors along non-planar crack fronts in three dimensions. *Engineering Fracture Mechanics* 2002; **69**(3):299–319.

21. González-Albuixech VF, Giner E, Tarancón JE, Fuenmayor FJ, Gravouil A. Domain integral formulation for 3-D curved and non-planar cracks with the extended finite element method. *Computer Methods in Applied Mechanics and Engineering* 2013; **264**(0):129–144.

22. Nakamura T, Parks DM. Determination of elastic t-stress along three-dimensional crack fronts using an interaction integral. *International Journal of Solids and Structures* 1992; **29**(13):1597–1611.

23. Xiao QZ, Karihaloo BL. Direct evaluation of accurate coefficients of the linear elastic crack tip asymptotic field. *Fatigue & Fracture of Engineering Materials & Structures* 2003; **26**(8):719–729.

24. Liu XY, Xiao QZ, Karihaloo BL. XFEM for direct evaluation of mixed mode sifs in homogeneous and bi-materials. *International Journal for Numerical Methods in Engineering* 2004; **59**(8):1103–1118.

25. Karihaloo BL, Xiao QZ. Accurate determination of the coefficients of elastic crack tip asymptotic field by a hybrid crack element with p-adaptivity. *Engineering Fracture Mechanics* 2001; **68**(15):1609–1630.

26. Xiao Q, Karihaloo B. Implementation of hybrid crack element on a general finite element mesh and in combination with XFEM. *Computer Methods in Applied Mechanics and Engineering* 2007; **196**(13–16):1864–1873.

27. Réthoré J, Roux S, Hild F. Hybrid analytical and extended finite element method (HAX-FEM): a new enrichment procedure for cracked solids. *International Journal for Numerical Methods in Engineering* 2010; **81**(3):269–285.

28. Dhia HB, Rateau G. The Arlequin method as a flexible engineering design tool. *International Journal for Numerical Methods in Engineering* 2005; **62**(11):1442–1462.

29. Passieux JC, Gravouil A, Réthoré J, Baietto MC. Direct estimation of generalized stress intensity factors using a three-scale concurrent multigrid X-FEM. *International Journal for Numerical Methods in Engineering* 2011; **85**(13):1648–1666.

30. Passieux JC, Réthoré J, Gravouil A, Baietto MC. Local/global non-intrusive crack propagation simulation using a multigrid X-FEM solver. *Computational Mechanics* December 2013; **52**(6):1381–1393.

31. Raju I, Newman J. Stress-intensity factors for a wide range of semi-elliptical surface cracks in finite-thickness plates. *Engineering Fracture Mechanics* 1979; **11**(4):817–829.

32. Rannou J, Gravouil A, Baietto-Dubourg MC. A local multigrid X-FEM strategy for 3-D crack propagation. *International Journal for Numerical Methods in Engineering* 2009; **77**(4):581–600.

33. Karihaloo BL, Xiao QZ. Linear and nonlinear fracture mechanics. *Comprehensive Structural Integrity: Fundamental Theories and Mechanisms of Failure* Milne I, Ritchie RO, Karihaloo BL (eds.), Vol. 2. Elsevier: London, 2003; 81–212.

34. Leblond JB, Torlai O. The stress field near the front of an arbitrarily shaped crack in a three-dimensional elastic body. *Journal of Elasticity* 1992; **29**(2):97–131.

35. Chaudhuri RA, Xie M. A novel eigenfunction expansion solution for three-dimensional crack problems. *Composites Science and Technology* 2000; **60**(12–13):2565–2580.

36. González-Albuixech VF, Giner E, Tarancón JE, Fuenmayor FJ, Gravouil A. Convergence of domain integrals for stress intensity factor extraction in 2-D curved cracks problems with the extended finite element method. *International Journal for Numerical Methods in Engineering* 2013; **94**(8):740–757.

37. Belytschko T, Black T. Elastic crack growth in finite elements with minimal remeshing. *International Journal for Numerical Methods in Engineering* 1999; **45**(5):601–620.

38. Sukumar N, Moës N, Moran B, Belytschko T. Extended finite element method for three-dimensional crack modelling. *International Journal for Numerical Methods in Engineering* 2000; **48**(11):1549–1570.

39. Bourgat JF, Glowinski R, Le Tallec P, Vidrascu M. Variational formulation and algorithm for trace operation in domain decomposition calculations. *Technical Report*, 1988.

40. Wang X. Elastic T-stress solutions for penny-shaped cracks under tension and bending. *Engineering Fracture Mechanics* 2004; **71**(16-17):2283–2298.

41. Li S, Mear ME, Xiao L. Symmetric weak-form integral equation method for three-dimensional fracture analysis. *Computer Methods in Applied Mechanics and Engineering* 1998; **151**(3-4):435–459.

42. Mi Y. *Three-Dimensional Analysis of Crack Growth*. WIT Press: Boston, 1996.

43. Qu J, Wang X. Solutions of T-stresses for quarter-elliptical corner cracks in finite thickness plates subject to tension and bending. *International Journal of Pressure Vessels and Piping* 2006; **83**(8):593–606.

44. Bažant ZP, Estenssoro LF. Surface singularity and crack propagation. *International Journal of Solids and Structures* 1979; **15**(5):405–426.

45. Karihaloo BL, Abdalla H, Xiao QZ. Coefficients of the crack tip asymptotic field for wedge splitting specimens. *Engineering Fracture Mechanics* 2003; **70**(17):2407–2420.

46. Lachambre J, Réthoré J, Weck A, Buffiere JY. Extraction of stress intensity factors for 3D small fatigue cracks using digital volume correlation and X-ray tomography. *International Journal of Fatigue* 2014. In press.