The hysteresis phenomenon of flow patterns due to thermal and solutal Marangonis convection in a liquid bridge under zero gravity was numerically simulated. A three-dimensional half-zone configuration of the Floating Zone (FZ), formed between hot and cold discs, was considered. The field equations; continuity, momentum, energy, and mass transport equations, were solved numerically by the PISO algorithm, and the open source software OpenFOAM was used. Three sets of initial conditions were used. The first is the conditions of no-flow and constant solutal Marangoni numbers, the other two are the gradually increasing or decreasing solutal Marangoni numbers. Results showed that the flow becomes 3-D unstable at a larger solutal Marangoni number than the critical number (Ma_C)_cr, and the (Ma_C)_cr depends on the initial condition. The flow structure clearly shows hysteresis behavior of the flow field with about 24% difference between the critical values in the hysteresis diagram.

Keywords: numerical simulation, marangoni convection, floating zone technique, half-zone, openfoam

FZ, floating zone; HZ, half zone; FVM, finite volume method; PISO, pressure-implicit split-operator; Pr, prandtl number; Sc, schmidt number

The floating zone (FZ) method is one of the growth techniques used to grow single crystal materials. FZ growth is achieved in a containerless environment (without a crucible) by the movement of a molten zone of a solid rod, which forms a free-standing liquid bridge between two solid sections. However, the surface tension along the free surface of the molten zone varies due to temperature and concentration gradients, and in turn gives rise to Marangoni convection in the melt. It is therefore essential for the growth of high quality crystals by FZ to have a better understanding and precise control of the Marangoni convection occurring in the melt.^1‒3 From this point of view, the role of thermal Marangoni convection in the FZ system has been investigated numerically and experimentally^4‒7 using a half-zone liquid bridge. The Thermal Marangoni convection in liquid bridge of low Prandtl number fluids becomes oscillatory via two-step flow transitions. The first transition occurs from axisymmetric steady to a 3-D steady flow, and then at a larger temperature difference, the second transition to an oscillatory flow takes place. A stability map of thermal Marangoni convection was proposed by Imaishi et al.⁴ In the growth of alloys such as Si_xGe_1-x, it was shown that it is necessary to consider not only the thermal Marangoni convection but also the solutal Marangoni convection due to surface tension differences of the components of the alloy.^8‒13 The flow transformation with the increase or decrease in the Marangoni number is accompanied by the hysteresis phenomena.¹⁴ In this study, the hysteresis phenomenon of flow patterns due to thermal and solutal Marangoni convections in a liquid bridge under zero gravity was numerically simulated.

Figure 1 shows a schematic description of the model domain and the coordinate system used. A half zone (HZ) model is considered as a liquid bridge between two hot and cold discs by establishing a temperature gradient in the volume.

Figure 1 Schematics of liquid bridge.

The model makes the following assumptions:

1. The fluid is incompressible and Newtonian
2. The solid/liquid interfaces are flat and
3. The system is under zero gravity. Since the gravity is considered absent, the liquid bridge was assumed to remain cylindrical (and thus the liquid/gas interface is flat) and the natural convection in the melt does not exist.

Under these assumptions, the governing equations of the liquid phase are obtained respectively from the overall mass conservation, the balance of momentum, the balance of energy and the conservation of mass of species in the following forms:

Continuity

$\nabla ·\overrightarrow{V}=0$ (1)

Momentum

$\frac{\partial \overrightarrow{V}}{\partial \text{t}}+\overrightarrow{V}·\nabla \overrightarrow{V}=-\frac{1}{\rho }\nabla p+\upsilon {\nabla }^2\overrightarrow{V}$ (2)

Energy

$\frac{\partial \text{T}}{\partial \text{t}}+\overrightarrow{V}·\nabla \text{T}=\alpha {\nabla }^2\text{T}$ (3)

Mass Transport

$\frac{\partial \text{C}}{\partial \text{t}}+\overrightarrow{V}·\nabla \text{C}=D{\nabla }^2\text{C}$ (4)

Equations (1)-(4) were discretized by the Finite Volume Method on a collocated grid system, and solved by the pressure-implicit split-operator (PISO) algorithm.¹⁵ The computation was carried out using the open source software OpenFOAM.¹⁶ Euler scheme, second-order quadratic upstream interpolation for convective kinematics scheme and second-order centered scheme were applied respectively to the time term, divergence term and Laplacian term in the governing equations. The residual of the sparse matrix solvers falls below the solver tolerance 1.0×10^-10.

Boundary conditions are:

On the upper-disc (z=L): non-slip condition on the velocity, prescribed temperature (T=T_c) and concentration (C=C_Si=1) values.

On the lower-disc (z=0): non-slip condition on the velocity, prescribed temperature (T=T_h) and concentration (C=C_Ge=0) values.

On the free surface (r=a): surface tension is assumed to be a function of temperature and Si_xGe_1-x concentration only.

$\text{μ}\frac{\partial V_z}{\partial r}=-\left(\frac{\partial \sigma }{\partial T}\frac{\partial T}{\partial z}+\frac{\partial \sigma }{\partial C}\frac{\partial C}{\partial z}\right)$ (5)

$\text{μ}\left\{r\frac{\partial }{\partial r}\left(\frac{V_{\theta }}{r}\right)\right\}=-\frac{1}{r}\left(\frac{\partial \sigma }{\partial T}\frac{\partial T}{\partial \theta }+\frac{\partial \sigma }{\partial C}\frac{\partial C}{\partial \theta }\right)$ (6)

Temperature and concentration conditions along the free surface are:

$\frac{\partial \text{T}}{\partial \text{r}}=\frac{\partial \text{C}}{\partial \text{r}}=0$ (7)

The associated dimensionless thermal and solutal Marangoni numbers are defined by:

$M{a}_T=-\frac{\partial \sigma }{\partial T}\frac{\Delta TL}{\mu \nu },M{a}_C=-\frac{\partial \sigma }{\partial C}\frac{\Delta CL}{\mu \nu }$ (8)

Initial conditions are using following three cases:

1. No-flow case: no-flow and constant values (T = T_c and C = C_Si = 1) in the melt
2. UP case: using result value at 1000s while the solutal Marangoni number is gradually increase $(\Delta M{a}_C=+35.7)$ .
3. DOWN case: using result value at 1000s while the solutal Marangoni number is gradually decrease $\left(\Delta M{a}_C=-35.7\right)$ .

The physical properties of Si_xGe_1-x melt¹⁷ and the operating conditions are listed in Table 1. Computations were performed using the N_R×N_Ө×N_Z=30×120×40 mesh. Non-uniform grids in the region near the upper and lower discs and free surface were used.

Solutal marangoni convection

Figure 2 shows the melt concentration at Ma_T=0 and Ma_C=1786 when the initial condition with No-flow case was considered. The melt concentration becomes time-depended and a 3-D pulsating oscillation of m=4 appears when the solutal Marangoni number is larger than the critical value (Ma_C)_cr. The concentration field oscillates at a period of 60seconds.

Figure 2 Time dependency of the melt concentration at Ma_T = 0 and Ma_C = 1786 (t_p = 60s).

Figures 3 & 4 present the results while the solutal Marangoni number is gradually increasing (UP case). The melt concentration appears a 3-D pulsating oscillation of m=4 as shown in Figure 3 when the solutal Marangoni number is larger than 1072. Figure 4 shows the time dependency of the z-direction velocity at sampling point (r, θ, z)=(a, 0, 0.5L). The flow velocity becomes faster with the increasing solutal Marangoni number. Furthermore, oscillation becomes notable when the solutal Marangoni number is larger than 1072.

Figure 3 Melt concentration distribution (z = 0.5L) at Ma_T=0 and Ma_C=1000(a), 1036(b), 1072 (c) and 1107(d).

Figure 4 Time dependency of z-direction velocity at sampling point (r, θ, z)=(a, 0, 0.5L).

(Figure 5) (Figure 6) present the results while the solutal Marangoni number is gradually decreasing (DOWN case). The flow pattern appears axisymmetric at all cases as shown in Figure 5. However, the flow fields oscillate as shown in Figure 6 when the solutal Marangoni number is smaller than 786.

Figure 5 Melt concentration distribution (z = 0.5L) at Ma_T=0 and Ma_C=857(a), 822(b), 786(c) and 750(d).

Figure 6 Time dependency of z-direction velocity at sampling point (r, θ, z)=(a, 0, 0.5L).

Figure 7 presents the variation of the flow velocity at a sampling point in which the red line (DOWN case) represents the values of the flow velocity while the solutal Marangoni number is gradually decreasing and the blue line (UP case) is for the flow velocity while the solutal Marangoni number is gradually increasing. The flow structure clearly shows hysteresis behaviour of the flow field with about 27% difference between the critical values in the hysteresis diagram.

Figure 7 The amount of change of the velocity at sampling point (r, θ, z)=(a, 0, 0.5L).

Thermal and solutal marangoni convections

Figure 8 shows the melt concentration for the case of thermal and solutal Marangoni convections at Ma_T=2679 and Ma_C=1365. The flow fields became three-dimensional and time-depended when the Marangoni number is larger than the critical value. Figure 9 presents the variation of the flow velocity at a sampling point. As seen in Figure 9, the critical values of the solutal Marangoni number are about 1357(UP case) and 1036 (DOWN case). The flow structure clearly shows hysteresis behaviour of the flow field with about 24% difference between the critical values in the hysteresis diagram. This implies that if we can control the steady flow before becoming unsteady, the thermal and solutal Marangoni convections may be kept under control by the application of weaker external forces such as applied low disc rotations or low magnetic fields.

Figure 8 Time dependency of the melt concentration at Ma_T=2679 and Ma_C=1365 (t_p=15s).

Figure 9 The amount of change of the velocity at sampling point (r, θ, z)=(a, 0, 0.5L).

The present numerical simulation study carried out to examine the hysteresis phenomena of flow patterns due to thermal and solutal Marangoni convections in a liquid bridge under zero gravity led to the following conclusions:

The critical value (Ma_C)_cr depends on the initial condition. The flow structure clearly shows hysteresis behaviour of the flow field with about 24% difference between the critical values in the hysteresis diagram.

None.

Author declares that there is no conflict of interest.

1. Kaiser Th, Benz KW. Floating-zone growth of silicon in magnetic fields III. Numerical simulation. Journal of Crystal Growth. 1998;183(4):564‒572.
2. Dold P, Cröll A, Lichtensteiger M, et al. Floating zone growth of silicon in magnetic fields: IV. Rotating magnetic fields. Journal of Crystal Growth. 2001;231(1‒2):95‒106.
3. Lan CW, Chian JH. Three-dimensional simulation of Marangoni flow and interfaces in floating-zone silicon crystal growth. Journal of Crystal Growth. 2001;230(1‒2):172‒180.
4. Imaishi N, Yasuhiro S, Akiyama Y, et al. Numerical simulation of oscillatory Marangoni flow in half-zone liquid bridge of low Prandtl number fluid. Journal of Crystal Growth. 2001;230(1‒2):164‒171.
5. Zeng Z, Mizuseki H, Simamura K, et al. Three-dimensional oscillatory thermocapillary convection in liquid bridge under microgravity. International Journal of Heat and Mass Transfer. 2001;44(19):3765‒3774.
6. Lappa M, Savino R, Monti R. Three-dimensional numerical simulation of Marangoni instabilities in non-cylindrical liquid bridges in microgravity. International Journal of Heat and Mass Transfer. 2001;44(10):1983‒2003.
7. Sakurai M, Ohishi N, Hirata A. Oscillatory thermocapillary convection in a liquid bridge: Part 1-1g Experiments. Journal of Crystal Growth. 2007;308(2):352‒359.
8. Campbell TA, Shweizer M, Dold P, et al. Float zone growth and characterization of Ge_1−xSi_x (x⩽10 at%) single crystals. Journal of Crystal Growth. 2001;226(2‒3):231‒239.
9. Minakuchi H, Okano Y, Dost S. A three-dimensional numerical simulation study of the Marangoni convection occurring in the crystal growth of Si_xGe_1−x by the float-zone technique in zero gravity. Journal of Crystal Growth. 2004;266(1‒3):140‒144.
10. Minakuchi H, Okano Y, Dost S. A three-dimensional numerical study of Marangoni convection in a floating full zone. Int J of Materials and Product Technology. 2005;22(1/2/3):151‒171.
11. Minakuchi H, Takagi Y, Okano Y, et al. Three-dimensional numerical simulation of thermal and solutal Marangoni convection in a liquid bridge under zero-gravity field. Transactions of JSASS, Aerospace Technology Japan. 2012;10(28):15‒20.
12. Minakuchi H, Takagi Y, Okano Y, et al. The relative contributions of thermo-solutal Marangoni convections on flow patterns in a liquid bridge. Journal of Crystal Growth. 2014;385:61‒65.
13. Minakuchi H, Okano Y, Dost S. Effect of thermo-solutal Marangoni convection on the azimuthal wave number in a liquid bridge. Journal of Crystal Growth. 2017;468:502‒505.
14. Lyubimova TP, Skuridin RV, Faizrakhmanova IS. Thermo- and soluto-capillary convection in the floating zone process in zero gravity conditions. Journal of Crystal Growth. 2007;303:274‒278.
15. Minakuchi H, Takagi Y, Okano Y, et al. A grid refinement study of half-zone configuration of the Floating Zone growth system. Journal of Advanced Research in Physics. 2012;3(1):011201.
16. Open source CFD software. Open FOAM, Open CFD Ltd, UK;
17. Abbasoglu S, Sezai I. Three-dimensional modelling of melt flow and segregation during Czochralski growth of Ge_xSi_1−x single crystals. International Journal of Thermal Sciences. 2007;46(6):561‒572.