Numerical solution of Eshelby's elastic inclusion problem in plane elasticity by using boundary integral equation

Engineering Analysis with Boundary Elements 68 (2016) 17–23

Contents lists available at ScienceDirect

Engineering Analysis with Boundary Elements journal homepage: www.elsevier.com/locate/enganabound

Numerical solution of Eshelby's elastic inclusion problem in plane elasticity by using boundary integral equation Y.Z. Chen Division of Engineering Mechanics, Jiangsu University, Zhenjiang, Jiangsu, 212013 PR China

art ic l e i nf o

a b s t r a c t

Article history: Received 23 December 2015 Received in revised form 18 March 2016 Accepted 22 March 2016

This paper provides a numerical solution for Eshelby's elastic inclusions in an plane elasticity based on the complex variable boundary integral equation (CVBIE) method. An inclusion with arbitrary shape is embedded in the infinite matrix. In the inclusion, the constant eigenstains are assumed. No remote loading is applied to the matrix. The displacements from the assumed eigenstrains are evaluated exactly, which in turn are the generalized loading in the problem. The proposed problem is reduced to solve interior and exterior boundary value problems simultaneously. For the elliptical inclusion case, the computed stresses along the interface of the inclusion side are nearly uniform. One more numerical example is devoted to a square-type inclusion with round corner. & 2016 Elsevier Ltd. All rights reserved.

Keywords: Eshelby's inclusion in plane elasticity Complex variable boundary integral equation method Stress distribution along interface Elliptic inclusion Square-type inclusion with round corners

1. Introduction In an earlier year, Eshelby [1] published a pioneering work in the field of inclusion problem in elasticity. In the problem, an inclusion is embedded in the infinite homogenous medium. The uniform eigenstrains are prescribed in the inclusion. In the 3D elasticity case, Eshelby 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. In addition, Mura [2] studied and summarized the eigenstrain problem in more detail. For the inclusion problems, Mura [2] proposed the following category, the problems for (1) inhomogeneous inclusions, (2) inhomogeneities and (3) homogeneous inclusions. The inhomogeneous inclusion is defined such that it has different elastic constants with surrounding matrix and contains eigenstarins. The inhomogeneity is defined such that it has different elastic constants with surrounding matrix. The homogeneous inclusion is defined such that it contains eigenstarins. In the Eshelby's inclusion problem, the explicit form of the Eshelby tensor was derived [3]. Properties of the Eshelby tensor for a rotational symmetrical inclusion are studied. A new integral equation formulation of two-dimensional infinite isotropic medium with various inclusions and cracks was presented [4]. The derivation depends on some expressions for displacements and E-mail address: [email protected] http://dx.doi.org/10.1016/j.enganabound.2016.03.014 0955-7997/& 2016 Elsevier Ltd. All rights reserved.

tractions at domain point. Eshelby's tensor fields in transport phenomena and anti-plane elasticity was studied with the purpose of finding their general properties and structure [5]. The null-field integral equation was suggested, which was used for the solution of an infinite medium containing circular holes and/or inclusions [6]. Using the suggested method, the multi-inclusion problem under antiplane shear was solved numerically. Wang considered the internal stress field of a three-phase elliptical inclusion bonded to an infinite matrix through an interphase layer when the matrix is subjected to remote uniform stresses [7]. A condition leading to internal uniform hydrostatic stresses was derived. Uniform stress state inside an inclusion of arbitrary shape in a three-phase composite was studied [8]. In the inclusion, the stresses are assumed to be hydrostatics. Complex variable method is used in the study. A comprehensive review for recent works on inclusions was provided [9]. The study of inclusions is of significance to the development of advanced materials. In addition, the review concludes with an outlook on future research directions. For the Eshelby's elliptic inclusion in plane elasticity with the constant distribution of the eigenstrains, several closed form solutions were provided [10]. Those solutions were obtained by using the complex variable, and they were expressed in an exact form. Lee et al. [11] provided a closed-form solution for the arbitrary polygonal inclusion problem with polynomial eigenstrains of arbitrary order in an anisotropic magneto-electro-elastic full plane. Physically, when the matrix and the inclusion are under different

18

Y.Z. Chen / Engineering Analysis with Boundary Elements 68 (2016) 17–23

temperature distributions, the relevant thermo-inclusion problem can be reduced to a particular Eshelby's inclusion problem. From those references we see that the inclusion problems have not been solved very well previously. For example, no general numerical method was provided to find out the stress distribution for an arbitrary shape of inclusion. This paper provides a numerical solution for Eshelby's elastic inclusions in plane elasticity based on the complex variable boundary integral equation (CVBIE) method. An inclusion with arbitrary shape is embedded in the infinite matrix. In the inclusion, the constant eigenstains are assumed. No remote loading is applied to the matrix. The displacements from the assumed eigenstrains are evaluated exactly, which in turn are the generalized loading in the problem. The proposed problem is reduced to solve interior and exterior BVPs (boundary value problem) simultaneously. For the elliptical inclusion case, the computed stresses along the interface of the inclusion side are nearly uniform. One numerical example is devoted to a square-type inclusion with round corner. Many computed results are provided.

∂ 2εy* ∂x2

+

∂ 2εx* ∂y2

−2

* ∂ 2εxy ∂x∂y

=0

(5)

* = c3( ci -constant), those strains Clearly, if εx* = c1, εy* = c2 , εxy satisfy the compatible condition (5) automatically. Now we study the Eshelby's inclusion problem in plane elasticity. Without losing generality, we take the elliptic inclusion as an example (Fig. 1). The inclusion portion is denoted by “ S1”, and the matrix portion is denoted by “ S2 ”. In the inclusion S1, the ei* are assumed. genstarins εx*, εy* and εxy The original problem can be decomposed into two BVPs (Fig. 2). In the interior BVP shown in Fig. 2(a), the applied boundary displacement is denoted by u1 + iv1 and the boundary traction is denoted by σN1 + iσNT 1. The complex potentials defined in the inclusion are denoted by ɸ1(z ) and ψ1 (z ). In the Eshelby's inclusion problem, we generally assume some * placed in the inclusion. It is assumed eigenstrains εx*, εy* and εxy * satisfy Eq.(5), and the relevant that those eigenstrains εx*, εy* and εxy displacement is denoted by u1* + iv1*. In this case, Eq. (3) will be modified to

2G (u1 + iv1) = κ ɸ1(z ) − z ɸ′1 (z ) − ψ1 (z ) + 2G (u1* + iv1* )

2. Analysis Analysis presented below mainly depends on two kinds of integral equation. Among them, one is for the interior BVP (boundary value problem) and other for the exterior BVP. Both of them are presented in the complex variable form [12]. For the problem of infinite plate with Eshelby's eigenstrain in inclusion, after discretization the relevant BIEs (boundary integral equation) are converted into algebraic equations, which are expressed in a matrix form. The way for the solution of the simultaneous algebraic equations is studied in details. 2.1. Some preliminary knowledge in complex variable method of plane elasticity The complex variable function method plays an important role in plane elasticity. Fundamental of this method is introduced below. 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 [13]

σx + σy = 4 Re ɸ′(z ),

σy − σx + 2iσxy = 2 [z¯ ɸ″(z ) + ψ ′(z )]

p = − Y + iX = ɸ(z ) + z ɸ′(z ) + ψ (z )

2G (u + iv) = κ ɸ(z ) − z ɸ′(z ) − ψ (z )

(1)

(2)

(3)

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 the Poisson's ratio. Sometimes, the displacements u and v are denoted by u1 and u2, the stresses σx , σy and σxy by σ1, σ 2 and σ12, the coordinates x and y by x1 and x2. It is assumed that, a particular displacement field ( u*, v*) is given beforehand. Therefore, we can define the relevant strains as

εx* =

∂u* * ∂v* * ∂u* ⎞ 1 ⎛ ∂v* + , εy = , εxy = ⎜ ⎟ ∂x ∂y ∂y ⎠ 2 ⎝ ∂x

Eq. (6) can be rewritten as

(u1 + iv1) = (u1e + iv1e ) + (u1* + iv1* )

(7)

where

(u1e + iv1e ) =

1 {κ ɸ1(z ) − z ɸ′1 (z ) − ψ1 (z )} 2G

(8)

In Eq. (7), (u1 + iv1) represents the total displacement, (u1* + iv1* ) represents the displacement caused by the eigenatrains, and (u1e + iv1e ) represents the modified part (or the elastic part) of the displacement in the solution. Note that only the portion (u1e + iv1e ) has the relevant stress in the inclusion. * taking constant value, we have In the case of εx*, εy* and εxy

* ) x + (εxy * + iεy* ) y u1* + iv1* = (εx* + iεxy

(9)

For the matrix portion, we will formulate the exterior BVP shown in Fig. 2(b). The applied boundary displacement is denoted by u2 + iv2 and the boundary traction is denoted by σN 2 + iσNT 2. The complex potentials defined in the matrix are denoted by ɸ2 (z ) and ψ 2 (z ) . In this case, Eq. (3) can be written as

2G (u2 + iv2 ) = κ ɸ2 (z ) − z ɸ′2 (z ) − ψ2 (z )

(10)

Finally, we can propose the following continuity conditions along the interface Γ (Figs. 1,2)

σ N1 + iσ NT 1 = σ N2 + iσ NT 2, (along the interface Γ )

(11)

u1 + iv1 = u2 + iv2, (along the interface), Γ

(12)

y

S2 b o * x

* y

S1 a

* xy

(4)

Clearly, those strain components must satisfy the following compatible condition

(6)

Fig. 1. Eshelby inclusion problem.

x

Y.Z. Chen / Engineering Analysis with Boundary Elements 68 (2016) 17–23

19

Fig. 2. Decomposition of the Eshelby's eigenstrain problem: (a) an interior BVP with assumed eigenstrains, (b) a usual exterior BVP.

2.2. Complex variable boundary integral equations for interior region and exterior region in the case of assumed eigenstrains in the inclusion In the present study, a CVBIE for the interior problem is introduced below [12] (Fig. 2(a)) U (to ) + B1i 2 = B2 i

∫Γ

∫Γ

⎛ κ−1 ⎞ U (t ) dt − L1 (t , t o )U (t ) dt + L2 (t , t o )U (t )dt ⎟ ⎜ ⎝ t − to ⎠

for the interior problem)

(14)

In Eq. (14), u(t) and v(t) take the real value and U(t) ¼u(t) þiv (t) is a complex value. Similarly, σN (t ) and σNT (t ) take the real value and Q (t ) = σN (t ) + iσNT (t ) is a complex value. Those notations have been indicated in Fig. 2(a). In addition, two elastic constants and two kernels are defined by

1 1 B1 = , B2 = 2π (κ + 1) 4πG (κ + 1)

(15)

⎛

⎞

∫ Γ ⎜⎝ tκ −− t1o U (t ) dt − L1 (t , t o)U (t ) dt + L2 (t , t o)U (t )dt ⎟⎠ ⎛

⎞

∫ Γ ⎜⎝ 2κ ln t − to Q (t ) dt + tt¯ −− tt¯oo Q (t )dt¯⎟⎠,

( to ∈ Γ , (13)

where Γ denotes the boundary of the interior region and the increase “dt” is defined in the anti-clockwise direction. Generally, the increase “dt” takes a complex value, which is indicated in Fig. 2 (a). In addition, dt¯ is a conjugate value with respect to the increase “dt”. In Eq. (13), U(t) and Q(t) denote the displacement and traction along the boundary Γ , which are defined by

U (t ) = u (t ) + iv (t ), Q (t ) = σ N (t ) + iσ NT (t ), (t ∈ Γ )

U (to ) − B1i 2 = − B2 i

⎛ ⎞ t − to Q (t )dt¯⎟ , ⎜ 2κ ln t − to Q (t ) dt + ⎝ ⎠ t¯ − t¯o

( to ∈ Γ ,

Generally, if one uses the BIE shown by Eq. (13) to an exterior boundary value problem, the increase “dt” should be going forward in a clockwise direction [12]. However, it is preferable to define increase “dt” in the anti-clockwise direction. In the case for the exterior boundary value problem, from Eq. (13) the relevant BIE for the exterior BVP shown in Fig. 2(b) should be written as

for the exterior problem)

(17)

As stated previously, for the interior problem with the imposed eigenstrains, only the part of displacement defined by (u1e + iv1e )( =(u1 + iv1) − (u1* + iv1* ), refer to Eq. (7) ) has the relevant contribution to the stress. Therefore, for the integral equation formulated to the interior BVP, we must input the (u1e + iv1e ) for U (t) in Eq. (13). In the derivation, the vector for modified part of displacement (u1e + iv1e ) is denoted by {U1e } ( ={U1} − {U *1}). Therefore, after discretization to Eq. (13), we have

1 1 {U1} + [H ]{U1}− {U1* }−[H ]{U1* } = [R ]{Q 1} , 2 2 (along Γ of the inclusion side)

(18)

If we define

[Hin ] {U1} =

1 {U1} + [H ] {U1} 2

(19)

Eq. (18) can be rewritten as

[Hin ] {U1} = [Hin ] {U *1} + [R ] {Q 1} , L1 (t , τ ) = − L2 (t , τ ) =

d dt

d dt

{

ln

t−τ t¯ − τ¯

1 1 dt¯ , + t−τ t¯ − τ¯ dt 1 t − τ dt¯ = − t¯ − τ¯ (t¯ − τ¯ )2 dt

t−τ t¯ − τ¯

{ }

}

(along Γ of the inclusion side)

= −

(20)

In (Eqs. (18) to 20), we define

(16)

where κ = 3 − 4ν (for plane strain condition), κ = (3 − ν ) /(1 + ν ) (for plane stress condition), G is the shear modulus of elasticity, and ν is the Poisson's ratio. In this paper, the plane strain condition and ν = 0.3 are assumed. In Eq. (16), τ denotes a domain point or a point on the boundary.

u1 (M ) v1 (M ) }T

{U1} = { u1 (1) v1 (1)

. .

u1 (j) v1 (j)

. .

. .

{U1* } = { u1*(1) v1*(1)

. .

u1*(j) v1*(j)

. .

. . u1*(M ) v1*(M ) }T

{Q1} = { σ N1 (1) σ NT 1 (1)

. . σ N1 (j) σ NT 1 (j)

. .

. .

σ N1 (M ) σ NT 1 (M ) }T

(21) (22)

(23)

Note that, for example, u1 (j ) and v1 (j ) represent the displacement

20

Y.Z. Chen / Engineering Analysis with Boundary Elements 68 (2016) 17–23

at the j-th node. In (Eqs. (21) to 23), the vector {U1} corresponds to (u1 + iv1), the vector {U1* }corresponds to (u1* + iv1* ), and the vector {Q 1} corresponds to (σN1 + iσNT 1) (Fig. 2(a)). * taking constant value, from Eq. (19) In the case of εx*, εy* and εxy the vector {U1* } can be evaluated from

* ) x + (εxy * + iεy* ) y u1* + iv*1= (εx* + iεxy * y) + i (εxy * x + εy* y) = (εx* x + εxy

(24)

* are some constant In Eq. (24), the eigenstrains εx*, εy* and εxy values given beforehand. In Eq. (18), the matrix [H ] is evaluated after discretization of the integral in the left hand side of Eq. (13), or the integral κ−1 B1i ∫ ( t − t U (t ) dt − L1 (t , t o )U (t ) dt + L2 (t , t o )U (t )dt ), and the matrix Γ

B2 i

⎛

Γ

o

In the matrix side, the displacement can be written as (u2 + iv2 ). Clearly, there is no eigenstain assumed in the matrix portion Therefore, after making discretization to the exterior BVP Eq. (17), we have

1 {U2 }−[H ] {U2 } = −[R ] {Q 2 } , (along Γ of the matrix side) 2

(25)

or

[Hex ]{U2 } = −[R ]{Q 2 } , (along Γ of the matrix side)

(26)

where the matrix [Hex ] is defined by

[Hex ] {U2 } =

1 {U2 }−[H ] {U2 } 2

(27)

In (Eqs. (26) and 27), we define {U 2 } = { u2 (1) v 2 (1)

. .

u 2 (j ) v 2 (j )

. .

. .

. . σ N2 (j) σ NT 2 (j) {Q 2 } = { σ N2 (1) σ NT 2 (1) . . σ N2 (M ) σ NT 2 (M ) }T

u2 (M ) v 2 (M ) }T

a dg1 =

⎛ 1−δ ⎞ ⎛ ⎞ 2 b exp ⎜ ⎟ , bdg1 = δa dg1, ⎜ with δ = , κ = 1.8⎟ ⎝ ⎠ 1+δ a ⎝ 2κ (1 + δ ) ⎠

(35)

a dg2 =

⎛ ⎛ ⎞ 2 1−δ ⎞ b exp ⎜ − ⎟ , b dg2 =δa dg2, ⎜ with δ = , κ = 1.8⎟ ⎝ ⎠ 1+δ a ⎝ 2κ(1 + δ) ⎠

(36)

(29)

(30)

Finally, the simultaneous algebraic equations from the BIE formulation for the Eshelby’ inclusion problem are composed of (Eqs. (20) and 30). From (Eqs. (19) and 27), we have the following property

(31)

From this property and Eqs. (20) and (30), we have

{U1} = [Hin ] {U *1} (along Γ)

(32)

From Eq., (31) we see that the vector {U1}( {U2 } = {U1}) can be directly obtained from the displacement vectors {U *1} caused by eigenstrain. In addition, from Eq. (30) we have

{Q 1} = −[R −1] [Hex ] {U1} (along Γ)

For the case of b/a ¼0.5, from (Eqs. (35) and 36), we have the following degenerate scale adg1 = 1.4626, bdg1 = 0.7313 and adg1 = 1.2154 , bdg1 = 0.6077. In this case, if we use a = 1.4626 and b = 0.7313 (or a = 1.2154 , b = 0.6077) in the computation, the obtained solution must be unstable or illness. In fact, we use a¼ 10 and b¼ 5 (for case b/a¼ 0.5) in the numerical computation, the relevant solution must be unique. In the second case, if the contour Γ is not a degenerate scale, the operator in the left hand side of (34) is invertible. Alternatively speaking, the BIE shown by Eq. (34) has a unique solution. If we explain this result after discretization of BIE, the inverse matrix [R −1] exists.

3. Numerical examples

In (Eqs. (28), 29), the vector {U2 } corresponds to (u2 + iv2 ), and the vector {Q 2 } corresponds to (σN 2 + iσNT 2 ) (Fig. 2(b)). Considering the continuity conditions along the interface shown by (Eqs. (11) and 12), Eq. (26) can be rewritten as

[Hin ] {U1} + [Hex ] {U1} = {U1}

(34)

(28)

. .

[Hex ] {U1} = −[R ] {Q 1} , (along Γ )

for the interior problem)

There are two particular cases for consideration. In the first case, if the contour Γ is a degenerate scale, the operator in the left hand side of (34) is not invertible. Alternatively speaking, the BIE shown by Eq. (34) has no unique solution. If we explain this result after discretization of BIE, the inverse matrix [R −1] does not exist. In the case of the elliptic configuration with two semi-axes “a” and b”, two sets of the degenerate scale were obtained previously [14]

o

[R] is evaluated after discretization of the integral in the right hand of Eq. (13) or the integral B2 i ∫ (2κ ln t − to Q (t ) dt + tt¯ −− tt¯o Q (t )dt¯ )

⎞

∫ Γ ⎜⎝ 2κ ln t−to Q (t ) dt + tt¯−−tt¯oo Q (t )dt¯⎟⎠ = f (to ), ( to ∈ Γ,

(33)

We know that, if the used scale does not reach the degenerate scale, the inverse matrix [R −1] exists. Finally, the problem is solved. The degenerate scale problem is generally arising from the Dirichlet boundary value problem. If substituting the boundary value for U(t) in the left hand side of Eq. (13), we will find the following integral equation

Example 1 In the first example, the infinite plate contains an elliptic inclusion with the applied eigenstrains εx* and εy*, which take constant value. The elliptic interface boundary Γ has two semi-axes “a” and “b”. ( Fig. 2). In computation, M¼ 96 divisions are used for the discretization of the contour Γ . The constant value for U(t) or Q (t) is assumed along all boundary elements. As stated previously, from Eqs. (32) and (33) we can evaluate two vectors {U1} = {U2 } and {Q 1} = {Q 2 }. In the example, after evaluating the vectors {U1} = {U2 } and {Q 1} = {Q 2 }, we can obtain the stresses σN1 = σN 2 and σNT 1 = σNT 2 along the interface (Fig. 2). In addition, we can also evaluate the stress components σ T1 at the inclusion side and σ T 2 at the matrix side (Fig. 2). As stated previously, in the inclusion side the vector for elastic part of displacement (u 1e+iv1e ) was denoted by {U 1e}, which can be obtained by the relation {U1e } = {U1} − {U *1}. From the two computed vectors {U 1e} and {Q 1}, we can evaluate the stress component σ T1 at the inclusion side. Similarly, from the two computed vectors {U2 } and {Q 2 }, we can evaluate the stress component σ T 2 at the matrix side. The technique for evaluating the component σ T1 and σ T 2 can be referred to [15]. It is known that, if some constant eigenstrians εx* and εy* are assumed in the elliptic inclusion, the relevant stresses and strains are still uniform in the inclusion [1,2]. A detailed solution was obtained as follows [10]

σx =

G G (H11εx* + H12 εy* ), σy = (H21εx* + H22 εy* ), κ+1 κ+1

where

(37)

In the example, for the following cases: b /a = 0.25, 0.5, 0.75 and 1.0, the computed stress components σx and σy along the contour Γ at the inclusion side are denoted by G (h11 (b/ a, θ)ε*x +h12 (b/ a, θ)ε *y ), κ+1 G σy = (h21 (b/ a, θ)ε*x +h22 (b/ a, θ)ε *y ), (at the point x = acosθ , y = b sin θ) κ+1 σx =

(39)

The computed non-dimensional stresses, or h11 (b /a, θ )and H11, h21 (b /a, θ ) and H21, h12 (b /a, θ ) and H12, and h22 (b /a, θ ) and H22 are plotted in Figs. 3–6, respectively. From plotted results we see that all those non-dimensional stresses take the negative value in general. It is seen that the deviation of the numerical solution shown by h11 (b /a, θ ) from the exact solution H11 is rather small. That is to say a rather higher accuracy has been achieved in the present method. In the case of εx* and εy* = 0, the h11 (b /a, θ )and H11 values are

b/a=0.25

in inclusion side

(38)

21

-0.5

-0.6

along

H11 = − (1 + m)(3 + m), H12 = H21 = − (1 − m2), ⎛ a − b⎞ H22 = − (1 − m)(3 − m), ⎜ with m = ⎟ ⎝ a + b⎠

Non-dimensional stresses h 21(b/a, ) , H21

Y.Z. Chen / Engineering Analysis with Boundary Elements 68 (2016) 17–23

-0.9

-0.7

at x=a cos y=b sin

h21(b/a, )

-0.8

- - numerical

H21

b/a=0.5

exact

b/a=0.75 -1.0

b/a=1 -1.1 0

15

30

45

60

75

90

(degree) Fig. 4. Non-dimensional stresses h21 (b/a, θ ) (numerical) and H21(exact) along Γ in inclusion side (see Eqs. (37),(38),(39) and Fig. 2).

dominant, and h21 (b /a, θ ) and H21 values are minor. We generally have H21/H11 = (1 − m) /(3 + m). Particularly, in the case of b/ a ¼0.25, we have H11 = − 5.760 and H21 = − 0.640. Note that the h11 (b /a, θ )and H11 values denote the influence to the component σx in the inclusion from the eigenstrain εx*. In addition, the h21 (b /a, θ ) and H21 values denote the influence to the component σy in the inclusion from the eigenstrain εx*. In the case of b/a ¼0.5, from Fig. 3 we see that a close agreement has been found between h11 (b /a, θ )values and H11( H11 = − 4.444 ). From Fig. 4, some difference can be found between h21 (b /a, θ ) and H21( H21 = − 0.888). However, if one comparing the value h21 (b /a, θ ) − H21 with the value H11( H11 = − 4.444 ), the computation error presented in Fig. 4 is small. Similarly, in the case of εy* and εx* = 0, the h22 (b /a, θ ) and H22

Fig. 5. Non-dimensional stresses h12 (b/a, θ ) (numerical) and H12 (exact) along Γ in inclusion side(see Eqs. (37),(38),(39) and Fig. 2).

-2.5

b/a=1

in inclusion side

-3.0

along

Non-dimensional stresses h11(b/a, ) , H11

values are dominant, and h12 (b /a, θ ) and H12 values are minor. We generally have H12/H22 = (1 + m) /(3 − m). Particularly, in the case of b/a¼ 0.25, we have H22 = − 0.960 and H12 = − 0.640. In the example, the computed stress components σN1 = σN 2, σNT 1 = σNT 2, σ T1 and σ T 2 along the contour Γ are expressed as (Fig. 2)

b/a=0.75 -3.5

-4.0

at x=a cos y=b sin

b/a=0.5

-4.5

h11(b/a, )

-5.0

- - numerical

H11

exact

b/a=0.25

-5.5

-6.0 0

15

30

45

60

75

90

(degree) Fig. 3. Non-dimensional stresses h11 (b/a, θ ) (numerical) and H11 (exact) along Γ in inclusion side (see Eqs. (37),(38),(39) and Fig.2).

Fig. 6. Non-dimensional stresses h22 (b/a, θ ) (numerical) and H22(exact) along Γ in inclusion side (see Eqs. (37),(38),(39) and Fig. 2).

22

Y.Z. Chen / Engineering Analysis with Boundary Elements 68 (2016) 17–23

G G ( f1 (b/a, θ) ε*x + g1 (b/a, θ) ε*y ) , σ NT1 = ( f2 (b/a, θ) ε*x + g2 (b/a, θ) ε*y ) κ+1 κ+1 G G σ T1 = ( f3 (b/a, θ) ε*x + g3 (b/a, θ) ε*y ) , σ T2 = ( f4 (b/a, θ) ε*x + g 4 (b/a, θ) ε*y ) κ+1 κ+1 (at the point x = a cosθ , y = b sin θ) σ N1 =

(40)

For the case of b/a ¼0.5 the computed non-dimensional stresses for (a) f1 (b /a, θ ), f2 (b /a, θ ), f3 (b /a, θ ) and f4 (b /a, θ ) (b) g1 (b /a, θ ), g2 (b /a, θ ), g3 (b /a, θ ) and g4 (b /a, θ )are plotted in Figs. 7 and 8, respectively. In the case of εx*, εy* = 0 and b/a ¼0.5, the f1 (b /a, θ )(for σN1 = σN 2) generally take a negative value (Fig. 7). In addition, the f3 (b /a, θ )(for σ T1) also take a negative value. The major portion for f4 (b /a, θ ) (for σ T 2in the range 15o ≤ θ ≤ 90o ) take a positive value, and f4 (b /a, θ ) max =f4 (b /a, θ ) θ= 90o ¼3.563. In the case of εy*, εx* = 0 and b/a¼ 0.5, the g1 (b /a, θ )(for σN1 = σN 2) generally take a negative value (Fig. 8). In addition, the g3 (b /a, θ )(for σ T1) also take a negative value. The major portion for g4 (b /a, θ )(for σ T 2in the range 0o ≤ θ ≤ 52o ) take a positive value, and g4 (b /a, θ ) max =g4 (b /a, θ ) θ= 0o ¼6.216. Example 2 The second example is devoted to find out the stress distribution in a square-type inclusion with the round corners when the constant eigenstains εx* and εy* are applied in the inclusion (Fig. 9). The square-type notch has a width “2a”, and the round corner has a radius “0.5a”. The square-type notch is approximated by many linear boundary elements with equal element length. On all the boundary elements, the constant density distributions for U(t) and Q(t) are assumed. The same method used in the first example can be used in the present example. It is known that, if some constant eigenstrians εx* and εy* are assumed in the square-type inclusion with round corners, the relevant stresses and strains are not uniform in the inclusion [1]. After some computation, in the cased of εx* = ε1and εy* = 0 the non-dimensional stress components σx and σy and σxy along the contour Γ at the inclusion side are denoted by (Fig. 9)

σx =

Gε1 Gε1 Gε1 h1 (θ ), σy = h2 (θ ), σxy = h3 (θ ) κ+1 κ+1 κ+1

they are located within the range 0.557 < h2 (θ ) < 0.797. Generally, h3 (θ )(for σxy along edge ABCD) values take a small one. In the example, the computed stress components σN1 = σN 2, σNT 1 = σNT 2, σ T1(at the inclusion side) and σ T 2 (at the matrix side) along the contour Γ are expressed as (Fig. 9)

σ N1 =

Gε1 Gε1 Gε1 f (θ ), σ NT 1 = f (θ ), σ T 1 = f (θ ), κ+11 κ+12 κ+13 Gε1 σ T2 = f (θ ) κ+14

(42)

The computed non-dimensional stresses, or f1 (θ ), f2 (θ ), f3 (θ ) and f4 (θ )are plotted in Fig. 11. The following results can be seen from the plotted results. In the case of εx* = ε1, εy* = 0, f1 (θ )(for σN1 along edge ABCD) generally take a negative value, for example, we have f1 (θ ) θ= 0o = − 2.384 and f1 (θ ) θ= 90o = − 0.577. In addition, all

f3 (θ )(for σ T1 along edge ABCD) values are negative, for example, we have f3 (θ )

θ= 0o

= − 0.581 and f3 (θ )

θ= 90o

= − 4.442. Only in the

o

condition of θ ≥ 37.5 , f4 (θ )(for σ T 2 along edge ABCD) values become positive, for example, we have

f4 (θ )

θ= 37.5o

= 0.688 and

(41)

The computed non-dimensional stresses, or h1 (θ ), h2 (θ ) and h3 (θ ) are plotted in Fig. 10. The following results can be seen from the plotted results. In the case of εx* = ε1 and εy* = 0, h1 (θ )(for σx along edge ABCD) generally take a negative value (Fig. 10), for example, we have h1 (θ ) θ= 0o = − 2.380 and h1 (θ ) θ= 90o = − 4.442. In addition, all h2 (θ )(for σy along edge ABCD) values are negative, and

Fig. 7. Non-dimensional stresses f1 (b/a, θ )(for σN1), f2 (b/a, θ ) (for σNT1), f3 (b/a, θ ) (for σ T1), f4 (b/a, θ ) (for σ T2 ) along Γ in inclusion side (see Eq. (40) and Fig. 2).

Fig. 8. Non-dimensional stresses g2 (b/a, θ ) (for g1 (b/a, θ ) (for σN1), σNT1), g3 (b/a, θ ) (for σ T1), g4 (b/a, θ ) (for σ T2) along Γ in inclusion side (see Eq. (40) and Fig. 2).

Fig. 9. A square-type elastic inclusion with applied eigenstrains εx* εy* embedded in an infinite plate.

Y.Z. Chen / Engineering Analysis with Boundary Elements 68 (2016) 17–23

23

inclusion, numerical examination proves that a higher accuracy has been achieved in the suggested method. The suggested method can be used to the case of non-elliptic configuration. However, the final physical stress and strain fields in the inclusion are no longer uniform.

References

Fig. 10. Non-dimensional stresses h1 (θ ) (for σx ), h2 (θ ) (for σy ), h3 (θ ) (for σxy ) along Γ in inclusion side (see Eq. (41) and Fig. 9).

Fig. 11. Non-dimensional stresses f1 (θ ) (for σN1), f2 (θ )(for σNT1), f3 (θ ) (for σ T1), along Γ in inclusion side and f4 (θ ) (for σ T2 ) in matrix side (see Eq. (42) and Fig. 9).

f4 (θ )

θ= 90o

= 3.558.

4. Conclusions This paper provides a general numerical method for evaluating the stresses in the Eshelby's inclusion problem, which is based on the usage of CVBIE. From the known solution for the elliptic

[1] Eshelby JD. The determination of the elastic field of an ellipsoidal inclusion and related problems. Proc Roy Soc Lond 1957;A 241:376–96. [2] Mura T. Micromechanics of defects in solids.Dordrecht: Martinus Nijhoff Publishers; 1987. [3] Wang MZ, Xu BX. The arithmetic mean theorem of Eshelby tensor for a rotational symmetrical inclusion. J Elast 2004;77:13–23. [4] Dong CY, Lee KY. A new integral equation formulation of two-dimensional inclusion–crack problems. Int J Solids Struct 2005;42:5010–20. [5] Le Quang H, He QC, Zheng QS. Some general properties of Eshelby's tensor fields in transport phenomena and anti-plane elasticity. Int J Solids Struct 2008;45:3845–57. [6] Chen JT, Wu AC. Null-field approach for the multi-inclusion problem under antiplane shears. ASME J Appl Mech 2007;74:469–87. [7] Wang X. Three-phase elliptical inclusions with internal uniform hydrostatic stresses in finite plane elastostatics. Acta Mech 2011;219:77–90. [8] Wang X, Gao XL. On the uniform stress state inside an inclusion of arbitrary shape in a three-phase composite. Z Angew Math Phys 2011;62:1101–16. [9] Zhou K, Hoh HJ, Wang X, Keer LM, Pang JHL, Song B. A review of recent works on inclusions. Mech Mater 2013;60:144–58. [10] Chen YZ. Closed form solution and numerical analysis for Eshelby's elliptic inclusion in plane elasticity. Appl Math Mech (Eng Ed ) 2014;35:863–74. [11] Lee YG, Zou WN, Pan E. Eshelby’s problem of polygonal inclusions with polynomial eigenstrains in an anisotropic magneto-electro-elastic full plane. Proc Math Phys Eng Sci 2011:471. http://dx.doi.org/10.1098/rspa.2014.0827. [12] Chen YZ, Lin XY. Dual boundary integral equation formulation in plane elasticity using complex variable. Eng Anal Bound Elem 2010;34:834–44. [13] Muskhelishvili NI. Some basic problems of the mathematical theory of elasticity.Groningen: Noordhoff; 1963. [14] Chen YZ, Lin XY, Wang ZX. Evaluation of the degenerate scale for BIE in plane elasticity and antiplane elasticity by using conformal mapping. Eng Anal Bound Elem 2009;33:147–58. [15] Chen YZ. Boundary integral equation method for two dissimilar elastic inclusions in an infinite plate. Eng Anal Bound Elem 2012;36:137–46.