# X-FEM a good candidate for energy conservation in simulation of brittle dynamic crack propagation

Alain Combescure, Anthony Gravouil, David Grégoire, Julien Réthoré

► To cite this version:

Alain Combescure, Anthony Gravouil, David Grégoire, Julien Réthoré. X-FEM a good candidate for energy conservation in simulation of brittle dynamic crack propagation. Computer Methods in Applied Mechanics and Engineering, 2008, 197 (5), pp.309-318. 10.1016/j.cma.2007.04.011. hal-00381473

HAL Id: hal-00381473

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

Submitted on 27 May 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

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

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# X-FEM a good candidate for energy conservation in simulation of brittle dynamic crack propagation

A. Combescure a,*, A. Gravouil a, D. Grégoire a, J. Réthoré b

a LaMCoS, INSA-Lyon, CNRS UMR5259, F69621, France

b Aerospace Engineering Department, Tu-Delft, Netherlands

This paper is devoted to the simulation of dynamic brittle crack propagation in an isotropic medium. It focuses on cases where the crack deviates from a straight-line trajectory and goes through stop-and-restart stages. Our argument is that standard methods such as element deletion or remeshing, although easy to use and implement, are not robust tools for this type of simulation essentially because they do not enable one to assess local energy conservation. Standard cohesive zone models behave much better when the crack's path is known in advance, but are difficult to use when the crack's path is unknown. The simplest method which consists in placing the cohesive segments along the sides of the finite elements leads to crack trajectories which are mesh-sensitive. The adaptive cohesive element for-mulation, which adds new cohesive elements when the crack propagates, is shown to have the proper energy conservation properties during remeshing. We show that the X-FEM is a good candidate for the simulation of complex dynamic crack propagation. A two-dimensional version of the proposed X-FEM approach is validated against dynamic experiments on a brittle isotropic plate.

Keywords: Dynamic crack propagation; X-FEM; Cohesive zone; Remeshing; Simulation; Crack arrest

## 1. Introduction

The calculation of dynamic crack propagation remains a difficult challenge. Many contributions have been made on this topic. For the mechanical part, one may mention the works of Freund [14], Bui [19] and Lemaitre [15]. For the computational aspects, many authors have also addressed the problem using different methods, such as local smeared cracking, which relies on material models which include damage [7,6,10], or cohesive zone models, which are clearly related to fracture mechanics concepts and have been proven effective for localized fracture Falk [30–32]. Cohesive zone models have been used extensively, especially in cases where the crack's trajectory is known in advance, and more recently have also been extended to adaptive calculations in which cohesive elements are inserted into the mesh progressively

as the crack travels or branches [33–35,28,29]. X-FEM simulation of dynamic crack propagation was first presented by Krysl and Belytschko [25]. The present paper focuses on the comparison of standard FEM dynamic crack propagation simulation with X-FEM simulation. Experimental results are used to assess the validity of the calculations. First, the paper presents the computational models commonly used for crack propagation. We introduce the global theory of dynamic rupture, based on the evaluation of stress intensity factors, followed by the local approach to rupture. Next, we present three usual calculation strategies for the simulation of dynamic crack propagation: element deletion, remeshing, and the use of cohesive zone elements. Then, we briefly present the X-FEM formulation and compare it to the other methods using the same DCB example. We explain the good quality of the dynamic crack propagations obtained using the X-FEM by applying the conservation of energy principle and proving mathematically that the X-FEM method

<footnote>* Corresponding author. Tel.: +33 4 72 43 64 26.
E-mail address: alain.combescure@insa-lyon.fr (A. Combescure).</footnote>
guarantees exact energy conservation when the crack propagates. This proof is also valid for the adaptive cohesive zone formulation of dynamic fracture problems using constant strain finite elements. Finally, we apply the X-FEM to the prediction of a crack's propagation in a simple experiment involving a complex crack path with kinks and a stop-and-restart history.

## 2. Mechanical modeling of dynamic crack propagation

### 2.1. Global and local approach of rupture

Even though the simulation of dynamic brittle crack propagation remains a difficult challenge, the underlying physical fracture mechanics model is relatively simple and based on three key concepts [19]:

(1) an equation which gives the crack propagation direction;

(2) a criterion for the initiation of crack propagation;

(3) an equation which gives the crack's velocity.

Usually, two approaches are in competition for this type of prediction: a global energy approach, which is often preferred for brittle rupture, and an approach based on local stresses, the latter more effective for ductile fracture. We will limit ourselves to crack propagation driven by the maximum hoop tensile stress alone. Let us recall the main features of the two approaches.

#### 2.1.1. The global energy approach to rupture

Brittle crack propagation is assumed to be governed by the maximum value of the hoop stress $\sigma_{00}$ near the crack's tip, which is evaluated using the stress intensity factor $k_{00}$:

$$k_{00} = \lim_{r \to 0} \sqrt{2\pi r \sigma_{00}}, \quad (1)$$

where $(r, \theta)$ are the local polar coordinates of the crack's tip.

The maximum hoop stress intensity factor and the corresponding local polar angle are denoted $K^*$ and $\theta^*$ respectively:

$$K^* = \max_{\theta \in [-\pi, \pi]} k_{00} = k_{0^*0^*} \quad (2)$$

The propagation begins when the maximum hoop stress intensity factor is greater than a critical value called the dynamic crack initiation toughness. The direction of propagation is that of the maximum hoop stress [27]. This criterion can be written as follows:

$$ \begin{aligned} K^* < K_{1d} \quad (\text{no initiation}), \\ K^* = K_{1d}, \quad \theta^* = \theta_c \quad (\text{initiation}). \end{aligned} \quad (3) $$

The dynamic crack initiation toughness is a material property which is obtained from experiments.

During the dynamic growth of the crack, the velocity of the crack's tip $\dot{a}$ adjusts itself so that the current maximum

hoop stress intensity factor $K^*$ remains equal to the dynamic crack growth toughness:

$$K^*(t, \dot{a}) \ge K_{1d} \Rightarrow K^*(t, \dot{a}) = K_{1D}(\dot{a}) \quad (\text{propagation}). \quad (4)$$

The evaluation of $K_{1D}(\dot{a})$ was given by Kanninen [26], who replaced the quasi-static toughness by the dynamic crack initiation toughness. Then, the dynamic crack growth toughness is assumed to be

$$K_{1D}(\dot{a}) = \frac{K_{1d}}{1 - \left(\frac{\dot{a}}{c_R}\right)}. \quad (5)$$

In Eq. (5), $c_R$ is the velocity of the Rayleigh waves. Bui [19] calculated the propagation direction $\theta^*$ analytically through the following equation:

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

The corresponding $K^*$ value is

$$K^* = \cos^3 \left( \frac{\theta^*}{2} \right) K_1^{\text{dyn}} - \frac{3}{2} \cos \left( \frac{\theta^*}{2} \right) \sin(\theta^*) K_2^{\text{dyn}}. \quad (7)$$

These physical crack propagation laws provide the key to the prediction of the change in the crack's length at each time step of a transient analysis as proposed by Bui, Freund, or Tuler [19,14,21]. These laws can also be used, when the crack is meshed explicitly or using cohesive zone models, to detect whether it should be remeshed. The calculation of the dynamic stress intensity factors $K_1^{\text{dyn}}$ and $K_2^{\text{dyn}}$ is necessary for this approach.

#### 2.1.2. Approach based on the local stress and damage

In local fracture models, one dismisses the previous concepts and uses an equivalent method based on the stress and damage fields. The stress and damage states at the crack's tip [15] define how the crack progresses:

(1) the crack propagates if the maximum hoop stress at the crack's tip is greater than a critical value or if damage reaches the critical value;

(2) the crack propagates in the direction of the maximum hoop stress;

(3) in general, the crack propagation velocity is not controlled.

However, this common method based on a local vision of fracture is mesh-dependent: the finer the mesh, the faster the crack's propagation. In order to model brittle failure, one introduces a simple elastic softening failure law.

Then, one calculates the principal stresses at all the integration points of all the elements. If one of these principal tensile stresses exceeds the failure stress $\sigma_c$, damage starts to grow. A typical stress-strain curve is shown in Fig. 1, where one can also observe the effect of damage on the
<figure><img src="image_2.png" /><figcaption>Fig. 1. Typical stress–strain curve for a brittle material.</figcaption></figure>

elastic properties. The Young's modulus decreases as soon as the strain is larger than $\epsilon_c$ and reaches zero when the strain is greater than $\epsilon_R$. As one can see, the model leads to softening of the material: this softening is known to induce a mesh dependant response if no special treatment is introduced such as non-local damage or a control of damage rate.

### 3. Numerical modeling of rupture

#### 3.1. Remeshing

An alternative approach consists in meshing the crack at its initial location, using the criterion given in Section 2.1.1 to predict the crack's new location, then remeshing the volume of the body with the crack's new geometry and performing another time step with this new geometry. This method seems rather simple, but presents a number of significant drawbacks. Let us consider the cracked structure at time $n$. Let $X_n^n$ and $X_{n+1}^{n+1}$ designate the geometry and state variables respectively before remeshing and at time $n+1$ after remeshing. In the process, the crack's length has increased by $\Delta a$.

One can define the state $X_n^{n+1}$, which describes the state at time $n$ on the new geometry (increased crack's length and new mesh). A first approximation of this state is obtained by a “projection” operation of state field $X_n^n$ onto the mesh at time $n+1$. This state must be in equilibrium with the applied loads, and the increased crack segment $\Delta a$ must be closed. This closing is achieved by applying the appropriate closing force to the extended crack's lips. Equilibrium iterations are necessary in order to achieve equilibrium of the projected field. To go from State $X_n^{n+1}$ to State $X_{n+1}^{n+1}$, one must increase the applied loads while simultaneously releasing the closing loads. These operations are very involved and require significant CPU time. There are several possible strategies for releasing the crack's closing loads, and the choice of a strategy influences the robustness of the method. It has been shown [24] that if

one uses a Newmark time integration scheme this type of approach cannot preserve local energy. The Newmark scheme is governed by the following equations:

$$ \underline{u}_{n+1} = \underline{u}_n + \Delta t \dot{\underline{u}}_n + \frac{\Delta t^2}{2} [(1 - 2\beta)\ddot{\underline{u}}_n + 2\beta \ddot{\underline{u}}_{n+1}], \quad (8) $$

$$ \dot{\underline{u}}_{n+1} = \dot{\underline{u}}_n + \Delta t [(1 - \gamma)\ddot{\underline{u}}_n + \gamma \ddot{\underline{u}}_{n+1}], $$

where $\beta$, $\gamma$ are the two integration parameters of the scheme, and $\underline{u}$, $\dot{\underline{u}}$, $\ddot{\underline{u}}$ designate the displacements, velocities and accelerations respectively. The subscript $n$ denotes the time step number, and $\Delta t$ is the time step. Let us now introduce the following notations:

$$ [X] = X_{n+1} - X_n, $$

$$ \langle X \rangle = \frac{1}{2} (X_n + X_{n+1}). $$

If one expresses the energy variation between state $n$ with Mesh $n$ and state $n+1$ with Mesh $n+1$ and uses Eq. (8), one gets [24]:

$$ \begin{aligned} & \frac{1}{2} [\dot{u}_{n+1}^{n+1T} M_{n+1} \dot{u}_{n+1}^{n+1} + \underline{u}_{n+1}^{n+1T} K_{n+1} \underline{u}_{n+1}^{n+1}] - \frac{1}{2} [\underline{u}_{n}^{nT} M_{n} \underline{u}_{n}^{n} + \underline{u}_{n}^{nT} K_{n} \underline{u}_{n}^{n}] \\ & = \left(\gamma - \frac{1}{2}\right) [\underline{u}]^T [\underline{F}]^{\text{ext}} + [\underline{u}]^T \langle \underline{F}^{\text{ext}} \rangle - \frac{\Delta t^2}{2} (2\beta - \gamma) [\ddot{\underline{u}}]^T M_{n+1} \langle \ddot{\underline{u}} \rangle \\ & \quad - \frac{\Delta t^2}{2} \left(\gamma - \frac{1}{2}\right) (2\beta - \gamma) [\ddot{\underline{u}}]^T M_{n+1} [\ddot{\underline{u}}] - \left(\gamma - \frac{1}{2}\right) [\underline{u}]^T K_{n+1} [\underline{u}] \\ & \quad - (\gamma - 1) [\underline{u}]^T ([M] \ddot{\underline{u}}_{n}^{n+1} + [K] \underline{u}_{n}^{n+1}). \end{aligned} \quad (9) $$

In Eq. (9), one finds the standard terms of common Newmark algorithms plus an additional last term $((\gamma - 1)[\underline{u}]^T([M]\ddot{\underline{u}}_n^{n+1} + [K]\underline{u}_n^{n+1}))$: this term is due to the change of discretization at time $n$. It cannot be proved to be zero in the general case. It can only be calculated, to check its relative value in each numerical experiment. This is the intrinsic reason why energy conservation cannot be ensured through remeshing. The remeshing operation itself is delicate: it can either introduce energy into the calculation (which leads to an unstable computation scheme), or dissipate energy in an uncontrolled manner. This remark also holds for the simulation of evolving cracks with remeshing. The use of remeshing in 3D problems is not very realistic because of the lack of robust remeshing tools for 3D situations.

#### 3.2. Local damage models

The material softening behavior is known to make the calculation dependent on the mesh size. The predicted failure load gets smaller as the mesh size decreases. It has been shown analytically that, for a one-dimensional bar subjected to a tensile step load, damage is always localized in the smallest element. More generally, the failure prediction is mesh-dependent. Non-local models enable such mesh dependency to be avoided. Two classes of models are currently available, mainly for quasi-static cases: non-local models, which have been widely used, and second-gradient
models, which require the resolution of an auxiliary problem. The two approaches have been applied successfully to many static 2D problems. Their extension to 3D situations is possible, but would involve high computational costs because the size of the elements must be small. Many authors have tried to overcome this difficulty, mainly by using static concepts such as non-local measures of damage or second-gradient theories. Others have introduced delayed damage models, which are essentially based on the idea that damage cannot extend at an infinite speed. The approach first introduced by Seaman and Curran [6] in their dynamic fracture constitutive model consists in limiting the damage growth rate. A similar idea, called “delayed damage” modeling, was proposed later [10,2]. Limiting the damage growth rate leads, as is the case in non-local damage models, to a converging solution when the mesh size decreases [11]. This was proved mathematically and numerically in Suffis [16,17]. The damage rate is bounded; consequently, when the mesh size tends to zero, localization cannot take place in a single element alone and the numerical localization effect disappears. An equivalence between geometrical non-local models and delayed damage models was proposed in Suffis [17]. The main advantages of this local approach to fracture are that it enables both initiation and propagation to be predicted indifferently, that it can be easily extended to ductile fracture, and that in a more basic sense fracture mechanics is not needed in order to predict failure: all the information is contained in the stress-strain representation of the material. One must also mention that one can establish a connection with the classical fracture mechanics energy release rate $G_c$ if one notes that the area of the triangle $\sigma_c \epsilon_c \epsilon_R$ is the released energy density of a fully damaged Gauss point in an element. This energy is simply transferred to the rest of the structure when the element disappears. Nevertheless, one can observe that the crack’s path and velocity are obviously mesh-dependent, and that this dependence decreases as the mesh size decreases Suffis [18]. Today, there is no model capable of defining the velocity of the crack’s tip. The local models rely on some simple basic ideas. One makes a fine mesh of the structure being considered, which may or may not present an initial crack. These models require very refined meshes in order to achieve correct crack propagation, and are certainly more appropriate for ductile fracture [7]. In practice, these models are not very powerful for extensive failure propagation in 3D dynamic situations.

### 3.3. X-FEM modeling of dynamic crack propagation

In this section, we will show that the X-FEM used appropriately to calculate the dynamic stress intensity factors provides a good model for the resolution of dynamic crack propagation problems. We will also show that exact energy conservation can be proved and achieved with this method, even when the crack propagates. Finally, we will

compare the predictions using the X-FEM and the FEM on the DCB specimen.

#### 3.3.1. Calculation of the stress intensity factors

The accurate calculation of the stress intensity factors $K_1^{dyn}$ and $K_2^{dyn}$ is an important requirement for this approach. The most efficient method consists in using interaction integrals, which are an efficient means of computing these quantities with results which are not too mesh-dependent as long as the calculated quantities are not too close to the crack’s tip. The interaction integral $I$ is

$$ I = - \int_{\Omega} q_{k,j} \left( \sigma_{ml}^{aux} u_{m,l} \delta_{jk} - \left[ \sigma_{ij}^{aux} u_{j,k} + \sigma_{ij} u_{j,k}^{aux} \right] \right) d\Omega + 2 \int_{\Omega} q_k \left( \sigma_{ij,j}^{aux} u_{i,k} + \rho \left[ \ddot{u}_i u_{j,k}^{aux} + \dot{u}_i^{aux} \dot{u}_{i,k} + \dot{u}_i \dot{u}_{i,k}^{aux} \right] \right) d\Omega. \quad (10) $$

In Eq. (10), $q$ is the virtual extension field, which must be tangent to the crack’s surface and kinematically admissible; $u$ and $\bar{u}$ are the displacement and stress fields of the cracked body; the superscript aux designates a virtual auxiliary displacement field and the corresponding stress field. One can prove that

$$ I = \frac{2}{E^*} \left( K_1^{dyn} K_1^{aux} + K_2^{dyn} K_2^{aux} \right). \quad (11) $$

In Eq. (11), $E^*$ is the Young’s modulus in the plane stress case and is equal to $E/(1 - v^2)$ in the plane strain case. If one chooses the known analytical Mode-1 or Mode-2 fields as the auxiliary displacement and stress fields, Eq. (11) enables the direct calculation of the stress intensity factors $K_1^{dyn}$ and $K_2^{dyn}$. The surface integral evaluation of $I$ [14] is generally preferred to the standard $J$-integral approach [20] because it is more compatible with the usual finite element integration schemes. The best evaluation is obtained with a $q$ field whose spatial derivatives are zero close to the crack’s tip [13]. This enables one to avoid the calculation of the $K_i^{dyn}$ in the vicinity of the crack’s tip, where the stress evaluation is rather poor.

#### 3.3.2. The X-FEM and energy conservation

The method is now well-known, having been described by many authors, such as Black [22], Belytschko [23], and Moes [5,8]. The crack is represented independently of the structure, either by using a discretized mesh (2-node segments in 2D or 3-node triangles in 3D) or by level set representation [9]. The latter constitutes the main practical advantage of the method, first because it enables one to introduce new cracks without remeshing the structure (and, therefore, leads to a very easy parametric analysis of cracked structures) and, second, because it enables one to propagate the cracks without remeshing the structure, which makes the crack propagation analysis automatic and easy. Usually, one uses the structural mesh to represent the level set values, but the simple explicit algorithms used to move the level set when the discontinuity evolves require
a fine, regular mesh of the finite-difference type. Hence, in the case of an unstructured mesh, one observes level set distortions. A simple, pragmatic way of avoiding this problem consists, as explained in Prabel [1], in using a specific, fine and regular mesh (distinct from that used for structural calculations and limited to the regions where the crack might propagate) to represent the level set values. The usual continuous displacement field is enriched with two sets of functions: discontinuous Heavyside-type functions for the nodes of the elements completely traversed by the crack, and singular enrichment for the nodes of the elements only partially cut.

The approximate displacement field in an element is

$$ \underline{u} = \sum_{i=1,\text{nodes}} N_i \underline{u}_i + \sum_{k=1,p} \Phi_k \underline{b}_{ik} \). \quad (12) $$

In this equation, $N_i$ denotes the standard shape function, $\underline{u}_i$ is the continuous displacement vector and $\underline{b}_{ik}$ is the value of the kth enrichment function at node i. When the crack traverses the element entirely, the enrichment function $\Phi_k$ can be chosen as the Heavyside function, while when the element is only partially traversed (i.e. the crack's tip lies somewhere within the element) the singularity is defined by the following four enrichment functions:

$$ \{\Phi_k(r, \theta)\}_{k \in \{1,4\}} = \sqrt{r} \left\{ \sin\frac{\theta}{2}, \cos\frac{\theta}{2}, \sin\frac{\theta}{2} \sin\theta, \sin\frac{\theta}{2} \cos\theta \right\}. \quad (13) $$

In Eq. (13), r designates the distance to the crack's tip and θ designates the angle between the point and the tangent at the crack's tip. This second part of the model makes the method very efficient and attractive because one introduces the analytical singular solutions in the vicinity of the crack's tip, thus suppressing the need for a very refined mesh near the crack's tip. This method has been used in the context of dynamic crack propagation [4]. Crack propagation is achieved either by adding a new segment to the representation of the crack or by updating the crack's level set representation and adding new degrees of freedom (because the set of elements traversed by the crack and the set of elements at the crack's tip both change).

### 3.3.3. Proof of energy conservation

Let us now assume, as before, that the crack is at its original position at time n. The mass matrix is denoted $M_n$, the stiffness $K_n$, the velocity field $\dot{u}_n = \dot{u}_{old}$, and the displacement $\underline{u}_n = \underline{u}_{old}$. The crack extends by $\Delta a$ and reaches its final position at time $n+1$. If one adds new degrees of freedom corresponding to the new position of the crack's tip, the velocity field $\dot{u}$ and the displacement field $\underline{u}$ are extended by two velocity and displacement fields as shown in the following equations:

$$ \begin{aligned} \underline{u}_{n+1} &= \left\{ \begin{array}{l} \underline{u}_{old} \\ \underline{u}_{new} \end{array} \right\}, \\ \dot{\underline{u}}_{n+1} &= \left\{ \begin{array}{l} \dot{\underline{u}}_{old} \\ \dot{\underline{u}}_{new} \end{array} \right\}. \end{aligned} $$

Consequently, the new mass matrix is

$$ M_{n+1} = \begin{bmatrix} M_n & \widetilde{M}_{no} \\ \widetilde{M}_{no}^T & \widetilde{M}_{new} \end{bmatrix}. \quad (14) $$

A similar equation holds for the stiffness matrix.

$$ K_{n+1} = \begin{bmatrix} K_n & \widetilde{K}_{no} \\ \widetilde{K}_{no}^T & \widetilde{K}_{new} \end{bmatrix}. \quad (15) $$

The following extension strategy is chosen: the crack is moved by $\Delta a$ at time n. New (extended) degrees of freedom are added at time n and initialized to zero. One obtains:

$$ \dot{u}_{n}^{n+1} = \left\{ \begin{array}{l} \dot{u}_{n}^{n} \\ 0 \end{array} \right\}, \quad (16) $$

$$ u_{n}^{n+1} = \left\{ \begin{array}{l} u_{n}^{n} \\ 0 \end{array} \right\}. \quad (17) $$

Then, the change in kinetic energy when one moves the crack at time n is

$$ 2[\mathbf{W}^{\text{kinetic}}] = \dot{u}_{n}^{n+1\mathrm{T}} M_{n+1} \dot{u}_{n}^{n+1} - \dot{u}_{n}^{n\mathrm{T}} M_{n} \dot{u}_{n}^{n}. \quad (18) $$

The change in potential energy for the extended crack at time n is

$$ 2[\mathbf{W}^{\text{def}}] = u_{n}^{n+1\mathrm{T}} K_{n+1} u_{n}^{n+1} - u_{n}^{n\mathrm{T}} K_{n} u_{n}^{n}. \quad (19) $$

First, let us observe that with this initialization choice the additional segment $\Delta a$ closes at time n. The forces applied to close the new crack segment are given by Eq. (20):

$$ \underline{F}_{n}^{n+1 \text { close }}=\left[\widetilde{K}_{\text { no }}^{\mathrm{T}}\right] \underline{u}_{\text { old }}. \quad (20) $$

Introducing Eqs. (14) and (16) into Eq. (18), one gets no change in kinetic energy. The same result holds for the potential energy by introducing Eq. (17) into Eq. (19). The work of the external loads at this instant is also unchanged. Thus, one can observe that the proposed strategy is such that the change of discretization at time n preserves the energy exactly, both locally and globally, even if the “mesh” is changed at this instant. It is impossible to prove this property for other types of “adaptive” meshes. Nevertheless, it is possible to set constraints on the adaptive meshing process so that the global energy remains constant. This, however, does not guarantee energy conservation at the crack's tip and makes the progress of the crack's tip difficult to control. Conversely, if one constrains the energy to remain in the vicinity of the crack's tip, it is very difficult to enforce conservation over the whole structure. In all cases, this conservation property is only approximate, whereas in the proposed X-FEM approach energy conservation holds exactly.
Then, Newmark's standard scheme is used to update time from $n$ to $n + 1$. During this time increment, the crack extends by $\Delta a$ and the closing forces 20 at time step $n$ are automatically released to zero at time step $n + 1$. The X-FEM formulation chosen, which includes the singular fields, ensures that the work of the closing forces during opening corresponds exactly to the energy release rate [4].

## 3.4. Cohesive zone models

### 3.4.1. Description of the family of models

The cohesive zone model is a good candidate for the simulation of crack propagation, even in 3D dynamics, when the crack's path is known *a priori*. Then, the method is very simple: one meshes the future crack path, even in 3D, using a set of cohesive zone elements. These elements are initially closed except at the location of the initial crack, and they progressively open as the (dynamic) loading is applied [12]. In these elements, one uses a cohesive law model similar to that of Fig. 2. The cohesive zone model relates the displacement jump across the interface to the traction at the interface. Again, with this model, there is no law giving the velocity of the crack's tip. The most stable solutions are obtained with $[u_n]_c = 0$.

If the crack's path is not known *a priori*, which is the case when one wants to predict the crack's propagation for a future experiment, the standard method for practical industrial applications is the following: one puts cohesive zone elements at all element side segments (or, in 3D, surfaces) connected to the crack's tip. Then, the crack's tip can progress. Nevertheless, it is natural to consider this method to be mesh-dependent, especially when the crack does not propagate in a straight line, and it has been very frequently observed that the final shape of the cracked structure depends on the chosen mesh. Therefore, the method should be used with an adaptive strategy which consists in inserting cohesive elements which “cut” the old mesh: thus, during

the simulation of a crack's propagation in the general case, the mesh changes with time.

### 3.4.2. energy conservation in case of adaptative cohesive zone elements

The same type of approach and demonstration also applies to the conservation of the kinetic energy when one adds a cohesive zone element, provided that the added nodes have no mass. In this case, the additional mass matrices $\tilde{M}_{\text{no}}$ and $\tilde{M}_{\text{new}}$ are null matrices. Now, let us make the following three assumptions:

(a) the basic finite elements are constant strain elements only,

(b) the mesh created after updating the cohesive zone elements contains constant strain elements only (e.g. in 2D, 3-node constant strain triangles),

(c) the displacements and velocities of the added nodes at time step $n$ are interpolated within the old elements through the shape functions.

Then, the strain energy is the same before and after remeshing at time $n$. This property is obvious because the strain is constant within an element and if this element is subdivided into any number of constant strain elements these inherit the same strain as the father element. The total strain energy is unchanged. Moreover, it is also obvious that if the material is history-dependent (for instance in the case of plasticity) and the internal variables as well as the stress states are copied to the Gauss points of the new elements, the irreversible energy is also preserved through this remeshing strategy.

Then, the implicit Newmark scheme can be applied as usual between time steps $n$ and $n + 1$. The cohesive zone formulation ensures that the energy released in the propagation is the fracture energy release rate.

If the mass matrix changes within the element cut by the new cohesive segment, the kinetic energy distribution within this element changes slightly, but the total kinetic energy within the element can be kept constant with appropriate velocity initialization at the added nodes. Energy conservation at time $n$ is then still preserved.

## 3.5. Application to the DCB virtual test

Let us now consider the example of a DCB specimen, which is rather simple because the crack's propagation follows a straight line and presents no stop-and-restart effects. The geometry is described in Fig. 3. The specimen is loaded by two opposite $\pm 0.0025$ m vertical step displacements at Nodes A and A'. The material's data are: Young's modulus 186 GPa, Poisson's ratio 0.3, specific mass 8000 kg/m$^3$ and fracture toughness $K_{\text{IC}} = 100$ MPa $\sqrt{m}$.

This example was calculated in plane strain using the remeshing method and cohesive zone modeling. The same calculation was also performed with the X-FEM: this will be discussed in the next section. A relatively coarse mesh

<figure><img src="image_3.png" /><figcaption>Fig. 2. Cohesive zone stress–strain law.</figcaption></figure>
<figure><img src="image_4.png" /><figcaption>Fig. 3. DCB geometry.</figcaption></figure>

was used in order to emphasize the difficulties associated with the remeshing phase of the operation. The initial mesh is shown in Fig. 4.

The remeshing strategy was applied using two different algorithms: the first consisted simply of remeshing and transferring the fields to the new geometry: this is the basic method in usual explicit dynamic codes, which generally ignore equilibrium iterations. The second strategy consisted in re-equilibrating the fields transferred to the new mesh and, in the same equilibrium iterations, applying forces to the crack's extension so that it closes at the beginning of the time step. This example was also calculated using the X-FEM. The mesh consisted of 400 4-node elements. As can be observed in Fig. 5, which compares the crack's evolutions with time given by 4 different approaches, the crack progresses much too fast when the

<figure><img src="image_5.png" /><figcaption>Fig. 4. Mesh of the DCB specimen.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Fig. 5. Compared crack evolution: I – FEM without equilibrium control; II – FEM with equilibrium iterations, III – X-FEM and cohesive elements.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Fig. 6. Energy conservation check on the DCB specimen test.</figcaption></figure>

equilibrium step is ignored during remeshing (Curve I). When using the remeshing algorithms, the crack ends up splitting the specimen into two parts, whereas when using a cohesive model or the X-FEM energy conserving scheme the crack stops, which is consistent with experimental observation.

The energy conservation check is presented in Fig. 6. An exact calculation should produce no energy. A huge amount of energy is introduced into the calculation if one does not re-equilibrate the structure after remeshing (FEM curve), but this input is drastically reduced if one re-equilibrates the structure (FEM+ curve). The cohesive zone model and the X-FEM simulation have much better energy conservation properties.

## 4. An original experiment

The objective of this section is to present the validation of the X-FEM in reference to an original dynamic crack growth experiment and to show that the X-FEM provides a good prediction of a crack's path and velocity and can even predict the interesting crack stop-and-restart phenomenon which is observed experimentally.

### 4.1. *Description of the experiment and main results*

The proposed experiment, described in Fig. 7, is based on a Hopkinson bar impact system. The specimen's material is PMMA and the cylindrical bars are made of Nylon in order for them to have the right impedance.

Fig. 8 shows the geometry of the test specimen, whose thickness is 15 mm.

The material properties for Nylon and PMMA are given in Table 1.

The experiments were performed on the Hopkinson bar system of the LMS Laboratory at Ecole Polytechnique, France. They are highly reproducible. A high speed camera was used to follow the crack's tip. Only four pictures could
<figure><img src="image_8.png" /><figcaption>Fig. 7. Hopkinson test setup.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Fig. 8. Geometry of the PMMA specimen tested.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Fig. 9. Experimental horizontal crack propagation history.</figcaption></figure>

<table><caption>Table 1
Material properties for the experimental test</caption><thead><tr><th>Property</th><th>Name</th><th>Unit</th><th>PMMA</th><th>Nylon</th></tr></thead><tbody><tr><td>Static Young's modulus</td><td><i>E</i></td><td>GPa</td><td>3.3</td><td>3.5</td></tr><tr><td>Dynamic Young's modulus</td><td><i>E<sub>D</sub></i></td><td>GPa</td><td>4.25</td><td>3.6</td></tr><tr><td>Poisson's ratio</td><td><i>ν</i></td><td>—</td><td>0.42</td><td>0.41</td></tr><tr><td>Density</td><td><i>ρ</i></td><td>kg m<sup>-3</sup></td><td>1180</td><td>1145</td></tr><tr><td>Rayleigh wave velocity</td><td><i>c<sub>R</sub></i></td><td>m s<sup>-1</sup></td><td>1064</td><td>996</td></tr><tr><td>Fracture toughness</td><td><i>K<sub>IC</sub></i></td><td>MPa √m</td><td>1.47</td><td>—</td></tr></tbody></table>

<figure><img src="image_11.png" /><figcaption>Fig. 10. The experimental crack's final path.</figcaption></figure>

be taken during the crack's propagation time. Reproducibility was used: by repeating the same experiment the desired number of times with different trigger times for the camera, we were able to catch about 20 different positions of the crack's tip. The tests were repeated five times for each experiment.

The evolution of the crack's length as a function of time is shown in Fig. 9.

One can observe a first propagation stage at 200 µs with a velocity close to 260 ms⁻¹. Then, the crack stopped for 50 µs, after which it restarted with a speed of 160 ms⁻¹ before coming to a final stop.

The crack's trajectory is shown in Fig. 10.

One can observe a first kink of about 38° at the beginning of the propagation, then a second kink after the first stop. This test enables one to test situations of dynamic crack propagation with mixed-mode loading and complex stop-and-restart. These experimental results are challenging for calculation tools. Yet, the test remains relatively simple because linear material behavior is sufficient to predict the response.

## 4.2. X-FEM predictions

Let us now present the X-FEM simulation of this experiment. A plane strain model was used. Had the crack's path been known, a cohesive zone approach could have been successful, but it could not be used in order to predict the crack's path prior to the experiment. A delicate point in setting up the X-FEM model is the modeling of the experimental Hopkinson bar system. One could mesh the two bars, but this would lead to long and useless calculations. We chose to model the test specimen alone. The mesh, consisting of 1500 4-node elements, is shown in Fig. 11. One can observe that the mesh also covers the hole. This is normal with the X-FEM, which is not required to follow the structure's boundaries exactly.
<figure><img src="image_12.png" /><figcaption>Fig. 11. Mesh and boundary conditions used for the simulation.</figcaption></figure>

The measured velocities were prescribed on the left-hand side of the model. They are shown in Fig. 12.

An impedance (or adsorbing) condition was prescribed on the part of the right-hand side of the specimen in contact with the nylon bar. The equation giving the associated load is

$$F_x(t) = -Z\dot{u}_x = -\sqrt{E_{nylon}\rho_{nylon}}\dot{u}_x.$$

The theory described in the previous section was applied. The calculation was carried out using Newmark's implicit mean acceleration scheme with 100 constant time steps of 5 µs each. The calculated crack's path is compared to the experimental path in Fig. 13: this result is rather good, which indicates that the angle model is reasonable and that static crack initiation values are sufficient to obtain a good approximation of the crack's trajectories. The crack's deviations (especially the two kinks after the stops) and the final length are well-reproduced. This indicates that the directional criterion given by Eq. (6) and the crack's initiation criterion 7 and values were correctly chosen.

Similarly, the calculated horizontal crack's velocity is compared to the measured results in Fig. 14 and the comparison is reasonable. One can observe the stop-and-restart behavior, but the calculated restart time is slightly different from the measured value. The calculation leads to the correct crack velocities in the two propagating phases as well as in the stop-and-restart: this inspires confidence in the equation proposed for the crack's propagation velocity.

<figure><img src="image_13.png" /><figcaption>Fig. 12. History of the prescribed velocity on the left-hand side of the model.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>Fig. 13. Comparison of calculated (continuous lines) vs. experimental (dotted lines) crack's paths.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>Fig. 14. Comparison of calculated vs. experimental horizontal crack tip velocities.</figcaption></figure>

The stopping phase is due to wave reflections within the specimen which “close” the crack's tip.

## 5. Conclusions

This paper showed that the X-FEM is an interesting tool for the simple and efficient simulation of dynamic crack propagation. Other methods based on local fracture mechanics require very refined meshes and have not yet been proven efficient and precise enough. Remeshing methods are not very robust. Cohesive zone models are interesting for the interpretation of known crack paths, but are difficult to use and not yet very robust for unknown crack paths and coarse meshes. Based on the cases presented here as well as other situations (such as those presented in Prabel [1] for 2D elastic-plastic simulation or in Menouillard
[3] for elastic crack propagation of other complex experiments using an explicit implementation of the X-FEM), the method clearly appears really promising and precise in predicting elastic and confined elastic–plastic dynamic crack propagation. The simple crack propagation laws based on an initiation, a directional criterion and a speed criterion seem sufficient to reproduce experimental observations. One can note that the results presented here concern only two-dimensional examples. The X-FEM seems to be a good candidate for a 3D extension of the proposed strategy to the case of elastic or confined plastic crack propagation. Of course, the computation time would be larger but, thanks to the limited mesh size constraints compared with other usual methods, one can reasonably expect that it would at least be possible to derive a good predictive tool based on 3D X-FEM. The determination of ad hoc, true 3D crack propagation criteria remains an open question, but the comparison of 3D X-FEM simulations with 3D dynamic crack propagation experiments should foster the development of interesting 3D dynamic crack propagation theories and models.

## Acknowledgements

This work was sponsored by ANR grant JC05_45254.

## References

[1] B. Prabel, A. Combescure, A. Gravouil, S. Marie, Level set X-FEM non-matching meshes: application to dynamic crack propagation in elastic-plastic media, Int. J. Numer. Methods Engrg. 69 (2007) 1533–1569.

[2] O. Allix, J.F. Deu, Delayed-damage modelling for fracture prediction of laminated composites under dynamic loading, Engrg. Trans. 45 (1) (1997) 29–46.

[3] T. Menouillard, J. Rethore, A. Combescure, H. Bung, Efficient explicit time stepping for the extended Finite Element Method (X-FEM), Int. J. Numer. Methods Engrg. 68 (9) (2006) 911–939.

[4] J. Rethore, A. Gravouil, A. Combescure, An energy conserving scheme for dynamic crack growth using extended finite element method, Int. J. Numer. Methods Engrg. 63 (2005) 631–659.

[5] N. Moes, J. Dolbow, T. Belytschko, A finite element method for crack growth without remeshing, Int. J. Numer. Methods Engrg. 46 (1) (1999) 133–150.

[6] D.R. Curran, L. Seaman, Dynamic failure of solids, Phys. Rep. 147 (1987) 253–388.

[7] A.L. Gurson, Continuum theory of ductile rupture by void nucleation and growth: Part I – Yield criteria and flow rules for porous ductile media, J. Mech. Phys. Solids 17 (1977) 201–217.

[8] N. Moes, A. Gravouil, T. Belytschko, Non-planar 3D crack growth by the extended finite element and level sets. Part I: Mechanical model, Int. J. Numer. Methods Engrg. 53 (11) (2002) 2549–2568.

[9] A. Gravouil, N. Moes, T. Belytschko, Non-planar 3D crack growth by the extended finite element and level sets—Part II: Level set update, Int. J. Numer. Methods Engrg. 53 (11) (2002) 2569–2586.

[10] P. Ladeveze, A damage computational method for composite structures, Comput. Struct. 44 (1992) 79–87.

[11] P. Ladeveze, O. Allix, J.F. Deu, A mesomodel for localization and damage computation in laminates, Comput. Methods Appl. Mech. Engrg. 183 (2000) 105–122.

[12] J. Remmers, R. de Borst, A. Needleman, A cohesive segments method for simulation of crack growth, Comput. Mech. 31 (2003) 69–77.

[13] P. Attigui, C. Petit, Mixed mode separation in dynamic fracture mechanics: new path independent integrals, Int. J. Fracture 84 (1) (1997) 19–36.

[14] L.B. Freund, Dynamic Fracture Mechanics, Cambridge Monographs on Mechanics and Applied Mathematics, 1990.

[15] J. Lemaitre, A Course on Damage Mechanics, Springer, Paris, 1996.

[16] A. Suffis, A. Combescure, Modèle d'endommagement à effet retard, étude numérique et analytique de l'évolution de la longueur caractéristique, Rev Européenne Eléments Finis 11 (5) (2002).

[17] A. Suffis, A. Combescure, A. Lubrecht, Damage model with delay effect: analytical and numerical studies of the evolution of the characteristics length, Int. J. Solids Struct. 40 (2003) 3463–3476.

[18] A. Suffis, Développement d'un modèle d'endommagement à taux de croissance contrôlée pour la simulation robuste de ruptures sous impact, Ph.D. thesis no. 04ISAL0033, 2004, INSA-Lyon.

[19] H.D. Bui, Mécanique de la rupture fragile, Masson, 1978.

[20] H.D. Bui, Associated path independent J integral for separation mixed modes, J. Mech. Phys. Solids 31 (1983) 439–448.

[21] F.R. Tuler, B.M. Butcher, A criterion for the time dependence of dynamic fracture, Int. J. Fracture 4 (1968) 431–437.

[22] T. Black, T. Belytschko, Elastic crack growth in finite element with minimal remeshing, Int. J. Numer. Methods Engrg. 45 (1999) 601–620.

[23] T. Belytschko, N. Moes, S. Usui, C. Parimi, Arbitrary discontinuities in finite elements, Int. J. Numer. Methods Engrg. 50 (4) (2001) 993–1013.

[24] J. Rethore, A. Gravouil, A. Combescure, A stable numerical scheme for the finite element simulation of dynamic crack propagation with remeshing, Comput. Methods Appl. Mech. Engrg. 193 (2004) 4493–4510.

[25] P. Krysl, T. Belytschko, Dynamic propagation of arbitrary 3D cracks, Int. J. Numer. Methods Engrg. 44 (6) (1999) 767–800.

[26] M. Kanninen, C.H. Popelar, Advanced Fracture Mechanics, Oxford University Press, 1985.

[27] H. Maigre, D. Rittel, Mixed-mode quantification for dynamic fracture initiation: application to the compact compression specimen, 1993, Int. J. Solids Struct. 30 (23) (1993) 3233–3244.

[28] F. Zhou, J.F. Molinari, Three dimensional dynamic crack propagation using cohesive elements: A methodology to address mesh dependency, Int. J. Numer. Methods Engrg. 59 (1) (2004) 1–24.

[29] F. Zhou, J.F. Molinari, T. Shioya, A rate-dependent cohesive model for dynamic crack propagation in brittle materials, Engrg. Fracture Mech. 72 (9) (2005) 1383–1410.

[30] X.P. Xu, A. Needleman, Numerical simulations of fast crack growth in brittle solids, J. Mech. Phys. Solids 42 (9) (1994) 1397–1407.

[31] G.T. Camacho, M. Ortiz, Computational modelling of impact damage in brittle materials, Int. J. Solids Struct. 33 (20) (1996) 2899–2938.

[32] M.L. Falk, A. Needleman, J.R. Rice, A critical evaluation of cohesive zone models of dynamic fracture, J. Phys. IV 5 (2001) 43–50.

[33] M. Ortiz, A. Pandolfi, Finite-deformation irreversible cohesive elements for three-dimensional crack-propagation analysis, Int. J. Numer. Methods Engrg. 44 (9) (1999) 1267–1282.

[34] A. Pandolfi, M. Ortiz, An efficient adaptive procedure for three-dimensional fragmentation simulations, Engrg. Comput. 18 (2002) 148–159.

[35] G. Ruiz, A. Pandolfi, M. Ortiz, Three-dimensional cohesive modeling of dynamic mixed-mode fracture, Int. J. Numer. Methods Engrg. 52 (2001) 97–120.