## Abstract

Eshelby’s problem of an ellipsoidal inclusion embedded in an infinite homogeneous isotropic elastic material and prescribed with a uniform eigenstrain and a uniform eigenstrain gradient is analytically solved. The solution is based on a simplified strain gradient elasticity theory (SSGET) that contains one material length scale parameter in addition to two classical elastic constants. The fourth-order Eshelby tensor is obtained in analytical expressions for both the regions inside and outside the inclusion in terms of two line integrals and two surface integrals. This non-classical Eshelby tensor consists of a classical part and a gradient part. The former involves Poisson’s ratio only, while the latter includes the length scale parameter additionally, which enables the newly obtained Eshelby tensor to capture the inclusion size effect, unlike its counterpart based on classical elasticity. The accompanying fifth-order Eshelby-like tensor relates the prescribed eigenstrain gradient to the disturbed strain and has only a gradient part. When the strain gradient effect is not considered, the new Eshelby tensor reduces to the classical Eshelby tensor, and the Eshelby-like tensor vanishes. In addition, the current Eshelby tensor for the ellipsoidal inclusion problem includes those for the spherical and cylindrical inclusion problems based on the SSGET as two limiting cases. The non-classical Eshelby tensor depends on the position and is non-uniform even inside the inclusion, which differ from its classical counterpart. For homogenization applications, the volume average of the new Eshelby tensor over the ellipsoidal inclusion is analytically obtained. The numerical results quantitatively show that the components of the newly derived Eshelby tensor vary with both the position and the inclusion size, unlike their classical counterparts. When the inclusion size is small, it is found that the contribution of the gradient part is significantly large. It is also seen that the components of the averaged Eshelby tensor change with the inclusion size: the smaller the inclusion, the smaller the components. Moreover, these components are observed to approach the values of their classical counterparts from below when the inclusion size becomes sufficiently large.

## 1. Introduction

Eshelby’s (1957, 1959) solution for an ellipsoidal inclusion in an infinite homogeneous isotropic elastic material has been playing a key role in micromechanics (e.g. Mura 1982; Nemat-Nasser & Hori 1999; Qu & Cherkaoui 2006; Li & Wang 2008), and Eshelby’s tensor for strain transformation is crucial for homogenization schemes including the Mori–Tanaka method (e.g. Mori & Tanaka 1973; Weng 1990; Genin & Birman 2009) and the self-consistent method (e.g. Budiansky 1965; Hill 1965; Huang *et al*. 1994; Le Quang & He 2007). The solution of the more complex dynamic Eshelby ellipsoidal inclusion problem has been obtained by Michelitsch *et al*. (2003), which reduces to the solution of Eshelby (1957) in the static limiting case.

However, Eshelby’s (1957, 1959) original formulation and most subsequent studies (e.g. Markenscoff 1998; Ju & Sun 1999; Wang *et al*. 2005; Liu 2008) are based on classical elasticity, and the resulting Eshelby tensor inside an ellipsoidal inclusion depends only on Poisson’s ratio of the material and the aspect ratios of the inclusion. As a result, homogenization methods based on this classical Eshelby tensor cannot capture the inclusion (particle) size effect on elastic properties exhibited by particle-matrix composites (e.g. Vollenberg & Heikens 1989; Cho *et al*. 2006; Marcadon *et al*. 2007; Gao 2008). This has motivated studies on Eshelby-type inclusion problems using higher order elasticity theories, which contain material length scale parameters and can describe size-dependent elastic deformations of polymers and other materials (e.g. Nikolov *et al*. 2007).

The higher order elasticity theories that have been used in studying Eshelby-type inclusion problems include a micropolar theory (Cheng & He 1995, 1997; Ma & Hu 2006), a microstretch theory (Kiris & Inan 2006; Ma & Hu 2007), a modified couple stress theory (Zheng & Zhao 2004) and a simplified strain gradient theory (Gao & Ma 2009, 2010; Ma & Gao 2010). Most solutions obtained in these studies are for the problem of a spherical or cylindrical inclusion embedded in an infinite elastic medium of non-Cauchy type. For the ellipsoidal inclusion problem, the two studies by Ma & Hu (2006, 2007) can be identified, which consider spheroidal inclusions and make use of an important result provided in Michelitsch *et al*. (2003). However, the micropolar theory used in Ma & Hu (2006) contains four material parameters in addition to the two Lamé constants, and the microstretch theory applied in Ma & Hu (2007) involves seven additional material parameters. In view of the difficulties in determining these microstructure-dependent material length scale parameters (e.g. Lakes 1995; Lam *et al*. 2003) and in dealing with fourth-order Eshelby tensors, there is a need to solve the Eshelby ellipsoidal inclusion problem using a higher order elasticity theory that contains only one additional material length scale parameter.

The aim of the current paper is to provide such a solution by using a simplified strain gradient elasticity theory (SSGET), which involves only one material length scale parameter in addition to two classical elastic constants (Gao & Park 2007). This theory has recently been used by Gao & Ma (2009, 2010) and Ma & Gao (2010) to derive closed-form expressions of the Eshelby tensor for the spherical and cylindrical inclusion problems.

The rest of this paper is organized as follows. In §2, the Eshelby tensor for a three-dimensional inclusion of arbitrary shape is presented in a general form using the Green’s function method. In §3, analytical expressions of the Eshelby tensor inside and outside an ellipsoidal inclusion with three distinct semi-axes are derived, and the volume average of the Eshelby tensor over the ellipsoidal inclu- sion is analytically evaluated. Numerical results are provided in §4 to quantita- tively illustrate the position dependence and the inclusion size dependence of the newly obtained Eshelby tensor for the ellipsoidal inclusion problem.

## 2. Eshelby tensor based on the simplified strain gradient elasticity theory

The SSGET, also known as the first gradient elasticity theory of Helmholtz type and the dipolar gradient elasticity theory (Gao & Ma 2010), uses the following strain energy density function, *w*, for an isotropic, linearly elastic material:
2.1
where *λ* and *μ* are the Lamé constants in classical elasticity, *L* is a material length scale parameter, and *ε*_{ij} and *κ*_{ijk} are, respectively, the components of the infinitesimal strain, , and the strain gradient, , given by
2.2
with *u*_{i} being the displacement components. Physically, the dependence of *w* on **∇ ε** included in equation (2.1) arises from the non-local nature of atomic forces (Kröner 1963).

It follows from equation (2.1) that the constitutive equations are
2.3and
2.4
where *τ*_{ij} are the components of the Cauchy stress, , that is energetically conjugated to the infinitesimal strain ** ε**, and

*μ*

_{ijk}are the components of the double stress, , that is energetically conjugated to the strain gradient

**κ**. Physically, the double stress

**μ**corresponds to double (or dipolar) forces per unit area (e.g. Lazar & Maugin 2005; Vardoulakis & Giannakopoulos 2006). Multipolar forces are anti-parallel forces acting in a microstructured continuum (Green & Rivlin 1964). Equations (2.3) and (2.4) can be combined to obtain the modified constitutive relations: 2.5 where

*σ*

_{ij}are the components of the total stress, , which is clearly symmetric with

*σ*

_{ij}=

*σ*

_{ji}. The total stress

**σ**is a balance (equilibrium) stress, and the expression of

*σ*

_{ij}as a combination of

*τ*

_{ij}and

*μ*

_{ijk}given in equation (2.5) arises naturally from the variation of the strain energy (Gao & Park 2007).

As shown in Gao & Park (2007), the Navier-like displacement equations of equilibrium read
2.6
and the boundary conditions have the form:
2.7
where *f*_{i} are the components of the body force vector (force per unit volume), and *t*_{i} and *q*_{i} are, respectively, the components of the Cauchy traction vector and double stress traction vector. In equations (2.6) and (2.7), *Ω* is the region occupied by the elastically deformed material, ∂*Ω* is the smooth bounding surface of *Ω*, and *n*_{i} is the outward unit normal vector on ∂*Ω*. In equations (2.7) and in the sequel, the overbar represents the prescribed value. Note that the standard index notation, together with the Einstein summation convention, is used in equations (2.1)–(2.7) and throughout this paper, with each Latin index (subscript) ranging from 1 to 3 unless otherwise stated.

Equations (2.6) and (2.7) define the displacement form of the SSGET. Clearly, the material length scale parameter *L* is explicitly involved in equation (2.6) in addition to the two Lamé constants *λ* and *μ*. This parameter *L* can be estimated from comparing the dispersion curves of Rayleigh waves obtained using the SSGET and those from lattice dynamics calculations. It can also be estimated by fitting the predicted results based on the SSGET to the measured data from bending and torsion tests of microstructured solids. The SSGET described by equations (2.6) and (2.7) and containing the length scale parameter *L* has been used to derive microstructure-dependent solutions of a number of important boundary value problems, which are reviewed in Gao & Ma (2010) along with other relevant issues. When the strain gradient effect is absent (i.e. *L*=0), *μ*_{ijk}=0 (see equation (2.4)) and *σ*_{ij}=*τ*_{ij} (see equations (2.5)), and equations (2.6) and (2.7) reduce to the governing equations and boundary conditions in terms of displacement in classical elasticity (e.g. Gao & Rowlands 2000).

Solving equation (2.6), subject to the boundary conditions of *u*_{i} and their first, second and third derivatives vanishing at infinity, gives the fundamental solution and Green’s function based on the SSGET. By using three-dimensional Fourier and inverse Fourier transforms, the fundamental solution of equation (2.6) has been obtained as (Gao & Ma 2009)
2.8
where **x** is the position vector of a point in the infinite three-dimensional space, **y** is the integration point, and *G*_{ij}(⋅) is the Green’s function (a second-order tensor) given by
2.9with
2.10
It can be readily shown that when *L*=0, equations (2.9) and (2.10) reduce to the Green’s function for three-dimensional problems in classical elasticity (e.g. Li & Wang 2008). In equations (2.9) and (2.10), *v* is Poisson’s ratio, which is related to the Lamé constants *λ* and *μ* through (e.g. Timoshenko & Goodier 1970)
2.11
where *E* is Young’s modulus.

By using the Green’s function method built on equations (2.8)–(2.10), the general expressions of the Eshelby tensor based on the SSGET can then be obtained, as summarized below.

Consider the problem of a three-dimensional inclusion of arbitrary shape embedded in an infinite homogenous isotropic elastic body, as shown in figure 1. The inclusion is prescribed with a uniform eigenstrain ** ε*** and a uniform eigenstrain gradient

**κ***, which can be independent of each other (Ma & Gao 2010). There is no body force or any other external force acting on the elastic body.

The disturbed strain, *ε*_{ij}, induced by ** ε*** and

**κ*** can be shown to be (Gao & Ma 2009) 2.12 where

*S*

_{ijlm}is the fourth-order Eshelby tensor having 36 independent components (or even fewer; Yin

*et al*. 2006), and

*T*

_{ijlmp}is a fifth-order Eshelby-like tensor with 108 independent components. Equation (2.12) shows that

**is solely linked to**

*ε**** in the absence of**

*ε***κ*** (i.e. the classical case) and is fully related to

**κ*** in the absence of

***. The fourth-order Eshelby tensor can be obtained as 2.13a 2.13b 2.13c where is the classical part, is the gradient part,**

*ε**δ*

_{ij}is the Kronecker delta, and 2.14a,b,c are three scalar-valued functions that can be obtained analytically or numerically by evaluating the volume integrals over the domain

*Ω*occupied by the inclusion, with

**y**(∈

*Ω*) being the integration variable. It can be shown that

*Φ*(

**x**),

**(**

*Λ***x**) and

*Γ*(

**x**) satisfy the following relations (see appendix A in the electronic supplementary material): 2.15a,b,c which have been used in obtaining equations (2.13

*b*,

*c*) that are simpler than those derived in Gao & Ma (2009) and Ma & Gao (2010).

Equations (2.14*a*–*c*) show that among the three functions, only the third one, *Γ*(**x**), involves the length scale parameter *L*. It then follows from equations (2.13*c*) and (2.14*b*,*c*) that depends on *L*, whereas , expressed in terms of only ** Φ**(

**x**) and

**(**

*Λ***x**) according to equations (2.13

*b*) and (2.14

*a*,

*b*), is independent of

*L*. Also, it is seen from equations (2.13

*c*) and (2.14

*b*,

*c*) that when

*L*=0 (and thus

*Γ*(

**x**)≡0), thereby giving (from equation (2.13

*a*)). That is, the newly obtained Eshelby tensor reduces to that based on classical elasticity when the strain gradient effect is not considered.

The fifth-order Eshelby-like tensor **T**, which relates the strain gradient **κ*** to the disturbed strain ** ε** in the elastic body (see equation (2.12)), can be shown to be
2.16
where
2.17
with the scalar-valued functions

**(**

*Φ***x**),

**(**

*Λ***x**) and

*Γ*(

**x**) defined in equations (2.14

*a*–

*c*). Note that use has been made of equations (2.15

*a*–

*c*) in reaching equation (2.16), which is simpler than that provided in Ma & Gao (2010). It is clear from equation (2.16) that

**T**vanishes when

*L*=0. Then, with as discussed above, equation (2.12) simply becomes when

*L*=0, which is the defining relation for the Eshelby tensor based on classical elasticity (Eshelby 1957), as expected.

Equations (2.13a–*c*) and (2.16) provide the general formulas for determining and *T*_{ijlmp} for an inclusion of arbitrary shape in terms of the scalar-valued functions ** Φ**(

**x**),

**(**

*Λ***x**) and

*Γ*(

**x**) defined in equations (2.14

*a*–

*c*). For the two simple cases of a spherical inclusion and a cylindrical inclusion, closed-form expressions have been obtained for

**(**

*Φ***x**),

**(**

*Λ***x**) and

*Γ*(

**x**) and thus for the Eshelby tensor (Gao & Ma 2009, 2010; Ma & Gao 2010). The more complex case of an ellipsoidal inclusion, for which

*Γ*(

**x**) is difficult to evaluate analytically, is examined below in this study.

## 3. Ellipsoidal inclusion

The problem of an ellipsoidal inclusion in an infinite elastic body has been studied analytically using classical elasticity (e.g. Eshelby 1957, 1959; Mura 1982; Markenscoff 1998; Ju & Sun 1999; Michelitsch *et al*. 2003; Liu 2008) or micropolar elasticity (Ma & Hu 2006, 2007). The Eshelby tensor for this problem is derived here using the general formulas based on the SSGET presented in the preceding section.

Consider an ellipsoidal inclusion of three semi-axes *a*_{1}, *a*_{2} and *a*_{3} centred at the origin of the coordinate system (*x*_{1}, *x*_{2}, *x*_{3}) in the physical space, as shown in figure 2. Then, the region *Ω* occupied by the inclusion is given by
3.1

For this ellipsoidal inclusion, *Φ*(**x**) and *Λ*(**x**) defined in equations (2.14*a*) and (2.14*b*) can be shown by integration and differentiation to satisfy (Mura 1982)
3.2a,b
3.2c,d
where
3.3a
3.3b
which are functions of *γ*, with *I*, *J*∈{1,2,3}. For **x**∈*Ω* (interior points) *γ*=0, and for **x**∉*Ω* (exterior points) *γ* is the largest positive root of the following equation:
3.4
which shows that *γ* is a function of **x**. Note that in equations (3.2*a*–*d*) and in the sequel each repeated lower-case index implies summation from 1 to 3, whereas each upper-case index takes the same value as its corresponding lower-case index but implies no summation.

For interior points with **x**∈*Ω*, it follows from equations (2.13b) and (3.2*a*,*c*) that the classical part of the Eshelby tensor is
3.5
which is the same as that provided in Li & Wang (2008). It is clear from equation (3.5) that the Eshelby tensor is uniform (i.e. independent of position **x**) inside the ellipsoidal inclusion, which is a well-known result (e.g. Markenscoff 1998) and has recently been elaborated in a broader context by Liu (2008). In fact, it can be shown that the components of given in equation (3.5) depend only on the two aspect ratios of the ellipsoidal inclusion, defined by *α*_{1}=*a*_{1}/*a*_{3} and *α*_{2}=*a*_{2}/*a*_{3}, and Poisson’s ratio of the matrix material, *v*.

For exterior points with **x**∉*Ω*, it can be shown that the use of equations (2.13b) and (3.2*b*,*d*) yields the classical part of the Eshelby tensor as (Mura 1982)
3.6
where
3.7
Clearly, equation (3.6) involves the first- and second-order derivatives of *I*_{I}(*γ*) and *I*_{IJ}(*γ*) defined in equations (3.3*a*,*b*), which are not explicit and will be replaced with more direct expressions. It can be shown from equations (3.3*a*,*b*), after some lengthy algebra, that
3.8a,b
3.8c
where
3.9a–c
3.9d–f
Using equations (3.7) and (3.8*a*–*c*) in equation (3.6) then gives
3.10
where
3.11
Note that the expression for the classical part of the Eshelby tensor outside the ellipsoidal inclusion (i.e. for exterior points with **x**∉*Ω*) given in equations (3.10) and (3.11) no longer contains derivatives of *I*_{I}(*λ*) and *I*_{IJ}(*λ*) and is therefore more convenient and more accurate to use (since differentiation tends to introduce more errors in numerical approximations). It can be readily shown that the expression of in equation (3.10) is the same as that derived earlier by Ju & Sun (1999) using a different notation. Clearly, it is seen from equations (3.6) and (3.7) that , being dependent on the position **x**, is not uniform outside the ellipsoidal inclusion, although (see equation (3.5)) is uniform inside the same inclusion.

The determination of the gradient part of the Eshelby tensor requires the evaluation of the integral defining *Γ*(**x**) in equation (2.14*c*). For the ellipsoidal inclusion described in equation (3.1) and shown in figure 2, a closed-form expression for *Γ*(**x**) can hardly be derived. However, the following results can be obtained (see appendix B in the electronic supplementary material for derivations):
3.12
for interior points **x**∈*Ω*, where
3.13a
3.13b,cand
3.14
for exterior points **x**∉*Ω*, where
3.15a
3.15band
3.15c
Clearly, equations (3.12) and (3.14) show that *Γ*(**x**)=0 when *L*=0 (i.e. when the strain gradient effect is not considered), as expected.

It should be mentioned that for interior points **x**∈*Ω*, evaluating *Γ*(**x**) defined in equation (2.14*c*) can also be reduced to the evaluation of one line integral on the interval by using an expression of the Helmholtz potential inside an ellipsoidal region derived in Michelitsch *et al*. (2003, see their eqn (3.18)), as was done in Ma & Hu (2006) for spheroidal inclusion cases with *a*_{1}=*a*_{2}.

For the special case of a spherical inclusion with *a*_{1}=*a*_{2}=*a*_{3}=*R*, thereby *s*=*R* and according to equations (3.13*a*–*c*), it can be readily shown by using equations (3.12) and (3.14), respectively, that
3.16
which are the same as those obtained in Gao & Ma (2009) by direct integration.

For the special case of a cylindrical inclusion with *a*_{1}=*a*_{2}=*a* and , it can be shown, after evaluating the integrals analytically, that equations (3.12) and (3.14) give
3.17a,b
where *I*_{n}(⋅) and *K*_{n}(⋅) (*n*=0,1) are, respectively, the modified Bessel functions of the first and second kinds of the *n*th order. Equations (3.17*a*,*b*) are the same as those derived in Ma & Gao (2010). In reaching equations (3.17*a*,*b*), use has been made of the identities:
3.17c,d

It can be shown that differentiating equation (3.12) or equation (3.14) leads to
3.18a
where
3.18b
3.18c
or
3.18d
Clearly, it is seen from equations (3.18*c*,*d*) that *f*(*P*)=0 in both the interior and exterior cases when *L*=0, as expected. Note that in reaching equation (3.18a) use has been made of the coordinate transformation:
3.19
which transforms the ellipsoidal region *Ω* defined in equation (3.1) into a unit sphere with the unit outward normal **n** on its surface |**X**|=1.

It then follows from equations (3.18a,*b*) and (3.19) that
3.20a
3.20b
3.20c
3.20d
where
3.21
with *n*_{i}=*x*_{i}/(*Pa*_{I}) according to equation (3.19).

Note that equations (3.20a–*d*) hold for both the interior and exterior cases, with *f* defined in equations (3.18c) and (3.18d), respectively. In equation (3.21), *f*′, *f*′′, *f*′′′ and *f*^{(4)} denote, respectively, the first-, second-, third- and fourth-order derivatives of *f* with respect to *P*. For the interior case with **x**∈*Ω*, *f*′, *f*′′, *f*′′′ and *f*^{(4)} can be obtained from equation (3.18c) by direct differentiation. For the exterior case with **x**∉*Ω*, the use of equation (3.18d) leads to
3.22
where *F*^{(1)} and *F*^{(2)} are defined in equations (3.15a) and (3.15b). In reaching equation (3.22) use has also been made of the result [*F*^{(1)}−*F*^{(2)}]|_{θ=α}=0, which enables the terms involving ∂*α*/∂*P* to vanish.

Using equations (3.20*a*,*c*) in equation (2.13c) then yields the gradient part of the Eshelby tensor as
3.23
where *n*_{i}=*x*_{i}/(*Pa*_{I}) from equation (3.19), and *P* is defined in equation (3.18b).

Equation (3.23) applies to both the interior and exterior cases, but the expressions for *f*(*P*) and *Λ*_{,ijlm} are different in each case. For the interior case with **x**∈*Ω*, *f*(*P*) is given in equation (3.18c) and *Λ*_{,ijlm}=0 (see equation (3.2*c*)), whereas for the exterior case with **x**∉*Ω*, *f*(*P*) is provided in equation (3.18d) and *Λ*_{,ijlm} is obtained from equation (3.2*d*), after some lengthy algebra, as
3.24
where
3.25a
3.25b
3.25c,d

From equations (3.23)–(3.25a–*d*), it is seen that the gradient part of the Eshelby tensor, , is position-dependent inside and outside the ellipsoidal inclusion, since *f*, *P*, *n*_{i} and *Λ*_{,ijlm} (for the exterior case) involved in are all functions of **x**. This differs from the classical part, , which is uniform inside the same inclusion.

Next, substituting equations (3.20*b*,*d*) into equation (2.16) gives the fifth-order Eshelby-like tensor as
3.26
where *n*_{i}=*x*_{i}/(*Pa*_{I}) from equation (3.19), *P* is defined in equation (3.18b) and *Λ*_{,ijp}, *Λ*_{,ijlmp}, *Φ*_{,ijlmp} can be obtained from equations (3.2*a*–*d*). Clearly, *T*_{ijlmp}=0 whenever *L*=0. That is, this fifth-order Eshelby-like tensor vanishes both inside and outside the ellipsoidal inclusion when the strain gradient effect is not considered.

Equation (3.26) is valid for both the interior and exterior cases, but the expressions for *f*(*P*), *Λ* and *Φ* are different in each case. For the interior case with **x**∈*Ω*, *f*(*P*) is given in equation (3.18c) and the derivatives of *Φ* and *Λ* involved are obtained from equations (3.2*a*,*c*) as
3.27
For the exterior case with **x**∉*Ω*, *f*(*P*) is provided in equation (3.18d) and the derivatives of *Φ* and *Λ* are determined from equations (3.2*b*,*d*), after some tedious derivations, as
3.28a
3.28band
3.28c
where
3.28d
and *M*_{IJK} and *N*_{IJKL} are given in equations (3.25*a*,*b*).

Considering that is position-dependent inside the ellipsoidal inclusion, its volume average over the ellipsoidal region occupied by the inclusion is examined next. This averaged Eshelby tensor is needed for predicting the effective elastic properties of a heterogeneous composite containing ellipsoidal inhomogeneities.

The volume average of a sufficiently smooth function *F*(**x**) over the ellipsoidal inclusion occupying region *Ω* is defined by, with the help of the coordinate transformation defined in equation (3.19),
3.29
where d, with *P*=|** X**| (see equation (3.18b)). It then follows from equations (2.13c), (3.2

*c*) and (3.29) that 3.30 where 3.31a,b with

*f*defined in equation (3.18c). Note that in reaching equations (3.31

*a*,

*b*) use has been made of the integral identities given in eqn (53) of Gao & Ma (2009). Using equations (3.31

*a*,

*b*) in equation (3.30) finally gives 3.32 as the average of the gradient part of the Eshelby tensor over the ellipsoidal inclusion

*Ω*. It can be readily shown that for the spherical inclusion case with

*a*

_{1}=

*a*

_{2}=

*a*

_{3}=

*R*, equation (3.32) recovers the closed-form expression of derived in Gao & Ma (2009).

Since is uniform inside the inclusion, the use of equations (3.5) and (3.29) gives . It then follows from equations (2.13a), (3.5) and (3.32) that
3.33
as the volume average of the Eshelby tensor inside the ellipsoidal inclusion based on the SSGET. In equation (3.33), *I*_{I}(0) and *I*_{IJ}(0) are constants obtainable from equations (3.3*a*,*b*), and *f*=*f*(*P*) is defined in equation (3.18c). Clearly, when *L*=0, equation (3.33) reduces to given in equation (3.5), since *f*(*P*)≡0 for any value of *P* (>0) when *L*=0 (see equation (3.18*c*)).

The volume average of *T*_{ijlmp} for **x** locating inside the ellipsoidal inclusion (i.e. **x**∈*Ω*) can be readily shown to vanish, i.e.
3.34
This is based on the fact that *T*_{ijlmp}(**x**) involved in equation (3.34), which is to be substituted from equations (3.26) and (3.27), contains the components of the unit normal vector on the unit sphere surface |**X**|=1 through *n*_{i}, *n*_{i}*n*_{j}*n*_{l} and *n*_{i}*n*_{j}*n*_{l}*n*_{m}*n*_{p}, which satisfy the following integral identities:
3.35

It then follows from equations (2.12), (3.29) and (3.34) that
3.36
where 〈*S*_{ijlm}〉_{V} is given in equation (3.33). Equation (3.36) shows that the average disturbed strain is only related to the eigenstrain ** ε*** even in the presence of the eigenstrain gradient

**κ***.

## 4. Numerical results

To quantitatively illustrate how the newly derived Eshelby tensor changes with the position and inclusion size, some numerical results are presented in this section.

Figure 3 shows the distribution of along the *x*_{3}-axis (with *x*_{1}=0=*x*_{2}) for four different values of *a*_{3}, where the values of and are, respectively, obtained from equations (3.5) and (3.23), with *f*(*P*) given in equation (3.18c) and *Λ*_{,ijlm}=0 from equation (3.2*c*). For comparison, the value of the counterpart component of the classical Eshelby tensor, which is the same as , is also displayed in figure 3, where *α*_{1}=*a*_{1}/*a*_{3} and *α*_{2}=*a*_{2}/*a*_{3} are the two aspect ratios of the ellipsoidal inclusion.

From figure 3, it is observed that *S*_{3333} varies with the position (with *x*_{1}=*x*_{2}=0, *x*=*x*_{3}) and depends on the inclusion size (*a*_{3}), unlike its classical part which, for the specified values of the aspect ratios *α*_{1} and *α*_{2}, is a constant independent of both *x*_{3} and *a*_{3} (i.e. from equation (3.5), as shown). When *a*_{3} is small (with *a*_{3}=*L*=17.6 μm here), *S*_{3333} is much smaller than , which indicates that the magnitude of is very large and the strain gradient effect is significant. However, when *a*_{3} is much greater than *L* (e.g. *R*=5*L*=88.0 μm shown here), *S*_{3333} is seen to be quite uniform and its value approaches from below, meaning that the magnitude of is very small. This indicates that for large inclusions the strain gradient effect is insignificant and may be neglected. Similar trends have been observed for other components of .

The variation of the component of the averaged Eshelby tensor inside the ellipsoidal inclusion, 〈*S*_{3333}〉_{V}, with the inclusion size *a*_{3} shown in figures 4 and 5. In figure 4, the spherical inclusion and the cylindrical inclusion cases are included as two limiting cases of the ellipsoidal inclusion problem solution with *α*_{1}=1 and , respectively. Note that 〈*S*_{3333}〉_{V} is obtained from equation (3.33) for all cases displayed in figures 4 and 5.

For the spherical inclusion (with *α*_{1}=1=*α*_{2}) and cylindrical inclusion (with and *α*_{2}=1) cases, the numerical results of 〈*S*_{3333}〉_{V} obtained using equation (3.33) and shown in figure 4 are almost identical to those determined using the closed-form formulas derived in Gao & Ma (2009) and Ma & Gao (2010). Moreover, it is clearly seen from figure 4 that the spherical inclusion case having *α*_{1}=1=*α*_{2} provides a lower bound, whereas the cylindrical inclusion case having and *α*_{2}=1 furnishes an upper bound for the ellipsoidal (spheroidal) inclusion cases with (and *α*_{2}=1), as expected. These facts verify and support the current analysis of the ellipsoidal inclusion problem.

It is observed from figure 4 that 〈*S*_{3333}〉_{V} is indeed varying with the inclusion size for all five cases considered: the smaller *a*_{3} is, the smaller 〈*S*_{3333}〉_{V} is. This size effect is seen to be significant when the inclusion is small (with *a*_{3}/*L*<10 or *a*_{3}<176 μm here). However, as the inclusion size increases, 〈*S*_{3333}〉_{V} in each case approaches from below the corresponding value of based on classical elasticity, which, for given values of the aspect ratios, is a constant independent of *a*_{3}, as discussed in §2. This comparison is further illustrated in figure 5, where it is shown that the classical Eshelby tensor component, as a constant (i.e. from equation (3.5) for *α*_{1}=3 and *α*_{2}=2), cannot explain the inclusion size effect.

## Acknowledgements

The work reported in this paper is supported by a grant from the US National Science Foundation, with Dr Clark V. Cooper as the programme manager. This support is gratefully acknowledged. The authors also thank one anonymous reviewer for the helpful comments.

## Footnotes

- Received November 28, 2009.
- Accepted February 11, 2010.

- © 2010 The Royal Society