Modelling Thermal Diffusivity of Heterogeneous Materials Based on Thermal Diffusivities of Components with Implications for Thermal Diffusivity and Thermal Conductivity Measurement

Effective thermal diffusivity models are useful for predicting thermal diffusivities of heterogeneous materials. The literature contains models that may be broadly categorised into four different types: (1) effective thermal diffusivity for highly specific applications (e.g. empirical curve fitting of measured data); (2) effective thermal diffusivity as a weighted averages of the components’ thermal diffusivities and volume fractions; (3) effective thermal diffusivity calculated from effective thermal conductivity, effective density and effective specific heat capacity known as the ‘lumped parameter’ approach (which is the most commonly employed method); (4) comparison of times for a fixed quantity of heat to be transferred to a composite material with the heat transfer time for a material with an effective thermal diffusivity. The latter three modelling methods were tested on theoretical composite materials, and none performed consistently better than the others, suggesting there is scope for further work in this area. Of the three methods, the least accurate on average was the lumped parameter method. Given that this relationship is often used to derive thermal conductivity data from thermal diffusivity data (or vice versa), it is possible that significant error is introduced to the derived property in addition to any measurement error, which is often not acknowledged.


Fo
Fourier number k Thermal conductivity (W·m −

Introduction
Modern technology often makes use of advanced materials that are composites or hybrids of two or more component materials [1]. The advantage of such materials is the ability to customise the material properties, within limits, according to a set of design criteria [2]. To this end, there is value in being able to model and predict properties of heterogeneous materials based on the properties of the component materials and some understanding of the structure of the composite, particularly in the initial stages of development when a range of specifications or design options need to be considered and evaluated.
The literature abounds with studies on the prediction of properties of heterogeneous or composite materials. Thermal conductivity prediction in particular has received much attention and numerous models have been developed [1,[3][4][5][6][7][8][9][10][11][12][13][14][15], and upper and lower bounds for model predictions have been established [2,[16][17][18]. By comparison, significantly less attention has been paid to predicting the effective thermal diffusivity of heterogeneous materials, perhaps because many researchers seem satisfied with deriving effective thermal diffusivity from effective thermal conductivity data (or vice versa) via the definition of thermal diffusivity, i.e.: where α is thermal diffusivity, k is thermal conductivity, ρ is density and c is specific heat capacity and the subscript e refers to the effective property (i.e. the property of the composite as a whole). However, despite making use of the formal definition of thermal diffusivity, it is claimed that this approach only provides approximations at best [19,20]. For composite materials comprised of layers of component materials, relatively simple formulae for effective thermal diffusivity prediction have been developed from first principles [20][21][22][23]; however, this is the only structure for which this is the case.
A number of researchers have produced empirical models to account for changes in effective thermal diffusivity with changes in composition and/or temperature (e.g. [24][25][26][27][28][29][30][31][32][33][34][35][36][37][38][39]). However, these types of models have limited value beyond data reduction or perhaps the identification of trends between variables. At the opposite end of the spectrum, there are highly elaborate models that may be used to obtain thermal diffusivity (and other properties), often by fitting temperature histories and profiles to measurement or simulation data (e.g. [40]); however, these types of modelling approaches rely on detailed structural information of the material being considered and are not suitable for quick calculations.
This paper describes a theoretical and numerical evaluation of different generic techniques (i.e. excluding purely empirical or highly material-specific models) for simple prediction of effective thermal diffusivity that are in current use; in particular the paper aims to test whether Eq. 1 (the Lumped Parameter method) is sufficient for use for every material, as appears to be assumed by many researchers. The implications from the evaluation exercise for the common practice of obtaining thermal diffusivity measurements for composites from thermal conductivity measurements (and vice versa) are also discussed.

Current Approaches to Modelling Effective Thermal Diffusivity
Of the types of modelling approaches described above, models that are simply empirical correlations of measured data, or are based on highly specific material structures will not be considered further in this paper, since they were never intended for generic application, unlike the ones outlined below.

Effective Thermal Diffusivity Models as Function of Composition and Components' Thermal Diffusivities: 'Averaging Method'
Two of the simplest, yet most widely employed thermal conductivity models are the Series and Parallel models [7], shown here for a 2-component material: where v is the volume fraction of the component. Inspection of Eqs. 2 and 3 reveals that they are simply the volume-weighted harmonic and arithmetic means, respectively, of the component materials' thermal conductivities. Similarly, the geometric mean of the components' thermal conductivities has been used to model effective thermal conductivities [7]. This simple approach of taking a weighted average of the components' physical properties has also been applied to thermal diffusivity, e.g. in the case of the Parallel model [26,41,42]: When applied to thermal conductivity, the Series and Parallel models do actually have theoretical bases, being derived from analysis of steady-state conduction resistances [43]. The Series model in particular provides exact thermal conductivity calculations for one-dimensional heat flow in layered materials (i.e. when the thermal resistances occur in series with respect to the temperature gradient); however, the same is not necessarily true when applied to thermal diffusivity.
Schimmel et al. [20] determined relationships for effective thermal diffusivity of multi-layered materials analytically for which one external surface was exposed to a heat flux and the other was insulated. For a two-component material their expression for effective thermal diffusivity was: Note than in Eq. 5, unlike Eqs. 2 to 4, the designation of component materials is not arbitrary and it is the external surface of Component 1 that experiences the temperature change, while the external surface of Component 2 is insulated. They compared predictions from their effective thermal diffusivity models as well as predictions using Eq. 1 (in which effective thermal conductivity was obtained using Eq. 2), to numerical simulations of temperature histories in the layered material. They concluded that while their effective thermal diffusivity models could not produce the same results as the numerical method (that they assumed provided the 'correct' solution), it was much closer than the predictions from Eq. 1. Similar analyses for layered materials have been performed allowing for the finite propagation time of the temperature front [21][22][23].
However, other than for layered materials there do not appear to be any effective thermal diffusivity models of this type that have been derived analytically, although a semi-empirical, modified form of an equation for layered materials derived by Mansanares et al. [21] has been proposed for use with materials having other structures [44].

Effective Thermal Diffusivities Determined from Effective Thermal Conductivity, Effective Density and Effective Specific Heat Capacities: 'Lumped Parameter Method'
The most common approach for predicting thermal diffusivity (sometimes referred to as the 'lumped parameter method' [20,38]), is to calculate it from effective thermal conductivity, effective density and effective specific heat capacities using Eq. 1 (e.g. [19,20,24,[45][46][47][48][49][50][51][52][53][54][55][56][57][58][59][60][61][62]). As mentioned in the Introduction, the number of different effective thermal conductivity models is large; however, the effective density and effective specific heat capacity are normally determined from the volume-weighted arithmetic mean according to Eq. 6 and the mass-weighted mean according to Eq. 7, respectively: Note that volume fractions and mass fractions are related by: Combining Eqs. 6 to 8: The rationale of this approach is that provided the correct model for effective thermal conductivity is employed Eq. 1 would be sufficiently accurate for predicting thermal diffusivity. However, it has been stated that Eq. 1 will only provide approximations for effective thermal diffusivity [19,20].
It is worth observing that the thermal diffusivities predicted using Eq. 1 are different from those predicted using the approach described in Sect. 2.1. For example, Eq. 10 shows an expanded version of Eq. 4 (i.e. the Parallel model applied to thermal diffusivity) while Eq. 11 is the application of Eq. 1 in which k e is modelled by the Parallel effective thermal conductivity model Eq. 3 and the effective density and specific heat capacity are modelled by Eqs. 6 and 7: It is clear by inspection that Eqs. 10 and 11 produce different results.

Ahmadi's Method
Ahmadi [19] defined effective thermal diffusivities not specifically in terms of the component materials' diffusivities or effective properties, but rather in terms of equivalent times for a transient heat transfer process to progress to a certain point. Using numerical methods, the time required for a heat conduction process in a representative micro-unit of a composite material to reach some fraction of completion defined in terms of total heat transferred was calculated, and compared to the analytical solution for the same process in a homogeneous material having an effective thermal diffusivity. For example, if the conduction process was to be some fraction β of completion, i.e.
then, according to his definition, the time required for the conduction process in the composite to reach β completion (t β , comp ) should be the same as the time required for an homogeneous material to reach β completion if the diffusivity of the homogeneous material is equal to the effective thermal diffusivity of the composite, i.e.: where the subscripts h and comp refer to the homogeneous and composite materials, respectively. His definition of α e is therefore: where t* is a dimensionless time variable (equivalent to the Fourier number multiplied by a factor of 4, according to the definition provided), and a is the characteristic dimension of the object.
Ahmadi claimed that no similar definition of effective thermal diffusivity exists in the open literature, and indeed no others were found as part of this study. However, it was acknowledged that the α e calculated using Eq. 10 has a slight dependence on the value of β, which is undesirable. What was not discussed was whether t and t* also depend on the boundary conditions imposed, i.e. whether if α e is determined from one set of boundary conditions for a given value of β it will necessarily be accurate if that α e value is applied to the same material under a different set of boundary conditions. In addition, the fact that it must be determined numerically makes this method less convenient than those outlined in Sects. 2.1 and 2.2.

Evaluation of Different Methods for Determining α e
In order to evaluate the different prediction methods a simple composite system is employed for which analytical solutions to transient heat transfer are available, following a similar approach to that employed by Schimmel et al. [20]. Consider a (12) composite material made up of two slabs with dissimilar physical properties where heat flows in one-dimension through the two-component materials in series (Fig. 1). For simplicity, it is assumed that contact resistance is negligible. The material as shown in Fig. 1 has top and bottom boundaries insulated (i.e. dT/dy = 0 for all x), while the left-hand boundary is maintained at T A and the right-hand boundary is at T B . The material is initially at uniform temperature T 0 , and then for t > 0 T A is suddenly increased to some new value, while T B is maintained at T B = T 0 . The origin of the x-coordinate is the interface between Components 1 and 2.
Within each material, the transient conduction is governed by the Diffusion equation [43]: The conditions at the interface between the two slabs are: Eventually a steady temperature distribution is attained across the material with a steady heat flux between T A and T B , in which case: Between T A and T B the steady heat flux is: The purpose of an effective property is to be able to account for the heat transfer behaviour with a single property rather than having to account for the actual spatial variation in the property, so that Eq. 15 becomes: Observing that in this case a/X = v 1 and b/X = v 2 , the heat flow evaluated by Eq. 19 will be the same as that calculated by Eq. 20 if k e is calculated according to Eq. 2; in other words, Eq. 2 may be derived by a comparison of Eqs. 19 and 20. It is important to observe that Eq. 19 applies regardless of whether T A or T B experiences the step change, with the only difference being the sign of the heat flux (i.e. positive or negative).
The solution of Eq. 15 with these boundary conditions for the transient heat transfer process is provided by Carslaw and Jaegar [63], with T B = T 0 = 0: where the λ n are the roots of: and: sin 2 b n sin n (a + x) n (asin 2 b n + bsin 2 a n ) sina n sin b n sin (b − x) n n (asin 2 b n + bsin 2 a n )    (Table 1, [64]) with a = b = 0.05 m, and X = Y = 0.1 m. The transition from the initial uniform temperature distribution to the steady-state distribution once the infinite series term on the right-hand side of Eqs. 21 and 22 has become negligible is clearly revealed. Figure 2b shows the temperature profiles when the positions of the component materials relative to the boundary that experiences the temperature change are reversed. It is clear that the temperature profile both during the transient process and under steady-state conditions depends on which component the heat flows through first [20,23,65].
The solution of Eq. 15 with the same external boundary conditions but for a homogeneous material is given by Carslaw and Jaegar [63] as: where in this case the coordinate origin is at the left-hand boundary (surface A). Unlike the steady-state problem for thermal conductivity where the relationship between k e and k 1 , k 2 was obvious, it is not a simple exercise to determine the proper relationship for α e with respect to α 1 and α 2 by comparing Eq. 26 with Eqs. 21 and 22). This demonstrates that attempting to obtain an effective thermal diffusivity is inherently a more challenging task than determining effective thermal conductivity from a steady-state analysis. Figure 3 shows the temperature history at the interface between the two-component materials (recalling that perfect thermal contact is assumed) according to the analytical solutions Eqs. 21 and 22) both with water as Component 1 (i.e. the component lying on the boundary which experiences the temperature change) and vice versa. Figure 3 also includes temperature histories plotted using Eq. 22 and the effective thermal diffusivities determined by the different methods outlined in Sect. 2, as follows. The first of these methods ('averaging method', Sect. 2.1) was to employ some volume-weighted average of the components' diffusivities. In this case the Series model (harmonic mean) seems logical since the components are in series with respect to the heat flow: The second method (Lumped Parameter method, Sect. 2.2) was to employ Eq. 1 with effective thermal conductivity modelled by the Series model and effective density and effective specific heat capacity modelled as per Eqs. 6 and 7: Ahmadi's method (Sect. 2.3) was implemented using finite element analysis based on a heat transfer completion fraction (β) of 0.9. From Eq. 14, when Component 1 had the properties of water the predicted effective thermal diffusivity was 1.92 × 10 -7 m 2 ·s −1 ; when Component 1 had the properties of sandstone it was 4.15 × 10 -7 m 2 ·s −1 . Note that in order to implement Ahmadi's method, the  Fig. 1 insulated rather than held at fixed temperature). Figure 3 shows that none of the three methods is able to model the exact solution for the temperature histories at the interface between the component materials, although Ahmadi's method is able to account for the fact that the temperature history at the interface is affected by which of the component materials the heat flows through first.
The material shown in Fig. 1 is highly non-homogeneous, so it is worth considering whether the effective thermal diffusivity prediction methods would be more accurate for a composite made up of a larger number of alternating layers. To this end numerical simulations using the Finite Element method were performed, since the analytical solutions become increasingly complex with the addition of extra layers [20,23,63]. The numerical methods were first checked against the analytical solutions for the 2-layer composite in Fig. 1. Simulations were performed for a composite having the same overall dimensions and   (Table 1). In this situation the positions of the components with respect to the boundary that experiences the temperature change has only a minor influence and the two temperature histories almost superimpose each other, so they are not displayed graphically. Table 2 lists the effective thermal diffusivities calculated by the three different methods applied to two different pairings of component materials: sandstone/water and sandstone/aluminium. Note that the averaging method and lumped parameter method effective thermal diffusivities are independent of the number of layers and so they are unchanged from the 2-layer material, while those obtained using Ahmadi's method are different. The mean square error (Eq. 29) between the temperature histories from the numerical simulations and those calculated using Eq. 26 with an effective thermal diffusivity during the transient stage are also shown in Table 2.
All three prediction methods perform much better for the 40-layer material than for the 2-layer material, and both the averaging method and Ahmadi's method provide reasonably accurate predictions. The lumped parameter method does not perform nearly as well as the other two methods.
It is worth considering a different geometry, such as the one illustrated in Fig. 4. This structure approximates the matrix + filler structure that is common for composite materials. The effective thermal conductivity of this structure is modelled well by the 2-dimensional form of the Maxwell-Eucken model (Eq. 30) [5,66], particularly for low volume fractions of the dispersed phase (i.e. the circles): where the subscripts c and d refer to the continuous and dispersed phases, respectively. Therefore, the lumped parameter effective thermal diffusivity is: In addition to Eqs. 4 and 27, it is also worth considering whether the Maxwell's model applied to thermal diffusivity provides reasonable accuracy using the 'averaging' approach, i.e.: Table 3 shows a comparison between the numerical results and the effective thermal diffusivity method for a system in which X = Y = 0.1 m, the radius of the Lumped parameter (Eq. 28  Of the three test cases considered (2-layered material, 40-layered material, matrix + filler material) none of the three modelling approaches consistently provided the most accurate predictions, although, if an MSE value of less than 0.01 can be considered 'sufficient accuracy', it is clear that sufficiently accurate predictions were obtained for the latter two cases by more than one method. However, it is worth observing that the lumped parameter approach, which appears to be employed most often for effective thermal diffusivity in the literature was the least accurate method on average.
The averaging method is simple, and has the potential for greater accuracy than the lumped parameter method; however, it may be difficult to determine the appropriate mathematical relationship between the components' thermal diffusivities, since it clearly will not always be either the arithmetic or harmonic mean. While it was logical to identify the harmonic mean (Eq. 27, Series model) for the layered materials, and it provided adequate predictions for the matrix + filler material when sandstone was the continuous phase, it was not suitable for the matrix-filler material when aluminium was the continuous phase. Simply 'translating' an effective thermal conductivity model to effective thermal diffusivity is not the solution, since Eq. 27 provided better predictions than Eq. 32 for the sandstone matrix/aluminium filler material.
Ahmadi's method [19] accounts for the transient nature of the problem by defining effective thermal diffusivity in terms of time; however, it is more cumbersome to implement than the lumped parameter or averaging methods, and for some materials it may be difficult to identify a representative sub-unit that is repeated throughout the volume of the material. While the method can provide very accurate effective diffusivities if they are applied under the same boundary conditions for which they were determined, they are not guaranteed to be more accurate than the averaging method if they are applied under different boundary conditions. It appears, then, that there is scope for more effort to be devoted to improve effective thermal diffusivity modelling. In particular it would be valuable if thermal diffusivity bounds analogous to those that have proven to be so useful in effective thermal conductivity analysis [16][17][18] could be determined, although this is likely to be a more challenging exercise than for thermal conductivity.

Implications for Thermal Conductivity and Thermal Diffusivity Measurement
It is worth analysing why the lumped parameter approach does not perform well compared to the other two approaches. Consider again the 2-layer material depicted in Fig. 1 with Component 1 having the properties of water. Figure 2a shows that during the early stages of the transport process the temperature gradient exists entirely within Component 1 and the effective properties of the composite are the same as the actual properties of Component 1. As the temperature 'front' crosses the interface between the two components, the effective properties now are dependent on the properties of Component 2 in addition to the properties of Component 1. However, shortly after the temperature front has crossed the interface, the temperature gradient is still mainly within Component 1 and so the effective properties are affected more by the properties of Component 1 than the properties of Component 2. For example, if the front has reached x = 0.03 m (Fig. 1) the effective volume fraction of Component 1 is 0.5/0.8 = 0.625, while that of Component 2 is 0.3/0.8 = 0.375, and the properties of Component 1 still influence the effective diffusivity significantly more than Component 2. It is only once the steady-state gradient has been fully established that the effective diffusivity is affected by the thermal properties of the two components in equal measure since v 1 = v 2 .
In order for the accuracy of Eq. 1 to be improved, the time-averaged values of effective thermal conductivity, specific heat capacity and density, rather than the steady-state values, should be used to estimate effective thermal diffusivity: However, usually it is the steady-state values that are employed, which causes α e either to be underestimated (e.g. in this case when Component 1 has the properties of sandstone), or overestimated (if Component 1 has the properties of water).
A large number of thermal conductivity measurements are performed under steady-state or quasi-steady-state conditions, and often there is subsequent calculation of thermal diffusivity within the same study. Alternatively, thermal conductivity data are presented that have been calculated from thermal diffusivity data that were measured using a transient method. Consider the 2-layer material in Fig. 1 being subjected to a transient heat transfer measurement (with boundary B insulated rather than held at fixed temperature) where the effective thermal diffusivity is fitted to the temperature history at the insulated boundary. If Component 1 has the properties of sandstone, and Component 2 has the properties of water, the effective thermal diffusivity is 4.18 × 10 -7 m 2 ·s −1 (from fitting the analytical solution for a homogeneous material to the numerical simulation history for the heterogeneous material). The effective thermal conductivity calculated from this thermal diffusivity using Eq. 1 would be 1.20 W·m −1 ·K −1 . However, if the same material was subjected to a steadystate measurement the thermal conductivity is 0.90 W·m −1 ·K −1 , and the thermal diffusivity calculated from it using Eq. 1 would be 3.15 × 10 -7 m 2 ·s −1 . For both thermal diffusivity and thermal conductivity there is a difference of 33 % between the property measured directly and the one derived using Eq. 1.
This example involving the 2-layered material is a case of extreme inhomogeneity, and the error would not be as significant in the case of the 40-layered material. However, Eq. 1 is commonly used to derive thermal conductivity data from thermal diffusivity data (or vice versa) for heterogeneous materials without any recognition or consideration of the potential for significant error being introduced by its use (for example, thermal conductivity data derived from thermal diffusivity data can be very different from directly measured thermal conductivity [67]). Given that the definitions of both thermal diffusivity and thermal conductivity are in macroscopic terms, defining the scale at which a multi-component material becomes 'homogeneous' is (33)  . a difficult task. However, if distinct 'phases' in a material can be identified by the naked eye then it is likely that error will be introduced if Eq. 1 is employed to derive thermal conductivity from measured thermal diffusivity data or vice versa. Clearly it is preferable for each property to be measured directly if possible.
It would be valuable to perform a study where direct measurements of thermal diffusivity and thermal conductivity are made on a range of different composite materials with varying degrees of inhomogeneity, and then to compare these to conductivities and diffusivities derived from data for the other property. In this way it may be possible to quantify the extent of any uncertainty introduced when data for either property is derived from measurements of the other.

Conclusions
Three categories of effective thermal diffusivity modelling methods (averaging, lumped parameter, and Ahmadi) were tested on three theoretical composite materials, and none performed consistently better than the others. Of the three methods, the least accurate on average was the lumped parameter method. Given that this relationship is often used to derive thermal conductivity data from thermal diffusivity data (or vice versa), it is probable that significant error is introduced to the derived property in addition to any measurement error, which is often not acknowledged.
Funding Open Access funding enabled and organized by CAUL and its Member Institutions.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.