Conservative differencing procedures for rotationally symmetric flow with swirl

COMPUTER METHODS IN APPLIED MECHANICS NORTH-HOLLAND PUBLISHING COMPANY

AND ENGINEERING

25 (1981) 315-333

CONSERVATIVE DIFFERENCING PROCEDURES FOR ROTATIONALLY SYMMETRIC FLOW WITH SWIRL W.E. LANGLOIS IBM Research Laboratory,

San Jose, California

Manuscript received 28 December

95193, U.S.A.

1979

When swirling motion is accompanied by a rapid radial component of flow toward the axis, numerical procedures for handling the coupling term in the azimuthal velocity equation become unstable. The difficulty is averted by replacing this equation with one which governs the angular momentum per unit mass. Moreover, since this is a physically conserved quantity, differencing schemes can be constructed directly from a balance law. Conservative schemes thus emerge automatically. In plane motion of an incompressible fluid, the vorticity is a pseudo-conserved quantity, i.e., it is governed by a differential equation which resembles one derived from a physical balance law. This is not the case for the vorticity of the meridional motion of a rotationally symmetric flow. However, a little known theorem due to the 19th century hydrodynamicist A.F. Svanberg demonstrates that the vorticity divided by the distance from the axis of symmetry is pseudo-conserved. This is useful for deriving differencing schemes governing the meridional motion. In the present paper these concepts are used to obtain prognostic procedures which are then applied to a problem of flow in a crystal-growing crucible. Crystal and crucible rotation provide the swirl and thermocapillary flow provides the inward radial motion which renders earlier procedures unsatisfactory.

1. Introduction

A series of recent papers [l-5] describes the results of numerical simulations on the fluid flow in a Czochralski crystal-growing crucible. As illustrated in fig. 1, the growing crystal is extracted slowly from the melt. To disperse inhomogeneities, the crystal is rotated as it is lifted; in some applications, the heated crucible is also rotated. Thermal convection in the differentially heated melt is treated with the Boussinesq approximation, i.e., fluid density changes are ignored except in the buoyancy term. In practice, the free surface is not really flat: there is a meniscus near the growing crystal and there may be ripples on the melt surface; these have been ignored (the Czochrulski bulk flow idealization). Also, the vertical motion induced by lifting the crystal is neglected, as it is typically three orders of magnitude smaller than the motion induced by the rotating crystal. For investigations that focus on the details of flow near the crystal [Ml, the vertical motion can be accounted for by treating the crystal face as a solid surface with suction [S]. The obvious mechanisms for generating a meridional flow in the melt are centrifugal pumping and buoyant convection. It has recently been demonstrated [9] that a third mechanism, thermocapillary flow, can also be significant because of the high thermal gradients, and concomitantly high surface-tension gradients, that obtain in Czochralski growth. The

316

W.E. Langlois, Conservative differencing procedures for rotationally symmem’c flow with swirl

Fig. 1. Geometry of the Czochralski technique.

importance of thermocapillary flow becomes overwhelming for crystal growth in space vehicles [5, lo] but some importance remains even under terrestrial conditions. Thermocapillary flow is accentuated by crucible rotation: the Taylor-Proudman effect gives it muscle. The primary intent of the present work was to study the coupling of thermocapillarity and rotation in a system of technological importance, viz., Czochralski growth of elemental silicon. In the course of the investigation, however, there arose a computational consideration which is of some interest in itself, and which is the main focus of this paper. References [l] through [5J successfully updated the azimuthal velocity v by using a scheme that was derived from the prognostic equation 2uv ~+$(uv)+-&v)+I=

vv2v - v ;,

(1)

where V2 signifies a’/&* + (l/r)a/& + a2/8z2 and v = g/p is the kinematic viscosity. In the present study that scheme fails. Near the top of the melt, thermocapillarity induces a strong inward component of radial flow. Eq. (1) then becomes dominated by the coupling term, i.e.,

*A!lv at

( > r

which is the equation

’

for a growing exponential.

In the earlier papers, its finite-difference

W.E. Langlois, Conservative differencing procedures for rotationally symmetric flow with swirl

317

approximation was unconditionally unstable, but the instability was damped by the action of the other terms in (1). The difficulty appears to be related to the fact that azimuthal velocity is not a physically conserved quantity. Eq. (1) is therefore somewhat removed from the basic balance laws. The called conserved variable is the product r~ = 0, the angular momentum per unit mass-often the swirl. Rearranging (1) yields a prognostic equation for this variable:

(2) As we shall see in section 2, the terms in this equation relate directly to the angular momentum changes produced in a grid cell by advection and viscous torque. The azimuthal velocity vanishes at the axis of symmetry, and hence LJ vanishes quite strongly. Nevertheless, the cylindrical grid cells lying along the axis may, because of their finite volume, contain a small residue of angular momentum. The implications of this are discussed in section 3. In view of the instability encountered with (1) it seems appropriate to take a closer look at the prognostic equation for the vorticity before it too starts giving trouble. The situation is somewhat different here because vorticity is a derived, rather than fundamental, variable. The primitive equations of viscous flow are obtained by starting with physical balance laws, which stipulate that mass, momentum, and energy are conserved. Gauss’ divergence theorem is invoked to convert these laws to partial differential equations, with various constitutive relationships and the constraint of incompressibility supplying certain key steps. Derivation of a finite-difference scheme in essence reverses the process: partial derivatives are replaced by their finite analogs, which leads to prognostic equations for the grid-cell averages of the dependent variables. In the case of fundamental variables, the differential equations can, if one desires, be by-passed completely: the finite-difference scheme can be derived by applying the integral balance laws directly to the grid cells. In fact, this is an excellent way to generate “conservative” differencing schemes. The vorticity transport equation is not one of the primitive equations, but rather is derived from them by cross-differentiating to eliminate the pressure, with well known computational advantages. It has been pointed out [ll (p. 421) 12 (p. 6811that the properties of this equation make it seem natural to regard vorticity as something carried about by the flow field, rather than as the purely kinematical quantity which it actually is. The image is especially strong in two-dimensional flow, for which the vorticity equation takes the exact form of the heattransfer equation. A conservative differencing scheme for the vorticity can then be obtained by re-tracing the analysis which led to a conservative scheme for the thermal energy. For rotationally symmetric flow, the analogy between vorticity-transfer and heat-transfer no longer holds. Finite-difference analogs of the vorticity equation are consequently obtained without any guidance from underlying conservation principles. However, in his numerical treatment of a flow with no swirl, Kyrazis [13] found that certain errors could be avoided by using o/r, rather than o, as a dependent variable. Other authors [14, 15] had previously used w/r, with varying degrees of success. A physical reason for the effectiveness of w/r as dependent variable is suggested in [ll (p. 421)], which completes the proof of a very old result known as Svanberg’s vorticity theorem.

318

W.E. Langlois, Conservative differencing procedures for rotationally svmmetric flow with swirl

For an incompressible fluid, this theorem states that a rotationally symmetric motion is circulation-preserving if and only if o/r is constant for each fluid particle. (It is the rotationally symmetric analog of D’Alembert’s vorticity theorem which establishes w as the corresponding key quantity in plane flow.) Thus, in a circulation-preserving rotationally symmetric flow, w/r behaves as if it measured some fundamental property of the fluid particles, moving about with them as they are advected through the flow field. Just as conservative differencing indirectly becomes meaningful for w in plane flow, so does it also for o/r in rotationally symmetric flow. The theorem also suggests a name for w/r. We shall term it the Svanberg vorticity and denote it by the symbol S. In section 4 we derive finite-difference equations governing S away from the axis of symmetry. The corresponding analysis for control volumes along the axis is developed in section 5. A somewhat surprising result emerges: when viscosity is present, the axis of symmetry acts as a source of S. This serves to remind us that the notion of S as a physical property of the fluid particles is a limited analogy. At bottom, the Svanberg vorticity remains merely a convenient kinematical variable. Section 6 presents a brief exposition of the results obtained for the effect of thermocapillary flow on Czochralski bulk flow in silicon melt. It should be acknowledged that the earliest work on digital simulation of the flow in a Czochralski crucible was carried out by N. Kobayashi and T. Arizumi [16-H], even though the present paper does not refer to any of their specific results.

2. Conservative

differencing for the angular momentum

- off the axis

Those grid cells which do not lie on the axis of symmetry are annular rings of rectangular section; a slice of a typical cell C is illustrated in fig. 2. Within the meridional section of C is the grid point (pi, rj), with r = (i - 1)Ar. The boundary 8C of C consists of the circular cylinders

and the plane annuli 2, = {(r, 2): ri- < r < ri+,

22).

2 =

4

42 &/%

Ri_

A\

,,,.’

-‘*

Crirzj)

‘\

‘,

\

\

\

\

Fig. 2. OR-axis grid cell.

W.E. Langlois, Conservative differencing procedures for rotationally symmetric flow with swirl

319

In these descriptions rik=rikAr/2=(i-l+)Ar, z- = zi - AZ-/~,

z+ = zi + A2+/2. The upward and downward z-increments may differ, because the study of Czochralski bulk flow is facilitated by using an axially stretched grid [3]. The height of the grid cell is AZj = (AZ- + Az+)/~

and its volume is V = TAZj(rf+ - rf_) = 2miArAzj.

(3)

We shall be concerned with the rate of change of angular momentum, and later of Svanberg vorticity, within C. Rather than repeat’mathematical displays, we shall define certain averages for a generic variable f, which may represent 0, S, or any other quantity whose averages are of interest. We use fi,j to denote the average off over the grid cell, i.e.,

At sufficiently high resolution, fi,j may be regarded as an approximation to the value off at the grid-point (ri, zj), but it is seldom necessary to invoke this interpretation. We also define averages over the bounding surface of C:

(f>z*= &

1:’ r,- rf(r, z,) dr.

(6)

For the purpose of computing advective fluxes, we introduce special weighted averages (f>$, (f)& defined by the equations

Since angular momentum is a physically conserved quantity, AOi,jy the change of fli,j during a time interval At, is determined by the angular momentum flux across 8C during At and by the viscous torque impulse applied to C during this same interval. Thus

320

W.E. Langlois,

The contribution

Conservative

differencing procedures for rotationally

from advection is determined

+

symmetric flow with swirl

by

2miAr((w)z-(a)k- - (w)z+(fl>%+).

The first term corresponds to radial advection (denoted henceforth the second term corresponds to axial advection (~a). With (3)

by the subscript ru) and

As expected, these expressions converge to the advective terms in (2) as the space and time differencing become infinitely fine. In finite-difference approximation, one must make estimates of the various averages. This can be done in many ways, and is the subject of an extensive literature. In the case of the unweighted velocity averages, however, use of the vorticity-streamfunction method to compute the meridional flow imposes a natural choice. If we introduce a Stokes stream-function $, the meridional velocity components are given by u=----,

1w

)$.I

=

r a.7

---*

1w ar

T

From (5) then,

I;-fU)R,=&

[44ri-,

z+)-

W-,

Z-)] A

&

(&,j+l

-

(l/i,j-1

+

+i-lj+l

-

+i-lj-I),

with analogous expressions for ri+(u)R+ and (w)~,. We now determine the torque impulse from viscous forces acting on K’. These arise from radial viscous stress (denoted by the subscript TU)and axial viscous stress (zu). Looking first at the cylindrical surfaces, we recall that the torque acting at Ri, on the fluid within C is

where r& denotes the viscous stress acting in the azimuthal direction upon surfaces oriented normal to the radial direction. Thus

p”@$L_ 2q&(rf+(T,e)tzi+ - rf-(Td)Ri_). For rotationally

symmetric flow, T& is constitutively

related to the azimuthal velocity v and

W.E. Langlois, Conservative differencing procedures for rotationally symmetric flow with swirl

321

its first radial derivative [ 11 (p. B-89)]. In terms of 0 = ru,

With (3) therefore,

The torque acting on the plane annu!ar surface 2, is given by

with

Thus, with (3)

As was the case for the advective terms, there are many ways to estimate the various averages. The grid cells lying just off the axis of symmetry may, however, pose difficulties. For example, if the derivative of O/r’ is estimated by second-order central differencing, then

Since O/r* is the local angular velocity of the fluid, a quantity which may well approach a finite value as r + 0, the second term is an indeterminate form when i = 2. A possible circumvention is to take

NOW, if one chooses to identify &

(-$($))R2_; -3.

as 0(0,

2,)= 0, using

this expression at i = 2 yields

(11)

Thus, if R,,j# 0, the annular grid cell loses angular momentum because of viscous torque at r = Ar/2. With the insistence that 0(0, zj) remain zero, this angular momentum simply

322

W.E. Langlois, Conservative differencing procedures for rotationally symmehic flow with swirl

disappears from the flow field. In other words, with the scheme (11) the cylinder r = Ar/2 acts as an absorbing barrier for viscous diffusion of angular momentum. For flows with intense swirling near the axis, this could be a significant error. There is of course nothing fallacious about product-rule differentiation. The fallacy lies in interpreting C&j as the (vanishing) angular momentum at r = 0 rather than as the average angular momentum in the core grid cell. A simple and accurate alternative to (11) is available. As we shall now show, it is acceptable to take =*

(12)

0.

Rt-

To establish (12), we examine the orders of the various kinematical quantities as r + 0. Since there can be no flow through the axis of symmetry, the radial velocity u vanishes at r = 0. Thus, for small r, u =

u;r + u:r*+O(r3),

(13)

where UL and Ub: may depend upon z and f. The axial velocity w can be finite at the axis. However its first radial derivative must vanish, for otherwise the three-dimensional flow field would have a conical singularity, which viscosity will not permit. Therefore, w = w(J +

wb’r’ + O(r3).

(14)

The coefficients in (13) and (14) are not independent.

The continuity equation

au u aw -$+;+-jg=O mandates that UA = -(1/2)8Wol& u =

U$ = 0. Thus (13) becomes

u;r + 0(r3).

The physical restriction value 0, as r + 0. Thus 0 =

and, more pertinently,

(1% on 0 is that the angular velocity O/r* approaches

AOr + A;r3+ A:r4 + O(r’),

a non-infinite

(16)

and establishing that (12) involves a small error is equivalent to showing that AA = 0.This is accomplished by substituting (14) through (16) into the angular momentum transport equation (2): 3vAL = r 9’+4U&4,+-$(W&,)-8vA;-

v $I+

0(r2).

Since A6 is independent of r, each side must vanish independently. Viscosity prevents a conical singularity in the angular velocity from occurring at the axis of symmetry.

W.E. Langlois, Conservative differencing procedures for rotationally symmetrk pOw with swirl

323

3. Conservative differencing for the angular momentum - the axial core Those grid cells with i = 1 are cylindrical volumes of radius R1+ = Ar/2, centered on the axis of symmetry, as shown in fig. 3. The volume of the cell Co about (0,21) is therefore Vo = ~(Ar)*Azj/4

(17)

and the average value fij of the generic variable f over this cell is given by

VOfl,j

=

The boundary

/fco f du = 2~ LA’* r

]*‘f dz dr = z-

27~ 1” IA”’ rf dr dt.

z- 0

(18)

X0 of Co consists of the cylinder R1+ and the disks

Z”, = {(I-, z): 0 <

r

Averages and flux-weighted

c Ar/2, z = 2,). averages over these disks are defined by

(wf)z: = (w>z$fE:.

(20)

If one accepts (12) as sufficiently accurate, evaluation of 0l.j is a diagnostic rather than prognostic calculation, because (12) is a statement that the angular velocity of the fluid within C is approximately the same as that within the nearest annular neighbor. This permits us to determine 0ij algebraically in terms of L?2j by using the definition of fli,j as the average of R over the grid cell. If we assume that Aor adequately represents 0 for 0 I r 5 r2+ = 3Ar/2, then 2m(Aor2)

dr = YA,,

(21) 2n-r(Aor2) dr =

so that a,,j

A

(22)

Lt*,jllO.

2,” R1+

(O.Zj 1 ,x--

@

--.

\zo-

Fig. 3. Axial grid cell.

324

W.E. Langlois, Conservative differencing procedures for rotationally symmetric flow with swirl

Use of (22) renders the remainder of this section superfluous. However, if one seeks to replace (12) by an approximation of higher order accuracy, then L!,,j must be determined prognostically. The angular momentum change within C,, due to advection across X,, during the time interval dt is determined by

With (17) the radial and axial contributions

are

(23)

Comparing (23) with (8), we see that -(Anl,j)JAt is eight times as large as the contribution to (AC!n2j),JAt from advection across RI+ -= R2-. This is correct for a conservative scheme, since the volume of the cell centered at (Ar, Zj) is eight times that of CO. The factor of eight also enters the calculation of (Afl,,j)JAt, since the torque acting at RI+ = R,_ on the fluid within COis the reaction to the torque which contributes to (AO(t,,j)JAt. We have (24) which brings us back to the problem discussed at the end of the preceding section. For the reasons pointed out there, it is quite acceptable to neglect the right side of (24) entirely. However, continuing this section beyond display (22) assumes that there is potential interest in higher-order schemes. A physically reasonable scheme for evaluating (24) and its counterpart in (9) is obtained by recomputing the integrals in (21), retaining the next non-zero term from (16): 27rr(A,$ + A;r4) dr = (d g

A o + 48 (Ar)4 A”o, (25)

2m(Aor2 + A;r4) dr =

Solving this pair of equations yields A;=------

l6 (J&j - lOi2,,j)* 27(Arr

325

W.E. Langlois, Conservative differencing procedures for rotationally symmetric flow with swirl

Since

=

a n

arf0

2Agr + O(r’),

(24) yields

(26) Finally, computing

the viscous torque on the planar surfaces 2: yields without difficulty

Parenthetical Remark

In principle, considerations similar to those which led to (25) could also be used away from the axis of symmetry, to obtain a differencing scheme which reflects the physics somewhat better than does (10): If we estimate the angular velocity as a linear function of r over the grid cells in question, i.e., if we take

O/r2 = A + A’r,

r(i-l)- I T5 ri+,

(27)

then A’ should replace the right side of (10). As above, A’ would be evaluated by averaging a, as given by (27) over the cells and then eliminating A from the resulting equations. After some algebra, this leads to

A’=&

[l

+ E+(i)] % - [l - E_(k)] +), I

i = 3,4,

(28)

. . . ,

1

where E,(i) =

4(i - 3)f20i2 - (30 & 17)i + 40]+ 240 + 131 4i(i - 3)(80i4 - 480i3 + 1220i2 - 1500i + 1067) + 3091’

Even for i = 3, the discrepancy between (10) and (28) is quite small, and it diminishes rapidly and monotonically as i increases. For i large,

both

E,(i) - 1/4i3.

4. Conservative

differencing for the Svanberg vorticity - off the axis

If the vorticity o is defined in terms of the meridional w =

awlar - &.dl&~,

velocity components

according to (29)

326

W.E. Langlois, Conservative differencing procedures for rotationally symmetric flow with swirl

then w satisfies the equation

~+$(luo)+-$wto)+~~= ,gF+

vv2w-v;;

(30)

the buoyancy term involves new symbols, viz., the volumetric expansion coefficient a, the gravitational acceleration g, and the local temperature T. From (30) one may derive a prognostic equation governing the Svanberg vorticity S = w/r: (31) In section 2, we derived procedures for advancing Q,j by consideration of angular momentum interchange between grid cells. Ultimately one must make estimates of various averages over the bounding surfaces of the cells, but the analysis itself dealt with finite quantities rather than infinitesimals. No analogous procedure applies here, because the Svanberg vorticity is defined differentially, i.e., locally. We must derive the differencing procedure from (31): there is no other route. Fortunately - and this is the advantage of using S rather than w - eq. (31) is of a form which is easily integrated over the grid-cell C. The resemblance between the advective and viscous terms in (31) to their counterparts in (2) suggests in advance that this should be so. Eq. (31) has two distinctive terms, the coupling (8/az)(J2’/r4) and the buoyancy (cuglr)aT/&, but these present no special difficulty. Thus our procedure will be to average (31) over C, obtaining a prognostic equation for advancing Si,j. Extending the notation of section 2 in an obvious fashion, we write

We now integrate (31) term by term, using either of the iterated forms in (18) as appropriate. First,

v

(AS>)ra At --

II

$(nrS)dV=

2?rj=+ [nrS];;dz. L-

With (3), (5) and (7) therefore, @+

=

& I

(ri-(U)R,-(S)&

- rj+(U)Ri+(S);ii+).

Next, (wS) dV = 27r j” r[wS]:; dr. ,iG

(32)

W.E. Langlois, Co~e#a~ue

~i~ere~ing procedures for ~o~tio~l~y symmetric ftow with swirl

327

Eqs. (3) (6) and (7) then yield (33) Applying one or the other of these procedures

to each of the remaining terms in (31) yields

(35) (36) (37) Each of the forms (32) through (37) is conseruatiue. That is, if JC lies entirely within the fluid and off the axis of symmetry, any increase in VSi+jis exactly offset by a decrease in the total Svanberg vorticity of the neighbo~ng grid cells. This applies even to the coupling and buoyancy terms: away from boundaries and the axis, these terms redistribute Svanberg vorticity rather than create it.

5. Conservative differencing for the Svanberg vorticity -the axial core Using (18) to integrate (31) over the cylindrical volume C, forces us to deal with indete~inate forms, since one of the integrations begins at r = 0. The orders of some relevant quantities as r + 0 were established near the end of section 2, but a few more are needed. With (14), (15), (29) and the definition of S as o/r, we obtain

s = S(0, z) + O(r) where

S(0, 2) = 2 w: - au~/az. When we integrate the viscous term in (31), it win be necessary to note that ; g (r2S) = 2S(O,2) + O(t). Since thermal diffusivity will not permit a conical singularity in the temperature T = T-j)+

TgrZ+ 0(r3).

(39) field, (40)

W.E. Langlois,

328

Conservative

differencing procedures for rotationally

symmetric flow with swirl

For any quantity f which does not vanish on the axis of symmetry, we define (41) which makes available the approximation

We now repeat the analysis of the preceding section, with appropriate changes to reflect that we are dealing with the cylindrical cell Co rather than the annular C. First, Vo (@A+

(42)

= 27r =+[ruS];r,z dz. I z-

With (5) (7) (17) (38) and (40), therefore,

Next,

vo(dslj),,= zT At

IAr/2

r[wS]=- z+ &

.

0

Eqs. (17) (19) and (20) then yield

i&!GL =-&((w),&& - (w)&&). At Applying this same calculation to the coupling term yields

which presents the problem of estimating grid-cell boundaries, must be estimated themselves, viz.,

((~/t2)2)z~. These quantities, which are defined on from analogous quantities defined for the cells

(43) For example, if centered estimates are used,

W.E. Lunglois, Conservative differencing procedures for rotationally symmetric ffow with swirl

329

It is consistent with (22) to take

The higher order procedure which led to (26) should be accompanied by a more complicated form here. Carrying out the integration (43) with the quadratic and quartic terms both retained in fl leads to . @a& + - 8 (a,j - lOfl~,j)(1102,j - 1163fi~,j)l.i = (Ar)4 (9Arr

iI*

(CH 7

Applying an integration

similar to (42) to the buoyancy term in (31) yields

(ASIJ)b = 8ag (A42

At

[fT

_ (T)o].

) RI+

The same technique, applied to the radial viscous term, involves the indeterminate resolved by (39). We have

form

Using (39) and (41) to evaluate the second term, we obtain

‘dsd:‘m= 8~ (Ar)2

[ (+

$

(r2S))Rl+

-

2(s)4

(45)

The first term, together with its counterpart in (36) represents an exchange of Svanberg vorticity between Co and its annular neighbor. As in section 3, the factor of eight reflects the volume-ratio of these two cells. The presence of the second term in (45) implies that, in a viscous fluid, Svanberg vorticity is created or destroyed at the axis of symmetry. If S were a physically conserved quantity, this would be a disturbing result because, as far as the three-dimensional flow field is concerned, x = y = 0 is just one more vertical line. One would expect physically conserved quantities to enter the flow field only at boundaries, or perhaps where external force fields act as sources, as in (44). The remaining term in (31) entails no difficulty. We obtain

6. Thermocapillary Procedures

flow in Czochralski

growth of silicon crystals

based on the foregoing analysis were applied to the flow in a typical Czochralski

330

W.E. Langlois, Conservative differencing procedures for rotationally symmetric flow with swirl

crucible for growing crystals of elemental were as follows:

silicon. The physical and geometrical

parameters

Melt density: 2.33 gm/cm3, Specific heat: 0.233 cal/gm-K, Shear viscosity: 0.0233 dyne-sec/cm2, Thermal diff usivity : 0.125 cm2/sec, Volumetric expansion coefficient: 0.000014 l/K, Surface emissivity: 0.318 Surface tension temperature coefficient: 0.43 dynes/cm-K, Crucible radius: 11.55 cm, Crystal radius: 3.75 cm, Melt depth: 10.91 cm, Crucible rotation rate: 1.57 radians/set, Crystal rotation rate: -2.31 radianslsec (counter-rotation), Crucible temperature: 1773.0 K, Crystal temperature: 1685.0 K. The simulation was carried out on a 41 x 60 grid, axially stretched [3] to concentrate cells near the top and bottom. Specifically, there were three uniform regions: Zj = (j - 1)X (0.05455 cm), Zj =

1.302 cm + (j - 13) X (0.2395 cm),

zj = 10.91 cm - (60 - j) x (0.04546 cm),

grid

j= 1,2 7 -*-, 6; j = 13,14, . . . ,48; j = 55,56, . . . ,60;

these were separated by two regions in which the grid spacing diminished by arithmetic progression from coarse to fine. The uniform radial grid spacing Ar was 0.2886 cm. Upstream estimates were used for the flux-weighted averages and centered estimates were used for the unweighted averages appearing in non-advective terms. Boundary conditions and heat transport were treated by the methods described in [3]. As frequently happens in simulations of Czochralski bulk flow [2, 41, the motion ultimately settled down to an oscillation rather than a steady state. The oscillation is more quantitative than qualitative, however, and fig. 4 is quite representative of the long-term flow field. Scanning the streamline plot from crucible wall to axis of symmetry reveals the following features of the meridional circulation: 1. An intense eddy driven by thermocapillarity. The surface flow is a rapid motion inward from the crucible wall. The return flow begins as a shallow layer but then forms a Taylor column along the crucible wall. 2. A weak back-eddy driven by the return flow of the thermocapillary circulation. 3. A broad region of near stagnation. 4. A columnar eddy driven mainly by crucible rotation, with some help from buoyant convection. The intensity of this eddy is reflected in the temperature plot, which shows isotherms displaced significantly downward by advective heat transfer in this region. 5. A weak but deep eddy generated by crystal rotation.

331

W.E. Langlois, Conservative differencing procedures for rotationally symmetric flow wrth swirl

a /

10.0

P 6.0

8.0 -

/: 6.0t z 4.0-

4.0

2.0

2.0

0.0 0.0

?I J 2.0

6.0 r-

4.0

/

0.0 0.0

1

6.0

2.0

Fig. 4. The flow field: a. Streamlines. b. Isotherms.

The simulation was then continued with the surface tension temperature coefficient set to zero, i.e., with thermocapillarity artificially turned off. The long-term flow then became a weak oscillation about that shown in fig. 5. With no thermocapillary flow, the region of near stagnation was the entire annulus of melt lying under the free surface. Somewhat surprisingly, the centrifugal eddies were significantly modified. The flow generated by crystal rotation was b -

a

6.1

f 6.0

t 6.0

2

z

4.0

2.0

0.0

4.0

6.0 r-_,

6.0

10.0

0.0 I 0.0

2.0

4.0

6.0 r--*

6.0

Fig. 5. The flow field with thermocapillarity turned off. a. Streamlines. b. Isotherms.

10.0

332

W.E. Langlois,

Conservative

differencing procedures for rotationally

symmetric flow with swirl

compressed into a narrow boundary layer under the crystal face. The eddy generated by crucible rotation (and buoyancy) expanded to fill the space, and became weaker in the process. The downward displacement of isotherms was correspondingly less pronounced. It seemed remarkable that thermocapillary flow, which occurs mainly near the crucible wall, could have such a strong effect on the flow beneath the crystal. Accordingly, a series of simulations was carried out with the surface tension temperature coefficient systematically reduced from the published [19] value for silicon value to zero. The flow patterns did indeed evolve smoothly from those of fig. 4 to those of fig. 5. The size of the thermocapillary eddy seems to exert a positioning effect on the other eddies within the melt.

References Langlois and C.C. Shir, Digital simulation of flow patterns in the Czochralskicrystal pulling process, Computer Methods in Appl. Mechanicsand Engineering 12 (1977)14%152.

[l] W.E.

[2] W.E. [3] [4] [5] [6] [7] [8] [9] [lo]

[ll] [12] [13] [14] [15] [16] [ 17] [18] [19]

Langlois, Digital simulation of Czochralski bulk flow in a parameter range appropriate for liquid semiconductors, J. of Crystal Growth 42 (1977) 386-399. W.E. Langlois, Czochralski bulk flow in the growth of garnet crystals, Proceedings of the Sixth International Conference on Numerical Methods in Fluid Dynamics, Springer-Verlag (1978) 339-344. W.E. Langlois, Effect of the buoyancy parameter on Czochralski bulk flow in garnet growth, J. of Crystal Growth 46 (1979) 743-746. W.E. Langlois, Digital simulation of Czochralski bulk flow in microgravity, J. of Crystal Growth 48 (1980) 25-28. Lynn 0. Wilson, On interpreting a quantity in the Burton, Prim and Slichter equation as a diffusion boundary layer thickness, J. of Crystal Growth 44 (1978) 247-250. Lynn 0. Wilson, A new look at the Burton, Prim and Slichter model of segregation during crystal growth from the melt, J. of Crystal Growth 44 (1978) 371-376. Lynn 0. Wilson, The effect of fluctuating growth rates on segregation in crystals grown from the melt, J. of Crystal Growth 48 (1980) 43w58. D. Schwabe and A. Scharmann, Thermocapillary convection in crystal growth melts, Letters in Heat and Mass Transfer 7 (1980) 283-292. V.G. Babskiy, I.L. Sklovskaya and Yu.B. Sklovskiy, Thermocapillary convection in weightless conditions, in: G.S. Pisarenko (ed.), Space Studies in the Ukraine 1 (Naukova Dumka, 1973) 62-68. Available in English Translation as NASA Technical Translation F 15,535. C. Truesdell and R.A. Toupin, The classical field theories, in: S. Flugge (ed.), Handbuch der Physik, III/l (Springer-Verlag, l%O) 226-793. W.E. Langlois, Slow Viscous Flow (Macmillan, 1964). D.T. Kyrazis, A numerical modeling of the fluid dynamic stability of Hagen-Poiseuille flow, Lawrence Livermore Laboratory Report UCRL-52289 (1977). H.Z. Barakat and J.A. Clark, Analytical and experimental study of the transient laminar natural convection flows in partially filled liquid containers, Proc. 3rd International Heat Transfer Conf. 2 (1966) 152. K.E. Torrance, Comparison of finite-difference computations of natural convection, J. of Research of the National Bureau of Standards-B. Mathematical Sciences 72B (1%8) 281-301. N. Kobayashi and T. Arizumi, Computational analysis of the.flow in a crucible, J. of Crystal Growth 30 (1975) 177-184. N. Kobayashi, Computational simulation of the melt flow during Czochralski growth, J. of Crystal Growth 43 (1978) 357-363. N. Kobayashi, Computer simulation of heat, mass and fluid flows in a melt during Czochralski crystal growth, Computer Methods in Appl. Mechanics and Engineering 23 (1980) 21-33. C.E. Chang and W.R. Wilcox, Inhomogeneities due to thermocapillary flow in floating zone melting, J. of Crystal Growth 28 (1975) 8-12.

W.E. Langlois, Conservative differencing procedures for rotationally symmetric j?ow with swirl

333

Note added in proof Professor Pieter Wesseling has pointed out that the angular momentum balance used in this paper is reminiscent of Isaac Newton’s attempt (Lib. II, Prop. LI of the Principia) to calculate what is now known as axisymmetric Couette flow, by examining the interactions between neighboring cylindrical laminae. Regrettably, Newton fell into the elementary error of balancing forces rather than torques, which led him to an incorrect result.