A stable numerical scheme for the finite element simulation of dynamic crack propagation with remeshing

This paper presents a general study of the stability of variable-mesh dynamic calculations using an energy approach. This study, whose scope is limited to the calculation of dynamic crack propagation with remeshing, enables us to establish the conditions which are necessary to ensure stability and allow control of energy transfers during the evolution of the mesh. The problem of the implicit calculation of the crack length is also presented. The results obtained on an sample problem are analyzed to illustrate the eﬀectiveness of the proposed methods


Introduction
In many fields, one may have to perform calculations during which the topology of the structure changes for various reasons. In some cases, it is necessary to apply remeshing procedures to represent large deformations of the mesh properly. In a dynamic calculation, because the spatial discretization of the problem varies with time and the time integration scheme requires the knowledge of the state vector at time t n in order to obtain the solution at time t nþ1 , the state vector at time t n must be projected onto the spatial discretization at time t nþ1 , for example to simulate the propagation of cracks. In 2D analysis, one can apply nodes splitting methods, which do not present the problem mentioned above but have the drawback of requiring a priori the knowledge of the path followed by the crack. The other methods involve an evolving mesh. Furthermore, numerical instabilities due to the successive projection operations have been reported in [1]. Therefore, later works focused on the quality of these projections. In [2], the authors formally described the need for a change of discretization. Meshless methods were suggested as an interesting alternative to finite elements for the simulation of dynamic propagation [3,4]. However, the authors did not comment on the stability conditions of the time integration schemes, which is a more complex problem in this case than in the case of finite elements methods.
In the present work, the problem of the instability of dynamic calculations with remeshing is taken from an energy point of view. Dynamic fracture mechanics is used as a basis for this study, which is of a completely general nature. First, we will redefine the continuous and discrete reference problems in order to derive a calculation strategy for the crack growth. Then, we will address the problem of evolving meshes in dynamics and attempt to express a stability condition, then to establish a discretized energy balance. To solve the problem, we propose a technique called the ''balance recovery method'', which guarantees both numerical stability and accuracy of the calculation for any type of projection used in dynamic simulations with varying meshes, particularly in the case of dynamic crack growth with remeshing. Finally, we will present an example to illustrate the results obtained with the tools we developed. The aim of this work is not to develop an efficient remeshing procedure with accurate projections. The problem we deal with, is the numerical stability of the scheme and the uncontrolled transfer of energy when projections are used. Consequently, we voluntary use a basic remeshing procedure and the projections are only linear interpolations to study the influence of such handlings on the stability and the energy balance of the simulation.

The continuous problem
Here, we are presenting a general definition of dynamic crack growth problem with remeshing. In this type of calculation, in addition to the usual unknowns, we are dealing with an additional unknown aðtÞ representing a measure of the extent of the crack (its length in two dimensions). One can write the reference problem as follows: • 8M 2 XðtÞ, 8t 2½0; T , hðtÞ¼2 arctan 1 4 where XðtÞ is the domain being considered, oX 1 the boundary on which displacements u d are prescribed and oX 2 the boundary on which forces F d are prescribed; f d are the prescribed volume forces; r and e are the symmetric stress and strain tensors; G and K i are respectively the dynamic energy release rate, dynamic stress intensity factors; C is the Hooke tensor, q the mass density and G c the critical energy release rate; c r is the celerity of Rayleigh waves in this material; _ a and hðtÞ are the cylindrical coordinates of _ a, the cracktip speed in the crack tip coordinate system; U and S are the function spaces associated with the problem and U 0 the vector space of the virtual fields defined by Remark 1. Here, the problem is presented in weak form. It is interesting to recall its underlying local equations in order to pinpoint the boundary conditions being considered ( Fig. 1): Remark 2. Although the domain was defined as time-dependent (XðtÞ), this evolution is carried out through aðtÞ alone. We shall make the simplifying hypothesis that the initial and deformed configurations are the same and assume small perturbations.
Remark 3. We assume the material to be homogeneous with linear isotropic behavior. 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 [5][6][7]. The dynamic energy release rate is equivalent to the dynamic Rice J integral [8] and can be related, in plane strain assumption, to dynamic stress intensity factors by Eq. (7) where f i are universal function of the crack tip speed (see [7]).
Here, the crack tip speed _ a is calculated with a criterion based on the critical energy release rate [4,7] considering that the material toughness is independent of the crack speed [9]. The crack growth direction calculated using Eq. (3) corresponds to the maximum hoop stress [10].

The discrete problem
Prior to discretizing the problem, let us look at the particular case of a two dimensional straight crack solicited in pure Mode 1. The space discretization is performed using finite elements. For the time discretization, we use a Newmark-type scheme with constants c and b for the second-order Eq. (10) and a generalized trapezoidal scheme with constant a for the first-order Eq. (11). As the discretization of the problem changes with time (because of crack growth), the time integration scheme requires that the components of the state vector at time t n be known in order to calculate the solution at time t nþ1 . Therefore, the state vector at time t n must be projected onto the space discretization at time t nþ1 . Let us denote X j i the discretization of a first-order tensor calculated at time t i onto the discretization at time t j . Then, one can define the projection operator P i;j as For the second-order tensors, let us denote the mass and stiffness matrices as As shown in Fig. 2, there are two possible strategies to calculate the solution at time t nþ1 . The remeshing step requires the determination of X nþ1 n . This stage will be detailed in the following section. Now, we can define the discretized reference problem: Discretized reference problem: Given U n n ; _ U n n ; € U n n a n ; _ a n ( find such that: a nþ1 ¼ a n þ 1 ðÀ aÞDt _ a n þ aDt _ a nþ1 : The quantities U nþ1 n , _ U nþ1 n and € U nþ1 n are unknowns of the problem at each time step. We will detail their calculation in Section 3. The resolution strategy presented here is implicit. Indeed, in order to determine a nþ1 (if a 6 ¼ 0), one must know the velocity _ a nþ1 at the end of the time step. To obtain this quantity, one iterates on Loop 1 defined in Fig. 2. Here, we consider first-order scheme coupled with the crack growth equation Eq. (14). This problem is highly nonlinear problem because of the dependency of G which is not known explicitly. The stability analysis of this scheme is delicate and we chose a P 0:5 based on the results for the linear case. a nþ1 ¼ a n þð1 À aÞDt _ a n þ aDt _ a nþ1 ; Remark 4. The discretized problem presented here leads to five equations with eight unknowns, which means that three projections are needed for its resolution.
Remark 5. If one considers only problems in which no crack growth is required, Eqs. (11) and (13) are irrelevant and the problem is much simpler.

Remark 6.
Here, boundary conditions are taken into account through the Lagrange multiplier technique, wherein F incorporates both Neumann and Dirichlet boundary conditions. Remark 7. Dynamic energy release rate and stress intensity factors are computed using Gh and Mh methods developed by [11,12]. Gh was first developed for static cases then extended for dynamic simulations. To evaluate numerically M and G integrals, both approaches convert them in path independent integrals thanks to h field (also called virtual crack extension). Crack auxiliary fields are used in the Mh method for mixed mode separation. We finally obtain (real fields are u and r u , auxiliary v and r v ): 3. The remeshing procedure: stability analysis and discretized energy balance Section 2 enabled us to develop a calculation strategy for the simulation of dynamic crack growth. As a consequence, the geometry and the discretization of the problem depend on time. In this context, we propose a stability study using the energy method (see [13,14]). Then, we will study its application to dynamic fracture mechanics and present a method to improve the stability and accuracy during changes of discretization. At first, we will present the theoretical tools involved for problems with no remeshing , then expand them to the case of problems with time-varying discretization.

Study of the numerical scheme
There are two possible approaches to the resolution of the equations of dynamics: modal superposition and time integration. We have chosen the second approach in which we use a scheme of the Newmark family to perform the integration of the equilibrium equations discretized in time and space. For a linear elastic material subjected to small perturbations, these equations can be written as follows: where U n designates the discretized displacement vector at time t n , _ U n the velocity, € U n the acceleration, M and K the mass and stiffness matrices respectively.
We define the properties of a numeric scheme by the concepts of stability, consistency and convergence. Stability ensures that a perturbation yields a non-increasing modification of the solution. The scheme is called consistent if the local truncation error is bounded by cDt k , where Dt ¼ t nþ1 À t n is the time step and k the convergence rate. Convergence enables to define that the limit of U n as Dt tends toward 0 is the true displacement field u at the time considered.
From now on, we will focus on the study of the Newmark scheme and their stability, as an instability is a sufficient condition of non-convergence. Thus, we will be interested in the energy method which enable us to evaluate the stability, then the quality of the results obtained. The method we are referring to is presented in [13,14]. First, we will review this method in order to detail the notations being used and the consequences they imply.
We consider the Newmark numerical scheme defined by two parameters b and c: Let us use the following notations for the average and the difference of a quantity between times t n and t nþ1 : with the property: Using these operators, one can rewrite Eq. (17) as follows: Let us consider first the equilibrium equations discretized at times t n and t nþ1 . Since the stability of the scheme is not dependent on external loads [15] (as the scheme must be stable regardless of the problem), we will consider a problem with no external loads: By premultiplying these two equations by ½ _ U , taking the difference of the two quantities obtained and using the relations defined by the scheme, we end up with where This leads to the following theorem [13]: Theorem 8. If c P 1 2 and A is positive definite, then € U nþ1 and _ U nþ1 are bounded.
Corollary 9. If, in addition, K À1 exists, then U nþ1 is bounded.
If these propositions are verified, then the scheme is stable. For this to be the case, the eigenvalues x of A must verify: Thus, we get back to the stability conditions of the Newmark scheme: In addition to this stability analysis, one can write a discretized energy balance. It involves numeric dissipation terms due to the scheme. To develop this equation, one uses the same approach as in the stability study to average the equilibrium equations and premultiply them by ½U . One gets the following equation: is the difference of total energy between times t n and t nþ1 , T ½F þ½U T hF i the work of external forces, the remainder being dissipation terms due to the scheme. Remark 11. If one uses the central difference explicit scheme (c ¼ 1 , what remains constant is not the total energy, but 1 2 Remark 12. Under no circumstance one can use the energy balance to prove the stability of a scheme. Nevertheless, it can be used to verify the consistency of a calculation from an energy point of view.

Stability study
Let us now consider the case of a calculation with remeshing. If one extends the method presented in [13,14] to this type of calculation, one can study the stability of the scheme. Let us use the notations defined in Part 2. Working with two different discretizations requires that we rewrite: • The equilibrium equations at times t n and t nþ1 on the same discretization: • The Newmark scheme: • The operators AEae and ½: Using the approach developed previously and observing that n and ½KU ¼K nþ1 nþ1 ½U þ½KU nþ1 n ; we obtain the stability condition: We recognize on the right-hand side the term in c À 1 2 which governs the stability of the scheme for a calculation with constant discretization, whose stability conditions are: However, the presence of a supplementary term involving ½M and ½K voids the stability conditions expressed in the previous case. Indeed, the stability depends on the sign of the right-hand side of Eq. (27). If this term is always negative or zero, then the calculation is stable, whereas if it is always positive the calculation is unstable.
In the particular cases of the central difference scheme c ¼ 1 2 ; b ¼ 0 ÀÁ or the mean acceleration scheme n Þ, whose sign is not known. Only a numerical study in a particular case would enable us to conclude about the stability of the calculation during remeshing. The method presented here enables us to perform the stability study when going from state X nþ1 n to state X nþ1 nþ1 . To treat the problem for the whole time step, one must study the transition X n n , X nþ1 nþ1 . For this purpose, one writes Remark 13. Here, the term which governs stability does not depend on the parameters of the scheme which, therefore, cannot be the culprit in case of instability. Such instability is caused by the evolutions of the mesh. It is for that reason that we referred earlier to an unstable calculation rather than an unstable scheme.

The discretized energy balance
The energy balance can be rewritten with the notations above in the same way as the stability study was. We get Again, these are the equations of a problem with no mesh evolution completed with a term involving ½M and ½K.
Eq. (29) enables one to express an energy balance between the energies of states X nþ1 n and X nþ1 nþ1 . Finally, one must calculate the balance between X n n and X nþ1 nþ1 , i.e.
Remark 14. In this expression, the variation of the total energy between states X n n and X nþ1 n is taken into account. Since this change of discretization takes place at time t n , the external forces produce no work and W ext designates the work done by these forces between X nþ1 n and X nþ1 nþ1 , whose expression is

Application to fracture mechanics: the balance recovery method
Let us now consider the calculation of dynamic crack growth. In this particular case, we know what physical phenomenon causes the change of geometry. Let us assume a crack growth Da 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 on the crack extension to close it. By making the residual defined by (31) vanish, the balance recovery method guarantees the balance of the state vector projected onto the new discretization (assuming that there are no external forces): Thus, we have the following relations:  À2GDa and À2G _ a (G being the energy release rate and _ a the crack velocity) [5]. In our algorithm, the balance recovery stage is performed after the state vector at time t n has been projected onto the new mesh. This enables one to verify the relations: Let us return to the stability study and energy balance described above. The calculation of ÞÀ2G _ a enables one to quantify the instability actually introduced at the transition between state X nþ1 n and state X nþ1 nþ1 . In what follows, 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 D X nþ1 nþ1 ; X n n ÀÁ denotes the complete energy balance obtained as the difference between the left-hand side and the right-hand side of Eq. (30) taking Eq. (34) into account. D X nþ1 nþ1 ; X nþ1 n ÀÁ is obtained as the difference between the left-hand side and the right-hand side of Eq. (29) taking (34) into account and D X nþ1 n ; X n n ÀÁ is defined by One can immediately write D X nþ1 nþ1 ; X n . D X nþ1 n ; X n n ÀÁ is a measure of the quantity of energy introduced or dissipated during the projection operations (from state X n n to state X nþ1 n ) and D X nþ1 nþ1 ; X nþ1 n ÀÁ when going from state X nþ1 n to state X nþ1 nþ1 .
Remark 15. This method enables us to set D 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 remeshing operations.
Remark 16. The balance takes into account the energy required to create new crack surfaces (2GDa). This term appears naturally because of its relation with the work of the force distribution F nþ1 þ [5].
Remark 17. In practice, the residual (31) is zeroed thanks to the use of an iterative method. Thus, the implementation of this balance recovery method into an explicit dynamic calculation code requires the introduction of iterations to ensure equilibrium prior to resuming the dynamic calculation. These iterations can be initialized by a projection of the displacements, velocities and accelerations of X n n . Subsequently, the residual defined by (31) is calculated, yielding the increments DU nþ1 n , D _ U nþ1 n and D € U nþ1 n which verify: A criterion based on the norm of the residual R is introduced in order to end the iterations. This procedure converges rapidly and only a few iterations suffice to verify the criterion, which makes the increase in calculation cost due to this method acceptable. This enables one to determine the three unknowns U nþ1 n , _ U nþ1 n and € U nþ1 n required for the complete solution of the discretized problem.

Dynamic crack growth under prescribed displacement
In order to illustrate the effectiveness of the methods presented above, we calculate the dynamic propagation of a crack in a plate with remeshing (see Fig. 4). The specimen consists of a homogeneous and isotropic linear elastic material. It is subjected to a prescribed vertical displacement of 0.0025 m at the ends. The evolution of this loading is of the Heaviside step function, with a rising time of 0.1 ms. The additional data are c ¼ 0:01 m, b ¼ 0:05 m, a b ¼ 0:4 for the geometry and E ¼ 186 GPa, m ¼ 0:3, q ¼ 8000 kg m À3 , for the material properties. With this example we want to show the effectiveness of the method. The remeshing scheme and the projection procedure are voluntary basic and not efficient if the balance recovery method is not used. This simulation is programmed into the finite element code CASTEM 2000. The remeshing scheme can be described as follows: when the criterion for crack advance (see Section 2) is checked a new crack front is created, the old crack front node is splited and a new mesh is automatically generated (see Fig. 5). Then fields (displacement, velocity and acceleration) are linearly approximated from the old mesh to the new one. Path independent integrals are used to calculate the energy release rate [11,12]. For time integration, the Newmark mean acceleration scheme is used (c ¼ 1 . In a first time we will use a very coarse mesh as shown in Fig. 11, and then we will compare with the results obtained with a fine mesh.

Influence of the balance recovery method
First, we studied the influence of the balance recovery method on the results obtained in the calculation defined above. Fig. 6 shows the evolutions of the crack length in calculations with and without recovery. Without balance recovery, the time at rupture was divided approximately by two (the calculation was terminated when a b ¼ 0:7). Energy balance can be plotted for each of these calculations to explain these differences. In Fig. 7, T is the kinetic energy, W the strain energy. Next the work of the external forces and balance designates the term DðX nþ1 nþ1 ; X n n Þ defined above (see Eq. (36)). Without balance recovery, the energy imbalance is so large that given the scale of the graph one cannot distinguish the energy levels involved at the onset of the calculation. To evaluate the effectiveness of the balance recovery method, we compared the cumulated lack of balance D ¼   13 double advantage of reducing DðX nþ1 n ; X n n Þ and zeroing DðX nþ1 nþ1 ; X nþ1 n Þ for a final energy imbalance which was negligible compared to that obtained with classical remeshing. Indeed, in the absence of remeshing, the amount of energy introduced through repeated remeshing impairs the accuracy of the calculation and raises doubts about its stability.
Regarding stability, we plotted the cumulated instability I ¼ P n i¼0 I i during the propagation of the crack. In the absence of balance recovery, even though the calculation does not diverge completely, one can see on Fig. 9 that instability indeed occurs.
This example clearly shows the quality of the results obtained with the balance recovery method in terms of both energy conservation and stability. Compared to a calculation without balance recovery, these results explain the differences in the evolution of the crack length observed in Fig. 6. Therefore, this method leads to a stable calculation during which energy conservation is respected.

Study of the time integration strategy for the propagation of the crack
We will now consider the time integration of the equation which yields the crack length given its velocity (see (14)). Let us analyze the results obtained using the algorithm defined by the first calculation strategy presented in Section 2. This algorithm is sketched in Fig. 3.
One can observe that if one chooses a ¼ 0 the calculation is identical to that which was made in the previous subsection, using only the balance recovery method and an explicit integration of crack velocity. To make the results even more significant, we stopped the calculations at a b ¼ 0:9. Fig. 10 shows the results obtained by a calculation without balance recovery, a calculation with a ¼ 0 and a calculation with a ¼ 0:6. One can observe that the calculation without balance recovery predicts a rapid failure of the specimen; with a ¼ 0 rupture takes about twice as long and with a ¼ 0:6 the test piece never breaks. Furthermore, we observe a bifurcation in the evolutions of the crack lengths predicted with a ¼ 0 and a ¼ 0:6 beyond 0.1 ms (i.e. when the loading reaches his final value). One can notice that in this approach (see Eqs. (12) and (13)), the time step is the same for the integration of the balance of momentum and the crack velocity. Consequently, we use an implicit scheme for the crack velocity integration because an explicit time scheme (a ¼ 0) can be instable: the time step is larger than the critical time step corresponding to the fastest crack velocity supposed to be the Rayleigh wave speed. When the time integration is implicit (a ¼ 0:6), the scheme is stable and the numerical crack behavior reveals a deceleration when the loading reaches his final value, until crack arrest. Now, let us focus on the influence of the value of a on the results obtained. Fig. 12 shows the evolution of the crack length for different values of a. We observe that the differences are slight as long as a is greater than 0.5. Finally, thanks to the tools developed above, we predict a crack arrest for the test case of a crack growth under prescribed displacement. This result has been theoretically established in by Freund [7] and has also been experimentally observed for DCB specimen by Kalthoff [16].

Dynamic crack growth under prescribed load
In this example, the crack propagation is theoretically unstable, which often happens in experimental configurations. The specimen is a three point bending specimen and also consists of a homogeneous and isotropic linear elastic material. It is subjected to a prescribed load of 25 N/m on the midspan of the beam. The evolution of this load is of the Heaviside step function, with a rising time of 0.05 ms. The additional  for the material properties (Fig. 13). On Fig. 15, we plot the evolution of the crack length for one simulation without the balance recovery method and an other with the balance recovery method. As soon as the crack propagation is unstable in this case, the difference between both results is not significant. The crack propagates at an approximately constant speed and rapidly breaks the beam. Here with the refined mesh (15,000 elements (Fig. 14)) and the use of a small time step, the numerical instability and the numerical energy introduced by remeshing and projections are negligible. Hence in this case, the balance recovery method gives similar results.

Conclusion
In the first section, we redefined the reference problem of dynamic fracture and identified an appropriate calculation strategy for dynamic crack growth. Then, we focused on the numerical stability of dynamic calculations with varying meshes. We proposed a theoretical study of the problem using an energy approach. This study led to the conditions required to guarantee the stability and accuracy of such calculations. Then, we considered the particular case of dynamic crack growth. We presented a balance recovery method and the results obtained on an example. This method provides a general framework to ensure the numerical stability of the calculation and control the transfers of energy during remeshing. It can be extended to any type of problem involving remeshing and used in any calculation code, implicit or explicit. To be even more general, one can think of extending this method to problems in which the type and/or the number of degrees of freedom depend on time. This would be the case, for example, with extended finite elements methods. The results obtained here suggest that the difficulties encountered to simulate crack arrest using the linear elastic fracture mechanics may come from numerical problems occurring during remeshing and calculating the crack extension. It has been shown with the first example that using a coarse mesh and a large time step, if projections are balanced and if the time integration scheme for the crack length is stable, we can obtained accurate results from an energetical stand point. The second example shows that, in most of usual applications using remeshing, the crack length history can be well predicted without any control of the energy balance if the mesh is refined and especially if the remeshing procedure is efficient.