A Multi-Scale, Semi-Analytical Model for Transient Heat Transfer in a Nano-Composite Containing Spherical Inclusions

This paper presents a semi-analytical model for transient heat conduction in a composite material reinforced with small spherical inclusions. Essential to the derivation of the model is the assumption that the size of the inclusions is much smaller than the length scale characterizing the macroscopic problem. An interfacial thermal resistance is also present between the two phases. During heating, the inclusions are treated as heat sinks within the matrix, with the coupling provided by the boundary conditions at the surface of the embedded particles. Application of Duhamel’s Theorem at the particle scale provides the local relationship between the temperature profile in a particle and the matrix that surrounds it. A simple spatial discretization at the macro-scale leads to an easily solvable system of coupled Ordinary Differential Equations for the matrix temperature, particle surface temperature and a series of ψ-terms related to the heat exchange between phases. The interfacial thermal resistance between the two phases can lead to the particle temperature lagging behind that of the surrounding matrix. The resulting transient response of the matrix temperature cannot be reproduced by a material with a single effective thermal conductivity. In the case where transient methods are used to determine effective thermal conductivity, this transient response may introduce errors into the measurement. Article info Received: 4 September 2018 Received in revised form: 20 October 2018 Accepted: 24 December 2018


Introduction
Prediction of the effective properties of composite materials is a prerequisite for their effective use. Derivations of effective properties of composites go back to Maxwell's initial calculations in 1873 for composites comprised of dispersed spheres [1]. In 1986, Davis proposed a method to calculate an effective thermal conductivity with spherical inclusions that extended Maxwell's work [2]. The next year, Hasselman proposed a new effective medium approach that accounted for a resistance to heat flow at the interface [3], before the approach was further generalized by Nan in 1997 [4]. More recent works have used micro mechanics-based solutions to determine effective thermal conductivity of materials [5] and the review by Zhai et al. provides a good overview of the relation between the different model approaches [6]. A commonality between these studies is that they consider steady state conduction through the multi-phase material and determine effective material properties that yield an identical heat flux for a given thermal gradient. In order to describe the transient response of such systems, the heat exchange between he phases must also be taken into account [7]. As techniques such as the Transient Hot Wire Method are commonly used to determine the effective thermal conductivity of a material [8], these effects may become a source of error in the measurement of effective properties.
Many of the works in the literature have been motivated by the desire to increase the thermal conductivity via conductive fillers for applications such as electronics' packaging [9]. In this work we focus on materials with low-conductivity inclusions for insulative applications [10,11]. We pay particular attention to the multi-scale aspect of composites and derive a semi-analytical model capable of linking the micro-scale heat exchange between the particles and the matrix with the macro-scale heat transport. Our model, being semi-analytical, does not require a complex/multi-scale mesh but instead requires only a simple discretization of the macro-scale.

Governing equations
Consider the composite slab, shown in Fig. 1, consisting of a matrix (phase m) reinforced with spherical inclusions (phase p) of radius R, where each phase has its own thermal conductivity For an averaging volume, V, large enough that it contains a representative sample of each phase, yet small enough that the temperature does not vary significantly within it, the volume averaged heat equation in phase m is given by Glatzmaier and Ramirez in Eq. 1 [12].
n mp is the unit vector pointing into the matrix and A mp is the total interfacial area contained within the volume. The integral term over n mp T represents the reduction in conduction due to the tortuosity of the matrix [12]. This is considered to result in the conductivity being reduced to that of a matrix containing voids, K m,voids = 2K m (1 -V p )/(2 + V p ) as evaluated by Hasselman's approach. The integral term over represents the heat exchange between the matrix and the particles. It can then be evaluated in Eq. 2, where the interfacial thermal resistance is γ (in ) and the temperature of the particle surface is T S = T p(r=R) .
Equation 1 now takes the form of Eq. 3 where the particles' effect on the heat transfer can be seen in both terms via the reduction matrix conductivity and their treatment as heat sources/sinks.
The heat equation inside an individual particle is given by Eq. 4. The surface boundary condition of heat exchange to the matrix is given in Eq. 5 and provides the coupling of Eqs. 3 and 4 via Eq. 2.
In this paper, the term is solved by first applying Duhamel's Theorem to solve Eq. 4 subject to the time-dependent boundary condition, T S . If Z(r,t), Eq. 6, is the time evolution of the particle temperature for an initial temperature of zero and a surface temperature of unity, then for a surface temperature of T S (t), the particle temperature is given by Eq. 7 by Carslaw and Jaeger [13].
The thermal diffusivity, κ p , is given by . (2) (K m or K p ), density (ρ m or ρ p ), specific heat capacity (c m or c p ) and volume fraction (V m or V p , where V p +V m = 1).
The functions ψ n (t) are described by Eq. 10, or in differential form in Eq. 11. Using Eq. 9, the partial derivative term in Eq. 5 can be evaluated to obtain Eq. 12.
From Eq. 11, it is evident that the ψ n functions are first order, subject to the forcing function .
For large enough values of the n 2 term, ψ n can be reasonably approximated by its quasi-steady state approximation. This allows the summation to be split into two parts, the first n 0 terms are given by the solutions to the ODEs in Eq. 11 and the second set of terms is approximated by the quasi-steady state approximation, leading to Eq. 13.
Equations 3, 11 and 13 can be rearranged into a set of n 0 +2 differential equations with respect to time, Eqs. 14-16. Equation 14 can then be discretized in the spatial dimension into N nodes, leading to a system of N(n 0 +2) ODEs that can be solved using standard numerical methods.

Numerical Application
The test-case involves a 1 cm thick slab, initially at T int where the surface temperature is T int + 1 K for t > 0 s. The properties of the particles and   In the limit as the thermal resistance tends to zero, the Maxwell and Hasselman estimates converge in Fig. 3 while the proposed model underestimates the conductivity of the composite. This is expected as the constant surface temperature at the particle-level used in this model implies no thermal gradient and therefore no heat flux through the particles themselves. When a small thermal resistance is present, a temperature lag appears between the particle and the matrix, as shown in Fig. 4. As the particles have a lower temperature than the matrix, they behave like heat sinks and slow the temperature increase at the center of the slab, leading to an apparent decrease in thermal conductivity. The temperature response of a homogenous slab having the K effective that minimizes the sum of squared error (SSE) between the predicted temperature and the matrix temperature is also included in the figure. Fig. 4. When a small thermal resistance is present, the particle temperature lags behind that of the matrix (γ = 10 -2 m² • K/W, V p = 0.1, R = 10 µm).
An interesting result of these simulations is that the particles also induce a spreading, or smearing in the temperature response. This is more clearly illustrated in the difference between the temperatures at the center of three different homogenous slabs and at the center of the composite in Fig. 5. The K effective that minimized the SSE does not perfectly replicate the transient response of composite as there is initially an underestimation and later an overestimation relative to the composite's temperature. The result is that the transient response of the composite materials cannot be precisely replicated using a single effective conductivity; the transient effects of the particles need to be taken into account. As the form of the matrix temperature rise does not follow exactly the analytical solution for an effective material, there may be errors introduced into the calculation of the effective conductivity. For example, the K effective that minimizes the SSE in Fig. 5 will depend on the range of time over which the SSE is calculated/minimized; using time ranges of 1 s, 0.25 s and 0.1 s leads to K effective values of 236, 241 and 249 W/(m•K), respectively.

Conclusions
A semi-analytical method has been proposed for evaluating transient heat flow in composite materials. These materials contain spherical inclusions of a lower thermal conductivity and size of which is much smaller than the macro-scale dimension of the composite. The heat exchange between the two phases is taken into account and the interfacial Eurasian Chemico-Technological Journal 21 (2019) 101-105 thermal resistance between the particles and matrix results in particle temperatures that lag behind those of the matrix. This leads to a smearing, or spreading out of the overall temperature evolution that cannot be accurately reproduced by treating the composite as a homogenous material having a single value of effective thermal conductivity. The example case discussed in Section 3 showed how this smearing can lead to fitted values of effective thermal conductivity that depend on the range of time used to obtain the fit. This may introduce a souce of error into characterization methods that measure only a part of the transient response.