Multiscale CFD Simulation of Multiphase Erosion Process in a Connecting Pipe of Industrial Polycrystalline Silicon Unit

Abstract

Severe erosion phenomena often occur in industrial polycrystalline silicon units, leading to hydrogen leakage accidents and affecting long-term operation. It is favorable to use a computational fluid dynamics (CFD) simulation with the dense discrete phase model (DDPM) and the sub-grid energy-minimization multi-scale (EMMS) drag model to improve the prediction accuracy of complex multiphase erosion phenomena in a connecting pipe of an industrial polycrystalline silicon unit. Furthermore, the effect of droplet the specularity coefficient on boundary conditions is thoroughly considered. The predicted erosion behaviors are consistent with industrial data. The effects of operations parameters were discussed with three-dimensional CFD simulation, including droplet size and hydrogen volume fraction on erosion behaviors. The results indicated that the non-uniform multiphase erosion flow behavior near the wall can be simulated accurately with the EMMS drag model in a coarse mesh. A suitable droplet specularity coefficient such as 0.5 can also improve the accuracy of erosion position. Small liquid droplets, such as those of 30 μm size, will follow the gas phase better and have a lower erosion rate. The inertia effect of large droplets, such as those of 150 μm size, plays a dominant role, resulting in obvious erosion on the elbow walls. The erosion range and thinning rate enlarge with the increase in hydrogen volume fraction. A few silicon solid particles, such as 0.01% volume fraction, change local flow behaviors and probably cause the variation of local erosion positions. The process of erosion deformation first circumferentially extended and then accelerated at the local center position deeper.

1. Introduction

Polycrystalline silicon is a critical raw material for new high-tech industries, such as semiconductors, electronic information, and solar photovoltaics [1]. Siemens technology of polycrystalline silicon production has become a primary process due to low cost and high product quality, the products of which occupied 80% of the market in recent years [2,3]. Here, an essential device of silicon tetrachloride recovery is directly related to the polysilicon yields and the amounts of saving energy consumption and cost, where silicon tetrachlorides are transformed into trichlorosilane via hydrogenation reaction [4]. However, the pipes often erode in the device, especially the connecting pipes from the gas–liquid static mixer to the silicon tetrachloride evaporator. The SA 106 Gr.B steel connecting pipes are eroded because of high-velocity (30 m/s to 60 m/s) gas droplet flow, including hydrogen and liquid silicon tetrachloride. Therein, the hydrogen volume fractions are more than 90%, and the liquid droplets’ diameters of silicon tetrachloride are in the range of 30 mm to 200 mm. Additionally, the erosion is further accelerated due to the invasion of a few silicon particles derived from the silicon tetrachloride fractionating tower (solid particles volume fraction about 0.01% and the diameters about 30–50 mm). These pipe erosions often cause pipe breaks, leading to hazardous medium leakage and even hydrogen gas explosion accidents [5,6]. That severely affects the profits of polycrystalline silicon units. Accurately predicting the erosion position and thinning rate of pipes in the device is an essential precondition of safety protection and long-time operation of the silicon tetrachlorides recovery device.

Currently, many researchers focus on the erosion simulations of gas–solid [7,8,9] or liquid–solid [10,11,12] in the pipelines of oil gas exploiting, transporting, and production processes. Many effect factors [13,14,15,16,17,18] were presented to compare the variation of erosion rate, including particle mass fraction [19,20], shape [21], and diameter distribution [22,23], flow velocity [24,25], material hardness [23], and so on. Therein, the erosion rate is in an exponential increase relationship with the flow velocity. The other factors are directly proportional to erosion rate to a certain degree. Some studies [26,27] have investigated the gas–liquid droplet erosion phenomenon only for the simple structures of elbow and tee, and few accounted for the influences of non-uniform gas–liquid droplet flow on the erosion position and rate under actual operation conditions in complex industrial pipe. Moreover, the gas–liquid–solid erosion phenomenon is less discussed and does not predict the actual erosion position and rate under complex industrial conditions.

It is difficult to disclose the erosion processes in the complex structure and the operation condition with experiments. Most studies to predict and analyze the complex erosion phenomenon used the Eulerian–Lagrangian methods, such as the discrete phase model (DPM) [28,29], dense discrete phase model (DDPM) [19,30], and discrete element method (DEM) [31,32,33]. DPM in Ansys Fluent^® is usually used to trace a low-volume fraction of particles but neglects particle–particle interactions. DDPM in Ansys Fluent^® further considers the particle volume fraction and particle–particle collision process described by the hard sphere model, but particle–particle collision forces are calculated by the kinetic theory of granular flow (KTGF). DEM can disclose the whole unsteady particle–particle collision process with the soft sphere model but needs massive computing resources [34]. The DEM is less used to simulate the erosion phenomenon in complex industrial pipes.

Despite many [35,36,37] having studied multiphase erosion behaviors using DPM or DDPM methods, fewer have paid attention to the effects of sub-grid inter-phase force and discrete phase boundary condition on the multiphase flow distributions and the erosion positions. In the multiphase non-uniform erosion process, the drag force plays a dominant role in the inter-phase force. The traditional drag models derive from experiment correlation based on uniform flow structure [38], such as Gidapsow [39,40], spherical [41], and non-spherical models [42], which are more used for the low-particle Reynolds numbers. However, gas–liquid droplet or gas–liquid–solid flow exists in the sub-grid heterogeneous structure of clusters due to droplets or particles aggregation, leading to an excessively uniform flow field and overestimated drag force in coarse mesh [43,44,45,46,47,48], especially in the dense particle zones. To disclose the non-uniform flow structure and consider the effects on multiphase flow, fine-grid simulation with a resolution down to less than 10 times the particle size [49] is required, so it is unaffordable for simulations of large-scale industrial pipelines. To resolve the above problems, many researchers proposed meso-scale drag models considering sub-grid structures in a coarse mesh, such as the filtered drag model [50,51,52] and the energy-minimization multi-scale (EMMS) drag models [38,46,53,54,55,56]. The former is based on the average statistics results with fine grid or direct numerical simulation approaches and a correlation model to consider the scale effects on flow behavior, such as the correlation between mesoscale and particle scale in a sub-grid. The latter is based on the governing principle of compromise in competition for mesoscale phenomena that includes three mechanisms [38,55], i.e., the fluid-dominant (FD), the particle-dominated (PD), and the compromise in competition between the dominant mechanism of particles (PFC). It is a coupling model to disclose the interactive relationship between mesoscale, microscale, and macroscale; furthermore, its effects of the mesoscale structures on the gas–solid or gas–liquid two-phase flow in large-scale fluidized bed reactors are investigated while demonstrating that it is insensitive to the grid size in previous literatures [43,57,58,59,60,61,62]. However, fewer studies have the effects of interphase drag force on gas-liquid-solid flow on wall erosion.

Additionally, the normal and tangent discrete phase reflection coefficients were used to describe the momentum change after particles rebounded off the wall boundary in previous literature reviews [63], while the droplet or particle specularity coefficient of the Johnson and Jackson boundary condition [64] in DDPM is also a considerable factor to controlling the amount of lateral momentum transfer between particles and the wall. Many factors, such as operating conditions, wall roughness, particle properties, and geometric configurations, are found to affect specularity coefficients [63,64]. Several studies have indicated that the flow field is sensitive to the value of specularity coefficient in the gas–solid two-fluid model (TFM) [43,65,66,67,68,69] and they have proposed a model to estimate the variable specularity coefficient and indicated that a constant value is also viable provided that it is properly set [70,71]. However, less discussion has taken place concerning its effects on the simulation with Eulerian–Lagrangian models, such as computational fluid dynamics with dense discrete phase models (CFD-DDPM).

In this study, to improve predicted accuracy, the CFD-DDPM method coupled with an energy minimization multi-scale (EMMS) drag model was proposed to simulate the gas–liquid and gas–liquid–solid non-uniform erosion process and the gas-droplet erosion phenomena under industrial operation conditions in a 3D industrial polycrystalline silicon production pipeline. We focused on the effects of interphase drag and wall specularity coefficient of droplets on erosion behaviors. Then, the operation conditions were analyzed and improved in terms of the simulation approach, while the erosion deformation process was predicted with a moving deforming mesh solver. Those results will be conducive to guiding the safe and long-term operation of industrial polycrystalline silicon units.

2. Simulation Description

2.1. Geometry of an Industrial Pipe of Polycrystalline Silicon Production

Figure 1 presents the geometry of a connecting pipe in an industrial silicon tetrachloride hydrogenation reactor, where severe erosion phenomena often appear in the A and B elbows and even cause hydrogen leakage accidents. The pipe material is SA106 Gr. B steel. This pipe consists of four 90° elbows and one 45° elbow with a bend diameter ratio of 1.5 at the C position, whose diameter and thickness are Φ300 and 8.8 mm, respectively. The lengths of every section are listed in Figure 1. To improve the quality of the mesh, we used the hexahedral mesh drawn with ANSYS ICEM CFD 2019 R1^® software, opted for an O-shape mesh of the circular cross-section, and refined the meshes at the elbows. To ensure accuracy and improve the efficiency of the calculation, 1,091,305 grids were selected after the optimization of grid independence, as shown in Appendix A.

2.2. Governing Equations and Constitutive Relations

The dense discrete particle model (DDPM) in Ansys Fluent^® software was adopted to solve the interaction between continuous (i.e., hydrogen gas) and discrete phases (i.e., liquid silicon tetrachloride droplets and solid silicon particles), whose relevant governing equations are summarized in

Table 1. Summary of governing equations for gas-particles flow in DDPM.

Continuity equation: $\frac{\partial }{\partial t}\left({\alpha }_g{\rho }_g\right)+\nabla \cdot \left({\alpha }_g{\rho }_g{\overrightarrow{u}}_g\right)=0$ $\frac{\partial }{\partial t}\left({\alpha }_{\mathrm{p}}{\rho }_{\mathrm{p}}\right)+\nabla \cdot \left({\alpha }_{\mathrm{p}}{\rho }_{\mathrm{p}}{\overrightarrow{u}}_{\mathrm{p}}\right)=0$ Momentum equation: $\frac{\partial }{\partial t}\left({\alpha }_{\mathrm{g}}{\rho }_{\mathrm{g}}{\overrightarrow{u}}_{\mathrm{g}}\right)+\nabla \cdot \left({\alpha }_{\mathrm{g}}{\rho }_{\mathrm{g}}{\overrightarrow{u}}_{\mathrm{g}}{\overrightarrow{u}}_{\mathrm{g}}\right)=-{\alpha }_{\mathrm{g}}\nabla p+\nabla \cdot {\overline{\overline{\tau }}}_{\mathrm{g}}+{\alpha }_{\mathrm{g}}{\rho }_{\mathrm{g}}\overrightarrow{g}+\beta \left({\overrightarrow{u}}_{\mathrm{p}}-{\overrightarrow{u}}_{\mathrm{g}}\right)$ $\frac{\partial }{\partial t}\left({\alpha }_{\mathrm{p}}{\rho }_{\mathrm{p}}{\overrightarrow{u}}_{\mathrm{p}}\right)+\nabla \cdot \left({\alpha }_{\mathrm{p}}{\rho }_{\mathrm{p}}{\overrightarrow{u}}_{\mathrm{p}}{\overrightarrow{u}}_{\mathrm{p}}\right)=-{\alpha }_{\mathrm{p}}\nabla p+\nabla \cdot {\overline{\overline{\tau }}}_{\mathrm{p}}+{\alpha }_{\mathrm{p}}{\rho }_{\mathrm{p}}\overrightarrow{g}+\beta \left({\overrightarrow{u}}_{\mathrm{p}}-{\overrightarrow{u}}_{\mathrm{g}}\right)$ Granular motion equation: $\frac{d{\overrightarrow{u}}_{\mathrm{p}}}{dt}=\frac{\left({\rho }_{\mathrm{p}}-{\rho }_{\mathrm{g}}\right)\overrightarrow{g}}{{\rho }_{\mathrm{p}}}+\frac{{\rho }_{\mathrm{g}}}{{\rho }_{\mathrm{p}}}{\overrightarrow{u}}_{\mathrm{g}}\nabla {\overrightarrow{u}}_{\mathrm{g}}+\beta \left({\overrightarrow{u}}_{\mathrm{g}}-{\overrightarrow{u}}_{\mathrm{p}}\right)+{\overrightarrow{F}}_{\mathrm{KTGF}}$ where ${\overrightarrow{F}}_{\mathrm{KTGF}}=-\frac{m_{\mathrm{p}}}{{\rho }_{\mathrm{p}}}\nabla {\overline{\overline{\tau }}}_{\mathrm{p}}$, ${\mathrm{Re}}_{\mathrm{p}}=\frac{{\rho }_{\mathrm{p}}{d}_{\mathrm{p}}\left|{\overrightarrow{u}}_{\mathrm{p}}-{\overrightarrow{u}}_{\mathrm{g}}\right|}{{\mu }_{\mathrm{g}}}$. Gidaspow drag coefficient [39]: $\begin{array}{l}\beta =\left\{\begin{array}{l}\frac{3}{4}{C}_{\mathrm{d}}\frac{{\alpha }_{\mathrm{p}}{\alpha }_{\mathrm{g}}{\rho }_{\mathrm{g}}\left|{\overrightarrow{u}}_{\mathrm{g}}-{\overrightarrow{u}}_{\mathrm{p}}\right|}{d_p}{\alpha }_{\mathrm{p}}^{-2.65},{\alpha }_{\mathrm{g}}\le 0.8\\ 150\frac{{\alpha }_{\mathrm{p}}^2{\mu }_{\mathrm{g}}}{{\alpha }_{\mathrm{g}}{d}_{\mathrm{p}}^2}+1.75\frac{{\rho }_{\mathrm{g}}{\alpha }_{\mathrm{p}}}{d_{\mathrm{p}}}\left|{\overrightarrow{u}}_{\mathrm{g}}-{\overrightarrow{u}}_{\mathrm{p}}\right|,{\alpha }_{\mathrm{g}}>0.8\end{array}\right.\\ C_{\mathrm{d}}=\left\{\begin{array}{l}\frac{24}{{\mathrm{Re}}_{\mathrm{p}}}\left[1+0.15{\mathrm{Re}}_{\mathrm{p}}{}^{0.687}\right],{\mathrm{Re}}_{\mathrm{p}}\le 1000\\ 0.44,{\mathrm{Re}}_{\mathrm{p}}>1000\end{array}\right.\end{array}$ EMMS drag coefficient [46]: $\begin{array}{l}\beta =\frac{3}{4}{C}_{\mathrm{d}}\frac{{\alpha }_{\mathrm{p}}{\alpha }_{\mathrm{g}}{\rho }_{\mathrm{g}}\left|{\overrightarrow{u}}_{\mathrm{g}}-{\overrightarrow{u}}_{\mathrm{p}}\right|}{d_{\mathrm{p}}}{\alpha }_{\mathrm{g}}^{-2.65}{H}_{\mathrm{d}}\\ C_{\mathrm{d}}=\left\{\begin{array}{l}\frac{24}{{\mathrm{Re}}_{\mathrm{p}}}\left[1+0.15{\mathrm{Re}}_{\mathrm{p}}^{0.687}\right],{\mathrm{Re}}_{\mathrm{p}}\ge 1000\\ 0.44,{\mathrm{Re}}_{\mathrm{p}}<1000\end{array}\right.\end{array}$ where the heterogeneity index, Hd, is provided in Appendix A. Spherical drag coefficient [41]: $\begin{array}{l}\beta =\frac{3{\mu }_{\mathrm{p}}{C}_{\mathrm{d}}{\mathrm{Re}}_{\mathrm{p}}}{4{\rho }_{\mathrm{p}}{d}_{\mathrm{p}}^2}\\ C_{\mathrm{d}}=\frac{24}{{\mathrm{Re}}_{\mathrm{p}}}\left[1+0.1806{\mathrm{Re}}_{\mathrm{p}}^{0.6459}\right]+\frac{0.4251{\mathrm{Re}}_{\mathrm{p}}}{{\mathrm{Re}}_{\mathrm{p}}+6880.95}\end{array}$ Granular temperature equation (Algebraic): $0=\left(-p_{p}I+{\tau }_p\right): \nabla u_p-{\gamma }_p-3\beta \cdot {\Theta }_p$ Stress–strain tensors: ${\overline{\overline{\tau }}}_{\mathrm{g}}={\mu }_{\mathrm{g}}\left[\nabla {\overrightarrow{u}}_{\mathrm{g}}+{\left(\nabla {\overrightarrow{u}}_{\mathrm{g}}\right)}^T\right]-\frac{2}{3}{\mu }_{\mathrm{g}}\nabla \cdot {\overrightarrow{u}}_{\mathrm{g}}\overline{I}$ ${\overline{\overline{\tau }}}_{\mathrm{p}}=-p_{\mathrm{p}}\overline{I}+{\alpha }_{\mathrm{p}}{\mu }_{\mathrm{p}}\left(\nabla {\overrightarrow{u}}_{\mathrm{p}}+\nabla {\overrightarrow{u}}_{\mathrm{p}}{}^T\right)+{\alpha }_{\mathrm{p}}\left({\lambda }_{\mathrm{p}}-\frac{2}{3}{\mu }_{\mathrm{p}}\right)\left(\nabla \cdot \nabla {\overrightarrow{u}}_{\mathrm{p}}\right)\overline{I}$ Granular shear viscosity using the Gidaspow model [73]: ${\mu }_{\mathrm{p}}={\mu }_{\mathrm{p},\mathrm{kin}}+{\mu }_{\mathrm{p},\mathrm{col}}+{\mu }_{\mathrm{p},\mathrm{fr}}$ Granular bulk viscosity using the Lun et al. model [74]: ${\lambda }_{\mathrm{p}}=\frac{4}{3}{\alpha }_{\mathrm{p}}{d}_{\mathrm{p}}{\rho }_{\mathrm{p}}{g}_0\left(1+e_{\mathrm{p},\mathrm{p}}\right){\left(\frac{{\Theta }_{\mathrm{p}}}{\pi }\right)}^{0.5}$ Granular pressure using the Syamlal–O’Brien model [75]: $p_{\mathrm{p}}={\alpha }_{\mathrm{p}}{\rho }_{\mathrm{p}}{\Theta }_{\mathrm{p}}+2{\rho }_{\mathrm{p}}\left(1-e_{\mathrm{p},\mathrm{p}}\right){\alpha }_{\mathrm{p}}^2{g}_0{\Theta }_{\mathrm{p}}$ Radial distribution function using the Lun et al. model [74]: $g_0={\left[1-{\left({\epsilon }_{\mathrm{p}}/{\epsilon }_{\mathrm{p},\max }\right)}^{1/3}\right]}^{-1}$ Collisional dissipation function using the Lun et al. model [74]: ${\gamma }_{\mathrm{p}}=\frac{12\left(1-e_{\mathrm{p}}^2\right){\epsilon }_{\mathrm{p}}^2{\rho }_{\mathrm{p}}{g}_0{\mathsf{\Theta }}_{\mathrm{p}}^{3/2}}{d_{\mathrm{p}}\sqrt{\pi }}$ Shear force of granular phase at wall: ${\overrightarrow{\tau }}_{\mathrm{p},\mathrm{wall}}=-\frac{\sqrt{3}\pi }{6}\phi \frac{{\alpha }_{\mathrm{p}}}{{\alpha }_{\mathrm{p},\max }}{\rho }_{\mathrm{p}}{g}_0\sqrt{{\Theta }_{\mathrm{p}}}{\overrightarrow{u}}_{\mathrm{p},\mathrm{wall}}$ Turbulent kinetic energy k and turbulent dissipation rate ε equation using the RNG k-ε model: $\frac{\partial }{\partial t}\left(\rho k\right)+\frac{\partial }{\partial x_i}\left(\rho k{u}_i\right)=\frac{\partial }{\partial x_i}\left[\left(\mu +\frac{{\mu }_t}{{\sigma }_k}\right)\frac{\partial k}{\partial x_j}\right]+G_k+G_b-\rho \epsilon -Y_M+S_k$ $\frac{\partial }{\partial t}\left(\rho \epsilon \right)+\frac{\partial }{\partial x_i}\left(\rho \epsilon u_i\right)=\frac{\partial }{\partial x_j}\left[\left(\mu +\frac{{\mu }_t}{{\sigma }_{\epsilon }}\right)\frac{\partial \epsilon }{\partial x_j}\right]+\rho C_{1}S\epsilon -\rho C_2\frac{{\epsilon }^2}{K+\sqrt{v\epsilon }}+C_{1\epsilon }\frac{\epsilon }{k}{C}_{3\epsilon }{G}_b+S_{\epsilon }$ Fluent generic erosion model [76]: $E_R={\displaystyle \sum _{p=1}^{N_p}\frac{m_{p}C\left(d_p\right)f\left(\theta \right)V^{n\left(u\right)}}{A_f}}$ where $C\left(d_p\right)=6.8\times {10}^{-9}$, $n\left(u\right)=2.41$, $f\left(\theta \right)=9.37-42.295{\theta }^2+110.864{\theta }^3-175.804{\theta }^4+170.137{\theta }^5-98.398{\theta }^6+31.211{\theta }^7-4.170{\theta }^8$. Discrete phase reflection coefficients of Grant and Tabakoff [63]: $e_{Nor}=0.993-1.76\theta +1.56{\theta }^2-0.49{\theta }^3$ $e_{\mathrm{Tan}}=0.998-1.66\theta +2.11{\theta }^2-0.67{\theta }^3$

. For better numerical convergence, the kinetic theory of granular flow [39] in its algebraic form was used to close the discrete phase stresses and ignore the coalescence process. It is possible to include the turbulence effect in the hydrogen gas phase viscosity. However, to constrain the number of adjustable parameters, we selected the RNG k-ε turbulence model for the gas phase as practiced by many other researchers for multiphase erosion flows. Furthermore, many previous studies [7,16,30,33,72] have proved that using the RNG k-ε turbulence model for the multiphase erosion in complex structures also has a reasonable agreement with experimental data.

As heterogeneous flow structures were found to significantly affect the non-uniform erosion process, the EMMS drag model [54] was adopted for modeling the gas–liquid droplet and gas–solid particle drag forces. The EMMS drag model reformulates the momentum equations of the original EMMS model by introducing inertial terms, thereby correlating the heterogeneity index of the drag coefficient, H_d, with both slip velocity and voidage. Simulation tests have shown that the EMMS drag is insensitive to the grid size and hence allows better applicability in the simulation of large-scale industrial units [43,46,57,60]. The heterogeneity index, H_d, for correction of the drag is provided in Appendix B. The generic erosion model of Ansys Fluent^® in Table 1 was adopted to predict the multiphase erosion process, in which the f(θ) correlation of erosion angle for SA106 Gr. B steel is derived from the experimental measurement in a jet testing machine.

2.3. Simulation Settings

Table 2. Summary of relevant parameter settings in simulation.

Types		Parameters	Values
inlet	velocity inlet	gauge pressure, MPa	2.68
		temperature, K	397.05
		turbulent intensity, %	5
		turbulent hydraulic diameter, m	0.3074
	hydrogen	velocity, uinlet,g m/s	47.71
		density, kg/m3	1.62
		volume flow rate, m3/s	3.54
		volume fraction, εhydrogen	99.26
		viscosity, kg/(m·s)	1.0862 × 10−5
	silicon tetrachloride droplet	velocity, uinlet,l m/s	0.36
		density, kg/m3	1248.31
		diameter, ddroplet m	0.0001
		viscosity, kg/(m·s)	3.40 × 10−4
		mass flux, kg/s	33.06
		volume fraction, εdroplet	0.74
	silicon particle	velocity, uinlet,s m/s	0.36
		density, kg/m3	2400
		diameter, ds m	0.00007
		mass flux, kg/s	0, 1
		volume fraction, εs	0, 0.01
outlet	outflow	discrete phase boundary condition type	escape
wall	material	density, kg/m3	7850
	boundary condition	hydrogen shear condition	no slip
		discrete phase shear condition	specularity coefficient, φ, 0.5
		wall roughness	standard

lists the industrial operation parameters in the simulation. In all cases, the hydrogen velocity is derived from calculating the measured gas volume flow rate divided by the cross-sectional area at the inlet. The gas and the discrete phases were assumed to flow uniformly into the top inlet and leave from the bottom outlet, where the atmospheric pressure boundary was prescribed. The initial velocities inside the pipe were assumed to be zero. The no-slip boundary condition was prescribed for the gas phase, whereas the partial-slip boundary condition developed by Johnson and Jackson [64] was used for the droplets and particles. The values for the specularity coefficient were chosen in line with the relevant literature [43,68,71] to analyze the influences of wall shear and lateral momentum transfer on erosion behaviors. The rebound behaviors between the discrete phase and wall were described by Grant and Tabakoff’s model [63] in accordance with the relevant literatures [7,16,30,33,72]. For the solution procedure, the SIMPLE algorithm for pressure–velocity coupling was used to obtain the solution for the continuous air phase. A second-order scheme was used for spatial discretization of the momentum and turbulence equations. PRESTO! and QUICK were used for spatial discretization of the pressure and volume fraction equations, respectively.

Furthermore, the calculation unit of erosion rate in Ansys Fluent^® is kg/(m²·s). The erosion thinning rates [77] were calculated with the expression below to directly compare the industrial measurement date.

$$R=\frac{E_R\times 8.76\times 3.6\times {10}^6}{{\rho }_w}$$

(1)

where E_R is the erosion rate with the Fluent generic erosion model, kg/(m²·s); R is erosion thinning rate, mm/a, whose unit means loss surface thickness of material per year; ρ_w is the pipe material density, g/cm³.

The erosion–dynamic mesh solver was used to further simulate the pipe erosion deformation process in Section 3.5. Therein, the variability time step was selected for the time-stepping method. The maximum node movement was set to 20% of the wall face length scale. The calculation time step was changed with iterations per flow simulation. The minimum and maximum time step sizes were set to 0.001 s and 10,000 s, respectively.

3. Results and Discussion

3.1. Model Validation

The designed thickness of the industrial pipe is 8.8 mm. Severe erosion damage often occurs in the local zones of the A and B elbows of this pipe after running for two years under industrial operation conditions, as shown in Figure 2, where the averaged thinning rates of SA106 Gr. B steel pipe walls are about 3 mm to 4 mm per year (i.e., 3–4 mm/a). Hydrogen gas and silicon tetrachloride droplets were assumed to flow uniformly into the pipeline at a high-speed gas velocity of 47.7 m/s and a lower liquid velocity of 0.3 m/s, respectively. Under the interphase drag force, the silicon tetrachloride droplets were quickly accelerated, and then heterogeneous flow following the hydrogen gas, appearing as the random clustering of droplets shown in Figure 3a. The clusters of silicon tetrachloride droplets are mesoscale structures expressing irregular strip line shapes, which influences the local gas–liquid interphase drag force. Additionally, a severe backflow appeared in the A elbow, accelerating the erosion, as shown in Figure 3b. Due to the turbulent flow, some backflow phenomena and inhomogeneous distribution of droplet volume fraction also occur in the B elbow. The turbulent flow along the pipeline flows to the bottom outlet.

Here, the hydrogen and liquid droplet erosion behaviors under the actual industrial running conditions were simulated with gas–liquid droplet EMMS drag models and a specularity coefficient φ = 0.5 for liquid droplet–wall partial slip. The predicted errors of erosion positions, such as bending angles and XYZ coordinates, are small compared to the industrial data in

Table 3. Comparison of predicted bending angle and XYZ coordinates of industrial erosion positions with measured in A and B elbows.

Type	Data	Bending Angle	Contact Ratio between Predicted and Industrial Data	Distance to Pipe Boundary along X	Distance to Pipe Boundary along Y	Distance to Pipe Boundary along Z
A elbow	Industrial data	26°~48°	/	316	472	191
	Predicted	25°~47°	0.917	340	485	175
B elbow	Industrial data	34°~68°	/	260	421	278
	Predicted	38°~69°	0.886	234	458	287

. Therein, a novel range contact ratio of C calculated with the below expression is introduced to quantitatively describe the deviation of the predicted bending angle range in comparison to industrialized measured data.

$$C=\frac{L_{B_P\cap B_i}}{L_{B_P\cup B_i}}$$

(2)

where, L is the length of the range; B_p is the range of predicted bending angle; B_i is the range of industrial measured bending angle. Here, ∩ and U symbols express the intersection and the union of the predicted and measured ranges, respectively. C closer to one means a smaller deviation.

At the A elbow, the predicted severe thinning rate at the back zone is 4.02 mm/a, close to the actual 3.2 mm/a calculated from industrial thickness measured data with an ultrasonic thickness gauge; the thickness value of every local measured point was averaged through arithmetic mean. Additionally, the whole tendencies of thinning rates at the erosion positions of B elbow are predicted and close to the industrial date, as shown in Figure 4. Differences exist for the predicted thinning rates at those points with industrial measured data because of a calculation error caused possibly by the current predicted boundary condition, some simplified models for discrete phases such as DDPM and KTGF, and the difference of predicted running conditions from the industrial parameters’ randomly fluctuation.

Above all, the erosion positions and rules are well predicted, except that the local wall thinning rate is overpredicted.

3.2. Effects of Key Simulation Parameters

To discuss single variable effects of the drag force and the liquid droplet specularity coefficient φ on the predicted erosion positions, different variables were compared with the simulation of hydrogen gas and liquid droplets of silicon tetrachloride, not considering the silicon particles due to very small solids volume fraction (0.01).

3.2.1. Drag Force between Hydrogen and Discrete Phase

The multiphase erosion simulations in the previous literature [16,19,78,79] often used the homogeneous drag models in the coarse mesh. To consider the effect of drag model on uneven multiphase flow, the multiscale flow erosion behaviors were contrasted respectively with homogeneous drag models and mesoscale drag models as shown in Figure 5. The worse prediction erosion distributions were found when using the Gidaspow or the spherical drag models. The simulated erosion positions and thinning rates with the EMMS drag model between hydrogen and discrete phase agree well with the industrial data in

Table 4. Comparison of predicted bending angle and XYZ coordinates of erosion positions industrially measured at different drag models.

Type	Data	Bending Angle	Contact Ratio between Predicted and Industrial Data	Distance to Pipe Boundary along X	Distance to Pipe Boundary along Y	Distance to Pipe Boundary along Z
A elbow	Industrial data	26°~48°	/	316	472	191
	Predicted with Gidaspow	13°~80°	0.338	584	131	157
	Predicted with spherical	34°~55°	0.500	381	369	274
	Predicted with EMMS	25°~47°	0.917	340	485	175
B elbow	Industrial data	34°~68°	/	260	421	278
	Predicted with Gidaspow	36°~60°	0.714	418	420	114
	Predicted with spherical	15°~73°	0.593	158	578	129
	Predicted with EMMS	38°~69°	0.886	234	458	287

. This is because the Gidaspow and spherical drag models are assumed to be homogeneous within a computational coarse-mesh multiphase simulation, ignoring the sub-grid structures of the particular droplet or particle cluster. This leads to an excessively uniform flow field and overestimated drag force [43,46,47,48]. However, the EMMS drag model is a heterogeneous approach derived from mesoscale-structure-based methods, and the mesoscale structures are broken down into a cluster phase and a dilute phase. It is applied with appropriate closure equations to address sub-grid structures and cluster formation [38,53,55], and the well describes the dynamic formation and dissolution of heterogeneous structures and significantly improves the accuracy of multiphase flow distribution.

Table 5. Comparison of predicted key flow parameters industrially measured at different drag models.

Data	Whole Pressure Drops, kPa	Hydrogen Velocity at Outlet, m/s	Mass Flux at Outlet, kg/s
Industrial data	80.00~100.00	/	0.20~0.30
Predicted with Gidaspow	42.69	1.79	0.21
Predicted with Spherical	76.55	3.61	0.25
Predicted with EMMS	98.09	5.55	0.28

further shows the comparison of predicted key flow parameters with industrial data at different drag models. It was found that lower hydrogen velocity, the mass flux at the outlet and total pressure drops (including the static and dynamic pressure drop) when using the Gidaspow and spherical uniform drag models. That is because over-predicted drag forces lead to overestimated local flow resistances in downward gas–liquid droplet flow. The flow resistance with the EMMS drag model is smaller; the mass flux at the outlet and pressure drops rises with the increase of hydrogen velocity as well as the industrial data.

Overall, the drag force between the hydrogen gas and the discrete phase is a key factor in the erosion behavior simulation.

3.2.2. Liquid Droplet Specularity Coefficient at wall Boundary Condition

Figure 6 shows the comparison of different liquid droplet specularity coefficients from 0.1 to 1 (i.e., no slip) on erosion behaviors. With respect to a perfectly specular collision at the wall results in larger solids transport, whereas the no-slip condition leads to stronger residence of particles near the wall and thus more pronounced, such a dependency is similar to that reported in the literature [43,62,68,69,70] that the wall specularity coefficient can significantly affect the predicted distribution.

With the decrease in specularity coefficient, it was found that the severe erosion position of A elbow deviates from the right side of the back bend to the left side and that the intensity of erosion thinning rate declines, but then the erosion range enlarges. That is because that droplet specularity coefficient changes the wall shear force and flow residence and then influences the velocity profiles of liquid droplets and cross-sectional slip velocity between droplets and hydrogen gas before entering into A elbow, as shown in Figure 6. This leads to different Reynolds numbers, causing the discrepancies of backflow range and the erosion positions of A elbow.

For the B elbow, the erosion positions are farther upward from the bottom due to the reduction in droplet specularity coefficient, as shown in Figure 6. The reason is that different wall flow residences affect the droplet accelerating process and its profiles of liquid droplet velocity, causing the increase in cross-sectional gas–liquid droplet slip velocity before entering into B elbow with the decrease of specularity coefficient, as shown in Figure 7.

Therein, the predicted droplet–wall partial slip condition of 0.5 is in better agreement with the industrial dates in terms of the bending angle and the coordinate values of X, Y and Z as compared in

Table 6. Comparison of predicted bending angle and XYZ coordinates of erosion positions with industrial measured at different specularity coefficients.

Type	Data	Bending Angle	Contact ratio between Predicted and Industrial Data	Distance to Pipe Boundary along X	Distance to Pipe Boundary along Y	Distance to Pipe Boundary along Z
A elbow	Industrial data	26°~48°	/	316	472	191
	Predicted with φ = 0.1	33°~53°	0.571	451	364	222
	Predicted with φ = 0.3	35°~51°	0.538	377	405	256
	Predicted with φ = 0.5	25°~47°	0.917	340	485	175
	Predicted with φ = 0.7	29°~45°	0.739	336	498	122
	Predicted with no slip (φ = 1.0)	32°~51°	0.654	387	432	81
B elbow	Industrial data	34°~68°	/	260	421	278
	Predicted with φ = 0.1	20°~42°	0.184	481	276	257
	Predicted with φ = 0.3	7°~46°	0.210	306	264	262
	Predicted with φ = 0.5	38°~69°	0.886	234	458	287
	Predicted with φ = 0.7	49°~87°	0.370	168	381	291
	Predicted with no slip (φ = 1.0)	53°~76°	0.372	295	504	216

. Thus, the liquid droplet specularity coefficient at the wall boundary condition is also a critical factor in the erosion behavior simulation.

3.3. Effects of Key Operation Parameters

3.3.1. Atomized Droplet Size

To analyze the effect of inlet-atomized droplet size derived from the mixer between hydrogen and silicon tetrachloride on the erosion behavior, 30 μm to 150 μm diameters according to the industrial measurement data were used to investigate, as shown in Figure 8. Therein, the 100 μm average atomized droplet size with Rosin–Rammler diameter is distributed from 30 μm to 150 μm, and its expression is given by [80]:

$$Y_d=e^{{-\left(d/\overline{d}\right)}^n}$$

(3)

where $\overline{d}$ is the size constant; n is the size distribution parameter and here is 3.5.

With the increase in atomized droplet size, the erosion range and thinning rate all gradually enlarge. Therein, smaller atomized droplets with such sizes of 30 mm, cause erosion less due to the droplets following with the hydrogen. The erosion ranges and intensities become more dispersible and weaker, respectively, in the A and B elbows when using the 100 mm average atomized droplet size. That is because the carrying capacities of hydrogen gas on different sizes of liquid droplets have large discrepancies. The gas flow easily entrains smaller droplets, while a bigger droplet means greater inertia and leads to difficulty in following the same gas flow and changing flow direction in time.

As for A elbow, the erosion maximum thinning rate enlarges, and the growth rate gradually trends to being smooth, but the bending angles at severe erosion positions are all located at about 30°–50°, with liquid droplet size increasing as shown in Figure 9. That is because the impact momentum of a single droplet has reached the peak by only increasing its mass under the carrying capacities of hydrogen gas unchanged. Despite further increasing atomized droplet size can enlarge single droplet mass, its velocity will decline. Additionally, more high-momentum droplets will impact the wall, causing an increase in erosion range. At the B elbow, the maximum erosion position and bending angle rates all rise under the action of gravity acceleration with the increase in droplet size shown in Figure 10. That indicates that the atomized droplet size easily affected the erosion behaviors, such as the erosion position, thinning rate, and bending angle.

3.3.2. Hydrogen Volume Fraction

To consider the effect of process fluctuation in the mixer on the hydrogen gas volume fraction of the inlet, the range of 93% to 99.2% according to the possible industrial conditions is selected to be discussed, as shown in Figure 11.

With the hydrogen gas volume fraction increasing, the flow velocity also enlarges from 4.74 m/s to 47.7 m/s, so the maximum thinning rates rise. Moreover, the bending angles at severe erosion positions also have few differences because of the small flow change.

At the A elbow, the erosion maximum thinning rate and bending angle exist few differences in the range of hydrogen volume fraction from 93% to 98.5% in Figure 12a, but obvious changes appear at 99.2% because the hydrogen velocity sharply increases to 47.7 m/s. For the B elbow, the thinning rate is approximately linear rise under the action of gravity acceleration as shown in Figure 12b, and the bending angle obviously enlarges at the 99.2% hydrogen volume fraction.

3.4. Coupling Effects of Liquid Droplet and Solid Particle Phases on Erosion

To consider further the effect of solid particles on industrial erosion behavior, a few particles in terms of the measured particle volume fraction (0.01%) are injected into the pipe, and the variation in erosion behavior is observed in Figure 13. Here, the solid particle specularity coefficient was set to 0.5, and the drag model between hydrogen gas and solid particles was also selected as the EMMS drag model. The other simulation settings are listed in Table 2.

The erosion area increases, but the position is unchanged at A elbow due to the existence of a larger vortex at the gas–solid–liquid droplets phase as shown in Figure 14a. As for the B elbow, the erosion position transfers upward to a positive Z direction. This is because more discrete phases (liquid droplet and solid particle) improve the flow resistance and then change the flow pattern as shown in Figure 14b.

3.5. Erosion Deformation Process in the Gas and Liquid Droplet Two-Phase Flow

To investigate the rule of erosion deformation and failure processes of this industrial pipe in the hydrogen gas and silicon tetrachloride liquid droplet two-phase flow, Figure 15 shows the extent of A and B elbows’ erosion deformation at different times with an erosion-dynamic mesh solver under the conditions of no solid particles and using the EMMS drag model and a 0.5 liquid droplet specularity coefficient. Other simulation parameters are shown in Table 2. To accurately solve erosion deformation, the calculation with erosion-coupled dynamic mesh starts after the steady erosion behavior using the static mesh approach. The residual curves converge in every time step.

In the beginning, the erosion ranges and thinning rates gradually enlarge, and the erosion positions remain unchanged with time; meanwhile, the wall surfaces appear to a certain extent of local wall-thinning deformation. After several hours, such as 4500 s, the erosion range reaches a stable state of range change at a slow velocity, but the thinning rates at the severe erosion positions further enlarge. More obvious thinning deformations appear at the positions of central local erosion points, leading to a local stagnated zone and easy-to-accelerate granular particles gathering and erosion toward the depth; meanwhile, the local range of erosion points also gradually expands to the surroundings with time.

4. Conclusions

The CFD-DDPM approach, coupled with the sub-grid EMMS drag model, was used to analyze the complex multiphase erosion phenomenon in a connecting pipe of an industrial polycrystalline silicon unit. The simulated erosion behaviors are very consistent with the industrial data. The main conclusions are as follows:

(1)

The drag force between the hydrogen and discrete phase and the droplet wall specularity coefficient are key effect factors on the predicted erosion position and thinning rate. The non-uniform multiphase erosion flow behavior near the wall in a coarse mesh can be modeled accurately with the sub-grid EMMS drag model. The droplet specularity coefficient influences the wall momentum exchange and flow distribution. A suitable partial slip boundary condition, such as a 0.5 specularity coefficient, could improve the accuracy of erosion position.

(2)

The variations of operation parameters directly influence the erosion positions and thinning rates. Small liquid droplets, such as those of 30 μm size, will follow the gas phase better and have a lower erosion rate. The inertia effect of large droplets, such as those of 150 μm size, plays a dominant role, resulting in obvious erosion on the elbow walls. The erosion range and thinning rate enlarge with the increase in hydrogen volume fraction. A greater hydrogen volume fraction will have a greater impact on kinetic energy and produce more serious erosion damage.

(3)

A few silicon solid particles (such as 0.01% volume fraction) can change local flow behaviors, such as flow velocity and residence, and probably cause the variation of local erosion positions.

(4)

In the gas–liquid droplet erosion thinning deformation process, the erosion ranges and thinning rates enlarge quickly. The erosion positions remain unchanged at the beginning, whereafter the pipe deformation accelerates to improve the thinning rates at the severe erosion positions, but the erosion range changes little.