- Email: [email protected]
Structure of temperature and velocity fields in the Si melt of a Czochralksi crystal growth system
,. . . . . . . .
ELSEVIER
CRYSTAL GROWTH
Journal of Crystal Growth 156 (1995) 383-392
Structure of temperature and velocity fields in the Si melt of a Czochralksi crystal growth system Kyung-Woo Yi a, Vicki B. Booker b, Minoru Eguchi a, Toshiyuki Shyo a, Koichi Kakimoto a,* a Fundamental Research Labs., NEC Corporation, 34 Miyukigaoka, Tsukuba 305, Japan b Department of Chemical Engineering, University of Pennsylvania, 220 S. 33rd Street, Philadelphia, Pennsylvania 19104-6393, USA
Received 17 November 1994;manuscript received in final form 2 June 1995
Abstract
Non-axisymmetric temperature and velocity profiles attributed to crucible rotation in the Si melt of a Czochralski (Cz) system were studied numerically by using three-dimensional (3D) numerical simulation and experimentally by measuring temperature. Two types of non-axisymmetric temperature and velocity profiles were found to exist in the melt. The type observed at low rotation rates occurs under conditions in which Rayleigh-B6nard instability occurs, and the other type of non-axisymmetric flow observed at high rotation rates occurs under conditions corresponding to the occurrence of the baroclinic instability in the rotating annulus. The 3D structure of the temperature field was studied experimentally to confirm the possibility of baroclinic instability in the melt. The predominant trend was a tilt of the temperature field in the same direction as crucible rotation. It was found that the temperature field tilted along the azimuthal direction under conditions in which baroclinic instability occurs.
1. I n t r o d u c t i o n
Many researchers are studying the instabilities in the melt of the Czochralski (Cz) system as a source of the impurity inhomogeneities and point defects in Cz-grown single crystals. Ristorcelli et al, [1] have pointed out that there are many possible causes of instabilities in the melt, and that the instabilities make a three-dimensional (3D) structure (i.e., non-axisymmetric flow and temperature profiles) in the melt even though the
* Corresponding author. Fax: +81 298 566137.
system is completely axisymmetric. The spoke pattern observed on the oxide melt [2-4] is one example of a non-axisymmetric temperature profile. This pattern also appears on the free surface of the silicon melt when the crystal and the crucible are stationary [5]. It is correlated with the roll structures in the velocity field, which are caused by the B6nard-type (Rayleigh-B6nard or Marangoni-B6nard) instability [5]. Kishita et al. [6] reported that when the rotation rate of the crucible increased, the temperature difference along the azimuthal direction showed different behavior: it decreased in the low rotation region while it increased in the high rotation region. The velocity field of the silicon
0022-0248/95/$09.50 © 1995 Elsevier Science B.V. All rights reserved SSDI 0022-0248(95)00277-4
384
K-IV. Y/et al. /Journal of Crystal Growth 156 (1995) 383-392
melt also changed as rotation rate increased [7]. This means that there are different types of nonaxisymmetric profiles in each region. In the case of rotating fluid, different types of instabilities, such as the T a y l o r - C o u e t t e [8], the baroclinic [9,10] and the rotating Rayleigh-B6nard [11] instability, were suggested as the cause of the 3D structure. Because the Cz system contains thermal configurations in which all of these instabilities can occur, it may not be easy to confirm that only one occurs in the Cz melt. The dominant instabilities in the melt depends on the temperature distribution of the crucible and rotation rates of the crystal and the melt. The first purpose of the work reported here is to clarify the origin of the non-axisymmetric pattern in the melt, especially for the case of high crucible rotation. To do this we used a numerical simulation. The results of the numerical simulation showed that two types of 3D structures exist in the Si melt when crucible rotation rate is changed. One occurs when rotation rates are low and disappears when they are high, and the other appears when the rotation rate is high. Analysis of the conditions enhancing or suppressing the instability showed that the profile when rotation rates are low is dominated by the rotating Rayleigh-B6nard instability. On the other hand, the conditions under which the non-axisymmetric profile occurs corresponds to the condition of the occurring baroclinic instability, which was obtained from the rotating annulus experiment. Experimental study is also needed to aid in the understanding of the numerous instabilities which may be generated in the Cz system [1], and the second purpose of the work reported in the present paper was to clarify the 3D structure of molten silicon flow by analyzing the temperature experimentally. The results presented focus on the spatial characteristics of the temperature field for non-axisymmetric case under high crucible rotation rate.
as that reported elsewhere [5]. The diameter of the melt is 70 mm with 35 mm diameter of a crystal, and melt depth is 30 mm. Crucible rotation rates were changed from 0 to 15 rpm, and the wall and the bottom temperatures were changed from Tm + 10 K to Tm + 60 K (Tin: melting temperature). The number of grids were set to 50, 30 and 30 for the azimuthal, radial and axial directions, respectively. The initial conditions for the transient calculation were obtained from the steady state calculation.
3. Experimental procedure
The experimental apparatus was the doublebeam X-ray radiography system described in detail in a previous article [12]. X-ray visualization of the melt allowed visual placement of the thermocouples and allowed their mechanical integrity to be monitored throughout the experiment. Fig. 1 shows a schematic of the system. Temperature was measured using type R; P t / P t - R h thermocouples placed in an AlzO 3 protective tube. The thermocouples were fixed in the inertial coordinate system. The results from four experiments are reported: two of the experiments were designed to detect the horizontal temperature field's wave profile (d12 = 0), and the other two were designed to detect any tilting phenomena (d12 0). A pre-grown single crystal of fixed diameter and length was placed at the melt surface in order to more closely model growth conditions. The crystal diameter was 35 mm, and aspect TC#I
TC#2
2. Calculation method and its results
The basic idea of the numerical simulation and Si properties used in the calculation are the same
Fig. 1. Schematicof experimentalparameters.
K.-W. Yi et al. /Journal of Crystal Growth 156 (1995) 383-392 Table 1 Parameter values of each experiment Experiment: A B C D H, Melt height (mm) 30 30 33 33 r, TCs radial position (mm) 22 22 24 24 Hrc#2, TC#2 position (mm) 10 9 16 18 d12, Vertical spacing (ram) 0 0 9 9 01z, Angle between TCs 12.0° 81.3° 8.3° 10.1°
385
guishing one reading from the other; and any tilt (i.e., incline) can be detected because the thermocouples are located in different vertical planes. Disturbance due to the thermocouples was minimized by locating the second thermocouple T C # 2 lower than the first thermocouple, and rotating the system in the direction from T C # 1 to T C # 2 .
4. Numerical results
ratios of 0.80 and 0.88 were studied. The crystal and crucible were rotated simultaneously in a clockwise direction as viewed from above. The choice of 6 and 8 r p m rotation rates was based on a flow instability m a p of a system with similar p a r a m e t e r s to obtain non-axisymmetric pattern with vortices referred to as rotating non-axismymmetric pattern in the present p a p e r [10]. T e m p e r ature data were collected for 20 min at a sampling rate of 20 Hz. The thermal boundary condition was determined with a constant heater power input: 18.7 kW for 6 r p m and 19.1 k W for 8 rpm. W h e n the rotation rates were changed to begin a new set of measurements, the bottom temperature was allowed to settle to a constant value before data collection began. T h e p a r a m e t e r s of each experiment are listed in Table 1. In experiments A and B, the length of the crystal, after melt-crystal contact was established, was 30 mm. The difference in the vertical position between the thermocouples d12, shown in Fig. 1, was set to zero so that both thermocouples were in the same horizontal plane. The azimuthal wave n u m b e r could then be determined by analyzing the phase difference of the traveling wave. Tracer particles were also placed in the melt and the particle path recorded. Experiments C and D were p e r f o r m e d to investigate possible tilting of the t e m p e r a t u r e field against the pulling axis. After m e l t / c r y s t a l contact was established; the crystal length was 22 ram. In these experiments, the thermocouples were set at fixed radii but different vertical and horizontal cross sections as shown in Fig. 1. If an azimuthal t e m p e r a t u r e wave is traveling with an incline relative to the pulling axis, t e m p e r a t u r e readings from different azimuthal planes will have the same s h a p e with only a phase shift distin-
4.1. Suppression o f the Rayleigh-B~nard instability by rotation W h e n the crucible is not rotating, the temperature profile shows a spoke pattern as reported elsewhere [5]. Ref. [5] also includes a precise description of the spoke pattern in the silicon melt. This pattern is caused by the B6nard-type instability, because cooling by radiation at the free surface of the melt makes a hydrodynamically unstable layer just beneath the free surface, as reported elsewhere [13]. The t e m p e r a t u r e o f the u p p e r boundary of the layer is lower than that of the lower boundary. This layer will be called the "unstable layer" in the present paper. Chandrasekhar [14] reported that the occurrence of the Rayleigh-B6nard instability is suppressed by the rotation of the fluid, so the critical Rayleigh n u m b e r for inducing the Rayleigh-B6nard instability, R a o increases as the rotation rate increases. Axisymmetric t e m p e r a t u r e profiles like these Chandrasekhar described: suppression of the R a y l e i g h - B 6 n a r d instability can be obtained in the present numerical simulation under the condition of a crucible rotation rate of 0.5 rpm. If the rotation rate increases further, a different type of non-axisymmetric t e m p e r a t u r e profile is obtained, as shown in Fig. 2. This profile will be discussed in the next section. Numerical simulations with different boundary conditions and rotation rates showed that when the t e m p e r a t u r e difference between the crystal and the crucible (ATw) increases, a higher crucible rotation speed is necessary to suppress the B6nard-type instability. By using the linear stability theory [14], Chandrasekhar obtained a relation between the rotation rates of the melt and the
K.-W. Yi et aL /Journal of Crystal Growth 156 (1995) 383-392
386
represents
the
Rayleigh number, Ra = g is the acceleration due to gravity, /3 is the thermal expansion coefficient, AT is the t e m p e r a t u r e difference across a fluid layer of height h, and a is the thermal diffusivity. T h e values of h and 2xT used for calculating the R a and 12 of the results of numerical simulations were obtained by taking into account the unstable layer. The values are difficult to determine, however, because h and AT of the unstable layer depend on the radial position [5]. For calculating, the data of the h and AT values o f the unstable layer at the position midway between the edge of the crystal and the inner wall of the crucible were used. The squares in Fig. 3 show that the nonaxisymmetric t e m p e r a t u r e profiles caused by the R a y l e i g h - B 6 n a r d instability a p p e a r in the results of the present numerical simulation under the same condition with 12 and Ra, whereas the circles show that the Rayleigh-B6nard instability does not exist. The line of the R a c for the numerical results should therefore be between the circles and squares. T h e r e are discrepancies between the results of the numerical simulations and the analytic solutions of Chandrasekhar (dotted lines). The R a c obtained from the numerical results is generally greater than that of Chandrasekhar's analytic solution at the same value of 12. In other words, lower rotation rates are sufficient to suppress the Rayleigh-B6nard instability in the Cz melt than those predicted by using the linear instability theory. Therefore, the present numerical simulation cannot reproduce the occurrence of oscillatory convection just above the analytical R a c. The R a c which was reported by Rossby [15] is the same as Chandrasekhar's result, although Rossby also found that there is another transition of fluid flow from oscillatory to irregular convection. The transition conditions are shown by the solid curve in Fig. 3. The numerical results of R a c almost correspond to this curve. This means that the non-axisymmetry of the azimuthal direction obtained in the present study corresponds to the irregular convection of the rotating R a y l e i g h - B 6 n a r d instability observed by Rossby [15]. W h e n the crucible rotation rates are low, the t e m p e r a t u r e and velocity profiles of the
g[3h3AT/ozv, where
(a)
(b)
Fig. 2. (a) Isothermal lines on a plane whose depth from the free surface is one third of the melt height when the crucible rotation rate is 8 rpm. (b) The isothermal lines on the same plane as that of (a) after 8 s.
critical Rayleigh n u m b e r (Ra c) in the rotating fluid within the infinite flat plane. The results of numerical simulations without Marangoni flow were p e r f o r m e d with different boundary conditions and rotation rates are c o m p a r e d in Fig. 3 with the analytic results. The results of the present numerical simulation are plotted along with the results of Chandrasekhar (dotted line) and Rossby's results obtained in an experiment with mercury in a cylinder (solid line) [15]. The horizontal axis of Fig. 3 is the dimensionless rotating rate, 12 = 12ch2/u, where 12c is the rotation rate of the crucible, h is the height of the fluid layer, and v is the kinematic viscosity. The vertical axis
10 6
. . . . . . . .
i
. . . . . . . .
i
. . . . . . . .
i
. . . . . . .
. . . . . . . .
i
. . . . . .
J
I-- 10 4
2 1010°
....
i
. . . . . . . .
i
102 f) (rch2/v )
104
Fig. 3. Relationship between the critical Rayleigh number and the crucible rotation rates: (solid line) the transition between oscillatory and irregular convection, determined experimentally by Rossby; (dotted line) the analytic results of Chandrasekhar; (squares) the Rayleigh-B6nard instability shown in numerical results; (circles) the Rayleigh-B6nard instability is not shown in numerical results.
K.-W. Yi et aL /Journal of Crystal Growth 156 (1995) 383-392
whole melt are dependent on the behavior of the non-axisymmetric profile originated by the Rayleigh instability. Therefore, it seems that the rotating Rayleigh-B6nard instability might play an important role in determining the temperature and velocity profiles in the melt when the crucible rotation rate is low. 4.2. Another asymmetry appearing at high rotation rates
At high crucible rotation rates, the temperature profile becomes non-axisymmetric again. Fig. 2a shows one example of this non-axisymmetric temperature profile on a plane one third of the melt height below the free surface when the crystal and crucible are both rotating at 8 rpm and in the same direction. The shape of the isothermal lines of a rotating non-axisymmetric pattern differs from that of the non-axisymmetric patterns at low rotation rates. Moreover, when the rotation rate was increased with the constant temperature boundary conditions, the rotating non-axisymmetric pattern started to appear not only after the temperature profile became axisymmetric but also when the patterns of the Rayleigh-B6nard instability remained. Therefore, the cause of the rotating non-axisymmetric pattern appears to differ from that of the Rayleigh-B6nard instability. The Taylor-Couette and the baroclinic instabilities are possible causes of the rotating non-axisymmetric pattern. The baroclinic instability can be considered the main cause of the rotating non-axisymmetric pattern in the two possible instabilities, since the silicon melt of the Czochralski crystal-growing system has temperature boundary conditions similar to those of the experimental equipment (the rotating annulus) for the baroclinic instability. The rotating non-axisymmetric pattern continuously changed its Shape even though boundary conditions and rotating rates were constant during the calculations. The temperature profile on the same plane was changed from that shown in Fig. 2a to that shown in Fig. 2b during 8 s in real time, for example. Fig. 4 shows the calculated temperature change half the radius and half the depth of the melt in the inertial coordinate sys-
387
143(
E
E
142
I
I
r
20
I
I
40
60
Time (sec) Fig. 4. Change in calculated temperature over a period of one minute at half the radius and half the height of melt depth.
tem. We can identify that temperature oscillates with several frequencies. The temperature oscillation is caused by the rotation of the non-axisymmetric pattern. The change in the amplitude of the temperature oscillation, on the other hand, is caused by the time-dependent deformation of the rotating non-axisymmetric temperature profile. Fig. 4 also shows that the temperature change is periodical. Therefore, the temperature fluctuation originates from the rotation of the non-
5 CM/SEC
Fig. 5. Velocity profiles on the same plane as that in Fig. 2a observed from the rotating view point with the same crucible rotation rate.
K.IW. Yi et aL /Journal of Crystal Growth 156 (1995) 383-392
388
axisymmetric profiles and the change in the rotating non-axisymmetric pattern. This non-axisymmetric temperature field is strongly coupled with the velocity field of the melt. Fig. 5 shows the velocity profiles on the same plane as that of Fig. 5a. The most interesting feature in this profile is the formation of two vortices shown in the periphery of the melt. The structure with two vortices in the silicon melt of the Cz crystal growing system was observed in an X-ray radiography experiment with a 3-inch crucible [7]. The vortices formation in the rotating melt was usually observed in the rotating annulus experiments. Hide and co-workers [16] observed velocity profiles similar to the present numerical results in a rotating-annulus experiment with a very small inner crucible radius. The above results show that the non-axisymmetric velocity and temperature profiles obtained from the 3D numerical simulation on the silicon melt in the Cz system are similar to those caused by baroclinic instability in rotating-annulus experiments.
4.3. Formation of the rotating non-axisymmetric profiles The occurrence of the baroclinic instability is well predicted using two dimensionless numbers: the thermal Rossby number, Ro t (=
I
.......
........
~
........
[]
d" 102 r'c
I
........
J
........
I
........
I
......
-,, b a r o c l i n i c f l o w ~,~ n o n - b a r o c l i n i c fic
•
LINE ]
•
<9 10 o BAROCLINIC
ION O
10-2 .......
i
........
106
[
........
i
........
i
........
I
........
108 1010 T a ( 4 ~ 2R5/u 2d)
i
,,,
i'012
Fig. 6. T h e diagram of the baroclinic instability with calculated results. Squares are the results of normal Cz conditions and triangles are values calculated assuming a thermally adiabatic free surface. T h e two lines are obtained from the experim e n t s of Fein and Pfeffer [20].
g/3ATh/f22R 2) and the Taylor number Ta ( = 4g-22RS/u2h) [9]. Here, R, g, AT, /3, and u are radius, acceleration of gravity, temperature difference between crucible and crystal, volume expansion coefficient, and kinematic viscosity, respectively. Fig. 6 shows a map of flow instability including the calculated results. The lines I and II show the upper and lower boundaries for the occurrence of baroclinic instability in rotating-annulus experiments [17]. The region between the lines in which the non-axisymmetric profile caused by the baroclinic instability was observed, All results of experiments on different liquids that covered a wide range of Prandtl numbers in various sizes of rotating annulus can be expressed by the lines in Fig. 6. The results of the numerical simulation are also plotted in Fig. 6. The filled squares show that rotating instability exists and the open square means that the rotating instability does not exist. If the baroclinic instability really exists in the melt, its occurrence should also show agreement with the boundaries expressed by line I. However, the numerical results for the Si melt in a crucible of 7 cm diameter, show that the transition to baroclinic instability occurs at a Ro t number much higher than that predicted from line I, as shown in Fig. 6. The experiment in the silicon melt of the Cz system also reported such a high value of the Ro t [18]. This m e a n s that the lowest crucible rotation rate at which the rotating non-axisymmetric profile is induced is lower than the value expected from the studies on the baroclinic instability. To explain this discrepancy, it is necessary that the definition of the Ro t be examined carefully. The difference in the critical Ro t number between the value of the rotating annulus with water and the value of the C z crucible with Si is mainly caused by the difference of the vertical temperature profile in the liquid. The horizontal temperature difference in the liquid, ATh, was usually used as the value of the AT term for calculating the R o t number for the diagram in Fig. 6. However, the vertical temperature difference ATv must originally be used as the value of the term AT in the Ro t number [19]. Under the assumption that the ATh and ATv are similar
K.-W. Yi et al. /Journal of Crystal Growth 156 (1995) 383-392
[18], the ATh was used instead of the ATv. This assumption might be acceptable in the rotatingannulus experiments performed at room temperature because the free surface of the fluid could be assumed to be thermally adiabatic. On the other hand, ATv in the Si melt is much smaller than ATh because the free surface of the Si melt is cooled by radiative heat transfer. Therefore, the actual Ro t number in the melt should be smaller than the calculated value. Because calculation of the actual Ro t of the melt is very difficult, the validity of this reasoning was confirmed by carrying out other calculations assuming that the free surface is thermally insulated. The results calculated under the assumption of adiabatic conditions are plotted by triangles in Fig. 6. The filled triangles show conditions under which rotating non-axisymmetric profiles exist, whereas the open triangles show conditions under which the results are axisymmetric. The occurrence of the baroclinic instability in these results shows good agreement with the occurrence of baroclinic instability in the upper boundary (line I), thus indicating that the radiation cooling at the free surface promotes the occurrence of the baroclinic instability by decreasing the temperature gradients along the axial direction. The most characteristic behavior of the baroclinic wave is tilting of the temperature field in the azimuthal direction [19]. Fig. 7 shows calculated temperature profile in the 0 - z plane when the vortex structure was formed in which the tilting of the temperature field can be recognized.
5. Experimental results
5.1. Data analysis As the numerical simulation predicted, tilting of the temperature profile along the azimuthal direction is possible to occur in the melt at a higher crucible rotation rate. Driving forces at the boundary of the fluid include the heating at the crucible walls inducing a buoyant flow, and the crucible and crystal co-rotation. T h e resulting fluid motion carries the temperature pattern between T C # 1 and T C # 2 . The azimuthal wavenumber N represents the ratio of the temperature wave's local frequency f to the phase velocity c(f) of the wave traveling between TC#1 and TC#2:
N= 2rrf/c(f).
L.o Fig. 7. T e m p e r a t u r e p r o f i l e in t h e 0 - z p l a n e .
(1)
This expression requires the frequency spectrum of each temperature signal. By direct Fourier transformation of the original time histories x(t) and y(t) obtained from TC#1 and TC#2, respectively, we obtained the expression of the signal in the frequency domain. A total of 20 480 sampled data points were analyzed using a record length T of 102.4 s (2048 data points). Hence, 10 data records were averaged, and the resulting minimum frequency-resolution bandwidth was 0.00976 Hz and the Nyquist cutoff frequency was 10 Hz. The 012 is the distance (in radians) between the thermocouples as shown in Fig. 1. Using the value of the phase difference /3(f) and the angle 012 between the thermocouples, we could calculate the degree of tilt. An angle 0t relative to a horizontal line through the measuring point at thermocouple # 2 Can be expressed as Ot--=tan-1[d12/r(/3'-012)],
z
389
(2)
where fl' is observed phase shift. The /3 is equal to fi'/N, where N is the wave number. Any tilt of the. temperature field relative to a vertical line is represented by the angle q~T, defined as 45T = 90° - 0t.
(3)
390
1C-W. Yi et al./Journal of Crystal Growth 156 (1995) 383-392
,i
'
'
,
,
t
. . . .
>.I--
1
,
15000
'
i
\
........ TC#1 -
-
TC#2
ii
+t H
&-
I
0.075
. . . .
i
i
0.100
t
i
i
0.75
LIJ
o O -1m 0.50 ~ m z 0-<
oF ~E I- = "n>~ O_ 1~-~ <
Expt. C
/
Expt. A
).25
gL
a
. . . .
r 6Sl
!
w w 0 n
. . . .
1,00
z w cn d < n-
13_
i
b '1
,
0.125
i
i
i
I
i
0,150
,
100
0
JL. . . .
0.00
,.!,~J
I
0.10
0,20
0.30
0.40
(Hz) F R E Q U E N C Y (Hz) Fig. 8. (a) A plot of the power spectral density function (left ordinate) of TC#1 and TC#2, and the corresponding coherency function (right ordinate) for experiment C. The crucible and crystal rotation rate was 6 rpm. (b) Amplitude results of experiments of A and C. The crystal and crucible co-rotation was 8 rpm. FREQUENCY
The tilting angle is defined such that + q~T represents a tilt of the temperature field in the same direction as crystal and crucible rotation. Temperature oscillations can be observed for crucible and crystal co-rotation at 6 rpm by using two thermocouples of T C # 1 and T C # 2 in Fig. 1. The result is from experiment C. The signal from T C # 2 lags behind a time lag of the signal from T C # 1 . This means the wave was traveling in the same direction as the crucible and crystal, since both rotated in the direction from T C # 1 to T C # 2 . The wave traveling direction was also determined by flow visualization. Fig. 8a shows the power spectral density of each signal. Both signals had a major frequency component of maximum amplitude at 0.10938 Hz. This was approximately 10% higher than the frequency of the periodic driving force applied by the crystal and crucible at the boundaries. No other peaks were detected outside of the range shown on the abscissa. The plot of the coherence function in the same figure shows good agreement between the two signals, since the value is equal to 1.0 at most frequencies. From each of the four experiments, comparison of the power spectral density and the resulting coherence function yields results similar to those shown in Fig. 8a. Hence, in each experi-
ment the thermocouples measured a signal of the same shape, even when the thermocouples were set at different vertical positions. Only differences in phase were observed. This implies that the temperature fields detected in different horizontal planes were the same. Fig. 8b shows the amplitude results for the cross-spectrum function for experiments A and C when the crucible and crystal rotate at 8 rpm. Experiment A was used to detect the azimuthal wave while experiment C was used to detect any temperature tilting against the pulling axis. The dashed line in the figure corresponds to the system rotation rate. It is most interesting to note the general behavior observed at both rotation rates. Although melt height and the vertical placement of the thermocouples differed between experiments in Table 1, the peak
Table 2 Calculated phase velocityand wave number Experiment A ExperimentB 6rpm c(f) (rpm)
2.98 Tracer average azimuthal velocity (rpm) N 2.33
8rpm
6rpm
8rpm
3.62 -
3.14 3.20
4.32 -
2.48
2.07
2.05
K.-W. Y~et al. /Journal of Crystal Growth 156 (1995) 383-392
was an average 8.9% higher than the system rotational frequency for the 6 rpm case, and 10.6% higher for the 8 rpm case. From experiments A and B, the azimuthal wave number of the maximum amplitude peak was found. The wave number was calculated from the phase velocity, by using Eq. (1). Table 2 lists results of the phase velocity and wave number calculations. The origin of the result of c ( f ) << w c is not clarified till now. However, azimuthal wave velocity should not always be larger than crucible rotation rate, since the azimuthal wave velocity is determined by momentum transfer at the m e l t crucible, melt-crystal and melt-gas interfaces. Therefore, the velocity of the traveling wave with smaller or larger angular velocity can be observed in the condition of baroclinic instability. In one experiment, the average azimuthal velocity of a tracer particle was calculated by tracking the particle's path in three dimensions. The value obtained, 3.2 rpm, was in good agreement with the 3.14 rpm calculated from F F T analysis. The degree of tilting in experiments C and D was calculated using Eqs. (2) and (3). Since the wave numbers found in the experiment A and B averaged 2.40 and 2.06, respectively, these two values were used to estimate error range of azimuthal temperature wave tilt. In three of the cases, the temperature profile was found to tilt in the azimuthal direction with 9° to 14°. In one case, however, tilting in the direction opposite to the direction of crucible and crystal rotation was observed. 5.2. Discussion o f experimental data
Although the wave number calculated in experiment B was very nearly an integer, the wave number calculated was not in experiment A. The non-integer value can be explained, roughly by taking experimental error into account. All the waves are thought to have wave number with integers, since a wave with a non-integer value cannot keep its stable oscillation as shown in Fig. 8a, Another characteristic of waves observed in the heated rotating-annulus geometry is a temperature field tilt against the pulling axis in the
391
direction of rotation. Pfeffer et al. [20] demonstrated that, in their annulus, waves in the temperature field tilt slightly along the azimuthal direction of rotation. This is consistent with the theoretical analysis of Eady [21]. While most of the cases we observed followed this pattern, we observed one case of slight tilting in t h e opposite direction. This seemingly errant case can provide insight to differences between the instability mechanisms of the rotating annulus geometry and the Czochralski system. Further work using a numerical model to investigate the observed tilting behavior and to make quantitative comparisons between the results of experiments and simulations will help elucidate the observed tilting phenomena.
6. Summary of the instabilities in the silicon melt of a Cz system
We showed that there are two types of instabilities in a rotating Si melt. At low rotation rates, the Rayleigh-B6nard instability appears and forms a spoke-pattern-like temperature profile in the melt. The instability is suppressed and finally disappeared by the rotation of the crucible. This behavior implies that the Rayleigh-B6nard instability is dominating for determining the temperature and velocity profiles at low rotation rates. When crucible rotation speed is increased, the rotating non-axisymmetric patterns appear in the melt. The profiles of these patterns are quite similar to those of the non-axisymmetric structure made by the baroclinic instability. Moreover, the conditions under which patterns occur show good agreement with those under which baroclinic instability occurs. Therefore, it seems that the baroclinic instability occurs in the silicon melt of a Cz system. Fluid flow in the melt is converted from axisymmetric to n o n - a x i s y m m e t r i c by the Rayleigh-B6nard instability at low crucible rotation rates and by the baroclinic instability at high crucible rotation rates. A tilt of the azimuthal temperature wave against the pulling axis was found at both 6 and 8 rpm crystal and crucible co-rotation rates. A sim-
392
K-W. Yi et al. /Journal of Crystal Growth 156 (1995) 383-392
pie frequency spectrum existed for each case, and 2-lobed azimuthal temperature waves seemed most probable.
Acknowledgements Part of this research was performed under a joint research project funded by the U.S. National Science Foundation Grant INT-9317251 "Dissertation Enhancement Award" granted to Dr. Stuart W. Churchill, University of Pennsylvania. We thank Dr. Roy Lang and Dr. Masashi Mizuta of the Fundamental Research Laboratories, NEC Corporation, for their support and continuous encouragement; and we greatly appreciate the fruitful discussions with Dr. Taketoshi Hibiya and Dr. Shin Nakamura. One of us (V.B.B.) thanks Dr. Stuart Churchill for his continuous encouragement and guidance.
References [1] J.R. Ristorcelli and J.L. Lumsby, J. Crystal Growth 116 (1992) 447.
[2] J.C. Brice, T.M. Bruton, O.F. Hill and P.A.C. Whiffin, J. Crystal Growth 24/25 (1974) 429. [3] A.D.W. Jones, J. Crystal Growth 63 (1983) 70. [4] D.C. Miller and T.L. Pernell, J. Crystal Growth 57 (1982) 253. [5] K.-W. Yi, K. Kakimoto, M. Eguchi, M. Watanabe, T. Shyo and T. Hibiya, J. Crystal Growth 144 (1994) 20, [6] Y. Kishita, M. Tanaka and H. Esaka, J. Crystal Growth 130 (1993) 75. [7] K. Kakimoto, M. Watanabe, M. Eguchi and T. Hibiya, J. Crystal Growth 126 (1993) 435. [8] J.C. Brice and P.A.C. Whiffin, J. Crystal Growth 38 (1977) 245. [9] C.D. Brandle, J. Crystal Growth 57 (1982) 65. [10] K. Kakimoto, M. Eguchi, H. Watanabe and T. Hibiya, J. Crystal Growth 102 (1990) 16. [11] A. Seidl, G. McCord, G. Miiller and H.-J. Leister, J. Crystal Growth 137 (1994) 326. [12] M. Watanabe, M. Eguchi, K. Kakimoto, and T. Hibiya, J. Crystal Growth 133 (1993) 28. [13] A.D.W. Jones, Phys. Fluids 28 (1985) 31. [14] S. Chandrasekhar, Hydrodynamics and Hydromagnetics Stability (Oxford, London, 1961) pp. 76-135. [15] H.T. Rossby, J. Fluid Mech. 36 (1969) 309. [16] J.S. Fein and R.L. Pfeffer, J. Fluid Mech. 75 (1976) 81. [17] R. Hide and P.J. Mason, Phil. Trans. R. Soc. A 268 (1970) 201. [18] M. Watanabe, M. Eguchi, K. Kakimoto, Y. Baros and T. Hibiya, J. Crystal Growth 128 (1993) 288. [19] R. Hide, Phys. Fluid 10 (1967) 56. [20] R.L. Pfeffer, G. Buzyna and R. Kung, J. Atm. Sci. 37 (1980) 2129. [21] E.T. Eady, Tellus 1 (1949) 33.
ELSEVIER
CRYSTAL GROWTH
Journal of Crystal Growth 156 (1995) 383-392
Structure of temperature and velocity fields in the Si melt of a Czochralksi crystal growth system Kyung-Woo Yi a, Vicki B. Booker b, Minoru Eguchi a, Toshiyuki Shyo a, Koichi Kakimoto a,* a Fundamental Research Labs., NEC Corporation, 34 Miyukigaoka, Tsukuba 305, Japan b Department of Chemical Engineering, University of Pennsylvania, 220 S. 33rd Street, Philadelphia, Pennsylvania 19104-6393, USA
Received 17 November 1994;manuscript received in final form 2 June 1995
Abstract
Non-axisymmetric temperature and velocity profiles attributed to crucible rotation in the Si melt of a Czochralski (Cz) system were studied numerically by using three-dimensional (3D) numerical simulation and experimentally by measuring temperature. Two types of non-axisymmetric temperature and velocity profiles were found to exist in the melt. The type observed at low rotation rates occurs under conditions in which Rayleigh-B6nard instability occurs, and the other type of non-axisymmetric flow observed at high rotation rates occurs under conditions corresponding to the occurrence of the baroclinic instability in the rotating annulus. The 3D structure of the temperature field was studied experimentally to confirm the possibility of baroclinic instability in the melt. The predominant trend was a tilt of the temperature field in the same direction as crucible rotation. It was found that the temperature field tilted along the azimuthal direction under conditions in which baroclinic instability occurs.
1. I n t r o d u c t i o n
Many researchers are studying the instabilities in the melt of the Czochralski (Cz) system as a source of the impurity inhomogeneities and point defects in Cz-grown single crystals. Ristorcelli et al, [1] have pointed out that there are many possible causes of instabilities in the melt, and that the instabilities make a three-dimensional (3D) structure (i.e., non-axisymmetric flow and temperature profiles) in the melt even though the
* Corresponding author. Fax: +81 298 566137.
system is completely axisymmetric. The spoke pattern observed on the oxide melt [2-4] is one example of a non-axisymmetric temperature profile. This pattern also appears on the free surface of the silicon melt when the crystal and the crucible are stationary [5]. It is correlated with the roll structures in the velocity field, which are caused by the B6nard-type (Rayleigh-B6nard or Marangoni-B6nard) instability [5]. Kishita et al. [6] reported that when the rotation rate of the crucible increased, the temperature difference along the azimuthal direction showed different behavior: it decreased in the low rotation region while it increased in the high rotation region. The velocity field of the silicon
0022-0248/95/$09.50 © 1995 Elsevier Science B.V. All rights reserved SSDI 0022-0248(95)00277-4
384
K-IV. Y/et al. /Journal of Crystal Growth 156 (1995) 383-392
melt also changed as rotation rate increased [7]. This means that there are different types of nonaxisymmetric profiles in each region. In the case of rotating fluid, different types of instabilities, such as the T a y l o r - C o u e t t e [8], the baroclinic [9,10] and the rotating Rayleigh-B6nard [11] instability, were suggested as the cause of the 3D structure. Because the Cz system contains thermal configurations in which all of these instabilities can occur, it may not be easy to confirm that only one occurs in the Cz melt. The dominant instabilities in the melt depends on the temperature distribution of the crucible and rotation rates of the crystal and the melt. The first purpose of the work reported here is to clarify the origin of the non-axisymmetric pattern in the melt, especially for the case of high crucible rotation. To do this we used a numerical simulation. The results of the numerical simulation showed that two types of 3D structures exist in the Si melt when crucible rotation rate is changed. One occurs when rotation rates are low and disappears when they are high, and the other appears when the rotation rate is high. Analysis of the conditions enhancing or suppressing the instability showed that the profile when rotation rates are low is dominated by the rotating Rayleigh-B6nard instability. On the other hand, the conditions under which the non-axisymmetric profile occurs corresponds to the condition of the occurring baroclinic instability, which was obtained from the rotating annulus experiment. Experimental study is also needed to aid in the understanding of the numerous instabilities which may be generated in the Cz system [1], and the second purpose of the work reported in the present paper was to clarify the 3D structure of molten silicon flow by analyzing the temperature experimentally. The results presented focus on the spatial characteristics of the temperature field for non-axisymmetric case under high crucible rotation rate.
as that reported elsewhere [5]. The diameter of the melt is 70 mm with 35 mm diameter of a crystal, and melt depth is 30 mm. Crucible rotation rates were changed from 0 to 15 rpm, and the wall and the bottom temperatures were changed from Tm + 10 K to Tm + 60 K (Tin: melting temperature). The number of grids were set to 50, 30 and 30 for the azimuthal, radial and axial directions, respectively. The initial conditions for the transient calculation were obtained from the steady state calculation.
3. Experimental procedure
The experimental apparatus was the doublebeam X-ray radiography system described in detail in a previous article [12]. X-ray visualization of the melt allowed visual placement of the thermocouples and allowed their mechanical integrity to be monitored throughout the experiment. Fig. 1 shows a schematic of the system. Temperature was measured using type R; P t / P t - R h thermocouples placed in an AlzO 3 protective tube. The thermocouples were fixed in the inertial coordinate system. The results from four experiments are reported: two of the experiments were designed to detect the horizontal temperature field's wave profile (d12 = 0), and the other two were designed to detect any tilting phenomena (d12 0). A pre-grown single crystal of fixed diameter and length was placed at the melt surface in order to more closely model growth conditions. The crystal diameter was 35 mm, and aspect TC#I
TC#2
2. Calculation method and its results
The basic idea of the numerical simulation and Si properties used in the calculation are the same
Fig. 1. Schematicof experimentalparameters.
K.-W. Yi et al. /Journal of Crystal Growth 156 (1995) 383-392 Table 1 Parameter values of each experiment Experiment: A B C D H, Melt height (mm) 30 30 33 33 r, TCs radial position (mm) 22 22 24 24 Hrc#2, TC#2 position (mm) 10 9 16 18 d12, Vertical spacing (ram) 0 0 9 9 01z, Angle between TCs 12.0° 81.3° 8.3° 10.1°
385
guishing one reading from the other; and any tilt (i.e., incline) can be detected because the thermocouples are located in different vertical planes. Disturbance due to the thermocouples was minimized by locating the second thermocouple T C # 2 lower than the first thermocouple, and rotating the system in the direction from T C # 1 to T C # 2 .
4. Numerical results
ratios of 0.80 and 0.88 were studied. The crystal and crucible were rotated simultaneously in a clockwise direction as viewed from above. The choice of 6 and 8 r p m rotation rates was based on a flow instability m a p of a system with similar p a r a m e t e r s to obtain non-axisymmetric pattern with vortices referred to as rotating non-axismymmetric pattern in the present p a p e r [10]. T e m p e r ature data were collected for 20 min at a sampling rate of 20 Hz. The thermal boundary condition was determined with a constant heater power input: 18.7 kW for 6 r p m and 19.1 k W for 8 rpm. W h e n the rotation rates were changed to begin a new set of measurements, the bottom temperature was allowed to settle to a constant value before data collection began. T h e p a r a m e t e r s of each experiment are listed in Table 1. In experiments A and B, the length of the crystal, after melt-crystal contact was established, was 30 mm. The difference in the vertical position between the thermocouples d12, shown in Fig. 1, was set to zero so that both thermocouples were in the same horizontal plane. The azimuthal wave n u m b e r could then be determined by analyzing the phase difference of the traveling wave. Tracer particles were also placed in the melt and the particle path recorded. Experiments C and D were p e r f o r m e d to investigate possible tilting of the t e m p e r a t u r e field against the pulling axis. After m e l t / c r y s t a l contact was established; the crystal length was 22 ram. In these experiments, the thermocouples were set at fixed radii but different vertical and horizontal cross sections as shown in Fig. 1. If an azimuthal t e m p e r a t u r e wave is traveling with an incline relative to the pulling axis, t e m p e r a t u r e readings from different azimuthal planes will have the same s h a p e with only a phase shift distin-
4.1. Suppression o f the Rayleigh-B~nard instability by rotation W h e n the crucible is not rotating, the temperature profile shows a spoke pattern as reported elsewhere [5]. Ref. [5] also includes a precise description of the spoke pattern in the silicon melt. This pattern is caused by the B6nard-type instability, because cooling by radiation at the free surface of the melt makes a hydrodynamically unstable layer just beneath the free surface, as reported elsewhere [13]. The t e m p e r a t u r e o f the u p p e r boundary of the layer is lower than that of the lower boundary. This layer will be called the "unstable layer" in the present paper. Chandrasekhar [14] reported that the occurrence of the Rayleigh-B6nard instability is suppressed by the rotation of the fluid, so the critical Rayleigh n u m b e r for inducing the Rayleigh-B6nard instability, R a o increases as the rotation rate increases. Axisymmetric t e m p e r a t u r e profiles like these Chandrasekhar described: suppression of the R a y l e i g h - B 6 n a r d instability can be obtained in the present numerical simulation under the condition of a crucible rotation rate of 0.5 rpm. If the rotation rate increases further, a different type of non-axisymmetric t e m p e r a t u r e profile is obtained, as shown in Fig. 2. This profile will be discussed in the next section. Numerical simulations with different boundary conditions and rotation rates showed that when the t e m p e r a t u r e difference between the crystal and the crucible (ATw) increases, a higher crucible rotation speed is necessary to suppress the B6nard-type instability. By using the linear stability theory [14], Chandrasekhar obtained a relation between the rotation rates of the melt and the
K.-W. Yi et aL /Journal of Crystal Growth 156 (1995) 383-392
386
represents
the
Rayleigh number, Ra = g is the acceleration due to gravity, /3 is the thermal expansion coefficient, AT is the t e m p e r a t u r e difference across a fluid layer of height h, and a is the thermal diffusivity. T h e values of h and 2xT used for calculating the R a and 12 of the results of numerical simulations were obtained by taking into account the unstable layer. The values are difficult to determine, however, because h and AT of the unstable layer depend on the radial position [5]. For calculating, the data of the h and AT values o f the unstable layer at the position midway between the edge of the crystal and the inner wall of the crucible were used. The squares in Fig. 3 show that the nonaxisymmetric t e m p e r a t u r e profiles caused by the R a y l e i g h - B 6 n a r d instability a p p e a r in the results of the present numerical simulation under the same condition with 12 and Ra, whereas the circles show that the Rayleigh-B6nard instability does not exist. The line of the R a c for the numerical results should therefore be between the circles and squares. T h e r e are discrepancies between the results of the numerical simulations and the analytic solutions of Chandrasekhar (dotted lines). The R a c obtained from the numerical results is generally greater than that of Chandrasekhar's analytic solution at the same value of 12. In other words, lower rotation rates are sufficient to suppress the Rayleigh-B6nard instability in the Cz melt than those predicted by using the linear instability theory. Therefore, the present numerical simulation cannot reproduce the occurrence of oscillatory convection just above the analytical R a c. The R a c which was reported by Rossby [15] is the same as Chandrasekhar's result, although Rossby also found that there is another transition of fluid flow from oscillatory to irregular convection. The transition conditions are shown by the solid curve in Fig. 3. The numerical results of R a c almost correspond to this curve. This means that the non-axisymmetry of the azimuthal direction obtained in the present study corresponds to the irregular convection of the rotating R a y l e i g h - B 6 n a r d instability observed by Rossby [15]. W h e n the crucible rotation rates are low, the t e m p e r a t u r e and velocity profiles of the
g[3h3AT/ozv, where
(a)
(b)
Fig. 2. (a) Isothermal lines on a plane whose depth from the free surface is one third of the melt height when the crucible rotation rate is 8 rpm. (b) The isothermal lines on the same plane as that of (a) after 8 s.
critical Rayleigh n u m b e r (Ra c) in the rotating fluid within the infinite flat plane. The results of numerical simulations without Marangoni flow were p e r f o r m e d with different boundary conditions and rotation rates are c o m p a r e d in Fig. 3 with the analytic results. The results of the present numerical simulation are plotted along with the results of Chandrasekhar (dotted line) and Rossby's results obtained in an experiment with mercury in a cylinder (solid line) [15]. The horizontal axis of Fig. 3 is the dimensionless rotating rate, 12 = 12ch2/u, where 12c is the rotation rate of the crucible, h is the height of the fluid layer, and v is the kinematic viscosity. The vertical axis
10 6
. . . . . . . .
i
. . . . . . . .
i
. . . . . . . .
i
. . . . . . .
. . . . . . . .
i
. . . . . .
J
I-- 10 4
2 1010°
....
i
. . . . . . . .
i
102 f) (rch2/v )
104
Fig. 3. Relationship between the critical Rayleigh number and the crucible rotation rates: (solid line) the transition between oscillatory and irregular convection, determined experimentally by Rossby; (dotted line) the analytic results of Chandrasekhar; (squares) the Rayleigh-B6nard instability shown in numerical results; (circles) the Rayleigh-B6nard instability is not shown in numerical results.
K.-W. Yi et aL /Journal of Crystal Growth 156 (1995) 383-392
whole melt are dependent on the behavior of the non-axisymmetric profile originated by the Rayleigh instability. Therefore, it seems that the rotating Rayleigh-B6nard instability might play an important role in determining the temperature and velocity profiles in the melt when the crucible rotation rate is low. 4.2. Another asymmetry appearing at high rotation rates
At high crucible rotation rates, the temperature profile becomes non-axisymmetric again. Fig. 2a shows one example of this non-axisymmetric temperature profile on a plane one third of the melt height below the free surface when the crystal and crucible are both rotating at 8 rpm and in the same direction. The shape of the isothermal lines of a rotating non-axisymmetric pattern differs from that of the non-axisymmetric patterns at low rotation rates. Moreover, when the rotation rate was increased with the constant temperature boundary conditions, the rotating non-axisymmetric pattern started to appear not only after the temperature profile became axisymmetric but also when the patterns of the Rayleigh-B6nard instability remained. Therefore, the cause of the rotating non-axisymmetric pattern appears to differ from that of the Rayleigh-B6nard instability. The Taylor-Couette and the baroclinic instabilities are possible causes of the rotating non-axisymmetric pattern. The baroclinic instability can be considered the main cause of the rotating non-axisymmetric pattern in the two possible instabilities, since the silicon melt of the Czochralski crystal-growing system has temperature boundary conditions similar to those of the experimental equipment (the rotating annulus) for the baroclinic instability. The rotating non-axisymmetric pattern continuously changed its Shape even though boundary conditions and rotating rates were constant during the calculations. The temperature profile on the same plane was changed from that shown in Fig. 2a to that shown in Fig. 2b during 8 s in real time, for example. Fig. 4 shows the calculated temperature change half the radius and half the depth of the melt in the inertial coordinate sys-
387
143(
E
E
142
I
I
r
20
I
I
40
60
Time (sec) Fig. 4. Change in calculated temperature over a period of one minute at half the radius and half the height of melt depth.
tem. We can identify that temperature oscillates with several frequencies. The temperature oscillation is caused by the rotation of the non-axisymmetric pattern. The change in the amplitude of the temperature oscillation, on the other hand, is caused by the time-dependent deformation of the rotating non-axisymmetric temperature profile. Fig. 4 also shows that the temperature change is periodical. Therefore, the temperature fluctuation originates from the rotation of the non-
5 CM/SEC
Fig. 5. Velocity profiles on the same plane as that in Fig. 2a observed from the rotating view point with the same crucible rotation rate.
K.IW. Yi et aL /Journal of Crystal Growth 156 (1995) 383-392
388
axisymmetric profiles and the change in the rotating non-axisymmetric pattern. This non-axisymmetric temperature field is strongly coupled with the velocity field of the melt. Fig. 5 shows the velocity profiles on the same plane as that of Fig. 5a. The most interesting feature in this profile is the formation of two vortices shown in the periphery of the melt. The structure with two vortices in the silicon melt of the Cz crystal growing system was observed in an X-ray radiography experiment with a 3-inch crucible [7]. The vortices formation in the rotating melt was usually observed in the rotating annulus experiments. Hide and co-workers [16] observed velocity profiles similar to the present numerical results in a rotating-annulus experiment with a very small inner crucible radius. The above results show that the non-axisymmetric velocity and temperature profiles obtained from the 3D numerical simulation on the silicon melt in the Cz system are similar to those caused by baroclinic instability in rotating-annulus experiments.
4.3. Formation of the rotating non-axisymmetric profiles The occurrence of the baroclinic instability is well predicted using two dimensionless numbers: the thermal Rossby number, Ro t (=
I
.......
........
~
........
[]
d" 102 r'c
I
........
J
........
I
........
I
......
-,, b a r o c l i n i c f l o w ~,~ n o n - b a r o c l i n i c fic
•
LINE ]
•
<9 10 o BAROCLINIC
ION O
10-2 .......
i
........
106
[
........
i
........
i
........
I
........
108 1010 T a ( 4 ~ 2R5/u 2d)
i
,,,
i'012
Fig. 6. T h e diagram of the baroclinic instability with calculated results. Squares are the results of normal Cz conditions and triangles are values calculated assuming a thermally adiabatic free surface. T h e two lines are obtained from the experim e n t s of Fein and Pfeffer [20].
g/3ATh/f22R 2) and the Taylor number Ta ( = 4g-22RS/u2h) [9]. Here, R, g, AT, /3, and u are radius, acceleration of gravity, temperature difference between crucible and crystal, volume expansion coefficient, and kinematic viscosity, respectively. Fig. 6 shows a map of flow instability including the calculated results. The lines I and II show the upper and lower boundaries for the occurrence of baroclinic instability in rotating-annulus experiments [17]. The region between the lines in which the non-axisymmetric profile caused by the baroclinic instability was observed, All results of experiments on different liquids that covered a wide range of Prandtl numbers in various sizes of rotating annulus can be expressed by the lines in Fig. 6. The results of the numerical simulation are also plotted in Fig. 6. The filled squares show that rotating instability exists and the open square means that the rotating instability does not exist. If the baroclinic instability really exists in the melt, its occurrence should also show agreement with the boundaries expressed by line I. However, the numerical results for the Si melt in a crucible of 7 cm diameter, show that the transition to baroclinic instability occurs at a Ro t number much higher than that predicted from line I, as shown in Fig. 6. The experiment in the silicon melt of the Cz system also reported such a high value of the Ro t [18]. This m e a n s that the lowest crucible rotation rate at which the rotating non-axisymmetric profile is induced is lower than the value expected from the studies on the baroclinic instability. To explain this discrepancy, it is necessary that the definition of the Ro t be examined carefully. The difference in the critical Ro t number between the value of the rotating annulus with water and the value of the C z crucible with Si is mainly caused by the difference of the vertical temperature profile in the liquid. The horizontal temperature difference in the liquid, ATh, was usually used as the value of the AT term for calculating the R o t number for the diagram in Fig. 6. However, the vertical temperature difference ATv must originally be used as the value of the term AT in the Ro t number [19]. Under the assumption that the ATh and ATv are similar
K.-W. Yi et al. /Journal of Crystal Growth 156 (1995) 383-392
[18], the ATh was used instead of the ATv. This assumption might be acceptable in the rotatingannulus experiments performed at room temperature because the free surface of the fluid could be assumed to be thermally adiabatic. On the other hand, ATv in the Si melt is much smaller than ATh because the free surface of the Si melt is cooled by radiative heat transfer. Therefore, the actual Ro t number in the melt should be smaller than the calculated value. Because calculation of the actual Ro t of the melt is very difficult, the validity of this reasoning was confirmed by carrying out other calculations assuming that the free surface is thermally insulated. The results calculated under the assumption of adiabatic conditions are plotted by triangles in Fig. 6. The filled triangles show conditions under which rotating non-axisymmetric profiles exist, whereas the open triangles show conditions under which the results are axisymmetric. The occurrence of the baroclinic instability in these results shows good agreement with the occurrence of baroclinic instability in the upper boundary (line I), thus indicating that the radiation cooling at the free surface promotes the occurrence of the baroclinic instability by decreasing the temperature gradients along the axial direction. The most characteristic behavior of the baroclinic wave is tilting of the temperature field in the azimuthal direction [19]. Fig. 7 shows calculated temperature profile in the 0 - z plane when the vortex structure was formed in which the tilting of the temperature field can be recognized.
5. Experimental results
5.1. Data analysis As the numerical simulation predicted, tilting of the temperature profile along the azimuthal direction is possible to occur in the melt at a higher crucible rotation rate. Driving forces at the boundary of the fluid include the heating at the crucible walls inducing a buoyant flow, and the crucible and crystal co-rotation. T h e resulting fluid motion carries the temperature pattern between T C # 1 and T C # 2 . The azimuthal wavenumber N represents the ratio of the temperature wave's local frequency f to the phase velocity c(f) of the wave traveling between TC#1 and TC#2:
N= 2rrf/c(f).
L.o Fig. 7. T e m p e r a t u r e p r o f i l e in t h e 0 - z p l a n e .
(1)
This expression requires the frequency spectrum of each temperature signal. By direct Fourier transformation of the original time histories x(t) and y(t) obtained from TC#1 and TC#2, respectively, we obtained the expression of the signal in the frequency domain. A total of 20 480 sampled data points were analyzed using a record length T of 102.4 s (2048 data points). Hence, 10 data records were averaged, and the resulting minimum frequency-resolution bandwidth was 0.00976 Hz and the Nyquist cutoff frequency was 10 Hz. The 012 is the distance (in radians) between the thermocouples as shown in Fig. 1. Using the value of the phase difference /3(f) and the angle 012 between the thermocouples, we could calculate the degree of tilt. An angle 0t relative to a horizontal line through the measuring point at thermocouple # 2 Can be expressed as Ot--=tan-1[d12/r(/3'-012)],
z
389
(2)
where fl' is observed phase shift. The /3 is equal to fi'/N, where N is the wave number. Any tilt of the. temperature field relative to a vertical line is represented by the angle q~T, defined as 45T = 90° - 0t.
(3)
390
1C-W. Yi et al./Journal of Crystal Growth 156 (1995) 383-392
,i
'
'
,
,
t
. . . .
>.I--
1
,
15000
'
i
\
........ TC#1 -
-
TC#2
ii
+t H
&-
I
0.075
. . . .
i
i
0.100
t
i
i
0.75
LIJ
o O -1m 0.50 ~ m z 0-<
oF ~E I- = "n>~ O_ 1~-~ <
Expt. C
/
Expt. A
).25
gL
a
. . . .
r 6Sl
!
w w 0 n
. . . .
1,00
z w cn d < n-
13_
i
b '1
,
0.125
i
i
i
I
i
0,150
,
100
0
JL. . . .
0.00
,.!,~J
I
0.10
0,20
0.30
0.40
(Hz) F R E Q U E N C Y (Hz) Fig. 8. (a) A plot of the power spectral density function (left ordinate) of TC#1 and TC#2, and the corresponding coherency function (right ordinate) for experiment C. The crucible and crystal rotation rate was 6 rpm. (b) Amplitude results of experiments of A and C. The crystal and crucible co-rotation was 8 rpm. FREQUENCY
The tilting angle is defined such that + q~T represents a tilt of the temperature field in the same direction as crystal and crucible rotation. Temperature oscillations can be observed for crucible and crystal co-rotation at 6 rpm by using two thermocouples of T C # 1 and T C # 2 in Fig. 1. The result is from experiment C. The signal from T C # 2 lags behind a time lag of the signal from T C # 1 . This means the wave was traveling in the same direction as the crucible and crystal, since both rotated in the direction from T C # 1 to T C # 2 . The wave traveling direction was also determined by flow visualization. Fig. 8a shows the power spectral density of each signal. Both signals had a major frequency component of maximum amplitude at 0.10938 Hz. This was approximately 10% higher than the frequency of the periodic driving force applied by the crystal and crucible at the boundaries. No other peaks were detected outside of the range shown on the abscissa. The plot of the coherence function in the same figure shows good agreement between the two signals, since the value is equal to 1.0 at most frequencies. From each of the four experiments, comparison of the power spectral density and the resulting coherence function yields results similar to those shown in Fig. 8a. Hence, in each experi-
ment the thermocouples measured a signal of the same shape, even when the thermocouples were set at different vertical positions. Only differences in phase were observed. This implies that the temperature fields detected in different horizontal planes were the same. Fig. 8b shows the amplitude results for the cross-spectrum function for experiments A and C when the crucible and crystal rotate at 8 rpm. Experiment A was used to detect the azimuthal wave while experiment C was used to detect any temperature tilting against the pulling axis. The dashed line in the figure corresponds to the system rotation rate. It is most interesting to note the general behavior observed at both rotation rates. Although melt height and the vertical placement of the thermocouples differed between experiments in Table 1, the peak
Table 2 Calculated phase velocityand wave number Experiment A ExperimentB 6rpm c(f) (rpm)
2.98 Tracer average azimuthal velocity (rpm) N 2.33
8rpm
6rpm
8rpm
3.62 -
3.14 3.20
4.32 -
2.48
2.07
2.05
K.-W. Y~et al. /Journal of Crystal Growth 156 (1995) 383-392
was an average 8.9% higher than the system rotational frequency for the 6 rpm case, and 10.6% higher for the 8 rpm case. From experiments A and B, the azimuthal wave number of the maximum amplitude peak was found. The wave number was calculated from the phase velocity, by using Eq. (1). Table 2 lists results of the phase velocity and wave number calculations. The origin of the result of c ( f ) << w c is not clarified till now. However, azimuthal wave velocity should not always be larger than crucible rotation rate, since the azimuthal wave velocity is determined by momentum transfer at the m e l t crucible, melt-crystal and melt-gas interfaces. Therefore, the velocity of the traveling wave with smaller or larger angular velocity can be observed in the condition of baroclinic instability. In one experiment, the average azimuthal velocity of a tracer particle was calculated by tracking the particle's path in three dimensions. The value obtained, 3.2 rpm, was in good agreement with the 3.14 rpm calculated from F F T analysis. The degree of tilting in experiments C and D was calculated using Eqs. (2) and (3). Since the wave numbers found in the experiment A and B averaged 2.40 and 2.06, respectively, these two values were used to estimate error range of azimuthal temperature wave tilt. In three of the cases, the temperature profile was found to tilt in the azimuthal direction with 9° to 14°. In one case, however, tilting in the direction opposite to the direction of crucible and crystal rotation was observed. 5.2. Discussion o f experimental data
Although the wave number calculated in experiment B was very nearly an integer, the wave number calculated was not in experiment A. The non-integer value can be explained, roughly by taking experimental error into account. All the waves are thought to have wave number with integers, since a wave with a non-integer value cannot keep its stable oscillation as shown in Fig. 8a, Another characteristic of waves observed in the heated rotating-annulus geometry is a temperature field tilt against the pulling axis in the
391
direction of rotation. Pfeffer et al. [20] demonstrated that, in their annulus, waves in the temperature field tilt slightly along the azimuthal direction of rotation. This is consistent with the theoretical analysis of Eady [21]. While most of the cases we observed followed this pattern, we observed one case of slight tilting in t h e opposite direction. This seemingly errant case can provide insight to differences between the instability mechanisms of the rotating annulus geometry and the Czochralski system. Further work using a numerical model to investigate the observed tilting behavior and to make quantitative comparisons between the results of experiments and simulations will help elucidate the observed tilting phenomena.
6. Summary of the instabilities in the silicon melt of a Cz system
We showed that there are two types of instabilities in a rotating Si melt. At low rotation rates, the Rayleigh-B6nard instability appears and forms a spoke-pattern-like temperature profile in the melt. The instability is suppressed and finally disappeared by the rotation of the crucible. This behavior implies that the Rayleigh-B6nard instability is dominating for determining the temperature and velocity profiles at low rotation rates. When crucible rotation speed is increased, the rotating non-axisymmetric patterns appear in the melt. The profiles of these patterns are quite similar to those of the non-axisymmetric structure made by the baroclinic instability. Moreover, the conditions under which patterns occur show good agreement with those under which baroclinic instability occurs. Therefore, it seems that the baroclinic instability occurs in the silicon melt of a Cz system. Fluid flow in the melt is converted from axisymmetric to n o n - a x i s y m m e t r i c by the Rayleigh-B6nard instability at low crucible rotation rates and by the baroclinic instability at high crucible rotation rates. A tilt of the azimuthal temperature wave against the pulling axis was found at both 6 and 8 rpm crystal and crucible co-rotation rates. A sim-
392
K-W. Yi et al. /Journal of Crystal Growth 156 (1995) 383-392
pie frequency spectrum existed for each case, and 2-lobed azimuthal temperature waves seemed most probable.
Acknowledgements Part of this research was performed under a joint research project funded by the U.S. National Science Foundation Grant INT-9317251 "Dissertation Enhancement Award" granted to Dr. Stuart W. Churchill, University of Pennsylvania. We thank Dr. Roy Lang and Dr. Masashi Mizuta of the Fundamental Research Laboratories, NEC Corporation, for their support and continuous encouragement; and we greatly appreciate the fruitful discussions with Dr. Taketoshi Hibiya and Dr. Shin Nakamura. One of us (V.B.B.) thanks Dr. Stuart Churchill for his continuous encouragement and guidance.
References [1] J.R. Ristorcelli and J.L. Lumsby, J. Crystal Growth 116 (1992) 447.
[2] J.C. Brice, T.M. Bruton, O.F. Hill and P.A.C. Whiffin, J. Crystal Growth 24/25 (1974) 429. [3] A.D.W. Jones, J. Crystal Growth 63 (1983) 70. [4] D.C. Miller and T.L. Pernell, J. Crystal Growth 57 (1982) 253. [5] K.-W. Yi, K. Kakimoto, M. Eguchi, M. Watanabe, T. Shyo and T. Hibiya, J. Crystal Growth 144 (1994) 20, [6] Y. Kishita, M. Tanaka and H. Esaka, J. Crystal Growth 130 (1993) 75. [7] K. Kakimoto, M. Watanabe, M. Eguchi and T. Hibiya, J. Crystal Growth 126 (1993) 435. [8] J.C. Brice and P.A.C. Whiffin, J. Crystal Growth 38 (1977) 245. [9] C.D. Brandle, J. Crystal Growth 57 (1982) 65. [10] K. Kakimoto, M. Eguchi, H. Watanabe and T. Hibiya, J. Crystal Growth 102 (1990) 16. [11] A. Seidl, G. McCord, G. Miiller and H.-J. Leister, J. Crystal Growth 137 (1994) 326. [12] M. Watanabe, M. Eguchi, K. Kakimoto, and T. Hibiya, J. Crystal Growth 133 (1993) 28. [13] A.D.W. Jones, Phys. Fluids 28 (1985) 31. [14] S. Chandrasekhar, Hydrodynamics and Hydromagnetics Stability (Oxford, London, 1961) pp. 76-135. [15] H.T. Rossby, J. Fluid Mech. 36 (1969) 309. [16] J.S. Fein and R.L. Pfeffer, J. Fluid Mech. 75 (1976) 81. [17] R. Hide and P.J. Mason, Phil. Trans. R. Soc. A 268 (1970) 201. [18] M. Watanabe, M. Eguchi, K. Kakimoto, Y. Baros and T. Hibiya, J. Crystal Growth 128 (1993) 288. [19] R. Hide, Phys. Fluid 10 (1967) 56. [20] R.L. Pfeffer, G. Buzyna and R. Kung, J. Atm. Sci. 37 (1980) 2129. [21] E.T. Eady, Tellus 1 (1949) 33.