Efficient explicit time stepping for the eXtended Finite Element Method (X‐FEM)

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


INTRODUCTION
An important domain of application of the finite element method is the dynamic analysis of impacts, crashes or explosions. When analysing the propagation of waves in a structure, one must use relatively small time steps in order to represent the physics correctly. This is also the case when one uses explicit time integration. Generally, one must go through a very large number of time steps, which sometimes leads to high computation costs. In order to reduce the computation cost of a transient analysis, one can use an explicit time integration scheme (usually the central difference method [1]). Moreover, if a lumping technique is used for the mass matrix, the solution of the problem does not require the resolution of a linear system. One drawback of the central difference method, however, is that its stability is conditional and that the time step must satisfy the Courant's condition. For wave propagation studies, the stability condition is not really burdensome. This technique is used in a great many commercial codes dedicated to transient analyses.
Unfortunately, a major improvement to the finite element method, the partition of unity method [2], appears to be incompatible with the stability condition of explicit time integrators. The eXtended Finite Element Method [3][4][5] uses a local partition of unity to incorporate a discontinuity in the displacement or strain field into the interpolation. Thus, the mesh does not need to describe the geometrical support of this discontinuity. This method enables one to simulate the propagation of arbitrary cracks, even in three-dimensional applications [6][7][8].I t was shown in Reference [9] that the use of a discontinuous enrichment results in a dramatic reduction of the critical time step corresponding to the Courant's condition when the crack's surface is close to the nodes of the finite element mesh. This can be easily understood by looking at a one-dimensional problem: let us consider the one-dimensional element enriched with a discontinuous function shown in Figure 1. The equivalent finite element problem consists of two elements with a double node at the location of the discontinuity. For both problems, when the support of the discontinuity approaches Node 1 or Node 2 , the critical time step tends to zero. In the equivalent FE problem, the critical time step tends to zero because of the presence of a very large term in the stiffness matrix (∼ 1/l) and a very small term in the mass matrix (∼ l) of the smallest element. For the X-FEM problem, the critical time step tends to zero because the mass matrix is singular. Several solutions have been proposed in the literature in order to overcome his difficulty. In a previous paper [10], the use of an implicit time integrator was proposed, which is expensive in terms of computation cost. In Reference [11], the authors used an implicit-explicit time integrator in order to treat the region which contains the enriched functions with an unconditionally stable implicit scheme. Using this technique, the mass matrix is not diagonal and the same time step is used both in the implicit and in the explicit regions. In Reference [12], numerical simulations were carried out using a special algorithm which prevents the discontinuity from crossing the element in the vicinity of a node. In addition to this modification of the crack's path, a large reduction factor was applied to the standard critical time step in order to ensure the stability of the explicit time integrator. Using such a procedure, the time step is often very small in relation to the minimum distance authorized by the algorithm between the discontinuity and the nodes. The solutions mentioned above are not efficient enough for an extensive use of the partition of unity concept in explicit dynamic codes.
In the present paper, a mass matrix lumping technique is proposed for enriched functions. The idea of this lumping technique is based on an exact kinetic energy representation for basic motions. A proof of stability for the 1D and 2D cases is given. We show that if one uses this technique, the critical time step does not tend to zero as the support of the discontinuity reaches the boundaries of the elements. This enables one to use an explicit time integrator with a time step of the order of the critical time step of the problem without enrichment.
The paper is organized as follows: in Section 2, the continuous and discrete formulations of the problem are presented; Section 3 is dedicated to the lumping technique and the proof of stability; numerical examples of dynamic crack propagation are presented in Section 4.

Continuous formulation
Let us consider the homogeneous body with a crack described in Figure 2. The material has a mass density and is assumed to have linear elastic isotropic behaviour. This behaviour can be described using Lamé's constants , along with the Young's modulus E and Poisson's ratio . The motion of the body is described by the displacement u(x,t), x being the location of a material point and t the time. We assume small perturbations. The body is subjected to prescribed displacements u d at the boundary * 1 and to external loads f d in and F d at the boundary * 2 .( * 1 and * 2 are such that * 1 ∪ * 2 = * and * 1 ∩ * 2 = л.) n is the outward normal to the material boundary. The crack's faces are traction-free. The strong form of the problem can be written as follows: where is the Cauchy stress tensor, the symmetric first-order strain tensor and C the constitutive law.˙=*/*t denotes the time derivative. The weak form of these equations becomes: where u is in the space of the kinematically admissible functions and v is an arbitrary function kinematically admissible to 0.

Discrete formulation
To capture the discontinuity and the singularity in the strain field at the crack's tip, the partition of unity properties of the finite element shape functions [2] are used to enrich the spatial interpolation as follows: In this equation, u h is the approximate displacement field, N the set of nodes used to discretize , and N i the finite element shape function associated with node i.Ū i are classical degrees of freedom, whereasŨ ij are additional degrees of freedom supported by the nodes in the set N j associated with the enrichment function j . LettingŨ ij = 0, u h is a standard finite element approximation. Letting U i = 0,Ũ ij o = a andŨ ij = 0 for j = jo,w eh a v e : Consequently, the enriched approximation can describe the function exactly. The approximate displacement field u h is: where N contains the standard shape functions and the complete basis of shape functions. U andŨ j contain the standard and additional degrees of freedom, respectively, and U all the degrees of freedom.
Following Reference [3], the basis of enriched functions contains the generalized Heaviside function H and the branch functions: The nodes enriched with the discontinuous function are those whose support is completely cut by the crack, whereas the nodes enriched with the branch functions are those whose support contains the crack's tip. This enrichment will be denoted X cla .
In the following discussion, we will introduce the enrichment functions by using additional degrees of freedom for the discontinuous function alone. This type of enrichment will be denoted X dis . The affected nodes are the nodes of the elements cut by the crack (see Figure 3). Using this technique, the element containing the crack's tip is assumed to be completely cut by the crack, as in References [12,13]. This leads to some imprecision regarding the location of the asymptotic behaviour of the displacement field, but such errors become negligible as the mesh is refined. In References [11,14], an enrichment technique which enables the effective positioning of the crack's tip with a discontinuous enrichment alone is proposed. Using one of the enrichment techniques presented above, the discrete balance of linear momentum in space is obtained from Equation (6) by choosing u, v in the complete basis of shape functions : U,U andÜ denote, respectively, the displacement, velocity and acceleration vectors discretized on the complete basis of shape functions. F is the external force vector. M and K are the mass and stiffness matrices: Then, we introduce the matrix: (K 11 , K 12 , K 21 , K 22 ,): Thus, the stiffness matrix of an enriched element is When enriched interpolation is used, the consistent mass matrix has standard terms ( , and coupling terms ( N i N j k l d ).

Time integration.
The Newmark method is chosen as the time integrator. The updating equations are given by an d are the two parameters of the Newmark scheme. In order to study the stability properties, we follow the energy method from Reference [15]. Thus, we get back to the stability conditions of the Newmark scheme: unconditionally stable scheme 1 2 and 2 stable scheme if t 1 max is the largest solution of det(K − 2 M) = 0. For the central difference method, = 0 and = 0.5, the problem is: given U n ,U n andÜ n , find U n+1 ,U n+1 andÜ n+1 such that: with the stability condition Here, we can see that if the mass matrix is lumped, the resolution of the problem does not require the use of a linear system solver. Since the consistent mass matrix contains standard terms, pure enriched terms and coupling terms, the question becomes: how can one lump such a mass matrix?

LUMPING TECHNIQUE FOR THE MASS MATRIX
The objective of this section is to construct a lumped mass matrix which enables one to perform explicit dynamic analyses in the X-FEM framework with time steps of reasonable size. A basic requirement for this lumped mass matrix is that for rigid body motions, the discrete kinetic energy ( 1 2U T MU) be exact. We propose to define the diagonal coefficients m diag of the lumped mass matrix corresponding to the enriched degrees of freedom as follows: where el is the element being considered, m its mass, mes( el ) its length (in 1D), area (in 2D) or volume (in 3D), n nodes the number of nodes of el , and the enriched function.

Proof for the 1D case
Following this idea, let us consider a one-dimensional element with two nodes. Each node has ordinary degrees of freedom corresponding to the shape functions N 1 ,N 2 and additional degrees of freedom corresponding to the enriched function 1 . The approximate displacement is A lumped mass matrix for this element is of the form: We want to find the coefficients m i such that First, we consider a motion described by a constant velocityu =ā. Thus, we setU i toā anḋ U i1 to 0. In this case and where m is the mass of the element. Therefore, if we choose m 1 = m 2 = m/2 we get T h = T . This constitutes a practical means of lumping the mass matrix. Now, let us consider a motion described byu =ã 1 (x) (this represents the separation of the element into two parts). We seṫ U i to 0 andU i1 toã, and the energies are and Therefore, an exact representation of the kinetic energy (T h = T ) is obtained if the terms of the lumped mass matrix are calculated using: Using this lumped mass matrix, the kinetic energy is represented exactly for all basic motions.

The one-dimensional case
In order to illustrate the properties of the lumping technique presented above, let us determine the critical time step for a one-dimensional element containing a discontinuity. The purpose of this discussion is to study the influence of the location of the discontinuity on the critical time step. Figure 1 presents the 1D element containing a discontinuity. The length of the element is l, and the distance between the discontinuity and the left node is s. Using the notations of Figure 1, the linear shape functions are First, in order to set a reference critical time step, let us determine the critical time step for a standard element (i.e. without discontinuity) of section S, length l, Young's modulus E and mass density . For this element, the consistent mass matrix M FE and the stiffness matrix K FE are The critical time step given by Equation (24) is Using the mass lumping technique for this element with standard shape functions alone, we obtain: and the corresponding stable time step is In the following discussion, t 0 c will denote t lumped c , the critical time step for a standard finite element of length l with a lumped mass matrix. Now, let us deal with the presence of a discontinuity in the displacement field by using the eXtended Finite Element Method to enhance the basis of shape functions. The discontinuity at location s is described with the generalized Heaviside function H centred at s: The approximate displacement is With this enhanced basis of shape functions, the consistent mass matrix and the stiffness matrix are Using the lumping technique presented in Section 3.1, the mass matrix becomes In both cases (lumped mass and consistent mass), the critical time step given by Equation (24) is a function of the location s of the discontinuity. The results are summarized in Table I. This table gives the critical time steps with consistent mass and lumped mass for three different problems: the first problem is a standard finite element of length l; the second is an extended finite element with a discontinuity centred at s; the third problem is the equivalent finite element problemFE with two elements (one of length s and the other of length l − s, as shown in Figure 1). A plot of the critical time steps for the first and third problems is also shown in     Table II for these problems are zero. Consequently, the stability of the central difference method cannot be guaranteed for an arbitrary location of the discontinuity.
Using the proposed lumping technique, the critical time step also decreases with the distance between the discontinuity and the closest node, but its minimum value is not zero. Indeed, one can observe in Figure 4 and in Table II that the limit of t c when s goes to 0 or l is 1 √ 2 t 0 c . This result enables one to perform explicit transient analyses with discontinuous enrichment using a time step of reasonable size regardless of the location of the discontinuity.

The two-dimensional case
In this section, we consider 2D elements cut by a crack. The kinematics can be written as where f is the signed distance to the crack's surface. We perform the same development as in the one-dimensional case. The interesting point is the evolution of the critical time step as a function of the crack's location for an extended finite element with lumped mass matrix. Figures 5 and 6 show results for three different crack orientations (0, 30 and 45 • ) for a four-and three-node element, respectively. t c is the critical time step of the extended cracked element using the proposed lumping technique. In Figures  5 and 6, the results are normalized by t 0 c , the critical time step of the same element without extended functions. The critical time step is plotted as a function of S/S 0 , where S 0 is the area of the element and S the smaller area obtained when the element is cut by the crack. Similarly, for the three-dimensional case, Figure 7 shows the critical time step as a function of V/V 0 , where V 0 is the volume of the element and V the smaller volume obtained when the element is cut by the crack.  The conclusions are the same as in the one-dimensional case. The proposed lumping technique leads to a reasonable, nonzero minimum value of the critical time step. In practice, according to Figures 5 and 6, explicit transient analysis in the X-FEM framework with discontinuous enrichment can be carried out with a time step of 2 3 t 0 c .

Implementation of the technique and choice of the discontinuous function
The use of the generalized Heaviside function as the enhanced function is a particular case. Indeed, this function has the property that H 2 is a constant function equal to 1. Therefore, in this case, the following two lumping techniques give the same results as that proposed above (provided the coupling terms of the consistent mass matrix are not taken into account): The implementation of the proposed technique using the H function reduces to the use of one of the two standard lumping techniques already mentioned on the pure enriched terms of the consistent mass matrix. This would not be the case for enriched functions whose square is not constant. We show in Appendix A that for other discontinuous enrichments functions, the key factor in getting efficient time steps for transient X-FEM calculations is the choice of the lumping technique.

APPLICATION TO DYNAMIC CRACK PROPAGATION
In this section, we illustrate the effectiveness of the proposed method. Here, we focus on dynamic crack propagation in the framework of linear elastic fracture mechanics. First, we explain the main principles of the theory. Then, we develop the fracture criteria and present three numerical examples. Through these examples, we aim to validate the lumping technique and to show that the enrichment strategies X cla and X dis give similar results. In order to do that, each numerical example will be solved using: 1. implicit time integration with enrichment strategy, X cla ; 2. implicit time integration with enrichment strategy, X dis ; 3. explicit time integration with enrichment strategy, X dis .
As the crack propagates, the enrichments must evolve to take this propagation into account. The procedure we follow, which was described in Reference [10], consists in keeping the previously enriched degrees of freedom, adding new enriched functions to model the crack's extension and initializing the new degrees of freedom to zero. This procedure is energy-consistent and preserves the stability properties of the time integrator.

Linear elastic fracture mechanics
In linear elastic fracture mechanics, the concept of dynamic stress intensity factor [16] can be used to deal with crack propagation in the static case. Dynamic stress intensity factors are (see Reference [16] for details and reference therein): Since K dyn i are local quantities defined as limit values, they cannot be estimated directly and accurately. A possible approach was described in Reference [10] (with more details in PhD Thesis [17], in French): a two-field problem (consisting of the actual field u and an auxiliary field u aux ) is constructed for mixed-mode separation, and a domain-independent integral I int is obtained from the Lagrangian conservation law using a virtual crack extension q.
[(div( aux ).∇u(q) + ü.∇u aux (q)) + ( u aux .∇u(q) + u.∇u aux (q)] dS (45) Figure 8 shows the virtual crack extension field defined on two surfaces S1 and S2. q is particularized as follows: q=0 outside surface S2 q =1 on surface S1 tangent to the crack's faces everywhere (46) The use of a virtual crack extension field as defined above (where the size of S1 is chosen to be half the size of S2) decreases the influence of the numerical errors in the estimation of the mechanical fields near the crack's tip on the stress intensity factors. If K aux i denotes the stress intensity factors of the auxiliary field, we obtain from the dynamic equivalent of Irwin's relation: where f i are universal functions of the speed of the crack's tipȧ: In these definitions, c 1 and c 2 are the dilatational and shear wave velocities. The positive root of D(ȧ) = 0 defines the Rayleigh wave speed c r , which is the theoretical maximum speed for a crack in a homogeneous medium. Equation (47) shows that the stress intensity factors K dyn i can be estimated through an appropriate choice of u aux (K aux 1 = 1, K aux 2 = 0 for the determination of K

Fracture criteria
Considering an arbitrary 2D crack under mixed-mode loading, one must check whether the crack is going to propagate and, if necessary, determine the crack's speed and direction. In the following examples, we will consider that the fracture phenomenon is governed by the intensity of the hoop stress ( in the co-ordinate system associated with the crack's tip). Consequently, the crack's direction c is the direction of the maximum hoop stress, given by Then, the intensity of the loading at the crack's tip is calculated using an equivalent stress intensity factor K dyn : Once the crack's direction and the intensity of the stress field in that direction are known, the speed of the crack's tipȧ is given by where K 1c is the quasi-static fracture toughness and K 1D (ȧ) the dynamic fracture toughness, which depends onȧ through the following equation (see Reference [18]): Here, we will consider that m = 1 and the speed of the crack's tip is obtained from the equation:

Mode 1
The first example is an infinite plate with a semi-infinite crack. A theoretical solution of this problem for a given crack velocity is known (see Reference [16]). Since this analytical solution was obtained under these assumptions (i.e. infinite plate with a semi-infinite crack and a given speed of the crack's tip) and a numerical one for the geometry described in Figure 9, those could be compared for time t 3t c = 3h/c d (where c d is the dilatational wave speed). Beyond that, the reflected stress wave reaches the crack's tip and the analytical solution is no longer valid. As the wave reaches the crack, the mode-1 stress intensity factor for the moving crack can be written as For the moving crack, one can write where k is a universal function of the velocity of the crack's tipȧ, which can be approximated by Explicit X dis Implicit X dis Implicit X cla Finally, one has: The following numerical results were obtained with the plate dimensions h = 2m, L = 10 m, l = 5 m and the material properties E = 210 GPa, = 0.3, = 8000 kg m −3 . The tensile stress 0 was 500 MPa. The stress intensity factors were computed over a J -domain of width 0.1 m. The numerical results are normalized by 0 √ h. Since the first enrichment strategy X dis gives less accurate fields in the vicinity of the crack's tip, the mesh used with that strategy was composed of 60 × 120 linear quadrangular elements. With the strategy X cla , the mesh was 40 × 80 linear quadrangular elements. We will focus on two cases: a stationary crack and a crack propagating at prescribed constant speed v 0 = 1500 m s −1 after time t = 1.5t c .

Stationary crack.
For this case, using the implicit mean acceleration method ( = 0.25 and = 0.5), the time T was computed in 100 time steps t imp . With the central difference method, the size of the time step was t exp = t imp /2. Figure 10 shows the comparison of the results obtained forK 1 . Using the explicit scheme and X dis , the results are less accurate than with the implicit scheme and X cla . The plot of the results obtained with the implicit scheme and X dis enables one to verify that the oscillations in the Explicit-X dis solution are not caused by the lumping of the mass matrix but by the loss of accuracy in the mechanical fields due to X dis .

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

Kalthoff's experiments
This example deals with the numerical simulation of Kalthoff's experiments of the failure mode transition under pure mode-2 loading [22]. A schematic description of the problem is shown in Figure 12. Kalthoff's experiments are modelled in plane strain. A plate with two symmetrical edge cracks is impacted by a projectile at speed V 0 . The two cracks are centred with respect to the specimen's geometry and their distance corresponds to the diameter of the projectile. At low speed, Kalthoff observed brittle fracture. The two cracks propagate simultaneously with an overall angle from 60 to 70 • . If one increases the projectile's speed, a transition between Crack path Explicit X dis Implicit X dis Implicit X cla Explicit X dis Implicit X dis Implicit X cla Figure 14. Evolution of the crack's length. l = 0.05 m, and the initial crack's length was a = 0.05 m. The mesh was regular and consisted of 80 × 80 linear quadrangular elements. For the implicit time integrator, the time step was t imp = 1 s. For the explicit scheme, following the results obtained in the previous section, the time step was t exp = 2 3 t 0 c = 0.125 s. Figure 13 shows the crack's paths obtained with three different simulations. These results are very similar and the overall angle agrees with the experimental results. Looking at the details of the crack's paths, we observe the same phenomena as in Reference [12]: the crack starts to propagate with an angle of 65 • at t ≈ 28 s; then, a small deviation is observed at t ≈ 50 s; finally, the crack continues at the initial angle of 65 • from t ≈ 65 s. These changes in the crack's direction could be related to the propagation of the initial compressive stress wave within the specimen. The initiation occurs at t ≈ 28 s, after the wave has travelled from the left side to the right side, has been reflected against the free surface as a tensile stress wave, and finally reaches the crack's tip. Then, this tensile wave travels from the crack's tip to the left side of the specimen and reflects again as a compressive wave which arrives at the moving crack's tip at t ≈ 50 s. This compressive wave induces the crack's deviation, but is subsequently reflected against the right side as a tensile wave which reaches the crack's tip at t ≈ 65 s, after which the crack continues with the initial angle.
These results enable us to confirm the effectiveness of the method. We can also verify that the crack's evolutions for the three types of simulations carried out are identical by looking at the plot of the evolution of the crack's extension in Figure 14. Although, in this example, the same mesh refinement was used with both enrichment strategies X dis and X cla , the results are, again, very similar. Figure 15 shows the norm of the displacement field plotted on the deformed mesh for the three types of simulations.

Compact compression specimen (CCS)
The interest of the CCS experiment in our case is that there is always a mixed-mode loading at the crack's tip. The CCS problem is described schematically in Figure 16. The material's properties are those of PMMA: E = 5.76 GPa, = 0.42, = 1180 kg m −3 . The force P 1 (t) is  due to an impact at velocity V 0 = 20 m s −1 . The CCS is assumed to be linear elastic. Although the CCS is symmetric, the deformation and, therefore, the crack's path are not. This is due to the non-symmetric loading and boundary conditions. To validate the lumping technique and to show that the enrichment strategies X cla and X dis give quite similar results, each numerical example is again analysed using Table III: 1. implicit time integration and enrichment strategy, X cla ; 2. implicit time integration and enrichment strategy, X dis ; 3. explicit time integration and enrichment strategy, X dis . Figure 17 shows the final deformed shape. Figure 18 shows the calculated crack's paths for the three numerical models. All these results agree with the experiments (see References [23,24]). We can observe that the results of the explicit simulation with a time step of 0.25 s are similar to those of the implicit simulation using the same mesh. Figure 19 shows the evolution of the crack's length as a function of time for the three cases. One can observe that the three methods give the same results. In this case, we showed that using a time step of 0.25 s enables us to perform a stable calculation of dynamic crack propagation using the eXtended Finite Element Method. Figure 17 shows the crack's path Figure 17. Deformed mesh.

Crack path
Explicit X dis Implicit X dis Implicit X cla  Explicit X dis Implicit X dis Implicit X cla Explicit X dis Implicit X dis Implicit X cla Figure 19. Evolution of the crack's length.
superimposed with the nodes: one can observe that although the crack sometimes runs very close to a node, the simulation remains stable.

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