Review Article 
Corresponding author: Vladimir S. Berdnikov ( berdnikov@itp.nsc.ru ) © 2021 Vladimir S. Berdnikov.
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): 151165. https://doi.org/10.3897/j.moem.7.4.51073

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 heatinduced 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
crystal growth, Czochralski technique, melt hydrodynamics, heat exchange, thermal gravity capillary convection, numerical simulation.
Below is a continuation of the overview of experimental and numerical results for crystal pulling from melt using the Czochralski technique [
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 selfconsistent 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, concentrationrelated 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 concentrationcapillary one. Analysis of joint heat and mass transport processes dramatically complicates the initial problem. The emergence thresholds and types of flow instabilities [
The ideal physicalmathematical 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 T, σ = σ_{0}(1 – β_{σ}ΔT), where β_{σ} = (1/β_{0})(–∂β/∂T); and β_{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 [
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 halfstep were taken from the previous halfstep. The equations were reduced to sets of equations with threediagonal matrices and solved by passing. The optimum grid parameters were determined preliminarily [
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 = (βg/ν^{2})ΔTR_{sc}^{3} or the Rayleigh number Ra = GrPr = (βg/αν)ΔTR_{sc}^{3} 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_{0}) /ΔT, T is the temperature, T_{0} 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.
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.
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^{1}; (b, f) 10^{2}; (c, g) 10^{3}; (d, h) 10^{4}.
Along with the isotherms, additional temperature field characteristics are temperature profiles in horizontal and vertical sections which are measured in the physical experiment [
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.
(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.
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^{2}; (4) 10^{3}; (5) 10^{4}; (6) 78274.
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.
Temperature profiles are used for calculating the radial distributions of the local heat flow at the CF (Fig.
Calculations for Pr = 0.05, 50 and 2700 were carried out for the same geometry (Figs
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: (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^{4}; (b, f) Pr = 16, Gr = 35.25; (c, g) Pr = 50, Gr = 10; (d, h) Pr = 2700, Gr = 0.185.
(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}; (2) 4.7 ∙ 10^{3}; (3) 4.7 ∙ 10^{4}; (4) 4.7 ∙ 10^{5}.
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}; (2) 4.7 ∙ 10^{3}; (3) 4.7 ∙ 10^{4}; (4) 4.7 ∙ 10^{5}.
Temperature distribution in layer height for TGC mode in subcrystal region at Pr = 0.05, H/R_{c} = 0.7, R_{c}/R_{sc} = 2.76 and Gr = 2.3 ∙ 10^{6}: (1) r = 0; (2) 0.25; (3) 0.5; (4) 0.75; (5) 1.
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^{2}; (3) 4.7 ∙ 10^{3}; (4) 4.7 ∙ 10^{4}; (5) 4.7 ∙ 10^{5}.
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^{5}; (d, e) 2.3 ∙ 10^{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.
Temperature distribution in layer height for TGC mode in subcrystal region at Pr = 0.05, H/R_{c} = 0.25, R_{c}/R_{sc} = 2.76 and Gr = 2.3 ∙ 10^{6}: (1) r = 0; (2) 0.25; (3) 0.5; (4) 0.75; (5) 1.0.
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.
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.
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^{3}; (c, g) 3 ∙ 10^{3}; (d, h) 10^{4}.
Temperature distribution in layer height for TGCC mode in subcrystal 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.
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^{3}; (5) 10^{4}; (6) 10^{5}; (b): Ma = 34240; (1) Gr = 0; (2) 10^{3}; (3) 10^{4}; (4) 41272.
TCC is theoretical zerogravity 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.
– circulation over the external contour and integration into the main superficial swirl having a thermocapillary origin;
– emergence of a bottomreflected 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.
(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
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 lowintensity 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 subcrystal 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 TCCcontrolled 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 ААААА171170228500213.