Numerical Simulation of Fluidization Behavior and Chemical Performance for Hydrochlorination of Silicon Tetrachloride in a Fluidized Bed Reactor

Abstract

Exploring the fluidization behaviors and chemical performance in silicon tetrachloride (SiCl₄) hydrochlorination processes within a fluidized bed reactor (FBR) poses significant challenges. In this study, we developed an Eulerian-granular model (EGM) by integrating the Eulerian–Eulerian two-fluid model with the kinetic theory of granular flow (KTGF). The effect of fluidization velocities on the flow regime, heat transfer, and chemical reaction performance were investigated. The applicability of the simulation method and the validity of the model were confirmed through comprehensive comparisons, including the simulated values of the maximum bed expansion height (H_max) with theoretical values derived from empirical formulas and the simulated gas temperature profile with Hsu’s experimental data. The results indicate that the present EGM can be feasible to describe the variation of the flow regime within the FBR. An increase in bed voidage over time, coinciding with transitions in the flow regime, can be observed. Particularly noteworthy was the attainment of a more uniform distribution of SiCl₄ under the bubbling fluidization state. Furthermore, the FBR possess high heat transfer characteristics, and the reaction gas can reach the set temperature of the bed after entering a small distance (about 10 mm). The presence of circulating bubbles within the FBR enhances the mixing uniformity of the SiCl₄ reaction gas and silicon particles, particularly in the central and upper regions of the bed under the bubbling fluidization state. As a result, the predicted highest concentration of SiHCl₃ was 13.08% and the conversion rate of SiCl₄ was 28.97% under the bubbling fluidization state. Our results can provide a theoretical basis for further understanding of the hydrochlorination process of SiCl₄ within the FBR.

1. Introduction

Polycrystalline silicon plays a crucial role as the primary raw material in the production of monocrystalline silicon substrates for semiconductor devices and as a substrate for solar cell manufacturing [1,2]. At present, the popular technique for producing polysilicon is the Siemens method [3], by which the high-purity polysilicon rods are prepared by chemical vapor deposition (CVD) of purified SiHCl₃ and H₂ in a Siemens reactor. The production process is well-established, stable, and consistently yields high-quality products [4]. However, it is worth noting that for every ton of polysilicon produced, approximately 15 to 20 tons of by-products, specifically silicon tetrachloride (SiCl₄), are generated [5]. If SiC_l4 is not effectively utilized, it not only represents a waste of silicon-rich chemical resources but also poses potential risks to human health and environmental quality [6,7]. Therefore, finding effective solutions for handling SiCl₄ has emerged as a significant challenge in this field.

Nowadays, the most effective approach for managing SiCl₄ involves the utilization of hydrogenation technology to transform SiCl₄ into SiHCl₃. This not only yields the desired SiHCl₃ feedstock but also mitigates the risk of secondary pollution stemming from SiCl₄. The central apparatus for SiCl₄ hydrogenation [8] is the fluidized bed reactor (FBR), and the primary chemical reaction is depicted in Equation (1). Furthermore, the hydrochlorination process [9], which represents an optimized hydrogenation technology, introduces HCl as an additional reactant. This process offers notable advantages, including high conversion rates and low energy consumption [10]. The principal chemical reaction in hydrochlorination is illustrated in Equation (2). As such, it holds immense significance to investigate the structural stability of the FBR, ensure a uniform flow distribution, assess the heat transfer efficiency, and evaluate the chemical performance. These endeavors are essential for enhancing the hydrochlorination process and maximizing the conversion rate of SiCl₄.

$$3{\mathrm{SiCl}}_4+{2\mathrm{H}}_2+\mathrm{Si}=4{\mathrm{SiHCl}}_3$$

(1)

$$2{\mathrm{SiCl}}_4+{\mathrm{H}}_2+\mathrm{Si}+\mathrm{HCl}=3{\mathrm{SiHCl}}_3$$

(2)

To date, the research on the chemical reaction process of converting SiCl₄ into SiHCl₃ has mainly focused on the thermodynamics and kinetics of the FBR at a laboratory scale. Lee et al. [11] investigated the hydrogenation of SiCl₄ to SiHCl₃ using a carbon-based catalyst with the experimental method, and the results showed that the yield of SiHCl₃ can be improved by retarding the reverse reaction rate of SiHCl₃ and HCl. Jain et al. [12] studied the hydrochlorination of Si to produce SiHCl₃ under experimental conditions of the FBR, and the reaction rate constants of the HCl reaction on the Si substrate were obtained. Becker et al. [13] employed a kinetic model to describe the hydrochlorination behavior of SiCl₄, and the effects of the operation mode on the SiHCl₃ yield were systemically analyzed. It is well known that the chemical reaction is also affected by transfer phenomena; thus, only focusing on the chemical thermodynamics and kinetics is insufficient.

Over the past few decades, many experimental and theoretical investigations of the gas–solid fluidization behavior in the FBR have been performed by researchers [14,15]. Shao et al. [16] systematically investigated the influences of temperature and pressure on the minimum fluidization velocity (U_mf) under different particle size distributions. Their results indicated that the U_mf decreased with increasing temperature and pressure for both the narrow and wide particle size distributions. Chein et al. [17] employed numerical simulations to examine the chemical performance of a combustion system in the traditional bubbling fluidized bed (BFB) and the BFB combined with internal particle circulation (ICBFB). The results showed that the ICBFB can increase the mixing time of reactants and enhance the chemical reaction performance. Pottimurthy et al. [18] constructed and analyzed 3D images of the mean solid volume fraction distribution in a slugging fluidized bed. The results demonstrated that three zones of bottom fluidization, gas slug, and solid slug co-existed in the bed, and the slugging fluidization velocity had a significant effect on the maximum length of the gas slug zone. Liu et al. [19,20] employed an Eulerian–Eulerian two-fluid model, coupled with population balance equations, to simulate the progression of silicon particle growth through the chemical vapor deposition of silane pyrolysis in a slugging fluidized bed reactor using FLUENT. The simulation of particle growth, which takes into account surface deposition, cluster scavenging, aggregation, and wall deposition, was conducted following the verification of the flow and heat transfer characteristics using widely accepted correlations. There have also been many other studies on fluidization characteristic parameters, including fluidization velocity [21,22], bed expansion height [23,24], solid volume fraction, and voidage [25,26], etc. Nevertheless, there have been limited numerical investigations into the fluidization behavior associated with SiCl₄ hydrochlorination for the production of SiHCl₃. Colomb et al. [27,28] developed a fluidized bed model based on the Kunii–Levenspiel fluidization framework, which can accurately predict the reaction performance of the hydrochlorination process. However, the effects of the fluidization behavior on the chemical reaction performance were ignored in the above-proposed model. Liu et al. [29] applied the computational fluid dynamics (CFD) method to investigate the influence of baffles on the hydrogenation process of SiCl₄ in an FBR. Their results showed that the louver baffles can improve the reaction performance of SiCl₄, and the predicted value of the conversion rate in the FBR with louver baffles was higher than that with channel baffles.

The aforementioned reports have made significant contributions to understanding the flow behaviors and heat transfer characteristics within the FBR. However, the previous studies have inadequately explored the fluidization behaviors specifically related to SiCl₄ hydrochlorination processes within the FBR. Additionally, there has been limited investigation into the impact of fluidization behaviors on chemical reaction performance. Hence, in this study, we established an EGM using the CFD approach. We systematically analyzed the effects of fluidization velocities on the flow regimes, heat transfer, and chemical reaction performance. Our primary objective was to offer theoretical insights that can support the enhancement of FBR fluidization stability and the improvement of SiCl₄ conversion rates.

2. Model Description

2.1. CFD Model

2.1.1. Eulerian-Granular Model

The Eulerian-granular model (EGM) was established by the Eulerian–Eulerian two-fluid model with the kinetic theory of granular flow (KTGF), which was employed to simulate the flow behaviors between the gas–solids in the FBR [30,31,32]. Considering the granular motion is primarily influenced by particle kinetic theory and direct collisions within the solid fluid, the detailed equations [33,34] of the EGM and KTGF are shown as follows.

The continuity equations are depicted in Equations (3) and (4), where the subscripts g and s represent the gas phase and solid phase, respectively; α is the volume fraction; the arrow represents the vector; ρ represents the density; and m_gs and m_sg are the mass transfer of the interphase mass transfer from gas to solid and solid to gas, respectively.

$$\frac{𝜕}{𝜕\mathrm{t}}({\alpha }_{\mathrm{g}}{\rho }_{\mathrm{g}})+\nabla \cdot ({\alpha }_{\mathrm{g}}{\rho }_{\mathrm{g}}{\overrightarrow{v}}_{\mathrm{g}})=m_{\mathrm{gs}}$$

(3)

$$\frac{𝜕}{𝜕\mathrm{t}}({\alpha }_{\mathrm{s}}{\rho }_{\mathrm{s}})+\nabla \cdot ({\alpha }_{\mathrm{s}}{\rho }_{\mathrm{s}}{\overrightarrow{v}}_{\mathrm{s}})=m_{\mathrm{sg}}$$

(4)

The momentum equations are depicted in Equations (5)–(8), where the τ is shear stress; p represents the pressure drop; μ represents the viscosity; K_gs and K_sg are the interphase exchange coefficients; λ_s is solid bulk viscosity and I represents identity matrix.

$$\frac{𝜕}{𝜕\mathrm{t}}({\alpha }_{\mathrm{g}}{\rho }_{\mathrm{g}}{\overrightarrow{v}}_{\mathrm{g}})+\nabla \cdot ({\alpha }_{\mathrm{g}}{\rho }_{\mathrm{g}}{\overrightarrow{v}}_{\mathrm{g}}{\overrightarrow{v}}_{\mathrm{g}})=-{\alpha }_{\mathrm{g}}{\rho }_{\mathrm{g}}g+\nabla {\tau }_{\mathrm{g}}-{\alpha }_{\mathrm{g}}\nabla p+K_{\mathrm{gs}}({\overrightarrow{v}}_{\mathrm{s}}-{\overrightarrow{v}}_{\mathrm{g}})-m_{\mathrm{gs}}\overrightarrow{v}{}_{\mathrm{g}}$$

(5)

$${\tau }_{\mathrm{g}}={\alpha }_{\mathrm{g}}{\mu }_{\mathrm{g}}(\nabla {\overrightarrow{v}}_{\mathrm{g}}+\nabla {\overrightarrow{v}}_g^T)$$

(6)

$$\frac{𝜕}{𝜕\mathrm{t}}({\alpha }_{\mathrm{s}}{\rho }_{\mathrm{s}}{\overrightarrow{v}}_{\mathrm{s}})+\nabla \cdot ({\alpha }_{\mathrm{s}}{\rho }_{\mathrm{s}}{\overrightarrow{v}}_{\mathrm{s}}{\overrightarrow{v}}_{\mathrm{s}})=-{\alpha }_{\mathrm{s}}{\rho }_{\mathrm{s}}g+\nabla {\tau }_{\mathrm{s}}-{\alpha }_{\mathrm{s}}\nabla p+K_{\mathrm{sg}}({\overrightarrow{v}}_{\mathrm{g}}-{\overrightarrow{v}}_{\mathrm{s}})-m_{\mathrm{sg}}\overrightarrow{v}{}_{\mathrm{s}}$$

(7)

$${\tau }_{\mathrm{s}}={\alpha }_{\mathrm{s}}{\mu }_{\mathrm{s}}(\nabla {\overrightarrow{v}}_{\mathrm{s}}+\nabla {\overrightarrow{v}}_{\mathrm{s}}^T)+{\alpha }_{\mathrm{s}}({\lambda }_{\mathrm{s}}-\frac{2}{3}{\mu }_{\mathrm{s}})\nabla {\overrightarrow{v}}_{\mathrm{s}}I$$

(8)

The energy equations are depicted in Equations (9)–(11), where the subscript m represents the gas phase and solid phase, respectively; h is the specific enthalpy; C_p represents the heat capacity; k is the turbulence kinetic energy tensor; and T is the temperature.

$$\frac{𝜕}{𝜕t}({\alpha }_{\mathrm{m}}{\rho }_{\mathrm{m}}{h}_{\mathrm{m}})+\nabla \cdot ({\alpha }_{\mathrm{m}}{\rho }_{\mathrm{m}}{\overrightarrow{v}}_{\mathrm{m}}{h}_{\mathrm{m}})={\alpha }_{\mathrm{m}}\frac{𝜕p_{\mathrm{m}}}{𝜕t}+{\tau }_{\mathrm{m}}\nabla \cdot {\overrightarrow{v}}_{\mathrm{m}}-\nabla \cdot q_{\mathrm{m}}-(Q_{\mathrm{gs}}+m_{\mathrm{gs}}{h}_{\mathrm{gs}}-m_{\mathrm{sg}}{h}_{\mathrm{sg}})$$

(9)

$$h_{\mathrm{m}}={\displaystyle {\int }_{T_{\mathrm{ref}}}^T{C}_{\mathrm{p},\mathrm{m}}}dT$$

(10)

$$q_{\mathrm{m}}=-{\alpha }_{\mathrm{m}}{k}_{\mathrm{m}}\nabla T_{\mathrm{m}}$$

(11)

The solid pressure and radial distribution function are described in Equations (12) and (13), where Θ_s is the granular temperature; e_ss is the particle–particle restitution coefficient; and g₀ is the radial distribution function.

$$P_{\mathrm{s}}={\alpha }_{\mathrm{s}}{\rho }_{\mathrm{s}}{\Theta }_{\mathrm{s}}+2{\rho }_{\mathrm{s}}(1+e_{\mathrm{ss}}){\alpha }_{\mathrm{s}}^2{g}_{0,\mathrm{ss}}{\Theta }_{\mathrm{s}}$$

(12)

$$g_0={\left[1-{(\frac{{\alpha }_{\mathrm{s}}}{{\alpha }_{\mathrm{s},\max }})}^{1/3}\right]}^{-1}$$

(13)

The solid shear stress, which includes the collisional viscosity, kinetic viscosity and frictional viscosity is shown in Equations (14)–(17), where μ_s, μ_s,col, μ_s,kin, and μ_s,fr represent the solid shear stresses, collisional viscosity, kinetic viscosity, and frictional viscosity, respectively; d_s represents the diameter of the solid particles.

$${\mu }_{\mathrm{s}}={\mu }_{\mathrm{s},\mathrm{col}}+{\mu }_{\mathrm{s},\mathrm{kin}}+{\mu }_{\mathrm{s},\mathrm{fr}}$$

(14)

$${\mu }_{\mathrm{s},\mathrm{col}}=\frac{4}{5}{\alpha }_{\mathrm{s}}{\rho }_{\mathrm{s}}{d}_{\mathrm{s}}{g}_{0,\mathrm{ss}}(1+e_{\mathrm{ss}})\sqrt{\frac{{\Theta }_{\mathrm{s}}}{\pi }}$$

(15)

$${\mu }_{\mathrm{s},\mathrm{kin}}=\frac{10{\rho }_{\mathrm{s}}{d}_{\mathrm{s}}\sqrt{{\Theta }_{\mathrm{s}}\pi }}{96{\alpha }_{\mathrm{s}}(1+e_{\mathrm{ss}})g_{0,\mathrm{ss}}}{\left[1+\frac{4}{5}{g}_{0,\mathrm{ss}}{\alpha }_{\mathrm{s}}(1+e_{\mathrm{ss}})\right]}^{{}^2}$$

(16)

$${\mu }_{\mathrm{s},\mathrm{fr}}=\frac{p_{\mathrm{s}}\sin \varphi }{2\sqrt{I_{2D}}}$$

(17)

The granular temperature is shown in Equation (18), the collision dissipation of energy is shown in Equation (19), and the transfer of the kinetic energy of random fluctuations is shown in Equation (20), where ${\gamma }_{{\mathsf{\Theta }}_{\mathrm{s}}}$ is the collisional dissipation of energy fluctuation; ${\varphi }_{\mathrm{gs}}$ is the exchange energy between gas and solid.

$$0=(-p_{\mathrm{s}}I+{\tau }_{\mathrm{s}})\nabla \cdot v_{\mathrm{s}}-{\gamma }_{{\mathsf{\Theta }}_{\mathrm{s}}}+{\varphi }_{\mathrm{gs}}$$

(18)

$${\gamma }_{{\mathsf{\Theta }}_{\mathrm{m}}}^{}=\frac{12(1-e^2{}_{\mathrm{ss}}{g}_{0,\mathrm{ss}})}{d_{\mathrm{s}}\pi }{\alpha }^2{}_{\mathrm{s}}{\rho }_{\mathrm{s}}{\Theta }_{\mathrm{s}}^{3/2}$$

(19)

$${\varphi }_{\mathrm{gs}}=-3K_{\mathrm{gs}}{\Theta }_{\mathrm{s}}$$

(20)

2.1.2. Turbulent Model

In this work, the RNG k-ε turbulent model and standard wall function were employed to solve the turbulent viscosity. The transient-state transport equations are expressed as follows: Equations (21)–(25), where k is the turbulent kinetic energy tensor, ε represents the turbulent dissipation rate, and G_k represents the generation of turbulence kinetic energy. C_1ε, C_2ε, C_ν and C_μ are constants, and C_1ε = 1.44, C_2ε = 1.92, C_ν = 100, and C_μ = 0.0845; μ represents the gas viscosity and μ_eff is the effective solid viscosity.

$$\frac{𝜕}{𝜕t}\left({\rho }_{\mathrm{m}}k\right)+\frac{𝜕}{𝜕x_i}({\rho }_{\mathrm{m}}k{\overrightarrow{v}}_{{}_{\mathrm{m}}})=\frac{𝜕}{𝜕x_j}\left[{\alpha }_k{\mu }_{\mathrm{eff}}\frac{𝜕k}{𝜕x_j}\right]+G_{k,\mathrm{m}}-{\rho }_{\mathrm{m}}\epsilon$$

(21)

$$\frac{𝜕}{𝜕t}\left({\rho }_{\mathrm{m}}\epsilon \right)+\frac{𝜕}{𝜕x_i}({\rho }_{\mathrm{m}}\epsilon {\overrightarrow{v}}_{{}_{\mathrm{m}}})=\frac{𝜕}{𝜕x_j}\left[{\alpha }_{\mathsf{\epsilon }}{\mu }_{\mathrm{eff}}\frac{𝜕\epsilon }{𝜕x_j}\right]+C_{1\mathsf{\epsilon }}\frac{\epsilon }{k}{G}_{k,\mathrm{m}}-C_{2\epsilon }{\rho }_{\mathrm{m}}\frac{{\epsilon }^2}{k}$$

(22)

$$d\left(\frac{{\rho }^{2}k}{\sqrt{\epsilon \mu }}\right)=1.72\frac{v}{\sqrt{v^3-1+C_v}}dv$$

(23)

$$v=\frac{{\mu }_{eff}}{\mu },C_v\approx 100$$

(24)

$${\mu }_t=\rho C_{\mu }\frac{k^2}{\epsilon },C_{\mu }=0.0845$$

(25)

2.1.3. Drag Model

For this case, the interphase exchange coefficient (K_gs) was described by the empirical drag model proposed by Gidaspow et al. [35], as defined by Equations (26)–(28). Here, C_D signifies the drag coefficient of individual spherical particles; Re_s denotes the Reynolds number of the solid particles.

$$K_{\mathrm{gs}}=\frac{3}{4}{C}_{\mathrm{D}}\frac{{\alpha }_{\mathrm{s}}{\alpha }_{\mathrm{g}}{\rho }_{\mathrm{g}}\left|{\overrightarrow{v}}_{\mathrm{g}}-{\overrightarrow{v}}_{\mathrm{s}}\right|}{d_{\mathrm{s}}}{\alpha }_{\mathrm{g}}^{-2.65}, {\alpha }_{\mathrm{g}}>0.8$$

(26)

$$C_{\mathrm{D}}=\frac{24}{{\alpha }_{\mathrm{g}}R{e}_{\mathrm{s}}}\left[1+0.15{({\alpha }_{\mathrm{g}}R{e}_{\mathrm{s}})}^{0.687}\right]$$

(27)

$$K_{\mathrm{gs}}=150\frac{{\alpha }_{\mathrm{s}}(1-{\alpha }_{\mathrm{g}}){\mu }_{\mathrm{g}}}{{\alpha }_{\mathrm{g}}{d}^2{}_{\mathrm{s}}}+1.75\frac{{\alpha }_{\mathrm{s}}{\rho }_{\mathrm{g}}\left|{\overrightarrow{v}}_{\mathrm{g}}-{\overrightarrow{v}}_{\mathrm{s}}\right|}{d_{\mathrm{s}}}, {\alpha }_{\mathrm{g}}\le 0.8$$

(28)

2.1.4. Heat Transfer Model

Assuming that the wall radiation was shielded by the erratic movement of particles and that the heating wall remained adiabatic, the heat transfer coefficient (h_w) for the wall-to-bed was determined using the following Equation (29), where the heat flux, q, was derived from a simulation, and T_w and T_f represent the temperature of the wall and fluid, respectively.

$$h_{\mathrm{w}}=\frac{q}{T_{\mathrm{w}}-T_{\mathrm{f}}}=\frac{q}{\frac{1}{H}{\displaystyle \underset{0}{\overset{\mathrm{H}}{\int }}\Delta Tdh}}$$

(29)

For the gas–solid heat transfer, the heat transfer coefficient (h_gs) was calculated by the Ranz–Marshall [36] and Guun [37] models. The Ranz–Marshall model equations were given as follows (Equations (30) and (31), respectively). The detailed Guun’s model is expressed as follows (Equations (32)–(35)). In those equations, Nu_s is the Nusselt number of the solid phase; Pr is the Prandtl number of the gas phase; k_g and k_s represent the gas and solid thermal conductivity, respectively, and C_p is the heat capacity.

$$h_{\mathrm{gs}}=\frac{k_{\mathrm{s}}N{u}_{\mathrm{s}}}{d_{\mathrm{s}}}$$

(30)

$$N{u}_{\mathrm{s}}=2.0+0.6R{e}_{\mathrm{s}}^{0.5}P{r}^{1/3}$$

(31)

$$h_{\mathrm{gs}}=\frac{6k_{\mathrm{g}}{\alpha }_{\mathrm{s}}{\alpha }_{\mathrm{g}}N{u}_{\mathrm{s}}}{d_{\mathrm{s}}^2}$$

(32)

$$N{u}_{\mathrm{s}}=(7-10{\alpha }_{\mathrm{g}}+5{\alpha }_{\mathrm{g}}^2)(1+0.7R{e}_{\mathrm{s}}^{0.2}P{r}^{1/3})+(1.33-2.4{\alpha }_{\mathrm{g}}+1.2{\alpha }_{\mathrm{g}}^2)R{e}_{\mathrm{s}}^{0.7}P{r}^{1/3}$$

(33)

$$Pr=\frac{C_{\mathrm{p}}^{\mathrm{g}}{\mu }_{\mathrm{g}}}{k_{\mathrm{g}}}$$

(34)

$$R{e}_{\mathrm{s}}=\frac{{\rho }_{\mathrm{g}}\left|{\overrightarrow{v}}_{\mathrm{g}}-{\overrightarrow{v}}_{\mathrm{s}}\right|d_{\mathrm{s}}}{\mu }$$

(35)

2.1.5. Species Transport Model

The mass transport behavior of the reactants within the FBR was solved using the species transport model, which is written as Equation (36), where the subscript n is the species; Y_n represents the fraction of each species; J_n represents the mass diffusion; R_n represents the net rate at which a chemical reaction produces a substance; and S_n is the source term.

$$\frac{𝜕}{𝜕t}(\rho Y_{\mathrm{n}})+\nabla \cdot (\rho \overrightarrow{v}{Y}_{\mathrm{n}})=-\nabla \cdot {\overrightarrow{J}}_{\mathrm{n}}+R_{\mathrm{n}}+S_{\mathrm{n}}$$

(36)

2.2. Chemical Reaction Kinetic Model

Finally, the kinetic model of Liu et al. [29] was coupled with the EGM model in our simulation. The chemical reaction rate of the hydrogenation (R_c1) and the hydrochlorination (R_c2) of SiCl₄ can be expressed by Equations (37)–(40), respectively, where k₀ is the apparent rate constant. The corresponding reaction kinetic parameters are shown in

Table 1. Reaction kinetics parameters [29].

Parameters	Values
Pre-exponential factor, A	0.22 mol/(g∙s∙atm2)
Activation energy, E	77.01 kJ/mol−1
Reference temperature, Tr	773 K
Equilibrium rate constant, Kp	0.2645 atm−1/3
Adsorption equilibrium constants, K1	0.2141 atm−1
Adsorption equilibrium constants, K2	8.92 atm−1
Enthalpy change, ΔH1	32.45 kJ/mol
Enthalpy change, ΔH2	−12.85 kJ/mol
Enthalpy change, ΔHP	−36.73 kJ/mol

.

$$\mathrm{Hydrogenation}:R_{\mathrm{c}1}=\frac{k_0{p}_{{\mathrm{SiCl}}_4}{p}_{{\mathrm{H}}_2}\left(1-\frac{p_{{\mathrm{SiHCl}}_3}^{4/3}}{K_{\mathrm{p}}{p}_{{\mathrm{SiCl}}_4}{p}_{{\mathrm{H}}_2}^{2/3}}\right)}{1+K_1{p}_{{\mathrm{SiCl}}_4}+\frac{P_{{\mathrm{SiHCl}}_3}^{1/3}{P}_{{\mathrm{H}}_2}^{1/3}}{K_2}}$$

(37)

$$\mathrm{Hydrochlorination}:R_{\mathrm{c}2}=\frac{k_0{p}_{{\mathrm{SiCl}}_4}{p}_{{\mathrm{H}}_2}{p}_{\mathrm{HCl}}\left(1-\frac{p_{{\mathrm{SiHCl}}_3}^{2/3}}{K_{\mathrm{p}}{p}_{{\mathrm{SiCl}}_4}{p}_{{\mathrm{H}}_2}^{1/2}{p}_{\mathrm{HCl}}^{1/2}}\right)}{1+K_1{p}_{{\mathrm{SiCl}}_4}+\frac{P_{{\mathrm{SiHCl}}_3}^{1/2}{P}_{{\mathrm{H}}_2}^{1/2}{P}_{\mathrm{HCl}}^{1/2}}{K_2}}$$

(38)

$$k_0=\mathrm{A}\exp [\frac{E}{R}(\frac{1}{T_r}-\frac{1}{T})]$$

(39)

$$K_n=K_{n,0}\exp [\frac{\Delta H_n}{R}(\frac{1}{T_{\mathrm{r}}}-\frac{1}{T})],n=p,1,2$$

(40)

2.3. Numerical Solution Algorithms

In this work, the simulation was carried out in commercial code ANSYS FLUENT 14.5. The computational domain was discretized using a structured mesh of 20,000 grids via the pre-processing module ANSYS Mesh 14.5. All the governing equations were solved using the SIMPLE algorithm, which was used for pressure–velocity coupling and correlation. The first-order upwind scheme was applied for the spatial discretization of the variables. Moreover, the inlet was defined as the velocity inlet and the outlet was the pressure outlet. On the wall, a non-slip boundary condition was set for the gaseous and solid phases. The basic parameter settings and boundary conditions in the present simulations are shown in

Table 2. Basic parameter settings [19,38].

Parameters	Values
Bed height	1.22 m
Bed width	0.1524 m
Solid density	2330 kg/m3
Particle diameter	235 μm
Solid–solid restitution coefficient	0.9
Gas density	1.225 kg/m3
Gas viscosity	1.72 μPa·s
Initial bed height	0.60 m
Initial solid volume fraction	0.60
Inlet gas velocity	0.15/0.30/0.65 m/s
Wall temperature	773 K
Inlet temperature	300 K
Operating pressure	2.0 MPa
SiCl4:H2:HCl	1:4:1

and

Table 3. Boundary conditions in present simulations.

Parameters	Values
Inlet	Velocity-inlet
Outlet	Pressure-outlet
Drag model	Gidaspow et al. [35]
Heat transfer model	Guun [37]
Granular viscosity	Gidaspow et al. [35]
Granular bulk viscosity	Lun et al. [32]
Restitution	0.9
Wall-condition	No-slip
Friction viscosity	Schaeffer [33]

.

2.4. Geometrical Configuration

In this study, the simulation focused on the FBR prototype as described in reference [38]. The experimental setup primarily consisted of a conveying and metering system, a fluidized bed body, and a pressure signal acquisition system. The schematic representation of the geometric model and meshing can be seen in Figure 1, and the essential parameters for the geometric configuration are detailed in Table 2. The FBR stands at a height of 1.22 m, with a fluidization zone diameter of 0.1524 m. The reaction gases, including SiCl₄, H₂, and HCl, are introduced from the bottom into the FBR, where silicon particles are fluidized, and the exhaust gases exit from the FBR. To optimize the computational efficiency, we simplified the cylindrical reaction zone into a two-dimensional model and employed a quadrilateral meshing approach [19].

2.5. Grid Independence Verification

The computational domain was meshed three times to verify the independence of the grids. The total elements present in the first, second, and third mesh are 11,400, 20,000, and 45,000, respectively. The pressure drop and mean solid volume fraction for all three meshes are shown in Figure 2. Figure 2 shows that the bed pressure drop oscillates in a certain range with the flow time, and the trend of the mean solid volume fraction is basically the same. The second mesh is adequate to predicate the parameters of the bed pressure drop and mean solid volume fraction. Hence, we adopted the mesh of 20,000 elements in the subsequent analysis.

3. Results and Discussion

3.1. Model Validation

In this section, we validate the EGM by comparing simulated values of the maximum bed expansion height (H_max) at various inlet gas velocities with theoretical values obtained from empirical formulas, as illustrated in Figure 3a. The empirical formulas used for comparison are listed in

Table 4. Empirical formulas of fluidization characteristics parameters [39,40].

Parameters	Formulas
$U_{\mathrm{mf}}$	$U_{\mathrm{mf}}=R{e}_{\mathrm{mf}}(\frac{{\mu }_{\mathrm{g}}}{{\rho }_{\mathrm{g}}{d}_{\mathrm{s}}})$	(41)
$U_{\mathrm{mb}}$	$U_{\mathrm{mb}}=0.35\sqrt{gD}$	(42)
$U_{\mathrm{ms}}$	$U_{\mathrm{ms}}=7U_{\mathrm{mf}}\sqrt{gD}+1.6\cdot {10}^{-5}{(60\times {(100D)}^{0.175}-100H_{\mathrm{mf}})}^2$	(43)
$H_{\mathrm{mf}}$	$H_{\mathrm{mf}}=\frac{H(1-\alpha )}{(1-{\alpha }_{\mathrm{mf}})}$	(44)
$\Delta P$	$\Delta P=H_{\mathrm{mf}}(1-{\alpha }_{\mathrm{mf}})({\rho }_{\mathrm{s}}-{\rho }_{\mathrm{g}})g$	(45)
$H_{\max }$	$H_{\max }=H_{\mathrm{mf}}+H_{\mathrm{mf}}\cdot \frac{U-U_{\mathrm{mf}}}{0.5\sqrt{gD}}$	(46)

, and they are sourced from references [39,40]. In these formulas, U_mf, U_mb, and U_ms represent the minimum fluidizing velocity, the minimum bubbling velocity, and the minimum slugging velocity, respectively. H_mf refers to the bed height at minimum fluidization, while H_max represents the maximum bed expansion height. Re_mf corresponds to the Reynolds number at minimum fluidization, and α_mf denotes the bed voidage at minimum fluidization. Figure 3a demonstrates that the simulated values of H_max predicted by the EGM exhibit a linear increase with increasing inlet gas velocities. When the inlet gas velocities are set at 0.15 m/s, 0.30 m/s, and 0.65 m/s, the discrepancies in the H_max between the simulated and theoretical values from the empirical formulas are 1%, 8%, and 9%, respectively. These results indicate that the predicted values from the EGM align well with the theoretical values derived from the empirical formulas. Moving on to Figure 3b, we present a comparison between the simulated and experimental gas temperature profiles within the FBR from the entrance. In accordance with the experimental parameters established by Hsu [38], the inlet temperature is maintained at 300 K, the wall temperature is set to 923 K, and the gas is uniformly introduced into the FBR at an inlet gas velocity of U/U_mf = 5. As depicted in Figure 3b, the simulated gas temperature profile shows a rapid increase from 400 K to approximately 910 K, which closely matches the experimental data. This successful alignment between the simulated and experimental results suggests that the proposed EGM is a viable model for describing the fluidization behavior at various stages within the FBR.

3.2. Flow Pattern and Regime

In Figure 4, the instantaneous solid volume fraction is depicted at various inlet gas velocities. It is important to note that, in this particular analysis, we have not considered chemical reactions and have focused solely on the fluidized bed behavior. It can be seen from Figure 4a that when the inlet gas velocity is 0.15 m/s (which is close to the theoretical initial fluidizing velocity, U_mf, of 0.13 m/s), the bed gradually forms a single “bubble” in the center from the state of static and fixed. The bubble gradually rises to the top of the bed and then disappears. After slight expansion, the top of the bed gradually falls back to a certain stable height. When the inlet gas velocity increases to 0.30 m/s (which is close to the theoretical initial bubbling velocity, U_mb, of 0.26 m/s), as shown in Figure 4b, the movement of the particle phase tends to be disordered. Small and dense bubbles are generated at the bottom of the bed, and semi-circular bubbles are generated near the wall region. The diameter of the bubbles exceeds about 0.6 times the bed diameter, which indicates that the flow pattern of the FBR is close to the state of bubbling. When the inlet gas velocity reaches 0.65 m/s (which is close to the theoretical initial slugging velocity, U_ms, of 0.69 m/s), as shown in Figure 4c, the solid particles flow more intensely and the bubbles disappear at the bottom of the bed. The formation of bubble cutoff layers in the middle of the bed indicates the flow pattern transition to a typical slugging regime. From the detailed analysis above, the developed EGM can be feasible to describe the variation of the flow regime.

Figure 5 illustrates the variation of the bed voidage with the flow time at different inlet gas velocities. It can be seen that the bed voidage increases from 0.40 to 0.60 within the range of the flow time from 0 to 1 s at an inlet gas velocity of 0.15 m/s. At the range of flow time from 1 s to 7 s, the fluctuation in the bed voidage nears the value of 0.60. At this stage, due to the effect of the single “bubble” in the center of the bed (as seen in Figure 4a), the fluidization behavior of the bed is not obvious, but the particle layer has begun to become loose. When the inlet gas velocity reaches 0.30 m/s, the entire bed transitions to the state of bubbling (as seen in Figure 4b), and the movement of the particle phase is obviously affected by the bubble fluctuation. The voidage of the bed increases with the aggregation of the bubble and decreases with the breakage of the bubble. When the inlet gas velocity increases to 0.65 m/s, the entire bed transitions to the state of slugging (as shown in Figure 4c), and the aggregation and breakage of the bubbles become more frequent. Due to the appearance of the bubble cutoff layers within the bed, the voidage of the bed increases rapidly to about 0.70 and then oscillates with the serious fluctuation of the bubbles. Therefore, owing to the transition of the flow pattern and regime, the predicated bed voidage varies with the flow time at different inlet gas velocities.

Figure 6 depicts the distribution of the mean mass fraction of SiCl₄ at different inlet gas velocities. When the inlet gas velocity is 0.15 m/s, there is no obvious fluidization of the bed because the drag force of the gas rising is only slightly greater than the particle gravity; thus, SiCl₄ is mainly distributed in the middle and lower part of the bed. The entire bed transitions to the bubble stage when the inlet gas velocity increases to 0.30 m/s, and the bubble fluctuations effectively promote the mixing of gas and solid particles. Therefore, SiCl₄ is distributed throughout the bed uniformly. However, when the gas velocity is further promoted to 0.65 m/s, SiCl₄ gas flows out of the bed with the solid particles owing to the slugging of the bubbles. This outcome results in less SiCl₄ distribution in the upper part of the bed. Therefore, it can be concluded that the SiCl₄ distribution is more even when the entire bed transitions to the bubbling stage.

3.3. Heat Transfer Characteristics

In this work, two different heat transfer models, namely Guun’s model and the Ranz–Marshall model, were considered to predict the gas–solid heat transfer coefficient. Figure 7 shows a simulated axial distribution of the gas–solid heat transfer coefficients and solid Reynolds number under a transient state. One can see that the axial variation trend of the gas–solid heat transfer coefficient predicted by the two models is consistent with the variation of the solid’s Reynolds number. This outcome demonstrates that the heat transfer between the gas and solids is related to the relative velocity of the two fluids, and the greater the drag force interphase, the higher the heat transfer efficiency. Moreover, the heat transfer coefficient simulated by Guun’s model is higher than that predicted by the Ranz–Marshall model because the Ranz–Marshall model is established for a single sphere particle and underestimates the degree of interphase heat transfer in a dense fluidization system.

Figure 8 displays the axial gas temperature profiles in the FBR from the entrance. It can be seen that when the axial height of the inlet region is Z = 0.05 cm, the average temperature value of the bed section is 369 K, which is close to the set value of the inlet temperature. When the axial height of the bed is Z = 0.2 cm, the average temperature value is 548 K. When the bed axial height is Z = 1 cm, the average temperature rapidly increases to 670 K, which is close to the set value of the FBR’s wall temperature. Our results indicate that the bed quickly reaches the wall temperature of the FBR through a small axial height of 10 mm, and the high heat transfer characteristics of the FBR will facilitate the hydrochlorination reaction performance within the bed.

3.4. Chemical Reaction Performance

In this section, the kinetic model of SiCl₄ hydrochlorination was coupled with the corresponding CFD model to explore the chemical reaction performance. Figure 9 shows the solid volume fraction of the bed during the hydrochlorination process at different inlet gas velocities. Notably, in this scenario, we introduce a comprehensive coupling of the flow field, heat transfer, and chemical reactions. It can be seen from Figure 9a that when the inlet gas velocity is 0.15 m/s, only a single small bubble is generated in the bed, and there is no significant fluctuation in the whole bed, indicating that the SiCl₄ hydrochlorination reaction at this stage mainly occurs in the middle and lower parts of the bed. When the gas velocity at the inlet reaches 0.30 m/s, there is an increase in the expansion height of the bed, leading to a transition of the bed into the bubbling stage. The presence of circulating bubbles within the FBR enhances the mixing uniformity of the SiCl₄ reaction gas and silicon particles, particularly in the central and upper regions of the bed under the bubbling fluidization state. When the gas velocity at the inlet further increases to 0.65 m/s, due to the high gas speed, the particles are easy to fluidize, and the particles are quickly brought into the top of the bed and blown out of the bed. In the slugging stage, the voidage of the bed increases due to the aggregation of large bubbles, which prevents the contact of gas–solid and, thus, affects the reaction performance of the hydrochlorination process.

Figure 10 further provides the mass fraction distribution of the SiHCl₃ products obtained from SiCl₄ hydrochlorination at different inlet gas velocities. One can see that when the inlet gas velocity is 0.15 m/s, the product, SiHCl₃, is mainly distributed in the middle and lower parts of the bed due to the entire bed at the minimum fluidization (as shown in Figure 9a), and the average mass fraction of SiHCl₃ is predicted to be 8.21%. When the inlet gas velocity is 0.30 m/s, the entire bed is in the bubbling stage (as shown in Figure 9b), and the product, SiHCl₃, is mainly distributed in the middle and upper regions of the bed. Due to the good mixing of the gas–solid phase, the SiCl₄ hydrochlorination reaction process is promoted; thus, the average mass fraction of SiHCl₃ obtained is the highest, which is 13.08%. When the inlet gas velocity is 0.65 m/s, the entire bed is in the slugging stage (as shown in Figure 9c), and the reactant of SiCl₄ quickly reaches the top of the bed. Owing to the short gas residence time, the hydrochlorination of SiCl₄ is insufficient; therefore, the average mass fraction of SiHCl₃ decreases to 10.07%.

The conversion rate of the SiCl₄ hydrogenation process (Equation (1)) and the hydrochlorination process (Equation (2)) at different inlet gas velocities was further studied, as shown in Figure 11. One can see that when the FBR is in the bubbling stage, the gas-solid phase exhibits excellent mixing. Consequently, this leads to the highest conversion rate of SiCl₄. By incorporating HCl into the hydrogenation process of SiCl₄, the chemical performance can be further enhanced, as noted in references [27,28]. This enhancement results in a predicted SiCl₄ conversion rate of 28.97% when the reactor is in the bubbling fluidization state.

4. Conclusions

In this study, the fluidization behavior and chemical performance for the hydrochlorination of SiCl₄ within the FBR were examined using the EGM. The effect of the fluidization velocities on the flow regime, heat transfer, and chemical reaction performance was studied. The detailed results are as follows:

• The established EGM can reasonably describe the variation of the flow regime for the hydrochlorination of SiCl₄ within the FBR. The effectiveness of the simulation method and the accuracy of the model were validated through comprehensive comparisons, encompassing a simulated H_max in comparison to theoretical values obtained from empirical formulas and a comparison of the simulated gas temperature profile with data from Hsu’s experiments.
• The bed voidage increased with the flow time, coinciding with transitions in the flow regime within the FBR. Particularly noteworthy was the attainment of a more uniform distribution of SiCl₄ under the bubbling fluidization state.
• The FBR exhibits remarkable heat transfer characteristics, as the reaction gas attains the preset bed temperature within a mere 10 mm from the entrance. Moreover, the presence of circulating bubbles within the FBR greatly enhances the uniform mixing of the SiCl₄ reaction gas with silicon particles.
• When the FBR is in the bubbling stage, it promotes the effective mixing of the gas–solid phases, resulting in the highest conversion rate of SiCl₄. Our predictions revealed a peak concentration of SiHCl₃ at 13.08%, and a SiCl₄ conversion rate of 28.97% was achieved under the conditions of bubbling fluidization.