The influence of two-point statistics on the Hashin–Shtrikman bounds for three phase composites

In this work we analyze the inﬂuence of the spatial distribution function, introduced by Ponte Castan˜eda and Willis [1], on the Hashin-Shtrikman bounds on the eﬀective transport properties of a transversely isotropic (TI) three-phase particulate composite, i.e. when two distinct materials are embedded in a matrix medium. We provide a straightforward mechanism to construct associated bounds, independently accounting for the shape, size and spatial distribution of the respective phases, and assuming ellipsoidal symmetry. The main novelty in the present scheme resides in the consideration of more than a single inclusion phase type. Indeed, unlike the two-phase case, a two-point correlation function is necessary to characterize the spatial distribution of the inclusion phases in order to avoid overlap between diﬀerent phase types. Moreover, once the interaction between two diﬀerent phases is described, the scheme developed can straightforwardly be extended to multiphase composites. The uniform expression for the associated Hill tensors and the use of a proper tensor basis set, leads to an explicit set of equations for the bounds. This permits its application to a wide variety of phenomena governed by Laplace’s operator. Some numerical implementations are provided to validate the eﬀectiveness of the scheme by comparing the predictions with available experimental data and other theoretical results.


Introduction
The prediction of the effective transport and elastic properties of multiphase materials has attracted the attention of scientists and engineers over many years now.Such materials, whether they be naturally occurring or synthetic, frequently exhibit enhanced physical and mechanical properties.The determination of such macroscopic or effective properties of reinforced materials, polymers, biomaterials and the exploration of the nature of hydrocarbon reservoirs are just a few of the many applications.From a mathematical point of view, the exact prediction of the effective properties of media characterized by a microstructure is generally an impossible assignment, since the associated physical phenomena are governed by partial differential equations with rapidly varying coefficients.
Before the early sixties, key results were given by Wiener [2] in 1912 for the effective conductivity and Voigt [3] and Reuss [4] in 1889 and 1928 respectively, in the case of elastic modulus tensor.These latter results were identified as upper and lower bounds (C R ≤ C * ≤ C V ) on the effective elastic properties C * , and for a composite whose rth constituent modulus is labelled C r , r = 0, ..., n they are , ξ r : rth phase volume fraction.
Wiener's results provided the same form of bounds but for the transport scenario.The rapid development of the aerospace industry initiated numerous contributions to the subject in the early sixties; especially as regards the understanding of the overall behaviour of more complicated geometries such as fibre-reinforced composites.In particular Hashin and Shtrikman [5] established a variational principle for elastostatics which they subsequently applied to multiphase (macroscopically isotropic) composites [6] and the resulting bounds have become known as the Hashin-Shtrikman (HS) bounds on multiphase media.The advantage of the HS bounds over the Reuss-Voigt bounds is that the former use information about the macroscopic anisotropy of the composite, and thus permit an improvement over the Reuss-Voigt bounds in almost all cases.Derivations of the HS bounds have been improved by many authors since they were originally devised, see e.g.[7], [8] and for a recent exposition of the nature of their construction [21].
Of particular note is the work by Ponte-Castañeda & Willis [1], who introduced a comparison material and included additional microstructural information represented by a two-point correlation function.This permitted the derivation of a more general methodology for n distinct types of inclusion phases that could be selected independently of their spatial distribution, although their applications were all associated with two phase materials.Indeed over the last few decades, although a large number of approaches have been proposed to predict the effective properties for the two phase case (see e.g.[2], [4], [9], [10], [3], [11]), the three-phase model (of special significance in the effective thermal conductivity of unsaturated soils for example) has been treated less frequently.
Amongst the authors who have treated such problems is the composite spheres model developed by Friedman [12] for permittivity.This notable example that does permit the study of this case, but then only with very special conditions on the microstructural information.Not all applications will possess this; an interesting application is the prediction of the overall properties of resins reinforced with one dimensional carbon nanotubes (fibres) or two dimensional graphene nanoplatelets (discs) [13], due to their light weight and good chemical resistance compared to more conventional materials.In fact, to study the associated properties and based on experimental analysis or well known expressions, analytical models can be employed ( [14], [15], [16], [17], [18]).However, these schemes have a limited applicability, cannot be used for more general geometries and none consider the spatial distribution effect between different phases.
In this work therefore, our aim in this work is to develop the methodology for three-phase HS bounds in order to accommodate geometries and microstructural parameters that incorporate the spatial distribution of the inclusions.It should be noted that the present scheme is certainly not a simple extension of that studied in [19] for the two phase case.Most importantly, the present case describes the interaction, not only between inclusions of the same phase, but also between different phases, to avoid their overlap, thanks to the probability density function introduced in [1].Therefore, as a consequence, the developed scheme can straightforwardly be employed and extended for the derivation of the HS bounds for multiphase composites.For the particular case of spheroidal inclusions, we analyze the influence of their aspect ratio on the volume fractions of any phase, something that as far as we know has not been clearly described in the literature.
It should be noted that although the general form for the HS bounds applicable to arbitrarily anisotropic composites can be written down in some cases, works concerning the construction of such bounds from first principles (volume fractions, elastic or physical properties, shapes of phases of the composite and their spatial distribution) of a given material are not easily found, if they exist at all, particularly for the three phase case.
For all these reasons, using a tensor-basis for transverse isotropy and exploiting the uniformity of the so-called Eshelby and Hill tensors [22], we construct explicit expressions for the HS bounds for the effective quasi-static transport properties (thermal or electrical conductivity, electrical resistivity or magnetic permeability) for transversely isotropic (TI) three phase media, incorporating information as regards the shape, relative size and distribution of the two filler phases.
Using the symmetrical way in which it relates two different types of inclusion phases, the fact that the Mori-Tanaka model cannot be realized from the Ponte Castanẽda& Willis method will be graphically exposed.The explicit dependence of the obtained formulae on the microstructural parameters provides the possibility of application to a wide range of geometries and not only the case of spherical inclusions and distributions, which is a frequent assumption.
We illustrate the implementation of the scheme with several examples where comparisons with other theoretical predictions confirm that the present model can predict and bound the effective transport properties with accuracy.

Problem statement
It is well known that Laplace's equation governs a significant range of applications, e.g.electrical and thermal conductivity and permittivity, magnetic permeability to name but a few.The mathematical formulation of such problems is therefore identical.In this work to fix notation, we restrict attention to the prediction of the macroscopical electrical conductivity of a bounded heterogeneous medium Ω ⊂ R N that in absence of internal current sources is governed by the Dirichlet problem In (2.1), ∇φ(x) = q(x) is the electric field, j(x) = Σ(x)q(x) the current, −φ(x) is the electrical potential and Σ(x) the conductivity tensor.We will assume that the heterogeneous medium consists of three homogeneous phases: the matrix Ω 0 and Ω k , k = 1, 2, the inclusion phases satisfying Ω = 2 k=0 Ω k , whose distributions are described by some characteristic functions χ k (x), k = 0, 1, 2, taking the value unity when x ∈ Ω k and zero otherwise.The heterogeneous electrical conductivity is defined as Σ(x) = 2 k=0 χ k (x)Σ k , where Σ k , k = 0, 1, 2, are the respective homogeneous conductivity tensors of the different phases, that are occupying constant volume fractions, defined by The problem (2.1) requires boundary conditions on the interface between the three phases, and we assume those of perfect continuity, i.e.
where ν k is the unit vector normal to the boundary ∂Ω k , k = 0, 1, 2.
Following Hill [20] the average behaviour of the composite is defined by its effective energy function where A is the set of admissible current fields A = {j : div j = 0, φ = φ 0 }, and f denotes the spatial volume average of the function f (x) over the domain Ω.Because of linearity, W * (q) = q • Σ * q, where Σ * is the effective conductivity tensor, which satisfies j = ∂W * (q) ∂ q = Σ * q.
In this work we assume: (i) For k = 1, 2, the kth-phase consists of n k ellipsoidal inclusions θ i k,α k whose shapes are defined by matrices E i k,α k ∈ M 3 (R) (the set of symmetrical matrices with real coefficients) in the following sense: Note that the superscript i refers to the shape associated with the inclusion rather than distribution for which we will use a superscript d, as defined shortly.
For every α k = 1, . . ., n k , the ellipsoidal inclusion θ i k,α k has the same shape (but is potentially of different size) as a reference ellipsoid θ i k given by a symmetric matrix E i k (see Fig. 1).consists of n k inclusions, all with the same shape (i.e.aspect ratio) as the associated reference ellipsoid θ i k but with potentially different sizes and orientations.
(ii) Ω is statistically homogeneous, i.e. the probability density p k (x) of finding an inclusion of type k centered at x is a constant, equal to the number of inclusions of type k per unit total volume, i.e.
(iii) Ω has ellipsoidal symmetry [8] for the distribution of the inclusions.This means that for k = 1, 2, the conditional probability density function p (k| ) (x, x ) for finding an inclusion of type k centered at x given that there exists an inclusion of type centered at point x , depends on x = x − x only through the expression which defines the ellipsoid (see Fig. 2) has the same shape (not necessarily size) as a reference ellipsoid θ d (k ) , defined by E d (k ) ∈ M(R 3 ).(iv) Ellipsoidal inclusions are distributed in a manner such that particle overlapping is excluded.This means, assuming that θ i k,α k and θ i ,α are two inclusions from phases k and respectively, centered at x and x respectively, then Remark 2.1.As remarked upon above, the superscripts i and d refer to the ellipsoidal inclusions and distributions respectively, characterizing the arrangement of both inclusion phases.Therefore, from hypothesis (iv) it follows that the ellipsoid θ d (k ),α (k 0 0 ) defining the two point correlation function of the distribution, must be chosen just large enough to guarantee p ( 0 |k0) = 0.Moreover, from the definition of conditional probability, one deduces the symmetry for the ellipsoidal distribution.Finally in the specific case θ d (kk),α (kk) , we write θ d k,α k for simplicity.

Hashin-Shtrikman bounds
Under the assumptions stated above, in the present section we will derive explicit HS bounds on the effective properties of TI three phase materials.Firstly, recall the following proposition (see [19]) for multiphase composites (see also [1], [21] for the elasticity setting).Proposition 3.2.Consider a homogeneous comparison material with uniform conductivity tensor Σ c .Assume that the composite is statistically homogeneous and the distribution of inclusion phases is defined by ellipsoidal symmetry.Then, Hashin-Shtrikman bounds Σ B for the effective conductivity of the heterogeneous media are given by the inequalities is the average of the optimal polarizations τ * k , which satisfy the relations The parameters M (k ) in (3.5) depend on Σ c and on the microstructure.They can be shown to be symmetric and to have the form Here, P k i and P (k ) d denote the uniform Hill-tensors of the inclusion from the kth phase and the spheroidal distribution associated with the interaction between the kth and th phases.See the Appendix for brief details associated with the Hill tensor and a full exposition is given in [22].
where Γ c is the second derivative of the Green's function [22] for an infinite medium with conductivity of the uniform comparison tensor Σ c .
, where Σ + and Σ − are obtained when the comparison material are the largest Σ c = max 0≤r≤n Σ r and the smallest Σ c = min 0≤r≤n Σ r conducting phase respectively.If the comparison material is neither the largest nor the smallest conducting phase, then Σ B is only an approximation to the effective properties and in the following this will be denoted by S.

Hashin-Shtrikman bounds for a three-phase anisotropic composite
Although some aspects of the above have been associated with the multiphase case, from hereon-in we shall restrict attention to the three-phase scenario and for simplicity, we will assume that Σ 2 and Σ 0 are the largest and the smallest conductivity tensors respectively, i.e.: for some arbitrary vector c.
Proof.Set n = 2 in (3.5) and as homogeneous comparison material take Σ c = Σ 0 .Then the polarization tensor τ 0 * = 0. From (3.4) and Remark 3.4 we deduce (3.8)where with (3.6), τ 1 * − and τ 2 * − satisfy which can be written in the form which when all inclusions have the same shape and their distributions have the same symmetry, i.e. (θ d k = θ i k = θ, k = 1, 2) leads to the corresponding expression given by Willis in [8].
Expanding expression (3.12) with respect to the volume fraction, we derive the following approximation for Σ * to second order in the volume fraction, In (3.13), the well-documented fact that the distribution of inclusions first appears at second order in the volume fraction is evident.
Remark 3.7.Expression (3.12) can be compared with the Mori-Tanaka approximation [23] , and S k i = P k i Σ k is Eshelby's tensor for the inclusion phase k = 1, 2 [22].The Mori-Tanaka scheme assumes that an inclusion behaves as an isolated inclusion in an infinite matrix, under some effective uniform boundary conditions in the far field [24], [23].In particular [23] established that the distribution spheroid associated with this scheme needs to be the same as the inclusion.This leads to inconsistencies in terms of microstructural interpretation and in particular here in the threephase (two-inclusion) case, one can see that it will cause an issue because it would give rise to a non-symmetric interaction spheroid θ d and then therefore p (k| ) = p ( |k) , something that is contradictory by definition.This situation therefore cannot be realized from Theorem 3.5.This explains past results for three-phase composites, which illustrated that the Mori-Tanaka method produced results residing outside the Hashin-Shtrikman bounds [25], unless all inclusions are spherical, in which case the conditions for microstructural construction are met.
Figure 3: The Mori-Tanaka method gives rise to problems for the multi-phase scenario when inclusions are spheroidal in the general case.The microstructural interpretation of the scheme is that the interaction spheroids must be the same shape as the associated inclusions.This means in particular that it would yield θ d in the general case, unless all inclusions are e.g.spherical.

A three-phase composite with aligned spheroidal inclusions
Consider now the restriction to the case of spheroidal inclusions with shapes defined by the reference spheroids θ i k , k = 1, 2 associated with the kth phase, where the long or short axis (prolate or oblate) of the spheroid is directed in the x 3 direction, see Figs. 4-6.The inclusion centres are randomly distributed throughout three dimensional space with the mechanism for distribution defined by the distribution spheroids, having aspect ratios ρ d k (for spheroids of phase k interacting with other spheroids from phase k) and ρ for (for spheroids of phase k interacting with spheroids from phase ).The associated reference spheroids for these distributions are defined by θ d k and θ d 12) .This arrangement will clearly lead to a transversely isotropic effective medium, with the x 1 x 2 plane being the plane of isotropy.
Recall that inclusion phase k consists of n k inclusions, each with specific spheroid (all with potentially different size) θ i k,α k and associated distribution spheroids 12),α (12) .If the semi-axis along the long/short axes are defined as a r α k , r = i, d for the α k th inclusion of phase k, the associated matrices are defined as 12) are the matrix representation of the reference spheroids  It should be noted that according to the non-overlapping assumption regarding the inclusions (iv), the volume fractions of the three materials as well as the aspect ratios the two families of spheroidal inclusions and their respective spheroidal distributions, must reside in so-called safe ranges related to how much of any inclusion phase can fit into a security spheroid (see Fig. 6).On one hand, this implies that the parameters ρ i k,α k , ρ d k,α k and ρ d (12),α (12) > 0, must be chosen large enough such that ∀x, At the same time, the feasible values of ξ r depend on the size of the inclusions and on whether In this sense, we must select Note that ξ k depend on the shape through the aspect ratio term multiplying the sum and on the inclusion size thanks to the terms under the summation.
As an example, consider phase 2 in Fig. 6, where it follows that given For phase 1 on the other hand, ρ i 1 < 1, ρ d 1 > 1 and of course ρ = 1 still.We therefore deduce that ξ 1 ≤ ρ i 1 /ρ.Alternatively, for k = 1, 2, fixing ρ i k , ξ k and ρ d k (or ρ), is possible to derive conditions on the maximum ρ (or ρ d k ) permitted.

Construction of the bounds
In this section we construct the Hashin-Shtrikman bounds for the threephase transversely isotropic (TI) composites described above, allowing them to be coded up explicitly in a direct manner in any mathematical package.We assume that the phases can also be TI with principal axes of anisotropy aligned with the axes of the spheroid and plane of isotropy being the x 1 x 2 plane.
The basic idea follows from observing that any second order TI tensor Σ ∈ R 3×3 , can be defined by using two arguments Σ r = (σ r 1 , σ r 2 ), r = 0, 1, 2, with σ r 2 /σ r 1 = σ r the coefficient of anisotropy, thanks to the proper definition of the following TI tensor basis set where With this notation, we can write down the operations of addition, contraction and the inverse (in short-hand notation) on two TI tensors Σ k = (σ k 1 , σ k 2 ) and Σ = (σ 1 , σ 2 ) respectively by Thanks to the explicit expressions (see in Appendix (5.18)-(5.20))for the P-tensor, (3.16) are employed when such expressions are introduced into (3.8)-(3.10)for the determination of the HS bounds, we deduce them in terms of the two components of the Hill tensor for phase k, P kr = (p 1 kr , p 2 kr ) where k = 1, 2 (phase) and r = i, d, m (inclusion, same phase distribution, mixed phase distribution (when no k will be present)).This leads to the following explicit expressions for the HS bounds (3.17)

Numerical results
In this section we illustrate the HS bounding scheme via its implementation for some specific examples, illustrating the mechanism of accounting for various geometries at the microscale.The resulting bounds are then compared with some experimental data and other predictive schemes.

Effective thermal conductivity of an epoxy matrix reinforced with carbon nanotubes (CNTs) and graphene nanoplatelets (GNPs)
Epoxy resin is a frequently used polymer matrix as it is considered a viable substitute for metal alloys due to its light weight and good chemical resistance.However, the low thermal and electrical conductivity constitute limits over its applications in certain areas.To overcome this problem, in recent studies CNTs and GNPs have attracted attention as fillers in epoxy resins since their 1D and 2D structure leads to an outstanding enhancement of the thermal conductivity of the resulting composites.
Furthermore, experimental data from [26] shows that their hybrid consideration leads to an enhancement of the properties with respect the use of one single filler.In this section, we apply our scheme to predict the effective properties of a three phase composite made up a host of epoxy resin with inclusions of carbon and graphene for different kinds of inclusions and spatial distributions of the fillers.Following [26] and [27] the (isotropic) thermal conductivities for the epoxy resin, CNTs (prolate inclusions) and GNPs (oblate inclusions) are respectively σ 0 j = σ e = 0.201 W/(Km), σ 1 j = σ c = 3000 W/(Km) and σ 2 j = σ g = 2000 W/(Km), j = 1, 2. In Fig. 7, we plot the bounds for the effective thermal conductivity on transverse (left) and axial (right) directions, plotting as a function of ξ = ξ 1 = ξ 2 .The significant enhancement of thermal conductivity is clearly visible.In Fig. 8 we plot the HS bounds on the effective transverse (on the left) and axial (on the right) conductivities as a function of ρ i k , for ξ k = 0.3, k = 1, 2. The limiting cases corresponding to carbon nanofibers (CNFs) (ρ i 1 → ∞) and GNPs (ρ i 2 → 0) are also considered, illustrating that fibers do not have to be particularly long to reach the conductivity limit: ρ i 2 = O (10).Note also, that bounds coincide in the layered case.

Effective thermal conductivity of a three-phase soil
The prediction of the effective thermal conductivity and electrical resistivity of a three-phase medium made up water, soil and air, are commonly employed in several applications such as artificial ground freezing, the food industry and for the prediction of the nature of sedimentary formation for petroleum exploration.In this example we present bounds for the transverse effective thermal conductivity of an unfrozen soil, comparing the obtained results with some experimental data and another predictive model based on coated spheres ( [28]).Following [28] the numerical values of thermal conductivities for the (isotropic) soil, water and air are respectively σ 0 j = σ s = 2.5 W/(Km), σ 1 j = σ w = 0.6 W/(Km) and σ 2 j = σ a = 0.026 W/(Km), j = 1, 2. In Fig. 9, we plot the Hashin-Shtrikman bounds for different kinds of inclusions with several geometrical distribution, showing also their improvement over the Voigt-Reuss bounds.In Fig. 10, we study how the effective property depends on the shape of the inclusion and/or distribution (on the right).Also, we compare the HS with some experimental data from [28], showing that these data are consistent with the bounds.

Concluding remarks
The Hashin-Shtrikman bounds are routinely employed to bound the effective properties of composites.In particular they are frequently written down for two-phase particulate composites where the distribution of inclusions is assumed to be the same as their shape.When this is not the case, as has been illustrated here, their construction is non-trivial and significant consideration has to be given to the manner of their construction and the associated distribution spheroids.Here we have illustrated how to construct the bounds in  the context of three-phase particulate composites, firstly considering the general case, before restricting attention to aligned spheroids, inducing transverse isotropy.By employing five parameters, the bounds take into account the shapes of the spheroidal inclusion phases and their corresponding spatial distributions [1].Before now, although thought has been given to this problem, it does not appear that the construction of such bounds for multi-phase materials with spheroidal distributions has been carried out, particularly as regards the mixed phase pair-correlation function.Incorporating this illustrates why, for example, the Mori-Tanaka method will reside outside the Hashin-Shtrikman bounds.
Here we have illustrated the implementation of the HS bounds for a number of examples.Proper use of the bounds ensures that correct HS bounds can be constructed on multiphase composites.

Appendix: The Hill (P) tensor
In this section, we discuss the nonzero components (P 11 , P 22 , P 33 ∈ R), of the second order P-tensor for spheroids with aspect ratio δ.This is employed for both the shape and distribution tensors for spheroids merely with an appropriate modification of the aspect ratio.The explicit expressions below, can be deduced following e.g. the details in [22].In particular consider the particular case of a spheroidal inclusion with semi-axes a 1 = a 2 = a = a 3 = c (see Fig. 5) embedded inside a transversely isotropic (TI) comparison phase with axis of symmetry along x 3 and with elastic modulus tensor Σ c = σ c 1 (Θ mn + σ c 2 σ c 1 δ m3 δ n3 ).In shorthand notation introduced in Section 3.3, we write this as Σ c = (σ c 1 , σ c 2 ) with σ c 2 /σ c 1 = σ c the coefficient of anisotropy and we also define the aspect ratio

Figure 1 :
Figure 1: Figure illustrating the inclusion phase geometries.The kth inclusion phase (k = 1, 2)consists of n k inclusions, all with the same shape (i.e.aspect ratio) as the associated reference ellipsoid θ i k but with potentially different sizes and orientations.

Figure 2 :
Figure2: Within the composite, inclusions comprising the kth phase and associated distributions all have the same shape but can be of different size and orientation, as illustrated here.

Figure 4 :
Figure 4: Aligned spheroidal distribution in a three-phase composite.

Figure 8 :
Figure 8: Plots of the upper and lower Hashin-Shtrikman bounds on the effective properties of an epoxy/CNTs&GNPs as a function of the CNTs and GNPs inclusion of aspect ratios ρ i k , k = 1, 2, for fixed volume fractions ξ k = 0.3, k = 1, 2. With dashed-lines, the case ρ i k = ρ d k = δ.With solid lines, 2δ = ρ d k = ρ i k = δ, k = 1, 2.

Figure 9 :
Figure9: Plots of the upper lower Hashin-Shtrikman bounds on the effective properties of an unfrozen soil with inclusions of water and air.In both figures, an identical spherical distribution (ρ i k = ρ d k = δ) is employed for both kinds of inclusions for the HS bounds (dashed-lines).The Voigt-Reuss bounds (solid lines) are also shown.Different symmetries are represented with dot-dashed lines, both with spherical mixed distribution (δ = 1).On the left spherical inclusions (ρ i 2 = 1) of air with oblate distribution (ρ d 2 = 0.5) and spherical inclusion of water (ρ i 1 = ρ d 1 = 1).On the right spherical (ρ i 1 = 1) inclusions of water with prolate distribution (ρ d 1 = 1.5).

Figure 10 :
Figure 10: On the left, plots of the upper and lower HS bounds (solid lines) on the effective properties of an unfrozen soil with spherical inclusions of water and air, spherically distributed (ρ i k = ρ d k = δ).Voigt-Reuss bounds (dot-dashed lines) and experimental data (circles marked).On the right, upper and lower HS bounds on the effective properties of soil/water&air as a function of the water and air inclusions aspect ratios ρ i k , k = 1, 2 for fixed volume fractions ξ k = 0.3, k = 1, 2. With solid lines, the case ρ i k = ρ d k = δ and with dashed-lines, when ρ i k = ρ d k = δ, k = 1, 2.