- Email: [email protected]
Three-dimensional convection and solute segregation in vertical Bridgman crystal growth
,. . . . . . . .
C R Y S T A L Q R O W T H
ELSEVIER
Journal of Crystal Growth 167 (1996) 320-332
Three-dimensional convection and solute segregation in vertical Bridgman crystal growth M.C. Liang, C.W. Lan * Chemical Engineering Department, National Central Universi~', Chung-Li 32054, Taiwan, ROC
Received 25 October 1995; accepted 30 November 1995
Abstract
Perfect axisymmetry is difficult to achieve in the vertical Bridgman crystal growth. Nonaxisymmetry is, to a certain degree, always present in a real crystal growth configuration. To study how fluid flow and solute fields in the process are affected by slight deviations from axisymmetry, a three-dimensional finite-volume/Newton method is developed. Pseudosteady-state heat and mass transfer, fluid flow, and the growth front are solved simultaneously. The calculated results for axisymmetric cases are compared with those by two-dimensional finite volume and finite element methods, and they are in excellent agreement. Calculated results also show that three-dimensional flow and solute fields are easily induced by imperfect axisymmetry. A slightly tilted ampoule or a small temperature gradient in the azimuthal direction could break axisymmetry and dramatically affect flow patterns and solute segregation.
1. Introduction
The vertical Bridgman (VB) technique has been widely used in growing high quality semiconductor crystals (e.g. Refs. [1-4]). With the fast development of optoelectronic devices and very large scale integrated circuits (VLSI), the need for crystal substrates with lower defects and higher dopant (solute) uniformity has increased dramatically. In order to meet the demand of higher quality substrate crystals, further process fine tuning and improvement for the VB technique are inevitable. During the VB crystal growth, heat and mass transfer dictates the growth interface and dopant segregation and thus plays a crucial role in crystal quality. To grow crystals with
* Corresponding author. Fax: +886 34227382; E-mail: lan @che730.che.ncu.edu.tw.
a better quality, precise control of heat and solute fields is necessary. Computer simulation provides an effective way to better understand the process and ultimately to develop a better control strategy. Numerous two-dimensional (2D) numerical studies on the VB crystal growth have been reported in the past decade (e.g. Refs. [5-12]). They were all based on the assumption of axisymmetry. The fluid flow, heat and mass transfer, and growth interface were simulated simultaneously by either the finite element method (FEM) [5-10] or the finite volume method (FVM) [11,12]. A detailed review of these 2D simulations can be found elsewhere [13]. In practice, however, perfect axisymmetry is difficult to achieve, and 3D flows and thus solute segregation can result. The axisymmetry can be broken by a slightly tilted ampoule, an imperfect furnace, or an anisotropic crystal or ampoule. For example, Neuge-
0022-0248/96/$15.00 Copyright © 1996 Elsevier Science B.V. All rights reserved PII S0022-0248(96)00209-6
M.C. Liang, C.W. Lan / Journal of Crystal Growth 167 (1996) 320-332
bauer and Wilcox [14] observed 3D flows being dominant in the VB crystal growth of salol through flow visualization. Even under perfectly axisymmetric conditions, symmetry breaking is possible due to flow bifurcations. Therefore, for a detailed simulation of the VB crystal growth, a 3D model is necessary. Recently, Jung [15] used a 3D model to simulate the melt flows in the vertical Bridgman crystal with axisymmetric boundary conditions (same as his 2D model). He found that the flows were 3D but did not deviate much from the axisymmetric ones. The reason for the slightly nonaxisymmetric flows was not explained. Since the solute transport, growth interface, and nonaxisymmetric boundary conditions were not considered in his calculations, the effect of slight nonaxisymmetry on the flow structures and hence solute segregation was not clear. Pulicani et al. [16] demonstrated the effect of azimuthal temperature gradients on the 3D flows in a vertical cylinder. They found that the degree of flow asymmetry increases with increasing Prandtl number and flow intensity; the solute segregation and the growth interface were not considered. With thermal boundary conditions and an interface shape from an axisymmetric system, Kuppurao et al. [17] further illustrated the effect of ampoule tilting on 3D flows and hence solute segregation. The radial solute segregation was found significantly affected by the induced 3D flows. However, the unknown growth interface was not considered in their calculations. In fact, both flow structures and solute segregation are affected by the growth interface [5,11]. To model the system in a self-consistent manner, the growth interface should be calculated simultaneously. However, to our best knowledge, a fully 3D crystal growth simulation of heat and mass transfer coupled with the unknown growth interface for the VB process has not been reported before. In crystal growth modeling, the computation of solute fields is usually very troublesome [12,18]. In addition to resolving thin solute boundary layers, a good global solute conservation is also necessary. In fact, the global conservation of solute is a major challenge to the VB simulation even in 2D calculations [18]. For some numerical schemes, the global error could be very sensitive to the mesh used. Accordingly, choosing an appropriate numerical scheme on the 3D solute calculations is crucial in
321
obtaining a reasonable result when a very fine mesh is not possible. The FVM is well known for its excellent global conservation properties. In a recent performance comparison of the FVM and the FEM on the segregation computation of the VB crystal growth [12], the FVM proves to be more efficient, and its global solute conservation is not affected by the mesh. On the contrary, a fine mesh is usually necessary for the FEM to achieve an acceptable global error, especially when the convection is strong. Furthermore, to simulate the process in a self-consistent manner and with fast convergence, the heat and mass transfer, fluid flow, and growth interface need to be solved simultaneously by a robust and fast solution scheme. According to previous studies [5,8,9,11,12,18], the global iteration scheme, which uses Newton's method to solve field and interface variables simultaneously, seems to be the best candidate for 3D computation. In the present report, a 3D FVM/Newton method is developed to study how slight deviations from axisymmetry affect flow patterns and radial solute segregation in VB crystal growth. Numerical solutions for axisymmetric cases are first qualified by 2D FVM and FEM codes. The nonaxisymmetry due to a slightly tilted ampoule and an imperfect furnace is then considered to illustrate the induced 3D flow and solute fields.
2. Model description
A generic VB crystal growth system is illustrated in Fig. 1. Since axisymmetry is no longer assumed here, the system is described by the Cartesian coordinate ( x , y , z ) . The nonaxisymmetry could be introduced by the tilted ampoule, which is represented by the tilted angle y, or the nonuniform ambient temperature profile T~(x, y , z ) . In Fig. 1, each region is characterized by a set of physical properties. When a pseudo-steady-state is achieved at the ampoule pulling speed Uamp, the growth rate can be set to be Uamp as well. The flow, temperature, and solute fields, as well as the melt/crystal interface, are also represented in the Cartesian coordinate (x, y, z). The dimensionless variables are defined by scaling length with the ampoule length L, velocity with O~m/L, pressure with Pm °~m/L 2 2, and solute concen-
M.C. Liang, C.W. Lan / Journal of Co'stal Growth 167 (1996) 320-332
322
the Prandtl number, Sc the Schmidt number, and 0m the melting temperature, which also serves as a reference temperature. Two important dimensionless variables, Ra s and Ra T, in the source term of the momentum equation are the solutal and thermal Rayleigh numbers, defined as follows:
i
/
g/3s CoL3 Ra s =
g/3T A TL3 ,
Ra T -
Ofm Pm
, O~m/'Pm
where g is the gravity acceleration, /3s and /3v are the solutal and thermal expansion coefficients, respectively, AT = Th - T~, and vm is the kinematic viscosity. The gravity direction eg can be decomposed into the Cartesian components (ex,ey,e ~)
,, °°"i/
eg = cos( a ) e x - c o s ( / 3 ) e y -
/
Fig. 1. Schematic sketch for the vertical Bridgman (VB) crystal growth system.
cos(y)e~,
where the o~, /3, and y are the angles between the gravity direction and the axes of x, y, and z, respectively. For a perfectly aligned ampoule, o~ = / 3 = ~-/2, and y = 0. In the present study, nonaxisymmetry is introduced only by tilting the ampoule in the y-z plane, i.e. a = ~-/2. In the crystal (c) and the ampoule (a), only heat transfer needs to be considered,
-Peie :. VO= KiV20 ( i = c , a ) , tration with its average concentration in the crystal, C 0, where a m is the thermal diffusivity and Prn the melt density. The dimensionless temperature ( 0 ) is defined as
O(x,y,z)=[T(x,y,z)-T~]/(Th-T~),
(1)
where T h and Tc are the top (hot) and bottom (cold) temperatures, as shown in Fig. 1. The pseudosteady-state governing equations describing convection and heat and solute transport in the melt (m) are as follows: V ' v = O, v" Vv = - V P + P r g 2 v + P r [ R a s ( C - R a T ( 0 - - 0m)] eg,
(2) 1) (3)
v" VO= V20,
(4)
Pr v-VC=--V2C, Sc
(5)
where v, P, and C are the dimensionless velocity, pressure, and solute concentration, respectively. Pr is
(6)
(7)
where Pe i =-PiCpiUampL/km is the Peclet number and Ki = k i / k m the dimensionless thermal conductivity of crystal or ampoule; k m is the thermal conductivity of the melt. Also, P i , Cpi, and k i a r e the density, specific heat, and thermal conductivity of the phase i (i = m, c, or a), respectively. The no-slip condition is used for the velocity on solid boundaries b' = -- Tct;ampez,
(8)
where Yc -= Pc/Pm and Uamp is the dimensionless ampoule pulling speed. Since the pseudo-steady-state is assumed, the upper open boundary is considered as an artificial boundary [5], and its velocity boundary condition is the same. The thermal and solutai boundary conditions at the melt/solid interfaces are set by the heat and mass flux balances. For example, at the growth front:
n.VOlm--n.K~VOL+PecSt(n.e:) = 0 , -n'VCIr.
(9)
Pe m Sc
--(1-K)C(n.e~)=O, Pr
(10)
M.C. Liang, C. W. Lan / Journal of Crystal Growth 167 (1996) 320-332
where n is the unit normal vector at the growth interface pointing to the melt. The Stefan number St = AH/Cp~AT scales the heat of fusion ( A H ) released during solidification to the sensible heat in the solid. The Pem = p~CpmUampL/km is the Peclet number and D the solute diffusivity in the melt. Also, K is the equilibrium segregation coefficient of the solute according to the phase diagram; K =- C J C at the growth interface, where C¢ is the solute concentration in the crystal. The temperature at the melt/solid interface is assumed to be the equilibrium liquidus temperature of the material. Solute diffusion in the solid is neglected here. Temperature at the top and bottom surface is set to the furnace temperature there. With the pseudo-steady-state assumption, an artificial boundary condition is used at the upper surface for the consistence of the overall solute balance [5]:
-n'VClm-
Pe m Sc
--(CPr
1)(n .e=).
ing. The complications of computing view factors from the sample to the furnace are neglected here.
3. Finite-volume/Newton method Owing to the unknown melt/crystal interface shape hc(x,y,z), the physical Cartesian coordinates (x, y, z) are transformed into general (nonorthogonal) curvilinear coordinates ( sc 1,~: 2, sc3), which fit all the interfaces. Fig. 2a is a sample converged mesh for computation; the mesh defines the edges of finite volumes. It should be pointed out that the interface hc(x,y,z) is unknown a priori and needs to be solved simultaneously with other field variables. During iteration, once the interface is updated, the mesh can be easily generated through an algebraic boundary-fitted coordinate transformation, and grid density control is through the stretch functions. The
(11)
The heat exchange between the ampoule and the furnace is by both radiation and convection according to the energy balance along the ampoule surface, - n . KaV0la= B i ( 0 - 0a) + Rad(0 4 - 04),
Ampoule
(12)
where n is the unit normal vector on the ampoule surface pointing outwards, B i - hL/k m the Biot number, and Rad - O'EaAT3L/km the radiation number; o- is the Stefan Boltzmann constant, while e a is the surface emissivity of the ampoule. For simplicity, the effective furnace temperature 0a is assumed to be 0a(x,y,z) = Z -- 0~( x, y,z),
Melt
(13)
3rowth nterface
where, Od(x,y,z) is a deviation temperature to describe thermal nonuniformity. A Gaussian distribution is used for 0d,
O,~(x,y,z)=AOdexp(-[(z-0.5)/a*]2},
323
Crystal
(14)
~ . ~
Y
x
and
= a0 p cos(
),
(15)
where A 0dp is the dimensionless peak (maximum) deviation temperature difference along the azimuthal angle q~ and a* a dimensionless parameter for the width of the distribution. In the present study, the nonaxisymmetry for a perfectly aligned ampoule is introduced by a nonzero A 0ap for nonuniform heat-
Ca)
(b)
Fig. 2. (a) 3D finite volume mesh for computation; a part of the mesh is removed to illustrate the finite volumes inside; (b) a typical finite volume.
M.C. Liang, C.W. Lan /Journal of Cr3,stal Growth 167 (1996) 320-332
324
transformation is similar to that used in the 3D thermal-capillary simulation of floating-zone crystal growth, and that has been described elsewhere [19]. The governing equations and boundary conditions can then be approximated by a finite volume method. As shown in Fig. 2a, the physical domain is divided into a number of finite volumes bounded by cell faces. A sample of the finite volume is illustrated in Fig. 2b. In general, each volume is an irregular hexahedron defined by eight vertices. The coordinate invariant strong conservative-law form of the transport equations (Eqs. (2)-(5) and (7)) can be written
where W is a nonzero weighting function. This isotherm approach without the weighting function has also been used by Ettouney and Brown [21] and Lan and Kou [22]. The purpose of using the weighting function here is to avoid the numerical breakdown during LU decomposition in the solution stage, since no pivoting is used there. A linear combination of the interface variables near the point of interest is used for the weighting function. Eq. (17) and the above weighted isotherm condition yield a system of nonlinear algebraic equations, which can be written symbolically as
as
f ( x ) =-f(v,P,r,C,hc)
V.I+S=0,
(16)
where I describes the flux vectors of the mass, momentum, energy, and solute. The source term S consists of the body forces for the momentum equation, but it is zero for other equations. If Eq. (16) is integrated over a general control volume AV in physical space and the Gauss divergence theorem is then applied, a discrete flux balance equation can be obtained, i.e.,
..~ I . A I [e + l - a l l w + l ' a 2 [ n + l ' a 2 1 ~ + / ' a 3 1 t
+I.A3Ib + SAV=O,
(17)
where A i, i = 1 to 3, is the area vector defined on the cell surfaces (see Fig. 2b). Working on the control volumes, one could derive the flux components using the central difference scheme. However, the mass fluxes used in the equation of continuity need to be taken care of separately to avoid the checkerboard pressure oscillation due to the velocity/pressure decoupling. A modified momentum interpolation scheme [20] is used in the present calculations. The detailed description of the numerical issues concerning the primitive variable formulation on the collocated grid for incompressible flow is described elsewhere [20]. The boundary conditions can be imposed easily during the finite volume approximation of the cells next to the boundaries. To the growth interface hc(x,y,z), a weighted melting-point isotherm approach is adopted
W[ O(x,y,hc) -
0re(C)] = 0,
(lS)
(19)
This nonlinear equation set is solved simultaneously by Newton's method for all variables. Starting from an initial approximation to vector of unknowns x °, successive updates are constructed as x , + l = x " + 6 n+l,
(20)
and the correction vector 6 "+1 is the solution of the linear equation set J6 "+1:
fAl'dA + f±vSdV
= 0.
-f(x").
(21)
The components of the Jacobian matrix J, formed by explicit differentiation as f j - Ofi/Oxj, represent the sensitivity of the residual vector to the solution vector. They are all determined by the forward difference [23,24]. Since Eq. (21) is large and sparse, the use of an efficient sparse matrix solver is necessary. The ILU(0) [25] preconditioned GMRES method [26] is used, and it is found effective and efficient in the computation here.
4. Results and discussion
For comparison purposes, we consider the gallium-doped germanium (GaGe) growth in the Grenoble furnace investigated by Adornato and Brown [5] in this study. The same system has also been studied by Lan and Chen [12] as a benchmark problem for comparing the performance of a 2D FVM (stream function-vorticity ( q - w ) formulation) and a 2D Galerkin FEM (primitive variable formulation). The physical properties of GaGe and some input parameters used in the study are listed in Table 1. The condition of Ra s = Rad = 0 used in Ref. [5] is also
M.C. Liang, C. W. Lan / Journal of Crystal Growth 167 (1996)320-332
[5,12] are tested again. Since the computation is performed in the ( x , y , z ) coordinates rather than the cylindrical coordinate (r,z), and no axisymmetry is assumed, fully 3D computation is still necessary. Fig. 3a-3c illustrate the 3D calculated results showing the effect of convection for the axisymmetric growth ( y = 0 and A0dp = 0). In each figure, the velocity field on the y - z plane of the system is illustrated on the left, and the solute field on the right. The solute distribution (on the x - y plane) in the grown crystal is illustrated underneath the solute field of the y - z plane. As shown, all solutions are perfectly axisymmetric with respect to the z axis as expected. Without the buoyancy force (Ra T = 0), as shown in Fig. 3a, the flow in the melt is due to ampoule translation, and the solute transport is mainly through diffusion. The radial segregation shown on the x - y
adopted here for comparison. In this system, the solute concentration is so low that the solutal effect is ignored. The solutal effect of a nondilute alloy can be considered if necessary [12]. The mesh used in the following calculations is similar to Fig. 2a. It contains 29 176 finite volumes leading to 111 572 nonlinear equations. All computations are performed in the IBM/ES9000 supercomputer, and one Newton iteration takes about 48 min CPU time (only one CPU is used for computation). The convergence of Newton's iterations is quadratic in all cases, and it takes only four to five iterations to get a convergent result with the L~-norm of 1 × 1 0 - 9 . 4.1. Axisymmetric solution and comparison with 2D cases To validate the code by different numerical approaches, the axisymmetric cases studied before RaT=
RaT = 0 [~-.-
325
RaT= 10 6
10 5 ....
! i
I
....
#Jt,ILJ.ttata.t.,lltili.~U JamUJ~:H,ILt~l.l.ti~.m~ laLm:;~,k~m~HhJl~:l,.illl I¸
. . . . . . . .
I
!
. . . . .
i
!
!:::: :,5: i;ii ii; :¸ ::ii
x Cm.~ = 1 2 . 3 9 *
C.~
+ Cc~ Cc...
* Cmt..~ = 1 . 0 2 7
1.122
+ CCm,~ = 1 . 2 2 9
- Cc,~., = 0 . 8 6 2
Cemt. = 0 . 8 5 4
*
= 1.078
+ Cc~=
0.880
Ca)
~
= 1.012
= 1.011
=
1 -- I 0 0 c m / s x Cm~ = 1 4 . 1 3
i ---~ = 5 0 0 c m / s x C,~¢ = 1 2 . 8 9
1
C.am
(b)
Co)
Fig. 3. Effect of convection on the flow and solute fields for axisymmetric growth; (a) Ra T = 0; (b) Ra w = 1 0 - 5 ; (c) Ra T = 10 6 in each figure, the left-hand side is velocity field and the right-hand side the solute fields. The upper part of the solute fields is taken from the ( y - z ) plane, while the lower one from the (x-y) plane.
M.C. Liang, C. W. Lan / Journal of Crystal Growth 167 (1996) 320-332
326
Table 1 P h y s i c a l p r o p e r t i e s a n d i n p u t p a r a m e t e r s [5,12]
GaGe p~ = 5.5 g c m - 3 Pm = 5.5 g c m - 3 Tm = 9 3 7 . 4 ° C ~H = 460 J g- i k~=0.17Wcm
I oc-I
k m=0.39Wcm
I °C
t
Cpc=Cpm=O.39Jg
i oC
/z= 0.00715 g cm-I s ]3T=5X10 4o C I
I
,8 s = 0 ( t o o l % G a ) D=2.1XI0 4cm
I
t 2/s
K = 0.087
Graphite (ampoule) p ~ = 1.8 g c m
3
k~ = 3.26 W c m - i ° C - i Cpa=l.814Jg
- l oC - I
ea = 0
Other input parameters L=7cm R c = 0.5 c m
R,
= 0.7 c m
T h = 1112.4°C T~ = 7 6 2 . 4 ° C Uam p = 4 X 10 - 4 c m / s h=46.571Wcm 2oc
I
here. Fig. 4 illustrates the comparison of the radial (v r) and axial (G) velocities on y - z plane with the ones from the FEM for Ra T = 1 0 6. As shown, they are in excellent agreement. The comparison of the calculated growth interfaces is shown in Fig. 5a for Ra T = 1 0 6. To facilitate the comparison, the scale in the vertical axis is enlarged, so that the small differences can be seen. As shown, they are also in good agreement. More importantly, an excellent agreement is also obtained when the calculated radial solute segregations by the three numerical approaches for different Ra T numbers are compared in Fig. 5b. The averaged solute concentration should be equal to one, which is often used for checking the global solute balance. For both 2D and 3D FVMs, the global balance is satisfied exactly. However, very small global errors still exist in the FEM. For an even larger Ra T number, it requires a finer mesh for the 3D computation to get an accurate result. Unfortunately, it is too costly to do so because of the limitation of present computing resources. From previous studies [5,12], the maximum radial solute segregation occurs at Ra T = 1 0 6 for the GaGe system. A further increase in Ra T leads to a better mixing and thus a smaller solute segregation. Accordingly,
a=lcm y = 0°_1.5 ° ATdp = 0 - 1 ° C
plane is caused by interface deflection. At Ra T = 105, small convection loops near the growth front are induced. Basically, they are induced by radial temperature gradients due to interface deflection. The solute field is also slightly distorted by the convection. As Ra T is further increased t o 1 0 6 , the convection becomes stronger, and the maximum velocity is now much larger than the growth rate (4/xm/s). As a result, the solute field is highly distorted, and the radial segregation becomes larger. For the convenience of further discussion, the degree of the radial solute segregation is defined as the maximum solute concentration difference in the crystal, i.e. (Ccm,x C ..... ) X 100%. Since the solutions from the 2D Galerkin FEM and the 2D FVM [12] are also available for the axisymmetric cases, a detailed comparison of our 3D computation with these two approaches is also made
3D FVM x maX(Vz}= 0 . 0 0 8 5 1 ( c m / s )
Vz
Vr
• min(vz)=-0.00372 + maX(Vr)= 0 . 0 0 1 2 3 - min(vr)=-0.00347
.......
2D FEM
x maX(Vz)= 0 . 0 0 8 7 4 ( c m / s )
f
* min(vz)=-0.00380 + max(vr}= 0 . 0 0 1 3 8 - min(vr}=-0.00368
f
F i g . 4. C o m p a r i s o n o f 3 D c a l c u l a t e d radial ( v r) a n d axial ( t , : ) v e l o c i t y f i e l d s w i t h the o n e s b y 2 D F E M .
M.C. Liang, C.W. Lan / Journal of Crystal Growth 167 (1996) 320-332
the effect of 3D flows on the solute segregation is also expected to be more significant in this region. The axisymmetric results for higher Ra T by the 2D FVM and FEM can be seen elsewhere [12]. 4.2. Tilted ampoule In practice, a perfect ampoule alignment with the gravity may not be possible. During crystal growth, a 3.5
¢
i
i
i
i
i
b
i
i
a
3.48
. . . . . . . . 2D-FVM(~¢-o~} -........ 2D-FEM
/ I
/
3.46
~
3.44
3.42 3.4 3.38 3.36
[
I
[
i
0.1
[
[
0.2
i
i
0.3
I
0.4
0.5
r (cm)
1.4
i
!
i
~
i
t
i
i
r
b
1.3 / R a T -- 1 0 6 1.2
~
~
/
RaT =105 Cc 1.1 RaT = 0
1
...... 09
~
_
2D_VVM¢ ,-,o) I
.......
0.8
i
0
i
0. I
i
i
i
i
0.2 0.3 r (cm)
l
l
0.4
l
0.5
Fig. 5. Comparison of 3D calculated results with those by the 2D FVM and 2D FEM; (a) interface shape (for Ra T = 106); (b) radial solute segregation.
327
small tilt (less than 1°) from the gravity could be very common, and this may break axisymmetry. Fig. 6 shows the effect of tilt angle for Ra T = 105. As shown in Fig. 6a, even though the tilt is only 0.5 °, the axisymmetry of lower cells is broken, and the convection becomes asymmetric. The lower right cell grows larger, and the lower left cell shrinks slightly, while a new big flow cell develops in the main body of the melt. The fluid flow also becomes faster near the right ampoule wall, and slower on the other side. As a result, the solute field is distorted slightly by the flow and becomes asymmetric. The maximum solute concentration (Ccmox) shifts slightly from the center to the left, and the radial solute segregation increases to 31.4% as compared with 26% in Fig. 3b. As the tilt angle is increased to 1.5 ° in Fig. 6b, the new flow cell becomes stronger. The lower left cell becomes much smaller, and the right cell grows much bigger and joins the big new flow cell. The fluid flow near the right wall is also enhanced. As shown, the flow is now fully 3D, and is not possible to be approximated by an axisymmetric model. It should be pointed out that the flow and solute fields in the x - z plane are still symmetrical. Due to the 3D fluid flow, solute transport becomes 3D as well. The solute field is distorted even more along the main flow direction, and the segregation is increased to 39.8%. Cc~x is also pushed to the left even more. As the convection becomes stronger (Ra T = 106), as shown in Fig. 7, the effect of tilt on the flow structure and hence the solute segregation becomes even more significant. Even a tilt of only 0.5 ° has a dramatic effect on the flow and solute fields (Fig.7a). As shown, the solute field is now highly distorted by the convection. The affected area now extends to the upper part of the melt as well. Consequently, the radial segregation is increased to 61.5% (39.5% in the axisymmetric case). As the tilted angle is further increased, as shown in Fig. 7b, the convection becomes even stronger, and Ccm,~ shifts to the left even more. However, the radial segregation is reduced to 54.5% due to better mixing. Apparently, the 3D flow promotes the better mixing of solute in the bulk melt and leads to a dramatic decrease in the concentration difference. Cmi n is increased to 4.175 as compared with 1.027 for the axisymmetric case in Fig. 3c. It is expected that to some extent the solute segregation
M.C. Liang, C. W. Lan /Journal of Crystal Growth 167 (1996) 320-332
328
may be reduced at higher Ra T and larger ampoule tilting [ 17]. Furthermore, in the flow visualization experiments by Neugebauer and Wilcox [14], it was observed that 3D flows were dominant in the melt. In their experiments, the furnace was specially designed to achieve a good thermal axisymmetry. From the implication of Figs. 6 and 7, the 3D flows might be caused by a slightly tilted ampoule used in the experiments. The melt height in their experiments happened to be much larger than the diameter, which was similar to the cases here. Such a high aspect ratio may also promote the formation of 3D flows
7=
0.5 °
RaT=
due to stronger convection. In other words, for VB crystal growth 3D flow and solute fields could be induced even more easily with increasing melt aspect ratio. Nevertheless, we do not exclude the possibility of the imperfect furnace as the cause of 3D convection. The effect of nonuniform heating due to an imperfect furnace is discussed next. 4.3. Nonuniform heating
In a VB furnace, no matter how good the design is, small temperature gradients in the azimuthal direction may still exist. If the ampoule is not rotated,
10511
7 = 1.5 ° !
:,T~11!,,,
I..~., ....
"~J;'il
I
lit!
::ji]~l~'iits,~ I..,,,:: .... 't~'~*lllii;it11*~l ! ir'~ll'll,
!
,r'H'il, ,,rrll[~
J I
,
i iij1111
,
' 'r11'~1
c
............
"
I
i
L~
i
., ILt~
i~
I
1 -~ = 500 cm/s
1
= 300 x Cm~ = * Cr~n = + Ccm.x = Ccmm =
x Cm,x = 1 2 . 9 7 * Cmtn ~- 1 . 0 1 4 ae
Ccmax
=
1.128
Ccml~ = 0 . 8 1 4
cm/s
13.42 1.048 I. 168 0.738
(a) Fig. 6. Effect of tilt angle for Ra T = 105: (a) 3' = 0.5°; (b) 3' = 1.5 °,
(b)
~
M.C. Liang, C.W. Lan/Journal of C~. stal Growth 167 (1996) 320-332
which is quite common in practice, the temperature nonuniformity may be even worse. To consider this effect, the deviation temperature distribution, Od(X,y,z), is used to introduce a small temperature nonuniformity in the azimuthal direction. Fig. 8 shows the effect of nonuniform heating for Ra T = 10 5. In Fig. 8a, the maximum temperature difference is I°C (or ATap = 0.5°C); the overall furnace temperature difference, T h - Tc, is 350 ° here. As shown, the flow and solute fields have become asymmetric near the growth interface. Since the thermal nonuniformity is only imposed near the growth front, the flow far away from the interface is not affected much,
7 = 0.50
RAT=
which is different from the cases with a tilted ampoule. In each figure, the ambient temperature Ta is higher on the left but lower on the right introducing thermal gradients to drive the flow clockwise. As a result, the right cell becomes stronger, while the left cell weaker. As shown in Fig. 8a, the solute segregation is increased to 35% as compared with 26% in the axisymmetric case (Fig. 3b) and 31.4% in the tilted case (Fig. 6a). Apparently, at Ra T = 105 the I°C nonuniformity seems to have a larger effect than a half-degree tilt. With a larger temperature difference (ATdp = I°C), as shown in Fig. 8b, the effect becomes even more significant. The solute segrega-
10 6 ]
7 = 1.5 °
11ilt llJ:~
,,qU~]TII,, ""'~tilillliLL'
il.E~II'l~
\,
,rTTIT~TTTI~ ' " ! J I I I I H I ] ,1~ ,~THITtI~,, ' " J l l ! l J l l l l J ' ,,,TIT;1~I~,.. i " ' ~ l ; l l l l l l ' l ' 1
i
\
1
= -cm/s 100 x Cm= = 15.01
--~ =--
1
cm/s
50 x Cm~ = 1 5 . 0 3
* Cm~ = 1 . 9 6 8
b-
* Cm~ = 4 . 1 7 5
+ Ccm= = 1 . 3 0 6
+ Ccm= = 1 . 3 0 7
Cc~n = 0 . 6 9 1
Ccm~ -- 0 . 7 6 2
(a)
329
{b}
Fig. 7. Effect of tilt angle for Ra T = 106: (a) Y = 0.5°; (b) 7 = 1-5°.
M.C. Liang, C.W. Lan / Journal of Crystal Growth 167 (1996) 320-332
330
tion is increased to 46%; however, the flow and solute fields in the upper half of the melt are still not affected much. Again, it should be pointed out that 1 or 2°C difference is quite trivial as compared with the overall ambient temperature difference, which is 350°C here. However, quite unexpected, its effect on flow structures and solute fields is significant. For stronger convection (Ra T = 106), as shown in Fig. 9, the effect of nonuniform heating becomes even more significant leading a larger radial solute segregation. In Fig. 9a (ATdp = 0.5°C), the fight flow cell becomes much stronger than the left one, and the solute field is highly distorted along the flow direc-
ATdp =
tion near the interface. As a result, the solute segregation is increased to 63%, as compared with 37.5% in the axisymmetric case (Fig. 3c) and 61.5% in the tilted case ( y = 0.5 ° in Fig. 7a). With a larger ATop of I°C, the flow becomes fully 3D at the bottom, and the left flow cell almost disappears. The flow becomes stronger as well. Due to the much stronger convection, the solute field near the interface is distorted even more; CC~a~ is also pushed to the left even more. Interestingly, the solute segregation is decreased slightly (57.8%) as compared with Fig. 9a. It should be pointed out that the maximum radial solute segregation usually occurs at an intermediate
IRa T- 10 5 ]]
0.5°C
ATdp= 1.0°C
~:llllt ] i I '. I ll:]~l:till I I ~ / I l:lltlt
--
"
~-
. . . . . . . . . . . i i i i i ! i i i i i ...........
I
I
\
1 = --- cm/s 300 x Cm~x = 1 3 . 5 5 * C~. = 1.016 + Cc~.~ = 1 . 1 7 9 - Cc~ = 0.725
I
--~ = 500 cm/s X Cmax ~- 1 3 . 0 9 *
Cmin
=
1.013
4"
Ccmax
=
1.138
Ccmm = 0 . 7 8 6
{a)
Co)
Fig. 8. Effect of nonuniform heating for Ra T = 105: (a) ATdp = 0.5°C; (b) ATdp = I°C.
M.C. Liang, C.W. Lan / Journal of Crystal Growth 167 (1996) 320-332
ATdp = 0.5°C iii I.
ii .
.
RAT=
10 6
331
= 1.0°C
ATdp
iii'. !
.
I
[ i L
i¸ - - - - - _ J ......... . . . . . . . . . . .
i............ ] .............
itlz ........ "'l
...........
1 ............. ....... L, ....
F''''XX,
,~hL,"! , ,,r ,, ,,1, t/ .......
~',~X~k~tL"
.... :r'e'elt//t/e'r~\\\\,
....',rYittt.,,...'? "X~II~L'k " , ,rll l; 'rt ,~. "'x\\'ll~llk~"]
"'\\\ ~l',ili~' I
,11[! I [ T I ~ lll]]/~.
' ii.
.,,'11 J~J]JJ ; "'~\~,'~1 I11~' ,l~,~ ~, "' i rl+ I I: ?,
1
--~" = ~ -
= 100 c m / s Cm~
1
*
em/s
x Cm~x = 15.22 • C~ = 1.172 + C e ~ , = 1.325 Ce~n = 0.74"7
15.24 * Croan = 1.092 + C e ~ = 1.326 Cc~, = 0.696
X
''C s
=
(a)
(b)
Fig. 9. Effect of n o n u n i f o r m heating for R a T = 106: (a) ATdp = 0.5°C; (b) AT@ = I°C.
flow intensity (or poor mixing) [5]. Too weak (without mixing) or too strong (well mixing) a convection results in less segregation. Therefore, at the maximum solute segregation, further increasing the flow intensity will result in a better mixing and thus less segregation. It is clear that the flow and solute fields become 3D, even though the maximum temperature deviation from uniform heating is only 1-2°C. The induced 3D flows are stronger than the 2D ones and enhance solute transport significantly, as shown by the much distorted solute fields in Figs. 9a and 9b. Such a significant effect has not been illustrated before. Nonaxisymmetric flows due to asymmetric heating
were also observed by Pulicani et al. [16]. However, because they used higher axial thermal gradients near the growth interface, their flow patterns were less sensitive to the azimuthal temperature deviation. In the present system, if Ra T is further increased, the system will be even more sensitive to the nonuniform heating. In other words, in practice a perfect axisymmetric growth could be difficult to achieve.
5. C o n c l u s i o n s
A three-dimensional finite-volume/Newton's method is developed to study the effect of slight
332
M.C. Liang, C. W. Lan / Journal of Crystal Growth 167 (1996) 320-332
nonaxisymmetry on the fluid flow and solute segregation in the vertical Bridgman crystal growth. The effects of a slightly tilted ampoule and a nonuniform thermal profile are studied. From the calculated results, it is found that slight deviations from axisymmetry can easily induce three-dimensional flow and solute fields, but their effect on the interface shape is small. As convection increases, the system becomes even more sensitive to such deviations. Compared with the axisymmetric cases, the induced three-dimensional convection is much stronger and causes a larger radial segregation.
Acknowledgements The authors are grateful for the support from the National Science Council and the National Center for High-performance Computing of the Republic of China under Grant No. NSC84-2215-E008-011.
References [1] W.A. Gault, E.M. Monberg and J.E. Clemans, J. Crystal Growth 74 (1986) 491. [2] K. Hoshikawa, H. Nakanishi, H. Kohda and M. Sasaura, J. Crystal Growth 94 (1989) 643. [3] E.M. Monberg, W.A. Gault, F. Simchock and F. Domingguez, J. Crystal Growth 83 (1987) 174. [4] E.M. Monberg, in: Handbook of Crystal Growth, Vol. 2a, Basic Techniques, Ed. D.T.J. Hurle (North-Holland, Amsterdam, 1994) p. 52. [5] P.M. Adomato and R.A. Brown, J. Crystal Growth 80 (1987)
155. [6] S. Motakef, J. Crystal Growth 97 (1989) 20l. [7] M.J. Crochet, F. Dupret, Y. Rychmans, F.T. Geyling and E.M. Monberg, J. Crystal Growth 97 (1989) 173. [8] R.A. Brown and D.J. Kim, J. Crystal Growth 109 (1991) 50. [9] S. Branton and J.J. Derby, J. Crystal Growth 121 (1992) 473. [10] D. Hofmann, T. June and G. Muller, J. Crystal Growth 128 (1993) 213. [11] C.W. Lan and C.C. Ting, J. Crystal Growth !49 (1995) 175. [12] C.W. Lan and F.C. Chen, Comput. Methods Appl. Mech. Eng., in press. [13] F. Dupret and N. Van den Bogaert, in: Handbook of Crystal Growth, Vol. 2b, Growth Mechanisms and Dynamics, Ed. D.T.J. Hurle (North-Holland, Amsterdam, 1994) p. 875. [14] G.T. Neugebauer and W.R. Wilcox, J. Crystal Growth 89 (1988) 143. [15] T. June, PhD Thesis, University of Erlangen-Nurnberg, 1993. [16] J.P. Pulicani, S. Krukowski, J.I.D. Alexander, J. Ouazzani and F. Rosenberger, Int. J. Heat Mass Transfer 35 (1992) 2119. [17] S. Kuppurao, Q. Xiao and J.J. Derby, AIChE 1995 Annual Meeting, Miami Beach, FL, USA. [18] P.M. Adornato and R.A. Brown, Int. J. Num. Methods Fluids 7 (1987) 761. [19] C.W. Lan and M.C. Liang, Int. J. Num. Methods. Eng., to be published. [20] M.C. Liang and C.W. Lan, J. Comp. Phys., accepted. [21] H.M. Ettouney and R.A. Brown, J. Comp. Phys. 49 (1983) 118. [22] C.W. Lan and S. Kou, J. Crystal Growth 108 (1991) 351. [23] A.R. Curtis, M.J.D. Powell and J.K. Reid, J. Inst. Maths. Appl. 13 (1974) 117. [24] T.F. Coleman, B.S. Garbow and J.J. More, ACM Trans. Math. Softw. 10 (1984) 329. [25] J.A. Meijerink and H.A. van der Vorst, Math. Comp. 31 (1977) 148. [26] Y. Saad and M.H. Schultz, SIAM J. Sci. Stat. Comput. 7 (1986) 856.
C R Y S T A L Q R O W T H
ELSEVIER
Journal of Crystal Growth 167 (1996) 320-332
Three-dimensional convection and solute segregation in vertical Bridgman crystal growth M.C. Liang, C.W. Lan * Chemical Engineering Department, National Central Universi~', Chung-Li 32054, Taiwan, ROC
Received 25 October 1995; accepted 30 November 1995
Abstract
Perfect axisymmetry is difficult to achieve in the vertical Bridgman crystal growth. Nonaxisymmetry is, to a certain degree, always present in a real crystal growth configuration. To study how fluid flow and solute fields in the process are affected by slight deviations from axisymmetry, a three-dimensional finite-volume/Newton method is developed. Pseudosteady-state heat and mass transfer, fluid flow, and the growth front are solved simultaneously. The calculated results for axisymmetric cases are compared with those by two-dimensional finite volume and finite element methods, and they are in excellent agreement. Calculated results also show that three-dimensional flow and solute fields are easily induced by imperfect axisymmetry. A slightly tilted ampoule or a small temperature gradient in the azimuthal direction could break axisymmetry and dramatically affect flow patterns and solute segregation.
1. Introduction
The vertical Bridgman (VB) technique has been widely used in growing high quality semiconductor crystals (e.g. Refs. [1-4]). With the fast development of optoelectronic devices and very large scale integrated circuits (VLSI), the need for crystal substrates with lower defects and higher dopant (solute) uniformity has increased dramatically. In order to meet the demand of higher quality substrate crystals, further process fine tuning and improvement for the VB technique are inevitable. During the VB crystal growth, heat and mass transfer dictates the growth interface and dopant segregation and thus plays a crucial role in crystal quality. To grow crystals with
* Corresponding author. Fax: +886 34227382; E-mail: lan @che730.che.ncu.edu.tw.
a better quality, precise control of heat and solute fields is necessary. Computer simulation provides an effective way to better understand the process and ultimately to develop a better control strategy. Numerous two-dimensional (2D) numerical studies on the VB crystal growth have been reported in the past decade (e.g. Refs. [5-12]). They were all based on the assumption of axisymmetry. The fluid flow, heat and mass transfer, and growth interface were simulated simultaneously by either the finite element method (FEM) [5-10] or the finite volume method (FVM) [11,12]. A detailed review of these 2D simulations can be found elsewhere [13]. In practice, however, perfect axisymmetry is difficult to achieve, and 3D flows and thus solute segregation can result. The axisymmetry can be broken by a slightly tilted ampoule, an imperfect furnace, or an anisotropic crystal or ampoule. For example, Neuge-
0022-0248/96/$15.00 Copyright © 1996 Elsevier Science B.V. All rights reserved PII S0022-0248(96)00209-6
M.C. Liang, C.W. Lan / Journal of Crystal Growth 167 (1996) 320-332
bauer and Wilcox [14] observed 3D flows being dominant in the VB crystal growth of salol through flow visualization. Even under perfectly axisymmetric conditions, symmetry breaking is possible due to flow bifurcations. Therefore, for a detailed simulation of the VB crystal growth, a 3D model is necessary. Recently, Jung [15] used a 3D model to simulate the melt flows in the vertical Bridgman crystal with axisymmetric boundary conditions (same as his 2D model). He found that the flows were 3D but did not deviate much from the axisymmetric ones. The reason for the slightly nonaxisymmetric flows was not explained. Since the solute transport, growth interface, and nonaxisymmetric boundary conditions were not considered in his calculations, the effect of slight nonaxisymmetry on the flow structures and hence solute segregation was not clear. Pulicani et al. [16] demonstrated the effect of azimuthal temperature gradients on the 3D flows in a vertical cylinder. They found that the degree of flow asymmetry increases with increasing Prandtl number and flow intensity; the solute segregation and the growth interface were not considered. With thermal boundary conditions and an interface shape from an axisymmetric system, Kuppurao et al. [17] further illustrated the effect of ampoule tilting on 3D flows and hence solute segregation. The radial solute segregation was found significantly affected by the induced 3D flows. However, the unknown growth interface was not considered in their calculations. In fact, both flow structures and solute segregation are affected by the growth interface [5,11]. To model the system in a self-consistent manner, the growth interface should be calculated simultaneously. However, to our best knowledge, a fully 3D crystal growth simulation of heat and mass transfer coupled with the unknown growth interface for the VB process has not been reported before. In crystal growth modeling, the computation of solute fields is usually very troublesome [12,18]. In addition to resolving thin solute boundary layers, a good global solute conservation is also necessary. In fact, the global conservation of solute is a major challenge to the VB simulation even in 2D calculations [18]. For some numerical schemes, the global error could be very sensitive to the mesh used. Accordingly, choosing an appropriate numerical scheme on the 3D solute calculations is crucial in
321
obtaining a reasonable result when a very fine mesh is not possible. The FVM is well known for its excellent global conservation properties. In a recent performance comparison of the FVM and the FEM on the segregation computation of the VB crystal growth [12], the FVM proves to be more efficient, and its global solute conservation is not affected by the mesh. On the contrary, a fine mesh is usually necessary for the FEM to achieve an acceptable global error, especially when the convection is strong. Furthermore, to simulate the process in a self-consistent manner and with fast convergence, the heat and mass transfer, fluid flow, and growth interface need to be solved simultaneously by a robust and fast solution scheme. According to previous studies [5,8,9,11,12,18], the global iteration scheme, which uses Newton's method to solve field and interface variables simultaneously, seems to be the best candidate for 3D computation. In the present report, a 3D FVM/Newton method is developed to study how slight deviations from axisymmetry affect flow patterns and radial solute segregation in VB crystal growth. Numerical solutions for axisymmetric cases are first qualified by 2D FVM and FEM codes. The nonaxisymmetry due to a slightly tilted ampoule and an imperfect furnace is then considered to illustrate the induced 3D flow and solute fields.
2. Model description
A generic VB crystal growth system is illustrated in Fig. 1. Since axisymmetry is no longer assumed here, the system is described by the Cartesian coordinate ( x , y , z ) . The nonaxisymmetry could be introduced by the tilted ampoule, which is represented by the tilted angle y, or the nonuniform ambient temperature profile T~(x, y , z ) . In Fig. 1, each region is characterized by a set of physical properties. When a pseudo-steady-state is achieved at the ampoule pulling speed Uamp, the growth rate can be set to be Uamp as well. The flow, temperature, and solute fields, as well as the melt/crystal interface, are also represented in the Cartesian coordinate (x, y, z). The dimensionless variables are defined by scaling length with the ampoule length L, velocity with O~m/L, pressure with Pm °~m/L 2 2, and solute concen-
M.C. Liang, C.W. Lan / Journal of Co'stal Growth 167 (1996) 320-332
322
the Prandtl number, Sc the Schmidt number, and 0m the melting temperature, which also serves as a reference temperature. Two important dimensionless variables, Ra s and Ra T, in the source term of the momentum equation are the solutal and thermal Rayleigh numbers, defined as follows:
i
/
g/3s CoL3 Ra s =
g/3T A TL3 ,
Ra T -
Ofm Pm
, O~m/'Pm
where g is the gravity acceleration, /3s and /3v are the solutal and thermal expansion coefficients, respectively, AT = Th - T~, and vm is the kinematic viscosity. The gravity direction eg can be decomposed into the Cartesian components (ex,ey,e ~)
,, °°"i/
eg = cos( a ) e x - c o s ( / 3 ) e y -
/
Fig. 1. Schematic sketch for the vertical Bridgman (VB) crystal growth system.
cos(y)e~,
where the o~, /3, and y are the angles between the gravity direction and the axes of x, y, and z, respectively. For a perfectly aligned ampoule, o~ = / 3 = ~-/2, and y = 0. In the present study, nonaxisymmetry is introduced only by tilting the ampoule in the y-z plane, i.e. a = ~-/2. In the crystal (c) and the ampoule (a), only heat transfer needs to be considered,
-Peie :. VO= KiV20 ( i = c , a ) , tration with its average concentration in the crystal, C 0, where a m is the thermal diffusivity and Prn the melt density. The dimensionless temperature ( 0 ) is defined as
O(x,y,z)=[T(x,y,z)-T~]/(Th-T~),
(1)
where T h and Tc are the top (hot) and bottom (cold) temperatures, as shown in Fig. 1. The pseudosteady-state governing equations describing convection and heat and solute transport in the melt (m) are as follows: V ' v = O, v" Vv = - V P + P r g 2 v + P r [ R a s ( C - R a T ( 0 - - 0m)] eg,
(2) 1) (3)
v" VO= V20,
(4)
Pr v-VC=--V2C, Sc
(5)
where v, P, and C are the dimensionless velocity, pressure, and solute concentration, respectively. Pr is
(6)
(7)
where Pe i =-PiCpiUampL/km is the Peclet number and Ki = k i / k m the dimensionless thermal conductivity of crystal or ampoule; k m is the thermal conductivity of the melt. Also, P i , Cpi, and k i a r e the density, specific heat, and thermal conductivity of the phase i (i = m, c, or a), respectively. The no-slip condition is used for the velocity on solid boundaries b' = -- Tct;ampez,
(8)
where Yc -= Pc/Pm and Uamp is the dimensionless ampoule pulling speed. Since the pseudo-steady-state is assumed, the upper open boundary is considered as an artificial boundary [5], and its velocity boundary condition is the same. The thermal and solutai boundary conditions at the melt/solid interfaces are set by the heat and mass flux balances. For example, at the growth front:
n.VOlm--n.K~VOL+PecSt(n.e:) = 0 , -n'VCIr.
(9)
Pe m Sc
--(1-K)C(n.e~)=O, Pr
(10)
M.C. Liang, C. W. Lan / Journal of Crystal Growth 167 (1996) 320-332
where n is the unit normal vector at the growth interface pointing to the melt. The Stefan number St = AH/Cp~AT scales the heat of fusion ( A H ) released during solidification to the sensible heat in the solid. The Pem = p~CpmUampL/km is the Peclet number and D the solute diffusivity in the melt. Also, K is the equilibrium segregation coefficient of the solute according to the phase diagram; K =- C J C at the growth interface, where C¢ is the solute concentration in the crystal. The temperature at the melt/solid interface is assumed to be the equilibrium liquidus temperature of the material. Solute diffusion in the solid is neglected here. Temperature at the top and bottom surface is set to the furnace temperature there. With the pseudo-steady-state assumption, an artificial boundary condition is used at the upper surface for the consistence of the overall solute balance [5]:
-n'VClm-
Pe m Sc
--(CPr
1)(n .e=).
ing. The complications of computing view factors from the sample to the furnace are neglected here.
3. Finite-volume/Newton method Owing to the unknown melt/crystal interface shape hc(x,y,z), the physical Cartesian coordinates (x, y, z) are transformed into general (nonorthogonal) curvilinear coordinates ( sc 1,~: 2, sc3), which fit all the interfaces. Fig. 2a is a sample converged mesh for computation; the mesh defines the edges of finite volumes. It should be pointed out that the interface hc(x,y,z) is unknown a priori and needs to be solved simultaneously with other field variables. During iteration, once the interface is updated, the mesh can be easily generated through an algebraic boundary-fitted coordinate transformation, and grid density control is through the stretch functions. The
(11)
The heat exchange between the ampoule and the furnace is by both radiation and convection according to the energy balance along the ampoule surface, - n . KaV0la= B i ( 0 - 0a) + Rad(0 4 - 04),
Ampoule
(12)
where n is the unit normal vector on the ampoule surface pointing outwards, B i - hL/k m the Biot number, and Rad - O'EaAT3L/km the radiation number; o- is the Stefan Boltzmann constant, while e a is the surface emissivity of the ampoule. For simplicity, the effective furnace temperature 0a is assumed to be 0a(x,y,z) = Z -- 0~( x, y,z),
Melt
(13)
3rowth nterface
where, Od(x,y,z) is a deviation temperature to describe thermal nonuniformity. A Gaussian distribution is used for 0d,
O,~(x,y,z)=AOdexp(-[(z-0.5)/a*]2},
323
Crystal
(14)
~ . ~
Y
x
and
= a0 p cos(
),
(15)
where A 0dp is the dimensionless peak (maximum) deviation temperature difference along the azimuthal angle q~ and a* a dimensionless parameter for the width of the distribution. In the present study, the nonaxisymmetry for a perfectly aligned ampoule is introduced by a nonzero A 0ap for nonuniform heat-
Ca)
(b)
Fig. 2. (a) 3D finite volume mesh for computation; a part of the mesh is removed to illustrate the finite volumes inside; (b) a typical finite volume.
M.C. Liang, C.W. Lan /Journal of Cr3,stal Growth 167 (1996) 320-332
324
transformation is similar to that used in the 3D thermal-capillary simulation of floating-zone crystal growth, and that has been described elsewhere [19]. The governing equations and boundary conditions can then be approximated by a finite volume method. As shown in Fig. 2a, the physical domain is divided into a number of finite volumes bounded by cell faces. A sample of the finite volume is illustrated in Fig. 2b. In general, each volume is an irregular hexahedron defined by eight vertices. The coordinate invariant strong conservative-law form of the transport equations (Eqs. (2)-(5) and (7)) can be written
where W is a nonzero weighting function. This isotherm approach without the weighting function has also been used by Ettouney and Brown [21] and Lan and Kou [22]. The purpose of using the weighting function here is to avoid the numerical breakdown during LU decomposition in the solution stage, since no pivoting is used there. A linear combination of the interface variables near the point of interest is used for the weighting function. Eq. (17) and the above weighted isotherm condition yield a system of nonlinear algebraic equations, which can be written symbolically as
as
f ( x ) =-f(v,P,r,C,hc)
V.I+S=0,
(16)
where I describes the flux vectors of the mass, momentum, energy, and solute. The source term S consists of the body forces for the momentum equation, but it is zero for other equations. If Eq. (16) is integrated over a general control volume AV in physical space and the Gauss divergence theorem is then applied, a discrete flux balance equation can be obtained, i.e.,
..~ I . A I [e + l - a l l w + l ' a 2 [ n + l ' a 2 1 ~ + / ' a 3 1 t
+I.A3Ib + SAV=O,
(17)
where A i, i = 1 to 3, is the area vector defined on the cell surfaces (see Fig. 2b). Working on the control volumes, one could derive the flux components using the central difference scheme. However, the mass fluxes used in the equation of continuity need to be taken care of separately to avoid the checkerboard pressure oscillation due to the velocity/pressure decoupling. A modified momentum interpolation scheme [20] is used in the present calculations. The detailed description of the numerical issues concerning the primitive variable formulation on the collocated grid for incompressible flow is described elsewhere [20]. The boundary conditions can be imposed easily during the finite volume approximation of the cells next to the boundaries. To the growth interface hc(x,y,z), a weighted melting-point isotherm approach is adopted
W[ O(x,y,hc) -
0re(C)] = 0,
(lS)
(19)
This nonlinear equation set is solved simultaneously by Newton's method for all variables. Starting from an initial approximation to vector of unknowns x °, successive updates are constructed as x , + l = x " + 6 n+l,
(20)
and the correction vector 6 "+1 is the solution of the linear equation set J6 "+1:
fAl'dA + f±vSdV
= 0.
-f(x").
(21)
The components of the Jacobian matrix J, formed by explicit differentiation as f j - Ofi/Oxj, represent the sensitivity of the residual vector to the solution vector. They are all determined by the forward difference [23,24]. Since Eq. (21) is large and sparse, the use of an efficient sparse matrix solver is necessary. The ILU(0) [25] preconditioned GMRES method [26] is used, and it is found effective and efficient in the computation here.
4. Results and discussion
For comparison purposes, we consider the gallium-doped germanium (GaGe) growth in the Grenoble furnace investigated by Adornato and Brown [5] in this study. The same system has also been studied by Lan and Chen [12] as a benchmark problem for comparing the performance of a 2D FVM (stream function-vorticity ( q - w ) formulation) and a 2D Galerkin FEM (primitive variable formulation). The physical properties of GaGe and some input parameters used in the study are listed in Table 1. The condition of Ra s = Rad = 0 used in Ref. [5] is also
M.C. Liang, C. W. Lan / Journal of Crystal Growth 167 (1996)320-332
[5,12] are tested again. Since the computation is performed in the ( x , y , z ) coordinates rather than the cylindrical coordinate (r,z), and no axisymmetry is assumed, fully 3D computation is still necessary. Fig. 3a-3c illustrate the 3D calculated results showing the effect of convection for the axisymmetric growth ( y = 0 and A0dp = 0). In each figure, the velocity field on the y - z plane of the system is illustrated on the left, and the solute field on the right. The solute distribution (on the x - y plane) in the grown crystal is illustrated underneath the solute field of the y - z plane. As shown, all solutions are perfectly axisymmetric with respect to the z axis as expected. Without the buoyancy force (Ra T = 0), as shown in Fig. 3a, the flow in the melt is due to ampoule translation, and the solute transport is mainly through diffusion. The radial segregation shown on the x - y
adopted here for comparison. In this system, the solute concentration is so low that the solutal effect is ignored. The solutal effect of a nondilute alloy can be considered if necessary [12]. The mesh used in the following calculations is similar to Fig. 2a. It contains 29 176 finite volumes leading to 111 572 nonlinear equations. All computations are performed in the IBM/ES9000 supercomputer, and one Newton iteration takes about 48 min CPU time (only one CPU is used for computation). The convergence of Newton's iterations is quadratic in all cases, and it takes only four to five iterations to get a convergent result with the L~-norm of 1 × 1 0 - 9 . 4.1. Axisymmetric solution and comparison with 2D cases To validate the code by different numerical approaches, the axisymmetric cases studied before RaT=
RaT = 0 [~-.-
325
RaT= 10 6
10 5 ....
! i
I
....
#Jt,ILJ.ttata.t.,lltili.~U JamUJ~:H,ILt~l.l.ti~.m~ laLm:;~,k~m~HhJl~:l,.illl I¸
. . . . . . . .
I
!
. . . . .
i
!
!:::: :,5: i;ii ii; :¸ ::ii
x Cm.~ = 1 2 . 3 9 *
C.~
+ Cc~ Cc...
* Cmt..~ = 1 . 0 2 7
1.122
+ CCm,~ = 1 . 2 2 9
- Cc,~., = 0 . 8 6 2
Cemt. = 0 . 8 5 4
*
= 1.078
+ Cc~=
0.880
Ca)
~
= 1.012
= 1.011
=
1 -- I 0 0 c m / s x Cm~ = 1 4 . 1 3
i ---~ = 5 0 0 c m / s x C,~¢ = 1 2 . 8 9
1
C.am
(b)
Co)
Fig. 3. Effect of convection on the flow and solute fields for axisymmetric growth; (a) Ra T = 0; (b) Ra w = 1 0 - 5 ; (c) Ra T = 10 6 in each figure, the left-hand side is velocity field and the right-hand side the solute fields. The upper part of the solute fields is taken from the ( y - z ) plane, while the lower one from the (x-y) plane.
M.C. Liang, C. W. Lan / Journal of Crystal Growth 167 (1996) 320-332
326
Table 1 P h y s i c a l p r o p e r t i e s a n d i n p u t p a r a m e t e r s [5,12]
GaGe p~ = 5.5 g c m - 3 Pm = 5.5 g c m - 3 Tm = 9 3 7 . 4 ° C ~H = 460 J g- i k~=0.17Wcm
I oc-I
k m=0.39Wcm
I °C
t
Cpc=Cpm=O.39Jg
i oC
/z= 0.00715 g cm-I s ]3T=5X10 4o C I
I
,8 s = 0 ( t o o l % G a ) D=2.1XI0 4cm
I
t 2/s
K = 0.087
Graphite (ampoule) p ~ = 1.8 g c m
3
k~ = 3.26 W c m - i ° C - i Cpa=l.814Jg
- l oC - I
ea = 0
Other input parameters L=7cm R c = 0.5 c m
R,
= 0.7 c m
T h = 1112.4°C T~ = 7 6 2 . 4 ° C Uam p = 4 X 10 - 4 c m / s h=46.571Wcm 2oc
I
here. Fig. 4 illustrates the comparison of the radial (v r) and axial (G) velocities on y - z plane with the ones from the FEM for Ra T = 1 0 6. As shown, they are in excellent agreement. The comparison of the calculated growth interfaces is shown in Fig. 5a for Ra T = 1 0 6. To facilitate the comparison, the scale in the vertical axis is enlarged, so that the small differences can be seen. As shown, they are also in good agreement. More importantly, an excellent agreement is also obtained when the calculated radial solute segregations by the three numerical approaches for different Ra T numbers are compared in Fig. 5b. The averaged solute concentration should be equal to one, which is often used for checking the global solute balance. For both 2D and 3D FVMs, the global balance is satisfied exactly. However, very small global errors still exist in the FEM. For an even larger Ra T number, it requires a finer mesh for the 3D computation to get an accurate result. Unfortunately, it is too costly to do so because of the limitation of present computing resources. From previous studies [5,12], the maximum radial solute segregation occurs at Ra T = 1 0 6 for the GaGe system. A further increase in Ra T leads to a better mixing and thus a smaller solute segregation. Accordingly,
a=lcm y = 0°_1.5 ° ATdp = 0 - 1 ° C
plane is caused by interface deflection. At Ra T = 105, small convection loops near the growth front are induced. Basically, they are induced by radial temperature gradients due to interface deflection. The solute field is also slightly distorted by the convection. As Ra T is further increased t o 1 0 6 , the convection becomes stronger, and the maximum velocity is now much larger than the growth rate (4/xm/s). As a result, the solute field is highly distorted, and the radial segregation becomes larger. For the convenience of further discussion, the degree of the radial solute segregation is defined as the maximum solute concentration difference in the crystal, i.e. (Ccm,x C ..... ) X 100%. Since the solutions from the 2D Galerkin FEM and the 2D FVM [12] are also available for the axisymmetric cases, a detailed comparison of our 3D computation with these two approaches is also made
3D FVM x maX(Vz}= 0 . 0 0 8 5 1 ( c m / s )
Vz
Vr
• min(vz)=-0.00372 + maX(Vr)= 0 . 0 0 1 2 3 - min(vr)=-0.00347
.......
2D FEM
x maX(Vz)= 0 . 0 0 8 7 4 ( c m / s )
f
* min(vz)=-0.00380 + max(vr}= 0 . 0 0 1 3 8 - min(vr}=-0.00368
f
F i g . 4. C o m p a r i s o n o f 3 D c a l c u l a t e d radial ( v r) a n d axial ( t , : ) v e l o c i t y f i e l d s w i t h the o n e s b y 2 D F E M .
M.C. Liang, C.W. Lan / Journal of Crystal Growth 167 (1996) 320-332
the effect of 3D flows on the solute segregation is also expected to be more significant in this region. The axisymmetric results for higher Ra T by the 2D FVM and FEM can be seen elsewhere [12]. 4.2. Tilted ampoule In practice, a perfect ampoule alignment with the gravity may not be possible. During crystal growth, a 3.5
¢
i
i
i
i
i
b
i
i
a
3.48
. . . . . . . . 2D-FVM(~¢-o~} -........ 2D-FEM
/ I
/
3.46
~
3.44
3.42 3.4 3.38 3.36
[
I
[
i
0.1
[
[
0.2
i
i
0.3
I
0.4
0.5
r (cm)
1.4
i
!
i
~
i
t
i
i
r
b
1.3 / R a T -- 1 0 6 1.2
~
~
/
RaT =105 Cc 1.1 RaT = 0
1
...... 09
~
_
2D_VVM¢ ,-,o) I
.......
0.8
i
0
i
0. I
i
i
i
i
0.2 0.3 r (cm)
l
l
0.4
l
0.5
Fig. 5. Comparison of 3D calculated results with those by the 2D FVM and 2D FEM; (a) interface shape (for Ra T = 106); (b) radial solute segregation.
327
small tilt (less than 1°) from the gravity could be very common, and this may break axisymmetry. Fig. 6 shows the effect of tilt angle for Ra T = 105. As shown in Fig. 6a, even though the tilt is only 0.5 °, the axisymmetry of lower cells is broken, and the convection becomes asymmetric. The lower right cell grows larger, and the lower left cell shrinks slightly, while a new big flow cell develops in the main body of the melt. The fluid flow also becomes faster near the right ampoule wall, and slower on the other side. As a result, the solute field is distorted slightly by the flow and becomes asymmetric. The maximum solute concentration (Ccmox) shifts slightly from the center to the left, and the radial solute segregation increases to 31.4% as compared with 26% in Fig. 3b. As the tilt angle is increased to 1.5 ° in Fig. 6b, the new flow cell becomes stronger. The lower left cell becomes much smaller, and the right cell grows much bigger and joins the big new flow cell. The fluid flow near the right wall is also enhanced. As shown, the flow is now fully 3D, and is not possible to be approximated by an axisymmetric model. It should be pointed out that the flow and solute fields in the x - z plane are still symmetrical. Due to the 3D fluid flow, solute transport becomes 3D as well. The solute field is distorted even more along the main flow direction, and the segregation is increased to 39.8%. Cc~x is also pushed to the left even more. As the convection becomes stronger (Ra T = 106), as shown in Fig. 7, the effect of tilt on the flow structure and hence the solute segregation becomes even more significant. Even a tilt of only 0.5 ° has a dramatic effect on the flow and solute fields (Fig.7a). As shown, the solute field is now highly distorted by the convection. The affected area now extends to the upper part of the melt as well. Consequently, the radial segregation is increased to 61.5% (39.5% in the axisymmetric case). As the tilted angle is further increased, as shown in Fig. 7b, the convection becomes even stronger, and Ccm,~ shifts to the left even more. However, the radial segregation is reduced to 54.5% due to better mixing. Apparently, the 3D flow promotes the better mixing of solute in the bulk melt and leads to a dramatic decrease in the concentration difference. Cmi n is increased to 4.175 as compared with 1.027 for the axisymmetric case in Fig. 3c. It is expected that to some extent the solute segregation
M.C. Liang, C. W. Lan /Journal of Crystal Growth 167 (1996) 320-332
328
may be reduced at higher Ra T and larger ampoule tilting [ 17]. Furthermore, in the flow visualization experiments by Neugebauer and Wilcox [14], it was observed that 3D flows were dominant in the melt. In their experiments, the furnace was specially designed to achieve a good thermal axisymmetry. From the implication of Figs. 6 and 7, the 3D flows might be caused by a slightly tilted ampoule used in the experiments. The melt height in their experiments happened to be much larger than the diameter, which was similar to the cases here. Such a high aspect ratio may also promote the formation of 3D flows
7=
0.5 °
RaT=
due to stronger convection. In other words, for VB crystal growth 3D flow and solute fields could be induced even more easily with increasing melt aspect ratio. Nevertheless, we do not exclude the possibility of the imperfect furnace as the cause of 3D convection. The effect of nonuniform heating due to an imperfect furnace is discussed next. 4.3. Nonuniform heating
In a VB furnace, no matter how good the design is, small temperature gradients in the azimuthal direction may still exist. If the ampoule is not rotated,
10511
7 = 1.5 ° !
:,T~11!,,,
I..~., ....
"~J;'il
I
lit!
::ji]~l~'iits,~ I..,,,:: .... 't~'~*lllii;it11*~l ! ir'~ll'll,
!
,r'H'il, ,,rrll[~
J I
,
i iij1111
,
' 'r11'~1
c
............
"
I
i
L~
i
., ILt~
i~
I
1 -~ = 500 cm/s
1
= 300 x Cm~ = * Cr~n = + Ccm.x = Ccmm =
x Cm,x = 1 2 . 9 7 * Cmtn ~- 1 . 0 1 4 ae
Ccmax
=
1.128
Ccml~ = 0 . 8 1 4
cm/s
13.42 1.048 I. 168 0.738
(a) Fig. 6. Effect of tilt angle for Ra T = 105: (a) 3' = 0.5°; (b) 3' = 1.5 °,
(b)
~
M.C. Liang, C.W. Lan/Journal of C~. stal Growth 167 (1996) 320-332
which is quite common in practice, the temperature nonuniformity may be even worse. To consider this effect, the deviation temperature distribution, Od(X,y,z), is used to introduce a small temperature nonuniformity in the azimuthal direction. Fig. 8 shows the effect of nonuniform heating for Ra T = 10 5. In Fig. 8a, the maximum temperature difference is I°C (or ATap = 0.5°C); the overall furnace temperature difference, T h - Tc, is 350 ° here. As shown, the flow and solute fields have become asymmetric near the growth interface. Since the thermal nonuniformity is only imposed near the growth front, the flow far away from the interface is not affected much,
7 = 0.50
RAT=
which is different from the cases with a tilted ampoule. In each figure, the ambient temperature Ta is higher on the left but lower on the right introducing thermal gradients to drive the flow clockwise. As a result, the right cell becomes stronger, while the left cell weaker. As shown in Fig. 8a, the solute segregation is increased to 35% as compared with 26% in the axisymmetric case (Fig. 3b) and 31.4% in the tilted case (Fig. 6a). Apparently, at Ra T = 105 the I°C nonuniformity seems to have a larger effect than a half-degree tilt. With a larger temperature difference (ATdp = I°C), as shown in Fig. 8b, the effect becomes even more significant. The solute segrega-
10 6 ]
7 = 1.5 °
11ilt llJ:~
,,qU~]TII,, ""'~tilillliLL'
il.E~II'l~
\,
,rTTIT~TTTI~ ' " ! J I I I I H I ] ,1~ ,~THITtI~,, ' " J l l ! l J l l l l J ' ,,,TIT;1~I~,.. i " ' ~ l ; l l l l l l ' l ' 1
i
\
1
= -cm/s 100 x Cm= = 15.01
--~ =--
1
cm/s
50 x Cm~ = 1 5 . 0 3
* Cm~ = 1 . 9 6 8
b-
* Cm~ = 4 . 1 7 5
+ Ccm= = 1 . 3 0 6
+ Ccm= = 1 . 3 0 7
Cc~n = 0 . 6 9 1
Ccm~ -- 0 . 7 6 2
(a)
329
{b}
Fig. 7. Effect of tilt angle for Ra T = 106: (a) Y = 0.5°; (b) 7 = 1-5°.
M.C. Liang, C.W. Lan / Journal of Crystal Growth 167 (1996) 320-332
330
tion is increased to 46%; however, the flow and solute fields in the upper half of the melt are still not affected much. Again, it should be pointed out that 1 or 2°C difference is quite trivial as compared with the overall ambient temperature difference, which is 350°C here. However, quite unexpected, its effect on flow structures and solute fields is significant. For stronger convection (Ra T = 106), as shown in Fig. 9, the effect of nonuniform heating becomes even more significant leading a larger radial solute segregation. In Fig. 9a (ATdp = 0.5°C), the fight flow cell becomes much stronger than the left one, and the solute field is highly distorted along the flow direc-
ATdp =
tion near the interface. As a result, the solute segregation is increased to 63%, as compared with 37.5% in the axisymmetric case (Fig. 3c) and 61.5% in the tilted case ( y = 0.5 ° in Fig. 7a). With a larger ATop of I°C, the flow becomes fully 3D at the bottom, and the left flow cell almost disappears. The flow becomes stronger as well. Due to the much stronger convection, the solute field near the interface is distorted even more; CC~a~ is also pushed to the left even more. Interestingly, the solute segregation is decreased slightly (57.8%) as compared with Fig. 9a. It should be pointed out that the maximum radial solute segregation usually occurs at an intermediate
IRa T- 10 5 ]]
0.5°C
ATdp= 1.0°C
~:llllt ] i I '. I ll:]~l:till I I ~ / I l:lltlt
--
"
~-
. . . . . . . . . . . i i i i i ! i i i i i ...........
I
I
\
1 = --- cm/s 300 x Cm~x = 1 3 . 5 5 * C~. = 1.016 + Cc~.~ = 1 . 1 7 9 - Cc~ = 0.725
I
--~ = 500 cm/s X Cmax ~- 1 3 . 0 9 *
Cmin
=
1.013
4"
Ccmax
=
1.138
Ccmm = 0 . 7 8 6
{a)
Co)
Fig. 8. Effect of nonuniform heating for Ra T = 105: (a) ATdp = 0.5°C; (b) ATdp = I°C.
M.C. Liang, C.W. Lan / Journal of Crystal Growth 167 (1996) 320-332
ATdp = 0.5°C iii I.
ii .
.
RAT=
10 6
331
= 1.0°C
ATdp
iii'. !
.
I
[ i L
i¸ - - - - - _ J ......... . . . . . . . . . . .
i............ ] .............
itlz ........ "'l
...........
1 ............. ....... L, ....
F''''XX,
,~hL,"! , ,,r ,, ,,1, t/ .......
~',~X~k~tL"
.... :r'e'elt//t/e'r~\\\\,
....',rYittt.,,...'? "X~II~L'k " , ,rll l; 'rt ,~. "'x\\'ll~llk~"]
"'\\\ ~l',ili~' I
,11[! I [ T I ~ lll]]/~.
' ii.
.,,'11 J~J]JJ ; "'~\~,'~1 I11~' ,l~,~ ~, "' i rl+ I I: ?,
1
--~" = ~ -
= 100 c m / s Cm~
1
*
em/s
x Cm~x = 15.22 • C~ = 1.172 + C e ~ , = 1.325 Ce~n = 0.74"7
15.24 * Croan = 1.092 + C e ~ = 1.326 Cc~, = 0.696
X
''C s
=
(a)
(b)
Fig. 9. Effect of n o n u n i f o r m heating for R a T = 106: (a) ATdp = 0.5°C; (b) AT@ = I°C.
flow intensity (or poor mixing) [5]. Too weak (without mixing) or too strong (well mixing) a convection results in less segregation. Therefore, at the maximum solute segregation, further increasing the flow intensity will result in a better mixing and thus less segregation. It is clear that the flow and solute fields become 3D, even though the maximum temperature deviation from uniform heating is only 1-2°C. The induced 3D flows are stronger than the 2D ones and enhance solute transport significantly, as shown by the much distorted solute fields in Figs. 9a and 9b. Such a significant effect has not been illustrated before. Nonaxisymmetric flows due to asymmetric heating
were also observed by Pulicani et al. [16]. However, because they used higher axial thermal gradients near the growth interface, their flow patterns were less sensitive to the azimuthal temperature deviation. In the present system, if Ra T is further increased, the system will be even more sensitive to the nonuniform heating. In other words, in practice a perfect axisymmetric growth could be difficult to achieve.
5. C o n c l u s i o n s
A three-dimensional finite-volume/Newton's method is developed to study the effect of slight
332
M.C. Liang, C. W. Lan / Journal of Crystal Growth 167 (1996) 320-332
nonaxisymmetry on the fluid flow and solute segregation in the vertical Bridgman crystal growth. The effects of a slightly tilted ampoule and a nonuniform thermal profile are studied. From the calculated results, it is found that slight deviations from axisymmetry can easily induce three-dimensional flow and solute fields, but their effect on the interface shape is small. As convection increases, the system becomes even more sensitive to such deviations. Compared with the axisymmetric cases, the induced three-dimensional convection is much stronger and causes a larger radial segregation.
Acknowledgements The authors are grateful for the support from the National Science Council and the National Center for High-performance Computing of the Republic of China under Grant No. NSC84-2215-E008-011.
References [1] W.A. Gault, E.M. Monberg and J.E. Clemans, J. Crystal Growth 74 (1986) 491. [2] K. Hoshikawa, H. Nakanishi, H. Kohda and M. Sasaura, J. Crystal Growth 94 (1989) 643. [3] E.M. Monberg, W.A. Gault, F. Simchock and F. Domingguez, J. Crystal Growth 83 (1987) 174. [4] E.M. Monberg, in: Handbook of Crystal Growth, Vol. 2a, Basic Techniques, Ed. D.T.J. Hurle (North-Holland, Amsterdam, 1994) p. 52. [5] P.M. Adomato and R.A. Brown, J. Crystal Growth 80 (1987)
155. [6] S. Motakef, J. Crystal Growth 97 (1989) 20l. [7] M.J. Crochet, F. Dupret, Y. Rychmans, F.T. Geyling and E.M. Monberg, J. Crystal Growth 97 (1989) 173. [8] R.A. Brown and D.J. Kim, J. Crystal Growth 109 (1991) 50. [9] S. Branton and J.J. Derby, J. Crystal Growth 121 (1992) 473. [10] D. Hofmann, T. June and G. Muller, J. Crystal Growth 128 (1993) 213. [11] C.W. Lan and C.C. Ting, J. Crystal Growth !49 (1995) 175. [12] C.W. Lan and F.C. Chen, Comput. Methods Appl. Mech. Eng., in press. [13] F. Dupret and N. Van den Bogaert, in: Handbook of Crystal Growth, Vol. 2b, Growth Mechanisms and Dynamics, Ed. D.T.J. Hurle (North-Holland, Amsterdam, 1994) p. 875. [14] G.T. Neugebauer and W.R. Wilcox, J. Crystal Growth 89 (1988) 143. [15] T. June, PhD Thesis, University of Erlangen-Nurnberg, 1993. [16] J.P. Pulicani, S. Krukowski, J.I.D. Alexander, J. Ouazzani and F. Rosenberger, Int. J. Heat Mass Transfer 35 (1992) 2119. [17] S. Kuppurao, Q. Xiao and J.J. Derby, AIChE 1995 Annual Meeting, Miami Beach, FL, USA. [18] P.M. Adornato and R.A. Brown, Int. J. Num. Methods Fluids 7 (1987) 761. [19] C.W. Lan and M.C. Liang, Int. J. Num. Methods. Eng., to be published. [20] M.C. Liang and C.W. Lan, J. Comp. Phys., accepted. [21] H.M. Ettouney and R.A. Brown, J. Comp. Phys. 49 (1983) 118. [22] C.W. Lan and S. Kou, J. Crystal Growth 108 (1991) 351. [23] A.R. Curtis, M.J.D. Powell and J.K. Reid, J. Inst. Maths. Appl. 13 (1974) 117. [24] T.F. Coleman, B.S. Garbow and J.J. More, ACM Trans. Math. Softw. 10 (1984) 329. [25] J.A. Meijerink and H.A. van der Vorst, Math. Comp. 31 (1977) 148. [26] Y. Saad and M.H. Schultz, SIAM J. Sci. Stat. Comput. 7 (1986) 856.