1 Introduction

Eshelby’s inclusion problem is important in solid mechanics and materials science,since the problem can model many phenomena that involve thermal expansion and structural transfor- mations in solids. In the three-dimensional (3D) elasticity case,Eshelby^[1] proved that the elastic field in an ellipsoidal inclusion is uniform when the uniform eigenstrains are applied in the inclusion which is embedded in an infinite homogenous matrix. Similar result holds in the two-dimensional (2D) elasticity. If the uniform eigenstrains are assumed for an inclusion in an infinite plate,the inclusion itself and the surrounding matrix will be stressed. Particularly,if the inclusion has an elliptic configuration,the real strains in the inclusion are also uniform. The pioneer work for the inclusion problem was carried out by Eshelby^[1] in an earlier year. Later,Mura^[2] summarized the eigenstrain problem in more detail.

It was proved that the stress field inside the m-pointed polygonal inclusion is uniform when m is an odd number^[3]. The elastic fields in an arbitrary convex polygon-shaped inclusion with uniform eigenstrains were investigated under the condition of plane strain^[4]. The elastic fields in an infinite elastic body containing a polyhedral inclusion with eigenstrains were investigated^[5]. It was proved that polyhedral shapes are impossible for inclusions with constant eigenstresses^[6].

With the aid of the techniques of analytic continuation and conformal mapping,a novel method is presented to obtain the analytic solution to Eshelby’s problem of an inclusion in a plane or a half-plane^[7]. Properties of Eshelby’s tensor for a rotational symmetrical inclusion in the plane elasticity were studied^[8].

Eshelby’s problem of non-elliptical inclusions was studied. Many explicit analytical solutions of Eshelby’s tensor field for various non-elliptical inclusions were provided. It was concluded that Eshelby’s tensor field inside a non-elliptical inclusion is quite non-uniform and cannot be replaced by its average^[9]. The stress field inside a 2D arbitrary-shape elastic inclusion bonded through an interphase layer to an infinite elastic matrix subjected to uniform stresses at infinity was analytically studied using the complex variable method in elasticity^[10]. Within the framework of 2D or 3D linear elasticity,a general approach based on the superposition principle is proposed to study the problem of a finite elastic body with an arbitrarily shaped and located inclusion^[11].

A solution was presented for determining the stresses in an infinite elastic plate containing a rectangular inclusion subjected to a uniform stress field^[12]. The 2D mathematical model of an interface crack which lies along an elastic inclusion embedded in an elastic matrix with dif- ferent elastic constants was considered^[13]. An analytical solution in infinite series form for two circular cylindrical elastic inclusions embedded in an infinite matrix with two circumferentially inhomogeneous imperfect interfaces interacting with a circular Eshelby inclusion in anti-plane shear was derived^[14]. Based on the transformation toughening theory,an approximate solution was developed for predicting the stress intensity factor for a crack interacting with an inclusion of arbitrary shape and size under I/II mixed mode loading conditions^[15].

A boundary-domain integral equation in elastic inclusion problems was introduced^[16]. In the formulation,the displacement integral equation was applied to all the boundary nodes,and the strain integral equation was used at all the collocation nodes of inclusions. A boundary integral equation approach was used to solve an infinite anisotropic elastic inclusion problem subjected to remote loading^[17]. Ting and Schiavone^[18] studied an anisotropic elastic inclusion of arbitrary shape embedded inside an infinite dissimilar anisotropic elastic medium(matrix) subjected to a uniform anti plane shear loading at infinity.

From the aforementioned publications,we see that many researchers paid attention to the structures of Eshelby tensor under different configuration of the inclusion,whereas less attention has been paid for the stress distribution in the matrix side.

This paper provides a closed form solution and numerical analysis for Eshelby’s elliptic inclusion in an infinite plate. In order to solve the problem and to formulate the conditions along the interface,the complex variable method and the conformal mapping technique are used. In the continuity condition for the displacement,one must consider the contribution from the assumed eigenstrains. In the present study,the eigenstrains take the constant values. The continuity conditions for the traction and displacement along the interface in the physical plane are reduced to the similar condition along the unit circle of the mapping plane. The properties of the complex potentials defined on the finite elliptic region are analyzed in detail. From the continuity conditions,one can separate and obtain the relevant complex potentials defined in the inclusion and the matrix.

This paper also provides a numerical analysis for Eshelby’s elliptic inclusion in plane elas- ticity. From the obtained complex potentials,the dependences of the real strains and stresses in the inclusion from the assumed eigenstrains are analyzed. In addition,the stress distribution on the interface along the matrix side is evaluated. These results are first obtained in the paper.

2 Analysis 2.1 Some preliminary knowledge in complex variable method of plane elasticity

The complex variable function method plays an important role in plane elasticity. The fun- damental of this method is introduced. In the method,the stresses (σx,σy,σxy),the resultant forces (X,Y ),and the displacements (u,v) are expressed in terms of complex potentials Φ(z) and (z) such that^[19]

where a bar over a function denotes the conjugated value for the function,G is the shear modulus of elasticity,κ = (3 − ν)/(1 + ν) in the plane stress problem,κ = 3 − 4ν in the plane strain problem,and ν is Poisson’s ratio. Sometimes,the displacements u and v are denoted by u1 and u2,the stresses σx,σy,and σxy are denoted by σ1,σ2,and σ12,and the coordinates x and y are denoted by x1 and x2,respectively.

.

The strain components are defined by

Clearly,for any assumed displacements u(x,y) and v(x,y),the relevant strain components must satisfy the following compatiblity condition: Clearly,if ε_x = c1,ε_y = c2,and ε_xy = c3,where ci are constant,those strains satisfy the compatiblity condition (5) automatically..

In addition,from Eq. (3),the relevant strain components can also be expressed as

In Eshelby’s inclusion problem,we generally assume some eigenstrains σ_x^*,σ_y^*,and σ_xy^* in the inclusion,which take constant values in the present study. Clearly,these eigenstrains ε_x^*, σ_y^* ,and ε_x^*y satisfy Eq. (5). In this case,from Eq. (4),the relevant displacements will be where the two coefficients c1 and c2 generally take complex values,which are

2.2 Closed form solution for Eshelby’s elliptic inclusion under assumed eigin- strains

For the solution to Eshelby’s elliptic inclusion problem,the following mapping function is used^[19]:

which maps the unit circle and its exterior region (in the &-plane) into the elliptic contour and its exterior region (in the z-plane) (see Fig. 1). Clearly,the mapping function has the following property:

The inverse of the mapping function of z = !(&) is denoted by

In addition,the function z = !(&) = R(& + m & ) also maps the ring region pm < |&| 6 1 (in the &-plane) into the interior region S2 with a cut line (−c,c) (see Fig. 1).

Clearly,if an analytic function f(z) is defined in the elliptic region S2 (see Fig. 1),we can obtain an analytic function

which is analytic in the ring region pm 6 |&| 6 1 (in the &-plane).

However,the inverse of the above-mentioned statement may not be correct. In fact,if an analytic function F(&) is defined in the ring region pm 6 |&| 6 1 (in the &-plane),we can formally define a function

Clearly,the defined function f(z) is analytic in the elliptic region S2 with a cut line (−c,c). However,the defined function f(z) may not be analytic in the elliptic region S2 without the cut line (−c,c)(in the z-plane). Clearly,the function f(z) may not satisfy the following condition among the cut line (−c,c):

Condition (15) can be rewritten as

From condition (16),we can get the following expression for F(&):

It is seen that if F(&) satisfies condition (16),F(&) must take the form of Eq. (17). To put it another way,if F(&) is defined in the form of Eq. (17),we can find a function f(z) = F(&)|&=Q(z), which is analytic in the interior region S2 without the cut line (−c,c) (in the z-plane).

We formulate Eshelby’s inclusion problem for an elliptic inclusion as follows (see Fig. 2). In the elliptic inclusion,the eigenstrains σ_x^*,"σ_y^*,and "σ_xy^* (=constant) are assumed beforehand. Inthe formulation,the complex potentials defined in the region z ∈ S1 and z ∈ S2 are denoted by Φ1o(z),1o(z) and Φ2o(z),2o(z),respectively.

The resultant force components and displacement components are continuous along the interface z ∈ ∑(see Fig. 2). Thus,from Eqs. (2) and (3),the continuity conditions along the interface z ∈ ∑ are as follows:

where ∑ denotes the interface on the z-plane (see Fig. 2). In Eq. (19),the term 2G(u∗ + iv∗) represents the influence from Eshelby’s eigenstrains.

Physically,the equality shown by Eq. (18) represents the continuity condition for the resul- tant force function p = −Y +iX = Φ(z)+zΦ′(z)+ (z) along the interface. Since this function, or p = −Y + iX = Φ(z) + zΦ′(z) + (z) is obtained after an integration for the traction com- ponents σN +iσNT along the interface (see Fig. 2),the component σN +iσNT in the matrix side must be equal to that in the inclusion side. However,the component σT in the matrix side may not be equal to that in the inclusion side (see Fig. 2).

Substituting Eq. (7) into Eq. (19) yields

where c1 and c2 have been defined by Eq. (8).

Considering Eqs. (18) and (20) as linear algebraic equations for Φ1o(z) and zΦ′1o(z)+ 1o(z), we have

Then,Eq. (22) can be rewritten as

In the present study,we only consider the case of the constant eigenstrains. Thus,we can assume Φ2o(z) takes the form Φ2o(z) = q1z,and 2o(z) takes the form 2o(z) = s1z. From thegeneral properties of the complex potentials shown by Eq. (17),those complex potentials must take the following expressions (taking the term for k=1 in Eq. (17)):

As for the complex potentials for z 2 S1,we assume

In Eqs. (24)-(27),q,s,pk,and rk (k=1,2,· · ·) are coefficients to be determined.

By the conformal mapping and the relations,x = a cosθ,y = b sinθ,2 cosθ = σ + σ−1, 2 sinθ = −i(σ − σ^-1),a = R(1 +m),and b = R(1 − m)(z ∈ ,or σ = e^iθ 2 ∈o,in which ∑o is the unit circle),Eqs. (21) and (23) can be rewritten as

where

Note that,for σ = e^iθ ∈∑_o,we use the relation σ = σ^-1. In addition,the term d₁σ + d₂σ^-1 is derived from (c₁x + c₂y)/(κ + 1) in Eq. (21).

By Eq. (28),Eq. (29) can be rewritten as

In the following analysis,we define the following two operators:

If one uses the integration operator I1(σ ! &) (pm < |&| < 1) to Eq. (28),we have^[19]

Thus,we can get

Substituting Eq. (18) into Eq. (35) yields

If one uses the integration operator I2(σ ! &) (for |&| > 1) to Eq. (28),we have^[19]

Substituting Eqs. (30) and (34) into Eq. (37) yields

For the solution from condition (31),we can define a function Q ∞ (&) as follows:

If one uses the integration operator I1(σ ! &) (for pm < |&| < 1) to Eq. (31),we have^[19]

Thus,from Eqs. (25),(30),and (40),we have

If one uses the integration operator I2(σ ! &) (for |&| > 1) to Eq. (31),we have^[19]

Substituting Eqs. (30),(39),and (40) into Eq. (42),we have

In the above mentioned derivation,the function 2(&) shown by Eq. (41) is regular in the region pm < |&| < 1. Thus,the function can be extended and defined on the region pm 6 |&| 6 1. Similarly,the function 1(&) shown by Eq. (43) is regular in the region |&| > 1. Thus,the function can also be extended and defined on the region |&| > 1.

Finally,we obtain the four complex potentials Φ1(&),Φ2(&),1(&),and 2(&) shown by Eqs. (38),(36),(43),and (41),respectively.

From the obtained complex potentials Φ1(&),Φ2(&),1(&),and 2(&),we can obtain the relevant complex potentials Φ1o(z),Φ2o(z),1o(z),and 2o(z) as follows:

)

Particularly,since Φ2(&) and 2(&) are expressed in the form of Eqs. (36) and (41),respec- tively,the relevant complex potentials Φ2o(z) and 2o(z) can be evaluated in a closed form as a function of z,which are as follows:

2.3 Influence on (ε_x,"ε_y, ε_xy)and( σ_x,"σ_y, σ_xyin inclusion caused by assumed eigenstrains (ε_x^*,ε_y^*,and ε_xy^*

In the following analysis,we shall derive the dependences of (ε_x,"ε_y, ε_xy) and (σ_x,"σ_y, σ_xy) in the inclusion on the assumed eigenstrains ("σ_x^*,"σ_y^*,and "σ_xy^*).

Note that 2G(u∗+iv∗) in the right hand term of Eq. (19) corresponds to relevant eigenstrains (σ_x^*,σ_y^*,and σ_xy^*) in the inclusion. In addition,the term corresponds to relevant strains in the inclusion.

Using Eqs. (6),(48),and (49),we can find the strains from the complex potentials Φ2o(z) and 2o(z) as follows:

where ε_xm,ε_ym,and ε_xym denote the modified portion of strains from the assumed eigenstrains σ_x^*,"σ_y^*,and "σ_xy^*.

Clearly,the real strains in the inclusion will be

Substituting Eqs. (50)-(52) into Eq. (53) yields

where

For b/a=0.1,0.2,· · ·,1.0 (note that m = (a − b)/(a + b)),the computed results for g11,g12, g21,g22,and g33 are listed in Table 1. Clearly,the coefficients g11,g12,g21,g22,and g33 are the counterparts of Eshelby’s tensor in literature.

From the tabulated results,we see that g11 and g22 are within the range 0 < g11 < 1 and 0 < g22 < 1. This phenomenon is easy to be explained. For example,in the case of b/a=0.5, ε_x^* = 1,ε_y^* = 0,and ε_x^*y = 0,we have ε_x = 0.492 1,ε_y = 0.127 0,and ε_xy = 0. Since the eigenstrian ε_x^* = 1 is an elongation in the x-direction and the elongation is interrupted by the surrounding matrix,we have the real strain ε_x = 0.492 1 < ε_x^* with ε_x^* = 1.

Similarly,by Eqs. (1),(48),and (49),the stress components in the inclusion can be expressed as

where

Note that those stress components really exist in the inclusion. For b/a=0.1,0.2,· · ·,1.0 (note that m = (a − b)/(a + b)),the computed results for h11,h12,h21,h22,and h33 are listed in Table 2.

From the tabulated results,we see that all the hij values are varied within the range hij < 0. This phenomenon is easy to be explained. For example,in the case of b/a=0.5,ε_x^* = 1,σ_y^* = 0, and ε_x^*y = 0,we have σx = −1.587 3G,σy = −0.317 5G,and σxy = 0. Since the eigenstrian ε_x^* = 1 is an elongation in the x-direction and the elongation is interrupted by the surrounding matrix,we have the real stresses σx = −1.587 3G,σy = −0.317 5G,and σxy = 0.

2.4 Influence at interface boundary along matrix side caused by assumed eigen- strains ε_x^*,"ε_y^*,and "ε_xy^*

For the matrix portion,we have obtained the complex potentials Φ1(&)and 1(&) shown by Eqs. (38) and (43). For any point z (z 2 S1,z = !(&)),from Eq. (2),we can evaluate the stress components as follows:

For the case of assumed eigenstrains ε_x^*,σ_y^*,and ε_x^*y,at the boundary point of matrix side x = a cosθ,y = b sinθ,the stress components σx,σy,and σxy can be evaluated using Eqs. (58)and (59). Clearly,the stress components σN,σNT,and σT at the boundary point of matrix side (x = a cosθ,y = b sinθ,and for details see Fig. 2) can also be evaluated,which are expressed as

The computed values for {f1(θ),g1(θ)},h1(θ),{f2(θ),g2(θ)},h2(θ),{f3(θ),g3(θ)},and h3(θ) are plotted in Figs. 3-8.

In the case of the applied eigenstrain ε_x^*y (ε_x^* = 0 and σ_y^* = 0),we find the following results from Fig. 7 and Fig. 8. Generally,σN and σNT take the values within the range |σN|<0.8Gε_x^*y and |σNT| < 0.8Gε_x^*y. For the σT components,we find the following maximum: (1) in the case of b/a=0.25,σT = −2.400Gε_x^*y atθ = 14◦; (2) in the case of b/a= 0.50,σT = −2.222Gε_x^*y at θ = 27◦; (3) in the case of b/a= 0.75,σT = −2.157Gε_x^*y atθ = 37◦; and (4) in the case of b/a= 1.00,σT = −2.143Gε_x^*y atθ = 45◦.

In the case of the applied eigenstrain ε_x^* (σ_y^* = 0 and ε_x^*y = 0),we find the following results from Fig. 3 and Fig. 4. In the case of b/a=0.25,we find the maximum σN = −2.057Gε_x^* atθ = 0◦, and the maximum σNT = 0.914Gε_x^* atθ = 14◦. In the case of b/a=1.00,we find the maximum σT = 1.786Gε_x^* atθ = 90◦.

In the case of the applied eigenstrain σ_y^* (ε_x^* = 0 and ε_x^*y = 0),we find the following results from Fig. 5 and Fig. 6. In the case of b/a=1.00,we find the maximum σN = −1.071Gσ_y^* atθ = 90◦, and maximum σNT = −0.357Gσ_y^* atθ = 46◦. In the case of b/a=0.25,we find the maximum σT = 2.514Gσ_y^* atθ = 0◦. That is to say,the crown point (x = a and y =0) is under a dangerous situation.

3 Conclusions

The conformal mapping function method plays an important role in the present study. The behaviors of the complex potentials Φ1(&) defined in |&| > 1 and Φ2(&) defined in pm 6 |&| 6 1 are studied exactly. In order that is an analytic function in the region z 2 S2,the function Φ2(&) must be expressed in the form of Eq. (17). This is an essential step in the present study. Similar property exists for the complex potentials ψ₂(&).

However,if one expands both sides in Eqs. (28) and (29) (including the complex potentials Φ2(σ) and φ₂(σ),σ = e^iθ ∈ ∑₀) in the Laurent series,the derivation is rather complicated. However,we use the integral operators (for ) and defined by Eqs. (32) and (33) to both sides of Eqs. (28) and (29),we can obtain all complex potentials in the closed form with higher efficiency.

Finally,the complex potentials Φ1(&) and 1(&) in the matrix and Φ2(&) and 2(&) in the inclusion are evaluated in closed form,which are shown by Eqs. (37),(42),(35),and (41), respectively. That is to say,the strain and stress distributions in the matrix and inclusion can be evaluated exactly.

Many previous researchers paid attention to the influence of the eigenstrains under different configurations of the inclusion. Even for the elliptic inclusion case,the detailed expressions for the influence on (ε_x,ε_y,ε_xy) and (σx,σy,σxy) in the inclusion caused by the assumed eigenstrains (σ_x^*,"σ_y^*,and "σ_xy^*) have not been reported. However,this paper provides a complete derivation for those influences,which are shown by Eqs. (54) and (56) and are listed in Tables 1 and 2.

In addition,many elastic responses caused by Eshelby’s eigenstrains,for example,the stress components σN,σNT,and σT (see Fig. 2) along the interface at the matrix side,are obtained in this paper for the first time.