INTERNATIONAL JOURNAL FOR NUMERICAL METHODS IN ENGINEERING Int. J. Numer. Meth. Engng (2010) Published online in Wiley InterScience (www.interscience.wiley.com). DOI: 10.1002/nme.2882

Smoothed nodal forces for improved dynamic crack propagation modeling in XFEM Thomas Menouillard∗, †, ‡ and Ted Belytschko§ Theoretical and Applied Mechanics, Northwestern University, 2145 Sheridan Road, Evanston, IL 60208-3111, U.S.A.

SUMMARY Improvements in numerical aspects of dynamic crack propagation procedures by the extended finite element method are described and studied. Using only the discontinuous enrichment function in XFEM gives a binary description of the crack tip element: it is either cut or not. We describe a correction force to modify the forces to smoothly release the tip element while the crack tip travels through the element. This avoids creating spurious stress waves and improves the accuracy of the stress intensity factors during propagation by decreasing the oscillations. Copyright q 2010 John Wiley & Sons, Ltd. Received 28 August 2009; Revised 26 January 2010; Accepted 28 January 2010 KEY WORDS:

dynamic; crack propagation; XFEM; extended finite element method; discontinuity; explicit

1. INTRODUCTION Several types of methods have been developed to deal with the dynamic crack propagation simulation which avoid remeshing. These include element separation methods [1–4] (often called cohesive surface methods), embedded crack methods [5] and element deletion methods [6]. In addition, meshless methods [7–11] also avoid this kind of problem but they are comparatively slow. One way to deal with fracture is to use FEM with few modifications: remeshing and projection of different fields between the different meshes [12]. However this approach, while suitable for statics, is quite ineffective for dynamics.

∗ Correspondence

to: Thomas Menouillard, Theoretical and Applied Mechanics, Northwestern University, 2145 Sheridan Road, Evanston, IL 60208-3111, U.S.A. † E-mail: [email protected] ‡ Post-Doctoral Research Fellow. § Walter P. Murphy Professor. Contract/grant sponsor: Office of Naval Research; contract/grant numbers: N00014-08-1-1191, N00014-09-1-0030 Copyright q

2010 John Wiley & Sons, Ltd.

T. MENOUILLARD AND T. BELYTSCHKO

In order to adapt the standard finite element methods to dynamic and static fracture, the extended finite element method (XFEM) has been developed. It completely avoids remeshing and projection [13, 14]. Belytschko et al. [15] and Stolarska et al. [16] combined this formulation with level sets to track the crack. In XFEM, a discontinuous enrichment function is used along the crack in order to describe the discontinuous displacement [14]. A cohesive zone is often used to deal with ductile or quasi-brittle fracture. This theory has been applied to finite element methods [17–19] as well as the XFEM [20]. De Borst [21] discussed the numerical aspects of using a cohesive zone model. Besides the cohesive law which delays the crack opening, a better description of the crack tip vicinity can improve the results. Zi and Belytschko [22] presented new crack tip elements for XFEM linked to a cohesive zone. Rabzcuk et al. [23] recently presented an improved crack tip element in the context of the rearranged XFEM basis [24], to describe a better kinematic of the crack tip element when the crack propagates through the element, for both static and dynamic cases. R´ethor´e et al. [25] studied the energy conservation in dynamic crack propagation with XFEM. He showed that an appropriate scheme for the initialization of the additional degrees of freedom ensures the stability in terms of energy. When a new cracked element appears, the additional degrees of freedom are injected with a zero initialization and it shows release stress wave emanating from the crack tip. The continuity in time of the space discretization is not ensured. Subsequently, R´ethor´e et al. [26] proposed a space–time method to deal with the new added degrees of freedom in the XFEM computation. Menouillard and Belytschko [27] presented an approach to deal with the releasing crack tip element when the crack propagation matches the mesh orientation. Belytschko and Chen [28] have studied dynamic cracks with singular enrichments. The results were quite accurate but the method required many quadrature points around the crack tip and is limited to elastodynamic fracture. Menouillard et al. [29] have presented an enriched formulation, including a moving crack tip enrichment. Substantial improvement was observed in the results, but this formulation is more complex. Added to the XFEM formulation, a Meshfree discretization near crack tip also improved the results in terms of stress intensity factors [30]. The elementwise description of the crack propagation in the interelement method and XFEM creates oscillations in the stress field due to the sudden release of the crack tip element. The smooth release of the elements due to the crack propagation has been first developed by Aberson and Anderson [31], Hardy [32], King and Malluck [33] and Malluck and King [34]. These developments focused on mode 1 crack propagation at constant velocity in a regular mesh within the interelement method. However they decreased the oscillations and their numerical results agreed well with the theory. Kobayashi et al. [35] compared the numerical and experimental results for a mode 1 crack propagation at constant velocity, and they showed remarkable agreement. Thesken and Gudmundson [36] applied a variable order singular element located near the crack tip in the interelement method still for a mode 1 crack propagation in a regular mesh, in order to deal with the relaxation of the nodes [37–39], which is based on the control of the stress inside the crack tip elements. This special higher order element had to move with the crack tip. Recently, Fan and Fish [40] proposed a method for propagating arbitrary failure modes using a local refined mesh. This paper describes an improvement in the XFEM methods for dynamic crack propagation based on element-by-element propagation. It introduces a slow release of the crack surfaces in the crack tip element when the crack propagates. It aims at characterizing the release of the crack tip element when additional degrees of freedom are injected due to the crack propagation. A correction force is introduced, and takes into account the fraction ratio of the crack tip element which is cut by the crack, and makes the new degrees of freedom continuous in time. The proposed method is Copyright q

2010 John Wiley & Sons, Ltd.

Int. J. Numer. Meth. Engng (2010) DOI: 10.1002/nme

IMPROVED DYNAMIC CRACK PROPAGATION MODELING IN XFEM

based on the control of the internal forces of the newly added degrees of freedom. Whereas no tip enrichment is involved to describe the crack tip vicinity, the release process decreases the spurious oscillatory release stress waves. The paper is organized as follows. Section 2 presents the governing equations, and introduces the cohesive zone model and the stress intensity factors computation. Section 3 describes the space discretization of the XFEM formulation using only the discontinuous enrichment and the time integration procedure. Then, Section 4 presents the new method for releasing the crack tip element by using correction force: it explains the link between the correction force and the position of the crack tip in the element and the smooth release process when the crack propagates. Section 5 presents the numerical examples to demonstrate the effect of the correction in the stress field and in the stress intensity factors.

2. GENERAL PROBLEM MODELING AND FRACTURE MECHANICS 2.1. Problem modeling The method is applicable to two and three dimensions, but only the two-dimensional (2D) case is described. Consider an initial domain 0 and its boundary as shown in Figure 1. As usual, the boundary is subdivided into displacement boundary *0u and traction boundary *0F , and *0u ∩*0F = ∅. The image of 0 is the current domain denoted by , and the motion is described by x = (X, t), where t is the time, X and x the material and spatial coordinates, respectively. The displacement at the material point X is denoted by u(X, t). Figure 1 shows the reference and current domains, and the different boundaries. The image of the crack in the reference configuration is denoted by 0c , whereas in the current domain  the crack is c . The momentum equation is written as: u¨ = div()+fd

in 

(1)

where  is the density, and fd a body force vector on the current domain ,  the Cauchy stress tensor and u¨ the acceleration. The global equation verified by the displacement solution u ∈ U is

Figure 1. A two-dimensional body with a discontinuity and its representation in the reference and the current domains. Copyright q

2010 John Wiley & Sons, Ltd.

Int. J. Numer. Meth. Engng (2010) DOI: 10.1002/nme

T. MENOUILLARD AND T. BELYTSCHKO

(see e.g. [41]): 



/c

¨ d+ u·v

 : ε(v) d





=

/c

/c

fd ·v d+

* F

 Fd ·v d+

c

f coh ·v d

∀v ∈ U0

(2)

where ε is the strain tensor. Fd is the applied traction and f coh the cohesive traction. The domains U and U0 are defined by: U = {u|u(x) = ud , U0 = {u|u(x) = 0,

∀x ∈ *u , ∀x ∈ *u ,

and regularity of u except on c } and regularity of u except on c }

(3) (4)

The boundary conditions are: u(x, t) = ud

∀x ∈ *u

∀t ∈ [1, T ]

(x, t)·n = Fd

∀x ∈ * F

∀t ∈ [1, T ]

(x, t)·n = fcoh

∀x ∈ c

(5)

∀t ∈ [1, T ]

where n denotes the outer normal of the domain , T is the total time. Moreover, the surface c depends on time as the crack propagates, hence c (t). 2.2. Fracture mechanics Here we review the physical laws that we used to drive crack propagation. First, considering a crack propagation from an energetic point of view, energy is dissipated during the propagation. One way to model this is by a cohesive law. Second, a crack is explicitly updated from the knowledge of the previous crack tip position and the crack velocity. This velocity along the crack front is simply a vector in the 2D case. As we only consider the 2D case in this paper, the complete updating process of the crack tip position requires the direction and norm of the crack velocity. In the framework of the linear fracture mechanics, the stress intensity factors give the velocity of the crack. 2.2.1. Cohesive law. The energy dissipation in fracture will in some cases be treated by a cohesive zone model, which was introduced first by Barenblatt [42] for elastic plastic fracture in ductile metals, and by Hillerborg et al. [43] for brittle materials. Figure 2 shows the crack opening and the tensile traction along the crack. It shows the closing effect of the cohesive traction near the crack tip with a limit opening max ; beyond this opening, the tensile traction vanishes. The relation between the traction and the crack opening is given by the cohesive law. Figure 3 presents a linear cohesive law [24, 44]; the area under the curve represents the fracture energy G F , defined by:  max 1 () d = max max (6) GF = 2 0 Copyright q

2010 John Wiley & Sons, Ltd.

Int. J. Numer. Meth. Engng (2010) DOI: 10.1002/nme

IMPROVED DYNAMIC CRACK PROPAGATION MODELING IN XFEM

Figure 2. Cohesive law described with the crack opening  as a function of the position x along the crack.

GF

GF 0 Pe n

alt

y

Cohesive traction:

max

0

max

Crack opening:

Figure 3. Linear cohesive law: the area in gray under the curve represents the fracture energy G F .

where  denotes the tensile stress along the crack, and max , max define the linear cohesive law described in Figure 3. 2.2.2. Stress intensity factors. To assign the direction and c (t) velocity a(t), ˙ we use the stress intensity factors. The initial value a(0) = a0 is the initial crack length, and the knowledge of the two functions a(t) ˙ and c (t) is sufficient to describe the history of the crack geometry in the structure. As we assume that the material is homogeneous with linear isotropic behavior, the problem falls within the framework of the linear dynamic fracture mechanics. Note that in XFEM, in contrast to interelement crack methods, a criterion for propagation of the crack must be given. The same is true in continuum mechanics. In interelement crack models, this requirement is bypassed by introducing a length scale (the element size) and a limited set of possible directions (the edges). We assume a K -dominant fracture behavior. The computation of the stress intensity factors is based on the auxiliary fields near the crack tip and the interaction integral which allows to write two scalar equations, and therefore determine the two unknown stress intensity factors K 1 and K 2 . Further developments of the interaction integral are presented in [25, 45]. Copyright q

2010 John Wiley & Sons, Ltd.

Int. J. Numer. Meth. Engng (2010) DOI: 10.1002/nme

T. MENOUILLARD AND T. BELYTSCHKO

S1 S2 S1

S2 Crack

Figure 4. Direction and norm of the virtual extension field q.

The stress intensity factors K 1 and K 2 are computed using a domain-independent integral I int and a virtual crack extension q by:   int aux aux ˙ u˙ )div(q) d+ aux : (∇u∇q)+ : (∇uaux ∇q) d I = − ( : ∇u−u· 



 +

aux



div(

)∇u(q)+¨u·∇u

where aux , uaux are auxiliary stress and is defined by: ⎧ 0 ⎪ ⎨ q = q = 1 ⎪ ⎩ q linear

aux

(q) d+

 

˙ ˙ u˙ aux (q) d u˙ aux ·∇ u(q)+ u·∇

(7)

displacement fields. The vector q, parallel to the crack, outside the domain S1 ∪ S2 inside the domain S1

(8)

inside the domain S2

where S1 and S2 are domains defined near crack tip. Figure 4 shows the domains S1 and S2 , and presents the direction and the norm of the virtual extension field q near the crack tip. Note that the shape of the domains S1 and S2 can be circular too. Then this integral is: I int =

2(1−2 ) ˙ 1 K 1aux +2 (a)K ˙ 2 K 2aux ) (1 (a)K E

(9)

where the Kolosov coefficient is 3−4 for the plane strain, or (3−)/(1+) for the plane stress and D, a function whose zero defines the Rayleigh wave celerity denoted by cr , is written as D(a) ˙ = 4 1 2 −(1+ 22 )2 . The parameters i are defined by 2i = 1−(a/c ˙ i )2 with i ∈ {1, 2} where c1 and c2 are, respectively, the dilatational and shear velocities which, √as a function of the Lam´e coefficients and the density, are given by: c1 = (+2 )/ and c2 = /. The different stress intensity factors are computed through an appropriate choice of uaux : i.e. the field K 1aux = 1 and K 2aux = 0 for the determination of K 1 , and the field K 1aux = 0 and K 2aux = 1 for the determination of K 2 . 2.2.3. Crack velocity law. The equivalent dynamic stress intensity factor K  is defined by:   3 c 3 c K  = cos (10) K 1 − cos sin c K 2 2 2 2 Copyright q

2010 John Wiley & Sons, Ltd.

Int. J. Numer. Meth. Engng (2010) DOI: 10.1002/nme

IMPROVED DYNAMIC CRACK PROPAGATION MODELING IN XFEM

Freund [46] developed a relation between the dynamic energy release rate and the crack velocity. The relation between the crack velocity and the equivalent dynamic stress intensity factor is [46]: ⎧ 0 if K 
2   a˙ = (11) K 1c ⎪ otherwise. ⎪ ⎩ cr 1− K  Freund and Douglas [47] and Rosakis and Freund [48] have verified this relation experimentally. Thus, a crack increment in the explicit algorithm is defined by a = at, ˙ where t is the time increment. In addition, we use the maximum hoop stress criterion to drive the direction of the crack: ⎞⎤ ⎡ ⎛   2 K1 1 K1 ⎠⎦ c = 2artan ⎣ ⎝ (12) −sign(K 2 ) 8+ 4 K2 K2

3. SPACE DISCRETIZATION AND TIME INTEGRATION 3.1. Space discretization: XFEM In XFEM, the discontinuous enrichments are added to the conventional finite element approximation. The discretized displacement u˜ is written using an enriched basis as:   ˜ t) = u(x, N I (x)u I (t)+ H ( f (x))N J (x)b J (t) (13) I ∈N

J ∈Ncut

where N is the set of all nodes of the mesh, Ncut the set of nodes belonging to the elements completely cut by the crack. f is the level set function whose isozero defines the crack surface, H is the discontinuous Heaviside function whose value, through the level set f , is +1 on one side of the crack and −1 on the other. The standard shape functions are denoted by N I . The standard degrees of freedom are denoted by u I ; b J corresponds to the additional degrees of freedom for the enrichment to represent the crack discontinuity. Figure 5 presents a mesh, a crack and the corresponding sets of enriched nodes for an XFEM model where the crack grows element-by-element. Note that the nodes corresponding to the edge on which the tip falls are not enriched. The discrete equations are obtained via the weak form, i.e. Equation (2). The test function v is taken to have the same shape than u˜ in Equation (13). Substituting the trial function into the weak form (2) gives the following equation: ¨ J = −f I MI J U ext coh int f I = fint I −f I −f I . f I

∀I ∈ N

(14)

where denotes the internal force at node I , the external force and fcoh I the cohesive force. Note that f denotes the level set function for the crack surface, whereas f is the force vector. The remaining terms are defined below, as standard and additional degrees of freedom at node I :    u,int  u¨ I fI ¨I = U , f Iint = (15) b¨ I f Ib,int Copyright q

2010 John Wiley & Sons, Ltd.

fext I

Int. J. Numer. Meth. Engng (2010) DOI: 10.1002/nme

T. MENOUILLARD AND T. BELYTSCHKO

Discontinuous enrichment

Figure 5. Different sets of the nodes along the crack: N defines all the nodes, and Ncut defines the circles.

 f Iext

=

f Iu,ext f Ib,ext



 f Icoh =

0

 (16)

f Ib,coh

where, ∀I ∈ N  f Iu,int i

=



N I, j  j i d

 f Iu,ext

=



(17)

 N I fd d+

* F

N I Fd d

(18)

and ∀I ∈ Ncut , one has  f Ib,int i

=



N I, j H ( f (x)) j i d

(19)

N I n d

(20)

 f Ib,coh

=

c

 f Ib,ext =



 N I H ( f (x))fd d+

* F

N I H ( f (x))Fd d

(21)

where  denotes the cohesive stress along the crack (see Figure 3). The consistent symmetric mass matrix is given by:  Muu IJ =  Mbb IJ =





N I N J d

N I N J H 2 ( f (x)) d

∀(I, J ) ∈ N×N

(22)

∀(I, J ) ∈ Ncut ×Ncut

(23)

∀(I, J ) ∈ N×Ncut

(24)

 Mub IJ = Copyright q



N I N J H ( f (x)) d

2010 John Wiley & Sons, Ltd.

Int. J. Numer. Meth. Engng (2010) DOI: 10.1002/nme

IMPROVED DYNAMIC CRACK PROPAGATION MODELING IN XFEM

However, the mass matrix was diagonalized as described in [49, 50]. Hence, the lumped mass matrix becomes  (25) MuII = N I d ∀I ∈ N 

 MbJJ =



N J H 2 ( f (x)) d ∀J ∈ Ncut

(26)

3.2. Time integration The well-known explicit central difference method [41] is used: ˙ t+t/2 = U ˙ t + 1 t U ¨t U 2

(27)

˙ t+t/2 Ut+t = Ut +t U

(28)

¨ t+t = f ext −f int +f coh M· U t+t t+t t+t

(29)

˙ t+t = U ˙ t+t/2 + 1 t U ¨ t+t U 2

(30)

˙ t , and U ¨ t ) denotes the displacement (respectively velocity, and accelerawhere Ut (respectively U tion) matrix at time t. t is the time step size, and can vary with the time instead of being constant. For the XFEM formulation, Menouillard et al. [49, 50] developed a mass lumping strategy for the discontinuous enrichment part. He found that the discontinuous enrichment does not decrease the stable time step much; indeed the XFEM critical time step is at least 70% of the finite element one, and 100% if an appropriate rearrangement of the shape function basis is performed [24]. 4. CORRECTION FORCE FOR RELEASING CRACK TIP ELEMENT In this section, we describe the proposed method to deal with the release of the forces from crack tip element enrichment when the crack propagates through. The aim is to avoid sudden tip element release during propagation, and thus avoid unphysical stress waves due to the crack propagation. 4.1. Kinematics of the crack tip element Here, we deal with the kinematics of the crack tip element that is only enriched with discontinuous enrichment as shown in Equation (13). Let us consider the element with nodes {1, 2, 3, 4} in Figure 6. The kinematics is described by Equation (13): it is a sum of a continuous part and the enrichment. Figure 6 shows the enriched nodes as the crack propagates. The newly added degrees of freedom are located on nodes 1 and 2 in Figure 6(b) and, nodes 1 and 4 in Figure 6(c). Figures 7(a) and (b) present the enriched part of the discretized displacement in the crack tip element shown in Figures 6(b) and (c). One can notice that the discretized displacement presents a discontinuity along the cut edge (i.e. edge 1–2), and no displacement discontinuity on any other edges (1–4, 3–4 and 2–3). The displacement field allows to represent a crack opening in the tip element with only the discontinuous enrichment. Thus the strategy to deal with the identification of the new enriched nodes is then finalized. Copyright q

2010 John Wiley & Sons, Ltd.

Int. J. Numer. Meth. Engng (2010) DOI: 10.1002/nme

T. MENOUILLARD AND T. BELYTSCHKO

2

3

2

3

1

4

1

4

(a)

(b)

2

1

3

4

(c)

Figure 6. Crack propagation possibilities in a 4-node element: (a) is the initial crack, (b) and (c) are the different cases at the next step; the circles represent the enriched nodes.

0

0 2

2 3

1

(a)

3 1

4

(b)

4

Figure 7. Kinematic of the enriched part of the discretized displacement: (a) corresponds to the crack tip located on the edge 3–4 (H (N1 b1 + N2 b2 )) and (b) to the edge (2–3) (H (N1 b1 + N4 b4 )).

Remark Another approach is to enrich all nodes of the crack tip element, and enforce continuous displacement along the three uncut edges. The constraint that the crack tip is located on an edge can be written as: ˜ Xtip = 0 [u]

(31)

where Xtip is the position of the crack tip. This constraint is given by enforcing Equation (31) on the field (13): 0 = b4 N4 (Xtip )−b3 N3 (Xtip )

for Xtip ∈ [3−4]

(32)

0 = b3 N3 (Xtip )−b2 N2 (Xtip )

for Xtip ∈ [2−3]

(33)

The minus sign is due to the Heaviside function. This constraint gives an explicit relation between two enriched degrees of freedom (i.e. b3 /b4 or b2 /b3 ). Figure 8 presents the displacement field with two additional degrees of freedom along the edge containing the crack tip. 4.2. The method with correction forces Figure 9 shows a generic situation where a crack tip passes from the edge between nodes 1 and 2 to the edge joining nodes 3 and 4. The configuration when the crack is injected is given in Figure 6(b). Copyright q

2010 John Wiley & Sons, Ltd.

Int. J. Numer. Meth. Engng (2010) DOI: 10.1002/nme

IMPROVED DYNAMIC CRACK PROPAGATION MODELING IN XFEM

0

0 2

2 3

1

3 1

(a)

(b)

4

4

Figure 8. Kinematic of the enriched part of the additional 2 nodes enriched of each side of the crack tip: (a) corresponds to the crack tip located on the edge 3–4 (H (N3 b3 + N4 b4 )) and (b) to the edge (2–3) (H (N2 b2 + N3 b3 )).

3

2

A

1

B

A

4

(a)

3

2

1

B

4

(b)

3

2

A

1

(c)

B

3

2

A

4

1

B

4

(d)

Figure 9. Evolution from the state (a)–(d) through the intermediate states (b) and (c). Mbnew denotes the mass matrix related to the new discontinuous degrees of freedom bnew denoted by the dotted circles. f b the forces related to the degrees of freedom bnew . The other degrees of freedom are only subjected to Equation (14). (b) Mbnew · b¨ new = −f b +f = 0; (c) M bnew · b¨ new = −f b +fcorr ; (d) Mbnew · b¨ new = −f b .

Note that the method is also valid for the case shown in Figure 6(c) but not illustrated in this section, i.e. the following can also be done for the joining of nodes 2 and 3. In the method, we consider a virtual crack that moves through the element with a constant velocity given by Equation (11) and a constant angle given by Equation (12). Before the crack reaches edge 1–2, nodes 1 and 2 are not enriched; when the enrichment is added to nodes 1 and 2, the crack tip passes through edge 1–2. Let the time at which the enrichment is injected at nodes 1 and 2 be denoted by tinj . The equation for the acceleration of the new enriched degrees of freedom, denoted by bnew , is then Mbnew · b¨ new = −f b Copyright q

2010 John Wiley & Sons, Ltd.

(34) Int. J. Numer. Meth. Engng (2010) DOI: 10.1002/nme

T. MENOUILLARD AND T. BELYTSCHKO

Figure 10. Correction force as a function of the fracture ratio of the crack tip element (i.e. the position of the crack tip in the element).

where

 bnew =

b1 b2



 f=

f1b



f2b

(35)

Mbnew is the mass matrix related to the degrees of freedom bnew only. The stress in elements A and B will generally be nonzero when t = tinj . As can be seen from Equations (14) and (19), the internal and external nodal forces associated with b1 and b2 , fb1 and fb2 , respectively, will also be nonzero, and in fact quite large. Therefore, the injection of the two enriched nodes is tantamount to introducing a step function in the nodal forces. Our proposed method aims to alleviate the effect of the step function (i.e. discontinuous) character of the nodal forces. This is accomplished by indirectly ramping up the effect of the existing stress, cohesive tractions and external forces. This is done through a correction force, so that Equation (34) is replaced by Mbnew · b¨ new = −f b +f corr

(36)

where the norm of the correction force f corr tends to zero when the crack tip reaches the next edge, i.e. when the element becomes completely cut by the discontinuity (see Figure 9(d)). All other degrees of freedom are subjected to Equation (14). The value of the correction force at t = tinj (when the crack tip is on the previous edge) is chosen to be f b such that the initial acceleration b¨ new vanishes. An example of the time history of f corr is given in Figure 10; the correction force goes from the initial force f b to zero when the crack tip travels from one edge to the other. Between these two crack tip positions, the correction force is taken to be linear. So the correction force corresponding to the new additional degrees of freedom bnew is given by  a(t) corr corr f (t) = 1− f (tinj ) (37) le f corr (tinj ) is such that: f corr (tinj ) = f b (tinj )

(38)

where f b (tinj ) = f int (tinj )−f ext (tinj )−f coh (tinj ). The decay of the correction force need not be linear; the initial value of this law (f corr , a/le ) is the only parameter of the method. Copyright q

2010 John Wiley & Sons, Ltd.

Int. J. Numer. Meth. Engng (2010) DOI: 10.1002/nme

IMPROVED DYNAMIC CRACK PROPAGATION MODELING IN XFEM

The continuity of the nodal force related to the new additional degrees of freedom gives the same property to the acceleration through the momentum equation. Therefore, the velocity and displacement remain quite continuous when the additional degrees of freedom are injected, and thus the property of continuity in time remains in the strain and stress field too, and a continuous progressive release of the tip element occurs.

5. NUMERICAL APPLICATIONS 5.1. Example Let us consider an example to underline the effect of the correction force on the computation of the response of the crack tip element. We consider a coarse regular mesh with a rectilinear crack. The velocity of the crack tip and its direction are prescribed; the computation is performed with a damping parameter = 2%. This example aims at showing the evolution of the response of the structure as a function of the crack length or time, as the crack length increases with the time. The aim of the correction force is to make the response as smooth as possible, compared with the case without the correction force. The displacement is zero at the bottom of the structure, and a constant tensile stress is applied at the top. The initial crack is on the left edge of the structure and is 2 elements long. Figure 11 presents the geometry and loading of the example: the length of the specimen L = 0.2 m and width H = 0.1 m. The tensile stress applied is 0 = 10 000 Pa. The initial crack length is a0 = 0.04 m. The media is considered linear elastic, and its material properties are: Young’s modulus = 1 MPa, the Poisson ratio = 0.3 and the density = 1 kg/m3 . The crack tip velocity is imposed to be constant = 1 m/s. The response of the structure is examined by the displacement of the top left corner point as a function of time. Two simulations are run with crack directions of 0 and 11◦ . Figure 12 presents the results for a crack propagating at 0◦ in the regular mesh as shown in Figure 12(a): the dashed line is the propagating crack. Figure 12(b) presents the vertical displacement of the top left corner of the structure as a function of time for both computations, with and without the correction force. One can notice that the displacement in the standard XFEM is usually discontinuous due to the propagation of the crack element by element. The correction force (Figure 12(b)) makes the response of the crack tip element continuous, so the displacement of the corner point becomes continuous. Figures 12(c) and (d) show the vertical velocity of the corner point as a function of time for the computation with and without the correction force. Note that the velocities are significantly

Figure 11. Geometry and loading of the illustrative example. Copyright q

2010 John Wiley & Sons, Ltd.

Int. J. Numer. Meth. Engng (2010) DOI: 10.1002/nme

T. MENOUILLARD AND T. BELYTSCHKO

0.14

0.08

0.12

With correction Without correction

0.07 0.1 0.06 Displacement (m)

0.08 0.06 0.04 0.02

0.05 0.04 0.03

0

0.02

-0.02

0.01

-0.04 -0.02

0 0

20

0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0.18 0.2 0.22

40

(b)

(a) 70

2

Without correction

60

80

100

120

Time (ms) With correction

1.8

60

1.6

50

Velocity (m/s)

Velocity (m/s)

1.4 40 30 20 10

1.2 1 0.8 0.6 0.4

0

0.2

-10

0

-20

-0.2 20

(c)

40

60

80

Time (ms)

100

120

20

(d)

40

60

80

100

120

Time (ms)

Figure 12. Test 1: a flat crack with imposed velocity and direction (0o ): (a) mesh and crack path; (b) displacement of the top left corner of the structure as a function of time; (c) velocity of the top left corner without the correction force; and (d) velocity of the top left corner with the correction force.

different. The correction force smoothes the velocity and makes it more physical. The velocity of the corner point without the correction force shows a large acceleration each time the crack tip reaches a new element. Figure 13 presents the same results (i.e. crack path, displacements and velocities) for a crack propagating with an angle of 11◦ so that the crack path does not match with the mesh orientation. The same observations as in the previous test are made in terms of the smoothness of the response of the structure. One can conclude that the correction makes the stiffness of the crack tip element smoother as a function of time. We also underline that the derivative of the stiffness is not continuous when the crack reaches a new element. However, the continuity of the stiffness directly affects the internal forces of the nodes of the crack tip element, and this continuity will bring this same property on the acceleration of these nodes. The release stress waves at the crack tip element during propagation are now expected to be significantly decreased. We remark that a different correction law (see Figure 10) may change the continuity of the velocity and stiffness. We show later that damping alone will not give smooth, accurate results. Copyright q

2010 John Wiley & Sons, Ltd.

Int. J. Numer. Meth. Engng (2010) DOI: 10.1002/nme

IMPROVED DYNAMIC CRACK PROPAGATION MODELING IN XFEM

0.14

0.08

0.12

With correction Without correction

0.07 0.1 0.06 Displacement (m)

0.08 0.06 0.04 0.02

0.05 0.04 0.03

0

0.02

-0.02

0.01

-0.04 -0.02

0 0

0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0.18 0.2 0.22

(a)

20

60

60

100

120

With correction

2.5

40

2 Velocity (m/s)

30 20 10

1.5 1 0.5

0

0

-10 -20

-0.5 20

(c)

80

Time (ms) 3

Without correction

50

Velocity (m/s)

40

(b)

40

60

80

100

Time (ms)

120

20

(d)

40

60

80

100

120

Time (ms)

Figure 13. Test 2: a mixed mode crack with imposed velocity and direction (11o ): (a) mesh and crack path; (b) displacement of the top left corner of the structure as a function of time; (c) velocity of the top left corner without the correction force; and (d) velocity of the top left corner with the correction force.

5.2. Moving semi-infinite mode 1 crack The aim of this example is to show that the proposed correction also has beneficial effects on the stress intensity factor. The example considered is an infinite plate with a semi-infinite crack [51] loaded as shown in Figure 14. A theoretical solution of this problem for a given crack velocity is given by Freund [45]. Since this analytical solution was obtained under the assumptions of an infinite plate with a semi-infinite crack and a given speed of the crack tip, according to the geometry described in Figure 14, the FEM model and the analytical solution are only comparable for time t3tc = 3h/c1 (where c1 is the dilatational wave speed). Beyond that, the reflected stress wave reaches the crack tip and the analytical solution is no longer applicable to the FEM model. The dimensions of the structure are the following: the length L = 10 m, the initial crack length a = 5 m, and the vertical position of the crack h = 2 m. The material properties of the linear elastic media are: Young’s modulus E = 210 GPa, the Poisson ratio  = 0.3 and the density  = 8000 kg/m3 . The tensile stress applied on the top surface is 0 = 500 MPa. The crack velocity is imposed to be zero until √ 1.5tc , and 1500 m/s thereafter. The mode 1 stress intensity factor is normalized by the factor 0 h. Copyright q

2010 John Wiley & Sons, Ltd.

Int. J. Numer. Meth. Engng (2010) DOI: 10.1002/nme

T. MENOUILLARD AND T. BELYTSCHKO

Figure 14. Geometry and loading of the semi-infinite plate example.

The analytical relation between the stress intensity factor [45] is: ⎧ 0 ⎪ ⎪ ⎪ ⎪  ⎪ ⎪ ⎪ ⎪ 20 c1 (t −tc )(1−2) ⎪ ⎪ ⎨ 1−

˙ t) = K 1 (a, a˙ ⎪  ⎪ 1− ⎪ ⎪ 2 c (t −t )(1−2) cr 0 1 c ⎪ ⎪ ⎪ ⎪ a˙ ⎪ 1−

⎪ ⎩ 1− 2cr

K 1 and the velocity a˙ of the crack if t
where the Rayleigh wave speed is cr = 2947 m/s and the dilatational wave speed is c1 = 5944 m/s. 5.2.1. Regular mesh. We first study the effect of a moving crack on the accuracy of the stress intensity factor with and without the correction force. Two regular meshes are used: 78×39 and 120×59 structured meshes with uniform 4-node elements with one point quadrature and stabilization [52]. Figures 15(a) and (b) present the normalized stress intensity factor as a function of time for the coarse mesh and fine mesh, respectively, with and without the correction force. Both figures show that the correction force improves the result during propagation by decreasing the magnitude of the oscillations. The results are definitely mesh dependent when the crack propagates, because the number of oscillation corresponds to the number of newly cracked elements. The error, as can be seen from Figure 15, is decreased from about 20 to 5% by adding the correction force. Figure 16 shows the stress state (i.e. yy ) in the fine mesh at the end of the computation, for the two cases: with and without the correction force. One can notice the circular waves centered on the crack tip; their number corresponds to the number of newly cracked elements and also corresponds to the number of oscillations observed in the graph showing the stress intensity factors as a function of time (see Figure 15). The correction force makes the stress smoother in the structure when the crack propagation occurs, and significantly attenuates the stress waves released due to the injection of additional degrees of freedom as the crack propagates. 5.2.2. Nonregular mesh. This problem was also run with an unstructured mesh specially constructed so that the orientation does not match the later part of the crack path. Figure 17 shows the nonregular mesh. Figure 18 presents the results obtained with and without the correction force in the unstructured mesh. One can notice that the oscillations are decreased in the computation using the correction force. The results are globally better than the regular mesh because the Copyright q

2010 John Wiley & Sons, Ltd.

Int. J. Numer. Meth. Engng (2010) DOI: 10.1002/nme

IMPROVED DYNAMIC CRACK PROPAGATION MODELING IN XFEM

1.2

1.2

Analytical With correction Without correction

1

0.8 K1 / σ0 √h

0.8 K1 / σ0 √h

Analytical With correction Without correction

1

0.6

0.6

0.4

0.4

0.2

0.2

0

0 0

(a)

0.5

1

1.5 t / tc

2

2.5

3

0

(b)

0.5

1

1.5 t / tc

2

2.5

3

Figure 15. Normalized stress intensity factor as a function of time in the (a) coarse and (b) fine mesh: in both cases, analytical and numerical results with and without the correction force are shown.

(a)

(b)

Figure 16. Stress field in the fine mesh: (a) without the correction force and (b) with the correction force.

Figure 17. Nonregular mesh for mode 1 crack propagation test.

nonregular mesh is finer. More oscillations appear after time 2.5tc because the crack tip reaches a region where the mesh is coarser. However, the correction force actually improves the results in terms of stress intensity factors independent of the mesh. Copyright q

2010 John Wiley & Sons, Ltd.

Int. J. Numer. Meth. Engng (2010) DOI: 10.1002/nme

T. MENOUILLARD AND T. BELYTSCHKO

1.2

Analytical Without correction

1

1

0.8

0.8 K1 / σ0 √h

K1 / σ0 √h

1.2

0.6

Analytical With correction

0.6

0.4

0.4

0.2

0.2

0

0 0

0.5

1

(a)

1.5

2

2.5

3

0

0.5

1

(b)

t / tc

1.5

2

2.5

3

t / tc

Figure 18. Stress intensity factor as a function of time for the unstructured mesh: (a) without the correction force and (b) with the correction force.

1.2

Analytical Without Correction and 1% Damping With Correction and 1% Damping

1

K1 / σ0 √h

0.8 0.6 0.4 0.2 0 0

0.5

1

1.5

2

2.5

3

t / tc

Figure 19. Stress intensity factor as a function of time for the unstructured mesh: analytical solution is compared with the numerical results obtained without and with the correction force.

Figure 19 presents the stress intensity factor as function of time computed with 1% damping. It shows that the damping dramatically decreases the stress intensity factor and results in significant error. Even very small damping does not improve the results compared with the correction force. 5.3. Infinite strip in mode 1 We consider a mode 1 crack propagating as shown in Figure 20. The bottom of the specimen is fixed, whereas the top is constrained by a vertical displacement. The initial crack is horizontal; thus mode 1 crack propagation occurs. The dimensions of the specimen shown in Figure 20 are the following: length L = 2 mm, height h = 0.2 mm and the initial crack length a = 0.09 mm. The imposed displacement U0 is equal to 0.007 mm. The material model is linear elastic. The material Copyright q

2010 John Wiley & Sons, Ltd.

Int. J. Numer. Meth. Engng (2010) DOI: 10.1002/nme

IMPROVED DYNAMIC CRACK PROPAGATION MODELING IN XFEM

Uo

a

h

L

Figure 20. Geometry and loading.

Kinetic Deformation Fracture SUM

1

Energy (Nm)

0.8 0.6 0.4 0.2 0 -0.2 0

0.2

0.4

0.6

0.8

1

1.2

t * cr / L

Figure 21. Energies as a function of time (Kinetic, Deformation and Fracture energies) for the Discretizations: with correction force and without. Both results are identical though.

properties are the following: Young’s modulus E = 3.24 GPa, the Poisson ratio  = 0.35, density  = 1090 kg/m3 . A cohesive law is used for the crack with fracture energy G F = 352.3 Pam. The propagation criterion is the maximum principal strain: εc = 0.02. The specimen is discretized with 90 by 9 4-node elements. Two computations are run: with and without the correction force. As expected, the crack path is strictly straight. Figure 21 presents the different energies for both computations (with and without the correction force); the kinetic, deformation and fracture energies are shown in the Figure 21 as a function of time. Both simulations give similar results in terms of energy, so that only one curve is seen in Figure 21 instead of two. Figure 21 also presents the total energy as a function of time during the crack propagation. Both curves show a constant total energy during the propagation. These results are important to underline the energy conservation in this simulation, whether or not the correction force is used. One concludes that there is no energy dissipation due to the correction force, and the correction force does not perturb the energy conservation or the energy partitioning. 5.4. Three points bending test Here, a beam in three points bending is considered as previously computed by Pandolfi et al. [3]. A precracked plate is supported at two bottom corners. The loading corresponds to a falling projectile which impacts the top center of the specimen at a given velocity; experimentally, a Copyright q

2010 John Wiley & Sons, Ltd.

Int. J. Numer. Meth. Engng (2010) DOI: 10.1002/nme

T. MENOUILLARD AND T. BELYTSCHKO

10m/s

304.8 127

27

Figure 22. Geometry and loading (dimensions are in mm). Table I. Material properties for the three points bending test. Material property Young modulus Poisson ratio Density Initial yield stress Hardening exponent Reference plastic strain Cohesive law Mode 1 static toughness Rayleigh wave velocity

Symbol

Value

E

200 GPa 0.3 7830 kg/m3 1 GPa 10 4.9E−6 15 000 Pa √ m 57.5 MPa m 2906 m/s

  0Y

n ε0p GF K Ic cr

drop tower is used to impose the impact velocity. The geometry and the boundary conditions are summarized in Figure 22. A regular mesh is used at the center region; the corresponding element size is 3 mm×3 mm. Two computations were performed: a standard XFEM without correction force and the same discretization with the correction force. The material is the brittle C300 Steel. The J2 plasticity with power hardening is used to model the material (see [41]). Table I presents the material model properties and the cohesive law parameters. The crack is driven by the mode 1 fracture toughness, and the direction is computed by Equation (12), i.e. the maximum hoop stress criterion. Both results (i.e. with and without the correction force) show a straight vertical crack path as expected according to the boundary conditions. Figure 23 presents the crack lengths as a function of time. The experimental data also appear in the graph. The two numerical computations give very similar results, and both agree quite well with the experiment. Figure 24 shows the evolution of the stress intensity factor as a function of time. One observes that the stress intensity factor is smoother with the correction force than without the correction force. However, the crack length as a function of time is the same with or without the correction force. 5.5. Kalthoff’s experiment This example is based on the experiment of Kalthoff [53] and by B¨ohme and Kalthoff [54] and deals with the crack propagation which initiates in mode 2. A plate with two symmetrical edge cracks is impacted by a projectile at speed V0 . The two cracks are centered with respect to the specimen’s geometry and their separation corresponds to the diameter of the projectile. The crack propagates Copyright q

2010 John Wiley & Sons, Ltd.

Int. J. Numer. Meth. Engng (2010) DOI: 10.1002/nme

IMPROVED DYNAMIC CRACK PROPAGATION MODELING IN XFEM

90

Without correction With correction Experiment

Crack length (mm)

80 70 60 50 40 30 20 0

20

40

60

80

100

120

140

160

180

200

Figure 23. Crack lengths as a function of time, with and without the correction force compared with the experimental results. 90

Without correction

80

80

70

70

Stress intensity factors (MPa√m)

Stress intensity factors (MPa√m)

90

60 50 40 30 20 10

60 50 40 30 20 10

0

0

-10

-10 0

20

40

60

80

100

120

140

(a)

160

180

200

With correction

0

20

40

60

80

100

120

140

160

180

200

(b)

Figure 24. Stress intensity factor as a function of time for the three points bending example: (a) without the correction force and (b) with the correction force.

with an overall angle from 60 to 70◦ . We chose V0 = 20 m/s as a typical velocity for brittle fracture. A schematic description of the problem and the geometry is given by Figure 25. The dimensions of the specimen are given by: L = 100 mm, l = 50 mm, a = 50 mm and the thickness is 16.5 mm. The material properties of the linear elastic media in plane strain are given by: Young’s modulus E = 190 GPa, √ the Poisson ratio  = 0.3 and the density  = 8000 GPa. The fracture toughness is K I c = 68 MPa m and the fracture energy for the cohesive model is 22 146 Pa m. This problem is widely used to validate numerical methods in dynamic crack propagation. For the XFEM, the interest is the fact that whereas the mesh is structured, the crack propagates with an angle of 60–70◦; so the crack path does not respect the orientation of the mesh, i.e. it is not parallel to any element edges. A regular mesh (38×38) with 4-node elements was employed in the two simulations (with and without the correction force). Copyright q

2010 John Wiley & Sons, Ltd.

Int. J. Numer. Meth. Engng (2010) DOI: 10.1002/nme

T. MENOUILLARD AND T. BELYTSCHKO

L

V0 2L

l

a

Figure 25. Geometry of Kalthoff’s experiment. With correction Without correction

Stress intensity factors (MPa√m)

140 120 100 80 60 40 20 0 30

40

50

60

70

80

90

Figure 26. Mode 1 stress intensity factors as a function of time for both computations: with and without the correction force.

Figure 26 presents the mode 1 stress intensity factors as a function of time for both computations (i.e. with and without the correction force). It shows that the computation using the correction force gives less oscillations in the stress intensity factor. Figure 27 shows the crack paths for both computations. One observes that both paths are similar. Figure 28 presents the crack lengths as a function of time. Both computations give similar results and they are almost identical in term of crack length too. Figure 29 presents snapshots of the stress xx for the two different computations; (a) corresponds to the computation without any correction force, whereas (b) corresponds to the one with the correction. The stress states are quite different in the two simulations. One observes that the stress obtained without the correction (Figure 29(a)) presents spurious waves centered on the current crack tip, whereas with correction the stress is smoother as shown in Figure 29(b). The Copyright q

2010 John Wiley & Sons, Ltd.

Int. J. Numer. Meth. Engng (2010) DOI: 10.1002/nme

IMPROVED DYNAMIC CRACK PROPAGATION MODELING IN XFEM

With correction Without correction

70˚

Figure 27. Crack paths: with and without the correction force. 90 80

Crack length (mm)

70 60 50 40 30 20 10

With correction Without correction

0 0

10

20

30

40

50

60

70

80

90

Figure 28. Crack lengths as a function of time for both computations: with and without the correction force.

stress waves observed in the results without the correction force correspond to the sudden release of the newly cracked elements during the propagation. These stress oscillations are numerical, and do not match any physical aspect of crack propagation, and they are mesh dependent.

6. CONCLUSION A method to smooth the release of the crack tip element with only the discontinuous enrichment in XFEM has been proposed. This method is based on enforcing the continuity of the forces Copyright q

2010 John Wiley & Sons, Ltd.

Int. J. Numer. Meth. Engng (2010) DOI: 10.1002/nme

T. MENOUILLARD AND T. BELYTSCHKO

Figure 29. Stress xx in the specimen for the two computations: (a) without the correction force and (b) with the correction force.

corresponding to the newly added enriched degrees of freedom, when the crack tip just reaches a new element; indeed it makes the acceleration corresponding to the enriched degrees of freedom continuous, and thus velocity and displacement once and twice continuously differentiable. This method is not restricted to a straight crack path, but applies to any direction of propagation in the mesh. The results of the Kalthoff experiment have shown the improvement in term of stress intensity factors and stress state by using this correction, even if the crack path does not follow the orientation of the mesh. The control of the acceleration is performed through the use of a correction force linked to the virtual position of the crack tip in the element. The restriction is that the crack tip velocity has to be known or estimated for the time required for the crack to cut an element. It has been shown that the release stress wave of the newly cracked element is significantly decreased by using the proposed correction. This has been shown through the computation of the stress intensity factors during propagation, and by observing the stress state in the structure. The use of this correction force improves the accuracy of the stress intensity factors during propagation and decreases its mesh dependence. It also smoothes the corresponding derivatives of the displacement, i.e. crack velocity and stress, because it decreases the spurious numerical stress due to the release of the mesh. However, it does not appear to change the global results (i.e. crack path and crack length) of the simulation significantly. Thus its relevance lies mainly in inverse problems where dynamic crack propagation is inferred from velocity or strain time histories and in problems where accurate stress intensity factors are desirable.

ACKNOWLEDGEMENTS

The support of the Office of Naval Research under Grants N00014-08-1-1191 and N00014-09-1-0030 is gratefully acknowledged. Copyright q

2010 John Wiley & Sons, Ltd.

Int. J. Numer. Meth. Engng (2010) DOI: 10.1002/nme

IMPROVED DYNAMIC CRACK PROPAGATION MODELING IN XFEM

REFERENCES 1. Pandolfi A, Ortiz M. Solid modeling aspects of three-dimensional fragmentation. Engineering with Computers 1998; 14(4):287–308. 2. Ortiz M, Pandolfi A. Finite-deformation irreversible cohesive elements for three-dimensional crack-propagation analysis. International Journal for Numerical Methods in Engineering 1999; 44(9):1267–1282. 3. Pandolfi A, Guduru PR, Ortiz M, Rosakis AJ. Three dimensional cohesive-element analysis and experiments of dynamic fracture in C300 steel. International Journal of Solids and Structures 2000; 37(27):3733–3760. 4. Xu XP, Needleman A. Numerical simulations of dynamic crack growth along an interface. International Journal of Fracture 1995; 74(4):289–324. 5. Linder C, Armero F. Finite elements with embedded strong discontinuities for the modeling of failure in solids. International Journal for Numerical Methods in Engineering 2007; 72:1391–1433. 6. Hallquist JO. LS-DYNA Theory Manual. Livermore software Technology corporation, 2006. 7. Belytschko T, Tabbara M. Dynamic fracture using element-free Galerkin methods. International Journal for Numerical Methods in Engineering 1996; 39(6):923–938. 8. Belytschko T, Krongauz Y, Organ D, Fleming M, Krysl P. Meshless methods: an overview and recent developments. Computer Methods in Applied Mechanics and Engineering 1996; 139(1–4):3–47. 9. Belytschko T, Krongauz Y, Dolbow J, Gerlach C. On the completeness of meshfree particle methods. International Journal for Numerical Methods in Engineering 1998; 43:785–819. 10. Fleming M, Chu YA, Moran B, Belytschko T. Enriched element-free Galerkin methods for crack tip fields. International Journal for Numerical Methods in Engineering 1997; 40(8):1483–1504. 11. Bordas S, Rabczuk T, Zi G. Three-dimensional crack initiation, propagation, branching and junction in non-linear materials by an extended meshfree method without asymptotic enrichment. Engineering Fracture Mechanics 2008; 75(5):943–960. 12. Swenson DV, Ingraffea AR. Modeling mixed-mode dynamic crack propagation using finite elements: theory and applications. Computational Mechanics 1988; 3(6):381–397. 13. Black T, Belytschko T. Elastic crack growth in finite elements with minimal remeshing. International Journal for Numerical Methods in Engineering 1999; 45:601–620. 14. Mo¨es N, Dolbow J, Belytschko T. A finite element method for crack growth without remeshing. International Journal for Numerical Methods in Engineering 1999; 46(1):131–150. 15. Belytschko T, Mo¨es N, Usui S, Parimi C. Arbitrary discontinuities in finite elements. International Journal for Numerical Methods in Engineering 2001; 50(4):993–1013. 16. Stolarska M, Chopp DL, Mo¨es N, Belytschko T. Modelling crack growth by level sets in the extended finite element method. International Journal for Numerical Methods in Engineering 2001; 51:943–960. 17. Wells GN, Sluys LJ. A new method for modelling cohesive cracks using finite elements. International Journal for Numerical Methods in Engineering 2001; 50(12):2667–2682. 18. Belytschko T, Chen H, Xu J, Zi G. Dynamic crack propagation based on loss of hyperbolicity and a new discontinuous enrichment. International Journal for Numerical Methods in Engineering 2003; 17(5):679–706. 19. Remmers JJC, de Borst R, Needleman A. A cohesive segments method for the simulation of crack growth. Computational Mechanics 2003; 31(1):69–77. 20. Mo¨es N, Belytschko T. Extended finite element method for cohesive crack growth. Engineering Fracture Mechanics 2002; 69:813–833. 21. de Borst R. Numerical aspects of cohesive-zone models. Engineering Fracture Mechanics 2003; 70(14): 1743–1757. 22. Zi G, Belytschko T. New crack-tip elements for XFEM and applications to cohesive cracks. International Journal for Numerical Methods in Engineering 2003; 57(15):2221–2240. 23. Rabczuk T, Zi G, Gerstenberger A, Wall WA. A new crack tip element for the phantom-node method with arbitrary cohesive cracks. International Journal for Numerical Methods in Engineering 2008; 75(5):577–599. 24. Song JH, Areias PMA, Belytschko T. A method for dynamic crack and shear band propagation with phantom nodes. International Journal for Numerical Methods in Engineering 2006; 67:868–893. 25. R´ethor´e J, Gravouil A, Combescure A. An energy-conserving scheme for dynamic crack growth using the extended finite element method. International Journal for Numerical Methods in Engineering 2005; 63:631–659. 26. R´ethor´e J, Gravouil A, Combescure A. A combined space–time extended finite element method. International Journal for Numerical Methods in Engineering 2005; 64:260–284. 27. Menouillard T, Belytschko T. Correction force for releasing crack tip element with xfem and only discontinuous enrichment. European Journal of Computational Mechanics 2009; 18(5/6):465–483. Copyright q

2010 John Wiley & Sons, Ltd.

Int. J. Numer. Meth. Engng (2010) DOI: 10.1002/nme

T. MENOUILLARD AND T. BELYTSCHKO

28. Belytschko T, Chen H. Singular enrichment finite element method for elastodynamic crack propagation. International Journal of Computational Methods 2004; 1(1):1–15. 29. Menouillard T, Song JH, Duan Q, Belytschko T. Time dependent crack tip enrichment for dynamic crack propagation. International Journal of Fracture 2009; DOI: 10.1007/s10704-009-9405-9. 30. Menouillard T, Belytschko T. Dynamic fracture with meshfree enriched XFEM. ACTA Mechanica 2009; DOI: 10.1007/s00707-009-0275-z. 31. Aberson JA, Anderson JM. Cracked finite elements proposed for NASTRAN. Third NASTRAN Users’ Colloquim. NASA TMX-2893, NASA Langley Research Center, Langley, VI, 1973; 531–550. 32. Hardy RH. A high order finite element for two dimensional crack problems. Ph.D. Thesis, Georgia Institute of Technology, August 1974. 33. King WW, Malluck F. Toward a singular element for propagating cracks. International Journal of Fracture 1978; 14(1):7–11. 34. Malluck JF, King WW. Fast fracture simulated by the conventional finite elements: a comparison of two energy-release algorithms. Crack Arrest Methodology and Applications 1980; 711:38–53. 35. Kobayashi AS, Emery AF, Mall S. Dynamic-finite-element and dynamic-photoelastic analyses of two fracturing homalite-100 plates. Experimental Mechanics 1976; 16(9):321–328. 36. Thesken JC, Gudmundson P. Application of a moving variable order singular element to dynamic fracture mechanics. International Journal of Fracture 1991; 52(1):47–65. 37. Rydholm G, Frederiksson B, Nilsson F. Numerical investigations of rapid crack propagation. In Numerical Methods in Fracture Mechanics; Proceedings of the First International Conference, Luxmoore AR, Owen DRJ (eds), Department of Civil Engineering, University College of Swansea, Swansea, January 1978; 660–672. 38. Owen DRJ, Hinton E. Finite Elements in Plasticity: Theory and Practice. Pineridge Press: Swansea, 1980. 39. Brickstad B, Nilsson F. Numerical evaluation by FEM of crack propagation experiments. International Journal of Fracture 1980; 16(1):71–84. 40. Fan R, Fish J. The rs-method for material failure simulations. International Journal for Numerical Methods in Engineering 2008; 73(11):1607–1623. 41. Belytschko T, Liu WK, Moran B. Nonlinear Finite Elements for Continua and Structures. Wiley: Chichester, 2000. 42. Barenblatt GI. The mathematical theory of equilibrium cracks formed in brittle fracture. Advanced Applied Mechanics 1962; 7:55–129. 43. Hillerborg A, Modeer M, Petersson PE et al. Analysis of crack formation and crack growth in concrete by means of fracture mechanics and finite elements. Cement and Concrete Research 1976; 6(6):773–782. 44. Rabczuk T, Song JH, Belytschko T. Simulations of instability in dynamic fracture by the cracking particles method. Engineering Fracture Mechanics 2009; 76:730–741. 45. Freund LB. Dynamic Fracture Mechanics. Cambridge University Press: Cambridge, 1990. 46. Freund LB. Crack propagation in an elastic solid subjected to general loading. Pt. 1. Constant rate of extension. Journal of the Mechanics and Physics of Solids 1972; 20(3):129–140. 47. Freund LB, Douglas AS. Influence of inertia on elastic–plastic antiplane-shear crack growth. Journal of the Mechanics and Physics of Solids 1982; 30(1):59–74. 48. Rosakis AJ, Freund LB. Optical measurement of the plastic strain concentration at a crack tip in a ductile steel plate. Journal of Engineering Materials and Technology (Transactions of the ASME) 1982; 104(2):115–120. 49. Menouillard T, R´ethor´e J, Combescure A, Bung H. Efficient explicit time stepping for the eXtended finite element method (X-FEM). International Journal for Numerical Methods in Engineering 2006; 68:911–939. 50. Menouillard T, R´ethor´e J, Mo¨es N, Combescure A, Bung H. Mass lumping strategies for X-FEM explicit dynamics: application to crack propagation. International Journal for Numerical Methods in Engineering 2008; 74(3):447–474. 51. Ravi-Chandar K. Dynamic Fracture. Elsevier Science: Amsterdam, 2004. 52. Flanagan DP, Belytschko T. A uniform strain hexahedron and quadrilateral with orthogonal hourglass control. International Journal for Numerical Methods in Engineering 1981; 58(12):1873–1905. 53. Kalthoff JF. On the measurement of dynamic fracture toughnesses—a review of recent work. International Journal of Fracture 1985; 27(3):277–298. 54. B¨ohme W, Kalthoff JF. The behavior of notched bend specimens in impact testing. International Journal of Fracture 1982; 20(4):139–143.

Copyright q

2010 John Wiley & Sons, Ltd.

Int. J. Numer. Meth. Engng (2010) DOI: 10.1002/nme