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.
Eshelby’s solution for one particle in an infinite domain under a far perturbation stress  has been a foundation for micromechanics of composite materials , 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  has proved that the Eshelby property does not hold for Jewish star or a rectangular inclusion, and Lubarda and Markenscoff  extended the conclusion to all inclusions bounded by a non-convex or polynomial order (higher than 2) surface. Ru  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  based on Waldvogel’s work of the Newtonian potential of homogeneous polyhedra . 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 , among others [20,21].
Rodin , 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  or surface energies of the nano-inclusions .
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  and Mura  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  and Mindlin , 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 , Zhang and Sharma  modified the classic Eshelby’s tensor to model the inclusion problems of small-scale elastic phenomena and isotropic chiral solids. Ma and Gao  and Liu and Gao  extended the concept for the plane strain problem with elliptical and polygonal inclusions, respectively. For the second approach, Nakasone et al.  discretized the inclusion domain and utilize the shape function to distribute the eigenstrain field without an explicit Eshelby’s tensor. While Zhou et al.  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  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
3 Eshelby’s Tensor for a Higher-Order Eigenstrains in a Polygonal Inclusion
3.1 Linear Eshelby’s Tensor for the Polygonal Inclusion.
3.2 Quadratic Eshelby’s Tensor for the Polygonal Inclusion.
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 NF-side polygon with an equal side length is considered in this subsection. When NF 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 E0 = 70 GPa and Poisson’s ratio v0 = 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 NF 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 , the comparison at the vertex ((0, 1) and (0, −1)) is not exactly showed for NF = 4, 6, 18, 720. In Eqs. (10b) and (14b), 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 (Ax2δij), and quadratic () thermal strain to the polygonal inclusion, respectively. Figures (2)–(4) show the variation of stresses along the x2 axis. The following features of stress distributions can be observed:
When NF 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 (NF = 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 E1 = 0)/stiffer fiber (Young’s modulus E2 = 210 GPa and Poisson’s ratio ν2 = 0.3) is embedded in a homogeneous infinite matrix (Young’s modulus E0 = 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) MPa and (ii) 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 x2-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 x2-axis in the range of 3 m higher/lower to the centroid of the triangle. To avoid the singularity issue at the top vertex , two close neighbor points are used for demonstration.
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 and . 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 x2 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 , 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 x2 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].
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.
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
Appendix B: Derivatives of the Potential Components
The derivatives with respect to are similar to Eq. (B6b).