Full text of DOI:10.1016/S0022-5096(01)00105-3

Journal of the Mechanics and Physics of Solids
50 (2002) 695–716
www.elsevier.com/locate/jmps
Self-consistent elastoplastic stress solutions
for functionally graded material systems
subjected to thermal transients
S. Pitakthapanaphong, E.P. Busso∗
Department of Mechanical Engineering, Imperial College of Science, Technology and Medicine,
Exhibition Road, London SW7 2BX, UK
Received 11 April 2001; received in revised form 23 July 2001
Abstract
In this work, a self-consistent constitutive framework is proposed to describe the behaviour
of a generic three-layered system containing a functionally graded material (FGM) layer sub-
jected to thermal loading. Analytical and semi-analytical solutions are obtained to describe the
thermo-elastic and thermo-elastoplastic behaviour of a three-layered system consisting of a metal-
lic and a ceramic layer joined together by an FGM layer of arbitrary composition proÿle. So-
lutions for the stress distributions in a generic FGM system subjected to arbitrary temperature
transient conditions are presented. The homogenisation of the local elastoplastic FGM behaviour
in terms of the properties of its individual phases is performed using a self-consistent approach.
In this work, power-law strain hardening behaviour is assumed for the FGM metallic phase.
The stress distributions within the FGM systems are compared with accurate numerical solutions
obtained from ÿnite element analyses and good agreement is found throughout. Solutions are
also given for the critical temperature transients required for the onset of plastic deformation
within the three-layered systems. ? 2002 Elsevier Science Ltd. All rights reserved.
Keywords: Functionally graded materials; Self-consistent; Homogenisation; Constitutive formulations
1. Introduction
In a large number of engineering applications, individual components are exposed
to local loading conditions which can vary greatly with location. In such cases, the
∗ Corresponding author. Tel.: +44-0171-594-7084; fax: +44-0171-823-8845.
E-mail address: e.busso@ic.ac.uk (E.P. Busso).
0022-5096/02/$ - see front matter ? 2002 Elsevier Science Ltd. All rights reserved.
PII: S0022-5096(01)00105-3

696
S. Pitakthapanaphong, E.P. Busso / J. Mech. Phys. Solids 50 (2002) 695–716
use of dissimilar materials, such as metals and ceramics, may be required. However,
the abrupt transitions in microstructure and composition between the di@erent materials
often result in high residual stress regions and local stress concentrations which can
lead to the subsequent nucleation of microcracks at or near the bimaterial interface.
The intensity and extent of the stress concentration e@ects due to the large mismatch in
properties can be substantially reduced if the microstructure is gradually changed from
that of the metal to that of the ceramic. A relatively new class of materials known as
functionally graded materials (FGMs) has recently been developed to minimise such
property mismatch e@ects.
The composition and microstructure of an FGM varies with location, resulting in
spatially dependent properties. The composition proÿle can take di@erent forms de-
pending on the required performance. Within an FGM, the di@erent material phases
have di@erent functions. For example, in a metal–ceramic FGM, the metal-rich side
is typically placed in regions where mechanical properties, such as toughness, need to
be high. In contrast, the ceramic-rich side, which has lower thermal conductivity and
can withstand high temperatures, is placed in regions where there are potentially sharp
temperature gradients.
The real advantages of using an FGM as an alternative to two dissimilar materials
joined directly together include: smoothing of thermal stress distributions across the
layers (Choules and Kokini, 1996; Wetherhold et al., 1996), minimisation or elimination
of stress concentrations and singularities at free edges (Lee and Erdogan, 1994=1995;
Yang and Munz, 1997), increase in the bimaterial bonding strength (Howard et al.,
1994a, b), and improved fracture toughness compared to that of monolithic ceramics as
a result of the plastic deformation of the metallic phase (Erdogan, 1995; Jin and Batra,
1996; Shaw, 1998). Typical applications of FGMs are in high temperature coatings
for combustion chambers and airfoils in the aerospace and power generation industries,
coatings in microelectronics and optoelectronics applications, orthopaedic implants and
wear-resistant coatings in bearings, gears, cams and machine tools (Miyamoto, 1996).
Analytical and numerical studies have been carried out to investigate the thermo-
mechanical behaviour of FGMs (e.g. Freund, 1993; Giannakopoulos et al., 1995; Finot
and Suresh, 1996; Dao et al., 1997; Weissenbek et al., 1997). The results of a micro-
scopic study, where the anisotropic properties of each individual metallic grain in the
FGM are described using crystal plasticity, have been reported in Dao et al. (1997).
Only limited amount of analytical work describing the thermo-elastic behaviour of
FGMs have been reported to date and most of it has been limited to FGMs with linear
compositional gradations. This includes the work of Freund (1993), where solutions
for the stress distributions in graded semiconductor layers were derived using linear
elastic plate theory. The resulting expressions for the in-plane strains and curvatures
were later applied to a Ni–Al2O3 system by Giannakopoulos et al. (1995). Since the
thermo-mechanical behaviour of an FGM strongly depends on its composition proÿle,
which in many instances is far from linear due to feasibility in manufacturing, the
ability to understand and predict the behaviour of an FGM system with a generic com-
position proÿle is of prime importance. Accurate FGM constitutive models can, in turn,
be used as a design tool for FGM manufacturing so that the optimal proÿle is tailored
to suit a speciÿc application.

S. Pitakthapanaphong, E.P. Busso / J. Mech. Phys. Solids 50 (2002) 695–716
697
In a metal–ceramic FGM, yielding in the metallic phase can greatly a@ect the stress
distribution within the FGM and thus a stress analysis based only on elastic material
behaviour will not be accurate. A similar issue applies to metal matrix composites
(e.g. see Olsson et al., 1995). In such cases, it is worth noting that plastic deforma-
tion can occur solely as a result of thermal loads. Numerical investigations have also
been carried out to study the thermo-elastoplastic behaviour of metal–ceramic FGMs
using ÿnite element (FE) techniques (Giannakopoulos et al., 1995; Finot and Suresh,
1996; Weissenbek et al., 1997). However, no analytical or semi-analytical solutions
have been reported till date. Finite element analyses, although accurate, involve exten-
sive and costly computations. Therefore there is a strong need for accurate analytical
formulations to predict nonlinear FGM behaviour.
The homogenisation of the local FGM material behaviour can be performed using
a self-consistent scheme. In the classical self-consistent approach, (e.g. see Kroner,
1961; Budiansky and Wu, 1962), the interaction of a single grain, or inclusion, with the
surrounding polycrystalline material is approximated by a spherical inclusion embedded
in a homogeneous polycrystalline matrix with the same elastic property as the inclusion.
For an elasto-plastic material, the term ‘self-consistent’ implies that the microscopic
stresses and strains satisfy equilibrium and compatibility at the grain=phase boundaries,
and that the average stress and strain over all the grains must yield the corresponding
macroscopic values. Eshelby’s solutions (Eshelby, 1957) are relied upon to link the
macroscopic stresses (or stress rates, ꢀ˙∗ij) with the microscopic or average stresses (or
stress rates, ꢀ˙ij) in each grain (Kroner, 1961; Budiansky and Wu, 1962),
ꢀ˙ij = ꢀ˙∗ + 2ꢁ∗(1 − ÿ) (ꢂ˙∗p − ꢂ˙p ); for i; j = 1; 2; 3:
(1)
ij
ij
ij
Here, ꢁ∗ is the polycrystal aggregate shear modulus, ÿ Eshelby’s elastic accommo-
dation factor (Eshelby, 1957), and ꢂ˙p and ꢂ˙∗p are the microscopic and macroscopic
ij
ij
plastic strain rates, respectively.
In Eq. (1), an elastic interaction between the grain and the polycrystal aggregate
is implicitly assumed, thus a high constraint is imposed on the inclusion by the sur-
rounding elastic aggregate. In reality, such high constraint is partially relaxed by the
plastic deformation of the polycrystalline aggregate. The work of Berveiller and
Zaoui (1979) addressed this problem by introducing an additional plastic accommoda-
tion factor into Eq. (1). In more recent work, Cailletaud and Pilvin (1995) and Busso
(1999) (see also Busso et al., 2001a, b) described the behaviour of multiphase materi-
als using tensorial interphase accommodation variables, which enable accurate uniaxial
self-consistent predictions to be obtained for inclusion volume fractions above which
self-consistent predictions cease to be accurate (e:g: ¿ 30%).
An alternative homogenisation method to the classical self-consistent approach is
that derived from a variational procedure. The early work of Hashin and Shtrikman
(1963) introduced a variational procedure to estimate the e@ective modulus of elastic
composites with statistically isotropic microstructures. This approach allowed rigorous
upper and lower bounds to be obtained for the e@ective modulus tensor of the com-
posite. A generalisation of Hashin–Shtrikman’s bounds for nonlinear composites was
presented by Talbot and Willis (1985), and an alternative variational formulation was
later on proposed by Ponte Castan˜eda (1991). The latter relies on the e@ective modulus

698
S. Pitakthapanaphong, E.P. Busso / J. Mech. Phys. Solids 50 (2002) 695–716
Fig. 1. Schematic representation of (a) a three-layered model, and (b) coordinate axes and dimensions
of the three-layered system.
tensor of ‘linear elastic comparison composites’, whereby the e@ective stress potentials
of nonlinear composites are expressed in terms of the corresponding potentials for lin-
ear composites with similar microstructural distributions, to generate the corresponding
bounds. Ponte Casten˜eda’s variational principles (1991) have also been used to derive
elastoplastic relations for metal matrix composites reinforced by spheroidal elastic par-
ticles (Li and Ponte Castan˜eda, 1994). The latter formulation will be applied to an
FGM system to provide a basis with which predictions obtained with the homogenisa-
tion formulation to be proposed here can be compared to.
In this work, the analytical thermo-elastic framework proposed by Freund (1993) is
ÿrst generalised to include various forms of nonlinear FGM compositional gradations.
This is then followed by the description of the elastoplastic behaviour of a generic FGM
using a self-consistent constitutive approach (Cailletaud and Pilvin, 1995; Busso, 1999).
The metallic phase within the FGM will be considered to be a generic power-law strain
hardening material. The calibration of the self-consistent formulation using embedded
unit cell FE techniques is next described, and the results compared with the variational
solutions of Li and Ponte Castan˜eda (1994). The self-consistent solutions for the stress
distributions within the three-layered FGM system are then compared with accurate FE
solutions. Finally, the critical temperature transients required for the onset of plastic
deformation within the three-layered systems are given, together with a comparison
between the magnitudes of elastic and elastoplastic curvatures.
2. Thermo-elastic analysis under an arbitrary temperature transient
2.1. Analytical formulation for a generic FGM composition proÿle
Consider a three-layered system made up of a compositionally graded layer in be-
tween two homogeneous materials, as shown in Fig. 1, subjected to a uniform temper-
ature which may vary with respect to time.
Let the total length, width and height of the multilayer system be given by l, w
and (h1 + h2), respectively, see Fig. 1(a). When lꢀ(h1 + h2)ꢀw, the geometry of the
model allows it to be idealised as a plane stress beam, and the variables of interest
depend only on the out-of-plane coordinate, z. Furthermore, a multi-layer system with

S. Pitakthapanaphong, E.P. Busso / J. Mech. Phys. Solids 50 (2002) 695–716
699
l = wꢀ(h1 + h2), corresponds to a system under equal in-plane or biaxial stresses,
herefrom to be referred to as the biaxial stress model. The in-plane geometry of the
layered structure is shown in Fig. 1(b). The FGM layer extends from z =−a to z =+a
and, for continuous property assumptions to be valid, the thickness of this layer must
be signiÿcantly larger than its dominant microstructural length scale (e.g. grain size).
The interfaces between the di@erent layers are assumed to be perfectly bonded at all
times and the multilayer system behaviour to be linear elastic. Moreover, the material
is assumed to be initially stress free.
By assuming small strain ∗kinematics, the total strain, ꢂ˜∗, applied in the plane of the
multilayer system, (where ꢂ˜ = ꢂx∗x for the plane stress case and ꢂ˜∗ = ꢂx∗x = ꢂy∗y for the
equal biaxial stress case), can be decomposed into an elastic, ꢂ˜∗e, and a thermal, ꢂ˜∗th
,
component as
ꢂ˜∗ = ꢂ˜∗e + ꢂ˜∗th
:
(2)
When the total strain ꢂ˜∗ is considered to be a function of the out-of-plane coordinate,
z, it can be shown that the small-strain compatibility equations lead to a linear relation
between the total strain and curvature, Ä,
ꢂ˜∗(z) = ꢂ0 + Äz;
(3)
where ꢂ0 is the strain at the mid-plane of the FGM layer (z = 0).
Under plane stress conditions, the only non-zero stress component, ꢀ˜∗(z), is given
by
ꢀ˜∗(z) = E∗(z) [ꢂ˜∗(z) − ꢃ∗(z)QT]:
(4)
Here, QT = T − T0, where T0 is the stress-free temperature and T the homogeneous
temperature in the multilayer system. For equal biaxial stress conditions,
ꢀx∗x = ꢀy∗y = ꢀ˜∗(z):
(5)
Substitution of Eq. (3) into Eq. (4) gives
ꢀ˜∗(z) = E∗(z) [ꢂ0 + Äz − ꢃ∗(z)QT]:
(6)
Expressions for ꢂ0 and Ä (given in Freund, 1993), can be derived on the basis that the
resultant force and moment arising from the stress distribution through the thickness,
i.e. along the z-axis, must equilibrate any applied in-plane force and bending moment,
respectively.
Let the volume fraction of the layer 2 material within the FGM vary as a function of
ˆ
the coordinate, z, and be arbitrarily deÿned by a generic function, V(z), which satisÿes
the following conditions at the homogeneous layers’ interfaces,
ꢀ
0
1
at z = − a;
at z = a:
ˆ
V(z) =
(7)
The elastic properties of the FGM, given by its Young’s modulus, E∗, and Poisson’s
ratio, ꢄ∗, together with the coeTcient of thermal expansion, ꢃ∗, are assumed to vary

700
S. Pitakthapanaphong, E.P. Busso / J. Mech. Phys. Solids 50 (2002) 695–716
Table 1
Thermo-elastic properties for the metallic (Ni) and ceramic (Al2O3)
phases (Giannakopoulos et al., 1995)
Material
E (GPa)
ꢄ
ꢃ (◦C−1
)
Al2O3
Ni
380
214
0.25
0.31
7:4 × 10−6
15:4 × 10−6
according to the rule of mixture. Then,
∗
ˆ
E (z) = E1 + (E2 − E1)V(z);
∗
ˆ
ꢄ (z) = ꢄ1 + (ꢄ2 − ꢄ1)V(z);
∗
ˆ
ꢃ (z) = ꢃ1 + (ꢃ2 − ꢃ1)V(z);
(8)
where the subscripts 1 and 2 refer to the respective homogeneous material layers (see
Fig. 1).
In the majority of practical applications involving functionally graded materials,
they are used as interlayers between two homogeneous materials. For such cases, the
through-thickness variation of the elastic modulus in the three-layered system is

E
for − h1 6 z 6 − a;

1


E∗(z) =
(9)
ˆ
E1 + QEV(z) for − a 6 z 6 a;
E2



for
a 6 z 6 h2;
where QE = E2 − E1. Similar relations can be obtained for ꢄ∗(z) and ꢃ∗(z).
Substituting Eq. (9) into Eq. (4) leaves the thermo-elastic plane stress solution for
ˆ
an arbitrary FGM characterised by a generic composition proÿle function V(z),

E [ꢂ + Äz − ꢃ QT]
for − h1 6 z 6 − a;

1
0
1


ꢀ˜∗(z) =
ˆ
ˆ
(E1 + QEV(z)) [ꢂ0 + Äz − (ꢃ1 + QꢃV(z))QT] for − a 6 z 6 a;
E2[ꢂ0 + Äz − ꢃ2QT]



for
a6z6h2:
(10)
The thermo-elastic solution for the biaxial stress case (i.e., l = w in Fig. 1(a)) can
be readily obtained from the plane stress solution by replacing Ei by Ei=(1 − ꢄi) in
Eq. (10), for i = 1; 2. Furthermore, for plane strain conditions, Ei and ꢃi in Eq. (10)
should be replaced by Ei=(1 − ꢄi2) and ꢃi(1 + ꢄi), respectively, for i = 1; 2.
2.2. Application to a Ni–FGM–Al2O3 system
Consider a system made up of Ni–FGM–Al2O3 layers, with all layers assumed to
be isotropic elastic materials, free of damage and having the temperature independent
properties given in Table 1.

S. Pitakthapanaphong, E.P. Busso / J. Mech. Phys. Solids 50 (2002) 695–716
701
Fig. 2. Di@erent composition proÿles used to describe the variation of ceramic volume fraction in the FGM
layer.
A thermal history consisting of a temperature drop from a typical processing (e.g.
stress free) temperature of 320◦C to 20◦C is considered. For illustration purposes, the
metal and ceramic layer thicknesses will be taken to be one-half that of the FGM layer,
viz. h = h1 = h2 = a. Although plane stress conditions are assumed, solutions to other
plane conditions can simply be obtained in the way described in Section 2.1.
The compositional gradation of the FGM layer is deÿned by the volume fraction of
ˆ
the ceramic phase. Here, the following functions of V(z) will be considered:
z + a
ˆ
1: Linear:
V(z) =
;
2a
ꢁ
ꢂ
2
z + a
2a
ˆ
2: Quadratic:
V(z) =
;
ꢁ
ꢂ
(11)
2
a − z
ˆ
3: Inverse quadratic: V(z) = 1 −
;
2a
ꢁ
ꢂ
2 − 2
ꢁ
ꢂ
3
z + a
2a
z + a
2a
ˆ
V(z) = 3
4: Cubic:
:
These functions are illustrated in Fig. 2.
To verify the accuracy of the analytical solutions, ÿnite element analyses of the
three-layered system were performed using the commercial ÿnite element code ABAQUS
(1999). Due to symmetry considerations, only half of the three-layered system was
modelled using 1000 equal sized linear plane stress elements. The length-to-height
ratio for the model was chosen to be l=(h1 + h2) = 2:5 so as to minimise edge e@ects.

702
S. Pitakthapanaphong, E.P. Busso / J. Mech. Phys. Solids 50 (2002) 695–716
Fig. 3. Analytical thermo-elastic residual stress distributions at room temperature and corresponding FE
predictions using an FGM with a (a) linear, (b) quadratic, (c) inverse quadratic and (d) cubic compositional
variation (T0 = 320◦C).
The FGM layer was modelled with continuously varying composition and properties
through the thickness. To implement that, an implicit user-deÿned material subroutine
was developed and used in the FE calculations.
Figs. 3(a)–(d) show the stress distributions through the thickness of the three-layered
system for the di@erent composition proÿles at room temperature. The analytical pre-
dictions are in good agreement with the FE results. It can be seen that the stress
distributions in the homogeneous metallic and ceramic layers are close to being linear
and hence consistent with Eq. (10), whereas the stress distribution in the FGM layer
depends on the type of compositional gradation. It is also worth noting that the stresses
at the ceramic=FGM interface are always compressive whereas those at the metal=FGM
interface are tensile in all cases. The reverse is true if the FGM system temperature
increases from the stress-free temperature (viz. 320◦C). In the linear proÿle case (Fig.
3(a)), the stress within the FGM varies parabolically and is close to zero at the FGM
mid-plane. Peak tensile and compressive stresses are observed at the ceramic surface
and the ceramic=FGM interface, respectively. The stresses across the metal=FGM in-
terface for the quadratic gradation (Fig. 3(b)) and across the ceramic=FGM interface

S. Pitakthapanaphong, E.P. Busso / J. Mech. Phys. Solids 50 (2002) 695–716
703
for the inverse quadratic gradation (Fig. 3(c)) are continuous. This is due to the very
gradual transition of the FGM composition near those interfaces.
Compared to the linear proÿle, the quadratic gradation gives a measurably smaller
(tensile) stress at the metal=FGM interface at the expense of a higher (compressive)
stress at the ceramic=FGM interface. A lower tensile stress is found in the ceramic layer
and the peak (tensile) stress develops near the FGM mid-plane. The inverse quadratic
gradation reduces the magnitude of the (tensile) stress at the metal surface by 65%
from the linear proÿle case. As can be seen in Fig. 3(d), by using the cubic gradation,
the stresses across both interfaces vary smoothly and the interfacial stress magnitudes
are comparable to the linear case.
Possible failure mechanisms in the ceramic include cleavage fracture from pre-existing
microcracks caused by high tensile stresses. Thus, the inverse quadratic composition
proÿle appears to be the least suitable as it yields the highest tensile stress in the
ceramic. In contrast, high triaxialities (ꢀm = ꢀ˜=3) in the metal layer or the metal-rich
region of the FGM can lead to ductile failure. Again, the quadratic composition proÿle
exhibits the highest triaxialities. Furthermore, since void growth is mainly a@ected by
the mean triaxial stress, it is therefore more likely to take place in the biaxial stress
case than in the plane stress case. Among all composition proÿles considered here, the
linear and the cubic proÿles appear to be the most benign.
3. Thermo-elastoplastic analysis
3.1. Self-consistent constitutive formulation
In this section, the elastoplastic behaviour of the metallic phase within the FGM in
the three-layered system illustrated in Fig. 1 is considered. Layer 1 is assumed to be an
elastoplastic metal and Layer 2 an elastic ceramic material. Plane stress conditions will
be considered as an example. The metal is assumed to have its constitutive behaviour
described by a power-law hardening law, namely

ꢀ˜1
E1




if ꢀ˜1 6 ꢀy ;
1
ꢂ˜1 =
ꢁ
ꢂ
(12)
n

ꢀ˜1



ꢂy
otherwise:
1
ꢀy
1
Here ꢀy and ꢂy (= ꢀy =E1) are the metal yield stress and strain, respectively, and
1
1
1
n is the strain hardening exponent.
The plastic strain rate in the metal is deÿned through an associated Vow rule as
˙
ꢀ˜1
p
˙
ꢂ˜1 =
;
(13)
ˆ
h(ꢀ˜1)
ˆ
where h(ꢀ˜1) is the plastic hardening modulus of the metalli∗cp phase (Phase 1).
˙
The uniaxial equivalent plastic strain rate in the FGM, ꢂ˜ , is simply given by the
volume average of that in the metallic phase, (ꢂ˜1 ),
p
˙
p
1
∗p
˙
˙
ꢂ˜ (z) = f (z)ꢂ˜
(14)
1

704
S. Pitakthapanaphong, E.P. Busso / J. Mech. Phys. Solids 50 (2002) 695–716
where f denotes the volume fraction of the metallic phase within the FGM. In
1
∗p
˙
Eq. (14), the dependency of ꢂ˜ on the out-of-plane coordinate has been included
to account for the spatial variation of the material behaviour.
In an analogous way as in the thermo-elastic case, the total strain in the three-layered
system, ꢂ˜∗, is considered to be made up of an elastic (ꢂ˜∗e), a plastic (ꢂ˜∗p) and a thermal
component (ꢂ˜∗th), and to depend also on the out-of-plane coordinate (z) (see Fig. 1).
∗
˙
The FGM stress rate at a generic location z; ꢀ˜ , is given by
∗
∗
∗p
˙
ꢀ˜ (z) = E (z) [ꢂ˜ (z) − ꢂ˜ (z) − ꢂ˜∗th(z)];
(15)
∗
˙
˙
˙
∗th
∗
∗
˙
˙
˙
where ꢂ˜ = ꢃ (z)T, with T being the temperature rate and E (z) is the FGM modulus
∗
˙
deÿned in Eq. (9). An expression for ꢂ˜ (z) can be readily found from the di@erential
form of Eq. (3),
∗
˙
ꢂ˜ (z) = ꢂ˙0 + Ä˙z;
(16)
where ꢂ˙0 and Ä˙ are their di@erential forms.
The homogenisation of the local elastoplastic properties within the FGM is done
using a self-consistent approach. Following the recent work of Cailletaud and Pilvin
(1995), Busso (1999) and Busso et al. (2001a, b), the homogenisation of the FGM
elastoplastic behaviour for ceramic volume fractions greater than 20% requires that the
s∗tress in each individual FGM phase, ꢀ˜i, be related to the macroscopic FGM stress,
˜
ꢀ˜ , in terms of interphase accommodation factors for each phase, Si. Then,
ꢀꢃ
ꢅ
ꢆ
m
ꢄ
∗
˙
˜
˙
˜
∗
˙
˙
ꢀ˜i(z) = ꢀ˜ (z) + 2ꢁ (1 − ÿ)
f Sk − Si
;
(17)
k
k=1
where ꢁ∗ is the FGM shear modulus, ÿ Eshelby’s elastic accommodation factor and m
the number of elastoplastic phases (m = 1 in this case). The evolution of the interphase
accommodation variable is given by
p
p
i
˙
˜
˙
˜
i
ˆ
Si = ꢂ˜ − H(f)Si|ꢂ˜ |;
(18)
i
ˆ
where H(f) is a dimensionless homogenisation function which needs to be calibrated
i
p
˙
from detailed ÿnite element calculations of the metal–ceramic system, and | ꢂ˜ | is the
i
p
˙
ˆ
magnitude of ꢂ˜ . The calibration of H using an embedded unit cell technique will be
i
described in Section 3.2.
Eqs. (12)–(18) constitute the complete set of constitutive and compatibility relations
for the three-layered system. In this work, they were numerically integrated using an
Euler forward integration scheme and the results will be shown in Section 3.4. For a
general three-dimensional form of the formulation, see Appendix A.
ˆ
3.2. Calibration of H(f) using an embedded unit cell ÿnite element approach
i
An embedded unit cell approach (e.g. Dong and Schmauder, 1996; Forest and Pilvin,
ˆ
1996) is relied upon to calibrate the function H(f) (Eq. (18)) so that self-consistency
i
in both strains and stresses is satisÿed under uniaxial loading conditions. The FE unit
cell model consists of a spherical inclusion, made up of a ceramic core embedded

S. Pitakthapanaphong, E.P. Busso / J. Mech. Phys. Solids 50 (2002) 695–716
705
Fig. 4. FE model for the embedded unit cell calculations with f = 0:4 consisting of a ceramic–metal inclusion
2
embedded in the HEM matrix. Only one-eighth of the unit cell is shown.
in an outer metallic shell which is, in turn, centred inside a cube of a homogeneous
equivalent material (HEM) (see Fig. 4). Note that embedded cell models of this type
are known to represent more accurately the behaviour of composites with randomly
arranged particles (Dong and Schmauder, 1996).
The behaviour of the HEM is described by the proposed FGM self-consistent consti-
tutive formulation and that of the metal and the ceramic phases within the spherical in-
clusion by their bulk elastoplastic and elastic properties, respectively. For self-consistency
ˆ
to be satisÿed, H(f) is calibrated so that, for a given ceramic volume fraction,
i
f (= 1 − f ), the stress–strain response of the HEM matches that of the spherical
2
1
two-phase inclusion. Due to symmetry considerations, only one-eighth of the model
was considered (see Fig. 4). A user-deÿned material subroutine was developed which
allows the material properties at each integration point of the FE model to be identi-
ÿed from its spatial coordinates and the ceramic phase volume fraction. Such approach
enables the use of a single relatively regular 3D FE mesh (see Fig. 4) to analyse the
response of the embedded unit cell model for a wide range of ceramic volume fractions.
A mesh sensitivity study revealed the 3D FE mesh shown in Fig. 7, which consists of
1000 3D quadratic isoparametric elements, to yield mesh-independent results. Symme-
try boundary conditions were prescribed on the X1 = 0, X2 = 0 and X3 = 0 planes and
periodic boundary conditions on the remaining faces of the unit cell. Uniaxial loading
was simulated by imposing appropriate displacement boundary condition along any one
of the three orthogonal directions.
The thermo-elastic properties used for the individual Ni and Al2O3 phases are those
given in Table 1. The Ni yield stress was taken to be temperature-independent and
equal to 100 MPa, and its strain hardening exponent, n = 5 (Frost and Ashby, 1982).

706
S. Pitakthapanaphong, E.P. Busso / J. Mech. Phys. Solids 50 (2002) 695–716
Fig. 5. Comparison between the uniaxial stress–strain response of the ceramic–metal composite calculated
with the FE unit cell models and the self-consistent predictions.
Eqs. (13)–(15) and (17)–(18) were numerically integrated using an Euler forward
ˆ
scheme. The dimensionless parameter H was calibrated for the range 0 ¡ f ¡ 1 and
2
the optimum values found to be linearly dependent on the ceramic volume fraction, f ,
2
ˆ
H(f ) = h0 + h1f for f 6 0:4;
(19)
2
2
2
ˆ
with h0 = 135 and h1 = − 340. If f ¿ 0:4, H = 0.
2
3.3. Uniaxial response of the Ni–Al2O3 composite: a comparison between the
proposed self-consistent model and an existing variational formulation
3.3.1. Self-consistent method
The self-consistent formulation derived in Section 3.1 is used together with the cal-
ˆ
ibrated homogenisation function H to predict the uniaxial stress–strain response of
Ni–Al2O3 composites with a range of reinforcement volume fractions. A comparison
between the self-consistent predictions of the uniaxial FGM stress–strain response and
accurate FE results obtained with the unit cell model is given in Fig. 5. It can be
seen that good agreement is generally found for the range of ceramic volume frac-
tions considered, except at f = 0:3 where the hardening rate is overpredicted at large
2
strains.
3.3.2. Alternative homogenisation procedure based on a variational principle
The e@ective behaviour of a nonlinear composite can alternatively be described using
a variational approach, such as that proposed by Li and Ponte Castan˜eda (1994). This
method will be brieVy summarised next and its predictive capabilities compared with
those of the proposed self-consistent approach for the Ni–Al2O3 system.

S. Pitakthapanaphong, E.P. Busso / J. Mech. Phys. Solids 50 (2002) 695–716
707
In the variational approach, the constitutive behaviour∗of nonlinear composite mate-
˜
rials relies on an e@ective stress potential function, U(ꢀij), such that,
∗
˜
9U(ꢀij)
9ꢀi∗j
ꢂi∗j =
;
(20)
where ꢂi∗j and ꢀi∗j represent the average composite strain and stress components, re-
spectively. For isotropic composites with a power-law strain hardening metallic matrix
reinforced by elastic particles, the uniaxial stress potential is given by (Li and Ponte
Castan˜eda, 1994),

∗2
ꢀ
11
∗



ˆ
ˆ
d(ꢁ0 = ꢁ1)
if "(ꢀ11) 6 0;
2
∗
˜
U(ꢀij) =
(21)
ꢇ
ꢈ
∗2
ꢀ



11
∗
ˆ
ˆ
ˆ
max
d(ꢁ0) − f V(ꢁ0)
if "(ꢀ11) ¿ 0:
1
2
∗
ˆ
Here, "(ꢀ11) is the yield function, ꢁ1 the elastic shear modulus of the metallic phase
1, and ꢁ0 the shear modulus of the linear elastic comparison composite. The function
V(ꢁ0) is, in turn, given by
ˆ
ꢉ
ꢊ
ꢁ
ꢂ
ꢁ
ꢂ
(n+1)=n−1 − 1 ;
(22)
2
#
n − 1
ꢁ1
ꢁ0
y1
ˆ
V(ꢁ0) =
n + 1 2ꢁ1
where #y and n are the metal yield stress in shear and hardening exponent, respectively.
1
ˆ
The function d(ꢁ0) is deÿned as
ꢉ
ꢊ
ꢇ
ꢈ
ˆ
1
9k1
ꢁˆ + f (Cs − 2Gs)
1
3ꢁ0
k + f (Cs + Gs)
1
1
ˆ
d(ꢁ0) =
1 + f
+
1 + f
;
(23)
2
2
C
C
ˆ
where f is the reinforcement volume fraction and f = 1 − f , k = k2=(k1 − k2) and
2
1
2
ꢁˆ = ꢁ2=(ꢁ0 − ꢁ2). Here ki and ꢁi are the bulk and shear moduli of the phase (i), for
i = 1 and 2.
The parameters C; Cs; Ds; Gs and Hs are given by
f
3
1
ˆ
ˆ
C = kꢁˆ + k(2Ds + Cs − 2Hs − 2Gs)
f
3
+
1 ꢁˆ(Ds + 2Cs + 2Hs + 2Gs) + f2(CsDs − 2GsHs);
1
Cs = 1 − (S2(12)22 + S(1) ); Ds = 1 − S(1)
;
1111
2233
Gs = −S(1)
;
Hs = − S(1)
:
2211
(24)
1122
(1)
Here S are the components of the Eshelby’s tensor (Eshelby, 1957) associated with
ijkl
∗
ˆ
ˆ
phase 1. The yield function, "(ꢀ11), deÿned so that " ¡ 0 for elastic deformation,

708
S. Pitakthapanaphong, E.P. Busso / J. Mech. Phys. Solids 50 (2002) 695–716
Fig. 6. Comparison between the uniaxial stress–strain response calculated with the FE unit cell models and
the variational method predictions.
ˆ
ˆ
" = 0 at yield, and " ¿ 0 for plastic loading, is given by
ꢉ
ꢊ
ꢁ
ꢂ
2
∗2
ˆ
ꢀ
9d(ꢁ0 = ꢁ1)
f
2
#
y1
ꢁ1
1
∗
11
ˆ
"(ꢀ11) = −
+
:
(25)
(26)
2
9ꢁ0
Finally, by using Eqs. (20) and (21), the uniaxial stress–strain relation is
ꢀ
∗
∗
ˆ
ˆ
ꢂ11=d(ꢁ0 = ꢁ1) if "(ꢀ11) 6 0;
ꢀ1∗1
=
∗
∗
ˆ
ˆ
ꢂ11=d(ꢁ0 = ꢁ˜0) if "(ꢀ11) ¿ 0;
where ꢁ˜0 is the optimal value of ꢁ0 obtained from Eq. (21).
In Fig. 6, the uniaxial responses predicted by Eqs. (26) for pure Ni reinforced with
di@erent Al2O3 volume fractions (i.e. 10%, 20%, 30% and 60%) are presented together
with the accurate reference solutions obtained from the FE unit cell analyses. It can
be seen that, for this particular composite, the variational model does not predict the
elastoplastic behaviour of the composites as accurately as the proposed self-consistent
method. This is due to the fact that this type of variational formulation is particularly
suited to composite materials with a large elastic property mismatch between its phases.
For the Ni–Al2O3 composite, the elastic modulus ratio is only 1.8 (see Table 1).
Next, the self-consistent approach will be used to study the thermo-elastoplastic
behaviour of the Ni–FGM–Al2O3 three-layered system.
3.4. Application to a Ni–FGM–Al2O3 system
The stress distributions within the Ni–FGM–Al2O3 system using the homogenised
constitutive formulation outlined in Section 3.1 were determined for the composition
proÿles illustrated in Fig. 2. A temperature drop from the stress-free temperature of

S. Pitakthapanaphong, E.P. Busso / J. Mech. Phys. Solids 50 (2002) 695–716
709
Fig. 7. Comparison between the self-consistent solution and FE results at room temperature for a linear
composition proÿle.
620◦C to room temperature is considered as the thermal transient. The cooling rates
are assumed to be such that time-dependent e@ects are negligible.
For comparison, reference solutions were obtained from ÿnite element analyses of
the three-layer system with the elastoplastic behaviour of the FGM given by the results
from the FE unit cell calculations.
Fig. 7 shows a comparison between the self-consistent predictions and the FE results
for the linear composition proÿle with continuously varying FE properties in the FGM
region. The stress distributions corresponding to an elastic, a perfectly plastic (n = ∞),
and a strain hardening (with n = 5) metallic phase are shown. It can be seen that
yielding takes place at the metal surface and at both the metal-rich and ceramic-rich
ends of the FGM. The magnitude of the peak (tensile) stress in the FGM is reduced
by 65% as a result of plasticity in the FGM metallic phase. A sharp inVection in
the stress distribution is observed at the ceramic=FGM interface due to the purely
elastic behaviour of the ceramic. Note that the locations where plasticity is predicted
to occur in Fig. 7 (e.g. at the metal-free surface, the metal=FGM interface and in
the ceramic-rich region) generally agree with the FE predictions of Weissenbek et al.
(1997), where the ceramic particles within the FGM were explicitly modelled as elastic
isotropic materials. However, they found that plastic deformation also occurs in the
middle of the graded FGM layer as a result of the local constraint exerted by the
ceramic particles.
In Figs. 8–10, both the elastic, and the self-consistent elastoplastic solutions with
n = 5 are shown together with the FE results for the nonlinear composition proÿles. It
can be seen that the self-consistent predictions for all the composition proÿles agree
well with the FE results. In the quadratic proÿle, Fig. 8, compressive yielding starts
from the metal surface and tensile yielding occurs within the FGM. The stresses in both
regions show a considerable reduction when compared with the elastic solution. Small
yielded regions are also seen in the ceramic-rich side of the FGM. The inverse quadratic

710
S. Pitakthapanaphong, E.P. Busso / J. Mech. Phys. Solids 50 (2002) 695–716
Fig. 8. Comparison between the self-consistent thermal stress distributions with the FE results for a quadratic
compositional variation (n = 5).
Fig. 9. Comparison between the self-consistent thermal stress distributions with the FE results for an inverse
quadratic compositional variation (n = 5).
proÿle (Fig. 9) shows yielding occurring mainly around the metal=FGM interface and
within the FGM, with the former region being in tension and the latter in compression.
A similar trend can be observed for the cubic proÿle case (Fig. 10).
As the self-consistent FGM predictions cannot account for the e@ects of the packing
arrangement of the individual ceramic particles within the FGM, it is of interest to
investigate their relative accuracy. To that purpose, the uniaxial FE results of Weis-
senbek et al. (1997) obtained from FE unit cell models with square and hexagonal
packings were compared with those obtained with the proposed self-consistent FGM
formulation. Fig. 11 shows the average stress–strain response of the three-layered sys-
tem obtained by loading it along its in-plane direction (i.e. x-axis in Fig. 1) under

S. Pitakthapanaphong, E.P. Busso / J. Mech. Phys. Solids 50 (2002) 695–716
711
Fig. 10. Comparison between the self-consistent thermal stress distributions with the FE results for a cubic
compositional variation (n = 5).
Fig. 11. Comparison between the average in-plane stress–strain response of the three-layered system predicted
by the proposed self-consistent model and by the FE results of Weissenbek et al. (1997) based on unit cell
computations with square and hexagonal packing of the ceramic particles.
plane stress. The dimensions of the three-layered model were chosen to be identical
to those used in Weissenbek et al. (1997). The predicted average in-plane stress–strain
response in the three-layered FGM system is shown in Fig. 11, together with the re-
sults reported by Weissenbek et al. (1997). It can be seen that the average response
given by the self-consistent model lies amongst the periodic unit cell FE results ob-
tained by Weissenbek et al. (1997) with square and hexagonal packing arrangements
for the ceramic particles. The overall hardening behaviour increases in the follow-
ing sequence: square–diagonal packing → hexagonal packing → self-consistent model →
square–edge packing. It is encouraging that the overall in-plane system response

712
S. Pitakthapanaphong, E.P. Busso / J. Mech. Phys. Solids 50 (2002) 695–716
Fig. 12. Radius of curvature of a plane stress three-layered FGM beam undergoing a temperature drop of
QT=QTmax from the stress-free temperature assuming elastic and elastoplastic behaviour (QTmax = 600◦C).
predicted by the self-consistent FGM formulation is well bounded by the periodic
unit cell FE results.
A variable which can be generally measured experimentally more easily than the
stress distributions is the overall curvature of the FGM system. An analytical solu-
tion for the elastic FGM system curvature induced by a temperature change from
a curvature-free reference temperature can be found in Freund (1993). However, an
elastoplastic curvature needs to be calculated numerically from ÿnite element results.
Fig. 12 shows the predicted radius of curvature of a plane stress three-layered FGM
beam subjected to di@erent temperature changes QT, up to a maximum of QTmax
=
600◦C. Here, the geometry of the FE model was chosen to ensure a uniform curvature
throughout the FGM system thickness, and an FGM with an inverse quadratic proÿle
was considered as it shows the most contrasting result. The elastic and elastoplastic
predictions shown in Fig. 12 reveal that plasticity only causes the curvature to devi-
ate a relatively small amount from the elastic curvature. This occurs at approximately
QT ≈ 0:3QTmax, which agrees with the onset of plasticity in the FGM metallic phase.
It is interesting to note that despite the small di@erences in the overall curvatures for
the elastic and elastoplastic cases, the corresponding rearrangement of the local elastic
stresses due to plastic deformation is considerable. For instance, the maximum stress
reduction as a result of plastic deformation was found to be as large as 200% (at
the metal=FGM interface, see Fig. 9). Therefore, it can be concluded that the curva-
ture alone cannot be used to infer accurately the local elastoplastic stresses within the
three-layered FGM system.
If failure is assumed to be controlled by yielding, it is important to determine the
critical temperature transient from the stress-free temperature, QTc, that causes the
onset of plasticity at either the metal surface or the metal=FGM interface. Fig. 13
shows the critical temperature at which yielding starts at either the metal surface or
the metal=FGM interface as a function of the relative FGM thickness, 2a=(h1 + h2),
for the plane stress case and the linear composition proÿle. These results show that,

S. Pitakthapanaphong, E.P. Busso / J. Mech. Phys. Solids 50 (2002) 695–716
713
Fig. 13. Critical temperature transient for the onset of plasticity at the metal surface and the metal=FGM
interface for the linear FGM compositional variation (n = 5).
in a three-layered system of the linear gradation where the FGM layer is relatively
thin, yielding will occur ÿrst at the metal=FGM interface. The critical relative FGM
thickness beyond which yielding switches from occurring ÿrst at the interface to the
metal surface is approximately 2a=(h1 + h2) = 0:46. This is roughly the same as that
of the cubic composition proÿle. Furthermore, these analyses revealed that, for the
quadratic composition proÿle, yielding also occurs ÿrst at the metal=FGM interface
and the corresponding critical relative FGM thickness is approximately 0.28. For the
case of the inverse quadratic composition proÿle, yielding always occurs ÿrst at the
metal=FGM interface regardless of the relative FGM thickness. Based on this type of
information, the geometry of the system or=and the FGM composition proÿle can be
chosen so that, for a given temperature transient, yielding at undesirable locations can
be avoided.
4. Conclusions
The proposed relations for the thermo-elastic stress distributions within a generic
metal–FGM–ceramic system can predict accurately complex stress distributions induced
by thermal transients. By choosing an appropriate FGM compositional gradation, the
stress distribution within the system can be controlled so that undesirable stresses at
critical locations are minimised or avoided. From the choices of the compositional
gradation considered here, the linear and the cubic proÿles o@er the most benign stress
distributions.
Elastoplastic solutions for the residual stress distribution in a generic metal–FGM–
ceramic system have been proposed and the results agree well with accurate ÿnite ele-
ment solutions. Two methods were pursued in order to homogenise the FGM elastoplastic

714
S. Pitakthapanaphong, E.P. Busso / J. Mech. Phys. Solids 50 (2002) 695–716
properties, namely a self-consistent method and a variational principle-based approach.
The uniaxial response of the Ni–Al2O3 composite obtained from both methods were
compared with ÿnite element solutions. The self-consistent method was found to pre-
dict the composite uniaxial response accurately over a wide range of ceramic volume
fractions, whereas the variational formulation, which is better suited to composites with
larger elastic property mismatch than what was considered in this work, is relatively
less accurate. The stress–strain response of a metal–FGM–ceramic system obtained with
the proposed self-consistent method was also in agreement with periodic unit cell FE
results reported in the literature.
Predictions of the critical temperature transients from the stress-free temperature
that causes yielding at either the metal surface or at the metal=FGM interface have
also been obtained. The analytical elastic and semi-analytical self-consistent elastoplas-
tic solutions presented here provide a simple, yet accurate tool for the prediction of
thermally induced stresses in an FGM layer sandwiched between two homogeneous
materials.
Acknowledgements
Support for this work has been provided by the Government of Thailand. The authors
are also very grateful to Prof. Pedro Ponte Castan˜eda, University of Pennsylvania, USA,
for his comments on the discussion of the variational approach. The ABAQUS was
provided under academic license by HKS Inc., Providence, Rhode Island, USA.
Appendix A. 3-Dgeneralisation of the self-consistent thermo-elastoplastic constitutive
formulation
The FGM uniaxial equivalent stress–strain relations expressed in terms of the out-of-
plane coordinate z, see Eq. (15), can be generalised into the following three-dimensional
form:
ꢀ˙i∗j(z) = 2ꢁ∗(z){ꢂ˙∗ (z) − ꢂ˙∗p(z)} + (ij)∗(z){ꢂ˙k∗k(z) − ꢂ˙∗p(z)}
ij
ij
kk
∗
∗
∗
˙
− (ij{2ꢁ (z) + 3) (z)}ꢃ (z)T; i; j; k = 1; 2; 3;
(A.1)
where ꢂ˙∗ is the total strain rate component due to the temperature-induced local curva-
ij
ture change (e.g., see Eq. (16)), ꢂ˙∗p the overall plastic strain rate component, ꢁ∗; )∗
ij
and ꢃ∗ are the FGM shear modulus, Lame constant, and coeTcient of thermal ex-
˙
pansion, respectively, (ij the Kronecker delta and T the temperature rate. The spatial
variations of ꢁ∗, )∗ and ꢃ∗ in the FGM are given by Eq. (8).
The FGM plastic strain rate is deÿned by
ꢂ˙∗p = f ꢂ˙p ;
(A.2)
1
ij
1ij

S. Pitakthapanaphong, E.P. Busso / J. Mech. Phys. Solids 50 (2002) 695–716
715
where f is the volume fraction of the metallic phase and ꢂ˙p its plastic strain rate
1
1ij
given by
ꢀ1ꢀ
ij
ꢂ˙p1 = ꢀ*˙1ꢁ
:
(A.3)
ij
ꢀe
1
Here ꢀ1ꢀ is the deviatoric stress component of the metallic phase, ꢀe the equivalent
1
ij
stress,
ꢁ
ꢂ
1=2
3
2
ꢀe =
ꢀ1ꢀ ꢀ1ꢀ
(A.4)
(A.5)
1
ij
ij
and ꢀ*˙1ꢁ a switching parameter so that

3 ꢀ˙e
1

if ꢀe ¿ ꢀy while loading;
1
1
2 h
ꢀ*˙1ꢁ =

0
if ꢀe ¡ ꢀy or ꢀe = ꢀy while unloading;
1
1
1
1
ˆ
where h = h(ꢀe ) is the plastic hardening modulus of the metallic phase, and ꢀy the
1
1
metal yield stress.
The stress within a generic elastoplastic FGM phase, r, can be expressed in terms
of the corresponding interphase accommodation factor, Sr , as,
ij
ꢀ
ꢆ
m
ꢄ
ꢀ˙r = ꢀ˙∗ij + 2ꢁ∗(1 − ÿ)
f Sk − Sr
;
(A.6)
˙
˙
k
ij
ij
ij
k=1
where ÿ is Eshelby’s elastic accommodation factor and m the number of elastoplastic
phases (= 1 in the Ni–Al2O3 system considered here). The evolution of Sr is given
ij
by
p
rij
p
˙
ˆ
Sr = ꢂ˙ − H(f )Sr |ꢂ˙ |;
(A.7)
r
ij
ij rij
ˆ
where H(f ) is the dimensionless homogenisation function.
r
References
ABAQUS, 1999. Version 5.8. HKS Inc., Providence, RI, USA.
Berveiller, B., Zaoui, A., 1979. An extension of the self-consistent scheme to plastically Vowing polycrystals.
J. Mech. Phys. Solids 26, 325–344.
Budiansky, B., Wu, T.Y., 1962. Theoretical prediction of plastic strains of polycrystals. Proceedings of the
Fourth U.S. National Congress of Applied Mechanics, p. 1175.
Busso, E.P., 1999. Oxidation-induced stresses in ceramic–metal interfaces. J. Phys. IV 9, 287–296.
Busso, E.P., Lin, J., Sakurai, S., Nakayama, M., 2001a. A mechanistic study of oxidation-induced degradation
in a plasma-sprayed thermal barrier coating system. Part I: model formulation. Acta Mater. 49=9,
1515–1528.
Busso, E.P., Lin, J., Sakurai, S., 2001b. A mechanistic study of oxidation-induced degradation in a
plasma-sprayed thermal barrier coating system. Part II: life prediction model. Acta Mater. 49=9,
1529–1536.
Cailletaud, G., Pilvin, P., 1995. Rev. Eur. EleÃments Finis 3, 515–541.

716
S. Pitakthapanaphong, E.P. Busso / J. Mech. Phys. Solids 50 (2002) 695–716
Choules, B., Kokini, K., 1996. Architecture of functionally graded ceramic coatings against surface thermal
fracture. J. Eng. Mater. Technol. Trans. ASME 118, 522–528.
Dao, M., Gu, P., Maewal, A., Asaro, R., 1997. A micromechanical study of residual stresses in functionally
graded materials. Acta Mater. 45 (8), 3265–3276.
Dong, M., Schmauder, S., 1996. Modeling of metal matrix composites by a self-consistent embedded cell
model. Acta Mater. 44 (6), 2465–2478.
Erdogan, F., 1995. Fracture mechanics of functionally graded materials. MRS Bull. 43–44.
Eshelby, J.D., 1957. The determination of the elastic ÿeld of an ellipsoidal inclusion and related problems.
Proc. Royal Soc. London A241, 398.
Finot, M., Suresh, S., 1996. Small and large deformation of thick and thin ÿlm multi-layers: e@ects of layer
geometry, plasticity and compositional gradients. J. Mech. Phys. Solids 44 (5), 683–721.
Forest, S., Pilvin, P., 1996. Modelling the cyclic behaviour of two-phase single crystal nickel-base superalloys.
In: Pineau, A., Zaoui, A. (Eds.), IUTAM Symposium on Micromechanics of Plasticity and Damage of
Multiphase Materials 1996. Kluwer Academic Publishers, Dordrecht, pp. 51–58.
Freund, L.B., 1993. Stress distribution and curvature of a general compositionally graded semiconductor
layer. J. Crystal Growth 132 (1–2), 341–344.
Frost, H.J., Ashby, M.F., 1982. Deformation-Mechanism Maps. Pergamon Press, Oxford.
Giannakopoulos, A.E., Suresh, S., Finot, M., Olsson, M., 1995. Elastoplastic analysis of thermal cycling;
layered materials with compositional gradients. Acta Metall. Mater. 43 (4), 1335–1354.
Hashin, Z., Shtrikman, S., 1963. A variation approach to the theory of the elastic behaviour of multiphase
materials. J. Mech. Phys. Solids 11, 127–140.
Howard, S.J., Tsui, Y.C., Clyne, T.W., 1994a. The e@ect of residual stresses on the debonding of coatings—I:
a model for delamination at a bimaterial interface. Acta Metall. Mater. 42 (8), 2823–2836.
Howard, S.J., Tsui, Y.C., Clyne, T.W., 1994b. The e@ect of residual stresses on the debonding of coatings—
II: an experiment study of a thermally sprayed system. Acta Metall. Mater. 42 (8), 2837–2844.
Jin, Z.H., Batra, R.C., 1996. Some basic fracture mechanics concepts in functionally graded materials. J.
Mech. Phys. Solids 44 (8), 1221–1235.
Kroner, E., 1961. Zur plastischen verformung des vielkristalls. Acta Metall. Mater. 9, 153.
Lee, Y.D., Erdogan, F., 1994=1995. Residual=thermal stresses in FGM and laminated thermal barrier coatings.
Int. J. Fract. 69, 145–165.
Li, G., Ponte Castan˜eda, P., 1994. Variational estimates for the elastoplastic response of particle-reinforced
metal–matrix composites. Appl. Mech. Rev. 47 (1), S77–94.
Miyamoto, Y., 1996. The applications of functionally graded materials in Japan—from structure to functional
uses. Mater. Technol. 11 (6), 230–236.
Olsson, M., Giannakopoulos, A.E., Suresh, S., 1995. Elastoplastic analysis of thermal cycling: ceramic
particles in a metallic matrix. J. Mech. Phys. Solids 43 (10), 1639–1671.
Ponte Castan˜eda, P., 1991. The e@ective mechanical properties of nonlinear isotropic composites. J. Mech.
Phys. Solids 39 (1), 45–71.
Shaw, L.L., 1998. The crack driving force of functionally graded materials. J. Mater. Sci. Lett. 17, 65–67.
Talbot, D.R.S., Willis, J.R., 1985. Variational principles for inhomogeneous nonlinear media. IMA J. Appl.
Math. 35, 39–54.
Weissenbek, E., Pettermann, H.E., Suresh, S., 1997. Elasto-plastic deformation of compositionally graded
metal–ceramic composites. Acta Mater. 45 (8), 3401–3417.
Wetherhold, R.C., Seelman, S., Wang, J., 1996. The use of functionally graded materials to eliminate or
control thermal deformation. Compos. Sci. Technol. 56, 1099–1104.
Yang, Y.Y., Munz, D., 1997. Stress analysis in a two-materials joint with a functionally graded material. In:
Shiota, I., Miyamoto, Y. (Eds.), Functionally Graded Materials 1996. Elsevier Science BV, Amsterdam,
pp. 41–46.