Elementary Micromechanics of

Heterogeneous Media

Konstantin Z. Markov
St. Kliment Ohridski University of Sofia
Faculty of Mathematics and Informatics
5 blvd J. Bourchier, P.O. Box 48
BG-1164 Sofia, Bulgaria
kmarkov@[Link]

Chapter 1 in the collection:

Heterogeneous Media: Modelling and Simulation,

edited by Konstantin Z. Markov and Luigi Preziosi, Birkhauser Boston, 1999, pp. 1162.

ii

Heterogeneous Media: Modelling and Simulation

CONTENTS
1.1 Introduction
1.1.1 The Aim of Micromechanics .
1.1.2 The Idea of Homogenization .
1.1.3 Brief Historical Remarks . . .
[Link] Navier and Cauchy . .
[Link] Poisson and Mossotti .
[Link] Maxwell . . . . . . . .
[Link] Einstein . . . . . . . .
[Link] Smoluchowski . . . . .
1.1.4 The Outline of the Survey . .

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

2
2
4
7
8
9
11
15
17
18

1.2 The Homogenization Problem

1.2.1 The Representative Volume Element . . . . . . . . .
1.2.2 Definition of the Effective Conductivity . . . . . . .
1.2.3 Definition of the Effective Resistivity . . . . . . . . .
1.2.4 Energy Definition of the Effective Constants . . . . .
1.2.5 Consistency of the Definitions of Effective Properties
1.2.6 The Concentration Factors . . . . . . . . . . . . .
[Link] Ideal Thermal Contacts . . . . . . . . . . . .
[Link] The Wu Representation . . . . . . . . . . . .
[Link] Nonideal Thermal Contacts . . . . . . . . . .
1.2.7 The Elastic Case . . . . . . . . . . . . . . . . . . . .
[Link] The Basic Equations . . . . . . . . . . . . . .
[Link] The Effective Elastic and Compliance Tensors
[Link] The Elastic Concentration Factors . . . . .
1.2.8 The Effective Absorption Coefficient . . . . . . . . .

.
.
.
.
.
.
.
.
.
.
.
.
.
.

21
22
25
28
32
35
39
39
42
44
45
46
47
49
50

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

1.3 Some Basic Results

1.3.1 Weakly Inhomogeneous Media . . . . . .
[Link] The Definition . . . . . . . . . .
[Link] The Perturbation Expansion and
order Approximation . . . . . . .
[Link] The Third-order Approximation

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

53
. . . . . . . . 54
. . . . . . . . 54
the Second. . . . . . . . 55
. . . . . . . . 58

Micromechanics of Heterogeneous Media

iii

1.3.2 The Hashin-Shtrikman Estimates . . . . . . . . . . . .

[Link] The Bergman Formula . . . . . . . . . . . . . .
[Link] The Matheron Inequalities . . . . . . . . . . . .
[Link] The Wiener Bounds . . . . . . . . . . . . . . .
[Link] The Hashin-Shtrikman Bounds . . . . . . . . .
1.3.3 The Hashin Assemblage . . . . . . . . . . . . . . . . .
1.3.4 The Variational Estimates . . . . . . . . . . . . . . . .
[Link] The Dirichlet Principle . . . . . . . . . . . . . .
[Link] The Thompson Principle . . . . . . . . . . . .
[Link] The Beran Bounds . . . . . . . . . . . . . . . .
[Link] The Torquato-Milton Patameter . . . . . . . .
1.3.5 The Bounds in the Elastic Case . . . . . . . . . . . . .
1.3.6 Cross-properties Relations: The Levin Formula . . . .
1.3.7 The Elementary Bounds on the Effective Absorption
Coefficient . . . . . . . . . . . . . . . . . . . . . . . . .

60
61
62
63
64
66
68
68
69
70
75
77
80
83

1.4 The Single Inclusion Problem

85
1.4.1 Scalar ConductivityIntegral Equation . . . . . . . . 85
1.4.2 Scalar ConductivityEllipsoidal Inhomogeneity . . . . 86
1.4.3 Singular Inclusions and Nonideal Contacts . . . . . . 86
[Link] Singular Inclusions . . . . . . . . . . . . . . . . 92
[Link] Nonideal Contacts . . . . . . . . . . . . . . . . 93
[Link] The Superconducting Spherical Inhomogeneity . . . . . . . . . . . . . . . . . . . . . . . . . 94
[Link] The Resisitive Spherical Inhomogeneity . . . 96
1.4.4 Elastic Case . . . . . . . . . . . . . . . . . . . . . . . . 99
1.5 One-particle Approximations
1.5.1 Dilute Filler FractionScalar Conductivity . . . .
1.5.2 Dilute Filler FractionElastic Case . . . . . . . .
1.5.3 Self-consistent Scheme . . . . . . . . . . . . . . . .
[Link] The Basic Idea . . . . . . . . . . . . . . . .
[Link] An Alternative Derivation . . . . . . . . . .
[Link] The Percolation Phenomenon . . . . . . . .
[Link] Some Particular Cases of the Self-consistent
Scheme . . . . . . . . . . . . . . . . . . . .
1.5.4 The Differential Scheme . . . . . . . . . . . . . . .
[Link] The Basic Idea . . . . . . . . . . . . . . . .
[Link] Some Particular Cases . . . . . . . . . . . .

.
.
.
.
.
.

.
.
.
.
.
.

105
106
110
113
113
113
117

.
.
.
.

.
.
.
.

118
120
121
122

iv

Heterogeneous Media: Modelling and Simulation

1.5.5 The Effective Field . . . . . . . . . . . . . . . . . . . . 125
[Link] The Basic Idea . . . . . . . . . . . . . . . . . . 126
[Link] Alternative Derivations and Interpretations . . 127
[Link] Some Particular Cases . . . . . . . . . . . . . . 128
1.5.6 The Conductivity of a Medium with Non-ideal Interphase Contacts . . . . . . . . . . . . . . . . . . . . . . 133
1.5.7 The Self-consistent Scheme for the Absorption Problem 135

1.6 Elastic Properties of Polycrystals

1.6.1 The Self-consistent Scheme for Polycrystals . . .
1.6.2 Formulation in Compliances and Consistency of
Scheme . . . . . . . . . . . . . . . . . . . . . . .
1.6.3 Example: Elastic Moduli of a Cubic Polycrystal .
1.7 References

139
. . . 139
the
. . . 142
. . . 143
146

1
Elementary Micromechanics of
Heterogeneous Media

Konstantin Z. Markov
St. Kliment Ohridski University of Sofia
Faculty of Mathematics and Informatics
BG-1164 Sofia, Bulgaria

Abstract. The introductory and more elementary ideas and results of micromechanics of heterogeneous media are collected in the
survey. The central problem under discussion is homogenization.
It replaces such media by homogeneous ones, which behave macroscopically in the same way and possess certain gross effective properties. These properties are related in a complicated manner to the
prescribed internal structure of the medium and their evaluation, in
general, represents a profound challenge in any specific situation. A
brief historical survey is given, underlying the reappearance of essentially the same homogenization quest in numerous guises and
contexts over the last two centuries. Within the framework of the
volume-averaging approach the basic notions are introduced and some
of the central, now classical, results are then derived and discussed
such as perturbation expansions, Hashin-Shtrikmans bounds, variational estimates and Levins cross-property relation. A general one1

Heterogeneous Media: Modelling and Simulation

particle scheme for approximate evaluation of the effective properties

(in the static case) is detailed in its various implementations like selfconsistency, iterated limits and effective field. Illustrations concern
conductivity, elasticity, and absorption phenomena in heterogeneous
particulate media, as well as a simple self-consistent model for polycrystals homogenization.

1.1 Introduction
1.1.1 The Aim of Micromechanics
Continuum mechanics deals with ideal homogeneous materials. Its
aim is to describe their response to external exertions using appropriate constitutive relations. The latter generally are specified by means
of macroscopical experiments without microstructural considerations.
It is a trivial statement, however, that any sample of material is microscopically inhomogeneous, even if it appears homogeneous at some
natural scale of observation. Inevitably, therefore, a description of
any material in terms of continuum mechanics is an approximation,
and any experimental determination of constitutive behavior yields,
in fact, a relationship between the overall properties measured in
the experiment. This observation leads us to a fundamental and
widespread problem of science and technology, concerning micromacro interconnection, i.e. a proper and reliable determination of
the macroscopic (or large scale) behavior of a medium which exhibits microscopic (or small scale) heterogeneity, on the base of
the appropriate and available microstructural information. (The exact meaning of large and small depends on the specific problems
and media under study.) The oldest problem of such a type, fundamental in statistical physics, is description of matter in terms of
its molecular constituents. Here small corresponds obviously to
molecular dimensions.
Micromechanics, in general, deals with heterogeneous media for
which small has a certain intermediate dimension ` which is large
compared to molecular dimensions, but is small in macroscale. The
length ` is connected with the characteristic size of the heterogeneities
in the medium, say, with mean radius of inclusions, voids, fibers, the

1.1 Micromechanics of Heterogeneous Media

size of a crystallite in polycrystalline aggregates, etc. Though the

molecular effects are not present, as a rule, in such a scale, many of
the methods used in the appropriate theories have their direct origin in statistical physics and statistical mechanics, as we shall point
out below. The aim of micromechanics is just to relate the gross
macroscopical behavior of heterogeneous media, characterized with
the above mentioned length scale `, to the details of their microscopical constitution. The basic idea is that of homogenization, which
consists in a replacement of a piece of a microheterogeneous solid
by a homogeneous one which, from a macroscopical point of view
behaves in the same manner as this piece. This fundamental idea
will be briefly illustrated in the next subsection, before considering it
more rigorously and in needed detail in Section 1.2.
Note that an important class of heterogeneous media, extensively
treated by micromechanics, are the compositesman-made mixtures
mixtures of two or more constituents, firmly (as a rule, but not always) bonded together. Among other reasons, technological demands
have stimulated a very extensive bulk of studies, devoted to their
gross mechanical behavior, like elasticity, strength, etc., as predefined
by a specific microstructural arrangement. A considerable part of the
more fundamental works along this line will be discussed in what follows. There are, however, plenty of micromechanical problems not
directly concerned with mechanical properties, but representing particular cases of the same micro-macro interconnection quest, in
which the same idea of homogenization applies.1 An example is furnished by the problem of describing propagation of electromagnetic
or acoustic waves through a heterogeneous solid or a turbulent fluid.
Another example is cooling of a heated body by means of a system of
tubes or holes, kept at a fixed temperature. Mathematically the latter problem is very closely connected to describing the gross behavior
of a diffusing species in a medium, containing absorbing sinks. The
1
It is more precise perhaps to speak about micro-physical problems, instead
of micromechanical since, say, diffusion, electromagnetic waves, etc., have little to
do with mechanics in its classical sense. However, introducing the term physics
would inevitably enlarge this survey far beyond the admissible length and (more
important)the authors competence. In this connection it seems appropriate to
recall the Constitution of the International Society for the Interaction of Mechanics
and Mathematics (ISIMM), where one reads that Mechanics is understood here
in the broad sense of the word, including relevant physical phenomena such as
electromagnetic and thermal fields.

Heterogeneous Media: Modelling and Simulation

third example to be specially mentioned concerns the viscous flow

through a porous solid (which itself could deform due to the flow).
1.1.2 The Idea of Homogenization
To illustrate the basic ideas of homogenization as simply as possible,
consider a (linear) elastic heterogeneous medium. Let, for example,
the tensile stress-strain behavior of the medium along the axis x1 be
under investigation. Imagine to this end that a large, say, cubical
specimen V with a side L is cut out from the medium; large means
here that L  `, where ` is the above mentioned microscale length.
Either the cube would be loaded to some level and its extension measured, or else it would be extended by some amount and the load
measured. Then the stress component 11 would be taken as load
divided by area of the cross-section L2 and the strain component 11
as extension divided by the original length L. These two ways to find
11 and 11 are obvious, if the cube were homogeneous; the heterogeneity results, however, in non-homogeneous and rapidly oscillating
in the microscale fields of both stress and strain. The latter quantities, calculated from the experiment, thus represent, in fact, averages
of the actual forces and displacements in the cube. More precisely,
they are just the so-called volume (or spatial) averages, to be denoted
by overbar:
11 =

1
V

Z
V

11 (x) dx ,

11 =

1
V

Z
V

11 (x) dx ;

(1.1)

hereafter V = vol V is the volume of the region V.

In turn, the proportionality coefficient
V
E11
= 11 / 11

(1.2)

defines the effective (or gross) Young modulus for the specimen (along
the axis x1 ). This means that the cube, through the relation (1.2),
is homogenized in the sense that its heterogeneity is smoothed out
and it is replaced by a homogeneous one, possessing the gross Young
V.
modulus E11
If we now repeat the same procedure for other cubical samples of
the material, of the same size and orientation as the first, we shall in
V , since the interior phase
general obtain slightly different values of E11
geometry will not be the same in full detail. This explains why the
V.
superscript V explicitly appeared in the notation E11

1.1 Micromechanics of Heterogeneous Media

To define a true material property from such experiments, i.e.

one that is connected with the heterogeneous medium and its internal
structure, independently of the specific choice of the subvolumes, two
natural procedures can be employed.
First, perform a great number, N , of experiments on differently
centered cubes (otherwise identical and identically oriented), and
V 0 , E V 00 , etc., for each one. Then,
measure the appropriate values E11
11
to suppress the specimens dependence, it is natural to define

E11
=


1  V0
V 00
E11 + E11
+ ,
N

(1.3)

which is already a true material property. The right-hand side of

(1.3) is the simplest example of the so-called ensemble averaging.
The meaning is that to obtain information about the expected gross
behavior of a heterogeneous medium, we should deal with average
reaction of a whole ensemble of specimens of identical shape and size,
and apply identical external influence. (In the foregoing reasoning
these are simply the cubes taken from different parts of a big piece
of the medium.) The ensemble averaging is one of the basic notion
in the theory of heterogeneous media of random constitution, see e.g.
the book of Beran [BEi]. An elucidating introduction, together with
basic methods and ideas concerning homogenization for such random
media, can be found as well in the lectures of Willis [WIe].
We shall be concerned in the present survey, however, with more
elementary aspects of the theory, without invoking the proper mathematical framework for random internal constitution. That is why
only the volume averaging (1.1) will be used throughout. It appears
in a natural way, if a second kind of mental experiment is performed.
Take a cube with a fixed center and start increasing its side L (assuming the medium unbounded). Then it is equally natural to expect
that in the limit L/` the sample dependence will disappear:
V

E11
= E11
(1 + O(`/L))

as

`/L 0 ,

(1.4)

, as in (1.3), shows up. Hence,

where the same effective modulus E11
with E V , provided the
we shall be on the safe side, equating E11
11
volume V is big enough compared to the typical inhomogeneity, i.e. if
L  `. But this was a basic assumption, already adopted for the
class of heterogeneous media under study at the beginning of this

Heterogeneous Media: Modelling and Simulation

subsection. The importance of this assumption will be once more

emphasized in Section 1.2.1, when formulating Hashins MMM principle.
Hereafter we shall consider only media for which the ensemble and
volume averages coincide; they are called ergodic.
A satisfactory proof of the statement that a medium is ergodic,
i.e. the two definitions (1.3) and (1.4) yield, in particular, one and the
is not a trivial matter. Such a statement amounts
same value of E11
to the so-called ergodic hypothesis. It can be rigorously justified
under certain additional conditions imposed on the microgeometry
of the media, the most important of which is the macroscopic (or
statistical) homogeneity, see again the book [BEi] for a brief, but
elucidating discussion.
Intuitively, the coincidence between the ensemble and volume av is to be specified by the
erages should not be a surprise, since E11
mediums microgeometry solely. Any natural definition of the homogenization properties, in which sample dependence is somehow
eliminated, should then result in a true and unique material characteristics.
Fully similar considerations can be repeated almost literally, using
different physical backgrounds, chosen and tailored by the personal
taste of the reader. For example, instead of elasticity, diffusion of a
solute through the same cube V can be considered, fixing the solute
concentration at a pair of opposite faces and, in the steady state,
measuring the total amount of the solute, passing through them, see
[HUa]. The dielectric context is also used very often, starting with
the classical works of Mossotti [MOc], [MOd], and Faraday [FAa];
the posing of the homogenization problem in this case is very clearly
explained and discussed, e.g. by Bergman [BEk]. A fully equivalent,
from a mathematical point of view, context concerns heat propagation which will be mainly used in what follows. That is why we should
recall here the well-known equations
q(x) = 0 ,

q(x) = (x)(x) ,

(1.5)

that govern the temperature field (x) in a medium at the absence of

body sources. In Eq. (1.5) q(x) is the flux vector, and (x) represents,
for the heterogeneous body under study, a rapidly fluctuating conductivity coefficient, taking different values when the point x scans

1.1 Micromechanics of Heterogeneous Media

the medium. In the case of a two-phase (binary) material, which very

often appear in various applications,


(x) =

1 , if x constituent 1,
2 , if x constituent 2.

(1.6)

It should be emphasized that all the results, obtained below in the

heat conduction context, are equally well applicable in the contexts of
electrical conductivity, dielectric behavior and magnetic permeability.
A list of analogous quantities are conveniently summarized in a table
by Hashin [HAi, p. 496]. The reason is the obvious mathematical
similarity between between these problems: in all of them a solenoidal
vector field (say, the dielectric polarizability, or the heat flux, etc.)
is proportional to a potential vector field (the electric field, the
temperature gradient, etc.) with a proportionality coefficient (which
may be a tensor) rapidly varying throughout the medium.
The second context to be more extensively used in the present paper is the already employed, for illustrative purposes, (linear) elasticity problem (see Section 1.2.7 below). The basic ideas and conclusions
are the same in both cases, but their implementation in the heat conductivity case is, however, simpler, since the basic quantitiesthe
heat flux and temperature gradientare vectors. In the elasticity
case the respective quantities are tensors (the stress and strain ones)
which makes the calculations more cumbersome technically. That
is why we shall concentrate in what follows primarily on the scalar
conductivity problem; afterward, the elasticity counterparts of the
results will be discussed more concisely.
1.1.3 Brief Historical Remarks
The homogenization quest, due to its utmost importance, has kept
reappearing in various guises and contexts in the last two centuries,
attracting the interest and the efforts of some of the most illustrious
names in science. Here we shall try to collect some historical details, without any attempt or claim to be full and exhaustive.2 One
of our aims is to underline the common micromechanical origin
2

We do not mention even the famous Darcys book [DAa] that initiated tremendous research on flow in porous media, and on Darcys law, in particular. A great
number of recent developments in this field can be found in other papers of the
present volume.

Heterogeneous Media: Modelling and Simulation

and nature of the appropriate problems, despite the widely different

scientific communities and contexts in which they have shown up.
[Link] Navier and Cauchy
It should be recalled above all that the origin and the base of the
modern continuum mechanics, as laid down by Navier and Cauchy,
is to a great extent micromechanical. The starting point of these
authors was the Newtonian picture of matter as an assemblage of
material molecules interacting by means of central forces. In a language, used today, this is nothing but a discrete model and their
main and great contribution in this field is the homogenization of
such a model. More precisely, imaging the molecules as material
points and the interaction forces depending upon their displacement
from the equilibrium position, first Navier and then, almost in the
same time and in a more general form, now classical, Cauchy derived3 the equations describing an elastic continuum, spread over
the discrete modelin the sense that it macroscopically behaves in
the same way as its underlying molecular skeleton. This discrete
model was then put aside for more than a century and the attention,
in continuum mechanics at least, was focused on the appropriate partial differential equations that formed, together with the equations of
hydrodynamics, the core of research in the Mathematical Physics of
the previous century.
In a new guise and on a new level, similar homogenization problems reappeared in continuum mechanics almost a century later when
Einstein [EIa] in 1905 and Bruggeman [BRg] in 1937 considered hosts
of inclusions in a fluid and in an elastic medium, respectively. In both
cases the problems were again micromechanical, in the sense that homogeneous continua, spread over discrete arrays of inclusions, should
be constructed. The interpretation of micro was of course different from the pioneering works of Navier and Cauchy, but the sense
of the homogenization quest was the same. However, the nature
of interaction between the inclusions exhibited a key new feature:
3
According to Love [LOd] (see the historical survey in his introduction), the
Navier Memoire was presented at the French Academy in May 1821, but was
published only in 1827 [NAa]. The Cauchy Memoire was presented at the French
Academy in September 1821, but also appeared in 1827 and 1928 in his Exercises
de mathematique. The Cauchy derivation revisited can be found in Appendix B
of the same Love book [LOd].

1.1 Micromechanics of Heterogeneous Media

its magnitude already was not prescribed, like a simple elastic one,
say, but instead resulted from the local disturbances of the fields like
displacement, velocity, etc., generated by each single particle in the
medium.
[Link] Poisson and Mossotti
With elastic bodies homogenized, the micromechanical studies of
heterogeneous media for almost a century turned toward the gross
behavior of charged materials.4
In 1824 Poisson presented his first Memoire on the theory of magnetism at a meeting of the French Academy [POa]. There he put
the basis of the theory of induced magnetism, assuming a model of
a nonconducting material (we would call it a matrix today) in which
small conducting particles are distributed. For magnetically isotropic
bodies Poisson took the particles spherical, isotropically arranged as
well, and noted that an adequate theory of magnetism can be constructed if one requires that each sphere is lying in a certain constant
external field. In his second Memoire [POb], he treated anisotropic
magnetic solids as well, representing them as an array of thin aligned
three-axial conducting ellipsoids in a matrix. In these two Poissons
papers one can see not only the origin but the rudiments of the basic
models and ideas that prevailed in the theory of heterogeneous media for almost a century after their appearance, and are still useful
and applicable. In particular, the model of the so-called dispersion of
nonoverlapping spheres and ellipsoids seems to have emerged for the
first time in this particular (magnetic) context; the model, in general, turned out very appropriate for many heterogeneous media of
particulate type (matrix-inclusion composites, say), and also in the
liquid-state theory.
It was Faraday in 1838 [FAa], who employed Poissons ideas to
dielectrics, proposing to treat them as an array of metallic spheres,
immersed into an insulating material (matrix) [FAa]. But a more
detailed analysis of this model was due to Mossotti [MOc] who specified the embedded particles as material molecules, modelled them
as little conducting spheres, and described the polarization of a single one. The main contribution of Mossotti was however a method
4
In the electric story we shall utilize a number of historical details, given by
Wiener [WIa], Brown [BRe], Landauer [LAc] and Trimarco [TRa].

10

Heterogeneous Media: Modelling and Simulation

of interconnecting the micro- and macroproperties of the dispersion

which turned out to be the first of the so-called cavity approaches,
and this was the heart of his second paper [MOd]. The basic, and
extremely fruitful, idea of Mossotti was to remove a material region C
around a point in order to evaluate the (electric) force exerted on the
point, due to the interactions with the molecules outside C. In this
way the first of the cavity methods in theory of heterogeneous media emerged. The field acting on the cavity is the local (or effective,
or internal) field that differs from the applied macroscopic one, as a
consequence of the presence of the other molecules. In Section 1.5.3
below, when discussing the so-called effective field approach for evaluating the effective properties, we shall specially outline its connection
with the cavity idea, see Section [Link]. As pointed out by Landauer [LAc], the Mossotti work [MOd] is based, however, on etherial
concepts and notions, typical for the epoch, which makes this classical paper difficult for a modern reader to follow. The idea of local
field was reintroduced and exploited by H. Lorentz around 1880. A
very clear and now classical exposition of this idea and its application
to dielectric and optical properties5 is given in his well-known book
[LOb, Section 117 and App. 54].
Almost at the same time Clausius in his book [CLa, pp. 6777]
revisited Mossottis internal field approach in a concise way, much
easier for a modern reader to follow [LAc]. He cut each molecule
out of the medium and replaced the rest with a material possessing
the effective dielectric constant, reaching the explicit and now classical formula for the latter, given in Eq. (1.19) below. This is, as a
matter of fact, the first appearance of the so-called effective medium
idea that much later was successfully exploited by Bruggeman, Skorohod, Hill, Budiansky et al. in the theory of heterogeneous solids,
see Section 1.5.2 below. As pointed out by Brown [BRe], Mossottis
paper does not contain explicitly the formula for the effective dielectric constant; Clausius book does, and that is why H. Lorentz [LOb,
Section 124] attributes the formula (1.19) to Clausius and Mossotti,
5
Recall that according to Maxwell, 2 = , where is the dielectric constant
of a transparent medium and = c/v is its refractivity index (i.e. the ratio c/v
between the velocity of the light in the vacuum and in the medium). This classical formula of the electromagnetic theory explains the appearance of equivalent
formulae for the effective properties, containing refractivity instead of dielectric
properties.

1.1 Micromechanics of Heterogeneous Media

11

and it has come to be known as the Clausius-Mossottis one. But

several years before Clausius publication, an analogous formula in
conductivity had been given by Maxwell [MAj], whose reasoning will
be presented below.
Note also that L. Lorenz in three papers between 1869 and 1880
considered optical problems, concerned with propagation of electromagnetic waves (optical vibrations as he called them) [LOc]. He
assigned a refractive index to each molecule that differs from that
of the surrounding medium and then employed a scheme which resembled very much what is now called the Coherent Potential Approximation. That is why the local field and its expression, to be
discussed in an appropriate context in Section [Link], is often called
Lorenz-Lorentzs field.
[Link] Maxwell
In his famous treatise [MAj] Maxwell addressed in passing the problem of predicting the overall behavior of a dispersion of spheres, immersed into a matrix of different conductivity. We shall repeat here
his elegant and simple reasoning, since it contains a number of basic
ideas and drawbacks that have kept reappearing in the later works
for a long period afterward.
The basic building block in Maxwells study is the single sphere
field. Let a single spherical inhomogeneity Va , located at the origin,
be immersed into an unbounded matrix. Let the applied far-field
have a constant and prescribed gradient, i.e. (x) = G x. The
temperature field, (x), is governed then by the equation:
((x)(x)) = 0 ,

(x) G x as |x| ,
(1.7)
[] = 2 1 ,

(x) = 1 + []ha (x) ,

where ha (x) is the characteristic function of the sphere Va , 1 and 2

are the conductivities of the matrix and of the sphere, respectively,
see Eqs. (1.5) and (1.6).
Since (x) is a step-constant function, the field (x) should be a
harmonic function both inside and outside the sphere; on the sphere
surface, |x| = a, it should be continuous together with the normal
component of the heat flux:




2
= 1
,

r r=a0
r  r=a+0

(1.8)

Heterogeneous Media: Modelling and Simulation

12

as it follows directly from (1.7); here r = |x| is the radial coordinate,

a is the spheres radius.
The solution of Eq. (1.7) depends linearly on the far-field gradient
G. Together with the obvious geometrical symmetry, this suggests
to look for the needed solution in the form
(x) = G (r) .

(1.9)

Since (x) is harmonic, Eq. (1.9) implies

(x) = G (r) = 0 ,

i.e.

(r) = const ,

and hence
(r) = 12 C 0 r2 C 00 a/r + const .

(1.10)

Obviously, the additive constant in the last formula does not influence
the solution (1.9). Inside the sphere C 00 = 0, in order to have a
bounded everywhere solution there, so that
(x) = C 0 G x ,

|x| < a .

(1.11)

Outside the sphere C 0 = 1, in order to meet the far-field condition in

Eq. (1.7), and hence


(x) = G x 1 C 00 (a/r)3 ,

|x| > a .

(1.12)

The continuity conditions for (x) and for the normal flux at |x| =
a, see Eq. (1.8), yield
C 0 = 1 + C 00 ,

2 C 0 = 1 (1 2C 00 ) ,

as it follows from Eqs. (1.11) and (1.12), and also from the simple
formula

= er (x) = G (r) er = 00 (r)G er ,

(1.13)
r
er = x/r, r = |x|, which holds for any function of the form (1.9).
The solution of the single inclusion problem (1.7) therefore reads
(

(x) = G x

if |x| < a ,

,
(a/r)3 ,

if |x| > a ,

=1 =

31
.
2 + 21

(1.14)

with the constants

=

[]
,
2 + 21

(1.15)

1.1 Micromechanics of Heterogeneous Media

13

Remark 1.1 In Section 1.4.2 we shall find the solution of Eq. (1.7)
inside the sphere in a much simpler way, which works equally well
for an ellipsoidal inhomogeneity also. In Maxwells considerations
one needs, however, the field (x) outside the sphere, so that a more
lengthy analysis is required.
Remark 1.2 The formula (1.14) shows that if the temperature gradient at infinity is constant, the same holds true for this gradient
inside a spherical inhomogeneity. This is a particular case of a more
general statement in which the sphere is replaced by an ellipsoidal inhomogeneity. A similar statement holds for the elastic inhomogeneity
problem, see Section 1.4.4.
With the single-sphere solution found, we can go back to Maxwells
analysis. Consider a large sphere VA of the radius A and conductivity 1 , containing N identical and nonoverlapping small spheres of
conductivity 2 . (Large means here that A  a, where a is the
radius of the small spheres.)
Immerse the sphere VA into a homogeneous medium possessing the
same conductivity 1 as that of the matrix. Sit down at the point
M , far away from the center O of the sphere, i.e. imagine R  A,
where R = |OM |. Assume again that the temperature gradient G at
infinity is constant. The presence of the inhomogeneous sphere VA
will disturb the field Gx that would have existed in the homogeneous
medium of conductivity 1 at its absence. We shall evaluate this
disturbance in two different ways.
First, let us consider each small sphere within VA as a single and
pretend there are no other small spheres around it. Then, according to Eq. (1.14), the disturbance, caused by any such sphere, is
a3 /R3 , and summing up all of them, gives
N

a3
a3 2 1
=
N
;
R3
R3 2 + 21

(1.16)

recall that the small spheres within VA are located practically at one
and the same distance R from our observation point M , since R  A.
Second, looking at the big sphere VA from a great distance, we
can homogenize it, neglecting its heterogeneous internal structure
and replacing it by a homogeneous sphere possessing a certain unknown effective conductivity . The disturbance, caused by such a

14

Heterogeneous Media: Modelling and Simulation

homogeneous sphere at the point M , is

A3
A3 1
=

,
R3
R3 + 21

(1.17)

as it follows again from Eq. (1.14).

It is natural now to equate the disturbances (1.16) and (1.17),
calculated in these two different ways (which just reflects the above
mentioned idea of homogenizationreplacement of the heterogeneous sphere VA by an equivalent homogeneous one; equivalent in the
sense that it behaves (i.e. disturbs) the temperature field in exactly
the same way). This results in the equation
N
or

a3 2 1
A3 1
=

,
R3 2 + 21
R3 + 21

1
= 2 ,
+ 21

2 = N

Va
a3
=N 3,
VA
A

(1.18)

so that 2 is the volume fraction of the spheres in VA ; Va = 43 a3 and

VA = 34 A3 are the volumes of the respective spheres. From (1.18) it
follows

1 + 22
=
,
(1.19)
1
1 2
and this is exactly the formula, known now as Clausius-Mossottis (in
dielectric context) or Maxwells (in conductivity context), and also
Lorenz-Lorentzs (in refractivity context, see the footnote on page 10).
Having derived (1.19), Maxwell remarked that it may be obtained in
other ways, but that here given involves only repetition of the result
already obtained for a single sphere [MAj, Section 314].
Due to the above explained reasons, and to the fact that we shall
deal mainly with conductivity problem in the sequel, we shall refer
to (1.19) as Maxwells formula. Its derivation possesses the obvious
drawback that each sphere is considered as single, so that the result
(1.19) is strictly valid only in the dilute limit 2  1:

= 1 + 32 + O(2 ) ,
1

(1.20)

when the spheres are sparse. This fact was entirely clear to Maxwell
himself (see the same Section 314 of his Treatise) and to his contemporaries. Lord Rayleigh [RAa], in order to check the exactness and

1.1 Micromechanics of Heterogeneous Media

15

applicability of (1.19) for higher values of 2 , considered a special kind

of internal constitution, assuming that the inclusions are spheres or
cylinders, forming a regular periodic array. (A vast number of studies concerning the properties of periodic media were initiated by this
classical paper.) The basic Rayleighs conclusion was that (1.19) does
hold for dilute to moderate concentrations 2 , but for closer packing
of spheres and cylinders, considerable deviations from the predicted
values are present.
The evaluation of the terms of higher orders in Eq. (1.20), e.g. the
one proportional to 22 , requires already that spheres interactions be
taken into account. If it is done properly, then the spatial statistics
of the sphere distribution should enter the result. The appropriate
problem appears very hard and had to wait exactly 100 years after
Maxwell, when Jeffrey [JEa], in 1973, derived the first exact formula
for the successive 22 -coefficient in the power expansion (1.20).
[Link] Einstein
The dawn of the 20th century witnessed the emergence of some new
micromechanical problems whose motivation was already different
from the mainly electrical preoccupation of the previous century
(a notable exception was the classical Voigt quest for the effective
elastic properties of a polycrystal [VOa]). In a pioneering paper on
viscosity of dilute suspensions Einstein [EIa], in his PhD dissertation,
addressed the question of experimental measurement of the size of
molecules. His simple and innovative reasoning may be summarized
as follows.
Consider a known amount of molecules, n, per unit volume of a
solvent. The molecules are identical, and each one is imagined as
a rigid sphere of an unknown radius a. Thus there exists a volume
concentration 2 = nVa = n 43 a3 of particles in the solvent, i.e. a
liquid suspension is to be studied. Let the liquid be Newtonian,
of viscosity 1 . The presence of the rigid spheres will undoubtedly
make the suspension more viscous, so that it will possess a certain gross macroscopical viscosity . The micromechanical problem
consists in evaluating , provided 1 and the volume fraction 2
of the molecules are known. Einstein solved this problem for slow
motions (in the Stokesian approximation, that is) and in the dilute
case, when 2  1 and hence the many-particle interactions can be

16

Heterogeneous Media: Modelling and Simulation

neglected. Each sphere (molecule) can be then considered as single,

immersed into a given (shear, to be more specific) flow of the fluid at
infinity. The presence of the sphere, perturbing the flow, will increase
the energy dissipation, by a certain quantity , proportional to the
spheres volume Va . Technically, the correct formula for was the
core of Einsteins calculations. Since the interactions are neglected,
the total dissipations increase will be n, i.e. n times higher. This
increase, from a macroscopical point of view, can be attributed to the
change of the gross viscosity of the suspension, which has the needed
value , and this fact led him to the now famous formula

= 1 + 2.52 + O(2 ) ,
1

(1.21)

named after him.6 Though Einstein reasoning is based on some simple solutions of the Stokes equations, his derivation is quite technical
and will be omitted here. The classical formula (1.21) will be derived
in Section 1.5.2 below as a particular case of a more general result
(see Remark 5.6).
The term O(2 ), added in the right-hand side of Eq. (1.21), underlines the fact that it is strictly valid, similarly to the Maxwells one
(1.19), only in the dilute case 2  1, when the spheres interactions
can be totally neglected.
Applying the formula Eq. (1.21) for sugar solutions and using
some viscosimetric data, Einstein was able to obtain experimentally
the volume concentration 2 = nVa of the sugar molecules; for a
given n, it led him to a certain value of their radius a which turned
out to represent reasonable approximation of the latter. For our purposes, however, much more interesting is the micromechanical problem, posed by Einstein, namely, how to evaluate a macroscopic property (viscosity), using the available microscopic information, i.e., in
the dilute case, the viscosity of the solvent and the volume fraction
of the particles.
[Link] Smoluchowski
Another classical micromechanical problem was introduced in 1916
by Smoluchowski in his Three papers on diffusion, Brownian motion
6
As a matter of fact, Einstein introduced and used the so-called energy
definition of the effective properties; see Section 1.2.4 below.

1.1 Micromechanics of Heterogeneous Media

17

and coagulation in colloid solutions [SMa]. There he considered, in

particular, a dilute array of fixed spheres that absorb a species diffusing around. The question was what is the mean rate of absorption in
such a system of spheres, provided the concentration is kept constant
at infinity.
Since the array is dilute, each sphere (of the radius a) can be
treated as alone, located at the origin. In the steady-state limit the
species concentration c(x) is governed by the equation:
c(x) = 0 ,

c(x) 

r=a

= 0,

c(x) c0 ,

(1.22)

r a.

(1.23)

whose solution is obvious



c(r) = c0

a
1
r

The flux of the diffusing species within the sphere is

Z

j(x) er dS ,

r=a

where j(x) = Dc(x) is the local flux vector (according to Ficks

law) and D is the diffusing coefficient. Employing Eq. (1.23) gives the
value 4aDc0 for this flux. If n is the number density of the spheres,
the total flux K due to all spheres, i.e. the total rate of absorption,
is n times bigger (recall the dilute assumption). Hence
K = ks2 c0 ,

ks2 = 4anD =

3
2 D ;
a2

(1.24)

similarly to the foregoing analysis of fluid dispersions, 2 = nVa is

the volume fraction of the spheres.
Thus there is a balance between creating diffusing species (in order
to keep its concentration, c0 , at infinity fixed) and its removal due
to the spherical traps or sinks. The proportionality coefficient
between c0 and this removal rate, K, is just the so-called effective
reaction-rate constant or the effective absorption coefficient, denoted
by k 2 :
K = k 2 c0 , k 2 = ks2 + O(2 ) .
(1.25)
(k 2 is also called trapping constant or sink strength.)
Similarly to the foregoing Einsteins formula (1.21), O(2 ) is added
in the right-hand side of Eq. (1.25) to underline the fact that the

18

Heterogeneous Media: Modelling and Simulation

Smoluchowskis value ks2 is the dilute approximation of the effective

sink strength, strictly valid only when the spheres interactions are
totally neglected.
Smoluchowskis paper initiated a broad and extensive research
activity, concerned with the so-called diffusion-controlled reactions
chemical events whose rate decisively depends on the transport processes, like diffusion, in a solution of reactant entities. For more details and different points of views and interpretations, concerning the
ideas, results and basic references in this important field the reader
is referred to the papers [FEa], [TAa], or to the survey of Calef and
Deutch [CAa]. For our aims here much more interesting, once again,
is the micromechanical problem, posed by Smoluchowski, namely,
how to evaluate a macroscopic property (the effective rate constant
now), using the available microscopic information, i.e. in the dilute
case, the volume fraction of the spheres.
1.1.4 Outline of the Survey
A number of the simplest ideas, methods and results, concerning homogenization procedures for heterogeneous media, are collected in the
present survey, with no pretense to give a full exposition. The aim
is much more modestto provide a certain basic reading and a general picture for a newcomer to gain some necessary knowledge, before
proceeding to more specific aspects and (extremely voluminous) specialized literature. Another aim lies also in underlying the common
nature and unity of the problems and approaches, spread through
many different communities, to mention only pure mathematics (homogenization of differential equations with periodic and random coefficients), solid mechanics (composite materials), hydromechanics
(liquid dispersion and emulsions, flows through porous media), geophysics (micromechanics of rocks), chemical physics (diffusion-absorption reactions), and so on and so forth. Though the context, the
terminology, the way of thinking and language may widely vary, behind many of the central, for these fields, problems lies one and the
same fundamental quest about the micro-macro interconnection.
The list of references, though extensive, by no means pretends to
be exhaustive. It includes some (but not all) of the basic works in the
fields and may provide a certain guide where to find further details
and references for interested readers. We have tried to mention as
well some of the earlier and perhaps not very well-known works in

1.1 Micromechanics of Heterogeneous Media

19

this field to demonstrate that the basic and classical ideas are old
and, typically for any classic, have been rediscovered many a time.
Only static behavior and static approaches will be discussed here.
A thorough and critical discussion of somewhat similar methods, using wave phenomena and, in particular, features of scattered fields by
single inclusions in a particulate medium, can be found in Berrymans
papers [BEm], [BEn].
The volume-average approach is used throughout, which may explain the word elementary in the title. Though this approach possesses certain inherent inconsistencies at a closer look (underlined in
due course below), it seems more appealing to the author for introductory purposes and for a broader audience. The word elementary underlines also the fact that in recent decades heterogeneous media have
been studied by means of new and much more refined methods. The
statistical modelling has already been mentioned, together with the
book [BEi] (see also [KRf]). The methods, based on various homogenization procedures (such as two-scale expansions, G-convergence and
so on), have received special attention; see, e.g. the books [SAa] or
[ZHa]. The modelling aspects, concerned with determination of optimal micro-structures, should be specially mentioned as well. The
appropriate problems, ideas, and results in this field can be found in
[MIg]. The author hence believes that in view of the recent extensive
research activity, combining new mathematical tools and approaches,
the reasoning and the results, collected in the present survey, indeed
deserve to be called elementary.
The basic notions and definitions are introduced in Section 1.2;
there, as in the rest of the paper, the scalar conductivity context is
first used for illustration. The reason is that, being simpler formally,
it allows one to concentrate on the basic problems, unshaded by technical details that may appear in a similar elastic context. There are
two central results of the Section: (i) the equivalence of the energy
definition of the effective properties to the usual ones, under homogeneous boundary conditions; (ii) the Willis result of consistency,
stating that the effective conductivity and resistivity are the inverse
of each other. Moreover, the Willis reasoning, reproduced in a simplified situation here, clarifies the special role of the homogeneous
boundary conditions and proclaims their extremum properties.
Some of the evergreens of the theory of heterogeneous media

20

Heterogeneous Media: Modelling and Simulation

are collected in Section 1.3. This is the perturbation theory for

a weakly inhomogeneous materials, as initiated by Brown in 1955
[BRd]. The Bergman formula is then derived [BEk], which interconnects the derivatives of the effective conductivity with respect to the
properties of the constituents with the mean energies, accumulated
within the latter. Only one of the simpler applications of this important formula is exploitedthe beautiful Matheron derivation [MAi] of
the Hashin-Shtrikman estimates on the effective conductivity, without invoking the original variational arguments of Hashin and Shtrikman [HAk]. These bounds hold under the sole assumption that the
material is macroscopically isotropic. The Hashin assemblage [HAe]
is then constructed to demonstrate that the bounds are sharp, in the
sense that they cannot be improved. The Berans bounds are then
derived and simplified, after Torquato and Milton; the appearance
of three-point statistical parameters there indicates that the effective
properties depend in general on all details of the internal structure,
as reflected by the appropriate multipoint characteristics (correlations) of a heterogeneous medium. The cross-property relations and
their philosophy are briefly discussed as well; the original derivation of the first such relation, due to Levin [LEa], is supplied: it
interconnects exactly the effective bulk modulus and thermal expansion coefficient.
In Section 1.4 the single particle problem is discussed in various
contexts. The well-known Eshelby result [ESa] is first proved, stating
that the temperature gradient (or the strain) within an ellipsoidal inhomogeneity immersed in a homogeneous matrix is constant provided
the applied gradient (or strain) at infinity is constant. The reasoning here is simplified, as compared to the original Eshelbys one, and
follows, as a matter of fact, Wu [WUa].
The aim of Section 1.5 is to show how various repetitions (to cite
again Maxwell [MAj, Section 314]) of the single-inclusion solution of
Section 1.4 can be systematically utilized in order to derive the basic
and most widely used approximate (static) models for predicting the
effective properties of particulate media, like self-consistent, differential scheme and effective field. The starting point are the concentration factors, introduced in Section 1.2, and the representations, due
to Reynolds and Hough [REb] and Hill [HIb] of the effective properties, through these factors. Each approximation then corresponds to
a certain qualitative and intuitively appealing replacement of a con-

1.2 Micromechanics of Heterogeneous Media

21

centration factor by a quantity that stems out from the appropriate

solution for a single inhomogeneity.
Section 1.6 demonstrates, after Hershey [HEb] and Kroner [KRd],
the application of the self-consistent arguments for predicting the effective elastic moduli of polycrystals. As an illustration, the simplest
case of cubical symmetry of grains is treated at some length.
The survey is a (considerably) extended version of the course on
composite and heterogeneous media, read by the author in the last
years for graduate students in the Faculty of Mathematics and Informatics at the St. Kliment Ohridski University of Sofia. (A number
of topics were included also in the authors lectures presented in 1994
at the Department of Engineering Sciences of the Istanbul Technical
University and, in 1997, in the Department of Mathematics of Torino
Polytechnics.) As a consequence, the exposition is quite detailed with
no incantations of the kind it is easily seen . . . . Besides a natural curiosity, nothing more than a basic knowledge of differential
equations, continuum mechanics and tensor calculus is required from
a reader.
As a rule, the bold-face lower-case letters e, x, y, etc., are used
for vectors. The upper-case such letters A, K, etc., are reserved for
second-rank tensors, the shadow symbols A, L, etc., denote fourthrank tensors. Dyadic notations are used, e.g. xy = xy is the dyadic
product of the vectors x and y. The dot stands for a contraction
with respect to one pair of indices, and the colonwith respect to
two pairs, say,
(A B)ij = Ai Bj ,

(A : B)ijkl = Aij Bkl ,

in a Cartesian system, with the summation convention adopted

throughout.

1.2 The Homogenization Problem

The simplest point of view, already introduced and briefly discussed
on the examples in the previous Section, is that a heterogeneous
medium behaves macroscopically in the same way, as do its constituents, but with different, effective, values of the appropriate material constants. In the present Section we shall try to put this point

22

Heterogeneous Media: Modelling and Simulation

of view on a clearer framework, introducing and analyzing the definition of the effective properties.
1.2.1 The Representative Volume Element
A basic notion in micromechanics is the representative volume element (RVE). Though it can be more rigorously defined, see e.g. [DRa],
it is appropriate to be imagined in the following manner [HIb]. This
is a volume which is small enough from a macroscopical point of view
and could be thus treated as a typical point of the heterogeneous
continuum under study. On the other hand, it should be large enough
in the microscopical scale, in order to contain a large number of single
inhomogeneities and therefore to be indeed typical representative
for the microstructure of the solid. A more detailed discussion of
RVEs, together with certain criteria how to identify them, can be
found in the book [NEa, Chapters 1.3 and 2.5.4].
It is noted that the notion of RVE (as well as the volume averaging
procedure to be used in the sequel) was clearly and vividly described
already in the classical Lorentz book [LOb, Section 113]. The author
introduces there the mean value , over a ball S, of a microscopically
fluctuating quantity , and adds that the size of the ball should be
neither too big nor too small, in order to get a meaningful macroscopic
characterization of . Lorentzs comments on the words too big and
too small, that followed in the same Chapter, is just what is now
understood as a RVE.
The foregoing definition of a RVE will undoubtedly leave unhappy any reader with a more rigorous taste. It should be pointed
out, however, that this notion is very convenient on a heuristic and
elementary level of consideration, pursued in this survey, and allows
us to introduce and develop the basic ideas and some of the simplest techniques of micromechanics, in general, and of mechanics of
composites, in particular. A more rigorous approach, say, for random
media, does not need the notion of RVE at all. Note that for periodic
media the RVE is unambiguously defined (this is simply the typical
unit cell of the structure). It is to be noted also that our RVE here
is just the well-known material point of a continuous mediumthe
basic notion, which lies in the foundation of continuum mechanics.
The only difference is the size of inhomogeneities: for the latter they
often correspond to molecular or atomistic scales. For heterogeneous
solids this size is defined by the dimensions of a typical inclusion or

1.2 Micromechanics of Heterogeneous Media

23

of a monocrystal (in a polycrystalline aggregate).

The heterogeneous media considered here can be thus characterized, after Hashin [HAi], in the following manner. It is possible to
distinguish within them representative volume elements which are
small macroscopically. Their typical dimension defines the so-called
mini-scale. A RVE consists of a big amounts of much smaller inhomogeneities (inclusions, monocrystals, etc.), whose dimension defines
the micro-scale. The composite medium itself consists of many RVE;
its size defines, in turn, the macro-scale. For the very idea of homogenization to be sensible and to produce widely applicable results, the
following inequalities between the above defined three length scales
should hold:
MICRO  MINI  MACRO .
(2.1)
After Hashin [HAi], we shall call (2.1) the MMM principle.
Note that though general, the violation of the MMM principle in
certain approximate homogenization schemes does not imply that the
resulting predictions are useless or should be immediately rejected.
It rather indicates that the internal structures, for which the schemes
may be applicable and useful, have little in common with the simplified models (like dispersions of identical particles), utilized when
deriving the appropriate approximations for the effective properties
(see Section [Link] below).
Remark 2.1 The MMM principle can be traced back as far as the
Poissons 1829 Memoire [POc, p. 149], where he wrote: The molecules are so small and so close to one other that a portion of a body
containing an extremely large number of them can be still supposed to
be extremely small and to have an insensibly small volume. (translation and citation after Arnold [ARb, p. 364]). As pointed out by
Arnold in the same paper, this statement is crucial to Poissons conception of matter. In our context, the latter is but a heterogeneous
medium (with the molecules treated as certain inhomogeneities),
which we macroscopically observe and study as homogeneous on this
level.
Remark 2.2 The fact that only two length scalesmicro and mini
are considered hereafter, besides the macroscopic one, means that one
more assumption is tacitly adopted: namely, that the applied external field is either macroscopically homogeneous or varies slowly. If

24

Heterogeneous Media: Modelling and Simulation

this field varies rapidly, as it is the case with the wave propagation
problems, when the wave-length is of the order of the micro-length
`, the situation drastically changes. The homogenization can still be
performed (though ensemble averaging should be used from the very
beginning), but the macroscopic equations will be nonlocal. Any
discussion of nonlocality goes, however, far beyond the scope of the
present survey. It suffices to say only that in this case the stress,
e.g. in a point would depend on the strain within the entire body
through a certain integral operator. For details, discussion and references, the reader is referred to the book [KUa], see also [WIe].
The position of a typical mini-element (a RVE) is identified by
X, so that all continuum quantities like temperature, flux, stress and
strain fields, etc., are functions of X. (Hereafter, stationary problems will be only discussed, and hence the time t is not present.) To
distinguish these fields from the same ones within the RVE which
exemplifies the internal structure of a typical continuums point,
the continuum fields are referred to as macro-fields, and those within
a RVE as micro-fields, respectively. That is, instead of saying, for example, the temperature within the continuum or the temperature
field within the RVE, that corresponds to the particle X, the expressions macrotemperature field and microtemperature field are used. In
a similar manner, the continuum displacement, mass-density, strain,
stress, and other physical quantities are identified by an appropriate
use of the prefix macro, and those within a RVE by the prefix micro. Again the analogy with the typical reasoning, contained at the
opening pages of textbooks on continuum mechanics, can be easily
observed: The only difference lies in the interpretation and in the size
of the material pointthe RVE in our case.
The physical quantities of interest on the micro-level, i.e. within
the RVE, depend also on the local coordinate x, e.g. the temperature
= (X; x) varies both on macro- and micro-levels. (The dependence of both macro- and micro-coordinates X and x indicates that
the internal fields in the medium vary, in general, in different ways
in different RVEs.) The connection between the macro- and microquantities is supplied by the volume averaging, with respect to the
micro-coordinate x, e.g.
1
(X) =
V

(X; x) dx ,
V

etc. ,

(2.2)

1.2 Micromechanics of Heterogeneous Media

25

where V is the RVE, attached to the macro-point X, and V =

vol V.
Since the point X spans the body, the quantities like (X) in
Eq. (2.2) are called moving averages. Such averages play a central
role in the elementary theory of effective properties, developed and
discussed below.
In what follows the analysis will be restricted to the important
class of statistically homogeneous materials. In the terminology of
moving averages it means, roughly speaking, that the macroscopic
properties of almost all RVEs (more precisely, except those near the
macroscopic boundary of the solid) are one and the same, see [HAf]
and [NEa] for more details and comments.
The assumption of statistical homogeneity allows us to deal in
what follows with a single RVE, V. The latter indeed will be then
representative, in the sense that it will possess the same gross properties as all of the rest of such elements. These properties will then
coincide with those of the medium as a whole. That is why the
words RVE, the volume V and the medium will be used hereafter
interchangeably as synonyms.
In most of the cases, treated below, we shall assume that the
medium is statistically isotropic as well; again, without entering a
longer discussion, we shall only point out that the latter means that
the macroscopic properties under study are independent of direction.
1.2.2 Definition of the Effective Conductivity
Let V be a RVE. Consider the heat propagation problem (1.5) in V,
with the boundary condition



(x) 

= G x,

(2.3)

V is the boundary of V. If the RVE is homogeneous, then (x)

G x, x V, and thus G. That is why (2.3) will be called a
homogeneous boundary condition.
For a heterogeneous RVE, the field (x) fluctuates inside V. However, the boundary conditions (2.3) yield that
= G ,

(2.4)

as if the volume were homogeneous. This fact is a consequence of a

more general statement, see [HAh, Section 3.1] or [NEa], which we
shall first formulate in the simpler scalar case:

Heterogeneous Media: Modelling and Simulation

26

Lemma 2.1 Whatever the composition and microstructure of the

RVE, the volume average of the temperature gradient is completely
determined by the temperature of the boundary V.
Proof. Let (x) be a temperature field in V, attaining prescribed
boundary values 0 (x). Using the Gauss theorem, we have
1
=
V

1
=
V

(x) dx

1
n(x) dS =
V
V

(2.5)
Z

n (x) dS ,
V

where n is the outward unit vector to the boundary V. Eq. (2.5)

shows that the boundary values 0 (x) do specify the volume average
uniquely.
In particular, if 0 (x) = G x, Eq. (2.5) implies
1
V

1
=
V

nG x dS

1
(G x) dx = G
V
V

(2.6)
Z

x dx = G ,

which proves (2.4) (having used once more the Gauss theorem and
the fact that x = I, where I is the unit second-rank tensor).
Consider now Eq. (1.5) with respect to the temperature field, i.e.
((x)(x)) = 0 ,

(2.7)

with the boundary condition (2.3). Since (x) > 0, it can be easily
shown that the BVP problem (2.7), (2.3) possesses a solution, (x),
which is unique. The latter allows us to evaluate the mean heat flux
over the RVE under study
Q = q = (x)(x) =

1
V

(x)(x) dx .

(2.8)

Due to the obvious linearity of the boundary-value problem (2.3),

(2.7) with respect to G, the vector Q is a linear function of G, i.e.
Q = K G ,

(2.9)

1.2 Micromechanics of Heterogeneous Media

27

where K is the second-rank tensor of effective conductivity of V. In

the statistically isotropic case, this tensor is spherical, so that (2.9)
reduces to
Q = G ,

(2.10)

with being the effective (or overall, or macroscopic) conductivity

of the volume V.
Eq. (2.9) means that we replace the micro-inhomogeneous volume
V by a homogeneous one, with conductivity K , which, from a microscopical point of view, reacts in the same way. Thus, if we are not
interested in detailed microscopic fields within the RVE, but rather in
the macroscopically measured quantities like the volume averages q,
, etc., we can homogenize V, i.e. treat it as a homogeneous entity
possessing the overall conductivity . The problem then consists
in evaluating K making use of appropriate information about the
microstructure of V. This is a typical and very difficult problem of
micromechanics of heterogeneous and composite media.
The fundamental reason that makes the evaluation of the effective properties so difficult, can be well seen from the very structure
of Eq. (2.7). At first glance, the BVP (2.7), (2.3) looks innocuous,
resembling the ones treated in the basic books on PDEs. In these
books, however, the emphasis lies on the dependence of the solution upon the boundary data, which is obviously linear. Here the
boundary condition is fixed, and one should find the solution (x)
for given coefficients (x); the dependence of (x) upon the latter
is already nonlinear, as first pointed out by Kraichnan [KRa] (in a
physically entirely different context though). This gives rise to very
serious mathematical complications, common for all problems concerning heterogeneous media and their homogenization. Typical and
very important among them is the fact that even the mean value of
the solution, like the one in Eq. (2.8), depends on all details of the
internal structures of the medium. More persuasive and specific arguments in support of this statement will be given in Section [Link]
below, on the example of a weakly inhomogeneous medium.
Remark 2.3 If, by chance, the temperature gradient is constant
within V, i.e. (x) G x, if x V, the effective conductivity
is exactly the mean of the constituents conductivities. Hence

28

Heterogeneous Media: Modelling and Simulation

for an N -phase medium

= ,

N
X

i i ,

(2.11)

i=1

as easily seen from Eqs. (2.8) and (2.10); here i = Vi /V is the volume
fraction of the ith constituent that occupies the volume Vi V; i is
its conductivity, i = 1, . . . , N .
The assumption that the temperature possesses a constant gradient within the medium was adopted, as a matter of fact, by Voigt
[VOa], in his study of elastic moduli of polycrystals (accordingly,
he proposed that the strain field is homogeneous within the polycrystalline sample). That is why the approximation (2.11) is called
Voigts and is designated by the superscript v, i.e. v = .
It is noted that a temperature field with a constant gradient does
appear in a heterogeneous solid, provided it is a fiber-reinforced or
layered material, and the applied macroscopic gradient G is along the
fibers or layers, see Section 1.5.2. If however the gradient G is perpendicular to the laminae in the layered material, though the volume
fractions of the constituents remain the same, its effective conductivity is already totally different. (In this case it coincides with
the Reuss value, see Eq. (2.18) and Remark 2.2 below.) This simple
observation demonstrates that is indeed sensitive to the internal
structure of the medium. As a consequence any theory that tries to
predict the effective behavior on the base of the volume concentrations and phases properties solely can serve, at the best, only as a
certain heuristic approximation.
1.2.3 Definition of the Effective Resistivity
Recall that the conduction through a solid can be alternatively studied, choosing as a basic variable the heat flux q(x) instead of temperature (this is the so-called dual formulation). Since q(x) is divergencefree, q(x) = 0, it can be represented by means of its vectorpotential (x):
q(x) = (x) ,
(2.12)
where () = curl () is the curl of the respective vector field.
In turn, the curl of any potential field and, in particular, that of
(x) = k(x)q(x), vanishes. Here k(x) = 1/(x) is the so-called

1.2 Micromechanics of Heterogeneous Media

29

resistivity or (inspired by the elasticity terminology), the compliance

field of the medium.
Hence, the vector potential (x) satisfies the equation
(k(x) (x)) = 0

(2.13)

which is the counterpart of the basic Eq. (2.7) for the temperature
field.
Let us prescribe the heat flux on the boundary V



qn 

= Q n,

qn = q n ,

(2.14)

where Q is a given vector. For a homogeneous V, q(x) Q, for

all x V, and thus q = Q. That is why the boundary condition
(2.14) will be also called homogeneous (with respect to flux); it is the
counterpart of the homogeneous condition (2.3) for the temperature.
For a heterogeneous V the field q(x) fluctuates. However, the
boundary condition (2.14) implies that
q = Q,

(2.15)

as if the volume V were homogeneous, similarly to Eq. (2.4). This

fact is again a consequence of a more general statement, dual to
Lemma 2.1, see [HAh, Section 3.1] or [NEa], which we shall formulate
once again in the simpler scalar case first:
Lemma 2.2 Whatever the composition and microstructure of RVE,
the volume average of the heat flux, at the absence of body sources,
is completely determined by its boundary values on V.
Proof. Let q(x) be a heat flux field in V, attaining the prescribed
boundary values. Due to the assumed absence of body sources
q(x) = 0. In the Cartesian system xi , i = 1, 2, 3, we have, using
the Gauss theorem,
1
V

qi =

1
V

qi (x) dx =

1
V

Z
V

ij qj (x) dx

1
V

1
xi qj,j (x) dx =
V
V

xi,j qj (x) dx =

(xi qj (x)),j dx
xi qn (x) dS ,

(2.16)

Heterogeneous Media: Modelling and Simulation

30

and hence q is indeed uniquely specified by the values qn | V .

In particular, let qn = Q n on V. To check the validity of (2.15),
we repeat almost literally the reasoning after Lemma 2.1, see (2.6).
Indeed, in virtue of Eq. (2.16),
q=

1
V

xqn dS =
V

1
V

xn Q dS =

1
Q
V

x dx = Q ,

which proves (2.15), since x = I.

Consider next Eq. (2.13) for the potential (x), together with
the boundary conditions (2.14). This BVP has a solution which is
unique in the class of solenoidal fields ( (x) = 0). Thus the field
(x) is uniquely defined, and we can evaluate the respective
mean temperature gradient
G = = k(x) (x) = k(x)q(x) = k Q .

(2.17)

The parameter k that enters (2.17) is called the effective (or

overall, or macroscopic) resistivity of the volume V. (Recall that we
have assumed the medium isotropic. Otherwise the constant k in
(2.17) should be replaced by the appropriate second-rank tensor of
effective resistivity.)
The interpretation of Eq. (2.17) is fully similar to that of Eq. (2.10);
it means that we homogenize the RVE, V, under study, replacing it
by a homogeneous one, with the resistivity k . The volume V, from
a macroscopical point of view, will behave then in the same way as
the heterogeneous one.
Remark 2.4 If, by chance, the heat flux is constant within V, i.e.
q(x) Q, x V, the effective resistivity is exactly the mean value
k of the constituents resistivities
k =

1
= k,

k=

N
X
i=1

i ki =

N
X
i
i=1

(2.18)

as easily seen from Eq. (2.17). Here, similarly to Remark 2.1, i and
i are the volume fraction and conductivity of the ith constituent,
i = 1, . . . , N .
The assumption of a constant flux is to be attributed to Reuss
[REa] who, similarly to Voigt, studied the elastic moduli of polycrystals (accordingly, he proposed that the stress field is homogeneous

1.2 Micromechanics of Heterogeneous Media

31

within the polycrystalline sample). That is why the approximation

(2.18) is called hereafter Reuss and is designated by the superscript
r, i.e. r = (1/)1 .
It is noted that such a constant heat flux field does appear in a
heterogeneous solid, provided it represents a layered material, and
the applied macroscopic gradient Q is perpendicular to the layers,
cf. Section 1.5.2.
We shall conclude this section with a comment concerning the
dual formulation (2.13) of the heat conduction problem. Comparing
(2.13) with (2.7), one notices that (x) is the exact counterpart of
the scalar potential (x) (the temperature), with divergence replaced
by the curl operator. One may feel tempted to extend this analogy to
the boundary conditions as well, and it seems at first glance natural
to combine Eq. (2.13) not with (2.14), but rather with the boundary
condition


(2.19)
(x)  = 0 (x) = 12 Q x
V

the exact counterpart of (2.3). Indeed, if the volume V is homogeneous, then q(x) = 12 (Q x) Q, x V, so that (2.19)
corresponds obviously to a constant heat flux Q throughout V. Then,
in particular, q Q.
For a heterogeneous RVE, both (x) and q(x) fluctuate inside V.
However, the boundary conditions (2.19) yield the validity of (2.15)
as if the volume were homogeneous. This fact is a consequence of the
statement which is a direct counterpart of Lemma 2.1:
Lemma 2.3 Whatever the composition and microstructure of the
RVE, the volume average of the heat flux is completely determined
by the values, 0 (x), of its vector potential (x) at the boundary V.
The proof literally repeats the one of Lemma 2.1, having replaced
(x) by (x). In particular, if the function 0 (x) on V be
the one, given in Eq. (2.19), that corresponds to a constant heat flux
in a homogeneous volume V, then q = Q. Hence, the homogeneous
boundary condition (2.19) also assures the validity of (2.15), whatever
the internal constitution of the medium.
Similarly to the foregoing reasoning, we can solve Eq. (2.13) with
the boundary condition (2.19) and find the appropriate mean value of

32

Heterogeneous Media: Modelling and Simulation

the temperature gradient G = , see Eq. (2.17). The so-obtained

G is undoubtedly proportional to the prescribed mean heat flux Q:
G = k Q .

(2.20)

We cannot claim however that the constant k , as defined in

Eq. (2.20), coincides with the effective resistivity k , introduced in
Eq. (2.17). The reason is that the appropriate temperature gradients
G 6= G, since they result from solutions of Eq. (2.13), corresponding
to different boundary conditions, respectively, to (2.14) and (2.19).
The question then arises why (2.14) was preferred when defining the
effective resistivity in Eq. (2.17). The answer will be given in Section 1.2.5 below, when discussing the consistency of the definitions
of effective conductivity and resistivity, in the sense that they should
be each others inverse. Moreover, it will be shown there (Theorem
2.2) that the homogeneous boundary conditions (2.3) and (2.14) minimize the appropriate energies of the RVE as compared to all plausible
boundary conditions that give rise to prescribed mean temperature
gradient and heat flux respectively. (This will imply in particular
that k k .)
1.2.4 Energy Definition of the Effective Constants
Consider the energy of the RVE
W = W (x; (x)) = 12 q(x) (x) ,

(2.21)

W (x; (x)) = 12 (x)|(x)|2 .

Lemma 2.4 If the temperature on the boundary V is G x, so that
(x) = G, then
(2.22)
W = 12 q = 12 G2 .
Proof. Let

(x) = G + 0 (x) ,
q(x) = Q + q 0 (x) ,

G = ,
Q = q,

so that 0 (x) and q 0 (x) are the fluctuations of the indicated quantities about their mean values G and Q respectively. Then
W = 12 Q G + 12 q 0 (x) 0 (x) .

(2.23)

1.2 Micromechanics of Heterogeneous Media

33

But
1
V

1
(x) q (x) dx+
V
V

q 0 (x) 0 (x) =
1
=
V

q 0 (x) 0 (x) dx

(2.24)
qn0 (x)0 (x) dS

= 0,

because
q 0 (x) = (q(x) Q) = q(x) = 0 ,
and



0 (x) 

= ((x) G x) 

= 0.

It follows, in virtue of Eq. (2.23), that

W = 12 Q G = 12 G2 ,

(2.25)

having used the definition (2.10) of the effective conductivity of

the RVE.
It is noted that the above energy definition of the effective properties, used in a particular situation already by Einstein [EIa] (see the
footnote on page 16), appeared first in Hashins papers [HAc] and
[HAf] and, independently, in [HIb].
Remark 2.5 It is important to point out that in the linear case
under study, the function W (x; (x)) plays the role of a micropotential, i.e.
W (x; (x))
q(x) =
(x)
in each point x V. Lemma 2.5 then means that
Q=q=

W (G)
= G ,
G

so that the average potential W (G) = 12 G2 plays the role of the

macropotential. In other words, averaging the micropotential we get
just the macropotential; in the latter, the rapidly fluctuating conductivity field (x) is smoothed out being replaced by the effective
conductivity of the RVE.

34

Heterogeneous Media: Modelling and Simulation

Consider now the energy of the RVE as a function of the dual

variablethe heat flux q(x), i.e.
= 12 (x; q(x)) = 12 q(x) (x) ,
(x; q(x)) =

(2.26)

2
1
2 k(x)|q(x)| .

Note that the energies W and , as defined in Eqs. (2.21) and

(2.26), coincide, provided the appropriate fields are interconnected by
the Fourier law (1.5) and k(x) = 1/(x). We shall keep however these
different notations for them in order to underline that the energy
is considered as a functional of different argumentsthe temperature
and the heat flux fields respectively.
Lemma 2.5 Let there be no body sources in the RVE. If the heat
flux qn on the boundary V equals Qn, so that q = Q, see Eq. (2.15),
then
= 12 q = 12 k Q2 .
(2.27)
Proof. We repeat the reasoning of the foregoing proof of Lemma 2.4.
In the case under study we have again q 0 (x) 0 (x) = 0, since now
qn0 (x) = qn (x) Q n = 0 on V. That is why
= 12 Q G = 12 k Q2 ,

(2.28)

having used the definition (2.17) of the effective resistivity k of

the RVE.
Remark 2.6 Similarly to Remark 2.1, it is noted that in the linear
case under study the function (x; q(x)) = 12 k(x)|q(x)|2 plays the
role of a micropotential, i.e.
(x) =

(x; q(x))
q(x)

in each point x V. Lemma 2.4 then means that

G = =

(Q)
= k Q ,
Q

so that the average potential (Q) = 12 k Q2 plays the role of the

macropotential. In other words, averaging the micropotential we get
just the macropotential; in the latter the rapidly fluctuating resistivity field k(x) is smoothed out and replaced by the effective resistivity k of the RVE.

1.2 Micromechanics of Heterogeneous Media

35

Remark 2.7 Let q(x) be an arbitrary divergence-free field (not necessarily proportional to (x)). The two Lemmas 2.4 and 2.5 imply
that
q(x) (x) = q ,
(2.29)
either for uniform temperature gradient (x) = G x or for uniform
flux qn (x) = Q n on the boundary V. These two boundary conditions, either of which implies the validity of (2.29), are often called
in the literature Hills.
Note also that Eqs. (2.22) and (2.27) can serve as alternative definitions of the effective conductivity and resistivity, respectively, of
the RVE. Their equivalence to the usual definitions (2.10) and (2.17)
has several important consequences. One of them will be explored
immediately in the next subsection.
1.2.5 Consistency of the Definitions of Effective Properties
The possibility to introduce the effective properties on the base of
energy notions allows, first of all, a demonstration of the consistency for the definitions of the effective conductivity and the effective
resistivity, in the sense explained in the following important theorem:
Theorem 2.1 The effective conductivity and resistivity are each
others inverse, i.e.
1
k = .
(2.30)

Proof. Note first of all that the relations (2.22) and (2.27) do not
suffice to claim the validity of Eq. (2.30). The reason is the following. When solving Eq. (2.7) with the boundary condition (2.3), we
shall obtain a certain heat flux field q(x) which is not uniform, in
general, on the boundary V, i.e. qn (x) 6= Q n, where Q = q. The
effective resistivity, on the other hand, is defined in Eq. (2.17) under the assumption of homogeneous heat flux on the boundary, see
Eq. (2.19). The same remark holds true when considering the solution of the dual BVP (2.13), (2.19)it generates a temperature field
whose temperature gradient does not conform in general with the
needed, in the definition (2.10) of the effective conductivity, boundary condition (2.3). That is why the proof of Theorem 2.1 should
be preceded by appropriate comparisons between the energies of a

Heterogeneous Media: Modelling and Simulation

36

RVE, corresponding to a homogeneous and to a certain inhomogeneous temperature gradients and heat fluxes on the boundary V.7
Let first the heat flux be homogeneous on V, i.e. qn = Q n,
x V, so that q = Q. The solution of the BVP (2.13), (2.14), generates a temperature field, denoted by (a) (x), whose mean gradient
is
G = (a) = k Q .
Consider next a second temperature field, (b) (x), corresponding
to the uniform temperature gradient G x on the boundary V, i.e.
(b) (x) is the solution of the BVP (2.7), (2.3). Obviously, (b) = G,
similarly to the field (a) (x). In order to compare the energies, stored
by these two fields, consider the difference:
W (b) W (a) = W (x; (b) (x)) W (x; (a) (x))
i W (x; (a) (x))

(b) (x) (a) (x)

i

= (b) (a) q (a) = 0 ,

(a) (x)

q (a) (x) =

(2.31)

W (x; (a) (x))

.
(a) (x)

Two facts have been used here. The first is Eq. (2.29) which holds
since in our case q (a) (x) corresponds to a uniform heat flux, see Remark 2.7. The second is that the energy function W (x; (x)) is
convex with respect to (x) provided (x) > 0; that is why one
of the well-known definitions of convexity is applicable. (Namely, a
smooth enough scalar function of a vector argument f (u) is convex,
iff f (v) f (u) (v u) f (u), u , v.)
However, in the simplest case under study, when the flux is linearly
connected to the temperature gradient, the inequality, used in the
second line of Eq. (2.31), can be easily checked, without even referring
to convexity arguments. Indeed, in this case it reads
h

1
2 (x)

|(b) (x)|2 |(a) (x)|2

h

(x) (b) (x) (a) (x) (a) (x) ,

7
In a much more general variational setting this problem has been comprehensively studied by Willis et al. [TAc], [TOa], [WIf].

1.2 Micromechanics of Heterogeneous Media

37

i.e.
|(b) (x)|2 |(a) (x)|2 2(a) (x) (b) (x) 2|(a) (x)| 2 ,
which immediately follows from the obvious inequality

2


 (a) (x) (b) (x)  0 .

Thus, Eq. (2.31) implies that W (b) W (a) . But W (a) = 12 k Q2 ,

according to the definition of the effective resistivity and of the field
(a) (x), see Eq. (2.17), while W (b) = 12 G2 , and therefore
1 2
2 G

12 k Q2 ,

i.e. k 1 .

(2.32)

Dually, let (x) = G x on V; through solving the BVP (2.7),

(2.3), we find the respective temperature field (x) such that = G;
this field generates the heat flux, denoted by q (a) (x). The mean value
of the latter is
Q = q (a) = G .
Consider next a second heat flux field, q (b) (x), corresponding to
the uniform heat flux Q n on the boundary V, i.e. the solution of
the BVP (2.13), (2.14). Obviously, q (b) = Q, similarly to the field
q (a) (x). In order to compare the energies, stored by these two fields,
consider again the difference:
(b) (a) = (x; q (b) (x)) (x; q (a) (x))
i (x; q (a) (x))

q (b) (x) q (a) (x)

q (a) (x)

(2.33)

(x; q (a) (x))

.
q (a) (x)
We have used, in a manner fully similar to Eq. (2.31), the convexity
of the energy function (x; q(x)) with respect to q(x) (which surely
holds if k(x) > 0). Once again Eq. (2.29) was utilized, since in our
case (a) (x) corresponds to a uniform temperature gradient on the
boundary, see Remark 2.7. Thus, (b) (a) . But (a) = 12 G2 ,
according to the definition of the effective conductivity and of the
field (a) (x), see (2.10), while (b) = 12 k Q2 , and therefore
h

= q (b) q (a) (a) = 0 ,

1 2
2k Q

12 G2 ,

(a) (x) =

i.e. k 1 .

Comparison of (2.32) and (2.34) proves (2.30).

(2.34)

38

Heterogeneous Media: Modelling and Simulation

Remark 2.8 In the multitude of approximate theories of heterogeneous media, a minimum requirement, that each reasonable one
should meet, is the consistency of its predictions in the sense of
Eq. (2.30). In other words, if we perform the appropriate reasoning, say, for the effective conductivity and for the effective resistivity,
the results should be each others inverse. Without going into details
(which are purely technical) we should point out that all the oneparticle schemes, discussed in Section 1.5 below, are consistent. An
explicit check will be performed, for illustrative purposes, only for
the self-consistent theory of a polycrystal (see Section 1.6.2). More
detailed discussion can be found in [HAh] and [NEa].
Note that in the proof of the Theorem 2.1 the only fact we have
used for the second (comparison) fields (b) (x) and q(b) (x) is that
they comply with the conditions (2.4) and (2.15) respectively. That
is why, we can claim the validity of the following minimum principles:
Theorem 2.2
(i) Among all boundary data for the temperature that produce a field
(x) with a fixed mean gradient G, cf. Eq. (2.4), the uniform boundary conditions (2.3) render the energy W an absolute minimum
Wmin ; moreover Wmin = 12 G2 .
(ii) Among all boundary data for a divergence-free (heat flux) field
that produce a fixed mean flux Q, cf. Eq. (2.15), the uniform boundary conditions (2.14) render the energy an absolute minimum
min ; moreover min = 12 k Q2 .
These are the two basic minimum principles in the theory of heterogeneous media, based upon the notion of volume averaging. We
shall return to them in Section 1.3.4 in connection with the classical
variational principles of Dirichlet and Thompson.
Theorems 2.1 and 2.2 clarify the special role of the homogeneous
boundary conditions (2.3) and (2.14) as well as the reason why they
have been adopted when defining the effective conductivity and resistivity. They also explain the reason why a preference was given to
to the boundary condition (2.14) when defining the effective conductivity in Eq. (2.17), and not to its rival (2.19) (which at first glance
looked equally appropriate, generating a homogeneous flux field in a
homogeneous medium).

1.2 Micromechanics of Heterogeneous Media

39

Remark 2.9 Obviously, the central point of the proof of Theorem 2.1
are the inequalities (2.31) and (2.33) which hold for any convex potential W and . Hence both theorems can be appropriately generalized
for mixtures of nonlinearly behaving media, provided the constitutive
relations for the constituents emerge from convex potentials (see the
book of Nemat-Nasser and Hori [NEa, pp. 4749] for details and a
proof). Here, for the sake of simplicity, the linear scalar case has been
only treated. As pointed out by the authors [NEa, pp. 4748], their
proof had been suggested by Willis in a private communication. It is
noted also that the above Theorems 2.1 and 2.2 are particular cases
of much more general results of Willis et al. [TAc], [TOa], [WIf],
which explains why they will be called Willis.
1.2.6 The Concentration Factors
The effective constants can be conveniently represented by means of
the so-called concentration factors. Determination of these factors, to be defined below, is a problem, equivalent to that of the said
constants. However, the appropriate representations provide a natural and useful guide when constructing approximate formulae for the
effective properties, as we shall see in Section 1.5.
[Link] Ideal Thermal Contacts
Let again V be a RVE of a two-phase medium. Then V = V1 V2 ,
where V1 and V2 are the subvolumes, occupied by the constituent i,
respectively, i = 1, 2.
We shall begin with the scalar conductivity context again. Let
us first consider the most common case when an ideal contact takes
place, so that both temperature and heat flux remain continuous in
the whole volume. Then obviously
= 1 hi1 + 2 hi2 ,

(2.35)

q = 1 hqi1 + 2 hqi2 ,
throughout the volume V. (If the temperature and heat flux have discontinuities, as it is the case with nonideal contacts, additional terms
should be added to the right-hand sides of (2.35), see Eqs. (2.48) and
(2.51) below.) In Eq. (2.35)
hii =

1
Vi

Z
Vi

dx ,

(2.36)

40

Heterogeneous Media: Modelling and Simulation

i.e. hii signify the spatial averages over the regions Vi , occupied by
the constituent i, Vi = vol Vi , i = 1, 2. According to the Fourier law
and the definition of the effective conductivity, see Eq. (2.10),
Q = q = K ,

hqii = i hii ,

i = 1, 2 .

(2.37)

From (2.35) and (2.37) we can exclude the mean temperature gradient
in one the phases, say, in 1:
Q = K = q = 1 hqi1 + 2 hqi2
= 1 1 hi1 + 2 2 hi2


= 1 2 hi2 + 2 2 hi2 ,
or
K = 1 + 2 [] hi2 ;

(2.38)

recall that the square brackets denote the appropriate jumps in material properties, e.g. [] = 2 1 , [L] = L2 L1 , etc.
Let us prescribe, similarly to Section 1.2.2, the average value
G = of the macroscopic gradient, having imposed the boundary
condition (2.3) for the temperature. The solution of the boundary
value problem (BVP) (2.7), (2.3) depends linearly on the vector G.
Then, in particular,
hii = Ai G ,

i = 1, 2 ,

(2.39)

with certain concentration factors A1 and A2 that represent, in the

case under study, second-rank tensors. (In the scalar and isotropic
case the Ai s reduce to constants, called in [REb] field factors.)
Note that the concentration factors Ai are simply interconnected:
1 A1 + 2 A2 = I ,

(2.40)

see their definition (2.39) and the first line of Eq. (2.35).
Together with Eq. (2.39) (at i = 2), Eq. (2.38) yields the following
formula for the effective conductivity tensor of the composite:
K = 1 I + 2 []A2 ,

(2.41)

or, in an equivalent, but more symmetric form,

1 (K 1 I) A1 + 2 (K 2 I) A2 = 0 ,
which follows from Eq. (2.40).

(2.42)

1.2 Micromechanics of Heterogeneous Media

41

Remark 2.10 An alternative and physically appealing starting point

can be the observation of Landau and Lifshitz [LAa] that the field
q(x) 1 (x) is non-zero only within the inclusions. That is why
its average value,
q 1 = Q 1 G = [K 1 I] G ,
is proportional to the number density n of the inclusions, multiplied
by the average value of the same quantity within a single inclusion W:
[ K 1 I ] G = 2 hq 1 iW = 2 [] hiW ,

(2.43)

having assumed the inclusions to be of the same shape W, 2 = nW ,

W = vol W. Eq. (2.43) obviously reproduces the basic formula (2.41),
if (2.39) is taken into account.
Both Eqs. (2.41) and (2.42) look exceptionally simple. But this
simplicity should not mislead the reader. The extreme difficulties in
the homogenization problem, already discussed in Section 1.1, are
not resolved by Eqs. (2.41) or (2.42); they are only moved from one
place to another. That is, instead of looking for the effective conductivity, the quest is now for one of the tensors A1 or A2 . The
specification of any one of them needs knowledge of the temperature
fields within the appropriate constituent, say, within all of the filler
particles. This means that, as a matter of fact, again the detailed solution of the BVP, Eqs. (2.7) and (2.3), is needed. Hence, it is little
wonder that the concentration factors can be rigorously determined
in a very few cases; moreover, any such case leads to a remarkable
formula that gives the effective properties exactly, under minimum
assumptions about the internal structure of the heterogeneous solids.
Two such famous formulae will be discussed in the next Section 1.3.
The first one is the Hill expression for the bulk modulus of a twophase elastic medium, whose constituents possess equal shear moduli,
see Eq. (3.66). The other is Levins result that interconnects the effective bulk modulus with the effective thermal expansion coefficient,
see (3.76).
The formula (2.41) is very convenient, however, for obtaining approximate expressions for the effective properties. The basic idea

Heterogeneous Media: Modelling and Simulation

42

is to replace the concentration factor A2 in Eq. (2.41) by an appropriate tensor, stemming appropriately from the solution of a single inhomogeneity problem, on the base of certain heuristic arguments. The approximations, devised in such a manner, can be called
one-particlethey will be discussed in Section 1.5 having, to this
end, considered in needed detail the single-inclusion problem in Section 1.4.
Remark 2.11 The simple and important formula (2.41) is usually
attributed to Hill [HIb]. However, in the scalar context it was derived
earlier by Reynolds and Hough [REb], who clearly understood and
demonstrated its key role in formulating some of the one-particle
approximations to be treated below.
[Link] The Wu Representation
An alternative expression of the effective properties through the concentration factors can be obtained, following a reasoning, sketched
by Wu [WUa].
Let q(x) be the flux in the heterogeneous volume V, under the homogeneous boundary conditions (2.14), and let (x) = k(x)q(x)
be the temperature gradient that appears. Consider a comparison
(or a reference) medium with the resistivity k1 and impose mentally the same flux field q(x) there; as a result, a certain fictious
temperature gradient g c (x) will show up:
g c (x) = k1 q(x) =

(x) ,

if x V1 ,

k2
(x) , if x V2 .
k1

(2.44)

By means of the field g c (x), rearrange the energy (2.26) as

1
=
2V
1
+
2V

1
q(x) (x) dx =
2V
V

Z h

Q g c (x) dx

(2.45)
i

q(x) (x) Q g (x) dx .

The Gauss theorem and the boundary conditions (2.14) yield

Z
V

q(x) (x) dx =

Z
V

Q (x) dx .

1.2 Micromechanics of Heterogeneous Media

43

Also
= 12 k Q2 ,

g c = k1 q = Q ,

see Eqs. (2.26) and (2.44). Hence from (2.45) it follows

1 2
2k Q

2
1
2 k1 Q

1
+
2V

Q p(x) dx ,

(2.46)

p(x) = (x) g (x) ,

where p(x) is the so-called polarization field. Obviously
p(x) = [(x) 1 ] q(x)
which vanishes in the phase 1, see Eq. (2.44). That is why (2.46)
can be recast as
1 2
2k Q

= 12 k1 Q2 + 12 2 [k] Q hqi2 .

But
hqi2 = 2 hi2 ,

hi2 = a2 G ,

where A2 = a2 I is the concentration factor (2.39) in the isotropic

case under study. In turn, G = k Q, and from the last formula the
simplest Wu representation emerges:
k = k1 + 2 [k] k 2 a2 .

(2.47)

The appropriate generalizations, concerning anisotropic constituents and/or internal alignment, as well as elastic behavior, are obvious. (Wu himself [WUa] considered the elastic case.) However, all
the obtained in such a way formulae will be unavoidably equivalent
to (2.41): if they are not, by chance, we would have two independent
relations for the effective tensor and for the concentration factor and
they would allow us to specify them both explicitly. The latter is
possible but for few simplest cases only such as laminate media, see
Remarks 2.1 and 2.2. (The equivalence of Eqs. (2.47) and (2.41) in
the isotropic case can be easily checked.)
The above simple considerations possess, however, the merit of
invoking the notions of reference media and polarization fields, and
that is why they have been included here. The polarization fields are
a useful tool in the theory of random heterogeneous media. One of the

Heterogeneous Media: Modelling and Simulation

44

reasons is that the basic BVP problems, like (2.7), (2.3) can be easily
recast as integral equations with respect to the polarizations. These
equations can be then reformulated variationally in an obvious way
and in the result the important and very useful variational principles
of the so-called Hashin-Shtrikman type emerge, see the survey [WIc]
or [WUb].
Simple arguments, employing ingeniously appropriate reference
media, can be also utilized to construct the basic approximate schemes
of Section 1.5 for the effective properties of heterogeneous media. Details can be found in the paper [BEo].
[Link] Nonideal Thermal Contacts
In the case of a nonideal contact the heat flux and/or the temperature are discontinuous at the interphase boundaries S12 within the
volume V. Accordingly, there are two main types of such contacts
superconducting and resistive, to be properly defined and briefly discussed in Section 1.4.3 below, having introduced the notion of a singular inclusion after Kanaun [KAc]. For the moment it suffices to
demonstrate only that the discontinuities of the flux or temperature
lead to the appearance of additional terms in Eqs. (2.35).
Indeed, let us start with a superconducting interphase, as introduced in Section 1.4.3. Its exact definition is not important here;
we need only the fact that the temperature is continuous, but the
heat flux has jumps on S12 in this case, see (4.37). Application of
Eq. (2.16) then leads to the formula, given recently by Miloh and
Benveniste [MIa]
1
q=
V

Z
V

q(x) dx = 1 hqi1 + 2 hqi2 + q (12) ,

q

(12)

1
=
V

(2.48)

Z
S12

[qn ] x dS ,

where q (i) is the flux in the phase i, i = 1, 2, and n is the unit

normal vector of the interphase surface, pointing from phase 2 to
phase 1; hence [qn ] = q (1) n q (2) n is the jump of the normal
component of the flux across S12 . At the same time, the first line of
Eq. (2.35) remains unchanged due to continuity of the temperature
in the case under study.

1.2 Micromechanics of Heterogeneous Media

45

Similarly to Eq. (2.39), introduce the interphase or surface concentration factor, A(12) , through the relation
q (12) = A(12) G .

(2.49)

Repeating literally the elementary reasoning, that led us to (2.41),

produces its superconducting counterpart
K = 1 I + 2 []A2 + A(12) .

(2.50)

A fully similar formula can be devised for the the other basic kind
of nonideal contact, called resistance. Again the exact definition,
given in Section 1.4.3 below, is not important in this moment; we
need only the fact that the heat flux is continuous in this case (so
that the second line of Eq. (2.35) remains now unchanged), but the
temperature field has a jump in this case. An application of the
Gauss theorem then yields
=

1
V

Z
V

(x) dx = 1 hi1 + 2 hi2 +

(12)

1
=
V

(12)

,
(2.51)

(1)

(2)

n dS ,

S12

where (i) is the temperature field in the phase i, i = 1, 2, and n is

again the unit normal vector on the interphase surface S12 , pointing
from phase 2 to phase 1, see [BEf] for a detailed derivation and
discussion. The counterpart of (2.49) is

(12)

= B (12) G ,

(2.52)

with B (12) representing the interphase concentration factor for the

resistance boundary under study. Repeating once again the foregoing reasoning yields now
K = 1 I + 2 []A2 1 B (12) .

(2.53)

1.2.7 The Elastic Case

The generalization of the foregoing scalar reasoning to the elastic
case is straightforward, and will be outlined now.

46

Heterogeneous Media: Modelling and Simulation

[Link] The Basic Equations

For a (linear) elastic solid the basic system Eq. (2.7) is replaced by
its tensorial counterpart
(x) = 0 ,
(x) = L(x) : (x) ,

= 21 (u + u) ,

(2.54)

so that (x) is the strain tensor, generated by the displacement

field u(x), (x) is the stress tensor, L(x) is the fourth-rank tensor
of elastic moduli and the colon stands for contraction with respect
to two pairs of indices. The analogy with the scalar conductivity
case, Eq. (1.5), is obvious: temperature is replaced by displacement,
temperature gradient by strain tensor, heat flux by stress tensor,
conductivity coefficient (x) by the elastic moduli tensor L(x) and
the resistivity k(x) = 1/(x) by the tensor of elastic compliance
M(x) = L1 (x).
All scalar conductivity considerations of Sections 1.2.2 to 1.2.5
are straightforwardly extended to the elastic case. For instance, the
counterparts of Lemma 2.1 and its corollary, Eq. (2.6), read:
Lemma 2.6 Whatever the composition and microstructure of the
RVE, the volume average of the strain tensor is completely determined by the displacement field on the boundary V.
Proof. Let u(x) be a displacement field in V, attaining the prescribed
boundary values u0 (x). Using the Gauss theorem, we have
1
=
V

1
(x) dx =
V
V

1
=
2V
u0 (x)

In turn, if
=
second-rank, then

since x = x = I.
Similarly we have

nu0 + u0 n dS .

x = x 0 , where 0 is a symmetric tensor of

n0 x + (0 x)n dS

Z 
1

2V

(u + u) dx
V

1
=
2V
=

0 x + 0 (x) dx = 0 ,

1.2 Micromechanics of Heterogeneous Media

47

Lemma 2.7 Whatever the composition and microstructure of RVE,

the volume average of the stress tensor, at the absence of body forces,
is completely determined by the boundary values of the traction t0 =
n on V.
Proof. Let (x) be a stress field in V, with the prescribed boundary
values of the traction. Due to the assumed absence of body forces,
(x) = 0 In the Cartesian system xi , i = 1, 2, 3, we have, using
the Gauss theorem and the equilibrium equation,
1
V

ij =
=

1
V

Z
V

Z
V

ij (x) dx =

1
V

(xi kj (x)),k dx

1
V

ik kj (x) dx =

1
V

xi kj,k (x) dx =

Z
V

1
V

xi,k kj (x) dx

xi t0j (x) dS ,

V
0
tn on

which indeed is uniquely specified by the traction

the boundary V.
In turn, if t0n = 0 n on V, i.e. the traction is uniform on the
boundary,
1
ij =
V

1
V

0
np pj
xi dS

0
(pj
xi ),p dx =

1 0

V pj

Z
V

0
0
xi,p pj
ip dx = ij
,

i.e. = 0 .
[Link] The Effective Elastic and Compliance Tensors
Similarly to the scalar case, consider first the linear surface displacements


u0 (x)  = 0 x ,
(2.55)
V

0 ,

which assures that =

due to Lemma 2.6.
The BVP (2.54), (2.55) possesses a solution which is unique (the
tensor field L(x) is positive-definite for all x V), which generates a
certain stress field (x). Then
=

1
V

Z
V

L(x) : (x) dx = L : 0 .

(2.56)

The tensor L that enters Eq. (2.56) is the effective elastic tensor of
the volume V.

Heterogeneous Media: Modelling and Simulation

48

Dually, the stress field (x) satisfies the equation

(x) = (M(x) : (x)) = 0 .

(2.57)

Assume uniform boundary tractions:




t0n (x) = (x) n 

= (x) n ,

(2.58)

which assures that = 0 , due to Lemma 2.7.

The BVP (2.57), (2.58) possesses a solution which is unique (the
compliance tensor M(x) = L1 (x) is positive definite for all x V,
once L(x) is such) which generates a certain strain field (x). Then
=

1
V

Z
V

M(x) : (x) dx = M : 0 .

(2.59)

The tensor M that enters Eq. (2.56) is the effective compliance tensor
of the volume V.
Let
W (x; (x)) = 12 (x) : L(x) : (x)
be the elastic energy at the point x V. If the linear surface displacements (2.55) are imposed, then
W = W (0 ) = 21 (x) : L(x) : (x) = 12 0 : L : 0 .

(2.60)

Dually, let
(x; (x)) = 21 (x) : M(x) : (x)
be the elastic energy at the point x V, considered as a function
of the local stresses. If the surface tractions are uniform, see (2.58),
then
= ( 0 ) = 12 (x) : M(x) : (x) = 12 0 : M : 0 .

(2.61)

Eqs. (2.60) and (2.61) are the elastic counterparts of the appropriate scalar equations (2.22) and (2.27), see Lemmas 2.4 and 2.5.
Making use of them, one can show that

M = L1 ,

(2.62)

1.2 Micromechanics of Heterogeneous Media

49

i.e. the effective tensors of elasticity and compliance are each others
inverse; the proof repeats almost literally that of the foregoing Theorem 2.1. The minimum properties of the homogeneous boundary
conditions (2.55) and (2.58) can be then formulated in a manner,
fully similar to the statements of Theorem 2.2.
[Link] The Elastic Concentration Factors
In the elastic case only ideal contacts will be treated for which both
displacement and normal tractions are continuous at the interphases.
In particular, the phases are assumed hereafter firmly bonded, hence
no slipping and/or delimination across the interphase boundaries are
allowed. The scalar considerations of Section [Link] can be then
literally repeated, with the only difference that the second-rank concentrations factors Ai are to be replaced with fourth-rank tensors Ai
such that
hii = Ai : 0 , i = 1, 2 .
(2.63)
These tensors satisfy the identity
1 A1 + 2 A2 = I ,
where I is the unit fourth-rank tensor. In Eq. (2.63) 0 is the prescribed uniform strain, that stems from the homogeneous boundary
condition (2.55).
The effective elastic tensor, for a two-phase medium, now reads

L = L1 + 2 [L] : A2 ,

(2.64)

which is the counterpart of the scalar formula (2.41); L i are the

tensors of elastic moduli of the constituents, i = 1, 2.
Remark 2.12 Recall that any realistic tensor of the elastic moduli should possess, besides the symmetry in the first and the second
pairs of indices, the external symmetry
Lijkl = Lklij ,

(2.65)

as a consequence of natural thermodynamical requirements. This

means that not any fourth-rank tensor can serve as a concentration
factor for a heterogeneous medium, even if it satisfies (2.65). The

50

Heterogeneous Media: Modelling and Simulation

reason is that the product [L] : A2 of two tensors with such a symmetry is not obliged to comply with (2.65). Hence, a care is needed
when constructing approximate theories, based on the formula (2.64),
since they can violate (2.65) for some special internal constitutions
(say, certain three-phase media), see Section [Link].
Remark 2.13 Nonideal contacts of the same superconducting and
resistive types can be considered in the elastic case as well, extending the reasoning of Section [Link]. The definition of such contacts due to Kanaun et al. [KAc], [KAd] is outlined below, see Section 1.4.3. An extensive survey of the approaches, results and references, concerning nonideal interphases in elasticity and solid mechanics can be found in the book of Theocharis [THa], or in [EBa],
where the emphasis is more on their material science aspects.
1.2.8 The Effective Absorption Coefficient
Here we shall consider the homogenization quest for Smoluchowskis
absorption problem in a heterogeneous medium (see Section [Link]).
Formally, it is more convenient and instructive to deal, following
Talbot and Willis [TAb], with the more general problem, described
by the equation
c(x) k 2 (x) c(x) + K = 0 ,

xV.

(2.66)

Additionally, the no-flux condition



c(x) 
=0
n  V

(2.67)

is imposed on the boundary V.

The problem (2.66), (2.67) describes, for example, the steadystate diffusion of a species, created at the constant rate K within the
volume V, and absorbed there with variable intensity (defined by the
absorption coefficient field k 2 (x)). For a binary medium, to be only
treated here,
(

k 2 (x) =

k12 , if x constituent 1 ,
k22 , if x constituent 2 ,

(2.68)

similarly to the conductivity field in Eq. (1.6). Hence the species

is absorbed at different rates by the two constituents. Note that

1.2 Micromechanics of Heterogeneous Media

51

the units are chosen so that in Eq. (2.66) the diffusion coefficient
D = 1. The condition (2.67) guarantees that the species, created in
V, remains within V.
A variety of different physical phenomena can be modelled by the
system (2.66), (2.67). First, the above mentioned Smoluchowskis
problem is recovered if one take the limits:
k12 0 ,

k22 .

(2.69)

This limit, called by Willis hard (due, in particular, to the hard technical difficulties it involves), corresponds obviously to the situation
when one of the constituents (the matrix 1) does not absorb at all,
while the inclusions (the phase 2) are perfect absorbers. This
very important case is thoroughly reviewed in Torquatos survey in
the present volume (Chapter 2), together with the most recent developments. Boundary conditions more general than (2.67) are discussed
there as well.
Note also that the limit (2.69) can be interpreted in terms of heat
removal of a solid in the steady-state limit, when a constant heating of
intensity K takes place, competing with a system of coolers (the
phase 2) kept at the constant (zero) temperature, see, e.g. [KOa].
In the case of a general binary field k 2 (x) an important application
of Eq. (2.66) is supplied by the field of irradiation damage: When
a solid is subjected to intense, say neutron radiation, defects like
vacancies and intersitials are created. They diffuse in the medium and
either recombine or are absorbed (trapped) by the sinks like voids,
dislocations or grain boundaries. This leads to a considerable change
of the microstructure and, in particular, to swelling of the voids and
ultimate rupture. A deep review of this important interpretation and
application of Eq. (2.66) can be found in the survey of Brailsford and
Bullough [BRa].
To homogenize Eq. (2.66) average both sides over the volume V:
c(x) k 2 (x)c(x) + K = 0 .
The Gauss theorem implies
1
c(x) =
V

Z
V

c
dS = 0 ,
n

Heterogeneous Media: Modelling and Simulation

52

due to the no-flux condition (2.67). Hence

2

k c=

k 2 (x)c(x)

=K,

1
c=
V

c(x) dx ,

(2.70)

which defines the effective absorption coefficient, k 2 , of the medium.

The sense of Eq. (2.70) is fully similar to that of the definition
(2.10) for the effective conductivity. It means, let us underline once
again, that from a macroscopic point of view the volume V looks homogeneous, with a certain effective absorption coefficient k 2 . The
latter specifies the rate of macroscopical absorption of the species,
i.e. the proportionality coefficient (the rate constant) between the
mean steady-state concentration c and the intensity K of their creation.
Remark 2.14 It is worth pointing out that a simple change of the
sign in the left-hand side of the Helmholtz equation (2.66), with a special form of the source term, drastically changes the physical context.
Namely, we recall that the equation
(x) + k02 n2 (x) (x) + (x) = 0 .

(2.71)

describes the amplitude of a scalar wave, generated by a harmonic

point source in the origin, in a lossless, isotropic, time-independent
medium. In Eq. (2.71) k0 is the free-space wave-number and n(x) is
the fluctuating index of refraction.
The homogenization of
Eq. (2.71) is an important and classical problem of micromechanics and physics which cannot be discussed here. (The relevant literature is very extensive, see e.g. the survey of Frisch [FRc].) We
shall only add that many of the basic ideas and approaches for an
approximate evaluation of the effective properties, like self-consistent
scheme, effective field, etc., are very well adapted for wave propagation in heterogeneous media and, in particular, for Eq. (2.71), as
discussed in detail in Kanauns paper in the present volume (Chapter 3). Moreover, a considerable number of them have their direct
origin in the study of wave phenomena in such media, as initiated by
the works of Foldy [FOa] and Lax [LAd].
Note finally that the rate constant k 2 can be expressed by means
of the appropriate concentration factors in a manner, fully similar to

1.3 Micromechanics of Heterogeneous Media

53

the conductivity and elastic cases. Indeed, from (2.70) it follows

c=

1
V

Z
V
2

c(x) dx = 1 hci1 + 2 hci2 ,

K=k c=

1 k12 hci1

(2.72)

+ 2 hci2 ,

with the same notations (2.36), cf. Eq. (2.35).

Eliminating hci1 from (2.72) gives
k 2 c = k12 c + 2 [k 2 ] hci2 ,

[k 2 ] = k22 k12 .

(2.73)

This exact relation is fully similar to Eq. (2.38). In the present context it defines the effective sink strength k 2 through the average
concentration hci2 within the inclusions (the phase 2). Hence, if
hci2 = A2 c ,

(2.74)

k 2 = k12 + 2 [k 2 ]A2 ,

(2.75)

then
which is the absorption counterpart of Eqs. (2.41) and (2.64) and,
obviously, A2 is just the concentration factor for the problem under
study.
Similarly to the scalar and elastic cases, the formula (2.75) can
serve as a convenient starting point when deriving one-particle approximations for the effective constant k 2 . One of them (the selfconsistent one) will be considered in more detail in Section 1.5.7.

1.3 Some Basic Results

In this Section we shall collect several basic results for the effective
properties of heterogeneous media. The so-called weakly inhomogeneous medium will be first treated. Though straightforward, its
study clarifies why the homogenization problem is so difficult in general: The main reason is that the effective parameters to be specified, though very simple-looking scalar (or tensor) quantities at first
glance, incorporate within them all details about the internal structure of the mediuma tremendous amount of information, that is.

54

Heterogeneous Media: Modelling and Simulation

In practice we however know, as a rule, only a very limited part of

the latter. Hence the best we can hope about is, first, to find a certain interval within which the effective parameters should be located
(the narrower the better, of course). Such an interval is provided in
many cases (but not always, as we shall see) by appropriate variational principles; from them the so-called variational estimates for the
effective properties result. A brief account will be given below for the
classical such estimates, like those of Wiener, Hill, Hashin-Shtrikman
and Beran.
Second, with a limited information at hand, we can try to construct certain approximate scheme, using heuristic arguments (like
the ones employed by Maxwell, see Section [Link], when deriving
his famous formula (1.19)). A systematic way of constructing a wide
class of such approximate formulae will be discussed in Section 1.5.
1.3.1 Weakly Inhomogeneous Media
The weakly inhomogeneous media consist of phases whose properties
do not differ much. For simplicitys sake, we shall additionally assume that they are intermixed in such a manner that the mixture is
macroscopically isotropic.
[Link] The Definition
More precisely, let (x) be the varying conductivity field of the
medium. Represent it in the form
(x) = + 0 (x) ,

(3.1)

where 0 (x) is the fluctuation of (x) about its mean value , so that
0 = 0. The medium is weakly inhomogeneous, if
max
xV

|0 (x)|
 1.

For the two-phase materials treated here, the equivalent condition,

to be used below, reads
||  1 ,

= []/ .

(3.2)

It is clear from the definition (3.2) that weakly inhomogeneous

media are of very limited practical importance. However, the brief

1.3 Micromechanics of Heterogeneous Media

55

digression toward their study now is warranted by several facts of

primary importance in the general theory, developed in the sequel:
(i) A simple formula for , valid to the second-order ()2 holds;
it will play a central role in deriving the famous Hashin-Shtrikmans
estimates on , already applicable for arbitrary two-phase media,
see Section 1.3.3 below.
(ii) The performed analysis will clearly demonstrate that the effective
conductivity is indeed a quantity that depends, in general, on all
details of the internal constitution of a heterogeneous medium.
(iii) More precisely, it will be seen that macroscopically the internal
constitution shows up through certain specific integral parameters
that incorporate the so-called correlation functions for the medium.
[Link] The Perturbation Expansion and the Second-order
Approximation
The appearance of the small parameter, , in the BVP (2.7), (2.3),
suggests to look, after Brown [BRd], for the solution as the formal
perturbation expansion
(x) = (0) (x) + (1) (x) + ,



(p) (x) 

(0) (x) = G x ,
(3.3)

(p) (x) O (()p ) ,

= 0,

p = 1, 2, . . . . Then
= +

Bp ,

Bp G = 0 (x)(p) (x) ,

(3.4)

p=1

p = 1, 2, . . . , as it follows from Eqs. (2.10), (3.1) and (3.3).

The terms of the series (3.3) are specified by means of the equations
(1) (x) + G 0 (x) = 0 ,


(1) (x) 


(2) (x) + 0 (x)(1) (x) = 0 ,

= 0,

(3.5)



(2) (x) 

= 0 , (3.6)

etc., derived in a fully standard way.

Truncation of the series (3.3) after the first term, i.e. taking (x)
G x, results in the already discussed Voigt approximation (2.11).

56

Heterogeneous Media: Modelling and Simulation

The simplest nontrivial case shows up, if the series (3.3) is truncated
after the second term. Then the perturbation (1) (x) is governed by
Eq. (3.5). In turn


q = G = G + 0 (x)(1) (x) + O ()2 ,

(3.7)

cf. Eq. (2.10), which means that the performed truncation

specifies


the effective conductivity exactly to the order O ()2 , once the
mean value B1 G = 0 (x)(1) (x) is found.
To find this value, it is noted that the solution of the BVP (3.5)
can be written as
(1) (x) = G

Z
V

0 (y)
dy .
4|x y|

(3.8)

This is clearly an approximation, since the Green function 1/(4|x|)

for the Laplace equation in the unbounded space is used; in other
words, we have tacitly assumed, when solving Eq. (3.5), that the
RVE coincides with the whole space R3 .
Eq. (3.8) allows us to represent now the quantity, needed in (3.7),
in the form
0 (x)(1) (x) = G

Z
V

0 (x)0 (x y) x x

1
dy .
4|y|

(3.9)

The assumptions of statistical homogeneity and isotropy, as discussed

in Section 1.2.1, imply in particular that the two-point average
0 (x)0 (y) depends only on the distance between x and y, excepts
for a thin layer near the boundary V . To be consistent with usage of
the infinite-body Green function, we should neglect this layer, taking
in (3.9) the integration domain as the entire R3 . Thus
B1 G = 0 (x)(1) (x)
=G

Z
R3

M2 (r)

1
M2 (0)
dx =
G,
4r
3

(3.10)

since the integral in the last formula represents an isotropic secondrank tensor and (1/(4|x|)) = (x), with (x) denoting the Dirac
delta-function. In Eq. (3.10)
M2 (r) = 0 (0)0 (z) ,

r = |z| .

(3.11)

1.3 Micromechanics of Heterogeneous Media

57

Note that for any binary function a(x) that takes the values a1
and a2 in the phases V1 and V2 , respectively, one has
a0 2 = (a a)2 = 1 22 [a]2 ,
a0 3

= (a

a)3

1 22 (1

(3.12)
3

2 )[a] ,

[a] = a2 a1 , as a simple check demonstrates.

Hence, for the two-phase medium under study,
M2 (0) = 0 2 = 1 2 []2 ,

B1 =

1 2 []2
.
3

(3.13)

Together with (3.7) and (3.10), this gives the simple, but very important formula,

[]
= 1 13 1 2

2

+ O ()2 ,

(3.14)

which, in particular, will be needed in Section [Link], when deriving

the Hashin-Shtrikmans estimates on the effective conductivity.
To the best of the authors knowledge, the formula (3.14) was
first given by Brown [BRd], and rederived by many authors afterward
(especially, see the book [LAa] for an alternative and more heuristic
derivation).
Remark 3.1 The central role of the formula (3.8), and of its consequence (3.14) within the frame of the volume averaging approach,
is to be particularly emphasized. From them one can already rigorously derive both Hashin-Shtrikmans and Berans estimates, as we
shall see in a moment.
Remark 3.2 The derivation of (3.10) demonstrates in passing
that for a reader with a more rigorous taste, the approach adopted
here has profound inherent problems, connected with the employed
notion of very big, but still finite RVEs V (see also the comments
of Willis [WIc, p. 15]). Much more consistent is the stochastic point of
view, employed in the pioneering Browns paper [BRd], in which the
fields like (x), (x), etc., are treated as random. But this approach,
as already pointed out, is not pursued here.

Heterogeneous Media: Modelling and Simulation

58

It is noted that the function M2 (r), defined in Eq. (3.11) is the

so-called two-point correlation for the medium. The term can be
explained by the fact that, while the one-point average = 1 1 +
2 2 provides information about the volume fractions i only (at fixed
values of i ), the two-point average 0 (x)0 (y) supplies already
much more detailed two-point information, concerning the spatial
distribution of the constituents within the heterogeneous medium.
It is a simple exercise to show that M2 (r) is directly connected to
the average amount of the phase, say, 1, located at the distance r
from the origin, provided the same phase is to be found there. In
stochastic terms M2 (r) allows evaluation of the probability that if a
rod of length r is thrown into the medium, both its ends fall within
one and the same constituent, see [FRb], [BEi], [TOe] et al.
Remark 3.3 For the integral (3.10) to be convergent, we have tacitly assumed that M2 (r) decays as r . The meaning of this
assumption is that there is no long-range order in the medium;
in the language of elementary probability theory the properties of
the medium at the points x and y and, in particular, the quantities
(x) and (y), become uncorrelated as the distance between them
increases, r = |x y| .
The formula (3.14) shows that to the order O ()2 the effective
conductivity of a macro-isotropic medium does not feel specific details of the internal structure. Even the exact form of the two-point
correlation function M2 (r) for r > 0 is not required. Only the value
of the latter at the origin, M2 (0), that is the volume fractions i ,
affects
to this order. To the higher orders however, starting with
O ()3 , this is already not true.

[Link] The Third-order Approximation

Indeed, to see this it suffices to look at the equations that govern
the higher-order terms in the expansion (3.3). For example, let us
truncate (3.3) after the third term. Then
!

G=

0 2

+ B2 G ,
3

B2 G = 0 (x)(2) (x) .

(3.15)

The function (2) (x) is governed by Eq. (3.6) and it can be easily
found, if we decide to employ Eq. (3.8) (replacing again the volume

1.3 Micromechanics of Heterogeneous Media

59

V with the entire R3 ). In particular,

B2 =

0 3
I ,
32 3

0 3 = 1 2 (1 2 )[]3 ,

(3.16)

cf. Eq. (3.12), so that

=



0 2
0 3
+ 2 I3 + O ()3 .
3
3

(3.17)

In the last formula

I3 =

1
0 3

ZZ

M3 (x, y)

1
1
:
dx dy
4|x|
4|y|

(3.18)

is a dimensionless parameter that first appeared (in different notation) in Browns paper [BRd]. (The notation adopted here aims at
preserving similarity between the second- and third-order terms in
the expansion (3.17), cf. Eq. (3.9) as well.) Also, in Eq. (3.18)
M3 (x, y) = 0 (0)0 (x)0 (y)

(3.19)

is the three-point counterpart of the function M2 (r), defined in

Eq. (3.11). Naturally, M3 (x, y) can be called three-point correlation.
It supplies more detailed information about the internal structure of
the medium, as compared to the volume fraction and the two-point
correlation M2 (r). (In probabilistic terms M3 (x, y) allows evaluation
of the probability that if a triangle, two sides of which coincide with
the vectors x and y, is thrown into the medium, all its vertices fall
within one and the same constituent, see [BEi], [TOe] et al.)
Note that I3 is a very important characteristic of a two-phase heterogeneous medium: as we shall see, the same quantity will reappear
in a certain simple and natural Ritzs type variational procedure for
bounding , proposed by Beran [BEg], see Section [Link] below.
Hence, the truncation of the series (3.3) after the third term de

termines the effective conductivity exactly to the order O ()3 ,
provided the parameter I3 is found. Its specification, however, needs
knowledge not only of the value M3 (0, 0) at the origin, but of the
whole function M3 (x, y) for all values of x and y (unlike the situation with the function M2 (r) and the corresponding to it second-order
term in the expansion (3.3)).

60

Heterogeneous Media: Modelling and Simulation

The evaluation of the higher-order terms Bp , p > 2, in (3.4) is

fully similar: If the validity of Eq. (3.8) is accepted and the RVE
is taken unbounded, they can be found in a form, tantamount to
(3.17), containing in the integrands the multi-point averages (correlation functions) 0 (x)0 (y) 0 (w). Explicit formulae and details
are supplied by Hori et al. [HOa], [HOb]. (In [HOb], [HOc] it is
shown that the series (3.3) can be summed under certain additional
(and physically reasonable) conditions for the multipoint correlations,
corresponding to special kind of heterogeneous media.) The general
problem of convergence of the expansions of the type of (3.4) is addressed by Gambin and Kroner [GAa].
The foregoing results show that the effective conductivity , like
any other macroscopical property, depends on the whole infinite set
of all multipoint correlations, i.e. of the full information about the
internal structure of a heterogeneous medium. This very important
fact was first demonstrated and clearly stated by Brown [BRd]. The
general underlying reason was already pointed out in Section 1.2.1
the nonlinear dependence of the solution of a PDE, like (2.7), upon
its coefficients.
As a consequence, any theory that tries to predict the macroscopic properties on the base of limited information is bound to be
only an approximation. Hence, a question of central importance for
such approximate theories is the existence and the type of realistic
heterogeneous media (if any), whose effective properties coincide with
those that the theories predict.
1.3.2 The Hashin-Shtrikman Estimates
With a limited information about the internal structure of a medium,
the best that one can do is to find an interval, spanned by the plausible effective properties. The more we know, the narrower this interval should be (collapsing to the exact value, in principle, if we
know everything, i.e. the full details of the microstructure). If
the medium is two-phase and the volume fractions are solely known,
this interval is supplied by the classical Wiener bounds. If we know
morethat the medium is also macroscopically isotropicthe Wiener
bounds can be narrowed to the Hashin-Shtrikman ones. These basic
and fundamental results in the theory of heterogeneous media will be
discussed here.

1.3 Micromechanics of Heterogeneous Media

61

[Link] The Bergman Formula

For a two-phase heterogeneous solid let us, at the fixed boundary
conditions (2.3), change the given conductivity field (x) to (x) +
(x), with a small (x), as compared to (x). (Or, in other words,
let us vary (x).) Note that
(x) = 1 (x)1 + 2 (x)2 ,

(3.20)

where 1 (x) and 2 (x) are the characteristic functions of the regions
V1 and V2 , occupied by the constituents 1 and 2 respectively, and
1 , 2 are the variations of their conductivity. As a result, the temperature field will change by (x), and the effective conductivity
by . Since both (x) and (x) + (x) satisfy the same boundary
condition (2.3),


(x)  = 0 ,
(3.21)
V

which implies, in particular, that = + = G.

According to the energy definition (2.21) of the effective conductivity
1
2 (

+ )G2 = 12 ((x) + (x))|(x) + (x)|2

= 12 (x)|(x)|2 + 12 (x)|(x)|2

(3.22)

+ (x)(x) (x) ,
having neglected quantities of higher orders. The first term in the
right-hand side of Eq. (3.22) equals 12 G2 , see (2.21). Applying the
Green formula for the last term there gives
(x)(x) (x) =

1
V

(x)(x) (x) dx

1
=
V

1
V

(x) ((x)(x)) dx

(3.23)

(x)(x)
S

(x)
dS ,
n

which vanishes, in virtue of (2.7) and (3.21). Hence

G2 = (x)|(x)|2 = 1 1 (x)|(x)|2 + 2 2 (x)|(x)|2 ,

Heterogeneous Media: Modelling and Simulation

62

which means that

1

= 2 i (x)|(x)|2 ,
i
G

i = 1, 2 .

(3.24)

The important formula (3.24) was first noticed (and exploited for
bounding the effective properties) by Bergman [BEk], and that is why
we shall refer to it as Bergmans. In the foregoing reasoning we have
followed the more rigorous derivation of Bobeth and Diener [BOa],
who rediscovered it later on. One more independent derivation was
given recently by Matheron [MAi].
The formula (3.24) has a simple interpretation. Imagine that,
for a given two-phase medium, we start varying the properties i
of the constituents. As a result, the effective conductivity will become a function of i , i.e. = (1 , 2 ). The formula (3.24) then
means that the mean energy, accumulated within each one of the constituents, i.e. i (x)|(x)|2 , is proportional to the partial derivatives
/i , i = 1, 2.
[Link] The Matheron Inequalities
Since i (x), as a characteristic function, takes only the values 0 and
1, 2i (x) = i (x). Then
Z
2 Z
Z





(x)(x)
dx

(x)
dx
i (x)|(x)|2 dx ,
i
i


V

as it follows from the Schwartz inequality. Obviously

Z
V

i (x) dx = Vi .

Recalling that i = Vi /V , we get from (3.24) the inequalities

2
1 

(x)|(x)|
,
i
i
i G2

i = 1, 2 ,

(3.25)

pointed out by Matheron [MAi].

Following now the reasoning of the same Matherons paper [MAi],
note that for the two-phase medium under study one has
1 (x) =

2 (x)
,
2 1

2 (x) =

(x) 1
,
2 1

(3.26)

1.3 Micromechanics of Heterogeneous Media

63

because (x), as a binary function, takes the values i in the phase

i, i = 1, 2. Using the definition (2.10) of the effective conductivity
together with (3.26) yields
1 (x)(x) =

2
G,
2 1

2 (x)(x) =

1
G,
2 1

in virtue of Matherons inequalities (3.25). Hence

1
1

2
2 1

2

2
2

1
2 1

2

(3.27)

[Link] The Wiener Bounds

Note next that is a homogeneous function (of first degree) of 1
and 2 , i.e. (1 , 2 ) = (1 , 2 ), 0. Then

1
= 1
+ 2

1
2
1

2
2 1

2

2
+
2

1
2 1

2

The latter inequality can be simply factorized to give

(v )( r ) 0 ,

(3.28)

where v and r are the Voigt and Reuss approximations, see (2.11)
and (2.18). Since r v for all i > 0, 1 + 2 = 1, Eq. (3.28) has
as a consequence
( 1/ )1 k .
(3.29)
Hence the Reuss and Voigt values supply an interval within which
the effective conductivity of the heterogeneous solid should always
lie, independently on the details of the internal (macroscopically
isotropic) constitution. In the scalar (dielectric, more precisely) context this important fact has been first proved by Wiener [WIa], by
means of certain algebraic arguments. That is why the bounds (3.29)
are often called Wieners.
It is important to point out that the bounds (3.29) are sharp, in
the sense that there exist heterogeneous solids whose conductivity
equals the limiting values r and v . (Such solids are the simple
laminate or fiber media, already discussed in Remarks 2.1 and 2.2.)

64

Heterogeneous Media: Modelling and Simulation

Remark 3.4 Both media whose effective conductivities (along certain directions) coincide with the Voigt and Reuss values are anisotropic. If, additionally, it is required that the medium be macroscopically isotropic, then there exists an interval of admissible values of ,
narrower than (3.29). This is one of the most remarkable facts in the
theory of heterogeneous media. The narrower interval is specified by
the so-called Hashin-Shtrikmans bounds, to be derived in a moment.
[Link] The Hashin-Shtrikman Bounds
To get more restrictive bounds consider, again after Matheron [MAi],
the function
1
1
F (1 ) =

.
(3.30)

2
1 (2 1 )
Here, for the given two-phase medium, the conductivity 2 of the
phase 2 has been fixed, so that depends on 1 solely. The derivative
"

 #
dF
1

1 2 2
=

0,
d1
(2 )2 1 1 2 1
as a consequence of the first of the inequalities (3.27). Therefore
F (1 ) is an increasing function of 1 and, consequently, if 1 > 2 ,
F (1 ) F0 ,

F0 =

lim

1 2 +0

F (1 ) .

(3.31)

(If 1 < 2 , the opposite inequality holds.) The evaluation of the

limit F0 is straightforward, using the formula (3.14) and the final
result reads
B1
2
F0 = 2
or F0 =
,
(3.32)
32 1
1 2
having recalled the value of B1 from (3.14).
Remark 3.5 This is the only place in the reasoning, yielding the
Hashin-Shtrikman estimates, where the approximation (3.8) and its
consequence (3.14) are to be used, together with the value of B1 from
(3.13). If this value is left unspecified, from the considerations below
the appropriate HS bounds for a finite body can be derived as well.
However, the finite-body Green function should be then invoked.
As a result the explicit evaluation of B1 will be extremely difficult
even for the simplest regions V (say a ball), when this function is

1.3 Micromechanics of Heterogeneous Media

65

known. Moreover, the constant B1 will then feel both the details
of the internal structure (through the two-point correlation) and the
shape of the RVE, V.
Inserting F0 from the second equality of (3.32) in (3.31) gives
F (1 ) =

1
1
2

2
1 (2 1 )
32 1

which can be recast in the more compact form

1 2 []2
,
32 2 []

if 1 > 2 .

(3.33)

The opposite inequality holds, if 1 < 2 .

In a fully similar manner, the function
F (2 ) =

1
1

1 2 (2 1 )

can be considered, instead of (3.30), with 1 being fixed. Application

of the second of the inequalities (3.27) will ensure that F (2 ) is decreasing. Hence, for 1 > 2 , a lower bound on shows upit has
the same form as (3.33), but with indices 1 and 2 interchanged.
Hence, together with (3.33), we have
+

HS HS ,

HS =

1 2 []2
,
32 2 []

if 1 > 2 ,
+
HS =

1 2 []2
.
31 + 1 []

(3.34)

For 1 < 2 , the same inequalities (3.34) hold, but with reversed
lower and upper values of the bounds.
The bounds (3.34) coincide with the ones derived by Hashin and
Shtrikman [HAl] by means of an original variational procedure. A
simple check demonstrates that they are more restrictive that the
Wiener bounds, given in (3.29).
Note that the upper bound in Eq. (3.34) can be easily recast in
a form, coinciding with the Maxwell formula (1.19). This means
that, whatever the internal structure of a (macroscopically isotropic)
medium, Eq. (1.19) always provides an upper estimate on the effective
conductivity, if the matrix is stiffer than the inclusions (1 > 2 )
and a lower one in the opposite case (1 < 2 ).

66

Heterogeneous Media: Modelling and Simulation

In principle, the generalization of the HS bounds to the elastic

case can be obtained in a similar manner, though the technical details
will be much more cumbersome. The original derivation, using the
Hashin-Shtrikman variational principle [HAk], can be found in the
paper [HAm], see also [WId]. For the bulk modulus, however, much
simpler derivation can be proposed, following certain Hills arguments
[HIb], as we shall see in Section 1.3.5 below (cf. Eq. (3.68)).
1.3.3 The Hashin Assemblage
The Hashin-Shtrikman (HS) bounds (3.34), apart from their extreme
simplicity, possess also the very important properties that they, similarly to the Voigt and Reuss bounds, are sharp [HAe]. In other words,
they are the narrowest bounds on the effective conductivity that can
be constructed, provided we know that the medium is two-phase,
macroscopically isotropic, and has prescribed volume fractions 1 , 2
of the constituents.
For a demonstration that the HS bounds are sharp it suffices to
construct special two-phase media whose conductivities coincide with
+
the values
HS and HS . These are the so-called Hashins assemblages
which look as follows [HAe].
The typical element of the assemblages is a composite sphere,
denoted by Va,b . It has an external radius b and contains a concentric
spherical inclusion of the radius a, whose conductivity is, say, 2 .
The coating around the latter, i.e. the region a < r < b, has the
conductivity 1 . If 2 is the given volume fraction of the phase 2,
we require that
a3 /b3 = 2 .
(3.35)
Next we fill up the whole space (the RVE, more precisely), with
such composite spheres of different external radius b, from finite down
to infinitesimally small. The obtained two-phase medium is just the
Hashin assemblage.8 The condition (3.35) assures that the volume
fractions of the constituents equal the prescribed values 2 and 1 =
1 2 .
Since the structural elements of the assemblage are fully similar,
and differ only in a scale factor, the temperature on their surfaces
will be the same, coinciding with the one imposed on the macroscopic
8
It is noted in passing that the assemblage is an excellent example of a selfsimilar or fractal-like structure which looks one and the same in all length scales.

1.3 Micromechanics of Heterogeneous Media

67

boundary, see Eq. (2.3). Thus, to find the effective conductivity of the
composite, it suffices to specify the temperature field within a single
composite sphere, i.e. to find the continuous function (x), such that


(x) 

(x) = 0 ,


= G x,
r=b

(3.36)



2
= 1
.

r r=a0
r  r=a+0

The solution of the problem (3.36) has the form (x) = G (r),
cf. (1.9), with the potential function
(

(r) =

1
2
2 C1 r ,
1
2
2 C2 r

if r < a,
C3 a3 /r,

(3.37)

if a < r < b.

The boundary conditions (3.36) elementary specify the unknown

constants:
C1 =

(2 1 )
,
1 1

C2 =

1
,
1 1

C3 =

,
1 1

(3.38)

with given in Eq. (1.15).

The mean flux in the RVE coincides with its mean value over the
typical composite sphere Va,b :
q = a G ,

q = hqiVa,b ,

1
=
Va,b

1
q(x) dx =
V
Va,b
a,b

hqiVa,b

1
=
2
Va,b


(x)(x) dx

(x) dx + 1

r<a

(3.39)

Va,b

(x) dx ,

a<r<b

where a is the effective conductivity of the assemblage under study.

Elementary calculations, eventually give
a = 1

1 + 22
,
1 2


(3.40)

having used that (x) = G (r), with (r) already found in

Eqs. (3.37) and (3.38), the isotropy of the integrals that enter
Eq. (3.39), as well as the condition (3.35).

68

Heterogeneous Media: Modelling and Simulation

Eq. (3.40) obviously coincides with the Maxwell formula (1.19).

The latter, as already pointed out, is always one of the HS bounds
(3.34). To get the other of the bounds one should only reverse the indices 1 and 2 of the phases (i.e. the phase 1 will be now surrounded
by an appropriate concentric layer of 2, such that a3 /b3 = 1 ,
cf. (3.35)). This observation completes the proof that, as far as
the scalar conductivity is concerned, the HS bounds (3.34) are indeed optimal for the given amount of information about the medium
(macroscopically isotropic two-phase mixtures with prescribed volume fractions).
1.3.4 The Variational Estimates
Consider now the standard variational principles for the transport
problems in heterogeneous media, again in the simplest context of
scalar conductivity. Their application for bounding the effective properties is a classical one, and will be briefly recalled and summarized.
Within the frame of the volume averaging approach, pursued here,
they turn out to be of very limited significance, since the boundary
values of the admissible fields are prescribed. Instead, the minimum
principle of Willis (Theorem 2.2), free of such boundary limitations,
will to be applied now.
[Link] The Dirichlet Principle
Let (x) be the actual temperature field in a heterogeneous solid,
i.e. the solution of Eq. (2.7) with the boundary condition (2.3). Consider the class A of admissible (or trial) fields which are continuous and almost everywhere differentiable, subjected only to the same
e
e
boundary condition (2.3). Let (x)
A, then both (x)
and (x)
satisfy the same condition (2.3), so that the variation (x) should
vanish on V:



(x) 

= 0,

e
(x)
= (x) + (x) .

(3.41)

The energy of an admissible field, as defined in Eq. (2.21), is a

e
functional over the class A. Its value, for a given field (x)
A, is
e
e
f = W [(x)]
2 = 1 (x)|(x)|2
W
= 12 (x)|(x)|
2

+ (x)(x) (x) +

1
2
2 (x)|(x)| .

(3.42)

1.3 Micromechanics of Heterogeneous Media

69

The middle term in the right-hand side of Eq. (3.42) vanishes, as

a consequence of the Gauss theorem, cf. Eq. (3.23). The last term
there is non-negative and therefore
f,
W W

W = 12 G2 ,

(3.43)

having recalled Eq. (2.25) as well. Hence the actual field in the
medium minimizes the energy functional W in the class A of admissible fields.
The statement (3.43) represents the first of the standard variational principles for the BVP (2.7), (2.3). It is often called Dirichlet
principle (since, for a homogeneous medium, when (x) = const, it
states that the solution of the Dirichlet problem (x) = 0, (x) | V =
0 (x), Rwith a prescribed function 0 (x), minimizes the Dirichlet integral V |(x)|2 dx).
[Link] The Thompson Principle
Dually, let q(x) be the actual heat flux in the medium, i.e. the solution of the BVP (2.13), (2.19). Consider the class B of admissible
e (x) which are continuous and almost everywhere differflux fields q
entiable, such that
e (x) = 0 ,
q

en
q

= Q n.

(3.44)

e (x) q(x) is divergence-free and vanThen, the variation q(x) = q

ishes on the boundary V:

q(x) = 0 ,

q n 

= 0.

(3.45)

The energy of an admissible field, as defined in Eq. (2.26), is

a functional over the class B. Its value, for an admissible flux field
e (x) B, is
q
e = [q
e (x)|2
e (x)|2 = 12 k(x)|q
e (x)] = 12 k(x)|q

(3.46)

+ k(x)q(x) q(x) + 12 k(x)|q(x)|2 .

The last term in the right-hand side of Eq. (3.46) is non-negative.
The middle term there vanishes, in virtue of Eq. (3.45):
Z
1
k(x)q(x) q(x) =
(x) q(x) dx
V V
1
=
V

1
(x) q(x) dx +
V
V

(x)qn (x) dS ,
V

70

Heterogeneous Media: Modelling and Simulation

having recalled that (x) = k(x)q(x), and employing the Gauss

theorem. Therefore,
e,

= 12 k Q2 .

(3.47)

Hence the actual flux field in the medium minimizes the energy
functional in the class B of admissible fields, defined by means of
(3.44).
The statement (3.47) represents the second of the standard variational principles for the heat propagation in a heterogeneous media.
Often called Thompson principle, it specifies the extremum property
of the actual heat flux in such media.
[Link] The Beran Bounds
The principles (3.43) and (3.47) allow us to derive variational estimates on the effective properties in a standard manner, using approe
e (x) B. Then (3.43) and (3.47)
priate trial fields (x)
A and q
yield
f
Q2
2W
2 ,
(3.48)
e
G
2
since k = 1/ , see (2.30).
Like many other disciplines, where variational principles are of
importance, the art of their applications lies, in general, in the convenient and skillful choice of trial fields. This means that the values
of the energy can be comparatively easily evaluated for them and
the resulting estimates (3.48) should be as close as possible.
As a simplest application of the variational principles (3.43) and
(3.47) let us choose as trial fields, respectively, (x) = G x and
q(x) = Q. They both are admissible and reproduce the Wiener
bounds (3.29), in which only the volume fractions of the constituents
show up.
To get more restrictive bounds, broader classes of trial fields should
be invoked. For a big, but still finite RVE, a very serious problem
immediately arises in this connection, namely, the necessity that the
boundary conditions (3.41) and (3.44) be met. Then, if mathematical rigour is to be preserved, the application of the appropriate finite body Green function seems unavoidable, together with a careful
limiting passage to an unbounded volume V. That is why, strictly
speaking, the only application of the classical variational principles of

1.3 Micromechanics of Heterogeneous Media

71

Dirichlet and Thompson here, within the volume averaging approach,

is a rederivation of the Wiener bounds. Instead, the Willis minimum
principle (Theorem 2.2) should be employed, since it imposes no restrictions on the boundary values of the admissible fields.
To this end, recall the idea of Beran [BEg] to employ as trial
fields the truncated perturbation series (3.3), multiplying the terms
with adjustable scalars. (This is obviously a Ritzs type variational
procedure, with basis functions ingeniously chosen as the terms of
the series (3.3).) In the simplest nontrivial case the class
e
(x)
= G x + (1) (x)

(3.49)

of trial fields shows up, where (1) (x) is defined in Eq. (3.8), and is
a certain adjustable constant.
For the trial fields (3.49), the energy (2.21) is a quadratic function
of


e
f = W [(x)]
W
= 12 A 2B + C2 G2 ,
(3.50)
with the coefficients
A = ,

B=

1
G 0 (x) (1) (x) ,
G2

i
1 h
C = 2 | (1) (x) |2 + 0 (x) | (1) (x) |2 .
G

(3.51)

The quantity | (1) (x) |2 that appears in the coefficient C can

be evaluated, using the equation (3.5) for (1) (x). Indeed, multiply both sides of Eq. (3.5) by (1) (x) and apply the Gauss theorem;
the boundary conditions for (1) (x) assure then disappearance of the
surface integral, thus yielding
| (1) (x) |2 = 0 (x)(1) (x) G = B1 G2 .

(3.52)

For the second term in the coefficient C, we employ in a similar

manner the equations for the functions (1) (x) and (2) (x), see (3.5)
and (3.6). The final result is
0 (x) | (1) (x) |2 = G 0 (x)(2) (x) ,
and hence
C = B2 B1 ,

B = B1 ;

(3.53)

Heterogeneous Media: Modelling and Simulation

72

the second equality is obvious from Eqs. (3.10) and (3.52).

Minimizing the right-hand side of Eq. (3.50) with respect to
gives the estimate
+
B,

+
B =A

B12
B2
=
.
C
B2 B1

(3.54)

on the effective conductivity. Eq. (3.54) is a consequence, let us

point out once again, of the Willis minimum principle. That is why
this estimate, as a result, does not need the explicit form of the
perturbation coefficients B1 and B2 from the series (3.4). It holds
therefore for a finite volume V as well. (But these coefficients will be,
in general, form-dependent, since the finite-body Green function,
corresponding to the RVE, V, should be exploited in their evaluation;
moreover, even in the macro-isotropic case B1 will depend on the
whole two-point correlation function.)
For an unbounded RVE, using the approximation (3.8) and, accordingly, the formulae (3.10) and (3.17) for B1 and B2 , respectively,
the bound (3.54) takes already the explicit (form-independent, so
to say) form, due to Beran [BEg]:
(

+
B,

+
B

0 2
= 1
2
1
3

1+

0 3

!)

I3
0 2

(3.55)

In a fully similar manner the lower Beran bound can be deduced,

using the Willis principle. The final result, skipping the purely technical details, reads
(

02
2k
3

1+

1
2

k0 3
k k0 2

!)1

(1 +

I3k )

(3.56)
where r = (1/)1 is the Reuss approximation (2.18) and k(x) =
1/(x) is the resistivity field for the medium. The parameter I3k has
exactly the same form as I3 in Eq. (3.18), but with (x) replaced by
k(x).
For a two-phase medium it is easy to show that I3 = I3k , i.e. the
upper and lower Beran bounds depend on a single parameter I3 . The
latter reflects in a certain integral way the internal structure of the
medium, by means of the three-point correlation function M3 (x, y),
see (3.19). In this sense the Beran bounds are three-point bounds.

1.3 Micromechanics of Heterogeneous Media

73

It is noted that the foregoing derivation of the Hashin-Shtrikman

bounds (3.34) implies that they are two-point because only the twopoint correlation M2 (r) is needed for their evaluation. (This is the
assumption of macro-isotropy that, happily, has left only the value
M2 (0) of the latter in the final result.)
Moreover, for a weakly inhomogeneous medium the Beran bounds,
+

effective conductivity,
B and B coincide with the exact value of the


3
as given in Eq. (3.17), to the order O () . This means that

+
B B =



0 2
0 3
+ 2 I3 + O ()3 .
3
3

(3.57)

The sign indicates

here that the appropriate quantities differ in the


3
order O () . In this sense the Beran bounds are third-order ones.
The formula (3.57) well illustrates a basic and well-known feature of variational estimates (that holds not only in the context of
heterogeneous media). Namely, the more structural information is
incorporated, the narrower the bounds become. Indeed, for the simplest Wieners bounds (3.29), when the volume fractions are only
known in a binary medium
e,
v r = 1 2 []2 /
e is defined in (3.62) below. Hence v and r coincide to the
where
order O (); they are therefore first-order. The HS bounds (3.34) incorporate the additional information that the medium is statistically
isotropic. Accordingly

+
HS HS =

21 2 []3
,
e)
(21 + )(21 +

so that the coincidence is to the order O ()2 ; they are therefore

second-order. Finally, incorporating some (but not full) three-point
information in the Beran bounds (3.55), (3.56) makes them thirdorder, as already pointed out.
To conclude the story about Berans bounds, let us return to the
series (3.4). It is obviously the Taylor expansion of the effective conductivity treated as a function of the parameter = []/, see
Eq. (3.2). (That is, with the geometry of the two-phase medium being fixed, we vary the conductivity of one of the constituents, say

74

Heterogeneous Media: Modelling and Simulation

2 , having fixed 1 .) This parameter need not be small, in general, for the series to converge. Moreover, it is not even necessary
to assume it realthe conductivities can well be complex numbers
which, physically, may account for the possible dielectric losses in the
medium. Then (3.4) defines a function in the complex plane, which
possesses a number of important and beautiful properties. Such
a complex-analytical viewpoint was introduced by Bergman [BEk],
[BEl], who, having recognized these properties, employed them for
deriving bounds on the effective conductivity that hold both in real
and complex cases, see also Miltons papers [MIc], [MId]. The reasoning of Bergman and Milton was put on a rigorous base by Golden
and Papanicolaou [GOa]they proved that the function = ()
in the complex region can be represented as a certain Stieltjes integral. This representation opened the way for systematic derivation of
bounds on , using truncation of the Taylor series and the so-called
Pade approximants, whose definition and a simple application will be
very briefly discussed below. Any details of this beautiful theory are
beyond the scope of the present elementary survey (see, e.g. [MIh] for
a very clear summary of the basic ideas). But, we have accumulated
until now some basic results, happily, that will suffice for a small
illustration.
Let us truncate the Taylor series (in the point z = 0) of a function f (z) after the z p+q+1 -term. The rational function P (z)/Q(z) is
called (p, q)-Pade approximant for f (z), if its Taylor expansions coincides to the order z p+q+1 with the said truncation; here P (z) and
Q(z) are polynomials of degrees p and q, respectively, and Q(0) = 1.
(The origin, importance, theory and numerous applications of such
approximations are clarified and collected, e.g. in the book [BAa].)
To apply a Pades approximant in the study of effective conductivity , observe that we already know the Taylor expansion (3.4)


of the latter, treated as a function of , to the order O ()3 ,
cf. Eq. (3.57). The simplest nontrivial approximant will be then of
the order (1,2) and it should have the form
= + ()2

a0
,
1 + b1

(3.58)

as it follows from (3.14). The coefficients a0 and b1 are specified

by the condition that the Taylor
expansion of the function (3.58)


coincides to the order O ()3 with that in (3.57). The result, as

1.3 Micromechanics of Heterogeneous Media

75

immediately seen, reproduces the Beran bound (3.55). Hence, this

bound is nothing but a simplest Pade approximant for the function
= ().
Remark 3.6 In a purely formal manner, consider = (2 ) as a
function of the volume fraction 2 of one of the constituents.9 For a
dispersion of spheres the Taylor expansion of this function is known to
the order 2 , due to Maxwell, see (1.20). We can try to approximate
(2 ) by a simple Pades like function:
(2 ) 1

1 + a1 2
.
1 + b1 2

(3.59)

To fix the constants in (3.59), we require that it coincides with (1.20)

to the order 2 and has the value 2 , if 2 = 1. A simple check shows
that the right-hand side of Eq. (3.59) reproduces then the Maxwell
(1.19) or which is the same, one of the Hashin-Shtrikman bounds
(3.34). The author cannot provide any explanation of this fact which
could be only a fortunate coincidence. At first glance it is strange as
well that starting with the formula (1.20), valid for spherical shape
of the inclusions in a dispersion, one gets the HS-bound out of the
approximation (3.59), already valid for arbitrary (macro-isotropic)
internal geometry.

[Link] The Torquato-Milton Parameter

The Beran bounds have been considerably simplified by Torquato
[TOb] and Milton [MIc]. They introduced the so-called -parameter
which, in term of the above defined parameter I3 , can be defined
through the relation:
3(2 1 )I3 = 21 + 31 2 ,

2 = 1 1 .

(3.60)

The function = (2 ) is very often considered and approximated in the

literature on composite media, together with its expansions in powers of 2 (called
virial, borrowing the terminology from the theory of dense gases); see Section 1.5
below. Despite this, the reader should be warned that it has in general no exact
meaning, if the internal structure is not carefully predefined. The point is that
for a given volume fraction 2 there exist (infinitely) many internal constitutions;
changing 2 by an infinitesimal amount only, this constitution can be drastically
changed and rearranged.
9

Heterogeneous Media: Modelling and Simulation

76

By means of s, the Beran bounds (3.55) and (3.56) adopt the

amazingly simple and symmetric form:

2k 0 2

e + hki
2hki

0 2
,
e i + 2 hi
h

(3.61)

hki = 1 k1 + 2 k2 ,

(3.62)

with the notations

e = 1 k2 + 2 k1 ,
hki

and similarly for any other quantity instead of k, see [MIc].

The bounds (3.61) should be more restrictive than the elementary
Wieners bounds (since the appropriate functionals are minimized
over broader classes of trial fields). This implies that the denominators in the both fractions in (3.61) should be non-negative, whatever
the properties of the constituents 1 , 2 0, which immediately
yields
0 1 , 2 1 .
(3.63)
The limiting values 1 or 2 = 0, 1 reproduce one of the HS-bounds
(3.34). Together with (3.63), this fact implies that the Beran bounds
are always more restrictive than the Hashin-Shtrikmans ones (something which becomes clear only after the simplification (3.61) of the
former). This sounds natural, but by no means is obvious (since the
Beran bounds make use of the three-point information only partially,
and this does not suffice to claim in general that they are better than
any two-point ones, like Hashin-Shtrikmans).
The evaluation of the -parameter is a nontrivial and complicated
task for a given heterogeneous medium; the formula (3.18) for I3
and the definition (3.60) imply that represents in general a six-fold
integral with a singular kernel. The latter should be first drastically
simplified before trying to apply numerical or analytical methods.
Considerable efforts have been dedicated to this end, starting perhaps with the series of Corsons papers, see [COa] and the ensuing
pages of the same issue of J. Appl. Phys. There the needed parameter
( as a matter of fact) has been evaluated numerically using the outcome of extensive experimental work, concerned with measurement
of two- and three-point correlations of a real two-phase alloy (based
on a Buffons kind of needle and triangle games).

1.3 Micromechanics of Heterogeneous Media

77

For special (and perhaps idealized) models of heterogeneous media

the parameter was analytically found. For example, one of the first
and simplest models due to Miller [MIb], can be recalled here. This
is the so-called cell material, constructed through dividing the space
R3 into cells. Each cell is supplied afterward with conductivity 1
and 2 with probability 1 and 2 , respectively, independently of the
properties of the surrounding cells. For this model both the two- and
three-point correlations are quite simple and it turns out, for example,
that 1 = 1 , in the case of spherical cells, see [MIb], [HOa] et al. For
other, more realistic and complicated models of heterogeneous media,
details and extensive references, concerning the evaluation of the parameter, are given in Torquatos survey [TOc], see also the recent
paper of Jeulin [JEc].
1.3.5 The Bounds in the Elastic Case
All scalar considerations, performed until now, can be generalized
to the elastic case in an obvious manner. The derivation of the appropriate expansions like (3.17) will be considerably more tedious,
due to the tensorial nature of the basic quantities (strain and stress
tensors).
There is however another, more important, complication for an
elastic heterogeneous solid, namely, the fact that there exist at least
two small parameters, say, k and in the isotropic case, instead of
the single parameter (3.2) of the scalar case. Both these parameters
will be small for a weakly inhomogeneous solid, k,  1, but they
may be small of different orders. Details can be found, e.g., in the
papers [MOa] and [LOa]. One of the central perturbation results
of these works is Eq. (3.67) below. It is however a simple particular
case of the general formula of Hill, see (3.66).
The Hashin-Shtrikmans bounds (3.34) can be generalized to the
elastic case extending the arguments of Section 1.3.2, once the perturbation expansions of the moduli (to the second order of k and
) are known. The original derivation of these bounds employs the
appropriate variational principle [HAl], [HAm], see also [WId]. It will
be omitted here, due to the number of technical details it employs.
Instead only the bound on the effective bulk modulus will be derived
below in a different (non-variational) way, following Hill [HIb].
The classical variational principles, based on minimization of the
energy of a RVE, V, are formulated in a straightforward manner for

Heterogeneous Media: Modelling and Simulation

78

an elastic medium. Following Hill [HIa] (who dealt with polycrystals,

to be precise), the simplest trial fields of homogeneous strain and
stress within the volume V produce the counterpart of the Wieners
bounds (3.29), namely:

Lr L Lv ,
Lv = L(x) ,

(3.64)

Lr = ( L1 )1 ,

where Lv and Lr are the Voigt and Reuss approximations for the
tensors of the effective moduli, respectively. That is why the bounds
(3.64) are often called Hills. For elastic (two-phase) composites their
derivation is due also to Paul [PAc], so that sometimes they are referred to as Pauls bounds as well.
The generalization of the Berans bounds to the elastic case is, in
principle, straightforward, though again considerably more cumbersome. For the effective bulk modulus k , when the applied macrostrain is spherical, the appropriate bounds are due to Beran and
Molyneux [BEj]. They contain the same statistical parameter , that
appeared in the scalar bounds (3.61), but this became clear only when
Milton [MIc] drastically simplified the original Beran-Molyneux estimates, recasting them as

1/k

4(1/k)0 2

g + 3 h1/i
4h1/ki

k k

3k 0 2
e + 4 hi
3hki

(3.65)

Obviously, (3.65) resembles very closely (3.61). (It is recalled that

for any binary quantity that takes the values, say, a1 and a2 in the
constituents, a0 2 = (a a)2 = 1 2 (a1 a2 )2 , cf. (3.13).)
Remark 3.7 We shall call Hills medium an elastic binary mixture
of isotropic constituents with the same shear moduli, 1 = 2 = ,
which is macroscopically isotropic. For such a medium hi = ,
h1/i = 1/ so that the parameter disappears from the bounds
(3.65). Moreover, a simple check shows that the bounds coincide,
yielding the exact value of the effective bulk modulus
k = k

k0 2
k + 43 + k 0 3 / k 0 2

(3.66)

1.3 Micromechanics of Heterogeneous Media

79

This is a remarkable result due to the fact that it holds independently

on the internal structure of the medium, under the only assumptions
of macroscopical isotropy and equal shear moduli.
The formula (3.66) is due to Hill [HIb]. His original reasoning,
which does not make use of any variational arguments, is based on
an ingenious guess for the real displacement field in the medium.
This field turns out to be a linear combination of a pure dilatation
and of the gravitational force generated by uniform distribution of
mass within one of the constituents; the dilatation it generates is
piece-wise constant.10
Remark 3.8 In the weakly inhomogeneous Hills medium (1 = 2 )
the perturbation expansion of the type (3.14) reads
k = k



k0 2
2
+
O (k)
,
k + 43

k =

[k]
,
k

(3.67)

and it immediately follows from the general formula (3.66). Moreover

Eq. (3.66) specifies easily the perturbation expansion to any order of
k. The particular case (3.67) was found independently in [MOa],
[LOa] directly, though the needed calculations are quite tedious.
For the effective shear modulus the evaluation of the Beran
type bounds is due to McCoy [MCa] and it is much more involved
than the ones for the conductivity and bulk modulus. Moreover, the
bounds include a second statistical parameter of the kind of (the
so-called -parameter of Torquato and Milton, which shows up as a
consequence of the tensorial structure of the Green function in the
elastic case). A much simpler form, similar to (3.61) and (3.65), of
the McCoy bounds on can be found again in [MIc].
It is to be noted finally in this subsection that the Hill formula
(3.66) allows to obtain easily the Hashin-Shtrikman bounds on the
bulk modulus of a binary mixture. Indeed, following Hill [HIb], con10
Strictly speaking, the Hill analysis only proves that the right-hand side of
Eq. (3.66) is an upper bound on k . The reason is that the displacement field,
devised by him, does not satisfy the homogeneous boundary condition (2.55)
exactly, and hence this can be only one of the trial displacement fields in the
minimum principle of Willis (Theorem 2.2). The mean dilatation is prescribed, but
nothing is known about the exact values of the displacement on the boundary V.

Heterogeneous Media: Modelling and Simulation

80

sider a more general case of an elastic mixture, when 2 < 1 and

k2 < k1 . (Note that this is a particular case of the so-called wellordered mixture, when (2 1 )(k2 k1 ) > 0; the badly-ordered
one corresponds to the inequality (2 1 )(k2 k1 ) < 0.) Keeping
the modulus k1 fixed, let us start increasing mentally 2 until it coincides with 1 . In this process the effective bulk modulus will increase
as well and its value will be the one, given in Eq. (3.66) when 2
reaches 1 . Hence (3.66) will provide an upper bound on k in
the case under study. In a similar way we can decrease 1 until it
coincides with 2 ; then (3.66) will give a lower bound on k , i.e.
k

k0 2
k + 43 2 + k 0 3 / k 0 2

k k

k0 2
k + 43 1 + k 0 3 / k 0 2 .

(3.68)

If 2 > 1 and k2 > k1 , the inequality signs are to be changed in

(3.68). The same inequalities (3.68) were derived by Hashin and
Shtrikman [HAm] by means of variational arguments; these authors,
however, were able to find estimates of a similar kind for the effective
shear modulus as well.
1.3.6 Cross-Properties Relations: The Levin Formula
Consider again a two-phase medium. Geometrically, this is simply a
division (in general, extremely complicated of course) of a region V
into two parts, V1 and V2 , that represent the subregions occupied by
the phases, V = V1 V2 . Such a division, purely geometrical in this
moment, can be materialized in many different ways and contexts,
in the sense that we can fill up the phases with materials of different
conductivities, elasticity constants, absorption abilities, etc. We can
as well let a fluid penetrate the volume through one of the phases,
assuming the other to be a rigid or deformable skeleton and so on
and so forth. In any of the cases, pursuing the appropriate homogenization problem, we can evaluate, in principle, the specific effective
properties. The fact that behind any one of these problems, a fixed
and common for all of them geometrical structure is hidden (the
given division of the volume V into two subregions V1 and V2 ) suggests that the effective properties, no matter how different in context
and interpretation, should be interconnected. Such interconnections
are often called cross-properties relations.
The simplest cross-properties relations appear due to an obvious
mathematical analogy between the equations that govern different

1.3 Micromechanics of Heterogeneous Media

81

phenomena. The already mentioned examples of dielectric permeability, magnetic permittivity, diffusion, heat conduction, etc., provide such examples of full formal analogy; that is why it suffices to
get results for one of these phenomena, for the rest only the interpretation should be changed. An analogy of this kind will be used
in Section 1.5.2, when deriving the Einstein formula (1.21) from a
study of an elastic dispersion of spheres. In some cases, however, the
mathematical analogy, though present, is not that obvious and some
efforts are needed to clarify it. An example is provided by the the heat
propagation in a fiber-reinforced medium, transverse to the axes of
fibers (assumed parallel cylinders) and the shear straining along the
fibers. As shown by Hashin [HAh, Chapter 5.5.1], the appropriate
effective parameters coincide, due to the same governing equations
for these two phenomena.
Much more nontrivial cross-properties relations are however those
where there is no direct mathematical analogy, and the relations result from the same underlying internal structure of the medium. Here
we shall present only the first and maybe the most famous relation of
this kindthe so-called Levin formula which ties the effective thermal expansion coefficient with the bulk modulus of a two-phase solid
[LEa].
Following Levin himself [LEa], consider two types of straining of
a thermoelastic solid: the first, due to given surface traction at fixed
temperature, the otherdue to temperature change with no such
tractions.
More specifically, let (x) and (x) be the strain and stress fields
that appear in the volume V at fixed temperature T = 0 and under
the homogeneous boundary conditions (2.58). These conditions imply
= 0 ,

= M : 0 ,

(3.69)

cf. Lemma 2.7, where M is the effective compliance tensor of the

medium.
Consider a second pair of strain and stress fields in the medium,
(t)
(x) and (t) (x), due solely to the change T of the temperature,
with no boundary tractions; let u(t) (x) be the appropriate displacement field. These fields appear as a consequence of the inhomogeneity of the thermal expansion coefficient (x) which takes the different
values, 1 and 2 , in the phases 1 and 2 respectively. According

Heterogeneous Media: Modelling and Simulation

82

to the basic thermoelastic law, we have

(t) (x) = M(x) : (t) (x) + (x)T I .

(3.70)

Moreover,
1
=
(t) (x) dx = T I ,
(3.71)
V V
which is just the definition of the effective thermal expansion coefficient of the medium (having assumed, for simplicitys sake,
isotropy of thermal expansion effects as well).
Note first the identity
Z

(t)

1
V

1
(x) : (x) dx =
V
V

(t)

1
=
V

(t)

(x) u(t) (x) dx

(3.72)
0

n (x) u (x) dS = :

(t)

= T Tr ,

which easily follows from Gauss theorem, the boundary conditions

(2.55), the self-equilibrium of the field (x) and the definition (3.71)
of the effective thermal expansion constant.
The second identity that we shall need reads
Z

(t) (x) : (x) dx =

Z
V

(t) (x) : M(x) : (x) dx = 0 ,

(3.73)

since the field (t) (x) is temperature-induced, with no surface tractions involved, cf. the derivation of (3.72).
Introduce now the Hooke law (3.70) into the first integral of (3.72)
1
V

Z
V

(x) : M(x) : (t) (x) + (x)T I

1
= T
V

dx
(3.74)

(x) Tr (x) dx = T Tr ,
V

having used (3.73). Since (x) is step-constant, Eq. (3.74) gives

1 1 hTr i1 + 2 2 hTr i2 = Tr 0 .

(3.75)

It remains to combine now (3.75) with the obvious formulae

Tr 0 = 1 hTr i1 + 2 hTr i2 ,
Tr 0 = 1 hTr i1 + 2 hTr i2 ,
Tr 0 = 3k Tr 0 ,

hTr ii = 3ki hTr ii ,

i = 1, 2 ,

1.3 Micromechanics of Heterogeneous Media

83

cf. Eq. (2.35), which allow to exclude all the traces of the stress and
strain tensors. The final result reads
=

2 (1/k 1/k1 ) 1 (1/k 1/k2 )

1/k2 1/k1

(3.76)

and this is just the Levin formula.

Remark 3.9 Almost at the same time Levins formula was rediscovered and generalized by Rosen [ROb], Shapery [SHa] et al., see also
[ROc], [KRb] and the book [CHd, Chapters 9.3 and 9.4]. The general
philosophy of deriving various cross-properties relations was recently developed and implemented by Grabovsky and Milton [GRa],
[MIf]. Torquato and co-authors have discovered another class of very
interesting and unexpected such relations, involving, e.g. fluid permittivity of a medium; details and references can be found again in
Torquatos survey in the present volume (Chapter 2).
1.3.7 The Elementary Bounds on the Effective Absorption
Coefficient
The variational principles of classical type for the absorption problem (2.66) can be easily formulated, together with their modification
[TAb]. Other more refined variational procedures in the hard limit
(2.69), that corresponds to the Smoluchowski problem, are presented
in Torquatos survey in this volume. Here we shall only derive the
elementary estimates of Wieners type on the effective rate constant
k 2 without engaging directly variational arguments. The reasoning
that follows is a reminiscent of a Berans one for the conductivity
case [BEh].
Let us multiply both sides of Eq. (2.66) by c(x) and average the
result over the volume V:
1
Kc =
V

Z h

|c(x)|2 + k 2 (x)c2 (x) dx ,

having used that

Z
V

c c dx =

Z
V

|c|2 dx ,

(3.77)

Heterogeneous Media: Modelling and Simulation

84

in virtue of the boundary condition (2.67). Let c(x) = c + c0 (x), so

that c0 (x) is the fluctuating part of c(x), c0 = 0. Then
1
Kc =
V

Z n

|c0 (x)|2 + k 2 (x) c2 + 2cc0 (x) + |c0 (x)|2

io

dx

k 2 c2

2c
+
V

k 2 (x)[c(x) c] dx = k 2 c2 + 2c k 2 (x)c(x)

i.e. k 2 c K = k 2 c. In virtue of Eq. (2.70), we thus have

k 2 c K = k 2 c ,
kv2

k2

i.e. k 2 kv2 ,

1 k12

(3.78)

2 k22 .

Fully similar manipulations, omitted here, lead to the conclusion

that the Reuss type approximation kr2 , defined as
1
1 2
= 2+ 2,
2
kr
k1
k2

(3.79)

provides a lower bound on the effective absorption coefficient. Together with Eq. (3.77), this means that
kr2 k 2 kv2 ,

(3.80)

whatever the (macroscopically isotropic) internal constitution of the

medium. This conclusion is obviously the exact absorption counterpart of the Wiener (3.29) or Hills (3.64) bounds in the scalar and
elastic contexts, respectively.
Similarly to these bounds, the estimates (3.80) are also sharp, in
the sense that they cannot be improved (see Section 1.5.7). However,
the particular constitution which realize the bounds (more precisely,
approximate them as closely as we wish), are here isotropic. Therefore, the fact that a binary mixture is macro-isotropic, does not allow
us to narrow (3.80) for given volume fractions, unlike the scalar conductivity case treated in Section 1.3.4. Bounds of Hashin-Shtrikmans
type here also exists but the full two-point correlation function enters them through a certain integral parameter, see [TAb] and also
the recent authors paper [MAg] for details.

1.4 Micromechanics of Heterogeneous Media

85

1.4 The Single Inclusion Problem

In this Section the so-called single inclusion problem is treated. It
concerns determination of temperature and/or displacement fields in
a solid, containing an inhomogeneity whose thermal and/or elastic
properties differ from those of the surrounding medium. Having in
mind the application to heterogeneous solids, the main interest for
us is the temperature and/or displacement gradient within the inhomogeneity. This gradient turns out to be constant for ellipsoidal
shape. The most important particular casessphere, fiber, layer or
diskare elaborated as examples and explicit results are listed first
in conductivity and then in elastic contexts.
The application of the obtained results will be postponed until
Section 1.5. There they will serve as a basic building block when
deriving approximate formulae for the effective properties of heterogeneous media.
1.4.1 Scalar ConductivityIntegral Equation
Consider an infinite body (matrix), with the conductivity tensor K 1 ,
containing an inhomogeneity W (a filler particle), with a different
conductivity tensor K 2 . For the moment both tensors K 1 and K 2 are
symmetric second-rank tensors, though the main attention in what
follows will be concentrated on the isotropic case when K i = i I,
i = 1, 2. The conductivity field of the inhomogeneous solid under
study is
K(x) = K 1 + [K]hw (x) ,

[K] = K 2 K 1 ,

(4.1)

where hw (x) is the characteristic function of the inhomogeneity W.

At the absence of body sources, the temperature field (x) in the
solid obeys the equation
{K 1 (x)} + {hw (x)[K] (x)} = 0 ,
as it follows from Eqs. (1.5) and (4.1).
Let G(x) be the Green function for the matrix, i.e.
{K 1 G(x)} + (x) = 0 ,

(4.2)

Heterogeneous Media: Modelling and Simulation

86

which represents the temperature field, generated by a (unit) point

source, located at the origin x = 0. The Green function allows us to
recast Eq. (4.2) as
Z

(x) = (x) + G(x x0 ) hw (x0 )[K] (x0 ) dx0

(4.3)

(the integral is over the whole R3 ), or

(x) = (x) +

(x x0 ) [K] (x0 ) dx0 ,

(4.4)

where
(x) = G(x) .

(4.5)

Eq. (4.4) is the basic integral equation which governs the temperature gradient in the solid with a single inhomogeneity, provided the
temperature gradient at infinity is prescribed. In the elasticity context it has been derived by a number of authors, see, e.g. [KUb]. Note,
however, that the kernel (x) has the singularity |x|3 at |x| 0,
so that the integration in Eq. (4.4) should be understood in the sense
of generalized functions, see [KUa].
1.4.2 Scalar ConductivityEllipsoidal Inhomogeneity
As seen from the analysis of Section 1.2.6 and the definition of the
concentration factors, the main interest when modelling effective conductivity of heterogeneous media, is the temperature gradient within
the inhomogeneity W. The latter can be easily found from the integral equation (4.4) for ellipsoidal inhomogeneities directly, without
using the elegant but longer considerations of Eshelby [ESa], concerned with the eigenstrains and eigenstresses in an inclusion. The
reasoning to be used follows the basic idea, sketched very concisely
by Wu [WUa]. (More details and many particular cases have been
explicitly treated by the author later on [MAd].)
Let G = (x) be a constant vector; it is easily seen that a
simple solution of Eq. (4.4) within W can be found, if the temperature gradient (x) is also constant in W. To simplify the technical
details, assume that the matrix is isotropic. Then
G(x) =

1
.
41 |x|

(4.6)

1.4 Micromechanics of Heterogeneous Media

87

The field (x) will be constant in W if and only if

P (x) =

(x x0 ) dx0 =

1
w (x)
1

(4.7)

is constant within the inhomogeneity W. Here

Z

w (x) =

dx0
4|x x0 |

(4.8)

denotes the Newtonian potential for the region W. Recall that the
latter solves the equation
w (x) + hw (x) = 0 ,

(4.9)

which will be used repeatedly in what follows.

Thus, to have a simple solution of Eq. (4.4), the potential w (x)
should be a quadratic function within W, which is true if W is ellipsoidal [ESa], [MUa]. Therefore, if the temperature gradient is constant at infinity and the matrix is isotropic, this gradient is also constant within any ellipsoidal inhomogeneity, whatever the anisotropy
(i.e. the tensor K 2 ) of the latter.11 That is why, Eq. (4.4) implies
(x) = h(x)iw =

1
W

Z
W

(x0 ) dx0 = Aw (1 , K 2 ) G

1

Aw (1 , K 2 ) = I + P [K]
[K] = K 2 1 I ,

(4.10)

xW.

The tensor Aw , that appeared in Eq. (4.10), plays a central role in

the elementary models of composite media, as it will be seen below.
That is why we shall specify it now for the ellipsoidal shape of W
and its particular casessphere, fiber, layer or disk.
11

This fact holds true for an anisotropic matrix as well. We shall not consider
however such a matrix, since the technical details are cumbersome, even in the
scalar context. Moreover, both in scalar conductivity and elasticity contexts,
the polynomial conservation property holds. This means that if the applied
temperature gradient or strain field at infinity is a polynomial of a degree n,
these fields are also polynomials, of the same degree n, within the ellipsoidal
inhomogeneity; see [KUb] or [KUa, Section 4.9].

Heterogeneous Media: Modelling and Simulation

88

Let W be a sphere. The potential w (x) is spherically symmetric then, so that w (x) = I is an isotropic tensor, x W.
Making a contraction and invoking Eq. (4.9), one finds = 1/3, i.e.
w (x) = 13 I ,

P =

1
I,
31

and thus

xW,

(4.11)

1

1
Aw (1 , K 2 ) = I +
[K]
31


(4.12)

If the inhomogeneity is also isotropic, then K 2 = 2 I and Eq. (4.12)

yields the already known result (1.14), (1.15),
Aw (1 , 2 ) = I ,

31
2 + 21

(sphere) ,

(4.13)

found directly in Section 1.1.4.

Let the inhomogeneity be a circular fiber along the axis x3 .
The potential w (x) then depends on the x1 , x2 coordinates only
and the same symmetry arguments, as those utilized above for the
sphere, give
w (x) = 12 (e1 e1 + e2 e2 ) ,

P = 21 1 (e1 e1 + e2 e2 ) ,

(4.14)

x W. Hence, for an isotropic fiber,

Aw =

21
(e1 e1 + e2 e2 ) + e3 e3
1 + 2

(fiber) .

(4.15)

Let the inhomogeneity be a layer h x3 h, perpendicular

to the x3 -axis. The potential w (x) then depends on x3 solely, so
that ,ij = 3i 3j , i.e.
w (x) = e3 e3 ,

P =

1
e3 e3 ,
1

xW,

(4.16)

and hence, for an isotropic layer,

Aw = e1 e1 + e2 e2 +

1
e3 e3
2

(layer) .

(4.17)

1.4 Micromechanics of Heterogeneous Media

89

Remark 4.1 The tensors Aw for the simplest shapes, treated until
now (sphere, fiber and layer), can be extracted as particular cases
of the more general result below, valid for an arbitrary spheroid W.
However, the evaluation of Aw in these cases is extremely simple,
utilizing obvious symmetry arguments for the potentials gradients
in the appropriate regions; the use of the classical, but not that elementary formula (4.20) below, is totally avoided in this way. The
same remark holds true for the elastic case, as we shall see in the
next subsection.
Let the inhomogeneity W be an ellipsoid with the semiaxes a,
b and c along the Cartesian axes x1 , x2 and x3 , respectively. Then
w (x) = (M1 e1 e1 + M2 e2 e2 + M3 e3 e3 ) ,

xW,

(4.18)

where ei are the unit vectors along the axes xi , i = 1, 2, 3.

The coefficients Mi in the latter formula are often called in the
physical literature depolarizations (the quantities 4Mi in magnetic
theory are referred to as demagnetizing factors). They have been
extensively tabulated, see, e.g. [BRe, p. 28]. Obviously
M1 + M2 + M3 = 1 ,
as it follows from Eq. (4.9). The factors Mi are expressed as the
classical elliptical integrals, cited here for the sake of completeness:
Z

M1 =

1
2 abc

(a2

du
,
+ u)(b2 + u)(c2 + u)

(a2

+ u)

(4.19)

and similarly for M2 and M3 with a cyclic change 1 2 3 1. A

clear and concise derivation of (4.19), using the appropriate ellipsoidal
coordinate system, can be found in the opening pages of Landau and
Lifshitzs book [LAa], see also the book of Stratton [STa, pp. 207215,
257258].
For our purposes it suffices however to consider only spheroids,
i.e. rotational ellipsoids with the semiaxes a = b 6= c. Then

M = M1 = M2 = 21 a2 c
Z

M = M3 = 12 a2 c

(a2

(a2

du

,
+ u)2 c2 + u

du
.
+ u)(c2 + u)3/2

Heterogeneous Media: Modelling and Simulation

90

Both integrals can be elementary evaluated:

!

1 e2
1

arcsin e , if c < a (oblate) ,

e2 1
e
M=



1 e2
1+e

ln

2e
,
if c > a (prolate) ,

2e3
1e
(

M =

1
2 (1

M) ,

c/a,

if c < a (oblate) ,

a/c,

if c > a (prolate) ,

(4.20)

e = 1 2 is the eccentricity of the ellipsoid, see again the book

[LAa]. The type of the spheroid is indicated as well in Eq. (4.20).
If both matrix and inhomogeneity are isotropic, Eqs. (4.10) and
(4.20) yield
1

Aw (1 , 2 ) = I + P [K]
h

= I+

1
[]
w
1

i1
[] 
M (e1 e1 + e2 e2 ) + M e3 e3
1



= I

[]
[]
M (e1 e1 + e2 e2 ) + 1 +
M e3 e3
1
1


1+

since I = e1 e1 + e2 e2 + e3 e3 . Thus
1
Aw (1 , 2 ) =
(e1 e1 + e2 e2 )
1 + []M
+

1
e3 e3
1 + []M

1

(4.21)

(spheroid) .

Remark 4.2 The above considered, in an elementary way, cases of

a sphere, fiber and layer, can be easily extracted from the general
formula (4.21). Indeed, in these cases = c/a = 1, and
0, respectively, i.e. the spheroid is a sphere or degenerates into
a fiber or a flat disk (layer). The depolarization factors M and M ,
defined in Eq. (4.20), then become
M = M =

1
3

(sphere) ,

M = 21 ,

M =0

(fiber) ,

M = 0 ,

M =1

(disk) .

(4.22)

1.4 Micromechanics of Heterogeneous Media

91

When inserting these values in (4.21), the particular cases (4.13),

(4.17), (4.15) will immediately appear one after another.

Let the spheroid be randomly (and uniformly) oriented, then

hAw i = I ,
=

1
3

21
1
+
1 + []M 1 + []M

(4.23)

where hAw i denotes the value of the tensor Aw averaged with respect to all possible orientations of the axes xi ; in deriving Eq. (4.23)
we have used that the tensor hAw i should obviously be isotropic.
Note that the above mentioned limiting cases = c/a and
0, when combined with the previous formula (4.23), can be
interpreted as the situations, in which the inclusion is either a randomly oriented needle or a randomly oriented disk, respectively. Using Eq. (4.22) in (4.23) yields the simple formulae in these two cases
hAw i =

51 + 2
I
3(1 + 2 )

1 + 22
hAw i =
I
2

(needle) ,
(4.24)
(disk) .

Remark 4.3 The foregoing analysis concentrated on the temperature gradient within the ellipsoidal inclusion, i.e. on the tensor Aw ,
defined in Eq. (4.10). However, once the tensor Aw is known the temperature field outside the inclusion can be immediately found, using
the integral equation (4.3). Indeed, for an isotropic matrix, with the
Green function given in (4.6), this equation yields
(x) = G x +

1
w (x) [K] Aw G
1

(4.25)

after an obvious integration by parts. If the gradients of both sides

of (4.25) are taken, the result will reproduce the formula (4.10) for
the tensor Aw , provided x W.

92

Heterogeneous Media: Modelling and Simulation

1.4.3 Singular Inclusions and Nonideal Contacts

The assumption of ideal contact between the phases, adopted until
now, is a starting point in the great majority of works on heterogeneous media. Under this idealization, the temperature field and the
normal component of the heat flux are to be continuous at phase interfaces. Imperfect interfaces are however a fact in many circumstances
and reflect themselves either in a discontinuity of the temperature
field or of the normal component of the heat flux. More details can
be found in the papers, cited below and in the references they contain.
[Link] Singular Inclusions
A convenient way to treat nonideal contacts is based on the notion
of singular inclusions, as introduced by Kanaun [KAc], see also the
earlier paper [SOa]. Though in [KAc] the elasticity context was used,
the underlying idea is simple, appealing and general.12
Let one of the inclusions dimensions be of the order , small
as compared to the other two. Let, at the same time, its material
property like conductivity s or resistivity ks = 1/s be of the order
1/. Then the limit 0 produces the two models of a singular
inclusion.
In the first model s , being of the order 1/, as 0. Then
the limit
1
CI =
lim s
(4.26)
a1 0

remains finite. Here, to make C I dimensionless, the conductivity 1

of the matrix is employed, together with a certain characteristic size,
a, of the inclusion under study. For simplicity sake, we assume that
C I = const, though in general C I can be an arbitrary function, defined
on the surface of the inclusion . (To this end it suffices to take the
inclusion inhomogeneous, with conductivity s = K(x)/ + O(1) as
0, where K(x) is a prescribed function on this surface.)
12

In a particular case, the same idea was formulated clearly by Pham Huy
and Sanchez-Palencia [PHa]. The authors treated conductivity problems in the
presence of a superconducting singular inclusion (in the terminology introduced
below, see (4.26)) and proved the appropriate existence theorems using functional
analysis arguments. Note that one of the motivation of [PHa] came, in particular,
from a diffraction problem, where the electromagnetic wave field is disturbed by
the presence of a very thin screen, possessing very high electric conductivity.

1.4 Micromechanics of Heterogeneous Media

93

The so-obtained singular inclusion can be called superconducting. The temperature remains continuous everywhere, but the heat
flux suffers a jump when crossing its surface . We shall illustrate
these facts below, when treating nonideal contact of a spherical inhomogeneity in a matrix.
In the elastic context the superconducting inclusion is somewhat
similar to a rigid flake (membrane) in the medium. The displacement is continuous, but the normal stress component has a jump on
the flake. An integral equation on the surface , which specifies the
stress field in an unbounded elastic medium with such an inclusion,
has been derived by Kanaun [KAc].
In the second model the resistivity ks = 1/s , being of the
order 1/, as 0. Then the limit
C II =

1
lim ks
ak1 0

(4.27)

ks

remains finite. Here, to make C II dimensionless, the compliance k1

of the matrix is employed; again a is a certain characteristic size
of the inclusion. Once more we have assumed, for simplicity, that
C II = const, though in general it can be an arbitrary function, defined
on the surface of the inclusion . (See the comment after Eq. (4.26).)
The so-obtained singular inclusion can be called resisting. Here
the heat flux remains continuous, but the temperature suffers a jump
when crossing the surface . Again, we shall illustrate these facts
below, when treating the appropriate nonideal contact of a spherical
inhomogeneity in a matrix.
In the elastic context the resisting inclusion is somewhat similar to
a crack. (The reader should be warned that the crack itself, as treated
in solid mechanics, is only a special particular case of the inclusion
under study, see [KAc].) The displacement has a jump, but the normal stress component is continuous on . An integral equation on
the surface , which specifies the strain field in an unbounded elastic
medium with such an inclusion, has been also derived by Kanaun
[KAc].
[Link] Nonideal Contacts
After Kanaun et al. [KAd], [KAe], we shall treat nonideal contacts as
singular inclusions, spread over the surface of a bulk inhomogeneity.

94

Heterogeneous Media: Modelling and Simulation

For an illustration consider, following Torquato and Rintoul [TOd],

the simplest case of a spherical inhomogeneity of radius a and conductivity 2 . It is coated with a thin concentric layer with thickness
and conductivity s , such that  a, s  1 . The so-obtained
composite sphere is then immersed into a homogeneous matrix of
conductivity 1 . The distinguished limit is just (4.26), with a finite
C = C I . This is the definition of the superconducting interphase,
as adopted in [TOd]. Obviously, such an interphase is nothing but a
superconducting singular inclusion, spread over the surface r = a of
the spherical inhomogeneity.
In the second case let the conductivity of the thin layer s be
vanishingly small. The distinguished limit is now (4.27), which we
shall write, after Torquato and Rintoul [TOd], in the equivalent form
R=

2
lim
,
a 0 s

(4.28)

s 0

with a finite R. This is the definition of the resistive interphase

[TOd]. Again it is clear that such an interphase is nothing but a
resisting singular inclusion, spread over the surface r = a of the
spherical inhomogeneity.
[Link] The Superconducting Spherical Inhomogeneity
Consider a spherical inhomogeneity with a superconducting interphase. At infinity, the same condition as in Eq. (1.7) applies,
i.e. (x) G x as |x| .
The solution of such a single inclusion problem is very simple
due to the assumed spherical shape. (For the ellipsoidal shape the
problem is much more complicated and the full analytical solution
has been only recently given by Miloh and Benveniste [MIa].) Indeed, consider, before passing to the limit (4.26), the above composite
sphere. The temperature field has the form (1.9), with

(r) =

1
2

2 C1 r ,

1
2
2 C2 r

1 2
2r

if r < a ,

C3

C4

a3
r

a3
r
,

, if a < r < b ,
if b < r < ,

(4.29)

1.4 Micromechanics of Heterogeneous Media

95

where b = a + . The continuity of (x) and d/dr at r = a and

r = b implies the following system for the unknowns C1 to C4 :
C1 = C2 + C3 ,
 3

a
b

C2 +

 3

C3 = 1 +

a
b

C4 ,
(4.30)

2 C1 = s (C2 2C3 ) ,
!

 3

a
C2 2
b

 3

C3

a
12
b

= 1

C4

cf. Eqs. (1.12) and (1.13).

In the superconducting case under discussion, the solution of (4.28)
has the form
Ck =

Ck0

 

+O

k = 1, . . . , 4 .

In the limit s , the first three of the equations (4.30) give

C20 + C30 = 1 + C40 ,

C20 = 2C30 ,

C10 = 3C30 .

(4.31)

Note next that

 3

a
b

+O
a
a

 

=13

(4.32)

(since = b a) which, when employed in the last of the equations

(4.30), yields
2 C10 + 61 CC30 = 1 (1 2C40 ) ,
(4.33)
having taken the limit (4.26), with (4.31) taken into account, and
recalling that C = C I .
The relations (4.31) and (4.33) specify the needed constants C10
to C40 in the superconducting limit:
C10 =
C30

= C1 /3 ,

31
,
2 + 2(1 + C)1
C40

C10

C20 = 2C1 /3 ,

1 (1 2C) 2
1=
.
2 + 2(1 + C)1

(4.34)

Heterogeneous Media: Modelling and Simulation

96

Hence, in the conductance case under study, the solution of the

spherical inclusion problem has the same form (1.14) as that for the
ideal contact, cf. Eq. (1.14), namely,
(

(x) = G x

s ,
1

if |x| < a ,
s (a/r)3 ,

if |x| > a ,

(4.35)

with the constants

s = C10 = 1 s =

31
,
2 + 2(1 + C)1

s = C40 .

(4.36)

The temperature field (4.35) is obviously continuous at the inclusion

surface r = a, but the heat flux suffers a jump in this case


[qn ] 

r=a



[qn ] 

6= 0 ,



= 1
2
.

r=a
r r=a+0
r  r=a0

Its value can be easily found, using Eqs. (4.35) and (4.36):



[qn ] 

r=a

= 1 (1 + 2 s ) 2 s G er
621 C
=
G er .
2 + 2(1 + C)1

(4.37)

The surface concentration factor A(12)

s , for a single sphere (we
underline this adding the subscript s), can be now easily evaluated:
q

(12)

1
=
Va

Z
r=a

A(12)
s (1 , 2 )

[qn ]x dS = A(12)
s (1 , 2 ) G ,

621 C
=
I,
2 + 2(1 + C)1

(4.38)

cf. Eq. (2.49).

[Link] The Resistive Spherical Inhomogeneity
The solution of the spherical inhomogeneity problem with a resisting interphase is fully similar to the superconducting considerations, presented in the previous subsection.
Indeed, the temperature field in this case, before passing to the
limit (4.28), has the same form (1.9), (4.29), so that the system (4.30)

1.4 Micromechanics of Heterogeneous Media

97

for the Ci s remains the same. When /a 0, the unknowns can be

expanded as
C1 =

C10

+O
a

 

C40

C4 =

C2 = D + C20 + O

 

+O

,
(4.39)

C3 = D + C30 + O

 

 

retaining only the terms that matter in the limit (4.28). Taking this
limit yields
C10 =
2 C10

3D
,
R

C10 + 3D = 1 + C40 ,

= 1 (1

2C40 ) ,

C10

C20

(4.40)
C30 ,

having recalled Eq. (4.32). The constants C20 and C30 remain undefined; they however describe the solution within the disappearingly thin coating. Moreover, only their sum, i.e. C10 , matters, see
Eq. (4.40). Solving the system (4.40) is elementary: the single sphere
field in the resistance case has again the familiar form (1.14), namely,
(

(x) = G x

r ,
1

if |x| < a ,
r (a/r)3 ,

if |x| > a ,

(4.41)

with the constants

31
,
2 + 2(1 + R)1

r = C10 =
r

C40

2 (1 + R)1
=
.
2 + 2(1 + R)1

(4.42)

The temperature field (4.41) jumps at the inclusion surface r = a:




[] 

r=a

= (1) (2) 
= RC10 G x =

r=a

3R1
G x.
2 + 2(1 + R)1

(4.43)

Hence the interphase concentration factor B (12)

s , for a single sphere

Heterogeneous Media: Modelling and Simulation

98

(we again underline this adding the subscript s), is

(12)

1
Va

B (12)
s (1 , 2 )

r=a

[]n dS = B (12)
s (1 , 2 ) G ,

3R1
=
I,
2 + 2(1 + R)1

(4.44)

cf. Eq. (2.52).

Note that at r = a the heat flux, being continuous, is proportional
to this jump:
1


(1)
(2)
2  (1)
= 2
=
(2) at r = a .
r
r
aR

(4.45)

The boundary condition (4.45) corresponds to the so-called skin effect, and the proportionality coefficient in it, 2 /(aR), is the skin constant.13 The skin effect reflects the phenomenon, known as Kapitzas
resistance of an interphase boundary, see, for example, the papers
[BEb], [BEf], [TOd], [MIa] et al. and the references therein, as well
as the original Kapitzas work of 1941, translated in [KAh, p. 561].
(This kind of resistance is well pronounced at very low temperatures,
increasing dramatically at temperatures T < 20 K.) Note also that
the more general ellipsoidal inclusion problem for the resistance case,
i.e. under the boundary condition (4.45), is much more complicated
than that for the ideal contact. The analytical solution is obtained
in the paper of Benveniste and Miloh [BEf]; our elementary formula,
Eqs. (4.41) and (4.42), for spherical shape is given there as a particular case (cf. their Eq. (82)); the skin constant is denoted in this paper
as (so that = 2 /(aR) in our notations).
Remark 4.4 It is worth noting that the exact counterpart of the
Kapitza resistance in elasticity context is the so-called Hooke law
for a thin inhomogeneity, spread along the 2-D surface [SOa]. More
precisely, for a crack-like inclusion, when the stress vector n is
continuous on , one has
n (x) = (x) [u(x)] ,

x ,

(4.46)

13
Often the condition (4.45) is associated with the interfacial thermal barrier
resistance; in composites this effect arises due to poor mechanical or chemical
adherence, mismatch of the coefficients of thermal expansions leading to interfacial
gaps (microcracks). For more details and references see the papers of Hasselman
et al. [HAo], [HAp].

1.4 Micromechanics of Heterogeneous Media

99

which means that the opening of the inclusion, [u(x)], is proportional to the normal stress vector n (x) on the surface . The secondrank tensor (x) in (4.46) is just the Hooke tensor for the inclusion,
see [SOa], [KAc] for more details and discussion.
1.4.4 Elastic Case
In this case the counterpart of Eq. (4.1) reads

L(x) = L1 + [L]hw (x) ,

[L] = L2 L1 ,

(4.47)

where hw (x) is again the characteristic function of the inhomogeneity W, L1 and L2 are the tensors of the elastic moduli of the matrix
and of the inhomogeneity, respectively. At the absence of body forces,
the displacement field u(x) in the solid obeys the equation
{L1 : u(x)} + {hw (x)[L] : u(x)} = 0 ,

(4.48)

see Eqs. (2.54) and (4.47).

Similarly to the scalar case (Section 1.4.2), the basic integral equation which governs the strain tensor field in the solid with a single
inhomogeneity, follows directly from Eq. (4.48). In a Cartesian system it has the form
in (x) =
in (x) +

Z
W

Gipjn (x x0 ) [L]pjkl kl (x0 ) dx0 ,

(4.49)

where

Gipjn (x) =

1
2

(Gip,jn (x) + Gnp,ji (x))

(4.50)

is the counterpart of the kernel (4.5), with Gij (x) denoting the component of the Green tensor for an unbounded matrix.
For an isotropic matrix the Green tensor is well-known
1
Gij (r) =
41

1
ij + 1 r,ij ,
r


1 =

1
,
4(1 1 )

(4.51)

r = |x|, where 1 and 1 are the shear modulus and the Poisson ratio
of the matrix.
Let = const; from Eq. (4.49) it is easily seen that the strain
within the inhomogeneity W will be also constant, if the tensor field

Pipjn (x) =

Z
W

Gipjn (x x0 ) dx0

i
1 h1
=
(,jn (x),ip + ,ji (x),np ) + 1 ,ipjn
2
1

(4.52)

100

Heterogeneous Media: Modelling and Simulation

is constant within the inhomogeneity W; here = w (x) is the

Newtonian potential (4.8) and
1
(x) = w (x) =
4

|x x0 | dx0

(4.53)

is the biharmonic potential for the region W; recall that the latter
solves the equations
w (x) = 2w (x) ,

w (x) + 2hw (x) = 0 .

(4.54)

In turn, P will be constant within W, if w (x) is a quadratic function

and w (x) is a polynomial of fourth degree there. The latter is again
true for an ellipsoidal inhomogeneity. We can thus conclude, similarly
to the scalar case, that if the strain tensor at infinity is constant, it
will be also constant within any such inhomogeneity, see e.g. [ESa],
[KUb], [CHd], etc. Moreover, as it follows from Eq. (4.49),


1

Aw = I + P : [L]

(4.55)

where I is the unit fourth rank tensor with the Cartesian components Iijkl = 12 (ik jl + il jk ).
The tensor Aw in Eq. (4.55), that interconnects the constant strain
at infinity with the constant strain in the ellipsoid, is often called in
the literature Wus. (Wu was the first who specified this tensor for
a spheroid, averaged over all possible orientations, and then treated
the limiting cases of a disk and fiber in order to clarify the influence
of the shape upon the effective elastic properties [WUa].)
Remark 4.5 Note, without going into details, that once the strain
tensor within the inclusion is found, i.e. the tensor Aw is specified,
the displacement and strain fields outside the inclusion can be immediately obtained by means of the basic integral equation (4.49),
cf. Remark 4.3, where the similar question was discussed in the scalar
case.
Consider first, similarly to the scalar case, the most important
particular shapes of the inclusion W.
Let the inhomogeneity W be a sphere. Both tensors w (x)
and w (x) are then constant and isotropic, if x W;
w (x) is given in Eq. (4.11) and
,ijkl = Hijkl ,

Hijkl = ij kl + ik jl + il jk ,

x W , (4.56)

1.4 Micromechanics of Heterogeneous Media

101

since ,ijkl should be completely symmetric with respect to its indices.

Making a full contraction in (4.56) and utilizing Eq. (4.54)2 , we find
= 2/15, so that
2
H,
w (x) = 15

xW.

(4.57)

Eqs. (4.11), (4.57) and (4.52) now specify the needed tensor Aw :

Aw =

k1
1
I0 +
I00
k1 + 1 [k]
1 + 1 []

I0ijkl = 31 ij kl ,

I00ijkl =

1
2

(sphere) ,


ij kl + il jk 23 ij kl ,

3k1
1 + 1
1 =
=
,
3k1 + 41
3(1 1 )
1 =

(4.58)

6 k1 + 21
2 4 51
=
,
5 3k1 + 41
15 1 1

with k1 = 1 + 23 1 denoting the bulk modulus of the matrix.

Let the inhomogeneity W represent a circular fiber along the
axis x3 . The second gradient w in this case has been found in
Eq. (4.14). We introduce the 2-D counterparts D = kDijkl k of the
fully symmetric tensor H from Eq. (4.56), i, j = 1, 2. Then ,ijkl =
Dijkl ; using Eqs. (4.54) and (4.53) gives = 1/4 and hence
w = 41 D .

(4.59)

Let us recall the short-hand notation of Hill [HId]. If a pair of

symmetric second-rank tensors , T with the Cartesian components
ksij k and ktij k, respectively, are interrelated in the same Cartesian
system as
1
2

s11 + s12 = a(t11 + t12 ) + b0 t33 ,

s33 = b00 (t11 + t12 ) + ct33 ,

s12 = 2dt12 ,

s13 = 2et13 ,

s11 s12 = 2d(t11 t22 ) ,

(4.60)

s23 = 2et23 ,

then the transversely isotropic fourth-rank tensor A (with the symmetry axis x3 ), such that = A : T , is symbolically written as

A = (2a, b0 , b00 , c, 2d, 2e) .

(4.61)

Heterogeneous Media: Modelling and Simulation

102

Using these notations and the already found gradients of the potentials w and w , see Eq. (4.59), we have eventually

Aw = (2A, B 0 , B 00 , C, 2D, 2E) (fiber),

1 + 21
,
2 + 2 + 1

2A =

B 00 = 0 ,

C = 1,

21
,
2 + 1

2E =

g=

2B 0 =
2D =
2
,
1

[]
,
2 + 2 + 1

4(1 1 )
,
1 + (3 41 )g

(4.62)

[] = 2 1 .

Let the inhomogeneity be a layer h x3 h, perpendicular

to the x3 -axis. The potential and will then depend on x3 solely
and utilizing once again simple symmetry arguments yields
ij = 3i 3j ,

ijkl = 23i 3j 3k 3l .

These expressions, together with Eqs. (4.55) and (4.52), give eventually
Aw = (2A, B 0 , B 00 , C, 2D, 2E) (layer) ,
B0 = 0 ,

2A = 1 ,
C0 =

1 + 21
,
2 + 22

B 00 =
2D = 1 ,

[]
,
2 + 22
2E =

(4.63)

1
.
2

Consider next the more general case of a spheroidal inhomogeneity W. The needed second gradient of the Newtonian potential is
given in Eq. (4.20). In turn, using Hills notations (4.60) and (4.61),
and appropriate results, listed in the books of Mura [MUa] or Muratov [MUb], we have, for x W,
w (x) = (4C 0 , C 00 , C 00 , 3D0 , 2C 0 , 2C 00 ) ,
C0 =

D0 =

1
4

1 2M + 38 D0 ,

C 00 = M 32 D0 ,
(4.64)

1 3M

if c < a ,

1 3M

(1 e2 )
,
2

if c > a ,

3e

3e

1.4 Micromechanics of Heterogeneous Media

103

where M , M , and the eccentricity e, let us recall, are defined in

Eq. (4.20).
In the isotropic case under study it remains to introduce Eqs. (4.20)
and (4.64) into Eq. (4.55) in order to find the needed fourth-rank tensor Aw that transforms the strain in infinity into the strain within
the spheroidal inhomogeneity. After some algebra, the eventual result
reads
Aw = (2A, B 0 , B 00 , C, 2D, 2E) (spheroid) ,

1
1
1+
M [] + 21 C 00 []
2
1


2A =

2B 0 =

o
1 n
M [] 1 C 00 [] ,
1

2B 00 = l

o
1 n
M [] 1 C 00 [] ,
1


1
1
C=
1+
2M [] + ([] + 1 C 00 )[]
2
1
1
2D =
,
1 + 2(M + 21 C 0 )[]
1
2E =
,
1 + (M + M + 41 C 00 )[]


2 = 1 +
+

(4.65)


 o
1 n
[k] + 31 C 00 + M + 13 []
1


3
00

C
+
2M
M
1
[k][] ,
12

1 =

1 + 1
1
=
,
1
1 21

i = i + 2i , i = 1, 2.
Remark 4.6 Once again, for the particular values of the parameters M and M , as listed in Eq. (4.22), the particular forms of the
tensor Aw are recovered from (4.65) for a sphere, fiber and layer,
cf. Eqs. (4.58), (4.62) and (4.63), respectively.
Let the spheroid be randomly (and uniformly) oriented in the
space. Then
hAw i = a0 I0 + a00 I00 ,
(4.66)
where hAw i denotes the value of the tensor Aw averaged with respect
to all possible orientations of the axes xi of the spheroid; the isotropic

Heterogeneous Media: Modelling and Simulation

104

tensors I0 and I00 are defined in Eq. (4.58). After some algebra, it
appears eventually
a0 =

1
3

4A + 2B 0 + 2B 00 + C

1
1
=
1+
M+
2
1


a00 =

1
15

1
3

00



+ 31 C []

(4.67)


2A 2B 0 2B 00 + 2C + 12D + 12E .

Unlike a0 , we were not able to find a simpler expression for the

coefficient a00 .
The quantities a0 and a00 have been first calculated by Wu [WUa],
and also by Kuster and Toksoz [KUc], some misprints are corrected
in [BEm]. A simple check shows that though different functions are
used (e.g. arcsinh ), the basic parameter , used in the papers [WUa],
[BEm] is simply 1 M , where M is the depolarization factor from
Eq. (4.20).
The Crack Limit. In Section 1.5, when dealing with the
effective properties of microcracked solids, we shall need information
about the asymptotic behavior of the tensor Aw , given in Eq. (4.65),
when the ratio = c/a 0. The reason is that the cracks will be
treated as degenerating spheroids, in the sense that while two of the
axes a = b are kept fixed, the third one, c, becomes negligibly small,
i.e. c/a 0. From the definitions of the respective parameters, see
Eqs. (4.20) and (4.65), it is easily checked that the following formulae
hold as 0, namely,
M =1
C 00 =

+ O() ,
2

+ O() ,
4

+ O() ,
4

M =
D0 =

2
3

= 0 + O() ,
A=
C=

1 1 21 1
+ O(1) ,
4 1 21
4 1 21 1
+ O(1) ,
1 21

+ O() ,
2
0 =

C 0 = 0.0625 + O() ,
=

2 1
+ O(1) ,
3

1 21
,
8 (1 21 )2

B 0 = 14 (1 1 ) + O() ,
B 00 =

(4.68)

4 1 (1 1 ) 1
+ O(1) ,
1 21

1.5 Micromechanics of Heterogeneous Media

D=

1
2

0.0125
+ 1
1 1


+ O() ,

a0 = 31 (4A + 2B 0 + 2B 00 + C) =
a00 =
=

1
15 (2A

E=

105

2 1 1 1
+ O(1) ,
2 1

4 1 12 1
+ O(1) ,
3 1 21

2B 0 2B 00 + 2C + 12D + 12E)

8 (1 21 )(5 1 ) 1
+ O(1) .
15
2 1

1.5 One-Particle Approximations

For simplicity, we shall consider the so-called particulate media, in
which one of the constituents, called matrix, is connected. Throughout the latter particles of a second constituent (filler), with a welldefined, in general, shapes are distributed. The properties of the
matrix and of the filler are denoted, similarly to the previous section,
by 1 and 2, respectively, so that, e.g. 1 and 2 are the appropriate
conductivities, etc. Such a material is called also a dispersion (to be
thus distinguished from the case of suspension, when the particles are
immersed into a liquid and hence their relative positions are not fixed
due to the flowa fact that tremendously complicates any theoretical
study).
Our aim here will be to demonstrate how various repetitions
(to cite again Maxwell [MAj, Section 314]) of the single-inclusion
solutions of Section [Link] can be systematically utilized in order
to derive the basic and most widely used approximate models for
predicting the effective properties of particulate media.
The basic idea is to invoke formulae like (2.41) that interconnect
the effective properties with the appropriate concentration factors, as
introduced in Sections [Link].8. The latter are in turn approximated by the appropriate solutions of the single-inclusion problem.
Two natural possibilities are then opened. They can be summarized as follows:
(i) Each inclusion is treated as single, immersed into a homogeneous medium. The latter, as a consequence of the presence

106

Heterogeneous Media: Modelling and Simulation

of the host of other inclusions, possesses the unknown effective
properties (the effective medium idea).

(ii) Each inclusion is again treated as single, but now immersed into
the matrix. The presence of the rest of inclusions is accounted
for by imaging that the inclusion is subjected to a certain external field that differs from the applied macroscopic one (the
effective field idea).
We shall consider below some of the simplest implementations of
the assumptions (i) and (ii), within the frame of the static behavior of
heterogeneous media. Their application and generalization for wave
propagation problems in such media are discussed and treated in
detail in Chapter 3 of this volume.
1.5.1 Dilute Filler FractionScalar Conductivity
Let us start with the simplest case, when there is no need either to
modify the mediums properties or to change the applied field.
Namely, let the dispersion contain identical particles, of a fixed
shape W, aligned throughout the matrix. We assume that the dispersion is dilute, in the sense that the volume fraction, 2 , of the
particles is small, 2  1, and hence their mutual influence can be
neglected. Each particle can be thus imagined as single, immersed
into an unbounded matrix, and the concentration factor A2 is simply the single-inhomogeneity tensor Aw (1 , 2 ), defined in Eq. (4.10).
Thus, as it follows from (2.41),
K = 1 I + 2 []Aw (1 , 2 ) + O (2 )

(5.1)

both constituents are assumed isotropic but, due to the alignment

and nonspherical, in general, shape of the particles, the effective tensor K may be anisotropic. As underlined in (5.1), the formula is
correct to the order O (2 ) only, providing however in this case an exact expression (since the particle interactions are totally neglected).
To get the successive terms of order 22 and higher, these interactions
should be somehow accounted for. But, naturally, all schemes considered below, that try to incorporate the interactions, agree in the
dilute limit. More precisely, they all reproduce the formula (5.1) in
this case.

1.5 Micromechanics of Heterogeneous Media

107

Remark 5.1 Note that the approximation (5.1) has been used in
many of the earlier works on conductivity of two-phase media also in
the non-dilute case, see [REb] for some references. Perhaps its origin
can be traced back to Laplace himself: According to Lorentz [LOb,
Section 123], considering the connection between the density, , of
a transparent solid and its refraction index, 2 , Laplace reached the
formula
2 1
= const .
(5.2)

But 2 is just the dielectric constant of the medium, as already

pointed out (see the footnote on page 10), and the density is proportional to the volume fraction of the dielectric spheres contained
in the medium. Then (5.2) can be interpreted, somewhat loosely of
course, as the dilute approximation of the type (5.1) for the effective
dielectric constant of a dispersion.
Remark 5.2 Dilute formulae like (5.1) can be derived for heterogeneous media in different physical context (see, e.g. Section 1.5.7
below where the absorption problem is addressed). The so-obtained
results will make sense, however, only if the solution of the appropriate single-inclusion problem exists and is finite, like the tensor Aw
in Eq. (5.1). There are, unfortunately, situations when this solution does not exist (say, for the same absorption problem in 2-D, see
again the same Section 1.5.7). Then the appropriate effective properties have no representation of the type (5.1); nonanalytic functions,
like ln 2 , may emerge for small 2  1 instead.
Consider in more detail the most important particular cases of the
formula (5.1).
Assume the inhomogeneities are spheres. From Eqs. (4.13) and
(5.1) it follows

= 1 + 32 + O(2 ) ,
1

[]
,
1 + 22

(5.3)

which, naturally enough, reproduces the dilute limit (1.20) of the

Maxwell prescription for .
Consider a fiber-reinforced material, in which inhomogeneities are parallel circular cylinders (along the axis x3 ). The effective

Heterogeneous Media: Modelling and Simulation

108

conductivity tensor is transversely isotropic

K = 11 (e1 e1 + e2 e2 ) + 33 e3 e3 ,

(5.4)

where
22

2[]
1+
2 + O(2 ) ,
1 + 2

11

33

= = 1 1 + 2 2 ,

= 1

(5.5)

see (5.1) and (4.15). The coefficients 11 = 22 represent the transverse conductivity of the fiber-reinforced medium, for dilute fiber
concentration 2  1, in the case when the applied temperature
gradient is perpendicular to fibers axes. The conductivity component 33 along the x3 -axis is exactly the arithmetic mean , which
explains why the higher order correction O(2 ) is missing in the second formula of (5.5). The reason is that the temperature field in
the medium is just (x) G3 x3 , if the applied gradient is along the
fibers, i.e. if G = G3 e3 a situation already discussed in Remark 2.1,
see Eq. (2.11).
Consider a layered material that consists of laminae, perpendicular to the axis x3 . The appropriate tensor Aw is given in (4.17),
so that Eq. (5.1) implies
11 = 22 = 1 + 2 [] + O(2 ) ,

33 = 1 + 2 []

1
+ O(2 ) .
2

However, for such an exceptionally simple internal constitution the

following exact formulae for the effective conductivities hold
11 = 22 = ,

i.e. 1/ = 1/ = 1 /1 + 2 /2 .

(5.6)

The reason is that the temperature gradient remains constant within

the volume V, if G is along the laminae. If G is along the axis x3 , i.e.
perpendicular to the laminae, the heat flux is constant within V; in
this case the effective compliance k = 1/ is the arithmetic mean
of the constituents compliances, see Remark 2.2 and Eq. (2.18).
Let the inhomogeneities be identical and isotropic spheroids,
aligned along the x3 -axis, with the semiaxes c and a, c 6= a = b.
In this case the tensor Aw is given in Eq. (4.21), the medium is

1.5 Micromechanics of Heterogeneous Media

109

transversely isotropic, with symmetry axis x3 , and its effective conductivity tensor K is
K = 11 (e1 e1 + e2 e2 ) + 33 e3 e3 ,
11
33

22

= 1

= 1

[]
1 + 2
1 + []M

[]
1 + 2
1 + []M

+ O(2 ) ,

(5.7)

+ O(2 ) .

For the particular values of the parameters M and M , as listed

in Eq. (4.22), the already given formulae (5.3), (5.5) and (5.6) for
the effective conductivities in the cases of spheres, fibers and layers,
respectively, show up.
Let the inhomogeneities be again identical, but randomly oriented and isotropic spheroids, with the semiaxes c and a. In this
case the tensor Aw = hAw i is given in Eq. (4.23), the medium is
macroscopically isotropic and its effective conductivity reads

[]
=1+
1
31

21
1
+
1 + []M 1 + []M

2 + O(2 ) .

(5.8)

The parameters M and M depend on the eccentricity of the spheroids, see Eq. (4.20).
Let the spheroids, again randomly oriented in the space, degenerate into penny-shaped cracks. In this case 2 = 0 and the term
in the brackets in the right-hand side of (5.8) has the asymptotic
m=

2
1
2 1
+
=
+ O() ,
1 M 1 M

= c/a 0 ,

(5.9)

as it follows from Eq. (4.68). Eq. (5.8) now yields

c

= 1 n 43 a2 cm = 1 (na3 ) 43 m
1
a
=1

43 m

+ O() 1
0

(5.10)

8
9 ,

in virtue of Eq. (5.9). Here n is number density of the spheroids

(i.e. their number per unit volume), so that 2 = nV , V = 43 a2 c is
the volume of the single spheroid, and
= na3

(5.11)

110

Heterogeneous Media: Modelling and Simulation

is the so-called crack density parameterthe counterpart of the volume fraction of the inclusions in the case of cracked solids (as introduced and very often employed for such solids, see the survey of
Kachanov [KAa]). Thus the effective conductivity of a solid, containing a dilute population of randomly oriented penny-shaped cracks, is

= 1 89 + O() ,
1

(5.12)

1.5.2 Dilute Filler FractionElastic Case

In the elastic case and for a dilute particle concentration we have, as
a consequence of Eq. (2.64),

L = L1 + 2 [L] : Aw (L1 , L2 ) + O(2 ) ,

(5.13)

where Aw is the fourth-rank tensor that comes out from the solution
of the elastic single inclusion problem, see Section 1.4.4.
Remark 5.3 Since the formula (5.13) is exact, the tensor
[L] : Aw (L1 , L2 )
possesses the needed external symmetry (2.65), whatever the shape
of the inclusions W.
Similarly to the previous subsection, consider the most interesting
particular cases of the formula (5.13).
For an isotropic dispersion of spheres the tensor Aw is given
in Eq. (4.58) and hence
[k]
k
=1+
2 + O(2 ) ,
k1
k1 + 1 [k]

[]
=1+
2 + O(2 ) ;
1
1 + 1 []

(5.14)

1 and 1 depend on the Poisson ratio, 1 , of the matrix, see the

same Eq. (4.58).
Remark 5.4 The formulae (5.14) have been first derived by Bruggeman [BRg]. (Though, as pointed out in [ESa], his result for the shear

1.5 Micromechanics of Heterogeneous Media

111

modulus is apparently incorrect, since the Poisson ratio of the matrix

does not appear.) The correct formula for was given by Dewey
[DEb]. The particular case of (5.14) for rigid particles was found
by Hashin [HAc]. Eshelby [ESa] obtained them as a consequence
of his well-known solution for the ellipsoidal inhomogeneity problem. Hashin [HAd] rederived independently the same formulae and,
on their base, proposed a method how the nondilute case could be
treated. In the wave propagation problem for a dilute dispersion,
the same formula (5.14) reappeared as a zero-frequency limit [CHa],
[MAb]. (In the Russian literature Eqs. (5.14) are often attributed to
Krivoglaz and Cherevko [KRc].) As it seems, Eqs. (5.14) are one of
the most often rederived formulae in the elasticity of two-phase media
for a period at least of forty years (19371957); the author suspects
that the foregoing list of references is far from exhaustive.
Remark 5.5 For a porous solid (k2 = 2 = 0), the formulae (5.14)
simplify
3 1 1
k
=1
2 + O(2 ) ,
k1
2 1 21
(5.15)

1 1
= 1 15
2 + O(2 )
1
7 51
in this case they are due to Mackenzie [MAa].
Remark 5.6 Let the matrix be incompressible, so that k1 = ,
1 = 0.5, and the spherical particles are rigid, k2 = 2 = . The
shear modulus of such a mixture is then

= 1 + 2.52 + O(2 ) ,
1
which is exactly the Einstein formula (1.21) for the viscosity of a
dilute dispersion. (Observe the full mathematical analogy between
the slow (Stokesian) flow of the dispersion, with fixed position of the
rigid particles, and the elastic deformation of a solid with the same
particles distribution, provided both fluid and solid are incompressible.) In the other limiting case of voids in an incompressible matrix,
the second of Eqs. (5.15) yields

= 1 53 2 + O(2 ) .
1

112

Heterogeneous Media: Modelling and Simulation

the counterpart of the Einstein formula, giving the viscosity of a

dilute bubbly liquid in the Stokesian approximation. This result is
attributed to Taylor [TAd].
Let the inhomogeneities be identical but randomly oriented
spheroids. Then the material is macroscopically isotropic and (5.13)
implies
k
[k] 0
=1+
a 2 + O(2 ) ,
k1
k1

[] 00
=1+
a 2 + O(2 ) ,
1
1

(5.16)

with the coefficients a0 and a00 , given in Eq. (4.67).

Let 2 = 0 and the spheroids, again randomly oriented in the
space, degenerate into penny-shaped cracks. Using the asymptotic
(4.68) of the coefficients a0 and a00 in the crack limit, together with
Eq. (5.16), one gets
k
16 1 12
=1
+ O() ,
k1
9 1 21

32 (1 1 )(5 1 )
=1
+ O() ,
1
45
2 1

(5.17)

having recalled the definition (5.11) of the crack density parameter .

The formulae (5.17) have been independently proposed by Bristow
[BRc] and Walsh [WAa].
Remark 5.7 Treating the crack as a vacuous oblate spheroid of eccentricity approaching zero is very convenient. To the best of the
authors knowledge it was first employed by Willis [WIb]. The reason
is that it allows one to derive straightforwardly the needed formula
for a microcracked solid, passing (with a due care) to the appropriate limits in the general formulae, concerning ellipsoidal inclusions.
The usefulness of this idea was already demonstrated in deriving
Eqs. (5.12) and (5.17). In the elastic context an alternative derivations employ specific notions and results from fracture mechanics, like
stress-concentration factors, energy release rate, etc. Such a derivation, employed in particular in the original Bristows paper [BRc], or
in [BEc], is often difficult for a fracture outsider to follow.

1.5 Micromechanics of Heterogeneous Media

113

1.5.3 Self-Consistent Scheme

When the volume fraction becomes higher (say, 2 0.05) it is not
already realistic to imagine the inclusions isolated. Somehow their
mutual interactions should be accounted for. The simplest way to
do this, at least qualitatively, is the application of the self-consistent
idea that appeared, in a certain form, already in the Clausius derivation [CLa] of the Maxwell formula, mentioned in Section 1.1. Next
it showed up as a very useful and natural tool in solid-state physics
and, more specifically, in the wave mechanics of atoms, as originated
by Hartree [HAb]. In the context of effective properties of heterogeneous media the self-consistent approach seems to be first employed
by Bruggeman [BRf], [BRg], see also [FRd], [LAb]. In the quest of
the elastic constants of polycrystals, the method is attributed to Hershey [HEb] and Kroner [KRd]. (In Section 1.6 below we shall present
in more details the reasoning and some results of these two authors.)
[Link] The Basic Idea
The basic idea of the self-consistent scheme (SCS) is to treat each
particle again as isolated, but immersed into a matrix possessing the
unknown effective properties of the composite. This means that we
should replace Eq. (5.1) and its elastic counterpart (5.13) by the
equations
K = 1 I + 2 []Aw (1 , K )

L = L1 + 2 [L] : Aw (L1 , L )

(scalar case) ,
(5.18)
(elastic case) .

In each particular case, Eqs. (5.18) will yield appropriate algebraic

equations for the unknown effective properties of the composite.
To demonstrate this, consider the simplest example when the inclusions are spheres. In the scalar context, Eq. (4.13) is to be introduced into the first equation of (5.18):

[]
3
=1+
2
(5.19)
1
1 2 + 2
which represents, obviously, a quadratic equation with respect to the
effective conductivity.
Eq. (5.19) was proposed by Bruggeman [BRf], see also [LAb], in
the equivalent form
1
2
1
+

= 0.
(5.20)
2
1 + 2
2 + 2

114

Heterogeneous Media: Modelling and Simulation

Note the symmetry of the latter with respect to the constituents

characteristics, that is under the simultaneous interchange 1 2
and 1 2 . That is why Eq. (5.19) is often called symmetric
effective medium approximation (EMA), or First Bruggemans symmetric EMAa term, introduced by Landauer [LAc]. Thus each
phase should be treated on the same footing, as pointed out by Davies
[DAb]; in the above reasoning however one of the phase was imagined
continuous, and the second in the form of distributed inclusions; this
already implies an inherent logical inconsistency of the SCS. More
fundamental, however, is the inconsistency, very clearly exposed by
Hashin [HAi], namely, the violation of the basic MMM principle, formulated in Section 1.2.1. The point is that a single sphere in the
dispersion is the basic micro element and as such it can only see
the other spheres (micro-elements) around, but not the mini-element
(the homogenized medium) as a whole, whose properties result from
averaging over the multitude of all the micro-elements within it. Or,
as Hashin put it: SCS assumes that a tree sees the forestbut a tree
sees only other trees.
The foregoing does not imply however that the SCS is useless and
should be rejected. It simply suggests that if the scheme is realistic,
the appropriate two-phase geometries cannot resemble a dispersion
of equi-sized spheres. That this is indeed the case follows from an
important result of Milton [MIe]. It states that the SCS (more precisely, Eq. (5.19) or Eq. (5.20)) is a realizable approximation in the
sense that there exists a special, ingeniously devised, internal geometry of a two-phase medium whose effective elastic moduli coincide
exactly with the SCS predictions of Eq. (5.19). The reader is referred
to the original Miltons work since the details and proofs cannot be
elementary summarized. We can only add that, like Hashins assemblage (Section 1.3.4), the SCS-geometry of Milton is of fractal-type
(self-similar in all length-scales, that is), comprising inclusions of all
possible sizes, from finite down to infinitesimally small. Due to the
foregoing remarks, it is a little wonder that such a geometry has
nothing in common with equi-sized particles like spheres, distributed
throughout a matrix.
Remark 5.8 The choice of the basic element, to be inserted into the
effective medium, is not unique in the SC scheme. In the foregoing,
the simplest possibility was chosen, when this is just the filler parti-

1.5 Micromechanics of Heterogeneous Media

115

cle. This is however not the only choice; more complicated elements
can be devised as well. For example one can take a composite
sphere, with a kernel of phase 1, coated by a layer of phase 2,
as first suggested by Kerner [KEa]. A detailed and critical analysis
of such a SCS modification can be found in [HAg]; for generalizations in the elastic case, see [CHe] or [CHd, Chapter 2.3]. This kind
of modification of the self-consistent method are usually referred to
as Generalized Self-Consistent-Schemes (GSCS), but their discussion
goes beyond the frame of the present survey. It should be only noted
that without detailed morphological information, no theoretical arguments can be given for any specific choice of a basic element in a
binary medium. (The composite sphere was used in [CHe] or [CHd]
because the appropriate model provided better agreement with certain precise experimental data; though no guarantee exists whether
it will work equally well for other sets of experiments.) For polycrystals, however, the basic element is naturally predetermined from the
very internal structurethis is just the single grain (see Section 1.6).
The latter perhaps explains why the SCS provides as a rule more
accurate predictions for the effective properties of polycrystals than
of particulate media.
Note finally that the self-consistent method can be employed in an
obvious manner to multiphase media as well. However, its predictions
will possess the needed symmetry (2.65) only when the dispersed
inclusions are aligned and of similar shape, see the paper [BEe] for
details, examples and comments on this important fact.
[Link] An Alternative Derivation
Eq. (5.20) can be also recast as
(x)
= 0,
(x) + 2

(5.21)

where the overbar indicates, let us recall, averaging over the RVE,
and (x) is the function that takes the values 1 or 2 depending on
whether x lies in the constituent 1 or 2, respectively. Eq. (5.21) has
a simple and clear interpretation. Namely, imagine that the RVE, V,
is filled with spherical particles, whose sizes vary from finite down to
infinitesimal. The spheres have conductivity 1 or 2 . Imagine each

116

Heterogeneous Media: Modelling and Simulation

sphere lying in a homogeneous medium with the effective conductivities ; at infinity the prescribed temperature gradient G is imposed.
The deviation from this gradient within the sphere, due to the fact
that its properties differ from those of the surrounding material, is
just
(x)
G,
(5.22)
(x) + 2
cf. Eqs. (1.14) and (1.15). (In the electrostatic context, the latter
value is proportional to the so-called polarization of the sphere.)
Since the volume V is considered as macroscopically homogeneous,
we should have, accordingly, no macroscopical polarization. The
mean value of the local polarizations (5.22) should then vanish,
as claimed in Eqs. (5.20) and (5.21).
The above reasoning can be interpreted as well in the following
manner, that represents a modification of the original Maxwell arguments of Section [Link]. Cut again a big sphere of the radius A (the
Maxwell sphere) from the heterogeneous medium and immerse it
into unbounded homogeneous host material whose conductivity +
will be appropriately chosen (or, so to say after Berryman [BEm],
tuned). Sit down at the point M far away, i.e. at the distance R
from the center of the sphere, R  A. Let the temperature field Gx
be applied at infinity. Each piece of the phases 1 or 2, located in
the vicinity of x, will create at the point M the disturbance
a3 (x) +
G,
R3 (x) + 2+

(5.23)

since we are far from the piece under consideration, it can be taken
as spherical with a radius a. Averaging (5.23) over the big heterogeneous sphere will give a quantity, proportional to
(x) +
.
(x) + 2+

(5.24)

Comparing (5.24) with (5.21) shows that if we tune the host medium
so that + = , the presence of the heterogeneous sphere will not
be felt at the point M ; in other words this sphere will become
transparent for us.
The above arguments represent a static (and scalar) implementation of Berrymans idea [BEm]: this author dealt with the longwave length limit for a harmonic wave, scattered by the Maxwell

1.5 Micromechanics of Heterogeneous Media

117

sphere, and tuned the properties of the host medium in such a

way that the total scattering is zeroin other words, the sphere becomes indeed invisible for us. Similar arguments in the language of
the so-called T -matrix were employed in [GUa], [KOa].
Remark 5.9 Sometimes the SC scheme is called the Coherent Potential Approximation due to the similarity of its implementation
(cf. Eq. (5.24)) to a certain approach, used in the theory of electronic properties of imperfect crystals and alloys, see [PAb] or the
excellent exposition in Zimans book [ZIa].
[Link] The Percolation Phenomenon
Going back to to the simplest SCS equation (5.19), note that it possesses the simple solutions

1
, if 2 /1 = ,
= 1 32
1
1 3 , if / = 0 ,
2
1
2 2

(5.25)

in the two limiting cases of ideally conducting particles (2 /1 = )

or non-conducting particles (2 /1 = 0), respectively.
Hence, a dispersion of ideally conducting spheres becomes ideally
conducting, if the volume fraction 2 1/3, whatever the conductivity of the matrix, the spatial distribution of the spheres and their
sizes. This is just an exhibition of the so-called percolation phenomenon and the value p2 (1/3 for the adopted SCS) is the percolation threshold. The essence of the phenomenon is the following:
when an ideally conducting constituent 2 is added to a poorly conducting matrix, it is most possible that it will be in the form of
separated particles, if the volume concentration 2 is small enough.
When increasing 2 , the particles start touching each other, forming
more complicated aggregates. At certain critical (threshold) value
p2 , the constituent 2 forms a continuous skeleton penetrating the
whole body of the mixture and as a result it becomes ideally conducting as well. A similar percolation phenomenon appears obviously in
the case when a nonconducting phase (2 = 0) is added to a conducting matrix. In this case the threshold concentration p2 = 2/3, again
independently of microstructural details, see Eq. (5.25).

118

Heterogeneous Media: Modelling and Simulation

It is noted that the percolation phenomena are very important

in solid state physics, in general, and in many of its applications.
The appropriate discussion goes however far beyond the scope of the
present work (see the survey [LAc] or [KIa] for basic ideas, results
and classical references).
Remark 5.10 In the literature the presence of a percolation threshold is sometimes considered as a virtue of an approximate model,
while its lack (in the differential scheme below, for example) is viewed
as a drawback (or vice versa, depending on authors taste). It should
be noted that the thresholds presence (or lack) has nothing in common with the validity of the model. It only informs us that for the
geometry that realizes the model (like Milton SCS-geometry), there
is (or there is not) a threshold. Hence this is a conclusion concerning above all a specific two-phase geometry, but not a model or its
applicability.

[Link] Some Particular Cases

Consider now, for illustration, several particular cases of the SCS
equations (5.18).
For randomly oriented and isotropic spheroids, in the scalar
conductivity case,
[]

=1+
1
31

+
+ (2 )M + (2 )M

2 , (5.26)

cf. Eqs. (4.21) and (5.18). When the spheroids are spheres, M =
M = 31 , and (5.26) reproduces (5.19).
Let 2 = 0, and the spheroids degenerate into penny-shaped
cracks. The effective conductivity is then eliminated from the
right-hand side of Eq. (5.26) and Eq. (5.12) is recovered, but without the correction term O(), cf. Eq. (5.9). Hence the self-consistent
scheme predicts the linear dependence

= 1 89
1
of upon the crack density parameter.

(5.27)

1.5 Micromechanics of Heterogeneous Media

119

The results in the elastic case are obtained in a fully similar way,
and we shall only list some of the basic ones.
For an isotropic elastic dispersion of spheres, the self-consistent
scheme yields
k = k1 + k
= 1 +
=

k2 k1
2 ,
+ (k2 k )

2 1
2 ,
+ (2 )

3k
,
3k + 4

(5.28)

6 k + 2
,
5 3k + 2

see Eqs. (5.18) and (4.58). This is an algebraic system with respect
to the unknown effective moduli k and which allows a convenient
numerical solution (by a straightforward iterative procedure).
The system (5.28) is usually attributed to Budiansky [BUa] and
Hill [HIe], see also [WUa]. (These authors were apparently not aware
of the earlier Skorohods paper [SKa], where the same classical selfconsistent reasoning was employed to derive (5.28).)
Very simple closed form solutions (5.28) can be found, if 1 = 2 =
0.2. Then, as pointed out by Budiansky [BUa], the Poisson ratio of
the medium is also 0.2, = 0.5.
Other simple solutions exist if the Poisson ratio of the matrix is
again 1 = 0.2, but the inclusions are rigid or represent voids. In
particular, for spherical voids in a matrix of Poisson ratio 1 = 0.2,
the solution of (5.28) is extremely simple

=
= 1 22
1
1

(k2 = 2 = ) ,

1
=
=
1
1
1 22

(k2 = 2 = 0) ,

(5.29)

and hence 0.2 in both cases.

The foregoing facts suggest that the value 0.2 of the Poisson ratio
is something special in the one-particle models of particulate media.14 A certain explanation can be sought perhaps in the fact that
14
It was Hashin [HAc] who first pointed out the curious fact that 0.2 in
a dilute dispersion of rigid spheres in a matrix for which 1 = 0.2.

120

Heterogeneous Media: Modelling and Simulation

at 1 = 0.2 the concentration factor for a single spherical inclusion

the tensor Aw , as given in Eq. (4.58), is proportional to the unit
fourth-rank tensor I.
Eq. (5.29) means that for a dispersion of spheres, in a matrix with
1 = 0.2, the percolation threshold, as predicted by SCS, is p2 = 0.5.
At sphere fractions, higher than 0.5, the solid either become rigid or
cannot sustain load. Again, this fact is neither a merit or a fault of
the SCS, but just a property of the special Miltons SCS-geometry
that lies behind the scheme.
Another simple solution of (5.28) emerges for the Hill medium,
when the shear moduli of the phases coincide, 1 = 2 = . Then,
from the first equation, one has
k k1
3k + 4
= 2
.
k2 k1
3k2 + 4

(5.30)

which reproduces the Hill formula (3.66). Hence the SCS gives the
exact value of the effective bulk modulus in the case of spherical
particles. (Recall, however, that the Hill results (5.30) holds for an
arbitrary two-phase geometry, which is macroscopically isotropic.)
For an isotropic elastic solid, containing randomly oriented
penny-shaped cracks, we get formally the same system (5.17), with
the only difference that the Poisson ratio of the matrix, 1 , is replaced
by this ratio of the composite medium, :
k
16 1 2
=1
,
k1
9 1 2

32 (1 )(5 )
=1
.
1
45
2

(5.31)

This is exactly the self-consistent system of Budiansky and OConnell

[Bub] for the effective moduli of a microcracked solid (cf. their Eqs.
(36) and (440 )).
1.5.4 The Differential Scheme
The differential scheme (DS) is a kind of an infinitesimal implementation of the self-consistent idea. It was initiated by Bruggeman
[BRf], [BRg] and used in fifties for media containing rigid particles
(by Brinkman [BRb] for the viscosity of a fluid suspension, and by

1.5 Micromechanics of Heterogeneous Media

121

Roscoe [ROa] in the elastic case). More recently it was elaborated by

McLaughlin [MCb], Norris [NOa], Zimmerman [ZIc] et al. A detailed
discussion can be found in the survey [CLb].
[Link] The Basic Idea
Imagine the small concentration of particles is added to an initially
homogeneous medium. Replace mentally the so obtained mixture by
a homogeneous medium of the effective conductivity K = K ().
In the latter a new concentration of particles is inserted; the
new medium again is replaced by a homogeneous one of conductivity K (2), and so forth. As a result, the effective conductivity
will be a certain function of the total amount of inserted particles,
K = K (). Since at each step a small concentration of particles is
added,  1, the dilute approximation (5.1) applies
K ( + ) = K () + (2 I K ()) Aw (K (), 2 ) .
In the limit 0, this yields the differential equation
dK
= (2 I K ()) Aw (K (), 2 )
d

(5.32)

for the unknown function K = K (), which should be solved together with the natural initial condition
K (0) = 1 I .
The reasoning in the elastic case is fully similar, and will be omitted.
It is very important to point out, after McLaughlin [MCb], that
the total concentration of particles, introduced into the matrix,
does not coincide with the volume concentration 2 of the second
phase. The reason is that when adding a new portion of particles,
some of them will fall within the region, occupied by the previously
deposited ones. Only the part 1 , that falls within the matrix,
will increase 2 . Thus, in the limit 0 again,
d2 = (1 2 ) d

or

= ln(1 2 ) ,

(5.33)

because, obviously, = 0 at 2 = 0. From (5.33) it follows that

= 2 + 12 22 ,

122

Heterogeneous Media: Modelling and Simulation

which means that 2 for 2  1. Hence only for small values of

2 one can choose = 2 .
Remark 5.11 In some of the earlier papers [HEa], [ZIb], the differential scheme was applied without distinguishing between the quantity and the fraction 2 . One of the motivation was Mackenzies
treatment of a porous solid [MAa], in which a new single pore is added
at the center of a sphere, cut from the homogenized medium, and subjected to hydrostatic pressure. As just pointed out, taking = 2
is safe only in the dilute case; for higher values of 2 the deviations
from the differential schemes predictions is already considerable.
Obviously, the phases do not enter symmetrically in the Eqs. (5.32)
and (5.33), unlike the self-consistent scheme. On the other hand again
the effective medium surrounding of the inclusion is imagined at each
step of the DS procedure. That is why the DS is often called, after
Landauer [LAc], the Second Bruggemans unsymmetrical Effective
Medium Approximation. Also the term iterated dilute limit approximation is sometimes used, which reflects well the reasoning that has
led us to the formulae (5.32) and (5.33).
The differential scheme, similarly to the self-consistent one, is a realizable approximation in the sense that there exist special two-phase
media whose effective behavior coincides exactly with the predictions
of the scheme. This important fact was first demonstrated by Norris
et al. [NOa], [NOb], [NOc], see also [AVa].
[Link] Some Particular Cases
Consider now in more detail the system (5.32) and (5.33) in several
particular cases.
For spherical inclusions in the scalar conductivity context the
tensor K = I is isotropic and Aw is given in (4.13). Eq. (5.32)
then becomes
df
3( f (2 ))f (2 )
=
, /1 = f (2 ) ,
(5.34)
d2
(1 2 )( + 2f (2 ))
= 2 /1 . Similarly to the self-consistent scheme, Eq. (5.34) allows
simple solutions

1
,
(1 2 )3
(1 2

)3/2 ,

if 2 /1 = ,
if 2 /1 = 0 ,

(5.35)

1.5 Micromechanics of Heterogeneous Media

123

in the two limiting cases of ideally conducting particles (2 /1 = )

or non-conducting particles (2 /1 = 0), respectively. Hence, there
is no percolation in the differential schemethe values of remain
finite and positive for all 2 (0, 1).
For spherical inclusions in the elastic case the DS supplies the
system
1 dk
(3k + 4 )(k2 k )
=
,

k d2
(1 2 )k (3k2 + 4 )
1 d
5(3k + 4 )(2 )
=
,
d2
(1 2 )(6 (k + 2 ) + (9k + 8 ) )

(5.36)

coupled with the obvious initial conditions k(0) = k1 , (0) = 1 .

The system (5.36) was derived by McLaughlin [MCb] who showed, in
particular, that its solution meets the HS bounds.15
A detailed analysis of the system (5.36) was performed by Zimmerman [ZIc]. He observed that (5.36) is autonomous, with the parameter , see Eq. (5.33), playing the role of time. This suggests
to divide the equations of (5.36) and the result will be an explicit,
though complicated, expression of 2 through the ratio k / (or
which is the same, through the effective Poisson ratio ). Simpler
solutions exist in the limiting cases of rigid and vacuous inclusions.
The equation for that follows from (5.36) has in these cases the
form
d
= F (1 )(1 51 ) .
d
Though the function F (1 ) is different for rigid and vacuous inclusions, the presence of the multiplier 1 51 in both cases assures that
0.2, if 1 = 0.2something that is identical with the predictions
of the SCS. Moreover, the value 0.2 of the Poisson ratio is an attractor or a fixed point, in the sense that the function = (2 ) is
monotonically increasing, if (0) = 1 < 0.2, and monotonically decreasing, if (0) = 1 > 0.2. In both cases (2 ) 0.2 as 2 1.
The same attractors type behavior has been observed in the SCS by
the author [MAd].
15
Once an approximate scheme is realizable, its predictions always satisfy the
HS bounds, due to the optimality of the latter; see Section 1.3.3.

124

Heterogeneous Media: Modelling and Simulation

Remark 5.12 It is curious to understand whether the peculiarity of

the value 0.2 for the Poisson ratio is something that indeed exists in
real heterogeneous solids, or it is a mathematical artifact produced
by the models. This intriguing question explains the considerable
interests in the last 15 years in studying the so-called Poisson ratio
flow diagrams i.e. the effective Poisson coefficient as a function of
the volume fraction 2 of the inclusions. It turned out that if the
inclusions are voids, the fixed points on such diagrams (like the above
value 0.2) do exist, and this was demonstrated both analytically (see,
e.g. [CHb] in 2-D) and numerically by a number of authors. Details,
discussion and many references can be found in the recent paper
[CHf].
For rigid inclusions it seems as if there are no fixed points for
the Poisson ratio, as numerical experiments suggest, see again [CHf].
The apparent discrepancy between voids and rigid inclusions may be
however a by-product of the numerical procedures, employed by the
authors. After Zimmerman [ZId], one can point out that calculations
for a medium with vacuous pores can be performed simply putting
the elastic moduli 2 = 2 = 0 in the appropriate equations and
numerical schemes. But the analogous calculations for rigid inclusions
are performed taking the moduli 2 , 2 very large, but still finite.
One should be very careful, dealing with limit of moduli, tending
both to zero and to infinity. The problem arises because, when taking
either limit 0 or , any references of the Poisson ratio may be lost. In
fact, one can construct paths 2 , 2 0 or in (2 , 2 )-space so that
the Poisson ratio can have any value. For more details, concerning
this subtle situation, the reader is referred to the original paper [Zid].
Another similarity with the SCS is the existence of the simple solutions of (5.36), when 1 = 0.2, for the two limiting types of inclusions,
namely
k

=
= (1 2 )2 (voids) ,
k1
1
(5.37)
k

1
=
=
(rigid inclusions) ,
k1
1
(1 2 )2
cf. Eq. (5.35).
Remark 5.13 Relations of the type (1 2 )m for the decrease of
a certain effective property of a porous solid as a function of the

1.5 Micromechanics of Heterogeneous Media

125

void ratio 2 are very popular and useful in some applications like
micromechanics of microporous rocks. (In this context m is known
as the cementation exponent and its value supplies additional microstructural information, see [SEa] for details and references.) They
are often called Archies after the work [ARa]. The reappearance
of the phenomenological Archies law within the DS, cf. (5.35) and
(5.37), cannot be interpreted, however, as a kind of justification of
the scheme or as an indication of its superiority as compared to other
approximate theories. It means solely that the specific two-phase
Norris microgeometry that lies behind the DS, exhibits Archies
type behavior of its effective conductivity and elasticity (for 1 = 0.2
in the latter case).
Another particular case, in which the system (5.36) possesses a
very simple solution, is when both the matrix and spheres are rigid.
Then k1 = k2 = k = and
1

=
1
(1 2 )5/2

(5.38)

which reproduces the Einstein formula (1.21) in the dilute limit (recalling once again the full mathematical analogy between an incompressible elastic solid and a fluid in Stokesian approximation). The
formula (5.38) was proposed by Brinkman [BRb] for the viscosity of
a fluid suspension, within his attempt of extending the Einstein formula (1.21) in nondilute case. To this end the author has utilized
literally the same DS reasoning of Section [Link].
We shall skip further study of the DS model and equations, referring to the already cited papers of McLaughlin, Norris et al. and
Zimmerman. Also the application of the DS to the microcracked
solids will not be discussed (in this case some care is needed, due
to the degenerated nature of the inclusions). Details and further
references can be found in the Hashin paper [HAj], see also [HEa].
1.5.5 The Effective Field
Here again each particle is treated as single. However, we imagine
it immersed now into a homogeneous medium possessing the properties of the matrix. The influence of the rest of the particles is
accounted for through the assumption that the particle lies within

126

Heterogeneous Media: Modelling and Simulation

a certain external fieldthe effective fieldthat differs from the applied macroscopic one. Similarly to the self-consistent scheme, this
is an old approach that goes back to the Mossotti cavity idea and to
the Lorenz-Lorentz local field, as we shall point out below.
[Link] The Basic Idea
We shall first implement the effective field idea again in the scalar
case. Let G be the local (effective, or Lorenz-Lorentz in the dielectric context) field that acts on each particle; in this simplest version
it is assumed that G is a constant.
Since each particle is lying as single in the field, it is natural also
to assume that
hi1 = G .
(5.39)
As a matter of fact, this is the basic assumption of the simplest
version of the effective field methodit states that each particle undergoes an external field that coincides with the mean field in the
matrix phase.
In turn, G should depend linearly on the applied macroscopic
gradient G =
G = B G ,
(5.40)
with a certain second-rank tensor B. The latter is specified by means
of (5.39) and the first equation of (2.35):
G = = 1 hi1 + 2 hi2
= 1 G + 2 Aw (1 , K 2 ) G

(5.41)

= (1 I + 2 Aw (1 , K 2 )) G ;
we have used that hi2 = Aw (1 , K 2 ) G since each particle is
treated as single in the far-field G , according to the basic assumption. Hence
B = (1 I + 2 Aw (1 , K 2 ))1 ,
(5.42)
which gives explicitly the tensor B that interconnects the applied
macroscopic field and the resulting local one.
The concentration tensor can be easily found now by means of
(5.41) and (5.42):
hi2 = Aw (1 , K 2 ) G
= Aw (1 , K 2 ) B G = A2 G ,

1.5 Micromechanics of Heterogeneous Media

127

so that, in virtue of (5.42),

K = 1 I + 2 Aw (1 , K 2 ) (1 I + 2 Aw (1 , K 2 ))1

(5.43)

which is the effective field approximation for the bulk conductivity of

the medium.
In the elastic case, fully similar arguments produce the formula

L = L1 + 2 Aw (L1 , L2 ) : (1 I + 2 Aw (L1 , L2 ))1 ,

(5.44)

having replaced the second-rank tensors in (5.43) by the appropriate

fourth-rank tensors.
[Link] Alternative Derivations and Interpretations
Note first that the formula (5.41) and its consequences (5.43) and
(5.44) for the effective properties has an appealing interpretation in
terms of the Mossotti cavity idea (see Section [Link]). Indeed, guided
by the heuristic derivation of the Clausius-Mossotti formula, as used
by Feynman at al. [FEc, Chapter 11], let us cut the inclusions, treated
as isolated, from the medium and replace them by cavities. Then
G = Gcav + Gincl ,

(5.45)

where Gcav is the field in the matrix with cavities and, correspondingly, Gincl is the field that appears in the inclusions. The cavity
field Gcav is Gcav = 1 G , where G is the above mentioned mean
(or effective) field that results from the interactions of the inclusions.
In turn, Gincl = 2 Aw G , since each one is treated as isolated in
the far-field G . Together with (5.45), these formulae reproduce the
basic relation (5.41) of the effective field method.
An approach which turns out to be closely connected, if not equivalent, to the effective field approximation is associated with the names
of Mori and Tanaka [MOb] in the literature on composites, see also
[BEc], [WEa], [WEb] et al. The original derivation of these authors
was expressed in terms of eigenstrains, equivalent inclusions, transformation energy consideration, etc., and at first glance looked totally
different in spirit from the rest of the models, used in mechanics of
heterogeneous media. (Details and further applications can be found
in the book of Mura [MUa].) The clarification of the nature of this
approximation is to be attributed to Benveniste who, in the paper

128

Heterogeneous Media: Modelling and Simulation

[BEd], pointed out that the basic idea of Mori-Tanakas method is

just the statement of Eq. (5.39), namely, that each particle is again
treated as single, but lying within the far-field hi1 that equals the
mean field in the matrix. This assumption was adopted also by the
author in the note [MAe], without any connection to Mori-Tanakas
work, as a starting point of the derivation of the approximations
(5.43) and (5.44) for the effective moduli, as given in Section [Link].
As we shall see below, the effective field method works well and
produces reasonable results in the cases when the inclusions are either
spherical or, if nonspherical, are aligned. For multi-phase systems
with different alignment and/or shape of the particles, the effective
field predictions fail to satisfy the necessary symmetry conditions
(2.65) and thus they are not acceptable, see [BEe], [FEb]. The same
shortcoming of the SCS was already pointed out in Remark 5.8.
This failure, in particular, suggests that the effective field method
should be modified for media with more complicated internal structure. A possible way was proposed by Kanaun and Levin in a number
of publications and elaborated in their book [KAf], see also [KAg] for
a detailed English exposition. Using the integral equation for the
displacement field in a particulate medium, they have derived an
approximate self-consistent equation for the effective field. The
solution of the latter, though simple, involves however two-point
correlation functions for the inclusions distribution. For spherical
particles, isotropically distributed, the Kanaun-Levins scheme yields
the same results (5.46) and (5.47), given below as particular cases of
the foregoing effective field (Mori-Tanaka) method. Details and some
of the latest results, concerning the application of this approach to
scalar waves phenomena in particulate media, can be found in Chapter 3 of this volume.
[Link] Some Particular Cases
Consider again some of the most common and interesting particular
cases, to illustrate the predictions of the effective field method.
Let the medium be a macroscopically isotropic random dispersion of spheres. The tensors Aw and Aw are given in Eqs. (4.13) and
(4.58) respectively, so that the effective field prescriptions (5.43) and
(5.44) yield

32
[]
=1+
, =
,
(5.46)
1
1 2
2 + 21

1.5 Micromechanics of Heterogeneous Media

129

in the scalar case, and

k
2 [k]
=1+
,
k1
k1 + 1 1 [k]

2 []
=1+
,
1
1 + 1 1 []

(5.47)

in the elastic case; the parameter 1 and 1 are defined in Eq. (4.58).
The formula (5.46) reproduces the Maxwell one, cf. Eq. (1.19),
which in turn coincides with the appropriate Hashin-Shtrikman
bound (3.34).
Similarly to the scalar case, the formulae (5.47) provide bounds on
the elastic moduli of a two-phase material in the well-ordered case
(see Section 1.3.5). Indeed, the expression for the bulk modulus in
(5.47) obviously coincides with one of the Hashin-Shtrikman bounds
on k , as given in (3.68). Similarly, the expression for in the same
formula (5.47) is the HS-bound [HAm], which is upper, if k2 > k1 ,
2 > 1 and lower, if both inequality signs are reversed.
It is noted that the formulae (5.47) were also derived by Kuster
and Toksoz in the paper [KUc] (which is most often cited in geophysical community, but almost unknown outside it). The method
employed there is a direct generalization of the Maxwell sphere,
treated in the scalar context in Section [Link]. Here the long-wavelength scattering by a big sphere, containing a host of small inhomogeneities, was evaluated in two different ways: (i) taking the
spheres as homogeneous with unknown effective moduli and, (ii) summing the scattering waves from all small particles. Equating these
quantities reproduces (for spherical shape of the latter) the formula
(5.47), see also Berrymans discussion in [BEm], [BEn].
For a fiber reinforced materialcircular cylinders, aligned
along the axis Ox3 , the tensors Aw and Aw , as given in Eqs. (4.15)
and (4.62), respectively, are utilized.
In the scalar case the effective conductivity tensor is transversely
isotropic and has the form (5.5). The transverse conductivity reads
11
22
=1+
,
1
1 2

[]
,
2 + 1

(5.48)

which is the 2-D counterpart of the Maxwell formula (1.19). Note that
the latter also provides a bound on for an arbitrary macroscopically
isotropic 2-D mixture.

130

Heterogeneous Media: Modelling and Simulation

For the conductivity along the fibers we have the same formula as
in (5.6), i.e. the Voigt value 33 = = v (and this is fully natural,
because the temperature gradient along the fibers is homogeneous).
In the elastic case the fiber-reinforced material under study is
again transversely isotropic and thus it is characterized macroscopically by means of five elastic moduli. Utilizing Hills notations (4.60)
and (4.61), with Eq. (4.62) for the tensor Aw inserted into the formula
(5.44) for L , we get, after some algebra:

L = (2k , ` , ` , n , 2m , 2 ) ,
k = 1 + 1 +

2 1
([] + []) ,

` = 1 + 1

n = 1 + 2 [] 1 2

m = 1

= 1

[]2
,

2 []
1+
1 + 1 1 []

22 []
1+
21 + 1 []

2 1 []
,

(5.49)

with the notations

= 1 (2 + 2 ) + 2 (1 + 1 ) + 1 ,
1 =

3 41
1 + 31
=
.
2(1 + 21 )
4(1 1 )

(5.50)

The above expressions (5.49) for the moduli k , ` , n , m and

coincide with those found by Hashin and Rosen [HAn], see also [HIc].
In their derivation, Hashin and Rosen used the exact 2-D counterpart
of the Hashin assemblage (see Section 1.3.4)i.e. a medium, consisting of parallel composite cylindrical elements (a cylindrical core of
one of the constituents, coated with a concentric shell of the other,
in such a way that the volume fractions 1 and 2 are fixed). They
also showed that the value of , as given in Eq. (5.49), lies within
the rigorous bounds, derived in the same paper [HAn]. The formulae (5.49) were rederived by Levin [LEb] as well, who employed, in
a quite complicated way, a version of the effective field method, see
also [KAf].

1.5 Micromechanics of Heterogeneous Media

131

In a similar manner one can treat a medium, containing identical, but randomly oriented spheroidal inhomogeneities. The scalar
conductivity of such a medium, according to Eqs. (4.23) and (5.43),
reads
2 []
= 1 +
,
(5.51)
1 + 2
where depends on the shape of the spheroid through its depolarization factors M and M , as seen from the formula (4.23).
Remark 5.14 The formula (5.51) in an equivalent form, was proposed by Fricke [FRa]. In this work, and in a series of papers that
followed, the author compared some accurate experimental results
with the prediction of (5.51), using data for conductivity of dogs
blood. The red corpuscles were approximated as identical, but randomly oriented spheroids, and Fricke found that the best fit to experiments corresponds to the choice of prolate spheroidal form, when
(in our notations) c/a = 4.25. These papers of Fricke are in a sense
certain natural extension of Einsteins idea, already discussed in Section [Link]: while the latter looked only for the size of the particles
in a solute by means of micromechanical experiment, the former was
interested in more specific details, like the particle shape (assuming
it spheroidal which, of course, is by no means obvious). A number
of new methods of such micromechanical nature, that are able to
provide information about the microstructure of a dispersion (connected, e.g. with sedimentation, diffusion, etc.), were developed later
on. Some ideas, results and references can be found in [JEb], where
polydisperse materials were treated by means of the notion of the
so-called effective diameter.
The elastic moduli of the same dispersion of spheroids are
k = k1 +

2 [k]a0
,
1 + 2 a0

= 1 +

2 []a00
,
1 + 2 a00

(5.52)

see (5.44) and (4.66), with the parameters a0 and a00 , defined in (4.67).
Consider finally the most interesting particular case of a microcracked solid, when the spheroids degenerate into randomly oriented
penny-shaped cracks (with their conductivity and elastic moduli vanishing at the same time).

132

Heterogeneous Media: Modelling and Simulation

In the scalar case Eq. (5.51) yields in this limiting case

1

=
,
1
1+

= lim 2 ,
0

where = c/a, see the end of Section 1.4. The quantity in the last
formula was evaluated in Eq. (5.10) and the result, let us recall, is
= 89 . Hence

1
=
,
(5.53)
1
1 + 89
where is the crack density parameter, cf. (5.11).
Remark 5.15 As shown recently by the author [MAf], the approximation (5.53) provides an upper bound on the effective conductivity
of the solid, if (i) there is no correlation between the location and the
orientation of the cracks and, (ii) The distance between the centers
of the cracks is greater than 2a (which means that overcrossing of
is forbidden, whatever the orientations, with cracks centers fixed).
We recall that a denotes the radius of the penny-shaped cracks under
consideration.
For an elastic microcracked solid the formula (5.52) yields
k
2 a0
=1
,
k1
1 + 2 a0

2 a00
=1
,
1
1 + 2 a00

in which the limit = c/a 0 should be taken once again. But

the appropriate limits were already evaluated in the dilute case,
cf. Eq. (5.17). Hence, the effective field predictions for the elastic
moduli of a microcracked solid read
k
=
k1

1
,
16 1 12
1+

9 1 21

=
1

1
.
32 (1 1 )(5 1 )
1+

45
2 1

(5.54)

The formulae (5.54) were proposed by Benveniste [BEd] within

the frame of Mori-Tanakas approximation (see also [BEc] for the 2-D

1.5 Micromechanics of Heterogeneous Media

133

casea microcracked plane). A thorough discussion of this approximation with many appropriate references, can be found in Kachanovs
survey [KAa]. Note that the foregoing coincidence of the Benveniste
results with Eq. (5.54) is fully natural, having in mind the already
pointed out equivalence of the Mori-Tanakas scheme and the effective field idea, in its simplest form treated here. It turned out also
that the values (5.54), similarly to the scalar case, cf. Remark 5.15,
provide upper bounds on the effective elastic moduli of the microcracked solid, under the conditions (i) and (ii) of the same remark, see
[MAf].
1.5.6 The Conductivity of a Medium with Nonideal
Interphase Contacts
In the case of thermally nonideal interphase contacts, the one-particle approximations can be constructed in a manner, fully similar
to the ideal contacts case, treated in the foregoing analysis. The only
difference is that the formulae (2.50) and (2.53) should be employed
instead of (2.41), thus taking into account the surface concentration
factors A(12) and B (12) , as defined in Eqs. (2.49) and (2.52).
The simplest case of spherical inclusions will be treated again for
illustration, within the frame of the self-consistent approach. The
formulae (4.38) and (4.44), valid for a single sphere, allow us to obtain
then the surface concentration factors

A(12) = 2 A(12)
s ( , 2 ) ,

B (12) = 2 B (12)
s ( , 2 ) ,

having imagined that each sphere lies in a matrix with unknown

conductivity , cf. (2.49) and (2.52). Recalling (4.13) and inserting
the explicit expressions for A(12)
and B (12)
from the same formulae
s
s
in (2.50) and (2.53), respectively, gives the self-consistent equations

= 1 + 32

[]
2C
+
2 + 2 2 + 2(1 + C)

(5.55)

in the superconducting case, and

= 1 + 32
in the resistance case.

2R1
[]

2 + 2
2 + 2(1 + R)

(5.56)

134

Heterogeneous Media: Modelling and Simulation

With the same ease the rest of the above discussed approximate
schemes of one-particle type can be implemented. Some other examples of similar schemes with more detailed exposition, discussion,
comparison with experiments, etc., can be found, e.g. in the papers
[HAo], [HAp], [BEb], [MIa] et al. Variational principles and estimates, accounting for nonideal contacts, are given, for instance, in
the papers [TOd], [LIa], which also include a number of appropriate
and useful references.
Here it suffices perhaps to demonstrate briefly only two specific
effects in two-phase media that stem out from the nonideal interphase
contacts.
First, imagine that the inclusions are highly conducting, =
2 /1 , but the thermal resistance of the contacts is very high
at the same time, R , with the ratio = R/ remaining finite.
The formula (5.56) then gives
1

=
.
1
1 32 (1 )
This relation, for small to moderate values of 2 , is in a reasonable
agreement with experimental data of De Araujo and Rosenberg [DEa]
and the appropriate lower bound of Torquato and Rintoul [TOd].
Second, assume that particles of higher conductivity are distributed in a matrix, i.e. > 1. Then
= 1 ,

if R = Rc = 16 ( 1)( + 2) ,

(5.57)

as it follows again from (5.56). The above means that at R = Rc

the thermal resistance of the interphases hides the inclusions, so
that macroscopically their presence is not felt at all. This fact is directly connected with the presence of the so-called critical radius
ac , see [CHc], [EVa] et al. The essence is that the interfacial effects
are size-dependent: the radius a enters the definitions of both interfacial characteristics C and R, cf. their definitions (4.26) and (4.28)
(or (4.27)) respectively. In particular, whatever be the fixed volume
fraction 2 , Eq. (4.28) implies that when decreasing the particles radius a, we shall reach a value a = ac the critical radiusfor which
R coincides with the value Rc , given in (5.57). Then the bad interphase will produce enough resistance to compensate for the higher
conductivity of the inclusions. The effective conductivity, as a result,
will not be able to exceed that of the matrix, when a ac .

1.5 Micromechanics of Heterogeneous Media

135

Torquato and Rintoul [TOd] have shown that the value of Rc =

1 exactly for a dispersion of nonoverlapping spheres. The foregoing self-consistent scheme gave a totally different value of Rc (see
Eq. (5.57)). This fact demonstrates once again that the self-consistent
reasoning, though based on consideration of a single inclusion in a homogeneous medium, has, in reality, little in common with dispersions
of equi-sized particles.
Note finally that in the papers [KAd], [KAe] the thermoelastic
properties of composites (dispersion of spheres) with nonideal contacts were studied in detail by means of the effective field method.
1.5.7 The Self-consistent Scheme for the Absorption
Problem
To illustrate the generality of the one-particle approximations, we
shall discuss here a self-consistent type scheme for evaluating the
effective sink strength. Recall that this problem, motivated in particular by the classical Smoluchowskis work [SMa], was discussed in
Section 1.2.7. It is specified by Eqs. (2.66) and (2.67). The scheme
to be applied is due to Brailsford and Bullough [BRa]. The implementation follows Talbot and Willis [TAb].
The self-consistent scheme provides an approximation for hci2 assuming that each inclusion is single and it is immersed into a medium
of the effective sink strength k 2 . Hence, choosing for simplicity the
inclusions shape spherical (of the radius a), one should solve the
single inclusion problem which, in the present context reads:
c(x) k22 c(x) + K = 0 ,
c(x) k 2 c(x) + K = 0 ,

if r < a ,
(5.58)
if r > a .

Due to the radial symmetry of the problem, the solution has the
form

K
A sinh k2 r

,
r < a,

k2 +
k2 r
2
c(x) =

K
Bek r

2 +
,
r > a.
k r
k2
The constants A and B are fixed by the conditions that both the
function c(x) and its normal derivative dc/dr on the sphere r = a

Heterogeneous Media: Modelling and Simulation

136

should be continuous. The calculations, though tedious, are straightforward, and the final result for c(x), together with (2.73), yields
eventually
k 2 k12
(5.59)
k22 k12
2
= 2
k2

(k 2 k 2 )(1 + k a)(k2 a cosh k2 a sinh k2 a)

+3 2
(k2 a)2 (k a sinh k2 a + k2 a cosh k2 a)

This cubic (with respect to the dimensionless quantity k a) equation

is the absorption counterpart of the self-consistent equations (5.19)
and (5.28) for the conductivity and the elastic moduli, respectively.
Similarly to them, Eq. (5.59) does not pretend to be exact, since
spheres interactions are not taken into account with due precision,
but only approximately. Nevertheless, it allows to reveal certain interesting features of the absorption problems homogenization, that
distinguish the latter from that of the more ordinary scalar and
elastic problems.
Consider first the limiting case, when
k1 a 0 ,

k2 a 0 .

(5.60)

The solution of Eq. (5.59) is then k 2 = kv2 = 1 k12 + 2 k22 . This

is just the Voigt type approximation, already observed in the scalar
and elastic cases, and also in the absorption problem under study, see
Eq. (3.78). Recall that it corresponds to the assumption that the concentration field c(x) remains constant in V. It is noted, after Talbot
and Willis [TAb], that the formula (3.78) is fully natural in the case
under study. Indeed, the limits (5.60) mean that, at fixed k1 and k2 ,
the radius of the inclusions decreases to zero. The fluctuation of the
field k 2 (x), as defined in (2.68), becomes then more and more rapid,
c(x) has, so to say, little chance to vary and remains approximately
constant. A rigorous justification of the foregoing reasoning, i.e. the
validity of the Voigts type formula (3.78) in the case (5.60) is due to
Papanicolaou [PAa].
In the other limiting case, when
k1 a ,

k2 a ,

(5.61)

Eq. (5.59) generates now the Reuss type formula k 2 = kr2 , see
Eq. (3.79). The physical explanation of this result is also transparent

1.5 Micromechanics of Heterogeneous Media

137

enough [TAb]. The limit (5.61) means that at fixed k12 and k22 the
spheres radius increases. The matrix-inclusion interrelations affect
then only thin boundary layers, elsewhere c K/k12 or c K/k22
within the phase 1 and 12 respectively. Thus
c=

K
1 K
2 K
= 2 + 2 ,
k 2
k1
k2

which does yield (3.79).

The foregoing results indicate that the elementary bounds (3.80)
are sharp, in the sense that they cannot be improved.
Hence in both limiting cases Eq. (5.59) is exact. The most interesting range is therefore k1 a, k1 a O(1), since the interactions
then become of crucial importance, and the exact spatial distribution of the sinks cannot be neglected. We shall illustrate this fact,
again qualitatively, on the classical Smoluchowskis problem (see Section 1.1.3) that corresponds to the hard limit (2.69) of Eq. (2.66).
In this case the self-consistent equation (5.59) becomes
(1 1 )(k a)2 = 31 (1 + k a)
whose elementary solution, when expanded as a series for small values
of 2 , reads
3/2
k 2 a2 = 32 + 3 32 + ,
or

p
k 2
=
1
+
32 + 52 2 + ,
ks2

(5.62)

where ks2 is the Smoluchowskis value of the effective absorption, valid

in the dilute limit, see Eq. (1.24) (at D = 1).
The formula (5.62) demonstrates an interesting and important feature of the Smoluchowski problem, namely, that the effective absorption coefficient k 2 represents a nonanalytic function of the concentration 2 of the trapping inclusions. (The expansion (5.62) suggests,
but does not
prove of course, that instead it maybe an analytical
function of 2 .) This fact sharply contrasts the scalar and elasticity counterparts, as well as the absorption problem Eq. (2.66) for any
other case, different from the hard limit.
The explanation of the nonanalytic dependence (5.62) lies in the
appearance of the so-called screening (or shielding) effect which

138

Heterogeneous Media: Modelling and Simulation

becomes of predominant importance in the Smoluchowski problem.

Roughly speaking, this means that a trap will absorb the particles
close to it, so that fewer of them can reach other traps behind to
be eaten by them. Such a screening effect is very important also
in systems of charged particles (the Debye shielding), see e.g. [VAa,
Chapter 8.1].
The first several terms of the expansion (5.62) have been obtained,
bounded by means of variational procedures and discussed by a number of authors [FEa], [TAa], [TOc], see also [CAa] and Torquatos
survey in the present volume. An interesting approach to the problem is the so-called point-particle approximation in which the traps
are replaced by appropriate point-sources whose intensity is looked
for from a natural consistency condition [FEa]. The same approach
is also useful for particulate media in studying their effective elastic
properties; see [KAf, Chapter 5.14].
A final remark in this Section is warranted. Everywhere in the
foregoing considerations we have considered 3-D heterogeneous media. The 2-D case can be treated as wella natural interpretation
will be provided by a medium containing an array of infinite cylindrical inclusions, parallel to a fixed axis, when the applied temperature
gradient or external loading is perpendicular to it. In the conductivity
and elastic contexts, all 2-D counterparts of the 3-D considerations
and results can be repeated straightforwardly which explains their
omission here.
There exists however a notable exception when the 2-D situation
is drastically different from its 3-D analog. It concerns the Smoluchowski problem in 2-D, when the inclusions represents a sparse system of aligned circular cylinders. In this case the solution of the
Laplace equation (1.23) in 2-D contains the function ln r/a, instead
of a/r, and hence the condition c(r) c0 as r cannot be satisfied. (Note that the situation here is fully similar to the well-known
Stokes paradox [HAa, Chapter 2.7], which assures us that there
exists no Stokesian flow, bounded at infinity, around a cylinder.)
The nonexistence of the single-cylinder solution in the case under
study does not mean of course that the effective absorption coefficient
k 2 cannot be defined. It only implies that the interactions are so
strong here, that a dilute formula, like the Smoluchowski one (1.25),

1.6 Micromechanics of Heterogeneous Media

139

cannot be found.16 The dependence k 2 upon 2 will be proportional,

say, to ln 2 even for 2  1. (Such logarithmic dependencies do
appear in homogenization problems with very strong interactionsan
example is provided by the very interesting and important problem,
concerned with the effective permeability of a randomly perforated
membrane, which however cannot be discussed here even briefly, see,
e.g. [MAc] and the references therein.

1.6 Elastic Properties of Polycrystals

We shall present here some basic ideas and simplest results concerning
the application of one of the approximate schemes of the previous
Section for predicting the effective properties of polycrystals. This
is the oldest scheme in the fieldthe self-consistentproposed by
Hershey [HEb] and Kroner [KRd]. As an illustration only the cubic
symmetry case will be treated in more detail.
1.6.1 The Self-Consistent Scheme for Polycrystals
The polycrystal is an assembly of monocrystalshomogeneous grains
with one and the same elastic properties, defined by the elastic tensor
L. The crystallographic axes of each grain vary. Hence the tensor
L is rotated in a complicated manner when one moves across the
solid and exactly this is what makes the polycrystal heterogeneous.
for simplicity sake and to avoid cumbersome details that may eclipse
the clarity of the basic ideas, we assume that there exist no preferable
orientations of grains, i.e. of the crystallographic axes. Hence there
is no texture presented so that the polycrystal is macroscopically
isotropic, with a tensor of effective moduli

L = 3k I0 + 2 I00 ,
where I0 and I00 are the two basic fourth-rank isotropic tensors, defined
in Eq. (4.58). Our aim is to develop a certain simple approximate
16

This explains the failure in 2-D of certain variational procedures for estimating
k for random dispersions; see the recent authors lecture [MAg] for details and
references.
2

140

Heterogeneous Media: Modelling and Simulation

scheme of self-consistent type for evaluating L by means of the given

elastic tensor L for a single grain.
Imagine each grain is a sphere, immersed into an unbounded matrix with the effective, but yet unknown, properties L . Fix one of
the grains, W. According to Eqs. (4.10) and (4.55) the strain within
such a grain is constant
gr = Aw (L, L ) : E ,

(6.1)

where E = is the prescribed macrostrain tensor, applied to the

polycrystalline RVE. Moreover, in virtue of Eqs. (4.7), (4.55) and
(4.58)
A(L, L ) = A = [I + P : (L L )]1 ,

P = p1 I0 + p2 I00 ,
p1 =

3k

1
,
+ 4

(6.2)
p2 =

3
+
.

5 3k + 4

(The subscript w will be suppressed from now on.) The asterisk

everywhere underlines that the appropriate quantities, like p1 , etc.,
are evaluated by means of the effective elastic moduli k , of the
polycrystal.
Note that in writing Eqs. (6.1) and (6.2) we have used the fact
the the formula (4.55) is applicable under the sole assumption that
the matrix is isotropic; the anisotropy of the inclusion can be at the
same time arbitrary, as pointed out in Section 1.4.
Let us now average Eq. (6.1) with respect to all possible crystallographic orientations of the axes of the grain. (This operation will
be denoted by hi , similarly to Section 1.4). Then
hgr i = hA i : ,

(6.3)

The key observation now is the identity

hgr i =

(6.4)

i.e. the average (with respect to all crystallographic orientations)

strain in a grain equals the macrostrain.
Remark 6.1 Assuming all grains spherical we have a kind of Hashins
assemblage (Section 1.3.3)a set of spheres of different sizes, from

1.6 Micromechanics of Heterogeneous Media

141

finite down to infinitesimal, filling the RVE, V. The spheres differ

(mechanically) only in crystallographic orientation so that, averaging
over the latter, we should get the macrostrain. (Observe that the
strain-concentration factor A in Eq. (6.1) does not depend on the
grain size.)
Inserting (6.3) into (6.4), one gets the basic self-consistent equation
hA i = I ,
(6.5)
through which the unknown effective tensor L can be found; I is the
unit fourth-rank tensor that appeared in Eq. (4.55).
An equivalent to Eq. (6.5) formulation follows, if one applies the
tensor L to both parts of Eq. (6.1), evaluating thus the stress tensor
within a grain:
gr = L : gr = L : A : ,
and average the result with respect to all possible orientations:
hi = hL : gr i = hL : A i : .

(6.6)

But, similarly to Eq. (6.4), observe that

h gr i = = L :
which, together with Eq. (6.6), yields the equation

L = hL : A i ,

(6.7)

equivalent to the self-consistent one (6.5).

Remark 6.2 The equivalence of Eqs. (6.5) and (6.7) can be formally
demonstrated in a simple way, using the expression for the tensor A ,
see (6.2). Indeed, from this expression we have
[ I + P : (L L ) ] : A = I
which, when averaged over orientations, yields
hA i + P : hL : A i P : L : hA i = I ,
or

I hA i = P : hL : A i L : hA i .

(6.8)

142

Heterogeneous Media: Modelling and Simulation

Now, if Eq. (6.5) holds, i.e. hA i = I, then Eq. (6.8) immediately

implies (6.7). Conversely, if Eq. (6.7) holds, then hL : A i = L and
(6.8) becomes
h

I hA i = P : L : I hA i ,
i.e. I hA i = 0 which is exactly the needed Eq. (6.5).
Other equivalent forms of the self-consistent equation (6.5) have
been proposed by Kroner [KRe]. He then ingeniously exploited them
to bound the effective properties for a polycrystal, using simple algebraic arguments.
1.6.2 Formulation in Compliances and Consistency of the
Scheme
Let
(6.9)
gr = B : , B = Bw (M, M ) ,
be the solution of the single sphere (i.e. single grain) problem in
stresses; M = L1 and M are the compliance tensors of the grain
and of the polycrystal respectively. Recall that both stress and strain
tensors gr and gr are constant within the (spherical) grain under
study. Upon comparing Eqs. (6.1) and (6.9) we find that

A = M : B : M1

or M : B = A M .

(6.10)

We employ again the reasoning of Section 1.6.1, but with respect

to stresses now. That is, average Eq. (6.9) with respect to all grains
orientations:
h gr i = hB : (M, M )i :
and note that h gr i = . Then the self-consistent equation for the
effective compliance tensor M emerges
hM i = I ,

(6.11)

which is the natural counterpart of Eq. (6.5). An equivalent formulation reads

M = hM : B i .
(6.12)
The latter can be obtained, say, if one applies the tensor M to both
parts of Eq. (6.9) to get the strain gr = M : gr in the grain, and then
average this result with respect to all crystallographic orientations.

1.6 Micromechanics of Heterogeneous Media

143

Theorem 6.1 Eqs. (6.5) and (6.11) are equivalent, in the sense that

L : M = I ,

(6.13)

i.e. the effective properties of the polycrystal, calculated by means

of the above described self-consistent procedure, are the same, no
matter whether the schemes, based on elastic moduli or compliances,
are utilized.
Proof. Let M be the effective compliance, as calculated from (6.11).
Denote L+ = M1 and average Eq. (6.10) over all orientations:

A(L, L+ )

: M = hM : B i = M ,

having used Eq. (6.12). Thus A(L, L+ ) = I; comparing the latter

with Eq. (6.5) that specifies L , one concludes that indeed L+ = L
= M1 .

1.6.3 Example: Elastic Moduli of a Cubic Polycrystal

The tensor of the elastic moduli L in this simplest case is

L = 3k I0 + 22 I00 + 2(1 2 )Oh ,

(6.14)

Oh = e41 + e42 + e43

(6.15)

where
is the basic fourth-rank tensor with cubic symmetry (whose axes are
along the orthonormal crystallographic basis ei , i = 1, 2, 3, see [SPa]).
For this simplest anisotropy, the effective bulk modulus coincides with
that of the grains:
k = k + 23 (1 2 ) ,
(6.16)
because if the applied stress tensor is spherical. the strain is a pure
dilatation, one and the same in all grains. To simplify the calculation
of the shear modulus, introduce three basic fourth-rank tensor with
cubic symmetry as follows
1 = I0 ,

2 = I0 + I00 Oh ,

2 = Oh I0 .

(6.17)

It is easily checked that they are orthogonal, in the sense that

i : j = 0 ,

i 6= j ,
(6.18)

i : j = i ,

i = 1, 2, 3 (no sum over i) .

144

Heterogeneous Media: Modelling and Simulation

In the sequel, the symbolic notation

L = (, , ) L = 1 + 2 + 3

(6.19)

will be useful. The isotropic tensors are characterized by the equality

= , i.e. they have the form (, , ). The most important formula, however, concerns the tensors of the form (6.19), averaged over
orientation; it reads, symbolically,
= 15 (3 + 2) .

hLi = h(, , )i = (, , ) ,

(6.20)

It is a straightforward consequence of the fact that hOh i = 15 H,

where H is the fully-symmetric isotropic fourth-rank tensor, defined
in (4.56); the multiplier 1/5 follows from a full contraction of both
sides of (6.20), having noticed that hOh i is proportional to H. Hence
hOh i =

1
5

2I + 3I0 = I0 + 25 I00 .


According to the definitions (6.17) of the basic tensors 2 and 3 ,

we have
h2 i = 35 I00 , h3 i = 25 I00 ,
and (6.20) immediately follows from the last two equations.
For the tensor L of a single grain, see (6.14), we have

L = (, , ) ,

= 3k , = 22 , = 21 .

(6.21)

Due to the orthogonal properties (6.18) of the basic tensors

i , i = 1, 2, 3, the multiplication and inversion of the tensors in the
symbolic form (6.19) is extremely simple. Namely,

L = (, , ) L1 = (1/, 1/, 1/) ,

L0 = (0 , 0 , 0 ) ,

L00 = (00 , 00 , 00 ) L0 : L00 = (0 00 , 0 00 , 0 00 ) .

Due to Eqs. (6.2) and (6.21) we thus have

I + P : (L L )
h

= (1, 1, 1) + (p1 , p2 , p2 ) (3k , 22 , 21 ) (3k , 2 , 2 )

= (1, 1 + 2p2 (2 ), 1 + 2p2 (1 ) ,

1.6 Micromechanics of Heterogeneous Media

145

so that,
hAi =

h

I + P : (L L )

i1 

3
1
2
1
=
+
.

5 1 + 2p2 (2 ) 5 1 + 2p2 (1 )
Here
2p2 =

3 k + 2
5 3k + 4

(6.22)

(6.23)

in virtue of Eq. (6.2) and k is the bulk modulus of the polycrystal,

given in (6.16). It remains to recall the self-consistent condition (6.5)
hA i = I = (1, 1, 1) which means that the unknown shear modulus
solves the equation = 1, with defined in (6.22). Using (6.23)
we recast this equation, after some algebra, to its final form:


9k + 41
a=
,
81

3

+a
1


2

b
1


32 (k + 41 )
b=
,
821

c = 0,
3k 2
c=
.
421

(6.24)

It is easy to see that if k, 1 , 2 > 0, this cubic equation has only one
real root, as it should.
Eq. (6.24) was first derived by Hershey [HEb] and, independently,
by Kroner [KRd], who developed the self-consistent theory of polycrystals. Details and further references can be found in the reviews
[WAb], [KRf], as well as in the book [SHb] (concerned mainly with
polycrystals and covering, in particular, the rich Russian literature
on the subject).
The experimental data reported and discussed in these references
support very favorably the predictions of the self-consistent theory for
polycrystals and, in particular, those of Eq. (6.24). This is natural
enough, since for polycrystals we have a well-defined basic element
to be embedded in an effective medium, namely the single crystallite.
For a particulate medium the choice of such an element is not that
obvious, and many choices may compete, as already pointed out in
Section [Link] (see Remark 5.8).

146

Heterogeneous Media: Modelling and Simulation

Acknowledgements
The author cordially thanks Z. Hashin, S. Kanaun, D.R.S. Talbot,
J. Willis and R. Zimmerman for a number of very useful critical remarks and suggestions when reading a draft of this survey. Many colleagues and friends decisively helped me with reprints of their works
and copies of the 19th century works, locally unavailable. Z. Hashin
sent me his rare and highly valuable 1972 NASA report [HAh], C. Trimarco provided a copy of original Mossottis papers, P. Germain and
D. Lhuilliercopies of the Poissons Memoires on magnetism from the
Archive of the French Academy. Special thanks are due to D. Jeulin
who turned authors attention to the very interesting Matherons paper [MAi]. The support of Bulgarian Ministry of Education and
Science under Grant No. MM 805/1998 is acknowledged as well.
The final form of this survey was completed during the Workshop
on Mathematical Developments in Modelling Microstructure and
Phase Transformation in Solids (September, 1999), organized by the
Isaac Newton Institute for Mathematical Sciences at the University of
Cambridge. The hospitality and the most stimulating atmosphere of
the Institute, together with discussions with many of the participants,
are gratefully acknowledged.

1.7 References
[ARa] Archie, G. E., The electrical resistivity log as an aid in determining some reservoir characteristics, Trans. AIME (Am.
Inst. Min. Metall. Eng.), 146 (1942), 5467.
[ARb] Arnold, D. H., The mecanique physique of Simeone Denis
Poisson: The evolution and isolation in France of his approach
to physical theory. VI. Elasticity: The crystallization of Poissons views on the nature of matter, Archive History Exact
Sci., 28 (1983), 343367.
[AVa] Avellaneda, M., Iterated homogenization, differential effective
medium theory and applications, Commun. Pure Appl. Math.,
40 (1987), 527554.

1.7 Micromechanics of Heterogeneous Media

147

[BAa] Baker, G. A., Jr., and Graves-Morris, P., Pad

e Approximants. Part I. Basic Theory. Part II. Extensions and
Applications, Addison-Wesley, Reading, Massachusetts
(1981).
[BEa] Bensoussan, A., Lions, J. L., and Papanicolaou, G., Asymptotic Analysis for Periodic Structures, North-Holland
(1978).
[BEb] Benveniste, Y., The effective mechanical behavior of composite materials with imperfect contact between the constituents,
Mech. Materials, 4 (1985), 197208.
[BEc] Benveniste, Y., On the Mori-Tanakas method in cracked bodies, Mech. Res. Commun., 13 (1986), 193201.
[BEd] Benveniste, Y., A new approach to the application of the MoriTanakas theory in composite materials, Mech. Materials, 6
(1987), 147157.
[BEe] Benveniste, Y., Dvorak, G. J., and Chen, T., On diagonal and
elastic symmetry of the approximate effective stiffness tensor
of heterogeneous media, J. Mech. Phys. Solids, 39 (1991), 927
946.
[BEf] Benveniste, Y., and Miloh, T., The effective conductivity of
composites with imperfect contact at constituent interfaces,
Int. J. Eng. Science, 24 (1986), 15371552.
[BEg] Beran, M., Use of a variational approach to determine bounds
for the effective permittivity of a random medium, Nuovo Cimento, 38 (1965), 771782.
[BEh] Beran, M., Bounds on field fluctuations in a random medium,
J. Appl. Phys., 39 (1968), 57125714.
[BEi] Beran, M., Statistical Continuum Theories, John Wiley,
New York (1968).
[BEj] Beran, M., and Molyneux, J., Use of classical variational principles to determine bounds for the effective bulk modulus in
heterogeneous medium, Q. Appl. Math., 24 (1966), 107118.
[BEk] Bergman, D. J., The dielectric constant of a composite materiala problem in classical physics. Phys. Reports, 43C
(1978), 377407.
[BEl] Bergman, D. J., Analytical properties of the complex dielectric constant of a composite medium with applications to the

148

Heterogeneous Media: Modelling and Simulation

derivation of rigorous bounds and to percolation problems, in
Electrical transport and optical properties of inhomogeneous media, AIP Conf. Proc. No. 40, Garland, J. C. and
Tanner, D. B., eds., Am. Inst. Phys., New York (1978), 4661.

[BEm] Berryman, J. G., Long-wavelength propagation in composite

elastic media. I. Spherical inclusions. II. Ellipsoidal inclusions,
J. Acoust. Soc. Am., 68 (1980), 18091819, 18201831.
[BEn] Berryman, J. G., Single-scattering approximations for coefficients in Biots equations of poroelasticity, J. Acoust. Soc.
Am., 91 (1992), 551571.
[BEo] Berryman, J. G., and Berge, P. A., Critique of two explicit
schemes for estimating elastic properties of multiphase composites, Mech. Materials, 22 (1996), 149164.
[BOa] Bobeth, M., and Diener, G., Field fluctuations in multicomponent mixtures. J. Mech. Phys. Solids, 34 117, 1986.
[BRa] Brailsford, A. D., and Bullough, R., The theory of sink
strengths, Philos. Trans. Roy. Soc. Lond. A, 302 (1981), 87
137.
[BRb] Brinkman, H. C., The viscosity of concentrated suspensions
and solutions, J. Chem. Phys., 20 (1952), 571.
[BRc] Bristow, J. R., Microcracks, and the static and dynamic elastic
constants of annealed and heavily cold-worked metals, British
J. Appl. Physics, 11 (1960), 8185.
[BRd] Brown, W. F., Solid mixture permittivities, J. Chem. Phys.,
23 (1955), 15141517.
[BRe] Brown, W. F., Dielectrics, in Encyclopedia of Physics,
vol. XII, Fl
ugge, D., ed., Springer-Verlag, Berlin-GotingenHeidelberg (1956), 1154.
[BRf] Bruggeman, D. A. G., Berechnung verschiedener physikalischer Konstanten von heterogenen Substanzen. I. Dielectrizitatkonstanten und Leitfahigkeiten der Mischkorper aus isotropen
Substanzen, Ann. Physik Leipzig, 24 (1935), 636679.
[BRg] Bruggeman, D. A. G., Berechnung verschiedener physikalischer Konstanten von heterogenen Substanzen. III. Die Elastische Konstanten der Quasiisotropen Mischkorper aus isotropen
Substanzen, Ann. Physik Leipzig, 29 (1937), 160178.

1.7 Micromechanics of Heterogeneous Media

149

[BUa] Budiansky, B., On the elastic moduli of some heterogeneous

materials, J. Mech. Phys. Solids, 13 (1965), 223227.
[BUb] Budiansky, B., and OConnell, R., Elastic moduli of a cracked
solid, Int. J. Solids Structures, 12 (1976), 8197.
[CAa] Calef, D. F., and Deutch, J. M., Diffusion-controlled reactions,
Ann. Rev. Phys. Chem., 34 (1983), 493523.
[CHa] Chaban, L. A., Self-consistent field approach to calculation
of effective parameters of microinhomogeneous media, Sov.
Phys.-Acous., 10 (1964), 298304.
[CHb] Cherkaev, A. V., Lurie, K. A., and Milton, G. W., Invariant properties of the stress in plane elasticity and equivalence
classes of composites, Proc. Roy. Soc. A, 438 (1992), 519529.
[CHc] Chiew, Y. C., and Glandt, E. D., Effective conductivity of dispersions: The effect of resistance at particle interfaces, Chem.
Eng. Sci., 42 (1987), 26772685.
[CHd] Christensen, R. C., Mechanics of Composite Materials,
John Wiley, New York (1979).
[CHe] Christensen, R. C., and Loo, L. K., Solution for effective shear
properties in three phase and cylinder models, J. Mech. Phys.
Solids, 27 (1979), 315330.
[CHf] Christiansson, H., and Helsing, J., Poissons ratio of fiberreinforced composites, J. Appl. Phys., 79 (1996), 75827585.
[CLa] Clausius, R., Die mechanische Behandlung der Elektricit
at, Vieweg, Braunshweig (1879).
[CLb] Cleary, M. P., Chen, I. W., and Lee, S. M., Self-consistent
techniques for heterogeneous solids, ASCE Journal Eng. Mechanics, 106 (1987), 861887.
[COa] Corson, P. B., Correlation functions for predicting properties
of heterogeneous materials. IV. Effective thermal conductivity
of two-phase solids, J. Appl. Phys., 45 (1974), 31803182.
[DAa] Darcy, H., Les Fountaines Publiques de la Ville de Dijon, Victor Dalmont, Paris (1856).
[DAb] Davies, W. E. A., The theory of elastic composite materials,
J. Phys. D, 4 (1971), 13251339.

150

Heterogeneous Media: Modelling and Simulation

[DEa] de Araujo, F. F.T., and Rosenberg, H. M., The thermal conductivity of epoxy-resin/metal-powder composite materials
from 1.7 to 300 K, J. Phys. D, 9 (1976), 665675.
[DEb] Dewey, J. M., The elastic constants of materials loaded with
non-rigid fillers, J. Appl. Phys., 18 (1947), 578586.
[DRa] Drugan, W. J., and Willis, J. R., A micromechanics-based
nonlocal constitutive equation and estimates of representative
volume element size for elastic composites, J. Mech. Phys.
Solids, 44 (1996), 497524.
[EBa] Ebert, L. J., and Wright, P. K., Mechanical aspects of interphases, in Composite Materials, Vol. 1, Metcalfe, A. G.
ed., Academic Press, New York (1974), 4276.
[EIa] Einstein, A., Eine neue Bestimmung der Molek
uldimensionen,
Ann. Physik, 19 (1905), 289306.
[ESa] Eshelby, J. D., The determination of the elastic field of an
ellipsoidal inclusion, and related problems, Proc. Roy. Soc.
London A, 241 (1957), 376396.
[EVa] Every, A. G., Tzou, Y., Hasselman, D. P. H., and Raj, R., The
effect of particle size on the thermal conductivity of ZnS/Diamond composites, Acta Metall. Mater., 40 (1992), 123129.
[FAa] Faraday, M., Experimental researches on electricity, Phil.
Trans. Roy. Soc. London, II Ser. (1838), page 1ff.
[FEa] Felderhof, B. U., and Deutch, J. M., Concentration dependence of the rate of diffusion-controlled reactions, J. Chem.
Phys., 64 (1976), 45514558.
[FEb] Ferrari, M., Asymmetry and the high concentration limit of
the mori-tanaka effective medium theory, Mech. Materials,
11:251256, 1991.
[FEc] Feynman, R. P., Leighton, R. B., and Sands, M., The Feynman lectures on physics, vol. 2, Addison-Wesley, Reading,
Massathussets (1964).
[FOa] Foldy, L. L., The multiple scattering of waves, Phys. Rev., 67
(1945), 107119.
[FRa] Fricke, H., A mathematical treatment of the electric conductivity and capacity of disperse systems. I. The electric conductivity and capacity of disperse systems, Phys. Revue, 24
(1924) 575587.

1.7 Micromechanics of Heterogeneous Media

151

[FRb] Frisch, H. L., Statistics of random media, Trans. Soc. Rheology, 9 (1965), 293312.
[FRc] Frisch, U., Wave propagation in random media, in Probabilistic Methods in Applied Mathematics, vol. 1,
Bharucha-Reid, A. T., ed., Academic Press, New York and
London, (1968), 75198.
[FRd] Frolich, H., and Sack, R., Theory of the rheological properties
of dispersions, Proc. R. Soc. London A, 185 (1946), 415430.
[GAa] Gambin, B., and Kroner, E., Convergence problems in the
theory of random elastic media, Int J. Eng. Sci., 19 (1981),
313318.
[GOa] Golden, K., and Papanicolaou, G., Bounds for effective properties of heterogeneous media by analytic continuation, Comm.
Math. Phys., 90 (1983), 473491.
[GRa] Grabovsky, Yu., and Milton, G. W., Exact relations for composites. Towards a complete solution, Documenta Math., Extra Volume ICM 1998 (1998), 110.
[GUa] Gubernatis, J. E., and Krumhansl, J. A., Microscopic engineering properties of polycrystalline materials: Elastic properties, J. Appl. Phys., 46 (1975), 18751883.
[HAa] Happel, J., and Brenner, H., Low Reynolds Number Hydrodynamics, Prentice-Hall (1965).
[HAb] Hartree, D. R., The wave mechanics of an atom with noncoulomb central field. II. Some results and discussion, Proc.
Camb. Phil. Soc., 24 (1928), 111ff.
[HAc] Hashin, Z., The moduli of an elastic solid reinforced by rigid
particles, Bull. Res. Counc. Israel, 5C (1955), 4659.
[HAd] Hashin, Z., The moduli of an elastic solid, containing spherical
particles of another elastic material, in Nonhomogeneity in
Elasticity and Plasticity, Olszak, W., ed., Pergamon Press,
Oxford (1959), 463478.
[HAe] Hashin, Z., The elastic moduli of heterogeneous materials,
ASME J. Appl. Mech., 29 (1962), 143150.
[HAf] Hashin, Z., Theory of mechanical behavior of heterogeneous
solids, Appl. Mech. Reviews, 17 (1963), 19.

152

Heterogeneous Media: Modelling and Simulation

[HAg] Hashin, Z., Assessment of the self-consistent scheme approximation: Conductivity of particulate composites, J. Composite
Materials, 2 (1968), 284300.
[HAh] Hashin, Z., Theory of Fiber Reinforced Materials, NASA
CR-1974 (1972).
[HAi] Hashin, Z., Analysis of composite materialsa survey, J. Appl.
Mech., 50 (1983), 481505.
[HAj] Hashin, Z., The differential scheme and its application to
cracked materials, J. Mech. Phys. Solids, 36 (1988), 719734.
[HAk] Hashin, Z., and Shtrikman, S., On some variational principles
in anisotropic and nonhomogeneous elasticity, J. Mech. Phys.
Solids, 10 (1962), 335342.
[HAl] Hashin, Z., and Shtrikman, S., A variational approach to the
theory of the effective magnetic permeability of multiphase
materials, J. Appl. Phys., 33 (1962), 31253131.
[HAm] Hashin, Z., and Shtrikman, S., A variational approach to the
theory of the elastic behavior of multiphase materials, J. Mech.
Phys. Solids, 11 (1963), 127140.
[HAn] Hashin, Z., and Rosen, B., The elastic moduli of fiber-reinforced materials, Trans. ASME, J. Appl. Mech., 31 (1964),
223235.
[HAo] Hasselman, D. P. H., and Johnson, L. F., Effective thermal
conductivity of composites with interfacial thermal barrier, J.
Composite Materials, 21 (1987), 505512.
[HAp] Hasselman, D. P. H., Donaldson, K. Y., and Thomas, J. R.
Jr., Effective thermal conductivity of uniaxial composite with
cylindrically orthtropic carbon fibers and interfacial thermal
barrier, J. Composite Materials, 27 (1993) 637644.
[HEa] Heney, F. S., and Pomphrey, N., Self-consistent elastic moduli
of a cracked solid, Geophys. Res. Lett., 9 (1982), 903906.
[HEb] Hershey, A. V., The elasticity of an isotropic aggregate of
anisotropic cubic crystals, J. Appl. Mech., 21 (1954), 236
240.
[HIa] Hill, R., The elastic behavior of a crystalline aggregate, Proc.
Phys. Soc. A, 65 (1952), 349354.

1.7 Micromechanics of Heterogeneous Media

153

[HIb] Hill, R., Elastic properties of reinforced solids: Some theoretical principles, J. Mech. Phys. Solids, 11 (1963), 357362.
[HIc] Hill, R., Theory of mechanical properties of fibre-strengthened
materials. I. Elastic behavior, J. Mech. Phys. Solids, 12
(1964), 199212.
[HId] Hill, R., Theory of mechanical properties of fibre-strengthened
materials. III. Self-consistent model, J. Mech. Phys. Solids, 13
(1965), 189198.
[HIe] Hill, R., A self-consistent mechanics of composite materials,
J. Mech. Phys. Solids, 13 (1965), 213222.
[HOa] Hori, M., Statistical theory of effective electrical, thermal,
and magnetic properties of random heterogeneous materials.
I. Perturbation expansions for the effective permittivity of cell
materials, J. Math. Phys., 14 (1973), 514523.
[HOb] Hori, M., and Yonezawa, F., Statistical theory of effective electrical, thermal, and magnetic properties of random heterogeneous materials. III. Perturbation treatment of the effective
permittivity in completely random heterogeneous materials,
J. Math. Phys., 15 (1974), 21772185.
[HOc] Hori, M., and Yonezawa, F., Statistical theory of effective electrical, thermal, and magnetic properties of random heterogeneous materials. IV. Effective-medium theory and cumulant
expansion method, J. Math. Phys., 16 (1975), 352364.
[HUa] Hughes, B. D., and Prager, S., Random processes and random systems. An introduction, in The Mathematics and
Physics of Disordered Media, Hughes, B. D., and Ninham, B. W., eds., Lecture Notes in Mathematics, vol. 1035
(1983), 1108.
[JEa] Jeffrey, D. J., Conduction through a random suspension of
spheres, Proc. Roy. Soc. London, A335 (1973) 355367.
[JEb] Jennings, B. R., and Parslow, K., Partcicle size measurement:
The equivalent spherical diameter, Proc. R. Soc. Lond. A, 419
(1988), 137149.
[JEc] Jeulin, D., Bounds of physical properties of some random
structures models, in Continuum Models and Discrete
Systems, Inan, E., and Markov, K. Z., eds., World Sci. (1998),
147154.

154

Heterogeneous Media: Modelling and Simulation

[KAa] Kachanov, M., Effective elastic properties of cracked solids:

critical review of some basic concepts, Appl. Mech. Review,
45 (1992), 304335.
[KAc] Kanaun, S. K., On singular models of a thin inclusion in a homogeneous elastic medium, Appl. Math. Mechanics (PMM),
48(1) (1984), 8191. (In Russian.)
[KAd] Kanaun, S. K., and Kudrjavtseva, L. T., Temperature stresses
in the composites with spherically layered inclusions, MTT
(Mechanics of Solids), No 4 (1987), 113121. (in Russian.)
[KAe] Kanaun, S. K., Kudrjavtseva, L. T., and Goldman, A. Y.,
Anomalous behavior of the viscoelastic properties of some
polymer composites, Mechanics of Composite Materials, No
3 (1988), 442448. (in Russian.)
[KAf] Kanaun, S. K., and Levin, V. M., Effective Field Method
in the Mechanics of Composite Materials, Petrozavodsk
University (1993). (in Russian.)
[KAg] Kanaun, S. K., and Levin, V. M., Effective field method in mechanics of matrix composite materials, in Recent Advances
in Mathematical Modelling of Composite Materials,
Markov, K. Z., ed., World Sci. (1994), 158.
[KAh] Kapitza, P. L., Collected Papers, Vol. 2, ter Haar, D., ed.,
Pergamon Press, Oxford, 2nd edition (1965).
[KEa] Kerner, E. H., The electrical conductivity of composite media,
Proc. Phys. Soc. B, 59 (1956), 802807.
[KIa] Kirkpatrick, S., Percolation and conductivity, Rev. Mod. Phys.,
45 (1973), 574588.
[KOa] Kohler, W. E., and Papanicolaou, G. C., Some applications
pf the coherent potential approximation, in Multiple Scattering and Waves in Random Media, Chow, P. L. et al.,
eds., North-Holland, New York and Oxford (1981), 199223.
[KRa] Kraichnan, R. H., The closure problem in turbulence theory,
in Proc. Symp. Appl. Math., AMS, Providence, Rhode
Island, 13 (1962), 199225.
[KRb] Kreher, W., Internal stresses and relations between effective
thermoelastic properties of stochastic solidssome exact solutions, Z. Angew. Math. Mech., 68(3) (1988), 147154.

1.7 Micromechanics of Heterogeneous Media

155

[KRc] Krivoglaz, M., and Cherevko, A., On the elastic moduli of a

two-phase solid. The Physics of Metals and Metallography,
8(2) (1959), 14. (Engl. Transl. from Phys. Metallov i Metallov., Moscow, 8(2) (1859) 161168.)
[KRd] Kroner, E., Berechnung der Elastischen Konstanten des Vielkristalls aus den Konstanten des Einkristalls, Z. Physik, 151
(1958), 504518,
[KRe] Kroner, E., Bounds for the self-consistent tensor of elastic
moduli, Nucl. Sci. Appl., B (1977), 5255,
[KRf] Kroner, E., Statistical modelling, in Modelling Small Deformation of Polycrystals, Gittus, J., and Zarka, J., eds.,
Elsevier, Barking, Essex, England (1986), 229291.
[KUa] Kunin, I. A., Elastic Medium with Microstructure, Vol. 2,
Springer-Verlag, New York-London-Heidelberg-Tokio (1983).
[KUb] Kunin, I. A., and Sosnina, E. G., An ellipsoidal inhomogeneity
in an elastic medium, J. Appl. Math. Mech. (PMM), 37 (1971)
501504. (Translated from Russian.)
[KUc] Kuster, G. C., and Toksoz, M. N., Velocity and attenuation of
seismic waves in two-phase media. I. Theoretical formulations.
II. Experimental results, Geophysics, 39 (1974), 587618.
[LAa] Landau, L. D., and Lifshitz, E. M., Electrodynamics of
Continuous Media, Pergamon Press, Oxford (1960).
[LAb] Landauer, R., The electrical resistance of binary metallic mixtures, J. Appl. Phys., 23 (1952), 779784.
[LAc] Landauer, R., Electrical conductivity in inhomogeneous media, in Electrical Transport and Optical Properties of
Inhomogeneous Media, Garland, J. C., and Tanner, D. B.,
eds., AIP Conf. Proc. No. 40, Am. Inst. Phys., New York
(1978), 243.
[LAd] Lax, M., Multiple scattering of waves. II. The effective field
in dense systems, Rev. Mod. Phys., 85 (1952), 621629.
[LEa] Levin, V. M., On the coefficients of thermal expansion of heterogeneous materials, Mekhanika Tverdogo Tela, No 1 (1967),
8893. (In Russian.) (English Transl. Mechanics of Solids, 2
(1967), 5861.)

156

Heterogeneous Media: Modelling and Simulation

[LEb] Levin, V. M., On determination of the elastic and thermoelastic moduli of composite materials, Mekhanika Tverdogo Tela,
No 6 (1976), 136145. (In Russian.)
[LIa] Lipton, R., Influence of interfacial conduction on the DC electrical conductivity of particle reinforced composites, Proc. Roy.
Soc. London Ser. A, 454 (1998), 13711382.
[LOa] Lomakin, V. A., On deformation of microinhomogeneous solids, Appl. Math. Mech. (PMM), 29 (1965), 888893. (In Russian.)
[LOb] Lorentz, H. A., The Theory of Electrons, B. U. Teubner,
Leipzig (1909); reprint Dover, New-York (1952).

[LOc] Lorenz, L., Uber

die Refraktionskonstante, Ann. Phys. Chem.,
11 (1880), 70ff.
[LOd] Love, A. E. H., A Treatise on the Mathematical Theory
of Elasticity, Cambridge University Press (1927).
[MAa] Mackenzie, J. K., Elastic constants of a solid containing spherical holes, Proc. Phys. Soc. London B, 63 (1950), 211.
[MAb] Mal, A., and Knopoff, L., Elastic waves velocities in twocomponents systems, J. Ins. Math. Appl., 3 (1967), 376387.
[MAc] Malone, G. H., Suh, S. Y., and Prager, S., Variational bounds
on the diffusive and hydrodynamics permeabilities of randomly perforated sheets, in The Mathematics and Physics
of Disordered Media, Hughes, B. D., and Ninham, B. W.,
eds., Lecture Notes in Mathematics, vol. 1035 (1983), 370390.
[MAd] Markov, K. Z., Mechano-mathematical modelling of microinhomogeneous deformable solids, Doctor of Sci. Dissertation,
Bulg. Acad. of Sciences and Sofia University (1981), 1243.
(In Bulgarian.)
[MAe] Markov, K. Z., One-particle approximations in mechanics
of composite materials, in Continuum Models and Discrete Systems, Brulin, O., and Hsieh, R. K. T., eds., NorthHolland (1981), 441448.
[MAf] Markov, K. Z., On the cluster bounds on the effective properties of microcracked solids, J. Mech. Phys. Solids, 46 (1998),
357388.
[MAg] Markov, K. Z., On a two-point correlation function in random
dispersions and an application, in Continuum Models and

1.7 Micromechanics of Heterogeneous Media

157

Discrete Systems, Inan, E., and Markov, K. Z., eds., World

Sci. (1998), 206215.
[MAi] Matheron, G., Quelques inegalites pour la permeabilite effective dun milieu poreux heterog`ene, Cahiers de Geostatistique,
Fasc. 3 (1993), 120.
[MAj] Maxwell, J. C., A Treatise on Electricity and Magnetism,
Dover, New York (1954). (Republication of 3rd edition of
1891.)
[MCa] McCoy, J. J., On the displacement field in an elastic medium
with random variations in material properties, in Recent Advances in Engineering Sciences, Vol. 5, Eringen, A. C.,
ed., Gordon and Breach, New York (1970), 235254.
[MCb] McLaughlin, R., A study of the differential scheme for composite materials, Int. J. Eng. Sci., 15 (1977), 237244.
[MIa] Miloh, T., and Benveniste, Y., On the effective conductivity
of composites with ellipsoidal inhomogeneities and highly conducting interfaces, Proc. Roy. Soc. London Ser. A, 455 (1999),
in press.
[MIb] Miller, M. N., Bounds for effective electrical, thermal and
magnetic properties of heterogeheous materials, J. Math.
Phys., 10 (1969), 19882004.
[MIc] Milton, G. W., Bounds on the electromagnetic, elastic and
other properties of two-component composites, Phys. Review
Lett., 46 (1981), 542545.
[MId] Milton, G. W., Bounds on the complex permittivity of a twocomponent composites, J. Appl. Phys., 52 (1981), 52365241.
[MIe] Milton, G. W., The coherent potential approximation is realizable effective medium scheme, Commun. Math. Phys., 99
(1985), 463500.
[MIf] Milton, G. W., Composites: A myriad of microstructure independent relations, in Theoretical and Applied Mechanics
1996, Tatsumi, T. et al., eds., Elsevier (1997), 443459.
[MIg] Milton, G. W., On characterizing the set of possible effective
tensors of composites: the variational method and the translation method, Comm. Pure Appl. Math., 43 (1990), 63125.

158

Heterogeneous Media: Modelling and Simulation

[MIh] Milton, G. W., and Golden, K., Thermal conduction in composites, in Thermal Conductivity 18, Ashworth, T., and
Smith, D. R., eds., Plenum Publ. Co. (1985), 571582.
[MOa] Molyneux, J., and Beran, M., Statistical properties of the
stress and strain fields in a medium with small variations in
elastic constants, J. Math. Mech., 14 (1965), 227351.
[MOb] Mori, T., and Tanaka, K., Average stress in matrix and average elastic energy of materials with misfitting inclusions, Acta
Metall., 21 (1973), 571574.
[MOc] Mossotti, O. F., Sur les forces qui regissent la constitution
interieur des corps apercu pour servir `a la determination de
la cause et des lois de laction moleculaire, Ac. Sci. Torino, 22
(1836), 136.
[MOd] Mossotti, O. F., Discussione analitica sulinfluenza che lazione
di un mezzo dielettrico ha sulla distribuzione dellelettricit`a
alla superficie di pi`
u corpi elettrici disseminati in eso, Mem.
Mat. Fis. della Soc. Ital. di Sci. in Modena, 24 (1850), 4974.
[MUa] Mura, T., Micromechanics of Defects in Solids, 2nd edition, Martinus Nijhoff, Dordrecht (1988).
[MUb] Muratov, R. Z., The Potentials of an Ellipsoid, Atomizdat, Moscow (1975). (In Russian.)
[NAa] Navier, Claude-L.-M.-H., Memoire sur les lois de lequilibre et
du movement des corps solides elastiques, Memoire de lAcademie royale des Sciences, VII (1827), 375394.
[NEa] Nemat-Nasser, S., and Hori, M., Micromechanics: Overall
Properties of Heterogeneous Solids, Elsevier (1993).
[NOa] Norris, A. N., A differential scheme for the effective moduli of
composites. Mech. Mat., 4 (1985), 116.
[NOb] Norris, A. N., Callegari, A. J., and Sheng, P., A generalized
differential effective medium theory, J. Mech. Phys. Solids, 33
(1985), 525543.
[NOc] Norris, A. N., Sheng, P., and Callegari, A. J., Effective medium
theories for two-phase dielectric media, J. Appl. Phys., 57
(1985), 19901996.
[PAa] Papanicolaou, G., Effective parameters and fluctuations for
boundary value problems, in Transactions of the 25th Con-

1.7 Micromechanics of Heterogeneous Media

159

ference of Army Mathematicians (1979), ARO Report

80-1 (1980), 733744.
[PAb] Parmenter, R. H., Energy levels of a disordered alloy, Phys.
Review, 97 (1955), 587598.
[PAc] Paul, B., Prediction of the elastic constants of multiphase materials, Trans AIME, 218 (1960), 3641.
[PHa] Pham Huy, H., and Sanchez-Palencia, E., Phenomenes de
transmission a travers des couches minces de conductivite
elevee, J. Math. Anal. and Appl., 47 (1974), 284309.
[POa] Poisson, S. D., Memoire sur la theorie du magnetisme, Mem.
de lAcad. roy. de France, V (1824), 247338. (copy from
Archives de lAcademie des Sciences, Paris.)
[POb] Poisson, S. D., Second Memoire sur la theorie du magnetisme,
Mem. de lAcad. roy. de France, V (1824), 488533. (copy from
Archives de lAcademie des Sciences, Paris.)
[POc] Poisson, S. D., Memoire sur la lequilibre et le mouvement des
corpes solides elastiques et des fluides (Extrait), Annales de
chimie et de physique, XLII (1829), 145171.
[RAa] Lord Rayleigh, On the influence of obstacles arranged in rectangular order upon the volume properties of a medium, Phil.
Mag., 34 (1892), 481502.
[REa] Reuss, A., Berechnung der Fliessgrenze von Mischkristallen
auf Grund der Plastizitatsbedingung f
ur Einkristalle,
Z. angew. Math. Mech., 9 (1929), 4958.
[REb] Reynolds, J. A., and Hough, J. M., Formulae for dielectric
constant of mixtures, Proc. Phys. Soc., 70 (1957), 369775.
[ROa] Roscoe, R., The viscosity of suspensions of rigid spheres, Brit.
J. Appl. Physics, 3 (1952), 267269.
[ROb] Rosen, B. W., Effective thermal expansion coefficient of composite materials, Ph. D. Dissertation, University of Pennsylvania, 1968.
[ROc] Rosen, B. W., and Hashin, Z., Effective thermal expansion
coefficient and specific heats of composite materials, Int. J.
Eng. Sci., 8 (1970), 157173.
[SAa] Sanchez-Palencia, E., Non-Homogeneous Media and Vibration Theory, Springer-Verlag, New York (1980).

160

Heterogeneous Media: Modelling and Simulation

[SEa] Sen, P. N., Kenyon, W., Takezaki, H., and Petricola, M. J.,
Formation factor of carbonate rocks with microporosity:
model calculations, G. Petroleum Sci. Eng., 17 (1997), 345
352.
[SHa] Shapery, R. A., Thermal expansion coefficients of composite materials, based on energy principles, J. Compos. Mat., 2
(1968), 380404.
[SHb] Shermergor, T. D., Theory of Elasticity of Microinhomogeneous Solids, Nauka, Moscow (1977). (In Russian.)
[SKa] Skorohod, V. V., Calculation of the effective isotropic moduli
of disperse solid systems, Poroshkovaja Metallurgija (Powder
metall.), No 1 (1961), 5055. (In Russian.)
[SMa] Smoluchowski, M. v., Drei Vortrage u
ber Diffusion, Brownsche Molekularbewegung und Koagulation von Kolloidteichen,
Physik. Zeitschr., 17 (1916), 557571, Schlu, 585599.
[SOa] Sotkilava, O. V., and Cherepanov, G. P., Some problems of
the nonhomogeneous theory of elasticity, Appl. Mat. Mech.
(PMM), 38(3) (1974), 539550. (In Russian.)
[SPa] Spencer, A. J. M., Theory of invariants, in Continuum Physics, vol. 1, Eringen, A. C., ed., Academic Press, New York
(1971), 239353.
[STa] Stratton, J. A., Electromagnetic Theory, Mc Graw-Hill,
New York (1941).
[TAa] Talbot, D. R. S., and Willis, J. R., The effective sink strength
of a random array of voids in irradiated material, Proc. R.
Soc. London A, 370 (1980), 351374.
[TAb] Talbot, D. R. S., and Willis, J. R., The overall sink strength of
an inhomogeneous lossy medium. I. Self-consistent estimates,
Mech. Materials, 3 (1984), 171181.
[TAc] Talbot, D. R. S., and Willis, J. R., Bounds and self-consistents
estimates for overall properties of nonlinear composites, IMA
J. Appl. Math., 39 (1987), 215240.
[TAd] Taylor, G. I., The viscosity of a fluid containing small drops
of another fluid, Proc. R. Soc. Lond., 138 (1932), 4148.
[THa] Theocharis, P. S., The Mesophase Concept in Composites, Springer-Verlag, Berlin (1987).

1.7 Micromechanics of Heterogeneous Media

161

[TOa] Toland, J. F., and Willis, J. R., Duality for families of natural variational principles in nonlinear electrostatics, SIAM J.
Math. Anal., 20 (1989), 12831292.
[TOb] Torquato, S., Microscopic approach to transport in two-phase
random media, Thesis, State University of New York at Stony
Brook, 1980.
[TOc] Torquato, S., Random heterogeneous media: Microstructure
and improved bounds on effective properties, Appl. Mech.
Rev., 44 (1991), 3776.
[TOd] Torquato, S., and Rintoul, M. D., Effect of the interface on
the properties of composite media, Physical Review Letters,
75 (1995), 40674070.
[TOe] Torquato, S., and Stell, G., Microstructure of two-phase random media. I. The n-point probability functions, J. Chem.
Phys., 77 (1982), 20712077.
[TRa] Trimarco, C., The Toupin-Mindlin theory of dielectrics in the
light of the Mossottis idea, Bulletin Polish Acad. Sci., Techn.
Sci., 42 (1994), 495504.
[VAa] Van Kampen, N. G., and Felderhof, B. U., Theoretical
Methods in Plasma Physics, North-Holland, Amsterdam
(1967).

[VOa] Voigt, W., Uber

die Beziehung zwischen den beiden Elastizitatskonstanten isotroper Korper, Wied. Ann., 38 (1889),
573587.
[WAa] Walsh, J. B., The effect of cracks on the compressibilty of
rocks, J. Geophys. Res., 70 (1965), 381389.
[WAb] Watt, J. P., The elastic properties of composite materials,
Review Geophys. Space Phys., 14 (1976), 541563.
[WEa] Weng, G. J., Some elastic properties of reinforced solids with
special reference to isotropic ones containing spherical inclusions, Int. J. Eng. Sci., 22 (1984), 845856.
[WEb] Weng, G. J., The theoretical connection between MoriTanakas theory and the Hashin-Shtrikman-Walpole bounds,
Int. J. Eng. Sci., 28 (1990), 11111120.
[WIa] Wiener, O., Die Theorie des Mischkorpers f
ur das Feld des stationaaren Stromung. Erste Abhandlung die Mttelswertsatze

162

Heterogeneous Media: Modelling and Simulation

f
ur Kraft, Polarisation und Energie. Abh. Math.-Physichen
Klasse Konigl. Sacsh. Gessel. Wissen., 32(6) (1912), 509604.

[WIb] Willis, J. R., Bounds and self-consistent estimates for the overall properties of anisotropic composites, J. Mech. Phys. Solids,
25 (1977), 185202.
[WIc] Willis, J. R., Variational and related methods for the overall
properties of composites, in Adv. Appl. Mechanics, vol. 21,
Yih, C., ed., Academic Press, New York (1981), 178.
[WId] Willis, J. R., Elasticity theory of composites, in Mechanics
of Solids, Hopkins, H. G., and Sewell, M. J., eds., Pergamon
Press, Oxford and New York (1982), 653686.
[WIe] Willis, J. R., Randomly inhomogeneous media, in Homogenization Techniques for Composite Media, Proc. Lect.
CISM, Udine, 1985, Lect. Notes Phys. 272 (1987), 279336.
[WIf] Willis, J. R., The structure of overall constitutive relations
for a class of nonlinear composites, IMA J. Appl. Math., 43
(1989), 237242.
[WUa] Wu, T. T., The effect of inclusion shape on the elastic moduli
of a two-phase material, Int. J. Sol. Structures, 2 (1966), 28.
[WUb] Wu, C.-T. D., and McGullough, R. L., Constitutive relationships for heterogeneous materials, in Developments in
Composite Materials, Holister, G. S., ed., Appl. Science
Publishers, London (1977), 119187.
[ZHa] Zhikov, V., Kozlov, S., and Oleinik, O., Homogenization of
Differential Operators, Springer, Berlin (1985).
[ZIa] Ziman, J. M., Models of Disorder. The Theoretical Physics of Homogeneously Disordered Systems, Cambridge
University Press, Cambridge (1979).
[ZIb] Zimmerman, R. W., Elastic moduli of a solid with spherical
pores: New self-consistent method, Int. J. Rock Mech., 21
(1984) 339343.
[ZIc] Zimmerman, R. W., Elastic moduli of a solid containing spherical inclusions, Mech. Materials., 12 (1991), 1724.
[ZId] Zimmerman, R. W., Behavior of the Poisson ratio of a twophase composite material in the high-concentration limit,
Appl. Mech. Rev., 47 (1994), S38S44.