# Mass lumping strategies for x-fem explicit dynamics: application to crack propagation

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

► To cite this version:

Thomas Menouillard, Julien Réthoré, Nicolas Moes, Alain Combescure, Harridh Bung. Mass lumping strategies for x-fem explicit dynamics: application to crack propagation. International Journal for Numerical Methods in Engineering, 2008, 74, pp.447-474. 10.1002/nme.2180. hal-00938642

HAL Id: hal-00938642

https://hal.science/hal-00938642

Submitted on 26 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
Mass lumping strategies for X-FEM explicit dynamics:
Application to crack propagation

T. Menouillard¹,², J. Réthoré¹, N. Moës³, A. Combescure¹,*,† and H. Bung²

¹LaMCoS, INSA-Lyon, CNRS UMR5259, 18-20 rue des Sciences, F69621 Villeurbanne Cedex, France

²CEA Saclay, DEN/DM2S/SEMT/DYN, 91191 Gif-sur-Yvette, France

³GEM, Ecole Centrale Nantes, CNRS UMR6183, 1 rue de la Noë, 44321 Nantes, France

This paper deals with the numerical modelling of cracks in the dynamic case using the extended finite element method. More precisely, we are interested in explicit algorithms. We prove that by using a specific lumping technique, the critical time step is exactly the same as if no crack were present. This somewhat improves a previous result for which the critical time step was reduced by a factor of square root of 2 from the case with no crack. The new lumping technique is obtained by using a lumping strategy initially developed to handle elements containing voids. To be precise, the results obtained are valid only when the crack is modelled by the Heaviside enrichment. Note also that the resulting lumped matrix is block diagonal (blocks of size 2 × 2). For constant strain elements (linear simplex elements) the critical time step is not modified when the element is cut. Thanks to the lumped mass matrix, the critical time step never tends to zero. Moreover, the lumping techniques conserve kinetic energy for rigid motions. In addition, tensile stress waves do not propagate through the discontinuity. Hence, the lumping techniques create neither error on kinetic energy conservation for rigid motions nor wave propagation through the crack. Both these techniques will be used in a numerical experiment.

KEY WORDS: extended finite element method; explicit dynamics; critical time step; dynamic crack propagation; lumped mass matrix; displacement discontinuity

# 1. INTRODUCTION

The extended finite element method (X-FEM) allows one to introduce a crack within an existing mesh without the need to modify it. Discontinuous enrichments are introduced within the elements cut by the crack using the partition of unity technique [1]. According to the original papers on the subject, the enrichment is composed of a tip enrichment [2, 3] and a Heaviside enrichment [4]

<footnote>*Correspondence to: A. Combescure, LaMCoS, INSA-Lyon, CNRS UMR5259, 18-20 rue des Sciences, F69621 Villeurbanne Cedex, France.</footnote>

<footnote>†E-mail: alain.combescure@insa-lyon.fr</footnote>
away from the crack tip. Initially developed in the 2D setting, the method was then extended to 3D in [5]. At about the same time, two level sets were introduced to conveniently store the crack as two finite element fields [6, 7]. This level set representation of cracks was also found to be extremely convenient to handle crack growth [7–9] (without any remeshing).

The implementation of the X-FEM requires some enhancements of the regular FEM implementation:

* A specific integration scheme on elements whose nodes are enriched must be used [4, 10].

* The number of degrees of freedom per element depends on the location of the element with respect to the crack. The use of an objected-oriented language such as C++ is quite convenient in dealing with the variable number of element degrees of freedom (in space and time as the crack is growing).

* Conditioning issues are raised when many layers of elements are enriched around the crack tip (a pre-conditioning technique for this event is provided in [10]).

When dynamic loading is applied, additional issues must be addressed for the X-FEM to work properly. For transient analysis, one usually distinguishes between implicit and explicit algorithms. The standard Newmark implicit approach is known to be unconditionally stable in the finite element context for stationary cracks. As cracks are growing and remeshing is used, stability may be obtained, provided that the energy release rate at the crack tip is properly computed [11]. A nice thing about the X-FEM is that the mesh is held fixed as the crack grows. The number of enriched degrees of freedom is, however, growing during computation. The stability of X-FEM with implicit algorithms was studied in [12] and with explicit algorithms in [13, 14]. The latter paper [14] introduces a special lumping technique leading to a critical time step smaller (factor $1/\sqrt{2}$) than the one in the FEM case. In this paper we shall introduce a new lumping technique allowing one to use the same critical time step as in FEM (considering only Heaviside enrichment).

The paper is organized as follows. In Section 2, the critical time step is theoretically studied with two different lumping techniques for simplex elements (1D, 2D and 3D). Section 3 deals with effects of lumping techniques on wave propagation through discontinuity, and Section 4 is dedicated to numerical experiments: a compact compression specimen with both different lumping techniques is studied using two different element types (triangular and quadrangular).

## 2. EXPLICIT DYNAMICS FOR X-FEM

To capture the discontinuity in the displacement field at the crack tip, the partition of unity property of the finite element shape functions [1] is used to enrich the spatial interpolation as follows:

$$ \mathbf{u} = \sum_{i \in \mathcal{N}} N_i(\mathbf{x}) \bar{\mathbf{U}}_i + \sum_{j \in \mathcal{N}_{\text{cut}}} N_j(\mathbf{x}) H(\mathbf{x}) \tilde{\mathbf{U}}_j \quad (1) $$

In this equation, $\mathbf{u}$ is the approximate displacement field, $\mathcal{N}$ the set of nodes used to discretize the domain, and $N_i$ the finite element shape function associated with node $i$. $\bar{\mathbf{U}}_i$ are classical degrees of freedom, whereas $\tilde{\mathbf{U}}_j$ are additional degrees of freedom supported by nodes $j$ in the set $\mathcal{N}_{\text{cut}}$ associated with the enrichment function $H$, see Figure 1. Only one enriched function is used here to describe the discontinuity in the structure. We choose the jump enrichment function as unit and
<figure><img src="image_2.png" /><figcaption>Figure 1. Selected nodes $N_{cut}$ for the Heaviside enrichment.</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>Figure 2. Classical ($f_I, f_{II}$) and enriched ($f_{I'}, f_{II'}$) approximation functions on a cracked 1D element.</figcaption></figure>

symmetric. The body $\Omega$ is separated by the discontinuity into two parts: $\Omega^+$ and $\Omega^-$. Then the chosen enriched function $H$ is defined as

$$H(\mathbf{x}) = \begin{cases} 1 & \text{if } x \in \Omega^+ \\ -1 & \text{if } x \in \Omega^- \end{cases} \qquad (2)$$

Other jump functions could be used [15] but the following study is valid only with $H$. It is well known that the critical time step of an explicit integration scheme for dynamics may be estimated by computing the minimum value of the critical time step for all elements in the mesh taken separately. Moreover, this estimate is an upper bound for the exact critical time step. This result was obtained in [16] as an extension of a theorem first given by Rayleigh.

The paragraph above remains valid if the approximation is obtained through the X-FEM. We can thus concentrate on a single cracked element.

## 2.1. One-dimensional element

### 2.1.1. The shape functions.
First of all we start with a 1D cracked element shown in Figure 2. The length of the element is $L$, its density $\rho$, its section $S$ and its Young modulus $E$. The crack cuts the element into two sub-segments of sizes $\varepsilon L$ and $(1-\varepsilon)L$, respectively. Using the Heaviside enrichment [4], four approximation functions are defined over the element (Figure 2): two classical and two enriched. We assume that the ‘Heaviside’-type enrichment function is +1 on the left side of the bar and -1 on the right side. Equivalently, the four functions shown in Figure 3 may be used. The latter base is in the spirit of Hansbo and Hansbo [17]. The equivalence of the basis depicted in Figures 2 and 3 was proven in [18].
<figure><img src="image_4.png" /><figcaption>Figure 3. Truncated function basis for a 1D cracked element (shadow node version). This basis is equivalent to the one depicted in Figure 2.</figcaption></figure>

Moreover, Song et al. [19] and Areais et al. [20] contain a proof of the equivalence between Hansbo and Hansbo [17] and X-FEM scheme and present numerous results for explicit/implicit dynamic X-FEM simulations. The two basis functions are related by

$$ \begin{pmatrix} f_I \\ f_{I'} \\ f_{II} \\ f_{II'} \end{pmatrix} = \begin{pmatrix} 1 & 1 & 0 & 0 \\ 1 & -1 & 0 & 0 \\ 0 & 0 & 1 & 1 \\ 0 & 0 & -1 & 1 \end{pmatrix} \cdot \begin{pmatrix} f_I \\ f_{I'} \\ f_{II} \\ f_{II'} \end{pmatrix} = \mathbb{P} \cdot \begin{pmatrix} f_I \\ f_{I'} \\ f_{II} \\ f_{II'} \end{pmatrix} \quad (3) $$

Using the approximation functions depicted in Figure 2, an approximated field $u$ on the 1D element reads

$$ u = u_I f_I + u_{I'} f_{I'} + u_{II} f_{II} + u_{II'} f_{II'} \quad (4) $$

Alternatively, one can use the equivalent basis depicted in Figure 3. The approximation reads with this basis as follows:

$$ u = u_1 f_1 + u_2 f_2 + u_3 f_3 + u_4 f_4 \quad (5) $$

Using link (3), the degrees of freedom of the approximations are related by

$$ \begin{pmatrix} u_1 \\ u_{1'} \\ u_2 \\ u_{2'} \end{pmatrix} = \begin{pmatrix} 1 & 1 & 0 & 0 \\ 1 & -1 & 0 & 0 \\ 0 & 0 & 1 & -1 \\ 0 & 0 & 1 & 1 \end{pmatrix} \cdot \begin{pmatrix} u_I \\ u_{I'} \\ u_{II} \\ u_{II'} \end{pmatrix} = \mathbb{P}^T \begin{pmatrix} u_I \\ u_{I'} \\ u_{II} \\ u_{II'} \end{pmatrix} \quad (6) $$

It can be noted that the basis functions used in Figure 3 are as if the two sub-segments were completely independent. The approximation on the left segment is as if the right one was void and vice-versa. On each sub-segment the approximation is exactly the one used to model voids with non-matching meshes as described in [21]. It is thus tempting to use the lumping strategy developed in [22] for voids leading to a critical time step independent of $\epsilon$ and value $h/c$, where $c$ is the wave velocity $c = \sqrt{E/\rho}$, i.e. the same critical time step as for a regular element.
2.1.2. The mass and stiffness matrices. Following [22], the lumped mass matrix in our case is expressed as

$$ M_{1,1',2,2'} = \frac{\rho SL}{2} \begin{pmatrix} \varepsilon & 0 & 0 & 0 \\ 0 & 1-\varepsilon & 0 & 0 \\ 0 & 0 & 1-\varepsilon & 0 \\ 0 & 0 & 0 & \varepsilon \end{pmatrix} \qquad (7) $$

The stiffness matrix is expressed as (see [14]):

$$ K_{1,1',2,2'} = \frac{ES}{L} \begin{pmatrix} \varepsilon & 0 & 0 & -\varepsilon \\ 0 & 1-\varepsilon & \varepsilon-1 & 0 \\ 0 & \varepsilon-1 & 1-\varepsilon & 0 \\ -\varepsilon & 0 & 0 & \varepsilon \end{pmatrix} \qquad (8) $$

The mass and stiffness matrices determine the critical time step of the element.

2.1.3. The critical time step. For an explicit scheme, the critical time step is computed as $2/\omega_{max}$, where $\omega_{max}$ is linked to the maximum eigenvalue of the generalized system:

$$ \det(K_{1,1',2,2'} - \omega^2 M_{1,1',2,2'}) = 0 \qquad (9) $$

Thus, we obtain in our case

$$ \varepsilon^2(1-\varepsilon)^2\alpha^2(\alpha-2)^2=0 \qquad (10) $$

where

$$ \alpha = \frac{\omega^2 \rho S L^2}{2 E S} \qquad (11) $$

Hence, the biggest value for $\omega$ is when $\alpha$ is 2, thus, $\omega_{max} = (2/L)\sqrt{E/\rho}$. Then the critical time step is $L\sqrt{\rho/E}$, which exactly corresponds to the critical time step of the 1D finite element problem.

The result obtained is promising: the critical time step of a 1D element whose nodes are enriched is exactly the same as that of a standard element. Moreover, this result does not depend on the shape function basis used: $(1, 1', 2, 2')$ or $(I, I', II, II')$. The proof is

$$ \begin{aligned} \det(K_{1,1',2,2'} - \omega^2 M_{1,1',2,2'}) &= \det(\mathbb{P}^T K_{I,I',II,II'} \mathbb{P}) - \omega^2 \mathbb{P}^T M_{I,I',II,II'} \mathbb{P}) \\ &= \det(\mathbb{P}^T (K_{I,I',II,II'} - \omega^2 M_{I,I',II,II'})) \mathbb{P}) \\ &= 16 \det(K_{I,I',II,II'} - \omega^2 M_{I,I',II,II'}) \end{aligned} \qquad (12) $$

where matrix $\mathbb{P}$ is defined in Equation (6) such that $\det(\mathbb{P} \cdot \mathbb{P}^T) = 16$ and $\mathbb{P} \cdot \mathbb{P}^T = \mathbb{P}^T \cdot \mathbb{P} = 2I_4$. Note that if the modified enrichment proposed in [3] (which was shown to be equivalent to the Hansbo basis in [18]) were used and the mass matrix were lumped in this enriched basis following Equation
(7), the result would be exactly what we obtain by lumping the mass matrix in the (1, 1', 2, 2') basis. Hence, going back to the X-FEM basis functions, we apply the transformation to the above mass matrix

$$ M_{I,I',II,II'} = \begin{pmatrix} 1 & 1 & 0 & 0 \\ 1 & -1 & 0 & 0 \\ 0 & 0 & 1 & -1 \\ 0 & 0 & 1 & 1 \end{pmatrix}^T \cdot M_{1,1',2,2'} \cdot \begin{pmatrix} 1 & 1 & 0 & 0 \\ 1 & -1 & 0 & 0 \\ 0 & 0 & 1 & -1 \\ 0 & 0 & 1 & 1 \end{pmatrix} \quad (13) $$

$$ = \frac{\rho SL}{2} \begin{pmatrix} 1 & 2\epsilon - 1 & 0 & 0 \\ 2\epsilon - 1 & 1 & 0 & 0 \\ 0 & 0 & 1 & 2\epsilon - 1 \\ 0 & 0 & 2\epsilon - 1 & 1 \end{pmatrix} \quad (14) $$

We did obtain a block-diagonal mass matrix, the blocks being of size 2.

Finally, had we decided that node 1 was on the negative side of the 'Heaviside', the results would have been different. The results would have been the one obtained by replacing $\epsilon$ by $1-\epsilon$ in (14). Result (14) is thus valid whatever the choice of sign for the Heaviside, provided that $\epsilon$ is always defined as the matter fraction on the positive Heaviside side.

### 2.1.4. Kinetic energy conservation.
This part aims at checking that the lumping techniques conserve the kinetic energy for rigid motions with a discontinuity.

In the classical finite element setting, the sum of the entries in a mass matrix is always equal to the total mass. The reason for this is that when all degrees of freedom are set to one, we obtain a rigid translation mode. In the X-FEM case, the sum of the entries of the matrix is not equal to the total mass because setting all degrees to 1 does not give a rigid mode. However, if we sum the entries corresponding to a rigid mode (first and third lines and column) we recover the total mass. Let us consider a 1D element with two nodes. Each node has ordinary degrees of freedom corresponding to the shape functions $f_I$ and $f_{II}$ and additional degrees of freedom corresponding to the enriched function $H$. The fraction ratio of this element is $\epsilon \in [0, 1]$. First we consider a rigid motion at velocity $V$ (as if there was no crack). The exact kinetic energy is

$$ E_c = \frac{1}{2} m V^2 \quad (15) $$

where $m = \rho SL$ is the mass of the element, $S$ its section, $L$ its length and $\rho$ its density. And the discretized velocity field corresponding to this motion is described by the vector (in the Hansbo basis (1, 1', 2, 2') as shown in Figure 3): $[U]=[V, V, V, V]^T$. Hence, the discretized kinetic energy reads

$$ \begin{aligned} E_{\text{dis}} &= \frac{1}{2}[U]^T \cdot M_{1,1',2,2'} \cdot [U] = \frac{1}{2} \frac{m}{2}(\epsilon + 1 - \epsilon + 1 - \epsilon + \epsilon)V^2 \\ &= \frac{1}{2}mV^2 = E_c \end{aligned} \quad (16) $$
Let us now consider the rigid motion of a cut element whose parts move away at velocity V. The discretized velocity vector is now $[U]=[V, -V, -V, V]^T$. And the discretized kinetic energy can be expressed as

$$ 
\begin{aligned}
E_{\text{dis}} &= \frac{1}{2}[U]^T \cdot \mathbb{M}_{1,1',2,2'}.[U] = \frac{1}{2} \frac{m}{2} (\varepsilon V^2 + (1-\varepsilon)(-V)^2 + (1-\varepsilon)(-V)^2 + \varepsilon V^2) \\
&= \frac{1}{2} m V^2 = E_c
\end{aligned} 
\quad (17) $$

The case corresponding to the use of the other lumped mass matrix has already been verified in [14].

Let us now consider a more complex structure. Figure 4 presents a bar composed of three elements; the middle one is cut. Each element has its length equal to L, its Young's modulus E, its section S and its density ρ. The length of the bar is l=3L; hence, its mass is m=ρSL. Here we check the conservation of kinetic energy for the whole bar. Figure 5 presents the X-FEM shape functions: four standard and two enriched functions. Figure 6 presents the six shape functions used for this 1D case: the Hansbo basis is chosen here. The discretized displacement vector is of size 6 and is expressed as

$$ [U]_{\text{diag}} = [u_1 \ u_2 \ a_2 \ u_3 \ a_3 \ u_4]^T \quad (18) $$

Hence, the approximate displacement field U reads

$$ U = u_1 \cdot N_1 + u_2 \cdot N_2 + a_2 \cdot H \cdot N_2 + u_3 \cdot N_3 + a_3 \cdot H \cdot N_3 + u_4 \cdot N_4 \quad (19) $$

The diagonal mass matrix [14] for the whole structure, in the standard basis (N₁, N₂, HN₂, N₃, HN₃, N₄), is as follows (with l=3L):

$$ \mathbb{M}_{\text{diag}} = \frac{\rho SL}{2} \begin{pmatrix} 1 & 0 & 0 & 0 & 0 & 0 \\ 0 & 2 & 0 & 0 & 0 & 0 \\ 0 & 0 & 2 & 0 & 0 & 0 \\ 0 & 0 & 0 & 2 & 0 & 0 \\ 0 & 0 & 0 & 0 & 2 & 0 \\ 0 & 0 & 0 & 0 & 0 & 1 \end{pmatrix} \quad (20) $$

The block diagonal mass matrix [22] is

$$ \mathbb{M}_{\text{block-diag}} = \frac{\rho SL}{2} \begin{pmatrix} 1 & 0 & 0 & 0 & 0 & 0 \\ 0 & 2 & 2\epsilon - 1 & 0 & 0 & 0 \\ 0 & 2\epsilon - 1 & 2 & 0 & 0 & 0 \\ 0 & 0 & 0 & 2 & 2\epsilon - 1 & 0 \\ 0 & 0 & 0 & 2\epsilon - 1 & 2 & 0 \\ 0 & 0 & 0 & 0 & 0 & 1 \end{pmatrix} \quad (21) $$
<figure><img src="image_5.png" /><figcaption>Figure 4. A three-element structure cut by a crack.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Figure 5. Shape functions for the three elements mesh: $N_1, N_2, HN_2, N_3, HN_3$ and $N_4$.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Figure 6. Shape functions (Hansbo basis for the enhanced element) for the three elements mesh.</figcaption></figure>

$M_{\text{diag}}$ is the diagonal mass matrix obtained using the technique developed in [14], and $M_{\text{block-diag}}$ the block-diagonal mass matrix obtained in [22]. As in the last paragraph let us consider two rigid motions.

*Motion 1:* First the whole structure moves as a rigid body in the same direction at velocity $V$. Hence, the velocity vector is $[v_1]=[V, V, 0, V, 0, V]^T$. The discretized kinetic energy is evaluated by $\frac{1}{2}[v_1]^T M_{\text{diag}} [v_1]$, and for both lumping techniques, one obtains

$$Ec_{\text{motion-1}}^{\text{diag}} = \frac{1}{2}[v_1]^T M_{\text{diag}} [v_1] = \frac{1}{2} \rho S I V^2 = Ec_{\text{motion-1}}^{\text{exact}} \quad (22)$$
$$E c_{\text{motion-1}}^{\text{block-diag}} = \frac{1}{2}[v_1]^T \cdot \mathbb{M}_{\text{block-diag}} \cdot [v_1] = \frac{1}{2} \rho S I V^2 = E c_{\text{motion-1}}^{\text{exact}} \quad (23)$$

**Motion 2:** For the second motion, the structure at the left of the crack moves at velocity $-V$, whereas the right part moves at velocity $V$. Thus, the two parts move away at velocity $V$. Hence, the velocity vector is expressed as $[v_2]=[-V, 0, -V, 0, -V, V]^T$. The discretized kinetic energy is evaluated by $\frac{1}{2}[v_2]^T \cdot \mathbb{M} \cdot [v_2]$ for both lumping techniques:

$$E c_{\text{motion-2}}^{\text{diag}} = \frac{1}{2}[v_2]^T \cdot \mathbb{M}_{\text{diag}} \cdot [v_2] = \frac{1}{2} \rho S I V^2 = E c_{\text{motion-2}}^{\text{exact}} \quad (24)$$

$$E c_{\text{motion-2}}^{\text{block-diag}} = \frac{1}{2}[v_2]^T \cdot \mathbb{M}_{\text{block-diag}} \cdot [v_2] = \frac{1}{2} \rho S I V^2 = E c_{\text{motion-2}}^{\text{exact}} \quad (25)$$

Hence, one has to note that the two lumped mass matrices allow to conserve the kinetic energy for both rigid body motions; moreover, the crack is taken into account in one motion. These verifications about kinetic energy conservation mean that both lumping techniques do not add mass.

To conclude with the 1D case, the lumping technique developed in [22] allows one to obtain the same critical time step for an element with additional degrees of freedom than for a standard element. This result improves the technique explained in [14] which was only $\sqrt{2}$ smaller than the standard one. However, the cost here to obtain a better result is that the mass matrix depends on the fraction ratio, whereas the mass matrix in [14] was simply diagonal constant. In addition, we show that these two lumping techniques allow one to preserve the kinetic energy for rigid motions.

Now the use of each technique can be motivated by different reasons: the use of block-diagonal matrix will be preferable, except for pure explicit codes which do not own a matrix (as LS DYNA, RADIOSS, EUROPLEXUS). For these codes, a constant diagonal mass matrix will be enough unless they were able to store 2×2 matrices for each enriched node.

## 2.2. Two-dimensional elements

This part studies the critical time step for triangular and quadrangular elements. Moving to 2D elements, we shall see that the mass matrix keeps the same topology as in (14).

### 2.2.1. Triangular element.
Consider a triangular element depicted in Figure 7: a triangular element cut by a crack and characterized by its fraction ratio $\epsilon$. As in the 1D case, we use roman indices to denote the classical shape functions and the roman prime indices to denote the enriched shape function. The arabic indices denote the truncated shape functions. For instance, the function $f_1$ is equal to the classical shape function $f_I$ on the sub-triangle containing node 1 and is zero on the other sub-triangle. Finally, the prime arabic indices indicate the ‘complementary’ truncated functions. For instance, $f_{1'}$ is such that $f_1 + f_{1'} = f_I$.

The truncated and enriched basis functions are related by

$$f_I = f_1 + f_{1'} \quad (26)$$

$$f_{I'} = H(x_I)(f_1 - f_{1'}) \quad (27)$$

$$f_{II} = f_2 + f_{2'} \quad (28)$$
<figure><img src="image_8.png" /><figcaption>Figure 7. A triangular element of area S cut by a crack. The fraction ratio is denoted by $\varepsilon$.</figcaption></figure>

$$f_{II'} = H(x_{II})(f_2 - f_{2'}) \qquad (29)$$

$$f_{III} = f_3 + f_{3'} \qquad (30)$$

$$f_{III'} = H(x_{III})(f_3 - f_{3'}) \qquad (31)$$

where $H(\cdot)$ indicates the sign of the Heaviside at the corresponding node. The mass matrix in the truncated base reads

$$M_{1,1',2,2',3,3'} = \frac{\rho S}{3} \begin{pmatrix} \varepsilon & 0 & 0 & 0 & 0 & 0 \\ 0 & 1-\varepsilon & 0 & 0 & 0 & 0 \\ 0 & 0 & 1-\varepsilon & 0 & 0 & 0 \\ 0 & 0 & 0 & \varepsilon & 0 & 0 \\ 0 & 0 & 0 & 0 & 1-\varepsilon & 0 \\ 0 & 0 & 0 & 0 & 0 & \varepsilon \end{pmatrix} \qquad (32)$$

Applying the proper transformations, based on Equations (26)–(31), we obtain the following mass matrix for the enriched basis:

$$M_{I,I',II,II',III,III'} = \frac{\rho S}{3} \begin{pmatrix} 1 & 2\varepsilon - 1 & 0 & 0 & 0 & 0 \\ 2\varepsilon - 1 & 1 & 0 & 0 & 0 & 0 \\ 0 & 0 & 1 & 2\varepsilon - 1 & 0 & 0 \\ 0 & 0 & 2\varepsilon - 1 & 1 & 0 & 0 \\ 0 & 0 & 0 & 0 & 1 & 2\varepsilon - 1 \\ 0 & 0 & 0 & 0 & 2\varepsilon - 1 & 1 \end{pmatrix} \qquad (33)$$

We have assumed that the material fraction $\varepsilon$ is computed on the positive side of the Heaviside. Otherwise, $\varepsilon$ must be replaced by $1-\varepsilon$. The obtained mass matrix is basically the same as in the
<figure><img src="image_9.png" /><figcaption>Figure 8. Normalized critical time step for a triangular element as a function of the position of the crack as explained in Figure 7.</figcaption></figure>

1D case. For 3D elements, the 2 × 2 block-diagonal matrix will not change, except for the factor in front of the matrix, which will be the element volume divided by the number of nodes of the element. The results for a triangular element are presented in Figure 8. Figure 8 shows the critical time step normalized by the finite element critical time step obtained with a lumped mass matrix. The two results are for the two different lumping techniques, which use a diagonal mass matrix [14] and a block diagonal one [22], respectively. It underlines the fact that with the block-diagonal mass matrix the critical time step of the enriched triangular element is the same as that of the triangular finite element, whatever the position of the crack in the element. Hence, the result is similar to the 1D case. Moreover, Figure 8 shows that the critical time step for the diagonal case is quite smaller than for the block-diagonal case. Note also that the block-diagonal 2 × 2 matrix is not the same for all the enriched nodes. It depends on the crack path.

To conclude, using a block-diagonal mass matrix, one obtains the same critical time step as in the finite element case, hence, around twice the value of the critical time step obtained using constant diagonal mass matrix. The block-diagonal mass matrix is composed of 2 × 2 matrices, each of which depends on the fraction ratios $\varepsilon$ of the elements around the considered nodes. The results shown in this part are quite interesting because the critical time step for an element with additional degrees of freedom is the same as for a standard one. It means that an enriched triangular elements and a standard triangular element have the same critical time step when the mass matrix is lumped by the appropriate technique.

2.2.2. *Quadrangular element.* Let us now consider a quadrangular element. Figure 9 describes the fraction ratio $\varepsilon$ on a quadrangular element cut by a crack. The critical time steps are shown in Figure 10 obtained by both diagonal and block-diagonal mass matrices. Both results are similar,
<figure><img src="image_10.png" /><figcaption>Figure 9. A quadrangular element of area S cut by a crack.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Figure 10. Normalized critical time step for a quadrangular element as a function of the position of the crack as explained in Figure 11</figcaption></figure>

but it has to be noted that the use of the block-diagonal mass matrix allows one to obtain a bigger critical time step. We must note that the critical time step of the standard element is never reached for any fraction ratio $\varepsilon$, whereas it was for the triangular case. Thus the biggest critical time step is reached for $\varepsilon=0.5$ and its normalized value is 0.93.

2.2.3. Comparison between triangular and quadrangular elements. This part compares the uses of one quadrangular and two triangular elements whose nodes are enriched. Figure 11 describes
the same 4-nodes structure with the two different element types. Figure 12 shows the critical time step for one quadrangular and two triangular elements whose nodes are enriched and cut by the same crack. We can note that the critical time step obtained using a block-diagonal mass matrix is still upper than the one by using a diagonal mass matrix. Moreover, the quadrangular element allows one to obtain a bigger critical time step than using two triangular elements. Nevertheless, the worst critical time step for a quadrangular element is obtained when the crack is close to a node, and its value is around the same as the one obtained using two triangular elements.

<figure><img src="image_12.png" /><figcaption>Figure 11. Two triangular elements, and one quadrangular element cut by a crack.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Figure 12. Normalized critical time step (by the quadrangular finite element one) for two triangular elements as a function of the position of the crack as explained in Figure 11. The two curves of Figure 10 are also reproduced.</figcaption></figure>
### 2.3. Three-dimensional elements

Let us now move to the 3D case and the studies of critical time step of tetrahedral and cubic elements.

#### 2.3.1. Tetrahedral element.
Figure 13 presents a tetrahedral element cut by a crack and Figure 14 the critical time step normalized by the finite element critical time step obtained with a lumped mass matrix. The two results are for the two different lumping techniques, which use a diagonal

<figure><img src="image_14.png" /><figcaption>Figure 13. Three-dimensional tetrahedral element of volume $V$ cut by a crack; $\varepsilon$ is the fraction ratio on the positive side of the Heaviside function.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>Figure 14. Normalized critical time step for a tetrahedral element as a function of the position of the crack as explained in Figure 13.</figcaption></figure>
<figure><img src="image_16.png" /><figcaption>Figure 15. Three-dimensional cubic element cut by a crack.</figcaption></figure>

mass matrix [14] and a block-diagonal mass matrix [22], respectively. For this 3D element, the critical time step obtained with a block-diagonal mass matrix is exactly the same as that of the non-enriched finite element. Hence, the results in this paper are similar for all dimensions: 1D linear element, triangular element and tetrahedral element allow one to obtain an X-FEM critical time step equal to the FEM critical time step.

### 2.3.2. Cubic element.
Figure 15 presents a cubic element cut by a crack and Figure 16 shows the critical time step of this cut element as a function of the position of the crack in the element. The time step is normalized by the critical time step of the classical finite element. The results of the critical time step of the enhanced element are available for the two different lumping techniques. These techniques give different results as shown in Figure 16. The critical time step results are quite similar to those of the quadrangular element. None of the two techniques are able to obtain the same critical time step as the one of the finite element problem. Nevertheless, the use of a block-diagonal mass matrix allows a bigger time step.

## 2.4. Results analysis

In this section, we take a deeper look at the reasons as to why we can obtain the same critical time step as for the finite element case for some elements: triangular and tetrahedral elements with the block-diagonal mass matrix give the same critical time step as for the finite element. The fact that the critical time step of a standard element can be obtained with an enriched one is directly related to the shape functions of the element. Using the Hansbo basis for an element whose nodes are enriched, we have for a 1D element (as shown in Figure 3)

$$ \frac{1}{L} \int_{0}^{L} \frac{\partial f_{1}(x)}{\partial x} dx = \int_{0}^{\epsilon} dX = \epsilon $$

$$ \frac{1}{L} \int_{0}^{L} \frac{\partial f_{1'}(x)}{\partial x} dx = \int_{0}^{\epsilon} dX = \epsilon $$
<figure><img src="image_17.png" /><figcaption>Figure 16. Normalized critical time step for a cubic element as a function of the fraction ratio $\epsilon$.</figcaption></figure>

$$ \frac{1}{L} \int_{0}^{L} \frac{\partial f_{2}(x)}{\partial x} dx = \int_{\epsilon}^{1} dX = 1 - \epsilon \qquad (34) $$

$$ \frac{1}{L} \int_{0}^{L} \frac{\partial f_{2'}(x)}{\partial x} dx = \int_{\epsilon}^{1} dX = 1 - \epsilon $$

Thus, the X-FEM stiffness matrix is proportional to the FEM stiffness matrix on each side of the discontinuity

$$ K_{1,2',1',2} = \frac{ES}{L} \begin{pmatrix} \epsilon & -\epsilon & 0 & 0 \\ -\epsilon & \epsilon & 0 & 0 \\ 0 & 0 & 1-\epsilon & \epsilon-1 \\ 0 & 0 & \epsilon-1 & 1-\epsilon \end{pmatrix} = \begin{pmatrix} \epsilon K_{\text{FEM}} & \circ \\ \circ & (1-\epsilon) K_{\text{FEM}} \end{pmatrix} \qquad (35) $$

where $K_{\text{FEM}}$ is given in the Appendix. The same may be said for the lumped mass matrix given in Equation (7). Thus, both stiffness and mass matrices are proportional to the FEM on each side of the discontinuity. Then, the resolution of Equation (9) is similar to the resolution in the classical finite element case, which explains why the X-FEM critical time step is the same as the FEM one.

For simplex elements (1D bar, triangular and tetrahedral), the standard shape functions are linear in space and are written as the following functions $\Phi_k$ (where $k \in \{1...n_{\text{nodes}}\}$, $n_{\text{nodes}}$ is the number of nodes in the element):

$$ \Phi_k(x, y, z) = I_0 + I_1 \cdot x + I_2 \cdot y + I_3 \cdot z \qquad (36) $$
<table><caption>Table I. Critical time steps for different elements: standard and enriched (Young's modulus $E=1$, length $L=1$, density $\rho=1$, the Poisson ratio $\nu=0,3$).</caption><thead><tr><th rowspan="2">Element</th><th colspan="2">Δt<sub>c</sub><sup>FEM</sup></th><th colspan="2">Δt<sub>c</sub><sup>XFEM</sup> normalized by Δt<sub>c</sub><sup>FEM lump</sup></th></tr><tr><th>Δt<sub>c</sub><sup>FEM</sup><br>M standard</th><th>Δt<sub>c</sub><sup>FEM lump</sup><br>M lumped</th><th>M block diagonal</th><th>M diagonal</th></tr></thead><tbody><tr><td>One-dimensional</td><td>0.5774</td><td>1</td><td>1</td><td>&gt; 0.707</td></tr><tr><td>Triangular</td><td>0.2819</td><td>0.5638</td><td>1</td><td>&gt; 0.707</td></tr><tr><td>Quadrangular</td><td>0.4354</td><td>0.7601</td><td>&gt; 0.77</td><td>&gt; 0.707</td></tr><tr><td>Tetrahedral</td><td>0.1919</td><td>0.4071</td><td>1</td><td>&gt; 0.707</td></tr><tr><td>Cubic</td><td>0.3651</td><td>0.6324</td><td>&gt; 0.76</td><td>&gt; 0.707</td></tr></tbody></table>

Hence, the gradients of the shape functions are constant as explained in Equation (34). This is the reason as to why the X-FEM stiffness matrix is proportional to the finite element one.

Note that this is not true for the quadrangular and cubic elements because the gradients of shape functions (used to compute the stiffness matrix) are not constant in space; thus, the X-FEM stiffness matrices are not proportional to the FEM matrices. The shape functions for these (quadrangular then cubic) elements are not linear and may represent the following polynomials:

$$ \Phi_k^{\text{quadrangular}}(x, y) = I_0 + I_1 \cdot x + I_2 \cdot y + I_3 \cdot x \cdot y \quad (37) $$

$$ \Phi_k^{\text{cubic}}(x, y, z) = I_0 + I_1 \cdot x + I_2 \cdot y + I_3 \cdot z + I_4 \cdot x \cdot y + I_5 \cdot x \cdot z + I_6 \cdot y \cdot z + I_7 \cdot x \cdot y \cdot z \quad (38) $$

These results are summarized in Table I, which presents the normalized critical time step for the five studied elements: 1D linear, triangular, tetrahedral, quadrangular and cubic elements, with the two lumping techniques (diagonal and block diagonal).

### 3. WAVE TRANSFER THROUGH THE DISCONTINUITY

First the mass matrix is used in the explicit algorithm as

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

$$ M \cdot \ddot{U}_{n+1} = F_{n+1}^{\text{ext}} - F_{n+1}^{\text{int}} \quad (40) $$

$$ \ddot{U}_{n+1} = \ddot{U}_n + \frac{\Delta t}{2}(\ddot{U}_{n+1} + \ddot{U}_n) \quad (41) $$

We now check that the dynamic wave propagation is not transmitted through the discontinuity. This has to be checked for the two different lumping techniques.

#### 3.1. One-dimensional elementary case

Figure 17 presents the first case for wave propagation: a 6-element structure in 1D (the third element is enriched and cut by a crack). The characteristics are Young's modulus $E=3.24$ GPa,
<figure><img src="image_18.png" /><figcaption>Figure 17. A 6-element structure cut by a crack.</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>Figure 18. Shape functions for the 6-element mesh: N₁, N₂, N₃, N₄, HN₃, HN₄, N₅, N₆ and N₇.</figcaption></figure>

density ρ=1190kg/m³, length L=l/6=10mm and section S=2mm². m=ρS6L=ρSl is the mass of the whole structure. The critical time step is given by Δt_c = L/√ρ/E; here Δt_c = 6.06μs. The time step for calculation is taken as 1μs, and the total time is 50μs. The fraction ratio ε is 0.2 in the third element. The loading is an initial velocity V₀=10m/s at node 7. The structure is assumed to be linear elastic. Figure 18 presents the standard shape functions of the 6-element mesh: there are seven standard shape functions (one per node) and two enriched ones (for the cut element nodes). The right extremity of the structure is loaded with an initial velocity. Afterwards the tension wave propagates and reaches the discontinuity. We would like to check whether or not the structure on the left part of the discontinuity is affected by the wave near the discontinuity. In the explicit algorithm, the important equation to know whether the wave will propagate through the discontinuity is as follows:

$$ \ddot{U}_{n+1} = M^{-1}(F_{n+1}^{\text{ext}} - F_{n+1}^{\text{int}}) \quad (42) $$

And in our case, there is no external force, and the material model is elastic linear; hence, this term is just written as

$$ \ddot{U}_{n+1} = -M^{-1}F_{n+1}^{\text{int}} \quad (43) $$
<figure><img src="image_20.png" /><figcaption>Figure 19. Wave propagation on the 6-element mesh at time 25 µs: displacement, velocity, acceleration, stress and strain fields.</figcaption></figure>

The fact that wave goes through the discontinuity would be due to a non-zero right-hand side for node 3 in Equation (43). The results of wave propagation are summed up in Figure 19, which presents five graphs for the two different mass matrices: diagonal and block diagonal. The five graphs in Figure 19 show the displacement, velocity, acceleration, stress and strain fields along the structure at time 25 µs. Note that the two different lumping techniques give similar results.
<table><caption>Table II. The induction proposition.</caption><thead><tr><th></th><th>At the left side of the discontinuity<br>(H=+1)</th><th>At the right side of the discontinuity<br>(H=-1)</th></tr></thead><tbody><tr><td>Displacement</td><td>u<sub>1</sub>=u<sub>2</sub>=u<sub>3</sub>+a<sub>3</sub>=0</td><td>u<sub>4</sub>=-a<sub>4</sub></td></tr><tr><td>Velocity</td><td>u&#x0307;<sub>1</sub>=u&#x0307;<sub>2</sub>=u&#x0307;<sub>3</sub>+a&#x0307;<sub>3</sub>=0</td><td>u&#x0307;<sub>4</sub>=-a&#x0307;<sub>4</sub></td></tr><tr><td>Acceleration</td><td>u&#x0307;<sub>1</sub>=u&#x0307;<sub>2</sub>=u&#x0307;<sub>3</sub>+a&#x0307;<sub>3</sub>=0</td><td>u&#x0307;<sub>4</sub>=-a&#x0307;<sub>4</sub></td></tr></tbody></table>

This is due to the mass matrix, which is not the same in the two cases. Figure 19 shows that there is no wave transfer through the discontinuity. Indeed nothing happens on the left side of the discontinuity. For this simple 1D example, let us now prove analytically that no wave transfer occurs with the two lumping techniques.

*Proof by mathematical induction*

On the basis of Figure 17, the displacement discretization is considered as follows:

$$ 
\begin{aligned}
U(x) &= [u_1, u_2, u_3, a_3, u_4, a_4, u_5, u_6, u_7] \cdot [N_1, N_2, N_3, HN_3, N_4, HN_4, N_5, N_6, N_7]^T \\
&= [U] \cdot [N]^T
\end{aligned} 
\quad (44) 
$$

and the *proposition to prove* is chosen as shown in Table II. One notes that the displacement of node 3 (respectively, node 4) is $u_3+a_3$ (respectively, $u_4-a_4$) because of the choice of the enrichment function $H$, which is +1 on the left side of the discontinuity and -1 on the right side. To initialize the proof, we consider that the statement described in Table II is verified at time 0. Now we suppose that it is also verified at time $n$, and we prove that it is also true at time $n+1$. First, to pass from time $n$ to $n+1$, Equation (39) gives the displacement at time $n+1$. We easily obtain that the displacement field $[U_{n+1}]$ preserves the properties of Table II at time $n+1$. Second, Equation (43) is explained for both lumping techniques:

$$ M \cdot [\ddot{U}_{n+1}] = -K \cdot [U_{n+1}] \quad (45) $$

where $[U_{n+1}]=[u_1, u_2, u_3, a_3, u_4, a_4, u_5, u_6, u_7]$; the stiffness matrix is expressed as follows:

$$ K = \frac{6ES}{l} \begin{bmatrix} 1 & -1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ -1 & 2 & -1 & 0 & -1 & 0 & 0 & 0 & 0 \\ 0 & -1 & 2 & -1 & 2s & 1-2s & 0 & 0 & 0 \\ 0 & 0 & -1 & 2 & 1-2s & 2s-2 & -1 & 0 & 0 \\ 0 & -1 & 2s & 1-2s & 2 & -1 & 0 & 0 & 0 \\ 0 & 0 & 1-2s & 2s-2 & -1 & 2 & 1 & 0 & 0 \\ 0 & 0 & 0 & -1 & 0 & 1 & 2 & -1 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & -1 & 2 & -1 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & -1 & 1 \end{bmatrix} \quad (46) $$
and the lumped mass matrix for both lumping techniques is expressed as

$$ M = \frac{m}{12} \begin{bmatrix} 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 2 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 2 & \delta(2s-1) & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & \delta(2s-1) & 2 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 2 & \delta(2s-1) & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & \delta(2s-1) & 2 & 1 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 2 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 2 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 \end{bmatrix} \quad (47) $$

with $\delta$ being a parameter selecting the lumping technique: $\delta=1$ corresponds to the block-diagonal mass matrix and $\delta=0$ to the diagonal mass matrix. Afterwards, we compute Equation (45) and subsequently obtain the following relationships for the acceleration field:

$$ \ddot{u}_1 = -\frac{ES}{ml}(u_1 - u_2) \quad (48) $$

$$ \ddot{u}_2 = -\frac{2ES}{ml}(-u_1 + 2u_2 - u_3 - a_3) \quad (49) $$

$$ \ddot{u}_3 + \delta(2s-1)\ddot{a}_3 = -\frac{2ES}{ml} \left( 2s(-u_2 + 2u_3 + 2sa_3 - u_4 + (1-2s)a_4) \right) \quad (50) $$

$$ \delta(2s-1)\ddot{u}_3 + \ddot{a}_3 = -\frac{2ES}{ml}(-u_2 + 2su_3 + 2a_3 + (1-2s)u_4 - a_4) \quad (51) $$

$$ \ddot{u}_4 + \delta(2s-1)\ddot{a}_4 = -\frac{2ES}{ml}(-u_3 + (1-2s)a_3 + 2u_4 + (2s-2)a_4 - u_5) \quad (52) $$

$$ \delta(2s-1)\ddot{u}_4 + \ddot{a}_4 = -\frac{2ES}{ml}((1-2s)u_3 - a_3 + (2s-2)u_4 + 2a_4 + u_5) \quad (53) $$

From the first two equations, $\ddot{u}_1$ and $\ddot{u}_2$ do satisfy their induction proposition. Taking the sum of Equations (50) and (51), (52) and (53), respectively, one obtains

$$ [1 + \delta(2\epsilon - 1)](\ddot{u}_3 - \ddot{a}_3) = -\frac{2ES}{ml}[-2u_2 + 2(s+1)(u_3 + a_3) - 2s(a_4 + u_4)] = 0 \quad (54) $$

$$ [1 + \delta(2\epsilon - 1)](\ddot{u}_4 + \ddot{a}_4) = -\frac{2ES}{ml}2s(-u_3 - a_3 + a_4 + u_4) = 0 \quad (55) $$
<figure><img src="image_21.png" /><figcaption>Figure 20. Geometry and loading conditions.</figcaption></figure>

Thus, the acceleration properties of Table II remain valid at time $n+1$. Third, Equation 41 allows one to verify that the velocity field at time $n+1$ preserves the properties of Table II. Finally, the initial properties described in Table II are conserved along time. The fact that zero values on the left side of the discontinuity will not change for both lumping techniques means that tensile waves do not propagate through the discontinuity.

## 3.2. Two-dimensional case

Let us now consider the plane strain wave propagation block described in Figure 20. The dimensions are $L=5$ mm, $W=10$ mm and materials properties are those of PMMA: Young's modulus $E=3.24$ GPa, the Poisson ratio $\nu=0.35$ and its density $\rho=1190$ kg/m³. The loads are just impact velocity in the tension direction (see Reference [15]). This impact velocity is defined as

$$V = \begin{cases} V_0 t / t_r & \text{for } t \le t_r \\ V_0 & \text{otherwise} \end{cases} \qquad (56)$$

where $t_r=0.1$ µs and the velocity $V_0=10$ m/s. The discretization of the structure is 78 × 39, with quadrangular elements. The time step was 0.1 µs, whereas Remmers et al. [15] used a time step 100 times smaller. Let us now simulate this test with both lumping techniques: diagonal and block-diagonal matrices. Figure 21 describes the profile of the stress $\sigma_{yy}$ at the centre of the plane strain wave specimen at $t=1.5$ and $2$ µs with both lumping techniques. The reason as to why one sees only one curve in Figure 21 is the fact that results are identical for both techniques. One notes that the tensile stress average is about 25 MPa and propagates from the top side of the specimen to reach the crack at around $t=2$ µs, which agrees with [15]. Figure 21 shows that the stress wave at $t=3$ µs is properly reflected, and the lower part of the specimen remains stress free for both lumping techniques. Hence, the no-wave propagation through the discontinuity using both lumping techniques is acceptable.
<figure><img src="image_22.png" /></figure>

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

<figure><img src="image_24.png" /><figcaption>Figure 21. Stress curves σyy versus vertical position y measured along the centre of the plane strain wave specimen at t = 1.5, 2 and 3 μs, for both lumping techniques.</figcaption></figure>

4. NUMERICAL APPLICATIONS

The compact compression specimen (CCS) problem is described schematically in Figure 22. The material properties are those of PMMA: E = 5.76 GPa; ν = 0.42; and ρ = 1180 kg/m³. The force
<figure><img src="image_25.png" /><figcaption>Figure 22. Modelling of the CCS: boundary conditions and geometry (specimen thickness: 16, 5 mm).</figcaption></figure>

<figure><img src="image_26.png" /><figcaption>Figure 23. CCS meshes: triangular and quadrangular elements.</figcaption></figure>

$P_1(t)$ is due to an impact at velocity $V_0 = 20$ m/s. The CCS is assumed to be linear elastic. Although the CCS geometry is symmetrical, the deformation and, therefore, the crack path are not. This is due to the non-symmetric loading and boundary conditions. We carry the computations with the two methods: first using a block-diagonal mass matrix and a time step close to the corresponding critical time step of this method, and second the technique using a diagonal mass matrix with its corresponding critical time step. Both the results agree with the experiments [23, 24]: crack path and velocity of the crack tip.

Moreover, the CCS experiment was simulated with two different meshes as shown in Figure 23 (using triangular or quadrangular elements) and with both lumping techniques for the mass matrix: these are the four simulations. The experiment time was 120 µs and the critical time steps are 0.63 and 0.85 µs for triangular and quadrangular classical elements. One notes that the different mesh strategies (quadrangular or triangular) give exactly the same crack paths for each lumping
<figure><img src="image_27.png" /><figcaption>Figure 24. Crack paths for both lumping techniques.</figcaption></figure>

<figure><img src="image_28.png" /><figcaption>Figure 25. Crack length for both lumping techniques.</figcaption></figure>
method; hence, Figure 24 represents only two paths corresponding to the two lumping techniques. Furthermore, the results obtained with the diagonal mass matrix was compared in [14] with those obtained using an implicit time integration that stands for a reference. Hence, Figures 24 and 25 present crack paths and lengths for both techniques compared with the implicit computation results. The critical time steps corresponding to each lumping technique are different as given in Table I: the one corresponding to the block-diagonal lumping technique is bigger. However, the block-diagonal approach will require more effort at each time step since a 2 × 2 matrix based on the fraction ratios of the cut elements needs to be built (and inverted). Note that if one uses the truncated basis directly, a truly diagonal mass matrix is obtained (but it still depends on the fraction ratio). This additional effort at each time step may be less of an issue for a non-linear problem for which the internal forces, computation is where most of the time is spent. However, both ways exposed to treat X-FEM in an explicit dynamic code can be interesting. The CPU time of both may be around the same if the integration of block-diagonal mass matrix is associated with effective strategies for the inversion of block 2 × 2 and evaluation of fraction ratio. Indeed fraction ratio could be computed only once for the cut elements, and the inversion of block could also be performed once, because the quantities of a cut element (fraction ratio and mass matrix) will not change after the propagation of the crack. In fact, the computation of fraction ratio will be performed only for the new cut elements, and the inversion of the block can be performed analytically to be integrated easily in an explicit code. Thus, X-FEM with the diagonal mass matrix technique has been integrated in EUROPLEXUS.

## 5. CONCLUSION

In this paper, we introduced a new lumping technique for the mass matrices of meshes enriched with Heaviside functions within the X-FEM framework. This lumping technique yields the same critical time step at the element level as the one for *standard* elements. The lumped mass matrix is not strictly diagonal but rather block diagonal. A 2 × 2 matrix needs to be stored at each enriched node. This additional storage provides a better critical time step than the one obtained using a true diagonal mass matrix for the following elements: 1D bar, triangular, quadrangular, tetrahedral and cubic. Numerical experiments demonstrate the robustness and stability of the approach.

Then the two lumping techniques (the one developed in this paper and the other in [14]) were compared with different meshes in the 2D case. It aims at showing that the lumped mass matrix does verify some simple validations: first, studies on critical time step are checked by the simulations, and second, numerical simulations prove that no wave transfer through the displacement discontinuity appears with both lumping techniques. Furthermore, they both preserve kinetic energy for different rigid motions.

In order to choose between the two lumping techniques, the following points must be taken into account: the use of block-diagonal mass matrix is slightly more complex to implement than the true diagonal mass matrix; the diagonal mass matrix does not depend on the element fraction ratios, whereas the block-diagonal matrix does. On the other hand, the bigger time step offered by the block-diagonal approach may be interesting, in particular, for non-linear applications for which the constitutive law integration is most likely predominant compared with mass matrix inversion.
# APPENDIX A: ONE-DIMENSIONAL FINITE ELEMENT PROBLEM

For the 1D finite element problem, the consistent mass matrix and the stiffness matrix are

$$ M_{FEM} = \rho S L \begin{bmatrix} \frac{1}{3} & \frac{1}{6} \\ \frac{1}{6} & \frac{1}{3} \end{bmatrix}, \quad K_{FEM} = \frac{ES}{L} \begin{bmatrix} 1 & -1 \\ -1 & 1 \end{bmatrix} \quad (A1) $$

The lumped mass matrix is

$$ M_{FEM}^{\text{lumped}} = \rho S L \begin{bmatrix} \frac{1}{2} & 0 \\ 0 & \frac{1}{2} \end{bmatrix} \quad (A2) $$

and the corresponding stable time step is

$$ \Delta t_c^{\text{lumped}} = L \sqrt{\frac{\rho}{E}} = \sqrt{3} \Delta t_c^{\text{consistent}} = \Delta t_c^0 \quad (A3) $$

This time step is used in the normalization.

## REFERENCES

1. Babuška I, Melenk I. Partition of unity method. *International Journal for Numerical Methods in Engineering* 1997; **40**(4):727–758.

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

3. Zi G, Belytschko T. New crack-tip elements for XFEM and applications to cohesive cracks. *International Journal for Numerical Methods in Engineering* 2003; **57**:2221–2240.

4. 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**:131–150.

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

6. Stolarska M, Chopp DL, Moës N, Belytschko T. Modelling crack growth by level sets in the extended finite element method. *International Journal for Numerical Methods in Engineering* 2000; **51**:943–960.

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**:2569–2586.

8. Sukumar N, Chopp DL, Moës N, Belytschko T. Modelling holes and inclusions by level sets in the extended finite-element method. *Computer Methods in Applied Mechanics and Engineering* 2001; **190**:6183–6200.

9. Zi G, Chen H, Xu J, Belytschko T. The extended finite element method for dynamic fractures. *Shock and Vibration* 2005; **1**:9–23.

10. Béchet É, Minnebo H, Moës N, Burgardt B. Improved implementation and robustness study of the x-fem method for stress analysis around cracks. *International Journal for Numerical Methods in Engineering* 2005; **64**:1033–1056.

11. Réthoré 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.

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

13. Menouillard T, Moës N, Combescure A. An optimal explicit time stepping scheme for cracks modeled with X-FEM. *IUTAM Symposium on Discretisation Methods for Evolving Discontinuities*. Kluwer Academic Publishers: Dordrecht, 2006.

14. Menouillard T, Réthoré J, Combescure A, Bung H. Efficient explicit time stepping for the extended finite element method (X-FEM). *International Journal for Numerical Methods in Engineering* 2006; **68**:911–939.
15. Remmers JJC, de Borst R, Needleman A. An evaluation of the accuracy of discontinuous finite elements in explicit calculations. *IUTAM Symposium on Discretisation Methods for Evolving Discontinuities*. Kluwer Academic Publishers: Dordrecht, 2006.

16. Belytschko T, Smolinski P, Liu W. Stability of multi-time step partitioned integrators for first order finite element systems. *Computer Methods in Applied Mechanics and Engineering* 1985; **49**:281–297.

17. Hansbo A, Hansbo P. A finite element method for the simulation of strong and weak discontinuities in solid mechanics. *Computer Methods in Applied Mechanics and Engineering* 2004; **193**:3523–3540.

18. Areias P, Belytschko T. A comment on the article a finite element method for the simulation of strong and weak discontinuities in solid mechanics [Computer Methods in Applied Mechanics and Engineering 2004; **193** 3523–3540]. *Computer Methods in Applied Mechanics and Engineering* 2006; **195**:1275–1276.

19. Song J-H, Areias PMA, Belytschko T. A method for dynamic crack and shear band propagation with phantom nodes. *International Journal for Numerical Methods in Engineering* 2006; **67**:868–873.

20. Areias PMA, Song J-H, Belytschko T. Analysis of fracture in thin shells by overlapping paired elements. *Computer Methods in Applied Mechanics and Engineering* 2006; **195**:5343–5360.

21. Daux C, Moës N, Dolbow J, Sukumar N, Belytschko T. Arbitrary branched and intersecting cracks with the extended finite element method. *International Journal for Numerical Methods in Engineering* 2000; **48**:1741–1760.

22. Rozycki P, Moës N, Béchet É, Dubois C. X-FEM explicit dynamics for constant strain elements to alleviate mesh constraints on internal or external boundaries. *Computer Methods in Applied Mechanics and Engineering* 2007; DOI:10.1016/j.cma.2007.05.011.

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.