Micro-mechanical model for the effective thermal conductivity of the multi-oriented inclusions reinforced composites with imperfect interfaces

International Journal of Heat and Mass Transfer 148 (2019) 119167

Contents lists available at ScienceDirect

International Journal of Heat and Mass Transfer journal homepage: www.elsevier.com/locate/hmt

Micro-mechanical model for the effective thermal conductivity of the multi-oriented inclusions reinforced composites with imperfect interfaces Wenlong Tian a, M.W. Fu a,∗, Lehua Qi b, Haihui Ruan a a b

Department of Mechanical Engineering, The Hong Kong Polytechnic University, Hong Kong, China School of Mechanical Engineering, Northwestern Polytechnical University, Xi’an 710072, China

a r t i c l e

i n f o

Article history: Received 29 September 2019 Revised 30 November 2019 Accepted 2 December 2019

Keywords: Imperfect interface Micro-mechanics Multi-oriented composites Thermal conductivity Two-step homogenization

a b s t r a c t For the accurate prediction of the Effective Thermal Conductivities (ETCs) of the Multi-Oriented Inclusions Reinforced Composites (MOIRCs) with imperfect interfaces, i.e., Weakly Conducting (WC) or Highly Conducting (HC) interfaces, this work develops a two-step mean-field homogenization method, in which the MOIRCs is virtually decomposed into a set of pseudo-grains, each of which is equivalent as a two-phase composite with imperfect interfaces and firstly homogenized by using the Mori-Tanaka (M-T) or DoubleInclusion (D-I) model, followed by the homogenization of all the pseudo-grains based on the assumption of the uniform intensity or heat flux and on the inclusion orientation distribution. The M-T and D-I models are derived by extending the solution of the Eshelby’s single inclusion problem to the heat transfer behaviors of the two-phase composites with imperfect interfaces. Through the comparison with the Finite Element (FE) homogenization method, the developed two-step mean-field homogenization method and the corresponding models are validated to accurately and efficiently predict the ETCs of the MOIRCs with imperfect interfaces. The simulation results show that the presence of WC interface decreases the ETCs of the MOIRCs, while the existence of HC interface has an opposite effect on the ETCs of the MOIRCs. In addition, the ETCs of the MOIRCs with imperfect interfaces show an asymptotic behavior with the variation of interfacial thermal conductance. © 2019 Elsevier Ltd. All rights reserved.

1. Introduction Since the 1950s, composites have been well developed and more widely used in the industries such that the prediction of the effective mechanical and physical properties of composites is critically required, as it is beneficial for both the application and design of composites. The heat transfer analysis of composites is non-trivial, and the modeling of the Effective Thermal Conductivities (ETCs) of composites thus has been widely studied. It is found that most of the modeling approaches assumes perfect interfaces between the matrix and inclusions of composites. However, these interfaces are generally imperfect [1,2] and have an immediate influence on the ETCs of composites. Concerning the prediction of the ETCs of the composites with imperfect interfaces, which specify the Weakly Conducting (WC) or Highly Conducting (HC) interfaces in this paper (unless otherwise indicated), many theoretical and numerical investigations

∗

Corresponding author. E-mail address: [email protected] (M.W. Fu).

https://doi.org/10.1016/j.ijheatmasstransfer.2019.119167 0017-9310/© 2019 Elsevier Ltd. All rights reserved.

have been conducted. The existing approaches can mainly be categorized into three classes: (1). The Laplace’s heat transfer equation based method. The solution of the Laplace’s heat transfer equation was initially solved by Rayleigh [3] and Maxwell [4], and was then extended to determine the ETCs of the oriented inclusions reinforced composites with perfect interfaces [5] and imperfect interfaces [6,7], and of the multi-oriented inclusions reinforced composites with perfect interfaces [8] and imperfect interfaces [9]. This method provides an analytical evaluation on the ETCs of the composites reinforced by the spherical, ellipsoidal and cylindrical inclusions with isotropic thermal conductivity. (2). The Eshelby’s single inclusion problem based method. This class of method originates from the solution of the famous Eshelby’s single inclusion problem initially developed for elasticity of heterogeneous material [10,11], which was then reformulated for the ETCs of the composites with perfect interfaces by Hatta and Taya [12,13]. The reformulated method was extended for the ETCs of the composites with imperfect interfaces by Quang et al. [14,15] (the case of spherical inclusion) and by Bon-

2

W. Tian, M.W. Fu and L. Qi et al. / International Journal of Heat and Mass Transfer 148 (2019) 119167

foh et al. [16,17] (the case of ellipsoidal inclusion). This approach can consider the inclusions with different shapes and anisotropic thermal conductivity, but is only available for the oriented inclusions reinforced two-phase composites. (3). The numerical modeling method. This class of methods employs numerical modeling methods, such as the Finite Element (FE) method [18,19], the extended finite element method [20] and the fourier transform method [21], to predict the ETCs of the composites with imperfect interfaces. Therefore, the complex micro-structures of composites can be fully considered. Meanwhile, the detailed field localization in each phase of composites can be presented. However, it might become extremely computational expensive, especially for multi-scale analysis. Considering the balance between capability, accuracy and efficiency of the above-mentioned methods, the mean-field homogenization method (one of the Eshelby’s single inclusion problem based methods) is adopted to predict the ETCs for the MultiOriented Inclusions Reinforced Composites (MOIRCs) with imperfect interfaces in this paper. Though the mean-field homogenization method has been well documented to predict the ETCs of the two-phase composites with imperfect interfaces, the corresponding works on the prediction of the ETCs for the MOIRCs with imperfect interfaces are seldom reported. Regarding the prediction of the effective mechanical properties of the MOIRCs, one available approach is the two-step meanfield homogenization method, initially proposed for the effective elastic properties of the MOIRCs [22,23] and then modified for the effective thermo-elastic, elasto-plastic and visco-elasto-plastic properties of the MOIRCs [24–26]. Analogously, the idea of the two-step mean-field homogenization method is adopted to predict the ETCs of the MOIRCs with imperfect interfaces in this research. Based on the pioneering works of Quang et al. [14,15] and Bonfoh et al. [16,17], the Mori-Tanaka (M-T) and Double-Inclusion (D-I) models for determining the ETCs of the two-phase composites with imperfect interfaces are firstly derived. The two-step mean-field homogenization method is then developed and the corresponding models are established to predict the ETCs of the MOIRCs with imperfect interfaces. The predicted results are compared with the “exact” predictions provided by the FE homogenization method to validate the established two-step mean-field homogenization models. Meanwhile, the effect of imperfect interface on the ETCs of the MOIRCs is discussed. This work would provide an accurate and efficient method for predicting the ETCs of the MOIRCs with imperfect interfaces and a comprehensive understanding of the effect of imperfect interface on the ETCs of composites.

inclusion I are linear and characterized by their thermal conductivities k0 and kI , respectively. The common interface S between the matrix 0 and the ellipsoidal inclusion I is imperfect. Denote q(x), e(x) (i.e., −∇ T (x )) and T(x) the local heat flux, intensity and temperature fields of the material point x in the RVE s , respectively. Here two types of imperfect interfaces, i.e., WC and HC interfaces, are considered. Regarding the WC interface, the Kapitza’s interfacial thermal resistance model [27], which assumes that the heat flux and temperature are continuous and discontinuous, respectively, at the interface S, is considered:

[q ( x ) ] · n ( x ) = 0

The Eshelby’s single inclusion problem initially proposed for elasticity of heterogeneous materials [10,11] is extended to determine the ETCs of the two-phase composites with imperfect interfaces. The topology (Representative Volume Element (RVE) s with the volume V) of the Eshelby’s single inclusion problem consists of an infinite homogeneous reference matrix 0 with the volume V0 and an ellipsoidal inclusion I with the volume VI . The heat transfer behaviors of the reference matrix 0 and the ellipsoidal

with

x ∈ S. (1)

where is the jump of the variable across the interface S. n(x) is the outward unit vector normal to the interface S. α is the Kapitza’s interfacial thermal resistance. For the HC interface, the interfacial thermal conductance model proposed by Miloh and Benveniste [28] is considered. In contrast to the Kapitza model, this HC interface model assumes that the heat flux and temperature are discontinuous and continuous, respectively, at the interface S: •

[q ( x ) ] · n ( x ) = β S T I

and

[T ( x ) ] = 0

with

x ∈ S.

(2)

where β represents the interfacial thermal conductance for the HC interface S, and the surface Laplacian S TI is defined as [28]:

S T I = −∇ · eI (x ) + n(x ) · (∇ · eI (x )) · n(x ) + n(x ) · eI (x ) ·(∇ · n(x ) ).

(3)

Given the boundary ∂ 0 of the infinite matrix 0 is remotely subjected to a homogeneous intensity E0 , the average intensity EI of the inclusion I can be related to the homogeneous intensity E0 by the following equation [16,17]:



E 0 = I + S I · r 0 · k I − H ( k S ) · α · ( S I − I ) · C · k I



+ ( 1 − H ( kS ) ) · β · S I · r 0 · D · E I ,

)−1 .

(4)

where kI = kI − k0 and r0 = (k0 S is the interior-point Eshelby’s thermal conductance tensor. C and D are the defined thermal conductance tensors for the WC and HC interfaces, respectively. The tensors S I , C and D are inclusion-size-dependent and their explicit expressions are given in Appendix A. I is the secondorder Kronecker tensor. H(kS ) is a Heaviside step function and defined as:



H ( kS ) =

1 for 0 for

I

kS  ≤ min (k0 , kI  ) kS  > max (k0 , kI  )

(5)

where kS denotes the ETC of the interface S and can be related to the interfacial thermal resistance α and the interfacial thermal conductance β as follows:

kS =

2.1. Eshelby’s single inclusion problem

[T (x )] = −α · q(x ) · n(x )

[•]

2. Micro-mechanical model The M-T and D-I models are derived to numerically estimate the ETCs of the two-phase composites with imperfect interfaces by extending the solution of the Eshelby’s single inclusion problem initially proposed for elasticity of heterogeneous materials in this section.

and

hS

α

·I

and kS =

β hS

· I,

(6)

where hS is the equivalent thickness of the interface S. Similar to the strain concentration tensor [29], an intensity concentration tensor Be is defined and introduced:

E I = Be · E 0 , where



Be

(7)

has the following formulation based on Eq. (4):

B e = I + S I · r 0 · k I − H ( k S ) · α · ( S I − I ) · C · k I + ( 1 − H ( kS ) ) · β · S I · r 0 · D

−1

.

Be

(8)

This intensity concentration tensor depends on the inclusion size, the thermal conductivity k0 of the reference matrix 0 , the thermal conductivity kI of the inclusion 1 and the interfacial thermal properties. It is noted that the intensity concentration tensor Be converges to the perfect interfacial bonding case when α → 0 or β → 0.

W. Tian, M.W. Fu and L. Qi et al. / International Journal of Heat and Mass Transfer 148 (2019) 119167

2.2. Mean-field homogenization model Based on the derived intensity concentration tensor Be (Eq. (8)), the mean-field homogenization models are developed to estimate the ETCs of the two-phase composites with imperfect interfaces in this section. Considering a RVE  of the two-phase composites, the RVE  are constituted by the aligned ellipsoidal inclusions I with the thermal conductivity kI , the finite matrix M with the thermal conductivity kM and the imperfect interfaces S between the matrix M and inclusions I , respectively. In the RVE , the volume fraction of the inclusions I is vI . The average intensity E and heat flux Q of the RVE  are computed from the local intensity e(x) and heat flux q(x) by the following equations:    1 1 1 T (x ) · n (x ) d x = e ( x ) d x − H ( kS ) · [T (x )] · n(x ) dx V ∂ V  V S = v I · E I + ( 1 − v I ) · E M + v I · H ( k S ) · α · C · kI · E I , (9)

E = −

 1 q (x ) · n (x ) · x d x V ∂   1 1 = q ( x ) d x + ( 1 − H ( kS ) ) · x · [ q ( x )] · n ( x ) d x V  V S



(17)

By substituting Eq. (17) to Eq. (13), the ETC k¯ of the two-phase composites is finally calculated. The D-I model [32] is another widely used mean-field homogenization model, and hypothesizes that each ellipsoidal inclusion I of the RVE  in the two-phase composites is wrapped in a hollow inclusion 0 with the identical aspect ratio, symmetry axis and center, and that the inclusion 0 is supported by the outer reference material R . Meanwhile, the volume ratio the inclusions I and M is the same with that of the actual two-phase composites. On one hand, by choosing kR = kM , the M-T model is retrieved:



−1

= Belower .

(18)

On the other hand, choosing kR = kI leads to an inverse M-T model with the following intensity concentration tensor Be :

vI · Q I + (1 − vI ) · Q M + vI · (1 − H (kS )) · β · D · Q I · (kI )−1 . (10)

Then, the ETC k¯ of the two-phase composites is related to the average intensity E and heat flux Q of the RVE based on the Fourier’s law [30]:

(11)

The average intensity EI of the inclusion I and the average intensity E of the RVE  can be related by another intensity concentration tensor BeI , which is expressed as:

E I = BeI · E .



BeI = Be · (1 − vI ) · I + vI · (I + H (kS ) · α · C · kI ) · Be .

+ ( 1 − H ( kS ) ) · β · S I · r M · D

Q =

Q = k¯ · E .

where rM = (kM )−1 . The Eshelby’s tensor S I is calculated by using kM instead of k0 . Based on Eqs. (7)–(10), (12), (15) and (16), the intensity concentration tensor BeI is explicitly related to the intensity concentration tensor Be :

B e = I + S I · r M · k I − H ( k S ) · α · ( S I − I ) · C · k I

and

=

3

(12)

By substituting Eqs. (9)–(11) to Eq. (12), the ETC k¯ of the two-phase composites is computed:

k¯ = kM + vI · {kI −kM −H (kS ) · α ·C · kI · kM + (1 −H (kS )) · β · D } · BeI .

Be = I − S I,r · rI · kI − H (kS ) · α · (S I,r − I ) · C · kM + (1 − H (kS ) ) · β · S I,r · rI · D = Beupper ,

(19)

where the Eshelby’s tensor S I,r is calculated by using kI instead of k0 , and rI = k−1 I . To more accurately predict the ETCs of the two-phase composites, the following Lielens’ type interpolative homogenization model, which is considered as a proper interpolation model from the M-T and inverse M-T models, is introduced [23]:



Be = (1 − ξ (v1 )) · (Belower )−1 + ξ (v1 ) · (Beupper )−1 and

1 2

ξ ( vI ) = vI ( 1 + vI ).

−1 (20)

Similarly, based on Eqs. (7)–(10), (12), (15) and (18)–(20), the intensity concentration tensor BeI is computed and then used to estimate the ETC k¯ of the two-phase composites based on Eq. (13).

(13) BeI

Generally, the intensity concentration tensor can be related to the intensity concentration tensor Be and is implicitly expressed by the following equation:

BeI = BeI (Be ).

(14)

The mean-field homogenization model can then be used to determine the intensity concentration tensors Be and BeI , and further determine the ETC k¯ of the two-phase composites. Note that various mean-field homogenization models only differ by the expression of the intensity concentration tensor Be . The M-T model [31] is one of the commonly used mean-field homogenization models and can provide the accurate prediction of effective properties of two-phase composites. The M-T model assumes that the matrix M of the RVE  of the two-phase composites plays the role of the infinite homogeneous reference matrix 0 of the RVE s in the Eshelby’s single inclusion problem and reads as:

k0 = kM

and E 0 = E M .

p = {cos θ , sin θ cos ϕ , sin θ sin ϕ } .

(16)

For all the inclusions in the RVE ω, their orientation distribution is characterized by using the Orientation Distribution Probability Function (ODPF) (p) [24], which defines the probability of the inclusions whose orientations are between p and p + d p.



B e = I + S I · r M · k I − H ( k S ) · α · ( S I − I ) · C · k I

−1

,

The two-step mean-field homogenization method is developed to predict the ETCs of the MOIRCs containing imperfect interfaces based on the micro-mechanical models for the ETCs of the twophase composites with imperfect interfaces in Section. 2. The RVE ω (the volume U) of the MOIRCs composes of the matrix (the volume U0 , the volume fraction v0 and the thermal conductivity of k0 ), a finite number of identical (except for orientation) ellipsoidal inclusions (the volume U1 , the volume fraction v1 and the thermal conductivity of k1 ) and the imperfect interfaces S between the matrix and inclusions. In the RVE ω, an inclusion orientation is characterized by a unit orientation vector p, which is related to two independent Euler angles: θ ∈ [0, π ] between p and the axis ox and φ ∈ [0, 2π ] between the projection of p on the plane oyz and the axis ox in the global coordinate system oxyz shown in Fig. 1:

(15)

Consequently, the intensity concentration tensor Be becomes:

+ ( 1 − H ( kS ) ) · β · S I · r M · D

3. Two-step mean-field homogenization method

T

(21)

4

W. Tian, M.W. Fu and L. Qi et al. / International Journal of Heat and Mass Transfer 148 (2019) 119167

Fig. 1. Illustration of an inclusion defined by two Euler orientation angles θ ∈ [0, π ] and φ ∈ [0, 2π ] with respect to the global and local cartesian coordinate systems: (a) Global cartesian coordinate system oxyz, and (b) Local cartesian coordinate system OXYZ.

Based on the orientation, the inclusions in the RVE ω can be categorized into a series of sets with the volume fraction of v1, i such that: N 

v1,i = v1 ,

(22)

Fig. 2. Two-step mean-field homogenization scheme for evaluation of the ETCs of the MOIRCs with imperfect interfaces.

i=1

where N represents the number of the categorized inclusion sets. The inclusions in each set (i) have the same orientation. The matrix is accordingly divided into N sets with the volume fraction of v0, i such that: N 

v0,i = v0 and

i=1

v1,i = v1 . v0,i + v1,i

(23)

By combining the inclusion set-i and the corresponding matrix set-i, N pseudo-grains are generated. Each pseudo-grain ωi (the volume of Vωi ) can be view as a two-phase composite containing the matrix in concentration v0 (the same as that of the matrix in the RVE ω), the oriented ellipsoidal inclusions in concentration v1 (the same as that of all the inclusions in the RVE ω) and the imperfect interfaces between the matrix and inclusions, and all the pseudo-grains can assemble the RVE ω. Therefore, the RVE ω is decomposed to a set of the pseudo-grains. It is obviously reasonable to hypothesize the perfect interfaces between the pseudo-grains. Therefore, the average intensity E of the RVE ω is related to the average intensity E ωi , p of the pseudo-grain ωi as:

E = e(x ) ω =

N N 1 1 Uωi · e(x ) ωi = Uωi · E ωi , p . U U i=1

(24)

i=1

According to the volume conservation condition Uωi = U · ( p), Eq. (24) is rewritten as:

E=

N  i=1

( p) · E ωi , p = E ωi , p ,

Q=

i=1





( p) · Q ωi , p = Q ωi , p .

(26)

Eqs. (25) and (26) show that the homogenization of the local intensity e(x) and heat flux q(x) in the RVE ω can be conducted in the two sequential steps following the decomposition of the RVE ω, as shown in Fig. 2: •

Second step: Homogenization (orientation average) of all the pseudo-grains.

The M-T and D-I models are generally used for homogenizing each pseudo-grain in the first step, while the Voigt [33] and Reuss [34] models can be used for the homogenization of all the pseudograins in the second step. Therefore, four two-step mean-field homogenization models, viz., the M-T/Voigt, M-T/Reuss, D-I/Voigt and D-I/Reuss models, are established. Regarding the Voigt model, the uniform intensity within all the pseudo-grains is assumed:

E = E ω1 , p = E ω2 , p = · · · = E ωi , p = · · · = E ωN , p . Thus, the ETC k¯ of the MOIRCs is derived as:







(27)







−1 ¯ k¯ = Q · E −1 = Q ωi , p · E −1 ω i , p = Q ω i , p · E ω i , p = kω i , p

N  =

( p) · k¯ ωi , p ,

First step: Homogenization of each pseudo-grain individually,

(28)

i=1

where k¯ ωi , p is the ETC of the pseudo-grain ωi with the orientation vector p. For the Reuss model, the uniform heat flux within all the pseudo-grains is assumed:

Q = Q ω1 , p = Q ω2 , p = · · · = Q ωi , p = · · · = Q ωN , p .

(29)

Through the Fourier’s law [30], the following expression is obtained:



−1 −1 k¯ = Q −1 · E = Q −1 ωi , p · E ωi , p = Q ωi , p · E ωi , p



=



k¯ ωi , p

−1

.

(30)

(25)

where • denotes the orientation-weighted average of the variable • over all the pseudo-grains. Similarly, the average heat flux Q of the RVE ω is related to the average heat flux Q ωi , p of the pseudo-grain ωi by the following equation: N 

•

Consequently, the ETC k¯ of the MOIRCs is derived as:

k¯ =



k¯ ωi , p

−1 −1



=

N 



( p) · k¯ ωi , p

−1

−1

.

(31)

i=1

Theoretically, the RVE ω can contain infinite ellipsoidal inclusions, and then the ODPF (p) of all the inclusions in the RVE ω is generally continuous and differentiable. In this case, the two-step mean-field homogenization method for the ETCs of the MOIRCs is still maintained, and the ETC k¯ of the MOIRCs has the following integral formulation instead of the sum formulation (Eqs. (28) and (31)):

   ψ ( p) · k¯ ωi , p d p for k¯ = k¯ ωi , p =

ω

Voigt model,

(32a)

W. Tian, M.W. Fu and L. Qi et al. / International Journal of Heat and Mass Transfer 148 (2019) 119167

5

Fig. 3. The introduced thin interphases between the matrix and inclusions for modeling imperfect interfaces in the MOIRCs.

k¯ =



k¯ ωi , p

−1 −1

=

 ω

−1 −1 ψ ( p) · k¯ ωi , p d p for Ruess model, (32b)

where ψ ( p) = d ( p)/d p is the inclusion orientation distribution probability density function [24]. It is noted that for the cases of the pure 3-D and 2-D random distribution, the inclusion orientation distribution probability density functions are ψ (θ , ϕ ) = 1/4π and ψ (θ , ϕ ) = 1/2π , respectively [35].

Eqs. (34)–(36), the heat flux q(x) and temperature T(x) of any material point x within the RVE ω are solved by adopting the FE method. Based on the solved local heat flux field q(x), the average heat flux Q = (Q1 , Q2 , Q3 )T of the RVE ω is computed by using Eq. (33) and given as follows:





 neint  1 1 Qi = qi ( x ) d x = Ve qi ( yI ) · W ( yI ) V ω V e I=1 with

4. Comparison with FE homogenization method The accuracy of the developed mean-field homogenization method is assessed by comparing with a full-field simulation method: the FE homogenization method [36], which is supposed to provide the “exact” prediction for the ETCs of the MOIRCs with imperfect interfaces. The imperfect interfaces can be viewed as limiting cases of a thin interphase layer [20,37,38] such that a thin interphase (see Fig. 3) between the matrix and an inclusion is introduced for modeling an imperfect interface. It is noted that the interfaces between the interphases and the others (the matrix and inclusions) are perfect. The ETC kS of the interphase is computed by using the interfacial thermal resistance α of the WC interface S (kS  < k1  and kS  < k0 ) or the interfacial thermal conductance β of the HC interface S (kS  > k1  and kS  > k0 )) through Eq. (6). Here the thickness of the interphase is selected as 50.0 nm. Consider a RVE ω (the volume of U) of the MOIRCs, in which the matrix ω0 (the volume of U0 ), the inclusions ω1 (the volume of U1 ) and the interphases ω2 (the volume of U2 ) are included. The matrix, inclusions and interphases are perfectly bound. Due to the absence of micro-defects and voids, the following average functions of the local heat flux q(x) within the RVE ω can be defined:

Q =

  1 1 q(x ) dx and Qωi = q(x ) dx, U ω Ui ωi i = 0, 1 and 2.

with (33)

To predict the ETCs of the MOIRCs, the steady state heat transfer equation is considered here:

div(q(x ) ) = 0.

(34)

The heat flux q(x) at the material point x in the RVE ω is related to its temperature gradient ∇ T(x) through the Fourier’s law [30]:

q(x ) = −k(x ) · ∇ T (x ).

(35)

The periodic boundary condition [39,40] is imposed on the boundary ∂ω of the RVE ω:

T ( x + ) − T ( x − ) = T x+

x−

on

∂ω,

(36)

where and are the corresponding material point pairs on the opposite boundaries ∂ω+ and ∂ω− of the RVE ω. Based on

i = 1, 2 and 3,

(37)

where neint and nint are the numbers of the integration points in the element e (with the volume of Ve ) and in the whole RVE ω, and W(yI ) is the integration point weight at the integration point positioned at yI in the element e [40]. The ETC k¯ of the MOIRCs is then derived as:

k¯ ii = −Qi /∇ Ti

with

i = 1, 2 and 3.

(38)

Here no Einstein’s summation convention is indicated. The temperature gradient ∇ Ti is computed by the following equation:

∇ Ti = T /Li with i = 1, 2 and 3,

(39)

where Li is the side length of the RVE ω in the direction xi . Note that the directions x1 , x2 and x3 refer to the directions x, y and z, respectively. 5. Numerical examples and validation 5.1. 1-D Randomly distributed inclusions reinforced composites In this section, the developed M-T and D-I models are validated to accurately estimate the ETCs of the 1-D randomly distributed inclusions reinforced composites with imperfect interfaces. The first investigated composites are the spherical SiC inclusions (with the volume fraction of 0.30 and the radius of 5.0 μm) reinforced aluminum alloy matrix composites with imperfect interfaces. The thermal conductivities of the aluminum alloy matrix and the SiC inclusions are isotropic: k0 = 173.0W · m−1 · K−1 and k1 = 393.0W · m−1 · K−1 , respectively [41]. The previous studies have shown that the ETCs of the randomly distributed spherical inclusions reinforced composites with imperfect interfaces are macroscopically isotropic such that the isotropic ETCs are used for characterizing the heat transfer behaviors of the composites [18,36]. The ETCs of the composites with the inclusion volume fraction of 0.30 are predicted by using the M-T model, the D-I model, the FE homogenization method and the Pavanello and Stefano (P-S) model [42,43], and plotted against the interfacial thermal conductance α¯ = 1/α and β in Fig. 4. The ETCs of the composites obtained by the M-T model agree well with those predicted by the P-S model and the FE homogenization method. On the other hand, the ETCs of the composites obtained by the D-I model agree well

6

W. Tian, M.W. Fu and L. Qi et al. / International Journal of Heat and Mass Transfer 148 (2019) 119167

Fig. 4. The ETCs of the spherical SiC inclusions reinforced aluminum alloy matrix composites with the WC or HC interfaces predicted by the M-T model, the D-I model, the P-S model and FE homogenization method: (a) WC interface and (b) HC interface.

Fig. 5. The normalized ETCs of the oriented dilute ellipsoidal inclusions reinforced composites with the WC interfaces as a function of the normalized interfacial thermal conductance α¯ : (a) Normalized longitudinal ETC κ L and (b) Normalized transverse ETC κ T .

with those predicted by the P-S model and the FE homogenization method in the case of α¯ ≥ 107 W · m−2 · K−1 and β ≤ 10−3 W · K−1 . The second example composites are the oriented dilute ellipsoidal inclusions (i.e., the volume fraction of 0.01) reinforced composites with the WC interfaces. The semi-minor axis a2 (a3 ) and the aspect ratio Ra of the inclusions are 2.5 μm and 5.0, respectively. The isotropic thermal conductivities of the matrix and inclusion are related as k1 = 5.0k0 . The ETCs of the oriented ellipsoidal inclusions reinforced composites with imperfect interfaces are transversely isotropic [16,17], and the longitudinal and transverse ETCs are adopted to characterize the heat transfer behaviors of the composites. Fig. 5 show the deviations of the normalized longitudinal and transverse ETCs κL = (k¯ L − k0 )/(v1 k0 ) and κT = (k¯ T − k0 )/(v1 k0 ) of the composites as a function of the normalized   interfacial thermal conductance α¯ = log

a21 −a23

α k0

. The difference

between the ETCs of the composites predicted by the M-T and D-I models and those reported in the work of Benveniste and Miloh [44] is not of significance. The third studied composites are identical with the second example composites except for the HC interfaces instead of the WC

interfaces and k1 = 0.2k0 . Fig. 6 presents the normalized longitudinal and transverse ETCs κL = (k¯ L − k0 )/(v1 k0 ) and κT = (k¯ T − k0 )/(v1 k0 ) of the composites regarding the normalized interfacial thermal conductance β¯ = log

β k0

a21 −a23

. The estimations of the

M-T and D-I models for the ETCs of the composites are identical with those of Ref. [28]. In addition, it is discovered that the ETCs of the composites with the WC interfaces are less than or equal to those of the composites with the perfect interfaces (i.e., α¯ → +∞ or β → −∞), which are less than or equal those of the composites with the HC interfaces. Therefore, it is concluded that the presence of the WC and HC interfaces decreases and increases the ETCs of the composites, respectively. Moreover, the ETCs of the composites with imperfect interfaces show an asymptotic behaviour, that is, behaving as a sigmoid curve [17,18]. Therefore, three zones can be identified for the curves of the ETCs of the composites: a transition zone and two plateaus zones. In the two plateaus zones, the effects of the interface thermal conductance α¯ and β on the ETCs of the composites are not significant. In the transition zone, the ETCs of the composites in-

W. Tian, M.W. Fu and L. Qi et al. / International Journal of Heat and Mass Transfer 148 (2019) 119167

7

Fig. 6. The normalized ETCs of the oriented dilute ellipsoidal inclusions reinforced composites with the HC interfaces regarding the normalized interfacial thermal conductance β¯ : (a) Normalized longitudinal ETC κ L and (b) Normalized transverse ETC κ T .

crease with the increase of the interfacial thermal conductance α¯ and β . 5.2. 2-D Randomly distributed inclusions reinforced composites The 2-D randomly distributed ellipsoidal inclusions (in the plane xoy of Fig. 1(a)) reinforced composites with imperfect interfaces are studied here. In the composites, the semi-minor axis a2 (a3 ), the aspect ratio Ra and the volume fraction v1 of the inclusions are 2.5 μm, 3.0 and 0.20, respectively. The matrix with the isotropic thermal conductivity k0 = 72.0W · m−1 · K−1 is considered, while the ellipsoidal inclusions with the isotropic and transversely isotropic thermal conductivities are considered, respectively. The thermal conductivity of the isotropic inclusion is k1 = 32.0W · m−1 · K−1 , and the longitudinal and transverse thermal conductivities of the transversely isotropic inclusion are kL1 = 32.0W · m−1 · K−1 and kT1 = 5.0W · m−1 · K−1 , respectively. Considering the WC interface with the interfacial thermal conductance α¯ (i.e., 1/α ) of 1.0 × 107 W · m−2 · K−1 , the ETCs of the composites evaluated by the M-T/Voigt model are given (in W · m−1 · K−1 ):



k¯ Iso =

0.00 58.97 0.00

0.00 0.00 55.49

57.06 0.00 0.00

0.00 57.06 0.00

0.00 0.00 . 51.67

 k¯ Tran =



58.97 0.00 0.00

and



On the other hand, considering the HC interface with the interfacial thermal conductance β of 1.0 × 10−3 W · K−1 , the predicted ETCs of the composites by the M-T/Voigt model are presented (in W · m−1 · K−1 ):



k¯ Iso =

0.00 129.15 0.00

0.00 0.00 100.85

128.83 0.00 0.00

0.00 128.83 0.00

0.00 0.00 . 100.37

 k¯ Tran =



129.15 0.00 0.00

and



It reveals that the ETCs of the composites with imperfect interfaces are transversely isotropic with the transverse plane xoy and the

symmetric axis z. The same conclusion can be obtained for the MT/Reuss, D-I/Voigt and D-I/Reuss models as well. Thus, the in-plane and out-of-plane ETCs are used for the 2-D randomly distributed inclusions reinforced composites with imperfect interfaces. The in-plane and out-of-plane ETCs of the composites with the WC interfaces are estimated by using the two-step mean-field homogenization models (i.e., the M-T/Voigt, M-T/Reuss, D-I/Voigt and D-I/Reuss models) and the FE homogenization method, and plotted against the interfacial thermal conductance α¯ in Figs. 7 and 8. The relative deviations between the estimations of the MT/Voigt, M-T/Reuss, D-I/Voigt and D-I/Reuss models and the estimations of the FE homogenization method are not more than 6.92%, and it thus indicates that all the M-T/Voigt, M-T/Reuss, DI/Voigt and D-I/Reuss models provide the accurate predictions of the ETCs of the composites with the WC interfaces. Moreover, it can be found that the M-T/Voigt model predicts the more accurate ETCs of the composites with the WC interfaces than the other models. Meanwhile, Figs. 7 and 8 demonstrate that the inplane and out-of-plane ETCs of the composites with the WC interfaces are less than or equal to those of the composites with the perfect interfaces (i.e., α¯ ≥ 1 × 1010 W · m−2 · K−1 ), which concludes that the presence of the WC interface decreases the ETCs of the composites. Regarding the interfacial thermal conductance β , Figs. 9 and 10 illustrate the in-plane and out-of-plane ETCs of the composites with the HC interfaces predicted by using the M-T/Voigt, MT/Reuss, D-I/Voigt, D-I/Reuss models and the FE homogenization method. It is shown that the M-T/Voigt, M-T/Reuss, D-I/Voigt and D-I/Reuss models provide the same predictions as the FE homogenization method when β ≤ 1.0 × 10−3 W · K−1 , and that only the predictions of the D-I/Voigt model agree well with those of the FE homogenization method when β > 1.0 × 10−3 W · K−1 . Meanwhile, the in-plane and out-of-plane ETCs of the composites with the perfect interfaces (i.e., β ≤ 1 × 10−7 × W · K−1 ) are less than or equal those of the composites with the HC interfaces, which reveals that the presence of the HC interface increases the ETCs of the composites. In addition, it is discovered that the ETCs of the composites with imperfect interfaces show an asymptotic behaviour with the variation of the interfacial thermal conductance and that a transition zone and two plateaus zones can thus be identified for the

8

W. Tian, M.W. Fu and L. Qi et al. / International Journal of Heat and Mass Transfer 148 (2019) 119167

Fig. 7. The ETCs of the 2-D isotropic inclusions reinforced composites with the WC interfaces predicted by using the M-T/Voigt, M-T/Reuss, D-I/Voigt, D-I/Reuss models and the FE homogenization method: (a) In-plane ETC k¯ 11 , and (b) Out-of-plane ETC k¯ 33 .

Fig. 8. The ETCs of the 2-D transversely isotropic inclusions reinforced composites with the WC interfaces predicted by using the M-T/Voigt, M-T/Reuss, D-I/Voigt, D-I/Reuss models and the FE homogenization method: (a) In-plane ETC k¯ 11 , and (b) Out-of-plane ETC k¯ 33 .

Fig. 9. The ETCs of the 2-D composite reinforced by the isotropic inclusions with the HC interfaces predicted by using the M-T/Voigt model, M-T/Reuss model, D-I/Voigt model, D-I/Reuss model and the FE homogenization method: (a) In-plane ETC k¯ 11 , and (b) Out-of-plane ETC k¯ 33 .

W. Tian, M.W. Fu and L. Qi et al. / International Journal of Heat and Mass Transfer 148 (2019) 119167

9

Fig. 10. The ETCs of the 2-D composite reinforced by the transversely isotropic inclusions with the HC interfaces predicted by using the M-T/Voigt model, M-T/Reuss model, D-I/Voigt model, D-I/Reuss model and the FE homogenization method: (a) In-plane ETC k¯ 11 , and (b) Out-of-plane ETC k¯ 33 .

Fig. 11. The isotropic ETCs of 3-D randomly distributed fibers reinforced composites with the WC interfaces predicted by using the M-T/Voigt, M-T/Reuss, D-I/Voigt, D-I/Reuss models and the FE homogenization method: (a) Isotropic fibers, and (b) Transversely isotropic fibers.

curves of the ETCs of the composites with imperfect interfaces, as mentioned in Section. 5.1.

:

 k¯ Iso =

5.3. 3-D Randomly distributed inclusions reinforced composites The composites reinforced by the 10.0 vol.% 3-D randomly distributed cylindrical fibers with the radius r of 2.5 μm and the aspect ratio Ra of 10.0 are investigated in this section. The imperfect interfaces between the matrix and fibers are assumed. The thermal properties of the matrix are isotropic, while the thermal properties of the fibers are isotropic and transversely isotropic, respectively. The thermal conductivities of the matrix and the fibers are identical with those of the matrix and ellipsoidal inclusions in Section. 5.2. In the case of the WC interface (the interfacial thermal conductance α¯ (1/α ) of 1.0 × 106 W · m−2 · K−1 ), the ETCs of the composites predicted by the M-T/Voigt model are given (in W · m−1 · K−1 )

0.00 62.67 0.00

0.00 0.00 62.67

62.55 0.00 0.00

0.00 62.56 0.00

0.00 0.00 . 62.57

 k¯ Tran =



62.65 0.00 0.00

and



In the case of the HC interface (the interfacial thermal conductance β of 1.0 × 10−2 W · K−1 ), the ETCs of the composites predicted by the M-T/Voigt model are given (in W · m−1 · K−1 ):

 k¯ Iso =

0.00 176.25 0.00

0.00 0.00 176.41

175.65 0.00 0.00

0.00 176.23 0.00

0.00 0.00 . 176.40

 k¯ Tran =



175.68 0.00 0.00

and



10

W. Tian, M.W. Fu and L. Qi et al. / International Journal of Heat and Mass Transfer 148 (2019) 119167

Fig. 12. The isotropic ETCs of 3-D randomly distributed fibers reinforced composites with the HC interfaces predicted by using the M-T/Voigt, M-T/Reuss, D-I/Voigt, D-I/Reuss models and the FE homogenization method: (a) Isotropic fibers, and (b) Transversely isotropic fibers.

It is found that the ETCs of the composites with imperfect interfaces are macroscopically isotropic, and the same conclusion can be obtained as well for the other two-step mean-field homogenization models (i.e. the M-T/Reuss, D-I/Voigt and D-I/Reuss models). Thus, the ETCs of the 3-D randomly distributed cylindrical fibers reinforced composites with imperfect interfaces are isotropic. The isotropic ETCs of the composites are predicted by using the M-T/Voigt, M-T/Reuss, D-I/Voigt, D-I/Reuss models and the FE homogenization method, respectively, as illustrated in Figs. 11 and 12. Regarding the composites with the WC interfaces, the maximum relative deviation between the predictions of the M-T/Voigt, M-T/Reuss, D-I/Voigt and D-I/Reuss models and those of the FE homogenization method is 3.6%, which thus validates the effectiveness of the M-T/Voigt, M-T/Reuss, D-I/Voigt and D-I/Reuss models to accurately predict the ETCs of the composites with the WC interfaces. Meanwhile, the M-T/Voigt model provides the most accurate predictions among all the two-step mean-field homogenization models. For the composites with the HC interfaces, all the M-T/Voigt, M-T/Reuss, D-I/Voigt and D-I/Reuss models can provide the reasonable predictions when β ≤ 1.0 × 10−3 W · K−1 . However, when β > 1.0 × 10−3 W · K−1 the M-T/Voigt and D-I/Voigt models provide the accurate predictions but the M-T/Reuss and D-I/Reuss models do not. In addition, the phenomena that the presence of the WC and HC interfaces decreases and increases, respectively, the ETCs of the composites and that the ETCs of the composites with imperfect interfaces show an asymptotic behaviour with an inflection point in the middle, can be found from Figs. 11 and 12, as aforementioned. In conclusion, the developed two-step mean-field homogenization method can accurately predict the ETCs of the MOIRCs with the WC interfaces or the HC interfaces when the interfacial thermal conductance is less than a certain value, and the M-T/Voigt model provides the best estimations. Therefore, the M-T/Voigt model will be adopted for predicting the ETCs for the MOIRCs with imperfect interfaces in future research.

faces. The MOIRCs with imperfect interfaces is decomposed into a set of pseudo-grains, and each pseudo-grain can be viewed as a two-phase composite with imperfect interfaces. The ETCs of the pseudo-grains are homogenized firstly by using the M-T or D-I model, which is derived based on the solution of the Eshelby’s single inclusion problem. The ETCs of the MOIRCs with imperfect interfaces are then predicted by homogenizing all the pseudo-grains based on the assumption of the uniform intensity or heat flux within all the pseudo-grains, respectively. The FE homogenization method is employed to provide the “exact” ETCs of the MOIRCs with imperfect interfaces. Comparing with the FE homogenization method, the two-step mean-field homogenization models, i.e. the M-T/Voigt, M-T/Reuss, D-I/Voigt and DI/Reuss models, are validated to be able to accurately predict the ETCs of the MOIRCs with the WC interface or with the HC interface when the interfacial thermal conductance is less than a certain value. Meanwhile, it is found that the M-T/Voigt model can provide the most accurate prediction of the ETCs of the MOIRCs with imperfect interfaces. The presence of the WC and HC interfaces decreases and increases the ETCs of the composites, respectively. Meanwhile, the ETCs of the composites with imperfect interfaces show an asymptotic behaviour with the variation of the interfacial thermal conductance, and three regions can be identified for the curves of the ETCs of the composites with imperfect interfaces, viz., a transition region and two plateaus regions. In the two plateaus regions, the ETCs of the composites almost maintain constant with the variation of the interfacial thermal conductance, while the ETCs of the composites are increased with the interfacial thermal conductance in the transition region.

Declaration of Competing Interest None.

6. Conclusions Acknowledgement In this research, the two-step mean-field homogenization method is developed and the corresponding models are established to accurately predict the ETCs of the MOIRCs with imperfect inter-

The authors thank the financial support from the project No.1YW3H from The Hong Kong Polytechnic University.

W. Tian, M.W. Fu and L. Qi et al. / International Journal of Heat and Mass Transfer 148 (2019) 119167

Appendix A. Determination of the tensors S I , C and D The explicit expressions of the Eshelby’s tensor S I , the tensor C and the tensor D for some specific cases are given in the following.

a. Generalized prolate spheroidal inclusion: λ1 > λ2 = λ3 (λi =  ai / k0,ii , ai is the semi-axis of the ellipsoidal inclusion in the direction xi ), a1 > a2 = a3 and k0,22 = k0,33 : I S11 = I S22 =

λ1 λ

√ 2

λ1

λ

λ

(λ1 /λ3 )−λ

(λ21 −λ23 )

3/2

λ λ

2− 2− 1 1 3

(

√

λ

2 3 arccosh 2 − 2 3/2 1 3

2 λ

I I S33 = S22 .

2 3

)

λ

λ

2− 2 1 3

,

(λ1 /λ3 )

,

(A.1)

=

I S22 =

λ23

√

λ23 −λ21 −λ1 λ23 arccos(λ1 /λ3 )

(λ23 −λ21 ) √ λ1 λ23 arccos(λ1 /λ3 )−λ21 λ23 −λ21 3/2

(

2 λ23 −λ21

I I S33 = S22 .

)

3/2

, (A.2)

,

c. Generalized isotropic infinite circular cylinder: λ1 → ∞, λ2 = λ3 , a1 → ∞, a2 = a3 and k0,22 = k0,33 : I S11 = 0,

1 2

I S22 =

I S33 =

and

1 . 2

(A.3)

d. Generalized isotropic sphere: λ1 = λ2 = λ3 , a1 = a2 = a3 = a and k0,11 = k0,22 = k0,33 : I S11 =

1 , 3

I S22 =

1 3

I S33 =

and

1 . 3

(A.4)

A2. Components of the tensor C a. Prolate spheroidal inclusion: a1 > a2 = a3 and k0,22 = k0,33 :

C11 = C22 =

3a3 a21 arccos(a3 /a1 )−3a23 3a1 a3

(

√

C33 = C22 .

√

a21 −a23

,

) (a21 −2a23 )arccos(a3 /a1 ) , 3/2 4a3 (a21 −a23 )

2a1 a21 −a23 a21 −a23 +3a1

3/2

(A.5)

b. Oblate spheroidal inclusion: a1 < a2 = a3 and k0,22 = k0,33 :

C11 =

3a23

√

a23 −a21 −3a3 a21 arccosh(a3 /a1 )

,

( ) √ 3a1 (2a23 −a21 )arccosh(a3 /a1 )−3a1 a3 a23 −a21 C22 = , 3/2 2 2 4a3 (a3 −a1 ) C33 = C22 .

2a1 a23 −a21

3/2

C11 = 0,

C22 =

C11 =

1 , a

C22 =

3π 8a3

and

C33 =

3π . 8a3

(A.7)

(A.6)

1 a

and

C33 =

1 . a

(A.8)

A3. Components of the tensor D The components of the tensor D are calculated by using the components of the tensor C :

D11 = C22 + C33 ,

b. Generalized oblate spheroidal inclusion: λ1 < λ2 = λ3 , a1 < a2 = a3 and k0,22 = k0,33 : I S11

c. Isotropic infinite circular cylinder: a1 → ∞, a2 = a3 and k0,22 = k0,33 :

d. Isotropic sphere: a1 = a2 = a3 = a and k0,11 = k0,22 = k0,33 :

A1. Components of the Eshelby’s tensor S I

2 3 arccosh

11

D22 = C11 + C33

and

D33 = C11 + C22 .

(A.9)

Appendix B. Details of the FE homogenization models The periodic RVEs containing the matrix, inclusions and interphases used for the FE homogenization of the ETCs of the composites with imperfect interfaces are generated by using the modified RSA algorithm [35,36], and the detailed generation procedures of the RVEs can be found in Ref. [36]. The generated RVEs of the composites reinforced by the 30.0 vol.% 3-D randomly distributed spherical inclusions, the 20.0 vol.% 2-D randomly distributed ellipsoidal inclusions and the 10.0 vol.% 3-D randomly distributed cylindrical fibers are, respectively, illustrated in Figs. B.1, B.2 and B.3. The sizes of the inclusions/fibers, interphases and the generated RVEs are given as follows: (1). The RVE with the 3-D randomly distributed spherical inclusions: the inclusion radius of 5.0 μm, the interphase thickness of 50.0nm and the RVE size of 40.2 μm × 40.2 μm × 40.2 μm; (2). The RVE with the 2-D randomly distributed ellipsoidal inclusions: the inclusion major length of 15.0 μm, the inclusion minor length of 5.0 μm, the interphase thickness of 50.0nm and the RVE size of 60.2 μm × 60.2 μm × 20.2 μm; (3). The RVE with the 3-D randomly distributed cylindrical fibers: the fiber length of 50.0 μm, the fiber radius of 2.5 μm, the interphase thickness of 50.0nm and the RVE size of 100.2 μm × 100.2 μm × 100.2 μm. The positions and orientation of the inclusions in the generated RVEs are provided in Appendix C: Supplementary material. The periodic boundary condition is introduced to homogenize the ETCs of the composites, and the numerical implementation algorithm of the periodic boundary condition is given in the previous work of the authors [40]. The mesh type of the FE models is selected as 4-node linear heat transfer tetrahedron finite element,

Fig. B1. The RVE of the composites reinforced by the 30.0 vol.% 3-D randomly distributed spherical inclusions: (a) Interphase, and (b) Inclusion and (c) Assembly of interphases and inclusions.

12

W. Tian, M.W. Fu and L. Qi et al. / International Journal of Heat and Mass Transfer 148 (2019) 119167

Fig. B2. The RVE of the composites reinforced by the 20.0 vol.% 2-D randomly distributed ellipsoidal inclusions: (a) Interphase, and (b) Inclusion and (c) Assembly of interphases and inclusions.

Fig. B3. The RVE of the composites reinforced by the 10.0 vol.% 3-D randomly distributed cylindrical fibers: (a) Interphase, and (b) Fiber and (c) Assembly of interphases and fibers.

Table B1 Mesh sizes of the inclusions/fiber, interphases and matrices of the RVEs of the composites. RVE #1, RVE #2 and RVE #3 refer to the RVE with the spherical, ellipsoidal and cylindrical inclusions/fibers, respectively. Mesh size (μm)

Matrix Inclusion/Fiber Interphase

RVE #1

RVE #2

RVE #3

1.00 0.75 0.50

1.00 0.75 0.50

2.00 1.00 1.50

and the mesh sizes of the inclusions/fibers, interphases and matrices of the RVEs of the composites are given in Table. B.1, which are sufficient to obtain convergent results. Supplementary material The positions (center points) and orientation (two Euler orientation angles) of the inclusions/fibers in the RVEs of the composites reinforced by the 30.0 vol.% 3-D randomly distributed spherical inclusions, the 20.0 vol.% 2-D randomly distributed ellipsoidal inclusions and the 10.0 vol.% 3-D randomly distributed cylindrical fibers, respectively, are given in this section. Supplementary material associated with this article can be found, in the online version, at 10.1016/j.ijheatmasstransfer.2019. 119167. References [1] X. Ren, Y. Dai, J. Gou, W. Tao, Numerical prediction of thermal contact resistance of 3D C/C-SiC needled composites based on measured practical topography, Int. J. Heat Mass Transf. 131 (2019) 176–188.

[2] Y. Koutsawa, A. Karatrantos, W. Yu, D. Ruch, A micromechanics approach for the effective thermal conductivity of composite materials with general linear imperfect interfaces, Compos. Struct. 200 (2018) 747–756. [3] L. Rayleigh, LVI. On the influence of obstacles arranged in rectangular order upon the properties of a medium, London Edinburgh Dublin Philos.Mag. J. Sci. 34 (1892) 481–502. [4] J.C. Maxwell, A Treatise on Electricity and Magnetism, vol. 1, Clarendon press, 1881. [5] A.S. Sangani, C. Yao, Bulk thermal conductivity of composites with spherical inclusions, J. Appl. Phys. 63 (1988) 1334–1341. [6] M.P. Lutz, R.W. Zimmerman, Effect of the interphase zone on the conductivity or diffusivity of a particulate composite using Maxwell’s homogenization method, Int. J. Eng. Sci. 98 (2016) 51–59. [7] D. Hasselman, L.F. Johnson, Effective thermal conductivity of composites with interfacial thermal barrier resistance, J. Compos. Mater. 21 (1987) 508–515. [8] H. Duan, B. Karihaloo, J. Wang, X. Yi, Effective conductivities of heterogeneous media containing multiple inclusions with various spatial distributions, Phys. Rev. B 73 (2006) 174203. [9] H. Duan, B. Karihaloo, Effective thermal conductivities of heterogeneous media containing multiple imperfectly bonded inclusions, Phys. Rev. B 75 (2007) 64206. [10] J.D. Eshelby, R.E. Peierls, The determination of the elastic field of an ellipsoidal inclusion, and related problems, Proc. R. Soc. London Ser.A 241 (1957) 376–396. [11] J.D. Eshelby, The elastic field outside an ellipsoidal inclusion, Proc. R. Soc. London Ser.A 252 (1959) 561–569. [12] H. Hatta, M. Taya, Effective thermal conductivity of a misoriented short fiber composite, J. Appl. Phys. 58 (1985) 2478–2486. [13] H. Hatta, M. Taya, Equivalent inclusion method for steady state heat conduction in composites, Int. J. Eng. Sci. 24 (1986) 1159–1172. [14] H.L. Quang, D.C. Pham, G. Bonnet, Q.-C. He, Estimations of the effective conductivity of anisotropic multiphase composites with imperfect interfaces, Int. J. Heat Mass Transf. 58 (2013) 175–187. [15] H.L. Quang, Q.-C. He, G. Bonnet, Eshelby’S tensor fields and effective conductivity of composites made of anisotropic phases with Kapitza’s interface thermal resistance, Philos. Mag. 91 (2011) 3358–3392. [16] N. Bonfoh, C. Dreistadt, H. Sabar, Micromechanical modeling of the anisotropic thermal conductivity of ellipsoidal inclusion-reinforced composite materials with weakly conducting interfaces, Int. J. Heat Mass Transf. 108 (2017) 1727–1739. [17] N. Bonfoh, H. Sabar, Anisotropic thermal conductivity of composites with ellipsoidal inclusions and highly conducting interfaces, Int. J. Heat Mass Transf. 118 (2018) 498–509.

W. Tian, M.W. Fu and L. Qi et al. / International Journal of Heat and Mass Transfer 148 (2019) 119167 [18] D. Marcos-Gómez, J. Ching-Lloyd, M. Elizalde, W. Clegg, J. Molina-Aldareguia, Predicting the thermal conductivity of composite materials with imperfect interfaces, Compos. Sci. Technol. 70 (2010) 2276–2283. [19] W. Yang, K. Peng, L. Zhou, J. Zhu, D. Li, Finite element simulation and experimental investigation on thermal conductivity of diamond/aluminium composites with imperfect interface, Comput. Mater. Sci 83 (2014) 375–380. [20] J. Yvonnet, Q.-C. He, C. Toulemonde, Numerical modelling of the effective conductivities of composites with arbitrarily shaped inclusions and highly conducting interface, Compos. Sci. Technol. 68 (2008) 2818–2825. [21] H.L. Quang, T.-L. Phan, G. Bonnet, Effective thermal conductivity of periodic composites with highly conducting imperfect interfaces, Int. J. Therm. Sci. 50 (2011) 1428–1444. [22] C.W. Camacho, C.L. Tucker III, S. Yalvaç, R.L. McGee, Stiffness and thermal expansion predictions for hybrid short fiber composites, Polym. Compos. 11 (1990) 229–239. [23] G. Lielens, Micro-Macro Modeling of Structured Materials, Universite Catholique de Louvain„ 1999 Ph.D. thesis. [24] O. Pierard, C. Friebel, I. Doghri, Mean-field homogenization of multi-phase thermo-elastic composites: a general framework and its validation, Compos. Sci. Technol. 64 (2004) 1587–1603. [25] I. Doghri, L. Tinel, Micromechanics of inelastic composites with misaligned inclusions: numerical treatment of orientation, Comput. Methods Appl. Mech. Eng. 195 (2006) 1387–1406. [26] C. Friebel, I. Doghri, V. Legat, General mean-field homogenization schemes for viscoelastic composites containing multiple phases of coated inclusions, Int. J. Solids Struct. 43 (2006) 2513–2541. [27] P. Kapitza, The study of heat transfer in helium ii, J. Phys. (Moscow) (1941) 181–210. [28] T. Miloh, Y. Benveniste, On the effective conductivity of composites with ellipsoidal inhomogeneities and highly conducting interfaces, Proc. R. Soc. London Ser.A 455 (1999) 2687–2706. [29] R. Hill, Elastic properties of reinforced solids: some theoretical principles, J. Mech. Phys. Solids 11 (1963) 357–372. [30] N. Bonfoh, A. Jeancolas, F. Dinzart, H. Sabar, M. Mihaluta, Effective thermal conductivity of composite ellipsoid assemblages with weakly conducting interfaces, Compos. Struct. 202 (2018) 603–614.

13

[31] T. Mori, K. Tanaka, Average stress in matrix and average elastic energy of materials with misfitting inclusions, Acta Metall. 21 (1973) 571–574. [32] M. Hori, S. Nemat-Nasser, Double-inclusion model and overall moduli of multi-phase composites, Mech. Mater. 14 (1993) 189–206. [33] W. Voigt, Ueber die beziehung zwischen den beiden elasticit–tsconstanten isotroper körper, Ann. Phys. 274 (1889) 573–587. [34] A. Reuss, Berechnung der fließgrenze von mischkristallen auf grund der plastizit–tsbedingung für einkristalle, ZAMM 9 (1929) 49–58. [35] W. Tian, L. Qi, J. Zhou, J. Liang, Y. Ma, Representative volume element for composites reinforced by spatially randomly distributed discontinuous fibers and its applications, Compos. Struct. 131 (2015) 366–373. [36] W. Tian, M. Fu, L. Qi, X. Chao, J. Liang, Interphase model for fe prediction of the effective thermal conductivity of the composites with imperfect interfaces, Int. J. Heat Mass Transf. 145 (2019) 118796. [37] Z. Hashin, Thin interphase/imperfect interface in conduction, J. Appl. Phys. 89 (2001) 2261–2267. [38] Y. Benveniste, Models of thin interphases and the effective medium approximation in composite media with curvilinearly anisotropic coated inclusions, Int. J. Eng. Sci. 72 (2013) 140–154. [39] W. Tian, L. Qi, C. Su, J. Liu, J. Zhou, Effect of fiber transverse isotropy on effective thermal conductivity of metal matrix composites reinforced by randomly distributed fibers, Compos. Struct. 152 (2016) 637–644. [40] W. Tian, L. Qi, X. Chao, J. Liang, M. Fu, Numerical evaluation on the effective thermal conductivity of the composites with discontinuous inclusions: periodic boundary condition and its numerical algorithm, Int. J. Heat Mass Transf. 134 (2019) 735–751. [41] A. Geiger, D. Hasselman, K. Donaldson, Effect of reinforcement particle size on the thermal conductivity of a particulate silicon carbide-reinforced aluminium-matrix composite, J. Mater. Sci. Lett. 12 (1993) 420–423. [42] F. Pavanello, F. Manca, P. Luca Palla, S. Giordano, Generalized interface models for transport phenomena: unusual scale effects in composite nanomaterials, J. Appl. Phys. 112 (2012) 84306. [43] F. Pavanello, S. Giordano, How imperfect interfaces affect the nonlinear transport properties in composite nanomaterials, J. Appl. Phys. 113 (2013) 154310. [44] Y. Benveniste, T. Miloh, The effective conductivity of composites with imperfect thermal contact at constituent interfaces, Int. J. Eng. Sci. 24 (1986) 1537–1552.