Mixed Oscillation Flow of Binary Fluid with Minus One Capillary Ratio in the Czochralski Crystal Growth Model

Abstract

This work presented a series of three-dimensional unsteady numerical simulations on the characteristics of the mixed oscillation flows of binary mixture in a Czochralski crystal growth model. The silicon-germanium melt is investigated and the capillary ratio is minus one. The simulation results showed that, for the special capillary ratio, the thermal and solutocapillary forces are imposed in opposite directions and counteract each other. With the effect of buoyancy, the balance between the capillary forces is disturbed. Mixed with the forced convection driven by rotation, the capillary-buoyancy convection is complex. The basic mixed flow streamlines are presented as various rolling cells. The directions of the rolls are dependent on the combinations of surface and body forces. With the increase of temperature gradient, the basic flow stability is broken, and the oscillations occur. The crucible rotation has an effective influence on the stability enhancement. However, affected by the crystal rotation, the critical condition experiences an increase to a turning point, and then undergoes a sharp reduction to zero. Once the instability is incubated, the surface oscillations are analyzed. For the three-dimensional steady flow, only spatial oscillations are observed circumferentially, and the surface patterns of spokes, rosebud, and pulsating ring are obtained. For the unsteady oscillation flow, the spiral hydrosoultal waves, rotating waves, and superimposition of spirals and spokes are observed, and the oscillation behaviors are also discussed.

1. Introduction

With the rapid development of information technology, the electronics and optoelectronic materials have attracted much attention in various research fields. As important representatives of photoelectric functional materials, the supplement of synthetic crystals increases continuously, and the improvement on the functions and qualities of the crystals are required [1,2]. Among the crystal growth methods, the Czochralski (Cz) technique is one of the most important methods employed for crystal production [3,4]. In the process of crystal growth by the Czochralski method, the temperature and concentration gradients induce the heat and mass transfer, and drive the diffusion and natural convections (including the thermo-solutal-capillary convection and capillary-buoyancy convections). Meanwhile, for the uniformity, the rotations of crystal and crucible are applied [5,6]. As a result, the forced convection is formed and simultaneously mixed with the natural convections. The mixed flow is very complex as the surface capillary forces and body forces interact on different scales, making it difficult for the flow stability to be well recognized and controlled [7,8]. However, the stability of the mixed convection plays a key role in the growth of high-quality crystals. When the flow instabilities occur, the oscillations in the flow field will lead directly to the crystallographic defects. Therefore, the characteristics of the complex convections during the Czochralski process are required to be investigated [9,10,11].

McTaggart [12] firstly studied the flow instabilities with vertical temperature and solute concentration gradients in an infinitely extended layer, through the linear stability analysis. It was found that the occurrence of instabilities is associated with the directions of the capillary forces. For the positive solutal Marangoni number, the convective flow cell is static. When the solutal effect is negative, the flow is oscillated after losing the stability. Ho and Chang [13] extended MaTaggart’s research to the nonlinear instabilities induced by the thermo-solutocapillary effect. It was reported that the steady rolls with nonlinear finite amplitude has a depression influence on the oscillations. Castillo [14,15] studied the Marangoni–Bénard flow in a mixture layer. The linear stability diagram was obtained. Bergeon [16,17] conducted two-dimensional (2D) and three-dimensional (3D) simulations to investigate the mixed thermo-solutal convections in a rectangular cavity. The Soret effect was considered and the results showed that solute concentration gradient drives the fluid flow. Afterwards, a great deal of studies have been performed to investigate the thermo-solutocapillary flows [18,19,20,21,22,23].

For the convections induced by horizontal concentration and temperature gradients, Bergman [24] defined the capillary ratio to express the relative magnitude of the solutal and thermo-capillary forces, which is analogous to the buoyancy ratio in the double diffusion system. Chen [25] and Li [26] performed linear stability analysis on the double-diffusive convections. They found that, with the variation of the Marangoni number, the flow transitions in the liquid pool belong to the Hopf bifurcation. Zhou and Huai [27,28] conducted numerical simulations on the thermo-solutocapillary flow in an opened cavity. The results showed that the deformation of the surface is closely linked with the capillary ratio. For the minus capillary ratio, an S-shaped surface is observed.

Meanwhile, rotation driven flow is a traditional and hot topic in the research field of fluid mechanics [29,30,31]. Zebib [32] pointed out the rotation effect is not ignorable in the investigation of thermal instabilities during the crystal growth process. Tian et al. [33] considered the coupled rotation and thermal flow in an annular pool through linear stability analysis. Energy balance analysis was explored to capture the flow instabilities, and a stability diagram was presented. It was shown that the counteraction of the rotation and thermocapillary effect induces three Hopf bifurcation points for the flow transition. Meanwhile, we have previously conducted several experiments [34] and simulations [35,36,37] to investigate the effect of rotation on the surface tension driven flow in pure liquid. The results showed that the combinations of the rotation and thermocapillary forces drive several kinds of unstable flow. The mechanisms of flow instabilities were revealed.

For the compound semiconductor crystal growth, the coupling thermal and solutocapillary-buoyancy effects in the mixture are rather more complex than that of pure fluid. Cröll et al. [38] conducted an experiment during the parabolic flight campaign to separate the buoyancy effect from the capillary convection. The solutocapillary convection was conclusively observed during the solidification of silicon-germanium. Also, Campbell et al. [39] reported that, during the crystal growth, the convection in a mixture is different from that of pure liquid melt. Especially for the silicon-germanium mixture, the large surface tension and density differences between silicon and germanium induce strong solutal and buoyancy convections. These convections are coupled with a thermal effect, which directly leads to the sharp change of the crystal interface. After that, a series of experiments and simulations were conducted to qualify the coupling effect of the thermal and solutal-capillary effect in the silicon-germanium mixture [40,41]. Abbasoglu et al. [42] pointed out that the concentration profiles is closely related with the temperature through the radial segregation. For the Czochralski growth of silicon-germanium crystal, the capillary ratio changes from −2.44 to 0.2. However, in their work, the effect of rotation was not considered. Once the rotation is taken into account, the centrifugal force is introduced and interacted with the capillary forces. For the mixture with small capillary ratio of −0.2 [41], the thermocapillary effect is dominant, and it was proved that rotation has similar a influence as that in pure fluid. However, for the mixture with minus one capillary ratio, the thermal and solutocapillary forces are counterbalanced. Once any disturbance, such as buoyancy or rotation, is introduced, the balance will be broken. The stability of such a mixed flow is still unclear. This work is devoted to the understanding of the mixed flow of binary fluid with minus one capillary ratio in a Czochralski crystal growth model. As the GeSi material [43] has many unique physical properties (such as high mobility), it has important application value in optoelectronics. Therefore, the GeSi melt is adopted as the test fluid in this work.

2. Model Formulation and Numerical Methodology

2.1. Physical Model and Basic Assumption

The sketch of the physical model investigated in this work is provided in Figure 1. For a comparison with the pure fluid, the same geometry is adopted [34,35]. The crucible with radius r_c is filled with GeSi melt, and the height of the melt is noted as h. The radius of the growth crystal is marked as r_s. The radius ratio (η = r_s/r_c) is 0.3. The temperature and concentration at the sidewall (T_c, C_c) are higher than that along the crystal/fluid interface (T_s, C_s). The bottom is thermally adiabatic, and no-slip and impermeable boundary conditions are applied. The flat surface is not deformable. The flow is assumed to be laminar and the test fluid is incompressible. The linear variations of surface tension σ and density ρ satisfy the following conditions:

$$\sigma (T,C)={\sigma }_{\mathrm{r}}-{\gamma }_{\mathrm{T}}\left(T-T_{\mathrm{r}}\right)-{\gamma }_{\mathrm{C}}\left(C-C_{\mathrm{r}}\right)$$

(1)

$$\rho (T,C)={\rho }_{\mathrm{r}}\left[1-{\beta }_{\mathrm{T}}\left(T-T_{\mathrm{r}}\right)-{\beta }_{\mathrm{C}}\left(C-C_{\mathrm{r}}\right)\right]$$

(2)

where β_T = −(∂ρ/∂T)_C/ρ_r, β_C = −(∂ρ/∂C)_T/ρ_r are the thermal and solutal expansion coefficients. γ_T =− (∂σ/∂T)_C and γ_C = −(∂σ/∂C)_T are the temperature and concentration coefficients of surface tension, respectively. The subscript r represents the reference value.

2.2. Mathematical Formulation

With the scale quantities of r_c, ν/r_c, r_c²/ν, and ρν²/r_c², the control equations are transformed to the non-dimensional form:

$$\nabla \cdot V=0$$

(3)

$$\frac{\partial \mathit{V}}{\partial \tau }+\mathit{V}\cdot \nabla \mathit{V}=-\nabla P+{\nabla }^2\mathit{V}+G{r}_{\mathrm{T}}\left(\Theta +R_{\rho }\Phi \right){\mathit{e}}_Z$$

(4)

$$\frac{\partial \Theta }{\partial \tau }+\mathit{V}\cdot \nabla \Theta =\frac{1}{Pr}{\nabla }^2\Theta$$

(5)

$$\frac{\partial \Phi }{\partial \tau }+\mathit{V}\cdot \nabla \Phi =\frac{1}{Le\cdot Pr}{\nabla }^2\Phi$$

(6)

The non-dimensional velocity vector is expressed as V, the temperature is in the form of Θ = (T−T_s)/(T_c−T_s), and the concentration is Φ = (C−C_s)/(C_c−C_s). In addition, Gr_T = gβ_T(T_c–T_s)r_c³/ν² is the thermal Grashof number, and R_ρ is the buoyancy ratio, which is defined as R_ρ = β_C(C_c–C_s)/[β_T(T_c–T_s)]. Pr = ν/a denotes the Prandtl number and Le = α/D is the Lewis number, where ν is kinematic viscosity, α is thermal diffusivity, and D is mass diffusivity of the working fluid.

The initial conditions are as follows:

$$V_{\theta }=V_{\mathrm{R}}=V_{\mathrm{Z}}=0$$

(7)

$$\Theta =\Phi =0{,}_{{}_{}}R<R_{\mathrm{s}}$$

(8)

$$\Theta =\Phi =-\frac{\ln \left(R/R_{\mathrm{s}}\right)}{\ln \left(R_{\mathrm{c}}/R_{\mathrm{s}}\right)}{,}_{{}_{}}R\ge R_{\mathrm{s}}$$

(9)

The following are the established boundary conditions:

at the crystal/fluid interface:

$$V_{\theta }=V_{\mathrm{Z}}=0,$$

(10)

$$V_{\mathrm{R}}=R{e}_{\mathrm{s}}\frac{R}{R_{\mathrm{s}}}$$

(11)

$$\Theta =\Phi =0$$

(12)

along the free surface:

$$\frac{\partial V_{\mathrm{R}}}{\partial Z}=-R{e}_T\frac{\partial \mathsf{\Theta }}{\partial R}-R{e}_{\mathrm{CON}}\frac{\partial \mathsf{\Phi }}{\partial R}$$

(13)

$$\frac{\partial V_{\theta }}{\partial Z}=-R{e}_{\mathrm{T}}\frac{\partial \mathsf{\Theta }}{R\partial R}-R{e}_{\mathrm{CON}}\frac{\partial \mathsf{\Phi }}{R\partial R}$$

(14)

$$V_Z=0$$

(15)

$$\frac{\partial \Theta }{\partial Z}=\frac{\partial \Phi }{\partial Z}=0$$

(16)

at the bottom:

$$V_R=V_Z=0$$

(17)

$$V_{\mathsf{\theta }}=R{e}_{\mathrm{c}}R$$

(18)

$$\frac{\partial \Theta }{\partial Z}=\frac{\partial \Phi }{\partial Z}=0$$

(19)

along the crucible’s sidewall:

$$V_R=V_Z=0$$

(20)

$$V_{\theta }={Re}_{\mathrm{c}}^{}$$

(21)

$$\Theta =\Phi =1$$

(22)

To describe the effect of thermal and solutal capillary effects, the corresponding Reynolds numbers are introduced:

$$R{e}_{\mathrm{T}}=\frac{{\gamma }_{\mathrm{T}}\left(T_{\mathrm{c}}-T_{\mathrm{s}}\right)r_{\mathrm{c}}}{\mu v}$$

(23)

$$R{e}_{\mathrm{CON}}=\frac{{\gamma }_{\mathrm{T}}\left(C_{\mathrm{c}}-C_{\mathrm{s}}\right)r_{\mathrm{c}}}{\mu v}$$

(24)

The rotation Reynolds numbers are defined as

$$R{e}_{\mathrm{s}}=\frac{2\pi r_{\mathrm{s}}{n}_{\mathrm{s}}{r}_{\mathrm{c}}}{60v}$$

(25)

$$R{e}_{\mathrm{c}}=\frac{2\pi n_{\mathrm{c}}{r}_{\mathrm{c}}^2}{60v}$$

(26)

The relative strength of the soluto and thermocapillary forces is expressed by the capillary ratio:

$$R_{\sigma }=\frac{R{e}_{\mathrm{CON}}}{R{e}_{\mathrm{T}}}=\frac{{\gamma }_C\left(C_{\mathrm{c}}-C_{\mathrm{s}}\right)}{{\gamma }_T\left(T_{\mathrm{c}}-T_{\mathrm{s}}\right)}$$

(27)

The physical properties of the test fluid are listed in

Table 1. Physical properties of silicon-germanium mixture.

Property	Symbol	Unit	Value
Kinematic viscosity	ν	m2/s	1.4 × 10−7
Thermal diffusivity	α	m2/s	2.2 × 10−5
Mass diffusivity	D	m2/s	1.0 × 10−8
Temperature coefficient of surface tension	γT	N/(m·k)	8.1 × 10−5
Concentration coefficient of surface tension	γC	N/m	−0.54
Prandtl number	Pr	-	6.4 × 10−3
Lewis number	Le	-	2197.8

[42].

2.3. Numerical Procedure and Method Validation

The control equations and boundary conditions are handled with the finite volume method with non-uniform grids. The mesh refinement is considered near the wall and free surface. The modified central-difference method is utilized for solving the diffusion term, and the convection term is approximated by QUICK format. Simultaneously, the coupling numeration of pressure and velocity is solved using the SIMPLE algorithm. The non-dimensional time step is between 2 × 10⁻⁶ and 2 × 10⁻⁵. The criteria for convergence is defined as follows:

$${\left|R^{\mathrm{i}}\right|}_{\mathrm{m}\mathrm{a}\mathrm{x}}\le {10}^{-12}\mathrm{and}{\left|\frac{{\xi }^{\mathrm{i}+1}-{\xi }^{\mathrm{i}}}{{\xi }^{\mathrm{i}+1}}\right|}_{\max }\le {10}^{-4}$$

where ${\left|R^{\mathrm{i}}\right|}_{\max }$ represents the maximum absolute residual in the iteration and ξ represents the velocity, temperature, and concentration. Mesh dependence is tested and typical results at four levels of grids with different densities are listed in

Table 2. Mesh dependence of oscillatory flow at Re_T = 9938 and Re_c = 936 (a).

Mesh	m	F1	F2
60R × 20Z × 80θ	8	650	1165
80R × 35Z × 120θ	8	665	1177
100R × 35Z × 160θ	8	657	1171
120R × 55Z × 200θ	8	659	1167

. It is shown that the four different grids produce almost the same surface patterns, the wave numbers are the same, and maximum deviations of the non-dimensional oscillation frequency at a monitoring point are less than 2%. Therefore, the moderate mesh of 100^R× 35^Z× 160^θ used in this work is appropriate for accurate simulations.

To validate the numerical method, the rotation driven flow experimentally investigated by Kanda [44] was numerically simulated, and the same surface pattern was reproduced. Further, 3D simulations on the double-diffusive flow in a rectangular enclosure were conducted under the same conditions as in the report by Zhan et al. [45]. It was found that the oscillation frequency and Nusselt number agree well with the results by Zhan et al. These sufficiently indicate the effectiveness of the simulation method.

3. Results and Discussions

3.1. Basic Two-Dimensional Flow and Stability

In this work, the capillary ratio is minus one and the thermal and soluto-capillary forces are counterbalanced. However, under a normal gravity environment, the appearance of buoyancy force will inevitably lead to the disturbance of the balance along the free surface, generating the thermo-solutal-buoyancy convection. When the system rotations are considered, the centrifugal and Coriolis forces are induced, and the rotating driven flow is coupled with thermo-solutal-buoyancy convection. When the driving force is small, the basic mixed flow is in two-dimensional (2D) steady and axisymmetric state. Owing to the good fluidity and thermal uniformity, this steady state is preferable for the crystal growth. The non-dimensional stream function ψ is introduced to present the velocity distribution:

$$V_{\mathrm{R}}=-\frac{1}{R}\frac{\partial \psi }{\partial Z}{,}_{}{V}_{\mathrm{Z}}=\frac{1}{R}\frac{\partial \psi }{\partial R}$$

(28)

As the thermo-solutocapillary Reynolds numbers are related by the introduction of capillary ratio, the values of thermal Grashof number and thermocapillary Reynolds number are also dependent on the temperature difference; thus, for a given value of Re_T, both Re_CON and Gr_T can be determined. For brevity, only the values of Re_T are given in Figure 2, which shows the typical flow structures for the basic two-dimensional steady flow. The marked positive or negative sign indicates the clockwise or counter-clockwise circulation directions. For the case of a mixture with minus one capillary ratio, the soluto and thermal-capillary forces are equal, but act in contrary directions, and the overall capillary effect is zero. However, even without rotation, the existence of buoyancy breaks the equilibrium between the capillary forces. The fluid always moves from the crucible sidewall to the crystal, producing a counterclockwise vortex, as shown in Figure 2a. Because of the large Lewis number of the melt, the heat diffusion speed is faster than that of mass, and the isotherms are almost kept as conductive state. For the small thermocapillary Reynolds number of Re_T = 394, the iso-concentration lines bent slightly, and the flow rate is small on the free surface, as illustrated in Figure 3.

When the crystal starts to rotate at Re_s = −561, as shown in Figure 2b, a clockwise cell forms under the crystal and dominates the flow field. Compared with the natural convection displayed in Figure 2a, the mixed natural and forced flow is stronger. When rotation rate is increased to −1400, as shown in Figure 2c, the clockwise cell becomes stronger and the maximum stream function increases from 0.12 to 0.80. The enhanced crystal rotation induces a uniform concentration distribution underneath the crystal–fluid interface. When the crucible rotation is introduced, as illustrated Figure 2d, a primary counterclockwise rolling cell appears. The surface flow rate near the crystal is greatly increased. For the case of counter rotation of crucible and crystal, two counter-circulations are formed and rival each other. The competition of the vortex creates the zigzag profile of radial velocity. Away from the crystal/fluid interface, thermal and soluto-capillary forces are imposed oppositely and counterbalance each other, causing the surface fluid flows to be rather sluggish.

Once the values of Re_T or rotation rate exceed threshold values, the basic flow rapidly loses the stability and transforms into 3D oscillatory state. For the crystal growth, the oscillations should be avoided as much as possible. The amplitudes of oscillations are calculated and the dichotomy method is utilized to obtain the critical values of thermocapillary Reynolds numbers Re_T,cri. Figure 4 gives the variations of critical thermocapillary Reynolds numbers at different rotation rates. Without rotation, Re_T,cri approximately equals 1010; this value is larger than that of fluid mixture with a capillary ratio of −0.2, but still smaller than that of pure silicon melt [35]. This indicates that the introduction of the solutocapillary effect makes it easier for the flow to lose stability.

With the crystal rotation, the directions of the thermocapillary and centrifugal forces are opposite, which reduces the intensity of natural convection flow. Thus, with the rising Re_s, the value of Re_T,cri increases. However, when the value of Re_s reaches 1123, the forced convection driven by crystal rotation is predominant and the enhanced disturbance stimulates flow instability. As a result, Re_T,cri undergoes a sharp drop. Particularly, when Re_s is increased to 1680, no stable state is observed. This indicates the pure rotating flow is unstable. On the other hand, when the crucible starts to rotate, the large contact area allows the melt to rotate with the crucible synchronically, inducing a stabilization effect on the mixed convection. Therefore, Re_T,cri experiences a monotonic increase with the increase of Re_c.

3.2. Three-Dimensional Steady Flow

The spatial fluctuation δξ is extracted to characterize the 3D flow:

$$\delta \xi (R,\theta ,Z,\tau )=\xi (R,\theta ,Z,\tau )-\frac{1}{2\mathsf{\pi }}{\displaystyle {\int }_0^{2\mathsf{\pi }}\xi }(R,\theta ,Z,\tau )d\theta$$

(29)

where ξ represents Θ, Φ, or velocity V.

Figure 5 shows the concentration fluctuations and the corresponding spatial-temporal diagram (STD) for the case of Re_T = 3943 without rotation. It is observed that the surface fluctuations are presented as straight spokes, and the STD is composed of four vertical strips. These surface patterns indicate the fluctuations are stationary in time, but oscillate in space. For the further understanding of the fluctuations, the streamlines are plotted in cross sections of θ = 0 and θ = π/4 (which correspond to the dark and bright strips, respectively, shown in Figure 5). As illustrated in Figure 6a, influenced by the thermocapillary effect, the fluid near the sidewall is pushed inward, forming a counterclockwise rolling cell. This convective vortex stirs the melt, inducing a local uniform concentration distribution near the crystal. Meanwhile, it can be seen that, away from the crystal, the deformation of iso-concentration lines reflects the existence of a reversed vertical concentration gradient. Therefore, the solutal buoyancy force is generated and the lighter fluid in the lower part is brought to the surface, forming a clockwise cell. At the cross section of θ = π/4, the rolling cell generated by the thermocapillary effect is squeezed by the small vortex near the crystal.

The surface fluctuations of the mixed flow with crystal rotation are shown in Figure 7. When Re_T = 3943 and Re_s = −561, the surface pattern is still shown as a spoke pattern, as seen in Figure 7a. However, an annular band is observed around the crystal. The concentration fluctuations in the annular band are almost zero. Compared with Figure 5, the amplitude of the surface fluctuations is decreased from 0.24 to 0.038, but the wave number increases greatly. The natural convection still dominates the flow field. The corresponding local capillary ratio $R_{\sigma }^{\mathrm{i}}$ on the free surface and the flow structures in the cross section are plotted in Figure 8 and Figure 9a, respectively. Near the crystal, the combination of rotation and solutocapillary effects drives the mixture with a low concentration to move outside, leading to a uniform concentration gradient in the range of 0.3 < R < 0.6. When 0.6 < R < 0.82, the sparse iso-concentration lines are observed, indicating the small concentration gradient. Meanwhile, the absolute value of $R_{\sigma }^{\mathrm{i}}$ is smaller than unity, which means the thermocapillary effect is stronger than the solutocapillary influence, thus the counterclockwise cell is dominant in this area. Near the sidewall, the combination of solutocapillary and buoyancy forces generates a weak clockwise circulation in the top right corner of the crucible. With the increase of Res, as shown in Figure 7b, the occupied area of the annular band is expanded and the surface fluctuations near the sidewall are presented as a rosebud pattern. The centrifugal force further pushes the fluid forward and squeezes the rolling cell driven by the thermocapillary effect, as described in Figure 9b. For the case of Re_T = 3943, Re_s = −1123, the surface fluctuations are shown as a pulsating ring, the oscillation amplitude is greatly deduced, and the wave number is lowered to 4, as illustrated in Figure 7c. From the variation of the local capillary ratio (Figure 8), it is shown that the region where $\left|R_{\sigma }^{\mathrm{i}}\right|$ is smaller than unity coincides with the position of the surface circular ring. This confirms that the thermocapillary force is dominant in the pulsating ring area.

3.3. Three-Dimensional Oscillatory Flow

With the increasing temperature and concentration gradients, the steady flow transits to the 3D unsteady state. This unsteady flow not only oscillates in space, but also fluctuates with time. When the crystal rotation rate is kept at Re_s = −1123, once Re_T is greater than 7098, the rotating wave appears on the free surface. For a typical case of Re_s = −1123 and Re_T = 7887, the surface fluctuations and the corresponding STD at different radial positions are displayed in Figure 10. It is noted that the surface pattern can be regarded as a combination of spokes and rotating curved arms. Near the crystal (R = 0.32), the STD is shown as oblique lines leaning to the left. This means that the curved arms move in a clockwise direction and the rotation rate is less than that of the crystal. At R = 0.75, the corresponding the STD is displayed as vertical lines. This verifies that the spokes can be considered as standing waves. Between the rotating and standing waves, the annular band is still observed; the concentration distribution at a monitoring point (R = 0.52) inside the annular space is exhibited in Figure 11a. Along the circumference, the concentration is almost uniform and no fluctuation is detected in the annular band. The characteristics of the time-dependent of concentration and azimuthal velocity are also analyzed in Figure 11b. It is noted that the concentration oscillation lags behind velocity. The hysteresis between the oscillations was also observed in the report by Smith and Davis [46], which is one of the typical characteristics of hydrothermal waves generated by temperature difference. Moreover, the propagating angle of 17° is another feature of the hydrothermal waves. As shown in Figure 10, the propagating angle of the concentration oscillations is about 16°. These two values of wave angles are very close. Therefore, the spiral concentration oscillation wave is considered as hydrosoultal waves. Meanwhile, the oscillation frequencies are obtained by fast Fourier transform (FFT), as shown in Figure 12, and the main dimensionless frequency F₀ is 1427 with a secondary frequency F₁ of 2853, which satisfies F₀ = 1/2F₁. This frequency relationship was also experimentally reported by Shen et al. [34].

When the effect of crucible rotation is considered, the evolutions of surface oscillations are presented in Figure 13. The standing spokes are replaced by the rotating waves. For these cases, the counteraction of the centrifugal solutocapillary effect produces an annular space surrounding the crystal. With the increase of Re_T, the natural convection is enhanced, and the oscillation amplitude is enlarged. The wave number also has a tendency to increase, as presented in Figure 13b, and the waves with weaker oscillations near the annular band appear. When Re_T is increased to 15,773, as shown in Figure 13a, the grown flow instabilities lead to the increase of the wave number in the azimuthal direction to dissipate energy. The FFT analysis corresponding to the three typical cases is displayed. It is noted that a multiple relationship is observed among the frequencies, which is f₁ = 1/2f₂ = 1/3f₃.

For the case of Re_T = 5442, Re_c = 1871, the forced flow driven by crucible rotation plays a decisive role in the occurrence of instability. From the STD shown in Figure 14a, it can be observed that synchronous rotation is achieved between the crucible and oscillation waves. The wave number is 4. Near the crystal, the absolute value of $\left|R_{\sigma }^{\mathrm{i}}\right|$ is less then unity, as displayed in Figure 15, the thermocapillary and centrifugal forces are superimposed to enhance the rotating forced flow. Thus, a strong counterclockwise circulation is formed and occupies the left-half part of the crucible, as illustrated in Figure 16a. When R is greater than 0.6, the local capillary ratio oscillates around its equilibrium value of −1. Near the crucible sidewall, the solutal-buoyancy effect drives a clockwise circulation at the corner. When the Re_T is increased to 19,717, the solutocapillary effect is enhanced, and $\left|R_{\sigma }^{\mathrm{i}}\right|$ is greater than one in the region of R < 0.5. As the solutocapillary force counteracts the centrifugal force, the strength of the rolling cell is weakened near the crystal, and a series of surface spiral waves appear and extrude the rotating wave outward. Through the comparison between Figure 14a,b, it is found that, with the increase of Re_T, the rotating wave number increases twice, and the oscillation amplitude is reduced by half. This reveals that the enhanced natural convection has a depression effect on the mixed 3D oscillation flow.

4. Conclusions

A series of 3D numerical simulations were performed to investigate the mixed forced convection driven by rotation and natural convection generated by concentration and temperature gradients in the Czochralski crystal growth model. The main conclusions are as follows:

(a)

For the mixture with capillary ratio of minus one, the total thermal and solutocapillary forces are counterbalanced. The introduction of buoyancy force leads to the disturbance of the balance, generating the mixed capillary-buoyancy convection. For a small concentration gradient, the flow is in two-dimensional steady state. Owing to the good fluidity and thermal uniformity, this state is important for the growth of high quality of crystals. The coupled capillary and buoyancy flow in the crucible is presented as a large counterclockwise circulation. When the rotations of the crystal and crucible are considered, the mixed natural and forced flow structures are more complex, and the directions of the rolling cells are associated with the competitions among the driving forces.

(b)

When the capillary force is greater than a certain value, the basic flow transits to three-dimensional state. The critical conditions for the mixed flow transitions at different rotation rates are obtained. Crucible rotation can obviously strengthen the flow stability. Influenced by the crystal rotation, the critical thermocapillary Reynolds number increases until it reaches a turning point. With the enhancement of crystal rotation driven flow, a dramatic decrease of the critical value is observed.

(c)

Once the instability is incubated, the basic mixed flow firstly bifurcates to the three-dimensional steady state, which oscillates spatially along the circumferential direction. Driven by the competition among the capillary-buoyancy forces, centrifugal, and Coriolis forces, the surface fluctuations are presented as spokes, rosebud, and pulsating ring. With the enhanced flow instabilities, three-dimensional unsteady oscillation occurs. Prosperous oscillation patterns are discussed, including the spiral hydrosoultal waves, superimposition of spirals and spokes, as well as rotating waves.