# Efficient explicit time stepping for the extended finite element method (X-FEM)

Thomas Menouillard, Julien Réthoré, Alain Combescure, Harridh Bung

► To cite this version:

Thomas Menouillard, Julien Réthoré, Alain Combescure, Harridh Bung. Efficient explicit time stepping for the extended finite element method (X-FEM). International Journal for Numerical Methods in Engineering, 2006, 68, pp.911-939. 10.1002/nme.1718. hal-00938641

HAL Id: hal-00938641

https://hal.science/hal-00938641

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 Attribution 4.0 International License
Efficient explicit time stepping for the eXtended Finite
Element Method (X-FEM)

T. Menouillard¹,², J. Réthoré¹, A. Combescure¹,*,† and H. Bung²

¹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

²Commissariat à l'Energie Atomique, CEA Saclay, DEN/DM2S/SEMT/DYN, F-91191 Gif-sur-Yvette, France

This paper focuses on the introduction of a lumped mass matrix for enriched elements, which enables one to use a pure explicit formulation in X-FEM applications. A proof of stability for the 1D and 2D cases is given. We show that if one uses this technique, the critical time step does not tend to zero as the support of the discontinuity reaches the boundaries of the elements. We also show that the X-FEM element's critical time step is of the same order as that of the corresponding element without extended degrees of freedom.

KEY WORDS: explicit time integration; extended finite element method; dynamic crack propagation

# 1. INTRODUCTION

An important domain of application of the finite element method is the dynamic analysis of impacts, crashes or explosions. When analysing the propagation of waves in a structure, one must use relatively small time steps in order to represent the physics correctly. This is also the case when one uses explicit time integration. Generally, one must go through a very large number of time steps, which sometimes leads to high computation costs. In order to reduce the computation cost of a transient analysis, one can use an explicit time integration scheme (usually the central difference method [1]). Moreover, if a lumping technique is used for the mass matrix, the solution of the problem does not require the resolution of a linear system. One drawback of the central difference method, however, is that its stability is conditional and that the time step must satisfy the Courant's condition. For wave propagation studies, the stability

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

<footnote>Contract/grant sponsor: Commissariat à l'Energie Atomique</footnote>
<figure><img src="image_2.png" /></figure>

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

<figure><img src="image_4.png" /><figcaption>Figure 1. One-dimensional enriched finite element and the equivalent finite element problem.</figcaption></figure>

condition is not really burdensome. This technique is used in a great many commercial codes dedicated to transient analyses.

Unfortunately, a major improvement to the finite element method, the partition of unity method [2], appears to be incompatible with the stability condition of explicit time integrators. The eXtended Finite Element Method [3–5] uses a local partition of unity to incorporate a discontinuity in the displacement or strain field into the interpolation. Thus, the mesh does not need to describe the geometrical support of this discontinuity. This method enables one to simulate the propagation of arbitrary cracks, even in three-dimensional applications [6–8]. It was shown in Reference [9] that the use of a discontinuous enrichment results in a dramatic reduction of the critical time step corresponding to the Courant’s condition when the crack’s surface is close to the nodes of the finite element mesh. This can be easily understood by looking at a one-dimensional problem: let us consider the one-dimensional element enriched with a discontinuous function shown in Figure 1. The equivalent finite element problem consists of two elements with a double node at the location of the discontinuity. For both problems, when the support of the discontinuity approaches Node [1] or Node [2], the critical time step tends to zero. In the equivalent FE problem, the critical time step tends to zero because of the presence of a very large term in the stiffness matrix (~ 1/l) and a very small term in the mass matrix (~ l) of the smallest element. For the X-FEM problem, the critical time step tends to zero because the mass matrix is singular. Several solutions have been proposed in the literature in order to overcome his difficulty. In a previous paper [10], the use of an implicit time integrator was proposed, which is expensive in terms of computation cost. In Reference [11], the authors used an implicit–explicit time integrator in order to treat the region which contains the enriched functions with an unconditionally stable implicit scheme. Using this technique, the mass matrix is not diagonal and the same time step is used both in the implicit and in the explicit regions. In Reference [12], numerical simulations were carried out using a special algorithm which prevents the discontinuity from crossing the element in the vicinity of a node. In addition to this modification of the crack’s path, a large reduction factor was applied to the standard critical time step in order to ensure the stability of the explicit time integrator. Using such a procedure, the time step is often very small in relation to the minimum distance authorized by the algorithm between the discontinuity and the nodes. The solutions mentioned above are not efficient enough for an extensive use of the partition of unity concept in explicit dynamic codes.
In the present paper, a mass matrix lumping technique is proposed for enriched functions. The idea of this lumping technique is based on an exact kinetic energy representation for basic motions. A proof of stability for the 1D and 2D cases is given. We show that if one uses this technique, the critical time step does not tend to zero as the support of the discontinuity reaches the boundaries of the elements. This enables one to use an explicit time integrator with a time step of the order of the critical time step of the problem without enrichment.

The paper is organized as follows: in Section 2, the continuous and discrete formulations of the problem are presented; Section 3 is dedicated to the lumping technique and the proof of stability; numerical examples of dynamic crack propagation are presented in Section 4.

## 2. GOVERNING EQUATIONS AND FORMULATIONS

### 2.1. Continuous formulation

Let us consider the homogeneous body $\Omega$ with a crack $\Gamma$ described in Figure 2. The material has a mass density $\rho$ and is assumed to have linear elastic isotropic behaviour. This behaviour can be described using Lamé's constants $\mu, \lambda$ along with the Young's modulus $E$ and Poisson's ratio $\nu$. The motion of the body is described by the displacement $\mathbf{u}(\mathbf{x}, t)$, $\mathbf{x}$ being the location of a material point and $t$ the time. We assume small perturbations. The body is subjected to prescribed displacements $\mathbf{u}_d$ at the boundary $\partial\Omega_1$ and to external loads $\mathbf{f}_d$ in $\Omega$ and $\mathbf{F}_d$ at the boundary $\partial\Omega_2$. ($\partial\Omega_1$ and $\partial\Omega_2$ are such that $\partial\Omega_1 \cup \partial\Omega_2 = \partial\Omega$ and $\partial\Omega_1 \cap \partial\Omega_2 = \emptyset$.) $\mathbf{n}$ is the outward normal to the material boundary. The crack's faces are traction-free. The strong form of the problem can be written as follows:

$$ \mathbf{u} = \mathbf{u}_d \quad \text{on} \quad \partial\Omega_1 \qquad (1) $$

$$ \underline{\underline{\sigma}}(\mathbf{n}) = \mathbf{F}_d \quad \text{on} \quad \partial\Omega_2 \qquad (2) $$

$$ \underline{\underline{\sigma}}(\mathbf{n}) = \mathbf{0} \quad \text{on} \quad \Gamma_+ \cup \Gamma_- \qquad (3) $$

$$ \mathrm{div}(\underline{\underline{\sigma}}) + \mathbf{f}_d = \rho\ddot{\mathbf{u}} \quad \text{in} \quad \Omega \qquad (4) $$

$$ \underline{\underline{\sigma}} = \mathcal{C} \underline{\underline{\varepsilon}}(\mathbf{u}) = \lambda \operatorname{tr}(\underline{\underline{\varepsilon}}(\mathbf{u})) \underline{\underline{Id}} + 2\mu \underline{\underline{\varepsilon}}(\mathbf{u}) \quad \text{in} \quad \Omega \qquad (5) $$

<figure><img src="image_5.png" /><figcaption>Figure 2. Material body $\Omega$.</figcaption></figure>
where $\underline{\sigma}$ is the Cauchy stress tensor, $\underline{\varepsilon}$ the symmetric first-order strain tensor and $\mathbb{C}$ the constitutive law. $^{\cdot} = \partial/\partial t$ denotes the time derivative. The weak form of these equations becomes:

$$ \int_{\Omega} \rho \underline{\ddot{u}} \cdot \mathbf{v} \, d\Omega + \int_{\Omega} \underline{\underline{\varepsilon}}(\mathbf{u}) : (\mathbb{C} \cdot \underline{\underline{\varepsilon}}(\mathbf{v})) \, d\Omega = \int_{\Omega} \mathbf{f}_d \cdot \mathbf{v} \, d\Omega + \int_{\partial\Omega_2} \mathbf{F}_d \cdot \mathbf{v} \, dS \quad (6) $$

where $\mathbf{u}$ is in the space of the kinematically admissible functions and $\mathbf{v}$ is an arbitrary function kinematically admissible to 0.

## 2.2. Discrete formulation

### 2.2.1. Space discretization.
To capture the discontinuity and the singularity in the strain field at the crack's tip, the partition of unity properties of the finite element shape functions [2] are used to enrich the spatial interpolation as follows:

$$ \mathbf{u}^h = \sum_{i \in \mathcal{N}} N_i(\mathbf{x}) \tilde{\mathbf{U}}_i + \sum_{j=1}^{m} \sum_{i \in \mathcal{N}_j} N_i(\mathbf{x}) \phi_j(\mathbf{x}) \tilde{\mathbf{U}}_{ij} \quad (7) $$

In this equation, $\mathbf{u}^h$ is the approximate displacement field, $\mathcal{N}$ the set of nodes used to discretize $\Omega$, and $N_i$ the finite element shape function associated with node $i$. $\tilde{\mathbf{U}}_i$ are classical degrees of freedom, whereas $\tilde{\mathbf{U}}_{ij}$ are additional degrees of freedom supported by the nodes in the set $\mathcal{N}_j$ associated with the enrichment function $\phi_j$. Letting $\tilde{\mathbf{U}}_{ij} = \mathbf{0}$, $\mathbf{u}^h$ is a standard finite element approximation. Letting $\mathbf{U}_i = \mathbf{0}$, $\tilde{\mathbf{U}}_{ij} = \mathbf{a}$ and $\tilde{\mathbf{U}}_{ij} = \mathbf{0}$ for $j \neq jo$, we have:

$$ \mathbf{u}^h = \sum_{i \in \mathcal{N}} N_i(\mathbf{x}) \phi_{jo}(\mathbf{x}) \mathbf{a} = \phi_{jo}(\mathbf{x}) \mathbf{a} \quad (8) $$

Consequently, the enriched approximation can describe the function $\phi$ exactly. The approximate displacement field $\mathbf{u}^h$ is:

$$ \mathbf{u}^h = N^T \left( \bar{\mathbf{U}} + \sum_{j=1}^{m} \phi_j \tilde{\mathbf{U}}_j \right) = \phi^T \mathbf{U} \quad (9) $$

where $\mathcal{N}$ contains the standard shape functions and $\phi$ the complete basis of shape functions. $\bar{\mathbf{U}}$ and $\tilde{\mathbf{U}}_j$ contain the standard and additional degrees of freedom, respectively, and $\mathbf{U}$ all the degrees of freedom.

Following Reference [3], the basis of enriched functions contains the generalized Heaviside function $H$ and the branch functions:

$$ [B_{\mathrm{el}}] = \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 (10) $$

The nodes enriched with the discontinuous function are those whose support is completely cut by the crack, whereas the nodes enriched with the branch functions are those whose support contains the crack's tip. This enrichment will be denoted $X_{\mathrm{cla}}$.

In the following discussion, we will introduce the enrichment functions by using additional degrees of freedom for the discontinuous function alone. This type of enrichment will be denoted
<figure><img src="image_6.png" /><figcaption>Figure 3. (a) Enrichment $X_{cla}$; and (b) enrichment $X_{dis}$.</figcaption></figure>

$X_{dis}$. The affected nodes are the nodes of the elements cut by the crack (see Figure 3). Using this technique, the element containing the crack’s tip is assumed to be completely cut by the crack, as in References [12, 13]. This leads to some imprecision regarding the location of the asymptotic behaviour of the displacement field, but such errors become negligible as the mesh is refined. In References [11, 14], an enrichment technique which enables the effective positioning of the crack’s tip with a discontinuous enrichment alone is proposed.

Using one of the enrichment techniques presented above, the discrete balance of linear momentum in space is obtained from Equation (6) by choosing **u**, **v** in the complete basis of shape functions $\varphi$:

$$ M\ddot{U} + KU = F \quad (11) $$

**U**, **\dot{U}** and **\ddot{U}** denote, respectively, the displacement, velocity and acceleration vectors discretized on the complete basis of shape functions. **F** is the external force vector. **M** and **K** are the mass and stiffness matrices:

$$ (M)_{ij} = \int_{\Omega} \rho \varphi_i \varphi_j \, d\Omega \quad (12) $$

Then, we introduce the matrix: ($K_{11}$, $K_{12}$, $K_{21}$, $K_{22}$):

$$ (K_{11})_{ij} = \int_{\Omega} (\lambda + \mu) \frac{\partial \varphi_i}{\partial x} \frac{\partial \varphi_j}{\partial x} + \mu \frac{\partial \varphi_i}{\partial y} \frac{\partial \varphi_j}{\partial y} \, d\Omega \quad (13) $$

$$ (K_{12})_{ij} = \int_{\Omega} \lambda \frac{\partial \varphi_i}{\partial y} \frac{\partial \varphi_j}{\partial x} + \mu \frac{\partial \varphi_i}{\partial x} \frac{\partial \varphi_j}{\partial y} \, d\Omega \quad (14) $$

$$ (K_{21})_{ij} = \int_{\Omega} \lambda \frac{\partial \varphi_i}{\partial x} \frac{\partial \varphi_j}{\partial y} + \mu \frac{\partial \varphi_i}{\partial y} \frac{\partial \varphi_j}{\partial x} \, d\Omega \quad (15) $$

$$ (K_{22})_{ij} = \int_{\Omega} (\lambda + \mu) \frac{\partial \varphi_i}{\partial y} \frac{\partial \varphi_j}{\partial y} + \mu \frac{\partial \varphi_i}{\partial x} \frac{\partial \varphi_j}{\partial x} \, d\Omega \quad (16) $$
Thus, the stiffness matrix of an enriched element is

$$ K = \begin{bmatrix} K_{11} & K_{12} \\ K_{21} & K_{22} \end{bmatrix} \qquad (17) $$

When enriched interpolation is used, the consistent mass matrix has standard terms ($\int_{\Omega} \rho N_i N_j d\Omega$), block-diagonal enriched terms ($\int_{\Omega} \rho N_i N_j \phi_k^2 d\Omega$), and coupling terms ($\int_{\Omega} \rho N_i N_j \phi_k \phi_l d\Omega$).

### 2.2.2. Time integration.
The Newmark method is chosen as the time integrator. The updating equations are given by

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

$$ \dot{U}_{n+1} = \dot{U}_n + \Delta t (1 - \gamma) \ddot{U}_n + \Delta t \gamma \ddot{U}_{n+1} \qquad (19) $$

$\beta$ and $\gamma$ are the two parameters of the Newmark scheme. In order to study the stability properties, we follow the energy method from Reference [15]. Thus, we get back to the stability conditions of the Newmark scheme:

$$ \begin{gathered} \frac{1}{2} \le \gamma \le 2\beta \quad \text{unconditionally stable scheme} \\ \frac{1}{2} \le \gamma \text{ and } 2\beta \le \gamma \quad \text{stable scheme if } \Delta t \le \frac{1}{\omega_{\max}\sqrt{(\gamma/2) - \beta}} \end{gathered} \qquad (20) $$

$\omega_{\max}$ is the largest solution of $\det(K - \omega^2 M) = 0$. For the central difference method, $\beta=0$ and $\gamma=0.5$, the problem is: given $U_n$, $\dot{U}_n$ and $\ddot{U}_n$, find $U_{n+1}$, $\dot{U}_{n+1}$ and $\ddot{U}_{n+1}$ such that:

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

$$ M \ddot{U}_{n+1} = F_{n+1} - K U_{n+1} \qquad (22) $$

$$ \dot{U}_{n+1} = \dot{U}_n + \frac{1}{2} \Delta t (\ddot{U}_n + \ddot{U}_{n+1}) \qquad (23) $$

with the stability condition

$$ \Delta t \le \Delta t_c = \frac{2}{\omega_{\max}} \qquad (24) $$

Here, we can see that if the mass matrix is lumped, the resolution of the problem does not require the use of a linear system solver. Since the consistent mass matrix contains standard terms, pure enriched terms and coupling terms, the question becomes: how can one lump such a mass matrix?

## 3. LUMPING TECHNIQUE FOR THE MASS MATRIX

The objective of this section is to construct a lumped mass matrix which enables one to perform explicit dynamic analyses in the X-FEM framework with time steps of reasonable size. A basic
requirement for this lumped mass matrix is that for rigid body motions, the discrete kinetic energy $(\frac{1}{2}\dot{U}^T M \dot{U})$ be exact. We propose to define the diagonal coefficients $m_{diag}$ of the lumped mass matrix corresponding to the enriched degrees of freedom as follows:

$$m_{diag} = \frac{m}{n_{nodes}} \frac{1}{\text{mes}(\Omega_{el})} \int_{\Omega_{el}} \phi^2 d\Omega \qquad (25)$$

where $\Omega_{el}$ is the element being considered, $m$ its mass, mes($\Omega_{el}$) its length (in 1D), area (in 2D) or volume (in 3D), $n_{nodes}$ the number of nodes of $\Omega_{el}$, and $\phi$ the enriched function.

### 3.1. Proof for the 1D case

Following this idea, let us consider a one-dimensional element with two nodes. Each node has ordinary degrees of freedom corresponding to the shape functions $N_1, N_2$ and additional degrees of freedom corresponding to the enriched function $\phi_1$. The approximate displacement is

$$\mathbf{u}^h = N_1 \bar{\mathbf{U}}_1 + N_2 \bar{\mathbf{U}}_2 + N_1 \phi_1 \tilde{\mathbf{U}}_{11} + N_2 \phi_1 \tilde{\mathbf{U}}_{12} \qquad (26)$$

A lumped mass matrix for this element is of the form:

$$M^{\text{lumped}} = \begin{bmatrix} m_1 & 0 & 0 & 0 \\ 0 & m_2 & 0 & 0 \\ 0 & 0 & m_3 & 0 \\ 0 & 0 & 0 & m_4 \end{bmatrix} \qquad (27)$$

We want to find the coefficients $m_i$ such that $T^h = \frac{1}{2}\dot{U}^T M^{\text{lumped}} \dot{U}$ equals $T = \frac{1}{2} \int_{\Omega_{el}} \rho \dot{\mathbf{u}}^2 d\Omega$. First, we consider a motion described by a constant velocity $\dot{\mathbf{u}} = \bar{\mathbf{a}}$. Thus, we set $\dot{\bar{\mathbf{u}}}_i$ to $\bar{\mathbf{a}}$ and $\dot{\bar{\mathbf{U}}}_{i1}$ to $\mathbf{0}$. In this case

$$T^h = \frac{1}{2} \dot{\mathbf{U}}^T M^{\text{lumped}} \dot{\mathbf{U}} = \frac{1}{2} (m_1 \dot{\bar{\mathbf{U}}}_1^2 + m_2 \dot{\bar{\mathbf{U}}}_2^2) = \frac{1}{2} (m_1 + m_2) \bar{\mathbf{a}}^2 \qquad (28)$$

and

$$T = \frac{1}{2} \int_{\Omega_{el}} \rho \dot{\mathbf{u}}^2 d\Omega = \frac{1}{2} m \bar{\mathbf{a}}^2 \qquad (29)$$

where $m$ is the mass of the element. Therefore, if we choose $m_1 = m_2 = m/2$ we get $T^h = T$. This constitutes a practical means of lumping the mass matrix. Now, let us consider a motion described by $\dot{\mathbf{u}} = \bar{\mathbf{a}}\phi_1(\mathbf{x})$ (this represents the separation of the element into two parts). We set $\dot{\bar{\mathbf{U}}}_i$ to $\mathbf{0}$ and $\dot{\bar{\mathbf{U}}}_{i1}$ to $\bar{\mathbf{a}}$, and the energies are

$$T^h = \frac{1}{2} \dot{\mathbf{U}}^T M^{\text{lumped}} \dot{\mathbf{U}} = \frac{1}{2} (m_3 \dot{\bar{\mathbf{U}}}_{11}^2 + m_4 \dot{\bar{\mathbf{U}}}_{12}^2) = \frac{1}{2} (m_3 + m_4) \bar{\mathbf{a}}^2 \qquad (30)$$

and

$$T = \frac{1}{2} \int_{\Omega_{el}} \rho \dot{\mathbf{u}}^2 d\Omega = \frac{1}{2} \rho \bar{\mathbf{a}}^2 \int_{\Omega_{el}} \phi_1^2 d\Omega \qquad (31)$$
Therefore, an exact representation of the kinetic energy ($T^h = T$) is obtained if the terms of the lumped mass matrix are calculated using:

$$m_3 = m_4 = \frac{m}{2} \frac{1}{\text{mes}(\Omega_{\text{el}})} \int_{\Omega_{\text{el}}} \phi_1^2 d\Omega \quad (32)$$

Using this lumped mass matrix, the kinetic energy is represented exactly for all basic motions.

### 3.2. The one-dimensional case

In order to illustrate the properties of the lumping technique presented above, let us determine the critical time step for a one-dimensional element containing a discontinuity. The purpose of this discussion is to study the influence of the location of the discontinuity on the critical time step. Figure 1 presents the 1D element containing a discontinuity. The length of the element is $l$, and the distance between the discontinuity and the left node is $s$. Using the notations of Figure 1, the linear shape functions are

$$N_1(x) = 1 - \frac{x}{l} \quad (33)$$

$$N_2(x) = \frac{x}{l} \quad (34)$$

First, in order to set a reference critical time step, let us determine the critical time step for a standard element (i.e. without discontinuity) of section S, length $l$, Young's modulus $E$ and mass density $\rho$. For this element, the consistent mass matrix $M_{FE}$ and the stiffness matrix $K_{FE}$ are

$$M_{FE} = \rho S l \begin{bmatrix} \frac{1}{3} & \frac{1}{6} \\ \frac{1}{6} & \frac{1}{3} \end{bmatrix} \qquad K_{FE} = \frac{E S}{l} \begin{bmatrix} 1 & -1 \\ -1 & 1 \end{bmatrix} \quad (35)$$

The critical time step given by Equation (24) is

$$\Delta t_c = \frac{2}{\omega_{\max}} = l \sqrt{\frac{\rho}{3E}} \quad (36)$$

Using the mass lumping technique for this element with standard shape functions alone, we obtain:

$$M_{FE}^{\text{lumped}} = \rho S l \begin{bmatrix} \frac{1}{2} & 0 \\ 0 & \frac{1}{2} \end{bmatrix} \quad (37)$$

and the corresponding stable time step is

$$\Delta t_c^{\text{lumped}} = l \sqrt{\frac{\rho}{E}} = \sqrt{3} \Delta t_c^{\text{consistent}} \quad (38)$$
In the following discussion, $\Delta t_c^0$ will denote $\Delta t_c^{\text{lumped}}$, the critical time step for a standard finite element of length $l$ with a lumped mass matrix.

Now, let us deal with the presence of a discontinuity in the displacement field by using the eXtended Finite Element Method to enhance the basis of shape functions. The discontinuity at location $s$ is described with the generalized Heaviside function $H$ centred at $s$:

$$ 
\begin{aligned}
H(x-s) &= -1 & \text{if } x < s \\
H(x-s) &= 1 & \text{if } x > s
\end{aligned} 
\qquad (39) 
$$

The approximate displacement is

$$ \mathbf{u}^h = N_1(x)\bar{\mathbf{U}}_1 + N_2(x)\bar{\mathbf{U}}_2 + H(x-s)N_1(x)\tilde{\mathbf{U}}_{11} + H(x-s)N_2(x)\tilde{\mathbf{U}}_{21} \quad (40) $$

With this enhanced basis of shape functions, the consistent mass matrix and the stiffness matrix are

$$ M_{\text{XFEM}} = \rho S l \begin{bmatrix} \frac{1}{3} & \frac{1}{6} & 2s^2 - 2s + \frac{1}{3} - \frac{2}{3}s^3 & \frac{1}{6} - s^2 + \frac{2}{3}s^3 \\ \frac{1}{6} & \frac{1}{3} & \frac{1}{6} - s^2 + \frac{2}{3}s^3 & \frac{1}{3} - \frac{2}{3}s^3 \\ 2s^2 - 2s + \frac{1}{3} - \frac{2}{3}s^3 & \frac{1}{6} - s^2 + \frac{2}{3}s^3 & \frac{1}{3} & \frac{1}{6} \\ \frac{1}{6} - s^2 + \frac{2}{3}s^3 & \frac{1}{3} - \frac{2}{3}s^3 & \frac{1}{6} & \frac{1}{3} \end{bmatrix} \quad (41) $$

$$ K_{\text{XFEM}} = \frac{ES}{l} \begin{bmatrix} 1 & -1 & 1 - 2s & 2s - 1 \\ -1 & 1 & 2s - 1 & 1 - 2s \\ 1 - 2s & 2s - 1 & 1 & -1 \\ 2s - 1 & 1 - 2s & -1 & 1 \end{bmatrix} \quad (42) $$

Using the lumping technique presented in Section 3.1, the mass matrix becomes

$$ M_{\text{XFEM}}^{\text{lumped}} = \frac{\rho S l}{2} \begin{bmatrix} 1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \end{bmatrix} \quad (43) $$

In both cases (lumped mass and consistent mass), the critical time step given by Equation (24) is a function of the location $s$ of the discontinuity. The results are summarized in Table I. This table gives the critical time steps with consistent mass and lumped mass for three different problems: the first problem is a standard finite element of length $l$; the second is an extended finite element with a discontinuity centred at $s$; the third problem is the equivalent finite element problem FE with two elements (one of length $s$ and the other of length $l-s$, as shown in Figure 1). A plot of the critical time steps for the first and third problems is also shown in
<table><caption>Table I. Normalized critical time steps.</caption><thead><tr><th>Δt<sub>c</sub>/Δt<sub>c</sub><sup>0</sup></th><th>FE length l</th><th>X-FEM</th><th>FE lengths s; l - s</th></tr></thead><tbody><tr><td>Consistent mass</td><td>1/√3</td><td>min(1/√3 l; 1/√3 (1 - s/l))</td><td>min(1/√3 l; 1/√3 (1 - s/l))</td></tr><tr><td>Lumped mass</td><td>1</td><td>min(1/√2 (1/s/l); 1/√2 (1/(1 - (s/l)))</td><td>min(s/l; 1 - s/l)</td></tr></tbody></table>

<figure><img src="image_7.png" /><figcaption>Figure 4. Critical time steps as functions of the location of the discontinuity in a 1D element, with $H$.</figcaption></figure>

<table><caption>Table II. Minimum values for the normalized critical time step.</caption><thead><tr><th>Δt<sub>c</sub>/Δt<sub>c</sub><sup>0</sup></th><th>X-FEM</th><th>FE</th></tr></thead><tbody><tr><td>Consistent mass</td><td>0</td><td>0</td></tr><tr><td>Diagonal mass</td><td>1/√2</td><td>0</td></tr></tbody></table>

Figure 4. As mentioned in Section 1, for the element with enriched shape functions (X-FEM) and consistent mass matrix as well as for the equivalent finite element problem $\tilde{FE}$, the critical time step decreases with the distance between the discontinuity and the closest node. Moreover, the minimum values of $\Delta t_c$ reported in Table II for these problems are zero. Consequently, the
<figure><img src="image_8.png" /><figcaption>Figure 5. Quadrangular element for three different angles $\alpha$: 0, 30 and 45°.</figcaption></figure>

stability of the central difference method cannot be guaranteed for an arbitrary location of the discontinuity.

Using the proposed lumping technique, the critical time step also decreases with the distance between the discontinuity and the closest node, but its minimum value is not zero. Indeed, one can observe in Figure 4 and in Table II that the limit of $\Delta t_c$ when $s$ goes to 0 or $l$ is $\frac{1}{\sqrt{2}} \Delta t_c^0$. This result enables one to perform explicit transient analyses with discontinuous enrichment using a time step of reasonable size regardless of the location of the discontinuity.

### 3.3. The two-dimensional case

In this section, we consider 2D elements cut by a crack. The kinematics can be written as

$$ \mathbf{u}^h = \sum_{i=1}^{n_{el}} N_i(\mathbf{x})(\bar{\mathbf{U}}_i + H(f(\mathbf{x}))\bar{\mathbf{U}}_{i1}) \quad (44) $$

where $f$ is the signed distance to the crack's surface.

We perform the same development as in the one-dimensional case. The interesting point is the evolution of the critical time step as a function of the crack's location for an extended finite element with lumped mass matrix. Figures 5 and 6 show results for three different crack orientations (0, 30 and 45°) for a four- and three-node element, respectively. $\Delta t_c$ is the critical time step of the extended cracked element using the proposed lumping technique. In Figures 5 and 6, the results are normalized by $\Delta t_c^0$, the critical time step of the same element without extended functions. The critical time step is plotted as a function of $S/S_0$, where $S_0$ is the area of the element and $S$ the smaller area obtained when the element is cut by the crack. Similarly, for the three-dimensional case, Figure 7 shows the critical time step as a function of $V/V_0$, where $V_0$ is the volume of the element and $V$ the smaller volume obtained when the element is cut by the crack.
<figure><img src="image_9.png" /><figcaption>Figure 6. Triangular element for three different angles $\alpha$: 0, 30 and 45°.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Figure 7. Three-dimensional cubic element for five different angles $\alpha$: 0, 30, 45, 60 and 90°.</figcaption></figure>

The conclusions are the same as in the one-dimensional case. The proposed lumping technique leads to a reasonable, nonzero minimum value of the critical time step. In practice, according to Figures 5 and 6, explicit transient analysis in the X-FEM framework with discontinuous enrichment can be carried out with a time step of $\frac{2}{3}\Delta t_c^0$.
### 3.4. Implementation of the technique and choice of the discontinuous function

The use of the generalized Heaviside function as the enhanced function is a particular case. Indeed, this function has the property that $H^2$ is a constant function equal to 1. Therefore, in this case, the following two lumping techniques give the same results as that proposed above (provided the coupling terms of the consistent mass matrix are not taken into account):

$$1. (M^{\text{lumped}})_{ii} = \sum_j (M^{\text{consistent}})_{ij}$$

$$2. (M^{\text{lumped}})_{ii} = m \frac{(M^{\text{consistent}})_{ii}}{\sum_j (M^{\text{consistent}})_{jj}}$$

The implementation of the proposed technique using the $H$ function reduces to the use of one of the two standard lumping techniques already mentioned on the pure enriched terms of the consistent mass matrix. This would not be the case for enriched functions whose square is not constant. We show in Appendix A that for other discontinuous enrichments functions, the key factor in getting efficient time steps for transient X-FEM calculations is the choice of the lumping technique.

## 4. APPLICATION TO DYNAMIC CRACK PROPAGATION

In this section, we illustrate the effectiveness of the proposed method. Here, we focus on dynamic crack propagation in the framework of linear elastic fracture mechanics. First, we explain the main principles of the theory. Then, we develop the fracture criteria and present three numerical examples. Through these examples, we aim to validate the lumping technique and to show that the enrichment strategies $X_{\text{cla}}$ and $X_{\text{dis}}$ give similar results. In order to do that, each numerical example will be solved using:

1. implicit time integration with enrichment strategy, $X_{\text{cla}}$;

2. implicit time integration with enrichment strategy, $X_{\text{dis}}$;

3. explicit time integration with enrichment strategy, $X_{\text{dis}}$.

As the crack propagates, the enrichments must evolve to take this propagation into account. The procedure we follow, which was described in Reference [10], consists in keeping the previously enriched degrees of freedom, adding new enriched functions to model the crack's extension and initializing the new degrees of freedom to zero. This procedure is energy-consistent and preserves the stability properties of the time integrator.

### 4.1. Linear elastic fracture mechanics

In linear elastic fracture mechanics, the concept of dynamic stress intensity factor [16] can be used to deal with crack propagation in the static case. Dynamic stress intensity factors are (see Reference [16] for details and reference therein):

$$K_1^{\text{dyn}} = \lim_{r \to 0} \sqrt{2\pi r \sigma_{22}} (\theta = 0)$$

$$K_2^{\text{dyn}} = \lim_{r \to 0} \sqrt{2\pi r \sigma_{12}} (\theta = 0)$$
Since $K_i^{dyn}$ are local quantities defined as limit values, they cannot be estimated directly and accurately. A possible approach was described in Reference [10] (with more details in PhD Thesis [17], in French): a two-field problem (consisting of the actual field $\mathbf{u}$ and an auxiliary field $\mathbf{u}^{aux}$) is constructed for mixed-mode separation, and a domain-independent integral $I^{int}$ is obtained from the Lagrangian conservation law using a virtual crack extension $\mathbf{q}$.

$$ I^{int} = - \int_{S2} [(\underline{\underline{\sigma}}^{aux} : \nabla \mathbf{u} - \rho \dot{\mathbf{u}} \cdot \dot{\mathbf{u}}^{aux}) \mathrm{div}(\mathbf{q}) - (\underline{\underline{\sigma}}^{aux} : (\nabla \mathbf{u} \cdot \nabla \mathbf{q}) + \underline{\underline{\sigma}}^{aux} : (\nabla \mathbf{u}^{aux} \cdot \nabla \mathbf{q}))] dS \\ + \int_{S2} [(\mathbf{div}(\underline{\underline{\sigma}}^{aux}) \cdot \nabla \mathbf{u}(\mathbf{q}) + \rho \ddot{\mathbf{u}} \cdot \nabla \mathbf{u}^{aux}(\mathbf{q})) + (\rho \dot{\mathbf{u}}^{aux} \cdot \nabla \dot{\mathbf{u}}(\mathbf{q}) + \rho \dot{\mathbf{u}} \cdot \nabla \dot{\mathbf{u}}^{aux}(\mathbf{q}))] dS \quad (45) $$

Figure 8 shows the virtual crack extension field defined on two surfaces S1 and S2. $\mathbf{q}$ is particularized as follows:

$\mathbf{q}=0$ outside surface S2

$\|\mathbf{q}\|=1$ on surface S1 \qquad (46)

tangent to the crack's faces everywhere

The use of a virtual crack extension field as defined above (where the size of S1 is chosen to be half the size of S2) decreases the influence of the numerical errors in the estimation of the mechanical fields near the crack's tip on the stress intensity factors.

<figure><img src="image_11.png" /><figcaption>Figure 8. Virtual crack extension field for a curved crack.</figcaption></figure>
If $K_i^{\text{aux}}$ denotes the stress intensity factors of the auxiliary field, we obtain from the dynamic equivalent of Irwin's relation:

$$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 (47)$$

where $f_i$ are universal functions of the speed of the crack's tip $\dot{a}$:

$$f_i(\dot{a}) = \frac{4\alpha_i(1 - \alpha_2^2)}{(\kappa + 1)D(\dot{a})}, \quad i \in \{1, 2\} \quad (48)$$

$$\alpha_i = \sqrt{1 - \left(\frac{\dot{a}}{c_i}\right)^2} \quad (49)$$

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

In these definitions, $c_1$ and $c_2$ are the dilatational and shear wave velocities. The positive root of $D(\dot{a}) = 0$ defines the Rayleigh wave speed $c_r$, which is the theoretical maximum speed for a crack in a homogeneous medium.

Equation (47) shows that the stress intensity factors $K_i^{\text{dyn}}$ can be estimated through an appropriate choice of $\mathbf{u}^{\text{aux}}$ ($K_1^{\text{aux}} = 1, K_2^{\text{aux}} = 0$ for the determination of $K_1^{\text{dyn}}$; $K_1^{\text{aux}} = 0, K_2^{\text{aux}} = 1$ for the determination of $K_2^{\text{dyn}}$). For an accurate numerical evaluation of $K_i^{\text{dyn}}$ in a finite element code, the interaction integral is calculated over a domain, called J-domain, using an additional mesh (the elements of the J-domain are used only for the calculation of the stress intensity factors and are independent of the mesh of the body.) This J-domain follows the crack's front and its orientation as it propagates (see Figure 8). All the numerical results were obtained using 36 square elements with 64 Gauss points for the J-domain.

## 4.2. Fracture criteria

Considering an arbitrary 2D crack under mixed-mode loading, one must check whether the crack is going to propagate and, if necessary, determine the crack's speed and direction. In the following examples, we will consider that the fracture phenomenon is governed by the intensity of the hoop stress ($\sigma_{\theta\theta}$ in the co-ordinate system associated with the crack's tip). Consequently, the crack's direction $\theta_c$ is the direction of the maximum hoop stress, given by

$$\theta_c = 2 \arctan \left[ \frac{1}{4} \left( \frac{K_1^{\text{dyn}}}{K_2^{\text{dyn}}} - \text{sign}(K_2^{\text{dyn}}) \sqrt{8 + \left(\frac{K_1^{\text{dyn}}}{K_2^{\text{dyn}}}\right)^2} \right) \right] \quad (51)$$

Then, the intensity of the loading at the crack's tip is calculated using an equivalent stress intensity factor $K_{\theta\theta}^{\text{dyn}}$:

$$K_{\theta\theta}^{\text{dyn}} = \cos^3\left(\frac{\theta_c}{2}\right) K_1^{\text{dyn}} - \frac{3}{2} \cos\left(\frac{\theta_c}{2}\right) \sin(\theta_c) K_2^{\text{dyn}} \quad (52)$$
Once the crack's direction and the intensity of the stress field in that direction are known, the speed of the crack's tip $\dot{a}$ is given by

$$ \begin{array}{l} \text{if } K_{00}^{\text{dyn}} < K_{1c} \quad \text{then} \quad \dot{a} = 0 \\ \text{else } \dot{a} > 0 \quad \text{and} \quad K_{00}^{\text{dyn}} = K_{1D}(\dot{a}) \end{array} \qquad (53) $$

where $K_{1c}$ is the quasi-static fracture toughness and $K_{1D}(\dot{a})$ the dynamic fracture toughness, which depends on $\dot{a}$ through the following equation (see Reference [18]):

$$ K_{00}^{\text{dyn}} = K_{1D}(\dot{a}) = \frac{K_{1c}}{1 - (\dot{a}/c_r)^m} \qquad (54) $$

Here, we will consider that $m = 1$ and the speed of the crack's tip is obtained from the equation:

$$ \dot{a} = \left( 1 - \frac{K_{1c}}{K_{00}^{\text{dyn}}} \right) c_r \qquad (55) $$

### 4.3. Mode 1

The first example is an infinite plate with a semi-infinite crack. A theoretical solution of this problem for a given crack velocity is known (see Reference [16]). Since this analytical solution was obtained under these assumptions (i.e. infinite plate with a semi-infinite crack and a given speed of the crack's tip) and a numerical one for the geometry described in Figure 9, those could be compared for time $t \le 3t_c = 3h/c_d$ (where $c_d$ is the dilatational wave speed). Beyond that, the reflected stress wave reaches the crack's tip and the analytical solution is no longer valid. As the wave reaches the crack, the mode-1 stress intensity factor for the moving 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 (56) $$

For the moving crack, one can write

$$ K_{1}^{\text{dyn}}(\dot{a}, t) = k(\dot{a}) K_{1}^{\text{dyn}}(0, t) \qquad (57) $$

where $k$ is a universal function of the velocity of the crack's tip $\dot{a}$, which can be approximated by

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

<figure><img src="image_12.png" /><figcaption>Figure 9. Geometry and loading for the infinite plate example.</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>Figure 10. Numerical and analytical solutions $\bar{K}_1$ for a stationary crack.</figcaption></figure>

Finally, one has:

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

The following numerical results were obtained with the plate dimensions $h=2$ m, $L=10$ m, $l=5$ m and the material properties $E=210$ GPa, $v=0.3$, $\rho=8000$ kg m$^{-3}$. The tensile stress $\sigma_0$ was 500 MPa. The stress intensity factors were computed over a J-domain of width 0.1 m. The numerical results are normalized by $\sigma_0\sqrt{h}$.

Since the first enrichment strategy $X_{\text{dis}}$ gives less accurate fields in the vicinity of the crack's tip, the mesh used with that strategy was composed of 60 × 120 linear quadrangular elements. With the strategy $X_{\text{cla}}$, the mesh was 40 × 80 linear quadrangular elements. We will focus on two cases: a stationary crack and a crack propagating at prescribed constant speed $v_0=1500$ m s$^{-1}$ after time $t=1.5t_c$.

### 4.3.1. Stationary crack.
For this case, using the implicit mean acceleration method ($\beta=0.25$ and $\gamma=0.5$), the time $T$ was computed in 100 time steps $\Delta t^{\text{imp}}$. With the central difference method, the size of the time step was $\Delta t^{\text{exp}} = \Delta t^{\text{imp}}/2$. Figure 10 shows the comparison of the results obtained for $\bar{K}_1$. Using the explicit scheme and $X_{\text{dis}}$, the results are less accurate than with the implicit scheme and $X_{\text{cla}}$. The plot of the results obtained with the implicit scheme and $X_{\text{dis}}$ enables one to verify that the oscillations in the Explicit-$X_{\text{dis}}$ solution are not caused by the lumping of the mass matrix but by the loss of accuracy in the mechanical fields due to $X_{\text{dis}}$.
<figure><img src="image_14.png" /><figcaption>Figure 11. Numerical and analytical solutions $\bar{K}_1$ for a stationary, then moving crack.</figcaption></figure>

### 4.3.2. Stationary, then moving crack.
This example has already been treated by many authors: see, e.g. References [8, 10, 11, 19–21]. The numerical solutions described in all these papers presented spurious numerical oscillations after the initiation of the crack’s propagation and during the propagation itself. The phenomenon observed is that the peaks in the oscillations occur just before each instant when the crack’s front crosses the boundary between two elements. In References [10, 21], the use of an implicit scheme or a high-order time integrator (Time eXtended Finite Element) allowed the use of large time steps which filtered the oscillations and provided solutions with very low levels of numerical noise. The results presented in Figure 11 for the Implicit-$X_{\text{cla}}$ case were obtained following this approach. The simulation throughout the time interval was completed in only 20 steps (instead of 100 steps in the previous case). The result is oscillation-free, but the large time discretization prevents the accurate positioning of the crack’s initiation in the time interval. In Figure 11, we also present the result obtained with Explicit-$X_{\text{dis}}$ using 200 time steps. We applied a five-point filter, as in Reference [11], in order to reduce numerical noise. The numerical solution fits the analytical solution quite well, but presents the peaks mentioned above.

### 4.4. Kalthoff’s experiments

This example deals with the numerical simulation of Kalthoff’s experiments of the failure mode transition under pure mode-2 loading [22]. A schematic description of the problem is shown in Figure 12. Kalthoff’s experiments are modelled in plane strain. A plate with two symmetrical edge cracks is impacted by a projectile at speed $V_0$. The two cracks are centred with respect to the specimen’s geometry and their distance corresponds to the diameter of the projectile. At low speed, Kalthoff observed brittle fracture. The two cracks propagate simultaneously with an overall angle from 60 to 70°. If one increases the projectile’s speed, a transition between
<figure><img src="image_15.png" /><figcaption>Figure 12. Geometry and loading for Kalthoff’s experiments.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>Figure 13. Crack’s paths.</figcaption></figure>

brittle fracture and shear band propagation (with a propagation angle of ≈ -10°) occurs. In our case, we chose $V_0 = 20 \text{ m s}^{-1}$ as a typical speed for brittle fracture.

The material’s properties were those of a 18Ni1900 maraging type steel: $E = 190 \text{ GPa}$, $v = 0.3$, $\rho = 8000 \text{ kg m}^{-3}$, $K_{1c} = 68 \text{ MPa} \sqrt{\text{m}}$. The dimensions of the specimen were $L = 0.1 \text{ m}$,
<figure><img src="image_17.png" /><figcaption>Figure 14. Evolution of the crack's length.</figcaption></figure>

$l=0.05$ m, and the initial crack's length was $a=0.05$ m. The mesh was regular and consisted of $80 \times 80$ linear quadrangular elements. For the implicit time integrator, the time step was $\Delta t^{\text{imp}} = 1$ µs. For the explicit scheme, following the results obtained in the previous section, the time step was $\Delta t^{\text{exp}} = \frac{2}{3}\Delta t_c^0 = 0.125$ µs.

Figure 13 shows the crack's paths obtained with three different simulations. These results are very similar and the overall angle agrees with the experimental results. Looking at the details of the crack's paths, we observe the same phenomena as in Reference [12]: the crack starts to propagate with an angle of 65° at $t \approx 28$ µs; then, a small deviation is observed at $t \approx 50$ µs; finally, the crack continues at the initial angle of 65° from $t \approx 65$ µs. These changes in the crack's direction could be related to the propagation of the initial compressive stress wave within the specimen. The initiation occurs at $t \approx 28$ µs, after the wave has travelled from the left side to the right side, has been reflected against the free surface as a tensile stress wave, and finally reaches the crack's tip. Then, this tensile wave travels from the crack's tip to the left side of the specimen and reflects again as a compressive wave which arrives at the moving crack's tip at $t \approx 50$ µs. This compressive wave induces the crack's deviation, but is subsequently reflected against the right side as a tensile wave which reaches the crack's tip at $t \approx 65$ µs, after which the crack continues with the initial angle.

These results enable us to confirm the effectiveness of the method. We can also verify that the crack's evolutions for the three types of simulations carried out are identical by looking at the plot of the evolution of the crack's extension in Figure 14. Although, in this example, the same mesh refinement was used with both enrichment strategies $X_{\text{dis}}$ and $X_{\text{cla}}$, the results are, again, very similar. Figure 15 shows the norm of the displacement field plotted on the deformed mesh for the three types of simulations.
<figure><img src="image_18.png" /><figcaption>Figure 15. Norm of the displacement field plotted on the deformed mesh (×2) at times $t = 40, 60$ and $80$ µs.</figcaption></figure>

### 4.5. *Compact compression specimen (CCS)*

The interest of the CCS experiment in our case is that there is always a mixed-mode loading at the crack's tip. The CCS problem is described schematically in Figure 16. The material's properties are those of PMMA: $E = 5.76$ GPa, $v = 0.42$, $\rho = 1180$ kg m$^{-3}$. The force $P_1(t)$ is
<figure><img src="image_19.png" /><figcaption>Figure 16. Modelling of the CCS: boundary conditions and geometry (specimen thickness: 16.5 mm).</figcaption></figure>

<table><caption>Table III. CCS experiment: time steps and meshes used in the three different methods.</caption><thead><tr><th></th><th>Implicit X<sub>cla</sub></th><th>Implicit X<sub>dis</sub></th><th>Explicit X<sub>dis</sub></th></tr></thead><tbody><tr><td>Mesh (entire structure)</td><td>39 × 60</td><td>39 × 60</td><td>39 × 60</td></tr><tr><td>Time step (10<sup>-6</sup> s)</td><td>1</td><td>1</td><td>0.25</td></tr></tbody></table>

due to an impact at velocity $V_0 = 20 \text{ m s}^{-1}$. The CCS is assumed to be linear elastic. Although the CCS is symmetric, the deformation and, therefore, the crack's path are not. This is due to the non-symmetric loading and boundary conditions.

To validate the lumping technique and to show that the enrichment strategies $X_{cla}$ and $X_{dis}$ give quite similar results, each numerical example is again analysed using Table III:

1. implicit time integration and enrichment strategy, $X_{cla}$;

2. implicit time integration and enrichment strategy, $X_{dis}$;

3. explicit time integration and enrichment strategy, $X_{dis}$.

Figure 17 shows the final deformed shape. Figure 18 shows the calculated crack's paths for the three numerical models. All these results agree with the experiments (see References [23, 24]). We can observe that the results of the explicit simulation with a time step of 0.25 µs are similar to those of the implicit simulation using the same mesh.

Figure 19 shows the evolution of the crack's length as a function of time for the three cases. One can observe that the three methods give the same results. In this case, we showed that using a time step of 0.25 µs enables us to perform a stable calculation of dynamic crack propagation using the eXtended Finite Element Method. Figure 17 shows the crack's path
<figure><img src="image_20.png" /><figcaption>Figure 17. Deformed mesh.</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>Figure 18. Crack's paths, 39 × 60 mesh.</figcaption></figure>
<figure><img src="image_22.png" /><figcaption>Figure 19. Evolution of the crack's length.</figcaption></figure>

superimposed with the nodes: one can observe that although the crack sometimes runs very close to a node, the simulation remains stable.

## 5. CONCLUSION

In this paper, we carried out a numerical study of the evolution of the critical time step for elements with displacement discontinuities. We studied the stability of this simulation in the case where the explicit central difference method is used in combination with an enriched space interpolation (taking advantage of the partition of unity properties of the finite element shape functions). We proposed a lumping technique in order to get a diagonal mass matrix. The presented energy-based technique enabled us to have a lower bound of the critical time step of an arbitrary cracked element of the same order of magnitude as the critical time step of a standard element (regardless of the location of the discontinuity). This very important feature enables us to define the time step in the X-FEM analysis irrespective of the crack's location. Furthermore, we compared our results (explicit time step and lumping technique) with those of implicit calculations for two mixed-mode cases: Kalthoff's experiment and the compact compression specimen experiment. The results of this comparison enable us to conclude that the explicit technique can be used with the eXtended Finite Element Method and a reasonable time step for numerical simulations. This important result demonstrates the potential of explicit time integration in dynamics with crack propagation simulation carried out with X-FEM.
APPENDIX A

First, we introduce the use of the Heaviside function in the 1D case.

Let us go back to the one-dimensional problem treated above and assume that the standard Heaviside function, denoted $\mathcal{H}$, is used instead of $H$:

$$ 
\begin{aligned}
\mathcal{H}(x-s) &= 0 && \text{if } x < s \\
\mathcal{H}(x-s) &= 1 && \text{if } x > s
\end{aligned} 
\qquad (A1) 
$$

The consistent mass matrix and the stiffness matrix are:

$$ M_{\text{XFEM}} = \rho S l \begin{bmatrix} \frac{1}{3} & \frac{1}{6} & s^2 - s + \frac{1}{3} - \frac{1}{3}s^3 & \frac{1}{6} - \frac{s^2}{2} + \frac{1}{3}s^3 \\ \frac{1}{6} & \frac{1}{3} & \frac{1}{6} - \frac{s^2}{2} + \frac{1}{3}s^3 & \frac{1}{3} - \frac{1}{3}s^3 \\ s^2 - s + \frac{1}{3} - \frac{1}{3}s^3 & \frac{1}{6} - \frac{s^2}{2} + \frac{1}{3}s^3 & s^2 - s + \frac{1}{3} - \frac{1}{3}s^3 & \frac{1}{6} - \frac{s^2}{2} + \frac{1}{3}s^3 \\ \frac{1}{6} - \frac{s^2}{2} + \frac{1}{3}s^3 & \frac{1}{3} - \frac{1}{3}s^3 & \frac{1}{6} - \frac{s^2}{2} + \frac{1}{3}s^3 & \frac{1}{3} - \frac{1}{3}s^3 \end{bmatrix} \qquad (A2) $$

$$ K_{\text{XFEM}} = \frac{ES}{l} \begin{bmatrix} 1 & -1 & 1-s & s-1 \\ -1 & 1 & s-1 & 1-s \\ 1-s & s-1 & 1-s & s-1 \\ s-1 & 1-s & s-1 & 1-s \end{bmatrix} \qquad (A3) $$

The first lumping technique leads to the following mass matrix (see Section 3.4):

$$ M_{\mathcal{H}}^{\text{lumped1}} = \rho S l \begin{bmatrix} \frac{1}{2} & 0 & 0 & 0 \\ 0 & \frac{1}{2} & 0 & 0 \\ 0 & 0 & \frac{1}{2} - s + \frac{s^2}{2} & 0 \\ 0 & 0 & 0 & \frac{1}{2} - \frac{s^2}{2} \end{bmatrix} \qquad (A4) $$
The second lumping technique yields:

$$ M_{\mathcal{H}}^{\text{lumped2}} = \rho S l \quad \left[ \begin{array}{cccc} \frac{1}{2} & 0 & 0 & 0 \\ 0 & \frac{1}{2} & 0 & 0 \\ 0 & 0 & \frac{\frac{1}{3} - s + s^2 - \frac{s^3}{3}}{\frac{2}{3} - s + s^2 - \frac{2s^3}{3}} & 0 \\ 0 & 0 & 0 & \frac{\frac{1}{3} - \frac{s^3}{3}}{\frac{2}{3} - s + s^2 - \frac{2s^3}{3}} \end{array} \right] \quad (A5) $$

Finally, our proposed technique yields:

$$ M_{\mathcal{H}}^{\text{lumped}} = \rho S l \quad \left[ \begin{array}{cccc} \frac{1}{2} & 0 & 0 & 0 \\ 0 & \frac{1}{2} & 0 & 0 \\ 0 & 0 & \frac{1}{2}(1-s) & 0 \\ 0 & 0 & 0 & \frac{1}{2}(1-s) \end{array} \right] \quad (A6) $$

Figure A1 shows the critical time steps obtained with these three lumped mass matrices. In this case where the square of the enriched function $\mathcal{H}$ depends on the location of the discontinuity, the results are different. Using standard lumping techniques (such as Techniques 1 and 2), the critical time step appears to depend strongly on $s$. Moreover, the critical time step tends to zero when the discontinuity reaches $s=l$. Conversely, in this configuration, the technique we propose leads to the same critical time step as for a standard element. Otherwise, for $s=0$, the three formulas lead to the same critical time step $\frac{1}{\sqrt{2}}\Delta t_c^0$.

Now, let us present the use of the two functions $H_1$ and $H_2$ in the 1D case (it corresponds to the use of an enrichment with the signed distance function [11, 14]):

$$ H_1(x) = 1 \quad \text{if } x < s \quad (A7) $$

$$ H_1(x) = 0 \quad \text{if } x > s \quad (A8) $$

$$ H_2(x) = 0 \quad \text{if } x < s \quad (A8) $$

$$ H_2(x) = 1 \quad \text{if } x > s \quad (A8) $$
<figure><img src="image_23.png" /><figcaption>Figure A1. Critical time steps as functions of the location of the discontinuity in a 1D element, with $\mathcal{H}$.</figcaption></figure>

where $s$ is the crack's location in the element. Thus, the stiffness matrix is

$$ K_{XFEM} = \frac{ES}{l} \begin{bmatrix} 1 & -1 & s & s-1 \\ -1 & 1 & -s & 1-s \\ s & -s & s & 0 \\ s-1 & 1-s & 0 & 1-s \end{bmatrix} \quad (A9) $$

The first lumping technique provides the following mass matrix:

$$ M_{\mathcal{H}}^{\text{lumped1}} = \rho S l \begin{bmatrix} \frac{1}{2} & 0 & 0 & 0 \\ 0 & \frac{1}{2} & 0 & 0 \\ 0 & 0 & \frac{s}{2} & 0 \\ 0 & 0 & 0 & \frac{1-s}{2} \end{bmatrix} \quad (A10) $$
The time step can be calculated as a function of the parameter s:

$$ \Delta t_c(s) = \sqrt{\frac{2}{3}} \Delta t_c^0 \qquad (A11) $$

We obtain the same time step with the following enrichment functions:

$$ H_1(x) = 0 \quad \text{if } x < s $$

$$ H_1(x) = 1 \quad \text{if } x > s $$

(A12)

$$ H_2(x) = -1 \quad \text{if } x < s $$

$$ H_2(x) = 0 \quad \text{if } x > s $$

ACKNOWLEDGEMENTS

The support of the French Commissariat à l'Energie Atomique is gratefully acknowledged.

REFERENCES

1. Newmark N. A method of computation for structural dynamics. *Journal of the Engineering Mechanics Division, Proceedings of the ASCE (Paper No. 2034)* 1959; **85**(EM3):67-94.

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

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

4. Dolbow J, Moës N, Belytschko T. Discontinuous enrichment in finite elements with a partition of unity method. *International Journal for Numerical Methods in Engineering* 2000; **36**(3):235-260.

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

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

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

8. Duarte C, Hamzeh O, Liska T, Tworzydlo W. 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.

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

10. Réthoré J, Gravouil A, Combescure A. An energy conserving scheme for dynamic crack growth with the extended finite element method. *International Journal for Numerical Methods in Engineering* 2005; **63**(5):631-659.

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

12. de Borst R, Remmers J, Needleman A. Mesh-independent discrete numerical representations of cohesive-zone models. *Engineering Fracture Mechanics* 2006; **73**(2):160-177.

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

14. Zi G, Belytschko T. New crack-tip elements for xfem and applications to cohesive cracks. *International Journal for Numerical Methods in Engineering* 2003; **57**(15):2221-2240.
15. Hughes T, Belytschko T. *Nonlinear Finite Element Analysis*. ICE Division, Zace Services Ltd., 2000.

16. Freund L. *Dynamic Fracture Mechanics*. Cambridge Monographs on Mechanics and Applied Mathematics, 1990.

17. Réthoré J. Méthode élément finis étendus en espace et temps: application à la propagation dynamique des fissures. *Ph.D. Thesis*, Institut National de Sciences Appliquées de Lyon, 2005.

18. Kanninen M, Popelar C. *Advanced Fracture Mechanics*. Oxford University Press: New York, 1985.

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

20. Belytschko T, Chen H. Singular enrichment finite element method for elastodynamic crack propagation. *International Journal of Computational Methods* 2001; **1**(1):1–15.

21. Réthoré J, Gravouil A, Combescure A. A combined space time extended finite element method. *International Journal for Numerical Methods in Engineering* 2005; **64**(2):304–321.

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

23. Rittel D, Maigre H. An investigation of dynamic crack initiation in pmma. *Mechanics of Materials* 1996; **23**:229–239.

24. Rittel D, Maigre H. A study of mixed-mode dynamic crack initiation in pmma. *Mechanics Research Communications* 1996; **23**:475–481.