D.-H Nguyen

H Le Quang
email: hung.lequang@univ-paris-est.fr

Q.-C He

A.-T Tran

Generalized Hill-Mendel lemma and equivalent inclusion method for determining the effective thermal conductivity of composites with imperfect interfaces

Keywords: Micromechanics, Composites, Thermal conductivity, Imperfect interfaces, Equivalent inclusion method

The present work aims at determining the effective thermal conductivity of twoor three-dimensional composites with imperfect interfaces between their constituent phases. These imperfect interfaces are described by the highly conducting, lowly conducting or general thermal imperfect model. To achieve the objective, the classical Hill-Mendel lemma is first extended to include the effects of imperfect interfaces and an equivalent inclusion method (EIM) is proposed. The basic idea of EIM is to replace an inclusion embedded in a matrix via an imperfect interface by an equivalent inclusion inserted in the same matrix via a perfect interface. Using EIM and applying the dilute distribution, Mori-Tanaka, self-consistent, generalized self-consistent and differential schemes, the effective thermal conductivities of layered composites and some particle-reinforced composites with imperfect interfaces are analytically and explicitly determined. These results are compared with the Voigt, Reuss and Hashin-Shtrikman bounds and checked against the numerical results provided by the fast Fourier transform (FFT) method. These comparisons and checks show that the methods proposed in this work are particularly efficient. The methods and results of the present work are directly transposable to other transport phenomena and anti-plane elasticity by their strict mathematical analogy with thermal conduction.

Introduction

Most of the classical microstructural models used in the micromechanics of composite materials to estimate their effective properties often adopt the hypothesis that the interfaces between the constituent phases are perfect. Within the context of thermal conduction, an interface is considered as perfect if and only if both temperature and normal heat flux are continuous across it. However, in many situations of practice, such as those in the presence of roughness, defect, damage or mismatch between the phases, the hypothesis of the perfect interface is no longer appropriate. Consequently, consideration of imperfect interfaces between the constituent phases of composites becomes indispensable for determining the effective properties of composite materials. Among all linear thermal imperfect interface models, the most widely used ones are the following three. First, the lowly contacting (LC) interface model, also called Kapitza's thermal resistance interface model, assumes that the normal component of the heat flux vector is continuous across an interface while the temperature across the interface presents a jump proportional to the normal heat flux component via the so-called Kapitza thermal resistance coefficient (see, e.g. [START_REF] Kapitza | The study of heat transfer in helium ii[END_REF]). The effect of Kapitza's thermal resistance on the effective thermal conductivity of composites has been studied in many works, i.e. [START_REF] Benveniste | A general interface model for a three-dimensional curved thin anisotropic interphase between two anisotropic media[END_REF][START_REF] Benveniste | The effective conductivity of composites with imperfect thermal contact at constituent interfaces[END_REF][START_REF] Benveniste | Effective thermal conductivity of composites with a thermal contact resistance between the constituents: Nondilute case[END_REF][START_REF] Hasselman | Effective thermal conductivity of composites with interfacial thermal barrier resistance[END_REF][START_REF] Lipton | Composites with imperfect interface[END_REF][START_REF] Lipton | Critical radius, size effects and inverse problems for composites with imperfect interface[END_REF][START_REF] Cheng | Effective conductivity of dispersions of spheres with a superconducting interface[END_REF][START_REF] Nan | Effective thermal conductivity of particulate composites with interfacial thermal resistance[END_REF][START_REF] Torquato | Effect of the interface on the properties of composite media[END_REF][START_REF] Le-Quang | Eshelby's tensor fields and effective conductivity of composites made of anisotropic phases with kapitza's interface thermal resistance[END_REF][START_REF] Le-Quang | Estimations of the effective conductivity of anisotropic multiphase composites with imperfect interfaces[END_REF][START_REF] Kushch | Effective conductivity of spheroidal particle composite with imperfect interfaces: Complete solutions for periodic and random micro structures[END_REF]. Second, viewed as being dual to the LC interface model, the highly contacting (HC) interface model or the coherent imperfect interface model, stipulates that the temperature field is continuous while the normal heat flux is discontinuous across an interface (see, e.g., [START_REF] Torquato | Effect of the interface on the properties of composite media[END_REF][START_REF] Lipton | Reciprocal relations, bounds, and size effects for composites with highly conducting interface[END_REF][START_REF] Lipton | Variational methods, bounds, and size effects for composites with highly conducting interface[END_REF][START_REF] Miloh | On the effective conductivity of composites with ellipsoidal inhomogeneities and highly conducting interfaces[END_REF][START_REF] Yvonnet | Numerical modelling of the effective conductivities of composites with arbitrarily shaped inclusions and highly conducting interface[END_REF][START_REF] Le-Quang | Size-dependent eshelby tensor fields and effective conductivity of composites made of anisotropic phases with highly conducting imperfect interfaces[END_REF][START_REF]Effective thermal conductivity of periodic composites with highly conducting imperfect interfaces[END_REF][START_REF] Le-Quang | Estimations of the effective conductivity of anisotropic multiphase composites with imperfect interfaces[END_REF]). According to the HC imperfect interface model, a surface intensity and a surface heat flux field defined on an interface are related by the surface Fourier law and satisfy the surface energy conservation equation. Third, according to the general imperfect (GI) model, proposed by Gu and his co-authors [START_REF] Gu | Coordinate-free derivation and weak formulation of a general imperfect interface model for thermal conduction in composites[END_REF][START_REF] Gu | Closed-form estimates for the effective conductivity of isotropic composites with spherical particles and general imperfect interfaces[END_REF], both the temperature and the normal heat flux are discontinuous across an interface and must satisfy some jump relations.

The general imperfect interface model was first proposed on the basis of some phenomenological arguments and next derived rigorously by applying asymptotic analysis. By considering a material interface as the limiting case of a very thin interphase between two bulk phases, it can be shown by asymptotic analysis that the LC model or the HC interface model can be derived from the general thermal imperfect model when the interphase is much less or more conducting than each of the constituents. This means that the general imperfect interface model includes the LC and HC models as particular cases. The asymptotic analysis applied to thermal conduction is closely related to some mathematical techniques of homogenization and was used in many works about interfaces [START_REF] Hashin | Thin interphase/imperfect interface in conduction[END_REF][START_REF] Benveniste | A general interface model for a three-dimensional curved thin anisotropic interphase between two anisotropic media[END_REF][START_REF] Gu | Interfacial discontinuity relations for coupled multifield phenomena and their application to the modeling of thin interphases as imperfect interfaces[END_REF][START_REF] Benveniste | The effective conductivity of composites with imperfect thermal contact at constituent interfaces[END_REF][START_REF] Miloh | On the effective conductivity of composites with ellipsoidal inhomogeneities and highly conducting interfaces[END_REF][START_REF] Sanchez-Palencia | Comportement limite d'un problème de transmission à travers une plaque mince et fablemenrt conductrice[END_REF][START_REF] Huy | Phénomènes de transmission a travers des couches minces de conductivité élevé e[END_REF][START_REF] Bövik | On the modelling of thin interface layers in elastic and acoustic scattering problems[END_REF].

The Hill-Mandel lemma plays a key role in analytical and numerical determinations of the effective properties of heterogeneous materials (see e.g. [START_REF] Hill | Elastic properties of reinforced solids: Some theoretical principles[END_REF][START_REF] Nemat-Nasser | Micromechanics: Overall Prosperities of Heterogeneous Materials[END_REF][START_REF] Qu | Fundamentals of Micromechanics of Solids[END_REF]). The Hill-Mandel lemma establishes the connection of the physical and mechanical energies at the microscopic scale with the ones at the macroscopic scale. This connection holds independently of any constitutive laws. The importance of the Hill-Mandel lemma to the mechanics and physics of heterogeneous materials is for two reasons. First, it allows us to check if the definitions of macroscopic quantities in terms of relevant microscopic ones comply with the micro-macro energy equivalence. Second, it allows us to know which are the boundary conditions compatible with the micro-macro energy equivalence. However, the classical Hill-Mandel lemma is established under the assumption that the interfaces between the phases of composites are perfect. In a variety of situations of theoretical or/and practical interest, this assumption is too idealized and a generalization of the classical Hill-Mandel lemma to the case where imperfect interfaces are involved becomes indispensable.

The present work is concerned with the determination of the effective thermal conductivity of two-or three-dimensional multiphase composites with imperfect interfaces. The aforementioned three imperfect interface models are used to describe and capture the effects of imperfect interfaces on the effective thermal conductivity of composites. By mathematical analogy, the methods elaborated and the results obtained in the present work on thermal conduction are directly transposable to other transport phenomena like electric conduction, dielectrics, magnetism, diffusion and flow in porous media. In addition, we can also exploit the fact that anti-plane elasticity is mathematically identical to two-dimensional thermal conduction. The objective of the present work is three-fold:

• First, it is to extend the classical Hill-Mandel lemma to the case where the interfaces between the constituent phases of composite materials are imperfect. In the generalized Hill-Mandel lemma, the total thermal energy at microscopic scale has to contain not only the thermal energy in the bulk phases but also the one of the imperfect interfaces. One of the important consequences of the generalized Hill-Mandel lemma is that when uniform intensity, uniform heat flux or periodic boundary conditions are prescribed, the micro-macro energy equivalence condition is ensured. • Second, it is to propose a novel approach, called equivalent inclusion method, in which each inclusion embedded in a host matrix phase via an imperfect interface is first replaced by an equivalent inclusion whose link with the matrix phase is perfect. The thermal conductivity of the equivalent inclusion can be determined with the aid of the generalized Hill-Mandel lemma.

Next, due to the fact that the interface between any inclusion and the ma-trix is perfect, classical schemes like the dilute distribution, Mori-Tanaka, self-consistent, generalized self-consistent and differential approximation models can be directly used to determine the effective thermal conductivity of composites. • Finally, it is to apply the equivalent inclusion method proposed to calculate the effective thermal conductivity of composites having different microstructure. The first application concerns a layered composite in which the interface between any two layers is imperfect. The second one is related to two-and three-dimensional composites consisting of circular or spherical inclusions embedded in a host matrix phase. Invoking the LC, HC and GI interface models, using the equivalent inclusion method, and applying the dilute distribution, Mori-Tanaka, self-consistent, generalized self-consistent and differential schemes, the effective thermal conductivities of layered composites and circular/spherical particle-reinforced composites are analytically and explicitly obtained. These results are compared with, and validated by, the corresponding numerical ones provided by the fast Fourier transform (FFT) method as well as the Voigt, Reuss and Hashin-Shtrikman bounds.

The paper is structured as follows. Section 2 is dedicated to specifying the phase constitutive laws of composites under investigation, the different thermal imperfect interface models and the general form of the effective thermal conduction behavior. In Section 3, the generalized Hill-Mendel lemma is given and proved. Section 4 presents the equivalent inclusion method and shows how to determine the thermal conductivity of composites with imperfect interfaces by using the generalized Hill-Mendel's lemma and the micro-macro energy equivalence. In Section 5, the effective thermal conductivities of layered composites and circular/spherical particle-reinforced composites with imperfect interfaces are analytically and explicitly determined by applying the equivalent inclusion method and some micromechanical schemes. The effects of imperfect interfaces on the effective conductivities of composites are discussed; comparisons with the corresponding numerical results provided by the FFT method and with the Voigt, Reuss and Hashin-Shtrikman bounds are made.

In section 6, a few concluding remarks are provided.

Problem setting

In a d-dimensional space R d with d = 2 or 3, the composite material under consideration consists of a matrix, called phase 0, in which inclusions, referred to as phase 1, are embedded. Both the matrix and inclusion phases are assumed to be individually homogeneous. Let Ω be the domain occupied by a representative volume element (RVE) of the two-phase composite investigated and let Ω (p) correspond to the subdomain occupied by phase p with p = 0 or 1. In what follows, we denote by ∂Ω the boundary of Ω and by ∂Ω (p) the boundary of Ω (p) .

Relative to the system of d-dimensional Cartesian coordinates {x 1 , . . . , x d } associated to a right-handed orthonormal basis {f 1 , . . . , f d }, the RVE Ω is subjected to general mixed boundary conditions as follows:

θ(x) = θ 0 (x), ∀x ∈ ∂Ω t , (1) q
(x) • m(x) = Q 0 (x), ∀x ∈ ∂Ω q (2)
where θ 0 and Q 0 are, respectively, the prescribed temperature and heat flux on the complementary parts, ∂Ω t and ∂Ω q , of ∂Ω; m is an outward unit vector normal to ∂Ω. For later use, three following particular boundary conditions are introduced:

• Uniform intensity boundary condition:

θ(x) = -E 0 • x, ∀x ∈ ∂Ω; (3)
• Uniform heat flux boundary condition:

q(x) • m(x) = Q 0 (x) • m(x), ∀x ∈ ∂Ω; (4)
• Periodic boundary condition:

θ(x) = -E 0 • x + θ(x), with θ(x) periodic on ∂Ω, q(x) • m(x) anti-periodic on ∂Ω. (5)
Here, E 0 and Q 0 are prescribed constant intensity and heat flux vectors.

In this work, the interface between phase 0 and phase 1, designated by Γ and defined by Γ = ∂Ω (0) ∩ ∂Ω (1) , is assumed to be imperfect in the sense that the temperature field or/and the normal component of heat flux field is/are discontinuous across it. More precisely, the following three imperfect interface models will be considered:

• Lowly conducting (LC) interface or Kapitza's interface thermal resistance model:

q(x) • n(x) = 0, θ(x) = -αq (1) (x) • n(x) = -αq (0) (x) • n(x), ∀x ∈ Γ; (6
) • Highly conducting (HC) interface or coherent imperfect interface model:

θ(x) = 0, q(x) • n(x) = -β∆ s θ (1) (x) = -β∆ s θ (0) (x), ∀x ∈ Γ; (7) • General imperfect (GI) interface model θ(x) = f q (1) (x) • n(x), q (0) (x) • n(x) , q(x) • n(x) = g θ (1) (x), θ (0) (x) , ∀x ∈ Γ. (8)
In these equations and hereafter, n designates a unit outward normal vector of Γ directed from Ω (1) to Ω (0) ; the symbol • = • (0) -• (1) presents the jump of • across Γ; ∆ s • denotes the surface Laplacian of •; α stands for the Kapitza thermal resistance while β is the surface thermal conductivity of Γ; f q (1) (x) • n(x), q (0) (x) • n(x) and g θ (1) (x), θ (0) (x) are two scalar functions whose expressions are omitted here. However, in the case where the imperfect interface is considered as replacing a thin interphase and the thermal constitutive law of this interphase is assumed to be linear, these two scalar functions f q (1) (x) • n(x), q (0) (x) • n(x) and g θ (1) (x), θ (0) (x) can be analytically and explicitly established by applying the Taylor's expansion.

For more details about these functions as well as their derivation, the reader can refer to the papers of Gu et al. [START_REF] Gu | Coordinate-free derivation and weak formulation of a general imperfect interface model for thermal conduction in composites[END_REF] and Le Quang [START_REF] Le-Quang | Estimations and bounds of the effective conductivity of composites with anisotropic inclusions and general imperfect interfaces[END_REF].

By definition, a temperature field θ(x) over Ω is said to be kinematically admissible if and only if it is continuous over subdomains Ω (1) and Ω (0) but not necessary to be continuous across the interface Γ between Ω (1) and Ω (0) , continuously differentiable over each subdomain Ω (p) and verifies the boundary condition (1) on ∂Ω t . Then, an almost everywhere continuous intensity field e(x) over Ω is qualified as kinematically admissible if and only if it is derived from a kinematically admissible temperature field θ(x) by e(x) = -∇θ(x).

On the other hand, a heat flux field q(x) over Ω is defined to be statically admissible if and only if it simultaneously satisfies the boundary condition [START_REF] Benveniste | A general interface model for a three-dimensional curved thin anisotropic interphase between two anisotropic media[END_REF] and the following energy conservation equation

∇ • q(x) = 0 (10)
in the case of steady thermal conduction and in the absence of heat source.

Next, we define the macroscopic intensity vector E and the macroscopic heat flux vector Q as follows:

E = - 1 |Ω| ∂Ω θ(x)m(x)dx, Q = 1 |Ω| ∂Ω [q(x) • m(x)]xdx (11)
where |Ω| denotes the volume or area of the domain Ω according as the tridimensional (3D) or bi-dimensional (2D) case is concerned.

It is important to notice that, when the normal heat flux field component and the temperature field are discontinuous across the matrix/inclusion interface Γ, the macroscopic heat flux field Q and the macroscopic intensity field E defined by [START_REF] Le-Quang | Eshelby's tensor fields and effective conductivity of composites made of anisotropic phases with kapitza's interface thermal resistance[END_REF] are not simply the volume averages of the local counterparts q(x) and e(x) as in the classical case with perfect interface. More precisely, from the definition [START_REF] Le-Quang | Eshelby's tensor fields and effective conductivity of composites made of anisotropic phases with kapitza's interface thermal resistance[END_REF], it can be shown that the macroscopic heat flux field Q is given by

Q = q + 1 |Ω| Γ q(x) • n(x)xdx, (12)
and the macroscopic intensity field E has the expression

E = e - 1 |Ω| Γ θ(x) n(x)dx. ( 13
)
Above • denotes the volume or area average of • defined over the domain Ω as

• = 1 |Ω| Ω (1)
• (1) dx +

Ω (0) • (0) dx . (14)
Finally, even in presence of imperfect interface, it can be shown with the help of the divergence theorem that, as in the classical case with perfect interface,

E = E 0 (15)
when the uniform intensity boundary condition (3) or of periodic boundary condition ( 5) is concerned, and

Q = Q 0 (16)
when the uniform heat flux boundary condition ( 3) is under consideration.

3 Generalized Hill-Mendel's lemma Lemma 3.1 Given the boundary conditions (1) and (2), for any statically admissible heat flux field q(x) and for any kinematically admissible temperature field θ(x) together with the associated kinematically admissible intensity field e(x), the following equation holds

q(x) • e(x) + 1 |Ω| Γ -θ(x)q(x) • n(x)dx = E • Q + 1 |Ω| ∂Ω {q(x) -Q} • m(x) {-θ(x) -E • x} dx ( 17
)
where E and Q are defined by Eq. [START_REF] Le-Quang | Eshelby's tensor fields and effective conductivity of composites made of anisotropic phases with kapitza's interface thermal resistance[END_REF].

Proof. Firstly, by introducing the fluctuation parts of the temperature, intensity and heat flux fields defined respectively by θ(x) = θ(x) + E • x, ẽ(x) = e(x) -E and q(x) = q(x) -Q, the volume or area average over Ω of the thermal energy q(x) • e(x) can be expressed in the following equivalent form

q(x) • e(x) = (q(x) + Q) • (ẽ(x) + E) = Q • E + Q • ẽ(x) + q(x) • E + q(x) • ẽ(x) . (18)
Next, by applying the divergence theorem and by using Eqs. [START_REF] Torquato | Effect of the interface on the properties of composite media[END_REF] and [START_REF] Le-Quang | Eshelby's tensor fields and effective conductivity of composites made of anisotropic phases with kapitza's interface thermal resistance[END_REF], it can be shown that

ẽ(x) = -∇ θ = -E -∇θ = 1 |Ω| Γ θ ndx, (19) q
(x) = q -Q = - 1 |Ω| Γ q • nxdx, (20) q
(x) • ẽ(x) = -(q -Q) • ∇ θ = - 1 |Ω| ∂Ω {q(x) -Q} • m(x) {θ(x) + E • x} dx - Q i |Ω| Γ θ n i dx + E i |Ω| Γ q • nx i dx + 1 |Ω| Γ θq • ndx. (21)
Finally, substitution of Eqs. ( 19)-( 21) into the right-hand side of Eq. ( 18) leads therefore to [START_REF] Yvonnet | Numerical modelling of the effective conductivities of composites with arbitrarily shaped inclusions and highly conducting interface[END_REF].

It is important to notice from Eq. ( 17) that:

• Firstly, the generalized Hill-Mendel's lemma is valid independently of the constitutive laws of the materials constituting the matrix and inclusion phases; • Secondly, if the interface Γ between matrix and inclusion phases is perfect, then θ(x)q(x) • n(x) = 0 for any x ∈ Γ and the Generalized Hill-Mendel's lemma reduces to the classical Hill-Mendel's lemma; • Thirdly, when a uniform intensity, uniform heat flux or periodic boundary condition as described in Eq. ( 3), (4) or ( 5) is prescribed on ∂Ω, it can be shown that the following condition, called also Hill-Mendel's condition,

1 |Ω| ∂Ω {q(x) -Q} • m(x) {θ(x) + E • x} dx = 0 ( 22
)
holds and the Generalized Hill-Mendel's lemma takes the following simple form q(x)

• e(x) + 1 |Ω| Γ -θ(x)q(x) • n(x)dx = E • Q. ( 23
)
Physically, with the help of the definitions [START_REF] Le-Quang | Eshelby's tensor fields and effective conductivity of composites made of anisotropic phases with kapitza's interface thermal resistance[END_REF], Eq. ( 23) guarantees the equivalence between the thermal energy at the macroscopic scale and the counterpart at the microscopic scale which is composed of the thermal energy in each phase and the interfacial energy of Γ.

Equivalent inclusion method

As mentioned in Section 2, the representative volume element (RVE) Ω of the composite under consideration is made of the matrix phase Ω (0) in which the inclusion phase Ω (1) is embedded via the imperfect interface Γ. The materials forming the matrix and inclusion phases are assumed to comply with the Fourier law:

q (p) (x) = -K (p) (x) • e (p) (x), x ∈ Ω (p) , (24)
where e (p) (x) = -∇θ (p) (x) is the local intensity field and K (p) (x) stands for the thermal conductivity tensor of the material constituting the inclusion phase with p = 1 or the matrix phase with p = 0. Recall that the imperfect interface Γ between the matrix and inclusion phases, Ω (0) and Ω (1) , can be described either by HC, LC or GI interface model. In particular, for GI interface model, owing to the linear thermal conduction behavior of material constituents, the temperature and normal heat flux jump relations can be expressed by

θ(x) = t 2 1 k (1) nn - 1 k (c-) nn q (1) n (x) + 1 k (0) nn - 1 k (c+) nn q (0) n (x) + t 2 s (1) -s (c-) ∇ s θ (1) (x) + s (0) -s (c+) ∇ s θ (0) (x) , (25)
q n (x) = t 2 ∇ s • S (c-) -S (1) • ∇ s θ (1) (x) + S (c+) -S (0) • ∇ s θ (0) (x) + t 2 ∇ s • s (c-) -s (1) q (1) n (x) + s (c+) -s (0) q (0) n (x) (26)
with x ∈ Γ. In these equations, the surface gradient and divergence operators for a quantity •, denoted by ∇ s (•) and ∇ s • (•), are defined as

∇ s (•) = ∇(•) • T, ∇ s • (•) = ∇(•) : T ( 27
)
where T = I (d) -n ⊗ n with n being the unit vector normal to Γ oriented from the inclusion phase to the matrix phase, I (d) denotes the d-dimensional second-order identity tensor and

k ( * ) nn = n•K ( * ) •n, s ( * ) = n • K ( * ) k ( * ) nn , S ( * ) = K ( * ) - (K ( * ) • n) ⊗ (K ( * ) • n) k ( * ) nn (28)
where * = 1, 0 or c± represents the quantity relative to the inclusion Ω (1) , matrix Ω (0) or interphase between the matrix and inclusion phases of thickness t associated with the inclusion and matrix side, respectively. In the following, we consider three particular and important cases where Γ is flat, circular and spherical and the matrix, inclusion and interphase are orthotropic, circularly or spherically transverse isotropic, respectively. Furthermore, the interphase between the matrix and inclusion is not necessarily heterogeneous but can exhibit functionally graded properties. Correspondingly, the temperature and normal heat flux jump relations [START_REF] Huy | Phénomènes de transmission a travers des couches minces de conductivité élevé e[END_REF] and ( 26) are given as follows:

• in the case of a flat interface Γ with the unit normal vector n = f 3 ,

K (c) (x 3 ) = K (c) γζ (x 3 )f γ ⊗ f ζ + K (c) 33 (x 3 )f 3 ⊗ f 3 and K (p) = K (p) γζ f γ ⊗ f ζ + K (p) 33 f 3 ⊗ f 3 with γ, ζ = 1 or 2, θ(x) = t 2 1 K (1) 33 - 1
K (c-) 33
q

(1)

3 (x) + 1 K (0) 33 - 1
K (c+) 33 q (0) 3 (x) , (29)
q 3 (x) = t 2 (K (c-) γζ -K (1)
γζ )

∂ 2 θ (1) (x) ∂x γ ∂x ζ + (K (c+) γζ -K (0) γζ ) ∂ 2 θ (0) (x) ∂x γ ∂x ζ ; (30)
• in the case of a circular interface Γ with radius R,

K (c) (r) = K (c) rr (r)f r ⊗ f r + K (c) φφ (r)f φ ⊗ f φ and K (p) = K (p) rr f r ⊗ f r + K (p) φφ f φ ⊗ f φ , θ(x) = t 2 1 K (1) rr - 1 K (c-) rr q (1) r (x) + 1 K (0) rr - 1 K (c+) rr q (0) r (x) , (31)
q r (x) = t 2   (K (c-) φφ -K (1) φφ ) R 2 ∂ 2 θ (1) (x) ∂φ 2 + (K (c+) φφ -K (0) rr ) R 2 ∂ 2 θ (0) (x) ∂φ 2   ; (32)
• in the case of a spherical interface Γ with radius R,

K (c) (r) = K (c) rr (r)f r ⊗ f r + K (c) φφ (r)f φ ⊗ f φ + K (c) ϕϕ (r)f ϕ ⊗ f ϕ with K (c) φφ = K (c) ϕϕ and K (p) = K (p) rr f r ⊗ f r + K (p) φφ f φ ⊗ f φ + K (p) ϕϕ f ϕ ⊗ f ϕ with K (p) φφ = K (p) ϕϕ , θ(x) = t 2 1 K (1) rr - 1 K (c-) rr q (1) r (x) + 1 K (0) rr - 1 K (c+) rr q (0) r (x) , (33)
q r (x) = t 2   (K (c-) φφ -K (1) φφ ) R 2 sin φ ∂ ∂φ sin φ ∂θ (1) (x) ∂φ + (K (c-) ϕϕ -K (1) ϕϕ ) R 2 sin φ ∂ 2 θ (1) (x) ∂ϕ 2 + (K (c+) φφ -K (0) φφ ) R 2 sin φ ∂ ∂φ sin φ ∂θ (0) (x) ∂φ + (K (c+) ϕϕ -K (0) ϕϕ ) R 2 sin φ ∂ 2 θ (0) (x) ∂ϕ 2   . (34)
In the above equations, (f r , f φ , f ϕ ) denotes a spherical orthonormal basis giving rise to a spherical coordinate system (r, φ, ϕ) whose origin coincides with the center of a spherical inclusion Ω (1) ; (f r , f φ ) is a polar orthonormal basis leading to the polar coordinate system (r, φ) with the origin situated at the center of a circular inclusion Ω (1) .

At the macroscopic scale, the composite medium under investigation is supposed to be statistically homogeneous. Furthermore, similar to the local thermal linear conduction law at the microscopic scale, the effective thermal conduction law is also assumed to be linear at the macroscopic scale and the corresponding effective thermal behavior is described by

Q = K eff • E (35)
where Q and E represent the macroscopic heat flux and intensity vectors defined in [START_REF] Le-Quang | Eshelby's tensor fields and effective conductivity of composites made of anisotropic phases with kapitza's interface thermal resistance[END_REF] and K eff corresponds to the effective thermal conductivity tensor of the composite to be determined.

In this section, attention is focused on the inclusion domain Ω (1) with imperfect interface Γ that is characterized by its thermal conductivity tensor K (1) and the interfacial jump conditions described by Eq. ( 6), ( 7) or ( 8). This inclusion Ω (1) with imperfect interface Γ is replaced now with an equivalent inclusion Ω(1) of same shape as Ω (1) and of unknown thermal conductivity tensor K(1) . At the same time, the imperfect interface Γ is substituted by a perfect interface, symbolized by Γ. For later notational convenience, the boundary of Ω( 1) is denoted by ∂ Ω(1) with ∂ Ω(1) = ∂Ω (1) . The determination of the unknown thermal conductivity tensor Ω(1) of the equivalent inclusion can be carried out by applying the generalized Hill-Mandel's lemma and by requiring that the thermal energy of the equivalent inclusion Ω(1) with perfect interface Γ be equal to be the one of the initial inclusion Ω (1) with the imperfect interface Γ.

First, the body Ω is assumed to be homogeneous and made of the matrix phase whose thermal conductivity tensor is K (0) . First, let Ω be subjected to either the uniform intensity (3) or uniform heat flux (4) or periodic boundary condition [START_REF] Hasselman | Effective thermal conductivity of composites with interfacial thermal barrier resistance[END_REF]. Due to the fact that Ω is homogeneous, the temperature field θ 0 (x), intensity field e 0 (x) and heat flux field q 0 (x) in Ω take the following simple forms:

θ 0 (x) = -E 0 • x, e 0 (x) = E 0 , q 0 (x) = K (0) • E 0 (36)
when the uniform intensity boundary condition (3) or periodic boundary condition ( 5) is concerned or

θ 0 (x) = -H (0) • Q 0 • x, e 0 (x) = H (0) • Q 0 , q 0 (x) = Q 0 (37)
when the uniform heat flux boundary condition ( 4) is under consideration. In Eq. ( 36) and hereafter, H (p) with p = 1 or 0 denotes the thermal resistivity tensor of phase p. We immediately obtain the thermal energy over Ω as follows:

U 0 = 1 2 |Ω| q 0 (x) • e 0 (x) = 1 2 |Ω|E 0 • K (0) • E 0 (38)
when the uniform intensity boundary condition (3) or periodic boundary condition ( 5) is prescribed or

U 0 = 1 2 |Ω| q 0 (x) • e 0 (x) = 1 2 |Ω|Q 0 • H (0) • Q 0 (39)
when the uniform heat flux boundary condition ( 4) is imposed.

Second, we cut the subdomain Ω (1) out of Ω and substitute back a inclusion phase Ω (1) with imperfect interface Γ and thermal conductivity tensor K (1) . Thus, after substituting the inclusion phase, the body Ω becomes heterogeneous. By applying the generalized Hill-Mandel's lemma described in section 2 to the particular case where the Hill-Mendel's condition ( 22) is satisfied, the macroscopic thermal energy over Ω after inserting the inhomogeneity Ω (1) with imperfect interface Γ is given by

U = 1 2 |Ω|E • Q = 1 2 |Ω| q(x) • e(x) - 1 2 Γ θ(x)q(x) • n(x)dx. (40)
Moreover, it can be shown that this macroscopic thermal energy U takes the following equivalent form

U = U 0 ± 1 2 Γ (+) q 0 • n(x)θ(x) -q(x) • n(x)θ 0 (x) dx + ∂Ω (1) \Γ q 0 • m(x)θ(x) -q(x) • m(x)θ 0 (x) dx (41)
where the positive sign " + " corresponds to the case where the uniform intensity boundary condition (3) or periodic boundary condition ( 5) is concerned while the negative sign " -" is relative to the case where the uniform heat flux boundary condition ( 4) is under consideration. In addition, Γ (+) denotes the side of Γ associated with the matrix phase Ω (0) ; θ(x) and q(x) are the temperature and heat flux solution fields inside the heterogeneous domain Ω with imperfect interface Γ.

Proof. Using Eqs. ( 38), ( 39) and ( 40), it implies that

2(U 0 -U ) = Ω (q 0 ± q) • (-∇θ 0 ± ∇θ)dx ± Ω (q • ∇θ 0 -q 0 • ∇θ)dx + Γ θ(x)q(x) • n(x)dx (42)
By taking into account Eq. ( 14) and by applying the divergence theorem together with the energy conservation equations ∇•q(x) = 0 and ∇•q 0 (x) = 0, we obtain

2(U 0 -U ) = ∂Ω (0) \Γ (q 0 ± q) • m(-θ 0 ± θ)dx - Γ (+) (q 0 ± q) • n(-θ 0 ± θ)dx + ∂Ω (1) \Γ (q 0 ± q) • m(-θ 0 ± θ)dx + Γ (-) (q 0 ± q) • n(-θ 0 ± θ)dx ± Ω (1)
(q • ∇θ 0 -q 0 • ∇θ)dx ±

Ω (0) (q • ∇θ 0 -q 0 • ∇θ)dx + Γ θ(x)q(x) • n(x)dx (43)
Next, it can be shown that:

• (i)
when Ω is subjected to either the uniform intensity (3) or uniform heat flux (4) or periodic boundary conditions (5), we have (q 0 ±q)•m(-θ 0 ±θ) = 0 on ∂Ω and 1) , we obtain

∂Ω (0) \Γ (q 0 ± q) • m(-θ 0 ± θ)dx + ∂Ω (1) \Γ (q 0 ± q) • m(-θ 0 ± θ)dx = ∂Ω (q 0 ± q) • m(-θ 0 ± θ)dx = 0; (44)
• (ii) since • = • (0) -• (
Γ (-) (q 0 ± q) • n(-θ 0 ± θ)dx - Γ (+) (q 0 ± q) • n(-θ 0 ± θ)dx = ± Γ ( q • n θ 0 -q 0 • n θ )dx - Γ θ(x)q(x) • n(x)dx; (45)
• (iii) due to the index symmetry of the thermal conductivity tensor K (0) such as

K (0) ij = K (0)
ji , it implies that

Ω (0) (q • ∇θ 0 -q 0 • ∇θ)dx = Ω (0) (∇θ • K (0) • ∇θ 0 -∇θ 0 • K (0) • ∇θ)dx = 0; (46)
• (iv) by using the divergence theorem with the energy conservation equations ∇ • q(x) = 0 and ∇ • q 0 (x) = 0, it reads

Ω (1) (q • ∇θ 0 -q 0 • ∇θ)dx = Γ (-) (q • nθ 0 -q 0 • nθ)dx + ∂Ω (1) \Γ (q • mθ 0 -q 0 • mθ)dx. (47)
Final, substitution of Eqs. ( 44)-( 47) into Eq. ( 43) leads to

2(U -U 0 ) = ± Γ (+) q 0 • n(x)θ(x) -q(x) • n(x)θ 0 (x) dx + ∂Ω (1) \Γ q 0 • m(x)θ(x) -q(x) • m(x)θ 0 (x) dx (48)
that is equivalent to Eq. ( 41).

Third, we consider the case where the inclusion Ω (1) with imperfect interface Γ is replaced with an equivalent inclusion Ω(1) of same shape as Ω (1) and of unknown thermal conductivity tensor K(1) and the imperfect interface Γ is substituted by a perfect interface Γ. Again by applying the simple form of the generalized Hill-Mandel's lemma [START_REF] Gu | Interfacial discontinuity relations for coupled multifield phenomena and their application to the modeling of thin interphases as imperfect interfaces[END_REF] where the Hill-Mendel's condition [START_REF] Hashin | Thin interphase/imperfect interface in conduction[END_REF] is verified, the corresponding macroscopic thermal energy, denoted by Ũ , is calculated by

Ũ = U 0 ± 1 2 Γ q 0 • n(x) θ(x) -q(x) • n(x)θ 0 (x) dx + ∂ Ω(1) \ Γ q 0 • m(x) θ(x) -q(x) • m(x)θ 0 (x) dx ( 49
)
where θ(x) and q(x) are the temperature and heat flux solution fields inside the heterogeneous domain Ω with perfect interface Γ. In addition, the local heat flux solution field q(x) is related to the local temperature field θ(x) by q(x) = -K(1) • ∇ θ(x) for x ∈ Ω (1) , q(x) = -K (0) • ∇ θ(x) for x ∈ Ω (0) . (50)

Finally, by demanding that the thermal energy Ũ after replacing the inclusion Ω (1) and imperfect interface Γ by the equivalent inclusion Ω(1) and perfect interface Γ be equal to the initial one U with the inclusion Ω (1) and imperfect interface Γ, it follows from Eqs. ( 49) and (50) that

Γ (+) q 0 • n(x)θ(x) -q(x) • n(x)θ 0 (x) dx + ∂Ω (1) \Γ q 0 • m(x)θ(x) -q(x) • m(x)θ 0 (x) dx = Γ q 0 • n(x) θ(x) -q(x) • n(x)θ 0 (x) dx + ∂ Ω(1) \ Γ q 0 • m(x) θ(x) -q(x) • m(x)θ 0 (x) dx (51)
holds for any E 0 or Q 0 . By substituting Eqs. ( 36), ( 39) and (50) into Eq. ( 51), the condition Ũ = U can be now detailed as follows:

Γ (+) K (0) • E 0 • n(x)θ(x)-K (0) • ∇θ(x) • n(x)E 0 • x dx + ∂Ω (1) \Γ K (0) • E 0 • m(x)θ(x)-K (1) • ∇θ(x) • m(x)E 0 • x dx = Γ K (0) • E 0 • n(x) θ(x) + q(x) • n(x)E 0 • x dx + ∂ Ω(1) \ Γ K (0) • E 0 • m(x) θ(x)-K(1) • ∇ θ(x) • m(x)E 0 • x dx (52)
which holds for any E 0 in the case where the uniform intensity boundary condition (3) or periodic boundary condition ( 5) is concerned, and

Γ (+) Q 0 • n(x)θ(x)-K (0) • ∇θ(x) • n(x)H (0) • Q 0 • x dx + ∂Ω (1) \Γ Q 0 • m(x)θ(x)-K (1) • ∇θ(x) • m(x)H (0) • Q 0 • x dx = Γ Q 0 • n(x) θ(x) + q(x) • n(x)H (0) • Q 0 • x dx + ∂ Ω(1) \ Γ Q 0 • m(x) θ(x)-K(1) • ∇θ(x) • m(x)H (0) • Q 0 • x dx (53)
which holds for any Q 0 in the case where the uniform heat flux boundary condition ( 4) is under consideration. The condition (52) or (53) allows us to calculate the thermal conductivity tensor K(1) of the equivalent inclusion Ω(1) .

In the particular case where the interface Γ is closed, i.e. ∂Ω (1) \Γ = ∂ Ω(1) \ Γ = ∅, the conditions (52) and (53) reduce therefore to

Γ (+) K (0) • E 0 • n(x)θ(x)-K (0) • ∇θ(x) • n(x)E 0 • x dx = Γ K (0) • E 0 • n(x) θ(x) + q(x) • n(x)E 0 • x dx (54)
for any E 0 when the uniform intensity boundary condition (3) or periodic boundary condition ( 5) is in question, or to

Γ (+) Q 0 • n(x)θ(x)-K (0) • ∇θ(x) • n(x)H (0) • Q 0 • x dx = Γ Q 0 • n(x) θ(x) + q(x) • n(x)H (0) • Q 0 • x dx (55)
for any Q 0 when the uniform heat flux boundary condition (4) is concerned.

In addition, when Ω is an infinite domain, it can be shown that the conditions (54) and ( 55) are equivalent to

Γ (+) K (0) • E 0 • n(x)θ(x)-K (0) • ∇θ(x) • n(x)E 0 • x dx = -E 0 • K (0) • Ω(1) [K (0) • (K (0) -K(1) ) -1 • K (0) -K (0) • SEsh ] -1 dx • K (0) • E 0 (56)
when the uniform intensity boundary condition (3) or periodic boundary condition ( 5) is concerned, or to

Γ (+) Q 0 • n(x)θ(x)-K (0) • ∇θ(x) • n(x)H (0) • Q 0 • x dx = -Q 0 • Ω(1) [K (0) • (K (0) -K(1) ) -1 • K (0) -K (0) • SEsh ] -1 dx • Q 0 (57)
when the uniform heat flux boundary condition ( 4) is under consideration.

In ( 56) and (57), SEsh denotes the Eshelby tensor field inside the inclusion Ω(1) . Specially, when the inclusion Ω(1) exhibits an ellipsoidal form for the three-dimensional case or an elliptic form for the bi-dimensional case and the matrix and inclusion phases are individually homogeneous, the Eshelby tensor field SEsh becomes uniform inside Ω(1) for any anisotropy of the material filling the infinite domain Ω. The expression of the Eshelby tensor SEsh of an ellipsoidal inclusion associated to an anisotropic media can be found in [START_REF] Le-Quang | Size-dependent eshelby tensor fields and effective conductivity of composites made of anisotropic phases with highly conducting imperfect interfaces[END_REF]. Consequently, the condition (56) and (57) take the following simple forms:

Γ (+) K (0) • E 0 • n(x)θ(x)-K (0) • ∇θ(x) • n(x)E 0 • x dx = -|Ω (1) |E 0 • K (0) • [K (0) • (K (0) -K(1) ) -1 • K (0) -K (0) • SEsh ] -1 • K (0) • E 0 , (58)
Γ (+) Q 0 • n(x)θ(x)-K (0) • ∇θ(x) • n(x)H (0) • Q 0 • x dx = -|Ω (1) |Q 0 • [K (0) • (K (0) -K(1) ) -1 • K (0) -K (0) • SEsh ] -1 • Q 0 (59)
with |Ω (1) | being the volume or area of Ω (1) .

5 Application of the equivalent inclusion method to estimating the effective thermal conductivity of composites

Layered composites with imperfect interfaces

As the first example of application, we consider a layered composites Ω consisting of n layers denoted by Ω (1) , Ω (2) , . . . , Ω (n) , respectively. The interface Γ (i) (i = 1, . . . , n -1) between two layers Ω (i) and Ω (i+1) is assumed to be imperfect and can be described by one of the three imperfect interface models specified in section 2. In addition, for simplicity, each layer Ω (i) is considered to be individually homogeneous, isotropic and of conductivity k i and thickness δ i (Fig. 1).

δ n n 1 n 2 1 K (1) K (n 2) K (n 1) K (n) Γ (1) Γ (n 2) Γ (n 1) K eff x 1 x 2 x 3 O Ω (1) Ω (n) Ω (n 1) Ω (n 2) δ δ δ i ∑δ i=1 n Imperfect interface Perfect interface 2 δ Γ (2) Ω (2) δ n n 1 n 2 1 K (1) K (n 2) K (n 1) K (n) Γ (1) Γ (n 2) Γ (n 1) Ω (1) Ω (n) Ω (n 1) Ω (n 2) δ δ δ 2 δ Γ (2) Ω (2) ~Ω (n) Ω (n 1) Ω (n 2) Ω (n) Ω (n 1) Ω (n 2) ~K (2) K (2) Fig 1.
Two-scale homogenization procedure for a layered composite in which each layer Ω (i) with imperfect interface Γ (i) is replaced in first with an equivalent layer Ω(i) with perfect interface Γ(i) , then this layered composite is homogenized by applying some classical schemes.

In order to determine the thermal conductivity tensor K(i) for the equivalent layer Ω (i) with imperfect interface Γ (i) , two new configurations of a layered composite Ω will be studied. In the first configuration, every layer except Ω (i) is replaced with a corresponding homogenized layer of the same thickness as Ω (i) and of unknown thermal conductivity tensor K eff (or thermal resistivity tensor H eff = (K eff ) -1 ). At the same time, each imperfect interface between them, except Γ (i) , is assumed to be perfect. The second configuration can be obtained from the first one by replacing the layer Ω (i) with imperfect interface Γ (i) by an equivalent layer Ω(i) of the same thickness as Ω (i) and an unknown thermal conductivity tensor K(i) (or unknown thermal resistivity tensor H(i) = ( K(i) ) -1 ). In the second configuration, the interface Γ(i) between Ω(i) and the remaining medium is perfect (Fig. 2).

K (n) Γ (n 1) x 3 O Ω (n) Γ (2) Ω (2) Ω (n) Ω (n 1) δ n n 1 1 K (i) K (n 1) x 2 δ δ 2 δ K (1) Γ (1) x 1 Ω (1) i δ Γ (n 2) K (n) Γ (n 1) K (2) Ω (i) i ∑δ i=1 n K eff K eff Γ (i) First configuration K (i) Ω (i) i ∑δ i=1 n K eff K eff Γ (i) K (i) Ω (i) Second configuration Γ (i) Imperfect interface Ω (i) Γ (i) K(i) Perfect interface
Fig2. First and second configurations used for determining the thermal conductivity tensor K(i) of the equivalent layer Ω(i)

We consider now two fundamental thermal conduction problems concerning a layered composite:

Problem 1: In-plane thermal conduction

In this problem, the layered composite Ω is be subjected to the uniform inplane intensity boundary condition (3) with E 0 1 = 0, E 0 2 = 0 but E 0 3 = 0 for both configurations. Under this boundary condition, it can be shown that the temperature solution fields θ(x) for the first configuration and θ(x) for the second configuration are given by

θ(x) = θ(x) = -E 0 1 x 1 -E 0 2 x 2 , x ∈ Ω. ( 60
)
It is interesting to remark from Eq. (60) that, even in the presence of the imperfect interface Γ (i) described either by HC, LC or GI model, the in-plane intensity components of the layered composite are also uniform and equal exactly to the ones applied on its boundary as in the classical case with perfect interface.

By substituting (60) into Eq. ( 52) and by replacing Γ, Γ, ∂Ω (1) , ∂ Ω(1) and K (0) with Γ (i) , Γ(i) , ∂Ω (i) , ∂ Ω(i) and K eff , respectively, we obtain from Eq. (52) that

K(i) γζ = K (i) γζ
for LC, HC or GI interface model.

It can be seen from Eq. (61) that, despite the fact that the interface Γ (i) is imperfect, the in-plane thermal conductivity tensor components K(i) γζ with γ, ζ = 1 or 2 of the equivalent layer Ω(i) are equal exactly to the ones of the layer Ω (i) .

Problem 2: Out-plane thermal conduction

In this case, the layered composite Ω undergoes the uniform out-plane heat flux boundary condition (4) with

Q 0 1 = 0, Q 0 2 = 0 but Q 0 3 = 0.
In the first configuration, under the boundary condition (4), by taking into account the imperfect interface conditions ( 6), ( 7) or ( 8) corresponding to the LC, HC or GI models, respectively, the temperature solution field is determined by

θ(x) =        -Q 0 3 H eff 33 x 3 for x ∈ Ω (1) ∪ Ω (2) ∪ • • • ∪ Ω (i-1) , -Q 0 3 H (i) 33 x 3 + C 1 for x ∈ Ω (i) , -Q 0 3 H eff 33 x 3 + C 2 for x ∈ Ω (i+1) ∪ Ω (i+2) ∪ • • • ∪ Ω (n) . (62)
Above, h i = 1/k i denotes the thermal resistivity of the layer Ω (i) ; H eff 33 designates the effective thermal resistivity tensor components of layered material under consideration; C 1 and C 2 are two constants determined from the interface conditions on Γ (i) as follows to within a constant:

C 1 = Q 0 3 H (i) 33 -H eff 33 i-1 p=1 δ p (63)
and

C 2 =          Q 0 3 δ i H eff 33 -δ i H (i) 33 -α i for LC interface model, Q 0 3 δ i H eff 33 -δ i H (i) 33
for HC interface model,

Q 0 3 δ i H eff 33 -δ i H (i) 33 + 1 2 t i H (i) 33 + 1 2 t i H eff 33 -t i H (ci) 33
for GI interface model, (64) with α i denoting the Kapitza thermal resistance of Γ (i) when the LC interface model is concerned, t i and H (ci) 33 being two interfacial parameters of the imperfect interface Γ (i) when the GI interface model is considered.

In the second configuration, when the layer Ω (i) is replaced with an equivalent layer Ω(i) of thermal conductivity tensor K(i) and thermal resistivity tensor H(i) and when all interfaces between layers become perfect, under the boundary condition (4), the temperature solution field takes the following simple expression to within a constant

θ(x) =        -Q 0 3 H eff 33 x 3 for x ∈ Ω (1) ∪ Ω (2) ∪ • • • ∪ Ω (i-1) , -Q 0 3 H(i) 33 x 3 + C1 for x ∈ Ω(i) , -Q 0 3 H eff 33 x 3 + C2 for x ∈ Ω (i+1) ∪ Ω (i+2) ∪ • • • ∪ Ω (n) , (65)
where C1 and C2 are given by

C1 = Q 0 3 H(i) 33 -H eff 33 i-1 p=1 δ p , (66)
C2 = Q 0 3 H eff 33 - H(i) 33 δ i . ( 67
)
Introducing the expressions of θ(x) and θ(x) given by Eqs. ( 62) and ( 65) together with Eqs. ( 63), ( 64), ( 66) and (67) into Eq. ( 53) and replacing Γ, Γ, ∂Ω (1) , ∂ Ω(1) and H (0) with Γ (i) , Γ(i) , ∂Ω (i) , ∂ Ω(i) and H eff , respectively, it follows that:

H(i) 13 = H(i) 23 = 0 for the LC, HC or GI interface model; (68)
• when the LC model is concerned,

H(i) 33 = H (i) 33 + α i δ i ; (69)
• when the HC model is under investigation,

H(i) 33 = H (i) 33 ; (70)
• when the GI model is considered,

H(i) 33 = H (i) 33 + H (ci) 33
t i δ i -H (i) 33 t i 2δ i -H eff 33 t i 2δ i . (71)
It is interesting to notice from Eqs. ( 69)-(71) that the thermal resistivity H(i)

33
in the layering direction of the equivalent layer Ω(i) is equal exactly to the one of Ω (i) for the HC interface model. However, for the LC and GI interface models, the thermal resistivity H(i) 33 of the equivalent layer Ω(i) depends not only on the thermal resistivity H (i) [START_REF] He | Homogenization of layered elastoplastic composites: Theoretical results[END_REF] of Ω (i) but also on the thickness δ i of the layer Ω (i) as well as on the interface parameters α i , t i and δ i . In particular, for the GI interface model, the thermal resistivity H(i) 33 of the equivalent layer Ω(i) depends also on the unknown effective resistivity H eff 33 in the layering direction of the homogenized layered material.

Finally, after replacing every layer Ω (i) and imperfect interface Γ (i) (i = 1, . . . , n -1) with an equivalent layer Ω(i) of thermal equivalent conductivity tensor K(i) or resistivity tensor H(i) , whose the tensor components given by Eqs. ( 61), ( 68)-(71), we obtain a new layered material in which all interfaces are perfect. As a consequence, this layered composite can be homogenized by applying the classical theory for laminates (see e.g. [START_REF] Milton | The Theory of Composites[END_REF] and [START_REF] He | Homogenization of layered elastoplastic composites: Theoretical results[END_REF]). The nonzero macroscopic thermal conductivity tensor components can be exactly and analytically determined as follows:

K eff γζ = n i=1 f i K (i) γζ
for the LC, HC or GI interface model;

• for the lowly conducting (LC) interface model,

K eff 33 = (H eff 33 ) -1 = n i=1 f i H (i) 33 + n-1 i=1 f i α i δ i -1 ; (73)
• for the highly conducting (HC) interface model,

K eff 33 = (H eff 33 ) -1 = n i=1 f i H (i) 33 -1 ; (74)
• for the general imperfect (GI) interface,

K eff 33 = (H eff 33 ) -1 = 1 + n-1 i=1 t i 2δ i n i=1 f i H (i) 33 + n-1 i=1 f i t i δ i H (ci) 33 - 1 2 H (i) 33 -1 ; ( 75
)
where f i = δ i / n i=1 δ i stands for the volume fraction of the i-th layer Ω (i) . The expression (73) of the effective thermal conductivity K eff 33 = (H eff 33 ) -1 in the layering direction for the LC interface model can be recovered from Eq. (75) by setting H (ci) 33 = α i /t i and by calculating the limit as t i → 0. Similarly, by setting H (ci) 33 = t i /β i into Eq. (75) and by computing the limit as t i → 0, the expression (74) of K eff 33 = (H eff 33 ) -1 for the HC interface model can be also found.

Composites with circular or spherical inclusions and imperfect interfaces

In the second example, the composite under consideration consists of a matrix in which circular inclusions in the two-dimensional case (d = 2) or spherical inclusions in the three-dimensional case (d = 3) are embedded. As before, we denote by Ω the d-dimensional (d = 2 or 3) domain occupied by a representative volume or area element of the composite material. We designate by Ω (0) and Ω (i) (with i = 1, . . . , n) the subdomains of Ω occupied by the matrix phase and i-th inclusion phase. The interface Γ (i) between the inclusion Ω (i) and matrix Ω (0) is assumed to be imperfect and modeled by the LC, HC or GI interface model. In addition, the matrix phase Ω (0) is considered to be individually homogeneous, isotropic and of conductivity k 0 while the i-th inclusion phase Ω (i) is supposed to be either circularly orthotropic with thermal conductivity tensor

K (i) = K (i) rr f r ⊗ f r + K (i)
φφ f φ ⊗ f φ when the 2D case is considered or spherically transversely isotropic with thermal conductivity tensor

K (i) = K (i) rr f r ⊗ f r + K (i) φφ f φ ⊗ f φ + K (i) ϕϕ f ϕ ⊗ f ϕ with K (i) φφ = K (i)
ϕϕ when the 3D case is concerned. The radius of the circular or spherical inclusion Ω (i) is denoted by R i (Fig. 3).

K eff K (0) K (i) R i ~K (i) K(2) ~K(1) K(0) K (2) K (1) Ω (n) Ω (n) Ω (0)
Ω( 1)

Ω (1) Ω (2) Ω (2) Ω (2) Ω (2) Ω (i) ~Ω (1) Γ (2) Γ (i) ~K (i) R i

Imperfect interface

Perfect interface Fig 3. Two-step homogenization procedure for a spherical particle reinforced composite: (i) each spherical inclusion Ω (i) with imperfect interface Γ (i) is replaced by an equivalent spherical inclusion Ω(i) with perfect interface Γ(i) ; (ii) the resulting spherical particle reinforced composite is homogenized by applying classical schemes.

As in section 4, in order to replace the circular or spherical inclusion Ω (i) with imperfect interface Γ (i) by an equivalent inclusion Ω(i) whose thermal conductivity is ki and whose interface Γ(i) with the matrix is perfect, we consider now the first configuration in which a circular or spherical inclusion Ω (i) of radius R i and thermal conductivity tensor K (i) is inserted into a d-dimensional infinitely extended matrix Ω of thermal conductivity k 0 . As before, the interface Γ (i) between the inclusion Ω (i) and matrix phase Ω (0) is supposed to be imperfect and described either by the LC, HC or GI interface model.

On the surface ∂Ω of Ω, the uniform intensity boundary condition (3) is prescribed with E 3 = E 0 = 0 for the three-dimensional case (d = 3). Relative to either the system of spherical coordinates (r, φ, ϕ) rela-tive to the spherical orthogonal basis (f r , f φ , f ϕ ) whose origin coincides with the center of Ω (i) for the 3D case or the system of polar coordinates (r, φ) relative to the polar orthogonal basis (f r , f φ ) whose origin is situated at the center of Ω (i) for the 2D case, the temperature solution field takes the following general form:

θ(x) =    A 1 r ξ i cos φ for x ∈ Ω (i) , A 2 r + B 2 r 1-d cos φ for x / ∈ Ω (i) (76)
where

ξ i = (d -2) 2 4 + (d -1) K (i) φφ K (i) rr - (d -2) 2 > 0
and r = x and three constants A 1 , A 2 and B 2 are determined by using the boundary condition (3) at r = ∞ and the interface conditions ( 6), ( 7) or (8) at r = R i . Precisely, they are given as follows:

A 2 = -E 0 for the LC, HC or GI interface model (77)

• for the LC interface model,

A 1 = - dk 0 E 0 R 2-ξ i i ξ i K (i) rr + (d -1)k 0 R i + (d -1)α i ξ i K (i) rr k 0 , B 2 = - (k 0 -ξ i K (i) rr )R i + α i k 0 ξ i K (i) rr E 0 R d i ξ i K (i) rr + (d -1)k 0 R i + (d -1)α i ξ i K (i) rr k 0 ; (78)
• for the HC interface model,

A 1 = - dk 0 E 0 R 2-ξ i i ξ i K (i) rr + (d -1)k 0 R i + (d -1)β i , B 2 = - (k 0 -ξ i K (i) rr )R i -(d -1)β i E 0 R d i ξ i K (i) rr + (d -1)k 0 R i + (d -1)β i ; (79)
• for the GI interface model,

A 1 = -E 0 K (ci-) rr d 4R 1-ξ i i k 0 K (ci+) rr + t i R i 2 R 1-ξ i i (d -1)(K (ci+) φφ -k 0 )(K (ci+) rr -k 0 ) × t i R i 2 (1 -d) (d -1)(K (i) φφ -K (ci-) φφ )K (ci-) rr (k 0 -K (ci+) rr ) +ξ i (k 0 -K (ci+) φφ )K (ci+) rr (K (i) rr -K (ci-) rr ) + 2(d -1) t i R i K (ci+) rr K (ci-) rr (K (ci-) φφ +K (ci+) φφ -K (i) φφ -k 0 (ξ i + 1) -K (i) rr ξ i ) + K (i) rr ξ i k 0 (K (ci+) rr + K (ci-) rr ) +4K (ci+) rr K (ci-) rr K (i) rr ξ i + k 0 (d -1) -1 , B 2 = E 0 R d i 4(K (i) rr ξ i -k 0 )K (ci+) rr K (ci-) rr + 2 t i R i (d -1)K (ci+) rr K (ci-) rr (K (ci+) φφ + K (ci-) φφ -k 0 -K (i) φφ ) -K (ci+) rr K (ci-) rr ξ i (k 0 + K (i) rr ) -k 0 K (i) rr ξ i (K (ci+) rr + K (ci-) rr ) +(d -1) t i R i 2 K (ci+) rr K (ci-) rr (K (ci-) φφ -K (i) φφ + ξ i k 0 -ξ i K (ci+) φφ ) +k 0 (K (i) φφ -K (ci-) φφ )K (ci-) rr + K (ci+) rr ξ i K (i) rr (K (ci+) φφ -k 0 ) × t i R i 2 (1 -d) (d -1)(K (i) φφ -K (ci-) φφ )K (ci-) rr (k 0 -K (ci+) rr ) +ξ i (k 0 -K (ci+) φφ )K (ci+) rr (K (i) rr -K (ci-) rr ) + 2(d -1) t i R i K (ci+) rr K (ci-) rr (K (ci-) φφ +K (ci+) φφ -K (i) φφ -k 0 (ξ i + 1) -K (i) rr ξ i ) + K (i) rr ξ i k 0 (K (ci+) rr + K (ci-) rr ) +4K (ci+) rr K (ci-) rr K (i) rr ξ i + k 0 (d -1) -1 . ( 80
)
The second configuration is obtained from the first one by replacing the circular or spherical inclusion Ω (i) with the imperfect interface Γ (i) by an equivalent inclusion Ω(i) of the shape identical to Ω (i) and with the equivalent thermal conductivity ki and the perfect interface Γ(i) . Under the same boundary condition on ∂Ω as in the first configuration, the temperature solution field can be obtained as follows:

θ(x) =    Ã1 r cos φ for x ∈ Ω(i) , Ã2 r + B2 r 1-d cos φ for x / ∈ Ω(i) . ( 81
)
In the above expression, Ã1 , Ã2 and B2 are three constants to be determined from the remote boundary condition (3) at r = ∞ and the continuity conditions of the temperature and normal heat flux field component across the perfect interface Γ(i) at r = R i . They are provided by

Ã2 = -E 0 , Ã1 = - dk 0 E 0 R i ki + (d -1)k 0 R i , B2 = - (k 0 -ki )E 0 R d i ki + (d -1)k 0 . (82) K (i) R i K (i) R i Imperfect interface Perfect interface Ω (i) K (i) Γ (i) K (0) Ω (0) Γ (i) Ω (i) K (i) K(0) Ω (0) ~Ω (i) Ω (i) ~First configuration
Second configuration

Γ (i) Γ (i)
Fig4.

First and second configurations used for determining the thermal conductivity tensor K(i) of the equivalent spherical inclusion Ω(i)

On the other hand, when Ω is an infinite domain and when the matrix/inclusion interface Γ(i) is closed and perfect, the determination of the thermal conductivity ki of the equivalent inclusion Ω(i) can be carried out via the Eshelby's tensor SEsh which takes the simple form

SEsh = 1 d I (d) (83)
under the condition that Ω(i) is circular or spherical.

By inserting the expressions of θ(x) and θ(x) or SEsh provided by Eqs. ( 76) and ( 81) or (82) together with Eqs. ( 77)-( 80) and (82) into Eq. ( 54) or (56) with K (0) = k 0 I (d) and q(x) = -ki ∇ θ(x) and by replacing Γ and Γ with Γ (i) and Γ(i) in Eq. ( 54) or (56), respectively, we obtain the following results:

• when the lowly conducting (LC) interface model is concerned,

ki = ξ i K (i) rr 1 + α i R i ξ i K (i) rr
or equivalently hi = 1

ξ i K (i) rr + α i R i ; ( 84
)
• when the highly conducting (HC) interface model is under investigation, ki = ξ i K (i) rr + (d -1)

β i R i ; ( 85
)
• when the general imperfect (GI) interface model is considered,

ki = (8 -2d)k 0 K (ci+) rr K (ci-) rr ξ i K (i) rr + 2k 0 K (ci+) rr K (ci-) rr t i R i (K (ci+) φφ +K (ci-) φφ -k 0 -K (i) φφ ) + (K (i) rr -K (ci-) rr )(K (ci+) φφ -k 0 )k 0 K (ci+) rr ξ i t i R i 2 × (8 -2d)k 0 K (ci+) rr K (ci-) rr + (4 -d) t i R i ξ i K (i) rr k 0 (K (ci+) rr + K (ci-) rr ) -K (ci+) rr K (ci-) rr (k 0 + K (i) rr ) + (K (ci-) φφ -K (i) φφ )(k 0 -K (ci+) rr )K (ci-) rr t i R i 2 -1 . (86)
As in the case of layered material, the thermal conductivity ki of the equivalent inclusion Ω(i) provided by (84) for the LC case and by (85) for the HC case can be recovered from Eq. ( 76) by setting

K (ci±) rr = K (ci±) φφ = t i /α i for the LC interface model and K (ci±) rr = K (ci±) φφ = β i /t i
for the HC interface model and by calculating the limit of Eq. ( 86) as t i → 0.

Next, owing to the fact that the inclusion/matrix interface is now perfect, we can apply any classical homogenization scheme to estimate the effective thermal conductivity tensor of the composite under consideration. Thus, this section, relative to the second-step homogenization procedure, consists in obtaining the closed-form expressions for the effective conductivity tensor by using some well-known classical homogenization schemes such as the dilute distribution, Mori-Tanaka, self-consistent, generalized self-consistent and differential approximation schemes. For more details about these schemes, the reader can refer to the works [START_REF] Hashin | Assessment of the self consistent scheme approximation: Conductivity of particulate composites[END_REF][START_REF] Hashin | Analysis of composite materials -a survey[END_REF][START_REF] Christensen | Solutions for effective shear properties in three phase sphere and cylinder models[END_REF][START_REF] Norris | A differential scheme for the effective moduli of composites[END_REF]. Denoting by k DD , k MT , k SC , k GSCS and k DA the effective thermal conductivity tensors derived from the dilute distribution, Mori-Tanaka, self-consistent, generalized self-consistent and differential approximation schemes, respectively, we obtain the following results.

• Dilute distribution (DD) scheme:

k DD = k 0 + n i=1 f i ki -k 0 dk 0 ki + (d -1)k 0 ; (87)
• Mori-Tanaka (MT) scheme or Generalized: self-consistent scheme (GSCS):

k M T = k GSCS = k 0 + n i=1 f i ki -k 0 dk 0 (1 -n i=1 f i ) ki + (d -1)k 0 + n i=1 f i dk 0 ; (88)
• Self-consistent (SC) scheme:

k SC = k 0 + n i=1 f i ki -k 0 dk SC ki + (d -1)k SC ; ( 89
)
The effective thermal conductivity k SC is calculated as the real positive root of Eq.( 89). • Differential approximation (DA) scheme:

Dk Dη = 1 1 -η n i=1 f i n i=1 f i ( ki -k)dk ki + (d -1)k with k(0) = k 0 . ( 90
)
The effective thermal conductivity k DA is determined as k DA = k(1) with k(η) being the solution of the differential equation (90). In particular case of a two-phase composite, or equivalently n = 1, the effective thermal conductivity k DA corresponds to the real positive root of the following d-th order equation

k 0 k1 -k DA d -k DA k1 -k 0 d (1 -f 1 ) d = 0. ( 91
)
The last part of this section is dedicated to numerically illustrating the results obtained above for the effective thermal conductivity of composites with imperfect interfaces. More precisely, in this numerical example, the two-dimensional composite is assumed to consist of a host isotropic matrix phase of conductivity k 0 in which mono-sized circular and isotropic inclusions of radius R 1 and conductivity k 1 are randomly distributed. In addition, the inclusion phase conductivity k 1 is chosen in such a way that k 1 /k 0 = 10 and k 1 /k 0 = 0.1 for two cases where the inclusion phase is more and less conducting than the matrix phase. The interfaces between the matrix and inclusions are described either by the HC, LC or GI interface model. In particular, for the GI interface model, the thickness t 1 is kept constant and is such that t 1 /R 1 = 0.001 whereas the conductivity

K (c1) rr (r) = K (c1) φφ (r) = k (c1) is set to vary.
Starting with the GI interface model, the normalized effective thermal conductivities with respect to the thermal conductivity of the matrix phase, obtained by applying the dilute distribution, Mori-Tanaka, self-consistent, generalized self-consistent and differential approximation schemes, are plotted versus the inclusion area fraction f 1 in Fig. 5 for the case where the matrix phase is less conducting than the inclusion and interphase and in Fig. 6 for the contrary case where the matrix phase is more conducting than the inclusion and interphase. In the same figures, these estimates for the effective thermal conductivity are then compared with the corresponding first-order upper and lower bounds (1stUB and 1stLB) as well as the generalized Hashin-Shtrikman lower and upper bounds which are derived and provided in [START_REF] Le-Quang | Estimations and bounds of the effective conductivity of composites with anisotropic inclusions and general imperfect interfaces[END_REF]. It can be seen from Figs. 5 and 6 that:

• The effective thermal conductivities k M T and k GSCS obtained by MT scheme and GSCS coincide with the Hashin-Shtrikman bound k HS . It is important to notice that this Hashin-Shtrikman bound k HS represents the lower one in the first case where the matrix phase is less conducting than the inclusion and interphase while it corresponds to the upper one in the second case where the matrix phase is more conducting than the inclusion and interphase. • Unlike the approximation of the effective thermal conductivity derived by the dilute distribution scheme, the values of the effective thermal conductivity based on the differential approximation scheme are situated always between the first-order upper and lower bounds as well as the second-order bounds. In other words, the effective thermal conductivity obtained from the differential approximation scheme never violates the first and secondorder upper and lower bounds. For this reason, except the MT scheme and GSCS, the estimation by applying the differential scheme can be also considered as an excellent approximation for the effective thermal conductivity of composites.

Concerning the LC interface model, we plot in Figs. 7 and8 the values of the normalized effective thermal conductivity provided by the dilute distribution, Mori-Tanaka, self-consistent, generalized self-consistent and differential approximation schemes in terms of the inclusion area fraction f 1 for a two-phase composite in which the inclusion/matrix interfaces are described by the LC interface model. At the same time, these values for the normalized effective thermal conductivities obtained with the LC interface model are compared with the corresponding ones derived with the GI interface model in which the conductivity k (c1) and thickness t 1 of the interphase are chosen in such a way that α 1 = t 1 /k (c1) . In addition, the results obtained and shown in Fig. 8 correspond to the case where the conductivity of interphase k (c1) is very lowly conducting with respect to the conductivities of the matrix and inclusion phases (k (c1) /k 0 = 0.001) while, in Fig. 7, the conductivity k (c1) of the interphase takes the same order of magnitude as the conductivities of the matrix and inclusion phases. It can be observed from Fig. 8 that the effective thermal conductivities obtained with the LC interface model has the same values as the one with the GI interface model when the interphase with thermal conductivity k (c1) is very lowly conducting. However, as can be seen from Fig. 7, this observation is no longer true when k (c1) is of the same order of magnitude as the conductivities of the matrix and inclusion phases.

In a similar way, the normalized effective thermal conductivities of the twophase composite with the HC interface model calculated by applying the dilute distribution, Mori-Tanaka, self-consistent, generalized self-consistent and differential approximation schemes are plotted in terms of the inclusion area fraction f 1 in Figs. 9 and 10. In addition, the values of the normalized effective thermal conductivities derived with the GI interface model are also represented in these figures for comparison with the ones obtained with the HC interface model. The surface conductivity β 1 for the HC interface model is related to the conductivity k (c1) and the thickness t 1 of the interphase for GI interface model by β 1 = k (c1) t 1 . It can be seen from Figs. 9 and 10 that the effective thermal conductivities with the HC interface model is recovered from the GI interface model when the interphase with thermal conductivity k (c1) is very highly conducting.

Thus, the GI model includes the LC and HC interface models as particular cases. More precisely, the HC interface model characterized by a surface thermal conductivity β 1 can be recovered from the GI model by considering that the interphase between the matrix and inclusions has a very small and uniform thickness t 1 but is a very conductivity given by k (c1) = β 1 /t 1 . Similarly, the LC interface model characterized by a Kapitza's interface thermal resistance α 1 can be obtained since the GI model by requiring the interphase of small thickness t 1 to have a very low conductivity k (c1) = t 1 /α 1 .

Inclusion area fraction Fig5. Estimates and bounds of the normalized effective thermal conductivity, with respect to the conductivity of the matrix phase k 0 , versus the inclusion area fraction of a two-phase composite with general imperfect (GI) interface model with k 1 /k 0 = 10, k (c1) /k 0 = 5 and t 1 /R 1 = 0.001. Fig6. Estimates and bounds of the normalized effective thermal conductivity, with respect to the conductivity of the matrix phase k 0 , versus the inclusion area fraction of a two-phase composite with general imperfect (GI) interface model with

k 1 /k 0 = 0.1, k (c1)
/k 0 = 0.5 and t 1 /R 1 = 0.001. Fig7. Estimates and bounds of the normalized effective thermal conductivity, with respect to the conductivity of the matrix phase k 0 , versus the inclusion area fraction of a two-phase composite with general imperfect (GI) and lowly conducting (LC) interface models with k 1 /k 0 = 0.1, k (c1) /k 0 = 10, t 1 /R 1 = 0.001 and Fig8. Estimates and bounds of the normalized effective thermal conductivity, with respect to the conductivity of the matrix phase k 0 , versus the inclusion area fraction of a two-phase composite with general imperfect (GI) and lowly conducting (LC) interface models with k 1 /k 0 = 0.1, k (c1) /k 0 = 0.001, t 1 /R 1 = 0.001 and Fig9. Estimates and bounds of the normalized effective thermal conductivity, with respect to the conductivity of the matrix phase k 0 , versus the inclusion area fraction of a two-phase composite with general imperfect (GI) and highly conducting (HC) interface models with k 1 /k 0 = 10, k (c1) /k 0 = 0.1, t 1 /R 1 = 0.001 and Fig10. Estimates and bounds of the normalized effective thermal conductivity, with respect to the conductivity of the matrix phase k 0 , versus the inclusion area fraction of a two-phase composite with general imperfect (GI) and highly conducting (HC) interface models with k 1 /k 0 = 10, k (c1) /k 0 = 1000, t 1 /R 1 = 0.001 and Estimates, bounds and comparaisons with results of Herve and Zaoui [START_REF] Herve | n-Layered inclusion-based micromechanical modelling[END_REF] and Garboczi and Bentz [START_REF] Garboczi | Analytical formulas for interfacial transition zone properties[END_REF] of the normalized effective thermal conductivity, with respect to the conductivity of the matrix phase k 0 , versus the inclusion volume fraction of a tridimensional two-phase composite with general imperfect (GI) interface model with k 1 /k 0 = 10, k (c1) /k 0 = 5 and t 1 /R 1 = 0.001. Fig14. Estimates and comparaisons with results of Sevostianov and Kachanov [START_REF] Sevostianov | Effect of interphase layers on the overall elastic and conductive properties of matrix composites. Applications to nanosize inclusion[END_REF] of the normalized effective thermal conductivity, with respect to the conductivity of the matrix phase k 0 , versus the inclusion volume fraction of a tridimensional two-phase composite with general imperfect (GI) interface model with k 1 /k 0 = 10, k

α 1 /R 1 = t 1 /(k (c1) R 1 ) = 0.0001k 0 .
α 1 /R 1 = t 1 /(k (c1) R 1 ) = k 0 .
β 1 /R 1 = t 1 k (c1) /R 1 = 0.0001k 0 .
β 1 /R 1 = t 1 k (c1) /R 1 = k 0 .
if /k 0 = 0.7 and t 1 /a 1 = 0.05. Fig17. Comparison between the effective thermal conductivities obtained by applying the EIM and three-phase model with their estimates and bounds for random and periodic composites whose k 1 /k 0 = 0.1, k (c1) /k 0 = 0.001 and t 1 /R 1 = 0.001.

L L t 1 k (c1) k 1 R 1 k 0 L L L L √3 t 1 k (c1) k 1 R 1 k 0 k 1 k 1 k 1 k 1 t 1 k (c1) k 1 R 1 k 0 L L k 1 R 1 k 0 L L √3 k 0 k 1 R 1 k 1 k 1 k 1 k 1 (a) (b) L L k 0 k 1 k 1 k R 1 1 k 1 k 1 k 1 Fig 15.
The second numerical example is related to a tridimensional composite consisting of a host isotropic matrix of conductivity k 0 in which identical spherical and isotropic inclusions of radius R 1 and conductivity k 1 are randomly inserted via interfaces described by the general imperfect (GI) interface model. Unlike the first example, in this GI interface model considered as replacing a thin interphase, two cases where the interphase situated between the matrix and inclusion phases possesses constant and functionally graded properties will be considered. In the first case, we choose the thickness t 1 and the conductivity

K (c1) rr (r) = K (c1) φφ (r) = K (c1) ϕϕ (r) = k (c1)
of the interphase to be such that t 1 /R 1 = 0.001 and k (c1) /k 0 = 5. The inclusion phase is assumed to be more conducing than the matrix phase with k 1 /k 0 = 10. By using the general imperfect (GI) interface model and by applying the dilute distribution, Mori-Tanaka, self-consistent, generalized self-consistent and differential approximation schemes, the normalized effective thermal conductivities with respect to the conductivity of the matrix phase k 0 are plotted in Fig. 11 versus the inclusion volume fraction. These values obtained for effective thermal conductivities are then compared with their first-and second-order upper and lower bounds as well as with the corresponding estimates provided by Herve and Zaoui [START_REF] Herve | n-Layered inclusion-based micromechanical modelling[END_REF] (denoted by kHZ) and Garboczi and Bentz [START_REF] Garboczi | Analytical formulas for interfacial transition zone properties[END_REF]. It can be observed from Fig. 11 that: (i) the estimate presented by Herve and Zaoui coincides exactly with our results obtained by the Mori-Tanaka or generalized self-consistent scheme; (ii) the estimate of Garboczi and Bentz is identical to our results derived with the dilute distribution scheme. These comparisons confirm that the method proposed in our work to estimate the effective thermal conductivity of composites with imperfect interfaces can recover the results obtained by Herve and Zaoui [START_REF] Herve | n-Layered inclusion-based micromechanical modelling[END_REF] and Garboczi and Bentz [START_REF] Garboczi | Analytical formulas for interfacial transition zone properties[END_REF]. In the second case, the interphase between the matrix and inclusion phases is assumed to exhibit functionally graded properties. More precisely, the conductivities of the interphase is supposed to be expressed in the following power-law function (see. e.g. Lutz and Zimmerman [START_REF] Lutz | Effect of the Interphase Zone on the Bulk Modulus of a Particulate Composite[END_REF][START_REF] Lutz | Effect of an inhomogeneous interphase zone on the bulk modulus and conductivity of a particulate composite[END_REF]):

K (c1) rr (r) = K (c1) φφ (r) = K (c1) ϕϕ (r) = k 0 + (k (c1) if -k 0 ) r a 1 ρ 1 . (92)
Here, a 1 = R 1 -t/2 is the radius of the inclusion domain, k

c1) if (
denotes the conductivity of the interphase related to the interface associated with the inclusion domain, i.e. k

(c1) if = K (c1-) rr = K (c1-) φφ = K (c1-)
ϕϕ , and the parameter ρ 1 controlling the thickness t of the interphase is chosen to be t 1 = 2.3a 1 /ρ 1 (see e.g. [START_REF] Lutz | Inhomogeneous interfacial transition zone model for the bulk modulus of mortar[END_REF]). We show, in Figs. 12 and13, the variation of the normalized effective thermal conductivity of the composite obtained by using the dilute distribution, Mori-Tanaka, self-consistent, generalized self-consistent and differential approximation schemes in terms of the inclusion volume fraction for the two cases of D = (k 0 -k

(c1) if )/k 0 = -0.75 and D = (k 0 -k (c1)
if )/k 0 = 0, respectively. In the same figures, the values obtained for the normalized effective thermal conductivity are compared with the ones derived by Lutz and Zimmerman [START_REF] Lutz | Effect of an inhomogeneous interphase zone on the bulk modulus and conductivity of a particulate composite[END_REF]. It is very interesting to notice from Figs. 12 and 13 that our values for the normalized effective thermal conductivity obtained by the self-consistent and differential approximation schemes are very close to the ones provided by Lutz and Zimmerman. These good agreements observed in Figs. 12 and 13 between our results and the ones of Lutz and Zimmerman confirm the validity of the analytical expressions (87)-(91). In addition, compared with the method proposed by Lutz and Zimmerman in which several complicated computations with hypergeometric functions are needed, the analytical expressions (87)-(91) derived in the present work allow the easy calculation of the effective thermal conductivity of both 2D and 3D composites.

By adopting the expression (92) for the thermal conductivity of the interphase between the matrix and inclusion phases and by setting k 1 /k 0 = 10, k

if /k 0 = 0.7 and t 1 /a 1 = 0.05, we compare in Fig. 14 the values of the effective thermal conductivity derived by the dilute distribution, Mori-Tanaka, self-consistent, generalized self-consistent and differential approximation schemes with the estimates, Hashin-Shtrikman upper and lower bounds obtained by Sevostianov and Kachanov [START_REF] Sevostianov | Effect of interphase layers on the overall elastic and conductive properties of matrix composites. Applications to nanosize inclusion[END_REF] with non-interaction assumption and with effective field method. It can be observed from Fig. 14 (i) when the interaction between the inclusions is omitted, the effective thermal conductivity obtained with dilute distribution scheme is very close to the one provided by Sevostianov and Kachanov; (ii) the values of the effective thermal conductivity with selfconsistent and differential approximation schemes are well situated between the Hashin-Shtrikman upper and lower bounds; (iii) there is a good accord between the estimate of the effective thermal conductivity by applying the differential approximation scheme and the one of Sevostianov and Kachanov with effective field method. This observation confirms that the differential approximation scheme is one of the best approximate schemes for estimating the effective thermal conductivity of composites with randomly distributed inclusions.

Next, in order to test the validity of the equivalent inclusion method proposed in sections 4 and 5 of the present work, we consider periodic composites as described by the two configurations of Fig. 15. In the first configuration (Fig. 15a), called three-phase one, three typical periodic two-dimensional composites in which circular inclusions of the same radius and thermal conductivity k 1 are squarely, hexagonally and randomly distributed in a matrix phase of thermal conductivity k 0 via an interphase of thickness t 1 and thermal conductivity k (c1) . The interphase between the inclusion and matrix phases is considered to be uniform and very thin with respect to the inclusion size. Denoting by R 1 the radius of the midline of the interphase, the thickness t 1 is chosen such as t 1 /R 1 = 0.001. In the second configuration (Fig. 15b), namely equivalent inclusion method (EIM), each circular inclusion coated with an interphase is replaced by an equivalent inclusion of radius R 1 and thermal conductivity k1 determined from Eq. ( 86) with d = 2, i = 1 and

K (c1)
rr (r) = K (c1) φφ (r) = k (c1) . The effective thermal conductivity of periodic composites with squarely, hexagonally and randomly distributed inclusions are then numerically computed by applying the method based on the fast Fourier transform (FFT) together with an iterative approach developed by Moulinec [START_REF] Moulinec | A fast numerical method for computing the linear and nonlinear mechanical properties of composites[END_REF], Moulinec and Suquet [START_REF] Moulinec | A numerical method for computing the overall response of nonlinear composites with complex microstructure[END_REF], Bonnet [START_REF] Bonnet | Effective properties of elastic periodic composite media with fibers[END_REF] for elastic problems and Le Quang et al. [START_REF] Le-Quang | Effective thermal conductivity of periodic composites with highly conducting imperfect interfaces[END_REF] for thermal conduction ones. We plot, in Figs. 16 and17, the values of the effective thermal conductivities of periodic composites with squarely, hexagonally and randomly distributed inclusions versus the inclusion area fraction for both configurations with the three-phase model and equivalent inclusion model in the case where k 1 /k 0 = 0.1, t 1 /R 1 = 0.001 and the ratio k (c1) /k 0 is set to vary from k (c1) /k 0 = 10 to k (c1) /k 0 = 0.001. In particular, for composites with randomly distributed inclusions, the number of inclusions is set to be equal to 10 per unit cell and the values obtained for effective thermal conductivities K ef f 11 and K ef f 22 correspond to the its average values over 10 realizations. The values given for the effective thermal conductivities of periodic composites with the three-phase model are compared first with the corresponding ones provided by the equivalent inclusion method and second with the estimates and bounds established in section 5. It is interesting to observe from Figs. [START_REF] Miloh | On the effective conductivity of composites with ellipsoidal inhomogeneities and highly conducting interfaces[END_REF] and 17 that:

• The computation of periodic two-phase composites with the equivalent inclusion method is very simple and fast with respect to the one with periodic three-phase composites. The effective thermal conductivities obtained with the three-phase model and EIM for all three typical periodically composites with squarely, hexagonally and randomly distributed inclusions are indistinguishable. This confirms the validity and advantage of the equivalent inclusion method proposed in this work. • The effective thermal conductivities derived for periodically composites with squarely and hexagonally distributed inclusions are very close to the ones obtained by applying the Mori-Tanaka model or GSCS. This means that the Mori-Tanaka model or GSCS is the best approximate scheme to estimate the effective thermal conductivity of periodically composites with squarely and hexagonally distributed inclusions. • As aforementioned, in the case of k 1 /k 0 = 0.1 and k (c1) /k 0 = 0.001 where the matrix phase is more conducting than the inclusion phase as well as the interphase, the effective thermal conductivities provided from the Mori-Tanaka model or GSCS correspond exactly to the Hashin-Shtrikman upper bound. It is clear from Fig. 17 that the effective thermal conductivities of periodically composites with randomly distributed inclusions obtained with three-phase model and EIM respect well the Hashin-Shtrikman upper bound. • The effective thermal conductivities obtained for periodically composites with randomly distributed inclusions are situated between the ones derived by the Mori-Tanaka model (or GSCS) and the differential approximation model. The best estimate proposed in this work for effective thermal con-ductivity of periodically composites with randomly distributed inclusions is the average value of the ones provided by Mori-Tanaka model (or GSCS) and differential approximation model.

Concluding remarks

In this work, the classical Hill-Mendel lemma, which plays a key role in the mechanics and physics of heterogeneous materials, has been extended to incorporating the effects of imperfect interfaces. Three imperfect interface models, namely the highly conducting, lowly conducting and general imperfect thermal models, have been adopted to describe imperfect interfaces. An inclusion embedded in a matrix via an imperfect interface has then been replaced by an equivalent inclusion of the initial shape inserted in the same matrix via the perfect interface. This replacement is based on the requirement that it do not change the total thermal energy. The equivalent inclusion method makes it possible to directly use appropriate micromechanical schemes to determine the effective conductivity of composites with imperfect interfaces. In the present work, the dilute distribution, Mori-Tanaka, self-consistent, generalized selfconsistent and differential approximation schemes have been applied. The validity of the approach proposed in this work has been confirmed through investigating layered composites and circular/spherical particle-reinforced composites made of isotropic constituent phases and by comparing the obtained results with those provided by the FFT method.

Owing to the fact that the thermal conduction phenomenon studied in this work is mathematically similar to other transport phenomena such as electric conduction, dielectrics, magnetism, diffusion and flow in porous media, the approach proposed and the results obtained in the present work are straightforwardly applicable to them. Additionally, it is also fruitful to exploit the mathematical analogy between anti-plane elasticity and 2D thermal conduction.

Fig 11 .

 11 Fig 11. Estimates, bounds and comparaisons with results of Herve and Zaoui[START_REF] Herve | n-Layered inclusion-based micromechanical modelling[END_REF] and Garboczi and Bentz[START_REF] Garboczi | Analytical formulas for interfacial transition zone properties[END_REF] of the normalized effective thermal conductivity, with respect to the conductivity of the matrix phase k 0 , versus the inclusion volume fraction of a tridimensional two-phase composite with general imperfect (GI) interface model with k 1 /k 0 = 10, k (c1) /k 0 = 5 and t 1 /R 1 = 0.001.

Fig 13 .

 13 Fig 13. Estimates and comparaisons with results of Lutz and Zimmerman[START_REF] Lutz | Effect of an inhomogeneous interphase zone on the bulk modulus and conductivity of a particulate composite[END_REF] of the normalized effective thermal conductivity, with respect to the conductivity of the matrix phase k 0 , versus the inclusion volume fraction of a tridimensional two-phase composite with general imperfect (GI) interface model withk 1 /k 0 = 10, D = (k 0 -k (c1)if )/k 0 = 0 and t 1 = 2.3a 1 /ρ 1 = 0.23.

  Fig 15. (a) Typical unit cell of 2D periodic three-phase composites with squarely, hexagonally and randomly distributed inclusions (three-phase model); (b) Typical unit cell of 2D periodic two-phase composites with squarely, hexagonally and randomly distributed equivalent inclusions (equivalent inclusion model)

  = -0.75 and t 1 = 2.3a 1 /ρ 1 = 0.23.

		2.6
			kDD kDD
		3	kSC kSC
		2.4	kMT or kGSCS kMT or kGSCS
			kDA kDA
			Approximation of Lutz and Zimmerman Approximation of Lutz and Zimmerman
		2.2
		2.5
		2
	0 0
	k eff /k k eff /k	2 1.8
		1.6
		1.5 1.4
		1.2
		1 1	0 0	0.1 0.1	0.2 0.2	0.3 0.3	0.4 0.4	0.5 0.5	0.6 0.6
	Inclusion volume fraction Fig 12. Estimates and comparaisons with results of Lutz and Zimmerman [40] of the normalized effective thermal conductivity, with respect to the conductivity of the matrix phase k 0 , versus the inclusion volume fraction of a tridimensional two-phase composite with general imperfect (GI) interface model with k 1 /k 0 = 10, D = (k 0 -k (c1) if )/k 0 Inclusion volume fraction

  Comparison between the effective thermal conductivities obtained by applying the EIM and three-phase model with their estimates and bounds for random and periodic composites whose k 1 /k 0 = 0.1, k (c1) /k 0 = 10 and t 1 /R 1 = 0.001.

		1
		0.9
		0.8
		0.7
	0 /k k eff	0.6 0.5		K eff 11 = K eff 22 First upper bound K eff 11 = K eff 22 First lower bound K eff 11 = K eff 22 DD model K eff 11 = K eff 22 SC model
				K eff 11 = K eff 22 MT model or kGSCS
		0.4		K eff 11 = K eff 22 DA model K eff 11 = K eff 22 Squared lattice with three-phase model
				K eff 11 = K eff 22 Squared lattice with EIM
		0.3		K eff 11 = K eff 22 Hexagonal lattice with three-phase model K eff 11 = K eff 22 Hexagonal lattice with EIM
				K eff 11 Random composite with three-phase model
		0.2		K eff 22 Random composite with three-phase model K eff 11 Random composite with EIM
				K eff 22 Random composite with EIM
		0.1	0	0.05	0.1	0.15	0.2	0.25	0.3	0.35	0.4	0.45	0.5
						Inclusion area fraction
	0 Fig16. Inclusion area fraction 0.05 0.1 0.15 0.2 0.25 0.3 k eff /k 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 K eff 11 = K eff 22 First upper bound K eff 11 = K eff 22 First lower bound K eff 11 = K eff 22 DD model K eff 11 = K eff 22 SC model K eff 11 = K eff 22 MT model or kGSCS K eff 11 = K eff 22 DA model K eff 11 = K eff 22 Squared lattice with three-phase model K eff 11 = K eff 22 Squared lattice with EIM K eff 11 = K eff 22 Hexagonal lattice with three-phase model K eff 11 = K eff 22 Hexagonal lattice with EIM K eff 11 Random composite with three-phase model K eff 22 Random composite with three-phase model K eff 11 Random composite with EIM K eff 22 Random composite with EIM	0.35	0.4	0.45	0.5

(i) K (i) K (n) Γ (n) K (n) Γ (n) Γ(1) Γ (2) Γ (i) Γ

Acknowledgmens

This research is funded by Vietnam National Foundation for Science and Technology Development (NAFOSTED) under grant number 107.02-2017.10