COMPUTER METHODS IN APPLIED MECHANICS AND ENGINEERING 12 (1977) 145-152 0 NORTH-HOLLAND PUBLISHING COMPANY
DIGITAL SIMULATION OF FLOW PATTERNS IN THE CZOCHRALSKI CRYSTAL-PULLING PROCESS W.E. LANGLOIS and C.C. SHIR IBM Research Laboratory, San Jose, CA 95193,
USA
Manuscript received 14 September 1976 Revised manuscript received 3 February 1977
A procedure is described for the numerical solution of the time-dependent equations governing the flow in a crystal-growth crucible. The molten crystal material is taken to be a Newtonian fluid with constant viscosity. The Boussinesq approximation is employed, i.e. compressibility is neglected except insofar as thermal expansion modifies the gravitational body force. Crucible rotation is the main driving force, but a series of numerical experiments suggests that buoyancy effects can modify appreciably the impurity-blending estimates obtained from an isothermal study.
1. Introduction In the Czochralski technique the growing crystal is pulled very gradually from the center of a rotating crucible of melt, as shown in fig. 1. The growing crystal is also rotated with a different angular velocity. The rotation may be steady or time-dependent, and the vertical motion of the crystal is negligible in comparison. The structure of the crystal is strongly influenced by flow patterns within the melt. These are generated and shaped by thermal convection as well as by the centrifugal pumping effect of the rotating solid boundaries. Hence, the rotary flow about the axis of symmetry is coupled to a secondary circulation in the meridional planes. Computational results for steady-state flow patterns have been reported by Kobayashi and Arizumi [ 11, [ 21. The present paper describes some early results obtained with a time-dependent model, motivated in part by current interest in the accelerated crucible-rotation technique [ 31, [ 41. So far, however, we have only examined flows driven by constant crucible and crystal rotation. These evolve to a steady state, and the final flow patterns can be compared with those obtained by Kobayashi and Arizumi. The melt is taken to be a Newtonian fluid with constant viscosity. We employ the Boussinesq approximation, i.e. compressibility is neglected except insofar as thermal expansion modifies the gravitational body force. Thus, in the interior of the melt the motion is governed by the Navier-Stokes equations specialized to the case of buoyancy-driven axisymmetric flow. Following well-known procedures, these can be transformed to prognostic equations for the temperature, secondary-flow vorticity and azimuthal velocity, along with a diagnostic equation for the secondary-flow streamfunction. On the solid boundaries the azimuthal velocity of the fluid matches that of the wall, the vorticity is adjusted such that the secondary flow vanishes, and the temperature is specified. Appropriate conditions for axial symmetry are specified along the common axis of rotation of crucible and crystal. The free surface of the melt requires special consideration. Its presence should, in principle, preclude use of the vorticity-streamfunction formulation because the normal stress boundary con-
146
W.E. Langlois and C.C. Shir, Digital simulation of flow patterns in the Czochralski crystal-pulling process
__
Fig. 1. Geometry of the Czochralski technique.
dition involves the pressure. However, in a Czochralski crucible the free surface remains nearly flat *. It has been shown [ 51 that the vorticity-streamfunction method can then be used after all with the vorticity set to zero on the nominal level of the free surface. To satisfy the tangential stress condition, the normal derivative of the azimuthal velocity is also set to zero. A temperature boundary condition is obtained by assuming that the free surface Ioses heat via Stefan-Boltzman radiation.
2. Equations of motion
To increase the generality of the results, we employ dimensionless variables. However, a detailed scale analysis will be deferred to a later paper since it entails two complicating factors. First, the proper velocity scale for the azimuthal velocity component is usually much larger than that of the meridional components. Second, the radiation boundary condition is expressed in terms of the absolute temperature, whereas the buoyancy depends on temperature gradients within the melt. * In practice, surface tension permits the growing crystal to be raised somewhat, giving the free surface a meniscus. Our slightly idealized geometry leads to a flow field - which we might term Czochraiski bufk flow - that is incorrect in the immediate vicinity of the meniscus.
W.E. Langlois and C.C. Shir, Digitalsimulation of flow patterns in the Czochralskicrystal-pullingprocess
147
Hence, for the present we simply use the units of measurement as scale factors. We denote the units of time, length and temperature by T, h and 0, respectively; it will not be necessary to introduce a unit of mass. To illustrate by example: suppose that the crucible radius R,h is 10 centimeters, and suppose further that the measuring system is MKS, so that A = 1 meter = 100 centimeters; then R, is the dimensionless number 0.1. All variables in subsequent formulas are dimensionless and are related to the physical variable according to the transformations
w w*
+B,
- K ‘KY x2/r
v+E, x2/r
0)
where (Yis the volumetric expansion matic viscosity, and K is the thermal In order to eliminate the pressure with Boussinesq approximation, we a Stokes streamfunction J/, 2.4=
(l/r)a+/az
)
w =
coefficient, g is the gravitational acceleration, v is the kinediffusivity. from the Navier-Stokes equations for axisymmetric flow express the meridional velocity components u, w in terms of
-(l/r)a$/ar
)
(1)
and introduce a vorticity w defined by w = aw/ar
- au/az
(2)
.
It follows directly from (1) and (2) that the streamfunction diagnostic equation
a s
(ia$ 1 Tar
ia2+
+;--+
and vorticity are related through the
(3)
=-cd.
The Navier-Stokes equation and the energy equation with viscous dissipation neglected yield prognostic equations for the vorticity, azimuthal velocity u, and absolute temperature T. With the dimensionless transformations we get
(44 (4b) gy+ia$
(ruT)+:
(wT)=KV2T,
(4c)
where V2 = a2/ar2 + (l/r)a/ar
+ a2/az2
.
148
W.E. Langlois and C.C. Shir, Digital sinudatiorz
of flow
patterns in the Czochralski crvstal-pulling process
Finite difference analogs of (4) are used to advance the vorticity, azimuthal velocity and tempotrature fields from one time step to the next. To describe these analogs, let the symbol 1 signify any of the dependent variables w, u, T or $. Further, let /,r, denote $[(i - l)Ar, Cj - I)Az, t] t” denote ,! [(i -- 1) Ar, (j -- 1) AZ, t + At] . The finite difference prognostic equations and let f used are ti;n
~kt’ 1,I = Co{ d[wfj
+ n(Ttj)
~ C(U;~: 1); 11 ;l> .
Ups’= cD{~[Ulli - c(~~j; 3); II : l} ) TiTit’ = CD{~[T~j;Ol with the operators
:O} ~
gQ, 9,
C, Cg (advective,
At sQ ‘[i,j:
cy3( f
i,j)
‘I) ={i.j
=
e
+ @ArAz
(fi+l,
j
-fi-1,
CfJ(fi,j;“)=fi,j +c:2
buoyant,
centrifugal,
diffusive)
defined
by
1 r__
plAr/zr
j) ’
[fi,j-1
-
2fi,j
+
fi,j+lI
At each time step the stream function is determined by an iterative solution of a finite difference analog to the diagnostic eq. (3). Applying the SOR (successive overrelaxation) method gives the following relation between successive iterates $” and $” +I :
where
W.E. Langlois and C.C. Shit+,Digital simulation of ,ilow patterns in the Czochralskicrystal-pullingprocess
149
i-l
i-l
(i-1)+1/2 + (i--1)--1/2 and 0 is a relaxation factor. At present we use fi = 1.6
3. Boundary conditions The region of computation is the rectangle 0 d YG R,, 0 < z < H. Since there is no flow across the periphery of this rectangle, it is a contour along which the streamfunction is constant. Since an additive constant does not affect its significance, we may take $ = 0 on all boundaries. On all solid boundaries the azimuthal velocity of the fluid equals that of the boundary. Since the free surface is stress-free, au/az vanishes there. Finally, u itself vanishes on the axis of symmetry. As shown in [S] , it is a reasonable approximation to take w = 0 on the almost-flat free surface. Since u and aw/ar both vanish on the axis of symmetry, so does w. On the solid boundaries, both meridional components of velocity vanish. The Taylor series expansion of the streamfunction then provides a boundary condition for w. For example, near the crucible bottom $4, AZ, f> = $(r, 0, t) + AzJ/,(r, 0, I) + (1/2)Az2 &(r,
0, t) + O(Az3).
By (1) and (2), along with the requirement that U, w and $ vanish along z = 0, neglecting the terms of order Az3 then yields 2 w(r, 0, t) = - $(r, AZ, 0 . r(Az)2 The temperature is specified on all solid surfaces. On the axis of symmetry, a T/h vanishes. The radiative heat loss from the free surface modifies the estimate of a T/at obtained from (4~) for those grid cells lying along the surface. In summary then the conditions obtaining on various segments of the computation region periphery are as follows: AXIS
r=O (i=l)
$=O,
o=O,
u=O,
aT/ar=O,
CRUCIBLE BOTTOM z = 0 (j = 1)
J/=0,
a=-2
$(r,Az,t),
u=d&,
T=T,,
r(Az)2
CRUCIBLE WALL r=Rc
J/=0,
a=-
2 R, CArI2
(i=M)
W,
- Ar, z, f),
v=R,S&,
T=T,,
150
W.E. Langlois
CRYSTALFACE
$=O,
and C.C. Shir, Digital simulation
of flow patterns
in the Czochralski
$(r,H-Az,t),
v = rfis ,
T= T s ’
am
$=o,
process
z=H,r~R,(j=N;i=1,2,...i,-11)
w=-2
FREESURFACE
crystal-pulling
z=H,R,
~=o,
(j=N;i=i,,i,+
av/az=o,
l,...M-
1)
aT/at=(arjat),,,.,-2e;~T4
1 ”
where E is the surface emissivity, and the dimensional parameters p, c, and u are, respectively, the melt density, the melt specific heat and the Stefan-Boltzman constant. The factor 2 arises in the radiation term since the grid lying in the surface correspond to grid cells of height AZ/~ rather than AZ.
4. Comparison with a result obtained by Kobayashi and Arizumi To test the program, we ran it using physical and geometrical parameters corresponding to the most complicated flow reported by Kobayashi and Arizumi [21, viz. the case they label “Rossby number-4”. If we choose the units of length and time so that R, = 1.0 and a2, = -1.0, then the other parameters for this case have the values R, = 2.5, H = 5.0, as = 3.0, E = 0.05, B = K = E = 0 (isothermal flow). We ran this case using a 39 X 97 grid (Ar = 0.06579, AZ = 0.05208) and a timestep of 0.001. Beginning with an initial guess for the J/ and o fields, we continued the computation until the flow settled down to a steady state. Fig. 2 illustrates the final $-field, i.e. the streamlines of the meridional flow at steady state. The pattern is quite similar to that obtained by Kobayashi and Arizumi using their technique for integration of the steady-state equations.
5. Nonisothermal
flows
We carried out a series of computations arbitrary units, we fixed R,=4.25,
T, = 1685.0 ,
R, = 9.84,
H= 15.24,
T, = 1773.0 ,
to ascertain the effects of buoyant convection. Using
!2,=
E = 20.0 ,
-2.31,
K =
10.0,
A-& = 1.57 )
837EO. AK”
7.64 x 10-13 )
and ran to steady state with B taking the values 0, 0.1, 0.2. In each case we used a 20 X 25 grid (Ar = 0.5 1789, A: = 0.63500) and a time step of 0.001. Changing B produced very little change in the distributions of temperature and azimuthal
W.E. Langlois and C.C. Shir, Digitalsimulationof flow patterns in the Czochralskicrystal-pullingprocess
Fig. 2. Steady-state streamlines for the computation
151
in section 4.
velocity. However, the meridional circulation was significantly altered, as can be seen by examining fig. 3. This series of numerical experiments was carried out using an artificially high viscosity (or, equivalently, an unrealistically small crucible) to avoid the necessity of a fine-mesh grid and the concomitant long computation times. Nevertheless, it suggests that buoyancy effects can modify appreciably the impurity-blending estimates obtained from an isothermal study.
References [ 1 ] N. Kobayashi and T. Arizumi, The numerical analyses of the solid-liquid interface shapes during the crystal growth by the Czochralski method, Japan. J. Appl. Phys. 9 (1970) 361-367. [2] N. Kobayashi and T. Arizumi, Computational analysis of the flow in a crucible, J. Crystal Growth 30 (1975) 177-184.
W.E. Langlois and C.C. Shiv. Digital simulation oj’fi’ow patterm in the Czochralski cr_vstal-pullirlgprocess
~
B = 0.
Fig. 3. Steady-state streamlines for the computation section 5.
in
H.J. Scheel, Accelerated crucuble rotation: a novel stirring technique in high-temperature solution growth, J. Crystal Growth 13/14 (1972) 560-565. E.O. Schulz-DuBois, Accelerated crucible rotation: hydrodynamics and stirring effect, J. Crystal Growth 12 (1972) 81-87. W.E. Langlois, Vorticity-streamfunction computation of incompressible fluid flow with an almost-flat free surface, Appl. Math. Modeliing, Vol. 1,4 (1977) 196-198.