The Eshelby, Hill, Moment and Concentration Tensors for Ellipsoidal Inhomogeneities in the Newtonian Potential Problem and Linear Elastostatics

Abstract

One of the most cited papers in Applied Mechanics is the work of Eshelby from 1957 who showed that a homogeneous isotropic ellipsoidal inhomogeneity embedded in an unbounded (in all directions) homogeneous isotropic host would feel uniform strains and stresses when uniform strains or tractions are applied in the far-field. Of specific importance is the uniformity of Eshelby’s tensor $\mathbf{S}$. Following Eshelby’s seminal work, a vast literature has been generated using and developing Eshelby’s result and ideas, leading to some beautiful mathematics and extremely useful results in a wide range of application areas. In 1961 Eshelby conjectured that for anisotropic materials only ellipsoidal inhomogeneities would lead to such uniform interior fields. Although much progress has been made since then, the quest to prove this conjecture is still not complete; numerous important problems remain open. Following a different approach to that considered by Eshelby, a closely related tensor $\mathbf{P}=\mathbf{S}\mathbf{D}^{0}$ arises, where $\mathbf{D}^{0}$ is the host medium compliance tensor. The tensor $\mathbf{P}$ is associated with Hill and is of course also uniform when ellipsoidal inhomogeneities are embedded in a homogeneous host phase. Two of the most fundamental and useful areas of applications of these tensors are in Newtonian potential problems such as heat conduction, electrostatics, etc. and in the vector problems of elastostatics. Knowledge of the Hill and Eshelby tensors permit a number of interesting aspects to be studied associated with inhomogeneity problems and more generally for inhomogeneous media. Micromechanical methods established mainly over the last half-century have enabled bounds on and predictions of the effective properties of composite media. In many cases such predictions can be explicitly written down in terms of the Hill tensor, or equivalently the Eshelby tensor and can be shown to provide excellent predictions in many cases.

Of specific interest is that a number of important limits of the ellipsoidal inhomogeneity can be taken in order to be employed in predictions of the effective properties of, for example, layered media and fibre reinforced composites and also to the cases when voids and cracks are present. In the main, results for the Hill and Eshelby tensors are distributed over a wide range of articles and books, using different notation and terminology and so it is often difficult to extract the necessary information for the tensor that one requires. The case of an anisotropic host phase is also frequently non-trivial due to the requirement of the associated Green’s tensor. Here this classical problem is revisited and a large number of results for problems that are felt to be of great utility in a wide range of disciplines are derived or recalled. A scaling argument leads to the derivation of the Eshelby tensor for potential problems where the host phase is at most orthotropic, without the requirement of using the anisotropic Green’s function. The Concentration tensor $\boldsymbol{\mathcal{A}}$ linking interior fields to those imposed in the far-field is derived for a wide variety of problems. These results can therefore be used in the various micromechanical schemes.

Directly related to the tensors of Eshelby and Hill is the so-called Moment tensor $\mathbf{M}$. As well as arising in the literature on micromechanics, this tensor is important in the vast area of research associated with inverse problems and specifically with the problem of identifying an object inside some domain given the application of a specific set of boundary conditions. Due to its fundamental importance and direct link to the Eshelby and Hill tensors, here we state the connection between $\mathbf{M}, \mathbf{P}$ and $\mathbf{S}$ in an effort to ensure that the work is of use to as wide a community as possible.

Both tensor and matrix formulations are considered and contrasted throughout. Appendices give various details that illustrate the implementation of both approaches.

1 Introduction

The canonical isolated inhomogeneity problem has been of fundamental importance in a number of materials modelling problems now for well over a century. This problem is the following: a single inhomogeneity or particle of general shape, with different material properties to that of the surrounding material is embedded inside an unbounded (in all directions, i.e., free-space) homogeneous host medium. Given some prescribed conditions in the far-field, what form do the fields take within the inhomogeneity? As well as being interesting in its own right, this problem is of utmost importance in homogenization, micromechanics and multiscale modelling.

The first to consider this kind of inhomogeneity problem was Poisson in 1826 [103] who studied the perturbed field due to an isolated ellipsoid in the context of the Newtonian potential problem. He showed that given a uniform electric polarization (or magnetization), the induced electric (or magnetic) field inside the ellipsoid is also uniform. In 1873 Maxwell [83] derived explicit expressions for this field. Early work in linear elasticity saw a number of studies determine the field inside and around inhomogeneities, including the important case of a cavity (since this was correctly recognized as a defect or flaw). Examples of these works are those associated with the case of spheres [43, 118], spheroids [29] and ellipsoids [107, 113, 114] but all considered specific loadings, usually of the homogeneous type in the far-field, meaning uniform tractions or displacements that are linear in the independent Cartesian variable, say $\mathbf{x}$.

The inhomogeneity problem is now usually associated with the name of Eshelby because in 1957 he showed that for general homogeneous conditions imposed in the far-field, the strain set up inside an isotropic homogeneous ellipsoid is uniform [31]. In 1961 Eshelby [32] conjectured that “…amongst closed surfaces, the ellipsoid alone has this convenient property…”. One could ask is this really true? In the sense of what it is thought that Eshelby meant when he made this conjecture (the so-called weak Eshelby conjecture, where the interior field must be uniform for any homogeneous far-field loading), this statement certainly is true although this was only proved in 2008, simultaneously by Kang and Milton [56] and Liu [71] in the case of isotropic media. There is a slightly different version (the so-called strong Eshelby conjecture), where the interior field must be uniform only for a specific, single uniform far-field loading. This strong conjecture has not been proven in the context of three dimensional isotropic linear elasticity, although significant progress has been made in the last decade, see [55] for a review. The results obtained in [4] go beyond the weak Eshelby conjecture but still do not fully prove the strong conjecture. Interestingly the weak conjecture for the associated Newtonian potential problem was proved some time before Eshelby’s 1957 elastostatics paper, by Dive in 1931 [24] and Nikliborc [97] in 1932, see also the discussion in [55, 56, 71]. In deriving these results, Dive and Nikliborc proved the converse of Newton’s theorem that if $V_{1}$ is an ellipsoid of uniform density, the gravitational force in $V_{1}$ is zero [58]. The strong conjecture in the context of the potential problem is true in two dimensions [111, 115] but is not true in dimensions greater than two: a non-ellipsoidal counterexample associated with a specific far-field loading (equivalently a specific eigenstress) was found by Liu [71].

It is important to note that the proofs of Eshelby’s conjectures in elastostatics referred to above correspond to simply connected, isotropic inhomogeneities with Lipschitz boundaries. Eshelby’s work was followed up by numerous researchers who considered the general anisotropic case [6, 8, 32, 61, 70, 121, 122, 129, 133]. In 1974 Cherepanov [22] proved that multiple inhomogeneities of non-ellipsoidal shape can interact in order to render the interior fields uniform; see also Kang and Milton [56] and Liu [71] who coined the term E-inclusions for such interacting inhomogeneities. Liu and co-workers have also considered the periodic Eshelby problem in two dimensions [72, 73]. Kang and Milton [56] used their approach to prove Eshelby’s weak conjecture in the context of the fully anisotropic potential problem. Most notably however, it is stressed again that the weak Eshelby conjecture for elasticity has not yet been proved in the context of anisotropic elasticity.

Interest in deriving the Eshelby tensor for non-ellipsoidal inhomogeneities has always been active in order to show that the conjecture holds for specific classes of inhomogeneities. Particular attention has been paid to polygonal and polyhedral inhomogeneities and the associated properties of Eshelby’s tensor [57, 74, 79, 80, 93, 94, 98, 108]. The supersphere case has been considered recently by [18] building on the work by [99, 100, 119]. A general method was developed by Ru [110] in order to obtain an analytical solution associated with a two dimensional inhomogeneity of arbitrary cross section and explicit forms of the stress inside hypotrochoidal and rectangular inhomogeneities were derived. Some analytical expressions have recently been derived for two-dimensional problems in the Newtonian potential and plane elastostatics problems where inhomogeneities are either polygonal or their shape can be described by finite Laurent expansions [140, 141]. Additional useful properties of the Eshelby tensor have been deduced, including the relationship of the averaged Eshelby tensor for non-ellipsoidal inhomogeneities to their ellipsoidal counterparts [126, 137].

More recently the inhomogeneity problem has been studied in the nonlinear elasticity context where results associated with Eshelby’s conjecture have been proved in two dimensions for so-called harmonic materials [59, 60, 112]. Although nonlinear problems are generally more difficult than linear elastostatics, the nonlinearity frees up a number of issues that are more constrained in linear problems. The study of nonlinear problems with dilatational eigenstrain was recently carried out in [135]. Giordano [41] considered the nonlinearly elastic inhomogeneity problem but where the constitutive behaviour is described via expansions in strain (Landau elasticity).

Here, due to the associated vast theoretical and practical interest, attention is restricted to the case of linear problems for ellipsoidal inhomogeneities and associated limits in the context of the Newtonian potential problem and elastostatics. A general approach to deriving the Hill tensor $\mathbf{P}$ and proving many of its properties is to use the integral equation form of the governing equations [132]. In fact Eshelby approached the problem in quite a different manner, using the concept of eigenstrain [31]. Hill [47] considered the so-called polarization (hence P) of an ellipsoid. The review articles of Walpole [123] and Willis [132], who developed the integral form of the Hill tensor, have been very influential and the text of Mura [92] describes the associated Green’s tensor and form of the Eshelby tensor $\mathbf{S}$ for elastostatics in detail. The consideration of isolated inhomogeneity problems allows the derivation of the so-called concentration tensor $\boldsymbol{\mathcal{A}}$, of use in dilute micromechanical schemes, where interactions between inhomogeneities are not important [134]. This tensor is related to the Hill tensor via the expression

$$\begin{aligned} \boldsymbol{\mathcal{A}} &= \bigl[\mathbf{I}+\mathbf{P} \bigl( \mathbf{C}^{1}-\mathbf{C}^{0} \bigr) \bigr]^{-1}, \end{aligned}$$

(1.1)

where $\mathbf{C}^{0}$ and $\mathbf{C}^{1}$ are the modulus tensors of the host and inhomogeneity respectively and $\mathbf{I}$ is the identity tensor. In the field of micromechanics a number of very ingenious approximations have been made that lead to rather excellent predictions of effective properties in the case where interactions amongst inhomogeneities are important (see [81, 127, 132] for broad overviews). Finally it is noted that variational bounds can be conveniently written down in terms of the Hill or Eshelby tensors [17, 45, 102, 105, 128, 130, 132], with those of the Hashin-Shtrikman type being particularly important.

There is no real preference for the direct integral equation approach leading to the Hill tensor, over the Eshelby eigenstrain approach. It is chiefly down to individual preference although it is important to note that Hill’s tensor possesses the major symmetries whereas Eshelby’s does not in general. Some find the notion of eigenstrain rather artificial, although in many cases it is a very useful concept as a means for solving harder problems such as determining fields inside multiple inhomogeneities [91, 138]. The simple relation

$$\begin{aligned} \mathbf{S}=\mathbf{P}\mathbf{C}^{0} \end{aligned}$$

(1.2)

between the Hill and Eshelby tensors, means that deriving one immediately yields the other.

The Hill and Eshelby tensors are of great utility in a number of micromechanical methods and what is quite astonishing is that they can be evaluated analytically in a large number of very important cases. However, results are distributed over a large number of articles, reviews and textbooks, and furthermore often in articles that span a wide range of scientific fields due to the wide ranging applicability of the theory. References dealing with derivations of specific results are those of [10, 14, 28, 81, 106, 123, 132] and [67]. The field is still very much alive, pushed forward by both unresolved theoretical issues as well as applications involving not only inhomogeneities but also cracks and dislocations [95, 139] and by the desire to fully resolve the open issues described above. Recent work has focused in more detail on inhomogeneities of general shape and how these can feed into models of inhomogeneous media with distributions of non-canonical inhomogeneities [14–16, 33–35, 138, 140, 141]. Such studies are important to understand how local stress fields develop in the medium under loading. This is highly dependent upon the inhomogeneity shape. It must always be stressed that the utility of the Eshelby tensor itself for general shapes in micromechanical methods is limited by the fact that fields interior to the inhomogeneity are not generally uniform in such cases and therefore the tensor does not arise as a natural quantity from the governing integral equations ((2.7) and (2.11) below) as it does in the case of ellipsoidal inhomogeneities.

Another important tensor that arises in inhomogeneity problems and has a direct link to the Eshelby and Hill tensors is the Moment tensor $\mathbf{M}$. In the case of ellipsoidal inhomogeneities with homogeneous boundary conditions in the far-field, this tensor is most easily written down in terms of the Concentration tensor as

$$\begin{aligned} \mathbf{M} &= |V_{1}| \bigl(\mathbf{C}^{1}- \mathbf{C}^{0} \bigr)\boldsymbol{\mathcal{A}}, \end{aligned}$$

(1.3)

where $|V_{1}|$ is the volume of the inhomogeneity. For further discussion on the broader nature of the Moment tensor, see the book by Ammari and Kang [1] and also [4]. It plays an important role in the vast area of research associated with inverse problems and specifically with the problem of identifying an object inside some domain given the application of a specific set of boundary conditions. The Moment tensor arises within the conjecture of Pólya-Szegö [104] from 1951 (pre-dating that of Eshelby), which states that the trace of the moment tensor is a minimum for circular (spherical) inhomogeneities in 2(3) dimensions. As should be obvious there are strong connections between this conjecture and the Eshelby conjecture. We refer to [4, 56] for further details.

Here the objective is to gather together important results associated with the Hill and Eshelby tensors for ellipsoidal inhomogeneities in a consistent notation, derive a number of important limiting cases such as those associated with cracks and cavities, derive compact results associated with the anisotropic potential problem and finally derive and state associated Concentration tensors. This should prove useful to many who frequently require the form of the $\mathbf{P}$- or $\mathbf{S}$-tensors in practice but who struggle to find the appropriate reference. The emphasis here is to derive the Hill, Eshelby and Concentration tensors but as is clear from (1.3) the Moment tensor follows straightforwardly from these.

An important point to note is that using the so-called invariant notation, potential and linear elastostatics problems can be considered simultaneously, only that the latter is a higher order tensor analogue of the former. Here however the applications are made distinct in order to stress the different results and mechanisms for deriving these expressions. In particular the results from potential theory feed into those from linear elastostatics. As a result, index notation shall be used almost entirely throughout.

In much of the literature on micromechanics the terms inclusion and inhomogeneity are used interchangeably. However in some cases they are used to make an important distinction. An inhomogeneity is defined as a particle of general shape having different material properties to those of the surrounding medium in which it is embedded. On the other hand the terminology inclusion is used to represent a general shaped region within some medium that has the same properties as the surrounding medium but where this finite inclusion region has been subject to some eigenstrain (e.g., thermal strain). This differentiation is used, e.g., in Mura [92] and Qu and Cherkaoui [106] and it is also adopted here.

In Sect. 2 the integral equation formulation of the inhomogeneity problem is stated, yielding integral equations for the potential gradient and strain inside an inhomogeneity. In Sect. 3 it is illustrated that such fields are uniform when the inhomogeneity is ellipsoidal and the general expressions for the associated Hill tensors are stated. The notion of Concentration tensors is also discussed. In Sects. 4 and 5 specific results are then stated and derived for the cases of the Newtonian potential problem and elastostatics respectively. A closing discussion is given in Sect. 6 describing how the results are used in micromechanical methods together with a summary of current areas of associated research. Numerous important details and results are stated in Appendices in order for this review to be comprehensive but also to aid the flow of reading, in particular the mechanism for representing tensors of certain symmetries with respect to appropriate tensor bases, tensor contractions, and the matrix formulation of operations between tensors are covered in Appendix C.

As many pertinent references are given as possible; it is important to stress that the focus is specifically on the formulation of the Eshelby, Hill and Concentration tensors rather than articles associated with micromechanical methods, of which there are thousands. For the latter the interested reader is referred to the many textbooks that have been written over the last decade, e.g., [14, 28, 54, 67, 106].

2 Integral Equation Formulation

Index notation shall be used for tensors throughout, working in Cartesian coordinates and using repeated subscripts to imply summation. The term unbounded will be used when referring to free-space, i.e., unbounded in all directions. Although a general invariant formulation can be employed to deal with problems in the potential and linear elastostatics context simultaneously [132], this approach can obfuscate details that are important when it comes to deriving specific Hill and Eshelby tensors for given anisotropies and inhomogeneity shapes, which is the main purpose of this article.

Notation is as defined in Fig. 1 for both the potential problem and linear elastostatics. A single isolated inhomogeneity $V_{1}$, for the time being of general shape and with surface $\partial V_{1}$ is embedded (perfectly) inside an unbounded homogeneous medium $V$ and the medium exterior to $V_{1}$ is denoted as $V\backslash V_{1}=V_{0}$. Both materials are considered generally anisotropic so that their material modulus tensors are

$$\begin{aligned} C_{ij}(\mathbf{x}) &= C_{ij}^{1} \chi^{1}(\mathbf{x}) + C_{ij}^{0} \bigl(1-\chi ^{1}(\mathbf{x}) \bigr) \end{aligned}$$

(2.1)

in the context of the potential problem and

$$\begin{aligned} C_{ijk\ell}(\mathbf{x}) &= C_{ijk\ell}^{1} \chi^{1}(\mathbf{x}) + C_{ijk\ell }^{0} \bigl(1- \chi^{1}(\mathbf{x}) \bigr) \end{aligned}$$

(2.2)

in the context of elastostatics. Here the so-called characteristic function associated with a domain $V_{1}$, has been employed, being defined as

$$\begin{aligned} \chi^{1}(\mathbf{x}) &= \textstyle\begin{cases} 1, & \mathbf{x}\in V_{1},\\ 0, & \mathbf{x}\notin V_{1}. \end{cases}\displaystyle \end{aligned}$$

(2.3)

Finally it is noted that the inhomogeneity ($C_{ij}^{1}$ and $C^{1}_{ijk\ell }$) and host ($C_{ij}^{0}$ and $C^{0}_{ijk\ell}$) modulus tensors are uniform tensors, meaning that each component of the tensor is constant but these constants can be different.

Fig. 1

An inhomogeneity $V_{1}$ of general shape and with boundary $\partial V_{1}$ is embedded perfectly inside the host medium $V_{0}$. The classical inhomogeneity problem is to determine the fields that arise within the inhomogeneity and host medium given some far-field condition

2.1 The Potential Problem

Since it is often useful to consider a specific physical problem, certainly in terms of language and terminology, the potential problem is described in the context of steady state thermal conductivity. The equation governing the steady state temperature distribution $T(\mathbf {x})$ in the medium described above and depicted in Fig. 1 is

$$\begin{aligned} \frac{\partial{}}{\partial{x_{i}}} \biggl(C_{ij}(\mathbf{x})\frac {\partial{T}}{\partial{x_{j}}} \biggr) &= 0, \end{aligned}$$

(2.4)

where the assumption is that no heat sources are present. The free-space Green’s function associated with the host phase satisfies

$$\begin{aligned} \frac{\partial{}}{\partial{x_{i}}} \biggl(C_{ij}^{0}\frac{\partial {G}}{\partial{x_{j}}}( \mathbf{x}-\mathbf{y}) \biggr) + \delta (\mathbf{x}-\mathbf{y}) &= 0, \end{aligned}$$

(2.5)

as well as the far-field condition $\lim_{\mathbf{x}\rightarrow \infty }G(\mathbf{x})=0$. Assuming continuity of temperature and normal flux across $\partial V_{1}$, the resulting temperature distribution may be straightforwardly derived in integral equation form as

$$\begin{aligned} T(\mathbf{y}) = T^{*}(\mathbf{y}) - \bigl(C_{kj}^{1}-C_{kj}^{0} \bigr) \int_{V_{1}}\frac{\partial{T}}{\partial{x_{k}}}(\mathbf{x})\frac {\partial{G}}{\partial{x_{j}}}( \mathbf{x}-\mathbf{y})\,d\mathbf {x}, \end{aligned}$$

(2.6)

which holds for all $\mathbf{y}$. Here $T^{*}(\mathbf{y})$ is the solution to the equivalent problem satisfying (2.4) with no inhomogeneity present (or equivalently with $C_{ij}^{1} = C_{ij}^{0}$). Upon taking derivatives of (2.6) with respect to $y_{i}$ and noting the property $\partial G/\partial x_{i}= -\partial G/\partial y_{i}$ it is found that for all $\mathbf{y}\in V$,

$$\begin{aligned} e_{i}(\mathbf{y}) &= e_{i}^{*}(\mathbf{y}) + \bigl(C_{kj}^{1}-C_{kj}^{0} \bigr) \frac{\partial^{2}{}}{\partial{y_{i}}\partial{y_{j}}} \int_{V_{1}}e_{k}(\mathbf{x})G(\mathbf{x}-\mathbf{y}) \,d\mathbf{x}, \end{aligned}$$

(2.7)

where the $i$th component of the temperature gradient has been defined as $e_{i}=\partial T/\partial x_{i}$.

2.2 Elastostatics

The origins of the $\mathbf{P}$-tensor reside in the context of elastostatics rather than in potential problems even though the theory is of course analogous. The $\mathbf{P}$-tensor originated with Hill [47] who also introduced the compact notation (now commonly referred to as Hill notation) for transversely isotropic fourth order tensors, which is summarized in Appendix C.2.3. Walpole [122], Willis [130–132] and Laws [62] amongst others followed this with influential work associated with inhomogeneities of specific shapes, paying particular attention in many cases to the scenarios of discs, fibres and cracks. A number of $\mathbf{P}$-tensors are also stated in the excellent concise review of micromechanics by Markov [81] although unfortunately, some typographical errors are present there and these are corrected here.

The integral equation formulation of the isolated inhomogeneity problem in elastostatics proceeds analogously to the potential problem with an expected increase in complexity. The equations governing the elastic displacement $\mathbf{u}$ in the medium described above and depicted in Fig. 1 are

$$\begin{aligned} \frac{\partial{}}{\partial{x_{j}}} \biggl(C_{ijk\ell}(\mathbf{x}) \frac{\partial{u_{k}}}{\partial{x_{\ell}}} \biggr) &= 0, \end{aligned}$$

(2.8)

where body forces have been neglected. The associated Green’s tensor of the host phase satisfies

$$\begin{aligned} \frac{\partial{}}{\partial{x_{j}}} \biggl(C_{ijk\ell}^{0} \frac{\partial{G_{kr}}}{\partial{x_{\ell}}} \biggr)+\delta _{ir}\delta( \mathbf{x}-\mathbf{y}) &= 0, \end{aligned}$$

(2.9)

as well as the far-field condition $\lim_{\mathbf{x}\rightarrow \infty }G_{ij}(\mathbf{x})=0$, noting that $G_{ij}=G_{ji}$. The resulting displacement field in the medium may be straightforwardly derived in integral equation form as

$$\begin{aligned} u_{i}(\mathbf{y}) = u_{i}^{*}(\mathbf{y}) - \bigl(C_{mnk\ell}^{1}-C_{mnk\ell}^{0} \bigr)\int _{V_{1}}^{}{e_{mn}(\mathbf{x})\frac{\partial{G_{ki}}}{\partial {x_{\ell}}}(\mathbf{x}- \mathbf{y})}\,d{\mathbf{x}}, \end{aligned}$$

(2.10)

which holds for all $\mathbf{y}$. Here $u_{i}^{*}(\mathbf{y})$ is the solution to the equivalent problem satisfying (2.8) with no inhomogeneity present, or equivalently $C^{1}_{ijk\ell}=C^{0}_{ijk\ell}$. As in the potential problem, take derivatives of both sides of (2.10) to form the strain tensor $e_{ij}=(\partial u_{i}/\partial x_{j} +\partial u_{j}/\partial x_{i})/2$ and use the property $\partial G_{ki}/\partial x_{j}= -\partial G_{ki}/\partial y_{j}$ so that, for all $\mathbf{y}\in V$,

$$\begin{aligned} e_{ij}(\mathbf{y}) &= e_{ij}^{*}(\mathbf{y}) + \bigl(C_{mnk\ell}^{1}-C_{mnk\ell }^{0} \bigr) \biggl[ \frac{\partial^{2}{}}{\partial{y_{\ell}}\partial{y_{j}}}\int _{V_{1}}^{}{e_{mn}( \mathbf{x})G_{ki}(\mathbf{x}-\mathbf{y})}\,d{ \mathbf{x}} \biggr] \bigg|_{(k\ell ),(ij)}. \end{aligned}$$

(2.11)

Here the notation $|_{(k\ell),(ij)}$ indicates symmetry with respect to the pair of indices grouped together inside parentheses, so that upon defining

$$\begin{aligned} Q_{mnki} &= \int _{V_{1}}^{}{e_{mn}( \mathbf{x})G_{ki}(\mathbf{x}-\mathbf{y})}\,d{\mathbf{x}}, \end{aligned}$$

(2.12)

one can write

$$\begin{aligned} \frac{\partial^{2}{Q_{mnki}}}{\partial{y_{\ell}}\partial{y_{j}}} \bigg|_{(k\ell),(ij)} &= \frac {1}{4} \biggl( \frac{\partial^{2}{Q_{mnki}}}{\partial{y_{\ell}}\partial {y_{j}}}+\frac{\partial^{2}{Q_{mn\ell i}}}{\partial{y_{k}}\partial {y_{j}}}+\frac{\partial^{2}{Q_{mnkj}}}{\partial{y_{\ell }}\partial {y_{i}}}+\frac{\partial^{2}{Q_{mn\ell j}}}{\partial{y_{k}}\partial {y_{i}}} \biggr). \end{aligned}$$

(2.13)

3 Uniformity of the Hill and Eshelby Tensors

3.1 The Potential Problem

Impose so-called homogeneous temperature gradient conditions (in the language of micromechanics: such conditions would lead to a homogeneous temperature gradient in a homogeneous medium) so that in the far-field when $|\mathbf{x}|\rightarrow\infty$,

$$\begin{aligned} T &\rightarrow\theta_{i} x_{i}, \end{aligned}$$

(3.1)

where $\boldsymbol{\theta}$ is uniform and therefore $T^{*} = \theta_{i} x_{i}$ and $e_{i}^{*}=\theta_{i}$. Referring to (2.7), one then asks, is there an inhomogeneity of any shape that can give rise to a uniform temperature gradient field inside the inhomogeneity, i.e., for $\mathbf{y}\in V_{1}$? If such an inhomogeneity does exist, then (2.7) is only consistent for $\mathbf{y}\in V_{1}$ if the tensor defined as

$$\begin{aligned} P_{ij}(\mathbf{y}) &= -\frac{\partial^{2}{}}{\partial{y_{i}}\partial {y_{j}}}\int _{V_{1}}^{}{G( \mathbf{x}-\mathbf{y})}\,d{ \mathbf{x}} \end{aligned}$$

(3.2)

is also uniform, i.e., is independent of $\mathbf{y}$. The tensor $\mathbf{P}$ with components $P_{ij}$ defined in (3.2) is known as Hill’s Polarization tensor for the potential problem and it possesses the symmetry $P_{ij}=P_{ji}$. If $\mathbf{P}$ is not uniform, it would mean that the assumption of a uniform temperature gradient field inside the inhomogeneity was incorrect.

It transpires that when the inhomogeneity region is ellipsoidal the $\mathbf{P}$-tensor defined in (3.2) is indeed uniform. This is proved in Appendix A.1, where it is also shown that the general form for the $\mathbf{P}$-tensor can be defined in terms of an integral over the surface of the unit sphere $S^{2}$.

Clearly since the integral is over the surface of the unit sphere $S^{2}$, it is sensible to resolve $\overline{\xi}_{i}$ into spherical polar coordinates for the purposes of evaluating this integral. The form of (3.3) illustrates the important general result that the $\mathbf{P}$-tensor is uniform for an arbitrarily anisotropic ellipsoidal inhomogeneity embedded inside an arbitrarily anisotropic host phase. That Eshelby’s (weak) conjecture is true for anisotropic potential problems [56, 71], means that the ellipsoid is the only shaped inhomogeneity for which the interior temperature gradient is uniform under any such far-field condition of the form (3.1). The fact that the strong conjecture is not true in the Newtonian potential problem means that there exists shapes where the interior temperature gradient field can be uniform for specifically chosen homogeneous far-field conditions [71].

To determine the appropriate $\mathbf{P}$-tensor in any circumstance then one can appeal to (3.3) and carry out the necessary integration. Alternatively, as shall be shown in Sect. 4, in many cases it is relatively straightforward to use symmetry arguments and results from potential theory in the isotropic host case together with scalings in some cases of host anisotropy, in order to derive explicit results, often in a more straightforward manner than by directly evaluating the general expression (3.3). In fact in the potential problem context, symmetry arguments and results from potential theory [58] are often sufficient to derive results for many special cases of ellipsoids in host media that are at most orthotropic. The general result (3.3) is thus suitable for more complex anisotropies than orthotropy or for example if the semi-axes of the ellipsoid are not aligned with the axes of symmetry of host anisotropy.

It should be recalled that the $\mathbf{P}$-tensor is independent of the anisotropy of the inhomogeneity and therefore arbitrary anisotropy can be retained for the inhomogeneity domain. The only aspects of the inhomogeneity that influence the $\mathbf{P}$-tensor are its shape and, for anisotropic host phases, its orientation with respect to the axes of anisotropy of the host phase. It is important to note the following three points:

• In the host region $V_{0}$ the temperature gradient is generally not uniform.

• For non-homogeneous temperature gradient conditions in the far-field, the temperature gradient field inside an ellipsoidal inhomogeneity is generally not uniform. However if the prescribed temperature gradient is a polynomial of order $n$, then so is the field inside an ellipsoidal inhomogeneity, see [6]. This is known as the polynomial conservation property for ellipsoids.

• Generally for non-ellipsoidal inhomogeneities in unbounded domains and general shaped inhomogeneities in bounded host domains $V$, the temperature gradient inside the inhomogeneities is not uniform, although interacting E-inclusions [71] can lead to uniform interior temperature gradients and for specific loadings, non-ellipsoidal inhomogeneities can yield uniform interior temperature gradients, e.g., the counterexample of the Strong Eshelby conjecture given by Liu [71].

Regarding the first point, once the interior temperature gradient field is known, it can be employed to determine the exterior field by using (2.6) so that for $\mathbf{y}\notin V_{1}$,

$$\begin{aligned} T(\mathbf{y}) &= \theta_{i} y_{i} - \bigl(C_{kj}^{1}-C_{kj}^{0} \bigr)\mathcal{A}_{k\ell }\theta_{\ell} \int _{V_{1}}^{}{ \frac{\partial{G}}{\partial{x_{j}}}(\mathbf{x}-\mathbf{y})}\,d{ \mathbf{x}}, \end{aligned}$$

(3.6)

where $\mathcal{A}_{ij}$ is the temperature gradient Concentration tensor linking the interior temperature gradient to that in the far-field, i.e., $\theta_{j}$, see Sect. 3.3. The gradient of (3.6) is not uniform since $\mathbf{y}$ now lies outside $V_{1}$.

3.2 Elastostatics

The symmetry relation $C_{ijk\ell}=C_{ij\ell k}$ was used in deriving (2.11) as this turns out to be preferable in various contexts. Analogously to the potential problem, let us take homogeneous displacement gradient conditions in the far-field, i.e., as $|\mathbf{x}|\rightarrow\infty$

$$\begin{aligned} u_{i} &\rightarrow \epsilon _{ij}x_{j}, \end{aligned}$$

(3.7)

where $\epsilon _{ij}$ is uniform and therefore $u^{*}_{i}=\epsilon _{ij}x_{j}$ and $e^{*}_{ij}=(\epsilon _{ij}+\epsilon _{ji})/2$. Note that $\epsilon _{ij}$ does not have to be symmetric but if it is then it is simply the strain in the far-field. As in the case of the potential problem, the aim is to determine if there exists an inhomogeneity of any shape that is consistent with the assumption of uniform interior strain. If such an inhomogeneity exists, (2.11) is only consistent for $\mathbf{y}\in V_{1}$ if the tensor defined as

$$\begin{aligned} P_{ijk\ell}(\mathbf{y}) &= - \biggl[\frac{\partial^{2}{}}{\partial {y_{j}}\partial{y_{\ell}}} \int _{V_{1}}^{}{G_{ik}(\mathbf{x}- \mathbf{y})}\,d{\mathbf{x}} \biggr] \bigg|_{(ij),(k\ell)} \end{aligned}$$

(3.8)

is uniform. The tensor defined here is the $\mathbf{P}$-tensor for elastostatics. It possesses the minor symmetries $P_{ijk\ell}=P_{ij\ell k}=P_{jik\ell }$ by construction. Furthermore, thanks to the symmetry of the free space Green’s tensor $G_{ij}=G_{ji}$ it also possesses the major symmetry $P_{ijk\ell}=P_{k\ell ij}$.

It transpires that when the inhomogeneity region is ellipsoidal the $\mathbf{P}$-tensor defined in (3.8) is indeed uniform. This is proved in Appendix A.1, where it is also shown that the general form for the $\mathbf{P}$-tensor can be defined in terms of an integral over the surface of the unit sphere $S^{2}$.

That Eshelby’s (weak) conjecture is true for isotropic elastostatics problems [56, 71], means that the ellipsoid is the only shaped inhomogeneity for which the interior temperature gradient is uniform under any such far-field condition of the form (3.7). It is stressed however that it is not yet clear whether the weak conjecture is true in the context of anisotropic problems.

To determine the $\mathbf{P}$-tensor for an ellipsoid for a given host anisotropy one merely has to evaluate the surface integral in (3.9) and this can be done numerically very efficiently. For host anisotropies more complex than transversely isotropic it is generally recommended that the form (3.9) be employed and integrals are evaluated numerically. In what follows here the $\mathbf{P}$-tensor shall be determined in the case of an isotropic host phase by appealing to various symmetries and potential theory. An important result derived by Withers [133] associated with a transversely isotropic host phase is also stated.

As in the potential problem, the only aspects of the inhomogeneity that influence the $\mathbf{P}$-tensor are its shape and, for anisotropic host phases, its orientation with respect to the axes of anisotropy of the host phase. Note also that the same three points described for the potential problem, preceding equation (3.6), also hold here in the elastostatics context with appropriate modifications to terminology. Furthermore, once the field is known inside the inhomogeneity region $V_{1}$ the exterior field can be determined in terms of the Green’s tensor, as

$$\begin{aligned} u_{i}(\mathbf{y}) &= \bar{\epsilon }_{ij}y_{j} - \bigl(C_{mnk\ell}^{1}-C_{mnk\ell }^{0} \bigr) \mathcal{A}_{mnpq}e^{*}_{pq}\int _{V_{1}}^{}{ \frac{\partial{G_{ki}}}{\partial{x_{\ell}}}(\mathbf{x}-\mathbf {y})}\,d{\mathbf{x}}, \end{aligned}$$

(3.12)

where $\mathcal{A}_{ijk\ell}$ are the components of the strain Concentration tensor (see Sect. 3.3), which links the strain inside the inhomogeneity to that in the far-field.

3.3 The Potential Gradient Tensor and Strain Concentration Tensor

Defining the volume average

$$\begin{aligned} \overline{f} &= \frac{1}{|V|} \int_{V} f(\mathbf{x}) \,d\mathbf{x} \end{aligned}$$

(3.13)

of the function $f$, it is straightforward to show that in the case of the conditions (3.1), the body averaged temperature gradient is [81]

$$\begin{aligned} \overline{e}_{i} &= \theta_{i} = e^{*}_{i}. \end{aligned}$$

(3.14)

Of immediate interest is the temperature gradient field $e^{1}_{i}$ inside an ellipsoidal inhomogeneity $V_{1}$, which from the theory developed above has been shown to be uniform so that it is equal to its phase average, $e^{1}_{i}=\overline{e}^{1}_{i}$ where the phase average is defined as

$$\begin{aligned} \overline{f}^{1} &= \frac{1}{|V_{1}|} \int_{V_{1}} f(\mathbf{x}) \,d\mathbf{x}. \end{aligned}$$

(3.15)

Therefore in the case of an isolated ellipsoidal inhomogeneity with homogeneous far-field conditions of the form (3.1), using (3.14), the expression in (2.7) becomes

$$\begin{aligned} e^{1}_{i} &= \overline{e}_{i} - \bigl(C^{1}_{kj}-C^{0}_{kj} \bigr)e^{1}_{k} P_{ij}. \end{aligned}$$

(3.16)

Using the symmetries $C_{ij}=C_{ji}$ and $P_{ij}=P_{ji}$ and re-arranging, (3.16) can thus be written in the form

$$\begin{aligned} \overline{e}_{i} &= \bigl(\delta_{ij}+P_{ik} \bigl(C_{kj}^{1}-C_{kj}^{0} \bigr) \bigr)e^{1}_{j}. \end{aligned}$$

(3.17)

Therefore one can relate the uniform temperature gradient inside the inhomogeneity to the average temperature gradient of the entire body via a second order tensor, which is thus identified as the temperature gradient Concentration tensor for this problem.

Note that $\mathcal{A}_{ij}$ is the Concentration tensor associated with an isolated inhomogeneity inside an unbounded host medium. The calligraphic notation $\mathcal{A}_{ij}$ has been used to stress the link with (and distinguish from) the exact Concentration tensor whose components are usually defined as $A_{ij}$, and which links the phase average of the true temperature gradient inside an inhomogeneity to that in the far-field in a complex inhomogeneous medium, which may consist of interacting inhomogeneities. For a dilute medium where interaction effects can be neglected, $A_{ij}=\mathcal{A}_{ij}$.

Moving on to the elastostatics case, it is straightforward to show that in the case of the conditions (3.7), the body averaged strain is

$$\begin{aligned} \overline{e}_{ij} &= \frac{1}{2}(\epsilon _{ij}+ \epsilon _{ji}) = e^{*}_{ij}. \end{aligned}$$

(3.20)

The (uniform) strain $e^{1}_{ij}$ inside an ellipsoidal inhomogeneity $V_{1}$ is thus equal to its phase average, $e^{1}_{ij}=\overline {e}^{1}_{ij}$. Therefore for an isolated ellipsoidal inhomogeneity with homogeneous far-field conditions (3.7), using (3.20), the expression in (2.11) can be used to determine the expression

$$\begin{aligned} \bar{e}_{ij} = \bigl(I_{ijk\ell}+P_{ijmn} \bigl(C_{mnk\ell}^{1}-C_{mnk\ell }^{0} \bigr) \bigr)e_{k\ell}^{1}, \end{aligned}$$

(3.21)

where $I_{ijk\ell}$ is the fourth order identity tensor defined in (C.9). Therefore the uniform strain inside the inhomogeneity can be related to the average strain of the entire body via a fourth order tensor, which is thus identified as the strain Concentration tensor for this problem.

4 The Potential Problem: Specific Cases

4.1 Isotropic Host Phase

Assume that the host phase is isotropic, so that $C_{ij}^{0}=k_{0}\delta_{ij}$ and therefore the associated free-space Green’s function is

$$\begin{aligned} G(\mathbf{x}-\mathbf{y}) &= \frac{1}{4\pi k_{0}}\frac{1}{|\mathbf {x}-\mathbf{y}|}. \end{aligned}$$

(4.1)

From (3.2) and (1.2) therefore

$$\begin{aligned} P_{ij}(\mathbf{x}) = \frac{1}{k_{0}}\frac{\partial^{2}{\varGamma }}{\partial{x_{i}}\partial{x_{j}}}, \qquad S_{ij}(\mathbf{x}) = \frac{\partial^{2}{\varGamma}}{\partial{x_{i}}\partial{x_{j}}}, \end{aligned}$$

(4.2)

where $\varGamma$ is the potential defined by

$$\begin{aligned} \varGamma(\mathbf{x}) &= -\frac{1}{4\pi} \int _{V_{1}}^{}{ \frac{1}{|\mathbf {x}-\mathbf{y}|}}\,d{\mathbf{y}}. \end{aligned}$$

(4.3)

Note that this is the negative of the Newtonian potential (see for example Kellogg [58]) associated with an ellipsoidal domain $V_{1}$. From potential theory therefore

$$\begin{aligned} \nabla^{2}\varGamma(\mathbf{x}) = \frac{\partial^{2}{\varGamma }}{\partial{x_{k}}\partial{x_{k}}} &= \chi^{1}(\mathbf{x}) \end{aligned}$$

(4.4)

and furthermore $\varGamma(\mathbf{x})$ is a quadratic function of the components of $\mathbf{x}$ (see Appendix B), illustrating the uniformity of the $\mathbf{P}$-tensor in this case.

As an aside, note that since the host phase is isotropic, the temperature field exterior to the inhomogeneity is determined via (3.6) so that

$$\begin{aligned} T(\mathbf{y}) &= \theta_{i} y_{i} + \bigl(C_{kj}^{1}-k_{0} \delta_{kj} \bigr)\mathcal{A}_{k\ell }\theta_{\ell} \frac{1}{k_{0}}\frac{\partial{\varGamma(\mathbf{y})}}{\partial {y_{j}}}. \end{aligned}$$

(4.5)

This solution tends to $\theta_{i} y_{i}$ in the far-field, as it should do.

Once $P_{ij}$ is determined for an isotropic host phase the associated Concentration tensor for an isolated inhomogeneity may then be found from (3.19) as is now illustrated in a number of special cases of specific inhomogeneities with given shape and anisotropy.

4.1.1 Sphere in an Isotropic Host Phase

When $V_{1}$ is a sphere, it is clear from (4.3) that $\varGamma (\mathbf{x})$ must be spherically symmetric and hence $\frac{\partial ^{2}\varGamma }{\partial x_{i}\partial x_{j}}$ must be isotropic (and uniform) and therefore

$$\begin{aligned} \frac{\partial^{2}{\varGamma}}{\partial{x_{i}}\partial{x_{j}}} &= \gamma \delta_{ij} \end{aligned}$$

(4.6)

for some constant $\gamma$. Note that the form (4.6) is a result of the spherical shape and not any assumption regarding isotropy of the inhomogeneity because such an assumption has not been made. Performing a contraction in (4.6) and using (4.4) with $\mathbf{x}\in V_{1}$ yields $\gamma =\frac{1}{3}$. Therefore from (4.2) and (1.2)

$$\begin{aligned} P_{ij} = \frac{1}{3k_{0}}\delta_{ij}, \qquad S_{ij} = \frac{1}{3}\delta_{ij}. \end{aligned}$$

(4.7)

For practical purposes, especially for use in micromechanical methods for bounds and estimates of effective material properties, it is useful to write down the associated Concentration tensors.

Isotropic Sphere

If the spherical inhomogeneity is isotropic with conductivity $C_{ij}^{1}=k_{1}\delta_{ij}$, one can show straightforwardly using (3.19) and properties of second order tensors (see Appendix C.1) that

$$\begin{aligned} \mathcal{A}_{ij} &= \frac{3k_{0}}{k_{1}+2k_{0}}\delta_{ij}. \end{aligned}$$

(4.8)

Anisotropic Sphere

Consider a transversely isotropic sphere where the plane of isotropy is the $x_{1} x_{2}$ plane. The conductivity tensor therefore takes the form $C_{ij}^{1} = k_{1}(\varTheta_{ij}+\upsilon\delta_{i3}\delta_{j3})$ where $\varTheta _{ij}$ is defined according to

$$\begin{aligned} \delta_{ij}=\varTheta_{ij}+\delta_{i3} \delta_{j3} \end{aligned}$$

(4.9)

and $\upsilon$ indicates the degree of anisotropy, with $\upsilon=1$ giving isotropy. Using the result derived in (4.7) and (4.9) together with properties from Appendix C.1, the Concentration tensor $\mathcal{A}_{ij}$ can be written down in the form

$$\begin{aligned} \mathcal{A}_{ij} &= \frac{3k_{0}}{k_{1}+2k_{0}}\varTheta_{ij} + \frac {3k_{0}}{\upsilon k_{1}+2k_{0}}\delta_{i3}\delta_{j3}. \end{aligned}$$

(4.10)

Setting $\upsilon=1$ recovers the isotropic result (4.8).

Averaging over all orientations of the anisotropy of the inhomogeneity uniformly will yield an isotropic Concentration tensor of the form

$$\begin{aligned} \underline{\mathcal{A}}_{ij} &= \gamma \delta_{ij}, \end{aligned}$$

(4.11)

where the underline denotes averaging over orientations. By performing this uniform orientation averaging (see Appendix C.3.3) on (4.10) it is straightforwardly shown that

$$\begin{aligned} \gamma &= \frac{2k_{0}}{k_{1}+2k_{0}} + \frac{k_{0}}{\upsilon k_{1}+2k_{0}}. \end{aligned}$$

(4.12)

It is also possible to weight the averaging should there be a preferential distribution of anisotropy. Details of how to implement this procedure are given in Appendix C.3.2.

4.1.2 Circular Cylinder in an Isotropic Host Phase

When $V_{1}$ is a circular cylinder with axis of symmetry in the $x_{3}$ direction, it is clear that $\varGamma(\mathbf{x})$ should be independent of $x_{3}$ and isotropic in the $x_{1}x_{2}$ plane so that

$$\begin{aligned} \frac{\partial^{2}{\varGamma}}{\partial{x_{i}}\partial{x_{j}}} &= \gamma \varTheta_{ij} \end{aligned}$$

(4.13)

for some constant $\gamma$, where $\varTheta_{ij}$ was defined in (4.9). Performing a contraction in (4.13) and using (4.4) the result $\gamma =\frac{1}{2}$ is obtained. From (4.2) therefore

$$\begin{aligned} P_{ij} = \frac{1}{2k_{0}}\varTheta_{ij}, \qquad S_{ij} = \frac{1}{2}\varTheta_{ij}. \end{aligned}$$

(4.14)

If the cylinder is isotropic with conductivity tensor

$$\begin{aligned} C^{1}_{ij} = k_{1}\delta_{ij}, \end{aligned}$$

(4.15)

it is straightforward to show that

$$\begin{aligned} \mathcal{A}_{ij} = \frac{2k_{0}}{k_{1}+k_{0}}\varTheta_{ij} + \delta_{i3}\delta_{j3}. \end{aligned}$$

(4.16)

Alternatively, suppose that the cylinder is transversely isotropic with conductivity tensor

$$\begin{aligned} C_{ij}^{1} &= k_{1}(\varTheta_{ij}+ \upsilon\delta_{i3}\delta_{j3}). \end{aligned}$$

(4.17)

Interestingly one can show that in this case the Concentration tensor is identical to the isotropic case, i.e., the form stated in (4.16): the parameter $\upsilon$ does not appear in the Concentration tensor. Of course if the axis of symmetry of transverse isotropy is not aligned with the cylinder axis then this Concentration tensor would then depend on $\upsilon$.

If the (uniform) orientation average of (4.16) is taken, the associated Concentration tensor is derived:

$$\begin{aligned} \underline{\mathcal{A}}_{ij} &= \frac{k_{1}+5k_{0}}{3(k_{1}+k_{0})} \delta_{ij}. \end{aligned}$$

(4.18)

This last result is often used when very long, thin needle-like inhomogeneities are uniformly distributed and oriented throughout some host medium.

4.1.3 Ellipsoid in an Isotropic Host Phase

Consider now the general case of an ellipsoidal inhomogeneity and as usual denote the semi-axes of the ellipsoid as $a_{j}$, $j=1,2,3$. It is straightforward to show, using the theory of the potential, as in Appendix B that for an ellipsoid in an isotropic host phase, the function $\varGamma(\mathbf{x})$ is quadratic in the components of $\mathbf{x}$ and can be written in the closed form

$$\begin{aligned} \varGamma(\mathbf{x}) &= \biggl(\frac{x_{1}^{2}}{a_{1}^{2}}+\frac {x_{2}^{2}}{a_{2}^{2}}+ \frac{x_{3}^{2}}{a_{3}^{2}}-1 \biggr)\varUpsilon- \sum_{j=1}^{3} \frac {x_{j}^{2}}{a_{j}}\frac{\partial{\varUpsilon}}{\partial{a_{j}}}, \end{aligned}$$

(4.19)

where

$$\begin{aligned} \varUpsilon&= \frac{1}{4}a_{1}a_{2}a_{3} \int_{0}^{\infty}\frac{dt}{\sqrt {(a_{1}^{2}+t)(a_{2}^{2}+t)(a_{3}^{2}+t)}}. \end{aligned}$$

(4.20)

In Appendix B it is then shown that

$$\begin{aligned} \frac{\partial^{2}{\varGamma}}{\partial{x_{i}}\partial{x_{j}}} &= \sum_{n=1}^{3} \mathcal{E}(\varepsilon _{n}; \varepsilon _{1},\varepsilon _{2}) \delta_{in}\delta_{jn}, \end{aligned}$$

(4.21)

where with $\varepsilon _{n}=a_{3}/a_{n}$,

$$\begin{aligned} \mathcal{E}(x;\varepsilon _{1},\varepsilon _{2}) &= \frac{x^{2}}{2} \int_{0}^{\infty}\frac {ds}{(1+sx^{2})\sqrt{(1+s\varepsilon _{1}^{2})(1+s\varepsilon _{2}^{2})(1+s)}}. \end{aligned}$$

(4.22)

Therefore

$$\begin{aligned} P_{ij} = \frac{1}{k_{0}}\sum_{n=1}^{3} \mathcal{E}(\varepsilon _{n};\varepsilon _{1},\varepsilon _{2}) \delta_{in}\delta_{jn}, \qquad S_{ij} = k_{0}P_{ij}. \end{aligned}$$

(4.23)

Finally note that using (4.4) it is easily shown that $\gamma _{1}+\gamma _{2}+\gamma _{3} = 1$.

4.1.4 Spheroid in an Isotropic Host Phase

Denote the semi-axes of the spheroid as $a_{1}=a_{2}=a\neq a_{3}$ and use this in (B.27) which becomes

$$\begin{aligned} \varUpsilon&= \frac{1}{2}a_{3}^{2} \int_{0}^{\pi/2}\frac{\cos\psi}{ \varepsilon ^{2}+(1-\varepsilon ^{2})\sin^{2}\psi}\,d\psi, \end{aligned}$$

(4.24)

where $\varepsilon =a_{3}/a$. Make the substitution $\beta =\sin\psi$ to find

$$\begin{aligned} \varUpsilon= \frac{1}{2}a_{3}^{2} \int_{0}^{1} \frac{d\beta }{\varepsilon ^{2}+(1-\varepsilon ^{2})\beta ^{2}} &= \frac{a_{3}^{2}}{2} \textstyle\begin{cases} \frac{\operatorname{arccosh}(\varepsilon )}{\varepsilon \sqrt{ \varepsilon ^{2}-1}}, & \varepsilon >1, \\ \frac{\arccos(\varepsilon )}{\varepsilon \sqrt{1-\varepsilon ^{2}}}, & \varepsilon < 1, \\ 1, & \varepsilon =1, \end{cases}\displaystyle \end{aligned}$$

(4.25)

noting that $\varepsilon =1$ is the case of a sphere.

Therefore from (4.19) it is clear that

$$\begin{aligned} \varGamma(\mathbf{x}) &= \frac{1}{2} \bigl(x_{1}^{2}+x_{2}^{2} \bigr)\mathcal{T}(\varepsilon ) +\frac{1}{2}x_{3}^{2} \mathcal{S}(\varepsilon )-\varUpsilon, \end{aligned}$$

(4.26)

where

$$\begin{aligned} \mathcal{S}(\varepsilon ) &= \frac{2}{a_{3}^{2}}\varUpsilon- \frac {2}{a_{3}} \frac{\partial{\varUpsilon}}{\partial{a_{3}}} \\ &= \frac{1}{1-\varepsilon ^{2}}-\frac{\varepsilon }{1-\varepsilon ^{2}} \textstyle\begin{cases} \frac{1}{\sqrt{ \varepsilon ^{2}-1}}\operatorname{arccosh}(\varepsilon ), & \varepsilon >1, \\ \frac{1}{\sqrt{1-\varepsilon ^{2}}}\arccos(\varepsilon ), & \varepsilon < 1 \end{cases}\displaystyle \end{aligned}$$

(4.27)

and $\mathcal{T}(\varepsilon )=\frac{1}{2}(1-\mathcal{S}(\varepsilon ))$. The function $\mathcal{S}(\varepsilon )$ has taken many forms in the literature but it is felt that this is a most clear, consistent and concise formulation. Note that $\mathcal{S}(\varepsilon )\rightarrow\frac{1}{3}$ as $\varepsilon \rightarrow1$ for the spherical case (see further details below).

Using (4.26) in (4.2), the resulting Hill and Eshelby tensors take the form

$$\begin{aligned} P_{ij} = \frac{1}{k_{0}} (\gamma \varTheta_{ij} + \gamma _{3} \delta_{i3}\delta_{j3} ), \qquad S_{ij} &= k_{0} P_{ij}, \end{aligned}$$

(4.28)

where $\gamma _{3}=\mathcal{S}(\varepsilon )$ and $\gamma =\mathcal{T}(\varepsilon )$.

Expressions for the Concentration tensors associated with the spheroidal inhomogeneity case can now be determined straightforwardly. For an isotropic spheroid,

$$\begin{aligned} \mathcal{A}_{ij} &= \dfrac{k_{0}}{k_{0}+(k_{1}-k_{0})\gamma }\varTheta_{ij} + \dfrac {k_{0}}{k_{0}+(k_{1}-k_{0})\gamma _{3}}\delta_{i3}\delta_{j3}. \end{aligned}$$

(4.29)

Averaging uniformly over orientations of the axes of the spheroid yields

$$\begin{aligned} \underline{\mathcal{A}}_{ij} &= \frac{1}{3} \biggl( \dfrac {2k_{0}}{k_{0}+(k_{1}-k_{0})\gamma }+\dfrac{k_{0}}{k_{0}+(k_{1}-k_{0})\gamma _{3}} \biggr)\delta_{ij}. \end{aligned}$$

(4.30)

One can take limits in the case of the spheroidal inhomogeneity in order to derive the following results, some of which confirm cases considered above. When $V_{1}$ is

• (i) a sphere, i.e., $\varepsilon \rightarrow1$, it is deduced that $\gamma =\gamma _{3} = \frac{1}{3}$,

• (ii) a cylinder, i.e., $\varepsilon \rightarrow\infty$, it is deduced that $\gamma _{3}= 0$, $\gamma =\frac{1}{2}$,

• (iii) a disc or layer, i.e., $\varepsilon \rightarrow0$, it is deduced that $\gamma _{3}=1$, $\gamma =0$.

When used in (4.29) (i) and (ii) confirm the results derived for the Concentration tensors for isotropic spheres and cylinders derived in Sects. 4.1.1 and 4.1.2 respectively. One has to be rather careful in taking these limits and for (i) use L’Hopital’s rule appropriately, noting that as $\varepsilon \rightarrow1$

$$\begin{aligned} \mathcal{S}(\varepsilon ) = \frac{1}{3} -\frac{4}{15}(\varepsilon -1) + \frac {6}{35}(\varepsilon -1)^{2}+O \bigl((\varepsilon -1)^{3} \bigr). \end{aligned}$$

(4.31)

In (ii) one has to use the fact that $\operatorname{arccosh}x\sim\log x$ as $x\rightarrow\infty$. The result for layers in (iii) can also be obtained via straightforward symmetry arguments.

4.1.5 Limiting Case of an Elliptical Cylinder

One can use the formulation for general ellipsoids above in order to obtain a result for an elliptical cylinder, unbounded in the $x_{3}$ direction with semi-axes $a_{1}$ and $a_{2}$ lying along the $x_{1}$ and $x_{2}$ axes respectively. Taking the limit $a_{3}\rightarrow\infty$ in (B.32), one can show that

$$\begin{aligned} \frac{\partial^{2}{\varGamma}}{\partial{x_{i}}\partial{x_{j}}} &= \sum_{n=1}^{2} \gamma _{n} \delta_{in} \delta_{jn}, \end{aligned}$$

(4.32)

where

$$\begin{aligned} \gamma _{1} = \frac{a_{1}a_{2}}{2} \int_{0}^{\infty}\frac{ds}{(a_{1}^{2}+s)^{\frac {3}{2}}(a_{2}^{2}+s)^{\frac{1}{2}}}, \qquad \gamma _{2} &= \frac{a_{1}a_{2}}{2} \int_{0}^{\infty}\frac{ds}{(a_{1}^{2}+s)^{\frac{1}{2}}(a_{2}^{2}+s)^{\frac{3}{2}}}. \end{aligned}$$

(4.33)

The integrals can be determined explicitly, noting that with $i,j=1,2$, $(i\neq j)$

$$\begin{aligned} \int\frac{ds}{(a_{i}^{2}+s)^{\frac{3}{2}}(a_{j}^{2}+s)^{\frac{1}{2}}} &= \frac {2}{(a_{i}^{2}-a_{j}^{2})} \biggl(\frac{a_{j}^{2}+s}{a_{i}^{2}+s} \biggr)^{\frac{1}{2}} \end{aligned}$$

(4.34)

and therefore

$$\begin{aligned} \gamma _{1} = \frac{a_{2}}{a_{1}+a_{2}}, \qquad \gamma _{2} = \frac{a_{1}}{a_{1}+a_{2}}. \end{aligned}$$

(4.35)

The Hill and Eshelby tensors therefore take the form

$$\begin{aligned} P_{ij} = \frac{1}{k_{0}(a_{1}+a_{2})} (a_{2}\delta_{i1} \delta_{j1} + a_{1}\delta_{i2}\delta_{j2}), \qquad S_{ij} = k_{0}P_{ij}. \end{aligned}$$

(4.36)

As regards the Concentration tensor for an isotropic cylinder with $C^{1}_{ij}=k_{1}\delta_{ij}$, this is determined in the form

$$\begin{aligned} \mathcal{A}_{ij} &= \frac{k_{0}(1+\epsilon )}{k_{0}+k_{1} \epsilon }\delta_{i1} \delta_{j1} + \frac{k_{0}(1+\epsilon )}{k_{1}+k_{0} \epsilon }\delta_{i2}\delta _{j2} + \delta_{i3}\delta_{j3}, \end{aligned}$$

(4.37)

where $\epsilon =a_{2}/a_{1}$ is the aspect ratio of the ellipse.

4.1.6 Limiting Cases of a Cavity, Penny-Shaped Crack and Ribbon-Crack

It does not make sense to define a temperature gradient Concentration tensor in the context of cracks or cavities because clearly there is no interior field. However it turns out that this concept is useful and can be interpreted as linking the far-field to the field on the surface of such inhomogeneities [50] with an appropriate definition of cavity temperature gradient. Here the results above are used in order to derive associated Concentration tensors for cracks and cavities.

Consider a spheroidal inhomogeneity and the limit $k_{1}\rightarrow0$ in (4.29). This yields

$$\begin{aligned} \mathcal{A}_{ij} &= \dfrac{1}{1-\gamma }\varTheta_{ij} + \dfrac{1}{1-\gamma _{3}}\delta_{i3}\delta_{j3}. \end{aligned}$$

(4.38)

This is the Concentration tensor for potential problems involving spheroidal cavities.

Next consider the so-called penny-shaped crack limit. The asymptotic form of $\gamma (\varepsilon )$ and $\gamma _{3}(\varepsilon )$ as $\varepsilon \rightarrow0$ is required. These are easily shown to be

$$\begin{aligned} \gamma _{3}(\varepsilon ) = 1-\frac{\pi}{2}\varepsilon + 2\varepsilon ^{2} + O \bigl(\varepsilon ^{3} \bigr), \qquad \gamma (\varepsilon ) = \frac{\pi}{4} \varepsilon -\varepsilon ^{2} + O \bigl(\varepsilon ^{3} \bigr). \end{aligned}$$

(4.39)

Therefore one can derive the form

$$\begin{aligned} \mathcal{A}_{ij} &= \varTheta_{ij} + \biggl( \frac{2}{\pi \varepsilon } + \frac {1}{2} \biggr)\delta_{i3} \delta_{j3} +O(\varepsilon ), \end{aligned}$$

(4.40)

where expansions have been taken for $\varepsilon \ll1$ and terms up to $O(1)$ have been retained since higher order terms will clearly vanish as $\varepsilon \rightarrow0$.

The coefficient of $\delta_{i3}\delta_{j3}$ in (4.40) involves an apparently singular limit as $\varepsilon \rightarrow0$. That this is not a problem arises from the fact that this expression is used in formulae for effective properties of cracked media where this term is always multiplied by a volume-fraction term (or rather a crack-density term) that is proportional to $\varepsilon$ [48, 49]. Note that taking the limits in the opposite order, i.e., $\varepsilon \rightarrow0$ and then $k_{1}\rightarrow0$ yields an inconsistent result, giving rise to singular effective material behaviour in the crack limit, which cannot be correct.

Finally, consider a different limit, the so-called ribbon-crack limit. Take $k_{1}=0$ in the elliptical cylinder result (4.37) to find

$$\begin{aligned} \mathcal{A}_{ij} &= (1+\epsilon )\delta_{i1}\delta_{j1} + \frac{(1+\epsilon )}{\epsilon }\delta_{i2}\delta_{j2} + \delta_{i3} \delta_{j3}. \end{aligned}$$

(4.41)

Therefore as $\epsilon \rightarrow0$

$$\begin{aligned} \mathcal{A}_{ij} &= \frac{1}{\epsilon }\delta_{i2} \delta_{j2} + \delta_{ij} +O(\epsilon ). \end{aligned}$$

(4.42)

As in the penny-shaped crack result above, the Concentration tensor for the ribbon-crack is singular.

4.2 Anisotropic Host Phase

The general form (3.3) for the $\mathbf{P}$-tensor associated with arbitrary host anisotropy requires the necessary surface integral to be evaluated. In the case of transversely isotropic and orthotropic media however, where principal axes are aligned with the semi-axes of the ellipsoid, the problem can be simplified significantly by employing a scaling of the Cartesian variables in order to reduce the isolated ellipsoidal inhomogeneity problem in an anisotropic medium to the case of an ellipsoidal inhomogeneity (with different semi-axes) in an isotropic medium. Therefore the results derived above for the isotropic host phase case can be used in the scaled domain and then mapped back to the physical domain in order to obtain the appropriate physical Hill and Eshelby tensors.

As usual consider the case of an ellipsoid with semi-axes $a_{j}, j=1,2,3$ but now embedded in an orthotropic host medium (with principal axes aligned along $x_{j}$, i.e., with the semi-axes of the ellipsoid) so that

$$\begin{aligned} C^{0}_{ij} &= k_{0} (\delta_{i1} \delta_{j1}+\upsilon_{2}\delta_{i2} \delta_{j2}+ \upsilon_{3}\delta_{i3} \delta_{j3} ), \end{aligned}$$

(4.43)

where $\upsilon_{2}=1$ (or $\upsilon_{3}=1$) for transverse isotropy. The governing partial differential equation is

$$\begin{aligned} \frac{\partial{}}{\partial{x_{i}}} \biggl(C^{0}_{ij}\frac{\partial {T}}{\partial{x_{j}}} \biggr) &= 0. \end{aligned}$$

(4.44)

Now employ the simple rescaling

$$\begin{aligned} x_{j} = \sqrt{\upsilon_{j}}\hat{x}_{j}, \qquad j =1,2,3, \end{aligned}$$

(4.45)

where $\upsilon_{1}=1$ is introduced for notational convenience (the conductivity along the $x_{1}$ axis is thus $k_{0}$), so that the semi-axes of the ellipsoid in the mapped domain become $\hat{a}_{j}=a_{j}/\sqrt {\upsilon_{j}}, j=1,2,3$ and the scaled ellipsoid is denoted as $\hat {V}_{1}$. The governing equation then becomes that governing isotropic media so that

$$\begin{aligned} P_{ij} &= \frac{1}{k_{0}}\sum_{n=1}^{3} \frac{\partial ^{2}\varGamma(\mathbf{x})}{\partial x_{n}^{2}}\delta_{in}\delta_{jn} = \frac{1}{k_{0}} \sum_{n=1}^{3} \frac{1}{\upsilon _{n}} \frac{\partial ^{2}\hat{\varGamma}(\hat{\mathbf{x}})}{\partial \hat {x}_{n}^{2}}\delta_{in}\delta_{jn}, \end{aligned}$$

(4.46)

where $\hat{\varGamma}(\hat{\mathbf{x}})$ is defined in terms of the isotropic (due to scaling) Green’s tensor as defined in (4.1) but now integrated over the scaled ellipsoid $\hat {V}_{1}$, i.e.,

$$\begin{aligned} \hat{\varGamma}(\hat{\mathbf{x}}) &= -\frac{1}{4\pi} \int_{\hat{V}_{1}}\frac {1}{|\hat{\mathbf{y}}-\hat{\mathbf{x}}|}d\hat{\mathbf{y}}. \end{aligned}$$

(4.47)

As a consequence for a general ellipsoid, the result (B.29) can be used but with $x_{j}$ replaced by $\hat{x}_{j}$ and $a_{j}$ replaced by $\hat{a}_{j}, j=1,2,3$. Therefore with reference to (4.22)

$$\begin{aligned} \frac{\partial^{2}{\hat{\varGamma}}}{\partial{\hat{x}_{i}}\partial {\hat{x}_{j}}} &= \sum_{n=1}^{3} \mathcal{E}(\hat{ \varepsilon }_{n};\hat{\varepsilon }_{1},\hat{ \varepsilon }_{2}) \delta_{in}\delta_{jn}, \end{aligned}$$

(4.48)

where $\hat{\varepsilon }_{n}=\hat{a}_{3}/\hat{a}_{n}$. The $\mathbf{P}$ and $\mathbf{S}$-tensors for an ellipsoid embedded inside an orthotropic host medium can then be written as

$$\begin{aligned} P_{ij} = \frac{1}{k_{0}}\sum_{n=1}^{3} \gamma _{n}\delta_{in}\delta_{jn}, \qquad S_{ij} = \sum_{n=1}^{3} \upsilon_{n} \gamma _{n}\delta_{in}\delta_{jn}, \end{aligned}$$

(4.49)

where $\gamma _{n} = \frac{1}{\upsilon_{n}}\mathcal{E}(\hat{\varepsilon }_{n};\hat {\varepsilon }_{1},\hat{\varepsilon }_{2})$.

Note that the above scaling approach is considerably simpler than carrying out the necessary integrals in the corresponding general expression (3.3) for the $\mathbf{P}$-tensor. If the inhomogeneity is embedded in a transversely isotropic host phase with conductivity tensor

$$\begin{aligned} C_{ij}^{0} &= k_{0}(\varTheta_{ij}+ \upsilon\delta_{i3}\delta_{j3}), \end{aligned}$$

(4.50)

then the (orthotropic) $\mathbf{P}$-tensor for an ellipsoid with semi-axes aligned with the principal directions of anisotropy is found by setting $\upsilon_{1}=\upsilon_{2}=1$ and $\upsilon_{3}=\upsilon$ in (4.49) above. Simplifications arise for a spheroid of course as shall now be illustrated.

4.2.1 Spheroid in a Transversely Isotropic Host Phase

Consider a spheroid in a transversely isotropic medium where the major/minor axis of the spheroid is aligned with the axis of transverse isotropy of the host phase. Denote the semi-axes of the spheroid as $a=a_{1}=a_{2}\neq a_{3}$ and the axis of transverse isotropy as $x_{3}$. The scaling argument above can be used to see immediately that the $\mathbf{P}$ and $\mathbf{S}$-tensors are given by

$$\begin{aligned} P_{ij} = \frac{1}{k_{0}} (\gamma \varTheta_{ij}+ \gamma _{3}\delta_{i3}\delta_{j3} ), \qquad S_{ij} = k_{0}P_{ij}, \end{aligned}$$

(4.51)

where with reference to (4.27)

$$\begin{aligned} \gamma _{3}(\varepsilon ) &= \frac{1}{\upsilon}\mathcal{S} \biggl( \frac{\varepsilon }{\sqrt{\upsilon}} \biggr), \end{aligned}$$

(4.52)

$\varepsilon =a_{3}/a$ and $\gamma =\frac{1}{2}(1-\upsilon \gamma _{3})$, the latter being derived by using

$$\begin{aligned} \frac{\partial^{2}{\hat{\varGamma}}}{\partial{\hat{x}_{1}}^{2}}+\frac {\partial^{2}{\hat{\varGamma }}}{\partial{\hat{x}_{2}}^{2}}+\frac {\partial^{2}{\hat{\varGamma}}}{\partial{\hat{x}_{3}}^{2}} &= 1 \end{aligned}$$

(4.53)

for $\hat{\mathbf{x}}\in V_{1}$.

Assuming the spheroid itself is isotropic with conductivity tensor $C^{1}_{ij} = k_{1}\delta_{ij}$, and using (4.50) together with the form of $\mathbf{P}$-tensor in (4.51), the Concentration tensor defined in (3.19) can be straightforwardly determined as

$$\begin{aligned} \mathcal{A}_{ij} &= \frac{k_{0}}{k_{0}+(k_{1}-k_{0})\gamma }\varTheta_{ij} + \frac {k_{0}}{k_{0}+(k_{1}-\upsilon k_{0})\gamma _{3}}\delta_{i3}\delta_{j3}, \end{aligned}$$

(4.54)

with $\gamma$ and $\gamma _{3}$ as defined above. Alternatively, supposing that the spheroid is now transversely isotropic with the same axis of symmetry as the host, i.e., $C^{1}_{ij} = k_{1}(\varTheta_{ij}+\zeta \delta_{i3}\delta _{j3})$, one finds that

$$\begin{aligned} \mathcal{A}_{ij} &= \frac{k_{0}}{k_{0}+(k_{1}-k_{0})\gamma }\varTheta_{ij} + \frac {k_{0}}{k_{0}+(\zeta k_{1}-\upsilon k_{0})\gamma _{3}}\delta_{i3}\delta_{j3}. \end{aligned}$$

(4.55)

4.2.2 Circular Cylinder in a Transversely Isotropic Host Phase

The circular cylinder limit can be taken in the spheroid case considered in Sect. 4.2.1 where the cross-section of the cylinder sits in the plane of isotropy of the host medium. It is then anticipated that the $\mathbf{P}$ and $\mathbf{S}$-tensors will be transversely isotropic. It has been discussed above that $\mathcal{S}(x)\rightarrow0$ as $x\rightarrow0$ and therefore as with the isotropic host case from (4.52) $\gamma _{3}\rightarrow0$. Therefore $\gamma =\frac {1}{2}(1-\nu \gamma _{3})= \frac{1}{2}$ and then

$$\begin{aligned} P_{ij} = \frac{1}{2k_{0}}\varTheta_{ij}, \qquad S_{ij} = \frac{1}{2}\varTheta_{ij}, \end{aligned}$$

(4.56)

so that in fact this tensor is unchanged from the case of an isotropic host phase as in (4.14). The Concentration tensor for an isotropic cylinder can be straightforwardly determined as

$$\begin{aligned} \mathcal{A}_{ij} &= \frac{2k_{0}}{k_{1}+k_{0}}\varTheta_{ij}+ \delta_{i3}\delta_{j3}. \end{aligned}$$

(4.57)

The Concentration tensor associated with a transversely isotropic cylinder is also given by that in (4.57).

An interesting non-standard example is the case of a spheroid embedded inside a transversely isotropic host phase where the axes of symmetry and semi-axes are not coincident. In this case the general (surface integral) form of the $\mathbf{P}$ and $\mathbf{S}$-tensors must be used with the semi-axes aligned with the $\mathbf{x}$ axes but with all components of the modulus tensor being generally non-zero.

4.2.3 Ellipsoid in an Orthotropic Host Phase

Consider an ellipsoid with semi-axes $a_{j}, j=1,2,3$ that are aligned with the axes of anisotropy of the host medium with orthotropic conductivity tensor as defined in (4.43). Analogous scaling arguments can be used as above in order to scale this problem into an ellipsoid in an isotropic host and then scale back to the physical domain, as described above to show that

$$\begin{aligned} P_{ij} = \frac{1}{k_{0}}\sum_{n=1}^{3} \gamma _{n}\delta_{in}\delta_{jn}, \qquad S_{ij} = \sum_{n=1}^{3} \gamma _{n} \delta_{in}\delta_{jn}, \end{aligned}$$

(4.58)

where

$$\begin{aligned} \gamma _{n} &= \frac{1}{\upsilon_{n}}\mathcal{E} (\hat{\varepsilon }_{n}; \hat{\varepsilon }_{1},\hat{\varepsilon }_{2} ), \end{aligned}$$

(4.59)

where $\upsilon_{1}=1$ and noting that $\hat{\varepsilon }_{n} = \sqrt{\frac {\upsilon_{n}}{\upsilon_{3}}}\varepsilon _{n}$ where $\varepsilon _{n}=\frac{a_{3}}{a_{n}}$.

For reference, $\mathbf{P}$-tensors for a variety of problems are summarized in Table 1.

Table 1 Summary of results for the $\mathbf{P}$-tensor associated with ellipsoidal inhomogeneities for potential problems

5 Elastostatics: Specific Cases

5.1 Isotropic Host Phase

For the case of an isotropic host phase case the elastic modulus tensor is defined as

$$\begin{aligned} C_{ijk\ell}^{0} &= 3\kappa _{0} I_{ijk\ell}^{1}+2 \mu_{0} I_{ijk\ell}^{2} \end{aligned}$$

(5.1)

in terms of the isotropic fourth order basis tensors as defined in (C.7) and (C.8). Here $\kappa _{0}$ and $\mu_{0}$ are the bulk and shear moduli of the host, noting the relation to Poisson’s ratio $\nu_{0}$

$$\begin{aligned} \kappa _{0} &= \frac{2\mu_{0}(1+\nu_{0})}{3(1-2\nu_{0})}. \end{aligned}$$

(5.2)

The appropriate isotropic Green’s tensor is

$$\begin{aligned} G_{ij}(\mathbf{x}-\mathbf{y}) &= \frac{1}{4\pi\mu_{0}} \frac{\delta _{ij}}{|\mathbf{x}-\mathbf{y}|}- \frac{1}{16\pi\mu_{0}(1-\nu _{0})}\frac{\partial^{2}{|\mathbf{x}-\mathbf{y}|}}{\partial {x_{i}}\partial{x_{j}}} \end{aligned}$$

(5.3)

and the expression for the $\mathbf{P}$-tensor in (3.8) therefore becomes

$$\begin{aligned} P_{ijk\ell} &= \dfrac{1}{4\mu_{0}} \biggl(\frac{\partial^{2}{\varGamma }}{\partial{x_{j}}\partial{x_{\ell }}} \delta_{ik}+\frac{\partial ^{2}{\varGamma}}{\partial{x_{j}}\partial{x_{k}}}\delta_{i\ell} + \frac {\partial^{2}{\varGamma}}{\partial{x_{i}}\partial{x_{\ell}}}\delta _{jk}+\frac{\partial^{2}{\varGamma}}{\partial{x_{i}}\partial {x_{k}}}\delta_{j\ell} \biggr) \\ &\quad{}+\dfrac{1}{4\mu_{0}(1-\nu_{0})}\frac{\partial ^{4} \varPsi}{ \partial x_{i}\partial x_{j} \partial x_{k} \partial x_{\ell}}. \end{aligned}$$

(5.4)

The potential $\varGamma$ is that already encountered and defined in (4.3) and $\varPsi$ is defined by

$$\begin{aligned} \varPsi(\mathbf{x}) &= \frac{1}{4\pi} \int _{V_{1}}^{}{|\mathbf{x}- \mathbf{y}|}\,d{\mathbf{y}}, \end{aligned}$$

(5.5)

which satisfies (see for example Kellogg [58])

$$\begin{aligned} \nabla^{4}\varPsi&= -2\nabla^{2}\varGamma= -2 \chi^{1}(\mathbf{x}). \end{aligned}$$

(5.6)

Once $P_{ijk\ell}$ is determined, the components of the Eshelby tensor can be calculated from (1.2) and the associated Concentration tensor can be found from (3.23). Recall that no assumptions have been made regarding the anisotropy of the inhomogeneity. This is not required in order for the $\mathbf{P}$-tensor to be determined. The only aspects of the inhomogeneity that influence the $\mathbf{P}$-tensor are its shape and, for anisotropic host phases, its orientation with respect to the axes of anisotropy of the host phase.

5.1.1 Sphere in an Isotropic Host Phase

Assume that the host phase is isotropic with elastic modulus tensor given in (5.1) and consider the case where $V_{1}$ is a sphere. The (uniform) tensors $\partial ^{2}\varGamma/\partial x_{i}\partial x_{j}$ and $\partial ^{4}\varPsi/\partial x_{i}\partial x_{j}\partial x_{k}\partial x_{\ell}$ will be spherically symmetric, i.e., isotropic and must possess full symmetry with respect to the interchange of any index. As with the potential problem therefore,

$$\begin{aligned} \frac{\partial^{2}{\varGamma}}{\partial{x_{i}}\partial{x_{j}}} &= \frac{1}{3}\delta_{ij} \end{aligned}$$

(5.7)

and for $\varPsi$ the general isotropic form, with the additional constraint regarding symmetry with respect to interchange of any indices, is

$$\begin{aligned} \frac{\partial ^{4}\varPsi}{ \partial x_{i}\partial x_{j}\partial x_{k}\partial x_{\ell}} &= \psi (\delta_{ij}\delta_{k\ell}+ \delta_{ik}\delta_{j\ell}+\delta_{i\ell} \delta_{jk}), \end{aligned}$$

(5.8)

where $\psi$ is a constant to be determined. Performing the contractions $j=i$, $\ell=k$ in (5.8) and utilizing (5.6), $\psi=-2/15$. Using (5.7) and (5.8) in (5.4), the components of the $\mathbf{P}$-tensor are

$$\begin{aligned} P_{ijk\ell} &= \frac{1}{6\mu_{0}} (\delta_{ik} \delta_{j\ell}+ \delta_{i\ell}\delta_{jk} ) - \frac{1}{30\mu_{0}(1-\nu_{0})} (\delta_{ij} \delta_{k\ell}+ \delta_{ik}\delta_{j\ell }+\delta_{i\ell} \delta_{jk} ). \end{aligned}$$

(5.9)

After simplification and writing in terms of the tensors $I_{ijk\ell }^{1}$ and $I_{ijk\ell}^{2}$ this becomes

$$\begin{aligned} P_{ijk\ell} &= p_{1} I_{ijk\ell}^{1} + p_{2} I_{ijk\ell}^{2}, \end{aligned}$$

(5.10)

where upon using the relationship (5.2) the components of the $\mathbf{P}$-tensor are determined as

$$\begin{aligned} p_{1} = \frac{1-2\nu_{0}}{6\mu_{0}(1-\nu_{0})} = \frac{1}{3\kappa _{0}+4\mu _{0}}, \qquad p_{2} = \frac{4-5\nu_{0}}{15\mu_{0}(1-\nu_{0})} = \frac{3(\kappa _{0}+2\mu _{0})}{5\mu _{0}(3\kappa _{0}+4\mu_{0})}. \end{aligned}$$

(5.11)

From this form the Eshelby tensor is easily determined via (1.2) as

$$\begin{aligned} S_{ijk\ell} &= s_{1} I_{ijk\ell}^{1} + s_{2} I_{ijk\ell}^{2}, \end{aligned}$$

(5.12)

where

$$\begin{aligned} s_{1} = \frac{1+\nu_{0}}{3(1-\nu_{0})} = \frac{3\kappa _{0}}{3\kappa _{0}+4\mu_{0}}, \qquad s_{2} = \frac{2(4-5\nu_{0})}{15(1-\nu_{0})} = \frac{6(\kappa _{0}+2\mu _{0})}{5(3\kappa _{0}+4\mu_{0})}. \end{aligned}$$

(5.13)

Isotropic Sphere

The $\mathbf{P}$-tensor derived in (5.10) holds for a spherical inhomogeneity of arbitrary anisotropy, embedded inside an isotropic host phase. If the inhomogeneity is also assumed isotropic with elastic modulus tensor of the form

$$\begin{aligned} C_{ijk\ell}^{1} &= 3\kappa _{1} I_{ijk\ell}^{1}+2 \mu_{1} I_{ijk\ell}^{2}, \end{aligned}$$

(5.14)

one can use (C.9) and the contraction properties of the tensors $I_{ijk\ell}^{1}, I_{ijk\ell}^{2}$ as defined in Appendix C.2 in order to write

$$\begin{aligned} \tilde{\mathcal{A}}_{ijk\ell} &= I_{ijk\ell}^{1}+I_{ijk\ell}^{2} + P_{ijmn} \bigl(C_{mnk\ell}^{1}-C_{mnk\ell}^{0} \bigr) \\ &= \bigl(1+3(\kappa _{1}-\kappa _{0})p_{1} \bigr)I_{ijk\ell}^{1} + \bigl(1+2(\mu_{1}-\mu _{0})p_{2} \bigr)I_{ijk\ell}^{2}. \end{aligned}$$

(5.15)

Finally the inversion properties of such a tensor are employed (as described in Appendix C.2.3) together with (5.11) in order to obtain the strain Concentration tensor:

$$\begin{aligned} \mathcal{A}_{ijk\ell} &= \biggl(\frac{3\kappa _{0}+4\mu_{0}}{3\kappa _{1}+4\mu _{0}} \biggr)I_{ijk\ell}^{1} + \frac{5\mu_{0}(3\kappa _{0}+4\mu_{0})}{3\kappa _{0}(3\mu_{0}+2\mu_{1})+4\mu _{0}(2\mu_{0}+3\mu_{1})}I_{ijk\ell}^{2}. \end{aligned}$$

(5.16)

Taking $\kappa _{1}, \mu_{1}\rightarrow0$ in (5.16) the Concentration tensor for a spherical cavity is obtained purely in terms of the host Poisson ratio:

$$\begin{aligned} \mathcal{A}_{ijk\ell} &= \frac{3(1-\nu_{0})}{2(1-2\nu _{0})}I^{1}_{ijk\ell} + \frac{15(1-\nu_{0})}{7-5\nu_{0}}I^{2}_{ijk\ell}. \end{aligned}$$

(5.17)

Cubic Sphere

Assume now that the sphere has cubic symmetry with elastic modulus tensor

$$\begin{aligned} C_{ijk\ell}^{1} &= 3\kappa _{1} I_{ijk\ell}^{1}+2 \mu_{1} I_{ijk\ell}^{2} + \eta_{1}\delta _{ijk\ell}, \end{aligned}$$

(5.18)

where

$$\begin{aligned} \delta_{ijk\ell} &= \textstyle\begin{cases} 1, & i=j=k=\ell,\\ 0, & \mbox{otherwise}. \end{cases}\displaystyle \end{aligned}$$

(5.19)

Instead of (5.15), it is found that

$$\begin{aligned} \tilde{\mathcal{A}}_{ijk\ell} = \bigl(1+3(\kappa _{1}- \kappa _{0})p_{1} \bigr)I_{ijk\ell}^{1} + \bigl(1+2(\mu_{1}-\mu_{0})p_{2} \bigr)I_{ijk\ell}^{2} + \eta_{1} P_{ijmn} \delta_{mnk\ell}. \end{aligned}$$

(5.20)

Next, using the form of the $\mathbf{P}$-tensor in (5.10) and the relationships (C.19), the expression (5.20) becomes

$$\begin{aligned} \tilde{\mathcal{A}}_{ijk\ell} = \tilde{\alpha }_{1} I_{ijk\ell}^{1} + \tilde{\alpha }_{2} I_{ijk\ell}^{2} + \tilde{\alpha }_{3}\delta_{ijk\ell}, \end{aligned}$$

(5.21)

where

$$\begin{aligned} \tilde{\alpha }_{1} = 1+3(\kappa _{1}-\kappa _{0})p_{1}+ \eta_{1}(p_{1}-p_{2}), \qquad\tilde{ \alpha }_{2} = 1+2(\mu_{1}-\mu_{0})p_{2}, \qquad\tilde{\alpha }_{3} = \eta_{1}p_{2}. \end{aligned}$$

(5.22)

Finally, appealing to the theory regarding cubic tensors in Appendix C.2.2 the Concentration tensor is

$$\begin{aligned} \mathcal{A}_{ijk\ell} &=\alpha _{1} I_{ijk\ell}^{1} + \alpha _{2} I_{ijk\ell}^{2} + \alpha _{3} \delta_{ijk\ell}, \end{aligned}$$

(5.23)

where

$$\begin{aligned} \alpha _{1} = \frac{\tilde{\alpha }_{2}^{2}+\tilde{\alpha }_{1}\tilde{\alpha }_{3}}{\tilde{\alpha }_{2}^{2}(\tilde{\alpha }_{1}+\tilde{\alpha }_{3})}, \qquad \alpha _{2} = \frac{1}{\tilde{\alpha }_{2}}, \qquad \alpha _{3} = -\frac{\tilde{\alpha }_{3}}{\tilde{\alpha }_{2}^{2}}. \end{aligned}$$

(5.24)

Note that when $\eta_{1}=0$ the case of an isotropic sphere is recovered as in (5.15).

5.1.2 Circular Cylinder in an Isotropic Host Phase

Now assume that $V_{1}$ is a circular cylinder with axis of symmetry in the $x_{3}$ direction. As for the potential problem, $\varGamma(\mathbf{x})$ must be independent of $x_{3}$ and isotropic in the $x_{1}x_{2}$ plane. Hence, as described in (4.13)

$$\begin{aligned} \frac{\partial^{2}{\varGamma}}{\partial{x_{i}}\partial{x_{j}}} &= \frac{1}{2}\varTheta_{ij}. \end{aligned}$$

(5.25)

Furthermore $\varPsi$ must also be independent of $x_{3}$ and be isotropic in the $x_{1}x_{2}$ plane and hence write

$$\begin{aligned} \frac{\partial \varPsi}{ \partial x_{i}\partial x_{j}\partial x_{k} \partial x_{\ell}} &= \psi (\varTheta_{ij}\varTheta_{k\ell}+ \varTheta_{ik}\varTheta_{j\ell}+\varTheta_{i\ell} \varTheta_{jk}), \end{aligned}$$

(5.26)

for some constant $\psi$, where the fact that this tensor should be fully symmetric with respect to interchange of any of its indices has been used. Performing the contractions $j=i$, $\ell=k$ and using (5.6) leads to the conclusion that $\psi=-\frac{1}{4}$. From (5.4) therefore

$$\begin{aligned} P_{ijk\ell} &= \dfrac{1}{8\mu_{0}} (\varTheta_{j\ell} \delta_{ik}+\varTheta_{jk}\delta_{i\ell}+ \varTheta_{i\ell}\delta_{jk}+\varTheta_{ik} \delta_{j\ell} ) \\ &\quad{}-\dfrac{1}{16\mu_{0}(1-\nu_{0})}(\varTheta_{ij}\varTheta _{k\ell}+ \varTheta_{ik}\varTheta_{j\ell}+\varTheta_{i\ell} \varTheta_{jk}). \end{aligned}$$

(5.27)

By writing $\delta_{ij}=\varTheta_{ij}+\delta_{i3}\delta_{j3}$ in the first term of (5.27) and then recognizing the appropriate Hill transversely isotropic basis tensors (the reader is referred to the discussion of this tensor basis in Appendix C.2.3) that arise as a result of the various contraction terms, the $\mathbf{P}$-tensor can be written in the form

$$\begin{aligned} P_{ijk\ell} &= \frac{1}{8\mu_{0}} \biggl[\frac{2(1-2\nu_{0})}{1-\nu _{0}}\mathcal {H}_{ijk\ell}^{1} + \frac{3-4\nu_{0}}{1-\nu_{0}}\mathcal{H}_{ijk\ell }^{5}+2 \mathcal{H}_{ijk\ell}^{6} \biggr] \end{aligned}$$

(5.28)

$$\begin{aligned} &= \frac{1}{4\mu_{0}} \biggl[\frac{2\mu_{0}}{\lambda _{0}+2\mu_{0}}\mathcal {H}_{ijk\ell}^{1} + \frac{\lambda _{0}+3\mu_{0}}{\lambda _{0}+2\mu_{0}}\mathcal{H}_{ijk\ell }^{5}+ \mathcal{H}_{ijk\ell}^{6} \biggr]. \end{aligned}$$

(5.29)

In order to derive the Eshelby tensor one contracts $\mathbf{P}$ with the host modulus tensor $\mathbf{C}^{0}$. To do this, write the isotropic basis tensors in terms of the Hill basis tensors using the expressions (C.26)–(C.28) and then use the contractions summarized in Table 2 in order to determine that

$$\begin{aligned} S_{ijk\ell} &= \frac{1}{4(1-\nu_{0})} \bigl(2\mathcal{H}^{1}_{ijk\ell} + 2\nu_{0} \mathcal{H}^{2}_{ijk\ell} +(3-4 \nu_{0})\mathcal{H}^{5}_{ijk\ell}+2(1-\nu_{0}) \mathcal{H}^{6}_{ijk\ell} \bigr). \end{aligned}$$

(5.30)

It may well be the case that a different tensor basis should be used if it transpires that the cylinder itself is more anisotropic than transversely isotropic, for use in Concentration tensors for example. However for computation, the matrix formulation, as described in Appendix C.4 can be of great utility when anisotropic basis tensors are becoming rather cumbersome.

Table 2 Short-hand notation for contractions of the basis tensors $\mathcal {H}_{ijkl}^{n}$

Suppose now that the cylinder is isotropic with elastic modulus tensor as defined in (5.14). Once again using (C.26)–(C.28), constructing $\tilde{A}_{ijk\ell}$ is then just a matter of exploiting the contractions in Table 2. It transpires that $\tilde{\mathcal{A}}_{ijk\ell}$ takes the form

$$\begin{aligned} \tilde{\mathcal{A}}_{ijk\ell} &= \sum_{n=1}^{6} \tilde{\alpha }_{n}\mathcal{H}^{n}_{ijk\ell}, \end{aligned}$$

(5.31)

where

$$\begin{aligned} \tilde{\alpha }_{1} &= \frac{\lambda _{1}+\mu_{1}+\mu_{0}}{(\lambda _{0}+2\mu_{0})}, \qquad\tilde{\alpha }_{2} = \frac{(\lambda _{1}-\lambda _{0})}{2(\lambda _{0}+2\mu_{0})}, \qquad\qquad\; \tilde{\alpha }_{3} = 0, \end{aligned}$$

(5.32)

$$\begin{aligned} \tilde{\alpha }_{4} &= 1, \qquad\qquad\qquad\quad\; \tilde{ \alpha }_{5} = \frac{\mu_{0}+\mu_{1}(3-4\nu_{0})}{4\mu_{0}(1-\nu_{0})}, \qquad \tilde{\alpha }_{6} = \frac{(\mu_{1}+\mu_{0})}{2\mu_{0}}. \end{aligned}$$

(5.33)

Using the inversion expression for transversely isotropic tensors stated in Appendix C.2.3, the appropriate Concentration tensor is thus derived as

$$\begin{aligned} \mathcal{A}_{ijk\ell} &= \sum_{n=1}^{6} \alpha _{n} \mathcal{H}_{ijk\ell}^{n}, \end{aligned}$$

(5.34)

where

$$\begin{aligned} \alpha _{1} &= \frac{\lambda _{0}+2\mu_{0}}{(\lambda _{1}+\mu_{1}+\mu_{0})}, \qquad \alpha _{2} = \frac{(\lambda _{0}-\lambda _{1})}{2(\lambda _{1}+\mu_{1}+\mu_{0})}, \qquad\ \; \alpha _{3} = 0, \end{aligned}$$

(5.35)

$$\begin{aligned} \alpha _{4} &= 1, \qquad\qquad \qquad \quad\ \ \ \; \alpha _{5} = \frac{4\mu_{0}(1-\nu_{0})}{\mu_{0}+\mu_{1}(3-4\nu _{0})}, \qquad \alpha _{6} = \frac{2\mu_{0}}{\mu_{1}+\mu_{0}}. \end{aligned}$$

(5.36)

Since $\alpha _{2}\neq \alpha _{3}$ it is seen that $\mathcal{A}_{1133}\neq \mathcal {A}_{3311}$. The circular cylindrical cavity limit can be obtained by setting $\mu_{1}=0$. Using the expression $\lambda =2\mu\nu/(1-2\nu)$,

$$\begin{aligned} \alpha _{1} &= \frac{2(1-\nu_{0})}{(1-2\nu_{0})}, \qquad \alpha _{2} = \frac{\nu_{0}}{1-2\nu _{0}}, \qquad\; \ \alpha _{3} = 0, \end{aligned}$$

(5.37)

$$\begin{aligned} \alpha _{4} &= 1, \qquad\qquad\quad\ \ \alpha _{5} = 4(1- \nu_{0}), \qquad \alpha _{6} = 2. \end{aligned}$$

(5.38)

Finally note that the fourth order tensor (5.34) can be represented in matrix form (see Appendix C.4 where the following “square bracket” notation is introduced) as

$$\begin{aligned}{} [\mathcal{A}] &= \left ( \textstyle\begin{array}{c@{\quad}c@{\quad}c@{\quad}c@{\quad}c@{\quad}c} a_{11} & a_{12} & a_{13} & 0 & 0 & 0 \\ a_{12} & a_{11} & a_{13} & 0 & 0 & 0 \\ a_{31} & a_{31} & a_{33} & 0 & 0 & 0 \\ 0 & 0 & 0 & a_{33} & 0 & 0 \\ 0 & 0 & 0 & 0 & a_{33} & 0 \\ 0 & 0 & 0 & 0 & 0 & a_{66} \end{array}\displaystyle \right ), \end{aligned}$$

(5.39)

where

$$\begin{aligned} a_{11} &= \frac{(1-\nu_{0})(3-4\nu_{0})}{(1-2\nu_{0})}, \qquad a_{12} = \frac{(\nu _{0}-1)(1-4\nu_{0})}{(1-2\nu_{0})}, \qquad a_{13} = \frac{\nu_{0}}{1-2\nu_{0}}, \end{aligned}$$

(5.40)

$$\begin{aligned} a_{31} &= 0, \qquad \qquad \qquad \qquad \ \ \;\, a_{33} = 1, \qquad\qquad \qquad \qquad\ \ \ a_{66} = 2(1- \nu_{0}). \end{aligned}$$

(5.41)

Clearly it is possible to write down explicit expressions for the Concentration tensor when the circular cylinder is anisotropic. This merely complicates the tensorial (or matrix) operations after the derivation of the $\mathbf{P}$-tensor in (5.27). Given this $\mathbf{P}$-tensor, perhaps the most important aspect is then to choose the tensor basis set correctly, given the anisotropy of the inhomogeneity. For example, when the cylinder itself is transversely isotropic (a common occurrence in applications) it is considered sensible to use a transversely isotropic tensor basis set. For practical purposes and especially for the sake of computation, using the matrix formulation of tensors is advantageous in cases where the tensor basis sets become rather cumbersome. The following procedure is used, in the usual notation, referring to Appendix C.4, and defining the $6\times6$ matrix $[P]$ associated with the $\mathbf{P}$-tensor, define

$$\begin{aligned}{} [\tilde{\mathcal{A}}] &= [I] + [P] [W] \bigl[C^{1}-C^{0} \bigr], \end{aligned}$$

(5.42)

where $[W]$ is defined in (C.56) and therefore

$$\begin{aligned}{} [\mathcal{A}] &= [W]^{-1}[\tilde{\mathcal{A}}]^{-1}[W]^{-1}. \end{aligned}$$

(5.43)

5.1.3 Spheroid in an Isotropic Host Phase

When $V_{1}$ is a spheroid, the potential theory outlined in Appendix B is once again of use. It is clear that the $\mathbf{P}$-tensor must be transversely isotropic and therefore will take the form

$$\begin{aligned} P_{ijk\ell} &= \sum_{n=1}^{6} p_{n} \mathcal{H}_{ijk\ell}^{n}. \end{aligned}$$

(5.44)

The separate contributions to the $\mathbf{P}$-tensor shall therefore first be decomposed into this form. Firstly, from the potential case described in Sect. 4.1.4

$$\begin{aligned} \frac{\partial^{2}{\varGamma}}{\partial{x_{i}}\partial{x_{j}}} &= \gamma \varTheta_{ij} + \gamma _{3} \delta_{i3}\delta_{j3}, \end{aligned}$$

(5.45)

where $\gamma _{3} = \mathcal{S}(\varepsilon )$ and $\gamma =\frac{1}{2}(1-\gamma _{3})$. Representing all terms in the Hill basis one can show that

$$\begin{aligned} &\frac{1}{4} \biggl(\frac{\partial^{2}{\varGamma}}{\partial {x_{j}}\partial{x_{\ell}}}\delta_{ik}+ \frac{\partial^{2}{\varGamma }}{\partial{x_{j}}\partial{x_{k}}}\delta_{i{\ell}} +\frac{\partial ^{2}{\varGamma}}{\partial{x_{i}}\partial{x_{\ell}}} \delta_{jk}+\frac {\partial^{2}{\varGamma}}{\partial{x_{i}}\partial{x_{k}}}\delta _{j{\ell}} \biggr) \\ &\quad= \gamma \bigl(\mathcal{H}^{1}_{ijk{\ell}}+ \mathcal{H}^{5}_{ijk{\ell}} \bigr) + \gamma _{3} \mathcal{H}^{4}_{ijk{\ell}}+ \frac{1}{2}(\gamma + \gamma _{3})\mathcal{H}_{ijk{\ell }}^{6}. \end{aligned}$$

(5.46)

This is seen by using $\varTheta_{ik}=\delta_{ik}-\delta_{i3}\delta _{k3}$ and writing for example,

$$\begin{aligned} \frac{\partial^{2}{\varGamma}}{\partial{x_{j}}\partial{x_{\ell }}}\delta_{ik} &= \frac{\partial^{2}{\varGamma}}{\partial{x_{j}}\partial{x_{\ell}}}( \varTheta_{ik}+\delta_{i3}\delta_{k3}) \\ &= (\gamma \varTheta_{j\ell} + \gamma _{3}\delta_{j3} \delta_{\ell3}) (\varTheta_{ik}+\delta_{i3} \delta_{k3}) \\ &= \gamma (\varTheta_{j\ell}\varTheta_{ik} + \varTheta_{j\ell} \delta_{i3}\delta_{k3}) + \gamma _{3}( \varTheta_{ik} \delta_{j3}\delta_{\ell3}+ \delta_{i3}\delta_{j3} \delta_{k3} \delta_{\ell3}). \end{aligned}$$

(5.47)

Doing this for each term on the left hand side of (5.46), combining and using the definitions of the transversely isotropic basis tensors in Appendix C.2.3 leads to the form on the right hand side of (5.46).

Further, after much algebraic manipulation using the simplifications of the integrals in Appendix B in the case of spheroids one can show that

$$\begin{aligned} \frac{1}{4}\frac{\partial ^{4}\varPsi}{ \partial x_{i}\partial x_{j}\partial x_{k}\partial x_{\ell}} &= \sum_{n=1}^{6} \psi_{n} \mathcal{H}_{ijk{\ell}}^{n} , \end{aligned}$$

(5.48)

where, upon using $\gamma _{3}=1-2\gamma$,

$$\begin{aligned} \psi_{1} &= \frac{\varepsilon ^{2}(4\gamma -1)-\gamma }{4(1-\varepsilon ^{2})}, \qquad \psi_{2} = \psi_{3} = \frac{\varepsilon ^{2}(1-2\gamma )-\gamma }{4(1-\varepsilon ^{2})}, \end{aligned}$$

(5.49)

$$\begin{aligned} \psi_{4} &= \frac{3 \gamma -1}{2(1-\varepsilon ^{2})}, \qquad\psi_{5} = \frac{1}{2}\psi_{1}, \qquad\psi_{6} = 2 \psi_{2}. \end{aligned}$$

(5.50)

Therefore

$$\begin{aligned} p_{1} &= \frac{1}{\mu_{0}} \biggl(\gamma +\frac{1}{(1-\nu_{0})} \psi_{1} \biggr), \qquad p_{2} = p_{3} = \frac{\psi_{2}}{\mu_{0}(1-\nu_{0})}, \end{aligned}$$

(5.51)

$$\begin{aligned} p_{4} &= \frac{1}{\mu_{0}} \biggl(1-2\gamma +\frac{1}{(1-\nu_{0})} \psi_{4} \biggr), \qquad p_{5} = \frac{1}{\mu_{0}} \biggl( \gamma +\frac{1}{2(1-\nu_{0})}\psi_{1} \biggr), \end{aligned}$$

(5.52)

$$\begin{aligned} p_{6} &= \frac{1}{\mu_{0}} \biggl(\frac{1}{2}(1-\gamma )+ \frac{2}{(1-\nu_{0})}\psi_{2} \biggr). \end{aligned}$$

(5.53)

A good check is to ascertain that the result for a sphere in an isotropic host phase is recovered by taking $\varepsilon \rightarrow1$ and using (4.31). This is easily done and yields the result derived in Sect. 5.1.1 associated with a sphere. Another useful limit is to take $\varepsilon \rightarrow0$. Results already derived can be employed, e.g., (4.39), which when used in (5.49)–(5.50) yield

$$\begin{aligned} &\psi_{1} = -\frac{\pi}{16}\varepsilon +O \bigl(\varepsilon ^{3} \bigr), \qquad\psi_{2} = \psi_{3} = -\frac{\pi}{16} \varepsilon + \frac{\varepsilon ^{2}}{2} + O \bigl(\varepsilon ^{3} \bigr), \end{aligned}$$

(5.54)

$$\begin{aligned} & \begin{aligned} &\psi_{4} = - \frac{1}{2}+\frac{3\pi}{8}\varepsilon - 2\varepsilon ^{2} + O \bigl( \varepsilon ^{3} \bigr), \qquad\psi_{5} = -\frac{\pi}{32} \varepsilon + O \bigl(\varepsilon ^{3} \bigr), \qquad \\ &\psi_{6} = -\frac{\pi }{8}\varepsilon + \varepsilon ^{2} + O \bigl(\varepsilon ^{3} \bigr). \end{aligned} \end{aligned}$$

(5.55)

For the components of the $\mathbf{P}$-tensor this then gives

$$\begin{aligned} p_{1} &= \frac{1}{\mu_{0}} \biggl(\frac{\pi(3-4\nu_{0})}{16(1-\nu _{0})}\varepsilon - \varepsilon ^{2} \biggr)+O \bigl(\varepsilon ^{3} \bigr), \end{aligned}$$

(5.56)

$$\begin{aligned} p_{2} &= p_{3} = \frac{1}{\mu_{0}} \biggl(- \frac{\pi}{16(1-\nu_{0})}\varepsilon + \frac{1}{2(1-\nu_{0})}\varepsilon ^{2} \biggr) + O \bigl(\varepsilon ^{3} \bigr), \end{aligned}$$

(5.57)

$$\begin{aligned} p_{4} &= \frac{1}{\mu_{0}} \biggl(\frac{1-2\nu_{0}}{2(1-\nu_{0})} - \frac{\pi (1-4\nu_{0})}{8(1-\nu_{0})}\varepsilon - \frac{2\nu_{0}}{1-\nu_{0}}\varepsilon ^{2} \biggr) + O \bigl(\varepsilon ^{3} \bigr), \end{aligned}$$

(5.58)

$$\begin{aligned} p_{5} &= \frac{1}{\mu_{0}} \biggl(\frac{\pi(7-8\nu_{0})}{32(1-\nu _{0})}\varepsilon - \varepsilon ^{2} \biggr) + O \bigl(\varepsilon ^{3} \bigr), \end{aligned}$$

(5.59)

$$\begin{aligned} p_{6} &= \frac{1}{\mu_{0}} \biggl(\frac{1}{2} - \frac{\pi(2-\nu_{0})}{8(1-\nu _{0})}\varepsilon + \frac{3-\nu_{0}}{2(1-\nu_{0})}\varepsilon ^{2} \biggr) + O \bigl(\varepsilon ^{3} \bigr). \end{aligned}$$

(5.60)

The Eshelby tensor with respect to a transversely isotropic basis, in the form

$$\begin{aligned} S_{ijk\ell} &= \sum_{n=1}^{6} s_{n} \mathcal{H}^{n}_{ijk\ell} \end{aligned}$$

(5.61)

for either the spheroid with components of $\mathbf{P}$-tensor (5.51)–(5.53) or the $\varepsilon \rightarrow0$ limit of the spheroid with components of the $\mathbf{P}$-tensor (5.56)–(5.60) has components that are related directly to the components of the $\mathbf{P}$-tensor via the expressions

$$\begin{aligned} s_{1} &= 2\mu_{0} \biggl(\frac{p_{1}+2 p_{2} \nu_{0}}{1-2\nu_{0}} \biggr), \qquad s_{2} = 2\mu_{0} \biggl(\frac{p_{1}\nu_{0}+(1-\nu_{0})p_{2}}{1-2\nu_{0}} \biggr), \end{aligned}$$

(5.62)

$$\begin{aligned} s_{3} &= 2\mu_{0} \biggl(\frac{p_{3} +\nu_{0} p_{4}}{1-2\nu_{0}} \biggr), \qquad \ \ s_{4} = 2\mu_{0} \biggl(\frac{2p_{3}\nu_{0}+(1-\nu_{0})p_{4}}{1-2\nu_{0}} \biggr), \end{aligned}$$

(5.63)

$$\begin{aligned} s_{5} &= 2\mu_{0}p_{5}, \qquad\qquad\qquad\quad \, s_{6} = 2\mu_{0}p_{6}. \end{aligned}$$

(5.64)

Finally, it is straightforward, but rather tedious, to show, using the $\mathbf{P}$-tensor derived in the previous example, that the strain Concentration tensor associated with an isotropic spheroid embedded in an isotropic host phase is

$$\begin{aligned} \mathcal{A}_{ijk\ell} &= \sum_{n=1}^{6} \alpha _{n} \mathcal{H}_{ijk\ell}^{n}, \end{aligned}$$

(5.65)

where with $\Delta=\frac{1}{2}q_{1}q_{4}-q_{2}q_{3}$,

$$\begin{aligned} \alpha _{1} &= \frac{q_{4}}{2\Delta}, \qquad \alpha _{2} = - \frac{q_{2}}{2\Delta}, \qquad \alpha _{3} = -\frac{q_{3}}{2\Delta}, \end{aligned}$$

(5.66)

$$\begin{aligned} \alpha _{4} &= \frac{q_{1}}{2\Delta}, \qquad \alpha _{5} = \frac{1}{q_{5}}, \qquad\quad\; \alpha _{6} = \frac {1}{q_{6}} \end{aligned}$$

(5.67)

and where

$$\begin{aligned} q_{1} &= 1 + 2p_{1} (\lambda _{d}+\mu_{d} ) + 2p_{2}\lambda _{d}, \qquad q_{2} = p_{1} \lambda _{d} + p_{2} (\lambda _{d}+2\mu_{d} ), \end{aligned}$$

(5.68)

$$\begin{aligned} q_{3} &= 2p_{3} (\lambda _{d}+\mu_{d} ) + p_{4}\lambda _{d}, \qquad\qquad\, q_{4} = 1 + 2p_{3} \lambda _{d} + p_{4} (\lambda _{d}+2 \mu_{d} ), \end{aligned}$$

(5.69)

$$\begin{aligned} q_{5} &= 1 + 2p_{5}\mu_{d}, \qquad\qquad\qquad \qquad\, \; \; q_{6} = 1 + 2p_{6}\mu_{d}, \end{aligned}$$

(5.70)

with

$$\begin{aligned} \lambda _{d}=\lambda _{1}-\lambda _{0}, \qquad\mu_{d} &= \mu_{1}-\mu_{0}. \end{aligned}$$

(5.71)

The spheroidal cavity result is simply (5.65)–(5.71) with $\lambda _{1}=\mu_{1}=0$ and so every occurrence of $\lambda _{d}$ and $\mu_{d}$ is simply replaced with $-\lambda _{0}$ and $-\mu_{0}$ respectively. The result can be obtained in terms of $\nu_{0}$ alone by using the expression $\lambda _{0}=2\mu_{0}\nu_{0}/(1-2\nu_{0})$.

The average of the Concentration tensor over uniform orientations of spheroids can be obtained using the result in (C.52) in order to derive an expression of the form

$$\begin{aligned} \underline{\mathcal{A}}_{ijk\ell} &= \alpha _{1} I^{1}_{ijk\ell} + \alpha _{2} I^{2}_{ijk\ell}. \end{aligned}$$

(5.72)

5.1.4 Elastic Layer

The result for the spheroid can be used in order to determine the $\mathbf{P}$-tensor and Concentration tensor for an elastic layer, taking $\varepsilon =0$ in (5.56)–(5.60),

$$\begin{aligned} P_{ijk\ell} &= \frac{1}{2\mu_{0}} \biggl(\frac{1-2\nu_{0}}{1-\nu_{0}} \mathcal{H}^{4}_{ijk\ell } + \mathcal{H}^{6}_{ijk\ell} \biggr), \end{aligned}$$

(5.73)

and Eshelby’s tensor easily follows as

$$\begin{aligned} S_{ijk\ell} &= \frac{\nu_{0}}{1-\nu_{0}}\mathcal{H}^{3}_{ijk\ell} + \mathcal{H}^{4}_{ijk\ell} + \mathcal{H}^{6}_{ijk\ell}. \end{aligned}$$

(5.74)

Using the $\mathbf{P}$-tensor the Concentration tensor for an isotropic layer is straightforwardly determined as

$$\begin{aligned} \mathcal{A}_{ijk\ell} = \mathcal{H}^{1}_{ijk\ell} + \biggl( \frac{\lambda _{0}-\lambda _{1}}{\lambda _{1}+2\mu_{1}} \biggr)\mathcal{H}^{3}_{ijk\ell} + \biggl( \frac{\lambda _{0}+2\mu_{0}}{\lambda _{1}+2\mu_{1}} \biggr)\mathcal{H}^{4}_{ijk\ell } + \mathcal{H}^{5}_{ijk\ell} + \frac{\mu_{0}}{\mu_{1}} \mathcal{H}^{6}_{ijk\ell}. \end{aligned}$$

(5.75)

5.1.5 Limiting Case of a Penny-Shaped Crack

For a penny shaped crack, terms up to $O(\varepsilon )$ are retained in (5.56)–(5.60) to obtain

$$\begin{aligned} P_{ijk\ell} &= \frac{1}{\mu_{0}} \biggl[\frac{\pi(3-4\nu _{0})}{16(1-\nu _{0})}\varepsilon \mathcal{H}^{1}_{ijk\ell} -\frac{\pi}{16(1-\nu_{0})}\varepsilon \bigl( \mathcal{H}^{2}_{ijk\ell}+\mathcal{H}^{3}_{ijk\ell} \bigr) \\ &\quad{}+ \biggl(\frac{1-2\nu_{0}}{2(1-\nu_{0})} - \frac{\pi(1-4\nu _{0})}{8(1-\nu _{0})}\varepsilon \biggr) \mathcal{H}^{4}_{ijk\ell} \\ &\quad{}+\frac{\pi(7-8\nu_{0})}{32(1-\nu_{0})}\varepsilon \mathcal {H}^{5}_{ijk\ell} + \biggl( \frac{1}{2} - \frac{\pi(2-\nu_{0})}{8(1-\nu_{0})}\varepsilon \biggr) \mathcal{H}^{6}_{ijk\ell} \biggr]. \end{aligned}$$

(5.76)

Using this to determine the Concentration tensor with $\mu_{1}=0$ as in the cavity limit, it is straightforwardly shown that

$$\begin{aligned} \mathcal{A}_{ijk\ell} &= (1-\nu_{0})\mathcal{H}^{1}_{ijk\ell} -\frac {1}{2}(1-\nu_{0})\mathcal{H}^{2}_{ijk\ell} + \biggl(\frac{4\nu_{0}(1-\nu_{0})}{\pi(1-2\nu_{0})}\frac{1}{\varepsilon }-\frac {1}{2}(1- \nu_{0}) (1+2\nu_{0}) \biggr)\mathcal{H}^{3}_{ijk\ell} \\ &\quad{}+ \biggl(\frac{4(1-\nu_{0})^{2}}{\pi(1-2\nu_{0})}\frac {1}{\varepsilon }+\frac {1}{2}(1+2 \nu_{0}) (1-\nu_{0}) \biggr)\mathcal{H}^{4}_{ijk\ell}+ \mathcal{H}^{5}_{ijk\ell} \\ &\quad{}+ \biggl(\frac{4(1-\nu_{0})}{\pi(2-\nu_{0})}\frac{1}{\varepsilon }+\frac{16(3-\nu _{0})(1-\nu_{0})}{\pi^{2}(2-\nu_{0})^{2}} \biggr) \mathcal{H}^{6}_{ijk\ell}, \end{aligned}$$

(5.77)

where terms of $O(\varepsilon )$ have been neglected in (5.77). This expression has singular behaviour as $\varepsilon \rightarrow0$ akin to the potential problem result (4.40) and when deriving effective properties for distributions of cracks, this singular nature is necessary to yield the correct effective behaviour [48]. In fact although $O(1)$ coefficients have been retained in (5.77) only the singular terms are required in order to determine effective properties. The expression (5.77) corrects the typographical errors given on p. 104 of [81].

A common requirement is the determination of the effective properties of a medium comprising penny shaped cracks that are uniformly distributed and uniformly oriented inside some host material. Using (C.52) the associated Concentration tensor is shown to be

$$\begin{aligned} \underline{\mathcal{A}}_{ijk\ell} &= \frac{4(1-\nu_{0}^{2})}{3\pi \varepsilon (1-2\nu_{0})}I^{1}_{ijk\ell} + \frac{8(1-\nu_{0})(5-\nu_{0})}{15\pi \varepsilon (2-\nu _{0})}I^{2}_{ijk\ell} + O(1). \end{aligned}$$

(5.78)

The case of an ellipsoid in an isotropic medium shall now be considered. In order to deal with this generally in a tensor setting, ideally an orthotropic tensor basis should be used. Although it is possible to write down such a basis, details are rather lengthy and in fact for practical computation, it is perhaps most sensible to write down the nine independent components of the $\mathbf{P}$-tensor and use matrix computations in the manner described after Sect. 5.1.2 above.

5.1.6 Ellipsoid in an Isotropic Host Phase

The nine independent components of the $\mathbf{P}$-tensor for an ellipsoid can be defined in terms of the function $\mathcal{E}(\varepsilon _{n};\varepsilon _{1},\varepsilon _{2})$ and the semi-axes ratios $\varepsilon _{n}$.

The nine independent components of the Eshelby tensor for an ellipsoid in an isotropic medium are usually stated in terms of the four components $S_{1111}, S_{1122}, S_{1133}$ and $S_{1212}$ together with cyclic properties of the indices, in terms of $I_{mn}$ and $I_{m}$ as defined in (B.42)–(B.45). In turn these lead to expressions in terms of the fundamental integral $\mathcal{E}(\varepsilon _{n};\varepsilon _{1},\varepsilon _{2})$ via (B.35) and (B.43)–(B.46). Now use (5.4) with (4.21) and (B.48) and employ the properties (B.43)–(B.46) to derive the following compact forms

$$\begin{aligned} P_{1111} &= \frac{3}{16\pi\mu_{0}(1-\nu_{0})}I_{11} + \frac{1-4\nu_{0}}{16\pi \mu_{0}(1-\nu_{0})}I_{1}, \end{aligned}$$

(5.79)

$$\begin{aligned} P_{1122} &= \frac{1}{16\pi\mu_{0}(1-\nu_{0})}(I_{21}-I_{1}), \end{aligned}$$

(5.80)

$$\begin{aligned} P_{1133} &= \frac{1}{16\pi\mu_{0}(1-\nu_{0})}(I_{31}-I_{1}), \end{aligned}$$

(5.81)

$$\begin{aligned} P_{1212} &= \frac{1}{32\pi\mu_{0}(1-\nu_{0})}(I_{12}+I_{21})+ \frac{(1-2\nu _{0})}{32\pi\mu_{0}(1-\nu_{0})}(I_{1}+I_{2}). \end{aligned}$$

(5.82)

All other non-zero components are obtained by cyclic permutation of the indices in the above equations. Those components that cannot be obtained via cyclic permutation are zero, e.g., $P_{1112}=P_{1223}=P_{1323}=0$.

For representations and calculations of the Concentration tensor it is convenient to use the matrix representation of the tensors. This is discussed in the next section by considering the elliptical cylinder and ribbon crack limits. First however the components of the Eshelby tensor are stated, using (1.2) and noting the slightly modified notation for $I_{mn}$ in (B.42) (i.e., the factor of $a_{m}^{2}$) as compared with the standard definition in Mura [92]. The components are expressed as

$$\begin{aligned} S_{1111} &= \frac{3}{8\pi(1-\nu_{0})}I_{11} + \frac{1-2\nu_{0}}{8\pi(1-\nu _{0})}I_{1}, \end{aligned}$$

(5.83)

$$\begin{aligned} S_{1122} &= \frac{1}{8\pi(1-\nu_{0})}\frac{\varepsilon _{1}^{2}}{\varepsilon _{2}^{2}}I_{12} - \frac{1-2\nu_{0}}{8\pi(1-\nu_{0})}I_{1}, \end{aligned}$$

(5.84)

$$\begin{aligned} S_{1133} &= \frac{1}{8\pi(1-\nu_{0})}\frac{\varepsilon _{1}^{2}}{\varepsilon _{3}^{2}}I_{13} - \frac{1-2\nu_{0}}{8\pi(1-\nu_{0})}I_{1}, \end{aligned}$$

(5.85)

$$\begin{aligned} S_{1212} &= \frac{1+\varepsilon _{1}^{2}/\varepsilon _{2}^{2}}{16\pi(1-\nu _{0})}I_{12}+\frac {1-2\nu_{0}}{16\pi(1-\nu_{0})}(I_{1}+I_{2}) \end{aligned}$$

(5.86)

and permutation rules follow as for the $\mathbf{P}$-tensor.

5.1.7 Elliptical Cylinder and Ribbon-Crack Limit

In Sect. 4.1.5 it was shown that in the limit as $a_{3}\rightarrow\infty$, $\mathcal{E}(1;\varepsilon _{1},\varepsilon _{2}) \rightarrow 0$ and

$$\begin{aligned} \mathcal{E}(\varepsilon _{1};\varepsilon _{1},\varepsilon _{2}) \rightarrow\frac {a_{2}}{a_{1}+a_{2}} = \frac{\epsilon }{1+\epsilon }, \qquad\mathcal{E}( \varepsilon _{2};\varepsilon _{1},\varepsilon _{2}) &\rightarrow \frac {a_{1}}{a_{1}+a_{2}} = \frac{1}{1+\epsilon }, \end{aligned}$$

(5.87)

where $\epsilon =a_{2}/a_{1}$. These are used in the expressions for $I_{mn}$ and $I_{n}$ in Appendix B and substituted into (5.79)–(5.82) to determine the associated $\mathbf{P}$-tensor components. Since the $\mathbf{P}$-tensor is still orthotropic there are nine independent components:

$$\begin{aligned} P_{1111} &= \epsilon \biggl(\frac{4(1+\epsilon )(1-\nu_{0})-(1+2\epsilon )}{4\mu _{0}(1-\nu _{0})(1+\epsilon )^{2}} \biggr), \end{aligned}$$

(5.88)

$$\begin{aligned} P_{2222} &= \frac{4(1+\epsilon )(1-\nu_{0})-(2+\epsilon )}{4\mu_{0}(1-\nu _{0})(1+\epsilon )^{2}}, \end{aligned}$$

(5.89)

$$\begin{aligned} P_{3333} &= 0, \end{aligned}$$

(5.90)

$$\begin{aligned} P_{1122} &= \frac{-\epsilon }{4\mu_{0}(1-\nu_{0})(1+\epsilon )^{2}},\qquad P_{1133} = 0, \end{aligned}$$

(5.91)

$$\begin{aligned} P_{2233} &= 0, \qquad\qquad\qquad\qquad\qquad P_{1313} = \frac {\epsilon }{4\mu_{0}(1+\epsilon )}, \end{aligned}$$

(5.92)

$$\begin{aligned} P_{2323} &= \frac{1}{4\mu_{0}(1+\epsilon )}, \qquad\qquad\quad\ \ \; P_{1212} = \frac{(1-\nu _{0})(1+\epsilon )^{2}-\epsilon }{4\mu_{0}(1-\nu_{0})(1+\epsilon )^{2}}. \end{aligned}$$

(5.93)

The Eshelby tensor components follow as

$$\begin{aligned} S_{1111} &= \epsilon \biggl(\frac{1+2(1+\epsilon )(1-\nu_{0})}{2(1-\nu _{0})(1+\epsilon )^{2}} \biggr), \qquad S_{2222} = \frac{\epsilon +2(1-\nu_{0})(1+\epsilon )}{2(1-\nu_{0})(1+\epsilon )^{2}}, \end{aligned}$$

(5.94)

$$\begin{aligned} S_{3333} &= 0, \qquad \qquad \qquad \qquad \qquad \qquad S_{1133} = \frac{\nu_{0}\epsilon }{(1-\nu_{0})(1+\epsilon )}, \end{aligned}$$

(5.95)

$$\begin{aligned} S_{3311} &= 0, \qquad \qquad \qquad \qquad \qquad \qquad S_{1122} = \frac{-\epsilon +2\epsilon (1+\epsilon )\nu_{0}}{2(1-\nu_{0})(1+\epsilon )^{2}}, \end{aligned}$$

(5.96)

$$\begin{aligned} S_{2211} &= \frac{-\epsilon +2\nu_{0}(1+\epsilon )}{2(1-\nu_{0})(1+\epsilon )^{2}}, \qquad \qquad \quad\, S_{2233} =\frac{\nu_{0}}{(1-\nu_{0})(1+\epsilon )}, \end{aligned}$$

(5.97)

$$\begin{aligned} S_{3322} &= 0, \qquad \qquad \qquad \qquad \qquad \qquad S_{1313} = \frac{\epsilon }{2(1+\epsilon )}, \end{aligned}$$

(5.98)

$$\begin{aligned} S_{2323} &= \frac{1}{2(1+\epsilon )}, \qquad \qquad \qquad \qquad \ \ \;\, S_{1212} =\frac{-\epsilon +(1-\nu _{0})(1+\epsilon )^{2}}{2(1-\nu_{0})(1+\epsilon )^{2}}, \end{aligned}$$

(5.99)

noting that Eshelby’s tensor does not possess the major symmetry, unlike Hill’s tensor.

Suppose that the host and elliptical cylinder are both isotropic. In order to determine the Concentration tensor orthotropic tensors are required. Although it is possible to use an orthotropic basis set, it is perhaps most convenient to work with the matrix formulation of the tensors and derive the Concentration tensor using a symbolic mathematical package such as Mathematica. In doing this the matrix formulation $[\mathcal{A}]$ of the tensor $\boldsymbol{\mathcal{A}}$ is employed as noted in (5.42)–(5.43). The components of the matrix are so long that to list these here would not be beneficial but two very useful limits shall be written down. The elliptical cylindrical cavity limit is obtained by taking $\mu _{1}\rightarrow0$ which yields a matrix form of the tensor as

$$\begin{aligned}{} [\mathcal{A}] &= \left ( \textstyle\begin{array}{c@{\quad}c@{\quad}c@{\quad}c@{\quad}c@{\quad}c} a_{11} & a_{12} & a_{13} & 0 & 0 & 0 \\ a_{21} & a_{22} & a_{23} & 0 & 0 & 0 \\ a_{13} & a_{23} & a_{33} & 0 & 0 & 0 \\ 0 & 0 & 0 & a_{44} & 0 & 0 \\ 0 & 0 & 0 & 0 & a_{55} & 0 \\ 0 & 0 & 0 & 0 & 0 & a_{66} \end{array}\displaystyle \right ), \end{aligned}$$

(5.100)

where

$$\begin{aligned} a_{11} &= \frac{(1-\nu_{0})(1+2\epsilon -2(1+\epsilon )\nu_{0})}{(1-2\nu_{0})}, \qquad a_{12} = \frac{(1-\nu_{0})(-1+2(1+\epsilon )\nu_{0})}{(1-2\nu_{0})}, \end{aligned}$$

(5.101)

$$\begin{aligned} a_{13} &= \frac{(2\epsilon -1+2(1-\epsilon )\nu_{0})\nu_{0}}{(1-2\nu_{0})}, \qquad \qquad \ \; a_{21} = \frac{(1-\nu_{0})(-\epsilon +2(1+\epsilon )\nu_{0})}{\epsilon (1-2\nu_{0})}, \end{aligned}$$

(5.102)

$$\begin{aligned} a_{22} &= \frac{(1-\nu_{0})(2+\epsilon -2(1+\epsilon )\nu_{0})}{\epsilon (1-2\nu _{0})}, \qquad\ \ a_{23} = \frac{(2-\epsilon +2(\epsilon -1)\nu_{0})\nu_{0}}{\epsilon (1-2\nu_{0})}, \end{aligned}$$

(5.103)

$$\begin{aligned} a_{33} &= 1, \qquad\qquad \qquad\qquad \qquad \qquad \quad \ \ \ \, a_{44} = \frac{1+\epsilon }{2\epsilon }, \end{aligned}$$

(5.104)

$$\begin{aligned} a_{55} &= \frac{1+\epsilon }{2}, \qquad \qquad \qquad \qquad \qquad \qquad a_{66} = \frac{(1+\epsilon )^{2}(1-\nu_{0})}{2\epsilon } \end{aligned}$$

(5.105)

and taking the limit as $\epsilon \rightarrow0$ yields the ribbon-crack limit, retaining terms up to $O(1)$ in $\epsilon$,

$$\begin{aligned} a_{11} &= 1-\nu_{0}, \qquad\qquad\qquad\quad\ \ \; a_{12} = -(1- \nu_{0}), \end{aligned}$$

(5.106)

$$\begin{aligned} a_{13} &= -\nu_{0}, \qquad\qquad\qquad\qquad\ \ \thinspace a_{21} = \frac{2(1-\nu_{0})\nu_{0}}{(1-2\nu_{0})\epsilon }-(1-\nu_{0}), \end{aligned}$$

(5.107)

$$\begin{aligned} a_{22} &= \frac{2(1-\nu_{0})^{2}}{(1-2\nu_{0})\epsilon }+1-\nu_{0}, \qquad a_{23} = \frac{2(1-\nu_{0})\nu_{0}}{(1-2\nu_{0})\epsilon }-\nu_{0}, \end{aligned}$$

(5.108)

$$\begin{aligned} a_{33} &= 1, \qquad\qquad\qquad\qquad\quad\ \ \ a_{44} = \frac {1}{2\epsilon }+ \frac{1}{2}, \end{aligned}$$

(5.109)

$$\begin{aligned} a_{55} &= \frac{1}{2}, \qquad\qquad\qquad\qquad\quad\ \ a_{66} = \frac{1}{2\epsilon }(1-\nu_{0})+(1- \nu_{0}). \end{aligned}$$

(5.110)

5.1.8 Flat Ellipsoid

Consider the case when $a_{1}> a_{2}\gg a_{3}$. It is straightforward to take this limit in (B.36)–(B.40) in order to obtain

$$\begin{aligned} I_{1} &= 4\pi \varepsilon _{2}\frac{(F(k)-E(k)}{((\varepsilon _{2}/\varepsilon _{1})^{2}-1)}, \end{aligned}$$

(5.111)

$$\begin{aligned} I_{2} &= 4\pi \varepsilon _{2}E(k) - I_{1}, \end{aligned}$$

(5.112)

$$\begin{aligned} I_{3} &= 4\pi- 4\pi \varepsilon _{2} E(k), \end{aligned}$$

(5.113)

where with reference to (B.39) and (B.40), $F(k)$ and $E(k)$ are introduced as the complete Elliptic integrals of the first and second kind, respectively

$$\begin{aligned} E(k) &= \int_{0}^{\pi/2}\frac{dx}{(1-k^{2}\sin^{2} x)^{1/2}}, \qquad F(k) = \int_{0}^{\pi/2} \bigl(1-k^{2} \sin^{2} x \bigr)^{1/2}\, dx \end{aligned}$$

(5.114)

and $k=1-\frac{\varepsilon _{1}^{2}}{\varepsilon _{2}^{2}}$. From (5.111)–(5.113), the $I_{mn}$ can be straightforwardly determined via (B.43)–(B.46) and thus the components of the $\mathbf{P}$-tensor from (5.79)–(5.82) and (5.83)–(5.86).

5.1.9 Spheroid Limit Check

One can straightforwardly take the spheroidal inhomogeneity limit $a_{1}=a_{2}=a\neq a_{3}$ in the ellipsoidal result above. In particular it is noted that in the limit as $a_{1}\rightarrow a_{2}=a$, referring to Sect. 4.1.4,

$$\begin{aligned} I_{1} &= I_{2} = 2\pi \bigl(1-\mathcal{S}(\epsilon ) \bigr), \end{aligned}$$

(5.115)

$$\begin{aligned} I_{3} &= 4\pi\mathcal{S}(\epsilon ), \end{aligned}$$

(5.116)

where $\varepsilon =a_{3}/a$. This then gives

$$\begin{aligned} I_{11} &= I_{22} = I_{12} = I_{21} = \pi- \frac{I_{1}-I_{3}}{4(\varepsilon ^{2}-1)}, \end{aligned}$$

(5.117)

$$\begin{aligned} I_{13} &= I_{23} = \frac{I_{1}-I_{3}}{\varepsilon ^{2}-1}, \end{aligned}$$

(5.118)

$$\begin{aligned} I_{33} &= \frac{4\pi}{3} - \frac{2}{3} \varepsilon ^{2} I_{13}, \end{aligned}$$

(5.119)

$$\begin{aligned} I_{31} &= I_{32} = \varepsilon ^{2} I_{13}. \end{aligned}$$

(5.120)

These can then be used in (5.79)–(5.82) together with the cyclic properties to derive the components of the $\mathbf{P}$-tensor for a spheroid. It is straightforward to check that this gives rise to the coefficients $p_{1}$ through $p_{6}$ as defined for a transversely isotropic tensor in (5.51)–(5.53).

5.2 Anisotropic Host Phase

In the potential problem case, scaling coordinate systems assisted in the derivation of results associated with anisotropic media. Although such methods can sometimes lead to modest simplifications in elasticity, the general theory does not lead to any significant advances, certainly for the problems that are of greatest interest in micromechanics. In particular such methods do not lead to significant simplifications for generally transversely isotropic media which is a material symmetry of great importance. Therefore to derive the $\mathbf{P}$-tensors associated with inhomogeneities in anisotropic host phases, it is best to work with the integral form of the $\mathbf{P}$-tensor as defined in (3.9) for an ellipsoid.

Few explicit results are available in general however since the Green’s tensor cannot generally be determined analytically. One of the few that can be determined is that associated with transversely isotropic media. Withers derived the associated Eshelby tensor for a spheroid [133] using the form of the Green’s function determined by Pan and Chou [101] and we state his result here where the semi-major or minor axis of the spheroid is aligned with the axis of transverse isotropy. This result shall then be validated by employing the general integral form (3.9). Only in the last decade have articles started to appear that compute effective properties of composite media with anisotropic inhomogeneities via micromechanical methods (e.g., [42, 116]). It is also important to note specific results for the Eshelby and Hill tensors associated with cracks in anisotropic media, e.g., Gruescu et al. [44] and Barthélémy [9].

5.2.1 Spheroid in a Transversely Isotropic Host Phase

Consider a spheroid with semi-axes $a=a_{1}=a_{2}\neq a_{3}$ embedded in a transversely isotropic host phase where the $x_{1}x_{2}$ plane is the plane of isotropy. The elastic modulus tensor of the host is

$$\begin{aligned} C_{ijk\ell}^{0} &= \sum_{n=1}^{6} c^{0}_{n} \mathcal{H}_{ijk\ell}^{n}, \end{aligned}$$

(5.121)

where

$$\begin{aligned} c_{1}^{0} &= 2K_{0}, \qquad c^{0}_{2} = \ell_{0}, \qquad c^{0}_{3} = \ell_{0}, \end{aligned}$$

(5.122)

$$\begin{aligned} c_{4}^{0} &= n_{0}, \qquad\ \ \; c^{0}_{5} = 2 m_{0}, \quad\; c^{0}_{6} = 2 g_{0}. \end{aligned}$$

(5.123)

Here $K_{0}$ and $m_{0}$ are the in-plane bulk and shear moduli and $g_{0}$ is the antiplane shear modulus (often $p_{0}$ is used for the anti-plane modulus but this is not employed here in order to avoid confusion associated with components $p_{j}$ of the Hill tensor).

Derivation from Withers’ Eshelby Tensor

Withers derived the Eshelby tensor for a spheroid in a transversely isotropic host medium. In order to state this result it is useful to first define the parameters

$$\begin{aligned} v_{1} &= \biggl(\displaystyle\frac{(\hat{\ell}_{0}-\ell_{0})(\hat{\ell }_{0}+\ell _{0}+2g_{0})}{4n_{0}g_{0}} \biggr)^{1/2}+ \biggl(\displaystyle\frac{(\hat {\ell }_{0}+\ell_{0})(\hat{\ell}_{0}-\ell_{0}-2g_{0})}{4n_{0}g_{0}} \biggr)^{1/2},\\ v_{2} &= \biggl(\displaystyle\frac{(\hat{\ell}_{0}-\ell_{0})(\hat{\ell }_{0}+\ell _{0}+2g_{0})}{4n_{0}g_{0}} \biggr)^{1/2}- \biggl(\displaystyle\frac{(\hat {\ell }_{0}+\ell_{0})(\hat{\ell}_{0}-\ell_{0}-2g_{0})}{4n_{0}g_{0}} \biggr)^{1/2},\\ v_{3} &= \biggl(\displaystyle\frac{m_{0}}{g_{0}} \biggr)^{1/2}, \end{aligned}$$

where $\hat{\ell}_{0}=(n_{0}(K_{0}+m_{0}))^{1/2}$. Note that for elastic materials $v_{3}\in\mathbb{R}$, the set of real numbers, but $v_{1},v_{2}\in\mathbb{C}$, the set of complex numbers, in general, with $v_{2}=\overline{v_{1}}$, where an overbar denotes the complex conjugate here (not to be confused with the body average!)¹ For $v_{i}\in\mathbb{R}$ define

$$\begin{aligned} v_{i}I_{3}(v_{i}) &= 4\pi\mathcal{S}(v_{i} \varepsilon ), \qquad I_{1}(v_{1}) = \frac{4\pi}{v_{i}}-2 I_{3}(v_{i}). \end{aligned}$$

(5.124)

When $v_{i}\in\mathbb{C}$, either of the cases in $\mathcal {S}(v_{i}\varepsilon )$ are valid since they are merely an analytic continuation of the function (of $v_{i}$) into the complex $v_{i}$-plane. Note that in the case of isotropy, $\ell_{0}=\lambda _{0}=\ell_{0}'$, $g_{0}=m_{0}=\mu_{0}$, $\hat{\ell}_{0}=n_{0}=K_{0}+m_{0}=\lambda _{0}+2\mu_{0}$ and thus $v_{1}=v_{2}=v_{3}=1$. The notation $I_{i}=I_{i}(1)$ is therefore appropriate for the isotropic case, as already introduced.

For a transversely isotropic host phase defined by elastic properties (5.121) Withers determined the components $S_{ijk\ell}$ in the form

$$\begin{aligned} S_{1111} &= \sum_{i=1}^{2} \bigl[2g_{0}(1+M_{i})v_{i}^{2}-m_{0} \bigr] L_{i}v_{i} I_{1}(v_{i}) + \frac{1}{2}D m_{0} I_{1}(v_{3}), \end{aligned}$$

(5.125)

$$\begin{aligned} S_{1122} &= \sum_{i=1}^{2} \bigl[2g_{0}(1+M_{i})v_{i}^{2} - 3m_{0} \bigr]L_{i} v_{i}I_{1}(v_{i})- \frac{1}{2}D m_{0} I_{1}(v_{3}), \end{aligned}$$

(5.126)

$$\begin{aligned} S_{3333} &= 2\sum_{i=1}^{2} \bigl[ \ell_{0}-n_{0} M_{i}v_{i}^{2} \bigr] v_{i}^{3} M_{i} L_{i} I_{3}(v_{i}), \end{aligned}$$

(5.127)

$$\begin{aligned} S_{1133} &= 2\sum_{i=1}^{2} \bigl[ \ell_{0}-n_{0} M_{i} v_{i}^{2} \bigr]v_{i} L_{i} I_{1}(v_{i}), \end{aligned}$$

(5.128)

$$\begin{aligned} S_{3311} &= 2\sum_{i=1}^{2} \bigl[g_{0} v_{i}^{2}(1+M_{i})-m_{0} \bigr] M_{i}L_{i}v_{i}^{3} I_{3}(v_{i}), \end{aligned}$$

(5.129)

$$\begin{aligned} S_{1313} &= \frac{1}{2}g_{0}\sum _{i=1}^{2} L_{i}v_{i}^{3}(1+M_{i}) \bigl(I_{3}(v_{i})-2M_{i} I_{1}(v_{i}) \bigr)+\frac{1}{4}D g_{0} I_{3}(v_{3})v_{3}^{2}, \end{aligned}$$

(5.130)

where

$$\begin{aligned} D &= \frac{1}{4\pi g_{0} v_{3}}, \qquad M_{i} = \frac {(K_{0}+m_{0})/v_{i}^{2}-g_{0}}{\ell_{0}+g_{0}}, \qquad L_{i} = (-1)^{i}\frac{g_{0}-n_{0}v_{i}^{2}}{8\pi n_{0}g_{0}(v_{1}^{2}-v_{2}^{2})v_{i}^{2}}. \end{aligned}$$

Note that slightly different notation has been used here from that in [133] and in particular the notation $I_{3}(v_{i})$ has been used whereas [133] used $I_{2}(v_{i})$ for this term in the corresponding equations. This is done here to preserve the symmetry with the isotropic case above so that as $v_{i}\rightarrow1$, $I_{3}(v_{i})\rightarrow I_{3}$.

Via straightforward contraction with the transversely isotropic compliance tensor $\mathbf{D}_{0}$, i.e., the inverse of (5.121), (1.2) then yields

$$\begin{aligned} P_{1111}+P_{1122} &= \frac{n_{0}(S_{1111}+S_{1122})-2\ell_{0} S_{1133}}{2\Delta}, \end{aligned}$$

(5.131)

$$\begin{aligned} P_{1111}-P_{1122} &= \frac{S_{1111}-S_{1122}}{2m_{0}}, \end{aligned}$$

(5.132)

$$\begin{aligned} P_{3333} &= \frac{K_{0} S_{3333}-\ell_{0} S_{3311}}{\Delta}, \end{aligned}$$

(5.133)

$$\begin{aligned} P_{1133} &= \frac{n_{0} S_{3311}-\ell_{0}S_{3333}}{2\Delta}, \end{aligned}$$

(5.134)

$$\begin{aligned} P_{1313} &= \frac{S_{1313}}{2g_{0}}, \end{aligned}$$

(5.135)

$$\begin{aligned} P_{1212} &= \frac{S_{1212}}{2m_{0}}, \end{aligned}$$

(5.136)

where $\Delta=K_{0} n_{0}-\ell_{0}^{2}$. Of course this calculation could also be done with the help of matrices rather than tensor forms. Since the $\mathbf{P}$-tensor is transversely isotropic however it is rather straightforward to write down the transversely isotropic tensor basis forms

$$\begin{aligned} P_{ijk\ell} &= \sum_{n=1}^{6} p_{n} \mathcal{H}^{n}_{ijk\ell}, \qquad S_{ijk\ell} = \sum_{n=1}^{6} s_{n} \mathcal{H}^{n}_{ijk\ell}, \end{aligned}$$

(5.137)

where

$$\begin{aligned} p_{1} &= 2P_{1111}-P_{1212}, \qquad p_{2} = p_{3} = P_{1133}, \end{aligned}$$

(5.138)

$$\begin{aligned} p_{4} &= P_{3333}, \qquad\qquad\quad\ \ \ p_{5} = P_{1212}, \qquad \qquad p_{6} = 2P_{1313} \end{aligned}$$

(5.139)

and similarly for the Eshelby tensor with $p_{n}\rightarrow s_{n}$ and $P_{ijk\ell}\rightarrow S_{ijk\ell}$ of course.

Derivation from the Direct Integral Form

As noted above, the $\mathbf{P}$-tensor will itself be transversely isotropic of the form given in (5.137). Using the direct integral formulation of the $\mathbf{P}$-tensor (3.9), let the unit vector $\overline{\boldsymbol{\xi}}$ pointing to the surface of the unit sphere be parametrized by the two angles $\varphi\in[0,2\pi)$ and $\vartheta\in[0,\pi)$, so that it takes the form

$$\begin{aligned} \overline{\xi}_{1} &= \cos\varphi\sin\vartheta, \qquad\overline{ \xi}_{2} = \sin\varphi\sin\vartheta, \qquad\overline{\xi}_{3} = \cos\vartheta. \end{aligned}$$

(5.140)

Therefore, together with (5.121), (3.9) becomes

$$\begin{aligned} P_{ijk\ell} &= \frac{\varepsilon }{4\pi} \int_{0}^{2\pi} \int_{0}^{\pi}\frac{\varPhi _{ijkl}}{(1+(\varepsilon ^{2}-1)\cos^{2}\vartheta)^{3/2}}\sin\vartheta \,d \varphi\, d\vartheta, \end{aligned}$$

(5.141)

where

$$\begin{aligned} \varPhi_{ijk\ell} &= (\overline{\xi}_{j}\overline{ \xi}_{\ell }N_{ik} ) |_{(ij),(k\ell)} \end{aligned}$$

(5.142)

with $N_{ij}$ defined via $N_{ik}\tilde{N}_{kj} = \delta_{ij}$. The components of $\tilde{\mathbf{N}}$ are defined as $\tilde {N}_{ij}=C_{ijk\ell}^{0}\overline{\xi}_{j}\overline{\xi}_{\ell}$. It is straightforward to implement this in a variety of mathematical packages or programming languages.

Here let us plot the five independent components of the $\mathbf{P}$-tensor, illustrating that the two approaches (Withers’ explicit form and the direct integral evaluation) described above agree. Let us take the elastic properties (associated with the host phase) to be transversely isotropic and choose the material PZT-7A.² Although realistically this material would normally be chosen as the reinforcing phase in a composite, it is appropriate to illustrate the calculations for a model material. Its properties are (all stated in GPa)

$$\begin{aligned} K &= 121.2, \qquad m = 35.8, \qquad\ell= 73, \qquad n= 175, \qquad g = 47.2. \end{aligned}$$

(5.143)

In Figs. 2–3 the components of the $\mathbf{P}$-tensor are plotted, using the explicit form arising from the Eshelby tensor and using the direct evaluation of the integral in order to confirm the results.

Fig. 2

Plot of the components $P_{1111}$ (solid black), $P_{3333}$ (dashed blue) and $P_{1313}$ (dotted red) of the Hill tensor associated with a spheroidal inhomogeneity (varying the aspect ratio $\varepsilon$) embedded in the transversely isotropic host phase PZT-7A, using the form of the $\mathbf{P}$-tensor derived from the explicit form of the Eshelby tensor for this problem. These explicit results are confirmed by calculating the components at discrete values of the aspect ratio by evaluating the direct integral form of the $\mathbf{P}$-tensor given in (5.141). Note that the limiting values as $\varepsilon \rightarrow \infty$ correspond to the circular cylinder result (5.29) and when $\varepsilon \rightarrow0$ the layer limit is obtained (Color figure online)

Fig. 3

As with Fig. 2 but for the components $P_{1122}$ (solid black) and $P_{1133}$ (dashed red) of the $\mathbf{P}$-tensor (Color figure online)

Concentration Tensor

Since the $\mathbf{P}$-tensor is known one can now go on to deduce the associated Concentration tensor. First assume that the spheroid is isotropic with elastic modulus tensor

$$\begin{aligned} C^{1}_{ijkl} &= 3\kappa _{1} I_{ijkl}^{1} + 2\mu_{1} I_{ijkl}^{2}. \end{aligned}$$

(5.144)

Since the Concentration tensor will be transversely isotropic, it is convenient to write $C^{1}_{ijkl}$ with respect to the transversely isotropic tensor basis so that

$$\begin{aligned} C^{1}_{ijkl} &= \sum_{n=1}^{6} c^{1}_{n} \mathcal{H}_{ijkl}^{n}, \end{aligned}$$

(5.145)

where the $c^{1}_{n}$ coefficients are defined in terms of the two independent elastic moduli $\kappa _{1}$ and $\mu_{1}$:

$$\begin{aligned} c^{1}_{1} &= 2\kappa _{1}+\frac{2}{3} \mu_{1}, \qquad c^{1}_{2} = \kappa _{1}- \frac{2}{3}\mu_{1}, \qquad c^{1}_{3} = \kappa _{1}-\frac{2}{3}\mu_{1}, \end{aligned}$$

(5.146)

$$\begin{aligned} c^{1}_{4} &= \kappa _{1}+\frac{4}{3} \mu_{1}, \qquad\ \ c^{1}_{5} = 2\mu_{1}, \qquad \qquad \,\, c^{1}_{6} = 2\mu_{1}. \end{aligned}$$

(5.147)

Let us employ (5.121), together with the form of $\mathbf{P}$-tensor defined in (5.137). The properties of the transversely isotropic basis tensors $\mathcal{H}_{ijkl}^{n}$ described in Appendix C.2.3 shall also be exploited (and in particular the contraction properties in Table 2), together with the expressions written down in (C.26)–(C.28). The inverse of the Concentration tensor defined in (3.23) can then be determined in the form

$$\begin{aligned} \tilde{\mathcal{A}}_{ijkl} &= \sum_{n=1}^{6} \tilde{\alpha}_{n} \mathcal{H}_{ijkl}^{n}, \end{aligned}$$

(5.148)

where upon defining $c_{n}=c^{1}_{n}-c^{0}_{n}$,

$$\begin{aligned} \tilde{\alpha}_{1} &= 1+p_{1} c_{1}+2 p_{2} c_{3}, \qquad \tilde{\alpha}_{2} = p_{1} c_{2}+p_{2} c_{4}, \qquad \tilde{ \alpha}_{3} = p_{3} c_{1}+p_{4} c_{3}, \end{aligned}$$

(5.149)

$$\begin{aligned} \tilde{\alpha}_{4} &= 1+2 p_{3} c_{2} + p_{4} c_{4}, \qquad \tilde{\alpha}_{5} = 1+p_{5} c_{5}, \qquad \quad \ \, \tilde{\alpha}_{6} = 1+p_{6}c_{6}. \end{aligned}$$

(5.150)

The tensor $\tilde{\mathcal{A}}_{ijkl}$ is then inverted, following the procedure in Appendix C.2.3, to yield the Concentration tensor

$$\begin{aligned} \mathcal{A}_{ijkl} &= \sum_{n=1}^{6} \alpha_{n} \mathcal{H}_{ijkl}^{n}, \end{aligned}$$

(5.151)

where

$$\begin{aligned} \alpha_{1} &= \frac{\tilde{\alpha}_{1}}{2\Delta}, \qquad\alpha_{2} = - \frac{\tilde{\alpha}_{2}}{2\Delta}, \qquad\alpha_{3} = -\frac {\tilde{\alpha}_{3}}{2\Delta}, \end{aligned}$$

(5.152)

$$\begin{aligned} \alpha_{4} &= \frac{\tilde{\alpha}_{4}}{2\Delta}, \qquad\alpha_{5} = \frac {1}{\tilde{\alpha_{5}}}, \qquad\quad\, \alpha_{6} = \frac{1}{\tilde {\alpha}_{6}}. \end{aligned}$$

(5.153)

Alternatively suppose that the inhomogeneity is transversely isotropic with the same symmetry axis as the host so that it possesses the elastic modulus tensor of the form (5.145) but where now the constants $c^{1}_{n}$ are defined in terms of the 5 independent components of this tensor. Then the Concentration tensor is again defined by (5.151) but of course now with the $c^{1}_{n}$ associated with the transversely isotropic cylinder. This indicates the merit of using the above notation since one can still use (5.151)–(5.153) in this case, merely modifying the $c^{1}_{n}$ to account for the transverse isotropy of the cylinder.

As usual, the matrix form of these fourth order tensors can be employed for computational efficiency when the problems lack simple symmetries.

5.2.2 Circular Cylinder in a Transversely Isotropic Host

Suppose now that the inhomogeneity is a circular cylinder with its cross-section residing in the plane of isotropy of the transversely isotropic host phase. One can arrive at the corresponding $\mathbf{P}$-tensor in two ways. The first is to take the limit $\varepsilon =a_{3}/a\rightarrow\infty$ in the prolate spheroid case in Sect. 5.2.1. The second way is to recognize that since the anisotropy of the host will not affect the in-plane components of the $\mathbf{P}$-tensor, the tensor will simply be the same as that for an isotropic host as derived in Sect. 5.1.2 but the elastic properties are modified via $\lambda _{0}+\mu_{0}\rightarrow K_{0}$ and $\mu_{0}\rightarrow m_{0}$ for in-plane components and $\mu _{0}\rightarrow g_{0}$ for the anti-plane component. Therefore, from (5.29)

$$\begin{aligned} p_{1} &= \frac{1}{2(K_{0}+m_{0})}, \qquad p_{2} = 0, \qquad\qquad \qquad\quad \ \ \; p_{3}=0, \end{aligned}$$

(5.154)

$$\begin{aligned} p_{4} &= 0, \qquad\qquad \qquad \ \ p_{5} = \frac{K_{0}+2m_{0}}{4m_{0}(K_{0}+m_{0})}, \qquad p_{6} = \frac {1}{4g_{0}}. \end{aligned}$$

(5.155)

The Concentration tensor may then be derived by using these coefficients in (5.149)–(5.150) and the expressions that follow.

6 Discussion

6.1 Association with Micromechanics

One of the primary reasons for deriving the Hill or Eshelby tensors and associated Concentration tensors is to understand a multitude of aspects of the behaviour of inhomogeneous media, including their macroscopic constitutive response and so-called effective properties. Following a relatively straightforward argument regarding volume averaging, the effective modulus tensor $\mathbf{C}^{*}$ of an $n+1$ phase medium with a distinguishable host phase (phase 0) can be stated as [81, 132]

$$\begin{aligned} \mathbf{C}^{*} &= \mathbf{C}^{0} + \sum_{r=1}^{n} \phi_{r} \bigl(\mathbf{C}^{r}-\mathbf{C}^{0} \bigr)\mathbf{A}^{r}, \end{aligned}$$

(6.1)

where $\phi_{r}$ is the volume fraction of phase $r$ and $\mathbf{A}^{r}$ is the exact Concentration tensor associated with embedded phase $r$. This is in contrast to the Concentration tensor $\boldsymbol {\mathcal{A}}$ introduced in Sect. 3.3 which is the Concentration tensor associated with an isolated inhomogeneity, i.e., the presence of other inhomogeneities is not accounted for in $\boldsymbol{\mathcal{A}}$. Consequently if the inhomogeneity phases are distributed dilutely then one can merely use the approximation $\mathbf{A}\approx\boldsymbol{\mathcal {A}}$ in (6.1). This links directly therefore to the moment tensor defined for (isolated) ellipsoidal inhomogeneities in (1.3), which permits (6.1) to be written as

$$\begin{aligned} \mathbf{C}^{*} &= \mathbf{C}^{0} + n_{r}\sum _{r=1}^{n} \mathbf{M}^{r}, \end{aligned}$$

(6.2)

where $n_{r}=N_{r}/|V|$ is the so-called number density of the inhomogeneities of phase $r$, with $N_{r}$ being the number of those inhomogeneities of that phase. We stress that (6.2) is valid here for dilute volume fractions only. For more discussion on the use of the moment tensor in micromechanics see [2, 3].

Many micromechanical methods use a more sophisticated approximation that can account, in an approximate manner at least, for interactions. One of the most commonly employed methods is the so-called classical self consistent method [54]. Interaction is approximated in this most simple self-consistent scheme by taking the host medium in the determination of $\boldsymbol{\mathcal{A}}$ to be the unknown effective medium. In general then (6.1) gives rise to a nonlinear system of equations for the determination of effective properties. In many cases these are not even algebraic equations. Furthermore for the self consistent method, one has to make an assumption in advance of the symmetry properties of the effective tensor. For example in the case of aligned spheroids $\mathbf{C}^{*}$ will be transversely isotropic.

The textbooks referred to at the end of Sect. 1 provide an excellent introduction to the numerous micromechanical methods, many of which are based on the form of effective modulus tensor defined in (6.1). A similar form can be deduced for media where there is no distinguishable host phase (for example the case of polycrystals) and also for media where multiphysics effects are important as described in the next section.

6.2 Beyond the Potential Problem and Elastostatics

A large number of explicit, compact results associated with the Hill and Eshelby tensors for ellipsoidal inhomogeneities, as well as their associated Concentration tensors have been collected, stated and in some cases derived. The intention is that this will be of great utility to a large number of researchers for implementation in micromechanical and bounding schemes. A thorough discussion of both matrix and tensor (where possible due to space limitations) formulations has been carried out. Typographical errors in past articles and reviews have been corrected and a common notation has been employed.

Although the general integral forms (3.3) and (3.9) are useful they should generally be avoided where explicit forms are available. Gavazzi and Lagoudas [40] described a numerical implementation for elasticity. It should be noted that recently Masson [82] derived a new form of the $\mathbf{P}$-tensor in terms of a single integral, although the integrand is inevitably more complex than that in the surface integral in (3.9).

In the literature many of the cases described above are considered as approximations to more complicated shaped inhomogeneities. In terms of the derivation of overall effective properties this is extremely useful, certainly as a first approximation, since it avoids complex computational simulations. However it must be stressed that more advanced analysis is required if detailed micromechanical information such as stress concentration calculations close to inhomogeneities of a complex shape is required [15]. For finite domains, provided the host phase is in some sense much larger than the inhomogeneity, if the inhomogeneity is ellipsoidal then the temperature gradient or strain field inside the ellipsoid is well approximated as being uniform. The inhomogeneity problem associated with bounded domains is described in the book by Li and Wang [67] which summarizes the work in [68, 69].

Continuing with discussion associated with the potential problem and elastostatics first, an important extension of the inhomogeneity problem is that of the coated inhomogeneity. This problem is popular, not least because it arises as a micromechanics problem in the generalized self-consistent method (GSCM) [23]. The so-called double inclusion problem dates back many decades and was solved approximately by Hori and Nemat-Nasser [50] although the approximations involved lead to some rather counter-intuitive predictions when used in the GSCM [51]. Exact solutions in the case of concentric spheroids or ellipsoids have been derived by Hatta and Taya [46] in the thermal context and Jiang et al. [52] in two-dimensional elasticity. See also [75, 76] for elasticity problems involving coated spheres and cylinders respectively. The case of inhomogeneities with radially dependent material properties has been considered in Chap. 3 of [54] amongst others. The coated inhomogeneity is also of interest due to its association with the neutral inclusion problem [89]. Associated with the coated inhomogeneity is the scenario when the interface of an inhomogeneity with the host phase is imperfect [37, 63]. This imperfection can itself be used as the basis for a neutral inclusion [12, 109]. Moment tensors have recently been discussed in the context of cloaking [5].

It is important to note that when the inhomogeneity becomes very small, i.e., the case of a nano-inhomogeneity, then surface energies become non-negligible. This problem has been considered by Sharmi and Ganti [117] and Duan et al. [25] for example. Including surface energies is important in order to incorporate size-dependent effects in effective properties. These are absent in classical micromechanical methods that use the standard Eshelby or Hill tensors.

Eshelby’s problem has also been considered in the context of micro-continuum elasticity models, which themselves were introduced in order to bridge the gap between continuum and atomistic/molecular models [30]. Micropolar (Cosserat) theory has been considered by Cheng and He [20, 21] and Ma and Hu [77]. Micro-stretch theory has been developed by Ma and Hu [78]. Strain gradient constitutive behaviour was studied by Gao and Ma [38, 39].

The dynamic problem was considered for spheres and cylinders by Mikata and Nemat-Nasser [87, 88] and more generally in [19, 84]. Rate dependence of the Hill and Eshelby tensors has been considered by Suvarov and Dvorak [120] and viscoelastic properties have been studied by Wang and Weng [125] by using transform techniques and correspondence principles. Nguyen et al. [96] studied cracked viscoelastic solids using the appropriate Eshelby tensor. Extensions to plasticity have been considered, see [36, 53, 64] for example.

The Newtonian potential and elastostatics problems are canonical problems that can assist with the development of coupled (multiphysics) problems. Dunn and Taya [26], Dunn and Wienecke [27] and Mikata [85, 86] considered the case of piezoelectricity and the prediction of the electroelastic moduli. Li and Dunn [66] and Zhang and Soh [136] considered full coupling and the resulting effective moduli associated with piezoelectromagnetic media. The theory associated with poroelastic and thermoelastic behaviour was developed by Berryman [11] and extended to the anisotropic case by Levin and Alvarez-Tostado [65].

Upon closing it should be noted that it is very fortuitous that such elegant and concise uniformity results hold for ellipsoidal inhomogeneities. These results allow a large number of expressions to be derived analytically and therefore the results have been utilized a great deal. Having said that there is much work to be done. As has been noted, the Eshelby conjecture is still not fully resolved [4], analysis for general shaped inhomogenities continues [16], specifically in the context of stress analysis and resulting effective properties and although computational methods are powerful, they are still only able to solve elasticity problems for inhomogeneous media with an order of 1000 inhomogeneities in “reasonable” times. For use in Monte-Carlo schemes this is therefore still computationally expensive. Nonlinear problems in the context of finite elasticity still require attention [135] and this applies to coupled problems as well.

Appendices

Appendix A: Uniform $\mathbf{P}$-Tensors for Ellipsoidal Inhomogeneities

Fourier transforms can be applied in a straightforward manner to derive forms of the Green’s tensors that are useful in the context of deriving properties of the Hill and Eshelby tensors. For arbitrary anisotropy in the potential problem the Green’s function takes the form [92]

$$\begin{aligned} G(\mathbf{z}) &= \frac{1}{16\pi^{3}} \int_{S^{2}} \int_{-\infty}^{\infty}\frac {1}{C_{ij}\bar{\xi}_{i}\bar{\xi}_{j}}\exp(i\xi\bar{ \boldsymbol{\xi}}\cdot\mathbf{z})\,d\xi\,dS(\bar{\boldsymbol{\xi}}), \end{aligned}$$

(A.1)

where the Fourier transform variable $\boldsymbol{\xi}=\xi\bar{\boldsymbol{\xi}}$ with $\xi=|\boldsymbol{\xi}|$ and where $S^{2}$ corresponds to $\xi=1$, the surface of the unit sphere. Next since

$$\begin{aligned} \delta(\bar{\boldsymbol{\xi}}\cdot\mathbf{z}) &= \frac{1}{2\pi} \int_{-\infty}^{\infty}\exp(i\xi\bar{\boldsymbol{\xi}}\cdot \mathbf{z})\, d\xi \end{aligned}$$

(A.2)

the form (A.1) becomes

$$\begin{aligned} G(\mathbf{z}) &= \frac{1}{8\pi^{2}} \int_{S^{2}}\delta(\bar{\boldsymbol{\xi}}\cdot\mathbf{z}) \frac{1}{C_{ij}\bar{\xi}_{i}\bar{\xi}_{j}} \,dS(\bar{\boldsymbol{\xi }}). \end{aligned}$$

(A.3)

In elastostatics an entirely analogous approach shows that

$$\begin{aligned} G_{ij}(\mathbf{z}) &= \frac{1}{8\pi^{2}} \int_{S^{2}}\delta(\bar{\boldsymbol{\xi}}\cdot \mathbf{z})N_{ij}( \bar{\boldsymbol{\xi}}) \,dS(\bar{\boldsymbol{\xi}}), \end{aligned}$$

(A.4)

where

$$\begin{aligned} \tilde{N}_{ik}(\boldsymbol{\xi})N_{kj}(\boldsymbol{\xi}) = \delta_{ij}, \qquad\tilde{N}_{ij}(\boldsymbol{\xi}) = C_{ijk\ell}\xi_{j}\xi_{\ell}. \end{aligned}$$

(A.5)

1.1 A.1 The Potential Problem

Substitute the general form (A.3) of the free-space Green’s function into (3.2) and with $\mathbf{x}\in V_{1}$ this leads to the following form of the Hill tensor:

$$\begin{aligned} P_{ij}(\mathbf{x}) &= -\frac{1}{8\pi^{2}}\frac{\partial ^{2}{}}{\partial{x_{i}}\partial{x_{j}}} \int_{S^{2}}{\frac{1}{C^{0}_{k\ell }\overline{\xi}_{k}\overline{\xi}_{\ell}}} J(\overline{\boldsymbol {\xi}} \cdot\mathbf{x})\,dS, \end{aligned}$$

(A.6)

where

$$\begin{aligned} J(p) &= \int_{\mathbf{y}\in V_{1}}\delta(p-\overline{\boldsymbol{\xi}}\cdot \mathbf{y}) \,d \mathbf{y}. \end{aligned}$$

(A.7)

Take the simplest case, where $V_{1}$ is a sphere of radius $a$. With $\mathbf{x}\in V_{1}$ then $p\leq a$ since $\overline{\boldsymbol{\xi}}$ is a unit vector. It is then recognized that the value of $J(p)$ is the area of the disc defined by the intersection of the plane $\overline{\boldsymbol{\xi}}\cdot\mathbf{y}=p$ with the sphere $V_{1}$. Since $|p|\leq a$, $J(p) = \pi(a^{2}-p^{2})$ and carrying out the necessary differentiation gives

$$\begin{aligned} P_{ij}^{\mathrm {sphere}} = \frac{1}{4\pi} \int_{S^{2}} \frac{\overline{\xi}_{i}\overline{\xi}_{j}}{C_{k\ell }^{0}\overline{\xi }_{k}\overline{\xi}_{\ell}}\,dS. \end{aligned}$$

(A.8)

Since this integral is over the unit sphere $S^{2}$ and involves only $\overline{\xi}_{i}$ then this shows that the $\mathbf{P}$-tensor is independent of $\mathbf{x}$ for a spherical inhomogeneity. This was evident from the fact that $J(p)$ is a quadratic function of its argument. In order to perform the integration over the unit sphere, introduce the parameters $\vartheta\in[0,\pi)$ and $\varphi\in[0,2\pi)$ via

$$\begin{aligned} \overline{\xi}_{1} = \cos\varphi\sin\vartheta, \qquad\overline{ \xi}_{2} = \sin\varphi\sin\vartheta, \qquad\overline{\xi}_{3} = \cos\vartheta \end{aligned}$$

(A.9)

and then

$$\begin{aligned} P_{ij}^{\mathrm {sphere}} &= \frac{1}{4\pi} \int_{0}^{2\pi} \int_{0}^{\pi} \varPhi_{ij}(\overline{ \boldsymbol{\xi}})\sin\vartheta d\vartheta d\varphi, \end{aligned}$$

(A.10)

where

$$\begin{aligned} \varPhi_{ij}(\overline{\boldsymbol{\xi}}) &= \frac{\overline{\xi }_{i}\overline{\xi}_{j}}{C_{k\ell}\overline{\xi }_{k}\overline{\xi}_{\ell}}. \end{aligned}$$

(A.11)

It is straightforward to extend this derivation in order to derive the corresponding result for an ellipsoid. Suppose that the ellipsoid is defined by $V=\{\mathbf{y}:\mathbf{y}^{T} (\mathbf{a}^{T}\mathbf{a})^{-1}\mathbf{y}<1\}$ where $\mathbf{a}$ is the second order tensor defined in (3.5) with $a_{j}$ being the semi-axes of the ellipsoid. With $\mathbf{x}\in V_{1}$, the only aspect that changes from the calculation for the sphere is that now the function $J(p)$ will be the area of the region defined by the intersection of the plane $\overline{\boldsymbol{\xi}}\cdot\mathbf{y}=p$ with the ellipsoid $V_{1}$. It transpires that

$$\begin{aligned} J(p) &= \frac{\det({\mathbf{a}})}{(\overline{\xi}_{k} a_{k\ell }a_{\ell m}\overline{\xi}_{m})^{3/2}}\pi \bigl(a^{2}-p^{2} \bigr). \end{aligned}$$

(A.12)

Importantly the integral is still only over the unit sphere and the subsequent result is (3.3).

1.2 A.2 Elastostatics

One can proceed entirely analogously to the potential problem case in order to derive the representation of the $\mathbf{P}$-tensor for an ellipsoid in the elastostatics context. Using (A.4) and following the same procedure as for the potential problem one obtains (3.9).

Appendix B: Potential Theory

Two integrals over ellipsoidal domains arise in potential theory, having important applications in micromechanics. Define the two functions

$$\begin{aligned} \varGamma(\mathbf{x}) = -\frac{1}{4\pi} \int_{V_{1}}\frac{1}{|\mathbf {x}-\mathbf{y}|}\,d\mathbf{y}, \qquad\varPsi( \mathbf{x}) = \frac{1}{4\pi} \int_{V_{1}}|\mathbf{x}-\mathbf{y}|\, d\mathbf{y}, \end{aligned}$$

(B.1)

where $V_{1}$ is the ellipsoidal domain defined by the inequality

$$\begin{aligned} \frac{y_{1}^{2}}{a_{1}^{2}} + \frac{y_{2}^{2}}{a_{2}^{2}} + \frac{y_{3}^{2}}{a_{3}^{2}} &\leq1 \end{aligned}$$

(B.2)

and of interest in (B.1) is the case when $\mathbf{x}\in V_{1}$.

Introduce local spherical polar coordinates via $y_{j}=x_{j}+z_{j}, j=1,2,3$ where

$$\begin{aligned} z_{1} = r\cos\varphi\sin\vartheta, \qquad z_{2} = r\sin \varphi\sin\vartheta, \qquad z_{3} = r\cos\vartheta, \end{aligned}$$

(B.3)

where $\vartheta\in[0,\pi), \varphi\in[0,2\pi)$ and $r\in [0,\infty)$. The surface of the ellipsoid is given by $r=R_{1}(\vartheta,\varphi)$ and the differential volume element is $dV=r^{2}\sin\vartheta drd\vartheta d\varphi$ so that

$$\begin{aligned} \varGamma(\mathbf{x}) &= -\frac{1}{8\pi} \int_{0}^{2\pi} \int_{0}^{\pi}R_{1}^{2}(\vartheta, \varphi)\sin\vartheta\,dr\,d\vartheta\, d\varphi \end{aligned}$$

(B.4)

and similarly

$$\begin{aligned} \varPsi(\mathbf{x}) &= \frac{1}{16\pi} \int_{0}^{2\pi} \int_{0}^{\pi }R_{1}^{4}(\vartheta, \varphi)\sin\vartheta d\vartheta d\varphi. \end{aligned}$$

(B.5)

Introduce the shifted variable $\psi=\vartheta-\pi/2$ and upon defining $r_{1}(\psi,\varphi)=R_{1}(\psi+\pi/2,\varphi)$ the integrals become

$$\begin{aligned} \varGamma(\mathbf{x}) &= -\frac{1}{8\pi} \int_{-\pi/2}^{\pi/2} \int_{0}^{2\pi }r_{1}^{2}(\psi, \varphi)\cos\psi\, d\varphi\,d\psi, \end{aligned}$$

(B.6)

$$\begin{aligned} \varPsi(\mathbf{x}) &= \frac{1}{16\pi} \int_{-\pi/2}^{\pi/2} \int_{0}^{2\pi }r_{1}^{4}(\psi, \varphi)\cos\psi\,d\varphi\,d\psi. \end{aligned}$$

(B.7)

At this point note that

$$\begin{aligned} \int_{-\pi/2}^{\pi/2} \int_{0}^{2\pi} f(\psi,\varphi)\,d\theta\, d\psi&= \int_{0}^{\pi/2} \int_{0}^{\pi} \bigl[f(\psi,\varphi)+f(\psi,\varphi+ \pi) \bigr] \\ &\quad{}+ \bigl[ f(-\psi,\varphi)+f(-\psi,\varphi+\pi) \bigr]\, d\varphi\,d\psi. \end{aligned}$$

(B.8)

This pairing is useful to argue that certain integrals below are zero. Evaluating the local spherical polar coordinates on the surface of the ellipsoid yields

$$\begin{aligned} A r_{1}^{2} + 2B r_{1} + C &= 0, \end{aligned}$$

(B.9)

where

$$\begin{aligned} A &= \frac{\cos^{2}\psi\cos^{2}\varphi}{a_{1}^{2}} + \frac{\cos^{2}\psi \sin ^{2}\varphi}{a_{2}^{2}} + \frac{\sin^{2}\psi}{a_{3}^{2}}, \end{aligned}$$

(B.10)

$$\begin{aligned} B &= \frac{x_{1}\cos\psi\cos\varphi}{a_{1}^{2}} + \frac{x_{2}\cos\psi \sin\varphi }{a_{2}^{2}} - \frac{x_{3}\sin\psi}{a_{3}^{2}}, \end{aligned}$$

(B.11)

$$\begin{aligned} C &= \frac{x_{1}^{2}}{a_{1}^{2}} + \frac{x_{2}^{2}}{a_{2}^{2}} + \frac{x_{3}^{2}}{a_{3}^{2}} - 1, \end{aligned}$$

(B.12)

noting that $A>0$ and $C<0$, and thus

$$\begin{aligned} r_{1} &= \frac{-B+\sqrt{B^{2}-AC}}{A}, \end{aligned}$$

(B.13)

where the positive root is chosen since $B^{2}-AC>B^{2}$ ($AC<0$). Hence

$$\begin{aligned} \varGamma(\mathbf{x}) &= \frac{1}{8\pi} \int_{-\pi/2}^{\pi/2} \int_{0}^{2\pi} \frac{AC-2B^{2}+2B\sqrt{B^{2}-AC}}{A^{2}}\cos\psi\, d\varphi \, d\psi \end{aligned}$$

(B.14)

and

$$\begin{aligned} \varPsi(\mathbf{x}) &= \frac{1}{16\pi} \int_{-\pi/2}^{\pi/2} \int_{0}^{2\pi} \biggl(\frac{8B^{4}}{A^{4}} - \frac{8B^{2} C}{A^{3}} + \frac{C^{2}}{A^{2}} -\frac{8B^{3}\sqrt{B^{2}-AC}}{A^{4}} \\ &\quad{} +\frac{4B C\sqrt{B^{2}-AC}}{A^{3}} \biggr) \cos\psi\, d\varphi\, d\psi. \end{aligned}$$

(B.15)

The radical contributions to both potentials can be shown to be zero by appealing to (B.8) since it transpires that the relevant integrand $f(\psi,\varphi)$ possesses the symmetry

$$\begin{aligned} f(\psi,\varphi) = -f(-\psi,\varphi+\pi), \qquad f(\psi,\varphi +\pi) = -f(-\psi, \varphi). \end{aligned}$$

(B.16)

The functions $\varGamma$ and $\varPsi$ thus reduce to

$$\begin{aligned} \varGamma(\mathbf{x}) &= C \varUpsilon- \varGamma_{1}, \end{aligned}$$

(B.17)

where

$$\begin{aligned} \varGamma_{1} = \frac{1}{4\pi} \int_{-\pi/2}^{\pi/2} \int_{0}^{2\pi}\frac {B^{2}}{A^{2}}\cos\psi\, d\varphi\, d\psi, \qquad\varUpsilon&= \frac{1}{8\pi} \int_{-\pi/2}^{\pi/2} \int_{0}^{2\pi}\frac{\cos \psi}{A}\, d\varphi\, d \psi \end{aligned}$$

(B.18)

and

$$\begin{aligned} \varPsi(\mathbf{x}) &= \varPsi_{1} - C \varPsi_{2} + C^{2} \varOmega, \end{aligned}$$

(B.19)

where

$$\begin{aligned} &\varPsi_{1} = \frac{1}{2\pi} \int_{-\pi/2}^{\pi/2} \int_{0}^{2\pi}\frac {B^{4}}{A^{4}}\cos\psi\,d\varphi\, d\psi, \qquad\varPsi_{2} = \frac{1}{2\pi} \int_{-\pi/2}^{\pi/2} \int_{0}^{2\pi}\frac {B^{2}}{A^{3}}\cos\psi\,d\varphi\,d \psi, \end{aligned}$$

(B.20)

$$\begin{aligned} &\varOmega= \frac{1}{16\pi} \int_{-\pi/2}^{\pi/2} \int_{0}^{2\pi}\frac {1}{A^{2}}\cos\psi\, d\varphi\, d\psi. \end{aligned}$$

(B.21)

2.1 B.1 Closed Integral Form for $\varGamma(\mathbf{x})$

Write $\varGamma_{1} = \varGamma_{11}+\varGamma_{12}$, where

$$\begin{aligned} \varGamma_{11} &= \frac{1}{4\pi} \int_{-\pi/2}^{\pi/2} \int_{0}^{2\pi} \biggl(\frac{\cos^{2}\psi\cos^{2}\varphi}{a_{1}^{2}} \frac{x_{1}^{2}}{a_{1}^{2}}+ \frac{\cos^{2}\psi\sin^{2}\varphi}{a_{2}^{2}}\frac {x_{2}^{2}}{a_{2}^{2}}+\frac{\sin ^{2}\psi}{a_{3}^{2}} \frac{x_{3}^{2}}{a_{3}^{2}} \biggr) \frac{\cos\psi}{A^{2}}\,d\varphi\,d\psi \end{aligned}$$

(B.22)

and

$$\begin{aligned} \varGamma_{12} &= \frac{1}{2\pi} \int_{-\pi/2}^{\pi/2} \int_{0}^{2\pi} \biggl(\frac{x_{1}x_{2}\cos^{2}\psi\sin\varphi\cos \varphi}{a_{1}^{2}a_{2}^{2}} - \frac{x_{2}x_{3}\cos\psi\sin\psi\sin\varphi}{a_{2}^{2}a_{3}^{2}} \\ &\quad{}- \frac{x_{3}x_{1}\cos\psi\sin\psi\cos\varphi}{a_{3}^{2}a_{1}^{2}} \biggr)\frac{\cos\psi}{A^{2}}\,d\varphi\, d\psi. \end{aligned}$$

(B.23)

The contribution from $\varGamma_{12}$ is zero—the first term due to $2\pi$ periodicity of the integrand in $\varphi$ and the second and third terms due to their being odd in $\psi$. Treating $\varUpsilon$ as a function of $a_{j}$, the form of $\varGamma_{1}$ can be exploited, writing

$$\begin{aligned} \varGamma(\mathbf{x}) &= C \varUpsilon- \sum_{j=1}^{3} \frac{x_{j}^{2}}{a_{j}}\frac{\partial{\varUpsilon}}{\partial{a_{j}}}. \end{aligned}$$

(B.24)

Therefore once $\varUpsilon$ is determined, $\varGamma(\mathbf{x})$ straightforwardly follows. Now introduce $A$ into the form of $\varUpsilon$ in (B.18) to obtain

$$\begin{aligned} \varUpsilon&= \frac{1}{8\pi} \int_{-\pi/2}^{\pi/2}\cos\psi \int_{0}^{2\pi} \frac{1}{M(\psi)\cos^{2}\varphi +N(\psi)\sin^{2}\varphi }\,d\varphi\,d\psi, \end{aligned}$$

(B.25)

where

$$\begin{aligned} M(\psi) = \frac{\cos^{2}\psi}{a_{1}^{2}} + \frac{\sin^{2}\psi}{a_{3}^{2}}, \qquad N(\psi) = \frac{\cos^{2}\psi}{a_{2}^{2}} + \frac{\sin^{2}\psi}{a_{3}^{2}}. \end{aligned}$$

(B.26)

Next, the evenness of the integrand is exploited in order to write it as

$$\begin{aligned} \varUpsilon&= \frac{1}{\pi} \int_{0}^{\pi/2}\cos\psi \int_{0}^{\pi/2}\frac{\sec ^{2}\varphi}{M+N\tan^{2}\varphi}\,d\varphi\, d\psi \\ &= \frac{1}{2} a_{1} a_{2} a_{3}^{2} \int_{0}^{\pi/2}\frac{\cos\psi}{\sqrt {(a_{1}^{2}\sin^{2}\psi+a_{3}^{2}\cos^{2}\psi)(a_{2}^{2}\sin^{2}\psi+a_{3}^{2}\cos ^{2}\psi)}}\, d\psi. \end{aligned}$$

(B.27)

Make the substitution $\sin\psi= a_{3} \big/\sqrt{a_{3}^{2}+t}$, where $t\in [0,\infty)$ so that

$$\begin{aligned} \varUpsilon&= \frac{1}{4} a_{1}a_{2}a_{3} \int_{0}^{\infty}\frac{dt}{\Delta(t)}, \end{aligned}$$

(B.28)

where $\Delta(t)=\sqrt{(a_{1}^{2}+t)(a_{2}^{2}+t)(a_{3}^{2}+t)}$. Finally therefore using (B.24)

$$\begin{aligned} \varGamma(\mathbf{x}) &= \frac{1}{4}a_{1}a_{2}a_{3} \int_{0}^{\infty}\frac{X(\mathbf{x},t)\,dt}{\sqrt {(a_{1}^{2}+t)(a_{2}^{2}+t)(a_{3}^{2}+t)}}, \end{aligned}$$

(B.29)

where

$$\begin{aligned} X(\mathbf{x},t) &= \sum_{n=1}^{3} \frac{x_{n}^{2}}{a_{n}^{2}+t}-1 \end{aligned}$$

(B.30)

is a quadratic polynomial in $\mathbf{x}$. It is then found that

$$\begin{aligned} \frac{\partial^{2}{\varGamma}}{\partial{x_{i}}\partial{x_{j}}} &= \sum_{n=1}^{3} \gamma _{n}\delta_{in} \delta_{jn}, \end{aligned}$$

(B.31)

where

$$\begin{aligned} \gamma _{n} &= \frac{a_{1}a_{2}a_{3}}{2} \int_{0}^{\infty}\frac{dt}{(a_{n}^{2}+t)\sqrt {(a_{1}^{2}+t)(a_{2}^{2}+t)(a_{3}^{2}+t)}}. \end{aligned}$$

(B.32)

Setting $t=a_{3}^{2} s$, $\gamma _{n}=\mathcal{E}(\varepsilon _{n};\varepsilon _{1},\varepsilon _{2})$ which is defined in (4.22), with $\varepsilon _{n}=a_{3}/a_{n}$. Mura [92] writes (B.29) in the form

$$\begin{aligned} \varGamma(\mathbf{x}) &= \frac{1}{8\pi} \Biggl(-I+\sum _{n=1}^{3} x_{n}^{2} I_{n} \Biggr), \end{aligned}$$

(B.33)

where

$$\begin{aligned} I = 2\pi a_{1}a_{2}a_{3} \int_{0}^{\infty} \frac{ds}{\Delta(s)}, \qquad I_{n} = 2\pi a_{1}a_{2}a_{3} \int_{0}^{\infty} \frac{ds}{(a_{n}^{2}+s)\Delta(s)}. \end{aligned}$$

(B.34)

The link between $I_{n}$ and ℰ is then clear:

$$\begin{aligned} I_{n} &= 4\pi\mathcal{E}(\varepsilon _{n};\varepsilon _{1}, \varepsilon _{2}). \end{aligned}$$

(B.35)

Finally, note that the integrals $I_{n}$ (or equivalently $\mathcal{E}(\varepsilon _{n};\varepsilon _{1},\varepsilon _{2})$) can be expressed in terms of elliptic integrals [92]. In particular assuming that $a_{1}>a_{2}>a_{3}$,

$$\begin{aligned} I_{1} &= \frac{4\pi \varepsilon _{2}}{(\varepsilon _{2}^{2}/\varepsilon _{1}^{2}-1)(1-\varepsilon _{1}^{2})^{1/2}} \bigl\{ F(\theta,k)-E(\theta,k) \bigr\} , \end{aligned}$$

(B.36)

$$\begin{aligned} I_{3} &= \frac{4\pi}{(1-\varepsilon _{2}^{2})(1-\varepsilon _{1}^{2})^{1/2}} \bigl\{ \bigl(1-\varepsilon _{1}^{2} \bigr)^{1/2}-\varepsilon _{2} E(\theta,k) \bigr\} , \end{aligned}$$

(B.37)

$$\begin{aligned} I_{2} &= 4\pi-I_{1}-I_{3}, \end{aligned}$$

(B.38)

where

$$\begin{aligned} F(\theta,k) &= \int_{0}^{\theta}\frac{dx}{(1-k^{2}\sin^{2} x)^{1/2}}, \qquad E(\theta,k) = \int_{0}^{\theta} \bigl(1-k^{2} \sin^{2} x \bigr)^{1/2} \,dx, \end{aligned}$$

(B.39)

$$\begin{aligned} \theta&= \sin^{-1} \bigl(1-\varepsilon _{1}^{2} \bigr)^{1/2}, \qquad\qquad\qquad\ \;\, k = \frac{1}{\varepsilon _{2}} \biggl( \frac{\varepsilon _{2}^{2}-\varepsilon _{1}^{2}}{1-\varepsilon _{1}^{2}} \biggr)^{1/2}. \end{aligned}$$

(B.40)

2.2 B.2 Closed Integral Form for $\varPsi(\mathbf{x})$

One can also derive an integral form for $\varPsi(\mathbf{x})$ although such a derivation is rather lengthy. Expression (11.38) of Mura [92] is employed, which establishes that (no sum over $i$ here, with sums being shown explicitly for clarity)

$$\begin{aligned} \frac{\partial{\varPsi(\mathbf{x})}}{\partial{x_{i}}} &= \frac {x_{i}}{8\pi} \Biggl( \bigl(I- a_{i}^{2} I_{i} \bigr) - \sum_{n=1}^{3} (I_{n}-I_{in} )x_{n}^{2} \Biggr), \end{aligned}$$

(B.41)

noting the additional factor of $1/(4\pi)$ here from our modified definition of potentials as compared with Mura. The integral $I_{mn}$ is defined as

$$\begin{aligned} I_{mn} &= 2\pi a_{m}^{2} a_{1}a_{2}a_{3} \int_{0}^{\infty} \frac {ds}{(a_{m}^{2}+s)(a_{n}^{2}+s)\Delta(s)}, \end{aligned}$$

(B.42)

where the slight modification to Mura’s notation should be noted: $I_{mn}=a_{m}^{2}I_{mn}^{M}$, where $I_{mn}^{M}$ is Mura’s definition of this integral as defined in Chap. 11 of [92]. This modification means that expressions are now defined in terms of non-dimensional quantities here. In particular it is possible to write $I_{mn}$ in terms of $I_{m}, I_{n}$ and $\epsilon _{m}, \epsilon _{n}$ as follows,

$$\begin{aligned} I_{mn} = \frac{(I_{n}-I_{m})}{(1-(\varepsilon _{m}/\varepsilon _{n})^{2})}, \quad m \neq n. \end{aligned}$$

(B.43)

Additional relations are noted as

$$\begin{aligned} I_{mn} &= \frac{\epsilon _{n}^{2}}{\epsilon _{m}^{2}}I_{nm} \end{aligned}$$

(B.44)

and

$$\begin{aligned} &I_{11} = \frac{4\pi}{3}-\frac{1}{3}(I_{12}+I_{13}), \qquad I_{22} = \frac{4\pi}{3}-\frac{1}{3}(I_{21}+I_{23}), \qquad I_{33} = \frac{4\pi }{3}-\frac{1}{3}(I_{31}+I_{32}), \end{aligned}$$

(B.45)

$$\begin{aligned} &3I_{1} = 3I_{11}+I_{21}+I_{31}, \qquad\quad3I_{2} = 3I_{22} + I_{12} + I_{32}, \qquad\quad\, 3I_{3} = 3I_{33} + I_{13}+I_{23}. \end{aligned}$$

(B.46)

Therefore, differentiating (B.41) with respect to $x_{j}, x_{k}$ and then $x_{\ell}$ (in that order) gives

$$\begin{aligned} \frac{\partial ^{4} \varPsi}{ \partial x_{i}\partial x_{j}\partial x_{k}\partial x_{\ell}} &= \frac {1}{4\pi }\delta_{ij} \delta_{k\ell}(I_{ik}-I_{k})+\frac{1}{4\pi}( \delta_{ik}\delta_{j\ell}+\delta_{jk} \delta_{i\ell}) (I_{ij}-I_{j}). \end{aligned}$$

(B.47)

Note that (B.47) is however not a fully symmetric fourth order tensor as it should be since derivatives should be able to be taken in any order. Therefore one can enforce in turn major, minor and then total symmetry [90] to show that (no sum over repeated coefficients)

$$\begin{aligned} \frac{\partial ^{4}\varPsi}{ \partial x_{i}^{4}} = \frac{3}{4\pi}(I_{ii}-I_{i}), \qquad\frac {\partial ^{4}\varPsi}{ \partial x_{i}^{2}\partial x_{j}^{2}} = \frac{1}{8\pi }(I_{ij}+I_{ji}-I_{i}-I_{j}), \quad i \neq j \end{aligned}$$

(B.48)

and odd derivatives are zero.

Appendix C: Cartesian Coordinates, Rotations and Tensors

Cartesian tensors are used throughout this article, referred to a fixed set of Cartesian coordinates $x_{i}, i=1,2,3$. Some of the properties of Cartesian tensors are summarized shortly, in particular those associated with higher order symmetries are discussed. After this a brief review of rotations of Cartesian coordinates is given for completeness and the matrix formulation of tensors is described. Tensors will be stated as $\mathbf{A}$ in invariant notation with components $A_{ij}$ and $A_{ijk\ell}$ for second and fourth order tensors (see below).

3.1 C.1 Second Order Cartesian Tensors

3.1.1 C.1.1 Isotropy

The second order identity tensor is $I_{ij}=\delta_{ij}$ and with $\alpha$ constant, the general second order isotropic tensor is therefore $A_{ij} = \alpha \delta_{ij}$. Its inverse, with components $\tilde{A}_{ij}$ is $\tilde{A}_{ij} = \frac{1}{\alpha }\delta_{ij}$.

3.1.2 C.1.2 Transverse Isotropy

Upon defining the tensor

$$\begin{aligned} \varTheta_{ij} &= \delta_{ij}-\delta_{i3} \delta_{j3}, \end{aligned}$$

(C.1)

a second order transversely isotropic tensor (with symmetry axis $x_{3}$) has the form

$$\begin{aligned} A_{ij} &= \alpha _{1}\varTheta_{ij} + \alpha _{3}\delta_{i3}\delta_{j3}. \end{aligned}$$

(C.2)

The tensor $\varTheta_{ij}$ defined in (C.1) possesses the following properties:

$$\begin{aligned} \varTheta_{ij}=\varTheta_{ji}, \qquad\varTheta_{ik} \varTheta_{kj} = \varTheta_{ij}, \qquad\varTheta_{ij} \varTheta_{ij} = 2. \end{aligned}$$

(C.3)

Using these properties, the inverse of $\mathbf{A}$ has components $\tilde{A}_{ij}$ that can be written

$$\begin{aligned} \tilde{A}_{ij} &= \frac{1}{\alpha _{1}}\varTheta_{ij} + \frac{1}{\alpha _{3}}\delta_{i3}\delta_{j3}. \end{aligned}$$

(C.4)

This concept can easily be generalized to other symmetry axes either by use of rotations of coordinate axes as will be described in Sect. C.3 or by use of the notation $n_{ij}=n_{i} n_{j}$, where $n_{i}$ are the components of the direction vector associated with the axis of symmetry (for the example described above $n_{i}=\delta_{i3}$). The associated generalization of (C.1) is therefore $\theta _{ij}=\delta_{ij}-n_{ij}$.

3.1.3 C.1.3 Orthotropy

A second order orthotropic tensor has components of the form

$$\begin{aligned} A_{ij} &= \alpha _{1}\delta_{i1}\delta_{j1}+ \alpha _{2}\delta_{i2}\delta_{j2} + \alpha _{3} \delta_{i3}\delta_{j3}, \end{aligned}$$

(C.5)

which has as its inverse

$$\begin{aligned} \tilde{A}_{ij} &= \frac{1}{\alpha _{1}}\delta_{i1} \delta_{j1}+\frac{1}{\alpha _{2}}\delta_{i2} \delta_{j2} + \frac{1}{\alpha _{3}}\delta_{i3} \delta_{j3}. \end{aligned}$$

(C.6)

3.2 C.2 Fourth Order Cartesian Tensors

The following tensors are used extensively in elasticity applications. See Walpole [124] for a comprehensive derivation of all associated theory.

3.2.1 C.2.1 Isotropy

First define the following tensors

$$\begin{aligned} I_{ijk\ell}^{1} &= \frac{1}{3}\delta_{ij} \delta_{k\ell}, \end{aligned}$$

(C.7)

$$\begin{aligned} I_{ijk\ell}^{2} &= \frac{1}{2} (\delta_{ik} \delta_{j\ell}+\delta_{i\ell}\delta_{jk} )- \frac{1}{3} \delta_{ij}\delta_{k\ell}, \end{aligned}$$

(C.8)

$$\begin{aligned} I_{ijk\ell} &= \frac{1}{2} (\delta_{ik} \delta_{j\ell}+ \delta_{i\ell}\delta_{jk} )= I_{ijk\ell}^{1}+I_{ijk\ell}^{2}. \end{aligned}$$

(C.9)

These have the following properties

$$\begin{aligned} I^{1}_{ijk\ell} = I^{1}_{k\ell ij}, \qquad I^{2}_{ijk\ell} = I^{2}_{k\ell ij}, \qquad I_{ijk\ell} = I_{k\ell ij} \end{aligned}$$

(C.10)

and

$$\begin{aligned} I^{1}_{ijmn}I^{1}_{nmk\ell} = I^{1}_{ijk\ell}, \qquad I^{2}_{ijmn}I^{2}_{nmk\ell } = I^{2}_{ijk\ell}, \qquad I^{1}_{ijmn}I^{2}_{nmk\ell} = 0, \qquad I^{2}_{ijmn}I^{1}_{nmk\ell} = 0. \end{aligned}$$

(C.11)

If $\boldsymbol{\sigma }$ is a second order tensor, whose components are written in the deviatoric/scalar form

$$\begin{aligned} \sigma _{ij} &= \sigma '_{ij}+\frac{1}{3}\sigma \delta_{ij}, \end{aligned}$$

(C.12)

where $\sigma =\sigma _{kk}$, then

$$\begin{aligned} I_{ijk\ell}^{1}\sigma _{k\ell} = \frac{1}{3}\sigma \delta_{ij}, \qquad I_{ijk\ell}^{2}\sigma _{k\ell} = \sigma '_{ij}. \end{aligned}$$

(C.13)

Given a fourth order isotropic tensor with components of the form

$$\begin{aligned} A_{ijk\ell} &= 3\alpha _{1} I_{ijk\ell}^{1} + 2 \alpha _{2} I_{ijk\ell}^{2}, \end{aligned}$$

(C.14)

then due to (C.13)

$$\begin{aligned} A_{ijk\ell} \sigma _{k\ell} &= 3\alpha _{1}\sigma \delta_{ij}+2 \alpha _{2}\sigma '_{ij}. \end{aligned}$$

(C.15)

Introduce a second fourth order isotropic tensor $B_{ijk\ell} = 3\beta _{1} I_{ijk\ell}^{1} + 2\beta _{2} I_{ijk\ell}^{2}$ and then

$$\begin{aligned} A_{ijmn}B_{nmk\ell}=B_{ijmn}A_{nmk\ell} &= 9 \alpha _{1}\beta _{1} I_{ijk\ell}^{1} + 4 \alpha _{2}\beta _{2} I_{ijk\ell}^{2} \end{aligned}$$

(C.16)

and the inverse of $A_{ijk\ell}$, denoted by $\tilde{A}_{ijk\ell}$, such that $A_{ijmn}\tilde{A}_{nmk\ell}=I_{ijk\ell}$ is

$$\begin{aligned} \tilde{A}_{ijk\ell} &= \frac{1}{3\alpha _{1}}I_{ijk\ell}^{1} + \frac{1}{2\alpha _{2}}I_{ijk\ell}^{2}. \end{aligned}$$

(C.17)

3.2.2 C.2.2 Cubic System

Define the cubic tensor

$$\begin{aligned} A_{ijk\ell} &= \alpha _{1} I_{ijk\ell}^{1} + \alpha _{2} I^{2}_{ijk\ell} + \alpha _{3} \delta _{ijk\ell}, \end{aligned}$$

(C.18)

where $\delta_{ijk\ell}$ has the property that $\delta_{ijk\ell}=1$ only if $i=j=k=\ell$ and is zero otherwise. Furthermore

$$\begin{aligned} I_{ijmn}^{1}\delta_{nmk\ell} = I_{ijk\ell}^{1}, \qquad I_{ijmn}^{2}\delta_{nmk\ell } &= \delta_{ijk\ell}-I_{ijk\ell}^{1}. \end{aligned}$$

(C.19)

One can write the inverse of $A_{ijk\ell}$ as

$$\begin{aligned} \tilde{A}_{ijk\ell} &= \frac{\alpha _{2}^{2}+\alpha _{1}\alpha _{3}}{\alpha _{2}^{2}(\alpha _{1}+\alpha _{3})} I_{ijk\ell}^{1} + \frac{1}{\alpha _{2}} I_{ijk\ell}^{2} -\frac{\alpha _{3}}{\alpha _{2}^{2}} \delta _{ijk\ell}. \end{aligned}$$

(C.20)

3.2.3 C.2.3 Transverse Isotropy

The Hill tensor basis shall be employed for transversely isotropic tensors. There are several slight variants on this but the Hill basis is used commonly in the micromechanics literature and so it appears sensible to adopt it here. This basis set enables a fourth order transversely isotropic tensor $A_{ijk\ell}$ to be written in the form

$$\begin{aligned} A_{ijk\ell} &= \sum_{n=1}^{6} \alpha _{n} \mathcal{H}_{ijk\ell}^{n}, \end{aligned}$$

(C.21)

where $\alpha _{n}$ are constants. In general it should be noted that $\alpha _{2}\neq \alpha _{3}$ since contraction of a transversely isotropic tensor with another transversely isotropic tensor does not result in a tensor with $\alpha _{2}=\alpha _{3}$. The basis tensors $\mathcal{H}^{n}_{ijkl}$ are defined by

$$\begin{aligned} &\mathcal{H}^{1}_{ijk\ell} = \frac{1}{2} \varTheta_{ij}\varTheta_{k\ell}, \qquad\mathcal{H}^{2}_{ijk\ell} = \varTheta_{ij}\delta_{k3}\delta_{\ell3}, \qquad \mathcal{H}^{3}_{ijk\ell} = \varTheta_{k\ell} \delta_{i3}\delta_{j3}, \end{aligned}$$

(C.22)

$$\begin{aligned} &\mathcal{H}^{4}_{ijk\ell} = \delta_{i3} \delta_{j3}\delta_{k3}\delta_{\ell3}, \qquad \mathcal{H}^{5}_{ijk\ell} = \frac{1}{2}( \varTheta_{ik}\varTheta_{\ell j}+\varTheta_{i\ell} \varTheta_{kj}-\varTheta_{ij}\varTheta_{k\ell}), \end{aligned}$$

(C.23)

$$\begin{aligned} &\mathcal{H}^{6}_{ijk\ell} = \frac{1}{2}( \varTheta_{ik}\delta_{\ell3}\delta_{j3}+ \varTheta_{i\ell}\delta_{k3}\delta_{j3}+ \varTheta_{jk}\delta_{\ell3}\delta_{i3}+ \varTheta_{j\ell}\delta_{k3}\delta_{i3}), \end{aligned}$$

(C.24)

where $\varTheta_{ij} = \delta_{ij}-\delta_{i3}\delta_{j3}$. The notation ℋ signifies the Hill basis.

Let us define the shorthand notation

$$\begin{aligned} \mathcal{H}^{m} \mathcal{H}^{n} &= \mathcal{H}_{ijpq}^{m} \mathcal{H}_{qpk\ell }^{n} \end{aligned}$$

(C.25)

for contraction between the basis tensors defined in (C.22)–(C.24). The contractions defined in (C.25) are then summarized in Table 2.

Note that it is often useful to write the fourth order isotropic identity tensor and basis tensors in the Hill transversely isotropic basis form:

$$\begin{aligned} I_{ijk\ell} &= \mathcal{H}_{ijk\ell}^{1}+ \mathcal{H}_{ijk\ell}^{4}+\mathcal{H}_{ijk\ell}^{5}+ \mathcal{H}_{ijk\ell}^{6}, \end{aligned}$$

(C.26)

$$\begin{aligned} I_{ijk\ell}^{1} &= \frac{1}{3} \bigl(2 \mathcal{H}_{ijk\ell}^{1}+ \mathcal{H}_{ijk\ell}^{2}+ \mathcal{H}_{ijk\ell}^{3}+ \mathcal{H}_{ijk\ell}^{4} \bigr), \end{aligned}$$

(C.27)

$$\begin{aligned} I_{ijk\ell}^{2} &= \frac{1}{3} \bigl(\mathcal{H}_{ijk\ell}^{1}- \mathcal{H}_{ijk\ell}^{2}-\mathcal{H}_{ijk\ell}^{3}+2 \mathcal{H}_{ijk\ell}^{4} +3\mathcal{H}_{ijk\ell}^{5}+3 \mathcal{H}_{ijk\ell}^{6} \bigr). \end{aligned}$$

(C.28)

This allows the definition of the inverse of the tensor $\mathbf{A}$ as $\tilde{\mathbf{A}}$ in a straightforward manner. Its components are

$$\begin{aligned} \tilde{A}_{ijk\ell} &= \sum_{n=1}^{6} \tilde{\alpha }_{n} \mathcal{H}_{ijk\ell}^{n}, \end{aligned}$$

(C.29)

where

$$\begin{aligned} \tilde{\alpha }_{1} &= \frac{\alpha _{4}}{2\Delta }, \qquad\tilde{ \alpha }_{2} = -\frac{\alpha _{2}}{2\Delta}, \qquad\tilde{\alpha }_{3} = - \frac{\alpha _{3}}{2\Delta }, \end{aligned}$$

(C.30)

$$\begin{aligned} \tilde{\alpha }_{4} &= \frac{\alpha _{1}}{2\Delta }, \qquad\tilde{ \alpha }_{5} = \frac {1}{\alpha _{5}}, \qquad\quad\, \tilde{\alpha }_{6} = \frac{1}{\alpha _{6}} \end{aligned}$$

(C.31)

and where $\Delta = \alpha _{1}\alpha _{4}/2-\alpha _{2}\alpha _{3}$.

3.2.4 C.2.4 Hill’s Shorthand Notation

Hill introduced a convenient short-hand notation regarding fourth order tensors. For a fourth order isotropic tensor $A_{ijk\ell}$ defined via

$$\begin{aligned} A_{ijk\ell} = \sum_{n=1}^{2} \alpha _{n} I^{n}_{ijk\ell}, \end{aligned}$$

(C.32)

Hill denoted it in shorthand notation as $A = (\alpha _{1},\alpha _{2} )$, i.e., the tensor basis is assumed from the outset so only the coefficients are to be prescribed. Commonly for linear isotropic elasticity, the case when $\alpha _{1}=3\kappa$, $\alpha _{2}=2\mu$ is assumed where $\kappa$ and $\mu$ are as usual the bulk and shear moduli. For a fourth order transversely isotropic tensor $A_{ijk\ell}$ defined via

$$\begin{aligned} A_{ijk\ell} = \sum_{n=1}^{6} \alpha _{n} \mathcal{H}^{n}_{ijk\ell}, \end{aligned}$$

(C.33)

Hill denoted it in shorthand notation as $A = (\alpha _{1},\alpha _{2},\alpha _{3},\alpha _{4},\alpha _{5},\alpha _{6} )$. Commonly for linear transversely isotropic elasticity, the case when $\alpha _{1}=2K$, $\alpha _{2}=\alpha _{3}=\ell$, $\alpha _{4}=n$, $\alpha _{5}=2m$, $\alpha _{6}=2g$ is assumed.

3.2.5 C.2.5 Orthotropy

Orthotropic basis tensors are described in the paper by Walpole [121] for example. For practical purposes the matrix formulation of fourth order tensors is extremely useful for higher order tensors, including the orthotropic case. This formulation is described in C.4 below.

3.3 C.3 Cartesian Coordinates and Tensors in Rotated Frames

3.3.1 C.3.1 Rotations of Cartesian Coordinates

Consider a fixed Cartesian coordinate system $x_{i}, i=1,2,3$ and an associated Cartesian coordinate system $x_{i}', i=1,2,3$ having the same origin, having been rotated arbitrarily in three dimensions. The general three dimensional rotation matrix is constructed as a product of three planar rotation matrices, each of which corresponds to a rotation of the axes about a given axis in three dimensional space, see Auld [7], p. 77. Begin by rotating clockwise about the $x_{3}$ axis (clockwise, when looking in the positive $x_{3}$ direction) by use of the matrix

$$\begin{aligned} \bigl[Q^{1} \bigr](\varphi) &= \left ( \textstyle\begin{array}{c@{\quad}c@{\quad}c} \cos\varphi& \sin\varphi& 0 \\ -\sin\varphi& \cos\varphi& 0 \\ 0 & 0 & 1 \end{array}\displaystyle \right ), \end{aligned}$$

(C.34)

having an associated (Cartesian) tensor with components $Q_{ij}^{1}$ and the $ij$th component is the $ij$th component of the matrix $[Q^{1}]$. This generates the rotated coordinate system $\bar{x}_{i} = Q_{ij}^{1}(\varphi)x_{j}$. This is followed by a clockwise rotation about the $\bar{x}_{2}$ axis by application of the matrix

$$\begin{aligned} \bigl[Q^{2} \bigr](\vartheta) &= \left ( \textstyle\begin{array}{c@{\quad}c@{\quad}c} \cos\vartheta& 0 & -\sin\vartheta\\ 0 & 1 & 0 \\ \sin\vartheta& 0 & \cos\vartheta \end{array}\displaystyle \right ), \end{aligned}$$

(C.35)

having an associated tensor with components $Q^{2}_{ij}$. This generates the rotated coordinate system $\hat{x}_{i} = Q_{ij}^{2}(\vartheta)Q_{jk}^{1}(\varphi)x_{k}$. Finally a clockwise rotation about the $\hat{x}_{3}$ axis is performed by application of the matrix $[Q^{1}](\psi)$, so that our required fully rotated system in three dimensions is derived as

$$\begin{aligned} x_{i}' &= Q_{ij}^{1}( \psi)Q^{2}_{jk}(\vartheta)Q^{1}_{k\ell}( \varphi)x_{\ell} = \mathcal{Q}_{i\ell}(\varphi,\vartheta, \psi)x_{\ell}. \end{aligned}$$

(C.36)

The domains of the angles of rotation are $\vartheta\in[0,\pi], \varphi\in[0,2\pi), \psi\in[0,2\pi)$. A $3\times3$ matrix $[\mathcal{Q}]$ is associated with the second order tensor $\boldsymbol{\mathcal{Q}}$, having components $\mathcal{Q}_{ij}$.

3.3.2 C.3.2 Cartesian Tensors in Rotated Frames and Orientational Averaging

Employing tensor product notation ⊗, a second order Cartesian tensor $\mathbf{A}$ with components $A_{ij}$ can be written

$$\begin{aligned} \mathbf{A} = A_{ij}\mathbf{e}_{i}\otimes \mathbf{e}_{j}, \end{aligned}$$

(C.37)

where $\mathbf{e}_{j}$ is the $j$th Cartesian unit basis vector. Similarly a fourth order tensor $\mathbf{A}$ is defined as

$$\begin{aligned} \mathbf{A} = A_{ijk\ell}\mathbf{e}_{i}\otimes \mathbf{e}_{j}\otimes\mathbf{e}_{k}\otimes \mathbf{e}_{\ell}. \end{aligned}$$

(C.38)

The relation between components of second and fourth order tensors in the rotated system $x_{i}'$, generated by application of the general rotation matrix $\boldsymbol{\mathcal{Q}}$ to the tensors in the original system $x_{i}$, i.e., $A_{ij}'$ and $A_{ijk\ell}'$ are defined via the following expressions

$$\begin{aligned} A_{ij}' &= \mathcal{Q}_{ik} \mathcal{Q}_{j\ell}A_{k\ell},\qquad\qquad\qquad\quad A_{ij} = \mathcal{Q}_{ki}\mathcal{Q}_{\ell j}A_{k\ell}', \end{aligned}$$

(C.39)

$$\begin{aligned} A_{ijk\ell}' &= \mathcal{Q}_{im} \mathcal{Q}_{jn}\mathcal{Q}_{kp}\mathcal{Q}_{\ell q}A_{mnpq}, \qquad A_{ijk\ell} = \mathcal{Q}_{mi}\mathcal{Q}_{nj} \mathcal{Q}_{pk}\mathcal{Q}_{q\ell}A_{mnpq}'. \end{aligned}$$

(C.40)

Often it is useful to determine the average of a general second order tensor, say $\mathbf{A}$ over all possible rotations of Cartesian axes (uniformly). The natural way to do this is to take a tensor with components $A_{ij}'$, diagonal in some coordinate system $x_{i}'$. This frame has been rotated from the fixed system $x_{i}$. However it is natural to work in a “fixed” coordinate system $x_{i}$, which can be considered as being obtained from the system $x_{i}'$ via a rotation; indeed $x_{i}= \mathcal{Q}_{ji}x_{j}'$. The components of $A_{ij}'$ are diagonal and therefore the components $A_{ij}$ are dependent on the angles $\varphi$, $\vartheta$ and $\psi$ as in (C.36). If one wishes to determine the average of the tensor $A_{ij}$ over all such orientations uniformly, it can be done by carrying out the following integration

$$\begin{aligned} \underline{A}_{ij} &= \frac{1}{8\pi^{2}} \int_{0}^{2\pi} \int_{0}^{2\pi} \int_{0}^{\pi} A_{ij}(\varphi,\vartheta, \psi)\sin\vartheta \,d\vartheta\, d\varphi\,d\psi \\ &= \frac{1}{8\pi^{2}} \int_{0}^{2\pi} \int_{0}^{2\pi} \int_{0}^{\pi} \mathcal{Q}_{ki}(\varphi, \vartheta,\psi)\mathcal{Q}_{\ell j}(\varphi,\vartheta, \psi)A_{k\ell}' \sin\vartheta \,d\vartheta\,d\varphi\,d \psi, \end{aligned}$$

(C.41)

where the underline denotes orientation averaging. Alternatively there may be some orientation distribution function, say $p(\varphi,\vartheta,\psi)$ that weights the importance of certain distributions. A weighted orientation average can then be defined as

$$\begin{aligned} \overline{\underline{A}}_{ij}(p) &= \frac{1}{8\pi^{2}} \int_{0}^{2\pi} \int_{0}^{2\pi} \int_{0}^{\pi}p(\varphi,\vartheta,\psi) \mathcal{Q}_{ki}(\varphi,\vartheta,\psi)\mathcal{Q}_{\ell j}( \varphi,\vartheta,\psi)A_{k\ell}' \sin\vartheta \,d \vartheta\, d\varphi\,d\psi \end{aligned}$$

(C.42)

with $\underline{A}_{ij}=\overline{\underline{A}}_{ij}(1)$. Note the normalization condition on the weighting distribution

$$\begin{aligned} \frac{1}{8\pi^{2}} \int_{0}^{2\pi} \int_{0}^{2\pi} \int_{0}^{\pi}p(\varphi,\vartheta,\psi) \sin\vartheta \,d\vartheta\, d\varphi\, d\psi&= 1. \end{aligned}$$

(C.43)

Analogous expressions to (C.41) and (C.42) hold for the fourth order tensor case of course.

3.3.3 C.3.3 Invariance Properties of Averaging for Second Order Tensors

Weighted orientation averaging of second order tensors can be done mechanically via rotation tensors as described above and as written explicitly in (C.42) or for uniform orientations (C.41). Alternatively for uniform orientation averaging, a simple aspect of tensor analysis associated with invariants can be exploited. Averaging uniformly would give rise to an isotropic tensor with components of the form

$$\underline{A}_{ij} = \alpha \delta_{ij}, $$

where the underline denotes averaging. Performing a contraction in the original tensor gives rise to a quantity that does not change with rotations, i.e., $a=A_{kk}$ is an invariant. Therefore

$$\underline{A}_{kk} = a = \alpha \delta_{kk} = 3\alpha , $$

so that $\alpha =a/3$. The components of the averaged tensor therefore take the form

$$\begin{aligned} \underline{A}_{ij} = \frac{1}{3}a\delta_{ij}. \end{aligned}$$

(C.44)

If averages need to be taken with respect to some weighting function, then the mechanical process of averaging over angles as in (C.42) needs to be followed.

In the case of a transversely isotropic second order tensor $A_{ij}$ with $n_{ij}=\delta_{i3}\delta_{j3}$, it is easily shown that

$$\begin{aligned} \underline{n}_{ij} = \frac{1}{3}\delta_{ij}, \qquad \underline{\varTheta}_{ij} = \frac{2}{3}\delta_{ij}, \end{aligned}$$

(C.45)

and so taking a uniform orientation average of $A_{ij}$ in (C.2) yields

$$\begin{aligned} \underline{A}_{ij} &= \frac{1}{3} (2\alpha _{1} + \alpha _{3} )\delta_{ij}. \end{aligned}$$

(C.46)

In the orthotropic case upon taking uniform averages of (C.5) it is shown that

$$\begin{aligned} \underline{A}_{ij} &= \frac{1}{3}(\alpha _{1}+ \alpha _{2}+\alpha _{3})\delta_{ij}. \end{aligned}$$

(C.47)

3.3.4 C.3.4 Invariance Properties of Averaging for Fourth Order Tensors

In the case of uniform averaging of fourth order tensors, this would give rise to an isotropic tensor with components of the form

$$\begin{aligned} \underline{A}_{ijk\ell} = \alpha I_{ijk\ell}^{1} + \beta I_{ijk\ell}^{2}. \end{aligned}$$

(C.48)

Here $A_{iikk}=a$ and $A_{ikik}=b$ are invariants of the original tensor. Therefore $\underline{A}_{iikk} = a = 3\alpha$ so that $\alpha =a/3$. Furthermore $\underline{A}_{ikik} = b = \alpha + 5\beta$ so that $\beta =b/5-a/15$. The components of the averaged tensor take the form

$$\begin{aligned} \underline{A}_{ijk\ell} = \frac{1}{3}\underline{A}_{ppqq} I_{ijk\ell}^{1} + \frac{1}{15} (3\underline{A}_{pqpq}- \underline{A}_{ppqq} ) I_{ijk\ell}^{2}. \end{aligned}$$

(C.49)

In the case of a transversely isotropic tensor of the form (C.21), taking contractions of the tensor and equating to the averaged tensor in (C.48) yields

$$\begin{aligned} \underline{A}_{ppqq} &= 3\alpha = 2\alpha _{1}+2\alpha _{2}+2 \alpha _{3} + \alpha _{4}, \end{aligned}$$

(C.50)

$$\begin{aligned} \underline{A}_{pqpq} &= \alpha +5\beta = \alpha _{1} + \alpha _{4}+2\alpha _{5}+2\alpha _{6}, \end{aligned}$$

(C.51)

which defines $\alpha$ and $\beta$ in terms of the coefficients $\alpha _{n}$, i.e.,

$$\begin{aligned} \underline{A}_{ijk\ell} = \frac{1}{3}(2\alpha _{1}+2 \alpha _{2}+2\alpha _{3} + \alpha _{4})I^{1}_{ijk\ell} + \frac{1}{15}(\alpha _{1}-2\alpha _{2}-2\alpha _{3}+2 \alpha _{4}+6\alpha _{5}+6\alpha _{6})I^{2}_{ijk\ell}. \end{aligned}$$

(C.52)

3.4 C.4 Matrix Formulation of Tensors

3.4.1 C.4.1 Second Order Tensors

Matrix representation and manipulation of second order Cartesian tensors is of course trivial. One can write the matrix $[A]$ associated with a second order tensor $\mathbf{A}$ with components $A_{ij}$ as

$$\begin{aligned}{} [A] &= \left ( \textstyle\begin{array}{c@{\quad}c@{\quad}c} A_{11} & A_{12} & A_{13} \\ A_{21} & A_{22} & A_{23} \\ A_{31} & A_{32} & A_{33} \end{array}\displaystyle \right ). \end{aligned}$$

(C.53)

Matrix operations can then be carried out, so that, for example, matrix multiplication takes the place of tensor contraction, etc. Furthermore defining the matrix $[\mathcal{Q}]$ associated with the second order general rotation tensor $\mathbf{\mathcal{Q}}$ as we have done below (C.36), the operation of general rotation defined in tensor form in (C.39) can be written in matrix form as

$$\begin{aligned} \bigl[A' \bigr] = [\mathcal{Q}] [A] [\mathcal{Q}]^{T}. \end{aligned}$$

(C.54)

3.4.2 C.4.2 Fourth Order Tensors

It is often of great utility to represent fourth order tensors in the form of six by six matrices, see the reference of Auld [7] for more details. In particular a fourth order tensor $\mathbf{A}$ (satisfying minor symmetries) with components $A_{ijk\ell}$ with respect to a Cartesian basis can be usefully written in matrix form $[A]$ as

$$\begin{aligned}{} [A] &= \left ( \textstyle\begin{array}{c@{\quad}c@{\quad}c@{\quad}c@{\quad}c@{\quad}c} A_{1111} & A_{1122} & A_{1133} & A_{1123} & A_{1131} & A_{1112} \\ A_{2211} & A_{2222} & A_{2233} & A_{2223} & A_{2213} & A_{2212} \\ A_{3311} & A_{3322} & A_{3333} & A_{3323} & A_{3313} & A_{3312} \\ A_{2311} & A_{2322} & A_{2333} & A_{2323} & A_{2313} & A_{2312} \\ A_{1311} & A_{1322} & A_{1333} & A_{1323} & A_{1313} & A_{1312} \\ A_{1211} & A_{1222} & A_{1233} & A_{1223} & A_{1213} & A_{1212} \end{array}\displaystyle \right ). \end{aligned}$$

(C.55)

Of specific importance is the ability to carry out the operations of tensor contraction and inversion with the matrix form stated here. Defining the matrix $[W]$ as

$$\begin{aligned}{} [W] &= \left ( \textstyle\begin{array}{c@{\quad}c@{\quad}c@{\quad}c@{\quad}c@{\quad}c} 1 & 0 & 0 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 & 0 & 0 \\ 0 & 0 & 1 & 0 & 0 & 0 \\ 0 & 0 & 0 & 2 & 0 & 0 \\ 0 & 0 & 0 & 0 & 2 & 0 \\ 0 & 0 & 0 & 0 & 0 & 2 \end{array}\displaystyle \right ), \end{aligned}$$

(C.56)

one can carry out the operation of tensor contraction between two fourth order tensors, i.e., $A_{ijmn}B_{nmk\ell}$, in matrix form as $[A][W][B]$. Furthermore, defining $[A^{-1}]$ as the matrix associated with the inverse of the tensor $\mathbf{A}$ it is straightforward to show that $[A^{-1}] = [W]^{-1}[A]^{-1}[W]^{-1}$.

The matrix $[W]$ can be used in the formulation of matrix forms of the linear elastic constitutive relations, i.e., given the tensor forms

$$\begin{aligned} \sigma _{ij} = C_{ijk\ell}e_{k\ell}, \qquad e_{ij} = D_{ijk\ell} \sigma _{k\ell}, \end{aligned}$$

(C.57)

the equivalent matrix forms are

$$\begin{aligned}{} [\sigma ]=[C] [W] [e], \qquad[e] = [D] [W] [\sigma ], \end{aligned}$$

(C.58)

where $[\sigma ]=(\sigma _{11}\, \sigma _{22} \,\sigma _{33} \,\sigma _{23}\,\sigma _{13}\,\sigma _{12})^{T}$ is the $1\times6$ column vector of stresses and similarly $[e]=(e_{11}\, e_{22} \,e_{33} \,e_{23}\,e_{13}\,e_{12})^{T}$. Multiplying the second of the equations in (C.58) from the left by $[W]$, it follows that

$$\begin{aligned}{} [\sigma ]=[\mathcal{C}] [\gamma ], \qquad[\gamma ] = [\mathcal{D}] [\sigma ], \end{aligned}$$

(C.59)

where

$$\begin{aligned}{} [\mathcal{C}] = [C], \qquad[\mathcal{D}]=[W] [D] [W] = [C]^{-1} \end{aligned}$$

(C.60)

and $[\gamma ]=[W][e]$ is the engineering strain.

General rotations of fourth order tensors are non-trivial operations to define in matrix form. Auld [7] goes into great detail regarding this using Bond’s method [13]. Here we state the result that the operation of rotation about the $x_{1}$ axis as defined above via the rotation matrix $[Q^{1}]$ associated with the fourth order tensor operation $A'_{ijk\ell} = Q^{1}_{im}Q^{1}_{jn}Q^{1}_{kp}Q^{1}_{\ell q}A_{mnpq}$ can be shown to be equivalent to the matrix operation

$$\begin{aligned} \bigl[A' \bigr] = [M](\varphi)[A] [M]^{T}(\varphi), \end{aligned}$$

(C.61)

where

$$\begin{aligned}{} [M](\varphi) &= \left ( \textstyle\begin{array}{c@{\quad}c@{\quad}c@{\quad}c@{\quad}c@{\quad}c} \cos^{2}\varphi& \sin^{2}\varphi& 0 & 0 & 0 & \sin2\varphi\\ \sin^{2}\varphi& \cos^{2}\varphi& 0 & 0 & 0 & -\sin2\varphi\\ 0 & 0 & 1 & 0 & 0 & 0 \\ 0 & 0 & 0 & \cos\varphi& -\sin\varphi& 0 \\ 0 & 0 & 0 & \sin\varphi& \cos\varphi& 0 \\ -\frac{1}{2}\sin2\varphi& \frac{1}{2}\sin2\varphi& 0 & 0 & 0 & \cos 2\varphi \end{array}\displaystyle \right ). \end{aligned}$$

(C.62)

The general rotation operation can be similarly defined via the operation

$$\begin{aligned} \bigl[A' \bigr] = [N](\varphi,\vartheta,\psi)[A] [N]^{T}(\varphi,\vartheta,\psi), \end{aligned}$$

(C.63)

where

$$\begin{aligned}{} [N](\varphi,\vartheta,\psi) &= [M](\psi) \bigl[M' \bigr]( \vartheta)[M](\varphi) \end{aligned}$$

(C.64)

with

$$\begin{aligned} \bigl[M' \bigr](\vartheta) &= \left ( \textstyle\begin{array}{c@{\quad}c@{\quad}c@{\quad}c@{\quad}c@{\quad}c} \cos^{2}\vartheta& 0 & \sin^{2}\vartheta& 0 & -\sin2\vartheta& 0 \\ 0 & 1 & 0 & 0 & 0 & 0 \\ \sin^{2}\vartheta& 0 & \cos^{2}\vartheta& 0 & \sin2\vartheta& 0 \\ 0 & 0 & 0 & \cos\vartheta& 0 & \sin\vartheta\\ \frac{1}{2}\sin2\vartheta& 0 & -\frac{1}{2}\sin2\vartheta& 0 & \cos 2\vartheta& 0 \\ 0 & 0 & 0 & -\sin\vartheta& 0 & \cos\vartheta \end{array}\displaystyle \right ). \end{aligned}$$

(C.65)