A Single-domain Approach to Simulate the Effect of Convective Flow on the Mushy-zone Structure in Czochralski Growth of Gadolinium Gallium Garnet Crystal

Abstaract A numerical study was carried out to describe the effect of the melt hydrodynamics on the crystallization front shape in the Czochralski growth of a semitransparent oxide crystal. In the present model calculation, the enthalpy-porosity method was used to solve the phase change problem with the convection due to buoyant, thermocapillary and centrifugal forces. It was shown that the rotationally-driven flow protrudes into the mushy zone when the crystal rotation rate was increased to a certain critical value corresponding to Gr/Re=0.89 as the ratio between the intensity of buoyancy and forced convection flow in the melt. The ratio between the vertical and horizontal temperature gradients beneath the mushy zone was found to be decreased by increasing the crystal rotation rate. It was shown that the shape of the zone deforms abruptly when the ratio between the axial and radial temperature gradients decreased to the values smaller than the unity. The Burger’s number condition was found to be violated in the case Gr/Re<0.89, at which the onset of geostrophic instability is expected.


Introduction
Oxide single crystals are conventionally grown by the Czochralski (Cz) method primarily characterized by hydrodynamics of the melt which is inextricably coupled to transport phenomena in the melt. The Cz melt is in fact a good example of situation in which the fluid motion is brought about by different coexisting mechanisms, that is, heating from side (the crucible wall) yields a convective flow of the Hadley type, whereas cooling from above (the melt/gas interface and the crystal/melt boundary) may lead a Rayleigh-Benard type convection in the system [1]. Motions relevant to the Cz melt can be classified by the principle of driving forces into the following main groups: (a) gravitational (natural convection), (b) mechanical (forced convection), and (c) surface-tension (Marangoni convection). These different kinds of the melt motion determine the structure of the flow, the heat and mass transport in the fluid, the crystal-melt interface shape, and consequently the quality of the growing crystal.
In a Cz melt model, the ratio between the buoyancy and the rotationally-driven forces, denoted by Gr/Re 2 , is generally assumed as the non-dimensional parameter governing the convective interactions, and the flow is supposed to become unstable whenever its value approaches to a critical one around the unity. It has been shown that refractory oxides may exhibit an abrupt change in the interface shape, from a convex-to-melt phase boundary to a flat or even concave one, during the growth process [2][3][4][5][6]. This sudden and uncontrollable change is usually assigned to a critical value of the ratio Gr/Re 2 in the melt. However, the manner in which the crystal-melt interface is represented is a central feature of bulk crystal growth models in which the phase change problems with convection play a significant role. A self-consistent growth model requires that the interface geometry be computed as part of the solution to the transport problem, and the main challenge is the presence of a moving phase boundary involving a strong coupling of mass and heat transfer. From an algorithmic point of view, two methods are primarily employed for solving phase change problems with convection and computing the position of the crystal-melt interface. The first one (front-tracking methods) is based on a multi-domain approach and defines a discrete moving surface to separate 74 A Single-domain Approach to Simulate the Effect of Convective Flow on the Mushy-zone Structure in Czochralski Growth of Gadolinium Gallium Garnet Crystal the interface between crystal and melt. The second category (diffuse-interface methods) is based on a single-domain approach where a system of momentum and energy equations is solved in the entire physical domain, and treats the interface as a region (mushy zone) of finite thickness across which physical properties vary rapidly but continuously from one bulk value to the other [7]. In contrast to the front-tracking methods, such as the widely used isotherm method, in a single-domain approach the thickness of the so-called mushy zone is very large compared to the width of a phase boundary of one or two layers of computational cells [8]. The main advantage of this approach is that the interface is not explicitly computed and the energy balance condition is automatically satisfied at the crystal-melt interface. Voller and Prakash [9] have argued that when the material is not a pure substance, the phase transition occurs over a temperature range around the melting point. That is, the evolution of latent heat has a functional relationship with temperature as opposed to the step change associated with an isothermal phase change. In such cases, the phase transition takes place through the so-called mushy region, in which the fluid flow plays an important role in the final solidified microstructure. This numerical approach to solve the convection-diffusion controlled mushy region phase change problems is called the enthalpy-porosity model which includes latent heat effects as a source term, and comprises a technique to ensure that the velocity field vanishes in the solid phase [10,11]. In the present study, the enthalpy-porosity model was used to describe the phase change problem involving convection and radiation in a Cz oxide crystal growth configuration. In this investigation, the finite volume approximation (FVM) is applied to discrete the governing equations on stationary grid mesh system incorporating enthalpy-porosity model to simulate the effect of rotationally-driven forces on the mushy region structure. Correlation between hydrodynamic stability of the flow and the mushy-zone structure is studied. We intended to show that the mushy-zone morphology would be considerably deformed at a critical value of the ratio between the buoyancy-and rotationally-driven forces.

Geometry, Thermal and Radiative Properties
The schematic in Figure 1 illustrate the geometry of the Cz growth configuration adopted in the present numerical simulation, where the crystal radius at the constant-diameter stage of the growth is r x = 30mm. The crystal to crucible radii ratio (r x /r c ) is 0.5, and the aspect ratio AR=h c /r c is 2.0, where h c =120mm is the height of the crucible side wall. The height of the tri-junction point J, is kept at 125mm, while the height of the melt free surface at the crucible wall is h l = 120mm. The crucible side wall is at a constant temperature T c =T mp +∆T max , where ∆T max = 67K is the applied temperature difference between the crucible wall and the melting point of the material. The crucible wall is assumed to be an opaque and gray surface diffusely emitting and reflecting but not absorbing the thermal radiation. The crucible bottom is thermally insulated, and its surface has the same radiation transfer properties ( =0.5) as the crucible side wall. Both the melt and the crystal have the same optical thickness, and they absorb and emit but not scatter radiation. The optical properties of the material (such as refractive index, n=1.8 and absorption coefficient, a=258m -1 [6]) are independent of temperature and wavelength. Throughout this work, the melt-gas interface is assumed to be a semitransparent diffuse gray surface, which is slightly curved to form a meniscus configuration at the crystal periphery. With the refractive index of the melt (and the crystal), the Spuckler-Siegel [12] approximation was used to estimate the transmissivites (τ ext =0.85, and τ in =0.23) on both sides of the melt free surface.
The emissivity, at the crystal surface is assumed to be 0.3. Energy is transmitted through the free surface of the melt (and the crystal) to the ambient walls ( =0.7) at T a ≅1650K. The ambient gas is assumed to be totally transparent. The crucible is at rest, and the no-slip condition is applied to all physical boundaries of the melt except the melt-gas interface on which thermocapillary forces are taken into account. The crystal rotates uniformly around the symmetry axis. The crystal rotation rate, 0.5<Ω (rad/s)<3.0, is the only variable external parameter in the present model. In contrast to the front-tracking method where the interface shape coincides with the melting point isotherm T mp =2023K, the crystal-melt mushy zone is treated as a region of finite thickness across which the phase change occurs. Figure 1. Geometry of the Cz growth model with a hypothetical mushy region between the crystal and melt. Dashed horizontal line, L r (5 mm below the melt free surface) and vertical line L z , crossing the point M(r,z)= (r x /2, 125mm ) are described. The point J is the tri-junction point, h m is the height of the melt meniscus region, and ∆ℎ is the convex to melt interface depth.

Dimensionless Parameters
The thermophysical properties of the material (gadolinium gallium garnet) employed for the present calculations, and the dimensionless parameters which describe the convective flow in the melt are listed in Table 1. The oxide (Gd 3 Ga 5 O 12 ) melt is characterized by the Prandtl number Pr≅4.69, the Grashof number Gr=1.141×10 3 ∆T max , the Reynolds number Re=1.271×10 2 Ω, and the Marangoni number Ma= 2.742×10 2 ∆T max , where Ω(rad/s) is the crystal rotation rate and ∆T max =67K is the driving temperature difference which ensures the growth of crystal at the tri-junction point, J(r,z)=(30mm, z=125mm). For the present mixed convection problem, the ratio between buoyancy and the rotationally-driven forces is given by Gr/Re 2 =4.732/Ω 2 , and the critical Reynolds number is approximated by Re * ≅ Gr 1/2 , at which (Re * ≅276.5) the intensity of forced convection flow is comparable to that of the buoyancy-driven one. In a Cz melt model, and particularly in the column of the flow beneath the phase boundary, both density variations (thermal stratification) and the crystal rotation rate are intrinsic to the hydrodynamic stability of flow. As a measure of the relative importance of thermal stratification and the rotationally-driven effects, the Burger number, Bu is used to predict the breakdown of axisymmetry with the onset of the geostrophic instability if the condition Bu < Bu * = 0.58 is fulfilled. In the present Cz model (AR=2.0), the Burger number can be expressed by Bu=(N/Ω) 2 , where N=( / ) 1/2 is the buoyancy frequency with which an adiabatically displaced fluid particle will oscillate about its equilibrium position [1]. The condition (Bu<Bu * ) is satisfied whenever the crystal rotation rate exceeds a critical value at which the convective flow is characterized by Gr/Re 2 ≤0.89. The most direct measure of the relative importance of the density and velocity gradients is provided by the gradient form of the Richardson number, Ri g , which can be expressed as Ri g =4Bu/(Ro T ) 2 [13], where Ro T is the thermal Rossby number. Fein and Pfeffer [13] have reported that, for all the experimental cases studied (Pr<1, and Pr>1), Ri g increases as the rotation rate is increasing.  Thermocapillary coefficient ( Thermal conductivity (

Mathematical Model
Assuming that the Cz/oxide melt is Newtonian and incompressible, the enthalpy-porosity governing equations [9][10][11] can be expressed as  Conservation of momentum:  Conservation of mass:  Conservation of energy: where (x,t) is the flow velocity vector at any point x(r,z) and at any time t, is the pressure, is the enthalpy, is the radiation heat flux through the melt. is defined so that the momentum equations are forced to mimic the Carman-Kozeny equation [10] which suggests the following form for the function in the momentum equation, where the factor is a small constant value introduced to avoid division by zero, is a constant accounting for the mushy-zone morphology, and is the porosity which takes the values, = 1 in the liquid phase, = 0 in the solid phase, and 0 < < 1 in the mushy region. The equation (5) implies that, vanishes in the fluid and has no influence in the liquid phase. In the solid phase, however, takes a large value and forces the velocity field to be vanished. Assuming the Boussinesq approximation to be valid, that is, density is constant in all terms except a gravity source term, that buoyancy source term can be expressed as where = (0,0, − ) is the acceleration vector due to gravity, is the thermal expansion coefficient, is the specific heat of the material, is a reference value for the enthalpy, and the enthalpy of the material (the total heat constant) is given by where ∆ is the latent heat of fusion.
To estimate the radiative heat flux, in Eq. (3), we have to solve the radiative transfer equation. In the present study, the S n -type (S4) discrete ordinates (DO) method [14], based on a discrete representation of the angular dependence of the radiation intensity, was used. The balance of energy passing in a specified direction through a small differential volume in an axisymmetric system is given by the following equation when the scattering is ignored: where , and are the direction cosines and is the intensity of radiation for the discrete direction . is the angle of revolution around the z-axis, and is the intensity of blackbody radiation at the temperature of medium. The absorption coefficient of the melt (and the crystal), is denoted by a. The radiative heat flux, can be expressed as follows: where is angular quadrature weight which sums to the surface area of unit sphere, and the index, m is used to sum over all directions.
We applied the finite volume technique (FVM) to compute steady and axisymmetric solutions to the fully coupled equations governing the convective flow involving solidification in the Cz growth system represented in Figure  (1). The SIMPLE algorithm was used to couple velocities and pressure on staggered grids, and second order upwind method was used for discretization of momentum and energy equations. The equations are integrated over each control volume (CV) and resulting system of algebraic equations was solved iteratively until convergence was reached. Non-uniform finite volume mesh structure used to perform the calculations in this work, consists of 36446 mixed (quadrilateral and triangular) cells. According to the investigations of Voller and Prakash [9], the constants b and C in the equations (4) and (5) are set at 10 -3 and 1.6×10 5 , respectively. The numerical results in the next section implies that the chosen value of C is small enough to allow for significant flow in the mushy region whereas the limiting value of the function , that is, = − / is large enough to suppress the flow velocity field in the solid phase.

Results and discussion
Steady-state simulations of the convective flow pattern and temperature field are represented in Fig. 2(a,b,c). These figures are relevant to three different crystal rotation rates, i.e., Ω= 1.5, 2.0, and 2.5 rad/s corresponding to Gr/Re 2 = 2.10, 1.18, and 0.75, respectively. The general feature of the fluid motion is described as following. The buoyancy-driven hot flow ascends along the crucible side wall and then, accompanied by the thermocapillary flow, travels along the melt free surface towards the crystal rim. Due to the convection and radiation at the melt free surface, the fluid is being cooled down along the path, and this creates a stream of cold fluid which descends along the melt centerline towards the crucible bottom. Due to the optical properties of the melt, the flow does not exhibit undulating structure observed in high Prandtl number opaque melts [3]. However, owing to the relatively large optical thickness of the fluid, the isotherms are more squeezed at the corner, and the vortex center of thermal convection flow is located at a position lower than the melt middle. When the rotation rate is smaller than ~2.3 rad/s, the flow is essentially a buoyancy-driven one with a large, unicellular counter-clockwise circulation. Increasing in the rotation rate to ~2.5 rad/s, small rotationally-driven, clockwise cells appear in the vicinity of the crystallization front and nearby the crucible bottom. As shown in Figure 2(a, b, c), in the cases with (Gr/Re 2 )>1.18, the flow pattern is associated by a convex-to-melt interface, which tends to be nearly flat when the rotation rate is increased to 2.5 rad/s, corresponding to Gr/Re 2 =0.75. In Figure 2(a), the isotherms close to the phase boundary are regular and more density of heat flux lines implies that the heat transfer is intensified along the interface. In this case the isotherms T(1)=2056.5K and T(2)=2060K both are found to be reached the crucible bottom. Figure 2(b) shows that increasing the rotation rate from 1.5 to 2.0, the regular shape of the isotherms in the vicinity of the interface is changed at around the interface mid-point, r=15mm, and more density of lines are located between meniscus region of the melt and the convex-to-melt interface mid-point. In this case, the isotherm T(1) turns toward the axis while T(2) reaches to the crucible bottom. Further increase in the rotation rate (from 2.0 to 2.5 rad/s) results in a nearly flat interface. As shown in Figure 2(c), the isotherms are more squeezed in the meniscus region, and they are strongly deformed along the rest of the interface. In this situation both of the isotherms T(1) and T(2) are found to be turned and reached to the symmetry axis. These results are consistent with the well-known [15] effect of the melt meniscus at the crystal periphery which allows for radial heat transfer (both convective and radiative) from the melt to surrounding near the crystallization front and reduces the net radial heat transfer from the melt to the interface. This effect leads to decrease in the convexity of the interface. The ratio between the axial and radial temperature gradients decreases with increasing rotation rate.
The effect of rotationally-driven forces on the morphology of mushy-zone is represented in Figure 3. Note that the ratio between the axial and horizontal temperature gradients, g z /g r decreases with increasing the crystal rotation rate, and the mushy region was found to be strongly deformed when (g z /g r ) beneath the phase boundary decreased to the values smaller than unity (i.e., in the cases with Ω ≥ 2.3rad/s). Figure 3 shows that the convex-to-melt interface at Ω=1.5 rad/s is changed to form a "gull-wing" geometry at Ω=2.0 rad/s. The effect of rotation rate on the radial velocity profile along the horizontal line, L r is illustrated in Figure 4. The line L r (5mm below the melt free surface) is crossing the mushy region (0<r(mm)<30). Figure 4 shows that the convective flow protrudes into the zone, and the velocity profile was found to be changed for the cases Ω ≥ 2.0 rad/s. In other words, when Ω ≤ 2.0 rad/s, the fluid motion within the mushy-zone is dominated by the buoyancy and thermocapillary forces (u r <0). The sign of the velocity component, which was altered at the rotation rates even slightly greater than 2.0 rad/s (e.g., Ω = 2.3 rad/s, corresponding to Gr/Re 2 =0.89), implies that the forced convection flow plays a significant role in deformation of the mushy-zone morphology. Correspondingly, the effect of rotation rate on the temperature profile along the same horizontal line (L r ) is depicted in Figure 5. It is shown that the crystallization front (0<r(mm)<30) tends to be flat at Ω=2.5 rad/s, and this can be inferred that the interface-inversion process occurs at Ω=2.8-3.0 rad/s. It is worthy to note (but not shown here) that, the flow field found to be dominated by the rotationally-driven forces when the rotation rate increased from 2.5 to 2.8 rad/s at which a large cold plume appears beneath the mushy-zone and, descending along the symmetry axis toward the crucible bottom, tends to break away from the axial position and the flow direction. The effect of rotation rate on the interface depth ( ∆ h), the thickness of the mushy region ( ∆ L z ) at r=r x /2, and the averaged axial and radial temperature gradients in the melt beneath the interface is illustrated in Figure 6. It is shown that the g r increases, and g z decreases first (0.5<Ω(rad/s)<2.0) slowly and then (Ω>2.0 rad/s) more rapidly with increasing rotation rate. The depth of the crystallization front was found to be sharply decreased for the cases Ω>2.0 rad/s (Re > 254.24 and Gr/Re 2 <1.18) which corresponds to (g z /g r )≤1.0. The morphology of mushy-zone depends on the rotationally-driven forces which affects g z /g r in the Cz melt. As expected [13], Figure 7 shows that the gradient Richardson number (Ri g ) increases and the thermal Rossby number (Ro T ) decreases as the forced convection flow, represented by the Reynods number (Re), is intensifying. It is shown that, increasing the rotation rate from 2.0 to 2.3 rad/s the governing parameter of mixed convection (Gr/Re 2 ) decreases from 1.18 to a value lower than unity (i.e., 0.89). Consistently, Ro T decreases from 1.59 to 0.81 in the case with Ω=2.3 rad/s at which the Burger number is equal to its critical value (Bu * =0.58) [1,13], and the morphology of the mushy region was found to be clearly modified.

Conclusions
The enthalpy-porosity method developed by Voller and Prakash [9,10] was used to solve the phase-change problem involving mixed convection and radiation in a Czochralski crystal growth model. A finite volume (FVM) method was applied to compute two-dimensional and axisymmetric solutions to the fully coupled equations governing the melt hydrodynamic and heat transfer in the model, and discrete ordinates method (DOM) was used to account for the influence of internal radiative heat transfer on the flow field and the temperature distribution in the material. The effect of rotationally-driven forces on the structure of the mushy region was studied. It was shown that the morphology of mushy-zone depends on the ratio between the axial and radial temperature gradients (g z /g r ), which decreases as the crystal rotation rate increased. The mushy-zone structure was found to be correlated with the hydrodynamic stability of the flow.