Modeling of Thermogasdynamics and Thermal Processes of Volatile Halides of Elements of III, IV and VI Groups in Thermal Plasma of RF-IC Discharge at Atmospheric Pressure

Abstract

A mathematical model has been proposed and verified, based on the combination of thermodynamic calculations of equilibrium plasma and CFD-modeling of the plasma-chemical reactor for hydrogen recovery of volatile halides. An experimental study of the hydrogen reduction processes of BCl₃ and SiCl₄ was carried out under conditions of RF-IC (RF inductively coupled) discharge at atmospheric pressure. Based on the experimental results, the developed mathematical model was validated. A series of numerical experiments on parametric study of gas dynamics of high-frequency plasmatron with vortex stabilization of flow for various hydrogen concentrations in argon-hydrogen plasma-forming mixture was carried out. Thermodynamic calculations and CFD modeling were made for input into the reaction chamber of chemically active gas mixtures (BCl₃ + H₂), (BF₃ + H₂), (SiCl₄ + H₂), (SiF₄ + H₂), (MoF₆ + H₂), (WF₆ + H₂), (BCl₃ + CH₄ + H₂), (BF₃ + CH₄ + H₂), (SiF₄ + CH₄ + H₂), (MoF₆ + CH₄ + H₂) and (WF₆ + CH₄ + H₂). The values for the predictive efficiency and the degree of conversion of target products were obtained.

Appendix

Mathematical Formulation of the Problem in the Zone of Plasmatron I (Fig. 12)

The flow of plasma-forming gas in plasmatron I (Fig. 12) while solving the first stage of the gas-dynamics problem of RF-IC-plasmatron with vortex generator is described by the system of Eqs. (6–10) and an additional relation (11) with the corresponding boundary and initial conditions. Due to the lack of physical symmetry of the flow in plasmatron with vortex generator, the thermogasdynamic component of the problem was solved in a three-dimensional formulation. The axial cross section of the calculated field of the problem is represented in Fig. 15a. Figure 15b shows the photo of the studied RF-IC-plasmatron. The field of solving the problem includes plasmatron (zone A Fig. 15a) with the walls of the final thickness (zone B Fig. 15a), the inductor (zone C Fig. 15a) and the outer part of the air space (zone D Fig. 15a) adjacent to the device. Excluding the bias current, the spiral inductor can be represented as a system of four cylindrically symmetric parallel rings.

$$ \frac{\partial \rho }{{\partial t}} + div(\rho \overline{V}) = 0, $$

(6)

$$ \frac{{\partial (\rho \overline{V})}}{\partial t} + div\left( {\rho \overline{V} \times \overline{V}} \right) = - gradp + div\left( {\mu_{t} \left( {grad\overline{V} + (grad\overline{V})^{T} - \frac{2}{3}\delta div\overline{V}} \right)} \right) + F_{e - m} , $$

(7)

$$ \frac{{\partial (\rho \,c_{p} T)}}{\partial t} + div(\rho \,c_{p} T \times \overline{V}) = div(\lambda_{t} \nabla T) + \frac{\partial p}{{\partial t}} + S_{J} - u_{rad} . $$

(8)

In the equation of continuity (6), movement (7) and energy (8) ρ is the density of the gas mixture Ar:H₂, V is the velosity, p is the pressure, $c_{p}$ is the specific heat capacity, T is the temperature, $\mu_{e} ,\lambda_{e}$ is the effective viscosity and thermal conductivity, $\mu_{e} = \mu + \mu_{t}$, $\lambda_{e} = \lambda + \lambda_{t}$, $\mu ,\lambda$ are the laminar viscosity and thermal conductivity, $\mu_{t} ,\lambda_{t}$ are the turbulent (vortex) viscosity and thermal conductivity, $F_{e - m} = 0.5 \cdot \sigma [\vec{E} \times \vec{B}]$ is the Lorentz power, k_b is the Boltzmann constant, e is the elementary electric discharge, $S_{J} = 0.5 \cdot \sigma \vec{E}^{2}$ is the power of joule heating, $u_{rad}$ is the radiation energy lost, σ is the plasma electrical conductivity, $\vec{E}$ is the electric field strength, $\vec{B}$ is the magnetic field strength, the upper index ^T is the transposition operator of the matrix.

While modeling the gasdynamics of plasmatrons with straightforward gas supply, as a rule, a laminar model is used which is justified by the results of physical experiments. However, preliminary calculations showed that for adequate description of the movement in the investigated plasmatron with vortex generator, the turbulence should be taken into account both near walls and in the reactor volume. The system of Eqs. (6–8) is closed by a semi-empirical k-ω SST model of turbulence of Menter [24]. The selection of the model is due to its popularity for solving a wide range of engineering problems and the potential possibility of joint use of the model of laminar-turbulent transition of γ-Reθ [25].

The values of density, viscosity, specific heat, thermal conductivity, electrical conductivity and the specific energy of the radiation of the gas mixture in the model are set by the temperature functions [26].

To find the values $\vec{E}$ и $\vec{B}$ it is necessary to find the solution to electromagnetic problem. For the process under consideration, the Maxwell equations system can be reduced to the Poisson equation for the complex vector of magnetic potential $\vec{A}$. In this case, in the cylindrical coordinates (rθ z) only the angular component of the vector $\vec{A}$ is different from zero, and then the equations for the real $A_{R\theta }$ and imaginary $A_{I\theta }$ parts of the magnetic potential are:

$$ \nabla^{2} A_{R\theta } + \mu_{0} \sigma *\omega A_{I\theta } - \frac{{A_{R\theta } }}{{r^{2} }} + \frac{{\mu_{0} J_{k} }}{{\pi a^{2} }} = 0, $$

(9)

$$ \nabla^{2} A_{I\theta } + - \mu_{0} \sigma *\omega A_{R\theta } - \frac{{A_{I\theta } }}{{r^{2} }} = 0. $$

(10)

where $\mu_{0}$ is the magnetic constant, $\omega = 2\pi f$, J_к is the inductor current (different from zero only in zone III and is determined from the conditions of the physical experiment), and a is the radius of the cross section of the inductor. The value of the electrical conductivity of the medium is:

$$ \sigma * = \left[ \begin{gathered} \sigma ,\;{\text{in zone I; }} \hfill \\ {\text{0, in zones II, III,IV}}{.} \hfill \\ \end{gathered} \right. $$

The strength of electrical and magnetic fields, included in Eqs. (7) and (8), have the form:

$$ B_{r} = - \frac{{\partial A_{\theta } }}{\partial z},\quad B_{\theta } = 0,\quad B_{z} = - \frac{1}{r}\frac{\partial }{\partial r}(rA_{\theta } ),\quad E_{r} = 0,\quad E_{\theta } = - i\omega A_{\theta } ,\quad E_{z} = 0. $$

(11)

The system of Eqs. (6–8) is solved in zone A of the computational domain (Fig. 15a) and is complemented by boundary conditions which are set in accordance with experimental data. The plasma-forming gas of room temperature with the given volume consumption and Ar/H₂ ratio is supplied to the input boundary of plasmatron, the zero excessive pressure is set on the output boundary, the adhesion condition is set on the solid walls while on the upper wall of plasmatron the temperature is 293 K, on the side wall 1450 K.

The boundary conditions of Eqs. (9) and (10) are: continuity $A_{\theta }$ on the internal boundaries of the computational domain; equality to zero of the vector of magnetic potential at the outer boundaries of the environment and walls of the channels; and equality to zero of the flow of magnetic potential at the outer boundaries of plasmatron.

It is assumed in the model that zones B, C, D of the computational domain (Fig. 15) are a fictitious solid medium with permanent properties. Due to the fact that the temperature of plasmatron walls is considered to be specified, there is no need to solve the thermodynamic problem in zones B, C and D. Only Eqs. (9) and (10) are solved in these zones. This assumption does not affect the correctness of the solution of the thermogasdynamic problem in zone A, but leads to a significant optimization of computational resources.

Mathematical Formulation of the Problem in Zones II and III of the Plasma Chemical Reactor (Fig. 12)

Consider, first of all, the case of input to the reactor of chemically active mixtures for the synthesis of heavy metals (molybdenum and tungsten) and their carbides. The investigation is carried out on the homogeneous flow of argon-hydrogen mixture, gas-carrier H₂, gas-reagent, powdered particles of target products, and, optionally, in the case of the synthesis of carbides, methane and carbon, formed in significant amount that is proved by thermodynamic calculations. Note that test calculations indicate a violation of the homogeneity of the flow under study with radius of the convertible particles more than 10 μm which is by two orders of magnitude larger of the assumed dimensions of fine particles (about 200 nm) obtained in the current RF-IC-plasmatron. The homogeneous turbulent flow of the mixture is described by the system of equations of continuity, momentum and energy. The heat exchange is carried out by convection and thermal conductivity. In the continuity equation for the corresponding components of the gas mixture, the mass sources have been added (for carbon, molybdenum or tungsten, carbides of molybdenum or tungsten) and discharge (for gas-reagent and methane). In this case, the mass source terms in each cell of the computational domain are determined on the basis of the dependence of the concentration of the corresponding substance on temperature obtained as a result of thermodynamic calculations. Based on test calculations, the description of the movement of media mixture in the reactor is carried out with the following assumptions:

• The mass transfer of molybdenum fluoride or tungsten fluoride, hydrogen and methane (in the case of the synthesis of carbides) into hydrogen fluoride and target reaction products–molybdenum or tungsten and, optionally, their carbides is considered. Note that according to thermodynamic calculations, the formation of other solid fluorides does not occur in the studied temperature range. The formation of atomic hydrogen, fluorine as well as carbon, in the case of the synthesis of carbides, occurs at temperatures above 3000 K (Figs. 5, 6, 10, and 11), and is not taken into account in the gas-dynamic model due to small volume of a part of the computational domain where temperature exceeds 3000 K (< 7% of the total volume);

• The formation of fine condensed phase of target products in the model of gas-dynamic flow is thermodynamically determined by the concentration in this cell of the computational domain of gas-reagents while the condensation process is not modeled in the problem;

• Excluding chemical kinetics, it is assumed that the thermodynamically determined composition of the mixtures, represented in Figs. 1, 2, 3, 4, 5, 6, 7,,8, 9, 10, and 11, is formed depending on the temperature in an infinitely small period of time.

The computational domain includes the reaction zone of plasma chemical reactor (Fig. 12, II) and the cooled channel (Fig. 12, III). The turbulent movement of viscous compressible perfect gases and fine particles in inhomogeneous temperature field is considered taking into account heat exchange by thermal conductivity and convection.

The flow of gases is described by the system of equations of momentum (12), energy (13) and continuity (14):

$$ \frac{{\partial (\varphi_{i} \rho \overline{V})}}{\partial t} + div(\rho \overline{V} \times \overline{V}) = - gradp + div(\mu_{e} (grad\overline{V} + (grad\overline{V})^{T} - \frac{2}{3}\delta div\overline{V})), $$

(12)

$$ \frac{{\partial (\varphi_{i} \rho \,c_{p} T)}}{\partial t}div(\rho \,c_{p} T \times \overline{V}) = div(\lambda_{e} \nabla T) + \frac{\partial p}{{\partial t}}, $$

(13)

$$ \frac{{\partial (\varphi_{i} \rho )}}{\partial t} + div(\varphi_{i} \rho_{i} \overline{V}) + S_{i} = 0,\quad i = \overline{1,8.} $$

(14)

In (12–14) index i = 1 corresponds to the argon–hydrogen plasma-forming mixture, i = 2–gas-carrier of hydrogen H₂; i = 3–gas-reagent (MoF₆ or WF₆); i = 4–methane CH4 (in the case of supplying it for hydrogen recovery of carbides), i = 5–hydrogen fluoride HF (g), i = 6–molybdenum Mo or tungsten W, i = 7–molybdenum carbide Mo₂C or tungsten carbide WC; i = 8–molybdenum carbide Mo₃C₂ (when applying gas-reagent MoF₆); ϕ_i, ρ_i,– volume fraction and density of the corresponding medium, while: $\sum\nolimits_{i = 1}^{8} {\varphi_{i} } = 1;\rho = \sum\nolimits_{i = 1}^{8} {\varphi_{i} \rho_{i} }.$

The system of Eqs. (12–14) is closed by the k–ω SST turbulence model which describes the vortex of the flow in the near-wall region and, according to the results of test calculations, provides a high speed of the task convergence.

In Eq. (14) S_i are the mass source terms of the corresponding substance in the unit of volume while in the case of supplying chemically active mixtures (MoF₆ + H₂) and (WF₆ + H₂) $S_{4} = S_{7} = S_{8} = 0.$ to the reactor.

As an example, in expressions (15), (16) and (17) the type of source member for the target reaction products is given –Mo and Mo₂C, Mo₃C₂, respectively, (while supplying CH₄ to the reactor):

$$ S_{6(Mo)} = \varphi_{3} \rho_{3} {\raise0.7ex\hbox{${\dot{n}_{6} (T)}$} \!\mathord{\left/ {\vphantom {{\dot{n}_{6} (T)} {n_{{_{3} }}^{0} }}}\right.\kern-\nulldelimiterspace} \!\lower0.7ex\hbox{${n_{{_{3} }}^{0} }$}} \cdot {\raise0.7ex\hbox{${M_{6} }$} \!\mathord{\left/ {\vphantom {{M_{6} } {M_{3} }}}\right.\kern-\nulldelimiterspace} \!\lower0.7ex\hbox{${M_{3} }$}}. $$

(15)

Hereinafter $\dot{n}(T) = \frac{\partial n(T(t))}{{\partial t}}$ is the change in the amount of substance per unit of time; n⁰_i(T) is the initial amount of the substance given in the thermodynamic calculation; M_i is the molar mass. The index i corresponds to the notations in Eqs. (12–14). For the convenience of perception the substance is added to the mass source index in brackets to which this source refers.

$$ \left[ \begin{gathered} if\;\varphi_{3} < 2\varphi_{4} :\quad S_{7(Mo2C)} = \varphi_{3} \rho_{3} \frac{{\dot{n}_{7} (T)}}{{n_{{_{3} }}^{0} }} \cdot \frac{{M_{7} }}{{(M_{3} + M_{4} )}} + \tfrac{1}{2}\varphi_{3} \rho_{4} \frac{{\dot{n}_{7} (T)}}{{n_{{_{4} }}^{0} }} \cdot \frac{{M_{7} }}{{(M_{3} + M_{4} )}}, \hfill \\ if\;\varphi_{3} > 2\varphi_{4} :\quad S_{7(Mo2C)} = 2\varphi_{4} \rho_{3} \frac{{\dot{n}_{7} (T)}}{{n_{{_{3} }}^{0} }} \cdot \frac{{M_{7} }}{{(M_{3} + M_{4} )}} + \varphi_{4} \rho_{4} \frac{{\dot{n}_{7} (T)}}{{n_{{_{4} }}^{0} }} \cdot \frac{{M_{7} }}{{(M_{3} + M_{4} )}}. \hfill \\ \end{gathered} \right. $$

(16)

$$ \left[ \begin{gathered} if\;\varphi_{3} < 1.5\varphi_{4} :\quad S_{8(Mo3C2)} = \varphi_{3} \rho_{3} \frac{{\dot{n}_{8} (T)}}{{n_{{_{3} }}^{0} }} \cdot \frac{{M_{8} }}{{(M_{3} + M_{4} )}} + \tfrac{2}{3}\varphi_{3} \rho_{4} \frac{{\dot{n}_{8} (T)}}{{n_{{_{4} }}^{0} }} \cdot \frac{{M_{8} }}{{(M_{3} + M_{4} )}}, \hfill \\ if\;\varphi_{3} > 1.5\varphi_{4} :\quad S_{8(Mo3C2)} = 1.5\varphi_{4} \rho_{3} \frac{{\dot{n}_{8} (T)}}{{n_{{_{3} }}^{0} }} \cdot \frac{{M_{8} }}{{(M_{3} + M_{4} )}} + \varphi_{4} \rho_{4} \frac{{\dot{n}_{8} (T)}}{{n_{{_{4} }}^{0} }} \cdot \frac{{M_{8} }}{{(M_{3} + M_{4} )}}. \hfill \\ \end{gathered} \right. $$

(17)

The values of density, viscosity, the specific heat and thermal conductivity of the plasma-forming gas, H₂ and CH₄ in the model are defined by the temperature functions [26,27,28]. For MoF₆, HF(g), as well as for Mo, the specific heat capacity was determined using the NIST database $(c_{p} (T) = \frac{1000}{M}(A + B\frac{T}{1000} + C(\frac{T}{1000})^{2} + D(\frac{T}{1000})^{3} + E(\frac{T}{1000})^{ - 2} )$, where A, B, C, D, E are the coefficients of the polynomial equation c_p(T) taken from the NIST database), viscosity and thermal conductivity were calculated by the formulas [28]:

$$ \mu (T) = 0.00063\frac{{M^{1/2} p_{cr}^{2/3} }}{{T_{cr}^{1/6} }} \cdot \frac{{\left( {{\raise0.7ex\hbox{$T$} \!\mathord{\left/ {\vphantom {T {T_{cr} }}}\right.\kern-\nulldelimiterspace} \!\lower0.7ex\hbox{${T_{cr} }$}}} \right)^{3/2} }}{{{\raise0.7ex\hbox{$T$} \!\mathord{\left/ {\vphantom {T {T_{cr} }}}\right.\kern-\nulldelimiterspace} \!\lower0.7ex\hbox{${T_{cr} }$}} + 0.8}}, $$

(18)

$$ \lambda (T) = \frac{{c_{p} (T) \cdot \mu (T)}}{\Pr }, $$

(19)

where p_cr is critical pressure (atm); T_cr is critical temperature (K); M is the molar mass (g/mol), $\Pr = \frac{{c_{p} }}{{1.204 \cdot c_{p} (T) + 1.47}}$. The target products of reaction are assumed as incompressible media, their densities are defined equal to the densities of the corresponding materials at room temperature. Temperature dependence of thermal conductivity Mo is taken from [29]. Temperature dependences of the specific heat and thermal conductivity coefficient for Mo₂C in the literature are absent; in the model the thermophysical properties of this material are set by constants corresponding to the room conditions [30]. The data on the thermophysical properties of Mo₃C₂ are not given in open sources, they are considered to be equal to the corresponding characteristics for Mo₂C. While developing the model, a parametric study was carried out on the effect of the physical properties of the reaction products on the gas-dynamic pattern in the reactor. The admissibility of applying these assumptions about the thermophysical properties of the fine condensed phase is established.

In the case of modeling the input of the mixtures of (BCl₃ + H₂), (BF₃ + H₂), (SiCl₄ + H₂), (SiF₄ + H₂), (SiF₄ + CH₄ + H₂), (BF₃ + CH₄ + H₂) into reactor, the simplified model is considered. It is assumed that the powdered particles B, Si, B₄C and SiC due to their low concentrations, expected on the basis of the thermodynamic calculations, do not affect the flow gasdynamics in the reactor. A homogeneous turbulent flow of the mixture is considered which includes only argon–hydrogen mixture, gas-carrier, gas reagent (one of the following compounds depending on the type mixture supplied to input 4 (Fig. 12): BCl₃, BF₃, SiCl₄, SiF₄) and, optionally, methane and carbon. To describe the distribution of powdered particles of target products in the gas mixture, the equation of diffusion was used with the addition of an appropriate source term that defines the formation of B, Si, B₄C and SiC according to thermodynamic calculations. The diffusion coefficients of the particles of boron, silicon and their carbides in the carrying gas stream are estimated by the Stokes–Einstein formula. The radius of particles was assumed equal to 200 nm.

$$ \frac{\partial (\rho c)}{{\partial t}} + div(\rho \overline{\user2{V}}c) = div(\sum\limits_{i = 1}^{5} {\rho_{i} D_{i} } \cdot grad\,c) + S_{reag} , $$

(18)

where c и D_i are the mass concentration and diffusion coefficient in the corresponding medium of the powder particles of the target product.

Boundary conditions are specified in accordance with the experimental conditions of the conversion process in the plasma-chemical reactor with the vortex stabilization of the flow. In the case of modeling the synthesis of carbides, the chemically active mixture of H₂, gas-reagent and CH₄ of room temperature is supplied to the inputs of tubes 4 (Fig. 12). While carrying out calculations for the reduction B, Si, Mo and W, the CH₄ is not applied to the computational domain. The values for the volumetric flow and the ratio of components at the inlet 4 (Fig. 12) for each studied chemically active mixture are set equal to the corresponding values for which thermodynamic calculations were carried out. The plasma-forming gas with the velocity fields, specified on the boundary, and temperature fields, defined from solving the problem in the plasmatron zone I is supplied to the inflow boundary separating plasmatron zone I and reaction chamber II. The zero excessive pressure is set at the outflow boundary 2 (Fig. 12). On all solid walls, the condition of adhesion is set and the temperatures of walls correspond to the experimental ones. The boundary conditions of diffusion equations for particles B, Si, B₄C and SiC are zero concentrations of particles at the inflow 1, 4 (Fig. 12) and the absence of particle flow through the solid walls of the computational domain. The physical properties of the components of the gas mixture are set by the functions of temperature. For BCl₃, BF₃, SiCl₄, SiF₄ specific heat capacity was determined using the NIST database, viscosity and thermal conductivity were calculated using formulas (18) and (19).