A micromechanics-based analytical solution for the effective thermal conductivity of composites with orthotropic matrices and interfacial thermal resistance

We obtained an analytical solution for the effective thermal conductivity of composites composed of orthotropic matrices and spherical inhomogeneities with interfacial thermal resistance using a micromechanics-based homogenization. We derived the closed form of a modified Eshelby tensor as a function of the interfacial thermal resistance. We then predicted the heat flux of a single inhomogeneity in the infinite media based on the modified Eshelby tensor, which was validated against the numerical results obtained from the finite element analysis. Based on the modified Eshelby tensor and the localization tensor accounting for the interfacial resistance, we derived an analytical expression for the effective thermal conductivity tensor for the composites by a mean-field approach called the Mori-Tanaka method. Our analytical prediction matched very well with the effective thermal conductivity obtained from finite element analysis with up to 10% inhomogeneity volume fraction.


Results and Discussion
Effective Thermal Conductivity in the Absence of Interfacial Resistance. In the steady state, the governing equation for heat conduction is written as, where g is a heat source, T is temperature field, and K 0 is a symmetric second order thermal conductivity tensor. The repeated small letter indices represent the dummy indices that imply summation over all the values from 1 to 3. The Green's function G(x − y) in the steady state heat conduction equation is defined as the temperature field at a position x in the presence of a unit heat source at another position y in an infinite medium, , the Greens' function reduces to a well-known function that is inversely proportional to the distance, .
π − x y k 1 4 0 In the absence of interfacial thermal resistance, the Green's function is introduced to derive the Eshelby tensor S ik of an eigen-intensity problem that relates the intensity field = −∇ e T and the eigen-intensity field e * within the inclusion as 30 where V is the volume of an inclusion whose thermal conductivity is identical to the matrix (see Fig. 1(a)). The eigen-intensity field e * can be considered as a fictitious temperature gradient field produced without external heat flux for an isolated inclusion, and the intensity field e is the temperature gradient within the inclusion when it is embedded in the matrix. With the Eshelby tensor, we can solve the eigen heat flux problem, which relates the heat flux within the inclusion q and eigen heat flux q * as, where the C is known as the conjugate Eshelby tensor,  is a diagonal matrix with its diagonal elements being the half the length of the principal axes of an ellipsoidal inclusion whose volume is defined as . For example, for the spherical inclusion, a can be expressed as δ = a R ij ij , where R is the radius of inclusion. In this study, we consider the Eshelby tensor for a spherical inclusion in the orthotropic matrix with three independent thermal conductivity coefficients. The thermal conductivity tensors of matrix (K 0 ) and inhomogeneity (K 1 ) are defined as follows, Now, we derive a closed form expression of S ij for a spherical inclusion by plugging the thermal conductivity tensor K 0 into the Eq. (7) with δ = a R ij ij as,  (see Supplementary Fig. S1). The three independent values (S 11 , S 22 , S 33 ) are plotted in terms of k k / 1 2 and k k / 1 3 , and we validate our analytical solutions against the numerical evaluation of Eq. (4) (see Supplementary Fig. S1).
It has been proven that the heat flux q inh within an ellipsoidal inhomogeneity embedded in an infinite matrix under the presence of a constant external far field heat flux q ext is uniform, and so is the Eshelby tensor S, regardless of the materials symmetry of the matrix 29,30 (see Fig. 1(b)). Heat fluxes, q inh and q ext , are related by the localization tensor B 22,30 Using the Eshelby tensor in the orthotropic matrix in Eq. (10), we predict the heat flux within the inhomogeneity with any arbitrary thermal conductivity tensor K 1 . For an inhomogeneity with isotropic or cubic symmetry whose thermal conductivity tensor is given as , the heat flux expression can be simplified as Our solution is validated against the numerical calculations based on FEA where we consider a single inhomogeneity surrounded by an orthotropic medium in a cubic shape. We set the edge of the medium to be 15 times larger than the diameter of the inhomogeneity to reasonably describe the infinite medium 22 as depicted in Fig. 2. The FEA is performed by COMSOL software with a total of 416,606 tetrahedral quadratic elements in the matrix and inhomogeneity. In this simulation, the unit heat flux boundary conditions are considered in the x, y, and z directions to study the effect of the anisotropy ( . We carry out the calcula- and κ = 10 (W/mK). As expected, the calculated heat flux within the inhomogeneity is uniform and dependent on the external heat flux direction (see Fig. 2(b)), and matches very well with the theoretical prediction (see Fig. 3).
The effective thermal conductivity of a composite with multiple inhomogeneities can be predicted by considering the interaction between the inhomogeneities. In a mean field approach, as in the Mori-Tanaka method, the heat flux within the inhomogeneity is related to the average heat flux within the matrix. The Mori-Tanaka model is known to predict effective properties well at a relatively low inhomogeneity volume concentration (<20%) and is more convenient than the self-consistent method, which relies on a nonlinear implicit equation. In the absence of the interfacial resistance, the effective thermal conductivity based on the Mori-Tanaka method can be obtained as 19,30

Effective Thermal Conductivity in the Presence of Interfacial Resistance.
We now turn our attention to the realistic system, where interfacial thermal resistance is present 31,32 . The interfacial thermal resistance α is defined as  where T out and T in refer to temperatures at the outer and inner surfaces of the interface, respectively, q is the heat flux at the interface, and the n is the outward surface normal vector (see Fig. 1(c,d)). The SI unit of interfacial thermal resistance α is [m 2 K/W]. The interfacial resistance augments an additional surface integration term in the eigen-intensity problem, as follows, It has been found that the heat flux within a spherical inclusion is uniform in the presence of interfacial thermal resistance 22 . Although a previous study 22 claims that the heat flux within an elliptical inclusion is also uniform in the presence of interfacial thermal resistance, our numerical tests reveal that it is non-uniform (See Supplementary Note 4). Because the intensity field in the spherical inclusion is uniform, we can simplify Eq. (14) as follows:  . After substituting the conjugate Eshelby tensor and the modified conjugate Eshelby tensor into the modified localization tensor equation, the heat flux within the inhomogeneity in the presence of interfacial thermal resistance can be obtained as follows, The theoretical predictions of heat fluxes along three directions match very well with the FEA calculation results with the same boundary conditions for a wide range of interfacial thermal resistance α, as shown in Fig. 4. At the limit of zero resistance, i.e., α → 0, the heat flux within the inhomogeneity converges to the zero interfacial resistance case depicted in Fig. 3. At the opposite limit, as α → ∞, the heat flux within the inhomogeneity reduces to zero, which implies that the infinite interfacial thermal resistance is equivalent to a void or a perfect thermal insulator as an inhomogeneity. Hence, the effective conductivity approaches that of a porous medium.
We then derive a closed form solution for the effective thermal conductivity based on a mean field micromechanics model, the Mori-Tanaka method. Following the previous study 22 , the effective thermal conductivity of a composite with interfacial thermal resistance can be determined by Although we consider a thermally isotropic inhomogeneity in the final solution for the sake of simplicity, the effective thermal conductivity with anisotropic inhomogeneity can be easily predicted by using an anisotropic K 1 .
We plot the effective thermal conductivities, K 1 eff , K 2 eff and K 3 eff , as a function of the inhomogeneity's volume fraction c 1 and the interfacial thermal resistance α, as shown in Fig. 5. At the limit of zero interfacial resistance, i.e., α → 0, Eq. (20) becomes identical to Eq. (12), where we assumed no interfacial resistance. At the opposite limit of α → ∞, because no heat flux is permitted within the inhomogeneity, the effective conductivity tensor converges to that of a porous medium,  and κ = 10 (W/mK), and the radius of the particle is 1 (mm).  , κ = 10 (W/mK), and = R mm 1 . Because the inhomogeneity is more conductive than the matrix (κ > k k k , , 3 ), the effective thermal conductivity increases with the volume fraction in the range where the interfacial thermal resistance is low. However, at high enough interfacial thermal resistance, the effective thermal conductivity decreases with the volume fraction because the heat flux through the inhomogeneity is significantly limited. We calculate the critical interfacial thermal resistance that makes the effective thermal conductivity of the composite K I eff identical to the thermal conductivity of the matrix k I , as follows: I critical When the interfacial thermal resistance is higher than the critical value, the effective thermal conductivity of the composite decreases as the volume fraction increases, even though the inhomogeneity is thermally more conductive than the matrix. We note that the critical interfacial resistance decreases for smaller particles, because the interface area increases with the decreasing size of the inhomogeneity for a given volume fraction. We validate our analytical solution presented in Eq. (20) by comparing it with the effective thermal conductivity calculated numerically by FEA, as depicted in Fig. 6. We obtain the effective conductivity by averaging the results from 10 independent RVEs, each containing multiple spherical inhomogeneities that are randomly distributed within a cube (see Fig. 7). We assign a uniform temperature boundary condition at two parallel surfaces while applying a periodic boundary conditions along the other two directions. We then compute the heat flux to obtain the effective thermal conductivity of each RVE. As shown in Fig. 6(a-c), the FEA results match very well with our solution for up to 10% of the inhomogeneity volume fraction for a wide range of interfacial thermal resistances. We also investigate the effect of the inhomogeneity's size at a fixed volume fraction of 5% and an interfacial thermal resistance (10 −3 m 2 K/W), as depicted in Fig. 6(d-f). When the interfacial resistance is absent, both theoretical predictions and numerical results find that the effective conductivity is independent of the size of inhogeneity. In contrast, in the presence of interfacial resistance, the effective thermal conductivity decreases as the radius of inhomogeneity decreases, because the interface fraction is bigger for smaller inhomogeneities for a fixed volume.

Conclusion
In conclusion, we have investigated the heat conduction problem of composites with orthotropic matrices and a spherical inhomogeneities in the presence of interfacial thermal resistance. We derive the modified Eshelby tensor of the eigen-intensity problem as well as the effective thermal conductivity based on a micromechanics approach by considering the interfacial thermal resistance effect, and validate our solution against FEA calculation results. We also demonstrate that the effective conductivity solution has the correct limiting behaviour at both the zero and infinite interfacial thermal resistance limit. The solution in the present paper is applicable to the composites with either transversely isotropic or isotropic matrices and an inhomogeneity with an arbitrary thermal conductivity tensor. We plan to extend the present study by considering the size effects of nanoscale inclusions for nanocomposites 49 in obtaining an analytic solution and by coupling molecular dynamics simulations of ceramic composites or polymer composites with the analytic solution. We believe that our study can provide an effective and accurate way of predicting the thermal conduction of composites, and it can be applied to better design technologically important materials such as polymer-based composite and thermoelectric materials. Mesh configuration of inhomogeneities in representative volume element for FEA. The volume fraction is 5% and the particle radius is 1 (mm).