Numerical simulation of crystallization on a stationary and rotating cooled disk

The process of crystallization of silicon melt on a monotonically cooled disk located on the free surface of the melt is studied numerically. The influence of thermogravitational, thermal gravitational-capillary and mixed convection on the shape of the crystal-melt interface is investigated. In mixed convection modes, the speed of uniform rotation of the disk is set. The calculations were carried out in an axisymmetric formulation of problems by the finite element method using an adaptive triangular mesh and taking into account the latent heat.


Introduction
Crystallization from melts occurs in systems that are inevitably non-isothermal and in a gravity field. Therefore, thermogravitational convection develops in melts [1][2][3]. In the presence of a free nonisothermal surface of the melt, the action of the thermocapillary effect (TCE) is added and convection has a thermal gravitational-capillary nature [1][2][3]. The spatial shape and intensity of flows depend on the characteristic temperature drops, geometry and absolute dimensions of the working section of the growth apparatus. In the Czochralski and Kyropoulos methods, crystals are pulled from the free surface of the melt layers [1][2][3]. In works dealing with the analysis of the processes of obtaining single crystals by directed crystallization methods, it is emphasized that free convection is poorly controlled. The Czochralski method uses crystal rotation to control the hydrodynamics of melts. In this case, the physical nature of convective motion is gravitational-centrifugal [2][3][4]. Despite a large number of studies of directed crystallization methods, due to the nonlinearity of the interaction of various flow generation mechanisms, there are no answers to many questions. Until now, there are few works in which convective heat transfer is studied for cases when the forces of buoyancy, TCE, and rotation effects act [4]. An important, but insufficiently studied stage is the initial stage of crystal growth. Technological practice shows that by choosing the mode of convective heat transfer, it is possible to correct the shape of the crystal-melt interface (CMI). The development of technologies for growing high-quality single crystals requires an understanding of the features of conjugate convective heat transfer for various combinations of thermal conductivity of melts, crystals, and crucible materials. This work is aimed at investigating this influence. For the silicon-graphite system, these studies are the development of works [2,3]. Direct experimental investigations of high-temperature technological processes are expensive and laborious. It is practically impossible to measure the characteristics of non-stationary temperature fields in the composite region of the crucible-melt-crystal of the growth apparatus. Therefore, it is advisable to numerically investigate the effect of conjugate convective heat transfer on the crystallization process and on the forms of CMI.

Statement of the problem
The problem statement corresponds to a simplified model of the initial stages of crystal growth in the Kyropoulos and Czochralski methods. In these methods, crystals are grown by pulling them out of the melt layer with a free surface [1]. Crystallization of a silicon melt in a graphite crucible is investigated. The crucible is a cylindrical container made of MPG-6 graphite with a wall and bottom thickness of 0.5 cm. The inner radius of the crucible is 5 cm. The height of the melt layer is 5 cm. The crystallization process begins on a monotonically cooled disk located at the upper boundary of the melt. The radius of the disc is 2 cm. The problems were solved in an axisymmetric setting in a composite computational domain, the diagram of which is shown in figure 1. Only the right side of the axisymmetric computational domain is shown here. At the initial moment of time t = 0, the crucible (Ω 3 ) is filled with melt (Ω 1 ), overheated relative to the crystallization temperature (1410 °C) by 10 K. The initial melt temperature is maintained on the outer vertical surface of crucible S 7 . The disk is cooled by lowering the temperature at the S 10 boundary at a rate of 2 K/min. The crystallization process begins when the temperature on the disk drops to the crystallization temperature. In the mixed convection mode, the angular velocity of uniform rotation of the disk is set (i.e., boundaries S 10 or S 9 and S 11 ). The upper free boundary of the melt (S 3 ) is set flat and non-deformable; either the frictionless condition is specified on it, or TCE is taken into account. The total volume of the crystallized substance (Ω 2 ) and the melt is assumed to be constant. Convective heat transfer in a melt is described by a system of dimensionless equations of unsteady thermogravitational convection in the Boussinesq approximation in the variables temperature, stream function, vortex of velocity and in a cylindrical coordinate system have the form: U,V,W are the components of the velocity vector, T is temperature, ω is a vortex of velocity, ψ is the stream function, The dimensionless equations of heat conduction in the crystal and crucible have the following form, respectively: 2 2 , .
In the equations being solved c f , c s , c c are the specific heat coefficients of the melt, crystal, and crucible, respectively, λ f , λ s , λ c are the thermal conductivity coefficients, ρ f , ρ s , ρ c are the densities, ν f is the kinematic viscosity of the melt, β f is the volumetric thermal expansion coefficient of the melt, g is acceleration due to gravity, R T is scale of geometric dimensions. The latent heat is accounted for through an internal heat source. Let for a time instant equal to с the  crystal-melt interface changes its position from B 0 (t) to B 1 (t) (figure 2), respectively, the source power per unit volume in the domain Ω 0 (t), bounded by the surfaces B 0 (t) and B 1 (1) and (2) the value of Q is defined as: is not known in advance. It is determined by solving discrete analogs of equations, by their multiple solution within one time step with a relaxation coefficient, while the position of the boundary B 1 (t), corresponding to the isothermal surface with a temperature value equal to the crystallization temperature is redefined, and the position of the boundary B 0 (t) remains unchanged.
At the boundaries between the «crucible-melt», «crucible-crystal», «melt-crystal», the conditions of ideal thermal contact are fulfilled, i.e., at the indicated boundaries, the temperature and heat flux are inseparable: , The following conditions are set on the axis of symmetry:  At the initial moment of time, the temperature in the entire system is constant and exceeds the crystallization temperature, and there is no convective flow: The initial temperature of the system is maintained on the outer vertical side surface of the crucible: 7 * .

S T T 
The temperature on the disk changes over time: where dT dt is the specified disk cooling rate. Azimuthal velocity at the disk and crystal surfaces (S 9 ,S 10 ) is a function of the radius: 9,10 * , W is the dimensionless azimuthal velocity at r = 1. The problems were solved by the finite element method [5]. To approximate the solution, linear functions on triangles were used. Mesh generation is based on the direct method of constructing the Delaunay triangulation over the largest angle with cellular acceleration. The grid is built for each time step, so as to track the shape of the CMI at the current step and to track the shape of the CMI at the next step with automatic determination of the volume of the crystallized substance. In addition, the grids are thickened to the CMI on both sides and to the boundaries of the calculated area.
The calculations were carried out at the values of thermophysical parameters characteristic of silicon crystallization in graphite crucibles [6][7][8]

Results and discussion
In the modes of nonstationary thermogravitational convection, the crystallization of silicon melt on the surface of a monotonically cooled disk after reaching the crystallization temperature. Then, similar studies were carried out taking into account the TCE in the mode of thermal gravitational-capillary convection. So the influence of unsteady natural convection on the hydrodynamics of the melt, heat transfer, and the shape of the CMI, excited either only by buoyancy forces or by the combined influence of buoyancy forces and the TCE, has been studied. At the next stage, these modes of natural convection were the initial ones in the calculations of mixed convection. In mixed convection modes, a forced flow occurs when the cooled disk rotates uniformly. Figures 3a and 4a show the combined fields of isotherms and isolines of the stream function in the modes of thermogravitational and thermal gravitational-capillary convection in the «silicon-graphite» system at the same times t = 750 s. Figures 3b and 4b show the combined fields of isotherms and isolines of the stream function in mixed convection modes at an angular rotation rate of the cooled disk ω К = 20 rpm at the same time moments t = 750 s. The step between the isolines of the stream function ψ h in figures 3a and 4a is 10 mm 2 /s, in figures 3b and 4b is 5 mm 2 /s. Isolines of the stream function with a value of 0 mm 2 /s are excluded. Figures 3, 4 show only the right-hand sides of the axisymmetric fields of the stream function isolines and isotherms. The directions and intensity of the flow can be understood from the data on the profiles of the radial and axial velocity components in    is directed from the crucible wall to the disk edge, and then to the CMI edge ( figure 7). In the mode of thermal gravitational-capillary convection, the speed along the free surface is significantly higher ( figure 7). This is also seen in figures 5 and 6 from a comparison of the velocity values on the free surface of the melt. A feature of the near-surface sections of the profiles in figures 5, 6 is the absence of a velocity gradient in the thermogravitational convection mode and its presence in the thermal gravitational-capillary convection. This is the influence of the TCE in the presence of a longitudinal temperature gradient along the free surface of the melt (figure 9).  A significant effect of TCE is also observed in mixed convection modes. This is especially clearly seen from the velocity profiles along the free surface of the melt in figure 7. Disregarding the TCE, the flow along the free surface is directed from the rotating surface to the crucible walls. In this case, a forced centrifugal flow from under the rotating surface displaces or even suppresses the free convective flow of the heated melt from the crucible walls to the CMI edge.
A clockwise vortex forms under the free surface. This can be seen in figures 5, 6 (curves 4). When TCE is taken into account along the free surface, the flow direction in the mixed convection regimes remains the same as in the thermal gravitational-capillary convection regime.