## Abstract

In this paper, we use approximate solutions of Nemat-Nasser *et al.* to estimate the effective conductivity of two-phase periodic composites with non-overlapping spherical inclusions. Systems with different inclusion distributions are considered: cubic lattice distributions (simple cubic, body-centred cubic and face-centred cubic) and random distributions. The effective conductivities of the former are obtained in closed form and compared with exact solutions from the fast Fourier transform-based methods. For systems containing randomly distributed spherical inclusions, the solutions are shown to be directly related to the static structure factor, and we obtain its analytical expression in the infinite-volume limit.

## 1. Introduction

Modelling composite materials and determining their effective properties from the microstructure have always been an active research area. Many analytical and numerical tools have been developed to solve the localization problem and predict the overall behaviour. Regarding the heat-conduction phenomena, analytical works are based on the potential theory and series of spherical harmonic functions [1–5]. Homogenization with numerical methods such as the finite-element method, the boundary-element method, the fast Fourier transform (FFT), etc. have been reported in numerous works [6–10]. For randomly heterogeneous materials, estimations and bounds for the effective properties are also established (see [11,12] and references therein).

In this paper, we estimate the effective conductivity for two-phase periodic composites with identical non-overlapping spherical inclusions. The estimations are based on the integral equation for **e*** and its dual form (integral equation on **q***),^{1} and the Nemat-Nasser–Iwakuma–Hejazi (NIH) approximation proposed by Nemat-Nasser *et al.* [13]. Previous works [14,15] have shown that the NIH approximation predicts the overall elastic constants for a large range of volume fractions *f* very well. The present work concerns the extension of the solutions to heat-conduction phenomena with applications in various practical cases. First, the cubic array arrangements (simple cubic (SC), body-centred cubic (BCC) and face-centred cubic (FCC)) are considered. The estimation of the conductivity involves infinite lattice sums that can be evaluated via continuous integrals and reduced to a fully closed form.

The second part of this work concerns the treatment of the random distribution case. Interestingly, the ingredient of the solutions is the static structure factor *S*(*ξ*) and the radial distribution function (RDF, *g*(*r*)). Two particular hard-sphere distributions are studied (the well-stirred and the Percus–Yevick (PY) distributions), and in the homogenization limit, the associated solutions converge towards a closed-form expression. Although the approximation is based only on the RDF, the results issued from this work are comparable with available solutions for suspensions of spherical particles by Bonnecaze & Brady [16,17]. In their works, a fully realistic particle distribution is adopted via Monte Carlo simulation and effective conductivity is then estimated using the moment expansion of the integral equation [16].

The paper is organized as follows. After the Introduction, §2 is dedicated to the formulation of the periodic thermal problem with integral equations, discussion of the NIH approximation and FFT-based numerical methods. Closed-form solutions and analytical expressions for effective conductivity of different composite structures are given in §3. Section 4 focuses on the application and comparison with FFT-based numerical solutions. Finally, conclusions and perspectives are given in §5.

## 2. Homogenization of two-phase periodic composites

### (a) Periodic thermal conduction problem and governing integral equations

We consider first a two-phase composite, periodic in three orthogonal directions *x*_{1},*x*_{2},*x*_{3} with periods *a*_{1},*a*_{2},*a*_{3}. Each phase is made of a homogeneous material whose thermal behaviour is governed by Fourier’s law,
2.1In (2.1), **q**(**x**) is the flux vector, **e**(**x**) the (minus) temperature gradient and **K**(**x**) the second-order local conductivity tensor. The latter can be either **K**^{M} or **K**^{I}, depending on which material, matrix or inclusion, is found at location **x**. For isotropic materials, the conductivity tensors have the following diagonal forms:
2.2with **I** being the identity tensor and *k*_{M},*k*_{I} the conductivities of both phases. In addition to (2.1), the flux **q** must satisfy the energy conservation equation
2.3To compute the effective conductivity **K**^{eff}, we apply a macroscopic temperature gradient **E** to the periodic cell and find its relation to the average flux **Q** of the cell volume *V* ,
2.4In (2.4) and from now on, we use the notation 〈*ϕ*〉_{V} to refer to the average of the quantity *ϕ* inside the bracket over the volume *V* ,
2.5The localization problem can be reduced to finding the *V* -periodic perturbation terms **e**^{per}, *T*^{per} and **q**^{per} given by the expressions
2.6Based on the reference material with conductivity **K**^{0}, we define a vector field **e*** by the formula
2.7Owing to the *V* -periodicity, **e**^{per}, *T*^{per}, **q**^{per}, **e*** admit the Fourier series representations
2.8in which **A** denotes **e**^{per}, *T*^{per}, **q**^{per},**e*** and denotes their Fourier transform , , , . The infinite sum in (2.8)_{1} involves all discrete wavevectors ** ξ**, whose components

*ξ*

_{i}satisfy the conditions 2.9Because

**e**

^{per},

*T*

^{per},

**q**

^{per}have zero averages in

*V*, the first term associated with

**=**

*ξ***0**in the Fourier series (2.8) vanishes, i.e 2.10Applying Fourier transform (2.8)

_{2}to equations (2.1), (2.3) and (2.7) results in the following relations between , , : 2.11where 1(

**) is the characteristic function such that 1(**

*ξ***0**)=1 and 1(

**≠**

*ξ***0**)=0. By some algebraic calculations, we obtain the expressions of and , which are valid for all

**≠**

*ξ***0**, 2.12The last equation of (2.12) can be recast into 2.13The spatial form of

**e**

^{per}can then be written as a Fourier series involving , 2.14or equivalently 2.15Finally, substituting (2.15) into (2.1) and (2.7) yields an integral equation of

**e***(

**x**) 2.16A dual integral equation can also be obtained based on the polarization

**q*** instead of

**e***. For thermal problems, the polarization field

**q*** is defined by 2.17where

**R**

^{0}is the reference resistivity tensor. By comparing (2.17) and (2.7), we must have 2.18

Applying Fourier transform to (2.17) and making use of (2.18) and (2.13), we obtain the following expressions that are valid for ** ξ**≠

**0**: 2.19

Next, we repeat the same steps as those used for deriving the **e***-based integral equation. First, we can compute **q**^{per}
2.20and then we obtain the dual integral equation for **q***
2.21where tensor **R**(**x**) in (2.21) is the local resistivity tensor.

### (b) Fast Fourier transform-based numerical methods

It is noted that tensor **S** in (2.13) is related to the periodic Green tensor associated with the reference material by
2.22and the infinite sum (2.14) is the Fourier series representation of the convolution involving the periodic Green operator ** Γ**,
2.23Combining (2.14), (2.23) and (2.6) yields an integral equation for

**e**, 2.24The dual formulation for flux

**q**can be deduced in a similar way from (2.21), 2.25with

**Δ**being the Green operator for fluxes. The representation of

**Δ**in Fourier space is 2.26Equations (2.24) and (2.25), and other related equations, are solved classically by iterative numerical schemes based on the FFT. Three of the main schemes with their features are summarized in the following.

— Primal iterative scheme (PIS) [9,14]: in the original method based on (2.24), the solution

**e**(**x**) is obtained by the recurrence process**e**^{i+1}=**E**−*(*Γ***K**(**x**)−**K**^{0})**e**^{i}starting from the initial value**e**^{1}=**E**. The PIS is suitable when the contrast ratio*k*_{I}/*k*_{M}is not too high.— Dual iterative scheme (DIS) [18]: the scheme is based on the dual form (2.25) and the recurrence process is applied to the flux

**q**, instead of**e**. The DIS can solve the problem related to the high ratio*k*_{I}/*k*_{M}, but diverges when the ratio tends to 0.— Polarization-based iterative scheme (PBIS): recently, Monchiet & Bonnet [10,19] proposed a PBIS based on a fictitious boundary-value problem where a uniform polarization field is applied to the periodic cell. The approach has proved to work at any

*k*_{I}/*k*_{M}ratio.

These iterative schemes will serve as numerical benchmarks for evaluating analytical solutions derived later in this paper.

### (c) Nemat-Nasser–Iwakuma–Hejazi approximation

The original approach of Nemat-Nasser *et al.* [13], abbreviated as NIH, was proposed for the effective elastic behaviour of composite materials. Numerical examples issued from Hoang [15] have shown that this approximation works well for composites with spherical inclusions of significant volume fractions *f*. Expecting similar results in the case of heat-conduction problems, we shall estimate the effective conductivity with the same procedure.

Taking the matrix as the reference material **K**^{0}=**K**^{M}, the effective conductivity can be computed from the average value of **e***. Indeed, averaging both sides of (2.7) over the cell volume *V* and making use of the fact that **e*** vanishes in the matrix, we obtain
2.27where the average is performed on the volume of inclusions, and *Ω* and *f* are, respectively, the volume occupied by the inclusion in the cell and its volume fraction. The effective conductivity defined by (2.4) is then calculated by
2.28The evaluation of 〈**e***(**x**)〉_{Ω} comes from the integral equation (2.16). Considering the matrix as the reference material and averaging both sides of (2.29) over *Ω* yields
2.29

The NIH approximation concerns the evaluation of 〈**e***(**x**) e^{−iξ⋅x}〉_{Ω} as
2.30

In the limit case where the inclusion is an ellipsoid and , this approximation is exact because **e*** is uniform inside the inclusion. This property is the consequence of the thermal Eshelby problem [20–25] concerning an ellipsoidal inclusion embedded in an infinite matrix. For finite *f*, the simplification (2.30) is expected to yield better prediction of **K**^{eff} than schemes based on Eshelby-type solutions. We define the two functions *I*(** ξ**) (the shape coefficient) and

*P*(

**) as 2.31Substituting (2.31) into (2.29) and accounting for (2.28) and (2.30) yields an explicit expression for the effective conductivity 2.32**

*ξ*A similar approximation can be applied to the dual formulation (2.21) with the matrix as the reference material. Instead of computing 〈**e***(**x**)〉_{Ω}, we compute 〈**q***(**x**)〉_{Ω} from (2.21),
2.33and find the effective resistivity, **R**^{eff}=(**K**^{eff})^{−1}, using the formula
2.34The final result for the effective resistivity is given by
2.35

## 3. Closed-form solutions for effective conductivity

### (a) Cubic cell with spherical inclusions

In this section, we exploit (2.32) and (2.35) by considering some particular structures: the periodic cell is a cube of dimension *a*, containing identical spherical inclusions of radius *R*. Both matrix and inclusions are made of isotropic thermal materials with conductivity *k*_{M} and *k*_{I}. To compute analytically the shape coefficient *I*(** ξ**), we shall make use of the following results:
3.1for a single sphere of volume

*V*

_{s}centred at

**x**

_{c}. As the matrix is thermally isotropic, the tensors

**S**(

**) and**

*ξ***L**(

**) defined by (2.13) and (2.19) admit simple forms, 3.2Let us denote by the set containing all wavevectors described by (2.9).**

*ξ*If *P*(** ξ**) and comply with the conditions

—

*P*() and are invariant under a reflection transformation with respect to any of the three planes*ξ**Oξ*_{i}*ξ*_{j}(*i*,*j*=1,2,3), 3.3—

*P*() and are invariant after a permutation between any two variables*ξ**ξ*_{i},*ξ*_{j}, 3.4

then the infinite sums and are reduced to the diagonal forms 3.5

The proof of the diagonal form is presented as follows. The first condition (3.3) implies that the non-diagonal terms vanish, 3.6whereas from the second condition (3.4), we can deduce that all the diagonal terms are identical, 3.7

As a result, the effective material is isotropic with the conductivity *k*^{eff} given by (2.32),
3.8or by (2.35),
3.9

The discrete vector set specified by (2.9) depends uniquely on the dimension of the periodic cell. As the latter is a cube, both conditions for are automatically satisfied. The validity of conditions for *P*(** ξ**) depends on the arrangement of the spheres inside the cell. However, for structures considered in later sections, both conditions for

*P*(

**) are also satisfied and allow us to obtain an expression for**

*ξ**k*

^{eff}in a closed form.

### (b) Cubic lattice structures

#### (i) Simple cubic

In SC structures, the inclusion sphere is located at the centre of the periodic cell (figure 1*a*). Owing to (3.1), both functions *I*(** ξ**) and

*P*(

**) are reduced to functions of its modulus**

*ξ**ξ*only, 3.10and thus verify both conditions (3.3) and (3.4). The remaining work is to evaluate in a closed form the infinite sum in (3.8).

First, we remark that the ** ξ**’s form a regular grid with spacing 2

*π*/

*a*as illustrated in figure 2. Taking a sufficiently large volume element d

*V*, the total number d

*N*of the grid points in d

*V*can be estimated by multiplying d

*V*with the grid density

*ρ*, 3.11

At large *ξ*, *P*(*ξ*) varies slowly, and we can calculate in the region d*V*
3.12or even in arbitrary remote domain *D* by the integral
3.13

At small *ξ*, *P*(*ξ*) fluctuates strongly, and the above approximation was found to be erroneous in previous tests. To obtain accurate approximations of the infinite sum, we keep several initial terms of the sum and estimate the remainder analytically with the integral. The infinite sum can be estimated by the formula
3.14in which *ξ*_{c} is the cut-off distance for ** ξ**, which defines the number of initial terms accounted for in this new expression. By changes of variable, the integral in (3.14) is rewritten as
3.15

The primitive integral is given by the analytical expression
3.16where *Si*(*η*) is the sine integral
3.17

We can now obtain the new expression for the infinite sum (3.14), 3.18

In most practical examples, keeping up to four initial terms can generate satisfactory results. For the SC lattice case, these four terms correspond to ** ξ** lying within a cut-off distance

*ξ*

_{c}=2

*ϵ*/

*R*(here,

*ϵ*=2

*πR*/

*a*). We can easily find that 3.19When substituting (3.19) back into (3.8) and (3.9), we obtain closed-form solutions for

*k*

^{eff}.

#### (ii) Body-centred cubic

In BCC structures, one sphere is located at the centre of the periodic cubic cell and eight others (each one counting for one eighth) at the eight corners (figure 1*b*). Making use of (3.1) and (2.31), function *P*(** ξ**) is given by the expression
3.20

This form of *P*(** ξ**) satisfies both conditions (3.3) and (3.4), and we continue to estimate the infinite sum of

*P*(

**) in (3.8) with the same procedure as the one established previously for the SC lattice. We remark that the term can take either one of values 0 or 4 with the same probability: it depends whether**

*ξ**n*

_{1}+

*n*

_{2}+

*n*

_{3}is even or odd. On average, we can estimate the sum of

*P*(

**) in domain d**

*ξ**V*, 3.21where is the average value 3.22

The infinite sum now becomes 3.23

Finally, by keeping four initial terms with ** ξ** lying within , the infinite sum can be computed with the following closed-form formula:
3.24

#### (iii) Face-centred cubic

In FCC lattice structures, eight spheres (each one counting for one eighth) are located at the corners and six others (each one counting for one half) at the centres of the cubic face (figure 1*c*). Applying the same procedure as the one described in §3*b*(ii)) and keeping four initial terms, we obtain
3.25

### (c) Randomly distributed spheres

#### (i) Relations to structure factor

In this subsection, we consider the case where a cube of size *a* contains a large number *N* of identical, non-overlapping spheres of radius *R* (figure 3). The shape functions *I*(** ξ**) and

*P*(

**) become 3.26where**

*ξ***x**

_{i}is the location of the sphere numbered

*i*in the cell. As long as the ergodic media hypothesis is valid, at the infinite-volume limit, we can replace the volume average 〈

**e***〉

_{Ω}(or 〈

**q***〉

_{Ω}) and the double sum in

*P*(

**) by their ensemble averages [12], and remark that 3.27**

*ξ*As suggested, the notation 〈〉_{ens} in (3.27) indicates the ensemble average and *S*(** ξ**) is the static structure factor, a statistical mechanics tool used for studying the local particle distribution [26–28]. The structure factor

*S*(

**) is related to the RDF**

*ξ**g*(

**r**) by 3.28where

*n*is the particle density. Regarding the particle density, it is related to the volume fraction

*f*by 3.29

When *g*(**r**) is isotropic and , *S*(** ξ**) is a function of

*ξ*only. By changes of variable while accounting for (3.29), (3.28) becomes 3.30

Finally, we can rewrite the new form for *P*(*ξ*) using the static structure factor *S*(*ξ*)
3.31

The infinite sum of *P*(*ξ*) can be evaluated either analytically with the procedure described in the previous section, i.e
3.32or computed numerically. Although equation (3.32) is obtained from the relation between *g*(*r*) and *S*(*ξ*), which is valid in the infinite-volume limit (i.e. ), it can be used to determine, as an approximation, the effective conductivity for finite-volume cases. Finally, for a given volume fraction *f*, *k*^{eff} depends on the cell size, via the parameter *R*/*a*.

#### (ii) Effective conductivity in the infinite-volume limit

To obtain the effective conductivity in this case, we compute the infinite sum of *P*(*ξ*) for a finite volume *V* and consider the following limit when ,
3.33

When the ratio *R*/*a* tends to 0, the ** ξ**’s grid becomes denser. For any given domain

*D*in

**space (figure 2), it is filled with an infinitely large number of grid points**

*ξ***with infinitesimally small spacing. Thus, equation (3.13) should hold for all**

*ξ**D*as long as , 3.34

The infinite sum of *P*(*ξ*) can now be evaluated by a continuous integral,
3.35or explicitly,
3.36

It can be noticed that for non-overlapping spheres, *g*(*r*)=0 when , and the following property holds:
3.37

Therefore, we can deduce that the infinite integral or the infinite sum of *P*(*ξ*) is given by
3.38

Substituting (3.38) back into (3.8) or (3.9), we obtain
3.39which is the Clausius–Mossotti equation (or Maxwell approximation formula, Hashin–Shtrikman bound) for conductivity [11]. It is interesting to remark that the NIH approximation in the homogenization limit gives results depending on *k*_{M},*k*_{I},*f*, but not on the local structure of the particles *g*(*r*),*S*(*ξ*). The accuracy of (3.39) for suspensions of spherical particles was discussed in Bonnecaze & Brady [16]. It was shown that (3.39) predicts very well the effective conductivities of these random composite materials when the contrast ratio *k*_{I}/*k*_{M} is not too high and the particles are more or less separated. In extreme situations, the interaction between the particles is much stronger and the NIH approximation can be erroneous. These effects can be viewed clearly from the percolation theory and should be studied by an adequate analysis [16,29].

To provide numerical proof to the infinite-volume solution (3.39), we carry out the convergence study of (3.33) as for two systems with different local structures (figure 4): the well-stirred distribution and the PY distribution for hard spheres.

#### (iii) Well-stirred distribution

The first case under consideration is a structureless system whose RDF function is expressed by the Heaviside function 3.40

This distribution is called the well-stirred distribution [30,31]. Combining (3.30) and (3.40) leads to the following analytical form of *S*(*ξ*) and *P*(*ξ*):
3.41
Markov & Willis [31] showed that the well-stirred distribution is only realistic for volume fraction . Otherwise, the positive-definiteness of the two-point correlation function is violated.

#### (iv) Percus–Yevick distribution for hard spheres

As shown previously, the well-stirred distribution is a simple/crude approximation and only valid for . We shall use now a more realistic and sophisticated distribution. In the framework of statistical mechanics, Percus & Yevick [32] proposed an integral equation to determine the pair correlation function of interacting particles. Wertheim [33] solved the PY equation for the case of hard spheres and Drugan & Willis [30] used Wertheim’s solution to determine the minimal size of a representative volume element. The following calculation is based on Wertheim’s solution presented by Drugan & Willis [30].

The solution of Wertheim [33] for *g*(*r*) is given as the Laplace transform of the RDF function,
3.42in which *L*(*t*) and *M*(*t*) are polynomials defined by
3.43

There is no explicit closed-form formula for *g*(*r*), but we can compute *S*(*ξ*) using (3.30) and (3.42), and the connection between the Laplace transform and the Fourier transform. Considering a special case where *s* is purely imaginary, say *s*=*iη*, we can prove that
3.44

Combining (3.42) and (3.44), we can obtain *S*(*ξ*) in closed form,
3.45

The variation of *S*(*ξ*) with respect to volume fraction *f* is shown in figure 5. The infinite sum of now becomes
3.46

## 4. Numerical applications

### (a) Cubic lattice structures

The closed-form solutions obtained previously are based on the approximations to the governing integral equations of **e*** (or **q***). In this section, these solutions will be compared with exact numerical solutions issued from the full FFT method using Neumann’s series. The inclusion/matrix conductivity ratio value *k*_{I}/*k*_{M} ranges from 10^{−3} to 10^{3}, and the volume fraction *f* from 0 to 90 per cent of the maximal value . Regarding numerical FFT iterative schemes, we adopt a resolution of 128×128×128 and control the difference between the consecutive steps at 10^{−3}. For the numerical evaluation of lattice sums in NIH approximations, a value *N*_{c}=128 is also used, i.e sum over all *n*_{i} such that |*n*_{i}|<*N*_{c}. Preliminary tests show that the given parameters generate satisfactory results with a reasonable computation time. The symmetry of the ** ξ**’s grid with respect to nine planes

*ξ*

_{i}=0 and

*ξ*

_{i}=±

*ξ*

_{j}(

*i*,

*j*=1,2,3) is also exploited to increase the computational efficiency.

From figures 6–9, we find an excellent agreement between NIH approximations and the closed-form expressions for all volume fractions and lattice structures at *k*_{I}/*k*_{M}=10. When compared with FFT-based solutions, figures 6–9 also showed a good agreement up to . Both FFT methods, PIS and PBIS, give the same results, while the deviation between the NIH-based and the FFT-based solutions is more noticeable at higher volume fractions. Despite the presence of different closed-form expressions, solutions based on **e*** and **q*** are very close to each other (figure 10).

Figure 10 shows a higher discrepancy for high *k*_{I}/*k*_{M} ratios (*k*_{I}/*k*_{M}=10^{3}) between the NIH-based solution and the FFT solution, whereas a very good agreement between them is observed for small *k*_{I}/*k*_{M}=10^{−3}.

### (b) Randomly distributed spherical inclusions

The ingredient of analytical solutions for these randomly heterogeneous systems *P*(*ξ*) in (3.41) and (3.46) depend on the ratio *R*/*a* and *f*. The former reflects the size-dependence effect, and we are interested in the convergence of the results in the infinite-volume limit, i.e. . In this paper, the following parameters are used: *k*_{I}/*k*_{M}=10, *f* ranging from 0 to 0.5, *R*/*a* from 0 to 0.3. For low *R*/*a* ratios, a higher resolution, *N*_{c}=1024, is required to ensure the convergence of the lattice sum.

Regarding the convergence at , from figures 11 and 12, we observe that the conductivity curves for both the PY and the well-stirred distribution become stable at *R*/*a*=0.1, and coincide with the analytical solution at *R*/*a*=0, the Hashin–Shtrikman bound (3.39).

## 5. Conclusions

In this paper, we have derived analytical expressions for effective conductivity of two-phase composites with spherical inclusions. Different arrangements of inclusions have been considered, including lattice structures and random distributions. Although the solutions are based on NIH approximations of solutions of integral equations, closed-form results for lattice structures agree very well with numerical FFT solutions for a large range of volume fractions. Regarding systems with randomly distributed spheres, the analytical solutions depend on the structure factor, and are size dependent. A convergence study at the infinite-volume limit was also carried out to find the effective conductivity.

Similar closed-form expressions can also be obtained in the case of anisotropic ellipsoidal inclusions. For random distribution cases, it will be of interest to compare with numerical or experiment results. These extensions will be the subjects of future works.

## Footnotes

↵1 Definitions of

**e*** and**q*** are given later in §2.

- Received June 4, 2012.
- Accepted December 11, 2012.

- © 2013 The Author(s) Published by the Royal Society. All rights reserved.