## Abstract

This paper presents the Eshelby’s tensor of a polygonal inclusion with a polynomial eigenstrain, which can provide an elastic solution to an arbitrary, convex inclusion with a continuously distributed eigenstrain by the Taylor series approximation. The Eshelby’s tensor for plane strain problem is derived from the fundamental solution of isotropic Green’s function with the Hadmard regularization, which is composed of the integrals of the derivatives of the harmonic and biharmonic potentials over the source domain. Using the Green’s theorem, they are converted to two line (contour) integrals over the polygonal cross section. This paper evaluates them by direct analytical integrals. Following Mura’s work, this paper formulates the method to derive linear, quadratic, and higher order of the Eshelby’s tensor in the polynomial form for arbitrary, convex polygonal shapes of inclusions. Numerical case studies were performed to verify the analytic results with the original Eshelby’s solution for a uniform eigenstrain in an ellipsoidal domain. It is of significance to consider higher order terms of eigenstrain for the polygon-shape inclusion problem because the eigenstrain distribution is generally non-uniform when Eshelby’s equivalent inclusion method is used. The stress disturbance due to a triangle particle in an infinite domain is demonstrated by comparison with the results of the finite element method (FEM). The present solution paves the way to accurately simulate the particle-particle, partial-boundary interactions of polygon-shape particles.

## 1 Introduction

Eshelby’s solution for one particle in an infinite domain under a far perturbation stress [1] has been a foundation for micromechanics of composite materials [2], including the Mori-Tanaka method, the self-consistent scheme, the differential scheme, etc. [3–5]. Based on the Eshelby’s equivalent inclusion method (EIM), the material mismatch between a particle and the matrix could be simulated by an eigenstrain field within the inclusion. For a single ellipsoid case, the eigenstrain is constant over the particle because the Eshelby’s tensor, derived from the volume integral of the Green’s function over an ellipsoidal source domain, is constant within the source as well. When the entire domain is finite or multiple particles exist or the shape of particle is arbitrary, the eigenstrain may not be constant any more although it is continuous. Although the classic Eshelby’s theory provided a beautiful, exact solution with the uniform eigenstrain, it is important to determine the Eshelby’s tensor for various shapes of inclusion for different non-uniform eigenstrain distribution for general cases [6,7], particularly when a particle exhibits other shapes or is close to the boundary or its neighbor particles.

This paper focuses on the Eshelby’s tensor of a polygonal inclusion with the polynomial-form eigenstrain. Several methods were developed to investigate the polygonal and polyhedral inclusions and its related properties of the Eshelby’s tensor. Chiu studied the stress field caused by an initial strain in both full and half-space for a cuboid inclusion [8–10]. Mura [11] has proved that the Eshelby property does not hold for Jewish star or a rectangular inclusion, and Lubarda and Markenscoff [12] extended the conclusion to all inclusions bounded by a non-convex or polynomial order (higher than 2) surface. Ru [13] provided the analytical solution of an arbitrary shaped inclusion and used the simple explicit expressions to obtain the internal stress in a plane or half-plane. The analytical Eshelby’s tensor for polygonal and polyhedral inclusions with uniformly distributed eigenstrain were first solved by Rodin [14] based on Waldvogel’s work of the Newtonian potential of homogeneous polyhedra [15]. Nozaki and Taya [16,17] proposed an explicit closed-form solution of the potentials, however, the method is limited to convex inclusions. Trotta et al. [7,18] derived an simplified analytical expression of the Eshelby’s tensor with the coordinates of vertices on the inclusion, which avoids the use of lengthy and complex coordinate variables. Subsequently, the special properties of the polygonal inclusions and its associated average Eshelby’s tensor were investigated by Xu and Wang [19], among others [20,21].

Rodin [14], Nozaki and Taya [16,17] and other pioneers [7,18] have developed the Eshelby’s tensor for arbitrary polygonal inclusions with a uniformly distributed eigenstrain field, which is consistent with the classic Eshelby’s problem. It provides elegant exact elastic solution for an infinite domain with an inelastic strain, namely eigenstrain, in the inclusion, which exhibits the same elasticity as the matrix. Besides the application to the physical problem for an inclusion with a uniform eigenstrain field, the most significant application of Eshelby’s tensor is to solve for inhomogeneity problems in which the particle with a different elasticity to the matrix is subjected to a uniform far field stress field. Eshelby’s equivalent inclusion method (EIM) simulates the material mismatch by an inclusion with an eigesntrain but the same elasticity as the matrix. For an ellipsoidal inhomogeneity, because the Eshelby’s tensor is a constant inside the particle, the EIM provides the exact solution with the Eshelby’s tensor. However, for a polygonal inhomogeneity, this feature does not exist due to the non-uniform Eshelby’s tensor. Therefore, the EIM with a uniform eigenstrain is not sufficient for an accurate solution. Similar problems also exist in the size effects caused by either microstructure [22] or surface energies of the nano-inclusions [23].

To improve the accuracy of solutions, two general approaches were proposed, (1) increase the order of eigenstrain field [24,25] or strain field through combination of the strain gradient theory with the classic Eshelby’s solution [6,22]; (2) discretize the inclusion domain with multiple small/basic elements [26,27]. As for the first approach, Moschovidis and Mura [24] and Mura [25] suggested the use of polynomial form eigenstrain by the Taylor expansion, which could not only address the size effect of the inclusion, but also the interactions between particles. The strain gradient theory, introduced by Toupin [28] and Mindlin [29], was developed to overcome the deficiencies of classic elasticity [30,31]. In different versions of the strain gradient theory, to accommodate the size effect, one (at least) or more characteristic length parameters are involved to describe the elastic fields. In the literature, Sharma [32], Zhang and Sharma [33] modified the classic Eshelby’s tensor to model the inclusion problems of small-scale elastic phenomena and isotropic chiral solids. Ma and Gao [6] and Liu and Gao [22] extended the concept for the plane strain problem with elliptical and polygonal inclusions, respectively. For the second approach, Nakasone et al. [26] discretized the inclusion domain and utilize the shape function to distribute the eigenstrain field without an explicit Eshelby’s tensor. While Zhou et al. [27] enforced a uniform distribution of eigenstrain in each of the element, which may arise problems of discontinuity of both eigenstrain and elastic fields in post-process. Eshelby’s tensor for polygonal inclusions with high-order terms of eigenstrain has not existed in the literature yet.

One goal of the paper is to extend the previous work of Rodin [14] to polynomial eigenstrain distributions following Mura’s work [24,25]. Using the polynomial eigenstrain could reach a tailorable accuracy for the solution of polygonal inhomogeneity problems. This is the first time in the literature to extend the EIM to a polygonal inhomogeneity although it is not the exact solution as the classic Eshelby’s problem for an ellipsoid. Section 2 provides the derivation of the explicit expression of harmonic potential through the direct integral [14,22]. Section 3 presents the scheme to derive Eshelby’s tensor with a high order polynomial form of eigenstrain and provides the explicit forms of Eshelby’s tensor for linear and quadratic cases. Section 4.1 conducts numerical case studies to verify the formulation with the existing results for circular inhomogeneity. Section 4.2 investigates the accuracy of solving the Eshelby’s EIM of a triangular cylindrical inclusion with different aspect ratio of width to height and various stiffness ratio of the particle to matrix.

## 2 Eshelby’s Tensor for a Uniform Eigenstrain on a Polygon Inclusion

*g*

_{ikl}) and strain (

*S*

_{ijkl}) [1,25] in 3-dimensional (3D) space is derived from the volume integral over the inclusion domain Ω in Eq. (1), and it could be split into biharmonic Ψ, harmonic Φ potentials.

*C*is an integral constant and it could be eliminated by partial differentiation in the process to obtain the Eshelby’s tensor.

*N*

_{F}-sided cross section of the cylindrical inclusion in the

*x*

_{1}−

*x*

_{2}plane. Let

**x**and

**x**′ denote the observing and source points, respectively. Following Rodin’s [14] work, transformed coordinates (TC) are constructed on each of the edge in Fig. 1. The base vectors of the

*f*

^{th}TC, namely the unit tangent, outward normal vectors of the

*f*

^{th}edge in the right-hand basis, are written as $\eta f0$ and $\lambda f0$, respectively. And the variables of the TC at the

*f*

^{th}edge are given as function of observing point

**x**and two vertices $vf\u2212$, $vf+$ [14,22],

*b*

_{f}is the perpendicular distance between the observing point to the edge. The vector between observing point

**x**and the source point

**x**′ is written as $x\u2212x\u2032=\u2212bf\lambda f0\u2212\eta \eta f0$, where $\eta =\eta f0(x\u2032\u2212x)$ is the position of the source point on the edge. To derive the Eshelby’s tensor, partial derivatives of the potentials are required; therefore, Green’s theorem could be utilized to simplify potentials as contour integrals. However, this work uses an alternative scheme to solve the potentials through direct integral over the source domain. Shown in Fig. 1(b), the integral of any arbitrary piece-wise smooth function $F(|x\u2212x\u2032|)$ could be performed as follows:

## 3 Eshelby’s Tensor for a Higher-Order Eigenstrains in a Polygonal Inclusion

*S*

_{ijklmpq···}is the polynomial series of the Eshelby’s tensor; $\Phi mpq\cdots =\u222b\Omega \u2212ln|x\u2212x\u2032|2xm\u2032xp\u2032xq\u2032\cdots dA(x\u2032)$; $\Psi mpq\cdots =\u222b\Omega \u2212|x\u2212x\u2032|2ln|x\u2212x\u2032|2\u22121/2xm\u2032xp\u2032xq\u2032\cdots dA(x\u2032)$. Since only the partial derivatives of the above potentials are of interest, the area integral can be transferred to the boundary by using the Green’s theorem. Shown in Fig. 1(b), the integral variable

**x**′ moves along each edge; thus, the distance between observing and source points, |

**x**−

**x**′| is written as $bf2+\eta 2$ when the source point is on the

*f*

^{th}edge. The detailed differentiation process and analytical solution to the uniform Eshelby’s tensor can be found in the previous work [7,22], which will be used in the following derivation.

### 3.1 Linear Eshelby’s Tensor for the Polygonal Inclusion.

Similarly to Eqs. (5) and (6), the components of the linear Eshelby’s tensor can be written as follows:

**x**′ exists in the area, and to apply the Green’s theorem, the distance vector at the

*f*

^{th}edge

**x**′ −

**x**is simplified as $bf\lambda f0+\eta \eta f0$. Thus, the derivatives of the two above linear potentials, which will be used in the fourth-rank linear Eshelby’s tensor, are written as follows:

*a*)

*b*)

*a*)

*b*)

*a*)

*b*)

*δ*

_{pi}Φ, which requires the analytical solution to the harmonic potential in Eq. (5). Finally, the substitution of Φ

_{p}and Ψ

_{p}potentials into Eq. (1) yields the linear Eshelby’s tensor.

### 3.2 Quadratic Eshelby’s Tensor for the Polygonal Inclusion.

*x*

_{p}+ (

*x*′

_{p}−

*x*

_{p}), to express the first-order source term

*x*

_{p}′ in the contour integrals. Similarly, the second-order term

*x*′

_{p}

*x*′

_{q}can be expressed as (

*x*

_{p}′ −

*x*

_{p})(

*x*

_{q}′ −

*x*

_{q}) −

*x*

_{p}

*x*

_{q}+

*x*

_{p}′

*x*

_{q}+

*x*

_{q}′

*x*

_{p}; thus, the potential components are written as

*a*)

*b*)

*x*

_{p}

*x*

_{q}′ and

*x*

_{q}

*x*

_{p}′ with the integrals are simplified to

*x*

_{p}(Φ/Ψ)

_{q}and

*x*

_{q}(Φ/Ψ)

_{p}, respectively. Since the uniform/linear potentials have been derived in the previous sections, the analytical solution to the last term in Eq. (13

*b*) lead to the closed-form expression of the quadratic Eshelby’s tensor. The derivatives of the two quadratic potentials, which will be used for the quadratic Eshelby’s tensor, are presented:

*a*)

*b*)

*a*)

*b*)

The same procedure can be extended to derive Eshelby’s tensor for higher-order eigenstrains.

## 4 Results and Discussion

The above analytical solution can be used to predict the elastic field caused by a continuously distributed inelastic strain or eigenstrain in the polynomial form over a polygonal inclusion, which exhibits the same stiffness of as the infinite domain. Using Eshelby’s equivalent inclusion method, it can be used to study the elastic field of a particle or inhomogeneity embedded in an infinite domain, in which the particle is simulated by an inclusion with an eigenstrain. In the following, we demonstrate some interesting results of the inclusion problem and inhomogeneity problem with the present analytical formulation.

### 4.1 Polygonal Inclusion Problem.

A *N*_{F}-side polygon with an equal side length is considered in this subsection. When *N*_{F} increases, the polygon can converge to a circular (radius *a* = 1 m) domain which recovers Eshelby’s classic solution. Consider a homogeneous infinite domain with the Young’s modulus *E*_{0} = 70 GPa and Poisson’s ratio *v*_{0} = 0.3 subjected to a local temperature change Δ*T* in the polygonal inclusion only, which can be represented by an eigenstrain or thermal strain *ε*_{ij}* = *α*Δ*Tδ*_{ij}, where *α* denotes the thermal expansion coefficient. Shown in Fig. 2(a), let *N*_{F} equal 3, 4, 6, 9, 18, 720 and centers of the polygonal inclusions locate at the origin point *O*. Fix one of the vertex at (0, 1), and the other vertices are evenly distributed on the circle with the radius *a* = 1 m. Beginning at the fourth-rank uniform Eshelby’s tensor, the strain field at the neighborhood of the vertices has logarithm singularity issues [16], the comparison at the vertex ((0, 1) and (0, −1)) is not exactly showed for *N*_{F} = 4, 6, 18, 720. In Eqs. (10*b*) and (14*b*), the linear and quadratic Eshelby’s tensors have the components of the uniform tensor; thus, the singularity issue also exists in linear, quadratic strain field too. Let *A* = *α*Δ*T* = 10^{−4} and assign uniform (*Aδ*_{ij}), linear (*Ax*_{2}*δ*_{ij}), and quadratic ($Ax22\delta ij$) thermal strain to the polygonal inclusion, respectively. Figures (2)–(4) show the variation of stresses along the *x*_{2} axis. The following features of stress distributions can be observed:

When

*N*_{F}increases, the results for polygonal inclusions approach the classic solution for the circular inclusion, which verifies the consistency and accuracy of the present formulation of polygonal inclusions.The stress variation is concentrated in the neighborhood of the inclusion with singularity on the vertices, and the far field stress approaches zero quickly.

Except the case of the triangle (

*N*_{F}= 3), the stress variations on the inner zone of the inclusions follow the similar trend of the temperature distributions, namely uniform, linear, and quadratic distributions, but exhibit bigger differences close to the boundary of the inclusion.By the combination the polynomial terms of eigenstrain distribution, the present formulation can be used to predict the elastic field caused by continuous eigenstrain on polygonal inclusions with the Taylor expansion of the eigenstrain.

### 4.2 Polygonal Inhomogeneity Problem.

In the Sec. 4.1, a verification of the linear and quadratic Eshelby’s tensors is provided by the comparison with a circular cylindrical inclusion. Although the Eshelby’s tensors are not uniform as the classic case for circular inclusions, because polynomial eigenstrains can be considered, Eshelby’s equivalent inclusion method can still be used to solve the polygonal inhomogeneity problem. This section will use triangular inhomogeneities to demonstrate the solution.

Consider one isosceles triangular cylindrical inhomogeneity embedded in an infinite domain with a uniform far field stress. The stress in the neighborhood of the inhomogeneity will be disturbed by its material mismatch with the matrix. In the following, an air void (Young’s modulus taken as *E*_{1} = 0)/stiffer fiber (Young’s modulus *E*_{2} = 210 GPa and Poisson’s ratio *ν*_{2} = 0.3) is embedded in a homogeneous infinite matrix (Young’s modulus *E*_{0} = 70 GPa and Poisson’s ratio *ν*_{0} = 0.3). The height of the triangle is 1 m, and the width of it varies from 0.5, 1, and 2 m. Consider two cases of the far-field uniform stress load: (i) $\sigma 220=\u22121$ MPa and (ii) $\sigma 110=\u22121$ MPa.

This paper provides the accuracy up to the quadratic term of eigenstrain but it can be straightforwardly extended to higher-order terms by adding higher order stress equivalent equations in the above. If the Eshelby’s tensor is uniform over the inhomogeneity, which is the case of circles, uniform eigenstrain and stress on the inhomogeneity can be obtained, and the linear or higher order eigenstrain terms will be zero. However, unlike the circular case, the stress field in the triangular inhomogeneity is not uniform. Using only the constant term will cause the loss of accuracy. When higher orders of eigenstrain are applied, it could better describe the actual solution. The finite element method (FEM) is used to evaluate the accuracy of the stress variation along the *x*_{2}-axis under the aforementioned two cases of load with different shapes of isosceles triangular inhomogeneities. To address the singularity, very fine mesh has been used in the neighborhood of the triangle’s vertices. Figure 5 shows one example of width *w* = 0.5 m. The total domain dimension is 20 times as the triangle to mimic an infinite domain.

Notice that in the following case studies, the EIM is applied at the centroid of the triangular inhomogeneity, which means the local Taylor expansion is based on the centroid. The accuracy will reduce with the distance from the centroid in general. Although other points on the inhomogeneity can be used to conduct the stress equivalent conditions, overall the centroid can be the most representative, so that this paper focuses on this case. To illustrate the stress distribution, the observing points are evenly distributed along the *x*_{2}-axis in the range of 3 m higher/lower to the centroid of the triangle. To avoid the singularity issue at the top vertex $(0,x2p)$, two close neighbor points $(0,x2p\xb110\u22123)$ are used for demonstration.

Figures 6 and 7 show the stress variations along *x*_{2} under two far-field stresses, respectively. The following features of stress distributions can be observed:

The stress variation is concentrated in the neighborhood of the inclusion, and the far field stress approaches the applied load quickly.

Stress singularity occurs at the top vertex but not at the mid-point of the bottom edge.

Overall, the uniform, linear, and quadratic cases asymptotically approach the finite element results although the trend is not 100 $%$ consistent. It confirms the present formulation can provide tailorable accuracy by using higher-order terms of eigenstrain.

When a horizontal compressive stress is applied in Fig. 6, both

*σ*_{11}and*σ*_{22}are in the compressive range; whereas the vertical compressive stress produces tensile horizontal stress at the bottom of the triangular void in Fig. 7, which may lead to open-mode cracking by compression.Although the stress vector around the air void surface is relaxed, the internal stress in the neighborhood of the void can be higher with the disturbed stress flow. However, no stress concentration factor can be defined due to the singularity.

When the inhomogeneity becomes stiffer, the stress on the inhomogeneity can be higher. Figure 8 shows the comparison of the stress distribution for an air void and a stiffer inhomogeneity (fiber). As Figs. 6 and 7 confirm the accuracy of the uniform, linear, and quadratic cases, and the quadratic order performs the best among them, in the following, only the quadratic case is shown.The following features can be observed:

Unlike the void case with zero stress, the stress on the stiffer inhomogeneity exhibits the same sign as the far-field load, i.e., compressive stress caused by compressive load for both $\sigma 110$ and $\sigma ss0$. Moreover, the stress value on the inhomogeneity is higher than the far-field load.

Although stress singularity still exists at the top vertex, the EIM uses two points to approximate the stress field and the results are close to that of FEM. The agreement between EIM and FEM is better for the case of a stiffer inhomogeneity because the deformation of the inhomogeneity can be approximated better with a quadratic eignestrain.

In both cases of the uniform horizontal and vertical far-field load conditions, the stress distribution along

*x*_{2}reverses the trend of the void case, which could be well explained with Eq. (16). Changing the material properties of the inhomogeneity affect the sign of $Cijkl1\u2212Cijkl0$, which produces series of eigenstrains with opposite signs.

When the aspect ratio (AR) (AR = width/height) of triangle changes, the stress distribution could also be significantly different. Figures 9(a) and 9(b) show the variation of the stress comparison along the *x*_{2} with different width 0.5, 1, and 2 m of a triangular void; however, Figs. 9(c) and 9(d) show the case for stiffer inhomogeneity (fiber). The following features can be observed:

Under a horizontal compressive far-field stress,

*σ*_{11}at the mid of the bottom edge increases with a shorter width, which could be caused by the higher influence by two bottom vertices.For the air void case, when a vertical compressive far-field stress is applied, the stress

*σ*_{22}around the top vertex increases significantly with the aspect ratio. It even causes a tensile stress for AR = 2.When AR is far from 1, the difference between the FEM and EIM results becomes larger. Because the selection of equivalent stress point is the centroid of the triangle, the irregularity of the triangle may lead to the loss of the accuracy of eigenstrain approximated by the Taylor expansion. Higher order terms of eigenstrain may provide better results.

The stiffer inhomogeneity cases still exhibit reverse trend of the void cases.

Overall, the EIM with quadratic eigenstrain on a single inhomogeneity provides very accurate results in comparison with the FEM that use a large number of elements, which will lead to a high-efficient, high-fidelity numerical method for simulation of multiple polygonal particles in a composite. This fundamental work also serve as a foundation of our software package, the inclusion-based boundary element method (iBEM) for virtual experiments with many arbitrary shaped particles in a finite domain with various loading conditions [35–37].

## 5 Conclusions

The integral scheme of the linear, quadratic, and higher order terms of eigenstrain for the isotropic elastic arbitrary-shaped polygonal inclusion has been presented. Rather than using the Green’s theorem converting the area integrals into the contour integrals, this paper uses the direct integral to derive the finite part of the biharmonic and harmonic potentials. The analytical closed-form solutions to the linear/quadratic Eshelby’s tensors are obtained by assembling the potential components for different orders of eigenstrains. Numerical verification of the circular inclusion was performed, and the present solution approaches the classic solution for the circular inclusions. The formulation has been used to study the stress concentration caused by an isosceles triangular hole, which approaches the FEM results when higher-order terms of eigenstrain are used. When a compressive load is applied along the height direction, tensile stress may be induced along the bottom; while a compressive load is applied along the width direction, the normal stresses along the symmetric axis are compressive. Parametric studies show the stress disturbance caused by the Young’s modulus, aspect ratio of the inhomogeneity. It is observed that when quadratic terms of eigenstrain were used, the EIM could provide fairly accurate solution for most cases of polygonal inhomogeneity.

## Acknowledgment

This work is sponsored by the National Science Foundation IIP#1738802 and CMMI#1762891, whose support is gratefully acknowledged.

## Conflict of Interest

There are no conflicts of interest.

## Data Availability Statement

The authors attest that all data for this study are included in the paper.

### Appendix A: Differentiation Rule of the Functions $P(bf,lf\u2212,lf+)$

### Appendix B: Derivatives of the Potential Components

#### Derivatives of $\Psi 0f$

*a*)

*b*)

*c*)

*d*)

*e*)

*f*)

#### Derivatives of $\Phi 0f$

*a*)

*b*)

*c*)

#### Derivatives of $\Psi If$

*a*)

*b*)

*c*)

*d*)

*e*)

*f*)

*g*)

#### Derivatives of $\Phi If$

*a*)

*b*)

*c*)

#### Derivatives of $\Psi IIf$

*a*)

*b*)

*c*)

*d*)

#### Derivatives of $\Psi IIf$

*a*)

*b*)

The derivatives with respect to $lf\u2212$ are similar to Eq. (B6*b*).