2D X-FEM simulation of dynamic brittle crack propagation

Summary. The application of X-FEM technique to the prediction of two dimensional dynamic brittle crack growth is presented in this paper. The method is known to guarantee exact energy conservation in case of crack propagation and it is applied to the simulation of one dynamic crack propagation experiment submitted to a mixed mode loading and showing stop and restart of a crack.


Introduction
Modeling and predicting dynamic crack propagation remain difficult challenges. There are at least five different methods to simulate dynamic crack propagations with a finite element technology.
The first one is the element deletion strategy which consists in removing the elements which have no remaining stresses or which are fully damaged. This method, currently used in explicit codes, is very simple to implement but it has two drawbacks: very fine meshes are needed and the results exhibit mesh dependencies.
The second one is the remeshing strategy: this type of method cannot preserve energy conservation both on the whole structure and in the crack tip region.
The third one is based on "element edge" cohesive segments: this method is not robust and leads to crack paths which are mesh dependent [1] because the crack paths are forced to follow the edges of the elements.
The fourth is often named the embedded discontinuity method [2]. It allows a discontinuity of the displacement field at the element level to represent the presence of a crack. The discontinuity amplitude is typically governed by a cohesive law.
The fifth is based on X-FEM method [3]. This method may be coupled or not to a level set technique which is used to represent the geometry of the crack [4]. We shall suppose here that the crack geometry is meshed with a set of 2 nodes geometrical segments. Those segments are totally independent on the structural finite element mesh. The crack growth can be either driven by a LEFM law (which can be used only for brittle fracture mechanics) or by a cohesive law approach which can rather easily be coupled with non linear behavior of the bulk material [5].
One of the main difficulties of dynamic fracture mechanics is the uncertainty upon the appropriate choice for models to be applied for dynamic crack propagation simulations and the method to obtain the corresponding material parameters. Some experiments are already available and used extensively for the validation of numerical predictions: the Kalthoff plate experiments [6] and the Rittel-Maigre ones [7,8]. These experiments show complex crack paths but no arrest and restart. The object of this paper is first to explain a new numerical scheme based on extended finite element method and second to apply it to the simulation of a very carefully designed experiment which provides precise and reliable data on dynamic brittle crack propagations. The paper is limited to bi-dimensional experiments and analysis. The numerical simulation of dynamic crack growth is a useful tool to understand the physics of the dynamic crack propagation. For instance it is difficult to explain from the measures alone the reasons why the crack changes direction or stops during the propagation. The numerical simulation helps to understand such complex questions: for instance the stop and direction changes can be explained by stresses wave reflections during the transient event, which are very difficult to obtain experimentally. The paper is organized as follows: The concepts used for brittle fracture mechanics numerical simulation are first presented.
The use of the X-FEM tool for dynamically evolving cracks is presented in a second section.
The experiment and its numerical simulation are then presented.

Dynamic brittle fracture mechanics
The dynamic crack propagation is a rather difficult subject because the experimental check of the theoretical models is difficult. We shall limit the thermodynamical presentation to a bi-dimensional case and only one crack of length a.

Theoretical background
Let us suppose that the crack propagates at a velocityȧ, the total energy variation writes: In equation (1) W is the total energy which includes kinetic and external loads energy. The second term of equation (1) has a dissipative nature and represents the energy necessary to advance the crack. A detailed and precise theory can be found in [9,10]. The definition of energy release rate G is the following in case of a load which is not applied on the crack lips: In equation (2) Γ is a closed contour surrounding the crack tip, w (resp. τ ) the strain (resp. kinetic) energy density, n 1 the normal to the contour Γ , σ n the traction vector normal to the contour Γ andu the velocity field. This definition can be used for any type of dynamic crack advance and any type of material. Under the hypothesis of linear material behavior and stationary crack tip fields the G integral can be proved to be independent of Γ [9,10] and its expression is: In equation (3) the convention of summation on repeated indices is supposed, δ 1j is the Kronecker symbol, 1 is the direction of the crack tip and n j the outward jth component of the normal to Γ contour. Let us now define the following stress intensity factors: In equation (4), 1 denotes the direction of the tangent direction of the crack at the crack tip, ν the Poisson's ratio and μ the second Lame coefficient. The intensity factors are related through the following relations: where , In equation (6), the indices i are I or II, κ is the Kolosov constant (1 for plane stress state and 3-4ν for plane strain case), c I is the compression elastic wave velocity and c II the shear wave one. The zero of D function gives the Rayleigh wave velocity which will be denoted c R . The two f i functions have a value of 1 for a null crack speed. Injecting equations (4), (5), (6) in equation (3) one gets: These equations are the extension of Irwin's equation [11] to the dynamic case. If we now want to predict dynamic crack propagation, we have to answer three questions: 1. In which direction will the crack propagate? 2. At what speed will it run? 3. Will the crack propagate?
We shall now limit the presentation to the case where the crack propagation direction is driven by the maximum principal hoop stress. The direction is then given by the following equation: Knowing the possible critical direction θ c , one defines an equivalent stress intensity factor K * I by the following equation: The crack is supposed to propagate at the speed which is such that: In equation (11) K I C is the static fracture toughness. Equation (11) is a nonlinear equation when the crack speed is not zero. These equations are sufficient to define brittle fracture crack propagation.

Numerical evaluation of fracture parameters
The precise numerical evaluation of the stress intensity factors is a rather difficult task. The line integral evaluation given in equation (3) is not easy to compute with finite elements because the best values of stresses are obtained on gauss points within the elements. The best numerical method consists to replace the line integral by a surface integral using the divergence theorem. The evaluation of dynamic stress intensity factors rely on the interaction integral presented in the book of Freund [10] and denoted M integral by Attigui [12]. This integral writes: In this equation S is the surface surrounded by the contour Γ , q is a virtual extension field which must be continuous and tangent to the crack lips and zero out of the Γ contour, (u aux ,u aux , σ aux ) are auxiliary displacement velocity and stress fields which are arbitrary but must be statically and kinematically admissible. Figure 1 shows a typical virtual extension field q and also displays the integration box around the crack tip.
If one chooses the q field module to be 1 at crack tip one has: It is now clear that if one chooses the auxiliary fields to be such that K aux I = 1, K aux II = 0 one will get directly K dyn I from the interaction integral. The same holds for the computation of K dyn II . Let us now compute the dynamic stress intensity factors: the simplest solution is to choose the analytical mode I and mode II solutions as auxiliary stress fields. The computation of these quantities implies the evaluation of integrals which are performed numerically using a domain S defined by a rectangular box filled with a large number of sub-elements using standard Gauss integration schemes. The auxiliary fields are computed using the analytical expressions and the other field values within the sub-element are calculated using the original mesh with its associated shape functions.

Remarks on extension to other cases like different brittle failure mechanism or ductile dynamic fracture
The preceding formulation is only valid for brittle fracture mechanics because of the use of energy release rate G concept. If one would like to use this type of approach for a different critical direction one has to replace equations (9) and (10) by different ones. For ductile fracture one has to come back to the basic definition given by equation (2) and hence to use the stress state close to the crack tip to define the properties of the propagation as it is often used in the cohesive zone approach of the simulation [13]. One may nevertheless observe that in existing cohesive zone models the velocity of the crack tip is governed by the de-cohesion rate.

X-FEM strategy: description and energy conservation properties
We shall describe in this section the special features associated with X-FEM simulation of dynamically evolving brittle cracks. The standard X-FEM method is used here [3] with the usual jump function added into the elements fully cut by the crack and the singular functions for the element partially cut. The following representation of displacement field within the elements is chosen: In equation (14) N i (ξ, η) are the standard shape functions, a i represent the amplitude of the degrees of freedom associated with the jump function H and b ij represent the amplitude of the jth singular functions at node i of the element containing the crack tip. Let us recall that the following 4 singular ψ ij functions are chosen: It must be observed that these extended functions are implicitly linked to the crack tip position (r is the distance to the crack tip and θ the angle with the tangent to the crack tip direction). When the crack propagates statically in an elastic medium, one can forget the previous singular function for the new position of the crack. This is not possible in time dependent problems, such as dynamic crack propagations or elasto-viscoplastic material behavior. For dynamic crack propagation it has been shown [14] that if one simply adds new extended degree of freedom associated with the new crack position and that these new extended displacements and velocities are initialized to 0. Consequently the energy is perfectly kept and the new crack segment is closed at the beginning of the time step. Within the next time step the crack opens and the opening work is exactly equal to the fracture energy necessary to propagate the crack. This is a very important result because it ensures that the change of geometry of the body (crack length change) does not affect the numerical results. The local energy is not altered by the geometry changes.

Time integration
The simplest implementation is to use a standard Newmark scheme to integrate the transient equations. The implicit integration is straight forward. One could prefer to use an explicit integrator because the crack propagation is generally very fast. But the method has to be adapted because the time step tends to zero when the crack passes very close to a node. This is simply due to the fact that the mass of the jump degree of freedom of the node very close to the interface is close to zero. A simple approximation has been proposed in [15] to overcome this difficulty: one replaces the added mass on the jump degrees of freedom (a i ) by the standard translational mass for the discontinuous enrichment. The mass matrix is diagonal and the time step of the explicit extended finite element constructed by this method has been shown analytically to be greater than one half of the standard time step of the non enriched element: this is valid for linear finite elements.

Specific crack propagation strategy
The crack propagation evaluation is a complex task because the equations defining the direction, the equivalent stress intensity and the crack speed are coupled and non linear (equations (9), (10), (11), (13)). The simplest way to proceed is to use an explicit time integration of crack advance with the known quantities at the beginning of the time step. In case of propagation the time step is generally chosen to be such that the crack does not progress of more than one element. Nevertheless the time step has to be large enough for the crack tip not to stay more than two time steps within a single element (because of the enrichment strategy choice that keeps the enriched functions corresponding to the old position of the tip). In case of dynamic explicit computation this leads to a constraint for the time steps which may not be compatible with the CFL condition. This explicit procedure works but sometimes induces instabilities in the crack propagation simulation. It needs to be improved in order to permit robust computations: a possible improvement is the implicitation of the prediction of crack increment.

Numerical interpretation of experimental test
This section describes a specific experiment for dynamic crack propagation. It uses a Hopkinson's bar system, and can be considered as a reference case. This experiment shows very interesting features as complex crack paths and crack arrest and restart.

Experimental set up
The experimental loading system is displayed in Figure 2.  The specimen is a simple rectangular PMMA plate with a wave concentrator (the circular hole) and an initial non symmetric crack. It is placed between two nylon bars with appropriate sections. The sections are close to that of the PMMA specimen (1050 mm 2 ). This permits to avoid the wave reflection on the interface. The thickness of the plates is 15 mm. Careful experiments are performed: they are found to be very repetitive, which proves the good quality of the loading system. The input and output velocities as well as the input and output stress waves have been measured and treated in order to have usable data for the numerical analysis. Table 1 describes the basic material data. The determination of Young's modulus for a PMMA under dynamic load is a difficult task because the elastic properties of this material are strain rate depend-ant. A mean value is proposed here. This value is chosen in order to get the best possible fit between experimental results and simple one dimensional wave propagation model of the experiment.

One hole specimen
The precise geometry is displayed in Figure 3.
The first modeling difficulty is the specification of appropriate boundary conditions. One could mesh the two nylon bars, but these bars are axisymetric and very long to avoid wave reflection. The resulting model would lead to axisymetric plane strain coupling, as well as a large model which would lead to heavy computations. A plane strain model of the specimen alone was used. The right hand side (output) bar is modeled by an impedance condition. This condition is applied on all nodes on a width of 40 mm. The impedance is a linear relationship between the velocity of the interface and the applied force on the interface. The following equation gives the impedance relationship: The second difficulty is the modeling of the moving crack path. Had the crack's path been known, an "element edge" cohesive zone approach could have been  successful, but it cannot be used in order to predict the crack's path prior to the experiment. The mesh, consisting of 1,500 4-node elements, is shown in Figure 4. The number of enriched X-FEM element is changing with time. The loading of the left impacting bar is simulated by an imposed speed on the left of the model. The measured velocity V (t) is applied to the left part of the specimen, on a 40 mm width. The input velocity time history is given by Figure 5.
The computation has been performed with a constant value of K I C (1.33 MPa √ m). The comparison of computed and experimental crack path is displayed in Figure 6. The comparison of crack velocity is given in Figure 7. The agreement is rather good except at the very beginning of the propagation: the crack path is slightly different from the measured one ( Figure 8). An explanation of this difference could be the radius of curvature of the notch, which is really smaller when the crack propagates. Figure 9 shows a photograph of the crack tip just after propagation. It is clear that the COD of the moving crack is much smaller of that of the arrested one. The computation has then been rerun with two different values of K I C . When the crack is at arrest (ȧ = 0) the value is taken as the static value of 1.42 MPa √ m whereas the previous value is kept for running cracks (1.33 MPa √ m). With these values the trajectories as well as crack velocities perfectly coincide with experimental data.

Conclusion
This paper has shown the interest of X-FEM to model dynamic brittle crack propagation. Simple macroscopic models seem to be sufficient to explain the crack complex propagations. Indeed, the results of an experimental dynamic brittle crack growth are compared successfully with numerical X-FEM simulation.