## 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  has been a foundation for micromechanics of composite materials , including the Mori-Tanaka method, the self-consistent scheme, the differential scheme, etc. . 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 . 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

Consider an infinite long arbitrary cylindrical inclusion Ω embedded in the infinite homogeneous isotropic elastic domain $D$, then the Eshelby-type problem could be simplified as a plane strain one. The classic Eshelby’s tensor for displacement (gikl) and strain (Sijkl) [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.
$gikl=18π(1−ν)[Ψ,ikl−2νδklΦ,i−2(1−ν)×(δilΦ,k+δikΦ,l)]Sijkl=18π(1−ν)[Ψ,klij−2νδklΦ,ij−(1−ν)×(δilΦ,jk+δ,ikΦ,jl+δjlΦ,ik+δjkΦ,il)]$
(1)
where the integrals of the potential functions are written as $Ψ=∫Ω|x−x′|dx′$ and $Φ=∫Ω1|x−x′|dx′$. Explicit form solution of the integrals over ellipsoidal inclusions can be found in Mura’s work . Based on the assumption for the plane strain problem, the volume integrals, Ψ, Φ, could be further reduced to the area integrals through the Hadmard regularization [6,22,34], thus the two potentials can be rewritten as:
$Ψ=∫ΩC−|x−x′|2ln|x−x′|2−12+CdA(x′)Φ=∫ΩC−ln|x−x′|2dA(x′)$
(2)
where C is an integral constant and it could be eliminated by partial differentiation in the process to obtain the Eshelby’s tensor.
In the following, the integrals over polygonal inclusions will be derived. Without the loss of any generality, consider a NF-sided cross section of the cylindrical inclusion in the x1x2 plane. Let x and x′ denote the observing and source points, respectively. Following Rodin’s  work, transformed coordinates (TC) are constructed on each of the edge in Fig. 1. The base vectors of the fth TC, namely the unit tangent, outward normal vectors of the fth edge in the right-hand basis, are written as $ηf0$ and $λf0$, respectively. And the variables of the TC at the fth edge are given as function of observing point x and two vertices $vf−$, $vf+$ [14,22],
$bf=((vf+)k−xk)(λf0)k,lf±=((vf±)k−xk)(ηf0)k$
(3)
where bf 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−x′=−bfλf0−ηηf0$, where $η=ηf0(x′−x)$ 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−x′|)$ could be performed as follows:
$I=∫ΩcF(|x−x′|)dA(x′)=∑f=1NF∫tan−1(lf−/bf)tan−1(lf+/bf)∫0bf1+tan2[θ]F(ρ)ρdρdθ$
(4)
Let $FΦ=ln|x−x′|2$, the analytical solution of the integral of the harmonic potential is obtained
$Φ=∑f=1NFbf{−32(lf−−lf+)+bf(tan−1(lf−bf)−tan−1(lf+bf))+12lf−ln(bf2+(lf−)2)−12lf+ln(bf2+(lf+)2)}$
(5)
Let $FΨ=−|x−x′|2ln|x−x′|2−1/2$, the integral of the finite part of the biharmonic potential is also obtained
$Ψ=∑f=1NF1144bf{−51bf2lf−−13(lf−)3+51bf2lf++13(lf+)3+24bf3tan−1(lf−bf)−24bf3tan−1(lf+bf)+6lf−(2bf2+(lf−)2)ln(bf2+(lf−)2)−6lf+(2bf2+(lf+)2)ln(bf2+(lf+)2)}$
(6)

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

Following Moschovdis and Mura’s  and Mura’s work , the continuous eigenstrain field ($ϵij*$) could be expanded in the polynomial form of a Taylor series: $ϵij*(x′)=ϵij0*+xp′ϵijp1*+xp′xq′ϵijpq2*+…$. The elastic field caused by the eigenstrain can be calculated with the superposition of the contribution of each order of the eigenstrains. Accordingly, the Eshelby’s tensor can be calculated for each order of the polynomial series. For example, similarly to Eq. (2) for constant eigenstrain, the eigenstrain term of $xm′xp′xq′ϵijmpq2*$ will lead to the following integrals:
$giklmpq⋯=18π(1−ν){Ψmpq⋯,ikl−2νδklΦmpq⋯,i−2(1−ν)×(δilΦmpq⋯,k+δikΦmpq⋯,l)}Sijklmpq⋯=18π(1−ν){Ψmpq⋯,klij−2νδklΦmpq⋯,ij−(1−ν)×(δilΦmpq⋯,jk+δikΦmpq⋯,jl+δjlΦmpq⋯,ik+δjkΦmpq⋯,il)}$
(7)
where the Sijklmpq··· is the polynomial series of the Eshelby’s tensor; $Φmpq⋯=∫Ω−ln|x−x′|2xm′xp′xq′⋯dA(x′)$; $Ψmpq⋯=∫Ω−|x−x′|2ln|x−x′|2−1/2xm′xp′xq′⋯dA(x′)$. 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, |xx′| is written as $bf2+η2$ when the source point is on the fth 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:

$Ψp=∫Ω(−|x−x′|2ln|x−x′|2−12)[xp+xp′−xp]dA(x′)=xpΨ+∫Ω(−|x−x′|2ln|x−x′|2−12)[xp′−xp]dA(x′)Φp=∫Ω(−ln|x−x′|2)[xp+xp′−xp]dA(x′)=xpΦ+∫Ω(−ln|x−x′|2)[xp′−xp]dA(x′)$
(8)
where Ψ, Φ are defined in Eq. (5). The integral limits are determined as $lf+=ηf0⋅(vf+−x)$ and $lf−=ηf0⋅(vf−−x)$. Following the chain rule of derivatives, the piece-wise smooth function $P(bf,lf+,lf−)$ could be differentiated as follows:
$P,i=−(λf0)i∂P∂bf−(ηf0)i[∂P∂lf++∂P∂lf−]$
(9)
and the higher order differentiation formulas are presented in Appendix  A similarly. Notice that, in Fig. 1(b), the source point x′ exists in the area, and to apply the Green’s theorem, the distance vector at the fth edge x′ − x is simplified as $bfλf0+ηη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:
$Ψp,ijkl=δpiΨ,jkl+δpjΨ,ikl+δpkΨ,ijl+δplΨ,ijk+∑f=1NF(λf0)i(λf0)p[−(λf0)j(Ψ0f),kl−(λf0)k(Ψ0f),jl−(λf0)l(Ψ0f),jk]+(λf0)i(ηf0)p(ΨIf),jkl$
(10a)
$Φp,ij=δpiΦ,j+δpjΦ,i+xpΦ,ij+∑f=1NF(λf0)i(λf0)p[−(λf0)jΦ0f+bf(Φ0f),j]+(λf0)i(ηf0)p(ΦIf),j$
(10b)
where $Ψ0f,Ψ1f$, $Φ0f$, and $Φ1f$ denote the contour integrals in the following:
$Ψ0f=∫lf−lf+(bf2+η2)ln(bf2+η2)−12dη=−7bf2(lf+−lf−)6−5[(lf+)3−(lf−)3]18+2bf33[tan−1(lf+bf)−tan−1(lf−bf))+12ln(bf2+(lf+)2)[bf2lf++(lf+)33]−12ln(bf2+(lf−)2)[bf2lf−+(lf−)33]$
(11a)
$ΨIf=∫lf−lf+η(bf2+η2)ln(bf2+η2)−12dη=116[2bf4+6bf2(lf−)2+3(lf−)4−2(bf2+(lf−)2)2ln(bf2+(lf−)2)]−116[2bf4+6bf2(lf+)2+3(lf+)4−2(bf2+(lf+)2)2ln(bf2+(lf+)2)]$
(11b)
and
$Φ0f=∫lf−lf+ln(bf2+η2)dη=2bf[tan−1(lf+bf)−tan−1(lf−bf)]+lf+ln(bf2+(lf+)2)−lf−ln(bf2+(lf−)2)$
(12a)
$Φ1f=∫lf−lf+ηln(bf2+η2)dη=12((lf−)2−(bf2+(lf−)2)ln(bf2+(lf−)2))−12((lf+)2−(bf2+(lf+)2)ln(bf2+(lf+)2))$
(12b)
The derivatives of the potentials Ψ and Φ (Appendix  B) could be derived through Eq. (5), i.e., $Ψ,ijkl=∑f=1NF(λf0)i(Ψ0f),jkl$ by the Green’s theorem. Notice that the linear Eshelby’s tensor for the displacement field contains δ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.

As previously shown in Sec. 3.1, the derivation of linear Eshelby’s tensor uses the observing points and the distance vector, xp + (xpxp), to express the first-order source term xp′ in the contour integrals. Similarly, the second-order term xpxq can be expressed as (xp′ − xp)(xq′ − xq) − xpxq + xpxq + xqxp; thus, the potential components are written as
$Ψpq=∫Ωc(−|x−x′|2ln|x−x′|2−12)×[(xp′−xp)(xq′−xq)−xpxq+xpxq′+xqxp′]dA(x′)=−xpxqΨ+xpΨq+xqΨp+∫Ωc(−|x−x′|2ln|x−x′|2−12)(xp′−xp)(xq′−xq)dA(x′)$
(13a)
$Φpq=∫Ωc(−ln|x−x′|2)[(xp′−xp)(xq′−xq)−xpxq+xpxq′+xqxp′]dA(x′)=−xpxqΦ+xpΦq+xqΦp+∫Ωc(−ln|x−x′|2)(xp′−xp)(xq′−xq)×dA(x′)$
(13b)
where the terms xpxq′ and xqxp′ with the integrals are simplified to xp(Φ/Ψ)q and xq(Φ/Ψ)p, respectively. Since the uniform/linear potentials have been derived in the previous sections, the analytical solution to the last term in Eq. (13b) 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:
$Ψpq,ijkl=[δpiΨq,jkl+δpjΨq,ikl+δpkΨq,ijl+δplΨq,ijk+xpΨq,ijkl]+[δqiΨp,jkl+δqjΨp,ikl+δqkΨp,ijl+δqlΨp,ijk+xqΨp,ijkl]−[δpiδqj+δqiδpj]Ψ,kl−[δpiδqk+δqiδpk]Ψ,jl−[δpiδql+δplδqi]Ψ,jk−[δpixq+xpδqi]Ψ,jkl−[δpjδqk+δpkδqj]Ψ,il−[δpjδql+δplδqj]Ψ,ik−[δpjxq+xpδqj]Ψ,ikl−[δpkδql+δplδqk]Ψ,ij−[δpkxq+xpδqk]Ψ,ij−[δplxq+xpδql]Ψ,ijk−xpxqΨ,ijkl$
$+∑f=1NF(λf0)i(λf0)p(λf0)q[2(λf0)k(λf0)jΨ,lf+2(λf0)l(λf0)jΨ,kf−2bf(λf0)jΨ,klf+2(λf0)l(λf0)kΨ,jf−2bf(λf0)kΨ,jlf−2bf(λf0)lΨ,jkf+bf2Ψ,jklf]+(λf0)i((λf0)p(ηf0)q+(λf0)q(ηf0)p)[−(λf0)j(ΨIf),kl−(λf0)k(ΨIf),jl−(λf0)l(ΨIf),jk+bf(ΨIf),jkl]+(λf0)i(ηf0)p(ηf0)q(ΨIIf)jkl$
(14a)
$Φpq,ij=[δpiΦq,j+δpjΦq,i+xpΦq,ij]+[δqiΦp,j+δqjΦp,i+xqΦp,ij]−δpiδqjΦ−[δpixq+δqixp]Φ,j−[δpjxq+xpδqj]Φ,i−xpxqΦ,ij+∑f=1NF(λf0)i(λf0)p(λf0)q[−2bf(λf0)jΦf+bf2Φ,jf]+(λf0)i[(λf0)p(ηf0)q+(λf0)q(ηf0)p][−(λf0)jΦf+bfΦ,jf]+(λf0)i(ηf0)p(ηf0)q(ΦIIf),j$
(14b)
where $ΨIIf$ and $ΦIIf$ denote the contour integrals as follows:
$ΨIIf=−12[−4bf4lf−15+19bf2(lf−)345+7(lf−)525+415bf5tan−1(lf−bf)−13bf2(lf−)3ln(bf2+(lf−)2)−15(lf−)5ln(bf2+(lf−)2)]+12[−4bf4lf+15+19bf2(lf+)345+7(lf+)525+415bf5tan−1(lf+bf)−13bf2(lf+)3ln(bf2+(lf+)2)−15(lf+)5ln(bf2+(lf+)2)]$
(15a)
$ΦIIf=19[−6bf2lf−+2(lf−)3+6bf2(lf+)2−2(lf+)3+6bf3tan−1(lf+bf)−6bf3tan−1(lf−bf)−3(lf−)3ln(bf2+(lf−)2)+3(lf+)3ln(bf2+(lf+)2)]$
(15b)

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* = αΔ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 (ij), linear (Ax2δij), and quadratic ($Ax22δij$) 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:

1. 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.

2. The stress variation is concentrated in the neighborhood of the inclusion with singularity on the vertices, and the far field stress approaches zero quickly.

3. 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.

Fig. 2

### 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) $σ220=−1$ MPa and (ii) $σ110=−1$ MPa.

Following Mura and Moschovidis , the origin will be set at the centroid of the triangle, and the stress caused by a polynomial eigenstrain can be written in the form of the Taylor expansion over the inhomogeneity, so that the stress equivalent condition can be used to derive the eigenstrain terms as follows:
$(Cijkl1−Cijkl0)(Sklabϵab0*+Sklabpϵabp1*+Sklabpqϵabpq2*)+Cijkl0ϵkl0*=0(Cijkl1−Cijkl0)(Sklabrϵab0*+Sklabprϵabp1*+Sklabpqrϵabpq2*)+Cijkl0ϵklr1*=012!(Cijkl1−Cijkl0)(Sklabrwϵab0*+Sklabprwϵabp1*+Sklabpqrwϵabpq2*)+Cijkl0ϵklrw2*=0$
(16)

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 $(0,x2p)$, two close neighbor points $(0,x2p±10−3)$ are used for demonstration.

Figures 6 and 7 show the stress variations along x2 under two far-field stresses, respectively. The following features of stress distributions can be observed:

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

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

3. 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.

4. 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.

5. 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:

1. 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 $σ110$ and $σss0$. Moreover, the stress value on the inhomogeneity is higher than the far-field load.

2. 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.

3. 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 $Cijkl1−Cijkl0$, 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:

1. 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.

2. 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.

3. 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.

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

Fig. 9 Fig. 9 Close modal

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 .

## 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−,lf+)$

Based on the Rodin and Liu’s work [14,22], the differentiation chain rule utilize the geometry construction of the variable $bf,lf−,lf+$ as shown in Eq. (9). Here, the second/third differentiation procedures are attached,
$P,ij=∂2P∂bf2(λf0)i(λf0)j+(∂2P∂bf∂lf++∂2P∂bf∂lf−)[(λf0)i(ηf0)j+(λf0)j(ηf0)i]+[∂2P∂(lf+)2+∂2P∂(lf−)2](ηf0)i(ηf0)j$
(A1)
$P,ijk=−∂3P∂bf3(λf0)i(λf0)j(λf0)k−(∂3P∂bf2∂lf++∂3P∂bf2∂lf−)[(λf0)i(λf0)j(ηf0)k+(λf0)i(ηf0)j(λf0)k+(ηf0)i(λf0)j(λf0)k]−[∂3P∂bf∂(lf+)2+∂3P∂bf∂(lf−)2][(λf0)i(ηf0)j(ηf0)k+(ηf0)i(λf0)j(ηf0)k+(ηf0)i(ηf0)j(λf0)k]−[∂3P∂(lf+)3+∂3P∂(lf−)3](ηf0)i(ηf0)j(ηf0)k$
(A2)
Since the upper and lower limits of the integral are $lf+$ and $lf−$, respectively. Hence, the derivatives with respect to $lf−$ and $lf+$ ($∂n+m/∂(lf−)n∂(lf+)m$) are discarded.

### Appendix B: Derivatives of the Potential Components

#### Derivatives of $Ψ0f$

$∂3Ψf∂(lf+)3=2bf4(lf+)2+4bf2(lf+)4+2(lf+)6+(bf6+3bf4(lf+)2+3bf2(lf+)4+(lf+)6)ln[bf2+(lf+)2](bf2+(lf+)2)3$
(B1a)
$∂3Ψf∂(lf−)3=−2bf4(lf−)2+4bf2(lf−)4+2(lf−)6+(bf6+3bf4(lf−)2+3bf2(lf−)4+(lf−)6)(bf2+(lf−)2)3ln[bf2+(lf−)2]$
(B1b)
$∂3Ψf∂(lf+)2∂bf=2bf5lf++4bf3(lf+)3+2bf(lf+)5(bf2+(lf+)2)3$
(B1c)
$∂3Ψf∂(lf−)2∂bf=−2bf5lf−+4bf3(lf−)3+2bf(lf−)5(bf2+(lf−)2)3$
(B1d)
$∂3Ψf∂lf+∂bf2=2bf6+4bf4(lf+)2+2bf2(lf+)4+(bf6+3bf4(lf+)2+3bf2(lf+)4+(lf+)6)ln[bf2+(lf+)2](bf2+(lf+)2)3$
(B1e)
$∂3Ψf∂lf−∂bf2=−2bf6+4bf4(lf−)2+2bf2(lf−)4+(bf6+3bf4(lf−)2+3bf2(lf−)4+(lf−)6)ln[bf2+(lf−)2](bf2+(lf−)2)3$
(B1f)

#### Derivatives of $Φ0f$

$∂Φf∂bf=2(tan−1[lf+bf]−tan−1[lf−bf])$
(B2a)
$∂Φf∂lf+=ln[bf2+(lf+)2]$
(B2b)
$∂Φf∂lf−=−ln[bf2+(lf−)2]$
(B2c)

#### Derivatives of $ΨIf$

$∂3ΨIf∂bf3=2bf3[(lf−)2−(lf+)2](bf2+(lf−)2)(bf2+(lf+)2)−3bfln[bf2+(lf−)2]+3bfln[bf2+(lf+)2]$
(B3a)
$∂3ΨIf∂bf2∂lf+=lf+[2bf2bf2+(lf+)2+ln[bf2+(lf+)2]]$
(B3b)
$∂3ΨIf∂bf2∂lf−=−lf−[2bf2bf2+(lf−)2+ln[bf2+(lf−)2]]$
(B3c)
$∂3ΨIf∂bf∂(lf+)2=bf[2(lf+)2bf2+(lf+)2+ln[bf2+(lf+)2]]$
(B3d)
$∂3ΨIf∂bf∂(lf−)2=−bf[2(lf−)2bf2+(lf−)2+ln[bf2+(lf−)2]]$
(B3e)
$∂3ΨIf∂(lf+)3=2(lf+)3bf2+(lf+)2+3lf+ln[bf2+(lf+)2]$
(B3f)
$∂3ΨIf∂(lf−)3=−2(lf−)3bf2+(lf−)2+3lf+ln[bf2+(lf−)2]$
(B3g)

#### Derivatives of $ΦIf$

$∂ΦPf∂bf=bf[−ln(bf2+(lf−)2)+ln(bf2+(lf+)2)]$
(B4a)
$∂ΦPf∂lf+=lf+ln[bf2+(lf+)2]$
(B4b)
$∂ΦPf∂lf−=−lf−ln[bf2+(lf−)2]$
(B4c)

#### Derivatives of $ΨIIf$

$∂3ΨIIf∂bf3=2bf(bf2+(lf−)2)(bf2+(lf+)2)[−4bf4lf−−3bf2(lf−)3+4bf4lf+⋅4bf2(lf−)2lf+−4bf2(lf+)2lf−−3(lf−)3(lf+)2+3bf2(lf+)3+3(lf−)2(lf+)3+4bf(bf2+(lf−)2)(bf2+(lf+)2)(tan−1[lf−bf]−tan−1[lf+bf])]$
(B5a)
$∂3ΨIIf∂bf∂(lf+)2=2bflf+((lf+)2bf2+(lf+)2+ln[bf2+(lf+)2])$
(B5b)
$∂3ΨIIf∂bf2∂lf+=(lf+)2(2bf2bf2+(lf+)2+ln[bf2+(lf+)2])$
(B5c)
$∂3ΨIIf∂(lf+)3=−bf4−2bf2(lf+)2+(lf+)4+(bf4+7bf2(lf+)2+6(lf+)4)ln[bf2+(lf+)2]bf2+(lf+)2$
(B5d)

The derivatives with respect to $lf−$ are similar to Eqs. (B5b), (B5c), and (B5d).

#### Derivatives of $ΨIIf$

$∂ΦIIf∂bf=2bf[−lf−+lf++bf(tan−1[lf−bf]−tan−1[lf+bf])]$
(B6a)
$∂ΦIIf∂lf+=(lf+)2ln[bf2+(lf+)2]$
(B6b)

The derivatives with respect to $lf−$ are similar to Eq. (B6b).

## References

1.
Eshelby
,
J. D.
,
1957
, “
The Determination of the Elastic Field of an Ellipsoidal Inclusion, and Related Problems
,”
Proc. R. Soc. London., A.
,
241
(
1226
), pp.
376
396
. 10.1098/rspa.1957.0133
2.
Zhou
,
K.
,
Hoh
,
H. J.
,
Wang
,
X.
,
Keer
,
L. M.
,
Pang
,
J. H.
,
Song
,
B.
, and
Wang
,
Q. J.
,
2013
, “
A Review of Recent Works on Inclusions
,”
Mech. Mater.
,
60
, pp.
144
158
. 10.1016/j.mechmat.2013.01.005
3.
Kanit
,
T.
,
Forest
,
S.
,
Galliet
,
I.
,
Mounoury
,
V.
, and
Jeulin
,
D.
,
2003
, “
Determination of the Size of the Representative Volume Element for Random Composites: Statistical and Numerical Approach
,”
Int. J. Solids. Struct.
,
40
(
13–14
), pp.
3647
3679
. 10.1016/S0020-7683(03)00143-4
4.
Yin
,
H. M.
,
Buttlar
,
W. G.
,
Paulino
,
G. H.
, and
Benedetto
,
H. D.
,
2008
, “
Assessment of Existing Micro-Mechanical Models for Asphalt Mastics Considering Viscoelastic Effects
,”
,
9
(
1
), pp.
31
57
. 10.1080/14680629.2008.9690106
5.
Jang
,
S.-H.
,
Hochstein
,
D. P.
,
Kawashima
,
S.
, and
Yin
,
H.
,
2017
, “
Experiments and Micromechanical Modeling of Electrical Conductivity of Carbon Nanotube/cement Composites With Moisture
,”
4Cement Concrete Composites
,
77
, pp.
49
59
. 10.1016/j.cemconcomp.2016.12.003
6.
Ma
,
H. M.
, and
Gao
,
X. L.
,
2009
, “
Eshelby’s Tensors for Plane Strain and Cylindrical Inclusions Based on a Simplified Strain Gradient Elasticity Theory
,”
Acta Mech.
,
211
(
1–2
), pp.
115
129
. 10.1007/s00707-009-0221-0
7.
Trotta
,
S.
,
Marmo
,
F.
, and
Rosati
,
L.
,
2017
, “
Evaluation of the Eshelby Tensor for Polygonal Inclusions
,”
Compos. Part B: Eng.
,
115
, pp.
170
181
. 10.1016/j.compositesb.2016.10.018
8.
Chiu
,
Y. P.
,
1977
, “
On the Stress Field Due to Initial Strains in a Cuboid Surrounded by An Infinite Elastic Space
,”
ASME J. Appl. Mech.
,
44
(
4
), pp.
587
590
. 10.1115/1.3424140
9.
Chiu
,
Y. P.
,
1978
, “
On the Stress Field and Surface Deformation in a Half Space With a Cuboidal Zone in which Initial Strains are Uniform
,”
ASME J. Appl. Mech.
,
45
(
2
), pp.
302
306
. 10.1115/1.3424292
10.
Chiu
,
Y. P.
,
1980
, “
On the Internal Stresses in a Half Plane and a Layer Containing Localized Inelastic Strains Or Inclusions
,”
ASME J. Appl. Mech.
,
47
(
2
), pp.
313
318
. 10.1115/1.3153661
11.
Mura
,
T.
,
1997
, “
The Determination of the Elastic Field of a Polygonal Star Shaped Inclusion
,”
Mech. Res. Communicat.
,
24
(
5
), pp.
473
482
. 10.1016/S0093-6413(97)00052-9
12.
Lubarda
,
V.
, and
Markenscoff
,
X.
,
1998
, “
On the Absence of Eshelby Property for Non-Ellipsoidal Inclusions
,”
Int. J. Solids. Struct.
,
35
(
25
), pp.
3405
3411
. 10.1016/S0020-7683(98)00025-0
13.
Ru
,
C. Q.
,
1999
, “
Analytic Solution for Eshelby’s Problem of An Inclusion of Arbitrary Shape in a Plane or Half-Plane
,”
ASME J. Appl. Mech.
,
66
(
2
), pp.
315
523
. 10.1115/1.2791051
14.
Rodin
,
G. J.
,
1996
, “
Eshelby’s Inclusion Problem for Polygons and Polyhedra
,”
J. Mech. Phys. Solids.
,
44
(
12
), pp.
1977
1995
. 10.1016/S0022-5096(96)00066-X
15.
Waldvogel
,
J.
,
1979
, “
The Newtonian Potential of Homogeneous Polyhedra
,”
Z. Angewandte Math. Phys. ZAMP
,
30
(
2
), pp.
388
398
. 10.1007/BF01601950
16.
Nozaki
,
H.
, and
Taya
,
M.
,
1997
, “
Elastic Fields in a Polygon-Shaped Inclusion With Uniform Eigenstrains
,”
ASME J. Appl. Mech.
,
64
(
3
), pp.
495
502
. 10.1115/1.2788920
17.
Nozaki
,
H.
, and
Taya
,
M.
,
2000
, “
Elastic Fields in a Polyhedral Inclusion With Uniform Eigenstrains and Related Problems
,”
ASME J. Appl. Mech.
,
68
(
3
), pp.
441
452
. 10.1115/1.1362670
18.
Trotta
,
S.
,
Zuccaro
,
G.
,
Sessa
,
S.
,
Marmo
,
F.
, and
Rosati
,
L.
,
2018
, “
On the Evaluation of the Eshelby Tensor for Polyhedral Inclusions of Arbitrary Shape
,”
Composit. Part B: Eng.
,
144
, pp.
267
281
. 10.1016/j.compositesb.2018.01.012
19.
Xu
,
B.
, and
Wang
,
M.
,
2005
, “
Special Properties of Eshelby Tensor for a Regular Polygonal Inclusion
,”
Acta. Mech. Sin.
,
21
(
3
), pp.
267
271
. 10.1007/s10409-005-0034-x
20.
Zou
,
W.
,
He
,
Q.
,
Huang
,
M.
, and
Zheng
,
Q.
,
2010
, “
Eshelby’s Problem of Non-Elliptical Inclusions
,”
J. Mech. Phys. Solids.
,
58
(
3
), pp.
346
372
. 10.1016/j.jmps.2009.11.008
21.
Kawashita
,
M.
, and
Nozaki
,
H.
,
2001
, “
Eshelby Tensor of a Polygonal Inclusion and Its Special Properties
,”
J. Elasticity
,
64
(
1
), pp.
71
84
. 10.1023/A:1014880629679
22.
Liu
,
M.
, and
Gao
,
X.-L.
,
2013
, “
Strain Gradient Solution for the Eshelby-Type Polygonal Inclusion Problem
,”
Int. J. Solids. Struct.
,
50
(
2
), pp.
328
338
. 10.1016/j.ijsolstr.2012.09.010
23.
Sharma
,
P.
, and
Ganti
,
S.
,
2004
, “
Size-dependent Eshelby’s Tensor for Embedded Nano-Inclusions Incorporating Surface/interface Energies
,”
ASME J. Appl. Mech.
,
71
(
5
), pp.
663
671
. 10.1115/1.1781177
24.
Moschovidis
,
Z. A.
, and
Mura
,
T.
,
1975
, “
Two-Ellipsoidal Inhomogeneities by the Equivalent Inclusion Method
,”
ASME J. Appl. Mech.
,
42
(
4
), pp.
847
852
. 10.1115/1.3423718
25.
Mura
,
T.
,
1987
,
Micromechanics of Defects in Solids
, 2nd ed.,
Springer Netherlands
,
Dordrecht
.
26.
Nakasone
,
Y.
,
Nishiyama
,
H.
, and
Nojiri
,
T.
,
2000
, “
Numerical Equivalent Inclusion Method: a New Computational Method for Analyzing Stress Fields in and Around Inclusions of Various Shapes
,”
Mater. Sci. Eng. A.
,
285
(
1–2
), pp.
229
238
. 10.1016/S0921-5093(00)00637-7
27.
Zhou
,
Q.
,
Jin
,
X.
,
Wang
,
Z.
,
Wang
,
J.
,
Keer
,
L. M.
, and
Wang
,
Q.
,
2014
, “
Numerical Implementation of the Equivalent Inclusion Method for 2D Arbitrarily Shaped Inhomogeneities
,”
J. Elasticity
,
118
(
1
), pp.
39
61
. 10.1007/s10659-014-9477-2
28.
Toupin
,
R.
,
1962
, “
Elastic Materials With Couple-Stresses
,”
Arch. Rational Mech. Anal.
,
11
, pp.
385
414
. 10.1007/BF00253945
29.
Mindlin
,
R.
, and
Eshel
,
N.
,
1968
, “
On First-Gradient Theories in Linear Elasticity
,”
Int. J. Solid Struct.
,
4
, pp.
109
124
. 10.1016/0020-7683(68)90036-X
30.
Delfani
,
M. R.
, and
Shahandashti
,
M. L.
,
2017
, “
Elastic Field of a Spherical Inclusion With Non-Uniform Eigenfields in Second Strain Gradient Elasticity
,”
Proc. R. Soc. A: Math., Phys. Eng. Sci.
,
473
(
2205
), p.
20170254
.
31.
Delfani
,
M.
, and
Sajedipour
,
M.
,
2018
, “
Spherical Inclusion With Time-Harmonic Eigenfields in Strain Gradient Elasticity Considering the Effect of Micro Inertia
,”
Int. J. Solids. Struct.
,
155
, pp.
57
64
. 10.1016/j.ijsolstr.2018.07.008
32.
Sharma
,
P.
,
2004
, “
Size-Dependent Elastic Fields of Embedded Inclusions in Isotropic Chiral Solids
,”
Int. J. Solids. Struct.
,
41
(
22–23
), pp.
6317
6333
. 10.1016/j.ijsolstr.2004.05.004
33.
Zhang
,
X.
, and
Sharma
,
P.
,
2005
, “
Inclusions and Inhomogeneities in Strain Gradient Elasticity With Couple Stresses and Related Problems
,”
Int. J. Solids. Struct.
,
42
(
13
), pp.
3833
3851
. 10.1016/j.ijsolstr.2004.12.005
34.
Yin
,
H.
, and
Zhao
,
Y.
,
2016
,
Introduction to the Micromechanics of Composite Materials
, 1st ed.,
CRC Press Inc.
,
Boca Raton, FL
.
35.
Song
,
G.
,
Wang
,
L.
,
Deng
,
L.
, and
Yin
,
H.
,
2015
, “
Mechanical Characterization and Inclusion Based Boundary Element Modeling of Lightweight Concrete Containing Foam Particles
,”
Mech. Mater.
,
91
, pp.
208
225
. 10.1016/j.mechmat.2015.07.014
36.
Song
,
G.
, and
Yin
,
H. M.
,
2018
, “
Stress Concentration of One Microvoid Embedded in An Adhesive Layer Under Harmonic Load
,”
J. Eng. Mech.
,
144
(
3
), p.
04018002
. 10.1061/(ASCE)EM.1943-7889.0001416
37.
Wu
,
C.
, and
Yin
,
H.
,
2021
, “
The Inclusion-Based Boundary Element Method (iBEM) for Virtual Experiments of Elastic Composites
,”
Eng. Anal. Boundary Elements
,
124
, pp.
245
258
. 10.1016/j.enganabound.2020.12.020