116
Journal of Crystal Growth 97(1989) 116—124 North-Holland, Amsterdam
SIMULATION OF STEADY AND TIME-DEPENDENT ROTATION-DRIVEN REGIMES IN A LIQUID-ENCAPSULATED CZOCHRALSKI CONFIGURATION J.-P. FONTAINE, A. RANDRIAMAMPIANINA, G.P. EXTREMET and P. BONTOUX Institut de Mécanique des Fluides, UM 34 du CNRS, I Rue Honnorat, F-13003 Marseille, France
The paper concerns the prediction and the analysis of flow regimes which can develop in liquid-encapsulated Czochralski (LEC) systems under forced convection. Axisymmetric solutions of the Navier—Stokes and energy equations are computed with a finite element method. The influence of the rotation of the crucible and of the crystal is investigated for a fixed temperature gradient. The flows and the temperature distributions are given in the two liquid phases and a number of steady and time-dependent regimes are identified.
1. Introduction In crystal growth from melt, natural convection flows are induced in the liquid phase by the effects of density gradients in the gravity field and/or by the effects of surface-tension gradients. Striations can result from periodic growth and remelts, coming from unsymmetrical heating or from oscillations in the melt. In order to reduce the effects of thermal convection, systems were developed, involving controlled forced-convection mechanisms as the rotations of the crystal and of the crucible in the Czochralski method [1—5].One main goal is to reduce the possible azimuthal variations of the heat transfer to the crystal. Some hydrodynamic aspects of LEC process are considered. The melt is contained in a cylindrical crucible maintained at temperature above the freezing point. The seed crystal is held above the melt surface on the axis of the crucible and is slowly pulled upward. The crystal and the crucible are sometimes rotated in opposite directions in the purpose of isolating the growth interface by providing a viscous shear layer in its vicinity [6—8]. Also for the growth of certain semiconductors as GaAs the melt is covered with a layer of immiscible liquid encapsulant as B203 in order to avoid
are devoted to the steady state and time-dependent axisymmetric regimes in a molten liquid layer (Prin 0.015) confined inside a crucible differentially rotating with the crystal and involving a liquid encapsulant (Pre 10). The rotational Reynolds numbers are varied up to Re ~ 15,000 [3,6—8,12]. =
=
2. Modelling and solution method The simulations concern a model of the actual system shown in fig. 1 assuming the interfaces to be flat. For the axisymmetric motion of a Boussinesq fluid, the governing equations are the coupled Navier—Stokes and energy equations. The thermocapillary flow and the pulling of the crystal are neglected in this study. No-slip conditions are imposed on the wall of the crucible and on the crystal which are maintamed at constant temperatures (hot and cold, respectively). The two immiscible liquids systems are coupled via viscous stress boundary conditions and by the continuity of velocity, temperature and heat flux at the interface. These conditions provide a mechanism for meridional circulations in
the melt. The upper eneapsulant-gas interface is the evaporation of volatil components as arsenic assumed to be insulated and a stress-free condi[9,10]. tion is applied. The numerical simulations are performed using The governing2/v, parameters are number, the Reynolds the Grashof Gr a finite element method (FIDAP [11]). The results number, Re 2R 0022-0248/89/$03.50 © Elsevier Science Publishers B.V. (North-Holland Physics Publishing Division) =
=
J. -P. Fontaine et al.
/ Simulation
ci~ii~
~22
~ V
0 (pulling
T
the heights of the melt and of the encapsulant, and the radii of the crucible and of the crystal. The subscripts m, e, c and s relate respectively to
speed)
R5
crystal. The characteristic scales for lengths, velocthe melt, the encapsulant, the crucible the 1~2 andand R/U, ities and time are R, U=(v/R)Gr respectively. The solution is based on a finite element technique [11]. The aspect ratios have been taken as Am 2, Ae 0.5 and R/RS 2. All the computations were lateral element carried mesh outinvolving with a non-uniform 21 x 21 nodes (confor centrated near the boundaries) four-node quadri-
He
li~~
117
of rotation-driven regimes in LEC configuration
Molten
______
=
MELT
_______________________
=
=
the encapsulant region and 24 x 36 for the melt region (along radial and axial directions, respectively).
R
I Fig. 1. Geometry of the Czochralski LEC model.
3. Hydrodynamicat regimes and effects of rotation ag.~TR3/v2, the Prandtl number, Pr = v/ic, and the aspect ratios Am = Hm/R, Ae = He/R and R/RS, where Hm, He, R and R~ are respectively
In a former paper we have studied the buoyancy driven convection without rotation up to Gr 2 =
-Re:
~
io2 5,
/
xx
I. I
•
/
101
100 100
.
0
101
/
A
___
+ +
0
tO2
io~
to~
Rec
10~
Fig. 2. Regime diagram (Re~,Ret) corresponding to the relative variation of the rotation of the crucible (Ret) and of the seed/crystal (Re~).Characteristic streamline patterns.
118
f-P. Fontaine et at.
/ Simulation
of rotaiion-drwen regimes in LEC configuration
x 10~conditions which are close to realistic LEC applications [12]. The attention is focused, in this study, on the effect of the rotation which is investigated up to actual rotation rates of about 50 rpm (see also ref. [61). However, the temperature gradient is yet taken arbitrarily small as Gr~ 25 and Gr~ 0.25 (the ratio of 100 corresponds to the different viscosities). The crystal and the crucible are rotated in opposite senses (Ret> 0, Re~< 0). The bulk flow is modified following the rota-
tion rates. The various characteristic flow patterns are identified in the regime diagram given in fig. 2. 3.1.
Steady slate regimes
The streamline and isotherm patterns corresponding to the different regimes are displayed in figs. 3a—3c. In these cases the flow is dominated by rotation and the contribution of buoyancy to the meridional flow remains minor. At the lower
=
=
crystal
~ Li
~
H
~ V
___
i
~
I
V
V
c
d
Fig. 3. Streamlines and isotherms for Gr~= 25, Gr~= 0.25. Pr,,, Re~= 37.5 (h); Re,
=
=
0.015, Pr
—900 and Re~= 20(c); Re,
=
3 =10: Re, = —5 and Re~= 20 (a); Re, = —75 and — 125 and Re~= 375 (d).
/ Simulation
J. -P. Fontaine et at.
of rotation -driven regimes in LEC configuration
values of Re~ and Re~(fig. 3a), the characteristic cell is induced by the rotation of the crucible: the flow is ascending along the hot crucible wall and descending along the axis (see refs. [6,12-14]). When I Re., increases with Ret. < Re, an anticlockwise vortex is generated in the neighbourhood of the crystal (fig. 3b). This second TaylorProudman cell is generated by the centrifugal force due to crystal rotation, by sucking the melt up beneath the crystal and by driving it radially
119
outwards near the interface [15]. Computed solutions reported in refs. [7.8,16] show a similar behaviour. We also observe near the meltencapsulant interface that the flow on the encapsulant side is opposed to the flow on the melt side. The primary cell disappears and the upper vortex spreads down into the melt when Re, is further increased up to 900 (figs. 2 and 3c). The streamlines are concentrated close to the encapsulant-melt interface which reveals large mass trans-
933
-.936
-.
936
a —. 940
—.943
—.945 I
16.5
16.7
I
16.7
I
16.6
16.6
I
16.9
~
Fig. 4. Time history for the streamfunction
(a) and streamlines, isotherms and isoazimuthal velocity lines at last time step (h) for Re, = —1675 and Re~= 5000.
i/~
120
/
f. -P. Fontaine et a!.
Simulation of rotation-driven regimes in LEC configuration
port. Large distortions of isotherms are then noticeable in the encapsulant.
region (fig. 4b). For Re~ 10,000 this counterrotating cell is stretched radially with its centre located near the side corner of the crystal (fig. 2). This configuration should give rise to inhomogeneities during the solidification of the crystal and in particular near the corner (the purpose should be to maintain a relative equilibrium in this area in order to keep the radius of the crystal constant during the growth). The isotherms are still weakly affected closely below the crystal for high Reynolds numbers. The vertical temperature gradient appears to remain nearly constant below the crystal, but the radial temperature gradient near the encapsulant layer can vary significantly with the crystal rotation and may then affect the uniformity of the growth interface. Simulations by Crochet et al. [6] (without liquid encapsulant) have shown that a counter-rotating cell develops along the crucible wall for Re~ 375 and Re~ —125, with a steady solution. With the =
3.2. Time-dependent regimes When Re~ is increased according to Re~ 3 Re,, a long transient stage is observed with a time-dependent scheme and weak oscillations establish for Re~ 100 (fig. 2). Fig. 3d presents the iso-q contours and isotherms at the latest time step computed for Re~ 375 and Re, 125. The primary cell (induced by the crucible rotation) fills the entire cavity hut is centered close below the crystal. A secondary cell grows in the encapsulant and dominates progressively, while a third one develops at the corner close to the interface. For Re~> 2,000, a counter-rotating cell appears in the vicinity of the crystal, while an Ekman layer starts to develop in the bottom of the crucible, as shown by deviations of azimuthal velocity pattern in this =
=
=
=
H
I
a
T..1
II
,II
TI(
~
: 1
b
~
L
__
I
.
Fig. 5. Streamlines and isotherms during the tr~in’.ient~iate C,i) •ind ,oom of the upper part of the cavity at last time step (b) for Re, = —5000 and Re~= 15,000.
J.-P. Fontaine et a.
/
1 21
Simulation of rotation-driven regimes in LEC configuration
Vr(z=2) max
I j~vr>O
~. .±~ ~T.T= 5
10_i
10
10
10
tO
10
10
Rec-Res
Fig. 6. Maxima of the radial velocity components along the melt—encapsulant interface (absolute values), of (Re~— Re,).
l,~1I1~(:
2). as a funciion
Vz(r=O) max
z,c,
~>‘
2
to
/‘ 101
,“
zr”
-
...
y
10
1
~O
~‘~1
io2
to~
io5
Rec-Res Fig. 7. Maxima of the axial velocity components along the axis (absolute values) v~~~ax(r = 0) as a function of(Re~— Re,).
122
f. -P. Fontaine et a!.
/ Simulation
of rotation-driven regimes in LEC configuration
Vz(r=O) max i0~
io2
too
10~
10
10
10
10
10
Vr(z=2) max Fig. 8. Maxima of the velocity components (absolute values) t’~”(r = 0) as a function of
present LEC model (using a finer grid close to the walls and the axis than Crochet et al. [61),we did not observe this secondary cell (fig. 3d). The result suggests that the encapsulant modifies the bulk flow in the melt and also the layer under the crystal interface. See the azimuthal velocity pattern and compare to Crochet et al. [6]. Also, for rotation rates as high as Re~~ 2000, the complex flow solution suggests also possible deformations of the shape of the melt—encapsulant interface due to opposite motions in the melt and in the encapsulant (figs. 2 and 4b). The distortions of the isotherms in the liquid encapsulant show that the flow is dominated by the convection, but not by the conduction, as mentioned in refs. [10,17]. Nevertheless, note the small Prandtl number considered, Pr,, 10, while Pre 3270 for B20, [17]. For =
=
Re.~ 5000 and Re, 1675. the time history for the streamfunction near the lower corner of the encapsulant and close to the crystal (fig. 4a) shows weak oscillations of about 0.7%. Oscillations exist in the melt, below the crystal interface and, also, outside the Ekman layer near the bottom. This oscillatory behaviour seems to be related to the =
=
~/i
—.
v,~”(z
=
2).
“small-scale periodicity” phenomenon mentioned by Hjellming and Walker [17]. For larger rotation rates (Ret 15,000, Re, —5000). the flow becomes more intricate: a third cell appears in the encapsulant near the hot wall and it is associated to large distortions of the isotherms; in the melt the stretching process of the cell amplifies radially close to the crystal (fig. 5a). The upper part of the cavity is magnified in fig. Sb. The largest velocities are located under the solid interface where a flattened counter-rotating cell develops. =
=
3.3. Variation of characteristic velocity scales
We have represented in terms of (Rev
—
Re,)
the maximum velocity components (absolute values) v~” along the interface between the melt and the encapsulant (fig. 6), and v~’ along the vertical axis of the cylinder (fig. 7). From the present amount of results, we discuss more precisely four situations: (i) the rotation of the crucible, Ret, when Re, 0; (ii) the rotation of the crystal, Re,, when Re~ 20; (iii) the rotation of the crucible =
=
J. -P. Fontaine et a!.
/ Simulation
of rotation-driven regimes in LEC configuration
when Re, —75; (iv) the relative rotation rates of the crystal and of the crucible as Re,, 3 Re.,. The variations are identified mainly with respect to (Re,, Re,) and (Re,, Re,)~dependencies. The (Re,, Re.,) trends correspond to the established regimes, while the (Re,, Re,)~trends correspond to the transition stages. When increasing (Re,, Re,) we have observed two breaks when the rotation of the crucible becomes dominant with respect to buoyancy for Re,, S and Re., 0. and when the rotation of the crystal dominates the flow for Re,, 20 and Re, 55. For Re, 75, an increase of Re,, results in reducing the speeds. If we consider the different results, the scattering is relatively important for (Re,, Re.,) < 100 when the cellular flow in the melt originates from the rotation of the crucible. For Re,, 3Re.,, the solution is time-dependent and it reaches the (Re,, Re,) behaviour after a wide transition (the mean values of maxima are plotted). The direction of the velocity components depend on the rotation dominated flow: it is downward (negative) for the crucible rotation and upward (positive) for the crystal rotation. The direction of the flow should be important, particularly on the axis below the crystal, because it may determine the shape (convex or concave) of the interface of the crystal. In order to determine the relative importance and the coupling between the vertical flow below the crystal and the radial flow at the melt-encapsulant interface we have represented in fig. 8 the variations of v2’5~ versus v~”1’. This figure suggests the following observations: when the rotation rates are low, the axial velocity, v~, is larger by about one order of magnitude than the radial velocity. ~ and ~ma\ varies linearly with v,~’5’: ~ 5.5 v~’; when the rotation rate is large, the axial and radial velocities are of the same order of magnitude, and also vm’~ 0.75 v~”~ When Re, is increased alone (Re,, 20). the transition between these two laws is clearly shown as corresponding to a variation as v2’~’ (v~~~~)05 with two breaks for small and large (Re,, Re,); when the two rotation rates are increased proportionally (Re,, —3 Re,), the evolution does not =
—
—
123
show any breaks but follows an asymptotical hehaviour.
—
—
4. Discussion
—
—
=
=
=
—
=
=
—
—
—
—
—
—
—
—
—
—
— =
—
A number of flow regimes (steady to time-dcpendent) have been investigated in a liquid-encapsulated Czochralski model when the rotational Reynolds numbers of the crucible and of the crystal vary. The first results already indicate a possible influence of the encapsulant fluid layer on the bulk convection in the melt. The large rotation rates of the crystal are shown to have stabilizing effects and also to reverse the basic flow in the melt (descending on the axis as generated by the buoyancy and the rotation of the crucible) into an ascending flow irrigating the crystal area. The regimes are steady and the investigation was made with respect to former numerical solutions [6,13.14.18]. The large rotation rates of the crucible exhibit time-dependent behaviour and steep shear layers are observed near the meltencapsulant interface and the crystal side. The onset of oscillations may be dependent on the Prandtl number of the encapsulant: in this study we took Pr,,. 10 which is much smaller than the actual values given in the experimental studies [10,17]. The stability of the time-dependent solution was confirmed by increasing the mesh resolution in the crystal zone. Further investigation is planned as the ones previously performed by the authors for a Bridgman type model [19.20]. The existence of axisymmetric solutions has been recently confirmed for situations without rotations [13,14.18]by two- and three-dimensional numerical solutions. The present study has not yet taken into account the influence of thermocapillary effects, neither moving shape interface effects. Recent applications have shown the stabilizing effects of the superimposition of a magnetic field on the oscillation regimes and its role on the good quality of the crystal [17]. Many situations of practical interest have not been yet analysed via numerical sirnulations. The global control of the various mechanisms would suggest further extensive numerical modellings, in particular all the hydrodynamic aspects. =
124
J. -P. Fontaine et a!.
/
Simulation of rotation-driven regimes in LEC configuration
Acknowledgements The authors wish to thank Drs. B. Roux, E. Crespo del Arco and A. Chaouche for their collaboration. They also want to acknowledge Professors F. Rosenberger and R.L. Sani for fruitful discussions and one referee for helpful comments. The computations were carried out on the Cray2 computer with support of the CCVR and using IBM facilities of CNUSC and P.A. St.-Charles. The research was supported by the DRET, the CNES and the CNRS.
References [1] M. Pimputkar and S. Ostrach, J. Crystal Growth 55 (1981) 614. [2] R. Balasubramaniam and S. Ostrach, Physico-Chem. Hydrodyn. 5 (1984) 3. [3] WE. Langlois, Ann. Rev. Fluid Mech. 17 (1985) 191. [4] R.A. Brown, in: Advanced Crystal Growth, Eds. P.M. Dryburgh, B. Cockayne and KG. Barraclough (Prentice Hall. New York, 1987) pp. 41—96. [5] R.A. Brown. J.J. Derby and P.M. Adornato, in: Proc. 1st Intern. Conf.on Electronic Materials Processing, 1986. [6] M.J. Crochet, P.L. Wouters, FT. Geyling and AS. Jordan, J. Crystal Growth 65 (1983) 153. [7] M. Mihel~it,C. Schrock-Pauli, K. Wingerath. H. Wenzl,
W. Uelhoff and A. van der Hart, Notes Numer. Fluid Mech. 5 (1982) 220. [8] [9] [10]
M. Mihelcic, K. 473. Wingerath and C. Pirron, J. Crystal Growth 69 (1984) P.S. Burggraaf, Semiconductor Intern. (June 1982) 44. J.J. Derby, R.A. Brown, FT. Geyling, AS. Jordan and
GA. Nikolakopoulou, J. Electrochem. Soc. 132 (1985) 470. [11] MS. Engelman, Users Manual (Fluid Dynamics Intern., Evanston,FIDAP IL, 1986). [12] J.-P. Fontaine, E. Crespo del Arco, A. Randriamampianina. G.P. Extremet and P. Bontoux, in; Advances in Space Research. Vol. 8, No. 12, Ed. J.C. Legros (Pergamon. Oxford, 1988) pp. 265—279. [13] A. Bottaro and A. Zebib, Phys. Fluids 31(3) (1988) 495. [14] A. Bottaro and A. Zebib, in: Bifurcation Phenomena in Thermal Processes and Convection, Eds. H. Bau, L.A. Bertram and S.A. Korpela (HTD 94-AMD 89. ASME, New York, 1987) p. 113. [15] R. Lamprecht, D. Schwabe. A. Scharmann and E. Schultheiss, J. Crystal Growth 65 (1983) 143. [16] N. Kobayashi, Comp. Methods AppI. Mech. Eng. 23 (1980) 21. [17] L.N. Hjellming and J.S. Walker, Physico-Chem. Hydrodyn. 10 (1988) 107. [18] A. Bottaro and A. Zebib, J. Crystal Growth 97 (1989) 50. [19] J.P. Pulicani, E. Crespo del Arco, A. Randriamampianina, P. Bontoux and R. Peyret, Intern. J. Numerical Methods Fluids, in press. [20] E. Crespo del Arco, A. Randnamampiamna and P. Bontoux, in: Synergetics, Order and Chaos, Ed. MG. Velarde (World Scientific, Singapore, 1988) pp. 244—257.