Heat transfer, fluid flow, and interface shapes in the floating-zone growth of tube crystals

NH CRYSTAL GROWTH ELSEVIER

Journal of Crystal Growth 141 (1994) 265—278

Heat transfer, fluid flow, and interface shapes in the floating-zone growth of tube crystals C.W. Lan Chemical Engineering Department, National Central Universily, Chung Li 32054, Taiwan, ROC Received 10 November 1993; manuscript received in final form 9 March 1994

Abstract A computer simulation is conducted to study crystal growth of tubes by the floating-zone (FZ) method. Because the inner and outer radii of the steadily growing tube are unknown a priori and are coupled with fluid flow and heat transfer in the melt, solid/melt interfaces, and inner and outer free surfaces, the modeling of this process becomes rather difficult. A robust computer model based on the body-fitted coordinate finite-volume method (BFCFVM) and Newton’s method with a preconditioned iterative linear equation solver is used to study the process. The effects of fluid flow and heat transfer on the solid/melt interfaces, free surfaces, and the dimension of the steadily grown tube crystal are demonstrated through computer simulations for both a high Prandtl-number material, i.e., NaNO 3, and a low Prandtl-number material, i.e., silicon. The convection induced by both thermocapillary and buoyancy forces is considered. The effects of key process parameters are also discussed, including feeding/pulling speeds and the internally applied pressure. 1. Introduction

ing tube crystals has not been widely used [1,61.

Floating-zone (FZ) crystal growth is unique in that it is a contamination-free process. Crystals with high purity can be grown or zone-refined by this process [1]. There are some special interests in the growth of crystal tubes, e.g., silicon and silicon carbide for high temperature reactors or furnaces employed in the fabrication of integrated circuit [21. Most of these tube crystals are grown by meniscus-defined processes, e.g. edgedefined film-fed growth (EFG), Stepanov process [3,41and crystal pulling [51.However, for the growth of high purity crystals, these methods may not be appropriate. The FZ process may be the best candidate to meet the high purity requirement. On the other hand, due to poor understanding of the process, the FZ process for grow-

Furthermore, thermal-capillary interaction, i.e., the interaction of heat transfer and free surfaces (and thus the size of the grown crystal) in tube system is complicated and poorly understood. In FZ crystal growth under normal gravity, convection in the melt zone can be induced by buoyancy forces, i.e., natural convection, and by surface-tension gradients on the free surface, i.e., thermocapillaiy convection. Convection in the melt can affect the melt zone length and hence the stability of the melt zone. It can also affect the shape of the melt/crystal interface and hence the quality of the resultant crystal. In FZ tube growth, the dimension of the grown tube crystal is unknown a priori and is determined by thermalcapillary interaction, which in turn is affected by the convection in the melt. The present study is

0022-0248/94/$07.00 © 1994 Elsevier Science B.V. All rights reserved SSDI 0022-0248(94)00101-0

C. W. Lan /Journal of Crystal Growth 141 (1994) 265—2 78

266

carried out mainly to help evaluate the effects of convection in the melt zone on the zone length, the shape of the melt/crystal interface, and the dimension of the steady grown crystal in the FZ tube crystal growth system. Since the extent of these effects can be affected significantly by the magnitude of the Prandtl number of the material being grown, both a material with a low Prandtl number (i.e., Si, Pr 0.01) and a material with a high Prandtl number (i.e., NaNO3, Pr 9.12) are =

=

considered in this report. Numerous computational studies, e.g., Refs. [7—12],have been carried out to investigate heat transfer and fluid flow in FZ growth of rod crystals. Especially, in the studies by Lan and Kou [11,121, fluid flow, heat and mass transfer, and interface shapes (non-isothermal melt/solid interfaces) can all be considered. However, no computational studies were carried out for FZ tube growth until recently. Lan [131 used a finitevolume/ Newton method to study the conduction-dominated heat transfer and the thermalcapillary interaction in the FZ tube growth. Thermal fields, melt/solid interfaces, melt/gas free surfaces, and the radii of the steadily grown tube crystal can all be predicted through the model, However, convection in the melt zone was not considered. In this report, convection in the melt under normal gravity will be considered, so that the model can be used for both conductiondominated (or for low Pr) and convectiondominated (or for high Pr) cases. In FZ rod growth, the size of the steadily growth crystal can be predetermined according to an overall mass balance. However, for FZ tube growth the inner and outer radii of the steadily growing crystal are unknown a priori and are strongly coupled with thermal fields and free surface shapes [13]. This, therefore, complicates the modeling of tube growth. Further, one more free surface, i.e., the inner melt/gas interface, needs to be considered in the tube growth than in the rod growth. Thus the modeling of tube growth is more difficult and relies on a more robust numerical method. Lan and Kou [11,12] have developed a numerical procedure based on the body-fitted coordinate finite-difference method with Gauss—Seidel iterations for computing the

fluid flow, heat and mass transfer, and interface shapes of FZ growth of rod crystals. During the iterations, the updates of the interfaces and the free surface are through outer iteration loops, i.e., the so-called decoupled approach. Their code has proven to be very versatile and useful in rod growth simulation [7,8,11,121. However, because of the much stronger interaction of thermal fields, interfaces, free surfaces, and the size of the grown crystal in the tube system, the simple decoupled iteration scheme used in FZ rod growth fails to converge, even for the cases without convection. Recently, this numerical procedure has been further improved [14] and has been applied, successfully, to the modeling of FZ tube crystal growth of low Prandtl-number materials, where the heat transport is dominated by conduction [13]. In this new approach, a global Newton iteration scheme that updates all variables simultaneously is implemented. In addition, to achieve a higher global conservation of energy, governing equations are written in conservative-law form, and discretized by the body-fitted coordinate finite-volume method (BFCFVM). Nearly quadratic convergence of this approach for the free surface, interfaces, and thermal and flow fields is observed [14]. Furthermore, the Newton’s method is optimal for computer-aided calculation of stability [15]. In this study, the new approach will be adopted to study the effects of convection on the FZ growth of tube crystals. The formulation of the model is described in the next section and the BFCFVM/ Newton scheme is presented in section 3. After the method of solution, the accuracy of the method is examined first in section 4, and then the effects of convections and some key process parameters on the growth of NaNO3 (Pr 9.12) and Si (Pr 0.01) are illustrated through calculated results. Finally, brief conclusions are drawn in section 5. =

=

2. Model description The steady-state FZ growth of tube crystals is simulated using a convective heat transfer model. The heating source considered here is a resistance heater, so that no electromagnetic force

C. W. Lan /Journal of Crystal Growth 141 (1994) 265—2 78

Feed Tube 1D

Rfo

i

feed and crystal are the same. Yet, the radii of the growing tube are both unknown a priori and coupled with capillary-thermal/flow field interaction. When a steady-state is achieved, the lengths of the feed and the crystal tube are usually many times their radii. In order to reduce computational load, while maintain accuracy at the same time, the system is divided into three regions. As illustrated in Fig. 1, there are an outer region in the crystal, an outer region in the feed, and an inner region which consists of the melt zone and the solid material near solid/melt interfaces. The inner region is made long enough so that the two-dimensional temperature distribution in this region becomes essentially one-dimensional (i.e.,

~

uniform in the radical direction) near its two

~

H—Rfi

I

1D zz

2

D1

~fi

aDf

fl

aDmj (Meniscus)

D

melt

m

0

ôD~ OD,~1

ZO 1D

C rySta1 sketch

Fig. 1. Schematic

1D —R~~

of floating-zone tube growth system.

and eddy-current heating due to induction coils need to be considered. If the heating is axisymmetric, the physical domain for feed, melt, and crystal tubes can be taken as shown in Fig. 1. As it is, it can be treated as a two-dimensional model; the flow instability, which may lead to three-dimensional convection, is not considered here. The domain consists of the feed tube (Df) being fed from above, the melt (Dm) being held by surface tension, and the steadily growing tube crystal (Do) that is being pulled downward from below. Each of these regions is characterized by a set of physical of the the subscript boundaryj have beenproperties. labeled asSegments ~ where corresponds to the parts shown in Fig. 1. The flow and temperature fields as well as the shapes of the feed front (hf(r)), the growth front (h~(r)), and the inner melt/gas (R 1(z)) and outer melt/gas (Rjz)) free surfaces are represented in cylindrical coordinate system (r,z). For a feed tube of inner radius Rf~and outer radius R~0at a feeding speed Uf, a steady-state growth can be obtained with the crystal tube being pulled at a speed U~with an inner radius ~ and an outer radius R~0,where (R~0 R~~)UC (R~ R~)U~,assuming that the densities of —

=

—

267

ends. Ascansuch, heat transfer in the two outer regions be treated as one-dimensional. The same strategy has been applied successfully in the numerical simulation of FZ rod growth [7,8,11,12]. In the melt (Dm), the governing equations in conservative-law form (or the so-called divergence form) for fluid flow and heat transfer can be described in terms of stream function ~i, vorticity w, and temperature T as follows: Equation of motion

a

(0)

~

ar

r

az

+

a

~

az

r

—

ar

a

i

a

—

—

—(~arw) +

r

ar

a

1

—

—

a —(ttrw)

az r az

0 PLIequation Stream 3S’~

1

— —

—

+

~

-~-

az pLr az

—~-— ~

+

ar pLr ar

—

~

0

2

—

Energy equation -~.-

ar

~

—

az

~°

a +

—

az

/

C T—~

az aT

i rk— ~,

—

az

~°

a +

—

ar

ar

a~ rk—

=

0.

(3)

ar

The energy equation is used to calculate T in both the melt phase (Dm) and the solid phase

268

C. W. Lan /Journal of Crystal Growth 141 (1994) 265—2 78

(Df ~), where k

kL and C~ CPL in the melt k5 and C~ ~ in the solid phase. kL and k5 denote the thermal conductivities of melt and solid, respectively, and CPL and C,~,5the specific heats of melt and solid, respectively, PL is the melt density, and f3 the thermal expansion coefficient. In Eq. (1), p. is the melt viscosity, which can be temperature dependent, but higher order terms related to the derivatives of p. are neglected here [16]. The stream function ~i and vorticity w in the above equations are defined in terms of radial (u) and axial (v) velocities as: phase and k

=

—

~—

w

=

az

v

~,

az

pLr

au

=

1

a~

(4)

—;

=

pLr

ar

av —

(5)

~.

ar

The thermal boundary conditions are as follows: (1) At the solid/melt interfaces (aDCf), the interfacial energy balance is k(n

‘

(4) Far away from the melt zone,

=

atp

1 u

=

=

VT) Is

—

k(n

VT)

IL +psi~HUfCn

.

0, (6) =

where n is the unit normal vector pointing into the solid. Uf~denotes U~and U~at the feed and growth fronts, respectively, Ps is the solid density, and ~H the heat of fusion. L and S mean the melt and solid sides, respectively. The temperature at the melt/solid interfaces (aDf~) is set to an equilibrium melting temperature of the material: TIf=TIC=Tm.

(7)

(2) Heat transfer from the system to the ambient is by both radiation and convection according to the energy balance along the material surface

(aDfo mu —k(n. VT) =h(T— 1~)+cr(T4— T~).

(8)

(3) Along the inner surface of the material

T= (Ti/2)f, (10) where (T1/2)f C is the temperature in the feed tube at z z2 or in the crystal at z 0 (see Fig. 1). These two points are located so far away from the melt zone that heat transfer becomes one-dimensional beyond them. The implementation of this boundary condition has been described elsewhere [17]. The heat exchange between the surrounding (heater and ambient) and the material is dictated =

=

by an effective ambient temperature Ta(Z) specified along the zone length. The effective ambient temperature distribution due to the heater, as in previous studies [11—13],is assumed Gaussian as follows: 2}+ i~, Ta(Z) (T~ Tax) exp{_ [(z -~z2)/a] (11) where T~ and Tao. are the peak and background =

—

—

temperatures, respectively, and the parameter a is the width (the standard deviation) of distribution. The complications of computing the view factors from points along the surface of the material to the ambient are neglected; the effect of the view factors is simply lumped into the effective ambient temperature Ta(Z). The fluid-flow boundary conditions are as follows: (1) Along the surfaces of the melt zone (aDmimo), fi=0

aDm,

on

(12) on

au ~

with

at

ay

—(s.VT),

t .

(13)

—

=

ns

aDmo,

(14)

=

(aDtjmicj), —k(n’ VT) 0, (9) where an adiabatic condition is assumed. Heat loss by both radiation and convection inside the tube can be considered if necessary. =

where n and s are respectively the unit normal vector and unit tangential vector at the free surfaces, t is the stress tensor, and ay/aT the surface-tension—temperature coefficient of the melt. Eq. (14) represents the shear stress balance at the

C.W Lan /Journal of Crystal Growth 141 (1994) 265—2 78

free surfaces. The stream function is set to zero at the inner free surface as a reference. (2) At the melt/solid interfaces (aDf ~),

au 2

~psU~,~(r

=

—

av

R~~ 1),

=

~—

—

-s—,

with

=

—

which is used as the equation for R~0.

u=0,

t)UfCps/pL,

(16)

where RfjC~denotes Rf1 at the melt/feed interface and R~1at the melt/crystal interface. The shapes of inner and outer free surfaces are governed by the normal stress balance. The normal-stress balance at the free surfaces can be written as: nfl: t

is also required. Eq. (19) is used for the unknown inner radius ~ of the grown tube. Finally, the overall mass balance also needs to be satisfied for the steady-state case, 2 R~,)U~ (R~ (R 0 R~)Uf, (20) —

(15)

1 =

~‘

1 +

~—

(17)

+

1

2

where R1ambient and R2pressure. are the radii of curvature, is the Outside the tube ~aand ~ set to zero as a reference, while inside the tube 1~a is an applied overpressure (i.e., the pressure difference between inside and outside of the tube). It should be pointed out that the pressure in the total normal stress (nn: t) needs to be integrated along the free surfaces starting from the melt/crystal/gas triple junctions. The pressure at the outer melt/crystal/gas triple junction 1~a

3. Finite volume analysis 3.1. Coordinate transformation

The above governing equations and boundary conditions are transformed into those in terms of general (nonorthogonal) curvilinear coordinates ~ ~) which fit all the interfaces and free surfaces, as shown in Fig. 2. treated In this accurately, way, all the boundary conditions can be as already described previously [7,8,11—14,17]. By following the procedure of Thompson et al.

I.LJ

I

0, however, is an unknown and is determined by specifying a growth angle constraint [18] as

dR0/dz

I z=h0(R~)

=

tan

~.

uZ

/ ‘~

z~h (R

) —

—

Lan

11 fl\ ~IY)

tII

-;

1

i

-

.

(18)

The growth angle ~ is approximately 2.5° for NaNO3 [7] and 11° for silicon growing on the <111) surface [18] This growth angle constraint is to ensure the crystal to grow upwards with a constant outer radius. The pressure at the inner triple junction is obtained by integrating the equation of motion along the melt/crystal interface from the outer triple junction (P0). Clearly, any setting of ~a to the outer free surface only affects the P0 value, and does not affect the calculated results. Again, to have a constant radius during crystal growth, a growth angleinner constraint for the inner melt/crystal/gas triple junetion, R

1

C—-

i

I

P

~‘ u

269

.~

:1.]

~

I .

I

I

.i~J~:’. ::. I

f.

.

.1

—

—~ I

i

(a)

I

-

______________ .

I

I

.

“ ..

I

~—.

. —-

I

.

-

1’tII~.

(b) o

.

—

:

:1

-

‘ I

~

Nx Physical Domain Compiilational [)omain Fig. 2. A nonorthogonal body-fitted grid system used for computation: (a) mark physical (b) interfaces. computational domain. The arrowheads the domain; melt/solid

C. W. Lan /Journal of Crystal Growth 141 (1994) 265—2 78

270

(2) At the outer surface of the feed tube, melt, and crystal ~ Nt),

Table 1 Coefficients a, b, c and d in Eq. (21)

~

=

a

b

c

d

0

l/pLr

1

~o

g22T~ g12~ 222 Jgi/

T w

C1/r l/r ~r —pp~gaT/ar 0 kr 1 0 ___________________________________________

Ta)+Eff(T4~T~).

(24) (3) At the inner surface of the feed tube, melt,

[19], Eqs. (1) through (3) can be transformed into the following general conservative-law form: a / ~ \ a / a~\

—la~—I-—la~—I a~k a~ a~\ a~

and crystal

e~1

=

0),

2 0. k(g22T,7 g12Tf)/Jg~ (4) Far away from the melt zone,

(25)

=

—

T= (Ti/

(26)

2)1 ~

a lb / a(c~) a(c4) g12 a~ j a~ a~ +—I—~g22 a [b I a(c4) a(c4) ~ —I a~ a~ a~ +Jd=o,

The transformed fluid-flow boundary condi-

—

tions are as follows: (1) Along the free surfaces of the melt zone

)j

.~

(~=0,

(21)

Ni), 1

a7 a~

2

+

w=+

(

az au

ar ac

~ where Iar\

gii=(—j

2

/

2

‘az~ +I—~ ~ g

ar

ar ar

az az

g12=——+———,

(az\ 2

2

(~)

22= ar az

+

where the

+

sign is for

i~=

ar az

=

(2) At the melt/solid interfaces ~ 2 —~ (28) ~psU~(r g 1’ar )2 +r—I~. a2r 11 (a2~ ~J2 ~ PSUf.C[~ ~ a~2~ (29) =

=

—

\

U

—

Along the free surface, i/i is constant and a~p~a~ 0. From this, it can be shown that n (UCr + vet) 0. This equation is the so-called kinematic condition. The normal stress balance equation, Eq. (17), is applied to both inner and outer free surfaces. The transformed form for inner and outer free surfaces is as follows: =

=

az d2R

a2z dR

g 22 az

g11T~—g12T~ Jg~/2 L

g~u/2

=

0~

(27) sign for

=

=

+

—

,

J=——————. a~a~ a~a~ a~a~equaThe coefficients a,b,c and da~a~ in the above tion are given in Table 1 for ~ ~, w, and T, respectively. Mathematically, the conservative-law form is exactly the same as the non-conservative one used in the previous papers [7,8,11,12]. However, the conservative-law form is easier to be implemented and less sensitive to grid roughness, because no second order derivatives of coordinates are needed in computation. The transformed thermal boundary conditions are as follows: (1) At the melt/solid interfaces (4 ~a,i,’ where ~a and ~b are explained in Fig. 2), T= Tm, (22)

~~11Tf_~12T~

0 and the

(23)

S

2=0, (30) — for the —total —g~ where S —is the term normal stress (including the unknown outer triple junction pressure P 3a’ static pressure due to gravity, 0,and overpressure dynamic pressure ‘ due to fluid

C.W. Lan /Journal of Crystal Growth 141 (1994) 265—2 78

flow) [7,8,11,12]. The boundary conditions for above the equation are R I T7=N0,f=f~ R~0, (31) =

RIN~,©..~b=Rfo,

(32)

RI ~ RI ~~=O,f=~b

(33)

=

=

R~1, Rf1,

(34)

and the growth angle constraints, Eqs. (18) and (19), can be transformed into aR ‘az

~

—‘I

where ‘e’ for example, represents the total flux of ~ across the face e. Eq. (37) involves no approximation, and represents the finite-volume equations of the differential equations. The numerical evaluation of the terms in Eq. (37) requires the calculation of geometrical factors for control volumes and an interpolation scheme for calculating quantities at the cell faces from their adjacent nodal values. The details of the discretization scheme is discussed elsewhere [14]. Eq. (37) is applied to all variables in the do-

±tan4~, (35) where the + sign is for the outer triple junction and the — sign for the inner one. Finally, the overall mass balance condition, Eq. (20), is the equation for R~. For convenience, ~1i~ z.1~ 1 is used in the computational domain, while N~+1 is the number of nodal points in the ~ direction. L10,Aç;

=

271

=

=

main. (30)—(35), For the normal stress balance equations, Eqs. a second order finite-difference discretization is used. Trapezoidal method is also adopted to integrate the pressure along boundaries. After the discretizations for both governing and boundary equations, a set of nonlinear equations can be obtained, f(x) =f(T, ~i, ~, R~), R0(~),h~(ij), h~(i~), P0)=0;

(38)

3.2. Method of solution

The equation of motion, the stream equation, and the energy conservation equation (Eqs. (1) (3)) are discretized by employing the finite-volume approach. The physical domain, in (r, z), is subdivided into a finite number of contiguous volumes (CVs) of volume V, which are bounded by cell faces located about halfway between consecutive nodal points. This domain can be transformed into the computation domain (~,~) with CV of volume V’, and rdV= rJdV’. All variables are stored at nodal points. Now, the integration of these equations over the CV can be proceeded as follows:

fv[Eqs. (1) to (3)] r dV f [Eq. (21)] di1 d~ 0. =

(36)

=

f (Jd) d~d~

=

V

where the correction vector ~5n+i is the solution of the linear equation set j~n+i —f(x~). (40) =

The components of the Jacobian matrix j, which are formed by explicit differentiation as J~

af1/ax1, represent the sensitivity of the residual

After the Gauss theorem is applied, the above equation can be transformed into surface integrals (fluxes) over the surface of CV. The resulting balance equations for each CV can then be expressed as le — 1w + I~— I, +

the unknowns R~1 and R~0are equal to Ri(~a) and Ro(~a), respectively. The nonlinear equation set is then solved by Newton’s method for all variables simultaneously. Starting from an initial guess to the vector of unknowns, x°,successive updates are constructed as ~ =..v~~ (39)

0,

(37)

vector to the solution vector, and are obtained by finite-difference approximations [201. The dependence of the residual equations on the unknowns leads J to a “fish-bone” structure, which is typical in free-boundary problems [21]. Eq. (40) is solved by the GMRES iterative linear equation solver with a preconditioner of the incomplete LU decomposition without fill-in, ILU(0) [22]. The ordering of variables and the choice of preconditioners are crucial to the success of iterative matrix methods here, because some diagonal ele-

272

C.W. Lan/Journal of Crystal Growth 141 (1994) 265—2 78

ments of j are zeros. In the implementation of Newton’s scheme, continuation can play an im-

Table 2

portant role, particularly when the solution behavior changes dramatically with small changes in parameters, pseudo continuation [23] isHence, used tothe locate thearc-length results from one to the other. Again, the detailed description of the numerical implementation is described elsewhere [14].

306.8°C z.1H = 182 J g~’

4. Results and discussion Due to the high Prandtl number of NaNO 3 (Pr 9.12), the effect of convection on heat transfer is very significant in the FZ process [7,8]. Therefore, the numerical results for NaNO 3 may be vulnerable to grid systems used in computation. The effect of the discretization on the accuracy of calculations is examined through a Sequence of mesh refinements. The meshes, which are similar to Fig. 2a, given in the order of feed, melt, and crystal, are 17 x 21, 17 x 21, and 17 x 21 for mesh I, 21 x 21, 21 x 31, and 21 x 21 for mesh II, and 27 X 21, 27 x 41, and 27 x 21 for mesh III. The physical properties of NaNO3 are listed in Table 2, and the input parameters are: R~1 0.7 cm, Rf0 0.9 cm, T0 612°C,Tao. 25°C, and a 0.3 cm. These parameters will also be used throughout this report for NaNO3 unless otherwise stated. The calculated results based on the three different meshes for a case of upward growth, while the feeding and pulling speeds are kept the same (Uf= U~, —2 X i0’~ cm/s), are summarized in Table 3. As shown, the calculated results based on meshes II and III are significantly close to each other. In order to have a reasonable balance between accuracy and CPU time, mesh II will be used for all calculations in =

=

=

=

=

Physical properties of NaNO

3 [24—26]

Tm

1

~ Wcm’ °C kL=5.65X10_3+44.7(T_Tm)X10_7Wcm_I0C_I

Cpç = CPL= 1.255 + 2.18(T

—

100) x 10

ly /IT = — 0.056 dyn cm °C -‘ y = 119.96+(T — Tm)ay/oT dyn cm ~=0.0302—1.533x104(T—T, 5)gcm~ = 6.6x iO~°C~ ES = 0.7 EL = 0.7 3 PS = 2.118 g cm3 p1=1.904gcm ~° —

J g’ °C1

si

this report. Furthermore, similar to Fig. 2a, the grid spacing is nonuniform, i.e., being finer near the melt/gas and melt/solid interfaces in order to treat the higher velocity and temperature gradients there. The nonuniform grid spacing can be generated easily through stretch functions during grid generation [14]. Similar mesh placements were used in the simulation of FZ rod growth [7,8,11,12]. In the following sections, calculated results based on the mesh II for NaNO3 and Si are discussed separately.

=

=

4.1. NaNO3

The effects of convection are demonstrated through three different sets of f~and ~7/~T. The calculated results are shown in Fig. 3. In each figure, the LHS shows the streamlines and the RHS isotherms. The arrowheads indicate the direction of flow and the asterisk indicates the zero streamline. The case of f3 0, E37/aT 0, i.e., conductive heat transfer only, is shown in Fig. 3a. =

=

Table 3 Effect of the grid mesh

Meshes

I 11 III

Total number of equations

Crystal inner radius (cm)

Crystal outer radius (cm)

Inner zone length

Outer zone length

~

I/max

Tmax

(10~g/s)

(102 g/s)

(°C)

1817 2859 4467

0.6676 0.6520 0.6471

0.8750 0.8633 0.8595

0.2166 0.2565 0.2715

0.5020 0.5140 0.517

—7.526 —7.284 —7.323

1.093 1.129 1.113

318.29 318.49 318.59

C. W Lan /Journal of Crystal Growth 141 (1994) 265—2 78

now reduced to and outer radii slightly. Fig. 3c shows ~7/~T= —0.056

~

Feed tube 0.2cm

‘1~=-6776xi0~

~i15880b0

4f==-7284xIo~

__~~

~ 1~ Lllhi11”

~

1289X1O~

II~

Cryxt/ ~

(a)

~Ibt~

Fig. 3. Calculated results for an NaNO3 tube based on: (a) conduction; (b) natural convection; (c) thermocapillary and natural convections. U~ = LI = —2x iO~cm/s; Rfi = 0.7 cm and Rf 0.9 cm. The left-hand side shows the streamlines and the right-hand side the isotherms. ‘~I~ = t”max /10 for positive i/i, ~i = ~‘mjfl /10 for negative ~I/;zlT,,, = (Tmax — Tm)/10 and L1T~= 1OLITm.

As shown, the streamlines show the material flow due to the feeding and the growth of the material. The feed front is almost flat, while the growth front is slightly convex. The maximum surface temperature is 352.8°C,which represents a rather high superheating (46°C).Due to the contractive surface tension force of free surfaces, the grown tube radii (R~1 0.595 cm and R00 0.821 cm) are smaller than the feed ones (R~~0.7 cm and R~0 0.9 cm) [13]. When the natural convection is considered, as 4°C~1and shown in Fig. 3b, where 13 6.6 X 10 a’y/aT= 0, there is one flow loop in the melt and it is clockwise in direction. The hotter and lighter melt near the outer free surface floats upward, while the cooler and heavier melt near the inner free surface sinks, thus producing the flow loop, Due to the action of the flow loop, heat transfer to the feed front is encouraged, while that to the growth front is discouraged. Consequently, the overall position of the melt zone and the location of maximum temperature shift upward as compared with those in Fig. 3a. Meanwhile, the melt/solid interfaces are still very flat. As a resuit of convection, the isotherms are distorted along the flow and Tmax (327.56°C) is significant lower than that in Fig. 3a; the superheating is =

=

=

=

=

273

20.76°C. In addition, the inner (0.607 and 0.820 cm) increase the result based on the case of

dyn cm~ °C’ and 13 6.6 X i0~ °C~,i.e., thermocapillary and natural convections. As shown, there are four flow loops in the zone. They are by even surface melt tension gradients at mainly the freeinduced surfaces, though the buoyancy force still exists. Since ~ 7/~T is negative, the surface flow is toward the melt/solid interfaces. Due to the action of the outer flow loops, heat transfer to the melt/solid interfaces is encouraged near the outer free surface, and thus results in a longer zone length =

there. On the other hand, since more heat is carried away by the outer flow loops, the inner zone length decreases significantly. As a result, the melt/solid interfaces become very convex. The inner flow loops are also smaller than the outer ones due to less space and smaller temperature gradients at the inner free surface. In addition, the shorter inner zone reduces the deformation of the inner free surface and thus restricts the further size reduction of the grown tube. Consequently, as compared with Figs. 3a and 3b, the resultant tube radii increase (R~1 0.652 cm and R~0 0.8633 cm). Again, due to the stronger convection in the melt zone, the maximum surface temperature (Tm.~ 318.49°C)is significantly lower than those in Figs. 3a and 3b. The ineffect of U~ feeding and feeding pulling speeds speedsare is shown Fig. 4; Uf. The —1 x iO~, —2 X i0~, and —3 x l0~ cm/s in Figs. 4a, 4b, and 4c, respectively. As shown, the overall position of the melt zone shifts downwards slightly with the increasing speed. Meanwhile, in spite of the three times increase in the speed, Tmax decreases only about 0.2°C. Due to the changes of zone positions, the space available for the upper flow loops decreases, while that for the lower ones increases. As a result, the upper flow loops become smaller and weaker, while the lower ones larger and stronger. The interface shapes and the dimension of the steadily grown tube at different growth speeds are further demonstrated in Fig. 5 for comparison. As shown, in addition to the changes in zone =

=

=

=

C. W~Lan /Journal of Crystal Growth 141 (1994) 265 —278

274

T~~HIHIT1~’T4~~’i~ij[‘~i~m1—-—] 4~~~

U

i -loin

I

U1=U0=lxlO

1=U =2xld~cm/x~

~ ~

Feed~ tube~

u1

U

‘I

3l0cm/x-~

—

~‘2~xi

I

~°

141

~

—-7.284xI0

~ I

W~~9.825lxIC°

I

_

I~IIl

I

~

I

‘.‘.,II’I’II.III

4

-

¶~~-L s’,

i~ III

x~.

~

I\2s9x1I32~~W,,,o.I.2625Xla2~

I(

~

I

W,.~,2’xIll

— xl

~

~

I~

I

._~)E, — —

0) ~

HEW

~I

I I i-t’ici till),

~ 3

1411=.8643x10)

‘I4~iii I

~ 2..Io°EEI/

)b) 4~=I08T,

Fig. 4. Effect of growth speed: (a) Uf = (J~= —1 (b) U~ = U~= —2x io~ cm/s; (c) U~ = U 0 = —3 NaNO3 tube of R~,= 0.7 cm and R10 = 0.9 cm.

~I ,-

I ~lI~I,l.l/

II

~

x i0~ cm/s; x i0~ cm/s.

(Ii)

Fig. 6. Effect of pulling speed: (a) b~ —1.75 x i0~ cm/s; (b) ~ = —2 x iO~cm/s; (c) ~ = —2.23 x i0~ cm/s. NaNO tube of R14 = 0.7cm and Re,,

positions, the shapes of melt/solid interfaces and the dimension of the tube crystal are not much affected by the growth speed. Nevertheless, as the growth speed increases, the growth front becomes andinthe become smaller.less Theconvex decrease the crystal crystal radii size, presumably, is due to the longer zone length, which may be caused by the stronger convection of lower flow loops, In FZ crystal growth, feeding and pulling

I ‘I

=

0.9 cm. Uf = —2x i0~ cm/s.

Fig. 6 shows the effect of the crystal pulling speed, while the feeding speed is fixed. The typical4case discussed earlier, where U,~ comparison. Uf —2 x cm/s, is shown in Fig. 6b for i0 When the pulling speed is reduced about 12% (U~ — 1.75 x iO~ cm/s), as shown in Fig. 6a, the tube can grow bigger (thicker wall) than the feed due to the mass continuity. As compared with Fig. 6b, the grown tube radii (R 01 0.681 cm and R00 0.911 cm) are much closer to the feed =

=

=

=

speeds usually can be controlled independently,

Feed tube

=

ones. As a result of the larger grown tube, the melt zone is shorter, while melt/solid interfaces are more convex. When the pulling speed is further reduced, the feed and growth fronts at the inner free surface tend to freeze together. On the contrary, as the pulling speed is increased to 2.23 X i0~ cm/s, the crystal radii decrease (R~~ 0.578 cm and R~0 0.788 cm) and the zone length increases significantly, even though the increase in the speed is only about 10% from Fig. 6b to Fig. 6c. As shown in Fig. 6, the inner convection loops =

II

Melt

/ 1

1 Fig. 5. Effect 4 of (—growth — —), speed —2X10~ on the( interface ) and shapes; —3x10~4 Uf = U~ cm/s =—1x10 (— - —) correspond to Figs. 4a, 4b and 4c, respectively,

=

are smaller than the outer ones due to smaller space for fluid flow, and the loops become larger with the increasing pulling speed. In addition, as the pulling speed increases, the zone length increases, leading to a highly necked zone shape. Accordingly, the upper flow cells become smaller and the lower cells larger.

C. W. Lan /Journal of Crystal Growth 141 (1994) 265—2 78

the melt and to counterbalance the contractive surface tension force (due to the curvature along the azimuthal direction) is required, which in turn is generated by a convex inner free surface. An additional force by the overpressure is simply to counterbalance the contractive force and to support the melt. Consequently, the inner free surface no longer needs to be convex. Due to the larger but thinner grown tube

0

14cm

Feed i.4cm

1.4cm

tube =

cm

0.2cm

hj 0

cm

=

cm

=

0075

02

3

W=-7.284x10°

~ 1~ V

=LI289xi05~

W,=~9.097xi0

W=-7.797x10~

‘~

W,,,~=i.l823x102

~

It

(bI ST~

1

H W,,,.~,el.3529xi02

i

‘°

275

Fig. 7. Effect of internally applied overpressure head: (a) = 0 cm; (b) h = 0.075 cm; (c) h. = 0.2 cm. NaNO

4 3cm/s. tube of

caused by the overpressure, the melt zone length and the maximum surface temperature increase. The longer zone also provides a larger space for convection. Therefore, the magnitudes of 1(1mm and il, also increase. 4.2. Si max

R~1O.7cm and R~0=0.9cm. U~=U0=—2x10

From the cases discussed so far, the grown tube is always smaller than the feed. This is mainly caused by the contractive force of surface tension on the free surfaces [13]. For smaller tubes, the principle curvature in the azimuthal direction is larger, and thus generates a larger contractive force, This not only causes the tube to grow smaller, but also reduces zone stability due to the highly distorted free surfaces [13]. In other words, the FZ growth of smaller tubes may be more difficult than that of bigger ones. In order to overcome the difficulties, and at the same time, to control crystal dimension and to enhance the zone stability, an overpressure can be introduced into the tube to counterbalance the contractive surface tension force [13]. This idea is illustrated in Fig. 7. Fig. 7a shows the case without the overpressure, and therefore the steadily grown tube is smaller than the feed. However, after pressure a small over headas h. P~/p~g), 0.075 cm (the headpressure is defined iS applied (Fig. 7b), a tube with radii close to the feed ones can be obtained (R~ 1=0.682cm and R~ 0 0.886 cm). By increasing the overpressure further (Fig. 7c), a larger tube can be obtain (R~ 1 0.749 cm and R~0 0.939 cm). As shown, the inner free surface is no longer convex toward the centerline. On the contrary, it bulges outwards. Indeed, the free surface shapes are determined by the normal stress balance. Without the overpressure, a surface tension force to support =

I

=

=

,

The calculated results for silicon (Si) are based on the physical properties shown in Table 4. According to the table, the Prandtl number for Si, i.e., 0.01, is much smaller than that of NaNO 3 (9.12). The growth angle of Si, i.e., 11°, is also significant larger than that of NaNO3 (2.5°). In the following discussion, the conditions of T~ 2097°C, Tam 64°C, a 1 cm, Uf U~ 1 X i0~ cm/s, Rf1 0.7 cm, and R~0 0.9 cm are used unless otherwise stated. The effect of convection is also evaluated through three different sets of 13 and ~7/~T. The calculated results are shown in Fig. 8. As shown, =

=

=

=

=

=

Table 4 Physical properties of Si [27,28]

Tm = 1410°C ~~H=1803Jg~

—3

—2o

—i

— . cm i... k =064Wcm~ °C1 k50.22(T +273)/(T+273)Wcm C~= 1.038 J g~°C~ 1°C1 Cp~=1.059Jg~ 87/8T = —0.43 dyn cm~°C~ t y=720+(T—Trn)ay/aT dyncm

=

~

=

= 0.3 3 Ps == 2.33 PL 2.55 gg cm cm3 ~° = 11

~L

t

-

‘°C

=

—

C. W. Lan /Journal of Crystal Growth 141 (1994) 265—278

276

shows the case with a lower pulling speed (U —0.8 x 10 ~ cm/s). As shown, the calculated

-~

1 4cm ~

Feed

~

=

4cm

1 4cm

tube~

result is similar to that of NaNO (Fig. 6a). The size of the grown tube increases significantly as

-~ H

4

Vooc5498X102

W,

W,~ 1c)728xI0

the pulling speed decreases. Since the inner and

10-0)28)

outer zone lengths are short, the deformation of

K 1

~. 1)1

0299)

H

-

Crystal Ib)

(a)

~

Ic)

~

Fig. 8. Calculated results for an Si tube based on: (a) conduction; (b) natural convection; (c) thermocapillary and natural convections. U1 = U4 = 1 X 10 cm/s; R1 = 0.7 cm and R10 =0.9 cm. iT —(T —T )/5and iT=lOziT . m

max

m

s

m

despite of the much stronger convection, the effeet of convection is similar to that for NaNO3 shown previously (Fig. 3), but less pronounced due to the much smaller Prandtl number of Si. For example, the difference of Tmax between Fig. 8a (conduction only) and Fig. 8c (natural and thermocapillary convection) is only 1.8°C. However, as compared with Fig. 3, the inner free surface is convex toward the axis more, while the necking of the outer free surface is reduced. Not surprisingly, because the Bond number, Bo pLgRfl/y, for Si is much smaller (1.7) than that for NaNO3 (8.48), the effect of contractive surface tension for Si is more significant. More interestingly, for Si, unlike NaNO3 (Fig. 3), when both natural and thermocapillary convections are considered (Fig. 8c), the grown tube radii do not increase but decrease. In fact, for Si although the thermocapillary convection enhances heat transfer toward the melt/solid interfaces at the outer surface, heat can still be brought into the inner free surface effectively by conduction. As it is, the outer zone length increases only slightly. On the contrary, since the inner free surface is convex toward the axis more space is available for flow to develop, which in turn enhances the convection, leading to a slightly longer zone length there. As a result, the grown tube radii decrease. The calculated results based on three different pulling speeds are illustrated in Fig. 9. Fig. 9a =

speed is ~:reased to ~ feeding speed, as demonstrated in Fig. 9c, the size of the grown tube decreases significantly while the outer free surface necks and the inner surface bulges inward severely. Clearly, to obtain a stable growth may be difficult if the pulling speed is further increased. Again, for all the cases the growth fronts are very convex. The effect of the internally applied overpres.

sure is demonstrated in Fig. 10. Without an over pressure (Fig. lOa) the inner free surface is convex inward significantly. After a small overpressure head (h1 0.5 cm) is applied (Fig. lob), the surface deformation decreases and the grown tube radii increase. When the applied pressure head is increased further to 1.1 cm (Fig. lOc), the inner free surface is no longer convex inward but bulges outward, while the grown tube radii increase significantly. Because of the increase of tube size, the wall thickness of the grown tube decreases as well. Owing to the slender melt zone, the radial temperature gradients decrease, while the axial =

~

[~Hp~ U11 OslO cm)s

Uf=-l OslO

U

1 OxlO

* -~

3285

~

-~

H

0.3283

~,

Voc~0347)

p

v—~=’

~V’

~ ~

H~

0)319

V,=o.299~ ~

W~o2867

Crystal AT,=OST. I~.

O.xxIs~s/s

(a)~

1II4

ioxioymis

(b)

Fig. 9. Effect of pulling speed: (a) (J~= —0.8x10 U~=

R1

=

(C)

cm/s; (b) cm/s; (c) U,~= — 1.1 X 10 ~ cm/s. Si tube of 0.7 cm and R10 = 0.9 cm. U~= — 1 x 10 cm/s.

—

1

x 10

C. W Lan /Journal of Crystal Growth 141 (1994) 265—2 78 -—--S

~

~

1.4cm h1=Ocm

1.4cm h1=O.5cm

(2) In both NaNO3 and Si, at least for the conditions examined in the present study, thermocapillary convection dominates in the melt zone. As compared with the case of no convec-

-

1.4cm h~=i1cm

.

0.2c

Feed tube W~,c~O 3283

W=-O.3446

P o~

0

I W~,o.c02993

H W,=O.3283

~

V,,~=-O5004 4889

-~

-~

(a)

Crystal)b)AT,100T

277

cl

Fig. 10. Effect of internally applied overpressure head: (a) = 0 cm; (b) h = 0.5 cm; (c) h = 1.1 cm. 3 Sicm/s. tube of ~ = 0.7 cm and R~0= 0.9 cm. U1 = U~= —lx i0

ones at the inner free surface increase, leading to the stronger thermocapillary convection there. Again, as shown, the control of the grown tube radii by an overpressure is very effective, even though the idea has not yet been used in practice. Finally, since the FZ method has not yet been used widely in the growth or the refining of tube crystals [1,6], quantitative comparison of the calculated results here with experiments is not available. However, the reported grown tube radii of tungsten crystals by Glebovsky et al. [6] could be fitted well by the thermocapillary model reported previously [13] and the present model, if the heat input parameters are adjusted properly. The simple comparison for the thermal-capillary model is discussed elsewhere [131. Again, the effect of convection for the tungsten system is also very small.

tion, the melt zone is lengthened at the outer free surface but shorten at the inner free surface, resulting in more convex melt/solid interfaces and a change in the steadily grown tube size. However, the effect of convection on the resultant tube radii is complicated and depends on the Prandtl number of grown materials. For example, with the natural and thermocapillary convections, the radii of the grown tube increase for NaNO3, but decrease for Si. force, (3) the Because inner of free thesurface contractive is convex surface toward tension the axis, and the grown tube is smaller than the feed. However, by controlling the overpressure inside the tube, the grown tube can be adjusted to be smaller, equal, or even larger than the feed. (4) When the pulling speed is increased, both the zone length and the grown tube size change significantly; the former increases while the latter decreases. However, the wall thickness is not much affected by the pulling speed. (5) The effect of thermocapillary convection is significantly more pronounced in NaNO 3 due to its much larger Prandtl number, though the extent of thermocapillary convection is much stronger in Si due to its much more negative a7/aT.

Acknowledgements The author is grateful for the support of National Research Council of ROC under Grant No. NSC84-0402-E008-002.

5. Conclusions (1) A computer model based on the BFCFVM/Newton method is used to study steady-state heat transfer, fluid flow, and interfaces in the FZ tube growth process, considering both natural and thermocapillary convections. The inner and outer free surfaces melt/solid interfaces and the size of the steadily grown tube crystal are calculated simultaneously.

References [1] W.G. Pfann, Zone Melting, 2nd ed. (Wiley, New York, 1966). [2] W. Dietze, L.P. Hunt and D.H. Sawyer, J. Electrochem. Soc. 121 (1974) 112. [31B.S. Redkin, V.N. Kurlov and V.A. Tatarchenko, J. Crystal Growth 82 (1987) 106.

[41D.S.

Harvey, J. Crystal Growth 104 (1990) 88.

278

C.W. Lan /Journal of C,ystal Growth 141 (1994) 265—278

[5] A.A. Alioshin, N.I. Bletscan, S.F. Bogatyriov and NV. Fedorenko, J. Crystal Growth 104 (1990) 130. 16] V.G. Glebovsky, V.N. Semenov and V.V. Lomeyko, J. Crystal Growth 98 (1989) 487. [7] C.W. Lan and S. Kou, J. Crystal Growth 108 (1991) 351. 18] C.W. Lan and S. Kou, J. Crystal Growth 114 (1991) 517. [9] B. Xiong and W.R. Hu, J. Crystal Growth 125 (1992) 149. 110] B. Xiong and W.R. Hu, J. Crystal Growth 133 (1993) 155. [11] C.W. Lan and S. Kou, J. Crystal Growth 132 (1993) 578. [12] C.W. Lan and S. Kou, J. Crystal Growth 133 (1993) 309. [13] C.W. Lan, J. Crystal Growth 135 (1994) 606. [14] C.W. Lan, mt. J. Numer. Methods Fluids, in press. [15] Y. Yamaguchi, Ci. Chang and R.A. Brown, Phil. Trans. Roy. Soc. London A 312 (1984) 519. [16] A.D. Gosman, W.M. Pan, A.K. Runchal, D.B. Spalding and M. Wolfshtein, Heat and Mass Transfer in Recirculating Flows (Academic Press, London 1969). [17] C.W. Lan and S. Kou, mt. J. Heat Mass Transfer 35 (1992) 433. [18] T. Surek and B. Chalmers, J. Crystal Growth 29 (1975) 1. [19] J.F. Thompson, F.C. Thames and C.W. Mastin, Boundary-fitted curvilinear coordinate systems for solution of

partial differential equations of field containing any number of arbitrary 2-dimensional bodies, NASA Report: NASA-CR-2729, July 1977. 120] T.F. Coleman and J.J. More, SIAM J. Numer. Anal. 108 (1983) 340. [21] H.L. Saito and L.E. Scriven, J. Comput. Phys. 42 (1981) 53. [22] M. Seager, A SLAP for the Masses, Lawrence Livermore National Laboratory Technical Report, UCRL-100267, December 1988. [23] H.B. Keller, Numerical Solution of Bifurcation and Nonlinear Eigenvalue Problems, in: Applications of Bifurcalion Theory, Ed. P.H. Rabinowitz (Academic Press, New York, 1977). [24] L.R. White and H.T. Davis, J. Chem. Phys. 47 (167) 5433. [25] G.J. Hanz, J. Phys. Chem. Reference Data 1 (1972) 583. [26] F. Preisser, D. Schwabe and A. Scharmann, J. Fluid Mech. 126 (1983) 545. [27] J.L. Duranceau and R.A. Brown, J. Crystal Growth 75 (1986) 367. [28] T.A. Kinney, D.E. Bornside, R.A. Brown and K.M. Kim, J. Crystal Growth 126 (1993) 413.