Progress in Materials Science 46 (2001) 379±405

www.elsevier.com/locate/pmatsci

Modeling the mechanical behaviour of

heterogeneous multi-phase materials
David S. Wilkinson a,*, Wolfgang Pompe b, Matthias Oeschner b
a
Department of Materials Science and Engineering, McMaster University, Hamilton, ON L8S 4L7, Canada
b
Technische UniveritaÈt Dresden, Dresden, Germany

Abstract
Many materials of engineering interest have highly heterogeneous microstructures. To a
®rst approximation, the response of multi-phase materials to external stimuli such as
mechanical loading depends on global parameters such as average particle size or phase
volume fraction. Most classical models of materials behaviour are based on such an assump-
tion. It is clear however that an accurate description must include parameters that character-
ize the distribution of phases. Moreover, some processes that we wish to model are inherently
stochastic in nature. This adds considerable complexity. First, the quantitative description of
microstructure containing higher order moments is fraught with diculties Ð both analytical
and experimental. Second, the inclusion of clustering into analytical models is prone to
assumptions and approximations. In this paper we will restrict ourselves to phenomena for
which a continuum approach is adequate. For these, self-consistent approaches are especially
valuable. The two examples that we discuss in some depth are related to (i) damage in porous,
brittle ®lms such as thermal barrier coatings and (ii) the simultaneous e€ects of damage and
particle clustering on the elasto-plastic response of metal matrix composites. # 2001 Elsevier
Science Ltd. All rights reserved.

Contents

1. Introduction......................................................................................................380

2. Approaches to self-consistent modeling............................................................381

3. Linear properties...............................................................................................383

* Corresponding author. Tel.: +1-905-525-9140; fax: +1-905-528-9295.

E-mail address: wilkinso@mcmail.cis.mcmaster.ca (D.S. Wilkinson).

0079-6425/01/$ - see front matter # 2001 Elsevier Science Ltd. All rights reserved.
PII: S0079-6425(00)00008-6
380 D.S. Wilkinson et al. / Progress in Materials Science 46 (2001) 379±405

4. Calculation of average phase stresses ...............................................................384

5. Damage development in porous, brittle solids under compression ..................386

5.1. The physical model ...................................................................................387
5.2. Non-linear macroscopic stress±strain relationship....................................389
5.3. Bifurcation phenomena.............................................................................393

6. Elastic±plastic ¯ow in two-phase solids ............................................................396

6.1. Modeling the role of spatial heterogeneity ...............................................396
6.2. Incorporation of damage in two-phase ¯ow models.................................399

7. Conclusions.......................................................................................................404

Acknowledgements................................................................................................404

References .............................................................................................................404

1. Introduction

With the possible exception of high purity, single crystals, the microstructure of a
real material is invariably heterogeneous. Even in single-phase materials the size and
orientation of grains is distributed nonuniformly. In multi-phase materials, the spatial
relationship between phases, and the size and orientation of particles can also be
distributed heterogeneously throughout the structure.
For most of its history, materials science has been content to describe micro-
structure in terms of average or global properties such as average grain size, overall
density or volume fraction of second phase. Moreover, most modeling of materials
behaviour has been based solely on such properties. These approaches have been
and continue to be useful, up to a point. Thus the density of a composite is de®ned
precisely by the average density of the constituent phases weighted by their volume
fraction. On the other hand the yield stress due to an array of small, hard particles
can be described, but only approximately, by the Orowan stress

Gb
ˆ 1†
l

using the average particle spacing l. To be more precise we must realize that the
spatial distribution of particles on the glide plane a€ects the resistance to dislocation
glide such that particle clustering for a ®xed average spacing lowers the yield stress.
To understand this one can consider what would happen if the particles were to
arbitrarily divide themselves into tight clusters, each of which contains say four
touching particles. These now behave as a single e€ective particle and the e€ective
average particle spacing increases by a factor of two, thus cutting the ¯ow stress in
two. This is of course an extreme example. (For a random distribution of particles of
 D.S. Wilkinson et al. / Progress in Materials Science 46 (2001) 379±405 381

uniform size, the ¯ow stress is reduced by about 15%). Clustering is rarely so severe,
and thus for properties such as this, which depend only weakly on higher order
moments of the dominant microstructural characteristics, the e€ect on the resulting
behaviour is modest.
At the other end of the spectrum are those material properties that clearly depend
on the extreme values of a distribution. These are generally related to material
breakdown emanating from defects or ¯aws in the microstructure. Common exam-
ples include the mechanical failure of brittle solids and the electrical breakdown of
dielectric materials. A well-developed methodology has been developed for dealing
with these properties based on both statistical approaches (Weibull analysis) and the
mechanics of ¯aw propagation.
Of greater concern in this paper are the wide varieties of properties that depend in
an intermediate way on the key microstructural parameters. There are many exam-
ples of this type. The elastic modulus of polycrystalline materials with a high degree
of elastic anisotropy depends on the degree of crystallographic texture. Similarly the
modulus of multi-phase materials containing elongated particles depends on the
orientation distribution and is generally anisotropic, even if those of the constituent
phases are not. This can be easily modeled for a unidirectional continuous ®bre
composite using bounding methods based on equal strain or equal stress in all con-
stituents (leading to the Voight and Reuss bounds). In general however, the actual
behaviour lies between these bounds. The prediction of yield stress in single-phase
polycrystals is another example, one with a long history, beginning with the seminal
papers by Taylor [1] and Bishop and Hill [2] who developed upper and lower bound
solutions. This work was later extended by KroÈner [3] and Hill [4] through the
development of self-consistent mechanics. In the period since then the self-consistent
method has been further re®ned and widely applied to a range of problems involving
both linear and non-linear phenomena. This method has proven to be extremely
valuable in estimating the mechanical and functional response of systems with het-
erogeneous microstructures. In the current age of ever-increasing computational
power, it is worth asking if this approach still has value, or whether problems
of sucient microstructural complexity can now be handled by large-scale ®nite
element approaches. It is our contention that analytical approaches will play an
important role in the modeling of materials behaviour for some time to come.
In the following we present a number of examples to support this hypothesis. We
start with a very brief survey of some classes of linear behaviour. However, we devote
the bulk of the paper to two non-linear phenomena. The ®rst concerns the development
of damage in a highly defected material while the other concerns incorporating the
e€ects of inhomogeneity and damage into the modeling of elasto-plastic behaviour.

2. Approaches to self-consistent modeling

The self-consistent method of continuum mechanics is grounded in solutions for

the behaviour of an ellipsoidal inclusion in a matrix. Such solutions were originally
developed by Eshelby [5], using what became known as the ``Eshelby equivalent
382 D.S. Wilkinson et al. / Progress in Materials Science 46 (2001) 379±405

inclusion method''. In Eshelby's original work the solutions are only valid if the
inclusion is surrounded by an in®nite matrix. Various approaches have since been
developed to incorporate the e€ect of a ®nite interparticle spacing. The earliest of these,
due to KroÈner [3], involved the suggestion that an equivalent homogeneous medium
could be used to represent the composite whole. The properties of this material were to
be found through an incremental self-consistent process. KroÈner later developed an
elastic-plastic method for polycrystals valid at small strain [6], which was further
re®ned by Hill [4]. This approach has been used to predict the response of both
polycrystal and multi-phase materials by Hutchinson [7] and Berveiller and Zaoui
[8]. Recent extensions of this ®eld have been proposed by Molinari et al. [9] and by
GonzaÂlez and Llorca [10], while a new approach based on variational principles has
been developed by Ponte-CastanÄeda and Suquet [11].
We are particularly interested in this paper in considering the behaviour of solids
containing heterogeneities, which may be a second phase of some form or a dis-
continuity associated with damage (microcracks, porosity, etc.). We will treat all such
heterogeneities as a second ``phase'' using a broader de®nition of this term that is the
norm. Depending on the nature of the system under study it may be possible to iden-
tify one phase as a matrix and the other as inclusions. Inclusions are often hard,
elastic particles but this is not necessarily so. In some cases the two phases are clearly
interconnected and there is no obvious matrix. This is important as two rather dif-
ferent approaches to self-consistency have been developed over time. The ®rst
approach is known as either the ``classical'' self-consistent method or the e€ective
medium approach (EMA). In this approach both phases are independently modeled
as inclusions sitting in an in®nite, homogeneous equivalent continuum (see Fig. 1).
A second approach known as the e€ective ®eld approach (EFA) or the ``general-
ized'' self-consistent method assigns a distinct inclusion phase to be embedded in a
matrix which then resides in the in®nite, homogeneous equivalent continuum [12]. This
is in fact an extension of Hashin's composite sphere model [13]. As one might expect,
the EMA model produces a sti€er response than the EFA model, especially when the
``inclusion'' phase is considerably harder or sti€er than the ``matrix''. This is because for
the EFA model the hard inclusions are assumed to be separated from one another by a
layer of the matrix phase at all volume fractions. On the other hand in the EMA
model some level of interconnectivity is assumed at all volume fractions. Neither one
of these assumptions is strictly correct. Thus, the choice of model should be made
with some understanding of the microstructure of the materials under study.
A second variant in modeling elastic-plastic behaviour follows from the treatment
of the rate of non-linear strain hardening. One can either use the secant modulus C  ˆ
=" or the tangent modulus T  ˆ @=@". The latter approach lends itself to incre-
mental methods, which is important for some problems, as we will see later. How-
ever, tangent modulus constructions tend to be rather sti€. Recently, GonzaÂlez and
Llorca [10] have shown however, that for proportional loading, the isotropic form of
the tangent modulus can be used and gives results that compare well with both
secant modulus and ®nite element calculations.
There have also been attempts to simplify some models by reducing the full stress
tensor to scalar form. This is especially valuable in the studies of inhomogeneous
 D.S. Wilkinson et al. / Progress in Materials Science 46 (2001) 379±405 383

Fig. 1. Two versions of self-consistent model, referred to as (a) the e€ective medium approach (EMA) and (b)
the e€ective ®eld approach (EFA), make di€erent assumptions about the connectivity of the two phases.

damage, as outlined below. There are however signi®cant concerns as to the validity
of such models. For this reason we have developed a three-dimensional model for an
elastic±plastic material based on a Mori±Tanaka approach [14]. For the case of
axisymmetric loading this gives good agreement with the simpler models.
To summarize, there are a variety of approaches to the development of self-con-
sistent models. However, these can generally be reconciled and shown to produce
similar results. In the following we present a number of studies in which this method
has been used to model the response of inhomogeneous solids involving both linear
and non-linear phenomena.

3. Linear properties

The most straightforward application of self-consistent methods is to problems in

which the response is linear. Examples include elasticity, thermal expansion and
viscosity. Consider for example, the elastic behaviour of porous solids. In this case
384 D.S. Wilkinson et al. / Progress in Materials Science 46 (2001) 379±405

we must use the EFA method since one of the component has zero sti€ness. The
e€ective modulus can be shown to equal
 
1 vp
E  ˆ Eo 2†
1 ‡ vp TE 1†

where vp is the pore volume fraction. TE is a parameter which depends only on

Poisson's ratio and pore shape:

1
TE ˆ ‰ 1 2†TK ‡ 1 ‡ †TG Š
3

where the functions TK and TG are given by Kreher and Pompe [15] as
  1   1
K K G G
TK ˆ 1 ‡ ; TG ˆ 1 ‡ 3†
KE GE

with

4 15  1 v
KE ˆ K ‡ G ; GE ˆ G 4†
3 2 4 5v

Note that the superscript * is used throughout to denote e€ective bulk properties.
TE is approximately equal to 2 for spherical pores.
Similar results can be obtained for composite materials containing a distribution
of hard particles. A good example involves aluminum reinforced with SiC particles.
Because of the large elastic modulus mismatch between these two materials (a factor
of close to 6) the e€ect is signi®cant. Fig. 2 shows a comparison of data for a range
of aluminum alloys containing up to 40 vol.% SiC particles, with a self-consistent
calculation. In this case an EMA model has been used. However, because of the
linear response such properties do not di€erentiate very strongly between models
and other approaches would provide a similar level of agreement.

4. Calculation of average phase stresses

Self-consistent methods can also be used to help understand the redistribution of

stress in a composite material as a result of elastic or elasto-plastic loading. Consider
for example the case of ceramic materials reinforced with platelets. It is generally
observed that the addition of platelets increases the fracture toughness of a material,
but only at the expense of a considerable loss of strength. This may occur for a
variety of reasons. For example, the increased complexity associated with processing
of composites can lead to an increase in defect density. In materials which fail by the
propagation of Grith cracks this can reduce the strength considerably. However,
we must also consider the e€ect of di€erences in thermal expansion coecient on the
 D.S. Wilkinson et al. / Progress in Materials Science 46 (2001) 379±405 385

Fig. 2. The tensile modulus of several SiC-reinforced aluminum alloys, normalized by the modulus of the
unreinforced alloy is plotted as a function of the SiC particulate content. An EMA self-consistent model
predicts the data rather well.

internal stress ®eld within the matrix. For example SiC has a CTE which is con-
siderably smaller than that of alumina but larger than that of silicon nitride. On
cooling from the ®ring temperature an internal stress ®eld is developed such that the
isotropic component of the stress tensor takes the form [15,16]

p m †
T vp
p ˆ 3vp ; m ˆ p 5†
vp vm 3 vm
‡ ‡
Km Kp 4Gm

where  is the stress, v the volume fraction, the coecient of thermal expansion, G
the shear modulus and K the bulk modulus, while subscripts p and m represent the
platelet and matrix, respectively. To this can be added a shape dependence which
results in a variation of the stress ®eld in the matrix around the (assumed elliptical)
platelet [17]. Consider the case of SiC-reinforced Si3N4, in which the matrix stresses
following cooling are tensile (Fig. 3). If one attempts to use the fracture strength and
toughness data to determine a critical ¯aw size in the material it appears to increase
by factor of 7 (from 30 to 204 mm) when up to 30 vol.% SiC is added. It seems
unlikely that such a large increase in processing defects should occur through the
addition of 25 mm platelets. Our analysis however [16], suggests that the residual
matrix stress on cooling is about 300 MPa. When this is added to the applied stress
to determine an e€ective stress the critical ¯aw size increases only by a factor of 2
upon the addition of platelets, a much more reasonable value. Similar results have
been obtained in other systems. This work suggests that the design of a suitable
386 D.S. Wilkinson et al. / Progress in Materials Science 46 (2001) 379±405

Fig. 3. The stress predicted for a SiC platelet-reinforced Si3N4 composite at di€erent positions, as a
function of the platelet aspect ratio.

residual stress ®eld is extremely important in the design of platelet-reinforced cer-

amics and similar composite materials.

5. Damage development in porous, brittle solids under compression

Non-linear stress±strain behaviour can be observed in brittle materials, even when

the material itself behaves elastically on the microscale. This non-linearity arises
from changes in the defect structure with increasing load, mainly due to micro-
cracking. One important example is the damage that develops in brittle coatings
such as thermal barrier coatings (TBCs). Usually the processing of a TBC results in
substantial porosity and a high density of microcracks which contributes to the high
compliance of these coatings. However with increasing load Ð in particular com-
pressive residual stress due to thermal cycling Ð damage evolution occurs which
®nally can cause spallation of the coating. In the following we show that a self-
consistent model incorporating a Weibull statistics approach for the single micro-
crack event can be used to describe the evolution of damage, as well as to derive a
criterion for materials instability.
TBCs show some obvious di€erences compared with structural ceramics that
exhibit high strength or high fracture toughness. The microstructure has been
developed to provide high compliance and low thermal conductivity, resulting in
high thermal shock resistance. Therefore, damage tolerant microstructures with a
high density of micropores and microcracks are favoured. The thermal conductivity
of the ceramics can be reduced by increasing the porosity of the coating. Thermal
conductivities less than 1 W/m K can be reached from thermal barrier coatings with
porosities up to 25%. This high defect density is connected with signi®cant changes
 D.S. Wilkinson et al. / Progress in Materials Science 46 (2001) 379±405 387

of the microstructure before failure as shown experimentally by Wesling and Socie

[18] and Johnsen et al. [19] for instance. An increase of the Young's modulus has
been observed under compressive load, which could be explained by a cumulative
conversion of micropores to microcracks. Therefore, we surmise that the increasing
accumulation of microcracks may lead to a macro instability, schematically shown
in Fig. 4.
One approach to analyzing the growth of pre-existing microcracks emanating
from pores has been presented by Ashby and Hallam [20] and by Sammis and Ashby
[21]. Their analysis is in an analytical formulation based on a homogeneous defect
distribution, using beam theory to describe the ligaments of solids between the
arrays of cracks. The constitutive equations derived by Sammis and Ashby [21] were
the starting point for a bifurcation analysis of porous, brittle solids under biaxial
compressive load given by Wimmer and Karr [22].

5.1. The physical model

The current model extends this analysis by incorporating some statistical features
of the microstructure into the damage evolution model. The material consists of
three ``phases'' Ð matrix, pores, and microcracks. As an idealization, the pores are
assumed to be spherical with a uniform radius, ap . They are randomly distributed. If
we apply biaxial compressive or tensile strain to the material, at some critical strain
the pores collapse and penny-shaped microcracks develop. Thus, in this model,
microcracks are assumed to develop from pre-existing pores. The orientation of the
microcracks is also assumed to be random (Fig. 4). The conversion of microcracks
to pores is considered to be a statistical phenomenon, which can be described by
means of Weibull statistics. The Weibull distribution relates the ratio of microcrack
density Nc to the initial pore density Np0. Because in TBC microstructures there is a
high density of such micropores and microcracks, the probability of the formation
of new microcracks is in¯uenced not only by the load concentration at a single

Fig. 4. Damage development due to microcrack evolution out of pre-existing pores under compressive
and tensile load.
388 D.S. Wilkinson et al. / Progress in Materials Science 46 (2001) 379±405

defect (micropore) but also by the interaction with neighbouring defects. This can be
included by applying an e€ective inclusion approach such as that proposed by pro-
posed by Pompe and Kreher [23]. As shown in Fig. 5, we calculate the fracture
probability of a spherical ceramic inclusion embedded in an e€ective medium con-
sisting of pores and microcracks randomly distributed in the ceramic matrix.
The mechanical behaviour of this e€ective medium is described by the macro-
scopic bulk modulus, K*, and the macroscopic shear modulus, G*, whereas the
ceramic inclusion is characterized by the corresponding values of the dense ceramic,
Km, Gm. The average local stress ®eld, ij ; in the matrix inclusion can be determined
by solving the elastic problem for a spherical, single inclusion embedded in the
e€ective medium [24]. Following Kreher and Pompe [15], the average local stress in
the matrix phase can be calculated by
  1
1 K Km 4G
ij ˆ  kk ij 1 ‡
3 Km 3K ‡ 4G
  6†
1 G Gm 9K ‡ 8G
‡  ij  kk ij †
1 ‡
3 Gm 5
3K ‡ 4G †

where  ij denotes the applied macroscopic stress, and ij is the Kronecker Delta.
Since the applied load to the solid is not uniaxial, an invariant representing the
applied load has to be developed for use in the Weibull theory. One option is the
elastic strain energy density, ; calculated with the average local inclusion stress due
to the applied macroscopic load. This has the form

1
 ˆ ij "ij 7†
2
with

ij ˆ Cm
ijkl "kl 8†

where Cm
ijkl denotes the elastic tensor of the inclusion phase.

Fig. 5. Spherical inclusion (ceramic) embedded in an e€ective medium (ceramic+pores+microcracks).

 D.S. Wilkinson et al. / Progress in Materials Science 46 (2001) 379±405 389

Therefore, the Weibull probability relates the microcrack density to the pre-exist-
ing pore density according to:
(   =2 !)

Nc ˆ Np0 1 exp 9†
0

where is the Weibull modulus and

1 2
0 ˆ  10†
2E0 0

is the characteristic strain energy density for fracture. Here E0 and 0 represent the
Young's modulus and the Weibull parameter under uniaxial loading of the inclusion
respectively. Note that Eq. (9) transforms to the well-known Weibull equation when
uniaxial load is applied.

5.2. Non-linear macroscopic stress±strain relationship

We assume a random distribution of pores and microcracks in the ceramic coat-

ing. Therefore, with increasing volume fractions of micropores, there is a growing
probability to ®nd interconnecting pores or microcracks. The constitutive equations
of such a structure consisting of ceramic matrix, pores, and microcracks, can be
developed through a self-consistent model based on the e€ective-medium approach
(EMA). In the vicinity of every defect (e.g. pores or cracks) all components of the
heterogeneous material are found with a probability corresponding to the various
volume fractions. Thus the formation of interconnecting micropores as well as
microcracks can be modelled. This percolation behaviour cannot be described using
the e€ective ®eld approach (EFA). In the latter model, micropores as well as
microcracks are assumed to be embedded always in the ceramic matrix. Since the
TBC microstructure exhibits microcrack coalescence, resulting in the formation of a
macro instability, we use in the following the EMA.
Following the notation of Kreher and Pompe [15], the e€ective Hooke's tensor C ,
can be determined by solving the following set of implicit equations for the e€ective
bulk modulus, K , and shear modulus, G , respectively:

< K K †TK K; G; ; K ; G † > ˆ 0 11†

vm
Km K †
TK Km ; Gm ; 1; K ; G † vp
K
TK 0; 0; 1; K ; G †
12†
vc
K
TK 0; 0; c ; K ; G † ˆ 0
and
< G G †TG G; K; ; G ; K † > ˆ 0 13†

vm
Gm G †
TG Gm ; Km ; 1; G ; K † vp
G
TG 0; 0; 1; G ; K †

14†
vc
G
TG 0; 0; c ; G ; K † ˆ 0
390 D.S. Wilkinson et al. / Progress in Materials Science 46 (2001) 379±405

where the tensor T is derived from the elastic polarization of the inclusions. To cal-
culate this, the Eshelby [5] assumption of an ellipsoidal shape has been used. The
shape parameter ˆ aa3 describes the ellipsoidal aspect ratio of the length of the
rotational axis of the inclusions ( =1: sphere; ! 1: oblate like crack; ! 0:
®ber like inclusion). TG and TK are given by Eqs. (3) and (4). v denotes the volume
fraction and the subscripts, m, p, and c indicate matrix, pore, and microcrack,
respectively. The volume fractions associated with each component are given by

4
vp ˆ  Np a3p
3

and

4 1
vc ˆ  Nc a3c
3 c

where N denotes the number density of each inclusion and the subscripts p, c refer to
pores and microcracks, respectively. For penny-shaped microcracks, ac3 ! 0 and
it is more convenient to substitute vc by the ®nite microcrack density, !:

4
! ˆ vc c ˆ Nc a3c
3

Following Kreher and Pompe, vc.TG and vc
TK can be derived by

 8 1  † 5  †
lim vc
TG 0; 0; c; G ; K † ˆ

! 15†
c !1 15 2 
vc ! 0

and

 4 1 2 †
lim vc
TK 0; 0; c; K ; G † ˆ

!; 16†
c!1 3 1 2
vc ! 0

respectively.
Furthermore, we introduce the parameter , which represents the ratio of the
e€ective volume dominated by a single microcrack to a single pore.
 3
ac
ˆ : 17†
ap

We assume the parameter always to be larger than 1 and to be dictated by the

microstructure of the given material.
Under hydrostatic pressure, microcracks are closed. Therefore, the third term in
Eq. (10) vanishes. The microcracks only in¯uence the macroscopic elastic constants
 D.S. Wilkinson et al. / Progress in Materials Science 46 (2001) 379±405 391

via load concentration under local shear. Thus, the two implicit equations can be
solved, the result being:

Km K
1‡ 2 3 1
4  
K ‡ G K 
!

K ˆ Km 3 K
6
41
7
5
vp0 18†
! 4
1 vp0 ‡ K ‡ G
3

and

Gm G
1‡ 20 1
5  3K ‡ 4G 
1

G G 
!
G ˆ Gm 6 K‡ 2G G
6B1 C
4@ A
vp0
1 5  3K ‡ 4G
1 vp0 !
1 G
6 K ‡ 2G

8 ‰1  Š
‰5  Š
‡
!
15 ‰2  Š
19†

with

1 3K 2G
 ˆ
2 3K ‡ G

Note that Np0 ˆ Np ‡ Nc ˆ constant:

Under tensile load, the hydrostatic stress also causes load concentrations at the
microcrack. Therefore, we get an asymmetric macroscopic stress±strain behaviour.
For the evaluation of these equations, for !c , the actual values due to the Weibull
equation [Eq. (9)] can be used [24].
Figs. 6 and 7 show the dependency of the e€ective bulk modulus, K , and the
e€ective shear modulus,G , respectively on the applied compressive biaxial strain at
di€erent values of the parameter . The e€ective bulk modulus increases up to that
of the matrix phase for all given values of since the pores vanish totally with
increasing applied load. The behaviour of the e€ective shear modulus, G* (Fig. 7) is
strongly dependent on . For small values the shear modulus increases, whereas for
larger values it decreases. The larger the ratio of microcrack to pore radius, the more
the volume ratio of the damaged material due to microcracks increases, as a result of
which the microcrack phase percolates more readily at a given number of micro-
cracks. Developing microcracks must have a certain radius related to the pore radius
in order to decrease the shear modulus. The pore density and the microcrack density
versus the applied strain are plotted in Fig. 8.
As stated above, the constitutive behaviour for the damaged material containing
cracks and pores is di€erent under tensile load. Although the e€ective shear modulus
392 D.S. Wilkinson et al. / Progress in Materials Science 46 (2001) 379±405

Fig. 6. E€ective bulk modulus versus the applied biaxial compressive strain at an initial porosity of 10%
( =6). Note that o= o/E.

Fig. 7. E€ective shear modulus versus the applied biaxial compressive strain at an initial porosity of 10%.
Note that o= o/E.

behaves equally under tension or compression, the e€ective bulk modulus di€ers.
Fig. 9 shows the e€ective tensile bulk modulus, KT , related to the biaxial tensile
strain at given initial porosity. The asymmetric stress-strain behaviour is well estab-
lished experimentally (e.g. [19]).
 D.S. Wilkinson et al. / Progress in Materials Science 46 (2001) 379±405 393

Fig. 8. Evolution of the porosity and microcrack density versus the applied biaxial compressive strain.
Note that o= o/E.

Fig. 9. Calculated values of the e€ective bulk modulus versus the biaxial tensile strain. Note that o= o/E.

5.3. Bifurcation phenomena

It is well known that non-linear stress±strain behaviour can be connected with a

change of the deformation mode under continuous loading. Microscopically, loca-
lized formation of interconnecting microcracks could occur when the microcrack
density approaches a critical value. This idea has been already employed by Wimmer
394 D.S. Wilkinson et al. / Progress in Materials Science 46 (2001) 379±405

and Karr [22] to describe compressive failure of microcracked porous brittle solids
with a regular micropore±microcrack array. A similar result can be derived from the
non-linear EMA approach of a random micropore±microcrack distribution.
A necessary condition for loss of material stability and any type of bifurcation is
given by be the loss of positive de®niteness of the rate of second order work [25±28].
This statement indicates that a necessary condition for loss of material stability is
: :

" ˆ 0; 20†

where
: :
 ij ˆ Dijkl " kl 21†

and the dot represents time derivation. The equivalent incremental form of Eq. (21) is:

@Cijkl
dij ˆ Cijkl d"kl ‡ d"kl
"mn ˆ Dijkl d"kl 22†
@"kl

When we apply this criterion for shear band formation (i.e. discontinuous bifur-
cation) the stability of an incremental shear deformation along an interface shown in
Fig. 10 has to be studied.
The shear band formation will occur ®rst when the ®rst eigenvalue of the
incremental sti€ness Dijkl goes to zero. To evaluate this criterion the whole set of
non-linear stress±strain equations including the Weibull failure condition has to
be evaluated. This has been done using one further approximation to simplify the
coupling between the Weibull criteria and the complete stress±strain curve by
the assumption that the e€ective shear modulus depends only on the porosity and the
microcrack density such that [24]

G ˆ Gm
1 2vp0 !  "††† 23†

Evaluating the bifurcation condition it has been shown that for a given -value,
initial porosities below a certain value do not lead to shear band formation. That means
no percolating microcrack clusters can be formed. In Fig. 11, the critical volume ratios,
, are plotted versus the initial porosity, when a shear band will be formed.

Fig. 10. Shear band geometry.

 D.S. Wilkinson et al. / Progress in Materials Science 46 (2001) 379±405 395

Fig. 12 shows the dependency of the shear band orientation angle  on the applied
compressive strain for di€erent values of initial porosity. At high porosity, a limited
density of cracked micropores will already form a percolating damage path along
the direction (/4) of maximum shear stress, whereas at lower porosities, only
microcracks with orientation more perpendicular to the loading direction can lead
to shear band formation (at higher compressive load).

Fig. 11. Critical values of the e€ective crack to pore volume ratio , for discontinuous bifurcation versus
the initial porosity vp0.

Fig.12. Biaxial compressive strain versus the shear angle for di€erent values of initial porosity vp0.
396 D.S. Wilkinson et al. / Progress in Materials Science 46 (2001) 379±405

6. Elastic±plastic ¯ow in two-phase solids

6.1. Modeling the role of spatial heterogeneity

We now turn our attention to a problem that has received much attention in
recent years Ð namely the plastic response of two-phase alloys (including metal
matrix composites). We are referring here in particular to systems containing relatively
large particles (greater than a few microns).1 One of the key features of this problem is
the importance of spatial inhomogeneity. As particles become clustered together the
yield stress and rate of initial work hardening increase signi®cantly. At large strains the
e€ect is lost. Thus to a certain extent the e€ect of clustering can be thought of as
extending the elastic±plastic transition. Physically, we envisage a process whereby
yielding starts to occur in regions of the microstructure that have a low volume fraction
of the reinforcing phase. Evidence for that has been found by studying the distribution
of slip bands following small strain deformation of Al±SiC(p) composites [29]. (This
is analogous to the problem of microyield in polycrystals which commences in
crystals that have a soft orientation for dislocation glide). In either case, as defor-
mation proceeds, the regions that deform ®rst begin to work harden and the other
regions are then forced to partake in the overall process of plasticity. Eventually all
regions are hardening at a more or less equal rate and the global hardening rate of
the inhomogeneous solid becomes parallel to that of the more uniform material.
Self-consistent modeling lends itself well to this kind of problem. The spatial
inhomogeneity in the distribution of the second phase can be modeled by treating
the material as a multi-phase alloy, in which each ``phase'' represents a material with
a speci®c volume fraction of the reinforcing phase. We have mostly dealt with
bimodal distributions (see Fig. 13). However, the extension to a larger number of
phases is not dicult and we have used this to handle problems related to particle
size distributions [30], as well as the incorporation of damage into these models as
outlined in the next section. In all of these problems one ®rst starts by developing a
constitutive relationship for each ``phase''. In the case of a bimodal composite then,
we develop a constitutive law for each volume fraction. This can be experimental or
it can be obtained from a model, using either self-consistent approaches or single-
cell ®nite element methods. In all cases the results are ®t to a constitutive law of
standard form, such as the Ramberg±Osgood equation. This two step process can be
used since the constitutive behaviour of each region is a function only of its local
stress±strain state. Thus interactions between the regions of high and low volume
fraction lead to stress partitioning, but they do not a€ect the shape of the stress±
strain response. We have used this method extensively to study the mechanical
behaviour of particulate±reinforced metal matrix composites [14,31,32]. The main
conclusions are summarized in Fig. 14, which shows the ratio of limit stress in a
clustered composite normalized by that of a uniform composite of the same overall

1
When particles are very small (4 0.1 mm) these problems are best treated using classical dislocation
mechanics. For particles of intermediate size the emerging ®eld of strain-gradient plasticity o€ers the best
approach.
 D.S. Wilkinson et al. / Progress in Materials Science 46 (2001) 379±405 397

Fig. 13. An idealized model of a clustered composite consisting of two ``phases'', one containing a high
volume fraction, the other with a low volume fraction.

Fig. 14. The limit stress of a particulate composite normalized by that of the unreinforced material plot-
ted as a function of the degree of clustering. The calculation has been performed for two di€erent rates of
work hardening N.

volume fraction of reinforcing phase, plotted as a function of the degree of cluster-

ing. As the degree of clustering increases so does the strength of the composite, at
least initially. If this becomes too severe, the strengthening ratio actually starts to
decrease. This is because the strengthening response is partly related to the shielding
of matrix material that occurs inside the clusters. The clusters behave in part as large
particles that undergo less plastic deformation than the regions with lower particle
398 D.S. Wilkinson et al. / Progress in Materials Science 46 (2001) 379±405

volume fractions. As the severity of clustering increases the e€ectiveness of this

shielding also increases. However, the e€ective volume fraction occupied by clusters
decreases. Eventually when the clustering severity V2 approaches one each cluster
behaves like a single large particle. Since the mechanical response is scale-indepen-
dent the strength returns to that of an unclustered material. Having noted this, it is
important to realize that the severity of clustering in most materials of interest puts
us on the left side of this curve, such that clustering is always expected to increase
the limit strength of a material.
Elasto-plastic models provide a good check on the di€erences between various
approaches to self-consistent models. We have done this by comparing an EMA
calculation that uses both an incremental tangent modulus algorithm and a secant
modulus approach with an EFA calculation based on a secant modulus construction
(the latter using software developed by Bornert and Zaoui) [32,33]. As shown in Fig.
15a for the case of a homogeneous two-phase composite (Al with SiC particles), there
is very little di€erence between these two calculations until the volume fraction of the
reinforcing phase approaches 30 vol.%. Thereafter the di€erences are signi®cant. The
EMA model is inherently sti€er, especially when used in conjunction with a secant
modulus formalism. This is due in part to the di€erent assumptions about micro-
structural connectivity inherent in these two models and may be re¯ective of second
order e€ects that can be realized in actual microstructures. In other words, we expect
a two-phase alloy in which percolation of the hard particles is suppressed to high
volume fractions to behave in a softer manner than one for which this is not the case.
There is, however, no simple way to correlate the extent of percolation found in a
real microstructure with the predictions of the model. We simply note that the dif-
ference between these models may re¯ect some of the true scatter observed in
materials made of the same constituents but with di€erent processing routes. Of
particular importance is how the model a€ects predictions about the role of clus-
tering. In Fig. 15b we show a calculation for a 20 vol.% composite in which all of
the particles are clustered in regions which have a local particle density of 40 vol.%.
Using the EFA model with a secant modulus construction there is very little e€ect of
clustering. The composite sphere model, in which the hard ``phase'' regions (i.e. the
particle-containing clusters) are all surrounded by the softer material, shields each
cluster and reduces the redistribution of load. In the EMA model, however, the
e€ect is much stronger, since this implicitly assumes that the regions of high and low
volume fraction are fully interpenetrating throughout the microstructure.
One area of concern in dealing with the models just discussed is that they use a
simpli®ed one-dimensional stress analysis which ignores the additional terms that
arise for a full tensorial analysis. In order to address this we have developed a three-
dimensional model, which extends the work of Estevez et al. [34]. In this work [14],
both 1-D and 3-D models use the same computational scheme and vary only in the
simplifying assumptions used in the 1-D version. For uniaxial loading, we ®nd that
there is essentially no di€erence between the models at all volume fraction levels
studies, as indicated in Fig. 16. This would appear to justify the use of the simpler
models when analyzing the tensile response of materials, although the more complex
models are certainly of value for assessing ¯ow behaviour during non-proportional
 D.S. Wilkinson et al. / Progress in Materials Science 46 (2001) 379±405 399

Fig. 15. (a) A comparison of results for three versions of self-consistent analysis using both EFA and
EMA approaches. The EMA calculation has been performed using both a tangent modulus and secant
modulus construction. These calculations assume a homogeneous AA2618(T4)±SiC composite. (b) A
comparison of self-consistent models for a clustered composite consisting of 20 vol.% SiC with all of the
particles in 40 vol.% clusters.

loading procedures as in forming operations, as well for incorporating e€ects due

non-spherical particles shapes and so on.

6.2. Incorporation of damage in two-phase ¯ow models

The incorporation of damage into models such as those described above rep-
resents an interesting challenge, particularly if one is also concerned with including
400 D.S. Wilkinson et al. / Progress in Materials Science 46 (2001) 379±405

Fig. 16. A comparison of models based on the full stress tensor (3-D model) and on only the uniaxial
stress (1-D model) for several di€erent levels of particulate loading.

e€ects due to clustering at the same time. The systems that have been most heavily
studied in this regard are particulate±reinforced metal matrix composites. This is largely
because damage is prevalent in these materials and central to their tensile behaviour.
Damage tends to develop at relatively low strains and accumulate as deformation pro-
ceeds. The loss of sti€ness associated with damage has a signi®cant impact on the
overall work hardening rate of the material and thus a€ects both tensile stability and
ductility. Clear evidence of this is provided by an experiment in which an Al±Si alloy
was damaged by pre-straining in the T4 condition, heat treated to T6, and then
subjected to a tensile test. The pre-damaged sample showed greater ductility than a
virgin sample in the same heat treatment condition, because of the destabilizing
e€ect of damage [35].
For the most part, damage takes the form of particle cracking, although particle
decohesion does occur, especially in materials processed by powder metallurgy
routes. Since the particles are brittle their fracture is stochastic and can be analyzed
using Weibull statistics. A number of aluminum alloys containing SiC particles have
been studied experimentally in which the density of cracked particles is measured as
a function of global strain. This is related to macroscopic stress and then, through a
stress concentration model, to the average particle stress. From this can be extracted
the Weibull parameters describing the distribution of fracture strengths in the SiC
particles. There is surprisingly good agreement amongst the di€erent studies
[32,36,37] with a characteristic fracture stress of about 1.4±1.6 GPa and a Weibull
modulus of 4±5. The Weibull modulus is quite low, especially given that following
thermomechanical working many of the larger defects in the SiC particles are lost
due to particle fracture. This is actually re¯ective of the analysis, which assumes that
 D.S. Wilkinson et al. / Progress in Materials Science 46 (2001) 379±405 401

all of the particles in the composite are subjected to the same stress. It is quite likely
that the distribution in such stress is broader than the distribution of particle
strengths. We must therefore be cautious is claiming overly good agreement between
the models we describe below and experimental data. The good ®t is due in part to
``hiding'' some of the diculties in the Weibull parameter measurements. If one
could independently measure the strength distribution of the particulate phase then
the models would need to be much more sophisticated to describe the true micro-
structural complexity of these systems. Despite this proviso, this ®eld has seen con-
siderable progress in recent years, which we now describe.
We can use the methodology described in the previous section to model damage
development also, by considering the undamaged and fully damaged composite as
separate ``phases''. The added complexity arises from the need to allow the volume
fraction of these ``phases'' to vary with time during a test. In this case the use of an
incremental model involving a tangent modulus construction is essential. Beyond
this the approach is relatively straightforward. There are several key elements. The
®rst is to gather a sucient body of experimental data on the system in question to
obtain parameters that are not available otherwise. This includes the Weibull dis-
tribution for particle strength, as well as the particle size distribution. The latter is
important since the stochastic nature of brittle particle fracture is inherently depen-
dent on the size of the specimen. This means that while in continuum mechanics
plasticity is scale-independent once we introduce damage this feature of the model is
lost. The second element is to obtain proper constitutive laws for each ``phase''. In
the case of the damaged composite single cell ®nite element models have led to the
best solutions [38]. These elements can then be assembled into a self-consistent
model. In a clustered composite the model has 2N phases in which N is the number
of classes we use to characterize the clustering while the factor 2 allows for an
undamaged and damaged version of each class. Thus

X
N

fi U ‡ fi D ˆ 1
iˆ1

where

fiU ‡ fiD ˆ fiO

where, f O U D
i is the volume fraction required by phase i, while f i and f i are the volume
fractions that are undamaged and damaged respectively. At the onset of defor-
mation fiD is zero for all phases. After each increment of macroscopic strain the
increase of damage level in each phase is reevaluated in terms of the local stress and
the Weibull distribution. The introduction of damage requires an additional process
of load relaxation. It is assumed that this process is fast and essentially elastic. This
allows it to be treated as a second independent step in the incremental loading
algorithm. Although these processes add bookkeeping complexity they do not
change the overall nature of the self-consistent model. The e€ective modulus of such
a composite can be found according to
402 D.S. Wilkinson et al. / Progress in Materials Science 46 (2001) 379±405

* +
T T
ˆ0 24†
3
T ‡ T
2
where the angle brackets represent an average over all phases in the solid, weighted
by the respective volume fraction. This formulation is based on the e€ective medium
approach and assumes a value for Poisson's ratio of 0.5. This is the fundamental
relationship, which lies at the heart of the model.
The results of a typical calculation are shown in Fig. 17. This has been calculated
for a 2618 T4 aluminum alloy containing 15 vol.% SiC particles. It shows the true
stress±true strain curve for four di€erent cases: with and without clustering, and with
and without damage. There are several conclusions one can draw from this: the ®rst is
that in the absence of damage clustering does increase the strength of a composite; the
second is that damage removes the bene®cial e€ect of clustering. This is because
damage occurs preferentially in the more highly stressed high volume fraction clusters.
This prediction is susceptible to experimental veri®cation. Since damage occurs
much more readily in tension than in compression the e€ect of microstructural
inhomogeneity should be stronger when testing is done in compression. However,
the complexity of most commercial materials makes this a dicult experiment to
interpret. We have therefore taken a di€erent approach of making materials that
exhibit a bimodal distribution as assumed in the models. These are made by a pow-
der metallurgy route using the Al±CuAl2 system [39,40]. Current studies are on the
Al±ZrO2 and Al±Al2O3 systems. When tested in tension materials which contain a
bimodal distribution exhibit almost identical stress±strain curves as compared with

Fig. 17. The prediction of EMA self-consistent calculations based on AA2618(T4) containing 15 vol.% SiC
particles. Calculations have been performed with and without damage, and with and without clustering.
 D.S. Wilkinson et al. / Progress in Materials Science 46 (2001) 379±405 403

materials that are homogeneous, having the same overall volume fraction of the
reinforcing phase. When tested in compression however the clustered materials are
clearly stronger, as shown in Fig. 18. Moreover, a self-consistent model using the
EMA formulation predicts the increase in strength rather well.
Despite this success, we once again need to address the impact of using a simpli®ed
scalar approximation in place of the full stress tensor. Thus, Fig. 19 shows the results
of a recent calculation using the full stress tensor, along the same lines as that outlined
in the previous section but now involving damage [14]. In this we plot the rate of
damage evolution, de®ned as the fraction of cracked particles for both the 1-D and
3-D models along with experimental data. The agreement is excellent.
There is however one important characteristic of these materials that the models
just discussed do not predict with any degree of accuracy. That is the ®nal ductility.
Extensive experimental work [35,41] has shown that particulate composites generally
start to neck once the ConsideÁre criterion has been met. However, when this criter-
ion is used with the models just described the ductility prediction is generally much
higher than observed. The reason appears to be due to details of the damage linkage
process. Prior to the onset of tensile instability, interparticle cracking commences.
This additional damage process accelerates the loss of work hardening thus moving
the onset of tensile instability to lower strains. None of the models developed to date
treats the problem of damage linkage.

Fig. 18. Experimental data for two Al±15 wt.% Cu alloys, one with a homogeneous distribution of the
CuAl2 particles, the other with a bimodal distribution of 10 and 24 wt.% regions. The EMA model ®ts the
data for the clustered composite well.
404 D.S. Wilkinson et al. / Progress in Materials Science 46 (2001) 379±405

Fig. 19. A comparison of models using the full stress tensor (3-D) and only the axial stress (1-D), in which
damage is plotted as function of the macroscopic plastic strain.

7. Conclusions

We have attempted in this paper to outline a range of problems in the mechanical

behaviour of materials that are amenable to modeling using self-consistent
mechanics. While there are always several approaches to modeling any phenom-
enon, we believe that this approach is particularly well suited to problems which
attempt to incorporate microstructural complexity. Two particular cases considered
here are spatial inhomogeneity or stochastic damage processes.

Acknowledgements

It is a pleasure to acknowledge the discussions we have had with many colleagues

on problems of the type discussed here over several years, in particular Professor
David Embury (McMaster) and Dr. Eric Maire (INSA de Lyon). One of us
(D.S.W.) is especially indebted to Professor Michael Ashby for nurturing and
inspiring a young scientist. Partial funding of this work by the Natural Sciences and
Engineering Research Council of Canada is gratefully acknowledged.

References

[1] Taylor GI. Journal of the Institute of Metals 1938;62:1307±24.

[2] Bishop JFW, Hill R. Philosophical Magazine 1951;42:1298±307.
[3] KroÈner E. Zeitschrift Physik 1958;151:504±18.
 D.S. Wilkinson et al. / Progress in Materials Science 46 (2001) 379±405 405

[4] Hill R. Journal of Mechanics and Physics of Solids 1965;13:89±101.

[5] Eshelby JD. Proceedings of the Royal Society 1957;A241:376±96.
[6] KroÈner E. Acta Metallurgica 1961;9:155.
[7] Hutchinson JW. Proceedings of the Royal Society 1970;A355:101±27.
[8] Berveiller M, Zaoui A. Journal of Mechanics and Physics of Solids 2000;26:325±44.
[9] Molinari A, Ahzi S, Kouddane R. Mechanics of Materials 1997;26:43±62.
[10] Gonzalez C, Llorca J. Journal of Mechanics and Physics of Solids 2000;48:675±92.
[11] Ponte CastanÄeda P, Suquet P. Advances in Applied Mechanics 1998;34:171±301.
[12] Christensen RM, Lo KH. Journal of Mechanics and Physics of Solids 1979;27:315.
[13] Hashin Z. Journal of Applied Mechanics 1962;29:143±50.
[14] Estevez R, Maire E, Francoisi P, Wilkinson DS. European Journal of Mechanics 1999;18:785±804.
[15] Kreher W, Pompe W. Internal stresses in heterogeneous solids. Berlin: Akademie-Verlag Berlin,
1989.
[16] Pompe W, Wilkinson DS. In: Pichoir R, editor. Proc. 3rd Intl. Conf. on Advanced Materials and
Processes Paris: L'Editions de Physique 1993. p. 1889±94.
[17] Kreher W, Janssen R. Journal of the European Ceramic Society 1992;10:167±73.
[18] Wesling KF, Socie DF. Journal of the American Ceramic Society 1994;77:1863.
[19] Johnsen BP, Cruse TA, Miller RA, Brindley WJ. ASME Journal of Engineering for Gas Turbines
and Power 1995;117:305.
[20] Ashby MF, Hallam SD. Acta Metallurgica 1986;34:497.
[21] Sammis CG, Ashby MF. Acta Metallurgica 1986;34:511.
[22] Wimmer SA, Karr DG. Mechanics of Materials 1996;22:265±77.
[23] Kreher W, Pompe W. Strength of ceramics. In: Kaldis E. editor. Current types in materials science,
vol. 12. Amsterdam: North-Holland, 1985. p. 362±91.
[24] Oechsner, M. PhD Thesis, U. Karlsruhe, 2000.
[25] Drucker DCQ. Appl Math 1950;7:411.
[26] Hill R. Journal of Mechanics and Physics of Solids 1958;6:103.
[27] Rice JR. In: Koiter WT, editor. Theoretical and applied mechanics. Amsterdam: North-Holland,
1976.
[28] Neilsen MK, Schreyer HL. International Journal of Solids and Structures 1993;30:521.
[29] Corbin SF, Wilkinson DS. Acta Metallurgica et Materialia 1994;42:1319±27.
[30] Wilkinson DS, Maire E, Fougeres R. Materials Science and Engineering A262 1999;264-70.
[31] Corbin SF, Wilkinson DS. Acta Metallurgica et Materialia 1994;42:1311±8.
[32] Wilkinson DS, Maire E, Embury JD. Materials Science and Engineering 1997;A233:145±54.
[33] Bornert M, Zaoui A. Private communication, 1996.
[34] Estevez R, Hoinard G, Francoisi P. Acta Metallurgica et Materialia 1997;45:1567±84.
[35] Kiser MK, Zok FW, Wilkinson DS. Acta Metallurgica et Materialia 1996;44:3465±76.
[36] Maire E, Lormand G, Gobin PF, Fougeres R. Journal de Physique IV 1993;C7:1849±52.
[37] Lewis CA, Withers PJ. Acta Metallurgica et Materialia 1995;43:3685±99.
[38] Brockenbrough JR, Zok FW. Acta Metallurgica et Materialia 1995;43:11±20.
[39] Conlon K, Maire E, Wilkinson DS, Henein H. Metallurgical and Materials Transactions
2000;A31:247±60.
[40] Maire E, Wilkinson DS, Embury JD, Henein H. Metallurgical Transactions 1998;29A:2613±20.
[41] Llorca J, Poza P. Materials Science and Engineering 1994;A185:25.