Hydrodynamics and heat exchange of crystal pulling from melts. Part II: Numerical study of free convection mode

Vladimir S. Berdnikov

doi:10.3897/j.moem.7.4.51073

Review Article

Hydrodynamics and heat exchange of crystal pulling from melts. Part II: Numerical study of free convection mode

Vladimir S. Berdnikov^‡

‡ Kutateladze Institute of Thermophysics, Siberian Branch of the Russian Academy of Sciences, Novosibirsk, Russia

Corresponding author: Vladimir S. Berdnikov ( berdnikov@itp.nsc.ru )

This is an open access article distributed under the terms of the Creative Commons Attribution License (CC BY 4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.

Citation: Berdnikov VS (2021) Hydrodynamics and heat exchange of crystal pulling from melts. Part II: Numerical study of free convection mode. Modern Electronic Materials 7(4): 151-165. https://doi.org/10.3897/j.moem.7.4.51073

Abstract

This work is a brief overview of numerical study results for hydrodynamics and convective heat exchange for the classic Czochralski technique obtained at the Institute of Thermophysics, Siberian Branch of the Russian Academy of Sciences. Flow structure evolution has been compared for melts with Prandtl numbers Pr = 0.05, 16, 50 and 2700. The regularities of local and integral convective heat exchange in the crucible/melt/crystal system for thermal gravity, thermocapillary and heat-induced gravity capillary convection have been studied. The calculations have been carried out using the finite and compact difference methods.

This work is a continuation of an article: Berdnikov VS (2019) Hydrodynamics and heat exchange of crystal pulling from melts. Part I: Experimental studies of free convection mode. Modern Electronic Materials 5(3): 91–100. https://doi.org/10.3897/j.moem.5.3.46647

Keywords

crystal growth, Czochralski technique, melt hydrodynamics, heat exchange, thermal gravity capillary convection, numerical simulation.

1. Introduction

Below is a continuation of the overview of experimental and numerical results for crystal pulling from melt using the Czochralski technique [1] obtained at the Institute of Thermophysics, Siberian Branch of the Russian Academy of Sciences. The physical origin of free convection in melts during crystal growth may vary. In gravity field affected systems the cause of buoyancy forces (Archimedes forces) can be the dependence of the density of fluids (melt, solution or solution in melt) on temperature or component concentrations. Providing the required the crystal growth conditions in the Czochralski technique implies maintaining the required temperature difference between the crystallization front (CF) and the crucible walls [2, 3]. As a result, by analogy with any other non-isothermal system containing a viscous media in a gravity field, thermal gravity convection (TGC) which emerges inevitably is fundamentally unavoidable and can hardly be controlled [1, 4–7]. The presence of a free melt surface portion between the CF and the crucible walls leads to a combined effect of the buoyancy forces (Archimedes forces) and the thermocapillary effect [1, 4–7]. It is therefore thermal gravity capillary convection (TGCC) in melt that is usually the initial mode for the search of optimum crystal growth conditions, and process engineers should understand its main features as they determine the quality of the growing crystals [8]. Importantly, if free convection mode is superimposed with crystal rotation, the relative contributions of free and forced convection are judged about based on the ratio of the Grashof and the Reynolds numbers Gr/Re². One fine point is that reference books on heat exchange provide typically heat exchange process parameters as a function of Gr/Re² but, along with this ratio, one should also know the absolute values of Gr and the Marangoni number Ma and hence the initial TGCC mode parameters which are strongly dependent on Gr and Ma.

One aim of this work is to illustrate the specific features of the expression of the two factors (buoyancy forces and thermocapillary effect) which initiate melt motion and the contribution of these factors in case of their joint action. In TGC mode and in flows of mixed origin the temperature and velocity fields are self-consistent under a significant effect of the buoyancy forces. The presence of the thermocapillary effect further strengthens the feedback between hydrodynamics and heat exchange. Separate experimental study of the action of the buoyancy forces and the thermocapillary effect in Earth conditions is an extremely complex or even impossible task. This greatly increases the importance of numerical studies. They open new opportunities for understanding the action and relative contributions of bulk and surface forces to convection flow generation and intelligent search for methods of controlling their contributions to the formation of the required flow structure. Furthermore, numerical studies based upon adequate physical and mathematical models and numerical methods allow studying the evolution of local parameters of the velocity and temperature fields at a higher resolution. On the one hand they significantly extend experimental results and on the other, provide fundamentally new information. For TGCC modes this primarily concerns understanding the relative contributions of the gravity and capillary effects. For sufficiently high solute concentrations in melt, concentration-related effects add to the thermal ones. This is also true for buoyancy force since density depends on temperature and solute concentration and therefore the thermocapillary effect will combine with the concentration-capillary one. Analysis of joint heat and mass transport processes dramatically complicates the initial problem. The emergence thresholds and types of flow instabilities [9] and the melt to solid transition conditions may change fundamentally. Below we will only consider the «simple» case of non-isothermal single-component systems. We will comparatively analyze the evolution of flow structures for fluids with Prandtl numbers Pr = 0.05, 16, 50 and 2700 with an increase in the characteristic temperature difference in the crucible/melt/crystal system. The calculations [10–12] just like the experiments [1] were carried out for Pr orders of magnitude covering the range from the values typical of liquid metals and semiconductors to those typical of highly viscous melts of oxides. Comparative analysis of the physical and numerical experimental data allows solving an important methodological problem – determining the range of parameters in which experimental results obtained for fluids with greater Pr can be safely used for the treatise of processes typical of liquid metals and specifically of silicon. Transparent media with 5 ≤ Pr ≤ 3000 allow hydrodynamic studies but differ greatly from liquid metals with Pr ~ 0.01. The work was aimed at revealing a correlation between the local features of spatial pattern and structure of melt flows and local heat flows. The effect of simulation area geometry and relative dimensions on hydrodynamics and heat exchange in melts was studied. The study was carried out for an axial symmetry problem statement and laminar flow modes. Free or mixed laminar convection are the initial states for studies of laminar-turbulent transition (LTT) in melt flow modes differing by physical origins. However they have not yet been studied to a sufficient extent for understanding the origins of LTT, i.e., the physical background of the development of different instability types is not completely clear. For the system in question the flow will be stable only in a limited range of dimensions and temperature differences, even for perfectly symmetrical boundary conditions and ideally elaborated control methods used (e.g. if the temperature control system does not produce disturbance). Determination of theoretical stability limits is another problem of physical and numerical simulation [1, 13–16].

2. Numerical study problem statement

The ideal physical-mathematical model of the processes studied for the crucible/melt/crystal system is TGCC in the vicinity of cooled discs with different diameters partially covering the free surface of fluid in a fixed crucible with uniformly heated walls. The crucible with the melt has a cylindrical symmetry and fixed rectilinear boundaries. The simulation area diagram is shown in Fig. 1. The problems are solved for an axial symmetry statement, and therefore Fig. 1 only shows the right-hand part of the melt: from the symmetry axis to the right-hand crucible wall. The calculated velocity and temperature field parameters are also shown for that region. Adhesion conditions are satisfied for the tangential velocity components at the crucible and CF boundaries in the numerical model which conforms to the model developed earlier [5, 17]. In a general case, tangential component balance condition is set at the free melt surface for the force caused by the surface tension gradient and for the friction forces:

Here s is the melt surface tension coefficient depending on the temperature T, σ = σ₀(1 – β_σΔT), where β_σ = (1/β₀)(–∂β/∂T); and β₀ is the surface tension coefficient at the characteristic system temperature. At all the boundaries of the area the normal velocity component is equal to zero, i.e., the rigid walls satisfy the impermeability conditions and the free melt surface meets the condition of no free melt surface deformation by convective flow. Furthermore the possibility of meniscus formation at the crystal edge and crucible walls is ignored at this stage [18–21]. These approximations mostly apply to thermocapillary convection (TCC) modes for which it is accepted that the area boundaries remain the same (planar) but the system is under theoretical zero-gravity conditions. The crystal T₀ and the crucible side wall T_w temperatures are set to be different but constant. Therefore the characteristic temperature difference in the system is ΔT = Т_w – Т₀. The heat insulation condition is set for the crucible bottom and the free melt surface.

The initial set of dimensionless convection equations in the Boussinesq approximation and in the motion and temperature field axial symmetry assumption was solved for the swirl ω, stream function ψ and temperature T variables using the implicit difference scheme. The stream function equation was solved using the fast Fourier transform variable splitting method with central difference approximation. The Peaceman–Rachford alternating direction difference scheme was used for approximating the heat and impulse transport equations. The convective equation terms for each next half-step were taken from the previous half-step. The equations were reduced to sets of equations with three-diagonal matrices and solved by passing. The optimum grid parameters were determined preliminarily [10]. The optimization criterion was disbalance between the integral dimensionless heat emission coefficients at the crucible side wall and at the simulated CF. The results were obtained mainly for 160 × 160 or 320 × 320 uniform grids. Iteration was stopped after the step swirl increment at nodes became sufficiently small: (|δω_ij|/max|ω_ij|) ≤ (0.1÷3.0) ∙ 10^-8. For liquid metal with the Prandtl number Pr = 0.05 the TGC set of equations was solved using the compact difference method in non-uniform structured grids with a high (fifth order) accuracy [22]. The minimum step for each of the spatial coordinates (r, z) was chosen at h = 0.004 or 0.002. This corresponds to 250 × 175 and 500 × 350 cell grids, respectively.

Laminar convection of different origins was studied: thermal gravity, thermocapillary and thermal gravity capillary [10–15].

The dimensionless set of equations for TGCC was as follows;

with the following boundary conditions:

– crucible bottom

– side surface

– free surface

– crystallization front

– symmetry axis

In TGC mode the no friction and deformation conditions are set for the free melt surface. Further simplification is the boundary heat insulation (adiabaticity) condition:

The identity criteria are introduced into the dimensionless set of equations and boundary conditions: the Grashof number Gr = (βg/ν²)ΔTR_sc³ or the Rayleigh number Ra = GrPr = (βg/αν)ΔTR_sc³ which determine TGC intensity, the Marangoni number Ma = [(–∂β/∂T)/aμ]ΔTR_sc which describes the contribution of the thermocapillary effect and the Prandtl number Pr = ν/a which describes the thermophysical properties of the melt. Here t is the time, U and V are the radial and axial velocity components, θ = (T – T₀) /ΔT, T is the temperature, T₀ is the temperature at the CF, ΔT is the temperature difference, H is the melt layer height, β is the volume expansion coefficient, μ and ν = μ/ρ are the dynamic and kinematic viscosity coefficients, respectively, a = λ/ρC_P is the temperature conductivity coefficient, λ is the heat conductivity coefficient, ρ is the density, C_P is the constant pressure heat capacity and R_sc and R_c are the crystal and crucible radii, respectively. The following scales were used for the transition to dimensionless equations: R_sc for length, ΔT for temperature and ν/R_sc for velocity. These scales allow all the results presented below in a dimensionless form to be brought to a dimensional form. The stream function and swirl are set by the following relationships: U = (1/r)(∂ψ/∂z); V = –(1/r)(∂ψ/∂r); ω = ∂U/∂z – ∂V/∂r.

Figure 1.

Simulation area diagram for ideal physical-mathematical crystal growth model: (1) crucible bottom, (2) crucible side wall, (3) free melt surface, (4) crystal and (5) symmetry axis.

3. Results and discussion

3.1. Thermal gravity convection mode

TGC was studied for no friction conditions at the free fluid surface. The fluid flow structure evolution with an increase in Gr was studied for preset geometry and Pr. Temperature fields were calculated for TGC in the initial heat conductivity mode the isotherms for which are presented in Fig. 2. Figure 3 shows calculation results for stream function isocountour fields and isotherms for a fluid with Pr = 16 and a constant geometry chosen in cooperation with the Ioffe Institute for Problems of Mechanics of the Russian Academy of Sciences [5, 17] for test experiments: H/R_sc = 0.7 and R_c/R_sc = 2.76. Figure 3 and similar figures below (Figs 8, 9, 14 and 17) only show the right-hand part of the axial symmetrical isocountour fields. Pr = 16 corresponds to ethyl alcohol which was used in the physical experiment as the melt simulator and is close to Pr for BGO (bismuth orthogermanate), paratellurite and garnet melts. The results for boundary conditions corresponding to the absence of friction at the free surface of fluids are of practical interest for oxides with weak surface tension coefficient temperature dependences (e.g. yttrium aluminum garnet with Pr = 8) or high viscosity. Figure 3 a–d shows stream function isocontours which for laminar stationary flow mode coincide with fluid motion trajectories. It is therefore safe to hypothesize about the evolution of the spatial flow pattern with an increase in Gr. One can therefore compare the calculated spatial flow pattern with the experimental one [1]. The flow pattern in Fig. 3 (and in similar figures below) is counterclockwise. Fig. 3 e–h shows isotherms whose shape as a function of Gr allows predicting CF curvature variation under real conditions. Comparison of Fig. 2 and Fig. 3 e shows that at Gr = 10, i.e., in low-intensity convection mode (low-velocity creep flow), the isotherm pattern already differs considerably from that for heat conductivity mode due to hot fluid entrainment to the top part of the fluid volume. At Gr = 100 the isotherms in the fluid core became horizontal. Thus at practically negligible temperature differences the system already clearly exhibits the formation of a stable stratified core. With further Gr growth this leads to a transition to boundary layer mode, i.e., a stagnant or low-mobility core forms, occupying large volume in the central crucible region, and fluid circulation around the perimeter of that region becomes increasingly intense with an increase in Gr. This entails a redistribution of the conductive heat flow component: as a result at Gr < 100 in heat conductivity mode the molecular heat flows are directed from the crucible walls toward the CF while at Gr ≥ 100 thermal interaction between radial branches of meridian flow starts: part of the molecular heat flow is directed downward (see Fig. 3 f–h). The quality criterion of temperature field reconstruction can be the magnitude and direction of the temperature gradient in the core [10].

Figure 2.

Isotherm field in heat conductivity mode (H/Rc = 0.7, R_c/R_sc = 2.76).

Figure 3.

Evolution of (a–d) stream function isocontours and (e–h) isotherms with increasing Gr (Pr = 16, Ma = 0, H/R_c = 0.7, R_c/R_sc = 2.76): (a, e) Gr = 10¹; (b, f) 10²; (c, g) 10³; (d, h) 10⁴.

Along with the isotherms, additional temperature field characteristics are temperature profiles in horizontal and vertical sections which are measured in the physical experiment [1]. The vertical temperature distribution in the r = (R_c + R_sc)/2 section at different Gr is presented in Fig. 4. These profiles provide additional information on the regularities of melt core stratification with an increase in Gr. First, comparison between Curves 1 and 2 show a dramatic change in the temperature distribution pattern from an unstable state in heat conductivity mode (the bottom melt portion is overheated relative to the average temperature) to a stable state in convection mode. Secondly, an increase in Gr leads to an increase in the core temperature. One can see this by comparing the profiles for the transition from Curve 2 to Curve 4. However, simple regularities only hold until Gr ≤ 400. At higher Gr the melt becomes colder near the crucible bottom relative to the average temperature while remaining overheated at the top. The flow structure becomes increasingly complex at Gr ≥ 500: two «rotation centers» form (Fig. 3 d) and the horizontal dimension of the poorly mixed melt core increases (Fig. 3 c and d). The change in the flow intensity can be seen from Fig. 5 showing dimension profiles of the radial and axial velocity components. For recalculation the crucible radius was taken to be 50 mm and the ethyl alcohol properties were taken for T = 20 °C. A specific feature of the velocity profiles is as follows: an important factor in the R_sc ≤ r ≤ R_c region at the top flow branch is that the top boundary is free and without friction (see Fig. 5 b, section r = (R_c + R_sc)/2), the velocity being the highest there. For 0 ≤ r ≤ R_sc the radial-convergent flow first impinges onto the crystal edge, and then a dynamic boundary layer forms downstream under a rigid wall (see Fig. 5 a, section r = R_sc/2). The hot melt impinges under the cold CF surface. An unstable stratified boundary layer forms downstream as can be seen from the isotherms in Fig. 3 f–h. A flow tear-off zone is at a certain distance from the crystal edge: a cold downward fluid stream forms there. Importantly, at all the stream development stages (tear-off, acceleration and deceleration) in the bottom region, fluid stratification in the melt core becomes increasingly expressed with an increase in Gr. The cold downward stream impinges onto the relatively cold fluid sublayer near the bottom, and radial stream oscillations may form as a result according to earlier experiments [6, 7].

Figure 4.

Temperature distribution in layer height in the r = (R_c + R_sc)/2 section for TGC mode at Pr = 16, H/R_c = 0.7 and R_c/R_sc = 2.76: (1) heat conductivity mode, (2) Gr = 10; (3) 20; (4) 100; (5) 500; (6) 1000.

Figure 5.

(a, b) radial U and (c) axial V velocity components at Pr = 16, H/R_c = 0.7 and R_c/R_sc = 2.76 for different Gr: (1) 10; (2) 20; (3) 50; (4) 100; (5) 200; (6) 300; (7) 500; (8) 1000.

Figure 6.

Radial distributions of local heat flows in TGC modes at Pr = 16, H/R_c = 0.7 and R_c/R_sc = 2.76: (1) heat conductivity mode, (2) Gr = 10; (3) 10²; (4) 10³; (5) 10⁴; (6) 78274.

Figure 7.

TGC structure dependence: (a–d) stream function isocontours and (e–h) isotherms as a function of Pr at preset Gr = 100: H/R_c = 0.7 and R_c/R_sc = 2.76: (a, e) Pr = 0.05; (b, f) 16; (c, g) 50; (d, h) 2700.

At the periphery of the cold downward fluid stream descending from the crystal to the crucible bottom the relatively overheated fluid is entrained (due to the viscous flow) downward from the core with an increase in Gr (see Fig. 3 g and h). In the outer parts of the upward flows at the crucible walls, on the contrary, the cold fluid is entrained upward. This generates a reversing force due to excess buoyancy which first generates a flow reflected from the crucible bottom (clearly seen in Fig. 3 d at Gr = 10⁴). Secondly the interaction of the oppositely directed viscous friction force and buoyancy force may develop «droplet» type instability (periodical regeneration of the outer part of the boundary layer) which was observed experimentally [1], or cause core oscillations. This will cause temperature oscillations under the CF.

Temperature profiles are used for calculating the radial distributions of the local heat flow at the CF (Fig. 6) which can be easily represented as dimensionless temperature gradients. The local heat flow has an expressed maximum at the crystal edge for all the above modes. Screening of the cold surface by the gradually cooling fluid which is typical of laminar flow mode occurs downstream initially at Gr ≤ 100 resulting in a decrease of the local heat flow from the CF edge to the center. This flow is less intense than that for heat conductivity mode (see Fig. 6, Curves 1 and 2). However with an increase in Gr and the flow velocity the thicknesses of the heat (see isotherms in Fig. 3 f and h) and dynamic (see Fig. 5 a) boundary layers decrease while the local temperature gradients and heat emission increase. It is typical that the local heat flows simultaneously become more homogeneous over increasingly large CF area (see Fig. 6, Curves 5 and 6). Along with convective heat transfer, there is a large contribution of the increasingly expressed stable temperature stratification of the melt which can be seen from Gr = 10⁴ isotherms in Fig. 3. Thus, as a result of an increase in the temperature difference in the system the CF will be located in the melt that is increasingly overheated relative to the average temperature.

Calculations for Pr = 0.05, 50 and 2700 were carried out for the same geometry (Figs 7 and 8). Figure 7 shows results of calculation for a preset Gr number which is the ratio of the buoyancy and viscous friction forces (or the respective work). Comparison of the meridian flow and isotherms for the same Gr and different Pr shows that at small Gr the stream function isocountour fields and isotherms already differ substantially. Heat exchange of the liquid metal with Gr = 100 is almost the same as for heat conductivity mode. The liquid with Pr = 2700 exhibits the largest change in the flow structure, i.e., with increasing Pr the system deviates increasingly from heat conductivity mode to well-developed boundary layer mode. With Pr growing in the 16 ≤ Pr ≤ 2700 range the isotherm pattern changes increasingly due to hot fluid entrainment to the top part of the volume and the formation of a stable stratified stagnant core. Figure 8 shows data for the same Rayleigh number but different Pr numbers. The Rayleigh number takes into account, along with the abovementioned parameters, that the floating unit volume of melt spends energy not only for overcoming viscous friction but also for equalizing its temperature with the ambient temperature due to molecular heat conductivity. Excess buoyancy force is then lost. Comparison of the data shown in Figs 7 and 8 suggests that an increase in Pr at the same Gr produces a flow structure reconstruction pattern which is to similar that for a Gr growth at a constant Pr. Figure 8 shows respective Gr for each Pr at a preset Ra. It can be seen that allowance for two dissipation factors leads to a far greater “identity degree” in system behavior. At a constant Rayleigh number, e.g. at Ra = 500, the flow patterns and isotherm fields are almost coincident for Pr = 50, 16 and 2700. Fluids with higher Pr have high viscosity and relatively low heat conductivity. At close Gr, floating fluid components in these fluids do not lose buoyancy and overcome friction forces almost as efficiently. The isotherm shapes at Pr = 0.05 and the same Ra = 500 are close to those at other Pr but the flow structure is different. In liquid metals with high heat conductivity, emergence of flows of similar patterns and generally similar temperature fields and local heat flow distributions require Gr numbers that are 2–4 orders of magnitude greater than those required for media with Pr ≤ 5 (Figs 8–12). Isotherms for Pr = 0.05 retain patterns close to those typical of heat conductivity mode almost up to Gr = 2 ∙ 10³. Melt core isotherms for 2 ∙ 10³ ≤ Gr ≤ 10⁴ acquire horizontal sections, this regularity being traced in radial temperature distributions (see Fig. 10 b), heat flow distributions (see Fig. 12) and integral heat flow vs Gr curves (Fig. 13). Radial heat flow distributions for Gr ≤ 10³ are almost coincident with that distribution for heat conductivity mode (see Fig. 12). Thus, convection largely affects the temperature distribution (see Fig. 10) and heat exchange in liquid metal beginning from Gr ≥ 2 ∙ 10³. Upon further Gr growth the ψ and T isocontour shapes are generally the same as those for large Pr. Velocity profiles vary in a similar manner with Gr growth for well-developed convection (see Figs 5 and 9) and a boundary heat layer forms in a similar manner as can be clearly seen from Fig. 11. Thus many behavior features of liquid metal systems including silicon melt can be simulated using transparent melt simulators, but taking into account the above regularities.

Lowering of the liquid metal melt level in the crucible leads to an increase in the contribution of molecular heat conductivity to heat exchange between the crucible walls and the CF (Figs 13 and 14). The reconstruction of local properties of the temperature and heat flow fields that is typical of convective modes occurs during melt level lowering at high Gr.

Figure 8.

TGC structure dependence: (a–d) stream function isocontours and (e–h) isotherms as a function of Pr at preset Ra = 500: H/R_c = 0.7 and R_c/R_sc = 2.76: (a, e) Pr = 0.05, Gr = 10⁴; (b, f) Pr = 16, Gr = 35.25; (c, g) Pr = 50, Gr = 10; (d, h) Pr = 2700, Gr = 0.185.

Figure 9.

(a) radial U and (b) axial V velocity components in TGC mode at Pr = 0.05, H/R_c = 0.7 and R_c/R_sc = 2.76: (1) Gr = 4.7 ∙ 10²; (2) 4.7 ∙ 10³; (3) 4.7 ∙ 10⁴; (4) 4.7 ∙ 10⁵.

Figure 10.

Temperature distribution in layer height in (a) r = (R_c + R_sc)/2 section and (b) radial direction for the z = H/2 section in TGC mode at Pr = 0.05, H/R_c = 0.7 and R_c/R_sc = 2.76: (1) Gr = 4.7 ∙ 10²; (2) 4.7 ∙ 10³; (3) 4.7 ∙ 10⁴; (4) 4.7 ∙ 10⁵.

Figure 11.

Temperature distribution in layer height for TGC mode in sub-crystal region at Pr = 0.05, H/R_c = 0.7, R_c/R_sc = 2.76 and Gr = 2.3 ∙ 10⁶: (1) r = 0; (2) 0.25; (3) 0.5; (4) 0.75; (5) 1.

Figure 12.

Radial distribution of local heat flows in TGC mode at Pr = 0.05, H/R_c = 0.7 and R_c/R_sc = 2.76 for different Gr: (1) 0; (2) 4.7 ∙ 10²; (3) 4.7 ∙ 10³; (4) 4.7 ∙ 10⁴; (5) 4.7 ∙ 10⁵.

Figure 13.

Evolution of (a, b, e) isotherms and (b, d) stream function isocontours for transition from (a) heat conductivity mode to (b, c) TGC at Pr = 0.05, H/R_c = 0.25 and R_c/R_sc = 2.76:(b, c) Gr = 2.3 ∙ 10⁵; (d, e) 2.3 ∙ 10⁶.

Integration of the local heat flows over the CF area yields the full heat flow from the crucible walls to the melt which is supplied to the CF due to heat conductivity and convection. Instead of local or integral heat flow, the dimensionless heat emission coefficient, or the Nusselt number, is conventionally calculated as a function of dynamic problem parameters (in the case in question, as a function of the Gr or Ra numbers, Fig. 15): Nu = αR_sc/λ = (∂T/∂z)_ctR_sc/ΔT (here α is the dimensional heat emission coefficient and (∂T/∂z)_ct is the temperature gradient near the CF). The Nusselt number characterizes the factor by which convective mode heat emission is higher than it would be at the same temperature difference in heat conductivity mode. It is based on the simple local heat flow expression q = –λ(∂T/∂z)_ct = αΔT. Integral parameters are determined by calculating average temperature gradients and heat emission coefficients. Figure 15 shows that local and integral TGC parameters calculated in the 16 ≤ Pr ≤ 2700 range can well be generalized if the Rayleigh number is used as the dynamic parameter. The convection contribution to heat exchange becomes noticeable at Pr = 0.05 and grows obeying almost the same law as for 16 ≤ Pr ≤ 2700 but at a sufficiently high convection flow intensities, i.e., at Grashof numbers being three orders of magnitude greater than those required for media with 16 ≤ Pr ≤ 2700. The Nu numbers for Ra = 500 in fluids with Prandtl numbers of 16, 50 and 2700 are almost the same and equal to 1.76, and for Pr = 0.05, Nu = 1.27. At Ra = 5 ∙ 10⁴, the differences in the Nu values become noticeable: Nu = 4.95, 5.07 and 5.15 for Pr = 16, 50 and 2700, respectively, and Nu = 3.21 for Pr = 0.05. The dependence Nu = (0.81 ± 0.08)Gr^0.227±0.01 holds with a high accuracy in the 0 ≤ Gr ≤ 2500 range at Pr = 16. Integral heat exchange data can be generalized with a high accuracy by the following Nu(Gr, Pr) dependence: Nu = (0.55 ± 0.035)Pr^0.206±0.009Gr^0.2 in the 0 ≤ Gr ≤ 2500 range. The data shown in Figs 9–14 for Pr = 0.05 were obtained using the compact difference method (high order accuracy schemes [20]) and those in Fig. 15, using the finite and compact difference methods. The heat flows calculated using the finite difference method for a 160 × 160 grid in the 0 ≤ Gr ≤ 2500 range are almost the same for the heat emitting and the heat receiving surfaces, the disbalance being less than 6% for Gr = 8 ∙ 10⁴.

Figure 14.

Temperature distribution in layer height for TGC mode in sub-crystal region at Pr = 0.05, H/R_c = 0.25, R_c/R_sc = 2.76 and Gr = 2.3 ∙ 10⁶: (1) r = 0; (2) 0.25; (3) 0.5; (4) 0.75; (5) 1.0.

Figure 15.

Dimensionless heat emission coefficient as a function of Rayleigh number: (1) Pr = 0.05, calculation using finite difference method; (2) Pr = 0.05, calculation using compact difference method; (3) Pr = 16, calculation using finite difference method; (4) Pr = 16, Ma = 34240; (5) Pr = 50; (6) Pr = 2700.

3.2. Thermal gravity capillary convection mode

Figures 2, 3, 7, 8 and 13 show that there is a large temperature gradient along the free melt surface at any Gr and Pr in both heat conductivity and TGC modes. Therefore if the melt surface tension coefficient is largely temperature-dependent then the formation of the flow structure is controlled both by the buoyancy (Archimedes) forces and the thermocapillary effect. Analysis of Figs 7 and 13 also suggests that an increase in Gr in TGC mode causes hot melt flow to the CF edge at a temperature gradually approaching the crucible wall temperature. The temperature of the most part of the free melt surface equalizes, the maximum temperature gradient being at the CF edge. It is natural to expect that the relative contribution of the thermocapillary effect will vary in accordance with this tendency with an increase in Gr.

First of all, TCC was studied at Pr = 16 in a system having the same geometry as the one for the TGC study (Fig. 3) with an increase in ΔT in the range corresponding to Marangoni numbers of 10 to 10⁵ [11, 14]. Then the basic mode was accepted to be TTC at Ma = 34240 corresponding to the experimentally implemented combined TGCC mode (Fig. 16 a). The initial thermocapillary flow (the system is in the initial theoretical zero-gravity state) was superimposed with the growing gravity field and the respective growth in Gr (Fig. 16 b–d). Thus the flow structure evolution (Fig. 16) and the local parameters (Figs 17 and 18) were studied for an increase in the relative contribution of buoyancy forces at constant Ma = 34240.

Figure 16.

Evolution of (a–d) stream function isocontours and (e–h) isotherms with increasing Gr at constant Ma = 34240, H/R_c = 0.7, R_c/R_sc = 2.76 and Pr = 16: (a, e) Gr = 0; (b, f) 10³; (c, g) 3 ∙ 10³; (d, h) 10⁴.

Figure 17.

Temperature distribution in layer height for TGCC mode in sub-crystal region at Pr = 16, H/R_c = 0.7, R_c/R_sc = 2.76, Gr = 41270 and Ma = 34240 with equal radial intervals (1–14) between (1) r = 1 mm and (14) r = 0.

Figure 18.

Radial distribution of local heat flows in (a) TGC and (b) TGCC modes at Pr = 16, H/R_c = 0.7 and R_c/R_sc = 2.76: (a): (1) heat conductivity mode, (2) Ma = 10; (3) 100; (4) 10³; (5) 10⁴; (6) 10⁵; (b): Ma = 34240; (1) Gr = 0; (2) 10³; (3) 10⁴; (4) 41272.

TCC is theoretical zero-gravity mode. The flow results from the force that acts intensely only along the free surface (due to the thermocapillary effect). The fluid flow is directed along the free surface from the crucible wall to the crystal edge. The fluid flows to below the CF by inertia and due to the fact that fluid from the volume is pulled up to the surface along the crucible walls in accordance with the continuity conditions. Under the crystal beginning from its edge the flow is decelerated by viscous friction. In TCC modes at Ma ≤ 1000 the rotation center formed by the closed stream function isocontours is approximately in the area with the coordinates r = (R_c + R_sc)/2, z = 2H/3 (see Fig. 4 [14]). As Ma increases the rotation center shifts gradually toward the top right corner (see Fig. 16 a). If compared with TGC mode the cooled melt stream from the CF to the crucible bottom becomes less intense and has a significantly greater diameter (see Figs 3 and 16 a). If even a weak gravity field is present (see Fig. 16 b) a compact cold downward stream forms flowing from under the cooled CF surface, by analogy with the data shown in Fig. 3. Hot melt from the crystal edge is pulled into the downward stream toward the center. As Gr approaches the experimental value Gr = 41270 the stream function isocontour shape becomes increasingly similar to the experimental [15] visualized fluid motion trajectories (Fig. 19):

– circulation over the external contour and integration into the main superficial swirl having a thermocapillary origin;

– emergence of a bottom-reflected flow from which part of fluid is entrained into the superficial swirl while passing through the central part of the melt core;

– a stagnant region clearly seen in Fig. 19 forms under the superficial swirl at Gr = 5000.

The flow pattern shown in Fig. 16 d is generally in a good agreement with the experimental one [1]. At Ma/Gr ratios close to those calculated from the temperature difference set in the experiment, the qualitative local parameters are also quite close to the experimental data, the absolute velocity approaching the experimental ones (see e.g. Figs 7–9 of earlier work [13]).

Figure 19.

(a) Experimental fluid motion trajectories and (b) calculated stream function isocontours for similar conditions: H/R_c = 0.7, R_c/R_sc = 2.76, Gr = 5000 and Ma = 34240.

Flow regularities affect the radial heat flow distribution at the CF. Figure 18 a shows change in flow regularities for Mg growth in TCC mode compared with heat conductivity mode. The melt flow concentrated at the free melt surface efficiently heats the crystal edge. The local heat flow decreases rapidly toward the center. In TGCC mode the growing contributions of buoyancy forces and TGC affect the most noticeably the growth of heat flows in the 0 ≤ r ≤ R_sc/2 range (see Figs 17 and 18 b). It is safe to affirm that TGC completely controls local heat exchange in this region. Hot melt moving from the crystal edge to the center is entrained into the downward stream formed by buoyancy forces. Therefore the local heat flow in the axial region grows with Gr while at the edge it remains almost the same as in TCC mode (see Fig. 18 b).

Studies of the temperature distribution regularities in the volume and along the free fluid surface in TCC, TGC and then TGCC modes [11–15] showed that the thermocapillary effect contribution is mainly resolvable against the TGC background in the vicinity of the cold crystal edge. Hot fluid entrainment toward the free surface from the heated crucible walls and temperature-stable fluid stratification in the volume at most part of the free surface eliminate the longitudinal temperature gradient, also due to the additional effect of TCC. This fact and hence the smaller thermocapillary effect contribution to the intensification of the general meridian flow affect the dependence of integral heat exchange on Gr at preset Ma. Curve 4 in Fig. 15 shows Nusselt number Nu as a function of Gr at preset Ma = 34240. At Gr = 0, Nu = 3.80 and is completely TCC-controlled. As can be seen from Fig. 15, Gr growth causes an asymptotic Nu tendency to the TCC-controlled value: Nu = 8.823 at Gr = 21000 and Ma = 34240 in TGCC mode, while Nu = 8.292 at the same Gr in TGC mode. The above data confirm that TCC has a clearly superficial pattern. The thermocapillary effect causes hot liquid “pulling” toward the crystal edge and a dramatic growth of the local heat flow at the edge (see Fig. 18 a). Unlike TGC mode the fluid is considerably decelerated in the axial region (see Fig. 16 a) and contributes considerably to the heat flow at the crystal edge, but affects the sub-crystal region only slightly. Thermocapillary flow has a high velocity along the free surface, i.e., it is higher than in TGC mode. Local U (r) maxima correlate with radial temperature gradients across the surface [10]. The influence of the thermocapillary effect shows itself against the TGC background in an increase in the velocity at the free surface, but this influence decreases because of temperature equalization over the free surface due to gravity convection and temperature stratification of fluid in the gravity field [10, 14]. This regularity of flow structure can be traced in the radial distributions of the local heat flow (see Fig. 18) and the integral heat flow as a function of Gr (see Fig. 15). The local heat flow velocity from the edge to the center decreases much faster (see Fig. 18 b, Curve 1 at Gr = 0) than in TGC mode. The superimposition of the thermocapillary effect with TGC causes an abrupt growth of the local heat flow at the crystal edge (almost threefold for media with Pr = 16 and 45.6 [14]). This local effect is observed in a wide range of Ma/Gr ratios. However with a decrease in the latter ratio the TCC contribution to the integral heat flow decreases monotonically as can be clearly seen from Fig. 15. Radial velocity and temperature distributions across the free surface from the crystal edge to the crucible wall were reported [10]. The patterns of the velocity and temperature component profiles in the radial and vertical sections are in a general agreement with the experimental results. The order of magnitude of the velocity components is also similar to the experimental one.

4. Conclusion

Thermal gravity, thermocapillary and gravity capillary modes of laminar convection were numerically studied. The calculations were carried out for the 0 ≤ Gr ≤ 10⁷ Grashof number range and the 0 ≤ Ma ≤ 10⁵ Marangoni number range. Regularities of the evolution of the general flow structure and local and integral heat exchange were analyzed for each convection type with an increase in the Grashof and Marangoni numbers. The spatial pattern of the axial symmetry laminar flow observed experimentally in combined TGCC mode was reproduced, testifying to a high authenticity of the results. Analysis of the data reported herein is of practical interest. At moderate and high Pr (10 ≤ Pr ≤ 3000), the isotherm shape changes noticeably in comparison with heat conductivity mode already at weak convection Gr = 10 due to hot fluid entrainment toward the top part of its volume. At practically negligible temperature differences the system already clearly exhibits the formation of a stable stratified core. This is a consequence of the transition from low-intensity creep melt flow at very low Gr ≤ 20 to boundary layer mode, i.e., a stagnant melt core forms and fluid circulation around the contour of the region becomes more intense with an increase in Gr. The experimentally observed delay in the reconstruction of the spatial flow shape at Pr = 0.05 with an increase in Gr shows itself in that the dimensionless heat emission coefficient values at the CF in the liquid metal system, which are similar to those values for Pr ≥ 10, are reached at three orders of magnitude greater Gr. Taking into account this delay one can simulate multiple qualitative behavior features of liquid metal systems including silicon melt using transparent melt simulators. For single crystal growth with predominant natural convection, the specific features of temperature field formation will make the CF shape close to conical (see Fig. 8 [14]). However with an increase in Gr the conical angle will grow and the volume of the solidifying material will decrease due to increasing temperature stratification and growth of heat flows at the CF. This was confirmed by physical experiments [1, 14, 15]. The shape and density of the isotherms allow predicting not only the CF curvature evolution tendency but even the solidifying material volume evolution tendency under established heat exchange conditions in the system. This is easy to prove by comparing the qualitative change of the isotherm field in Fig. 3 with the Gr growth and the established CF shape shown in Fig. 7 of an earlier report [1]. Along with the isotherm shape information the data obtained in this work allow one to give a general estimate of the radial distribution of the thermal stresses in the crystal close to the CF. Highest stresses can be expected at the crystal edge onto which the hot flow across the free surface impinges.

Studies of the temperature and velocity distribution regularities in the volume and across the free fluid surface in TCC, TGC and then TGCC modes show that the thermocapillary effect contribution is mainly visible against the TGC background in the vicinity of the cold crystal edge. TCC gives a large contribution to the heat flow at the crystal edge but affects the sub-crystal region only slightly. This fact and hence the smaller thermocapillary effect contribution to the intensifying of the general meridian flow affect the dependence of the integral heat exchange on Gr at preset Ma. Gr growth causes an asymptotic Nu tendency to the TCC-controlled value. Physical prerequisites of the development of two instability types were clarified. The origin of «droplet» type instability that was experimentally observed earlier [4, 6, 7] can be entrainment of the relatively overheated fluid downward from the core at the periphery of the cold downward fluid stream from the crystal to the crucible bottom. This generates a reversing force due to excess buoyancy. This instability type may be the cause of banded crystal inhomogeneity even for the so-called low-gradient technology. The calculated T (z, ri) temperature profiles in the sub-crystal region (see Figs 12, 15 and 18) allow tracing the formation of the unstable stratified boundary layer and estimating the formation threshold of the experimentally observed instability (Rayleigh–Benard type instability [6, 7]) from the local Rayleigh number, its critical value for a horizontal layer with one free boundary being 1100.

The great practical value of numerical simulation is the possibility of studying the main evolution tendencies of local parameters with good resolution and completeness that are superior to the experiment’s capabilities. The range of the parameters Gr, Ma ≤ 5000 is not practically achievable for experimental studies with acceptable control accuracy, maintaining of boundary conditions and temperature measurements in the volume of fluids having e.g. Pr = 16. In the physical experiments the results of which were presented earlier [1] the minimum temperature difference was ΔT = 0.41 K, the thermocouple sensitivity being 41.7 μV/K and the measurement accuracy being 1 μV. Thus it is hard to overestimate the role of numerical experiments in this range of parameters.

Acknowledgments

This research was realized under the project III.18.2.5, number of state registration АААА-А17-117022850021-3.

References

1. Berdnikov V.S. Hydrodynamics and heat exchange of crystal pulling from melts. Part I: Experimental studies of free convection mode. Modern Electronic Materials. 2019; 5(3): 91–100. https://doi.org/10.3897/j.moem.5.3.46647

2. Konakov P.K., Verevochkin G.E., Goryainov L.A. et al. Teplo- i massoobmen pri poluchenii monokristallov [Heat and mass transfer upon receipt of single crystals]. Moscow: Metallurgiya, 1971, 240 p. (In Russ.)

3. Shashkov Yu.M. Vyrashchivanie monokristallov metodom vytyagivaniya [Single-crystal growth by the extrusion method]. Moscow: Metallurgiya, 1982, 312 p. (In Russ.)

4. Berdnikov V.S., Borisov V.L., Panchenko V.I. Experimental modeling of hydrodynamics and heat transfer during the growth of single crystals by the Czochralski method. In: Teplofizicheskie yavleniya pri kristallizatsii metallov [Thermophysical phenomena during crystallization of metals]. Novosibirsk, 1982: 77–92. (In Russ.)

5. Polezhaev V.I. Hydrodynamics, heat and mass transfer during crystal growth. In: Freyhardt H.C., Müller G. (eds) Growth and Defect Structures. Crystals (Growth, Properties, and Applications). Berlin; Heidelberg: Springer, 1984; 10: 87–147. https://doi.org/10.1007/978-3-642-69866-8_3

6. Berdnikov V.S., Borisov V.L., Markov V.A., Panchenko V.I. Laboratory simulation of macroscopic transport processes in the melt during the growth of single crystals by drawing. In: Gidrodinamika i teploobmen v tekhnologii polucheniya materialov [Hydrodynamics and heat transfer in the technology of materials]. Moscow: Nauka, 1990: 68–88. (In Russ.)

7. Berdnikov V.S., Panchenko V.I., Soloviev S.V. Teplofizika kristallizatsii i vysokotemperaturnoi obrabotki materialov [Thermophysics of crystallization and high-temperature processing of materials]. Novosibirsk, 1990: 162–199. (In Russ.)

8. Milvidskii M.G. Poluprovodnikovye materialy v sovremennoi elektronike [Semiconductor materials in modern electronics]. Moscow: Nauka, 1986, 144 p. (In Russ.)

9. Turner J.S. Buoyancy effects in fluids. Cambridge University Press, 1973, 408 p. https://doi.org/10.1017/CBO9780511608827

10. Berdnikov V.S., Gaponov V.A., Kovrizhnykh L.S. Thermal gravitational-capillary convection in a cavity with a longitudinal temperature gradient. J. Engineering Physics and Thermophysics. 2001; 74(4): 999–1006. https://doi.org/10.1023/A:1012323810381

11. Berdnikov V.S., Vinokurov V.V., Panchenko V.I., Solovev S.V. Heat exchange in the classical Czochralski method. J. Engineering Physics and Thermophysics. 2001; 74(4): 1007–1014. https://doi.org/10.1023/A:1012375827219

12. Berdnikov V.S., Vinokourov V.A., Vinokourov V.V., Gaponov V.A., Markov V.A. General regularities of convective heat transfer in the crucible melt-crystal system of Czochralski method and their influence on the solidification front shape. Vestnik of the Lobachevsky NSU. 2011; (4-3): 641–643. (In Russ.)

13. Prostomolotov A.I., Milvidskii M.G. Modeling of thermal processes and defect formation during the growth and thermal annealing of dislocation-free silicon single crystals and wafers. Izvestiya Vysshikh Uchebnykh Zavedenii. Materialy Elektronnoi Tekhniki = Materials of Electronics Engineering. 2008; (3): 49–53. (In Russ.)

14. Bergfelds K., Sabanskis A., Virbulis J. Validation of mathematical model for CZ process using small-scale laboratory crystal growth furnace. IOP Conf. Ser.: Mater. Sci. Eng. 2018; 355: 012004. https://doi.org/10.1088/1757-899X/355/1/012004

15. Yunzhong Zhu, Feng Tang, Xin Yang, Mingming Yang, Decai Ma, Xiaoyue Zhang, Yang Liu, Shaopeng Lin, Biao Wang. In-situ detection of convection and rotation striations by growth interface electromotive force spectrum. J. Cryst. Growth. 2018; 487: 120–125. https://doi.org/10.1016/j.jcrysgro.2018.02.010

16. Nikitin N.V., Nikitin S.A., Polezhaev V.I. Convective instabilities in the hydrodynamic model of crystal growth by the Czochral method. Uspekhi mekhaniki. 2003; 2(4): 63–105. (In Russ.)

17. Polezhaev V.I., Bune A.V., Verezub N.A., Glushko G.S., Gryaznov B.L., Dubovik K.G., Nikitin C.A., Prostomolotov A.I., Fedoseev A .I., Cherkasov S.G. Matematicheskoe modelirovanie konvektivnogo teplomassoobmena na osnove uravnenii Nav’e-Stoksa [Mathematical Modeling of Convective Heat Transfer on the Basis of the Navier-Stokes Equations]. Moscow: Nauka, 1987, 274 p. (In Russ.)

18. Tatarchenko V.A. Ustoichivyi rost kristallov [Sustainable crystal growth]. Moscow: Nauka, 1988, 240 p. (In Russ.)

19. Babsky V.G., Kopachevsky N.D., Myshkis A.D., Slobozhanin L.A., Tyuptsov A.D. Gidromekhanika nevesomosti [Hydromechanics of zero gravity]. Moscow: Nauka, 1976, 504 p. (In Russ.)

20. Prostomolotov A.I., Verezub N.A., Mezhennii M.V., Reznik V.Ya. Thermal optimization of CZ bulk growth and wafer annealing for crystalline dislocation-free silicon. J. Cryst. Growth. 2011; 318(1): 187–192. https://doi.org/10.1016/j.jcrysgro.2010.11.080

21. Prostomolotov A.I. Coupled modeling of thermal processes and defect formation during the growth of dislocation-free silicon single crystals. Izvestiya Vysshikh Uchebnykh Zavedenii. Materialy Elektronnoi Tekhniki = Materials of Electronics Engineering. 2009; (2): 66–77. (In Russ.)

22. Gaponov V.A. Kompaktnye raznostnye approksimatsii povyshennogo poryadka tochnosti v zadachakh vychislitel’noi gidrodinamiki [Compact difference approximations of higher accuracy in computational fluid dynamics problems]. Novosibirsk, 1994. (Prep. RAS. Institute of Thermophysics; No. 272-94). (In Russ.)

﻿Abstract