Effect of a rigid toroidal inhomogeneity on the elastic properties of a composite

An analytical solution is obtained for the problem of an infinite elastic medium containing a rigid toroidal inhomogeneity under remotely applied uniform strain. The traction on the torus surface is determined as a function of torus parameters and strain components applied at infinity. The results are utilized to calculate components of the stiffness contribution tensor of the rigid toroidal inhomogeneity that is required for calculation of the overall elastic properties of a material containing multiple toroidal inhomogeneities. The analytical results are verified by comparison with finite element model calculations.


Introduction
The paper focuses on the problem of a rigid inhomogeneity of toroidal shape embedded in an elastic matrix.Inhomogeneities of this kind occur in both natural and man-made materials.Figure 1 provides several examples.Barium titanate nanotori are used as nonvolatile memory devices, transducers, optical modulators, sensors, and more recently possible energy storage in supercapacitors [1].Toroidal particles represent the preferred morphology of Li 2 O 2 deposition on porous carbon electrodes in lithium-oxygen batteries [2,3].Polymeric ''microdonuts'' are used in bioengineering [4].The toroidal shape of nanoparticles is reported to be preferred for microwave absorption properties of BaTiO 3 [5].Onaka et al. [6] 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 modeling of materials with such microstructure has not been well developed.In the homogenization schemes, the inhomogeneities are usually assumed to be of ellipsoidal shape.This unrealistic assumption is largely responsible for insufficient linkage between methods of micromechanics and materials science applications.While many analytical and numerical results have been obtained for two-dimensional non-elliptical inhomogeneities [10][11][12][13], only a limited number of numerical results and approximate estimates are available for non-ellipsoidal three-dimensional (3D) shapes.Most of them are related to pores and cracks.In the context of compliance of irregularly shaped cracks, certain results have been obtained [14][15][16][17][18]. Various concave pores have also been analyzed [19][20][21][22].
For toroidal shapes, Argatov and Sevostianov [23] used asymptotic methods to evaluate the contribution of a thin rigid toroidal inhomogeneity into overall stiffness.Onaka and colleagues [6,24,25] derived analytical expressions for components of the Eshelby tensor for a toroidal inclusion.We emphasize, however, that the Eshelby tensor for non-ellipsoidal inhomogeneities is irrelevant to the problem of effective properties of a heterogeneous material (as clearly seen from the formulation of the two Eshelby problems; see, for example, Chen et al. [21]).The problem on the effective conductivity (thermal or electric) of a material containing toroidal insulating inhomogeneities has been addressed by Radi and Sevostianov [26].
The present work focuses on the evaluation of the contribution provided by a rigid toroidal inhomogeneity to effective elastic properties.We first consider a homogeneous elastic material (matrix), with the stiffness tensor C 0 assumed to be isotropic.It contains an inhomogeneity, of volume V 1 ð Þ , of a different elastic material with the compliance and stiffness tensor C 1 .The contribution of the inhomogeneity to the overall stress, per representative volume V (the extra stress, as compared to the homogeneous matrix) is given by the fourth-rank tensor Nthe stiffness contribution tensor of the inhomogeneity -defined by where e ' is the ''remotely applied'' strain.For a rigid inhomogeneity, this additional stress field Ds can be represented as integral over the inhomogeneity boundary [27] where n is the outward unit normal to the inhomogeneity surface S at point x.
To calculate components of the stiffness contribution tensor of a rigid toroidal inhomogeneity, a displacement boundary value problem has to be solved for 3D elastic space containing such and inhomogeneity.This problem has been addressed previously [28][29][30].In the text to follow, we modify the previous solution to calculate the components of the stiffness contribution tensor.

Formulation of the problem in toroidal coordinates
We consider a rigid circular torus embedded in an infinite elastic medium subject to remotely applied homogeneous strain field.Following Morse and Feshbach [31] and Lebedev and Silverman [32], we introduce a toroidal coordinate system (a,b,u) (Figure 2) defined by the following relations where c denotes the distance of the poles from the origin, a .0, b2 [2p, p], and j2 [0, 2p].Constant values of the bipolar coordinate a describe a family of toroidal surfaces.R 0 and r 0 denote medium radius of the torus and the radius of its circular cross-section, respectively.Then With reference to an arbitrary meridian section, the poloidal angle c 2 Àp, p ½ around the circle defined by a constant value of a is associated with the toroidal coordinate b by the following relations: Then a torus surface a = a 0 can be defined parametrically as follows: denote the associated Legendre function of semi-integer index and order k of the first or the second kind, respectively [32][33][34].They satisfy the following equations where the operator D k is defined by

Analytical solution for a rigid toroidal inclusion in an infinite elastic medium
The total displacement field u in the infinite elastic medium containing a rigid toroidal inclusion can be represented as a sum u = u ' + u 0 where u ' corresponds to the remotely applied homogeneous strain field and u 0 is the correction due to the presence of the toroidal inhomogeneity.Due to the geometrical symmetry of the problem, the stiffness contribution tensor N is expected to be transversely isotropic.Therefore, it is sufficient to consider the following Cartesian components of the remotely applied strains: e ' zz , e ' yy , and e ' yz .Then, the fundamental displacement field admits the following Cartesian components The corresponding displacement components in cylindrical coordinates can be written as For the general case of arbitrary load, the inhomogeneity may be subjected to rigid body translation and rotation.In this case the boundary conditions on the inhomogeneity surface could be represented in the form of the Robin problem [34] as u = d + v × x, where d and v are vectors of small translation and rotation of the rigid inhomogeneity, and x is a position vector on the interface.In the case of tensile strain e ' zz and e ' yy , the rigid motion of the torus can be ignored due to the symmetry of the applied loads and geometry.In contrast, simple shear deformation e ' yz leads to the rigid rotation of the toroid with respect to its center with d = 0 and v = v e x , where v is the rotation angle of the torus (to be determined) and e x is a unit vector along the x-axis.Then, in the cylindrical coordinate system, In the study by Krokhmal [35], it was suggested to seek for the corrective term u 0 in the following form where g = 1=½4 (1 À n), and the functions q k , f k , c k , and x k satisfy the following relations and the differential constraint with c 21 = 0.The general solutions to equations ( 16) matching the displacement field of equations ( 12) and ( 13) have the following form for k = 0, 1, 2, 3.The primed summation in equation (18) denotes that the first term of the sum corresponding to n = 0 is halved.The unknown coefficients A n,k , B n,k , C n,k , and D n,k can be determined from the boundary condition at a = a 0 , together with the differential constraint (equation ( 17)).By using equations ( 11)-( 14), the boundary conditions at a = a 0 require Introducing the unknowns x n , k for the coefficients B n , k and using equations ( 15) and ( 19), all the coefficients in the Fourier series of equation ( 18) can be written in terms of x n , k as follows where being d 0 = 1/2, e 0 = 0, and d k = e k = 1, for k! 1. Coefficients p n,k , q n,k , and s n,k are the Fourier coefficients for the following functions Using the decomposition given in the Appendix, we finally obtain the following relations and obtain the following infinite system of algebraic equations for the unknown variables x n,k : where for n! 0. Note that for k = 0, 2 equation ( 25) holds for n! 0, while for k = 1 it holds for n! 1.
A necessary condition for the series in equation ( 18) to be convergent is that the coefficients A n,k , B n,k , C n,k and D n,k tend to zero as n !'.Thus, the unknown coefficients x n,k have to converge to zero.In this case, the infinite system of equations ( 24) can be solved explicitly using the procedure for tridiagonal matrix where From the conditions for n = 0 and k = 0, 2 From the conditions for n = 0 and k = 1 Once the unknown coefficients x n,k have been determined, the total displacement field in cylindrical coordinates is given by the sum of equations ( 12), (13), and (14).Then, the corresponding strain field can be calculated from the compatibility condition: where we used connection between derivatives in cylindrical and toroidal coordinate systems and expressions for derivatives of the associated Legendre functions given in the Appendix.Finally, the stress field is obtained from the Hooke's law where E and n are the Young's modulus and Poisson ratio of the isotropic matrix and I is the secondrank unit tensor.The Cartesian components of the stress tensor s can be easily expressed via cylindrical components s xx = s rr cos 2 u + s uu sin 2 u À 2s ru sinu cosu, s yy = s rr sin 2 u + s uu cos 2 u À 2s ru sinu cosu, The traction vector at the torus surface a = a 0 is t =2e a Á s with the outer unit normal The principal moment produced by torus surface tractions is where x = c cosh a 0 À cos b fsinh a 0 cos u, sinh a 0 sin u, sin bg ð 33Þ is the position-vector of a point on the torus surface.The infinitesimal element of torus surface is given as [31] dA = The solution obtained of the boundary value problem for an elastic medium with rigid torus is specified below for three load cases: e ' zz , e ' yy and e ' yz .
(i) Tension along z-axis.The nonzero component of remote elastic strain is e ' zz .In this case the total displacement field has axial symmetry about the z-axis and vanishes on the surface of the torus.The maps for the displacement components u y and u z are shown in Figure 3 for the meridional cross-section of the matrix near the inhomogeneity.
The distribution of Cartesian component of stress vector on the torus surface (a = a 0 ) is shown in Figure 4.
(ii) Tension along y-axis.The nonzero component of remote elastic strain is e ' yy .The total displacement vanishes on the surface of the rigid toroidal inclusion.The distribution of the displacement components u y and u z is illustrated in Figure 5 for the meridional cross-section of the matrix near the inhomogeneity.The distribution of Cartesian component of stress vector on the torus surface (a = a 0 ) is shown in Figure 6.It is seen that the boundary condition are satisfied on rigid torus surface for both (i) and (ii) load cases.(iii) Simple shear in yz plane.The nonzero component of remote elastic strain is e ' yz .Under shear, the toroidal inhomogeneity has rigid rotation on angle v to be determined from the moment equilibrium condition with respect to the inhomogeneity where M v is resultant moment required to produce rotation on angle v and M g is the reactive moment caused by simple shear e ' yz .Moments could be calculated using the equation ( 32) Due to the linearity of the problem, forces and moments are proportional to the displacements and rotations, respectively.Then, accounting for the equilibrium equation ( 34), the rotation angle can be calculated as The dependence of v on the aspect ratio of torus is shown in Figure 7.
The distribution of the displacement components u y and u z under simple shear e ' yz is shown in Figure 8 for the meridional cross-section of the matrix near the torus.
The distribution of Cartesian component of the traction vector on torus surface is shown in Figure 9.We now can calculate the components of the stiffness contribution tensor N for a rigid toroidal inhomogeneity using equations ( 1) and ( 2), namely where V Ã = 2p 2 r 2 0 R 0 is the volume of the torus and dA is the infinitesimal element of torus surface (see equation (34)).In the toroidal coordinate system traction t and position vector x is given by equations ( 29)-( 31) and (34).Then the normalized dependence of the stiffness contribution tensor N on torus aspect ratio l = r 0 /R 0 can be obtained from equation (38).

Finite element model
To verify the analytical solution, we also calculate components of the stiffness contribution tensors using finite element analysis (FEA).To produce the necessary 3D mesh of the considered shape we start by generating surface mesh in a custom MATLAB script [36,37].A torus surface is defined parametrically using equation (6).Since the coordinates are stored in the form of an ordered list they can easily be connected into the continuous mesh with triangular elements.Note that each mesh generated using our script is composed of approximately 30,000 surface elements (see Figure 10).
Torus surface mesh is then placed into a large reference volume V that has cubic shape with sides at least five times bigger than the largest linear dimension of the inhomogeneity to reduce boundary effects and simulate remote loading.This setup is auto meshed with non-linear tetrahedral 3D elements due to higher accuracy of results compared to linear elements.An example of the mesh is given in Figure 11.To simulate perfectly rigid inhomogeneity we assume that the ratio of matrix and inhomogeneity Young's moduli is E i /E = 10 10 .The boundary ∂V of the reference volume I is subject to the kinematic conditions At the next step, to find all nonzero components of tensor N we perform a set of six load cases in FEA package MSC Marc Mentat: three normal loadings in the directions of three global coordinate axes and three shear loads according to equation (39).Once the numerical simulations are completed for the given shape, the result files are processed to calculate all volume-averaged stress components within V from each load case: where s ij m is the volume average of the stress component ij calculated from the mth load case, V is the reference volume, (s (l) ij ) m is the stress component ij at the centroid of the finite element l calculated  from the mth load case, V (l) is the volume of the element l, and N e is the total number of elements in the model.Given the average stress components we then calculate the stiffness contribution tensor from: where (e ' kl ) m are the components of the prescribed strain and (s ' ij ) m are the stress components inside V in the absence of the inhomogeneity.Components of the stiffness contribution tensors are than normalized by particle volume fraction and matrix Young's modulus E as follows, N ijkl = V EV Ã N ijkl , where V * is torus inhomogeneity volume.Figure 12 illustrates comparison of the nonzero components of stiffness contribution tensor calculated analytically using equation (38) with ones obtained by FEA.

Approximation by spheroid
We now compare the obtained results with ones for stiffness contribution tensor of an oblate rigid spheroid of the aspect ratio x components of this tensor can be written as follows: where and the following notation is used: Argatov and Sevostianov [23], using asymptotic methods, showed that in-plane components of the stiffness contribution tensor of a thin rigid 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 l = r 0 =R 0 as follows: Equation (44) gives us a spheroid that has the same volume and the same radius as the torus.These components are shown in Figure 12(a), (c), and (e) as functions of l.
Out-of-plane components of the stiffness contribution tensor of a rigid torus can be approximated by the in-plane components of a rigid prolate spheroid that has the same radius and the same volume.In this case, to preserve the radius and the volume, the aspect ratio of the spheroid has to be Figure 12(b), (d), and (f) illustrate a comparison between the components of the stiffness contribution tensors for rigid inhomogeneities of toroidal and prolate spheroidal pores.It is seen that the lines are quite close to each other with the exception of N 1313 .

Discussion and conclusions
We solved an elasticity boundary-value problem for an infinite elastic medium containing a rigid toroidal inhomogeneity subjected to remotely applied uniform strains.The solution is obtained in explicit form of associated Legendre function's series using the approach proposed by Eroshkin and Tsukrov [27].Displacement and traction is determined as a function of coordinates, torus parameters and strain components applied at infinity.We managed to combine the procedures for solving the axisymmetric (k = 0) and asymmetric (k ! 1) case and write the solution in the most general form.The boundary-value problems were reduced to infinite systems of algebraic equations with three-diagonal matrices.We succeeded in extending the classical method of solving finite algebraic equations with a tridiagonal matrix to the case of an infinite number of equations.
Considering the toroidal inhomogeneity subject to simple shear in the meridional plane, we have taken into account the rigid rotation of the inhomogeneity.The angle of rotation can be calculated from the equation of moment equilibrium.The angle decreases when the torus aspect ratio increases.
The obtained solution is used to get explicit expressions for the components of stiffness contribution tensor of the rigid toroidal inhomogeneity.The results are verified by comparison with finite element simulating.The components of the stiffness contribution tensor for a torus are compared with ones for a rigid spheroid having the same radius and the same volume as torus.It is shown that the components of the stiffness contribution tensor for an oblate spheroid serve as good approximations for corresponding in-plane components of the torus while two out-plane components can be approximated by components of the stiffness contribution tensor for a prolate spheroid.
The expression for stiffness contribution tensor yields the possibility to calculate effective elastic properties of a material containing multiple toroidal inhomogeneities.Indeed, if interaction between the inhomogeneities is neglected, each inhomogeneity can be assumed to be subjected to the same remotely applied strain, their contributions into the change in the stiffness tensor can be treated separately.In none interaction approximation the extra stiffness due to inhomogeneities of the same shape can always be represented in the form where summation over inhomogeneities may be replaced by integration over their orientation.The noninteraction approximation, in turn, serves as the basic building block for various homogenization schemes -self-consistent, differential, Mori-Tanaka, Maxwell, etc. [38].

Figure 2 .
Figure 2. Toroidal coordinate system (a, b, j).Associated Cartesian (x, y, z) and cylindrical (r, u, z) coordinate systems are also shown.The coordinate surfaces a 0 = const are the eccentric family of tori and b 0 = 6const are spherical caps having there centers along the z-axis.c is the so-called poloidal angel of the torus ring.

Figure 3 .
Figure 3.The distribution of the (a) u x = u y and (b) u z displacement components near rigid toroidal inhomogeneity (R 0 /r 0 = 2) in infinite media under tension along the z-axis, i.e. for e ' zz = 1 and e ' yy = g ' yz = 0.The maps are plotted for a meridian plane y = 0.The displacement values are given in units of R 0 for n = 0.3.

Figure 4 .
Figure 4.The distribution of the (a) t x , (b) t y , and (c) t z traction components at the torus surface a = a 0 , R 0 / r 0 = 2 with respect to azimuth u and poloidal angel c.The nonzero component of remote strain is e ' zz = 1.The stress values are given in units of E for n = 0.3.

Figure 5 .
Figure 5.The distribution of the (a) u x , (b) u y , and (c, d) u z displacement components near rigid toroidal inhomogeneity (R 0 / r 0 = 2) in infinite media under tension along the y-axis, i.e. for e ' yy = 1 and e ' zz = g ' yz = 0.The maps are plotted for the meridian planes: (a) and (c) for y = 0; (b) and (d) for x = 0.The displacement values are given in units of R 0 for n = 0.3.

Figure 6 .
Figure 6.The distribution of the (a) t x , (b) t y and (c) t z traction components at the torus surface a = a 0 , R 0 / r 0 = 2 with respect to azimuth u and poloidal angel c.The nonzero component of remote strain is e ' yy = 1.The stress values are given in units of E for n = 0.3.

Figure 7 .
Figure 7. Dependences rotation angle v of inclusion with respect to torus aspect ratio r 0 /R 0 for g ' yz = 1 and different values of Poisson ratio 0.00, 0.30, and 0.45.

Figure 8 .
Figure 8.The distribution of the (a) u y and (b) u z displacement components near rigid toroidal inhomogeneity (R 0 /r 0 = 2) in infinite media under simple shear in the yz plane, i.e. for g ' yz = 1, v ' 0.794, and e ' yy = e ' zz = 0.The maps are plotted for the meridian plane x = 0.The displacement values are given in units of R 0 for n = 0.3.

Figure 9 .
Figure 9.The distribution of the (a) t x , (b) t y , and (c) t z traction components at the torus surface a = a 0 , R 0 / r 0 = 2 with respect to azimuth u and poloidal angel c.The nonzero component of remote strain is g ' yz = 1.The stress values are given in units of E for n = 0.3.

Figure 11 .
Figure 11.Example of a mesh density of the matrix and torus with r 0 /R 0 = 0.5: (a) general view; (b) close-up view of the highlighted region.

Figure 12 .
Figure 12.Normalized nonzero components of the stiffness contribution tensor of inhomogeneity (a) N 1111 = N 2222 , (b) N 3333 , (c) N 1122 , (d) N 1133 , (e) N 1212 , and (f) N 1313 = N 2323 as function of the aspect ratio r 0 /R 0 .The solid curves correspond to the rigid toroidal inclusion; The dashed curves correspond to the rigid ellipsoidal inclusion; the FEM results are noted by marks.The stiffness values are given in units of e ij E for different values of Poisson ratio 0.0, 0.3, and 0.45. x