D.-H Nguyen

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

Q.-C He

T.-T Nguyen

A.-T Tran

The effective thermal conductivity of composites with interfaces oscillating in two directions around a curved surface

Keywords: micromechanics, composites, thermal conductivity, asymptotic analysis, rough interfaces

In many situations of practical or/and theoretical interest, the assumption that the interfaces between constituent phases of a composite are smooth is no longer appropriate, and the consideration of rough interfaces at microscopic scale is necessary. However, in micromechanics, when the interfaces between the constituent phases of composites become rough, all classical well-known micromechanical schemes resorting to Eshelby's formalism cannot be applicable and the problem of determining the effective properties of composites become largely open. The present work aims to determine the effective thermal conductivity of a composite in which the interfaces between its constituent phases are perfectly bonded but oscillate quickly around a curved surface and along two directions. To achieve this objective, a two-scale homogenization method is proposed. In the first-scale homogenization, or microscopic-to-mesoscopic upscaling, the interfacial zone in which the interface oscillates is homogenized as an equivalent interphase by applying an asymptotic analysis. The thermal properties of the equivalent interphase can generally be determined by using a numerical approach based on the fast Fourier transform (FFT) method. In particular case where the equivalent interphase is very thin, this interphase is then replaced with a general imperfect interface situated at its middle surface. By applying the equivalent inclusion method, every inclusion with imperfect interface is further substituted by an equivalent inclusion with perfect interface. In the second scale homogenization, or mesoscopic-to-macroscopic upscaling, due to the fact that the interfaces are perfect, the effective thermal conductivity can be analytically obtained by using some well-known classical micromechanical schemes. To illustrate the two-scale homogenization method proposed in this work, the case of a layered

Introduction

In physics and mechanics of composite materials, most of investigations dedicated to determining the effective properties of composite materials in terms of local phase properties and their microstructure often adopt the hypothesis that the interfaces between their constituent phases are smooth. When this assumption about the interfaces of composite materials holds, classical approximation schemes, such as diluted, self-consistent, Mori-Tanaka, differential approximation models, resort to Eshelby's tensor for estimating the effective physical and mechanical properties. However, in many situations of practice, the assumption of smooth interfaces is too idealistic and the consideration of rough interfaces is unavoidable. Another necessity of accounting for rough interfaces in studying composite materials resides in the fact that an interface considered as nominally smooth at a given scale may be rough at a finer scale. Consequently, the problem of determining the effective properties of composites while taking account into the roughness of their interfaces between constituent phases is an important issue, which has captured attention of some scientists. In particular, a number of studies in the fields of physics and mechanics of solids have been dedicated to homogenizing the interfacial zone in which a rough interface oscillates. We can cite the studies [START_REF] Kohler | Boundary and interface problems in regions with very rough boundaries[END_REF][START_REF] Nevard | Homogenization of rough boundaries and interfaces[END_REF][START_REF] Kristensson | Homogenization of corrugated interfaces in electromagnetics[END_REF][START_REF] Vinh | Homogenized equations of the linear elasticity in two-dimensional domains with very rough interfaces[END_REF][START_REF] Vinh | Homogenization of rough two-dimensional interfaces separating two anisotropic solids[END_REF][START_REF] Vinh | Homogenized equations of the linear elasticity theory in two-dimensional domains with interfaces highly oscillating between two circles[END_REF] in which the interfacial zone is homogenized by applying asymptotic analysis. In the work of Le-Quang et al. [START_REF] Le-Quang | Multiscale homogenization of elastic layered composites with unidirectionally periodic rough interfaces[END_REF], it was shown that, when the rough interface oscillates quickly and periodically along only one direction, the effective properties of the equivalent interphase obtained by homogenizing an interface zone correspond exactly to the ones of a two-phase layered composite material. In addition, the determination of effective properties of composite materials while accounting for interfacial roughness has been carried out in the contexts of elasticity and thermal conduction [START_REF] Le-Quang | Multiscale homogenization of elastic layered composites with unidirectionally periodic rough interfaces[END_REF][START_REF] Le | The effective elastic moduli of columnar composites made of cylindrically anisotropic phases with rough interfaces[END_REF][START_REF] Nguyen | Determination of the effective conductive properties of composites with curved oscillating interfaces by a two-scale homogenization procedure[END_REF][START_REF] Le-Quang | Two-scale homogenization of elastic layered composites with interfaces oscillating in two directions[END_REF]. the results obtained in these works are limited to the cases where rough interfaces oscillate either along one direction around a curved surface [START_REF] Le-Quang | Multiscale homogenization of elastic layered composites with unidirectionally periodic rough interfaces[END_REF][START_REF] Le | The effective elastic moduli of columnar composites made of cylindrically anisotropic phases with rough interfaces[END_REF][START_REF] Nguyen | Determination of the effective conductive properties of composites with curved oscillating interfaces by a two-scale homogenization procedure[END_REF] or along two directions but around a plane surface [START_REF] Le-Quang | Two-scale homogenization of elastic layered composites with interfaces oscillating in two directions[END_REF]. The present work aims to solve, in the situation of thermal conduction, the problem of determining the effective conductivity of composite materials in which the interface between any constituent phases oscillates quickly and periodically about an arbitrarily curved surface and along two directions. Thus, the present work can be considered as a continuation and an extension of our previous ones [START_REF] Le-Quang | Multiscale homogenization of elastic layered composites with unidirectionally periodic rough interfaces[END_REF][START_REF] Le | The effective elastic moduli of columnar composites made of cylindrically anisotropic phases with rough interfaces[END_REF][START_REF] Nguyen | Determination of the effective conductive properties of composites with curved oscillating interfaces by a two-scale homogenization procedure[END_REF][START_REF] Le-Quang | Two-scale homogenization of elastic layered composites with interfaces oscillating in two directions[END_REF].

To achieve the objective mentioned above, a two-scale homogenization method is proposed in the present work. The first-scale homogenization, or microscopic-to-mesoscopic upscaling, is dedicated to homogenizing, via an asymptotic analysis, the interfacial zone in which the interface undulates as an equivalent interphase whose thermal conductivity matrix components are expressed and computed with respect to a system of orthogonal curvilinear coordinates. Unlike the case of unidirectionally periodic rough interfaces, the thermal properties of the homogenized interphase cannot be analytically and explicitly determined. To overcome this difficulty, a transformation from an orthogonal curvilinear coordinates system to a Cartesian coordinates system is first realized and a numerical approach based on the fast Fourier transform (FFT) method is then applied to numerically compute all thermal conductivity matrix components of the homogenized interphase. Particularly, for the case where the equivalent interphase obtained by the first-scale homogenization is very thin, this equivalent interphase can be then replaced by a general imperfect interface situated at its middle surface. By applying the equivalent inclusion method proposed recently by Nguyen et al. [START_REF] Nguyen | Generalized Hill-Mendel's lemma and equivalent inclusion method for determining the effective thermal conductivity of composites with imperfect interfaces[END_REF], an inclusion with imperfect interface is further substituted by an equivalent inclusion with perfect interface. The thermal conductivity of the equivalent inclusion is determined in such a way that the thermal energy does not change before and after this substitution. The second scale homogenization, or mesoscopic-to-macroscopic upscaling, concerns the determination of the effective (or macroscopic) properties of composites by employing some well-known micromechanical schemes like diluted, self-consistent, Mori-Tanaka, differential approximation models. It is important to notice that all classical micromechanical schemes are actually applicable because the interfaces between the matrix and equivalent inclusion phases of composites are now not only perfect but also smooth.

Layered composite and cylindrical composites are two important classes of composites for the practical and theoretical points of view. Indeed, many natural and man-made composites, such as geological stratified strata and artificial laminates fall within the class of layered composites while fibrous composites consisting of a homogeneous matrix reinforced by aligned parallel continuous homogeneous fibers, porous media containing parallel cylindrical pore or polycrystalline aggregates formed of columnar monocrystals are considered as some examples of cylindrical composites. From the theoretical standpoint, layered composites are homogeneous in the plane of layers but heterogeneous along the layering direction. The simple microstructure of layered composites leads to the fact that their effective properties can be analytically and exactly determined (see e.g. [START_REF] Milton | The Theory of Composites[END_REF][START_REF] He | Homogenization of layered elastoplastic composites: Theoretical results[END_REF]). Concerning cylindrical composites, their microstructures are such that their properties are homogeneous along the axial direction but heterogeneous in its transverse plane. In view of the importance of layered and cylindrical composites aforementioned and in order to illustrate the efficiency and accuracy of the two-scale homogenization method proposed in this work, two important examples of application in which a layered composite with rough interfaces oscillating around a plane surface and a composite cylinder with rough interface undulating around a circumferential surface are studied in details. In the two cases of layered and cylindrical composites, the analytical or semi-analytical results obtained by the two-scale homogenization method for the effective thermal conductivities are compared with the numerical ones provided by the finite element method (FEM) as well as with the Reuss, Voigt and Hashin-Shtrikman bounds.

The paper is organized as follows. Section 2 is dedicated to specifying the setting of the problem. In section 3, the general two-scale homogenization method is presented. In the first-scale homogenization process, an asymptotic analysis is carried out at the mesoscopic scale to homogenize a rough interfacial zone as an equivalent interphase. Then, this equivalent interphase is replaced by a general imperfect interface for the case where the equivalent interphase is very thin. Next, each inclusion with its general imperfect interface is be substituted by an equivalent inclusion with perfect interface. In the second-scale homogenization process, by applying some classical micromechanical schemes, the effective thermal conductivities of composites are obtained. Section 4 is dedicated to studying in detail a layered composite with rough interfaces oscillating around a plane surface and in two directions. Section 5 is related to the second example of application in which a composite cylinder possesses rough interfaces oscillating around the circumferential surface and in two directions. Finally, a few conclusions and remarks are given in section 6.

Problem setting

In a three-dimensional (3D) euclidian space R 3 , we consider a composite Ω consisting of two subdomains, Ω (1) and Ω (2) , formed of an inclusion phase 1 and a matrix phase 2, respectively. In addition, these two phases, Ω (1) and Ω (2) , are assumed to be perfectly bounded together at their interface Γ which oscillates around a curved surface. In order to obtain a mathematical characterization of Γ , we denote by (x 1 , x 2 , x 3 ) a system of Cartesian coordinates associated to an orthogonal basis {j 1 , j 2 , j 3 }, and by (y 1 , y 2 , y 3 ) a system of orthogonal curvilinear coordinates. The position x of any point of R 3 can be defined by x = x(y 1 , y 2 , y 3 ) = [x 1 (y 1 , y 2 , y 3 ), x 2 (y 1 , y 2 , y 3 ), x 3 (y 1 , y 2 , y 3 )].

(1)

The tangent vector to the y i -curved coordinate is defined by

t i = ∂x ∂y i = h i f i with h i = ∂x ∂y i , (2)
where the summation convention does not apply, h i is a metric coefficient and f i is the unit tangent vector to the y i -curved coordinate. Since (y 1 , y 2 , y 3 ) is an orthogonal curvilinear coordinate system, then f 1 , f 2 and f 3 satisfy f i • f j = δ ij with δ ij being the Kronecker symbol defined as δ ij = 1 for i = j and δ ij = 0 for i = j. The transformation from Cartesian coordinates to curvilinear coordinates is illustrated in Fig. 1. In this work, we are interested in the case where the interface Γ is periodically oscillated along the y 1 -and y 2 -directions with respective periods 1 and 2 . The values of 1 and 2 are j ( )

x 1 1 j ( ) x 2 2 j ( ) x 3 3 O f ( ) y 1 1 f ( ) y 2 2 f ( ) y 3 3 y =c 3 3 y =c 2 2
y =c

1 1
Fig. 1: Transformation from the Cartesian to curvilinear coordinates.

assumed to be of the same order of magnitude. For simplicity, the periods 1 and 2 can be expressed as 1 = η 1 and 2 = η 2 . In the orthogonal curvilinear coordinate system (y 1 , y 2 , y 3 ), the interface Γ between Ω (1) and Ω (2) is defined by

Γ = y = y(y 1 , y 2 , y 3 ) ∈ R 3 | y 3 = γ(ξ 1 , ξ 2 ), ξ 1 = y 1 1 , ξ 2 = y 2 2 (3)
where γ(ξ 1 , ξ 2 ) is a periodic function of period 1 in both ξ 1 and ξ 2 . Without loss of generality, this function γ(ξ 1 , ξ 2 ) is supposed to have the minimal value γ min and maximal value γ max satisfying γ max = -γ min = δ 2 with δ standing for the thickness of the rough interface zone where the interface Γ oscillates. In addition, we assume also that 0 < max{ 1 , 2 } δ. This means that Γ is a very rough interface in both y 1 -and y 2 -directions. The rough interface zone denoted by ω (c) and characterized by

O f ( ) y 1 1 δ/2 -/2 δ δ f ( ) y 2 2 f ( ) y 3 3
Fig. 2: Curved surface oscillating along two curvilinear coordinate directions

ω (c) = y = y(y 1 , y 2 , y 3 ) ∈ Ω | - δ 2 ≤ y 3 ≤ δ 2 , (4)
is situated between two smooth planes Γ 1 and Γ 2 specified by

Γ 1 = y = y(y 1 , y 2 , y 3 ) ∈ ω (c) | y 3 = - δ 2 , Γ 2 = y = y(y 1 , y 2 , y 3 ) ∈ ω (c) | y 3 = δ 2 . ( 5
)
It is convenient for later use to introduce a generic smooth surface Π(y 3 ) defined by

Π(y 3 ) = y = y(y 1 , y 2 , y 3 ) ∈ ω (c) | y 3 = y 3 , - δ 2 < y 3 < δ 2 . (6)
Clearly, it can be deduced from ( 5) and ( 6) that Γ 1 = Π -δ 2 and Γ 2 = Π δ 2 (see Fig. 2). At the same time, we define as follows two sub-domains, denoted by ω (1) and ω (2) , which are outside ω (c) but belong to Ω (1) and Ω (2) , respectively, ω (1) = Ω (1) \(Ω (1) ∩ ω (c) ), ω (2) = Ω (2) \(Ω (2) ∩ ω (c) ). [START_REF] Le-Quang | Multiscale homogenization of elastic layered composites with unidirectionally periodic rough interfaces[END_REF] With respect to the curvilinear coordinate system {y 1 , y 2 , y 3 } associated to the orthonormal curvilinear basis {f 1 , f 2 , f 3 }, the local thermal behavior in Ω is described by the following Fourier's law

q(y) = -K(y) • ∇θ(y) (8)
where q(y) and θ(y) stand for the heat flux and temperature fields, respectively, and K(y) denotes the second-order local thermal conductivity tensor that is characterized by

K(y) = Φ (1) (y)K (1) + Φ (2) (y)K (2) . (9)
In this equation, Φ (α) (y) is the characteristic function of the subdomain Ω (α) such that Φ (α) (y) = 1 if y ∈ Ω (α) and Φ (α) (y) = 0 otherwise. Relative to the curvilinear coordinate system {y 1 , y 2 , y 3 }, the matrix and inclusion phases are assumed to be homogeneous in the sense that K (α) are independent of the curvilinear coordinates {y 1 , y 2 , y 3 } but they are in fact heterogeneous in terms of the Cartesian coordinates {x 1 , x 2 , x 3 }. Above and hereafter, as a general rule, Latin superscripts and subscripts, such that i, j and k range from 1 to 3 while Greek superscripts and subscripts, like α and β, takes only the value 1 or 2. The temperature gradient ∇θ (α) is calculated in the orthonormal curvilinear basis {f 1 , f 2 , f 3 } by

∇θ (α) = θ (α) ,1 f 1 + θ (α) ,2 f 2 + θ (α) ,3 f 3 ( 10
)
where the derivatives (•) ,1 , (•) ,2 and (•) ,3 are, hereafter, defined as

(•) ,1 = 1 h 1 ∂(•) ∂y 1 , (•) ,2 = 1 h 2 ∂(•) ∂y 2 , (•) ,3 = 1 h 3 ∂(•) ∂y 3 (11)
with h 1 , h 2 and h 3 being the metric coefficients of the curvilinear coordinates y 1 , y 2 and y 3 , respectively. The local heat flux q (α) must satisfy the following energy conservation equation

∇ • q (α) + r (α) = (α) ζ (α) ∂θ (α) ∂t (12)
where r (α) , (α) and ζ (α) are the heat source, mass density and specific heat capacity of phase α, respectively; the divergence of heat flux ∇ • q (α) is expressed in the orthogonal curvilinear coordinates by

∇ • q (α) = 1 h 1 h 2 h 3 ∂(h 2 h 3 q (α) 1 ) ∂y 1 + ∂(h 1 h 3 q (α) 2 ) ∂y 2 + ∂(h 1 h 2 q (α) 3 ) ∂y 3 . ( 13
)
Owing to the fact that the interface Γ between Ω (1) and Ω (2) is assumed to be perfect, so that both the temperature and the normal heat flux component are continuous across Γ :

θ = 0 and K ij θ ,j n i = 0 ( 14
)
where 1) denotes the jump of • and n i stands for the components of the outward normal vector n of Γ directed from Ω (1) into Ω (2) and given by n = -1

• = • (2) -• (
( -1 γ ,ξ1 ) 2 + ( -1 γ ,ξ2 ) 2 + h -2 3 γ ,ξ1 , γ ,ξ2 , -h -1 3 ( 15
)
where the derivatives (•) ,ξ1 and (•) ,ξ2 are, hereafter, defined as

(•) ,ξ1 = ∂(•) η 1 h 1 ∂ξ 1 and (•) ,ξ2 = ∂(•) η 2 h 2 ∂ξ 2 . ( 16
)
Finally, at the macroscopic scale, the composite under consideration is assumed to be statistically homogeneous and Ω is considered as a representative volume element (RVE) of this composite. Relative to the Cartesian coordinate system {x 1 , x 2 , x 3 }, the corresponding effective thermal behavior is characterized by

Q = K eff • E ( 17
)
where K eff is the effective thermal conductivity second-order tensor, Q and E denote the macroscopic heat flux and intensity vectors, respectively. More precisely, these macroscopic intensity vector E and heat flux vector Q are defined as

E = - 1 |Ω| ∂Ω θ(x)ν(x)dS, Q = 1 |Ω| ∂Ω (q • ν)xdS (18)
where ν(x) is the outward unit normal vector to ∂Ω and |Ω| denotes the volume of Ω.

Let the RVE Ω undergo on its boundary ∂Ω either a uniform intensity boundary condition

θ(x) = -E 0 • x, ∀x ∈ ∂Ω (19)
where

E 0 = (E 0 1 , E 0 2 , E 0
3 ) T is a prescribed constant intensity vector, or a uniform heat flux boundary condition

q(x) • ν(x) = Q 0 • ν(x), ∀x ∈ ∂Ω ( 20
)
where

Q 0 = (Q 0 1 , Q 0 2 , Q 0 3 )
T is a prescribed constant heat flux vector or a mixed boundary condition

θ(x) = -E 0 • x, ∀x ∈ ∂Ω E , q(x) • ν(x) = Q 0 • ν(x), ∀x ∈ ∂Ω Q (21)
with ∂Ω E and ∂Ω Q being the complementary parts of ∂Ω. In particular, when [START_REF] Moulinec | A numerical method for computing the overall response of nonlinear composites with complex microstructure[END_REF] corresponds therefore to the thermal insulated one. Moreover, it can be shown that the macroscopic intensity vector E defined by Eq. ( 18) 1 is directly given by

Q 0 1 = Q 0 2 = Q 0 3 = 0, the boundary condition
E = E 0 (22)
when the boundary condition [START_REF] Moulinec | A fast numerical method for computing the linear and nonlinear mechanical properties of composites[END_REF] is used; the macroscopic heat flux Q specified by Eq. ( 18) 2 is immediately determined by

Q = Q 0 ( 23
)
when the boundary condition ( 20) is considered.

General two-scale homogenization method

The homogenization procedure proposed in this work to determine the effective conductivity of composites with interfaces oscillating in two directions around a curved surface corresponds to a general two-scale homogenization method summarized and illustrated in Fig. 3. The details of each step of this homogenization method are represented in the following subsections.

Ω ( 1)

ω ( ) ω ( ) (1)
(2)

c c Ω (1) Ω (2) Ω (1) Ω (2) ̭ ̭ Γ1 Ω (1) Ω (2) ω ω K (1)

K

(1)

K (2) K ( ) c K (1)
K

K

(2)

K (2) K (1)

K

(1) K

K

(2)

Step 1

Step 2

Step 3 In this section, attention is focused on the rough interface zone ω (c) that is geometrically specified by Eq. ( 4) and can be obtained from the periodic distribution of a unit cell defined by

K (1) K (2)
Y = [0, 1 [×[0, 2 [× -δ 2 , δ 2
along both y 1 -and y 2 -directions. The corresponding effective mesoscopic thermal conductivity tensor K (c) of the rough interface zone ω (c) is determined by applying the asymptotic analysis method presented in [START_REF] Papanicolaou | Asymptotic Analysis for Periodic Structures[END_REF][START_REF] Sanchez-Palencia | Non-Homogeneous Media and Vibration Theory[END_REF][START_REF] Bakhvalov | Averaging of Processes in Periodic Media: Mathematical Problems of the Mechanics of Composite Materials[END_REF][START_REF] Persson | The Homogenization Method: An Introduction[END_REF][START_REF] Jikov | Homogenization of differential operators and integral functionals[END_REF]. According to this method, the temperature field is first considered as a function of the macroscopic (or slow) variables y 1 , y 2 , y 3 and the microscopic (or fast) variables ξ 1 and ξ 2 :

θ(y 1 , y 2 , y 3 , t) = θ (y 1 , y 2 , y 3 , ξ 1 , ξ 2 , t). ( 24
)
The temperature field θ (y 1 , y 2 , y 3 , ξ 1 , ξ 2 , t) can be expressed in terms of an asymptotic expansion as follows:

θ (y 1 , y 2 , y 3 , ξ 1 , ξ 2 , t) = T + (N (1) T + N (1i) T ,i ) + 2 (N (2) T + N (2i) T ,i + N (2ij) T ,ij ) + 0( 3 ). ( 25
)
In this expression, T = T (y 1 , y 2 , y 3 , t) is a scalar function independent of the microscopic variables ξ 1 and ξ 2 and corresponds therefore to the overall temperature field at mesoscopic scale;

N (•) = N (•) (ξ 1 , ξ 2 , y 3
) are localization scalar functions that depend only on ξ 1 , ξ 2 and y 3 , and are determined from the energy conservation equation ( 12) and the continuity conditions [START_REF] Papanicolaou | Asymptotic Analysis for Periodic Structures[END_REF].

Since ξ 1 = y 1 / 1 and ξ 2 = y 2 / 2 , the derivatives of the temperature field with respect to the spatial variables y 1 , y 2 and y 3 are given by

θ ,1 = θ ,1 + -1 θ ,ξ1 , θ ,2 = θ ,2 + -1 θ ,ξ2 , θ ,3 = θ ,3 . (26)
In this equation, the derivatives of the temperature field with respect to y 1 , y 2 , y 3 and ξ 1 , ξ 2 are defined by [START_REF] Nguyen | Generalized Hill-Mendel's lemma and equivalent inclusion method for determining the effective thermal conductivity of composites with imperfect interfaces[END_REF] and [START_REF] Bakhvalov | Averaging of Processes in Periodic Media: Mathematical Problems of the Mechanics of Composite Materials[END_REF] where the metric coefficients h 1 , h 2 and h 3 are functions of y 1 , y 2 and y 3 . By introducing Eq. ( 25) into Eq. ( 26), we obtain

θ ,α = T ,α + (N (1)
,ξα T + N

,ξα T ,i ) + [N (1) T ,α + N (1i) T ,iα + N

(2)

,ξα T + N (2i) ,ξα T ,i + N (2ij) ,ξα T ,ij ] + 2 (N (2) T ,α + N (2i) T ,iα + N (2ij) T ,ijα ) + 0( 3 ), ( 27
)
θ ,3 = T ,3 + [N (1) T ] ,3 + [N (1i) T ,i ] ,3 + 2 [N (2) T ] ,3 + [N (2i) T ,i ] ,3 + [N (2ij) T ,ij ] ,3 + 0( 3 ).( 28
)
Substituting Eqs. ( 27) and ( 28) into Eq. ( 12) and requiring the coefficient associated with -1 to be equal to zero, we get

K αβ N (1) ,ξ β ,ξα T + K αβ N (1k) ,ξ β + K αk ,ξα T ,k = 0. ( 29
)
In a similar way, by introducing ( 27) and ( 28) into ( 14) 2 and letting the coefficient associated with -1 be equal to zero, yields

K αβ N (1)
,ξ β γ ,ξα T + K αβ N (1k) ,ξ β + K αk γ ,ξα T ,k = 0. ( 30
)
It can be deduced from ( 29) and ( 30) that

K αβ N (1) ,ξ β ,ξα = 0, ξ ∈ S, ξ / ∈ L; (31)
K αβ N (1k) ,ξ β + K αk ,ξα = 0, ξ ∈ S, ξ / ∈ L; (32)
K αβ N (1)
,ξ β γ ,ξα = 0, ξ ∈ L;

K αβ N (1k) ,ξ β + K αk γ ,ξα = 0, ξ ∈ L. (33)
Above, k = 1, 2 or 3, and the surface S(y 3 ) and the line L(y 3 ) are defined by

S(y 3 ) = Π(y 3 ) ∩ Y, L(y 3 ) = S(y 3 ) ∩ Γ. (35)
In addition, the continuity condition (14) 1 of the temperature field across the interface Γ and the periodic condition demand that

N (1) = 0, N (1k) = 0, ξ ∈ L; (36)
N (1) (ξ 1 , 0, y 3 ) = N (1) (ξ 1 , 1, y 3 ), N (1) (0, ξ 2 , y 3 ) = N (1) (1, ξ 2 , y 3 ); (37)
N (1k) (ξ 1 , 0, y 3 ) = N (1k) (ξ 1 , 1, y 3 ), N (1k) (0, ξ 2 , y 3 ) = N (1k) (1, ξ 2 , y 3 ). ( 38
)
Remark that Eqs. ( 31)- [START_REF] Kerner | The elastic and thermo-elastic properties of composite media[END_REF] constitute the governing equations of a localization problem. The solution of this problem will provide the first-order localization matrix N (1) and N (1k) with k = 1, 2 and 3. It can be shown that the function N (1) corresponds to an integration constant in ξ 1 and ξ 2 , so that N (1)

,ξ β = 0 in [START_REF] Bonfoh | Anisotropic thermal conductivity of composites with ellipsoidal inclusions and highly conducting interfaces[END_REF], [START_REF] Gu | Closed-form estimates for the effective conductivity of isotropic composites with spherical particles and general imperfect interfaces[END_REF] and [START_REF] Hashin | Analysis of composite materials -a survey[END_REF]. Moreover, unlike what happens in the case where the rough interface Γ oscillates periodically only in one direction (see Le-Quang et al. [START_REF] Le-Quang | Multiscale homogenization of elastic layered composites with unidirectionally periodic rough interfaces[END_REF]), the expressions of N (1k) cannot be analytically determined in the present situation where Γ oscillates periodically in two directions. Therefore, a numerical method or an approximate analytical scheme is needed to calculate or estimate N (1k) .

Next, by substituting ( 27) and ( 28) into the energy conservation equation [START_REF] Milton | The Theory of Composites[END_REF] and requiring the coefficient associated with 0 to be equal to zero, we have

K βα N (2) ,ξα + K β3 N (1) ,3 ,ξ β T + K βα N (2k) ,ξα + K βk N (1) + K β3 N (1k) ,3 ,ξ β T ,k +    1 h 1 h 2 h 3 ∂ h 1 h 2 (K 3α N (1k) ,ξα + K 3k ) ∂y 3 + 1 h 1 h 2 h 3 ∂ h 2 h 3 (K 1α N (1k) ,ξα + K 1k ) ∂y 1 + 1 h 1 h 2 h 3 ∂ h 1 h 3 (K 2α N (1k) ,ξα + K 2k ) ∂y 2    T ,k + K βα N (2ij) ,ξα + K βj N (1i) ,ξ β + K ji + K jβ N (1i) ,ξ β T ,ij + r = ζ ∂T ∂t . (39)
In this equation, the scalar localization functions N (2) , N (2k) and N (2ij) depend only of ξ 1 , ξ 2 and y 3 , and they are periodic with period 1 in both ξ 1 and ξ 2 . They are determined through the energy conservation equation and continuous conditions. By substituting ( 27) and ( 28) into ( 14) 2 and requiring the coefficient associated to 0 to be zero, it follows that

K βα N (2)
,ξα + K β3 N

(1) ,3

= 0,

K βα N (2ij) ,ξα + K βj N (1i) = 0, γ ,ξ β K βα N (2k) ,ξα + K βk N (1) + K β3 N (1k) ,3 = h -1 3 K 3α N (1k) ,ξα + K 3k . (40)
By integrating equation [START_REF] Der Poel | On the rheology of concentrated suspension[END_REF] over the surface S of the cell Y while accounting for [START_REF] Smith | Correction and extension of van der poel's method for calculating the shear modulus of a particulate composite[END_REF], and by applying convergence theorem, we arrive at

1 h 1 h 2 h 3    ∂ h 1 h 2 K 3α N (1k) ,ξα + K 3k ∂y 3 S + ∂ h 2 h 3 (K 1α N 1k ,ξα + K 1k ) ∂y 1 S + ∂ h 1 h 3 (K 2α N (1k) ,ξα + K 2k ) ∂y 2 S    T ,k + K ji + K jβ N (1i) ,ξ β S T ,ij + |S| -1 L h -1 3 (K 3α N (1k) ,ξα + K 3k ) dl |∇γ(ξ)| T ,k + r S = ζ S ∂T ∂t (41)
where |∇γ(ξ)| = γ 2 ,ξ1 + γ 2 ,ξ2 and • S denoting the area average of a quantity (•) over S is defined by

• S = |S| -1 S (•)dS = |S| -1 1 0 1 0 (•)dξ 1 dξ 2 (42)
with |S| being the area of S. By using Lemma A.1 proven by Kristensson [3], it can be shown that

∂ h 1 h 2 (K 3k + K 3α N (1k) ,ξα ) ∂y 3 S = ∂ h 1 h 2 K 3k + K 3α N (1k) ,ξα S ∂y 3 - h 1 h 2 h 3 |S| L h -1 3 (K 3α N (1k) ,ξα + K 3k ) dl |∇γ(ξ)| . (43)
Finally, combining [START_REF] Christensen | Solutions for effective shear properties in three phase sphere and cylinder models[END_REF] with [START_REF] Rayleigh | On the influence of obstacles arranged in rectangular order upon the properties of a medium[END_REF] delivers

[K (c) ik T ,k ] ,i + r S = ζ S ∂T ∂t (44)
where

K (c) ik (y 3 ) = K ik + K iβ N (1k) ,ξ β S . (45)
Thus, at the mesoscopic scale, the rough interface zone ω (c) can be replaced with an equivalent layer or interphase, whose effective thermal conductivity tensor is determined by [START_REF] Hashin | The elastic moduli of fiber-reinforced materials[END_REF]. It is seen from ( 45) that the equivalent layer is curvilinearly anisotropic. It is important to note that the equivalent layer characterized by its thermal conductivity tensor K (c) is homogeneous along the y 1 and y 2 -directions but generally heterogeneous in the y 3 -direction.

Unlike the simple and particular case where the interface oscillates only in one direction, it is in general impossible to analytically determine the exact expression for the scalar function N (1k) in the case of an interface oscillating in two directions. To solve this problem, either an approximate method is needed first to estimate N (1k) and then to calculate K (c) ik (y 3 ) of corrugated interface zone or a numerical method can be used to numerically compute N (1k) and K (c) ik (y 3 ). The calculations of N (1k) and K (c) ik (y 3 ) by using a numerical method will be presented in following paragraph; the estimations of N (1k) and K (c) ik (y 3 ) by analytical and semi-analytical methods for a simple case will be presented in the next section of this paper.

General numerical method for determining the effective thermal conductivity tensor of a rough interface zone By taking the same origin 0 as in the Cartesian coordinate system (x 1 , x 2 , x 3 ) and by using the curvilinear coordinate system (y 1 , y 2 , y 3 ), we introduce a new Cartesian coordinate system (ζ 1 , ζ 2 , ζ 3 ) such that

ζ 1 = h 1 η 1 ξ 1 = h 1 y 1 , ζ 2 = h 2 η 2 ξ 2 = h 2 y 2 , ζ 3 = h 3 y 3 . (46)
Relative to the new Cartesian coordinate system (ζ 1 , ζ 2 , ζ 3 ), the interface Γ is transformed into Γ defined by

Γ = ζ = ζ(ζ 1 , ζ 2 , ζ 3 ) ∈ R 3 | ζ 3 = h 3 γ ζ 1 h 1 η 1 , ζ 2 h 2 η 2 = γ (ζ 1 , ζ 2 ) (47)
and the smooth surface Π(y 3 ) into a planar surface Π (ζ 3 ) characterized by

Π (ζ 3 ) = ζ = ζ(ζ 1 , ζ 2 , ζ 3 ) | ζ 3 = ζ 3 , - h 3 δ 2 ≤ ζ 3 ≤ h 3 δ 2 . ( 48
)
For a given value of ζ 3 ∈ -h3δ 2 , h3δ 2 , we consider now a fictitious three-dimensional domain M defined as

M = Π (ζ 3 ) × -h3δ 2 , h3δ 2 .
From this definition, it is clear that this domain M exhibits a columnar microstructure which is invariant in the ζ 3 -direction, but periodic in ζ 1 -and ζ 2 -directions with periods h 1 η 1 and h 2 η 2 . As a consequence, the study on M can be focused only on a unit cell Y of M defined, for example, by Y = [0,

h 1 η 1 ] × [0, h 2 η 2 ] × -h3δ 2 , h3δ 2 .
In addition, the boundary conditions prescribed on the surface ∂Y of Y are linear in ζ 3 but periodic in ζ 1 and ζ 2 . Owing to the geometry of Y , given a mesoscopic intensity vector E, we are interested in the temperature solution field over Y of the following form

θ(ζ 1 , ζ 2 , ζ 3 ) = -E α ζ α -E 3 ζ 3 -χ α (ζ 1 , ζ 2 )E α -χ 3 (ζ 1 , ζ 2 )E 3 ( 49
)
where

χ k (ζ 1 , ζ 2 )
is the localization function to be determined. On the right-hand side of Eq. ( 49), the two first terms correspond to the overall part of the temperature field while the two last ones are relative to its fluctuation part. The resulting intensity field components are then determined by

e α = E α + χ k,ζα E k , e 3 = E 3 ( 50
)
where the subscript α following a comma denotes derivative with respective to ζ α . The localization functions χ k (ξ 1 , ξ 2 ) must verify the equilibrium equation in the case of absence of heat source,

K αβ χ k,ζ β + K αk ,ζα = 0, ζ ∈ S ≡ Π (ζ 3 ) ∩ Y , ζ / ∈ L ≡ S ∩ Γ , (51)
and the continuity and periodicity conditions,

χ k (ζ 1 , ζ 2 ) = 0, ζ ∈ L , (52) K αβ χ k,ζ β + K αk γ ,ζα = 0, ζ ∈ L , (53)
χ k (h 1 η 1 , ζ 2 ) = χ k (0, ζ 2 ), χ k (ζ 1 , 0) = χ k (ζ 1 , h 2 η 2 ), ζ ∈ ∂Y . (54)
By comparing (51)-( 54) with ( 31)- [START_REF] Kerner | The elastic and thermo-elastic properties of composite media[END_REF], it can be proven that, for a given value of y 3 = ζ 3 /h 3 , the functions N (1k) (ξ 1 , ξ 2 , y 3 ) and χ k (ζ 1 , ζ 2 ) have the same form. Consequently, the localization function N (1k) can be determined via the temperature solution field θ(ζ 1 , ζ 2 , ζ 3 ) over Y by using (49). The homogenized thermal conductivity tensor of the periodic rough interfacial zone ω (c) given by Eq. ( 45) in the curvilinear coordinate system (y 1 , y 2 , y 3 ) can be calculated by

K (c) ik (ζ 3 ) = K ik + K iβ χ k,ζ β S ( 55
)
in the Cartesian coordinate system (ζ 1 , ζ 2 , ζ 3 ). Here the bracket • S designating the area average of a quantity (•) over the planar surface S ≡ Π (ζ 3 ) ∩ Y is defined by

• S = |S | -1 S (•)dS = |S | -1 h1η1 0 h2η2 0 (•)dζ 1 dζ 2 (56)
with |S | being the area of S .

FFT-based numerical method

This paragraph presents a numerical approach based on the fast Fourier transform (FFT) together with an iterative method allowing to compute the temperature solution field θ(ζ 1 , ζ 2 , ζ 3 ) of the form (49) satisfying Eqs. ( 31)-( 38) of the boundary thermal conduction problem described above. The corresponding numerical algorithm is summarized as follows:

-Iteration i = 1 : e 1 (ζ 1 , ζ 2 ) = E, q 1 (ζ 1 , ζ 2 ) = K(ζ 1 , ζ 2 )e 1 (ζ 1 , ζ 2 ), (57)
-Iteration i > 1 :
Considering that e i (ζ 1 , ζ 2 ), and

q i (ζ 1 , ζ 2 ) are known, Computing qi (τ 1 , τ 2 ) = F q i (ζ 1 , ζ 2 ) ,
Checking the convergence test :

-the iterative process will be stopped when qi

(τ 1 , τ 2 ) -qi-1 (τ 1 , τ 2 ) qi (τ 1 , τ 2 ) < tol, -otherwise êi+1 (τ 1 , τ 2 ) = êi (τ 1 , τ 2 ) -Γ0 (τ 1 , τ 2 )q i+1 (τ 1 , τ 2 ), e i+1 (ζ 1 , ζ 2 ) = F -1 êi+1 (τ 1 , τ 2 ) , q i+1 (ζ 1 , ζ 2 ) = K(ζ 1 , ζ 2 )e i+1 (ζ 1 , ζ 2 ). ( 58
)
In this algorithm, F(•) denotes the discrete Fourier transform of a quantity • while F -1 (•) corresponds to its inverse; τ = (τ 1 , τ 2 ) stands for the discrete wave vector and Γ0 represents the Fourier transform of the Green operator associated to a reference homogeneous medium of thermal conductivity tensor K 0 ; tol is a tolerance value which is set to be equal to 0.001 in our calculations. For more details about this method, the reader can refer to 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.

Finally, it is important to emphasize that the foregoing FFT-based method to determine the temperature solution field θ(ζ 1 , ζ 2 , ζ 3 ) of the boundary thermal conduction problem characterized by ( 31)- [START_REF] Kerner | The elastic and thermo-elastic properties of composite media[END_REF] and to calculate the homogenized thermal conductivity tensor K (c) for the rough interface zone ω (c) is valid for any anisotropic constituent phases. More remarkably, for a given value of ζ 3 ∈ -h3δ 2 , h3δ 2 , the method is applicable to any periodic microstructure of the plane Π (ζ 3 ).

Second step: replacement of the homogenized interphase by an imperfect interface

Consider in this subsection the case where the constant thickness δ of the rough interface zone ω (c) , after the homogenization procedure called also equivalent interphase with homogenized thermal conductivity tensor K (c) (y 3 ), is very small. This equivalent interphase ω (c) situated between two surfaces Γ 1 and Γ 2 is now replaced with an imperfect interface Γ of zero thickness located at the middle surface of ω (c) . At the same time, the two sub-domains ω (1) and ω (2) limited by two surfaces Γ 1 and Γ 2 are extended up to the middle surface Γ of ω (c) . Correspondingly, the new sub-domains occupied by phases 1 and 2 are denoted by Ω(1) and Ω(2) , respectively. By using asymptotic expansions and by demanding that the jumps of both temperature and normal component of heat flux across two surfaces Γ 1 and Γ 2 have the same values in two configurations with and without equivalent interphase ω (c) , we can establish the interfacial jump conditions that the imperfect interface Γ has to satisfy. This approach was first proposed by Sanchez-Palencia [START_REF] Sanchez-Palencia | Comportement limite d'un problème de transmission à travers une plaque mince et faiblement conductrice[END_REF], Pham-Huy and Sanchez-Palencia [START_REF] Pham-Huy | Phénomènes de transmission à travers des couches minces de conductivité élevée[END_REF] in some particular situations, then extended and completed by Miloh and Benveniste [START_REF] Miloh | On the effective conductivity of composites with ellipsoidal inhomogeneities and highly conducting interfaces[END_REF], Hashin [START_REF] Hashin | Thin interphase/imperfect interface in conduction[END_REF] and 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] in the general case. In these works, it was shown that there are three imperfect interface models. When a linearly thermal interphase between two linearly thermal phases 1 and 2 is very thin, the Kapitza's thermal resistance imperfect interface model or the highly conducting imperfect interface model is applicable according as the thermal conductivity of the interphase is much lower or much higher than the ones of the phases 1 and 2 (see e.g. [START_REF] Le-Quang | Estimations of the effective conductivity of anisotropic multiphase composites with imperfect interfaces[END_REF][START_REF] Le-Quang | Determination of the effective conductivity of composites with spherical and spheroidal anisotropic particles and imperfect interfaces[END_REF][START_REF] Bonfoh | Micromechanical modeling of the anisotropic thermal conductivity of ellipsoidal inclusion-reinforced composite materials with weakly conducting interfaces[END_REF][START_REF] Bonfoh | Anisotropic thermal conductivity of composites with ellipsoidal inclusions and highly conducting interfaces[END_REF]). In the Kapitza's thermal resistance imperfect interface model, the temperature field is discontinuous while the normal component of the heat flux is continuous across the imperfect interface Γ . In the highly conducting imperfect interface model, the temperature is continuous while the normal component of heat flux field is discontinuous across the imperfect interface Γ . In the "intermediary" case where the thermal conductivity of the interphase is neither much lower nor much higher than those of the phases 1 and 2, a general imperfect interface model in which both the temperature and the normal heat flux component are discontinuous across the interface Γ is more appropriate (see e.g. [START_REF] Le-Quang | Estimations and bounds of the effective conductivity of composites with anisotropic inclusions and general imperfect interfaces[END_REF][START_REF] Gu | Closed-form estimates for the effective conductivity of isotropic composites with spherical particles and general imperfect interfaces[END_REF]).

In the present work, due to the fact that the value of the homogenized thermal conductivity tensor K (c) (y 3 ) given by Eq. ( 45) for equivalent interphase ω (c) is always situated between the thermal conductivity tensors K (1) and K (2) of phases 1 and 2, or equivalently to min K (1) , K (2) ≤ K (c) ≤ max K (1) , K (2) , the imperfect interface model obtained by applying the asymptotic expansion approach corresponds therefore to the general one whose the jump conditions for the temperature and the normal heat flux component have the expressions

qn = δ 2 [Q 1 ( θ(-) , q(-) n , K(c) ) + Q 1 ( θ(+) , q(+) n , K(c) )] - δ 2 [Q 1 ( θ(-) , q(-) n , K (1) ) + Q 1 ( θ(+) , q(+) n , K (2) )] + 0(δ 2 ), (59) θ
= δ 2 [P 1 ( θ(-) , q(-) n , K(c) ) + P 1 ( θ(+) , q(+) n , K(c) )] - δ 2 [P 1 ( θ(-) , q(-) n , K (1) ) + P 1 ( θ(+) , q(+) n , K (2) )] + 0(δ 2 ). ( 60
)
In these equations, K(c) represents K (c) (ζ 3 ) at the middle surface of ω (c) and takes therefore the value of K (c) (ζ 3 ) for ζ 3 = 0; • denotes the value of a quantity • in the configuration with presence of the imperfect interface Γ ; • designates the interfacial jump operator defined as • = •(+) -• (-) with •(-) and •(+) denoting the values of a quantity • evaluated at Γ on the sides of phase 1 and phase 2, respectively; ∇ s (•) and ∇ s • (•) stand for the surface gradient and divergence of •, respectively, defined as

∇ s (•) = ∇(•) • T, ∇ s • (•) = ∇(•) : T (61)
where T = In ⊗ n with n being the unit vector normal to Γ oriented from the phase 1 to the phase 2. In addition, the two operators 59) and (60) are specified by

P 1 ( θ(±) , q(±) n , K (i) ) and Q 1 ( θ(±) , q(±) n , K (i) ) in Eqs. (
P 1 ( θ(±) , q(±) n , K (i) ) = - q(±) n n • K (i) • n -s (i) • ∇ s θ(±) , (62)
Q 1 ( θ(±) , q(±) n , K (i) ) = ∇ s • (S (i) • ∇ s θ(±) ) -∇ s • (s (i) q(±) n ) (63)
with

s (i) = n • K (i) n • K (i) • n , S (i) = K (i) - (K (i) • n) ⊗ (n • K (i) ) n • K (i) • n (64)
where i = 1, 2 or c denote the quantity relative to the phase 1, phase 2 or interphase, respectively.

Third step: replacement of the inclusion with imperfect interface by an equivalent inclusion with perfect interface

As described in Section 2, the interface Γ between the inclusion and matrix is assumed to be perfect and oscillates quickly and periodically around a curved surface. This interface Γ is not necessarily closed but can be open. It has been shown that, after the first-step of the first-scale homogenization procedure, the rough interfacial zone ω (c) in which the interface Γ oscillates is replaced by an equivalent interphase whose thermal conductivity tensor K (c) is determined by Eq. ( 45). At the second step of the first-scale homogenization procedure, by applying the asymptotic expansion approach to the case where the equivalent interphase ω (c) is very thin, the equivalent interphase can be replaced by a general imperfect interface Γ whose jump conditions for the temperature and the normal heat flux component are characterized by Eqs. ( 59) and (60).

In this subsection, the inclusion phase Ω( 1) is replaced by an equivalent inclusion Ω(1) having the same shape as Ω(1) and the thermal conductivity tensor K(1) of the material forming it is determined. At the same time, the general imperfect interface Γ between the inclusion and matrix phases, Ω(1) and Ω(2) , is substituted by a perfect interface Γ across which both the temperature and the normal heat flux component are continuous. The problem of determining the unknown thermal conductivity tensor K(1) of the equivalent inclusion Ω(1) has been solved in the recent work of Nguyen et al. [START_REF] Nguyen | Generalized Hill-Mendel's lemma and equivalent inclusion method for determining the effective thermal conductivity of composites with imperfect interfaces[END_REF] by developing an approach based both on a micro-macro energy equivalence condition and the generalized Hill-Mendel's lemma. More precisely, according to the results of [START_REF] Nguyen | Generalized Hill-Mendel's lemma and equivalent inclusion method for determining the effective thermal conductivity of composites with imperfect interfaces[END_REF], by prescribing on the boundary ∂Ω of Ω a boundary condition satisfying an integral condition on ∂Ω, the unknown thermal conductivity tensor K(1) of the equivalent inclusion is determined in such a way that the following condition

Γ (+) (q (0) • n θ(2) -q(2) • nθ (0) )dS + ∂ Ω(1) \ Γ (q (0) • ν θ(1) -q(1) • νθ (0) )dS = Γ (q (0) • n θ(α) -q(α) • nθ (0) )dS + ∂ Ω(1) \ Γ (q (0) • ν θ(1) -q(1) • νθ (0) )dS (65)
holds for any E 0 or Q 0 . In Eq. (65), θ(α) and q(α) with α = 1 or 2 are the temperature and heat flux solution fields of phase α in the configuration with the imperfect interface Γ ; θ(α) and q(α) designate the counterparts in the configuration with the perfect interface Γ and equivalent inclusion Ω(1) ; θ (0) and q (0) are the temperature and heat flux solution fields for the problem in which the domain Ω consisting of a homogeneous material with thermal conductivity K (2) is subjected to the same forgoing boundary condition; in addition, q(2) = -K (2) • ∇ θ(2) , q(2) = -K (2) • ∇ θ(2) , q(1) = -K (1) • ∇ θ(1) , q(1) = -K(1) • ∇ θ(1) ; ∂ Ω(1) and ∂ Ω(1) are the boundary of Ω(1) and Ω [START_REF] Kohler | Boundary and interface problems in regions with very rough boundaries[END_REF] . The derivation of this equation is omitted here. However, the reader can refer to the work of Nguyen et al. [START_REF] Nguyen | Generalized Hill-Mendel's lemma and equivalent inclusion method for determining the effective thermal conductivity of composites with imperfect interfaces[END_REF] for more details.

In the particular case where the interface Γ is closed, i.e. ∂ Ω(1) \ Γ = ∂ Ω(1) \ Γ = ∅, and Ω is assumed to be an infinite domain, it can be shown that the condition (65) reduces to

Γ (+) (q (0) • n θ(2) -q(2) • nθ (0) )dS = E 0 • K (2) Ω(1) [K (2) • (K (2) -K(1) ) -1 • K (2) -K (2) • SEsh ] -1 dV • K (2) • E 0 (66)
when the boundary condition ( 22) is concerned and

Γ (+) (q (0) • n θ(2) -q(2) • nθ (0) )dS = Q 0 Ω(1) [K (2) • (K (2) -K(1) ) -1 • K (2) -K (2) • SEsh ] -1 dV • Q 0 (67)
when the boundary condition ( 23) is under consideration. Here SEsh denotes the Eshelby tensor field inside the inclusion Ω(1) . Specially, when the inclusion Ω(1) exhibits an ellipsoidal form, the Eshelby tensor field SEsh becomes uniform inside Ω(1) for any material anisotropy. It can be shown also that if the inclusion Ω(1) processes a cylindrical from of elliptic section, then the Eshelby tensor field SEsh is uniform inside Ω(1) when the material is orthotropic and one of its privileged directions coincides with the axial direction of the cylindrical inclusion. The expression of the Eshelby tensor SEsh of an ellipsoidal inclusion and for a general 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 (66) or (67) takes the following simple form:

1 | Ω(1) | Γ (+) (q (0) • n θ(2) -q(2) • nθ (0) )dS = E 0 • K (2) • [K (2) • (K (2) -K(1) ) -1 • K (2) -K (2) • SEsh ] -1 • K (2) • E 0 (68) 1 | Ω(1) | Γ (+) (q (0) • n θ(2) -q(2) • nθ (0) )dS = Q 0 • [K (2) • (K (2) -K(1) ) -1 • K (2) -K (2) • SEsh ] -1 • Q 0 (69)
where | Ω(1) | is the volume of Ω(1) .

3.2 Second-scale homogenization: determination of the effective conductivity of composites

Case of thin interphase

As shown in the previous section, when the equivalent interphase ω (c) is very thin, after the firstscale homogenization procedure, the inclusion phase Ω (1) and the rough matrix/inclusion interface Γ can be replaced with the equivalent inclusion Ω(1) of the same form as Ω(1) and of the thermal conductivity tensor K(1) determined by ( 65), (66), ( 67), (68) or (69). At the same time, the imperfect interface Γ is substituted by a perfect interface Γ . Next, owing to the fact that the inclusion/matrix interface is now perfect, we can therefore apply any classical homogenization schemes to estimate the effective thermal conductivity tensor of the composite under consideration. For this reason, this section, relative to the second-scale homogenization procedure, consists in obtaining the closed-form expressions for the effective conductivity tensor by using some well-known classical estimation schemes such as dilute, Mori-Tanaka, self-consistent and differential approximation ones. For more details about these schemes, the reader can refer to the review papers [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] Norris | A differential scheme for the effective moduli of composites[END_REF]. Denoting by K DD , K MT , K SC and K DA the effective thermal conductivity tensors derived from the dilute, Mori-Tanaka, self-consistent and differential approximation schemes, respectively, we have the following expressions for them:

-Dilute distribution (DD) model

K DD = K (2) + c1 [ K(1) -K (2) ] • [I + SEsh • H (2) • ( K(1) -K (2) )] -1 (70)
-Mori-Tanaka (MT) model

K M T = K (1) -c2 [K (1) -K (2) ] • {c 1 [I + SEsh • H (2) • ( K(1) -K (2) )] -1 + c2 I} -1 (71) -Self-consistent (SC) model K SC = K (2) + c1 [ K(1) -K (2) ] • [I + SEsh • H SC • ( K(1) -K SC )] -1 (72)
-Differential approximation (DA) model

DK Dκ = c1 1 -c1 κ [ K(1) -K] • [I + SEsh • H • ( K(1) -K)] -1 (73)
with K (0) = K (2) and K (DA) = K (1) .

In these equations, I stands for the identity second-order tensor;

H (2) = (K (2) ) -1 , H (SC) = (K (SC)
) -1 and H = K -1 denote the thermal resistivity second-order tensors; c1 and c2 are the volume fractions of the inclusion and matrix phases, respectively.

Case of thick interphase

For the case of thick equivalent interphase ω (c) , after the first-scale homogenization procedure, we have obtained at mesoscopic scale a three-phase inclusion/interphase/matrix model. At the macroscopic scale, after calculating the macroscopic heat flux and intensity vectors Q and E by Eq. ( 18), ( 22) or [START_REF] Sanchez-Palencia | Comportement limite d'un problème de transmission à travers une plaque mince et faiblement conductrice[END_REF], the effective thermal conductivity tensor K eff will be directly identified by using Eq.( 17).

Homogenization of layered composites with rough interfaces oscillating around a planar surface

In this section, as a first example of application of the general homogenization procedure elaborated above, we consider a layered composite Ω consisting of two layers Ω (1) and Ω (2) whose interface Γ is perfect and oscillates at the microscopic scale quickly and periodically in both the x 1 -and x 2 -directions between the parallel planes Γ 1 and Γ 2 defined by x 3 = -δ 2 and x 3 = δ 2 (see Fig. 4). For simplicity, the oscillating interface Γ in the x 1 -and x 2 -directions are assumed to have the same period, i.e. 1 = 2 = and η 1 = η 2 = 1. In addition, Ω (1) and Ω (2) are assumed to be made of two isotropic materials whose thermal conductivity tensors are given by K (1) = k 1 I and K (2) = k 2 I with k 1 /k 2 = 10. In this example of application, owing to the fact that the interface Γ oscillates in two directions around a planar surface, the orthogonal curvilinear coordinates (y 1 , y 2 , y 3 ) used to describe the interface Γ coincide with the Cartesian ones (x 1 , x 2 , x 3 ), i.e. y 1 = x 1 , y 2 = x 2 , y 3 = x 3 and h 1 = h 2 = h 3 = 1. The rough interface Γ is assumed to possess a cone profile (see Fig. 4). For a given value of y 3 ∈ -δ 2 , δ 2 , the line L(y 3 ) = S(y 3 ) ∩ Γ with S(y 3 ) = Π(y 3 ) ∩ Y becomes a circle of radius r(y 3 ) given by

r(y 3 ) = R δ 2 -y 3 δ ( 74
)
where δ is the thickness of the corrugated zone ω (c) and R is the radius of the bottom of a cone. In the present example, the value of R is chosen to be such that R/ = 0.48039. In addition, the thicknesses of the two layers ω (1) and ω (2) defined by Eq. ( 7) are denoted by δ 1 and δ 2 , respectively.

ζ 1 = y 1 , ζ 2 = y 2 , ζ 3 = y 3 . ( 75
)
For a given value of ζ 3 ∈ -δ 2 , δ 2 , the unit cell Y has a square prism form whose base is a square of unit length and whose height is equal to δ. Moreover, this square prism medium consists of phase 2 in which a circular cylinder of radius r(

ζ 3 ) = R( δ 2 -ζ3) δ
made of phase 1 is embedded. Owing to the isotropy of the two constituent materials, the columnar microstructure of Y and the boundary conditions prescribed on the surface ∂Y of Y with mesoscopic intensity tensor E in such a way that they are periodic in the plane ζ 1 -ζ 2 but linear in ζ 3 -direction, it can be shown that the temperature solution field over Y takes the form of (49) with

χ 3 (ζ 1 , ζ 2 ) = 0. ( 76
)
The

f 1 (ζ 3 ) = πr 2 (ζ 3 ) = πR 2 δ 2 -ζ 3 2 δ 2 2 = π 0.240195 -0.48039 ζ 3 δ 2 . ( 77
)
Generalized self-consistent estimate scheme In parallel to the forgoing FFT-based numerical method, to analytically estimate the homogenized thermal conductivity matrix components K (c) ik (ζ 3 ) of the rough interfacial zone ω (c) as defined by Eq. (55), an analytical approach based on the generalized self-consistent scheme (GSCS) is now proposed. This scheme was presented in the first time by Kerner [START_REF] Kerner | The elastic and thermo-elastic properties of composite media[END_REF], then improved and completed by Van der Poel [START_REF] Der Poel | On the rheology of concentrated suspension[END_REF], Smith [START_REF] Smith | Correction and extension of van der poel's method for calculating the shear modulus of a particulate composite[END_REF][START_REF] Smith | Correction and extension of van der poel's method for calculating the shear modulus of a particulate composite[END_REF] and Christensen and Lo [START_REF] Christensen | Solutions for effective shear properties in three phase sphere and cylinder models[END_REF].

For a given value of ζ 3 ∈ -δ 2 , δ 2 , the GSCS is constructed by considering a double-coated cylinder of height δ, composed of a core of phase 1 with conductivity k 1 coated by a layer of phase 2 of conductivity k 2 , embedded in a homogeneous and transversely isotropic with respect to the axis ζ 3 medium M * of the same height δ but infinite in the transverse plane ζ 1 -ζ 2 . The radii of the core and the coating, denoted by r 1 and r 2 , are chosen in such a way that they are compatible with

ζ 3 δ/2 δ/2 r 1 r 2 K (1) K (2) M* ζ 2 ζ 1 ( , ,δ/2) ζ ζ 1 2 (ζ ζ , ,δ/2) 1 2 ( , ,δ/2) ζ ζ 1 2 O ( , ,-δ/2) ζ ζ 1 2 ( , ,-δ/2) ζ ζ 1 2 ( , ,-δ/2) ζ ζ 1 2
Fig. 5: Sketch of the generalized self-consistent scheme their phase volume fractions, i.e. (r 1 /r 2 ) 2 = f 1 = 1 -f 2 with f 1 given by Eq. (77) (Fig. 5). The external homogeneous medium M * surrounding the composite cylinder and consisting of the effective homogeneous material with unknown transversely isotropic thermal conductivity tensor is subjected to the following uniform intensity boundary conditions on its lateral surface ∂M * :

θ(ζ) = E 1 ζ 1 + E 2 ζ 2 , ζ ∈ ∂M * (78)
and to the thermal insulated boundary condition on the top and bottom surfaces of M * . According to the GSCS, the unknown thermal conductivity tensor of the effective homogeneous material can be determined by requiring that the presence of the circular composite cylinder do not change the initial energy without the circular composite cylinder. By taking into account the expression (49) for the temperature solution field, this requirement can be shown to be equivalent to

χ k,ζ β S = f 1 χ k,ζ β S 1 + f 2 χ k,ζ β S 2 = 0 ( 79
)
where • S 1 and • S 2 denote the area averages of • over the phases 1 and 2, respectively. By substituting Eq. (79) into Eq. ( 55), the expression of the homogenized thermal conductivity tensor K (c) (ζ 3 ) can be recast to

K (c) ik (ζ 3 ) = f 1 k 1 ∆ ik + f 2 k 2 ∆ ik + f 1 (k 1 -k 2 )B ik (80)
where ∆ ik is the Kronecker symbol and B ik are the components of the second-order tensor B defined by

B ik = χ i,ζ k S 1 . (81)
) that K (c) 13 (ζ 3 ) = K (c) 31 (ζ 3 ) = K (c) 23 (ζ 3 ) = K (c) 32 (ζ 3 ) = 0 and K (c) 33 (ζ 3 ) = f 1 k 1 + f 2 k 2 . 80
Under the uniform intensity boundary conditions (78), we seek the temperature solution field in the following form:

θ(ζ) =          C i1 r1 r 2 + C i2 (E 1 ζ 1 + E 2 ζ 2 ), for 0 ≤ r = ζ 2 1 + ζ 2 2 ≤ r 1 C m1 r1 r 2 + C m2 (E 1 ζ 1 + E 2 ζ 2 ), for r 1 ≤ r = ζ 2 1 + ζ 2 2 ≤ r 2 C e1 r1 r 2 + C e2 (E 1 ζ 1 + E 2 ζ 2 ), for r 2 ≤ r = ζ 2 1 + ζ 2 2 ( 82
)
where C i1 , C i2 , C m1 , C m2 , C e1 and C e2 are constants to be determined from the boundary, interface inclusion/matrix condition together with a condition avoiding the temperature singularity in the core of the composite cylinder and a self-consistency condition that the energies with and without composite cylinder have the same value. We obtain from these conditions the following expressions of the constants:

C i1 = C e1 = 0, C e2 = 1, C i2 = 2k 2 k 2 (1 + f 1 ) + k 1 (1 -f 1 ) , (83)
C m1 = (k 2 -k 1 ) k 2 (1 + f 1 ) + k 1 (1 -f 1 ) , C m2 = (k 2 + k 1 ) k 2 (1 + f 1 ) + k 1 (1 -f 1 ) . ( 84
)
Combining Eqs. ( 82)-(84) with Eqs. ( 49) and (81) yields the expressions of the non-zero components of the tensor B as follows:

B 11 = B 22 = - (k 1 -k 2 )(1 -f 1 ) k 2 (1 + f 1 ) + k 1 (1 -f 1 ) . ( 85
)
Finally, introducing the expressions of B given by ( 85) into (80) gives rise to the following formulas for the non-zero components of the homogenized thermal conductivity tensor K (c) (ζ 3 ) of the rough interfacial zone

K (c) 11 (ζ 3 ) = K (c) 22 (ζ 3 ) = f 1 k 1 + f 2 k 2 - (k 1 -k 2 ) 2 f 1 f 2 k 2 (1 + f 1 ) + k 1 (1 -f 1 ) , K (c) 33 (ζ 3 ) = f 1 k 1 + f 2 k 2 . ( 86
)
We recall that f 1 = 1 -f 2 in Eq. ( 86) is determined by Eq. (77).

To summarize, by applying the GSCS, for a given value of ζ 3 ∈ -δ 2 , δ 2 , all the components of the homogenized thermal conductivity tensor K (c) (ζ 3 ) of the rough interfacial zone are obtained in an explicit analytical way. The numerical values obtained by applying the GSCS for the homogenized thermal conductivity tensor components K (c) ik (ζ 3 ) are plotted in Fig. 6 and compared with the ones provided by the FFT-based numerical method.

Chen and Kuo's method In addition to the analytical results obtained by applying the GSCS and the numerical ones derived by using the FFT-based numerical method, in this subsection, we use a semi-analytical method proposed by Chen and Kuo [START_REF] Chen | Transport properties of composites consisting of periodic arrays of exponentially graded cylinders with cylindrically orthotropic materials[END_REF] to solve a boundary value problem of heat conduction. In this method, a unit cell consisting of a circular cylindrical inclusion of height δ and of radius

r 1 = R( δ 2 -ζ3) δ
coated by a square prism matrix with height δ and square base of unit length is subjected to a periodic boundary conditions in-plane ζ 1 -ζ 2 . Taking account the periodicity of the composite in question, the following boundary conditions are adopted:

θ | ζ1=-1/2 = 1 2 E 1 , θ | ζ1=1/2 = - 1 2 E 1 , ∂θ ∂ζ 2 | ζ2=±1/2 = 0 ( 87
)
where E 1 represents the macroscopic intensity along the ζ 1 -direction. The symmetry of unit cell yields:

θ(ζ 1 , ζ 2 , ζ 3 ) = θ(ζ 1 , -ζ 2 , ζ 3 ), θ(ζ 1 , ζ 2 , ζ 3 ) -θ(0, 0, ζ 3 ) = -[θ(-ζ 1 , ζ 2 , ζ 3 ) -θ(0, 0, ζ 3 )] . (88)
Relative to a cylindrical coordinates system (r, φ, ζ 3 ), by applying the conservation energy equation ∇ • [K(ζ)∇θ] = 0 to an isotropic material and by assuming that the temperature field is independent of ζ 3 , we have

∂ 2 θ ∂r 2 + 1 r ∂θ ∂r + 1 r 2 ∂ 2 θ ∂φ 2 = 0. ( 89
)
The solution to this equation takes the following general form

θ (1) = C 0 + ∞ n=1 C 2n-1 r (2n-1) cos(2n -1)φ ( 90
)
in the inclusion phase and

θ (2) = A 0 + ∞ n=1 A 2n-1 r 2n-1 + B 2n-1 r -2n+1 cos(2n -1)φ. ( 91
)
in the matrix phase. In these expressions, the coefficients A 2n-1 , B 2n-1 and C 2n-1 are unknown constants to be determined while the coefficients A 0 and C 0 representing the reference values of the temperature fields in the matrix and inclusion phases are set to be equal to zero, i.e. A 0 = C 0 = 0. Since the continuity condition of the heat flux normal component at the matrix/inclusion interface r = r 1 , we obtain

A 2n-1 = k 2 + k 1 k 2 -k 1 r 2(1-2n) 1 B 2n-1 , C 2n-1 = k 2 + k 1 k 2 -k 1 + 1 r 2(1-2n) 1 B 2n-1 . (92)
The temperature fields in the inclusion and matrix phases, θ (1) and θ (2) , must also satisfy the periodic boundary conditions (87). In order to achieve this, the consistency condition of Rayleigh's identity [START_REF] Rayleigh | On the influence of obstacles arranged in rectangular order upon the properties of a medium[END_REF] is used in this work. According to this identity, the consistency condition is equivalent to

∞ n=1 A 2n-1 r 2n-1 cos(2n -1)φ = -E 1 ζ 1 + j =0 ∞ n=1 B 2n-1 r -2n+1 j cos(2n -1)φ j . (93)
Here (r, φ) are the polar coordinates for a point measured from the centre of the central cylinder while (r j , φ j ) correspond to the polar coordinates for the same point but measured from the centre of the j-th cylinder with index j running over all cylinders except the one at the origin. By applying the method proposed by Nicorovici et al. [START_REF] Nicorovici | Transport properties of a three-phase composite material: the square array of coated cylinders[END_REF] in which we replace first r j and φ j in the right-hand of (93) with r and φ, we equate then the (2n -1)-th partial derivative with respect to ζ 1 on both sides of (93) at the point corresponding to centre of the central cylinder, and finally obtain

A 2n-1 + ∞ m=1 2m + 2n -3 2n -1 Σ 2m+2n-2 B 2m-1 = -E 1 ∆ n1 , (94)
where Σ m = j =0 r -m j cos mφ j with (r j , φ j ) stands for the polar coordinates for the central cylinder centre measured from the centre of the j-th cylinder and m n denotes the binomial coefficient. By substituting Eq. (92) 1 into Eq. (94), we get a system of linear equations which can be rewritten in the following matrix form:

([T] + [W]){X} = {F} (95)
where [T] is a diagonal matrix whose diagonal components are equal to (k 2 + k 1 )/(k 2 -k 1 ), {F} is a column vector with only one non-zero component F 1 = -1, {X} and [W] are a column vector and a square matrix, respectively, whose components are specified by

X 2m-1 = B 2m-1 E 1 r 2(2m-1) 1 , W mn = (2m + 2n -3)!Σ 2m+2n-2 r 2(2m-1) i (2m -2)!(2n -1)! . (96)
First solving Eq.( 96) to obtain B 2n-1 and then invoking Eq. (92) to calculate A 2n-1 and C 2n-1 , we derive through Eq. (90) the temperature field inside the inclusion phase. It is important to note from (90) that, unlike the case of GSCS where the intensity field is uniform inside the inclusion, the resulting intensity field obtained in this case with Chen and Kuo's method is not uniform inside the inclusion. Correspondingly, by using Eq. ( 49), the scalar localization function χ 1 (ζ 1 , ζ 2 ) can be determined. In addition, since the symmetry of the unit cell in ζ 1 -and

ζ 2 -directions, it is immediate that χ 1 (ζ 1 , ζ 2 ) = χ 2 (ζ 1 , ζ 2 )
. By substituting these functions into Eqs. ( 81) and (80), we obtain the values of

K (c) 11 (ζ 3 ) = K (c) 22 (ζ 3
). Moreover, with the temperature field inside the inclusion expressed as in Eq. (90), it can be shown with the help of (81) and (80) that B 12 = B 21 = 0 and

K (c) 12 (ζ 3 ) = K (c)
21 (ζ 3 ) = 0. Now, we show, in Figs. 6 and7, the variation of the non-zero normalized homogenized thermal conductivity matrix components K (c) ij /k 2 of the rough interfacial layer ω (c) in the terms of both the ratio ζ 3 /δ and the area fraction f 1 of phase 1 that is defined by (77). In these figures, the homogenized thermal conductivity tensor components are not only obtained by using FFT-based numerical method but also derived by applying the GSCS and Chen and Kuo's methods. It can be observed from Figs. 6 and 7 that the values of the homogenized thermal conductivity matrix components provided by FFTbased numerical method are very close to the ones derived by Chen and Kuo's method and have a good agreement with the ones obtained by GSCS when the area fraction of phase 1 is lower than 0.5. However, when the area fraction of phase 1 is greater than 0.5, the values of the homogenized thermal conductivity matrix components obtained by FFT-based numerical method and Chen and Kuo's method are slightly different with the ones provided by GSCS. This is because, the intensity solution field inside the inclusion obtained by FFT-based numerical method is generally not uniform inside the inclusion as in the case of Chen and Kuo's method while the intensity solution field becomes uniform in the case of GSCS. Furthermore, by comparing the values obtained for K

(c) (FFT) = K 33 (c) (CK) = K 33 (c) (GSCS) K 11 (c) (GSCS) K 11 (c) (CK) K 11 (c) (FFT)
Fig. 6: Normalized values of the homogenized thermal conductivity matrix components for the rough interfacial zone, obtained by GSCS, Chen and Kuo's method and FFT; comparisons with the Reuss and Voigt bounds.

Replacement of the homogenized interphase with an imperfect interface

Following the process presented in section 3.1.2, when the rough interfacial layer ω (c) is very thin, the interphase layer obtained by applying the homogenization process can be replaced with an imperfect interface Γ of zero thickness situated at the middle surface of ω (c) and characterized by y 3 = 0, or equivalently to x 3 = 0. Immediately, the two layers ω (1) and ω (2) occupied by two phases 1 and 2 are prolonged up to the middle surface Γ . The two corresponding extended subdomains obtained from ω (1) and ω (2) become Ω(1) and Ω(2) . The jump conditions (59) and (60) for both temperature and the normal heat flux component across the imperfect interface Γ can be recast into

q3 = δ 2 ( K(c) 11 -k 1 ) ∂ 2 θ(-) ∂x 2 1 + ( K(c) 22 -k 1 ) ∂ 2 θ(-) ∂x 2 2 +( K(c) 11 -k 2 ) ∂ 2 θ(+) ∂x 2 1 + ( K(c) 22 -k 2 ) ∂ 2 θ(+) ∂x 2 2 + 0(δ 2 ), (97) θ
= δ 2 1 k 1 - 1 K(c) 33 q(-) 3 + 1 k 2 - 1 K(c) 33 q(+) 3 + 0(δ 2 ). ( 98
)
Here, we recall that K(c) ij takes the value of K In this problem, the two-phase layered composite Ω with imperfect interface Γ is assumed to be subjected to the uniform in-plane intensity boundary condition [START_REF] Moulinec | A fast numerical method for computing the linear and nonlinear mechanical properties of composites[END_REF] with E 0 1 = 0, E 0 2 = 0 but E 0 3 = 0. Under this boundary condition, with the help of the imperfect interface conditions (97) and (98), the temperature, intensity and heat flux fields in the phases 1 and 2 are given by θ

(1) (x) = θ(2) (x) = -E 0 1 x 1 -E 0 2 x 2 , x ∈ Ω(1) ∪ Ω(2) , (99) ê
(1) (x) = E 0 1 j 1 + E 0 2 j 2 , q(1) (x) = k 1 (E 0 1 j 1 + E 0 2 j 2 ), x ∈ Ω(1) , (100) ê
(2) (x) = E 0 1 j 1 + E 0 2 j 2 , q(2) (x) = k 2 (E 0 1 j 1 + E 0 2 j 2 ), x ∈ Ω(2) . (101)
It is interesting to remark from Eqs. ( 99)-( 101) that, even with the presence of the imperfect interface Γ between phases 1 and 2, 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. Problem 2: Out-plane thermal conduction

The second problem concerns the two-phase layered composite Ω with imperfect interface Γ undergoing the uniform out-plane heat flux boundary condition [START_REF] Moulinec | A numerical method for computing the overall response of nonlinear composites with complex microstructure[END_REF] with Q 0 1 = 0, Q 0 2 = 0 but Q 0 3 = 0. Under this boundary condition, by taking into account the imperfect interface conditions (97) and (98), the temperature, intensity and heat flux fields in phases 1 and 2 take the following form:

θ(1) (x) = - Q 0 3 k 1 x 3 - δQ 0 3 2 1 k 1 + 1 k 2 - 2 K(c) 33 , ê(1) (x) = Q 0 3 k 1 j 3 , q(1) (x) = Q 0 3 j 3 , x ∈ Ω(1) , (102) θ(2) (x) = - Q 0 3 k 2 x 3 + θ ref , ê(2) (x) = Q 0 3 k 2 j 3 , q(2) (x) = Q 0 3 j 3 , x ∈ Ω(2) . ( 103
)
Similarly again to the classical case of layered composite with perfect interface, we notice from Eqs. ( 102) and (103) that, notwithstanding the presence of the imperfect interface Γ , the heat flux component field in the direction normal to the layer plane of the layered composite is uniform and identical to the one prescribed on its boundary.

Replacement of the layer with imperfect interface by an equivalent layer with perfect interface

Relative to the third step of the first-scale homogenization procedure, the layer Ω( 1) is now replaced with an equivalent layer Ω(1) processing the same shape as Ω(1) and the thermal conductivity tensor of the form K(1) = K(1)

11 j 1 ⊗ j 1 + K(1) 22 j 2 ⊗ j 2 + K(1) 33 j 3 ⊗ j 3 in which K(1) 11 , K (1)
22 and K(1) 33 are unknown. In addition, the general imperfect interface Γ between two layers Ω(1) and Ω(2) is substituted by a perfect interface Γ across which both the temperature field and the normal component of heat flux field are continuous.

In order to determine the unknown thermal conductivity tensor component K(1) 11 , K(1) 22 and K(1)

33

of K(1) of the equivalent layer, the condition ( 65) is now applied. The temperature field θ(α) and the heat flux field q(α) are obtained by using Eqs. ( 99)-(103) while the temperature field θ(α) and the heat flux field q(α) in the classical layered composite with perfect interface as well as the temperature field θ (0) and the heat flux field q (0) in the homogeneous layered composite of conductivity tensor K (2) , are as follows: Problem 1: In-plane thermal conduction

θ(1) (x) = θ(2) (x) = θ (0) (x) = -E 0 1 x 1 -E 0 2 x 2 , x ∈ Ω(1) ∪ Ω(2) , (104)
e (0) (x) = E 0 1 j 1 + E 0 2 j 2 , q (0) (x) = k 2 (E 0 1 j 1 + E 0 2 j 2 ), x ∈ Ω(1) ∪ Ω(2) , (105) ẽ
(1) (x) = E 0 1 j 1 + E 0 2 j 2 , q(1) (x) = K(1) 11 E 0 1 j 1 + K(1) 22 E 0 2 j 2 , x ∈ Ω(1) , (106) ẽ
(2) (x) = E 0 1 j 1 + E 0 2 j 2 , q(2) (x) = k 2 (E 0 1 j 1 + E 0 2 j 2 ), x ∈ Ω(2) . ( 107
)
Problem 2: Out-plane thermal conduction

θ (0) (x) = - Q 0 3 k 2 x 3 , e (0) (x) = Q 0 3 k 2 j 3 , q (0) (x) = Q 0 3 j 3 , x ∈ Ω(1) ∪ Ω(2) , (108) θ
(1) (x) = - Q 0 3 K(1) 33 x 3 , ẽ(1) (x) = Q 0 3 K(1) 33 j 3 , q(1) (x) = Q 0 3 j 3 , x ∈ Ω(1) , (109) θ
(2) (x) = - Q 0 3 k 2 x 3 , ẽ(2) (x) = Q 0 3 k 2 j 3 , q(2) (x) = Q 0 3 j 3 , x ∈ Ω(2) . (110)
Finally, we obtain, from the condition (65), the following expressions for the non-zero thermal conductivity matrix components K(1) ij of the equivalent layer Ω(1) :

K(1) 11 = K(1) 22 = k 1 , K(1) 33 = 2k 1 k 2 K(c) 33 δ 1 + δ 2 2δk 1 k 2 + 2 K(c) 33 k 2 δ 1 + δ 2 -δk 2 K(c) 33 -δk 1 K(c) 33 . (111)
It is interesting to remark from Eq. ( 111) that the in-plane thermal conductivity K(1)

11 = K(1)
22 of the equivalent layer Ω(1) is exactly equal to the one of Ω (1) . In other words, the in-plane thermal conductivity K(1)

11 = K(1)
22 of the equivalent layer Ω(1) does not depend on the homogenized thermal conductivity tensor K (c) as well as the thickness δ of the rough interface zone ω (c) . However, unlike the in-plane thermal conductivity K(1)

11 = K(1)
22 , the thermal conductivity in the layering direction K(1) 33 of the equivalent layer Ω(1) given by Eq. ( 111) is a function of both the homogenized thermal conductivity K(c) [START_REF] Gu | Closed-form estimates for the effective conductivity of isotropic composites with spherical particles and general imperfect interfaces[END_REF] and the thickness ratio δ/δ 1 of the rough interface zone ω (c) .

4.2 Second-scale homogenization procedure: homogenization of layered composite

Case of thin interphase

When the equivalent interphase is very thin, after replacing the layer Ω(1) made of phase 1 with imperfect interface Γ by an equivalent layer Ω(1) with perfect interface Γ , the layered composite under consideration becomes therefore a laminated composite consisting of two layers Ω(1) and Ω(2) whose the interface Γ is now perfect. As a consequence, this two-layered composite can be homogenized by applying the classical theory for laminates (see e.g. [START_REF] Milton | The Theory of Composites[END_REF]). The non-zero macroscopic thermal conductivity tensor components can be exactly and analytically determined by

K eff 11 = K eff 22 = c1 k 1 + c2 k 2 , K eff 33 = c1 k 1 + c2 k 2 + cc 1 K(c) 33 - 1 2k 1 - 1 2k 2 -1 (112)
where c1 , c2 and cc corresponding to the volume fractions of Ω(1) , Ω(2) and ω (c) are defined by

c1 = δ 1 + δ 2 δ 1 + δ 2 + δ , c2 = δ 2 + δ 2 δ 1 + δ 2 + δ , cc = δ δ 1 + δ 2 + δ . ( 113
)
For later use, the expression (112) of K eff 11 , K eff 22 and K eff 33 can be recast into the following form:

K eff 11 = K eff 22 = c1 k 1 + c2 k 2 + k 2 δ F k (k 1 /k 2 , c1 , c2 , δ/ ) + 0(δ 2 ), K eff 33 = (H eff 33 ) -1 = c1 h 1 + c2 h 2 + h 2 δ F h (k 1 /k 2 , c1 , c2 , δ/ ) -1 + 0(δ 2 ) ( 114
)
where h 1 = 1/k 1 and h 2 = 1/k 2 denote the thermal resistivities of phases 1 and 2, respectively, and H eff 33 = (K eff 33 ) -1 designates the effective thermal resistivity of the layered composite along the layering direction; the two scalar dimensionless functions F k (k 1 /k 2 , c1 , c2 , δ/ ) and F h (k 1 /k 2 , c1 , c2 , δ/ ), that depend on the contrast k 1 /k 2 and the volume fractions c1 and c2 of phases 1 and 2 as well as on the ratio δ/ , are called fluctuating parts of the effective thermal conductivity K eff 11 and resistivity H eff 33 , respectively.

Case of thick interphase

When the equivalent interphase obtained from the homogenization process of the rough interfacial zone is not sufficiently thin, the replacement of the homogenized interphase with an imperfect interface as well as the replacement of the layer made of phase 1 with imperfect interface by an equivalent layer with perfect interface are no longer to apply. For this reason, the direct method based on the definition [START_REF] Persson | The Homogenization Method: An Introduction[END_REF] of the effective thermal conductivity tensor K eff for three-layered (phase 1/equivalent interphase/phase 2) composite together with the definitions (18) of the macroscopic intensity vector E and heat flux vector Q is now used to determine the effective thermal conductivities. It can be shown that the corresponding effective thermal conductivity tensor components can be calculated by (see also [START_REF] Milton | The Theory of Composites[END_REF] for more details)

K eff 33 = K -1 33 -1 Ω , (115)
K eff 22 = K 22 Ω + K 23 K -1 33 2 Ω K -1 33 -1 Ω -K 2 23 K -1 33 Ω , (116)
K eff 11 = K 11 Ω + K 13 K -1 33 2 Ω K -1 33 -1 Ω -K 2 13 K -1 33 Ω , (117)
K eff 12 = K eff 21 = K 12 Ω + K 13 K -1 33 Ω K -1 33 -1 Ω K -1 33 K 23 Ω -K 13 K -1 33 K 23 Ω , (118)
K eff 13 = K eff 31 = K 13 K -1 33 Ω K -1 33 -1 Ω , (119)
K eff 23 = K eff 32 = K -1 33 -1 Ω K -1 33 K 23 Ω (120)
where the average operator • Ω over the layered composite domain Ω is defined by

• Ω = 1 δ + δ 1 + δ 2 -δ 2 -δ 2 -δ1
• (1)

dx 3 + δ 2 -δ 2 • (c) dx 3 + δ 2 +δ2 δ 2 • (2) dx 3 . (121)
Next, in order to numerically illustrate the approach developed as well as the results obtained in sections 3 and 4 and to study the influence of interface zone to the effective thermal properties of the two-phase laminated composite under consideration, we consider the case where the thickness δ of the corrugated zone ω (c) is chosen to be equal to the thickness δ 1 and δ 2 of two layers ω (1) and ω (2) , i.e. δ = δ 1 = δ 2 , while the period of rough interface is set to be vary from δ/10 to δ. The normalized values of the effective thermal conductivity components K eff 11 /k 2 , K eff 22 /k 2 and K eff 33 /k 2 of the laminated composite are plotted in Fig. 9. In the latter, the effective thermal conductivity components K eff 11 , K eff 22 and K eff 33 are calculated by applying the analytical formula ( 115)-( 121) wherein the thermal effective conductivity of the rough interface zone is determined by using FFT, GSCS and the method of Chen and Kuo. Then, the values obtained for K eff 11 , K eff 22 and K eff 33 are compared with the ones provided by employing the FEM with Matlab software. The computation of the effective thermal conductivity with FEM is carried out over a unit cell with periodic boundary conditions in the plane x 1 -x 2 .

It can be seen from Fig. 9 that: (i) the values of the effective thermal conductivities K eff 11 and K eff 33 obtained by using either the numerical method FFT or the semi-analytical method proposed by Chen and Kuo almost coincide with the approximate ones derived by GSCS. Furthermore, these values of K eff 11 and K eff 33 are very close to the corresponding ones computed by FEM which will be considered as benchmarks in our study; (ii) with the same volume fractions of phases 1 and 2, when the ratio δ/ increases, the values of the effective thermal conductivity K eff 11 decrease significantly while the values of the effective thermal conductivity K eff 33 decrease very slightly; (iii) the higher the ratio δ/ , closer the values of the effective thermal conductivities K eff 11 and K eff 33 obtained by FFT, CK method or GSCS are to the ones provided by FEM. This means that the numerical and approximate values obtained for versus the ratio δ 1 /δ = δ 2 /δ for a two-phase layered composite with rough interface of cone profile the effective thermal conductivities K eff 33 with FFT, CK method or GSCS are more precise when the interface between two layers of the laminated composite is very rough.

Concerning the case where the equivalent interphase in which the rough interface oscillates is very thin, we consider the second example where the thickness δ 1 and δ 2 of two layers ω (1) and ω (2) take the same value but the thickness δ of the corrugated zone ω (c) is chosen in such a way that the ratio δ 1 /δ = δ 2 /δ varies from 10 to 100. At the same time, the value of the ratio /δ is set to be equal to 0.1 to guarantee that the interface between two phases is very rough. The fluctuating functions F k and F h of the effective thermal conductivity K eff 11 and resistivity H eff 33 , defined by Eq. ( 114), are computed and plotted in Figs. 10 and 11 in terms of the ratio δ 1 /δ = δ 2 /δ. It can be observed from Figs. 10 and11 that the values of the fluctuating functions F k and F h obtained by FFT, CK method and GSCS are very close to the ones provided by FEM. In particular, when the value of δ 1 /δ = δ 2 /δ is lager enough, or equivalently when the corrugated zone is very thin, the values of the fluctuating functions F k and F h calculated by FFT, CK method, GSCS and FEM converge to the one derived from the two-scale homogenization approach described in subsection 4.2.1.

Homogenization of two-phase composite cylinder with a rough interface oscillating around a circumferential surface

The second example of application is related to a composite cylinder consisting of a host matrix phase Ω (2) in which inclusions Ω (1) are inserted. The inclusion/matrix interface Γ is assumed to be perfect and oscillates quickly and periodically at the microscopic scale around a cylindrical surface in the longitudinal and angular directions (see Fig. 12). As a consequence, the curvilinear coordinate system (y 1 , y 2 , y 3 ) used in section 2 to describe the rough interface Γ corresponds to a cylindrical coordinate system (r, φ, z), i.e. Fig. 11: Fluctuating part of the effective thermal resistivity H eff 33 versus the ratio δ 1 /δ = δ 2 /δ for a two-phase layered composite with rough interface of cone profile where R denotes the radius of the circumferential surface around which the interface Γ oscillates. Relative to the cylindrical coordinate system (r, φ, z), the interface Γ between Ω (1) and Ω (2) is described by

y 1 = z = x 3 , y 2 = Rφ = R arctan(x 2 /x 1 ), y 3 = r -R = x 2 1 + x 2 2 -R (122)
Γ = y = y(z, φ, r) ∈ R 3 | r -R = γ(ξ 1 , ξ 2 ), ξ 1 = z 1 , ξ 2 = Rφ 2 . (123)
Above, γ(ξ 1 , ξ 2 ) is defined by γ(ξ 1 , ξ 2 ) = δ 2 sin(2πξ 2 ) sin(2πξ 1 ) with δ denoting the thickness of the rough interfacial zone; 1 and 2 standing for the periods of Γ along the longitudinal and angular directions are given by 1 = η 1 and 2 = η 2 with = H, η 1 = 1, η 2 = Rφ 0 /H and φ 0 = 2π/N . Here, N denotes the number of oscillations in the angular direction. In addition, the metric coefficients defined by (2) are given by

h 1 = 1, h 2 = r R , h 3 = 1. (124)
5.1 First-scale homogenization of a rough interface

FFT-based numerical method for determining the effective thermal conductivity tensor of a rough interfacial zone

As presented in section 3.1.1, to numerically compute the scalar localization function N (1k) and the homogenized thermal conductivity matrix components K (c) ik (y 3 ) of the periodic rough interfacial zone ω (c) , a linear transformation from the coordinate system (y 1 , y 2 , y 3 ) into the Cartesian coordinate system (ζ 1 , ζ 2 , ζ 3 ) is realized as follows: 125), the rough interface oscillates therefore around a planar surface (see Fig. 13).

ζ 1 = z H , ζ 2 = rφ H , ζ 3 = r -R. (125)
Thus, for a given value of (1) and Ω (2) , the columnar microstructure of Y and the periodic boundary condition in-plane ζ 1 -ζ 2 and linear boundary condition in ζ 3 -direction, it can be proven that the temperature solution field over Y takes the form of (49) with

ζ 3 ∈ -δ 2 , δ 2 corresponding to a value of r = R + ζ 3 ∈ R -δ 2 , R + δ
χ 3 (ζ 1 , ζ 2 ) = 0. (126)
The calculation of the two remaining scalar localization functions χ 1 (ζ 1 , ζ 2 ) and χ 2 (ζ 1 , ζ 2 ) that allow us to determine the homogenized conductivity matrix components K (c) ik (ζ 3 ) can be carried out by applying the FFT-based numerical method as presented in section 3.1.1. Recall that the thermal conductivity tensor K (c) of the equivalent interphase obtained by homogenizing the rough interface zone ω (c) is cylindrically anisotropic and heterogeneous in the radial direction but homogeneous in the azimuthal and vertical directions.

Finally, by keeping the ratio k 1 /k 2 constant, for example k 1 /k 2 = 10, the values obtained for the non-zero homogenized conductivity matrix components K

Replacement of the homogenized interphase by an imperfect interface

When the rough interfacial layer ω (c) is very thin, as mentioned in section 3.1.2, the first emplacement process in which the equivalent or interphase layer obtained by applying the homogenization process is substituted by an imperfect interface Γ of zero thickness situated at the middle surface of ω (c) specified by y 3 = 0, or equivalently to r = R, is now realized. The two subdomains ω (1) and ω (2) made of two phases 1 and 2 are enlarged to the middle surface Γ and denoted now by Ω(1) and Ω(2) , respectively.

Relative to the cylindrical coordinate system, the jump conditions (59) and (60) established for the temperature and normal heat flux component across the imperfect interface Γ can be rewritten in the following simple form:

Ratio ζ 3 /δ -0.5 -0.4 -0.3 -0.2 -0.1 0 0.1 0.2 0.3 0.4 0.5 K (c) ij /k
qr = δ 2 ( K(c) φφ -k 1 ) R 2 ∂ 2 θ(-) ∂φ 2 + ( K(c) zz -k 1 ) ∂ 2 θ(-) ∂z 2 + ( K(c) φφ -k 2 ) R 2 ∂ 2 θ(+) ∂φ 2 + ( K(c) zz -k 2 ) ∂ 2 θ(-) ∂z 2 + 0(δ 2 ), (127) θ
= δ 2 1 k 1 - 1 K(c) rr q(-) r + 1 k 2 - 1 K(c) rr q(+) r + 0(δ 2 ) ( 128
)
where K(c) rr , K(c) φφ and K(c) zz take the values of

K (c) rr (ζ 3 ), K (c)
φφ (ζ 3 ) and K (c) zz (ζ 3 ) when ζ 3 = 0, or equiva- lently to r = R.
In the following paragraphs, we consider two fundamental thermal conduction problems about the composite cylinder:

Problem 1: In-plane thermal conduction

The first problem is concerned with a two-phase composite cylinder Ω consisting of a cylindrical inclusion phase Ω (1) embedded via the imperfect interface Γ into a cylindrical matrix phase Ω (2) which is assumed to be infinite in the plane x 1 -x 2 . The composite cylinder Ω is subjected to the mixed boundary conditions [START_REF] Bonnet | Effective properties of elastic periodic composite media with fibers[END_REF] in which a uniform in-plane intensity boundary condition with E 0 1 = 0, E 0 2 = 0 but E 0 3 = 0 is prescribed on the lateral surface of Ω and a thermal insulated boundary condition is imposed on the bottom and top bases of Ω. Accounting for these boundary conditions together with the imperfect interface conditions (127) and (128), the temperature, intensity and heat flux fields in phases 1 and 2 are given by θ

(1) (x) = -Â1 (E 0 1 r cos φ + E 0 2 r sin φ), x ∈ Ω(1) , (129) θ(2)
(x) = -Â2 + B2 R r 2 (E 0 1 r cos φ + E 0 2 r sin φ), x ∈ Ω(2) , (130) ê(1)
(x) = Â1 (E 0 1 cos φ + E 0 2 sin φ)f r + Â1 (-E 0 1 sin φ + E 0 2 cos φ)f φ , x ∈ Ω(1) , (131) ê(2)
(x) = Â2 + B2 R r 2 (E 0 1 cos φ + E 0 2 sin φ)f r + Â2 -B2 R r 2 (-E 0 1 sin φ + E 0 2 cos φ)f φ , x ∈ Ω(2) , (132) q
(1) (x) = k 1 Â1 (E 0 1 cos φ + E 0 2 sin φ)f r + k 1 Â1 (-E 0 1 sin φ + E 0 2 cos φ)f φ , x ∈ Ω(1) , (133) q
(2) (x) = k 2 Â2 + B2 R r 2 (E 0 1 cos φ + E 0 2 sin φ)f r + k 2 Â2 -B2 R r 2 (-E 0 1 sin φ + E 0 2 cos φ)f φ , x ∈ Ω(2) , (134)
where three constants Â1 , Â2 and B2 are given by

Â2 = 1, B2 = B21 B22 , Â1 = Â11 Â12 (135) with B21
= 4 K(c) rr (k 1 -k 2 ) + 4 δ R ( K(c) rr K(c) φφ -k 1 k 2 ) + δ R 2 ( K(c) rr - K(c) φφ )(k 2 -k 1 ), (136)
B22 = Â12 = 4(k 1 + k 2 ) K(c) rr + 4 δ R (k 1 k 2 + K(c) rr K(c) φφ -k 1 K(c) rr -k 2 K(c) rr ) + δ R 2 (k 1 + k 2 )( K(c) rr + K(c) φφ ) -2(k 1 k 2 + K(c) rr K(c) φφ ) , (137)
Â11 = 8 K(c) rr k 2 + 2 δ R 2 k 2 2 -k 2 ( K(c) rr + K(c) φφ ) + K(c) rr K(c) φφ . (138)
Problem 2: Axial thermal conduction In the second problem, mixed boundary conditions [START_REF] Bonnet | Effective properties of elastic periodic composite media with fibers[END_REF], i.e. a uniform axial intensity boundary condition with E 0 1 = E 0 2 = 0 but E 0 3 = 0 prescribed on the bottom and top bases of Ω and a thermal insulated boundary condition applied on the lateral surface of Ω, are now imposed. By accounting for this boundary condition and the imperfect interface conditions (127) and (128), the temperature, intensity and heat flux fields in phases 1 and 2 are provided by

θ(1) (x) = -E 0 3 z, ê(1) (x) = E 0 3 f z , q(1) (x) = k 1 E 0 3 f z , x ∈ Ω(1) , (139) θ(2) (x) = -E 0 3 z, ê(2) (x) = E 0 3 f z , q(2) (x) = k 2 E 0 3 f z , x ∈ Ω(2) . (140)

Replacement of the circular cylindrical core with imperfect interface by an equivalent circular cylinder with perfect interface

The circular cylindrical core Ω(1) with imperfect interface Γ is replaced with an equivalent core Ω(1) of the same shape as Ω(1) and unknown thermal conductivity tensor K(1) = K(1)

11 j 1 ⊗ j 1 + K(1) 22 j 2 ⊗ j 2 + K(1)
33 j 3 ⊗ j 3 . Owing to the geometry of Ω(1) , it is clear that K(1) 11 = K(1) 22 . At the same time, the general imperfect interface Γ is now replaced with a perfect interface Γ across it both the temperature and the normal heat flux component are continuous. The determination of the unknown thermal conductivity tensor K(1) of the equivalent layer is carried out by using the condition (65) together with the temperature field θ(α) and the heat flux field q(α) given by Eqs. ( 129)-(149) while the temperature field θ(α) and the heat flux field q(α) in the case of perfect interface are provided by: Problem 1: In-plane thermal conduction

θ(1) (x) = -Ã1 (E 0 1 r cos φ + E 0 2 r sin φ), x ∈ Ω(1) , (141) θ
(2) (x) = -Ã2 + B2 R r 2 (E 0 1 r cos φ + E 0 2 r sin φ), x ∈ Ω(2) , (142) ẽ
(1) (x) = Ã1 (E 0 1 cos φ + E 0 2 sin φ)f r + Ã1 (-E 0 1 sin φ + E 0 2 cos φ)f φ , x ∈ Ω(1) , (143) ẽ
(2) (x) = Ã2 + B2 R r 2 (E 0 1 cos φ + E 0 2 sin φ)f r + Ã2 -B2 R r 2 (-E 0 1 sin φ + E 0 2 cos φ)f φ , x ∈ Ω(2) , (144) q
(1) (x) = K(1) 11 Ã1 (E 0 1 cos φ + E 0 2 sin φ)f r + K(1) 11 Ã1 (-E 0 1 sin φ + E 0 2 cos φ)f φ , x ∈ Ω(1) , (145) q(2) (x) = k 2 Ã2 + B2 R r 2 (E 0 1 cos φ + E 0 2 sin φ)f r + k 2 Ã2 -B2 R r 2 (-E 0 1 sin φ + E 0 2 cos φ)f φ , x ∈ Ω(2) , (146)
where three constants Ã1 , Ã2 and B2 are given by

Ã2 = 1, B2 = K(1) 11 -k 2 K(1) 11 + k 2 , Ã1 = 2k 2 K(1) 11 + k 2 . ( 147
)
Problem 2: Axial thermal conduction

θ(1) (x) = -E 0 3 z, ẽ(1) (x) = E 0 3 f z , q(1) (x) = K(1) 33 E 0 3 f z x ∈ Ω(1) , (148) θ(2) (x) = -E 0 3 z, ẽ(2) (x) = E 0 3 f z , q(2) (x) = k 2 E 0 3 f z , x ∈ Ω(2) . (149)
By applying condition (65) or (66) and (67), the expressions for the non-zero thermal conductivity tensor components K(1) ij of the equivalent layer Ω(1) are obtained as follows:

K(1) 33 = k 1 , K(1) 11 = K(1) 22 = 4k 1 k 2 K(c) rr + 2k 2 K(c) rr δ R (2 K(c) φφ -k 1 -k 2 ) -(k 1 -K(c) rr )(k 2 - K(c) φφ )k 2 δ R 2 × 4k 2 K(c) rr + 2 δ R (2k 1 k 2 -k 1 K(c) rr -k 2 K(c) rr ) -(k 1 - K(c) φφ )(k 2 -K(c) rr ) δ R 2 -1 . (150)
5.2 Second-scale homogenization of composite cylinder

Case of thin interphase

If the equivalent interphase is very thin, after replacing the circular cylindrical core consisting of phase 1 and the imperfect interface Γ by an equivalent circular cylindrical inclusion of conductivities given by (150) and with perfect interface Γ , the composite cylinder under consideration consists now of a matrix phase containing a circular cylindrical inclusion with perfect matrix/inclusion interface. Correspondingly, Eshelby's tensor for a circular cylindrical inclusion takes the following simple form:

SEsh = 1 2 (j 1 ⊗ j 1 + j 2 ⊗ j 2 ). (151)
The expressions (70)-( 73) for non-zero effective thermal conductivity tensor components obtained from the dilute, Mori-Tanaka, self-consistent and differential approximation schemes reduce to:

-Dilute distribution (DD) model

K DD 11 = K DD 22 = k 2 + 2c 1 k 2 ( K(1) 11 -k 2 ) K(1) 11 + k 2 , K DD 33 = c1 K(1) 33 + c2 k 2 ; (152)
-Mori-Tanaka (MT) model

K MT 11 = K MT 22 = k 2 + 2c 1 k 2 ( K(1) 11 -k 2 ) 2k 2 + c2 ( K(1) 11 -k 2 ) , K MT 33 = c1 K(1) 33 + c2 k 2 ; (153)
-Self-consistent (SC) model

K SC 11 = k 2 + 2c 1 K SC 11 ( K(1) 11 -k 2 ) K(1) 11 + K SC 11 , K SC 22 = K SC 11 , K SC 33 = c1 K(1) 33 + c2 k 2 (154)
where the effective thermal conductivity K SC 11 is calculated as the real positive root of (154) 1 ; c1 and c2 are the volume fractions of Ω(1) and Ω(2) , respectively; -Differential-approximation (DA) model

K DA 33 = c1 K(1) 33 + c2 k 2 (155)
and the effective thermal conductivities K DA 11 and K DA 22 are given as K DA 11 = K DA 22 = K(1) where K(κ) with 0 ≤ κ ≤ 1 is the solution to the following differential equation

DK Dκ = 1 1 -κc 1 2c 1 ( K(1) 11 -K) K(1) 11 + K with K(0) = k 2 . ( 156
)
It can also be shown that the effective thermal conductivities K DA 11 and K DA 22 correspond to the real solution root of the following second-order equation

k 2 ( K(1) 11 -K) 2 -K( K(1) 11 -k 2 ) 2 (1 -c1 ) 2 = 0. ( 157
)
It is very interesting to demonstrate that, when the equivalent cylindrical inclusions of circular section are assumed to be randomly distributed in a host matrix phase and the effective thermal behaviour of the composite at the macroscopic scale is supposed to be transversely isotropic, the effective thermal conductivity obtained by Eq.( 153) with Mori-Tanaka model coincides exactly with the one derived by applying the generalized self-consistent scheme (GSCS) or coated cylinder assemblage (CCA) model. Recall that the GSCS has been described in section 4.1.1. The CCA model, proposed for the first time by Hashin and Rosen [START_REF] Hashin | The elastic moduli of fiber-reinforced materials[END_REF] and considered as the two-dimensional version of the well-known coated sphere assemblage (CSA) of Hashin [START_REF] Hashin | The elastic moduli of heterogeneous materials[END_REF], will be described in details in the next section.

By considering the thin interphase case in which the thickness δ of the interphase is chosen to be such that δ/R = 0.001, we plot in Fig. 16 the variation of the non-zero normalized effective thermal conductivity tensor components k eff 11 /k 2 , k eff 22 /k 2 and k eff 33 /k 2 in terms of the equivalent inclusion volume fraction of a two-phase composite cylinder whose rough interface between these two phases oscillates around a circumferential surface with k 1 /k 2 = 10, H/R = 0.1 and φ 0 = π/20. The determination of the effective thermal conductivities of this composite is carried out by combining the equivalent inclusion method (EIM) with the DD, MT (equivalent to GSCS et CCA), DA and SC models. In order to validate the values obtained for these effective thermal conductivities, we compare in Fig. 16 them with the corresponding Voigt and Reuss bounds. It can be observed in Fig. 16 that the effective axial thermal conductivity k eff 33 of the composite cylinder coincides exactly with its Reuss bound. Except the DD model, the values of the effective in-plane thermal conductivity obtained by SC, MT (or equivalently to GSCS et CCA) and DA models respect well the Voigt and Reuss bounds. Moreover, in the recent work of Nguyen et al. [START_REF] Nguyen | Generalized Hill-Mendel's lemma and equivalent inclusion method for determining the effective thermal conductivity of composites with imperfect interfaces[END_REF], it is shown that the MT (or GSCS or CCA) model is the best approach to estimating the effective in-plane thermal conductivity of periodic composites where cylinder inclusions are either squarely or hexagonally distributed into a host matrix phase. However, for composite cylinders with random distribution of cylinder inclusions, the best estimation for the effective in-plane thermal conductivity is the average value of the ones obtained by MT (equivalent to GSCS or CCA) and differential approximation models.

Case of thick interphase

In the case of thick interphase, the coated cylinder assemblage (CCA) model introduced by Hashin and Rosen [START_REF] Hashin | The elastic moduli of fiber-reinforced materials[END_REF] is now applied to compute the effective macroscopic conductivity of the composite (see Fig. 17). It is important to notice that the effective thermal conductivity obtained by applying CCA model coincides exactly with the one provided by using the GSCS or MT model. The CCA model is built by successively and completely filling up the domain Ω, consisting of a homogeneous and transversely isotropic material of unknown effective macroscopic conductivity tensor K eff = K eff 11 j 1 ⊗ j 1 + K eff 22 j 2 ⊗ j 2 + K eff 33 j 3 ⊗ j 3 with K eff 11 = K eff 22 , with double-coated circular inclusions of all sizes in such a way that: (i) each double-coated cylinder inclusion in the assemblage is identical up to a scale factor to any other double-coated cylinder inclusion; (ii) the presence of each double-coated cylinder inclusion in Ω does not disturb the original fields in the outside medium. Each double-coated cylinder inclusion in the assemblage is composed of a circular cylinder core coated by two concentric layers. The core of the double-coated circular inclusion is made of phase 1 (inclusion phase) whereas the outer coating layer is formed of phase 2 (matrix phase) and the inner coating layer (equivalent interphase c) consists of the rough interfacial zone that has been homogenized and replaced with an equivalent interphase. The radii of the core and of the inner and outer coating layers 1 and 2, denoted by ρ 1 , ρ c and ρ 2 , are chosen in such a way that they are compatible with the area fractions of the core, inner and outer coating layers, symbolized by c 1 , c c and c 2 , and specified as follows:

c 1 = ρ 2 1 ρ 2 2 , c c = ρ 2 c -ρ 2 1 ρ 2 2 , c 2 = ρ 2 2 -ρ 2 c ρ 2 2 . ( 158
)
ρ 1 Ω ρ 2 ρ c K (1) K (2) K ( ) c
Fig. 17: Sketch of the coated cylinder assemblage (CCA) model By using the first-scale homogenization procedure shown in section 5.1.1, the effective mesoscopic thermal conductivity tensor K (c) (r) of the equivalent layer ω (c) can be numerically determined for the case where the interface Γ between the matrix phase and each inclusion is periodically rough along the circumferential and axial directions. The effective mesoscopic thermal conductivity tensor K (c) (r) of the equivalent layer is cylindrically anisotropic and heterogeneous in the radial direction but homogeneous in the circumferential and axial directions. More precisely, relative to the cylindrical coordinate system, K (c) (r) takes the following form:

K (c) (r) = K (c) rr (r)f r ⊗ f r + K (c) φφ (r)f φ ⊗ f φ + K (c) zz (r)f z ⊗ f z (159)
Since all double-coated circular cylindrically inclusions in the CCA are similar, the volume average fields over each double-coated circular cylindrically inclusion are equal to the ones over any other doublecoated circular cylindrically inclusion. Consequently, the volume average fields over the composite under consideration Ω are also identical to the ones calculated in each double-coated circular cylindrically inclusion. Thus, each double-coated circular cylindrically inclusion acts as the whole composite sample, and can be taken to be a representative volume element (RVE). For this reason, we consider now a typical RVE with

ρ 1 = R - δ 2 , ρ c = R + δ 2 . ( 160
)
Here we recall that R represents the radius of the cylindrical surface around which the rough interface oscillates, and δ is the thickness of the equivalent interphase. Problem 1: In-plane thermal conduction Under the mixed boundary conditions (21) composed of a uniform in-plane intensity boundary one with E 0 1 = 0 but E 0 2 = 0 and E 0 3 = 0 prescribed on the lateral surface and a thermal insulated boundary one applied on the bottom and top surfaces of Ω, by taking into account the geometric and material symmetries, we seek the temperature solution fields in the core, inner and outer layer coatings of a double-coated circular cylindrical inclusion in the following form

θ (i) (x) = f (i) (r) cos φ (161)
where i = 1, c or 2 and f (i) (r) is a scalar function while the temperature solution field outside the double-coated cylindrically circular inclusion θ (e) (x) is given by

θ (e) (x) = -E 0 1 r cos φ. (162)
By using the energy conservation equation in the case of the stationary thermal conduction without heat source, it can be shown that: (i) the scalar functions f (1) (r) and f (2) (r) for the core (phase 1) and outer layer (phase 2) of the double-coated circular inclusion take the following general form

f (1) (r) = A 1 r, f (2) (r) = A 2 r + B 2 r -1 (163)
where A 1 , A 2 and B 2 are unknown constants to be determined from the boundary and interface conditions; (ii) the function f (c) (r) of the equivalent interphase must satisfy the following equation

d dr K (c) rr df (c) dr + K (c) rr r df (c) dr -f (c) K (c) φφ r 2 = 0. (164)
In general, it is impossible to obtain an analytical and exact solution for the differential equation (164). For this reason, a numerical method is used to find its solution. A numerical approach based on the finite difference method is proposed in this work to solve Eq. ( 164). According to this method, we define first

∆r = δ N iter , r i = i∆r + ρ 1 (165)
where i = 0, 1, ..., N iter and N iter is the number of iterations chosen in such a way that it is large enough to guarantee the convergence of computations. By using an approximation as follows

df (c) (r) dr r=ri f (c) i+1 -f (c) i-1 2∆r , d 2 f (c) (r) d 2 r r=ri f (c) i+1 -2f (c) i + f (c) i (∆r) 2 (166)
where f (c) i = f (c) (r i ), Eq. ( 164) can be rewritten in the following form:

K (c) rr (r j ) f (c) j+1 -2f (c) j + f (c) j-1 (∆r) 2 - K (c) φφ (r j ) r 2 j f (c) j + K (c) rr (r j+1 ) -K (c) rr (r j-1 ) 2∆r + K (c) rr (r j ) r j f (c) j+1 -f (c) j-1 2∆r = 0 (167)
with j = 1, 2, ..., N iter -1 and r 0 = ρ 1 , r Niter = ρ c . The perfect bonding across the interfaces at r = ρ 1 , r = ρ c and r = ρ 2 implies the following continuity conditions for the temperature field and the normal component of the heat flux vector:

-At r = ρ 1 : A 1 ρ 1 = f (c) 0 , A 1 k 1 = K (c) rr (ρ 1 ) f (c) 1 -f (c) 0 ∆r . ( 168
) -At r = ρ c : f (c) Niter = A 2 ρ c + B 2 ρ -1 c , K (c) rr (ρ 2 ) f (c) Niter -f (c) Niter-1 ∆r = k 2 (A 2 -B 2 ρ -2 c ). (169)
-At r = ρ 2 :

A 2 + B 2 ρ -2 2 = E 0 1 , k 2 (A 2 -B 2 ρ -2 2 ) = K eff 11 E 0 1 . (170)
In these equations, K eff 11 is the unknown effective macroscopic conductivity of the composite cylinder under consideration.

Finally, Eqs. ( 167)-(170) constitute a system of N iter + 5 homogeneous linear equations for

N iter + 5 unknowns A 1 , A 2 , B 2 , f (c) 0 , f (c) 1 ,..., f (c)
Niter and K eff 11 which can be recast in the matrix form:

[Y] {X} = 0 (171)
where [Y] is a (N iter + 5) × (N iter + 5) matrix and {X} is a vector defined by

{X} = A 1 , f (c) 0 , f (c) 1 , ..., f (c) Niter , A 2 , B 2 , E 0 1 T . (172)
A non-trivial solution to this system exists if and only if the determinant of the relevant (N iter + 5) × (N iter + 5) matrix [Y] is equal to zero, or equivalently to

F (K eff 11 ) = det[Y] = 0. ( 173
)
This necessary and sufficient condition yields a linear equation for K eff 11 , symbolized by F (K eff 11 ) = 0. This equation allows us to determine the effective macroscopic conductivity K eff 11 as

K eff 11 = F (0) F (0) -F (1) . ( 174
)
When the effective macroscopic conductivity K eff 11 has been obtained,

N iter + 4 remaining unknown constants A 1 , A 2 , B 2 , f (c) 0 , f (c) 1 ,..., f (c)
Niter are then expressed in terms of the prescribed intensity constant E 0 1 by using the system of homogeneous liner equations ( 167)-(170). Then, the temperature, intensity and heat flux fields in the core, inner and outer coating layers of the double-coated circular cylindrically inclusion can be determined.

Problem 2: Axial thermal conduction When the mixed boundary condition ( 21) is applied, i.e. a uniform axial intensity boundary condition with E 0 1 = E 0 2 = 0 but E 0 3 = 0 prescribed on the top and bottom surface of Ω and a thermal insulated boundary condition imposed on the lateral surface of Ω, we find the temperature, intensity and heat flux solution fields in the core, inner and outer layer coatings of a double-coated circular cylindrically inclusion in the form as follows:

θ (1) (x) = θ (c) (x) = θ (2) (x) = -E 0 3 x 3 , (175)
e (1)

(x) = e (c) (x) = e (2) (x) = E 0 3 j 3 , (176)
q (1) (x) = k 1 E 0 3 j 3 , q (c) (x) = K (c) zz E3 0 j 3 , q (2) (x) = k 2 E 0 3 j 3 . (177)
Owing to the fact that each double-coated circular cylindrically inclusion can be taken to be a representative volume element of Ω, we can evaluate the macroscopic intensity vector and macroscopic heat flux vector over any double-coated circular cylindrically inclusion. Combining (175)-( 177) with [START_REF] Jikov | Homogenization of differential operators and integral functionals[END_REF] implies the following expressions for macroscopic intensity vector and macroscopic heat flux vector

E = E 0 3 j 3 , Q = 2 ρ 2 2 ρ1 0 k 1 rdr + ρc ρ1 K (c) zz rdr + ρ2 ρc k 2 rdr E 0 3 j 3 . (178)
The effective macroscopic conductivity K eff 33 is therefore given by

K eff 33 = 2 ρ 2 2 ρ1 0 k 1 rdr + ρc ρ1 K (c) zz (r)rdr + ρ2 ρc k 2 rdr (179)
or equivalently by

K eff 33 = c 1 k 1 + c c K (c) zz ω (c) + c 2 k 2 (180)
Next, in order to numerically illustrate the approach developed above to the case of thick interphase, we consider a two-phase composite cylinder where the interface between two phases oscillates periodically around a circumferential surface of radius R in the longitudinal direction with a period H = R and in the angular direction with a period φ 0 = π/20. The conductivity ratio of matrix/inclusion is kept constant with k 1 /k 2 = 10 while the thickness of the rough interface zone where the interface oscillates is set to vary from 0.001R to R. The radius R of the middle surface of the rough interface zone is chosen in such a way that the ratio c1 = (R/ρ 2 ) 2 takes successively four values, i.e. c1 = (R/ρ 2 ) 2 = 0.1, 0.2, 0. It can be observed from Figs. 18-21 that, in comparison with smooth interface, the effect of rough interface on the effective thermal conductivities becomes very impor-tant when the (homogenized) interphase associated to the oscillating zone is thick, say δ/R ≥ 0.1. However, when the interphase is very thin, the effect of interface roughness on the effective thermal conductivities is insignificant. Futher, it can be also seen from Figs. 18-21 that the effect of interface roughness depends not only on the thickness of the interphase but also on the inclusion volume fraction. The larger the inclusion volume fraction, the greater the effect of interface roughness.

Finally, we compare in Figs. [START_REF] Le-Quang | Effective thermal conductivity of periodic composites with highly conducting imperfect interfaces[END_REF] and 23 the values obtained for the normalized effective in-plane and axial thermal conductivities k eff 11 /k 2 = k eff 22 /k 2 and k eff 33 /k 2 by applying the CCA model with the ones derived by using the equivalent inclusion method (EIM), MT, DD and DA models. As before, the thickness of the rough interface zone is set to vary from 0.001R to R while the conductivity ratio of matrix/inclusion k 1 /k 2 as well as the periods H and φ 0 along the longitudinal and angular directions of the oscillating interface are kept constant with k 1 /k 2 = 10, H = R and φ 0 = π/20. It can be seen from Fig. 22 that when the thickness of the rough interface zone is very small, the normalized effective in-plane thermal conductivity k eff 11 /k 2 = k eff 22 /k 2 obtained by CCA model is very close to the ones derived with EIM and MT models. This means that, for the case of thin interface, the MT model is the best one to estimate the effective in-plane thermal conductivity. It can be observed in Fig. 23 that the normalized effective axial thermal conductivity k eff 33 /k 2 take the same value when the EIM and MT, DD and DA models are used. Moreover, these values of k eff 33 /k 2 are very close to the ones obtained from CCA model.

Concluding remarks

In the present work, a two-scale homogenization method has been elaborated to solve the problem of determining the effective thermal conductivity of a composite in which the interfaces between constituent phases oscillate in two directions around a curved surface. The method proposed in this work is general and can be applied to any composite material with rough interfaces oscillating around an arbitrary curved surface along two directions. The effectiveness and validity of the proposed method are shown through two examples of application. The first one consists of a layered composite in which the interface between two neighboring layers oscillates fast about a plane and along two directions and the second one concerns a two-phase composite cylinder with rough interface oscillating around a circumferential surface in the axial and angular directions. In these two cases, thin and thick interfacial zones have been studied in this work. The results obtained for effective thermal conductivities by applying the two-scale homogenization method are shown to agree well with the numerical results provided by the finite element method (FEM) and comply with the corresponding Reuss, Voigt and Hashin-Shtrikman bounds.

Since the thermal conduction phenomenon studied in this work is mathematically similar to other transport phenomena like electric conduction, dielectrics, magnetism, diffusion and flow in porous media, the general homogenization procedure elaborated in the present work for the thermal conduction phenomenon can be directly applicable to other transport phenomena.

Finally, the extension of the results obtained in this work for thermal conduction phenomenon to elasticity and coupled multifield phenomena a complete and systematic study of the dependance of the effective properties of composites on the ratio of interphase thickness/inclusion size, the inclusion volume fraction and the roughness geometry will be carried out in a forthcoming work.

Fig. 3 :

 3 Fig. 3: Two-scale homogenization procedure

4. 1 Fig. 4 :

 14 Fig. 4: Two-phase layered composite with periodically rough interface of cone profile

  following paragraphs are therefore dedicated to compute only two remaining scalar localization functions χ 1 (ζ 1 , ζ 2 ) and χ 2 (ζ 1 , ζ 2 ) needed to calculate the homogenized conductivity tensor components K (c) ik (ζ 3 ).FFT-based numerical method By applying the numerical method based on the FFT and the iterative method which have been developed in the section 3.1.1, we compute the scalar localization functions χ 1 (ζ 1 , ζ 2 ) and χ 2 (ζ 1 , ζ 2 ) and the homogenized conductivity tensor components K (c) ik (ζ 3 ). The values obtained for K (c) ik (ζ 3 ) are plotted in Fig.6 versus the area fraction f 1 of the inclusion phase 1 in the plane ζ 1 -ζ 2 of Y , given by

Moreover, due to the fact that χ 3

 3 (ζ 1 , ζ 2 ) = 0 and the scalar functions χ 1 and χ 2 are independent of ζ 3 , it is clear from (81) that B 13 = B 23 = B 31 = B 32 = B 33 = 0. Consequently, it follows from Eq. (

  ij with the corresponding Voigt and Reuss bounds, it can be seen from Figs.6 and 7that the values of K (c) 33 computed by FFT-based numerical method, GSCS and Chen and Kuo's method coincide completely with the Voigt bound while the values of K situated between the Voigt and Reuss bounds as expected.

Ratio ζ 3 Fig. 7 :

 37 Fig.7: Normalized values of the homogenized thermal conductivity matrix components for the rough interfacial zone obtained by GSCS, Chen and Kuo's method and FFT; comparisons with its Reuss and Voigt bounds.

ij (ζ 3 )

 3 when ζ 3 = 0. In the following paragraphs, we consider two fundamental thermal conduction problems for layered composites:Problem 1: In-plane thermal conduction

Fig. 8 :Fig. 9 :Fig. 10 :

 8910 Fig. 8: Mesh of a unit cell used in FEM

Fig. 12 : 2 OFig. 13 :

 12213 Fig.12: Two-phase composite cylinder with rough interface periodically oscillating around a circumferential surface

2

 2 and relative to the Cartesian coordinate system (ζ 1 , ζ 2 , ζ 3 ), the unit cell Y obtained by a transformation from Y possesses a prism form whose base is a rectangle of dimensions 1 and rφ 0 /H in ζ 1 -and ζ 2directions, respectively, and the height is equal to δ. As an example, by setting H/R = δ/R = 0.1 and φ 0 = π/20, we show in Fig. 14, the microstructure of the base for the rectangular prism unit cell Y with three different values of ζ 3 /δ, i.e. ζ 3 /δ = 1/4, ζ 3 = -1/4 and ζ 3 /δ = 0. Due to the isotropy of two materials consisting Ω

Fig. 14 :

 14 Fig. 14: Base of the rectangular prism unit cell Y with H/R = δ/R = 0.1, φ 0 = π/20 and (a) ζ 3 /δ = 1/4; (b) ζ 3 /δ = 0; (c) ζ 3 /δ = -1/4

  φφ (ζ 3 ) and K (c) zz (ζ 3 ) associated to the cylindrical coordinate system (r, φ, z) are plotted in Fig.15 versus the values of ζ 3 /δ for different values of δ/R and φ 0 . It can be seen from Fig.15 that, for any value of the ratio ζ 3 /δ, the homogenized conductivity K (c) rr (ζ 3 ) coincides exactly with the Reuss bound while the normalized homogenized conductivities K (c) φφ (ζ 3 )/k 2 and K (c) zz (ζ 3 )/k 2 take the same values as K (c) φφ /k 2 = 3.7522 and K (c) zz /k 2 = 2.7614 when ζ 3 = 0.

Fig. 15 :

 15 Fig.15: Normalized values of the homogenized thermal conductivity matrix components for the rough interfacial zone obtained by using the FFT method when k 1 /k 2 = 10; comparisons with its Reuss and Voigt bounds for different values of δ/R and φ 0 .

33 Fig. 16 :

 3316 Fig.16: Normalized values obtained by the equivalent inclusion method (EIM) for effective thermal conductivity matrix components in terms of the volume fraction of equivalent inclusions for a two-phase composite cylinder with rough interface oscillating around a circumferential surface; comparisons with its Reuss and Voigt bounds in the case of k 1 /k 2 = 10, δ/R = 0.001, H/R = 0.1 and φ 0 = π/20.

1 k. 4 Fig. 18 : 2 K

 14182 Fig. 18: Normalized values for the in-plane effective thermal conductivity K eff 11 = K eff 22 in terms of ratio log 10 (R/δ) for a two-phase composite cylinder with k 1 /k 2 = 10 and rough interface oscillating around a circumferential surface for different values of inclusion volume fraction c1 .

3 and 0. 4 .

 4 The normalized effective in-plane thermal conductivity k eff 11 /k 2 = k eff 22 /k 2 and the normalized effective axial thermal conductivity k eff 33 /k 2 are plotted in Figs. 18 and 19 versus the logarithm of the ratio R/δ. Similarly, we show in Figs. 20 and 21 the variations of the normalized effective in-plane and axial thermal conductivities k eff 11 /k 2 = k eff 22 /k 2 and k eff 33 /k 2 in terms of the ratio c1 = (R/ρ 2 ) 2 with the values of δ/R equal to 0.001, 0.1 or 1.

4 Fig. 19 :

 419 Fig. 19: Normalized values for the axial effective thermal conductivity K eff 33 in terms of ratio log 10 (R/δ) for a two-phase composite cylinder with k 1 /k 2 = 10 and rough interface oscillating around a circumferential surface for different values of inclusion volume fraction c1 .

1 Fig. 20 : 1 Fig. 21 :

 120121 Fig. 20: Normalized values for the in-plane effective thermal conductivity K eff 11 = K eff 22 in terms of inclusion volume fraction c1 for a two-phase composite cylinder with k 1 /k 2 = 10 and rough interface oscillating around a circumferential surface for different values of ratio δ/R.

1 CCA/ 20 Fig. 22 :

 12022 Fig. 22: Comparison of the in-plane effective thermal conductivity values, K eff 11 = K eff 22 , obtained by EIM and CCA for a two-phase composite cylinder with k 1 /k 2 = 10, c1 = 0.2 and rough interface oscillating around a circumferential surface for different values of ratio δ/R.

20 Fig. 23 :

 2023 Fig. 23: Comparison of the in-plane effective thermal conductivity values, K eff 33 , obtained by EIM and CCA for a two-phase composite cylinder with k 1 /k 2 = 10, c1 = 0.2 and rough interface oscillating around a circumferential surface for different values of ratio δ/R.

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