On the computation of the Hashin–Shtrikman bounds for transversely isotropic two-phase linear elastic fibre-reinforced composites

Although various forms for the Hashin–Shtrikman bounds on the effective elastic properties of inhomogeneous materials have been written down over the last few decades, it is often unclear how to construct and compute such bounds when the material is not of simple type (e.g. isotropic spheres inside an isotropic host phase). Here, we show how to construct, in a straightforward manner, the Hashin–Shtrikman bounds for generally transversely isotropic two-phase particulate composites where the inclusion phase is spheroidal, and its distribution is governed by spheroidal statistics. Note that this case covers a multitude of composites used in applications by taking various limits of the spheroid, including both layered media and long unidirectional composites. Of specific interest in this case is the fact that the corresponding Eshelby and Hill tensors can be derived analytically. That the shape of the inclusions and their distribution can be specified independently is of great utility in composite design. We exhibit the implementation of the computations with several examples.


Introduction
Fibre-reinforced composites (FRCs) are commonly employed in numerous applications in science and engineering, one of their main uses being to provide improved tensile strength. It is therefore of great interest to possess knowledge of the overall (effective) properties of such materials; of specific interest in this article will be their effective linear elastic properties and the ability to determine useful bounds on these for a variety of fibre microstructures. FRCs typically occur as one of two types: long unidirectional fibre materials where the fibres are aligned (parallel) to some common axis and they are so long that end effects can be neglected, or short fibre materials where the short fibres may or may not be aligned, depending upon the application. Such FRCs are frequently then used as ply phases in order to build up layered composites, see e.g. Tsai [1]. Furthermore, such layered and FRC materials have the potential to be of great use in modern metamaterial applications as was described by Torrent and Sanchez-Dehesa [2] and Amirkhizi et al. [3] for example. The main objective of this article is to pull together ideas relating to FRCs from a variety of diverse publications, placing particular emphasis on the explicit construction and computation of the Hashin-Shtrikman bounds on effective elastic constants, including information regarding the two-point correlation functions associated with the distribution of inclusion phases.
From a mathematical viewpoint, modelling inhomogeneous materials such as FRCs, which possess a microstructure (i.e. a lengthscale that is much smaller than the characteristic lengthscale of loading, e.g. a propagating wavelength) is difficult, because the partial differential equations that arise possess rapidly varying coefficients (in space). Before the early 1960s, there was a paucity of theoretical work regarding the prediction of the effective behaviour of elastic composite materials, i.e. elastic materials composed of two or more so-called phases that are mixed together in order to improve and/or optimize the overall material behaviour in some sense. In particular, very limited information was available regarding the prediction of the fourth-order effective linear elastic modulus tensor C * . Key results were those of Voigt [4] and Reuss [5] in 1889 and 1929, respectively, whose approximations of uniform strain and stress (respectively) throughout the composite are straightforwardly identified as upper and lower bounds (respectively) on the effective elastic properties. For a composite with phases labelled r = 0, . . . , n and with associated elastic modulus tensors C r , these bounds are written Here the bounds are interpreted component-wise and the Reuss (lower) bound C R is defined by where φ r is the volume fraction of the r th phase, and I is the fourth-order identity tensor with components I i jk = (δ ik δ j + δ i δ jk )/2. We have used A −1 to denote the inverse of the tensor A. Later, we discuss technical aspects for the inversion of transversely isotropic tensors. Finally, the Voigt (upper) bound C V is defined by In addition to these bounds, there also existed Eshelby's dilute estimate for particulate composites derived in his seminal paper of 1957 [6]. Such dilute estimates are only valid for small volume fractions φ r (asymptotically correct to O(φ r )) and furthermore, although the establishment of the Reuss/Voigt bounds was an important step, in many cases, the distance between these bounds can be large, meaning that little can be inferred as regards the effective elastic properties themselves. Numerous important contributions to the subject occurred in the early 1960s, the main motivation being the understanding of the overall behaviour of FRCs. At that time, such materials were used in a variety of applications but principally within the aerospace industry because of their high strength-to-weight ratio compared with more conventional structural materials. In particular, Hashin and Shtrikman [7] established a variational principle for elastostatics which they subsequently applied to multiphase (macroscopically isotropic) composites [8]. The resulting bounds on the effective bulk and shear moduli are the celebrated Hashin-Shtrikman bounds. Hashin [9] extended the principle in order to derive bounds for long FRCs where the macroscopic anisotropy is that of transverse isotropy and phases are isotropic. Hashin [10] also derived results for fibre reinforced composites where phases are transversely isotropic but not via the classical Hashin-Shtrikman variational scheme. Derivations of the Hashin-Shtrikman bounds have been improved and revised by many authors since they were originally devised. In particular, we note the works [11] and [12][13][14][15][16]. The analogous problem for nonlinear heterogeneous media has been treated less frequently. In this regard, we can refer to Willis [17] where the approach of Hill [18] was applied to nonlinear dielectrics, and Talbot and Willis [19] where theoretical aspects of the extension of the Hashin-Shtrikman variational principles to nonlinear heterogeneous systems were formulated. Other variational structures for nonlinear media, based on comparisons with linear homogeneous media that allow the estimation of the effective energy function of nonlinear composites, in terms of the corresponding properties for linear problems with the same microscopic distribution, were obtained by Ponte Castañeda [20,21].
The articles just described were instrumental in enabling general forms of the Hashin-Shtrikman bounds to be written down for arbitrarily anisotropic composites.
Working in Cartesian coordinates x 1 , x 2 , x 3 , let us denote differentiation of a second-order tensor σ = σ (x 1 , x 2 , x 3 ) with respect to x k by σ ,k , and we adopt the Einstein summation convention for summation over repeated indices. We also introduce the notation x = (x 1 , x 2 , x 3 ). The key step in a most succinct derivation of the Hashin-Shtrikman bounds in the structure developed by Willis [12] is to write the standard governing equations of elastostatics for an inhomogeneous material in an alternative form. Upon neglecting body forces and supposing that the elastic modulus tensor of the inhomogeneous material is C = C(x), the equations of elastostatics are div σ T = 0, σ = C(x)e, (1.2) where σ is the Cauchy stress tensor, div is the divergence operator, e is the linear strain tensor, and the superscript T here denotes transpose. In common with much of the existing literature, we use notation AB to denote contraction over the final two and first two indices of the tensors A and B, respectively. This is rather convenient as it can be applied to tensors of arbitrary order. Here the i jth component of Ce is therefore C i jk e k = C i jk e k and also later, we will use this notation extensively for the contraction of two fourth-order tensors A and B so that the i jk th component of AB is A i jmn B nmk . The alternative form of (1.2) used by Willis is derived by writing Hooke's law as σ = C c e + τ (x) where C c is a homogeneous linear elastic comparison material (the choice of which will be discussed shortly), and τ is the so-called polarization stress, which is necessarily spatially dependent. The equations of elastostatics thus reduce to div(σ c ) T = −div τ T , where σ c = C c e. The divergence of the polarization stress can therefore be interpreted as a spatial distribution of body forces in a homogeneous medium with elastic modulus tensor C c . Upon introducing the second-order elastostatic Green's tensor G c associated with the (anisotropic) comparison material, a variational principle on the strain energy is then found in terms of the polarization stress τ (x) [12,16]. Next, we choose the polarization stress to be piecewise uniform (uniform in each phase), usually written as τ (x) = N k=0 χ k (x)τ r where χ k (x) is unity when x resides in the kth phase and is zero otherwise, and each component of τ r is constant.
Optimizing with respect to τ r and choosing the distribution function to have the same spheroidal statistics as the shape of the corresponding inclusion phase, we are led to the classical (and much quoted) form of the Hashin-Shtrikman bounds ((3.13) below) derived by Willis [12]. This form can also be deduced from the more general formulation (3.8) developed by Ponte Castañeda and Willis [16] which allows for more general distributions. The key advantage of the Hashin-Shtrikman bounds over the Reuss-Voigt bounds is that the former uses information about the macroscopic anisotropy; this permits an improvement over the Reuss-Voigt bounds in almost all cases.
It is worth noting at this point that improved bounds, using more general polarizations leading to higher-order statistics, can also be derived [22][23][24]. However, it is often the case that such higher-order statistical correlation information for a given material is not known or difficult to determine accurately.
Although the general form for the Hashin-Shtrikman bounds applicable to arbitrarily anisotropic composites can be derived in a straightforward manner in some cases, only very recently have explicit bounds for generally transversely isotropic materials been written down [25]. Furthermore, as far as the authors are aware, works concerning the construction and computation of such bounds for a given material, using tensor bases (thus permitting fast implementation), are not available. Indeed it appears that the Hashin-Shtrikman bounds have often taken on a mysterious air in the literature. These bounds always appear to be merely stated (not derived constructively) and almost no information can be found in the literature regarding bounds on the effective properties of composites that are not simply either macroscopically isotropic or transversely isotropic of the long FRC variety. Therefore, the step taken in [25] was a useful and important one since explicit bounds were stated and compared with numerous homogenization and micromechanical methods. However, in [25], very little discussion of the required tensor-basis for transverse isotropy was given, and furthermore the uniformity of the so-called P-tensor (directly related to the Eshelby tensor and sometimes known as the Hill tensor in the literature) was not exploited fully since expressions for the P-tensor were given in integral form, rather than explicit expressions which have been previously derived for spheroidal inclusions and distributions.
For all of the reasons above, the construction and computation of the Hashin-Shtrikman bounds for those "not in the know" are far from straightforward. This is important specifically for engineers, materials scientists and industrialists who may wish to construct the Hashin-Shtrikman bounds for a variety of such media. They are often restricted from doing so by lack of information regarding the detail of the construction. It would be convenient to have a systematic and prescriptive way of constructing the Hashin-Shtrikman bounds from the first principles. That is, given the volume fractions, elastic properties, shapes of phases of the composite and their spatial distribution, it would be convenient to have a mechanism by which the Hashin-Shtrikman bounds could be constructed in a straightforward manner using the correct tensor basis set and the appropriate expressions for the Eshelby (S) tensors and/or Hill (P) tensors. In particular in this respect, it appears that although knowledge of explicit bounds is very useful, we would argue that it is less important than being able to construct the bounds from the first principles.
The principal objective of this article is therefore the discussion of such a construction and computation of the Hashin-Shtrikman bounds for transversely isotropic composites restricting attention to two-phase fibre-reinforced media. We consider cases that incorporate separately the influence of the fibre (or "inclusion") phase and the distribution of these inclusions. Initially, we shall consider the general case of an arbitrary number of phases, but considering general statistics makes the study of media with more then one inclusion phase more complicated, and therefore that case shall be considered in a follow-up article.
The paper is organized as follows. In Sect. 2, we introduce notation and the basics concerning Hill's tensor. In Sect. 3, we describe the general formulation of the Hashin-Shtrikman bounds, initially for the general multiphase case before restricting attention to two phases, the primary concern of this article. We specialize in Sect. 4 to the case of macroscopically transversely isotropic materials. This specialization therefore motivates the discussion of transversely isotropic fourth-order tensors in Sect. 5 and in particular we describe a convenient basis set for such tensors which was introduced in this context by Hill [26]. In Sect. 6, we describe how such a basis set, together with knowledge of the appropriate Hill and Eshelby tensors (stated in the appendices), allows us to construct the Hashin-Shtrikman bounds in a straightforward, procedural manner. We discuss the specific limiting cases of layered and long cylindrical fibre-reinforced media in Sect. 7 before illustrating the implementation of the construction via some examples in Sect. 8. The main focus is the illustration of the simple implementation via the basis set employed as well as the study of distributions of inclusions that have a different spheroidal spread to their shape, the latter being associated with the appropriate Hill tensor. We conclude in Sect. 9 indicating required areas for future study.

Notation and preliminaries
Let us first define some language: we shall say that a tensor is uniform if all of its components are constant (and these components can be different constants of course). In this article, we speak generally at first with regard to multi-phase particulate composite materials occupying the domain Ω with closed boundary ∂Ω and with n types of inclusion phase that could be chosen independently of their spatial distribution [16]. We denote the elastic modulus tensor of the r th inclusion phase by C r , r = 1, . . . , n. Inclusion phases are embedded inside a host phase with elastic modulus tensor C 0 . Although results can in principle be obtained for ellipsoidal inclusions, let us restrict attention to the case of spheroidal inclusions; this is convenient for reasons that will become clear shortly. We note that this case is extremely useful since many composites are of this class (including the limiting cases of long FRCs and layered media). Let us denote Ω r as the total domain of the r th phase, i.e. it is the collection of all spheroidal inclusions which constitute that phase, each aligned and having the shape defined by a domain V r . Each spheroid constituting the r th phase is located at a different point in space; in terms of the computation of the Hashin-Shtrikman bounds, only their shape and distribution is important. The total volume fraction of this phase is therefore φ r = |Ω r |/|Ω| where | · | denotes a volume. The volume fraction of the host phase is therefore φ 0 = 1 − n r =1 φ r . In micromechanics and homogenization, a specific tensor known as the P-tensor frequently arises [13]. This fourth-order tensor, also often known as Hill's polarization tensor, has components defined by being linked with the free-space Green's tensor having components G c i j (x − y), i.e. that associated with the linear comparison phase having elastic modulus tensor C c . We note that the integral in (2.1) is over the domain of a spheroid in the r th phase, i.e. V r (hence the superscript r on P) and dependent only on the ratio of semi-axes, and not the size of the spheroid. It is also independent of the material properties of the inclusion so that in the particular case where all inclusion phases are of the same shape, the P-tensor considered here is identical for any phase. The notation | (i j),(k ) indicates symmetrization with respect to these indices, i.e. for a tensor Q i jk such a symmetrization would be If the domain V r is spheroidal as here (ellipsoidal even) then the P-tensor defined in (2.1) is uniform (i.e. each component is constant) [6,27,28]. Note that P c i jk , Hill's tensor associated with the comparison phase, is related to the Eshelby tensor S c i jk via the expression where the second equality is due to the inclusions being ellipsoids so that S r i jmn = S r jimn . The P-tensor associated with an ellipsoidal inclusion of the kth phase, say V k = {x : |A (k) s x| 2 < 1}, for some second-order tensor A (k) s (which is usually diagonal with diagonal components corresponding to the semi-axes of the ellipsoid) may be computed from the expression where S 2 is the unit sphere andξ is a unit vector pointing from the centre of S 2 to a point on its surface. Furthermore, For spheroidal inclusions, explicit forms may be obtained as is detailed in Appendix B. The subscript 's' is associated with shape, to distinguish it from an analogous quantity (with subscript d) that we will encounter shortly associated with distributions.

The Hashin-Shtrikman variational principle
In his paper of 1977, Willis [12] developed a variational structure to derive the Hashin-Shtrikman bounds. Introducing a polarization stress relative to a so-called comparison material, both the inclusion shape and their distribution (from integrals of the associated two-point correlation functions) are incorporated through appropriate influence tensors. The distribution tensors were assumed to be the same as the associated shape tensor for simplicity. Later, Ponte Castañeda and Willis [16] generalized the structure of Willis [12,14], whereby the derived bounds contain general distribution tensors. It makes use of the alternative variational representation for the effective energy function given by Talbot and Willis [19], which also applies for nonlinear composites.
Let us now consider the linear composite described above whose stress-strain constitutive relation is given by Under homogeneous boundary conditions, the fourth-order-effective moduli tensor C * can be defined as and the average strain energy density (see [29]) is determined as Following [7,11], a comparison material is introduced with uniform modulus C c and a polarization field is defined by Certainly, choosing C c −C r ≤ (≥) 0 componentwise achieves this. Finally select a piecewise constant polarization stress Such an assumption subsequently yields at most two-point statistics regarding phase distribution. Under the hypothesis of ellipsoidal symmetry for the distribution of the inclusions (possibly with a different ellipsoidal shape to the shape of the inclusion), the following bounds on the effective energy are derived in [16]: is the average of the so-called optimal polarizations τ * k . By linearity, using (3.1), the following expression for the optimal bounds is derived from (3.2): where ≤ is here interpreted componentwise, and where C B = C + refers to an upper bound and C B = C − to a lower bound. These bounds depend on the choice of comparison modulus tensor and the optimal polarizations which satisfy the relations: The tensor parameters M (k ) in (3.4) depend on C 0 and on the microstructure. They can be shown to be symmetric and to have the form (see [16]): Here P k s is the P-tensor associated with the shape of the inclusion in the kth phase as introduced in (3.6). On the other hand, the notation P (k ) d is the (uniform) P-tensor associated with the ellipsoidal distribution tensor concerning interaction between the kth and th phases. It is associated with the ellipsoid V and may be computed from the expression: Recall that the subscripts s and d refer to the shape and distribution of inclusions, respectively. By definition, we have P d , and for conciseness, we will write P k d , to denote P (kk) d , k = 0, . . . , n. Given the above, let us therefore recall the Hashin-Shtrikman variational principle which states the following regarding the choice of the elastic modulus tensor C c of the comparison material and the way that this leads to bounds, denoted by C B = C + for the upper and C B = C − for the lower bounds. The theorem that we state has been stated in a number of different ways depending upon the context and application. See e.g. [16].

The Hashin-Shtrikman bounds are thus
Therefore, the lower bound C B in (3.7) is given whenever C r − C c is positive semi-definite for all r = 0, . . . , n. If C r − C c is negative semi-definite, the upper bound is given.
In some instances, instead of obtaining bounds, authors have chosen to obtain an approximation to the effective properties by choosing C c = C 0 , the host elastic modulus tensor, regardless of its relationship to other moduli [16]. Let us consider this case first.

The Hashin-Shtrikman estimates
Let us take C c = C 0 . This immediately yields τ * 0 = 0 and due to the form of (3.4), we are able to write down a simple system for the polarizations, i.e. for k = 1, 2, . . . , n, we have and for this general case, one simply has to solve this linear system for the optimal polarizations, which are then fed into the expression for C B = C * to obtain an estimate (note that it is not a bound unless C 0 is either minimal or maximal over all phase modulus tensors). Various cases enable clean simplifications and explicit expressions. Firstly suppose that although we allow the shape and elastic moduli to vary between phases, the distributions are the same for all so that P (k ) It is then straightforward to show that the estimate simplifies to A little work enables the form (3.8) to be written as If all of the inclusions have the same shape, then P r s = P s and furthermore, if it transpires that P s = P d , then the form (3.9) simplifies significantly to

The Hashin-Shtrikman bounds
As the title suggests we shall restrict attention to two-phase particulate media so that there is a single inclusion phase. The bounding technique for greater than one inclusion phase is not well established. Suppose that the inclusion phase is a distribution of (possibly different sized) aligned spheroids but where each spheroid has the same aspect ratio δ = a 3 /a and where the long/short axis (same direction as the semi-axis a 3 of the spheroid) is aligned with x 3 . Note that this spheroidal shape is taken into account thanks to the tensor P s . Their distribution is accounted for by virtue of the tensor P d which we shall consider to be governed by spheroidal statistics, of aspect ratio . We shall consider this notion further shortly. Let us first consider the case when the comparison phase can be chosen as either the host or inclusion phase. Note that this may not always be possible, however.

Comparison phase can be identified as either host or inclusion phase
Let us first suppose that we are able to identify C c = C 0 . In this case, τ 0 = 0 and it is straightforward to determine that and this forms the lower bound if C 0 ≤ C 1 or the upper bound if C 0 ≥ C 1 . Analogously, if we are able to identify C c = C 1 , then τ 1 = 0 and Of course, it may not always be the case that the comparison material can be identified as exactly one of the inclusion phases. As such, it is instructive to determine general bounds in terms of the comparison modulus tensor.

Comparison phase cannot be identified as either host or inclusion phase
In this case, both polarization stresses have non-zero components. As such, the linear system is which we shall conveniently write as so that we can solve to find that Other forms can also be found.

The distribution P-tensor is identical to the shape P-tensor
If P s = P d = P, one can write down a clean form involving the comparison phase. We have This also generalizes to the multi-phase case, so that we can write

Distribution P-tensor
Let us now consider the P-tensor P d which as noted above is associated with a spheroid of aspect ratio whose long or short axis is in the x 3 direction. We must consider the relation between the aspect ratio of the inclusions δ and that of P d , i.e. . Let us refer to Fig. 1) where we note that the distribution spheroid is a safety spheroid, containing a single spheroidal inclusion and which is not intersected by any other security spheroid (hence the terminology). The details of the distribution tensors for more than one inclusion phase become quite complex and will be discussed in a future article. Here we focus on a single inclusion phase. Assuming that P s and P d are given, this construction of the composite means that there is a maximal volume fraction associated with how much of the inclusion can fit into the security spheroid. This depends on whether δ > or δ < (see Fig. 1). Simple calculations show that when δ > , we have Alternatively, consider that the inclusion aspect ratio δ is fixed as well as the volume fraction. This then gives a condition on the maximum permitted. In particular, when < δ, we can determine that 0 ≤ ≤ δ √ φ ≤ δ whereas when > δ, we have δ ≤ ≤ δ/φ. A special case is the example when inclusions are spherical so that δ = 1, when φ max = 1.

Transversely isotropic composites
We have spoken of a general construction of bounds in terms of the elastic modulus tensors C r . Of course, above when it comes to implementation, one must consider a specific form and in general, this can cause great difficulty, particularly for anisotropic phases. Let us now restrict attention to two-phase composites where phases are (at most) transversely isotropic (TI), all with x 3 as the axis of transverse symmetry so that the x 1 x 2 plane is the plane of isotropy. This is a very commonly occurring case in practice. Many fibre-reinforced media are of this type with fibres being TI in order to provide high tensile strength. Elastic modulus tensors that are TI are in fact symmetric TI tensors, a concept that we will discuss shortly. Therefore, in the r th phase, the Cauchy stress σ r is linked to the linear elastic strain e r via σ r = C r e r where C r is (symmetrically) transversely isostropic. We use the following notation (following Hill [26]) in component form: Here k r and m r are known as the in-plane bulk and shear moduli of phase r = 1, 0, respectively, and p r is the anti-plane (or longitudinal) shear modulus. We note here that C r 1111 = k r + m r , C r 1133 = r , C r 3333 = n r , C r 1122 = k r − m r , C r 1313 = p r and C r 1212 = (C r 1111 − C r 1122 )/2 = m r . Finally, we note that the engineering notation for the elastic modulus tensor is often used, which expresses stress in terms of strain via multiplication by a six by six matrix C r i j . The coefficients of this matrix are related to the elastic modulus tensor by the expressions C r 1111 = C r 11 , C r 1133 = C r 13 , C r 3333 = C r 33 , C r 1212 = C r 66 , C r 1313 = C r 55 . Note that in the above, we have worked with the elastic modulus tensor and the Hill notation for this. We also note the following strain-stress relationship in terms of the so-called engineering constants where E T and E A are the transverse and axial Young's moduli and ν T and ν A are the transverse and axial Poisson ratios, respectively. In particular, ν A = ν 31 corresponds to the lateral contraction when a force is applied in the axial direction, and we note the relationship ν A 31 /E A = ν A 13 /E T . Note also that there are still only five independent constants due to the relationship E T r = 2m r (1 + ν T r ). The relationships between the Engineering moduli and the constants appearing in the Hill notation of the elastic modulus tensor are given in appendix.
As detailed above, we restrict attention to inclusions (and distributions) that are spheroidal and with the axis of TI of the phases being the x 3 axis. This yields macroscopically TI materials. Note that the limits δ → ∞ and δ → 0 correspond to the long cylindrical fibre and layer (or penny-shape) limits, respectively.
The macroscopic stress-strain law associated with the effective material will be that in (4.1)-(4.3) with r replaced by * everywhere (denoting the effective material). In what follows then we shall describe how to construct the Hashin-Shtrikman bounds on the effective elastic properties k * , * , n * , m * and p * and subsequently on the Engineering moduli E T * , E A * , ν T * and ν A * .

Transversely isotropic tensors
A fourth-order isotropic tensor can be defined with respect to the tensor basis set {I (1) i jk , I i jk } where i jk = I i jk − I i jk .
Therefore, a fourth-order isotropic tensor H i jk can be written as Given the basis set (5.1), an isotropic tensor H i jk may therefore be specified by a list of two parameters, in shorthand notation: H = (X 1 , X 2 ). In order to specify an isotropic elastic modulus tensor C i jk , frequently the short-hand notation C = (3κ, 2μ) = (X 1 , X 2 ) is adopted where κ = λ + 2μ/3 corresponds to the bulk modulus and μ to the shear modulus. The parameters λ and μ are the Lamé constants of the material. Analogously, suppose now that instead of isotropy we consider the case of transverse isotropy, where the axis of transverse symmetry is the x 3 axis. In this case, a transversely isotropic (TI) tensor H i jk can also be defined with respect to a tensor basis set. Several of these have been proposed; all slight variants of each other, see for example [30,31]. However we shall use the basis associated with the notation introduced for TI tensors by Hill [26] since this is a frequently used notation in the composite materials community, particularly when the Hashin-Shtrikman bounds are discussed. Although it does not appear that Hill wrote down this specific basis set, it is clear that this was what his notation referred to. This basis set enables a TI tensor H i jk to be written in the form: where X n are constants and the basis tensors H (n) i jk are defined by where Θ i j = δ i j − δ i3 δ j3 . The notation H has been used to signify the Hill basis. The set of symmetric TI tensors (H i jk = H k i j ) is defined by setting X 2 = X 3 (e.g. a TI elastic modulus tensor falls into this class). The set of symmetric TI tensors is not closed under multiplication in the sense that if we take two symmetric TI tensors A i jk and B i jk and perform the double contraction A i jmn B nmk = D i jk , the resulting tensor D i jk is not symmetric TI because D 3311 = D 1133 , and so X 2 = X 3 in this case; see Eq. (5.8) below. Hence in general, when considering transversely isotropic tensors, we must consider them with respect to the six basis tensors H (n) i jk defined above; this basis set is closed with respect to double contraction. Note this important difference between transversely isotropic tensors and symmetric transversely isotropic tensors. For more discussion of this point, see [32], although the terminology symmetric TI tensors is specific to this article.
Analogous to the short-hand notation introduced for isotropic tensors above, in order to specify a TI tensor H i jk , the shorthand notation H = (2k H , H , H , n H , 2m H , 2 p H ) = (X 1 , X 2 , X 3 , X 4 , X 5 , X 6 ) is adopted, i.e. the tensor is defined by these six constants (a symmetric TI tensor corresponds to H = H ). This notation is motivated by the fact that Hooke's law with respect to this basis (4.1)-(4.3) links stress to strain via the (symmetric) transversely isotropic tensor C i jk which in short-hand notation is written C = (2k, , , n, 2m, 2 p) with = . Before we proceed, let us define the shorthand notation:  (6) for contraction between the basis tensors defined in (5.3)-(5.5). The contractions defined in (5.6) are summarized in Table 1.
With the notation above, using the contractions in Table 1, we can therefore write down the operations of addition and contraction (in short-hand notation) on the two TI tensors H 1 i jk and H 2 i jk defined in short-hand notation by H 1 = (2k 1 , 1 , 1 , n 1 , 2m 1 , 2 p 1 ) and H 2 = (2k 2 , 2 , 2 , n 2 , 2m 2 , 2 p 2 ), respectively. We define the operation of addition H i jk = H 1 i jk + H 2 i jk in short-hand notation by and the operation of double contraction H i jk = H 1 i jmn H 2 nmk by H = (4k 1 k 2 + 2 1 2 , 2k 1 2 + 1 n 2 , 2k 2 1 + 2 n 1 , n 1 n 2 + 2 1 2 , 4m 1 m 2 , 4 p 1 p 2 ). (5.8) Note that even if 1 = 1 and 2 = 2 (corresponding to symmetric TI tensors H 1 and H 2 ), the second and third elements on the right-hand side of (5.8) are not equal in general, i.e. contraction between two symmetric TI tensors produces a non-symmetric TI tensor in general. This is why we must write out the theory for the evaluation of the Hashin-Shtrikman bounds for TI materials by using the general TI tensor basis set. We note that (5.8) is a correction to the typographical error in the corresponding expression in Appendix A of [25].
The fourth-order identity tensor I i jk and isotropic basis tensors I (1) i jk and I (2) i jk , written with respect to this TI basis are Finally the inverse of the tensor H 1 = (2k 1 , 1 , 1 , n 1 , 2m 1 , 2 p 1 ) is

Construction and computation of the bounds
Here we describe the procedure for constructing the Hashin-Shtrikman bounds for TI composites of the type considered above. This procedure can be easily coded in any standard mathematical package or alternatively in widely used programming languages, in a straightforward manner. In particular, we have described the construction in a way that should be simple to implement in such languages, by defining functions that represent the operations of taking the inverse and double contraction of TI tensors.
The construction is achieved by using the short-hand notation for TI tensors described in the previous section. In particular, we specify a TI tensor by writing this short-hand list of six constants as a 6-vector (2k, , , n, 2m, 2 p), with symmetric transverse isotropy if = . We then exploit the properties of the Hill basis defined above by defining the appropriate tensor operations of contraction and inversion as in (5.8) and (5.10). This vector notation is particularly convenient for computational implementation of the computation of the HS bounds.
We also define the fourth-order identity tensor and isotropic basis tensors in 6-vector forms as in (5.9). Two TI comparison materials are defined (as 6-vectors) as in Remark 1: This means that we take the minimal and maximal elements over all of the phases. For a multi-phase material therefore, it could be that the comparison phase does not correspond to any of the specific phases constituting the composite, since, for example, the modulus m 1 of phase 1 may be maximal but the modulus n 2 of phase 2 could be maximal.
We now define the operations of double contraction and inversion in terms of the 6-vector notation introduced here. With reference to (5.8) and (5.10), we introduce the functions C[H 1 , H 0 ] : R 6 × R 6 → R 6 and I[H 1 ] : R 6 → R 6 operating on the 6-vectors H 1 = (2k 1 , 1 , 1 , n 1 , 2m 1 , 2 p 1 ) and H 0 = (2k 0 , 0 , 0 , n 0 , 2m 0 , 2 p 0 ) in the following way: Such operations can therefore be straightforwardly coded by introducing such functions with arguments as 6-vectors. Therefore, given the (non-symmetric) TI Eshelby tensor (defined in Appendix A) corresponding to the appropriate comparison phase, written in 6-vector form as S δ = (2k δ , δ , δ , n δ , 2m δ , 2 p δ ) the super or subscript δ (or ) refers to the aspect ratio of the spheroid in question involved in the Eshelby tensor. The P-tensor (also written as a 6-vector) can thus be defined via (2.2) as (6.4) and we note that P c is a symmetric transversely isotropic tensor. We also recall that the P tensor is associated with the comparison material. Let us now consider the different computations, according to the choice of comparison material. Suppose first that we are able to choose the comparison material as either the host or inclusion as in Sect. 3.2.1.

Construction when comparison phase is host or inclusion phase
Start by letting N = I[C 1 − C 0 ]. In the case when C c = C 0 , define and then we have Similarly, in the case when C c = C 1 then we define and then we have which is a lower bound if C 0 − C 1 is positive semi-definite or an upper bound if C 0 − C 1 is negative semi-definite. If we have P d = P s = P, then one can simply use this in the procedure above, or alternatively one could use the contraction and inversion procedures to code up the approach written down in Sect. 3.2.3. In this two phase case (a single inclusion phase), re-writing in this manner is no more advantageous. However, in the multi-phase case, this alternative approach does make the construction more simple because we may exploit the form (3.13). We will describe this further in future work.
Since we can write the effective elastic modulus tensor in the short-hand notation form C * = (2k * , * , * , n * , 2m * , 2 p * ), (6.5) we have therefore constructed the bounds on each of the properties, i.e. from (3.7) and (6.5), we have As regards the Engineering constants, we can follow Hill [26] who noted the expression Therefore according as to whether P is positive or negative we can derive bounds on ν A * and E A * from the bounds on k * in (6.6). Furthermore from [10], we can then use their expressions to state the following bounds on the transverse properties: where

Construction in the general case
Let us now consider the most general case when the host or inclusion phases cannot be identified as the comparison material. Let us take into account the formulation in Sect. 3.2.2 and using the contractions and inversion operations described in (6.2) and (6.3), we introduce the following fourth-order tensors: Having defined (6.7), from (3.3) we can therefore write This provides lower (upper) bound whenever C r − C c , r = 0, 1 is positive (negative) semi-definite.

Limiting cases of transversely isotropic composites
We note that the two limits δ = → 0 and δ = → ∞ correspond to a layered medium and a long fibre-reinforced medium, respectively. As regards the former, the effective properties of such media are known exactly. They are the well-known [33] results, which for a material with TI phases, of the type considered here take the form: where the superscript L denotes a layered medium. It is straightforward (but time consuming) to derive the results (7.1)-(7.3) analytically using the construction of the Hashin-Shtrikman bounds in the form (3.13). In this specific limiting case, the Hashin-Shtrikman bounds coincide for all five effective properties of course, thus giving the exact results (7.1)-(7.3) for a layered medium (or a medium consisting of aligned uniformly distributed penny-shaped inclusions).
In the limit δ → ∞, the Hashin-Shtrikman bounds do not coincide; they become the form of bounds given in [9] on long FRCs for k * , m * and p * . Bounds on n * and * were not given in [9] but they may be determined by appealing to the so-called universal properties of composites [18]). In [9], Hashin dealt with only isotropic phases, although since he wrote the paper in terms of (de-coupled) in-plane and anti-plane problems and wrote the expressions in terms of shear moduli and in-plane bulk moduli, it is straightforward to derive the bounds for TI phases by simple replacement of the corresponding moduli with their TI counterparts. As in the layered medium case above, it is straightforward to derive these bounds analytically in this limiting case using the details regarding transversely isotropic tensors in Sect. 5 (once again noting the use of (5.8) and (5.10) in particular) in the construction of the Hashin-Shtrikman bounds described above, together with the exact form of the Eshelby tensor in this limit (given in 1). These results for an n phase composite of the type considered in this article are the following, in terms of k B , B , n B , m B and p B , specified in terms of the comparison phase: and therefore lower (upper) bounds follow from setting the comparison material equal to the smallest (largest) moduli, as described above. Therefore in this limit we are not as fortunate as the layered case where the bounds coincide. Many theories have been proposed to model the effective elastic properties of long FRC materials and a first check on the efficacy such theories is that the prediction sits within the bounds (7.4)-(7.8).

Implementation
We now implement the construction above in various scenarios. Table 2 states the phase properties (taken from [25]). In particular, note that traditionally the Hashin-Shtrikman bounds are plotted as a function of volume fraction φ = φ 1 of the inclusion phase, for a fixed aspect ratio δ of inclusion. Here we not only plot this case, but also consider the bounds plotted as a function of δ for a fixed volume fraction of inclusion. This was also done in [25] but note that here the P-tensor is computed analytically, and therefore it is easy to plot such results for a wide Table 2 Elastic material properties of the phases in GPa  In the last case, the bounds coincide as can also be seen clearly in Fig. 6. These correspond to the Voigt (m * ) and Reuss ( p * ) bounds in those cases range of δ (we plot results for a logarithmic scale in δ). Moreover, here we also study the dependence of the bounds on the distribution aspect ratio ( ). Such plots are very informative particularly since the Reuss-Voigt bounds are independent of δ and and so they cannot give this information.
Example 1: Isotropic glass inclusions in the isotropic epoxy host phase Let us here suppose that P s = P d so that = δ. Consider the example of isotropic glass inclusions uniformly distributed in the epoxy host phase. First we consider the case of spherical glass inclusions, and therefore the effective material is isotropic. In Fig. 2, we plot bounds on the effective bulk modulus κ * and shear modulus μ * of the composite material. In particular, we plot the Reuss (lower) and Voigt (upper) bounds on these effective properties together with the Hashin-Shtrikman bounds. Note the improvement of the Hashin-Shtrikman bounds over the Reuss/Voigt bounds.
Alternatively, we can consider the case of spheroidal glass inclusions with general aspect ratio δ so that the effective material is transversely isotropic. In Fig. 3, we plot the Hashin-Shtrikman bounds on the effective shear moduli m * /m 0 and p * / p 0 . Note that the bounds coincide in the (layered) limit when δ → 0. In the case of the in-plane shear modulus, the bounds become equal to the Voigt bound, whereas in the anti-plane shear modulus they become equal to the Reuss bound; this can be seen from the exact results in (7.3). Note also that the anti-plane shear modulus does not change significantly when we increase δ from δ = 1 to the δ → ∞ limit. This can also be seen on the right of Fig. 6 where we plot the bounds on the effective shear moduli as a function of δ for φ = 0.3. In Fig. 4, we plot the effective properties k * /k 0 and n * /n 0 , noting that the layered limit bounds coincide. We also k * /k 0 n * /n 0 note that in this instance, the bounds on n * /n 0 for δ → ∞ are very tight. The different regimes of * / 0 are plotted in Fig. 5.
In Fig. 6, we illustrate the dependence of the bounds as a function of δ for fixed volume fraction φ = 0.5. In particular, we note the recovery of the (exact) layered expressions (7.1)-(7.3) as δ → 0 (denoted by disks) and also the recovery of the classical Hashin-Shtrikman bounds (7.4)-(7.8) for long FRCs when δ → ∞ (also denoted by disks). Note in particular that the fibres do not have to be particularly long before they reach this limit: δ = O (10) suffices for this limit to be reached, with the exception of n * which requires a larger δ, perhaps δ = O(100). For the k * /k 0 case, we have only plotted the upper bound for the δ = 1 case since δ has to be extremely small before it departs from this significantly. This can be seen in Fig. 10. In the case when δ → 0, the bounds coincide as can also be seen clearly in Fig. 10 Example 2: Transversely isotropic PZT-7A inclusions in the isotropic Epoxy II host phase Let us once again suppose that P s = P d so that = δ, but now consider the effect of a TI inclusion phase (PZT-7A). We shall embed such inclusions in an epoxy phase with properties designated as Epoxy II in Table 2. We plot Hashin-Shtrikman bounds as a function of φ, the inclusion volume fraction in Figs. 7, 8 and 9 for different values of δ, the inclusion aspect ratio. In Fig. 10, we plot the bounds as a function of δ for a fixed volume fraction, φ = 0.3.
In Fig. 7, we plot the bounds on the effective shear moduli of this composite. Note in particular the high contrast achieved by addition of the inclusion phase (particularly the in-plane shear modulus). As in example 1, the effective shear moduli are not affected greatly by increasing δ from unity, particularly the anti-plane shear modulus. This is also seen in the plot on the right of Fig. 10. In Fig. 8, we plot corresponding bounds on k * /k 0 and n * /n 0 . For k * /k 0 , we plot only the upper bound for δ = 1 since this does not change significantly until δ is very small, when the upper and lower bounds coincide in the layered limit as can be seen in the plot on the left. This can also be seen in Fig. 10 (noting the logarithmic scale for δ). As in the isotropic case, bounds on n * /n 0 become relatively tight in the long cylindrical fibre limit, δ → ∞. In Fig. 9, we plot * / 0 , on the left as we decrease δ below unity and on the right as we increase δ above unity. On the left, we see that bounds become progressively tighter as δ is decreased until the layered limit is reached. On the right, we see that the bounds do not become modified greatly by increasing δ above unity. For all properties, as δ → 0, the layered limit is achieved, as should be expected. We note from Fig. 10 that this limit is not reached as quickly as in the isotropic phase case in example 1. In particular even for δ = 10 −3 , this limit has not yet been reached for the in-plane shear modulus m * /m 0 .

Example 3: Spheroidal glass inclusions with spheroidal distribution in the epoxy host phase
Let us now consider a case when the distribution spheroid has a different shape to that of the inclusion phase. Let us first take spherical glass inclusions (δ = 1) and consider both prolate ( > 1) and oblate ( < 1) distributions, referring to Fig. 1. We shall take the case of the isotropic epoxy host phase. Following the discussion in Sect. 3.3 for a given volume fraction, we take the distribution spheroid to be = √ φ < 1 in the oblate case and = 1/φ in the prolate case and plot the resulting bounds in Figs. 11 and 12, noting how tight the Hashin-Shtrikman bounds are in these cases and that the distribution has an apparently greater effect on m * and n * than on p * and * . Of course, each point on these curves represents a different type of composite in the sense that the distribution has different spheroidal distribution statistics. If we were to fix (i.e. fix the type of distribution), then we would have a maximum volume fraction φ max < 1. As an example, in Fig. 13, we therefore take = 1/ √ 2, 1, 2 and plot bounds on the shear moduli up to a maximal (in the non-spherical distribution cases) volume fraction of φ = 0.5. Bounds on the in-plane shear modulus m * are very tight, particularly in the case of the case = 1/ √ 2. The bounds for m * tend to shift upwards with a decrease in and vice versa, whereas with p * , the bounds move down with a decrease in .
Let us now consider the effect when the inclusion becomes a spheroid with aspect ratio δ possibly distinct from the aspect ratio of the spheroidal distribution. In particular, we refer to Fig. 14 13 Bounds on the effective in-plane shear modulus m/m 0 (left) and effective anti-plane shear modulus p/ p 0 (right) for a glass/epoxy composite with spheroidal inclusions with a spheroidal distribution. We plot Hashin-Shtrikman with = 2 (dashed lines) and = 1/ √ 2 (dotted lines) and additionally for spherical distributions, = 1 (dot-dashed lines). Voigt-Reuss bounds (solid lines) are also represented. The maximum volume fraction in the spheroidal distribution cases is φ = 0.5 to the plot in Fig. 6 when the spheroidal distribution is = δ √ φ = δ/ √ 2 < δ. Note that the plot has log-linear scaling and as should be expected, the most significant effect is felt away from the limiting cases when δ → 0 and δ → ∞ corresponding to the layer and long fibre limits.
Finally with reference to Fig. 15, we study the glass/epoxy composite with fixed volume fraction φ = 0.2 of glass inclusions. We then fix the inclusion aspect ratio and study the dependence of the shear moduli m * /m 0 and p * / p 0 on the spheroidal distribution aspect ratio . In particular, we present results for the four values δ = 1.5, 1  (spheres), δ = 0.6 and 0.2. In plotting the range of , we take into account the inequalities given in Sect. 3.3. In particular from subfigure to subfigure, we see the shift in bounds as δ is decreased. Note in particular that there is little variation in the bounds as increases above unity. More variation is seen for < 1.

Concluding remarks
We have shown how to construct, in a straightforward manner, the Hashin-Shtrikman bounds for transversely isotropic composites focusing in particular on the case of two-phase FRCs. The shape of the inclusion phase and their corresponding distributions can be chosen independently by use of the appropriate Hill and Eshelby tensors and TI tensor basis set. Note in particular that for TI materials, the Eshelby tensor can be derived analytically, and we summarized its various forms for spheroids in the Appendix. We implemented different computations for two specific composite materials, exhibiting the improvement of the Hashin-Shtrikman bounds over the Reuss-Voigt bounds, also showing how bounds behave for different inclusion (δ) and distribution ( ) aspect ratios, including the limiting cases of layered media (δ → 0) where bounds coincide, agreeing with the Backus expressions [33] and also long fibre-reinforced media (δ → ∞) where the bounds are then the classical Hashin-Shtrikman bounds (some of which were derived by Hashin [9]) for such media. Entirely analogous approaches may be followed for phases and materials of arbitrary anisotropy, although in general, the corresponding Green's tensor (and hence P-tensor) cannot be derived analytically.
Future work needs to consider the construction and computation of the Hashin-Shtrikman bounds for multi-phase composites, a case that has not been studied sufficiently in the literature. 4mk E A mn + k E A , and , ,

Eshelby tensor for an isotropic comparison phase
For an isotropic comparison phase defined by the Lamé moduli λ c and μ c , S c i jk takes the form (see [6,35], although note the slightly different notation involving the I i j ): , I 13 = δ 2 (I 1 − I 3 ) δ 2 − 1 (10.9) and thus the Eshelby tensor for the isotropic case is defined purely in terms of I 1 = I 1 (1) from (10.4). In the limiting case of an isotropic (v i = 1) sphere (δ = 1), the appropriate limit of (10.4) gives I 1 = 4π/3 which completely defines the Eshelby tensor via (10.5)-(10.9). Only in this spherical, isotropic case does S 1133 = S 3311 .

Eshelby tensor for a transversely isotropic comparison phase
For a TI comparison phase defined by elastic properties k c , c , n c , m c and p c , the components S i jk take the form [34] (although note some typographical errors there that are corrected here): 14) where D = 1 4π p c v 3 , Note that we have used the notation I 3 (v i ) = −2I 1 (v i ) + 4π/v i in (10.10)-(10.15) whereas [34] used I 2 (v i ) instead of I 3 (v i ) in the corresponding equations. This is to preserve the symmetry with the isotropic case above. Finally, let us define two very useful limits of the Eshelby tensor for TI comparison materials.
Layer (penny-shaped) limit (strongly oblate) If we take the limit as δ → 0, we find that Long cylindrical fibre limit (strongly prolate) If we take the limit as δ → ∞, we find that The P-tensor The P-tensor is defined from (2.1), (2.2) and (6.4) and is a symmetric transversely isotropic tensor, which we write in short-hand notation as