Euler-Euler CFD simulation of high velocity gas injection at pool scrubbing conditions

Pressure relief by blowdown is one of the most important measures to prevent excessive pressures in the primary circuit or containment in severe nuclear accidents. Pool scrubbing can significantly reduce the release of radioactive materials, e.g., aerosols, to the environment during the pressure relief. The decontamination factor indicating the particle retention efficiency depends, among other factors, on the hydrodynamic conditions of the gas-liquid two-phase flow inside the pool. In the present work, the hydrodynamics in two typical pool scrubbing experiments is investigated with the two-fluid model, and the influence of some key factors including bubble diameter, nozzle submergence as well as interaction models are analysed. One case is a rectangular pool and the other is a cylindrical column, and their injection Weber number is around 2×103 and 4×105, respectively. The numerical results show that as the distance from the nozzle exit increases, the void fraction and velocity field expand from the central region, where the nozzle is located, to the whole cross section. The profile and its development depends largely on the bubble size and the interaction force model. It reveals that in the monodisperse simulation, the tuning of bubble diameter is necessary for achieving good agreement, although it is difficult for high velocity gas injection. More information is required to properly describe the bubble size distribution as well as its evolution in pool scrubbing conditions. Furthermore, the experimental data show clear drag reduction in the bubble swarm generated by the gas jet, and the mechanism and model improvement possibilities need to be investigated.


Introduction
During a severe accident in nuclear power plants, large amount of fission products may be released as aerosols and escape from the containment, which has a serious impact on the environment as well as the health of nearby residents. It is therefore vital to find reliable ways to eliminate these fission products effectively. To ensure the integrity of the containment and reduce the risk of radioactivity, the containment is equipped with a pressure relief vent. When the temperature and pressure inside the containment becomes too high, the venting facility is activated to reduce the pressure and minimize the release of radioactivity to the environment (Swiderska-Kowalczyk et al., 1996). The wet filter discharge system is one of the most internationally used filter discharge systems. The gases and aerosols in the containment are discharged through pipes into a water pool, which acts as a filter for the fission products and retains a part of them inside the water. This process is referred to as pool scrubbing. Besides the filtered containment venting system, pool scrubbing scenarios can be witnessed in some other components, e.g., the suppression pool in boiling water reactors and the secondary side of steam generators in pressurized water reactors (He et al., 2021;Ouallal et al., 2021).
It is of utmost importance to understand the aerosol removal efficiency of pool scrubbing as well as its influence factors. The retention efficiency is often characterized by the decontamination factor (DF), which is defined as the ratio between the aerosol particle mass that is injected and that escapes from the pool. DF depends on many parameters such as particle size and concentration, bubble size and distribution, nozzle submergence, pool depth, gas composition, etc. Numerous experimental studies on pool scrubbing have been carried out since the 1970s, including the ACE experiment at the American Electric Power Research Institute (EPRI) (Ramsdale et al., 1992), the LACE-Espana experiment at the Center for Energy, Environmental and Technological Research (CIEMAT) (Marcos et al., 1994), the POSEIDON experiment at the Paul Schell Institute (PSI), Switzerland (Hillary et al., 1966). More recently, Li et al. (2021) suggested that the DF varies exponentially with pool depth in particular for the retention of large particle size aerosols. For particle diameters in the range of 0.2-0.52 μm, it has almost no effect on the DF. For low-pool-depth scrubbing, the steam condensation mechanism dominates and particle diameter has little effect on the DF. Sun et al. (2019) measured the relationship between DF and aerosol particle number concentration in a 2.4 m deep water pool. Under constant thermal hydraulic conditions, the DF increased monotonically with decreasing concentration, and the dependence became more significant in deeper nozzle submergence. Kanai et al. (2016) studied the decontamination characteristics for BaSO 4 aerosol with air as the carrier gas, and derived an empirical formula through relevant experimental data. The results showed that the decontamination efficiency of aerosol was gradually enhanced with the increase of particle size.
The above experimental research mainly focused on the decontamination efficiency as well as the influence of some design parameters. The hydrodynamic conditions such as flow velocity, bubble size, and void fraction were not the center of concern, which are however key influencing factors. The pool scrubbing is a three-phase flow process involving complex interaction between aerosols, bubbles, and the liquid. The influence of bubble hydrodynamics on the retention of aerosol particles in the pool has gained widespread attention. With the increase in computational capability and algorithmic efficiency, the high-resolution computational fluid dynamics (CFD) method is able to provide detailed information on the flow field in technical multiphase systems. Although most numerical studies on pool scrubbing so far are based on system codes, bubble columns have been widely studied with three-dimensional CFD simulations (Ma et al., 2015;Ziegenhein et al., 2015). Compared to the Euler-Lagrange framework (Gouesbet and Berlemont, 1999;Zuzio et al., 2018), the Euler-Euler or two-fluid model (Ferry and Balachandar, 2001;Laborde-Boutet et al., 2009;Azadbakhti et al., 2019) is more commonly used to investigate the hydrodynamic behaviour in bubble columns, particularly at high phase fractions due to its low computational resource consumption. The key to using the Euler-Euler method is to describe the gas-liquid interphase forces (Swearingen et al., 2023), the liquid turbulence, and the bubble distribution (Tabib et al., 2008;Mühlbauer et al., 2019;Fan et al., 2021). Despite effort and progress over past decades, precise description of the interfacial physical processes by closure models is still greatly limited (Oey et al., 2003;Liao et al., 2018), and validation against experimental data is often necessary. While certain confidence has been obtained in low void fraction pipe flows and homogeneous bubble columns (Liao et al., 2019), columns with oscillating bubble plumes and high injection rates still pose challenges. In the presence of particles, hydrodynamics and bubble dynamics become even more complicated (Mühlbauer et al., 2021). In pool scrubbing, a liquid column or pool with submerged nozzle is the key component where jet injection regime and high momentum scenarios take place. Compared to analyses with system codes, CFD simulations on bubble column under pool scrubbing conditions are scarce. Recently,  studied the bubble hydrodynamics based on the volume of fluid mixture model at an unaffordable computational cost.
The objective of the present work is to evaluate the capability of two-fluid model and its closures in analysing the hydrodynamic conditions of pool scrubbing with the help of experimental data from the literature (Abe et al., 2018;Sun et al., 2019). The simulations are carried out in the open source CFD software OpenFOAM v8.0 with Helmholtz-Zentrum Dresden-Rossendorf (HZDR) addons . After a detailed description of the numerical method in Section 2, Section 3 briefly summarizes the experimental and numerical setup of the investigated cases. Section 4 discusses the major results of the simulations and compares them with the experimental data. Finally, a summary of the major findings and recommendations for future research concludes the paper.

Numerical methods
In the Euler-Euler two-fluid model, both discrete and continuous phases are considered as inter-penetrating continuum, and the interface is not resolved. The transportation and motion of each phase is described by a set of conservation equations, and the information on interfacial area density and transfer rates is reconstructed by means of closure models. The conservation equations of the two-fluid model have been discussed at length in a number of books (Drew and Passman, 1998;Ishii and Hibiki, 2010;Yeoh and Tu, 2019) and a broad consensus has been reached, so only a brief summary of the pertinent equations will be given here. Closure relations, in contrast, are still subject to considerable variations between researchers. In this work, the same set of models as a series of previous studies (Liao et al., 2018(Liao et al., , 2019 are adopted in line with the HZDR baseline closure concept (Lucas et al., 2016;Rzehak et al., 2017), and their reliability in describing the pool scrubbing scenarios is assessed with the aid of experimental data.

Two-fluid model equations
For adiabatic gas-liquid flows, the equations for mass and momentum conservation are expressed as where the subscripts G and L denote the gas and liquid phase, respectively, and α, , ρ p, g, and u denote the volumetric phase fraction, density, pressure, acceleration of gravity, and phase velocity, respectively.
The stress tensor T including molecular and turbulent stresses is defined as where eff μ and k represent the effective viscosity and the turbulent kinetic energy, respectively, and I is a 3×3 identity matrix. eff mol turb where mol μ and turb μ are the molecular and turbulent eddy viscosity, respectively. Note that in the present study the gas phase is assumed laminar. Turbulence effects are considered only for the liquid phase, including the bubble-induced turbulence (BIT).
The term F inter represents the sum of all forces between the phases as shown in Eq. (7): and describes the overall interphase momentum transfer. In addition to the closures for turbulence and interphase interactions, the gas phase is assumed to be present in the form of discrete bubbles and their sizes are constant and uniform, i.e., a monodisperse approach is adopted. The effect of the prescribed bubble diameter on the simulated void fraction and velocity distribution is investigated additionally.

Interphase forces
In the present study, the interphase drag force, lift force, wall lubrication force, turbulent dispersion force, and virtual mass force are considered, i.e., The modelling of each force is explained below.

Drag force
The form and skin-friction drag affects the rising velocity of the bubbles in the liquid phase, and consequently the magnitude of the void fraction, the residence and scrubbing of particles in the pool. The evolution of gas-liquid interfacial topology during the pool scrubbing is complex, especially at high injection velocities. Both gas and liquid may be present in either a continuous or a dispersed form. To be consistent with the physical picture, the drag force is modelled according to local flow regimes. For bubbly flow, where d B is the equivalent size of bubbles, which are assumed spherical. The drag coefficient C D is given by the correlation from Ishii and Zuber (1979), which distinguishes between different bubble shape regimes. D,ellipse D,cap D,ellipse D,sphere D D,sphere D,ellipse D,sphere min( , ), , where σ denotes the surface tension. For droplet flow regime, the drag force is modelled in a similar way, but the coefficient is evaluated with the Schiller-Naumann (Schiller and Naumann, 1935) correlation. In the mixed region, the drag force is computed by using the segregated model of Marschall (2011). The overall interface drag between gas and liquid is obtained by blending the various regimes linearly.

Lift force
The shear lift force describes the interaction between the rising bubble and the liquid shear field. The difference in velocity between the bubble and the liquid phase causes an asymmetric pressure distribution on the bubble surface, which results in a lateral migration of the bubble. It therefore has an effect on the distribution of void fraction and axial liquid velocity in the radial direction. In upward flows, it drives large bubbles to the center and small ones to the wall. In other words, the coefficient changes its sign at a critical bubble size.
The lift coefficient is calculated using the empirical model in Eq. (13), recently proposed by Hessenkemper et al. (2021).
ln 1 exp 4 5.6 0.11 0.14 5.2 0.44 4 where the major axis of the ellipsoidal bubbles is calculated according to the correlation of Ziegenhein and Lucas (2017) in Eq. (16): Note that the coefficients are slightly different from those in the widely used Wellek correlation (Wellek et al., 1966).

Wall lubrication force
The effect on the discrete phase due to the presence of the wall is referred to as wall lubrication force. For tube flows, the effect of wall on the void distribution and flow pattern cannot be ignored, in particular in the presence of small bubbles, which are accumulated in the near-wall region due to the effect of lift force introduced above.
u is the wall tangential component of the relative velocity between phases, ŷ is the unit vector normal to the wall pointing into the fluid. The coefficient W C is calculated with reference to Hosokawa et al. (2002).
where y is the distance between the bubble and the wall.

Turbulent dispersion force
The dispersion of bubbles due to liquid velocity fluctuation is described by the turbulent dispersion force. Burns et al. (2004) derived an explicit expression by Favre averaging the drag force as in Eq. (20): where TD σ is the turbulent Schmidt number and is set to 0.9.

Virtual mass force
The virtual mass force describes that a portion of liquid mass is accelerated due to the flow acceleration and velocity difference between the phases. According to Drew and Passman (1998), it may be formulated as Eq. (21): where the substantial derivative of a scalar or vector field  is defined as in Eq. (22), and the virtual mass coefficient vm C is set to 0.5.

Turbulence model
The k-ω SST is a blending model combining the standard k-ε model and the k-ω model, which combines the advantages of the k-ω model near the wall and the k-ε formulation in the bulk (Menter, 2009). The model solves the equations for the turbulent kinetic energy k and the turbulent eddy frequency ω. For the liquid phase, It transforms automatically between the k-ω model and the k-ε model via the blending function F 1 . In the area near the wall, it activates the k-ω model, while in the freestream as well as in the core it recovers the k-ε model. The blending function is shown in Eq. (25): It is also used to interpolate the model constants μ C , (1 ) ωP ωP ωP ,2 ωP C , 1 ,2 k σare coefficients of the k-ω model and the k-ε model, respectively. The default values of these constants are summarized in Table 1. Table 1 Values for k-ω and k-ε model constants The production term for the shear-induced turbulence: The strain rate tensor: There is a broad consensus in the literature (Kataoka et al., 1992;Troshko and Hassan, 2001;Liao and Lucas, 2012) regarding the source term in the k-equation describing the bubble effects. A reasonable assumption is that a part of the energy lost by the bubble due to drag is converted into turbulent kinetic energy at the wake of the bubble. Hence, the k-source becomes For the ε-source term, we use a heuristic model in analogy to the single-phase theory, i.e., the k-source divided by a time scale τ (Yao and Morel, 2004), so that Regarding the determination of coefficients k C , ε C , and time scale τ, there is still controversy (Rzehak and Krepper, 2013). In this work, the model of Ma et al. (2017) is adopted: which has been validated for vertical pipe flow and homogeneous bubble columns in the previous work (Liao et al., 2019;Liao and Ma, 2022). For the use of the SST model, the ε-source is transformed to an equivalent ω-source, which gives The turbulent viscosity is evaluated from the standard relation: where γ 1 / 0.31 C = is a further model constant, and F 2 is a second blending function:

Pool scrubbing experiments
In order to qualitify the two-fluid model and closures described above for investigating pool scrubbing hydrodynamics, the experiments reported by Abe et al. (2018) and Sun et al. (2019Sun et al. ( , 2021 are simulated, considering that comprehensive data are available for validation of the numerical results. In Abe et al.'s case, the gas is injected through a circular nozzle with an inner diameter of 6 mm into a transparent rectangular pool of 500 mm × 500 mm × 3000 mm. The material system is air-water and the nozzle is submerged to a depth height of 1100 mm. A high-speed camera is used to record the motion and shape of bubbles, while two-layer wire mesh sensors (WMS) are installed at heights of z = 100, 300, 500, 700, and 900 mm (distance to the nozzle exit) in order to measure the evolution of void fraction and gas velocity distribution profiles. Four groups of injection flow rates, Q  = 8.5, 42.4, 84.8, and 254.0 L/min, were investigated, corresponding to an injection velocity of 5.0, 25.0, 50.0, and 150.0 m/s, respectively. In the present work we choose the first case with Q  = 8.5 L/min for simulation, since the previous work of  shows that bubble coalescence and breakup is non-negligible at high flow rates. Furthermore, the experimental data show considerable deformation of bubbles. They mostly have an ellipsoidal shape, and the aspect ratio can reach 3.5. Figure 1 shows the comparison between the data and several correlations. The Ziegenhein and Lucas (2017) correlation (Eq. (16)), which is adopted in this study, provides a better agreement than that of Wellek et al. (1966) and Okawa et al. (2003), particularly for d B in the range of 2-8 mm (Fig. 1). Another case is taken from Sun et al. (2019Sun et al. ( , 2021, which is different from the last case in both geometry and injection velocity. Sun et al. investigated the relationship between DF and inlet number concentration of aerosol particles. The experiment was carried out using a pool scrubbing test apparatus named PONTUS (pool scrubbing test unit on separate effect). The test section, which is a circular pipe with 0.2 m inner diameter, is filled with deionized water. Using dry air as the carrier gas, the solid test particles are injected upwards through a 1.0×10 -2 m diameter nozzle at the bottom of the test section. The experiments analysed in the papers Sun et al. (2019Sun et al. ( , 2021 were conducted with insoluble, monodisperse spherical SiO 2 particles (Nippon Shokubai Co., KE-S50, 0.5 μm in diameter, amorphous silica) without any surface modification. The particle number density and diameter distribution in the airflow entering and leaving the test section were measured by an aerosol spectrometer and DF derived as their ratio. In addition to DF, the distributions of void fraction and Sauter bubble mean diameter are obtained by a four-sensor optical probe and a high-speed camera, both measured at 1.6 m above the nozzle and a water depth of 1.9 m. The injection flow rate at the nozzle exit was kept at 1.3×10 -3 m 3 /s, which corresponds to an injection velocity of nearly 17 m/s.  Sun et al. (2019Sun et al. ( , 2021 and Abe et al. (2018). Note that the present study focuses on gas-liquid hydrodynamics under pool scrubbing conditions, while its coupling with particles will be investigated in the near future.

Numerical settings
The computational domain for the Abe et al.'s case is a 500 mm × 500 mm × 1500 mm rectangular pool with the initial water height of 1100 mm. In order to reduce computational time, the domain is downsized by first performing a comparison study between a quarter and a full pool. As shown in Fig. 2, the time-averaged void fraction and velocity distributions obtained with a full pool and a quarter are nearly identical, and thus the domain with a quarter section of the pool is taken for further studies (see Fig. 3(a)).
For Sun et al.'s case, the circular pool height is 4500 mm, the water depth is 1900 mm, and the void fraction is monitored at H = 1600 mm. The flow field is assumed axisymmetric and numerical simulations are carried out for a wedge with an opening angle of 1° instead of the whole cylinder, as shown in Fig. 3(b).
To ensure gird-independent results, a comparison study among three grids is performed for each case, as shown in Table 2. In Abe et al.'s case, starting from the jet inlet, the mesh expands with a scaling ratio of 2 in the horizontal direction and 10 in the vertical direction. The scaling ratio is defined as the ratio between the last (largest) and the first (smallest) cell size. In Sun et al.'s case, the cell size is uniform in the radial direction while expands with a scaling ratio of 10 in the vertical direction. The transient simulations are run for 20 and 40 s for the two cases, respectively. Adjustable time step is used with an initial time step of 1×10 -6 s and the maximum Courant number is limited to 0.5.

Figures 4(a) and 4(b)
show the results of the grid independence test for the two cases, respectively. From the void-fraction distribution, it can be seen that there is no significant difference in Abe et al.'s case, while in Sun et al.'s case small deviation is observable particularly in the near wall and central regions, but the change from grid II to grid III is sufficiently small. Therefore, as a good compromise between accuracy and efficiency, grid II is used for further studies.
The multiphaseEulerFoam solver with HZDR addons, referred to as HZDRMultiphaseEulerFoam , was chosen to perform the simulations. The types of boundary conditions applied in the two simulations are shown in Table 3. In this work, since the numerical calculations were performed for a quarter of a rectangular pool for Abe et al.'s case and a wedge for Sun et al.'s case (see Fig. 4), symmetryPlane boundary condition was applied to the two cut planes in the former case, while wedge condition (axis-symmetry) was applied to the latter one. The injection of pure air through the nozzle was realized by fixing G 1 α = at the inlet, and the injection velocity was determined by the injection flow rate and nozzle cross sectional area, which is G (0,0,5.0) = u m/s and G (0,0,16.552) = u m/s for Abe et al.'s case and Sun et al.'s case, respectively. A pressure boundary was applied at the outlet and the static pressure at the outlet boundary was prescribed 0 Pa with the boundary condition prghPressure. The wall was assumed adiabatic and smooth, at which the liquid and gas phases have no-slip and free-slip conditions, respectively. Due to the lack of reliable models and data for evaluating the bubble size distribution, we assumed a fixed monodisperse size distribution for both cases and analysed the effect of the bubble size on the results additionally. The modelling of bubble coalescence or breakup processes as well as morphology transfer in pool scrubbing will be added in the next step.  I Ⅱ Ⅲ Total

Numerical results
Void fraction is a key factor reflecting the hydrodynamic behaviour in pool scrubbing. The position and motion of bubbles cause a density and velocity difference in the radial direction of the pool, which drives the formation of an overall circulation flow. It further affects the residence time of aerosol bubbles and particles inside the pool and the retention efficiency. Because of the coupling between the void fraction and momentum transfer, the bubble diameter has a significant impact on the distribution. In Abe et al.'s case, the experiment revealed that due to coalescence and breakup phenomena, there is a spectrum of bubble sizes present in the pool at high injection conditions. Unfortunately, no data were available for determining the mean bubble diameter or size distribution properly. In the present study, bubble coalescence and breakup were assumed in equilibrium and a constant bubble diameter was prescribed. According to Akita and Yoshida (1974), the fully-developed Sauter mean bubble diameter far away from the nozzle in a bubble column depends on column diameter, gas superficial velocity, liquid properties (density and viscosity), surface tension, and gravitational acceleration.
where D c and J G denote the column diameter and superficial gas velocity, respectively. In the case of a rectangular pool, c D is set to hydraulic diameter. For Abe et al.'s case with an injection rate of 8.5 L/min, B d ≈ 6.75 mm, while for Sun et al.'s case it was approximately 5.18 mm, which is close to the measurement shown in Fig. 3(b) in Sun et al. (2021). Nevertheless, for Abe et al.'s case the recent study of  gives a value of 13.7 mm, which was obtained by averaging all bubbles in the whole domain. i.e., including large globules formed near the nozzle. One can see that the definition of a single bubble diameter has large uncertainty in pool scrubbing cases, particularly in the developing zone. The influence of the bubble diameter is investigated below. Figure 5 shows the gas-jet structure and the comparison of instantaneous and time-averaged void fraction and gas velocity. It can be seen that the jet is confined in the central region, and the largest void fraction and gas velocity is located near the nozzle at the bottom. Because of momentum exchange with the liquid phase, the jet spreads gradually in the lateral direction, and the peak of void fraction and velocity decreases rapidly. A relatively stable plume is formed above the nozzle region up to the free surface, where a bulge of liquid is observed due to strong interaction.

Abe et al. 's case
The comparison between the instantaneous and timeaveraged void fraction and gas velocity reveals that the gas jet at t = 20 s is almost in a steady state, since they are  similar in most regions. Nevertheless, obvious difference exists near the nozzle. The instantaneous void fraction distribution shows that the gas enters the pool in forms of large discrete structures, which is consistent with the experimental observation. The high-speed images reveal that the globules experience breakup shortly after their detachment and form a swarm of smaller bubbles. Since the size of these globules is much larger than the mean diameter prescribed in the monodisperse study, precise description of their hydrodynamics is difficult for two-fluid model. In the present work, a multi-scale concept was considered in the interphase drag modelling to overcome this difficulty. In the future, a hybrid framework, e.g., the one recently developed by Meller et al. (2021), can be adopted to improve the capture of morphology transfer in the pool scrubbing. As shown in Fig. 6, one can see that the void-fraction distribution is closer to the experimental observation by assuming a constant bubble diameter of 4.30 mm. The void-fraction distribution has an obviously too high central peak if a larger value is used for B d , in particular at the vertical position z = 100 mm (see Fig. 6(a)). This is understandable since the lift force pushes the large bubbles  to the center and the effect increases with the bubble diameter. It is worth noting that in the injection region, where the bubble is not yet detached from the nozzle, the two-fluid model assuming spherical bubbles may cause large uncertainty.
As the height z increases, the sharp void fraction gets flattened and wider toward downstream, which is consistent with the overall trend of void-fraction distribution observed in the experiment. Note that for better representation the horizontal coordinate is limited to 100 mm, although the pool has a half width of 250 mm, because in the region outside x = 100 mm the void fraction is nearly zero.
The radial distribution of the vertical velocity of bubbles is approximately the same as that of the void-fraction distribution, as shown in Fig. 7. Because of the injection at the center, most of the bubbles are accumulated in the central region, although the distribution width increases as they rise away from the nozzle exit. The measured velocity in the majority of the peripheral region is zero, since there are no bubbles present there. The predicted profile is similar to the measured one, but the peak is obviously over-predicted even with d B = 4.3 mm. At z = 100 mm height position, the velocity peaks of d B = 6.75 mm and d B = 8 mm are significantly higher than that of d B = 4.3 mm. From z = 100 mm to z = 900 mm, the velocity distribution of d B = 4.3 mm is more consistent with the experimental data.
The reduction of bubble diameter helps to improve the results, but the agreement with the data is still unsatisfactory. Furthermore, d B = 4.3 mm is much smaller than the value provided by the empirical relation of Akita and Yoshida (1974) and the volume-of-fluid prediction of . Apart from the bubble size, drag coefficient is another key parameter affecting the void fraction and velocity distribution. Based on the experimental data for bubble size and rise velocity, an average drag coefficient of around 0.5 was evaluated, which is much smaller than the value given by the Ishii and Zuber (1979) drag model (see Fig. 8).
For a comparison study, the Schiller-Naumann (1935) drag model was adopted in a simulation for bubbly flow regime. Although limited for spherical bubbles or particles, it predicted a constant drag coefficient of around 0.44 in this case, which was close to the measured mean value (see Fig. 8). The effect of drag models is shown in Fig. 9. Compared with the Ishii-Zuber model, the void-fraction profile obtained using the Schiller-Naumann model had a smaller peak and wider distribution, which was more consistent with the experimental data. Correspondingly, the velocity distribution of gas phase in the radial direction was broader. One can speculate that specifying a bubble size larger than 6.75 mm will improve the agreement further, which conforms to the results of . Nevertheless, it is known that the applicability of the Schiller-Naumann model is limited   Ishii and Zuber (1979) and Schiller-Naumann (1935) to the experimental data.
to spheres, so the mechanism of drag reduction in the swarm of ellipsoidal bubbles under pool scrubbing conditions requires further investigations.
The effect of turbulence modelling was investigated as well, by performing an additional simulation using the standard k-ε model. The comparison in Fig. 10 shows a very small difference, except for the void fraction at the pipe center, where the turbulence is the strongest due to central injection. This is because bubble-induced turbulence plays a dominant role in the pool scrubbing case, which is accounted for in the two simulations with the same model from Ma et al. (2017) (see Eqs. (28)-(30)).

Sun et al. 's case
Compared to the inlet conditions in Abe et al.'s case where the gas injection velocity was 5 m/s, the gas injection velocity in Sun et al.'s case was as high as 17 m/s. As can be seen in Fig. 11, the increase in gas injection velocity leads to more large gas structures. The transient void fraction distribution at t = 40 s shows that the globules detached from the nozzle, merged with the preceding bubble, and formed even larger bubbles, which rose up to the middle of the column and then destructed and smeared. Although the time-averaged distribution seems continuous, the instantaneous void fraction and velocity fields display clear morphology change, Fig. 9 (a) Time-averaged void-fraction distribution and (b) velocity distribution of gas phase at z = 500 mm for different drag models (dB = 6.75 mm).

Fig. 10
Effect of turbulence models on (a) time-averaged void-fraction distribution and (b) velocity distribution of gas phase at z = 500 mm (dB = 6.75 mm, Schiller-Naumann drag model).

Fig. 11
Comparison of instantaneous and time-averaged void fraction (left) and gas velocity (right) (note that the sketch is not to scale).
interface, and discontinuity. These features need to be revisited when choosing numerical methods and closure models for future studies. In the present work, a multi-scale strategy was applied to drag modelling. Furthermore, because the height of the column was much higher than that of Abe et al.'s pool, the void profile in the radial direction had enough time to develop. Before reaching the free surface, the bubbles had been redistributed over the whole section and a fully developed profile was developed, and consequently the interaction between the gas jet and the free surface was gentler compared to Abe et al.'s case. Sun et al.'s DF experiments (Sun et al., 2019) were carried out for different nozzle submergences, while the hydrodynamic parameters (void fraction and Sauter mean diameter) were measured only for a submergence of 1900 mm (Sun et al., 2021). The simulation was set up according to the test geometrical and injection conditions. In order to find out whether the liquid level has an effect on the void-fraction distribution and velocity distribution, two simulations with liquid levels of H = 1900 and 2400 mm were performed. As shown in Fig. 12, the void-fraction and velocity distribution at position z = 1600 mm was almost the same for the two liquid levels. Therefore, the initial liquid level of H = 1900 mm was assumed in further studies. The shape of the velocity curve conforms to a turbulent pipe flow, which is essentially zero at the solid wall and increases toward the center. A rather flat and uniform distribution was obtained for both velocity and void fraction with d B = 5 mm according to the data (see Fig. 12).
The experimental data (from Fig. 3(b) in Sun et al. (2021)) showed a cross-section averaged Sauter mean diameter around 5 mm. Since the void-fraction distribution obtained with a constant bubble diameter of 5 mm in the present simulation differs significantly from the experimental data, the uncertainty brought by d B was investigated. We increased the bubble diameter to 7 mm and found out that the void-fraction distribution became consistent with the experiment, which took a core-peak instead of a uniform or wall-peak profile as shown in Fig. 13. It indicates that the mean size of the bubbles in the column is larger than the critical value, at which the lift force coefficient changes its sign. As mentioned above, the lift force resulting from the interaction of rising bubbles with the liquid shear field caused the lateral accumulation of bubbles. It drove small bubbles towards the wall and large bubbles to the center, forming a wall-peak and core-peak of void fraction distribution, respectively. If the size of the bubbles is intermediate, they will distribute uniformly over the cross section. The velocity profile is correspondingly much steeper at d B = 7 mm than that at d B = 5 mm, and validation against experimental data is desirable. The study revealed that the critical bubble size for lift force to change its sign is smaller than the one predicted by the Hessenkemper et al. (2021) model. In the majority of references describing the experiments, detailed information about the nozzle, e.g., the wall thickness, the distance between the exit and the pool bottom, is missing, which inhibits the evaluation of these effects. As a result, the nozzle is routinely simplified as a hole in the pool bottom in numerical studies such as in Okagaki et al. (2020), Liao and Lucas (2020), Bicer et al. (2021), and . In this work, the difference between considering a nozzle part below the pool and simplifying it to a hole in the pool bottom was analysed. As shown in Fig. 14, there is slight difference in the two configurations, especially in the central region, where the void fraction is higher while the velocity is lower without considering the nozzle. The radial distribution of the void fraction when considering the nozzle is closer to the experimental values. More geometrical details on the nozzle in the experiment are required for further analyses.

Conclusions
Because of its function in alleviating the release of aerosol particles, pool scrubbing has gained increasing interest in the past years. Compared to laboratory experiments, CFD investigations are scarce and the capability of numerical approaches still requires assessment. The understanding of the gas jet structure and hydrodynamics in pool scrubbing is far from satisfactory. The purpose of the present work is to gain insight into submerged gas jets under pool scrubbing conditions using the Euler-Euler two-fluid model, which is mostly suitable for practical scenarios due to high computational efficiency. The evolution of the void fraction and velocity distribution along the jet was investigated in OpenFoam for two cases, taken from the experiment of Abe et al. (2018) and Sun et al. (2019Sun et al. ( , 2021, respectively. In both cases, the gas was injected from the middle of the bottom plane of the pool. The detached bubbles built an inverted conical plume above the nozzle exit, and diffused from the center to the wall of the pool. After certain developing length, a stable core-peak or wall-peak profile would be formed depending on the bubble size. Based on the comparison with experimental data, the following conclusions and perspectives were obtained. (1) In Abe et al.'s case, both numerical and experimental results showed that the diffusion of void fraction from the pool center to the wall is not accomplished up to a distance of 900 mm away from the nozzle. The speed of redistribution  depends on the bubble size and gas-liquid interaction forces. A smaller bubble size indicates a smaller negative lift coefficient or a positive coefficient if it exceeds the turning point, which helps to speed up the void redistribution. A smaller drag coefficient results in a higher bubble rise velocity, and likewise the development of the radial profile is accelerated.
(2) In Sun et al.'s case, data on the void distribution are only available for the position of 1600 mm above the nozzle exit, which revealed a fully developed profile. The CFD simulation is able to reproduce the core-peak profile by assuming a bubble diameter of 7 mm, while the measured Sauter mean diameter around 5 mm gave a rather flat profile. It indicated that the turning point between positive and negative lift force was not satisfactorily predicted. In addition to bubble size, data on the velocity distribution are desirable for further assessment of the numerical models and assumptions.
(3) It was found that the development of hydrodynamic parameters like void fraction and velocities depends on the distance from the nozzle, while the influence of the distance to the free surface is negligible. However, the nozzle submergence is still an important parameter affecting the DF value by increasing the bubble residence time in the pool and alleviating the interaction with the free surface. Furthermore, the geometrical simplification of the nozzle made in the simulation was shown to have an effect on the results.
(4) The data on bubble size and rise velocity available in Abe et al.'s experiment showed clear drag reduction compared to the prediction of the Ishii-Zuber model, which is widely used for turbulent bubbly flows. The mechanism of drag reduction and swarm effect should be a focus of future studies.
(5) The bubble size was shown to be a key parameter in predicting the pool scrubbing hydrodynamics. The monodisperse approach by tuning the bubble size is incapable of fully describing the interphase interactions. A spectrum of bubble sizes may be envisaged in the process of pool scrubbing, and poly-dispersed simulations with appropriate bubble coalescence and breakup models are necessary.
(6) In addition, high velocity injection generates large gas structures near the nozzle, which destruct into a swarm of bubbles during their rise through the pool. This feature rationalizes the use of a morphology-adaptive hybrid approach, which resolves the gas elements that are much larger than the grid size.