An efficient approach for numerical simulation of flows in Czochralski crystal growth

N

j. . . . . . . .

ELSEVIER

CRYSTAL GROWTH

Journal of Crystal Growth 181 (1997)427 436

An efficient approach for numerical simulation of flows in Czochralski crystal growth C. Shu*, Y.T. Chew, Y. Liu Department of Mechanical & Production Engineering, National Universi~ of Singapore, 10 Kent Ridge Crescent, Singapore 119260, Singapore

Received 3 February 1997; accepted 28 April 1997

Abstract

The global method of generalized differential quadrature (GDQ) is applied to simulate the flows in Czochralski crystal-growth configuration. The GDQ solver is validated by its application to the benchmark problem suggested by Wheeler. It is demonstrated in this study that the GDQ method can obtain accurate numerical solutions using just a few points and needing small computational resources.

I. I n t r o d u c t i o n

It is well known that single crystals form the foundations of modern technology. G o o d single crystals are essential for a variety of scientific and commercial purposes. They are needed for scientific appraisal of crystallography, topography and tensor properties of all crystalline materials. Basically, there are three major prototypical systems for melt-crystal growth, namely, Czochralski (CZ) growth; vertical Bridgman growth; and floatingzone growth. Among them, the CZ growth has received the most attention due to its ability to provide large single crystals. Unfortunately, the large size of a typical crystal puller and the presence

*Corresponding author. Fax: + 65 779 1459; e-mail: [email protected].

of the meniscus make the simulation of this process perhaps the most difficult. So far, the simulation of flows in CZ crystal growth is usually conducted by low-order finite difference, finite element and finite volume methods [1-11]. A comprehensive review for these solvers has been given by Langlois [4] and Brown et al. [9]. Recently, Raspo et al. [12] successfully simulated the problem by using the spectral multidomain technique. In general, the low-order methods can obtain accurate numerical results by using a large number of grid points. This would lead to the use of lot of computational effort and virtual storage. As will be shown in this paper, the global method of generalized differential quadrature (GDQ) offers an efficient way to obtain accurate numerical results by using a considerably small number of grid points and thus requiring very small computational effort and virtual storage. The G D Q method was

0022-0248/97/$17.00 @ 1997 ElsevierScienceB.V. All rights reserved PIl $ 0 0 2 2 - 0 2 4 8 ( 9 7 ) 0 0 2 9 6 - 0

428

C. Shu et h i . / J o u r n a l

o f C r y s t a l G r o w t h 181 (1997) 4 2 7 436

presented by Shu and Richards [13, 14] to improve the differential quadrature (DQ) method [15] in computing the weighting coefficients. G D Q approximates any spatial derivative at a discrete point by a linear weighted sum of all the functional values in the whole domain. Under the analysis of a highorder polynomial approximation and the analysis of a linear-vector space, G D Q computes the weighting coefficients of the first-order derivative by a simple algebraic formulation without any restriction on choice of grid points, and the weighting coefficients of the second- and higher-order derivatives by a recurrence relationship. Previous applications of G D Q method to solve the fluid flow problems [16-18] and the structural and vibration problems [19 23] demonstrated that accurate numerical results can be obtained by using very few grid points with tiny requirement of computational effort and virtual storage. In this paper, the benchmark problem suggested by Wheeler [24] will be used to validate the performance of G D Q method in the simulation of flows in CZ crystal growth. The G D Q results will be compared with available data in the literature.

(xi, y j) as N

f~(")(xi, y j) = ~ W!~)f(Xk, y j),

..... N - l ,

m-l,2

..... M - l ,

M k=l

(2) for i = 1, 2 . . . . . N;j - 1, 2 ..... M, where N, M are the number of grid points in the x and y direction, respectively,,,Vik, (") ,-(m) Wjk are the weighting coefficients to be determined as follows: Weighting coefficients for the first-order derivative

(xi A(1)(xi) - xj)A(1)(xJ) '

,.,(1) "iv =

w h e n j :# i,

v

(3)

V /,

w h e n j = i,

,rVik (1)

k = 1,k¢i

~ r i , j = 1,2 ..... N, li

2. Generalized differential quadrature (GDQ) ,7(1)

rvij

G D Q approach was developed by Shu and Richards [13, 14] based on the differential quadrature (DQ)technique [15]. It approximates the spatial derivative of a function with respect to a space coordinate at a given grid point as a weighted linear sum of all the functional values at all grid points in the whole domain of that space coordinate. The computation of weighting coefficients by G D Q is based on the analysis of a high-order polynomial approximation and the analysis of a linear-vector space. The weighting coefficients of the first-order derivative are calculated by a simple algebraic formulation, and the weighting coefficients of the second- and higher-order derivatives are given by a recurrence relationship. The details of G D Q method can be found in Ref. [13]. Some two-dimensional results are described as follows. For a smooth functionf(x, y), G D Q discretizes its nth order derivative with respect to x, and the ruth order derivative with respect to y, at the grid point

n=l,2

k-1

z

B(1)(Yi) Yi - Yj)B(1)(Yj) '

when j =fi i, (4)

M

--

~'

,7,(1)

/" ~

~

ik

w h e n j = i,

,

k=l.k¢i

for i,j = l, 2 ..... M, where N

A('(x3 =

[[ j

(x~- xj),

1,j~i M

Bm(yi) =

I~

(Yi -- Y j),

j= 1,j¢i

Weighting coefficients for the second- and higher-order derivatives

=

,~ ( , ,(1), (n ¢t ~v~,ij N vv ii

i

1,

Wij ~, ,(n 1 ) X i -- X j/

when i ¢ j, when i = j,

(5)

429

C. Shu et al. / Journal of Oystal Growth 181 (1997) 427 436

f o r i , . j = 1,2 ..... N , n = 2 , 3 WII ) :

-(1) 7,(m W i • ~&'ii

i(; ITl

[t

2

k = 1,kvLi

1)

..... N - - l ,

}vij ~i =

, Yjj

when i ¢.j, when i = j,

(6)

for i,j = 1,2 .... , M , m = 2,3 . . . . . M - 1. It is obvious from above equations that the weighting coefficients of the second- and higherorder derivatives can be completely determined from those of the first-order derivatives. Thus, the key computation is the weighting coefficients of the first-order derivative. It should be indicated that the application of G D Q is very easy and straightforward. Actually, when the coordinates of grid points are known, the weighting coefficients for the discretization of derivatives can be easily calculated from Eqs. (3) (6). Then using Eqs. (1) and (2), all the spatial derivatives can be discretized using a similar form. The difference for the respective derivatives is to use different weighting coefficients, which are usually computed in advance. This avails as an easily implementable scheme on the computer with greatly simplified code-editing features. G D Q is really a computer method, which makes the code-editing easier. The low-order finite difference method is very familiar to many researchers in CFD. As we know, the form for discretizing different derivatives at different grid points on a nonuniform mesh are quite different. Even for the simplest case where the mesh is uniform and only the first-order derivative is considered, the discretization form at the boundary point is different from that at the interior point. One should consider all these cases in the computer code. In contrast, G D Q uses the same form for discretizing different derivatives at different grid points on a nonuniform mesh. One standard "do loop" can cover the discretization at all the grid points including the interior and boundary points. As a result, G D Q makes the computer code shorter and clearer.

The configuration of Wheeler's benchmark problem is shown in Fig. 1, where Rc and H are the radius and height of a vertical cylindrical crucible. The melt is bounded above by a coaxial crystal of radius Rx < Re. The crucible rotates with an angular velocity Qc while the crystal rotates with an angular velocity Qx. The nondimensional governing equations can be written in the cylindrical coordinate system as [5]

~(ru)

~w + ~ z = 0,

(7a)

~U ~b/ OU ~t -~- U ~ -t- W ~Z

U2 F

~p U ~r -~- V2u -- --r2,

Ov Ov ~v uv & +u~+W~z -+ r = v2v

a 7 -~- b/ ~-F -]- ~' 82"

v 1"2,

(7c)

aZ -~- V2W ~- G r T,

(7d)

8T 8T ~T 1 ~ t + u ~ + w ~z - Pr V 2 T '

(7e)

where V 2 = (1/r)(~/Or)(r ~/~r) + (~Z/~z2) is the Laplacian operator, (u, v, w) are radial, azimuthal and

Crystal rotation z

Phase b

o

u

n

Free surface

~

/

/

Gravity H

Melt

• 3. Governing equations and numerical discretization To validate the efficiency of G D Q method, the Wheeler's benchmark problem is chosen for study.

(7b)

I

~r

Rc

Crucible rotation

Fig. I. The configuration of Czochralski crystal growth.

430

C. Shu et al. /Journal of C~stal Growth 181 (1997) 4 2 ~ 4 3 6

axial velocity components, (r, z) are polar coordinates, t is the time, p pressure, T temperature, v kinematic viscosity, g gravitational acceleration, and p is the density of the melt. The boundary conditions are given by

03--

8w u = v - 0r

u-

~T ~r

u=w=0,

0,

v=R%,

for r = 0, 0 ~< z ~< e, T=I,

(8a)

f o r r = 1,

0 ~< z ~< ~,

(8b)

0T u=w=--=0, 5z

v=rRe~,

for0~
(8c) T=0,

u=w=

v=rRe~,

forO<~r<~fl, z = : ¢

By introducing the vorticity co 5u Oz

5w Or'

(9)

and the stream function 0 1 ~0

r 5z'

1 00 r &"

w-

(10)

Equation system [Eqs. (Ta), (7b), (7c), (Td) and (Te)] can be converted into the vorticity-stream function formulation. In order to make the computation more stable, Langlois [4] suggested using f2 = rv to replace the azimuthal component in the equation governing v, and using S = co/r rather than o) as a dependent variable in the vorticity equation. So the final equations can be written as 5Q 10 0 1 O [ r 3 ~ {Q'~q ©~ + (ruf2) + (wf2) =

r~

(8d)

~z

-

02~c~

r~L ~\;~)1 + az2' (lla)

r--fl 1 - fi'

u=v=w=O,T-

forfi<~r~l,z=:¢.

~7 +

(8e) The nondimensional parameters that occur in Eqs. (8a), (8b), (8c), (Sd) and (8e) are the aspect ratios

r ~r

Oz

=_Grl~T

(wS)

-- ~z

1 11!0

,7~r + 7 Ur

7

]

52S

~ (r2s) + ~~z'

5T 1 5 0 ~ + -r~r (ruT) + ~z ( w r )

Rx

H

~=L,

~=~,

Pr\rer r ~ / + ez2j, the Reynolds numbers, Re~

-

(llb)

8(100) ~rr 7 ~ r

R2f20

2 RoY2~ ,

R%

-

V

Gr

(lid)

y2

Accordingly, the boundary conditions are changed to 530 5T 0 -- 0r 3 -- O = ~0t = 0,

gfl(T~ - T~)R~ =

1020 + r 5z 2 - r S .

,

and the Prandtl and Grashof numbers v Pr = ~,

(llc)

O2W S-

~r2,

(12a)

f2 = r2Re~,

(12b)

'

on the symmetry lines, where fi is the coefficient of volumetric thermal expansion of the melt, ]Pc is the temperature of the crucible wall, T~ is the crystal temperature. In this study, the aspect ratios and the Prandtl number are fixed at = 1.0, f l = 0 . 4 ,

Pr=0.05.

00 0--5z - T=O,

S-r5

1 0u z,

on the growth interface,

0 - ~ 0 = s=°--°=o, 5z

5z

T = (r - fi)/(1 - fi), (12c)

C. Shu et al.JJournal cfCysta1 Growth 181 (1997) 427-436

on the free surface,

on the crucible $ = g

Eqs. (12a), (12b), (12c), (12d) and (12e) gives two boundary conditions for $ at each boundary. After numerical discretization by GDQ method, these two boundary conditions can be combined to give two-layer conditions for $ as

bottom,

= S = 0,

fi = r2Re,,

T=

$‘l,j =

1.

dQ..

1 N

dt

+-’‘ik=l

$N,j

=

‘hi.1

=

‘/‘i,M

=

0,

(144

(12e)

on the crucible wall. Applying the GDQ method to discretize the spatial derivatives in Eqs. (1 la), (11 b), (1 lc) and (1 Id) gives rJ

431

*2,j

X

&

=

k=

1 k$2

.

ry_

1 (C\%#,‘N-

1 -

C!$,‘kC\:bp

l)$k,j,

(14b)

..

M

cl:‘rkUkj~2kj+ C

~~:‘Wiks2ik

k=l

(14c) k=l

k=l

k=l

(134

dSi.

M

1 N

(14e) for i = 1, 2, . . , N; j = 1,2, . . . , M, where

k=l

(13b) dT,, dt

M

1 N C{,$‘rkUkjTkj

+-=rik-l

+

1

AYM

CS:‘M?ikTik

(13c)

k$l

C!kz)$kj

= c’&

1Cg!2 - C\l,icc!M- 1,

k=l

-

i

t tk-1

Ci:)$fkj

+

z k=l

$)$/ik

=

rf’Sij,

and c\:k is the weighting coefficient related to the approximation of a3/dr3 at r = 0. Eqs. (13a), (13b) and (13~) can be further reduced to algebraic equations after discretizing the time derivative by Euler implicit difference scheme. The resultant algebraic equations are then solved by using SOR iterative method.

(134 where cjj’, cjj”’ are the weighting coefficients related to a/& and a2jar2 while c!:‘, 513’ represent the weighting coefficients related to a/az and a2/az2, N, M are the number of grid points in the r and z direction, respectively. Similarly, the derivatives in Neumann boundary conditions can also be approximated by GDQ scheme. It is indicated that

4. Results and discussion To validate the efficiency of GDQ method in numerical simulation of Czochralski melt flows, three cases of Wheeler’s problem are chosen for study. These three cases have the same geometry which is shown in Fig. 1, but they have different

C. Shu et al. /Journal of Crystal Growth 181 (1997) 427 436

432

flow conditions which are given as follows: Case A: Gr = 0, Rec = 0, Rex = 100; Case B: G r = 0, Rec = - 25, Rex = 100; Case C: G r = 105, Ree = 0, Rex = 0. It is noted that cases A and B are forced convection problems, case C is a natural convection problem. The flow in case A is driven by rotation of the crystal while the flow in case B is generated by opposite rotations of the crystal and the crucible. In the following, the minimum and maximum values of stream function denoted by ~min and ~max will be used to compare the G D Q results with available data in the literature for the above three cases. Since the analytical expression of stream function is unknown, the values of @rain and ~max are usually determined through the comparison of functional values at mesh points. Thus, to obtain accurate values of @min and ~ . . . . the mesh should be very fine (very large mesh size). As for the G D Q simulation, the numerical results are generated on a very coarse mesh which may not give accurate values of @rain and 0 . . . . To improve this, the following Lagrange interpolation scheme is applied to give the functional values on a fine mesh after the numerical solution on the coarse mesh is obtained, N

f(r, z j) = ~, f(rk, Zj)Sk(r),

(15a)

k-1 M

f(rl, z) = y" f(ri, zk)tk(z),

(lSb)

k-1 N

f(r, z) = y" i--1

M

~ f(ri, za)si(r)tj(z),

(15c)

j--I

where si(r), tj{z) are the Lagrange interpolation polynomials in the r and z direction, respectively. It should be indicated that the use of Eqs. (15a), (15b) and (15c) will not lose the accuracy of numerical solution since G D Q method is also based on the high-order polynomial approximation. Table 1 shows the comparison of computed minimum and m a x i m u m stream function for cases A, B, and C. The mesh size used for G D Q solution is 15 x 15. Also included in Table 1 are the results of Buckle et al. [-5] using a mesh size of 65 × 65 and the results o f X u et al. [11] using a mesh size of 80 x 80. The results of Buckle et al. [5] are obtained by the second-order central difference scheme while the results o f X u et al. [11] are given from the secondorder upwind Q U I C K scheme. It is noted that in the present study, the value of stream function is directly given from the solution of vorticity-stream function formulation, whereas in the work of Buckle et al. [-5] and Xu et al. [11], the value of stream function is computed from the solution of primitive variable formulation. From the table, it can be seen that G D Q results agree well with the second-order difference results even though a very coarse mesh is used. Actually, the m a x i m u m absolute value of stream function (I]/mi n for case A, I//ma x for cases B and C) computed by the G D Q method agrees very well with that computed by the second-order difference scheme. However, there are some deviations between the computed minimum absolute values of stream function (~ma~ for case A, //./mi n for cases B and C). These deviations can be considered to be negligible since the minimum absolute values of stream function are very small.

Table 1 Comparison

of computed

minimum

and maximum

stream function

M e s h size

I/./....

Case A

GDQ B u c k l e et al. [ 5 ] X u e t al. [ 1 1 ]

CaseB

GDQ

15×15

6.8088×10

2

1.1722×10-1

B u c k l e et al. [ 5 ]

65 × 65

-- 5.0203 × 10

2

1.1796 × 1 0 - 1

Xuet

80×80

-4.4332×10

z

GDQ

15 × 15

- 7.4951 × 10 - 3

2.8316 × 101

B u c k l e et al. [ 5 ]

65 x 65

-

2.8437 × 101

X u et al. [ 1 1 ]

80 × 80

- 5.7979 × 10 - 4

Case C

al. [ 1 1 ]

15 × 15 65 × 65 80×80

@min -- 2 . 2 2 3 4 × 10 - t -- 2.3447 × 10 1 2.1724×10 i

1.1936 × 10

3

5.4623 × 10 ~' 1.5642 × 10 o 4.0630×10 0

1.1722x10

2.8409 × 101

1

C. Shu et al. /Journal o f Crystal Growth 181 (1997) 427 436

The G D Q results in Table 1 are obtained by using SOR iterative method. In the present study, the relaxation factors are chosen as 1.35, 1.75, 1.15, 1.0 for S, Q, T, and ~p equation, respectively. The convergence criterion is set as Res =

N/

takes more iterations to meet the convergence criterion. For cases A and B, the residual of vorticity equation is reduced by about 9 orders (from l0 s to 10 4) within 570 iterations. For case C, the initial residual of vorticity is very high (107) due to large temperature gradient from initial guess. Thus, to satisfy Eq. (16), case C requires much more iterations than cases A and B. It should be indicated that the reduction of vorticity residual for case C is actually 11 orders (from 107 to 10 -4) rather than 9 orders for cases A and B. The efficiency of G D Q method in numerical simulation of above three cases is very obvious. Figs. 5 7 show the streamlines and isotherms of G D Q results for cases A, B, and C. For case A, the flow is generated by rotation of the crystal, and there is a primary vortex. For case B, the crystal and crucible are rotated in opposite directions. As a result, there are two vortices with opposite directions appearing in the upper left corner just below the crystal and the lower right corner. Case C is a natural convection problem since the flow is

I N-I M-I (N - 2)(M - 2) ~ ~ ReszJ < 10-4' i=2

j=2

(16) where Res represents residual of S, Y2, T, and equation, respectively. Starting from zero value of initial guess, it was found that the convergence of G D Q solution for all three cases is very fast. To satisfy Eq. (16) for S, t2, T, and ~ equations, cases A and B only need 570 iterations while case C requires 2550 iterations. As a result, the CPU times required on CRAY J916 are 13.2, 13.2, 56.5 s, respectively for case A, B, and C. The convergence histories of G D Q solution for cases A, B, and C are displayed in Figs. 2 4. It can be seen from these figures that for all the cases, the vorticity equation

t .OE+6

433

-7

| 1.0E+5 7 /

~

I.OE+4 ~-, 10E+3 H

1.0E-6

--

1.0E-7

--

1.0E-8

--

1.0E-9

--

x,\

,

~

-. . . .

Temperature

..............

Azimuthal

........

Stream

""

velocity function

-

"-....

"-.i ' '

• ...

I

'

1.0E-10 - 1.0E-11 '

I 1 O0

'

I 200

'

I 300

'

I 400

Iteration Numbers Fig. 2. Convergence histories for case A.

'

500

I

600

434

C. Shu et al. /Journal of Crystal Growth 181 (1997) 42~436 1.0E+6 1.0E+5

t

Vorticity

1.0E+4

~\

1.0E+3 1.0E+2

--

\.\

~

~'-.....,.~ \ \

"~

......

Temperature

..............

Azimuthal velocity

...........

Stream function

1.0E+I 1.0E+O 1.0E-1 1.0E-2 1.0E-3 1.0E-4 1.0E-5

1.0E-6 1.01=-7

1.0E-8 1.0E-9

1.0E-10

'

I 100.00

0.00

I 200.00

I ' ' T300.00 400.00 500.00

I ' L 600.00 700.00

Iteration Numbers Fig. 3. Convergence histories for case B. 1.0E+8

--

1.0E+7

--

~

1.0E+6

I$

1.0E+5

--

1.0E+4

--

1.0E+3

--

1.0E+2

--

1.0E+I

--

1.0E+O

--

1.0E-1

--

1.0E-3

--

I

.OE~

1.0E-5

,

Vo~dty

~

-

~

,"

~

,

,

,

......

Tenr~eral'ure

...........

Stream function

~zx, ,

\F\'

---

1.0E-6

-

1.0E-7

--

-

1.0E-8 -1.0E-9 -1.0E-10

'

I

500

'

1

1000

'

I

'

1500

I

2000

Iteration Numbers Fig. 4. Convergence histories for case C.

'

I

2500

'

I

3000

C. Shu et al. / Journal of Crystal Growth 181 (1997) 42~436

purely induced by the b u o y a n c y force (the crystal and crucible are at rest). There is a p r i m a r y vortex for this case. F o r the t e m p e r a t u r e field, it can be seen from Figs. 5 and 6 that although the streamlines are quite different for cases A and B, the

contours of t e m p e r a t u r e are very similar. This m a y indicate that for forced convection p r o b l e m s (Gr = 0), the t e m p e r a t u r e fields are similar. The contours of t e m p e r a t u r e for case C shows the effect of b u o y a n c y force on the t e m p e r a t u r e field.

S//ll!

a) Streamlines

b) Isotherm patterns Fig. 5. Flow configuration for case A.

a) Streamlines

435

b) Isotherm patterns Fig. 6. Flow configuration for case B.

436

C Shu et al. /Journal of C~stal Growth 181 (1997) 42~436

a)

Streamlines

b) Isotherm patterns

Fig. 7. Flow configuration for case C.

5. Conclusions The application of G D Q method to simulate flows in Czochralski crystal growth demonstrated that GDQ method is a very efficient global method. As compared to low-order methods, it can acquire the same order of accuracy by using much smaller number of grid points, and hence requiring very small computational resources. The implementation of GDQ method is very easy and straightforward. It is expected that GDQ approach can be widely used in future numerical simulations.

References Eli N. Kobayashi, J. Crystal Growth 52 (1978) 425. [2] N. Kobayashi, Comput. Meth. Appl. Mech. Eng. 23 (1980) 21. 1-3] M.J. Crochet, P.J. Wouters, J. Crystal Growth 65 (1983) 153. [4] W.E. Langlois, Ann. Rev. Fluid Mech. 17 (1985) 191. E5] U. Buckle, M. Schafer, J. Crystal Growth 126 (1993) 682. [6] J.R. Ristorcelli, J.L. Lumley, J. Crystal Growth 129 (1993) 249. [7] Q. Xiao, J.J. Derby, J. Crystal Growth 129 (1993) 593.

[8] Q. Xiao, J.J. Derby, J. Crystal Growth 139 (1994) 147. [9] A. Brown, T.A. Kinney, W. Zhou, D.E. Bornside, Proc. 10th Int. Heat Transfer Conf., vol. 1, 1994, pp. 189 203. [10] C.J. Kim, S.T. Ro, Numer. Heat Transfer, Part B 28 (1995) 385. [11] D. Xu, C. Shu, B.C. Khoo, J. Crystal Growth 173 (1997) 123. [12] 1. Raspo, J. Ouazzani, R. Peyret, Int. J. Numer. Meth. Heat Fluid Flow 6 (1996) 31. [13] C. Shu, Generalized Differential-lntegral Quadrature and Application to the Simulation of Incompressible Viscous Flows Including Parallel Computation, PhD Thesis, University of Glasgow, UK, 1991. [14] C. Shu, B.E. Richards, Int. J. Numer. Meth. Fluids 15 (1992) 791. 1-15] R. Bellman, B.G. Kashef, J. Casti, J. Comput. Phys. 10 (1972). EI6] C. Shu, B.E. Richards, Comput. Systems Eng. 3 (1992) 271. E17] C. Shu, Y.T. Chew, B.E. Richards, Int. J. Numer. Meth. Fluids 21 (1995) 723. [18] C. Shu, B.C. Khoo, Y.T. Chew, K.S. Yeo, Comput. Meth. Appl. Mech. Eng. 135 (1996) 229. 1-19] H. Du, M.K. Lim, R.M. Lin, Int. J. Numer. Meth. Eng. 37 (1994) 1881. [20] H. Du, M.K. Lim, R.M. Lin, J. Sound Vib. 181 (1995) 279. [21] H. Du, K.M. Liew, M.K. Lim, J. Eng. Mech. 122 (1996) 95. [22] C. Shu, Int. J. Mech. Sci. 38 (1996) 935. [23] C. Shu, J. Sound Vib. 194 (1996) 587. [24] A.A. Wheeler, J. Crystal Growth 102 (1990) 691.