# An energy-conserving scheme for dynamic crack growth with the extended finite element method

Julien Réthoré, Anthony Gravouil, Alain Combescure

► To cite this version:

Julien Réthoré, Anthony Gravouil, Alain Combescure. An energy-conserving scheme for dynamic crack growth with the extended finite element method. International Journal for Numerical Methods in Engineering, 2005, pp.631-659. 10.1002/nme.1283. hal-00373827

HAL Id: hal-00373827

https://hal.science/hal-00373827v1

Submitted on 27 Apr 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 CC BY 4.0 - Attribution - International License
# An energy-conserving scheme for dynamic crack growth using the eXtended finite element method

J. Réthoré, A. Gravouil and A. Combescure¹,*,†

LaMCoS, Laboratoire de Mécanique des Contacts et des Solides, UMR 5514, INSA Lyon, Bat. Jean d'Alembert, 18,20 rue des Sciences 69621 Villeurbanne, France

This paper proposes a generalization of the eXtended finite element method (X-FEM) to model dynamic fracture and time-dependent problems from a more general point of view, and gives a proof of the stability of the numerical scheme in the linear case. First, we study the stability conditions of Newmark-type schemes for problems with evolving discretizations. We prove that the proposed enrichment strategy satisfies these conditions and also ensures energy conservation. Using this approach, as the crack propagates, the enrichment can evolve with no occurrence of instability or uncontrolled energy transfer. Then, we present a technique based on Lagrangian conservation for the estimation of dynamic stress intensity factors for arbitrary 2D cracks. The results presented for several applications are accurate for stationary or moving cracks.

KEY WORDS: dynamic fracture mechanics; numerical stability; energy balance; extended finite element method; dynamic stress intensity factors

## 1. INTRODUCTION

The modelling of crack growth is a problem of great importance in the simulation of failure. Many techniques have been developed to take into account the discontinuity produced by a crack. First, Benzley [1] presented an enrichment to finite elements using asymptotic solutions of static fracture problems. This idea was also developed by Fleming et al. [2] for the element-free Galerkin (EFG) method. The partition of unity method (PUM) introduced by Babuska and Melenk [3] gives a theoretical framework to new techniques such as hp-clouds, generalized finite elements (GEM) or eXtended finite elements (X-FEM). For static fracture problems, Black and Belytschko [4], Moës et al. [5] and Dolbow et al. [6] developed the

Contract/grant sponsor: Office National d'Etudes et de Recherches Aérospatiales, centre de Lille; contract/grant number: F30075

<footnote>*Correspondence to: A. Combescure, LaMCoS, INSA Lyon, Bat. Jean d'Alembert, 18,20 rue des Sciences 69621 Villeurbanne, France.</footnote>

<footnote>†E-mail: alain.combescure@insa-lyon.fr</footnote>
X-FEM method. This technique was used in combination with the level set method by Gravouil et al. [7] for arbitrary 3D cracks. This method was also generalized to arbitrary discontinuities by Belytschko et al. [8]. It was also applied, then modified, in References [9, 10]. In the field of dynamic fracture, EFG and GEM have already been used to simulate the propagation of arbitrary 3D cracks (see References [11, 12]). Chen et al. [13] and, more recently, Belytschko et al. [14] developed a new enrichment technique to avoid the difficulties encountered with the original X-FEM in time-dependent problems [4–6]. In Reference [14], instead of a singular enrichment of the original X-FEM, the authors use a discontinuous enrichment combined with a cohesive law and damage model. The present paper proposes a methodology for using a singular enrichment for time-dependent problems.

A study of energy conserving methods for classical FEM with remeshing was introduced in Reference [15]. The present paper proposes a generalization of the X-FEM to model dynamic fracture and time-dependent problems in a general sense, and gives a proof of the stability of the numerical scheme in the linear case. This stability study is an extension of the balance recovery method presented in Reference [15] for a specific implementation of the X-FEM. Our approach enables us to simulate the dynamic propagation of arbitrary 2D cracks using an enrichment strategy for time-dependent problems. This strategy is justified by a theoretical study of Newmark-type time integrators for problems with evolving discretizations. This study yields stability conditions which also ensure energy conservation. Using this approach, as the crack propagates, the enrichment can evolve with no occurrence of instability or uncontrolled energy transfer.

A technique for the estimation of dynamic stress intensity factors for arbitrary 2D cracks, based on the law of conservation of the Lagrangian, is also presented. This technique can be viewed as an extension of the work by Attigui and Petit [16]. The method uses the concept of virtual crack extension to obtain a domain-independent integral and auxiliary fields for mixed-mode separation. Previous works can be found in References [7, 11, 17–24]. The development presented below takes into account the capabilities of the X-FEM whereby the implicit description of the crack’s geometry enables one to construct a virtual crack extension field tangent to the crack’s faces everywhere. This specificity makes it possible to write a domain-independent integral for an arbitrary crack.

In Section 2, the reference problem is presented and discretized using the X-FEM. In Section 3, the stability of this scheme is studied in the linear case and the discretized energy balance is written for dynamic problems with evolving discretizations. The enrichment technique is also presented in this section. The approach based on the Lagrangian conservation law is developed in Section 4 in order to calculate dynamic energy release rates and stress intensity factors using a domain-independent integral. In the last section, several examples of calculations of stationary or moving cracks are presented to show the accuracy of the method.

## 2. THE REFERENCE PROBLEM IN DYNAMIC FRACTURE MECHANICS

### 2.1. *The continuous problem*

Here, we present a general definition of the dynamic crack growth problem. For this type of calculation, we add an additional unknown $a(t)$, representing a measure of the crack (its length
in two dimensions), to the usual displacement and stress fields. One can express the reference problem as follows:

Reference problem: Given
$$
\begin{cases}
\underline{u}(\underline{x}, 0) \\
\underline{u}(\underline{x}, 0), \text{ find} \\
a(0)
\end{cases}
\quad
\begin{cases}
\underline{u}(\underline{x}, t) \in U \\
\underline{\sigma}(\underline{x}, t) \in S, \underline{x} \in \Omega(t), t \in [0; T] \\
a(t)
\end{cases}
$$
such that:

$$
\bullet \quad \forall \underline{x} \in \partial \Omega_1, \forall t \in [0; T]
$$

$$
\underline{u}(\underline{x}, t) = \underline{u}_d \tag{1}
$$

$$
\bullet \quad \forall t \in [0; T] \quad \forall \underline{v} \in U_0
$$

$$
\int_{\Omega(t)} \rho \ddot{u} \cdot \underline{v} \, d\Omega + \int_{\Omega(t)} \underline{\sigma} : \underline{\varepsilon}(\underline{v}) \, d\Omega = \int_{\Omega(t)} \underline{f}_d \cdot \underline{v} \, d\Omega + \int_{\partial \Omega_2} \underline{F}_d \cdot \underline{v} \, dS \quad (2)
$$

$$
\bullet \quad \forall \underline{x} \in \Omega(t), \forall t \in [0; T]
$$

$$
\underline{\sigma}(\underline{x}, t) = \underline{C}_{\underline{\varepsilon}}(\underline{u}(\underline{x}, t)) \quad (3)
$$

$$
\bullet \quad \forall t \in [0; T]
$$

$$
\begin{align}
\dot{a} &= \dot{a}(\mathcal{G}(t, \dot{a}), \mathcal{G}_c(\dot{a})) \tag{4} \\
\theta &= \theta(K_1(t, \dot{a}), K_2(t, \dot{a})) \nonumber
\end{align}
$$

where $\Omega(t)$ is the domain being considered (Figure 1), $\partial\Omega_1$ the boundary along which dis-
placements $u_d(t)$ are prescribed and $\partial\Omega_2$ the boundary along which forces $F_d(t)$ are prescribed.
$f_d(t)$ are the prescribed volume forces; $u$, $\dot{u}$ and $\ddot{u}$ the displacement, velocity and acceleration
fields; $\underline{\sigma}$ and $\underline{\varepsilon}$ the symmetric stress and strain tensors; $\mathcal{G}$ and $K_i$ the dynamic energy release
rate and stress intensity factor; $\mathbf{C}$ the Hooke's tensor; $\rho$ the mass density; $\mathcal{G}_c$ the critical energy
release rate; $\dot{a}$ and $\theta(t)$ the velocity of the crack's tip and its orientation in the cylindrical
co-ordinate system of the crack's tip; $U$ and $S$ the function spaces associated with the problem,
and $U_0$ the vector space of the virtual fields defined by

$$
U_0 = \{ \underline{v} / \underline{v}(\underline{x}) = 0 \quad \forall \underline{x} \in \partial \Omega_1 + \text{regularity} \} \quad (5)
$$

In numerical simulations, the first step consists in calculating the dynamic energy release
rate and the dynamic stress intensity factors; the second step consists in taking the crack’s
extension and the new geometry into account.

<figure><img src="image_2.png" /><figcaption>Figure 1. Notations used for the different parts of domain Ω.</figcaption></figure>
Remarks

(i) Here, a weak form of the problem is presented. It is interesting to review its underlying local equations in order to pinpoint the boundary conditions being considered:

$$
\begin{align*}
\underline{u} &= \underline{u}_d \quad \text{on } \partial\Omega_1 \\
\underline{\sigma}.\underline{n} &= \underline{F}_d \quad \text{on } \partial\Omega_2 \\
\underline{\sigma}.\underline{n} &= 0 \quad \text{on } \Gamma^+ \text{ and } \Gamma^- \\
\underline{\mathrm{div}}(\underline{\sigma}) + \underline{f}_d &= \rho\ddot{\underline{u}} \quad \text{in } \Omega
\end{align*}
\tag{6}
$$

(ii) Although the domain was defined as time-dependent ($\Omega(t)$), we consider that the crack's length alone ($a(t)$) changes. We make the simplifying hypothesis that the initial and deformed configurations are the same and we assume small perturbations.

(iii) We assume the material to be homogeneous with linear isotropic behaviour. Therefore, the problem falls within the framework of linear dynamic fracture mechanics, in which one can define the dynamic energy release rate $\mathcal{G}$ and the dynamic stress intensity factors $K_i$ (see References [18, 22, 25]).

2.2. The discrete problem

In our approach, we use the X-FEM first introduced in References [4, 5]. In this method, an enrichment to the classical finite element approximation is defined using the partition of unity method developed in Reference [3]. For a static problem (see References [5, 26]), the displacement field can be written using an enriched basis of shape functions:

$$
U = \sum_{i \in \mathcal{N}} N_i(x)U_i + \sum_{i \in \mathcal{N}_{\text{cut}}} N_i(x)H(x)a_i + \sum_{i \in \mathcal{N}_{\text{branch}}} \sum_{\alpha} N_i(x)B_{\alpha}(x)b_{i, \alpha} \quad (7)
$$

where $\mathcal{N}$ is the set of all the nodes in the mesh, $\mathcal{N}_{\text{cut}}$ the set of the nodes which belong to elements completely cut by the crack and $\mathcal{N}_{\text{branch}}$ the set of the nodes which belong to elements partially cut by the crack. $N_i$ are the classical shape functions and $H$ is a discontinuous function whose value is 1 if $x$ is above the crack's surface, -1 otherwise. [$B_\alpha$] are branch functions:

$$
[B_{\alpha}] = \left[ \sqrt{r} \sin\left(\frac{\theta}{2}\right), \sqrt{r} \cos\left(\frac{\theta}{2}\right), \sqrt{r} \sin\left(\frac{\theta}{2}\right) \sin(\theta), \sqrt{r} \cos\left(\frac{\theta}{2}\right) \sin(\theta) \right] \quad (8)
$$

In our approach, all the fields (displacement, velocity and acceleration) are discretized using Equation (7). Consequently, the discrete enriched problem can be written as a classical dynamic problem formulated according to the finite element method: we define the mass matrix and the stiffness matrix in the usual manner. The discrete momentum balance at time $t_{n+1}$ is

$$
M_{n+1} \ddot{U}_{n+1} + K_{n+1} U_{n+1} = F_{n+1} \quad (9)
$$

where $M_{n+1}$ and $K_{n+1}$ are the mass and stiffness matrices at time $t_{n+1}$. Since we intend to represent a moving crack, we have to solve a discretized problem whose number of degrees of freedom (size of the state vector) and crack geometry are evolving. Let us consider a first-order tensor $V$. Its value at time $t_i$ can be written using the set of shape functions of the
<figure><img src="image_3.png" /><figcaption>Figure 2. Calculation strategy.</figcaption></figure>

discretization either at time $t_i$ (denoted $V_i^i$) or at time $t_j$ (denoted $V_i^j$). Then, one can define the projection operator $\Pi_{ij}$ as

$$V_i^i = \Pi_{ij} V_i^j \qquad (10)$$

For the second-order tensors, let us designate the mass and stiffness matrices by

$$M_i^j = \Pi_{ij}^T M_i^i \Pi_{ij} \qquad (11)$$ $$K_i^j = \Pi_{ij}^T K_i^i \Pi_{ij}$$

We will now look at what happens during a time step $\Delta t$ from time $t_n$ to time $t_{n+1}$. As shown in Figure 2, one has to perform two basic steps: change the discretization and solve the momentum balance equation. Changing the discretization (step a) requires the determination of $X_n^{n+1}$ (the state vector at time $t_n$ written on the set of shape functions at time $t_{n+1}$), which will be detailed in Section 3. In step b, the momentum balance equation is solved from time $t_n$ to time $t_{n+1}$ in order to obtain $X_{n+1}^{n+1}$. This equation is discretized using a Newmark-type scheme with constant $\gamma$ and $\beta$ (Equations (13) and (14)). Working with two different discretizations also requires that we rewrite the Newmark scheme:

$$\begin{aligned} U_{n+1}^{n+1} &= U_n^{n+1} + \Delta t \dot{U}_n^{n+1} + \Delta t^2 (\frac{1}{2} - \beta) \ddot{U}_n^{n+1} + \Delta t^2 \beta \ddot{U}_{n+1}^{n+1} \\ \dot{U}_{n+1}^{n+1} &= \dot{U}_n^{n+1} + \Delta t (1 - \gamma) \ddot{U}_n^{n+1} + \Delta t \gamma \ddot{U}_{n+1}^{n+1} \end{aligned} \qquad (12)$$

Now, we can define the discretized reference problem:

Discretized reference problem: Given $$ \begin{cases} U_n^n, \dot{U}_n^n, \ddot{U}_n^n \\ a_n, \dot{a}_n \end{cases} $$ , find $$ \begin{cases} U_n^{n+1}, \dot{U}_n^{n+1}, \ddot{U}_n^{n+1} \\ U_{n+1}^{n+1}, \dot{U}_{n+1}^{n+1}, \ddot{U}_{n+1}^{n+1} \\ a_{n+1}, \dot{a}_{n+1} \end{cases} $$ such that:

$$M_{n+1}^{n+1} \ddot{U}_{n+1}^{n+1} + K_{n+1}^{n+1} U_{n+1}^{n+1} = F_{n+1}^{n+1} \qquad (13)$$

$$U_{n+1}^{n+1} = U_n^{n+1} + \Delta t \dot{U}_n^{n+1} + \Delta t^2 (\frac{1}{2} - \beta) \ddot{U}_n^{n+1} + \Delta t^2 \beta \ddot{U}_{n+1}^{n+1} \qquad (14)$$

$$\begin{aligned} \dot{U}_{n+1}^{n+1} &= \dot{U}_n^{n+1} + \Delta t (1 - \gamma) \ddot{U}_n^{n+1} + \Delta t \gamma \ddot{U}_{n+1}^{n+1} \\ \dot{a}_{n+1} &= \dot{a}_{n+1} (\mathcal{G}(\dot{a}_{n+1}), \mathcal{G}_c(\dot{a}_{n+1})) \end{aligned} \qquad (15)$$

$$\theta_{n+1} = \theta(K_1(\dot{a}_{n+1}), K_2(\dot{a}_{n+1})) \qquad (16)$$
The quantities $U_n^{n+1}$, $\dot{U}_n^{n+1}$ and $\ddot{U}_n^{n+1}$ are additional unknowns of the problem at each time step. We will detail their calculation in Section 3.

### Remark

The discretized problem presented here leads to five equations with eight unknowns, which means that its resolution requires three projections.

## 3. STABILITY ANALYSIS AND DISCRETIZED ENERGY BALANCE FOR DYNAMIC CRACK GROWTH SIMULATION

In this section, we will extend the method presented in Reference [15] to the case where the X-FEM is used. This method enables us to study the stability and the discretized energy balance of dynamic analyses with evolving meshes. When one uses the X-FEM for dynamic crack growth analysis, the mesh does not change but the basis of shape functions varies as the crack grows. Consequently, we can use the same approach and the same notations. First, we develop the stability analysis using the energy method (see References [27, 28]). Then, we study its application to dynamic fracture mechanics and present a method to improve the stability and accuracy for an evolving discretization. Finally, we present the application of this method, called the balance recovery method, in the context of the X-FEM.

### 3.1. Stability study

Let us now consider the case of a calculation with an evolving discretization. One can study the stability of this scheme by extending the method presented in References [27, 28] to this type of calculation. Let us use the notations defined in Section 2 and define the following operators:

$$ \langle V \rangle = \frac{1}{2} (V_{n+1}^{n+1} + V_n^{n+1}) \quad (17) $$

$$ [V] = V_{n+1}^{n+1} - V_n^{n+1} \quad (18) $$

Working with two different discretizations requires that we rewrite the equilibrium equations at times $t_n$ and $t_{n+1}$ on the same discretization. Since the scheme must be stable in the absence of external loads, these are not considered in the stability study (see Reference [29]):

$$ M_n^{n+1} \ddot{U}_n^{n+1} + K_n^{n+1} U_n^{n+1} = 0 \quad (19) $$

$$ M_{n+1}^{n+1} \ddot{U}_{n+1}^{n+1} + K_{n+1}^{n+1} U_{n+1}^{n+1} = 0 $$

Using the approach developed in Reference [27] and observing that:

$$ [M\ddot{U}] = M_{n+1}^{n+1}[\ddot{U}] + [M]\ddot{U}_n^{n+1} $$

$$ [KU] = K_{n+1}^{n+1}[U] + [K]U_n^{n+1} $$
we obtain the stability equation:

$$ \langle \ddot{U} \rangle^T A_{n+1}^{n+1} [\ddot{U}] + \langle \dot{U} \rangle^T K_{n+1}^{n+1} [\dot{U}] = - \left( \gamma - \frac{1}{2} \right) [\ddot{U}]^T A_{n+1}^{n+1} [\ddot{U}] - \frac{1}{\Delta t} [\dot{U}]^T ([M] \dot{U}_n^{n+1} + [K] U_n^{n+1}) \quad (20) $$

with

$$ A_{n+1}^{n+1} = M_{n+1}^{n+1} + \frac{\Delta t^2}{2} (2\beta - \gamma) K_{n+1}^{n+1} $$

On the right-hand side, we recognize the term in $\gamma - \frac{1}{2}$ which governs the stability of the scheme for a calculation with constant discretization. However, the presence of a supplementary term involving $[M]$ and $[K]$ destroys the stability conditions of the standard case. The stability depends on the right-hand side of Equation (20).

The method presented here enables one to perform the stability study when going from $X_n^{n+1}$ to $X_{n+1}^{n+1}$ (step b). In order to treat the problem for the whole time step, one must study the transition $X_n^n, X_{n+1}^{n+1}$ (step a followed by step b). For this purpose, one writes:

$$ \begin{aligned} & \frac{1}{2} [\ddot{U}_{n+1}^{n+1\mathrm{T}} A_{n+1}^{n+1} \ddot{U}_{n+1}^{n+1} + \dot{U}_{n+1}^{n+1\mathrm{T}} K_{n+1}^{n+1} \dot{U}_{n+1}^{n+1}] - \frac{1}{2} [\ddot{U}_n^{n\mathrm{T}} A_n^n \ddot{U}_n^n + \dot{U}_n^{n\mathrm{T}} K_n^n \dot{U}_n^n] \\ & = -\left(\gamma - \frac{1}{2}\right) [\ddot{U}]^T A_{n+1}^{n+1} [\ddot{U}] - \frac{1}{\Delta t} [\dot{U}]^T ([M] \ddot{U}_n^{n+1} + [K] U_n^{n+1}) \end{aligned} \quad (21) $$

## 3.2. The discretized energy balance

With the notations above, the energy balance can be written in the same way as the stability study. We get:

$$ \begin{aligned} & [\dot{U}]^T M_{n+1}^{n+1} \langle \dot{U} \rangle + [U]^T K_{n+1}^{n+1} \langle U \rangle \\ & = \mathcal{W}_{\mathrm{ext}} - \frac{\Delta t^2}{2} (2\beta - \gamma) [\ddot{U}]^T M_{n+1}^{n+1} \langle \ddot{U} \rangle - \frac{\Delta t^2}{2} \left(\gamma - \frac{1}{2}\right) (2\beta - \gamma) [\ddot{U}]^T M_{n+1}^{n+1} [\ddot{U}] \\ & \quad - \left(\gamma - \frac{1}{2}\right) [U]^T K_{n+1}^{n+1} [U] - (\gamma - 1) [U]^T ([M] \ddot{U}_n^{n+1} + [K] U_n^{n+1}) \end{aligned} \quad (22) $$

Again, these are the equations of a problem whose discretization does not evolve, with an additional term involving $[M]$ and $[K]$.

Equation (22) enables one to express an energy balance between $X_n^{n+1}$ and $X_{n+1}^{n+1}$ (step a). Finally, one must calculate the balance between $X_n^n$ and $X_{n+1}^{n+1}$ (step a followed by step b), which is

$$ \frac{1}{2} [\ddot{U}_{n+1}^{n+1\mathrm{T}} M_{n+1}^{n+1} \ddot{U}_{n+1}^{n+1} + U_{n+1}^{n+1\mathrm{T}} K_{n+1}^{n+1} U_{n+1}^{n+1}] - \frac{1}{2} [\dot{U}_n^{n\mathrm{T}} M_n^n \dot{U}_n^n + U_n^{n\mathrm{T}} K_n^n U_n^n] $$
$$
= \mathcal{W}_{\text{ext}} - \frac{\Delta t^2}{2} (2\beta - \gamma)[\ddot{U}]^T M_{n+1}^{n+1} \langle \ddot{U} \rangle - \frac{\Delta t^2}{2} \left(\gamma - \frac{1}{2}\right) (2\beta - \gamma)[\ddot{U}]^T M_{n+1}^{n+1} [\ddot{U}] \\
- \left(\gamma - \frac{1}{2}\right) [U]^T K_{n+1}^{n+1} [U] - (\gamma - 1)[U]^T ([M]\ddot{U}_n^{n+1} + [K]U_n^{n+1}) \quad (23)
$$

$\mathcal{W}_{\text{ext}}$ is the work of the external loads, which can be written as follows:

$$
\mathcal{W}_{\text{ext}} = (\gamma - \frac{1}{2}) [U]^T [F] + [U]^T \langle F \rangle
$$

3.3. Application to fracture mechanics: the balance recovery method

Let us now consider the case of dynamic crack growth analysis. In this particular case, we
know what physical phenomenon causes the change of geometry. Indeed, let us assume a crack
growth $\Delta a$ between times $t_n$ and $t_{n+1}$. The state vector at time $t_n$ can be in equilibrium on this
new geometry only if a distribution of forces $F_+$ is applied along the crack's extension to close
it. By making the residual defined by (24) vanish, the balance recovery method guarantees that
the state vector projected onto the new discretization is in equilibrium (assuming that there are
no external forces):

$$
R = M_{n+1}^{n+1} \ddot{U}_n^{n+1} + K_{n+1}^{n+1} U_n^{n+1} - F_+^{n+1} \quad (24)
$$

Thus, we have the following relations:

$$
[U]^T ([M] \ddot{U}_n^{n+1} + [K] U_n^{n+1}) = [U]^T F_+^{n+1} \quad (25)
$$

$$
[\dot{U}]^T ([M] \ddot{U}_n^{n+1} + [K] U_n^{n+1}) = [\dot{U}]^T F_+^{n+1} \quad (26)
$$

where $[U]^T F_+^{n+1}$ and $[\dot{U}]^T F_+^{n+1}$ correspond to the power and the work of the force distribu-
tion $F_+^{n+1}$, respectively. We also know that these two terms are, from a physical standpoint,
equal to $-2\mathcal{G}\Delta a$ and $-2\mathcal{G}\dot{a}$, respectively ($\mathcal{G}$ being the energy release rate and $\dot{a}$ the crack's
velocity, see Reference [18]). This enables us to verify the relations:

$$
[U]^T ([M] \ddot{U}_n^{n+1} + [K] U_n^{n+1}) = -2 \mathcal{G} \Delta a \quad (27)
$$

$$
[\dot{U}]^T ([M] \ddot{U}_n^{n+1} + [K] U_n^{n+1}) = -2 \mathcal{G} \dot{a} \quad (28)
$$

Let us return to the stability study and energy balance described above. Concerning the energy
balance, we will focus on the influence of the balance recovery on the two stages performed
during the time step $[t_n, t_{n+1}]$. The term $\Delta(X_{n+1}^{n+1}, X_n^n)$ denotes the complete energy balance
defined as the difference between the left-hand side and the right-hand side of Equation (23)
taking Equation (27) into account. $\Delta(X_{n+1}^{n+1}, X_n^{n+1})$ is defined as the difference between the
left-hand side and the right-hand side of Equation (22) taking Equation (27) into account.
Furthermore, $\Delta(X_n^{n+1}, X_n^n)$ is defined by

$$
\Delta(X_n^{n+1}, X_n^n) = \frac{1}{2} [\dot{U}_n^{n+1 \mathrm{T}} M_n^{n+1} \dot{U}_n^{n+1} + U_n^{n+1 \mathrm{T}} K_n^{n+1} U_n^{n+1}] - \frac{1}{2} [\dot{U}_n^{n \mathrm{T}} M_n^n \dot{U}_n^n + U_n^{n \mathrm{T}} K_n^n U_n^n] \quad (29)
$$
One can immediately deduce that

$$ \Delta(X_{n+1}^{n+1}, X_n^n) = \Delta(X_{n+1}^{n+1}, X_n^{n+1}) + \Delta(X_n^{n+1}, X_n^n) \quad (30) $$

One can observe (Equation (30)) that $\Delta(X_{n+1}^{n+1}, X_n^n)$ is the sum of $\Delta(X_n^{n+1}, X_n^n)$, which is the quantity of energy introduced or dissipated during the projection operations (step a), and $\Delta(X_{n+1}^{n+1}, X_n^{n+1})$, which is the energy introduced or dissipated when one goes from $X_n^{n+1}$ to $X_{n+1}^{n+1}$ (step b).

### Remarks

(i) The balance takes into account the energy required to create new crack surfaces ($2\mathcal{G}\Delta a$). This term appears naturally because of its relation to the work of the force distribution $F_+^{n+1}$ (see Reference [18]).

(ii) This method enables us to set $\Delta(X_{n+1}^{n+1}, X_n^{n+1})$ equal to zero a priori and, therefore, to minimize the amount of energy introduced by the successive discretization change operations.

(iii) Using the classical FEM, the discretization change step (remeshing) cannot be controlled and can introduce numerical energy into the simulation.

## 3.4. Use of the balance recovery method in the X-FEM

In the method presented above, we considered dynamic calculations with evolving discretizations. When one uses the X-FEM to simulate dynamic crack propagation (as presented in Reference [15]), the mesh does not change between time $t_n$ and $t_{n+1}$, but new space shape functions are added to simulate the crack's growth. Consequently, the mass and stiffness matrices at times $t_n$ and $t_{n+1}$ are different and one can apply the balance recovery method. If we choose to retain all enrichments from time $t_n$ to time $t_{n+1}$ (see Figure 3), the basis of shape functions at time $t_{n+1}$ is greater than that at time $t_n$ and one can choose the new degrees of freedom to be initially zero. The projected displacement field is

$$ [U_n^{n+1}] = \begin{bmatrix} \frac{U_n^n}{0} \\ \vdots \\ 0 \end{bmatrix} \quad (31) $$

The new stiffness matrix is:

$$ [K_{n+1}^{n+1}] = \begin{bmatrix} K_n^n & \tilde{K}_{n, n+1} \\ \tilde{K}_{n, n+1}^T & \tilde{K}_{n+1, n+1} \end{bmatrix} \quad (32) $$

The force distribution $F_+$ needed to close the crack's extension is replaced by setting new degrees of freedom which represent the crack's extension to zero. Consequently, the state vector calculated at time $t_n$ with the discretization at time $t_{n+1}$ verifies the momentum balance on the geometry at time $t_{n+1}$ (see Equation (24)). Moreover, the total energy associated with
<figure><img src="image_4.png" /><figcaption>Figure 3. Enrichment strategy when the crack propagates.</figcaption></figure>

this state vector is the same for both discretizations because the initialization reduces the
contribution of the new shape functions to zero

$$
U_n^{n+1} K_{n+1}^{n+1} U_n^{n+1} = [U_n^n | \dots | 0] \left[ \frac{K_n^n}{\tilde{K}_{n,n+1}} \middle| \frac{\tilde{K}_{n,n+1}}{\tilde{K}_{n+1,n+1}} \right] \left[ \frac{U_n^n}{0} \middle| \dots \middle| 0 \right] = U_n^{n\mathrm{T}} K_n^n U_n^n \quad (33)
$$

The same proof holds for the kinetic energy term. The equilibrium recovery method for the X-FEM is easy and consists simply in initializing properly the new degrees of freedom which appear as the crack grows. This method guarantees that the numerical scheme is stable and that no numerical energy is introduced when the crack propagates.
### 3.5. Enrichment strategy

The previous section showed how to introduce new degrees of freedom when using the X-FEM. We will now develop the enrichment strategy. As shown in Figure 3, a new singular enrichment is added onto the new set $N_{\text{branch}}^{n+1}$. These new shape functions are associated with the new front, i.e. the local cylindrical co-ordinate of the crack’s tip for these functions is centered on the new crack’s tip. The old singular enrichments are kept and always associated with the old tip. For a discontinuous enrichment, new shape functions are added only on the set $N_{\text{cut}}^{\text{new}} = N_{\text{cut}}^{n+1} \setminus N_{\text{cut}}^n$. If the crack propagates but remains within an element, singular shape functions are added onto the same set of nodes $N_{\text{branch}}^{n+1} = N_{\text{branch}}^n$, but there is no new discontinuous shape function. With such an enrichment strategy, the number of degrees of freedom increases as the crack propagates and the balance recovery method can be used.

As shown in Reference [30], for the partition of unity method with an explicit Newmark-type scheme, the stable time step of the enriched problem is a small fraction of the stable time step of the problem which contains no enriched shape function. For this reason, we now choose the mean acceleration scheme, which is unconditionally stable.

## 4. CALCULATION OF THE DYNAMIC ENERGY RELEASE RATE AND STRESS INTENSITY FACTORS

### 4.1. Interaction integral

In dynamic fracture mechanics, the dynamic energy release rate of a two-dimensional problem can be deduced from the variation of the Lagrangian corresponding to a variation of the crack’s length. In this section, we merely outline the method (see Appendix A for full details). We define the Lagrangian $L$ by the following expression, where $W$ is the strain energy, $T$ the kinetic energy and $W_{\text{ext}}$ the external work:

$$L = T - (W_{\text{int}} - W_{\text{ext}}) \qquad (34)$$

For mixed-mode separation, we will consider a problem with two fields $\underline{u}$ and $\underline{u}^{\text{aux}}$. $\underline{u}$ is the actual displacement field and $\underline{u}^{\text{aux}}$ an auxiliary displacement field, both of which are kinetically admissible, and we can decompose the Lagrangian of the total displacement field $\underline{u}^{\text{tot}} = \underline{u} + \underline{u}^{\text{aux}}$ as

$$L(\underline{u}^{\text{tot}}) = L(\underline{u}) + L(\underline{u}^{\text{aux}}) + L^{\text{int}}(\underline{u}, \underline{u}^{\text{aux}}) \qquad (35)$$

In this equation, $L^{\text{int}}(\underline{u}, \underline{u}^{\text{aux}})$ is the interaction Lagrangian, whose expression can be developed into:

$$L^{\text{int}}(\underline{u}, \underline{u}^{\text{aux}}) = - \int_{\Omega} W^{\text{int}} - T^{\text{int}} d\Omega \qquad (36)$$

where

$$W^{\text{int}} = \underline{\underline{\sigma}} : \nabla \underline{u}^{\text{aux}} = \underline{\underline{\sigma}}^{\text{aux}} : \nabla \underline{u} = [\lambda u_{k,k} \delta_{ij} + \mu (u_{i,j} + u_{j,i})] u_{i,j}^{\text{aux}}$$

$$T^{\text{int}} = \rho \dot{\underline{u}} \cdot \dot{\underline{u}}^{\text{aux}} = \rho u_i u_i^{\text{aux}}$$
<figure><img src="image_5.png" /><figcaption>Figure 4. Different contours for the integration of $\delta L^{\text{int}}$.</figcaption></figure>

### Remarks

(i) In the above equation, , or ,t denote time partial derivation, whereas ,i denotes space partial derivation.

(ii) The external work does not appear in the expression of the interaction Lagrangian because of its linear dependency with respect to the displacement field.

Let us consider the contour and domain represented in Figure 4. Assuming that the Lagrangian conservation law holds, that the virtual crack's extension $\delta l$ is collinear to the tangent to the crack's face and that the crack's faces are traction-free in both the actual and the auxiliary solutions, one obtains the following path-independent integral ($n_j$ being the components of the outward normal vector to the closed domain $A(\Gamma)$ defined by its contour $\Gamma$ enclosing the crack's tip):

$$ \Phi(\Gamma, \underline{\delta l}) = \int_{\Gamma} P_{kj}^{\text{int}} \delta l_k n_j \, ds + \int_{A(\Gamma)} Q_k \delta l_k \, dS \quad (37) $$

where

$$ P_{mj}^{\text{int}} = \left[ (\sigma_{kl}^{\text{aux}} u_{k,l} - \rho \dot{u}_k \dot{u}_k^{\text{aux}}) \delta_{mj} - (\sigma_{ij}^{\text{aux}} u_{i,m} + \sigma_{ij} u_{i,m}^{\text{aux}}) \right] $$

$$ Q_m = \left[ (\sigma_{ij,j}^{\text{aux}} u_{i,m} + \sigma_{ij,j} u_{i,m}^{\text{aux}}) + (\rho \dot{u}_i^{\text{aux}} \dot{u}_{i,m} + \rho \dot{u}_i \dot{u}_{i,m}^{\text{aux}}) \right] $$

In the case of a straight crack, $\dot{a}$ denotes the velocity of the crack's front in the direction $x_1$ defined by the crack. In the above development, ,t denotes the time partial derivative. The front's velocity $\dot{a}$ is obtained by applying the total time derivative in the local co-ordinates of the crack's tip:

$$ \begin{aligned} (\rho \dot{u}_i^{\text{aux}} u_{i,k} + \rho \dot{u}_i u_{i,k}^{\text{aux}})_t &= \frac{\partial}{\partial t} (\rho \dot{u}_i^{\text{aux}} u_{i,k} + \rho \dot{u}_i u_{i,k}^{\text{aux}}) \\ &= \frac{d}{dt} (\rho \dot{u}_i^{\text{aux}} u_{i,k} + \rho \dot{u}_i u_{i,k}^{\text{aux}}) - \dot{a} \frac{\partial}{\partial x_1} (\rho \dot{u}_i^{\text{aux}} u_{i,k} + \rho \dot{u}_i u_{i,k}^{\text{aux}}) \end{aligned} \quad (38) $$
Using Equations (37) and (38), we retrieve the M-integral presented in Reference [16]. Now, if we choose the auxiliary field to be the actual field, we obtain:

$$ \Phi(\Gamma, \delta l_1) = \int_{\Gamma} [(\sigma_{kl}u_{k,l} - \rho\dot{u}_k\dot{u}_k)\delta_{1j} - 2\sigma_{ij}u_{i,1} - 2\dot{a}\rho\dot{u}_i u_{i,1}\delta_{1j}] n_j \delta l_1 \, ds \\ + \int_{A(\Gamma)} 2\frac{d}{dt}[\rho\dot{u}_i\dot{u}_{i,1}]\delta l_1 \, dS \quad (39) $$

With these assumptions, taking the limit of $\Phi$ when $\Gamma$ goes to the crack's tip, the result is $2J$, where $J$ is Rice's $J$-integral (see Reference [31]), which is equal to the dynamic energy release rate $\mathcal{G}$ and, for plane strain conditions, can be related to the dynamic stress intensity factors by the following dynamic equivalent of Irwin's relation (see Reference [25]):

$$ \mathcal{G} = \frac{1 - \nu^2}{E} (f_1(\dot{a})(K_1^{\text{dyn}})^2 + f_2(\dot{a})(K_2^{\text{dyn}})^2) \quad (40) $$

In the above equation, $K_i^{\text{dyn}}$ are the dynamic stress intensity factors defined by References [18, 22] and $f_i$ are universal functions defined (with $c_1$ and $c_2$, respectively, the di-latational and shear wave velocities, $k = 3 - 4\nu$ for plane strain) by

$$ f_i(\dot{a}) = \frac{4\alpha_i(1 - \alpha_j^2)}{(k+1)D(\dot{a})}, \quad (i, j) \in \{1, 2\} \quad (41) $$

$$ \alpha_i = \sqrt{1 - \frac{\dot{a}}{c_i}} \quad (42) $$

$$ D(\dot{a}) = 4\alpha_1\alpha_2 - (1 + \alpha_2^2)^2 \quad (43) $$

The interaction integral is defined as the limit of $\Phi$ when $\Gamma_1$ goes to the crack's tip.

$$ I^{\text{int}} = \lim_{\Gamma_1 \to 0} \Phi(\Gamma_1, q) $$

Choosing the virtual crack's extension field $q$ as defined above enables us to write the interaction integral $I^{\text{int}}$ as a domain-independent integral. If we define $q$ by

$$ q = \begin{cases} \underline{0} & \text{on } \Gamma_2 \text{ and outside } \Gamma_2 \\ \underline{x_1} \delta l & \text{on the crack's tip} \\ \text{tangent to the crack's faces everywhere} & \end{cases} \quad (44) $$

we obtain:

$$ I^{\text{int}} = - \int_{A(\Gamma_2)} q_{k,j} \left[ (\sigma_{ml}^{\text{aux}} u_{m,l} - \rho \dot{u}_l \dot{u}_l^{\text{aux}}) \delta_{kj} - (\sigma_{ij}^{\text{aux}} u_{i,k} + \sigma_{ij} u_{i,k}^{\text{aux}}) \right] dS \\ + \int_{A(\Gamma_2)} q_k \left[ (\sigma_{ij,j}^{\text{aux}} u_{i,k} + \ddot{u}_i u_{i,k}^{\text{aux}}) + (\rho \dot{u}_i^{\text{aux}} \dot{u}_{i,k} + \rho \dot{u}_i \dot{u}_{i,k}^{\text{aux}}) \right] dS \quad (45) $$
## 4.2. Calculation of stress intensity factors

The objective of this section is to develop a domain-independent integral enabling us to estimate stress intensity factors. The expression obtained above does contain a domain-independent integral. Thus, in this subsection, we will focus on the estimation of stress intensity factors based on this integral. In the previous subsection, the virtual crack's extension field was rewritten in order to transform the mixed path, domain-independent integral $\Phi$ into a domain-independent integral $I^{\text{int}}$. In the particular case where the auxiliary fields are taken to be the real fields, we also related $\Phi$ to Rice's $J$-integral. Then, the $J$-integral was related to stress intensity factors using the dynamic energy release rate and the dynamic equivalent of Irwin's relation (Equation (40)). Here, if $K_i^{\text{aux}}$ denotes the Mode-$i$ stress intensity factors corresponding to the auxiliary fields, we can write:

$$I^{\text{int}} = \frac{2(1 - v^2)}{E} (f_1(\dot{a}) K_1^{\text{dyn}} K_1^{\text{aux}} + f_2(\dot{a}) K_2^{\text{dyn}} K_2^{\text{aux}}) \quad (46)$$

Consequently, $K_1^{\text{dyn}}$ and $K_2^{\text{dyn}}$ are calculated by choosing $K_1^{\text{aux}} = 1$, $K_2^{\text{aux}} = 0$ and $K_1^{\text{aux}} = 0$, $K_2^{\text{aux}} = 1$, respectively.

For numerical implementation purposes, the interaction integral is calculated using a $J$-domain (see Reference [24]). This $J$-domain is an additional mesh (see Figure 5) onto which the fields are projected, which is used only for the numerical space integration of the

<figure><img src="image_6.png" /><figcaption>Figure 5. The $J$-domain used to calculate $\mathcal{G}$ and the $K_i^{\text{dyn}}$.</figcaption></figure>
interaction integral. The virtual crack's extension is taken to be $x_1$ at the crack's tip and decreases to $0$ on the boundary of the $J$-domain. For example:

$$ \underline{q} = \left(1 - \left|\frac{x}{d}\right|\right) \left(1 - \left|\frac{y}{d}\right|\right) \underline{x}_1 \qquad (47) $$

where $x$ and $y$ are the local co-ordinates of the crack's tip and $d$ is the size of the $J$-domain. For curved cracks, $\underline{q}$ is curved with respect to the cracks' faces according to the assumptions made in the previous development.

# 5. EXAMPLES

## 5.1. Three-point bending

First, let us consider the problem of a stationary crack in a three-point bending specimen whose geometry is shown in Figure 6. The loading is pure Mode-1 and the dynamic stress intensity factor oscillates between 0 and twice the static value given in Reference [18]. This value is

$$ K_{1s} = \frac{6Sl\sigma_0}{4BW^2} \sqrt{\pi a \Phi\left(\frac{a}{W}\right)} \qquad (48) $$

where $S(0.04 \text{ m})$ is the distance between the supports, $a(0.005 \text{ m})$ the crack's length, $L$ the total length of the beam, $W(0.01 \text{ m})$ its height and $B(1.0 \text{ m})$ its thickness. The loading consists in a constant stress $\sigma_0$ (400 Pa) applied over a length $l$ (0.0025 m). The function $\Phi$ is given by

$$ \Phi\left(\frac{a}{W}\right) = 1.09 - 1.735\frac{a}{W} + 8.2\left(\frac{a}{W}\right)^2 - 14.18\left(\frac{a}{W}\right)^3 + 14.57\left(\frac{a}{W}\right)^4 \qquad (49) $$

The material's properties are $E = 200 \text{ GPa}$, $v = 0.3$ and $\rho = 7860 \text{ kg m}^{-3}$. The numerical results were obtained using the mesh shown in Figure 8, which consists of 20 × 100

<figure><img src="image_7.png" /><figcaption>Figure 6. Three-point bending geometry.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Figure 7. Different domains used for the calculation of the interaction integral.</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>Figure 8. The mesh and its subdivisions for the three-point bending specimen.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Figure 9. $\bar{K}_1$ calculated with different J-domain sizes for the three-point bending specimen.</figcaption></figure>

quadrangles. The time window was 250 µs divided into 200 time steps. In this figure, the subdivided elements (whether cut by the crack or including the crack's tip) are also displayed. Subelements were used only for the numerical integration of singular and discontinuous shape functions (see Reference [5]). The results obtained for different sizes of the J-domain (Figure 7) are shown in Figure 9. $K_1$ is normalized by the static value and oscillates between 0 and 2 with good accuracy. One can also observe that the influence of the size of the J-domain over the calculation of $K_1$ is small.

## 5.2. Semi-infinite plate with a stationary edge crack under mixed-mode loading

A schema of the problem is shown in Figure 10. An analytical solution was obtained by Lee and Freund [32]. To model this theoretical configuration, a 4 m × 6 m mesh with a crack of length $a_0 = 1$ m was used. Because of this mesh's finite dimension, the analytical solution is not valid when the compressive reflected wave arrives at the crack's tip. The prescribed velocity was taken to be $v_0 = 16.5$ m s$^{-1}$ and symmetry boundary conditions were applied at the bottom of the mesh. The material's properties were $E = 200$ GPa, $v = 0.25$ and $\rho = 7833$ kg m$^{-3}$. The results were obtained with 40 time steps for a total computation time of 542 µs (i.e. time $3t_c$ where $t_c = c_d/a_0$). Figure 11 shows the results obtained for different mesh refinements. $\bar{K}$ are the dynamic stress intensity factors normalized by $-Ev_0/2(1-v^2)c_d(\sqrt{\pi a_0})$.
<figure><img src="image_11.png" /><figcaption>Figure 10. Semi-infinite plate with a stationary edge crack under mixed-mode loading.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Figure 11. Numerical and analytical solution $\bar{K}$ for different mesh refinements.</figcaption></figure>

This example was calculated using two different meshes. The coarse mesh used 40 × 60 linear quadrangular elements and the fine mesh 80 × 120 elements. The accuracy of the results was satisfactory with both meshes. The fine mesh seemed to be more sensitive to oscillations associated with the time integration scheme. Nevertheless, we can conclude that the method does not depend on the mesh size.

### 5.3. Semi-infinite crack in an infinite plate subjected to a tensile stress wave

We choose as our third example the case of an infinite plate with a semi-infinite crack because its theoretical solution is known (see Reference [22]). Several authors already treated this example and, therefore, we can also compare our results with theirs. Under the given assumptions (an infinite plate with a semi-infinite crack) and for the geometry described in Figure 12, the analytical solution is valid only at time $t \le 3t_c = 3H/c_1$ when the reflected stress wave reaches
<figure><img src="image_13.png" /><figcaption>Figure 12. Geometry and loading for the infinite plate example.</figcaption></figure>

the crack's tip. When this situation occurs, the Mode-I stress intensity factor for a stationary crack can be written as

$$K_1^{\text{dyn}}(0, t) = \frac{2\sigma_0}{1-\mu} \sqrt{\frac{c_1 t (1-2\mu)}{\pi}} \qquad (50)$$

For a moving crack, we can write:

$$K_1^{\text{dyn}}(\dot{a}, t) = k(\dot{a}) K_1^{\text{dyn}}(0, t) \qquad (51)$$

where $k$ is a universal function of the velocity of the crack's tip $\dot{a}$ which we approximate by

$$k(\dot{a}) = \frac{1 - \dot{a}/c_r}{1 - \dot{a}/2c_r} \qquad (52)$$

Finally, we have:

$$K_1^{\text{dyn}}(\dot{a}, t) = \frac{2\sigma_0}{1-\mu} \sqrt{\frac{c_1 t (1-2\mu)}{\pi}} \frac{1-\dot{a}/c_r}{1-\dot{a}/2c_r} \qquad (53)$$

This theoretical solution will be compared with our numerical results. We will be looking at three cases: the crack does not propagate; the crack starts propagating at a prescribed constant velocity $v_0$ as soon as the wave reaches the crack ($t = t_c$); and the crack propagates at a prescribed constant velocity $v_0$ after time $t = 1.5t_c$. The following numerical results were obtained with the plate dimensions $H=2$ m, $L=10$ m and $l=5$ m and the material properties $E=210$ GPa, $v=0.3$ and $\rho=8000$ kg m$^{-3}$. The tensile stress $\sigma_0$ was 500 MPa. The dimension $d$ of the J-domain was 0.5 m and it contained 36 elements with 16 Gauss points each. $v_0$ was $1500$ m s$^{-1}$ and the stress intensity factors were normalized by $\sigma_0\sqrt{H}$. The solutions were calculated using 40 × 80 quadrangular elements with linear approximation. We stopped the computation at time $T = 3t_c$ when the stress wave reflected from the boundary reached the crack's tip and the analytical solution was no longer valid. For the stationary crack case, $T = 3t_c$ was reached after 200 time steps. The results are plotted in Figure 13. These results match the analytical solution with great accuracy. For the second case, the number of time steps was 20 and the results are not as accurate as in the first case. Oscillations appear at time $t = 2.5t_c$. These oscillations were also observed by the authors of several References [12, 13, 23]. We may assume that although the analytical solution was valid until
<figure><img src="image_14.png" /><figcaption>Figure 13. Numerical and analytical solutions $\bar{K}_1$ for a stationary and a moving crack.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>Figure 14. Numerical and analytical solutions $\bar{K}_1$ for a stationary, then moving crack.</figcaption></figure>

time $T = 3t_c$ the reflected wave reached the J-domain before that and affected the calculation of $K_1^{dyn}$.

In the third case, the crack remained stationary until time $t = 1.5t_c$, then propagated at the prescribed constant velocity $v_0$. The results are presented in Figure 14. The number of time steps was still 20; the slope was reproduced with great accuracy, after which the oscillations previously observed appeared again. In this case, the solution was calculated using a small time step when the crack does not propagate. This can be justified by the fact that one can assume a steady-state situation around the crack's tip. Consequently, the dynamic effects in this region are less than in the rest of the structure. The crack's velocity is approximately ten times less than the dilatational wave speed. Moreover, as shown in Reference [15], dynamic crack propagation can be simulated effectively using a large time step. This technique enabled us to predict the crack's loading prior to the initiation and, consequently, the initiation itself with
<figure><img src="image_16.png" /><figcaption>Figure 15. Definition of the initial geometry.</figcaption></figure>

good accuracy, thus avoiding the difficulties due to the high-frequency numerical solicitations when the crack grows. Let us mention that new, enriched degrees of freedom were added and the stiffness properties of the elements were enriched with new shape functions at each time step. This discontinuity (in time) of the stiffness produced a high-frequency response of the structure around the crack’s tip, which can explain the oscillations observed when a small time step was used. Here, we are reaching the limits of Newmark-type schemes, which do not deal with velocity or displacement discontinuities very well.

### 5.4. Dynamic crack propagation under prescribed displacement

In order to illustrate the effectiveness of the methods presented above from an energy point of view, we calculated the dynamic propagation of a crack in a plate subjected to a prescribed displacement (see Figure 15). The specimen was made of a homogeneous, isotropic, linear elastic material. It was subjected to a prescribed vertical displacement $u_d = 0.0025$ m at its ends. The evolution of this loading followed the Heaviside step function with a rising time of 0.1 ms. The additional data were $c = 0.01$ m, $b = 0.05$ m, $a/b = 0.4$ for the geometry and $E = 186$ GPa, $v = 0.3$, $\rho = 8000$ kg m$^{-3}$, $K_{1c} = 110$ MPa$\sqrt{m}$ for the material’s properties. The criterion used for calculating the velocity of the crack’s tip can be written as follows:

$$ \dot{a} = \begin{cases} 0 & \text{if } K_1 \le K_{1c} \\ c_r \left(1 - \frac{K_{1c}}{K_1}\right) & \text{otherwise} \end{cases} \qquad (54) $$

This example had already been calculated in a classical FEM framework in Reference [15]. In Figure 16, we compare the total energy introduced by the X-FEM scheme when the crack grows (i.e. when the discretization evolves) with classical FEM simulation. Using the FEM, the balance recovery method used to control step b (see Figure 2) was very useful in controlling the energy, but step a still produced some numerical energy due to remeshing and projections. Using the X-FEM and the strategy defined in Section 3, the energy balance was ensured during both steps a and b, as shown clearly in Figure 16. The energy introduced using the X-FEM was non-zero only because the discretized energy balance was calculated using the predicted dynamic energy release rate (i.e. using the interaction integral presented in Section 4), which is quite different from the dynamic energy release rate estimated with the energy definition (see Reference [18]).
<figure><img src="image_17.png" /><figcaption>Figure 16. Evolution of the cumulated unbalanced numerical energy.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>Figure 17. Geometry and loading used to model Kalthoff’s experiments.</figcaption></figure>

5.5. Edge-cracked plate under mixed-mode impulse loading

This example deals with a numerical simulation of Kalthoff’s experiments presented in Reference [33]. By modifying the projectile’s velocity, one observes a transition in the type of failure. At low velocity, i.e. under low strain rate, brittle failure is observed with a global angle of 70°. The experimental configuration modelled is presented in Figure 17. The material’s properties were those of a maraging steel type 18Ni1900: E = 190 GPa, ν = 0.3, ρ = 8000 kg m⁻³, K₁c = 68 MPa√m. For brittle failure, the characteristic velocity of the projectile was v₀ = 16.5 m s⁻¹, L was 0.1 m and the initial crack’s length was l = 0.05 m.

The results were obtained using the maximum hoop stress criterion [34] for the propagation
angle and a constant velocity $\dot{a}=750\,\text{m}s^{-1}$. The overall angle of the final crack's path obtained
was about 65° (Figures 18 and 19) which agrees well with the experimental value of 70°.
These results also agreed with those obtained by Belytschko *et al.* [14] when the authors used
the same criterion.
<figure><img src="image_19.png" /><figcaption>Figure 18. Final crack’s path.</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>Figure 19. Evolution of the crack path for time $t = 25, 50, 75$ and $100$ µs.</figcaption></figure>
## 6. CONCLUSION

The objective of the approach presented above is to simulate dynamic crack propagation using the X-FEM. The main idea consists in using the same singular and discontinuous enrichment as that presented in the scope of linear elastic fracture mechanics for static problems in Reference [5]. First, we generalized the work of Attigui and Petit [16] to the case of propagating arbitrary 2D cracks. Then, we studied Newmark-type schemes for problems with time-dependent discretizations. The stability conditions obtained from this study show how time-dependent problems can be dealt with using the X-FEM in a general sense. The application to dynamic crack propagation leads to an enrichment strategy which ensures numerical stability and energy conservation. The results obtained for several examples were accurate compared to theoretical or experimental data. For stationary cracks under dynamic loading, the method has the same advantages as in the static case (see previous works about the X-FEM). For moving cracks, this approach is proved to be stable and to conserve the energy. The comparison from an energy point of view of the X-FEM and the FEM with remeshing provides a good illustration of the great interest of the approach presented. This method avoids the difficulties encountered with remeshing and projections, which sometimes generate uncontrolled energy transfers. The main effort in our study was devoted to the improvement of the spatial representation of discontinuities. A method was also provided to assess the stability and energy conservation properties of the time integration scheme with time-dependent discretizations. This demonstration was based on a Newmark-type scheme. In dynamic fracture mechanics, a moving crack also implies time discontinuities which are not properly modelled by Newmark integrators. Another time integration scheme could be preferred. The enrichment strategy could be generalized to arbitrary moving discontinuities and the method for arbitrary 2D cracks is expected to be applicable to three-dimensional applications using the level set concept.

# APPENDIX A: INTERACTION INTEGRAL

## A.1. Interaction Lagrangian

Using the notations of Section 4 and considering that each displacement field ($u$, $u^{aux}$ and $u^{tot}$) is kinetically admissible, the variation $\delta L^{int}$ of the interaction Lagrangian for a variation $\delta a_i$ of the crack's length in each direction $x_i$ is zero. Let us define three configurations and their associated descriptions: Lagrangian, Eulerian and arbitrary Lagrangian Eulerian (ALE). The Lagrangian configuration $\Omega_0$ is the initial configuration, whose description (co-ordinates $X_{i0}$) corresponds to the initial location of the crack's tip. The Eulerian configuration $\Omega$ corresponds to an arbitrary extension $\delta l$ of the crack and its description (co-ordinates $X_i$) moves along with the crack's tip. The ALE configuration $\Omega$ and its description are chosen arbitrary between the Lagrangian and Eulerian configurations. Following the notations for the variations of the displacement field between the three configurations defined in Figure A1, one can write:

$$ \delta \tilde{u} = \delta u + \delta \hat{u} \qquad (A1) $$

First, the method is developed in the ALE configuration using the Eulerian description. In Equation (A2), the variation of the interaction Lagrangian is expressed as a function of the
<figure><img src="image_21.png" /><figcaption>Figure A1. Lagrangian, Eulerian and ALE configurations (Ω₀, Ω, Ω̄) and their descriptions (Xᵢ₀, Xᵢ, X̄ᵢ).</figcaption></figure>

true and auxiliary displacement and velocity fields ($l^{\text{int}} = W^{\text{int}} - T^{\text{int}}$):

$$ -\delta L^{\text{int}} = \int_{\Omega} \left[ \frac{\partial l^{\text{int}}}{\partial u_{i,j}} \delta \tilde{u}_{i,j} + \frac{\partial l^{\text{int}}}{\partial u_{i,t}} \delta \tilde{u}_{i,t} + \frac{\partial l^{\text{int}}}{\partial u_{i,j}^{\text{aux}}} \delta \tilde{u}_{i,j}^{\text{aux}} + \frac{\partial l^{\text{int}}}{\partial u_{i,t}^{\text{aux}}} \delta \tilde{u}_{i,t}^{\text{aux}} \right] d\Omega \quad (A2) $$

Using the relation among the three configurations, we obtain:

$$ -\delta L^{\text{int}} = \int_{\Omega} \left[ \frac{\partial l^{\text{int}}}{\partial u_{i,j}} \delta u_{i,j} + \frac{\partial l^{\text{int}}}{\partial u_{i,t}} \delta u_{i,t} + \frac{\partial l^{\text{int}}}{\partial u_{i,j}^{\text{aux}}} \delta u_{i,j}^{\text{aux}} + \frac{\partial l^{\text{int}}}{\partial u_{i,t}^{\text{aux}}} \delta u_{i,t}^{\text{aux}} \right] d\Omega \\ + \int_{\Omega} \left[ \frac{\partial l^{\text{int}}}{\partial u_{i,j}} \delta \hat{u}_{i,j} + \frac{\partial l^{\text{int}}}{\partial u_{i,t}} \delta \hat{u}_{i,t} + \frac{\partial l^{\text{int}}}{\partial u_{i,j}^{\text{aux}}} \delta \hat{u}_{i,j}^{\text{aux}} + \frac{\partial l^{\text{int}}}{\partial u_{i,t}^{\text{aux}}} \delta \hat{u}_{i,t}^{\text{aux}} \right] d\Omega \quad (A3) $$

If $\delta$ denotes the variation between the Eulerian and Lagrangian configurations, we can write in the Eulerian description $\delta = \partial/\partial x_k \delta a_k$.

$$ -\delta L^{\text{int}} = \int_{\Omega} \frac{\partial l^{\text{int}}}{\partial x_k} \delta a_k \, d\Omega + \int_{\Omega} \left[ \frac{\partial l^{\text{int}}}{\partial u_{i,j}} \delta \hat{u}_{i,j} + \frac{\partial l^{\text{int}}}{\partial u_{i,t}} \delta \hat{u}_{i,t} + \frac{\partial l^{\text{int}}}{\partial u_{i,j}^{\text{aux}}} \delta \hat{u}_{i,j}^{\text{aux}} + \frac{\partial l^{\text{int}}}{\partial u_{i,t}^{\text{aux}}} \delta \hat{u}_{i,t}^{\text{aux}} \right] d\Omega \quad (A4) $$
Then, we can use the Gauss-Ostrogradsky theorem and the derivation rule for a product to deduce from Equation (A4):

$$ -\delta L^{\text{int}} = \int_{\partial\Omega} l^{\text{int}} n_k \delta a_k \, dS + \int_{\Omega} \left[ \left( \frac{\partial l^{\text{int}}}{\partial u_{i,j}} \delta \hat{u}_i \right)_{,j} + \left( \frac{\partial l^{\text{int}}}{\partial u_{i,j}^{\text{aux}}} \delta \hat{u}_i^{\text{aux}} \right)_{,j} \right] \delta a_k \, d\Omega \\ - \int_{\Omega} \left[ \left( \frac{\partial l^{\text{int}}}{\partial u_{i,j}} \right)_{,j} \delta \hat{u}_i + \left( \frac{\partial l^{\text{int}}}{\partial u_{i,j}^{\text{aux}}} \right)_{,j} \delta \hat{u}_i^{\text{aux}} \right] \delta a_k \, d\Omega \\ + \int_{\Omega} \left[ \left( \frac{\partial l^{\text{int}}}{\partial u_{i,t}} \delta \hat{u}_i \right)_{,t} + \left( \frac{\partial l^{\text{int}}}{\partial u_{i,t}^{\text{aux}}} \delta \hat{u}_i^{\text{aux}} \right)_{,t} \right] \delta a_k \, d\Omega \\ - \int_{\Omega} \left[ \left( \frac{\partial l^{\text{int}}}{\partial u_{i,t}} \right)_{,t} \delta \hat{u}_i + \left( \frac{\partial l^{\text{int}}}{\partial u_{i,t}^{\text{aux}}} \right)_{,t} \delta \hat{u}_i^{\text{aux}} \right] \delta a_k \, d\Omega \quad (A5) $$

where $\partial\Omega$ is the boundary of the closed domain $\Omega$ and $\bar{n}$ the outward normal to that domain. Now, we define the ALE configuration by choosing the arbitrary configuration to be the Eulerian configuration. Thus, we have $\delta\tilde{u} = 0$ and $\delta\hat{u} = -\delta u$. If we use the Gauss-Ostrogradsky theorem again, we obtain the following expression:

$$ -\delta L^{\text{int}} = \int_{\partial\Omega} l^{\text{int}} n_k \delta a_k \, dS - \int_{\partial\Omega} \left[ \frac{\partial l^{\text{int}}}{\partial u_{i,j}} \delta \hat{u}_i + \frac{\partial l^{\text{int}}}{\partial u_{i,j}^{\text{aux}}} \delta \hat{u}_i^{\text{aux}} \right] n_j \delta a_k \, dS \\ - \int_{\Omega} \left[ \left( \frac{\partial l^{\text{int}}}{\partial u_{i,j}} \right)_{,j} \delta \hat{u}_i + \left( \frac{\partial l^{\text{int}}}{\partial u_{i,j}^{\text{aux}}} \right)_{,j} \delta \hat{u}_i^{\text{aux}} \right] \delta a_k \, d\Omega \\ + \int_{\Omega} \left[ \left( \frac{\partial l^{\text{int}}}{\partial u_{i,t}} \delta \hat{u}_i \right)_{,t} + \left( \frac{\partial l^{\text{int}}}{\partial u_{i,t}^{\text{aux}}} \delta \hat{u}_i^{\text{aux}} \right)_{,t} \right] \delta a_k \, d\Omega \\ - \int_{\Omega} \left[ \left( \frac{\partial l^{\text{int}}}{\partial u_{i,t}} \right)_{,t} \delta \hat{u}_i + \left( \frac{\partial l^{\text{int}}}{\partial u_{i,t}^{\text{aux}}} \right)_{,t} \delta \hat{u}_i^{\text{aux}} \right] \delta a_k \, d\Omega \quad (A6) $$

Then, we can use the definition of $l^{\text{int}}$ to express the derivatives:

$$ \frac{\partial l^{\text{int}}}{\partial u_{i,j}} = \frac{\partial W^{\text{int}}}{\partial u_{i,j}} = \sigma_{ij}^{\text{aux}}, \quad \frac{\partial l^{\text{int}}}{\partial u_{i,j}^{\text{aux}}} = \frac{\partial W^{\text{int}}}{\partial u_{i,j}^{\text{aux}}} = \sigma_{ij} $$

$$ \frac{\partial l^{\text{int}}}{\partial u_{i,t}} = -\frac{\partial T^{\text{int}}}{\partial u_{i,t}} = -\rho u_{i,t}^{\text{aux}}, \quad \frac{\partial l^{\text{int}}}{\partial u_{i,t}^{\text{aux}}} = -\frac{\partial T^{\text{int}}}{\partial u_{i,t}^{\text{aux}}} = -\rho u_{i,t} $$
<figure><img src="image_22.png" /><figcaption>Figure A2. Different contours for the integration of $\delta L^{\text{int}}$.</figcaption></figure>

Finally, we obtain the following expression:

$$ -\delta L^{\text{int}} = \int_{\partial \Omega} \left[ (\sigma_{kl}^{\text{aux}} u_{k,l} - \rho u_{k,t} u_{k,t}^{\text{aux}}) \delta_{mj} - (\sigma_{ij}^{\text{aux}} u_{i,m} + \sigma_{ij} u_{i,m}^{\text{aux}}) \right] n_j \delta a_m \, dS \\ + \int_{\Omega} \left[ (\sigma_{ij,j}^{\text{aux}} u_{i,m} + \sigma_{ij,j} u_{i,m}^{\text{aux}}) + (\rho u_{i,t}^{\text{aux}} u_{i,tm} + \rho u_{i,t} u_{i,tm}^{\text{aux}}) \right] \delta a_m \, d\Omega \quad (A7) $$

Let us define $\underline{P}^{\text{int}}$ and $\underline{Q}$ in order to simplify the expression of $\delta L^{\text{int}}$.

$$ P_{mj}^{\text{int}} = \left[ (\sigma_{kl}^{\text{aux}} u_{k,l} - \rho u_{k,t} u_{k,t}^{\text{aux}}) \delta_{mj} - (\sigma_{ij}^{\text{aux}} u_{i,m} + \sigma_{ij} u_{i,m}^{\text{aux}}) \right] $$

$$ Q_m = \left[ (\sigma_{ij,j}^{\text{aux}} u_{i,m} + \sigma_{ij,j} u_{i,m}^{\text{aux}}) + (\rho u_{i,t}^{\text{aux}} u_{i,tm} + \rho u_{i,t} u_{i,tm}^{\text{aux}}) \right] $$

## A.2. Path-independent integral

The problem is now considered to be two-dimensional, so that $d\Omega = b \, dS$ and $dS = b \, ds$. The closed domain $S$ represented in Figure A2 has an infinitesimal thickness and its contour is $\partial S = \Gamma_1 \cup \Gamma_2 \cup S_+ \cup S_-$. Then, considering the contours and the domain represented in Figure A2 and assuming that $\delta L^{\text{int}}$ is zero, that $\underline{\delta l}$ is collinear to the tangent at the crack's face along $M_+$ or $M_-$ and that the crack's faces are traction-free in the actual and auxiliary solutions, we obtain the following equation (where $n_j$ are the components of the outward normal vector to domain $S$):

$$ \int_{\Gamma_1} P_{kj}^{\text{int}} \delta l_k n_j \, ds - \int_{\Gamma_2} P_{kj}^{\text{int}} \delta l_k n_j \, ds + \int_{A(\Gamma_1)} Q_k \delta l_k \, dS - \int_{A(\Gamma_2)} Q_k \delta l_k \, dS = 0 \quad (A8) $$

This equation shows that $\Phi(\Gamma, \underline{\delta l}) = \int_{\Gamma} P_{kj}^{\text{int}} \delta l_k n_j \, ds + \int_{A(\Gamma)} Q_k \delta l_k \, dS$ does not depend on $\Gamma$.

### Remarks

(i) Since $\Phi$ does not depend on $\Gamma$, using three independent vectors $\underline{\delta l}$, one can write a local equation linking $\underline{P}^{\text{int}}$ to $\underline{Q}$:

$$ \underline{\text{div}} \, \underline{P}^{\text{int}} = \underline{Q} \quad (A9) $$
(ii) Equation (A9) enables us to simplify the final expression and is a consequence of the linear elastic material behaviour.

### A.3. Interaction integral

In this section, the interaction integral is defined as the limit of $\Phi$ when $\Gamma$ goes to the crack's tip.

$$I^{\text{int}} = \lim_{\Gamma_1 \to 0} \Phi(\Gamma_1, q)$$

Then, the virtual crack's extension field $q = \delta l$ is chosen such that we can write the interaction integral $I^{\text{int}}$ as a domain integral.

$$q = \begin{cases} \underline{0} & \text{on } \Gamma_2 \text{ and outside } \Gamma_2 \\ \underline{x_1} \delta l & \text{at the crack's tip} \\ \text{tangent to the crack's faces everywhere} & \end{cases} \qquad (\text{A10})$$

Now, we can replace $\delta l$ by $q$ in the expression of $\Phi$:

$$\Phi(\Gamma_1, q) = \int_{\Gamma_1} P_{kj}^{\text{int}} q_k n_j \, ds + \int_{A(\Gamma_1)} Q_k q_k \, dS \qquad (\text{A11})$$

Considering the properties of $q$, we obtain (with the thickness of domain $S$ no longer being necessarily infinitesimal):

$$\Phi(\Gamma_1, \underline{q}) = - \int_{\partial S} P_{kj}^{\text{int}} q_k n_j \, ds + \int_{A(\Gamma_1)} Q_k q_k \, dS = - \int_{S} (P_{kj}^{\text{int}} q_k)_j \, dS + \int_{A(\Gamma_1)} Q_k q_k \, dS \qquad (\text{A12})$$

If we denote $A = A(\Gamma_2)$, we can write the limit of the first term in the equation above as

$$\left(-\int_{S} (P_{kj}^{\text{int}} q_k)_j \, dS\right) \rightarrow -\int_{A} (P_{kj}^{\text{int}} q_k)_j \, dS \qquad (\text{A13})$$

Now, considering that $q$ is $\underline{0}$ on $\Omega - (S \cup A(\Gamma_1))$, the second term becomes:

$$\int_{A(\Gamma_1)} Q_k q_k \, dS = \int_{\Omega} Q_k q_k \, dS - \int_{S} Q_k q_k \, dS \rightarrow \int_{\Omega} Q_k q_k \, dS - \int_{A} Q_k q_k \, dS \qquad (\text{A14})$$

But $q$ is also $\underline{0}$ on $\Omega - A$; therefore:

$$\int_{A(\Gamma_1)} Q_k q_k \, dS \rightarrow \int_{A} Q_k q_k \, dS - \int_{A} Q_k q_k \, dS = 0 \qquad (\text{A15})$$

Finally, $I^{\text{int}}$ is written as a domain-independent integral:

$$I^{\text{int}} = - \int_{A} (P_{kj}^{\text{int}} q_k)_j \, dS \qquad (\text{A16})$$
Using the definitions of $P^{\text{int}}$ and $Q$, and the local equation (Equation (A9)), the expression of $I^{\text{int}}$ becomes:

$$ I^{\text{int}} = - \int_A q_{k,j} \left[ (\sigma_{ml}^{\text{aux}} u_{m,l} - \rho \dot{u}_l \dot{u}_l^{\text{aux}}) \delta_{kj} - (\sigma_{ij}^{\text{aux}} u_{i,k} + \sigma_{ji} u_{j,k}^{\text{aux}}) \right] dS \\ + \int_A q_k \left[ (\sigma_{ij,j}^{\text{aux}} u_{i,k} + \ddot{u}_i u_{i,k}^{\text{aux}}) + (\rho \dot{u}_i^{\text{aux}} \dot{u}_{i,k} + \rho \dot{u}_i \dot{u}_{i,k}^{\text{aux}}) \right] dS \quad (A17) $$

REFERENCES

1. Benzley SE. Representation of singularities with isoparametric finite elements. *International Journal for Numerical Methods in Engineering* 1974; **8**:537–545.

2. Fleming M, Chu YA, Moran B, Belytschko T. Enriched element free Galerkin methods for singular fields. *International Journal for Numerical Methods in Engineering* 1997; **40**:1483–1504.

3. Babuska I, Melenk JM. The partition of unity method. *International Journal for Numerical Methods in Engineering* 1997; **40**:727–758.

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

5. Moës N, Dolbow J, Belytschko T. A finite element method for crack growth without remeshing. *International Journal for Numerical Methods in Engineering* 1999; **46**(1):133–150.

6. Dolbow J, Moës N, Belytschko T. Discontinuous enrichment in finite elements with a partition of unity method. *Finite Elements in Analysis and Design* 2000; **36**(3):235–260.

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. Belytschko T, Moës N, Usui S, Parimi C. Arbitrary discontinuities in finite elements. *International Journal for Numerical Methods in Engineering* 2001; **50**(4):993–1013.

9. Wells GN, de Borst R, Sluys LJ. A consistent geometrically non-linear approach for delamination. *International Journal for Numerical Methods in Engineering* 2002; **54**:1333–1355.

10. Remmers JJC, de Borst R, Needleman A. A cohesive segments method for the simulation of crack growth. *Computational Mechanics* 2003; **31**:69–77.

11. Krysl P, Belytschko T. Dynamic propagation of arbitrary 3-d cracks. *International Journal for Numerical Methods in Engineering* 1999; **44**(6):767–800.

12. Duarte CA, Hamzeh ON, Liska TJ, Tworzydlo WW. A generalized finite element method for the simulation of three-dimensional dynamic crack propagation. *Computer Methods in Applied Mechanics and Engineering* 2001; **190**:2227–2262.

13. Chen H, Gerlach C, Belytschko T. Dynamic crack growth with X-FEM. *6th USACM 2001*, Dearborn.

14. Belytschko T, Chen H, Jingxiao X, Goangseup Z. Dynamic crack propagation based on loss of hyperbolicity and a new discontinuous enrichment. *International Journal for Numerical Methods in Engineering* 2003; **58**:1873–1905.

15. Rethore J, Gravouil A, Combescure A. A stable numerical scheme for the finite element simulation of dynamic crack propagation with remeshing. *Computer Methods in Applied Mechanics and Engineering* 2004; **193**:4493–4510.

16. Attigui M, Petit C. Mixed-mode separation in dynamic fracture mechanics: new path independent integrals. *International Journal of Fracture* 1997; **84**(1):19–36.

17. Noether E. Invariant variations-problem. *Transport Theory and Statistical Physics* 1918; **1**(3):183–207.

18. Bui HD. *Mécanique de la rupture Fragile*. Masson: New York, 1978.

19. Bui HD. Associated path independent J-integral for separating mixed modes. *Journal of Mechanics and Physics of Solids* 1983; **31**:439–448.

20. Destuynder P. Quelques remarques sur la mécanique de la rupture élastique. *Journal de Mécanique Théorique et Appliquée* 1983; **2**:113–135.
21. Chen FHK, Shield RT. Conservation laws in elasticity of the J-integral type. *Journal of Applied Mechanics and Physics* 1977; **28**:1–22.

22. Freund LB. *Dynamic Fracture Mechanics*. Cambridge Monographs on Mechanics and Applied Mathematics, 1990.

23. Organ D. Numerical solutions to dynamic fracture problems using the element-free Galerkin methods. *Ph.D. Thesis*, Northwestern University, 1996.

24. Goz M, Dolbow J, Moran B. Domain integral formulation for stress intensity factor computation along curved three-dimensional interface cracks. *International Journal of Solids and Structures* 1998; **35**(15):1763–1783.

25. Irwin GR. Analysis of stress and strains near the end of a crack traversing a plate. *Journal of Applied Mechanics* 1957; **24**(3):361–364.

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

27. Hughes TJR, Belytschko T. *Nonlinear Finite Element Analysis*. ICE Division, Zace Services Ltd., 2000.

28. Gravouil A. Méthode multi-échelles en temps et en espace avec décomposition de domaines pour la dynamique non-linéaire des structures. *Ph.D. Thesis*, LMT-Cachan, 2000.

29. Hughes TJR, Liu WK. Implicit-explicit finite element transient analysis: stability theory. *Journal of Applied Mechanics* 1978; **45**:371–374.

30. Gerlach C. Computational methods for the dynamic response of cracked specimens. *Ph.D. Thesis*, 1999.

31. Rice JR. A path independent integral and the approximate analysis of strain concentration by notches and cracks. *Journal of Applied Mechanics* 1968; **35**:379–386.

32. Lee YJ, Freund LB. Fracture initiation due to asymmetric impact loading of an edge cracked plate. *Journal of Applied Mechanics* 1990; **57**:104–111.

33. Kalthoff JF. Modes of dynamic shear failure in solids. *International Journal of Fracture* 2000; **101**:1–31.

34. Broek D. *Elementary Engineering Fracture Mechanics*. Martinus Nijhoff Publishers: Dordrecht, 1982.