See discussions, stats, and author profiles for this publication at: https://www.researchgate.

net/publication/303289895

Explicit effective elasticity tensors of two-phase periodic composites with

spherical or ellipsoidal inclusions

Article in International Journal of Solids and Structures · April 2016

DOI: 10.1016/j.ijsolstr.2016.05.005

CITATIONS READS

8 189

3 authors:

Quy Dong To Guy Bonnet

University of Paris-Est Université Gustave Eiffel
75 PUBLICATIONS 423 CITATIONS 179 PUBLICATIONS 2,146 CITATIONS

SEE PROFILE SEE PROFILE

H. Hoang-Duc
University of Paris-Est
4 PUBLICATIONS 27 CITATIONS

SEE PROFILE

Some of the authors of this publication are also working on these related projects:

Homogenization of transport properties of porous media View project

All content following this page was uploaded by Guy Bonnet on 14 March 2018.

The user has requested enhancement of the downloaded file.

Explicit effective elasticity tensors of two-phase periodic
composites with spherical or ellipsoidal inclusions
Quy-Dong Toa,b,∗, Guy Bonneta , Duc-Hieu Hoanga
a
Université Paris-Est, Laboratoire Modelisation et Simulation Multi Echelle, UMR 8208
CNRS, 5 Boulevard Descartes, 77454 Marne-la-Vallée Cedex 2, France
b
Duy Tan University - Institute of Research & Development, K7/25 Quang Trung,
Danang, Vietnam

Abstract

The effective elasticity tensors of two-phase composites are estimated by

solving the localization problem in the wave-vector domain for the case of non
overlapping spherical or ellipsoidal inclusions. With previous works showing
that the effective properties can be computed from lattice sums, we propose a
method to compute the sums analytically and obtain the explicit expressions
for the effective tensors. In the case of different periodic cells leading to cubic
or orthotropic elasticity tensors, the effective elasticity tensors are obtained
in closed forms that are in good agreement with the exact solutions for a
large range of physical parameters. In the random distribution cases, the
statistical connection of the effective tensor to the structure factor is shown
and a closed-form expression is obtained in the infinite volume limit.
Keywords: Closed-form expression; Fourier transform; effective elasticity
tensors; non overlapping spherical inclusions; ellipsoidal inclusions; periodic
problem; random distribution; orthorhombic lattice; structure factor; form
factor

∗
Email: quy-dong.to@u-pem.fr

Preprint submitted to Elsevier May 17, 2016

1. Introduction

One objective of micromechanics is to model the overall behavior of com-

posites by studying physical problems at the scale of heterogeneities, fol-
lowed by averaging the physical quantities on volumes of interest. The usual
homogenization procedure is to construct a representative volume element
(RVE) containing the distribution of different constituting phases and then
to use continuum mechanics to analyze the localization problem. For elastic
periodic media, it was shown that the effective elasticity tensor is rigorously
defined from the average stress/strain linear relation (Sanchez-Palencia, 1974,
1980) obtained by studying only a unit cell.

Most closed-form estimations of the effective elasticity tensor are based on

simplifications of RVE containing assemblages of coated spheres or ellip-
soids (Eshelby, 1957; Christensen and Lo, 1979; Benveniste and Milton, 2003;
Tsukrov and Kachanov, 2000). For periodic problems, the solutions are usu-
ally obtained from numerical methods such as the finite element method
(FEM), boundary element method (BEM) and fast Fourier transform (FFT)
(Kaminski, 1999, 2005; Liu et al., 2005; Michel et al., 1999; Eyre and Milton,
1999; Monchiet and Bonnet, 2012). Generally, it is more difficult to apply
analytical techniques to periodic microstructures due to the special bound-
ary conditions and the geometry of the unit cell. The most sophisticated
complete semi-analytical methods use an expansion of the solution on the
basis of periodic functions and expand the basis functions on spherical har-
monics or spheroidal functions (Nunan and Keller, 1984; Sangani and Lu,
1987; Kushch, 1997, 2013), whereas estimates can be obtained by perturba-
tion methods (Cohen, 2004).

In this paper, we use the Nemat-Nasser-Iwakuma-Hejazi (NIH) approxima-

tion (Nemat-Nasser et al., 1982) to treat the specific case of two-phase pe-
riodic composites with spherical or ellipsoidal inclusions. This method has

2
the same level of approximation as provided by the Clausius-Mossotti (CM)
assumption, that is, it provides good estimates for a large range of inclusion
concentrations. It takes into account the distribution of inclusions, but it is
not more accurate for very large concentrations of inclusions. In this con-
text, it rests on Fourier expansion and is very versatile, thus accounting for
the interaction between inclusions or cracks of different kinds: such as ellip-
soids, and cuboids. Further, it can address a fully anisotropic behavior. This
method is much easier to use than fully numerical methods. Nevertheless,
it necessitates the computation of lattice sums, which requires obviously a
minimum of numerical calculations: storing a large number of Fourier com-
ponents, checking the convergence of the lattice sums, etc. In some cases, like
for some cubic cases, alternative methods provided, under similar assump-
tions, simple expressions of the effective elasticity tensors (Cohen, 2004).
For further applications, the computation of lattice sums for further appli-
cations has been to interpolate can also be avoided by interpolating some
results coming from the NIH method, like for fiber composites (Luciano and
Barbero, 1994). The main contribution of our work is also to avoid the com-
putation of the lattice sums, but by estimating these sums analytically in
more general periodic cells containing spheres or ellipsoids. By considering
such different arrangements inside the unit cell (e.g., simple orthorhombic
(SO), body-centered orthorhombic (BCO), and face-centered orthorhombic
(FCO)), we derive explicit closed-form expressions of the effective elasticity
tensors.

The next part of this work deals with the random distribution case. As an
extension of our previous work on conduction phenomena (To et al., 2013; To
and Bonnet, 2015), we establish the statistical connections between the effec-
tive tensor and local arrangement factors such as the form factor P(ξ), the
structure factor S(ξ), and the scattering intensity I(ξ) used in the scatter-
ing theory (Hunter and White, 2001; Bohren and Huffman, 1998) and more

3
generally in solid-state physics (Rössler, 2009). These factors provide useful
insight into the local structure of the particles, and it is worth noting that
they can be obtained experimentally. Finally, a closed-form solution is also
obtained in this case.

The present paper is composed of five sections. After the Introduction, Sec-
tion 2 is dedicated to the homogenization theory of periodic media using the
integral equation approach and NIH estimation. In Section 3, closed-form
solutions for elasticity tensors related to different arrangements of spherical
or ellipsoidal inclusions are derived. Section 4 deals with the random cases;
some numerical applications and comparisons with different solutions of the
literature are reported in Section 5 and finally, the summary and remarks
are presented in Section 6.

The usual notations of tensor algebra are adopted throughout the paper.
For example, tensors are in bold characters, tensor products are denoted by
” ⊗ ” (tensor products) and ”., : ” (inner products). The Einstein summation
convention is used for repeated indices.

2. Homogenization of elastic periodic composites

First, the governing integral equations of the homogenization problem

are derived along classical means (Christensen, 1979; Milton, 2002), how-
ever with more attention given to the periodic case and the use of Fourier
transforms. Next, the estimation of effective elasticity tensors based on the
integral equations (Nemat-Nasser et al., 1982) is introduced.

2.1. Governing integral equations

We consider an infinite elastic composite the fourth-order elasticity tensor
C of which is a periodic function of the local coordinates x(x1 , x2 , x3 ) with

4
periods a1 , a2 , a3 :

C(x1 , x2 , x3 ) = C(x1 + n1 a1 , x2 + n2 a2 , x3 + n3 a3 ), ∀n1 , n2 , n3 ∈ Z (1)

The homogenization procedure of the periodic material was rigourously es-

tablished from the asymptotic development of the involved quantities, stress
σ, strain ϵ, and displacement u in terms of the scaling parameter ε (Sanchez-
Palencia, 1974, 1980). Matching the powers of ε in the elasticity equations
yields different relations between the quantities, including the definition
of the effective elasticity tensor Ce . To summarize, the following periodic
boundary value problem in a unit cell V :

σ(x) = C(x) : ϵ(x) ∀x ∈ V

1
ϵ(x) = (∇u(x) + ∇T u(x)), ∀x ∈ V
2
∇ · σ(x) = 0, ∀x ∈ V
u(x) − E.x periodic,
σ(x).n antiperiodic (2)

must be solved, allowing finally the computation of the effective elasticity

tensor Ce from the relation between macroscopic strain E and stress Σ:

Σ = Ce : E, Σ = ⟨σ⟩V , E = ⟨ϵ⟩V . (3)

Here, we adopt the notation ⟨⟩V to refer to the average over volume V of the
quantity inside the brackets, for example,
∫
1
⟨ϕ⟩V = ϕdV. (4)
V

Due to the periodicity of the problem, it is useful to express the periodic

quantities in the form of Fourier series and apply Fourier analysis to the
elasticity equations (2). For example, if ϕ is a periodic function of x1 , x2 , x3 ,

5
it can be expressed as a Fourier series:
∑ ⟨ ⟩
ϕ(x) = ϕ(ξ)eiξ.x , ϕ(ξ) = ϕ(x)e−iξ.x V , (5)
ξ̸=0

where ϕ can be stress σ, strain ϵ, elasticity tensor C, etc. For the sake of
simplicity, we differentiate quantities in Fourier space and in real space by
adding the variable after the same symbol: ϕ(ξ) is the Fourier transform of
the real function ϕ(x) defined by (5)2 . We note that the infinite sum in (5)1
involves all discrete wave vectors ξ with components ξ1 , ξ2 and ξ3 satisfying

2πni
ξi = , i = 1, 2, 3, n1 , n2 , n3 ∈ Z. (6)
ai

The periodic boundary value problem (2) can be solved by an integral equa-
tion approach. By introducing a constant reference elasticity tensor C0 and
polarization tensor σ ∗ (or eigenstress tensor) defined by

σ(x) = C0 : ϵ(x) + σ ∗ (x), (7)

it can be shown that ϵ, σ can be determined via the expressions

ϵ(ξ) = −Γ0 (ξ) : σ ∗ (ξ). (8)

In (8), Γ0 (ξ) is the Green operator for strain in Fourier space. When the
reference material is isotropic with Lamé constants λ0 and µ0 , for example,

0
Cijkl = λ0 δij δkl + µ0 (δik δjl + δil δjk ), (9)

6
tensor Γ0 admits the following form:

1
Γ0ijkl (ξ) = (δik ξ¯j ξ¯l + δil ξ¯j ξ¯k + δjk ξ¯i ξ¯l + δjl ξ¯i ξ¯k )
4µ0
λ0 + µ0 ¯ ¯ ¯ ¯
− ξi ξj ξk ξl . (10)
µ0 (λ0 + 2µ0 )

In (9,10), δij is the usual delta Kronecker symbol and ξ¯i is the direction cosine
of the wave vector ξ. By using the definition of the polarization tensor (7)
and relation (8), we obtain the integral equation in ϵ(x):

ϵ(x) = E − Γ0 ∗ (C(x) − C0 ) : ϵ(x). (11)

The transformation of (11) into an equation for eigenstrain ϵ∗ is straightfor-

ward by making use of the relations defining the eigenstrain ϵ∗ :

σ = C0 : [ϵ(x) − ϵ∗ (x)],
or equivalently (C0 − C(x)) : ϵ(x) = C0 : ϵ∗ (x) = −σ ∗ . (12)

Finally, we obtain the same equation as that reported by Nemat-Nasser et al.

(1982):
[ ]
∑
C0 : ϵ∗ (x) = (C0 − C(x)) : E + eiξ.x Γ0 (ξ) : C0 : ϵ∗ (ξ) . (13)
ξ̸=0

2.2. Estimations of effective elasticity tensors based on integral equations

From now on, we consider the specific case of a matrix-inclusion composite
where each phase is isotropic with elasticity tensors and Lamé constants being
Cm , λm , µm (matrix), and Ci , λi , µi (inclusion). Following the NIH procedure
(Nemat-Nasser et al., 1982), we shall estimate the effective elastic properties
on the basis of integral equation (13). Taking the matrix as the reference
material, that is, C0 = Cm and averaging both sides of (13) over the inclusion

7
volume Ω with fraction f = Ω/V yields
[ ]
∑
Cm : ⟨ϵ∗ (x)⟩Ω = (Cm − Ci ) : E + ⟨eiξ.x ⟩Ω Γm (ξ) : Cm : ϵ∗ (ξ) . (14)
ξ̸=0

Nemat-Nasser et al. (1982) noted that the eigenstrain vanishes outside the
inclusion and thus proposed the following approximate evaluation of ϵ∗ (ξ):

ϵ∗ (ξ) = f ⟨ϵ∗ (x)e−iξ.x ⟩Ω ≃ f ⟨ϵ∗ (x)⟩Ω ⟨e−iξ.x ⟩Ω . (15)

By defining the following shape functions I(ξ) and P (ξ)

⟨ ⟩ f
I(ξ) = Ω eiξ.x Ω , P (ξ) = I(ξ)I(−ξ), (16)
Ω2

and applying the approximation (15) to (14), we obtain

Cm : ⟨ϵ∗ (x)⟩Ω ≃ (Cm − Ci ) : [E + Υ : Cm : ⟨ϵ∗ (x)⟩Ω ] , (17)

where the tensor Υ is the following lattice sum in the reciprocal space:
∑
Υ= P (ξ)Γm (ξ). (18)
ξ̸=0

Inverting (17) yields the average eigenstrain ⟨ϵ∗ (x)⟩Ω

[ ]−1
⟨ϵ∗ (x)⟩Ω = (Cm − Ci )−1 : Cm − Υ : Cm : E. (19)

Next, we average (12) with the matrix as the reference material, and we find
the relation between E, Σ, and ⟨ϵ∗ (x)⟩Ω :

Σ = Cm : (E − ⟨ϵ∗ (x)⟩V ), or Σ = Cm : (E − f ⟨ϵ∗ (x)⟩Ω ). (20)

8
Finally, comparing (20) with (3) and accounting for (19), we can derive the
overall tensor Ce :
[ ]−1
Ce = Cm − f (Cm − Ci )−1 − Υ . (21)

The success of the NIH estimation relies on the accuracy of Eq. (15). The-
oretically speaking, if the inclusions have an ellipsoidal shape and occupy
a sufficiently small volume fraction f , approximation (15) is valid. Indeed,
Eshelby (1957) proved that the stress/strain fields inside an ellipsoidal inclu-
sion embedded in an infinite matrix and subject to homogeneous stress/strain
boundary conditions at infinity are also homogeneous. For interacting inclu-
sions, the inclusion stress/strain fields are no longer uniform, although (15)
is still expected to yield accurate results for a large range of volume fractions.

2.3. Tensor Υ and its relation with the periodic Eshelby tensor
The periodic Eshelby problem can also be addressed, where the inclusions
in the above study are made of the same material as the matrix (i.e., Ci =
Cm ) but are subjected to a uniform eigenstrain ϵ∗ = E∗ inside their domains
Ω. Different from the classical Eshelby problems, the strain field inside each
inclusion is generally not uniform. However, if we average this strain field
over each inclusion, we can still define the periodic Eshelby tensor Sp as
follows:
⟨ϵ(x)⟩Ω = Sp : E∗ , (22)

As strain field ϵ can be computed from ϵ∗ via (8,12) with relation (15) be-
coming exact, this tensor can be determined by
∑
Sp = P (ξ)Γm (ξ) : Cm . (23)
ξ̸=0

It is clear that if the domain V is sufficiently large compared with the size
of the inclusion, Sp should be equal to the classical Eshelby tensor S∞ . Nev-

9
ertheless, there is always a simple connection between Sp and Υ:

Sp = Υ : Cm . (24)

This leads to an alternative interpretation of NIH approximation: the NIH

approximation corresponds to an estimation of the mean strain over the in-
clusions produced by the periodic Eshelby tensor. From (24), we find that the
tensor Υ plays the same role as the classical Hill tensor, but here in the pe-
riodic setting. This ”periodic Hill tensor” provides the average deformation
inside the inclusions due to a constant polarization field.1

On returning to our original problem, combining (74) with (10), this ten-
sor Υ can be written as follows:

1 λm + µm
Υ= W− U (25)
2µm µm (λm + 2µm )

where W and U are fourth-order tensors that depend only on the geometry.

In the following, the unit cell is assumed to be symmetrical with respect

to three orthogonal planes. Under this condition, the matrix representations
of W and U are given in Kelvin’s notation by:
 
2S1 0 0 0 0 0
 
 0 2S2 0 0 0 0 
 
 
 0 0 2S3 0 0 0 
[W] =   (26)
 0 0 0 S2 + S3 0 0 
 
 
 0 0 0 0 S3 + S1 0 
0 0 0 0 0 S1 + S2

1
The classical Hill tensor, frequently used in micromechanics, is defined as the product
between the compliance tensor (Cm )−1 and the classical Eshelby tensor S∞ (Hill, 1965)

10
and  
S4 S9 S8 0 0 0
 
 S9 S5 S7 0 0 0 
 
 
 S8 S7 S6 0 0 0 
[U] =   (27)
 0 0 0 2S7 0 0 
 
 
 0 0 0 0 2S8 0 
0 0 0 0 0 2S9
where Si with i = 1, 2, .., 9 are the lattice sums given by
∑ ∑ ∑
Si = P (ξ)ξ¯i2 , Si+3 = P (ξ)ξ¯i4 , Si+6 = P (ξ)ξ¯j2 ξ¯k2 ,
ξ̸=0 ξ̸=0 ξ̸=0
i, j, k = 1, 2, 3, i ̸= j ̸= k ̸= i. (28)

These lattice sums are not independent because of the relation ξ¯22 +ξ¯22 +ξ¯32 = 1,
which leads to three relations:

Si+3 + Sj+6 + Sk+6 = Si , i, j, k = 1, 2, 3, i ̸= j ̸= k ̸= i. (29)

It implies that only six of these lattice sums are independent and that not
all of the elastic constants are independent, as shown, for example, by Cohen
(2004) in the cubic case. It is worth noting that the symmetry of the effective
tensor is the same as that of tensor Υ and that this tensor is orthotropic,
which corresponds to the symmetry assumptions on the periodic cell. To
obtain the overall elasticity tensor, we need to invert (Cm − Ci )−1 − Υ in
(21). These tensors can be inverted through lengthy expressions for a general
distribution with orthotropic symmetry (Luciano and Barbero, 1994). Fur-
ther, in some special cases (for example, cubic symmetries), we can obtain
simple analytical expressions, as recalled thereafter.

11
3. Closed-form expressions for effective tensors in the case of in-
clusions centered on periodic lattices

3.1. Lattice sum computation and integral approximation

Before obtaining the effective properties, the method for analytical esti-
mation of the lattice sums Si with i = 1, 2, .., 9 is introduced in the general
form
∑
η̄i2α η̄j2β H(η), i, j = 1, 2, 3 (30)
η̸=0

with
η
η = Rξ, H(η) = P (ξ), η̄ = , η = |η|. (31)
η
Here, α and β are non negative integers and R is the inclusion characteristic
length used to generate the dimensionless wave vectors η from ξ. It is noted
that the η points form a rectangular lattice in three-dimensional (3D) space
with density:
a1 a2 a3
ρ= . (32)
(2πR)3
If H(η) is a decaying function of η, we can reasonably estimate the sum
∑ ∫
η̄i2α η̄j2β H(η) ≃ρ ξ¯i2α ξ¯j2β H(η)dη, (33)
ξ∈D D

for a sufficiently remote and large region D where H(η) varies slowly. As
a result, the infinite lattice sum can be estimated by a finite sum of several
leading terms and a continuous integral for longer wave vectors, for example,
∑ ∑ ∫ ∞
η̄i2α η̄j2β H(η) ≃ η̄i2α η̄j2β H(η) +ρ H(η)η̄i2α η̄j2β dη. (34)
η̸=0 |η|<ηc ηc

The cutoff radius ηc in (34) determines the number of leading terms from
the series that are retained in the new formula. In many cases where H(η)
depends on its modulus η alone, we can calculate the integral in the right-

12
 Figure 1: Grid η and η ′ in Fourier space

hand side (RHS) of (34):

∫ ∞ ∫ ∞ ∫
ρ H(η)η̄i2α η̄j2β dη =ρ 2
H(η)η dη η̄i2α η̄j2β dS (35)
ηc ηc |η|=1

The second integral of the RHS involving the unit sphere surface |η| = 1
can be evaluated analytically in numerous cases. The first integral can also
be computed if the explicit expression of H(η) is known. For example, the
simple cubic distribution of spherical inclusions studied in the next section
corresponds to the following expression of H(η):

9[sin(η) − η cos(η)]2
H(η) = f φ(η), φ(η) = . (36)
η6

The associated integrals involving H(η) can be computed by

∫ ∞
ρ H(η)η 2 dη = 9ρf J(ηc ) (37)
ηc

where

π cos 2η 1 1 1 sin 2η cos 2η

J(η) = − − Si(2η) + + 3− − (38)
6 6η 3 2η 6η 3η 2 6η 3

13
and Si(η) is the sine integral
∫ η
sin η ′ ′
Si(η) = dη . (39)
0 η′

The above summation method can be extended to the more general cases
where H(η) is the function of the modulus |L−1 .η| with L−1 being any linear
transformation. Indeed, by simply posing

η ′ = L−1 .η, H(η) = H ′ (η ′ ), (40)

we can calculate the sum using grid η ′ as follows:

∑ (Lik ηk′ )2α (Ljl ηl′ )2β H ′ (η ′ )

(41)
η ′ ̸=0
[(L1k ηk′ )2 + (L2k ηk′ )2 + (L3k ηk′ )2 ]α+β

We can keep several leading terms by using the radius ηc′ and estimating the
remaining series with an integral. The integral estimation method introduced
previously can be applied as the grid η ′ is obtained from η by a homogeneous
deformation L−1 ; thus, the grid η ′ is also uniform. The only difference now
is to compute the integral on the unit sphere surface. Taking the example of
an ellipsoidal inclusion studied in the next section, we have
√
′ ′
H(η) = f φ(η ), η = η12 /χ21 + η22 /χ22 + η22 /χ22 (42)

where χi are dilatation coefficients. The associated surface integral becomes

∫
(χi ηi′ )2α (χj ηj′ )2β
dS, i, j not summed (43)
|η ′ |=1 [(χ1 η1′ )2 + (χ2 η2′ )2 + (χ3 η3′ )2 ]α+β

which can yield closed-form solutions in numerous cases.

Finally, we note that if the ratio R/ai → 0 ∀i = 1, 2, 3, the η grid becomes

infinitely dense. As a result, the relation (33) must hold for all domain D,

14
and the infinite lattice sum is equivalent to the integral over the whole space:
∑ ∫
η̄i2α η̄j2β H(η) = η̄i2α η̄j2β H(η)ρdη. (44)
η̸=0

Although ρ → ∞, the product H(η)ρ is generally bounded for the cases

examined in the later section. As a result, the integral in RHS of (44) is
well defined. Further, one can also expect to recover analytical results cor-
responding to the infinite volume domain limit.

3.2. Cubic lattice arrangements

In this subsection, we consider microstructures composed of identical
spheres of radius R arranged in a cubic lattice of period a (see Fig. 2).
We note that the Fourier transform of the indicator function of a spherical
inclusion of radius R located at xc admits the simple closed form
∫
sin η − η cos η iξ.xc 4πR3
eiξ.x dx = 3Vs e , Vs = , η = ξR. (45)
Vs η3 3

Using the elementary result (45), one can determine I(ξ) and P (ξ) for any
distribution of spheres inside the unit cell. For cubic crystal system arrange-
ments (see Fig. 2), To et al. (2013) derived the following expressions for P (ξ)
(or equivalently H(η)):

P (ξ) = H(η) = αf φ(η), (46)

with φ(η) being defined in (36). The coefficient α depends on the geometry,
equal to 1 for simple cubic arrangement. For other arrangements, it can be
computed as follows:

• Body-centered cubic (BCC)

α = 1 if n1 + n2 + n3 is even, otherwise α = 0, (47)

15
 • Face-centered cubic (FCC)

α = 1 if n1 , n2 , n3 are all even or odd, otherwise α = 0. (48)

Due to the symmetry of the cell and the definition of P (ξ) in (16), the lattice
sums present the equalities

S1 = S2 = S3 , S4 = S5 = S6 , S 7 = S8 = S9 ,

whereas S1 , S4 , and S7 are linked via the property

S4 + 2S7 = S1 . (49)

As a result, Υ is a cubic tensor and one can easily compute the effective
tensor Ce from the results of Appendix A.

• The effective bulk modulus κe is determined by

f
κe = κm − (50)
1
κm −κi
− 3κm9S3
+4µm

• The first and the second effective shear modulus µe and µe∗ are deter-
mined by

f (µi − µm ) f (µi − µm )
µe = µm − , µe∗ = µm − (51)
1−β 1 − β∗

where the coefficients β and β ∗ are defined as

2(µi − µm ) (3κm + 4µm )S1 − 2(3κm + µm )S7

β= .
µm 3κm + 4µm
6(µ − µ ) µ S + 3 (κ m + µm ) S7
β∗ =
i m m 4
. (52)
µm 3κm + 4µm

16
Figure 2: Unit cell of cubic lattice structures (from left to right: simple cubic, body-
centered cubic, and face-centered cubic)

It is worth noting that these expressions are identical in form to those given by
Cohen (2004). Using this latter work, the terms depending on the geometry
and concentration are given by

1−f
S1 = , (53)
3
1
S4 = [1 − f − 2(A1 f + A2 f 5/3 )], (54)
5
1
S7 = [1 − f + 3(A1 f + A2 f 5/3 )], (55)
15

where A1 and A2 are coefficients depending on the geometry that character-

ize the kind of lattice up to second order. They were computed by Cohen
(2004). This shows the very similarity between the present approach and
that of Cohen (2004) despite the difference in nature. It is clear that the
result of Cohen (2004) is easier to use than ours, because the expressions for
S1 , S4 and S7 are the same for any cubic lattice of the same kind. However,
as shown subsequently, the previously described results and their expressions
in closed forms that are described thereafter can be extended readily to other
geometries of the microstructure.

Equations (50–51) can be further simplified by evaluating the lattice sums

S1 , S4 , and S7 analytically using the method described in subsection 3.1 with

17
grid density ρ = a3 /8π 3 R3 . It has been previously reported that S7 can be
easily computed from S1 and S4 via relation (49). For BCC and FCC ar-
rangements, due to the fluctuation of α (see Eqs. 47 and 48), H(η) in the
integral (34) should be replaced by its average H̄(η). Generally, we have

H̄(η) = ᾱf φ(η) (56)

where ᾱ = 1(SC), 21 (BCC) and 14 (FCC). However, we note that the composite
coefficient
1
ρᾱf = 2 (57)
2π
is independent of the microstructure, leading to a unique expression for the
long wave integral ∫ ∞
3
ρ H̄(η)η 2 dη = 2 J(ηc ) (58)
ηc 2π
Assuming ϵ = 2πR/a, the expressions of S1 and S4 for different cubic sys-
tems are finally presented by keeping the first terms of the lattice sums. The
numerical evidence shows that keeping four to six leading terms produces
satisfying results. The terms that are kept corresponding to the value of ηc
are given as follows:

- Simple cubic system (SC) with ηc = 2ϵ

2 √ 8 √
S1 = J(2ϵ) + f [2φ(ϵ) + 4φ( 2ϵ) + φ( 3ϵ) + 2φ(2ϵ)],
π 3
6 √ 8 √
S4 = J(2ϵ) + f [2φ(ϵ) + 2φ( 2ϵ) + φ( 3ϵ) + 2φ(2ϵ)]. (59)
5π 9

18
 √
- BCC with ηc = 12ϵ

2 √ √ √ √
S1 = J( 12ϵ) + f [4φ( 2ϵ) + 2φ(2ϵ) + 8φ( 6ϵ) + 4φ( 8ϵ) +
π
√ 8 √
+8φ( 10ϵ) + φ( 12ϵ)],
3
6 √ √ √ √
S4 = J( 12ϵ) + f [2φ( 2ϵ) + 2φ(2ϵ) + 4φ( 6ϵ) + 2φ( 8ϵ) +
5π
164 √ 8 √
+ φ( 10ϵ) + φ( 12ϵ)]. (60)
25 9

- FCC with ηc = 4ϵ

2 8 √ √ √
S1 = J(4ϵ) + f [ φ( 3ϵ) + 2φ(2ϵ) + 4φ( 8ϵ) + 8φ( 11ϵ)
π 3
8 √
+ φ( 12ϵ) + 2φ(4ϵ)],
3
6 8 √ √ 664 √
S4 = J(4ϵ) + f [ φ( 3ϵ) + 2φ(2ϵ) + 2φ( 8ϵ) + φ( 11ϵ)
5π 9 121
8 √
+ φ( 12ϵ) + 2φ(4ϵ)]. (61)
9

The maximal relative difference between the analytical expressions and the
full lattice sums is 2% (see Fig. 3). For FCC and BCC, more leading terms
than in the SC case are retained to achieve the high precision level, due to
the fluctuation of the coefficient α. Using four leading terms for FCC and
BCC cases can result in a simpler expression but a higher maximal error (up
to 8%). Even if this error was found to have a small impact on the effective
properties in the case of thermal conduction, we shall use the accurate ex-
pressions (60,61) in the remaining part of the paper.

As mentioned earlier (see Eq. 44), when the ratio R/a → 0 and f → 0,
we obtain the limit that is independent of the microstructure and volume

19
 0.45
S1 (full sum)
Figure 0.4
3: Comparison between the full lattice sums withS4resolution
(full sum) −128 < ni ≤ 128
and Eq.
0.35 (59) for the simple cubic case. At f = 0, the lattice
S 7
sums
(full sum) take the limit values
S1 = 1/3, S4 = 1/5, and S7 = 1/15. S (present)
1
0.3
S4 (present)
0.25 S (present)
Sum

fraction f :
0.2

0.15
2 1 6 1
0.1 S1 = J(0) = , S4 = J(0) = . (62)
π 3 5π 5
0.05

This interesting
0 property can be explained in relation to the Eshelby tensor
0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5
for periodic spheres discussed in fraction
Volume Section f 2.3 (which becomes equal to the
classical Eshelby tensor in the dilute case), as also observed in Fig. 3.

3.3. Orthorhombic arrangements of spheres

As will be seen thereafter in the numerical examples, the previously de-
scribed closed-form solution for the effective properties of cubic lattices leads
to results that are very similar to those of Cohen (2004). In addition, this
earlier solution is expressed by simple expressions of the effective elasticity
tensor. Thus, in this section and the following, we show that the solution
described in the previous subsection can be readily extended to other kinds

20
of periodic cells.

First, the case of an array of spheres with orthorhombic symmetry is consid-

ered, the center of the spheres being lattice points of an orthorhombic lattice
characterized by a1 , a2 , a3 , with a1 ̸= a2 ̸= a3 ̸= a1 . In this case, the effective
elasticity tensor is again given by expression (21). It can be seen that all
tensors appearing in this expression are isotropic except the orthotropic ten-
sor Υ that is characterized by the symmetry of the lattice. Expressions for
all effective elasticity components are explicit, but cumbersome, as a 3 × 3
matrix needs to be inverted, except for the shear components that are listed
below:

f µm
e
Cjkjk = µm − ,
µm
µm −µi
− (Sj + Sk ) + 4 λλmm+2µ
+µm
m
S i+6

i, j, k = 1, 2, 3, i ̸= j ̸= k ̸= i. (63)

Regarding the computation of the lattice sum Si , it is clear that the continu-
ous integral estimation of the remaining series is identical to the cubic lattice
case, except for the density that is now given by the general expression (32).
This implies that the contribution of longer wave vectors to the sum does
not contribute to the orthotropic anisotropy of the material. The related
anisotropic terms are contained in the complementary finite sum, and the
number of independent terms to keep in the finite sum is usually higher than
in the case of the cubic lattice.

The infinite sum in (63) can be computed as follows:

2 ∑ 2 ∑
Si = J(ηc ) + f η̄i2 αφ(η), Si+6 = J(ηc ) + f η̄j2 η̄k2 αφ(η),
π η<η
5π η<η
c c

i, j, k = 1, 2, 3, i ̸= j ̸= k ̸= i (64)

21
The lattice sums S4 , S5 , and S6 contributing to the effective components
e e e
other than C2323 , C1313 , and C1212 are given as follows for completeness:

6 ∑
Si+3 = J(ηc ) + f η̄i4 αφ(η), i = 1, 2, 3. (65)
5π η<η c

The coefficient α is a function of n1 , n2 , n3 with the same meaning as in

(47,48). It can be translated to the SO, BCO, and FCO cases.

3.4. Orthorhombic arrangements of ellipsoids

For an ellipsoid of principal radii R/χ1 , R/χ2 , and R/χ3 located at xc ,
we know the elementary result
∫
sin η ′ − η ′ cos η ′ iξ.xc 4πR3
eiξ.x dx = 3Ve e , Ve = ,
Ve η3 3χ1 χ2 χ3
η ′ = (η12 /χ21 + η22 /χ22 + η32 /χ23 )1/2 , η = Rξ. (66)

Again, the procedure described in subsection 3.1 is used to obtain the ana-
lytical solution. Assuming that the ellipsoid is located at the center of the
rectangular cuboid of dimension a1 × a2 × a3 and performing the transfor-
mation η1′ = η1 /χ1 ,η2′ = η2 /χ2 and η3′ = η3 /χ3 , the lattice sums can now be
evaluated by

3 ∑ ηi′2 χ2i
′
Si = 2 Ti J(ηc ) + f ′2 2 ′2 2 ′2 2
αφ(η ′ ),
2π η χ + η2 χ2 + η3 χ3
η ′ <ηc′ 1 1

3 ∑[ ηi′2 χ2i
]2
′
Si+3 = 2 Ti+3 J(ηc ) + f ′2 2 ′2 2 ′2 2
αφ(η ′ ),
2π η ′ <η ′
η χ
1 1 + η2 2 χ + η3 3χ
c

3 ∑ ηj′2 χ2j ηk′2 χ2k

′
Si+6 = 2 Ti+6 J(ηc ) + f ′2 2 ′2 2 ′2 2 2
αφ(η ′ ),
2π η ′ <η ′
[η χ
1 1 + η 2 2χ + η3 3χ ]
c

i, j, k = 1, 2, 3, i ̸= j ̸= k ̸= i. (67)

22
The surface integrals on the unit sphere surface Ti are defined by
∫
(χi ηi′ )2 dS
Ti =
|η ′ |=1 [(χ1 η1′ )2 + (χ2 η2′ )2 + (χ3 η3′ )2 ]
∫
(χi ηi′ )4 dS
Ti+3 =
|η ′ |=1 [(χ1 η1′ )2 + (χ2 η2′ )2 + (χ3 η3′ )2 ]2
∫
(χj ηj′ )2 (χk ηk′ )2 dS
Ti+6 = ′ 2 ′ 2 ′ 2 2
,
|η ′ |=1 [(χ1 η1 ) + (χ2 η2 ) + (χ3 η3 ) ]
i, j, k = 1, 2, 3, i ̸= j ̸= k ̸= i. (68)

Regarding the cutoff radius ηc′ in function of ϵ, it can be set to the same
value as ηc in (59) with ϵ = max(2πR/χ1 a1 , 2πR/χ2 a2 , 2πR/χ3 a3 ). It is
interesting to note that if χ1 a1 = χ2 a2 = χ3 a3 , that is, if the ellipsoid and
the unit cell have the same aspect ratio, the grid η ′ is a cubic grid and the
leading terms kept in the series are as simple as for the cubic lattice cases
(59–61). Figure 4 is an example showing that the analytical expression for
the FCO arrangement is highly accurate, compared to the full computation
of the lattice sums.
It can be shown that T1 , T2 , and T3 are expressed as functions of elliptical
integrals (see, e.g., Eshelby, 1957; Mura, 1987). Assuming that χ1 < χ2 < χ3 ,
these integrals read

4πχ31 χ22 χ3
T1 = √ [F (θ, k) − E(θ, k)]
(χ22 − χ21 ) χ23 − χ21
[ √ ]
4πχ1 χ22 χ33 1
T3 = 2 √ χ3 − χ1 − E(θ, k)
2 2
(χ3 − χ22 ) χ23 − χ21 χ2
T2 = 4π − T1 − T3 (69)

23
 S (full sum)
1
S4 (full sum) −128 < ni ≤ 128 and
Figure 4: Comparison between the full lattice sums with resolution
0.25
Eqs. (67,68) for the FCO case. The spheroids and the unitS cell have the same dimension
(full sum)
7
ratio a × a × 0.5a and R × R × 0.5R. The cutoff radius is ηc′S =(present)
4ϵ and ϵ = 2πR/a. For the
1
sake of0.2
clarity, only the sums S1 , S4 and S7 are presented. S (present)
4
S (present)
Sum

7
0.15

where0.1
∫ θ ∫ θ √
dw
0.05
F (θ, k) = √ , E(θ, k) = 1 − k 2 sin2 wdw
0 1− k2 2
sin w 0
√
0 √ 0.15 0.2 0.25 0.3 20.35 2 0.4 2 0.45
0 0.05 0.1
χ (χ − χ1 ) 0.5
θ = arcsin 1 − χ21 /χ 2
Volume
3, k = f 32 22
fraction . (70)
χ2 (χ3 − χ21 )

In the general case, there are no available explicit analytical expressions for
the remaining integrals Ti , i = 4, 5, .., 9 except for spheroidal inclusions. Ob-
viously, only Ti , i = 4, 5, 6 are necessary, due to the general relation between
Si . It can be shown that each of these coefficients can be obtained from one
simple scalar integral, which is in any case easier to compute than the full
lattice sum. The related scalar integrals are provided in Appendix C.

For the case of spheroids where the axes of cylindrical symmetry are ori-

24
ented along the x3 axis, we can set χ1 = χ2 = 1 and χ3 = χ and evaluate
analytically Ti in spherical coordinates (see Appendix B). The final results
are given as follows:

2π 4πχ2
T1 = T2 = (δχ2
− γ), T3 = (γ − δ)
γ3 γ3
3π
T4 = T5 = 5 (χ2 (χ2 − 4)δ + γ(χ2 + 2)),
4γ
2
2πχ
T6 = (γ(2χ2 + 1) − 3χ2 δ)
γ5
πχ2
T7 = T8 = 5 ((χ2 + 2)δ ′ − 3γ),
γ
π
T9 = 5 (χ2 (χ2 − 4)δ + γ(χ2 + 2)), (71)
4γ
√
with γ = χ2 − 1 and δ = arctan γ for oblate inclusions (χ > 1) and
√
γ = 1 − χ2 and δ = tanh−1 γ for prolate inclusions (χ < 1).

4. Closed-form expressions for random distributions of spheres or

ellipsoids

In the previous sections, we considered only the cases of lattice distri-

butions of spheres and ellipsoids. This section is devoted to random distri-
butions of spheres and ellipsoids. Therefore, the periodic cell now contains
a random distribution of aligned inclusions. To proceed, the ergodic media
hypothesis (Torquato, 2001) is adopted, that is, the ensemble average results,
notation ⟨...⟩ens , are identical results for one sample in the infinite-volume
limit. This assumption guarantees the existence and uniqueness of the effec-
tive tensor Ce .

Obviously, the ergodic assumption is not sufficient to estimate the effective

properties. It is necessary to introduce also probabilistic assumptions on the
distributions of inclusions to reach this objective. These assumptions can be

25
based on the introduction of probability density functions, correlation func-
tions in the real space, but also on quantities defined in Fourier space. Let
us consider in a first step the notions defined in the real space that sustain
our work. Following Torquato (2001), it is convenient to introduce the prob-
ability density function PN such that PN (r1 , r2 , ...)dr1 dr2 ... represents the
probability of finding the center of inclusion 1 in dr1 , the center of inclusion
2 in dr2 ,... By a convenient partial integration of this probability function,
the pair correlation function can be simplified into the radial distribution
function g2 (r12 ) (with r12 = ||r2 − r1 ||) under the assumption of statistical
isotropy.

However, our method uses quantities defined in Fourier space and it is pos-
sible to define also ensemble average on functions defined in Fourier space.
So, in a first step, such functions will be defined. In a second step the link
between such ensemble averages and the radial distribution function g2 de-
fined by using the full probabilistic machinery will be provided.

As a result of the ergodic assumption, from (20), we can find that the fol-
lowing relation must hold true:

⟨ϵ∗ (x)⟩Ω = ⟨⟨ϵ∗ (x)⟩Ω ⟩ens if V → ∞ (72)

Taking the ensemble average of (17) and accounting for (72), we find that

Cm : ⟨ϵ∗ (x)⟩Ω ≃ (Cm − Ci ) : [E + ⟨Υ⟩ens : Cm : ⟨ϵ∗ (x)⟩Ω ] , (73)

Consequently, the final expression of the effective tensor (21) for random
media is slightly modified, Υ now being replaced by its ensemble average
⟨Υ⟩ens . In what follows, we shall examine the properties of ⟨Υ⟩ens given by

26
the expression
∑
⟨Υ⟩ens = ⟨P (ξ)⟩ens Γm (ξ). (74)
ξ̸=0

Considering the integral over the volume of one particle Vp of arbitrary shape
located at xc and the definition of two functions F(ξ) and P(ξ)
∫
1
eiξ.x dx = F(ξ).eiξ.xc , P(ξ) = F(ξ)F(−ξ) (75)
Vp Vp2

Here, F(ξ) is the integral on the same volume Vp located at the origin. In the
scattering theory, F(ξ) and P(ξ) are called form factors and are equivalent
to the definitions I(ξ) and P (ξ) for the case of one particle. For random
distribution of particles, P (ξ) in (14–74) can be replaced by the ensemble
average ⟨P (ξ)⟩ens . The latter is linked to the form factor P(ξ), structure
factor S(ξ), and scattering intensity I(ξ) of the system via the relations

Vp
⟨P (ξ)⟩ens = I(ξ), I(ξ) = P(ξ)S(ξ),
⟨ VN ⟩
1 ∑ iξ·xi ∑ −iξ·xi
N
S(ξ) = e e . (76)
N i i ens

In (76), we assume that there are N particles of characteristic dimension R

in a cubic cell of dimension a and xi is the location of the particle number i in
the cell. It is interesting to note that I(ξ) can be obtained experimentally by
scattering techniques. Conversely, theoretical results of I(ξ) for some ideal
systems are known, for example, those governed by Ornstein–Zernike (OZ)
equation and Percus–Yevick (PY) closure approximation. As long as I(ξ)
is known, the infinite lattice sums Si can also be computed. At the infinite
volume limit R/a → 0, taking R as the radius of the sphere having the same

27
volume as Vp (i.e., Vp = 4πR3 /3), those sums can be replaced by the integrals
∫ ∫
1 1
Si = 2 η̄i I(ξ)dη, Si+3 = 2 η̄i4 I(ξ)dη,
2
6π 6π
∫
1
Si+6 = 2 η̄j2 η̄k2 I(ξ)dη, η = Rξ,
6π
i, j, k = 1, 2, 3, i ̸= j ̸= k ̸= i. (77)

In the case where the particles are identical spheres of radius R with isotropic
distribution (see Fig. 5a), the associated isotropic quantities (iso as super-
script) I iso (ξ), P iso (ξ), and S iso (ξ) are all functions of the modulus ξ (or η).
In addition, the structure factor S iso (ξ) is related to the radial distribution
function g(r) via the relation
∫ ∞
sin ηr̄ r
S iso
(ξ) = 1 + 3f [g(r) − 1]r̄dr̄, r̄ = , (78)
0 η R

For nonoverlapping spheres g(r) = 0 when r̄ < 2, the following property

holds:
∫ ∞
[η cos η − sin η]2
sin(ηr̄)dη = 0, ∀r̄ ≥ 2. (79)
0 η5

As a result, the integrals Si admit the simple expressions

(1 − f ) (1 − f )
S1 = S2 = S3 = , S 4 = S5 = S6 = ,
3 5
(1 − f )
S7 = S8 = S9 = . (80)
15

28
Figure 5: Randomly distributed spheres (left) and ellipsoids (right). The distribution of
spheres of radius R is isotropic. The latter is obtained by scaling the former with ratios
1/χ1 , 1/χ2 , and 1/χ3 along directions 1, 2, and 3.

Substituting (80) back into (50, 51, 52), we obtain:

f
κe = κm − ,
1
κm −κi
− 3(1−f )
3κm +4µm
f
µe = µe∗ = µm − , (81)
1
µm −µi
− 6(κm +2µm )(1−f )
5µm (3κm +4µm )

Equations (81) show that, at the infinite volume limit, the effective mate-
rial is isotropic and that its effective elasticity tensor corresponds to the
Hashin–Shtrikman bound (Hashin and Shtrikman, 1963), or equivalently to
the Mori–Tanaka estimation for spherical inclusions (Mori and Tanaka, 1973).
Interestingly, this equivalence has been noted previously in the issue of heat
conduction (To et al., 2013) and has recently been rediscovered in linear elas-
ticity.

It is worth noting that the Mori–Tanaka estimate has been recovered along
the previous lines using statistical information, including the ergodic hypoth-
esis, although this estimate is usually presented by using different ad hoc as-
sumptions. Indeed, following Benveniste (1987), such an ad hoc assumption

29
can be expressed as ”if introducing a single inclusion into a homogeneous ma-
trix under boundary conditions corresponding to an overall strain ϵ0 results
in an average strain in the inclusion given by ϵI = T.ϵ0 , then introducing
a single inclusion into a deformed matrix having an average strain ϵM will
result in an average strain in the inclusion given by ϵI = T.ϵM .”

Our results can be compared with full results on random distributions of

particles. For example, Segurado and Llorca (2002) modeled the behaviour
of composites containing randomly distributed spherical inclusions using fi-
nite elements. For low concentrations of inclusions, their results are very close
to the Mori–Tanaka estimates, as predicted previously by our developments.
However, for higher concentrations where the inclusions strongly interact,
the present approach is not sufficiently accurate. In this case, higher-order
correlation functions can be used to address these issues, as done by Nguyen
et al. (2016) for the heat conduction phenomenon.

For anisotropic distributions of identical ellipsoids, there exists an analyt-

ical solution for the special case where the distribution of inclusions can be
obtained from an isotropic distribution by uniform dilatation transforma-
tion with coefficients 1/χ1 , 1/χ2 , and 1/χ3 along the three directions 1, 2,
and 3 (see Fig. 5b). In other words, not only are the spherical inclusions
transformed into ellipsoids of dimensions R/χ1 , R/χ2 , and R/χ3 but their co-
ordinates are scaled with the same ratios as well. In doing so, we also obtain a
system containing nonoverlapping ellipsoids with the same volume fraction f .
To respect the definition of the characteristic length R defined in the begin-
ning of the present subsection, χ1 χ2 χ3 = 1 is assumed. The anisotropic (ani
as superscript) and isotropic statistical quantities associated to two systems

30
are interconnected, say
(√ )
P (ξ) = P
ani iso 2 2 2 2 2 2
ξ1 /χ1 + ξ2 /χ2 + ξ3 /χ3 ,
(√ )
ani iso 2 2 2 2 2 2
S (ξ) = S ξ1 /χ1 + ξ2 /χ2 + ξ3 /χ3 ,
(√ )
I (ξ) = I
ani iso 2 2 2 2 2 2
ξ1 /χ1 + ξ2 /χ2 + ξ3 /χ3 . (82)

Again, by changing the variable, we obtain

∫ ∞
Ti
Si = 2 I iso (ξ ′ )dη ′ ,
6π 0
√
ξ ′ = ξ12 /χ21 + ξ22 /χ22 + ξ32 /χ23 , η ′ = Rξ ′ . (83)

Combined with the results for the isotropic case, Si admits the simple form

Ti
Si = (1 − f ), (84)
4π

which results again in a closed-form expression of the effective constants.

5. Numerical examples and comparisons

To illustrate the method described previously, we first consider numerical

applications in the case of cubic microstructures. The closed-form solutions
derived in this work are used to compute the effective properties, and these
results will be compared with exact semi-analytical, numerical, and approx-
imate solutions from the literature including

- Numerical results with NIH approximation (Iwakuma and Nemat-Nasser,

1983; Nemat-Nasser et al., 1982)
- CM-type approximation in electrostatics, which was applied to linear elas-
ticity and carefully compared with other types of solutions by Cohen (2004)

31
- Exact series expansion of elastic polarizabilities and approximation up to
third-order correlation (To) (Torquato, 1997, 1998)
- Semi-analytical solution based on periodic singular distribution (SL) (San-
gani and Lu, 1987), using an expansion onto a basis of periodic functions and
projection onto a basis of spherical harmonics.
- The multipole expansion (ME) method, which can provide accurate solu-
tions by keeping a high number of terms in the series (Kushch, 2013).

As an example, we choose a specific composite with the same Poisson’s ratio

in both the matrix and inclusions νi = νm = 0.3. Next, the dimensionless
e
effective coefficient C1111 /µm is computed at different contrast ratios µi /µm .
The results tabulated in Table 2 show that the closed-form solution, which is
based on NIH approximation, agrees very well with CM and ME solutions in
the range µi /µm ≤ 100 and f ≤ 0.5. At higher contrasts, some discrepancies
are observed between the approaches. However, as noted in previous works
(Hoang and Bonnet, 2013; Cohen, 2004), both methods (NH and CM) fail in
this range. In this case, only the ME solution provides an accurate solution
(Kushch, 2013).

Table 1: Ratio of the effective elastic constants C1111 /µm of a simple cubic array of
sphere, νi = νm = 0.3. Comparison between the closed-form solution of the present work
(PR), Clausius–Mossotti (CM)-type solution from Cohen (2004), and solutions based on
multipole expansion (ME) from Kushch (2013).

µi /µm f = 0.1 f = 0.3 f = 0.5

PR CM ME PR CM ME PR CM ME
0 2.794 2.799 2.799 1.852 1.836 1.818 1.137 1.148 1.048
10 4.097 4.101 4.102 5.913 5.880 5.930 8.580 8.645 9.646
100 4.226 4.248 4.248 6.720 6.758 6.887 10.94 11.453 17.71

The first extreme case is related to void inclusions. By making µi /µm → 0,

32
the limit of (50–51) yields the following simple expressions:

κe f (1 − νm )
=1− , (85)
κm 1 − νm − (1 + νm )S3
µe f (1 − νm )
=1− , (86)
µm 2(S2 − (1 − νm )S3 ) + 1 − νm
µe∗ f (1 − νm )
=1− . (87)
µm 1 − νm − (1 − 2νm )S1 − (5 − 4νm )S2

In the other extreme case of rigid inclusions, that is, µi /µm → ∞, the effective
properties can be determined as follows:

κe f (1 − νm )
=1+ , (88)
κm (1 + νm )S3
µe f (1 − νm )
=1+ , (89)
µm 2((1 − νm )S3 − S2 )
µe∗ f (1 − νm )
=1+ . (90)
µm (1 − 2νm )S1 + (5 − 4νm )S2

In these examples, the Poisson’s ratio of the matrix is fixed at νm = 0.3. The
resulting effective properties are plotted in Fig. 6 for spherical voids and Fig.
7 for rigid inclusions.

All figures show that the results obtained by our closed-form solution com-
pare well with Cohen’s results. For voids, a discrepancy with Sangani and
Lu’s results is observed, as with Cohen’s results. However, in this case, Cohen
(2004) expressed doubts with Sangani and Lu’s results. For rigid inclusions,
our results coincide in the range of concentrations up to 0.5 with Cohen’s
results. Cohen compared these results with a few semi-analytical or approxi-
mate models, which were found to be accurate in this range of concentrations
of inclusions.

In the case of orthotropic materials, one considers distributions of ellipsoids

33
 1

Effective incompressibility modulus κe / κm

Present
0.9 Cohen (2004)
Torquato (1997)
0.8 Sangani and Lu (1987)

0.7

0.6

0.5

0.4

0.3

0.2
0 0.1 0.2 0.3 0.4 0.5
Volume concentration of inclusion

Figure 6: Effective bulk modulus κe /κm for a simple cubic array of void spheres.

6
Second effective shear modulus µe*/ µ m

Present
Cohen (2004)
5

1
0 0.1 0.2 0.3 0.4 0.5
Volume concentration of inclusion

Figure 7: Effective shear modulus µe∗ /µm for a simple cubic array of rigid spheres.

located at central positions of a cubic lattice. Table 2 shows the compari-

son between our estimation of Si , i = 1..6 and the results obtained by the
computation of the full lattice sums (Iwakuma and Nemat-Nasser, 1983) for

34
two examples characterized by χ1 = χ2 and the values of χ3 /χ1 = 4/3 or
χ3 /χ1 = 2. This comparison proves that the closed-form expressions obtained
in the previous section accurately reproduce the values of Si computed from
the full lattice sum.

Table 2: Comparison between Si results of the present work (present) and those from
Table 3 in Iwakuma and Nemat-Nasser (1983) (IN). The unit cell is cubic of dimension a
and the aspect ratio χ1 : χ2 : χ3 of the spheroids in SO arrangement are 3 : 3 : 4 (case 1)
and 1 : 1 : 2 (case 2). Note that results presented by Iwakuma and Nemat-Nasser (1983)
are obtained from full lattice sums with resolution −50 < ni < 50.

f 0.05 0.10 0.15 0.20 0.25 0.30 0.35

Case 1
S1 = S2 Present 0.268 0.240 0.230 0.216 0.196 0.171 0.148
IN 0.273 0.255 0.236 0.217 0.196 0.175 0.154
S3 Present 0.397 0.387 0.376 0.363 0.353 0.345 0.340
IN 0.392 0.380 0.369 0.359 0.351 0.343 0.337
S4 = S5 Present 0.172 0.172 0.170 0.164 0.153 0.139 0.123
IN 0.170 0.171 0.168 0.161 0.152 0.139 0.124
S6 Present 0.279 0.286 0.292 0.294 0.295 0.297 0.299
IN 0.275 0.281 0.286 0.290 0.293 0.295 0.297
Case 2
S1 = S2 Present 0.218 0.196 0.174 0.148 0.120
IN 0.216 0.195 0.172 0.147 0.119
S3 Present 0.514 0.505 0.502 0.502 0.508
IN 0.507 0.500 0.496 0.497 0.504
S4 = S5 Present 0.129 0.126 0.120 0.107 0.088
IN 0.127 0.126 0.119 0.107 0.089
S6 Present 0.397 0.407 0.422 0.439 0.458
IN 0.391 0.403 0.418 0.436 0.455

Finally, we consider examples concerning orthorhombic and random arrange-

ments of spheroids. For orthorhombic arrangements, the unit cell has the
dimensions a × a × 0.5a and the spheroids has axes R × R × 0.5R. Re-
garding the constitutive materials, the Poisson’s ratios of both materials are

35
νi = νm = 0.3 and the Young modulus ratio is Ei /Em = 10. Figures 8 and
9 display the results of analytical solutions for the shear constant obtained
in this paper. It is noted that all the curves are relatively close at a small
volume fraction (f ≤ 0.15), which can be explained by the independent be-
havior of each inclusions. In this range, the dilute estimation is valid, that
is, the effective properties only depend on the volume fraction f . At higher
f , the interaction of the inclusions is significant and their relative positions
in the matrix are not negligible. Although curves start to deviate from each
other, the FCO curves are still close to the random curves. This observation
suggests that the structure of the FCO arrangements is comparable to the
random structures proposed previously.

3.5

FCO
3 BCO
SO
Random
2.5
Ce1212/µm

2
e
Figure 8: Dimensionless effective shear moduli C1212 /µm for FCO, BCO, SO, and random
arrangements of spheroids.
1.5

1
0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5
6. Concluding remarks Volume fraction f

Effective elasticity tensors of two-phase matrix inclusion composite mate-

rials have been obtained when the inclusions are spherical or ellipsoidal and
distributed either along the sites of an orthorhombic crystal system or ran-
domly. The basis of the analysis is the eigenstrain integral equation and the

36
 2.5
FCO
BCO
SO
2 Random
Ce1313/µm

e
Figure 9: Dimensionless effective
1.5 shear moduli C1313 /µm for FCO, BCO, SO, and random
arrangements of spheroids.

NIH-type estimation (Nemat-Nasser

1 et al., 1982). We have shown that this
0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5
Volumetensor”
estimation is related to the use of the ”periodic Eshelby fraction f which allows

to recover the mean strain within the inclusions when a constant eigenstrain
is applied over the periodic set of inclusions.

A first validation is provided in the specific case of cubic materials: all our
results fully reproduce those of Cohen (2004), which were obtained by a per-
turbation method, in relation to the CM approximation of electrostatics. We
have shown that the expressions of the components of the elasticity tensor
reproduce those of Cohen (2004), in the sense that the lattice sums Si of
the NIH approximation can be obtained using the coefficients obtained by
Cohen’s perturbation method.

We have shown that our method can also estimate the orthotropic elasticity
tensor in the case of arrangements obtained by parallelelipedic lattice cells of
spheres or ellipsoids. All closed-form solutions are fully extended to the case
of cells containing spheroids. For the case of fully ellipsoidal inclusions, only

37
lattice sums S1 , S2 , S3 can be obtained in a closed form, whereas S4 , S5 , S6 are
provided by simple scalar integrals, which are in any case easier to compute
than full lattice sums. Our results compare well with those obtained by the
full computation of lattice sums as shown by Iwakuma and Nemat-Nasser
(1983).

In the case of a random distribution of spheres, the statistical connection

of the effective tensor to the structure factor is shown. A closed-form expres-
sion is obtained in the infinite volume limit for an isotropic distribution of
spheres, which reproduces the Mori–Tanaka estimate. Although this approx-
imation is usually obtained by ”ad hoc” assumptions, our work shows that
this estimate can be based on precise statistical information, thus extending
a result obtained similarily in the case of conduction through a random dis-
tribution of spheres (To et al., 2013).

As a final remark, the present method can also be extended to a more general
case where the principal axes of ellipsoids (or spheroids) do not necessarily
coincide with the axes x1 , x2 , and x3 of a rectangular unit cell, as considered
in this paper. In such situations, it is sufficient to use the linear transforma-
tion L as discussed in section 2; that is, the grid η ′ is obtained by rotation
and dilatation of the grid η. However, in this case, the effective elasticity
tensor is no longer orthotropic, and it leads to a larger number of lattice
sums Si to be computed.

References

Benveniste, Y., 1987. A new approach to the application of Mori-Tanaka’s

theory in composite materials. Mech. Mater. , 147–157.

Benveniste, Y., Milton, G., 2003. New exact results for the effective electric,
elastic, piezoelectric and other properties of composite ellipsoid assem-
blages. J. Mech. Phys. Solids 51, 1773–1813.

38
Bohren, C.F., Huffman, D.R., 1998. Absorption and scattering of light by
small particles. John Wiley & Sons, New York.

Christensenn, R.M., 1979. Mechanics of composite materials. Wiley, N.Y..

Christensen, R., Lo, K., 1979. Solutions for effective shear properties in three
phase sphere and cylinder models. J. Mech. Phys. Solids 27, 315–330.

Cohen, I., 2004. Simple algebraic approximations for the effective elastic
moduli of cubic arrays of spheres. J. Mech. Phys. Solids 52, 2167–2183.

Eshelby, J., 1957. The determination of the elastic field of an ellipsoidal

inclusion, and related problems. Proc. R. Soc. A 241, 376–396.

Eyre, D.J., Milton, G.W., 1999. A fast numerical scheme for computing the
response of composites using grid refinement. Eur. Phys. J. Appl. Phys. 6,
41–47.

Hashin, Z., Shtrikman, S., 1963. A variational approach to the theory of the
elastic behaviour of multiphase materials. J. Mech. Phys. Solids 11, 127 –
140.

Hill, R., 1965. A self-consistent mechanics of composite materials. J. Mech.

Phys. Solids 13, 213 – 222.

Hoang, D.H., Bonnet, G., 2013. Effective properties of viscoelastic heteroge-

neous periodic media: An approximate solution accounting for the distri-
bution of heterogeneities. Mech. Mater. 56, 71–83.

Hunter, R.J.,White, L.R., 2001. Foundations of colloid science. Oxford

University Press, Oxford.

Iwakuma, T., Nemat-Nasser, S., 1983. Composites with periodic microstruc-

ture. Comput. Struct. 16, 13–19.

39
Kaminski, M., 1999. Boundary element method homogenization of the peri-
odic linear elastic fiber composites. Eng. Anal. Bound. Elem. 23, 815–823.

Kaminski, M., 2005. Multiscale homogenization of n-component composites

with semi-elliptical random interface defects. Int. J. Solids Struct. 42,
3571–3590.

Kushch, V., 1997. Microstresses and effective elastic moduli of a solid rein-
forced by periodically distributed spheroidal particles. Int. J. Solids Struct.
34, 1353–1366.

Kushch, V., 2013. Micromechanics of composites. Multipole expansion ap-

proach. Butterworth-Heinemann, New York..

Liu, Y., Nishimura, N., Otani, Y., Takahashi, T., Chen, X., Munakata, H.,
2005. A fast boundary element method for the analysis of fiber-reinforced
composites based on a rigid-inclusion model. J. Appl. Mech.-Trans. ASME
72, 115–128.

Luciano, R., Barbero, E.J., 1994. Formulas for the stiffness of composites
with periodic microstructure. Int. J. Solids Struct. 31, 2933–2944.

Michel, J., Moulinec, H., Suquet, P., 1999. Effective properties of composite
materials with periodic microstructure: a computational approach. Com-
put. Meth. Appl. Mech. Eng. 172, 109–143.

Milton, G.W., 2002. The theory of composites. Cambridge Univ. Press,

Cambridge.

Monchiet, V., Bonnet, G., 2012. A polarization-based FFT iterative scheme

for computing the effective properties of elastic composites with arbitrary
contrast. Int. J. Numer. Meth. Eng. 89, 1419–1436.

40
Mori, T., Tanaka, K., 1973. Average stress in matrix and average elastic
energy of materials with misfitting inclusions. Acta. Metall. Mater. 21,
571–574.

Mura, T., 1987. Micromechanics of defects in solids. Kluwer Academic

Publisher, New York.

Nemat-Nasser, S., Iwakuma, T., Hejazi, M., 1982. On composites with peri-
odic structure. Mech. Mater. 1, 239–267.

Nguyen, M.T., Monchiet, V., Bonnet, G., To, Q.D., 2016. Conductivity esti-
mations of spherical particle suspension based on triplet structure factor.
Phys. Rev. E 93, 022105.

Nunan, K., Keller, J., 1984. Effective coefficients of a composite containing

spherical inclusions in a periodic array. J. Mech. Phys. Solids 32, 259–280.

Rössler, E., 2009. Solid state theory: an introduction. Springer Verlag,

Berlin.

Sanchez-Palencia, E., 1974. Comportements local et macroscopique d’un

type de milieux physiques hétérogènes. Int. J. Eng. Sci. 12, 331–351.

Sanchez-Palencia, E., 1980. Non-homogeneous media and vibration theory.

Lecture Notes in Physics, Volume 127. Springer Verlag, Berlin.

Sangani, A.S., Lu, W., 1987. Elastic coefficients of composites containing

spherical inclusions in a periodic array. J. Mech. Phys. Solids 35, 1–21.

Segurado, J., Llorca, J., 1987. A numerical approximation to the elastic

properties of sphere-reinforced composites. J. Mech. Phys. Solids 50, 2107–
2121.

To, Q.D., Bonnet, G., 2015. Conductivity of periodic composites made of

matrix and polydispersed aggregates. Phys. Rev. E 91, 023206.

41
To, Q.D., Bonnet, G., To, V.T., 2013. Closed-form solutions for the effective
conductivity of two-phase periodic composites with spherical inclusions.
Proc. R. Soc. A 469, 1471–2946.

Torquato, S., 1997. Effective stiffness tensor of composite media. I. Exact

series expansions. J. Mech. Phys. Solids 45, 1421–1448.

Torquato, S., 1998. Effective stiffness tensor of composite media : II. Appli-
cations to isotropic dispersions. J. Mech. Phys. Solids 46, 1411 – 1440.

Torquato, S., 2001. Random heterogeneous materials: microstructure and

macroscopic properties. Springer Verlag, Berlin.

Tsukrov, I. , Kachanov, M., 2000. Effective moduli of an anisotropic material

with elliptical holes of arbitrary orientational distribution. Int. J. Solids
Struct. 37, 3571–3590.

Appendix A. Algebraic operations between cubic tensors

Algebraic operations between cubic tensors can be made simple using the
base B = {I, i ⊗ i, N} and Table A.3 where I is the fourth-order identity
tensor and i the second-order identity tensor. Tensor N in this case is equal
to

N = e1 ⊗ e1 ⊗ e1 ⊗ e1 + e2 ⊗ e2 ⊗ e2 ⊗ e2 + e3 ⊗ e3 ⊗ e3 ⊗ e3 ,
(A.1)

in which e1 , e2 , and e3 are the unit vectors along directions 1, 2, and 3 of the
coordinate system. The inversion of cubic tensors is also simple enough to

42
 Table A.3: Basic algebraic operations between tensors in base B

: I i⊗i N
I I i⊗i N
i ⊗ i i ⊗ i 3i ⊗ i i ⊗ i
N N i⊗i N

work with symbolic notations, for example,

λ 1
[λi ⊗ i + 2µI + αN]−1 = − i⊗i+ I−
(2µ + α)(3λ + 2µ + α) 2µ
α
− N (A.2)
2µ(2µ + α)

Thus, when combined with (10), the tensor Υ can be rewritten in base B as

Υ = AI + Bi ⊗ i + CN (A.3)

where the coefficients A, B, C are defined by the formulae:

2 (λm + µm ) 1
A=− S7 + S1 ,
µm (λm + 2µm ) µm
(λm + µm )
B=− S7 , (A.4)
µm (λm + 2µm )
(λm + µm ) 3 (λm + µm )
C=− S4 + S7 .
µm (λm + 2µm ) µm (λm + 2µm )

Only three lattice sums appear, due to the cubic symmetry (49). In addition,
due to (29), these sums are interdependent via the relation S4 + 2S7 = S1 .
As each phase is an isotropic material, the tensors Cm and (Ci − Cm )−1 , etc.
in (21) can be expressed in base B. This evidence leads to the final cubic

43
effective tensor Ce

Ce = λe i ⊗ i + 2µe I + αe N (A.5)

where the coefficients λe , µe , and αe can be used to determine the three

effective elasticity coefficients κe , µe∗ by

3λe + 2µe + αe 2µe + αe

κe = , µe∗ = (A.6)
3 2

The final explicit expressions of κe , µe , and µe∗ are given in (50) and (51).

Appendix B. Computation of Ti coefficients for spheroidal inclu-

sions

In the spherical coordinate system (r, θ, ϕ), the surface integrals have the
following forms:
∫ ∫ 2
sin2 ϕ sin2 θ sin θdθdϕ χ cos2 θ sin θdθdϕ
T1 = T2 = , T 3 = ,
sin2 θ + χ2 cos2 θ sin2 θ + χ2 cos2 θ
∫ ∫ 4
sin4 ϕ sin4 θ sin θdθdϕ χ cos4 θ sin θdθdϕ
T4 = T5 = , T6 = ,
[sin2 θ + χ2 cos2 θ]2 [sin2 θ + χ2 cos2 θ]2
∫ 2
χ cos2 θ sin2 θ sin2 ϕ sin θdθdϕ
T7 = T8 = ,
[sin2 θ + χ2 cos2 θ]2
∫
sin4 θ sin2 ϕ cos2 ϕ sin θdθdϕ
T9 = , θ ∈ [0, π], ϕ ∈ [0, 2π]. (B.1)
[sin2 θ + χ2 cos2 θ]2

44
Eliminating ϕ and making variable change t = cos θ, the integrals Ti can be
reduced to more tractable forms:
∫ 1 ∫ 1
(1 − t2 )dt χ2 t2 dt
T1 = T2 = π , T 3 = 2π ,
−1 (1 − t ) + χ t −1 (1 − t ) + χ t
2 2 2 2 2 2
∫ ∫ 1
3π 1 (1 − t2 )2 dt χ4 t4 dt
T4 = T5 = , T 6 = 2π ,
4 −1 [(1 − t2 ) + χ2 t2 ]2 −1 [(1 − t ) + χ t ]
2 2 2 2
∫ 1 ∫
χ2 (1 − t2 )t2 dt π 1 (1 − t2 )2 dt
T7 = T8 = π , T9 = .
−1 [(1 − t ) + χ t ] 4 −1 [(1 − t2 ) + χ2 t2 ]2
2 2 2 2

(B.2)

Finally, the integration leads to the closed-form expressions given in (71).

Appendix C. Computation of tensors Ti , i = 4..6 for fully ellip-

soidal inclusions

For fully ellipsoidal inclusions, coefficients Ti , i = 4..6 are given by the

scalar integrals
∫
χ4i π/2
Ti = 4 4 Fi (ϕ)dϕ (C.1)
χ3 0
where Fi are the following functions:

M (1 − 4ζ) + ζ(2ζ + 1)Q 4

F4 = cos (ϕ)
ζ 4 Q5
F5 = tan4 ϕF4
P ζ(ζ + 2) − 3ζM
F6 = (C.2)
ζ 3 Q5

cos2 ϕχ21 +sin2 ϕχ22

where ζ = χ23
,
√
ζ −1
Q= (C.3)
ζ
M = tanh−1 (Q) (C.4)

45
 for ζ > 1,
√
1−ζ
Q= (C.5)
ζ
M = arctan(Q) (C.6)

for ζ < 1, whereas

16
F4 = cos4 (ϕ) (C.7)
5
F5 = tan4 ϕF4
2
F6 = (C.8)
5

for ζ = 1.

46

View publication stats