Computation of thermal properties via 3D homogenization of
multiphase materials using FFT-based accelerated scheme
Sophie Lemaitre, Vladimir Salnikov, Daniel Choı̈, and Philippe Karamian
Nicolas Oresme Mathematics Laboratory
University of Caen Lower Normandy

arXiv:1504.07499v1 [] 28 Apr 2015

CS 14032, Bd. Maréchal Juin, BP 5186
14032, Caen Cedex, France

Abstract
In this paper we study the thermal effective behaviour for 3D multiphase composite material
consisting of three isotropic phases which are the matrix, the inclusions and the coating media.
For this purpose we use an accelerated FFT-based scheme initially proposed in Eyre and Milton
(1999) to evaluate the thermal conductivity tensor. Matrix and spherical inclusions media are
polymers with similar properties whereas the coating medium is metallic hence better conducting.
Thus, the contrast between the coating and the others media is very large. For our study, we
use RVEs (Representative volume elements) generated by RSA (Random Sequential Adsorption)
method developed in our previous works, then, we compute effective thermal properties using
an FFT-based homogenization technique validated by comparison with the direct finite elements
method. We study the thermal behaviour of the 3D-multiphase composite material and we show
what features should be taken into account to make the computational approach efficient.
Keywords: composite material, coated medium, 3D-multiphases, stochastic homogenization, FFT

1

I.

INTRODUCTION AND MOTIVATIONS

In this paper, we study the thermal conductivity tensor using the 3D homogenization
method for materials constituted of 3 phases. Since composite materials are widely used
in modern industry, there has been a number of works studying their effective properties:
they concern descriptions of modelling approaches as well as results of experimental measurements. The purpose of this paper is to explore the range of application and features
for existing modelling approaches, and at the same time the work is motivated by precise
application within the framework of an industrial project.
The importance of modelling for the analysis of composite materials results from various
challenges in carrying out experimental works: such problems include purely technical difficulties as well as financial optimization of the applied activities. It is thus necessary to
master modelling approaches that are reliable, comparable with experiment, and still affordable from the point of view of computational efficiency. Our strategy in this paper is
related to the idea of stochastic homogenization. We consider samples of a composite material containing coated inclusions; and in order to take into account possible imperfections
and random factors, we average the result for a series of tests representing the same macro
characteristics. For this work we have chosen to study relatively simple geometry of inclusions: spherical shape with uniform coating. On the one hand this permits to concentrate
on the details of application of the adopted modelling technique, on the other hand this
can still provide some insight on the influence of morphology on the effective properties of
composite materials.
The paper is organized as follows. In the next section, we explain the RSA (Random Sequential Adsorption) method chosen for the generation of RVEs (Representative volume
elements) proposed in [1] and we give the context of the study. Next, we recall the model
used with the computational method. In the third section we compare the FFT (Fast Fourier
Transform)-based homogenization technique and the finite elements method (FEM) for simple geometries and we motivate our choice of the FFT-based ‘accelerated scheme’. The last
section describes the numerical results. We conclude by recapitulating the main points of
this work and describing the work in progress.

2

II.

SAMPLE GENERATION AND SETTING OF THE STUDY

In this section, we describe the method used to generate RVEs with three phases. As
we have outlined in the introduction, we consider a material with the matrix, the spherical
inclusions and the coating media.
We proceed in two steps. First, we generate RVEs using the tool developed by V. Salnikov
et al. in [1]. With this tool using the RSA algorithm generation method, we are able to
generate three phased RVEs containing up to 100 identical spherical inclusions and up to
40% in volume fraction. In the cited paper, it has been shown that the method is very
efficient being able to generate such RVEs in fractions of a second. The described algorithm
produces a list of inclusions in a ”vector” form namely for each sphere coordinates of the
center (in 3D) and the radius. The spheres do not overlap. In the second step, we use a
voxelization tool ([2]) to obtain a voxelized sample with three phases. We add a coating
defined by the user for each spherical inclusion. The figure 1 shows a scheme representing
a section of a sphere with the notations used throughout the article. In the sequel, we note
l the layer which is the ratio between the coating width and the radius of a sphere rs , the
layer is thus a number between 0 and 1. On the figures 2, we see one 3D RVE voxelization
and two of theses sections.

FIG. 1: Sphere made of a spherical inclusion with a coating parametrized by a layer l

The setting of the study is the following: we fix the parameters like the volume fraction of
spheres fsp equal to 30% (which takes into account spherical inclusion and coating) and the
phases conductivity tensor for each phase equals to kI where k is a positive number (it takes
3 different values for the 3 materials and I stands for the identity matrix in dimension 3.
In this setting, we study how the thermal properties of a composite vary depending on a
geometrical parameter denominated as the layer keeping in mind that the spherical volume
3

FIG. 2: Example of one RVE meshed with 3 phases generated by RSA method and two sections
in 1283 voxels representing the top face and one of these center sections

fraction is fixed here at 30%. For this purpose, we generated samples with different quantity
of spheres and different resolutions to observe the behaviour of the composite.

III.

COMPUTATIONAL TECHNIQUES

Consider a representative volume element V with periodic boundary conditions, and
denote θ(x) the temperature, φ(x) the heat flux and L(x) the thermal conductivity at any
point x ∈ V .
Fourier’s law in the linear case states that: φ(x) = −L(x) ∇θ(x). Moreover, we suppose that we are in stationary phase and without heat source so the heat flux verifies
div(φ(x)) = 0. Notice that for a composite material, the thermal conductivity tensor
depends on the point x: the dependence is governed by microscopic geometry of the sample,
namely which phase (matrix, inclusion or coating) the point belongs to.
The volume element V is subjected to an average temperature gradient h∇θ(x)i = ∇Θ
which induces local temperature gradient ∇θ(x) and heat flux φ(x) inside the RVE. The
effective constitutive relations of the composite are the relations between the (spatial)
average flux Φ and average temperature gradient ∇Θ.
It is a common practice to introduce a constant reference tensor L0 to solve the problem
of recovering the local fields and gradients. Introducing this homogeneous linear tensor, we

4

can rewrite the problem as:


0


 φ(x) = −L (x)∇θ(x) + τ (x), ∀x ∈ V
div(φ(x)) = 0, ∀x ∈ V



 τ (x) = −(L(x) − L0 )∇θ(x), ∀x ∈ V

(1)

where τ denotes the polarization tensor.
The solution of (1) can be expressed in real and Fourier spaces, where b
• denotes the Fourier
image of • and ξj ’s are the coordinates in the Fourier space, respectively, with the help of
the periodic Green operator Γ0 associated with L0 we can also write:
∇θ(x) = Γ0 (x)∇τ (x) ∀x ∈ V

(2)

c
c
b0 (ξ)b
∇θ(ξ)
=Γ
τ (ξ) ∀ξ 6= 0, ∇θ(0)
=0

(3)

or in Fourier space:

The equations (2) and (3) give the Lippmann-Schwinger integral equations in real and
b0 is easily expressed and computed in the
Fourier spaces respectively. The Green operator Γ
Fourier space by:
b0 (ξ) = P
Γ
ij

ξi ξj
L0mn ξ m ξ n

m,n

The Lippmann-Schwinger equation can be solved iteratively using the algorithm based on
the accelerated scheme proposed by Eyre and Milton in [3] and written as algorithm 3 in [4]
for the elastic case. For the thermal conductivity case, we obtain with adapted notations
the following algorithm:

5

Algorithm: FFT-based homogenization scheme for thermal conductivity case
Initialize ∇θ0 (x) ≡ ∇Θ, fix the convergence criterion acc.
while (not converged)
1. Convergence test:
n
n
n
bn
if ( comp
q < acc) compute φ (x) = −L(x)∇θ (x), φ = F F T (φ ),
b n (ξ)k2 i/kφ(0)k
b
eq = hkξ φ

if (eq < acc) then converged, stop
2. τ n (x) = −(L(x) + L0 )∇θ(x)
3. τb n = F F T (τ n )
c n (ξ) = −Γ
c n (0) = ∇Θ
b0 (ξ)b
4. ∇θ
τ n (ξ), ξ 6= 0, ∇θ
comp
comp
n

n
c
= F F T −1 (∇θ
5. ∇θcomp
comp )
q
n
k2 i/k∇Θk
6. comp = hk∇θn − ∇θcomp
n
(x) − ∇θn (x))
7. ∇θn+1 (x) = ∇θn (x) − 2(L(x) − L0 )−1 L0 (∇θcomp

When this algorithm converges, we compute hφ(x)i and we obtain Lhom using the following
equality: hφ(x)i = −Lhom h∇θ(x)i = −Lhom ∇Θ.
There are several FFT-based numerical schemes like the ’basic scheme’, the ’dual scheme’,
the augmented Lagrangian scheme, the ’polarization scheme’ and the ’accelerated scheme’.
We have chosen the last one for its computational efficiency like in the elastic case [4].
We justify this choice by relying on Moulinec and Silva in [5]; they have shown that the
accelerated scheme is the optimal compromise except for composite with infinite contrast.
In this paper, the contrast of the materials studied is between 5 and 2.103 .

IV.

COMPARISON FEM AND FFT METHOD

We briefly describe how we evaluate the homogenized thermal properties of the composite using the FEM. We use a double-scale method described by Sanchez-Palencia [6] and
Bensoussan, Lions and Papanicolaou [7]. It deals with setting a multi-scale problem via
an asymptotic expansion of the equations describing the behaviour of the material. In our

6

case, we consider a double-scale problem in the context of the thermal conductivity. The
composite material is associated to the macroscopic scale and the RVE to the microscopic
one. We impose periodic boundary conditions.
In order to compare both approaches, we study simple cases namely: one sphere or four
spheres centred as it is shown on figure 3.

FIG. 3: Two cuts representing RVE with one coated sphere and four coated spheres in resolution
643 voxels with sufficient coating

We generate RVEs in resolution 643 and 1283 voxels and with several thickness of the layer.
For each generated RVE, we construct a mesh based on the voxels then each voxel defines
a cub8 element and we use a python code to compute the thermal homogenized tensor.
The RVEs with a thin coating in low resolution namely 643 voxels or lower, exhibits some
voids on the coating as one can notice on figure 4. In these cases, the results obtained by
the two methods show some deviation which can reach 20%.

FIG. 4: A cut representing RVE with one sphere with thin coating

However, for the RVEs with sufficient coating such in figure 3, the results obtained by the
7

different methods are very similar, the discrepancy between both methods does not exceed
1%. An important remark is that with more complex geometries, the 3D homogenization
requires a lot of memory resources and takes a lot of CPU time with the FEM whereas the
FFT method has a lower cost in time and it allows the computation in 2563 voxels and even
in 5123 voxels (using a computation cluster). Thus we extend the FFT method to more
complex geometries namely with at least 30 spheres located randomly without intersection.

V.

NUMERICAL SIMULATIONS

Numerical computations are made for 3 resolutions: 643 , 1283 and 2563 voxels. According
to these resolutions, for the same geometries, the number of voxels representing each phase
increases. We also take 30 up to 100 spheres for nsp . We recall that the layer is a real
number between 0 and 1. We make the tests with l equal to 0.02 up to 0.08 and to 0.1 up
to 0.9. In fact, from a layer equal to 0.2, the coating reaches almost the half of the sphere
volume as we see in the table I. The coating volume fraction is evaluated by the equality :
fcoating = (1 − (1 − l)3 ) × 0.3, where 0.3 denotes the spheres volume fraction.

For each set of (Res, nsp , l), we generate 10 different samples in order to make the average
of obtained macroscopic conductivity tensor.
To compute the thermal conductivity tensor, we recall that the materials are supposed to
be isotropic namely with a tensor equal to kI, where k is a positive number and I stands
for the identity matrix in dimension 3. We take:
Lmatrix phase = I , Linclusion phase = 0.2I , Lcoating phase = 400I
With theses values, the higher contrast is 2.103 and the lower is 5. We deliberately choose
these coefficients for industrial applications with matrix and inclusion phases poorly conducting and coating phase highly conducting. According to [3] and [8] the constant reference
thermal tensor L0 is set to:
r
√
min(x) × max(x) I = −4 5 I
L =−
0

x∈{1,0.2,400}

x∈{1,0.2,400}

8

layer fcoating
0.02 0.0176
0.04 0.0346
0.06 0.0508
0.08 0.0664
0.1 0.0813
0.2 0.1464
0.3 0.1971
0.4 0.2352
0.5 0.2625
0.6 0.2808
0.7 0.2919
0.8 0.2976
0.9 0.2997
TABLE I: Coating volume fraction depending on the thickness of the layer

We remark that the computed thermal homogenized tensor is close to an isotropic material tensor, so we take as the homogenized thermal conductivity the mean of the diagonal
elements of the tensor for the series of the 10 geometries generated as discussed above.
The curves below on figures 8, 9, 10 give these homogenized values. First, we naturally
observe that the homogenized thermal conductivity increases with the thickness of the layer
and for all the three resolutions tested. We can also note that the larger the resolution is,
the less the dispersion is. We conclude that resolution 64 is not accurate enough.
We observe for very thin layers that the homogenized thermal conductivity decreases when
the nsp increases because the coating is thin and the number of voxels representing it decrease
and can be less than one pixel, i.e. it can be missed: this is shown in the tables II, III, IV
where we use the equality:
Number of voxels = l × Res ×



3×0.3×nsp
4×π

1/3

, with Res ∈ {64, 128, 256}.

The values in the tables are deliberately not rounded up.

9

FIG. 5: Resolution 643 : homogenized thermal conductivity evolution vs layer for nsp from 30 up
to 100

nsp

30

40

50

60

70

80

90 100

l
0.02

0.17 0.16 0.14 0.14 0.13 0.12 0.12 0.11

0.04

0.34 0.31 0.29 0.27 0.26 0.25 0.24 0.23

0.06

0.51 0.47 0.43 0.41 0.39 0.37 0.36 0.34

0.08

0.68 0.62 0.58 0.54 0.52 0.49 0.47 0.46

0.1

0.86 0.78 0.72 0.68 0.64 0.62 0.59 0.57

0.2

1.71 1.55 1.44 1.36 1.29 1.23 1.19 1.15

0.3

2.57 2.33 2.16 2.04 1.93 1.85 1.78 1.72

0.4

3.42 3.11 2.89 2.72 2.58 2.47 2.37 2.29

TABLE II: Voxels coating thickness for the resolution 643 voxels

10

FIG. 6: Resolution 1283 : homogenized thermal conductivity evolution vs layer for nsp from 30 up
to 100

We also note for resolution 2563 voxels that the number of spheres almost does not influence
the computed values, which agrees with the definition of acceptable RVE size and means we
no longer have the ’artefact’ due to the voxelization. Moreover, for a layer l equal to 0.08,
there is a kind of stagnation, it is shown on figures 8, 9, 10. On these figures, we clearly
observe two behaviours, one before the stagnation and the other one during the stagnation.
This effect is less visible in resolution 128 since it is difficult to capture the coating when
the layer is too low which appears in figures 11. From the point of view of applications we
11

FIG. 7: Resolution 2563 : homogenized thermal conductivity evolution vs layer for nsp from 30 up
to 100

nsp

30

40

50

60

70

80

90 100

l
0.02

0.34 0.31 0.29 0.27 0.26 0.25 0.24 0.23

0.04

0.68 0.62 0.58 0.54 0.52 0.49 0.47 0.46

0.06

1.03 0.93 0.87 0.81 0.77 0.74 0.71 0.69

0.08

1.37 1.24 1.15 1.09 1.03 0.99 0.95 0.92

0.1

1.71 1.55 1.44 1.36 1.29 1.23 1.19 1.15

0.2

3.42 3.11 2.89 2.72 2.58 2.47 2.37 2.29

0.3

5.13 4.66 4.33 4.07 3.87 3.70 3.56 3.44

TABLE III: Voxels coating thickness for the resolution 1283 voxels

12

nsp

30

40

50

60

70

80

90 100

l
0.02

0.68 0.62 0.58 0.54 0.52 0.49 0.47 0.46

0.04

1.37 1.24 1.15 1.09 1.03 0.99 0.95 0.92

0.06

2.05 1.87 1.73 1.63 1.55 1.48 1.42 1.37

0.08

2.74 2.49 2.31 2.17 2.06 1.97 1.90 1.83

0.1

3.42 3.11 2.89 2.72 2.58 2.47 2.37 2.29

TABLE IV: Voxels coating thickness for the resolution 2563 voxels

can say that it is not necessary to have a big coating to increase the homogenized thermal
conductivity. There is a limit due to the sphere volume fraction and due to the type of
spherical inclusions.

FIG. 8: Homogenized thermal conductivity evolution vs coating volume fraction for 30 spheres

VI.

CONCLUSIONS AND OUTLOOK

To conclude, we can summarize the main ideas of this paper. We have shown that
the FFT-based iterative accelerated scheme is a good tool for computing thermal effective
13

FIG. 9: Homogenized thermal conductivity evolution vs coating volume fraction for 70 spheres

FIG. 10: Homogenized thermal conductivity evolution vs coating volume fraction for 100 spheres

properties of a composite generated randomly in 3D. We have restricted our study to RVE
composed with only spheres but it can be made with cylinders or both spheres and cylinders.
In [1], RVEs generation algorithms in MD (Molecular dynamics) or RSA methods have

14

FIG. 11: Example of sections with lack of coating when the layer is too low

been presented using both spheres and cylinders and in [4] these RVEs have been used for
calculations in elastic case.
We have observed that the resolution 2563 with 50 − 70 spheres and an sufficient thickness
of the layer gives good results. We have also noted that the coating must be thick enough
in order to be ”captured”. In our work in progress, we explore the case with a thin coating,
where the methods are more subtle and some preparation of the samples is needed. We
will address this issue as well as more advanced morphological parameters in a separate
publication.

Acknowledgements

Most of the computations described in this paper have been carried out at the cluster of
the Centre of Informatics Resources of Higher Normandy (CRIHAN - Centre de Ressources
Informatiques de Haute-Normandie).
This work has been supported by ACCEA project selected by the ”Fonds Unique Interministriel (FUI) 15 (18/03/2013)” program.

15

[1] V. Salnikov, D. Choi, Ph. Karamian-Surville On efficient and reliable stochastic generation
of RVEs for analysis of composites within the framework of homogenization, Computational
Mechanics, Volume 55, Issue 1, 2015. Preprint: http://arxiv.org/abs/1408.6074, 2014
[2] S. Lemaitre, V. Salnikov, D. Choi, P. Karamian Approche par la dynamique molculaire pour
la conception de VER 3D et variations autour de la pixellisation, accepted for publication in
the proceedings of the CSMA 2015 https://csma2015.csma.fr/)
[3] D.J. Eyre, G.W. Milton A fast numerical scheme for computing the response of composites
using grid refinement, Journal of Physique III 1999; 6; 41-47
[4] V. Salnikov, S. Lemaitre, D. Choi, Ph. Karamian-Surville Measure of combined effects of
morphological parameters of inclusions within composite materials via stochastic homogenization to determine effective mechanical properties, To appear in Composite Structures.
DOI:10.1016/j.compstruct.2015.03.076, 2015, Preprint: http://arxiv.org/abs/1411.4037
[5] H. Moulinec, F. Silva Comparison of the three accelerated FFT-based schemes for computing
the mechanical response of composite materials, Int. J. Numer. Meth. Engng. 2014; 97:960-985
[6] E. Sanchez-Palencia Non-homogeneous media and vibration theory, vol. 127, Berlin, 1980)
[7] A. Bensoussan, J.L. Lions, G.C. Papanicolaou Asymptotic Analysis for Periodic Structures,
Springer Verlag, North Hollland, Amsterdam Edition, 1978)
[8] J.C. Michel,H. Moulinec, P. Suquet A computational scheme for linar and non-linear composites with arbitrary phase contrast , Int. J. Numer. Meth. Engng. 2001; 52:139-160(DOI:
10.1002/nme.275)

16

