Virtual crack extension method for calculating the second order derivatives of energy release rates for multiply cracked systems

Engineering Fracture Mechanics 74 (2007) 1468–1487 www.elsevier.com/locate/engfracmech

Virtual crack extension method for calculating the second order derivatives of energy release rates for multiply cracked systems C.G. Hwang a

a,*

, A.R. Ingraffea

b

Department of Ubiquitous Technology in Infrastrure, Seoul University of Venture and Information, 1603-54 Seocho-Dong #413, Seocho-Gu, Seoul 137-070, Republic of Korea b 643 Frank Rhodes Hall, Cornell Fracture Group, Cornell University, Ithaca, NY 14853, USA Received 29 November 2005; received in revised form 23 July 2006; accepted 17 August 2006 Available online 20 October 2006

Abstract In this paper, we further generalize the work of Lin and Abel [Lin SC, Abel JF. Variational approach for a new directintegration form of the virtual crack extension method. Int J Fract 1988;38:217–35.] to the case of higher order derivatives of energy release rates for two-dimensional, multiply cracked systems. The direct integral expressions are presented for the energy release rates and their first and second order derivatives. The salient feature of this numerical method is that the energy release rates and their first and second order derivatives can be computed in a single analysis. It is demonstrated through a set of examples that the proposed method gives expectedly decreasing, but acceptably accurate results for the energy release rates and their first and second order derivatives. The computed errors were approximately 0.5% for the energy release rates, 3–5% for their first order derivatives and 10–20% for their second order derivatives for the mesh densities used in the examples. Potential applications of the present method include a universal size effect model and a probabilistic fracture analysis of cracked structures. Ó 2006 Elsevier Ltd. All rights reserved. Keywords: Virtual crack extension method; Second order derivative of energy release rates; Universal size effect model; Probabilistic fracture mechanics analysis

1. Introduction The derivatives of energy release rate in cracked structures provide useful information for the prediction of stability and arrest of a single crack [1], the growth pattern analysis of a system of interacting cracks [2,3], configurational stability analysis of evolving cracks [4], probabilistic fracture mechanics analysis [5] and universal size effect model [6]. In the case of multiple crack systems, for example, the variation of energy release rate at one crack tip due to the growth of any other crack must be calculated to determine the strength of the *

Corresponding author. Tel.: +82 2 3470 5281; fax: +82 2 523 6767. E-mail addresses: [email protected], [email protected] (C.G. Hwang).

0013-7944/$ - see front matter Ó 2006 Elsevier Ltd. All rights reserved. doi:10.1016/j.engfracmech.2006.08.009

C.G. Hwang, A.R. Ingraffea / Engineering Fracture Mechanics 74 (2007) 1468–1487

1469

Nomenclature ai B cf D E fx(x) f Gi G0 Gf GIc g g0 g00 J K KI k0 L N PF P S S0 Sb W u X ~e D j h m l r0 r22 rN x

length of crack i strain-nodal displacement matrix effective length of fracture process zone elastic constitutive matrix Young’s modulus joint probability density function applied nodal force vector energy release rate at crack i normalized energy release rate fracture energy mode-I fracture energy dimensionless energy release rate the first order derivative of dimensionless energy release rate the second order derivative of dimensionless energy release rate jacobian structural stiffness matrix mode-I stress intensity factor stress intensity factor length of plate standard shape function probability of failure potential energy of a cracked body structure size constant related to structure size S constant related to structure size S width of plate nodal displacement vector N-dimensional random vector with crack geometry virtual strain-like matrix geometry change of the elements due to virtual crack extension (3  4m) for plane strain, (3  m)/(1 + m) for plane stress angle of inclined radial crack Poisson’s ratio shear modulus remote tensile stress y-directional tensile stress nominal strength constant

interaction. In probabilistic fracture mechanics analysis of linear-elastic cracked structures, the first and second order reliability methods (FORM/SORM) require accurate estimates of energy release rates, their first and second order derivatives. Another use of the higher order derivatives is for size effect models that relate nominal strength to the structure size. The universal size effect model by Bazant [6] requires the first and second order derivatives of energy release rate. There have been several numerical methods for calculating the first order derivative of energy release rate. To the authors’ best knowledge, however, there is no method able to calculate the second order derivatives of energy release rates in the multiple crack systems. Sumi and his co-workers [7] presented a combined analytical and finite-element solution method in which special singular crack-tip elements were used. Those singular

1470

C.G. Hwang, A.R. Ingraffea / Engineering Fracture Mechanics 74 (2007) 1468–1487

elements are generated as a result of the second differentiation of potential energy and produce the higher order singularity at the crack-tip. De Koning and Lof [8] calculated the first order derivatives of crack-opening displacement and stress intensity factor due to virtual crack front perturbations for a three-dimensional (3-D) planar crack of arbitrary shape, based on finite elements and finite difference approximation. Meade and Keer [9] analyzed a half-plane crack with a slightly wavy crack front subjected to mode-I loading using an asymptotic expansion. They presented the first order variation of stress intensity factor due to variation in crack geometry. Rice [10] used 3-D weight function theory to develop a linear perturbation scheme for calculating the first order variation in stress intensity factor and crack opening displacement due to small changes in 3-D planar crack geometry. Nguyen and his co-workers [11] introduced an explicit expression for the matrix of the second order derivatives of energy with respect to the crack lengths in terms of path-independent integrals. Suo and Combescure [12] presented a double virtual crack extension method for the calculation of the second variation of energy for interacting two-dimensional (2-D) linear cracks. Chen and his co-workers [13] developed a numerical method for the first order derivatives of individual stress intensity factors for a 2-D single crack, using the concepts of continuum mechanics, M-integral and direct differentiation. Lin and Abel [14] introduced an analytical virtual crack extension method that uses variational theory in the finite element formulation. The method provides the direct integral forms of energy release rate and its higher order derivatives for a structure containing a single 2-D crack. The authors generalized Lin and Abel’s approach to a general case of a system of multiple cracks under arbitrary loading conditions in 2-D [15,16] and 3-D [17]. In this paper, the analytical virtual crack extension method is further extended for the calculation of the second order derivatives of energy release rates in a 2-D multiply cracked structure. In Section 2, the general formulation for the second order derivatives of energy release rates in the multiple crack systems is presented. In Section 3, several example problems are solved to demonstrate the accuracy of the present method. Section 4 discusses the potential application to the universal size effect model and probabilistic fracture analysis of linear-elastic cracked structures. 2. Formulation In this section, the analytical expressions for derivatives of energy release rates for a multiply cracked system modeled using finite element method are derived. The general formulation, first shown in [15], is first y

1

1

1

1

1

1

x

1

1

Fraction of δ a 0

0.75

1 Crack-tip node First ring of quarter-point nodes Second ring of nodes

Fig. 1. The uniform rosette of standard, quarter-point elements and mesh perturbation of the first ring of crack-tip singular elements.

C.G. Hwang, A.R. Ingraffea / Engineering Fracture Mechanics 74 (2007) 1468–1487

1471

repeated here for convenience. Next, the new formulation for the second order derivatives of energy release rates in multiply cracked systems is presented. For all the developments reported herein, it is assumed that each crack tip will be surrounded by a symmetric rosette of standard, quarter-point singular elements, as shown in Fig. 1, which can be augmented by surrounding rings of non-singular elements as shown in Fig. 2. 2.1. General formulation for the second order derivatives of energy release rates for multiple crack systems The potential energy P of a cracked body with multiple cracks is given by 1 P ¼ uT Ku  uT f ; 2

ð1Þ

where u, K and f are the nodal displacement vector, the structural stiffness matrix and the applied nodal force vector, respectively. The energy release rate at crack tip i can be expressed as Gi ¼ 

dP 1 dK df ¼  uT u þ uT ; dai 2 dai dai

ð2Þ

where ai is the length of crack i, and nonzero contributions to dK/dai and df/dai occur only over elements adjacent to the crack tip. It is noted that whenever crack-face, thermal and body force loadings are applied, the variations of loading must be taken into account to reflect the local load change on the crack-face or in the crack tip vicinity as a result of virtual crack extension. The variation of Gi, Eq. (2), with respect to the growth of any other crack, j, is,

Fig. 2. The uniform rosette of standard, quarter-point elements and mesh perturbation of the first and second rings of elements.

1472

C.G. Hwang, A.R. Ingraffea / Engineering Fracture Mechanics 74 (2007) 1468–1487

dGi dK du 1 T d2 K duT df d2 f ¼ uT  u uþ þ uT : dai daj 2 dai daj dai daj daj daj dai

ð3Þ

The second variation of Gi, Eq. (2), with respect to the growth of any other cracks, j and k, is, d2 Gi duT dK du d2 K du dK d2 u 1 duT d2 K 1 d3 K 1 ¼  uT  uT  u  uT u  uT dai dak daj dai daj dak 2 dak dai daj 2 dai daj dak 2 daj dak dak dai daj 

d2 K du d2 uT df duT d2 f duT d2 f d3 f þ þ þ þ uT : dai daj dak daj dak dai daj dai dak dak dai daj dai daj dak

ð4Þ

The variation of the displacement can be obtained from the variation of the global equilibrium equation Ku = f with respect to aj,   dK du df du df dK uþK ¼ or ¼ K 1  u : ð5Þ daj daj daj daj daj daj The second variation of the displacement is  2  d2 u df d2 K dK du dK du 1 ¼K  u  : daj dak daj dak daj dak daj dak dak daj

ð6Þ

We assume that the elements influenced by each crack-tip in a multiply cracked body comprise disjoint sets. Therefore, if i 5 j, then the second order variations of stiffness and loading with respect to two different crack extensions, ai and aj, vanish, d2 K d2 f ¼ ¼ 0: dai daj dai daj

ð7Þ

Hence, when i 5 j, the first order derivative of energy release rate with respect to two different crack extensions ai and aj is given as dGi dK du duT df ¼ uT þ : dai daj daj dai daj

ð8Þ

For i = j, dGi dK du 1 T d2 K duT df d2 f ¼ uT  u uþ þ uT 2 : 2 dai dai 2 dai dai dai dai dai

ð9Þ

When i 5 j 5 k, the second and third order variations of stiffness and loads are null, d2 K d2 K d3 K d2 f d2 f d3 f ¼ ¼ ¼ ¼ ¼ ¼ 0; dai daj dai dak dai daj dak dai daj dai dak dai daj dak and the second order derivatives of displacements are   d2 u dK du dK du ¼ K 1   : daj dak daj dak dak daj

ð10Þ

ð11Þ

Therefore, the second order derivatives of energy release rates are given by d2 Gi duT dK du dK d2 u d2 uT df ¼  uT þ : dai daj dak daj daK dai daj dak dak dai daj

ð12Þ

In the same manner, for i = j 5 k, d2 K d3 K d2 f d3 f ¼ ¼ ¼ ¼ 0; dai dak dai daj dak dai dak dai daj dak

ð13Þ

C.G. Hwang, A.R. Ingraffea / Engineering Fracture Mechanics 74 (2007) 1468–1487

1473

  d2 u dK du dK du 1 ¼K   ; dai dak dai dak dak dai

ð14Þ

d2 Gi duT dK du dK d2 u d2 K du d2 uT df duT d2 f ¼  uT  uT 2 þ þ : dai dai dak dai dak dai dak dai daK da2i dai dak dai dai dak

ð15Þ

For i = k 5 j, d2 K d3 K d2 f d3 f ¼ ¼ ¼ ¼ 0; dai daj dai daj dak dai daj dai daj dak   d2 u dK du dK du 1 ¼K   ; dai daj dai daj daj dai

ð16Þ ð17Þ

d2 G i duT dK du d2 K du dK d2 u d2 uT df duT d2 f ¼  uT 2  uT þ þ : dai daj dai dai daj dai daj dai daj da2i dai daj dai dai daj

ð18Þ

For i 5 j = k, d2 K d2 K d3 K d2 f d2 f d3 f ¼ ¼ ¼ ¼ ¼ ¼ 0; dai daj dai dak dai daj dak dai daj dai dak dai daj dak ! 2 d2 u d2 K dK du 1 d f ¼K  2 u2 ; da2j da2j daj daj daj d2 Gi duT dK du dK d2 u d2 uT df ¼  uT þ 2 : 2 dai da2j daj daj dai daj daj dai

ð19Þ ð20Þ ð21Þ

For i = j = k,  2  d2 u d2 K dK du 1 d f ¼K  2 u2 ; dai dai da2i da2i dai

ð22Þ

d2 Gi 1 d3 K d2 K du duT dK du dK d2 u d2 uT df duT d2 f d3 f ¼  uT 3 u  2uT 2   uT þ 2 þ2 þ uT 3 : 2 2 2 2 dai dai dai dai dai dai dai dai dai dai dai dai dai dai

ð23Þ

Element stiffness variations dke/da, d2ke/da2 and d3Ke/da3 are assembled to produce the global stiffness variations dK/da, d2K/da2, and d3K/da3. The element stiffness variations are given by, Z  T  dk e ¼ dB DB þ BT DdB þ Trð~eÞBT DB dV ; ð24Þ m

d2 k e ¼

Z



 d2 BT DB þ 2dBT DdB þ BT Dd2 B dV þ

m

d3 k e ¼

Z

Z

  2j~ejBT DB þ 2Trð~eÞðdBT DB þ BT DdBÞ dV ;

ð25Þ

m



 d3 BT DB þ 3d2 BT DdB þ 3dBT Dd2 B þ BT Dd3 B dV m Z  2 T  þ3 d B DB þ 2dBT DdB þ BT Dd2 B Trð~eÞ dV Z m   þ 2j~ejBT DB þ 6j~ejðdBT DB þ BT DdBÞ þ 2j~ejTrð~eÞBT DB dV ;

ð26Þ

m

where ~e is the virtual strain–like matrix, B is the strain-nodal displacement matrix, and D the elastic constitutive matrix. ~e is defined as, 9 8 oN > > >  = < 1 >  ~e11 ~e12 on 1 1 2 ~e ¼ J ; ð27Þ D Dn ¼ > oN > n ~e21 ~e22 > ; : 2;> on

1474

C.G. Hwang, A.R. Ingraffea / Engineering Fracture Mechanics 74 (2007) 1468–1487

where D’s are the geometry changes of the elements due to virtual crack extension and N is the standard shape function. The described formulation has been implemented in the fracture analysis code, FRANC2D, a 2-D finite-element-based code for simulating crack propagation [18]. This code performs automatic crack propagation simulation under LEFM conditions. 3. Numerical examples 3.1. Example 1. A center cracked infinite plate subjected to a uniform remote tensile stress: verification for the single crack case The first numerical example investigates a central crack in a large plate subjected to uniform remote tensile stress, r0. As shown in Fig. 3, the initial crack length to width ratio a/W is 0.01 to approximate a central crack in the infinite plate under a plane stress condition. For a direct comparison with analytical solutions for a single crack, only one-half of the plate was modeled, using the symmetry in the problem, with about 500 linear strain triangular elements including quarter-point elements at the crack-tip (Fig. 4a and b). It should be noted that df/da and d2f/da2 terms are null for a virtual crack extension under this loading. This model is analyzed to demonstrate the capability of the proposed method for evaluating the second derivative of mode-I energy release rate and stress intensity factor. The exact KI, dKI/da and dK 2I =da2 solutions for a mode-I crack under uniform remote tensile stress, r0, in an infinite plate can be expressed analytically as [19], pffiffiffiffiffiffi K I ¼ r0 pa;

dK I r0 ¼ da 2

rffiffiffi p ; a

d2 K I r0 ¼ da2 4a

rffiffiffi p : a

ð28Þ

Tables 1 and 2 show that the best computed solutions differ from the exact by about 0.1% for KI, 2–3% for dKI/da, and 5–10% for dK 2I =da2 , respectively, for the particular value of crack-tip element size used. Table 2 also shows that the solution accuracy for dK 2I =da2 is affected considerably by the number of rings of elements surrounding the crack tip that are involved in the mesh perturbation due to the virtual crack extension. When additional rings of surrounding elements from the crack tip are used in the perturbation, increasingly more accurate solutions for dK 2I =da2 are obtained. This observation also holds for computed values of KI and dKI/da.

σ =1.0

2H=20 2a

2W=20

Fig. 3. A central crack in a simulated infinite plate under remote unit tensile stress. Not to scale: initial a/W = 0.01. Example 1.

C.G. Hwang, A.R. Ingraffea / Engineering Fracture Mechanics 74 (2007) 1468–1487

1475

Fig. 4. (a) Finite element mesh for a central crack, Example 1, 2H = 2W = 20. (b) The detail of mesh near the crack tip, Example 1.

Table 1 Comparison of computed with exact solutions for KI and dKI a

KI exact

KI computed (error %)

dKI exact

dKI computed (error %)

0.10 0.11 0.12 0.13 0.14

0.5605 0.5879 0.6140 0.6391 0.6632

0.5610 0.5884 0.6145 0.6395 0.6627

2.802 2.672 2.558 2.458 2.369

2.739 2.617 2.561 2.404 2.313

(0.09) (0.08) (0.08) (0.06) (0.08)

(2.25) (2.06) (1.17) (2.20) (2.36)

Example 1, r0 = 1, 2W = 2H = 20, crack-tip element size = a/8.

The effect of crack tip element size on solution accuracy was also investigated. That is, the crack tip element size was decreased by repeatedly dividing the crack-tip element at the mid-point to create a new, regular 8-noded element, and a singular element. As the crack tip element size decreases, the radii of the second and third rings are also decreased. In Tables 3–5, it is shown that, as the crack tip element size is decreased, the solutions for dKI/da and dK 2I =da2 deteriorate while the accuracy of the computed values of KI is retained, unless at least 2 rings of elements are involved in the mesh perturbation. Therefore, it is recommended based

1476

C.G. Hwang, A.R. Ingraffea / Engineering Fracture Mechanics 74 (2007) 1468–1487

Table 2 Comparison of computed with exact solutions for d2KI for three different mesh perturbations a

Exact

1st Ring (error %)

1st + 2nd Rings (error %)

1st + 2nd + 3rd (error %)

0.10 0.11 0.12 0.13 0.14

14.012 12.145 10.659 9.454 8.459

52.373 58.599 40.801 46.023 22.607

17.465 15.386 13.447 12.394 10.941

15.044 11.637 11.915 9.881 8.911

(273.8) (382.5) (282.8) (386.8) (167.3)

(24.6) (26.7) (26.2) (31.1) (29.3)

(7.37) (4.18) (11.8) (4.52) (5.34)

Example 1, r0 = 1, 2W = 2H = 20, crack-tip element size = a/8.

Table 3 The computed values of KI for various crack-tip element sizes in three different mesh perturbation schemes Element size

1st Ring (error %)

1st + 2nd Rings (error %)

1st + 2nd + 3rd (error %)

a/8 a/16 a/32 a/64 a/128

0.5625 0.5622 0.5620 0.5617 0.5616

0.5617 0.5613 0.5610 0.5608 0.5606

0.5610 0.5611 0.5607 0.5604 0.5602

(0.36) (0.30) (0.27) (0.21) (0.20)

(0.21) (0.14) (0.09) (0.05) (0.02)

(0.09) (0.11) (0.04) (0.02) (0.05)

Example 1, r0 = 1, 2W = 2H = 20 (exact value of KI = 0.5605 for a = 0.1).

Table 4 The computed values of dKI for various crack tip element sizes in three different mesh perturbation schemes Element size

1st Ring (error %)

1st + 2nd Rings (error %)

1st + 2nd + 3rd (error %)

a/8 a/16 a/32 a/64 a/128

2.698 2.618 2.501 2.250 1.788

2.745 2.744 2.736 2.732 2.732

2.739 2.741 2.747 2.747 2.765

(3.71) (6.57) (10.7) (19.7) (36.2)

(2.03) (2.07) (2.36) (2.50) (2.50)

(2.25) (2.18) (1.96) (1.96) (1.32)

Example 1, r0 = 1, 2W = 2H = 20 (exact value of dKI = 2.802 for a = 0.1).

Table 5 The computed values of d2KI for various crack tip element sizes in three different mesh perturbation schemes Element size

1st Ring (error %)

1st + 2nd Rings (error %)

1st + 2nd + 3rd (error %)

a/8 a/16 a/32 a/64 a/128

52.373 (273.8) 67.957 (384.9) 225.25 (1507.4) 856.94 (6015.3) 3404.7 (24196.7)

17.465 (24.6) 18.990 (35.5) 29.614 (111.3) 79.044 (464.1) 281.20 (1906.7)

15.044 14.967 14.264 13.151 17.647

(7.4) (6.8) (1.8) (6.2) (25.9)

Example 1, r0 = 1, 2W = 2H = 20 (exact value of d2KI = 14.013 for a = 0.1).

on these results that at least one ring of nonsingular elements for dKI/da and two rings for dK 2I =da2 should be used in the mesh perturbation, along with the first ring of crack tip elements, for a more accurate analysis. However, adding additional rings alone does not guarantee uniform convergence as shown, for example, in the last column of Table 5. These data indicate that if the crack tip elements get too small, insufficient information from the singular field, which this element only can reproduce, is available for an accurate solution. In summary, for this particular example, it appears that optimal accuracy obtains with singular element size of a/32, and with two additional rings of non-singular elements participating in mesh perturbation. It is expected, of course, that this optimal arrangement would not necessarily hold for other structural configurations, but it is an initial guideline.

C.G. Hwang, A.R. Ingraffea / Engineering Fracture Mechanics 74 (2007) 1468–1487

1477

3.2. Example 2. A center cracked infinite plate subjected to a uniform remote tensile stress: the full model of plate with two crack tips The second example considers a full model of center-cracked plate under the mode-I uniform tensile stress, where the central crack is modeled with both crack tips. Three cases are considered for crack lengths 2a = 0.6, 0.8 and 1.0. The height (2H) and width (2W) of the plate are taken to be 40, respectively. The cracked plate is modeled with about 560 quadratic displacement elements, including quarter-point elements at the crack-tip, as shown in Fig. 5. In this example, crack-tip element size is a/24 and three rings of elements are used around the crack-tip for the mesh perturbation. Note that the analytical solutions in Eq. (28) represent mode-I stress intensity factors and their derivatives due to self-similar, co-linear virtual crack extensions at both crack tips. In order to estimate the derivatives of mode-I stress intensity factor at one crack tip due to the extension of any other crack tip, a 2-by-2 derivative matrix of mode-I stress intensity factors is calculated. Hence, we adopt a finite difference approximation to estimate the first order derivatives of mode-I stress intensity factors for the purpose of comparison with the present numerical solutions. The finite difference calculations for the first order derivative of stress intensity factor require three independent finite element analyses: one for the original crack length, a second with a finite, small, co-linear extension of crack tip 1, and the last with a finite, small, co-linear extension of crack tip 2. The crack extension increment used in these analyses is 0.01. In contrast, only a single finite element analysis is required to create the 2-by-2 derivative matrix of stress intensity factors by using the present numerical method. Table 6 compares the computed values of KI for three crack lengths, 2a = 0.6, 0.8, 1.0, with the exact solutions given in Eq. (28). The results are in a good agreement with exact solutions, giving maximum errors of 0.3% for those crack lengths, respectively.

Fig. 5. Finite element discretization for an inclined crack in a simulated infinite plate, Example 2, 2H = 2W = 40.

Table 6 The computed values of stress intensity factors for crack in a finite plate for crack lengths 2a = 0.6, 0.8 and 1.0 Crack tip

2a = 0.6

2a = 0.8

2a = 1.0

1 2

0.9719 (0.1%) 0.9719 (0.1%)

1.1234 (0.2%) 1.1234 (0.2%)

1.2579 (0.3%) 1.2579 (0.3%)

Example 2, 2W = 2H = 40 (analytical solutions for infinite plate: KI = 0.9708, 1.1210, 1.2533 for 2a = 0.6, 0.8, 1.0).

1478

C.G. Hwang, A.R. Ingraffea / Engineering Fracture Mechanics 74 (2007) 1468–1487

Table 7 compares present numerical results with finite difference calculations for the first order derivatives of stress intensity factors. In these cases, the present numerical results are in good agreement with finite difference solutions, giving a maximum difference of 1.4%. Table 8 lists the present numerical solutions for the second order derivatives of energy release rates in comparison with the finite difference solutions using Eqs. (29) and, (30), 2 3 2 3 dðK I Þ1 dðK I Þ1 dðK I Þ1 dðK I Þ1  2  6 d ðK I Þi 1 6 da2 7 da2 7 6 da1 7ða1 þ da1 ; a2 Þ  1 6 da1 7ða1 ; a2 Þ; ¼ ð29Þ 4 5 4 da1 dðK I Þ2 dðK I Þ2 5 daj da1 FDM da1 dðK I Þ2 dðK I Þ2 da1 da2 da1 da2 2 3 2 3 dðK I Þ1 dðK I Þ1 dðK I Þ1 dðK I Þ1  2  6 d ðK I Þi 1 6 da2 7 da2 7 6 da1 7ða1 ; a2 þ da2 Þ  1 6 da1 7ða1 ; a2 Þ; ¼ ð30Þ da2 4 dðK I Þ dðK I Þ 5 da2 4 dðK I Þ dðK I Þ 5 daj da2 2

FDM

da1

2

2

da2

da1

2

da2

where 2-by-2 derivative matrices, d(KI)i/daj, used in these calculations are calculated by the present numerical method for three different crack geometries (a1,a2), (a1 + da1, a2) and (a1, a2 + da2). It is shown that the results are in good agreement with finite difference solutions with differences of 10–20%. Table 7 Comparison between the present numerical solution and a finite difference method solution for the first order derivatives of stress intensity factors Present solutions

FDM solutions

Difference (%)

2a = 0.6 d(KI)1/da1 d(KI)2/da1

d(KI)1/da2 d(KI)2/da2

0.814 0.815

0.814 0.814

0.805 0.814

0.814 0.807

1.1 0.1

0.0 0.9

2a = 0.8 d(KI)1/da1 d(KI)2/da1

d(KI)1/da2 d(KI)2/da2

0.698 0.708

0.708 0.699

0.708 0.708

0.708 0.709

1.4 0.0

0.0 1.4

2a = 1.0 d(KI)1/da1 d(KI)2/da1

d(KI)1/da2 d(KI)2/da2

0.633 0.638

0.638 0.633

0.634 0.640

0.634 0.641

0.2 0.3

0.6 1.2

Example 2, 2W = 2H = 40.

Table 8 Comparison between the present numerical solution and FDM solution for the second order derivatives of stress intensity factors Present solutions

FDM solutions

d2(KI)1/da2da1 d2(KI)2/da2da1 d2 ðK I Þ1 =da22 d2 ðK I Þ2 =da22

0.669 0.675 0.668 0.667

0.668 0.667 0.673 0.629

0.597 0.647 0.603 0.647

0.646 0.670 0.647 0.675

12.1 4.3 10.8 3.1

3.4 0.4 4.0 6.8

Crack length, 2a = 0.6 d2 ðK I Þ1 =da21 d2 ðK I Þ2 =da21 d2(KI)1/da1a2 d2(KI)2/da1a2

d2(KI)1/da2da1 d2(KI)2/da2da1 d2 ðK I Þ1 =da22 d2 ðK I Þ2 =da22

0.464 0.435 0.429 0.429

0.429 0.429 0.435 0.480

0.385 0.417 0.389 0.416

0.416 0.382 0.416 0.372

17.0 4.1 9.3 3.0

3.0 11.0 4.4 22.5

Crack length, 2a = 0.6 d2 ðK I Þ1 =da21 d2 ðK I Þ2 =da21 d2 ðK I Þ1 =da1 da2 d2(KI)2/da2da2

d2(KI)1/da2da1 d2(KI)2/da2da1 d2(KI)1/da22 d2 ðK I Þ2 =da22

0.283 0.303 0.301 0.301

0.301 0.301 0.303 0.281

0.305 0.296 0.305 0.296

0.292 0.269 0.292 0.268

7.2 2.4 1.3 1.7

3.1 11.9 3.8 4.9

Crack length, 2a = 0.6 d2 ðK I Þ1 =da21 d2 ðK I Þ2 =da21 d2(KI)1/da1a2 d2(KI)2/da1a2

Example 2, 2W = 2H = 40.

Difference (%)

C.G. Hwang, A.R. Ingraffea / Engineering Fracture Mechanics 74 (2007) 1468–1487

1479

To assess the accuracy of the present numerical method with the analytical solution for the first order derivative of stress intensity factors, the first row sum in a derivative matrix of stress intensity factors, d(KI)I/ da1 + d(KI)1/da2 is used. The row sum represents dKI/da in Eq. (28) that is the first order derivative of mode-I stress intensity factors due to unit virtual crack extension at both crack tips 1 and 2. To compare with the analytical solutions for d2KI in Eq. (28), the sum of the first four components in the second order derivatives of stress intensity factors in Table 8 is used, which can be written by 2 X 2 X d2 ðK 1 Þi d2 ðK 1 Þ1 d2 ðK 1 Þ1 d2 ðK 1 Þ2 d2 ðK 1 Þ2 ¼ þ þ þ : daj da1 da21 da2 da1 da21 da2 da1 i¼1 j¼1

ð31Þ

Table 9 shows a direct comparison between the present solutions and analytical solutions for the first and second order derivatives of stress intensity factors. A good agreement between the numerical and analytical solutions was obtained. 3.3. Example 3. Radial multiple cracks embedded in the infinite plate The third example problem is a set of radial cracks embedded in an infinite plate subjected to the normal stress ryy = 1, as shown in Fig. 6. Two cracks are located along two different local radial axes at distance d = 1 from the origin and the length of each crack is taken to be 2a = 2. The orientation of crack 1 is kept constant at 30° while the orientation of the second crack is changed from 45° to 90°. Fig. 7 shows a finite element discretization for the case of h1 = 30° and h2 = 60°. The numerical studies are presented for the normalized pffiffiffi energy release rates, and their derivatives, where l is a shear modulus, G0 ¼ k 0 p; ð1 þ jÞ=8l; k 0 ¼ r22 a and j is defined as j ¼ 3  4m for plane strain;

j¼

3m for plane stress: 1þm

ð32Þ

Table 9 Comparison of the present numerical solutions with the analytical infinite plate solutions Eq. (28) for the first and second order derivatives of stress intensity factors P2 P2 d2 ðK I Þi P2 dðK I Þ1 dðK I Þ1 dðK I Þ1 Crack length j¼1 daj ¼ da1 þ da2 i¼1 i¼1 daj da1 (2a) Present solution Analytical solution Present solution Analytical solution (2W = 2H = 40) (2W = 2H = infinite) (2W = 2H = 40) (2W = 2H = infinite) 0.6 0.8 1.0

1.628 (0.6%) 1.407 (0.4%) 1.271 (1.4%)

2.679 (0.7%) 1.796 (2.5%) 1.187 (5.3%)

1.618 1.401 1.253

2.697 1.752 1.253

Example 2, 2W = 2H = 40.

y

Crack 2 2a 3 d

1

4

2a

θ2 2

Crack 1

θ1 x

Fig. 6. Radial multiple cracks embedded in the infinite plate, d = 1, Example 3.

1480

C.G. Hwang, A.R. Ingraffea / Engineering Fracture Mechanics 74 (2007) 1468–1487

Fig. 7. Finite element discretization for radial multiple cracks embedded in the infinite plate (2H = 2W = 40), Example 3.

Table 10 Comparison between numerical and analytical solutions for the normalized energy release rate, Gi/G0, h1 = 30°, Example 3 Crack tip

1 2 3 4

Present method

Analytical solution [21]

h2 = 45°

h2 = 60°

h2 = 90°

h2 = 45°

h2 = 60°

h2 = 90°

0.246 0.805 0.541 0.083

0.539 0.802 0.292 0.076

0.799 0.769 0.007 0.016

0.255 0.796 0.540 0.085

0.536 0.795 0.289 0.076

0.787 0.759 0.003 0.012

Shear modulus and Poisson’s ratio are chosen to be 0.5 and 0, respectively, for this calculation. Therefore, k0 = 1 and G0 = p. In this example, crack-tip element size is a/6 and three rings of elements are used around the crack-tip for the mesh perturbation. Table 10 compares the normalized energy release rates calculated by the present method with analytical solutions [20,21]. The analytical solutions are for multiple cracks in an infinite non-homogeneous plate, but they contain solutions for homogeneous material, which can be used as reference solutions. Table 10 indicates that a good agreement between numerical and analytical solutions was obtained for normalized energy release rate, with about 2–3% differences for the mesh density used in this analysis. For the first order derivatives of energy release rates, finite difference solutions are adopted as reference values for the case of h1 = 30° and h2 = 60°. Table 11 compares the first order derivatives of normalized energy release rates by the present method and finite difference approximation. The crack extension increment dak used in the finite difference solution was 1% of the half crack length, a. A difference of 5% was observed between the present and finite difference solutions in major diagonal terms. We obtain finite difference solutions for the second order derivatives of energy release rates by  2   d Gi 1 dGi dGi ffi ða þ dak ÞVCE  ðaÞVCE ; ð33Þ daj dak FDM dak daj daj where 4-by-4 derivative matrices dGi/daj are calculated by the present numerical method for five different crack geometries (a1, a2, a3, a4), (a1 + da1, a2, a3, a4), (a1, a2, + da2, a3, a4), (a1, a2, a3 + da3, a4) and (a1, a2, a3, a4 + da4). When k = 1, for example, Eq. (33) can be rewritten by

C.G. Hwang, A.R. Ingraffea / Engineering Fracture Mechanics 74 (2007) 1468–1487

1481

Table 11 The calculated values of the first order derivatives of normalized energy release rates, dGi/(dajG0), by the present virtual crack extension method and the finite difference approximation, h1 = 30° and h2 = 60°, Example 3 Present numerical solution

1 2 3 4

Finite difference solution

1

2

3

4

Row sum

0.418 0.350 0.289 0.122

0.350 0.442 0.083 0.043

0.289 0.083 0.454 0.070

0.122 0.043 0.070 0.052

0.357 0.832 0.318 0.043

2 dG1 

d2 Gi daj da1

da1

6 dG2 1 6 6 da1 ¼ 6 dG da 16 3 FDM 4 da1



2

3

4

Row sum

0.420 0.350 0.277 0.115

0.347 0.462 0.089 0.041

0.290 0.086 0.455 0.070

0.124 0.041 0.073 0.051

0.353 0.857 0.340 0.035

2 dG1

7

da1

dG1 da4

dG2 dG2 da2 da3

dG2 7 da4 7

dG4 dG4 dG4 da1 da2 da3

1

3

dG1 dG1 da2 da3

dG3 dG3 da2 da3

1 2 3 4

6 dG2 1 6 6 da1 7ða þ da1 ;a2 ;a3 ;a4 ÞVCE  6 dG3 7 1 da1 6 dG3 da4 5 4 da1 dG4 da4

dG1 dG1 dG1 da2 da3 da4

3 7

dG2 dG2 dG2 7 da2 da3 da4 7

7ða1 ;a2 ;a3 ; a4 ÞVCE :

dG3 dG3 dG3 7 da2 da3 da4 5

dG4 dG4 dG4 dG4 da1 da2 da3 da4

ð34Þ It is shown from Table 12 that the present numerical solutions for the second order derivative of normalized energy release rate are in good agreement with finite difference solutions. 3.4. Example 4. A system of 3 inclined cracks in a finite plate under remote tensile stress The next example problem is a system of 3 inclined cracks (6 crack tips) in a finite plate under a remote unit tensile stress, Fig. 8. The finite element model used is shown in Fig. 9. The crack tip element size is taken to be a/6. A plane stress condition is assumed and Young’s modulus, E, is 1.0. In this example, crack-tip element size is a/6 and three rings of elements are used around the crack-tip for the mesh perturbation. Table 12 The calculated values of the second order derivatives of normalized energy release rates, dG2i =ðdaj dak G0 Þ, by the present numerical method and the finite difference approximation, h1 = 30° and h2 = 60°, Example 3 Present numerical method 1

2

Finite difference approximation 3

3

4

0.087 0.042 0.072 0.040

k ¼ 1 : dG2i =ðdaj da1 G0 ÞFDM 1 1.299 0.047 2 0.047 0.038 3 1.484 0.111 4 0.093 0.043

1.484 0.111 0.863 0.073

0.093 0.043 0.073 0.035

0.111 0.037 0.149 0.025

0.042 0.258 0.025 0.249

k ¼ 2 : dG2i =ðdaj da2 G0 ÞFDM 1 0.064 0.050 2 0.050 0.020 3 0.112 0.039 4 0.041 0.255

0.112 0.039 0.119 0.023

0.041 0.255 0.023 0.245

0.111 0.037 0.149 0.025

0.930 0.149 0.553 0.017

0.072 0.025 0.017 0.008

k ¼ 3 : dG2i =ðdaj da3 G0 ÞFDM 1 1.556 0.112 2 0.112 0.038 3 0.959 0.148 4 0.071 0.025

0.959 0.148 0.576 0.016

0.071 0.025 0.016 0.010

0.042 0.258 0.025 0.249

0.072 0.025 0.017 0.008

0.040 0.249 0.008 0.193

k ¼ 4 : dG2i =ðdaj da4 G0 ÞFDM 1 0.092 0.042 2 0.042 0.260 3 0.081 0.024 4 0.039 0.252

0.081 0.024 0.017 0.009

0.039 0.252 0.009 0.218

k ¼ 1 : dG2i =ðdaj da1 G0 ÞVCE 1 1.736 2 0.059 3 1.480 4 0.087

0.059 0.064 0.111 0.042

1.480 0.111 0.930 0.072

k ¼ 2 : dG2i =ðdaj da2 G0 ÞVCE 1 0.059 2 0.064 3 0.111 4 0.042

0.064 0.483 0.037 0.258

k ¼ 3 : dG2i =ðdaj da3 G0 ÞVCE 1 1.480 2 0.111 3 0.930 4 0.072 k ¼ 4 : dG2i =ðdaj da4 G0 ÞVCE 1 0.087 2 0.042 3 0.072 4 0.040

4

1

2

1482

C.G. Hwang, A.R. Ingraffea / Engineering Fracture Mechanics 74 (2007) 1468–1487

σ =1.0

Tip 2

2a

Tip 4

Tip 6

θ

2H

Tip 3

Tip 1

s

Tip 5

s

2W

Fig. 8. A series of inclined cracks in a finite plate under remote tensile stress (2W = 2H = 20, 2a = 2.0, S = 5.0, r = 1.0, h = 45°). Not to scale. Example 4.

Fig. 9. Finite element discretization for the series of inclined cracks in a plate, Example 4.

Table 13 lists the values of energy release rates computed by the present virtual crack extension method in comparison with other numerical methods, including J-integral method, displace correlation method and modified crack closure integral method. As showin in the table, the comparable results are obtained. Table 14 compares the first order derivatives of energy release rates calculated by the present method with finite difference solutions. The crack extension increment used in the finite difference solutions was one hundredth of the half crack length, a. Good agreement was obtained with 2–3% differences between major diagonal terms clustered along the diagonal of the derivative matrix.

C.G. Hwang, A.R. Ingraffea / Engineering Fracture Mechanics 74 (2007) 1468–1487

1483

Table 13 The calculated values of the energy release rates, Gi, by the present virtual crack extension method and other numerical methods, Example 4

1 2 3 4 5 6

Present VCE method (using 1st + 2nd + 3rd rings)

J-integral (using 1st ring)

Displacement correlation

Modified crack closure integral

2.674 2.829 2.943 2.943 2.828 2.674

2.682 2.837 2.951 2.951 2.836 2.681

2.679 2.839 2.952 2.952 2.839 2.679

2.631 2.782 2.894 2.894 2.781 2.632

(0.3%) (0.3%) (0.3%) (0.3%) (0.3%) (0.3%)

(0.2%) (0.4%) (0.3%) (0.3%) (0.4%) (0.2%)

(1.6%) (1.7%) (1.7%) (1.7%) (1.7%) (1.6%)

Table 14 The calculated derivatives of energy release rate, dGi/daj, by the present virtual crack extension method and a finite difference approximation, Example 4 Present numerical method

1 2 3 4 5 6

Finite difference approximation

1

2

3

4

5

6

Row sum

1.008 1.114 0.254 0.146 0.055 0.026

1.114 1.150 0.309 0.253 0.091 0.055

0.254 0.309 1.097 1.160 0.252 0.146

0.146 0.253 1.160 1.097 0.309 0.254

0.055 0.091 0.252 0.309 1.150 1.114

0.026 0.055 0.146 0.254 1.114 1.008

2.603 2.972 3.218 3.219 2.971 2.603

1 2 3 4 5 6

1

2

3

4

5

6

Row sum

1.037 1.123 0.254 0.146 0.055 0.026

1.123 1.179 0.311 0.253 0.092 0.056

0.254 0.311 1.130 1.169 0.253 0.148

0.146 0.253 1.170 1.127 0.311 0.255

0.054 0.091 0.253 0.310 1.182 1.124

0.025 0.055 0.147 0.254 1.125 1.037

2.639 3.012 3.265 3.259 3.018 2.646

The second order derivatives of energy release rates are calculated by the present method and compared with finite difference solutions for k = 1, 2, 3 in Table 15. The finite difference calculations have been done by Eq. (33), where the 6-by-6 derivative matrices dGi/daj are computed by the present numerical method

Table 15 The second order derivatives of energy release rates (Young’s modulus E = 1), dG2i =ðdaj dak Þ, Example 4 Present numerical method 1 k ¼ 1 : dG2i =ðdaj da1 ÞVCE 1 0.108 2 0.138 3 0.053 4 0.024 5 0.004 6 0.004

2

Finite difference approximation

3

4

5

6

1

2

3

4

5

6

0.025 0.092 0.089 0.012 0.011 0.013

0.005 0.032 0.021 0.011 0.035 0.018

0.004 0.018 0.013 0.013 0.018 0.004

0.138 0.140 0.124 0.092 0.031 0.018

0.053 0.124 0.162 0.089 0.020 0.013

0.024 0.092 0.089 0.012 0.011 0.013

0.004 0.031 0.020 0.011 0.035 0.018

0.004 0.018 0.013 0.013 0.018 0.004

k ¼ 1 : dG2i =ðdaj da1 ÞFDM 1 0.187 0.152 0.056 2 0.152 0.166 0.125 3 0.056 0.125 0.163 4 0.025 0.092 0.089 5 0.005 0.032 0.021 6 0.004 0.018 0.013

k ¼ 2 : dG2i =ðdaj da2 ÞVCE 1 0.138 0.140 2 0.140 0.044 3 0.124 0.052 4 0.092 0.167 5 0.031 0.051 6 0.018 0.035

0.124 0.052 0.050 0.123 0.021 0.011

0.092 0.167 0.123 0.039 0.021 0.020

0.031 0.051 0.021 0.021 0.051 0.031

0.018 0.035 0.011 0.020 0.031 0.004

k ¼ 2 : dG2i =ðdaj da2 ÞFDM 1 0.162 0.153 0.125 2 0.153 0.045 0.055 3 0.125 0.055 0.051 4 0.092 0.169 0.123 5 0.031 0.052 0.021 6 0.018 0.035 0.011

0.092 0.169 0.123 0.039 0.021 0.021

0.031 0.052 0.021 0.021 0.052 0.031

0.018 0.035 0.011 0.021 0.031 0.004

0.162 0.050 0.062 0.109 0.039 0.012

0.089 0.123 0.109 0.109 0.123 0.089

0.020 0.021 0.039 0.123 0.167 0.092

0.013 0.011 0.012 0.089 0.092 0.024

k ¼ 3 : dG2i =ðdaj da3 ÞFDM 1 0.053 0.125 0.165 2 0.125 0.052 0.055 3 0.165 0.055 0.021 4 0.090 0.124 0.124 5 0.021 0.021 0.041 6 0.013 0.011 0.014

0.090 0.124 0.124 0.136 0.124 0.090

0.021 0.021 0.041 0.124 0.168 0.092

0.013 0.011 0.014 0.090 0.092 0.024

k ¼ 3 : dG2i =ðdaj da3 ÞVCE 1 0.053 2 0.124 3 0.162 4 0.089 5 0.020 6 0.013

0.124 0.052 0.050 0.123 0.021 0.011

1484

C.G. Hwang, A.R. Ingraffea / Engineering Fracture Mechanics 74 (2007) 1468–1487

for four different crack geometries, (a1, a2, a3, a4, a5, a6), (a1,+da1, a2, a3, a4, a5, a6), (a1, a2,+da2, a3, a4, a5, a6), (a1, a2, a3,+da3, a4, a5, a6). It should be noted that agreement between two solutions is excellent, except for the terms with i = j = k. It is due to the fact that the third order variation of stiffness is involved in the calculations of the diagonal terms because the higher order variations possess more errors. The results for k = 4, 5, 6 are not presented for brevity, but the similar solution accuracy is obtained. It is noted that the errors produced in the computation of the second order derivatives of energy release rate appear to vary from 5% to 20% for the example problems in the above. It is partly because the comparison is made with respect to crude finite difference solutions that are significantly influenced by the mesh configurations used for the example problems and the choice of a ‘finite’ virtual crack extension. The comparison with exact solutions only can give a reliable estimate on the accuracy, but it is only available for a single crack problem in Example 1. For the multiple crack systems, finite difference solutions are used for comparison due to the lack of exact solutions. 4. Applications In this section, potential application problems of the present method are discussed. The problems include the universal size effect model [6,22] and probabilistic fracture mechanics analysis of linear-elastic cracked structures [5]. 4.1. Universal size effect model It is well known that fracture of quasi-brittle materials exhibits significant size effects. The size effect is the change of nominal strength rN as a function of structure size S in geometrically similar structures with similar cracks. While the energy dissipated in geometrically similar structures with similar cracks is proportional to the crack length a and to the structure size S, the energy release caused by fracture in structures under the same nominal stress grows with structure size S faster than proportionally – hence the size effect [6,22]. In the universal size effect model by Bazant [22], the nominal strength rN is given as a function of structure size S, fracture energy Gf, effective length of fracture process zone cf, dimensionless energy release rate g and its higher order derivatives, g 0 and g00 . The model represents a smooth transition from the case of plasticity, for which there is no size effect, to the case of LEFM, for which the size effect is the strongest possible. In the universal size effect model, the nominal strength of the structure, rN, is given by  ð1=2Þ S 2S b rN ¼ r0 þ1 ; ð35Þ S0 S where sffiffiffiffiffiffiffiffi EGf r0  ; cf g 0 cf g0 ; g hg00 i S b ¼ cf : 4g0 S0 ¼

ð36Þ ð37Þ ð38Þ

As shown in these equations the universal size effect model requires the dimensionless energy release and its higher order derivatives, where the appropriate notations are g ¼ gða0 Þ; dgðaÞ at a ¼ a0 ; g0 ¼ g0 ða0 Þ ¼ da d2 gðaÞ g00 ¼ g00 ða0 Þ ¼ at a ¼ a0 : da2

ð39Þ ð40Þ ð41Þ

C.G. Hwang, A.R. Ingraffea / Engineering Fracture Mechanics 74 (2007) 1468–1487

1485

Dimensionless current and initial crack lengths are defined as a¼

a S

and

a0 ¼

a0 : S

ð42Þ

In addition, the effective length of the fracture process zone is given by, cf ¼ cf for a0 P cf ; cf ¼ xcf for a0 ¼ 0;

ð43Þ ð44Þ

where x is a certain constant. As shown in Eqs. (35)–(41), the energy release rate and its higher order derivatives are indispensible for the determination of the nominal strength rN according to the universal size effect model, and those values can be accurately calculated by the present numerical method, as shown in Section 3. 4.2. First and second order reliability methods (FORM/SORM) for probabilistic fracture mechanics analysis This subsection shows that derivatives of energy release rates are also required for the probabilistic fracture-mechanics analysis of linear-elastic cracked structures. Rahman and Rao [5] presented a stochastic meshless method for probabilistic fracture-mechanics analysis. In this stochastic approach, the performance of the cracked structure is evaluated using the probability of failure, PF, defined as Z def P def ¼ P ½H ðX Þ < 0 ¼ fx ðxÞ dx; ð45Þ r F gðxÞ<0

where fx(x) is the joint probability density function of X. In a cracked structure, X denotes an N-dimensional random vector with crack geometry (a), crack orientation, mode-I fracture energy (GIc) at crack initiation, and loads. If the propagation of a mode-I crack constitutes a failure condition, then the performance function is H ðxÞ ¼ GIc  G;

ð46Þ

where H(x) = 0 means crack propagation. To compute the failure probability in Eq. (45), Rahman and Rao [5] used the FORM which leads to nonlinear constrained optimization, requiring the derivative of the performance function with respect to crack size (a) as oH oG ¼ ; oa oa

ð47Þ

if fracture energy GIc is a constant. If the SORM is used, then the probabilistic analysis requires the second order derivatives of energy release rate. Then, the second order derivative of the performance function will be given as o2 H o2 G ¼  : oa2 oa2

ð48Þ

If the method is extended to multiple crack systems, then the first and second order derivatives of the performance function must be evaluated using the derivative matrix of energy release rates by oH i oGi ¼ ; oaj oaj

ð49Þ

o2 H i o2 G i ¼ : oaj oak oaj oak

ð50Þ

By Eqs. (49) and (50), we can estimate the failure probability of the cracked structure due to an interaction between different cracks i and j, using FORM and SORM, respectively. It is notable that for mixed-mode crack propagation the probabilistic analysis requires the higher order derivatives of individual stress intensity factors instead of energy release rates [5,16]. Again, values of those derivatives can be accurately calculated by the present numerical method, as shown in Section 3.

1486

C.G. Hwang, A.R. Ingraffea / Engineering Fracture Mechanics 74 (2007) 1468–1487

5. Conclusions In this paper, the work of Lin and Abel [14] is further extended to the general case of multiple crack systems. Analytical expressions are presented for energy release rates, and their first and second order derivatives for a multiply cracked body. The salient feature of this approach is that all of these values can be computed in a single finite element analysis. The present method has been implemented and verified in FRANC2D, fracture analysis software developed by the Cornell Fracture Group. It is demonstrated through numerical examples that the proposed method gives expectedly decreasing, but acceptably accurate results for the energy release rates, their first and second order derivatives. The computed errors are about 0.5% for the energy release rates, 3–5% for their first order derivatives and 10–20% for the second order derivatives for the mesh density used in these examples. Sensitivity of the accuracy of the proposed method to crack tip element size and to the number of rings of surrounding elements participating in the local mesh perturbation was also evaluated. The solution accuracy of the present method is influenced by several factors: the number of rings of additional non-singular elements around the crack-tip used for the mesh perturbation, the size of the mesh perturbation zone, and the mesh configuration in the neighborhood of crack-tips. Potential application problems of the present method include probabilistic fracture mechanics analysis of linear-elastic cracked structures and the universal size effect model both of which require calculation of derivatives of energy release rates. Acknowledgements The first author acknowledges the financial support from the Kyo-Hak publishing company and the YongOk Foundation in Seoul, South Korea. The authors thank Drs. Shbeeb and Binienda for providing analytical solutions of Example 3. The second author acknowledges the support of the Northrop Grumman Corporation and the Cornell Theory Center. References [1] Alvarado, Shah SP, John R. Mode-I fracture in concrete using center-cracked plate specimens. J Engng Mech ASCE 1989;115(2):366–83. [2] 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. [3] 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. [4] Gao H, Rice JR. A first-order perturbation analysis of crack trapping by arrays of obstables. J Appl Mech 1989;56:828–36. [5] Rahman S, Rao BN. Probabilistic fracture mechanics by Galerkin meshless methods – part II: reliability analysis. Comput Mech 2002;28(5):365–74. [6] Bazant ZP. Scaling theories for quasi-brittle fracture: recent advances and new directions. In: Wittmann Folker H, editor. Fracture mechanics of concrete structures, Proceedings FRAMCOS-2. D-79014 Freiburg: AEDIFICATIO Publishers; 1995. p. 515–34. [7] Sumi Y, Nemat-Nasser S, Keer LM. A new combined analytical and finite-element solution method for stability analysis of the growth of interacting tension cracks in brittle solids. Int J Engng Sci 1980;18:211–24. [8] DeKoning AU, Lof CJ. K-distributions extrapolated on the basis of stress intensity rates. In: Luxmoore AR, Owen DRJ, editors. Proceedings of the third international conference on numerical methods in fracture mechanics. Swansea, UK: Pineridge Press; 1984. p. 195–203. [9] Meade KP, Keer LM. Stress intensity factors for semi-infinite plane crack with a wavy front. J Elasticity 1984;14:79–92. [10] Rice JR. First-order variation in elastic fields due to variation in location of a planar crack front. J Appl Mech 1985;52:571–9. [11] Nguyen QS, Stolz C, Debruyne G. Energy methods in fracture mechanics: stability, bifurcation and second variations. Eur J Mech A/ Solids 1990;9(2):157–73. [12] Suo X, Combescure A. Double virtual crack extension method for crack growth stability assessment. Int J Fract 1992;57:127–50. [13] Chen G, Rahman S, Park YH. Shape sensitivity analysis in mixed-mode fracture mechanics. Comput Mech 2001;27:282–91. [14] Lin SC, Abel JF. Variational approach for a new direct-integration form of the virtual crack extension method. Int J Fract 1988;38:217–35. [15] Hwang CG, Warwzynek PA, Tayebi AK, Ingraffea AR. On the virtual crack extension method for calculating rates of energy release rate. Engng Fract Mech 1998;59(4):521–42. [16] 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. [17] Hwang CG, Warwzynek PA, Ingraffea AR. On the virtual crack extension method for calculating the derivatives of energy release rates for a three-dimensional planar crack of arbitrary shape under mode-I loading. Engng Fract Mech 2001;68:925–47.

C.G. Hwang, A.R. Ingraffea / Engineering Fracture Mechanics 74 (2007) 1468–1487

1487

[18] Wawrzynek PA, Ingraffea AR. Interactive finite element analysis of fracture processes: an integrated approach. Theoret Appl Fract Mech 1987;8:137–50. [19] Paris PC, Sih GC. Stress analysis of cracks. ASTM-STP 381; 1965. p. 65. [20] Shbeeb NI, Binienda WK, Kreider KL. Analysis of the driving forces for multiple cracks in an infinite nonhomogeneous plate, Part I: theoretical analysis. ASME J Appl Mech 1999;66:492–500. [21] Shbeeb NI, Binienda WK, Kreider KL. Analysis of the driving forces for multiple cracks in an infinite nonhomogeneous plate, Part II: numerical solutions. ASME J Appl Mech 1999;66:501–6. [22] Bazant ZP. Scaling of quasibrittle fracture: asymptotic analysis. Int J Fract 1997;83:19–40.