Evaluation of J1 and J2 integrals for curved cracks using an enriched boundary element method

Engineering Fracture Mechanics 78 (2011) 623–637

Contents lists available at ScienceDirect

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

Evaluation of J1 and J2 integrals for curved cracks using an enriched boundary element method R. Simpson ⇑, J. Trevelyan ⇑ School of Engineering & Computing Sciences, Durham University, South Road, Durham DH1 3LE, UK

a r t i c l e

i n f o

Article history: Received 29 July 2009 Received in revised form 1 December 2010 Accepted 18 December 2010 Available online 25 December 2010 Keywords: BEM Fracture J-integral J1 J2 enrichment

a b s t r a c t This paper introduces an enriched Boundary Element Method in which functions are introduced that are known to model singularities or discontinuities from a priori knowledge of the solution space. Additional fundamental solutions are introduced to solve for the additional unknowns created by enrichment and a numerical integration routine is outlined for the evaluation of strongly singular and hypersingular enriched boundary integrals. The solution of a curved crack in an infinite domain by Muskhelishvili is used to assess the accuracy of the method. Using an appropriate technique to evaluate J1 and J2 integrals, it is found that very good agreement with the exact solution is seen with improvements in accuracy over similar FEM implementations. Ó 2011 Elsevier Ltd. All rights reserved.

1. Introduction The design of engineering structures is often governed by the presence of cracks that can, under certain conditions, propagate rapidly leading to catastrophic failure. The goal of engineers is to prevent such failures by understanding the mechanisms that control the initiation and growth of cracks. By then incorporating this knowledge within computational methods, it is then possible to offer practising engineers efficient and accurate tools that are beneficial in the design process. In the context of linear elastic fracture mechanics (LEFM) it is generally accepted that the stress intensity factors (SIFs), components which give a measure of the singularity of a crack in defined orientations, are the key parameters that govern crack growth. Much research has been carried out to develop computational methods that provide accurate SIFs. It was discovered early on [1] that the presence of a singularity affects convergence significantly. Since then a great number of methods have been developed to overcome the difficulties imposed by the presence of a crack tip with varying degrees of success. In some early studies using the Finite Element Method (FEM) detailed mesh refinements were common. Such approaches have been augmented by recent advances in remeshing that appear in Miranda et al. [2] and Bouchard et al. [3]. Further improvements were made by incorporating analytical expressions known to capture the asymptotic nature of a crack within elements surrounding the crack tip where SIFs were output directly [4], but a significant development was the introduction of quarter-point elements that provided a much more readily implemented solution. The work can be attributed to Henshell and Shaw [5] and Barsoum [6] who, independently, realised that by moving the midnodes of elements to quarter-point posipffiffiffiffi tions the desired q (where q is the distance to the crack tip) variation in displacements could be achieved. However, certain limitations were shown to be inherent with the use of these elements; Harrop [7] described the reasons for size dependence and Ingraffea and Manu [8] showed that for 3D cases the results were dependent on Poisson’s ratio. Much more recently, the Partition of Unity Method (PUM) [9] has been used extensively in conjunction with the FEM to model fracture problems with ⇑ Corresponding authors. Tel.: +44 1913342487 (R. Simpson). E-mail addresses: [email protected] (R. Simpson), [email protected] (J. Trevelyan). 0013-7944/$ - see front matter Ó 2011 Elsevier Ltd. All rights reserved. doi:10.1016/j.engfracmech.2010.12.006

624

R. Simpson, J. Trevelyan / Engineering Fracture Mechanics 78 (2011) 623–637

Nomenclature Cij Dkij e na e na Ena kij ; E kijl ; F kij F, F1, F2 I1I ; I2I ; I1II ; I2II Jn J1, J2 KI, KII e I; K e II K M Na ni na e na Pna ij ; P ijl ; Q ij R Skij Tij tj t j Uij uj uaj uj na e na V na Ij ; V Ijl ; W Ij W x 0 x (x, y) dij

C C0 Cc c Ccþ int ; Cint

j

K

l m wuIj ; wuIIj (q, h)

rij n np

scalar function of surface geometry derivative fundamental solution residual error in auxiliary boundary integral equation boundary integrals terms used for singular integral evaluation boundary integrals for auxiliary equations Jacobian of coordinate transformation on element n J-integrals mode I and II stress intensity factors amplitudes of the enrichment functions number of nodes on the element shape function associated with node a component, in direction given by index j, of outward pointing unit normal boundary integrals radius of region around crack tip for J2 integral evaluation derivative fundamental solution traction fundamental solution traction component in direction given by index j auxiliary traction component in direction given by index j displacement fundamental solution displacement component in direction given by index j coefficient describing displacement at node a in direction given by index j auxiliary displacement component in direction given by index j boundary integrals strain energy ‘field’ point on boundary collocation point Cartesian coordinates Kronecker delta domain boundary portion of J-integral contour within the domain crack boundary portions of upper and lower crack surfaces used in J-integral Kosolov constant, function of Poisson’s ratio strain energy jump term at crack tip shear modulus Poisson’s ratio enrichment functions for displacement in direction given by index j polar coordinate system with origin at crack tip ij-component of the Cauchy stress tensor parametric coordinate defining element parametric coordinate value at collocation point

great success. The landmark paper by Moës et al. [10] introduced the Extended Finite Element Method (XFEM) that not only implements enrichment of nodes by using the PUM with a basis that captures the singular field, but represents cracks independently of the mesh preventing the need for remeshing under crack growth. The method has since been extended to 3D where use has been made of level sets to represent the crack surface aiding in crack propagation and SIF evaluation [11]. Bordas and Duflot [12–14] presented various a posteriori error estimators for XFEM. Other attempts to improve fracture modelling include the hybrid element of Karihaloo and Xiao [15], and the development of meshless methods, which also automatically overcome the problem of remeshing in problems involving crack growth (since there is no mesh). Fracture mechanics computations have been considered using meshfree methods from the early papers on such methods [16]; Nguyen et al. [17] provide a useful recent review and include a discussion of the practical use of the methods for cracked bodies. Meshfree algorithms for fracture mechanics remain a subject of considerable research activity, with recent work focussing on locally enriched approximations for problems containing material and geometric non-linearity [18,19]. Liew et al. [20] extended enriched meshless schemes to a meshless boundary integral formulation, using similar enrichment to that of the current paper. Although not enjoying as much attention as the FEM, the Boundary Element Method (BEM) has seen success when applied to fracture problems. The inherent advantage of a reduction in dimensionality when using the BEM is often

R. Simpson, J. Trevelyan / Engineering Fracture Mechanics 78 (2011) 623–637

625

countered by the argument that BEM matrices are fully populated and unsymmetric (in comparison to the symmetric and sparsely populated FEM matrices). But in terms of fracture, there are some distinct features of the BEM that make it a competitive numerical tool for engineering analysis. For crack problems we are mostly concerned with boundary parameters making the BEM particularly suitable, and for the simulation of crack propagation, where normally with a domain discretised method complete remeshing would be required, simply including a small number of additional boundary elements for each crack advance is all that is needed. Of course, with the use of quadratic interpolation for most BEM implementations, the method suffers from the same problems encountered by the FEM due to singularities at crack tips. Special crack tip elements have been used where Yamada et al. [21] developed special shape functions to interpolate displacements. In addition, due to the independent representation of tractions in the BEM, special traction crack-tip shape functions were developed by Tanaka and Itoh [22] to model the singularity of O(q1/2). But, as described earlier for the FEM, a simple technique is to use quarterpoint elements for those elements positioned adjacent to the crack tip. Implementation of quarter-point elements in the BEM has shown that their practical use is restricted to 2D examples with flat cracks since their construction is greatly complicated for curved cracks and, as described by Aliabadi [23] the geometry restrictions for quarter-point elements in 3D are such that meshing becomes awkward. Further developments in the BEM include work by Watson [24] who included enrichment by incorporating the Williams expansion within elements surrounding the crack. To account for the additional unknowns introduced by the method, auxiliary boundary integral equations were used requiring the use of rather complicated procedures to evaluate singular integrals, including an additional requirement for Hermitian elements. Perhaps the most popular method at present to model fracture problems using boundary integral equations (BIEs) is the Dual Boundary Element Method (DBEM) developed by Portela et al. [25]. This uses an additional BIE to allow collocation at coincident nodes on crack surfaces overcoming the problem of a singular system when applying the conventional BEM to fracture problems. In addition, to allow the accurate evaluation of SIFs, a circular J-integral path centred at the crack tip can be used with an appropriate routine to allow decomposition into separate fracture modes. The DBEM has become established also in the commercial software industry (BEASY [26]). The J-integral (more precisely the J1-integral) [27] is perhaps the most popular method within computational fracture mechanics to evaluate SIFs. It precludes the need to evaluate explicitly singular terms seen at a crack tip by calculating a path independent integral surrounding the crack tip that can be related to the strain energy release rate and, if required, stress intensity factors. For flat cracks with tractions equal to zero on the crack faces the implementation is straightforward. However, for curved cracks the J1-integral only becomes path independent with the inclusion of an additional boundary term and, to allow a full description of the crack-tip behaviour, it is necessary to evaluate the J2-integral. This involves the evaluation of a singular strain energy term that becomes problematic, especially in the FEM. Various methods have been proposed to overcome this problem; an excellent overview of the various methods applied to XFEM is given by Legrain et al. [28]. In the current work curved cracks using the BEM are considered and evaluation of J1 and J2-integrals is given. Enrichment is applied by incorporating the analytical solution for crack-tip displacements given by Williams [29] and using the DBEM to allow modelling of coincident crack surfaces. Then, by employing the technique first proposed by Eischen [30] and later implemented for curved cracks using the FEM [31] accurate evaluation of J1 and J2 integrals is achieved. This paper is organised in the following fashion. In Section 2 we describe our strategy for enrichment of the boundary elements local to a crack tip. In Section 3 we present some notes on implementation while in Section 4 the procedure for evaluating the J1 and J2 integrals is outlined. Finally, an example of a curved crack in an infinite domain is given where results are compared against an exact solution. 2. Enrichment of the BEM Use of standard quadratic interpolation of displacements in the BEM is insufficient at capturing the crack tip singularity unless refined meshes are used. Instead, by incorporating the analytical expression for displacements around a crack tip within the approximation, more accurate results are to be expected. This can be achieved by taking the first order terms of the Williams expansion allowing the following expression for crack tip displacements to be written:

uj ¼ K I wuIj ðq; hÞ þ K II wuIIj ðq; hÞ

ð1Þ

where KI and KII are defined as the mode I and mode II stress intensity factors, and the terms wIj(q, h) and wIIj(q, h) are given by the following functions obtained from the Williams solution

wuIx ¼

1 2l

rffiffiffiffiffiffiffi

q

2

cosðh=2Þ½j  1 þ 2 sin ðh=2Þ 2p rffiffiffiffiffiffiffi

1 q sinðh=2Þ½j þ 1 þ 2 cos2 ðh=2Þ 2l 2p rffiffiffiffiffiffiffi 1 q wuIy ¼ sinðh=2Þ½j þ 1  2 cos2 ðh=2Þ 2l 2p rffiffiffiffiffiffiffi 1 q 2 wuIIy ¼  cosðh=2Þ½j  1  2 sin ðh=2Þ 2l 2p wuIIx ¼

ð2aÞ ð2bÞ ð2cÞ ð2dÞ

626

R. Simpson, J. Trevelyan / Engineering Fracture Mechanics 78 (2011) 623–637

Here, q and h define the polar coordinates centred at the crack tip, such that h = ±p on the upper and lower crack surfaces, l is the material shear modulus and j is the parameter defined as

8 < 3  4m plane strain j¼ 3m : plane stress 1þm

ð3Þ

These expressions can then be used to approximate displacements around a crack by interpolating in a similar fashion to that of Benzley [4]

e II wu þ uj ¼ f K I wuIj þ K IIj

M X

Na uaj

ð4Þ

a¼1

where uaj represents a nodal coefficient (but usually represents nodal displacements for conventional interpolation), Na is the conventional Lagrangian shape function for node a and M is the number of nodes on the element. The last term of (4) is included to allow represention of rigid body motion of the crack tip while the first and second terms are used to represent the e I and K e II no longer represent SIFs but are coefficients that instead, when they are combined with the singular crack-tip field. K functions of (2) and the conventional interpolation (given by the last term of (4)) return enriched crack-tip displacements. To return the SIFs directly, shape functions would be required to interpolate the functions of (2) while the value of the enrichment functions at nodal points would need to be subtracted. It would be entirely possible to include this within the enriched formulation, but in the present work the expression (4) is used throughout. Since the expression for displacements is now proposed, it is possible to include enrichment within the BEM formulation by substituting Eq. (4) into the displacement boundary integral equation (DBIE). To aid understanding, this is first stated in its unenriched form (where boundary collocation is used) as

Z Z C ij ðx0 Þuj ðx0 Þ þ -- T ij ðx0 ; xÞuj ðxÞdCðxÞ ¼ U ij ðx0 ; xÞt j ðxÞdCðxÞ; C

i; j ¼ x; y

ð5Þ

C 0

where Tij and Uij are the traction and displacement fundamental solutions, x and x are the source and field points that lie on 0 the surface C, uj and tj are displacements and tractions and Cijuj represents a jump term due to the singular integral at x . The R Einstein summation convention is assumed throughout. The integral -- denotes that the integral is taken in the Cauchy Principal Value sense. For numerical implementation, Eq. (5) is then discretised to give

C ij ðx0 Þuj ðx0 Þ þ

Ne X M X n¼1 a¼1

na Pna ij uj ¼

Ne X M X

na Q na ij t j

ð6Þ

n¼1 a¼1

where

Pna ij ¼ Q na ij ¼

Z

1

1 Z 1

Na ðnÞT ij ½x0 ; xðnÞJ n ðnÞdn

ð7aÞ

N a ðnÞU ij ½x0 ; xðnÞJ n ðnÞdn

ð7bÞ

1

Ne is the number of elements and Jn(n) is the Jacobian of the transformation (x, y) ? n for element n. If Eq. (4) is now substituted into (5) the enriched form of the boundary integral is given by Ne X Ne X Ne X M M M   X X X na na e na K el ¼ e l wu ðx0 Þ þ P C ij ðx0 Þ uj ðx0 Þ þ K P na Q na lj ij uj þ ij t j ; ijl n¼1 a¼1

n¼1 a¼1

l ¼ I; II

ð8Þ

n¼1 a¼1

e na is expressed as The additional term P ijl

e na ¼ P ijl

Z

1

1

Na ðnÞT ij ½x0 ; xðnÞwulj ðnÞJ n ðnÞdn

ð9Þ

e na ¼ 0 otherwise. However, the problem of modelling coincident crack surfaces (which is a if element n is enriched and P ijl well-known problem for the conventional BEM) is present also in Eq. (8). To eliminate this problem the dual boundary element method is adopted to allow the modelling of any general crack problem. The method uses an independent BIE known as the traction boundary integral equation (TBIE) used for collocation on one of the crack surfaces (see Fig. 1) which 0 is derived by differentiating the DBIE with respect to the position of the source point x . The result is a BIE that contains 2 hypersingular terms of O(r ) in 2D which require special consideration for numerical implementation. The TBIE is stated here as

1 tj ðx0 Þ þ ni ðx0 Þ 2

Z Skij ðx0 ; xÞuk ðxÞdCðxÞ ¼ ni ðx0 Þ -- Dkij ðx0 ; xÞt k ðxÞdCðxÞ;

C

i; j; k ¼ x; y

ð10Þ

C

where Skij and Dkij represent fundamental solutions that are derived by differentiating Tij and Uij and the integral sents a Hadamard finite-part integral. In the same manner as before, the TBIE is discretised to give

repre-

R. Simpson, J. Trevelyan / Engineering Fracture Mechanics 78 (2011) 623–637

627

Fig. 1. DBEM crack modelling strategy. Ne X Ne X M M X X 1 na 0 na t j ðx0 Þ þ ni ðx0 Þ Ena F na kij uk ¼ ni ðx Þ kij t k 2 n¼1 a¼1 n¼1 a¼1

ð11Þ

where

Ena kij ¼ F na kij ¼

Z

1

1 Z 1

Na ðnÞSkij ½x0 ; xðnÞJ n ðnÞdn

ð12aÞ

Na ðnÞDkij ½x0 ; xðnÞJ n ðnÞdn

ð12bÞ

1

To apply enrichment to the TBIE, Eq. (4) is substituted into Eq. (11) to give Ne X Ne X Ne X M M M X X X 1 na 0 0 na e e E na t j ðx0 Þ þ ni ðx0 Þ Ena F na kijl K l ¼ ni ðx Þ kij uk þ ni ðx Þ kij t k 2 n¼1 a¼1 n¼1 a¼1 n¼1 a¼1

where

e E na kijl ¼

Z

1 1

Na ðnÞSkij ½x0 ; xðnÞwulk ðnÞJ n ðnÞdn

ð13Þ

ð14Þ

e na ¼ 0. Then, by applying the strategy of the DBEM outlined in [25] where the TBIE is if element n is enriched, otherwise E kijl applied for collocation on one of the crack surfaces while the DBIE is used for all other boundaries, it is possible, through the use of Eqs. (8) and (13) to apply enrichment to any general crack model. 3. Implementation 3.1. Additional fundamental solutions It is clear, after inspecting the enriched BIEs of (8) and (13) that enrichment introduces additional degrees of freedom (DOF) which, in the present method, correspond to each mode of fracture. Several methods can be used to overcome this problem such as additional BIEs derived by differentiating Eq. (5) and the use of additional collocation points (which the present authors have investigated). Instead, a much simpler scheme was used where the crack-tip displacements and tractions derived by Williams [29] can be used to formulate an independent BIE. This is illustrated by first considering Betti’s reciprocal theorem which, if body forces are ignored, can be expressed in terms of displacements and stresses as

Z C

tj uj dC ¼

Z C

t j uj dC

ð15Þ

where (uj, tj) and ðuj ; t j Þ correspond to two systems. The system ðuj ; tj Þ is taken as the crack-tip solution where displacements are given by Eq. (1) and tractions are expressed as

tj ¼ K I wtIj ðq; hÞ þ K II wtIIj ðq; hÞ

ð16Þ

wtlj ð

where the functions q; hÞ are constructed using the crack-tip stresses given by the Williams’ solution [29] and the relation ti = rijnj. If we arbitrarily choose KI = 1 and KII = 0 and substitute Eqs. (1) and (16) for u and t, then we arrive at the BIE

Z

C

wtIj ðq; hÞuj dC ¼

Z

C

wuIj ðq; hÞt j dC

ð17Þ

628

R. Simpson, J. Trevelyan / Engineering Fracture Mechanics 78 (2011) 623–637

This can be discretised as Ne X M X

na V na Ij uj ¼

n¼1 a¼1

Ne X M X

na W na Ij t j

ð18Þ

n¼1 a¼1

where

V na Ij ¼

Z

W na Ij ¼

1

Na ðnÞwtIj ðnÞJ n ðnÞdn

1 Z 1

1

ð19aÞ

Na ðnÞwuIj ðnÞJ n ðnÞdn

ð19bÞ

Fig. 2. Flow chart of enriched Dual BEM process.

R. Simpson, J. Trevelyan / Engineering Fracture Mechanics 78 (2011) 623–637

629

Fig. 3. Enrichment of crack-tip elements.

Enrichment can now be applied by substituting Eq. (4) into (18) to give Ne X M X

na V na Ij uj þ

n¼1 a¼1

Ne X M X

e na K el ¼ V Ijl

n¼1 a¼1

Ne X M X

na W na Ij t j

ð20Þ

n¼1 a¼1

e na is expressed as where the additional term V Ijl

e na ¼ V Ijl

Z

1 1

Na ðnÞwtIj ðnÞwulj ðnÞJ n ðnÞdn

ð21Þ

where it should be noted that I is not a tensorial subscript and the Einstein summation convention should not be applied to the repeated index j. To derive the second additional BIE required to form a square system, Eqs. (1) and (16) are used in the same manner as before with KI = 0 and KII = 1 leading to a simple change of replacing wtIj and wtIj with wtIIj and wtIIj in the integrals (19a), (19b) and (21). A flow chart showing the enriched Dual BEM process is presented in Fig. 2. 3.2. Numerical integration of strongly singular and hypersingular boundary integrals In any BEM implementation the evaluation of singular integrals is necessary. The use of rigid-body motion is often made to circumvent the problem. But, as outlined in [25], explicit singular integration is required for the DBEM implementation where, if flat elements are used to model the crack surfaces, simple analytical expressions can be used. The introduction of enrichment terms no longer permits the use of these equations but instead requires a general singular integration routine for strongly singular and hypersingular boundary integrals. Another need for explicit singular integration comes from the use of curved elements on crack surfaces since the analytical expressions given in [25] are no longer valid. The numerical integration method used in the present work is based on Guiggiani et al. [32] where the integrand is first expressed as a Taylor series allowing the singular component to be subtracted. This then leaves a regular integral that can be evaluated by conventional quadrature routines while the singular component can be found analytically. If the singular integrand is expressed as F(np, n) where np denotes the local coordinate of the source point, then the integral can be written as

I¼

Z

þ1

"

Fðnp ; nÞ 

1

F 2 ðnp Þ ðn  np Þ

þ 2

F 1 ðnp Þ n  np

!#

    np  1  þ F 2 ðnp Þ  1  1 dn þ F 1 ðnp Þ ln  np 2 ð1; 1Þ 1  np 1 þ np n p þ 1

ð22Þ

where F2(np) and F1(np) are regular functions determined by a Taylor series expansion of the integrand about the source point. Full details of these expressions can be found in Appendix A. The method is valid for both hypersingular and strongly singular integrals where it is found, for the case of strongly singular integrals of O(1/r), that the term F2(np) simply equals zero. 3.3. Enrichment strategy The current method allows the enrichment of multiple elements surrounding the crack tip therefore allowing the singular approximation to be extended to a user-defined area. This is a distinct advantage over singular crack-tip elements (such as quarter-point) which confine the singular region to the size of the crack tip element while in reality, the singularity dominated zone may extend beyond this. Fig. 3 illustrates enrichment of a crack where four elements (two for each crack surface) have enrichment applied. The current work included a study of the effect of varying the enrichment region by increasing the number of enriched elements and observing the accuracy of the results. 4. Evaluation of J1 and J2 integrals The J-integral, whose components are given by

Jk ¼

  Z  @ui dC; Wnk  rij nj @xk C

k ¼ 1; 2

ð23Þ

630

R. Simpson, J. Trevelyan / Engineering Fracture Mechanics 78 (2011) 623–637

Fig. 4. Definitions of integration paths used for J1 and J2 evaluation.

is often chosen to evaluate fracture paramaters where, in the case of flat cracks, the J1 component is path-independent and can be decomposed to produce individual SIFs [33]. For a curved crack the situation is complicated since the J1-integral is no longer path independent due to a non-zero contribution from the crack surfaces. It is also necessary to evaluate the J2-integral which presents difficulties due to a singular strain energy jump term. Many FEM implementations circumvent the problem of evaluating the singular boundary integral by transforming the term into a domain integral but the method used in the present implementation is based on work by Chang and Wu [31]. By specifying a radius R that is centred at the crack tip the J-integral components can be approximated by

Jk ¼

   Z @ui dC þ Wnk  rij nj Wnk dC þ 2Kdk2 R1=2 @xk C0 ðCcþ RÞþðCc RÞ int

Z

int

Fig. 5. Biaxially loaded curved crack in an infinite domain.

ð24Þ

R. Simpson, J. Trevelyan / Engineering Fracture Mechanics 78 (2011) 623–637

631

c where the integral paths C0 ; Ccþ int  R, Cint  R and the radius R are defined in Fig. 4. The variable denoted by K represents the singular strain energy jump term seen across the crack faces and d represents the Kronecker delta function. Then, by noting that J2 and K are invariant for a given problem, several values of R can be used to give a set of simultaneous equations which can be solved by a least-squares scheme to give J2 and K. In the case of J1 the formulation is simplified greatly by the omission of the third term in Eq. (24).

Fig. 6. Boundary element mesh used for curved crack analysis.

Fig. 7. J-integral integration paths (path 3 shown).

632

R. Simpson, J. Trevelyan / Engineering Fracture Mechanics 78 (2011) 623–637

5. Examples To illustrate the accuracy of the method the analytical solution of a curved crack in an infinite domain (see Fig. 5) given by Muskhelishvili was used where exact expressions for J1 and J2 are given [34]. Due to symmetry only half the problem was modelled where an example mesh of six elements per line is shown in Fig. 6. This is in keeping with the FEM implementation of [31]. In addition, to assess the accuracy gained by modelling the crack surface with curved elements in comparison with flat elements, tests were carried out with both element types. The J-integration routine outlined in Section 4 was implemented by using integration paths centred at the crack tip which progressively increased in size (see Fig. 7) while the radius R was chosen to be 0.01l, 0.015l, 0.02l, 0.025l and 0.03l to allow the evaluation of K in Eq. (24). To assess convergence the number of elements per line was increased while uniform element spacing was used throughout. It should be noted that grading was not used in any of the meshes. Discontinuous elements were used throughout, but it would be entirely feasible to use continuous elements on all non-crack surfaces using the strategy as illustrated by Portela et al. [25]. The geometry was defined with Rc = 0.9 mm, b = 45°, w/2 = 20 mm and L = 40 mm. The loadings were chosen as rx = ry = 100 MPa with material properties E = 210 GPa and m = 0.3.

Fig. 8. J1 convergence study results.

Fig. 9. J2 convergence study results.

Table 1 Comparison of J1 and J2 values, with % errors in parentheses.

Exact Unenriched (flat) Enriched (flat) QP elements (flat) Unenriched(curved) Enriched (curved) FEM [31]

J1 (Pa m)

J1 % error

J2 (Pa m)

J2 % error

0.06592 0.06719 0.06698 0.06708 0.06695 0.06632 0.06684

– 1.92 1.61 1.75 1.56 0.61 1.40

0.04661 0.04557 0.04595 0.04593 0.04749 0.04679 0.04716

– 2.24 1.43 1.46 1.89 0.39 1.18

R. Simpson, J. Trevelyan / Engineering Fracture Mechanics 78 (2011) 623–637

633

We note that Liu and Xu [35] present quarter-point BEM results for this problem for varying angles b, but direct comparison is difficult because details of the mesh were not provided. 6. Results As mentioned in Section 3.3 the effect of increasing the size of the enrichment region was investigated. In fact, what was found for all meshes was an increase in accuracy, over the unenriched case, when enrichment was applied only to the elements adjacent to the crack tip; any further increase in enrichment had an insignificant effect on the accuracy of the results. It was therefore decided to enrich only crack-tip elements for the convergence study. Similar results were shown for the

(a) J2 path independence for flat unenriched elements

(b) J2 path independence for flat enriched elements

(c) J2 path independence for QP elements

(d) J2 path independence for curved unenriched elements

(e) J2 path independence for enriched elements Fig. 10. J2 path independence.

634

R. Simpson, J. Trevelyan / Engineering Fracture Mechanics 78 (2011) 623–637

XFEM by Laborde et al. [36] and Béchet et al. [37], and Duflot and Bordas [14] have presented convergence results for different values of the XFEM enrichment radius. The results of the convergence study are shown in Figs. 8 and 9 while comparisons are made with the FEM solution of [31] in Table 1. There is a clear improvement in accuracy due to enrichment, and we note that the results for the unenriched formulations have still not converged to their final values over the range of model sizes considered. In fact, at convergence, the error is reduced by more than half in comparison to the unenriched analysis with curved elements. By comparing the relative errors of unenriched curved and flat elements it can be seen, especially for coarse meshes, that the accuracy suffers due to the approximation of a curved surface with straight lines. Much of the same is true for the J2 integrals where a clear improvement in accuracy is seen after introducing enrichment. This is even more significant than the improvement in J1 and a reduction in error by a factor of more than three is seen. Interestingly, by observing the unenriched element results for J2, an improvement in accuracy is realised for coarse meshes by using curved elements over flat elements. But at convergence, where the meshes are so fine that the flat elements very closely approximate a curve, the accuracies are very similar and, in fact, the flat elements give slightly better results.

(b) Comparison of uniform and graded mesh results for J2

(a) Graded mesh Fig. 11. Convergence in J2 for uniform and graded meshes.

Fig. 12. Residual error in auxiliary BIEs.

R. Simpson, J. Trevelyan / Engineering Fracture Mechanics 78 (2011) 623–637

635

It is noted from Table 1 that the enriched dual BEM method presented in the current paper outperforms the FEM results in [31] while using considerably fewer degrees of freedom. Using the integration paths defined in Fig. 7, J2 results are shown (see Fig. 10) for enriched, unenriched with curved elements, unenriched with flat elements and quarter-point elements for each of the paths. In general, path independence is seen in all results while a few observations can be made. For the analyses using curved elements the largest variation is seen at the lowest mesh density but this quickly reduces as the mesh density increases. In contrast, the analysis using flat elements shows a relatively uniform path-dependence as the mesh density is increased but, for meshes other than the coarsest one used for analysis, is consistently larger than that seen for curved elements. By observing path 2 in Fig. 10d and e it can be seen that the results for this path differ somewhat from those for all other paths at higher mesh densities. The reason for this can be explained by the numerical integration inaccuracies introduced by the close proximity of this path to the crack tip. Finally, a comparison between uniformly graded elements and elements graded towards the crack surfaces (see Fig. 11a for an example graded mesh) was made while uniformly graded elements were used on the crack surfaces. As can be seen in Fig. 11b for J2 values, a gain in accuracy is achieved using grading for coarse meshes while both meshes converge to the same result. Interestingly, for the graded meshes convergence is approached from below in contrast to uniform elements. The reason for this can only be speculated upon, and further study would be required. Readers may notice that the algorithm does not appear to converge to the exact Muskhelishvili solution. This is a result of the additional fundamental solutions used to provide the auxiliary equations as described in Section 3.1; these fundamental solutions are the Williams expressions for the flat crack. As a result there is a residual error in the boundary integral Eq. (17) when applied to curved cracks. Let us denote using I1I and I2I the integrals on the left and right hand sides (respectively) of (17), and similarly I1II and I2II the corresponding integrals in the BIE derived using the mode II solution as the fundamental solution. Let the residual error be defined as e ¼ ðI1N  I2N Þ; N ¼ I; II. The error term e is found to be small in comparison with the integrals I1N ; I2N even for the high curvature in the example studied; this is demonstrated by plotting the variation in e with the radius of curvature of the crack as shown in Fig. 12. The small value of e allows us still to obtain results markedly more accurate than the other formulations studied. One method of overcoming this small error, at the cost of slightly degraded matrix conditioning, is to use extra collocation points in preference to the mode I and mode II solutions to provide the auxiliary equations, as described in Simpson and Trevelyan [38]. 7. Conclusions This paper outlines a method for enriching the BEM and DBEM for analysis of general crack geometries where use is made of the Williams solution for crack-tip displacements and tractions. Additional BIEs are derived to solve for the additional unknowns introduced through enrichment while a numerical integration scheme is outlined for the evaluation of strongly singular and hypersingular enriched integrals. Results show a marked improvement in accuracy when enrichment is applied while tests also showed an improvement with the use of unenriched curved crack elements over flat elements. In addition the method compares favourably with FEM results obtained with higher DOF. Since there are far fewer elements in a BEM model than a corresponding FEM model, direct comparison of the formulation against FEM (both the unenriched FEM and the enriched XFEM) is difficult. It remains to undertake an exhaustive study based on accuracy obtained for equivalent CPU time, which is probably the most meaningful comparison. However, preliminary results show that the enriched DBEM described herein is very competitive. The accuracy and efficiency of the scheme show promise that similar benefits are likely to be found in 3D, and justify further research in implementing this type of enrichment using the approach (of Sukumar et al. [39]) that has been used in XFEM. Appendix A The procedure outlined in [32] is used to allow the evaluation of the hypersingular integrals arising in the present method where, for illustration purposes, the kernel Skij is used. For 2D elastostatics, Skij is given by

Skij ¼



@r  ð1  2mÞdij r ;k þ mðr ;j dik þ r ;i djk Þ  4r ;i r;j r;k þ 2mðni r ;j r ;k þ nj r ;i r ;k Þ @n  þ ð1  2mÞð2nk r ;i r ;j þ nj dik þ ni djk Þ  ð1  4mÞnk dij

l

2pð1  mÞr2

2

ð25Þ

All hypersingular integrals involving this term are multiplied by the shape function Na(n) and the Jacobian of transformation Jn(n). We are then left with an integral of the form

Z

þ1

1

Na ðnÞSkij wUl ðnÞJ n ðnÞdn

ð26Þ

which is of O(1/r2) when the source and field point coincide. The method is based on expressing the integrand seen in (26) in a Taylor series form where definitions are made to simplify later expressions.

636

R. Simpson, J. Trevelyan / Engineering Fracture Mechanics 78 (2011) 623–637

If the components of the field and source point locations are expressed as xi and yi respectively (in keeping with the notation of [32]), then the following Taylor series expansion about the point np can be written

  2 dxi  d xi  xi  yi ¼ ðn  np Þ þ 2  dn n¼np dn 

n¼np

ðn  np Þ2 þ  2

¼ Ai ðn  np Þ þ Bi ðn  np Þ2 þ    ¼ Ai d þ Bi d2 þ Oðd3 Þ;

ð27Þ

which defines the constants Ai and Bi along with the term d: = n  np. The constants A and C are also defined as 2 X

A :¼

!1=2 A2k

ð28Þ

Ak Bk

ð29Þ

k¼1

C :¼

2 X k¼1

However, to determine Ai and Bi (and therefore A and C), the first and second derivatives about the source point must be found. This is achieved by utilising the relevant shape functions and the nodal coordinates in the following way

dxi dNa a ¼ x dn dn i 2

d xi dn

2

ð30aÞ

2

¼

d Na dn2

xai

ð30bÞ

Now the derivative r,i can then be expressed as

r;i ¼

  xi  y i A i Bi Ak Bk ¼ þ  Ai 3 d þ Oðd2 Þ ¼: di0 þ di1 d þ Oðd2 Þ r A A A

ð31Þ

while the term 1/r2 can also be rewritten as

1 1 2C ¼  þ Oð1Þ r 2 A2 d2 A4 d S2 S1 þ Oð1Þ ¼: 2 þ d d

ð32Þ ð33Þ

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi It is also useful to express the Jacobian of transformation in terms of its components Ji(n) where J n ðnÞ ¼ J 1 ðnÞ2 þ J 2 ðnÞ2 and

J 1 ¼ A2 þ 2B2 d þ Oðd2 Þ

ð34aÞ

J 2 ¼ A1  2B1 d þ Oðd2 Þ

ð34bÞ

As a generalisation, these are written as

J k ¼ J k0 þ J k1 d þ Oðd2 Þ

ð35Þ

Finally, we express the shape functions Na and the enrichment functions

wUl

as Taylor expansions

 dNa  Na ðnÞ ¼ Na ðnp Þ þ ðn  np Þ þ    ¼ Na0 þ Na1 d þ Oðd2 Þ dn n¼np

ð36Þ

and

wUl ðnÞ ¼ wUl ðnp Þ þ

 dwUl   dn 

ðn  np Þ þ    ¼ wUl0 þ wUl1 d þ Oðd2 Þ:

ð37Þ

n¼np

The integrand in (26) can now be expressed as a Taylor series by substituting in expressions (31), (32), (35), (36) and (37) while also noting that Ji = niJn. By collecting all the terms that contain 1/d2 and 1/d where, due to the use of quadratic shape functions, any higher order terms are zero, the following expression can be written for the integrand

Na ðnÞSkij wUl ðnÞJ n ðnÞ

( ) 

S2 ðnp ÞNa0 hðnp ÞwUl0  U U U ¼D þ S2 ðnp Þ Na0 hðnp Þwl1 þ wl0 Na1 hðnp Þ þ gðnp ÞNa0 þ S1 N a0 hðnp Þwl0 =d d2 ð38Þ

where the constant D is defined as l/2p(1  m), and the terms h(np) and g(np) are given by

R. Simpson, J. Trevelyan / Engineering Fracture Mechanics 78 (2011) 623–637

637

hðnp Þ ¼ 2mðJ i0 dj0 dk0 þ J j0 di0 dk0 Þ þ ð1  2mÞð2J k0 di0 dj0 þ J j0 dik þ J i0 djk Þ  ð1  4mÞJ k0 dij

ð39Þ

 gðnp Þ ¼ 2ðdl1 J l0 þ dl0 J l1 Þ ð1  2mÞdk0 dij þ mðdj0 dik þ di0 djk Þ  4di0 dj0 dk0  þ 2m J i0 ðdj1 dk0 þ dj0 dk1 Þ þ J i1 dj0 dk0 þ J j0 ðdi1 dk0 þ di0 dk1 Þ þ J j1 di0 dk0  þ ð1  2mÞ 2ðJ k1 di0 dj0 þ J k0 ðdi1 dj0 þ di0 dj1 ÞÞ þ J j1 dik þ Ji1 djk  ð1  4mÞJ k1 dij

ð40Þ

where the summation rule is applied to the first two terms in expression (40). References [1] Tong P, Pian THH. On the convergence of the finite element method for problems with singularity. Int J Numer Methods Engng 1973;9:313–21. [2] Miranda ACO, Martha LF, Wawrzynek PA, Ingraffea AR. Surface mesh regeneration considering curvatures. Engng Comput 2009;25:207–19. [3] Bouchard PO, Bay F, Chastel Y. Numerical modelling of crack propagation: automatic remeshing and comparison of different criteria. Comp Methods Appl Mech Engng 2003;192:3887–908. [4] Benzley SE. Representation of singularities with isoparametric finite elements. Int J Numer Methods Engng 1974;8:537–45. [5] Henshell RD, Shaw KG. Crack-tip finite elements are unnecessary. Int J Numer Methods Engng 1975;9:495–507. [6] Barsoum RS. On the use of isoparametric finite elements in linear fracture mechanics. Int J Numer Methods Engng 1976;10:25–37. [7] Harrop LP. The optimum size of quarter-point crack tip elements. Int J Numer Methods Engng 1982;18:1101–3. [8] Ingraffea AR, Manu C. Stress-intensity factor computation in three dimensions with quarter-point elements. Int J Numer Methods Engng 1980;15:1427–45. [9] Babuška I, Melenk JM. The partition of unity method. Int J Numer Methods Engng 1996;40:727–58. [10] 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. [11] Duflot M. A study of the representation of cracks with level sets. Int J Numer Methods Engng 2007;70:1261–302. [12] Bordas S, Duflot M. Derivative recovery and a posteriori error estimate for extended finite elements. Comp Meth Appl Mech Engng 2007;196:3381–99. [13] Bordas S, Duflot M, Le P. A simple error estimator for extended finite elements. Commun Numer Methods Engng 2008;24:961–71. [14] Duflot M, Bordas S. A posteriori error estimation for extended finite elements by an extended global recovery. Int J Numer Methods Engng 2008;76:1123–38. [15] Karihaloo BL, Xiao QZ. Accurate determination of the coefficients of elastic crack tip asymptotic field by a hybrid crack element with p-adaptivity. Engng Fract Mech 2001;68:1609–30. [16] Belytschko T, Lu YY, Gu L. Element-free Galerkin methods. Int J Numer Methods Engng 1994;37:229–56. [17] Nguyen VP, Rabczuk T, Bordas S, Duflot M. Meshless methods: a review and computer implementation aspects. Math. Comput. Simul. 2008;79:763–813. [18] Rabczuk T, Belytschko T. A three-dimensional large deformation meshfree method for arbitrary evolving cracks. Comp Methods Appl Mech Engng 2007;196:2777–99. [19] Bordas S, Rabczuk T, Zi G. Three-dimensional crack initiation, propagation, branching and junction in non-linear materials by an extended meshfree method without asymptotic enrichment. Engng Fract Mech 2008;75:943–60. [20] Liew KM, Cheng Y, Kitipornchai S. Analyzing the 2D fracture problems via the enriched boundary element-free method. Solids Struct 2007;44:4220–33. [21] Yamada Y, Ezawa Y, Nishiguchi I. Reconsiderations on singularity or crack tip elements. Int J Numer Methods Engng 1979;14:1525–44. [22] Tanaka M, Itoh H. New crack elements for boundary element analysis of elastostatics considering arbitrary stress singularities. Appl Math Modell 1987;11:357–63. [23] Aliabadi MH. The boundary element method. Applications in solids and structures, vol. 2. Wiley Publishing; 2002. [24] Watson JO. Singular boundary elements for the analysis of cracks in plane-strain. Int J Numer Methods Engng 1997;38:2389–411. [25] Portela A, Aliabadi MH, Rooke DP. The dual boundary element method – effective implementation for crack problems. Int J Numer Methods Engng 1992;6:1269–87. [26] Mellings S, Baynham J, Adey RA. Fully automatic 3D crack growth with BEM. Boundary elements XXII. WIT Press; 2000. ISBN: 1-85312-824-4 [27] Rice JR. A path independent integral and approximate analysis of strain concentraion by notches and cracks. J Appl Mech 1968;35:379–86. [28] Legrain G, Moës N, Verron E. Robust and direct evaluation of J2 in linear elastic fracture mechanics with the X-FEM. Int J Numer Methods Engng 2008;76:1471–88. [29] Williams ML. Stress singularities resulting from various boundary conditions in angular corners of plates in extension. J Appl Mech 1952;19:526–8. [30] Eischen JW. An improved method for computing the J2 integral. Engng Fract Mech 1987;26:691–700. [31] Chang JH, Wu DJ. Computation of mixed-mode stress intensity factors for curved cracks in anisotropic elastic solids. Engng Fract Mech 2007;74:1360–72. [32] Guiggiani M, Krishnasamy G, Rudolphi TJ, Rizzo FJ. A general algorithm for the numerical-solution of hypersingular boundary integral equations. J Appl Mech – Trans ASME 1992;59:604–14. [33] Ishikawa H, Kitagawa H, Okamura H. J integral of a mixed mode crack and its application. Proc 3rd Int Conf Mech Behav Mater 1980;3:447–55. [34] Sih GC, Paris PC, Erdogan F. Crack-tip, stress-intensity factors for plane extension and plate bending problems. J Appl Mech 1962;29:306–12. [35] Liu YJ, Xu N. Modeling of interface cracks in fiber-reinforced composites with the presence of interphases using the boundary element method. Mech Mater 2000;32:769–83. [36] Laborde P, Pommier J, Renard Y, Salaün M. High-order extended finite element method for cracked domains. Int J Numer Methods Engng 2005;64:354–81. [37] Béchet E, Minnebo H, Moës N, Burgardt B. Improved implementation and robustness study of the X-FEM for stress analysis around cracks. Int J Numer Methods Engng 2005;64:1033–56. [38] Simpson R, Trevelyan J. A partition of unity enriched dual boundary element method for accurate computations in fracture mechanics. Comp Meth Appl Mech Engng. in press. [39] Sukumar N, Moës N, Moran B, Belytschko T. Extended finite element method for three-dimensional crack modelling. Int J Numer Methods Engng 2000;48:1549–70.