Computation of dynamic stress intensity factors in cracked functionally graded materials using scaled boundary polygons

Engineering Fracture Mechanics 131 (2014) 210–231

Contents lists available at ScienceDirect

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

Computation of dynamic stress intensity factors in cracked functionally graded materials using scaled boundary polygons Irene Chiong a,⇑, Ean Tat Ooi b, Chongmin Song a, Francis Tin-Loi a a b

School of Civil and Environmental Engineering, The University of New South Wales, Sydney, NSW 2031, Australia School of Science, Information Technology and Engineering, Federation University, Ballarat, VIC 3353, Australia

a r t i c l e

i n f o

Article history: Received 17 June 2014 Accepted 26 July 2014 Available online 23 August 2014 Keywords: Scaled boundary finite element Polygon element Functionally graded materials Fracture Dynamic stress intensity factors

a b s t r a c t In this paper, the recently developed scaled boundary polygons formulation for the evaluation of stress intensity factors in functionally graded materials is extended to elasto-dynamics. In this approach, the domain is discretized using polygons with arbitrary number of sides. Within each polygon, the scaled boundary polygon shape functions are used to interpolate the displacement field. For uncracked polygons, these shape functions are linearly complete. In a cracked polygon, the shape functions analytically model the stress singularity at the crack tip. Therefore, accurate dynamic stress intensity factors can be computed directly from their definitions. Only a single polygon is necessary to accurately compute the stress intensity factors. To model the material heterogeneity in functionally graded materials, the material gradients are approximated locally in each polygon using polynomial functions. This leads to semi-analytical expressions for both the stiffness and the mass matrices, which can be integrated straightforwardly. The versatility of the developed formulation is demonstrated by modeling five numerical examples involving cracked functionally graded specimens subjected to dynamic loads. Ó 2014 Elsevier Ltd. All rights reserved.

1. Introduction Functionally graded materials (FGMs) are a new class of customizable material arising from advancements in materials science and technology. Heat, oxidization, and fracture resistances can be increased by tailoring the volume fractions of the individual material constituents in a predetermined profile. This allows engineers to optimize the performance of FGM components under critical and extreme conditions for safety–critical structures in key modern industries such as microelectronics, aerospace and nuclear energy. As the application of FGMs become more common in engineering practice, the extent to which these materials can be tailored against damage becomes more important. Therefore, the capability to analyze the fracture behavior of structures and components made of FGMs is crucial to their design and development. For this purpose, both static and dynamic analyses can be used. While static analyses provide engineers and designers with an indication of the critical state of stress in a cracked body, real world structures are invariably loaded dynamically. As such, there is an emphasis for the accurate analyses of fractured FGM structures and components that are subjected to dynamic loads.

⇑ Corresponding author. E-mail address: [email protected] (I. Chiong). http://dx.doi.org/10.1016/j.engfracmech.2014.07.030 0013-7944/Ó 2014 Elsevier Ltd. All rights reserved.

I. Chiong et al. / Engineering Fracture Mechanics 131 (2014) 210–231

211

Nomenclature a B cd D

 g Ei E F I JðgÞ KðhÞ K L M N

X U W W ðsÞ WrL Wr ðsÞ Wr Wu

q Sn

r Dt t €b u u_ b ub

m n Z

crack length strain–displacement matrix dilatational wave speed constitutive material matrix strain circumferential coordinate coefficient matrices Young’s modulus force vector identity matrix Jacobian matrix vector of generalized stress intensity factors stiffness matrix characteristic length mass matrix circumferential shape function matrix domain polygon shape function matrix Schur transformation matrix strain modes singular stress modes at characteristic length stress modes singular stress modes displacement modes mass density block diagonal Schur matrix stress time step time vector of nodal acceleration vector of nodal velocity vector of nodal displacements Poisson’s ratio radial coordinate Hamiltonian coefficient matrix

The stress intensity factors (SIFs) play an important role in characterizing the fracture of FGMs. To compute the SIFs, both analytical and numerical methods can be used. The analytical solutions of the dynamic SIFs in FGMs have been reported by many researchers e.g. [1–10]. These studies are, however, limited to a finite crack in an infinite medium subjected to simple load cases. Although analytical methods provide direct closed-form solutions, the trade-off is a loss of generality and the need for complex solution techniques. Alternatively, numerical methods are more versatile and can be used to model a wider class of problems. The finite element method (FEM) is the most widely used numerical method for structural analysis. When standard FEM procedures are applied to fracture analyses, a very fine mesh is required at the crack tip in order to obtain reliable results. Otherwise, special elements such as the quarter-point element [11], or elements with embedded asymptotic expansions [12], have to be used. Moreover, the FEM requires additional post-processing methods such as the J-integral [13] or the M-integral [14]. These techniques need to be specifically formulated for nonhomogeneous materials such as FGMs. Wu et al. [15] reported the formulation of the J-integral for dynamic fracture analysis of FGMs. Song and Paulino [16] derived the M-integral for the dynamic response of FGMs, and found the method to be superior to the J-integral and the direct correlation technique in calculating the dynamic SIFs. The extended finite element method (XFEM) resolves some of the issues of mesh refinement about the crack tip in the FEM by enriching the shape functions of the finite elements in the vicinity of the crack with asymptotic crack tip functions [17]. Application of the XFEM to compute the dynamic SIFs in FGMs has been reported by Motamedi and Mohammadi [18], Singh et al. [19], Bayesteh and Mohammadi [20] and Liu et al. [21]. To compute the stress intensity factors, additional postprocessing techniques e.g. the J- and M-integrals are still required.

212

I. Chiong et al. / Engineering Fracture Mechanics 131 (2014) 210–231

The meshless methods [22] are able to analyze the response of structures without the need for a mesh. Like the XFEM, the meshless methods also require singular enrichment functions so that the stress field in the vicinity of the crack tip can be captured more accurately. Application of the meshless methods to calculate the dynamic SIFs in cracked FGMs have been reported by Zhang et al. [23], and Sladek et al. [24]. In the boundary element method [25–27], the computational burden is significantly reduced because only the boundary of the domain has to be discretized. However, this advantage is offset by the need for a fundamental solution, and also the complex mathematical procedures involved; e.g. integration of temporal shape functions [28], and convoluted fundamental solutions [29]. Application of the boundary element method to calculate the dynamic SIFs in FGMs was reported by Zhang et al. [30], and Sladek et al. [31]. While good results can be obtained using these methods, special singular elements are still required to compute accurate dynamic SIFs [26]. The scaled boundary finite-element method (SBFEM) [32] is a semi-analytical numerical technique that has emerged as an attractive alternative to model many types of engineering problems. In fracture mechanics, a study by Deeks and Wolf [33,34] showed that compared with the FEM, the SBFEM is able to reduce the computational effort considerably in problems involving stress singularities. This unique feature arises from the semi-analytical nature of the SBFEM solution [32,35,36]. In the radial direction emanating from the crack tip, the stress singularities are analytically represented [37,38]. Local mesh refinement at the crack tip is therefore avoided. Local enrichment with analytical asymptotic expansion is also not required. Furthermore, the SIFs can be computed directly from their definitions, as has been demonstrated in various studies on fracture in isotropic-, anisotropic- and bi-materials [33,35,37,39–45]. In problems arising in structural dynamics, a solution in the frequency domain for bounded media was originally derived by Song and Wolf [35,46]. An alternative formulation was developed by Yang et al. [47] using the Frobenious solution technique, which was successfully applied to the analysis of fracture in isotropic and orthotropic materials. In the time-domain, a SBFEM based super-element formulation was developed by Song [40] and was applied to compute the dynamic SIFs in homogeneous and simple bi-materials [40]. However, this formulation could only capture the inertial effects at low frequencies unless the domain was substructured into many subdomains. To improve capability of the SBFEM in capturing the inertial effects at high frequencies, Song [48] developed a continued fraction based SBFEM formulation. This approach has been shown to be accurate in capturing the dynamic SIFs and T-stresses in homogenous and bi-material specimens [43], and multi-material cracks and notches [49,50]. Recently, Chiong et al. [51] developed a novel SBFEM formulation that is capable of accurate and efficient computation of static stress intensity factors in FGMs. The formulation was implemented on polygons with arbitrary numbers of sides, and is based on the scaled boundary polygon formulation developed by Ooi et al. [52]. In order to model the response of FGMs, scaled boundary shape functions were introduced. These shape functions are obtained analytically within a polygon and have been shown to naturally include the strain singularities at a crack or a notch [51]. Higher order shape functions can be more easily constructed than other polygon elements reported in the literature e.g. [53]. The scaled boundary shape functions enable a polygon element to be formulated using standard finite element procedures. The material variation is locally approximated in each polygon using a polynomial function. The stiffness matrix is evaluated semi-analytically. In this paper, the SBFEM formulation developed by Chiong et al. [51] will be extended to model transient dynamic fracture problems in FGMs. This paper is organized as follows. In Section 2, geometrical requirements of the scaled boundary polygons are over viewed, and the scaled boundary polygon shape functions are summarized. Section 3 describes the scaled boundary polygon formulation for functionally graded materials, including the derivations for the stiffness and the mass matrices. In Section 4, procedures used to evaluate the dynamic SIFs in cracked FGMs are shown. In Section 5, the developed formulation is validated using five numerical examples. The major conclusions of this study are summarized in Section 6. 2. Overview of scaled boundary polygons 2.1. Scaled boundary transformation of geometry The SBFEM can be formulated on polygons with arbitrary number of sides. The geometry of the polygon has only to satisfy the SBFEM scaling requirement Song and Wolf [32]. This condition can be easily met for any convex polygon and some concave polygons. For more complex shapes, this requirement can always be achieved by sub-division. A polygon mesh can be generated from Delaunay triangulation, following the approach reported by Ooi et al. [52]. The polygons generated by this approach automatically satisfy the scaling requirement. Fig. 1a shows a polygon modeled by the SBFEM. The scaling center is chosen at a point where the scaling requirement is satisfied, and is usually the geometric center of the polygon. At the scaling center, the SBFEM coordinate system (n; g) is defined. The radial coordinate, n, varies from 0 at the scaling center to 1 at the polygon boundary. On the boundary, each edge is discretized by one dimensional finite elements. These elements can be of any order. Each edge can also be discretized using more than one element. The local coordinate, g, is defined for each line element on each polygon edge and has an interval of 1 6 g 6 þ1. The Cartesian coordinates of a point on a line element with M nodes are

xb ðgÞ ¼ NðgÞxb

ð1Þ

yb ðgÞ ¼ NðgÞyb

ð2Þ

I. Chiong et al. / Engineering Fracture Mechanics 131 (2014) 210–231

213

Fig. 1. Scaled boundary polygon representation: (a) Standard polygon; (b) cracked (a  0) or notched (a > 0) polygon.

where xb ¼ ½x1 x2 ; . . . ; xM  and yb ¼ ½y1 y2 ; . . . ; yM  are their respective nodal coordinates. NðgÞ ¼ ½ N 1 ðgÞ N 2 ðgÞ . . . N M ðgÞ  is the one-dimensional shape function vector of the line elements on the polygon boundary. The geometry of the domain is described by scaling the boundary along n. The Cartesian coordinates of a point inside the domain are given by the scaled boundary transformation equations

xðn; gÞ ¼ nNðgÞxb

ð3Þ

yðn; gÞ ¼ nNðgÞyb

ð4Þ

When a crack or a notch is modeled by a polygon, as shown in Fig. 1b, the scaling center is chosen as the crack or notch tip. The crack surfaces are defined by scaling the crack mouth nodes to the crack tip, and are not discretized. The coordinate transformation between the Cartesian coordinates ðx; yÞ and the scaled boundary coordinates ðn; gÞ can be performed in a manner similar to standard isoparametric finite elements. The Jacobian matrix on the boundary JðgÞ required for this transformation is Song and Wolf [32]

" JðgÞ ¼

xb ðgÞ xb ðgÞ;g

yb ðgÞ yb ðgÞ;g

# ð5Þ

where jJðgÞj is the determinant of jJðgÞj. For two-dimensional problems, an infinitesimal area in the domain, dX, is given in [54] as

dX ¼ jJðgÞjndndg

ð6Þ

On the boundary, an infinitesimal length dC is mapped as

dC ¼ DC ðgÞdg

ð7Þ

where

DC ¼

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi  2  2ffi NðgÞ;g xb þ NðgÞ;g yb

ð8Þ

2.2. Scaled boundary polygon shape functions For a sector covered by a line element on the boundary, the displacement field uðn; gÞ is interpolated using the scaled boundary shape functions [51] as

uðn; gÞ ¼ Uðn; gÞub ¼ Nu ðgÞWu nSn W1 u ub

ð9Þ

where Nu ðgÞ is

Nu ðgÞ ¼



N1 ðgÞ

0

N2 ðgÞ

0

0

N1 ðgÞ

0

N2 ðgÞ

...

0

NM ðgÞ

0

...

0

NM ðgÞ

 ð10Þ

214

I. Chiong et al. / Engineering Fracture Mechanics 131 (2014) 210–231

and Sn and Wu are the block diagonal Schur matrix and the corresponding transformation matrix related to the modal displacements, respectively. These are obtained from a block-diagonal Schur decomposition of the Hamiltonian matrix Z Song [36]

" Z¼

#

T E1 0 E1

E1 0

T E2 þ E1 E1 0 E1

E1 E1 0

ð11Þ

where E0 ; E1 and E2 are coefficient matrices that are related to the geometry and material properties in each polygon

E0 ¼ E1 ¼ E2 ¼

Z

þ1

1 Z þ1 1 Z þ1 1

BT1 ðgÞDB1 ðgÞjJðgÞjdg

ð12Þ

BT2 ðgÞDB1 ðgÞjJðgÞjdg

ð13Þ

BT2 ðgÞDB2 ðgÞjJðgÞjdg

ð14Þ

and D is the material constitutive matrix. They can be assembled element by element on the polygon boundary. B1 ðgÞ and B2 ðgÞ are the strain–displacement matrices Song [36]

B1 ðgÞ ¼ b1 ðgÞNu ðgÞ

ð15Þ

B2 ðgÞ ¼ b2 ðgÞNu ðgÞ;g

ð16Þ

with

2 1 6 b1 ðgÞ ¼ 4 jJðgÞj

yb ðgÞ;g 0 xb ðgÞ;g

0

3

xb ðgÞ;g 7 5

yb ðgÞ;g 2 3 0 yb ðgÞ 1 6 7 xb ðgÞ 5 b2 ðgÞ ¼ 4 0 jJðgÞj xb ðgÞ yb ðgÞ

ð17Þ

ð18Þ

For uncracked polygons, the scaled boundary shape functions can be shown to be linearly complete i.e. they can reproduce rigid body motions and constant strain modes [51]. For polygons modeling cracks and notches, the derivatives of these shape functions naturally include the strain singularity at the crack/notch tip. This allows problems with cracks and notches to be modeled with higher accuracy without any enrichment functions or fine crack/notch tip meshes. The strain field can be derived from the displacement field in Eq. (9) as

ðn; gÞ ¼ Bðn; gÞub

ð19Þ

where Bðn; gÞ is the scaled boundary strain displacement matrix

Bðn; gÞ ¼ W ðgÞnSn I W1 u

ð20Þ

and the strain mode W ðgÞ is

W ðgÞ ¼ B1 ðgÞWu Sn þ B2 ðgÞWu

ð21Þ

3. Scaled boundary polygon formulation for elastodynamics of functionally graded materials 3.1. Modeling of material heterogeneity in functionally graded materials In FGMs, the material properties i.e. Young’s modulus E, Poisson’s ratio m and mass density q vary according to some predetermined profile. In this paper, the heterogeneity of E; m and q is accounted for by assuming that their variations within each polygon are polynomial functions of the forms

Dðx; yÞ ¼ D0 þ D1 x þ D2 y þ D3 x2 þ D4 xy þ D5 y2 þ . . .

ð22Þ

qðx; yÞ ¼ q0 þ q1 x þ q2 y þ q3 x2 þ q4 xy þ q5 y2 þ . . .

ð23Þ

and

where D is the standard elastic constitutive matrix.The coefficient matrices Di in Eq. (22) and constants the qi in Eq. (23) are determined from a least squares fit over each polygon. The Gaussian/Gauss–Lobatto integration points are used for fitting. To incorporate Eqs. (22) and (23) into the scaled boundary finite element formulation, they are first written in terms of scaled boundary coordinates ðn; gÞ by substituting Eqs. (3) and (4) for x and y, resulting in

215

I. Chiong et al. / Engineering Fracture Mechanics 131 (2014) 210–231

X ðkÞ D ðgÞnk

Dðn; gÞ ¼ Dð0Þ ðgÞn0 þ Dð1Þ ðgÞn1 þ Dð2Þ ðgÞn2 þ . . . ¼

ð24Þ

k¼0

and

qðn; gÞ ¼ qð0Þ n0 þ qð1Þ ðgÞn1 þ qð2Þ ðgÞn2 þ . . . ¼

X

qðkÞ ðgÞnk

ð25Þ

k¼0

where k is the order of the complete polynomial used to approximate the material variation in each polygon. The assumptions made in Eqs. (24) and (25) imply that the size of the polygons should be sufficiently small and k should be sufficiently large to accurately model the material gradient. The choice of k depends on the number of fitting points available in a polygon, and also the number of polygons in mesh. It was shown in Chiong et al. [51] that a second order polynomial k ¼ 2 is adequate for most problems when there is a sufficient number of polygons in the mesh. 3.2. Governing equations for elastodynamics In linear elastodynamics, the equilibrium condition of a polygon can be formulated using the virtual work principle

Z

dT rdX þ

Z

X

€ dX ¼ duT qu

Z

duT tdC þ

C

X

Z

duT bdX

ð26Þ

X

€ is the acceleration field, r is the stress field, and t where d is the virtual strain field, du is the virtual displacement field, u and b are surface tractions and body loads respectively. For two-dimensional problems, dC is an infinitesimal line on the polygon boundary and dX is an infinitesimal area in the polygon. Using Hooke’s law, r ¼ D, and substituting Eqs. (9) and (19) for u and , in Eq. (26) results in

duTb

Z

BT ðn; gÞDBðn; gÞdXub þ

X

Z

€b UT ðn; gÞqUðn; gÞu

X



¼ duTb

Z

UT ðn; gÞtdC þ

Z

C

 UT ðn; gÞbdX

ð27Þ

X

Invoking the arbitrariness of the virtual displacement dub , Eq. (27) becomes

Z

 Z  Z Z €b ¼ BT ðn; gÞDBðn; gÞdX ub þ UT ðn; gÞqUðn; gÞdX u UT ðn; gÞtdC þ UT ðn; gÞbdX

X

X

C

ð28Þ

X

3.2.1. Stiffness matrix for functionally graded materials The first term on the left hand side of Eq. (28) is the stiffness matrix of the polygon

Kp ¼

Z

BT ðn; gÞDBðn; gÞdX

ð29Þ

X

Substituting Eq. (20) for Bðn; gÞ, and Eq. (6) for dX, the stiffness matrix becomes

Kp ¼ WT u

Z

1

Z

1 0

1

T

nSn I WT ðgÞDW ðgÞnSn I njJðgÞjdndgW1 u

ð30Þ

Eq. (30) can be first integrated numerically in the g direction and then analytically in the n direction. Substituting the expression for the elastic constitutive matrix from Eq. (24) into Eq. (30), leads to the expression

Kp ¼ WT u

XZ 0

k¼0

1

 T nSn I YðkÞ nSn þkI dn W1 u

ð31Þ

where

YðkÞ ¼

Z

1 1

WT ðgÞDðkÞ ðgÞW ðgÞjJðgÞjdg

ð32Þ

YðkÞ is integrated numerically for each term k using Gauss/ Gauss–Lobatto techniques. For a single term k, the n integral in Eq. (31) is defined as

XðkÞ ¼

Z

1

T

nSn I YðkÞ nSn þkI dn

ð33Þ

0

Using integration by parts, it can be shown that XðkÞ is the solution of the Lyapunov equation T

ðSn þ 0:5kIÞ XðkÞ þ XðkÞ ðSn þ 0:5kIÞ ¼ YðkÞ

ð34Þ

As the coefficient matrix Sn þ 0:5kI is in Schur form, only a back substitution is required. The stiffness matrix, Kp can be computed once XðkÞ has been determined for all the k terms. Eq. (31) therefore simplifies to

216

I. Chiong et al. / Engineering Fracture Mechanics 131 (2014) 210–231

! X ðkÞ W1 X u

Kp ¼ WT u

ð35Þ

k¼0

3.2.2. Mass matrix for functionally graded materials The second term in Eq. (28) is the mass matrix of the polygon

Mp ¼

Z

UT ðn; gÞqUðn; gÞdX

ð36Þ

X

Substituting the scaled boundary shape functions in Eq. (9), and dX from Eq. (6), into Eq. (36), results in the expression

Mp ¼ WT u

Z

1

0

Z

þ1

T

nSn WTu NTu ðgÞqNu ðgÞWu nSn njJðgÞjdndgW1 u

1

ð37Þ

Substituting the approximation of q ¼ qðn; gÞ in Eq. (25) leads to the expression

Mp ¼ WT u

n Z X k¼0

 T nSn HðkÞ nSn ndn W1 u

ð38Þ

NTu ðgÞqðkÞ ðgÞNu ðgÞjJðgÞjdgWu

ð39Þ

1

0

where HðkÞ is

HðkÞ ¼ WTu

Z

þ1

1

and is integrated numerically for each term k using Gauss/Gauss–Lobatto quadrature. For each of these terms, the integral with respect to n in Eq. (38) is

Z

GðkÞ ¼

1



T nSn þI HðkÞ nSn þkI dn

ð40Þ

0

Using integration by parts, it can be shown that GðkÞ is the solution of the following Lapunov equation T

GðkÞ ðSn þ ð1 þ 0:5kÞIÞ þ ðSn þ ð1 þ 0:5kÞIÞ GðkÞ ¼ HðkÞ

ð41Þ

The mass matrix in Eq. (38) now simply reduces to

M¼

WT u

! X ðkÞ W1 G u

ð42Þ

k¼0

3.2.3. Polygon load vector The terms on the right-hand-side of Eq. (28) are the equivalent nodal force vector, Fp acting on the polygon

Fp ¼

Z

UT ðn; gÞtdC þ

Z

C

UT ðn; gÞbdX

ð43Þ

X

The first term corresponds to the distributed load on the polygon boundary, and the second corresponds to the body loads. The first term in Eq. (43) can be simplified by considering that at the polygon boundary, n ¼ 1. Therefore, Uðn; gÞ ¼ Nu ðgÞ. Using the relation in Eq. (7), the vector of distributed loads is thus

Z

UT ðn; gÞtdC ¼

Z

C

þ1

1

NTu ðgÞtDC ðgÞdg

ð44Þ

For a constant body load, the second term on the right-hand-side of Eq. (43) can be integrated numerically in g, and analytically in n. Substituting Eqs. (9) and (6), for Uðn; gÞ and X, the body load vector becomes

Z X

UT ðn; gÞbdX ¼

Z 0

1

Z

þ1

1



T Nu ðgÞWu nSn W1 jJðgÞjdgdnb u

ð45Þ

Integrating analytically in the n direction yields

Z X

1 T UT ðn; gÞbdX ¼ WT u ðSn þ 2IÞ Wu

Z

þ1

1

 NTu ðgÞjJðgÞjdg b

ð46Þ

3.2.4. Assembly of system of equations Substituting Eq. (35) for Kp , Eq. (42) for Mp , and Eq. (43) for Fp into Eq. (28) leads to a set of coupled second-order linear differential equations

I. Chiong et al. / Engineering Fracture Mechanics 131 (2014) 210–231

€ b þ Kp ub ¼ Fp Mp u

217

ð47Þ

Mp ; Kp , and Fp can be assembled polygon by polygon by similar means to the FEM. The assembled system of equations can be solved using standard time integration techniques. In this study, the Newmark method [55] with c ¼ 0:5 and b ¼ 0:25, is used. 4. Evaluation of dynamic stress intensity factors In the SBFEM, a crack or a notch is modeled by a single polygon, such as that shown in Fig. 2. Under this condition, some of the real parts of the eigenvalues in Sn , kðSn Þ, are between 0 and 1. These eigenvalues lead to singular strain fields in Eq. (20), as n ! 0. Denoting these diagonal blocks in Sn as SðsÞ , with corresponding displacement modes WðsÞ u , the singular strain modes of the strain field in Eq. (21) are ðsÞ ðsÞ ðsÞ WðsÞ  ðgÞ ¼ B1 ðgÞWu S þ B2 ðgÞWu

ð48Þ

ðsÞ

The singular stress modes Wr ðgÞ are related to the singular strain modes by the relationship ðsÞ WðsÞ r ðgÞ ¼ Dtip W ðgÞ

ð49Þ

were Dtip is the material constitutive matrix evaluated at the crack tip. At the crack tip, the singular stress field rðsÞ ðn; gÞ is therefore S rðsÞ ðn; gÞ ¼ WðsÞ r ðgÞn

ðsÞ

I

ðsÞ

ðW1 u Þ ub

ðsÞ ðW1 u Þ

ð50Þ ðsÞ

W1 u

where are the rows of corresponding to S . To compute the generalized SIFs, the singular stress field rðsÞ ðn; gÞ in Eq. (50) is first transformed into polar coordinates using the relations [38]

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ^rðn; gÞ ¼ n x2 ðgÞ þ y2 ðgÞ

ð51Þ

hðgÞ ¼ arctanðyðgÞ=xðgÞÞ

ð52Þ

Using Eqs. (51) and (52), and introducing a characteristic length L; n is expressed at an angle h as

n¼

^r L ^r ¼ rðhÞ rðhÞ L

ð53Þ

where rðhÞ is the distance from the scaling center to the boundary along the radial line at angle h. The singular stress field in Eq. (50) can be expressed in the polar coordinates following standard transformation procedures as

^ S rðsÞ ð^r; hÞ ¼ WðsÞ rL ðgðhÞÞðr =LÞ

ðsÞ

I

ðsÞ

ðW1 u Þ ub

ð54Þ

where the stress modes at the characteristic length L are ðsÞ

ðsÞ

S WrL ðhÞ ¼ WðsÞ r ðgðhÞÞðL=rðhÞÞ

I

The generalized stress intensity factors KðhÞ at angle h are defined in

Fig. 2. Scaled boundary cracked polygon (the scaling center is chosen at the crack tip).

ð55Þ

218

I. Chiong et al. / Engineering Fracture Mechanics 131 (2014) 210–231

1

S rðsÞ ð^r; hÞ ¼ pffiffiffiffiffiffiffiffiffi WðsÞ r =LÞ rL ðhÞð^ 2pL

ðsÞ

I



1 ðsÞ WrL ðhÞ KðhÞ

ð56Þ

Comparing Eq. (56) with Eq. (54), the generalized SIFs can be evaluated as

KðhÞ ¼

pffiffiffiffiffiffiffiffiffi ðsÞ ðsÞ 2pLWrL ðhÞðW1 u Þ ub

ð57Þ

5. Numerical examples To illustrate the accuracy of the present method in computing the dynamic SIFs in FGMs, five numerical examples are modeled. The computational domain is discretized by a mesh of polygons having arbitrary numbers of sides. Each polygon edge is discretized using one linear element. As was substantiated by Chiong et al. [51], a cracked polygon requires approximately 50 nodes to model the angular variation of the stress singularity accurately. As the polygons in the mesh used in this study has typically 5 or more edges, each edge on the cracked polygon is discretised using ten linear elements. For the FGM specimens, the material variations are locally approximated by the polynomial function, Eqs. (24) and (25). For all the examples presented, a complete second order polynomial, k ¼ 2, has been used. For time integration, the implicit Newmark scheme, which is unconditionally stable, is employed. Parametric studies show that accurate results are obtained when the size of time step Dt is selected to satisfy the Courant–Friedrichs–Lewy (CFL) condition which states that the distance travelled by the dilatational wave speed in one time step is at most the size of a line element in the mesh i.e.

c d Dt 6 h

ð58Þ

where h is the representative length of the edges of the polygons in the mesh, and cd is the dilatational wave speed. In plane strain the dilatational wave speed is given by

sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Eð1  mÞ cd ¼ qð1 þ mÞð1  2mÞ

ð59Þ

and in plane stress the dilatational wave speed is

sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi E cd ¼ qð1  m2 Þ

ð60Þ

5.1. Center cracked tension specimen A center cracked tension specimen with a horizontal central crack, as shown in Fig. 3a, is considered. The dimensions of the specimen are 2W ¼ 20 mm and 2H ¼ 40 mm. The crack length is 2a ¼ 4:8 mm. The Young’s modulus E and mass density q are graded according to the relations

E ¼ E0 expðbyÞ

Fig. 3. Center crack tension problem: (a) geometry; (b) coarse mesh; (c) medium mesh; (d) fine mesh.

ð61Þ

219

I. Chiong et al. / Engineering Fracture Mechanics 131 (2014) 210–231

and

q ¼ q0 expðbyÞ

ð62Þ

where E0 ¼ 199:992 GPa, and q0 ¼ 5000 kg=m3 . b is a parameter with units of mm1 . In this example three cases are considered, b ¼ 0 (homogeneous), 0:05, and 0:1. The Poisson’s ratio m ¼ 0:3 and is assumed to be constant throughout the specimen. Plane strain conditions are assumed. The specimen is discretized by polygon meshes of arbitrary number of sides. Fig. 3b shows the coarse mesh with 60 polygons and 202 nodes, Fig. 3c shows the medium mesh with 188 polygons and 482 nodes, and Fig. 3d shows the fine mesh with 642 polygons and 1442 nodes.

5.1.1. Modal analysis A modal analysis is first performed to validate the present formulation. For this analysis, the lower edge the specimen shown in Fig. 3a is clamped. The convergence of the first ten natural modal frequencies is investigated for the meshes shown in Fig. 3b–d. The result obtained from a very fine mesh of 328,571 PLANE183 elements with 988,291 nodes generated using ANSYS is used as a reference solution. Table 1a shows the first ten modal frequencies for the homogeneous specimen ðb ¼ 0Þ. The results obtained using the present method are in good agreement with the reference values calculated using ANSYS. The modal frequencies obtained using the coarse polygon mesh shown in Fig. 3b are within 2% of the reference solution. Results from the fine polygon mesh, shown in Fig. 3d, converge to within 0.1% of the reference solution. The natural frequencies for the FGM specimen are shown in Table 1b for the FGM with b ¼ 0:05, and in Table 1c for the FGM with b ¼ 0:1. The convergence trend for both FGM cases is similar to the homogeneous case. For both b ¼ 0:05 and b ¼ 0:1, the results obtained using the coarse and fine polygon meshes are within 2% and 0.1% of the reference solution, respectively.

Table 1 Convergence of the first ten natural modal frequencies of the center cracked tension problem. Mode

ANSYS

Scaled boundary polygons

(Converged Result)

Coarse

Medium

Fine

(a) Homogeneous b ¼ 0 1 2 3 4 5 6 7 8 9 10

11.347 40.757 42.716 87.185 112.30 115.96 144.34 147.24 155.81 157.11

11.404 40.836 42.946 88.072 113.31 116.98 146.15 149.76 158.82 159.67

11.367 40.790 42.787 87.440 112.60 116.27 144.89 147.97 156.74 157.86

11.354 40.769 42.739 87.265 112.38 116.06 144.50 147.45 156.07 157.33

(b) FGM b ¼ 0:05 1 2 3 4 5 6 7 8 9 10

6.2738 26.315 32.736 80.939 111.19 112.02 139.41 146.61 153.00 154.63

6.3802 26.614 33.061 81.560 111.87 111.91 141.72 149.92 156.21 157.31

6.2988 26.371 32.818 81.179 111.51 112.19 140.07 147.42 153.93 155.41

6.2795 26.324 32.758 81.018 111.29 112.11 139.56 146.84 153.21 154.86

(c) FGM b ¼ 0:1 1 2 3 4 5 6 7 8 9 10

3.1961 15.396 23.180 76.713 112.05 112.84 133.70 149.23 153.41 155.24

3.2615 15.733 23.580 76.531 110.63 113.63 136.27 152.02 156.81 157.30

3.2180 15.455 23.263 76.855 112.05 113.23 134.45 150.05 154.25 155.89

3.2006 15.404 23.198 76.783 112.12 112.94 133.85 149.47 153.55 155.45

220

I. Chiong et al. / Engineering Fracture Mechanics 131 (2014) 210–231

5.1.2. Transient analysis A transient dynamic fracture analysis of the center cracked tension problem is now performed. The specimen shown in Fig. 3a is loaded along the top and bottom edges with a time dependent uniform traction that varies according as PðtÞ ¼ PHðtÞ, where HðtÞ is the Heaviside step function. The homogeneous plate is considered first. The coarse and fine polygon meshes shown in Fig. 3b and d, are used for the analysis. For the homogeneous and three FGM cases studied, the dilatational wave speed throughout the plate, according to Eq. (59), is 7.338 mm/ls. The value h of the coarse and fine meshes are 2.5 mm and 0.6 mm respectively. To satisfy the condition in Eq. (58), the maximum allowable time step is Dt ¼ 0:08 ls. A time step of Dt ¼ 0:05 ls was selected for the analyses. pffiffiffiffiffiffi The time history of K I =P pa is shown in Fig. 4. The FEM solution reported in Song and Paulino [16] and the SBFEM solution reported in Song and Vrcelj [43] are used for comparison. The present results are smooth and do not display the small amplitude oscillations reported in Song and Paulino [16]. Overall, good agreement is observed between the present results, obtained using both the coarse and the fine mesh, and the reference solutions. When compared with the finite element results of Song and Paulino [16], which were obtained from a mesh of 816 Q8 (eight-noded quadrilateral) and 146 T6 (six-noded triangular) elements, the present method achieves similar accuracy using significantly coarser meshes. The FGM specimens are now considered. As the material varies in the y-direction, the dynamic SIFs of the left and right crack tips are equal. Therefore, only the normalized dynamic SIFs computed from the right crack tip are presented. The pffiffiffiffiffiffi dynamic SIFs are normalized by P pa. Two polygon meshes are employed for the analysis: The coarse mesh in Fig. 3b, and the fine mesh in Fig. 3d. As the ratio of q to E remains constant throughout the specimen, the dilatational wave speed and time step chosen is the same as it is in the homogeneous case. Fig. 5a and b shows the normalized values of K I and K II of the FGM plate for the case of b ¼ 0:05 over a period of 14 ls. The normalized dynamic SIFs for the FGM plate with b ¼ 0:1, are shown in Fig. 6a and b. The finite element solution obtained with a mesh of 816 Q8 and 146 T6 elements reported in Song and Paulino [16] is used for comparison. The present results are smooth, and do not display the small amplitude oscillations observed in the finite element results of Song and Paulino [16]. Overall the results are in good agreement. It is noted that the magnitude of K II is much smaller than K I . To further verify the results, an additional analysis was performed with a finer polygon mesh of 2426 polygons and 12,652 nodes. No appreciable difference with the results obtained from the fine polygon mesh is observed and results of this analysis are not shown.

5.1.3. Sensitivity of numerical results with respect to time step size The effect of the time step size on the accuracy of the result is now investigated. To evaluate the effect of the size of time step on the accuracy of the results the transient analyses is performed with four different sized time steps. The FGM case where b ¼ 0:1 is chosen for this analysis. The fine mesh shown in Fig. 3d is used. h for this mesh is approximately 1 mm and the dilatational velocity cd ¼ 7:338 mm=ls. The CFL condition Eq. (58) requires a time step of at least 0.14 ls for this mesh. Four time steps were trialled. Compared with the time step based on the CFL condition (0:14 ls), two time steps (Dt ¼ 0:4 ls and 0:2 ls) were larger and two (Dt ¼ 0:1 and 0.07 ls) were smaller. The results for K I at the right crack tip were pffiffiffiffiffiffi recorded. The values, normalized using by P pa, are presented in Fig. 7. As evidenced in the plotted results, the transient response obtained using the two larger time steps results in a loss of accuracy. Significant differences to the reference result are observed, and are especially pronounced near the peaks. The responses obtained with the two smaller time steps are very close to each other and in good agreement with the reference solution. This result indicates that the size of time step determined by Eq. (58) leads to accurate results without excessive number of time steps.

2.5

Normalized K

I

2 1.5 1 0.5

Present (Coarse) Present (Fine) Song and Paulino (2006) Song and Vrcelj (2008)

0 0

2

4

6

8

10

12

14

Time (μs) Fig. 4. Normalized K I for the homogeneous center cracked tension problem.

I. Chiong et al. / Engineering Fracture Mechanics 131 (2014) 210–231

Normalized KI

2.5 2 1.5 1 0.5 Present (Coarse) Present (Fine) Song and Paulino (2006)

0 0

2

4

6

8

10

12

14

Time (μs) 0.15

Normalized KII

0.1 0.05 0 −0.05 Present (Coarse) Present (Fine) Song and Paulino (2006)

−0.1 0

2

4

6

8

10

12

14

Time (μs) Fig. 5. Normalized dynamic SIFs for the FGM b ¼ 0:05 of the center cracked tension problem: (a) normalized K I ; (b) normalized K II .

2.5

Normalized K

I

2 1.5 1 0.5 0

Present (Coarse) Present (Fine) Song and Paulino (2006)

−0.5 0

2

4

6

8

10

12

14

10

12

14

Time (μs) 0.3

Normalized K

II

0.2 0.1 0 −0.1 −0.2

Present (Coarse) Present (Fine) Song and Paulino (2006)

−0.3 0

2

4

6

8

Time (μs) Fig. 6. Normalized dynamic SIFs for the FGM b ¼ 0:1 of the center cracked tension problem: (a) normalized K I ; (b) normalized K II .

221

222

I. Chiong et al. / Engineering Fracture Mechanics 131 (2014) 210–231

2.5

Normalized K I

2 1.5 1

Δt=0.4μs Δt=0.2μs Δt=0.1μs Δt=0.07μs Song and Paulino(2006)

0.5 0 −0.5 −1 0

2

4

6

8

10

12

14

Time (μs) Fig. 7. Normalized mode I dynamic SIFs at the right crack tip for 4 different time steps: 0.4 ls, 0.2 ls, 0.1 ls and 0.07 ls for b ¼ 0:1.

Fig. 8. Inclined center cracked tension problem: (a) geometry and boundary conditions; (b) coarse mesh; (c) fine mesh.

5.2. Inclined center crack tension specimen A specimen with an inclined crack at its center shownpin ffiffiffi Fig. 8a is considered. The dimensions of the problem are W ¼ 15 mm and H ¼ 30 mm, and the crack length is a ¼ 10= 2 mm. The crack is inclined at an angle h ¼ 45 from the horizontal axis. An external uniform traction, which varies according to the Heaviside step function PðtÞ ¼ PHðtÞ, is applied to the top and bottom edges of the specimen. The materials properties, Young’s modulus E and mass density q, vary exponentially in the x-direction resulting in a constant ratio E=q and are given as

E ¼ E0 expðbxÞ

ð63Þ

q ¼ q0 expðbxÞ

ð64Þ

and

where E0 ¼ 199:992 GPa and q0 ¼ 5000 kg=m3 . b is a material parameter with units of mm1 . Three cases: b ¼ 0:05; b ¼ 0:1, and b ¼ 0:15 are considered. The Poisson’s ratio m ¼ 0:3 is constant throughout the specimen. Plane strain conditions are assumed. For this example, two polygon meshes are considered. Fig. 8b shows the coarse mesh with 169 polygons and 514 nodes and Fig. 8c shows the fine mesh with 589 polygons and 1332 nodes. For all three cases, the dilatational wave speed throughout the specimen is cd ¼ 7:338 mm=ls. In the polygon meshes used, the value h is 2 mm and 1 mm for the coarse and fine mesh respectively. To satisfy the condition in Eq. (58), the maximum time step Dt ¼ 0:14 ls is required. For the present analysis, Dt ¼ 0:1 ls is selected.

Normalized Dynamic SIFs (Left)

I. Chiong et al. / Engineering Fracture Mechanics 131 (2014) 210–231

1 KI

0.5

Present (Coarse) Present (Fine) Song and Paulino (2006)

0 −0.5

K

II

−1 0

Normalized Dynamic SIFs (Right)

223

5

10 Time (μs)

15

20

1.5 1

K

I

Present (Coarse) Present (Fine) Song and Paulino (2006)

0.5 0 −0.5 K

II

−1 −1.5 0

5

10 Time (μs)

15

20

Normalized Dynamic SIFs (Left)

Fig. 9. Normalized dynamic SIFs K I and K II for the FGM b ¼ 0:05 case of the inclined center cracked tension problem: (a) left crack tip; (b) right crack tip.

1 K

I

0.5

Present (Coarse) Present (Fine) Song and Paulino (2006)

0 −0.5

K

II

−1 0

5

10

15

20

Normalized Dynamic SIFs (Right)

Time (μs) 2 1.5 K

I

1 0.5

Present (Coarse) Present (Fine) Song and Paulino (2006)

0 −0.5 KII

−1 −1.5 0

5

10

15

20

Time (μs) Fig. 10. Normalized dynamic SIFs K I and K II for the FGM b ¼ 0:1 case of the inclined center cracked tension problem: (a) left crack tip; (b) right crack tip.

224

I. Chiong et al. / Engineering Fracture Mechanics 131 (2014) 210–231

pffiffiffiffiffiffi Figs. 9–11 show the dynamic SIFs, normalized by P pa, for b ¼ 0:05, b ¼ 0:1 and b ¼ 0:15, respectively. The results obtained using the 208 Q8 and 198 T6 finite element mesh, reported in Song and Paulino [16], are used for comparison. Excellent agreement between the present numerical results and the reference solution is observed. 5.3. Circular hole

Normalized Dynamic SIFs (Left)

A rectangular FGM specimen with two cracks emanating from a central circular hole, shown in Fig. 12a, is investigated. The geometry of the plate is W ¼ 15 mm, and H ¼ 30 mm. The hole has a radius of r ¼ 3:75 mm. Two cracks emanate from

1 KI 0.5 Present (Coarse) Present (Fine) Song and Paulino (2006)

0 −0.5

K

II

−1 0

5

10

15

20

Normalized Dynamic SIFs (Right)

Time (μs)

2 K

I

1

Present (Coarse) Present (Fine) Song and Paulino (2006)

0 K

−1

II

0

5

10

15

20

Time (μs) Fig. 11. Normalized dynamic SIFs K I and K II for the FGM b ¼ 0:15 case of the inclined center cracked tension problem: (a) left crack tip; (b) right crack tip.

Fig. 12. Circular hole problem: (a) geometry and boundary conditions; (b) coarse mesh; (c) fine mesh.

I. Chiong et al. / Engineering Fracture Mechanics 131 (2014) 210–231

225

the hole, at an angle h ¼ 30 from the horizontal axis. The distance between the two crack tips is 2a ¼ 15 mm. An external uniform distributed load is applied instantaneously to the top and bottom edges according to the Heaviside step function PðtÞ ¼ PHðtÞ. The Young’s modulus E (MPa), and mass density q (kg=m3 ) vary in the x-direction according to

EðxÞ ¼ 7470:5 þ 3659:5

x

ð65Þ

W

and

x

qðxÞ ¼ 1380 þ 432

ð66Þ

W

The Poisson’s ratio is 0.3 and is constant throughout the specimen. Plane strain conditions are assumed. Two polygon meshes are used to discretize the specimen. Fig. 12b shows the coarse mesh with 176 polygons and 570 nodes, and Fig. 12c shows the fine mesh with 348 polygons and 814 nodes. The reference result in [16] was obtained using a mesh of 1350 Q8 and 204 Q6 elements. The dilatational wave speed varies throughout the specimen. The maximum value of cd ¼ 2:8755 mm=ls occurs along the far right edge of the specimen, and the minimum value of cd ¼ 2:3263 mm=ls occurs along the far left edge of the specimen. h is 3 mm and 1.5 mm for the coarse and fine mesh, respectively. In order to satisfy the condition in Eq. (58), the time step should not exceed 0:52 ls. A time step of Dt ¼ 0:5 ls is used for pthe ffiffiffiffiffiffiffiffiffiffianalysis. pffiffiffiffiffiffi In Fig. 13, the normalized dynamic SIFs (K I =P pa and K II =P ðpaÞ) computed using the present method are compared pffiffiffiffiffiffi with the reference values in [16]. The dynamic SIFs presented are normalized by P pa. Excellent agreement is observed between the present numerical results and the reference solutions. 5.4. Edge cracked specimen subjected to impact load

Normalized Dynamic SIFs

The edge cracked specimen shown in Fig. 14a is considered. The dimensions are W ¼ 200 mm, H ¼ 300 mm and a ¼ 50 mm. A constant velocity v ¼ 6:5 m/s is imposed on the upper half of the left most boundary above the crack tip.

2 KI Right tip KI Left tip

1

0 KII Left tip

Present (Coarse) Present (Fine) Song and Paulino (2006)

−1

0

5

10

15

20

25

K Right tip II

30

35

40

Time (μs) Fig. 13. Normalized K I and K II for the left and right crack tips of the FGM circular hole problem.

Fig. 14. Edge cracked plate problem: (a) geometry and boundary conditions; (b) coarse mesh; (c) fine mesh.

I. Chiong et al. / Engineering Fracture Mechanics 131 (2014) 210–231

0.5

Dynamic SIFs (MPa m )

226

Present (Coarse) Present (Fine) Song and Paulino (2006)

30 20

KII

10 0 KI

−10 0

5

10

15

20

25

30

Time (μs) Fig. 15. Normalized K I and K II for the FGM the edge cracked specimen subjected to impact load.

The Young’s modulus of the specimen is characterized by the relation

x

EðxÞ ¼ E1 exp b W

ð67Þ

  E2 E1

ð68Þ

where

b ¼ log

E1 ¼ 300 GPa and E2 ¼ 100 GPa are the Young’s moduli along the left and right edges of the plate, respectively. The mass density q ¼ 7850 kg=m3 and Poisson’s ratio m ¼ 0:25 are constant throughout the specimen. Plane strain conditions are assumed. The specimen is discretized using the coarse mesh shown in Fig. 14b (137 polygons and 392 nodes), and the fine mesh shown in Fig. 14c (346 polygons and 793 nodes). The dilatational wave speed varies throughout the specimen, and has a maximum value of cd ¼ 6:771 mm=ls along the far left edge of the specimen, and a minimum value of cd ¼ 3:9098 mm=ls on the far right edge. h is 13 mm and 7 mm for the

Fig. 16. V-notched problem: (a) geometry and boundary conditions; (b) mouth opening and crack tip nodes at which displacements are calculated; (c) coarse mesh; (d) fine mesh.

I. Chiong et al. / Engineering Fracture Mechanics 131 (2014) 210–231

227

Table 2 Material properties along the left and right edge of the V-notched plate.

Left edge: Iron Right edge: Alumina

E (GPa)

q (kg=m3 )

cd (mm/ls) Eq. (60)

100 350

7150 3800

3.817 9.795

Present (Coarse) Present (Fine) ANSYS

0.5 0.4

M1

u (mm) x

0.3 0.2 0.1

M2 0 −0.1 0

5

10

15

20

25

30

35

40

45

50

Time (μs) Present (Coarse) Present (Fine) ANSYS

0.7 0.6

y

u (mm)

0.5

M1

0.4 0.3 0.2

M2

0.1 0 0

5

10

15

20

25

30

35

40

45

50

Time (μs) Fig. 17. Crack mouth displacement of the V-notched problem: (a) X-displacement; (b) Y-displacement.

coarse and fine mesh, respectively. In order to satisfy the condition in Eq. (58), the time step should not exceed 1.03 ls. The time step of Dt ¼ 0:8 ls is chosen for the analysis. Fig. 15 shows the predicted dynamic SIFs histories of the specimen. The finite element results computed using the M-integral Song and Paulino [16] are used as a reference. The results are in good agreement. To further validate the results of the present formulation, an additional analysis was performed using a finer polygon mesh with 2821 polygons and 5887 nodes. The K I and K II histories computed from this finer mesh showed no appreciable difference to the result presented in Fig. 15, and are not shown. 5.5. V-notched plate The V-notched iron-alumina FGM plate, shown in Fig. 16a, with height 2H ¼ 30 mm, width W ¼ 20 mm and opening angle 2a ¼ 90 is considered. The origin of the coordinate system is at the center of the plate. The plate is clamped along the lower edge. A uniform tension is applied to the top edge according to the Heaviside function PðtÞ ¼ PHðtÞ. Table 2 shows the material properties of the respective individual FGM constituents. The Young’s modulus E and mass density q are assumed to vary exponentially in the x-direction according to

E ¼ E0 exp bE

x

W x

q ¼ q0 exp bq W where E0 ¼ 187:08 GPa and q0 ¼ 5212:5 kg=m3 . The constants bE and bq are

ð69Þ ð70Þ

I. Chiong et al. / Engineering Fracture Mechanics 131 (2014) 210–231

0.3

Present (Coarse) Present (Fine) ANSYS

u (mm) x

0.2

T1

0.1

T2

0

−0.1

−0.2 0

5

10

15

20

25

30

35

40

45

50

35

40

45

50

Time (μs) Present (Coarse) Present (Fine) ANSYS

0.35 0.3

T1

0.2

y

u (mm)

0.25

0.15 0.1 0.05 0

T2 0

5

10

15

20

25

30

Time (μs) Fig. 18. Near tip displacement of the V-notched problem: (a) X-displacement; (b) Y-displacement.

Present (Coarse) Present (Fine) ANSYS

x

σ and τ

xy

(GPa)

15

σ

10

y

5 τ

xy

0 0

5

10 15 20 25 30 35 40 45 50

Time (μs) Present (Coarse) Present (Fine) ANSYS

7

σ and τ (GPa) x xy

228

6 5

σ

y

4 3 2 1

τ

xy

0 0

5

10 15 20 25 30 35 40 45 50

Time (μs) Fig. 19. Stresses,

ry and sxy , of the V-notched problem at: (a) Q(0.9, 0) and; (b) R(0.43431, 0).

I. Chiong et al. / Engineering Fracture Mechanics 131 (2014) 210–231

x 10 12

229

6

Present (Coarse) Present (Fine)

Dynamic SIFs

10 8

K

I

6 4 2

K

II

0 0

5

10

15

20

25

30

35

40

45

50

Time (μs) Fig. 20. Normalized dynamic SIFs for the FGM edge notched example.

bE ¼ ln bq ¼ ln

  EA EI  

qA qI

ð71Þ ð72Þ

where EI and EA are the Young’s modulus of iron and alumina, respectively; and qI and qA are the mass density of iron and alumina. Poisson’s ratio m ¼ 0:2 is assumed to be constant throughout the plate. Plane stress conditions are assumed. Two polygon meshes are used to discretize the plate. Fig. 16c shows the coarse mesh (116 polygons and 330 nodes), and Fig. 16d shows the fine mesh (378 polygons and 859 nodes). To validate the solutions obtained for displacements and stresses, a very fine ANSYS mesh containing 136,906 nodes was used to obtain reference values. The dilatational wave speed varies throughout the specimen, and has a maximum value according to Eq. (60) of cd ¼ 9:795 mm=ls along the far right edge of the specimen. h is 2 mm and 1 mm for the coarse and fine mesh, respectively. To satisfy the condition in Eq. (58), the time step should not exceed 0:1 ls. The time step Dt ¼ 0:05 ls was used in the analyses. 5.5.1. Validation of displacements and stresses The displacement responses at points M1, M2, T1 and T2, shown in Fig. 16b are computed. Figs. 17 and 18 show the time histories of the displacements at the crack mouth nodes (M1 and M2), and two points close to the crack tip (T1 and T2), respectively. The predicted displacements agree very well with the reference solution. The stress components ry and sxy are evaluated at two points: Q ð0:9; 0Þ, and Rð0:43431; 0Þ. Their respective time histories are shown in Fig. 19a for point Q, and Fig. 19b, for point R. It is observed that the present results and the reference solutions are in good agreement. 5.5.2. Dynamic stress intensity factors Next, the dynamic SIFs of the V-notched problem computed using the coarse (Fig. 16c), and fine (Fig. 16d), polygon meshes are presented. These results are shown in Fig. 20. As expected, the time history of K I is similar to the variation of ry near to the crack tip, whereas K II follows a similar trend to that of the shear stresses sxy (see Fig. 19a). To confirm the present result, additional analyses were performed using increasingly finer polygon meshes. No appreciable difference to the present solution was observed, and these results are not shown. 6. Conclusions A polygon based SBFEM formulation for elastodynamic analysis of FGMs has been developed. The displacement field in a polygon is approximated by scaled boundary shape functions. For uncracked polygons, the scaled boundary shape functions are linearly complete. Information on any kind of stress singularity is automatically included in the shape functions when modeling a crack or a notch. This enables accurate calculation of dynamic SIFs directly from their definitions. The formulation does not require local mesh refinement like the FEM, or singular enrichment functions like the XFEM. By using the scaled boundary shape functions, the stiffness and mass matrices in each polygon can be formulated using standard finite element procedures. To model the material heterogeneity in a typical FGM, it is assumed that the material gradients vary locally within each polygon as smooth polynomial functions. These functions are determined from a least squares fit sampled at the Gaussian/Gauss–Lobatto integration points. This assumption lead to semi-analytical expressions for the stiffness and mass matrices of each polygon in the form of matrix power functions and can be integrated analytically in the radial n direction.

230

I. Chiong et al. / Engineering Fracture Mechanics 131 (2014) 210–231

The efficiency of the developed formulation in calculating the dynamic SIFs in FGMs was demonstrated by modeling five numerical benchmarks. It was found that the results of the developed formulation compared well with the numerical results reported in the literature. The numerical examples indicate that present formulation is capable of achieving similar accuracy to the FEM using coarser meshes with significantly fewer DOFs. The results of this study indicate that the scaled boundary polygon formulation is an attractive alternative for the elastodynamic fracture analyses of FGMs.

References [1] Wang BL, Han JC, Du SY. Multiple crack problem in nonhomogeneous composite materials subjected to dynamic anti-plane shearing. Int J Fract 2000;100(4):343–53. [2] Li C, Weng GJ, Brunswick N, Duan Z. Dynamic stress intensity factor of a functionally graded material under antiplane shear loading. Acta Mech 2001;149(1–4):1–10. [3] Guo L, Wu L-Z, Zeng T, Ma L. Fracture analysis of a functionally graded coating-substrate structure with a crack perpendicular to the interface – part II: transient problem. Int J Fract 2004;127(2004):39–59. [4] Ma L, Wu L-Z, Guo L-C, Zhou Z-G. Dynamic behavior of a finite crack in the functionally graded materials. Mech Mater 2005;37(11):1153–65. [5] Ma L, Li J, Abdelmoula R, Wu L-Z. Dynamic stress intensity factor for cracked functionally graded orthotropic medium under time-harmonic loading. Eur J Mech A Solids 2007;26(2):325–36. [6] Monfared M, Ayatollahi M. Dynamic stress intensity factors of multiple cracks in an orthotropic strip with FGM coating. Engng Fract Mech 2013;109:45–57. [7] Bai X, Guo L, Wang Z, Zhong S. A dynamic piecewise-exponential model for transient crack problems of functionally graded materials with arbitrary mechanical properties. Theor Appl Fract Mech 2013;66:41–61. [8] Ayatollahi M, Bagheri R. Dynamic behaviour of several cracks in functionally graded strip subjected to anti-plane time-harmonic concentrated loads. Acta Mech Solida Sinca 2013;26:691–705. [9] Ding S, Li X. The fracture analysis of an arbitrarily orientated crack in the functionally graded material under in-plane impact loading. Theor Appl Fract Mech 2013;66:26–32. [10] Zhang H, Ma G. Fracture modelling of isotropic functionally graded materials by the numerical manifold method. Engng Anal Bound Elements 2014;38:61–71. [11] Barsoum RS. Trianglar quarter-point elements as elastic and perfectly-plastic crack tip elements. Int J Numer Meth Engng 1977;11:85–98. [12] Karihaloo BL, Xiao QZ. Accurate determination of the coefficients of elastic crack tip asymptotic field by a hybrid crack element with p-adaptivity. Eng Fract Mech 2001;68:1609–30. [13] Rice JR. A path independent integral and the approximate analysis of strain concentration by notches and cracks. J Appl Mech 1968;35:379–86. [14] Freund LB. Stress intensity factor calculations based on a conservation integral. Int J Solids Struct 1978;14:241–50. [15] Wu C, He P, Li Z. Extension of J integral to dynamic fracture of functional graded material and numerical analysis. Comput Struct 2002;80(5–6):411–6. [16] Song SH, Paulino GH. Dynamic stress intensity factors for homogeneous and smoothly heterogeneous materials using the interaction integral method. Int J Solids Struct 2006;43(16):4830–66. [17] Sukumar N, Moes N, Moran B, Belytschko T. Extended finite element method for three-dimensional crack modelling. Int J Numer Meth Engng 2000;48:1549–70. [18] Motamedi D, Mohammadi S. Dynamic analysis of fixed cracks in composites by the extended finite element method. Engng Fract Mech 2010;77(17):3373–93. [19] Singh IV, Mishra BK, Bhattacharya S. XFEM simulation of cracks, holes and inclusions in functionally graded materials. Int J Mech Mater Des 2011;7(3):199–218. [20] Bayesteh H, Mohammadi S. XFEM fracture analysis of orthotropic functionally graded materials. Compos Part B: Engng 2013;44(1):8–25. [21] Liu P, Yu T, Bui TQ, Zhang C. Transient dynamic crack analysis in non-homogeneous functionally graded piezoelectric materials by the X-FEM. Comput Mater Sci 2013;69:542–58. [22] Belytschko T, Lu Y, Gu L, Tabbara M. Element-free Galerkin methods for static and dynamic fracture. Int J Solids Struct 1995;32(17/18):2547–70. [23] Zhang C, Sladek J, Sladek V. Effects of material gradients on transient dynamic mode-III stress intensity factors in a FGM. Int J Solids Struct 2003;40(20):5251–70. [24] Sladek J, Sladek V, Zhang C, Solek P, Pan E. Evaluation of fracture parameters in continuously nonhomogeneous piezoelectric solids. Int J Fract 2007;145(4):313–26. [25] Banerjee P, Butterfield R. Boundary element methods in engineering science; 1981. [26] Dominguez J, Gallego R. Time domain boundary element method for dynamic stress intensity factor computations. Int J Numer Meth Engng 1992;33(3):635–47. [27] Albuquerque E, Sollero P, Aliabadi M. The boundary element method applied to time dependent problems in anisotropic materials. Int J Solids Struct 2002;39(5):1405–22. [28] Israil aSM, Banerjee P. Two-dimensional transient wave-propagation problems by time-domain BEM. Int J Solids Struct 1990;26(8):851–64. [29] Fedelinski P, Aliabadi M, Rooke DP. A single-region time domain BEM for dynamic crack problems. Int J Solids Struct 1995;32:3555–71. [30] Zhang C, Savaidis a, Savaidis G, Zhu H. Transient dynamic analysis of a cracked functionally graded material by a BIEM. Comput Mater Sci 2003;26:167–74. [31] Sladek J, Sladek V, Zhang C. An advanced numerical method for computing elastodynamic fracture parameters in functionally graded materials. Comput Mater Sci 2005;32(3–4):532–43. [32] Song C, Wolf JP. The scaled boundary finite-element method – alias consistent infinitesimal finite-element cell method – for elastodynamics. Comput Methods Appl Mech Engng 1997;147:329–55. [33] Deeks AJ, Wolf JP. Stress recovery and error estimation for the scaled boundary finite-element method. Int J Numer Meth Engng 2002;54(4):557–83. [34] Deeks AJ, Wolf JP. An h-hierarchical adaptive procedure for the scaled boundary finite-element method. Int J Numer Meth Engng 2002;54(4):585–605. [35] Song C, Wolf JP. Semi-analytical representation of stress singularities as occurring in cracks in anisotropic multi-materials with the scaled boundary finite-element method. Comput Struct 2002;80(2):183–97. [36] Song C. A matrix function solution for the scaled boundary finite-element equation in statics. Comput Methods Appl Mech Engng 2004;193(23– 26):2325–56. [37] Song C. Evaluation of power-logarithmic singularities, T-stresses and higher order terms of in-plane singular stress fields at cracks and multi-material corners. Engng Fract Mech 2005;72(10):1498–530. [38] Song C, Tin-Loi F, Gao W. Transient dynamic analysis of interface cracks in anisotropic bimaterials by the scaled boundary finite-element method. Int J Solids Struct 2010;47:978–89. [39] Bird G, Trevelyan J, Augarde C. A coupled BEM/scaled boundary FEM formulation for accurate computations in linear elastic fracture mechanics. Engng Anal Boundary Elem 2010;34(6):599–610. [40] Song C. A super-element for crack analysis in the time domain. Int J Numer Meth Engng 2004;61(8):1332–57.

I. Chiong et al. / Engineering Fracture Mechanics 131 (2014) 210–231

231

[41] Chidgzey SR, Deeks AJ. Determination of coefficients of crack tip asymptotic fields using the scaled boundary finite element method. Engng Fract Mech 2005;72(13):2019–36. [42] Mittelstedt C, Becker W. Semi-analytical computation of 3D stress singularities in linear elasticity. Commun Numer Methods Engng 2005;21(5):247–57. [43] Song C, Vrcelj Z. Evaluation of dynamic stress intensity factors and T-stress using the scaled boundary finite-element method. Engng Fract Mech 2008;75(8):1960–80. [44] Song C, Tin-Loi F, Gao W. A definition and evaluation procedure of generalized stress intensity factors at cracks and multi-material wedges. Engng Fract Mech 2010;77(12):2316–36. [45] Goswami S, Becker W. Computation of 3-D stress singularities for multiple cracks and crack intersections by the scaled boundary finite element method. Int J Fract 2012;175(1):13–25. [46] Song C, Wolf JP. The scaled boundary finite-element method: analytical solution in frequency domain. Comput Methods Appl Mech Engng 1998;164(1–2):249–64. [47] Yang Z, Deeks A, Hao H. Transient dynamic fracture analysis using scaled boundary finite element method: a frequency-domain approach. Engng Fract Mech 2007;74(5):669–87. [48] Song C. The scaled boundary finite element method in structural dynamics. Int J Numer Meth Engng 2009;77(8):1139–71. [49] Song C, Tin-Loi F, Gao W. Transient dynamic analysis of interface cracks in anisotropic bimaterials by the scaled boundary finite-element method. Int J Solids Struct 2010;47(7–8):978–89. [50] Chiong I, Sun Z, Xiang T, Song C. Evaluation of dynamic generalised stress intensity factors at cracks and multi-material wedges using the scaled boundary finite element method. Aust J Struct Engng 2012;12(3):197–210. [51] Chiong I, Ooi ET, Song C, Tin-Loi F. Scaled boundary polygons with application to fracture analysis of functionally graded materials. Int J Numer Meth Engng 2014;98:562–89. [52] Ooi ET, Song C, Tin-loi F, Yang Z. Polygon scaled boundary finite elements for crack propagation modelling. Int J Numer Meth Engng 2012;91:319–42. [53] Sukumar N, Tabarraei A. Conforming polygonal finite elements. Int J Numer Meth Engng 2004;61(12):2045–66. [54] Deeks A, Wolf JP. A virtual work derivation of the scaled boundary finite-element method for elastostatics. Comput Mech 2002;28(6):489–504. [55] Newmark N. A method of computation for structural dynamics. ASCE J Engng Mech 1959;85:67–94.