Revisiting the problem of a 2D infinite elastic isotropic medium with a rigid inclusion or a cavity

International Journal of Engineering Science 126 (2018) 68–96

Contents lists available at ScienceDirect

International Journal of Engineering Science journal homepage: www.elsevier.com/locate/ijengsci

Revisiting the problem of a 2D infinite elastic isotropic medium with a rigid inclusion or a cavity W.-N. Zou a,∗, Q.-C. He b a

Institute for Advanced Study, Nanchang University, Nanchang 330031, China Université Paris-Est, Laboratoire de Modélisation et Simulation Multi Echelle, UMR 8208 CNRS, 5 bd Descartes, Marne-la-Vallée 77454, France

b

a r t i c l e

i n f o

Article history: Received 18 June 2017 Revised 1 January 2018 Accepted 24 January 2018

2010 MSC: 00-01 99-00 Keywords: K–M potentials Rigid inclusion Cavity Relative rigid-body displacement Laurent polynomial Stress concentration

a b s t r a c t The problem of analytically finding the elastic fields inside a 2D infinite elastic isotropic medium containing a rigid inclusion or a cavity and subjected to uniform remote loading is a classical elasticity problem of theoretical and practical interest. In the present work, we revisit the Kolosov–Muskhelishvili potential theory which is a powerful tool for solving the problem in question. In particular, a novel strategy is proposed to deal with the rigid-body displacements that the rigid inclusion or cavity may undergo. When the shape of the rigid inclusion or cavity is described by a Laurent polynomial, a general method is elaborated to solve the problem. The results obtained by using our method include as special cases all the relevant results reported in the literature. In light of our results, some errors in the literature are corrected. Finally, the cases of a rhombus rigid inclusion and a pentagonal cavity saturated with a fluid are studied in detail and some general properties are discussed. © 2018 Elsevier Ltd. All rights reserved.

1. Introduction Soft-hard materials integration, namely doping a hard phase into a soft, ductile matrix, is not only a well-established industrial technology for the production of ultra high strength materials (Michler, Adhikari, & Henning, 2004; Öztürk, Mirmesdagh, & Ediz, 1994), but also an ubiquitous method in biological materials and structures in nature and bio-inspired design (Naleway, Porter, McKittrick, & Meyers, 2015; Tang, Kotov, Magonov, & Öztürk, 2003; Wang, Yang, McKittrick, & Meyers, 2016; Wegst, Bai, Saiz, Tomsia, & Ritchie, 2015). In engineering structures, different types of holes and openings are made to satisfy some service requirements, resulting in strength degradation even failure of the structure (Penzien, 20 0 0; Rezaeepazhand & Jafari, 2010; Sharma, 2012; Zappalorto & Lazzarin, 2011). Pervasive porous microstructures in both naturally occurring and man-made materials received much attention of researchers studying the influence of their features on the local stress concentration and the effective properties of materials (Duan, Wang, & Karihaloo, 2006; Fan & Yang, 2017; Haller, Monerie, & Pagano, 2016; Mindeguia, Pimienta, Noumowe, & Kanema, 2010; Phan, 2008). Eshelby (1957, 1959) proposed an elegant method to link the contribution of inhomogeneities with different properties, now known to be suitable only for ellipsoidal/elliptical shapes (Eshelby, 1961; Kang & Milton, 2008; Liu, 2008). Dundurs (1989) showed that the stress field for cavity problem in plane elasticity can be obtained from that of rigid inclusion by setting in the solution the Kolosov ∗

Corresponding author. E-mail address: [email protected] (W.-N. Zou).

https://doi.org/10.1016/j.ijengsci.2018.01.001 0020-7225/© 2018 Elsevier Ltd. All rights reserved.

W.-N. Zou, Q.-C. He / International Journal of Engineering Science 126 (2018) 68–96

69

constant κ = −1. In the present paper, a novel procedure is developed to find the stress functions for the remotely loaded 2D infinite plate with a rigid inclusion or a cavity characterized by a conformal mapping in the form of Laurent polynomial. Even though powerful computer methods are now available for numerically solving various kinds of linear and nonlinear problems in mechanics and physics, it remains quite useful and helpful to obtain analytical solutions to some special classes of problems. This is mainly for two reasons. First, problems which admit analytical solutions can serve as benchmarks for testing the correctness and accuracy of numerical methods. Second, probably more importantly, analytical solutions even for some particular or idealized cases can often provide meaningful physical insight into the problems under investigation and may also, sometimes, reveal the mathematical structure of their general solutions. Bearing these two facts in mind, we propose, in this work, to revisit the problem of a two-dimensional (2D) infinite elastic isotropic medium containing a rigid inclusion or a cavity and subjected to remote uniform loading. This classical elastic problem is of great theoretical significance and practical value (Kachanov, Tsukrov, & Shafiro, 1994; Savin, 1961; Zhao & Yang, 2015). It is well-known that the potential theory or complex variable method of Muskhelishvili (1953) is a very efficient and elegant way for solving the problem under consideration. For the rigid inclusion or cavity which is a simply-connected region ω in a 2D space, the Riemann theorem (see, e.g., Zou, He, Huang, & Zheng, 2010) asserts that there exists a unique conformal mapping



z = h + r φ (w ) = h + r w +

∞  n=1



bn w

−n

= rw

∞   n=1

1−



wn , w

(1)

normalized by f (∞ ) = ∞ and f  (∞ ) = r > 0, that maps the exterior of the unit disk with the origin as its center on to the exterior of ω. This mapping is such that |wn | < 1 for all h ∈ ω. An alternative form can be used to map the interior of the unit disk onto the exterior of ω by the transformation w → w−1 . In particular, as recalled in Appendix B, for a polygon, the Schwartz–Christoffel transformation (143) holds. It admits the Laurent series expansion (1). The degree of accuracy of the finite description of a polygon by a Laurent polynomial is controllable; in other words, the corners of the approximated polygon can have any specified radius of curvature and its sides can be rendered as straight as desired. The problem in question was completely solved when the rigid inclusion or cavity, ω, is elliptic and has been investigated in the case where the function f (w ) is a polynomial or rational function. For example, in his monograph, Muskhelishvili (1953) collected various relevant solutions. Later, Savin (1961) decomposed the potentials into basic parts and disturbed parts, and give solutions for a number of cavity problems, many of which are important from the applicative point of view. In particular, Chang and Conway (1968) observed that all existing investigations were confined to mapping functions consisting only of three terms, so that the effect of the geometric approximation of ω, say a polygonal cavity, on stress concentration could not be studied. The effort made by Chang and Conway (1968) to improve this situation in the case where the rigid inclusion ω is characterized by a Laurent polynomial confirmed that one of the perturbation potentials can be represented by a polynomial with the coefficients obtained by solving a system of algebraic equations but failed in obtaining the Kolosov–Muskhelishvili (K–M) potentials with finite number of terms. A careful study of the relevant literature shows that, even though the problem in question is a very classical one, a general efficient solution strategy for it is still lacking when the number of terms in the mapping function is more than three. This is the origin of the present situation where one has still to solve rigid-inclusion/cavity problems case by case, sometimes in an incomplete numerical way (Ekneligoda & Zimmerman, 2008; Tsukrov & Novak, 2002; 2004). It is important to note that a one-to-one correspondence exists between the rigid inclusion problem and the cavity problem in question. This correspondence, first discovered by Dundurs (1989) and widely discussed later (see, e.g., Jasiuk, 1995; Markenscoff & Paukshto, 1995), plays an important role in, for example, doing experiments or designing composite materials and makes it possible to limit our study to the rigid inclusion or cavity problem. In the present work, we first revisit the K–M potential theory in relation to the problem under investigation. Note that relative rigid-body displacements appear generally between the rigid-inclusion/cavity and the matrix. For instance, a rigidbody rotation takes place when the matrix containing an elliptic rigid inclusion undergoes a uniaxial tension or compression loading not along with one of the principal axes of the elliptical inclusion (Muskhelishvili, 1953); a rigid-body translation occurs in the case where the matrix with a triangular cavity is under uniaxial traction (Savin, 1961). To account for these rigid-body displacements efficiently, a novel strategy is elaborated. Next, we deal in detail with the case where the shape of a rigid-inclusion or cavity is described by a Laurent polynomial with any finite number of terms. It will be shown that, in this case, the ordinary power expansions of the perturbation potentials must consist of infinite terms and admit no exact solution in a truncated way. As a solution to this difficulty, we shall suggest a novel general method to construct the perturbation potentials such that they have a finite number of terms and the homogeneous displacement or traction boundary conditions can be exactly satisfied. The rest of this paper is arranged as follows. In Section 2, we recapitulate the complex potential theory of Kolosov– Muskhelishvili, suggest a novel treatment of possible rigid-body displacements and pore pressure, and write the boundary balance relations for the rigid inclusion problem and the cavity problem in a unified form; then we elaborate the general finite analytical solution when the shape of the inclusion is characterized by a Laurent polynomial while postponing the detailed derivations to Appendix A. In Section 3, some comments are presented: the effect of a finite rigid-body rotation is considered, the solution of an elliptical inhomogeneity with eigenstrains and under remote loadings is derived based on the interior uniform stresses/strains; two kinds of correspondence for the inhomogeneity problems, including the rigid

70

W.-N. Zou, Q.-C. He / International Journal of Engineering Science 126 (2018) 68–96

inclusion and cavity problems, are discussed; finally, the transformation of potential coefficients with the rotation of the coordinate system is given. Applications are expanded in Section 4: first some reported results for the rigid inclusion and cavity problems are reexamined in light of the results extracted from our general solution; the problem of a rhombus rigid inclusion with sharp tips is further analyzed by considering the effect of finite rigid-body rotation and different levels of remote loading; In addition, the problem of a pentagonal cavity saturated with a pore pressure is studied in detail to illustrate the stress concentration along the contour of the cavity and the curvature change of the contour due to the deformation. Some concluding remarks are drawn in Section 5. Appendices B and C provide the Laurent polynomial approximation for rectangle, rhombus and pentagram by using the Schwarz–Christoffel mapping, and give the formulae for the curvature of the boundary and its change. 2. General solutions 2.1. Formulation with the K–M potentials The Kolosov–Muskhelishvili (K–M) potentials γ and ψ , two arbitrary functions of the complex variable z, have been proved to be a powerful formulation in the context of 2D isotropic elasticity (England, 1971; Lu, 1995; Muskhelishvili, 1953; Sadd, 2005). Let  be the 2D space excluding a subdomain ω which is either a rigid body or a cavity. The material forming  is linearly elastic and isotropic, so that it is characterized, for example, by its Young’s modulus E and Poisson’s √ ratio ν . Using a 2D system of Cartesian coordinates (x1 , x2 ) and introducing the complex variable z = x1 + ιx2 with ι = −1, the displacement and stress/strain components are expressed in terms of γ and ψ as

u 1 + ιu 2 =



1  κγ (z ) − zγ  (z ) − ψ (z ) , 2μ

(2)

 σ11 + σ22 = 2 γ  (z ) + γ  (z ) ,  σ22 − σ11 + 2ισ12 = 2 z¯γ  (z ) + ψ  (z ) ,

⎧ ⎪ ⎨

(3)

κ − 1  γ (z ) + γ  (z ) , 2μ  ⎪ ⎩ε22 − ε11 + 2ιε12 = 1 z¯γ  (z ) + ψ  (z ) , μ ε11 + ε22 =

(4)

where ( ) indicates the conjugate of a complex variable () and

⎧ ⎨ 3−ν E 1+ν μ= , κ= 2 (1 + ν ) ⎩3 − 4ν

in the case of plane stress;

(5)

in the case of plane strain.

The traction vector fn = f1n + ι f2n on the boundary is given by

fn = −ι

d γ ( z ) + zγ  (z ) + ψ ( z ) , ds

(6)

where s is the arc length along the boundary with n as its outward unit normal. The stress field thus obtained satisfies the equilibrium equations without body forces and the compatibility relation automatically. So, solving a 2D elastic problem amounts to finding the potentials γ and ψ such that the displacement and stress fields derived from them satisfy the prescribed  boundary∞conditions.  Let σx∞ , σy∞ , τxy be the prescribed remote uniform stresses on . From (3), we have for z → ∞,

⎧ ⎪ ⎨

1  ∞ = Cσ ,1 , γ  (z ) = σ11∞ + σ22 2   ⎪ ⎩z¯γ  (z ) + ψ  (z ) = 1 σ ∞ − σ ∞ + ισ ∞ = Cσ ,2 , 22 11 12 2Re



(7)

2

resulting in the basic potentials

1 2

γo (z ) = Cσ , 1 (z − h ),

1 2

ψo (z ) = Cσ , 2 (z − h ) − Cσ , 1 h¯ ,

(8)

without rigid-body displacement in the infinity, where the points in the matrix are indicated by (1). The perturbation potentials, denoted by γ p (z) and ψ p (z), giving rise to the stresses due to the presence of a rigid body or a cavity ω, will be determined by the homogeneous conditions on the boundary of ω. For the conformal mapping (1), the perturbation potentials γ p and ψ p , being single-valued due to the assumption that no resultant is exerted on the boundary of ω, admit in general the series expansions:

W.-N. Zou, Q.-C. He / International Journal of Engineering Science 126 (2018) 68–96

γ p (z (w ) ) = r

∞ 

αk w−k ,

ψ p (z (w ) ) = r

k=1

∞ 

71

βk w−k , |w| ≥ 1,

(9)

k=1

because they vanish for w → ∞ or z → ∞. Above, r is the size parameter of ω. Under the remote loading, the subdomain ω as a rigid inclusion may undergo an unknown rigid-body displacement, say a rigid-body translation U0 plus a rigid-body rotation α which are in general assumed to be small. Then, the displacement condition U = U0 + ια rφ (η ) on ∂ω (|η| = 1) can be achieved from (2) as

U=

κγ p (η ) −

z (η ) z¯ η−1

(

)

    γ¯ d η−1 − ψ¯ d η−1

2μ

+

  1 κ −1 Cσ ,1 rφ (η ) − C¯σ ,2 rφ¯ η−1 . 4μ 2μ

(10)

It is convenient to combine the real parameters Cσ , 1 and α so as to yield a complex:

Cσα,1 = −

κ −1 2

Cσ , 1 + 2μαι.

(11)

Thus, the basic relation (10) is simplified as

      z (η ) κγ p (η ) −   −1  γ¯ d η−1 − ψ¯ d η−1 = Cσα,1 rφ (η ) + C¯σ ,2 rφ¯ η−1 + 2μU0 . z¯ η

(12)

Otherwise, if the subdomain ω corresponds to a cavity, we start from the traction f on the boundary, whose complex expression in terms of K–M potentials is given by (6). Under the assumption that there is a pressure P applied on the boundary ∂ω such that

fn = P n = −ιP

dz , ds

(13)

the force balance on the boundary |η| = 1 can be formulated as

γ (z ) + zγ  (z ) + ψ (z ) = −Pz − f0

(14) f0 2μ .

where the counterpart of the integration constant f0 is the rigid-body translation defined by U0 = Due to the property of the surface traction on the boundary, the term associated with the rigid-body rotation naturally disappears, namely α ≡ 0. Substituting the basic potentials (8) into (14) yields

      z (η ) γ p (η ) +   −1  γ¯ p η−1 + ψ¯ p η−1 + Cσp,1 rφ (η ) + C¯σ ,2 rφ¯ η−1 + f0 = 0 z¯ η

(15)

p

with Cσ ,1 = Cσ , 1 + P . Comparing Eq. (15) for a cavity with Eq. (12) for a rigid inclusion allows establishing the equivalence between the cavity and rigid-inclusion problems upon setting (Dundurs, 1989)

α → 0, 2μU0 → f0 ; κ → −1, Cσα,1 → Cσp,1 .

(16)

Thus, we can write (12) and (15) in a unified form

      z (η ) κγ p (η ) −   −1  γ¯ d η−1 − ψ¯ d η−1 = Cσ∗ ,1 rφ (η ) + C¯σ ,2 rφ¯ η−1 + 2μU0 . z¯ η

(17)

2.2. Main results Now we are concerned with an infinite elastic body  with a rigid inclusion or a cavity defined by a Laurent polynomial, namely a truncation of (1), as:

z = h + rφN (w ),

φN ( w ) = w +

N  bk w−k , |w| ≥ 1.

(18)

k=1

Making use of a novel solution procedure we are able to derive the explicit perturbation potentials of (17). In order to highlight the results from many lengthy formulae, we postpone the detailed derivations to the Appendix A. Combining the obtained perturbation potentials with the basic potentials (8), we finally derive the total potentials from (17) in the form:

1 2

γ (w ) = Cσ , 1 rφ (w ) + r

N 

αk w−k ,

(19)

k=1



 r ψ (w ) =  C w+ βk w−k φ (w ) σ , 2 N+2 k=1

 −

1 Cσ , 1 h¯ 2

(20)

72

W.-N. Zou, Q.-C. He / International Journal of Engineering Science 126 (2018) 68–96

with N+2 

βk w

−k

k=1

2μU¯ 0  = 1 − φ  (w ) − C¯σ∗ ,1 w−1 φ  (w ) + r

 w

−1

h¯ + r



N 

kαk w−k−1 +

k=1

N 

λk w−k .

(21)

k=1

In the above solution, the rigid-body rotation α and coefficients {α1 , . . . , αN } are solved from (117) and (136) by an iterative procedure, and then the parameter U0 and {λ1 , . . . , λN } are explicitly given by (118) and (123), respectively. Back to the special case of a cavity, the iterative procedure is confined only to the coefficients {α1 , . . . , αN }, since the rigid-body p rotation vanishes constantly, and there is an identity βN+2 = αN + Cσ ,1 bN ≡ 0 because of κ = −1. These results imply that the perturbation potentials of a rigid inclusion or a cavity in an infinite elastic body subjected to a uniform remote loading have finite representations when the shape of the excluded subdomain is characterized by a Laurent polynomial with finite terms. Finally, the displacement, stress and strain fields are calculated through (2)–(4). 3. Some comments 3.1. Consideration of finite rigid-body rotation Softhard materials integration has broad applications in a large number of situations of practical interest. Large rigid-body rotation of hard phases in company with the large deformation of soft matrix poses a challenge to the theory of elasticity. As an easy part of it, we consider the effect of finite rigid-body rotation first, using a rigid inclusion in the shape of sharp rhombus as an example. We will express the elastic fields in the coordinate system associated with the deformed matrix. When the rigid inclusion ω undergoes a translation U0 and a rotation α due to the remote loading, the deformed matrix can be defined by

z = H + reια φ (w ), H = h + U0 , |w| ≥ 1

(22)

if  = {z = h + rφ (w ), |w| ≥ 1} represents the matrix without deformation. Thus, the rigid-body displacement is described by





Urigid = 1 − e−ια (z − H ) + U0 = (eια − 1 )rφ (η ) + U0 ,

|η| = 1,

(23)

and the potentials can be derived in the form

1 2

γ (w ) = Cσ , 1 reια φ (w ) + reια

N 

αk w−k ,

(24)

k=1



 re−ια ψ (w ) =  e2ια Cσ , 2 w + βk w−k φ (w ) N+2



k=1

−

1 Cσ , 1 H¯ , 2

(25)

where the translation U0 , the angle α of the rigid-body rotation and the coefficients {α k , β k } are given in Appendix A2. For the case where the inclusion is elliptical, say φ (w ) = w + b1 w−1 , the rigid-body rotation can be obtained explicitly as

α=

(1 + κ )Im(Cσ ,2 b1 )   2μ κ + b1 b¯ 1

(26)

for the small rigid-body rotation, or implicitly as

  (1 + κ )Im e2ια Cσ ,2 b1 sin α =   2μ κ + b1 b¯ 1

(27)

for α being finite. These results coincide with those given by Muskhelishvili (1953) and Yang et al. (2017), respectively. Two points must be mentioned herein: first, the so-called explicit solution for the problem of an arbitrary rigid inclusion obtained by Yang et al. (2017) is quite limited since the analogical replacement of (14) by (15) in their paper is unsuitable for shapes other than ellipse; second, the above treatment for finite rotation is not easily extended to problem of inhomogeneity other than rigid inclusion. 3.2. Solutions for elliptical inhomogeneities, including rigid inclusions and cavities The global displacement of a rigid inclusion or a cavity relative to the matrix, as used in the above section, can be extended in general to Eshelby’s inclusion/inhomogeneity problems (Eshelby 1957, 1959). The additional rigid-body displacement of the inclusion or the inhomogeneity is also relative and must meet the requirement of self-consistency. The equations for Eshelby’s inhomogeneity problem in two-dimensional space are formulated with the K–M potentials as follows.

W.-N. Zou, Q.-C. He / International Journal of Engineering Science 126 (2018) 68–96

73

Using the subscript ‘M’ for the matrix and ‘I’ for the inhomogeneity, and employing the superscript ‘+’ for the fields in the inhomogeneity and ‘−’ for the fields out of the inhomogeneity, if the basic potentials are chosen to be

1 2

γo (z ) = Cσ ,1 (z − h ) +

2 μI C (z − h )χ ω , κI − 1 ε,1

(28)

1 2

ψo (z ) = Cσ ,2 (z − h ) − Cσ ,1 h + 2μICε,2 (z − h )χ ω −

2 μI C hχ ω , κI − 1 ε,1

(29)

where χ ω is the characteristic function of the subdomain ω, namely χ ω = 1 if z ∈ ω and 0 if z ∈ ω, then the continuity conditions of displacement and traction on the interface can be expressed with the perturbation potentials by

  κM γ p− (t ) − t γ˙ p− (t ) − ψ p− (t ) κI γ p+ (t ) − t γ˙ p+ (t ) − ψ p+ (t ) = + C1 rφ (η ) + C¯2 r φ¯ η−1 + 2[U0 + ια rφ (η )], μM μI

(30)

γ p− (t ) + t γ˙ p− (t ) + ψ p− (t ) = γ p+ (t ) + t γ˙ p+ (t ) + ψ p+ (t ).

(31)

Above, t is indicated by (18) but w → η with |η| = 1, the matrix is fixed through the basic potentials (8), the parameters C1 and C2 are given by the remote loadings Cσ , 1 and Cσ , 2 and/or the eigenstrains Cε, 1 and Cε, 2 by

D2 Cσ ,1 , C2 = −2Cε,2 + D1Cσ ,2 , 2

C1 = 2Cε,1 − where

D1 =

1

−

μM

1

μI

, D2 =

(32)

κM − 1 κI − 1 − . μM μI

(33)

Let us consider an elliptical shape: φ (η ) = η + b1 η−1 . It is well-known that the interior perturbation elastic stresses/strains of an elliptical inhomogeneity due to uniform remote loadings or uniform eigenstrains are uniform, say

γ p+ (z ) = s+1 (z − h ), ψ p+ (z ) = s+2 (z − h ) − s+1 h¯ , s+ 1

(34)

s+ . 2

with a real coefficient and a complex coefficient Substituting (31) and (34) into (30), and posing C1α = C1 + 2ια , we have

 κ + 1       κM + 1 − M α r η + b η−1 + D s¯+ + C¯ r η−1 + b ¯ 1 η + 2U0 γ p (η ) = − D2 s+ + C 1 1 2 2 1 1 μM μM

resulting in U0 = 0,

γ ( w ) = μM − p

and

 κM +1

(35)



+ + α ¯ μM − D2 s1 b1 + C1 b1 + D1 s¯2 + C2

κM + 1

−1 rw−1 = s− 1 rw

κ + 1   +  M α ¯ ¯ ¯ − D2 s+ 1 + C1 + D1 s2 + C2 b1 = 0. μM

(36)

(37)

Substituting (36) into (31) yields

  −1   h¯ + rη−1 + r b¯ 1 η − −2 + −1 r s1 η − r s¯− + b¯ 1 η + r s+ 2 η + b1 η 1 η + 2r s1 η 1 − b1 η−2

ψ p− (η ) = such that



ψ (w ) = − p



h¯ w−2 + 1 + b1 b¯ 1 rw−3 1 − b1 w−2





−¯ + + −1 s− 1 + s1 b1 + 2s1 + s2 b1 rw

(38)

and + ¯ ¯− 2s+ 1 b1 + s2 − s1 = 0.

(39)

Solving the algebraic equations (37) and (39), s+ and s+ can be worked out as follows 1 2

1









M +1 1 − b1 b¯ 1 D1C¯1α − κμ C¯ α + C2 b1  κM +1  M 1 κI +1 κM +1 , ¯ μM + D 2 1 − b 1 b 1 D 1 + μI μM

s+ 1

=

s+ 2

κM +1 κM +1 ¯ ¯α ¯ μM 2 C1 b1 + C2 − μM + D2 1 − b1 b1 C2 = .  κM +1   κI +1 κM +1 ¯ μM μM + D 2 1 − b 1 b 1 D 1 + μI μM

μM 1









(40)



(41)

Since the coefficient s+ is assumed to be a real, we get 1

α=

κM +1 ImC1α 1 μ Im(C2 b1 ) = κ +1 M  , M 2 2 − D1 1 − b1 b¯ 1 μM

(42)

74

W.-N. Zou, Q.-C. He / International Journal of Engineering Science 126 (2018) 68–96

and so

μM γ p− (w ) 1 − b1 b¯ 1

=

 κM +1 μM

μM ψ p− (w ) 1 − b1 b¯ 1

=

 κM +1

 M +1  D1C1α b1 + κμ − D2 C¯2 M + D2



rw  I +1 κM +1 1 − b1 b¯ 1 D1 + κμ μM I

 M +1  D1C1α b1 + κμ − D2 C¯2 M

μM + D 2



 

,

(43)



κM +1



h¯ w−2 + 1 + b1 b¯ 1 rw−3

 I +1 κM +1 1 − b1 b¯ 1 D1 + κμ μM I 

−1

1 − b1 w−2





D1 1 + b1 b¯ 1 − μM 2C1 − D2 C2 b1 + C¯2 b¯ 1 + rw−1 .  κM +1   ¯ 1 D1 + κI +1 κM +1 b + D 1 − b 2 1 μM μI μM

(44)

From (135) and (38), together with (39), it is easy to verify that the resultant moment of exterior potentials on the boundary vanishes constantly:



 −  −1 + ψ p− (z )dz = r2 Re s1 b¯ 1 + 2s+ 1 + s2 b1 η d η ∂ω |η|=1   +¯ + + ¯ ¯ b b = −2π r 2 Im 2s+ b + s 1 1 1 2 1 + 2s1 + s2 b1 ≡ 0

Mo = Re

(45)

s+ 1

since the coefficient is real. For the inclusion problem, κI = κM , μI = μM , the exterior perturbation potentials become

γ p− (w ) = −

1 − b1 b¯ 1 2μMC¯ε,2 rw−1 , κM + 1

1 − b1 b¯ 1 ψ p− (w ) = − 2 μM κM + 1



(46)







h¯ w−2 + 1 + b1 b¯ 1 rw−3 1 − b1 w−2

C¯ε,2 + 4Cε,1 rw−1 ,

(47)

which coincide with the results derived for the inclusion problem by Zou and Lee (2017) in a different way, except that the rigid-body rotation is hidden in the image part of the linear coefficient of the perturbation potential γ p+ in the previous consideration. When the inhomogeneity is elliptical and rigid, the perturbation potentials are obtained directly from (12) as

γ p− (w ) =

Cσα,1 b1 + C¯σ ,2

κM

rw−1 ,

(48)



ψ (w ) = − d



Cσα,1 b1 + C¯σ ,2 h¯ w−2 + 1 + b1 b¯ 1 rw−3

κM

1 − b1 w−2



  κM − 1  κM − b1 b¯ 1 − Cσ ,1 + Re Cσ ,2 b1 rw−1 . κM κM

(49)

For the cases of elliptical cavities, the perturbation potentials can be obtained from the rigid inclusion solutions by setting

κM = −1 and Cσα,1 → Cσp,1 , namely

  γ p− (w ) = − Cσp,1 b1 + C¯σ ,2 rw−1 , ψ

− d



p

(w ) = − Cσ ,1 b1 + C¯σ ,2

(50)

   h¯ w−2 + 1 + b1 b¯ 1 rw−3 1 − b1 w−2

−







1 + b1 b¯ 1 Cσp,1 + 2Re Cσ ,2 b1



rw−1 .

(51)

It is easy to verify that the solutions of the elliptical rigid inclusion and cavity problems cannot be deduced from (43) and (44) by any parameter reduction, say μI → ∞, since Eqs. (30) and (31) are not equivalent to (12) or (15). 3.3. Two kinds of correspondence According to the equivalent eigenstrain idea (Eshelby, 1957), if the interior stresses are uniform, there is a correspondence between the exterior perturbation fields of an inhomogeneity under remote loadings and those of an inclusion undergoing eigenstrains. Now we have proved that such a correspondence exists only when the shape is ellipse (Kang and Milton, 2008; Liu, 2008). From the derivation of the last subsection, we see that in the correspondence the eigenstrains can be calculated from the remote stresses by









D2 M +1 C − 2αι b¯ 1 + κμ − D2 Cσ ,2 2μMCε,2 2 σ ,1 M = −  κ +1 ,   M (κM + 1 )D1 + D2 1 − b1 b¯ 1 D1 + κI +1 κM +1

μM



μI



μM

(52)



κM +1 ¯ 8μMCε,1 μ Cσ ,1 + D1 2Re(Cσ ,2 b1 ) − 1 + b1 b1 Cσ ,1 = − M  κ +1   κI +1 κM +1 M ¯ (κM + 1 )D2 μM + D 2 1 − b 1 b 1 D 1 + μI μM

(53)

W.-N. Zou, Q.-C. He / International Journal of Engineering Science 126 (2018) 68–96

75

with

α=

κM +1 D1 Im(Cσ ,2 b1 ) 2 μM



κM +1 ¯ μM − D 1 1 − b 1 b 1



(54)

to carry out the same perturbation fields of inhomogeneity problem from the inclusion problem. Detailed investigation shows that the perturbation potentials of elliptical rigid inclusions and cavities have the same function form as the perturbation potentials of elliptical inclusions with eigenstrains, which implies that there also exist equivalent eigenstrains for the rigid-inclusion and cavity problems. Certainly, it comes true only for elliptical shapes. The correspondence of the exterior perturbation fields cannot be extended to the shapes other than ellipse. The unified formulation (17) of the cavity and rigid inclusion problems for arbitrary shapes presents another valuable correspondence that mainly focuses on the function form instead of the function value. This kind of study can be traced back to the pioneering work of Dundurs (1989). From the above derivations, we may give two comments: (1) there is a mathematical correspondence between the rigid inclusion problem and the cavity problem by setting κM = −1 when their shapes are characterized by the same type of Laurent polynomial; (2) there is no physical correspondence, say the same perturbation stresses, since κM = −1 is physically unreal, though in the case of ellipse the perturbation potentials of both the rigid inclusion and the cavity have counterparts of inclusions with proper eigenstrains. The first kind of respondence is applicable to the 3D ellipsoidal problems. But the second kind of correspondence is in general unsuitable for the 3D elasticity problems of cavities and rigid inclusions. Markenscoff and Paukshto (1995) tried to extend such a correspondence to three-dimensional elasticity, and obtained an additional condition that the dilatation must be constant or linear function of position. For the same shape, the solution for a rigid inclusion and a cavity have the same polynomial structure, and the latter could be regarded as a subset of the former in some way. The viewpoints that may not be clarified before for the cavity problem include: (1) the rigid-body translation correspondence in derivation by 2μU0 → f0 and the real rigid-body transκf lation U0 = − 2μ0 in calculation, (2) the null rigid-body rotation, (3) the always vanishing of the coefficient βN+2 , (4) the degeneration of hydrostatic pressure. Then, we can focus on the study of the rigid inclusion problem instead of the cavity one. Finally, from the derivation of Appendix A, we see that the rigid-body translation is just an extraction as (118) but the rigid-body rotation is deeply involved in the perturbation potential. 3.4. Transformation with the coordinate system In the general derivations of perturbation potentials, we purposely maintain the position parameter h and the size parameter r of the subdomain. From the obtained results, we find that bk , ck and gk , k = 1, . . . , N, are independent of h and r, and so the coefficients α k of the perturbation potential γ p and the rigid-body rotation α are irrelevant to these two parameters. On the other hand, the coefficients β k seem to be apparently connected with the position h, but careful investigation shows that the terms involving h identically disappear in the combination zγ  (z ) + ψ (z ), implying that they actually have nothing to do with the elastic fields and behave just as the compensation of h occurring in the complex variable z. Therefore, with no loss of generality, the position of the inclusion can be set to be the origin of the physical plane, namely h = 0. Beside this property, all elastic fields are directly proportional to r or the size of the subdomain. Almost all involved variables will change their values when the coordinate system undergoes a rotation of angle ϕ . In order to use the original conformal mapping, the new mapping between the complex variable Z = (z − h )e−ιϕ in the physical plane and the complex variable W = we−ιϕ in the image plane can be written by



Z = rφ (w )e

−ιϕ

∞  =r W+ BnW −n



= r(W ),

(55)

n=1

with

Bn = bn eι(n+1 )ϕ .

(56)

Further considering the transformation of displacement and stress, we have the expressions of the potentials in new coordinate system as

(W ) = eιϕ γ (w ), (W ) = e−ιϕ ψ (w ),

(57)

resulting in the perturbation potentials in the new coordinate system written by

p = r

N k=1

AkW −k ,  p =

2 N+1 r BkW −k .  (W ) k=1

Thus, we obtain the transformation relations

Ak = αk eι(k+1 )ϕ , Bk = βk eι(k−1 )ϕ due to  (W ) =

φ  (w ).

(58)

76

W.-N. Zou, Q.-C. He / International Journal of Engineering Science 126 (2018) 68–96

4. Examples and illustration In this section, some reported explicit solutions are collected in order to be compared with the corresponding solutions obtained by the general method of the present work. Then, two typical examples for a rhombus rigid inclusion and a pentagram cavity are analyzed and illustrated. For the convenience of comparison, the complementary terms, say the rigid-body displacement of rigid inclusion or the hydrostatic pressure and the global translation of cavity, are added into the basic potentials, namely using

γo =

1 2

Cσ , 1 −



2μαι (z − h ), κ +1

1 2

ψo = Cσ , 2 (z − h ) − Cσ , 1 h¯ + 2μU¯

(59)

for the rigid inclusion problem, and

1 2

1 2

γo = (Cσ , 1 + P )(z − h ), ψo = Cσ , 2 (z − h ) − Cσ , 1 h¯ + f¯0

(60)

for the cavity problem. Then we can verify the results by checking the displacement condition

κγ (η ) −

z (η ) z (η )

γ  (η ) − ψ (η ) = 0, |η| = 1

(61)

on the boundary of the rigid inclusion, or the traction condition

γ (η ) +

z (η ) z (η )

γ  (η ) + ψ (η ) = 0, |η| = 1

(62)

on the boundary of the cavity. 4.1. Some examples for polygons approximated For a polygon, the Schwarz–Christoffel formula (Appendix B) provides a precise description with all parameters exactly determined. When its series expansion is truncated to obtain an approximate description, its area becomes a little larger: the less the terms are, the larger it is. In order to make the approximation more efficiently, we introduce the area equivalence in which the size parameter r is calculated according to the area of the approximate shape. Assuming that |ω0 | is the area of a polygon and N is the highest index in the truncation, the approximate size parameter ra is calculated by



ra =

|ω0 |   < r.  π 1 − Nk=1 kbk b¯ k

(63)

Concerning a rectangle, a rhombus, a pentagram and a regular polygon, their recursive formulae for the shape parameters including bk are provided in Appendix B. In all examples of this subsection, the remote loading is a uniaxial tension or compression along with the x-direction, 1 ∞ 1 ∞ = 1 σ, C ια namely Cσ ,1 = 12 σ11 σ ,2 = − 2 σ11 = − 2 σ . If the body changes its orientation, say z = h + re φ , it suffices to change 2 1 2 ια the coefficient Cσ , 2 by setting Cσ ,2 = − 2 σ e . In the first two examples, the body is in plane stress with Poisson ratio ν = 0.49 like (Misseroni, Dal Corso, Shahzad, & Bigoni, 2014). Other examples in this subsection are taken from Savin (1961) for various cavities. (1) The first example concerns a rectangular rigid inclusion with aspect ratio ly /lx = 0.25, which is approximated by taking a truncation with 8 terms. The shape expression (150) is the same as that given by Misseroni et al. (2014), while the area-equivalent size parameter r/lx = 0.3538 is a little smaller than the theoretical value. But the potentials are quite different even though the equivalence should be established absolutely between the interior-to-exterior mapping and the exterior-to-exterior mapping by the transformation w → w−1 . We also check the balance relation (61) with the potentials given by Misseroni et al. (2014). The residuals compared with those calculated from the new derived potentials listed in the following are shown in Fig. 1(a), which shows that the previous results cannot well satisfy the boundary relation. In our results, the denominator polynomial in the perturbation potential ψ p is taken to be φ  and the coefficients α i and β i in (129) and (130) for this situation are such that

γ (w ) =

ψ (w ) =

σ



z (w ) + rσ −0.5453w−1 + 1.464 × 10−3 w−3 + 3.605 × 10−3 w−5 + 3.345 × 10−3 w−7 + 1.480 × 10−3 w−9

4  − 1.327 × 10−4 w−11 − 5.139 × 10−4 w−13 + 8.710 × 10−5 w−15 ,

(64)



1 rσ − w + 0.3649w−1 − 0.9546w−3 − 5.797 × 10−2 w−5 + 3.821 × 10−2 w−7 + 7.429 × 10−2 w−9 φ  (w ) 2



+ 4.917 × 10−2 w−11 + 3.192 × 10−3 w−13 − 2.002 × 10−2 w−15 + 3.507 × 10−3 w−17 .

(65)

W.-N. Zou, Q.-C. He / International Journal of Engineering Science 126 (2018) 68–96

77

Fig. 1. Displacement residual on the boundary of the rigid inclusion: (a) the present potentials (64) and (65) vs. the Misseroni et al. potentials for the rectangular rigid inclusion, (b) the present potentials (66) and (67) vs. the Misseroni et al. potentials for the rhombus rigid inclusion.

(2) The second example is about a rhombus rigid inclusion with axis ratio ly /lx = 2/15 approximated by a truncation with 8 terms. The shape expression (160) is also identical to that in Misseroni et al. (2014), and the area-equivalent size parameter r/lx = 0.2659 is almost the same as the theoretical value. The comparison of the residuals is similar to the first example, as shown in Fig. 1(b). Our results for such a rhombus rigid inclusion are

γ (w ) =

ψ (w ) =



σ

z (w ) + rσ −0.5087w−1 − 8.106 × 10−3 w−3 + 3.117 × 10−3 w−5 − 1.167 × 10−3 w−7 + 1.102 × 10−3 w−9 4  − 4.193 × 10−4 w−11 + 5.676 × 10−4 w−13 − 1.753 × 10−4 w−15 , (66)



1 rσ − w + 0.3611w−1 − 0.6214w−3 − 7.639 × 10−2 w−5 + 5.075 × 10−2 w−7 − 2.787 × 10−2 w−9 φ  (w ) 2



+ 3.030 × 10−2 w−11 − 1.659 × 10−2 w−13 + 2.178 × 10−2 w−15 − 7.058 × 10−3 w−17 .



(67)



(3) For a triangle cavity characterized by z (w ) = r w + 13 w−2 , the reported results (Kachanov, Shafiro, & Tsukrov, 2003; Savin, 1961) can be recast in the form

γ (w ) =





1 rσ w + 2w−1 − w−2 , 4 3

ψ (w ) = −

rσ 2



w−



3 − 11w−1 + 9w−3 . 9 − 6w−3

(68)

According to our formulae, the simple parameters

g1 = b 1 = 0 , g2 = b 2 =

1 3

(69)

in this case lead to the explicit solutions

α1 = −c1 = −C¯σ ,2 = and

σ 2

,

σ rσ α2 = −c2 = −Cσp,1 b2 = − , f0 = rα¯ 1 g2 = 6

6

(70)

    σ λ1 = α1 b¯ 1 + b1 α¯ 1 + Cσp,1 b1 b¯ 1 + 2 α2 b¯ 2 + b2 α¯ 2 + Cσp,1 b2 b¯ 2 = − ,

(71)

  σ λ2 = 2 α2 b¯ 1 + b2 α¯ 1 + Cσp,1 b2 b¯ 1 = 2b2 α¯ 1 = .

(72)

9

3

Then, the potentials are obtained to be

⎧ ⎪ ⎪ ⎨

 1 1 1 γ (w ) = rσ w + w−1 − w−2 , 4 2 12   −1 1 1 11 w − w−3 ⎪ ⎪ . ⎩ψ (w ) = −rσ 2 w − 6 + 18 φ  w ( )

(73)

78

W.-N. Zou, Q.-C. He / International Journal of Engineering Science 126 (2018) 68–96

This is the same as (68). But the constant term from f0 is hidden in (68)2 and incompatible with the expansion reported by Chang and Conway (1968), and Savin (1961) added a constant terms in the expansion of perturbation potential ψ d , say in (II.5)2 , to contain it, which is somehow unusual.   1 If one more coefficient is involved so that z (w ) = r w + 13 w−2 + 45 w−5 , Savin (1961) reported the potentials (by setting α = 0 in equation (II.59)) as

γ (w ) = r σ

1 4

w+

 ψ (w ) = −rσ



1 −2 1 −5 675 −1 15 w − w + w−3 − w , 1348 12 1348 180

(74)

31004w−4 − 546675w−5 − 31725w−6 − 1 909989w−1 − 19035w−2 − 101475w−3     w+ + 2 181980 9 − 6w−3 − w−6 181980 9 − 6w−3 − w−6

112 −7 w 3

 . (75)

Using the non-zero parameters in this case

g5 = b 5 =

1 47 , g2 = b2 + b¯ 1 g4 + 2b¯ 2 g5 = , 45 135

c1 = C¯σ ,2 + Cσp,1 b1 = −

(76)

1 1 1 σ , c2 = Cσp,1 b2 = σ , c5 = Cσp,1 b5 = σ , 2 6 90

(77)

we have algebraic equations

α5 = −c5 = −

1 σ, 90

α3 = −c3 + g5 α¯ 1 = 1 6

α2 = −c2 + g4 α¯ 1 + 2g5 α¯ 2 = − σ +

1 α¯ 1 , 45

2 α¯ 2 , 45 1 2

α1 = −c1 + g3 α¯ 1 + 2g4 α¯ 2 + 3g5 α¯ 3 = σ +

1 α¯ 3 15

to yield

α5 = −

1 σ, 90

α3 =

15 σ, 1348

α2 = −

15 σ, 86

α1 =

675 σ 1348

(78)

and

f0 = r α¯ 1 g2 =

235 rσ . 1348

(79)

The determination of the additional coefficients

    4273 1 λ1 = 2 α2 b¯ 2 + b2 α¯ 2 + σ b2 b¯ 2 + 5 α5 b¯ 5 + b5 α¯ 5 + Cσp,1 b5 b¯ 5 = − σ, λ2 λ4

2 34830 465 5 ¯ = 2b2 α¯ 1 + 3α3 b2 = σ , λ3 = 5b5 α¯ 3 = σ, 1348 4044   5 1 75 = 5 α5 b¯ 2 + b5 α¯ 2 + σ b5 b¯ 2 = − σ , λ5 = 5b5 α¯ 1 = σ. 2 258 1348

finally yields the potentials

1

(80)



1 −5 675 −1 15 −2 15 w − w + w−3 − w 4 1348 86 1348 90  1 47 −2 1 −5 675 −1 15 = rσ w+ w − w + w−3 − w , 4 1348 516 1348 180

γ (w ) = r σ

φ (w ) +

235rσ ψ (w ) = − 1348 =

rσ φ  (w )

235rσ − rσ 1348





(81)

1395w−2 625w−3 30w−5 235w−6 w 21688w−1 3w−4 + − − + − − 2 34830 4044 1011 86 337 12132

w + 2

21688w−1 34830

−

47w−2 4044

−

625w−3 1011

−4

+ 3w − 86  φ (w )

203w−5 6066

−

235w−6 12132





.

(82)

By comparison, we find that, besides many differences in the coefficients of the polynomials, there are structural flaws in (75): lack of constant term and impossible highest-negative-power term, and consequently the old result cannot satisfy the balance condition (62), as shown in Fig. 2(a).

W.-N. Zou, Q.-C. He / International Journal of Engineering Science 126 (2018) 68–96



(4) For a square cavity described by z (w ) = r w − 16 w ishes. Our result



1

3 1 −3 w + w−1 + w , 4 7 24

γ (w ) = r σ

79

 −3

, due to the symmetry, the relative rigid-body translation van-



ψ (w ) = −rσ

1 91w−1 − 78w−3 w+ 2 84(2 + w−4 )

 (83)

is the same as (II.35) in Savin (1961); there is a missing negative sign in the potential ψ in Kachanov et al. (2003). When the number of the terms in z (w ) increases, the results given by Savin (1961) deviate from ours more and more, and should be incorrect because they do not satisfy (62). We believe that these mistakes, as shown in the following,   should come mainly 1 from the informal derivation. The reported solutions for the case where z (w ) = r w − 16 w−3 + 56 w−7 are

1

γ (w ) = r σ

4



w + 0.426w−1 + 0.046w−3 + 0.008w−5 − 0.004w−7 ,

 ψ (w ) = −rσ

1 0.548w−1 − 0.457w−3 − 0.026w−5 − 0.029w−7 w+ 2 φ  (w )

(84)

 (85)

from (II.41) while those of us are

1

γ (w ) = r σ

4



w + 0.4259w−1 + 4.638 × 10−2 w−3 + 7.605 × 10−3 w−5 − 4.464 × 10−3 w−7 ,

 ψ (w ) = −rσ

1 w+ 2

(86)



 1  0.5475w−1 − 0.4576w−3 − 2.516 × 10−2 w−5 − 2.876 × 10−2 w−7 . φ  (w ) 

The reported solutions for z (w ) = r w − 16 w−3 +

1 −7 56 w

−

1 −11 176 w



(87)

(Note: Savin (1961) presented only the solution for

α = π /2 in (II.44), the solution for α = 0 can be achieved by changing the signs of terms associated with w−1 , w−5 , w−9 in γ (w ) and terms associated with w, w−3 , w−7 , w−11 in ψ (w ). During the comparison, there could be typos for the signs of the terms associated with w−3 , w−7 , w−11 in ψ (w ), so we show two solutions in Fig. 2(c, d) for demonstration): the old potentials are

1

γ (w ) = r σ

4



w + 0.4254w−1 + 0.0476w−3 + 0.0086w−5 − 0.0060w−7 − 0.0024w−9 + 0.0014w−11 )

1 2

ψ (w ) = − r σ w −

(88)

 rσ  0.479w−1 − 0.457w−3 + 0.269w−5 − 0.037w−7 + 0.073w−9 − 0.017w−11 − 0.031w−13 ). φ  (w )

(89)

whereas the newly obtained ones are

γ (w ) = r σ

1 4

w + 0.4254w−1 + 4.763 × 10−2 w−3 + 8.562 × 10−3 w−5

− 5.986 × 10−3 w−7 − 2.417 × 10−3 w−9 + 1.421 × 10−3 w−11 1 2

ψ (w ) = − r σ w −





(90)

 rσ  0.5493w−1 − 0.4564w−3 − 3.147 × 10−2 w−5 φ  (w )  −2 −7 −2 −9 −2 −11

− 3.658 × 10

w

+ 1.624 × 10

w

+ 1.709 × 10

w

.

(91)

We find that the first potentials γ (w ) given by (90) and (88) are almost the same, but the second potentials ψ (w ) are obviously different, even though the possible typos have been corrected. In (89), the impossible highest-negative-power term happens again.   (5) For a hypocycloid cavity characterized by z (w ) = r w + bN w−N with N > 1, due to the symmetry, there are

gN = bN , c1 = C¯σ ,2 , cN = Cσp,1 bN .

(92)

Then, we obtain

α2 = −c2 = −Cσp,1 b2 ; α1 = −c1 = −C¯σ ,2

(93)

for N = 2,

α3 = −c3 = −Cσp,1 b3 , α1 = −c1 + g3 α¯ 1 = −C¯σ ,2 + b3 α¯ 1 or α1 = −

C¯σ ,2 +Cσ ,2 b3 1−b b¯ 3 3

for N = 3, and the equations for N > 3

αN = −cN = −Cσp,1 bN ; αN−2 = gN α¯ 1 , α1 = −c1 + (N − 2 )gN α¯ N−2 ;

(94)

80

W.-N. Zou, Q.-C. He / International Journal of Engineering Science 126 (2018) 68–96

Fig. z (w ) =  2. Traction1 residual  on the boundary of the cavity: (a) the present potentials (81) and (82) vs. the Savin potentials for thetriangular cavity 1 w−5 , (b) the present potentials (86) and (87) vs. the Savin potentials (84) and (85) for square cavity z (w ) = r w − 16 w−3 + 56 w−7 , (c) r w + 13 w−2 + 45





1 1 the present potentials (90) and (91) vs. the Savin potentials (88) and (89) for square cavity z (w ) = r w − 16 w−3 + 56 w−7 − 176 w−11 , (d) the same as the former but without considering the typos. Fig. 2(b) and Fig. 2(c) are similar to each other, and for Savin’s solutions Fig. 2(d) without considering the typos have bigger residuals than Fig. 2(c) from the potentials (88) and (89) correcting the typos.

αN−k−1 = kgN α¯ k , αk = (N − k − 1 )gN α¯ N−k−1 ,

N − 1 2

≥k>1

(95)

resulting in three non-zero coefficients

αN = −Cσp,1 bN ; αN−2 = −

Cσ ,2 bN 1 − (N − 2 )bN b¯ N

,

α1 = −

C¯σ ,2 1 − (N − 2 )bN b¯ N

.

(96)

The constant term f0 = r (N − 1 )α¯ N−1 gN always vanishes except that

f0 = r α¯ 1 b2 = −rCσ ,2 b2

(97)

for N = 2. The derivation of the additional coefficients

  λ1 = N αN b¯ N + bN α¯ N + Cσp,1 bN b¯ N = −NCσp,1 bN b¯ N , λ3 = NbN α¯ N−2 = −

yields the potentials

NC¯σ ,2 bN b¯ N 1 − (N − 2 )bN b¯ N



,

λN = NbN α¯ 1 = −

N bNCσ ,2 1 − (N − 2 )bN b¯ N

(98)



r C¯σ ,2 w−1 + Cσ ,2 bN w2−N 1 γ (w ) = Cσ ,1 z(w ) + f0 − − Cσp,1 r bN w−N , 2 1 − (N − 2 )bN b¯ N

(99)

W.-N. Zou, Q.-C. He / International Journal of Engineering Science 126 (2018) 68–96

81

Fig. 3. A rhombus rigid inclusion in the matrix, with Young module E = 22 MPa and Poisson ratio 0.49, under the remote loading σ along with an angle inclined to the x-axis. The in-plane principal stress difference on the doted line is shown in Fig. 5.

Fig. 4. Stresses along the boundary of a rhombus rigid inclusion, where the remote loading σ = 0.25 MPa is parallel to the x-axis and the boundary is indicated by z(η ) = z(eιθ ).



ψ (w ) =



 r N bN f¯0 −N−1  C¯σ ,2 w−3 Cσ ,2 w + w − 1 + N bN b¯ N Cσp,1 w−1 +  φ (w ) r 1 − (N − 2 )bN b¯ N

 −

 (2N − 2 )Cσ ,2 bN −N w . 1 − (N − 2 )bN b¯ N (100)

1 3,

These potentials coincide with: (i) the formulae (68) by setting N = 2 and b2 = and (83) by posing N = 3 and b3 = − 61 , p 1 1 under uniaxial tension Cσ ,1 = 2 σ , Cσ ,2 = − 2 σ ; (ii) the result for polygonal hole in page 66 of Kachanov et al. (2003) by setting bN = m, Cσ ,1 = 12 p and Cσ ,2 = − 12 pe−2ια , (iii) the formulae (3.10) and (3.18) in Ekneligoda and Zimmerman (2008) by setting bN = m1 , Cσ ,1 = 0 and Cσ ,2 = ιτ , where the remote shear loading is transformed into the force on the boundary of the pore. In summary, many known potentials for plane rigid inclusion and cavity problems can be treated as special cases of the results presented in the present work. In addition, in light of our results, many mistakes in the reported results can be found, especially when the shape of the rigid inclusion or the cavity involves terms more than three. p

4.2. Analysis of the solution of a rhombus rigid inclusion Here we report some interest features of the elastic fields given by the solution of the second example in the last subsection, wherein the elastic matrix is in plane stress, with Young module E = 22 MPa and Poisson ratio equal to 0.49. We use the Laurent polynomial (160) to approximate the rhombus as shown in Fig. 3, whereby the curvature radius at the tip reaches 0.1136μm while the size parameter r = 7.976 mm. When the remote loading parallels the x-axis, namely ϕ = 0, there is no rigid-body displacement of the rigid inclusion relative to the matrix, but we find that the potentials given by

82

W.-N. Zou, Q.-C. He / International Journal of Engineering Science 126 (2018) 68–96

Fig. 5. Toriodal stress and surface traction at the tip for different loading angles ϕ , where the remote loading σ = 0.25 MPa.

Fig. 6. Rigid-body rotation of the rigid inclusion, where two levels of remote loading are considered, namely σ = 0.25 MPa and σ = 2.5 MPa. The solutions based on the small deformation exhibit the symmetry with respect to the loading angle ϕ while the symmetry breakup emerges when considering the effect of finite rigid rotation, especially for the case of large loading.

Misseroni et al. (2014) cannot satisfy the displacement condition on the boundary, as shown in Fig. 1(b), and should be replaced by (66) and (67). We calculate the toroidal stress σ θ along boundary and the surface traction on the rigid inclusion, with the normal component σ n and the tangent component τ n , and illustrate them in Fig. 4. The prominent concentration around the tip is very serious, while the biggest concentration of toroidal stress is about 168σ . When the loading direction deviates from the x-axis, the stress on the boundary changes. From Fig. 5, we find that the normal stresses at the tip weaken with the loading angle ϕ increasing, until they switch from stretch to compression at the loading angle about 55°. An interesting property is the symmetry of shear stress τ with respect to the loading angle ϕ (Fig. 5), though there is no such a geometric symmetry actually. We further check and find the same symmetry of the rigidbody rotation α of the rigid inclusion (Fig. 6), and its breakup when the effect of finite rigid-body rotation is considered by using the solution presented in Section 3.1. From Fig. 6, we also see that the breakup of symmetry becomes much obvious with the remote loading increasing. This phenomenon has been investigated in the solutions of elliptical rigid inclusions

W.-N. Zou, Q.-C. He / International Journal of Engineering Science 126 (2018) 68–96

83

Fig. 7. In-plane principal stress difference σ on a line, as shown in Fig. 3, coincides with the experimental data, where the remote loading σ = 0.25 MPa parallel to the x-axis.

Fig. 8. A pentagram cavity in an infinite space, where the dashed circle and colored dot define the domains whose fields are illustrated in Fig. 10. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article)..

(Yang et al., 2017). Finally, in Fig. 7, we check the in-plane principal stress difference, defined by σ = 2|z¯γ  (z ) + ψ  (z )| , on a line that coincides well with the experimental data provided by Misseroni et al. (2014). We have checked carefully the results for the rhombus rigid inclusion with lx /ly = 7.5 as shown in Fig. 5, and found that the strain at the tip could not be small while the rigid-body rotation (or the remote loading) becomes large. For example, the minimal tip stress with respect to the loading angle ϕ is (σn , σθ , στ ) = (0.355, 0.725, −8.71 )σ with σ = 0.25 MPa at 11 ϕ = 36 π , and we obtain the rigid-body rotation angle α = 0.432◦ = 0.007537rad, not greater than the maximal shear strain

84

W.-N. Zou, Q.-C. He / International Journal of Engineering Science 126 (2018) 68–96

Fig. 9. Toroidal stress along the contour of the cavity, curvature of the contour and its change due to the deformation.

εn = 0.194. A careful observation shows that both the strain and the rigid-body rotation are proportional to the remote loading (the latter is a little different if the rigid-body rotation is treated as being finite, as shown by Eq. (55)). Therefore, using the assumption of small strains under the finite rotation of the rigid inclusion is unsuitable. However, following the work of Yang et al. (2017), we have still investigated the effect of finite rigid-body rotation in Section 3.1, regardless of the synchronous change of strain due to the stress concentration around the tip. According to Fig. 6, the rigid-body rotation 10 rises to α = 4.58◦ = 0.07985rad for the loading angle ϕ = 36 π and the effect of finite rigid-body rotation comes true when the remote loading is relatively large, say σ = 2.5 MPa in comparison with E = 22 MPa. 4.3. A poromechanics example The elastic isotropic material forming the solid phase of a porous medium has the shear modulus equal to 10 GPa and the Poisson ratio equal to 0.15, typically for concrete, and the plane strain assumption (κ = 3 − 4ν ) is used. The porous medium under consideration is saturated of a fluid with internal pressure P = 2 MPa (cf. Phan, 2008). The pore is taken to

W.-N. Zou, Q.-C. He / International Journal of Engineering Science 126 (2018) 68–96

85

Fig. 10. Hydrostatic pressure and maximal shear stress (MPa) due to a pentagram cavity with inner pressure are shown. The fields are enlarged  gradually to show the detailed structures, with the hydrostatic pressure σ0 =

1 2

(σ11 + σ22 )in the cases (a1)–(a3) and the maximal shear stress τmax =

 σ11 −σ22 2 2

2 + σ12

in the situations (b1)–(b3).

be a regular pentagram which is centered at the origin and whose circumcircle has the radius R = 5 mm. This pentagram is approximated by a polynomial with 10 coefficients, namely

4ι −14 214 −19 2908ι −34 3ι −4 13 −9 1231ι −24 20974 −29 ζ + ζ − ζ − ζ + ζ + ζ − ζ 10 225 125 11875 93750 2265625 390625 441199 −39 95763ι −44 6890609 − ζ + ζ + ζ −49 76171875 19531250 1708984375

φ (ζ ) = ζ +

(101)

and r = 0.7864R = 3.932 mm, as shown in Fig. 8. The maximal curvature reaches 287.17 mm−1 and the minimal radius of curvature is about 3.48 μm. Following our method, the potentials are determined as follows:

γ (w ) = rP (−0.2604ιw−4 − 3.515 × 10−2 w−9 + 2.018 × 10−2 ιw−14 + 1.044 × 10−2 w−19 − 6.436 × 10−3 ιw−24 − 5.980 × 10−3 w−29 + 2.516 × 10−3 ιw−34 + 4.449 × 10−3 w−39 − 7.027 × 10−4 ιw−44 − 4.032 × 10−3 w−49 ), (102)

86

W.-N. Zou, Q.-C. He / International Journal of Engineering Science 126 (2018) 68–96

ψ (w ) =

rP

φ  (w )

(−1.277w−1 + 0.2421ιw−6 + 0.2771w−11 − 0.2160ιw−16 − 0.1833w−21 + 0.2044ιw−26 + 0.1161w−31

− 0.2116ιw−36 − 6.005 × 10−2 w−41 + 0.2363ιw−46 ).

(103)

In particular, note that these potentials allow satisfying (14) due to the boundary condition of the cavity

σn = −P, τn = 0.

(104)

The stress state along the contour of the cavity is characterized by the toroidal normal stress σ θ . From σn + σθ = σ11 +

σ22 and (3)1 , we have (Savin, 1961)

!  N 1  γ  (η ) −k−1 σθ = 4Re  − σn = P − 4 Re kα η , z (η ) φ  (η ) k=1 k 

(105)

as shown in Fig. 9, together with the curvature of the boundary and its change. Their calculation formulae are derived in Appendix C. Apart from the local compression around the corner, most of the contour is in tension. The tension stress is about 6 MPa, while the compression stress is highly concentrated, about 80 MPa, 40 times the internal pressure, near the places with maximal curvature (Fig. 9(a)). It is interesting to observe from Fig. 9(b) that: (i) the top of the corner is quite flat; (ii) the change of curvature has a dramatic jump around the corner, namely the curvature increases for the contour outside the corner but decreases in the flat top, and even changes from positive curvature to negative curvature. This, together with the double peak characteristic of stress concentration, implies the possibility of bifurcation extension of the crack if the material strength is broken through. In Fig. 10, the fields are enlarged gradually to show the detailed structures, with the hydrostatic pressure σ 0 in the cases (a1)–(a3) and the maximal shear stress τ max in the situations (b1)–(b3). Besides the flat tops and the associated double peak characteristic of stress concentration, the illustrations exhibit a remarkable features of the stress fields due to the pore pressure of the cavity with multiple cusps, namely the concentration of hydrostatic pressure is separated while the maximal shear stress intends to circumferentially penetrate into each other. 5. Concluding remarks In this work, the fundamental problem of a 2D infinite elastic isotropic medium containing a rigid inclusion or a cavity and undergoing uniform remote loading has been studied by revisiting the powerful Kolosov–Muskhelishvili (K–M) potential theory. The global displacement of the rigid inclusion or the cavity relative to the matrix plays an important role in determining the K–M potentials. When the shape of subdomain is characterized by a Laurent polynomial with a finite number of terms, the K–M potentials are proved to possess also a finite polynomial expression, and a general procedure is elaborated to work out the solutions. All of them are explicit apart from one iterative step. In particular, in the situation where the shape of the rigid inclusion or the cavity is described by a Laurent polynomial with no more than three terms, the problem admits an explicit analytical solution. The present solution includes all the relevant results reported in the literature as its special cases. In light of the present method, some mistakes in the literature have been identified and corrected. Finally, besides their own usefulness, the results obtained in the present work can be applied as the benchmarks of numerical methods. Via some micromechanics schemes, they can also be applied to estimating the effective elastic moduli of 2D composites with rigid inclusions or cavities. Acknowledgments This work is supported by the National Science Foundation of China (NSFC), grant no. 11372124. The main work of this paper was done during the research stay of WNZ at the Université Paris-Est, Laboratoire de Modélisation et Simulation Multi Echelle, as a visiting scholar through QCH and under the financial supporting of China Scholarship Council (CSC), file no. 201506825048. Appendix A. Detailed derivations of explicit solutions A1. General derivation From (17), when the matrix extracting a subdomain ω is indicated by (18), the perturbation potentials, as the analytic functions in the domain including infinity, must take the forms

γ p (w ) = r

N 

αk w−k ,

k=1

ψ p (w ) =

2 N+1 r βk w−k , |w| ≥ 1, φ  (w ) k=1

(106)

W.-N. Zou, Q.-C. He / International Journal of Engineering Science 126 (2018) 68–96

87

and will be explicitly determined. To the best of our knowledge, such an algebraic structure of potentials has never been clarified before. The basic idea to solve the perturbation potentials is to expand all terms of (17) into the series of η, which can be regarded as the Fourier expansions with respect to the angular variable θ by introducing the expression η = eιθ , and the require every sum of the coefficients of the same power of η to be null. In order to make the derivation smooth, we introduce the equivalence relation ’ ∼ ’ between two functions F(η) and G(η) in the following. Precisely, F(η) ∼ G(η) means that the two functions have the same non-positive power part in their expansions with respect to η. Obviously,

F (η ) = G(η ) ⇒ F (η ) ∼ G(η ).

(107)

We proceed first to figure out γ p (η). Accounting for (107), it follows from (12) that

  z (η ) κγ p (η ) −   −1  γ¯ d η−1 ∼ C¯σ ,2 rη−1 + Cσ∗ ,1 r[φ (η ) − η] + 2μU0 . z¯ η −1

In the case where [z (η )]

(108)

can be expanded into a non-positive polynomial of η, namely any one of the roots of the −1

wN+1 z

algebraic equation (w ) = 0 has modulus less than the unit, [z¯ (η−1 )] powers of η. Thus, the non-positive polynomial

G

can be expanded in terms of non-negative



 z (η ) h φ (η ) η−1 ∼   ∼  −  −1  −1 ¯ ¯ r φ η rφ η

(109)

has the same maximal negative order N as that of φ (η). For |w| ≥ 1, we have dz = z (w )dw, meaning that dw = 0 ↔ dz = 0, and so that no root of wN+1 z (w ) = 0 has module equal to or greater than 1. Therefore, the above precondition certainly holds for all exterior-to-exterior conformal mappings. It is easy to see

G



η

−1



=

N 

 gk η

−k

∼

1−

k=0

N−1 

−1 l b¯ l η

l+1

N 

bk η−k ,

(110)

k=1

l=1

resulting in N 

gk η

k=0

−k

−

N 

gk η

k=0

−k

N−1 

l b¯ l ηl+1 ∼

l=1

N 

gk η

−k

−

k=0

N 



k=0

N−k −1

 l gk+l+1 b¯ l

η−k =

N 

bk η−k .

(111)

k=1

l=1

The balance of (111) gives rise to the explicit recursive formulae

gN−k = bN−k +

k−1 

l b¯ l gN−k+l+1 , k = 0, 1, . . . , N − 1

(112)

l=1

besides g0 =

N−1 ¯ l bl gl+1 . For instance, from gN = bN , gN−1 = bN−1 , we have l=1

gN−2 = bN−2 + b¯ 1 gN , gN−3 = bN−3 + b¯ 1 gN−1 + 2b¯ 2 gN ; . . . .

(113)

Substituting (9)1 and (109) into (108), and making use of the definition N 

ck η−k = C¯σ ,2 η−1 + Cσ∗ ,1 [φ (η ) − η] +

k=0

2μU0 , r

(114)

namely

c0 =

2μU0 , c1 = C¯σ ,2 + Cσ∗ ,1 b1 , ck = Cσ∗ ,1 bk , k = 2, . . . , N, r

we have ∞  n=1

καn η−n +

N  k=2

gk η−k

∞  l=1

l α¯ l ηl+1 ∼

∞  n=1

καn η−n +

N−2  n=0



N−n −1 l=1

(115)

 l α¯ l gn+l+1

η−n =

N  k=0

ck η−k .

(116)

88

W.-N. Zou, Q.-C. He / International Journal of Engineering Science 126 (2018) 68–96

The balance of both sides of (116) yields αn = 0 for all n > N, that is to say, the potential γ p (η) , when expanded as a polynomial of η, has the same maximal negative order as ω(η) (Chang & Conway, 1968). In detail, the relations derived from (116) consist of a system of linear equations N−k −1

καk = ck −

pgk+ p+1 α¯ p , k = 1, . . . , N

(117)

p=1

for {α1 , . . . , αN }, besides

2μU0 = rc0 = r

N−1 

kα¯ k gk+1 .

(118)

k=1

The linear equations (117) can be solved by an iterative process (n = 0, 1, 2, . . .)



(n+1 )

αk

=κ

−1

ck − κ

−1

N−k −1

"

(n )

l gk+l+1 α¯ l , k = 1, . . . , N ,

(119)

l=1

starting from

αk(0) = κ −1 ck , k = 1, . . . , N − 2.

(120)

Up to now, the only undetermined parameter is the rigid-body rotation α in σ ,1 . For the potential ψ p (η), multiplying the conjugate of (17) by φ  (η), we have C∗

      κφ  (η )γ¯ d η−1 − r−1 z¯ η−1 γd  (η ) − φ  (η )ψ p (η ) ∼ Cσ ,2 rφ  (η )φ (η ) + C¯σ∗ ,1 rφ  (η )φ¯ η−1 + 2μU¯ 0 φ  (η ).

Accounting for (18), (106)1 and (9)2 , and using the derivation like

φ

(η )

N 



ak η = k

1−

N 

k=1



kbk η

−k−1

k=1

N 

al η ∼ − l

N−1 

kbk ak+1 −

k=1

l=1



N N   p=1

(121)

 kbk ak−p+1

η−p ,

(122)

k= p

we obtain

φ  (η )ψ p (η ) ∼ φ  (η )ψ p (η ), φ  (η )φ (η ) ∼ φ  (η )φ (η ) − η,   N−1 N N      φ  (η )γ¯ d η−1 ∼ −r kbk α¯ k+1 − r kbk α¯ k−p+1 η−p , k=1

φ z¯



(η )φ¯ η



η

−1



r2

−1



∼η

−1

φ

p=1

(η ) −

N−1 

k= p

kbk b¯ k+1 −

k=1



N N   p=1

 kbk b¯ k−p+1

η−p ,

k= p

    N−1 N N   h¯ γ p  (η )   −1 ¯ ¯ γ p (η ) ∼ η + − kαk bk+1 − kαk bk−p+1 η−p . r

r

k=1

p=0

k= p

Introducing the notation

λk =

N 



p

 α p b¯ p−k+1 − κ b p α¯ p−k+1 + C¯σ∗ ,1 b p b¯ p−k+1 , k = 1, . . . , N,

(123)

p=k

the balance of all non-positive-power terms in (121) and the use of (118) yield

C=

N−1 

k



 αk b¯ k+1 − κ bk α¯ k+1 + C¯σ∗ ,1 bk b¯ k+1 − c¯0

k=1

=

N−1 

k



 αk b¯ k+1 − κ bk α¯ k+1 + bk c¯k+1 − αk g¯ k+1 = 0

(124)

k=1

for the constant term, and

φ  (η )ψ p (η ) = r

    h¯ λk η−k − η−1 + γ p  (η ) − C¯σ∗ ,1 rη−1 φ  (η ) + c¯0 r 1 − φ  (η ) + Cσ , 2 r η − φ  (η )φ (η )

N 

r

k=1

(125)

for the negative-power terms, where the highest order of the negative powers is 2N + 1 from the last term. Making use of the relations

gk+1 = bk+1 +

N−k −2 l=1

l b¯ l gk+l+2 ,

καk+1 = ck+1 −

N−k −2 p=1

pgk+ p+2 α¯ p ,

W.-N. Zou, Q.-C. He / International Journal of Engineering Science 126 (2018) 68–96

89

and noting that

καN = cN , καN−1 = cN−1 ,

gN = bN , gN−1 = bN−1 ;

(126)

we can derive N−1 

C=

k



αk b¯ k+1 − αk g¯ k+1 − κ bk α¯ k+1 + bk c¯k+1



k=1 N−3 

=

kbk

N−k −2

k=1 N−3 

=

pα p g¯ p+k+2 −

N−k −2

kbk g¯ p+k+2 −

k=1

p=1

kαk

N−k −2

k=1

p=1

pα p

N−3 

pb p g¯ p+k+2

p=1

N−3 

kαk

N−k −2

k=1

pb p g¯ p+k+2 ≡ 0,

(127)

p=1

meaning that the relation (124) is an identity, and the remaining relations (125) result in the explicit expression of ψ p (w ):

φ  (w ) r

ψ p (w ) =

N 

λk w−k +

k=1

N 



kαk w−k−2 +

k=1

+ Cσ ,2

N 2 N+1  kbk w−k + k=1

k=1



k=3



N  h¯ k α + c¯0 bk w−k−1 r k



min (N,k−2 )



pb p bk−p−1 w

p=max (1,k−N−1 )

! −k

 ¯∗

− Cσ ,1 w

−1

N  − kbk w−k−2

 (128)

k=1

with the highest negative exponent 2N + 1. Combining with the basic potentials (8), we finally obtain the total potentials for the rigid inclusion problem with remote uniform stresses as:

1 2

γ (w ) = Cσ , 1 rφ (w ) + r

N 

αk w−k ,

(129)

k=1

1 2

ψ (w ) = Cσ , 2 rφ (w ) − Cσ , 1 h¯ +

2 N+1 r βk w−k . φ  (w )

(130)

k=1

It is interesting to remark that the expansion

1 2

φ  (w )ψ (w ) = Cσ , 2 rw − C¯σ∗ ,1 rw−1 φ  (w ) − Cσ , 1 h¯ φ  (w ) + r 

h¯ + r w−1 + r



N 

N 

λk w−k + 2μU0

k=1

N 

kbk w−k−1

k=1

kαk w−k−1

(131)

k=1

associated with the total potential ψ is more compact than φ  (w )ψ p (w ), having the highest negative exponent N + 2 of w. This feature may induce the misunderstanding to some extent concerning the structure of ψ p (w ) (Misseroni et al., 2014). For simplicity, it is preferable to write the potential ψ (w ) in the form



 r ψ (w ) =  C w+ βk w−k φ (w ) σ , 2 N+2



−

k=1

with N+2 

1 Cσ , 1 h¯ , 2

(132)

  N N   h¯  βk w−k = c¯0 1 − φ  (w ) − C¯σ∗ ,1 w−1 φ  (w ) + w−1 + kαk w−k−1 + λk w−k . r

k=1

k=1

(133)

k=1

It remains to determine the rigid-body rotation α . To this end, we calculate the resultant moment applied on the interface of ω by





Mo = Re

∂ω

d

 ψ (z )dz − zψ (z ) − zz¯γ  (z ) .

(134)

According to the aforementioned derivations, the potentials γ and ψ are single-valued, requiring the resultant moment to be null yields (Muskhelishvili, 1953, Page 364)



0 = Mo = Re = −2π r 2 Im

∂ω

ψ (z )dz = r2 Re

N  k=1

k







|η|=1

 λ1 − C¯σ∗ ,1 η−1 dη !  ∗ ∗

αk b¯ k − κ bk α¯ k + C¯σ ,1 bk b¯ k − C¯σ ,1 .

(135)

90

W.-N. Zou, Q.-C. He / International Journal of Engineering Science 126 (2018) 68–96

Reminding (11), we obtain



2μα =

1−

N 

−1

kbk b¯ k

Im

N 

k=1

k



κ bk α¯ k − αk b¯ k



k=1

= π r 2 |ω|−1 Im

N 

k



 κ bk α¯ k − αk b¯ k ,

(136)

k=1

   ¯ where use is made of the area formula |ω| = π r 2 1 − N k=1 kbk bk of the rigid inclusion ω . However, the above formula is an implicit relation for the relative rigid-body rotation since the iterative process (119) used to solve {αk , k = 1, . . . , N } actually involves the parameter α . Therefore, the relation (136) has to be combined with (119) so as to constitute an extended iterative system, starting from α (0 ) = 0. In our practice, this iterative procedure turns to be very efficient. p In the case of the cavity, due to κ = −1 and Cσ∗ ,1 = Cσ ,1 being real, the vanishing of the resultant moment holds identically, and so the solution of iterative process without α becomes simpler. In addition, such a property has the surprising p consequence for N > 1: the highest negative exponent reduces to N + 1 since βN+2 = αN + Cσ ,1 bN ≡ 0 due to (115) and (117). A careful investigation reveals that, when the mapping function (18) has few terms, say two or three, the algebraic equations (117) for {α1 , . . . , αN }, plus (136) in the case of rigid inclusion, become easy to be solved and so the completely explicit solutions of potentials are available. In summary, for the elasticity problem of a remotely loaded infinite matrix extracting a subdomain specified by (18), as a rigid inclusion or a cavity, we list the solution steps as follows: calculate the coefficients gk , solve the parameters {α1 , . . . , αN } and α from (117) and (136) with the iterative procedure, get U0 from (118), and obtain the potentials γ and ψ from (129) and (132). Finally, the displacement, stress and strain fields are calculated through (2)–(4). A2. Consideration of finite rigid-body rotation Based on the descriptions of (22) and (23), where w is the same complex variable of the image plane for whatever the matrix is deformed or not, the basic balance equation (12) on the boundary of the rigid inclusion for the perturbation potentials is changed to be

      z (η ) κγ p (η ) −   −1  γ¯ d η−1 − ψ¯ d η−1 = Cσα,1 reια φ (η ) + C¯σ ,2 re−ια φ¯ η−1 + 2μU0 z¯ η

(137)

where the basic potentials are assumed the same as (8) but using H instead of h and

Cσα,1 = −

κ −1 2





Cσ , 1 + 2μ 1 − e−ια .

(138)

The parameters in the potentials (25) are then determined simply by some replacements as follows. First, {gk } and {λk } share the same definitions as (110) and (123), respectively. The coefficients {α k , α } are iteratively solved from (119) and

2μ sin α = π r 2 |ω|−1 Im

N 

k



 κ bk α¯ k − αk b¯ k ,

(139)

k=1

where {ck } are modified to be

c1 = e−2ια C¯σ ,2 + Cσα,1 b1 , ck = Cσα,1 bk , k = 2, . . . , N

(140)

with Cσα,1 defined by (138). U0 is calculated by

2μU0 = rc0 eια = reια

N−1 

kα¯ k gk+1 ,

(141)

k=1

whilst {β k } are determined by N+2 

βk w

−k







= c¯0 1 − φ  (w ) − C¯σα,1 w−1 φ  (w ) +

k=1

w

−1

H¯ + eια r



N  k=1

kαk w−k−1 +

N 

λk w−k .

(142)

k=1

Appendix B. Laurent polynomials approximating rectangle, rhombus and pentagram The conformal map from the exterior of the unit disk to the exterior of a polygon having vertices z1 , z2 , . . . , zL and exterior angles α1 π , α2 π , . . . , αL π in counterclockwise order has the Schwarz–Christoffel formula (Savin, 1961)

z (w ) = h + r



w 1

L  

k=1

1−

 ak αk −1 ζ

dζ , |ak | = 1, k = 1, . . . , L

(143)

W.-N. Zou, Q.-C. He / International Journal of Engineering Science 126 (2018) 68–96

91

satisfying L 

L 

αk = L + 2,

k=1

(αk − 1 )ak = 0, zk = z(wk ) for k = 1, 2, . . . , L.

(144)

k=1

For the polygon being a rectangle with vertices

z1 =

1 1 1 1 1 1 1 1 l x − ιl y , z 2 = l x + ιl y , z 3 = − l x + ιl y , z 4 = − l x − ιl y , 2 2 2 2 2 2 2 2

(145)

we have

L = 4,

3 2

α1 = α2 = α3 = α4 = ,

(146)

and it is reasonable by assuming

a1 = e−ιπη , a2 = eιπ η , a3 = eιπ (1−η ) , a4 = eιπ (1+η ) ,



0<η<

1 2

 (147)

to maintain the symmetry. Using the above parameters, we can calculate the series expansion of the conformal map within a constant. From the relation

φ (w ) =



w



1−

2 cos 2ηπ 1 + 4 t2 t

 12

 dt =

∞ w

ck t −2k dt =

k=0

∞  k=0

ck t 1−2k , 1 − 2k

(148)

we obtain the recursive formulae of the coefficients

c0 = 1, c1 = − cos 2ηπ , ncn = (2n − 3 )cn−1 cos 2ηπ − (n − 3 )cn−2 .

(149)

For the special case ly /lx = 0.25, the truncation with 8 terms is

φ (w ) = w + 0.56325w−1 − 1.1379 × 10−1 w−3 − 3.8456 × 10−2 w−5 − 7.1478 × 10−3 w−7 + 4.1618 × 10−3 w−9 + 5.1512 × 10−3 w−11 + 2.2116 × 10−3 w−13 − 6.0667 × 10−4 w−15 . If the shape becomes square by letting lx = ly , we must have η =

r/lx =

 





1 1 1 1 RF , 0, 1 − Rd , 0, 1 2 2 3 2

−1

1 4

(150)

due to the symmetry, and so

 0.59017

(151)

and

φ (w ) = w +

∞  (− )k (2k − 3 )!! 1−4k w . (4k − 1 )(2k )!!

(152)

k=1

In order to obtain η and r in general, we can calculate the side lengths from as

lx = r

 2π −πη  π +π η

= 4r cos π η

1 − e2ι(πη−θ )

 

 2π −πη  π +πη

1 − e−2ι(πη+θ )

1/2

ιeιθ dθ

1/2  π /2−πη   π 2 sin θ 1 1− d θ = 4 r cos π η E − π η , ( ) 2 cos π η cos2 π η 0

= 4r cos2 π η RF sin

ly = r

1/2 

2

 1   π η, 0, 1 − Rd sin2 π η, 0, 1 , 3

1 − e2ι(πη−θ )

1/2 

1 − e−2ι(πη+θ )

1/2

eιθ dθ

1/2   θ 1 π η = 4r sin π η 1− d θ = 4 r sin π η E , ( ) 2 sin π η 0 sin π η       1 = 4r sin 2 π η RF cos2 π η, 0, 1 − Rd cos2 π η, 0, 1  πη 

2

sin

3

where E(φ , k) is the Legendre elliptic integral of the second kind in Legendre’s notation, while RF /RD are the Legendre elliptic integrals of the first/second kinds in Carlson’s notation (Press, Teukolsky, Vetterling, & Flannery, 1992). So the scalar η can be solved from



tan

2

πη





RF cos2 π η, 0, 1 − 13 Rd cos2 π η, 0, 1



2

RF sin



π η , 0, 1 −

1 R 3 d



2

sin

π η , 0, 1



 =

ly lx

(153)

92

W.-N. Zou, Q.-C. He / International Journal of Engineering Science 126 (2018) 68–96 Table B.1 Parameters η and r for different aspect ratios ly /lx of rectangles. ly /lx

1

1/2

1/4

η

0.25 0.59017032

0.20033011 0.4373786

0.1547744572 0.3538574

r/lx

Table B.2 Parameter r and β for different aspect ratios ly /lx of rhombuses. ly /lx

9/15

4/15

2/15

r/lx

0.33882660 0.65595826

0.28410851 0.83409536

0.26588766 0.91561507

β

and r is determined by

r/lx =





 2  1  2  1 cos−2 π η/ RF sin π η, 0, 1 − Rd sin π η, 0, 1 . 4 3

(154)

Table B.1 lists some numerical results for rectangles. The second important shape is rhombus. Letting the vertices located on the axes, namely

z1 =

1 1 1 1 l x , z 2 = ιl y , z 3 = − l x , z 4 = − ιl y , 2 2 2 2

(155)

we have the completely explicit parameters for L = 4

α1 = α3 = 2 −

2

π

arctan

ly ≡ 1 + β, lx

α2 = α4 = 2 − α1 = 1 − β

(156)

with β ∈ (−1, 1 ), and the only choice

a1 = 1, a2 = ι, a3 = −1, a4 = −ι

(157)

due to the symmetry. Thus, from the relation

    1 β ι 1−β  1 β ι 1−β 1− 1+ 1− dt t t t t  w ∞ ∞  ck = ck t −2k dt = w1−2k 1 − 2k

φ (w ) =



w



1−

k=0

(158)

k=0

we can derive the recursive formulae

c0 = 1, c1 = 1 − 2β , ncn = (n − 3 )cn−2 + (1 − 2β )cn−1 ,

(159)

and the truncation with 8 terms for ly /lx = 2/15 is

φ (w ) = w + 0.83123w−1 + 5.1509 × 10−2 w−3 − 8.5632 × 10−3 w−5 + 6.7899 × 10−3 w−7 − 2.7809 × 10−3 w−9 + 2.4756 × 10−3 w−11 − 1.3489 × 10−3 w−13 + 1.2561 × 10−3 w−15 .

(160)

The size parameter r can be determined from





1 1 1 π − lx + ιly = − lx 1 − ι tan (1 − β ) 2 2 2 2  π /2  β  1−β  β  1−β ιθ =r 1 − e−ιθ 1 − ιe−ιθ 1 + e−ιθ 1 + ιe−ιθ ιe d θ 0

= 2ιr eιβπ /2 to be



r/lx =



1 0



−β /2

xβ 1 − x2

1 −1 βπ 1+β β sin /B ,1 − 2 2 2 2



1 1+β β dx = 2ιr eιβπ /2 B ,1 − 2 2 2



 ,

(161)

where is the Beta function. Some numerical results for rhombuses are listed in Table B.2. k−1 For the L-side regular polygons with vertices zk = Reι L 2π , k = 1, 2, . . . , L, we have the geometric parameters as

αk =

L+2 k−1 π , ak = eι L 2π , k = 1, 2, . . . , L, L

(162)

W.-N. Zou, Q.-C. He / International Journal of Engineering Science 126 (2018) 68–96

93

and so from

φ (w ) =



w

L  

1 − t −1 eι

2kπ L

 2L

 dt =

w



1 − t −L

 2L

dt

k=0

=

∞  k=0

ck w1−Lk , 1 − Lk

(163)

one can derive the recursive relations between the coefficients (Jasiuk, 1995)

# 2 kp=1 ( pL − 2 ) L Lk − 2 c0 = 0, c1 = − and ck+1 = c = ··· = − , k > 0. 2 (k + 1 )L k (k + 1 )!Lk+1 From the relation

z2 − z1 π L+2 = Reι 2L π sin =r 2 L



π L



0

1−e

−ιLθ

 2L

2 L+2 ιeιθ dθ = 2 L reι 2L π



π L



0

Lθ sin 2

(164)

 2L

dθ

2  − 12   r ι L+2 π 1 Lθ L Lθ Lθ 2 =2 e 2L sin 1 − sin d sin L 2 2 2 0  1   1 L+2 1 2 r L+2 2−L 2 r L+2 = 2 L eι 2L π x 2L (1 − x )− 2 dx = 2 L eι 2L π B , L L 2L 2 0 L+2 L

we derive the formula of the size parameter

r/R = 2−2/L L sin

π L

/B

L + 2 1 2L

,

.

2

(165)

For example, we have r/R = 0.7305 for the equilateral triangle and r/R = 0.92037 for the regular hexagon. The final shape concerned in this paper is pentagram. We define its vertices as

z2k+1 = Rιeι2kπ /5 , z2k+2 = −Rι

cos 2π /5 ι2(k+3 )π /5 e , k = 0, 1, 2, 3, 4, sin 3π /10

(166)

and consequently have the following parameters for L = 10

9 5

3 5

α1 = α3 = α5 = α7 = α9 = π , α2 = α4 = α6 = α8 = α10 = π .

(167)

According to the symmetry, the imaging points of the vertices on the unit circle can be chosen to be

a2k+1 = ιeι2kπ /5 , a2k+2 = −ιeι2(k+3 )π /5 , k = 0, 1, 2, 3, 4, so from

φ (w ) =

 

=

4 w



k=0 w



1−

1−

ιe ι 2 k π / 5

4/5

t

ι 4/5  t5

4 

 1+

ιe ι 2 k π / 5

l=0

1+

ι −2/5 t5

−2/5 dt

t

dt = w +

(168)

∞  k=1

ck w1−5k , 1 − 5k

(169)

one can derive the recursive relations between the coefficients

c0 = 1 , c1 = −





6 6 12 ι and ncn = − ιcn−1 − n − cn−2 , n > 1. 5 5 5

(170)

The truncation with 10 terms is

4ιw−14 214w−19 2908ιw−34 441199w−39 3ιw−4 13w−9 1231ιw−24 20974w−29 + − − + + − − 10 225 125 11875 93750 2265625 390625 76171875 95763ιw−44 6890609w−49 + + . (171) 19531250 1708984375

φ (w ) = w +

The size parameter r can be evaluated from

z3 − z2 = Rιeι2π /5 + Rι =r

cos 2π /5 ι6π /5 π e = −R tan sin 3π /10 5

 π /2+2π /5  π /2+π /5

1 − ιe−5ιθ

4/5 

1 + ιe−5ιθ

−2/5

ιeιθ dθ

94

W.-N. Zou, Q.-C. He / International Journal of Engineering Science 126 (2018) 68–96

= −r

 2π /5  π /5

1 − e−5ιθ

4/5 

1 + e−5ιθ

−2/5

eιθ dθ

4/5  −2/5  2π /5  5θ 5θ sin cos dθ 2 2 π /5  π /2 7/5 2 =− r (cos θ )4/5 (sin θ )−2/5 dθ 5 0  1 3 9 2 2/5 2 2/5 =− r x−7/10 (1 − x )−1/10 dx = − rB , 5 5 10 10 0 = −22/5 r

to be



r/R =

3 9 22/5 B 10 , 10 π 5 tan 5

 = 0.790707.

(172)

Appendix C. Curvature of the boundary and its change For a plane curve given by the parametric representation

x1 = x1 (t ), x2 = x2 (t ),

(173)

its curvature has expression

K0 =

x˙ 1 x¨2 − x˙ 2 x¨1



x˙ 21 + x˙ 22

3/2 .

(174)

When the curve is materialized, the new curvature becomes

Ku =

(x˙ 1 + u˙ 1 )(x¨2 + u¨ 2 ) − (x˙ 2 + u˙ 2 )(x¨1 + u¨ 1 ) ,  3/2 (x˙ 1 + u˙ 1 )2 + (x˙ 2 + u˙ 2 )2

(175)

where ui (i = 1, 2 ) are the Cartesian components of the displacement. Dundurs (1989) derived the first-order approximation of the curvature change for small strain as

K = Ku − K0      1 3 ∂ ε12 ∂ ε11 ≈ K0 ε11 + ε22 − 2 x˙ 21 ε11 + 2x˙ 1 x˙ 2 ε12 + x˙ 22 ε22 + 3 x˙ 31 2 − ∂ x1 ∂ x2 R R       ∂ ε22 ∂ ε11 ∂ ε11 ∂ ε22 ∂ ε12 ∂ ε22 + x˙ 21 x˙ 2 2 − − x˙ 1 x˙ 22 2 − − x˙ 32 2 − , ∂ x1 ∂ x1 ∂ x2 ∂ x2 ∂ x2 ∂ x1

(176)

where R2 = x˙ 21 + x˙ 22 . In this paper, using the complex variable z = x1 + ιx2 , and denoting by

τ=

dz z˙ = ; ds |z˙ |

1 2

1 2

ε0 = (ε11 + ε22 ), ε2 = (ε22 + ε11 ) + ιε12 ,

(177)

we can rewrite (174) and (176) as

K0 = Im(τ¯ z¨ )/|z˙ |

2

and

(178)

    ∂ε ∂ε ∂ε K ≈ K0 Re 3τ 2 ε2 − ε0 + Im 2τ 0 + τ 3 2 + 3τ 2 , ∂z ∂z ∂ z¯

(179)

respectively. In the context of linear elasticity, by virtue of (4), the strain in the above formula can be calculated from the potentials as















2μK ≈ K0 Re 3τ 2 z¯γ  (z ) + ψ  (z ) − 2(κ − 1 )γ  (z ) + Im (2κ + 1 )τ γ  (z ) + τ 3 z¯γ  (z ) + ψ  (z )



For the material curve to be the boundary of the cavity, the mapping z = rφ (η ) = rφ e (173) can be replaced by θ , and so

τ=

ιηφ  (η ) . |φ  (η )|

 ιθ

.

(180)

means that the parameter t in

(181)

W.-N. Zou, Q.-C. He / International Journal of Engineering Science 126 (2018) 68–96

95

The curvature and its change have formulae

r K0 =

1



Re 1 +

|φ  (η )|

 ηφ  (η ) , φ  (η )

(182)

  3 ηφ  (η )  ηφ (η )  γ z − r φ (η )γ  (z ) + ψ  (z ) ( )   φ η φ η | ( )| | ( )| " η 2 φ  (η )  − 3K0 r φ (η )γ  (z ) + ψ  (z ) − 2(κ − 1 )K0 γ  (z ) φ  (η )

2μK ≈ Re (2κ + 1 )

(183)

with

r γ  (z ) = r 3 γ  (z ) =

γ  (η ) 2  γ  (η ) φ  (η )γ  (η ) , r ψ (z ) = − ,  2 3 φ (η ) [φ  (η )] [φ  (η )]   2  φ (η ) γ (η ) γ  (η ) φ  (η )γ  (η ) [φ  (η )]

3

−3

[φ  (η )]

4

+3

[φ  (η )]

5

−

φ  (η )γ  (η ) 4 [φ  (η )]

(184)

and similar formulae of rψ  (z) and r2 ψ   (z). References Chang, C. S., & Conway, H. D. (1968). A parametric study of the complex variable method for analyzing the stresses in an infinite plate containing a rigid rectangular inclusion. The International Journal of Solids and Structures, 4(11), 1057–1066. Duan, H. L., Wang, J., Karihaloo, B. L., et al. (2006). Nanoporous materials can be made stiffer than non-porous counterparts by surface modification. Acta Materialia, 54, 2983–2990. Dundurs, J. (1989). Cavities vis-a-vis rigid inclusions and some related general results in plane elasticity. Journal of Applied Mechanics, 56, 786790. Ekneligoda, T. C., & Zimmerman, R. W. (2008). Shear compliance of two-dimensional pores possessing n-fold axis of rotational symmetry. Proceedings of the Royal Society A, 464, 759–775. England, A. H. (1971). Complex variable methods in elasticity. London; New York: Wiley-Interscience. Eshelby, J. D. (1957). The determination of the elastic field of an ellipsoidal inclusion and related problems. Proc. Roy. Soc. London A, 241, 376–396. Eshelby, J. D. (1959). The elastic field outside an ellipsoidal inclusion. Proc. Roy. Soc. London A, 252, 561–569. Fan, T., & Yang, L. (2017). Effective youngs modulus of nanoporous materials with cuboid unit cells. Acta Mechanica, 228, 21–29. doi: 10.1007/ s00707- 016- 1682- 6. Haller, X., Monerie, Y., & Pagano, S. (2016). Elastic behavior of porous media with spherical nanovoids. The International Journal of Solids and Structures, 84, 99–109. Jasiuk, I. (1995). Cavities vis-a-vis rigid inclusions: Elastic moduli of materials with polygonal inclusions. The International Journal of Solids and Structures, 32(3/4), 407–422. Kachanov, M., Shafiro, B., & Tsukrov, I. (2003). Handbook of elasticity solutions. Springer-Science + Business Media, B.V. Kachanov, M., Tsukrov, I., & Shafiro, B. (1994). Effective modulus of solids with cavities of various shapes. Applied Mechanics Reviews, 47, 151–174. Kang, H., & Milton, G. W. (2008). Solutions to the PólyaSzegö conjecture and the weak Eshelby conjecture. Archive for Rational Mechanics and Analysis, 188(1), 93–116. Liu, L. P. (2008). Solutions to the Eshelby conjectures. Proc. Roy. Soc. A, 464, 573–594. Lu, J. K. (1995). Complex variable method in plane elasticity. Singapore: World Scientific. Markenscoff, X., & Paukshto, M. (1995). The correspondence between cavities and rigid inclusions in three-dimensional elasticity and the cosserat spectrum. The International Journal of Solids and Structures, 32, 431–438. Michler, G. H., Adhikari, R., & Henning, S. (2004). Micromechanical properties in lamellar heterophase polymer systems. The Journal of Materials Science, 39, 3281–3292. Mindeguia, J. C., Pimienta, P., Noumowe, A., & Kanema, M. (2010). Temperature, pore pressure and mass variation of concrete subjected to high temperature experimental and numerical discussion on spalling risk. Cement and Concrete Research, 40, 477–487. Misseroni, D., Dal Corso, D., Shahzad, S., & Bigoni, D. (2014). Stress concentration near stiff inclusions: Validation of rigid inclusion model and boundary layers by means of photoelasticity. Engineering Fracture Mechanics, 121-122, 87–97. Muskhelishvili, N. I. (1953). Some basic problems of the mathematical theory of elasticity. Dordrecht: Springer. Naleway, S. E., Porter, M. M., McKittrick, J., & Meyers, M. A. (2015). Structural design elements in biological materials: application to bioinspiration. Advanced Materials, 27, 5455–5476. Öztürk, T., Mirmesdagh, J., & Ediz, T. (1994). Strain partitioning and plastic flow in some metal/metal laminates. Materials Science and Engineering A, 175, 125–129. Penzien, J. (20 0 0). Seismically induced raking of tunnel linings. Earthquake and Engineering Structural Dynamics, 29, 683–691. Phan, L. T. (2008). Pore pressure and explosive spalling in concrete. Materials and Structures, 41, 1623–1632. Press, W. H., Teukolsky, S. A., Vetterling, W. T., & Flannery, B. P. (1992). Numerical recipes in Fortran 77: The art of scientific computing. (2nd). Press Syndicate of the University of Cambridge. Rezaeepazhand, J., & Jafari, M. (2010). Stress concentration in metallic plates with special shaped cutout. The International Journal of Mechanical Sciences, 52, 96–102. Sadd, M. H. (2005). Elasticity: Theory, applications, and numerics. Elsevier Butterworth–Heinemann. Savin, G. N. (1961). Stress concentration around holes. Pergamon Press. Sharma, D. S. (2012). Stress distribution around circular/elliptical/triangular holes in infinite composite plate. Engineering Letters, 20, 1–9. Tang, Z., Kotov, N. A., Magonov, S., & Öztürk, B. (2003). Nanostructured artificial nacre. Nature Materials, 2, 413–418. Tsukrov, I., & Novak, J. (2002). Effective elastic properties of solids with defects of irregular shapes. The International Journal of Solids and Structures, 39, 1539–1555. Tsukrov, I., & Novak, J. (2004). Effective elastic properties of solids with two-dimensional inclusions of irregular shapes. The International Journal of Solids and Structures, 41, 6905–6924. Wang, B., Yang, W., McKittrick, J., & Meyers, M. A. (2016). Keratin: structure, mechanical properties, occurrence in biological organisms, and effort s at bioinspiration. Progress in Materials Science, 76, 229–318.

96

W.-N. Zou, Q.-C. He / International Journal of Engineering Science 126 (2018) 68–96

Wegst, U. G., Bai, H., Saiz, E., Tomsia, A. P., & Ritchie, R. O. (2015). Bioinspired structural materials. Nature Materials, 14, 23–36. Yang, Q. Z., Liu, Q. C., Yue, Z. F., Li, X. D., & Xu, B. X. (2017). Rotation of hard particles in a soft matrix. Journal of the Mechanics and Physics of Solids, 101, 285–310. Zappalorto, M., & Lazzarin, P. (2011). Stress fields due to inclined notches and shoulder fillets in shafts under torsion. The Journal of Strain Analysis for Engineering Design, 46, 187–199. Zhao, G. P., & Yang, S. L. (2015). Analytical solutions for rock stress arounds quare tunnels using complex variable theory. International Journal of Rock Mechanics and Mining Sciences, 80, 302–307. Zou, W. N., & Lee, Y. G. (2017). Completely explicit solutions of Eshelby’s problems of smooth inclusions embedded in a circular disk, full- and half-planes. Acta Mech. doi:10.10 07/s0 0707- 017- 2058- 2. Eshelby, J. D. (1961). Elastic inclusion and inhomogeneities. In I. N. Sneddon, & R. Hill (Eds.), Progress in Solid Mechanics, vol. 2 (pp. 89–140). Amsterdam: North-Holland. Zou, W. N., He, Q. C., Huang, M. J., & Zheng, Q. S. (2010). Eshelby’s problem of non-elliptical inclusions. Journal of the Mechanics and Physics of Solids, 58(3), 346–372.