Generalized Effective Medium Theory for Particulate Nanocomposite Materials

The thermal conductivity of particulate nanocomposites is strongly dependent on the size, shape, orientation and dispersion uniformity of the inclusions. To correctly estimate the effective thermal conductivity of the nanocomposite, all these factors should be included in the prediction model. In this paper, the formulation of a generalized effective medium theory for the determination of the effective thermal conductivity of particulate nanocomposites with multiple inclusions is presented. The formulated methodology takes into account all the factors mentioned above and can be used to model nanocomposites with multiple inclusions that are randomly oriented or aligned in a particular direction. The effect of inclusion dispersion non-uniformity is modeled using a two-scale approach. The applications of the formulated effective medium theory are demonstrated using previously published experimental and numerical results for several particulate nanocomposites.


Introduction
Typically, the process of heat conduction is treated using the classical Fourier's law. Although Fourier's law is widely applied, its application to systems with characteristic lengths comparable to or lower than the mean-free-path of the energy carriers (phonons or electrons) leads to large errors in any or all variables in the system such as the thermal conductivity, temperature and the temperature gradient [1]. Examples of such systems include nanoparticles, nanowires and thin films. The reason for the inapplicability of Fourier's law for nanostructures is that during heat conduction in such systems, equilibrium conditions are not achieved and therefore, a local temperature gradient is not established. Without a temperature gradient, Fourier's law is not applicable and thermal conductivity, which relates heat flux to the temperature gradient, has no meaning. However, if one is interested in heat conduction in a domain with a characteristic length much greater than the characteristic length of the nanostructures, Fourier's law may still be applied using an effective value of thermal conductivity. An example of such a case is a nanocomposite where the domain under consideration is considerably larger than the characteristic length of the nanoparticles.
The problem of estimation of the effective thermal conductivity of composite materials has been widely studied. Early works in the area were done by Maxwell [2] and Lord Rayleigh [3] who studied the thermal conductivities of composites with low concentrations of inclusions. Their works were later extended by Hasselman and Johnson [4] and Benveniste [5] who studied the effects of thermal boundary conductance on the effective thermal conductivity of the composite. Bruggeman [6] derived a model for the effective thermal conductivity of a composite when the inclusion concentration was high. His model was later extended by Every and coworkers [7] to include the effect of thermal boundary conductance. Modifications to the model by Every and coworkers to include the effect of particle include the effect of particle shape have also been presented [8]. An effective medium theory (EMT) for the estimation of thermal conductivity of composites with dilute concentrations of inclusions of different shapes was presented by Nan and coworkers [9]. The major drawbacks of the effective medium theory approach discussed above include the limitation of using inclusions of regular shapes, inability to handle non-uniformly dispersed and percolating inclusions and inability to deal with nanocomposites [1]. To overcome the inability to handle percolating networks of inclusions, Prasher and coworkers [10,11] and Evans and coworkers [12] presented a three-level homogenization methodology capable of handling clustering of inclusions.
The application of the models mentioned above to nanocomposites can lead to significant errors in the predicted effective thermal conductivity. Many times, the addition of nanoparticles in a matrix can result in an effective thermal conductivity which is significantly lower than the thermal conductivities of both the matrix and the inclusion [13,14]. The use of conventional modeling approaches for such composites results in an overestimation of the effective thermal conductivity even when using very high interface thermal resistance [15,16]. Wu and coworkers [15] used Nan and coworkers' model [17] to estimate the thermal conductivity of aluminum matrix composites with carbon nanotube (CNT) inclusions and found that the model over-predicted the results even when the thermal interface resistance was set to infinity. Ahmad and coworkers [16] carried out a similar study for alumina-CNT composite and arrived at the same conclusion.
To reduce the errors in effective medium theory predictions for nanocomposites, a modified effective medium theory was presented by Minnich and Chen [18] for spherical inclusions and was extended by Ordonez-Miranda and coworkers [19] for spheroidal inclusions. In their approach, modified thermal conductivities of the matrix and inclusions are first calculated and then used in the effective medium theory. Minnich and Chen used the modified values of matrix and inclusion thermal conductivities in Nan and coworkers' EMT for spherical inclusions and found good agreement between the effective thermal conductivities predicted by the modified EMT and Monte Carlo simulations. The shortcomings of the modified effective medium theory approach of Minnich and Chen and Ordonez-Miranda and coworkers include the inability to handle multiple inclusions, randomly oriented inclusions, and non-uniformly distributed inclusions. Other approaches used for the estimation of effective thermal conductivity of nanocomposites include the Monte Carlo simulation method [14,20,21], molecular dynamics [22] and the Boltzmann transport equation [23][24][25][26][27].
The current work shows the formulation of a generalized effective medium theory which can be used to determine the effective thermal conductivity of particulate nanocomposites. The proposed effective medium theory overcomes several shortcomings in the effective medium theory approaches reported in the literature. These include capabilities to include the effect of multiple nanometer-sized inclusions, the effect of oriented (randomly or at any angle relative to heat flow direction) spheroidal, cylindrical or platelet inclusions and the effect of non-uniformly dispersed inclusions. Figure 1 graphically shows the capability of the formulated generalized effective medium theory in comparison to those proposed by Minnich and Chen [18] and Ordonez-Miranda and coworkers [19]. The formulation of the generalized effective medium theory is presented in Section 2 while Section 3 shows the application of the formulated EMT.   [18]; (b) Ordonez-Miranda and coworkers [19]; and (c) generalized EMT-current work; (d) A schematic representation nonuniformly distributed inclusions-current work (inclusion sizes are not to scale).

Framework of the Generalized Effective Medium Theory
The formulation of the generalized effective medium theory is presented in this section. First, the effective medium theory formulation for multiple inclusions is presented in Section 2.1. Second, the effect of multiple nanometer-sized inclusions on the matrix and inclusion thermal conductivities is derived in Section 2.2. Lastly, a two-scale approach to handle non-uniformly dispersed inclusions is presented in Section 2.3.

Effective Medium Theory for Composites with Multiple Inclusions
To derive the effective medium theory for composites with multiple inclusions, a two-phase composite was first considered. The thermal conductivity of the composite was assumed to vary from point to point according to the function thermal conductivity function and   n K r is the fluctuation in thermal conductivity due to inclusion particle n. The effective thermal conductivity of the composite was then defined by Equation (1) [28].
where I is the unit tensor, G is the Green's tensor [29] and denotes volumetric average.
The tensor T, defined by Equation (2) is the transition matrix and it describes the effect of individual inclusion particles on the effective thermal conductivity of the composite.
If the inclusion volume fraction is low, K 0 can be taken as the matrix thermal conductivity Kmat, the fluctuating part becomes Kinc − Kmat and T can be approximated as, Using this approach, Nan and coworkers [9] formulated the effective medium theory of composites with spheroidal inclusions of one type. In the current work, their approach was used to extend the methodology for multiple inclusions. First, Equation (3) was extended to Equation (4) for multiple types of inclusions. The effective thermal conductivity of the composite was then calculated using Equations (5)- (13) Figure 1. Inclusion geometries and orientation for: (a) Minnich and Chen [18]; (b) Ordonez-Miranda and coworkers [19]; and (c) generalized EMT-current work; (d) A schematic representation non-uniformly distributed inclusions-current work (inclusion sizes are not to scale).

Framework of the Generalized Effective Medium Theory
The formulation of the generalized effective medium theory is presented in this section. First, the effective medium theory formulation for multiple inclusions is presented in Section 2.1. Second, the effect of multiple nanometer-sized inclusions on the matrix and inclusion thermal conductivities is derived in Section 2.2. Lastly, a two-scale approach to handle non-uniformly dispersed inclusions is presented in Section 2.3.

Effective Medium Theory for Composites with Multiple Inclusions
To derive the effective medium theory for composites with multiple inclusions, a two-phase composite was first considered. The thermal conductivity of the composite was assumed to vary from point to point according to the function K (r) = K 0 + ∑ n δK n (r) where K 0 is the constant part of thermal conductivity function and δK n (r) is the fluctuation in thermal conductivity due to inclusion particle n. The effective thermal conductivity of the composite was then defined by Equation (1) [28].
where I is the unit tensor, G is the Green's tensor [29] and denotes volumetric average. The tensor T, defined by Equation (2) is the transition matrix and it describes the effect of individual inclusion particles on the effective thermal conductivity of the composite.
If the inclusion volume fraction is low, K 0 can be taken as the matrix thermal conductivity K mat , the fluctuating part becomes K inc − K mat and T can be approximated as, Using this approach, Nan and coworkers [9] formulated the effective medium theory of composites with spheroidal inclusions of one type. In the current work, their approach was used to extend the methodology for multiple inclusions. First, Equation (3) was extended to Equation (4) for multiple types of inclusions. The effective thermal conductivity of the composite was then calculated using Equations (5)-(13) which are valid for spheroidal inclusions of multiple types.
K e f f , 11 = K e f f ,22 = K mat where ϕ i is the volume fraction, a i 1 and a i 3 are the particle radii, p i is the aspect ratio, R i TB is the interfacial thermal resistance, K inc.i is the thermal conductivity of inclusion of type i, K mat is the thermal conductivity of the matrix, and cos 2 θ i is a factor defining the orientation of inclusion of type i.

Effect of Nanometer-Sized Inclusions on Matrix and Inclusion Thermal Conductivities
The effect of multiple nanometer-sized inclusions oriented in any random direction on the thermal conductivities of the matrix and inclusions was calculated by extending the approach of Minnich and Chen [18]. According to their approach, the addition of the nanometer-sized particles modifies the thermal conductivities of the matrix and the inclusions. The modified thermal conductivities are calculated by first calculating the effective mean-free-path of the energy carriers (phonons or electrons) in the matrix or the inclusions using Matthiessen's rule.
where x can be p or e for phonons or electrons respectively and y can be mat for the matrix or inc.i for the ith inclusion.
The bulk mean-free-paths, Λ y x,bulk for a material are known quantities. The unknown quantity in Equation (14) is the collision mean-free-path, Λ y coll , defined as the average distance traveled by the energy carriers between collisions. Once the effective mean-free-paths of the energy carriers have been calculated, the modified thermal conductivities of the matrix and inclusions can be calculated using, x,e f f (15) To calculate the collision mean-free-paths of the matrix, Λ mat coll , the densities of the inclusions in the matrix, n i , were calculated using Equation (16).
where N is the number of inclusions in a sample of volume V and V i is the volume of a single particle of inclusion type i.
To determine the collision mean-free-path of the matrix, it was noted that the number of inclusions of type i that an energy carrier (phonon or electron) encounters in a volume A ⊥1 L is n i A ⊥1 L assuming is the collision cross-section area of the ith inclusion type. The matrix collision mean-free-path is therefore, By replacing the expressions for n i in Equation (17), we get, where Equation (18) reduces to Ordonez-Miranda and coworkers' formulation [19] when only a single inclusion was considered; that is when ϕ i = 0 for i = 2, 3, . . . , N. It is also interesting to note that in Equation (18), the relative sizes of the inclusions, represented by r Vi , determines their contribution to the collision mean-free-path for the matrix.
Applying Matthiessen's rule, the effective mean-free-path of energy carriers in the matrix was estimated as, Λ mat x,e f f = Λ mat x,bulk The modified thermal conductivity of the matrix material was calculated using Equation (15), which led to Equation (20).
The total matrix thermal conductivity is given by Equation (21).
For the calculation of the modified matrix thermal conductivity, the collision cross-sectional area of the inclusion is required. For the general case of a spheroidal inclusion oriented at any angle, Equation (22) was derived for the calculation of the collision cross-sectional area. The detailed derivation for Equation (22) is presented in Appendix A.
a 1 1 sin 2 θ 1 + a 1 3 cos 2 θ 1 (23) and θ 1 is the average orientation angle of inclusion 1. The dependence of the spheroidal cross section area A ⊥1 on θ 1 is shown in Figure 2. Equation (22) was derived for the calculation of the collision cross-sectional area. The detailed derivation for Equation (22) is presented in appendix A.
where,  is the average orientation angle of inclusion 1.
The dependence of the spheroidal cross section area A  on 1  is shown in Figure 2. For the inclusions, the collision mean-free-path depends only on the size and orientation of the inclusion itself. For the special cases of inclusion aligned parallel ( 0 to the direction of heat conduction, the collision mean-free-path for inclusion particles can be calculated using Equations (25) or (26), respectively [19]. For the general case of particles oriented at any angle between 0° and 90°, Equation (24) were derived to calculate the collision mean-free-path of inclusion particles. Figure 3 shows the reasoning behind Equation (24). For randomly oriented inclusions, the average orientation of the inclusions was estimated from Figure 4 for spheroids of aspect ratios 0.1, 1 and 10.  For the inclusions, the collision mean-free-path depends only on the size and orientation of the inclusion itself. For the special cases of inclusion aligned parallel (θ i = 0 • ) or perpendicular (θ i = 90 • ) to the direction of heat conduction, the collision mean-free-path for inclusion particles can be calculated using Equations (25) or (26), respectively [19]. For the general case of particles oriented at any angle between 0 • and 90 • , Equation (24) were derived to calculate the collision mean-free-path of inclusion particles. Figure 3 shows the reasoning behind Equation (24). For randomly oriented inclusions, the average orientation of the inclusions was estimated from cos 2 The variation of Λ inc.i coll with angle θ i is shown in Figure 4 for spheroids of aspect ratios 0.1, 1 and 10.
where ε is the eccentricity given by, and K ellip ( ) and E ellip ( ) are elliptic integrals of first and second kind.
, for p 1 3 1 arcsinh where  is the eccentricity given by, ,for 1   The effective mean-free-path of the energy carriers in the inclusion particles was calculated by applying the Matthiessen's rule.
, for p 1 3 1 arcsinh where  is the eccentricity given by, ,for 1   The effective mean-free-path of the energy carriers in the inclusion particles was calculated by applying the Matthiessen's rule.  The effective mean-free-path of the energy carriers in the inclusion particles was calculated by applying the Matthiessen's rule.
x,e f f = By replacing the effective mean-free-path in Equations (15) for thermal conductivity, the reduced electron or phonon thermal conductivities of the inclusions were derived using Equation (29).
The total thermal conductivity for inclusion i is the sum of the phonon and electron thermal conductivities given by Equation (30).

Two-Scale Approach for Non-Uniformly Distributed Inclusions
A major drawback of the effective medium theory approach is its inability to handle non-uniform dispersion of inclusions. In the current work, a two-scale approach was proposed that combines the EMT approach at the lower scale with computational homogenization approach at the upper scale [30,31] to introduce the effect of non-uniform dispersion of inclusions in the matrix on the effective thermal conductivity of the composite.
The application of the method requires a quantitative knowledge of inclusion distribution in the matrix. This could be a statistical parameter like standard deviation of inclusion volume fraction inside the matrix that can be used to develop a representative volume element (RVE) for computational homogenization. The distribution could also be experimentally determined using material characterization techniques. The key point in the application of the two-scale methodology is the requirement of scale separation between the size scale of inclusions in the composites and the scale at which the homogenized properties are being estimated. This means that the distribution needs to be measured at a length scale larger than the inclusion size. As an example, if the inclusions are nanometer-sized, the distribution needs to be determined at the micrometer length scale. The two-scale methodology is presented in Figure 5.
By replacing the effective mean-free-path in Equations (15) for thermal conductivity, the reduced electron or phonon thermal conductivities of the inclusions were derived using Equation (29).
The total thermal conductivity for inclusion i is the sum of the phonon and electron thermal conductivities given by Equation (30).

Two-Scale Approach for Non-Uniformly Distributed Inclusions
A major drawback of the effective medium theory approach is its inability to handle nonuniform dispersion of inclusions. In the current work, a two-scale approach was proposed that combines the EMT approach at the lower scale with computational homogenization approach at the upper scale [30,31] to introduce the effect of non-uniform dispersion of inclusions in the matrix on the effective thermal conductivity of the composite.
The application of the method requires a quantitative knowledge of inclusion distribution in the matrix. This could be a statistical parameter like standard deviation of inclusion volume fraction inside the matrix that can be used to develop a representative volume element (RVE) for computational homogenization. The distribution could also be experimentally determined using material characterization techniques. The key point in the application of the two-scale methodology is the requirement of scale separation between the size scale of inclusions in the composites and the scale at which the homogenized properties are being estimated. This means that the distribution needs to be measured at a length scale larger than the inclusion size. As an example, if the inclusions are nanometer-sized, the distribution needs to be determined at the micrometer length scale. The two-scale methodology is presented in Figure 5. In the first two steps in the methodology, the variation of inclusion volume fraction at different points in the composite is determined and represented in the form of a grid of points for each inclusion. As a result, the distribution of inclusions in the composite can be shown as in Figure 6, which shows the distribution of a nanometer sized Si inclusion in Ge matrix. Because of the scale In the first two steps in the methodology, the variation of inclusion volume fraction at different points in the composite is determined and represented in the form of a grid of points for each inclusion. As a result, the distribution of inclusions in the composite can be shown as in Figure 6, which shows the distribution of a nanometer sized Si inclusion in Ge matrix. Because of the scale separation requirement, each point in the grid represents the composite material having a certain volume fraction of inclusions. In the example considered in the figure, the average inclusion volume fraction in the Ge-Si RVE was 5% and the standard deviation of the volume fraction in the microscale domain was 1%. The distribution of the effective thermal conductivity in the RVE determined using the generalized effective medium theory is shown in Figure 7.
Materials 2016, 9, 694 9 of 21 separation requirement, each point in the grid represents the composite material having a certain volume fraction of inclusions. In the example considered in the figure, the average inclusion volume fraction in the Ge-Si RVE was 5% and the standard deviation of the volume fraction in the microscale domain was 1%. The distribution of the effective thermal conductivity in the RVE determined using the generalized effective medium theory is shown in Figure 7.  The final step in the methodology is to apply computational homogenization to estimate the overall effective thermal conductivity of the composite. An important aspect of the accuracy of computational homogenization results is the correct determination of the RVE size required for homogenization. In the current work, the validity of the RVE size used was determined using the methodology presented by Gitman and coworkers [32]. Gitman and coworkers suggested the calculation of a variation coefficient, known as the chi-square criterion (χ 2 ) using Equation (31). If the chi-square criterion is below a threshold value, the homogenization results are considered to be independent of RVE size. A threshold of 0.1 was selected in the present study as suggested by Gitman and coworkers    separation requirement, each point in the grid represents the composite material having a certain volume fraction of inclusions. In the example considered in the figure, the average inclusion volume fraction in the Ge-Si RVE was 5% and the standard deviation of the volume fraction in the microscale domain was 1%. The distribution of the effective thermal conductivity in the RVE determined using the generalized effective medium theory is shown in Figure 7.  The final step in the methodology is to apply computational homogenization to estimate the overall effective thermal conductivity of the composite. An important aspect of the accuracy of computational homogenization results is the correct determination of the RVE size required for homogenization. In the current work, the validity of the RVE size used was determined using the methodology presented by Gitman and coworkers [32]. Gitman and coworkers suggested the calculation of a variation coefficient, known as the chi-square criterion (χ 2 ) using Equation (31). If the chi-square criterion is below a threshold value, the homogenization results are considered to be independent of RVE size. A threshold of 0.1 was selected in the present study as suggested by Gitman and coworkers    The final step in the methodology is to apply computational homogenization to estimate the overall effective thermal conductivity of the composite. An important aspect of the accuracy of computational homogenization results is the correct determination of the RVE size required for homogenization. In the current work, the validity of the RVE size used was determined using the methodology presented by Gitman and coworkers [32]. Gitman and coworkers suggested the calculation of a variation coefficient, known as the chi-square criterion (χ 2 ) using Equation (31). If the chi-square criterion is below a threshold value, the homogenization results are considered to be independent of RVE size. A threshold of 0.1 was selected in the present study as suggested by Gitman and coworkers where a i is the homogenized property under consideration for RVE realization i, n is the total number of random RVE realizations of the same size and a is the mean value of the property under consideration. For the computational homogenization problem of the current work, the parameter χ 2 was determined for multiple average inclusion volume fractions and multiple dispersion non-uniformity values (modeled using standard deviation of inclusion volume fraction). The results of this study are presented in Figure 8. From the results, it was concluded that the use of RVEs of 10 µm edge length will result in homogenized thermal conductivity values which are independent of the RVE size. Therefore, this RVE size was used.
values (modeled using standard deviation of inclusion volume fraction). The results of this study are presented in Figure 8. From the results, it was concluded that the use of RVEs of 10 μm edge length will result in homogenized thermal conductivity values which are independent of the RVE size. Therefore, this RVE size was used.
Computational homogenization was carried using COMSOL/MATLAB. A mesh convergence study was conducted to determine the number of elements required for results to be independent of the element size. It was found that increasing the number of elements beyond 50 × 50 did not result in any significant change in the results. The final mesh having 50 × 50 elements used in the current work to carry out computational homogenization is shown in Figure 9.
The effect of distribution non-uniformity was incorporated into the mesh using an interpolation function which applied the thermal conductivity distribution, as shown in Figure 7, to the integration points in the finite element mesh. It is important to note here that due to the condition of scale separation, the inclusion particles were not explicitly modeled in the geometry of the finite element domain. A temperature gradient was applied along the x-and y-directions of the RVE and Equation (32) was used to determine the effective thermal conductivity of the composite.
where q is the heat flux, T  is the temperature gradient, and  Computational homogenization was carried using COMSOL/MATLAB. A mesh convergence study was conducted to determine the number of elements required for results to be independent of the element size. It was found that increasing the number of elements beyond 50 × 50 did not result in any significant change in the results. The final mesh having 50 × 50 elements used in the current work to carry out computational homogenization is shown in Figure 9. For the example of 5% Ge-Si composite considered above, the overall effective thermal conductivity of the composite was determined to be 21.41 W/m·K. For the case when inclusion distribution was uniform, the effective thermal conductivity was determined to be 21.06 W/m·K.

Applications of the Generalized Effective Medium Theory
In this section, the generalized effective medium theory formulated in the previous section was applied to different particulate composites to show its various capabilities. Before the studies were carried out, the model was validated against two experimental datasets for Al2O3-SiC platelet The effect of distribution non-uniformity was incorporated into the mesh using an interpolation function which applied the thermal conductivity distribution, as shown in Figure 7, to the integration points in the finite element mesh. It is important to note here that due to the condition of scale separation, the inclusion particles were not explicitly modeled in the geometry of the finite element domain. A temperature gradient was applied along the x-and y-directions of the RVE and Equation (32) was used to determine the effective thermal conductivity of the composite.
K e f f = q 11 ∇T 11 + q 22 ∇T 22 + q 33 ∇T 33 where q is the heat flux, ∇T is the temperature gradient, and • = 1 |V| V •dV and V is the domain volume.
For the example of 5% Ge-Si composite considered above, the overall effective thermal conductivity of the composite was determined to be 21.41 W/m·K. For the case when inclusion distribution was uniform, the effective thermal conductivity was determined to be 21.06 W/m·K.

Applications of the Generalized Effective Medium Theory
In this section, the generalized effective medium theory formulated in the previous section was applied to different particulate composites to show its various capabilities. Before the studies were carried out, the model was validated against two experimental datasets for Al 2 O 3 -SiC platelet composite [33] and SiO 2 -CNT composite [34]. The results of these validations are presented in Figures 10 and 11 for Al 2 O 3 -SiC composite and SiO 2 -CNT composite, respectively. For both validations, the matrix thermal conductivity was set as the experimentally determined value and the inclusion dimensions were taken from the original experimental data. The model predictions were calculated by varying the thermal interface resistance in the range of 1 × 10 −8 to 8 × 10 −8 m 2 ·K/W. This was done to take into account the variation in thermal interface resistance that can occur due to change in process parameters, the source of raw materials and can even vary greatly from sample to sample [35]. Model predictions for both composites showed good agreement with experimental measurements. For the example of 5% Ge-Si composite considered above, the overall effective thermal conductivity of the composite was determined to be 21.41 W/m·K. For the case when inclusion distribution was uniform, the effective thermal conductivity was determined to be 21.06 W/m·K.

Applications of the Generalized Effective Medium Theory
In this section, the generalized effective medium theory formulated in the previous section was applied to different particulate composites to show its various capabilities. Before the studies were carried out, the model was validated against two experimental datasets for Al2O3-SiC platelet composite [33] and SiO2-CNT composite [34]. The results of these validations are presented in Figure  10 and Figure 11 for Al2O3-SiC composite and SiO2-CNT composite, respectively. For both validations, the matrix thermal conductivity was set as the experimentally determined value and the inclusion dimensions were taken from the original experimental data. The model predictions were calculated by varying the thermal interface resistance in the range of 1 × 10 −8 to 8 × 10 −8 m 2 ·K/W. This was done to take into account the variation in thermal interface resistance that can occur due to change in process parameters, the source of raw materials and can even vary greatly from sample to sample [35]. Model predictions for both composites showed good agreement with experimental measurements.

Effect of Nanometer-Sized Inclusions
The generalized EMT formulated in the current work was applied to three particulate nanocomposites, Ge-Si composites, alumina-CNT composites and aluminum-CNT composites, to study the effect of nanometer-sized inclusions on the effective thermal conductivity of the composites.
The Ge-Si composite was studied for a case reported in the literature for spherical Si inclusions in Ge matrix. The thermal interface resistance, R, was calculated using the diffuse mismatch model given by Equation (33) and Si and Ge material properties were taken from [14] and are given in Table 1. Figure 12 shows a comparison of the effective thermal conductivities of Ge-Si composite predicted by the generalized EMT formulation with those predicted by Monte Carlo simulations [14]. The results showed that the effective thermal conductivity of Ge-Si nanocomposites is significantly lower than the thermal conductivities of Si and Ge. The effective thermal conductivity falls to less than 5 W/m·K for Si inclusions with a diameter equal to 10 nm and a volume fraction of 20%. This drastic reduction occurs because both Silicon and Germanium have mean-free-paths much greater than the size of the Si inclusion.  The results showed that the effective thermal conductivity of Ge-Si nanocomposites is significantly lower than the thermal conductivities of Si and Ge. The effective thermal conductivity falls to less than 5 W/m·K for Si inclusions with a diameter equal to 10 nm and a volume fraction of 20%. This drastic reduction occurs because both Silicon and Germanium have mean-free-paths much greater than the size of the Si inclusion.
The root-mean-squared error (RMSE) in model predictions for the nanocomposite with Si particles of radii 5 nm, 25 nm and 100 nm were 1.68 W/m·K, 2.18 W/m·K, and 4.45 W/m·K, respectively. These translated to normalized errors of 12.72%, 10.22%, and 12.25%, respectively. The RMSE for the case with inclusion size of 100 nm was higher than the other cases because the thermal conductivity of the composite was higher for this case. On the other hand, the normalized root mean squared error all three cases were almost the same. That being said, the difference in the results of the Monte Carlo simulation and the generalized EMT can be attributed to a number of reasons including the difference in the details included in the models and model domain size. Monte Carlo simulations involve the detailed simulation of the transport of a large number of energy carriers through the material whereas the generalized EMT only models the average effect of transport of energy carriers. But because of their complexity, Monte Carlo simulations are extremely computationally intensive and therefore are usually solved over a unit cell containing a single inclusion particle.
The generalized EMT was also applied to the case to alumina and aluminum matrix composites with randomly oriented multi-walled carbon nanotubes (MWCNT). A comparison between the generalized EMT and Nan and coworkers' model [36] for the two cases is shown in Figure 13. Nan and coworkers' model can take into account the effect of size and orientation of CNTs on the effective thermal conductivity of the composite but ignores the effect of CNT size on the thermal conductivity of matrix and CNT inclusion itself. The properties of alumina, aluminum, and MWCNTs used in the current work are shown in Table 2.  To investigate the reason for the difference in the effect of nanometer inclusion sizes on the effective thermal conductivity of alumina-CNT and aluminum-CNT composites, the variation of the  For comparison between the two models, the value of the thermal interface resistance, R, was varied from 0 to 1 × 10 −7 . For high values of thermal interface resistance, a very small difference was observed between the generalized EMT and Nan and coworkers' model for alumina-MWCNT composite. In this case, the effect of interfacial resistance dominated the effect of inclusion size and therefore, the difference is more prominent at lower values of R. On the other hand, the aluminum-MWCNT composites showed dependence on CNT size for all interfacial resistance values. A significant result of the comparison is that the theoretical maximum effective thermal conductivity possible when R = 0 is reduced from 85 W/m·K to 78 W/m·K for alumina matrix composite and from 298 W/m·K to 278 W/m·K for aluminum matrix composite when the CNT volume fraction is 5%.
To investigate the reason for the difference in the effect of nanometer inclusion sizes on the effective thermal conductivity of alumina-CNT and aluminum-CNT composites, the variation of the reduced thermal conductivities of alumina and aluminum matrices with CNT diameter was analyzed. The results are shown in Figure 14. As shown in the figure, the CNT diameter had a minimal effect on the alumina matrix thermal conductivity, which dropped by around 1 W/m·K when the CNT diameter was reduced from 20 nm to 5 nm. For the same change in CNT diameter, the thermal conductivity of the aluminum dropped from 236 W/m·K to 224.5 W/m·K. This large difference between aluminum and alumina in the sensitivity to CNT size can be attributed to the difference in the mean-free-paths of energy carriers for aluminum and alumina. For aluminum, the dominant energy carriers are the electrons whose mean-free-path is 14.1 nm compared to the 4.6 nm phonon mean-free-path for alumina. Since the CNT diameters were larger than the phonon mean-free-path for alumina, the addition of CNTs had little effect on the thermal conductivity of alumina.  To investigate the reason for the difference in the effect of nanometer inclusion sizes on the effective thermal conductivity of alumina-CNT and aluminum-CNT composites, the variation of the reduced thermal conductivities of alumina and aluminum matrices with CNT diameter was analyzed. The results are shown in Figure 14. As shown in the figure, the CNT diameter had a minimal effect on the alumina matrix thermal conductivity, which dropped by around 1 W/m·K when the CNT diameter was reduced from 20 nm to 5 nm. For the same change in CNT diameter, the thermal conductivity of the aluminum dropped from 236 W/m·K to 224.5 W/m·K. This large difference between aluminum and alumina in the sensitivity to CNT size can be attributed to the difference in the mean-free-paths of energy carriers for aluminum and alumina. For aluminum, the dominant energy carriers are the electrons whose mean-free-path is 14.1 nm compared to the 4.6 nm phonon mean-free-path for alumina. Since the CNT diameters were larger than the phonon meanfree-path for alumina, the addition of CNTs had little effect on the thermal conductivity of alumina.

Effect of Inclusion Aspect Ratio
The generalized EMT was also used to study the effect of inclusion aspect ratio on the effective thermal conductivity of randomly oriented spheroidal inclusions cos 2 θ = 1/3 using the example of Ge-Si nanocomposite. The results are shown in Figure 15. To compare the inclusions of various aspect ratios, it was assumed that the volume of a single inclusion is same for all cases. The results in Figure 15 show that for inclusions having same volume, the best thermal conductivity reduction was achieved for spherical inclusions. The reason for this was that spherical inclusions had the largest cross-sectional area to obstruct the path of the energy carriers in the matrix (1963.5 nm 2 for the inclusion of radius 25 nm as compared to 1119.2 nm 2 and 1392.5 nm 2 for inclusions of aspect ratios 0.2 and 5, respectively). This resulted in greater reduction in thermal conductivity for spherical inclusions.

Effect of Inclusion Aspect Ratio
The generalized EMT was also used to study the effect of inclusion aspect ratio on the effective thermal conductivity of randomly oriented spheroidal inclusions   2 cos 1 / 3   using the example of Ge-Si nanocomposite. The results are shown in Figure 15. To compare the inclusions of various aspect ratios, it was assumed that the volume of a single inclusion is same for all cases. The results in Figure 15 show that for inclusions having same volume, the best thermal conductivity reduction was achieved for spherical inclusions. The reason for this was that spherical inclusions had the largest cross-sectional area to obstruct the path of the energy carriers in the matrix (1963.5 nm 2 for the inclusion of radius 25 nm as compared to 1119.2 nm 2 and 1392.5 nm 2 for inclusions of aspect ratios 0.2 and 5, respectively). This resulted in greater reduction in thermal conductivity for spherical inclusions.

Effect of Inclusion Orientation
The effective thermal conductivity of aligned spheroidal nanometer-sized inclusions was estimated by Ordonez-Miranda and coworkers [19]. The generalized EMT presented in the current work extends their approach to composites with randomly oriented spheroidal inclusions. A comparison of the axial direction effective thermal conductivity of Ge-Si nanocomposites with aligned and randomly oriented inclusion is presented in Figure 16, for two different aspect ratios of inclusions. For oblate inclusions (p < 1), the alignment of inclusions in a specific direction resulted in reduced effective thermal conductivity, which is the main purpose of the development of Ge-Si nanocomposites. On the other hand, for prolate inclusions (p > 1), the effect was reversed and randomly oriented prolate inclusions were more effective in reducing the thermal conductivity of Ge-Si composites.
The difference in the behavior of oblate and prolate shaped inclusions was because of two reasons. First, the cross-sectional area of the inclusion obstructing the motion of energy carriers in the matrix was reduced due to random orientation of particles in case of oblate inclusions while it increased for prolate inclusions. This resulted in the matrix thermal conductivity to increase due to random orientation of oblate inclusion and reduce due to the random orientation of prolate inclusions. Second, the collision mean-free-path of the energy carriers inside the inclusion increased due to random orientation of oblate inclusions while it reduced due to the random orientation of prolate inclusions. This resulted in an increase in the inclusion thermal conductivity for oblate inclusions and a decrease in the inclusion thermal conductivity for prolate inclusions due to random orientation. Figure 15. Dependence of Ge-Si effective thermal conductivity on inclusion aspect ratio p for inclusion size a 1 equal to: (a) 5 nm; and (b) 25 nm (R TB = 6.80 × 10 −9 m 2 ·K/W).

Effect of Inclusion Orientation
The effective thermal conductivity of aligned spheroidal nanometer-sized inclusions was estimated by Ordonez-Miranda and coworkers [19]. The generalized EMT presented in the current work extends their approach to composites with randomly oriented spheroidal inclusions. A comparison of the axial direction effective thermal conductivity of Ge-Si nanocomposites with aligned and randomly oriented inclusion is presented in Figure 16, for two different aspect ratios of inclusions. For oblate inclusions (p < 1), the alignment of inclusions in a specific direction resulted in reduced effective thermal conductivity, which is the main purpose of the development of Ge-Si nanocomposites. On the other hand, for prolate inclusions (p > 1), the effect was reversed and randomly oriented prolate inclusions were more effective in reducing the thermal conductivity of Ge-Si composites.

Effect of Non-Uniform Dispersion of Inclusions
The two-scale approach presented in the current was used to study the effect of non-uniform dispersion on the effective thermal conductivity of Ge-Si nanocomposites. Without loss of generality, Ge-Si nanocomposites with 2.5%, 5%, 7.5% and 10% Si inclusions with a diameter equal to 10 nm were considered in the study. The results, presented in Figure 17, show the effective thermal The difference in the behavior of oblate and prolate shaped inclusions was because of two reasons. First, the cross-sectional area of the inclusion obstructing the motion of energy carriers in the matrix was reduced due to random orientation of particles in case of oblate inclusions while it increased for prolate inclusions. This resulted in the matrix thermal conductivity to increase due to random orientation of oblate inclusion and reduce due to the random orientation of prolate inclusions. Second, the collision mean-free-path of the energy carriers inside the inclusion increased due to random orientation of oblate inclusions while it reduced due to the random orientation of prolate inclusions. This resulted in an increase in the inclusion thermal conductivity for oblate inclusions and a decrease in the inclusion thermal conductivity for prolate inclusions due to random orientation.

Effect of Non-Uniform Dispersion of Inclusions
The two-scale approach presented in the current was used to study the effect of non-uniform dispersion on the effective thermal conductivity of Ge-Si nanocomposites. Without loss of generality, Ge-Si nanocomposites with 2.5%, 5%, 7.5% and 10% Si inclusions with a diameter equal to 10 nm were considered in the study. The results, presented in Figure 17, show the effective thermal conductivity of the Ge-Si nanocomposite normalized with effective thermal conductivity of Ge-Si nanocomposite with uniformly dispersed inclusions plotted against the non-uniformity in inclusion distribution (shown in the figure as the ratio of standard deviation, σ ϕ , and average inclusion volume fraction, ϕ Si,avg ). To analyze the variation of results due to randomness, each case was repeated five times by randomly generating an inclusion distribution in the RVE. As can be seen from the figure, the effectiveness of Si inclusions in reducing the effective thermal conductivity of the composite is adversely affected by non-uniformity in inclusion dispersion.

Effect of Non-Uniform Dispersion of Inclusions
The two-scale approach presented in the current was used to study the effect of non-uniform dispersion on the effective thermal conductivity of Ge-Si nanocomposites. Without loss of generality, Ge-Si nanocomposites with 2.5%, 5%, 7.5% and 10% Si inclusions with a diameter equal to 10 nm were considered in the study. The results, presented in Figure 17, show the effective thermal conductivity of the Ge-Si nanocomposite normalized with effective thermal conductivity of Ge-Si nanocomposite with uniformly dispersed inclusions plotted against the non-uniformity in inclusion distribution (shown in the figure as the ratio of standard deviation, σφ, and average inclusion volume fraction, φSi,avg). To analyze the variation of results due to randomness, each case was repeated five times by randomly generating an inclusion distribution in the RVE. As can be seen from the figure, the effectiveness of Si inclusions in reducing the effective thermal conductivity of the composite is adversely affected by non-uniformity in inclusion dispersion. Figure 17. Effect of non-uniform dispersion on Ge-Si effective thermal conductivity (RTB = 6.80 × 10 −9 m 2 ·K/W). Figure 17.

Sensitivity Analysis of the Generalized Effective Medium Theory
To identify the relative importance of the material properties and other parameters used in the generalized effective medium theory, a sensitivity analysis was carried. The model parameters considered in the sensitivity analysis were the matrix and inclusion thermal conductivities, mean-free-path of the energy carriers in the matrix and inclusion, the inclusion volume fraction and the thermal interface resistance. Using the case of Ge-20% Si nanocomposite as the nominal case, normalized sensitivity coefficients (NSCs) were calculated for each parameter. Normalized sensitivity coefficients, calculated using Equation (34) express the order of magnitude change in the analyzed function that will result from one order of magnitude change in the concerned parameter [37,38]. Using normalized sensitivity coefficients, one-to-one comparison between the model parameters can be carried out and parameters to which the model is more sensitive to can be determined.
where Y is nominal value of the function at nominal model parameters X i and ∆Y is the change in the function value to due to a change of ∆X i in the model parameter X i . For the case of Ge-Si nanocomposite, the results of the sensitivity analysis are presented in Table 3 for three inclusion sizes and the variation of the normalized sensitivity coefficients with inclusion size is presented in Figure 18. The results of the study showed that the NSCs of three model parameters, bulk matrix thermal conductivity, matrix phonon mean-free-path and inclusion volume fraction, were at least an order of magnitude higher than the NSCs of the remaining three parameters. This implies that the effective thermal conductivity of the Ge-Si nanocomposite is an order of magnitude more sensitive to changes in these three parameters. An analysis of Figure 18 also revealed several important points regarding the sensitivity of the generalized effective medium theory to the model parameters. First, the NSC of the bulk matrix thermal conductivity remained almost constant for inclusion radii from 5 nm to 100 nm. This means the sensitivity of the effective thermal conductivity of the composite is independent of the inclusion size. Second, with an increasing inclusion size, the NSCs of matrix phonon mean-free-path and the inclusion volume fraction decreased while the NSC of the thermal interface resistance showed a significant increase. Both inclusion volume fraction and matrix phonon mean-free-path affect the reduction in the matrix thermal conductivity due to nanometer sized inclusion. This suggests that with an increase in the inclusion size, the effect of interfacial thermal resistance starts to dominate the effect of reduced matrix thermal conductivity.  Figure 18. Variation of normalized sensitivity coefficients for: (a) bulk matrix thermal conductivity, matrix phonon mean-free-path and inclusion volume fraction; and (b) bulk inclusion thermal conductivity, inclusion phonon mean-free-path and thermal interface resistance.

Conclusions
In this paper, we presented a generalized effective medium theory for the estimation of the effective thermal conductivity of particulate nanocomposites. The formulated EMT has the capability of incorporating the effects of size, shape, orientation and dispersion non-uniformity of multiple inclusions on the estimated thermal conductivity of particulate composites. Several applications of the formulated EMT were also presented.
For the Ge-Si nanocomposite, it was found that spherical Si inclusions result in better effective thermal conductivity reduction in the nanocomposite. It was also found that aligned oblate inclusions result be better thermal conductivity reduction than randomly oriented oblate inclusions. The effect is reversed for the case of prolate inclusions for which randomly oriented prolate inclusions show better thermal conductivity reduction. Finally, the effective thermal conductivity was found to be strongly dependent on the dispersion uniformity of the inclusion particles for Ge-Si nanocomposites. The effectiveness of nanometer-sized Si particles in reducing the thermal conductivity of Ge matrix reduced with increasing non-uniformity in Si dispersion.
For alumina-MWCNT and aluminum-MWCNT nanocomposites, the effect of high scattering of energy carriers on the effective thermal conductivity increased with reducing interface thermal resistance. Between alumina and aluminum, aluminum showed a greater sensitivity to CNT size due to its relativity large electron mean-free-path.
A sensitivity analysis carried out to determine the relative effective of the model parameters on the effective thermal conductivity of Ge-Si nanocomposites showed that the model was an order of magnitude more sensitive to changes in the matrix thermal conductivity, matrix phonon mean-freepath and inclusion volume fraction than changes in inclusion thermal conductivity, inclusion phonon Figure 18. Variation of normalized sensitivity coefficients for: (a) bulk matrix thermal conductivity, matrix phonon mean-free-path and inclusion volume fraction; and (b) bulk inclusion thermal conductivity, inclusion phonon mean-free-path and thermal interface resistance.

Conclusions
In this paper, we presented a generalized effective medium theory for the estimation of the effective thermal conductivity of particulate nanocomposites. The formulated EMT has the capability of incorporating the effects of size, shape, orientation and dispersion non-uniformity of multiple inclusions on the estimated thermal conductivity of particulate composites. Several applications of the formulated EMT were also presented.
For the Ge-Si nanocomposite, it was found that spherical Si inclusions result in better effective thermal conductivity reduction in the nanocomposite. It was also found that aligned oblate inclusions result be better thermal conductivity reduction than randomly oriented oblate inclusions. The effect is reversed for the case of prolate inclusions for which randomly oriented prolate inclusions show better thermal conductivity reduction. Finally, the effective thermal conductivity was found to be strongly dependent on the dispersion uniformity of the inclusion particles for Ge-Si nanocomposites. The effectiveness of nanometer-sized Si particles in reducing the thermal conductivity of Ge matrix reduced with increasing non-uniformity in Si dispersion.
For alumina-MWCNT and aluminum-MWCNT nanocomposites, the effect of high scattering of energy carriers on the effective thermal conductivity increased with reducing interface thermal resistance. Between alumina and aluminum, aluminum showed a greater sensitivity to CNT size due to its relativity large electron mean-free-path.
A sensitivity analysis carried out to determine the relative effective of the model parameters on the effective thermal conductivity of Ge-Si nanocomposites showed that the model was an order of magnitude more sensitive to changes in the matrix thermal conductivity, matrix phonon mean-free-path and inclusion volume fraction than changes in inclusion thermal conductivity, inclusion phonon mean-free-path and thermal interface resistance. The analysis also showed that as inclusion size increases, the model becomes increasingly sensitive to variation in thermal interface resistance.