An energy‐conserving scheme for dynamic crack growth using the eXtended finite element method

This paper proposes a generalization of the eXtended finite element method (X‐FEM) to model dynamic fracture and time‐dependent problems from a more general point of view, and gives a proof of the stability of the numerical scheme in the linear case. First, we study the stability conditions of Newmark‐type schemes for problems with evolving discretizations. We prove that the proposed enrichment strategy satisfies these conditions and also ensures energy conservation. Using this approach, as the crack propagates, the enrichment can evolve with no occurrence of instability or uncontrolled energy transfer. Then, we present a technique based on Lagrangian conservation for the estimation of dynamic stress intensity factors for arbitrary 2D cracks. The results presented for several applications are accurate for stationary or moving cracks. Copyright © 2005 John Wiley & Sons, Ltd.


INTRODUCTION
The modelling of crack growth is a problem of great importance in the simulation of failure. Many techniques have been developed to take into account the discontinuity produced by a crack. First, Benzley [1] presented an enrichment to finite elements using asymptotic solutions of static fracture problems. This idea was also developed by Fleming et al. [2] for the element-free Galerkin (EFG) method. The partition of unity method (PUM) introduced by Babuska and Melenk [3] gives a theoretical framework to new techniques such as hp-clouds, generalized finite elements (GEM) or eXtended finite elements (X-FEM). For static fracture problems, Black and Belytschko [4], Moës et al. [5] and Dolbow et al. [6] developed the 1 X-FEM method. This technique was used in combination with the level set method by Gravouil et al. [7] for arbitrary 3D cracks. This method was also generalized to arbitrary discontinuities by Belytschko et al. [8]. It was also applied, then modified, in References [9,10]. In the field of dynamic fracture, EFG and GEM have already been used to simulate the propagation of arbitrary 3D cracks (see References [11,12]). Chen et al. [13] and, more recently, Belytschko et al. [14] developed a new enrichment technique to avoid the difficulties encountered with the original X-FEM in time-dependent problems [4][5][6]. In Reference [14], instead of a singular enrichment of the original X-FEM, the authors use a discontinuous enrichment combined with a cohesive law and damage model. The present paper proposes a methodology for using a singular enrichment for time-dependent problems.
A study of energy conserving methods for classical FEM with remeshing was introduced in Reference [15]. The present paper proposes a generalization of the X-FEM to model dynamic fracture and time-dependent problems in a general sense, and gives a proof of the stability of the numerical scheme in the linear case. This stability study is an extension of the balance recovery method presented in Reference [15] for a specific implementation of the X-FEM. Our approach enables us to simulate the dynamic propagation of arbitrary 2D cracks using an enrichment strategy for time-dependent problems. This strategy is justified by a theoretical study of Newmark-type time integrators for problems with evolving discretizations. This study yields stability conditions which also ensure energy conservation. Using this approach, as the crack propagates, the enrichment can evolve with no occurrence of instability or uncontrolled energy transfer.
A technique for the estimation of dynamic stress intensity factors for arbitrary 2D cracks, based on the law of conservation of the Lagrangian, is also presented. This technique can be viewed as an extension of the work by Attigui and Petit [16]. The method uses the concept of virtual crack extension to obtain a domain-independent integral and auxiliary fields for mixed-mode separation. Previous works can be found in References [7,11,[17][18][19][20][21][22][23][24]. The development presented below takes into account the capabilities of the X-FEM whereby the implicit description of the crack's geometry enables one to construct a virtual crack extension field tangent to the crack's faces everywhere. This specificity makes it possible to write a domain-independent integral for an arbitrary crack.
In Section 2, the reference problem is presented and discretized using the X-FEM. In Section 3, the stability of this scheme is studied in the linear case and the discretized energy balance is written for dynamic problems with evolving discretizations. The enrichment technique is also presented in this section. The approach based on the Lagrangian conservation law is developed in Section 4 in order to calculate dynamic energy release rates and stress intensity factors using a domain-independent integral. In the last section, several examples of calculations of stationary or moving cracks are presented to show the accuracy of the method.

The continuous problem
Here, we present a general definition of the dynamic crack growth problem. For this type of calculation, we add an additional unknown a(t), representing a measure of the crack (its length in two dimensions), to the usual displacement and stress fields. One can express the reference problem as follows: Reference problem: Given T ] such that: where (t) is the domain being considered (Figure 1), * 1 the boundary along which displacements u d (t) are prescribed and * 2 the boundary along which forces F d (t) are prescribed. f d (t) are the prescribed volume forces; u,u andü the displacement, velocity and acceleration fields; and the symmetric stress and strain tensors; G and K i the dynamic energy release rate and stress intensity factor; C the Hooke's tensor; the mass density; G c the critical energy release rate;ȧ and (t) the velocity of the crack's tip and its orientation in the cylindrical co-ordinate system of the crack's tip; U and S the function spaces associated with the problem, and U 0 the vector space of the virtual fields defined by In numerical simulations, the first step consists in calculating the dynamic energy release rate and the dynamic stress intensity factors; the second step consists in taking the crack's extension and the new geometry into account.

Remarks
(i) Here, a weak form of the problem is presented. It is interesting to review its underlying local equations in order to pinpoint the boundary conditions being considered: (ii) Although the domain was defined as time-dependent ( (t)), we consider that the crack's length alone (a(t)) changes. We make the simplifying hypothesis that the initial and deformed configurations are the same and we assume small perturbations. (iii) We assume the material to be homogeneous with linear isotropic behaviour. Therefore, the problem falls within the framework of linear dynamic fracture mechanics, in which one can define the dynamic energy release rate G and the dynamic stress intensity factors K i (see References [18,22,25]).

The discrete problem
In our approach, we use the X-FEM first introduced in References [4,5]. In this method, an enrichment to the classical finite element approximation is defined using the partition of unity method developed in Reference [3]. For a static problem (see References [5,26]), the displacement field can be written using an enriched basis of shape functions: where N is the set of all the nodes in the mesh, N cut the set of the nodes which belong to elements completely cut by the crack and N branch the set of the nodes which belong to elements partially cut by the crack. N i are the classical shape functions and H is a discontinuous function whose value is 1 if x is above the crack's surface, −1 otherwise.
[B ] are branch functions: In our approach, all the fields (displacement, velocity and acceleration) are discretized using Equation (7). Consequently, the discrete enriched problem can be written as a classical dynamic problem formulated according to the finite element method: we define the mass matrix and the stiffness matrix in the usual manner. The discrete momentum balance at time t n+1 is where M n+1 and K n+1 are the mass and stiffness matrices at time t n+1 . Since we intend to represent a moving crack, we have to solve a discretized problem whose number of degrees of freedom (size of the state vector) and crack geometry are evolving. Let us consider a firstorder tensor V . Its value at time t i can be written using the set of shape functions of the For the second-order tensors, let us designate the mass and stiffness matrices by We will now look at what happens during a time step t from time t n to time t n+1 .A s shown in Figure 2, one has to perform two basic steps: change the discretization and solve the momentum balance equation. Changing the discretization (step a) requires the determination of X n+1 n (the state vector at time t n written on the set of shape functions at time t n+1 ), which will be detailed in Section 3. In step b, the momentum balance equation is solved from time t n to time t n+1 in order to obtain X n+1 n+1 . This equation is discretized using a Newmark-type scheme with constant and (Equations (13) and (14)). Working with two different discretizations also requires that we rewrite the Newmark scheme: Now, we can define the discretized reference problem: Discretized reference problem: Given U n n ,U n n ,Ü n n a n ,ȧ n , find such that: 5 The quantities U n+1 n ,U n+1 n andÜ n+1 n are additional unknowns of the problem at each time step. We will detail their calculation in Section 3.

Remark
The discretized problem presented here leads to five equations with eight unknowns, which means that its resolution requires three projections.

STABILITY ANALYSIS AND DISCRETIZED ENERGY BALANCE FOR DYNAMIC CRACK GROWTH SIMULATION
In this section, we will extend the method presented in Reference [15] to the case where the X-FEM is used. This method enables us to study the stability and the discretized energy balance of dynamic analyses with evolving meshes. When one uses the X-FEM for dynamic crack growth analysis, the mesh does not change but the basis of shape functions varies as the crack grows. Consequently, we can use the same approach and the same notations. First, we develop the stability analysis using the energy method (see References [27,28]). Then, we study its application to dynamic fracture mechanics and present a method to improve the stability and accuracy for an evolving discretization. Finally, we present the application of this method, called the balance recovery method, in the context of the X-FEM.

Stability study
Let us now consider the case of a calculation with an evolving discretization. One can study the stability of this scheme by extending the method presented in References [27,28] to this type of calculation. Let us use the notations defined in Section 2 and define the following operators: Working with two different discretizations requires that we rewrite the equilibrium equations at times t n and t n+1 on the same discretization. Since the scheme must be stable in the absence of external loads, these are not considered in the stability study (see Reference [29]): Using the approach developed in Reference [27] and observing that: we obtain the stability equation: On the right-hand side, we recognize the term in − 1 2 which governs the stability of the scheme for a calculation with constant discretization. However, the presence of a supplementary term involving [M] and [K] destroys the stability conditions of the standard case. The stability depends on the right-hand side of Equation (20).
The method presented here enables one to perform the stability study when going from X n+1 n to X n+1 n+1 (step b). In order to treat the problem for the whole time step, one must study the transition X n n , X n+1 n+1 (step a followed by step b). For this purpose, one writes:

The discretized energy balance
With the notations above, the energy balance can be written in the same way as the stability study. We get: Again, these are the equations of a problem whose discretization does not evolve, with an additional term involving [M] and [K].
Equation (22) enables one to express an energy balance between X n+1 n and X n+1 n+1 (step a). Finally, one must calculate the balance between X n n and X n+1 n+1 (step a followed by step b), which is W ext is the work of the external loads, which can be written as follows:

Application to fracture mechanics: the balance recovery method
Let us now consider the case of dynamic crack growth analysis. In this particular case, we know what physical phenomenon causes the change of geometry. Indeed, let us assume a crack growth a between times t n and t n+1 . The state vector at time t n can be in equilibrium on this new geometry only if a distribution of forces F + is applied along the crack's extension to close it. By making the residual defined by (24) vanish, the balance recovery method guarantees that the state vector projected onto the new discretization is in equilibrium (assuming that there are no external forces): Thus, we have the following relations: correspond to the power and the work of the force distribution F n+1 + , respectively. We also know that these two terms are, from a physical standpoint, equal to −2G a and −2Gȧ, respectively (G being the energy release rate andȧ the crack's velocity, see Reference [18]). This enables us to verify the relations: Let us return to the stability study and energy balance described above. Concerning the energy balance, we will focus on the influence of the balance recovery on the two stages performed during the time step t n ,t n+1 . The term X n+1 n+1 ,X n n denotes the complete energy balance defined as the difference between the left-hand side and the right-hand side of Equation (23) taking Equation (27) into account.
X n+1 n+1 ,X n+1 n is defined as the difference between the left-hand side and the right-hand side of Equation (22) taking Equation (27) One can immediately deduce that X n+1 n+1 ,X n n = X n+1 n+1 ,X n+1 n + X n+1 n ,X n n (30) One can observe (Equation (30)) that X n+1 n+1 ,X n n is the sum of X n+1 n ,X n n , which is the quantity of energy introduced or dissipated during the projection operations (step a), and X n+1 n+1 ,X n+1 n , which is the energy introduced or dissipated when one goes from X n+1 n to X n+1 n+1 (step b).

Remarks
(i) The balance takes into account the energy required to create new crack surfaces (2G a). This term appears naturally because of its relation to the work of the force distribution F n+1 + (see Reference [18]). (ii) This method enables us to set X n+1 n+1 ,X n+1 n equal to zero a priori and, therefore, to minimize the amount of energy introduced by the successive discretization change operations.
(iii) Using the classical FEM, the discretization change step (remeshing) cannot be controlled and can introduce numerical energy into the simulation.

Use of the balance recovery method in the X-FEM
In the method presented above, we considered dynamic calculations with evolving discretizations. When one uses the X-FEM to simulate dynamic crack propagation (as presented in Reference [15]), the mesh does not change between time t n and t n+1 , but new space shape functions are added to simulate the crack's growth. Consequently, the mass and stiffness matrices at times t n and t n+1 are different and one can apply the balance recovery method. If we choose to retain all enrichments from time t n to time t n+1 (see Figure 3), the basis of shape functions at time t n+1 is greater than that at time t n and one can choose the new degrees of freedom to be initially zero. The projected displacement field is The new stiffness matrix is: The force distribution F + needed to close the crack's extension is replaced by setting new degrees of freedom which represent the crack's extension to zero. Consequently, the state vector calculated at time t n with the discretization at time t n+1 verifies the momentum balance on the geometry at time t n+1 (see Equation (24)). Moreover, the total energy associated with this state vector is the same for both discretizations because the initialization reduces the contribution of the new shape functions to zero The same proof holds for the kinetic energy term. The equilibrium recovery method for the X-FEM is easy and consists simply in initializing properly the new degrees of freedom which appear as the crack grows. This method guarantees that the numerical scheme is stable and that no numerical energy is introduced when the crack propagates.

Enrichment strategy
The previous section showed how to introduce new degrees of freedom when using the X-FEM. We will now develop the enrichment strategy. As shown in Figure 3, a new singular enrichment is added onto the new set N n+1 branch . These new shape functions are associated with the new front, i.e. the local cylindrical co-ordinate of the crack's tip for these functions is centered on the new crack's tip. The old singular enrichments are kept and always associated with the old tip. For a discontinuous enrichment, new shape functions are added only on the set N new cut = N n+1 cut \N n cut . If the crack propagates but remains within an element, singular shape functions are added onto the same set of nodes N n+1 branch = N n branch , but there is no new discontinuous shape function. With such an enrichment strategy, the number of degrees of freedom increases as the crack propagates and the balance recovery method can be used.
As shown in Reference [30], for the partition of unity method with an explicit Newmark-type scheme, the stable time step of the enriched problem is a small fraction of the stable time step of the problem which contains no enriched shape function. For this reason, we now choose the mean acceleration scheme, which is unconditionally stable.

Interaction integral
In dynamic fracture mechanics, the dynamic energy release rate of a two-dimensional problem can be deduced from the variation of the Lagrangian corresponding to a variation of the crack's length. In this section, we merely outline the method (see Appendix A for full details). We define the Lagrangian L by the following expression, where W is the strain energy, T the kinetic energy and W ext the external work: For mixed-mode separation, we will consider a problem with two fields u and u aux . u is the actual displacement field and u aux an auxiliary displacement field, both of which are kinetically admissible, and we can decompose the Lagrangian of the total displacement field u tot = u + u aux as In this equation, L int u,u aux is the interaction Lagrangian, whose expression can be developed into: where W int = :∇u aux = aux :∇u = u k,k ij + u i, j + u j,i u aux Let us consider the contour and domain represented in Figure 4. Assuming that the Lagrangian conservation law holds, that the virtual crack's extension l is collinear to the tangent to the crack's face and that the crack's faces are traction-free in both the actual and the auxiliary solutions, one obtains the following path-independent integral (n j being the components of the outward normal vector to the closed domain A( ) defined by its contour enclosing the crack's tip): where In the case of a straight crack,ȧ denotes the velocity of the crack's front in the direction x 1 defined by the crack. In the above development, ,t denotes the time partial derivative. The front's velocityȧ is obtained by applying the total time derivative in the local co-ordinates of the crack's tip: Using Equations (37) and (38), we retrieve the M-integral presented in Reference [16].N o w , if we choose the auxiliary field to be the actual field, we obtain: ( , l 1 ) = kl u k,l − u kuk 1j −2 ij u i,1 −2ȧ u i u i,1 1j n j l 1 ds With these assumptions, taking the limit of when goes to the crack's tip, the result is 2J , where J is Rice's J -integral (see Reference [31]), which is equal to the dynamic energy release rate G and, for plane strain conditions, can be related to the dynamic stress intensity factors by the following dynamic equivalent of Irwin's relation (see Reference [25]): In the above equation, K dyn i are the dynamic stress intensity factors defined by References [18,22] and f i are universal functions defined (with c 1 and c 2 , respectively, the dilatational and shear wave velocities, k = 3 − 4 for plane strain) by The interaction integral is defined as the limit of when 1 goes to the crack's tip.
Choosing the virtual crack's extension field q as defined above enables us to write the interaction integral I int as a domain-independent integral. If we define q by we obtain:

Calculation of stress intensity factors
The objective of this section is to develop a domain-independent integral enabling us to estimate stress intensity factors. The expression obtained above does contain a domain-independent integral. Thus, in this subsection, we will focus on the estimation of stress intensity factors based on this integral. In the previous subsection, the virtual crack's extension field was rewritten in order to transform the mixed path, domain-independent integral into a domain-independent integral I int . In the particular case where the auxiliary fields are taken to be the real fields, we also related to Rice's J -integral. Then, the J -integral was related to stress intensity factors using the dynamic energy release rate and the dynamic equivalent of Irwin's relation (Equation (40)). Here, if K aux i denotes the Mode-i stress intensity factors corresponding to the auxiliary fields, we can write: Consequently, K For numerical implementation purposes, the interaction integral is calculated using a J -domain (see Reference [24]). This J -domain is an additional mesh (see Figure 5) onto which the fields are projected, which is used only for the numerical space integration of the interaction integral. The virtual crack's extension is taken to be x 1 at the crack's tip and decreases to 0 on the boundary of the J -domain. For example: where x and y are the local co-ordinates of the crack's tip and d is the size of the J -domain. For curved cracks, q is curved with respect to the cracks' faces according to the assumptions made in the previous development.

Three-point bending
First, let us consider the problem of a stationary crack in a three-point bending specimen whose geometry is shown in Figure 6. The loading is pure Mode-1 and the dynamic stress intensity factor oscillates between 0 and twice the static value given in Reference [18]. This value is where S(0.04 m) is the distance between the supports, a(0.005 m) the crack's length, L the total length of the beam, W(0.01 m) its height and B(1.0m) its thickness. The loading consists in a constant stress 0 (400 Pa) applied over a length l( 0.0025 m). The function is given by The material's properties are E = 200 GPa, = 0.3 and = 7860 kg m −3 . The numerical results were obtained using the mesh shown in Figure 8, which consists of 20 × 100   quadrangles. The time window was 250 s divided into 200 time steps. In this figure, the subdivided elements (whether cut by the crack or including the crack's tip) are also displayed. Subelements were used only for the numerical integration of singular and discontinuous shape functions (see Reference [5]). The results obtained for different sizes of the J -domain ( Figure 7) are shown in Figure 9. K 1 is normalized by the static value and oscillates between 0 and 2 with good accuracy. One can also observe that the influence of the size of the J -domain over the calculation of K 1 is small.

Semi-infinite plate with a stationary edge crack under mixed-mode loading
A schema of the problem is shown in Figure 10. An analytical solution was obtained by Lee and Freund [32].  This example was calculated using two different meshes. The coarse mesh used 40 × 60 linear quadrangular elements and the fine mesh 80 × 120 elements. The accuracy of the results was satisfactory with both meshes. The fine mesh seemed to be more sensitive to oscillations associated with the time integration scheme. Nevertheless, we can conclude that the method does not depend on the mesh size.

Semi-infinite crack in an infinite plate subjected to a tensile stress wave
We choose as our third example the case of an infinite plate with a semi-infinite crack because its theoretical solution is known (see Reference [22]). Several authors already treated this example and, therefore, we can also compare our results with theirs. Under the given assumptions (an infinite plate with a semi-infinite crack) and for the geometry described in Figure 12 the crack's tip. When this situation occurs, the Mode-I stress intensity factor for a stationary crack can be written as For a moving crack, we can write: where k is a universal function of the velocity of the crack's tipȧ which we approximate by Finally, we have: This theoretical solution will be compared with our numerical results. We will be looking at three cases: the crack does not propagate; the crack starts propagating at a prescribed constant velocity v 0 as soon as the wave reaches the crack (t = t c ); and the crack propagates at a prescribed constant velocity v 0 after time t = 1.5t c . The following numerical results were obtained with the plate dimensions H = 2m, L = 10 m and l = 5 m and the material properties E = 210 GPa, = 0.3 and = 8000 kg m −3 . The tensile stress 0 was 500 MPa. The dimension d of the J -domain was 0.5 m and it contained 36 elements with 16 Gauss points each. v 0 was 1500 m s −1 and the stress intensity factors were normalized by 0 √ H . The solutions were calculated using 40 × 80 quadrangular elements with linear approximation. We stopped the computation at time T = 3t c when the stress wave reflected from the boundary reached the crack's tip and the analytical solution was no longer valid. For the stationary crack case, T = 3t c was reached after 200 time steps. The results are plotted in Figure 13. These results match the analytical solution with great accuracy. For the second case, the number of time steps was 20 and the results are not as accurate as in the first case. Oscillations appear at time t = 2.5t c . These oscillations were also observed by the authors of several References [12,13,23]. We may assume that although the analytical solution was valid until time T = 3t c the reflected wave reached the J -domain before that and affected the calculation of K dyn 1 . In the third case, the crack remained stationary until time t = 1.5t c , then propagated at the prescribed constant velocity v 0 . The results are presented in Figure 14. The number of time steps was still 20; the slope was reproduced with great accuracy, after which the oscillations previously observed appeared again. In this case, the solution was calculated using a small time step when the crack does not propagate. This can be justified by the fact that one can assume a steady-state situation around the crack's tip. Consequently, the dynamic effects in this region are less than in the rest of the structure. The crack's velocity is approximately ten times less than the dilatational wave speed. Moreover, as shown in Reference [15], dynamic crack propagation can be simulated effectively using a large time step. This technique enabled us to predict the crack's loading prior to the initiation and, consequently, the initiation itself with a c b ud Figure 15. Definition of the initial geometry. good accuracy, thus avoiding the difficulties due to the high-frequency numerical solicitations when the crack grows. Let us mention that new, enriched degrees of freedom were added and the stiffness properties of the elements were enriched with new shape functions at each time step. This discontinuity (in time) of the stiffness produced a high-frequency response of the structure around the crack's tip, which can explain the oscillations observed when a small time step was used. Here, we are reaching the limits of Newmark-type schemes, which do not deal with velocity or displacement discontinuities very well.

Dynamic crack propagation under prescribed displacement
In order to illustrate the effectiveness of the methods presented above from an energy point of view, we calculated the dynamic propagation of a crack in a plate subjected to a prescribed displacement (see Figure 15). The specimen was made of a homogeneous, isotropic, linear elastic material. It was subjected to a prescribed vertical displacement u d = 0.0025 m at its ends. The evolution of this loading followed the Heaviside step function with a rising time of 0.1 ms. The additional data were c = 0.01 m,b = 0.05 m,a/b = 0.4 for the geometry and E = 186 GPa, = 0.3, = 8000 kg m −3 ,K 1c = 110 MPa √ m for the material's properties. The criterion used for calculating the velocity of the crack's tip can be written as follows: This example had already been calculated in a classical FEM framework in Reference [15].
In Figure 16, we compare the total energy introduced by the X-FEM scheme when the crack grows (i.e. when the discretization evolves) with classical FEM simulation. Using the FEM, the balance recovery method used to control step b (see Figure 2) was very useful in controlling the energy, but step a still produced some numerical energy due to remeshing and projections.
Using the X-FEM and the strategy defined in Section 3, the energy balance was ensured during both steps a and b, as shown clearly in Figure 16. The energy introduced using the X-FEM was non-zero only because the discretized energy balance was calculated using the predicted dynamic energy release rate (i.e. using the interaction integral presented in Section 4), which is quite different from the dynamic energy release rate estimated with the energy definition (see Reference [18]).

Edge-cracked plate under mixed-mode impulse loading
This example deals with a numerical simulation of Kalthoff's experiments presented in Reference [33]. By modifying the projectile's velocity, one observes a transition in the type of failure. At low velocity, i.e. under low strain rate, brittle failure is observed with a global angle of 70 • . The experimental configuration modelled is presented in Figure 17. The material's properties were those of a maraging steel type 18Ni1900: E = 190 GPa, = 0.3, = 8000 kg m −3 , K 1c = 68 MPa √ m. For brittle failure, the characteristic velocity of the projectile was v 0 = 16.5ms −1 , L was 0.1 m and the initial crack's length was l = 0.05 m.
The results were obtained using the maximum hoop stress criterion [34] for the propagation angle and a constant velocityȧ = 750m s −1 . The overall angle of the final crack's path obtained was about 65 • (Figures 18 and 19) which agrees well with the experimental value of 70 • . These results also agreed with those obtained by Belytschko et al. [14] when the authors used the same criterion.

CONCLUSION
The objective of the approach presented above is to simulate dynamic crack propagation using the X-FEM. The main idea consists in using the same singular and discontinuous enrichment as that presented in the scope of linear elastic fracture mechanics for static problems in Reference [5]. First, we generalized the work of Attigui and Petit [16] to the case of propagating arbitrary 2D cracks. Then, we studied Newmark-type schemes for problems with time-dependent discretizations. The stability conditions obtained from this study show how time-dependent problems can be dealt with using the X-FEM in a general sense. The application to dynamic crack propagation leads to an enrichment strategy which ensures numerical stability and energy conservation. The results obtained for several examples were accurate compared to theoretical or experimental data. For stationary cracks under dynamic loading, the method has the same advantages as in the static case (see previous works about the X-FEM). For moving cracks, this approach is proved to be stable and to conserve the energy. The comparison from an energy point of view of the X-FEM and the FEM with remeshing provides a good illustration of the great interest of the approach presented. This method avoids the difficulties encountered with remeshing and projections, which sometimes generate uncontrolled energy transfers. The main effort in our study was devoted to the improvement of the spatial representation of discontinuities. A method was also provided to assess the stability and energy conservation properties of the time integration scheme with time-dependent discretizations. This demonstration was based on a Newmark-type scheme. In dynamic fracture mechanics, a moving crack also implies time discontinuities which are not properly modelled by Newmark integrators. Another time integration scheme could be preferred. The enrichment strategy could be generalized to arbitrary moving discontinuities and the method for arbitrary 2D cracks is expected to be applicable to three-dimensional applications using the level set concept.

A.1. Interaction Lagrangian
Using the notations of Section 4 and considering that each displacement field (u,u aux and u tot ) is kinetically admissible, the variation L int of the interaction Lagrangian for a variation a i of the crack's length in each direction x i is zero. Let us define three configurations and their associated descriptions: Lagrangian, Eulerian and arbitrary Lagrangian Eulerian (ALE). The Lagrangian configuration 0 is the initial configuration, whose description (co-ordinates X i0 ) corresponds to the initial location of the crack's tip. The Eulerian configuration corresponds to an arbitrary extension l of the crack and its description (co-ordinates X i ) moves along with the crack's tip. The ALE configuration¯ and its description are chosen arbitrary between the Lagrangian and Eulerian configurations. Following the notations for the variations of the displacement field between the three configurations defined in Figure A1, one can write: First, the method is developed in the ALE configuration using the Eulerian description. In Equation (A2), the variation of the interaction Lagrangian is expressed as a function of the Figure A1. Lagrangian, Eulerian and ALE configurations 0 , ,¯ and their descriptions X i0 ,X i ,X i . true and auxiliary displacement and velocity fields l int = W int − T int : Using the relation among the three configurations, we obtain: If denotes the variation between the Eulerian and Lagrangian configurations, we can write in the Eulerian description = */*x k a k .
Then, we can use the Gauss-Ostrogradsky theorem and the derivation rule for a product to deduce from Equation (A4): where * is the boundary of the closed domain and n the outward normal to that domain. Now, we define the ALE configuration by choosing the arbitrary configuration to be the Eulerian configuration. Thus, we have ũ = 0 and û =− u. If we use the Gauss-Ostrogradsky theorem again, we obtain the following expression: Finally, we obtain the following expression: − L int = * aux kl u k,l − u k,t u aux k,t mj − aux ij u i, m + ij u aux i, m n j a m dS + aux ij, j u i, m + ij, j u aux i, m + u aux i, t u i, tm + u i, t u aux i, tm a m d (A7) Let us define P int and Q in order to simplify the expression of L int .
P int mj = aux kl u k,l − u k,t u aux k,t mj − aux ij u i, m + ij u aux i, m Q m = aux ij, j u i, m + ij, j u aux i, m + u aux i, t u i, tm + u i, t u aux i, tm

A.2. Path-independent integral
The problem is now considered to be two-dimensional, so that d = b dS and dS = b ds. The closed domain S represented in Figure A2 has an infinitesimal thickness and its contour is *S = 1 ∪ 2 ∪ S + ∪ S − . Then, considering the contours and the domain represented in Figure A2 and assuming that L int is zero, that l is collinear to the tangent at the crack's face along M + or M − and that the crack's faces are traction-free in the actual and auxiliary solutions, we obtain the following equation (where n j are the components of the outward normal vector to domain S): 1 P int kj l k n j ds − 2 P int kj l k n j ds + This equation shows that ( , l) = P int kj l k n j ds + A( ) Q k l k dS does not depend on .

Remarks
(i) Since does not depend on , using three independent vectors l, one can write a local equation linking P int to Q: div P int = Q (A9) (ii) Equation (A9) enables us to simplify the final expression and is a consequence of the linear elastic material behaviour.

A.3. Interaction integral
In this section, the interaction integral is defined as the limit of when goes to the crack's tip.
Then, the virtual crack's extension field q = l is chosen such that we can write the interaction integral I int as a domain integral. If we denote A = A ( 2 ), we can write the limit of the first term in the equation above as − S P int kj q k ,j dS −→ − A P int kj q k ,j dS (A13) Now, considering that q is 0 on − (S ∪ A ( 1 )), the second term becomes: But q is also 0 on − A; therefore: Finally, I int is written as a domain-independent integral: Using the definitions of P int and Q, and the local equation (Equation (A9)), the expression of I int becomes: