## Abstract

An analytic solution for the steady-state temperature distribution in an infinite conductive medium containing an insulated toroidal inhomogeneity and subjected to remotely applied uniform heat flux is obtained. The temperature flux on the torus surface is then determined as a function of torus parameters. This result is used to calculate the resistivity contribution tensor for the toroidal inhomogeneity required to evaluate the effective conductive properties of a material containing multiple inhomogeneities of this shape.

## 1. Introduction

In this paper, we discuss a pore or an insulating inhomogeneity of toroidal shape. Such inhomogeneities occur in both natural and man-made materials. Figure 1 provides several examples: toroidal particles represent the most desirable morphology of Li_{2}O_{2} deposition on a porous carbon electrode in lithium–oxygen batteries [5,6]; polymeric ‘microdonuts’ are used in bioengineering [2]; the toroidal shape of nanoparticles is reported to be preferred for the microwave absorption properties of BaTiO_{3} [3]; Onaka *et al.* [4] reported the formation of toroidal particles of SiO_{2} in a Cu matrix due to internal oxidation of a Cu–Si solid-solution polycrystal.

Analytical modelling of materials with such micro- structure has not been well developed. The inhomogeneities are typically assumed to be ellipsoids of identical aspect ratios. This unrealistic assumption is largely responsible for insufficient linkage between the methods of micromechanics and materials science applications.

While for two-dimensional non-elliptical inhomo- geneities many analytical and numerical results have been obtained [7–10], only a limited number of numerical results and approximate estimates are available for non-ellipsoidal three-dimensional shapes. Several examples of pores of irregular shape that are typical for carbon–carbon composites have been analysed in the context of their contribution to overall compliances by Drach *et al.* [11] using the finite element method. The authors, however, did not discuss the effect of any particular irregularity. In the narrower context of compliance of irregularly shaped cracks, certain results were obtained by Fabrikant [12], Sevostianov & Kachanov [13] (planar cracks), Grechka *et al*. [14] (intersecting planar cracks), Mear *et al*. [15] (non-planar cracks) and Kachanov & Sevostianov [16] (cracks growing from pores and cracks with partial contact between the faces). The effect of irregularly shaped cracks on the overall material resistivity has been discussed by Kachanov *et al*. [17]. Concave inhomogeneity shapes have been analysed by Sevostianov *et al*. [18], Sevostianov & Giraud [19], Chen *et al*. [20] and Sevostianov *et al*. [21] (the last two works analyse the effect of the pores’ concavity on the overall elastic and conductive properties).

For the toroidal shape, certain results have been obtained by Argatov & Sevostianov [22], who found the contribution of a thin rigid toroidal inhomogeneity to overall stiffness; and by Onaka *et al*. [4] and Onaka [23,24], who found the analytical expression for the Eshelby tensor for a toroidal inclusion. We emphasize, regarding the latter works, that the Eshelby tensor for non-ellipsoidal inhomogeneities is irrelevant for the problem of the effective properties of a heterogeneous material—as is clear from the formulation of the two Eshelby problems, of an inhomogeneity and of an eigenstrain (e.g. [20]). Unfortunately, the confusion between the two problems (attempts to calculate the effective properties of materials with non-ellipsoidal inhomogeneities using the Eshelby tensor) can be found in a number of publications.

This work is motivated by the problem of the effective conductivity (thermal or electric) of a material containing toroidal insulating inhomogeneities. The key quantity in this problem is the resistivity contribution tensor that gives the extra temperature gradient produced by the introduction of the inhomogeneity into a material subjected to an otherwise uniform heat flux. This tensor was introduced by Sevostianov & Kachanov [25] in the context of the cross-property connection between the elastic and conductive properties of heterogeneous materials. Following this work, we assume that the background material of volume *V* having the isotropic thermal conductivity *k*_{0} contains an isolated inhomogeneity of volume *V*_{1} of the isotropic thermal conductivity *k*_{1}. The limiting cases *k*_{1}=0 and *T* and the heat flux vector **q** per reference volume (Fourier law) for both constituents, the change in ∇*T* is required to maintain the same heat flux if the inhomogeneity introduced is given by
** R** is the

*resistivity contribution tensor*of the inhomogeneity.

For an insulating inhomogeneity (e.g. a pore), the additional temperature gradient due to the presence of such an inhomogeneity can be represented as an integral over the inhomogeneity boundary
*T* and ** n** are the temperature and outward unit vector normal to the boundary, respectively. Equation (1.2) shows that a Neumann boundary-value problem for the Laplace equation has to be solved for three-dimensional space containing a toroidal insulating inhomogeneity in order to find the resistivity contribution tensor. To the best of our knowledge, this problem has not been solved, although certain steps in the analyses of the toroidal shape have been made. Chapko [26] obtained the numerical solution for the Dirichlet boundary-value problem for the heat equation in the case of a torus. In his PhD thesis, Lauria [27] found the temperature distribution around an insulated torus. Payne & Peel [28], Wakiya [29] and Majumdar & O’Neill [30] found exact solutions to the axisymmetric problem of Stokes flow past a torus. Andrews [31] proposed an alternative separation of the Laplace equation in toroidal coordinates which finds application in electrostatics. In the following text, we obtain an analytic solution for the Neumann problem on temperature distribution in an infinite conductive medium containing an insulated toroidal inhomogeneity and subjected to a remotely applied uniform heat flux. The solution is first obtained in the form of a series expansion of the associated Legendre functions and then we calculate the components of the resistivity contribution tensor. The temperature flux on the torus surface is then determined as a function of torus parameters. The methodology of solving the problem is similar to one used by Radi

*et al*. [32]. The results are used to calculate the resistivity contribution tensor for the toroidal inhomogeneity. As an illustration of the application, we show how to formulate typical homogenization schemes for materials containing multiple toroidal insulating inhomogeneities.

## 2. Formulation of the problem in toroidal coordinates

We consider an insulating circular torus embedded in an infinite conductive medium subject to a constant heat flux at infinity. The corresponding boundary condition on the torus surface requires the normal component of the heat flux to vanish therein. A Cartesian coordinate system (0, *x*, *y*, *z*) is introduced, as shown in figure 2. Following Morse & Feshbach [33] and Lebedev [34], we introduce a toroidal coordinate system (*α*, *β*, *γ*) defined by the following relations:
*a* denotes the distance of the poles from the origin, *α*>0, *β*∈[−*π*,*π*] and *γ*∈[0,2*π*].

Constant values of the bipolar coordinate *α* describe a family of toroidal surfaces. In particular, if *R* and *r* denote the medium radius of the toroidal surface and the radius of its circular section, respectively, then the corresponding coordinate *α*_{0} and length *a* are defined by
*θ* around the circle defined by a constant value of *α* is associated with the toroidal coordinate *β* by the following relations:
*T* under steady-state heat flux satisfies the Laplace equation
*k* is the heat conduction coefficient. The boundary condition on the torus surface then requires
** n** is the outer unit normal on the torus surface (figure 2), coinciding with the unit vector −

*e*_{α}. Therefore, the temperature field must obey the following Neumann boundary condition on the torus surface obtained from equations (2.5) and (2.6):

*α*and

*β*tend to zero), the heat flux must obey the condition

**=**

*q*

*q*_{0}.

## 3. Flux along the symmetry axis

First, we consider the axisymmetric problem corresponding to the heat flux at infinity directed along the symmetry axis *z*:*q*_{0}=(0,0,*q*_{0}). Following Lauria [27], let us split the temperature field into the fundamental contribution *T*_{0} induced by the assigned heat flux in a homogeneous medium and the correction *T*_{1} due to the presence of the toroidal inhomogeneity, namely *T*=*T*_{0}+*T*_{1}. According to (2.1) and (2.5), the contribution *T*_{0} is then given by
*T*_{0} is clearly harmonic, as it depends linearly on the coordinate *z* only; from (2.4) it follows that the field *T*_{1} is also harmonic. By using toroidal coordinates, the most general harmonic function under axisymmetric symmetry conditions has the following representation:
*P*_{n−1/2} is the Legendre function of the first kind [34–36] and the unknown coefficients *B*_{n} can be found by imposing the boundary condition (2.7). Using (A 1), the basic temperature field *T*_{0} in (3.1) can be expanded in Fourier series as
*Q*_{n−1/2} denotes the Legendre function of the second kind. Substitution of the representations (3.2) and (3.3) of the temperature fields into the boundary condition (2.7) yields
*B*_{0}=0 has been assumed. Multiplying equation (3.9) by *p*′_{n} and then summing from *n*=1 to *n*=*m*, we reduce the problem to the following bidiagonal infinite system of algebraic equations:
*m*=1 to *m*=*n*−1 in each term of (3.13), we obtain

A necessary condition for the series (3.2) to be convergent is that the *n*th term of this series goes to zero as *n* goes to *q*′_{n}/*p*′_{n}→e^{−2nα0} as *n* goes to *ξ* follows from (3.16) as
*B*_{1} can be retrieved from (3.14). Correspondingly, the introduction of (3.17) in (3.15) yields

Once the unknown coefficients *B*_{n} have been determined, the total temperature field *T* is given by the sum of equations (3.1) and (3.2), and the corresponding heat flux can be calculated from the Fourier law (2.5) as

The corresponding temperature distribution and heat flux are illustrated in figure 3 for *r*/*R*=0.5. Figure 4 shows the variation in temperature and tangential flux on the torus surface (*α*=*α*_{0}) with the angular coordinate *θ* for various torus aspect ratios *r*/*R*.

## 4. Applied flux orthogonal to the symmetry axis

Let us consider the non-axisymmetric problem corresponding to the heat flux directed along the *y*-axis, and thus orthogonal to the symmetry axis: *q*_{0}=(0,*q*_{0},0). Similar to the previous section, we split the temperature field into the basic contribution *T*_{0} induced by the assigned heat flux *q*_{0} in a homogeneous medium and the correction *T*_{1} due to the presence of the toroidal inhomogeneity *T*=*T*_{0}+*T*_{1}. In this case, the contribution *T*_{0} is given by
*T*_{0} is harmonic, as it depends linearly on the coordinate *y* only; therefore, according to (2.4), the field *T*_{1} is also harmonic. In toroidal coordinates, the harmonic function *T*_{1} can be written in the following form:
*C*_{n} can be found by imposing the boundary condition (2.7). By using the result (A 2), the basic temperature field *T*_{0} in (4.1) can be expanded in Fourier series as
*p*′_{n} and *q*′_{n} are given by (3.6) and

Equation (4.4) can be rearranged as

By using the orthogonality property of cosinusoidal functions, the following tridiagonal infinite system of algebraic equations can be derived from (4.5):
*C*_{n} by using the procedure suggested by Kutsenko & Ulitko [38], which yields the following analytical solution:
*C*_{n} have been determined, the total temperature field *T* is given by the sum of equations (4.1) and (4.2), and the corresponding heat flux can be calculated from the Fourier law (2.5) as
*r*/*R*=0.5. Figure 6 shows the variation in temperature and tangential flux on the torus surface (*α*=*α*_{0}) with the angular coordinate *θ* for various torus aspect ratios *r*/*R*.

## 5. Resistivity contribution tensor for a toroidal pore

We now can calculate both components of the resistivity contribution tensor for a toroidal pore using equation (1.2).

Component *R*_{33} can be calculated from the axisymmetric solution given in §3. The flux of temperature on the torus surface *S* is given by the surface integral
** n**=−

*e*_{α}is the outer unit normal on the torus surface:

*T*=

*T*

_{0}+

*T*

_{1}defined by equations (3.1) and (3.2) into (5.4), accounting for relation (2.2), yields

*V**=2

*π*

^{2}

*r*

^{2}

*R*is the volume of the torus and λ=

*r*/

*R*is its aspect ratio.

Component *R*_{11} can be calculated from the non-axisymmetric solution given in §4. In this case, the non-vanishing component of the flux vector (5.1) on the torus surface *S* is
*T*=*T*_{0}+*T*_{1} from (4.1) and (4.2) into (5.10) yields
*r*/*R* are illustrated in figure 7.

We now compare the obtained results with those for the resistivity contribution tensor of an oblate spheroidal insulator with aspect ratio *γ* which has the following form (e.g. [39]):
*stiffness contribution tensor* of a rigid thin torus can be approximated with good accuracy by the corresponding components calculated for a rigid spheroid that has the same radius and the same volume. For an oblate spheroid, it gives the aspect ratio of the spheroid in terms of the aspect ratio of the torus λ=*r*/*R* as follows:
*b*. Substitution of (5.13) into (5.12) and (5.11) gives us components of the resistivity tensor of a spheroid that has the same volume and the same radius as the torus. These components are shown in figure 8 as functions of λ. Comparison with figure 7*b* shows that component *R*_{11} for a torus can be approximated by *R*_{11} for a spheroid only if the torus is very thin (in agreement with the results of [22]). By contrast, *R*_{33} components for the two shapes behave completely differently from each other. In particular, note that for spheroidal inhomogeneity the coefficient *R*_{33} does not vanish as λ decreases, as in this limit the inhomogeneity approaches the shape of thin circular discs, whereas for toroidal inclusion the shape of a thin circular ring is approached as λ is vanishing small and such a shape does not contribute to the three-dimensional effective properties. Generally, we can state that it is impossible to approximate a torus by an oblate spheroid for any finite interval of variation of λ.

Note, however, that the resistivity contribution tensor for a torus is a symmetric second-rank transversely isotropic tensor and, therefore, one-to-one correspondence with a spheroid must exist for such a tensor. Taking this into account, let us try to approximate it by one for a *prolate*spheroid. In this case, to preserve the radius and the volume, the aspect ratio of the spheroid has to be

Figure 8 illustrates a comparison between the components of the resistivity contribution tensors for insulating inhomogeneities of toroidal and prolate spheroidal pores. It can be seen that the lines are quite close to each other for the entire interval of the variation of parameter λ (the difference is less than 2%). This means that the toroidal inhomogeneity can be approximated by the prolate spheroid of aspect ratio (5.14) in the context of the calculation of the effective conductive properties.

## 6. Effective conductivity of a material containing multipletoroidal inhomogeneities

The resistivity contribution tensor serves as the basic building block for the calculation of the overall thermal properties of heterogeneous materials. In the following sections, we illustrate this with an example of randomly oriented toroidal inhomogeneities. For non-randomly distributed inhomogeneities the approach has to incorporate an orientation distribution function. This generalization is straightforward (e.g. [40]).

### (a) Non-interaction approximation

This approximation is accurate at low concentrations of inhomogeneities (dilute limit). However, its applicability is much wider than those of the dilute limit [41]. As shown by Sevostianov & Sabina [42], the non-interaction approximation (NIA) may remain accurate at a volume fraction of inhomogeneities up to 0.15.

If the interaction between the inhomogeneities is neglected, each inhomogeneity can be assumed to be subjected to the same remotely applied thermal flux, their contributions to the change in the temperature gradient can be treated separately and the total ∇*T* is
*c* is the volume fraction of the inhomogeneities and the overall thermal conductivity is

The NIA is used for various approximate schemes that place non-interacting inclusions into some sort of ‘effective environment’ (effective field or effective matrix) so that the effective conductivities predicted by various homogenization schemes are expressed in terms of the same parameter *η*(λ). This can be done in two different ways, as follows.

(i)

*Effective media approaches*, where a representative inhomogeneity, treated as a single one, is placed into a homogeneous material possessing unknown effective properties.(ii)

*Effective field approaches*, where each inhomogeneity, treated as an isolated one, is placed into the unaltered matrix material and the interactions between different inhomogeneities are accounted for by assuming that the external field acting on them differs from the remotely applied one.

The following modifications of this approach have been developed in the literature.

In the *self-consistent* scheme, the inhomogeneity is embedded into material with effective properties that can be found from the solution of a single inclusion problem. This scheme was first employed probably by Clausius [43] and further developed in the works of Bruggeman [44] for effective conductivity; by Kröner [45] for the elastic properties of polycrystals and by Skorohod [46], Hill [47] and Budiansky [48] for the effective elastic properties of matrix composites. For the isotropic distribution of insulating inhomogeneities, this method yields the following formula for the effective electrical conductivity:
*generalized self-consistent* scheme, first proposed by Kerner [49] and then developed by Christensen & Loo [50], differs from the self-consistent one in that it introduces an intermediate layer with the properties of the matrix being between the inhomogeneity and the background homogeneous material. Analytical results for this scheme can be obtained in the case of spherical inhomogeneities only.

The *differential scheme* can be considered as an infinitesimal version of the self-consistent scheme. It was first proposed by Bruggeman [44,51] and was further developed by Vavakin & Salganik [52], McLaughlin [53] and Zimmerman [54]. The scheme assumes that inhomogeneities are incrementally added to the material until the final volume fraction is reached. On each increment, a set of non-interacting inhomogeneities is added to the homogenized material with the properties determined at the previous step. As pointed out by McLaughlin [53], the total concentration of the inhomogeneities introduced into the matrix does not coincide with the volume concentration of the dispersed phase as a certain fraction of the volume where ‘new’ inhomogeneities are placed is already occupied by the ‘old’ ones. For the isotropic distribution of insulating inhomogeneities, the differential scheme gives the effective thermal conductivity as follows:
*K* (thermal or electrical conductivity, elastic compliance, etc.) behaves as a power of 1−*c*,
*K*_{0} is the corresponding property of the dense material. Comparison of formulae (6.7) and (6.8) allows one to extract information about an ‘average’ pore shape using the fitting parameter that may be obtained experimentally.

In the simplest version of the effective field method, the Mori–Tanaka–Benveniste (MTB) scheme, the above-mentioned field is taken as the average over the matrix and is the same for all inhomogeneities [56,57]. In the more advanced Kanaun–Levin (KL) method [58], the effective field can reflect the statistics of the mutual positions of inhomogeneities. For multi-phase composites, fields acting on different families of inhomogeneities are different [59]. For the isotropic distribution of insulating inhomogeneities, the macroscopic thermal conductivity according to the MTB method is

### Remark .

Maxwell’s [60] scheme, which is probably the oldest method of homogenization, equates the far field produced by the considered set of inhomogeneities to the far field produced by a fictitious domain of a certain shape that possesses unknown effective properties. In the case of the isotropic distribution of inhomogeneities, prediction of Maxwell’s scheme coincides with (6.9).

Figure 9 compares the predictions according to the above-mentioned methods of homogenization. As expected, the difference between the various schemes is substantial as the volume fraction of the inhomogeneities increases. Assessment of the accuracy of different schemes is, however, beyond the scope of the present work—it can be done by comparison with experimental or numerical data only.

## 7. Concluding remarks

We solved analytically the problem of the distribution of the temperature field and heat flux around an insulated toroidal inhomogeneity embedded in an isotropic three-dimensional space subjected to an arbitrarily oriented steady heat flux at infinity. The solution is obtained in the form of an infinite series which rapidly converges everywhere except for the singularity point *r*→*R*. The obtained solution is used to construct the resistivity contribution tensor for an insulating toroidal inhomogeneity—the quantity that describes the extra temperature gradient due to the presence of the inhomogeneity. It is shown, in particular, that this tensor can be approximated by one for a *prolate spheroidal* insulator. The resistivity contribution tensor is the key quantity for the calculation of the effective properties of the material. We illustrated this with the simplest example of isotropic distribution of toroidal inhomogeneities. The obtained results can be reformulated for electric conductivity and for the diffusion process (in the electric conductivity problem, temperature and heat flux should be replaced by electric potential and electric current; in the diffusion problem they should be replaced by concentration and diffusion flux).

## Authors' contributions

E.R. and I.S. developed the mathematical models, interpreted the results and wrote the paper. E.R. performed most of the derivations on temperature distribution in consultation with I.S. I.S. performed most of the derivations on the overall resistivity contribution tensor and the overall resistivity in consultation with E.R. Both authors gave final approval for publication.

## Competing interests

The authors do not have any competing interests.

## Funding

This work was supported by the National Group of Mathematical Physics GNFM-IndAM (Italy), FP7 Project TAMER IRSES-GA-2013-610547 and the New Mexico Space Grant Consortium contained in the NASA Cooperative Agreement NNX13AB19A to New Mexico State University.

## Appendix A

The following useful definite integrals involving the Legendre functions and associated Legendre functions of the first and second kinds, denoted as

- Received November 9, 2015.
- Accepted January 29, 2016.

- © 2016 The Author(s)

Published by the Royal Society. All rights reserved.