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

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.

Below is a continuation of the overview of experimental and numerical results for crystal pulling from melt using the Czochralski technique [^{2}. One fine point is that reference books on heat exchange provide typically heat exchange process parameters as a function of Gr/Re^{2} 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 [

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.

Here s is the melt surface tension coefficient depending on the temperature _{0}(1 – β_{σ}Δ_{σ} = (1/β_{0})(–∂β/∂_{0} 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 [_{0} and the crucible side wall _{w} temperatures are set to be different but constant. Therefore the characteristic temperature difference in the system is Δ_{w} – _{0}. 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 _{ij}|/max|ω_{ij}^{-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 [

Laminar convection of different origins was studied: thermal gravity, thermocapillary and thermal gravity capillary [

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 = (β^{2})Δ_{sc}^{3} or the Rayleigh number Ra = GrPr = (β_{sc}^{3} which determine TGC intensity, the Marangoni number Ma = [(–∂β/∂_{sc} which describes the contribution of the thermocapillary effect and the Prandtl number Pr = ν/_{0}) /Δ_{0} is the temperature at the CF, Δ_{P}_{P}_{sc} and _{c} are the crystal and crucible radii, respectively. The following scales were used for the transition to dimensionless equations: _{sc} for length, Δ_{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:

Simulation area diagram for ideal physical-mathematical crystal growth model: (

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. _{sc} = 0.7 and _{c}/_{sc} = 2.76. Figure

Isotherm field in heat conductivity mode (_{c}/_{sc} = 2.76).

Evolution of (_{c} = 0.7, _{c}/_{sc} = 2.76): ^{1}; (^{2}; (^{3}; (^{4}.

Along with the isotherms, additional temperature field characteristics are temperature profiles in horizontal and vertical sections which are measured in the physical experiment [_{c} + _{sc})/2 section at different Gr is presented in Fig. _{sc} ≤ _{c} region at the top flow branch is that the top boundary is free and without friction (see Fig. _{c} + _{sc})/2), the velocity being the highest there. For 0 ≤ _{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. _{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.

Temperature distribution in layer height in the _{c} + _{sc})/2 section for TGC mode at Pr = 16, _{c} = 0.7 and _{c}/_{sc} = 2.76: (

(_{c} = 0.7 and _{c}/_{sc} = 2.76 for different Gr: (

Radial distributions of local heat flows in TGC modes at Pr = 16, _{c} = 0.7 and _{c}/_{sc} = 2.76: (^{2}; (^{3}; (^{4}; (

TGC structure dependence: (_{c} = 0.7 and _{c}/_{sc} = 2.76: (

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. ^{4}). 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 [

Temperature profiles are used for calculating the radial distributions of the local heat flow at the CF (Fig. ^{4} isotherms in Fig.

Calculations for Pr = 0.05, 50 and 2700 were carried out for the same geometry (Figs ^{3}. Melt core isotherms for 2 ∙ 10^{3} ≤ Gr ≤ 10^{4} acquire horizontal sections, this regularity being traced in radial temperature distributions (see Fig. ^{3} are almost coincident with that distribution for heat conductivity mode (see Fig. ^{3}. Upon further Gr growth the ψ and

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

TGC structure dependence: (_{c} = 0.7 and _{c}/_{sc} = 2.76: (^{4}; (

(_{c} = 0.7 and _{c}/_{sc} = 2.76: (^{2}; (^{3}; (^{4}; (^{5}.

Temperature distribution in layer height in (_{c} + _{sc})/2 section and (_{c} = 0.7 and _{c}/_{sc} = 2.76: (^{2}; (^{3}; (^{4}; (^{5}.

Temperature distribution in layer height for TGC mode in sub-crystal region at Pr = 0.05, _{c} = 0.7, _{c}/_{sc} = 2.76 and Gr = 2.3 ∙ 10^{6}: (

Radial distribution of local heat flows in TGC mode at Pr = 0.05, _{c} = 0.7 and _{c}/_{sc} = 2.76 for different Gr: (^{2}; (^{3}; (^{4}; (^{5}.

Evolution of (_{c} = 0.25 and _{c}/_{sc} = 2.76:(^{5}; (^{6}.

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. _{sc}/λ = (∂_{ct}_{sc}/Δ_{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 _{ct} = αΔ^{4}, 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.009}Gr^{0.2} in the 0 ≤ Gr ≤ 2500 range. The data shown in Figs ^{4}.

Temperature distribution in layer height for TGC mode in sub-crystal region at Pr = 0.05, _{c} = 0.25, _{c}/_{sc} = 2.76 and Gr = 2.3 ∙ 10^{6}: (

Dimensionless heat emission coefficient as a function of Rayleigh number: (

Figures

First of all, TCC was studied at Pr = 16 in a system having the same geometry as the one for the TGC study (Fig. ^{5} [

Evolution of (_{c} = 0.7, _{c}/_{sc} = 2.76 and Pr = 16: (^{3}; (^{3}; (^{4}.

Temperature distribution in layer height for TGCC mode in sub-crystal region at Pr = 16, _{c} = 0.7, _{c}/_{sc} = 2.76, Gr = 41270 and Ma = 34240 with equal radial intervals (

Radial distribution of local heat flows in (_{c} = 0.7 and _{c}/_{sc} = 2.76: (^{3}; (^{4}; (^{5}; (^{3}; (^{4}; (

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 _{c} + _{sc})/2,

– 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.

The flow pattern shown in Fig.

(_{c} = 0.7, _{c}/_{sc} = 2.76, Gr = 5000 and Ma = 34240.

Flow regularities affect the radial heat flow distribution at the CF. Figure _{sc}/2 range (see Figs

Studies of the temperature distribution regularities in the volume and along the free fluid surface in TCC, TGC and then TGCC modes [

Thermal gravity, thermocapillary and gravity capillary modes of laminar convection were numerically studied. The calculations were carried out for the 0 ≤ Gr ≤ 10^{7} Grashof number range and the 0 ≤ Ma ≤ 10^{5} 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.

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 [

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 [

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