Mass lumping strategies for X‐FEM explicit dynamics: Application to crack propagation

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


INTRODUCTION
The extended finite element method (X-FEM) allows one to introduce a crack within an existing mesh without the need to modify it. Discontinuous enrichments are introduced within the elements cut by the crack using the partition of unity technique [1]. According to the original papers on the subject, the enrichment is composed of a tip enrichment [2,3] and a Heaviside enrichment [4] 1 away from the crack tip. Initially developed in the 2D setting, the method was then extended to 3D in [5]. At about the same time, two level sets were introduced to conveniently store the crack as two finite element fields [6,7]. This level set representation of cracks was also found to be extremely convenient to handle crack growth [7][8][9] (without any remeshing).
The implementation of the X-FEM requires some enhancements of the regular FEM implementation: • A specific integration scheme on elements whose nodes are enriched must be used [4,10].
• The number of degrees of freedom per element depends on the location of the element with respect to the crack. The use of an objected-oriented language such as C++ is quite convenient in dealing with the variable number of element degrees of freedom (in space and time as the crack is growing). • Conditioning issues are raised when many layers of elements are enriched around the crack tip (a pre-conditioning technique for this event is provided in [10]).
When dynamic loading is applied, additional issues must be addressed for the X-FEM to work properly. For transient analysis, one usually distinguishes between implicit and explicit algorithms. The standard Newmark implicit approach is known to be unconditionally stable in the finite element context for stationary cracks. As cracks are growing and remeshing is used, stability may be obtained, provided that the energy release rate at the crack tip is properly computed [11].A nice thing about the X-FEM is that the mesh is held fixed as the crack grows. The number of enriched degrees of freedom is, however, growing during computation. The stability of X-FEM with implicit algorithms was studied in [12] and with explicit algorithms in [13,14]. The latter paper [14] introduces a special lumping technique leading to a critical time step smaller (factor 1/ √ 2) than the one in the FEM case. In this paper we shall introduce a new lumping technique allowing one to use the same critical time step as in FEM (considering only Heaviside enrichment).
The paper is organized as follows. In Section 2, the critical time step is theoretically studied with two different lumping techniques for simplex elements (1D, 2D and 3D). Section 3 deals with effects of lumping techniques on wave propagation through discontinuity, and Section 4 is dedicated to numerical experiments: a compact compression specimen with both different lumping techniques is studied using two different element types (triangular and quadrangular).

EXPLICIT DYNAMICS FOR X-FEM
To capture the discontinuity in the displacement field at the crack tip, the partition of unity property of the finite element shape functions [1] is used to enrich the spatial interpolation as follows: In this equation, u is the approximate displacement field, N the set of nodes used to discretize the domain, and N i the finite element shape function associated with node i.Ū i are classical degrees of freedom, whereasŨ j are additional degrees of freedom supported by nodes j in the set N cut associated with the enrichment function H , see Figure 1. Only one enriched function is used here to describe the discontinuity in the structure. We choose the jump enrichment function as unit and symmetric. The body is separated by the discontinuity into two parts: + and − . Then the chosen enriched function H is defined as Other jump functions could be used [15] but the following study is valid only with H .I ti sw e l l known that the critical time step of an explicit integration scheme for dynamics may be estimated by computing the minimum value of the critical time step for all elements in the mesh taken separately. Moreover, this estimate is an upper bound for the exact critical time step. This result was obtained in [16] as an extension of a theorem first given by Rayleigh.
The paragraph above remains valid if the approximation is obtained through the X-FEM. We can thus concentrate on a single cracked element.

One-dimensional element
2.1.1. The shape functions. First of all we start with a 1D cracked element shown in Figure 2. The length of the element is L, its density , its section S and its Young modulus E. The crack cuts the element into two sub-segments of sizes L and (1− )L, respectively. Using the Heaviside enrichment [4], four approximation functions are defined over the element ( Figure 2): two classical and two enriched. We assume that the 'Heaviside'-type enrichment function is +1 on the left side of the bar and −1 on the right side. Equivalently, the four functions shown in Figure 3 may be used. The latter base is in the spirit of Hansbo and Hansbo [17]. The equivalence of the basis depicted in Figures 2 and 3 was proven in [18]. Moreover, Song et al. [19] and Areais et al. [20] contain a proof of the equivalence between Hansbo and Hansbo [17] and X-FEM scheme and present numerous results for explicit/implicit dynamic X-FEM simulations. The two basis functions are related by Using the approximation functions depicted in Figure 2, an approximated field u on the 1D element reads Alternatively, one can use the equivalent basis depicted in Figure 3. The approximation reads with this basis as follows: Using link (3), the degrees of freedom of the approximations are related by It can be noted that the basis functions used in Figure 3 are as if the two sub-segments were completely independent. The approximation on the left segment is as if the right one was void and vice-versa. On each sub-segment the approximation is exactly the one used to model voids with non-matching meshes as described in [21]. It is thus tempting to use the lumping strategy developed in [22] for voids leading to a critical time step independent of and value h/c, where c is the wave velocity c = √ E/ , i.e. the same critical time step as for a regular element.

The mass and stiffness matrices.
Following [22], the lumped mass matrix in our case is expressed as The stiffness matrix is expressed as (see [14]): The mass and stiffness matrices determine the critical time step of the element.

The critical time step.
For an explicit scheme, the critical time step is computed as 2/ max , where max is linked to the maximum eigenvalue of the generalized system: Thus, we obtain in our case Hence, the biggest value for is when is 2, thus, max = (2/L) √ E/ . Then the critical time step is L √ /E, which exactly corresponds to the critical time step of the 1D finite element problem. The result obtained is promising: the critical time step of a 1D element whose nodes are enriched is exactly the same as that of a standard element. Moreover, this result does not depend on the shape function basis used: (1, 1 ′ , 2, 2 ′ ) or (I, I ′ , II, II ′ ). The proof is where matrix P is defined in Equation (6) such that det(P·P T ) = 16 and P·P T = P T ·P = 2I 4 .Note that if the modified enrichment proposed in [3] (which was shown to be equivalent to the Hansbo basis in [18]) were used and the mass matrix were lumped in this enriched basis following Equation (7), the result would be exactly what we obtain by lumping the mass matrix in the (1, 1 ′ , 2, 2 ′ ) basis. Hence, going back to the X-FEM basis functions, we apply the transformation to the above mass matrix We did obtain a block-diagonal mass matrix, the blocks being of size 2.
Finally, had we decided that node 1 was on the negative side of the 'Heaviside', the results would have been different. The results would have been the one obtained by replacing by 1− in (14). Result (14) is thus valid whatever the choice of sign for the Heaviside, provided that is always defined as the matter fraction on the positive Heaviside side.

Kinetic energy conservation.
This part aims at checking that the lumping techniques conserve the kinetic energy for rigid motions with a discontinuity.
In the classical finite element setting, the sum of the entries in a mass matrix is always equal to the total mass. The reason for this is that when all degrees of freedom are set to one, we obtain a rigid translation mode. In the X-FEM case, the sum of the entries of the matrix is not equal to the total mass because setting all degrees to 1 does not give a rigid mode. However, if we sum the entries corresponding to a rigid mode (first and third lines and column) we recover the total mass. Let us consider a 1D element with two nodes. Each node has ordinary degrees of freedom corresponding to the shape functions f I and f II and additional degrees of freedom corresponding to the enriched function H . The fraction ratio of this element is ∈[0, 1]. First we consider a rigid motion at velocity V (as if there was no crack). The exact kinetic energy is where m = SL is the mass of the element, S its section, L its length and its density. And the discretized velocity field corresponding to this motion is described by the vector (in the Hansbo basis (1, 1 ′ , 2, 2 ′ ) as shown in Figure 3): Hence, the discretized kinetic energy reads Let us now consider the rigid motion of a cut element whose parts move away at velocity V . The discretized velocity vector is now And the discretized kinetic energy can be expressed as The case corresponding to the use of the other lumped mass matrix has already been verified in [14].
Let us now consider a more complex structure. Figure 4 presents a bar composed of three elements; the middle one is cut. Each element has its length equal to L, its Young's modulus E, its section S and its density . The length of the bar is l = 3L; hence, its mass is m = Sl.Herewe check the conservation of kinetic energy for the whole bar. Figure 5 presents the X-FEM shape functions: four standard and two enriched functions. Figure 6 presents the six shape functions used for this 1D case: the Hansbo basis is chosen here. The discretized displacement vector is of size 6 and is expressed as Hence, the approximate displacement field U reads The diagonal mass matrix [14] for the whole structure, in the standard basis (N 1 , N 2 , HN 2 , N 3 , HN 3 , N 4 ), is as follows (with l = 3L): The block diagonal mass matrix [22] is  Figure 4. A three-element structure cut by a crack. Figure 5. Shape functions for the three elements mesh: N 1 , N 2 , HN 2 , N 3 , HN 3 and N 4 . M diag is the diagonal mass matrix obtained using the technique developed in [14],andM block-diag the block-diagonal mass matrix obtained in [22]. As in the last paragraph let us consider two rigid motions.
Motion 1: First the whole structure moves as a rigid body in the same direction at velocity V . Hence, the velocity vector is , and for both lumping techniques, one obtains Motion 2: For the second motion, the structure at the left of the crack moves at velocity −V , whereas the right part moves at velocity V . Thus, the two parts move away at velocity V . Hence, the velocity vector is expressed as for both lumping techniques: Hence, one has to note that the two lumped mass matrices allow to conserve the kinetic energy for both rigid body motions; moreover, the crack is taken into account in one motion. These verifications about kinetic energy conservation mean that both lumping techniques do not add mass.
To conclude with the 1D case, the lumping technique developed in [22] allows one to obtain the same critical time step for an element with additional degrees of freedom than for a standard element. This result improves the technique explained in [14] which was only √ 2s m a l l e rt h a n the standard one. However, the cost here to obtain a better result is that the mass matrix depends on the fraction ratio, whereas the mass matrix in [14] was simply diagonal constant. In addition, we show that these two lumping techniques allow one to preserve the kinetic energy for rigid motions.
Now the use of each technique can be motivated by different reasons: the use of block-diagonal matrix will be preferable, except for pure explicit codes which do not own a matrix (as LS DYNA, RADIOSS, EUROPLEXUS). For these codes, a constant diagonal mass matrix will be enough unless they were able to store 2×2 matrices for each enriched node.

Two-dimensional elements
This part studies the critical time step for triangular and quadrangular elements. Moving to 2D elements, we shall see that the mass matrix keeps the same topology as in (14). Figure 7: a triangular element cut by a crack and characterized by its fraction ratio . As in the 1D case, we use roman indices to denote the classical shape functions and the roman prime indices to denote the enriched shape function. The arabic indices denote the truncated shape functions. For instance, the function f 1 is equal to the classical shape function f I on the sub-triangle containing node 1 and is zero on the other sub-triangle. Finally, the prime arabic indices indicate the 'complementary' truncated functions. For instance,

Triangular element. Consider a triangular element depicted in
The truncated and enriched basis functions are related by  Figure 7. A triangular element of area S cut by a crack. The fraction ratio is denoted by .
where H (·) indicates the sign of the Heaviside at the corresponding node. The mass matrix in the truncated base reads Applying the proper transformations, based on Equations (26)-(31), we obtain the following mass matrix for the enriched basis: We have assumed that the material fraction is computed on the positive side of the Heaviside. Otherwise, must be replaced by 1− . The obtained mass matrix is basically the same as in the 1D case. For 3D elements, the 2×2 block-diagonal matrix will not change, except for the factor in front of the matrix, which will be the element volume divided by the number of nodes of the element. The results for a triangular element are presented in Figure 8. Figure 8 shows the critical time step normalized by the finite element critical time step obtained with a lumped mass matrix. The two results are for the two different lumping techniques, which use a diagonal mass matrix [14] and a block diagonal one [22], respectively. It underlines the fact that with the block-diagonal mass matrix the critical time step of the enriched triangular element is the same as that of the triangular finite element, whatever the position of the crack in the element. Hence, the result is similar to the 1D case. Moreover, Figure 8 shows that the critical time step for the diagonal case is quite smaller than for the block-diagonal case. Note also that the block-diagonal 2×2 matrix is not the same for all the enriched nodes. It depends on the crack path.
To conclude, using a block-diagonal mass matrix, one obtains the same critical time step as in the finite element case, hence, around twice the value of the critical time step obtained using constant diagonal mass matrix. The block-diagonal mass matrix is composed of 2×2 matrices, each of which depends on the fraction ratios of the elements around the considered nodes. The results shown in this part are quite interesting because the critical time step for an element with additional degrees of freedom is the same as for a standard one. It means that an enriched triangular elements and a standard triangular element have the same critical time step when the mass matrix is lumped by the appropriate technique.

Quadrangular element.
Let us now consider a quadrangular element. Figure 9 describes the fraction ratio on a quadrangular element cut by a crack. The critical time steps are shown in Figure 10   . but it has to be noted that the use of the block-diagonal mass matrix allows one to obtain a bigger critical time step. We must note that the critical time step of the standard element is never reached for any fraction ratio , whereas it was for the triangular case. Thus the biggest critical time step is reached for = 0.5 and its normalized value is 0.93.

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

Three-dimensional elements
Let us now move to the 3D case and the studies of critical time step of tetrahedral and cubic elements. Figure 13 presents a tetrahedral element cut by a crack and Figure    mass matrix [14] and a block-diagonal mass matrix [22], respectively. For this 3D element, the critical time step obtained with a block-diagonal mass matrix is exactly the same as that of the non-enriched finite element. Hence, the results in this paper are similar for all dimensions: 1D linear element, triangular element and tetrahedral element allow one to obtain an X-FEM critical time step equal to the FEM critical time step. Figure 15 presents a cubic element cut by a crack and Figure 16 shows the critical time step of this cut element as a function of the position of the crack in the element. The time step is normalized by the critical time step of the classical finite element. The results of the critical time step of the enhanced element are available for the two different lumping techniques. These techniques give different results as shown in Figure 16. The critical time step results are quite similar to those of the quadrangular element. None of the two techniques are able to obtain the same critical time step as the one of the finite element problem. Nevertheless, the use of a block-diagonal mass matrix allows a bigger time step.

Results analysis
In this section, we take a deeper look at the reasons as to why we can obtain the same critical time step as for the finite element case for some elements: triangular and tetrahedral elements with the block-diagonal mass matrix give the same critical time step as for the finite element. The fact that the critical time step of a standard element can be obtained with an enriched one is directly related to the shape functions of the element. Using the Hansbo basis for an element whose nodes are enriched, we have for a 1D element (as shown in Figure 3) Thus, the X-FEM stiffness matrix is proportional to the FEM stiffness matrix on each side of the discontinuity where K FEM is given in the Appendix. The same may be said for the lumped mass matrix given in Equation (7). Thus, both stiffness and mass matrices are proportional to the FEM on each side of the discontinuity. Then, the resolution of Equation (9) is similar to the resolution in the classical finite element case, which explains why the X-FEM critical time step is the same as the FEM one. For simplex elements (1D bar, triangular and tetrahedral), the standard shape functions are linear in space and are written as the following functions k (where k ∈{1 ...n nodes }, n nodes is the number of nodes in the element): k (x, y, z) = I 0 + I 1 .x + I 2 .y + I 3 .z (36) Hence, the gradients of the shape functions are constant as explained in Equation (34). This is the reason as to why the X-FEM stiffness matrix is proportional to the finite element one. Note that this is not true for the quadrangular and cubic elements because the gradients of shape functions (used to compute the stiffness matrix) are not constant in space; thus, the X-FEM stiffness matrices are not proportional to the FEM matrices. The shape functions for these (quadrangular then cubic) elements are not linear and may represent the following polynomials: These results are summarized in Table I, which presents the normalized critical time step for the five studied elements: 1D linear, triangular, tetrahedral, quadrangular and cubic elements, with the two lumping techniques (diagonal and block diagonal).

WAVE TRANSFER THROUGH THE DISCONTINUITY
First the mass matrix is used in the explicit algorithm as We now check that the dynamic wave propagation is not transmitted through the discontinuity. This has to be checked for the two different lumping techniques. density = 1190 kg/m 3 , length L =l/6 = 10 mm and section S = 2mm 2 . m = S6L = Sl is the mass of the whole structure. The critical time step is given by t c = L E ;h e r e t c = 6.06 s. The time step for calculation is taken as 1 s, and the total time is 50 s. The fraction ratio is 0.2 in the third element. The loading is an initial velocity V 0 = 10 m/s at node 7. The structure is assumed to be linear elastic. Figure 18 presents the standard shape functions of the 6-element mesh: there are seven standard shape functions (one per node) and two enriched ones (for the cut element nodes). The right extremity of the structure is loaded with an initial velocity. Afterwards the tension wave propagates and reaches the discontinuity. We would like to check whether or not the structure on the left part of the discontinuity is affected by the wave near the discontinuity.

One-dimensional elementary case
In the explicit algorithm, the important equation to know whether the wave will propagate through the discontinuity is as follows:Ü And in our case, there is no external force, and the material model is elastic linear; hence, this term is just written asÜ The fact that wave goes through the discontinuity would be due to a non-zero right-hand side for node 3 in Equation (43). The results of wave propagation are summed up in Figure 19, which presents five graphs for the two different mass matrices: diagonal and block diagonal. The five graphs in Figure 19 show the displacement, velocity, acceleration, stress and strain fields along the structure at time 25 s. Note that the two different lumping techniques give similar results. At the left side of the discontinuity At the right side of the discontinuity (H =+1)( H =−1) This is due to the mass matrix, which is not the same in the two cases. Figure 19 shows that there is no wave transfer through the discontinuity. Indeed nothing happens on the left side of the discontinuity. For this simple 1D example, let us now prove analytically that no wave transfer occurs with the two lumping techniques.

Proof by mathematical induction
On the basis of Figure 17, the displacement discretization is considered as follows: and the proposition to prove is chosen as shown in Table II. One notes that the displacement of node 3 (respectively, node 4) is u 3 +a 3 (respectively, u 4 −a 4 ) because of the choice of the enrichment function H ,whichis+1 on the left side of the discontinuity and −1 on the right side.
To initialize the proof, we consider that the statement described in Table II is verified at time 0. Now we suppose that it is also verified at time n, and we prove that it is also true at time n +1. First, to pass from time n to n +1, Equation (39) gives the displacement at time n +1. We easily obtain that the displacement field [U n+1 ] preserves the properties of Table II where [U n+1 ]=[u 1 , u 2 , u 3 , a 3 , u 4 , a 4 , u 5 , u 6 , u 7 ]; the stiffness matrix is expressed as follows: and the lumped mass matrix for both lumping techniques is expressed as with being a parameter selecting the lumping technique: = 1 corresponds to the block-diagonal mass matrix and = 0 to the diagonal mass matrix. Afterwards, we compute Equation (45) and subsequently obtain the following relationships for the acceleration field: Thus, the acceleration properties of Table II remain valid at time n +1. Third, Equation 41 allows one to verify that the velocity field at time n +1 preserves the properties of Table II. Finally, the initial properties described in Table II are conserved along time. The fact that zero values on the left side of the discontinuity will not change for both lumping techniques means that tensile waves do not propagate through the discontinuity.

Two-dimensional case
Let us now consider the plane strain wave propagation block described in Figure 20. The dimensions are L = 5 mm, W = 10 mm and materials properties are those of PMMA: Young's modulus E = 3.24 GPa, the Poisson ratio = 0.35 and its density = 1190 kg/m 3 . The loads are just impact velocity in the tension direction (see Reference [15]). This impact velocity is defined as where t r = 0.1 s and the velocity V 0 = 10 m/s. The discretization of the structure is 78×39, with quadrangular elements. The time step was 0.1 s, whereas Remmers et al. [15] used a time step 100 times smaller. Let us now simulate this test with both lumping techniques: diagonal and block-diagonal matrices. Figure 21 describes the profile of the stress yy at the centre of the plane strain wave specimen at t = 1.5a n d2 s with both lumping techniques. The reason as to why one sees only one curve in Figure 21 is the fact that results are identical for both techniques. One notes that the tensile stress average is about 25 MPa and propagates from the top side of the specimen to reach the crack at around t = 2 s, which agrees with [15]. Figure 21 shows that the stress wave at t = 3 s is properly reflected, and the lower part of the specimen remains stress free for both lumping techniques. Hence, the no-wave propagation through the discontinuity using both lumping techniques is acceptable.

NUMERICAL APPLICATIONS
The compact compression specimen (CCS) problem is described schematically in Figure 22.  P 1 (t) is due to an impact at velocity V 0 = 20 m/s. The CCS is assumed to be linear elastic. Although the CCS geometry is symmetrical, the deformation and, therefore, the crack path are not. This is due to the non-symmetric loading and boundary conditions. We carry the computations with the two methods: first using a block-diagonal mass matrix and a time step close to the corresponding critical time step of this method, and second the technique using a diagonal mass matrix with its corresponding critical time step. Both the results agree with the experiments [23,24]: crack path and velocity of the crack tip. Moreover, the CCS experiment was simulated with two different meshes as shown in Figure 23 (using triangular or quadrangular elements) and with both lumping techniques for the mass matrix: these are the four simulations. The experiment time was 120 s and the critical time steps are 0.63 and 0.85 s for triangular and quadrangular classical elements. One notes that the different mesh strategies (quadrangular or triangular) give exactly the same crack paths for each lumping Crack path diagonal quadrangular diagonal triangular Reference (Implicit) Figure 24. Crack paths for both lumping techniques. method; hence, Figure 24 represents only two paths corresponding to the two lumping techniques. Furthermore, the results obtained with the diagonal mass matrix was compared in [14] with those obtained using an implicit time integration that stands for a reference. Hence, Figures 24 and 25 present crack paths and lengths for both techniques compared with the implicit computation results. The critical time steps corresponding to each lumping technique are different as given in Table I: the one corresponding to the block-diagonal lumping technique is bigger. However, the block-diagonal approach will require more effort at each time step since a 2×2 matrix based on the fraction ratios of the cut elements needs to be built (and inverted). Note that if one uses the truncated basis directly, a truly diagonal mass matrix is obtained (but it still depends on the fraction ratio). This additional effort at each time step may be less of an issue for a non-linear problem for which the internal forces, computation is where most of the time is spent. However, both ways exposed to treat X-FEM in an explicit dynamic code can be interesting. The CPU time of both may be around the same if the integration of block-diagonal mass matrix is associated with effective strategies for the inversion of block 2×2 and evaluation of fraction ratio. Indeed fraction ratio could be computed only once for the cut elements, and the inversion of block could also be performed once, because the quantities of a cut element (fraction ratio and mass matrix) will not change after the propagation of the crack. In fact, the computation of fraction ratio will be performed only for the new cut elements, and the inversion of the block can be performed analytically to be integrated easily in an explicit code. Thus, X-FEM with the diagonal mass matrix technique has been integrated in EUROPLEXUS.

CONCLUSION
In this paper, we introduced a new lumping technique for the mass matrices of meshes enriched with Heaviside functions within the X-FEM framework. This lumping technique yields the same critical time step at the element level as the one for standard elements. The lumped mass matrix is not strictly diagonal but rather block diagonal. A 2×2 matrix needs to be stored at each enriched node. This additional storage provides a better critical time step than the one obtained using a true diagonal mass matrix for the following elements: 1D bar, triangular, quadrangular, tetrahedral and cubic. Numerical experiments demonstrate the robustness and stability of the approach. Then the two lumping techniques (the one developed in this paper and the other in [14])w e r e compared with different meshes in the 2D case. It aims at showing that the lumped mass matrix does verify some simple validations: first, studies on critical time step are checked by the simulations, and second, numerical simulations prove that no wave transfer through the displacement discontinuity appears with both lumping techniques. Furthermore, they both preserve kinetic energy for different rigid motions.
In order to choose between the two lumping techniques, the following points must be taken into account: the use of block-diagonal mass matrix is slightly more complex to implement than the true diagonal mass matrix; the diagonal mass matrix does not depend on the element fraction ratios, whereas the block-diagonal matrix does. On the other hand, the bigger time step offered by the block-diagonal approach may be interesting, in particular, for non-linear applications for which the constitutive law integration is most likely predominant compared with mass matrix inversion.