On dependence of the stress intensity factor and T-stress from imposed boundary conditions in a rectangular cracked plate

Available online at www.sciencedirect.com

Computational Materials Science 42 (2008) 149–155 www.elsevier.com/locate/commatsci

On dependence of the stress intensity factor and T-stress from imposed boundary conditions in a rectangular cracked plate Y.Z. Chen *, X.Y. Lin Division of Engineering Mechanics, Jiangsu University, Zhenjiang, Jiangsu 212013, PR China Received 1 February 2007; received in revised form 18 June 2007; accepted 28 June 2007 Available online 23 August 2007

Abstract This paper investigates the dependence of the stress intensity factor (abbreviated as SIF) and T-stress from the imposed boundary conditions. A rectangular cracked plate is subject to a tension in the vertical direction. In the first boundary condition, the tension loading is applied at the tops of plate. In the second boundary condition, the relevant assumed displacement is applied on the edges of plate. In the study, the complex potentials are assumed in the series form with undetermined coefficients, and those coefficients can be evaluated by using the variational principle. The aim of study is to obtain the SIF and T-stress at the crack tips for two different boundary conditions. Finally, numerical results are presented. The obtained results may have their application in the fatigue test of a cracked component. Ó 2007 Elsevier B.V. All rights reserved. Keywords: Stress intensity factor; T-stress; Dependence of SIF and T-stress; Influence of boundary conditions

1. Introduction In the stress expansion form near a crack tip, the first term is singular which has a relation with the stress intensity at crack tip. The second term is finite and bounded. In the notation of Rice [1], the second term is denoted as the T-stress and can be regarded as the stress acting parallel to the crack flanks. The T-stress evaluation may have engineering application in the following fields. For example, The T-stress can considerably influence the plastic zone near crack tip in the case of small scale yielding [2]. In addition, the T-stress plays an important role in directional stability for the crack growth path [1,3]. Therefore, the T-stress evaluation got much attention by many researchers [2–5]. The traction boundary condition and the assumed displacement boundary condition are two kinds of condition, which are generally encountered in the analysis as well as in *

Corresponding author. E-mail address: [email protected] (Y.Z. Chen).

0927-0256/$ - see front matter Ó 2007 Elsevier B.V. All rights reserved. doi:10.1016/j.commatsci.2007.06.014

the test of a cracked component. The conditions must significantly influence behavior of solution for the stress intensity factor and the T-stress at the crack tip. The relevant problem is important in the fatigue test of a cracked component. In the test for evaluating the propagating rate of crack for a cracked plate with respect to the various boundary conditions, researcher needs to preserve the SIF value on a constant value. However, there are two cases encountered. In the first case, the tension po is applied on the upper and lower edges of the cracked plate (Fig. 1a). In the second case, the assumed displacements in Fig. 1b are uo = mpob/E and vo = poh/E. If the plate is perfect with no crack, the mentioned two states are actually in the same level of tension. Since the loading is a tension in the first case, the SIF value must become larger when the crack is growing. In this case, one must adjust the tension loading a little lower to preserve the SIF value on a constant (Fig. 1a). On contrary, the boundary displacements are assumed on the boundaries of the cracked plate in the second case (Fig. 1b). In this case, SIF value must become smaller when the crack is growing.

150

Y.Z. Chen, X.Y. Lin / Computational Materials Science 42 (2008) 149–155

a

y

u = −u o x / b

b

po

y

v = vo

u = uo v = yvo / h

o 2h

2a

o

x 2h

x u = −u o

2a

v = yvo / h 2b

2b u = −u o x / b v = −vo

In this paper, the dependence of the SIF and T-stress from the imposed boundary conditions is investigated. A rectangular cracked plate is subject to a tension in vertical direction. In the first boundary condition, the tension loading is applied at the tops of plate. In the second boundary condition, the relevant assumed displacement is applied on the edges of plate. In the study, the complex potentials are assumed in the series form with undetermined coefficients, and those coefficients can be evaluated by using the variational principle. The aim of study is to obtain the SIF and T-stress at the crack tips for two different boundary conditions. Finally, numerical results are presented. 2. Analysis

Fig. 1. (a) A rectangular cracked plate with tension po, (b) A rectangular cracked plate with an equivalent boundary displacements to case (a), uo = mpob/E and vo = poh/E.

In this case, one must adjust the applied displacement at the grip a little larger to preserve the SIF value on a constant (Fig. 1b). Therefore, it is natural to study the SIF and T-stress solution under two kinds of the boundary condition. The boundary collocation method (abbreviated as BCM) was widely used in the crack problems [6–8]. A particular advantage of this method is straightforward in derivation. The results obtained in the method depend on the collocation scheme used [6]. If the boundary collocation method (BCM) is used in the crack problems, one may be concerned about the influence of the used boundary collocation scheme to the computed result. For example, it was said, ‘‘the reason for the difference in stability of various boundary collocation procedure is not clear’’ [6]. A unified method of analysis is developed for asymmetric problems of a crack in a plate [7]. A set of general complex stress functions is proposed. Various problems wherever the crack is in the plate can be calculated by the method. For the particular cases, accurate values of the stress intensity factors (SIF) for the symmetric problems are obtained and they are used to show the validity of the method. For different eccentric and slant crack cases, the SIF have been calculated effectively. An edge crack problem in finite plate subjected to wedge forces is solved by the superposition of two solutions. In one solution of them, the boundary collocation method was used. It was shown that if sufficient terms are adopted in the expansion form, the stability in computation is not a problem for BCM any more [8]. The eigenfunction expansion variational method (abbreviated as EEVM) is developed to evaluate the SIFs and T-stress in this paper [9]. In the traction boundary value problem, EEVM is equivalent to the theorem of least potential energy in elasticity. Therefore, EEVM possesses a clear physical meaning. Meantime, EEVM does not need to design any boundary collocation scheme.

The following analysis depends on the complex variable function method in plane elasticity [10]. In the method, the stresses (rx, ry, rxy), the resultant forces (X,Y) and the displacements (u,v) are expressed in terms of complex potentials /(z), w(z) and x(z) such that rx þ ry ¼ 4ReUðzÞ ry  irxy ¼ 2ReUðzÞ þ zU0 ðzÞ þ WðzÞ ¼ UðzÞ þ ðz  zÞU0 ðzÞ þ XðzÞ

ð1Þ 0

0

f ¼ Y þ iX ¼ /ðzÞ þ z/ ðzÞ þ wðzÞ ¼ /ðzÞ þ ðz  zÞ/ ðzÞ þ xðzÞ

ð2Þ 0

0

2Gðu þ ivÞ ¼ j/ðzÞ  z/ ðzÞ  wðzÞ ¼ j/ðzÞ  ðz  zÞ/ ðzÞ  xðzÞ

ð3Þ

where z = x + iy denotes complex variable, G is the shear modulus of elasticity, j = (3  m)/(1 + m) is for the plane stress problems, j = 3  4m is for the plane strain problems, and m is the Poisson’s ratio. In the present study, the plane stress condition is assumed thoroughly and m = 0.3 is adopted. In Eqs. (1)–(3), the functions x(z), U(z), W(z) and X(z) are defined [10]  0 ðzÞ þ wðzÞ  xðzÞ ¼ z/ ð4Þ UðzÞ ¼ /0 ðzÞ;

WðzÞ ¼ w0 ðzÞ;

XðzÞ ¼ x0 ðzÞ

ð5Þ

 Note that, in Eq. (4) the function wðzÞ is an analytic func tion, which is defined by wðzÞ ¼ wðzÞ [10]. For a rectangular cracked plate, the two kinds of boundary condition can be formulated in the following manner. If the plate width 2b and height 2h is stretched in the y-direction with tension po, the traction boundary condition is indicated in Fig. 2a. For a perfect rectangular plate under the tension loading po, the strains (ex, ey, exy) and displacements (u, v) are as follows: mpo p ; ey ¼ o ; E E mpo x po y ; v¼ u¼ E E

ex ¼ 

exy ¼ 0

ð6Þ ð7Þ

where E (E = 2G(1 + m)) denotes Young’s modulus of elasticity.

Y.Z. Chen, X.Y. Lin / Computational Materials Science 42 (2008) 149–155

a

y

b

po

y* d B A

y* d x* B A

v = yvo / h

α

c o

2a

u uox / b v = vo

y

u = uo x*

c o

x

2h

2a

α x u = −u o v = yvo / h

2h 2b

151

2b u = −u o x / b v = −vo

c

y

po

A 2h

o 2a 2b

x

Fig. 2. (a) An inclined line crack in a rectangular cracked plate with the traction boundary condition, (b) an inclined line crack in a rectangular cracked plate with the assumed displacement boundary condition, (c) a central crack in a rectangular plate.

From Eq. (7), we can formulate the following assumed displacement boundary condition to the cracked plate (Fig. 2b) vo y h ðat the right edge x ¼ b; h 6 y 6 hÞ uo x ; v ¼ vo u¼ b ðat the upper edge  b 6 x 6 b; y ¼ hÞ vo y u ¼ uo ; v ¼ h ðat the left edge x ¼ b; h 6 y 6 hÞ uo x ; v ¼ vo u¼ b ðat the lower edge  b 6 x 6 b; y ¼ hÞ

u ¼ uo ;

v¼

rij nj ¼ pi ; ð8Þ

ð9Þ

ð10Þ

ðon cp Þ

where pi are the given fractions on the boundary cp, and nj denote the direction cosines. It is assumed that the complex potentials can be expressed in the expansion form /ðzÞ ¼

K X

X k /ðkÞ ðzÞ;

xðzÞ ¼

k¼1

a p h vo ¼ o E

ð13Þ

K X

X k xðkÞ ðzÞ

ð14Þ

k¼1

ð11Þ

where mp b uo ¼ o ; E

The eigenfunction expansion variational method (abbreviated as EEVM) is developed to solve a crack problem of finite region [9]. The cracked finite plate is shown in Fig. 3a. The traction boundary condition is assumed in the form

y

cp

ð12Þ

Clearly, for a perfect rectangular plate, the traction boundary condition shown by Fig. 2a and the assumed displacement boundary condition shown by Fig. 2b must provide the same elastic solution. However, for a cracked plate, two kinds of the boundary condition must give different solutions. This will be a topic in the present study.

b

y

cu

Imposed displacement

o 2a

x

o 2a

x

Fig. 3. (a) A finite cracked plate with the traction boundary condition, (b) a finite cracked plate with the assumed displacement boundary condition.

152

Y.Z. Chen, X.Y. Lin / Computational Materials Science 42 (2008) 149–155

where each pair /(k)(z) and x(k)(z) (k = 1, 2, . . ., K) satisfies the traction free condition along the crack faces and the single-valuedness condition of displacements. In Eq. (14), Xk(k = 1,2, . . ., K) are some undetermined coefficients. By using the variational principle in elasticity [9,11], the coefficients Xk(k = 1,2, . . ., K) can be determined by K X

x(z) = idM+1z, represents a rigid motion of plate. This term should be excluded in the summation. Meantime, the complex potentials are defined in the local coordinates ox*y* (Fig. 2). After the algebraic equation is solved, the stress intensity factors at the right crack tip A can be evaluated by pffiffiffiffiffiffiffiffiffiffiffi 1=2 ðK I  iK 2 ÞA ¼ 2ð2pÞ lim z  a/0 ðzÞ z!a

Amk X k ¼ Bm ;

ðm ¼ 1; 2; . . . ; KÞ

ð15Þ

M X ¼ 2ðpaÞ1=2 ðak þ ibk Þak1

k¼1

where

Z ðmÞ ðkÞ Amk ¼ Akm ¼ ui rij nj ds; ðm; k ¼ 1; 2; . . . ; KÞ cp Z ðmÞ ui  pi ds; ðm ¼ 1; 2; . . . ; KÞ Bm ¼

Similarly, at the left crack tip B we have ðK I  iK 2 ÞB ¼ 2ðpaÞ

ðkÞ (16), ui (k)

ðkÞ rij (k)

In Eq. and are derived from the complex potentials / (z) and x (z) using Eqs. (1)–(3). Similarly, the assumed displacement boundary condition is assumed in the form (Fig. 3b) ui ¼ ui ;

ðon cu Þ;

ð17Þ

Similarly, by using the variational principle in elasticity [9,11], the coefficients Xk(k = 1, 2, . . ., K) can be determined by Amk X k ¼ Bm ;

ðm ¼ 1; 2; . . . ; KÞ

ð18Þ

k¼1

where

Z

ðmÞ

ðkÞ

Amk ¼ Akm ¼ rij nj ui ds; ðm; k ¼ 1; 2; . . . ; KÞ cu Z ðmÞ Bm ¼ rij nj  ui ds; ðm ¼ 1; 2; . . . ; KÞ

1=2

ð19Þ

M X k1 ðak þ ibk ÞðaÞ

ð23Þ

k¼1

ð16Þ

cp

K X

ð22Þ

k¼1

Meantime, the T-stress at the cracks tip A and B can be evaluated by TA ¼ 4

M X

kck ak1 ;

TB ¼ 4

k¼1

M X

kck ðaÞk1

ð24Þ

k¼1

3. Numerical examples Example 1. In the first example (Fig. 2), we assume that the rectangular cracked plate with width 2b and height 2h has a traction boundary condition (Fig. 2a), or the assumed displacement boundary condition (Fig. 2b). Meantime, the center of crack is at the position with coordinates (c, d), and a denotes the inclined angle of crack. In computation, the inclined angle is subject to change from a = 0°, 10°, 20°, . . ., 90°. In the case of tension loading po, the stress intensity factors and the T-stresses at the crack tip A and B are expressed as

cu ðkÞ

ðkÞ

In Eq. (19), ui and rij are derived from the complex potentials /(k)(z) and x(k)(z) using Eqs. (1)–(3). In the cracked plate case, the complex potentials can be expressed in the expansion form 4M X k¼1

X k /ðkÞ ðzÞ;

xðzÞ ¼

4M X

X k xðkÞ ðzÞ

ð20Þ

k¼1

Two sets of appropriate expansion form in Eq. (20) can be expressed as follows [9]: pffiffiffiffiffiffiffiffiffiffiffiffiffiffi /ðkÞ ðzÞ ¼ xðkÞ ðzÞ ¼ ðak þ ibk Þ z2  a2 zk1 ; ð1 6 k 6 MÞ /ðkÞ ðzÞ ¼ xðkÞ ðzÞ ¼ ðck þ id k ÞzkM ;

ðM þ 1 6 k 6 2MÞ ð21Þ

It is easy to prove that the displacements, strains and stresses derived from the complex potentials /(k)(z) and x(k)(z) (k = 1, 2, . . ., 2M) satisfy all governing equations of plane elasticity and the traction free condition on the crack faces [9]. In Eq. (21), the coefficients ak, bk, ck, dk (k = 1, 2, . . ., M) stand for the coefficients Xk(k = 1, 2, . . ., 4M) in Eq. (20). Note that one complex potential in Eq. (21), or /(z) =

Non-dimensional SIF

/ðzÞ ¼

1.9764 2.0

h/b=1 c=0 d=0 a/b=0.8 F1A

1.5

1.0

0.5656

F2A

0.5

f1A f2A

0.0

0

20

40

60

80

α (degree) Fig. 4. Non-dimensional stress intensity factors F1A (a), F2A(a) (under the traction boundary condition) and f1A(a), f2A(a) (under the assumed displacement boundary condition), in the case of h/b = 1, c = d = 0 (see Fig. 1a and b and Eqs. (25) and (26)).

Y.Z. Chen, X.Y. Lin / Computational Materials Science 42 (2008) 149–155

pffiffiffiffiffiffi K 1A ¼ F 1A ðaÞpo pa; pffiffiffiffiffiffi K 1B ¼ F 1B ðaÞpo pa;

pffiffiffiffiffiffi K 2A ¼ F 2A ðaÞpo pa; T A ¼ GA ðaÞpo pffiffiffiffiffiffi K 2B ¼ F 2B ðaÞpo pa; T B ¼ GB ðaÞpo

ð25Þ

In the case of the assumed displacement boundary condition, the stress intensity factors and the T-stress at the crack tip A and B are expressed as Evo pffiffiffiffiffiffi pa; h Evo T A ¼ gA ðaÞ h Evo pffiffiffiffiffiffi pa; K 1B ¼ f1B ðaÞ h

K 1A ¼ f1A ðaÞ

K 2A ¼ f2A ðaÞ

K 2B ¼ f2B ðaÞ

Evo pffiffiffiffiffiffi pa; h

Evo pffiffiffiffiffiffi pa; h

T B ¼ gB ðaÞ

Evo h

153

In the case of the angle a being changed, three sets of computation, (1) c = d = 0, a/b = 0.8 (2) c/b = d/b = 0.25, a/b = 0.6, and (3) c/b = d/b = 05, a/b = 0.4 are presented. The computed results for SIFs and T-stresses are plotted in Figs. 4–9, respectively. From calculated results we see that the imposed boundary conditions can considerably influence the computed results. For example, in the case of c = d = 0, a/b = 0.8 and a = 0°, we have (1) F1A(a) = 1.9764 (for SIF), G(a) = 2.3359 (for T-stress) (for the imposed tension po), (2) f1A(a) = 0.5656 (for SIF), g(a) = 0.5304 (for T-stress) (for the assumed displacement boundary condition), respectively.

ð26Þ 1.0

1.5

h/b=1

1.0

c=0 d=0 a/b=0.8

0.5

gA

0.0

-0.5304

-0.5

GA

-1.0 -1.5

-2.3359

-2.0

h/b=1 c/b=d/b=0.25 a/b=0.6

0.5

Non-dimensional T-stress

Non-dimensional T-stress

where Evo/h = po.

gA

0.0

gB

-0.5

GB

-1.0 -1.5

GA

-2.0 -2.5

-2.5

0

0

20

40

60

Fig. 5. Non-dimensional T-stresses GA(a) (under the traction boundary condition) and gA(a) (under the assumed displacement boundary condition), in the case of h/b = 1 c = d = 0 (see Fig. 1a and b and Eqs. (25) and (26)).

40

60

80

α (degree)

80

α (degree)

Fig. 7. Non-dimensional T-stresses GA(a), GB(a) (under the traction boundary condition) and gA (a), gB(a) (under the assumed displacement boundary condition), in the case of h/b = 1 c/b = d/b = 0.25 (see Fig. 1a and b and Eqs. (25) and (26)).

2.0

2.0

1.7327

h/b=1

F1A

h/b=1

c/b=d/b=0.25 a/b=0.6

F1A

1.5

Non-dimensional SIF

1.5

Non-dimensional SIF

20

F1B 1.0

f1B f1A

0.5

F2A f2B

F2B

f2A

d/b=c/b=0.5 a/b=0.4

1.4855

F1B

1.0

0.7514

0.5

f1B

0.6621 f1A

F2A

F2B

f2B

0.0

f2A

0.0 -0.5 0

20

40

60

80

α (degree) Fig. 6. Non-dimensional stress intensity factors F1A(a), F2A(a), F1B(a), F2B (a) (under the traction boundary condition) and f1A(a), f2A(a), f1B(a), f2B(a) (under the assumed displacement boundary condition), in the case of h/b = 1, c/b = d/b = 0.25 (see Fig. 1a and b and Eqs. (25) and (26)).

0

20

40

60

80

α (degree)

Fig. 8. Non-dimensional stress intensity factors F1A(a), F2A(a), F1B(a), F2B(a) (under the traction boundary condition) and f1A(a), f2A(a), f1B(a), f2B(a) (under the assumed displacement boundary condition), in the case of h/b = 1, c/b = d/b = 0.5 (see Fig. 1a and b and Eqs. (25) and (26)).

154

Y.Z. Chen, X.Y. Lin / Computational Materials Science 42 (2008) 149–155

K 1A ¼ f1A ðh=b; a=bÞ

h/b=1

1.0

Non-dimensional T-stress

c/b=d/b=0.5 a/b=0.4

Evo pffiffiffiffiffiffi pa; h

T A ¼ gA ðh=b; a=bÞ

ð28Þ

0.5

where Evo/h = po.

gA

0.0

-0.4647 -0.5

-0.8477

gB

-1.0 -1.5

GB GA

-1.1232 -2.0079

-2.0 0

20

Evo h

40

60

80

α (degree)

Fig. 9. Non-dimensional T-stresses GA(a), GB(a) (under the traction boundary condition) and gA(a), gB(a) (under the assumed displacement boundary condition), in the case of h/b = 1 c/b = d/b = 0.5 (see Fig. 1a and b and Eqs. (25) and (26)).

Example 2. In the second example (Fig. 2c), we assume that the rectangular cracked plate with width 2b and height 2h has a traction boundary condition (refer to Fig. 2a), or the assumed displacement boundary condition (refer to Fig. 2b). Meantime, the central crack has a length 2a. In computation, the ratio a/b is subject to change from a/b = 0.1, 0.2, . . ., 0.9. In the case of tension loading po, the stress intensity factors and the T-stresses at the crack tip A and B are expressed as pffiffiffiffiffiffi K 1A ¼ F 1A ðh=b; a=bÞpo pa; T A ¼ GA ðh=b; a=bÞpo ð27Þ In the case of the assumed displacement boundary condition, the stress intensity factors and the T-stress at the crack tip A and B are expressed as

In the case of the ratio a/b being changed, two sets of computation, (1) h/b = 1 and (2) and h/b = 2 are presented. The computed results for SIFs and T-stresses are listed in Table 1. From calculated results we see that in the case of a/b being larger the imposed boundary conditions can considerably influence the computed results. For example, in the case of h/b = 1 and a/b = 0.6, we have (1) F1A(a) = 1.4805 (for SIF), G(a) = 1.5472 (for T-stress) (for the imposed tension po), (2)f1A(a) = 0.6857 (for SIF), g(a) = 0.7183 (for T-stress) (for the assumed displacement boundary condition), respectively. However, in the case of a/b being small the influence is minor. For example, in the case of h/b = 1 and a/b = 0.1, we have (1) F1A (a) = 1.0139 (for SIF), G(a) = 1.0188 (for T-stress) (for the imposed tension po), (2)f1A(a) = 0.9846 (for SIF), g(a) = 1.0072 (for T-stress) (for the assumed displacement boundary condition), respectively. 4. Conclusions Theoretically, the EEVM provides a reasonable solution, which can give a least potential energy for the cracked plate. However, in some other methods, one cannot judge if the least potential energy is achieved in the solution. It is also significant to use the computed results to the fatigue test of a cracked component. For example, if one wants to preserve the same level of SIF in the fatigue test of a cracked plate under tension, one must adjust the loading as the crack is growing.

Table 1 Non-dimensional stress intensity factors F1A(h/b,a/b), f1A(h/b,a/b) and T-stresses GA(h/b,a/b), gA(h/b,a/b) for a central crack in the rectangular plate (see Fig. 2c and Eqs. (27) and (28)) h/b = 1 case a/b 0.1 F1A 1.0139 f1A 0.9846 a/b 0.1 GA 1.0188 gA 0.9843

0.2 1.0554 0.9422 0.2 1.0726 0.9427

0.3 1.1232 0.8823 0.3 1.1554 0.8873

0.4 1.2160 0.8153 0.4 1.2600 0.8295

0.5 1.3337 0.7486 0.5 1.3849 0.7743

0.6 1.4805 0.6857 0.6 1.5472 0.7183

0.7 1.6750 0.6262 0.7 1.8104 0.6474

0.8 1.9764 0.5656 0.8 2.3359 0.5304

0.9 2.6273 0.4897 0.9 3.5888 0.2839

h/b = 2 case a/b 0.1 1.0061 F1A 0.9946 f1A a/b 0.1 GA 1.0072 gA 0.9940

0.2 1.0248 0.9786 0.2 1.0301 0.9758

0.3 1.0582 0.9527 0.3 1.0723 0.9452

0.4 1.1098 0.9178 0.4 1.1416 0.9009

0.5 1.1864 0.8748 0.5 1.2518 0.8402

0.6 1.3001 0.8238 0.6 1.4270 0.7573

0.7 1.4754 0.7639 0.7 1.7092 0.6414

0.8 1.7752 0.6912 0.8 2.1928 0.4687

0.9 2.4656 0.5944 0.9 3.3492 0.1747

Y.Z. Chen, X.Y. Lin / Computational Materials Science 42 (2008) 149–155

Acknowledgement The project was supported by National Natural Science Foundation of China. References [1] J.R. Rice, J. Mech. Phys. Solids 22 (1974) 17–26. [2] C. Betegon, J.W. Hancock, ASME J. Appl. Mech. 58 (1991) 104–110.

[3] [4] [5] [6] [7] [8] [9] [10] [11]

155

S. Melin, Int. J. Fract. 114 (2002) 259–265. T.L. Sham, Int. J. Fract. 48 (1991) 81–102. T. Fett, G. Rizzi, Int. J. Fract. 132 (2005) L9–L16. W.K. Wilson, J. Basic Eng. 93 (1971) 685–690. Y. Wang, Eng. Fract. Mech. 40 (1991) 133–143. Q.Z. Xiao, B.L. Karihaloo, Eng. Frac. Mech. 69 (2002) 959–981. Y.Z. Chen, Eng. Fract. Mech. 17 (1983) 387–394. N.I. Muskhelishvili, Noordhooff, Groningen, 1953. W.C. Wang, J.T. Chen, Com. Meth. Appl. Mech. Eng. 73 (1989) 153– 171.