# Optimal convergence analysis for the extended finite element method

Serge Nicaise, Yves Renard, Elie Chahine

► To cite this version:

Serge Nicaise, Yves Renard, Elie Chahine. Optimal convergence analysis for the extended finite element method. International Journal for Numerical Methods in Engineering, 2011, 86, pp.528-548. 10.1002/nme.3092 . hal-00339853v2

HAL Id: hal-00339853

https://hal.science/hal-00339853v2

Submitted on 6 Jun 2018

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 CC BY 4.0 - Attribution - International License
# Optimal convergence analysis for the eXtended Finite Element Method

Serge Nicaise<fnref n="1" />, Yves Renard<fnref n="2" />, Elie Chahine<fnref n="3" />

## Abstract

We establish some optimal *a priori* error estimate on some variants of the eXtended Finite Element Method (Xfem), namely the Xfem with a cut-off function and the standard Xfem with a fixed enrichment area. The results are established for the Lamé system (homogeneous isotropic elasticity) and the Laplace problem. The convergence of the numerical stress intensity factors is also investigated. We show some numerical experiments which corroborate the theoretical results.

**Keywords:** extended finite element method, error estimates, stress intensity factors.

## 1 Introduction

Inspired by the Pufem [26], the Xfem (extended finite element method) was introduced by Moës *et al.* in 1999 [28, 27] for plane linear isotropic elasticity problems (Lamé system) in cracked domains. The main advantage of this method is the ability to take into account the discontinuity across the crack and the asymptotic displacement at the crack tip by addition of special functions into the finite element space. It allows the use of a mesh which is independent of the geometry of the crack. This avoids the remeshing operations when the crack propagates and the corresponding re-interpolation operations which can cause numerical instabilities. In the original method, the asymptotic displacement is incorporated into the finite element space multiplied by the shape function of a background Lagrange finite element method. However, we deal also with a variant, introduced in [12], where the asymptotic displacement is multiplied by a cut-off function. This variant is similar to the classical singular enrichment method introduced in 1973 by Strang and Fix [32] but it additionally preserves the independence of the mesh to the geometry of the crack which is indeed the essential contribution of Xfem.

Another classical method to take into account a singular behavior of the solution is the dual singular function method introduced by M. Dobrowolski *et al.* in [5] (see also [19, 10]) or a more recent variant the singular complement method introduced by P. Ciarlet Jr. *et al.* in [17] (for a *L*-shape domain, see [29]). These methods require the use of dual singular functions which can be difficult to obtain in some situations (even for the Lamé system) or quite impossible to obtain when just the asymptotic behavior is known (for non-linear elasticity [2] or Mindlin plate model for instance).

The Xfem strategy can be adapted to various situations. See among many other references [3, 6, 7, 8, 23, 25, 36, 37, 35, 38]. In particular, a fictitious domain method can be

<footnote n="1">¹Université de Valenciennes et du Hainaut Cambrésis, LAMAV, FR CNRS 2956, Institut des Sciences et Techniques de Valenciennes, F-59313 - Valenciennes Cedex 9 France, email: Serge.Nicaise@univ-valenciennes.fr</footnote>

<footnote n="2">²Corresponding author. Université de Lyon, CNRS INSA-Lyon, ICJ UMR5208, LaMCoS UMR5259, F-69621, Villeurbanne, France, Yves.Renard@insa-lyon.fr</footnote>

<footnote n="3">³Laboratory for Nuclear Materials, Nuclear Energy and Safety Research Department, Paul Scherrer Institute OVGA/14, CH-5232 Villigen PSI, Switzerland, elie.chahine@psi.ch.</footnote>
derived from the principle of Xfem (see [24, 4]) and it is possible to adapt some strategies when the asymptotic behavior is unknown or only partially known (see [11, 13, 14]).

In the present paper, we improve the results given in [12] concerning the variant which uses a cut-off function. We also give some additional error estimates concerning the stress intensity factors and the standard Xfem. The theoretical results are established for both the Lamé system and the Laplace problem. Some numerical tests that illustrate and confirm the theoretical results are presented.

## 2 The model problems

The analysis will be performed on a cracked domain $\Omega \subset \mathbb{R}^2$ for two model problems: The Laplace equation and the Lamé system. The crack $\Gamma_C$ is assumed to be straight. In both cases, the boundary $\partial\Omega$ of $\Omega$ is partitioned into $\Gamma_D$, $\Gamma_N$ and $\Gamma_C$ (see Fig. 1). A Dirichlet condition is prescribed on $\Gamma_D$, a Neumann one on $\Gamma_N$ while on the crack $\Gamma_C$ we consider an homogeneous Neumann condition.

<figure><img src="image_2.png" /><figcaption>Figure 1: The cracked domain $\Omega$.</figcaption></figure>

The weak formulation of the (scalar) Laplace equation on this domain reads as follows:

$$ egin{cases} 
\text{Find } u \in V \text{ such that } a(u,v) = l(v) \quad \forall v \in V, \\ 
a(u,v) = \int_{\Omega} \nabla u \cdot \nabla v dx, \\ 
l(v) = \int_{\Omega} f v dx + \int_{\Gamma_N} g v \, d\Gamma, \\ 
V = \{v \in H^1(\Omega); v = 0 \text{ on } \Gamma_D\}. 
\end{cases} 
\quad (1) $$

While one of the Lamé (vectorial) system (linear elasticity problem on this domain for an isotropic material) is:

$$ egin{cases} 
\text{Find } u \in V \text{ such that } a(u,v) = l(v) \quad \forall v \in V, \\ 
a(u,v) = \int_{\Omega} \sigma(u) : \varepsilon(v) \, dx, \\ 
l(v) = \int_{\Omega} f \cdot v \, dx + \int_{\Gamma_N} g \cdot v \, d\Gamma, \\ 
\sigma(u) = \lambda \operatorname{tr}(\varepsilon(u)) I + 2\mu \varepsilon(u), \\ 
V = \{v \in H^1(\Omega; \mathbb{R}^2); v = 0 \text{ on } \Gamma_D\}, 
\end{cases} 
\quad (2) $$

where $\sigma(u)$ denotes the stress tensor, $\varepsilon(u) = \frac{1}{2}(\nabla u + \nabla u^T)$ is the linearized strain tensor,
f and g are some external load densities on Ω and Γ$_N$ respectively, and λ > 0, µ > 0 are the Lamé coefficients.

In both cases, we suppose Ω, f and g smooth enough for the solution u of Problem (1) or (2) to be written as a sum of a singular part u_s and a regular part u − u_s (see [22, 21]) satisfying:

$$
u - u_s \in H^2(\Omega; \mathbb{R}^d), \qquad (3)
$$

with

$$
d = 1 \quad \text{and} \quad u_s = K_L u_L, \tag{4}
$$

for the solution to the Laplace equation (1), and

$$
d = 2 \quad \text{and} \quad u_s = K_I u_I + K_{II} u_{II}, \tag{5}
$$

for the solution to the Lamé system (2). The scalars K_L, K_I and K_{II} are the so-called stress intensity factors and the functions u_L, u_I and u_{II} are given in polar coordinates relatively to the crack tip (Fig. 2) by:

$$
u_L(r, \theta) = \sqrt{r} \sin \frac{\theta}{2}, \qquad (6)
$$

$$
u_I(r, \theta) = \frac{1}{E} \sqrt{\frac{r}{2\pi}} (1 + \nu) \begin{pmatrix} \cos \frac{\theta}{2}(\delta - \cos \theta) \\ \sin \frac{\theta}{2}(\delta - \cos \theta) \end{pmatrix}, \quad (7)
$$

$$
u_{II}(r, \theta) = \frac{1}{E} \sqrt{\frac{r}{2\pi}} (1 + \nu) \begin{pmatrix} \sin \frac{\theta}{2}(\delta + 2 + \cos \theta) \\ \cos \frac{\theta}{2}(\delta - 2 + \cos \theta) \end{pmatrix}, \quad (8)
$$

where $\nu = \frac{\lambda}{\lambda + 2\mu}$ denotes the Poisson ratio, $E = \frac{4\mu(\lambda + \mu)}{\lambda + 2\mu}$ the Young modulus and $\delta = 3 - 4\nu$ (plane stress approximation). Note that $u_L, u_I$ and $u_{II}$ belong to $H^{3/2-\eta}(\Omega)$ for any $\eta > 0$ (see [22]) which limits the order of the convergence rate of a classical finite element method to $O(h^{1/2})$ where $h$ is the mesh parameter.

<figure><img src="image_3.png" /><figcaption>Figure 2: Polar coordinates respectively to the crack tip Ω.</figcaption></figure>

3 Xfem with a cut-off function

The Xfem variant which uses a cut-off function was proposed in [12]. The principle of the
standard Xfem (see [28, 27]) is to consider a mesh independent of the crack geometry. An
Heaveaside type function is used to represent the discontinuity across the straight crack:

$$
H(x) = \begin{cases} +1 & \text{if } (x - x^*) \cdot n \ge 0, \\ -1 & \text{elsewhere,} \end{cases} \tag{9}
$$
where $x^*$ denotes the crack tip and $n$ is a given normal to the crack. Moreover, the nonsmooth functions $u_L$, $u_I$ and $u_{II}$ are integrated to the discrete space to take into account the asymptotic behavior at the crack tip.

<figure><img src="image_4.png" /><figcaption>Figure 3: Enrichment strategy.</figcaption></figure>

We consider an affine Lagrange finite element method defined on a regular triangulation $\mathcal{T}_h$ (in the sense of the Ciarlet [16]) of the non-cracked domain $\bar{\Omega}$, $h$ being the mesh parameter i.e. the largest diameter of the elements of $\mathcal{T}_h$. The piecewise $P_1$ basis functions are denoted $\varphi_i$. In this section, We consider the variant of Xfem proposed in [12] for which a whole area around the crack tip is enriched by using a cut-off function denoted $\chi$. The approximation of the Laplace equation reads as

$$ 
\left\{ 
\begin{array}{l}
\text{Find } u^h \in V^h \quad \text{such that} \quad a(u^h, v^h) = l(v^h) \quad \forall v^h \in V^h, \\
a(u^h, v^h) = \int_{\Omega} \nabla u^h \cdot \nabla v^h dx, \\
l(v^h) = \int_{\Omega} fv^h dx + \int_{\Gamma_N} gv^h d\Gamma, \\
V^h = \left\{ v^h = \sum_{i \in I} a_i \varphi_i + \sum_{i \in I_H} b_i H \varphi_i + K_{L,h} \chi u_L; \ a_i, b_i, K_{L,h} \in \mathbb{R} \right\}.
\end{array}
\right. 
\qquad (10) 
$$

where $I$ is the set of node indices of the $P_1$ finite element method, $I_H$ is the sub-set of node indices whose corresponding shape functions have their supports completely cut by the crack (see Fig. 3) and $\chi$ is a $W^{2,\infty}(\bar{\Omega})$ cut-off function verifying for fixed $0 < r_0 < r_1$

$$ 
\left\{ 
\begin{array}{l}
\chi(r) = 1 \text{ if } r < r_0, \\
0 < \chi(r) < 1 \text{ if } r_0 < r < r_1, \\
\chi(r) = 0 \text{ if } r_1 < r.
\end{array}
\right. 
\qquad (11) 
$$

Note that $I_H$, the set of finite element shape functions enriched by the Heaviside function $H$, cannot be larger. In particular, it cannot contain the shape functions having the element which contains the crack tip in their support. Otherwise, the geometry of the crack would not be well represented since it is prolonged one element forward. In fact, the representation of the discontinuity across the crack inside the element containing the crack tip is ensured by the nonsmooth functions.
Concerning now the Lamé system, we consider two different ways to incorporate the asymptotic displacement. The first one is directly based on a vectorial enrichment with $u_I$ and $u_{II}$:

$$ egin{cases} 
\text{Find } u^h \in V^h \text{ such that } a(u^h, v^h) = l(v^h) \quad \forall v^h \in V^h, \\
a(u^h, v^h) = \int_{\Omega} \sigma(u^h) : \varepsilon(v^h) \, dx, \\
l(v^h) = \int_{\Omega} f \cdot v^h \, dx + \int_{\Gamma_N} g \cdot v^h \, d\Gamma, \\
\sigma(u^h) = \lambda \operatorname{tr} \varepsilon(u^h) I + 2\mu \varepsilon(u^h), \\
V^h = \left\{ v^h = \sum_{i \in I} a_i \varphi_i + \sum_{i \in I_H} b_i H \varphi_i + K_{I,h} \chi u_I + K_{II,h} \chi u_{II}; \ a_i, b_i \in \mathbb{R}^2, K_{I,h}, K_{II,h} \in \mathbb{R} \right\}. 
\end{cases} 
\quad (12) $$

The second one corresponds to a more classical Xfem approximation with a scalar enrichment of each component:

$$ egin{cases} 
\text{Find } u^h \in V^h \text{ such that } a(u^h, v^h) = l(v^h) \quad \forall v^h \in V^h, \\
a(u^h, v^h) = \int_{\Omega} \sigma(u^h) : \varepsilon(v^h) \, dx, \\
l(v^h) = \int_{\Omega} f \cdot v^h \, dx + \int_{\Gamma_N} g \cdot v^h \, d\Gamma, \\
\sigma(u^h) = \lambda \operatorname{tr} \varepsilon(u^h) I + 2\mu \varepsilon(u^h), \\
V^h = \left\{ v^h = \sum_{i \in I} a_i \varphi_i + \sum_{i \in I_H} b_i H \varphi_i + \sum_{j=1}^4 c_j F_j \chi; \ a_i, b_i, c_j \in \mathbb{R}^2 \right\}, 
\end{cases} 
\quad (13) $$

where the set of functions $\{F_j(x)\}_{1 \le j \le 4}$ is defined by

$$ \{F_j(x)\}_{1 \le j \le 4} = \{ \sqrt{r} \sin \frac{\theta}{2}, \sqrt{r} \cos \frac{\theta}{2}, \sqrt{r} \sin \frac{\theta}{2} \cos \theta, \sqrt{r} \cos \frac{\theta}{2} \cos \theta \}. \quad (14) $$

Note that the nonsmooth functions $u_I$ and $u_{II}$ can be decomposed on this set of functions.

## 4 Optimal Error estimate for the Xfem with a cut-off function

We use the notation $a \lesssim b$ to signify that there exists a constant $C > 0$ independent of the mesh parameter of the solution and of the crack-tip position such that $a \le Cb$. For a non negative real number $s$ let $H^s(D)$ denote the standard Sobolev space of order $s$ in $D$ of norm (resp. semi-norm) denoted by $\| \cdot \|_{s,D}$ (resp. $| \cdot |_{s,D}$), see for instance [1].

The aim of this section is to establish the following result which is the optimal version of Theorem 1 in [12]:

**Theorem 1** Assume that the displacement field $u$, solution to Problem (1) (resp. Problem (2)), satisfies Condition (3). Then, the following estimate holds

$$ \|u - u^h\|_{1,\Omega} \lesssim h \|u - \chi u_s\|_{2,\Omega}, \quad (15) $$

where $u^h$ is the solution to Problem (10) (resp. to Problem (12) or to Problem (13)), $u_s$ is the singular part of $u$ (see (3)) and $\chi$ is the $W^{2,\infty}(\Omega)$ cut-off function introduced before.
The outline of the proof globally follows the one of Theorem 1 in [12]. Some sub-optimal
intermediary results are here replaced by optimal ones.

We recall the definition of the adapted interpolation operator Π^h. The interpolation error estimates are then computed locally over every different type of triangles: triangles totally enriched by the Heaveaside function, triangles partially enriched by the Heaviside function and the triangle containing the crack-tip.

<figure><img src="image_5.png" /><figcaption>Figure 4: Domain decomposition.</figcaption></figure>

The domain $\Omega$ is divided into $\Omega_1$ and $\Omega_2$ according to the crack and a straight extension of the crack (Fig.4) such that the value of $H$ is $(-1)^k$ on $\Omega_k$, $k = 1, 2$. Let us denote $u_r = u - \chi u_s$, and $u_r^k$ the restriction of $u_r$ to $\Omega_k$, $k \in \{1, 2\}$. Then, there exists in $H^2(\bar{\Omega}; \mathbb{R}^d)$ an extension $\tilde{u}_r^k$ of $u_r^k$ on $\bar{\Omega}$ such that (see [1])

$$
\|\tilde{u}_r^1\|_{2,\bar{\Omega}} \lesssim \|u_r^1\|_{2,\Omega_1}, \quad (16)
$$

$$
\| \tilde{u}_r^2 \|_{2, \bar{\Omega}} \lesssim \| u_r^2 \|_{2, \Omega_2}. \quad (17)
$$

**Definition 1** (from [12]) Given a displacement field $u$ satisfying (3) and two extensions $\tilde{u}_r^1$ and $\tilde{u}_r^2$ respectively of $u_r^1$ and $u_r^2$ in $H^2(\bar{\Omega}; \mathbb{R}^d)$, we define $\Pi^h u$ as the element of $V^h$ such that

$$
\Pi^h u = \sum_{i \in I} a_i \varphi_i + \sum_{i \in I_H} b_i H \varphi_i + \chi u_s, \quad (18)
$$

where $a_i$, $b_i$ are given as follows ($x_i$ denotes the node associated to $\varphi_i$):

$$
\begin{align}
& \text{if } i \in \{I \setminus I_H\} \text{ then } a_i = u_r(x_i), \nonumber \\
& \text{if } i \in I_H \text{ and } x_i \in \bar{\Omega}_k \text{ then } (k \in \{1,2\}, l \neq k) \quad \left\{ \begin{array}{@{}l@{}}
      a_i = \frac{1}{2} (u_r^k(x_i) + \tilde{u}_r^l(x_i)), \\
      b_i = \frac{1}{2} (u_r^k(x_i) - \tilde{u}_r^l(x_i)) (-1)^k.
      \end{array} \right. \tag{19}
\end{align}
$$

From this definition, the following result holds:

**Lemma 1** (from [12]) The function Π^h u satisfies

(i) Π^h u = I^h u_r + χu_s over a triangle non-enriched by H,

(ii) Π^h u|_{K∩Ω_k} = I^h \tilde{u}_r^k + χu_s over a triangle K totally enriched by H,

where I^h denotes the classical interpolation operator for the associated finite element method.
<figure><img src="image_6.png" /><figcaption>Figure 5: (a) Totally enriched triangle and (b) partially enriched triangle (Fig.3).</figcaption></figure>

For K a subset of Ω, we denote $h_K = \text{diam}(K) = \max_{x_1,x_2 \in K}|x_1 - x_2|$ and $\rho_K = \{\text{sup}(\text{diam}(B)); B \text{ ball of } \mathbb{R}^2, B \subset K\}$. The following lemma, established in [12], derives simply from the classical interpolation of the extensions of $u_r^1$ and $u_r^2$.

**Lemma 2** (from [12]) Let $\mathcal{T}_h^H$ be the set of triangles totally enriched by H (Fig.3) and $\sigma_K = h_K \rho_K^{-1}$. For all K in $\mathcal{T}_h^H$, and for all u satisfying (3) we have the estimates

$$ \|u - \Pi^h u\|_{1,K \cap \Omega_1} \lesssim h_K \sigma_K \|\tilde{u}_r^1\|_{2,K}, \quad (20) $$

and

$$ \|u - \Pi^h u\|_{1,K \cap \Omega_2} \lesssim h_K \sigma_K \|\tilde{u}_r^2\|_{2,K}. \quad (21) $$

The optimal convergence is of course obtained for non enriched triangles. It remains to treat the partially enriched triangles and the triangle containing the crack tip. We will now detailled the optimal intermediary results which are original in this paper.

Let us start with the Laplace equation and recall that in that case $u_r^1$ and $u_r^2$ satisfy

$$ \partial_n u_r^1 = \partial_n u_r^2 = 0 \text{ on } \Gamma_C. \quad (22) $$

Since the extension $\tilde{u}_r^1$ and $\tilde{u}_r^2$ are $H^2(\bar{\Omega})$ extension, they also satisfy the Neumann boundary condition on $\Gamma_C$, namely

$$ \partial_n \tilde{u}_r^1 = \partial_n \tilde{u}_r^2 = 0 \text{ on } \Gamma_C. \quad (23) $$

We now give the main technical result.

**Lemma 3** Assume that $x_1 \in \Omega_1$ is a node belonging to a triangle K containing the crack tip. Then

$$ |u_r^1(x_1) - \tilde{u}_r^2(x_1)| \lesssim h_K |\tilde{u}_r^1 - \tilde{u}_r^2|_{2, B(0, h_K)}. \quad (24) $$

**Proof:** For shortness write $v = \tilde{u}_r^1 - \tilde{u}_r^2$. Using a Taylor expansion, we have

$$ v(x_1) = \int_0^1 (x_1^{(1)} \partial_1 v(tx_1) + x_1^{(2)} \partial_2 v(tx_1)) dt, $$

where $x_1 = (x_1^{(1)}, x_1^{(2)})$. Without loss of generality and modulo an orthonormal change of coordinates we assume that the position of the crack tip is (0,0) and the crack $\Gamma_C$ is a part of $(\mathbb{R}_-, 0)$. By setting $v^{(1)} = \partial_1 v$ and $v^{(2)} = \partial_2 v$, and making the change of variable $s = tx_1$, the above identity is equivalent to

$$ v(x_1) = \int_e (n_2 v^{(1)}(s) - n_1 v^{(2)}(s)) ds, \quad (25) $$
where e is the edge joining the crack tip and $x_1$ and $n = (n_1, n_2)$ is (one of) the normal vector to e. Denote by C the truncated sector determined by e and the crack:

$$C = \{(r \cos \theta, r \sin \theta) : 0 < r < h_1 \quad \theta_0 < \theta < \pi\},$$

when $x_1 = (h_1 \cos \theta_0, h_1 \sin \theta_0)$, see Fig. 6. Now setting $e_2 = \{(h_1 \cos \theta, h_1 \sin \theta) : \theta_0 < \theta < \pi\}$, by Green's formula we remark that

$$\int_C \partial_1 v^{(2)} dx = \int_{\partial C} n_1 v^{(2)} ds = \int_e n_1 v^{(2)} ds + \int_{e_2} n_1 v^{(2)} ds,$$

because $n_1 = 0$ on $\Gamma_C$. Hence

$$\int_e n_1 v^{(2)} ds = \int_C \partial_1 v^{(2)} dx - \int_{e_2} n_1 v^{(2)} ds. \qquad (26)$$

The first term of this right-hand side will be estimated by a simple Cauchy-Schwarz inequality. For the second term by a scaling argument, we show that

$$\int_{e_2} |v^{(2)}| ds \lesssim h_1 \|\nabla v^{(2)}\|_{B(0,h_1)}. \qquad (27)$$

<figure><img src="image_7.png" /><figcaption>Figure 6: The truncated sector C.</figcaption></figure>

Indeed by construction $v^{(2)}$ satisfies

$$v^{(2)} = 0 \text{ on } \Gamma_C.$$

Therefore the change of variable $x = h_1 \hat{x}$ maps $B(0, h_1)$ to the unit ball. By setting $\hat{v}^{(2)}(\hat{x}) = v^{(2)}(x)$, we deduce that

$$\begin{align*}
\int_{e_2} |v^{(2)}| ds &\le \int_{\partial B(0, h_1)} |v^{(2)}| ds \\
&= h_1 \int_{\partial B(0, 1)} |\hat{v}^{(2)}| d\hat{s} \\
&\lesssim h_1 \left( \int_{B(0, 1)} |\nabla \hat{v}^{(2)}|^2 d\hat{x} \right)^{\frac{1}{2}}.
\end{align*}$$

This last estimate follows from the property

$$\hat{v}^{(2)}(\hat{x}) = 0 \text{ on } \{(x_1, 0) : -1 < x_1 < 0\},$$
and the compact embedding of $H^1(B(0, 1))$ into $L^2(B(0, 1))$. Coming back to $B(0, h_1)$, we obtain (27).

Using the estimate (27) into (26) and Cauchy-Schwarz inequality, we have shown that

$$ \left| \int_e n_1 v^{(2)} ds \right| \lesssim h_1 \| \nabla v^{(2)} \|_{B(0, h_1)}. \qquad (28) $$

Let us now pass to the estimate of $\int_e n_2 v^{(1)}(s) ds$: Denote by $D$ the truncated sector determined by $e$ and the extended crack:

$$ D = \{(r \cos \theta, r \sin \theta) : 0 < r < h_1 \quad 0 < \theta < \theta_0\}, $$

when we recall that $x_1 = (h_1 \cos \theta_0, h_1 \sin \theta_0)$ (see Fig. 6). As before setting $e_3 = \{(h_1 \cos \theta, h_1 \sin \theta) : 0 < \theta < \theta_0\}$, we remark that

$$ \int_D \partial_2 v^{(1)} dx = \int_{\partial D} n_2 v^{(1)} ds = \int_e n_2 v^{(1)} ds + \int_{e_3} n_2 v^{(1)} ds, $$

because $v^{(1)} = 0$ on $\tilde{\Gamma}_C = \{(x_1, 0) : x_1 > 0\}$, the extension of the crack $\Gamma_C$ to $x_1 > 0$. Hence

$$ \int_e n_2 v^{(1)} ds = \int_D \partial_2 v^{(1)} dx - \int_{e_3} n_2 v^{(1)} ds. \qquad (29) $$

It suffices to estimate the second term of this right-hand side. Again using a scaling argument, we show that

$$ \int_{e_3} |v^{(1)}| ds \lesssim h_1 \|\nabla v^{(1)}\|_{B(0, h_1)}. \qquad (30) $$

Indeed by construction $v^{(1)}$ satisfies

$$ v^{(1)} = 0 \text{ on } \tilde{\Gamma}_C. $$

Therefore the same scaling arguments as before lead to (30).

Using the estimate (30) into (29) and Cauchy-Schwarz inequality, we have shown that

$$ \left| \int_e n_2 v^{(1)} ds \right| \lesssim h_1 \| \nabla v^{(1)} \|_{B(0, h_1)}. \qquad (31) $$

The estimates (28) and (31) into the identity (25) lead to the estimate (24) because $h_1 \le h_K$. $\square$

Let us go on with the Lamé system and recall that in that case $u_r^1$ and $u_r^2$ satisfy

$$ \sigma(u_r^1) \cdot n = \sigma(u_r^2) \cdot n = 0 \text{ on } \Gamma_C. \qquad (32) $$

Since the extension $\tilde{u}_r^1$ and $\tilde{u}_r^2$ belong to $H^2(\bar{\Omega}; \mathbb{R}^2)$, they also satisfy the traction free boundary condition on $\Gamma_C$, namely

$$ \sigma(\tilde{u}_r^1) \cdot n = \sigma(\tilde{u}_r^2) \cdot n = 0 \text{ on } \Gamma_C. \qquad (33) $$

**Lemma 4** Assume that $x_1 \in \Omega_1$ is a node belonging to a triangle $K$ containing the crack tip. Then the estimate (24) holds.
**Proof:** The proof starts as before with the identity (25).

We notice that by (33) and since $n = (0, 1)^{\top}$ on $\Gamma_C$, $v = \tilde{u}_r^1 - \tilde{u}_r^2$ satisfies

$$ \lambda(\partial_1 v_1 + \partial_2 v_2) + \mu(\partial_1 v_2 + \partial_2 v_1) = 0 \text{ on } \Gamma_C, \quad (34) $$

$$ (\lambda + 2\mu)\partial_2 v_2 + \lambda\partial_1 v_1 = 0 \text{ on } \Gamma_C, \quad (35) $$

where $v_1, v_2$ are the two components of $v$, i.e., $v = (v_1, v_2)^{\top}$. Note that by construction, we also have $v = 0$ on $\tilde{\Gamma}_C$ and therefore

$$ \partial_1 v_1 = \partial_1 v_2 = 0 \text{ on } \tilde{\Gamma}_C. \quad (36) $$

Since $v^{(1)}$ still satisfies

$$ v^{(1)} = 0 \text{ on } \tilde{\Gamma}_C, $$

the arguments of the previous lemma show that (31) is valid.

For the estimate of the term involving $v^{(1)}$, since $n_1 = 0$ on $\Gamma_C$, as before the identity (26) holds. To estimate the second term of the right-hand side of (26), we again use a scaling argument: The change of variable $x = h_1 \hat{x}$ maps $B(0, h_1)$ to the unit ball and by setting $\hat{w}(\hat{x}) = \nabla v(x)$, where

$$ \nabla v = \begin{pmatrix} \partial_1 v_1 & \partial_2 v_1 \\ \partial_1 v_2 & \partial_2 v_2 \end{pmatrix}, $$

we deduce that

$$ \begin{aligned} \int_{e_2} |v^{(2)}| ds &\le \int_{\partial B(0, h_1)} |v^{(2)}| ds \\ &\le \int_{\partial B(0, h_1)} |\nabla v| ds \\ &\le h_1 \int_{\partial B(0, 1)} |\hat{w}| ds. \end{aligned} $$

Now we notice that the conditions (34), (35) and (36) satisfied by $v$ lead to

$$ \lambda(\hat{w}_{11} + \hat{w}_{22}) + \mu(\hat{w}_{21} + \hat{w}_{12}) = 0 \text{ on } \{(x_1, 0) : -1 < x_1 < 0\}, \quad (37) $$

$$ (\lambda + 2\mu)\hat{w}_{22} + \lambda\hat{w}_{11} = 0 \text{ on } \{(x_1, 0) : -1 < x_1 < 0\}, \quad (38) $$

$$ \hat{w}_{11} = \hat{w}_{12} = 0 \text{ on } \{(x_1, 0) : 0 < x_1 < 1\}. \quad (39) $$

Hence the compact embedding of $H^1(B(0, 1))$ into $L^2(B(0, 1))$ and a contradiction argument lead to

$$ \int_{\partial B(0, 1)} |\hat{w}| ds \lesssim \|w\|_{1, B(0, 1)} \lesssim |w|_{1, B(0, 1)}. $$

This last estimate holds since otherwise we would find a vector field $v \in H^1(B(0, 1))^{2 \times 2}$ satisfying (37) to (39) such that

$$ |w|_{1, B(0, 1)} = 0 \text{ and } \|w\|_{0, B(0, 1)} = 1. $$

Such a matrix field does not exist because $w$ would be a constant matrix and by (37) to (39), it would be zero.

This estimate leads to (27) and we conclude as in the previous Lemma.

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

These lemmas allow to treat the non-optimal cases from [12] as follows:
**Corollary 1** Let K be a triangle partially enriched and let $K^* = K \setminus \Gamma_C$. Then

$$ \|u - \Pi^h u\|_{1,K^*} \lesssim h_K(|\tilde{u}_r^1|_{2,B(0,2h_K)} + |\tilde{u}_r^2|_{2,B(0,2h_K)}). \quad (40) $$

**Proof:** It is sufficient to estimate $\|u_r - \Pi^h u_r\|_{1,K^*}$ since the singular part of $u - \Pi^h u$ vanishes. We treat the situation of Fig. 5 (b). Other situations can be treated similarly. We have

$$ \Pi^h u_r = u_r^1(x_1)\varphi_1 + u_r^2(x_2)\varphi_2 + \tilde{u}_r^2(x_3)\varphi_3 \text{ on } K_2 = K \cap \Omega_2, $$

or equivalently

$$ \begin{aligned} \Pi^h u_r &= \tilde{u}_r^2(x_1)\varphi_1 + u_r^2(x_2)\varphi_2 + \tilde{u}_r^2(x_3)\varphi_3 + (u_r^1(x_1) - \tilde{u}_r^2(x_1))\varphi_1 \text{ on } K_2 \\ &= \Pi^h \tilde{u}_r^2 + (u_r^1(x_1) - \tilde{u}_r^2(x_1))\varphi_1 \text{ on } K_2. \end{aligned} $$

By the triangular inequality, we may write

$$ \begin{aligned} \|u_r - \Pi^h u_r\|_{1,K_2} &\le \|u_r^2 - \Pi^h \tilde{u}_r^2\|_{1,K_2} + |u_r^1(x_1) - \tilde{u}_r^2(x_1)|\|\varphi_1\|_{1,K_2} \\ &\lesssim \|\tilde{u}_r^2 - \Pi^h \tilde{u}_r^2\|_{1,K} + |u_r^1(x_1) - \tilde{u}_r^2(x_1)|. \end{aligned} $$

By a standard interpolation error estimate and Lemma 3 (or 4), we conclude that

$$ \|u_r - \Pi^h u_r\|_{1,K_2} \lesssim h_K(|\tilde{u}_r^2|_{2,K} + |\tilde{u}_r^1 - \tilde{u}_r^2|_{2,B(0,h_K)}). $$

For the part on $K_1 = K \cap \Omega_1$, we remark that

$$ \begin{aligned} \Pi^h u_r &= u_r^1(x_1)\varphi_1 + u_r^2(x_2)\varphi_2 + u_r^1(x_3)\varphi_3 \text{ on } K_1 \\ &= \Pi^h \tilde{u}_r^1 + (\tilde{u}_r^1(x_1) - u_r^2(x_1))\varphi_2 \text{ on } K_2. \end{aligned} $$

And we conclude as before because $\tilde{u}_r^1 - u_r^2$ satisfies the same conditions than $u_r^1 - \tilde{u}_r^2$ on $\Gamma_C$ and $\tilde{\Gamma}_C$. $\square$

**Corollary 2** Let K be the triangle containing the crack tip. Then

$$ \|u - \Pi^h u\|_{1,K^*} \lesssim h_K(|\tilde{u}_r^1|_{2,B(0,h_K)} + |\tilde{u}_r^2|_{2,B(0,h_K)}). \quad (41) $$

**Proof:** In this case we have

$$ \Pi^h u_r = u_r^1(x_1)\varphi_1 + u_r^2(x_2)\varphi_2 + u_r^2(x_3)\varphi_3 \text{ on } K. $$

Without loss of generality, we may assume that K has one vertex $x_1$ in $\Omega_1$ and the two other ones $x_2, x_3$ in $\Omega_2$. In this case on $K_1 = K \cap \Omega_1$, we have

$$ \begin{aligned} \|u_r - \Pi^h u_r\|_{1,K_1} &\le \|u_r^1 - \Pi^h \tilde{u}_r^1\|_{1,K_1} + |u_r^1(x_2) - \tilde{u}_r^1(x_2)|\|\varphi_2\|_{1,K_1} + |u_r^1(x_3) - \tilde{u}_r^1(x_3)|\|\varphi_3\|_{1,K_1} \\ &\lesssim |\tilde{u}_r^1 - \Pi^h \tilde{u}_r^1|_{1,K} + \|u_r^1(x_2) - \tilde{u}_r^1(x_2)\| + |u_r^1(x_3) - \tilde{u}_r^1(x_3)|. \end{aligned} $$

We then conclude as in the previous Corollary. The estimate on $K_2 = K \cap \Omega_2$ is treated similarly. $\square$

As in [12], these two Corollaries and Lemma 2 lead to the global error estimate of Theorem 1.
## 5 Optimal error estimate for the standard Xfem

We give now an a priori error estimate for the standard Xfem with a fixed enrichment area. In the original method proposed in [28] the enrichment with the asymptotic displacement at the crack tip is done only on the element containing the crack tip. The rate of convergence of this method is the same than the one without the enrichment (i.e. $O(\sqrt{h})$, see [31, 25]) since the area of enrichment tends to vanish when the mesh parameter decreases. Of course, this rate of convergence is not difficult to establish. Instead, we prove here an optimal error estimate for the strategy introduced independently in [25] and [3] and called “Xfem with a fixed enrichment area” in the first reference and “Xfem with geometrical enrichment” in the second one and consisting in an enrichment area for the asymptotic displacement whose size is independent of the mesh parameter. The approximation of the Laplace equation with this method reads as

$$ egin{cases} 
\text{Find } u^h \in V^h \text{ such that } a(u^h, v^h) = l(v^h) \quad \forall v^h \in V^h, \\ 
a(u^h, v^h) = \int_{\Omega} \nabla u^h \cdot \nabla v^h dx, \\ 
l(v^h) = \int_{\Omega} f v^h dx + \int_{\Gamma_N} g v^h d\Gamma, \\ 
V^h = \left\{ v^h = \sum_{i \in I} a_i \varphi_i + \sum_{i \in I_H} b_i H \varphi_i + \sum_{i \in I_F} c_i \varphi_i F_1; \ a_i, b_i, c_i \in \mathbb{R} \right\}, 
\end{cases} 
\quad (42) $$

and the one of the Lamé system is:

$$ egin{cases} 
\text{Find } u^h \in V^h \text{ such that } a(u^h, v^h) = l(v^h) \quad \forall v^h \in V^h, \\ 
a(u^h, v^h) = \int_{\Omega} \sigma(u^h) : \varepsilon(v^h) \, dx, \\ 
l(v^h) = \int_{\Omega} f \cdot v^h \, dx + \int_{\Gamma_N} g \cdot v^h \, d\Gamma, \\ 
\sigma(u^h) = \lambda \operatorname{tr} \varepsilon(u^h) I + 2\mu \varepsilon(u^h), \\ 
V^h = \left\{ v^h = \sum_{i \in I} a_i \varphi_i + \sum_{i \in I_H} b_i H \varphi_i + \sum_{i \in I_F} \sum_{j=1}^{4} c_{i,j} \varphi_i F_j; \ a_i, b_i, c_{i,j} \in \mathbb{R}^2 \right\}. 
\end{cases} 
\quad (43) $$

where $I_F$ is the set of finite element nodes which are inside a disk centered on the crack tip and of a fixed radius $r_2$ independent of the mesh parameter. Let us prove now the optimality of this method.

**Theorem 2** Assume that the displacement field $u$, solution to Problem (1) (resp. Problem (2)), satisfies Condition (3). Then, the following estimate holds

$$ \|u - u^h\|_{1,\Omega} \lesssim h (\|u - u_s\|_{2,\Omega} + \|u_s\|_{1,\Omega} + \|u_s\|_{2,\Omega \setminus B(x^*, \frac{r_2}{2})}), \quad (44) $$

where $u^h$ is the solution to Problem (42) (resp. to Problem (43)).

**Proof:** Let $\chi$ be a $W^{2,\infty}$ cut-off function satisfying (11) such that $r_1 < r_2$ and $r_0 > \frac{r_2}{2}$. Let

$$ \chi^h = I^h \chi, $$

be the interpolate of $\chi$ on the $P_1$ finite element method. Using the notation of Section 4, the following interpolation operator

$$ \Pi_S^h u = I^h u_r + \chi^h u_s, $$
clearly satisfies $\Pi_S^h u \in V^h$ for $V^h$ defined by (42) (resp. by (43)) at least for $h$ sufficiently small since $r_1 < r_2$. Then,

$$ \begin{aligned} \|\Pi^h u - \Pi_S^h u\|_{1,\Omega} &= \|(\chi - \chi^h)u_s\|_{1,\Omega} \\ &\le \|(\chi - \chi^h)\|_{W^{1,\infty}} \|u_s\|_{1,\Omega} \\ &\lesssim h\|\chi\|_{W^{2,\infty}} \|u_s\|_{1,\Omega}, \end{aligned} $$

using a classical error estimate on the interpolation of $\chi$ with a $P_1$ finite element method (see for instance [20]). Thus, using the estimates established in section 4, one has

$$ \begin{aligned} \|u - \Pi_S^h u\|_{1,\Omega} &\le \|u - \Pi^h u\|_{1,\Omega} + \|\Pi^h u - \Pi_S^h u\|_{1,\Omega} \\ &\lesssim h\|u - \chi u_s\|_{2,\Omega} + h\|u_s\|_{1,\Omega} \\ &\lesssim h(\|u - u_s\|_{2,\Omega} + \|1 - \chi\|_{W^{2,\infty}} \|u_s\|_{1,\Omega} + \|u_s\|_{2,\Omega \setminus B(x^*, \frac{r_2}{2})}) \end{aligned} $$

Which ends the proof thanks to Céa's lemma. □

An interpretation of the proof of Theorem 2 is that the standard Xfem performs a better approximation than the Xfem with a cut-off function because the cut-off function used in the proof is arbitrary. As a consequence, the error bound of the standard Xfem is less than the infimum taken on all the $W^{2,\infty}$ cut-off functions satisfying (11). This is also corroborated with the result on Fig. 8. Nevertheless, the standard Xfem is more expensive than the Xfem with a cut-off function since the number of enrichment degrees of freedom can be greatly higher (see next section).

# 6 Numerical experiments

The analysis presented in the two previous sections is corroborated by the numerical tests also presented in [12]. We reproduce on Fig. 8 the convergence curves obtained in this paper for the approximation (13).

These numerical tests were done on a non-cracked domain defined by $\bar{\Omega} = [-0.5; 0.5] \times [-0.5; 0.5]$ and the crack was the line segment $\Gamma_C = [-0.5; 0] \times \{0\}$. The cut-off function $\chi \in C^2(\bar{\Omega})$ was defined such that

$$ \left\{ \begin{array}{ll} \chi(r) = 1 & \text{if } r < r_0 = 0.01, \\ \chi(r) = 0 & \text{if } r > r_1 = 0.49, \end{array} \right. \qquad (45) $$

and $\chi$ was identical to a fifth degree polynomial for $r_0 \le r \le r_1$.
<figure><img src="image_9.png" /><figcaption>Figure 7: The reference solution, a mixed mode (from [12]).</figcaption></figure>

The exact solution was a combination of a regular solution to the elasticity problem, the mode I and the mode II analytical solutions and a higher order mode (for the deformed configuration, see Fig. 7 with the Von Mises stress). Fig. 8 shows a comparison of the convergence curves of the non-enriched classical method, the standard Xfem and the cut-off strategy. The optimal rate is obtained for both the cut-off enrichment and the standard Xfem with a fixed enrichment area.

<figure><img src="image_10.png" /><figcaption>Figure 8: $H^1$ error with respect to the number of cells in each direction (ns) for a mixed mode and different enrichment strategies of a P1 elements (from [12]).</figcaption></figure>

Table 6 is also reproduced from [12]. It shows a comparison between the number of
the degrees of freedom for different refinements of the classical method, the XFEM with a fixed enrichment area and the cutoff method. In accordance with the theoretical analysis presented in Section 5, the standard Xfem is more accurate but needs more degrees of freedom than the variant with the cut-off function. One might also remark that the choice between the two variants should account for the fact that the stiffness matrix is sparser for the standard Xfem.

<table><caption>Table 1: Number of degrees of freedom</caption><thead><tr><th>Number of cells<br>in each direction</th><th>Classical FEM</th><th>XFEM<br>(enrichment radius =0.2)</th><th>Cutoff enrichment</th></tr></thead><tbody><tr><td>40</td><td>3402</td><td>4962</td><td>3410</td></tr><tr><td>60</td><td>7502</td><td>11014</td><td>7510</td></tr><tr><td>80</td><td>13202</td><td>19578</td><td>13210</td></tr></tbody></table>

Fig. 9 and 10 present some new numerical tests on the comparison between strategies (12) and (13) (i.e. between a scalar and a vectorial enrichment) for the same experimental situation. The discrete space corresponding to the scalar enrichment (13) strictly includes the one for the vectorial enrichment (12). However, the gain in $H^1(\Omega)$ norm for the error is rather small (Fig. 9). Consequently, the vectorial enrichment appears to be a better choice since the number of additional degrees of freedom is lower and the condition number of the linear system obtained is better (Fig. 10).

<figure><img src="image_11.png" /><figcaption>Figure 9: $H^1$ error with respect to the number of cells in each direction (ns). Comparison of strategies (12) and (13).</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>Figure 10: Condition number of the linear system. Comparison of strategies (12) and (13).</figcaption></figure>

# 7 Error estimate on the stress intensity factor

In this section, we show an error estimate between the exact stress intensity factors and the approximated ones. Let us start with the Laplace equation. Recall that we write

$$u = u_r + K_L \chi u_L,$$

and that our Galerkin solution $u_h \in V_h$ solution to (12) admits the splitting

$$u_h = u_{rh} + K_{L,h} \chi u_L,$$

where $u_{rh} \in S^h$, the space $S^h$ being defined by

$$S^h = \left\{ v^h = \sum_{i \in I} a_i \varphi_i + \sum_{i \in I_H} b_i H \varphi_i; a_i, b_i \in \mathbb{R} \right\},$$

so that our approximation space $V_h$ is spanned by $S^h$ plus the singular function $\chi u_L$.

Adapting the arguments from Theorem 9.1 of [9] we have the next error estimate:

**Theorem 3** Assume that the triangulation is quasi-uniform in the sense that

$$h \lesssim h_K \quad \forall K \in \mathcal{T}_h.$$

Then we have

$$|K_L - K_{L,h}| \lesssim h^{\frac{1}{2}}. \qquad (46)$$

**Proof:** As in Theorem 9.1 of [9], we have

$$K_L - K_{L,h} = - \frac{a((I - G_h)u_r, (I - G_h)(\chi u_L))}{a((I - G_h)(\chi u_L), (I - G_h)(\chi u_L))},$$

where $G_h u$ is the Galerkin approximation of $u$ on $S^h$, namely $G_h u \in S^h$ is the unique solution of

$$a(G_h u, v_h) = a(u, v_h) \quad \forall v_h \in S^h.$$
By Cauchy-Schwarz's inequality, we deduce that

$$|K_L - K_{L,h}| \le \frac{\|(I - G_h)u_r\|_{1,\Omega}}{\|(I - G_h)(\chi u_L)\|_{1,\Omega}}, \quad (47)$$

Since $u_r$ belongs to $H^2(\Omega)$ by Theorem 1, we have

$$\|(I - G_h)u_r\|_{1,\Omega} \lesssim h|u_r|_{2,\Omega}, \quad (48)$$

and it remains to estimate from below the denominator of (47). For that purpose, we need to adapt the arguments from Lemma 7.1 of [9] because here the triangulation is not aligned with the crack. The main point is to find a small truncated cone $C_\rho$ included into the triangle $K$ containing the crack tip with $\rho$ equivalent to $h$. Let us denote by $x_i, i = 1, 2, 3$ the three nodes of $K$. First we remark that by a scaling argument we have

$$\max_{i=1,2,3} |x_i| \ge \frac{\rho_K}{\sqrt{2}} \max_{i=1,2,3} |\hat{x}_i - \hat{O}|,$$

where $|x_i|$ is the Euclidean norm of $x_i$, $\hat{O}$ is the pull back of the crack tip $O$ by the affine transformation $F_K$ that sends the standard reference element $\hat{K}$ to $K$. Simple calculations show that

$$\max_{i=1,2,3} |\hat{x}_i - \hat{O}| \ge \frac{1}{4},$$

and therefore since the triangulation is regular we have

$$\max_{i=1,2,3} |x_i| \gtrsim h_K.$$

We now fix $j \in \{1, 2, 3\}$ such that

$$|x_j| = \max_{i=1,2,3} |x_i| \gtrsim h_K.$$

Let $e_1$ and $e_2$ be the two edges of $K$ having $x_j$ as vertex and denote by $\gamma_\ell, \ell = 1, 2$, the angle between $e_\ell$ and the segment joining $x_j$ to $O$. Without loss of generality we may assume that $\gamma_1 \ge \gamma_2$, and therefore

$$\gamma_1 \ge \frac{\alpha_0}{2},$$

where $\alpha_0 \in (0, \frac{\pi}{3})$ is the minimal angle of all triangles of $\mathcal{T}_h$ (equivalent to the regularity of the mesh thanks to Zlamal's result [39]).

<figure><img src="image_13.png" /><figcaption>Figure 11: Sub-triangle $K_1$.</figcaption></figure>
We now consider the sub-triangle $K_1$ of $K$ of vertices $O, x_j, m$, where $m$ is the mid-point of the edge $e_1$. Denote that $\alpha$ and $\beta$ the angle of $K_1$ at $O$ and $m$ respectively (see Fig. 11). Now if $2l$ is the length of the edge $e_1$, by the sinus formula, we notice that

$$ \frac{\sin \alpha}{\sin \beta} = \frac{l}{|x_j|} \sim 1. $$

This property and the fact that

$$ 2\alpha_0 \le \alpha + \beta = \pi - \gamma_1 \le \pi - \frac{\alpha_0}{2}, $$

leads to the existence of a minimal angle $\alpha_1 > 0$ (independent of $h$) such that

$$ \alpha > \alpha_1. $$

Denoting by $\rho$ the distance from $O$ to $m$, again by the sinus formula, we have

$$ \rho = \frac{\sin \gamma_1}{\sin \alpha} l \sim h_K, $$

due to the previous property and the fact that $\alpha \le \alpha + \beta \le \pi - \frac{\alpha_0}{2}$.

We now denote by $\theta_O$ the angle of the half-line containing the segment joining $O$ to $x_j$ and consider the truncated cone:

$$ C_{\rho} = \{(r \cos \theta, r \sin \theta) : 0 < r < \rho \quad \theta_O < \theta < \theta_O + \alpha_1\}. $$

By construction, $C_\rho$ is included into $K$, and degenerates only in the radial direction. Indeed, by setting $C = \{(s \cos \theta, s \sin \theta) : 0 < s < 1 \quad 0 < \theta < \alpha_1\}$ we can introduce the change of variables

$$ F : C \to C_{\rho} : (s, \omega) \to (\rho s, \theta_O + \omega). $$

Then for every $w_h \in V^h$ we see that

$$ \|\chi u_L - w_h\|_{1,\Omega} \ge \|u_L - w_h\|_{1,C_\rho} \ge |u_L \circ F - w_h \circ F|_{1,C}. $$

Since $w_h \circ F$ belongs to $\mathbb{P}_1(C)$, we deduce that

$$ \|\chi u_L - w_h\|_{1,\Omega} \ge |(I - P)(u_L \circ F)|_{1,C}, $$

where $P$ is the projection on $\mathbb{P}_1(C)$ with respect to the inner product of $H^1(C)/\mathbb{P}_0(C)$. Since

$$ u_L(r, \theta) = r^{\frac{1}{2}} \sin \frac{\theta}{2}, $$

we have

$$ \begin{aligned} u_L \circ F(s, \omega) &= \rho^{\frac{1}{2}} s^{\frac{1}{2}} \sin\left(\frac{\theta_O + \omega}{2}\right) \\ &= \rho^{\frac{1}{2}} \left(\sin\frac{\theta_O}{2} S_D(s, \omega) + \cos\frac{\theta_O}{2} S_N(s, \omega)\right), \end{aligned} $$

where we have set

$$ S_N(s, \omega) = s^{\frac{1}{2}} \sin \frac{\omega}{2} \quad \text{and} \quad S_D(s, \omega) = s^{\frac{1}{2}} \cos \frac{\omega}{2}. $$

With these notations, we have

$$ (I - P)(u_L \circ F) = \rho^{\frac{1}{2}} \left(\sin \frac{\theta_O}{2}(I - P)S_D + \cos \frac{\theta_O}{2}(I - P)S_N\right), $$
and therefore

$$ \|\chi u_L - w_h\|_{1,\Omega} \gtrsim \rho^{\frac{1}{2}} |\sin \frac{\theta_O}{2}(I-P)S_D + \cos \frac{\theta_O}{2}(I-P)S_N|_{1,C}. $$

If we can show that

$$ |\sin \frac{\theta_O}{2}(I-P)S_D + \cos \frac{\theta_O}{2}(I-P)S_N|_{1,C} \gtrsim 1, \quad (49) $$

then

$$ \|\chi u_L - w_h\|_{1,\Omega} \gtrsim \rho^{\frac{1}{2}} \gtrsim h^{\frac{1}{2}}. \quad (50) $$

This estimate with (48) in (47) then lead to the conclusion.

It remains to prove (49). For that purpose, we introduce the function $g$ from $[-\frac{\pi}{2}, \frac{\pi}{2}]$ into $\mathbb{R}$ defined by

$$ g(\gamma) = |\sin \gamma(I - P)S_D + \cos \gamma(I - P)S_N|_{1,C}^2. $$

We first notice that $g(\gamma) > 0$ for all $\gamma \in [-\frac{\pi}{2}, \frac{\pi}{2}]$ simply because $\sin \gamma S_D + \cos \gamma S_N$ is not a polynomial. Moreover, $g$ is clearly continuous. Therefore

$$ \min_{\gamma \in [-\frac{\pi}{2}, \frac{\pi}{2}]} g(\gamma) = g(\gamma_0) > 0, $$

for some $\gamma_0 \in [-\frac{\pi}{2}, \frac{\pi}{2}]$. The main point is that this minimum is now independent of $\theta_O$ and therefore the estimate (49) is proved, and the Theorem follows. $\square$

In the same manner for the Lamé system approximated by (12) we recall that

$$ u = u_r + K_I \chi u_I + K_{II} \chi u_{II}, $$

and that our Galerkin solution $u_h \in V_h$ admits the splitting

$$ u_h = u_{rh} + K_{I,h} \chi u_I + K_{II,h} \chi u_{II} $$

where $u_{rh} \in (S^h)^2$.

As before, we can prove the

**Theorem 4** Assume that the triangulation is quasi-uniform in the sense that

$$ h \lesssim h_K \quad \forall K \in \mathcal{T}_h. $$

Then we have

$$ |K_I - K_{I,h}| + |K_{II} - K_{II,h}| \lesssim h^{\frac{1}{2}}. \quad (51) $$

**Proof:** Following [9], we introduce

$$ V_I^h = (S^h)^2 \oplus \operatorname{Span}\{\chi u_{II}\} \text{ and } V_{II}^h = (S^h)^2 \oplus \operatorname{Span}\{\chi u_I\}, $$

and denote by $G_{I,h}u$ and $G_{II,h}u$ the Galerkin approximation of $u$ on $V_I^h$ and $V_{II}^h$ respectively. By Theorem 9.1 of [9] we know that

$$ K_I - K_{I,h} = - \frac{a((I - G_{I,h})u_r, (I - G_{I,h})(\chi u_I))}{a((I - G_{I,h})(\chi u_I), (I - G_{I,h})(\chi u_I))}, $$

$$ K_{II} - K_{II,h} = - \frac{a((I - G_{II,h})u_r, (I - G_{II,h})(\chi u_{II}))}{a((I - G_{II,h})(\chi u_{II}), (I - G_{II,h})(\chi u_{II}))}. $$
Therefore by Cauchy-Schwarz's and Korn's inequalities, we have

$$ 
\begin{aligned}
|K_I - K_{I,h}| & \lesssim \frac{\|(I - G_{I,h})u_r\|_{1,\Omega}}{\|(I - G_{I,h})(\chi u_I)\|_{1,\Omega}}, \\
|K_{II} - K_{II,h}| & \lesssim \frac{\|(I - G_{I,h})u_r\|_{1,\Omega}}{\|(I - G_{I,h})(\chi u_{II})\|_{1,\Omega}}.
\end{aligned}
$$

The remainder of the proof is the same as the one of the previous Theorem. □

<figure><img src="image_14.png" /><figcaption>Figure 12: Numerical convergence of stress intensity factors with respect to the mesh parameter h.</figcaption></figure>

Let us now present some numerical experiments obtained on the Lamé system with the same reference solution as the one on Fig. 7. The implementation of the discrete problem (12) uses Getfem++, the freely available C++ finite element library developed by our team (see [30]). The two stress intensity factors have the same value. The approximation of the stress intensity factor given by $K_{I,h}$ and $K_{II,h}$ in (12) is presented on Fig. 12. Different values of the radiuses $r_0$ and $r_1$ corresponding to the definition of the cut-off function (11) are tested in order to show the crucial influence of the shape of the cut-off function. The optimal rate of convergence is reached in the two cases $(r_0, r_1) = (0.01, 0.4)$ and $(r_0, r_1) = (0.01, 0.2)$. The sharper is the cut-off function, the worst is the approximation of the stress intensity factors. In the case $(r_0, r_1) = (0.2, 0.4)$, the optimal rate of convergence is not reached in the range of values of h studied. This important sensitivity to the shape of the cut-off function suggests to investigate in the future the variant with a pointwise matching [25] or an integral matching [15, 11] which avoid the use of a cut-off function.

The convergence rate obtained is lower than the one obtained by J-integral and interaction integral (see [18, 28] for the principle and [33, 34, 25] for some numerical tests). Such methods require a postprocessing but are superconvergent, for instance in [25] the
order of convergence numerically observed for a $P_1$ finite element method is close to $O(h^2)$.
However, the advantage of the coefficients $K_{I,h}$ and $K_{II,h}$ of (12) is that they are directly
given by the approximation without any postprocessing. Moreover, there is no particular
difficulty when the crack tip is near a boundary of the domain.

Concluding remarks

In this paper we obtained new advances in the analysis of Xfem methods. First, in contrast
with [12] we provide optimal *a priori* error estimates. We also provide an *a priori* error
estimate on the standard Xfem with fixed enrichment area which shows the optimality of
this method. An error estimate on the stress intensity factors computed by the variant
with a cut-off function is also established. We prove that the convergence order is $O(h^{1/2})$
which is confirmed by numerical experiments. This order is rather low compared to the one
obtained with the J-integral (see [18, 33, 34, 28, 25]). However, it permits to have a first
estimate without post-treatment of the solution.

Another interesting perspective is the generalization to 3D cracks where the computation
of the stress intensity factors is more complex. However, a variant with an integral matching
or a cut-off function could be adapted.

References

[1] R. A. Adams. Sobolev Spaces. Academic Press, 1975.

[2] M. Arfaoui, K. Mansouri, A. Rezgui. An Asymptotic finite plane deformation analysis of the elastostatic fields at a notch vertex of an incompressible hyperelastic material. *C.R. Mécanique*, 336(9):737-743, 2008.

[3] E. Béchet, H. Minnebo, N. Moës and B. Burgardt. Improved implementation and robustness study of the X-FEM for stress analysis around cracks. *Int. J. Numer. Meth. Engng.*, 64:1033–1056, 2005.

[4] E. Béchet, N. Moës and B. Wohlmuth. A stable Lagrange multiplier space for stiff interface conditions within the extended finite element method. *To appear*.

[5] H. Blum, M. Dobrowolski. On finite element methods for elliptic equations on domains with corners. *Computing* **28**(1) (1982) 53–63.

[6] S. Bordas, J.G.C. Conley, B. Moran, J. Gray, and E. Nichols. A simulation-based design paradigm for complex cast components. *Engineering with Computers*, 23 Issue 1:25–37, 2007.

[7] S. Bordas and B. Moran. Enriched Finite Elements and Level Sets for Damage Tolerance Assessment of Complex Structures. *Engng. Fract. Mech.*, 73:1176–1201, 2006.

[8] S. Bordas, V.P. Nguyen, C. Dunant, H. Nguyen-Dang, and A. Guidoum. An extended finite element library. *Int. J. Numer. Meth. Engng.*, 71(6):703-732, 2007.

[9] M. Bourlard, M. Dauge and S. Nicaise. Error estimates on the coefficients obtained by the singular function method. *Numer. Funct. Anal. and Optimiz.*, 10(11 & 12):1077–1113, 1989.
[10] M. Bourlard, M. Dauge, M.-S. Lubuma, S. Nicaise. Coefficients of the singularities for elliptic boundary value problems on domains with conical points. III. Finite element methods on polygonal domains. *SIAM J. Numer. Anal.* **29**(1) (1992) 136-155.

[11] E. Chahine. Étude mathématique et numérique de méthodes d'éléments finis étendues pour le calcul en domaines fissurés. PhD thesis. INSA Toulouse, France. 2008.

[12] E. Chahine, P. Laborde and Y. Renard. Crack-tip enrichment in the Xfem method using a cut-off function. *Int. J. Numer. Meth. Engng.*, 75(6):629-646, 2008.

[13] E. Chahine, P. Laborde and Y. Renard. Spider-xfem, an extended finite element variant for partially unknown crack-tip displacement. *European Journal of Computational Mechanics*, 15(5-7):625-636, 2008.

[14] E. Chahine, P. Laborde and Y. Renard. A reduced basis enrichment for the extended finite element method. *Math. Model. Nat. Phenom.*, to appear.

[15] E. Chahine, P. Laborde and Y. Renard. A non-conformal eXtended Finite Element approach: Integral matching XFEM. Submitted.

[16] P. G. Ciarlet. *The Finite Element Method For Elliptic Problems*. North Holland Publishing Company, 1979.

[17] P. Ciarlet, Jr. and J. He. The singular complement method for 2d scalar problems. *C. R. Math. Acad. Sci. Paris*, 336(4):353-358, 2003.

[18] Ph. Destuynder, M. Djaoua. Sur une interprétation mathématique de l'intégrale de Rice en théorie de la rupture fragile. *Math. Meth. in the Appl. Sci.*, 3:70-87, 1981.

[19] M. Dobrowolski. Numerical approximation of elliptic interface and corner problems. Habilitationsschrift, Bonn 1981.

[20] A. Ern, J.L. Guermond. *Eléments finis: théorie, applications, mise en œuvre*. Springer, 2001.

[21] P. Grisvard. Problèmes aux limites dans les polygones - Mode d'emploi. *EDF Bull. Directions Etudes Rech. Sér. C. Math. Inform.* 1, MR 87g:35073, 21-59, 1986.

[22] P. Grisvard. * Singularities in boundary value problems*. Masson, 1992.

[23] A. Hansbo, P. Hansbo. A Finite Element Method for the Simulation of Strong and Weak Discontinuities in Solid Mechanics. *Comput. Methods Appl. Mech. Engrg.*, 193:3523-3540, 2004.

[24] J. Haslinger, Y. Renard. A new fictitious domain approach inspired by the extended finite element method. *Siam J. on Numer. Anal.*, 47(2):1474-1499, 2009.

[25] P. Laborde, Y. Renard, J. Pommier, M. Salaun. High Order Extended Finite Element Method For Cracked Domains. *Int. J. Numer. Meth. Engng.* 64:354-381, 2005.

[26] J.M. Melenk and I. Babuška. The partition of unity finite element method: Basic theory and applications. *Comput. Meths. Appl. Mech. Engrg.*, 139:289-314, 1996.

[27] N. Moës, T. Belytschko. X-FEM: Nouvelles Frontières Pour les Eléments Finis. *Revue européenne des éléments finis* (Calcul des structures GIENS'01), 11:305-318, 01/2002.
[28] N. Moës, J. Dolbow, and T. Belytschko. A finite element method for crack growth without remeshing. *Int. J. Numer. Meth. Engng.*, 46:131–150, 1999.

[29] M.-A. Moussaoui. Sur l'approximation des solutions du problème de Dirichlet dans un ouvert avec coins. In *Singularities and constructive methods for their treatment (Oberwolfach, 1983)*, pages 199–206. Springer, Berlin.

[30] Y. Renard, J. Pommier. An open source generic C++ library for finite element methods. http://home.gna.org/getfem/

[31] F.L. Stazi, E. Budyn, J. Chessa, and T. Belytschko. An extended finite element method with higher-order elements for curved cracks. *Computational Mechanics*, 31:38–48, 2003.

[32] G. Strang, G. Fix. *An Analysis of the Finite Element Method*. Prentice-Hall, Englewood Cliffs, 1973.

[33] B. A. Szabó, Z. Yosibash. Numerical analysis of singularities in two dimensions. II. Computation of generalized flux/stress intensity factors. *Internat. J. Numer. Methods Engrg.* **39**(3) (1996) 409–434.

[34] B. A. Szabó, Z. Yosibash. Superconvergent extraction of flux intensity factors and first derivatives from finite element solutions. *Comput. Methods Appl. Mech. Engrg.* **129**(4) (1996) 349–370.

[35] G. Ventura. On the elimination of quadrature subcells for discontinuous functions in the eXtended Finite-Element Method. *Int. J. Numer. Meth. Engng.*, 66:761–795, 2006.

[36] E. Wyart, D. Coulon, P. Martiny, T. Pardoen, Jean-François Remacle and Frédéric Lani. A substructured FE/XFE method for stress-intensity factors computation in an industrial structure. *Revue Européenne de Mécanique Numérique*, vol. 16/2:199–212, 2007.

[37] E. Wyart, D. Coulon, M. Duflot, T. Pardoen, JF Remacle and F. Lani, A substructured FE/shell/XFE-3D method for crack analysis in thin walled structures. *Int. J. Numer. Meth. Engng.*, Vol. 72:757–779, 2007.

[38] Q.Z. Xiao, B.L. Karihaloo. Improving the accuracy of Xfem crack tip fields using higher order quadrature and statically admissible stress recovery. *Int. J. Numer. Meth. Engng.*, 66:1378–1410, 2006.

[39] M. Zlámal. On the finite element method. *Numer. Math.*, 12:394–409, 1968.