JOURNAL OF THEORETICAL

AND APPLIED MECHANICS

44, 3, pp. 603-636, Warsaw 2006

EVALUATION OF FRACTURE PARAMETERS FOR CRACK

PROBLEMS IN FGM BY A MESHLESS METHOD

Jan Sladek

Vladimir Sladek

Institute of Construction and Architecture, Slovak Academy of Sciences, Bratislava, Slovakia

e-mail: usarslad@savba.sk

Chuanzeng Zhang

Department of Civil Engineering, University of Siegen, Germany

Choon-Lai Tan

Department of Mechanical and Aerospace Engineering, Carleton University, Ottawa, Canada

Ameshless method based on the local Petrov-Galerkin approach is proposed
for crack analysis in two-dimensional (2D), anisotropic and linear elastic so-
lids with continuously varyingmaterial properties. Both quasi-static thermal
and transient elastodynamic problems are considered. For time-dependent
problems, the Laplace transform technique is utilized. The analyzed doma-
in is divided into small subdomains of circular shapes. A unit step function
is used as the test function in the local weak form. It leads to Local Inte-
gral Equations (LIE) involving a domain-integral only in the case of transient
dynamic problems. TheMoving Least Squares (MLS) method is adopted for
approximating the physical quantities in theLIE.Efficient numericalmethods
are presented to compute the fracture parameters, namely, the stress intensity
factors and the T-stress, for a crack inFunctionallyGradedMaterials (FGM).
The path-independent integral representations for stress intensity factors and
T-stresses in continuously non-homogeneous FGM are presented.

Key words: stress intensity factors, T-stress,meshlessmethods, thermoelasti-
city, nonhomogenity, orthotropic materials

1. Introduction

Functionally GradedMaterials (FGM) possess continuously nonhomogeneous
material properties. These materials have been introduced in recent years to



604 J. Sladek et al.

benefit from the ideal performance of their constituents, e.g. high heat and
corrosion resistances of ceramics on one side, and large mechanical strength
and toughness of metals on the other side. In FGM, the composition and the
volume fraction of their constituents vary continuously with spatial coordi-
nates. A review on various aspects of FGM can be found in the monograph
by Suresh and Mortensen (1998) and the review chapter by Paulino et al.
(2003). FGM may exhibit isotropic or anisotropic material properties depen-
ding on the processing technique and the practical engineering requirements.
In the present paper, anisotropic material properties of FGM are considered
since isotropic FGM have been investigated by several authors previously, see
e.g., Erdogan (1995), Sladek et al. (2000), Rao andRahman (2003), Kim and
Paulino (2002), Dolbow andGosz (2002), and Yue et al. (2003).
Cracking of structures is an important phenomenon in many engineering

applications. The use of linear elastic fracture mechanics to predict the crack
behaviour in elastic solids underdifferent loading conditions iswell established
in engineering applications. In conventional linear elastic fracture mechanics
analysis in two dimensions, the primary focus is on the Stress IntensityFactors
(SIF), KI and KII, which characterise the near-tip stress fields. A number of
studies performed over the years have suggested the need to consider leading
non-singular term in the Williams’ series expansion in order to offer better
predictions of the crack path direction and the stability (Cotterell and Rice,
1980), and fracture toughness in elastic solids under conditions of low crack
tip stress triaxiality (Williams and Ewing, 1972; Sumpter, 1993; Betegon and
Hancock, 1991; Sumpter andHancock, 1991). This leading non-singular term
is often referred to in the literature as the T-stress. It represents the stress
acting parallel to the crack plane. Both crack parameters, the stress intensity
factor and the T-stress, are important for a more precise characterisation of
crack fields in the crack tip vicinity. Accurate computationalmethods for their
evaluation are thus needed. Many efficient techniques for evaluation of these
fracture parameters in homogeneous solids exist in the literature. However,
due to the highmathematical complexity of the boundary or initial-boundary
value problem, most investigations on cracked FGM known in the literature
are restricted to isotropic materials (Noda and Jin, 1993a,b; Jin and Noda,
1993a,b; Nemat-Alla and Noda, 1996). Gu and Asaro (1997) studied ortho-
tropic FGM considering a four-point bending specimen with varying Young’s
modulus and varying Poisson’s ratio. Ozturk and Erdogan (1997, 1999) used
the singular integral equation method to investigate mode I andmixed-mode
crack problems in an infinite nonhomogeneous orthotropic medium, with a
crack aligned to one of the principal material axes and a constant Poisson’s
ratio. Recently, Kim and Paulino (2003b) have proposed the domain interac-



Evaluation of fracture parameters for... 605

tion integralmethod for the evaluationof stress intensity factors inorthotropic,
functionally graded materials. Less attention is devoted to the evaluation of
the T-stress. Larsson and Carlsson (1973), Leevers and Radon (1982), Kfo-
uri (1986), Sherry et al. (1995), Nakamura and Parks (1992), Olsen (1994),
J. andV. Sladek (1997a,b), Sladek et al. (1997), Smith et al. (2001), Kim and
Paulino (2003a), and Shah et al. (2005) are among the authorswho have inve-
stigated the T-stress but their works were restricted to homogeneous bodies.
Recently Kim and Paulino (2004) presented a computational method for the
evaluation of the T-stress in orthotropic, functionally gradedmaterials under
a static load. The method is based on the interaction integral expressed by a
domain integral. It requires an accurate evaluation of quantities occurring in
the domain integral at the crack tip vicinity. It is well known that the accura-
cy of the computed stress and displacement fields at the crack tip vicinity is
lower than for those located far away from the crack tip. Therefore, a contour-
domain integral formulation for the evaluation of stress intensity factors and
T-stresses in orthotropic FGMunder thermal andmechanical impact loads is
proposed in this paper.Only the inertial termand the termwith the gradients
of material parameters appear in the domain integral.
The solution of the boundary or initial boundary value problems for conti-

nuously nonhomogeneous solids requires advanced numerical methods due to
the highmathematical complexity. Beside the well-established Finite Element
Method (FEM), the Boundary Element Method (BEM) provides an efficient
and popular alternative to the FEM for solving certain classes of boundary
or initial boundary value problems (Aliabadi, 1997). The conventional BEM
is accurate and efficient for many engineering problems. However, it requires
the availability of the fundamental solutions orGreen’s functions to the gover-
ning equations. Material anisotropy increases the number of elastic constants
in Hooke’s law, and hence makes the construction of the fundamental solu-
tions cumbersome. For 2D elastostatic problems in homogeneous, anisotropic
and linear elastic solids, the fundamental solution is available in closed form
(Eshelby et al., 1953; Schclar, 1994) and it is given in a complex variable
space. Closed-form elastostatic fundamental solutions for 3D anisotropic ela-
sticity exist only for special cases, such as for transversally isotropic or cubic
materials (Ding et al., 1997). In contrast to the static case, very few appli-
cations of the BEM to elastodynamic problems in homogeneous, anisotropic
and linear elastic solids can be found in the literature (Wang andAchenbach,
1996; Albuquerque et al., 2002, 2004; Kgl andGaul, 2000), although theBEM
has been successfully applied to elastodynamic problems in homogeneous, iso-
tropic and linear elastic solids for many years. The main reason lies in the
elastodynamic fundamental solutions for anisotropic and linear elastic solids,



606 J. Sladek et al.

which cannot be given in simple and closed forms and thus make their nu-
merical implementation somewhat cumbersome. Time-domain elastodynamic
fundamental solutions for 2D anisotropic and linear elastic solids have been
applied byWang et al. (1996) to transientwave scattering analysis by a cavity.
The dual reciprocity BEM has also been used by Albuquerque et al. (2002)
andKögl andGaul (2000), where the corresponding elastostatic fundamental
solutions have been utilized.

In this paper, a new computational method is presented to analyze the
boundaryvalue problems in anisotropic FGMwith cracks. The governing equ-
ations for nonhomogeneous, anisotropic and linear elastic solids aremore com-
plex than those for the homogeneous counterpart. The material nonhomoge-
neity gives rise to an additional complication in the derivation of elastostatic
and elastodynamic fundamental solutions. For general nonhomogeneous, ani-
sotropic and linear elastic solids, elastostatic and elastodynamic fundamental
solutions are, hitherto, to the best of the authors’ knowledge, not available.
One possibility to obtain aBEM formulation is based on the use of fundamen-
tal solutions for a fictitious homogeneous medium (Sladek et al., 1993). This
approach, however, leads to a boundary-domain integral formulation with a
domain-integral containing the gradients of the primary fields. Owing to the
singularities, special care is required in the numerical integration when the
standard boundary-domain formulation is employed. Moreover, the system
matrix is relatively large and fully populated. To overcome such difficulties, a
local integral formulation can be used for general nonhomogeneous solids (Sla-
dek et al., 2000;Mikhailov, 2002). The application of Local Integral Equations
(LIE) requires the use of a domain approximation of the physical fields in the
numerical implementation.

Meshless formulations have become increasingly popular in recent years
due to their higher adaptivity and lower cost for preparing input data in the
numerical analysis. Several meshless methods have already been proposed in
literature (Belytschko et al., 1994; Atluri and Shen, 2002; Atluri, 2004). Many
of them are derived from a weak form formulation on a global domain or a
set of local subdomains. The global formulation requires background cells for
the integration of the weak form. In contrast, the local weak form formulation
needs no cells and therefore, the correspondingmethods are often called truly
meshless methods. If a simple form for the geometry of the subdomains is
chosen, numerical integrations over them can be easily carried out. The Me-
shless Local Petrov-Galerkin (MLPG) method is a fundamental base for the
derivation ofmanymeshless formulations, since trial and test functions can be
chosen from different functional spaces. If a unit step function is used as the



Evaluation of fracture parameters for... 607

test function in the local weak form to derive LIE, the form of LIE is much
simpler than that provided by utilizing the singular fundamental solutions.
Such an approach has been recently applied to problems in homogeneous, ani-
sotropic and linear elastic solids by Sladek et al. (2004). It is extended in this
paper to continuously nonhomogeneous, anisotropic and linear elastic solids.
It yields a pure contour or boundary integral formulation on local boundaries
for static problems in anisotropic linear elasticity. In anisotropic elastodyna-
mics an additional domain-integral containing the inertial terms is involved.
The Laplace transform is applied to eliminate the time variable in the gover-
ning equations and theboundary conditions of elastodynamic problems.Then,
the local boundary integral equations are derived in the Laplace transformed
domain. Several quasi-static boundary value problems have to be solved for
various values of the Laplace transformparameter. The integral equations ha-
ve a very simple nonsingular form. Moreover, both the contour and domain
integrations can be easily carried out on circular subdomains. The Stehfest’s
inversion method (Stehfest, 1970) is applied to obtain the time-dependent
solutions. The spatial variation of the displacements is approximated by the
moving least-squares (MLS) scheme. Several numerical examples for crack pro-
blems in nonhomogeneous orthotropic and linear elastic solids are presented
and discussed.

2. Evaluation of the first fracture parameter – the SIF

In stationary thermoelasticity for non-homogeneous anisotropic bodies, the
thermal and the stress fields aredescribedby the following governing equations
(Nowacki, 1975)

(λijθ,j ),i=−Q σij,j +Xi =0 (2.1)

where θ is the temperature measured in the scale with its origin at the equ-
ilibrium state, Xi and Q are the body force vector and the heat source re-
spectively, and σij is the stress tensor. The thermal conductivity tensor λij
is a continuous function of Cartesian coordinates. A subscript preceded by
a comma denotes differentiation with respect to the corresponding Cartesian
coordinate.
The strain tensor εij is given by the displacement gradients as follows

εij =
1

2
(ui,j +uj,i) (2.2)



608 J. Sladek et al.

In the presence of a temperature gradient, the total strain tensor can be de-
composed into its elastic part εeij and another one accounting for the free
thermal expansion of the medium. Thus,

εij = ε
e
ij +αijθ (2.3)

where αij is the tensor of thermal expansions.

According to the Neumann hypothesis, the stress tensor is related to the
elastic strain in the usual way, viz

σij = cijklε
e
kl (2.4)

where cijkl is the elasticity tensor. For an isotropic continuum, it is given by

cijkl = µ
( 2ν
1−2ν

δijδkl + δikδjl + δilδjk
)

(2.5)

where µ is the shearmodulus and ν is Poisson’s ratio. Note that thematerial
is non-homogeneous in general, since the elasticity tensor and the coefficients
of thermal expansions depend on the spatial position.

The elastic strain energy density W = W(εij,xi) can be written in the
following form

W =
1

2
cijkl(εij −αijθ)(εkl −αklθ) (2.6)

Then, the gradient of strain energy density is given as

W,m=
∂W

∂εij

∂εij
∂xm
+
(∂W
∂xm

)

exp l
=
∂W

∂εeij

∂εij
∂xm
+
(∂W
∂xm

)

exp l
(2.7)

where the term for the ”explicit” derivative of the strain energy density for
non-homogeneous materials becomes

(∂W
∂xm

)

exp l
=
1

2
cijkl,mε

e
ijε

e
kl −σklαkl,mθ−σklαklθ,m (2.8)

Utilizing Eqs. (2.2) and (2.4), the gradient of the strain energy density can be
rewritten in the form

W,m=(σijui,m),j −σij,jui,m +(W,m )exp l (2.9)

Then, from Eqs. (2.9) and (2.1)2, it follows that

(Wδjm −σijui,m),j = Xiui,m +(W,m )exp l (2.10)



Evaluation of fracture parameters for... 609

An integral form of Eq. (2.10) may be obtained upon application of the di-
vergence theorem. If Ω is a regular bounded region enclosed by a surface Γ
whose unit outward normal vector is n, it follows that

∫

Γ

(Wδjm −σijui,m)nj dΓ =
∫

Ω

Xiui,m dΩ +

∫

Ω

(W,m )exp l dΩ (2.11)

The integral identity (2.11) is valid in a region where no field irregularities
prevail. In a crack problem, the stresses at the crack tip are singular and the
displacements are discontinuous across both crack-surfaces. Therefore, a cut-
off along the crack with a small region in the vicinity of a crack tip Ωε has to
be excluded. This region is surrounded by Γε as shown in Fig.1. The global
Cartesian coordinate system is defined in such away that the principal axes of
orthotropy are aligned with the global coordinates. Both fields σij and ui are
regular in the region Ω−Ωε. The contour Γ = Γ0+Γ+c −Γε+Γ−c is a closed
integration path in the counter-clockwise direction. The radius ε is considered
to be very small and shrunk to zero in the limiting process. The crack surfaces
Γ+c and Γ

−
c are assumed to be free of tractions, i.e., ti = σijnj =0, and the

crack is parallel to the x1-axis of the local Cartesian coordinate system.Then,
Eq. (2.11) can be written as

∫

Γε

(Wδjm −σijui,m)nj dΓ =
∫

Γ0

(Wδjm −σijui,m)nj dΓ +

(2.12)

+

∫

Γ
+
c

(W+−W−)δ2m dΓ − lim
ε→0

∫

Ω−Ωε

Xiui,m dΩ − lim
ε→0

∫

Ω−Ωε

(W,m )exp l dΩ

The left-hand side of Eq. (2.12) is identical to the definition of the Ĵ-integral
(Kishimoto et al., 1980) for m =1, which has the following form

Ĵ =

∫

Γ0

(Wn1−σijnjui,1) dΓ − lim
ε→0

∫

Ω−Ωε

Xiui,1 dΩ − lim
ε→0

∫

Ω−Ωε

(W,1 )exp l dΩ

(2.13)
Consider two independent equilibrium states in an orthotropic, functional-

ly graded material. Let the state 1 correspond to the actual state for given
boundary conditions, and let the state 2 correspond to an auxiliary state,
which can be either the mode I or the mode II near-tip displacement and
stress fields. The asymptotic distribution of the stress and the displacement
fields in the vicinity of the crack tip in a continuously nonhomogeneous me-
dium is the same as in the homogeneous one.We select the auxiliary stress and



610 J. Sladek et al.

Fig. 1. Integration paths and coordinate definitions

displacement fields, denoted by superscript (2), as the crack tip asymptotic
fields in the orthotropic elastic solid given by (Sih et al., 1965)

σ
(2)
11 =

K
(2)
I√
2πr
Re
[ µtip1 µ

tip
2

µ
tip
1 −µ

tip
2

(µtip2
A2
− µ

tip
1

A1

)]
+

+
K
(2)
II√
2πr
Re
[ 1
µ
tip
1 −µ

tip
2

((µtip2 )
2

A2
−
(µ
tip
1 )
2

A1

)]

σ
(2)
22 =

K
(2)
I√
2πr
Re
[ 1
µ
tip
1 −µ

tip
2

(µtip1
A2
−
µ
tip
2

A1

)]
+

+
K
(2)
II√
2πr
Re
[ 1
µ
tip
1 −µ

tip
2

( 1
A2
− 1
A1

)]

σ
(2)
12 =

K
(2)
I√
2πr
Re
[ µtip1 µ

tip
2

µ
tip
1 −µ

tip
2

( 1
A1
−
1

A2

)]
+

+
K
(2)
II√
2πr
Re
[ 1
µ
tip
1 −µ

tip
2

(µtip1
A1
− µ

tip
2

A2

)]
(2.14)

u
(2)
1 = K

(2)
I

√
2r

π
Re
[ 1
µ
tip
1 −µ

tip
2

(µ
tip
1 P12A2−µ

tip
2 P11A1

)]
+

+K
(2)
II

√
2r

π
Re
[ 1
µ
tip
1 −µ

tip
2

(P12A2−P11A1
)]



Evaluation of fracture parameters for... 611

u
(2)
2 = K

(2)
I

√
2r

π
Re
[ 1
µ
tip
1 −µ

tip
2

(µ
tip
1 P22A2−µ

tip
2 P21A1

)]
+

+K
(2)
II

√
2r

π
Re
[ 1
µ
tip
1 −µ

tip
2

(P22A2−P21A1
)]

and

A1 =
√
cosϕ+µ

tip
1 sinϕ A2 =

√
cosϕ+µ

tip
2 sinϕ

where (r,ϕ) are polar coordinates with the origin at the crack tip and related
to the local Cartesian coordinate system (x1,x2), Re denotes the real part of
a complex function, µ

tip
i are material parameters at the crack tip, which are

roots of the following characteristic equation in termsof the elastic compliances
βmn (m,n =1,2 and 6) of the anisotropic material (Lekhnitskii, 1963)

β11µ
4−2β16µ3+(2β12+β66)µ2−2β26µ+β22 =0 (2.15)

and

Pik =

[
β11µ

2
k +β12−β16µk

β12µk +β22/µk −β26

]
(2.16)

The auxiliary stress fields in Eq. (2.14) are in equilibrium in absence of the

inertial effects, i.e., σ
(2)
ij,j =0. The auxiliary strain field is chosen as

ε
(2)
ij = Sijkl(x)σ

(2)
kl (2.17)

which differs from S
tip
ijklσ

(2)
kl , where Sijkl(x) is the compliance tensor of the

actual FGM.Thus, the auxiliary strain field inEq. (2.17) is incompatible with
the symmetric part of the auxiliary displacement gradients (Kim andPaulino,
2003b)

ε
(2)
ij 6=

1

2
(u
(2)
i,j +u

(2)
j,i )

Although this incompatibility of the strainfieldvanishes as the contour shrinks
to the crack tip, it gives finite contributions on the contour Γ0 and the domain
Ω −Ωε.
Superposition of the actual and the auxiliary fields leads to another equ-

ilibrium state (state s) for which the J-integral is given as



612 J. Sladek et al.

J(s) =

∫

Γ0

[1
2
(σij +σ

(2)
ij )(εij +ε

(2)
ij )n1− (σij +σ

(2)
ij )nj(ui,1+u

(2)
i,1)
]
dΓ +

− lim
ε→0

∫

Ω−Ωε

Xi(ui,1+u
(2)
i,1) dΩ − lim

ε→0

∫

Ω−Ωε

1

2
Cijkl,1(εkl +ε

(2)
kl )(εij +ε

(2)
ij ) dΩ +

− lim
ε→0

∫

Ω−Ωε

(σkl +σ
(2)
kl )(αkl,1θ+αklθ,1 ) dΩ (2.18)

which is conveniently decomposed into

J(s) = J +J(2)+M (2.19)

where

J(2) =

∫

Γ0

[1
2
σ
(2)
ij ε
(2)
ij n1−σ

(2)
ij nju

(2)
i,1

]
dΓ − lim

ε→0

∫

Ω−Ωε

1

2
Cijkl,1ε

(2)
kl ε
(2)
ij dΩ (2.20)

The interaction integral M is then given by

M =

∫

Γ0

[1
2
(σijε

(2)
ij +σ

(2)
ij εij)n1− (σijnju

(2)
i,1 +σ

(2)
ij njui,1)

]
dΓ +

− lim
ε→0

∫

Ω−Ωε

Xiu
(2)
i,1 dΩ + (2.21)

− lim
ε→0

∫

Ω−Ωε

[1
2
Cijkl,1(εklε

(2)
ij +ε

(2)
kl εij)−σ

(2)
kl (αkl,1θ+αklθ,1 )

]
dΩ

Using the following identities

σijε
(2)
ij = σijSijkl(x)σ

(2)
kl = εklσ

(2)
kl

(2.22)

Cijkl,1εklε
(2)
ij = Cklij,1εklε

(2)
ij = Cijkl,1εijε

(2)
kl

the M integral can be rewritten in the form

M =

∫

Γ0

[σijε
(2)
ij n1−σijnju

(2)
i,1 −σ

(2)
ij njui,1] dΓ − lim

ε→0

∫

Ω−Ωε

Xiu
(2)
i,1 dΩ +

(2.23)

− lim
ε→0

∫

Ω−Ωε

[Cijkl,1εklε
(2)
ij −σ

(2)
kl (αkl,1θ+αklθ,1 )] dΩ



Evaluation of fracture parameters for... 613

For mixed-mode crack problems, the relationship of J-integral and stress in-
tensity factors is given by Kim and Paulino (2003b)

J = c11K
2
I + c12KIKII + c22K

2
II (2.24)

where

c11 =−
a
tip
22

2
Im
(µtip1 +µ

tip
2

µ
tip
1 µ

tip
2

)

c12 =−
a
tip
22

2
Im
( 1
µ
tip
1 µ

tip
2

)+
a
tip
22

2
Im(µ

tip
1 µ

tip
2 ) (2.25)

c22 =
a
tip
11

2
Im(µ

tip
1 +µ

tip
2 )

For two admissible fields (actual and auxiliary) one obtains

J(s) = c11(KI +K
(2)
I )
2+ c12(KI +K

(2)
I )(KII +K

(2)
II )+ c22(KII +K

(2)
II )
2 =

= J +J(2)+M

where
J(2) = c11(K

(2)
I )
2+ c12K

(2)
I K

(2)
II + c22(K

(2)
II )
2

and

M =2c11KIK
(2)
I + c12(KIK

(2)
II +K

(2)
I KII)+2c22KIIK

(2)
II (2.26)

The mode I and the mode II stress intensity factors are evaluated by solving
the system of linear algebraic equations

2c11KI + c12KII = M
I c12KI +2c22KII = M

II (2.27)

resulting from Eq. (2.26) by taking K
(2)
I = 1, K

(2)
II = 0 for M

I, and

K
(2)
I = 0, K

(2)
II = 1 for M

II, respectively. The values MI and MII are
computed numerically by Eq. (2.23) with an adequate choice of the auxiliary
solutions according to Eq. (2.14).
The advantage of the presented method for the computation of stress in-

tensity factors through the M-integral technique is high accuracy, because the
contour integral is evaluated at points far away from the crack tip. Numerical
errors of computed quantities at the crack tip vicinity in a direct computation
of SIF from the asymptotic expansion formulae are hence reduced in this ap-
proach. The results are also insensitive to the distance of the evaluation point
from the crack tip.



614 J. Sladek et al.

Similarly, one can derive the M-integral representation for elastodynamic
problems in continuously nonhomogeneous orthotropic materials. The M in-
tegral in this case has the following form

M =

∫

Γ0

(σijε
(2)
ij n1−σijnju

(2)
i,1 −σ

(2)
ij njui,1) dΓ − lim

ε→0

∫

Ω−Ωε

Xiu
(2)
i,1 dΩ +

(2.28)

− lim
ε→0

∫

Ω−Ωε

(Cijkl,1εklε
(2)
ij −ρüiu

(2)
i,1) dΩ

It should be remarked that the terms in Eq. (2.23) or (2.28) are given in the
local coordinate system with the origin at the crack tip as shown in Fig.1.
Recall that the derived expressions (2.23) and/or (2.28) for the M-integral
are related to the local Cartesian coordinate system, where M1(local) = M.
Due to material orthotropy, it is convenient to perform all the numerical cal-
culations in the global Cartesian coordinate system. Thus, in the numerical
computation, the M-integrals (M1 and M2) are obtained in global coordina-
tes and then transformed according to the law for the transformation of vector
components to get the local component M1(local) or equivalently M.

The transformation is given by Kim and Paulino (2004)

M = M1(local) = M1(global) cosω+M2(global) sinω (2.29)

where

Mm(global) =

∫

Γ0

(σijε
(2)
ij nm −σijnju

(2)
i,m −σ

(2)
ij njui,m) dΓ +

(2.30)

− lim
ε→0

∫

Ω−Ωε

Xiu
(2)
i,m dΩ − lim

ε→0

∫

Ω−Ωε

(Cijkl,mεklε
(2)
ij −ρüiu

(2)
i,m) dΩ

and ω is the angle between the x1-axes of the local and the global coordinate
systems.

3. Evaluation of the second fracture parameter – the T-stress

The T-stress is very frequently considered as the second fracture parame-
ter in fracture mechanics analysis. It is the leading non-singular term in the



Evaluation of fracture parameters for... 615

Williams’ (1957) asymptotic expansion of stresses. The T-stress can be com-
puted directly from the above-mentioned asymptotic expansion of stresses or
displacements, if the stress intensity factors and the stress or the displacement
values at nodal points close to the crack tip are obtained from a numerical
analysis. In the literature, such a method where T-stress is computed from
the known stress component σ11, is called the boundary layermethod (Sherry
et al., 1995). A drawback of this method the sensitivity of the results with re-
spect to the distance of the evaluation point from the crack tip. It is therefore
suggested that amethod shpuld be adopted inwhich the T-stress is expressed
in terms of the solution at points far away from the crack tip. In the following,
an integral representation of the T-stress will be derived for a cracked body
analysed in the framework of stationary thermoelasticity.

The stresses and the displacements in the crack tip vicinity of a material
with a continuous non-homogeneity have the same singularity and angular
distributions as those in a homogeneous material. The displacements and the
stresses close to the crack tip can be written in an asymptotic form as (Shah
et al., 2005)

u1(r,ϕ) =
KI
4

√
r

2π
gI1(ϕ,µ

tip
1 ,µ

tip
2 )+

KII
4

√
r

2π
gII1 (ϕ,µ

tip
1 ,µ

tip
2 )+a11Trcosϕ

(3.1)

u2(r,ϕ) =
KI
4

√
r

2π
gI2(ϕ,µ

tip
1 ,µ

tip
2 )+

KII
4

√
r

2π
gII2 (ϕ,µ

tip
1 ,µ

tip
2 )+a12Trsinϕ

σij(r,ϕ) =
KI√
2πr

fIij(ϕ,µ
tip
1 ,µ

tip
2 )+

KII√
2πr

fIIij (ϕ,µ
tip
1 ,µ

tip
2 )+Tδi1δj1 (3.2)

where the functions f
I,II
ij (ϕ,µ

tip
1 ,µ

tip
2 ) and g

I,II
i (ϕ,µ

tip
1 ,µ

tip
2 ) are given inEqs.

(2.14).

An attempt will be made to find an appropriate auxiliary field to obtain
the integral expression of the T-stress. In order to obtain a finite contribution
of this term in the interaction integral M given byEq. (2.23), the stress tensor

σ
(2)
ij should be proportional to r

−1. However, such auxiliary fields multiplied
by the first part of the asymptotic expansion of stresses and displacements,
which contain the stress intensity factors, give the singular integrand on Γε.
These singular terms have to be eliminated only by the angular variation of
the auxiliary fields.

The application of a concentrated point-force f at the crack tip in the
plane of the crack will give the stress and the displacement fields which will
meet the conditions mentioned above. The auxiliary stress field solution due
to a concentrated point-force acting at the crack tip of an anisotropic, ho-



616 J. Sladek et al.

mogeneous, and infinite wedge-shaped beam with a unit thickness is given as
(Lekhnitskii, 1963)

σ(2)rr =
C cosφ+D sinφ

rL(φ)
σ
(2)
φφ = σ

(2)
rφ =0 (3.3)

where

L(φ)= β
tip
11 cos

4φ+(2β
tip
12 +β

tip
66 )sin

2φcos2φ+β
tip
22 sin

4φ

with β
tip
ij being the material compliances at the crack tip, and the constants

C and D may be determined from the equilibrium forces in global coordinate
system

C

ψ2∫

−ψ1

cos2φ

L(φ)
dφ+D

ψ2∫

−ψ1

sinφcosφ

L(φ)
dφ =−f cosω

(3.4)

C

ψ2∫

−ψ1

sinφcosφ

L(φ)
dφ+D

ψ2∫

−ψ1

sin2φ

L(φ)
dφ =−f sinω

InEq. (3.4), the integration bounds ψ1 and ψ2 are givenby the angles between
the X1-axis and the lower and upper crack faces, respectively.

The auxiliary strains are defined by using the actual material compliances
at the considered point of the non-homogeneous medium as

ε(2)rr (x)= β11(x)σ
(2)
rr ε

(2)
φφ(x)= β12(x)σ

(2)
rr (3.5)

The auxiliary stress fields in global coordinate system can be expressed as

σ
(2)
11 = σ

(2)
rr cos

2φ =
C cosφ+D sinφ

rL(φ)
cos2φ

σ
(2)
22 = σ

(2)
rr sin

2φ =
C cosφ+D sinφ

rL(φ)
sin2φ (3.6)

σ
(2)
12 = σ

(2)
rr sinφcosφ =

C cosφ+D sinφ

rL(φ)
sinφcosφ

Thedisplacementderivatives for theauxiliaryfieldderivedbyKimandPaulino
(2004) are



Evaluation of fracture parameters for... 617

u
(2)
1,1 =(β

tip
11 cos

2φ+β
tip
12 sin

2φ)
C cosφ+D sinφ

rL(φ)

u
(2)
2,2 =(β

tip
12 cos

2φ+β
tip
22 sin

2φ)
C cosφ+D sinφ

rL(φ)
(3.7)

u
(2)
1,2 =

CL(φ)sinφ−CH1(φ)−DH2(φ)
rL(φ)

u
(2)
2,1 =

2β
tip
66 (C cosφ+D sinφ)sinφcosφ−CL(φ)sinφ+CH1(φ)+DH2(φ)

rL(φ)

where

H1(φ)= β
tip
11 cos

4φsinφ+β
tip
12 sinφcos

2φ−βtip22 sin
3φcos2φ+

−(βtip12 +β
tip
66 )sinφcos

4φ

H2(φ)= β
tip
11 cos

3φ+β
tip
12 sin

2φcosφ

As it is seen in the previous section, the auxiliary strain field, which follows

fromEq. (3.5), has an incompatibility since ε
(2)
ij 6=(u

(2)
i,j +u

(2)
j,i )/2. The interac-

tion integral M given by Eq. (2.23) is path-independent since it is expressed
in terms of J-integrals. Thus, the contour can be chosen arbitrarily, say, as a
circle with the radius ε, which is shrunk to zero. In such a case one can write

M = lim
ε→0

∫

Γε

(σijε
(2)
ij n1−σijnju

(2)
i,1 −σ

(2)
ij njui,1) dΓ (3.8)

The asymptotic displacements and stresses can be split into singular and non-
singular parts as follows

ui = u
s
i +u

T
i σij = σ

s
ij +σ

T
ij (3.9)

where the terms with superscript s contain the stress intensity factors in Eq.
(3.2) and superscript T is related to the T-stress term

σTij = Tδi1δj1 (3.10)

The elastic strains at the crack tip corresponding to the uniform T-stress in
thermoelasticity are given as

εT11 = β
tip
11 T +α

tip
11θ

tip

(3.11)

uTi,1 = u
T
1,1δi1 = ε

T
11δi1 = β

tip
11 Tδi1+α

tip
11 θ

tipδi1



618 J. Sladek et al.

Thus, the M-integral reduces to

M =− lim
ε→0

∫

Γε

σ
(2)
ij nju

T
i,1 dΓ =−(β

tip
11 T +α

tip
11 θ

tip) lim
ε→0

∫

Γε

σ
(2)
ij nj dΓ =

(3.12)

= (β
tip
11 T +α

tip
11 θ

tip)f

The final result of Eq. (3.12) may be rearranged and take the form

T =
M −αtip11 θtipf

fβ
tip
11

(3.13)

The unknown displacements, traction vector and temperature, required along
the integrationpathandwithin thedomainenclosedby the integration contour
in Eq. (2.30), can be obtained from a numerical or experimental analysis. In
thispaper, themeshless local integral equationmethod isused for thispurpose.

4. Meshless local Petrov-Galerkin method in continuously

non-homogeneous solids

4.1. Uncoupled thermoelasticity

Let us consider a boundary value problem defined in the stationary un-
coupled thermoelasticity for a continuously non-homogeneous anisotropicme-
dium, which in 2D is described by the governing equations (2.1). Since Eqs.
(2.1) are uncoupled, they can be solved separately. In the first step we solve
the heat conduction equation (2.1)1.

Instead of writing the global weak form for the above governing equation,
the MLPG methods construct the weak form over local subdomains such
as Ωs, which is a small region taken for each node inside the global domain
(Atluri and Shen, 2002). The local subdomains overlap each other, and cover
the whole global domain Ω; they could be of any geometric shape and size.
In the present paper, the local subdomains are taken to be of circular shape
(Fig.2). The local weak form of the governing equation (2.1)1 can be written
as ∫

Ωs

[(λij(x)θ,j (x)),i+Q(x)]θ
∗(x) dΩ =0 (4.1)

where θ∗(x) is a weight (test) function.



Evaluation of fracture parameters for... 619

Fig. 2. Local boundaries for weak formulation, the domain Ωx forMLS
approximation of the trial function, and support area of weight function around

node xi

Applying the Gauss divergence theorem to Eq. (4.1) one can write

∫

∂Ωs

q(x)θ∗(x) dΓ −
∫

Ωs

λij(x)θ,j (x)θ,
∗

i (x) dΩ +

∫

Ωs

Q(x)θ∗(x) dΩ =0 (4.2)

where ∂Ωs is the boundary of the local subdomain and

q(x)= λij(x)θ,j (x)ni(x)

The localweak form inEq. (4.2) is a startingpoint to derive the local boundary
integral equation if an appropriate test function is selected. If aHeaviside step
function is chosen as the test function θ∗(x) in each subdomain

θ∗(x)=

{
1 at x∈ Ωs
0 at x /∈ Ωs

the local weak form, Eq. (4.2), is transformed into a simple local integral
equation ∫

∂Ωs

q(x) dΓ =−
∫

Ωs

Q(x) dΩ (4.3)

Equation (4.3) is recognized as the flow balance condition of the subdomain.
In theMLPGmethod, the test and trial functions are not necessarily from the
same functional spaces. For internal nodes, the test function is chosen as the



620 J. Sladek et al.

Heaviside step function with the support on the local subdomain. The trial
function, on the other hand, is chosen to be theMoving Least-Squares (MLS)
interpolation over a number of nodes randomly distributedwithin the domain
of influence. While the local subdomain is defined as the support of the test
function on which the integration is carried out, the domain of influence is
defined as a region where the weight function is not zero, i.e., all nodes lying
inside that region influence the interpolation. The approximate function can
be written as (Atluri and Shen, 2002)

θh(x)=Φ>(x) · θ̂=
n∑

a=1

φa(x)θ̂a (4.4)

where θ̂a are fictitious parameters and φa(x) is the shape function associated
with the node a. The number of nodes n used for the approximation of θ(x)
is determined by the weight function wa(x). A spline-type weight function is
considered in the present work

wa(x)=




1−6

(da

ra

)2
+8
(da

ra

)3
−3
(da

ra

)4
0¬ da ¬ ra

0 da ­ ra
(4.5)

where da = ‖x−xa‖ and ra is the radius of the supportdomain for theweight
function wa. The directional derivatives of θ(x) are approximated in terms of
the same nodal values as

∂θ
h

∂n
(x,p)= nk(x)

n∑

a=1

θ̂a(p)φ,ak (x) (4.6)

Making use of theMLS-approximations (4.4) and (4.6) for θ(x) and

q(x)= λij(x)ni(x)
n∑

a=1

θ̂aφ,aj (x) (4.7)

the local boundary integral equation, Eq. (4.3), for all subdomains yields the
following set of equations

n∑

a=1

θ̂a
∫

∂Ωs

λij(x)nj(x)φ,
a
j (x) dΓ =−

∫

Ωs

Q(x) dΩ (4.8)

It should be noted that neither the Lagrange multipliers nor penalty parame-
ters are introduced into the local weak form in Eq. (4.1), because the essential



Evaluation of fracture parameters for... 621

boundary conditions on Γθ can be imposed directly using the interpolation
approximation in Eq. (4.4)

n∑

a=1

φa(x)θ̂a = θ̃(x) for x∈ Γθ (4.9)

The same procedure to derive the local integral equation for the heat conduc-
tion equation can be applied to obtain the corresponding LIE for mechanical
fields described by Eq. (2.1)2. In the case of elastic materials, the relation
between the stress and the strain is given by Hookes law for an anisotropic
body

σij(x, t)= Cijkl(x)εkl(x, t)= Cijkl(x)uk,l(x, t) (4.10)

where Cijkl is the elasticity tensor which exhibits the symmetries

Cijkl =Cjikl = Cklij

The Cauchy formula has been employed for strains in the last equality in Eq.
(4.10). Then, the traction vector ti = σijnj is given by

ti(x, t)= Cijkl(x)uk,l(x, t)nj(x) (4.11)

where nj denotes a unit outward normal vector.

For the plane stress problem of orthotropic materials, one can write



σ11
σ22
σ12


=D(x)



ε11
ε22
2ε12


 (4.12)

where

D(x)=



E1/e E2ν12/e 0

E2ν12/e E2/e 0
0 0 G12


 with e =1−

E2
E1
(ν12)

2

The local weak form of the governing equation (2.1)2 can be written as

∫

Ωs

[σij,j(x)+Xi(x)]u
∗

ik(x) dΩ =0 (4.13)

where u∗ik(x) is a test function.



622 J. Sladek et al.

Using the same procedure as for the above thermal analysis, and assuming
the test function as

u∗ik(x)=

{
δik at x∈ (Ωs ∪∂Ωs)
0 at x /∈ Ωs

the local weak form, Eq. (4.13), leads to the local boundary integral equations
∫

∂Ωs

ti(x) dΓ +

∫

Ωs

Xi(x) dΩ =0 (4.14)

Note here that a pure contour integral formulation is obtained under the as-
sumption of vanishing body forces. Analogous to the approximation for the
temperature, one can approximate the displacements by

uh(x)=Φ>(x) · û=
n∑

a=1

φa(x)ûa (4.15)

where the nodal values ûa are fictitious nodal parameters for displacements.
The traction vectors ti(x) at a boundarypoint x∈ ∂Ωs are approximated

in terms of the same nodal values ûa as

t
h(x)=N(x)D

n∑

a=1

B
a(x)ûa (4.16)

where the matrix N(x) is related to the normal vector n(x) on ∂Ωs by

N(x)=

[
n1 0 n2
0 n2 n1

]

and the matrix Ba is represented by the gradients of the shape functions as

B
a =



φ,a1 0
0 φ,a2
φ,a2 φ,

a
1




Satisfying the boundary conditions at those nodal points on the global bo-
undary, where displacements are prescribed, and making use of the approxi-
mation (4.15), one obtains the discretized form of the displacement boundary
conditions given as

n∑

a=1

φa(ζ)ûa = ũ(ζ) for ζ ∈ Γu (4.17)



Evaluation of fracture parameters for... 623

With the MLS-approximations, Eqs.(4.15) and (4.16) for the unknown fields
in the local boundary integral equations (4.14), we obtain the discretized LIE

n∑

a=1

û
a
∫

Ls+Γsu

N(x)D(x)Ba(x) dΓ =−
∫

Γst

t̃(x) dΓ −
∫

Ωs

X(x) dΩ (4.18)

which are considered on the sub-domains adjacent to interior nodes as well as
to the boundary nodes on Γst.

4.2. Elastodynamics

Now, we consider a linear elastodynamic problem in an anisotropic, conti-
nuously nonhomogeneous and linear elastic domain Ω bounded by the boun-
dary Γ . The governing equations can be expressed as

σij,j(x, t)−ρ(x)üi(x, t)=−Xi(x, t) (4.19)

where ρ is the mass density, and double dots indicate the second derivati-
ve with respect to time. The following boundary and initial conditions are
assumed

ui(x, t)= ũi(x, t) on Γu

ti(x, t)= t̃i(x, t) on Γt

ui(x, t)
∣∣
t=0
= ui(x,0)

u̇i(x, t)
∣∣
t=0
= u̇i(x,0)



 in Ω

where Γu and Γt are the parts of the global boundarywith prescribed displa-
cements and tractions, respectively.
Applying the Laplace transform to the governing equations, Eq. (4.19), we

obtain
σij,j(x,p)−ρ(x)p2ui(x,p)=−F i(x,p) (4.20)

where p is the Laplace transform parameter, and

F i(x,p)= Xi(x,p)+pui(x,0)+ u̇i(x,0)

The local weak form of Eq. (4.20) can be written as
∫

Ωs

[σij,j(x,p)−ρp2ui(x,p)+F i(x,p)]u∗ik(x) dΩ =0 (4.21)

Using the same procedure as for the uncoupled thermoelasticity discussed
above, one obtains the local integral equations
∫

Ls+Γsu

ti(x,p) dΓ −
∫

Ωs

ρp2ui(x,p) dΩ =−
∫

Γst

t̃i(x,p) dΓ −
∫

Ωs

F i(x,p) dΩ

(4.22)



624 J. Sladek et al.

Equation (4.22) is recognized as the overall equilibrium of forces including the
inertial ones acting on the subdomain. The approximate Laplace transforms
of the displacements can be written as

uh(x,p)=Φ>(x) · û(p)=
n∑

a=1

φa(x)ûa(p) (4.23)

Similarly, the traction vector ti(x,p) at a boundary point x∈ ∂Ωs is appro-
ximated in terms of the same nodal values ûa(p) as

t
h
(x,p)=N(x)D

n∑

a=1

B
a(x)ûa(p) (4.24)

With theMLS approximations (4.23) and (4.24) for the unknownfields in the
local boundary integral equations, Eq. (4.22), we obtain the discretized LIE

n∑

a=1

û
a(p)
( ∫

Ls+Γsu

N(x)D(x)Ba(x) dΓ −ρp2
∫

Ωs

φa(x) dΩ
)
=

(4.25)

=−
∫

Γst

t̃(x,p) dΓ −
∫

Ωs

F(x,p) dΩ

which are considered on the sub-domains adjacent to the interior nodes aswell
as to the boundary nodes on Γst. Collecting the discretized LIE together with
the discretized boundary conditions for displacements, we get the complete
systemof linearalgebraic equations for the computationof thenodalunknowns
whichare theLaplace transformsof fictitiousparameters ûa(p).Theboundary
and domain integrals are evaluated numerically by the Gaussian quadrature
formula since both integrands are regular.

The time-dependent values of the transformed variables can be obtained
by an inverse transform. There are many inversion methods available for the
Laplace transform.As theLaplace transform inversion is an ill-posedproblem,
small truncation errors can be greatly magnified in the inversion process and
lead to poor numerical results. In the present analysis the Stehfest’s algorithm
(Stehfest, 1970) is used. An approximate value fa of the inverse f(t) at a
specific time t is given by

fa(t)=
ln2

t

N∑

i=1

vif
(ln2
t
i
)

(4.26)



Evaluation of fracture parameters for... 625

where

vi =(−1)i+N/2
min(i,N/2)∑

k=(i+1)/2

kN/2(2k)!

(N/2−k)!k!(k −1)!(i−k)!(2k − i)!

The number N = 10 with single precision arithmetic has been found to be
optimal for obtaining accurate results. It means that at each time t, N bo-
undary value problems for the corresponding Laplace transform parameters
p = (i ln2)/t, with i = 1,2, . . . ,N, need to be solved. If M is the number of
the time instants at which the quantity f(t) has to be found, the number of
the Laplace transform solutions f(pj) is then M ×N.

5. Numerical results

5.1. A finite plate with a central crack

In the first numerical example, a rectangular orthotropic plate with a cen-
tral crack is analyzed.Theplate considered is subjected to a thermal loadwith
different prescribed temperatures at its bottom and top sides, which are con-
strained in x2-direction such that it is equivalent to the casewith a prescribed
strain ε0 = α(x1)θ(h), as shown in Fig.3.

Fig. 3. A finite plate with a center crack parallel to the material gradation

The following geometry is considered: w =10m, a/w =0.1, 0.2, 0.3 and
h = w. First, isotropic material properties with an exponential variation of



626 J. Sladek et al.

the Young’s modulus parallel to the crack line are considered

α(x)= α0 exp(δx1) E(x)= E0 exp(γx1) (5.1)

with α0 = 1.67 · 10−5deg−1, E0 = 105MPa, and a constant Poisson’s ratio
ν = 0.3. In this case, the temperature distribution is a function of x2 since
the thermal conductivity is considered to be uniform. By virtue of symmetry,
only one half of the cracked plate is numerically analyzed. A regular node
distributionwith 61×30=1830 nodes (61 nodes along each line x2 = const)
is used in our numerical calculations. The integration path Γ0 for evaluation
of the M-integrals has a rectangular shape, as shown in Fig.3 (dashed line,
5× 5m). To test the accuracy of the proposed method, an additional inte-
gration path with a size 20% larger than the dashed line (i.e. 6× 6m) has
also been used for evaluating the fracture parameters. The discrepancies of
the numerical values obtained for these fracture parameters from the two pa-
ths were less than 1%. The stress intensity factors obtained are normalized as
fI(±a) = KI(±a)/[E0α(a)θ

√
πa]. The numerical results of the variations of

the normalized stress intensity factor with the crack length at both crack tips
in an isotropic FGMwith γ =0.25m−1 are presented inFig.4 for δ =0.1m−1

and δ = −0.1m−1. The normalized stress intensity factor at the right crack
tip increases with increasing crack length; the opposite trend is observed at
the left crack tip, however. For an increasing gradation of the thermal expan-
sionwith the x1-coordinate, the normalized stress intensity factor at the right
crack tip is also enhanced with respect to the homogeneous case. For a de-
creasing gradation of the thermal expansion, the normalized stress intensity
factor at x1 = a is reduced, while at x1 =−a, it is enhanced.

Fig. 4. Variation of stress intensity factors with crack length at both crack tips in an
isotropic FGMwith γ =0.25m−1 and δ =0.1m−1 (a), and δ =−0.1m−1 (b)



Evaluation of fracture parameters for... 627

Next, orthotropic material properties with a constant Poisson’s ratio
ν12 = 0.3 and a constant shear modulus G12 = 4 · 104MPa are considered
for the crack problem analyzed above. The Young’s moduli are expressed as
functions of the parameter R = E1/E2 with E1 = G12(R +2ν12 +1) and
E2 = E1/R. Two different ratios R =0.5 and R =4.5 are considered in our
numerical analyses, and Ei have an exponential variation in the xi-direction
as follows

Ei(x)= Ei0exp(γx1) (5.2)

Numerical results for the normalized stress intensity factors fI(±a) =
KI(±a)/[α(a)θE20

√
πa] are given in Fig.5. The gradation exponent γ is the

same as in the previous isotropic case. One can observe that the orthotropy
parameter R has a relatively small influence on the value of the normalized
stress intensity factors at both crack tips, at least for the case considered here.

Fig. 5. Variation of stress intensity factors with crack length at both crack tips in an
anisotropic FGMwith δ =0.1m−1 and γ =0.25m−1

The T-stresses arealso analyzed for the sameboundaryvalueproblemas in
the previous examples. First, an isotropic FGM is considered.Variations of the
normalized T-stresses T

√
πa/KI with the crack length at x1 = a and x1 =

−a are presented inFig.6a andFig.6b, respectively. The gradation of thermal
expansion coefficient has a small influence on T-stress values.With increasing
crack length, the absolute value of the T-stresses is reduced, andat both crack
tips the T-stresses are almost the same. Amuchmore significant influence of
orthotropic properties on the T-stresses is observed in Fig.7. Increasing the
ratio R increases the T-stresses. It is in agreementwith the analytical solution
for a crack in an infinite plane (Gao andChiu, 1992) and the FEM results for
mechanical loading (Kim and Paulino, 2004).



628 J. Sladek et al.

Fig. 6. Variation of T-stress with crack length at the right crack tip x1 = a (a) and
x1 =−a (b) in an isotropic FGMwith γ =0.25m−1

Fig. 7. Variation of T-stress with crack length at both crack tips in an anisotropic
FGMwith δ =0.1m−1 and γ =0.25m−1

5.2. A finite plate with an edge crack

In the next example, a rectangular orthotropic and linear elastic FGM
plate with an edge crack subjected to an impact of mechanical loading is
analyzed. The plate length is 2h = 30cm, the width w = 10cm, and the
crack length a/w =0.4 (see Fig.8). At the top and the bottom of the plate, a
uniform impact tensile stress σ22(t)= σH(t−0) with the Heaviside step time
variation is applied.Orthotropicmaterial propertieswith a constant Poisson’s
ratio ν12 = 0.3, a constant shear modulus G12 = 0.785 · 1010N/m2, and a
mass density ρ =5 ·103kg/m3 are considered. The Young’s moduli Ei have



Evaluation of fracture parameters for... 629

the same exponential variation (5.2) as in the previous static case. A regular
node distribution with 930 nodes is used for theMLS approximation.

Fig. 8. An edge-cracked orthotropic plate with thematerial gradation in x1-direction

The dynamic stress intensity factor is normalized by the static value
KstatI = σ

√
πa for convenience. The time variations of the mode I dynamic

stress intensity factor are presented in Fig.9. The gradient parameter γ cor-
responds to the ratio of Young’s moduli E1w/E10 = 5.0 in the FGM pla-
te. The influence of the material anisotropy, characterized by the parameter
R = E1/E2, on the dynamic stress intensity factor is presented in Fig.9. The
numerical results for the corresponding isotropic case have been obtained by
Sladek et al. (2005). It can be seen that if the Young’s modulus in the x1-
direction is lower than in the direction perpendicular to the crack line, e.g.
R = 0.5, the wave velocity in the x1-direction is lower too, and the peak
values of the normalized mode I dynamic stress intensity factor are reached
at larger time instants than in the isotropic case R = 1.0. For R > 1.0 the
effect of the material anisotropy on the position of the peak KI(t)-values is
expected to be opposite.

The time variation of the normalized T-stress T/T0 is given in Fig.10,
where the corresponding static equivalent is T0/σ = −0.32 for the isotropic
case. In this case the time variations for the SIF and the T-stresses are very
similar. On the basis of the asymptotic expansion of stresses, one can say that
the time variation of the stress component σ22 ahead of a crack tip follows
that of σ11 at a small distance behind a crack tip for the considered boundary
value problem. For the orthotropic plate with the ratio R =0.5, the T-stress
is lower than that for an isotropic plate. The position of the peak T-value



630 J. Sladek et al.

Fig. 9. Normalizedmode I dynamic stress intensity factors for an FGMplate with
an edge crack

is shifted to larger time instants than in the variation of the stress intensity
factor, since thewave propagation velocity in the x1-direction is reduced.The
influence of thematerial orthotropy on the peak values of the T-stress ismore
pronounced than that for the stress intensity factors.

Fig. 10. Time variation of T-stresses for an FGMplate with an edge crack



Evaluation of fracture parameters for... 631

6. Conclusions

This paper presents efficient numerical methods for the evaluation of stress
intensity factors and T-stresses for crack problems in orthotropic, functio-
nally graded materials under thermal and impact mechanical loads. Path-
independent integral representations for both fracture parameters are derived.
Thepresent integralmethodsare numericallymore expedient than thosebased
on direct computation of the fracture parameters from the asymptotic expan-
sion of the stresses and/or displacements. The integral approach is well suited
for elastic meshless analyses. A local boundary integral equation formulation
based on theMLPGwith ameshless approximation has been successfully im-
plemented to solve 2D boundary and initial-boundary value problems. A unit
step function is used as the test function in the local symmetric weak form on
the local subdomains. The analyzed domain is divided into small overlapping
circular sub-domains on which the local boundary integral equations are ap-
plied. The derived local boundary-domain integral equations are non-singular.
The proposed method is a truly meshless method, wherein no elements or
background cells are involved in either the interpolation or the integration.
For elastodynamic crack problems, the Laplace transform technique is imple-
mented. Themain difficulty in the application of the classical boundary inte-
gral equation formulations for nonhomogeneous, anisotropic and linear elastic
solids is the absence of well-established fundamental solutions. This difficulty
is circumvented by using the present local integral equation method.

Acknowledgements

The authors acknowledge the supports offered by the Slovak Research andDeve-

lopment Support Agency registered under the project number APVT-20-035404, by

the Slovak Grant Agency under the project number VEGA-2/6109/6, by the Ger-

man Research Foundation (DFG), and the project supported jointly by the German

Academic Exchange Service (DAAD) and the Ministry of Education of the Slovak

Republic under the project number D/04/25722.

References

1. Albuquerque E.L., SolleroP., AliabadiM.H., 2002, The boundary ele-
mentmethod applied to time dependent problems in anisotropicmaterials, Int.
Journal Solids and Structures, 39, 1405-1422



632 J. Sladek et al.

2. Albuquerque E.L., Sollero P., Aliabadi M.H., 2004, Dual boundary
elementmethod foranisotropicdynamic fracturemechanics, Int. J.Num.Meth.
Eng., 59 1187-1205

3. Aliabadi M.H., 1997, Boundary element formulations in fracture mechanics,
Mech. Rev., 50, 83-96

4. Aoki S., KishimotoK., SakataM., 1981,Crack tip stress and strain singu-
larity in thermally loaded elastic-plasticmaterials, J. Appl. Mech., 48, 428-429

5. Atluri S.N., Shen S., 2002, The Meshless Local Petrov-Galerkin (MLPG)
Method, Techn. Science Press

6. Atluri S.N., 2004, The Meshless Local Petrov-Galerkin Method for Domain
& BIE Discretizations, Techn. Science Press

7. Belytschko T., LuY., Gu L., 1994, Element freeGalerkinmethods, Int. J.
Num. Meth. Eng., 37, 229-256

8. Betegon C., Hancock J.W., 1991, Two parameter characterization of
elastic-plastic crack tip fields, J. Appl. Mech., 58, 104-110

9. Cotterell B., Rice J.P., 1980, Slightly curved or kinked cracks, Int. J.
Fracture, 16, 155-169

10. Ding H., Liang J., Chen B., 1997, The unit point force solution for both
isotropic and transversaly isotropicmedia,Communications in Numerical Me-
thods in Engineering, 13, 95-102

11. Dolbow J., Gosz M., 2002, On the computation of mixed-mode stress in-
tensity factors in functionally graded materials, Int. J. Solids Structures, 39,
2557-2574

12. ErdoganF., 1995,Fracturemechanics of functionally gradedmaterials,Com-
pos. Eng., 5, 753-770

13. Eshelby J.D, Read W.T., Shockley W., 1953, Anisotropic elasticity with
applications to dislocations,Acta Metallurgica, 1, 251-259

14. Gao H., Chiu C.H., 1992, Slightly curved or kinked cracks in anisotropic
elastic solids, Int. J. Solids Structures, 29, 947-972

15. Gu P., Asaro R.J., 1997, Cracks in functionally graded materials, Int. J.
Solids and Structures, 34, 1-17

16. JinZ.H.,NodaN., 1993a,Minimization of thermal stress intensity factor for a
crack in ametal-ceramicmixture,Trans. American Ceramic Soc. Functionally
Gradient Materials, 34, 47-54

17. Jin Z.H., Noda N., 1993b, An internal crack parallel to the boundary of
a nonhomogeneous half plane under thermal loading, Int. J. Eng. Sci., 31,
793–806.



Evaluation of fracture parameters for... 633

18. Jin Z.H., Noda N., 1994, Crack tip singular fields in nonhomogeneousmate-
rials, J. Appl. Mech., 64, 738-739

19. Kfouri A.P., 1986, Some evaluations of the elastic T-term using Eshelby’s
method, Int. J. Fracture, 30, 301-315

20. Kim J.H., Paulino G.H., 2002, Finite element evaluation of mixed mode
stress intensity factors in functionally graded materials, Int. J. Num. Meth.
Eng., 53, 1903-1935

21. KimJ.H.,PaulinoG.H., 2003a, T-stress,mixedmode stress intensity factors,
and crack initiation angles in functionally gradedmaterials:A unified approach
using the interaction integralmethod,ComputerMeth. Appl. Mech. Eng., 192,
1463-1494

22. Kim J.H., Paulino G.H., 2003b, The interaction integral for fracture of or-
thotropic functionally graded materials: Evaluation of stress intensity factors,
Int. J. Solids and Structures, 40, 3967-4001

23. Kim J.H., Paulino G.H., 2004, T-stress in orthotropic functionally graded
materials: Lekhnitskii and Stroh formalisms, Int. J. Fracture, 126, 345-384

24. KishimotoK.,Aoki S., SakataM., 1980,On thepath independent integral-
Ĵ, J. Eng. Fracture Mech., 13, 841-850

25. KöglM.,GaulL., 2000,A3-Dboundaryelementmethod fordynamicanaly-
sis of anisotropic elastic solids,CMES: Computer Modeling in Eng. and Scien-
ces, 1, 27-43

26. Larsson S.G., Carlsson A.J., 1973, Influence of nonsingular stress terms
and specimen geometry on small-scale yielding at crack tips in elastic-plastic
materials, J. Mech. Phys. Solids, 21, 263-277

27. Lekhnitskii S.G., 1963,Theory of Elasticity of an Anisotropic Body, Holden
Day

28. LeeversP.S.,RadonJ.C., 1982, Inherent stressbiaxiality in various fracture
specimen geometries, Int. J. Fracture, 19, 311-325

29. Mikhailov S.E., 2002, Localized boundary-domain integral formulations for
problemswith variable coefficients,Eng. Analysis with Boundary Elements,26,
681-690

30. Nakamura T., Parks D.M., 1982, Determination of elastic T-stresses along
three-dimensional crack frontsusingan interaction integral, Int. J. Solids Struc-
tures, 29, 1597-1611

31. Nemat-Alla M., NodaN., 1996, Thermal stress intensity factor for functio-
nally gradient half space with an edge crack under thermal load,Archive Appl.
Mech., 66, 569-580

32. Noda N., Jin Z.H., 1993a, Thermal stress intensity factors for a crack in a
functionally gradientmaterial, Int. J. Solids Structures, 30, 1039-1056



634 J. Sladek et al.

33. Noda N., Jin Z.H., 1993b, Steady thermal stresses in an infinite nonhomoge-
neous elastic solid containing a crack, J. Thermal Stresses, 16, 181-196

34. Noda N., Jin Z.H., 1995, Crack tip singularity fields in non-homogeneous
body under thermal stress fields, JSME Int. Jour. Ser. A, 38, 364-369

35. Nowacki W., 1975,Dynamic Problems of Thermoelasticity, PWN,Warsaw

36. Olsen P.C., 1994, Determining the stress intensity factors KI, KII and the
T-term via the conservation laws using the boundary element method, Eng.
Fracture Mech., 49, 49-60

37. Ozturk M., Erdogan F., 1997,Mode I crack problem in an inhomogeneous
orthotropic medium, Int. J. Eng. Sci., 35, 869-883

38. Ozturk M., Erdogan F., 1999, The mixedmode crack problem in an inho-
mogeneous orthotropic medium, Int. J. Fracture, 98, 243-261

39. Paulino G.H., Jin Z.H., Dodds R.H., 2003, Failure of functionally graded
materials, In:Comprehensive Structural Integrity, Vol.2, KarihalooB., Knauss
W.G. (edit.), Elsevier Science, 607-644

40. PaulinoG.H.,Kim J.H., 2004,A new approach to compute T-stress in func-
tionally gradedmaterials bymeans the interaction integralmethod,Eng. Frac-
ture Mechanics, 71, 1907-1950

41. Rao B.N., Rahman S., 2003, Mesh-free analysis of cracks in isotropic func-
tionally gradedmaterials,Eng. Fracture Mechanics, 70, 1-27

42. Schclar S.N., 1994,Anisotropic Analysis Using Boundary Elements, Compu-
tationalMechanics Publications

43. ShahP.D.,TanC.L.,WangX., 2005,T-stress solutions for two-dimensional
crack problems in anisotropic elasticity using the boundary element method,
Fatigue Fract. Eng. Mater. Struc., accepted for publication

44. Sherry A.H., France C.C., Goldthorpe M.R., 1995, Compendium of
T-stress solutions for two and three dimensional cracked geometries, Fatigue
Fract. Eng. Mater. Struct., 18, 141-155

45. SihG.C.,ParisP.C., IrwinG.R., 1965,Oncracks in rectilinearlyanisotropic
bodies, Int. J. Fracture Mechanics, 1, 189-203

46. Sladek V., Sladek J., Markechova I., 1993, An advanced boundary ele-
ment method for elasticity problems in nonhomogeneous media, Acta Mecha-
nica, 97, 71-90

47. Sladek J., SladekV., 1997a,Evaluations of the T-stress for interface cracks
by the boundary element method, Eng. Fracture Mech., 56, 813-825

48. Sladek J., Sladek V., 1997b, Evaluation of T-stresses and stress intensity
factors in stationary thermoelasticity by the conservation integralmethod, Int.
J. Fracture, 86, 199-219



Evaluation of fracture parameters for... 635

49. Sladek J., SladekV., Fedelinski P., 1997, Integral formulation for elasto-
dynamic T-stresses, Int. J. Fracture, 84, 103-116

50. Sladek J., Sladek V., Atluri S.N., 2000, Local boundary integral equ-
ation (LBIE) method for solving problems of elasticity with nonhomogeneous
material properties,Computational Mechanics, 24, 456-462

51. Sladek J., Sladek V., Van Keer R., 2003a,Meshless local boundary inte-
gral equationmethod for 2Delastodynamicproblems, Int. J. Num.Meth. Eng.,
57, 235-249

52. Sladek J., Sladek V., Zhang Ch., 2003b, Application of meshless local
Petrov-Galerkin (MLPG) method to elastodynamic problems in continuously
nonhomogeneous solids, CMES: Computer Modeling in Eng. & Sciences, 4,
637-648

53. Sladek J., Sladek V., Atluri S.N., 2004, Meshless local Petrov-Galerkin
method in anisotropic elasticity,CMES: Computer Modeling in Eng. & Scien-
ces, 6, 477-489

54. SladekJ., SladekV., ZhangCh., 2005,Anadvancednumericalmethod for
computing elastodynamic fracture parameters in functionally gradedmaterials,
Computational Materials Science, 32, 532-543

55. Smith D.J., Ayatollahi M.R., Pavier M.J., 2001, The role of T-stress in
brittle fracture for linear elastic materials under mixed-mode loading, Fatigue
Fracture Eng. Material Structure, 24, 137-150

56. Stehfest H., 1970, Algorithm 368: numerical inversion of Laplace transform,
Comm. Assoc. Comput. Mech., 13, 47-49

57. Sumpter J.D.G., 1993, An experimental investigation of the T-stress appro-
ach, In:Constraint Effects in Fracture, ASTM STP, 1171, 429-508

58. Sumpter J.D.G., Hancock J.W., 1991, Shallow crack toughness of HY80
welds – an analysis based on T, Int. J. Press. Vess. Piping, 45, 207-229

59. Suresh S., Mortensen A., 1998,Fundamentals of Functionally Graded Ma-
terials, Institute of Materials, London

60. YueZ.Q.,XiaoH.T.,ThamL.G., 2003,Boundary element analysis of crack
problems in functionally graded materials, Int. J. Solids and Structures, 40,
3273-3291

61. WangC.Y., Achenbach J.D., 1996,Twodimensional time domainBEMfor
scattering of elastic waves in solids of general anisotropy, Int. Journal of Solids
and Structures, 33, 3843-3864

62. Williams M.L., 1957, On the stress distribution at the base of stationary
crack, J. Appl. Mech., 24, 109-114

63. Williams J.G., Ewing P.D., 1972, Fracture under complex stress – the an-
gled crack problem, J. Fracture Mech., 8, 441-446



636 J. Sladek et al.

Wyznaczanie parametrów pękania dla szczelin w materiałach

funkcjonalnie gradientowych metodą bezsiatkową

Streszczenie

Przedstawiono bezsiatkowąmetodę analizy szczelin opartą na podejściu Petrova-
Galerkina dla dwuwymiarowych liniowo-sprężystych i anizotropowych ośrodków
o zmieniających się własnościachmateriałowych.Rozważono zarównokwazistatyczne
problemynaprężeń cieplnych, jak i zagadnienia elastodynamiki, w których zastosowa-
no aparat transformacji Laplace’a.Badanyobszarpodzielono namałe podobszaryko-
łowe. Jako funkcję testowąw lokalnej, słabej postaci zastosowano jednostkową funkcję
schodkową, co prowadzi do lokalnych równań całkowych (LIE). Metodę ruchomych
najmniejszych kwadratów (MLS) zastosowano do przybliżenia wielkości fizycznych
w LIE. Przedstawiono efektywne metody numeryczne wyznaczania parametrów pę-
kania, a w szczególności współczynników koncentracji naprężeń oraz naprężeń T dla
szczelinwmateriałach funkcjonalnie gradientowych(FGM).Przedstawiononiezależne
od drogi całkowania reprezentacje tych parametrów w materiałach FGM o kontynu-
alnie zmieniającej się niejednorodności.

Manuscript received January 4, 2006; accepted for print April 4, 2006