Hydrodynamics of an Airlift Column for Aeration in Molten Sulfur

The airlift column is a promising technology for the removal of volatile gas from high-viscosity molten sulfur. However, a detailed analysis is lacking on the hydrodynamic properties inside the column, due to the difficulty in flow behavior detection in the opaque molten sulfur. In this work, we adopted the computational fluid dynamics simulation to understand the hydrodynamic behaviors in an airlift column for molten sulfur aeration. In addition, we analyzed the impacts of the superficial gas velocity (UGr) and column height on the hydrodynamic characteristics, such as gas holdup, average bubble diameter, and liquid circulation velocity (ULr) in the column. The simulation shows that at a constant column height of 15 m, an increase on gas holdup can be obtained with the increase of the superficial gas velocity, while the bubble diameter remains almost constant. Once the superficial gas velocity exceeded 0.333 m/s, the liquid circulation velocity increased slowly. With a variation on the column height from 5 to 25 m, a negligible change on gas holdup, but an obvious increase on liquid circulation velocity and bubble diameter is observed at the given superficial gas velocity of 0.0389 m/s. Furthermore, the simulation shows a similar trend, but with considerably more detailed information, on the relationship between the gas holdup and liquid circulation velocity when compared to the predictions from the Chisti correlation (1988) and an optimized correlation proposed in this work.


Introduction
The Claus process is widely used in the treatment of hydrogen sulfide found in raw natural gas and by-product gases derived from oil refineries [1]. However, one of the main challenges in elemental sulfur production is the dissolution of hydrogen sulfide into molten sulfur, which can escape into the surrounding environment gradually in the later stage [2]. Currently, aeration or agitation is the most commonly used method for the removal of hydrogen sulfide from molten sulfur. Air stripping and injected oxygen components contribute to the separation of dissolved hydrogen sulfide from molten sulfur. A suitable flow field of molten sulfur facilitates the gas-liquid mass transfer. However, agitation involves rotating equipment with a complex structure, which are not suitable for a long-term operation. Therefore, a passive mixture column has been proposed.
An aeration airlift column exhibits liquid circulation and a gas upflow pattern, which provide good mass transfer and mixing performances. The airlift column consists of an upward zone and a downcomer. Pressurized air in the form of dispersed bubbles enters the upward zone through the gas distribution device located at the lower part of the airlift zone. Since the gas holdup in the riser is higher than in the downcomer, there is a gas holdup difference between the riser and the downcomer, forming a liquid circulation driving force. Therefore, the liquid in the downcomer exhibits a continuous downward trend. The liquid flows upward in the upward zone and downward in the downcomer to form a liquid circulation pattern in the entire column.
The aeration airlift column is suitable for the removal of dissolved hydrogen sulfide from molten sulfur. Compared with the prior art, the airlift column significantly changes the fluid flow state in the column and enhances the mixing effect due to the existence of the draft tube. The hydrodynamics of the airlift column affect the oxidation and removal of hydrogen sulfide. However, there are few studies on the hydrodynamics of aeration airlift columns with the industrial-scale physical properties of molten sulfur. This work focuses on the effects of superficial gas velocity and height on the internal flow field of the airlift column. Three main parameters of airlift column were investigated: Gas holdup, liquid circulation velocity, and bubble diameter. Some empirical correlations are optimized by considering the influence of bubble diameter, and through a comparison with the Chisti correlation and simulation results. The schematic of the airlift column is shown in Figure 1.
holdup difference between the riser and the downcomer, forming a liquid circulation driving force. Therefore, the liquid in the downcomer exhibits a continuous downward trend. The liquid flows upward in the upward zone and downward in the downcomer to form a liquid circulation pattern in the entire column.
The aeration airlift column is suitable for the removal of dissolved hydrogen sulfide from molten sulfur. Compared with the prior art, the airlift column significantly changes the fluid flow state in the column and enhances the mixing effect due to the existence of the draft tube. The hydrodynamics of the airlift column affect the oxidation and removal of hydrogen sulfide. However, there are few studies on the hydrodynamics of aeration airlift columns with the industrial-scale physical properties of molten sulfur. This work focuses on the effects of superficial gas velocity and height on the internal flow field of the airlift column. Three main parameters of airlift column were investigated: Gas holdup, liquid circulation velocity, and bubble diameter. Some empirical correlations are optimized by considering the influence of bubble diameter, and through a comparison with the Chisti correlation and simulation results. The schematic of the airlift column is shown in Figure 1.
The outer column of the airlift column had an inner diameter of 1600 mm, and different heights of the column were considered. The inner diameter of the concentric draft tube was 1150 mm. Elliptical heads with a short axis length of 400 mm were placed at both ends of the outer column. Nine gas injection pipes were inserted into the riser. The detailed parameters of the airlift column are shown in Table 1.   D1  D2  D3  D4  D5  D6  H1  H2  H3  H4  H5  mm  mm  mm  mm  mm  mm  mm  mm  mm  mm  mm  1600  1150  200  50  690  230  600  200  150 400 100 The outer column of the airlift column had an inner diameter of 1600 mm, and different heights of the column were considered. The inner diameter of the concentric draft tube was 1150 mm. Elliptical heads with a short axis length of 400 mm were placed at both ends of the outer column. Nine gas injection pipes were inserted into the riser. The detailed parameters of the airlift column are shown in Table 1.

Mathematical Modeling
Mathematical modeling provides correlations for hydrodynamic parameters. Using the energy balance principle, a predictive model for the liquid circulation velocity was Appl. Sci. 2022, 12, 117 3 of 15 developed by Chisti and Moo-Yong [3]. The model provides a simplified correlation to simultaneously predict the superficial liquid velocity and void fraction in a riser using an iterative procedure [4]. The results obtained by this correlation show a better agreement with the measurements, as compared with the results obtained using the correlations proposed by Bello and Choi et al. [5].
Only the physically relevant root of the polynomial was used. A steady-state value of U Lr was obtained. For the column discussed in this work, the value of K B was 11.7.
The Chisti correlation proposed a correlation that considered the circulation of bubbles in the downcomer and it does not consider the bubble diameter. Generally, a narrow bubble diameter distribution is expected at a high temperature since the liquid viscosity and surface tension decrease with the increase of temperature [6]. However, the bubble size distribution is considerably wider in the molten sulfur condition with a high viscosity, giving the requirement on the consideration of bubble diameter in an optimized correlation.
Equations (5) and (6) are the two well-known correlations for the bubble diameter proposed by Akita and Wilkinson et al., respectively. The correlation obtained by Wilkinson et al. accurately predicts the Sauter mean diameter since it considers gas density [7]. Therefore, this correlation is adopted to assess the bubble diameter [8].
Based on the two-phase drift flux model deduced by Zuber and Findlay [9], the gas velocity in the downcomer can be expressed as follows: where C d is the distribution parameter. According to [10], the value of C d is 1.065. According to Gao et al., the gas holdup in an airlift column is expressed as follows [11]: where U Gr includes U Gr1 and U Gr2 [12]: One for gas injection, i.e., U Gr1 = Q G,in /A r , and the other for the downcomer, U Gr2 = U Gr − U Gr1 . Therefore, the downcomer gas velocity can be calculated as: Jamialahmadi proposed new correlations to predict the terminal velocity of bubbles over a wide range of gas and liquid properties [13]. These correlations are suitable for high temperatures [14].

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 − ε 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 Zuber 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-ε 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 × 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 • 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.  To investigate the local hydrodynamic parameters, the riser was divided into six equal parts with the height of ℎ, starting from the top of , 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, , and the superficial gas velocity, = , / , due to the gas injection on the hydrodynamics was discussed.

Gas Holdup
The area-averaged value of the gas holdup at a cross section of the riser is defined as . The total gas holdup is the volume-averaged gas holdup in the riser zone. Figure 4 shows the contours of the gas holdup with at a total column height of 15 m. The gas holdup increases with the superficial gas velocity . Apparently, the gas holdup at the center of the riser is larger than the sidewalls of the riser. The growth of gas volume increases the gas holdup difference between the riser and the downcomer, accelerates the liquid circulation, and brings more bubbles into the annulus. With the increase of , a small amount of bubbles gradually circulates from the downcomer into the riser.  To investigate the local hydrodynamic parameters, the riser was divided into six equal parts with the height of ℎ, starting from the top of , 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, , and the superficial gas velocity, = , / , due to the gas injection on the hydrodynamics was discussed.

Gas Holdup
The area-averaged value of the gas holdup at a cross section of the riser is defined as . The total gas holdup is the volume-averaged gas holdup in the riser zone. Figure 4 shows the contours of the gas holdup with at a total column height of 15 m. The gas holdup increases with the superficial gas velocity . Apparently, the gas holdup at the center of the riser is larger than the sidewalls of the riser. The growth of gas volume increases the gas holdup difference between the riser and the downcomer, accelerates the liquid circulation, and brings more bubbles into the annulus. With the increase of , a small amount of bubbles gradually circulates from the downcomer into the riser.

Gas Holdup
The area-averaged value of the gas holdup at a cross section of the riser is defined as ε a . The total gas holdup is the volume-averaged gas holdup in the riser zone. Figure 4 shows the contours of the gas holdup with U Gr at a total column height of 15 m. The gas holdup increases with the superficial gas velocity U Gr . Apparently, the gas holdup at the center of the riser is larger than the sidewalls of the riser. The growth of gas volume increases the gas holdup difference between the riser and the downcomer, accelerates the liquid circulation, and brings more bubbles into the annulus. With the increase of U Gr , a small amount of bubbles gradually circulates from the downcomer into the riser. Figure 5 shows the gas holdup distribution with the total column height H at a constant U Gr of 0.0389 m/s. The gas holdup exhibits similar distributions at different values of column heights. It can be seen that even if the column height is 5 m, a large number of bubbles cannot circulate through the annular gap into the riser, indicating that the column height is not a leading factor of bubble circulation. Figure 5 shows the gas holdup distribution with the total column height at a constant of 0.0389 m/s. The gas holdup exhibits similar distributions at different values of column heights. It can be seen that even if the column height is 5 m, a large number of bubbles cannot circulate through the annular gap into the riser, indicating that the column height is not a leading factor of bubble circulation.  The gas is enriched at the top region of the downcomer, and the bottom region of the downcomer contains a negligible amount of the gas. This is in agreement with flow regime II according to Heijnen et al. [22]. The bubble-enriched zone in the downcomer in the  Figure 5 shows the gas holdup distribution with the total column height at a constant of 0.0389 m/s. The gas holdup exhibits similar distributions at different values of column heights. It can be seen that even if the column height is 5 m, a large number of bubbles cannot circulate through the annular gap into the riser, indicating that the column height is not a leading factor of bubble circulation.  The gas is enriched at the top region of the downcomer, and the bottom region of the downcomer contains a negligible amount of the gas. This is in agreement with flow regime II according to Heijnen et al. [22]. The bubble-enriched zone in the downcomer in the The gas is enriched at the top region of the downcomer, and the bottom region of the downcomer contains a negligible amount of the gas. This is in agreement with flow regime II according to Heijnen et al. [22]. The bubble-enriched zone in the downcomer in the column with a height of 5 m is smaller than the column with a height of 25 m. The displacement of the bubbles in the downcomer tends to increase with the total column height. Figure 6 shows the variation of the area-averaged gas holdup at different cross sections (local heights) of the riser when the total column height is 15 m. It is interesting that there is a downturn region when the gas holdup increases with the increment of local height in a global trend. Moreover, the area-averaged gas holdup is almost constant in the riser region between the second and sixth cross sections. It shows that the bubbles were ejected from the bubble generator, after full development and uniformity, it belongs to the homogeneous flow, with only a small amount of bubble coalescence and breakup.
column with a height of 5 m is smaller than the column with a height of 25 m. The displacement of the bubbles in the downcomer tends to increase with the total column height. Figure 6 shows the variation of the area-averaged gas holdup at different cross sections (local heights) of the riser when the total column height is 15 m. It is interesting that there is a downturn region when the gas holdup increases with the increment of local height in a global trend. Moreover, the area-averaged gas holdup is almost constant in the riser region between the second and sixth cross sections. It shows that the bubbles were ejected from the bubble generator, after full development and uniformity, it belongs to the homogeneous flow, with only a small amount of bubble coalescence and breakup.  Figure 7 shows the area-averaged gas holdup at the cross sections of the riser at different total column heights with a fixed of 0.0389 m/s. The area-averaged gas holdup increases with the , demonstrating a similar trend for different total column heights. Figures 6 and 7 show that the gas holdup is the lowest at the first cross section, which indicates the undeveloped multiphase flow at this location.   Figures 8 and 9 show the comparison of the total gas holdup calculated by the Chisti correlation [23], the newly proposed correlation, and the CFD simulation results in the present study. Figure 8 implies that the total gas holdup calculated by the simulation is lower than the gas holdup calculated by the Chisti correlation, while the proposed  Figure 7 shows the area-averaged gas holdup at the cross sections of the riser at different total column heights with a fixed U Gr of 0.0389 m/s. The area-averaged gas holdup increases with the U Gr , demonstrating a similar trend for different total column heights. Figures 6 and 7 show that the gas holdup is the lowest at the first cross section, which indicates the undeveloped multiphase flow at this location.
column with a height of 5 m is smaller than the column with a height of 25 m. The displacement of the bubbles in the downcomer tends to increase with the total column height. Figure 6 shows the variation of the area-averaged gas holdup at different cross sections (local heights) of the riser when the total column height is 15 m. It is interesting that there is a downturn region when the gas holdup increases with the increment of local height in a global trend. Moreover, the area-averaged gas holdup is almost constant in the riser region between the second and sixth cross sections. It shows that the bubbles were ejected from the bubble generator, after full development and uniformity, it belongs to the homogeneous flow, with only a small amount of bubble coalescence and breakup.  Figure 7 shows the area-averaged gas holdup at the cross sections of the riser at different total column heights with a fixed of 0.0389 m/s. The area-averaged gas holdup increases with the , demonstrating a similar trend for different total column heights. Figures 6 and 7 show that the gas holdup is the lowest at the first cross section, which indicates the undeveloped multiphase flow at this location.   Figures 8 and 9 show the comparison of the total gas holdup calculated by the Chisti correlation [23], the newly proposed correlation, and the CFD simulation results in the present study. Figure 8 implies that the total gas holdup calculated by the simulation is lower than the gas holdup calculated by the Chisti correlation, while the proposed  Figures 8 and 9 show the comparison of the total gas holdup calculated by the Chisti correlation [23], the newly proposed correlation, and the CFD simulation results in the present study. Figure 8 implies that the total gas holdup calculated by the simulation is lower than the gas holdup calculated by the Chisti correlation, while the proposed correlation is higher. Figure 9 shows the results for different total column heights at a constant U Gr of 0.0389 m/s. The total gas holdup calculated by the simulation is slightly lower than the two correlations. The simulation results exhibit a good agreement. Although the proposed correlation considers the bubble diameter to extend the validity of the calculation results, the assumption that bubbles enter from the downcomer into the riser at the bottom of the column may cause some deviations, therefore, the value is larger. The proposed correlation may be more accurate for other gas-liquid conditions since the gas holdup is also affected by physical and operating parameters, such as liquid viscosity [24], liquid surface tension, and pressure. correlation is higher. Figure 9 shows the results for different total column heights at a constant of 0.0389 m/s. The total gas holdup calculated by the simulation is slightly lower than the two correlations. The simulation results exhibit a good agreement. Although the proposed correlation considers the bubble diameter to extend the validity of the calculation results, the assumption that bubbles enter from the downcomer into the riser at the bottom of the column may cause some deviations, therefore, the value is larger. The proposed correlation may be more accurate for other gas-liquid conditions since the gas holdup is also affected by physical and operating parameters, such as liquid viscosity [24], liquid surface tension, and pressure.

Height (m)
Correlation-Chisti Correlation of this work Simulation Figure 9. Comparison of total gas holdup calculated by the correlations and CFD simulation at different total column heights.

Liquid Circulation Velocity
The liquid circulation is caused by the gas holdup difference between the riser and downcomer. The liquid circulation velocity is a significant hydrodynamic parameter of the airlift column, which is associated with the liquid axial velocity and gas holdup in riser , = (1 − ). correlation is higher. Figure 9 shows the results for different total column heights at a constant of 0.0389 m/s. The total gas holdup calculated by the simulation is slightly lower than the two correlations. The simulation results exhibit a good agreement. Although the proposed correlation considers the bubble diameter to extend the validity of the calculation results, the assumption that bubbles enter from the downcomer into the riser at the bottom of the column may cause some deviations, therefore, the value is larger. The proposed correlation may be more accurate for other gas-liquid conditions since the gas holdup is also affected by physical and operating parameters, such as liquid viscosity [24], liquid surface tension, and pressure.

Height (m)
Correlation-Chisti Correlation of this work Simulation Figure 9. Comparison of total gas holdup calculated by the correlations and CFD simulation at different total column heights.

Liquid Circulation Velocity
The liquid circulation is caused by the gas holdup difference between the riser and downcomer. The liquid circulation velocity is a significant hydrodynamic parameter of the airlift column, which is associated with the liquid axial velocity and gas holdup in riser , = (1 − ).

Liquid Circulation Velocity
The liquid circulation is caused by the gas holdup difference between the riser and downcomer. The liquid circulation velocity U Lr is a significant hydrodynamic parameter of the airlift column, which is associated with the liquid axial velocity V Lr and gas holdup in riser ε r , U Lr = V Lr (1 − ε r ). Figure 10 shows the contours of the liquid axial velocity at different superficial gas velocities. The V Lr increases with the U Gr . The growth rate of the V Lr decreases when the U Gr is larger than 0.0333 m/s. Figure 11 shows the contours of the V Lr at different total column heights with U Gr of 0.0389 m/s. The V Lr increases with the total column height. This is due to the fact that as the total column height increases, the hydrostatic pressure increases, and consequently, the driving force generated by the differential pressure between the riser and the downcomer increases. Finally, the enhanced driving force is an increase in the liquid circulation velocity.
velocities. The increases with the . The growth rate of the decreases when the is larger than 0.0333 m/s. Figure 11 shows the contours of the at different total column heights with of 0.0389 m/s. The increases with the total column height. This is due to the fact that as the total column height increases, the hydrostatic pressure increases, and consequently, the driving force generated by the differential pressure between the riser and the downcomer increases. Finally, the enhanced driving force is an increase in the liquid circulation velocity.  When the superficial gas velocity increases, the liquid circulation velocity gradually increases when the gas holdup is high, as reported by previous studies [25][26][27][28]. This is velocities. The increases with the . The growth rate of the decreases when the is larger than 0.0333 m/s. Figure 11 shows the contours of the at different total column heights with of 0.0389 m/s. The increases with the total column height. This is due to the fact that as the total column height increases, the hydrostatic pressure increases, and consequently, the driving force generated by the differential pressure between the riser and the downcomer increases. Finally, the enhanced driving force is an increase in the liquid circulation velocity.  When the superficial gas velocity increases, the liquid circulation velocity gradually increases when the gas holdup is high, as reported by previous studies [25][26][27][28]. This is Figure 11. Axial liquid velocity distribution at different total column heights with U Gr of 0.0389 m/s. When the superficial gas velocity increases, the liquid circulation velocity gradually increases when the gas holdup is high, as reported by previous studies [25][26][27][28]. This is due to the fact that the displacement of the bubbles in the downcomer tends to increase with the superficial gas velocity. As seen from Figure 4, the gas holdup difference between the riser and downcomer decreases as the U Gr increases. As a result, the liquid axial velocity no longer changes evidently.
The comparison of the liquid circulation velocity calculated by the Chisti correlation, the proposed correlation, and the present simulation are shown in Figures 12 and 13. As shown in Figure 12, when the U Gr is lower than 0.0389 m/s, the simulation results are larger than the results calculated by the two correlations.
with the superficial gas velocity. As seen from Figure 4, the gas holdup difference between the riser and downcomer decreases as the increases. As a result, the liquid axial velocity no longer changes evidently.
The comparison of the liquid circulation velocity calculated by the Chisti correlation, the proposed correlation, and the present simulation are shown in Figures 12 and 13. As shown in Figure 12, when the is lower than 0.0389 m/s, the simulation results are larger than the results calculated by the two correlations.  Figure 13 illustrates the results for different total column heights with of 0.0389 m/s. The results calculated by the Chisti correlation fit in well with the simulation results, while the results calculated by the proposed correlation are slightly larger than the simulation results. The results calculated by the two correlations exhibit the same variation trends. When the total column height is less than 10 m and more than 15 m, the simulation results are larger and smaller than the results obtained using the two correlations, respectively.   Figure 13 illustrates the results for different total column heights with of 0.0389 m/s. The results calculated by the Chisti correlation fit in well with the simulation results, while the results calculated by the proposed correlation are slightly larger than the simulation results. The results calculated by the two correlations exhibit the same variation trends. When the total column height is less than 10 m and more than 15 m, the simulation results are larger and smaller than the results obtained using the two correlations, respectively.

Height (m)
Correlation-Chisti Correlation of this work Simulation Figure 13. Comparison of liquid circulation velocity calculated from the correlations and CFD simulation at different total column heights. Figure 13. Comparison of liquid circulation velocity calculated from the correlations and CFD simulation at different total column heights. Figure 13 illustrates the results for different total column heights with U Gr of 0.0389 m/s. The results calculated by the Chisti correlation fit in well with the simulation results, while the results calculated by the proposed correlation are slightly larger than the simulation results. The results calculated by the two correlations exhibit the same variation trends. When the total column height is less than 10 m and more than 15 m, the simulation results are larger and smaller than the results obtained using the two correlations, respectively.

Bubble Diameter
The equilibrium value of the average bubble diameter is determined by the local hydrodynamic conditions. The well-known correlation for the bubble diameter is as follows [29,30]: The bubble diameter distribution is an important parameter of hydrodynamics, which depends on the coalescence and breakup of bubbles [31]. Figure 14 shows the average bubble diameter distribution with superficial velocity at a total column height of 15 m. The bubble diameter gradually increases from the nozzle to the top of the column, and the large bubbles in the downcomer are mainly concentrated at the top. The bubble diameter at the middle and bottom of the downcomer is smaller than the riser. Figure 15 shows the average bubble diameter distribution in the riser for different total column heights with U Gr of 0.0389 m/s. As the total column height increases, the average bubble diameter gradually increases up to a height of 20 m and decreases when the height is 25 m. It shows that the column height has a significant effect on the average bubble diameter. The bubble diameter grows gradually from the bottom to the top. As the bubble reaches the maximum critical stable size, the bubble begins to breakup. Herein, the breakup of the bubbles will increase the turbulence in the surrounding flow field, which causes additional bubbles to break. dynamic conditions. The well-known correlation for the bubble diameter is as follows [29,30]: The bubble diameter distribution is an important parameter of hydrodynamics, which depends on the coalescence and breakup of bubbles [31]. Figure 14 shows the average bubble diameter distribution with superficial velocity at a total column height of 15 m. The bubble diameter gradually increases from the nozzle to the top of the column, and the large bubbles in the downcomer are mainly concentrated at the top. The bubble diameter at the middle and bottom of the downcomer is smaller than the riser. Figure 15 shows the average bubble diameter distribution in the riser for different total column heights with of 0.0389 m/s. As the total column height increases, the average bubble diameter gradually increases up to a height of 20 m and decreases when the height is 25 m. It shows that the column height has a significant effect on the average bubble diameter. The bubble diameter grows gradually from the bottom to the top. As the bubble reaches the maximum critical stable size, the bubble begins to breakup. Herein, the breakup of the bubbles will increase the turbulence in the surrounding flow field, which causes additional bubbles to break.   Figure  16 shows the average bubble diameter at different cross sections with at a total column height of 15 m. The average bubble diameter gradually increases from 10 to 22 mm  Figure 16 shows the average bubble diameter at different cross sections with U Gr at a total column height of 15 m. The average bubble diameter gradually increases from 10 to 22 mm from the first to the seventh cross section. However, the average bubble diameter decreases in the upper part of the column when the U Gr is more than 0.0333 m/s. As the U Gr increases, the turbulence degree of the column increases, which results in the decrease of the maximum stable bubble diameter and the average bubble diameter [32,33]. A few researchers argued that the average bubble diameter slightly increases with the increasing U Gr [34]. A maximum stable bubble diameter exists at a critical value of the U Gr . When the U Gr is less than the critical value, the average bubble diameter slightly increases with the U Gr . However, when the U Gr is larger than the critical value, a severe bubble breakup occurs and the average bubble diameter decreases.  Figure  16 shows the average bubble diameter at different cross sections with at a total column height of 15 m. The average bubble diameter gradually increases from 10 to 22 mm from the first to the seventh cross section. However, the average bubble diameter decreases in the upper part of the column when the is more than 0.0333 m/s. As the increases, the turbulence degree of the column increases, which results in the decrease of the maximum stable bubble diameter and the average bubble diameter [32,33]. A few researchers argued that the average bubble diameter slightly increases with the increasing [34]. A maximum stable bubble diameter exists at a critical value of the . When the is less than the critical value, the average bubble diameter slightly increases with the . However, when the is larger than the critical value, a severe bubble breakup occurs and the average bubble diameter decreases.  Figure 17 shows the average bubble diameter at different cross sections for various total column heights with the U Gr of 0.0389 m/s. The average bubble diameter increases with the local height up to a total column height of 20 m. In addition, the growth rate of the bubble diameter from the first to the third section increases with the total column height. At the total column heights of 20 and 25 m, the bubble diameter decreases from the fifth and the third section onwards, respectively. This is investigated on the basis of the turbulent dissipation in each section.  Figure 17 shows the average bubble diameter at different cross sections for various total column heights with the of 0.0389 m/s. The average bubble diameter increases with the local height up to a total column height of 20 m. In addition, the growth rate of the bubble diameter from the first to the third section increases with the total column height. At the total column heights of 20 and 25 m, the bubble diameter decreases from the fifth and the third section onwards, respectively. This is investigated on the basis of the turbulent dissipation in each section.  Figure 18 shows the turbulent eddy dissipation at different cross sections for various total column heights with the of 0.0389 m/s. When the gas holdup is constant, the bubble diameter has a negative correlation with turbulent dissipation. In this paper, the unstable top and bottom sections are not considered. When the total column height is 25 m, the turbulent dissipation significantly increases from the third section onwards. This  Figure 18 shows the turbulent eddy dissipation at different cross sections for various total column heights with the U Gr of 0.0389 m/s. When the gas holdup is constant, the bubble diameter has a negative correlation with turbulent dissipation. In this paper, the unstable top and bottom sections are not considered. When the total column height is 25 m, the turbulent dissipation significantly increases from the third section onwards. This is an important reason for the decrease in the average bubble diameter.  Figure 18 shows the turbulent eddy dissipation at different cross sections for various total column heights with the of 0.0389 m/s. When the gas holdup is constant, the bubble diameter has a negative correlation with turbulent dissipation. In this paper, the unstable top and bottom sections are not considered. When the total column height is 25 m, the turbulent dissipation significantly increases from the third section onwards. This is an important reason for the decrease in the average bubble diameter.

Conclusions
The hydrodynamic characteristics of an industrial-scale internal loop airlift column were studied based on a molten sulfur system with high viscosity using the three-dimensional transient Eulerian-Eulerian model coupled with the population balance model. The

Conclusions
The hydrodynamic characteristics of an industrial-scale internal loop airlift column were studied based on a molten sulfur system with high viscosity using the three-dimensional transient Eulerian-Eulerian model coupled with the population balance model. The effects of the U Gr and column height on the gas holdup, liquid circulation velocity, and average bubble diameter were examined.
As the total height of the column remains unchanged, with the U Gr mainly affecting the gas holdup, there was a critical value for the variation trend of liquid circulation velocity and bubble diameter. Under the determined superficial gas velocity, the column height had an obvious impact on the liquid circulation velocity, while the gas holdup had no evident change. When the total column height and superficial gas velocity were large, turbulent dissipation was a significant factor that affects the bubble diameter.
The CFD simulation results were slightly different from the results obtained using the Chisti and proposed correlations. The fitting accuracy of the CFD simulation results and correlations should consider the influences of the column geometry size and physical parameters simultaneously. The CFD simulation results provide detailed information regarding the hydrodynamics, which are beneficial for understanding the aeration behavior of an airlift column in molten sulfur.