## Abstract

Obtaining an accurate higher order statistical description of heterogeneous materials and using this information to predict effective material behaviour with high fidelity has remained an outstanding problem for many years. In a recent letter, Gillman & Matouš (2014 *Phys. Lett. A* 378, 3070–3073. ()) accurately evaluated the three-point microstructural parameter that arises in third-order theories and predicted with high accuracy the effective thermal conductivity of highly packed material systems. Expanding this work here, we predict for the first time effective thermo-mechanical properties of granular Platonic solid packs using third-order statistical micromechanics. Systems of impenetrable and penetrable spheres are considered to verify adaptive methods for computing *n*-point probability functions directly from three-dimensional microstructures, and excellent agreement is shown with simulation. Moreover, a significant shape effect is discovered for the effective thermal conductivity of highly packed composites, whereas a moderate shape effect is exhibited for the elastic constants.

## 1. Introduction

Quantifying the morphology of many-body systems has applications in many scientific fields at a variety of length scales from molecular configurations [1,2] up to composite microstructures [3,4] and celestial bodies [5–7]. Often, these systems have varying degrees of randomness and can only be described though use of statistics. Moreover, the connection of a statistical description to observed physical phenomena has long been studied [3,8]. Bernal [9] first motivated the importance of higher order statistics when investigating the molecular arrangement of liquids using the radial distribution function. *n*-point probability functions were introduced by Frisch & Stillinger [10] to describe radiation scattering patterns in packs of voided media. Statistical characterization of many-body systems has remained prevalent as evidenced by more recent studies of the background cosmic radiation [6,11] and characterization of molecules in amorphous condensed matter (e.g. glasses [12]). Note that many-body systems are often spatially complex, and accurately obtaining higher order statistical descriptions in three dimensions have proved difficult.

The focus of this work is to predict effective transport and mechanical properties of particulate composites from higher order statistical data. This field has applications in many scientific and engineering disciplines including flow in fractured rock [13], design of structural composites [14] and use of conductive adhesives in electronics [15], just to name a few. Specifically, a primary objective of this work is to predict thermo-mechanical properties of highly filled crystalline packs (Platonic solids) and investigate the effect of the inclusion shape for a wide range of microstructures. Experimental observations for particular systems have shown significant shape effects. For example, Timofeeva *et al.* [16] show that alumina nanofluids have higher thermal conductivities for alumina particles with larger aspect ratios, and linear models dependent on the particle shape were developed. However, application of these linear models is limited to low particle concentrations. Direct numerical modelling of a shape effect was performed by Yvonnet *et al*. [17,18] for both mechanical and thermal properties. However, studies were limited to single particle simulations at low volume fractions. Chauhan *et al.* [19] developed a micromechanics model for the effective conductivity from a chain of single particle unit cells for various shapes. In [19], good agreement was obtained for some experiments at moderate particle volume fractions, but large errors were present for other tests.

Statistical micromechanics theories developed over the past half century have proved accurate for wide ranging particle concentrations. Brown [20] first used *n*-point probability functions to connect the microstructural details to the effective material behaviour. Hashin & Shtrikman [21] derived the most accurate second-order bounds for the permittivity/conductivity and elastic constants [22] for isotropic two-phase materials. However, the sole microstructural information included in these bounds is the volume fractions (one-point probability functions) of the respective constituents. Three-point bounds for the effective permittivity were introduced by Beran [23], and these bounds were later simplified independently by Torquato [24] and Milton [25] as a function of the volume fraction *c*_{q} and a microstructural parameter *ζ*_{q} of each material phase *q*. Beran & Molneux [26], McCoy [27] and Milton & Phan-Thien [28] derived three-point bounds for the effective shear and bulk moduli incorporating *c*_{q}, *ζ*_{q} and an additional microstructural parameter *η*_{q}. These microstructural parameters *ζ*_{q} and *η*_{q} depend on three-dimensional integrals involving one-, two- and three-point probability functions. While additional theories have been introduced incorporating nonlinear material response [29] and imperfect interface behaviour [30], demonstration of these theories has been severely restricted. Strong microstructural assumptions are an impediment to accurately computing *ζ*_{q} and *η*_{q} for complex microstructures. Therefore, application of third-order models has been primarily directed to systems of spheres.

Many have studied systems of overlapping spheres [31–34], where the probability functions can be described by a simple exponential decay function. Numerical methods have also been proposed for computing the *n*-point probability functions directly from microstructures [32,35]. However, application of these methods has been limited to computationally generated morphologies of overlapping spheres [36] or imaging approaches of real systems like sandstone, which closely resemble the overlapping sphere model [37]. More complex configurations of impenetrable spheres were considered by Torquato and co-workers [38–41], where the microstructural parameters *ζ*_{q} and *η*_{q} were computed using a statistical mechanics or a single-body operator concept. These methods unfortunately bypassed the computation of the probability functions and have thus proved inaccurate. Despite years of progress in computing these microstructural parameters, direct computation of the statistical functions from three-dimensional microstructural domains with resolution required in third-order models has been elusive until the recent work of Gillman & Matouš [42].

The adaptive numerical method introduced in [42] to compute statistical information directly from three-dimensional microstructures is designed to extend beyond systems of spheres and allows for investigation of an inclusion shape effect. In this work, we extend the computational framework presented in [42] to the microstructural parameter *η*_{q}, and a unique study is conducted on the effective thermal conductivity and elastic constants for highly filled isotropic systems of crystalline inclusions (Platonic solids). First, we present an overview of the mathematical theory of third-order statistical micromechanics. Next, we describe and rigorously verify the adaptive numerical methods for computing *ζ*_{q} and *η*_{q} in third-order models. For the first time, novel predictions of the effective thermo-mechanical properties for isotropic systems of Platonic solids are obtained.

## 2. Mathematical theory

In this section, a brief summary of higher order statistical micromechanics is presented for computing the effective thermal conductivity and elastic constants of statistically isotropic microstructures.

### (a) *n*-point probability functions

In this work, *n*-point probability functions are used to describe the morphology of a heterogeneous material. First, a phase indicator function of material phase *q* at a position *α* of an ensemble space *p*(*α*) is the probability density function of *α* in *n*-point probability function, *S*_{qs…t}(*x*_{1},*x*_{2},…,*x*_{n}), reads
*q*,*s*,…,*t* existing at points *x*_{1},*x*_{2},…,*x*_{n}, simultaneously.

In this work, materials containing statistically homogeneous (translationally invariant) and isotropic (rotationally invariant) phases are considered. For statistically homogeneous systems, it is useful to define volume averages. When assuming ergodicity of homogeneous systems, ensemble averaging corresponds to volume averaging in the infinite volume limit
** l** is a translation vector, and

*Ω*is the volume of the domain. For ergodic, homogeneous and isotropic systems, the one-, two- and three-point probability function reduce to

*S*

_{q}(

*x*_{1})=

*c*

_{q},

*S*

_{qs}(

*x*_{1},

*x*_{2})=

*S*

_{qs}(

*r*

_{1}=|

*x*_{2}−

*x*_{1}|) and

*S*

_{qst}(

*x*_{1},

*x*_{2},

*x*_{3})=

*S*

_{qst}(

*r*

_{1}=|

*x*_{2}−

*x*_{1}|,

*r*

_{2}=|

*x*_{3}−

*x*_{1}|,

*θ*), where

*θ*is the angle between the vectors

*x*_{2}−

*x*_{1}and

*x*_{3}−

*x*_{1}. Moreover, long range order is also absent for randomly configured microstructures, and the following limits hold for the two-point probability functions

Analytical expressions of the *n*-point probability functions do not generally exist for random morphologies. Therefore, Monte Carlo (MC)-based statistical sampling algorithms are often used. When computing function values of a *n*-point probability function, many random samples or throws are considered and averaged. The accuracy of MC methods is *Stat3D*. The *N*_{r} random samples are decomposed on *N*_{r}, consists of a random translation (described by three position values) and random rotation (described by three angles) of a (*n*–1)-simplex within the three-dimensional microstructure as illustrated in figure 1. For example, computing a function value of the isotropic three-point probability function *S*_{qst}(*r*_{1},*r*_{2},*θ*) requires random placement of a triangle (a two-simplex) described by two side lengths, *r*_{1} and *r*_{2}, and the angle between the sides *θ*. In a similar manner, the two-point probability function *S*_{qs}(*r*_{1}) requires random samples of a line segment (one-simplex) with length *r*_{1}.

Figure 2 shows selected planes of *θ* for the isotropic three-point probability function, *S*_{ppp}(*r*_{1},*r*_{2},*θ*), for a system of impenetrable (non-overlapping) spheres considered by Gillman & Matouš [42]. Note that this pack was generated using a packing algorithm [43] based on the Lubachevsky–Stillinger method [44,45]. Others have computed one-, two- and three-point probability functions for systems of spheres for the purpose of quantifying the configuration [46,47]. However, computing these complex functions for the entire function domain, *Ω*, with the fidelity required by third-order statistical micromechanics models presented in §2b has not been previously accomplished. Therefore, we introduce the adaptive interpolation strategy described in §3 to overcome these difficulties. See [42] for brief information on the adaptive interpolation scheme.

### (b) Third-order models of effective material behaviour

For completeness of presentation, this section summarizes the statistical micromechanics theories. Books by Milton [48] and Torquato [3] present the detailed mathematical theory for determining effective material properties.

#### (i) Thermal conductivity

Regarding effective thermal properties, the steady-state heat conduction problem is considered. At the local scale, the conservation of energy is
** q**(

**) denotes the heat flux at a spatial point**

*x**Ω*. The heat flux is related to the temperature field,

*T*(

**), assuming Fourier’s Law**

*x***(**

*κ***) is the second-order thermal conductivity tensor and**

*x***(**

*Q***)=−∇**

*x**T*(

**). The local thermal conductivity for a two-phase particulate composite is defined as**

*x*Using variational principles, Beran [23] derived third-order bounds for the effective thermal conductivity, *κ*^{L}≤*κ*_{e}≤*κ*^{U}, of homogeneous and isotropic heterogeneous materials. The superscripts L and U denote the lower and upper bounds, respectively. The bounds were later simplified by Milton [25] and Torquato [24] for two-phase materials assuming isotropic material behaviour and are given as
*c*_{q}, the local conductivity, *κ*_{q}, and the microstructural parameter, *ζ*_{q}, of the particle (*q*=*p*) and matrix (*q*=*m*) material phases. The microstructural parameter, *ζ*_{q} is defined as

#### (ii) Elastic constants

When considering the effective mechanical response, the conservation of linear momentum neglecting inertia and body forces reads
** σ**(

**) is the symmetric second-order stress tensor at a point**

*x***(**

*u**x*) is the local displacement field, and

**(**

*L***) is the local elasticity tensor given as**

*x*

*L*_{q}=3

*K*

_{q}

*Λ*_{h}+2

*G*

_{q}

*Λ*_{s}, where

*K*

_{q}and

*G*

_{q}are the bulk and shear moduli of material phase

*q*, and the hydrostatic and shear projection tensors are defined as

**1**and

**are the second and symmetric fourth-order identity tensors.**

*I*Using variational principles, third-order bounds of the effective bulk (*K*^{L}≤*K*_{e}≤*K*^{U}) and shear moduli (*G*^{L}≤*G*_{e}≤*G*^{U}) have been derived by Beran & Molyneux [26], McCoy [27], and were later simplified and improved upon by Milton [25], Milton & Phan-Thien [28] and Gibiansky & Torquato [51]. The bounds of the effective moduli considered in this work are given in Milton & Phan-Thien [28] as
*Ξ*_{L} and *Ξ*_{U} are presented in appendix A. The bounds are functions of *K*_{m}, *K*_{p}, *G*_{m}, *G*_{p}, *c*_{m}, *c*_{p}, *ζ*_{m}, *ζ*_{p}, *η*_{m} and *η*_{p}. Note that in addition to the microstructural parameter *ζ*_{q} that also appeared in the bounds of the effective thermal conductivity, the bounds of the shear modulus are functions of an additional microstructural parameter *η*_{q}. This parameter is defined as

Similar to the TPA presented for the effective thermal conductivity (see equation (2.13)), Torquato [52,53] has more recently derived a TPA for the effective bulk and shear modulus (see appendix A).

## 3. Adaptive interpolation/integration method for computing microstructural parameters

The integrals presented in equations (2.11) and (2.20) for computing the microstructural parameters *ζ*_{q} and *η*_{q} involve the interaction of the one-, two- and three-point probability functions. Moreover, the statistics are multiplied by the term 1/(*r*_{1}*r*_{2}) in the integral kernel. When the values for *r*_{1} and *r*_{2} are small, the error in the statistics values can be substantially amplified. These interactions, especially for systems with complex statistics functions (figure 2) have prevented the accurate evaluation of *ζ*_{q} and *η*_{q} for particulate microstructures until recent work of Gillman & Matouš [42]. In this section, we expand on the numerical methods introduced in [42] for constructing an adaptive interpolant of

Others [46,54] have attempted an MC sampling strategy to compute probability functions, but on regular structured grids. An adaptive triangulation technique is used to overcome inefficiencies of structured grids. To simply explain our adaptive scheme, an illustration is shown for a two-dimensional product peak function introduced by Genz [55] in figure 3. This function is typical for testing multi-dimensional numerical integration schemes and is defined as
*a*_{1}=5, *a*_{2}=10, *w*_{1}=0.25 and *w*_{2}=0.4 for this example.

In order to capture the important features of this function, a Delaunay triangulation is iteratively constructed with *C*^{0} continuity using the CGAL library [56,57]. The triangulation is denoted as *c*). This triangulation and associated function values, *f*(*x*_{1},*x*_{2},…,*x*_{n}), define the initial interpolant *l* is the adaptive level). For a given iteration, a bisection method is used to refine the triangulation based on the local error of the interpolant. For all triangles/tetrahedrons in a given iteration of the interpolant *b*, where the line segment along *x*_{1}=0.5 is considered. The solid black line is the smooth function, *f*_{G}(*x*_{1}=0.5,*x*_{2}), and the dashed black line is the initial interpolant,

In this work, the adaptive triangulation technique is used to create an interpolant of *r*_{1}=*r*_{2}). This triangulation and associated function values, *ζ*_{q} and *η*_{q} in equations (2.11) and (2.20), are computed via simplex integration. MC integration of each tetrahedron is performed [58], and a convergence study determined that using *N*_{int}=1000 random integration points per tetrahedron is sufficient for all particulate systems considered. Given that there are

## 4. Numerical results

Using the methods described in the prior sections, bounds and estimates of the effective thermal conductivity and elastic constants were computed with the well-resolved (WR) statistical description. First, verification studies were conducted on commonly studied microstructures of overlapping and impenetrable spheres. After verification of the methods, results are presented for statically isotropic systems of Platonic solids. Note that throughout this section, we denote our results with WR due to use of WR statistics.

### (a) Verification

In order to assess the accuracy of our numerical method, a verification study was conducted. The interpolation/integration scheme is verified using the analytical form of the overlapping sphere model. Then, the methods are used with MC sampling for systems of impenetrable spheres. The resulting predictions of the effective thermal conductivity and elastic constants are compared to simulation data.

#### (i) Penetrable sphere model

For the purpose of verification, microstructures composed of overlapping spheres are studied. This is one of the few systems, where the *n*-point probability functions can be defined analytically. The probability function for the matrix phase is defined as
*ρ* is the number density of spheres and *V*_{n} is the union volume of *n* spheres with diameter *d*.

Using this analytical expression for the *n*-point probability function, an interpolant is created for *ζ*_{m} and *η*_{m} are presented in table 1. Note that results for *ζ*_{m} (or *ζ*_{p}) were first presented in [42] and are provided here for completeness. In order to assess the accuracy of these methods, the values computed in this work are compared to the most accurate result presented in the literature [34]. Error measures introduced to quantify the accuracy of the microstructural parameters are defined as

#### (ii) Random systems of impenetrable spheres

Following verification of the penetrable sphere model, systems of impenetrable (non-overlapping) monodisperse spheres are considered. Such systems have been studied extensively by Torquato and co-workers [38–40]. However, despite significant progress, computing the microstructural parameters *ζ*_{q} and *η*_{q} for three-dimensional systems directly has been unachievable (due to the complexity of the *n*-point probability functions) until the recent work on thermal conductivity of Gillman & Matouš [42].

The same monodisperse sphere systems presented in [42] for verification of the effective thermal conductivity are considered for verification of the third-order estimates of the effective bulk and shear moduli. These systems were generated using a packing algorithm [43] based on the Lubachevsky–Stillinger method [44,45]. Cubic domains consisting of *d* were generated, and a high expansion rate was prescribed to ensure statistical isotropy. After a convergence study, it was determined that *ζ*_{p} and *η*_{p} with acceptable accuracy. Furthermore, the convergence study uncovered that it is necessary to use *N*_{r}=10^{7} random samples for each *n*-point probability function evaluation in equation (2.12) to satisfy the tolerance tol. The computed microstructural parameters *ζ*_{p} and *η*_{p} are presented in table 2. Comparison is made to approximations given in [40,59]. The improvement of *ζ*_{p} from the result presented in [40], denoted as *ζ*^{R2}_{p}, is quantified by introducing an error metric *ζ*_{p} is shown for lower volume fractions, but significant improvement has been achieved for higher volume fractions. For the parameter *η*_{p}, the results are compared to the work of Torquato *et al.* [59]. The differences from the results in [59] are quantified by the following metric *ε*^{η}_{IS} is non-monotonic with respect to *c*_{p}. However, this is not unexpected given the simple inaccurate relation, *η*_{p}=0.4827*c*_{p}, derived in [59], as *η*_{p} is only first-order accurate with respect to *c*_{p}.

The computed microstructural parameters, *ζ*_{p} and *η*_{p}, are then used within the third-order statistical micromechanics theories for the effective bulk and shear moduli as described in §2b(ii). We use our in-house software package *Prop3D* to compute all thermo-mechanical properties. In the remainder of this work, an infinite contrast ratio between particle and matrix phases is considered (*ν*_{m}=0.25). The WR three-point lower bounds *ζ*^{R2}_{p} and *ζ*_{p} and *η*_{p} (table 2). However, for *c*_{p}=0.6, where improvements in *ζ*_{p} and *η*_{p} were 46.8% and 18.0%, respectively, the lower bounds of *K*_{e} and *G*_{e} are improved by 2.78% and 3.54%.

In figure 4, we show results for the well-resolved three-point approximation (TPA-WR, solid line with circles) of *K*_{e} and *G*_{e} (see equations (A 2) and (A 5)). These results are compared to FEM data from Segurado & Llorca [60] (filled circles), the TPA using Torquato’s approximations [40,59] (TPA, dash-dotted line) and the widely used Hashin–Shtrikman bound (HS, dashed line). Note that the HS bounds (presented in appendix B) are unable to capture details of the morphology other than the volume fraction, *c*_{p}, and is the least accurate prediction. When comparing TPA-WR to Torquato’s TPA, the largest improvement in the estimate of *K*_{e} and *G*_{e} is seen for higher volume fractions. For the volume fraction *c*_{p}=0.6, the estimates of the bulk modulus and shear modulus have been improved by 5.25% and 6.4%, respectively. Note that the improvement in the elastic constants is less pronounced than those of the effective thermal conductivity [42] (17.4% improvement in *κ*_{e} achieved for *c*_{p}=0.6). This smaller improvement shows that the third-order estimates for the elastic constants are less sensitive to the morphology. When comparing the TPA-WR to finite-element method (FEM) data in [60] (filled black circles), all finite-element method (FEM) data points are within 4.3% of the TPA-WR for the bulk modulus (*K*_{e}) and within 9% of the TPA-WR for the shear modulus (*G*_{e}). We point out that FEM data are not error free. Therefore, small differences between our statistical micromechanics and FEM results [60] for *K*_{e} and *G*_{e} must be interpreted more qualitatively. It is also important to note that the FEM data presented in [60] is only reported for particle volume fractions up to *c*_{p}=0.5. Direct numerical modelling methods like FEM can become challenging for highly filled systems due to numerical issues, especially for high material contrast. However, the computational complexity of higher order statistical micromechanics presented in this work does not increase for highly filled packs.

### (b) Random systems of Platonic solids

As mentioned in the Introduction, predicting the effective material behaviour of particulate materials by third-order statistical micromechanics has proved elusive due to numerical difficulties. Therefore, in this seciton, we elucidate the importance of our WR computations with the novel application of these theories to packs of monodisperse Platonic solids. For the first time, the statistical micromechanics is applied to crystalline packs (Platonic solids) with fidelity and resolution needed for meaningful thermo-mechanical engineering analysis.

#### (iii) Generation of microstructures

To study the effect of inclusion shape on the effective material properties, the packing code *Rocpack* [61–63] is used to generate statistically isotropic systems of Platonic solids. The algorithm is a hybrid of the Lubachevsky–Stillinger [44] and Adaptive Shrinking Cell [64,65] packing algorithms, where infinitesimal particles are randomly distributed within a volume and allowed to grow until a desired volume fraction is reached. Collisions are handled to prevent particle intersection and maintain statistical isotropy. In this work, packs with volume fractions of *c*_{p}=0.2, 0.4, 0.5 and 0.6 containing *N*_{p}=10 313, 20 626, 25 783 and 30 940 particles, respectively, are packed in a unit cubic volume. This corresponds to a side length of approximately 30*d*, where *d* is the equivalent sphere diameter (diameter computed from volume of a particle assuming spherical shape).

Geometric quantities associated with the Platonic solids are summarized in table 4. The number of faces (*N*_{F}), vertices (*N*_{V}) and edges (*N*_{E}) are presented along with the dihedral angle (*α*), inradius (*R*_{in}), midradius (*R*_{mid}) and circumradius (*R*_{circ}). Note that the inradius (*R*_{in}) is the radius of an inscribed sphere for the polyhedra, the midradius (*R*_{mid}) is the radius of a sphere intersecting the midpoints of the edges and the circumradius (*R*_{circ}) is the radius of a circumscribed sphere. The inradius and dihedral angle (*α*) increase with the number of faces *N*_{F}. To aid further analysis, the percolation threshold volume fraction, *c*^{per}_{p}, for systems of overlapping shapes from [66] is also presented. Note that *c*^{per}_{p} increases with *N*_{F}, and the percolation results for overlapping systems assist in explaining the structures and their effect on the effective material response.

#### (iv) Statistical analysis

After generating the packs of Platonic solids, the systems are statistically characterized. The isotropic two-point probability function *S*_{pp}(*r*) for *c*_{p}=0.6 packs of tetrahedra, octahedra and icosahedra is shown in figure 5. The functions are compared to the system of spheres considered in §4b(ii). Note that packs with polyhedra that are composed of more faces (*N*_{F}) approach the statistical function for a pack of spheres. At the origin (*r*=0), the slope of the function is steepest for the system of tetraheda. The three-point probability functions for packs with *c*_{p}=0.6 are presented for selected values of *θ* in figure 6. It can be seen that the system of icosahedra closely resembles the function for spheres (cf. figure 2). Meanwhile, the probability function for the systems of tetraheda has less complex features. In general, the packs of polyhedra with smaller *N*_{F} result in less complex *n*-point probability functions as compared with ones with more faces. Note that the one-point and two-point probability functions are contained within the three-point probability function (see equations (2.5) and (2.6)). The function *S*_{ppp}(*r*_{1},*r*_{2},*θ*) degenerates to the volume fraction *c*_{p}=0.6 at the point *r*_{1}=*r*_{2}=0, and the function degenerates to the two-point probability function along the line *r*_{1}=*r*_{2} (*θ*=0) and along the planes *r*_{1}=0 and *r*_{2}=0. These degeneracies are used to evaluate the MC sampling error in the statistical functions.

It is particularly important to understand the error related to the interactions of the probability functions in *c*_{p}=0.2, both of these error metrics are below 1.3×10^{−4} when the number of random samples is *N*_{r}=10^{7}. For all packs with *c*_{p}=0.6, both of these metrics are below 2.7×10^{−4}. These errors are at least 2× less than the tolerance tol in the adaptive interpolation scheme.

Following the statistical analysis and error quantification, the microstructural parameters *ζ*_{p} and *η*_{p} are computed. The same numerical parameters reported for the systems of impenetrable spheres (*N*_{r}=10^{7}) are used. The parameters *ζ*_{p} and *η*_{p} for all packs of Platonic solids are shown in figure 7*a*,*b*, respectively. Note that there is a monotonic increase of *ζ*_{p} as the number of faces, *N*_{F}, decreases. However, the value of *η*_{p} does not follow this simple trend.

Considering the local configurations of the shapes aids in understanding the values obtained for *η*_{p}. A measure *β*_{IJ} is introduced to quantify the local arrangement of particles. *β*_{IJ}, as illustrated in figure 8, is the distance between a particle *I* and a neighbouring particle *J*. The minimum distance between particle *I* and all neighbouring particles is denoted *β*_{I}=min(*β*_{IJ}). Finally, *β*=mean(*β*_{I}) is defined as the average *β*_{I} for all particles in a pack. In figure 9*a*, *β* is presented for each pack and is normalized by the equivalent sphere diameter, *d*. In order to complement the values of *β*/*d*, the distances, *β*_{12}/*d*, associated with the various types of contacts between only two single Platonic solids are shown in figure 9*b*. The distances associated with a vertex–vertex (V–V, *β*_{12}=*R*_{circ}+*R*_{circ}), vertex–edge (V–E, *β*_{12}=*R*_{circ}+*R*_{mid}), vertex–face (V–F, *β*_{12}=*R*_{mid}+*R*_{in}), edge–edge (E–E, *β*_{12}=*R*_{mid}+*R*_{mid}), edge–face (E–F, *β*_{12}=*R*_{mid}+*R*_{in}) and face–face (F–F, *β*_{12}=*R*_{in}+*R*_{in}) contact between two particles are presented. The F–F contact represents the lower limit of *β*/*d* and corresponds to the inscribed radius, *R*_{in} (table 4). *β*_{12}/*d* less than one signifies particles in close proximity and assembling mostly in a E–F of F–F contact. For illustration, note that two tetrahedrons come in V–V contact at larger distances apart and are closest together with a F–F contact. In figure 9*a*, *β*/*d* is less than one at the smallest volume fractions for the system of tetrahedra. Moreover, when the measure *β*/*d* becomes less than one, a significant slope decrease in the function *η*_{p} is observed (figure 7*b*). With tetrahedra for example (dotted line with triangular markers), *β*/*d* becomes less than one between *c*_{p}=0.2−0.3, and a change in slope of *η*_{p} is observed near *c*_{p}=0.2. A similar change in slope can be also observed for *ζ*_{p} as *β*/*d* approaches one (figure 7*a*).

#### (v) Third-order estimates of effective material behaviour

With *ζ*_{p} and *η*_{p} accurately computed, the TPA of the thermal conductivity and bulk and shear moduli are considered. Figure 10 presents the TPA-WR estimates of the thermal conductivity, *κ*_{e}, using equation (2.13) for an infinite contrast ratio (*κ*_{e}/*κ*_{m}) for *c*_{p}=0.6 is 16.81. In comparison to the pack of spheres (*κ*_{e}/*κ*_{m}=8.42), the effective thermal conductivity has increased by 1.99 times (99.6% increase). In general, as the number of faces, *N*_{F}, decreases, the effective conductivity increases. Moreover, when revisiting the local morphology measure *β*, tetrahedra came to proximity frequently through F–F contact (*β*/*d*<1) (figure 9*a*) and are capable to contact at large distances, *β*_{12}/*d*>1, through V–V contact (figure 9*b*). This contact behaviour suggests that tetrahedra are near percolation compared to other packs, which leads to greater effective conductivity. This trend is also consistent with the percolation threshold, *c*^{per}_{p}, for systems of overlapping Platonic solids listed in table 4. This is the first time to the authors’ knowledge that a shape effect has been predicted for isotropic systems with varying shapes up to high packing fractions using third-order statistical micromechanics. The resulting predictions of the thermal conductivity are also compared to the HS lower bound (see equation (B 1)) that is unable to capture any shape effect.

The TPA of the effective bulk modulus (equation (A 2)) and shear modulus (equation (A 5)) are presented in figure 11*a*,*b*, respectively. The same shape effect trend seen for the effective thermal conductivity is exhibited for both the shear and bulk moduli. However, the magnitude is less pronounced. With tetrahedra for example (dotted line with triangular markers), the normalized effective bulk modulus (*K*_{e}/*K*_{m}) for *c*_{p}=0.6 is 5.00. In comparison to the pack of spheres (*K*_{e}/*K*_{m}=4.29), the effective bulk modulus increases by 1.16 times (16.6% increase). Meanwhile, the normalized effective shear modulus (*G*_{e}/*G*_{m}) at *c*_{p}=0.6 for the pack of tetrahedra increases by 1.13 times (12.6% increase) above the system of spheres. This indicates that the effective elastic properties are less sensitive to the inclusion shape. These predictions are also compared to the HS lower bounds (see equations (B 2) and (B 3)).

## 5. Conclusion

In this work, we predict with high accuracy thermal conductivity and elastic constants of isotropic packs of Platonic solids (crystalline materials). Verification studies are conducted for systems of overlapping and hard monodisperse spheres, and numerical approaches are found very accurate. Good agreement is shown between the third-order models and finite-element simulations for rigid particles in a deformable matrix, and the TPA using the WR microstructural parameters *ζ*_{p} and *η*_{p} is improved.

For the first time, TPAs of the thermal–mechanical properties are computed for isotropic systems of Platonic solids at various volume fractions. A significant particle shape effect is predicted for thermal conductivity, whereas the effective elastic moduli are less sensitive to the microstructural configuration. Based on our statistical framework, a large class of materials with arbitrary inclusion shapes can now be easily studied. Moreover, image-based modelling, using micro-computed tomography for example, can now be successfully employed for real material systems.

## Funding statement

A.G. and K.M. acknowledge support from the Department of Energy, National Nuclear Security Administration, under award no. DE-NA0002377, and from Multiscale Design Systems, LLC under the contract no. FA8651-14-C-0045 (Eglin Air Force Base, STTR Phase II project) by the Department of Defense. T.L.J. was supported in part by the U.S. Department of Energy, National Nuclear Security Administration, Advanced Simulation and Computing Program, as a Cooperative Agreement under the Predictive Science Academic Alliance Program, under contract no. DE-NA0002378. G.A. and T.L.J. were also supported in part by the Defense Threat Reduction Agency, Basic Research award no. HDTRA1-13-0010 to University of Illinois, and T.L.J. in part by HDTRA1-14-1-0031 to University of Florida; Dr Suhithi Peiris, program manager.

## Author contributions

A.G. developed the numerical integration/interpolation and statistical sampling methods, conducted numerical experiments and performed subsequent analysis of results. G.A. and T.J. developed the packing algorithm for Platonic solids. K.M. conceived and coordinated the studies, and contributed to analysis of results. A.G. and K.M. drafted the manuscript. All authors approved manuscript for publication.

## Conflict of interests

The authors have no competing interests.

## Appendix A. Third-order estimates for elastic constants

The parameters *Ξ*_{L} and *Ξ*_{U} presented in equation (2.19) are defined as

The TPA for the bulk and shear modulus of a two-phase isotropic material were derived by Torquato [52,53]. These approximations were simplified assuming stiff particles in a deformable matrix (

## Appendix B. Hashin–Shtrikman bounds

The HS lower bound for the effective thermal conductivity [21] of an isotropic two-phase composite is given as

- Received January 28, 2015.
- Accepted March 11, 2015.

- © 2015 The Author(s) Published by the Royal Society. All rights reserved.