# Hydrodynamics of an Airlift Column for Aeration in Molten Sulfur
## 3 CFD Modeling  
The commercial CFD software package ANSYS FLUENT19.0 was used to model the hydrodynamics of airlift column in this work. Considering the flow characteristics, the \(k-\varepsilon\) turbulence model was adopted in the simulation due to its stability and convergence [15]. It is a two-equation model that includes the transport equations for the turbulent kinetic energy and dissipation rate [16]. The model is suitable for complex secondary flows [17]. The two-fluid Euler-Euler model [18; 19] was used to solve the governing equations for each phase under the condition of adiabatic process. The interfacial momentum forces in the momentum equation were considered to simulate the interactions between the continuous and dispersed phases. The drag, lift, and turbulent dispersion forces are important for gas-liquid flows in the present simulation. The Ishii and [PERSON] drag coefficient model [20] that is valid for the modelling of individual bubble motion under a wide range of Reynolds number was employed in this paper. The population balance model was used to describe the coalescence and breakup of bubbles [21]. The population balance method discretizes the transport equation into different groups based on the predefined classification of bubble size distribution. The equal-diameter model was used for the size group distribution. The bubble diameter ranged from 0.5 to 25 mm. The CFD model and calculation formulas used in this work: Equations (S1)-(S24), and the values used in k-\(\varepsilon\) model: Table S1 were described in the Supplementary Materials.  
The \"velocity inlet\" boundary condition was used at the gas sparger and the \"pressure outlet\" condition was used at the outlet. The specification of inlet velocity was based on the superficial gas velocity ranging from 0.0056 to 0.05 m/s. Considering the precision and computational costs, the entire computational domain was divided into around 2 million unstructured tetrahedral cells, as shown in Figure 2. Mesh convergence tests were conducted to ensure the independence of the solution on the mesh size. The number of grids was 1, 1.6, and 2.2 million, respectively. The average gas holdup in the rising area and downward area of different height sections of the column was investigated, and the column height was selected as 15 m. As shown in Figure 3, compared with the number of 2.2 million grids, the maximum error of simulation results of 1 and 1.6 million grids is no more than 2%, which can be considered grid independent. In order to ensure the accuracy of simulation, 2.2 million grids were selected. The time step was set as 2 \(\times\) 10\({}^{-2}\) s and the tests showed that the results were not sensitive to the time step. The column was filled with the liquid during the initialization. The liquid density and gas density were 1090 and 20 kg/m\({}^{3}\), respectively. The liquid viscosity was 0.18 Pa-s at the temperature of 160 \({}^{\circ}\)C.  
To investigate the local hydrodynamic parameters, the riser was divided into six equal parts with the height of \(h\), starting from the top of \(H_{2}\), as shown in Figure 1. Therefore, seven cross sections along the riser were obtained. The bottom cross section was the first cross section.  
The contours of the gas holdup, liquid axial velocity, and bubble diameter were used to illustrate the hydrodynamics. The influence of the total height of the column, \(H\), and the superficial gas velocity, \(U_{Gr1}=Q_{G,in}/A_{r}\), due to the gas injection on the hydrodynamics was discussed.