An analytical stiffness derivative extended finite element technique for extraction of crack tip Strain Energy Release Rates

Engineering Fracture Mechanics 77 (2010) 3204–3215

Contents lists available at ScienceDirect

Engineering Fracture Mechanics journal homepage: www.elsevier.com/locate/engfracmech

An analytical stiffness derivative extended finite element technique for extraction of crack tip Strain Energy Release Rates Haim Waisman Department of Civil Engineering & Engineering Mechanics, Columbia University, New York, NY 10027, United States

a r t i c l e

i n f o

Article history: Received 19 January 2010 Received in revised form 16 August 2010 Accepted 17 August 2010 Available online 22 August 2010 Keywords: Extended finite element method XFEM Stiffness derivative J-integral Mixed mode fracture Stress intensity factors Strain Energy Release Rates

a b s t r a c t An analytical approach combined with the extended finite element method (XFEM) is proposed to extract the Strain Energy Release Rates within the classical stiffness derivative technique. The proposed idea hinges on the following two XFEM properties: (i) the crack is mesh independent, i.e. there is no need for mesh perturbations in the vicinity of the crack and (ii) the asymptotic crack tip field is embedded in the mathematical formulation of the stiffness matrix. By employing these properties we show that the derivative of the stiffness matrix with respect to the crack extension can be computed in a closed form and on the fly during the analysis. Thus the virtual crack extension, and the error inherent in the finite difference scheme of the classical stiffness derivative technique is completely avoided. Numerical results on few benchmark problems show that this method is comparable to the J-integral method. Ó 2010 Elsevier Ltd. All rights reserved.

1. Introduction Several methods have been proposed in the literature to extract the Stress Intensity Factors (SIFs) or equivalently the Strain Energy Release Rates (SERRs) in linear elastic fracture mechanics. These methods can be grouped as direct and indirect methods. Direct methods are the simplest methods to compute SIFs, which are mainly based on correlation of crack opening displacements, obtained directly from the FE solution [1,2]. Indirect methods are based on integral or energy quantities to calculate the SERRs and often require additional steps during the FE analysis process or a special postprocessing, however yield more accurate results [3,4]. Some of the most well known and used indirect methods are the J-integral method [5], its domain (area) variant [4,6,7], the stiffness derivative method [8,9] and the virtual crack closure techniques which are based on Irwin’s integral [10,11]. These methods have received tremendous amount of attention in the literature in the past few decades and some have also been implemented in commercial software. Nevertheless, due to the global character of these energetic methods it is sometimes not easy to extract independent estimations of the SIF associated with each crack opening mode, KI, KII and KIII in mixed mode problems, e.g. complicated crack fronts in three dimensions or cracks in nonhomogeneous bi materials. The classical stiffness derivative method proposed by Parks [8] and extended by Hellen [9] employs finite difference schemes to compute the derivative of the stiffness matrix with respect to a virtual extension in crack length. The method has been found very accurate, within 1% of analytical solutions, with a comparable accuracy in two and three dimension to the J-integral method [3,12]. Nevertheless, it approximates the stiffness derivative by subtracting two stiffness matrices divided by a finite infinitesimal crack increment which may introduce numerical truncation errors. To improve this method Lin and Abel [13] introduced an approach for evaluation of stiffness derivatives directly from the finite element formulations. E-mail address: [email protected] 0013-7944/$ - see front matter Ó 2010 Elsevier Ltd. All rights reserved. doi:10.1016/j.engfracmech.2010.08.015

H. Waisman / Engineering Fracture Mechanics 77 (2010) 3204–3215

3205

Nomenclature x spatial coordinates a, Da crack length, and virtual crack extension r the distance from a point to the crack tip h angle from a point and crack tip to the crack line u, d coefficient vectors corresponding to standard and all degrees of freedom a, b coefficient vectors corresponding to jump and tip enriched degrees of freedom N and B shape functions and strain–displacement matrix J Jacobian of the transformation D stress–strain constitutive matrix E, m and b Young’s modulus, Poisson ratio and a constant Fj four Branch functions H Heaviside function W and U Levelset functions G Strain Energy Release Rate P total potential energy K stiffness matrix M mutual energy release rate KI and KII stress intensity factors

This later served as the basis for higher order derivatives calculations in multiple cracks analysis [14,15], and in a design sensitivity framework which yielded more accurate results [16]. However, these methods are based on standard finite element formulations where the implementation requires perturbation of elements surrounding the crack tip. As illustrated in Fig. 1, there are many possibilities to obtain a virtual crack extension (two are shown in the figure) and it is therefore not straight forward to asses which one yields better results. In addition, modeling cracks with standard FEM is quite limited due to special mesh considerations which requires the mesh to conform with the crack surfaces. Furthermore, for better accuracy, quarter point elements should be placed at the crack tip [17,18]. Thus if a crack is propagating in a structure in an arbitrary direction, mesh generation is not always trivial. Recently, an extended finite element method (XFEM) has been proposed by Belytschko and co-workers [19–21] and a generalized finite element method (GFEM) by Strouboulis et al. [22,23]. These methods have similar features and in this paper we will focus on the XFEM. The main outcome and attractiveness of XFEM is that it enables to model cracks without the need for special meshes, i.e. cracks can cut through elements, which implies that mesh generation may be avoided. The key idea of XFEM is to locally enrich the standard FE shape functions with Heaviside functions to enable opening displacepffiffiffi ments, and asymptotic functions at the tip element, which incorporate a r term in the displacement field, and thus result in

(a)

active elements

Δa

(b)

active elements

Δa

Fig. 1. Surrounding crack elements that contribute to the stiffness derivative matrix due to a virtual crack extension Da [3]. The crack is illustrated by red line and active elements are colored in gray. Two different possibilities are shown: (a) crack tip elements and (b) first ring of elements. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

3206

H. Waisman / Engineering Fracture Mechanics 77 (2010) 3204–3215

a stress singularity at the crack tip. We will review the XFEM in a later section. Other methods based on partition of unity concepts are given in [24–26]. In the current paper we propose an analytical method combined with the XFEM to compute the stiffness derivative without the need for virtual crack extensions. The proposed approach employs the formulations proposed by Lin and Abel [13], which avoids the classical finite difference approximation of the stiffness derivative matrix [8,9,12]. Specifically, the idea hinges on the following two properties of the XFEM: (i) the underlying mesh is fixed and no perturbation of mesh is required for the virtual extension and (ii) the crack is embedded in the mathematical formulation of the stiffness matrix and thus derivatives may be computed directly. Combining these features yield a simple and straight forward algorithm for extracting the Strain Energy Release Rates. To extract the mixed mode SERRs we employ Yau’s mutual potential energy representation [27], which is based on Betti’s reciprocal theorem. This method employs an analytical solution in the vicinity of crack tip to decompose the SERRs into components. A similar decomposition is also employed by the interaction integral which decomposes the J-integral and extracts the SIF’s in mixed mode problems [28,7,29]. The reminder of the paper is organized as follows. In Section 2, we present a brief introduction to the XFEM. In Section 3 we describe the stiffness derivative method and the proposed method that is based on the XFEM. Section 4 discusses Yau’s mutual potential energy to extract the mixed mode SERRs followed by several benchmark examples of 2D pure and mixed mode problems. 2. Modeling cracks by the extended finite element method The key idea of XFEM is to locally enrich the standard finite element approximation with local partitions of unity enrichment functions which are chosen according to the physics of the problem at hand. It follows that for crack problems, the mesh is independent of the crack orientation [30–32]. An excellent review on the XFEM can be found in [33]. Similar enrichment methods for modeling cracks are based on the generalized finite element method [34,35]. Let uh 2 Uh be an extended finite element approximation to the descretized weak form of elasticity, where Uh is the appropriate Sobolev space [36]. The XFEM enriches the conventional shape function space with a set of functions H(x) and Fj(x), such that

uh ðxÞ ¼

n X

N I ðxÞuI þ

I¼1

nJ X

N I ðxÞHðxÞaI þ

I¼1

nT X

" N I ðxÞ

I¼1

4 X

# F j ðxÞbjI

ð1Þ

j¼1

where x = {x, y, z}Tare the spatial coordinates. n, nJ and nT are the number of total nodal points, jump-enriched nodes and tip enriched nodes, respectively. NI(x) are the standard finite element shape functions associated with standard degrees of freedom uI, whereas aI and bjI are the degrees of freedom associated with the enriched nodes. For linear elastic fracture problems, the crack tip zone is enriched with the classical analytical solution for the near tip field [37], given by

F j ðrðxÞ; hðxÞÞ ¼

  pffiffiffi pffiffiffi h pffiffiffi h pffiffiffi h h r sin ; r cos ; r sin cos h; r cos cos h 2 2 2 2 j¼1;...;4

ð2Þ

while the rest of the element nodes cut by the crack are enriched with the Heaviside function

( HðxÞ ¼

þ1 above Cþc

ð3Þ

1 below Cc

pffiffiffi  where Cþ r term in c and Cc defines the edges of the discontinuity line that splits the element into two parts. Note that the 1 the displacement field is directly built into the equations and hence the stress singularity of pffir appears in the solution. The enriched nodes of an inclined 2D crack and the corresponding deformed shape under tension traction loads is illustrated in Fig. 2. In the illustration the Young’s modulus and Poisson ratio are set to E = 1  107 [N/m2] and m = 0.3, respectively. The plate is fixed at the bottom and pulled upwards by a traction of t ¼ 100 ½N=m2 . Typically, these enriched nodes (and the corresponding elements) are easily obtained with the aid of the Levelsets method. The Levelset method is a numerical technique for representing discontinuities and boundaries, in particular when they are moving [38,39]. In the context of finite elements, this method has been found very useful as information of the crack geometry is only stored as additional nodal variables and then interpolated to gauss (integration) points employing standard shape functions. It is also straight forward to update this stored data as the crack propagates in the structure. Moreover, the Levelset method allows an automatic detection of enriched nodes and the representation of the enrichment functions in Eq. (1) by means of Levelset values. For crack problems one may define two Levelset functions W and U which are simply the distance in local crack coordinate from any point to the crack front and crack surface, respectively. Thus the angle h and distance r given in the enrichment functions in Eq. (2), may be obtained by

r¼

  qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi U W2 þ U2 and h ¼ tan1

W

For more details the reader is referred to [32,40].

ð4Þ

H. Waisman / Engineering Fracture Mechanics 77 (2010) 3204–3215

3207

Fig. 2. XFEM modeling of an inclined crack. Traction B.C (tension) are applied on the top edge and the bottom edge is fixed. (a) Enrichments visualization: Crack line is illustrated in red, black circles are jump enriched nodes and green squares are tip enriched nodes. (b) Deformed configuration (1.5  104). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

3. The classical stiffness derivative method and the proposed approach The Stiffness derivative method [8,9], which belongs to the category of virtual crack extension methods, is an energy approach to extract the energy release rates and is considered as one of the most accurate and efficient methods [3,12]. The idea is to compute the rate of change in total potential energy (the Strain Energy Release Rate under LEFM) of the system due to a virtual (infinitesimal) extension of the crack in a self similar direction. For 2D unit thickness problems, in the absence of body forces, the SERR (G) is given by

G¼

  @P @ 1 T T d Kd  d f ¼ @a 2 @a

ð5Þ

where P is the total potential energy, a is the crack length, K is the global stiffness matrix, d is the global nodal displacement and f is the global nodal force vector. Assuming that the applied forces do not change due to the virtual crack extension reduces Eq. (5) to

G¼

@P 1 T @K d ¼ d 2 @a @a

ð6Þ

Parks [8] proposed to approximate the derivative @K/@a by a forward difference, such that

›K 

@K K aþDa  K a  @a Da

ð7Þ

where Da is an infinitesimal crack extension and the resulting matrix is defined as the symbol ›K. The global matrix ›K is obtained by assembly of local element contributions

›K ¼ Ae2g ›K e  Ae2g

Ke  K ea @K e  Ae2g aþDa @a Da

ð8Þ

where A is the assembly operator over a set of elements {g} and superscript e denotes individual elements. Note that the element set {g} consists only those elements, in the vicinity of the crack tip, which will have some nonzero contribution to ›K due to the virtual crack extension. Two different possibilities for the element set {g} were illustrated in Fig. 1. This finite difference approach proposed by Parks [8] is a simple and straight forward way to compute the SERRs. Nevertheless since it is based on a difference between two matrices, it may introduce errors due to the non unique way to determine the increment Da. Furthermore, it is sometimes important to compute higher order derivatives (first and second) of the SERRs or SIF’s and numerical differentiation may not be accurate [14,15,41]. Higher order derivatives are important quantities for the prediction of the stability and arrest of a single crack [42], propagation and interaction of system of cracks [43,44], perturbation analysis of crack trapping [45], probabilistic fracture mechanics analysis [46,47] and the universal size effect model that relates nominal strength to the Fracture Process Zone [48]. In the current paper we follow the explicit derivation given in [49] and [13], where it was shown that for isoparametric element mappings with the following stiffness matrix

Ke ¼

Z Xe

BT DBjJjdXe

ð9Þ

3208

H. Waisman / Engineering Fracture Mechanics 77 (2010) 3204–3215

the derivatives with respect to the crack length a can be obtained directly using the chain rule, i.e.

" # @BT @jJj T @B T dXe jJj þ B DB ›K ¼ DBjJj þ B D @a @a @a Xe e

Z

ð10Þ

where B is the strain–displacement matrix, D is the stress–strain constitutive matrix, J is the Jacobian of the transformation between the natural to physical coordinate system and Xe is the element integration domain in the natural system. It is assumed that the constitutive relations do not depend on the crack length. Several researchers have used these formulations to develop improved stiffness derivative methods. For instance Barbero and Reddy [50] developed a Jacobian derivative method as a postprocessing tool to evaluate the SERRs and Giner et al. [16] proposed an interesting approach based on shape design sensitivity. Nevertheless, these methods are inherently mesh dependent as they are based on the standard FEM, and hence have to compute the term @jJj/@a in Eq. (10) by some means of mesh perturbation around the crack tip, which may not be trivial. The key point of the proposed approach is that the term @jJj/@a  0, or

Z X

BT DB e

@jJj e dX  0 @a

ð11Þ

Eq. (11) holds due to the mesh independence in XFEM (see Section 2). In other words, the mesh is left unchanged due to a virtual extension of Da. Thus the stiffness derivative matrix reduces to

›K exfem

" # @BT T @B jJj dXe ¼ DBjJj þ B D @a @a Xe Z

ð12Þ

Moreover, the strain–displacement relation (B matrix) in XFEM is obtained by differentiating (1), i.e.

 nJ nT X n 4  X X @uh ðxÞ X N I ðxÞ @N I ðxÞ @N I ðxÞ @F j ðxÞ bjI ¼ uI þ HðxÞaI þ F j ðxÞ þ N I ðxÞ @x @x @x @x @x I¼1 I¼1 j¼1 I¼1 h

tip jump ¼ Bfem u þ Bjump xfem a þ Bxfem b ¼ Bfem Bxfem

2 3 i u 6 7 e Btip xfem 4 a 5 ¼ Bd

ð13Þ

ð14Þ

b where B is the overall strain–displacement matrix and de is the element vector with all the unknown coefficients. Bfem is the part of the B matrix and u is the part of the coefficient vector de corresponding to conventional degrees of freedom. Bjump xfem is the part of the B matrix and a is the part of the coefficient vector d corresponding to the jump enriched degrees of freedom. Similarly, Btip xfem and b are the parts associated with the tip enriched degrees of freedom. Next, taking a derivative of Eq. (14) with respect to the virtual crack extension leads to

" # i @Btip @B @ h xfem jump tip ¼ Bfem Bxfem Bxfem ¼ 0 0 @a @a @a

ð15Þ

where the only contribution in Eq. (15) comes from the tip enriched part, as shown in Fig. 3 (since the mesh do not change as the crack extends). Note however that one may enrich a bigger region (several element layers around the crack tip) with tip enriched nodes in which all will contribute to the derivative computation. Thus to obtain the derivative @Btip xfem =@a, one has to differentiate the tip functions in Eq. (2) with respect to @a, i.e. @Fj/@a and @ 2Fj/(@x@a). These functions may be expended by the chain rule, which is given by

@F j ðrðxÞ; hðxÞÞ @F j @r @F j @h ¼ þ @a @r @a @h @a     2 @ F j ðrðxÞ; hðxÞÞ @ @F j @ @F j @r @F j @h ¼ þ ¼ @x@a @x @a @x @r @a @h @a

ð16Þ ð17Þ

The derivatives in Eqs. (16) and (17) may be computed numerically or analytically. Either way these derivatives are computed on the fly during the element formulation. In the semi-analytical (numerical) approach one may employ, e.g. a central difference scheme, which leads to

@F j F j ðx þ DaÞ  F j ðx  DaÞ  2 Da @a     2 @ Fj @ @F j 1 @F j @F j ¼  ðx þ DaÞ  ðx  DaÞ 2Da @x @x@a @a @x @x

ð18Þ ð19Þ

We have chosen the central difference approach due to the higher accuracy (compared to forward difference) and also to differentiate these function evaluations from the classical stiffness derivative method [8]. Nonetheless, the choice of the

H. Waisman / Engineering Fracture Mechanics 77 (2010) 3204–3215

(a)

active elements

(b)

3209

active elements

Δa

Δa

Fig. 3. Crack tip elements which contribute to the stiffness derivative matrix due to a virtual crack extension Da. The crack is illustrated by a red line and active elements are colored in gray. Two different possibilities are shown: (a) closest tip elements enriched and (b) larger domain around the crack is enriched. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

increment Da is quite sensitive as illustrated in Fig. 8a and b. Thus, the analytical approach is preferable and we proceed with the closed form derivations. To compute the derivatives analytically, we first note that the derivatives @Fj/@r,@Fj/@h, @h/@x and @r/@x are already computed and available as part of the extended finite element method. Hence, one only needs to obtain closed form solutions for the derivatives @r/@a, @h/@a, @ 2Fj/@a@h, @ 2Fj/@a@r, which can be done through the definition of the Levelset functions in Eq. (4). The key point in the derivation is to observe that for a self similar crack extension in local coordinates

@W ¼ 1 and @a

@U ¼0 @a

ð20Þ

hence the derivatives @r/@a and @h/@a are obtained by

 

1=2  1

1=2  @r @ 2 @W @U W ¼ ¼ ¼ W þ U2 W2 þ U2 2W þ 2U @a @a 2 @a @a r

ð21Þ

#      "       @h @ U 1 @ U W2 @ 1 1 @U W2 U @W U 1 ¼ ¼ þ ¼ tan U  ¼ ¼ 2 @a @a @a W W r2 W @a r2 r W2 @a 1 þ ðU=WÞ2 @a W

ð22Þ

Plugging that back to Eqs. (16) and (17), along with the following notation simplification

@F j ¼ fjr @r

and

@F j ¼ fjh @h

ð23Þ

gives

@F j W U ¼ fjr  fjh 2 @a r r     @2Fj @ @F j @ U W ¼ fjh 2  fjr ¼ @x @x@a @x @a r r

ð24Þ ð25Þ

Without the loss of generality, but due to more simple relations, we will now focus on a 2D case. Noting that

@W ¼ 1; @x

@W ¼ 0 and @y

@U ¼ 0; @x

@U ¼1 @y

ð26Þ

the derivatives in x or y in Eq. (25) can further be expended by using the chain rule, and we arrive at

@fj @2Fj ¼ @x@a @x @fjh @2Fj ¼ @y@a @y h

  @fjr W W @r r 1   f j r r2 @x r r 3 @x @x r   @fjr W U U @r h 1 r W @r þ f  2 þ f  j j r2 r2 r 3 @y r2 @y @y r

U

 2fjh 2

U @r



ð27Þ ð28Þ

3210

H. Waisman / Engineering Fracture Mechanics 77 (2010) 3204–3215 @fjr

Here the chain rule is employed once more to compute the derivatives

@fjr

@fjr @r @fjr @h þ ¼ @x @r @x @h @x

@x

, that is

@fjh

@fjh @r @fjh @h þ ¼ @x @r @x @h @x

and

ð29Þ

Note that the derivatives @fjr =@r; @fjr =@h; @fjh =@r and @fjh =@r are essentially the second derivatives of the tip enrichment functions in Eq. (2) with respect to r and h, which in general are not computed in standard XFEM but can easily be computed (see Appendix A for complete derivations). To summarize, the stiffness derivative may be computed analytically following the outlined procedure in this Section. Hence the virtual crack extension proposed in the classical technique by Parks [8] may completely be avoided. Finally, we illustrate in Fig. 4a and b the sparsity pattern of a fully assembled stiffness derivative matrix. The examples correspond to the mixed mode case in Section 5.2. In these figures the total number of Equations are 628 (Fig. 4a) and 2380 (Fig. 4b), and the tip enriched degrees of freedom are always ordered last as given in Eq. (14). Hence, it can be seen that the terms associated with the tip degrees of freedom are the only one that are active, i.e. populate this matrix. As expected, the number of nonzero terms in the matrix is 2304, for both matrices, since only one tip element (four nodes) is enriched no matter how fine the discretization is. 4. Yau’s mutual potential energy for mixed mode extraction To extract the mixed mode SERRs we employ Yau’s mutual potential energy representation [27], which is based on Betti’s reciprocal theorem. This method employs the modes I and II analytical solution fields, in the vicinity of the crack tip, to decompose the SERRs into components. Yau’s method was further extended in [49]. Following [13] the mutual potential energy for equilibrium states (1) and (2) may be written in the following way:

1 2

ð2Þ

Pð1;2Þ ¼ ðd ÞT Kd

ð1Þ

1 ð1Þ ð2Þ ð2Þ þ ðd ÞT Kd  ðd ÞT f 2

ð1Þ

ð1Þ

 ðd ÞT f

ð2Þ

ð30Þ

The mutual energy release rate is therefore obtained by

Mð1;2Þ ¼ 



T @ Pð1;2Þ ð2Þ ð1Þ ›Kd ¼ d @a

ð31Þ

Note that ›K is defined in Eq. (8) and (12), and is nonzero only in the enriched tip elements. The mutual potential energy release rate can also be represented as

h i ð1Þ ð2Þ ð1Þ ð2Þ Mð1;2Þ ¼ 2b K I K I þ K II K II

ð32Þ

1m2

where b ¼ E for plain strain and b ¼ 1E for plane stress. The first equilibrium state (1) is the state solved by the finite element method and d(1) corresponds to the discrete solution. To extract the mixed mode stress intensity factors, two auxiliary states (2a) and (2b) corresponding to pure modes I and II are defined

(a)

0

(b)

0

100

200

Number of Eq.

Number of Eq.

500

nonzero terms in the stiffness derivative matrix (nz=2304)

300

400

1000

nonzero terms in the stiffness derivative matrix (nz=2304)

1500

500 2000 600 0

100

200

300

400

Number of Eq.

500

600

0

500

1000

1500

2000

Number of Eq.

Fig. 4. Structure and nonzero terms of the stiffness derivative matrix for the mixed mode problem in Section 5.2: (a) discretization leading to 628 equations and (b) discretization leading to 2380 equations.

H. Waisman / Engineering Fracture Mechanics 77 (2010) 3204–3215

pure mode I : pure mode II :

ð2aÞ

KI

ð2bÞ KI

ð2aÞ

¼ 1 and K II ¼ 0 and

ð2bÞ K II

¼0

3211

ð33Þ

¼1

For more details the readers are referred to reference [13]. 5. Numerical examples The proposed method is studied on pure modes I and II, and a mixed mode benchmark problems. The Strain Energy Release Rates results are converted to stress intensity factors and compared to the results obtained by the interaction integral [28,7,29]. In the comparison we evaluate the stiffness derivative in an analytical manner and hence the virtual crack extension is not required. Nonetheless, the semi-analytical derivation given by Eqs. (18) and (19) is also employed to illustrate the sensitivity of this approach to the selection of virtual crack extensions, which makes the fully analytical method much more robust and practical. We plot the rates of convergence for decreasing h (element size). In the current examples only nine elements contribute to the stiffness derivative computation, since only the tip element is enriched with the branch functions. Future work will consider enrichment of more layers around the tip element. For all problem considered we have used the trapezoidal integration rule for the tip element with 30  30 equally spaced integration points. Hence integration points do not change due to a virtual extension in the semi-analytical method. Similar type of integration strategy has been previously tested in [51,52]. 5.1. Pure modes I and II The first two examples show the convergence rates of the proposed approach (with the fully analytical derivation) for pure modes I and II. The problem consists of a square plate with dimensions 10  10 units, with an edge crack of 5 units in length, shown in Fig. 5. Displacement boundary conditions with normalized stress intensity factors KI = 1 and KII = 1 [53] are applied to generate modes I and II cracks openings, respectively. These essential boundary conditions (assuming a plane strain state) are given by

( Mode I :

  cosð2hÞ 2  2m  cos2 2h

  mÞ pffiffiffiffi uy ¼ 2ð1þ K I 2rp sinð2hÞ 2  2m  cos2 2h E

mÞ ux ¼ 2ð1þ KI E

pffiffiffiffi r 2p

ð34Þ

8

  pffiffiffiffi mÞ < ux ¼ 2ð1þ K II 2rp sin 2h 2  2m þ cos2 2h E h Mode II :  i pffiffiffiffi : uy ¼ 2ð1þmÞ K II r cos h 1  2m  sin2 h 2p 2 2 E

ð35Þ

where r is the distance from an enriched boundary node to the crack tip; h is the angle between the positive x axis and any boundary node. The boundary conditions are applied to all four sides of the rectangular domain. The Young’s modulus and Poisson ratio are chosen as E = 1  107 and m = 0.3, respectively. We employ the proposed stiffness derivative method to extract the Strain Energy Release Rates and then convert them into stress intensity factors. The proposed method is also compared to the interaction integral. The domain of influence and the active elements used to extract the SERRs are shown in Fig. 5. Note that in the stiffness derivative approach these elements are associated with tip enriched nodes in which the derivative ›K is not zero. Modes I and II deformed shapes obtained by XFEM are shown in Fig. 6a and b, respectively. The convergence results for Modes I and II are shown in Fig. 7a and b, respectively. It can be seen that

(a)

(b)

Fig. 5. Pure modes I and II edge crack example problem. Enriched nodes, boundary conditions and active elements are shown in both figures. (a) Active elements in the exact stiffness derivative method and (b) J-integral domain and contour.

3212

H. Waisman / Engineering Fracture Mechanics 77 (2010) 3204–3215

Fig. 6. Modes I and II deformed shapes: (a) mode I opening displacement and (b) mode II opening displacement.

(a)

KSDM I

1.6 1.4 1.2 1 0.8 0.6 0.4 0.2 0

0.9 KSDM II Jint

0.8

Jint

KI

relative error in KII [%]

I

relative error in K [%]

(b)

1.8

KII

0.7 0.6 0.5 0.4 0.3 0.2 0.1

0

0.2

0.4

0.6

0.8

1

1.2

0

1.4

0

0.2

0.4

0.6

h

0.8

1

1.2

1.4

h

Fig. 7. Convergence results for pure modes I and II. circles denote the SIFs obtained by the fully analytical stiffness derivative method and x denote the SIFs obtained from the interaction integral. (a) Convergence to mode I and (b) convergence to mode II.

while the interaction integral yields slightly more accurate results the convergence rates and overall behavior of both methods are quite similar. In the following Fig. 8a and b we study the sensitivity of the semi-analytical approach to a virtual crack extension, i.e. when computing the derivatives in Eqs. (18) and (19) for different increments Da. As it can be seen, the most accurate results are obtained when choosing an increment sizes of 106 6 Da 6 108 (depending on the density of the mesh). Hence the fully analytical approach is more attractive as it completely avoids the virtual extension philosophy.

(b)

0

10

9 × 9 elements 19 × 19 elements 39 × 39 elements 59 × 59 elements

−1

10

9 × 9 elements 19 × 19 elements 39 × 39 elements 59 × 59 elements

−2

10

II

relative error: Kfully −Ksemi II

−2

10

−4

10

I

relative error: Kfully−K

semi I

(a)

−6

10

−8

10

−10

10

−3

10

−4

10

−5

10

−6

10

−7

10

−8

10

−9

10

−10

10 −12

10

−11

−12

10

−10

10

−8

−6

10

10

Δa

−4

10

−2

10

10

−12

10

−10

10

−8

−6

10

10

−4

10

−2

10

Δa

Fig. 8. Sensitivity of the semi-analytical approach as function of different increments in Da. (a) Convergence to mode I and (b) convergence to mode II.

3213

H. Waisman / Engineering Fracture Mechanics 77 (2010) 3204–3215

τ

(a)

(b)

h a

w Fig. 9. Mixed mode edge crack example problem: (a) geometric definition of the problem with boundary conditions and (b) enriched jump and tip nodes.

(a)

SDM

KI

6

SDM

KII

Jint

5 4 3 2 1 0

4

KJint I

relative error in KII [%]

I

relative error in K [%]

(b)

7

0

0.1

0.2

0.3

0.4

h

0.5

0.6

0.7

KII

3.5

3

2.5

2

1.5

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

h

Fig. 10. Convergence results for mixed mode problem. circles denote the SIFs obtained by the fully analytical stiffness derivative method and x denote the SIFs obtained from the interaction integral. (a) Convergence to mode I and (b) convergence to mode II.

5.2. Mixed mode problem We study the performance of the proposed stiffness derivative method on a benchmark mixed mode problem [54,55]. The problem considered consists of a rectangular domain with dimensions 7  16 units subjected to a 1 unit stress shear traction on its upper edge and is fixed on the bottom, as shown in Fig. 9a. The crack length is 3.5 units and the enriched jump and tip nodes are illustrated in Fig. 9b. We employ the proposed stiffness derivative method to compute the stress intensity factors KI and KII, and compare the results to those obtained by the interaction integral. The analytical solution for this problem, given in [54,55], is:

K I ¼ 34:00;

K II ¼ 4:55

ð36Þ

Fig. 10a and b show the convergence studies of both methods to the SIFs KI and KII, respectively. Yau’s mutual potential energy (see Section 4) is used to decompose the SERRs into components and extract KI and KII for the proposed stiffness derivative method. It can be observed the results from both methods are comparable. 6. Conclusions A new approach to extract the Strain Energy Release Rates based on the extended finite element method (XFEM) is proposed. The idea is to take advantages of the tip functions in XFEM to compute the stiffness derivative in a closed form. Hence the forward difference approach of the classical Parks method [8], which introduces mesh perturbations due to a virtual crack extension, is no longer needed. The paper also investigates a semi-analytical approach which employs numerical derivative calculations (as opposed to the close form derivations). However sensitivity analysis shows that the semi-analytical approach is sub optimal as it depends on the increment Da taken, and hence may lead to larger errors. The implementation

3214

H. Waisman / Engineering Fracture Mechanics 77 (2010) 3204–3215

of the method is straight forward and is done on the fly during the analysis process. On the benchmark problems presented, the method is shown to perform close to the J-integral method. Future work will investigate how different enrichment schemes (more enriched layers around the tip) effect the accuracy of the method. Acknowledgements The financial support of the National Science Foundation under Grant CMMI-0856161 and the Department of Energy under Grant DE-SC0002137 is gratefully acknowledged. The author is also thankful to the outstanding suggestions made by reviewer #2, contributing to the current state of the paper. Appendix A The set of derivatives given in Section 3, Eqs. (16) and (17) are detailed in this Appendix. The crack-tip functions given in Eq. (2) are explicitly given by

8 F1 > > > F3 > > : F4

pffiffiffi ¼ r sin 2h pffiffiffi ¼ r cos 2h pffiffiffi ¼ r sin 2h cos h pffiffiffi ¼ r cos 2h cos h

ð37Þ

Their first derivatives with respect to r and h are given as

8 1ffi > f1r ¼ @F@r1 ¼ 2p sin 2h > r > > > 1ffi < f r ¼ @F 2 ¼ p cos 2h 2 @r 2 r 1ffi > f3r ¼ @F@r3 ¼ 2p sin 2h cos h > r > > > : f r ¼ @F 4 ¼ p 1ffi cos 2h cos h 4 @r 2 r

8 h f1 > > > > > f3h > > : h f4

¼ @F@h1 ¼ ¼ ¼ ¼

@F 2 @h @F 3 @h @F 4 @h

¼ ¼ ¼

ð38Þ

pffi r 2

cos 2h pffi  2r sin 2h pffi pffiffiffi r cos 2h cos h  r sin 2h sin h 2 pffi pffiffiffi  2r sin 2h cos h  r cos 2h sin h

ð39Þ

Using the relations given in Eqs. (21) and (22) and plugging the above derivatives into Eq. (16), one may obtain the analytical form for @Fj/@a. Next, to obtain the analytical form for @ 2Fj/@x@a in Eq. (17), one needs to obtain the second and mixed derivatives of the crack-tip functions (these derivatives are typically not used in standard extended finite element formulations). These are given by

8 @f r 2 1 > ¼ @@rF21 > @r > > > r > < @f2 ¼ @ 2 F22 @r @r @f3r @2 F3 > > ¼ > @r @r 2 > > > : @f4r @ 2 F 4 ¼ @r @r 2

3

¼  14 r 2 sin 2h 3

¼  14 r 2 cos 2h 3

¼  14 r 2 sin 2h cos h

ð40Þ

3

¼  14 r 2 cos 2h cos h

8 h pffi 2 @f1 > ¼ @@hF21 ¼  4r sin 2h > @h > > > pffi h > > < @f2 ¼ @ 2 F22 ¼  r cos h 4 2 @h @h pffi  @f3h > @2 F3 r 5 h > ¼ ¼  sin cos h þ 2 cos 2h sin h > 2 2 2 2 @h > @h > > pffi  > : @f4h @ 2 F 4 ¼ @h2 ¼  2r 52 cos 2h cos h  2 sin 2h sin h @h

ð41Þ

8 @f r @f h @2 F1 1ffi h 1 1 > p > > @h ¼ @r ¼ @h@r ¼ 4 r cos 2 > > r h > 2 > 1ffi < @f2 ¼ @f2 ¼ @ F 2 ¼  p sin 2h @h @r @h@r 4 r r h  2 @f @f > 3 ¼ 3 ¼ @ F3 ¼ p 1ffi 1 > cos 2h cos h  sin 2h sin h > @h @r @h@r 2 r 2 > > > >  : @f4r @f4h @ 2 F 4 1ffi 1 ¼ @r ¼ @h@r ¼  2p sin 2h cos h þ cos 2h sin h @h r 2

ð42Þ

Plugging the above derivatives into Eqs. (27) and (28) will give the analytical form of the derivatives @ 2Fj/@x@a in Eq. (17).

H. Waisman / Engineering Fracture Mechanics 77 (2010) 3204–3215

3215

References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] [25] [26] [27] [28] [29] [30] [31] [32] [33] [34] [35] [36] [37] [38] [39] [40] [41] [42] [43] [44] [45] [46] [47] [48] [49] [50] [51] [52] [53] [54] [55]

Chan SK, Tuba IS, Wilson WK. On the finite element method in linear fracture mechanics. Engng Fract Mech 1970;2:1–17. Shih CF, Delorenzi HG, German MD. Crack extension modeling with singular quadratic isoparametric elements. Int J Fract 1976;1:647–51. Banks-Sills L, Sherman D. Comparison of methods for calculating stress intensity factors with quarter-point elements. Int J Fract 1986;32:127–40. Li FZ, Shih CF, Needleman A. A comparison of methods for calculating energy release rates. Engng Fract Mech 1985;21(2):405–21. Rice JR. A path independent integral and the approximate analysis of strain concentration by notches and cracks. J Appl Mech 1968;35:379–86. Destuynder P, Lescure S, Djaoua M. Some remarks on elastic fracture mechanics. J Mécaniq Théoriq Appliq 1983;2(1):113–35. Moran B, Shih CF. A general treatment of crack tip contour integrals. Int J Fract 1987;35:295–310. Parks DM. A stiffness derivative finite element technique for determination of crack tip stress intensity factors. Int J Fract 1974;10(4):487–502. Hellen TK. On the method of virtual crack extensions. Int J Numer Methods Engng 1975;9(1):187–207. Rybicki EF, Kanninen MF. A finite element calculation of stress intensity factors by a modified crack closure integral. Engng Fract Mech 1977;9:931–8. Krueger R. The virtual crack closure technique: history, approach and applications. Technical report, NASA. NASA/CR-2002-211628, ICASE Report No. 2002-10; 2002. Banks-Sills L, Sherman D. On the computation of stress intensity factors for three-dimensional geometries by means of stiffness derivative and jintegral methods. Int J Fract 1992;53:1–20. Lin SC, Abel J. Variational approach for a new direct-integration form of the virtual crack extension method. Int J Fract 1988;38:217–35. Hwang CG, Wawrzynek PA, Tayebi AK, Ingraffea AR. On the virtual crack extension method for calculation of the rates of energy release rate. Engng Fract Mech 1998;59(4):521–42. Hwang CG, Wawrzynek PA, Ingraffea AR. On the calculation of derivatives of stress intensity factors for multiple cracks. Engng Fract Mech 2005;72:1171–96. Giner E, Fuenmayor FJ, Besa AJ, Tur M. An implementation of the stiffness derivative method as discrete analytical sensitivity analysis and its application to mixed mode in LEFM. Engng Fract Mech 2002;69:2051–71. Henshell RD, Shaw KG. Crack tip finite elements are unnecessary. Int J Numer Methods Engng 1975;9:495–507. Barsoum RS. One the use of isoparametric finite elements in linear fracture mechanics. Int J Numer Methods Engng 1976;10:25–37. Moës N, Dolbow J, Belytschko T. A finite element method for crack growth without remeshing. Int J Numer Methods Engng 1999;46:131–50. Belytschko T, Black T. Elastic crack growth in finite elements with minimal remeshing. Int J Numer Methods Engng 1999;45:601–20. Sukumar N, L Chopp D, Moës N, Belytschko T. Modeling holes and inclusions by level sets in the extended finite element method. Comput Methods Appl Mech Engng 2000;190(46–47):6183–200. Strouboulis T, Babuska I, Copps K. The design and analysis of the generalized finite element method. Comput Methods Appl Mech Engng 2000;181(1– 3):43–69. Strouboulis T, Babuska I, Copps K. The generalized finite element method. Comput Methods Appl Mech Engng 2001;190(32-33):4081–193. Fish J, Yuan Z. Multiscale enrichment based on the partition of unity. Int J Numer Methods Engng 2005;62(10):1341–59. Gracie R, Wang HW, Belytschko T. Blending in the extended finite element method by discontinuous Galerkin and assumed strain methods. Int J Numer Methods Engng 2008;74(11):1645–69. Waisman H, Belytschko T. Parametric enrichment adaptivity by the extended finite element method. Int J Numer Methods Engng 2008;73:1671–92. F Yau J, Wang SS, Corten HT. A mixed-mode crack analysis of isotropic solids using conservation-laws of elasticity. J Appl Mech 1980;47(2):335–41. Shih CF, Moran B, Nakamura T. Energy release rate along a three-dimensional crack front in a thermally stressed body. Int J Fract 1986;30(2):79–102. Nakamura T. Three-dimensional stress fields of elastic interface cracks. Trans ASME, J Appl Mech 1991;58(4):939–46. Moës N, Dolbow J, Belytschko T. A finite element method for crack growth without remeshing. Int J Numer Methods Engng 1999;46:131–50. Belytschko T, Black T. Elastic crack growth in finite elements with minimal remeshing. Int J Numer Methods Engng 1999;45:601–20. Sukumar N, L Chopp D, Moës N, Belytschko T. Modeling holes and inclusions by level sets in the extended finite element method. Comput Methods Appl Mech Engng 2000;190(46–47):6183–200. Belytschko T, Gracie R, Ventura G. A review of the extended/generalized finite element methods for material modelling. Modell Simulat Mater Sci Engng 2009;17(043001):24. Duarte CA, Babuska I, Oden JT. Generalized finite element methods for three-dimensional structural mechanics problems. Comput Struct 2000;77(2):215–32. Duarte CA, Hamzeh ON, Liszka TJ, Tworzydlo WW. A generalized finite element method for the simulation of three-dimensional dynamic crack propagation. Comput Methods Appl Mech Engng 2001;190(15–17):2227–62. Hughes TJR. The finite element method: linear static and dynamic finite element analysis. Prentice-Hall, Inc.; 1987. Williams ML. On the stress distribution at the base of a stationary crack. ASME J Appl Mech 1957;24:111–4. Osher S, Sethian JA. Fronts propagating with curvature-dependent speed: algorithms based on Hamilton–Jacobi formulations. J Comput Phys 1988;79:12–49. Sethian JA. Level set methods and fast marching methods: evolving interfaces in computational geometry. Fluid mechanics, computer vision, and materials science. Cambridge University Press; 1999. 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 Engng 2002;53:2549–68. Hwang CG, Ingraffea AR. Virtual crack extension method for calculating the second order derivatives of energy release rates for multiply cracked systems. Engng Fract Mech 2007;74(9):1468–87. Alvarado MA, Shah SP, John R. Mode-I fracture in concrete using center-cracked plate specimens. ASCE J Engng Mech 1989;115(2):366–83. Bazant ZP, Ohtsubo H, Aho K. Stability and post-critical growth of a system of cooling or shrinkage cracks. Int J Fract 1979;15(5):443–56. Nemat-Nasser S, Keer LM, Parihar KS. Unstable growth of thermally induced interacting cracks in brittle solids. Int J Solids Struct 1978;14:409–30. Gao H, Rice JR. A first-order perturbation analysis of crack trapping by arrays of obstacles. J Appl Mech 1989;56:828–36. Rao BN, Rahman S. Probabilistic fracture mechanics by Galerkin meshless methods. Part I: rates of stress intensity factors. Comput Mech 2002;28(5):351–64. Rahman S, Rao BN. Probabilistic fracture mechanics by Galerkin meshless methods. Part II: reliability analysis. Computat Mech 2002;28(5):365–74. Bazˇant ZP. Scaling of quasibrittle fracture: asymptotic analysis. Int J Fract 1997;83:19–40. Haber RB, Koh HM. Explicit expressions for energy release rates using virtual crack extensions. Int J Numer Methods Engng 1985;21:301–15. Barbero EJ, Reddy JN. The Jacobian derivative method for three-dimensional fracture mechanics. Commun Appl Numer Methods 1990;6:507–18. Elguedj T, Gravouil A, Combescure A. Appropriate extended functions for X-FEM simulation of plastic fracture mechanics. Computer Methods Appl Mech Engng 2006;195(7-8):501–15. Prabel B, Combescure A, Gravouil A, Marie S. Level set X-FEM non-matching meshes: application to dynamic crack propagation in elastic–plastic media. Int J Numer Methods Engng 2007;68(8):1553–69. Ewalds HL, Wanhill RJH. Fracture mechanics. New York: Edward Arnold; 1985. Wilson WK. Combined mode fracture mechanics. PhD thesis, University of Pittsburgh; 1969. Liu XY, Xiao QZ, Karihaloo BL. XFEM for direct evaluation of mixed mode SIFs in homogeneous and bi-materials. Int J Numer Methods Engng 2004;59:1103–18.