Dynamic crack tip fields and dynamic crack propagation characteristics of anisotropic material

Theoretical and Applied Fracture Mechanics 51 (2009) 73–85

Contents lists available at ScienceDirect

Theoretical and Applied Fracture Mechanics journal homepage: www.elsevier.com/locate/tafmec

Dynamic crack tip fields and dynamic crack propagation characteristics of anisotropic material X. Gao a,*, X.W. Kang a,b, H.G. Wang a a b

501 Section of Xi’an High-Tech Institute, Xi’an 710025, PR China MOE Key Laboratory for Strength and Vibration, Xi’an Jiaotong University, Xi’an 710049, PR China

a r t i c l e

i n f o

Article history: Available online 23 January 2009 Keywords: Fracture dynamics Stress intensity factor Crack propagation Strain energy density Anisotropy

a b s t r a c t Based on mechanics of anisotropic material, the dynamic crack propagation problem of I/II mixed mode crack in an infinite anisotropic body is investigated. Expressions of dynamic stress intensity factors for modes I and II crack are obtained. Components of dynamic stress and dynamic displacements around the crack tip are derived. The strain energy density theory is used to predict the dynamic crack extension angle. The critical strain energy density is determined by the strength parameters of anisotropic materials. The obtained dynamic crack tip fields are unified and applicable to the analysis of the crack tip fields of anisotropic material, orthotropic material and isotropic material under dynamic or static load. The obtained results show Crack propagation characteristics are represented by the mechanical properties of anisotropic material, i.e., crack propagation velocity M and fiber direction a. In particular, the fiber direction a and the crack propagation velocity M give greater influence on the variations of the stress fields and displacement fields. Fracture angle is found to depend not only on the crack propagation but also on the anisotropic character of the material. Ó 2009 Elsevier Ltd. All rights reserved.

1. Introduction The dynamic fracture problems by inertia force have been broadly classified as stationary and propagating crack. Knowledge of the dynamic crack tip fields and dynamic crack propagation characteristic are important to understanding the process and nature of dynamic fracture. Several problems have been examined in dynamic fracture mechanics for isotropic materials. The earliest theoretical work about dynamic fracture mechanics was provided in [1,2]. They presented the relationship between displacement components and the crack propagation velocity, and the relationship between stress components and the crack propagation velocity. Subsequently, the studies [3–5] considered the stress and displacement fields of dynamic crack problems in isotropic material. Considered in [6,7] are the moving crack problems of orthotropic materials. These investigations solved the Griffith crack problem by using integral transform techniques. In dynamic fracture of orthotropic material, refer to the studies in [8–10] are the stress and displacement fields of a propagating crack tip with constant velocity. But their dynamic stress and displacement components of a propagating crack tip in orthotropic material do not clearly represent the stress intensity factor. Studied in [11] is the propagating crack problems of ortho* Corresponding author. E-mail address: [email protected] (X. Gao). 0167-8442/$ - see front matter Ó 2009 Elsevier Ltd. All rights reserved. doi:10.1016/j.tafmec.2009.01.006

tropic material under the dynamic plane mode. He derived the dynamic stress components and dynamic displacement components around the crack tip of orthotropic material with propagating crack under the dynamic load and the steady state, and studied the crack propagation characteristics of orthotropic material. Studied in [12] is the dynamic propagation problem of mode III crack in an orthotropic infinite body and presented the simple approximate expressions of the stress components and displacement components around crack tip. Presented in [13] is a method for analyzing the elastodynamic response of an infinite orthotropic material with a semi-infinite crack under impact loads. A closed form solution for the dynamic stress intensity factor is obtained, and he found expressions for the stress ahead of the crack tip and the displacements on the crack faces. Studied in [14] is the propagating crack problem of an anisotropic material. The dynamic crack tip fields were obtained, and the crack propagation characteristics of the anisotropic materials are also studied. These former works were important. Unfortunately, their works only dealt with the dynamic crack problem of isotropic and orthotropic material and did not investigate the direction of crack propagation. More than two decades ago, the volume energy density was proposed [15] as a fracture criterion in contrast to the idea of Griffith’s energy release rate. This provided an alternative approach to failure prediction for the same stress solution. The distinctions were emphasized in the series Mechanics of Fracture [16]. The theory requires no calculation on the energy release rate and thus

74

X. Gao et al. / Theoretical and Applied Fracture Mechanics 51 (2009) 73–85

possesses the inherent advantage of being able to treat all mixed mode crack extension problems. Unlike the conventional theory of G and K which measures only the amplitude of the local stress, the fundamental parameter in this theory, the strain energy density factor S which is defined as the coefficient of 1/r singular behavior of the volume energy density dW/dV, is also direction sensitive. The stationary values of the density factor S can predict the direction of crack growth in any condition. The direction of crack initiation is assumed to occur when Smin reaches a critical value Scri that can be used as an intrinsic materials parameter and is independent of the crack geometry and loading. This theory gained momentum and credibility in engineering. A review on the use of the volume energy density can be found in [17,18]. Recently, Shen and Nishioka [19] used this theory to develop a fracture criterion for piezoelectric materials. Their theoretical result agrees qualitatively with the empirical evidence by assuming impermeable crack. In the present work, a method for analyzing the propagating crack problem of anisotropic plate under the dynamic plane mode is presented based on mechanics of anisotropic material. The expressions of dynamic stress intensity factors for modes I and II crack are obtained, the dynamic stress fields and dynamic displacement fields around the crack tip are derived, and the energy density theory is applied as a fracture criterion to study the direction of crack propagation.

ð1Þ

where t is time, q is the mass density of the material, and U and V are the displacement components along the X and Y-axis, respectively. When the crack is moving at velocity c along the X-axis. One introduce the moving coordinate system oxy (see Fig. 1), it is convenient to introduce the transformation

y ¼ Y:

ð2Þ

Therefore, Eq. (1) can be transformed into

@ rx @ sxy @2U þ ¼ qc 2 2 ; @x @x @y

@ sxy @ ry @2V þ ¼ qc 2 2 : @x @x @y

Substituting Eqs. (4) and (5) into Eq. (3), Eq. (3) becomes

@ rx @ ry @ sxy @ rx @ ry @ sxy A þB þC ¼ 0; D þE þF ¼ 0; @z @z @z @z @z @z 2 2 2 A ¼ 1  qc a11 ; B ¼ qc a12 ; C ¼ m  qc a16 ; D ¼ qc2 a12 =m;

E ¼ m  qc2 a22 =m;

F ¼ 1  qc2 a26 =m:

ð6Þ

Using the Airy stress function w in Eq. (6), stress components can be obtained as follows:

@2w @2w @2w ; ry ¼ fy 2 ; sxy ¼ fxy 2 ; 2 @z @z  @z  f x ¼ m3  mqc2 ða22 þ ma16  a12 Þ þ q2 c4 ða16 a22  a12 a26 Þ =m;   f y ¼ m  qc2 ða26 þ ma11  ma12 Þ þ q2 c4 ða11 a26  a12 a16 Þ =m;   ð7Þ f xy ¼ m2 þ qc2 ða22 þ m2 a11 Þ  q2 c4 ða11 a22  a212 Þ =m:

rx ¼ fx

The strain-compatibility equation is 2 @ ex @ 2 ey @ cxy þ ¼ : @y2 @x2 @x@y 2

ð8Þ

Substituting Eqs. (4) and (7) into Eq. (8), yields

@4w ¼ 0; @z4

a ¼ 2a16 =a11 ;   b ¼ 2a12 þ a66 þ qc2 ða212 þ a216  a11 a22  a11 a66 Þ =a11 ;  2  c ¼ 2 qc ða22 a16 þ a11 a26  a12 a16  a12 a26 Þ  a26 =a11 ;  d ¼ a22 þ qc2 ða212 þ a226  a11 a22  a22 a66 Þ

 þq2 c4 ð2a12 a16 a26 þ a11 a22 a66  a11 a226  a22 a216  a66 a212 Þ =a11 : ð9Þ

The solution of Eq. (9) is

wðzÞ ¼ w1 ðz1 Þ þ w2 ðz2 Þ þ w1 ðz1 Þ þ w2 ðz2 Þ ¼ 2Re½w1 ðz1 Þ þ w2 ðz2 Þ;

ð10Þ

where z1 = x + m1y, z2 = x + m2y. Here m1 and m2 are the roots of characteristic Eq. (11)

ð3Þ

Under plane strain or plane stress conditions, the relations of the strain and stress are

ex ¼ a11 rx þ a12 ry þ a16 sxy ; ey ¼ a12 rx þ a22 ry þ a26 sxy ; cxy ¼ a16 rx þ a26 ry þ a66 sxy :

ð5Þ

2

The equilibrium equations of stress with inertia force in fixed Cartesian coordinate system OXY are

x ¼ X  ct;

z ¼ x þ my:

ðm4 þ am3 þ bm þ cm þ dÞ

2. Stress and displacement fields of propagating crack

@ rx @ sxy @ 2 U @ sxy @ ry @2V þ ¼q 2 þ ¼q 2 @X @Y @Y @t @X @t

The complex variable z is

2

m4 þ am3 þ bm þ cm þ d ¼ 0: Set

/ðz1 Þ ¼ ð4Þ

ð11Þ

dw1 ; dz1

/0 ðz1 Þ ¼

Fig. 1. Notation for bipolar coordinates (r, h), (r1, h1) and (r2, h2).

d/ ; dz1

uðz2 Þ ¼

dw2 ; dz2

u0 ðz2 Þ ¼

du : dz2

ð12Þ

75

X. Gao et al. / Theoretical and Applied Fracture Mechanics 51 (2009) 73–85

Substituting Eqs. (10) and (12) into Eq. (7), it is found that

rx ¼ 2Re½fx ðm1 Þ/0 þ fx ðm2 Þu0 ; ry ¼ 2Re½fy ðm1 Þ/0 þ fy ðm2 Þu0 ; sxy ¼ 2Re½fxy ðm1 Þ/0 þ fxy ðm2 Þu0 : ð13Þ The angle h between the crack tip and crack plane is p and p, respectively, because the crack plane is discontinuous in the upper and lower parts. In such a case, h = h1 = h2 = ±p and r = r1 = r2. The normal and shear stress on the crack plane do not exist, that is

ry ¼ 0; sxy ¼ 0: 0

pffiffiffiffiffiffiffiffiffi K I ¼ lim 2pr ry ðh ¼ 0Þ; r!0

0

D1 / þ D2 / þ D3 u þ D4 u ¼ 0; 0

0

0

0

D5 / þ D6 / þ D7 u þ D8 u ¼ 0; 0

D1 ¼ fy ðm1 Þ; D5 ¼ fxy ðm1 Þ;

D2 ¼ fy ðm1 Þ;

D3 ¼ fy ðm2 Þ;

D6 ¼ fxy ðm1 Þ;

D7 ¼ fxy ðm2 Þ;

pffiffiffiffiffiffiffi pffiffiffiffiffiffiffi 2pRe½D1 A1 þ D3 B1  ¼ 2pðD1 A1 þ D3 B1 Þ; pffiffiffiffiffiffiffi pffiffiffiffiffiffiffi K II ¼ 2pRe½D5 A1 þ D7 B1  ¼ 2pðD5 A1 þ D7 B1 Þ:

D8 ¼ fxy ðm2 Þ:

1 D7 K I  D3 K II A1 ¼ pffiffiffiffiffiffiffi ; 2p D1 D7  D3 D5

/ðz1 Þ ¼

uðz2 Þ ¼

X

An zk1n ¼

n X

Bn zk2n ¼

n

X n X

Bn r k2n eikn h2 ;

ð17Þ

n

ð18aÞ

(2) h1 = h2 = p

ðD1 An þ D3 Bn Þeiðkn 1Þp þ ðD2 An þ D4 Bn Þeiðkn 1Þp ¼ 0; ð18bÞ

To get the solution of Eq. (18) regardless of (D1An + D3Bn) and ðD2 An þ D4 Bn Þ, and (D5An + D7Bn) and ðD6 An þ D8 Bn Þ, the coefficient matrix of Eq. (18) should be zero, that is

 eiðkn 1Þp   ¼ 0; eiðkn 1Þp 

sin 2ðkn  1Þp ¼ 0; n kn ¼ ; ðn ¼ 1; 2; 3 . . .Þ; 2

ð26Þ

 n2 n D7 K nI  D3 K nII pffiffiffiffiffiffiffi ½rðcos h þ m1 sin hÞ 2 ; 2 2p D1 D7  D3 D5 n  X  n D1 K n  D5 K n n2 II I pffiffiffiffiffiffiffi u0 ðz2 Þ ¼ ½rðcos h þ m2 sin hÞ 2 : 2 2p D1 D7  D3 D5 n

ð27aÞ

X

ð27bÞ

When n = 1, Eq. (27) becomes

ðD1 An þ D3 Bn Þeiðkn 1Þp þ ðD2 An þ D4 Bn Þeiðkn 1Þp ¼ 0;

  eiðkn 1Þp   iðk 1Þp e n

1 D1 K nII  D5 K nI Bn ¼ pffiffiffiffiffiffiffi : 2p D1 D7  D3 D5

n 1 X D7 K nI  D3 K nII /ðz1 Þ ¼ pffiffiffiffiffiffiffi ½rðcos h þ m1 sin hÞ2 ; 2p n D1 D7  D3 D5

(1) h1 = h2 = p

ðD5 An þ D7 Bn Þeiðkn 1Þp þ ðD6 An þ D8 Bn Þeiðkn 1Þp ¼ 0:

ð25Þ

From Eqs. (17), (23) and (26), we get

/0 ðz1 Þ ¼

where An and Bn are complex constants, respectively, kn are real eigenvalues which will be determined later. Substituting Eq. (17) into Eq. (15), there results

ðD5 An þ D7 Bn Þeiðkn 1Þp þ ðD6 An þ D8 Bn Þeiðkn 1Þp ¼ 0;

1 D1 K II  D5 K I B1 ¼ pffiffiffiffiffiffiffi : 2p D1 D7  D3 D5

n 1 X D1 K nII  D5 K nI uðz2 Þ ¼ pffiffiffiffiffiffiffi ½rðcos h þ m2 sin hÞ2 ; 2p n D 1 D 7  D 3 D 5

An r k1n eikn h1 ;

ð24Þ

For general complex constants An and Bn, we have

ð16Þ

Analytical complex functions /(z1) and u(z2) can be expressed as power series as follows:

ð23Þ

From Eq. (24), we get

1 D7 K nI  D3 K nII An ¼ pffiffiffiffiffiffiffi ; 2p D1 D7  D3 D5

D4 ¼ fy ðm2 Þ;

r!0

KI ¼

ð15Þ where

pffiffiffiffiffiffiffiffiffi K II ¼ lim 2pr sxy ðh ¼ 0Þ:

When n = 1, by substituting Eqs. (13), (17) and (22) into Eq. (23), the relationships between the complex constants and stress intensity factor are

ð14Þ

From Eqs. (13) and (14), there results 0

The dynamic stress intensity factors are defined as

ð19Þ ð20Þ ð21Þ

where kn is positive, otherwise it has not physical meaning. Substituting Eq. (21) into Eq. (18), we get:

1 1 D7 K I  D3 K II /ðz1 Þ ¼ pffiffiffiffiffiffiffi ½rðcos h þ m1 sin hÞ2 ; 2p D 1 D 7  D 3 D 5 1 1 D K  D5 K I uðz2 Þ ¼ pffiffiffiffiffiffiffi 1 II ½rðcos h þ m2 sin hÞ2 ; 2p D1 D7  D3 D5

ð28aÞ

1 1 D7 K I  D3 K II ½rðcos h þ m1 sin hÞ2 ; /0 ðz1 Þ ¼ pffiffiffiffiffiffiffi 2 2p D1 D7  D3 D5 1 1 D K  D5 K I u0 ðz2 Þ ¼ pffiffiffiffiffiffiffi 1 II ½rðcos h þ m2 sin hÞ2 : 2 2p D 1 D 7  D 3 D 5

ð28bÞ

From Eq. (4), we get

@U ¼ a11 rx þ a12 ry þ a16 sxy ; @x @V ¼ a12 rx þ a22 ry þ a26 sxy : ¼ @y

ex ¼ ey

ð29Þ

The dynamic displacement fields can be obtained from Eqs. (13) and (29) as follows:

U ¼ 2Re½p1 /ðz1 Þ þ p2 /ðz2 Þ;

V ¼ 2Re½q1 /ðz1 Þ þ q2 /ðz2 Þ;

p1 ¼ a11 fx ðm1 Þ þ a12 fy ðm1 Þ þ a16 fxy ðm1 Þ; (1) When n is odd number

p2 ¼ a11 fx ðm2 Þ þ a12 fy ðm2 Þ þ a16 fxy ðm2 Þ;

D1 An þ D3 Bn ¼ D2 An þ D4 Bn ¼ D1 An þ D3 Bn ; D5 An þ D7 Bn ¼ D6 An þ D8 Bn ¼ D5 An þ D7 Bn :

q1 ¼ ½a12 fx ðm1 Þ þ a22 fy ðm1 Þ þ a26 fxy ðm1 Þ=m1 ; ð22aÞ

(2) When n is even number

D1 An þ D3 Bn ¼ ðD2 An þ D4 Bn Þ ¼ ðD1 An þ D3 Bn Þ; D5 An þ D7 Bn ¼ ðD6 An þ D8 Bn Þ ¼ ðD5 An þ D7 Bn Þ:

ð22bÞ

q2 ¼ ½a12 fx ðm2 Þ þ a22 fy ðm2 Þ þ a26 fxy ðm2 Þ=m2 :

ð30Þ

Substituting Eq. (27) into Eqs. (13) and (30), the general mixed mode displacement components and stress components of anisotropic material for constant velocity propagating crack can be represented as

76

X. Gao et al. / Theoretical and Applied Fracture Mechanics 51 (2009) 73–85

"

# n p1 X D7 K nI  D3 K nII 2 U ¼ 2Re pffiffiffiffiffiffiffi ½rðcos h þ m1 sin hÞ 2p n D1 D7  D3 D5 " # n p2 X D1 K nII  D5 K nI 2 þ 2Re pffiffiffiffiffiffiffi ½rðcos h þ m2 sin hÞ ; 2p n D1 D7  D3 D5 " # n q1 X D7 K nI  D3 K nII V ¼ 2Re pffiffiffiffiffiffiffi ½rðcos h þ m1 sin hÞ2 2p n D1 D7  D3 D5 " # n q2 X D1 K nII  D5 K nI ffi þ 2Re pffiffiffiffiffiffi ½rðcos h þ m2 sin hÞ2 ; 2p n D1 D7  D3 D5

ð31Þ

" # o X n2 D7 K nI  D3 K nII ½rðcos h þ m1 sin hÞ 2 nfx ðm1 Þ D1 D7  D3 D5 n " #  o X n2 1 D1 K nII  D5 K nI 2 þ pffiffiffiffiffiffiffi Re ; ½rðcos h þ m2 sin hÞ nfx ðm2 Þ D1 D7  D3 D5 2p n " # o X n2 1 D7 K nI  D3 K nII ½rðcos h þ m1 sin hÞ 2 ry ¼ pffiffiffiffiffiffiffi Re nfy ðm1 Þ D1 D7  D3 D5 2p n " #  o X n2 1 D1 K nII  D5 K nI 2 þ pffiffiffiffiffiffiffi Re ; nfy ðm2 Þ ½rðcos h þ m2 sin hÞ D1 D7  D3 D5 2p n " # o X n2 1 D7 K nI  D3 K nII sxy ¼ pffiffiffiffiffiffiffi Re nfxy ðm1 Þ ½r ðcos h þ m1 sin hÞ 2 D D  D D 1 7 3 5 2p n " #  n o X n2 1 D1 K II  D5 K nI þ pffiffiffiffiffiffiffi Re : nfxy ðm2 Þ ½rðcos h þ m2 sin hÞ 2 D1 D7  D3 D5 2p n 1

rx ¼ pffiffiffiffiffiffiffi Re 2p

ð32Þ When n < 0 and r approaches the crack tip in Eq. (31), the displacement is infinite, thus negative eigenvalues are not suitable. When n = 0, the stress components of Eq. (32) are zero, and the displacement components of Eq. (31) represent rigid displacement. Therefore the eigenvalues n are greater than 0 and are integer in the stress, strain and displacement components of elastic body with a crack. Substituting Eq. (28) into Eqs. (13) and (30), the mixed mode displacement components and stress components around crack tip for constant velocity propagating crack can be represented as

rffiffiffiffiffi 2r

Re

 pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 p D cos h þ m1 sin h p D1 D7  D3 D5 1 7 pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi  p2 D5 cos h þ m2 sin h rffiffiffiffiffi  pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2r 1 Re  K II p D cos h þ m1 sin h p D1 D7  D3 D5 1 3 pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi  p2 D1 cos h þ m2 sin h ; rffiffiffiffiffi  pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2r 1 Re q1 D7 cos h þ m1 sin h V ¼ KI p D1 D7  D3 D5 pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi  q2 D5 cos h þ m2 sin h rffiffiffiffiffi  pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2r 1 Re q1 D3 cos h þ m1 sin h  K II p D1 D7  D3 D5 pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi  q2 D1 cos h þ m2 sin h ; U ¼ KI

rx

" KI fx ðm1 Þ D7 ¼ pffiffiffiffiffiffiffiffiffi Re pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2pr cos h þ m1 sin h D1 D7  D3 D5 # fx ðm2 Þ D5  pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi cos h þ m2 sin h D1 D7  D3 D5

" K II fx ðm1 Þ D3  pffiffiffiffiffiffiffiffiffi Re pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi D D 2pr cos h þ m1 sin h 1 7  D3 D5 # fx ðm2 Þ D1  pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ; cos h þ m2 sin h D1 D7  D3 D5 " KI fy ðm1 Þ D7 ffi Re pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ry ¼ pffiffiffiffiffiffiffiffi 2pr cos h þ m1 sin h D1 D7  D3 D5 # fy ðm2 Þ D5  pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi cos h þ m2 sin h D1 D7  D3 D5 " K II fy ðm1 Þ D3  pffiffiffiffiffiffiffiffiffi Re pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2pr cos h þ m1 sin h D1 D7  D3 D5 # fy ðm2 Þ D1  pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ; cos h þ m2 sin h D1 D7  D3 D5 " KI fxy ðm1 Þ D7 p ffiffiffiffiffiffiffiffi ffi Re pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi sxy ¼ 2pr cos h þ m1 sin h D1 D7  D3 D5 # fxy ðm2 Þ D5  pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi cos h þ m2 sin h D1 D7  D3 D5 " K II fxy ðm1 Þ D3  pffiffiffiffiffiffiffiffiffi Re pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi D D 2pr cos h þ m1 sin h 1 7  D3 D5 # fxy ðm2 Þ D1  pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi : cos h þ m2 sin h D1 D7  D3 D5

ð34Þ

When KII = 0 or KI = 0, Eqs. (33) and (34) are mode I or II displacement fields and stress fields around the crack tip of anisotropic material under dynamic load. When a16 = a26 = 0, Eqs. (33) and (34) coincide with the displacement and stress components of Lee [8] which are those of the crack tip for orthotropic material. When c = 0, Eqs. (33) and (34) coincide with the displacement and stress components of Sih et al. which are those of the crack tip for orthotropic material under the static load [20]. Modifying the constants aij of Eq. (4), Eqs. (33) and (34) are the displacement and stress components around the crack tip of isotropic material under dynamic load, respectively. Similarity, The dynamic strain fields can be obtained from Eqs. (4), (32) and (34). 3. Crack propagation characteristics Carbon fiber epoxy composites T300/5208 is used to survey the characteristics of stress and displacement distribution at the propagating crack tip. Table 1 shows the mechanical properties of T300/ 5208. The influence of crack propagation velocity on stress and displacement components of plane stress is studied when the fiber orientation angle a = 0, p/6, p/3, p/2, respectively. The crack propagation velocity used in this study is M = c(q a66)1/2. The velocity of shear stress wave velocity is Cs = (1/qa66)1/2, and M is c/Cs. Because the crack propagation velocity c cannot be faster than the shear stress wave velocity Cs, therefore M < 1, and here M = 0, 0.7, 0.9. 3.1. ModeI crack

ð33Þ

Fig. 2 shows the relationship between rx =K I and M with h when a = 0, p/6, p/3, p/2. The rx =K I value of a propagating crack is much greater than that of a stationary crack, but the inclinations of stress

Table 1 Mechanical properties of T300/5208. Material

E1 (GPa)

E2 (GPa)

G12 (GPa)

v21

X (MPa)

Y (MPa)

S (MPa)

T300/ 5208

181

10.3

7.17

0.28

1347

68.9

91

77

X. Gao et al. / Theoretical and Applied Fracture Mechanics 51 (2009) 73–85

Fig. 2. Mode I

rx =K I versus angle h.

distribution at crack tip are almost the same. Distributions of rx =K I vary with a but the inclinations of stress distribution are almost the same. The faster the crack propagation velocity is, the greater the stress value of rx =K I around crack tip for any fiber direction

Fig. 3. Mode I

ry =K I versus angle h.

a is. The maximum of rx =K I occurs at h = 0° when a = 0. As a increases, the maximum value of

rx =K I decreases gradually when

78

X. Gao et al. / Theoretical and Applied Fracture Mechanics 51 (2009) 73–85

M = 0, 0.7. But when M = 0.9, the value rx =K I increases gradually with the increment of a, the maximum of the value occurs at a = p/6, then the value rx =K I decreases gradually with the increment of a, the minimum of the value occurs at a = p/2. When a = 0, p/2, the stress value of rx =K I has a symmetrical distribution to the crack plane at any velocity of crack propagation. Fig. 3 shows the relationship between ry =K I and M with h when a = 0, p/6, p/3, p/2. The faster the crack propagation velocity is, the greater the stress value of ry =K I around crack tip for any fiber direction a is. As a increases, the stress value ry =K I increases gradually when M = 0, 0.7, the maximum of ry =K I occurs at a = p/2. When M = 0.9, the stress value ry =K I increases gradually with the increment of a, the maximum of ry =K I occurs at a = p/3; then the stress value ry =K I decreases gradually with the increment of a, the minimum of ry =K I occurs at a = 0. The variations of ry =K I with h are greater when a = p/6, p/3 than when a = 0, p/2. When a = 0, p/2, the stress value of ry =K I has a symmetrical distribution to the crack plane at any velocity of crack propagation. Distributions of ry =K I vary with a but the inclinations of stress distribution are almost the same. The maximum of ry =K I occurs at h = ±90° when a = p/2. Fig. 4 shows the relationship between sxy =K I and M with h when a = 0, p/6, p/3, p/2. The faster the crack propagation velocity is, the greater the stress value of sxy =K I around crack tip is. For any crack propagation velocity M, the stress value sxy =K I increases gradually with the increment of a, the maximum of sxy =K I occurs at a = p/6; then the stress value sxy =K I decreases gradually with the increment of a, the minimum of sxy =K I occurs at a = p/2. Distributions of sxy =K I vary with a but the inclinations of stress distribution are almost the same. When a = 0, p/2, the stress value of sxy =K I has a dissymmetrical distribution to the crack plane at any velocity of crack propagation. The maximum of sxy =K I occurs at h = ±90° when a = p/2. It is known from above results that stresses rx, ry, sxy around the propagation crack tip increase with the increment of crack propagation velocity and depend on the fiber direction angle a. The stress value of rx and sxy are minimums when a = p/2. The stress value of ry is minimum when the direction of the stress component is normal to the fiber direction (a = 0). The maximum of sxy at a = 0 are greater than the maximum of sxy at a = p/2. When M increases from 0 to 0.9, the magnitude order influenced by the velocity of crack propagation among rx, ry and sxy is rx, sxy, ry in any fiber directions. When the fiber direction angle a is changed from 0 to p/2, the order influenced by the fiber direction among rx, ry and sxy is rx, ry, sxy at any velocity of crack propagation. Fig. 5 shows the relationship between U=K I and M with h when a = 0, p/6, p/3, p/2. U=K I of a = 0 is greater than that of a = p/2 when M = 0 and M = 0.7, but U=K I of a = p/2 is greater than that of a = 0 when M = 0.9, and the maximum values of U=K I , respectively, occur at h = ±85° when a = 0 and h = 0° when a = p/2. The faster the crack propagation velocity is, the greater the displacement value of U=K I around crack tip is. For any crack propagation velocity M, the displacement value U=K I increases gradually with the increment of a, the maximum of U=K I occurs at a = p/6; then the displacement value U=K I decreases gradually with the increment of a, the minimum of U=K I occurs at a = p/2. Distributions of U=K I vary with a but the inclinations of displacement distribution are almost the same. When a = 0, p/2, the displacement value of U=K I has a symmetrical distribution to the crack plane at any velocity of crack propagation. Fig. 6 shows the relationship between V=K I and M with h when a = 0, p/6, p/3, p/2. The faster the crack propagation velocity is, the greater the displacement value of V=K I around crack tip is. For any crack propagation velocity M, the displacement value V=K I increases gradually with the increment of a, the maximum of V=K I

Fig. 4. Mode I

sxy =K I versus angle h.

occurs at a = p/6; then the displacement value V=K I decreases gradually with the increment of a, the minimum of V=K I occurs

X. Gao et al. / Theoretical and Applied Fracture Mechanics 51 (2009) 73–85

at a = p/2. The maximum of V=K I occurs at h = ±180° for any M and a. The influence of M and a on the displacement value of V=K I is great at any M and a when h 2 [180°, 90°] [ [90°, 180°], on the contrary, the influence is slight when h 2 [90°, 90°] at any M and a. When a = 0, p/2, the displacement value of V=K I has a dissymmetrical distribution to the crack plane at any velocity of crack propagation. We know from Figs. 5 and 6 that displacements U and V around the propagation crack tip increase with the increment of crack propagation velocity M and depend on the fiber direction angle a. Displacement V, which is normal to the crack direction at the crack tip, is zero at h = 0° and the maximums of V occur at h = ±180° at any crack propagation velocity. The crack propagation velocity gives greater influence on the displacement U than on the displacement V when M increases from 0 to 0.9. When the fiber direction angle a is changed from 0 to p/2, the order influenced by the fiber direction among U and V is U, V. 3.2. Mode II crack Fig. 7 shows the relationship between rx =K II and M with h when a = 0, p/6, p/3, p/2. The faster the crack propagation velocity is, the greater the stress value of rx =K II around crack tip for any fiber direction a is. As a increases, the stress value rx =K II decreases gradually when M = 0, 0.7. When M = 0.9, the stress value rx =K II increases gradually with the increment of a, the maximum of rx =K II occurs at a = p/6; then the stress value rx =K II decreases gradually with the increment of a, the minimum of rx =K II occurs at a = p/2. When a = 0, p/2, the stress value of rx =K II has a dissymmetrical distribution to the crack plane and the maximum of rx =K II occurs at h = ±180° at any velocity of crack propagation. Fig. 8 shows the relationship between ry =K II and M with h when a = 0, p/6, p/3, p/2. The faster the crack propagation velocity is, the greater the stress value of ry =K II around crack tip for any fiber direction a is. As a increases, the stress value ry =K II increases gradually when M = 0, 0.7, the maximum of ry =K II occurs at a = p/2. When M = 0.9, the stress value ry =K II increases gradually with the increment of a, the maximum of ry =K II occurs at a = p/3; then the stress value ry =K II decreases gradually with the increment of a, the minimum of ry =K II occurs at a = 0. The variations of ry =K II with h are greater when a = p/6, p/3 than when a = 0, p/2. For any M and a, the stress value of ry =K II are zero at h = ±180°. When a = 0, p/2, the stress value of ry =K II has a dissymmetrical distribution to the crack plane at any velocity of crack propagation. Fig. 9 shows the relationship between sxy =K II and M with h when a = 0, p/6, p/3, p/2. The maximum stress value of sxy =K II occurs at h = 0 at any crack propagation velocity M when a = 0. For any M and a, the stress value of ry =K II are zero at h = ±180°. Distributions of sxy =K II vary with a but the inclinations of stress distribution are almost the same. When a = 0, p/2, the stress value of sxy =K II has a symmetrical distribution to the crack plane at any velocity of crack propagation. It is known from above results that stresses rx, ry, sxy around the propagation crack tip increase with the increment of crack propagation velocity and depend on the fiber direction angle a. The stress value of rx is minimum when the direction of the stress component is normal to the fiber direction (a = p/2). The stress value of ry and sxy are minimums when a = 0. Fig. 10 shows the relationship between U=K II and M with h when a = 0, p/6, p/3, p/2. The faster the crack propagation velocity is, the greater the displacement value of U=K II is. For any M and a, the displacement value of U=K II is zero at h = 0°. When a = 0, p/2, the displacement value of U=K II has a dissymmetrical distribution to the crack plane at any velocity of crack propagation. The displacement value of U=K II is greater when a = p/2 than when a = 0 at any crack propagation velocity M.

Fig. 5. Mode I U=K I versus angle h.

79

80

X. Gao et al. / Theoretical and Applied Fracture Mechanics 51 (2009) 73–85

Fig. 11 shows the relationship between V=K II and M with h when a = 0, p/6, p/3, p/2. The faster the crack propagation velocity is, the greater the displacement value of V=K II around crack tip is. The influence of M on the displacement value of V=K II is slight when a = p /2 and h 2 [180°, 90°] [ [90°, 180°], on the contrary, the influence is obvious when h 2 [90°, 90°] and a = p/2. When a = 0, p/2, the displacement value of V=K II has a symmetrical distribution to the crack plane at any velocity of crack propagation. The maximum displacement value of V=K II occurs at h = 0 at any crack propagation velocity M when a = 0, p/2. It is known from Figs. 10 and 11 that displacements U and V around the propagation crack tip increase with the increment of crack propagation velocity M and depend on the fiber direction angle a. Displacement V is zero at h = ±180° and the maximums of V occur at h = 0° at any velocity of crack propagation when a = 0, p/2. The crack propagation velocity gives greater influence on the displacement V than on the displacement U when M increases from 0 to 0.9. When the fiber direction angle a is changed from 0 to p/ 2, the order influenced by the fiber direction among U and V is V, U. 4. The strain energy density theory Carbon fiber epoxy composites T300/5208 is used to survey the fracture angle of the propagating crack. The influence of crack propagation velocity on fracture angle is studied when the fiber orientation angle a = 0 and M = 0–0.7. 4.1. The strain energy density ratio criterion The strain energy density (SED) criterion can be applied to predict the crack propagation in isotropic materials [21–23]. The strain energy density can be written as

dW 1 ¼ rij eij : dV 2

ð35Þ

Zhang and Cai [24] proposed a strain energy density ratio (SEDR) criterion for predicting cracking direction in composites. For orthotropic materials and plane problems, the SEDR can be written as

dW 1 ¼ rij eij : dV 2

ð36Þ

Define SR as the strain energy density factor such that

dW SR ¼ : dV r

ð37Þ

The stress ratio rij and the strain ratio eij were, respectively, defined as

rx ¼ rx =X; ry ¼ ry =Y; sxy ¼ sxy =S; ex ¼ ex =exc ; ey ¼ ey =eyc ; cxy ¼ cxy =cxyc ; exc ¼ X=E1 ; eyc ¼ Y=E2 ; cxyc ¼ S=G12 :

ð38Þ ð39Þ ð40Þ

Substituting Eqs. (38)–(40) into Eq (37), Eq. (37) takes the specific form

SR ¼

Fig. 6. Mode I V=K I versus angle h.



 2  r rx 2 ry 2 m21 E E s þ  rx ry 12 þ 22 þ xy ; 2 X Y E1 S X Y

ð41Þ

where X and Y are the tensile strength in the principal direction, and S is shear yield strength. ex, ey are the normal strain in the principal direction, and cxy is shear strain, and exc, eyc, cxyc are the critical values of them, respectively. E1 and E2 are Young modulus, G12 is shear modulus and v21 is in-plane Poisson’s ratio. The subscript 1 and 2 represent the principal axes of material in the longitudinal and transverse directions, respectively.

X. Gao et al. / Theoretical and Applied Fracture Mechanics 51 (2009) 73–85

81

Two fundamental hypotheses of crack extension in compose materials can be stated that Hypothesis 1. Crack propagation will start in a radial direction along which the inner energy density is a minimum. Hypothesis 2. The crack is assumed to take place when Smin reaches a critical value Sc. Hence, the necessary and sufficient conditions for SR to be a minimum value Smin are

@SR ¼ 0; @h

@ 2 SR @h2

> 0:

ð42Þ

From Eq. (42), the direction of crack propagation h0 can be determined. The driving force Smin can then be obtained by substituting h0 into Eq. (37). 4.2. The influence of crack propagation velocity Fig. 12 shows the variations of h0 with crack propagation velocity M for various values of m. The magnitude of m dominates the fracture behavior. The greater the value of m is, the greater the fracture angle h0 is. For a specific mode crack, the fracture angle h0 decreases gradually with the increment of the crack propagation velocity M. When m is as small as 0.6, the crack behaves similarly to the case of m = 0. The influence of M on the fracture angle h0 can be ignored and the fracture angle h0 almost holds the line. Note that when m = 0 the crack propagation velocity effect disappears. The fracture angle h0 = 0° for I mode crack, and the angle is about 65° < h0 < 69° for the different M in case of II mode. It is known from above results that the crack extension depends on the crack propagation velocity. 4.3. The influence of anisotropic constants In order to investigate the influence of anisotropy constants, one define b1 = E2/G12, b2 = E1/E2, and analyze the influence of b1 and b2 on the fracture angle when M = 0. For fixed value of E1, E2, m21, X, Y, S, and b1 = 0.5, 0.75, 1, 1.25, 1.5. Fig. 13 shows the influence of b1 on the fracture angle h0 for various values of m. The greater the value of m is, the greater the fracture angle h0 is. When 0 6 m 6 1, the influence of anisotropic constant b1 disappears and the fracture angle h0 holds the line. The fracture angle h0 increases gradually as b1 increases when 1 < m 6 5. However, the fracture angle h0 decreases gradually as b1 increases when m = 20 and m = 50. It indicates that the fracture angle decreases gradually with the increment of b1 for II mode crack. For fixed value of E2, G12, m21, X, Y, S, and b2 = 15, 17.5, 20, 22.5, 25. Fig. 14 shows the influence of b2 on the fracture angle h0 for various values of m. The magnitude of m dominates the fracture behavior. The greater the value of m is, the greater the fracture angle h0 is. When m = 0, i.e. for I mode crack, the influence of anisotropic constant b2 disappears and the fracture angle h0 holds the line. The fracture angle h0 decreases gradually as b2 increases for 0 < m 6 5. On the contrary, the fracture angle h0 increases gradually as b2 increases when m = 20 and m = 50. It indicates that the fracture angle increases gradually with the increment of b2 for II mode crack. It is known from Figs. 13 and 14 that the crack extension also depends on elastic properties of the composite materials. 4.4. The SED factor Fig. 7. Mode II

rx =K I versus angle h.

While there are many minima of the strain energy density function and hence SED factor, it is the maximum of the Smin that is assumed to initiate crack propagation. As it is to be expected, the results are sensitive to the mix mode coefficient m = KII/KI. The like-

82

X. Gao et al. / Theoretical and Applied Fracture Mechanics 51 (2009) 73–85

Fig. 8. Mode II

ry =K II versus angle h.

Fig. 9. Mode II

sxy =K II versus angle h.

X. Gao et al. / Theoretical and Applied Fracture Mechanics 51 (2009) 73–85

Fig. 10. Mode II U=K II versus angle h.

lihood dynamic instability such as bifurcation is increased with increasing m and M. The effects of the degree of the material

83

Fig. 11. Mode II V=K II versus angle h.

anisotropy as reflected by the parameters b1 and b2 are shown respectively in Figs. 16 and 17.

84

X. Gao et al. / Theoretical and Applied Fracture Mechanics 51 (2009) 73–85

Fig. 12. Fracture angle h0 versus crack velocity M.

Fig. 15. Normalized Smin versus M with varying m.

Fig. 13. Variations of h0 with moduli ratio b1.

Fig. 16. Normalized Smin versus b1 with varying m(M = 0).

Fig. 14. Variations of h0 with moduli ratio b2.

Fig. 17. Normalized Smin versus b2 with varying m(M = 0).

X. Gao et al. / Theoretical and Applied Fracture Mechanics 51 (2009) 73–85

Fig. 15 shows the variations of minimum SED factor Smin with crack propagation velocity M for various values of m. The magnitude of m dominates the value of Smin. The greater the value of m is, the greater the minimum SED factor Smin is. For a specific mode crack, the minimum SED factor Smin decreases gradually with the increment of the crack propagation velocity M when m > 1. The influence of M on the minimum SED factor Smin disappears and the minimum SED factor Smin holds the line when m = 1. The minimum SED factor Smin increases slightly with the increment of the crack propagation velocity M when m < 1. Fig. 16 shows the influence of b1 on the minimum SED factor Smin for various values of m. The greater the value of m is, the greater the minimum SED factor Smin is. When 0 6 m < 1, the influence of anisotropic constant b1 disappears and the minimum SED factor Smin holds the line. The minimum SED factor Smin increases gradually as b1 increases when 1 6 m. Fig. 17 shows the influence of b2 on the minimum SED factor Smin for various values of m. The magnitude of m dominates the value of Smin. The greater the value of m is, the greater the minimum SED factor Smin is. The minimum SED factor Smin increases slightly as b2 increases for m 6 2.5. On the contrary, the minimum SED factor Smin decreases gradually as b2 increases when m > 3. Note that the influence of b2 on the minimum SED factor Smin disappears and the minimum SED factor Smin holds the line when m = 3. It is known from above results that the minimum SED factor Smin depend not only on the crack propagation velocity M but also on elastic properties of the composite materials. Compared with the anisotropic constants, the crack propagation velocity M gives greater influence on the minimum SED factor Smin. When m < 1, the influence of M, b1 and b2 on the minimum SED factor Smin is slight. On the contrary, the influence is great when m P 10. It indicates that the crack propagation velocity and anisotropic constants effect on Smin almost disappears for I mode crack (m = 0), and the influence is greater for II mode crack (m P 10) compared with I mode crack.

5. Results The following results are obtained on the dynamic fracture analysis of anisotropic cracked plates: (1) The expression of dynamic stress intensity factors for mode I and II crack are obtained, and the dynamic stress components and dynamic displacement components around the crack tip are derived. The obtained crack tip fields are unified and applicable to the analysis of the stress fields, strain fields and displacement fields of mode I, mode II and I/II mixed mode for anisotropic material, orthotropic material and isotropic material under dynamic or static load. (2) The velocity of crack propagation M obviously affects the stress fields and displacement fields. The faster the crack velocity is, the greater the dynamic stress components and dynamic displacement components around crack tip is. (3) The fiber direction a gives greater influence on the variations of the stress fields and displacement fields. The variations of the stress fields and displacement fields with h are greater when a = p/6, p/3 than when a = 0, p/2. (4) Distributions of the stress fields and displacement fields vary with a but the inclinations of stress distribution are almost the same. When a = 0, p/2, the crack tip stress fields and displacement fields have a symmetrical or dissymmetrical distribution to the crack plane at any velocity of crack propagation.

85

(5) The fracture angle is found to depend not only on the velocity of crack propagation M but also on the anisotropic character of the material. The faster the crack velocity is, the less the dynamic crack extension angle is. The anisotropic material constants b1 and b2 give a little influence on the fracture angle. (6) The crack extension angle of mode I crack is 0 for the different crack propagation velocity and anisotropic material parameters. (7) The minimum SED factor Smin depends not only on the crack propagation velocity M but also on elastic properties of materials. Compared with the anisotropic constants, the crack propagation velocity M gives greater influence on Smin. When m < 1, the influence of M, b1 and b2 on the minimum SED factor Smin is slight. The influence is great when m P 10. Acknowledgement This work has been supported by National Advance Research Foundation of China (No. 513270301). References [1] N.F. Mott, Fracture of metals: theoretical consideration, Engineering 165 (1948) 16–18. [2] E.H. Yoffe, The moving Griffith crack, Philos. Mag. 42 (1951) 739–750. [3] E.P. Chen, Sudden appearance of a crack in a stretched finite strip, J. Appl. Mech. 45 (1978) 270–280. [4] B.R. Baker, Dynamic stresses created by a moving crack, J. Appl. Mech. 29 (1962) 449–458. [5] L.B. Freund, Dynamic crack propagation, in: F. Erdogan (Ed.), Mechanics of Fracture, vol. 19, ASME, 1976, pp. 105–134. [6] M.K. Kassir, S. Tse, Moving Griffith crack in an orthotropic material, Int. J. Eng. Sci. 21 (1983) 315–325. [7] M. Arcisz, G.C. Sih, Effect of orthotropy on crack propagation, Theor. Appl. Fract. Mech. 1 (1984) 225–238. [8] J.D. Achenbach, Z.P. Bazant, Elastodynamic near-tip stress and displacement fields for rapidly propagating crack in orthotropic materials, J. Appl. Mech. 42 (1975) 183–189. [9] A. Piva, E. Viola, Crack propagation in an orthotropic medium, Eng. Fract. Mech. 29 (1988) 535–548. [10] E. Viola, A. Piva, E. Radi, Crack propagation in an orthotropic medium under general loading, Eng. Fract. Mech. 34 (1990) 1155–1174. [11] K.H. Lee, J.S. Hawong, S.H. Choi, Dynamic stress intensity factors KI, KII and dynamic crack propagation characteristics of orthotropic material, Eng. Fract. Mech. 53 (1996) 119–140. [12] H.M. Xu, X.F. Yao, X.Q. Feng, Dynamic stress intensity factor KIII and dynamic crack propagation characteristics of orthotropic material, Eng. Mech. 23 (2006) 68–72. [13] C. Rubio-Gonzalez, J.J. Mason, Dynamic stress intensity factors at the tip of a uniformly loaded semi-infinite crack in an orthotropic material, Int. J. Mech. Phys. Solids 48 (2000) 899–925. [14] X. Gao, H.G. Wang, X.W. Kang, Dynamic stress intensity factor KIII and dynamic crack propagation characteristics of anisotropic materials, Appl. Math. Mech. 29 (2008) 1119–1129. [15] G.C. Sih, Some basic problems in fracture mechanics and new concepts, J. Eng. Fract. Mech. 5 (1973) 365–377. [16] G.C. Sih (Ed.), Mechanics of Fracture, vols. I–VII, Martinus Nijhoff, Netherlands, 1973. [17] E.E. Gdoutos, Problems of mixed mode crack propagation, Engineering Application of Fracture Mechanics, vol. II, Martinus Nijhoff, Netherlands, 1984. [18] A. Carpinteri, Crack growth and material damage in concrete: plastic collapse and brittle fracture, Engineering Application of Fracture Mechanics, vol. V, Martinus Nijhoff, Netherlands, 1986. [19] S. Shen, T. Nishioka, Fracture of piezoelectric materials: energy density criterion, Theor. Appl. Fract. Mech. 33 (2000) 57–65. [20] G.C. Sih, P.C. Paris, G.R. Irwin, On crack in rectilinearly anisotropic bodies, Int. J. Fract. Mech. 1 (1965) 189–203. [21] G.C. Sih, A special theory of crack propagation: methods of analysis and solutions of crack problems, Mechanics of Fracture I, Noordhoff, Leyden, 1973. pp. 21-45. [22] G.C. Sih, Strain energy factor applied to mixed mode crack problems, Int. J. Fract. 10 (1974) 305–321. [23] G.C. Sih, Mechanics of Fracture Initiation and Propagation, Kluwer Academic Publisher, 1991. [24] S.Y. Zhang, L.W. Cai, The strain energy density ratio criterion for predicting cracking direction in composite materials, Appl. Math. Mech. 14 (1993) 79– 87.