Thomas Menouillard · Ted Belytschko

Dynamic fracture with meshfree enriched XFEM

Received: 1 December 2009 / Published online: 21 January 2010 © Springer-Verlag 2010

Abstract The enrichment of the extended finite element method (XFEM) by meshfree approximations is studied. The XFEM allows for modeling arbitrary discontinuities, but with low order elements the accuracy often needs improvement. Here, the meshfree approximation is used as an enrichment in a cluster of nodes about the crack tip to improve accuracy. Several numerical examples show that this leads to more accuracy for stress intensity factors computations, and to the capability to capture the branching point of a propagating crack from the stresses.

1 Introduction Finite element methods have been widely used for computational mechanics because of their efficiency. However, they have some difficulties to simulate crack propagation mainly due to the use of a mesh. Crack propagation requires an update of the mesh during the computation, and thus projections of various fields. Whereas it is possible to use remeshing in 2D, it can become a real challenge in 3D, particularly for dynamic crack propagation. The simplest way is to use element deletion, but it brings a significant mesh dependence. Some techniques can avoid using a mesh and are particularly adapted to crack propagation, i.e., meshfree methods [1–4]. A second way is the extended finite element method (XFEM) [5,6]. An enriched basis of shape function is added into the standard finite element discretization via the partition of unity [7], and this allows to enhance the description of the physics (i.e., crack) independently of the mesh [8,9]. The crack propagation phenomenon is modeled through an enrichment function; this is discontinuous to characterize the crack discontinuity. The description of the crack tip vicinity can be also enhanced by a tip enrichment. The main advantage of XFEM for crack propagation remains in the possibility to compute 3D crack propagation, where finite elements show their limits due to the remeshing process. For 3D problems, the coupling of XFEM with the level set method to deal with the tracking process of the crack independently of the mesh [10,11] is promising. Dynamic crack propagation by XFEM has been developed in [12–14]. The first two restricted the method to element-by-element propagation of the crack. In [14], a special element was developed for dynamics with partially cracked element. Belytschko and Chen [15] studied dynamic fracture with singular enrichments. One difficulty with XFEM is that it sometimes does not have the desired accuracy. Recent work on the time-dependent tip enrichment [16] and on the release of the crack tip element [17,18] showed some improvement of the accuracy of XFEM in the stress intensity factor computation. An alternative way to The authors gratefully acknowledge the support of the Army Research Office through Grant W911NF-08-1-0212 and the Office of Naval Research through Grant N00014-08-1-1191. T. Menouillard (B) · T. Belytschko Theoretical and Applied Mechanics, Northwestern University, Evanston, IL 60208-3111, USA Tel.: +1-847-4914029 E-mail: [email protected]; [email protected]

54

T. Menouillard, T. Belytschko

enhance the quality of the stress field near the crack tip is adaptivity driven by the error estimation developed in [19–22]. In this paper we study whether the enrichment of an XFEM approximation by a meshless approximation can improve the accuracy of stresses in the vicinity of the crack tip and stress intensity factor computations. In particular, we consider cohesive cracks, where nearfield asymptotic solutions are not known for dynamics. We also examine whether with improved stresses the branching of cracks can be driven by the stress criteria. The enrichment of finite element approximations by meshfree approximations was reported by Liu et al. [23] and Huerta and Fernandez-Mendez [24]. Convergence studies are given in [25]. Coupled finite element/meshless methods have been reported in [26,27]. To our knowledge, this is the first paper which studies the XFEM enriched with meshfree approximations. The main advantage of the meshfree enrichment is that nodes can easily be added in the vicinity of the crack tip and then removed when they are not needed. In contrast, finite element adaptivity requires rather complex data structures. We enrich XFEM with the meshfree method near the crack tip to enhance the stress state and to capture a better description to drive the crack or to catch physical phenomena such as crack branching. The paper is organized as follows. Section 2 presents the variational equations. Section 3 introduces and develops the space discretization, i.e., the coupling of the XFEM and the meshfree method. Section 4 presents the time integration and the discrete equations. Section 5 explains the fracture mechanics and crack law. Section 6 describes the numerical examples underlining the effect of the coupling in the different results. Section 7 presents the remarks and conclusion.

2 Variational equations The method is applicable to two and three dimensions, but only the 2D case is described. The formulation of the 3D case is identical, but the limitation comes from computer memory to run simulations. Consider an initial domain Ω 0 and its boundary ∂Ω F0 and ∂Ωu0 as shown in Fig. 1 with ∂Ωu0 ∩ ∂Ω F0 = ∅. 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 describes the reference and current domains, and the different boundaries with their corresponding conditions. The initial crack is denoted by Γc0 and its image in the current domain Ω is Γc . The momentum equation is written as: ρ u¨ = div σ + fd , (1) where ρ is the density, 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 (see e.g., [28]): ρ u¨ · v dΩ + σ : ε(v) dΩ = fd · v dΩ + Fd · v dΓ + Fcoh · v dΓ ∀ v ∈ U0 , Ω/Γc

Ω/Γc

Ω/Γc

∂Ω F

Γc

(2)

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

Dynamic fracture with meshfree enriched XFEM

55

where Fd is the cohesive force along the crack, ε the strain tensor, and v ∈ U0 . The different domains U and U0 are defined by: U = {u | u (x) = ud , ∀ x ∈ ∂Ωu , and regularity of u} , U0 = {u | u (x) = 0, ∀ x ∈ ∂Ωu , and regularity of u} .

(3) (4)

The boundary of the current domain Ω is partitioned into ∂Ωu on which displacements ud are prescribed, ∂Ω F on which tractions Fd are prescribed, and then Γc which corresponds to the displacement discontinuity, i.e., crack. One can sum up the boundary conditions by: u(x, t) = ud , σ (x, t) · n = Fd , σ (x, t) · n = Fcoh ,

∀ x ∈ ∂Ωu , ∀ x ∈ ∂Ω F , ∀ x ∈ Γc ,

(5)

where n denotes the outer normal of the boundary ∂Ω, t is the time. Moreover, the latter boundary Γc also depends on time as the crack propagates: Γc (t). Note that ∂Ω = ∂Ωu ∪ ∂Ω F ∪ Γc , ∂Ωu ∩ Γc = ∅ and ∂Ω F ∩ Γc = ∅. 3 Space discretization: coupling XFEM and meshfree The coupling between the two discretizations (i.e., XFEM and meshfree) arises because the strain is common ˜ to both, so the stress is too. The discretization of the displacement field shows that the total displacement u(x, t) is the sum of the XFEM discretized displacement denoted by u˜ XFEM (x, t) plus the meshfree part denoted by u˜ MF (x, t), as follows: ˜ u(x, t) = u˜ XFEM (x, t) + u˜ MF (x, t).

(6)

Thus, the deformation gradient is the sum of both deformation gradients, and the two algorithms can run independently. 3.1 XFEM discretization In XFEM, the enrichments are added to the conventional finite element approximation according to the partition of unity method. These functions, denoted by H and φ in this paper, take into account the physical phenomena in their formulation. Different functions are used to enhance different parts of the structure (i.e., along the crack and the crack tip vicinity). A discontinuous function H is typically used to describe a displacement discontinuity currently characterized by the crack; f is the level set function whose isozero defines the crack surface. Thus, the discretized displacement u˜ XFEM is written as: N I (x) u I + H ( f (x)) N J (x) b J + φ(x)N K (x) e K , (7) u˜ XFEM (x, t) = I ∈N

J ∈Ncut

K ∈Ntip

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, and Ntip the set of nodes belonging to the crack tip element or zone (e.g., 3 × 3 elements near crack tip). The standard shape functions are denoted by N I , and H denotes 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 tip enrichment function, denoted by φ, is related to the crack tip position and the direction of the crack, and enhances the behavior of the structure at the vicinity of the crack tip. The standard degrees of freedom are denoted by u I ; b J correspond to the enriched degrees of freedom due to the discontinuous function H , and e K correspond to the tip enrichment degrees of freedom. Figure 2 presents a mesh, a crack, and the corresponding sets of enriched nodes. For linear elastic fracture mechanics, the most commonly used tip enrichment is defined by four functions with a square root singularity, which builds a basis allowing to represent the displacement near the crack tip. The discontinuity ahead of the crack tip is represented by the expression sin θ2 . One has to notice that this popular basis is similar to the asymptotic displacement field near the crack tip, but is not exactly the same. Nevertheless, the asymptotic displacement field of modes 1 and 2 can be generated by the function basis.

56

T. Menouillard, T. Belytschko

Fig. 2 Different sets of the nodes along the crack: N , Ncut defined by squares, and Ntip defined by circles

Recently, Cox [29] used an analytical displacement field near the crack tip based on a static solution for the cohesive stresses to enrich the shape functions. Our developments aim at describing more accurately the displacement near the crack tip. We chose to use φ as tip enrichment: φ ∈ {r sin(θ/2); r 2 sin(θ/2)}.

(8)

This enrichment was used by Moës and Belytschko [9] for static crack propagation. It underlines the behavior in front of the crack tip with a larger displacement gradient [29], and the discontinuity in the term sin(θ/2). The motivation to use such enrichments [16] is that the discontinuity along the crack is taken into account by the Heaviside function, and the displacement behavior in front of this discontinuity is taken into account by the tip enrichment with a higher order function. Only the current crack tip zone is enriched with the tip enrichment. When the crack tip reaches a new element, the enrichment type of some nodes of this zone changes, from φ to H , and the degree of freedom corresponding to the Heaviside function has to be initiated. With this strategy, the initial value of the new discontinuous enrichment will not be taken as zero as it is effected [30]. The main point is that the sum on the old crack tip enrichment is not kept during the propagation. Menouillard et al. [16] used this same strategy in their simulations of dynamic crack propagation. Whereas the energy conservation is not exact, the results in term of crack propagation (i.e., stress intensity factors, crack velocity and length) are similar to the strategy strictly conserving the energy.

3.2 Meshfree discretization An additional term u˜ MF (x, t) of meshfree enrichment is added to the XFEM discretization as shown in Eq. (6). This term is written as follows: u˜ MF (x, t) =

Ψ L (x) d L (t),

(9)

L∈Nmf

where Nmf is the set of meshfree nodes, Ψ L (x) the meshfree shape function at node L, and d L (t) the corresponding degree of freedom of the meshfree node L. The shape functions Ψ L evaluated at point XG are given by: Ψ L (XG ) = w L (XG ) P(X L )T · α(XG ),

(10)

where P(X) is the basis functions vector, usually [1, x, y] for linear interpolation, [1, x, y, x2 , xy, y 2 ] for quadratic interpolation. Here, the linear basis is described by the linear finite elements shape functions, while the meshfree method carries only the quadratic basis, i.e.,[x2 , y 2 ]. L| where r0 is the radius of the influence Let us introduce the kernel function w L (XG ) = w |XGr−X 0 domain in meshfree methods. We use the following kernel function:

Dynamic fracture with meshfree enriched XFEM

57

⎧2 2 3 ⎪ ⎨ 3 − 4s + 4s w(s) = 43 − 4s + 4s 2 − 43 s 3 ⎪ ⎩ 0

s < 21 , 1 2

≤ s < 1,

(11)

1 ≤ s.

The matrix α has the same dimension as the vector P, and has to be determined so that the partition of unity is verified; so for all points XG , it must verify

Ψ I (XG ) P(X I ) = P(XG ), ∀ XG ∈ Ω.

(12)

I ∈Nmf

By injecting Eq. (10) into Eq. (12), the resolution of α becomes the solution of a linear system whose dimension is the same as the matrix P : ⎛ ⎝

⎞ w I (XG ) P(X I )T · P(X I )⎠ · α(XG ) = A · α(XG ) = P(XG ), ∀ XG ∈ Ω,

(13)

I ∈Nmf

where the matrix A is defined by the weight function w I and the matrix P at all points. Once α is known, the shape function Ψ I can be determined at all points XG by injecting the solution into Eq. (10). Figure 3 presents a summary of the discretization coupling XFEM and meshfree; it shows the FEM mesh and the crack, and also the enriched nodes of the XFEM formulation, the meshfree nodes and the quadrature points. As the quadrature points remain the same for both discretizations, it makes the stress and strain update between both quite easy. Note that the meshfree approximation functions are continuous. The discontinuity in the XFEM approximation suffices to capture the physics of the crack.

Fig. 3 Space discretization: meshfree and XFEM formulation

58

T. Menouillard, T. Belytschko

4 Time integration and discrete equations 4.1 Time integration The central difference method [28] is used: ˙ t + 1 Δt U ˙ t+Δt/2 = U ¨ t, U 2 ˙ t+Δt/2 , Ut+Δt = Ut + Δt U ¨ t+Δt = M ·U

ext Ft+Δt

int − Ft+Δt

(14) (15)

coh + Ft+Δt ,

(16)

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

(17)

˙ t , and U ¨ t ) denotes the displacement (respectively, velocity, and acceleration) vector where Ut (respectively U as degrees of freedom at time t indicated at the bottom right of the variables. Δt is the time step size, M the mass matrix, and Fext (respectively, Fint and Fcoh ) the external forces (respectively, internal and cohesive forces) as a vector of degrees of freedom. Equation (16) is described in the next paragraph. The time step Δt defines the time between the current time t and the next t + Δt; it can vary with time instead of being constant. A diagonal mass matrix is used. The advantage is to speed up the computation, and use less memory by storing only vectors instead of matrices. The time step of the computation is given by the critical time step of the meshfree part because it is smaller than the critical time step of the corresponding XFEM element. We used a time step of 0.5 h M/c1 , where h M is the distance between meshfree particles. 4.2 Discrete equations The discrete equations are obtained via the weak form, i.e., Eq. (2). The test function v is taken to be N I (x) u I (t)+ H ( f (x)) N J (x) b J (t)+ φ(x) N K (x) e K (t) + Ψ L (x) d L (t). v(x, t) = I ∈N

J ∈Ncut

K ∈Ntip

L∈Nmf

(18) Substituting the trial function into the weak form (2) gives the following equation: ∀ I ∈ N ∪ Nmf int coh ¨ J = Fext MI J U I − FI + FI ,

and the remaining terms are defined below: ⎧ u,int ⎫ ⎧ u,ext ⎫ ⎧ ⎫ ⎧ ⎫ 0 ⎪ FI ⎪ FI ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ u¨ J ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨¨ ⎪ ⎨ Fb,int ⎪ ⎨ Fb,ext ⎪ ⎨ Fb,coh ⎪ ⎬ ⎬ ⎬ ⎬ b I I I J int ext coh ¨J = , FI = , F , F , = = U I I e,coh e,ext e,int ⎪ e¨ J ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ F F F ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ I I ⎩ d¨ ⎭ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ I ⎪ ⎪ ⎪ ⎩ d,int ⎪ ⎩ d,ext ⎪ ⎩ d,coh ⎪ ⎭ ⎭ ⎭ J FI FI FI where the internal forces defined on different degrees of freedom are: u,int FI i = N I, j σ ji dΩ ∀ I ∈ N ,

(19)

(20)

(21)

Ω

Fb,int Ii

=

N I, j H ( f (x)) σ ji dΩ ∀ I ∈ N cut ,

(22)

(N I φ(x)), j σ ji dΩ ∀ I ∈ N tip ,

(23)

Ψ I, j σ ji dΩ ∀ I ∈ Nmf .

(24)

Ω

= Fe,int Ii Ω

Fd,int Ii

= Ω

Dynamic fracture with meshfree enriched XFEM

59

The external forces defined on the different degrees of freedom are: u,ext = N I fd dΩ + N I Fd dΓ ∀ I ∈ N , FI Ω

∂Ω F

Fb,ext I

=

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

N I H ( f (x)) Fd dΓ ∀ I ∈ N cut ,

(26)

∂Ω F

= Fe,ext I

(25)

N I φ(x) fd dΩ ∀ I ∈ N tip ,

(27)

Ψ I fd dΩ ∀ I ∈ Nmf .

(28)

Ω

Fd,int Ii

= Ω

The cohesive forces are defined by: Fb,coh = I

N I n τ dΓ ∀ I ∈ N cut ,

(29)

N I φ( f (x)) n τ dΓ ∀ I ∈ N tip ,

(30)

Γc

= Fe,coh I Γc

where τ denotes the opening tensile stress along the crack. N I, j denotes the derivative of N I by x j . The consistent mass matrix is given by uu (31) M I J = ρ N I N J dΩ, Ω

Mbb IJ

=

ρ N I N J H 2 ( f (x)) dΩ,

(32)

ρ N I N J H ( f (x)) dΩ,

(33)

ρ N I N J φ 2 (x) dΩ,

(34)

ρ N I N J φ(x) dΩ,

(35)

ρ N I N J H ( f (x)) φ(x) dΩ.

(36)

Ω

Mub IJ = Ω

Mee IJ

= Ω

Mue IJ

= Ω

Mbe IJ = Ω

However, the mass matrix was diagonalized as described in [31–33]. Hence, the lumped mass matrix becomes u M I I = ρ N I dΩ, (37) Ω

MbI I

=

ρ N I H 2 ( f (x)) dΩ, Ω

MeI I =

e

ne 2 J =1 φ (xe J )

(38)

n e

ρ φ 2 (x) dΩe , Ωe

(39)

60

T. Menouillard, T. Belytschko

where xe J is the coordinate of node J in the element e. The element e owns n e nodes. The consistent mass matrix related to the meshfree part is written as: 2 Mdd = ρ Ψ I Ψ J dΩ ∀ (I, J ) ∈ Nmf . (40) IJ Ω

The corresponding lumped mass matrix is: Mdd II

=

ρ Ψ I dΩ ∀ I ∈ Nmf .

(41)

Ω

We assumed there are no external forces on the meshfree part, because it is always inside the structure. 5 Fracture mechanics and crack law Fracture mechanics deals with the physical law which drives the crack propagation. First, considering a crack propagation from an energetic point of view, energy is dissipated within the propagation. One way to model this is a cohesive law. Second, the crack position 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. Thus, the complete updating process of the crack tip position requires the direction and norm of the crack velocity. Moreover, the stress intensity factors give the velocity of the crack according to the framework of the linear fracture mechanics. 5.1 Cohesive law The energy dissipation in this fracture zone needs to be taken into account. Since no singularity was included at the crack tip, the energy dissipation is accomplished by a cohesive traction. A good way is to use cohesive-zone models, which were introduced first by Barenblatt [34] for elastic plastic fracture in ductile metals, and by Hillerborg et al. [35] for brittle materials. Figure 4a shows the crack opening and the tensile traction along the crack. It shows the closing effect near the crack tip with a limit opening δmax ; beyond this opening, there is no more tensile stress. The link between the tensile stress and the crack opening along the crack lips defines the cohesive law. Figure 4b presents a linear cohesive law [36,37]; the area under the curve represents the fracture energy G F , defined by: δmax 1 τ (δ) dδ = τmax δmax , GF = 2

(42)

0

where τ denotes the tensile stress along the crack, and τmax , δmax define the linear cohesive law described in Fig. 4b. (a)

(b)

Fig. 4 Cohesive law: a crack opening δ as a function of the position x along the crack, b linear cohesive law where the area in gray under the curve represents the fracture energy G F

Dynamic fracture with meshfree enriched XFEM

61

5.2 Stress intensity factors The main step consists in relating the propagation of the crack to the state of the structure, to evaluate the direction and velocity of the crack as functions of time: denoted, respectively, by θc (t) and a(t). ˙ 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 whole history of the crack geometry in the structure. Rice [38] defined the path independent integral which leads to the energy release rate. The dynamic stress intensity factors K 1 and K 2 are defined by the asymptotic behavior of the stress near the crack tip as [39]: √ K 1 = lim 2πr σyy , (43) r →0 √ (44) K 2 = lim 2πr σxy . r →0

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 [40,30]. 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 σ div(q)dΩ + σ aux : ∇u∇q + σ : ∇uaux ∇q dΩ : ∇u − ρ u˙ u˙ I =− Ω

Ω

¨ aux dΩ + div σ aux ∇u(q) + ρ u∇u

+ Ω

˙ (q)dΩ, ˙ ˙ uaux ρ u˙ aux ∇ u(q) + ρ u∇

(45)

Ω

where σ aux , uaux are auxiliary stress and displacement fields. The vector q, parallel to the crack, is defined by: ⎧ ⎪ ⎨0 q = q = 1 ⎪ ⎩q linear

outside the domain S1 ∪ S2 , inside the domain S1 , inside the domain S2 ,

(46)

where S1 and S2 are domains defined near the crack tip. Figure 5 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. The sides of the domain S1 are about six element sizes, and the lengths of the domain S2 are twice as large. The domain S1 is chosen to be sufficiently large so that the crack surface outside this domain has no cohesive tractions. Then, this integral is: I int =

2(1 − ν 2 ) ˙ 1 K 1aux + β2 (a)K ˙ 2 K 2aux , β1 (a)K E

Fig. 5 Direction and norm of the virtual extension field q

(47)

62

T. Menouillard, T. Belytschko

where βi are the universal functions defined by:

4αi 1 − α22 ˙ = i ∈ {1, 2}, βi (a) (κ + 1)D(a) ˙

(48)

where the Kolosov coefficient κ is 3 − 4ν for the plane strain, or (3 − ν)/(1 + ν) for the plane stress, and D, a 2 ˙ = 4α1 α2 − 1 + α22 . function whose zero defines the Rayleigh wave celerity denoted by cr , is written as D(a) 2 The parameters αi are defined by αi2 = 1 − ca˙i with i ∈ {1, 2} where c1 and c2 are, respectively, the dilatational which, as a functions of the Lamé coefficients and the density are given by: and shear velocities λ+2µ µ and c2 = ρ . c1 = ρ The different stress intensity factors are estimated through an appropriate choice of uaux , i.e., K 1aux = 1 and K 2aux = 0 for the determination of K 1 , and K 1aux = 0 and K 2aux = 1 for the determination of K 2 . The asymptotic displacement fields are actually used to perform the mixed mode decomposition. 5.3 Crack velocity law From the knowledge of the dynamic stress intensity factors, the equivalent dynamic stress intensity factor K θ θ is defined by: 3 θc 3 θc K θ θ = cos (49) K 1 − cos sin θc K 2 . 2 2 2 Freund [41] developed a relation between the dynamic energy release rate and the crack velocity. Freund and Douglas [42] and Rosakis and Freund [43] have linked experimentally the crack velocity to the stress intensity factors through the relation explained by Freund [41]: ⎧ if K θ θ < K 1c , ⎨0 2 a˙ = (50) K 1c otherwise. ⎩ cr 1 − K θ θ Although this law was not derived for cohesive laws, we assume they apply when the dissipation due to the 2 2 cohesive law corresponds to the critical energy dissipation 1−ν E K 1c . In addition, the direction of the crack is driven by the maximum hoop stress criterion as: ⎡ ⎛ ⎞⎤ 2 1 K1 K 1 ⎠⎦ θc = 2 artan ⎣ ⎝ . (51) − sign(K 2 ) 8 + 4 K2 K2 So a crack increment in the explicit algorithm is defined by Δa = a˙ Δt, where Δt is the time step size. 6 Numerical examples 6.1 Mode 1 stationary and moving crack The example considered in this Section is a semi infinite plate with a crack [44] loaded as shown in Fig. 6. A theoretical solution of this problem for a given crack velocity is given in Freund [40]. On this problem no cohesive force is used, and we are strictly interested in the ability of the meshfree approximations to improve stress intensity factors. 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 Fig. 6, those could be compared for time t ≤ 3tc = 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 valid. The dimensions of the structure are the following: the length is L = 10 m, the initial crack length a = 5 m, and the vertical position of the crack is h = 2 m. A regular mesh of 78 × 39 four-node elements is used. The material properties of the linear elastic media are: Young’s modulus E = 210 GPa, Poisson’s ratio ν = 0.3, and the density

Dynamic fracture with meshfree enriched XFEM

63

Fig. 6 Geometry and loading of the semi-infinite plate example

ρ = 8,000 kg/m3 . The tensile stress applied on the top surface is σ0 = 500 MPa. The mode 1 stress intensity √ factor is normalized by the factor σ0 h. The Rayleigh wave speed is cr = 2,947 m/s, and the dilatational wave speed is c1 = 5,944 m/s. First, we consider a stationary crack, i.e., the velocity is imposed to be zero during the simulation. We compare to the analytical relation for the stress intensity factor K 1 of the stationary crack [40]: K 1 (a, ˙ t) =

0 2σ0 1−ν

if t < tc , c1 (t−tc )(1−2ν) π

if tc ≤ t,

(52)

Second, let us now consider a moving crack, whose velocity is imposed to be zero until 1.5tc , and 1, 500 m/s after. The analytical relation linking the stress intensity factor K 1 and the velocity a˙ of the crack to the time t is [40]:

K 1 (a, ˙ t) =

⎧ ⎪ 0 ⎪ ⎪ ⎨ 2σ0 c1 (t−tc )(1−2ν) π

⎪ ⎪ c1 (t−tc )(1−2ν) 2σ0 ⎪ ⎩ 1−ν π 1−ν

if t < tc , if tc ≤ t < 1.5tc , 1− ca˙r

1− 2ca˙r

(53)

if 1.5 tc ≤ t.

Figure 7 presents the normalized mode 1 stress intensity factor as a function of time for the stationary and moving crack. Four results are illustrated: the analytical, the one obtained with the standard Heaviside function without any tip enrichment, and two results obtained with different tip enrichments (i.e., r sin(θ/2) and r 2 sin(θ/2)). In addition, Figs. 7c, d show the results obtained with tip enrichment and meshfree enrichment near the crack tip, whereas there is no meshfree enrichment in the result presented in Fig. 7a, b. For the stationary crack, one notices that the results obtained with the additional meshfree enrichment are improved compared to the one with or without tip enrichment and without meshfree. Both results obtained with meshfree and tip enrichment agree well with the analytical one. However, a slightly better result is obtained with the quadratic tip enrichment (i.e., r 2 sin(θ/2)), which in the later stages is almost exact. Figure 8 presents the relative error on the stress intensity factor as a function of time corresponding to the results shown in Figs. 7b, d. It can be seen that the meshfree part somewhat improves the results, but the improvement is limited. For the moving crack, Menouillard et al. [16] previously computed the stress intensity factor with only the quadratic tip enrichment. The result was slightly improved. In another way to deal with the oscillations occurring when the crack propagates, Menouillard and Belytschko [17,18] proposed a method to release smoothly the crack tip element using XFEM and only discontinuous enrichment. The results were significantly improved too. Here, the meshfree enrichment handles the description of the crack tip vicinity during propagation. One can notice again that the results obtained with the tip enrichment and meshfree enrichment are improved when the crack is moving, i.e., after t = tc (see Fig. 7d). The relative error is much lower for the computations with tip and meshfree enrichment than for the computation without or with tip enrichment without meshfree. To conclude, the use of a meshfree part in the discretization improves the accuracy of the stress intensity factor whether or not the crack is moving. The description is enhanced near the crack tip, thus the stress wave due to releasing tip element is decreased when propagation occurs.

64

T. Menouillard, T. Belytschko

(a)

(b)

(c)

(d)

Fig. 7 Normalized stress intensity factor as a function of time: a stationary crack with XFEM, b moving crack with XFEM, c stationary crack with XFEM + meshfree, and d moving crack with coupled XFEM + meshfree

(b) 50

No tip Zone + r sin(θ/2) MeshFree + r sin(θ/2)

Relative Error ΔK/K (in %)

Relative Error ΔK/K (in %)

(a) 50 40

30

20

40

30

20

10

10

0 1.6

No tip Zone + r2 sin(θ/2) MeshFree + r2 sin(θ/2)

0 1.7

1.8

1.9

2

2.1

t / tc

2.2

2.3

2.4

2.5

1.6

1.7

1.8

1.9

2

2.1

2.2

2.3

2.4

2.5

t / tc

Fig. 8 Relative error on the mode 1 stress intensity factor as a function of time for the use of the tip enrichment (a) r sin(θ/2) and b r 2 sin(θ/2), with and without the meshfree corresponding to Fig. 7b, d

6.2 Mode 2 stationary crack This example deals with a stationary mixed mode crack in a semi-infinite body, and focuses on the computation of the two dynamic stress intensity factors K 1 and K 2 . Figure 9 presents the geometry and the loading. The dimensions are the following: H = 6 m, a = 1 m, and L = 4 m. The prescribed velocity on the boundary is V0 = 16.5 m/s. Two meshes are used: 40 × 59 and 60 × 119 four-node elements. The material properties of the linear elastic media are: Young’s modulus E = 200 GPa, Poisson’s ratio ν = 0.25 and the density ρ = 7,833 kg/m3 . An analytical solution for both stress intensity factors as a function of time is available for this problem [44,45]. The numerical results are valid only until the compressive stress wave reaches again the crack tip after having been reflected on the right edge of the structure. The time needed for the wave to first

Dynamic fracture with meshfree enriched XFEM

65

Fig. 9 Geometry and loading 3

(b)

Analytical Standard Heaviside Meshfree + r2 sin(θ/2)

2.5

Normalized Stress Intensity Factors

Normalized Stress Intensity Factors

(a)

K2

2 1.5 1 K1 0.5 0 -0.5

3

Analytical Standard Heaviside Meshfree + r2 sin(θ/2)

2.5

K2

2 1.5 1 K1 0.5 0 -0.5

0

0.5

1

1.5

t / tc

2

2.5

3

0

0.5

1

1.5

2

2.5

3

t / tc

Fig. 10 Normalized stress intensity factors as functions of time, analytical and numerical results. a 40 × 59 element mesh and b 60 × 119 element mesh

reach the crack tip is given by tc = a/c1 where c1 is the dilatational wave speed. Therefore, the computation √ V0 a/π is run until 3 tc . The stress intensity factors are normalized by the factor −E . Belytschko and Tabbara 2 c1 (1−ν 2 ) [46] simulated this problem with the element free Galerkin method, and obtained a good agreement with the analytical result. Figure 10a, b shows both stress intensity factors K 1 and K 2 as functions of time for the coarse mesh (i.e., 40 × 59 elements) and the fine mesh (i.e., 60 × 119 elements), respectively. The three different numerical results globally agree with the analytical one. However, one can notice significant oscillations appearing for the computation performed on the coarse mesh. There are still some for the finer mesh. These are due to the time integration (explicit scheme) and the boundary condition (velocity imposed). Second, the computations for both meshes underline that the results are improved with the tip enrichment zone. In fact, the most accurate results are obtained using the quadratic meshfree enrichment. The computation without any tip enrichment tends to overestimate the stress intensity factors. This is due to the fact that the mechanical crack tip is always further than the virtual crack tip in the center of the J-integral domain. Once again, the accuracy in terms of stress intensity factors is slightly improved using the proposed discretization.

6.3 Crack branching problem Crack branching is a common phenomenon in dynamic crack propagation, and has been experimentally investigated [47–54]. However, only a few numerical results have been reported [14,36,55–57]. Experiments [50,51] led to the conclusion that branching occurred when the stress intensity factor was around two to three times the critical fracture toughness. Figure 11 presents the geometry and the loading of the specimen. The dimensions of the structure are the following: the length is L = 100 mm, the initial crack length a = 50 mm, and the vertical position of the crack is h = 40 mm. The material properties of the linear elastic media are: Young’s modulus E = 32 GPa, Poisson’s

66

T. Menouillard, T. Belytschko

Fig. 11 Geometry and loading of the crack branching problem

(a) 0.7

(b)

XFEM XFEM + Meshfree

1

XFEM XFEM + Meshfree [Song2009]

0.6 0.8

0.4

Vtip / cr

K (MPa√m)

0.5

0.3

0.6

0.4

0.2 0.1

0.2 0 -0.1

0 0

10

20

30

40

50

60

Time (μs)

0

10

20

30

40

50

60

Time (μs)

Fig. 12 a Mode 1 Stress intensity factors √ as a function of time for the numerical computations using only XFEM and the coupling XFEM/meshfree (K 1c = 0.33 MPa m), and b crack tip velocities as a function of time for the numerical computations using only XFEM and the coupling XFEM/meshfree, compared to [58]

ratio ν = 0.2 and the density ρ = 2,450 kg/m3 . The Rayleigh wave speed is cr = 2,125 m/s. The cohesive law from Sect. 5.1√was used with a fracture energy G F = 3.0 N/m, which corresponds to a fracture toughness K I c = 0.33 MPa m. The tensile stress applied on the top surface is σ¯ = 1 MPa. The mesh is 166 × 79 and consists of four-node quadrilaterals. As before, the crack velocity was driven by Eq. (50). The criterion for branching was first chosen to be Döll’s [49] criterion G > 2Gc .

(54)

We also monitored the stress before and after Döll’s criterion was met. Linder and Armero [57] have used a critical velocity to trigger branching. We note that the cohesive dissipation needs to be added here because in the absence of a singularity at the crack tip, there is no dissipation (in the previous two problems the crack tip was advanced arbitrarily). Although the Freund formula is not strictly applicable in this case, it provides a reasonable law for crack speed since it has experimental support [43]. Figure 12a presents the mode 1 stress intensity factor as a function of time for both computations, i.e., XFEM and XFEM + meshfree. One can notice that both results are similar, and present a maximum value for time 30 µs, when the crack branching occurs. The mode 1 stress intensity factor begins to increase when√the tensile stress wave reaches the crack tip. It increases until reaching the√ fracture toughness K 1c = 0.33 MPa m. Then the crack propagates keeping a stress intensity around 0.38 MPa m. At the time 25 µs, the stress intensity √ factor increases suddenly to reach almost √ 0.57 MPa m at time 32 µs. The ratio between the stress intensity factor and the fracture toughness is 2 at 30–31 √ µs, and crack branching occurs. Subsequently, the stress intensity factor decreases to around 0.40 MPa m and keeps this value. Figure 13a presents the position of the quadrature points near the crack tip at which the stress is monitored. Figure 13b shows the hoop stress at the quadrature point near crack tip as a function of its angular position for four different times (5.7, 12.7, 19.7 and 28.5 µs) for the computation using the XFEM + meshfree discretization. Note that the branching occurs after time 28.5 µs. One notices that the shape of the hoop stress near

Dynamic fracture with meshfree enriched XFEM

(a)

2

67

(b)

Quadrature points

6

Time 5.7μs Time 12.7μs Time 19.7μs Time 28.5μs

1.5 5 1 4

σθθ (MPa)

Y (mm)

0.5 0

3

-0.5 2 -1 1 -1.5 -2 50.5

0 51

51.5

52

52.5

53

53.5

54

-80

-60

-40

-20

0

20

40

60

80

Angle (deg)

X (mm)

Fig. 13 a Quadrature points near the crack tip at which the stress is monitored, and b Hoop stress at the quadrature points as a function of the polar angle from the crack tip for different time of the propagation before branching

(a)

Branching (XFEM) Branching (XFEM+Meshfree)

30

(b)

20

Y (mm)

10 0 -10 -20 -30 0

20

40

60

80

100

X (mm)

Fig. 14 a Crack paths for the numerical computations using only XFEM and the coupled XFEM/meshfree, and b the final deformed specimen (amplified ×200)

the crack tip has two major maxima, which lead to the crack branching. One objective of this computation is to see whether a more accurate computation of the stresses in the vicinity of the crack tip will enable crack branching criteria to depend on the stresses and whether the stresses will give a reasonable angle of branching. We also examined how closely a criterion based on Gc agrees with the appearance of two stress maxima. Figure 12b shows the crack speed as a function of time for both computations, i.e., XFEM and XFEM + meshfree. The result of Song and Belytschko [58] is also represented in Fig. 12b. Note that the crack speed is related to the stress intensity through Eq. (50) in our computations. The branching occurs when the crack speed reaches around 90% of the Rayleigh wave speed. Figure 14 presents the crack path for both computations. Both results are quite similar, and show a good agreement with [58].

7 Conclusion A method to enrich the XFEM by meshless approximations has been described. The advantage is that the meshfree approximation can easily be added and removed, so it lends itself to adaptivity. As the crack moves, it can be added around the new areas around the crack tip, and removed in the subdomains the crack has passed. The meshfree part of the discretization aims at smoothing the stress state near the crack tip during the propagation, and thus decreasing the unphysical oscillations in the stress due to the propagation of the discontinuity.

68

T. Menouillard, T. Belytschko

By the numerical examples, we showed that the accuracy of the stress intensity factor has been improved for stationary and moving crack, although the improvement is moderate. Furthermore, the computations show that the improved accuracy in stresses enables the branching point to be identified from the stresses. Branching criteria in terms of the fracture toughness are based on experiments, and it is not known how applicable they are to many materials. The results obtained here clearly show that with increasing crack velocity, the stresses develop two maxima. Such results have also analytically been demonstrated by Yoffe [47]. In the √ case of our computations, this also coincides with the time when the stress intensity factor equals 0.52 MPa m, which is a criterion that has been deduced from experiment. References 1. Belytschko, T., Krongauz, Y., Organ, D., Fleming, M., Krysl, P.: Meshless methods: an overview and recent developments. Comput. Methods Appl. Mech. Eng. 139, 3–47 (1996) 2. Belytschko, T., Krongauz, Y., Dolbow, J., Gerlach, C.: On the completeness of meshfree particle methods. Int. J. Numer. Methods Eng. 43, 785–819 (1998) 3. Fleming, M., Chu, Y.A., Moran, B., Belytschko, T.: Enriched element-free Galerkin methods for crack tip fields. Int. J. Numer. Methods Eng. 40, 1483–1504 (1997) 4. Nguyen, V.P., Rabczuk, T., Bordas, S., Duflot, M.: Meshless methods: a review and computer implementation aspects. Math. Comput. Simul. 79, 763–813 (2008) 5. Moës, N., Dolbow, J., Belytschko, T.: A finite element method for crack growth without remeshing. Int. J. Numer. Methods Eng. 46, 131–150 (1999) 6. Dolbow, J., Moës, N., Belytschko, T.: Discontinuous enrichment in finite elements with a partition of unity method. Int. J. Numer. Methods Eng. 36, 235–260 (2000) 7. Melenk, J.M., Babuška, I.: The partition of unity finite element method: basic theory and applications. Comput. Methods Appl. Mech. Eng. 139, 289–314 (1996) 8. Belytschko, T., Moës, N., Usui, S., Parimi, C.: Arbitrary discontinuities in finite elements. Int. J. Numer. Methods Eng. 50, 993–1013 (2001) 9. Moës, N., Belytschko, T.: Extended finite element method for cohesive crack growth. Eng. Fract. Mech. 69, 813–833 (2002) 10. Moës, N., Gravouil, A., Belytschko, T.: Non-planar 3D crack growth by the extended finite element and level sets. Part I : Mechanical model. Int. J. Numer. Methods Eng. 53, 2549–2568 (2002) 11. Duan, Q., Song, J.H., Menouillard, T., Belytschko, T.: Element-local level set method for three-dimensional dynamic crack growth. Int. J. Numer. Methods Eng. 80, 1520–1543 (2009) 12. Wells, G.N., Sluys, L.J., De Borst, R.: Simulating the propagation of displacement discontinuities in a regularized strainsoftening medium. Int. J. Numer. Methods Eng. 53, 1235–1256 (2002) 13. Réthoré, J., Gravouil, A., Combescure, A.: A stable numerical scheme for the finite element simulation of dynamic crack propagation with remeshing. Comput. Methods Appl. Mech. Eng. 193, 4493–4510 (2004) 14. Belytschko, T., Chen, H., Xu, J., Zi, G.: Dynamic crack propagation based on loss of hyperbolicity and a new discontinuous enrichment. Int. J. Numer. Methods Eng. 58, 1873–1905 (2003) 15. Belytschko, T., Chen, H.: Singular enrichment finite element method for elastodynamic crack propagation. Int. J. Comput. Methods 1, 1–15 (2004) 16. Menouillard, T., Song, J.H., Duan, Q., Belytschko, T.: Time dependent crack tip enrichment for dynamic crack propagation. Int. J. Fract. doi:10.1007/s10704-009-9405-9 (2009) 17. Menouillard, T., Belytschko, T.: Correction Force for releasing crack tip element with XFEM and only discontinuous enrichment. Eur. J. Comput. Mech. 18, 465–483 (2009) 18. Menouillard, T., Belytschko, T.: Improvement with smoothly discontinuous enrichment in XFEM for dynamic crack propagation. Int. J. Numer. Methods Eng. (2009, accepted) 19. Xiao, Q.Z., Karihaloo, B.L.: Improving the accuracy of XFEM crack tip fields using higher order quadrature and statically admissible stress recovery. Int. J. Numer. Methods Eng. 66, 1378–1410 (2006) 20. Bordas, S., Duflot, M.: Derivative recovery and a posteriori error estimate for extended finite elements. Comput. Methods Appl. Mech. Eng. 196, 3381–3399 (2007) 21. Duflot, M., Bordas, S.: A posteriori error estimation for extended finite elements by an extended global recovery. Int. J. Numer. Methods Eng. 76, 1123–1138 (2008) 22. Ródenas, J.J., González-Estrada, O.A., Tarancón, J.E., Fuenmayor, F.J., Valenciana, G.: A recovery type error estimator for the extended finite element method based on singular+ smooth stress field splitting. Int. J. Numer. Methods Eng. 76, 545–571 (2008) 23. Liu, W.K., Uras, R.A., Chen, Y.: Enrichment of the finite element method with the reproducing kernel particle method. Trans. ASME J. Appl. Mech. 64, 861–870 (1997) 24. Huerta, A., Fernandez-Mendez, S.: Enrichment and coupling of the finite element and meshless methods. Int. J. Numer. Methods Eng. 48, 1615–1636 (2000) 25. Fernandez-Mendez, S., Diez, P., Huerta, A.: Convergence of finite elements enriched with mesh-less methods. Numer. Math. 96, 43–59 (2003) 26. Rabczuk, T., Xiao, S.P., Sauer, M.: Coupling of meshfree methods with finite elements: basic concepts and test results. Commun. Numer. Methods Eng. 22, 1031–1065 (2006) 27. Gu, Y.T., Zhang, L.C.: Coupling of the meshfree and finite element methods for determination of the crack tip fields. Eng. Fract. Mech. 75, 986–1004 (2008) 28. Belytschko, T., Liu, W.K., Moran, B.: Nonlinear Finite Elements for Continua and Structures. Wiley, Chichester (2000)

Dynamic fracture with meshfree enriched XFEM

69

29. Cox, J.V.: An extended finite element method with analytical enrichment for cohesive crack modeling. Int. J. Numer. Methods Eng. 78, 48–83 (2009) 30. Réthoré, J., Gravouil, A., Combescure, A.: An energy-conserving scheme for dynamic crack growth using the extended finite element method. Int. J. Numer. Methods Eng. 63, 631–659 (2005) 31. Menouillard, T., Réthoré, J., Moës, N., Combescure, A., Bung, H.: Mass lumping strategies for X-FEM explicit dynamics: application to crack propagation. Int. J. Numer. Methods Eng. 74, 447–474 (2008) 32. Menouillard, T., Réthoré, J., Combescure, A., Bung, H.: Efficient explicit time stepping for the eXtended Finite Element Method (X-FEM). Int. J. Numer. Methods Eng. 68, 911–939 (2006) 33. Elguedj, T., Gravouil, A., Maigre, H.: An explicit dynamics extended finite element method. Part 1: mass lumping for arbitrary enrichment functions. Comput. Methods. Appl. Mech. Eng. 198, 2297–2317 (2009) 34. Barenblatt, G.I.: The mathematical theory of equilibrium cracks formed in brittle fracture. Adv. Appl. Mech. 7, 55–129 (1962) 35. Hillerborg, A., Modeer, M., Petersson, P.E.: Analysis of crack formation and crack growth in concrete by means of fracture mechanics and finite elements. Cement Concrete Res 6, 773–782 (1976) 36. Song, J.H., Areias, P.M.A., Belytschko, T.: A method for dynamic crack and shear band propagation with phantom nodes. Int. J. Numer. Methods Eng. 67, 868–893 (2006) 37. Rabczuk, T., Song, J.H., Belytschko, T.: Simulations of instability in dynamic fracture by the cracking particles method. Eng. Fract. Mech. 76, 730–741 (2009) 38. Rice, J.R.: A path independent integral and the approximate analysis of strain concentration by cracks and notches. J. Appl. Mech. (Trans. ASME) 35, 379–386 (1968) 39. Bui, H.D.: Mecanique de la rupture fragile. Masson, Paris (1978) 40. Freund, L.B.: Dynamic Fracture Mechanics. Cambridge University Press, London (1990) 41. Freund, L.B.: Crack propagation in an elastic solid subjected to general loading. Pt. 1. Constant rate of extension. J. Mech. Phys. Solids 20, 129–140 (1972) 42. Freund, L.B., Douglas, A.S.: Influence of inertia on elastic-plastic antiplane-shear crack growth. J. Mech. Phys. Solids 30, 59–74 (1982) 43. Rosakis, A.J., Freund, L.B.: Optical measurement of the plastic strain concentration at a crack tip in a ductile steel plate. J. Eng. Mater. Technol. (Trans. ASME) 104, 115–120 (1982) 44. Ravi-Chandar, K.: Dynamic Fracture. Elsevier, Amsterdam (2004) 45. Lee, Y.J., Freund, L.B.: Fracture initiation due to asymmetric impact loading of an edge cracked plate. J. Appl. Mech. 57, 104 (1990) 46. Belytschko, T., Tabbara, M.: Dynamic fracture using element-free Galerkin methods. Int. J. Numer. Methods Eng. 39, 923–938 (1996) 47. Yoffe, E.H.: The moving Griffith crack. Philos. Mag. 42, 739–750 (1951) 48. Schardin, H.: Velocity effects in fracture. In: Averbach, B., et al. (eds.) Fracture, pp. 297–330. Wiley, New York (1959) 49. Döll, W.: Investigations of the crack branching energy. Int. J. Fract. 11, 184–186 (1975) 50. Kobayashi, A.S., Wade, B.G., Bradley, W.B.: Fracture dynamics of Homalite-100. In: Kansch, H.H., Haseele, J.A., Jaffee, R.L. (eds.) Deformation and Fracture of High Polymers, pp. 487–500. Plenum Press, New York (1973) 51. Kobayashi, A.S.: Photoelastic Techniques. Exp. Tech. Fract. Mech., 126–145 (1973) 52. Sumi, Y., Nemat-Nasser, S., Keer, L.M.: On crack branching and curving in a finite body. Int. J. Fract. 21, 67–79 (1983) 53. Ramulu, M., Kobayashi, A.S.: Mechanics of crack curving and branching dynamic fracture analysis. Int. J. Fract. 27, 187– 201 (1985) 54. Ravi-Chandar, K.: Dynamic fracture of nominally brittle materials. Int. J. Fract. 90, 83–102 (1998) 55. Xu, X.P., Needleman, A.: Numerical simulations of fast crack growth in brittle solids. J. Mech. Phys. Solids 42, 1397–1434 (1994) 56. Seelig, T., Gross, D.: On the interaction and branching of fast running cracks—a numerical investigation. J. Mech. Phys. Solids 47, 935–952 (1999) 57. Linder, C., Armero, F.: Finite elements with embedded branching. Finite Elem. Anal. Des. 45, 280–293 (2009) 58. Song, J.H., Belytschko, T.: Cracking node method for dynamic fracture with finite elements. Int. J. Numer. Methods Eng. 77, 360–385 (2009)