Mori-Tanaka Based Estimates of Effective Thermal Conductivity of Various Engineering Materials
SSubmitted to
Micromachines . Pages 1 - 21.
OPEN ACCESS micromachines
ISSN 2072-666X
Article
Mori-Tanaka Based Estimates of Effective ThermalConductivity of Various Engineering Materials
Jan Str´ansk´y , Jan Vorel , Jan Zeman and Michal ˇSejnoha , ,(cid:63) Department of Mechanics, Faculty of Civil Engineering, Czech Technical University in Prague,Prague, Czech Republic Centre for Integrated Design of Advanced Structures, Faculty of Civil Engineering, Czech TechnicalUniversity in Prague, Prague, Czech Republic (cid:63)
Author to whom correspondence should be addressed; E-Mail: [email protected], Tel.:+420-2-2435-4494 and Fax:+420-2-2435-4494
Version August 8, 2011 submitted to
Micromachines . Typeset by L A TEX using class file mdpi.cls
Abstract:
The purpose of this paper is to present a simple micromechanics-based model toestimate the effective thermal conductivity of macroscopically isotropic materials of matrix-inclusion type. The methodology is based on the well-established Mori-Tanaka method forcomposite media reinforced with ellipsoidal inclusions, extended to account for imperfectthermal contact at the matrix-inclusion interface, random orientation of particles and parti-cle size distribution. Using simple ensemble averaging arguments, we show that the Mori-Tanaka relations are still applicable for these complex systems, provided that the inclusionconductivity is appropriately modified. Such conclusion is supported by the verification ofthe model against a detailed finite-element study as well as its validation against experimen-tal data for a wide range of engineering material systems.
Keywords:
Engineering materials; homogenization; Mori-Tanaka method; thermal conduc-tivity; imperfect interface; finite element simulations
1. Introduction
There has been a clear trend over the last decade to exploit ever greater detail of the material structuretowards better predictions of its response from simulations. Hierarchical modeling strategy, regardlesswhether coupled or uncoupled but mostly of the bottom-up type, has served to provide estimates of the a r X i v : . [ c ond - m a t . m t r l - s c i ] A ug ersion August 8, 2011 submitted to Micromachines macroscopic response. In this process, geometric details decisive for a given scale are first quantifiedemploying various statistical descriptors [1], but eventually smeared via homogenization to render largerscale property. Greater precision is expected when introducing the results of microstructure evaluationinto the homogenization step. However, the actual gain when compared to the cost of this analysis is stillin question. Obviously, description of evolving microstructures or rigorous representation of deformationmechanisms would require to account for almost every detail of the microstructure on a given scale. Buthow 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 var-ious engineering materials with a significant degree of heterogeneity whereas focusing on imperfectthermal 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 infor-mation 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 bystatistically isotropic distribution of spheres. Figure 1 shows micro-images of selected material repre-sentatives which seem to admit this classification. Note that whatever material phase embedded into thematrix (the basic material) is henceforth termed the heterogeneity in real material systems while it istermed inclusion in approximations adopted for calculations.
Figure 1.
Examples of micro-graphs of real engineering materials taken in back scatteredelectrons: a) Alkali-activated fly ash, b) Alumino-silicate ceramics with Fe and silicium par-ticles (dark phase), c) Superspeed - alloying ingredient into crude iron for cast iron workingwith silicon particles (dark phase). Reproduced with permission of L. Kopeck´y (CTU inPrague). (a) (b) (c)Strong motivation for this seemingly swingeing simplification is supported by experimental measure-ments presented in [2] for cement matrix based mixture of rubber particles and air voids. Comparisonbetween experimental data and predictions provided by the Mori-Tanaka averaging scheme under thepremise of random distribution of spherical inclusions, the method in this particular format promotedherein, 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 approxi-mation of particle phase in Figure 1(c) will yield rather erroneous predictions. Note, however, that thisattempt is not hopeless providing the microstructure can still be considered as macroscopically isotropic,ersion August 8, 2011 submitted to
Micromachines a) Evolution of effective thermal conductivity χ H as a function of volume fractionof rubber in solid phase, b) correlation of measured and calculated values; ρ is the correlationcoefficient. Volume fraction of rubber [%] χ H [ W / m K ] Measured value [W/mK] C a l c u l a t ed v a l ue [ W / m K ] ρ = 0.987 (a) (b)ensured for statistically isotropic distribution of heterogeneities having isotropic material symmetry. Inthat case, it can be shown that the previously mentioned Mori-Tanaka method written out for sphericalinclusions is adequate providing the material properties of the inclusions are suitably modified. Al-though this step requires information beyond that of volume fractions of phases, the benefit of gatheringadditional data will become particularly appreciable once turning our attention to material systems withimperfect interfaces, the principal objective of this study.The problem of quantifying the influence of imperfect thermal contact on the overall thermal conduc-tivity has been under intense study in the past. Hasselman and Johnson [3] provided estimates for diluteconcentration of mono-disperse spherical and cylindrical heterogeneities. Successful application of thissimple model to Al/SiC porous composites is presented in [4]. The Hasselman-Johnson results werethen extended by Benveniste and Miloh [5] to spheroidal particle shapes with imperfect interfaces andsubsequently applied in the framework of the Mori-Tanaka method [6]. These early developments werelater generalized by Nogales and B¨ohm, who proposed in [7] a simple method for dealing with polydis-perse systems of spherical particles. In addition, rigorous third-order bounds for effective conductivityof macroscopically isotropic distribution of particles with imperfect interfaces were derived by Torquatoand Rintoul [8]. Alternatively, as demonstrated by Hashin [9], the material systems with imperfect inter-faces can be accurately approximated by the coated inclusion model due to Dunn and Taya [10], whichalso accounts for different orientation of the inclusions. If limiting attention to spherical inclusions theresults presented in [7] can be obtained in a very elegant way by simple extension of one-dimensionalanalysis. This is demonstrated in Appendix B.To exploit this result in practical applications of the Mori-Tanaka method to a heat conduction prob-lem, prediction of effective thermal conductivity in particular, we adopt the analysis scheme graphicallypresented in Figure 3. We start from the assumption of multidisperse system of randomly orientedspheroidal inclusions with possibly imperfect thermal contact (non-zero temperature jump along the in-terface). To arrive at the desired approximation of multidisperse system of spherical inclusions withperfect interfaces (temperature continuity along the interface) we proceed in five consecutive steps.ersion August 8, 2011 submitted to Micromachines
Mori-Tanaka based scheme: Strategy of derivation ≈ ⇒ ⇒⇒ ⇒ ⇒ (i) Single inclusion (ii) Multiple inclusions(iii) Orientation averaging (iv) Imperfect interface (v) PolydispersivityOriginal system
These steps are mathematically described in Section 2. Section 3 is devoted to both validation andverification of the proposed scheme against available experimental data and finite element simulationsperformed for several representatives of statistically isotropic random microstructures. The crucial re-sults and principal recommendations are finally summarized in Section 4.
2. Theoretical background of the Mori-Tanaka method
In this section attention is accorded to essential theoretical details of the Mori-Tanaka method in viewof the five steps in Figure 3. In the first step, we consider a single inclusion with perfect interface sub-ject to far-field loading. This step is theoretically elaborated in Section 2.1. Solution of this problem isthen employed in Section 2.2 to estimate the overall conductivity of a composite consisting of multipleellipsoidal inclusions bonded to a matrix phase. The third step addressed in Section 2.3 is reserved forsystems with randomly oriented inclusions with a uniform distribution over the hemisphere. Here, a sim-ple orientation averaging argument is shown to demonstrate that the effective conductivity of the systemcoincides with the conductivity of a system reinforced by spherical inclusions, thermal conductivity ofwhich is appropriately modified. Following [7], an analogous argument is employed in the next stepoutlined in Section 2.4 to account for imperfect thermal contact along the matrix-inclusion interface. Itis shown that in this case the modified conductivity becomes size-dependent. This eventually allows usto extend the scheme to polydisperse systems in Section 2.5.
Let us consider an ellipsoidal inclusion Ω i , with semi-axes a ≤ a ≤ a embedded in the matrixdomain Ω m . We attach to the inclusion a Cartesian coordinate system with the origin at the inclusioncenter and axes aligned with the semi-axes. The distribution of local fields follows from the problem q ( x ) = − χ ( x ) h ( x ) , ∇ T q ( x ) = 0 for x ∈ R , (1)ersion August 8, 2011 submitted to Micromachines where q ∈ R denotes the heat flux, h ∈ R denotes gradient of the temperature θ (i.e. h ( x ) = ∇ θ ( x ) )and χ designates the × symmetric positive-definite matrix of thermal conductivity given by χ ( x ) = (cid:40) χ i for x ∈ Ω i , χ m otherwise. (2)Eq. (1) is completed by the far-field boundary conditions, cf. [5, Eq. (14)], θ ( x ) = H T x for (cid:107) x (cid:107) → ∞ , (3)with H ∈ R denoting the overall (macroscopic) temperature gradient. Due to linearity of the problem,we can introduce the temperature gradient concentration factor A ∈ R × in the form: h ( x ) = A ( x ) H for x ∈ R . (4)As shown first by Hatta and Taya [11], the concentration factor is constant inside the inclusion and admitsthe expression A − ( x ) = (cid:0) A i (cid:1) − = I − S ( χ m ) − (cid:0) χ m − χ i (cid:1) for x ∈ Ω i , (5)where I denotes the unit matrix and S ∈ R × is the Eshelby-like tensor which depends only on thematrix conductivity χ m and the ratios of semi-axes lengths β = a : a and β = a : a , see alsoAppendix 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 A i = A isph I with A isph = 3 χ m χ m + χ i . (6) In the next step, we adopt the results of the previous section to estimate the overall behavior of acomposite material consisting of distinct phases r = 0 , , . . . , N . The value r = 0 is reserved for thematrix phase Ω m and the r -th phase ( r > ) corresponds to the ellipsoidal inclusion Ω ( r ) , characterizedby its semi-axes a ( r )1 , a ( r )2 and a ( r )3 , volume fraction c ( r ) and conductivity χ ( r ) . Following Benveniste’sreformulation [12] of the original Mori-Tanaka scheme [13], the interaction among phases is approxi-mated 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 , , . . . , N. (7)Here, H ( r ) = 1 | Ω ( r ) | (cid:90) Ω ( r ) h ( x ) d x , (8)and T ( r ) is the × partial temperature gradient concentration factor of the r -th phase, given by T ( r ) = (cid:40) I for r = 0 , R ( r ) T ( r ) (cid:96) ( R ( r ) ) T for r = 1 , , . . . , N, (9)ersion August 8, 2011 submitted to Micromachines where the × rotation matrix R ( r ) accounts for the difference in the global and local coordinatesystems, see Section 2.3 for additional details, and T ( r ) (cid:96) equals to A i in Eq. (5) with χ i = χ ( r ) and S determined from values a ( r )1 , a ( r )2 and a ( r )3 .The volume consistency of the overall temperature gradient H and local averages H ( r ) requires H = N (cid:88) r =0 c ( r ) H ( r ) = (cid:32) N (cid:88) r =0 c ( r ) T ( r ) (cid:33) H m . (10)Inverting the above equation gives the average temperature gradient in the matrix as H m = (cid:32) N (cid:88) r =0 c ( r ) T ( r ) (cid:33) − H = A m H , (11)which, when substituted into Eq. (7), yields the explicit expression for the phase temperature gradientsin the form H ( r ) = T ( r ) (cid:32) N (cid:88) r =0 c ( r ) T ( r ) (cid:33) − H = A ( r ) H , (12)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 Q ( r ) = 1 | Ω ( r ) | (cid:90) Ω ( r ) q ( x ) d x for r = 0 , , . . . , N, (13)equals to Q ( r ) = − χ ( r ) H ( r ) for r = 0 , , . . . , N. (14)This allows us to express the macroscopic heat flux in the form Q = N (cid:88) r =0 c ( r ) Q ( r ) = − (cid:32) N (cid:88) r =0 c ( r ) χ ( r ) T ( r ) (cid:33) (cid:32) N (cid:88) r =0 c ( r ) T ( r ) (cid:33) − H , (15)from which we obtain the effective conductivity Q = − χ H H in its final form χ H = (cid:32) c m χ m + N (cid:88) r =1 c ( r ) χ ( r ) T ( r ) (cid:33) (cid:32) c m I + N (cid:88) r =1 c ( r ) T ( r ) (cid:33) − . (16)Finally, assuming that the composite consists of isotropic matrix with conductivity χ m and sphericalisotropic inclusions with conductivities χ ( r ) , Eq. (16) becomes χ H = χ H I , with χ H = c m χ m + N (cid:88) r =1 c ( r ) χ ( r ) T ( r )sph c m + N (cid:88) r =1 c ( r ) T ( r )sph and T ( r )sph = 3 χ m χ m + χ ( s ) . (17)ersion August 8, 2011 submitted to Micromachines
We are now in a position to provide estimates of the effective thermal conductivity for compositeswith M (with M (cid:28) N ) inclusion classes indexed by s = 1 , , . . . , M . Each class is characterized bya single Eshelby-like matrix S in Eq. (5) and represents the reference ellipsoidal inclusion randomlyoriented over the unit hemisphere with an independent uniform distribution of orientation angles.To this goal, consider a quantity X (cid:96) ∈ R × , expressed in a local coordinate system aligned with acertain reference inclusion. Its form in the global coordinate system follows from X ( ϕ, φ, ψ ) = R ( ϕ, φ, ψ ) X (cid:96) R T ( ϕ, φ, ψ ) , (18)where ϕ, φ and ψ denote the Euler angles and the transformation matrix R is provided by R ( ϕ, φ, ψ ) = cos ψ − sin ψ ψ cos ψ
00 0 1 cos φ − sin φ φ φ cos ϕ − sin ϕ ϕ cos ϕ
00 0 1 . (19)The orientation average of X is denoted by double angular brackets: (cid:104)(cid:104) X (cid:105)(cid:105) = 18 π (cid:90) π (cid:90) π (cid:90) π X ( ϕ, φ, ψ ) sin φ d ϕ d φ d ψ. (20)Straightforward calculation, presented e.g. in [14, Appendix A.2.3], reveals that the orientation averag-ing of an arbitrary X (cid:96) ∈ R × yields (cid:104)(cid:104) X (cid:105)(cid:105) = (cid:104)(cid:104) X (cid:105)(cid:105) I with (cid:104)(cid:104) X (cid:105)(cid:105) = 13 (cid:88) i =1 ( X (cid:96) ) ii . (21)Repeating the steps of the previous section with partial temperature gradient concentration factorsreplaced with their orientation averages, we obtain, after some manipulations presented e.g. in [15,Section B], the scalar homogenized conductivity in the form χ H = c m χ m + M (cid:88) s =1 c ( s ) (cid:104)(cid:104) χ ( s ) T ( s ) (cid:105)(cid:105) c m + M (cid:88) s =1 c ( s ) (cid:104)(cid:104) T ( s ) (cid:105)(cid:105) . (22)Therefore, it follows from the comparison of Eq. (22) with Eq. (17) that the system of randomly orientedinclusions embedded in an isotropic matrix is indistinguishable, from the point of view of homogenizedconductivity, from the system of spherical inclusions with an apparent conductivity (cid:101) χ ( s ) = (cid:104)(cid:104) χ ( s ) T ( s ) (cid:105)(cid:105)(cid:104)(cid:104) T ( s ) (cid:105)(cid:105) , (23) Note that so-called ” x convention” is used, in which a conversion into a new coordinates system follows three consecu-tive steps. First, the rotation of angle ϕ around the original X axis is done. Then, the rotation of angle φ around the new x axis is followed by the rotation of angle ψ around the new x axis to finish the conversion. ersion August 8, 2011 submitted to Micromachines which yields χ H = c m χ m + M (cid:88) s =1 c ( s ) (cid:101) χ ( s ) (cid:101) T ( s )sph c m + M (cid:88) s =1 c ( s ) (cid:101) T ( s )sph with (cid:101) T ( s )sph = 3 χ m χ m + (cid:101) χ ( s ) . (24) The presence of imperfect thermal contact at the matrix-inclusion interface ∂ Ω i results in temperaturejump (cid:74) θ (cid:75) , the magnitude of which is provided by Newton’s law, e.g. [16, Section 1.3]: n T ( x ) q ( x ) = k ( x ) (cid:74) θ ( x ) (cid:75) for x ∈ ∂ Ω i , (25)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. Thisrelation, together with Eqs. (1)–(3), defines the single inclusion problem accounting for the presence ofimperfect interface. Its solution is, however, substantially more involved as the temperature gradient in-side an ellipsoidal inclusion becomes position-dependent; the concentration factor is then available onlyin the form of complicated infinite series expansion for spheroidal inclusions [5] or ellipsoidal coatedinclusions [10]. Nevertheless, when restricting the attention to spherical inclusion of radius a , it can beshown that the temperature gradient within inclusion recovers the constant value and the concentrationfactor becomes [5] (cid:98) A sph ( χ i , k, a ) = 3 χ m χ m + (cid:98) χ i ( χ i , k, a ) where (cid:98) χ i ( χ i , k, a ) = χ i akak + χ i , (26)see also Appendix B for simple derivation of this result. Hence, analogously to the previous section, theinterfacial effect can be modeled by replacing the “true” conductivity χ i by the size-dependent apparentvalue (cid:98) χ i provided by Eq. (26) . Assuming in addition that each class of inclusions is characterized byidentical semi-axes lengths a ( s )1 , a ( s )2 and a ( s )3 and interfacial conductance k ( s ) , we propose to extend therelation (24) into the form, cf. [7, Section 2] χ H = c m χ m + M (cid:88) s =1 c ( s ) (cid:98) χ ( s ) (cid:98) T ( s )sph c m + M (cid:88) s =1 c ( s ) (cid:98) T ( s )sph , (cid:98) T ( s )sph = 3 χ m χ m + (cid:98) χ ( s ) , (cid:98) χ ( s ) = (cid:101) χ ( s ) a ( s ) k ( s ) a ( s ) k ( s ) + (cid:101) χ ( s ) , (27)with a ( s ) = (cid:113) a ( s )1 a ( s )2 a ( s )3 and the apparent conductivity (cid:101) χ ( s ) given by Eq. (23). Even though Eq. (27) is applicable to very general material systems, in practice we typically assumesingle inclusion family, with the particle size distribution characterized by a probability density function p ( a ) satisfying p ( a ) ≥ for − ∞ < a < ∞ , (cid:90) ∞−∞ p ( a ) d a = 1 . (28)ersion August 8, 2011 submitted to Micromachines
In this context, the effective conductivity finally becomes χ H = c m χ m + c (1) (cid:110)(cid:98) χ (1) (cid:98) T (1)sph (cid:111) c m + c (1) (cid:110) (cid:98) T (1)sph (cid:111) , (29)where, for an arbitrary function g ( a ) , { g } denotes its expected value given by { g } = (cid:90) ∞−∞ g ( a ) p ( a ) d a. (30)Following e.g. [7,17], the log-normal distribution with the probability density function p ( a ) = 1 √ πaσ exp (cid:32) − (cid:20) ln( a ) − µ √ σ (cid:21) (cid:33) , a > , (31)will be employed to characterize materials’ polydispersity. The parameters µ and σ are provided by [7,Eq. (17)] µ = ln( a ) , σ = 11 . (cid:32) S + √ S + 44 (cid:33) , S = a − a a , (32)where a x denotes the x -th percentile of the particle radii and S is the span of the size distribution.
3. Mori-Tanaka estimates - example results
Two particular examples of real engineering materials are examined in this section to show applica-bility of Eq. (27) and its extension for polydisperse distribution of heterogeneities, Eq. (29), even whendisregarding their actual shape and simply accepting a spherical representation of the inclusions in theMori-Tanaka predictions. The results provided by these two equations are corroborated by availableexperimental data.Random dispersion of copper particles in the Epoxy matrixIn this first example we compare the single-phase Mori-Tanaka predictions with the experimentalresults of de Araujo and Rosenberg [18] who measured the effective thermal conductivities of systemsconsisting of random dispersion of metal particles in the epoxy matrix for several values of the interfacialresistance arising due to acoustic mismatch at the particle-matrix interface particularly for low temper-atures below K. To enhance credibility of the Mori-Tanaka predictions we focus on one particularsystem made from epoxy resin filled with copper particles also examined in [8] in view of the three-point lower bound assuming a random array of superconducting hard spheres (i.e. χ (1) → ∞ ). It hasbeen shown, see [18,19], that for metal-filled composites with ratio α = χ (1) /χ m > the macroscopicconductivity does not depend on the thermal conductivity of the metallic filler, but only on its volumefraction together with the properties of matrix and matrix-particle interface. This becomes evident whenrewriting Eq. (27) in terms of αχ H χ m = 1 + 3 c (1) [ α − (1 + R )] c m [ α + 2 (1 + R )] + 3 c (1) (1 + R ) , (33)ersion August 8, 2011 submitted to Micromachines
10 of 21 where we introduced the dimensionless quantity R adopted in [8] R = χ (1) ka = αχ m ka , (34)where a is the sphere radius. It follows from Eq. (33) that for the inclusion size equal to the Kapitzaradius [20] a K = αα − χ m k ≈ χ m k (35)the effective conductivity equals to that of the matrix, thus the effect of inclusions becomes completelyshielded by the interface. For a < a K , the overall properties are dominated by interfaces and the effectiveconductivity decreases with increasing c (1) even when α > . For a > a K , the bulk properties of phasesbecome dominant; see also [20] for further discussion.In addition to experimental measurements, we also present a comparison with the Torquato-Rintoulthree-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 inter-facial properties can be estimated by measuring the ratio of temperature jump to the applied heat fluxacross a thin bi-material layer, by acoustic mismatch model [4] or, indirectly, from an inverse approachas discussed below. The results are presented for two different temperatures θ = 4 K ( a = 14 . a K ) and θ = 3 K ( a = 4 . a K ). Figure 4.
Evolution of the normalized effective thermal conductivity χ H /χ m as a functionof a) phase contrast α and relative particle radius a/a K ( c (1) = 30% ) and b) volume fraction c (1) of copper particles, c) correlation of measured and calculated values; ρ is the correlationcoefficient. a/a K χ H / χ m α = 10 α = 50 α = 100 α = 1,000 α = 10,000 c (1) [%] χ H / χ m Measured value C a l c u l a t ed v a l ue ρ = 0.999 (a) (b) (c)The influence of parameter α together with expected particle size dependence, now hidden in param-eter R through Eq. (34) , is evident from Figure 4(a) plotted for θ = 4 K ( a K = 3 . µ m). Clearly, theMori-Tanaka predictions confirm negligible influence of χ (1) observed for particulate composites with χ (1) (cid:29) χ m ( α > ) as well as decreasing trend for χ H with decreasing particle size caused by imper-fect thermal contact. Figure 4(b) then displays evolution of normalized effective thermal conductivityas a function of volume fraction of copper particles. Predictive capability of the Mori-Tanaka methodis 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) andersion August 8, 2011 submitted to Micromachines
11 of 21 experimental results. The solid line was derived by linear regression of measured and calculated ef-fective conductivities using the least square method. Another indication of the quality of numericalpredictions can be presented in the form of Pearson’s correlation coefficient written as ρ = n (cid:104) xy (cid:105) − (cid:104) x (cid:105)(cid:104) y (cid:105) (cid:112) n (cid:104) x (cid:105) − (cid:104) x (cid:105) (cid:112) n (cid:104) y (cid:105) − (cid:104) y (cid:105) , (36)where (cid:104) a (cid:105) = (cid:80) ni =1 a i , n is the number of measurements, x stands for the experimental and y for thecorresponding theoretical values. Note that for ρ = 1 the correspondence is exact. In this case thecorrelation coefficient equals to 0.999 suggesting almost perfect match between measured and predictedvalues, also evident from graphical presentation.Al/SiC compositeIn [4] the authors studied the effect of imperfect thermal contact on the macroscopic response ofAl/SiC porous composites. The paper presents the results of a thorough experimental investigation andtraces of an inverse approach in material mechanics for inferring material properties of unknown com-ponents of the composite by matching numerical and experimental results. This approach was firstexploited 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 distributionof spherical particles of identical size to estimate the particle thermal conductivity and interfacial ther-mal conductance for pore-free specimens and subsequently utilized in the two-step application of theHasselman-Johnson model to address the influence of pores. Note that the material data used in thesepredictions 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. Inaddition, we take advantage of available characteristics of the SiC particles, the span S and the 50-thpercentile a , to extend the analysis to polydisperse systems as presented in Section 2.5. The inputmaterial data are listed in Table 1. Table 1.
Material properties [4].Al matrix SiC particles Interface[ Wm − K − ] [ Wm − K − ]187 252.5 72.5 × The Mori-Tanaka predictions stored in Table 2 were provided by Eq. (29). The integral (30) wasevaluated such that the entire interval was split into 1,000 segments, thus considering 1,000 differentparticle sizes of spherical shape. Within each segment s the probability function p ( a ) was approximatedby 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. Standardtrapezoidal integration rule was then used to sum over all 1,000 segments. Examples of probabilitydistribution 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 sensitivityof the effective thermal conductivity to the mean particle size distribution. Note that individual pointsin Figure 5(b) correspond to slightly different volume fraction, see Table 2.ersion August 8, 2011 submitted to Micromachines
12 of 21Table 2.
Characteristics of SiC particles (cf. [4, Table 1]), and comparison of effectivethermal conductivity χ H [ Wm − K − ] between experimental measurements [4] and MT pre-dictions Eq. (29). Sample Radius [ µ m] SiC ResultsNo. a a S vol. Exp. MT1 55 114.5 0.71 0.58 219 217.82 23 65.5 1.02 0.58 210 212.33 19.5 37.5 0.66 0.60 208 208.54 11.5 25 0.79 0.59 198 199.95 7 17 0.86 0.58 195 190.86 5 12 0.82 0.55 184 182.57 2.4 7 1.05 0.53 160 161.3 Figure 5. a) Examples of probability distribution functions of particle radii a for specimensNo. 1 and 7, b) evolution of effective thermal conductivity χ H as a function of particle radius a , c) correlation of measured and calculated values; ρ denotes the correlation coefficient. ln( a ) p ( a ) -2 0 2 4 6 800.20.40.60.81 No. 1No. 7 Median particle radius a (50) [ µ m] χ H [ W / m K ] Measured value [W/mK] C a l c u l a t ed v a l ue [ W / m K ]
160 180 200 220160170180190200210220 Exp [4]Regression ρ = 0.993 (a) (b) (c)Almost negligible deviation from experimental results measured as E = (cid:118)(cid:117)(cid:117)(cid:117)(cid:117)(cid:117)(cid:116) (cid:88) i ( χ exp − χ MT ) (cid:88) i ( χ MT ) = 1 . is, however, not surprising owing to the used material parameters, which were not measured but ratherfitted to a micromechanical model. Even when comparing the Pearson correlation coefficients, . forthe Hasselman-Johnson model and . for polydisperse MT model, the improvement when accountingfor more accurate representation of particle size distribution is marginal. This can be explained by avery small variance associated with adopted distributions, recall Figure 5(a). Nevertheless, it is fair topoint out that unlike the Hasselman-Johnson model the Mori-Tanaka approach is not limited to sphericalparticles providing the transformations given by Eqs. (24) and (27) are admissible. The influence ofshape of particles on the macroscopic response has been put forward, e.g. by J¨ackel in [19], and isersion August 8, 2011 submitted to Micromachines
13 of 21 numerically investigated in the next section suggesting increasing thermal conductivity of the compositewith transition from spherical to needle-like particles, the trend also observed experimentally [18,19].
It has been argued in the previous sections that even very limited information about microstructureamounted to phase properties and corresponding volume fractions might be sufficient to provide a rea-sonable estimate of macroscopic response of various engineering material systems generally classifiedas being macroscopically isotropic. This naturally invites the assumption of spherical representation ofotherwise irregular heterogeneities. Although supported by several practical examples discussed in theprevious 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 influenceof shape and size of inclusions or mismatch of phase material properties on the predicted macroscopicresponse.All numerical results reported bellow are obtained in the two-dimensional setting, hence the argu-ments presented in Section 2 need to be translated to the planar case. In particular, inclusions are mod-eled as elliptic cylinders of aspect ratio β = a /a with a → ∞ . The corresponding Eshelby-liketensor is given by Eq. (42), which leads to the two-dimensional thermal concentration factor for circularinclusion in the form A i = A icir I (2 × with A icir = 2 χ m χ m + χ i . (37)The apparent conductivity due to random orientation is provided by (cid:101) χ ( s ) = (cid:104)(cid:104) χ ( s ) T ( s ) (cid:105)(cid:105) (cid:104)(cid:104) T ( s ) (cid:105)(cid:105) . (38)where (cid:104)(cid:104)•(cid:105)(cid:105) 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 asin the three-dimensional setting, compare Eq. (26) with Eq. (50). Figure 6.
Examples of random macroscopically isotropic microstructures: a) circular cylin-ders ( β = 1 ), b) elliptical cylinders with aspect ratio β = 3 , c) elliptical cylinders withaspect ratio β = 9 .(a) (b) (c)ersion August 8, 2011 submitted to Micromachines
14 of 21
Three particular representatives, generated such as to approximately resemble the real microstruc-tures in Figure 1, appear in Figure 6. To comply with general assumptions put forward in the previoussections, 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 homogenizationstrategy, see e.g. [22,23], to provide estimates of the macroscopic response. The results are plottedin Figure 7.
Figure 7.
Variation of the normalized effective conductivity χ H /χ m for three microstructuresin Figure 6 with perfect interfaces, determined by periodic FEM homogenization for phasecontrast a) α = 3 , b) α = 10 and c) α = 20 ; β denotes the inclusion aspect ratio and c (1) the inclusion volume fraction. c (1) [%] χ H / χ m
20 25 30 35 401.31.41.51.6 β =1 β =3 β =9 c (1) [%] χ H / χ m
20 25 30 35 4011.522.53 β =1 β =3 β =9 c (1) [%] χ H / χ m
20 25 30 35 4011.522.533.5 β =1 β =3 β =9 (a) (b) (c)These results clearly indicate not only the influence of the shape of inclusions on the macroscopicresponse but also a strong dependence of these predictions on the contrast of material properties of indi-vidual phases. Thus drawing from the plots presented in Figure 7(a) one may suggest that the proposedcircular representation of generally non-circular heterogeneities is still acceptable when their shapes onlymoderately deviate from a circle and when the mismatch of phase properties is not too severe, which cer-tainly is the case of a number of real materials as demonstrated in the previous section. This expectationis quite important particularly when dealing with imperfect thermal contact in which case only sphericaland circular inclusions can be easily handled analytically.If the circular approximation of heterogeneities is no longer acceptable or the contrast of phase prop-erties is excessive, one needs to look for more details about microstructure. In such a case, even two-dimensional 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 hetero-geneities, see e.g. [24]. Then, being given the elliptical shape of the inclusion allows us to appropriatelymodify 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 k corr = (cid:101) T (1)cir T (1)cir , (39)now determined for circular inclusions.ersion August 8, 2011 submitted to Micromachines
15 of 21Figure 8. a) Evolution of correction factor k corr for systems with perfect interfaces as a func-tion of the ratio of semi-axes β and variation of normalized effective conductivity χ H /χ m obtained by FEM and MT with modified conductivity (cid:101) χ (1) for phase contrast b) α = 3 andc) α = 20 . β k c o rr α =3 α =10 α =20 c (1) [%] χ H / χ m
20 25 30 35 401.31.41.51.6 FEM β =1FEM β =9MT β =1MT β =9 c (1) [%] χ H / χ m
20 25 30 35 4011.522.533.5 FEM β =1FEM β =9MT β =1MT β =9 (a) (b) (c)Variation of this parameter as a function of the shape of inclusion is seen in Figure 8(a) further con-firming quite strong influence of the phase properties mismatch. The modified conductivity when intro-duced successively into Eq. (37) and Eq. (24) then renders the estimate of effective conductivity almostidentical 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 theMori-Tanaka method is indeed well suited for statistically isotropic random microstructures and secondto promote applicability of this simple transformation from elliptical to circular representations even forshapes markedly distinct from circles. Small but evident deviation of the results observed in Figure 6for the third needle-like microstructure and large mismatch of conductivities of the inclusion and ma-trix equal to can be attributed to the finite size, although infinite in the sense of periodicity, of therepresentative model not large enough to yield statistically isotropic microstructure. Figure 9.
Statistically isotropic distribution of circular cylinders with variable radius of theircross-section. (a) (b) (c) Note that the parameter k corr is determined for two dimensional systems and thus is not applicable to ellipsoidal inclusionsof the same aspect ratio. ersion August 8, 2011 submitted to Micromachines
16 of 21
The second set of numerical simulations addresses the theoretically derived dependence of macro-scopic predictions on the size of inclusions in cases with imperfect interfaces generating jumps in thelocal temperature field. For simplicity, we limit our attention to statistically isotropic distributions ofcircular cylinders with a radius varying from sample to sample. Three such microstructures are shownin Figure 9. Note that the same volume fraction and the same number of inclusions was maintained in allsimulations. 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 im-perfect interface ( k = 60 × Wm − K − and χ (1) = 100 Wm − K − , χ m = 10 Wm − K − and c (1) = 40% ) as a function of inclusion radius a and b) variation of effective conductivityfor systems weakened by cylindrical voids as a function of inclusion volume fraction c (1) ; β denotes inclusion aspect ratio. a [ µ m] χ H / χ m
20 40 60 80 100 1201.11.21.31.41.51.61.71.8 MTFEM β =9 β =3 β =1 c (1) [%] χ H / χ m
20 25 30 35 400.10.20.30.40.50.60.7 β =1 β =3 β =9 (a) (b)The relevant results appear in Figure 10(a). Both the results found from finite element simulations andcorresponding 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 theresults generated for elliptical microstructures from Figure 6 with cross-sectional area equal to the areaof the circle. These are indicated by x-symbol and the ratio of semi-axes of the corresponding ellipticalcylinder. The Mori-Tanaka estimates found from the application of Eqs. (38) and (26) are reasonablyclose further supporting the proposed approach for the modeling of real materials with an imperfectthermal contact. Intuitively, it can be expected that the value of interfacial conductance k may alsoshow some effect as to the estimates of effective conductivities for non-circular inclusions. This notionis supported by the results presented in Figure 10(b) showing variation of effective conductivity of anisotropic matrix weakened by voids with a very low conductivity. Clearly, the influence of shape of theinclusions is quite pronounced.
4. Conclusions
The Mori-Tanaka micromechanical model has been often the primary choice among engineers toprovide quick estimates of the macroscopic response of generally random composites. Motivated by anearly theoretical as well as experimental works on this subject, the Mori-Tanaka method was examinedhere in the light of the solution of a linear steady state heat conduction problem allowing us to estimateersion August 8, 2011 submitted to
Micromachines
17 of 21 the effective thermal conductivity of variety of real engineering materials experiencing an imperfectthermal contact along the constituents interfaces.Adhering to the only limitation, an assumption of macroscopically isotropic composite, it was shownthat the method originally proposed by B ¨ohm and Nogales [7] for a spherical representation of particlesstill applies even to non-spherical particles providing their shape can be suitably quantified, e.g. by anellipsoidal inclusion. In this particular case the Mori-Tanaka predictions were partially corroborated bytwo-dimensional numerical simulations confirming experimentally observed considerable sensitivity ofmacroscopic conductivities to the shape of particles.The fact that for composites with imperfect thermal contact the macroscopic predictions depend onparticle size can be effectively handled by introducing the particle size probability density function di-rectly 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 stan-dard deviation of particle sizes from its mean value. This is particularly appealing, since grading curvesare 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 handledvery effectively with simple micromechanical models with no need for laborious finite element simula-tions of certain representative volumes of real microstructures.
Acknowledgments
The authors are thankful to two anonymous referees for their constructive remarks and sugges-tions on the original version of the manuscript. The financial support provided by the GA ˇCR grantsNo. 106/08/1379 and P105/11/0224 and partially also by the research project CEZ MSM 6840770003 isgratefully acknowledged.
References
1. Torquato, S.
Random heterogeneous materials: Microstructure and macroscopic properties . Springer-Verlag, 2002.2. Benazzouk, A.; Douzane, O.; Mezreb, K.; Laidoudi, B.; Quneudek, M. Thermal conductivity ofcement composites containing rubber waste particles: Experimental study and modelling.
Con-struction and Building Materials , , 573–579.3. Hasselman, D.P.H.; Johnson, L.F. Effective Thermal Conductivity of Composites with InterfacialThermal Barrier Resistance. Journal of Composite Materials , , 508–515.4. Molina, J.M.; Prieto, R.; Narciso, J.; Louis, E. The effect of porosity on the thermal conductivity ofAl-12 wt.% Si/SiC composites. Scripta Materialia , , 582–585.5. Benveniste, Y.; Miloh, T. The Effective Conductivity Of Composites with Imperfect Thermal Con-tact at Constituents Interfaces. International Journal of Engineering Science , , 1537–1552.6. Benveniste, Y. On the effective thermal conductivity of multiphase composites. Journal of AppliedMathematics and Physics , , 696–713.ersion August 8, 2011 submitted to Micromachines
18 of 21
7. B¨ohm, H.J.; Nogales, S. Mori-Tanaka models for the thermal conductivity of composites with inter-facial resistance and particle size distribution.
Composites Science and Technology , , 1181–1187.8. Torquato, S.; Rintoul, M.D. Effect of the Interface on the Properties of Composite Media. PhysicalReview Letters , , 4067–4070.9. Hashin, Z. Thin interphase/imperfect interface in conduction. Journal of Applied Physics , , 2261–2267.10. Dunn, M.L.; Taya, M. The effective thermal conductivity of composites with coated reinforcementand the application to imperfect interfaces. Journal of Applied Physics , , 1711–1722.11. Hatta, H.; Taya, M. Equivalent inclusion method for steady state heat conduction in composites. International Journal of Engineering Science , , 1159–1170.12. Benveniste, Y. A new approach to the application of Mori-Tanaka theory in composite materials. Mechanics of Materials , , 147–157.13. Mori, T.; Tanaka, K. Average stress in matrix and average elastic energy of materials with misfittinginclusions. Acta Metallurgica , , 571–574.14. Str´ansk´y, J. Micromechnical models for thermal conductivity of composite materials with imperfectinterface. Bachelor thesis, Faculty of Civil Engineering, Czech Technical University in Prague,2009. 57 pages (in Czech), http://mech.fsv.cvut.cz/ ∼ zemanj/teaching/theses/09 stransky.pdf.15. Benveniste, Y.; Chen, T.; Dvorak, G.J. The effective thermal conductivity of composites reinforcedby coated cylindrically orthotropic fibers. Journal of Applied Physics , , 2878–2884.16. Lienhard IV, J.H.; Lienhard V, J.H. A heat transfer textbook - 3rd edition . Phlogiston press, Cam-bridge, Massachusetts, 2008.17. Molina, J.; Pinero, E.; Narciso, J.; Garciacordovilla, C.; Louis, E. Liquid metal infiltration intoceramic particle preforms with bimodal size distributions.
Current Opinion in Solid State and Ma-terials Science , , 202–210.18. de Araujo, F.; Rosenberg, H. The thermal conductivity of epoxy-resin/metal-powder compositematerials from 1.7 to 300K. Journal of Physics D: Applied Physics , , 665–675.19. J¨ackel, M. Thermal properties of polymer/particle composites at low temperatures. Cryogenics , , 713–716.20. A.G. Every, A.; Tzou, Y.; Hasselman, D.; Raj, R. The effect of particle size on the thermal conduc-tivity of ZnS/diamond composites. Acta Metallurgica et Materialia , , 123–129.21. Miller, C.A., T.S. Effective conductivity of hard-sphere dispersions. Journal of Applied Physics , , 5486–5493.22. Michel, J.C.; Moulinec, H.; Suquet, P. Effective properties of composite materials with periodicmicrostructure: A computational approach. Computer Methods in Applied Mechanics and Engi-neering , , 109–143.23. Zeman, J.; ˇSejnoha, M. Numerical evaluation of effective properties of graphite fiber tow impreg-nated by polymer matrix. Journal of the Mechanics and Physics of Solids , , 69–90.24. Tsukrov, I.; Piat, R.; Novak, J.; Schnack, E. Micromechanical Modeling of Porous Carbon/CarbonComposites. Mechanics of Advanced Materials and Structures , , 43–54.ersion August 8, 2011 submitted to Micromachines
19 of 21
25. Chen, T.; Yang, S.H. The problem of thermal conduction for two ellipsoidal inhomogeneities in ananisotropic medium and its relevance to composite materials.
Acta Mechanica , , 41–58. A. Eshelby-like tensor
The Eshelby-like tensor for the solution of thermal conductivity problem was introduced by Hattaand Taya in [11]. For an ellipsoidal inclusion with semi-axes a , a , a found in an isotropic matrix itreceives the form S ij = a a a ∂∂x i ∂x j (cid:90) ∞ (cid:18) x a + s + x a + s + x a + s (cid:19) s d s, (39) ∆ s = (cid:113) ( a + s ) ( a + s ) ( a + s ) . (40)Closed form solutions of integral (39) for some special cases of ellipsoidal shapes of the inclusion canbe found in [11]. For a general ellipsoid the solution was introduced by Chen and Yang in [25]. Forcircular and spherical shapes needed in the present study the S tensor reads • Sphere ( a = a = a ) S = S = S = 13 and S ij = 0 for i (cid:54) = j. (41) • Elliptic cylinder ( a → ∞ ) S = a a + a , S = a a + a , S = 0 and S ij = 0 for i (cid:54) = j. (42) B. Single spherical homogeneity with imperfect interface
This section outlines derivation of the replacement conductivity (cid:98) χ i and the concentration factor (cid:98) A sph introduced in Eq. (26). Since used in numerical calculations in Section 3.2, its two dimensional formatis presented as well. It it shown that both 2D and 3D concentration factors can be recovered from thesolution of a 1D problem using a simple geometrical argument. Figure 11.
Imperfect interface and temperature progress for: a) 1D, b) 2D. k k c m c m c i c i realitysubstitution ║ q ║ L n j c m c i k x (a) (b)To that end, consider one-dimensional heat conduction problem depicted in Figure 11(a). Assum-ing imperfect thermal contact, the temperature drop across an infinitely thin interface layer is given byersion August 8, 2011 submitted to Micromachines
20 of 21
Eq. (25). The local temperature gradient for perfect interface between a solitary inclusion embedded intoan infinite matrix follows from Eq. (4) H i = A i H = χ m χ i H. (43)To arrive at similar format of Eq. (43) for imperfect contact, we imagine the interface temperature jumpbeing smeared over the inclusion. Since the heat flux Q associated with the macroscopic temperaturegradient H is constant throughout the composite, we obtain the total temperature change in the substituteinclusion in the form ∆ (cid:98) θ i = ∆ θ i + 2 (cid:74) θ (cid:75) = − Q (cid:18) Lχ i + 2 k (cid:19) , (44)where L stands for the inclusion length. Next, defining the local temperature gradient in the substituteinclusion as H i = ∆ (cid:98) θ i /L yields the local constitutive law in terms of the replacement conductivity (cid:98) χ i Q = − (cid:98) χ i H i = − (cid:98) χ i ∆ (cid:98) θ i L = − (cid:98) χ i L (cid:20) − Q (cid:18) Lχ i + 2 k (cid:19)(cid:21) , (45)so that (cid:98) χ i = χ i LkLk + 2 χ i , (46)and in analogy with Eq. (43) we finally get H i = (cid:98) A i H = χ m (cid:98) χ ( i ) H. (47)The two-dimensional problem of a solitary circular inclusion is treated similarly. We build up on theassumption that the temperature gradient in the inclusion is constant and collinear with the prescribed farfield gradient parallel to the local x -axis, see Figure 11(b). To draw similarity with 1D case we dividethe inclusion into parallel filaments with the length L = 2 a cos ϕ . Next, define a unit vector normal tothe inclusion surface n = (cos ϕ, sin ϕ ) T and in analogy to Eq. (44)) write the total temperature changein each filament for the constant heat flux q i = (cid:0) q i , (cid:1) T as ∆ (cid:98) θ i = ∆ θ i + 2 (cid:74) θ (cid:75) = − q i (cid:18) d cos ϕχ i + 2 cos ϕk (cid:19) , (48)where d is the inclusion diameter. The equivalent conductivity has to fulfill the condition q i = − (cid:98) χ i ∆ (cid:98) θ i d cos ϕ = − (cid:98) χ i d cos ϕ (cid:20) − q i (cid:18) a cos ϕχ i + 2 cos ϕk (cid:19)(cid:21) , (49)which yields (cid:98) χ i = χ i akak + χ i . (50)Consequently, the concentration factor of the substitute inclusion attains the form (cid:98) A i = 2 χ m χ m + (cid:98) χ i . (51)ersion August 8, 2011 submitted to Micromachines
21 of 21
The analysis of a spherical inclusion follows identical steps. The replacement thermal conductivityfor constant heat flux q i = (cid:0) q i , , (cid:1) T thus receives the same form as in Eq. (50) rendering the searchedconcentration factor as, recall Eqs. (27), (cid:98) A i = 3 χ m χ m + (cid:98) χ i . (52)c (cid:13) August 8, 2011 by the authors; submitted to