Comparison of three turbulence models for simulation of melt convection in Czochralski crystal growth of silicon

Journal of Crystal Growth 205 (1999) 71}91

Comparison of three turbulence models for simulation of melt convection in Czochralski crystal growth of silicon Aleksey Lipchin, Robert A. Brown* Department of Chemical Engineering, Massachusetts Institute of Technology, Cambridge, MA 02139-4307, USA Received 15 June 1998; accepted 22 March 1999 Communicated by J.J. Derby

Abstract Accurate simulation of heat transfer and oxygen segregation in large-scale Czochralski growth of silicon requires modeling of turbulent convection in the melt. We test three models for turbulent convection for the calculation of the #ow, temperature "eld and oxygen transport in a prototypical model of the Czochralski melt. Each of the models is based on the k}e turbulence model, but di!ers in the form of model and the description of the #ow near solid boundaries. The three formulations are: (I) a traditional k}e model using wall functions at solid boundaries, (II) a k}e model with a one-equation model for the #ow near solid boundaries, and (III) a low-Reynolds number k}e model that does not require an independent description of the #ow near the wall. Calculations using a "nite-volume method with collocated grids shows large di!erence in the predictions of these models, especially for the highly separated #ows driven by buoyancy. In the limit of week turbulence only the low-Reynolds number model gives solutions which are consistent with the laminar regime #ow.  1999 Elsevier Science B.V. All rights reserved. PACS: 47.27.!i; 47.27.E; 81.10.Aj; 81.10.Eq Keywords: Turbulence models; Czochralski; Silicon

1. Introduction Large-scale Czochralski (CZ) crystal growth is the dominant method for producing the singlecrystalline silicon that is the mainstay of the microelectronic industry. The demand for increasing larger silicon wafers has driven Czochralski crystal growth from 150 mm (6) crystal diameters typical in 1985 to 200 mm (8) crystals today and to the

* Corresponding author. Fax: #1-617-253-5726. E-mail address: [email protected] (R.A. Brown)

deployment of 300 mm (12) crystals in the very near future [1]. As the size of the crystal has increased, so has the volume of melt, the size of the crucible and the size and cost of the entire growth apparatus. This increased cost is coupled with the ever more stringent demands on crystal quality, measured in terms of the oxygen concentration, uniformity, and the type, size and distribution of microdefects in the crystal [2]. The cost and complexity of large-scale Czochralski systems has made numerical simulation an important tool for optimizing the design and operation of these systems. Beginning with Derby and

0022-0248/99/$ - see front matter  1999 Elsevier Science B.V. All rights reserved. PII: S 0 0 2 2 - 0 2 4 8 ( 9 9 ) 0 0 2 4 2 - 0

72

A. Lipchin, R.A. Brown / Journal of Crystal Growth 205 (1999) 71}91

Nomenclature g k Nu p r Sh ¹ u u   u O v w y z l  e p q  W X

gravitational constant turbulent kinetic energy Nusselt number normalized by its value for pure conductive (*"0) regime pressure radial cylindrical coordinate Sherwood number (analog of Nu for mass transfer) temperature radial velocity component velocity component parallel to the wall friction velocity azimuthal velocity component vertical velocity component distance from the nearest wall vertical cylindrical coordinate turbulent viscosity dissipation rate of turbulent energy Stefan}Boltzmann constant shear stress at the wall stream function rv, swirl

Subscripts bt bottom surface (crucible bottom) cr crystal surface (crystal/melt interface) fr free surface (gas/melt interface) sd side surface (crucible side)

Brown [3}6], thermal-capillary models of Czochralski growth have been developed that model conductive and radiative heat transfer in the crystal, melt, crucible and the surroundings components of the furnace and that include the self-consistent calculation of the melt/crystal and melt/gas interfaces, as well as the crystal shape; see Refs. [7}9] for a review of these models and the "nite element analysis that is typically used for their solution. Simulations of heat transfer using thermal-capillary models are routinely used for prediction of temperature "eld in the crystal and for design of the heat transfer elements [7]. The accuracy of these calculations depends on the level of detail included in the heat-transfer model, on the accuracy of the thermophysical properties and on the accuracy of the

description for convection in the melt. In this last feature that is the focus of this paper. Modeling convection in even moderate-scale Czochralski systems for silicon growth is a very di$cult problem. Although early extensions of thermal-capillary analyses to include convection modeled the #ow driven by buoyancy force, crystal and crucible rotation, and surface-tension gradient as laminar and steady state [10], there is ample evidence that the #ows in these systems are at least chaotic and, most probably turbulent. There is growing direct evidence of the form of convection in large-scale Czochralski silicon melts. The scaling analysis in Refs. [11,12] give dimensionless driving forces which are indicative of turbulent convection. Even so, there is insu$cient

A. Lipchin, R.A. Brown / Journal of Crystal Growth 205 (1999) 71}91

direct numerical simulation of intensive convection in large-scale Czochralski melts to prove this assumption. The axisymmetric time-dependent simulations of Kim and Langlois [13] on relatively small-scale crucibles do show chaotic #ow patterns. Another evidence of turbulent convection in Czochralski melts is the experimental data of Kim and Smetana [14,15] on oxygen segregation where chaotic #uctuations of oxygen concentration were observed along the center line of crystal. The measurements of temperature #uctuations in CZ systems with crucible diameters up to 600 mm) by Seidl et al. [16] show that the #ow has the character of the transition from soft to hard turbulence. Incorporating turbulent convection into practical simulations of convection and heat transfer in CZ growth must rely on the use of single-point turbulence models. These models are based on approximations for the Reynolds average equations of momentum, mass, heat and species transport * velocity, pressure, temperature, and concentration * that are solved by using approximations for the couplings between the time #uctuating components of these variables [17]. The well-known family of k}e models describes the Reynolds stress tensor for the pairwise coupling between #uctuating velocity components through the solution of two additional model equations for the kinetic energy k and dissipation function e [18]. In this context, the calculation of turbulent convection in CZ system is reduced to solution of equations for the time-average "eld variables in addition to the equations for k and e. With the usual quasisteady-state crystal growth approximation [19] and appropriate thermal boundary conditions it may be reasonable to imagine these #ows as steady and axisymmetric. Even so, accurate prediction of turbulent convection in a CZ melt is extremely di$cult. Most importantly the combinations of driving forces for the #ow and the e!ects of the boundaries make the #ows highly anisotropic and lead to #ow separation where large #ow structures attach to solid surface. Both conditions are notoriously di$cult for k}e models, because the description of the Reynolds stress tensor in terms of k and e models these terms as isotropic. The occurrence of #ow separation points to perhaps the most di$cult feature for k}e models;

73

accurate #ow modeling requires the coupling between this representation for the fully turbulent core #ow and the region of the #ow near walls where the #ow must transit to a viscous sublayer. At the Reynolds numbers appropriate for fully developed turbulence the transition and viscous layers are very thin and are the focus of much of the modeling e!ort based on k}e models. Three types of approaches have been used. These include the classical k}e model formulation in which boundary conditions along solid surfaces are modeled using wall functions constructed from the structure of turbulent shear #ow [18], two-layer models that replace the use of wall function near solid surfaces with a one-equation model for the #ow near the wall [20,21], and, "nally, low-Reynolds number k}e models that are designed to model the transition from the turbulent core to the wall without the use of special near-wall models [22}25]. The purpose of this paper is to compare the applications of representative k}e turbulence models of these three types for the calculation of the #ow in a prototype model for melt convection in the Czochralski system. The model is based on a cylindrical melt volume with "xed surface with convection driven by buoyancy di!erence and crystal and crucible rotation. The turbulent melt #ow is assumed to be axisymmetric and steady. However, there are experimental and numerical evidence in the literature that the #ow in large-scale CZ system does have large-scale, three-dimensional and timedependent structures. The limitations of the axisymmetric and steady-state approach to CZ #ow modeling for estimating dopant transport and heat transfer are still open questions. Our goal is not to address these issues, but to look into the applicability of speci"c turbulence models for simulating #ow in CZ systems. This study reports one of the "rst [26] direct comparisons of turbulence models for simulating convection in this geometry. There have been few calculations of turbulent convection in Czochralski growth. Previously, Kobayashi [27,28] has used a standard k}e model to compute convection in a prototype Czochralski con"guration and compared these calculations to measurements of temperature #uctuations in a silicon melt. Although direct detailed comparisons between simulations and experiments were not

74

A. Lipchin, R.A. Brown / Journal of Crystal Growth 205 (1999) 71}91

reported, the simulations did reproduce the decrease in the temperature of melt surface seen with increasing crucible rotation. Chung et al. [29] applied the standard k}e model with wall functions for a similar prototype Czochralski system and reported good agreement of calculated oxygen distribution in the crystal with experimental data. The only comprehensive simulations of heat transfer in Czochralski crystal growth to include modeling turbulent convection were performed by Kinney et al. [30] and were based on a lowReynolds number k}e model proposed by Launder and Sharma [23,31]. Because the simulations included a complete thermal-capillary model for heat transfer in an experimental Czochralski system, results for temperature measurements for a thermocouple probe and for the shape of the melt/crystal interface could be compared directly with measurements. The agreement was reasonably good. Ristorcelli and Lumley [32] have used a secondmoment di!erential closure Reynolds stress model for calculations of melt convection in a prototype Czochralski melt. This model is extremely complicated * including tensorial second-order di!erential moment equations and third-rank algebraic models for the third-order moments * and was not used in extensive calculations. Finally, Zhang et al. [26] compared performance of four turbulence models for CZ #ow driven by buoyancy and thermocapillary force: the standard k}e model with wall function, the low-Reynolds number k}e model of Abe and Kondoh, the renormalization group k}e model of Yakhot and the algebraic stress model by Gatski and Speziale (see Ref. [26] for references). All the models gave similar results for velocity and temperature distributions. The later two models incorporate nonisotropic and nonlinear Reynolds stress tensor. However, rotation of the crystal and crucible were not considered, which are important contributions to the nonisotropic nature of the #ow. Both the renormalization group formulation and the algebraic stress model are high-Reynolds formulations and are used with wall functions and, therefore, are limited by the drawbacks of the wallfunction approach; see Section 2.2. In this paper we compare three models for turbulent convection in Czochralski growth of silicon

for various rotation rates of the crystal and the crucible: (I) A standard k}e model using wall functions along solid surfaces [33]. (II) A two-layer k}e model that includes a oneequation model for the #ow near walls [21]. (III) The low-Reynolds number k}e model proposed by Jones and Launder [24]. Simulations are based on a transient "nite volume method with collocated grids, as "rst developed by PericH [34]. The turbulence models are described in Section 2 along with the prototype geometry for Czochralski bulk #ow. The three models are compared in Section 3 in a series of calculations for crystal and crucible rotation alone, for buoyancy-driven convection, and for the combination of these driving forces. Predictions are reported for the #ow, for the temperature in the melt and for the species transfer for oxygen.

2. Turbulence models 2.1. General formulation Three turbulence models are compared using the bulk #ow model [35,36] for melt convection in the CZ system. The idealized con"guration is shown in Fig. 1 and consists of a vertical cylindrical crucible,

Fig. 1. Schematic diagram of prototype Czochralski con"guration and mesh used for the JL and TL model computations.

A. Lipchin, R.A. Brown / Journal of Crystal Growth 205 (1999) 71}91

75

Table 1 Geometrical parameters, physical properties of the melt, operational parameters, and dimensionless groups Variable

Symbol

Value

Source

Crucible radius Crucible height Crystal radius Heat capacity Thermal conductivity Density Thermal di!usivity Viscosity Kinematic viscosity Coe$cient of thermal expansion Surface emissivity Melting point temperature Crucible temperature Ambient temperature Crystal rotation rate Crucible rotation rate Prandtl number Grashof number Radiation parameter

R  H R  C N k 2 o a 2 k l "k/o

b e 2 ¹

¹ #*¹

¹ X  X  Pr"l /a

2 Gr"bg¹ R/l

 R"R ¹ e p/l oC  2 N

20 cm R  0.4R  1.0 Jg\ K\ 0.64 W cm\ K\ 2.42 g cm\ 0.26 cm s\ 0.7 cP 2.9;10\ cm s\ 1.4;10\ K\ 0.05 1683 K *¹"50 K ¹

0, 20, 100 rpm 0, !8, !20 rpm 0.011 2;10 3.9

[37] [38,36] [38,36] [10] [10] [10] [10] [10] [10] [36] [36] [10] [37] [35,32] [30,36] [30,36]

of height H and radius R , rotating at angular  velocity X and "lled by melt. The melt is bounded  above partially by a crystal, of radius R , rotating at  angular velocity X and by a #at free surface. The  temperature of the crystal surface is set to the melting temperature ¹ and the temperature of

the crucible is set to ¹ #*¹. The heat #ux from

the melt free surface is governed by Stefan} Boltzmann equation, i.e. by radiation to an imposed ambient environment at temperature ¹ . Other assumptions made in the bulk #ow model are that all the thermophysical properties of melt are independent of temperature and have the values listed in Table 1 and that the #ow is axisymmetric. Because the dependence on temperature of the surface tension is neglected, thermocapillary motion is not included in our analysis. Other calculations including this e!ect for large diameter CZ system show little e!ect of thermocapillary motion on the intensely turbulent #ow that is predicted near the meniscus in these systems. The Reynolds-averaged equations of motion and the turbulence models are written in a cylindrical coordinate system centered at the bottom of the crucible, as shown in Fig. 1. Dimensionless vari-

ables are formed by scaling length with the radius of the crucible, R , time with R/l (Gr, velocity   with l (Gr/R , pressure with ol (Gr/R and



 temperature with ¹ , where o is the melt density,

l is the melt molecular kinematic viscosity and Gr

is the Grashof number de"ned in Table 1. In the k}e turbulence models, the dimensionless kinetic energy is scaled with l Gr/R and the dissipation rate

 with l Gr/R.

 The dimensionless Reynolds-average equations of motion, continuity and energy transport are




(Gr



Ru Rp X #* ) u " ) (l u)! #(Gr Rt Rr r Rl Ru Rl Rw u #  #  !l , (1) Rr Rr Rz Rr r




(Gr



Rw Rp #* ) w " ) (l w)! Rt Rz Rl Ru #(Gr(¹!1)#  Rr Rz Rl Rw #  , Rz Rz

(2)

76

(Gr

A. Lipchin, R.A. Brown / Journal of Crystal Growth 205 (1999) 71}91






RX #* ) ) " ) (l X) Rt






2 RX Rl ! l #X  , r Rr Rr

) *"0,

(3) (4)















 

1 l #  ¹ , (5) Pr p 2 where *2"(v , v , v )"(u, w, v), l,1#l and X, P X F  vr is the swirl function [35]. The equations of k}e model are written in the form R¹ (Gr #* ) ¹ " ) Rt

 

l 1#  k p I #P #G !e#D, I I Re l (Gr #* ) e " ) 1#  e Rt p C #(C f (P #C G ) C  I C I e !C f e) #E, C  k (Gr

Rk #* ) k " ) Rt








(6)

 

(7)

where the turbulent production rates due to dissipation P and buoyancy G are de"ned as I I l l R¹ P , ( *# *2) : ( *# *2), G ,!  (8) I 2 I p Rz 2 and the turbulent kinematic viscosity l is  k l "Gr f C . (9)  I Ie Closure of Eqs. (1)}(7) requires setting the remaining constants and functions in the model and gives some of the distinctions between the three k}e models tested here. In all the models we use the generally agreed values of C "0.09, p "1.00, I I p "1.30, C "1.44, C "1.92, and p "0.9. C C C 2 There is no agreement on the value of C [39}43]. C We use the approach of Henkes and collaborators [42,44,45] and set C "tanh("w"/(u#v). The C choice of this value only weakly a!ects the results. Some researchers [33] have neglected the term G in Eq. (7). I

The damping functions ( f , f , f , D, E) are used   I in the low-Reynolds number models to describe near-wall turbulence. To properly make the transition from the wall to the turbulent core, these terms should satisfy the high-Reynolds number limit in the core f "f "f "1, D"E"0. (10)   I In the standard k}e model the values listed in Eq. (10) are used everywhere and special models are used for the near-wall region of the #ow. The boundary conditions at the crystal/melt interface, z"H/R , 0)r)R /R , are    ¹"1, u"0, w"0, X"X r. (11)  The conditions on the crucible surface, z"0, 0)r)1 and r"1, 0)z)H/R , are  ¹"1#*¹/¹ , u"0, w"0, X"X r.

 (12) Along the meniscus surface R /R (r(1) the conditions are   Ru RX Rk Re w"0, " " " "0, Rz Rz Rz Rz




(z"H/R ,  (13a)



l R¹ 1 #  "R(¹!(¹ /¹ )), (13b)

Pr p Rz 2 where R is the radiation number (see Table 1) and ¹ is the dimensionless ambient temperature. We set ¹ "1, i.e. the ambient temperature is set to the melting temperature; the choice of ¹ has little qualitative e!ect on the results because of the relatively low emissivity of the silicon melt. The conditions in Eq. (13a) on k and e set the normal #ux to zero, but allow the #ow to remain turbulent at the melt/gas surface. Conversely, one can argue that the #ow is constrained to be two-dimensional at the surface and so a viscous sublayer should be present. Calculations performed replacing Eq. (13a) with the conditions on k and e appropriate for a rigid surface showed little quantitative variation in the bulk of the #ow. The line (r"0, 0(z(H/R ) is a symmetry  axis:

!

u"X"0,

R¹ Rw Rk Re " " " "0. Rr Rr Rr Rr

(14)

A. Lipchin, R.A. Brown / Journal of Crystal Growth 205 (1999) 71}91

2.2. The standard k}e model with the wall functions: WF Using wall functions with the standard k}e model is meant to avoid having to resolve the transition from the turbulent core to the viscous wall layer by modeling the transition with simple models valid for a parallel shear #ow and matching these models to the core #ow computed with the k}e model. This approach is presented in many texts [18,33] and is only described brie#y here. The model near the wall is described in the wall coordinates de"ned as y>,(Gru y, u ,Gr\(q , O O 

(15)

where q is the wall shear stress and u is the friction  O velocity. u is computed from the transcendental O relationship between the tangential velocity u , u   O and y: u 1  " ln(9y>). (16) u 0.41 O This relationship is valid outside the viscous sublayer, or for y>'11.5, where the molecular viscosity is negligible compared to l . At this point  (closest grid point to the wall) the boundary conditions on k and e are u u k" O , e"(Gr O . (17) 0.41y (C I No-slip boundary conditions on velocity are applied on the wall. In the "nite volume method all boundary conditions are implemented by calculating #uxes through the faces of "nite volumes at boundaries. In the wall-function formulation the momentum #ux at the solid boundary, e.g. shear stress q , is calculated by Eqs. (15)}(17). Calculating  the shear stress from the normal derivative of the tangential velocity at the wall (u /y) is inconsist  ent with the wall-function approach. The shear stress q computed in the traditional way can be  several time less than q . In the stagnation regions  (at the corners of the crucible or at the central line) q is less than q and we use q instead of q in these     regions. It is not surprising that where q "q co 

77

incides to points where y>"11.5, i.e. boundaries of regions where the wall functions are valid. The boundary conditions on a solid surface are completed by assuming that for a low Prandtl number #uid, conductive heat transfer dominates near the wall so that the classical boundary conditions on ¹ at a solid surface are unchanged. This approach avoids adding a wall function for temperature as is usually done for high Prandtl number #uids [33]. There are many di$culties with applying wall functions to the CZ #ow problem; even the applicability of wall functions for convective #ows has been the topic of much discussion [33,41,42,46]. Comparison of predictions of the k}e model with wall functions for turbulent natural convection of air (Pr"0.71) and water (Pr"7.0) for thermal boundary layers and enclosures are described in several references [33,42,43,45], where simulations with wall functions are shown to over-predict the wall heat transfer by as much as 30}50%. Moreover, the numerical solution is sensitive to the distance of the closest grid point to the wall. Researchers [40,43,44,47] have tried replacing Eq. (17) with e<1 and k"0 at the wall; however this approach does not lead to better predictions. The strongly rotating #ow near the wall in CZ system is an additional source of concern, because the derivation of Eq. (17) assumes a one-dimensional, not a two-dimensional, averaged motion near the wall [20]. 2.3. The two-layer model: TL Iacovides, Launder and collaborators [20,21] have addressed several of these issues by developing a two-layer k}e model for simulation of rotating and curvilinear #ows. The two-layer models consist of an algebraic stress model or the standard k}e model for the core #ow and a low-Reynolds, one-equation model or a simple mixing length hypothesis for the near-wall regions. Calculations for orthogonally rotating U-bend presented in Refs. [20,21] gave promising results compared to results based on using wall functions. Following Refs. [20,21] we have used the standard k}e model for the core #ow and Wolfshtein's one-equation model [48] for the near-wall regions.

78

A. Lipchin, R.A. Brown / Journal of Crystal Growth 205 (1999) 71}91

In Wolfshtein's approach Eq. (6) is used to compute k near the wall and e and l are obtained as  k e"(Gr , l "C l (Grk (18a)  II l C by estimating the length scales l and l from C I l "2.55y(1!exp(!0.263yH)), (18b) C l "2.55y(1!exp(!0.016yH)), (18c) I yH"2.55(Gryk,

(18d)

2.4. Low-Reynolds model of Jones}Launder: JL Low-Reynolds number forms of k}e models have been introduced to give smooth transitions between the turbulent core #ow and the wall [18,22] and have shown improved performance relative to standard k}e models for isothermal near-wall #ows [23], for the convective boundary layer caused by a heater plate [33,43,45] and for natural convection in enclosures [25,40]. For example, in the calculations of natural convection for air in an enclosure, the predicted Nusselt numbers for the Rayleigh number Ra"4;10 di!ered from measurements by only 0.5}1.6%, giving the best predictions [25]. The damping functions used in the Jones}Launder form of the low-Reynolds model [24] are




k Re "Gr ,  e



(19a)



 

D"!2

R(k  R(k  # , Rr Rz (19b)








1 Ru u  1 Rv v  E"l 3 ! #3 !  r Rr r r Rr r #

2.5. Oxygen transport model A simple model for oxygen transport in the melt is used as a vehicle for comparison of the #ow predicted by turbulent transport. Here oxygen is treated as a dilute species with nondimensional concentration C(r, z) governed by the turbulent transport equation




1 Rw  1 Ru  1 Rv  #2 #2 , (19c) r Rr r Rz r Rz

where the derivatives Ru/Rr, Rw/Rz have been neglected in E. The boundary conditions used at solid





 

1 l #  C , (20) Sc p 2 where Sc,l /D is the Schmidt number and

D"3;10\ cm s\ is the molecular di!usivity of oxygen. We follow Kinney et al. [30] and use the eddy di!usivity hypothesis; the turbulent Schmidt number is equal to the turbulent Prandtl number p "0.9. Very idealized boundary conditions are 2 used on the oxygen concentration. The crucible is assumed to be dissolved so that it is a source of oxygen. This is approximated by setting (Gr

respectively. At the wall, k"0 is set.

f "1.0, f "1!0.3 exp (!Re),    !2.5 , f "exp I 1#Re /50 

walls are the classical no-slip conditions on velocity and k"e"0.

RC #* ) C " ) Rt

C"1.

(21)

Oxygen transport to the gas in the form SiO is modeled by C"0,

(22)

at the planar melt/gas surface. The condition for incorporation of oxygen at the melt/crystal interface is RC ! "(k!1)PeC, Rz

(23)

where k is the equilibrium segregation coe$cient, Pe,
3. Finite volume method Finite volume methods have proved to be the most robust for solving the nonlinear partial di!erential equations governing k}e turbulence models [34]. We apply a "nite-volume method based on

A. Lipchin, R.A. Brown / Journal of Crystal Growth 205 (1999) 71}91

velocity}pressure splitting and calculation of steady-state solution by transient simulation. Discretized equations are formed on a collocated grid using the SIMPLE-like algorithm [34] for calculation of the pressure. Three schemes of di!erence approximations for convective terms are tested: central di!erence approximation, upwind di!erences and the hybrid power-law scheme developed by Kim and Ro [49] for cylindrical coordinates. The grids for application of the "nite-volume discretizations were generated slightly di!erently for the TL and JL models and for WF discretization. For the TL and JL models a grid with equidistant nodes was used everywhere, except near the melt/crystal/meniscus tri-junction and near the solid surfaces, where the grids were re"ned. Re"ned radial grids were used near the tri-junction and the grid were decreased in volume geometrically as the solid surface were approached. The mesh used in the calculations with the JL and TL models is presented in Fig. 1; the mesh had 49 control volumes in the r-direction and 69 in the z-direction. More mesh re"nement is required in the vertical direction because of the need for mesh re"nement both on the crystal surface and the bottom of the crucible. The mesh used with the WF model was similar to the one used in the bulk region with JL and TL models, but omitted the mesh re"nement near the boundaries; this mesh had 49 control volumes in the r-direction and 59 in the z-direction. Regarding the solution for the doubly re"ned mesh (99;139) as exact, errors on mesh 49;69 are computed for each quantity. All errors are within 5% for the three schemes of di!erence approximations for convective terms. The upwind di!erence

79

scheme was used in all other computations. One of the most awkward issues in the WF approach is the solution's sensitivity to the distance of the "rst grid point from the wall. We have analyzed this sensitivity and for the calculations reported in this work increasing this distance by a factor of three increased the maximum variation in the turbulent viscosity by 20% and maximum value of the stream function by 7%. We have developed a strategy of choosing this distance based on physical reasoning; however, it is doubtful that any strategy can eliminate the intrinsic shortcomings of the wall-function formulation. For the TL model some uncertainty is created by the choice of the size of the near-wall layer where the one-equation model is applied. To keep all advantages of the structured grid the boundary between the layers was taken to be planar and parallel to the wall. For each near-wall layer the value of y> (see Eq. (15)) averaged along the layer interface was calculated. The value y> monotonically increases with increasing y or tangential velocity near the surface. The standard k}e model is valid for y>'11.5 [18]. In the present study the one-equation layer size was chosen so that average y> satis"ed the condition 20)1y>2)80. Each one-equation layer (near the crystal, side crucible wall and crucible bottom) has its own thickness. However, in some calculations the melt motion near a wall was so weak that 1y>2 was less than 20; otherwise, the one-equation region spread unrealistically far from the wall. We carried out a series of computations with di!erent sizes of y for the oneequation layers adjacent to the crystal surface and the side crucible wall. The variations of y, 1y>2, and other quantities are listed in Table 2. These

Table 2 In#uence of the size of one-equation layer region on solution characteristics computed with the TL turbulence model Regime X "0, X "0,   Side wall X "20 rpm, X "!8 rpm,   Side wall X "20 rpm, X "!8 rpm,   Crystal surface

1y>2

y

l  

Nu 

W



20}60

0.014}0.062

5.8%

0.5%

1.7%

20}100

0.03}0.10

2.0%

2.6%

2.0%

20}100

0.004}0.021

0.3%

0.6%

2.6%

80

A. Lipchin, R.A. Brown / Journal of Crystal Growth 205 (1999) 71}91

results suggest that the solution becomes less sensitive to the width of the one-equation layer with increasing #ow intensity near the corresponding wall. A possible explanation is that, with increasing rotation rate, the #ow becomes more similar to a pure shear driven motion in the near-wall region. The one-equation and the standard k}e models become more consistent for pure shearing motions. Selecting C "0.8 instead of C " C C tanh((u#v/"w") in Eq. (7) a!ected slightly the results; e.g., for the TL model, with X "20 rpm  and X "!8 rpm, the maximum turbulent viscos ity varied by 3% and the maximum stream function changed by 0.6% with this change.

4. Results 4.1. Rotationally driven yows Calculations without buoyancy-driven #ow (*¹"0) were used to derive a better understanding of the #ow "eld and of the di!erences between the predictions of the models for #ows driven by crucible and crystal rotation. The contours of the stream function W(r, z) for X "20 rpm and  several crucible rotation rates (X ) are shown in  Fig. 2 for JL and TL models. The maximum and minimum values of W are also listed. Here a cellular motion with W'0 corresponds to a clock-wise motion and W(0 corresponds to a counterclock-wise #ow. The WF model predicted similar results for X "!8 rpm: (W "2.01;10\, 

 W "!0.14;10\, l "259.2). However, be    cause of the low intensity of the motion near the crucible surface the distance of the "rst grid point from the wall (the size of the control volumes nearest to the wall) had to be set to be too large to provide a reasonable value for y>. This problem leads to wiggles in the velocity "eld near the crucible. Di$culties in predicting physically smooth velocity "eld near walls is a computational problem intrinsic to the WF model and the problem is worst with decreasing crucible rotation rate. For this reason results are shown only for JL and TL models. For su$ciently high crucible counter-rotation rates the #uid body is in rigid rotation with angular

velocity close to X . The Ekman #ow at the crucible  bottom forms a vortex near the bottom. Sharp changes in the angular velocity near the crystal interface disturb the balance between centrifugal force and the pressure gradient. This leads to the radial motion towards the center line under the crystal surface and to the formation of a counterclockwise rotating eddy near the crystal surface. It is the commonly known Taylor}Proudman cell [12]. The Ekman #ow under the crystal is bounded by a small region near the center line and too weak to be resolved in the upper pair of "gures in Fig. 2. The #ow driven by crucible counter-rotation dominates over the motion caused by crystal rotation because of the di!erences in torque. Slowing crucible rotation leads to stronger involvement of crystal rotation with the #ow and to the formation of a clockwise rotating vortex under the crystal. Further decreasing the crucible rotation causes this vortex to migrate towards the free surface and grow until it "lls the whole melt. The di!erences in the predictions of JL and TL models increases with decreasing crucible rotation rate. At X "!8 rpm and X "!4 rpm the #ow   structures predicted are in qualitative agreement, although the #ow intensities and values of maximum turbulent viscosity di!er by 25%. For lower crucible rotation rates the predictions qualitatively disagree: at X "!0.5 rpm the TL predicts  a single #ow cell while the JL model predicts two cells. The di!erences in the predictions are linked to the di!erences in the turbulent viscosity distributions. The corresponding turbulent viscosity "elds are shown in Fig. 3. The TL model predicts much large values of l (r, z) near the crystal and meniscus  than the JL model; the larger value of l causes the  more intense melt motion. For X "!0.5 rpm, turbulence is su$ciently  weak that the models give completely di!erent distributions of turbulent viscosity. To understand which model better reproduces the transition to the laminar regime we computed the laminar #ow for the same set of parameters. The stream function for the laminar #ow is shown in Fig. 4. The calculations are performed with a laminar version of our "nite volume code using time integration and the same (49;69) mesh. A steady-state solution is not reached; small oscillations of the #ow persist near

A. Lipchin, R.A. Brown / Journal of Crystal Growth 205 (1999) 71}91

81

Fig. 2. Contours of the stream function computed for varying crucible rotation rate (X ) without buoyancy-driven convection for the JL  and TL turbulence models: X "20 rpm. Stream function is scaled with l (Gr. 

82

A. Lipchin, R.A. Brown / Journal of Crystal Growth 205 (1999) 71}91

Fig. 3. Contours of the turbulent viscosity l computed for varying crucible rotation rate (X ) without buoyancy-driven convection for   the JL and TL turbulence models: X "20 rpm. The turbulent viscosity is scaled with l . 

A. Lipchin, R.A. Brown / Journal of Crystal Growth 205 (1999) 71}91

Fig. 4. Contours of the stream function computed for laminar #ow with X "20 rpm and X "!0.5 rpm without buoy  ancy-driven convection.

the wall. Interestedly, the stream function of the laminar solution is smooth in the region of the melt where the JL model predicts nearly laminar #ow, i.e. the turbulent viscosity is close to zero, and there are wiggles in the laminar #ow in the region where the JL model predicts weak turbulence. 4.2. Buoyancy-driven yow For pure buoyancy-driven #ow all models predict the same patterns for the stream function and temperature "elds. Contours of these variables are shown in Fig. 5. Besides the large rotating cell there is a small counter-rotating eddy in the bottom corner of the crucible. This second eddy is too weak to be resolved in Fig. 5, but its existence is con"rmed by "nite-element computations [50]. The turbulent viscosity "elds are compared in Fig. 6 for the three models. The location of the maximum in

83

the turbulent viscosity corresponds approximately to the location of the maxima of the turbulent kinetic energy k (see Fig. 6) and the dissipation rate e. These variables are set by turbulence generation by the vertical temperature gradient; see Eqs. (6) and (7). The distributions of the turbulent viscosity for the JL and WF models are much closer to each other than the results for the TL model. The same conclusion can be drawn from various solution characteristics listed in Table 3. It is also clear that the WF model predicts a more turbulent #ow, compared to the JL model and as measured by the turbulent viscosity. As was mentioned above, the same conclusion has been reported for buoyancydriven convection in two-dimensional enclosures [40,43]. 4.3. Combine rotation and buoyancy driven yow The combination of crucible rotation and buoyancy-driven #ow generate a counter-clock-wise vortex while crystal rotation causes clockwise rotation. Even slow crucible rotation suppresses the motion induced by crystal rotation; see Fig. 2. To illustrate the interaction of the #ow induced by crystal rotation and buoyancy di!erence calculations were carried out for X "100 rpm and  X "0. The results for the JL model are shown in  Fig. 7. The strong azimuthal motion leads to values of the turbulence viscosity that are much larger than the ones predicted for buoyancy-driven convection alone. This leads to a steep thermal boundary layer adjacent to the crystal surface. The heat #ux from the free surface is less than the values predicted for conduction alone (Nu (1) because 

Fig. 5. Contours of the stream function and temperature computed for buoyancy-driven convection without rotation for the JL turbulence model.

84

A. Lipchin, R.A. Brown / Journal of Crystal Growth 205 (1999) 71}91

Fig. 6. Contours of the turbulent viscosity l for the (a) JL, (b) TL and (c) WF turbulence models and the (d) turbulence kinetic energy  for the JL model computed for buoyancy-driven convection without rotation.

Table 3 Comparison of predictions of the three turbulence models with buoyancy-driven convection and varying crystal and crucible rotation rates Rotation model

Nu  Nu  Nu
 Nu  W ;10

 !W ;10

 l   k ;10

 Sh  Sh ;10\  Sh ;10\
 Sh ;10\  C ;10 

X "0 X "0  

X "20 X "!8 (rpm)  

X "100 X "!20 (rpm)  

JL

WF

TL

JL

WF

TL

JL

JL

WF

TL

JL

9.63 2.89 27.8 1.26 69.7 0.70 459 1.85 1.18 1.16 3.74 2.21 1.08

9.19 2.88 26.1 1.27 72.2 5.48 590 1.60 } } } } }

7.23 1.92 22.0 1.27 62.9 5.81 253 0.95 0.24 0.86 2.44 1.54 0.26

6.62 4.96 8.71 1.17 2.21 0.23 2909 1.30 4.09 3.98 6.39 5.74 2.22

6.43 4.68 8.99 1.17 2.35 3.50 3207 1.29 } } } } }

5.84 4.34 7.89 1.11 2.87 2.75 2467 1.19 0.69 0.53 1.08 0.82 0.31

} } } } 1.94 1.19 250 0.27 } } } } }

7.07 4.29 12.8 1.13 9.04 7.03 1342 4.53 3.15 1.25 2.30 1.78 1.15

6.64 3.96 12.3 1.14 9.34 7.84 1267 2.48 } } } } }

6.84 4.51 11.1 1.10 10.1 14.6 1270 7.06 3.23 1.17 2.15 1.67 1.15

} } } } 7.45 15.9 841 4.76 } } } } }

Results for JL model without buoyancy-driven #ow (*¹"0).

of the outward motion along the relatively cold melt/crystal interface. Two #ow regimes are investigated: #ow with moderate rotation rates (X "20 rpm, X "  

!8 rpm) and #ow with high rotation rates (X "  100 rpm, X "!20 rpm). Flows computed for  moderate rotation are shown in Figs. 8 and 9; the corresponding values of #ow parameters are

A. Lipchin, R.A. Brown / Journal of Crystal Growth 205 (1999) 71}91

85

Fig. 7. Contours of the stream function and temperature computed for the JL model with buoyancy-driven convection and X "100 rpm and X "0; for this solution l "2279.    

Fig. 8. Contours of the stream function for the (a) JL, (b) TL and (c) WF turbulence models and the (d) temperature for the JL model calculated for buoyancy-driven convection and rotation with X "20 rpm and X "!8 rpm.  

reported in Table 3. As was shown above, the counter-rotating eddy dominates for large crucible rotation rates. The #ow patterns for all models are similar to each other because the models di!er only in near-wall #ow treatment. Turbulence generation mostly occurs in the near-wall regions and leads to some quantitative di!erence in the predictions. As was for the case for buoyancy-driven #ow, the TL model predictions are farther from the results of the JL and WF models than results of the JL and WF are from each other. However, for moderate rotation rates the #ow pattern near the crucible side

surface predicted by the WF model is qualitatively di!erent from the ones predicted by the JL and TL models. The predicted temperature distributions are nearly the same for all models because of the low value of the Prandtl number. The distributions of azimuthal velocity and turbulent energy also are similar for all models. The #ow pattern is very similar to the one computed for pure rotation. Buoyancy-driven convection only slightly enhances this motion; compare Figs. 2 and 8. However, the turbulent viscosity for the combine rotation and buoyancy-driven #ow is

86

A. Lipchin, R.A. Brown / Journal of Crystal Growth 205 (1999) 71}91

Fig. 9. Contours of the turbulent viscosity for the (a) JL, (b) TL and (c) WF turbulence models and the (d) azimuthal velocity for the JL model calculated for buoyancy-driven convection and rotation with X "20 rpm and X "!8 rpm.  

larger by an order-of-magnitude than the one for the pure rotation driven #ow. A possible perspective on this result is to consider the balance of di!erent forces in the momentum equation. The ratio of viscous to Coriolis forces is evaluated in terms of the Ekman number as Ek,l/2XR. If the magnitude of the turbulence viscosity is approximately l +2;10l , X"8 rpm and R"R , then 

 Ek+10\; so that that the #ow is determined by a balance between Coriolis, inertia and pressure forces in Eq. (1). Buoyancy forces a!ect the #ow directly through the Eq. (2) and indirectly by generating the turbulence energy in Eqs. (6) and (7). This e!ect is taken into account through the eddy viscosity and does not directly alter the bulk #ow. Viscous forces are important only in Ekman layers near boundaries where the turbulent viscosity is small. To evaluate the balance of the buoyancy and inertial forces in Eq. (2) consider the Burgers number Bu,gb*¹/4XR, i.e. the ratio of buoyancy to Coriolis forces. The Coriolis force is of the same order-of-magnitude as the inertial force because these terms balance in Eq. (1). For the parameter values listed above, Bu+10\ and the buoyancy

force should not play a signi"cant role in the determining the #ow for these rotation rates. The #ows for higher rotation rates are shown in Figs. 10 and 11 and characteristic values are listed in Table 3. The results for pure rotation (*¹"0) also are reported in Table 3 for comparison. The value of the turbulent viscosity is less than the result for more moderate rotation rates, while the motion intensity and turbulent energy are larger. For high rotation rates buoyancy plays an even smaller role than for moderate rotation rates and the motion is almost completely determined by the Coriolis force. The distribution of turbulence viscosity for high rotation rates (Fig. 11a) is qualitatively similar to the results for pure rotation with either moderate (Fig. 3) or high (Fig. 11d) rotation rates. The di!erences between the temperature distributions for moderate and high rotation rates are caused by the di!erences in the turbulent viscosity distribution. For high rotation rates there is a region of low l near the center line under the crystal  surface. The changes in the heat #uxes through the surfaces of the melt with rotation rate are listed as

A. Lipchin, R.A. Brown / Journal of Crystal Growth 205 (1999) 71}91

87

Fig. 10. Contours of the stream function for the (a) JL, (b) TL and (c) WF turbulence models computed for buoyancy-driven convection and rotation with X "100 rpm and X "!20 rpm. The contours computed with the JL model without buoyancy-driven   convection, but with the same rotation rates, are shown in (d).

Fig. 11. Contours of the (a) turbulent viscosity, (b) temperature, and (c) azimuthal velocity computed with the JL model with buoyancy-driven convection and rotation with X "100 rpm and X "!20 rpm. The contours of the turbulent viscosity computed   with the JL model without buoyancy-driven convection, but with the same rotation rates, are shown in (d).

88

A. Lipchin, R.A. Brown / Journal of Crystal Growth 205 (1999) 71}91

Nusselt numbers in Table 3 and are connected with corresponding changes in the intensity of motion in the (r, z) plane. 4.4. Oxygen concentration The transport of oxygen from the dissoluting crucible to the crystal is determined by di!usion, by convection driven by the mean velocity "eld and by turbulent transport modeled by the turbulent diffusivity. The oxygen concentration distribution was computed using the velocity and turbulent viscosity "elds. Comparison of the oxygen distribution predicted for the three turbulent models is another means of illustrating the di!erences in the predictions. The mesh used in the calculations with the WF model did not allow calculation of a solution for the oxygen concentration because the size of the mesh near to the walls is too large to describe the steep concentration boundary layers there. Wall functions on the concentration have to be introduced to avoid this problem. The oxygen distributions computed with the JL and TL models are shown in Fig. 12. In both cases, there are boundary layers near the crucible surface that arise because of the high Schmidt number. Near the wall the e!ective di!usion coe$cient, a sum of molecular di!usion (1/Sc) and turbulent di!usion (l /p ), drops by several orders-of-magni 2 tude. This e!ect does not occur in the calculation of the temperature distribution because the low Prandtl number makes conduction dominant near the wall. Large discrepancies exist between the predictions of the JL and TL models for nonrotational #ow and moderate rotation rates. These di!erences are largest in the near-wall regions where turbulent transport plays a signi"cant role, because Sc<1. For high rotation rates the #ow is highly turbulent (l +10) and the di!erences in the prediction for  the two model are smallest. The Sherwood numbers (the mass #ux normalized by the value for pure di!usion) are reported in Table 3 for all surfaces. The Sherwood numbers increase as the rotation rates are increased from no rotation to moderate rotation rate. A similar increase was observed in the experiments mentioned in Ref. [30]. For high rotation rates the Sherwood numbers are less than

Fig. 12. Contours of oxygen concentration computed with the JL and TL turbulence models for (a) (X "0, X "0),   (b) (X "20 rpm, X "!8 rpm), and (c) (X "100 rpm,    X "!20 rpm). Buoyancy-driven convection is included in the  calculations.

the values for moderate rotation rates. This nonmonotonic behavior is opposite to the prediction for the Nusselt numbers and is explained by the di!erence in the Prandtl and Schmidt numbers; Pr;1 and Sc<1. The computed radial distribution of oxygen in the crystal is shown in Fig. 13. Experiment data is available for the oxygen distribution; see Ref. [30]. Unfortunately, the crudeness of the computational model makes comparison useless. The concentration value on the center line at the crystal surface is reported in Table 3. There are a large discrepancy between radial concentration distributions in the crystal predicted by the JL and TL models with and without moderate rotation. Again, for high rotation rates the predictions are closer.

A. Lipchin, R.A. Brown / Journal of Crystal Growth 205 (1999) 71}91

Fig. 13. Pro"les of oxygen concentration in the crystal computed radially along the melt/crystal interface with the JL and TL models for the rotation rates (a) (X "0, X "0),   (b) (X "20 rpm, X "!8 rpm), and (c) (X "100 rpm,    X "!20 rpm). Buoyancy-driven convection is included in the  calculations.

5. Conclusions Comparison of the predictions for the three k}e turbulence models for the variety of the #ow conditions considered here leads to some general conclusions. First, the discrepancies in the predictions of the three models become smaller with increasing rotation rates. For intensive enough #ow driven by rotation, CZ bulk #ow is governed by the balance of inertial and Coriolis forces, which is independent of the turbulence model. For high rotation rates, buoyancy force does not play a signi"cant role, and the #ow becomes highly turbulent and more tractable for all three models. However, near the boundaries viscous force is important and the turbulence is weak. The heat and oxygen #uxes across surfaces change di!erently with rotation rates, because of

89

the di!erence in Pr and Sc. The conditions Pr;1 and Sc<1 cause the oxygen transport to be more sensitive to the the turbulent viscosity than is the heat transport. The results of these calculations put us in a position to recommend using the low Reynolds number model of Jones}Launder (JL) preferentially over the other two k}e models for simulating CZ bulk #ow. This recommendation is based in part on the intuitive connection between the JL results and the predictions of laminar #ow calculation for condition where the #ow is weak. As shown above, the JL model predicts #ow states that approach the laminar #ow with reducing intensity while neither the TL or WF model have this qualitative feature. The shortcomings of the TL and WF models are coupled to the intrinsic problem of the representations of the boundary conditions along solid boundaries in both the wall function and two-layer models. For weakly turbulent #ows or near turning or stagnation regions of the #ow the viscous sublayer become large and using a uniform size for either wall function in the WF model or the wall layer in the TL model leads to poor solution quality adjacent to the wall. In the case of the predictions of the oxygen concentration the situation is compounded by Sc<1 which forces the use of wall function so close to the wall that the calculation fails altogether. The validity of the JL model for quantitative prediction of heat and species transport in the CZ system cannot be assessed from these results, but must await quantitative comparison with experiment in actual CZ system. This requires the coupling of the #ow calculation with a thermal-capillary heat transport model. Following Kinney et al. [30], we have accomplished this for both a generalization of the "nite volume simulation of melt convection described in Ref. [51] and a "niteelement-based simulation of convection [52]; both implementations use the "nite element heat transfer analysis of Bornside et al. [7]. These analysis and careful experiment, such as those of Seidl et al. [16], will lead to quantitative assessment of these predictions and, hopefully, to simulations of engineering utility in design of CZ #ow system.

90

A. Lipchin, R.A. Brown / Journal of Crystal Growth 205 (1999) 71}91

Acknowledgements We thank MEMC Corporation, Sematech, Shin-Etsu Hanotai Company, and Wacker Siltronic for "nancial support of this research.

References [1] K.M. Kim, Solid State Technol. 39 (11) (1996) 70. [2] T. Abe, Silicon materials science and technology, Electrochemical Society Proceedings, Vol. 98-1, 1998, p. 157. [3] J.J. Derby, R.A. Brown, J. Crystal Growth 74 (1986) 605. [4] J.J. Derby, R.A. Brown, J. Crystal Growth 75 (1986) 227. [5] A.B. Crowley, IMA J. Appl. Math. 30 (1983) 173. [6] F. Dupret, Y. Ryckmans, P. Wouters, M.J. Crochet, J. Crystal Growth 79 (1986) 84. [7] D.E. Bornside, T.A. Kinney, R.A. Brown, G. Kim, Int. J. Numer. Meth. Eng. 30 (1990) 133. [8] R.A. Brown, A.I.Ch.E. J. 34 (1988) 881. [9] F. Dupret, N. Van den Bogaert, in: D.T.J. Hurle (Ed.), Handbook of Crystal Growth, Vol. 2B, Modelling Bridgman and Czochralski Growth, North-Holland, Amsterdam, 1994 (Chapter 15). [10] P.A. Sackinger, R.A. Brown, J.J. Derby, Int. J. Numer. Meth. Fluids 9 (1989) 453. [11] J.R. Ristorcelli, J.L. Lumley, J. Crystal Growth 116 (1992) 447. [12] T.A. Kinney, D.E. Bornside, R.A. Brown, K.M. Kim, J. Crystal Growth 126 (1993) 413. [13] K.M. Kim, W.E. Langlois, J. Electrochem. Soc. 138 (1991) 1850. [14] K.M. Kim, P. Smetana, J. Appl. Phys. 58 (1985) 2731. [15] K.M. Kim, P. Smetana, J. Electrochem. Soc. 133 (1986) 1682. [16] A. Seidl, G. MuK ller, E. Dornberger, E. Tomzig, B. Rexer, W. von Ammon, Silicon materials science and technology, Electrochemical Society Proceedings, Vol. 98-1, 1998, p. 417. [17] H. Tennekes, J.L. Lumley, A First Course in Turbulence, MIT Press, Cambridge, MA, 1974. [18] D.C. Wilcox, Turbulence Modeling for CFD, DCW Industries, La Cana da, CA, 1994. [19] J.J. Derby, R.A. Brown, J. Crystal Growth 87 (1987) 251. [20] H. Iacovides, B.E. Launder, Int. J. Heat Fluid Flow 16 (1995) 454. [21] H. Iacovides, B.E. Launder, H.-Y. Li, Int. J. Heat Fluid Flow 17 (1996) 22. [22] K. HanjalicH , S. VasicH , Int. J. Heat Mass Transfer 36 (14) (1993) 3603. [23] V.C. Patel, W. Rodi, G. Scheuerer, AIAA J. 23 (9) (1985) 1308. [24] W.P. Jones, B.E. Launder, Int. J. Heat Mass Transfer 15 (1972) 301.

[25] R.A.W.M. Henkes, C.J. Hoogendoorn, in: R.A.W.M. Henkes, C.J. Hoogendoorn (Eds.), Turbulent Natural Convection in Enclosures. A Computational and Experimental Benchmark Study, EETI, Paris, 1993, p. 64. [26] T. Zhang, F. Ladeinde, H. Zhang, V. Prasad, Proceedings of the 1996 31st ASME National Heat Transfer Conference. Part 1. August 3-6 1996, Vol. 323, no. 1, p. 17. [27] S. Kobayashi, J. Crystal Growth 99 (1990) 692. [28] S. Kobayashi, S. Miyahara, T. Fujiwara, T. Kubo, H. Fujiwara, J. Crystal Growth 109 (1991) 149. [29] H.-T. Chung, S.-Ch. Lee, J. Yoon, J. Crystal Growth 163 (1996) 249. [30] T.A. Kinney, R.A. Brown, J. Crystal Growth 132 (1993) 551. [31] B.E. Launder, B.I. Sharma, Lett. Heat Mass Transfer 1 (1974) 131. [32] J.R. Ristorcelli, J.L. Lumley, J. Crystal Growth 129 (1993) 249. [33] R.A.W.M. Henkes, C.J. Hoogendoorn, Int. J. Heat Mass Transfer 32 (1) (1989) 157. [34] J.H. Ferziger, M. PericH , Computational Methods for Fluid Dynamics, Springer, Berlin, 1996. [35] W.E. Langlois, Ann. Rev. Fluid Mech. 17 (1985) 191. [36] Q. Xiao, J.J. Derby, J. Crystal Growth 129 (1993) 593. [37] A. Anselmo, V. Prasad, J. Koziol, K.P. Gupta, J. Crystal Growth 134 (1993) 116. [38] A.A. Wheeler, J. Crystal Growth 102 (1990) 691. [39] W. Rodi, Turbulence Models and their Application in Hydraulics, a State of the Art Review, International Association for Hydraulic Research, Delft, The Netherlands, 1980. [40] E. Nobile, in: R.A.W.M. Henkes, C.J. Hoogendoorn (Eds.), Turbulent Natural Convection in Enclosures, a Computational and Experimental Benchmark Study, EETI, Paris, 1993, p. 214. [41] K. HanjalicH , Proceedings of 10th International Heat Transfer Conference, Brighton, UK, 1994, Vol. 1, p. 1. [42] R.A.W.M. Henkes, F.F. Van Der Vlugt, C.J. Hoogendoorn, Int. J. Heat Mass Transfer 34 (2) (1991) 377. [43] R.A.W.M. Henkes, C.J. Hoogendoorn, in: R.A.W.M. Henkes, C.J. Hoogendoorn (Eds.), Turbulent Natural Convection in Enclosures, a Computational and Experimental Benchmark Study, EETI, Paris, 1993, p. 185. [44] R.A. Kuyper, Th.H. Van der Meer, C.J. Hoogendoorn, R.A.W.M. Henkes, Int. J. Heat Mass Transfer 36 (11) (1993) 2899. [45] R.A.W.M. Henkes, C.J. Hoogendoorn, Numer. Heat Transfer Pt. B 28 (1995) 59. [46] N.Z. Ince, B.E. Launder, Int. J. Heat Mass Transfer 10 (2) (1989) 110. [47] S. VasicH , K. HanjalicH , in: R.A.W.M. Henkes, C.J. Hoogendoorn (Eds.), Turbulent Natural Convection in Enclosures, a Computational and Experimental Benchmark Study, EETI, Paris, 1993, p. 133.

A. Lipchin, R.A. Brown / Journal of Crystal Growth 205 (1999) 71}91 [48] M. Wolfshtein, Int. J. Heat Mass Transfer 12 (1969) 301. [49] C.-J. Kim, S.T. Ro, Numer. Heat Transfer Pt. B 28 (1995) 385. [50] T. Mori, R.A. Brown, private communication, 1996. [51] A. Lipchin, R.A. Brown, Computational technology for #uid/thermal/chemical system and industrial applications,

91

Proceedings of the 1998 Joint ASME/JSME Pressure Vessels and Piping Conference, July 26}30, 1998, San Diego, CA. [52] T. Mori, T.R. Sinno, R.A. Brown, K.M. Kim, Abstracts of the 193rd Meeting of The Electrochemical Society, May 3}8, 1998, San Diego, CA, p. 496.