# X-FEM explicit dynamics for constant strain elements to alleviate mesh constraints on internal or external boundaries

Patrick Rozycki, Nicolas Moes, Eric Bechet, Céline Dubois

► To cite this version:

Patrick Rozycki, Nicolas Moes, Eric Bechet, Céline Dubois. X-FEM explicit dynamics for constant strain elements to alleviate mesh constraints on internal or external boundaries. Computer Methods in Applied Mechanics and Engineering, 2008, 197 (5), pp.349-363. 10.1016/j.cma.2007.05.011. hal-01007433

HAL Id: hal-01007433

https://hal.science/hal-01007433

Submitted on 12 Mar 2017

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.

Public Domain
# X-FEM explicit dynamics for constant strain elements to alleviate
mesh constraints on internal or external boundaries

P. Rozycki<sup>a</sup>, N. Moes<sup>a</sup>, E. Bechet<sup>c</sup>, C. Dubois<sup>a,b</sup>

<sup>a</sup> Institut de Recherche en Génie Civil et Mécanique, UMR CNRS 6183, Ecole Centrale de Nantes, 1 rue de la Noë, BP 92101, F-44321 Nantes Cedex 3, France

<sup>b</sup> PSA – Peugeot Citroën, DRIA/SARA/STEO/MAST, Route de Gisy, Bât. 62 – 1, 78943 Velizy, Villacoublay, France

<sup>c</sup> Laboratoire de Physique et Mécanique des Matériaux, UMR CNRS 7554, Université de Metz, Ile du Saulcy, F-57045 Metz Cedex 1, France

This paper deals with the use of the extended Finite Element Method (X-FEM) for rapid dynamic problems. To solve the equations of motion, a common technique is the explicit direct integration with a Newmark scheme. Since this temporal scheme is only conditionally stable, the critical time step must be determined. It is generally induced by mesh constraints. The idea of the paper is to weaken con- straints on mesh generation algorithms so that the critical time step is as large as possible. Using the X-FEM one allows a non-confor- mity between mesh and discontinuities such as cracks, holes or interfaces. In a first part, we present a summary about direct integration schemes and about the eXtended Finite Element Method. Then, we focus on the theoretical description of a 1D X-FEM finite element and its generalization to 2D and 3D finite elements. Then, dynamic numerical simulations are shown. They concern structures under impact with holes or external boundaries not exactly matched by the mesh. Comparisons are made with numerical results coming from the ABAQUS software. It shows that developments are satisfactory. We conclude with some outlooks concerning this work.

Keywords: X-FEM; Explicit schemes; Dynamics; Holes; Non-conforming meshes

## 1. Introduction

The complexity of current numerical simulations in rapid dynamics leads to problems with very high CPU time requirements and often cumbersome pre- or post-process- ing steps. Direct integration methods that are used to solve the equations of motion require relatively small time steps in order to keep the stability of the temporal scheme. The- oretically, the critical time step is based on the highest pul- sation of the structure which, in practice, depends on the smallest finite element within the discretized model (and its constitutive law). Often, mesh constraints related to the rigorous respect of complex surface geometries produce “small” finite elements. Even if the mesh can be optimized

later on in order to make it more uniform, the procedure can be difficult and may not succeed despite of the progress of meshing tools.

The objective of this paper is therefore to avoid such meshing constraints and to have an easier control on the size of the elements. For this reason, the work is based on developments of the eXtended Finite Element Method (X-FEM), which has been successfully applied to static problems exhibiting discontinuities or heterogeneities such as cracks, holes or material interfaces. The governing idea of this method is to enrich the classical FEM approxima- tion thanks to the Partition of Unity technique with specific functions representing surfaces of discontinuities or hetero- geneities. Level sets are used to locate the physical surfaces on the mesh. Their sign indicate the side on which a point is located. Level sets use node-valued functions and are inter- polated with the basis functions of the finite element. This description allows to release the underlying mesh from the
description of surfaces of discontinuity or external boundaries.

As previously described, the X-FEM allows one to treat problems showing discontinuities or heterogeneities. Some works have been done for dynamic problems relative to crack propagations: different techniques of enrichment have been developed to take into account the discontinuity produced by a crack. Some discontinuous enrichments can be used such as in [5] for time-dependent problems. Other authors purpose a generalization of X-FEM by introducing enrichment strategies for time-dependent problems [16]. The discontinuities can also be enriched in time as in [15,7]. In this paper, we restrict our work to the treatment of holes and external boundaries. The X-FEM approach is easier to manage in this case since no specific enrichment functions are necessary. A direct use of enrichment functions is problematic for explicit schemes because a lumping technique is not easily available for such approximation. A work in this direction for cracks may be found in [11].

In the first part of this paper, the explicit direct integration method is first summarized: we recall the definition of the critical time step and the definition of stability and dispersion properties. Then we present the X-FEM for static problems. The second part of the paper introduces a 1D X-FEM finite element: we propose here a theoretical approach to treat the presence of void volume within this 1D finite element as well as a specific technique of mass lumping. Results from this simple 1D model are then generalized for triangular (2D) and tetrahedral (3D) finite elements. In the last part of this work, some numerical simulations in 1D, 2D and 3D cases are carried out. They are related to pierced structures under impact loads and are compared to simulations carried out on conforming meshes with the ABAQUS software. We show several results about convergence and stability of the technique. The conclusions drawn from these results finally guides the outlook concerning the X-FEM for rapid dynamics.

## 2. Foreword about the explicit Newmark scheme and the X-FEM

### 2.1. Explicit direct integration method

The following set of relations (1) gives the equations of dynamics, under matrix form, for a solid body at the finite element level

$$M\ddot{U} + KU = F^{\text{ext}}. \qquad (1)$$

The most general approach to solve (1) in time consists in using direct integration methods because the size of the system is usually very large. Among many techniques, the Newmark scheme [3,4,13] is commonly used for rapid dynamics: for each time step $\Delta t$, the accelerations (2) are computed and then allow to update the displacement and velocity fields

$$\ddot{U}_{n+1} = M^{-1} \left( F_{n+1}^{\text{ext}} - K \left[ U_n + \Delta t \dot{U}_n + \frac{\Delta t^2}{2} \ddot{U}_n \right] \right). \quad (2)$$

Eq. (2) shows that the explicit Newmark method possesses the advantage in the fact that during the resolution of equilibrium equations, only the mass matrix requires to be inverted. Furthermore, it is possible to apply mass lumping techniques which render this mass matrix diagonal. This leads to a very efficient numerical scheme.

However the explicit Newmark scheme also possesses the disadvantage of the necessity to impose a critical time step in order to maintain its stability. This constraint is also called CFL (Freidrich-Lévy-Courant) condition [3,4,9,13].

Stability is insured when the norm of the greatest eigenvalue of system is strictly inferior to one. Rather than identifying this previous eigenvalue which is quite time consuming, one can use the highest eigenvalue of each finite element taken separately. It is easier and much faster to find the smallest finite element within the most disadvantageous material to impose the critical time step. It can be shown [9] that this time step (3) constitutes a lower bound of the one issued from the complete structure, and that therefore respects the CFL condition. It is also a good approximation if the mesh is almost uniform

$$\Delta t \le \frac{2}{\omega_{\max}} = \Delta t_{C}. \qquad (3)$$

This reveals a potential problem with the explicit approach: if the discretized structure has just one very small finite element, the latter will constrain the critical time step for the simulation of the whole structure. That is the reason why some freedom is offered to the users of commercial software in choosing an “appropriate” time step (possibly bigger) with stabilization techniques. Nevertheless, the choice of the time step is generally limited by instability as upper bound and by the phenomenon of dispersion as lower bound:

* Instability occurs if the time step is greater than the critical time step. Overestimated forces are transmitted and reflected within some finite elements. The process repeats itself at each time step and it quickly leads to the divergence of the numerical scheme.

* Dispersion occurs if the time step is much lower than the critical time step. Small disturbances are generated and they are transmitted with a wave speed higher than physical wave speed. It implies an attenuation of the real efforts. It should be noticed that even if the quality of the solution is somewhat degraded, it is generally accurate. The smallest dispersion is however obtained for a time step equal to the critical time step.

### 2.2. The X-FEM

Within the finite element method the presence of discontinuities or heterogeneities such as cracks, inclusions or material interfaces constrains the mesh to be in conformity
<figure><img src="image_1.png" /><figcaption>Fig. 1. Example of conforming and non-conforming meshes to model an inclusion: (a) conforming mesh classical FEM; (b) non-conforming mesh X-FEM.</figcaption></figure>

with the geometry of these entities (Fig. 1a). In this case, surface management has to be completely integrated into the meshing process and can be very time consuming, specifically for the three-dimensional simulations.

For a couple of years now, a new approach has been developed: the eXtended Finite Element Method. Its main characteristic is to separate the previous problem into two parts:

* The first part corresponds to the discretization of whole domain, which does not include some or all the surfaces related to discontinuities or boundary conditions. The approximation of the displacement field is consequently the classical approximation used in FEM, and mesh generation is generally straightforward.

* The second part consists in supplementary shape functions added to some nodes of the previous approximation. The goal of these additional functions is to enrich the basic approximation of the existing displacement field with less regular functions able to model jumps on the surfaces which are not meshed. In Fig. 1b, black circles indicate these nodes. The determination of these nodes can be made with the help of a level set representation of the surface [17]. When nodes have been selected, specific enrichment functions are then associated to these and offer for instance the possibility to take the deformation discontinuity at a material interface into account.

Due to the Partition of Unity method [10], the approximation of the displacement field, described in Fig. 1b, follows Eq. (4), where N represents all the nodes of the global domain and D represents the nodes belonging to elements which are cut by the material interface

$$u(x) = \sum_{i \in N} u_i \varphi_i(x) + \sum_{i \in D} a_i \varphi_i(x) \psi(x). \quad (4)$$

The additional degrees of freedom generally have no point wise significance, but they allow one to model discontinu-

ities. This is an important point because lumping techniques for the mass matrix are commonly used to get a more efficient resolution. These techniques are based on some physical representation of phenomena: a straight application of enrichment functions is problematic for explicit schemes because a lumping technique is not directly available for such approximation. In this paper, we restrict our work to the treatment of holes and external boundaries. The X-FEM approach is easier to manage in this case since no specific enrichment functions are necessary. The linear and bilinear weak forms are simply integrated on the matter part of the elements [8,12].

## 3. X-FEM elements for explicit dynamic in 1D, 2D and 3D cases

### 3.1. The X-FEM bar element

Fig. 2 shows the elementary finite element bar: it is composed of two domains. The first one of length $\epsilon l_e$ contains matter while the second one of length $(1 - \epsilon)/l_e$ corresponds to an empty space. The total length of the finite element is $l_e$.

A classical sub-element integration gives the consistent stiffness and mass matrices (5). In order to respect the nature of each domain constituting the finite element, integration is only done inside the matter, i.e. using a range going from 0 to $\epsilon l_e$

$$ \begin{align} (K_{e4}) &= \frac{\epsilon ES}{l_e} \begin{bmatrix} 1 & -1 \\ -1 & 1 \end{bmatrix}; \notag \\ (M_{e4}) &= \frac{1}{3} \rho S \epsilon l_e \begin{bmatrix} \epsilon^2 - 3\epsilon + 3 & -\epsilon(2\epsilon - 3)/2 \\ -\epsilon(2\epsilon - 3)/2 & \epsilon^2 \end{bmatrix}. \tag{5} \end{align} $$

As said previously, the explicit direct integration generally uses a lumped mass matrix. Indeed, for systems with numerous degrees of freedom, lumping the mass matrix leads to an efficient algorithm inverting the matrix is straightforward. Among all techniques of lumping, the most common consists in replacing the diagonal matrix
<figure><img src="image_2.png" /><figcaption>Fig. 2. Description of the X-FEM bar element.</figcaption></figure>

terms by the sum of all terms from the corresponding row. That is equivalent to consider that the total mass of element is equally concentrated on the finite element nodes. By applying this approach to the mass matrix from (5), one obtains

$$ (K_e) = \frac{\varepsilon ES}{l_e} \begin{bmatrix} 1 & -1 \\ -1 & 1 \end{bmatrix}; \quad (M_e) = \frac{1}{2} \rho S \varepsilon l_e \begin{bmatrix} 2 - \varepsilon & 0 \\ 0 & \varepsilon \end{bmatrix}. \quad (6) $$

The critical time step associated to this system which satisfies the condition of Freidrich-Lévy-Courant, is therefore

$$ \Delta t_C \le \frac{2}{\omega_{\max}} = \frac{l_e \sqrt{\varepsilon(2 - \varepsilon)}}{c}. \quad (7) $$

Eq. (7) shows that the critical time step depends on $\varepsilon$. If $\varepsilon$ tends to zero, the critical time step also tends to zero making the algorithm unpractical. One idea to solve this problem is to avoid distributing the mass to nodes according to the domain length, but to equally distribute it on each node regardless of $\varepsilon$. For instance, if the total mass of the material domain is $\rho\varepsilon Sl_e$, each node receives one half of this mass; then the X-FEM mass and stiffness matrices are given by

$$ (K_e) = \frac{\varepsilon ES}{l_e} \begin{bmatrix} 1 & -1 \\ -1 & 1 \end{bmatrix}; \quad (M_e) = \frac{1}{2} \rho S \varepsilon l_e \begin{bmatrix} 1 & 0 \\ 0 & 1 \end{bmatrix}. \quad (8) $$

The critical time step associated to this X-FEM bar element is then

$$ \Delta t_C \le \frac{2}{\omega_{\max}} = \frac{l_e}{c}. \quad (9) $$

The difference between (8) and (6) lies in the fact that in the X-FEM context, summing the elements of each row does not equally distribute the mass over the nodes.

The latter lumping technique yields a critical time step independent of the material fraction within the element. This fact hints at a possible smoothing of critical time steps within a structure. For instance, we suppose that the mesh is quite regular except at a boundary where they are small finite elements. Regular finite elements possess a mass and stiffness matrices defined by (10) and a critical time step defined by (11). The rest of finite elements located near the boundary induces a smaller critical time step because of their size. Classical simulation will then use this time step that causes a high computational cost. One possible solution is to replace these small finite elements with X-FEM ones. In this case, the critical time step can be smoothly

computed and this fact can offer an additional flexibility concerning the numerical simulation

$$ (K_{e1,2,3}) = \frac{ES}{l_e} \begin{bmatrix} 1 & -1 \\ -1 & 1 \end{bmatrix}; \quad (M_{e1,2,3}) = \frac{1}{2} \rho S l_e \begin{bmatrix} 1 & 0 \\ 0 & 1 \end{bmatrix}, \quad (10) $$

$$ \Delta t_C \le \frac{2}{\omega_{\max}} = \frac{l_e}{c}. \quad (11) $$

To fully describe the X-FEM bar element, we precise how to take loads into account. One simple way is to distribute the force as in Fig. 3, using interpolation functions (12). This distribution will be justified in the numerical applications

$$ \left\{ \begin{array}{l} \text{node 1: } N1(\varepsilon) * F = (1 - \varepsilon)F, \\ \text{node 2: } N2(\varepsilon) * F = \varepsilon F. \end{array} \right. \quad (12) $$

As the characteristics of the X-FEM element are defined, we can now focus on the convergence of the explicit Newmark method. The main objective is to prove that the use of this element provides at least results that are similar to those of the FEM.

Therefore, lets take the dynamic problem of a bar which is clamped on one side and submitted to a step load $F$ on the other side. The bar possesses a length $L$, a Young's modulus $E$ and a section $S$. This problem admits an analytical solution given by (13) where $x$ denotes the distance to the clamped point and $t$ the current time

$$ u(x, t) = \frac{F}{ES} x + \frac{8L}{\pi^2} \sum_{i=1}^{\infty} \frac{(-1)^i}{(2i-1)^2} \sin\left(\frac{(2i-1)\pi x}{L}\right) \times \cos\left(\frac{(2i-1)\pi ct}{L}\right). \quad (13) $$

<figure><img src="image_3.png" /><figcaption>Fig. 3. Load distribution for the X-FEM bar element.</figcaption></figure>
In the case of a classical numerical simulation with a conforming regular mesh, the explicit Newmark method gives the exact solution for every nodes and at every discrete time step. However, if the solution is required at an intermediate time or position, an error appears due to the interpolation that is used. This error is summarized in Fig. 4: for a given critical time step, displacements issued from numerical simulation are exact at both nodes. Owing to the interpolation functions of finite element, the displacement field is linearly approximated between nodes. To treat the approximation in time (dashed curve), we decided to also use a linear interpolation because the critical time step is linked to mesh size by material celerity (11). Let us now calculate the L2 norm error (14) between the approximation (dashed curve) and the analytical solution (solid curve)

$$err = \sqrt{\int_{0}^{T} \int_{0}^{L} (U_{\text{analytical}} - U_{\text{calculated}})^2 dx dt.} \quad (14)$$

The analytical expression of the two curves are

$$\left\{ \begin{array}{ll} U_{\text{sim}} = \frac{F_c}{E_s} t & \text{for } 0 \le t \le \frac{h}{c}, \\ U_{\text{ref}} = \begin{cases} 0 & \text{for } 0 \le t \le \frac{\alpha h}{c}, \\ \frac{LF}{h(\alpha-1)ES} [-ct + \alpha h] & \text{for } \frac{\alpha h}{c} < t \le \frac{h}{c}. \end{cases} & (15) \end{array} \right. $$

The displacement error is then

$$\begin{aligned} err &= \sqrt{\frac{F^2 h (\alpha^2 L - \alpha L^2 + L \alpha h + L^2 + h^2 - 2Lh)}{3cE^2 S^2}} \\ &\approx \frac{Fh^{3/2}}{3\sqrt{cE_s}} + \dots \end{aligned} \quad (16)$$

This result shows that the order of convergence for the L2-norm in displacement, in the case of a bar impacted by a step load, is 3/2. This result is due to the fact that the impacted bar solution is not smooth. Indeed, for more general cases, under smoothness conditions, the order of convergence is 2 [2,6,18]. We keep in mind that the convergence results shape those obtained for a uniform mesh and a time step equal to the critical time step.

Fig. 5 represents some convergence curves coming from numerical simulations using the X-FEM element. The material ratio $\varepsilon$ is respectively 1%, 50% and 99%. The X-FEM finite element utilization does not modify the convergence rate of 3/2.

## 3.2. Generalization to 2D and 3D constant strain finite elements

The extension of the technique for 2D and 3D constant strain elements is easily carried out. Indeed, in the 1D case,

<figure><img src="image_4.png" /><figcaption>Fig. 4. Exact and numerical displacements over a time step.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Fig. 5. Convergence curves for different ratio of matter.</figcaption></figure>
$ε$ represents a ratio between the length occupied by the matter and the total length of bar. In the 2D case, this quantity will represent a ratio between the material surface and the total surface of the finite element. The same is true for the 3D case with volumes.

As for the X-FEM bar element, the stiffness matrix is obtained by consistent integration on the material domain using interpolation functions of classical finite elements. As for the mass matrix, it is built by the adapted lumping method, i.e. each node will receive an identical quantity of the mass issued from the material domain. Thus, just like in the 1D case, the critical time step does not depend on the material fraction within finite elements but only on the dimensions of the X-FEM finite element, as in standard FEM.

To illustrate the generalization of the method, we propose to identify the eigenvalue of a triangular T3 FEM and X-FEM finite element (Fig. 6). Both finite elements are defined in plane stresses and they possess a unit Young's modulus, a unit thickness, a Poisson's ratio equal to 0.5 and a unit density. The interpolation functions in reduced coordinates of the element are given in Eq. (17). In this example, the material fraction is equal to 8%

$$ \left\{ \begin{array}{l} N1(s,t) = 1 - s - t, \\ N2(s,t) = s, \\ N2(s,t) = t. \end{array} \right. \qquad (17) $$

After calculating the FEM and X-FEM stiffness matrices, one can see the proportionality factor of 12.5 (i.e. material fraction of 8%) between these two matrices

$$ Ke_{T3 \times FEM} = \frac{Ke_{T3}}{12.5} = \frac{1}{6 * 12.5} \begin{bmatrix} 5 & 3 & -4 & -1 & -1 & -2 \\ & 5 & -2 & -1 & -1 & -4 \\ & & 4 & 0 & 0 & 2 \\ & & & 1 & 1 & 0 \\ & & & & 1 & 0 \\ & & & & & 4 \end{bmatrix}_{\text{sym}}. \qquad (18) $$

The special lumping technique gives the following mass matrix:

<figure><img src="image_6.png" /><figcaption>Fig. 6. T3 FEM and X-FEM finite elements.</figcaption></figure>

$$ Ke_{T3 \times FEM} = \frac{0.04}{3} \begin{bmatrix} 1 & 0 & 0 & 0 & 0 & 0 \\ & 1 & 0 & 0 & 0 & 0 \\ & & 1 & 0 & 0 & 0 \\ & & & 1 & 0 & 0 \\ & & & & 1 & 0 \\ & & & & & 1 \end{bmatrix}_{\text{sym}}. \qquad (19) $$

As a consequence, the eigenvalue of this T3 X-FEM finite element, which will be used to satisfy the CFL condition, is 3.6457 and is identical to the pulsation of the T3 FEM finite element.

To conclude it is worth noting that stiffness and mass matrices for the X-FEM finite elements can be directly obtained from the “full material” regular FEM matrices, just multiplying them by the factor $ε$; this fact is always true for linear interpolation on simplex elements even in the case of nonlinear constitutive laws.

## 4. Numerical examples

### 4.1. Numerical simulation of an impact on a bar

This section deals with the numerical simulation of an impacted bar shown in Fig. 7. The length of the bar is denoted $L$. The left side, located at $x = 0$, is clamped and the other side, located at $x = L$, is subjected to a step load $F$. The bar is homogeneous and is composed of an isotropic elastic material. The bar is discretized using four finite elements; each of length $l_e$. This mesh is not conforming on the right side of bar ($x = L$): indeed, the last finite element only contains material over a length $εl_e$.

The material is steel for which Young's modulus $E$ is equal to $2.1 \times 10^{11}$ Pa and the density $ρ$ is equal to $7800$ kg/m$^3$. The bar's length $L$ is $1.5$ m and its section $S$ is $0.1$ m$^2$. The simulation time lasts $1.5$ ms. The impact load ($F$) on the right side of bar is $10,000$ N. The last finite element has a material fraction $ε$ of $10\%$.

Four types of simulations (Fig. 8) are carried out:

* Case 1: A fine uniform conforming discretization of the bar is done. It respects the physical geometry of bar.

* Case 2: The mesh is not uniform but still conforming. The critical time step is given by the shortest element.

<figure><img src="image_7.png" /></figure>
<figure><img src="image_8.png" /><figcaption>Fig. 7. Description of the impact of bar.</figcaption></figure>

* Case 3: The mesh is uniform but not conforming. The right finite element implements the new developed X-FEM technique. Therefore, the critical time step is the same for all finite elements.

* Case 4: An artificial mass [1] is added on nodes of the last finite element so that the critical time step of case 4 is identical to the case 3. This additional mass allows one to increase the critical time step of the last element to the value of the other elements (this simulation's purpose is informative).

Figs. 9 and 10 show numerical results for the displacements at points located at $x = L$ (after interpolation of displacements) and at $x = L - \varepsilon L_e$. These results are compared to the analytical solution.

The first conclusion with regard to these figures are very satisfactory: X-FEM results agree closely to the analytical results. It establishes that the developments of the X-FEM

<figure><img src="image_9.png" /><figcaption>Fig. 8. Different numerical simulations carried out.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Fig. 9. Displacement of the node, located at $x = (1 - \varepsilon)L$ for the bar under impact.</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>Fig. 10. Displacement of the point, located at $x = L$ for the bar under impact.</figcaption></figure>

are successful. Moreover, the behavior of the X-FEM is quite similar to the one that of a classical FEM numerical simulation (case 2). For both, the approximation is acceptable. However, it can be seen that the mass scaling option does not yield satisfactory results. Note that more evolved mass scaling strategies exist in the literature [14] and might not yield such unphysical results.

Table 1 summarizes some calculation characteristics. It points out that X-FEM finite element also allows one to decrease strongly the number of time step while offering solution in agreement with the analytical solution.

To investigate the good approximation properties of the new finite element, a set of numerical simulation were done for $\epsilon$ varying from 0 to 1. A comparison between numerical and analytical results is made using the L2 norm error defined by (14). The numerical solution is interpolated linearly in time between the computed instants. The results are plotted in Fig. 11. It shows that the error remains acceptable when compared to the analytical solution.

In addition, Fig. 11 shows also another interesting aspect: the dispersion is weaker with the new finite element than with classical FEM, and this is particularly true for very small value of $\epsilon$. Indeed, when the mesh is conforming, very small elements are needed. This induces small time step in order to respect the CFL condition. This small time step increases the dispersion in the FEM results. Therefore, if there exists small elements, which can be seen as distur-

bance of the discretization of the structure, the use of the X-FEM finite element may yield better results than conventional FEM.

## 4.2. Numerical simulation of dynamic loading of pierced plates

This example deals with the impact of a plate which is clamped at one side and subjected to a sudden pressure of 100,000 Pa, on its other side, which remains constant after. The dimensions of the plate are 500 mm × 250 mm × 2 mm. The plate possesses three holes (50 mm in diameter) located respectively at 125 mm, 250 mm and 375 mm in X-coordinate and at 125 mm in Y-coordinate. The constitutive material is steel ($E = 210,000$ MPa, $\nu = 0.3$, $\rho = 7800$ kg/m$^3$). The total time of the simulation is 5 ms in order to observe stress waves propagating over several periods.

Two triangular meshes were used (under the assumption of plane stress). The first mesh is in conformity with the geometry and uses standard T3 FEM elements while the second one uses T3 X-FEM finite elements and does not conform to the three holes (Fig. 12). For both meshes, the average element size is about 15 mm. The simulation for the conforming elements is made with ABAQUS software.

Fig. 13 summarizes the vertical displacements observed for nodes denoted “hole” and “border” (Fig. 12). It shows

<table><caption>Table 1
Summary of the results for the bar under impact</caption><thead><tr><th></th><th>Case 1: Uniform fine conforming mesh</th><th>Case 2: Non-uniform conforming mesh</th><th>Case 3: Non-conforming mesh – X-FEM option</th><th>Case 4: Conforming mesh – mass scaling option</th></tr></thead><tbody><tr><td>Number of elements</td><td>31</td><td>4</td><td>4</td><td>4</td></tr><tr><td>Time step</td><td>9.325395 ms</td><td>9.325395 ms</td><td>93.25395 ms</td><td>93.25395 ms</td></tr><tr><td>Number of time steps</td><td>161</td><td>161</td><td>16</td><td>16</td></tr><tr><td>Maximum displacement (x = L)</td><td>1.41 µm</td><td>1.46 µm</td><td>1.38 µm</td><td>0.36 µm</td></tr></tbody></table>
<figure><img src="image_12.png" /><figcaption>Fig. 11. L2 displacement error norm as a function of the mesh size and epsilon.</figcaption></figure>

that the results are very close for the FEM and the X-FEM results. Indeed, for the “hole” node as well as for the “border” node, the displacements agree very well for both simulations.

In addition, the critical time step identified by ABAQUS code is 1.00 ms and it is induced by a small finite element around the upper hole. In the X-FEM code, the critical time step is 1.60 ms. Releasing the mesh from constraints is shown here to improve the critical time step and thus increases performance.

As previous results seem to be reliable, we aim to study the convergence of the simulation. Therefore, we focus again on the dynamic problem of a plate clamped at one side and subjected, on the other side, to a sudden pressure of 100,000 Pa. The dimensions of the plate are now of 50 mm × 50 mm × 2 mm. The plate, in this case, possesses

<figure><img src="image_13.png" /><figcaption>Fig. 12. Modeling of a plate subjected to a dynamic compression.</figcaption></figure>

only one hole of variable diameter (from 15 to 30 mm by step of 0.5 mm). The hole is located at the center of plate. The constitutive material is steel ($E = 210,000$ MPa, $v = 0.3$, $\rho = 7800$ kg/m$^3$). The total time of the simulation is also equal to 0.5 ms in order to observe the stress waves propagating over several periods.

Two triangular meshes were used under the assumption of plane stress. The first mesh (Fig. 14) is in conformity with the geometry and uses standard T3 FEM elements while the second one (Fig. 15) uses T3 X-FEM finite elements. For both meshes, the size of the elements varies from one time the thickness (2 mm) to five times the thickness by steps of 0.5 mm.

The variable size of the hole and the element size aim at studying the behavior T3 X-FEM elements compared to classical FEM elements. Results are collected for a point located at $x = 50$ mm and $y = 25$ mm and relate to hori-

<figure><img src="image_14.png" /><figcaption>Fig. 13. Comparisons between ABAQUS and X-FEM codes for the displacement field at two different locations.</figcaption></figure>
<figure><img src="image_15.png" /><figcaption>Fig. 14. Conforming mesh for the plate with one hole.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>Fig. 15. Non-conforming mesh for the plate with one hole.</figcaption></figure>

zontal displacements. Fig. 16 shows, for a given case (size
of the hole of 17.5 mm), that results about displacements
agree well with explicit classical FEM computations.

In order to investigate the behavior of the new
approach, we computed the L2 norm (16) for each simula-
tion, taking as reference a classical FEM calculation made
with a very refined mesh. Fig. 17 shows that the X-FEM
method in rapid dynamics does not modify the overall
numerical behavior when compared to classical FEM.
The trends are nearly identical. A remark, on the other
hand, is to be raised concerning the absolute value of the
error: for a sufficiently discretized mesh, the error is less
in the X-FEM case.

As the study relative to the displacement fields leads to
good conclusions, we now focus on stresses within the
plate. The results in this section are presented for a plate
with a hole of 17.5 mm of diameter. The mesh size varies
from 1 mm to 5 mm for both FEM and X-FEM cases.
Moreover, in order to have a good comparison, we impose
to the mesh processing to create two “fixed” finite elements
on specific locations. The first one is depending on the mat-
ter ratio that is included in the range from 0.01% to
99.99%. It belongs to the internal boundary, i.e. it is near
the hole (location around the point of coordinates
x = 25 mm ; y = 33.75 mm). The second finite element is
located near the external boundary (near the point of coor-
dinates x = 25 mm and y = 50 mm). This latter possesses
characteristic length imposed by mesh size.

For these two elements and for each numerical simula-
tion, we extract the von Mises stress along time. To com-
pare their evolution, we take as reference the results
coming from refined uniform meshes which are based on
the finite element size near the hole. It must be noted that
we cannot compute a reference for the case of 0.1% and
0.01% of matter due to the too large number of degree of
freedom.

Fig. 18 summarizes, for the finite element located near
the hole, the evolution of the stresses according to the ratio
of matter. The stresses are collected at the final time of sim-
ulation which is 25 µs. We represent on the graph only the

<figure><img src="image_17.png" /><figcaption>Fig. 16. Comparisons of the displacement results for different meshes.</figcaption></figure>
<figure><img src="image_18.png" /><figcaption>Fig. 17. Error evolution for conforming and non-conforming meshes.</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>Fig. 18. Evolution of the longitudinal stress for an element near the hole for different mesh sizes and matter ratio e.</figcaption></figure>

curves coming from the sizes of mesh of 1, 3 and 5 mm in order to clarify the figure. The first conclusion is that FEM and X-FEM behaviors are quite similar: indeed, the stress of von Mises increases as the mesh size decreases. However, that raises the following problem: if a level set imposes that an element is cut very near to a node, the stress will be overestimated within this X-FEM finite element. Nevertheless, it should be noticed that compared to FEM simulation, the time step used will be much less important with X-FEM: as comparison, for the numerical simulation where the mesh size is about 1 mm and for which the matter ratio is 0.01%, the time step used by X-FEM is approximately $5.75 \times 10^{-8}$ s whereas it is $9.6 \times 10^{-10}$ s for FEM.

Another important advantage is shown in Figs. 19 and 20: the curves represent the evolution along time of von Mises stresses in the two finite elements (near the hole and on external boundary) for FEM, X-FEM and reference simulations. The results come from the case where

the size of mesh is 1 mm and where the percentage of matter is 0.01%. Fig. 19 again shows that evolution of stress near the hole is similar in FEM and X-FEM cases. Fig. 20 shows on the other hand that von Mises stress in the outside element (which is a “standard” finite element in both cases) is better evaluated for X-FEM than for FEM: indeed, in this last case, the time step employed is imposed by the small finite element located near the hole. As we said previously, this last one is approximately hundred times smaller than for X-FEM case. Some phenomena of dispersion are then observed which slightly disturb the FEM results. This fact is of course not observed for X-FEM since we use the critical time step, which is close to one for all finite elements.

In conclusion, the use of X-FEM finite elements gets the double advantage of not being dependent on the mesh generation and on using a time step more “homogeneous”. Even if one sometimes observes some finite elements for
<figure><img src="image_20.png" /><figcaption>Fig. 19. Comparisons of the longitudinal stresses in an element near the hole for FEM and X-FEM.</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>Fig. 20. Comparisons of the longitudinal stresses in an element near the external boundary for FEM and X-FEM.</figcaption></figure>

which the ratio of matter is small and which could induce an abusive increase of stresses, the preceding advantages mask this disadvantage.

### 4.3. Numerical simulation of a diabolo under dynamic loading

In order to show the relevance of the approach for 3D cases, we propose in this part a simulation on a “diabolo” subjected to dynamic loading. The results obtained with the X-FEM code are compared at various locations of the part with the results coming from the ABAQUS code.

The part (Fig. 21) consists of a unit cube whose matter has been extracted in a torus whose center is located at the center of cube. The ring’s radius is equal to half-diagonal (0.5¹/² m) of the cube and the width radius of ring section is equal to 0.3 m.

The material characteristics are those of steel ($E = 2.1 \times 10^{11}$ Pa, $\rho = 7800$ kg/m³ and $v = 0.3$). The boundary

<figure><img src="image_22.png" /><figcaption>Fig. 21. Diabolo geometry.</figcaption></figure>
<figure><img src="image_23.png" /><figcaption>Fig. 22. ABAQUS mesh of the diabolo.</figcaption></figure>

<figure><img src="image_24.png" /><figcaption>Fig. 23. X-FEM mesh for the diabolo (6 elements per edge).</figcaption></figure>

conditions are a complete clamping of the lower face and a
sudden traction of 10 MPa applied on the upper face
(Fig. 21).

<figure><img src="image_25.png" /><figcaption>Fig. 24. Iso-zero level set used in X-FEM code to represent the boundary.</figcaption></figure>

The mesh consists of tetrahedral elements. For the con-
forming mesh shown in Fig. 22 (used in ABAQUS code),
there are 10 elements along the edges of the model resulting
in 5063 elements. For X-FEM simulations, a coarser mesh
(Fig. 23) is used and has only 6 elements along the edges,
resulting in 2456 finite elements. For the X-FEM mesh,
the geometry of the diabolo is represented by a level set
(Fig. 24).

The comparisons between displacements in X-FEM and
the ABAQUS codes are carried out for points located at
the following coordinates: (1, 0, 1) (point no. 1) and (–.289,
0, .289) (point no. 2).

Figs. 25 and 26 show the results for the displacement at
points 1 and 2. The comparison shows that the results com-
ing from both simulations are equivalent: it should be
noted here that the mesh complexity is higher than in pre-
vious 2D cases, because the volume of “empty space” is
about half that of the initial cube.

The time steps used in calculations are respectively
8.77 ms for the X-FEM (6 elements per edge) and
8.27 ms for ABAQUS (10 elements per edge). An X-
FEM computation with a coarser mesh, which offers a time

<figure><img src="image_26.png" /><figcaption>Fig. 25. X-FEM/ABAQUS comparison for the point no. 1 of measurement.</figcaption></figure>
<figure><img src="image_27.png" /><figcaption>Fig. 26. X-FEM/ABAQUS comparison for the point no. 2 of measurement.</figcaption></figure>

step equivalent to that of an ABAQUS computation, allows one to obtain a solution of very good quality.

## 5. Conclusions

We exposed in this paper the work carried out in connection with an X-FEM type finite elements in its application to rapid dynamics. One of the major problems in this case relates to the critical time step, which is controlled by both the smallest characteristic size of finite element and the most unfavorable material celerity. It is possible to identify two cases which can strongly penalize an explicit numerical simulation: the first case is related to a material heterogeneity within the structure and second case is related to mesh constraints imposed by a precise description of defects, interfaces or holes. The latter is considered in this study: avoiding constraints in mesh generation for complex models.

We first propose a 1D bar X-FEM finite element. The objective is to be able to model the behavior of an element composed of both material and empty space. Thanks to the adapted lumping technique, we could get back a critical time step that is identical to that of a completely filled finite element. We then proposed the generalization of this technique to 2D and 3D constant strain finite elements. This generalization is done in a straightforward manner: as in the 1D case, one represents the percentage of material according to the ratio of surfaces (in 2D) or volumes (in 3D). The technique of adapted lumping proposed in the 1D case remains valid.

Next, some numerical results are shown. The results of numerical and analytical comparisons in the 1D case demonstrate that errors remain acceptable when compared to the present existing solutions (mass scaling, imposition of small time steps, etc.). Moreover, these results also highlight another interesting aspect: the dispersion with the X-FEM approach is smaller than that of the classic FEM approach, especially when the material fraction $\varepsilon$ tends to zero. Numerical simulations of a perforated plate subjected

to compression or tensile loads show the quality of the results obtained by the X-FEM when compared to the classic FEM in the 2D case. The displacements that are measured in various points of the plate are coherent and relative errors remain small when compared to reference cases. Finally, we propose a 3D simulation of a diabolo. As for the 1D and 2D cases, the results agree well with the ABAQUS software.

The outlooks of this study is to extend the technique for quadrangular and hexahedron finite elements. In this case, the interpolation functions are not linear and there are some potential issues with the adapted lumping technique. It should be noted that quadrangular and hexahedral finite elements are very popular among structural analysts.

## Acknowledgements

We would like to thank PSA Peugeot-Citroën S.A. for their financial support and more particularly MM. Laurent ROTA and Malek ZARROUG.

## References

[1] ABAQUS/Explicit User's Manual, ABAQUS 6.5, 2005.

[2] T. Belytschko, N. Holmes, R. Mullen, Explicit integration—stability, solution properties, cost, ASME AMD 14 (1975) 1–21.

[3] T. Belytschko, J.L. Lin, Explicit algorithms for the nonlinear dynamics of shells, Comput. Methods Appl. Mech. Engrg. 42 (1984) 225–251.

[4] T. Belytschko, T.J.R. Hughes, Computational Methods for Transient Analysis, North-Holland, 1986.

[5] T. Belytschko, H. Chen, X. Jingxiao, Z. Goangseup, Dynamic crack propagation based on loss of hyperbolicity and a new discontinuous enrichment, Int. J. Numer. Methods Engrg. 58 (2003) 1873–1905.

[6] H.C. Chan, C.W. Cai, Y.K. Cheung, Convergence studies of dynamic analysis by using the finite element method with lumped mass matrix, J. Sound Vib. 165 (2) (1993) 193–207.

[7] J. Chessa, T. Belytschko, Arbitrary discontinuities in space–time finite elements by level sets and X-FEM, Int. J. Numer. Methods Engrg. 61 (2004) 2595–2614.

[8] C. Daux, N. Moës, J. Dolbow, N. Sukumar, T. Belytschko, Arbitrary branched and intersecting cracks with the eXtended Finite Element Method, Int. J. Numer. Methods Engrg. 48 (2000) 1741–1760.
[9] T.J.R. Hughes, J.S. Pister, R.L. Taylor, Implicit-explicit finite elements in non-linear transient analysis, Comput. Methods Appl. Mech. Engrg. 17/18 (1979) 159-182.

[10] J.M. Melenk, I. Babuska, The partition of unity finite element method: basic theory and applications, Comput. Methods Appl. Mech. Engrg. 139 (1996) 289-314.

[11] T. Menouillard, J. Réthoré, A. Combescure, H. Bung, Efficient explicit time stepping for the eXtended Finite Element Method (X-FEM), Int. J. Numer. Methods Engrg. 69 (2006) 911-939.

[12] N. Moës, M. Cloirec, P. Cartraud, J.F. Remacle, A computational approach to handle complex microstructure geometries, Comput. Methods Appl. Mech. Engrg. 192 (2003) 3163-3177.

[13] M.N. Newmark, A method of computation for structural dynamics, ASCE Journal of the Engineering Mechanics Division 85 (EM3) (1959).

[14] L. Olovsson, K. Simonsson, M. Unosson, Selective mass scaling for explicit finite element analyses, Int. J. Numer. Methods Engrg. 63 (10) (2005) 1436-1445.

[15] J. Réthoré, A. Gravouil, A. Combescure, An energy-conserving scheme for dynamic crack growth using the eXtended finite element method, Int. J. Numer. Methods Engrg. 63 (5) (2005) 631-659.

[16] J. Réthoré, A. Gravouil, A. Combescure, A combined space-time extended finite element method, Int. J. Numer. Methods Engrg. 64 (2) (2005) 260-284.

[17] N. Sukumar, D.L. Chopp, N. Moës, T. Belytschko, Modeling holes and inclusions by level sets in the extended finite element method, Comput. Methods Appl. Mech. Engrg. 190 (2001) 6183-6200.

[18] S.R. Wu, Lumped mass matrix in explicit finite element method for transient dynamics of elasticity, Comput. Methods Appl. Mech. Engrg. 195 (44-47) (2006) 5983-5994.