Dispersive waves in composites, a comparison between various laminated plate theories

Composite Structures 25 (1993) 449-457

Dispersive waves in composites, a comparison between various laminated plate theories W. J. N. Lima & A. M. B. Braga Department of Mechanical Engineering, Pontificia Universidade Cat6lica do Rio de Janeiro, CEP 22453, Rio de Janeiro, Brazil The exact Rayleigh-Lamb wave dispersion spectrum of laminated plates is obtained by applying a recursive algorithm developed for the analysis of stress waves in anisotropic layered media. Results from this benchmark problem are used in comparisons between the classical laminated plate, first-order shear deformation, and two higher-order shear deformation theories. Attention is paid to the effects of the number of layers, direction of propagation, and lay-up (symmetric or asymmetric), on the range of frequencies where these theories apply.

give evidence of the boundaries between different approximate theories. This approach has been pioneered by Mindlin 7 in the evaluation of shear correction factors for his first order theory of homogeneous, isotropic plates. Nelson and Lorch 8 and Whitney and S u n 9 have also relied on comparisons between exact and approximated dispersion curves for the calculation of correction factors needed by their higher-order laminated plate theories. Valisetty and Rehnfield l° have compared the first branch of the flexural wave dispersion spectrum obtained via a ply level analysis and other approximated theories with exact results reported by Kulkarni and Pagano. 5 These authors have studied laminates with up to five layers, and with a variety of lay-ups. Due to numerical stability problems, inherent in the method applied in the exact analysis presented in Ref. 5, the results from Kulkarni and Pagano are limited to the low-frequency range (wavelengths not shorter than twice the plate thickness). Here, the exact Rayleigh-Lamb wave dispersion spectrum is obtained by applying a recursive algorithm developed for the analysis of stress waves in anisotropic layered mediaJ 1,12 This algorithm is based on the discrete invariant imbedding technique, and has been shown to be stable over a wide range of frequencies, even in the presence of evanescent wave modes. 11 Due to the numerical stability of the algorithm, it is possible to obtain exact results for the whole frequency spectrum and for a variety of laminates.

INTRODUCTION Over the last two decades, considerable effort has been directed towards the development of approximate theories for the deformation of laminated composite plates. Thorough literature reviews are found in articles by Reddy, 1 Noor and Burton 2 and Kapania and Raciti, 3 to name only a few. In most cases, in order to assess the range of validity in the frequency domain of a given theory, comparisons are made with available exact results obtained from three-dimensional elasticity, such as those presented, for instance, in Refs 4-6. With a few exceptions, most of these results are restricted to simply-supported, rectangular, crossply laminates with a limited number of layers, and to the low-frequency/long-wavelength range. This paper presents a study on the limits of application in the frequency domain of laminated plate theories. This investigation is made via a wave dispersion analysis. The dispersion spectrum of waves propagating in a composite plate calculated via a higher-order theory, has as many real branches as the number of kinematic variables used in the approximation of the displacement field. Within the range of frequencies where the theory is valid, each of these branches must have a counterpart in the Rayleigh-Lamb wave dispersion spectrum obtained from three-dimensional elasticity. A comparison between dispersion curves calculated through both methods allows one to assess the limits of validity of the laminated plate model, as well as to 449

IV. J. N. Lima, A. M. B. Braga

450

These results are used in comparisons between the classical laminated plate theory, first-order shear deformation theory, ~3 Reddy's refined third-order theory ~4 and Lo et al. ~5 first-order theory. The governing equations of the approximate theories are written in matrix form using a state variable approach, and the dispersion curves are constructed from the solution of the associated eigenvalue problem. Attention is paid to the effects of the number of layers, direction of propagation, and lay-up (symmetric or asymmetric), on the range of frequencies where these theories apply.

the field variables into two parts, corresponding to the contributions of the upgoing and downgoing waves to the total fields. Using the subscripts 1 or 2 to represent, respectively, the upgoing or downgoing wave motions, we may write:

fi(z)=fit(z)+fi2(z),andi(z)=it(z)+L(z)

(4)

It may be shown that the fields fi~(z) and ia(z) assume the following form lt,~: (no sum over repeated index) fi~,(z) = W~,(z)fia(0), and/~(z) = - i~oZ~fi,(z) (5)

EXACT ANALYSIS The field equations governing the linear, elastic motions of homogeneous anisotropic media are div 0 = / 9 02u/Ot 2, and o = C:Vu

(1)

where o represents the Cauchy stress tensor, u the displacement field, and C the rank-4 elasticity tensor. Since we are concerned with laminated planar structures, it is convenient to employ as a field variable the traction t = ae z, which must be continuous across planes parallel to the layering. Assuming plane harmonic motions in the form

f(x, z, t)=f(z) e ''~..... '! (2) where f(x, z, t) represents any of the field variables, while w is the angular frequency and k x the wave-number in the x-direction, eqns (1) may be recast into a six-dimensional system of ordinary differential equations, written in matrix form as t 6

Iu/z)l

ddx~(z)--iN~(z)'where~(z)=|ii(z)J

(3)

while the linear operator N is a real valued sixtensor depending upon the anisotropic material properties, and on the pair (w, kx). The matrix N for orthotropic media may be found in Ref. 16. In order to obtain a solution for eqn (3), one must find the eigenvalues of N. These six eigenvalues are the wave-numbers in the z-direction, and occur in pairs representing waves propagating or decaying in opposite directions of the z-axis. 16 We use the terms upgoing and downgoing when referring to these partial waves. The eigenvectors of N are the polarization vectors of these six waves. A common practice in the description of the wave motion in layered media is the separation of

where the Wa(z) and Z , ( a = l , 2 ) are ( 3 x 3 ) matrices constructed from the eigenvalues and eigenvectors of the six-dimensional system matrix N. Equations (4) and (5) completely describe the plane harmonic motions in homogeneous anisotropic media. The operator W~(z) is the propagator matrix that relates the upgoing (a = 1 ), or downgoing (a = 2), displacement field at a plane of coordinate z with its value at z = 0. The matrix Z~ is the local impedance tensor of the upgoing (a = 1) or downgoing (a =2) waves. The vectors fit(0) and fia(0) are the only unknown terms in the formalism represented by (4)and (5), and must be obtained from boundary conditions applied at planes normal to the z-axis. Surface impedance tensor of a homogeneous plate We now consider an infinite, homogeneous anisotropic plate of thickness hi, which is bonded along its lower surface ( z = 0 ) to a substrate characterized by the impedance tensor G 0, i.e. along z = 0 - (the plane immediately below the interface) i(0 )= -itoG0fi(O-)

(6)

We assume that G Ois given. If the layer is tractionfree at the lower surface, we take G0=0. In general, the impedance Go will depend on both w and k x. Since the layer is perfectly bonded to the substrate, both the displacement and traction must be continuous along the interface, i.e. fi(O-)= fi(O+ ), and i(O- )= i(O +)

(7)

In this case, tractions and displacements at z = h I , the top surface of the homogeneous layer, will be

Dispersive waves in composites related through the expression 11,12 /(h~) = -itoG(hl)li(h~)

(8)

where

G(hl)=[Z1H(hl)+ Z2] [I+/-/(hi)] -~

(9)

451

solving the eigenvalue problem for the N matrix of each homogeneous layer. A detailed analysis of this algorithm may be found elsewhere, 11 where its stability over a wide range of frequencies has been verified. Rayleigh-Lamb waves in a composite plate

and H(hl)=

W l ( h l ) [ Z 1 - G 0 ] - l [ G 0 - Z2] W~-'(h,)

(10) Equation (8) expresses a linear relation between the traction and velocity fields at the plane z = hi. The (3 x 3) operator G(hl) represents the surface impedance tensor of an anisotropic layer of thickness h 1bonded to a substrate with impedance G 0. Surface impedance tensor of a multilayered composite plate

We now consider a laminated composite plate which is in contact along its lower surface with a substrate of given surface impedance Go. The n homogeneous layers have thicknesses hi, h2,..., h~, and the total thickness of the laminated plate is h. The goal is to find an expression which is a generalization of (8) to the case of a multilayered plate. That is, we need to find the tensor Gn such that i = - itoG, ti

(11)

where i and ti are, respectively, the traction and displacement on the top of the plate. We start by using the results of the last section to evaluate the impedance G1 on the top of the first layer (for the sake of simplicity, we write Gj = G(hj); likewise, we write Hj= H(hi) ). We then proceed by employing this result to calculate G2, the impedance on the top of the second layer, and so on, up to the nth layer. This procedure is summarized in the following algorithm:

In this section, we consider an infinite, laminated composite plate which is traction-free on both surfaces. We are interested in generating the dispersion curves of plane waves that can propagate freely in the plate. In order to apply the algorithm (23), we choose Go = 0, which, according to eqn (6), guarantees that the lower surface of the plate is free of tractions. For each (to, kx) pair, we apply (23) to evaluate the impedance tensor G, on the upper surface of the plate. From eqn (11), we conclude that the traction-free boundary condition will be satisfied only for those values of to and k x such that the following equation holds: det Gn=0

This equation provides an implicit relation between the frequency to and the wave-number kx. It represents the dispersion equation of free waves that propagate in the plate in the absence of an external disturbance. The numerical results of eqn (13) will be compared with those predicted by different approximate plate theories. APPROXIMATE THEORIES

The highest order laminated plate model considered in this paper is the third-order shear deformation theory introduced by Lo et aL 15 (TSDT), based on an assumed displacement field which, for time-harmonic motions, has the form

ux(x, y, z, t)=[u(x, y)+ Z¢l(X, y) + z2v;l(X, y) + z301(x, y)] e -iwt

tly(X, y, z, t) =[v(x, y) + zfJ2(x, y) + z2g'2(x, y)

Given Go, to and k~;

+z302(x,y)]e -i~''

For j = 1 to n repeat HI .= W (()(hi) [Z(( ) -

(13)

Gj_1]-1[ G]_I

(14)

Uz(X, y, z, t) =[w(x, y) + ZCa(X,y)

-

+ z2~)3(x, y)] e -i~t

x(Wj/(hj)) -1 n,.]- '

(12)

where 14~l(hj) and Z~) are the propagator matrix and the local impedance of the upgoing (a = 1 ), or downgoing (a = 2), waves propagating in the jth layer. These (3 x 3) matrices are evaluated by

Other theories may be considered particular cases of (14). The classical laminated theory (CLT), is obtained from expression (14) by assuming

¢l=-Ow/ax, and Oi=O for i= 1,2,3.

¢2=-Ow/ay,

~3 --" 0,

~)i --'--0,

(15)

452

W. J. N. Lima, A. M. B. Braga

The first-order shear (FSDT) corresponds to

V)i=O,

¢2=0,

deformation

theory j3

and

(16)

0 i=O

while Reddy's refined third-order (RSDT) is obtained when we let

4 ( ¢, +aw/

v,,=o, o,-

¢~ = 0,

and

theory 14

(]7)

02 =

1.2

4[ / ] 3h 2 ~¢2 +Ow Oy/

1.6 II

L) r"

3

0.0

- I

I

I

I

I

I

I

I

I

I

I

I

1

I

i

l

I

I

I

0.0

T

I

I

I

I

I

I

I

1

I

I

(a) 2.0

I

I

I

I

I

I

I

2.0 o m

-

o| |

•

no

-u

I

(d)

oono °u

--

oo o.

°•

•

t-

3

0.0

T-]

1

I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

o.o

I

-i

i

i

i

i

i

i

i

i

(b)

~

2.0

K c-

°°

o °°

i

i

i

i

i

i

i

i

i

(e) 2.0

f

._,:."

i

3

0.0

0.0

0.0

Fig. 1.

0.5

1.0

7 0.0

I

I

I

I

I

I

I

I

I 0.5

x:., h/'n"

x,, h/"rr

(c)

(f)

I

I

I

I

I

I

I

I 1.0

Dispersion curves for a symmetric, cross-ply laminate with 15 layers: (a) CLT, ~f=0°; (b) FSDT, cp=0°; (c) RSDT, cf = 0°; (d) CLT, q9= 60°; (e) FSDT, q~= 60°; (f) RSDT, q~= 60 °.

Dispersive waves in composites Inertia coefficients and generalized forces are defined, respectively, as

Ik =

f hl2 d-h~2

Ch/2

PZ k- 1 dz,

and

FSDT and RSDT. The matrix M for TSDT may be found in Ref. 17. From the symmetries of $ and L follows

R~k/=] ojZ k dz 3 -h/e (18)

where oj(] 1, 2,..., 6) are the six components of the stress tensor. Stiffness coefficients are defined by =

./M3"r=-M r, where

f h/2 -h/2

Qij(1, z, z2, z3, z4, zS, z6)dz

;I]

(26)

4.0

(19)

where Q0 are the elastic constants in the plate coordinates -- for CLT, FSDT and RSDT, they represent the plane-stress-reduced constants. In the case of cylindrical bending, when all derivatives in the y direction vanish, the balance of linear momentum -- in the absence of forcing terms -- and constitutive equations for all four theories may be written in matrix form as 17 d d--x ~ = M ~

j=[Ol

From (26), we conclude that the eigenvalues of M, which we denote by ik x, occur in pairs of opposite

(Aij, Bij, Dij, Eij, Fij, Gij, Hi]) =

453

>

0.0

4.0

(20)

where the state vector, ~, is defined as

~T=[u, v, - d w / d x , w, RI°), R~6°t,R~'), R~°)] for CLT (21)

'///

(a)

13 ¢-

3

CT=[U, v, ~bl, ¢2, w, Rt°), RC6°),Rql), R~611,R~°1] for FSDT (22) ~T= [U, V, O1, qk2, - 4( dw/dx)/3 h2, w, Rt °1, R(6°1,

R~I) - 4RI3)/3h 2, R(61)- 4R(63)/3h 2, R~3), R~°)] for RSDT (23)

0.0 u

I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

(b)

4.0

and U = [u, v, ~,, ¢2, W,, ~2, 0,, 02, w, ~3, ~03, RId'l, R~6°/, Rt 1), R(61),Rt2/, R(62),R~3), R(63), R~°/, R~'/, R~2/] for TSDT

(24)

Even though its size depends on the theory that is being applied, the system matrix Min eqn (20) can in all cases be partitioned into four real, square submatrices: 0.0

0.00

3.00

,cxh/n (c)

where both S and L are symmetric. These matrices are presented in the Appendix for CLT,

Fig. 2.

Dispersion curves for a (0°/90°/0°/90°) laminate: (a) CLT, ~0= 00; (b) FSDT, q0 = 0°; (c) RSDT, q0 = 0 °.

454

W. J. N. Lima, A. M. B. Braga

sign. 17 It is clear that these eigenvalues are the wave-numbers associated with partial waves that propagate -- or decay when kx is complex -along the unloaded infinite plate. When solving a boundary value problem, these partial waves must be linearly combined in order to satisfy a proper set of boundary conditions. Next, we compare pairs formed by these eigenvalues and frequencies with the dispersion curves of Rayleigh-Lamb waves evaluated through the exact method presented in the preceding section.

makes an angle q) with the stiffer direction of the top layer in the laminate. In all plots, dotted Fines correspond to approximate theories, while continuous lines are reserved for exact solutions. The following material properties have been used in all the calculations: E I = 16E:, GI2= Gl.~=0-57E2, G23 = 0 . 3 7 E 2 and vj2=0.25. The curves drawn in Fig. 1 are for a cross-ply, symmetric laminate with 15 layers. Figures 1(a-c) have been obtained for cp--0 °. In this case, the lower branch of the dispersion spectrum corresponds to a flexural wave, the second to an antiplane shear wave, and the third to an extensional wave. All three start from the origin. The fourth branch is associated with another antiplane shear wave and the fifth corresponds to the second flexural wave. Plots in Figs l(d,e) are for 9 = 60 °. As expected, the flexural branch from CLT (Figs

RESULTS AND DISCUSSION The dispersion curves presented in this section have been obtained by plotting the nondimensional frequency, o)h/ng, versus the nondimensional wave-number kxh/n, where ~= , J Q S p is a reference wave speed. The wavelength is ,~ = 2n/k~. Except for the slowness plots in Fig. 4 below, the direction of propagation (x-direction) 4.0

"~

~

I

~

I

0.0 4.0

0.0 0,00

3.00

(xh/~ (c) Fig. 3.

0.00

3.00

(.h/TT (d)

Dispersion curves obtained via TSDT: (a) (0°/90°/0°), ~ = 0°; plane motions; (b) (0°/90°/0°), (p = 0 °, antiplane motions; (c) (0°/90"/0"/90"), 9 = 0°, plane motions; (d)(0°/900/0°/900), to = 00, anfiplane motions.

Dispersive waves in composites

the one obtained from RSDT and shown in Fig. 2(c), whose first flexural branch matches the exact curve up to the nondimensional frequency 0.45, corresponding to a ratio 2/h = 2. The first two plots in Fig. 3 present the dispersion spectrum of a (00/900/0 °) plate calculated from TSDT. The results are for q~= 0 °. In this case, plane and antiplane motions uncouple, and are plotted respectively in Figs 3(a) and 3(b). The first three branches of Fig. 3(a) and the first two of Fig. 3(b) can be identified with those obtained from both FSDT and RSDT. The results from these two theories for the three-ply laminate are not shown here, but we have observed that the agreement of the lower branches of TSDT with exact dispersion curves is better in this case. As the frequency increases, the matching of the higher branches of TSDT with those obtained from the exact analysis is very poor though, even for low wave-numbers (long wavelengths). The

l(a) and l(d)) matches the exact result only at very low frequencies. The shear-correction factor used for the FSDT results presented in Figs l(b) and l(e) has been calculated following Mindlin's approach7 We observe from Figs 1(c) and l(f) that, as the wave-number becomes larger, the RSDT will predict higher frequencies for the lower flexural mode than the exact analysis or FSDT. The results in Fig. 2 are for an antisymmetric laminate with four plies (0°/90°/0°/90°). The angle is 0 °. Again, except for the first antiplane shear mode, results from CLT match the exact dispersion spectrum only in the low-frequency region. The shear correction factor :r2/12 was used for FSDT results in Fig. 2(b). The matching of the first three branches in Fig. 2(b) to exact results is remarkably good up to the nondimensional wavenumber 3.0, which corresponds to a wavelength equal to 2 h / 3 . This is a much better result than

2.0

O.

--

2.0

--

0.0

0.0

-2.0

-2.0

--

(a) 2,0

CL

0.0

(b) 2.0

--

0.0

- -

-2.0 --'

-2.0

-2.0

Fig. 4.

455

0.0

2,0

-2.0

0.0

px

px

(c)

(d)

2.0

Slowness curves for a (0°/45°/0°/45°/0 °) laminate at cohln(= 0-8: (a) CLT: (b) FSDT; (c) RSDT; (d) TSDT.

456

W. J. N. Lima, A. M. B. Braga

same conclusions regarding the TSDT results are reached for the plots in Figs 3(c) and 3(d), corresponding to a (00/900/00/90 °) plate. Finally, we consider a (0°/45°/0°/45°/0°) laminate and observe the matching with exact results as we vary the angle of propagation go. The plots in Fig. 4 correspond to slowness curves for a fixed frequency. They are obtained by dividing the nondimensional wave-number by the frequency, which yields the inverse of the phase velocity. Results are for wh/ar(=0"8. The slownesses Px and Pr represent the projections of the inverse of the phase velocity, kxd/0), on the x- and y-axes. The x-axis is parallel to the stiffer direction of the top layer in the laminate. The outermost and innermost curves correspond, respectively, to the first and second flexural modes. For all four theories considered, the agreement with exact results depends on the angle of propagation. Poorest results are again from the CLT. At this frequency, both the FSDT and RSDT provide very good estimates of the first flexural phase velocity for angles of propagation close to 90 °. The best match is obtained from TSDT. The matching of the second flexural mode is poor for all three theories, except for angles go close to 0 °.

REFERENCES 1. Reddy, J. N., A review of the literature on finite-element modeling of laminated composite plates and shells. Shock Vib. Digest, 17 (4)(1985) 3-8. 2. Noor, A. K. & Burton, W. S., Assessment of shear deformation theories for multilayered composite plates. Appl. Mech. Rev., 42 (1)(1989) 1-13.

3. Kapania, R. K. & Raciti, S., Recent advances in analysis of laminated beams and plates, part lI: vibrations and wave propagation. AIAA J., 27 (1989) 935-46. 4. Srinvas, S., Rao, J. & Rao, A. K., An exact analysis for vibration of simply-supported homogeneous and laminated thick rectangular plates. J. Sound Vib., 12 (1970) I87-99. 5. Kulkarni, S. V. & Pagano, N. J., Dynamic characteristics of composite laminates. J. Sound Vib., 23 (1972) 127-43. 6. Dong, S. B. & Nelson, R. B., On natural vibrations and waves in laminated orthotropic plates. J. Appl. Mech., 39 (1972) 739-45. 7. Mindlin, R. D., Influence of rotatory inertia and shear on flexural motions of isotropic elastic plates. J. Appl. Mech., 18 (1951) 31-8. 8. Nelson, R. B. & Lorch, D. R., A refined theory for laminated orthotropic plates. J. Appl. Mech., 41 (1974) 177-83. 9. Whitney, J. M. & Sun, C. T., A higher order theory for extensional motion of laminated composites, J. Sound l/ib., 30 (1973) 85-97. 10. Valisetty, R. R. & Rehnfield, L. W., Application of ply level analysis to flexural wage propagation. J. Sound Vib., 126 (1988) 183-94. 11. Braga, A. M. B., Wave propagation in anisotropic layered composites. PhD thesis, Stanford University, 1990. 12. Braga, A. M. B. & Herrmann, G., Free waves at a fluid/ layered-composite interface. In Elastic Wave Propagation and Ultrasonic Nondestructive Evaluation, ed. S. K. Datta, J. D. Achenbach & Y. S. D. Rajapakse. Elsevier, Amsterdam, 1990, pp. 171-6. 13. Whitney, J. M. & Pagano, N. J., Shear deformation in heterogeneous anisotropic plates. J. Appl. Mech.. 37 (1970) 1031-6. 14. Reddy, J. N., A simple higher-order theory for laminated composite plates. J. Appl. Mech., 51 (1984) 745-52. 15. Lo, K. H., Christensen, R. M. & Wu, E. M., A high-order theory of plate deformation, Part 2: laminated plates. J. Appl. Mech., 44 (1977) 669-76. 16. Braga, A. M. B. & Herrmann, G., Floquet waves in anisotropic periodically layered media. J. Acoust. Soc. Amer., 92 (1992) 1211-27. 17. Lima, W. J. N., Dispersive waves in laminated composite plates (in Portuguese). MSc thesis, Pontiffcia Universidade Cat61ica do Rio de Janeiro, 1993.

APPENDIX: MATRIX M

For CLT: Pi/= 0 (i, j = 1,..., 4), except for/'43 = - 1 Lll = L22 = L44 --- _ 0)211, L33 = - 0)213, LI3 --- L31 = _ 0)212 L12

=

L21 = L14

=

L41 = L23 = L32 = L24 = L42 = L34 = L43 = 0

$44 = Sa4 = Saa = 0, Sa~ = K,~a for a, fl = 1, 2, 3, where

K-l=

I

-All

Al6

Bll ]

AI6

A66

BI6 •

[BII

B16

D11J

Dispersive waves in composites

457

For FSDT:

Po=O (i,j = 1,..., 5), except for P53 = - 1 and P54 = -A45/A55 Lil = L22 =L55 = - 09211 , L33 = - 09213, L4a = _ 0)213 + A 4 4 -AZ5/A55 L13 = L31 = L24 = L42 = - 09212, L12 = L21 = L14 = L41 = L15 = L51 = 0 L23 = L32 = L25 = L52 = L34 = L43 = L35 = L53 = L45 = L54 = 0 and

S-1=

-All

AI6

Bll

BI6

0 -

AI6

A66

BI6

B66

0

B11

BIn

Dll

DI6

0

B16

B66

DI6

D66

0

0

0

0

0

A55

For RSDT:

P~/= 0 (i,j= 1, ..., 6), except for P65 = - 3h2/4 LI1 = L22 = L 6 6 = _ 09211 , L33 = _ 092[3 -k-~t55 , L44 = - 09213 + A 4 4 L55 = - 09217 + 3h:fi.55/4, L13 = L31 = L24 --- L42 = -- 092[2 Ll5 = L51 = _ 09214 , L34 = L 4 3 =A45 , L35 =L53 = - 092(s - 3hZA55/4 L12 = L21 = L14 = L41 = L16 = L61 = L23 = L 3 2 = L 2 5 =L52 = 0 L26 = L62 = L36 = L63 = L46 = L64 = L56 = L65 = 0 $66 = Sa6 = S6a = 0, Sai1 =

Ka~f o r

a , fl = 1,..., 6, w h e r e

- A 11

A 16

Bll

/~16

E~ l-

A16

A66

/~16

B66

El6

Bll

K-i=

m

1311

FI1

B16

B66

/~16

/)66

El6

Ell

El6

El1

El6

Hl l w

and Bij = Bi/- 4DiJ 3h2, 1~/= D i / - 8Fii/3h2 + 16HiJ9h 4, ~ / = Fij - 4Hr/3h 2, for i,j = 1, 6, fta~ = Aa~ - 8 D a J h 2 + 16Fa~/h 4, for a, f l = 4, 5, [2 = 12 -414/3hz, [~ = 13 - 8Is/3h ~ + 1617/9h 4, is = 15 -417/3h 2.