The extended finite element method for dynamic fractures

A method for modelling arbitrary growth of dynamic cracks without remeshing is presented. The method is based on a local partition of unity. It is combined with level sets, so that the discontinuities can be represented entirely in terms of nodal data. This leads to a simple method with clean data structures that can easily be incorporated in general purpose software. Results for a mixed-mode dynamic fracture problem are given to demonstrate the method.


Introduction
Modelling of arbitrary dynamic fracture still poses a significant challenge.Recently, considerable success has been achieved with what are called cohesive surface models, cf., [1][2][3].These models could more accurately be called edge-separation models, for they achieve simplicity by allowing the cracks to propagate only along element edges.Consequently, the amount of the dissipated energy during fracture can be substantially greater than in reality because of the limitations on the crack path.
In this paper, we describe the application of a new method that allows arbitrary crack paths.This method is an example of a general framework for approximating discontinuities independent of the finite element mesh.We call the framework the extended finite element method (XFEM).An initial form of the method is reported in [4,5]; the methodology has recently been generalized in [6].The approach is based on a local partition of unity [7,8].
The paper is organized as follows.Section 2 describes the governing equations, the crack representation and the motion in the presence of a crack.Section 3 describes the finite element approximation of the motion by XFEM.Section 4 describes the cohesive laws.In Section 5 we give the weak form and the discretized equations.Then we will describe numerical studies in Section 6.Finally, Section 7 provides a summary and some concluding remarks.

Governing equations and motion
Consider a body Ω 0 in the reference configuration as shown in Fig. 1.The material coordinates are denoted by X and the motion is described by x = ϕ(X, t) where x are the spatial coordinates.In the current configuration the image of Ω 0 is Ω.We define a crack surface implicitly by f 1 (x) = 0. To specify the edge of the crack we construct Fig. 1.A two dimensional body with a crack and its representation in the initial and the current domain.another set of implicit functions f 2 (x) of which gradient is orthogonal to the gradient of f 1 (x) (Fig. 2).The crack edge is given by f 1 (x) = f 2 (x) = 0; f i (x) is only defined in a subdomain about the crack denoted by Ω c .We also define f 0 1 (X) and f 0 2 (X, t) where f 0 i (X) is the image of f i (x) in the initial configuration.The surface f 0 1 (X) = 0 corresponds to the surface where the material has failed, i.e. the crack surface, and is also denoted by Γ 0 c .We arbitrarily choose one side of the crack to be the direction of the normal to the crack, m 0 and choose the sign of f 0 1 so that The function f 0 2 (X, t) is defined so that f 0 2 (X, t) > 0 on the crack surface and vanishes on the crack edge.The discontinuous approximation is constructed in term of the level set functions.Any implicit definition of the surface can be used, but generally we employ the signed distance function defined by f (x) = min x∈Γc x − x , where Γ c is the surface of discontinuity.The level set field is approximated by where N I (x) are finite element shape functions or other interpolants.Generally, the same shape functions are used as for the approximation of the dependent variable.However, other functions can also be employed.
As a consequence of the nodal definition of the geometry of the discontinuity surface, the surface of discontinuity and everything that needs to be known about the discontinuity can be represented by nodal values.This provides great simplicity to the data structure.By contrast, in surface representations of cracks in three-dimensional problems, it is necessary to construct a surface mesh.The task of matching a three-dimensional mesh to a growing surface mesh is quite overwhelming.These new methods, by relying only on nodal values to describe the discontinuity, simplify the data and make updates of evolving surfaces more straightforward.
The strong form consists of the following where P is the nominal stress, ρ 0 the initial density and b the body force per unit mass, τ 0 is the cohesive traction across the crack, t 0 is the applied traction on the Neumann boundary Γ t and u is the applied displacement on the Dirichlet boundary Indicial notation is used for any lower case indices, so repeated subscripts imply summations.We can also write the momentum equation and the traction boundary condition in terms of the Cauchy stress as follows ) where ti is the prescribed traction on the current surface, and Γ c is the image of Γ 0 c in the current configuration.We consider only rate-independent materials.The Truesdell rate of the Cauchy stress can then be related to the rate of deformation D kl by where the tangent modulus C may depend on the stress σ and other state variables.
A key step of this new method is that the discontinuous part of the displacement field can be expressed by where H is the Heaviside function given by and Ψ(X, t) is a function satisfying the following conditions.
The first condition enforces continuity of the crack opening displacement at the cracktip.The second condition precludes overlap of the crack surfaces.The function f 0 1 (X) is not a function of time.This construction of f 0 1 is motivated by the fact that once a crack develops, the image of that crack configuration in the reference domain is fixed, i.e. any crack surface that has developed is frozen.In the computational method, this implies that the crack surface must be considered separately from the continuum for the remainder of the analysis.The function f 0 2 (X, t) is a function of time to reflect the growth of the crack surface; as the crack advances, f 0 2 (X, t) changes.
The discontinuous part of deformation gradient F is obtained from its definition where δ(•) is the Dirac delta function.Since Ψ(X, t) = 0 at the crack edge, the second term of Eq. ( 17) vanishes.Therefore, the deformation gradient is given by Thus the deformation gradient is infinite on the crack surface.
To treat the mechanics on this surface, we then switch from a continuum description of the constitutive behavior Eq. ( 11) to jump relations.We write these relations in rate form in terms of the tractions on the current surface of discontinuity: where D c ij is a constitutive matrix for the interface and τ ∇ i is an objective rate of the traction [see 14; P. 585].In the numerical examples in this paper, we neglect the tangential components of the traction, i.e.only the normal component is considered.

Finite element approximation
In this section, we describe the enriched finite element approximation that models the discontinuity across the crack interface.The crack progresses element-by-element, i.e. the crack advances one complete element at a time.The displacement approximation, i.e. the approximation to Eqs ( 12) and ( 13), is given by where q I are the enrichment parameters and Φ I is where f (X) ≡ f 0 1 (X) = 0 defines the crack interface as shown in Figs 1 and 2 (we drop the subscript and superscript on f ), S is the set of enriched nodes, and n is the total number of nodes.The shifting of the enrichment function limits the enrichments to the cracked elements [11,13].
We will also use the following compact expression for the displacement instead of Eq. ( 20): where The above is very similar to a partition of unity.However it is not strictly a partition of unity, but is instead a local partition of unity [7]: only the set of elements around the crack are enriched.On the elements adjacent to these, the set of functions Ψ I (X) vanish.We next describe the procedure for selecting the set of nodes S to be enriched for a mesh of constant strain triangular elements.The procedure depends on whether we are modelling the nucleation of a discontinuity in an element not contiguous to an existing crack, or the propagation of a discontinuity.We assume here that a law that models both crack nucleation and crack propagation is available.
Consider a mesh of three-node, constant strain triangular elements.The nodes to be enriched are selected as follows Fig. 3(a).Let the nodes of a generic element e 1 be I, J, K. Let the normal to side JK be denoted n I ; it is the normal to the side opposite to node I. Suppose that the material in element reaches to a certain fracture  criterion -for example, the stress exceeds the tensile strength of the material -and element e 1 is not contiguous to any existing cracked elements.The specific node I to be enriched is then given by I = arg(max i.e. the node I is the node opposite to the side that is most closely aligned with the normal to the incipient crack m.
When the enrichment is added to node I, the displacement fields in the adjacent two elements also become discontinuous, as illustrated in Fig. 3.In the two elements denoted as e 2 and e 3 , the discontinuity varies linearly and vanishes at the edge of the elements as shown in Fig. 3.
The procedure for selecting the subsequent nodes to be enriched depends on whether or not the fracture criterion is triggered in an element adjacent to an element with a discontinuity.When the element is contiguous to a failed element, the node to be enriched depends on the direction of the crack, i.e. the normal obtained by the maximum hoop stress criterion [15].The next node to be enriched is selected as follows.From the normal to the crack, we know the direction of propagation and hence which edge of the adjacent element the crack will impinge.The node to be enriched, i.e. the node added to the set S, is then the node opposite to that edge; this is illustrated for two cases in Fig. 3(b,c).

Cohesive law
The mechanics of the crack will be treated by a cohesive model.In a cohesive model, the energy dissipated due to crack propagation is modelled by a traction-displacement jump relation at the crack surface [16][17][18][19][20].
We only consider cohesive laws involving the normal component of the traction.The displacement jump δ N is defined by The normal traction is defined by We will use the linear cohesive law of the form depicted in Fig. 5(a) with a radial return algorithm as in [1,13,21,22].When the current value of the traction in an element gets equal to or greater than the cohesive strength τ max the material switches from continuous to discontinuous.
Given the shape of a cohesive law, the cohesive crack model is completely defined by the fracture energy G F and the cohesive strength τ max .Therefore, the critical crack opening δ max is calculated from G F and τ max .For example, for the linear cohesive model shown in Fig. 5(a), the critical crack opening δ max is The traction law for the case shown in Fig. 5(a) is formulated as an incremental law with radial return by the following algorithm Note that if δ n+1 N 0, i.e. whenever the crack surfaces close, the tractions τ N are computed by the equations of motion; the cohesive law is only active when the crack surfaces are separated.

Variational principle
The trial and test functions reside in the following spaces, respectively: The weak form of the momentum equation is: for u ∈ U and ∀u ∈ U 0 , Here, we have both integrals in the current and reference configurations for convenience.The discrete equations are developed from a weak form.For the case of continua, the weak form is the principle of virtual work.The discrete equations are obtained by standard Bubnov-Galerkin procedures.We let the trial displacement field Eq. ( 22) be where Ψ I (X) have been defined in Eq. ( 23).The test displacement field is In taking any derivatives, we consider Ω 0 the open set with Γ 0 c excluded.Furthermore, we carefully distinguish between Γ + c and Γ − c in evaluating the surface integrals of the crack.The derivatives of the trial displacement field are then given by A similar expression holds for the derivatives of the test displacement field ∂δui ∂Xj .Substituting the above into the weak form Eq. ( 30) gives for all δu iI and δq iI that are not constrained by displacement boundary conditions or the inequality Eq. ( 16); the two parentheses give the discrete equations.In the above, Note that the cohesive traction only effects the nodal force Q ext iI , and that it is evaluated in the current configuration.As noted in [14], the nodal forces can be evaluated in either the current or the reference configurations; it is often convenient to evaluate surface integrals over the current configuration, volume integrals over the initial configuration.
It might be noted that the change of the enrichment during a time step is neglected for simplicity.Otherwise, the inertia term in Eq. (34) would be nonlinear.Although this may introduce some inaccuracy, this simplification was also adopted in meshless methods [23].
The mass matrix corresponding to the standard degrees-of-freedom u iI will be lumped, i.e. diagonalized, by a row-sum procedure.The consistent mass matrix is retained for any terms linked to the enrichment parameters q iI , because we have been unable to determine an effective diagonal mass matrix for these degrees-of-freedom.

Quadrature
One of the major drawbacks of the method is that the integrands in elements cut by discontinuities are discontinuous.For elements containing the crack, the integrals in Eqs (35a) to (37c) can not be evaluated by standard quadrature methods since the integrand is discontinuous.The following quadrature procedure was used, similar to [5].For an element completely cut by the crack, we first divided the element into three triangular subdomains by a Delaunay triangulation for five points -the three nodes and the two crossing points of the crack.Then for each triangular subdomain, we used a standard quadrature method to evaluate the integrals over the subdomain.The integrals were then obtained by summing over the three subdomains.Recently, It has been tried to use the approximated step functions representing the cracks to avoid the subdivision for the numerical integration [24,25].

Dynamic mixed-mode fracture experiment
By using a modified Charpy instrumented impact testing system [26,27], examined are the crack paths and crack initiation angles of three point bending specimens under impact loading.The geometry of the specimen is shown in Fig. 7.The dimensions of the specimen are as follows: height h = 76.2 mm (3 in); length L = 228.6 mm (9 in); distance between supports L s = 2 l 1 = 203.2mm (8 in).
To induce a mixed-mode fracture, a notch n 0 of 19.0 mm (0.75 in) was made at an offset l 2 from the midspan in the specimen.A dimensionless parameter γ was defined to describe the location of the offset notch.
They reported the crack paths and initial crack propagation directions for test specimens with the offset notch at various locations ranging from γ = 0.5 to γ = 0.75.
The initial crack propagation angles observed in the experiments vary from 22 • at γ = 0.5 to 30 • at γ = 0.72.These values roughly agree with values given by linear elastic fracture mechanics.They found a transition point in the notch location a t which changes the mode.For γ > γ t , failure occurs at the midspan; a crack is initiated at the midspan and then propagates.For γ < γ t , the specimen fails at the offset notch.Near the transition point, both midspan crack and crack growth at the notch occur.The experimentally obtained value of a t is about 0.7.
For the loading rates, the material properties are practically strain-rate independent; the change of the Young's modulus and the fracture energy due to strain-rate effects was less than 10%.

Simulation
The experiments of [27] were simulated to see whether we could reproduce the observed crack paths and the transition point γ t , where failure changes from growth at the notch to two-crack growth.For these simulations, we used a cross-triangular mesh constructed from a 24 × 73 rectangular mesh.Constant strain triangular elements were used.To initiate crack growth at the center, a small initial crack (representing a small flaw) was placed at the center of the bottom surface.The length of the flaw n m was taken as 4.76 mm (3/16 = 0.1875 in), which is one fourth of the maximum aggregate size of the concrete.We allowed the cracks to grow only from the notch and the flaw for simplicity, although the method does not restrict the crack growth in this way.The bulk of the concrete is modelled as elastic while the behavior of the crack is modelled by a cohesive crack model.The elastic material properties are chosen as: mass density ρ = 2400 kg/m 3 ; Young's modulus E = 31.37Gpa; Poisson's ratio ν = 0.2.Therefore, the dilatational wave speed is c d = 3810.9m/s, the shear wave speed c s = 2333.7 m/s and the Rayleigh wave speed c R = 2119.8m/s.Here, the rate effect to the material properties is neglected because of the experimental observation that it was insignificant.
The tensile strength τ max is assumed to be 1/3000 of Young's modulus (10.45 Mpa) and the critical crack opening δ max = 0.00375 mm, so the fracture toughness K Ic = 0.8 Mpa √ m which corresponds to a fracture energy G F = 19.58J/m 2 .The direction of the crack advance is determined by the maximum hoop stress criterion.
The velocity boundary condition at the loading point is chosen to be the same as in the EFG simulation by [28].The velocity was imposed as: where v 1 = 0.06 m/s and t ramp = 196 µs.Three runs with different values of the edge notch location γ are made with the cross-triangle mesh.The total simulation time is 1500 µs with Courant number of 0.1.

Results and discussion
The crack evolutions are shown in Figs 8 to 10.As can be seen from the figures, the feature that the transition of the failure mode depends on the location of the notch is captured quite well.As the notch moves away from the center of the beam, the crack pattern obtained by the simulation changes as in the experiment.After several simulations, the transition location γ t was to be γ t = 0.635, which is close to the experimentally observed value 0.7.To show how well the method models the transition of the failure mode, the crack patterns γ = 0.63, 0.635, 0.64,  angle reported in the experiment.The crack from the notch is a little bit toward the loading point during the final stage of crack growth, which is also observed in the experiment (as shown in Fig. 8).
The crack propagation speeds obtained in the simulation are shown in Fig. 11.The ratio of the crack speeds to the Rayleigh wave speed c R ranges from 0.30 to 0.50.Because the crack propagation speed was not measured in the experiment, we can not verify the speed quantitatively.Our computed crack speed is quite high compared to the crack speed 75-115 m/s −1 observed by [29].On the other hand, the computed crack speed ratio to the Rayleigh wave speed is of the order observed in many materials (cf., [30]).
The possible reasons for the discrepancy in the crack speeds might be that we used a linear constitutive law for the concrete bulk.The method is not limited to linear constitutive models, but a nonlinear law has not been added yet.
Although the model still needs some refinement, the simulations show good qualitative agreement with the important features, such as the crack patterns and the transition of the failure mode from one to two crack growth.

Conclusions
A method for simple treatments of arbitrarily growing cracks has been developed.This method does not need remeshing to model growing cracks; the crack paths are independent of the mesh.It is based on the local partition of unity, which makes the enrichment local and increases computational efficiency.
In this method, the extended finite element method is coupled with the level set method.This enables to represent the location of the crack implicitly.The level set function is discretized by the standard finite element method.Therefore the method requires only nodal values of the level set function to represent a crack.
The discretized equation of the method is obtained by the standard Galerkin procedure.The crack is taken into account by a discontinuous enrichment function.Because the discontinuity occurs only in the level of the element which is enriched, the equations are sparse.However, the discontinuity does require the subdivision of the element for the numerical integration of Galerkin weak form.The subdivision does not increase the cost significantly.
An example of dynamic fracture has demonstrated the capability of the method to handle the growing cracks without remeshing.The important experimental observations of the mixed mode fracture under impact loading are well captured by XFEM.

Fig. 2 .
Fig. 2. A two dimensional crack represented by two sets of implicit function f 1 and f 2 .

Fig. 3 .
Fig. 3. (a) Crack formation in an element, in which nodes to be enriched depend on the direction of the crack as shown in (b) and (c)..

Fig. 4 .
Fig. 4. Enrichment scheme for crack initiating at the edge of a body.

Fig. 5 .
Fig. 5. Cohesive models, in which the area under a cohesive law is the fracture energy G F ; (a) a linear and (b) a general nonlinear cohesive law.