L E Quang

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

Q.-C He

G Bonnet

Eshelby's tensor fields and effective conductivity of composites made of anisotropic phases with Kapitza's interface thermal resistance

Keywords: anisotropy, composite materials, conductivity, interfaces, micromechanics, thermal transport Eshelby problem, Eshelby tensor, Kapitza interface thermal resistance Eshelby problem, Eshelby tensor, Kapitza interface thermal resistance, Interfaces, Anisotropy, Composite materials, Conductivity, Micromechanics, Thermal transports

Eshelby's tensor fields and effective conductivity of composites made of anisotropic phases with Kapitza's interface thermal resistance

Hung Le Quang, Qi-Chang He, Guy Bonnet

Introduction

Most of the studies dedicated to estimating the effective properties of inhomogeneous materials from their local phase properties and their microstructure adopt the assumption that the interfaces between the phases are perfect. In the context of thermal conduction, a perfect interface is a material surface across which the temperature and the normal component of the heat flux are both continuous. In practice, many situations occur where the assumption of perfect interface is inappropriate. Frequently, a thermal resistance appears at the interface between the constituent phases, so that the heat flux remains continuous but the temperature suffers a jump.

The thermal interface resistance in liquid helium was first observed in 1936 by Kurti et al. [27] and Keesom and Keesom [26] and then systematically studied in 1941 by Kapitza [25]. The presence of an interfacial thermal resistance, called hereafter Kapitza's thermal resistance, was also shown experimentally at solid/solid interfaces by many other researchers (see e.g. Swartz and Pohl [43,44], Challis [6], Stoner [40] and Cahill et al [5]). In fact, in physics and materials science, the Kapitza thermal interface resistance has been recognized to be of a great theoretical and practical importance. Applications concern, for example, the development of microelectronic semiconductor devices which need an efficient heat dissipation. Another example is related to the interface of high thermal resistance which can be applied to elements requiring a good thermal isolation such as those in jet engine turbines. To achieve a high thermal resistance, metal-ceramic composites and multilayered systems are currently used for these applications.

The effect of thermal resistance interfaces on the effective conductivity of composites was the subject of many studies (see, e.g., Benveniste and Miloh [4], Benveniste [2], Hasselman and Johnson [21], Lipton and Vernescu [30,31], Torquato and Rintoul [41], Cheng and Torquato [8], Nan et al. [START_REF] Nan | [END_REF], Hashin [18], Dunn and Taya [11], Hatta and Taya [22,23]). However, all of these works are limited to the case where the phases constituting the composite are isotropic.

The present work is concerned with the effective conductivity of composites consisting of a matrix in which circular inclusions in the two-dimensional (2D) case or spherical inclusions in the three-dimensional (3D) case are embedded via interfaces exhibiting Kapitza's thermal resistance. It has the following two objectives: (i) First, it aims at deriving the closed-form expressions for the generalized Eshelby's interior and exterior conduction tensor fields and for the localization tensor fields in the general case where the materials constituting the matrix and the inclusions have the most general thermal anisotropy, i.e., orthotropy. In contrast with the results presented by Hashin [19], Qu [36] and Zhong and Meguid [46] for Eshelby's elastic inhomogeneity problem with the linear spring-layer interface which is the elastic counterpart of Kapitza's thermal resistance interface model, our results show that Eshelby's thermal conduction tensor fields inside the circular or spherical inclusion are uniform. This uniformity makes it possible to transform the circular or spherical inhomogeneity problem into an equivalent circular or spherical inclusion problem and thus to preserve Eshelby's formalism even in the presence of Kapitza's thermal interface resistance. The temperature gradient and heat flux localization tensors can then be obtained and expressed compactly in terms of Eshelby's tensors. (ii) Next, it has the purpose of estimating the effective conductivity of composites under consideration by using the derived expressions of the localization tensor fields for extending some well-established micromechanical schemes to situations where the interface resistance cannot be neglected.

The paper is structured as follows. In Section 2, the constitutive laws of the constituent phases of composites under investigation, the interface model with Kapitza's thermal resistance and the general form of the effective thermal conduction behavior are specified. Section 3 is dedicated to deriving the size-dependent Eshelby's and localization tensor fields in the context of thermal conduction in the presence of Kapitza's interface thermal resistance. In Section 4, closed-form expressions are obtained for the effective conductivity moduli by using the dilute, Mori-Tanaka, self-consistent and generalized self-consistent schemes, respectively. In Section 5, the Kapitza thermal interface resistance and inhomogeneities size effect on the effective conductivity of composites are numerically discussed and illustrated; in addition, the closed-form solutions obtained are compared with the results provided by the finite element method (FEM). In section 6, a few concluding remarks are provided.

Local constitutive laws

The composite under investigation consists of circular or spherical inclusions embedded in a matrix via imperfect interfaces. Let Ω be the domain of a 2D or 3D Euclidean space R d (d = 2, 3) occupied by a representative volume element (RVE) of the composite and let ∂Ω be the boundary of Ω. The subdomains of Ω inhabited by the inclusions and the matrix are denoted by Ω (1) and Ω (2) , respectively. The interface between Ω (1) and Ω (2) and its unit vector directed from Ω (1) toward Ω (2) are designed by Γ and n (see Fig. 1b). In a system of Cartesian coordinates {x i } associated to a right-handed orthonormal basis {f i } (1 ≤ i ≤ d), the matrix, referred to as phase 2, and inclusions, denoted as phase 1, are assumed to be individually homogeneous and have the linear thermal conduction behavior described by an anisotropic Fourier's law

q (i) = K (i) • e (i) or e (i) = H (i) • q (i) , (1)
where K (i) and H (i) , with i = 1 or 2, stand for the thermal conductivity and resistivity second order tensors of phase i, which are symmetric, positive definite and in general orthotropic. The intensity field e (i) is related to the temperature field T (i) by

e (i) = -∇T (i) . (2)
In the case of stationary thermal conduction without heat source, the heat flux field q (i) must fulfill the energy conservation equation

∇ • q (i) = 0. (3)
In addition, the interface Γ between the matrix and inclusions is described by an imperfect interface model with Kapitza's thermal resistance. According to the latter, the normal component of the heat flux q(x) is continuous across Γ, i.e. holding on Γ:

q (1) (x) • n(x) = q (2) (x) • n(x)
T (2) -T (1) = -α(n)q (1) • n = -α(n)q (2) • n, (4)
β(n)(T (2) -T (1) ) = -q (1) • n = -q (2) • n, (5)
where α and β = 1/α stand for the Kapitza thermal resistance and joint conductivity of Γ, respectively, and depend in general on n.

Next, to clarify the physical background and validity domain of the interface model with Kapitza's thermal resistance used here to describe Γ, we consider two configurations in Fig. 1. In the three-phase one (Fig. 1a), a representative volume element Ω consists of a finite number of inclusions Ω (1) embedded in the matrix phase Ω (2) via the interphase Ω (s) . The interface between Ω (1) and Ω (s) and the one between Ω (2) and Ω (s) are assumed to be perfect, and the thickness h of Ω (s) is uniform and very small in comparison with the minimum size of the inclusions. In the two-phase configuration (Fig. 1b), Ω (s) is replaced by a thermally resistive surface of zero thickness located at the middle surface Γ of the interphase, and the inclusions and matrix are extended up to the middle surface Γ. By using asymptotic expansions, the jump conditions that the imperfect interface Γ has to verify for the two configurations to be physically equivalent between themselves within an error of order 0(h) were derived first by Sanchez-Palencia [37] and Pham Huy and Sanchez-Palencia [35] in a particular case and then completed by Miloh and Benveniste [32], Hashin [18] and Benveniste [3] in the general situation. More precisely, when the interphase with thermal conductivity tensor K (s) is assumed to be weakly conducting, namely K (s)

K (1) and K (s) K (2) , then the conditions that the imperfect interface Γ must satisfy are exactly those characterizing the imperfect interface model with Kapitza's thermal resistance described above. The corresponding Kapitza thermal resistance of Γ can be expressed in terms of K (s) and h as follows (see Benveniste [3]):

α(n) = h n • K (s) • n . (6)
In what follows, we make the assumption that K (s) is circularly or spherically orthotropic (see e.g. Dryden [10]; Chen [7]), so that

k s = n • K (s) • n is constant.
Then, it is immediate from Eq. ( 6) that

α = h k s (7)
is independent of n.

Finally, at the macroscopic scale, the composite under consideration is assumed to be statistically homogeneous. The corresponding effective thermal behavior is characterized by

Q = K ef f • E or E = H ef f • Q (8)
where K ef f and H ef f are the effective thermal conductivity and resistivity secondorder tensors, Q and E denote the macroscopic heat flux and intensity fields. More precisely, the macroscopic intensity field E and the heat flux field Q are defined as

E = - 1 |Ω| ∂Ω T (x)ν(x)dx, (9)
Q = 1 |Ω| ∂Ω (q • ν)xdx, (10)
where ν(x) is the outward unit normal vector to ∂Ω and |Ω| denotes the volume or area of the domain Ω according to the case (3D or 2D) which is concerned. As in the classical case of perfect interfaces, the macroscopic heat flux field Q defined by ( 9) is equal to the volume or surface average of the local counterpart q(x) over Ω, i.e.,

Q = 1 |Ω| Ω q(x)dx = q , ( 11
)
where • is the volume or surface average of quantity • over the domain Ω. However, due to the fact that the temperature field suffers a jump across Γ, the macroscopic intensity field E defined by ( 9) is not simply the volume or surface average of the local counterpart e(x) over Ω but given by

E = e - 1 |Ω| Γ [T (2) (x) -T (1) (x)]n(x)dx, (12)
or equivalently by

E = e + 1 |Ω| Γ α[q (1) (x) • n(x)]n(x)dx. (13)
Thus, compared with the classical case, new terms due to the Kapitza interface thermal resistance appear in both Eqs (12) and (13).

3. Generalized Eshelby's formalism

Integral formulation

Now, let V be a subdomain of a 2D or 3D domain Ω made of a homogeneous material with thermal conductivity and resistivity tensors denoted by K 0 and H 0 , respectively. The boundary of V and its unit outward normal vector are denoted by ∂V and n. The interface Γ between the subdomain V and the remaining part of Ω is not necessarily perfect. If a heat-free (or an eigen) intensity field e * (x) is prescribed over V , the induced heat flux field q(x) is related to the intensity field e(x) produced by e * (x) through Fourier's law

q(x) = K 0 • e(x) + q * (x) (14)
with x ∈ Ω \ Γ and

q * (x) = -K 0 • e * (x). ( 15
)
With the help of the thermal result analogous to the Maxwell-Betti reciprocal theorem of elasticity, we can derive an integral representation of the temperature field T (x) due to the prescribed heat-free intensity field e * (x) as follows:

K 0 ij ∂V [T (y)G ,j (y -x) -T ,j (y)G(y -x)]n i (y)dy -K 0 ij V e * j,i (y)G(y -x)dy = T (x) for x ∈ V , 0 for x ∈ Ω \ V . (16)
Hereafter, the subscript i following a comma denotes the derivative with respect to y i ; G(yx) represents the Green's function giving the temperature at a point x produced by a unit point heat source at a point y in an infinite domain. More precisely, G(yx) is the solution to the equation

K 0 ij G ,ij (y -x) = δ(y -x), ( 17
)
where δ is the Dirac delta function.

Generalized Eshelby's conduction tensor fields

In this subsection, we are interested in the thermal counterpart of the well-known Eshelby's elastic inclusion problem [12,[START_REF] Eshelby | Elastic inclusion and inhomogeneities[END_REF] while accounting for imperfect interface effects. Indeed, the domain Ω described above in the last subsection is now assumed to be infinitely extended and the heat-free intensity e * (x), which plays a role similar to the eigenstrain tensor in elasticity, is distributed within a subdomain ω of Ω and vanishes outside ω (Fig. 2). The interface Γ between the inclusion ω and the external medium has the outward normal vector n and is assumed to exhibit the Kapitza thermal resistance. When the domain V is taken to be the inclusion ω, applying the integral formulation (16) and using the divergence theorem yield

T (x) = K 0 ij Γ - {T (y)G ,j (y -x) -G(y -x)[T ,j (y) + e * j (y)]}n i (y)dy + K 0 ij ω e * j (y)G ,i (y -x)dy, for x ∈ ω. (18)
Here Γ -denotes the inside part of the interface between the inclusion and the external medium, i.e. the limit of a surface inside ω tending to Γ. For the points outside the inclusion ω, it follows from Eq. ( 16) together with the divergence theorem that

0 = -K 0 ij Γ - {T (y)G ,j (y -x) + G(y -x)[T ,j (y) + e * j (y)]}n i (y)dy -K 0 ij ω e * j (y)G ,i (y -x)dy, for x ∈ Ω \ ω. (19)
Next, if V is taken to be the exterior domain of ω, i.e. V = Ω\ω, from the integral formulation (16) and by taking into account the fact that the Green function G(yx) and the temperature field T (y) are equal to zero for y ∈ ∂Ω, we can show that

0 = K 0 ij Γ + [T (y)G ,j (y -x) -G(y -x)T ,j (y)] n i (y)dy (20)
for x ∈ ω. Here, Γ + is the outside part of the interface Γ as opposed to Γ -. For the points exterior to the inclusion ω, we derive from Eq. ( 16) the temperature field outside the inclusion as follows:

T (x) = -K 0 ij Γ +
[T (y)G ,j (yx) -G(yx)T ,j (y)] n i (y)dy (21) for x ∈ Ω \ ω.

Due to the fact that the normal component of the heat flux field q(x) is continuous across Γ, subtracting Eqs. (20) and (19) from Eqs. (18) and (21), respectively, leads to the temperature field both inside and outside the inclusion

T (x) = -K 0 ij Γ [[T (y)]]G ,j (y -x)n i (y)dy + K 0 ij ω e * j (y)G ,i (y -x)dy, (22)
where [[T (y)]] = T (y + ) -T (y -) with y + ∈ Γ + and y -∈ Γ -. Correspondingly, the resulting intensity field components both inside and outside the inclusion are given by

e m (x) = - ∂T (x) ∂x m = -K 0 ij Γ [[T (y)]]G ,jm (y -x)n i (y)dy + K 0 ij ω e * j (y)G ,im (y -x)dy. (23)
Moreover, by using the interface condition Eq.( 4), the intensity field components

Now, let us show that if the heat-free intensity e * is uniformly distributed within a spherical or circular inclusion, i.e. e * (x) = χ ω (x)e 0 with the characteristic function χ ω (x) of ω defined by

χ ω (x) = 1 for x ∈ ω, 0 for x / ∈ ω, (25)
then the intensity field given by Eq.( 24) inside the inclusion ω is, as in the case of perfect interfaces, independent of the position x. To this end, we first compute the right members of Eq. ( 22) and Eq. ( 24) for a constant intensity field e(x) = e ω within the inclusion.

In the important special case where ω is a spherical or circular inclusion of radius R, substituting n i (y) = y i /R for any y ∈ Γ into Eq.( 24), using the divergence theorem to transform the second integral of Eq. ( 24) and accounting for the fact that

∂G(y -x) ∂x i = - ∂G(y -x) ∂y i , (26)
the resulting intensity components are given by

e m (x) = K 0 ij e 0 j ω G ,im (y -x)dy + α R K 0 ij K 0 sr M ijms (x)(e ω r -e 0 r ), (27)
for x ∈ ω. Here,

M ijms (x) = - ∂N ijs (x) ∂x m + δ is D mj (x), (28)
with

N ijs (x) = ∂ 2 ∂x i ∂x j ω G(y -x)y s dy, (29)
D mj (x) = ∂ 2 ∂x m ∂x j ω G(y -x)dy. (30)
Next, it is shown in the appendix that when ω is an ellipsoidal or elliptical inclusion and when x ∈ ω, then N mji (x) and D mj (x) are polynomial functions of degree 1 and 0 in terms of the position vector x. Thus, the tensor M(x) with the components M ijms (x) given by Eq.( 28) is uniform inside ω and equal to a constant tensor M ω , i.e. M(x) = M ω when x ∈ ω. As a consequence, the intensity field is constant within the inclusion and the components e m are identical to e ω m for any x inside ω. Thus, for any x inside ω, we

(x) = K 0 ij e 0 j ω G ,im (y -x)dy + α R K 0 ij K 0 sr M ω ijms (e ω r -e 0 r ). (31)
This last equation can be recast into

e ω = Sω • e 0 (32)
with

Sω = [I - α R K 0 : M ω • K 0 ] -1 • [S ω - α R K 0 : M ω • K 0 ]. (33)
Here, S ω stands for the classical Eshelby's conduction tensor inside ω without interface effect while Sω is defined as the Eshelby's conduction tensor inside ω with Kapitza interface thermal resistance; I is the second-order identity tensor. Indeed, the classical Eshelby's conduction tensor fields both inside and outside ω with perfect interface are, in general, given by (see, e.g. Le Quang et al. [28])

S(x) = ω ∇ ⊗ ∇G(y -x)dy • K 0 . ( 34
)
Using the definition (30), the classical Eshelby's conduction tensor field S(x) can be determined via the tensor D(x):

S ij (x) = D im (x)K 0 mj . (35)
Thus, the general solution to equation ( 24), both inside and outside ω, can be written as

e(x) = S(x) • e 0 (36)
where the tensor S(x), called thereafter Eshelby's conduction tensor field with Kapitza interface thermal resistance, is defined as:

S(x) = S(x) + α R K 0 : M(x) • K 0 • [ Sω -I], (37)
where Sω is given by Eq. [START_REF] Mura | Mircomechanics of Defects in Solids[END_REF]. Eq. (37) shows that the Eshelby's conduction tensor S(x) with Kapitza's interface thermal resistance is uniform inside the inclusion ω but varies outside the inclusion ω. Note that S(x) depends not only on the geometry of ω and the material thermal anisotropy as in the case of perfect interfaces but also on the interfacial thermal resistance parameter α and the inclusion size R. When α = 0 or when the size R of the heterogeneity is very large, the imperfect interface effect disappears or becomes negligible so that the expression (37) of S(x) reduces to the formula [START_REF] Nan | [END_REF] of the classical Eshelby's conduction tensor.

The elastic counterpart of Kapitza's interface thermal resistance model is the well-known spring-layer imperfect interface model. When the latter is adopted in studying a spherical inclusion embedded in an isotropic matrix, Hashin [19], Qu [36] and Zhong and Meguid [46] have shown that the Eshelby's elastic tensor S(x) is in general a quadratic function of the position vector, so that it is not uniform inside the spherical inclusion. The results derived above for Eshelby's thermal conduction tensor S(x) with Kapitza's interface thermal resistance hold for any thermal anisotropy. In order to obtain the explicit expression of S(x), we now consider the most important special case where the 2D or 3D infinite body Ω consists of an isotropic material whose thermal conductivity tensor takes the form K 0 = k 0 I with k 0 being a positive scalar.

By using (35) with D im (x) given by Eqs. (A-43) and (A-44) for the 3D case or by Eqs. (A-55) and (A-56) for the 2D case, the classical Eshelby's conduction tensor field reads

S(x) = 1 d Iχ ω (x) + (1 -χ ω (x))ρ d 1 d I -x ⊗ x ( 38
)
where

ρ = R r = R |x| , x = x r . (39)
By substituting Eq. ( 38) into Eq. ( 37) with D mj (x) and N mij (x) determined in the appendix, Eshelby's conduction tensor field with Kapitza's interface thermal resistance in the case of an isotropic medium is obtained as

S(x) = 1 + (d -1)α d + (d -1)α Iχ ω (x) + (1 -χ ω (x)) dρ d d + (d -1)α 1 d I -x ⊗ x , (40)
where the dimensionless parameter α is defined by

α = αk 0 R . ( 41
)
From a dual point of view, we now consider the conjugate Eshelby's conduction problem where a heat-free flux is uniformly distributed within a spherical or circular inclusion, or equivalently q * (x) = χ ω (x)q 0 (see Fig. 3). The general heat flux field solution, both inside and outside the inclusion ω, can be determined via the conjugate Eshelby's conduction tensor field with Kapitza's interface thermal resistance, denoted by C(x), such as

q(x) = C(x) • q 0 . ( 42
)
It can be shown that C(x) is related to S(x) by

C(x) = -K 0 • S(x) • H 0 + χ ω (x)I. ( 43
)
This equation indicates that, as expected, the conjugate Eshelby's conduction tensor C(x) with Kapitza's interface thermal resistance is also uniform inside ω, i.e. C(x) = Cω for x ∈ ω with

Cω = -K 0 • Sω • H 0 + I. (44)
In particular, when the interface between the inclusion and the exterior medium is perfect, i.e. α = 0 or when the size R of the inclusion ω is very large, we obtain from Eq. ( 43) the classical conjugate Eshelby's conduction tensor

C(x) = -K 0 • S(x) • H 0 + χ ω (x)I (45)
and the classical conjugate Eshelby's conduction tensor inside the inclusion reads

C ω = -K 0 • S ω • H 0 + I. (46)
In the particular case where the 2D or 3D infinite body Ω is made of an isotropic material, by substituting the expressions of S(x) and S(x) given by Eqs. (40) and (38) into Eqs. ( 43) and (45), respectively, we obtain the expressions for the conjugate Eshelby's conduction tensors with and without Kapitza's thermal resistance effect:

C(x) = d -1 d + (d -1)α Iχ ω (x) -(1 -χ ω (x)) dρ d d + (d -1)α 1 d I -x ⊗ x , (47)
C(x) = d -1 d Iχ ω (x) -(1 -χ ω (x))ρ d 1 d I -x ⊗ x . ( 48
)

Generalized localization tensor fields

Now, consider a 2D or 3D infinitely extended matrix Ω with the thermal conductivity and resistivity tensors, K 0 and H 0 , in which an inclusion (or more precisely an inhomogeneity) ω made of another material with the thermal conductivity and resistivity tensors denoted by K (1) and H (1) , is embedded. As before, the interface Γ between matrix and inclusion is assumed to exhibit the Kapitza thermal resistance. On the boundary ∂Ω of Ω, the following homogeneous heat flux boundary condition is prescribed

q(x) • ν(x) = q 0 • ν(x), (49)
where x ∈ ∂Ω, ν(x) is the outward unit normal to ∂Ω and q 0 is a constant heat flux field.

Similarly to the results obtained in the previous subsection for the Eshelby's conduction tensor with Kapitza interface thermal resistance, we proceed to show that the heat flux and intensity fields are uniform within the spherical or circular inhomogeneity i.e. q(x) = q ω and e(x) = H (1) • q ω for any x ∈ ω. For this, we decompose the inhomogeneity problem under consideration into three Sub-problems, denoted as A, B and C (see Fig. 4). Each Sub-problem and its solution are presented in detail as follows. heat flux solution field q (A) (x) is provided directly by

q (A) (x) = q 0 (50) for any x ∈ Ω. b) Sub -problem B
Sub-problem B consists in computing the heat flux field within the infinite heterogeneous matrix medium Ω with the thermal conductivity and resistivity tensors, K 0 and H 0 , in which the sub-domain ω made of the same material as the matrix and surrounded by the interface Γ with Kapitza's thermal resistance is subjected to a uniform heat-free flux field q * = ( Cω ) -1 • q ω . Applying the conjugate Eshelby's conductivity tensor fields obtained in Section 3, the heat flux field is provided by

q (B) (x) = C(x) • ( Cω ) -1 • q ω (51)
where

x ∈ Ω. c) Sub -problem C
In this problem, the sub-domain ω embedded in the infinite homogeneous matrix medium Ω is formed of the same material as the matrix, with the thermal conductivity and resistivity tensors, K 0 and H 0 , and the interface between ω and the external medium is assumed to be perfect. Then, we subject ω to uniform heat-free flux field

q * = -K 0 • (H (1) -H 0 ) • q ω -( Cω ) -1 • q ω .
Due to the fact that the interface is perfect, the heat flux solution field q (C) (x) can be determined via the classical conjugate Eshelby's conductivity tensor field by

q (C) (x) = -C(x) • K 0 • (H (1) -H 0 ) • q ω -C(x) • ( Cω ) -1 • q ω (52) with x ∈ Ω.
Finally, by applying the superposition principle, the heat flux field of the inhomogeneity problem is given by

q(x) = q (A) (x) + q (B) (x) + q (C) (x) = q 0 + C(x) • ( Cω ) -1 • q ω -C(x) • K 0 • (H (1) -H 0 ) • q ω -C(x) • ( Cω ) -1 • q ω (53)
where x ∈ Ω.

When x ∈ ω, we have q(x) = q ω . Introduction of this relation into Eq. (53) allows us to calculate q ω as follows:

q ω = Bω • q 0 , ( 54
)
where Bω , called thereafter the localization tensor for the inhomogeneity ω, has the expression

Bω = [C ω • ( Cω ) -1 + C ω • K 0 • (H (1) -H 0 )] -1 , (55)
with C ω and Cω given by Eqs. ( 46) and (44), respectively. After obtaining the expression of q ω , the heat flux field q(x) in Eq. ( 53) for both inside and outside ω can be calculated by combining equations ( 54) and ( 55) with (53), so that

q(x) = B(x) • q 0 (56)
where B(x), named the localization tensor field with Kapitza's interface thermal resistance, is given by:

B(x) = I -{C(x) • K 0 • (H (1) -H 0 ) + [C(x) -C(x)] • ( Cω ) -1 } • Bω (57)
where Bω , C(x), C(x) and Cω are provided by Eqs. ( 55), ( 45), ( 43) and (44).

In the important particular case where the matrix and inclusion phases are both isotropic, their thermal conductivity and resistivity tensors are given by K 0 = k 0 I, H 0 = 1/k 0 I, K (1) = k 1 I and H (1) = 1/k 1 I with k 0 and k 1 being two positive scalars. Then, the localization tensor field B(x) has the simple expression

B(x) = dk 1 k 1 + (d -1)(k 0 + k 1 α) Iχ ω (x) + (1 -χ ω (x)) I + dρ d [k 0 + k 1 (α -1)] k 1 + (d -1)(k 0 + k 1 α) 1 d I (d) -x ⊗ x ( 58
)
where α is defined by Eq. ( 41). Remark that, if α = 0, B(x) reduces to the wellknown classical localization tensor field B(x).

Effective thermal behavior

Starting from the results obtained in the previous two sections, we will derive in this section the closed-form expressions for the effective conductivity of composites accounting for the effects of the Kapitza interface thermal resistance between the matrix and inclusion phases following the dilute, Mori-Tanaka, self-consistent and generalized self-consistent schemes. For more details about these schemes in the case of perfect interfaces, the reader can refer to the review papers of Hashin [15,17].

Dilute distribution model

Assuming that the interaction between the inhomogeneities of the composite described in Section 2 is neglected, a dilute distribution of inhomogeneities is adopted to estimate the effective thermal conductivity or resistivity tensors, K ef f or H ef f , respectively. For this purpose, we first let Ω be subjected to the homogeneous heat flux boundary conditions given by Eq. ( 49). Under boundary condition (49) and owing to the fact that the normal component of heat flux is continuous across the interface between the matrix and inclusions phases, the macroscopic heat flux defined by Eq. ( 10) is obtained by

Q = q = q 0 . ( 59
)
Within the framework of the dilute distribution model, the interaction between the inhomogeneities is omitted and, accounting for the results obtained for the

(1) = [C ω • ( Cω ) -1 + C ω • K (2) • (H (1) -H (2) )] -1 • q 0 , (60)
ē(1) = H (1) • q(1) .

Here, K (i) and H (i) are thermal conductivity and resistivity second-order tensors of phase i, respectively. The conjugate Eshelby's conduction tensors without and with interface effects, C ω and Cω , inside ω are respectively given by Eqs. ( 46) and ( 44) taking into account Eqs. ( 33) and ( 34) with K 0 = K (2) and H 0 = H (2) . From Eq. ( 59), the volume or surface averages of the heat flux and intensity fields in the matrix phase are given by

q(2) = c -1 2 (q 0 -c 1 q(1) ), (62)
ē(2) = H (2) • q(2) , (63)
where c 1 and c 2 = 1c 1 are the volume or surface fractions of the inclusion and matrix phases, respectively. The macroscopic intensity field can be obtained from Eq. ( 13) and is given by

E = c 1 ē(1) + c 2 ē(2) + c 1 α R q(1) . (64)
Unlike the classical case where the matrix-inclusion interface is perfect, the macroscopic intensity field given by Eq. ( 64) is not simply the volume or surface average of the different contributions over the domain Ω. More precisely, a new term due to the discontinuity of both temperature field and normal component of the intensity field across the Kapitza interface Γ is added in Eq. (64). By introducing Eqs.

(60)-(63) into Eq. ( 64) and by identifying the resulting expression with Eq. ( 8), we obtain the expression of the effective thermal conductivity tensor

K DD = {H (2) + c 1 [H (1) -H (2) + α R I] • [C ω • ( Cω ) -1 + C ω • K (2) • (H (1) -H (2) )] -1 } -1 . (65)
It is interesting to note that when α = 0, the expression (65) of K DD reduces to the well-known formula without interface effects. In the particular case where the matrix and inclusion phases are isotropic and have the thermal conductivities k 2 and k 1 , respectively, Eq. (65) leads to

k DD = k 2 1 + c 1 d[k 2 -k 1 + αk 1 k 2 /R] k 1 + (d -1)k 2 (1 + αk 1 /R) -1 (66)
where k DD is the effective thermal conductivity of the composite for the dilute distribution model. The expression (66) for the effective thermal conductivity allows us to recover the formula derived by Benveniste and Miloh [4].

Mori-Tanaka model

To account for the interaction between the heterogeneities, the Mori-Tanaka model is now applied. In this model, under the boundary conditions written as in Eq. ( 49), the heat flux field in the matrix and far from the inclusions is assumed to be uniform and equal to Q 0 . Consequently, the heat flux and intensity fields inside the inhomogeneities ω are obtained by using Eqs. ( 60) and (61) with q 0 being replaced by Q 0 . Thus, this yields

q(1) = [C ω • ( Cω ) -1 + C ω • K (2) • (H (1) -H (2) )] -1 • Q 0 , (67)
ē(1) = H (1) • q(1) ,

where the conjugate Eshelby's conduction tensors without and with interface effects, C ω and Cω , inside ω are respectively provided by Eqs. ( 46) and ( 44), taking into account Eqs. ( 33) and ( 34) with K 0 = K (2) and H 0 = H (2) . Concerning the matrix phase, the volume or surface average intensity and heat flux fields are given by

q(2) = Q 0 , ē(2) = H (2) • Q 0 . ( 69
)
Due to the fact that Q = q = q 0 , we can express Q 0 in terms of the macroscopic heat flux field as follows

Q 0 = L • q 0 . (70)
with

L = {c 1 [C ω • ( Cω ) -1 + C ω • K (2) • (H (1) -H (2) )] -1 + c 2 I} -1 . (71)
By substituting Eqs. (67)-(71) into Eq. ( 64) and by identifying the resulting equation with Eq. ( 8), the effective thermal conductivity tensor is given by

K M T = {H (1) + α R I -c 2 [H (1) -H (2) + α R I] • L} -1 . (72)
In the special case where the matrix and inclusion phases are isotropic and have the thermal conductivities k 2 and k 1 , respectively, we obtain from Eq. (72) the effective thermal conductivity of the composite for the Mori-Tanaka model

k M T = k 2 1 + c 1 d[k 2 -k 1 + αk 1 k 2 /R] k 1 c 1 d + c 2 [(d -1)k 2 (1 + αk 1 /R) + k 1 ] -1 . ( 73
)
We notice that this expression for the effective thermal conductivity is identical to the one given by Benveniste [2]. Moreover, it is easy to check that when the volume or surface fraction of the inclusion is small enough (c 1 1), Eq. (73) giving the effective thermal conductivity of the composite for the Mori-Tanaka model reduces to the one obtained by Eq. (66) for the dilute distribution model. As for all estimate schemes for the effective thermal behavior of the composite, it is necessary to determine the volume or surface average of the intensity and heat flux fields inside the inhomogeneities. In the self-consistent model, the intensity and heat flux fields inside each inhomogeneity are estimated by embedding this inhomogeneity in a fictitious infinite medium which has unknown effective properties. By using the expression for the solution of the localization problem for thermal conduction described in Section 4, we can calculate the heat flux and intensity fields inside the inhomogeneity ω as follows:

q(1) = [C ω • ( Cω ) -1 + C ω • K SC • (H (1) -H SC )] -1 • q 0 , (74)
ē(1) = H (1) • q(1) .

In these expressions, C ω and Cω are given by Eqs. ( 46) and ( 44) together with Eqs. ( 33) and [START_REF] Nan | [END_REF] in which K 0 and H 0 are replaced by the unknown effective thermal conductivity and resistivity tensors, K SC and H SC , respectively. As in the dilute distribution model, owing to the fact that Q = q = q 0 , the volume or surface averages of the heat flux and intensity fields in the matrix phase are given by q

(2) = c -1 2 (q 0 -c 1 q(1) ), (76) ē
(2) = H (2) • q(2) . (77)
Introducing Eqs. ( 74)-(77) into Eq. ( 64) and taking into account the macroscopic behavior expressed by Eq. ( 8), the effective thermal conductivity tensor is obtained from a quadratic equation in K SC which can be written as

K SC = {H (2) + c 1 [H (1) -H (2) + α R I] • [C ω • ( Cω ) -1 + C ω • K SC • (H (1) -H SC )] -1 } -1 . (78)
This second order polynomial matrix equation provides the effective conductivity tensor K SC . If the matrix and inclusion phases are isotropic and characterized by the thermal conductivities k 2 and k 1 , Eq. ( 78) is reduced to a quadratic equation

k SC = k 2 1 + c 1 d[k 2 -k 1 + αk 1 k 2 /R] k 1 + (d -1)k SC (1 + αk 1 /R) -1 , (79)
and the effective thermal conductivity k SC is given by the real positive root of Eq. (79).

Generalized self-consistent model

The model proposed in this section does not use the Eshelby's tensor fields and localization tensor fields but is presented here to provide results allowing a comparison with those coming from the three models described above. It can be considered as an extension to the thermal conduction problem with interface effects of the classical GSCM for elasticity which was initiated by Van der Poel [42], improved and completed by Smith [38,39] and Christensen and Lo [9]. According to this model, we first consider an infinite three-dimensional medium M consisting of the effective homogeneous and isotropic medium whose thermal behavior is characterized by Eq. ( 8) in which the effective thermal conductivity tensor is denoted by K ef f = k GSCM I (3) . As before, let M be subjected to the uniform boundary conditions Eq. ( 49) in which the constant heat flux vector q 0 is chosen such as q 0 = (0, 0, q 0 ) T where q 0 is a constant heat flux field (Fig. 2). In the system of spherical coordinates (r, θ, ϕ) corresponding to the spherical orthogonal basis (f r , f θ , f ϕ ), this uniform boundary condition takes the equivalent form

q 0 (x) = q 0 (cos θf r -sin θf θ ), x ∈ ∂M. ( 80
)
This boundary condition gives rise to the following temperature, intensity and heat flux fields at point x of M T 0 (x) = -q 0 k GSCM r cos θ, e 0 (x) = q 0 k GSCM (cos θf rsin θf θ ), q 0 (x) = q 0 (cos θf rsin θf θ ).

(81)

Let us introduce the following energy U 0 (q 0 ) of M (see e.g. Hashin [16])

U 0 (q 0 ) = M q 0 (x) • e 0 (x)dx = vol(M )(q 0 ) 2 /k GSCM . ( 82
)
Next, we cut a sphere out of the foregoing infinite effective medium and substitute into its place a composite sphere Ω, while applying the same boundary condition on ∂M as before. The interface between the composite sphere and the outside medium is assumed to be perfect. The core of this composite sphere is made of the inclusion phase, referred to as phase 1 and surrounded by a concentric shell consisting of the matrix phase, denoted by phase 2. The core and outer coating consist of two isotropic materials whose thermal conductivities are k 1 and k 2 . The radii of the core and coating, symbolized by r 1 and r 2 are chosen so as to be compatible with the prescribed phase volume fraction

c 1 = 1 -c 2 = r 3 1 r 3 2 . ( 83
)
The 3D spherical interface Γ between the matrix and inclusion is modeled by the Kapitza interface thermal resistance as described in Section 2. Under boundary condition (80), the expressions of the temperature field, nonzero intensity and heat flux field components are given by (see, e.g. Hervé [24])

T (i) = -a i r + b i r 2 cos θ, (84)
e (i) r = a i - 2b i r 3 cos θ, e (i) θ = -a i + b i r 3 sin θ, (85)
q (i) r = k i a i - 2b i r 3 cos θ, q (i) θ = -k i a i + b i r 3 sin θ. ( 86
)
with i = 1 referring to the core, i = 2 to the outer coating and i = e to the external effective medium. In these expressions, a i and b i are constants to be determined from the boundary and interface conditions together with a condition avoiding the displacement singularity in the core of the composite sphere. More precisely, the requirement of the temperature in the core at r = 0 to be finite implies that b 1 = 0. The value of a e , determined by using the boundary condition (80) with r → ∞, is given by a e = q 0 /k GSCM . At the interface between the core and the coating, the normal component of heat flux field is continuous, so that

k 2 a 2 - 2b 2 r 3 1 = k 1 a 1 . ( 87
)
The surface condition (4) at the interface between the core and the coating can be specified by :

a 2 + b 2 r 3 1 -a 1 = α a 2 - 2b 2 r 3 1 , (88)
with α = αk 2 /r 1 . The interface at r = r 2 between the coating matrix and outside effective medium is perfectly bonded. Thus, the continuity conditions of both temperature field T and normal component heat flux field q r across the interface at r = r 2 are expressed as

q 0 k GSCM + b e r 3 2 = a 2 + b 2 r 3 2 , ( 89
)
k GSCM q 0 k GSCM - 2b e r 3 2 = k 2 a 2 - 2b 2 r 3 2 . ( 90
)
As in the GSCM of Christensen and Lo [9] in the context of elasticity, the effective thermal conductivity is required to be such that the energy U (q 0 ) after introducing the sphere composite is equal to the initial one U 0 (q 0 ) which is given by Eq. ( 82).

On the other hand, we can show that the actual energy U (q 0 ) can be expressed in terms of U 0 (q 0 ) as follows

U = U 0 - Γ (q 0 T (e) -q (e) T 0 ) • ndx. (91)
Thus, the self-consistency condition U 0 = U is reduced to

Γ (q 0 T (e) -q (e) T 0 ) • ndx = 0. ( 92
)
Substituting Eqs. ( 80), ( 81), ( 84) and ( 84) with i = e into Eq. ( 92), we obtain the simple equation

b e = 0. ( 93
)
Finally, substituting b e = 0 into Eqs. ( 87), ( 88), ( 89) and ( 90), we obtain a system of four homogeneous linear equations for the four unknowns a 1 , a 2 , b 2 and q 0 . A non-trivial solution to this system exists if and only if the determinant of the relevant 4 × 4 matrix is equal to zero. This necessary and sufficient condition yields the expression for the effective thermal conductivity as follows:

k GSCM = k 2 1 + 3c 1 [k 2 + (α -1)k 1 ] 3k 1 c 1 + c 2 [2(k 2 + αk 1 ) + k 1 ] -1 . ( 94
)
We remark that the effective thermal conductivity obtained by Eq. ( 94) depends not only on the phase thermal properties and volume fractions but also on the interface thermal conductivity and on the size of the inhomogeneities through α.

When setting α = 0, the expression (94) of k GSCM reduces to the one obtained by Hashin and Shtrikman [20] or Hervé [24] without interface effects.

In the 2D case, the corresponding effective thermal conductivity is given by

k GSCM = k 2 1 + 2c 1 [k 2 + (α -1)k 1 ] 2k 1 c 1 + c 2 [k 2 + (α + 1)k 1 ] -1 . ( 95
)
The details of the derivation are omitted here. Finally, it is convenient to express (94) and ( 95) for k GSCM by the following compact one :

k GSCM = k 2 1 + c 1 d[k 2 + (α -1)k 1 ] k 1 c 1 d + c 2 [(d -1)(k 2 + αk 1 ) + k 1 ] -1 , (96)
where d (= 3 or 2) is the dimension of the problem. Compared Eq. (96) with Eq. (73) for the case where the matrix and inclusion phases are isotropic, we notice that k GSCM = k M T . This remark has been also mentioned in the work of Benveniste ( [2]) for the tridimensional case.

Numerical examples

This section is dedicated to numerically illustrate the features of the results obtained above. The first example is related to the plane problem of a circular inhomogeneity having the radius R = 1µm located within a 2D host matrix phase. The matrix and inclusion phases are both assumed to be anisotropic. Additionally, the thermal conductivity tensors, K (1) [W m -1 K -1 ] of the inclusion and K (2) [W m -1 K -1 ] of the matrix, are provided by :

K (1) = 1.5 -0.5 -0.5 1.5 , K (2)
= 0.15 0.05 0.05 0.15 (97) relative to the orthonormal basis {f 1 , f 2 } and by

K (1) = 1 0 0 2 , K (2) = 0.2 0 0 0.1 (98)
according to the orthonormal basis

{f 1 = 1 √ 2 (f 1 + f 2 ), f 2 = 1 √ 2 (f 2 -f 1 )
} which is obtained from the orthonormal basis {f 1 , f 2 } by an anti-clockwise rotation of angle π/4.

Concerning the thermal parameters of the interface, we choose different values of the Kapitza interface thermal resistance α ranging from 10 -11 m 2 K/W (nearly perfect interface) to 10 -3 m 2 K/W (nearly insulating interface).

Next, by applying a uniform heat flux vector q 0 = (1, 0) T [W/m 2 ] at the boundary of the heterogeneous medium consisting of a circular inclusion with radius R = 1µm embedded in a square matrix domain 10µm × 10µm and 20µm × 20µm and by assuming the interface between the matrix and inclusion to be imperfect with Kapitza thermal resistance α = 10 -3 , 10 -5 , 10 -6 and 10 -11 m 2 K/W , we calculate first the heat flux field in the inclusion and matrix phases from the closed-form solution of section 3 corresponding to an infinite matrix. Then, these values obtained for the heat flux field are compared with the ones given by using the Finite Element Method (FEM) with Comsol Multiphysics 3.5 for the meshes described in Fig. 6. A temperature compatible with a constant intensity field is applied at the boundary of the model and Fig. 7 provides the temperature distribution along the axis parallel to the given intensity field obtained from the FEM results provided by using the mesh described in Fig. 6a. Owing to the fact that the phase configurations are dilute, it is meaningful and efficient to compare the analytical estimates obtained for homogeneous flux boundary conditions with the numerical results delivered by the FEM for homogeneous intensity boundary. It can be seen from Fig. 7 that any modelling provides a resulting heat flux field which is uniform within the inclusion. Due to the higher value of conductivity taken within the inclusion, it can be observed that the heat flux inside the inclusion can be higher (for small interface effect) or smaller (for high interface effect) than the heat flux at infinity. However, we observe also from Fig. 7 that there is a difference for the heat flux field inside the inclusion and matrix phases between both solutions. This difference can be explained by the fact that the closed-form solution obtained in Section 3 is related to a matrix extending to infinity while the numerical solution obtained by FEM is calculated for a finite square domain. A second computation was performed for the mesh defined in Fig. 6b where the size of the square domain is increased from 10µm × 10µm to 20µm × 20µm. As for the first example, we observe from Fig. 8 that the heat flux field is uniform in the inclusion. It can be seen now from Fig. 8 that the proposed analytical solution presents a complete agreement with the numerical solution by FEM both inside and outside the inclusion, meaning that the size of the square domain is now such that the solution is no more affected by the limited extension of the meshed domain.

A second numerical example is provided in order to study the size-dependence of the effective thermal conductivity of the composite. In this example, the matrix and inclusion phases are assumed to be both isotropic and the inhomogeneity radius R takes different values between 1 and 50µm while the inhomogeneity volume fraction is kept constant. Additionally, the inclusion phase is supposed to be less conducting than the matrix one. The thermal conductivities of the matrix and of the inhomogeneities are :

k 1 = 0.1W m -1 K -1 , k 2 = 1W m -1 K -1 . (99)
The ratio k ef f /k ef f c , where k ef f c is computed without taking into account the interface and k ef f is obtained from various models of section 4, is plotted in Fig. 9 in terms of the inhomogeneity radius R. k ef f is the effective thermal conductivity The ratio k ef f /k ef f c in terms of the inhomogeneity volume fraction c 1 is presented in Fig. 10 for an inhomogeneity radius R = 1µm and a Kapitza thermal resistance α = 10 -5 m 2 K/W . It can be seen from Fig. 10 that when the inhomogeneity volume fraction increases, the values of the ratio k ef f /k ef f c obtained from closed-form solutions decrease with the the concentration of inclusions. This result depicts a size effect which increases with the concentration. In addition, all models give the same results for low concentrations (below, say, 0.1), but differ for higher concentrations.

Concluding remarks

Being complementary with the recent work of Le Quang et al. [29] devoted to the highly conducting imperfect interface model, the present work is based on Kapitza's concept of thermal resistance and contributes to the development of the micromechanics of inhomogeneous materials with imperfect interfaces. However, both problems have quite different solutions because, in the present work, it is noteworthy that the solution has been obtained by using an extension to conduction problems of the classical Eshelby's inclusion problem in the case of a polynomial heat free intensity field.

In the context of elasticity, Eshelby's problem for a spherical inclusion embedded in an infinite matrix via an interface characterized by the linear spring-layer model was solved, in the case where the constituent material of the matrix is isotropic, by Hashin [19] and Zhong and Meguid [46] who showed that Eshelby's strain tensor field inside the spherical inclusion is generally no longer uniform. In the case with anisotropic elastic material constituting the matrix, Qu [36] obtained a semi-analytical solution by using an iterative computation, leading again to a non-uniform strain field inside the inclusion. Thus, Eshelby's formalism allowing the transformation of an inhomogeneity problem into an inclusion problem is no more valid. This is because Eshelby's formalism is based on the uniformity of the strain field inside an ellipsoidal inclusion embedded in an infinite matrix and subjected to a prescribed uniform eigenstrain. Unlike the results of Hashin [19], Zhong and Meguid [46] and Qu [36] for elasticity, the results given in the present paper for thermal conduction show that the temperature gradient field inside a circular or spherical inclusion remains uniform even in the presence of a Kapitza interface and holds for any thermal anisotropy of materials constituting the matrix and the inclusions. This fact allows to preserve the Eshelby's formalism and to apply some well-known micromechanical schemes to obtain closed-form expressions for the effective conductivity of composites with Kapitza thermal interface resistance.

Finally, it is well known that due to the physical analogy existing between the different transport phenomena, the approach elaborated and results obtained by the present work for the thermal conduction phenomenon are straightforwardly transposable to other classes of problems such as electric conduction, dielectrics, magnetism, diffusion and flow in porous media. Moreover, owing to the close relation existing between the depolarization tensor in electrostatics and Eshelby's conduction tensor under consideration (see e.g. [45] and [START_REF] Giordano | [END_REF]), all results previously obtained for the thermal conduction phenomenon can be directly generalized to the electromagnetic phenomenon. Additionally, since a mathematical correspondence between anti-plane elasticity and 2D thermal conduction exists, new results for anti-plane anisotropic elasticity can be deduced directly from the 2D anisotropic results derived above. with k * (z) = (K 0 ij z i z j ) -1 . The Fourier inversion of Eq. (A-4) , using Eq. (A-5) is given by

e p (x) = -(2π) -3 Θ z p z j k * (z)dz ∞ 0 |ξ| 2 exp(-i |ξ| z • x)d |ξ| × ω q * j (y) exp(i |ξ| z • y)dy (A-6)
where Θ is the surface of the unit sphere.

Next, owing to the fact that

∇ 2 x cos {|ξ| z • (y -x)} = -|ξ| 2 cos {|ξ| z • (y -x)}, (A-7) and ∞ 0 cos {|ξ| z • (y -x)}d |ξ| = πδ[z • (y -x)], (A-8)
where

∇ 2 x • = ∂ 2 • ∂x s ∂x s , (A-9)
Eq. (A-6) can be rewritten in the following equivalent form

e p (x) = 1 8π 2 Θ z p z j k * (z)dz∇ 2 x ω q * j (y)δ(z • y -z • x)dy. (A-10)
We remark that this equation holds for an arbitrary thermal anisotropy of material forming Ω and for an arbitrary closed inclusion ω. However, we are interested now in the important case where ω is an ellipsoidal inclusion and can be expressed by the following equation of its surrounding surface:

3 γ=1 x 2 γ /a 2 γ = 1, (A-11)
where a γ are semi principal axes of the ellipsoidal inclusion.

Owing to the fact that a polynomial of degree M in {x i } remains a polynomial of same degree M under any change of a system of Cartesian coordinates {x i } by rotation or translation, without loss of generality, we can choose a system of Cartesian coordinates {x i } such that the {x i } axes coincide with the principal axes of the ellipsoidal inclusion.

Next, in order to calculate the intensity field e p (x) given by Eq. (A-10), the polynomial q * j (y) of degree M can be expressed in the following general form where [ 1 2 r] is the integer part of 1 2 r. First, we remark from Eq. (A-28) that I (r) j (s • t) with x ∈ ω is a polynomial of degree r + 2 in x. Thus, the resulting intensity field e(x) given by Eq. (A-18) and the corresponding heat flux field q(x) = K 0 • [e(x)e * (x)] inside the inclusion ω are polynomials of degree M in x.

q * j (y) = M r=0 A ji1i2
Next, in order to calculate integral (A-18), it can be seen from Eq. (A-14) that g pj (z) = g pj (z 1 , z 2 , z 3 ) with

z 1 = µs 1 a 1 , z 2 = µs 2 a 2 , z 3 = µs 3 a 3 (A-29)
However, since g pj (z 1 , z 2 , z 3 ) are rational fractions with the numerator and denominators being both homogeneous polynomials of degree 2 of z 1 , z 2 and z 3 then g pj (z) can be computed alternatively from:

z 1 = s 1 a 1 , z 2 = s 2 a 2 , z 3 = s 3 a 3 , (A-30)
Moreover, the components of the unit vector s can be expressed in the following form

s 1 = (1 -s 2 3 ) 1/2 cos ζ, s 2 = (1 -s 2 3 ) 1/2 sin ζ, s 3 = s 3 . (A-31)
Thus, by substituting ds = ds 3 dζ into Eq. (A-18), we obtain for x / ∈ ω. Next, for M = 1, in addition to the polynomial of degree 0, the polynomial q * j (y) contains a polynomial which is homogeneous and linear in y, i.e. q * jα (y) =

e p (x) = 1 8π 2 1 -1 ds 3 2π 0 µ 3 g pj (z) M r=0 ∇ 2 x I (r) j (s • t) dζ, (

Application to the case with elliptic or circular inclusion

The 2D case where ω is an elliptic or circular inclusion is now considered. This case can be considered as a special case of the one described above related to an ellipsoidal inclusion which is now degraded to a cylindrical one, which corresponds to a 3 → ∞. Next, for a 3 → ∞ Eq. (A-30) implies that z 3 = 0. Therefore, the homogeneous functions g pj (s 1 /a 1 , s 2 /a 2 , 0) depend no more on s 3 and the values of s 1 and s 2 can be computed using s 3 = 0. As a consequence s 1 and s 2 have now the simple form

s 1 = cos ζ, s 2 = sin ζ, (A-52)
The intensity fields e(x) for interior points given by Eq. (A-41) with x = (x 1 , x

Figure 1 .

 1 Figure 1. Two-and three-phase configurations: (a) matrix/interphase/inclusion composite; (b) matrix/inclusion composite with imperfect interface exhibiting the Kapitza thermal resistance.

  Quang, Q.-C. He and G. Bonnet

Figure 2 .

 2 Figure2. The Eshelby's conduction problem for the thermal conduction phenomenon consisting of a subdomain ω embedded in an infinite medium Ω , subjected to a uniform heat-free intensity e 0 within ω and vanishing outside ω.

11 Figure 3 .

 113 Figure3. The conjugate Eshelby's conduction problem concerned with a subdomain ω embedded in an infinite medium Ω and subjected to a uniform heat-free flux equal to q 0 within ω and vanishing outside ω.

  MS-HLQ-QCH-GB-RV 12 H. Le Quang, Q.-C. He and G. Bonnet

Figure 4 .

 4 Figure 4. Decomposition of the initial heterogeneity problem into three Sub-problems (A), (B) and (C) with q ω being the uniform heat flux solution field of the initial heterogeneity problem within ω.

  a) Subproblem AThis Sub-problem is concerned with the determination of the heat flux solution field within the infinite homogeneous matrix Ω subjected to the homogeneous flux boundary condition (49) on its external surface ∂Ω. Clearly, the corresponding MS-HLQ-QCH-GB-RVPhilosophical Magazine Letters 13

  Quang, Q.-C. He and G. Bonnet

  MS-HLQ-QCH-GB-RVPhilosophical Magazine Letters localization tensor fields in Section 3, it can be shown that the heat flux and intensity fields inside the inhomogeneity ω are uniform and given by q

  Quang, Q.-C. He and G. Bonnet

4 . 3 .

 43 Self-consistent model

  Quang, Q.-C. He and G. Bonnet

Figure 5 .

 5 Figure 5. Generalized self-consistent model applied to thermal conduction problem.

  Quang, Q.-C. He and G. Bonnet

Figure 6 .

 6 Figure 6. Meshing and temperature solution field distributed in the matrix and inclusion phases for the heterogeneity problem undergoing a uniform heat flux vector q 0 = (1, 0) T [W/m 2 ] at the boundary with α = 10 -5 m 2 K/W and solved by FEM: (a) Case 1 with inclusion radius R = 1µm embedded in a square matrix domain 10µm × 10µm; (b) Case 2 with inclusion radius R = 1µm embedded in a square matrix domain 20µm × 20µm.

2 ]Figure 7 . 2 ]

 272 Figure 7. The heat flux fields in the matrix and inclusion obtained by the present analytical method and by FEM for Case 1 with inclusion radius R = 1µm embedded in a square matrix domain 10µm × 10µm.

Figure 8 .Figure 9 .

 89 Figure 8. The heat flux fields in the matrix and inclusion obtained by the present analytical method and by FEM for Case 2 with inclusion radius R = 1µm embedded in a square matrix domain 20µm × 20µm.calculated by Eqs. (66), (73), (79) and (96) with α = 10 -5 m 2 K/W according to the dilute distribution, Mori-Tanaka, self-consistent or generalized self-consistent models, respectively. It can be seen from Fig.9that: (i) k ef f depends on the inhomogeneity radius R while k ef f c (corresponding to the "surface C" case) is obviously independent of R; (ii) the difference between k ef f and k ef f c decreases when R in-

Figure 10 .

 10 Figure 10. The ratio k ef f /k ef f c versus the inclusion volume fraction c 1 with inclusion radius R = 1µm and Kapitza thermal resistance α = 10 -5 m 2 K/W .

j

  (s • t) given by Eq. (A-28). (ii) For exterior points of the inclusion ω When x is outside the inclusion ω, or equivalently |t| > 1, the delta functionδ[µ(s • τs • t)] = 1 µ δ(s • τs • t) is non-zero only if |s • t| ≤ 1 and s • τs • t = 0.These conditions imply that τ must take the same form as in Eq. (A-19) under the condition |s • t| ≤ 1.

Figure 11 .|t| - 1 g|t| - 1 g

 1111 Figure 11. The subdomains Θ * and L of Θ satisfying the condition |s • t| ≤ 1 and s • t = 1 and the chosen new coordinate system {s 1 , s 2 , s 3 } for a given vector t

3 α=1A 3 α=1A 3 α=1A jα Ξ αγ a α 2π 0 |t| - 1 s|t| - 1 s

 333011 jα y α . The values of f (r) (ψ) and I (r) j (x, z) defined by Eqs. (A-23) and (A-28) become:f (0) (ψ) = 1, f (0) α (ψ) = s α , f (1) jα s α a α (s • t) 1 -(s • t) 2 , (A-46)Owing to the fact that∇ 2 x (s • t) 3 = 6 µ 2 (s • t) , (A-47) γ g pj (s ) s 3 =1/|t| dζ (A-49)outside the inclusion. The comparison of Eqs. (A-48) and (A-49) with Eq. (A-3) implies thatN pji (x) = ∂ 2 ∂x p ∂x j ω G(yx)y i dy = s • t) g pj (z)dζ (A-50)when x ∈ ω andN pji (x) = ∂ 2 ∂x p ∂x j ωG(yx)y i dy γ g pj (s ) s 3 =1/|t| dζ (A-51) when x / ∈ ω.

  ...ir y i1 y i2 ...y ir , (A-12) where A ji1i2...ir are constant.

	April 3, 2011 7:17 Philosophical Magazine Letters MS-HLQ-QCH-GB-RV April 3, 2011 7:17 Philosophical Magazine Letters MS-HLQ-QCH-GB-RV
	F o r F o r Philosophical Magazine Letters F o r
	P P P
	e e e
	e r e r e r
	R R R
	e e e
	v i e v i e v i e
	w w w
	O n l O n l O n l
	y y y

  Le Quang, Q.-C. He and G. BonnetFor exterior points, Eq. (A-42) leads to: When a given heat-free intensity within an elliptic or cylindrical inclusion is linear, the intensity fields e(x) inside and outside the inclusion are determined by (s 1 t 1 + s 2 t 2 )g pj (s 1 /a 1 , s 2 /a 2 , 0)dζ (A-57) Finally, we obtain from Eqs. (A-57) and (A-58) the expressions of N pji (x) with p, j and i = 1 or 2 for the case where the inclusion is elliptic or cylindrical (s 1 t 1 + s 2 t 2 )g pj (s 1 /a 1 , s 2 /a 2 , 0)dζ (A-59)

	April 3, 2011 7:17 Philosophical Magazine Letters MS-HLQ-QCH-GB-RV April 3, 2011 7:17 Philosophical Magazine Letters MS-HLQ-QCH-GB-RV
	34 when x ∈ ω and H. e p (x) = -1 4π A j N pji (x) = 3 4π Ξ iγ a i + 1 2π A j 2π Philosophical Magazine Letters 1/|t| -1/|t| ds 3 2π 0 g pj ( Ξ 1m s m a 1 , 1/|t| -1/|t| ds 3 2π 0 s γ s 3 |t|g pj ( Ξ 1m s m Ξ 2m s m a 1 , Ξ 2m s m , 0)dζ a 2 a 2 0 |t| -1 g pj ( Ξ 1m s m a 1 , Ξ 2m s m a 2 , 0) s 3 =1/|t| -1 2π Ξ iγ a i 2π 0 |t| -1 s γ g pj ( Ξ 1m s m a 1 , Ξ 2m s m s 3 =1/|t| , 0)dζ dζ . a 2 , 0) dζ , (A-60) 35 (A-54)
	Correspondingly, D pj (x) with p and j = 1 or 2 takes the following form when x / ∈ ω.
	2 , 0) g (A-53) ds 3 2π 0 g pj (s 1 /a 1 , s 2 /a 2 , 0)dζ 2π 1 -1 0 e e p (x) = -1 4π A j = -1 2π A j F reduce now to o r P e r R e v i e w O n l y D pj (x) = ∂ 2 ∂x p ∂x j ω G(y -x)dy = 2π 1 g pj (s 1 /a 1 , s 2 /a 2 , 0)dζ (A-55) 2π 0 when x ∈ ω and D pj (x) = ∂ 2 ∂x p ∂x j ω G(y -x)dy = 1 4π 1/|t| -1/|t| ds 3 × 2π 0 g pj ( Ξ 1m s m a 1 , Ξ 2m s m a 2 , 0)dζ -1 2π 2π 0 |t| -1 g pj ( Ξ 1m s m a 1 , Ξ 2m s m a 2 , 0) s 3 =1/|t| dζ (A-56) when x / ∈ ω. e p (x) = -1 π 2 α=1 A jα a α 2π 0 s α when x ∈ ω and e p (x) = -3 4π 2 α=1 A jα Ξ αγ a α 1/|t| -1/|t| ds 3 2π 0 s γ s 3 |t|g pj ( Ξ 1m s m a 1 , Ξ 2m s m a 2 , 0)dζ + 1 2π 2 α=1 A jα Ξ αγ a α 2π 0 |t| -1 s γ g pj ( Ξ 1m s m a 1 , Ξ 2m s m a 2 , 0) s 3 =1/|t| dζ , (A-58) when x / ∈ ω. N pji (x) = a i 0 π 2π o r P e e r R e v i e w O n l s i F y

pj (s 1 /a 1 , s 2 /a 2 , 0)dζ.

http://mc.manuscriptcentral.com/pm-pml Philosophical Magazine & Philosophical Magazine Letters

Appendix : temperature field induced by a polynomial heat-free intensity inside an ellipsoidal inclusion

The intensity field induced by a constant heat-free intensity inside an ellipsoidal inclusion is obtained from tensor N mij defined in Eq. (29) which is related to a linearly varying heat-free intensity inside an ellipsoidal inclusion. In this appendix, similarly to the generalized Eshelby's results in the case of elasticity (see e.g. Asaro and Barnett [1] and Mura [START_REF] Mura | Mircomechanics of Defects in Solids[END_REF]) we show in the thermal conduction Eshelby's problem with perfect interfaces that: when an ellipsoidal region in an infinite anisotropic medium undergoes a heat-free intensity e * (x) which is a polynomial of degree M in the position coordinates x, then the resulting intensity and heat flux fields in this ellipsoidal region are also polynomials of degree M in x.

The results described in [START_REF] Mura | Mircomechanics of Defects in Solids[END_REF] and [1] are completely developed in the case of elasticity, more specifically in the description of the strain field outside inclusions produced by polynomial eigenstrains of order 0 and 1. These results cannot however be easily extended to the case of conduction problems. However, the process leading to the solution described in the previously cited references can be used again. The main difference is now that, instead of using harmonic and biharmonic potentials as in [START_REF] Mura | Mircomechanics of Defects in Solids[END_REF] and [1], only harmonic potentials are now necessary.

As described above in Section 3, we consider now a subdomain ω embedded in a 3D infinitely extended domain Ω constituting of a homogeneous material whose thermal conductivity and resistivity tensors are K 0 and H 0 , respectively. The interface Γ between inclusion ω and external medium is assumed to be perfect.

By using the Green function, G(yx) and by applying the divergence theorem, the temperature field T (x) produced by e * (x) can be expressed in terms of G(y-x) and e * (x) by (see e.g. Le Quang et al. [28])

where the heat flux field q * (x) relates to the heat-free intensity e * (x) by

It is clear from the last relation that if e * (x) is a polynomial of degree M in the position coordinates x then the heat flux field q * (x) is also a polynomial of degree M in x. By using the property (26) of G(yx), the resulting intensity field components are given by

Now, this equation can be written in Fourier space as follows

where i = √ -1, ξ is the discrete wave vector and

. The Fourier transform Ĝ(ξ) of G(yx) can be derived from Eq. ( 17) and is given by Ĝ

Substituting Eq. (A-12) into Eq. (A-10), we have

with

It can be shown, by the change of variables

where the summation convention rule does not apply for index γ, that the integral (A-15) reduces to an integral over the unit sphere, |τ| ≤ 1, i.e.,

, s γ = z γ a γ µ (no sum on γ).

By transformation of variable z by s in Eq. (A-13), we obtain the following equivalent form

with I (r) j (x, z) given by Eq. (A-17).

Explicit computation of the intensity field e(x)

In the following, the previous results, and more specifically (A-17) and (A-18) are provided in a more explicit form. They are studied in a first step for interior points and next for exterior points.

(i) For interior points of the inclusion ω We note that the delta function δ

When x is inside the inclusion ω, or equivalently when |t| ≤ 1, then |s • t| ≤ 1 and s • τs • t = 0 if and only if τ takes the following form:

where 0 ≤ η ≤ 1 -(s • t) 2 and u and v are two orthonormal basis vectors of the plane normal to s, while ψ is a polar angle in this plane. Thus, the volume integral (A-15) reduces to an integral over a plane circular region of radius 1 -(s • t) 2 as follows

where

with h = cos ψu + sin ψv.

Applying the binomial theorem for J (r) (s, t, ψ) yields

i1,i2,...,ir (ψ) = h i1 s i2 ...s ir + s i1 h i2 ...s ir + ... + s i1 s i2 ...h ir , . . .

Moreover, it is easy to show that

Finally, by combining Eq. (A-20) with Eqs. (A-22)-(A-27), we obtain the following expression of

i1,i2,...,ir (ψ)dψ

By using a process similar to the one described above for interior points of inclusion ω, we obtain the following expression of the intensity field e(x) for exterior points of the inclusion ω (see, e.g. [START_REF] Mura | Mircomechanics of Defects in Solids[END_REF])

where Θ * and L are subdomains of Θ satisfying the condition |s • t| ≤ 1 and s•t = 1, respectively; I (r) j (s • t) is given by Eq. (A-28). More precisely, the integral on Θ * and L can be calculated by introducing a new coordinate system {s 1 , s 2 , s 3 } in which s 3 -axis is taken in the direction of t, and s 1 -and s 2 -axis are given along two perpendicular directions in the plane normal to t (see Fig. 11). We can then express the relation between {s 1 , s 2 , s 3 } and {s 1 , s 2 , s 3 } by

where Ξ ij are the direction cosines between the s i -and s j -axis. Since

and by setting (A-37)

Application to the cases of constant and linear polarization tensors

The following particular cases with M = 0 and 1 are studied in detail, because they are explicitly used within the paper.

First, for M = 0, the polynomial q * j (y) takes the simple form q * j (y) = A j . By using Eqs. (A-23) and (A-28), we obtain the expressions of f (0) (ψ) and I (0) j (x, z) as follows f (0) (ψ) = 1, (A-38)

(A-39)