Steady motion of a crack parallel to a bond-plane

hf. 1. Engng SC;. Vol. 18, pp. 225-138 Pergamon Press Ltd., 1980 Printed in Great Britain

STEADY MOTION OF A CRACK PARALLEL TO A BOND-PLANE S. H. CHEN, L. M. KEER and J. D. ACHENBACH Department of Civil Engineering, Northwestern University, Evanston, IL 60201,U.S.A. combined effects of high crack-tip speedand the proximityof a bond-plane on the elastodynamic stress-intensity factorare investigated in thispaper.The model-problem that is considered

Abstract-The

concerns the steady propagation of a crack of length 2a, parallel to a bond-plane with a half-plane of different material properties. By using a moving coordinate system and applying Fourier transform techniques and superposition methods, the mixed boundary-value problem is reduced to a dual singular integral equation with Cauchy-type kernels, which is solved numerically. Stress intensity factors and cleavage angles are computed with the material constants, crack-interface distance and the crack-tip speed as variables. It is found that the effect of the bond-plane is more intense for smaller crack-interface distance and larger crack-tip speed. The propagation speed is subsonic for both materials.

INTRODUCTION

A BONDPLANE betwen two materials of different mechanical properties often provides a natural path for crack propagation. It is not uncommon that a crack originates and continues to propagate in such an interface. The rather pathological behavior predicted by linearized elasticity for a crack in a plane of otherwise perfect bond is well known, see e.g. EnglandLll. In[l] it was pointed out that oscillating singularities in the stress fields near the crack tips are present and that a peculiar overlap of the materials in this region will also occur, in contradiction to the physical boundary conditions. Achenbach and Bazant[2] and AtkinsonL31 showed that a propagating interface crack will have similar oscillating singularities in the elastodynamic near-tip stress field. In recent work two different approaches to the question of overlapping of materials near the crack tip have been proposed [4,5] to obtain a physically reasonable solution within the range of linear elasticity. Alternatively, a crack may have its origin at a position near the interface, follow a path close to the interface, and then either be attracted or repelled. For static problems cracks near and at interfaces have been studied by Erdogan[6, 73 and Erdogan and Gupta[& 91. For dynamic problems the mode of crack propagation parallel to an interface has apparently not been studied in any detail. Dynamic effects generated by rapid crack propagation tend to affect the angular variatidn of the near-tip stress fields, and thus possibly affect the directional stability of a rapidly propagating crack [ lo]. In the present paper the additional effect of an interface on these near-tip stress fields is studied for a crack propagating rapidly parallel to the interface. An interesting model to study the influence of a bond-plane would be a crack of timeincreasing width, with the crack tips propagating in different direction parallel to the bondplane. Unfortunately this problem is very diacult. However, after such a crack has been propagating for a short time at a constant velocity a relatively steady state situation near the crack tips may be assumed, with the instantaneous crack length entering as a parameter. This suggests a problem which is amenable to mathematical analysis, namely, the steady-state problem where both crack tips propagate in the same direction. On fist thought this model would seem rather unrealistic since the trailing crack tip closes, but then if the interest is focussed on the leading crack tip, the finite length of the crack introduces a relevant instantaneous length dimension in the analysis, and the fact that the trailing crack tip moves in the wrong direction may be assumed of minor importance for the field near the leading crack tip. The problem that is to be solved then, is for internal pressure applied to the faces of a crack of finite length, 2a moving with a steady velocity c. Yoffe[lO] considered a similar model for a crack moving with a steady velocity in a homogeneous elastic material. It was found in [ 10, 1l] that the maximum of the cleavage stress shifts out of the plane of the crack when the propagation velocity exceeds a certain value. Since the stress intensity factors for a crack near an 22s

S. H. CHEN et al.

226

interface will not only depend on the speed but also on the relative depth below the interface and the material constants, one can anticipate that the shift of the maximum cleavage stress will also depend upon these geometrical and material parameters. To develop a solution a coordinate system attached to the moving crack is used and the elastodynamic equations are written for the case of steady motion. By the use of Fourier integral transforms the boundary value problem is reduced to two simultaneous singular integral equations of the first kind with Cauchy kernels. The numerical scheme of Erdogan and Gupta[ 121is used to solve these equations. Material constants, geometrical parameters, and the crack velocity are the parameters that affect the stress and displacement fields. The tendency of the crack to move out of its plane, either towards or away from the interface, can be inferred from the cleavage angle which is obtained from the stress intensity factors and the near-tip stress fields. For the considered analysis the propagation speed is subsonic for both materials. GOVERNING

EQUATIONS

The geometrical configuration and the coordinate system for the problem are shown in Fig. la. A crack of constant length 2a, moving at a uniform speed c to the right, is located at a fixed distance h below the bond-plane of two elastic half planes. The symbols IL, v, and p denote the shear modulus, Poisson’s ratio, and the mass density, respectively, and the subscripts 1 and 2 refer to materials 1 and 2. A moving coordinate system (x, y), which is attached to the moving crack [ 133is introduced as x = x’ -

ct,

(1) (2)

Y = Y',

where (x’, y’) is a fixed coordinate system. To solve the governing elastodynamic equations of plane elasticity, two displacement potentials cp and $ are introduced such that

ux= Q,x +

+,yr

(3) (4)

u, = (P,y- &I, where a comma denotes partial differentiation. The expressions for the stresses in terms of cp and 4 follow from Hooke’s law as ax = AV2Q + 2/4Qm+ &,xy),

(5)

“y = AV2Q + %(Q,yy

(6)

cxy = /42Q,x,

- &xx +

- &,xy), ‘+&y

1,

(7)

where A and p are Lame’s elastic constants and V2is Laplace’s operator. Using eqns (3)-(7), the equations of motion are satisfied if

0,

(8)

PT%xx + hy = 0.

(9)

/k2Q,xx + Q,yy=

E PI

PI’vt,

l-20-4 ---cX

(a) Fig. 1. Geometry and coordinate

,1 -c (b)

system of a crack moving at a constant speed c parallel to a bond-plane.

Steady motion of a crack parallel to a bond-plane

227

sr* = 1- c*/cr*,

(10) (11)

In eqns (81, (9),

/3L’= 1 - c*/cL*,

cL and cr are the speeds for longitudinal and transverse waves,

FORMULATION

CL2 = (A + 2p)/p,

(12)

CT’= p/p.

(13)

AND METHOD OF SUPERPOSITION

The faces of the moving crack are subjected to a uniform internal pressure uo; on y = 0, the stresses are continuous everywhere, and the displacements are continuous for [xl> a. The boundary conditions on y = 0 are: UY

ux

@Y

(1)= gI” =

(1)= (1)=

u;2’,

u;:) = g$ =

-go,

u(1) Y

=

(1) a(2) Y 9 gxy

0,

1x1
a.

u(*) Y 7

JxI>

(2) (Jxy,

1x1> 0.

(14) (15) (16)

where superscripted field quantities refer to boundary values taken from below (1) or above (2) and evaluated on y = 0. One method to obtain a solution is to consider representations of cp and $ in the three different regions, namely, y c-h, -h 5 y ~0 and y > 0; then match the continuous quantities at y = -h and apply the mixed boundary conditions at y = 0. In this paper the manipulations required in the procedure above are avoided by taking advantage of the superposition method in linear elasticity. The solution of the considered problem is obtained by combining the solutions of (a) a line crack which moves at a constant speed c in a homogenous material 1 and is subjected to given loadings (this is in fact Yoffe’s problem) and (b) two different undamaged materials bonded together at y = -h, subject to prescribed loading at the bond-plane. The loadings in (a) and (b) are superposed so that the stresses and the displacements are continuous at y = -h and the stresses on the crack surfaces are 0; = -uo, and uXy= 0. (a) A moving crack in an infinite homogeneous elastic material The coordinate system for this part is shown in Fig. lb. Using Fourier transform techniques, the displacement potentials satisfying eqns (8) and (9) for regions (1) and (2) can be expressed as m

where i = 1, 2, A#), defined as

cp%,Y)

=

_mAi

I

e-i& e-@lkllYldt,

(174

l//“‘(X,Y)

=

1-1 Bj(5)

e+@ e-~2l~llyI d,$

(17b)

&(r) are yet unknowns, and i = d/(-l).

c*/c;,p

p, = (I-

The constants 8, and p2 are

p2 = (1 - c*/&)“*.

(18)

Since the tractions on the crack surfaces are the only external loads and since the geometry possesses a symmetry with respect to x = 0, the symmetric and antisymmetric parts of the problem are separated. The boundary conditions on y = 0 are #)

= a$? =

XY

0,

Up = up’,

-$U$I)- u$“l = g(x)H(a

1x1> 0,

(19)

1x1),

cw

-

228

S. H. CHEN et al.

for symmetric loading and (T(l) Y = c7y = 0, a$ = a$,

-j$p- u:‘] = f(x)H(a

1x1> 0

(21)

7

- Ix),

(22)

for antisymmetric loading. The functions g(x) and f(x) defined in eqns (20) and (22) are the normal and tangential dislocation densities, and H(x) denotes the Heaviside function. If eqns (17), (18) are substituted into eqns (6) and (7), the expressions for stresses in terms of Ai and Bj are obtained. To satisfy the boundary conditions (19) and (20), we find

Bl(5) =

g(t) eiB dt )/27r(1 - /3:), I-2( 1” -a

AI(~)= C2isgn(Z)(l +

Pz?(~’ g(t) e@ dt)/4@,(1 -0

(23)

-p:),

(24)

where sgn(t) = f 1 for t P 0. With eqns (23) and (24), the normal stress on y = 0 is expressed as

where

(26) and Ml,

P2)

=

U+

(27)

b2>‘-4P~P2

is the well-known Rayleigh function[ 141. Similarly, the expression for the shear stress in terms of the tangential dislocation density f(x), satisfying boundary conditions (21) and (22) is

(28) where

k2=j$$$$.

(2%

As c + 0, with the substitutions p, i= 1 - c2/2cL2,/32= 1 - c2/2cT2eqns (25) and (28) reduce to the corresponding equations for the static case: fly(x)_

PI

1 2r(l-v1)

I“go

-at-x’

e_Yw= CL1

1

I“f(t)

27r(l- VI) -_(It-x’

Whb)

The solution of eqn (25) reproduces Yoffe’s result. The displacements and stresses at any point (x, y) can be expressed in terms of the dislocation densities f(x) and g(x); for example, the normal stress is

(31)

229

Steady motion of a crack parallel to a bond-plane

or in alternate form,

where the superscript I denotes stress due to part (a) (later II will denote the contribution from part (b)). The form of eqn (31) is more desirable when considering the continuity conditions at Y = -h. (b) Two materials without a crack bonded together at y = -h. In Fig. lc, which shows the geometry and coordinate system, the subscripts denote materials 1 and 2. The displacement potentials are chosen similar to those in eqns (17) and (18) except for a shift in the y-exponential term. Ix:

d$yf=

I

_-m (J(s)

e-blSilY+hle-i@d&

where j = 1, 2 correspond to regions 1 and 2, respectively, (p, 4) = (1, 2) and (3,4) for j = 1, 2, respectively. By substituting eqns (33) and (34) into eqns (3X7), the stresses and displacements at any point (x, y) can be written in terms of Ci and Q, j = 1, 2. Some useful expressions are recorded here: u;?(x,y = -h + 0+) = pI _: t2f(1 + ~2)2C,(~)- 2i s~(~)~2~~(~)1e-‘@ d5,

(35)

a(‘)(~ XY ) y = -h +0+) = /A,

(1+ P?PdS)l e-‘@ d&

(36)

Bd2Cd5)t 2i sgn(M&D&)] e-@ d&

(37)

I

_L

t2t2i

%n(S)BICl(t)

= ~2 I -1

152K1

+

I

o$%Y = -h -0

d2)(x XY 7 y = -h - 0”) = p2

t2[

-

2i

+

s@@)b3c2(8

+

(I+

&I?&(@]

e-“’

d&

(3%

cc

I I

f&&y

=

0) =

agx,y

=

0) = p1 _mt2[2iSgn(!%% e-81"'hCl([) t (1 + 02) e-B~SlhD,(~)]e-‘fx d&

/.Lt

_-m t21(l + fl?) e-B1’6%(5.) - 2&i sgn(e) e-~2’~‘~~,(~)]e-%Xd&

(39)

m

(4N

where 83

=

(I

-

C%ip

84

=

(1 - c2&)“2

(4la,b)

(c) S~~e~osit~o~ of (a) and (b) The unknowns C#> and L+(l) are chosen such that ou, a,,, u,, U, are continuous at the interface y = -h. Thus, for example u$~‘(x,~= -h -O+) = LT;‘)(x,Y = -h + 0+)+ q’(x,y = -h),

(42)

where oi2), cr$‘),a,’ are from eqns (37), (351, and (31). Three other equations can be similarly obtained by replacing a, in eqn (42) by o;,, uX,U, respectively. These four equations result in a set of four equations for the four unknowns Cj([)t Q(s), j = 1, 2 in terms of f(x) and g(x). Substituting solutions of C,(t) and D,(f) into eqn (39), the expressions of a: are obtained in

230

S. H. CHEN et al.

terms of f(x) and g(x) for part (b). The superposition of solutions in part (a) and (b) yield

uy=uy’+uy’r=-u~;

IxI
(43)

similarly, =ufiy+ulxly= 0; IxI
(44)

or explicitly, for y = 0,

(45) W21(~,f)f(f)

z-z,

+

Mx,~k(~)l dt, Id> 0

IxI
W)

The kernels k,,(x,t), k12(x,f), kZ1(x,f), k22(x,f) in the above singular integral equations are defined in the Appendix. In addition to the above eqns (45), (46), the following conditions for crack closure are to be satisfied. +j;

0

f(t)dt =0,

(47)

-1 a g(t) dt = 0. lr I -0

(4)

Equations (45) and (46), together with the conditions, given by (47) and (48), can be solved numerically by the method of Erdogan and Gupta [ 121. NUMERICAL

SOLUTION

AND RESULTS

Equations (45-48) are normalized by using new variables defined by:

t*=t/U,x*=x/a,h*=h/a.

(49)

The equations cited above remain unchanged except that the integration interval becomes from -1 to 1. Thus, the original notations can be retained without ambiguity. The nature of the singular behavior of f(t) and g(t) is as (1- t?-1’2. It is convenient to introduce f(t) and I(t) by

f(O

=

E(l

-

P)-l’2f(t),g(t)

By using the method described in[12], eqns (454)

= EJl-

t2)-“2&)*

(50)

may then be written as (51)

ff&[~(ti)k2l(x,,ti)t~(li)[~+ kU(xk, ti)]}=-1,

(52) (53) (54)

Steady motion of a crack parallel to a bond-plane

231

The integration and collocation points in eqns (51-54) are chosen as the zeros of Tn(ti) = 0, i=l,..., n, U”_,(&) = 0, k = 1,. . . ) r~- 1. Here, T,(x) and U,_i(x) are the nth and (n - 1)th order Chebyshev polynomials of the fist and second kind. The number of integration points, n, has to be selected properly to ensure efficient calculation and accurate solutions. As observed, n has to be fairly large for small h/a, although n can be halved by taking advantage of the symmetry and antisymmetry of functions f(x) and g(x). It is found that n/2 = 20 and 40, give sufficient accuracy for h/a = 0.1, 0.05 respectively. We choose pi L 0, j = 1,2,3,4, which ensures that the propagation speed is subsonic in both materials and hence

c < min (CTI,CTZ)

(55)

and cT1 and cT2 are the transverse wave speeds in the materials 1 and 2, respectively. The stress intensity factors for mode I and mode II are defined as Kr = t/(2) lim V(x - a)u,(x, y = 01, x+(I+

(56)

KII = V(2) lim .\/(x - a)oJx, x+4+

(57)

y = 0).

They can be expressed in terms of f(l) and g(l) as Kr = -C7&r&l)~(a),

(58)

&I =

(59)

-mkzdl)d/(4,

where p(l), g(l) may be calculated by using the following extrapolation formula[l5], after f(ti), g(ti), i = 1,. . . n are obtained, n sin

i(l)=AZ

[

y(2i .

- 1)~ 2i-1

sln [ 4n*

I

j(fi).

I

For several values of h/a and c/err and for pI = p2 and rq = y2= 0.3, the stress intensity factors K, and Kh are plotted with respect to ~21~~in Figs. 2-5. For a fixed crack speed, K,

KI

QoJ;j Y,‘Ya 5.3 2.0

Fig. 2. Stress intensity factor KI vs pJpl (&I

< 1) for various values of C/CT,and h/a.

232

S. H. CHEN et al.

(k,, 1.5

-

9=(.3,.l)

ry

P, = Pz Y, - Y2 = .3

1

Fig. 3. Stress intensity factor Ku vs p2/pI @Z/F, < 1) for various values of C/CT,and h/a.

and &I increase for decreasing h/a when ~21~1< 1, and K1 and KII decrease for decreasing h/a when ~2/pi > 1. For ~2Ipi = 1 (identical materials) Kr and KII are independent of h/a and C/Q,, (K,, KII) = (a&a, 0). It is also observed that the stress-intensity factors at large crack speeds are qualitatively similar to those at small h/a. This observation indicates that the larger the crack speed, the more intense is the effect of the top half-plane upon the stress intensity factors. In addition to the stress intensity factors, the cleavage angle at which uee is a maximum is of interest. After KI and Krr have been computed this angle can be obtained from meee = (~w)-“‘K = (2~rr)-“2[K,T’Be(0,c) + KJ’:$(B, c)]

(61)

where the expressions for I’:@ was obtained from[16], and Tii was obtained explicitly by

Fig. 4. Stress intensity factor KI vs ~d~I (p2/pLI> 1) for various values of C/Q, and h/a.

233

Steady motion of a crack parallel to a bond-plane

Fig. 5. Stress intensity factor KU vs ~z//.LI(&MI > 1) for various values of C/CT,and h/a.

following the derivations in[ 111.We note that the derivations in[l I, 161are calculated for a stress free crack. For our case this corresponds to a tensile loading, uo, at y = -+m. Thus,

L2

(62) (2 - (YT*)* (2-(Y7-3&1 COS28 - 2(1 _ (y:)&T2 + 2( 1 -

gin 28

(uL2)“*fiL2 sin 28 - [(2 - CU.*) cos 28 + (IL2- (YT*& ,

(63)

where fiTI

=

(1 - (Ye*sin*@“*- cos 1 - (YT2sin28

I/2

8 I

64

’

1’

(1 - (YT2sin*@“*+ cos 8 “* nT2=

(YT =

1 - (YT* sin*8

c/CT, (YL

=

(65)

w

c/CL.

LEGEND Q+

&T’.2

+Y

.4

*

.6

+

.o

Fig. 6. Function TL(0, c) vs 6 for various values of C/CT.

234

S. H. CHEN el al.

The expressions for &I and &Z are obtained from eqns (64) and (65) by replacing subscripts T by subscripts L. For various values of C/Q and for v = 0.3, Tg(& c) is plotted vs 8 in Fig. 6.

8, DEG Fig. 7. Function K vs 6 for C/CT,= 0, h/a = 0.4 and various values of pJpI.

-180

-120

-60

60

0

I20

8, DEG Fig. 8. Function K vs 0 for C/CT,= 0.4, h/a = 0.4 and various values of p2Ik1.

2.500

1.900 K I 300

,700

,100

-.5GO -110

-120

-60

0

60 ‘3,

120

II

DEG

Fig. 9. Function K vs tl for C/CT,= 0.8, h/a = 0.4 and various values of pp./p,.

235

Steady motion of a crack parallel to a bond-plane 1,360

S.020 K .6*0

LEGEND

,340

“El0.000

-.340 -fOO

-60

-120

60

0 e,

Fig. 10. Function X vs 0 for p&l=

I20

g-y0 r

-f9

.2

*

.4

+

.6

*

.6

/

OEG

10, h/a = 0.4 and various values of cjcr,.

K .osc

.360 t270 _ .020 -160

-I20

-60

60

0 e,

120

DEG

Fig. II. Function K vs B for clcrl = 0 and various values of hia and &J&t.

8, DEG

Fig. 12. Function K vs 8 for C/CT,=0.4 and various values of hia and &L,.

236

S. H. CHEN et al.

Fig. 13. Function K vs 0 for c/cT, = 0.8 and various values of

h/a

and gLzIp,.

Figures 7-13 show the intensity factor K, which is related to the near-tip stress uee by eqn (61), for various values of c/crI, h/a and PJ~L~. All cases correspond to pI = p2 and vI = u2= 0.3. It is observed that the cleavage angle 8 is positive for p2/p1 > 1 and negative for ~LJP, < 1. Furthermore, as the crack speed is increased the magnitude of the cleavage angle tends to increase as shown in Fig. 10. It can be noted from Figs. 11-13 that the upper material has little effect for h/a > 5. At zero speed the stress intensity factors and the cleavage angles obtained in this paper are consistent with the corresponding results of [6], provided that eqn (17) in [6] is corrected as follows:

The corrected values of the stress intensity factors corresponding to Table 1 [17], for a crack parallel to the free boundary of the cracked half plane are given in Table 1. Table 1. Stress intensity factors for a stationary crack parallel to the free surface h/a 0.2 0.4 0.6 0.8 1.0 1.2 1.5 2.0 3.0 oc

K$uoda 5.9499 2.9055 2.0810 1.7138 1.5110 1.3848 1.2677 1.1633 1.0779 1.0

KI,lUO~/a 3.0300 0.9940 0.4936 0.2890 0.1849 0.1252 0.0750 0.0368 0.0124 0.

Acknowledgement-The work reported here was carried out in the course of research sponsored by the U.S. Army Office IXuham under Grant DAAG29-78-Ga031. REFERENCES [l] A. H. ENGLAND, J. Appl. Mech. 32,400 (1%5). (21 J. D. ACHENBACH, Z. P. BAZANT and R. P. KHETAN, Inl. J. Engng Sci. 14,797 (1976). [3] C. ATKINSON, In Mechanics Fracture, Elastodynamic Crack Problems (Edited by G. C. Sib), Vol. IV. Noordhoff, L.eyden (1976). [4] M. COMNINOU, 1. Appl. Mech. 44,631 (1977). [5] J. D. ACHENBACH, L. M. KEER, R. P. KHETAN and S. H. CHEN, J. Elasticity. To be published (1979). [6] F. ERDOGAN, Engng Fracture Mech. 3,231 (1971). [7] F. ERDOGAN, Engng Fracture Mech. 4,811 (1972).

237

Steady motion of a crack parallel to a bond-plane

[8] F. ERDOGAN and G. D. GUPTA, Jnt. J. Solids Struclures 7, 39 (1971). 191F. ERDCGAN and G. D. GUPTA, Jnl. J. Solids Structures 7, 1089(1971). [ iOj E. YOFFE, Phil. Msg. 42,739 (1951). [ll] J. D. ACHENBACH and Z. P. BAZANT, 3. Appl. Mech. 42, 183(1975). [12] F. ERDGGAN and G. GUPTA, Q. J. Appl. Math. 30,525 (1972). [13] I. N. SNEDDGN, Fourier Transforms. McGraw-Hill, New York/Toronto/Londin. pp. 444 (1951). [14] J. D. ACHENBACH, Wave Propagation in Elastic Solids, North-Holland/American Elsevier, Amsterdam/New York (1973). [lS] S. KRENK, Q. J. Appl. Math. 32,489 (1975). [16] J. D. ACHENBACH, Wave Propagation, Elastodynamic Stress Singularities and Fracture, Theoretical and Applied Mechanics (Edited by W. T. Koiter). North-Holland, pp. 71 (1976). [17] F. ERDOGAN, SIAMJ. Appl. Math. 17, 1041(1%9). (Received for publication

3 July 1979)

APPENDIX In eqns (51) and (52), k,,(x,t) = 0 -xl

I2

14

I6

4~,2h2+(t-x)Z+(~,+~~)2hZ+(f-x)2+4/32th2+(~-x)i

Whl,

MB, +

1’

2jhhls

B2)h

(Al)

k'2(x~')=4~,2h2+(t-~)2+(~,+~~)2h2+(t-~)2+4~~h2t(~-x)2'

(“a

m(B~+/%)h 'V&m6 Whm kz~(x~'~=4~,2h2t(~-~)2+(~,+~~)2h2t(f-x)2+4~~~2+(t-x)*'

(W

ml

k&,t) = 0 -xl

m3

m5

4~,‘h2+(r-x)Z+(~,+,92)zh2+(t-x)2+4~~2~2+(r-x)i

1

(A4)

To obtain these results two Fourier transform identities are used,

I-

e-hltl e-4

d[

-I

= I -I

fa f

f(t)@

df = 2f:0

--(I

i sgn(l) e-h’ele-‘o dt

hZt (:_x)j(r)df,

h2:(I!X)j(r)dt’ f_-(I

’ _yf(t)eit’dr=-2

a

The coefficients Ii, mi, i = 1 - 6,aredefined as I,= 4p,D-'(fjlfj1)'.3 =4~,D-'(-f,J,,+li*Lf2,-f3J31+~41f41) I*= 4j?,D-‘(ijJj3)3.4

(A@

I, = 2D-‘{2,Y,(&fj~)‘~3+

(1 +&?($jj#*}

l4 = m-‘{2fi,(ij,fj4)3’4 15= 2(1t ,9:)D-‘(4 16= 2(1t

(.47)

+ (i +

(fw

&%4~fi3)~‘3>

(AlO)

jj#*,

(All)

/32*)0-‘(4~fj4)~.~.

(A0

m,=2(1+~~)D-‘(~~Jj~)“3~

(A13)

mz = 2(1t &*)D-‘(9fi3)‘~*,

(A14)

m3 = 2D-‘{(l +

fit)($Ji2)“‘-

2&(9 j/Z”b

W5)

m4 = 2D-‘{( 1t

&*)(ij Jj4)"* -

2&(.$

(‘416)

jj3)“99

mS= -4&D-‘(fr jrz)3’4,

(.417)

m6 = -4fizD-‘($

(‘418)

jj4)“‘.

where D is the determinant of the matrix X defined below (eqn (A19)),and lilj represents the real (or imaginary) part of the determinant of the minor matrix of X (i.e. the matrix obtained by deleting the ith row and jth column). In eqns (A7-A18), the dummy index j is summed over j = 1, 2, 3.4, and the two superscripts p. q indicate that the pth and qth terms in the summation are negative, see eqn (A7). The matrix, X is defined as

x=

It822

I

1+g*

- i2&

- 31

-1 [ iSi where i = v(- 1) in eqn (A19). IJFS Vol. 18, No. I-P

9’--Ba -p+s49

i2b

PI

iS2 1

t fi4’) - icfi4 I

$3

ifi4 -1

1 ,

(A19)

S. H. CHEN et al.

238

The matrix, F, whose elements are fii, i, j = 1, 2, 3,4 is defined as -(1+/322) _y

1+022 2p2

F=jGii&

28, -ut/%?

It& zf3, -- It822 2

-- (l+/.%? 282 1+/3z

1+/3t2 --b

1

1 61

2

-- 1+/h* 282