Digital simulation of Czochralski bulk flow in a parameter range appropriate for liquid semiconductors

Journal of Crystal Growth 42 (1977) 386—399 © North-Holland Publishing Company

DIGITAL SIMULATION OF CZOCHRALSKI BULK FLOW IN A PARAMETER RANGE APPROPRIATE FOR LIQUID SEMICONDUCTORS W.E. LANGLOIS IBM Research Laboratory, San Jose, California 95193, USA

Six numerical simulations of time-dependent flow in a crystal-growth crucible are described. The geometrical and physical parameters are in realistic ranges for growth from liquid semiconductors. Except when crucible rotation is absent, a well-defined Taylor column forms under the crystal face. Outside of the Taylor column, the meridional flow oscillates indefinitely between two qualitatively different configurations. The Taylor column itself is vertically stratified into four distinct regions.

1. Introduction

Clearly this situation is unsatisfactory. Quite apart from computer expense, the long turnaround time makes it difficult to carry out the parameter studies which are potentially a most useful application of digital simulation. Moreover, as work progresses the model will become more sophisticated. Physically, we may want to incorporate the variation of viscosity with temperature (which is relatively easy) and radiative heat transfer within the melt (bad!). Geometrically, numerical studies thus far have been limited to what might be called Czochralski bulk flow: the top boundary is assumed flat. In practice, the face of the growing crystal may be concave or convex, and is somewhat elevated above the nominal level of the free surface. This discrepancy may not be too important as far as the flow in the bulk of the melt is concerned. However, the properties of the growing crystal are presumably influenced most strongly by the flow in its immediate vicinity. This complicated region bounded by the crystal face, the free-surface meniscus, and the nominal level of the melt (where the local flow couples to the bulk flow) will presumably be the subject of detailed attention at some future time. To accommodate these generalized models, more efficient numerical schemes must be evolved. Some limited improvement can undoubtedly be achieved with better computational algorithms, but a more promising approach is the use of nonuniform grids. If the mesh-points can be densely clustered in the regions of steep gradients, and sparsely distributed

A procedure for numerical solution of the timedependent equations governing the flow in a Czochralski crystal-growth crucible was described and demonstrated in a recent article [1]. For one of the cases treated, the steady-state flow had previously been computed by Kobayashi and Arizumi [2]. When the time-dependent computation of [1] settled down to a steady state, the streamlines of the rather cornplicated flow in the meridional planes matched quite well with those ifiustrated in [2]. This lends credibility to both techniques. Nevertheless, [1J must be regarded as a pilot study because all the cases cornputed assumed an artifically high viscosity (or, equivelently, an unrealistically small crucible), at least as far as application to liquid semiconductor Czochraiski growth is concerned. When the viscosity is lower, steep gradients of the flow variables occur within the fluid. To resolve these layers of rapid change requires a fine-mesh grid, and hence long computation time. If the simulation timestep is limited only by considerations of truncation error, the run time is approximately proportional to the number of grid-cells. More commonly, however, refming the spatial grid requires shortening the timestep in order to avoid computational instability. The run time is then approximately proportional to the 3/2 power of the number of cells in a two-dimensional grid. The simulations described in the present paper each require more than 25 mm of computation on an IBM 360/195 per simulated second.

—

—

386

WE. Langlois / Digital simulation of Czochralski bulk flow

387

elsewhere, the required resolution can be achieved at a considerable saving in computer time. However, as Schulz-Dubois [3] has pointed out, “The phenomena of hydrodynamic flow in rotating systems are complex and liable to defeat intuitive prediction even for the experienced specialist.” Carruthers and Nassau [4] have illustrated several intricate possibilities afforded by Czochralski bulk flow. Thus it seems prudent to defer the use of nonuniform grids until we have a semi-quantitative understanding of the locations and sizes of the steep-gradient regions. The present paper is intended as a first step toward that understanding. Six related cases were simulated

are set out in the appendix. The numerical procedures are those employed in [1], where they are described in detail, except that “upstream differencing” is used instead of centered differencing. The lower viscosity necessitates this modification.

on a 50 X 150 uniform grid, with physical and geometrical parameters in the appropriate ranges for Czochralski growth of semiconductorcrystals. Section 2 describes what might be termed the basic isothermal case of this article. The intent was to repeat, with realistic parameters this time, the approac o [1]. That is, to continue the simulation until the flow field evolved from an initial guess to a steady state. never did. after the transients died out, the It secondary flowLong in the meridional planes slowly oscillated between two qualitatively different configurations. Section 3 illustrates the strikingly different behavior of rotating and non-rotating fluids. The simulation is carried out with the same parameters except that the crucible does not rotate. This case does approach an asymptotic steady state which is completely different from every phase of section 2’s basic case. The effect of changing the viscosity coefficient is examined via the two cases presented in section 4. Section 5 describes a simulation which differs from the basic case only in that the crucible rotation rate is reduced ten percent. Section 6 adds heat transfer to the basic case. The neglect of heat transfer in five of the six cases is in no way intended to suggest that buoyancy is relatively unimportant compared with the centrifugal pumping. The actual reason is far more prosaic: at the time this series of numerical experiments was begun, the portion of the code which handles heat transfer along the upper boundary wasn’t working properly. Hence, priority was given to isothermal cases until the problem was corrected. The governing equations and boundary conditions

cible radius Rc, melt depth H, and crystal rotation rate &2~.Specifically, R5 = 4.18, Rc = 6.40, H~3.94, Il~= —2.31

2. An isothermal simulation

In the notation of the appendix, the parameters which were kept the same for all six cases discussed in this paper are the normalized crystal radius R5, cru-

A 50 X 150 network of grid points was used, so that —

—

—

—

The present section describes a simulation the 2c and thewith dimennormalized crucible rotation rate ~ sionless viscosity E specified as ~ 1 E c 0.01 Heat transfer was neglected i.e. the dimensionless buoyancy B and conductivity K, and the surface emissivity were all set to zero. Computation was begun from an initial guess for the azimuthal velocity, stream-function and vorticity fields. The original plan was to proceed as in [1], i.e., to continue the simulation until the flow settled down to a steady state. After 35 simulated seconds, the azimuthal velocity field was as shown in fig. 1. This exhibits a rather pronounced “Taylor column” effect the influence of the crystal face is felt all the way to the bottom of the crucible. In the annular region under the free surface, the flow is very nearly a rigid body rotation. Within the Taylor column, however, there are four distinct regions, stratified vertically. Directly under the crystal there is a region of counter-rotation in which the fluid turns with the crystal rather than with the crucible. Below this, there is a transition region where the flow changes to iso-rotation. Most of the Taylor column consists of a region where the azimuthal velocity is essentially independent of depth. This is not, however, a region of solid body —

—

—

—

—

W.E. Langlois / Digital simulation of Czochralski bulk flow

388

0.00

0.88

1.75

2.63

3.50

4.38

5.25

6.13 ~

°0.00

0.88

1.75

2.63

3.50

4.38

5.25

6.13

____-

°

Fig. 1. Azimuthal velocity field for the basic case described in section 2.

rotation: the angular velocity of the fluid, which is considerably slower than that of the crucible, depends on the distance from the axis. Finally, there is another transition region at the bottom of the column a boundary layer across which the angular velocity of the fluid adjusts to that of the crucible. Note that in both transition regions the velocity profile “overshoots”: within these regions, the azimuthal velocity is not a monotone function of depth. The secondary flow, i.e., the circulation in the meridional planes, behaved unexpectedly. After more than 40 simulated seconds, even its qualitative appearance failed to reach a steady state. The ultimate time-variation within the Taylor column is not great and could have been attributed to numerical noise. Outside the column, however, the secondary flow oscillates slowly between two strikingly different patterns, as ifiustrated in the squence of fig. 2. In fig. 2a, there is a pronounced back-eddy, which is considerably weakened in fig. 2b. In fIgs. 2c and 2d, it has disappeared completely; the circulation everywhere has the same sense as that of the Taylor —

column. The back-eddy then begins to reform. In fig. 2e, the secondary flow is almost stagnant outside the Taylor column, and in fig. 2f the back-eddy has nearly regained its full strength of fig. 2a. This sequence was the sixth repetition of a cycle which had a period of about 2 1/3 simulated time units. Although visually dramatic, the appearance and disappearance of the back-eddy is a relatively minor aspect of the Czochralski bulk flow. Even within the Taylor column, the magnitudes of the secondary flow velocity components are only about 10% of the typical azimuthal speeds. Outside the column they are weaker still less than 2% of the azimuthal flow. —

3. Stationary crucible If the crucible is not rotated, the situation is cornpletely different. A simulation with all parameters the same as in section 2, except with ~ = 0, evolved uneventfully to a steady state. Since fluid rotation is confined mainly to a boundary layer near the crystal

~—

C.

/1 )

0.00

2.19

3.06

3.94

9’

9’

\~

I

I

I~II’/”’~

IIii7 UI)’

0.00

0.00

9’

(TI

~ ~

~

0

~ (Si

(Ti

~

~

0

.,) (Ti

r.i

I)

\U

~ t ~

9’

1.31

/

/

/ ~~/\U

394

r...

2.63

9’

~--

~—~~—————-~

~

~

I

~

1.31

F.J (Si

9’

~

~

~

0

ui

~

;0

0.00

~_--~----

—---~

1.31

~

-~—

--

1.31

2.63

2.19

___

_j_~__~__~

-

z—*

\

3.06

//

~

‘~

/

/

/

~

(Si

91

3.94

/~

-

~II~I#H(“ ~

~

//‘~iIfl / I HI

7~it\~

7!

3.94

00

‘~‘~

I

/

~

I

~

0)

~

~

~---~

—

N.

-.—~

a

0.00

I

N

N

~ I~~~--.:

11/

i///

II ~ II1/ Il (

(71

~

0

0

b

0

0.00

~—--~

1.31

—

~_--—--—~--.-—-—

1.31

--

- --

--

-

-~

2.19

I

-

~

co

/

\t1~iI’

‘

/

1

c~)

9)

9’

co

°

3.94

I

i~i ~

‘I

/~~U!1~i.~ ~ ~

—~-

3.06

——----—--

--

0

~

9’

/_~

~

/

~-

~----~-

0.00

,

,~

7/

N N ~-

1.31

~-

I~

-~

1.31

I\-~~N~

~

~

1/7

(/~~

Ii II 7

/

I/

I

/

0.00

co

~ ~

~ ~

~ o’

p

b

0 0

0

3.94 0 b

/ /~I\Li~

\i

/~~‘(J1

/

---.1

.~:iiiiiI_~~_ ~ -.———

—-

--

2.63

—

-~

-

-

-

~

-

—

-

- -

----------

2.19

-

—

-

--

2.63

--—-------~/‘\iU1)~l

z

~\II~hi(tII ~

3.06

/

//

-.-—7--~~

~

i

,,,,l~Ii(

‘Iii

I Ii

\iiI~

co

OD

p

0

/1

3.94

9’

~

9’

co

--4~~~

-

——_~

—

~

7\IiL~°’

/,-

/

1

I

3.94 ~1p

W.E. Langlois / Digital simulation of Czochralski bulk flow

0.00

0.88

1.75.

2.63

3.50

391

4.38

5.25

6.13

~

c c c c c ççc~-.~ 2 /iHi /i/1’I I I

II

I

~-~.—/

I j II I i~’)i’~\ ft I I

I)

I

/

I

/

I I I I I I I

N

I

(e)

\~ \ \

II

~

//I~\l\~\

~

~

___ ________

Q 0

7/

~I

I I H ~/ I

\

I

0

I I

I

3

I

~

/

0

q I

I

0.00

0.88

1.75

2.63

3.50

—

4.38

5.25

6.13

0.00

0.88

1.75

2.63

3.50

4.38

5.25

6.13

5.25

6.13

/

/

I

/ I I I I \, I

1

/ //

~---~

‘\

\

i//Hr / I I I /1! II /

\

//

N~/

//~

H1

I

tI\\

0.00

0.88

1.75

2.63

Fig. 2. Meridional circulation for the basic case: (a) t

=

~-

/

~

c 7/ c ~ c c c c ~

0

3.50

35.0; (b) t

=

35.4;

4.38 (c)

~

t = 35.8; (d) t = 36.2; (e)

r = 36.6; (1) t

37.0.

~

I

ii

I

~

1<

0) C’)

0.00

(ii

C.

:..~j

(~

0 0

0.00

~

a

0

~

8~

~ I

1.31

I

1.31 I

~

I

2.19

z—* I

2.63

3.06

I

I

0.00

~/

~/

-.---~

___

3.94

3.94

0) C’)

(71

.~

0) C’)

01

1.31

(

2.19

7~N

3.06

3.94

0)

~

a,

.~

0

I

P

-

2.63

0

/1

z—* 0 0

/

1.31

0 0

-~

0.00

0 0

3.94.

2

t’J

W.E. Langlois / Digital simulation of Czochralski bulk flow

0.00

0.88

1.75

II~

2.63

3.50

~

4.38

~

I

393

5.25 I

I

I

6.13

~-

I

0)

~TIII~IIH 0.00

0.88

1.75

2.63

3.50

4.38

5.25

6.13

5.25

6.13

I

I

r Fig. 4. Meridional ciculation for the high viscosity case.

0.00

0.88

1.75 I

2.63 ____

3.50

4.38 1

__________

~.

9)

~~~/7II

I

I

0.00

I

0.88

I

C

I

1.75

2.63

3.50

4.38

r Fig. 5. Meridional circulation for the low viscosity case.

0

I

5.25

6.13

W.E. Langlois / Digital simulation of Czochralski bulk flow

394

face, there is no Taylor column. The steady-state meridional flow consists instead of a “tip vortex” located just off the edge of the rotating crystal face, as illustrated in fig. 3a. This vortex transports some of

section are all in the regime for which there is a boundary layer on the crystal face, it is of interest to see if the thicknesses correlate. Fig. 6 is a plot for each of the three cases of L/v’l2E/I~25lversus r,

the angular momentum imparted by the rotating crystal out of the boundary layer into the main body of the melt. The isotachs of fig. 3b are logarithmically spaced to manifest this effect,

where L(r) is the dimensionless vertical distance from the crystal face to the dividing streamline of the meridional flow. This dividing streamline corresponds quite closely to the boundary between the counterrotation and iso-rotation regions. From this figure it appears that ~J12v/I (~)~I is a representative estimate of the boundary layer thickness for each of the three simulations. This cannot, of course, be taken as a general conclusion on the basis of so few runs with no variation at all of &2~,nor of other factors (such as crucible rotation rate, melt aspect ratio, thermal effects) which might modify the flow pattern near the crystal face.

4. Different viscosity Two simulations were carried out with all parameters as in Section 2 except for the coefficient of viscosity. In one case it was increased toE = 0.015; in the other it was lowered to E = 0.003. Qualitatively, both runs were similar to the basic case with E = 0.01, including the ephemeral back-eddy. The meridional flow for the high viscosity case is shown in fig. 4 and for the low viscosity case in fig. 5. Analytical and semi-analytical results [5] suggest that the thickness of the_boundary layer on a rotating disc is proportional to ~Jv/I w~ I, the square root of the kinematic viscosity divided by the rotation rate. Since the flows described in section 2 and the present

1. SO

1.25

E

=

E

0.03

=

—

5. Slower crucible rotation The next simulation used the parameters of sections 2 except for the crucible rotation rate, which was reduced ten percent i.e., &~= 1.4 13. The results were again qualitatively similar to those of sections 2 and 4. Fig. 7 illustrates a late stage of the flow. Fig. 8 indicates that the boundary layer on the crystal face is somewhat thicker than would be expected from the \/1 2z’/I w~estimate.

0.10

6. Thermal effects L

—

1.00

12 E/I~25I

ff

E

=

\/

0,15

In the final simulation, buoyant convection and heat transfer were included by taking

0.75

B=0.0l382,

K=0.l25, X i0—°~.

3rea/Xpc~= 7.94

O

The crucible and crystal-face temperatures were specified, respectively, as

0.25

T~= 1773.0, 0.00 I

0

I

1

I

2

r

3

T~= 1685.0.

I

‘

Fig. 6. Comparison of the boundary layer thicknesses for, three values of viscosity. The function L(r) is defined such that the dividing streamline between the two vortices of the Taylor column is the locus z = H — Lfr).

All other parameters were as in section 2. Since the buoyant force is distributed throughout the melt, the Taylor column is somewhat less clearly delineated than in the isothermal cases. In particular, it tends to mushroom outward at its base, as can be

W.E. Langlois / Digital simulation of Czochralski bulk flow

0.00

0.88

1.75

II

tN

cc~ C’, C\J 0

I’

0.00

I

I

I

2.63

I

I

1.75

I

2.63

Fig. 7. Meridional

1.413

7 L 12

1.00

EfIs~I

2~’\

-~

~~-.-J~-~ /~

~2

-

1.57

‘~

j

0.75

0.50

Fig. 8. Comparison of the boundary layer thicknesses for two

values of the crucible rotation rate,

5.25

I

I

3.50

I

I

I

4.38

6.13

I

I

~-

0)

I

I

I

5.25

6.13

—

c’i 0 0) 0

circulation fox the reduced rotation case.

1.50

1.

4.38 I

:I]~

I

0.88

3.50

~

395

seen by comparing fig. 9 with figs. 1 and 2. The final temperature distriubution, shown in fig. 10, exhibits several interesting features. First of all there is the small region along the free surface, just beyond the crystal boundary where the temperature is below T~.This is a region of exceptionally feeble convective heat transfer. Also, the model has no phase-change mechanism to warm the region by release of latent heat. Consequently, the heat lost by radiation from the free surface is replaced almost exclusively by the relatively weak mechanism of conduction. The net heat transfer comes to equilibrium at a temperature slightly below T~.In practice, the —

—

melt may instead freeze, leading to the inconvenient phenomenon called “flaring”. Next, there is a rather large region of relatively cool melt beneath the crystal face. The region of counter-rotation provides a layer of “insulation”, which blocks convective heat transfer from the main body of the melt.

WE. Langlois / Digital simulation of Czochralski bulk flow

396

0.00

0.88

I

1.75

I

I

I

2.63

3.50

I

I

4.38 I

I

5.25 I

6.13

I

I

~. 0)

0

c-.j 0)

‘1~ N

C’)

(a) 0

P

____

0.00

0.88

I

I

1.75 I

_____________

2.63

r 0.00

I

0.88

1.75

~

3.50

4.38

3.50

4.38

I

5.25

6.13

I

—#

2.63 ___

5.25 I

I

6.13 I\

I

~

0) :

0

I

(0

0)

N

~

I)

(b) 0 0 I

0.00

I

I

0.88

1.75

2.63

3.50

r

4.38

5.25

6.13

-+

Fig. 9. Flow field for the non-isothermal case: (a) meridonial circulation; (b) azimuthal velocity field.

W.E. Langlois I Digital simulation of Czochralski bulk flow

0.00

0.88

1.75

2.63

3.50

C)

4.38

397 5.25

6.13 ~ C) ~vj

___

CO

0 c’~i C’,

t N

C,.~

0 0

—~

__________

I

—

0.00

~

I

I

0.88

1.75

I

I

0 0

709

________j I

I

2.63

3.50 r

I

I

4.38

I

I

5.25

I

I

6.13

—~

Fig. 10. The temperature distribution.

Finally, along the outer boundary of the Taylor column, the isotherms show a tendency to bulge upward as a result of the pronounced upwelling in this part of the flow field. Thermal effects cause the boundary layer on the

1.50

Isothermal

crystal face to become slightly thinner, as illustrated in fig. 11. The thickness of the transition region at defithe bottom of the crucible is also ofinterest, but its

,7-___.__~_~~,_/

I

am-

Heat transfer and buoyant convection included

nition is more subjective. Near the top of the Taylor column, there is a dividing streamline between the two principal eddies of the meridional flow, and its distance from the crystal face gives a convenient definition of the boundary layer thickness. Near the

a. so0.

2&

0.00 0

I

I

I

I

I

1

2

3

4

S

Fig. 11. Comparison of the boundary layer thicknesses for isothermal and nonisothermal flow.

in the same sense. Forcolumn, each of fiveis everywhere simulations bottom of the Taylor thethe flow with rotation, however, it appears meaningful to defme the thickness of the lower transition region as \/12v/I Wc I.

W.E. Langlois / Digital simulation of Czochralski bulk flow

398

7. Conclusions and program

Appendix

It would be inadvisable to draw far-reaching conclusions on the sole basis of the six simulations described above. (Actually only five are relevant because section 2 was included only to contrast rotating and non-rotating fluids.) The high-viscosity simulations described in [I], for example, bear only a superficial resemblance to those of the present paper. Also, one might anticipate a more-or-less gradual transition out of the regime of section 5 into that of section 3 as the crucible rotation rate is reduced further. Nevertheless, the results presented here do exhibit some trends which seem to persist as we move about the region of parameter space appropriate for Czochralski growth of semiconductor crystals. The well-defined Taylor column, with its four vertically-stratified regions, was present in each of the five simulations having ~‘2~ 0, as was the oscillatory character of the flow at large time. The thickness of the boundary layer on a crystal face or crucible bottom rotating with angular velocity

Equations ofmotion and boundary conditions

1

As in ref. [1], we normalize the dependent and independent variables with respect to the physical units of time, length, and temperature, denoted by r, X, and 0, respectively. If the measuring system is CGS (r = I sec, X = 1 cm, 0 = 1 K) then the physical and geometrical parameters used in this paper are representative of those which might well obtain in Czochralski growth from a liquid semiconductor. The radial and axial velocity components, u and w, respectively, are derived from a Stokes streamfunction according to i~i

u=(l/r)a~i/az

w=—(l/r)a~i/ar.

Defining a vorticity w = aw/ar nostic for i,L’: 2~~ a ~l a~\+ 1 a —i- = —w. r r r T az

-s-)

~—~--

—

au/az yields a diag-

—

W 0 was, in each case, not far from ‘IV/l2P/kooI. This estimate, however, does not reflect other influences, such as the rotation rate of the opposite surface or thermal effects. Figs. 8 and 11 show definitively that these influences matter, but they do not provide enough information for quantitative modification of the estimate. More simulations of theofsort reportedat here indicated. The persistence oscillation largeseem simulation time, as well as possible interest in accelerated crucible rotation techniques [6], suggest that timedependent simulation is mandatory for the parameter range of liquid semiconductor Czochralski growth. In other ranges, of course, steady-state techniques may be quite useful. Improving the computational efficiency should be given high priority. The pronounced vertical stratification of the flow can be used to immediate advantage. A straightforward transformation, involving only the vertical coordinate, would allow one to concentrate the grid points in the steep-gradient regions at the top and bottom of the flow field. Roache [7] has pomted out that coordmate transformation is, m general, the preferable method of increasing resolution locally,

Prognostic equations for the vorticity, the azimuthal velocity v, and the absolute temperature T are provided by +

~

3r

aT +

—

~

E r

a (uu) +

+

~

r az

B _T ar + EV

a

av

(ww) +

-~-~

az 2w

a

=

at

(uw) +

~-

at

2uv (wu)

+—

u =

EV2v

—

E

i a a — 2 r ar (i~T)+ az (wT) — KV T,

where ~2

=

a2/ar2 + (l/r)a/ar + a2/az2

the dimensionless coefficients are defined in terms of the dimensional physical parameters of the melt according to 2

B—ogr 0/A,

E=vr/X

2

2

,

K=~r/X

wnere a is the volumetric expansion coefficient, g is the gravitational acceleration, v is the kinematic

W.E. Langlois / Digital simulation of Czochralski bulk flow

viscosity, and ,~is the thermal diffusivity. The region of computation is the rectangle 0 ~ r s~ R~,0~z~H,whichiscoveredbyanMXNnetwork of evenly spaced grid points. The grid spacing is ~r X z~.z,where =

RC/(M

—

1)

and

~.Z

=

H/(N

—

1).

The crucible rotation rate is w~ &2c/T and the crystal rotation rateis w~ c15/r. The boundary conditions are those used in [1]: axial symmetry; no-slip velocity condition and specified temperature all solid surfaces; Boltzmann radiationon from the free surface,Stefan— which remains atregs-free and almost flat. As shown in [1], this leads to the following: Axis, r = 0 o~=0, w=0,

u’0,

Crucible bottom, z =

0,

w

=

—

=

aT/ar=o;

v=r~’l~, T=T5 Free surface, z = H, R5
o,

aT/at =

w = 0,

au/az

=

0,

ar/at

from the prognostic equation 3reoT4IXpcvA.z

—2O where e is the surface emissivity and the dimensional parameters p, ci,, and a are, respectively, the melt density, the melt specific heat, and the Stefan—

correspond to grid cells of height ~z/2 rather than References

t),

~,

Crystal face, z = H, r ~ R~ 2 0, W=_(~)2 o,(i(r,H— L~.z,t),

Boltzmann constant. The factor 2 arises in the radiation term since the grid-points lying in the surface

0

r(~z)2~(r,

i,D

399

a = rc2~, T = T~

[11 W.E. Langlois and C.C. Shir, Computer Methods in Applied Mechanics and Engineering 12 (1977) 143.

Crucible wall, r = Rc

[2] N. Kobayashi and T. Arizumi, J. Crystal Growth 30 (1975) 177. [3] E.O. Schulz-Dubois, J. Crystal Growth 12(1972)81. [41J.R. Carruthers and K. Nassau, J. Appi. Phys. 39(1968) 5205. [51H.P. Greenspan, The Theory of Rotating Fluids (Cambridge Univ. Press, 1968). [6] H.J. Scheel, J. Crystal Growth 13/14 (1972) 560. [7] P.J. Roache, Computational Fluid Dynamics (Hermosa Publ., 1972).

=

0

W

a = Rc~lc,

2 R(~r)2I/I(R~ ~.r,z, t) —

T Tc

‘