DYNAMIC CRACK PROPAGATION BETWEEN TWO BONDED ORTHOTROPIC PLATES

(1)

DYNAMIC CRACK PROPAGATION BETWEEN TWO BONDED ORTHOTROPIC PLATES

M. S. MATBULY

Received 30 June 2003 and in revised form 30 August 2003

The problem of crack propagation along the interface of two bonded dissimilar orthotropic plates is considered. Using Galilean transformation, the problem is reduced to a quasistatic one. Then, using Fourier transforms and asymptotic analysis, the problem is reduced to a pair of singular integral equations with Cauchy-type singularity.

These equations are solved using Gauss-Chebyshev quadrature formulae. The dynamic stress intensity factors are obtained in closed form expressions. Furthermore, a paramet- ric study is introduced to investigate the eﬀect of crack growth rate and geometric and elastic characteristics of the plates on values of dynamic stress intensity factors.

1. Introduction

Composite materials have been extensively employed in many engineering fields such as mechanical and aerospace structures. When the material used as member of such structures contains a crack, it is seriously necessary to study the stress field distribution at the immediate vicinity of crack tips. The inertia action must be considered when the applied loads or crack length depend on time. Also, the most frequently observed phenomenon in the experiments shows that the crack growth rate is constant during the extending his- tory except in the final unstable or arresting stage [11]. So, the elastodynamic analysis of a moving crack, with constant velocity, is one of the most important problems in fracture mechanics. The dynamic stress intensity factors (DSIF) play a key role in characteriz- ing the fracture behaviour of such problems. Thus, analytical determination of DSIF in predicting the fracture cannot be overemphasized.

In general, there are two approaches for analytical determination of DSIF. The first one employs the integral transforms and asymptotic analysis to reduce the problem to that of a system of singular integral equations [1,3,10,15,20,22,23]. The second approach employs complex analysis to reduce the problem to that of a system of Riemann-Hilbert problems [12,13,14,19].

The present work is concerned with elastodynamic stress disturbance problem of a moving Griﬃth crack with constant velocity. The crack is located at the interface of two

(2)

σ₁₂¹ =σ₂₂¹ =0

σ₁₂² =σ₂₂² =0 a0+ct a0+ct

X2

σ₂₂¹ =F1(X1)

σ₁₂¹ =F2(X1)

σ₁₂² =F2(X1)

σ₂₂² =F1(X1)

X1

=1

=2 H1

H2

Interfacial line

Plate 1

Plate 2

Figure 1.1. Two bonded dissimilar plates containing moving interfacial crack.

bonded dissimilar orthotropic plates, as shown inFigure 1.1. Each plate possesses a finite width and is subjected to a static stress distribution along crack surfaces. This is the main diﬀerence between the present work and previous ones [1,3,4,6,7,8,10,12,13,14,15, 19,20,22,23], which were concerned only with plates of infinite widths. The governing equations of the problem are described. Then, Fourier transforms and asymptotic analysis are employed to reduce the solution of the problem to that of a system of first-kind singular integral equations with Cauchy-type singularity. These are solved numerically according to the algorithm in [17]. Then, closed-form expressions for the asymptotic stress field distribution at the immediate vicinity of crack tips are obtained.

2. Governing equations

Assuming that the Cartesian coordinate axes are the axes of symmetry of the elastic materials, the displacement equations of motion for orthotropic plates are [3]

C₁₁∂²U₁

∂X₁² +C₆₆ ∂²U₁

∂X₂² +C₁₂+C₆₆ ∂²U₂

∂X1∂X2 =m∂²U₁

∂t² , C₆₆∂²U₂

∂X₁² +C₂₂ ∂²U₂

∂X₂² +C₁₂+C₆₆ ∂²U₁

∂X1∂X2 =m∂²U₂

∂t² ,

(2.1)

where is a superscript (=1 for orthotropic material in X2>0, while =2 for orthotropic material inX2<0), as shown inFigure 1.1;U_j (j=1, 2) are the displacement components in direction ofX1andX2, respectively;C_{i j}(i,j=1, 2) andC₆₆ are the elastic constants of orthotropic plate materials; andm andtare the material mass density and time, respectively.

The boundary conditions along the interface of plates are

Xlim2→0⁺σ₂₂¹X1,X2,t= lim

X2→0⁻σ₂₂²X1,X2,t=F1

X1

,

Xlim2→0⁺σ₁₂¹X1,X2,t= lim

X2→0⁻σ₁₂²X1,X2,t=F2

X1

, X1< a0+d,

(3)

Xlim2→0⁺U₁¹X1,X2,t= lim

X2→0⁻U₁²X1,X2,t,

Xlim2→0⁺U₂¹X1,X2,t= lim

X2→0⁻U₂²X1,X2,t,

X1> a0+d

Xlim2→0⁺σ₂₂¹X1,X2,t= lim

X2→0⁻σ₂₂²X1,X2,t,

Xlim2→0⁺σ₁₂¹X1,X2,t= lim

X2→0⁻σ₁₂²X1,X2,t,

X1> a0+d,

(2.2) wherea0is the initial half crack length, andd= |c|t, wherecis the magnitude of crack propagation velocity. Moreover,F_i(X1) (i=1, 2) are known functions. They represent the applied static stress along crack surfaces;F1(−X1)=F1(X1) andF2(−X1)= −f2(−X1).

The boundary conditions along the external boundaries are σ12¹

X1,H1,t=σ22¹

X1,H1,t=σ12²

X1,−H2,t

=σ₂₂²X1,−H2,t=0, X1<∞, (2.3) whereH1andH2are as shown inFigure 1.1.

Using Galilean transformation:x=(X1−ct)/a0, y=X2/a0, andt=t, the governing equations (2.1), (2.2), and (2.3) can be reduced to a quasistatic form as follows:

∂²u

∂x² +a₁∂²ν

∂x∂y+a₂∂²u

∂y² ⁼0,

∂²ν

∂x² +b₁∂²u

∂x∂y+b₂∂²ν

∂y² ⁼0,

(=1, 2) (2.4)

where

u(x,y)=U1

X1,X2,t

a0 , ν(x,y)=U2

X1,X2,t

a0 , (2.5a)

a₁= C12+C66

C₁₁ 1−

M₁², a₂= C66

C₁₁1−

M₁², (2.5b)

b₁= C₁₂+C₆₆ C66

1− M2

2, b₂= C₂₂ C66

1− M2

2, (2.5c)

M_j = c

V_j (j=1, 2), V₁=

C₁₁

m, V₂=

C₆₆

m. (2.5d)

Furthermore, the values of Mach numbersM_j (,j=1, 2) are assumed to be less than unity for subsonic crack propagation.

(4)

The boundary conditions for the reduced quasistatic problem can be expressed as follows:

ylim→0⁺σ_{y y}¹ (x,y)= lim

y→0⁻σ²_{y y}(x,y)= f1(x),

ylim→0⁺σ_xy¹(x,y)= lim

y→0⁻σ_xy²(x,y)=f2(x), ^|x|<1, (2.6)

ylim→0⁺u¹(x,y)₌ lim

y→0⁻u²(x,y),

ylim→0⁺ν¹(x,y)= lim

y→0⁻ν²(x,y), ^|x|>1, (2.7)

ylim→0⁺σ_{y y}¹ (x,y)= lim

y→0⁻σ²_{y y}(x,y),

ylim→0⁺σ_xy¹(x,y)= lim

y→0⁻σ_xy²(x,y), ^|x|>1, (2.8) σ¹_{y y}x,h1

=σ_xy¹x,h1

=σ_{y y}² x,−h2

=σ_xy² x,−h2

=0, |x|<∞, (2.9)

where

σ_xy =σ₁₂

Q, σ_{y y}=σ₂₂

Q , Q= C¹66C²66, h=H

a0 (=1, 2).

(2.10)

3. Solution of the problem

Decoupling (2.4) then employing Fourier sine (cosine) transforms with respect tox, one can find that

u(x,y)= 1 π

4 j=1

∞

0 A_j(α)e^αr^j^ysinαx dα

, ν(x,y)= 1

π 4 j=1

∞

0 k_jA_j(α)e^αr^j^ycosαx dα

,

(=1, 2) (3.1)

whereαis the transform variable,A_j(α) (=1, 2 andj=1, 4) are unknown functions, andr_j(=1, 2 andj=1, 4) are the real distinct roots of

f₁r⁴−2f₂r²+ 1=0 (=1, 2), (3.2) provided that

f₁=a₂b₂, 2f₂=a₂+b₂−a₁b₁, f₂>

f₁, f₁>0, (3.3) k_j=a₂r_j²−1

a₁r_j (=1, 2, j=1, 4). (3.4)

(5)

From the stress-displacement relationship of orthotropic materials [9], one can get σ_xy (x,y)=C66

Q ∂ν

∂x +∂u

∂y

=1 π

4 j=1

∞

0 αp_jA_j(α)e^αr^j^ysinαx dα

(=1, 2), σ_{y y} (x,y)= 1

Q

C₁₂ ∂u

∂x +C₂₂∂ν

∂y

=1 π

4 j=1

∞

0 αo_jA_j(α)e^αr^j^ycosαx dα

(=1, 2),

(3.5)

where

p_j=C₆₆ r_j−k_j

Q , o_j=

C₁₂ +C₂₂r_jk_j

Q . (3.6)

On suitable substitution from (2.9) into (3.5), one can find that A_m(α)=

4 j=3

L_jme^αβ^jmA_j(α) (m,=1, 2), (3.7) where

β_jm=



 h1

r_j−r_m for=1 h2

r_m −r_j for=2 (m=1, 2, j=3, 4), (3.8)

L_jm= 2

n=1(−1)ⁿ1−δ_nmp_no_j−p_jo_n

p₁o₂−p₂o₁ (,m=1, 2, j=3, 4), (3.9) δnmis the Kronecker delta.

Substituting (3.5) and (3.7) into (2.8), one can find that A²₃(α)=E₄²D₃¹−D²₄E¹₃

∆ A¹₃(α) +E²₄D¹₄−D²₄E₄¹

∆ A¹₄(α), A²₄(α)=D²3E¹3−E²3D¹3

∆ A¹₃(α) +D3²E¹4−E3²D4¹

∆ A¹₄(α),

(3.10)

where

∆=D²₃E²₄−D₄²E₃², D_j=p_j+

2 m=1

p_mL_jme^αβ^jm, E_j=o_j+ 2 m=1

o_mL_jme^αβ^jm (=1, 2, j=3, 4). (3.11)

(6)

The boundary conditions (2.6) and (2.7) in conjunction with (3.7), (3.8), (3.9), (3.10), and (3.11) lead to

1 π

_∞

0 αE¹₃A¹₃(α) +E¹₄A¹₄(α)cosαx dα=f1(x), 1

π _∞

0 αD¹₃A¹₃(α) +D¹₄A¹₄(α)sinαx dα=f2(x),

|x|<1, (3.12) 1

π _∞

0 ατ1(α)A¹₃(α) +τ2(α)A¹₄(α)sinαx dα=0, 1

π _∞

0 αξ1(α)A¹₃(α) +ξ2(α)A¹₄(α)cosαx dα=0,

|x_|>1, (3.13)

where τ1(α)=

1 +

2 m=1

L¹_3me^αβ¹^3m−

1 + 2 m=1

L²_3me^αβ²^3m

∆1−

1 + 2 m=1

L²_4me^αβ²^4m

∆2

, τ2(α)=

1 +

2 m=1

L¹_4me^αβ¹^4m−

1 + 2 m=1

L²_3me^αβ²^3m

∆3−

1 + 2 m=1

L²_4me^αβ²^4m

∆4

, ξ1(α)=

k¹₃+

2 m=1

k_m¹L¹_3me^αβ¹^3m−

k²₃+ 2 m=1

k²_mL²_3me^αβ²^3m

∆1

−

k²₄+ 2 m=1

k_m²L²_4me^αβ²^4m

∆2

, ξ2(α)=

k¹₄+

2 m=1

k_m¹L¹_4me^αβ¹^4m−

k²₃+ 2 m=1

k²_mL²_3me^αβ²^3m

∆3

−

k²₄+ 2 m=1

k_m²L²_4me^αβ²^4m

∆4

,

∆1=E²₄D¹₃−D₄²E¹₃

∆ , ∆2=D₃²E¹₃−E₃²D₃¹

∆ ,

∆3=E²₄D¹₄−D₄²E¹₄

∆ , ∆4=D₃²E¹₄−E₃²D₄¹

∆ .

(3.14)

The unknown functionsA¹_j(α) (j=3, 4) can be determined through solving the integral equations (3.12) and (3.13) as follows [21].

Let

ylim→0⁺

∂ν¹(x,y)

∂x ⁻ lim

y→0⁻

∂ν²(x,y)

∂x ⁼φ1(x)1−H(x−1),

ylim→0⁺

∂u¹(x,y)

∂x ⁻lim

y→0⁻

∂u²(x,y)

∂x ⁼φ2(x)1−H(x−1),

(3.15)

whereH(x−1) is the unit step function [2], whileφ_j(x) (j=1, 2) are unknown odd and even functions ofx, respectively.

(7)

From (3.13) and (3.15), one can deduce that A¹₃(α)=τ2Φ1(α) +ξ2Φ2(α)

(αΨ) , A¹₄(α)=−

τ1Φ1(α) +ξ1Φ2(α)

(αΨ) ,

(3.16)

where

Φ1(α)=2 1

0φ1(x) sinαx dx, Φ2(α)=2

1

0φ2(x) cosαx dx, Ψ=τ1ξ2−τ2ξ1.

(3.17)

Substituting (3.16) into (3.12), the problem is reduced to the following system of integral equations:

2 j=1

₁

−1

H˘i j(x,t)φj(t)dt= fi(x), |x|<1,i=1, 2, (3.18)

where the kernels ˘Hi j(x,t) are H˘11(x,t)= 1

π _∞

0

1 Ψ

E₃¹τ2−E₄¹τ1

sinα(t−x)dα, H˘12(x,t)= 1

π _∞

0

1 Ψ

E₃¹ξ2−E¹₄ξ1

cosα(t−x)dα, H˘21(x,t)= 1

π _∞

0

1 Ψ

D¹₃τ2−D¹₄τ1

cosα(t−x)dα, H˘22(x,t)= −1

π _∞

0

1 Ψ

D¹₃ξ2−D₄¹ξ1

sinα(t−x)dα.

(3.19)

Since the integrands of (3.19) are continuous functions ofα, then it is clear that any possible singularity of the kernels must result from the asymptotic analysis of the integrands ast→xand α→ ∞[16,21]. Then, by adding and subtracting the asymptotic expressions of these integrands under the integral sign, the problem can be reduced to the following pair of singular integral equations with Cauchy-type singularity:

1 π

2 j=1

δ_{i j}G_{i j} 1

−1

φ_j(t) t−xdt+

1

−1H_{i j}(x,t)φ_j(t)dt

=f_i(x), |x|<1,i=1, 2, (3.20)

(8)

where

H11(x,t)= _∞

0

1 Ψ

E₃¹τ2−E₄¹τ1

−G11

sinα(t−x)dα, H12(x,t)=

_∞

0

1 Ψ

E₃¹ξ2−E¹₄ξ1

−G12

cosα(t−x)dα, H21(x,t)=

_∞

0

1 Ψ

D¹₃τ2−D¹₄τ1

−G21

cosα(t−x)dα, H22(x,t)=

_∞

0

−1 Ψ

D¹₃ξ2−D¹₄ξ1

−G22

sinα(t−x)dα, G11=

o¹3

B−s2

/B−o¹4

B−s1

/B

Z ,

G12=

o¹₃k¹₄B−s4

/B−o¹₄k¹₃B−s3 /B

Z ,

G21=

p₃¹B−s2

/B−p¹₄B−s1

/B

Z ,

G22=

p₄¹k₃¹B−s3

/B−p¹₃k₄¹B−s4

/B

Z ,

Z=

k¹₄−k¹₃+s3−s4+k3¹s2−k4¹s1

B +s1s4−s2s3

B² , B=

L²₃₁L²₄₂−L²₄₁L²₃₂p²₁o²₂−p²₂o²₁, s1=

B13+B24

L²₃₁L²₄₂+B12+B27 L²₄₁L²₃₂, s2=

B33+B44

L²₃₁L²₄₂+B32+B47

L²₄₁L²₃₂, s3=

k²₁B13+k₂²B24

L²₃₁L²₄₂+k₂²B12+k²₁B27

L²₄₁L²₃₂, s4=

k²1B33+k2²B44

L²31L²42+k2²B32+k²1B47

L²41L²32, B13=p₃¹o²₂−p²₂o¹₃, B24=p²₁o¹₃−p¹₃o²₁,

B12=p₃¹o²₁−p²₁o¹₃, B27=p²₂o¹₃−p¹₃o²₂, B33=p₄¹o²₂−p²₂o¹₄, B44=p²₁o¹₄−p¹₄o²₁, B32=p₄¹o²₁−p²₁o¹₄, B47=p²₂o¹₄−p¹₄o²₂.

(3.21)

From (3.15) and (2.5d), one can deduce the single-valuedness conditions ensuring the uniqueness ofφ_i(α) (i=1, 2) as follows:

₁

−1φ1(t)dt=0, ₁

−1φ2(t)dt=2 c

C²₁₁ m² ⁻

C¹₁₁

m¹

. (3.22)

Equations (3.20) and (3.22) can be solved as follows [17].

(9)

Assume that

φ_i(t)= qi(t)

√1−t² (i=1, 2), (3.23)

whereqi(t) (i=1, 2) are bounded continuous functions for allt∈[−1, 1].

By substituting from (3.23) into (3.20) and (3.22) then employing Gauss-Chebyshev integration formulae, the solution of the problem can be reduced to the following system of linear algebraic equations [17]:

2 j=1

n1

k=1

wk

δi jGi j

tk−xl+Hi j xl,tk

qj tk

=fi xl

(i=1, 2)xl=1,n1−1,

n1

k=1

wkq1

tk

=0,

n1

k=1

w_kq2

t_k= 2 πc

C²11

m² ⁻

C¹11

m¹

,

(3.24)

wheren1is the number of collocation points in the interval [−1, 1],

w1=wn1= 1

2n1−2, wk= 1 n1−1

fork=2,n1−1,

tk=cos π(k−1)

n1−1 k=1,n1

,

x_l=cos π(2l−1)

2n1−2 l=1,n1−1.

(3.25)

For the concerned problem, one can deduce that the DSIF at the left and right crack tips are equal. Then, by making use of the following asymptotic relations, asα→ ∞, [5]:

Φ1(α)= ₁

−1

q1(t)

√1−t²sinαt dt≈q1(1) 2π

α sin α−π 4

+ϑ 1 α

, Φ2(α)=

₁

−1

q2(t)

√1−t²cosαt dt≈q2(1) 2π

α cos α−π 4

+ϑ 1 α

, _∞

0

√1 αe⁻^bα

sinhα coshα

dα=

√π b²+h²^0.25

sin

cos 0.5 tan⁻¹ h b

, b >0,

(3.26)

(10)

the leading terms of the asymptotic stress field distribution (3.5) at the immediate vicinity of the right crack tip (x→1 andy→0^±) can be obtained as follows:

σxy(x,y)≈q1(1) ζ1

2ρ1

sinθ1

2 ⁻ ζ2

2ρ2

sinθ2

2

+q2(1) ζ3

2ρ1

cosθ1

2 ⁻ ζ4

2ρ2

cosθ2

2

, σy y(x,y)≈q1(1) γ1

2ρ1cosθ1

2 ⁻ γ2

2ρ2cosθ2

2

+q2(1) γ3

2ρ1sinθ1

2 ⁻ γ4

2ρ2sinθ2

2

,

(3.27)

where

ρ1=

(x−1)²+r3¹y², ρ2=

(x−1)²+r4¹y², θ1=tan⁻¹ r3¹y

x−1

, θ2=tan⁻¹ r4¹y x−1

, ξ1= p¹₃1−s2/B

Z , ξ2= p¹₄1−s1/B

Z ,

ξ3= p¹₃k¹₄−s4/B

Z , ξ4= p¹₄k¹₃−s3/B

Z ,

γ1=o¹₃1−s2/B

Z , γ2=o¹₄1−s1/B

Z ,

γ3=o¹3

k4¹−s4/B

Z , γ4=o¹4

k¹3−s3/B

Z .

(3.28)

Therefore, the DSIF,K_I, andK_IIcan be determined as follows [18]:

K_I+iK_II=Q2a0 lim

x→1,y→0⁺

√x−1σ_{y y}(x,y) +iσ_xy(x,y) (3.29)

such that by substituting from (3.27) and (3.28) into (3.29) then evaluating the limits, one can find that

KI=Q^√a0 γ1−γ2

q1(1), KII=Q^√a0

ξ3−ξ4

q2(1). (3.30)