The effect of aspect ratio and axial magnetic field on thermocapillary convection in liquid bridges with a deformable free-surface

ABSTRACT

Three-dimensional numerical simulations are performed to analyze the effect of the aspect ratio Ar and axial magnetic field on thermocapillary convection in liquid bridges with deformable free-surface under microgravity, in which the volume of fluid (VOF) method is adopted to track the free-surface movement. The simulation results elucidate that the oscillation wave number m and frequency of temperature fluctuation decrease with increasing Ar, while the amplitude of temperature fluctuation increases with increasing Ar. The deformation ratio ξ of the free-surface increases as Ar increases. The numerical results also reveal that the axial magnetic field causes a concentration of convection vortexes near the free-surface and effectively suppresses the flow in both the radial and axial directions. Moreover, the axial magnetic field effectively damps free-surface deformation, so that the deformation ratio ξ decreases as the Hartmann number Ha increases. The temperature displays a uniform and linear distribution along the free-surface and the axis of molten zone when B_a = 0.3 T.

KEYWORDS:

• thermocapillary convection
• deformable free-surface
• microgravity
• VOF
• axial magnetic field

Nomenclature

Ar	=	aspect ratio, Ar = L/r1
A	=	average oscillation amplitude,
b	=	induced magnetic field vector (T)
B	=	magnetic field vector (T)
Ba	=	intensity of the applied axial magmatic fields (T)
Ca	=	Capillary number,
d	=	distance of the cell center apart from the midline EF
fl	=	Lorentz force (N·m−3)
f	=	frequency at monitoring point N
g	=	gravitational acceleration (m2·s−1)
Ha	=	Hartmann number,
J	=	induced current (A·m−2)
k	=	unit tangential vector
L	=	length of the liquid bridge (cm)
Ma	=	Marangoni number,
m	=	hydrothermal wave number
n	=	unit normal vector
N	=	monitoring point
P	=	pressure (Pa)
Pr	=	Prandtl number,
Prm	=	magnetic Prandtl number,
Q	=	convergence criteria
r	=	radius of the liquid bridge (cm)
rmax	=	maximum radius of the free-surface (cm)
rmin	=	minimum radius of the free-surface (cm)
Re	=	Reynolds number,
Rem	=	magnetic Reynolds number,
t	=	time (s)
T	=	temperature (K)
Tc	=	temperature of cold disk (K)
Th	=	temperature of hot disk (K)
Tav	=	average temperature of half-zone (K)
TNav	=	average temperature at the monitoring point N (K)
ΔT	=	temperature difference, ΔT = Th−Tc (K)

Greek Symbols

μ	=	dynamic viscosity (kg·m−1·s−1)
υ	=	kinematic viscosity (m2·s−1)
χ	=	thermal diffusivity (K−1)
ρ	=	density (kg·m−3)
Cp	=	specific heat (J·kg−1·K−1)
λ	=	thermal conductivity (W·m−1·K−1)
δ	=	thickness of the air layer
σ	=	surface tension (N·m−1)
ξ	=	deformation ratio,
	=	temperature coefficient of surface tension (N·m−1·K−1)
	=	electric conductivity
u, v, w	=	the velocity components in the r, z and θ directions (cm·s−1)
r, θ, z	=	radial, azimuthal and axial coordinates

Subscripts

a	=	axial magmatic field
av	=	average
c	=	cold wall
h	=	hot wall
max	=	maximum
min	=	minimum
1	=	inner melt fluid
2	=	outer air fluid

1. Introduction

Single crystal silicon is an important semiconductor material that is in high demand in electronics, information technologies and optoelectronic devices, etc. due to its high performance and low cost (Lappa & Savino, 2002; Takagi, Okano, & Dost, 2012). Therefore, the compositional uniformity of and absence of defects in the grown crystals is extremely important. The floating zone (FZ) method for growing higher-quality single crystals is considered very promising under microgravity. However, the effect of thermocapillary convection induced by the surface tension gradient on the crystal quality is an issue. Furthermore, the instability of thermocapillary convection is the major cause of macroscopic defects during the FZ crystal growth process (Sankar, Venkatachalappa, & Do, 2011). In order to produce homogeneous bulk silicon crystals under microgravity, the oscillatory flow must be theoretically analyzed and its instability must be well controlled.

Schwabe and Frank (1999) found that oscillatory thermocapillary convection was one of the sources of crystal growth striations. Since the above-mentioned discovery, many studies have been conducted on thermocapillary flow in liquid bridges with a non-deformed free-surface. Lappa, Savino, and Monti (2001) found that the azimuthal wave number is a monotonically decreasing function of the geometrical aspect ratio. Yano, Nishino, Kawamura, Ueno, and Matsumoto (2015) experimentally investigated the instability and associated roll structures of Marangoni convection in liquid bridges on the International Space Station. The roll structures traveled in the same direction as the surface flow (co-flow direction) for 1.00 ≤ AR ≤ 1.25 (where AR = height/diameter), while they traveled in the opposite direction (counter-flow direction) for AR ≥ 1.50. Melnikov, Shevtsova, Yano, and Nishino (2015) investigated the Marangoni convection in liquid bridges for a wide range of aspect ratios by experiment and simulation. Their results illustrated that the heat loss has a destabilizing effect on aspect ratios (ratio of radius to height) below 2.4.

Up to now, the oscillatory thermocapillary flow in liquid bridges with a deformable free-surface has only been reported on by a few researchers. Sim, Kim, and Zebib (2004a, 2004b) investigated the dynamic free-surface deformations in a two-dimensional (2D) axisymmetric liquid bridge. The free-surfaces were convex near the cold wall, while near the hot wall they changed from concave to convex with an increasing Reynolds number (Re). A static deformation free-surface was analyzed by Lappa (2005), in which the free-surface was assumed to be a convex or concave shape. The simulation results demonstrated that the stability of the flow field was very sensitive to the geometrical aspect ratio A_F (length/diameter) and free-surface shape (convex or concave). Zhou and Huai (2015a, 2015b) investigated the thermo-solutocapillary convection in an open rectangular cavity and a liquid bridge with a dynamic free-surface. The free-surface bulged out at the two ends of the zone and bulged in at the central zone when the thermal Marangoni number was identical to the solutal Marangoni number. Yang, Liang, Yan, Gao, and Feng (2015) investigated thermocapillary convection and surface fluctuation in a liquid bridge under lateral vibrations which inhibited the surface flow.

Due to the high conductivity of silicon-melt, the application of a magnetic field is the most appropriate method for controlling and damping the deleterious flow during the crystal growth process. Many kinds of magnetic fields, such as the axial magnetic field, transverse magnetic field, rotating magnetic field and non-uniform magnetic field, have been performed. Substantial theoretical and numerical works have been dedicated to the magnetic damping of thermocapillary convection (Yao, Zeng, Chen, & Li, 2012; Yao, Zeng, & Li, 2011; Zhou & Huang, 2010). For the axial magnetic field, Liang, Yang, and Li (2014) found that the Lorentz force strongly damped the fluid flow not only in the radial direction inside the liquid bridge but also in the circumference direction near free-surface. The influence of a static axial magnetic field was also performed by Jaber, Saghir, and Viviani (2009), and they found that the convective flow was efficiently damped by the axial magnetic field. However, most simulations neglected the effect of free-surface deformation.

Although few works have been performed which have considered free-surface deformation, most of the studies have focused on a 2D axisymmetric liquid bridge or static deformable surface. The flow and free-surface deformation behaviors in 3D liquid bridges and the impact of the axial magnetic field on free-surface deformation are not clearly illustrated. The very existence of these issues result in an intellectual curiosity that needs to be satisfied. Thus, the purpose of present works is to investigate and actively control the transport behaviors (flow field, temperature field and free-surface deformation) which occur during the crystal growth process. Both the effect of Ar and the applied axial magnetic field are analyzed with regard to the instability of thermocapillary convection.

2. Physical and mathematical models

The physical configuration of the present study consists of two concentric columns (inner and outer) of immiscible fluids, as shown in Figure 1(a). The inner column fluid, the characteristic length of which is L and the radius of which is r₁, limiting the aspect ratio of Ar = L/r₁, is suspended between two coaxial disks at a temperature difference of ΔT = T_h−T_c, and the fluid material is silicon-melt. The lower disk (hot wall) is maintained at a constant temperature T_h, while the upper disk (cold wall) is at a lower temperature T_c. The outer column fluid of height H and thickness δ = r₂−r₁ has air as its fluid material. A cylindrical coordinate system (r, θ, z) is adopted, with the z-axis along the vertical centerline of the cylindrical liquid bridge and the origin located at the center of lower disk. The magnetic lines of the external axial magnetic field are assumed to permeate the whole molten zone without change, as shown in Figure 1(b).

Figure 1. Sketch of the liquid bridge and magnetic lines of the axial magnetic field: (a) coordinate system, and (b) magnetic lines.

The following assumptions are introduced in our models: (1) the two kinds of fluids are incompressible Newtonian fluids; (2) the flow is laminar; (3) the system is under microgravity (set at 0g in the present investigation, where g denotes gravitational acceleration); and (4) the interface of the melt and the air is deformable. The surface tension of the interface is assumed a linear function of temperature (1) where σ₀ is the surface tension for T = T₀, T₀ indicates the reference temperature of the initial surface tension (which is equal to the melting point temperature; Lappa, 2005; Zeng et al., 2001), and is the negative rate of change of the surface tension with temperature.

2.1. Governing equations

With the above assumptions, the non-dimensional governing equations are obtained respectively from the overall continuity, momentum, energy (i = 1, 2; where 1 denotes the inner fluid and 2 denotes the outer fluid): (2) (3) (4) where the variables are non-dimensionalized as , , , , , , , , . And , , , and are the dimensionless dynamic viscosity, thermal conductivity, density and specific heat, respectively.

The induced magnetic field b can be derived from Maxwell's equations and Ohm's law as follows: (5) where B* is the magnetic field intensity, which includes both the applied and the induced magnetic field b* component. The induced current density is and the Lorentz force is expressed as f_l = J × B.

The dimensionless parameters Reynolds, magnetic Reynolds, Prandtl, magnetic Prandtl, Marangoni, capillary and Hartmann numbers are given by:

2.2. Free-surface deformation

Figure 2 presents the sketch of the free-surface and boundary conditions in the vertical section. The position of the interface is described by a function F(z), where F(0) = F(L) = r₁ = R₀, and the initial shape of the interface is assumed to be a cylindrical surface. The stress balance between the silicon-melt and the air on the free-surface, R = F(z), is controlled by (6) where S_ik, R_i (i = 1, 2) are the tensor of viscous stress and normal radii of curvature, P and P₀ are the pressure of the inner and outer fluid, and k and n are the unit tangential vector and outward unit normal vector, respectively. The mean curvature is defined as (7) where the subscript z denotes the derivative F_z = dF/dz. The normal and tangential projections of Equation (7) are given by (8) (9)

Figure 2. Sketch of the free-surface and boundary conditions in the vertical section.

The location of the interface is determined by the normal projection Equation (8). The tangential projection Equation (9) defines the driving thermocapillary force. Since the liquid is assumed to be incompressible, its total volume must remain constant: (10)

2.3. Boundary and initial conditions

The boundary conditions are written as follows: (11a) (11b) (11c) (11d)

The initial condition is expressed as (12)

In the present works, the thermophysical properties of the two kinds of fluids used are readily available in Table 1, the properties of which are in accordance with Dold, Cröll, Lichtensteiger, Kaiser, and Benz (2001).

Table 1. Thermophysical properties of the two kinds of fluids.

Fluids	μ Pa·s	ρ kg·m^−3	Cp J·kg^−1·K^−1	λ W·m^−1·K^−1	σm Ω^−1·m^−1	σ0 N·m^−1	σT N·m^−1·K^−1
Silicon-melt	8.89 × 10^−4	2530	1000	67	1.29 × 10^6	0.4868	2.8 × 10^−4
Air	1.78 × 10^−5	1.225	1006.43	0.0242	0

3. Numerical method

The volume of fluid (VOF) method is widely adopted to track the free-surface movement (Chokri, Zouhaier, Mohamed, & Khlifa, 2010; Zhou, Liu, & Ou, 2011), and it is also applied in the present numerical strategy. The governing equations are discretized by the finite volume method. The central difference approximation is applied to the diffusion terms. The first-order upwind scheme for the convection term is used, and the flow field is solved using the SIMPLE algorithm. The surface axial velocity distributions at B_a = 0.5 T obtained from present work (using the same physical properties and parameters as Liang et al., 2014) are compared with the work of Liang et al. (2014), and an expectantly close agreement is observed for this case, as shown in Figure 3.

Figure 3. Surface axial velocity distribution at B_a = 0.5 T, Ar = 1 and Pr = 61.

The grid system employed in the present computations is an axisymmetric and non-uniform grid. A finer mesh gradation is utilized in the regions where strong velocity and temperature gradients are anticipated, which includes the regions near the upper and lower disks, as well as the interface. In order to check the grid independence, simulations with four sets of different meshes, as shown in Table 2, are performed at Ma = 1.64 × 10⁴, B_a = 0.5 T, Ar = 1.5. The last two sets of grid (M3 and M4) show almost the same average axial velocity and temperature at the monitoring point N(r/r₁ = 0.5, θ = 0°, z/L = 0.5). The choice of grid 80^r × 120^θ × 80^z (M3) corresponds to a reasonable compromise between the minimum dependence of the threshold of the instability upon the mesh size and the computational time. For each time step (which has been validated using the adapted calculation time step), the quotient is calculated for all dependent variables in all grid points. If Q ≤ 10⁻⁵ is valid for all variables in all grid points, the solution is interpreted as converged. All calculations were performed using six CPUs on the local Beowulf-Linux cluster.

Table 2. Main characteristics of the different meshes for Ma = 1.64 × 10⁴, B_a = 0.5 T, and Ar = 1.5.

Mesh	M1	M2	M3	M4
Number of nodes	70^r × 100^θ × 70^z	75^r × 110^θ × 75^z	80^r × 120^θ × 80^z	85^r × 130^θ × 85^z
vav (cm·s^−1)	1.04 × 10^−2	1.73 × 10^−2	1.85 × 10^−2	1.86 × 10^−2
TNav (K)	1682.279	1682.438	1682.534	1682.535

4. Results and discussion

4.1. The effect of the aspect ratio on the oscillatory flow

In order to explore the spatio-temporal details of the flow pattern and free-surface deformation of a liquid bridge, and also understand the impact of Ar on the oscillatory thermocapillary flow, numerical simulations were conducted based on liquid bridge models with values of Ar in the range of 1 to 2. Here, the radius (r₁) of the two disks is 1 cm and the height (L) of the molten zone is from 1 cm to 2 cm. Other parameters adopted are as follows: ΔT = 45 K, Ma = 1.64 × 10⁴, and Ca = 0.026.

Any temperature difference (ΔT > 0) may produce a surface tension gradient at the free-surface, which is called the Marangoni phenomenon. The thermocapillary induced by the Marangoni effect flows from the hot wall towards the cold wall along the interface and returns back to the center (see Figure 4). As can be seen from the streamlines in Figure 4(a, b, and c), the whole molten zone appears to be almost entirely occupied by two larger convection vortexes with a clockwise circulation in the left half of the domain and a counterclockwise circulation in the right half of the domain. The streamlines all show a narrow ‘neck-shaped’ structure strongly concentrated near the interface. However, the differences in flow pattern are also detectable for different values of Ar. At Ar = 1 (see Figure 4(a)), two vortexes (C and D) form in the molten zone. Due to the small deformation of the free-surface, the vortexes C and D are located close to the midline of EF; with Ar increasing to 1.25 (see Figure 4(b)), the two vortexes stretch along the axial direction. Owing to the increase of free-surface deformation, the distance d of the vortexes (C and D) away from the midline of EF increases; when increasing the value of Ar further to 1.5 (see Figure 4(c)), the distance d also increases and the main vortexes break up into two pairs of vortexes. The vortexes of G and H at the upper region are the newly emerged vortexes. Near the interface and in the center of the liquid bridge, the isotherms all perform a non-uniform and distortion distribution (see Figure 4 (d, e, f)). This is because the surface flow is strong near the interface and the return flow is also strong at the center of molten zone, which makes an obvious distortion of the isotherms. Near the hot and cold walls, the isotherms perform in a relatively flat and uniform manner because the convection heat transfer is weak and the heat conductivity is predominant. Comparing the isotherms of the three cases, they also exhibit similar structures. Nevertheless, with the increasing values of Ar, the narrow low temperature area becomes sharper at the center (indicated by the arrow zone in Figure 4(d, e, f)) and the extent of the bulge decreases at the free-surface. In addition, the extent of the contraction of the isotherms increases with increasing values of Ar, because the degree of inward shrinking of free-surface increases with increasing values of Ar.

Figure 4. Streamlines (a, b, c) and isotherms (d, e, f) on the θ = 0° plane for (a, d) Ar = 1, (b, e) Ar = 1.25, and (c, f) Ar = 1.5.

For a better understanding of the effect which Ar has on the instability melt flow, the radial velocity distributions for different values of Ar on the z/L = 0.5 plane are shown in Figure 5. It was found that the radial velocity distribution presents a ‘petal shape’ on the middle plane of the liquid bridge. At Ar = 1, the ‘petal’ number is 12, indicating that the oscillation wave number m is 12. But some ‘petals’ are pretty blurry in the region from θ = 90° to θ = 150° (see Figure 5(a)) because of the oscillatory and deformable free-surface. At Ar = 1.25, the oscillation wave number m decreases to 8 (see Figure 5(b)) and as Ar increases to 1.5, the oscillation wave number m decreases to 6 (see Figure 5(c)). Thus, at Ma = 1.64 × 10⁴ and Ca = 0.026, hydrothermal wave instability of the thermocapillary convection occurs in the whole liquid bridge and is weakened as Ar increases, i.e., the wave number m decreases from 12 to 6 as Ar increased from 1 to 1.5 (see Figure 6), demonstrating that the wave number decreases as Ar increases. Moreover, there is a negative velocity zone at the interface because the free-surface shrinks inward in the middle section of the molten zone.

Figure 5. Radial velocity distributions on the z/L = 0.5 plane for (a) Ar = 1, (b) Ar = 1.25, and (c) Ar = 1.5.

Figure 6. Oscillatory wave number m varying with Ar.

In order to provide more insight into the above observations, the axial and radial velocity distribution along the circumferential direction at r/r₁ = 0.5 on the z/L = 0.5 plane for different values of Ar are exhibited in Figure 7. From examining the axial velocity in Figure 7(a), obvious oscillation behaviors are found on the three axial velocity oscillatory curves. The axial velocity of the three curves is negative, indicating that the flow direction is from the upper disk towards the lower disk. At Ar = 1, the velocity oscillation is disorder, ad hoc in the region from θ = 40° to θ = 150°. As Ar increases, the amplitude of the axial velocity increases and the peak number decreases. Strong oscillation behaviors are also found on the radial velocity oscillatory curves of Figure 7(b). At Ar = 1, the peak number is 12 and the amplitude of the radial velocity is small, but the wave crests are pretty blurry in the region from θ = 90° to θ = 150°. As Ar increases, the amplitude of radial velocity increases and the peak number decreases from 12 to 6 because, for small values of Ar, the surface tension gradient is higher, inducing stronger instability in the thermocapillary convection in the liquid bridge. It is worth mentioning that the peak number of the radial velocity of Figure 7(b) is consistent with the oscillation wave number m of Figure 5. For example, the peak number of radial velocity is 8 at Ar = 1.25 (see Figure 7(b)), and the oscillation wave number m is also equal to 8 (see Figure 5(b)).

Figure 7. Velocity distribution along the circumferential direction at r/r₁ = 0.5 on the z/L = 0.5 plane for different values of Ar for (a) the axial velocity, and (b) the radial velocity.

Figure 8 presents the transient temperature at the monitoring point N(r/r₁ = 0.5, θ = 0°, z/L = 0.5), varying with time for different values of Ar. Strong oscillatory behaviors emerge on the temperature fluctuation curves. At Ar = 1 (see Figure 8(a)), the frequency of temperature fluctuation is highest while the period average amplitude of the transient temperature is lowest in that of the three values of Ar. The period average amplitude of transient temperature obviously changes in different periods, indicating that a second transition of thermocapillary flow occurs. When Ar is 1.25 (see Figure 8(b)), the frequency of temperature fluctuation decreases to 1.367 Hz while the period average amplitude of the transient temperature increases to 1.2 K; when Ar increases to 1.5 (see Figure 8(c)), the frequency of temperature fluctuation further decreases to 1.2 Hz and the period average amplitude of the transient temperature further increases to 1.5 K, indicating that the instability of thermocapillary convection decreases. The period average temperature and frequency at the monitoring point N varying with different values of Ar are shown in Figure 9. As Ar increases from 1 to 1.5, the frequency of the monitoring point N decreases from 1.67 Hz to 1 Hz, while the average temperature of the monitoring point N increases from 1679.1 K to 1682.73 K.

Figure 8. Temperature fluctuations at monitoring point N for (a) Ar = 1, (b) Ar = 1.25, and (c) Ar = 1.5.

Figure 9. Average temperature and frequency at monitoring point N varying with Ar.

Figure 10 illustrates variations of the free-surface deformation F(z) (Figure 10(a)) and deformation ratio ξ (Figure 10(b)) varying with values of Ar (the deformation ratio ξ is defined as ξ = (r_max−r_min)/r₁, where r_max is the maximum radius of the free-surface and r_min is the minimum radius of the free-surface). As can be seen in Figure 10(a), the free-surfaces all bulge outward at the two ends of the half-zone and shrink inward in the middle section of the molten zone. In other words, it is convex near the two ends of the half-zone and concave in the middle region of the liquid bridge. As a result, the free-surface presents a typical narrow ‘neck-shaped’ structure, which shows the same behavior as the streamlines in Figure 4(a, b, c). All four lines keep the same tendency, the maximum bulge points whose axial positions are z/L = 0.15 and z/L = 0.95 and the maximum contraction point whose axial positions is z/L = 0.63, while the extent of the maximum bulge and constriction points increases as Ar increases in the radial position. For instance, at Ar = 1, the radial position of the maximum convex point is r_max/r₁ = 1.058 and the maximum concave point is r_min/r₁ = 0.93. However, at Ar = 1.75, the radial position of the maximum convex point shifts outward to r_max/r₁ = 1.08 and the maximum concave point shrinks inward to r_min/r₁ = 0.878. This shows that the extent of the free-surface deformation increases as Ar increases, which can also clearly be seen in deformation ratio ξ from Figure 10(b). The deformation ratio ξ increases from 0.14 to 0.22 as Ar increases from 1 to 2. Moreover, the volume of convex at two ends of the zone is equal to that of the concave in the central region for different values of Ar, which makes the total volume of deformed liquid bridge equal to that of the initial cylindrical one (see Figure 10(a)). That is to say, the volume is conserved. Similarly, at θ = 90°, the free-surface shows the same behavior as that when θ = 0°.

Figure 10. (a) Free-surface deformation F(z) and (b) deformation ratio ξ varying with Ar on the θ = 0° plane.

4.2. The effect of the axial magnetic field on the oscillatory flow

Since the thermocapillary convection and its instability are the major causes of some macroscopic defects in the crystal produced, it is desirable to suppress this deleterious flow. One of the effective methods of controlling convection flow is magnetic damping, in which the applied magnetic field generates a Lorentz force to inhibit the unstable flow. In order to investigate the effect of the axial magnetic field on the melt convection and temperature fields as well as the instability of the thermocapillary flow in a liquid bridge with a deformable free-surface, the liquid bridge models are kept with the constants of Ar = 1.5, Ma = 1.64 × 10⁴, ΔT = 45 K, and Ca = 0.026. The intensity of the applied axial magnetic field B_a varies from 0 to 0.3 T, and the corresponding Hartmann number Ha is 0 to 171.8.

Figure 11 presents the streamlines and isotherms on the θ = 0° plane at different intensities of the magnetic field B_a. Since the damping effect of Lorentz forces of axial magnetic fields reduces the volume transport of the thermocapillary vortex along the radial direction, and the damping effect increases with increasing B_a, the thermocapillary vortexes are confined near the vicinity of the free-surface – that is, the axial magnetic field causes a concentration of convection cells in the vicinity of the free-surface. Comparing Figure 11(a, b, c) to Figure 4(c), in the absence of a magnetic field, the main vortexes occur in the lower region and the newly formed vortexes occur in the upper region of the molten zone. At B_a = 0.1 T, the vortexes C and D move gradually towards the upper region and the broken-off vortexes G and H merge together with C and D (see Figure 11(a)). At B_a = 0.2 T, the distance d between the vortexes (C and D) and the midline of EF increases. The vortexes further concentrate towards the interface (see Figure 11(b)). When B_a = 0.3 T, the distance d further increases and the main vortexes C and D are broken into two pairs of new vortexes again by further compression of the axial magnetic field. The vortexes are maintained at the interface and the narrow ‘neck-shaped’ structure gradually disappears (see Figure 11(c)). Similarly, an obvious distortion of isotherms is found near the free-surface and also at the center of molten zone in the case when B_a = 0.1 T (see Figure 11(d)), which is similar to that without axial magnetic field (Figure 4(f)). However, the distortion of the isotherms diminishes gradually as B_a increases (see Figure 11(e, f)). Because the flow structure and heat transport are dominated by Lorentz force, the thermocapillary convection is also weakened. Thus, the isotherms in the liquid bridge tend to be flatter and become a uniform and regular distribution.

Figure 11. Streamlines (a, b, c) and isotherms (d, e, f) on the θ = 0° plane for (a, d) B_a = 0.1 T, (b, e) B_a = 0.2 T, and (c, f) B_a = 0.3T.

Figure 12 illustrates the axial and radial velocity distribution along the circumferential direction at r/r₁ = 0.5 on the z/L = 0.5 plane for different values of B_a. From Figure 12(a) showing the axial velocity oscillatory curves, a strong oscillation of axial velocity is found and the peak number is 7 in the case without the axial magnetic field. However, the amplitude of axial velocity decreases as B_a increases. When the intensity of the magnetic field is at 0.3 T, the axial velocity is fully damped, the wavelike state does not appear, and the velocity is almost 0. The radial velocity distribution of Figure 12(b) shows similar behavior to the axial velocity distribution of Figure 12(a). The radial velocity is also efficiently damped when the axial magnetic field is applied. The physical mechanism occurring here is that of the Lorentz force strongly damping thermocapillary convection under the effect of the axial magnetic field. Therefore, the oscillatory thermocapillary convection is effectively suppressed along both the radial and axial directions.

Figure 12. Velocity distribution along the circumferential direction at r/r₁ = 0.5 on the z/L = 0.5 plane with different values of B_a for (a) the axial velocity, and (b) the radial velocity.

The velocity vector of thermocapillary convection along the axial direction does not produce a velocity component along the radial direction when the flow is axisymmetric. When the flow is oscillation and non-axisymmetric, thermocapillary flow produces the velocity component along the radial direction and the radial velocity u produces inhomogeneity along the circumference. The average oscillation amplitude A along the circumference is calculated as (Chen & Li, 2008): where N is the number of control volumes along the circumference direction. Figure 13 gives the values of A along the circumference at r/r₁ = 0.5 on the z/L = 0.5 plane under different values of B_a. It can be seen that the value of A decreases linearly with as Ha increases. The critical magnetic field strength B_a can be determined using a linear regressive extrapolation method and the point of A = 0 is the steady state of thermocapillary convection. Therefore, the critical magnetic field strength is Ha = 172.5, and the corresponding B_a is 0.301 T.

Figure 13. Average amplitude of the radial velocity varying with Ha at r/r₁ = 0.5 on the z/L = 0.5 plane.

In order to investigate the effect of the intensity of magnetic field on free-surface deformation, the free-surface deformation F(z) and deformation ratio ξ varying with Ha are shown in Figure 14. In the case without the axial magnetic field, the axial positions of the maximum convex and maximum concave points are z/L = 0.13 and z/L = 0.6, respectively. When the axial magnetic field is applied, the extent of both the bulge and the contraction becomes smaller, and the convex region moves inward while the concave region moves towards cold wall, as shown in Figure 14(a). The extent of the free-surface deformation decreases as B_a increases and the free-surface tends to be a straight line, which can be confirmed by the deformation ratio ξ in Figure 14(b). When Ha increases from 0 to 114.5, the deformation ratio ξ decreases rapidly, and as Ha increases from 114.5 to 171.8, the deformation ratio ξ increases slowly. As Ha = 171.8 (B_a = 0.3 T), the slope of the deformation ratio curve is almost 0 and the deformation ratio ξ is 0.03, as shown in Figure 14(b). Therefore, the critical magnetic field strength B_a is approximately 0.3 T, which is almost consistent with the critical magnetic field strength in Figure 13.

Figure 14. (a) Free-surface deformation F(z) and (b) deformation ratio ξ varying with Ha on the θ = 0° plane.

The influence of the axial magnetic field on the temperature distributions along the free-surface and axis of half-zone is also investigated. Figure 15 elucidates the free-surface and axis temperature distributions under a magnetic field strength B_a from 0 to 0.3 T. In the absence of the magnetic field, both the free-surface and axis temperature present ‘bow-shaped’ distributions. The ‘bow’ shape of temperature distribution on the free-surface bends upward (see Figure 15(a)), while that along the axis bends downward (see Figure 15(b)). With increasing values of B_a, the ‘bow’ shape gradually disappears. As B_a = 0.3 T, the temperature displays a uniform and linear distribution along the free-surface and axis of molten zone, indicating that the thermocapillary flow is efficiently damped by the axial magnetic field.

Figure 15. Temperature distributions along (a) the free-surface of the liquid bridge and (b) the axis of the liquid bridge for different values of B_a.

5. Conclusion

In the present paper, 3D numerical simulations are carried out to analyze the effect of the aspect ratio Ar and the axial magnetic field B_a on oscillatory thermocapillary convection and free-surface deformation during the single crystal process, in which the VOF method was employed to track the free-surface movement. From the results, the following conclusions can be drawn:

1. The effect of Ar plays an important role with the flow and temperature fields. There are strong oscillations in both the axial and radial velocities in the whole molten region, and the oscillation wave number of the radial velocity m decreases from 12 to 6 as Ar increases from 1 to 1.5. Moreover, as Ar increases, the frequency of temperature fluctuation decreases and the amplitude increases.

2. The free-surface is convex near the two ends of the liquid bridge and concave in the middle section, and the deformation ratio ξ increases as Ar increases.

3. The axial magnetic field causes a concentration of convection cells near the free-surface and effectively suppresses the occurrence of oscillatory convection. It also damps the convection flow along both the radial and axial directions.

4. The axial magnetic field is helpful for damping free-surface deformation. The deformation ratio ξ decreases as Ha increases, and the critical magnetic field strength is B_a = 0.301 T in this case. In addition, the temperature displays a uniform and linear distribution along the free-surface and axis of molten zone when B_a = 0.3 T.

Our present works are limited to the investigations, and related experimental work will be considered in the future. In addition, the effect of Marangoni number, rotating magnetic field and non-uniform magnetic field on thermocapillary convection in deformable liquid bridges will also be investigated in our future works. We will aim to determine the critical Marangoni number of the critical temperature difference.

Acknowledgements

This work was supported by the National Nature Science Foundation of China [grant number 51276089].

Disclosure statement

No potential conflict of interest was reported by the authors.