Numerical evaluation on the effective thermal conductivity of the composites with discontinuous inclusions: Periodic boundary condition and its numerical algorithm

International Journal of Heat and Mass Transfer 134 (2019) 735–751

Contents lists available at ScienceDirect

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

Numerical evaluation on the effective thermal conductivity of the composites with discontinuous inclusions: Periodic boundary condition and its numerical algorithm Wenlong Tian a,b, Lehua Qi a,⇑, Xujiang Chao a, Junhao Liang c, M.W. Fu b,** a b c

School of Mechanical Engineering, Northwestern Polytechnical University, Xi’an 710072, PR China Department of Mechanical Engineering, The Hong Kong Polytechnic University, Hung Hom, Kowloon, Hong Kong, China School of Materials Science and Engineering, Northwestern Polytechnical University, Xi’an 710072, PR China

a r t i c l e

i n f o

a b s t r a c t

Article history: Received 24 May 2018 Received in revised form 21 December 2018 Accepted 16 January 2019

Boundary condition plays an important role in prediction of the effective thermal conductivity of composites. In this research, the periodic boundary condition and the representative volume element (RVE) based finite element (FE) homogenization method are adopted to evaluate the effective thermal conductivities of the composites reinforced by the spherical, ellipsoidal and cylindrical inclusions, and the emphases are on the numerical implementation algorithm and validation of the periodic boundary condition. The heat flux continuity of the node pairs on the opposite surfaces of the RVEs of the composites is analyzed and the effective thermal conductivity of the composites are homogenized. The results show that the heat flux continuity of the node pairs on the opposite surfaces of the RVEs of the composites can be guaranteed by the proposed numerical implementation algorithm for the periodic boundary condition, and that the predicted effective thermal conductivities of the composites agree well with those determined by the Lewis-Nielsen model and the experimental tests. Therefore, the RVE based FE homogenization method with the periodic boundary condition can accurately evaluate the effective thermal conductivity of the composites with discontinuous inclusions. Ó 2019 Elsevier Ltd. All rights reserved.

Keywords: FE homogenization Heat flux continuity Numerical algorithm Periodic boundary condition Thermal conductivity

1. Introduction Composites have been widely used in the automotive, aviation, aerospace and national defense industries, etc., because of their excellent properties. Therefore, it claims an essential requirement of the accurate characterization for the properties of composites. Due to the advantages of acceptable accuracy, time-saving and cost efficiency, the theoretical and numerical methods are becoming more and more attracting to investigate the effective (or macroscopic) properties of composites, such as thermo-elastic properties [1,2], elasto-plastic behaviors [3,4] and thermal conductivities [5,6], as an alternative or at least a complement to experimental test. Over the last decades, the theoretical models have been well developed to evaluate the effective thermal conductivity of composites, ranging from the simple weight average models, such as

⇑ Corresponding author. ⇑⇑ Corresponding author. E-mail addresses: (M.W. Fu).

[email protected]

(L.

Qi),

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

[email protected]

series model [7] and parallel model [7], that only account for the volume fraction and thermal conductivities of the constituents to the complicated formulations, for example Lewis-Nielsen model [7,8], which consider the spatial and orientation variation of inclusions. These theoretical models are appropriate for the composites with the simple micro-structures or with regular inclusions (i.e. spherical [9,10], ellipsoidal [11,12] and cylindrical [13,14]), but not for the composites with very complicated micro-structures. Moreover, these theoretical models just provide the evaluations on the effective thermal conductivity of composites and cannot provide the detailed thermal flux and temperature fields, which are very crucial under certain conditions. Therefore, this research will adopt the numerical method to predict the effective (or macroscopic) thermal conductivity of composites. For studying the effective thermal conductivity of composites, a preferable numerical method is the RVE based FE homogenization method [15,16], which not only provides the accurate predictions for the effective thermal conductivity for the composites but presents the detailed thermal flux and temperature fields within the composites as well. Marcos-Gómez et al. [17] adopted the RVE based FE homogenization method to predict the effective thermal

736

W. Tian et al. / International Journal of Heat and Mass Transfer 134 (2019) 735–751

conductivities of the composite with spherical and cylindrical inclusions. The effects of inclusion anisotropy, inclusion orientation distribution, thermal interface conductance and inclusion dimensions on the effective thermal conductivity of the composites have been considered. Jiang et al. [18] numerically evaluated the effective thermal conductivity and temperature distribution of 3-D braided composites by using the modified RVE based FE homogenization model. The predicted results were compared against those obtained by the experimental test to demonstrate the accuracy of the modified RVE based FE homogenization model. In the RVE based FE homogenization method, it is unnecessary and not efficient to simulate the whole micro-structures of composites. Instead, a region of interest generally named RVE [19,20] is selected. The RVE has a certain boundary with the surrounding environment, and therefore the numerical simulations have to consider the physical processes in the boundary regions of the RVE by choosing the appropriate boundary conditions. Different boundary conditions may cause the quite different simulation results, and improper boundary conditions even may introduce the nonphysical influences on the simulation results. Thus, the boundary condition becomes critical when the RVE based FE homogenization method is introduced to evaluate the effective thermal conductivity of composites. The existing investigations mainly focus on the periodic boundary condition and demonstrate that the periodic boundary condition can accurately predict the effective thermal conductivity of composites. Gou et al. [21] studied the temperature distribution in a representative volume element of 3-D four-directional braided composites to predict the effective thermal conductivity of the composite by using the periodic boundary condition and the FE homogenization method. Dong et al. [22] adopted a multi-scale FE homogenization method with the periodic boundary condition to characterize the thermal conductive behaviors of 3-D braided composites along in-plane and out-of-plane directions. Although there exists a large quantity of work concerning the predictions of effective thermal conductivity of composites by using the periodic boundary condition, to the best of authors’ knowledge, there are less literatures on the validation of the numerical algorithm of the periodic boundary condition for the evaluation of the effective thermal conductivity of composites. In this paper, the RVE based FE homogenization method and the periodic boundary condition are adopted to evaluate the effective thermal conductivities of the composites with discontinuous inclusions, and the emphases are on the periodic boundary condition, its numerical algorithm and validation. The RVEs of the composites with the spherical, ellipsoidal and cylindrical inclusions are generated by using the modified RSA algorithm [23,24]. The heat flux continuity of the node pairs on the opposite surfaces of the generated RVEs for the composites is analyzed to validate the proposed numerical algorithm for imposing the periodic boundary condition. The predicted effective thermal conductivities of the composites are compared against those obtained by the Lewis-Nielsen model [7,8] and by the experimental tests to verify the RVE based FE homogenization method with the periodic boundary condition.

and hear flux) within the RVE x, the following average functions are defined,




Z
 1 lðxÞ dx and l xi V x Z 1 ¼ lðxÞ dxi ; with i ¼ 0; 1: V i xi

l ¼

ð1Þ

For any material point x within the RVE x, the following heatconduction equation coming from the first law of thermodynamics (also called the conservation of energy) should be satisfied,

qðxÞ  cðxÞ 

@T ðxÞ þ div ðqðxÞÞ  Q ðxÞ ¼ 0; @t

ð2Þ

where qðxÞ and cðxÞ are the density and the specific heat of the material point x. T ðxÞ and qðxÞ are the temperature and heat flux of the material point x. Q ðxÞ is the rate of internal energy generation per unit volume of the material point x. By using the Fourier’s law [25], the heat flux qðxÞ can be expressed as,

qðxÞ ¼ kðxÞ  rT ðxÞ;

ð3Þ

where k is the thermal conductivity of the material point x. Combining with Eq. (2) and considering the following boundary conditions,

T ðxÞ ¼ T ðxÞ on @ XT


ðxÞ on @ Xq ; and qðxÞ ¼ q

ð4Þ

the heat flux qðxÞ and temperature T ðxÞ of any material point x within the RVE x can be solved. The presence of internal heat generation and boundary conditions may strongly affect the temperature and heat flux fields in the RVE x. The effective thermal conductivity, however, is intrinsic to materials and should not depend on such external effects. Therefore, the steady state heat transfer with the prescribed temperature boundary condition is considered to evaluate the effective thermal conductivity of the RVE x,

div ðk  rT Þ ¼ 0 with T ¼ T

on @ XT :

ð5Þ

By introducing the FE method, the heat flux qðxÞ and temperature T ðxÞ of any material point x within the RVE x are solved, and the effective heat flux vector hqi ¼ ðhqi1 ; hqi2 ; hqi3 ÞT can be obtained by the post-processing of the FE analysis results and is given as follows, neint X 1 X Ve qi ðyI Þ  W ðyI Þ hqii ¼ V RVE e I¼1

! with i ¼ 1; 2 and 3;

ð6Þ

where neint and nint are the numbers of the integration point in the element e and in the whole RVE, respectively. W ðyI Þ is the integration point weight at an integration point positioned at yI in the element e, whose volume is V e . Then, the effective thermal conductivity tensor hkiij is derived based on the effective heat flux vector hqii and the temperature gradient vector rT j ,

hkiij ¼ hqii =rT j :

ð7Þ

where the temperature gradient vector rT j can be calculated from the imposed temperature boundary condition, 2. FE homogenization scheme

rT j ¼ T=Lj ;

The RVE based FE homogenization method is adopted in this paper to predict the effective thermal conductivity of the composites with discontinuous inclusions. Given a RVE x of the composites, the volume of the RVE x is V. The RVE x consists of the matrix x0 (the volume of V 0 and volume fraction of v 0 ¼ V 0 =V) and the inclusions x1 (the volume of V 1 and volume fraction of v 1 ¼ 1  v 0 ¼ V 1 =V). For any micro-field l (such as temperature

where T and Lj are the imposed temperature boundary condition and the side length of the RVE x in the direction xj , respectively. Note that in this research the preprocessing (including the generation of RVE and application of periodic boundary condition), FE analysis and postprocessing (including the generation of heat flux contour and derivation of effective thermal conductivity) concerning the evaluations on the effective thermal conductivities of the

ð8Þ

W. Tian et al. / International Journal of Heat and Mass Transfer 134 (2019) 735–751

composites with discontinuous inclusions are conducted in the commercial FE package ABAQUS [26]. 3. Periodic geometry model This research emphasizes the periodic boundary condition to evaluate the effective thermal conductivity of the composites with discontinuous inclusions, and thus the periodic RVEs characterizing the micro-structures of the composites with discontinuous inclusions are preferred. To generate the periodic RVEs of the composites with discontinuous inclusions, the modified RSA algorithm [23,24] is adopted. The modified RSA algorithm for RVE generation consists of adding the inclusions with the random positions (and orientations) sequentially into the matrix until the predefined inclusion volume fraction or a predefined number of attempts is reached. Physically, the intersection of any two inclusions is not acceptable in the RVE and the inclusions penetrating the RVE boundary should be cut and shifted to the opposite parallel surfaces of the RVE. Meanwhile, a minimum separation distance (taking 5:0% of the inclusion diameter as an example [23,27]) between the surfaces of any two inclusions and between the surfaces of the inclusions and the RVE is recommended to avoid that the meshed finite elements become distorted. Based on the modified RSA algorithm, the periodic RVEs for the composites reinforced by the randomly distributed spherical, the randomly distributed ellipsoidal and the randomly distributed and orientated cylindrical inclusions are generated and given in Figs. 1–3, respectively. The sizes of the spherical, ellipsoidal and cylindrical inclusions are listed in Table 1, and the sizes of the RVEs for the composites are critically selected as 50:0 lm  50:0 lm  50:0 lm; 75:0 lm  37:5 lm  37:5 lm and 100:0 lm  100:0 lm  100:0 lm, respectively, based on the previous researches [28–30]. The volume fractions of the inclusions in the RVEs are 30:0%; 30:0% and 12:5% for the composites reinforced by the randomly distributed spherical, the randomly distributed ellipsoidal and the randomly distributed and orientated cylindrical inclusions, respectively. The meshing of these periodic RVEs for the composites is performed by using ABAQUS built-in free meshing algorithm, and the mesh type is selected as 4-node linear heat transfer tetrahedron (DC3D4) finite element. The investigations on the effect of RVE mesh density on the prediction of the effective thermal conductivities of the composites reinforced by the randomly distributed spherical, ellipsoidal and cylindrical inclusions have been conducted and find that the mesh size of 1:25 lm for the RVEs of the composites reinforced by the spherical inclusions, the mesh size of 0:94 lm for the RVEs of the composites reinforced by the ellipsoidal inclusions and the mesh size of 2:5 lm for the RVEs of the composites reinforced by the cylindrical inclusions are sufficient to obtain convergent results.

737

To impose the periodic boundary condition, the periodic meshes on the surfaces of the generated periodic RVEs for the composites are required. In the ABAQUS, the function ”Copy Mesh Pattern” [26] is used to periodically mesh the surfaces of the RVEs for the composites, whose procedure is given as follows, (1) Copy the original matrix and inclusions and then obtain an identical matrix and inclusions, (2) Mesh the selected surface of the generated identical matrix and inclusions, (3) Copy the meshes of the selected surface of the generated identical matrix and inclusions to the corresponding surfaces of the original matrix and inclusions. 4. Periodic boundary condition and numerical algorithm The better approximations for the effective properties of composites are obtained when the periodic boundary condition instead of other boundary conditions is adopted, and the periodic boundary condition is thus considered as the first choice to predict the effective thermal conductivity of the composites with discontinuous inclusions in this paper. Given a periodic RVE x of the composites, the boundary of the RVE x is firstly decomposed into two opposing parts @ xþ and @ x under the constraint of @ x ¼ @ xþ [ @ x and 0 ¼ @ xþ \ @ x . Each material point xþ on the part @ xþ is associated with a unique material point x on the part @ x and the normal vectors at these two material points satisfy n ¼ nþ . The periodic boundary condition on the boundary @ x of the RVE x is then given as follows,

T þ  T  ¼ DT:

ð9Þ

In the ABAQUS, the periodic boundary condition can be imposed on the nodes located on the surfaces of the RVE x through a linear multi-point constraint equation [26], which requires a linear combination of nodal degree of freedom (DOF) equals zero,

a1 DPi þ a2 DQj þ    þ aN DRk ¼ 0;

ð10Þ

where DPi is the DOF i (i ¼ 11 for temperature) at the node P, and aN is the coefficient that defines the relative motion of the nodes. The periodic boundary condition in Eq. (9) is introduced to the linear multi-point constraint equation system through a Reference Point, which actually is not linked to any element in the RVE x but simply to provide the necessary DOFs to control the thermal transfer of the RVE x. The general form of the complete set of the linear multi-point constraint equation system is then represented as,


T kþ  T k  T RP ¼ 0 with T RP ¼ DT;

ð11Þ

where kþ and k are superscripts representing the relative periodic node sets on the opposite surfaces of the RVE x, and T RP is the tem-

Fig. 1. Periodic RVEs for the composites reinforced by randomly distributed spherical inclusions and their meshes: (a) Spherical inclusions, (b) Matrix and (c) Meshes of RVE.

738

W. Tian et al. / International Journal of Heat and Mass Transfer 134 (2019) 735–751

Fig. 2. Periodic RVEs for the composites reinforced by randomly distributed ellipsoidal inclusions and their meshes: (a) Ellipsoidal inclusions, (b) Matrix and (c) Meshes of RVE.

Fig. 3. Periodic RVEs for the composites reinforced by randomly distributed and orientated cylindrical inclusions and their meshes: (a) Cylindrical inclusions, (b) Matrix and (c) Meshes of RVE.

Table 1 Sizes of the spherical, ellipsoidal and cylindrical inclusions. Inclusion

Size (lm)

Spherical

Ellipsoidal

Cylindrical

Diameter

Major length

Minor length

Diameter

Length

10.4

20.76

10.38

5

50

perature perturbation carried on the Reference Point. This form of constraint equations is applied to all nodes in relevance to each other on the opposite surfaces of RVE x. Here the over-constraint should be avoided because that once a DOF has been used in a constraint equation, it cannot be used in another constraint equation in the finite element analysis. To impose the periodic boundary condition constraint equation (Eq. (11)) on the boundary of the RVE x, one of the prerequisites is to choose and categorize the nodes on the opposite surfaces of the RVE x, and the nodes on the opposite surfaces of the RVE x are categorized into three sets (as shown in Fig. 4(a)):  Set-I (Inner-Face): Face-BCC1B1, Face-ADD1A1, Face-ABCD, FaceA1B1C1D1, Face-ABB1A1 and Face-DCC1D1.  Set-II (Inner-Edge): Edge-AA1, Edge-BB1, Edge-CC1, Edge-DD1, Edge-AB, Edge-A1B1, Edge-D1C1, Edge-DC, Edge-AD, Edge-BC, Edge-B1C1 and Edge-A1D1.  Set-III (Vertex): Vertex-A, Vertex-B, Vertex-C, Vertex-D, VertexA1, Vertex-B1, Vertex-C1 and Vertex-D1. The nodes in Set-I are those nodes that belong to the faces excluding the nodes on the common edges and corners of the RVE x. The nodes in Set-II are those nodes that belong to the edges excluding the common corners of the RVE x. And the nodes in SetIII are eight corners of the RVE x (see Fig. 4(b)).

The detailed implementation equations of the periodic boundary condition for predicting the effective thermal conductivity of the composites are given as follows,  Set-I: Inner-face (1) Face-BCC1B1 and Face-ADD1A1:

T FaceBCC 1 B1  T FaceADD1 A1  T Rp1 ¼ 0:

ð12Þ

(2) Face-ABCD and Face-A1B1C1D1:

T FaceABCD  T FaceA1 B1 C 1 D1  T Rp2 ¼ 0:

ð13Þ

(3) Face-ABB1A1 and Face-DCC1D1:

T FaceABB1 A1  T FaceDCC 1 D1  T Rp3 ¼ 0:

ð14Þ

 Set II: Inner-edge (1) Edge-AA1, Edge-BB1, Edge-CC1 and Edge-DD1:

T EdgeBB1  T EdgeAA1  T Rp1 ¼ 0; T EdgeAA1  T EdgeDD1  T Rp3 ¼ 0; T EdgeCC 1  T EdgeDD1  T Rp1 ¼ 0: ð15Þ (2) Edge-AB, Edge-CD, Edge-C1D1 and Edge-A1B1:

T EdgeAB  T EdgeCD  T Rp3 ¼ 0; T EdgeCD  T EdgeC 1 D1  T Rp2 ¼ 0; T EdgeA1 B1  T EdgeC 1 D1  T Rp3 ¼ 0: ð16Þ

739

W. Tian et al. / International Journal of Heat and Mass Transfer 134 (2019) 735–751

Fig. 4. Planar and vertex notations of the RVE x and the categorized sets for the nodes on the opposite surfaces of the RVE x: (a) Planar and vertex notations and (b) Categorized sets.

(3) Edge-BC, Edge-DA, Edge-D1A1 and Edge-B1C1: EdgeBC

EdgeDA

Rp1

EdgeDA

T T T ¼ 0; T T EdgeB1 C 1  T EdgeD1 A1  T Rp1 ¼ 0:

T

EdgeD1 A1

T

Rp2

¼ 0; ð17Þ

T VerticeA  T VerticeD  T Rp3 ¼ 0; T VerticeC  T VerticeD  T Rp1 ¼ 0; T VerticeB  T VerticeC  T Rp3 ¼ 0: ð18Þ (2) Vertice-A1, Vertice-B1, Vertice-C1 and Vertice-D1: T VerticeB1  T VerticeC 1  T Rp1 ¼ 0; T VerticeA1  T VerticeD1  T Rp3 ¼ 0; T VerticeC 1  T VerticeD1  T Rp1 ¼ 0: ð19Þ

(3) Vertice-D and Vertice-D1:

ð20Þ

To determine the components kij (ki1 ; ki2 and ki3 with i ¼ 1; 2 and 3) of the effective thermal conductivity tensor of the composites, the RVE x has to be numerically simulated three times by choosing the proper magnitudes of T Rp1 ; T Rp2 and T Rp3 , which can be expressed as, (1) The components ki1 of the effective thermal conductivity tensor:


T Rp1 ¼ DT;

T Rp2 ¼ 0 and T Rp3 ¼ 0:

ð21Þ

(2) The components ki2 of the effective thermal conductivity tensor:

T Rp1 ¼ 0;

T Rp2 ¼ DT


and T Rp3 ¼ 0:

ð22Þ

(3) The components ki3 of the effective thermal conductivity tensor:

T Rp1 ¼ 0;


T Rp2 ¼ 0 and T Rp3 ¼ DT:

5. Numerical examples 5.1. Composites with spherical inclusions

 Set III: Vertex (1) Vertice-A, Vertice-B, Vertice-C and Vertice-D:

T VerticeD  T VerticeD1  T Rp2 ¼ 0:

preprocessing. The python script has been provided in the supplementary materials.

ð23Þ

It is noted that the periodic boundary condition constraint equations are imposed on the boundary of the RVE x in the ABAQUS by implementing a python script as a part of the automated

The firstly investigated example is the spherical SiC inclusions reinforced aluminum alloy matrix, and the inclusion volume fraction is 30:0%. The thermal conductivities of SiC inclusions and aluminum

alloy

matrix 1

are

393.0 W  m1  K1

and

173.0 W  m1  K , respectively. The interfaces between the inclusions and matrix are assumed to be prefect and can totally transfer the heat flux. First of all, the continuity of the thermal flux fields on the boundaries of the RVE for the composites is checked. The thermal flux components q11 ; q22 and q33 of the nodes on the corresponding opposite surfaces (i.e. Surface pair BCC1B1-ADD1A1, Surface pair ABCD-A1B1C1D1 and Surface pair ABB1A1-DCC1D1) of the RVE yielded from the predefined periodic boundary condition and their relative deviations are contoured in Fig. 5. It is found that the relative deviations of the thermal flux components q11 ; q22 and q33 of the node pairs located on the regions without the inclusion concentration of the surfaces of the RVE (see Fig. 6) are less than 10.0%, while those of the node pairs located on the regions with the inclusion concentration of the surfaces of the RVE are more than 30.0%. The percentages of the nodes on the corresponding opposite surfaces of the RVE with the thermal flux differences of 6 5:0%; 5:0  10:0% and P 10:0% are calculated and given in Fig. 7, which shows that the percentages for the thermal flux components q11 ; q22 and q33 of the nodes with the thermal flux difference of 6 10:0% are 52:5%; 53:7% and 50:0%. These relative deviations of the thermal flux components of the corresponding node pairs on the opposite surfaces of the RVE mainly result from the inclusion concentration on or near the surfaces of the RVE (in other words, the non-symmetric distribution of inclusions in the RVE), and partly from the numerical error due to the nodal extrapolation of the thermal flux (Note that the thermal flux of a node is extrapolated from the respective integration points using the shape functions of the element or elements sharing the nodes.). The thermal flux continuity can also be quantified by the calculation of the relative difference and average difference of the thermal flux for the node pairs that belong to the opposite surfaces of the RVE, respectively. The relative difference is the difference of

740

W. Tian et al. / International Journal of Heat and Mass Transfer 134 (2019) 735–751

Fig. 5. The thermal flux components q11 ; q22 and q33 of the nodes on the surface pairs BCC1B1-ADD1A1, ABCD-A1B1C1D1 and ABB1A1-DCC1D1 of the RVE and their relative deviations: (a)–(c) q11 on surface pair ABCD-A1B1C1D1 and their relative deviations, (d)–(f) q22 on surface pair ABCD-A1B1C1D1 and their relative deviations and (g)–(i) q33 on surface pair ABB1A1-DCC1D1 and their relative deviations.

Fig. 6. Inclusion concentration on the corresponding surfaces of the RVE for the composites with the spherical inclusions: (a) Surface BCC1B1 or ADD1A1, (b) Surface ABCD or A1B1C1D1 and (c) Surface ABB1A1 or DCC1D1.

741

W. Tian et al. / International Journal of Heat and Mass Transfer 134 (2019) 735–751

Fig. 7. Percentage of the node pairs with the thermal flux differences of 6 5:0%; 5:0  10:0% and P 10:0% on the opposite surfaces of the RVE. Set-1, Set2 and Set-3 represent q11 on surface-pair: BCC1B1-ADD1A1, q22 on surface-pair: ABCD-A1B1C1D1, q33 on surface-pair: ABB1A1-DCC1D1, respectively.

two maximum thermal flux of all the nodes on two opposite surfaces of the RVE, and the average difference is calculated by averaging the differences of the thermal flux of all the corresponding node pairs on the opposite surfaces of the RVE [30],

     Max qþ  Max q  i     i  E ¼ Min Max qþi ; Max qi r

 kþ  n q  qk  1X i i  ; and E ¼ k n k¼1 Max qkþ i ; qi a

ð24Þ where qkþ and qk represent the heat flux component qii of the cori i responding node pair k. qkþ and ql i represent the heat flux compoi nent qii of the nodes on the opposite surfaces of the RVE. Functions MaxðxÞ and MinðxÞ stand for the maximum and minimum of the absolute value of x, respectively. Fig. 8 shows the relative differences and average differences of the thermal flux components q11 ; q22 and q33 of the node pairs that belong to the opposite surface pairs of the RVE. Regarding the thermal flux components q11 ; q22 and q33 , these node-pairs have an average difference of 12:4%; 11:4% and 12:1%, respectively, and have a relative difference of 8:9%; 11:6% and 11:2%, respectively. To alleviate the effect of the non-symmetric distribution of the inclusions in the RVE, here the RVE with the symmetrically distributed inclusions (as shown in Fig. 9, marked RVE#2) is generated and investigated. The thermal flux components q11 ; q22 and q33 of the nodes on the corresponding opposite surfaces of the symmetrical RVE yielded from the predefined periodic boundary conditions and their relative deviations are contoured in Fig. 10. The relative deviations of the thermal flux components q11 ; q22 and q33 of the most node pairs located on the surfaces of the symmetrical RVEs are less than 5.0%. The relative differences and average

differences of the thermal flux components q11 ; q22 and q33 of the node pairs that belong to the opposite surfaces of the symmetrical RVE and the percentages of the nodes on the corresponding opposite surfaces of the symmetrical RVE with the thermal flux differences of 6 5:0%; 5:0  10:0% and P 10:0% are calculated and given in Fig. 11. Regarding the thermal flux components q11 ; q22 and q33 , these node-pairs have an average difference of 2:87%; 2:93% and 1:63%, respectively, and have a relative difference of 0:0348%; 0:0259% and 0:0159%, respectively. The percentages of the nodes with the thermal flux differences of 6 10:0% are 93:2%; 93:7% and 98:7%. These relative deviations of the thermal flux components of the corresponding node pairs on the opposite surfaces of the symmetrical RVE mainly result from the numerical error due to the nodal extrapolation of the thermal flux. The relative deviations, relative differences and average differences of the thermal flux components q11 ; q22 and q33 of the nodes on the corresponding opposite surfaces of those two types of RVE can be neglected. Therefore, the thermal flux continuity on the boundaries of the RVEs for the composites reinforced by the spherical inclusions can be considered as satisfied, and consequently the periodic boundary condition imposed by the numerical algorithm is validated. Imposing the periodic boundary condition with rT ¼ 0:50, the corresponding heat flux fields within the matrices and inclusions of the RVEs are shown in Figs. 12 and 13, which imply that heat flux fields within the matrices and inclusions of the unsymmetrical and symmetrical RVEs are completely different under the same boundary condition. The effective thermal conductivities of the composites reinforced by the spherical inclusions are homogenized (in W  m1  K1 ),

2

220:27

6 kc;1 ¼ 4 0:1923 2 kc;2

0:1152 223:53

6 ¼ 4 0:0170 0:0005

0:2037 0:2081

3

220:36

7 0:2113 5 and

0:1717

220:15

0:0169 0:0007 223:53 0:0002

3

7 0:0003 5: 223:53

Since the computed values of the off-diagonal entries in the 1st-3rd rows and columns of thermal conductivity tensors kc;1 and kc;2 are five orders of magnitude smaller than the other entries and, thus, are considered as null. In spite of the different heat flux fields within the matrices and inclusions of the RVEs, the predicted thermal conductivities of the composites reinforced by the spherical inclusions are identical. Because of the random and symmetrical distributions of the spherical inclusions, the homogenized thermal conductivities of the unsymmetrical and symmetrical RVEs are supposed to be isotropic. To investigate the isotropy of the effective thermal conductivity of the composites, the parameter defined in Ref. [30] is adopted here,

  min kii ; kP   Iii ¼ max kii ; kP

with

kP ¼ ðk11 þ k22 þ k33 Þ=3 and i ¼ 1; 2 and 3;

Fig. 8. Relative differences and average differences of the thermal flux components q11 ; q22 and q33 on the node pairs that belong to the opposite surfaces of the RVE.

ð25Þ

where Iii takes the value 1:0 for the composites with the isotropic thermal conductivity. Table 2 presents the values of Iii by substituting the entries in the thermal conductivity tensors kc;1 and kc;2 , and illustrates that the effective thermal conductivities of the composites reinforced by the randomly and symmetrically distributed spherical inclusions can be approximated to be perfectly isotropic. The homogenized thermal conductivities of the composites are firstly compared against those obtained by using the Lewis-Nielsen

742

W. Tian et al. / International Journal of Heat and Mass Transfer 134 (2019) 735–751

Fig. 9. Periodic RVEs for the composites reinforced by symmetrically distributed spherical inclusions and their meshes: (a) Spherical inclusions, (b) Matrix and (c) Meshes of RVE.

Fig. 10. The thermal flux components q11 ; q22 and q33 of the nodes on the surface pairs BCC1B1-ADD1A1, ABCD-A1B1C1D1 and ABB1A1-DCC1D1 of the symmetrical RVE and their relative deviations: (a)–(c) q11 on surface pair ABCD-A1B1C1D1 and their relative deviations, (d)–(f) q22 on surface pair ABCD-A1B1C1D1 and their relative deviations and (g)–(i) q33 on surface pair ABB1A1-DCC1D1 and their relative deviations.

model [7,8] to validate the effectiveness of the RVE based FE homogenization method with the periodic boundary condition. For the prediction of the effective thermal conductivity of compos-

ites, Lewis and Nielsen [7,8] derived a semi-theoretical model by a modification of the Halpin-Tsai equation [31] to include the effects of the inclusion shape and orientation,

W. Tian et al. / International Journal of Heat and Mass Transfer 134 (2019) 735–751

743

Fig. 11. (a) The relative differences and average differences of the thermal flux components q11 ; q22 and q33 on the node pairs that belong to the opposite surfaces of the symmetrical RVE, (b) Percentage of the node pairs with the thermal flux differences of 6 5:0%; 5:0  10:0% and P 10:0% on the opposite surfaces of the RVE.

Fig. 12. Heat flux fields within the matrix and spherical inclusions of the unsymmetrical RVE: (a), (b) and (c) q11 ; q22 and q33 in the matrix, (d), (e) and (f) q11 ; q22 and q33 in the spherical inclusions.

kc ¼ k0

1 þ v 1 AB ; 1  v 1 Bw

ð26Þ

metrical RVE and 205.56 W  m1  K1 for the symmetrical RVE) of the spherical SiC inclusions reinforced aluminum alloy matrix by the RVE based FE homogenization method with the periodic

where,

ðk1 =k0 Þ  1 B¼ ðk1 =k0 Þ þ A

and w ¼ 1 þ v 1

1  /m /2m

obtained from the experimental test [33]. The homogenized effective thermal conductivities (205.94 W  m1  K1 for the unsym-

! :

ð27Þ

where k0 and k1 are the thermal conductivities of the matrix and inclusions, and v1 is the inclusion volume fraction. A is 1.50, and /m is 0.637 for random close or 0.524 for simple cubic close, respectively, of the spherical inclusions [32]. The effective thermal conductivities of the composites with the different inclusion volume fractions predicted by the RVE based homogenization method with the periodic boundary condition and obtained by the Lewis-Nielsen model are listed in Table 3, which implies that the errors between the estimations of the Lewis-Nielsen model and the RVE based homogenization method with the periodic boundary condition are not of significance. The predicted thermal conductivity of the spherical SiC inclusions (the volume fraction of 0.20) reinforced aluminum alloy matrix by using the RVE based FE homogenization method with the periodic boundary conditions is compared against that

boundary condition agree well with that (203.25 W  m1  K1 ) obtained from the experimental test. Therefore, the RVE based on homogenization method with the periodic boundary condition can be used for the accurate prediction of the effective thermal conductivity of the composites reinforced by the spherical inclusions. 5.2. Composites with ellipsoidal inclusions The second example composite is the ellipsoidal SiC inclusions (with the volume fraction of 30:0%) reinforced aluminum alloy matrix. The thermal conductivities of SiC inclusions and aluminum alloy matrix are 393.0 W  m1  K1 and 173.0 W  m1  K1 , respectively. Analogously, the RVE with the symmetrically distributed ellipsoidal inclusions (as shown in Fig. 14, marked as RVE #2) is generated and investigated, and the interfaces between the ellipsoidal inclusions and matrix are assumed to totally transfer the heat flux.

744

W. Tian et al. / International Journal of Heat and Mass Transfer 134 (2019) 735–751

Fig. 13. Heat flux fields within the matrix and spherical inclusions of the symmetrical RVE: (a), (b) and (c) q11 ; q22 and q33 in the matrix, (d), (e) and (f) q11 ; q22 and q33 in the spherical inclusions.

Table 2 Isotropic parameters of the effective thermal conductivity of the composites reinforced by the randomly distributed spherical inclusions. RVE #1

Isotropic parameter

RVE #2

I11

I22

I33

I11

I22

I33

0.9996

0.9995

0.9995

1.0000

1.0000

1.0000

Table 3 Effective thermal conductivities of the composites with the randomly and uniformly distribtued spherical inclusions, respectively, predicted by the RVE based homogenization method and obtained by the Lewis-Nielsen model, respectively. RVE #1 (W  m1  K1 )

V 1 = 0.10 V 1 = 0.20 V 1 = 0.30

RVE #2 (W  m1  K1 )

FE method

L-N model

Error #1

FE method

L-N model

Error #2

188.91 205.94 220.26

188.68 206.95 228.57

0.12% 0.49 % 3.64%

188.70 205.56 223.53

189.24 209.53 235.45

0.29% 1.89% 5:05%

Fig. 14. Periodic RVEs for the composites reinforced by symmetrically distributed ellipsoidal inclusions and their meshes: (a) Ellipsoidal inclusions, (b) Matrix and (c) Meshes of RVE.

Figs. 15–17 contour the relative deviations of the thermal flux components q11 ; q22 and q33 of the node pairs on the corresponding opposite surfaces of the RVEs with the randomly and symmetrically distributed ellipsoidal inclusions yielded from the predefined periodic boundary condition. It is seen that the most node pairs

located on the regions without the inclusion concentration of the surfaces of the unsymmetrical RVE have the relative deviation of  10:0% for the thermal flux components q11 ; q22 and q33 , while the node pairs located on the regions with the inclusion concentration of the surfaces of the unsymmetrical RVE have the relative

W. Tian et al. / International Journal of Heat and Mass Transfer 134 (2019) 735–751

745

Fig. 15. The thermal flux components q11 ; q22 and q33 of the nodes on the surface pair BCC1B1-ADD1A1 of the unsymmetrical and symmetrical RVEs and their relative deviations: (a) Unsymmetrical RVE and (b) Symmetrical RVE.

Fig. 16. The thermal flux components q11 ; q22 and q33 of the nodes on the surface pair ABCD-A1B1C1D1 of the unsymmetrical and symmetrical RVEs and their relative deviations: (a) Unsymmetrical RVE and (b) Symmetrical RVE.

deviation of P 30:0% for the thermal flux components q11 ; q22 and q33 . Regarding the symmetrical RVE, the most node pairs located on its surfaces have the relative deviation of  5:0% for the thermal flux components q11 ; q22 and q33 . The relative differences and average differences of the thermal flux components q11 ; q22 and q33 of the node pairs that belong to the opposite surfaces of the RVEs are given Fig. 18, and the percentages of the nodes on the corresponding opposite surfaces of the RVEs with the thermal flux differences of 6 10:0%; 10:0  20:0% and P 20:0% are calculated and shown in Fig. 19. For the unsymmetrical RVE, the average differences and relative differences of the thermal flux components q11 ; q22 and q33 of the node pairs are not great than 16.7% and 15.8%, respectively, and the relative deviations of the thermal flux components q11 ; q22 and q33 for half the nodes pairs are less than 10:0%. For the symmetrical RVE, the average differences and relative differences of the thermal flux components q11 ; q22 and q33 of the node pairs are not great than 5.3% and 4.8%, respectively, and the relative deviation of the thermal flux components q11 ; q22 and q33 for 80.0% of the nodes pairs are less than 10:0%. Therefore, based on the analyses of the relative deviations, relative differences and average differences of the thermal flux

components q11 ; q22 and q33 of the nodes on the corresponding opposite surfaces of the unsymmetrical and symmetrical RVEs, the thermal flux continuity on the boundaries of the RVEs for the composites reinforced by the ellipsoidal inclusions can be considered as satisfied, and consequently the numerical algorithm for imposing the periodic boundary condition is verified. The corresponding heat flux fields within the matrices and inclusions of the RVEs under the identical periodic boundary condition (with rT ¼ 0:50) are shown in Figs. 20 and 21, which shows that the heat flux fields within the matrices and inclusions of the unsymmetrical and symmetrical RVEs are completely different. The resulted effective thermal conductivities of the composites reinforced by randomly distributed ellipsoidal inclusions are given (in W  m1  K1 ),

2

kc;1

kc;2

3 228:40 0:3431 0:2127 6 7 ¼ 4 0:3952 221:10 0:0630 5 and 0:1396 0:0419 221:11 2 3 225:46 0:0200 0:0027 6 7 ¼ 4 0:0237 221:36 0:0048 5: 0:0149 0:0082 221:36

746

W. Tian et al. / International Journal of Heat and Mass Transfer 134 (2019) 735–751

Fig. 17. The thermal flux components q11 ; q22 and q33 of the nodes on the surface pair ABB1A1-DCC1D1 of the unsymmetrical and symmetrical RVEs and their relative deviations: (a) Unsymmetrical RVE and (b) Symmetrical RVE.

Fig. 18. Relative differences and average differences of the thermal flux components q11 ; q22 and q33 on the node pairs that belong to the opposite surfaces of the unsymmetrical and symmetrical RVEs: (a) Unsymmetrical RVE and (b) Symmetrical RVE. Set-1, Set-2 and Set-3 represent q11 on surface-pair: BCC1B1-ADD1A1, q22 on surfacepair: ABCD-A1B1C1D1, q33 on surface-pair: ABB1A1-DCC1D1, respectively.

Fig. 19. Percentage of the node pairs with the thermal flux differences of 6 5:0%; 5:0  10:0% and P 10:0% on the opposite surfaces of the unsymmetrical and symmetrical RVEs: (a) Unsymmetrical RVE and (b) Symmetrical RVE.

747

W. Tian et al. / International Journal of Heat and Mass Transfer 134 (2019) 735–751

Fig. 20. Heat flux fields within the matrix and ellipsoidal inclusions of the unsymmetrical RVE: (a), (b) and (c) q11 ; q22 and q33 in the matrix, (d), (e) and (f) q11 ; q22 and q33 in the spherical inclusions.

Fig. 21. Heat flux fields within the matrix and ellipsoidal inclusions of the symmetrical RVE: (a), (b) and (c) q11 ; q22 and q33 in the matrix, (d), (e) and (f) q11 ; q22 and q33 in the spherical inclusions.

Table 4 Effective axial and transverse thermal conductivities of the composites with the ellipsoidal inclusions predicted by the RVE based homogenization method and obtained by the Lewis-Nielsen model, respectively. Axial conductivity (W  m1  K1 )

V 1 = 0.10 V 1 = 0.20 V 1 = 0.30

Transverse conductivity (W  m1  K1 )

FE method

L-N model

Error #1

FE method

L-N model

Error #2

190.96 209.92 228.40

191.01 210.04 230.22

0.03% 0.44% 0.79%

187.79 203.59 221.10

185.72 200.30 217.21

1.11% 1.64% 1.79%

In spite of the different heat flux fields within the matrices and inclusions of the RVEs, the difference of the predicted thermal conductivities of the composites reinforced by the ellipsoidal inclusions can be ignored, and the predicted thermal conductivities of the composites reinforced by the ellipsoidal inclusions can be approximated to be transversely isotropic. Hereafter, we thus just choose the RVE with the unsymmetrically distributed ellipsoidal inclusions to evaluate the effective thermal conductivities of the composites.

To validate the RVE based homogenization method with the periodic boundary condition, the homogenized effective thermal conductivities of the composites are compared against those obtained by using the Lewis-Nielsen model (Eqs. (26) and (27), in which A are 2L=D and 0.5 for predicting the axial and transverse thermal conductivities of the composites and /m is 0.82 for random close of the ellipsoidal inclusions [32].). The effective axial and transverse thermal conductivities of the composites with the inclusion volume fraction of 0.10, 0.20 and 0.30 predicted by the RVE

748

W. Tian et al. / International Journal of Heat and Mass Transfer 134 (2019) 735–751

Fig. 22. The relative deviations for the thermal flux components q11 ; q22 and q33 of the nodes on the surface pairs BCC1B1-ADD1A1, ABCD-A1B1C1D1 and ABB1A1-DCC1D1 of the RVE: (a) q11 on surface pair ABCD-A1B1C1D1, (b) q22 on surface pair ABCD-A1B1C1D1 and (c) q33 on surface pair ABB1A1-DCC1D1.

Fig. 23. Inclusion concentration on the corresponding surfaces of the RVE for the composites with the randomly distributed fibers: (a) Surface BCC1B1 or ADD1A1, (b) Surface ABCD or A1B1C1D1 and (c) Surface ABB1A1 or DCC1D1.

based homogenization method with the periodic boundary condition and obtained by the Lewis-Nielsen model are given in Table 4. The errors between the estimations of the Lewis-Nielsen model and the RVE based homogenization method are less than 1:80%, which validates the effectiveness of the RVE based homogenization method with the periodic boundary condition to predict the effective thermal conductivity of the composites reinforced by the ellipsoidal inclusions. 5.3. Composites with cylindrical inclusions The investigated composites here consist of short random cylindrical fibers (volume fraction of 12.5%) with the thermal conductivity of 32.0 W  m1  K1 and metal alloy matrix with the thermal conductivity of 72.0 W  m1  K1 . In the composites, the interfaces between the inclusions and matrix are assumed to be perfect.

The contours of the relative deviations for the thermal flux components q11 ; q22 and q33 of the corresponding node pairs on the opposite surfaces of the RVE yielded from the predefined periodic boundary condition are given in Fig. 22. It is found that the thermal flux differences of the node pairs located on the regions without the inclusion concentration of the surfaces of the RVE (see Fig. 23) are less than 10.0% and those of the node pairs located on the regions with the inclusion concentration of the surfaces of the RVE are more than 30.0%. As aforementioned, these thermal flux differences of the corresponding node pairs on the opposite surfaces of the RVE mainly result from the unsymmetrical microstructures of the RVE, and partly from the numerical errors due to the extrapolation of the thermal flux at the nodes on the surfaces of the RVE. The percentages of the node pairs on the corresponding opposite surfaces of the RVE with the thermal flux differences of 6 5:0%; 5:0  10:0% and P 10:0% are given in Fig. 24(a), which

Fig. 24. (a) Percentage of the node pairs with the thermal flux differences of 6 5:0%; 5:0  10:0% and P 10:0% on the opposite surfaces of the RVE, (b) Relative differences and average differences of the thermal flux components q11 ; q22 and q33 on the node pairs that belong to the opposite surfaces of the RVE. Set-1, Set-2 and Set-3 represent q11 on surface-pair: BCC1B1-ADD1A1, q22 on surface-pair: ABCD-A1B1C1D1, q33 on surface-pair: ABB1A1-DCC1D1, respectively.

749

W. Tian et al. / International Journal of Heat and Mass Transfer 134 (2019) 735–751

Fig. 25. Heat flux fields within the matrix and inclusions of the RVE: (a), (b) and (c) q11 ; q22 and q33 in the matrix, (d), (e) and (f) q11 ; q22 and q33 in the inclusions.

implies that the percentage of the node pairs with the thermal flux difference of 6 10:0% is great than 75:0%. Fig. 24(b) shows the relative differences and average differences of the thermal flux components q11 ; q22 and q33 of the node pairs that belong to the opposite surfaces of the RVE. Regarding the thermal flux component q11 , these node-pairs have an average difference of 7:90% yielding a relative difference of 5:09%. Regarding the thermal flux component q22 , these node-pairs have an average difference of 7:55% yielding a relative difference of 5:08%. And regarding the thermal flux component q33 , these node-pairs have an average difference of 7:30% yielding relative difference of 3:25%. Therefore, it is concluded that the thermal flux continuity on the boundaries of the RVE for the composites reinforced by the cylindrical inclusions can be considered to be satisfied, and consequently the periodic boundary condition imposed by the numerical algorithm in Section 4 is verified. By imposing the periodic boundary condition (with rT ¼ 0:50), the corresponding heat flux fields within the matrices and inclusions of the RVE are obtained and shown in Fig. 25, and the effective thermal conductivities of the composites are homogenized and given (in W  m1  K1 ) as follows,

2

66:52

6 kc ¼ 4 0:0331 0:0149

0:0309 66:55

0:0140

3

7 0:0502 5:

0:0498

66:59

The values of the off-diagonal entries in the 1st-3rd rows and columns of the thermal conductivity tensor kc are considered as null because that they are four orders of magnitude smaller than the other entries. Because of the random distribution of the short inclusions, the effective thermal conductivities of the composites are expectedly isotropic (i.e. I11 ¼ 0:9995; I22 ¼ 1:0000 and I33 ¼ 0:9994) [23,24,30]. The effective thermal conductivities of the composites evaluated by the RVE based FE homogenization method with the periodic boundary condition are compared against those obtained based on the Lewis-Nielsen model [7,8] (Eqs. (26) and (27)) to validate the effectiveness of the RVE based FE homogenization method with the periodic boundary conditions. Here, A and /m in Eqs. (26) and (27) are 4.93 and 0.52, respectively, for randomly distributed cylindrical inclusions [32]. Fig. 26 illustrates the effective thermal conductivities of the composites with the different inclusion volume fractions predicted by the RVE based homogenization method with the periodic boundary condition and by the LewisNielsen model. It is seen that the predictions of the RVE based homogenization method with the periodic boundary condition agree well with those predicted by using the Lewis-Nielsen model. Meanwhile, the RVE based homogenization method with the periodic boundary condition is compared against the experimental test [34], as listed in Table 5, which shows the error between the predictions of the RVE based homogenization method with the periodic boundary condition and the Lewis-Nielsen model is 5:28%.

Table 5 Comparison of the effective thermal conductivities of the composites obtained by the RVE based homogenization method and the experimental test [34] (Here, the effective thermal conductivity of the composites is calculated by averaging the effective thermal conductivities of the composites parallel to the pressed plane and perpendicular to the pressed plane.). Thermal conductivity (W  m1  K1 ) Fig. 26. The homogenized thermal conductivities of the composites with the different inclusion volume fraction obtained by the RVE based homogenization method with the periodic boundary condition and the Lewis-Nielsen model.

Value

FE homogenization

Experimental test

Error

102.24

97.11

5:28%

750

W. Tian et al. / International Journal of Heat and Mass Transfer 134 (2019) 735–751

The effectiveness of the RVE based homogenization method with the periodic boundary condition to predict the effective thermal conductivity of the composites reinforced by the cylindrical inclusions is thus validated.

5.4. Some remarks (1) The inclusion distribution in the RVE has a significant impact on the continuity of the heat flux of the corresponding node pairs on the opposite surfaces of the RVE. For the RVEs with the symmetrical inclusion distribution, the relative differences and average differences of the thermal flux components of the node pairs on their opposite surfaces are lower than other types of RVE. However, the inclusion distribution in the RVE has a neglectable effect on the effective thermal conductivity of composites. (2) The ratio between the thermal conductivities of matrix and inclusions affect the continuity of the heat flux of the corresponding node pairs on the opposite surfaces of the RVE, and the ratio between the thermal conductivities of matrix and inclusions close to 1.0 provides the better continuity of the heat flux of the corresponding node pairs on the opposite surfaces of the RVE. (3) The proposed periodic boundary condition and its numerical algorithm can be used for numerically predicting the effective thermal conductivities of other kinds of composites, such as 3-D four directional braided composites and 3-D woven composites.

6. Conclusions In this paper, the RVE based FE homogenization method and the periodic boundary condition are introduced to evaluate the effective thermal conductivity of the composites with discontinuous inclusions, and the research focus is on the periodic boundary condition, its numerical implementation algorithm and validation. The heat flux continuity of the corresponding node pairs on the opposite surfaces of the RVE is checked, and the validation of the RVE based FE homogenization method with the periodic boundary condition is conducted by comparing with the Lewis-Nielsen model and the experimental test. The main conclusions are summarized as follows: (1) The heat flux continuity of the corresponding node pairs on the opposite surfaces of the RVE of composites can be guaranteed, which validates that the periodic boundary condition is imposed properly by the proposed numerical implementation algorithm. (2) Through comparing with the Lewis-Nielsen model and the experimental test, the RVE based FE homogenization method with the periodic boundary condition is verified in terms of accurate evaluation of the effective thermal conductivity of the composites with discontinuous inclusions. (3) For the composites reinforced by randomly distributed spherical and cylindrical inclusions, their effective thermal conductivities can be approximated as macroscopically isotropic, while the effective thermal conductivity of the composites reinforced by the unidirectional ellipsoidal inclusions can be approximated as transversely isotropic.

Conflict of interest The authors declare here that there are no known conflicts of interest with other people or organizations.

Data availability The raw/processed data required to reproduce these findings cannot be shared at this time due to technical or time limitations. Acknowledgement The authors would like to thank the financial support from the National Natural Science Foundation of China (Nos. 51821091 and 51772245) and the project from The Hong Kong Polytechnic University (No.1-YW3H) to support the post-doc research of Dr. Wenlong Tian in the university. Appendix A. Supplementary material Here, the python codes for the auto-preprocessing of FE simulation in the ABAQUS to predict the effective thermal conductivities of composites are attached, which mainly contain some modules for the material properties assignment, analysis step creation, mesh generation, node sets generation and boundary condition application. Supplementary data associated with this article can be found, in the online version, at https://doi.org/10.1016/j.ijheatmasstransfer.2019.01.072. References [1] W. Tian, L. Qi, C. Su, J. Zhou, Mean-field homogenization based approach to evaluate macroscopic coefficients of thermal expansion of composite materials, Int. J. Heat Mass Transf. 102 (2016) 1321–1333. [2] W. Tian, L. Qi, X. Chao, L. Ju, S. Li, J. Liang, Experimental and multi-scale numerical evaluations for effective mechanical properties of 2-D Cf/Mg composites, Compos. Struct. 189 (2018) 1–8. [3] I. Doghri, M.I. El Ghezal, L. Adam, Finite strain mean-field homogenization of composite materials with hyperelastic-plastic constituents, Int. J. Plast. 81 (2016) 40–62. [4] F. Naya, M. Herráez, C. Lopes, C. González, S.V. der Veen, F. Pons, Computational micromechanics of fiber kinking in unidirectional FRP under different environmental conditions, Compos. Sci. Technol. 144 (2017) 26–35. [5] J. Zhai, S. Cheng, T. Zeng, Z. Wang, D. Fang, Extended multiscale FE approach for steady-state heat conduction analysis of 3D braided composites, Compos. Sci. Technol. 151 (2017) 317–324. [6] I.L. Ngo, C. Byon, A generalized correlation for predicting the thermal conductivity of composites with heterogeneous nanofillers, Int. J. Heat Mass Transf. 90 (2015) 894–899. [7] R.C. Progelhof, J.L. Throne, R.R. Ruetsch, Methods for predicting the thermal conductivity of composite systems: a review, Polym. Eng. Sci. 16 (1970) 615– 625. [8] H. Yu, D. Heider, S. Advani, Prediction of effective through-thickness thermal conductivity of woven fabric reinforced composites with embedded particles, Compos. Struct. 127 (2015) 132–140. [9] S. Lee, J. Lee, B. Ryu, S. Ryu, A micromechanics-based analytical solution for the effective thermal conductivity of composites with orthotropic matrices and interfacial thermal resistance, Sci. Rep. 8 (2018) 7266. [10] H.L. Quang, Determination of the effective conductivity of composites with spherical and spheroidal anisotropic particles and imperfect interfaces, Int. J. Heat Mass Transf. 95 (2016) 162–183. [11] S.Y. Kim, Y.J. Noh, J. Yu, Thermal conductivity of graphene nanoplatelets filled composites fabricated by solvent-free processing for the excellent filler dispersion and a theoretical approach for the composites containing the geometrized fillers, Compos. Part A: Appl. Sci. Manuf. 69 (2015) 219–225. [12] 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. [13] M. Kulkarni, R. Brady, A model of global thermal conductivity in laminated carbon/carbon composites, Compos. Sci. Technol. 57 (1997) 277–285. [14] Y. Benveniste, T. Chen, G.J. Dvorak, The effective thermal conductivity of composites reinforced by coated cylindrically orthotropic fibers, J. Appl. Phys. 67 (1990) 2878–2884. [15] A.M. Thiele, A. Kumar, G. Sant, L. Pilon, Effective thermal conductivity of threecomponent composites containing spherical capsules, Int. J. Heat Mass Transf. 73 (2014) 177–185. [16] I. Ioannou, A. Hodzic, I.M. Gitman, Numerical investigation of thermal and thermo-mechanical effective properties for short fibre reinforced composite, Appl. Compos. Mater. 24 (2017) 999–1009. [17] 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 (16) (2010) 2276–2283.

W. Tian et al. / International Journal of Heat and Mass Transfer 134 (2019) 735–751 [18] L. Jiang, G. Xu, S. Cheng, X. Lu, T. Zeng, Predicting the thermal conductivity and temperature distribution in 3D braided composites, Compos. Struct. 108 (2014) 578–583. [19] S. Bargmann, B. Klusemann, J. Markmann, J.E. Schnabel, K. Schneider, C. Soyarslan, J. Wilmers, Generation of 3D representative volume elements for heterogeneous materials: a review, Prog. Mater Sci. 96 (2018) 322–384. [20] W.M. Harris, W.K. Chiu, Determining the representative volume element size for three-dimensional microstructural material characterization. Part I: predictive models, J. Power Sources 282 (2015) 552–561. [21] J. Gou, H. Zhang, Y. Dai, S. Li, W. Tao, Numerical prediction of effective thermal conductivities of 3D four-directional braided composites, Compos. Struct. 125 (2015) 499–508. [22] K. Dong, J. Zhang, L. Jin, B. Gu, B. Sun, Multi-scale finite element analyses on the thermal conductive behaviors of 3D braided composites, Compos. Struct. 143 (2016) 9–22. [23] 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. [24] K. Babu, P. Mohite, C. Upadhyay, Development of an RVE and its stiffness predictions based on mathematical homogenization theory for short fibre composites, Int. J. Solids Struct. 130–131 (2018) 80–104. [25] H. Thoemen, T. Walther, A. Wiegmann, 3D simulation of macroscopic heat and mass transfer properties from the microstructure of wood fibre networks, Compos. Sci. Technol. 68 (3) (2008) 608–616.

751

[26] ABAQUS/Standard, ABAQUS scripting/analysis user’s manual, Version 6.14-3, RI: ABAQUS Inc., 2014. [27] H. Bohm, A. Eckschlager, W. Han, Multi-inclusion unit cell models for metal matrix composites with randomly oriented discontinuous reinforcements, Comput. Mater. Sci. 25 (2002) 42–53. [28] C. González, J. Segurado, J. LLorca, Numerical simulation of elasto-plastic deformation of composites: evolution of stress microfields and implications for homogenization models, J. Mech. Phys. Solids 52 (7) (2004) 1573–1593. [29] O. Pierard, C. González, J. Segurado, J. LLorca, I. Doghri, Micromechanics of elasto-plastic materials reinforced with ellipsoidal inclusions, Int. J. Solids Struct. 44 (21) (2007) 6945–6962. [30] 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. [31] J. Halpin, Stiffness and expansion estimates for oriented short fiber composites, J. Compos. Mater. 3 (1969) 732–734. [32] D. Kumlutas, I.H. Tavman, A numerical and experimental study on thermal conductivity of particle filled polymer composites, J. Thermoplast. Compos. Mater. 19 (2006) 441–455. [33] 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. [34] K. Asano, H. Yoneda, Y. Agari, Microstructure and thermal properties of squeeze cast aluminum alloy composite reinforced with short potassium titanate fiber, Mater. Trans. 49 (2008) 2664–2669.