1. Introduction

Chemical vapor deposition (CVD) is the critical process of making compound semiconductor optical devices, including blue LED, high-power lasers, solar cells and so on. This technology uses organic compounds of Ⅱ, Ⅲ group elements and hydride of Ⅳ, Ⅴ group elements as the growth source material and conducts the vapor phase epitaxy on the substrate, depositing Ⅲ-Ⅴ, Ⅱ-Ⅳ group compound semiconductors and their multiple solid thin layer of single crystal^[1]. According to the flow direction of the main gas flow relative to the substrate, the reactor can be divided into two categories^[2]: one is a horizontal reactor, the main gas flow direction parallel to the substrate; the other is a vertical reactor, the main gas flow perpendicular to the substrate. Because the structure of the horizontal reactor system is simple and easy to operate, so it has been favored by a large number of scientists conducting laboratory research. In the design and application of CVD, the key issue is to obtain uniform film, which is closely related to the temperature and concentration fields. In view of this fact, previous study^[3] showed that when the Rayleigh number Ra > 2000, it is easy to see longitudinal rolls in the horizontal reactor. The longitudinal rolls would destroy the uniformity of the temperature and concentration fields, resulting in transverse non-uniform deposition of the film. In order to reduce this non-uniformity, Nicolas et al. ^[4] proposed and analyzed a method which made the parallel longitudinal rolls unsteady, periodic and sinuous in order to get a uniform time-average deposition in a laminar reactor. However, these studies focused on the laminar reactor and deliberately avoided the development of turbulence. Actually, the turbulence can be easily obtained with the increase of Ra. In 1985, an experimental result through interference holography from Giling^[5] pointed out the coexistence of a stable, laminar boundary layer near the susceptor, and an unstable turbulent layer above it was discerned. According to the fact that turbulence can enhance the heat transfer, Santen et al. ^[6] proposed to use turbulent vertical CVD reactors as a deposition model, showing that the deposition rate and the uniformity are far better in a large part of the susceptor than in a laminar state. In 2011, Tian et al.^[7] obtained a similar conclusion that turbulent CVD can get a better film growth. Then Liu et al. ^[8] studied the effects of varying the Reynolds number and Rayleigh number on the heat transfer in horizontal reactors, finding that when the flow pattern is within the soft turbulent stage, the time-averaged heat transfer is more uniform in the entire substrate. Depending on the limited reports of turbulence reactors, it was found that previous studies have mainly focused on the effects of turbulence on flow and temperature fields, not on the quality of the film deposition. In this study, numerical simulation applying the Reynolds Averaged Navier-Stokes (RANS) equations was undertaken to investigate the cases in the horizontal reactor, with longitudinal and transversal aspect ratios of 10 and 4, respectively. The simulation is based on simple reaction models that get silicon from silane (SiH$_{4})$, and directly investigates the effects of different Ra numbers and Re numbers on the Si deposition rate and uniformity, further exploring the feasibility of the turbulent CVD reactor.

2. Numerical model and boundary conditions

2.1 Physical model

A horizontal rectangular cold wall reactor is our numerical simulation and analysis case, and its schematic diagram is presented in Fig. 1. The longitudinal and transversal aspect ratios of this reactor are $A=L/H=$ 10 and $B=l/H=$ 4, and $l_{\rm r}/H=l_{\rm e}/H=$ 1, where $l$, $L$ and $H$ are the width, length and height of the reactor; $l_{\rm r}$ and $l_{\rm e}$ are the distances of the reserved developing and exit extension regions. Note that the capital $X, Y, Z$ are $x$/$H$, $y$/$H$ and $z$/$H$, respectively.

2.2 Computing method

Because the CVD process usually occurs at a temperature between 700-1300 K and the flow rate is low, so it can be considered as incompressible gas with properties depending on the local temperature. The governing equations of this model involve a continuity equation, momentum equation, energy equation and component equation, listed below:

Continuity equation:

$\begin{equation} \frac{\partial u_i }{\partial x_i }=0. \end{equation}$

(1)

Momentum equation:

$\begin{equation} \begin{split} \frac{\partial }{\partial t}(\rho u_i ) {}& +\frac{\partial }{\partial x_j }(\rho u_i u_j) \\& =-\frac{\partial p}{\partial x_i }+\frac{\partial }{\partial x_j}\left(\mu \frac{\partial u_i }{\partial x_j }-\rho \overline {u_i ^\prime u_j ^\prime } \right)+\rho g. \end{split} \end{equation}$

(2)

Energy equation:

$\begin{equation} \frac{\partial }{\partial t}(\rho {\rm T})+\frac{\partial }{\partial x_j }(\rho u_j T)=\frac{\lambda }{C_{\rm P} }\frac{\partial ^2T}{\partial x_j \partial x_j}. \end{equation}$

(3)

Component equation:

$\begin{equation} \frac{\partial }{\partial t}(\rho Y_i )+\frac{\partial }{\partial x_j }(\rho u_j Y_i )=\rho D_i \frac{\partial ^2Y_i }{\partial x_j \partial x_j }+\omega _i. \end{equation}$

(4)

Here $u$ is the gas velocity, $\mu$, $C_{\rm P}$ and $\lambda$ are the viscosity, specific heat and thermal conductivity, respectively, all of which are determined through kinetic theory. $T$ and $p$ denote the temperature and the pressure. $Y_i$, $D_i$, $\omega _i$ represent the mass fraction, diffusion coefficient and reaction rate of component $i$, respectively. The temperature dependence of $C_{\rm P}$, $\lambda$ and $\mu$ is determined through kinetic theory. The main work of the numerical simulation is to solve the governing equations. Because the RANS equations are not closed, so a realizable $k$-$\varepsilon$ turbulence model is introduced to ensure the governing equations can be calculated. The $k$, $\varepsilon$ transport equations are as follows^{[9, 10]}:

$\begin{equation} \begin{split} \frac{\partial (\rho k)}{\partial t}+\frac{\partial (\rho ku_i )}{\partial x_i } {}& =\frac{\partial }{\partial x_j }\left[(\mu +\frac{\mu _t }{\sigma _k })\frac{\partial k}{\partial x_j }\right] \\& +\mu _t \left(\frac{\partial u_i }{\partial x_j }+\frac{\partial u_j }{\partial x_j }\right)\frac{\partial u_i }{\partial x_j }-\rho \varepsilon, \end{split} \end{equation}$

(5)

$\begin{equation} \begin{split} \frac{\partial (\rho \varepsilon )}{\partial t}+\frac{\partial (\rho \varepsilon u_i )}{\partial x_i } {}& =\frac{\partial }{\partial x_j }\left[\left(\mu +\frac{\mu _t }{\sigma _\varepsilon }\right)\frac{\partial \varepsilon }{\partial x_j }\right] \\& +\rho C_1 E\varepsilon -\rho C_2 \frac{\varepsilon ^2}{k+\sqrt {\nu \varepsilon } }, \end{split} \end{equation}$

(6)

where $\sigma _k$ $=$ 1.0, $\sigma _\varepsilon$ $=$ 1.2, $C_2$ $=$ 1.9, $C_1 =\max (0.43, \frac{\eta }{\eta +5})$, $\eta =(2E_{ij} E_{ij} )^{1/2}\frac{k}{\varepsilon }$, $E_{ij} =\frac{1}{2}(\frac{\partial u_i }{\partial x_j }+\frac{\partial u_j }{\partial x_i })$, $\mu _t=\rho C_\mu \frac{k^2}{\varepsilon }$. In the realizable $k$-$\varepsilon$ turbulence model, $C_\mu$ is not a constant, and the $C_\mu$ expression is much more complicated, which can be found in Refs. [9, 10].

In the CVD reactor, what one is most concerned about is the thin film deposition rate and the uniformity. This paper focuses on the impact of the turbulent complex flow field on the Si deposition rate and the uniformity compared with a laminar reactor. The deposition model of silicon from silane includes one gaseous reaction (reaction 1) and two surface reactions (reactions 2 and 3), presented as follows:

Reaction 1: ${\rm SiH}_{4} ({\rm g})\to {\rm SiH}_{2} {\rm (g)+H}_{2} {\rm (g)}$,

Reaction 2: ${\rm SiH}_{4} {\rm (g)}\to {\rm Si(s)+2H}_2 {\rm (g)}$,

Reaction 3: ${\rm SiH}_2 {\rm (g)}\to {\rm Si(s)+H}_2 {\rm (g)}$.

Reaction rate meets the Arrhenius equation:

$\begin{equation} m_{\rm f, r} =A_{\rm r} T^{\beta _{\rm r}}{\rm e}^{-E_{\rm r}/RT}, \end{equation}$

(7)

where $A_{\rm r}$ is the pre-exponential factor, $\beta _{\rm r}$ is the temperature exponent, $E_{\rm r}$ is the reaction activation energy, and $R$ is the gas constant.

2.3 Boundary condition

The boundary conditions applied to the computational domain are as follows:

(1) At the inlet, a mixture gas of SiH$_{4}$ and N$_{2}$ (the mass fraction of SiH$_{4}$ is 2%, nitrogen is the carrier gas) is injected from the inlet with 300 K and flows out from the right (Fig. 1). The spatial distribution of the velocity is assumed to be uniform for this analysis.

(2) At the side walls, a no-slip and an adiabatic condition is imposed on the velocity and temperature, respectively.

(3) At the upper wall, a no-slip condition is imposed on the velocity for the entire wall. An adiabatic condition is imposed on the temperature except for the reaction section where the temperature is 300 K.

(5) At the lower wall, a no-slip condition is imposed on the velocity for the entire wall. An adiabatic condition is imposed on the temperature except for the substrate where the temperature is 1300 K.

(4) At the exit, the free outflow is adopted.

2.4 Solution method and model validation

A 3-D structured grid is used with a control volume formulation in accordance with the SIMPLE algorithm. Pressure-correction and velocity-correction schemes are implemented in this algorithm to arrive at a converged solution. The 2nd-order upwind scheme is employed for the discretization of the momentum equation. Different under-relaxation factors are used for the pressure, velocity and temperature solutions. In order to verify the grid-independence, the cases of Ra $=$ 1.60 $\times$ 10$^{4}$ are taken as examples. Different meshing and relative average deposition rates of the substrate $v_{\rm m}$ are listed in Table 1. The result shows that the meshing of $X \times Y \times Z=$ 250 $\times$ 100 $\times$10 can meet the grid-independence, and the meshing of subsequent computations are based on corresponding proportions.

In an attempt to verify the accuracy of the simulation model, a case with longitudinal and transversal aspect ratios $A$ $=$ 10, $B$ $=$ 200 was simulated by RANS model and the result was compared with that from the direct numerical simulation (DNS) extracted from Ref. [4] and presented as follows. A close examination of the temperature field shown in Fig. 2 indicates that the number, the spanwise gap and the onset points of longitudinal rolls modeled by RANS basically agree with the result given by DNS. Accordingly, the RANS numerical method can be used to simulate the cases involved in this study.

3. Results and discussion

In the process of CVD, the deposition rate and uniformity of thin film are the key issues. The turbulent flow can cause the instantaneous deposition to be highly non-uniform, however, the time-average of that may be uniform in some certain part. It is not a necessary problem in CVD for an instantaneously non-uniform deposition, and obtaining a uniform averaged deposition over a required growth time is essential. It is well known that there are two characteristic parameters affecting the turbulence. The Reynolds number, Re $=$ ($u_0 H)/v$, represents the strength of the forced flow, thus, increasing the Re number can be expected to improve the effect of the forced flow. The Rayleigh number, Ra $=$ ($g\beta \Delta TH^3)/(vk)$, can be used to depict the complex buoyancy induced vortex flows. Here, $u_0$ is the inlet velocity, $\beta$, $v$ and $k$ are thermal expansion coefficient, kinematic viscosity and thermal diffusivity, respectively. $H$ is the height of the reactors and $\Delta T$ denotes the temperature difference between the top and bottom wall. In order to study the deposition rate and uniformity of silicon growth from silane in turbulent reactors, varying Ra and Re numbers are explored below.

3.1 Effect of the Rayleigh number

It is well known that the buoyancy induced convection can be strengthened by increasing Ra numbers. Specifically, the convection can be divided into several regions^[11], such as the onset region, transition region (Ra $=$ 2.5 $\times$ 10$^{5}$-5 $\times$ 10$^{5})$ and soft turbulence region (Ra $=$ 5 $\times$ 10$^{5}$-4 $\times$ 10$^{7})$, in terms of the Ra number. In the onset region, the mixed parameter Ra/(PrRe$^{2})$ was not high enough, therefore, stable longitudinal rolls may emerge in horizontal reactors. With Ra/(PrRe$^{2})$ increasing, the buoyancy is improved by the ascending Ra number. As Ra increased to the soft turbulence stage, the regular rolls become unstable and break up eventually. The shape and location of the vortex cell also affect the internal mass transfer, thereby affecting the deposition rate and uniformity.

From Table 2, it has been found that the deposition rate of the turbulent reactor (Ra $=$ 2.24 $\times$ 10$^{6})$ is faster than the laminar reactor (Ra $=$1.6 $\times$ 10$^{4})$. In the view of industrial production, the increase of the deposition rate can improve the production efficiency and bring better economic benefits. It is well known that the chemical reaction rate depends on the Arrhenius equation that are influenced by the pre-exponential factor, temperature index and activation energy. By examining the given values of these parameter, it can be easily found that reaction 3 is much faster than reaction 2, so the SiH$_{2}$ concentration distribution is the main factor affecting the deposition rate. The instantaneous concentration distribution of SiH$_{2}$ above the substrate 5 mm at Re $=$ 50 is shown in Fig. 3. Clearly, it can be seen that the SiH$_{2}$ concentration distribution at Ra $=$ 2.24 $\times$ 10$^{6}$ is higher and more uniform than that at Ra $=$ 1.6 $\times$ 10$^{4}$, which can be explained by the turbulence strengthening the mixture of reactant gas, making more SiH$_{2}$ molecules reach the substrate.

In this study, different Ra numbers were obtained by changing the height of the reactors. In order to compare conveniently with various cases, a modified deposition rate $D$ was given and a dimensionless symbol, $D=(v_{\rm i} -v_{\rm m})/v_{\rm m}$, was proposed to denote the deposition rate. It can be discerned from Figs. 4 and 5 that the distribution of $D$, i.e. $v$, averaged over a long time along the Y direction varied with the Ra numbers at different cross sections $X=$ 4 and $X=$ 6. Firstly, we noticed from Fig. 4 that the deposition of Si fluctuates strongly along the $Y$ direction at Ra $=$ 1.6 $\times$ 10$^{4}$. With the increase of Ra, the natural convection gradually enhanced, and the internal plumes move more fiercely and randomly, so the fluctuation of deposition becomes more and more mild and uniform. When the Ra increases to the soft turbulent stage, at Ra $=$ 2.24 $\times$ 10$^{6}$, the $D$ in this region fluctuated more weakly than in any of the other cases studied in this paper. The same conclusion can also be obtained by observing Fig. 5. In short, when the flow is located in the soft turbulent stage, the internal heat and mass transfer have a better distribution, which is conducive to deposition.

3.2 Effect of the Reynolds number

In horizontal reactors, it is common to see longitudinal and transversal rolls whose rotation axes are parallel and perpendicular to the main flow direction, respectively, but the latter only emerges at low Re numbers^[12]. The changes of Re number have a significant effect on the formation of the longitudinal rolls in the reactor, which affect the deposition rate and uniformity of Si. The distribution of the deposition rate for varying Re are shown in Figs. 6 and 7 at Ra $=$ 2.24 $\times$ 10$^{6}$. It can be seen from Fig. 6 that when Re $=$ 50, the fluctuation of the deposition rate is slightly milder than in the other two situations. It can be explained that the higher Re number can be expected to boost the influence of the forced flow, and then prevent the hot plumes, ascending from the hot walls, from distributing randomly and cause non-uniformity in a range of times. Moreover, it can also be seen in Fig. 7 that the deposition rate increases significantly at Re $=$ 100 compared to that at Re $=$ 50. It indicates that when Re is in the small range, the main reason affecting the deposition rate is the reactants concentration. Increasing the inlet velocity will compensate the reactants loss along the duct, which greatly improves the downstream deposition rate. But increasing the Re number will damage the transversal distribution of the deposition rate, resulting in a pronounced fluctuation. When Re $=$ 150, it has been found that the deposition rate of Si does not increase, which indicates that the loss of reactants along the way has been fully compensated.

4. Conclusion

Generally, it is too difficult to attain uniform films in CVD processes due to instantaneous non-uniform heat transfer, resulted by turbulence. So, in industry, people do not dare to exceed the critical Ra number in case of turbulence, which occurs easily in many applications. Practically, turbulent CVD can also achieve uniform time-averaged heat transfer in a large part of the susceptor^[8]. In this study, the effect of turbulence on the deposition of silicon in a horizontal reactor has been numerically studied using RANS.

In the present simulation, the effect of different Ra numbers and Re numbers on silicon deposition were obtained. Because the turbulence induced by free-convection can enhance the mixture of reactant gases, the deposition rate of the turbulent reactor is faster than the laminar reactor. When the Ra number increases to the soft turbulence stage, a fairly uniform deposition also can be obtained. Increasing the Re number within a certain range will compensate the reactants along the path losses, which can significantly increase the downstream deposition rate. But the enhancement of forced convection can prevent the plume from distributing randomly, resulting in transverse non-uniformity of the deposition rate and \linebreak deteriorate the deposition quality. Therefore, in practical applications of turbulent CVD, it is very important to control the Re number and Ra number in order to get a relatively uniform film deposition.