Hydrodynamics and heat exchange of crystal pulling from melts. Part I: Experimental studies of free convection mode

This work is a brief overview of experimental study results for hydrodynamics and convective heat exchange in thermal gravity capillary convection modes for the classic Czochralski technique setup obtained at the Institute of Thermophysics, Siberian Branch of the Russian Academy of Sciences. The experiments have been carried out at test benches which simulated the physics of the Czochralski technique for 80 and 295 mm diameter crucibles. Melt simulating fluids with Prandtl numbers Pr = 0.05, 16, 45.6 and 2700 have been used. Experiments with transparent fluids have been used for comparing the evolution of flow structure from laminar mode to well-developed turbulent mode. Advanced visualization and measurement methods have been used. The regularities of local and integral convective heat exchange in the crucible/melt/crystal system have been studied. The experiments have shown that there are threshold Grashof and Marangoni numbers at which the structure of the thermal gravity capillary flow undergoes qualitative changes and hence the regularities of heat exchange in the melt change. The effect of melt hydrodynamics on the crystallization front shape has been studied for Pr = 45.6. Crystallization front shapes have been determined for the 1 × 105 to 1.9 × 105 range of Grashof numbers. We show that the crystallization front shape depends largely on the spatial flow pattern and the temperature distribution in the melt.

Below we will present the first part of this four-part overview of experimental and numerical results for convection in various fluids simulating various setups of the Czochralski technique obtained at the Institute of Thermophysics, Siberian Branch of the Russian Academy of Sciences. By now the most complete hydrodynamic and convective heat exchange simulation data have been obtained for the classic fixed crucible setup of the Czochralski technique. The aim of the present series of studies is to simulate the crystal growth conditions at different process stages and develop a fundamental basis for improving the methods of controlling the melt hydrodynamics and the heat exchange conditions. When designing the benches for the simulation of physical models of various directional crystallization techniques, we used proprietary methods for the visualization and parameter measurement of velocity and temperature fields developed at the IT SB RAS for the fundamental investigations of near-wall turbulence [42]. The experiments dealt with laminar-turbulent transition (LTT) processes in free-convention boundary layers and heat exchange at free liquid surfaces [42] which can be considered to be models of the processes occurring at the melt/ambience boundary. Similar boundary layer LTT and heat exchange studies were also carried out for regions at vertical walls [31,42] which can be considered as models of gas dynamic processes at the side surfaces of crystals.
Under real process conditions, melt hydrodynamics are determined by the cooperative action of a set of bulk forces (buoyancy, centrifugal, Coriolis, shear and electromagnetic forces) and surface forces (thermo-and concentration-capillary, Laplacian and friction forces). If the effect of the buoyancy forces is significant, the temperature and velocity fields are self-consistent. The presence of the thermocapillary effect further increases the hydrodynamics/heat exchange feedback. The development of hydrodynamic control methods requires understanding of the specific influence mechanisms (or the specific expression regularities) and the relative contributions of each of the forces in question to the formation of the melt flow structure and analysis of the results of their nonlinear interaction. This problem cannot be solved using solely experimental tools since the buoyancy forces which are inevitable in Earth conditions and other bulk forces act jointly with surface forces. Controlling their relative contributions or separately studying their individual action under the conditions of a physical experiment is an extremely complex or even impossible task. Therefore there is a need for combined complementary experimental and numerical investigations. Results of numerical investigations into free convection flow modes will be presented in the second part of this overview, and experimental and numerical results for the effect of the abovementioned forces on the convective processes in the melt during crystal and crucible rotation will be discussed in the third and the fourth parts.

Experimental setup
Natural convection was experimentally studied at different stages from laminar mode to well-developed turbulent flow mode [12][13][14]30]. The experiments were carried out at benches which simulated the physics of the Czhochralski technique for 80 to 660 mm diameter crucibles using the visualization and measurement methods described in detail earlier [12]. Digital videos were recorded for illustrating the flow of visualized transparent melt-simulating fluids at different magnifications. Computer processing of the videos provides qualitative and quantitative information on the shape and dimensions of the main and secondary flows and quantitative information on the velocity fields.
Below we present the results of studies conducted using two physical models of the classic Czhochralski technique. These models had almost similar designs of the working zones, but the crucibles differed considerably in size: the crucible diameter was D c = 80 mm for the first section and D c = 295 mm for the second section. The working media were gallium/indium/tin melt (Prandtl number Pr = 0.05), ethyl alcohol (96 %, Pr = 16), PES-5 siloxane fluid (Pr = 2700), and saturated hydrocarbons, i.e., hexadecane and heptadecane (Pr = 45.6). The two latter fluids have the solidification points at 18 and 22 °C, respectively, and allow studying crystallization processes for a transparent fluid. Transparent fluids were visualized using specially selected fraction of scale-shaped aluminum particles 10 to 15 µm in size. Fluid movement was observed in the central section (the r-z plane). The observation area was continuously illuminated with 1 to 2 mm thick planar light beams. The observation results were also photographed and video recorded in the r-z plane at specific preset levels relative to the bottom. In these cases the working zones had transparent bottoms and were observed through a rotating mirror installed underneath at a 45 arc deg angle to the bottom plane. The thermocouples for temperature measurement at the crucible surface and at the model crystals were made from 80 µm diam. copper and 60 µm diam. constantan wires. The temperature distribution in the working fluid volume was measured with the thermocouple probes. The measurement increment was controlled accurate to 0.01 mm along both coordinates. The nichrome-constantan probe thermocouples (0.1 mm wire diameter) were electrochemically etched off to 10-25 µm diam. at the working solders. The solders of the three thermocouple probes were arranged one under another in the same vertical plane so the temperature distribution could be measured in the CF boundary layers, at the crucible side wall and at the free melt surface. The two top thermocouple probes were installed with a constant spacing (δz ≤ 0.3 mm); this arrangement allowed measuring local (instantaneous and average) heat flows in the fluid for differential switching of the thermocouple probes and in simultaneous signal recording mode. The thermocouple sensitivity was 42 µV/K. The accuracy of the measurements with Sch300 or F30 microvoltmeter/ ammeters was 1 µV.

Results and discussion
Thermal gravity convection caused by the temperature difference between the CF and the crucible walls in the Czochralski technique is fundamentally unavoidable and can hardly be controlled, and this is the case for any non-isothermal system in the gravity field. Its physical nature, contribution and intensity are determined by the temperature dependence of density and by the effect of buoyancy (Archimedes) forces and are described by the Grashof number Gr = (βg/ν 2 )ΔTR sc 3 or the Rayleigh number Ra = GrPr = (βg/aν)ΔTR sc 3 , where g is the gravity acceleration, β is the volumetric expansion coefficient, ν = µ/ρ is the kinematic viscosity coefficient, µ is the dynamic viscosity coefficient, ρ is the density, a = λ/ρC p is the thermal diffusivity coefficient, λ is the thermal conductivity coefficient, C p is the constant pressure heat capacity and R sc is the crystal radius. The Prandtl number Pr = ν/a describes the thermophysical properties of the melt. The presence of a free melt surface portion and the temperature drop ΔT along the free surface between the CF and the crucible walls leads to a combined effect of the buoyancy forces and the thermocapillary effect. The thermocapillary effect originates from the temperature dependence of the surface tension coefficient σ and the forces acting along the melt surface in radial directions. As a rule the surface tension coefficient decreases with an increase in temperature. The resultant tangential force directed from the heated melt surface portion to the cold one causes thermocapillary convection. The intensity of thermocapillary convection is described by the Marangoni number Ma = (-δσ/δT)ΔTR sc /aµ. Thermal gravity capillary convection (TGCC) in the melt is always the initial mode for the search of optimum crystal growth conditions, or the working mode for the growth of some oxide crystals. Depending on the free surface state and the conditions of heat emission to the ambience, this initial mode may have variable relative contributions of the buoyancy forces and the thermocapillary effect. The two limit cases are thermal gravity convection mode for which the surface tension coefficient depends but slightly on temperature and gravity capillary convection mode with a free adiabatic boundary for which the contribution of the thermocapillary effect is the greatest in the presence of a radial temperature gradient.
For free convection which always has a gravity capillary origin in a physical experiment, the flow arrangement is as follows. In the center (Fig. 1a-c) a cooled melt downward stream forms which is directed from the bottom boundary of the crystal model toward the crucible bottom. Then the melt spreads over the adiabatic bottom towards the hot crucible walls, followed by a upward flow along the crucible wall and a flow toward the crystal edge along the free surface. A large superficial swirl forms in the top melt portion (Fig. 1a, b) in which the fluid flow direction coincides with that of the global meridian flow. This swirl has a predominantly thermocapillary origin. The reverse flow occurs under the free surface due to an excess buoyancy of the melt. Figure 1a shows a strictly axially symmetrical steadystate laminar flow. In laminar mode the relative height of the melt layer (in the 0.4 ≤ H/R c ≤ 1.5 range) and the ratio of the model crystal radii (in the 7.0 ≤ R c /R sc ≤ 1.3 range) have a negligible effect on the qualitative aspect of the spatial pattern of the flow. This mode occurs in media with medium Pr numbers.
In highly viscous media with Pr = 2700 the thermocapillary effect has but a little effect on the spatial pattern of the laminar flow, and no superficial swirl forms. In non-transparent gallium/indium/tin melt with Pr = 0.05 the temperature distributions normal to the CF were measured with microprobe thermocouples. The temperature profiles were similar to those observed in media with Pr = 16 ÷ 2700 having clearly delimited thermal boundary layer portions and a radial distribution of the local heat fluxes that is typical for free laminar convection.
In TGCC mode an increase in the temperature gradient and hence a growth of Gr and Ma cause an unstable and non-steady-state flow which transits gradually to random turbulent mode. The qualitative parameters of the transition to turbulent flow mode depend largely on the relative and absolute dimensions of the thermo-hydrodynamic system in question. For a 80 mm crucible diameter the volume of fluid in turbulent mode is filled with variable size swirls (Fig. 1b). The pattern of the fluid flows shown here corresponds to an instantaneous state of the system. For a 295 mm crucible diameter and Gr having the same or comparable order of magnitude, the flow in the boundary layers at the contour of the working area is turbulent, but there is also a large core of stagnant fluid (Fig. 1c).
The Gr values are almost the same. Laminar fluid convection mode with Pr = 16 shown in Fig. 1a can be implemented for models with up to 100 mm crucible diameter at temperature gradients of within 2 K. Some details of the 3D non-steady-state flow transition scenarios are universal and depend neither on the size (in a certain size range) nor on Pr. During ΔT growth for crystal models with R c /R sc ≥ 3 radii ratios, the cold fluid downward stream loses stability at an initial stage: low-frequency axial-symmetric single-harmonic velocity and temperature oscillations develop; then the second harmonic appears after the next threshold ΔТ (or Gr) value is reached, and at the next stage the oscillations become random. The second harmonic has a lower frequency. The analog signal of the thermocouple located in the sub-crystal melt region acquires a harmonic oscillation form with low-frequency amplitude modulation. Instability of this type also develops with an increase in ΔT in fluids with Pr = 2700 for any crystal radius. Interpretation of the video recording for the sub-crystal melt region shows that this instability shows itself as a periodical increase in the thickness of the boundary layer of the fluid cooled at the CF, combined with a periodical entrainment of the outer part of the boundary layer by the downward stream after the threshold thickness is reached. What is observed resembles the formation of a water droplet on the tip of an icicle with its periodical tear-off. The transition to random velocity and temperature oscillations is caused by the transition of the melt in the crucible to a stable stratified state (i.e. the melt becomes divided in layers by temperature and hence density) and the formation of a bottom melt layer the temperature of which is lower than the bulk-average one. As a result the cooled melt stream impinging onto a relatively cold layer starts oscillating randomly in radial and azimuth directions. The stream eventually becomes turbulent and the oscillations penetrate into the sub-crystal melt region.
Rayleigh-Benard instability develops under the CF in large-diameter crystals: cells forming in the boundary layer are entrained by the main flow from the crystal edge toward the center. Figure 2 shows general results of the studies for TGCC mode in fluid with Pr = 16 including visual observations, measurements and analytical parameters of the temperature field. At H/R c = 0.7 but for different crystal model radii, one can trace the following temperature oscillation regularities during an increase in the temperature difference. In Region I to the left of Curve 1 the global convective flow is laminar at low Gr. The boundary layers are laminar under the entire model CF, and temperature oscillations are absent. Region II between Curves 1 and 2 corresponds to laminar-turbulent transition mode in the free convective boundary layer at the CF. Region III to the right of Curve 2 corresponds to free turbulent convection mode. Each of these regions has qualitatively different temperature oscillation spectrum patterns (Figs. [3][4][5]. Region IIa between Curves 1 and 3 is distinguished by low-frequency low-amplitude irregular temperature oscillations. The oscillation spectrum patterns are typical of noise signals (see Fig. 3).
Once threshold conditions are reached (Region IIb between Curves 3 and 4) when space-ordered secondary flows emerge, the temperature oscillations under the CF at a certain distance from the front edge are regular, the respective spectra becoming discrete. Oscillation power spectra typical of this region are shown in Fig. 4. Discrete spectra of this type correspond to the paired swirls which are observed under the cold surface of the crystal model and entrained by the main flow directed from the edge to the downward stream. These secondary swirls have a Rayleigh-Benard origin and emerge once the threshold conditions are reached, i.e. the threshold thickness of the fluid layer cooled down from above and the threshold temperature difference in the boundary layer. The flows in the swirls of such pairs have opposite directions. A stationary thermocouple probe detects quasi-periodical signals corresponding to the alternating up and down flows in the entrained swirls. The vertical and horizontal dimensions of the swirls are comparable with the boundary layer thickness. The temperature oscillation amplitude at the outer edge of the hot boundary layer is approximately ΔT/2.
An increase in Gr causes a shrinkage of the ordered secondary flow existence region and noise contamination of the temperature oscillation spectra. The outer edge of the boundary layer exhibits irregular outbursts of the melt cooled under the CF to far beyond the boundary layer. In Region IIc (Fig. 2) the ordered secondary flow zone shrinks to r ≤ R sc /2, the spectra becoming increasingly noise-contaminated as one moves down the flow, i.e., their quality depends on r. Curve 2 corresponds to the transition to irregular and relatively high-frequency oscillations under the entire CF (Fig. 2). The oscillation spectra for the turbulent boundary layer become continuous in the entire frequency range typical of free convection, the highest frequencies being about 1 Hz.
The temperature oscillation behavior features noted above correlate with the regularities of heat boundary layer development at the CF and the local heat emission regularities. Figure 6 shows radial distributions of time-averaged local heat flows. The curves in Fig. 6 a correspond to different modes (the case of R c = 80 mm, H/R c = 0.7 and R c /R sc = 2.76) for an increase in the temperature difference from 0.41 to 10.3 °C and Gr from 0.62 × 10 4 to 3.13 × 10 5 : (1) stable laminar boundary layer; (2) irregular temperature oscillations in Region IIa of Fig. 2; (3) emergence of regular secondary flow in the boundary layer; (4-6) modes with Gr above which the boundary layer is completely turbulent. Curves 3 and 4 in Fig. 6a clearly  show laminar boundary layer sections between the front edge r/R sc = 1 and the local minimum in the q(r) curve. The minimum position does not coincide with the secondary flow emergence point. The rise in the q(r) curve down the flow is caused by the development of secondary flow and an intensification of the heat exchange between the local flows normal to the CF. The drop in the q(r) curve after the local maximum corresponds to the transition to turbulent boundary layer. With an increase in Gr this local maximum shifts toward the front edge. Curves 5 and 6 correspond to modes in which the flows impinging onto the crystal edge are non-steady-state. Curves 2-6 were obtained by measuring time-averaged temperature profiles normal to the CF, i.e., time-averaged temperature gradients. The radial heat flux distributions at the CF with increasing Gr were obtained in a wide range of Gr on physical models with crucible diameters of 80 and 290 mm with relative radii for crystal models R c /R sc ranging from 1.29 to 7.0 [12].   For a large-scale model with a 295 mm diameter crucible the laminar-turbulent transition processes in the boundary layer at the CF occur at as small ΔT as several tenths of a Kelvin. Typical q(r) curves for this case are shown in Fig. 6b in order of increasing temperature difference from 0.78 °C (curves 1-4) to 2.38 °C (curve 5). These curves clearly distinct the classic laminar boundary layer regions: the region between the front edge x/R c = 0 and the q(r) curve minimum, the transition region between the minimum and the maximum at x/R c ≈ 0.4 and the turbulent boundary layer region at 0.4 ≤ x/R c ≤ 1. Figure 6c also shows the RMS deviations of the local heat flow from the time-averaged value. These deviations are uniquely related to the amplitude of local temperature oscillations. It can be seen from Fig. 6a that the temperature oscillations at large-diameter CFs emerge before the heat-induced laminar-turbulent transition boundary. The secondary flows achieve saturation before the transition to the turbulent boundary layer. The temperature oscillation amplitudes in the region of well-developed quasi-2D secondary flows and in the region of 3D irregular secondary flow are similar, i.e., they do not change in the direction down the flow. The results presented in Fig. 6b, c were obtained at H/R c = 0.33 and R c / R sc = 3.69. Results of quite detailed investigations into the evolution of temperature oscillation spectra in radial directions with an increase in ΔT were presented in an earlier work [12] which dealt with the ΔT dependence of the main harmonic frequency of the discrete spectra. The growth of the main harmonic frequency is mainly caused by an increase in the rate of ordered secondary swirl entrainment by the main meridian flow relative to the stationary probe.
Practically important feature of free convective flow modes is that the abovementioned CF boundary layer evolution stages correlate with the bends in the curve of the integral heat flow Q at the CF as a function of temperature or Gr. This feature is typical of flow modes in various -scale setups of the classic Czochralski technique with significantly different absolute dimensions of the crucible/ melt/crystal system. This integral heat flux vs temperature difference curve pattern is typical of convection in horizontal fluid layers heated from below [12,42] (these are the so-called discrete transitions when the integral heat flux depends on the temperature difference during convection in horizontal layers of liquid heated from below). The abovementioned bends are caused by a stepwise pattern of space-time convective flow arrangement complexity evolution. New instabilities develop sequentially at certain intervals with an increase in ΔT (once the next ΔT i or Gr i threshold is reached). Then either new secondary swirls or a new oscillation harmonic emerge against the background of the already existing flow. Thus the initial mode is molecular heat conductivity through laminar boundary layers following which the effective heat conductivity coefficient increases in a discrete manner due to convective mixing in the normal direction to the CF. Effective diffusion ratios should also change in a similar manner. Understanding of these regularities will probably be useful in analyzing the origins of crystal inhomogeneity.
Another important problem is the dependence of melt flow structure and heat exchange regularities at the CF on the CF shape. Free convection modes can be naturally expected to produce a close to conical front shape. This is confirmed by a real hexadecane crystallization experiment (Pr = 45.6). Figure 7 shows CF shape for established convection at steady-state temperature boundary conditions as a function of the Gr number. It can be seen from Fig. 7 that the volume and mass of the solidifying material decrease with an increase in Gr which is equivalent to an increase in the temperature difference and local heat flows for constant dimensions. As a result the crystallization front becomes increasingly rounded. This can be accounted for by the steady-state melt stratification and the accumulation of the melt overheated relative to the bulk-average temperature in the top part of the fluid volume.
One more problem is associated with the effect of temperature boundary conditions at the crucible bottom. Figure  8 shows flow pattern as a function of the CF curvature and the ratio of the crucible bottom T b and the crucible side wall T w temperatures. The experiments were carried out with independent bottom and wall heating. The left-hand photos show the flow and the right-hand ones present the main flow trajectories providing qualitative and numerical characterization of the meridian flow dimensions and direction. The general flow pattern is initially close to the one that is typical of the abovementioned flat CF shape with side heating only. Crucible bottom heating at T b < T w only increases the fluid circulation rate at the meridian flow contour. An increase in the crucible bottom temperature with small overheating relative to the crucible wall temperature produces a flow structure that is typical of Rayleigh-Benard convection. A cell forms consisting of two torus-shaped swirls with an upward flow at a certain distance from the crucible walls (see Fig. 8d, e). The heat supplied to the crucible bottom is partially transferred to the side wall, not to the CF. For a flat CF the spatial flow pattern development and evolution scenario is similar but the ring-shaped upward stream loses stability, the flow becoming time-dependent.

Conclusion
Buoyancy and thermocapillary driven convection was studied experimentally for melt simulating fluids in thermal hydrodynamic systems similar to the classic Czochralski technique, for laminar to well-developed turbulent flow modes, with 80 and 295 mm crucible diameter models. The experiments proved that an increase in the Grashof and Marangoni numbers induces multiple stepwise qualitative changes in the structure of the initial buoyancy and thermocapillary driven flow. Once the next threshold temperature difference (or the threshold Gr and Ma numbers) is reached the space-time flow arrangement becomes more complex (new secondary swirl flow develops against the background of the main flow or a new oscillation harmonic emerges). The heat exchange regularities change accordingly. The scale factor proved to affect the hydrodynamics, laminar-turbulent transition and heat exchange regularities, yet further studies for fullscale models are required for practical purposes. TGCC modes are distinguished by a significant radial inhomogeneity of local heat flows in laminar, transition and turbulent modes. Model fluid experiments (Pr = 45.6) with real crystallization showed that the CF shape depends largely on the spatial melt flow structure and confirm the possible crystallization front shape hypotheses made on the basis of single-phase convection studies.