Numerical Simulation of Thermal-Solutal Capillary-Buoyancy Flow of Ge_1–xSi_x Single Crystals Driven by Surface-Tension and Rotation in a Czochralski Configuration

Abstract

A series of three-dimensional numerical simulations were performed to understand the thermal-solutal capillary-buoyancy flow of Ge_1-xSi_x melts during Czochralski crystal growth with a rotating crystal or crucible. The crystal and crucible rotation Reynolds numbers in this work are 0∼3.5 × 10³ (0∼4.4 rpm) and 0∼−2.4 × 10³ (0∼−1.5 rpm), respectively. Simulation results show that if the thermal capillary Reynolds number is relatively low, the flow will be steady and axisymmetric, even though the crystal or crucible rotates at a constant rate. The critical thermal capillary Reynolds number for the initiation of the three-dimensional oscillatory flow is larger than that of pure fluids. As the crystal or crucible rotation rate increases, the critical thermal capillary Reynolds number first increases and then decreases. The dominant flow pattern after the flow destabilization is azimuthal traveling waves. Furthermore, a reversed evolution from the oscillatory spoke pattern to traveling waves appears in the melt. Once the crystal or crucible rotation rate is relatively large, the traveling waves respectively evolve to rotating waves at the crystal rotation and a spindle-like pattern at the crucible rotation. In addition, the maximum amplitude of solute concentration oscillation on the free surface initially decreases, but finally rises with the crystal or crucible rotation rate increasing.

1. Introduction

Single crystals of oxides and semiconductors mostly grow from melt, which takes on an important role in microelectronic and optoelectronic devices [1]. More than 80% of single crystals are grown by the Czochralski (Cz) method [2]. This method is particularly successful for the growth of germanium (Ge), silicon (Si), gallium-arsenide (Ga_1–xAs_x), germanium-silicon (Ge_1–x Si_x) and so on. Here, x is the mass fraction of arsenide or silicon.

A classical Cz crystal growth system consists of four important parts, crucible, crystal, melt, and heater. The crucible is heated by a heater to keep the materials in molten state. The crystal is kept in a lower temperature and placed so as to contact the surface of the melt. The crystal and crucible are usually rotated to improve the uniformity of thermal distribution and get a high quality of single crystals. Several forces, including surface tension, buoyancy, Coriolis, and centrifugal forces interact with each other on different scales through the crystal growth process, which makes the process difficult to characterize and control [3,4]. Due to the interaction between aforementioned forces, different flow instabilities always occur and significantly affect the transport process in the melt [3,4,5,6]. This transport process, in turn, considerably affects the quality of single crystals. Therefore, there is an urgent demand to fully understand flow phenomena in melts during the Cz crystal growth process.

The thermal capillary effect generated by the temperature difference between the crystal and the crucible during the Cz crystal growth process, as well as the buoyancy, have great effects on the flow field. Galazka and Wilke [7] numerically predicted that the coupling action of thermal capillary and buoyancy forces leads to the interface being much more convex to the melt, compared to the results without the thermal capillary flow. Schwabe [8] carried out experiments on the thermal capillary flow in a Cz configuration with various gravities. The thermal capillary flow penetrates fully into the liquid pool under conditions of microgravity, whereas it separates into a layer-like flow close to the free surface on the earth. Li et al. [9,10] studied the flow bifurcation of the thermal capillary flow in a shallow molten silicon pool. They observed that the steady two-dimensional (2D) flow bifurcates to a steady three-dimensional (3D) flow, and then evolves to a 3D oscillatory flow as the thermal capillary Reynolds number rises. However, a steady 2D flow is directly invaded by a 3D oscillatory flow in a deep pool filled with molten silicon. Subsequently, Shen et al. [11,12] observed the flow patterns and the transitions of the thermocapillary-buoyancy flow of a pure fluid in a Cz configuration with rotating crystal and crucible.

The solute/impurity concentration shows an uneven distribution, generated by segregation during the Cz crystal growth process, which induces the solutal capillary force and the solutal buoyancy force in the melts [13,14,15]. It is obvious that the coupling effect of various driving forces makes the melt flow become much more complex. Bergman [16] first reported that, even if the solutal and thermal capillary forces are equal but opposite due to horizontal solute concentration and temperature gradients, the thermal-solutal capillary flow can be observed in the absence of a buoyancy force, because of the difference between mass and thermal diffusivities. Soon after, Chen et al. [17,18] systematically investigated flow instability with different Prandtl (Pr) numbers, Lewis (Le) numbers, and aspect ratios by 3D numerical simulations. Kamotani et al. [19] reported complex flow patterns and the bifurcations of the thermal-solutal capillary-buoyancy convection at different buoyancy ratios. If the buoyancy ratio is relatively small, a unicellular flow pattern, accompanied by secondary cells, was observed. However, a three-layered pattern with secondary cells appears with a large buoyancy ratio. Arafune et al. [20,21] observed that the classical surface velocities in solutal capillary flows are about 3–5 times larger than those in thermal capillary flows.

Recently, thermal-solutal capillary-buoyancy flows in annular pools have attracted much more attention. The annular pool is a classic simplified model of a Cz configuration. Firstly, Li et al. [22] conducted 2D simulations to study the thermal-solutal capillary-buoyancy flow of binary mixtures. They found that the development from a steady flow to an oscillatory flow experiences a typical Hopf bifurcation. Due to the limit of a 2D geometric model in the azimuthal direction, Yu et al. [23,24,25] pursued a similar direction by conducting experiments and 3D numerical simulations. They reported that the solutal capillary and buoyancy forces, which are generated by the solute gradient as a result of the Soret effect, promote the destabilization of the thermal-solutal capillary-buoyancy flow in binary mixtures. However, they have almost no impact on flow patterns and their transition. It is surprising that the thermal-solutal capillary-buoyancy flow can always be observed, even if the solutal and thermal capillary effects are equal but opposite. The critical thermal Reynolds number for the generation of a 3D flow decreases with the aspect ratio of the annular pool. However, it increases alongside the increase of the capillary ratio. Flow patterns strongly depend on the capillary ratio, aspect ratio, and thermal capillary Reynolds number.

Although annular pools can fulfill fundamental requirements for the study of flow instability during the crystal growth process, it still has some limitations, such as the lack of the crystal face touching the melt free surface, its inability to realize crystal and crucible rotations and so on. Therefore, it is necessary to explore thermal-solutal capillary-buoyancy flow in a Cz configuration. Abbasoglu and Sezai [26] reported that the approximate range of capillary ratios in a Ge_1-xSi_x melt, due to the segregation during Cz crystal growth, varies from −2.44 to −0.20. In addition, the radial solute concentration difference decreases along the rising of the aspect ratio and radial temperature difference. This is because the strong flow suppresses the solute concentration difference. Unfortunately, the rotations of crystals and crucibles were not taken into consideration, which is an essential industrial practice for improving the quality of crystals in the Cz crystal growth process. A lot of previous works focusing on the melt flow and the flow instability have not yet paid attention to solute diffusion in binary alloys [27,28,29,30].

Therefore, it is indispensable to carry out numerical simulations to precisely understand the influences of rotating crystals and crucibles on the thermal-solutal capillary-buoyancy convections in a Cz configuration. The present work is beneficial for improving crystal qualities in Cz crystal growth.

2. Physical and Mathematical Model

2.1. Basic Assumptions and Governing Equations

Figure 1 gives the physical model, which is a classical simplified model of the Cz growth system [31,32,33]. A disc with the radius of r_s = 23 mm contacts the free surface of Ge_1–xSi_x melt. The Ge_1–xSi_x melt is in a cylindrical crucible with the radius of r_c = 46 mm (hereafter are respectively called as “crystal” and “crucible”). The crystal and crucible rotate at constant rates n_s and n_c. The melt/crystal interface and sidewall of the crucible are maintained at constant temperatures and solute concentrations are T_s, C_s, T_c, and C_c, where T_s = T_m < T_c and C_s < C_c. T_m = 1271.3 K, is the melting point of the Ge_1–xSi_x melt. The depth of the Ge_1–xSi_x melt is h. Furthermore, the Ge_1–xSi_x melt is considered as an incompressible Newtonian fluid with constant properties, except for the surface tension and density in the buoyancy term of the momentum equation. The flow is laminar. The melt/crystal interface and free surface are non-deformable and flat, where the solute and thermal capillary forces are considered. A no-slip condition is adopted in all solid-liquid boundaries. In addition, the free surface and crucible bottom are assumed to be impermeable and adiabatic.

The surface tension and density are linear functions of temperature and solute concentration, respectively, which are defined as

$$\sigma \left(T,C\right)={\sigma }_0-{\gamma }_{\mathrm{T}}\left(T-T_0\right)-{\gamma }_{\mathrm{C}}\left(C-C_0\right)$$

(1)

$$\rho \left(T,C\right)={\rho }_0\left[1-{\beta }_{\mathrm{T}}\left(T-T_0\right)-{\beta }_{\mathrm{C}}\left(C-C_0\right)\right]$$

(2)

where γ_T = −(∂σ/∂T)_C,γ_C = −(∂σ/∂C)_T, β_T = −(∂ρ/∂T)_C,β_C = −(∂ρ/∂C)_T. Here, T₀ and C₀ are the reference temperature and solute concentration, respectively.

Using r_c, r_c²/ν, ν/r_c, and μν /r_c² as the scale quantities for length, time, velocity, and pressure, the governing equations are expressed by following dimensionless forms,

$$\nabla \cdot \mathit{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}}(\Theta +R_{\rho }\Phi ){\mathit{e}}_{\mathit{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)

where V (U, V, W) is the dimensionless velocity vector and e_Z is the unit vector in the Z direction. P, τ, Θ = (T − T_s)/(T_c−T_s), and Φ = (C − C_s)/(C_c − C_s) are the dimensionless dynamic pressure, time, temperature, and solute concentration, respectively. R_ρ = β_CΔC/(β_TΔT), Gr_T = gβ_TΔT(r_o-r_i)³/ν², Pr = ν/α, and Le = α/D are the ratio of the solute to thermal buoyancy forces, the thermal Grashoff number, Prandtl number, and Lewis number. Where ΔT = T_c − T_s and ΔC = C_c − C_s. Thus, the last term in Equation (4) shows the buoyancy effect. In addition, D, a, ν, and μ are the mass diffusivity of the silicon, thermal diffusivity, kinematic viscosity, and dynamic viscosity.

Based on these assumptions, the boundary conditions can be expressed as follows.

At the free surface

$$R{\mathrm{s}}_{}<R<R_{\mathrm{c}},0<\theta \le 2\pi ,Z=Г$$

(7)

$$W=0,\frac{\partial U}{\partial Z}=-R{e}_{\mathrm{T}}\frac{\partial \Theta }{\partial R}-R{e}_{\mathrm{C}}\frac{\partial \Phi }{\partial R}$$

(8)

$$\frac{\partial V}{\partial Z}=-R{e}_{\mathrm{T}}\frac{\partial \Theta }{R\partial \theta }-{Re}_{\mathrm{C}}\frac{\partial \Phi }{R\partial \theta }$$

(9)

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

(10)

Here, Г is the aspect ratio, which is defined as Г = H/R_c. Re_T and Re_c are the thermal and solutal capillary Reynolds numbers. They are respectively defined as

$$R{e}_{\mathrm{T}}=\frac{{\gamma }_T{r}_{\mathbf{c}}(T_c-T_s)}{\nu \mu },R{e}_{\mathrm{C}}=\frac{{\gamma }_C{r}_{\mathbf{c}}(C_{\mathrm{c}}-C_{\mathrm{s}})}{\nu \mu }$$

(11)

The ratio between the solutal and thermal capillary effects is denoted as the capillary ratio R_σ, which is defined as:

$$R\sigma =\frac{R{e}_{\mathrm{C}}}{R{e}_{\mathrm{T}}}=\frac{{\gamma }_{\mathrm{C}}\Delta C}{{\gamma }_{\mathrm{T}}\Delta T}$$

(12)

At the bottom (0 < R < R_c, $0<\theta \le 2\pi$, Z = 0)

$$U=W=0,V_{}={Re}_{\mathrm{c}}\cdot R,\frac{\partial \Theta }{\partial Z}=\frac{\partial \Phi }{\partial Z}=0$$

(13)

Re_c is the crucible rotation Reynolds number, which is expressed as

$$R{e}_{\mathrm{c}}=\frac{2\mathsf{\pi }{n}_{\mathrm{c}}{r}_{\mathrm{c}}{}^2}{60\nu }$$

(14)

At the melt/crystal interface (R ≤ R_s, $0<\theta \le 2\pi$, Z = Г)

$$U=W=0,V={Re}_{\mathrm{s}}\cdot \frac{R}{R_s},\Theta =\Phi =0$$

(15)

Here, Re_s is the crystal rotation Reynolds number, which is defined as

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

(16)

At the sidewall of the crucible (R = R_c, $0\le \theta \le 2\pi$,$0\le Z\le \Gamma$)

$$U=W=0,V={Re}_{\mathrm{c}},\Theta =\Phi =1$$

(17)

The initial conditions ($\tau =0$)

$$U=W=V=0$$

(18)

$$\Theta =\Phi =0\left(R\le R_{\mathrm{s}}\right)$$

(19)

$$\Theta =\Phi =\ln \left(R/R_{\mathrm{s}}\right)/\ln \left(R_{\mathrm{c}}/R_{\mathrm{s}}\right),(R_{\mathrm{s}}<R<1)$$

(20)

2.2. Calculations’ Conditions and Numerical Method

The radius ratio and aspect ratio of the Cz configuration are fixed at R_s = Г = 0.5. The Ge_1−xSi_x melt is chosen as the working fluid. Its initial mass fraction of silicon is 1.99% (x = 1.99%), which corresponds to the molar fraction of 5%. This mass fraction agrees well with the mean solute concentration in the melt, due to segregation during Cz crystal growth [26,34]. The physical properties of the Ge_1–xSi_x melt at 1271.3 K are shown in

Table 1. Physical parameters of the Ge_0.98Si_0.02 melt.

Property	Symbol	Unit	Value
Density	ρ	kg/m3	5246.00
Thermal diffusivity	α	m2/s	2.20 × 10−5
Viscosity	μ	kg/(m·s)	7.34 × 10−4
Mass diffusivity of species	D	m2/s	1.00 × 10−8
Temperature coefficient of surface tension	γT	N/(m·k)	8.10 × 10−5
Concentration coefficient of surface tension	γC	N/m	−0.536
Prandtl number	Pr	-	6.37 × 10−3
Lewis number	Le	-	2197.80

[26,34]. As reported in Reference [26], the approximate range of capillary ratios in the Ge_xSi_1–x melt due to segregation during Cz crystal growth is from −2.44 to −0.20. Therefore, R_σ is fixed at −0.8, which ensures that much more attention is paid to the influence of the rotating crystal and crucible on the thermal-solutal capillary-buoyancy flow. The crystal and crucible rotation Reynolds numbers range from 0 to 3.5 × 10³ and −2.4 × 10³. The thermal capillary Reynolds number ranges from 0 to 4 × 10⁴.

The governing equations are handled by the finite volume method (FVM). The QUICK scheme is applied for the convective terms. The second order central-difference approximation is introduced for the diffusion terms. Meanwhile, the SIMPLEC algorithm is adopted to deal with the coupling between the velocity and pressure. The dimensionless time step is from 1.2 × 10⁻⁵ to 1.1 × 10⁻⁴. When the largest relative error of these governing equations in the computational domain gets below 10⁻⁵ at each time step, the solution is taken for converging.

In order to verify the grid convergence, simulations on the thermal-solutal capillary-buoyancy flow of Ge_1–xSi_x melt with three different meshes at Re_c = 9.4 × 10², Re_s = −1.2 × 10³, and Re_T = 2 × 10⁴ were conducted.

Table 2. Mesh dependence tests at Re_c = 9.4 × 10², Re_s = −1.2×10³_, and Re_T = 2 × 10⁴.

Grids	m	f (Hz)
42R × 40Z × 60θ	4	0.1879
62R × 50Z × 80θ	4	0.1887
80R × 60Z × 104θ	4	0.1848

gives the corresponding wave numbers and oscillation frequencies. It is obvious that the 62^R × 50^Z × 80^θ grid is a suitable choice for accurate simulations in this work.

In order to validate the numerical method, simulations on the double-diffusive capillary convection in rectangular cavities were conducted at the same conditions with Refs. [17,35]. Results show that the maximum deviations of the Nusselt number and oscillation frequency are 1.1% [35] and 1.5% [17], respectively, as listed in

Table 3. Comparisons of the Nusselt numbers at the left vertical wall.

ReT	Nu		Deviation (%)
	Present	Reference [35]
10	1.010	1.007	0.30
120	1.332	1.328	0.30
200	1.496	1.497	0.07
340	1.770	1.762	0.45
380	1.832	1.821	0.60
450	1.933	1.953	1.02

and

Table 4. Comparisons of the oscillation frequencies.

ReT	f/Hz		Deviation (%)
	Present	Reference [17]
350	6.94	6.84	1.46
410	7.89	7.82	0.90
500	9.35	9.26	0.97
540	10.14	10.02	1.20

. In addition, the present numerical method was also checked by experiments. The results from the numerical simulations agree well with those from the experiments [23]. Therefore, the thermal-solutal capillary-buoyancy flow can be solved well by the present numerical method.

3. Results and Discussion

3.1. Basic Flow

When the thermal capillary Reynolds number is relatively low, the flow is axisymmetric and steady, even though the crystal and crucible rotate at constant values. The steady axisymmetric flow is an ideal flow condition for the Cz growth process, due to its thermal and crystal uniformities. This typical flow was also reported in binary mixtures or pure fluids in a crystal or crucible rotating system, which is usually called “basic flow” [36,37,38]. The dimensionless stream function ψ is applied to show the flow field and defined as

$$U=-\frac{1}{R}\frac{\partial \psi }{\partial Z},{{\displaystyle }}_{}W=\frac{1}{R}\frac{\partial \psi }{\partial R}$$

(21)

Figure 2 shows the streamlines and iso-concentration lines at Re_T = 2 × 10³ with different rotation Reynolds numbers. When Re_c = Re_s = 0, the flow is only moved by the buoyancy and capillary forces. When the capillary ratio is fixed at −0.8, the thermal capillary force is dominant, which induces flow movement from the crucible sidewall to the crystal edge along the free surface, as displayed in Figure 2a. The iso-concentration lines mainly concentrate near the crucible sidewall and below the crystal. Once the crucible rotates, for example, at Re_c = 46.8 and Re_s = 0, the flow is slightly suppressed. The maximum stream function of the basic flow decreases from 16.4 to 16.1. However, the streamlines and iso-concentration lines have slight changes, as shown in Figure 2b. When only the crystal rotates, at Re_c = 0 and Re_s = −467.5, the melt near the crystal rotates along with the crystal. In addition, the maximum stream function decreases from 16.4 to 15.8 (see Figure 2c). It implies that the flow strength is slightly weakened by the rotation.

3.2. Critical Conditions for the Flow Destabilization

As the thermal capillary Reynolds number increases, the basic flow is clearly enhanced. Once the thermal capillary Reynolds number is over a threshold value, the basic flow rapidly destabilizes and then bifurcates into a 3D oscillatory flow. Figure 3 gives the variations of critical thermal capillary Reynolds numbers along the crystal or crucible rotation Reynolds number. This dichotomy is successfully applied for obtaining the critical thermal capillary Reynolds numbers [24,25], and maximum deviation is controlled below 100. The critical thermal capillary Reynolds number at a non-rotation Cz configuration is 4.1 × 10³, which is much larger than that of pure fluids [9,10]. This discrepancy is because of the opposite directions of the thermal and solutal capillary forces at R_σ = −0.8 [25]. Thus, the quality of a pure component crystal is much more essential to the temperature difference between the crystal and crucible than that of a compound.

The critical thermal capillary Reynolds number rises first, and then reduces along the crystal or crucible rotation rate. When R_σ = −0.8, the thermal capillary effect plays a dominant role on the flow. The melt on the free surface moves from the crucible sidewall to crystal edge. Once the crystal rotates only at a relatively low rate, the centrifugal force suppresses the surface melt flow. Therefore, the crystal rotation enhances the flow stability until |Re_s| = 470, which is also reported in pure fluids [9,39] or binary mixtures [37]. With the crystal rotation rate increasing, a clockwise rotating cell under the crystal is formed for the crystal rotation, which shears with a counter-clockwise rotating cell, driven by the coupling action between the thermal and solutal capillary effects. Thus, the fluctuation of the melt flow is enhanced, while the flow stability is weakened, as displayed in Figure 3.

At Re_c ≤ 230 and Re_s = 0, almost all melt moves azimuthally with the crucible, which leads to the melt flow being much more uniform and more stable [36]. As the crucible rotation Reynolds number rises, the disturbance formed by the rotation is enhanced, which promotes the steady flow towards bifurcating to a 3D oscillatory flow. Thus, the critical thermal capillary Reynolds number initially grows, and then decreases, alongside the crucible rotation Reynolds number rising. When the crystal or crucible rotation rate is smaller than a certain value, the rotation benefits the quality of a crystal. However, the crystal or crucible rotation can have negative effects on the quality of a crystal at a larger rotation rate. Once the crystal and crucible rotation Reynolds numbers respectively exceed about 1.2 × 10³ and 1.4 × 10³, the 3D oscillatory flow is dominated by the crystal or crucible rotation. Furthermore, the flow is almost independent of the thermal Reynolds number [11,36].

3.3. Characteristics of the 3D Oscillatory Flow

The fluctuation δX of physical quantity X is usually used to extract the spatial disturbance characteristics,

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

(22)

It should be pointed out that X is the dimensionless velocity, solute concentration, or temperature. Since the Schmidt number is much greater than the Prandtl number of a Ge_1–xSi_x melt, the solute concentration fluctuation attracts more attention here. Figure 4 gives the overview of surface solute concentration fluctuations of 3D oscillatory flows in a non-rotation Cz system. The flow pattern in Figure 4 looks like spokes of a wheel. It was named the “spoke pattern” in Reference [40]. The spoke pattern, with three traveling waves along a clockwise direction, is observed at Re_T = 5 × 10³, as shown in Figure 4a. The wave number m is obtained from the bright and dark spokes. This flow pattern does not appear in pure fluids [9]. In pure fluids, the thermal capillary and buoyancy are the only driving forces. However, in a Ge_1–xSi_x melt, there is a nonuniform solute concentration distribution along the axial direction, close to the free surface. This unique solute concentration distribution in Ge_1–xSi_x melts will result in the Rayleigh-Marangoni-Benard instability.

The wave number of spoke pattern increases as the thermal capillary Reynolds number increases. For example, once the thermal capillary Reynolds number grows from 5 × 10³ to 10⁴, the wave number increases from 3 to 9. Furthermore, the azimuthal traveling spoke pattern transits to an oscillatory spoke pattern. The distribution of spokes on the free surface becomes nonuniform. The enhanced buoyancy-driven flow, due to the nonuniform solute concentration distribution in the axial direction, carries a lot of Ge_1–xSi_x melt from the crucible bottom to the free surface, close to the sidewall at a large thermal capillary Reynolds number. Thus, the oscillatory spoke pattern swings in the azimuthal direction and merges, which is reported in Reference [24]. On the other hand, the amplitude of the surface solute concentration fluctuation decreases slightly, because of the increase of the wave number.

3.3.1. The Influences of Crystal Rotation

Figure 5 shows the effects of the crystal rotation on 3D oscillatory flow at Re_c = 0 and Re_T = 10⁴. When the crystal rotates at Re_s = −4.7 × 10², the oscillatory spoke pattern (see Figure 4b) develops to the azimuthal traveling spoke pattern with three waves, as shown in Figure 5a. The traveling waves at Re_c = 0, Re_s = −4.7 × 10² and Re_T = 10⁴ are the same, with the flow pattern at Re_T = 5 × 10³ in a non-rotation Cz system. As reported in Ref. [25], the capillary effect is a key factor of flow destabilization for the binary mixtures at −0.8 ≤ R_σ ≤ −0.2. Compared with the flow transition in a non-rotation Cz system, a reversed evolution from the oscillatory spoke pattern to traveling waves is observed when the crystal rotation rate increases. Thus, the crystal at a relatively low rotation rate plays a suppression effect on the flow instability. When Re_s = −1.8 × 10³, the traveling waves are disturbed by the rotating waves, which is induced by centrifugal force instability [11,12], as displayed in Figure 5b. The traveling waves still take a dominant role, while the rotating waves will affect the traveling waves. The wave number reduces from 3 to 2. Meanwhile, the traveling waves begin to exhibit disorder, especially in the region close to the crucible sidewall. The rotating waves are further enhanced with the increased crystal rotation rate. The rotating waves are observed near the crucible sidewall (see Figure 5c).

Figure 6 gives the time dependence of azimuthal velocity at the monitoring point P (R = 0.75, Z = 0.5, θ = 0). At Re_c = 0, Re_s = −1.8 × 10³, and Re_T = 10⁴, the azimuthal velocity appears as quasi-periodical variation with a period of 0.1, as shown in Figure 6a, which means that the traveling wave shows a quasi-periodic oscillatory flow. If the crystal rotation Reynolds number increases to 2.3 × 10³, the average amplitude of azimuthal velocity rapidly increases, while its periodicity diminishes, as displayed in Figure 6b. When the crystal rotation Reynolds number rises from 1.8 × 10³ to 2.3 × 10³, the rotating waves are enhanced significantly, as shown in Figure 5b,c. Meanwhile, the traveling waves lose the dominant position and then compete with the rotating waves. The flow loses its periodicity, since the competition between the traveling waves and rotating waves becomes fiercer alongside the crystal rotation Reynolds number increasing. Thus, the maximum azimuthal velocity is not observed at the crystal rotation Reynolds number of 2.3 × 10³ (see Figure 6b). Figure 7 shows that the maximum amplitude of solute concentration fluctuation initially decreases, but finally increases as the crystal rotation rate grows. If the crystal rotation rate is relatively low, the capillary effect plays a leading role. Meanwhile, the crystal rotation suppresses the flow instability generated by the capillary effect. Thus, the amplitude of solute concentration fluctuation reduces with the rise of the crystal rotation rate. The rotating waves are dominant when the crystal rotates at a large rate. The centrifugal force instability related to the rotation results in the increase of the amplitude of the solute concentration fluctuation along the crystal rotation.

3.3.2. The Influence of Crucible Rotation

Once the crucible rotates, the oscillatory spoke pattern also converts to traveling waves, as shown in Figure 8. Here, the traveling waves are mainly located close to the crystal. This flow bifurcation with the crucible rotation is same as that with the crystal rotation (compare Figure 5 and Figure 8). Thus, the crucible rotation with a low rotation rate also has a suppression effect on flow instability. Figure 9 gives the streamlines and iso-concentration lines of 3D oscillation flow at the R-Z plane of θ = 0 when Re_s = 0, Re_c = 9.4 × 10² and Re_T = 1.4 × 10⁴. Two counter-clockwise rotating cells are observed. These two rotating cells are the result of the coupling action of the thermal and solutal capillary effects. This has also been reported in non-rotating annular pools [25]. Due to the suppression role of the crucible rotation, the other one close to the crucible sidewall is weakened. Therefore, the oscillations of the solute concentration close to the crucible sidewall are much smaller than those near the crystal. The wave number m rises from 3 to 4 as the crucible rotation Reynolds number is raised from 9.4 × 10² to 1.4 × 10³. In addition, the traveling waves gradually extend to the sidewall. If the crucible rotates at a large rate, the traveling waves develop a spindle-like pattern, as displayed in Figure 8c. As reported in Refs. [11,12], the spindle-like pattern should be a result of the circular shear instability.

Figure 10 gives the variations of the maximum amplitude of solute concentration fluctuation along the crucible rotation Reynolds number. Obviously, the tendency of the maximum amplitude of solute concentration fluctuation with the crucible rotation is similar to that of the crystal rotation. Specifically, the maximum amplitude of solute concentration fluctuation first decreases, then increases with the growth of the crucible rotation Reynolds number. The crucible rotation effectively suppresses the solute fluctuation of traveling waves at first. After that, the flow is gradually affected by the circular shear instability. Thus, the maximum amplitude of solute concentration fluctuation rises with the crucible rotation.

4. Conclusions

3D numerical simulations were performed to investigate the influence of a rotating crystal and crucible on the thermal-solutal capillary-buoyancy flow of a Ge_1–xSi_x melt in a Cz configuration. The critical condition for the initiation of a 3D oscillation flow was obtained, and the mechanism of flow destabilization was analyzed. Meanwhile, the basic characteristics of steady and 3D oscillatory flows were discussed. This work contributes to the development of the Cz crystal growth method. Based on these above results, some conclusions are drawn as follows.

1. If the thermal capillary Reynolds number is relatively low, the flow appears as a steady axisymmetric flow, even though the crystal and crucible rotate at constant values. A steady axisymmetric flow is an ideal flow condition for the Cz growth process, due to its good thermal and crystal uniformity. This steady axisymmetric flow is suppressed by the crystal or crucible rotation. The strength of the steady axisymmetric flow exists as a decreasing trend as the crystal or crucible rotation rate increases.

2. The critical thermal capillary Reynolds numbers for the onset of a 3D oscillatory flow in a Cz configuration are much greater than those of pure fluids. Furthermore, the critical thermal capillary Reynolds number first rises, and then falls as the crystal or crucible rotation rate increases. Clearly, the quality of a pure component crystal is much more essential to the temperature difference than that of a compound. If the crystal or crucible rotation rate is smaller than a certain value, the rotation benefits the quality of a crystal. However, the crystal or crucible rotation can have negative effects on the quality of a crystal at a large rotation rate.

3. Due to the suppression effect of the crystal or crucible rotation at a low rotation rate, a reversed evolution from an oscillatory spoke pattern to traveling waves is found as the crystal or crucible rotation rate increases. Once the crystal or crucible rotates at a large rate, the traveling waves transition to rotating waves in a crystal rotation system and a spindle-like pattern in a crucible rotation system. In addition, the maximum amplitude of the solute concentration fluctuation of the 3D oscillation flow initially decreases, but finally increases alongside the growth of the crystal or crucible rotation rate.