Mori-Tanaka Based Estimates of Effective Thermal Conductivity of Various Engineering Materials

The purpose of this paper is to present a simple micromechanics-based model to estimate the effective thermal conductivity of real-world macroscopically isotropic materials of matrix-inclusion type. The methodology is based on the well-established Mori-Tanaka method for composite media reinforced with ellipsoidal inclusions, extended to account for imperfect thermal contact at the matrix-inclusion interface, random orientation of particles and particle size distribution. Using simple ensemble averaging arguments, we show that the original Mori-Tanaka relations are still applicable for these complex systems, provided that the inclusion conductivity is appropriately modified. Such conclusion is supported by the verification of the model against a detailed finite-element study as well as its validation against experimental data for a wide range of engineering material systems.


Introduction
There has been a clear trend over the last decade to exploit ever greater detail of the material structure towards better predictions of its response from simulations. Hierarchical modeling strategy, regardless whether coupled or uncoupled but mostly of the bottom-up type, has served to provide estimates of the macroscopic response. In this process, geometric details decisive for a given scale are first quantified employing various statistical descriptors [1], but eventually smeared via homogenization to render larger scale property. Greater precision is expected when introducing the results of microstructure evaluation into the homogenization step. However, the actual gain when compared to the cost of this analysis is still in question. Obviously, description of evolving microstructures or rigorous representation of deformation mechanisms would require to account for almost every detail of the microstructure on a given scale. But how deep do we have to go if only the effective macroscopic response (i.e. linear macroscopic properties) is of the primary interest? Such a goal is addressed in this contribution.
Here, the modeling effort concentrates on the evaluation of effective thermal conductivities of various engineering materials with a significant degree of heterogeneity whereas focusing on imperfect thermal contact along constituents interfaces. We shall argue, shielded by available experimental data, that reasonably accurate predictions of macroscopic response can be obtained with very limited information about actual microstructure such as volume fractions and local properties of material phases. Consequently, we lump the entire analysis on the assumption of representing true material structures by statistically isotropic distribution of spheres. Figure 1 shows micro-images of selected material representatives which seem to admit this classification. Note that whatever material phase embedded into the matrix (the basic material) is henceforth termed the heterogeneity in real material systems while it is termed inclusion in approximations adopted for calculations. Strong motivation for this seemingly swingeing simplification is supported by experimental measurements presented in [2] for cement matrix based mixture of rubber particles and air voids. Comparison between experimental data and predictions provided by the Mori-Tanaka averaging scheme under the premise of random distribution of spherical inclusions, the method in this particular format promoted herein, appears in Figure 2. The match is almost remarkable.
Going back to Figure 1 one may object that while admissible for Figure 1(a) the spherical approximation of particle phase in Figure 1(c) will yield rather erroneous predictions. Note, however, that this attempt is not hopeless providing the microstructure can still be considered as macroscopically isotropic, ensured for statistically isotropic distribution of heterogeneities having isotropic material symmetry. In that case, it can be shown that the previously mentioned Mori-Tanaka method written out for spherical inclusions is adequate providing the material properties of the inclusions are suitably modified. Although this step requires information beyond that of volume fractions of phases, the benefit of gathering additional data will become particularly appreciable once turning our attention to material systems with imperfect interfaces, the principal objective of this study. The problem of quantifying the influence of imperfect thermal contact on the overall thermal conductivity has been under intense study in the past. Hasselman and Johnson [3] provided estimates for dilute concentration of mono-disperse spherical and cylindrical heterogeneities. Successful application of this simple model to Al/SiC porous composites is presented in [4]. The Hasselman-Johnson results were then extended by Benveniste and Miloh [5] to spheroidal particle shapes with imperfect interfaces and subsequently applied in the framework of the Mori-Tanaka method [6]. These early developments were later generalized by Nogales and Böhm, who proposed in [7] a simple method for dealing with polydisperse systems of spherical particles. In addition, rigorous third-order bounds for effective conductivity of macroscopically isotropic distribution of particles with imperfect interfaces were derived by Torquato and Rintoul [8]. Alternatively, as demonstrated by Hashin [9], the material systems with imperfect interfaces can be accurately approximated by the coated inclusion model due to Dunn and Taya [10], which also accounts for different orientation of the inclusions. If limiting attention to spherical inclusions the results presented in [7] can be obtained in a very elegant way by simple extension of one-dimensional analysis. This is demonstrated in Appendix B.
To exploit this result in practical applications of the Mori-Tanaka method to a heat conduction problem, prediction of effective thermal conductivity in particular, we adopt the analysis scheme graphically presented in Figure 3. We start from the assumption of multidisperse system of randomly oriented spheroidal inclusions with possibly imperfect thermal contact (non-zero temperature jump along the interface). To arrive at the desired approximation of multidisperse system of spherical inclusions with perfect interfaces (temperature continuity along the interface) we proceed in five consecutive steps. These steps are mathematically described in Section 2. Section 3 is devoted to both validation and verification of the proposed scheme against available experimental data and finite element simulations performed for several representatives of statistically isotropic random microstructures. The crucial results and principal recommendations are finally summarized in Section 4.

Theoretical background of the Mori-Tanaka method
In this section attention is accorded to essential theoretical details of the Mori-Tanaka method in view of the five steps in Figure 3. In the first step, we consider a single inclusion with perfect interface subject to far-field loading. This step is theoretically elaborated in Section 2.1. Solution of this problem is then employed in Section 2.2 to estimate the overall conductivity of a composite consisting of multiple ellipsoidal inclusions bonded to a matrix phase. The third step addressed in Section 2.3 is reserved for systems with randomly oriented inclusions with a uniform distribution over the hemisphere. Here, a simple orientation averaging argument is shown to demonstrate that the effective conductivity of the system coincides with the conductivity of a system reinforced by spherical inclusions, thermal conductivity of which is appropriately modified. Following [7], an analogous argument is employed in the next step outlined in Section 2.4 to account for imperfect thermal contact along the matrix-inclusion interface. It is shown that in this case the modified conductivity becomes size-dependent. This eventually allows us to extend the scheme to polydisperse systems in Section 2.5.

Single inclusion with perfect interface
Let us consider an ellipsoidal inclusion Ω i , with semi-axes a 1 ≤ a 2 ≤ a 3 embedded in the matrix domain Ω m . We attach to the inclusion a Cartesian coordinate system with the origin at the inclusion center and axes aligned with the semi-axes. The distribution of local fields follows from the problem where q ∈ R 3 denotes the heat flux, h ∈ R 3 denotes gradient of the temperature θ (i.e. h(x) = ∇θ(x)) and χ designates the 3 × 3 symmetric positive-definite matrix of thermal conductivity given by Eq. (1) is completed by the far-field boundary conditions, cf. [5,Eq. (14)], with H ∈ R 3 denoting the overall (macroscopic) temperature gradient. Due to linearity of the problem, we can introduce the temperature gradient concentration factor A ∈ R 3×3 in the form: As shown first by Hatta and Taya [11], the concentration factor is constant inside the inclusion and admits the expression where I denotes the unit matrix and S ∈ R 3×3 is the Eshelby-like tensor which depends only on the matrix conductivity χ m and the ratios of semi-axes lengths β 2 = a 2 : a 1 and β 3 = a 3 : a 1 , see also Appendix A for additional details. For the spherical inclusion with isotropic conductivity χ i = χ i I embedded in the isotropic matrix with χ m = χ m I, Eq. (5) simplifies as

Multiple inclusions with perfect interface
In the next step, we adopt the results of the previous section to estimate the overall behavior of a composite material consisting of distinct phases r = 0, 1, . . . , N . The value r = 0 is reserved for the matrix phase Ω m and the r-th phase (r > 0) corresponds to the ellipsoidal inclusion Ω (r) , characterized by its semi-axes a 3 , volume fraction c (r) and conductivity χ (r) . Following Benveniste's reformulation [12] of the original Mori-Tanaka scheme [13], the interaction among phases is approximated by subjecting each inclusion separately to the mean temperature gradient in the matrix phase H m in Eq. (3). As a result, the temperature gradient inside the r-th phase remains constant and reads H (r) = T (r) H m for r = 0, 1, . . . , N. Here, and T (r) is the 3 × 3 partial temperature gradient concentration factor of the r-th phase, given by where the 3 × 3 rotation matrix R (r) accounts for the difference in the global and local coordinate systems, see Section 2.3 for additional details, and T (r) equals to A i in Eq. (5) with χ i = χ (r) and S determined from values a Inverting the above equation gives the average temperature gradient in the matrix as which, when substituted into Eq. (7), yields the explicit expression for the phase temperature gradients in the form where A m , A (r) are the matrix and inclusion temperature gradient concentration factors, respectively. As each phase is assumed to be homogeneous, the average heat flux in the r-th phase equals to This allows us to express the macroscopic heat flux in the form from which we obtain the effective conductivity Q = −χ H H in its final form Finally, assuming that the composite consists of isotropic matrix with conductivity χ m and spherical isotropic inclusions with conductivities χ (r) , Eq. (16) becomes χ H = χ H I, with

Orientation averaging
We are now in a position to provide estimates of the effective thermal conductivity for composites with M (with M N ) inclusion classes indexed by s = 1, 2, . . . , M . Each class is characterized by a single Eshelby-like matrix S in Eq. (5) and represents the reference ellipsoidal inclusion randomly oriented over the unit hemisphere with an independent uniform distribution of orientation angles.
To this goal, consider a quantity X ∈ R 3×3 , expressed in a local coordinate system aligned with a certain reference inclusion. Its form in the global coordinate system follows from where ϕ, φ and ψ denote the Euler angles 1 and the transformation matrix R is provided by The orientation average of X is denoted by double angular brackets: Straightforward calculation, presented e.g. in [14,Appendix A.2.3], reveals that the orientation averaging of an arbitrary X ∈ R 3×3 yields Repeating the steps of the previous section with partial temperature gradient concentration factors replaced with their orientation averages, we obtain, after some manipulations presented e.g. in [15, Section B], the scalar homogenized conductivity in the form Therefore, it follows from the comparison of Eq. (22) with Eq. (17) that the system of randomly oriented inclusions embedded in an isotropic matrix is indistinguishable, from the point of view of homogenized conductivity, from the system of spherical inclusions with an apparent conductivity which yields

Imperfect interface
The presence of imperfect thermal contact at the matrix-inclusion interface ∂Ω i results in temperature jump θ , the magnitude of which is provided by Newton's law, e.g. [16, Section 1.3]: where k denotes the interfacial conductance (with k → ∞ corresponding to perfect interface and k = 0 to ideal insulation) and n denotes the normal vector at the interface oriented outside the inclusion. This relation, together with Eqs. (1)-(3), defines the single inclusion problem accounting for the presence of imperfect interface. Its solution is, however, substantially more involved as the temperature gradient inside an ellipsoidal inclusion becomes position-dependent; the concentration factor is then available only in the form of complicated infinite series expansion for spheroidal inclusions [5] or ellipsoidal coated inclusions [10]. Nevertheless, when restricting the attention to spherical inclusion of radius a, it can be shown that the temperature gradient within inclusion recovers the constant value and the concentration factor becomes [5] A see also Appendix B for simple derivation of this result. Hence, analogously to the previous section, the interfacial effect can be modeled by replacing the "true" conductivity χ i by the size-dependent apparent value χ i provided by Eq. (26) 2 . Assuming in addition that each class of inclusions is characterized by identical semi-axes lengths a 3 and the apparent conductivity χ (s) given by Eq. (23).

Polydisperse systems
Even though Eq. (27) is applicable to very general material systems, in practice we typically assume single inclusion family, with the particle size distribution characterized by a probability density function p(a) satisfying In this context, the effective conductivity finally becomes where, for an arbitrary function g(a), {g} denotes its expected value given by Following e.g. [7,17], the log-normal distribution with the probability density function will be employed to characterize materials' polydispersity. The parameters µ and σ are provided by [7, Eq. (17)] where a x denotes the x-th percentile of the particle radii and S is the span of the size distribution.

Validation against available experimental data
Two particular examples of real engineering materials are examined in this section to show applicability of Eq. (27) and its extension for polydisperse distribution of heterogeneities, Eq. (29), even when disregarding their actual shape and simply accepting a spherical representation of the inclusions in the Mori-Tanaka predictions. The results provided by these two equations are corroborated by available experimental data.

Random dispersion of copper particles in the Epoxy matrix
In this first example we compare the single-phase Mori-Tanaka predictions with the experimental results of de Araujo and Rosenberg [18] who measured the effective thermal conductivities of systems consisting of random dispersion of metal particles in the epoxy matrix for several values of the interfacial resistance arising due to acoustic mismatch at the particle-matrix interface particularly for low temperatures below 20K. To enhance credibility of the Mori-Tanaka predictions we focus on one particular system made from epoxy resin filled with copper particles also examined in [8] in view of the threepoint lower bound assuming a random array of superconducting hard spheres (i.e. χ (1) → ∞). It has been shown, see [18,19], that for metal-filled composites with ratio α = χ (1) /χ m > 10 2 the macroscopic conductivity does not depend on the thermal conductivity of the metallic filler, but only on its volume fraction together with the properties of matrix and matrix-particle interface. This becomes evident when rewriting Eq. (27) 1 in terms of α where we introduced the dimensionless quantity R adopted in [8] where a is the sphere radius. It follows from Eq. (33) that for the inclusion size equal to the Kapitza radius [20] the effective conductivity equals to that of the matrix, thus the effect of inclusions becomes completely shielded by the interface. For a < a K , the overall properties are dominated by interfaces and the effective conductivity decreases with increasing c (1) even when α > 1. For a > a K , the bulk properties of phases become dominant; see also [20] for further discussion.
In addition to experimental measurements, we also present a comparison with the Torquato-Rintoul three-point lower bound in resistance case derived in [8,Eqs. (8)], evaluated for the statistical parameter ζ(c (1) ) given for the model of impenetrable spheres in [21,Table II (simulations)]. Note that the interfacial properties can be estimated by measuring the ratio of temperature jump to the applied heat flux across a thin bi-material layer, by acoustic mismatch model [4] or, indirectly, from an inverse approach as discussed below. The results are presented for two different temperatures θ = 4K (a = 14.8a K ) and θ = 3K (a = 4.93a K ). Lower bound for perfect interface MT for perfect interface Exp [18] The influence of parameter α together with expected particle size dependence, now hidden in parameter R through Eq. (34) 1 , is evident from Figure 4(a) plotted for θ = 4K (a K = 3.38 µm). Clearly, the Mori-Tanaka predictions confirm negligible influence of χ (1) observed for particulate composites with χ (1) χ m (α > 10 2 ) as well as decreasing trend for χ H with decreasing particle size caused by imperfect thermal contact. Figure 4(b) then displays evolution of normalized effective thermal conductivity as a function of volume fraction of copper particles. Predictive capability of the Mori-Tanaka method is supported here by a very good agreement with both experimental measurements and the Torquato-Rintoul bounds [8]. Finally, Figure 4(c) shows correlation between theoretical (MT predictions) and experimental results. The solid line was derived by linear regression of measured and calculated effective conductivities using the least square method. Another indication of the quality of numerical predictions can be presented in the form of Pearson's correlation coefficient written as where a = n i=1 a i , n is the number of measurements, x stands for the experimental and y for the corresponding theoretical values. Note that for ρ = 1 the correspondence is exact. In this case the correlation coefficient equals to 0.999 suggesting almost perfect match between measured and predicted values, also evident from graphical presentation.

Al/SiC composite
In [4] the authors studied the effect of imperfect thermal contact on the macroscopic response of Al/SiC porous composites. The paper presents the results of a thorough experimental investigation and traces of an inverse approach in material mechanics for inferring material properties of unknown components of the composite by matching numerical and experimental results. This approach was first exploited to derive the matrix thermal conductivity from known electrical conductivity of the composite. Next, the Hasselman and Johnson model [3] was employed under the premise of random distribution of spherical particles of identical size to estimate the particle thermal conductivity and interfacial thermal conductance for pore-free specimens and subsequently utilized in the two-step application of the Hasselman-Johnson model to address the influence of pores. Note that the material data used in these predictions can be thought as optimal with respect to the adopted Hasselman-Johnson model.
Hereinafter, we compare the results presented in [4] for seven specimens with pore-free matrix. In addition, we take advantage of available characteristics of the SiC particles, the span S and the 50-th percentile a 50 , to extend the analysis to polydisperse systems as presented in Section 2.5. The input material data are listed in Table 1. The Mori-Tanaka predictions stored in Table 2 were provided by Eq. (29). The integral (30) was evaluated such that the entire interval was split into 1,000 segments, thus considering 1,000 different particle sizes of spherical shape. Within each segment s the probability function p(a) was approximated by a straight line and the volume fraction c (1) of a given set of particles, given by the segment area, was then associated in a logarithmic way with the mean radius a (1) of this set of particles. Standard trapezoidal integration rule was then used to sum over all 1,000 segments. Examples of probability distribution functions for specimens No. 1 and 7 are plotted in Figure 5.
Graphical representation of the results is further seen in Figure 5(b),(c) confirming the sensitivity of the effective thermal conductivity to the mean particle size distribution. Note that individual points in Figure 5(b) correspond to slightly different volume fraction, see Table 2.  Almost negligible deviation from experimental results measured as is, however, not surprising owing to the used material parameters, which were not measured but rather fitted to a micromechanical model. Even when comparing the Pearson correlation coefficients, 0.98 for the Hasselman-Johnson model and 0.993 for polydisperse MT model, the improvement when accounting for more accurate representation of particle size distribution is marginal. This can be explained by a very small variance associated with adopted distributions, recall Figure 5(a). Nevertheless, it is fair to point out that unlike the Hasselman-Johnson model the Mori-Tanaka approach is not limited to spherical particles providing the transformations given by Eqs. (24) and (27) are admissible. The influence of shape of particles on the macroscopic response has been put forward, e.g. by Jäckel in [19], and is numerically investigated in the next section suggesting increasing thermal conductivity of the composite with transition from spherical to needle-like particles, the trend also observed experimentally [18,19].

Verification against finite element simulations
It has been argued in the previous sections that even very limited information about microstructure amounted to phase properties and corresponding volume fractions might be sufficient to provide a reasonable estimate of macroscopic response of various engineering material systems generally classified as being macroscopically isotropic. This naturally invites the assumption of spherical representation of otherwise irregular heterogeneities. Although supported by several practical examples discussed in the previous section, we should expect and even identify, at least qualitatively, limitations to such perception. In doing so, this section presents numerical investigation of some specific issues such as the influence of shape and size of inclusions or mismatch of phase material properties on the predicted macroscopic response.
All numerical results reported bellow are obtained in the two-dimensional setting, hence the arguments presented in Section 2 need to be translated to the planar case. In particular, inclusions are modeled as elliptic cylinders of aspect ratio β 2 = a 2 /a 1 with a 3 → ∞. The corresponding Eshelby-like tensor is given by Eq. (42), which leads to the two-dimensional thermal concentration factor for circular inclusion in the form The apparent conductivity due to random orientation is provided by where • 2D stands for the orientation averaging for the uniform distribution of the Euler angle φ. Finally note that the apparent conductivity due to imperfect interface is provided by the same relation as in the three-dimensional setting, compare Eq. (26) with Eq. (50). Figure 6. Examples of random macroscopically isotropic microstructures: a) circular cylinders (β 2 = 1), b) elliptical cylinders with aspect ratio β 2 = 3, c) elliptical cylinders with aspect ratio β 2 = 9.
Three particular representatives, generated such as to approximately resemble the real microstructures in Figure 1, appear in Figure 6. To comply with general assumptions put forward in the previous sections, we consider locally isotropic phases with variable contrast of material properties. Additionally, we assume the above microstructures being periodic and adopt the classical first-order homogenization strategy, see e.g. [22,23], to provide estimates of the macroscopic response. The results are plotted in Figure 7. Figure 7. Variation of the normalized effective conductivity χ H /χ m for three microstructures in Figure 6 with perfect interfaces, determined by periodic FEM homogenization for phase contrast a) α = 3, b) α = 10 and c) α = 20; β 2 denotes the inclusion aspect ratio and c (1) the inclusion volume fraction.
These results clearly indicate not only the influence of the shape of inclusions on the macroscopic response but also a strong dependence of these predictions on the contrast of material properties of individual phases. Thus drawing from the plots presented in Figure 7(a) one may suggest that the proposed circular representation of generally non-circular heterogeneities is still acceptable when their shapes only moderately deviate from a circle and when the mismatch of phase properties is not too severe, which certainly is the case of a number of real materials as demonstrated in the previous section. This expectation is quite important particularly when dealing with imperfect thermal contact in which case only spherical and circular inclusions can be easily handled analytically.
If the circular approximation of heterogeneities is no longer acceptable or the contrast of phase properties is excessive, one needs to look for more details about microstructure. In such a case, even twodimensional images of real systems, at present almost standard input for any material based analysis, may play an important role in assessing better approximations of shapes, say elliptical, of these heterogeneities, see e.g. [24]. Then, being given the elliptical shape of the inclusion allows us to appropriately modify its material data, recall Eq. (38), and define a certain indicator of the real microstructure k corr , e.g. as a ratio of the modified and original partial temperature gradient concentration factors now determined for circular inclusions.
Variation of this parameter as a function of the shape of inclusion is seen in Figure 8(a) further confirming quite strong influence of the phase properties mismatch. 2 The modified conductivity when introduced successively into Eq. (37) 2 and Eq. (24) 1 then renders the estimate of effective conductivity almost identical to actual microstructure with non-circular inclusions as evident from plots in Figure 8(b),(c). Note that only the first and the third microstructure in Figure 6 were examined to first confirm that the Mori-Tanaka method is indeed well suited for statistically isotropic random microstructures and second to promote applicability of this simple transformation from elliptical to circular representations even for shapes markedly distinct from circles. Small but evident deviation of the results observed in Figure 6 for the third needle-like microstructure and large mismatch of conductivities of the inclusion and matrix equal to 20 can be attributed to the finite size, although infinite in the sense of periodicity, of the representative model not large enough to yield statistically isotropic microstructure. The second set of numerical simulations addresses the theoretically derived dependence of macroscopic predictions on the size of inclusions in cases with imperfect interfaces generating jumps in the local temperature field. For simplicity, we limit our attention to statistically isotropic distributions of circular cylinders with a radius varying from sample to sample. Three such microstructures are shown in Figure 9. Note that the same volume fraction and the same number of inclusions was maintained in all simulations. Zero thickness interface elements were introduced to account for imperfect thermal contact. Figure 10. a) Variation of normalized effective conductivity χ H /χ m for a system with imperfect interface (k = 60 × 10 6 Wm −2 K −1 and χ (1) = 100 Wm −1 K −1 , χ m = 10 Wm −1 K −1 and c (1) = 40%) as a function of inclusion radius a and b) variation of effective conductivity for systems weakened by cylindrical voids as a function of inclusion volume fraction c (1) ; β 2 denotes inclusion aspect ratio.
The relevant results appear in Figure 10(a). Both the results found from finite element simulations and corresponding Mori-Tanaka predictions are displayed to clearly identify the mentioned size dependence. Proper modifications in the sense of Eq. (38) now become even more important as indicated by the results generated for elliptical microstructures from Figure 6 with cross-sectional area equal to the area of the circle. These are indicated by x-symbol and the ratio of semi-axes of the corresponding elliptical cylinder. The Mori-Tanaka estimates found from the application of Eqs. (38) and (26) are reasonably close further supporting the proposed approach for the modeling of real materials with an imperfect thermal contact. Intuitively, it can be expected that the value of interfacial conductance k may also show some effect as to the estimates of effective conductivities for non-circular inclusions. This notion is supported by the results presented in Figure 10(b) showing variation of effective conductivity of an isotropic matrix weakened by voids with a very low conductivity. Clearly, the influence of shape of the inclusions is quite pronounced.

Conclusions
The Mori-Tanaka micromechanical model has been often the primary choice among engineers to provide quick estimates of the macroscopic response of generally random composites. Motivated by an early theoretical as well as experimental works on this subject, the Mori-Tanaka method was examined here in the light of the solution of a linear steady state heat conduction problem allowing us to estimate the effective thermal conductivity of variety of real engineering materials experiencing an imperfect thermal contact along the constituents interfaces.
Adhering to the only limitation, an assumption of macroscopically isotropic composite, it was shown that the method originally proposed by Böhm and Nogales [7] for a spherical representation of particles still applies even to non-spherical particles providing their shape can be suitably quantified, e.g. by an ellipsoidal inclusion. In this particular case the Mori-Tanaka predictions were partially corroborated by two-dimensional numerical simulations confirming experimentally observed considerable sensitivity of macroscopic conductivities to the shape of particles.
The fact that for composites with imperfect thermal contact the macroscopic predictions depend on particle size can be effectively handled by introducing the particle size probability density function directly into the Mori-Tanaka estimates. Although not confirmed for material systems studied in the paper, this may considerably improve final predictions especially for grading curves showing significant standard deviation of particle sizes from its mean value. This is particularly appealing, since grading curves are one of the few information supplied by the manufacturer.
To conclude, it is interesting to point out that there exist many material systems that can be handled very effectively with simple micromechanical models with no need for laborious finite element simulations of certain representative volumes of real microstructures. 25. Chen, T.; Yang, S.H. The problem of thermal conduction for two ellipsoidal inhomogeneities in an anisotropic medium and its relevance to composite materials. Acta Mechanica 1995, 111, 41-58.

A. Eshelby-like tensor
The Eshelby-like tensor for the solution of thermal conductivity problem was introduced by Hatta and Taya in [11]. For an ellipsoidal inclusion with semi-axes a 1 , a 2 , a 3 found in an isotropic matrix it receives the form S ij = a 1 a 2 a 3 4 ∂ ∂x i ∂x j ∆s = (a 2 1 + s) (a 2 2 + s) (a 2 3 + s).
Closed form solutions of integral (39) for some special cases of ellipsoidal shapes of the inclusion can be found in [11]. For a general ellipsoid the solution was introduced by Chen and Yang in [25]. For circular and spherical shapes needed in the present study the S tensor reads • Sphere (a 1 = a 2 = a 3 ) S 11 = S 22 = S 33 = 1 3 and S ij = 0 for i = j.

B. Single spherical homogeneity with imperfect interface
This section outlines derivation of the replacement conductivity χ i and the concentration factor A sph introduced in Eq. (26). Since used in numerical calculations in Section 3.2, its two dimensional format is presented as well. It it shown that both 2D and 3D concentration factors can be recovered from the solution of a 1D problem using a simple geometrical argument.  To that end, consider one-dimensional heat conduction problem depicted in Figure 11(a). Assuming imperfect thermal contact, the temperature drop across an infinitely thin interface layer is given by Eq. (25). The local temperature gradient for perfect interface between a solitary inclusion embedded into an infinite matrix follows from Eq. (4) To arrive at similar format of Eq. (43) for imperfect contact, we imagine the interface temperature jump being smeared over the inclusion. Since the heat flux Q associated with the macroscopic temperature gradient H is constant throughout the composite, we obtain the total temperature change in the substitute inclusion in the form where L stands for the inclusion length. Next, defining the local temperature gradient in the substitute inclusion as H i = ∆ θ i /L yields the local constitutive law in terms of the replacement conductivity χ i so that and in analogy with Eq. (43) we finally get The two-dimensional problem of a solitary circular inclusion is treated similarly. We build up on the assumption that the temperature gradient in the inclusion is constant and collinear with the prescribed far field gradient parallel to the local x 1 -axis, see Figure 11(b). To draw similarity with 1D case we divide the inclusion into parallel filaments with the length L = 2a cos ϕ. Next, define a unit vector normal to the inclusion surface n = (cos ϕ, sin ϕ) T and in analogy to Eq. (44)) write the total temperature change in each filament for the constant heat flux q i = q i , 0 T as where d is the inclusion diameter. The equivalent conductivity has to fulfill the condition which yields Consequently, the concentration factor of the substitute inclusion attains the form