Modeling of the Free ‐ Surface Vortex ‐ Driven Bubble Entrainment into Water

: The recently developed GENTOP (Generalized Two Phase Flow) concept, which is based on the multifield Euler ‒ Euler approach, was applied to model a free ‐ surface vortex—a flow situation that is relevant for hydraulic intake. A new bubble entrainment model has been developed and implemented in the concept. In general, satisfactory agreement with the experimental data can be achieved. However, the gas entrainment can be significantly affected by several parameters or models used in the CFD (Computational Fluid Dynamics) simulation. The scale of curvature correction 𝐶 (cid:3046)(cid:3030)(cid:3028)(cid:3039)(cid:3032) in the turbulence model, the coefficient in the entrainment model 𝐶 (cid:3032)(cid:3041)(cid:3047) , and the assigned bubble size to be entrained have a significant influence on the gas entrainment rate. The gas entrainment increases with higher 𝐶 (cid:3046)(cid:3030)(cid:3028)(cid:3039)(cid:3032) values, which can be attributed to the stronger rotation captured by the simulation. A smaller bubble size gives higher gas entrainment, while a larger bubble size leads to a smaller entrainment. The results also show that the gas entrainment can be controlled by adjusting the entrainment coefficient 𝐶 (cid:3032)(cid:3041)(cid:3047) . Based on the modeling framework presented in this paper, further improvement of the physical modeling of the entrainment process should be done.


Introduction
A free-surface vortex (see Figure 1) may exist in a wide range of scales; it can be as small as a "bathtub vortex" [1][2][3] or can be as big as an ocean whirlpool [4]. The topic of a free-surface vortex is often found in the discussion of hydropower plants, nuclear reactors and other applications using pumps. The supply of water for irrigation, domestic, industry, and power generation is usually taken from rivers or reservoirs through an intake that is located near the surface [5]. Insufficient submergence (a short distance between the water surface and an intake) may lead to the formation of a free-surface vortex that can induce gas entrainment into the intake [5].
A free-surface vortex and its associated gas entrainment may lead to several operational and safety problems [5][6][7][8][9]. They may cause mechanical damage and loss of performance in fluid machinery such as turbines and pumps [6,7]. A swirl in a sump leads to rotational flow in a pipe, which may reduce the performance of the pump [8,9]. If such flow is unsteady, it may also cause fluctuating loads on pump bearings [8]. The gas entrainment induced by a free-surface vortex will reduce the delivery of a pump (1% air reduces the efficiency of a centrifugal pump by 5-15%) [6]. This reduction may cause a severe problem such as the overtopping of a dam [3], which may lead to a safety hazard and cause loss of life [7]. Vortex-induced gas entrainment is also an important issue for more specific applications such as nuclear reactors. In sodium-cooled fast reactors (SFRs), an inert cover gas such as argon is used and maintained above the sodium surface to accommodate the volume changes of sodium and prevent the contact of sodium with air [11,12]. Several mechanisms could lead this cover gas to entrain the sodium region, including the entrainment driven by a free-surface vortex [12]. The gas entrainment may cause safety and operational problems and, for this reason, becomes an essential issue in SFR safety analyses [13][14][15][16]. This issue has been intensively investigated by the Japan Atomic Energy Agency (JAEA) and its institutional partners [16][17][18][19][20][21][22][23][24][25][26][27][28][29][30]. The problems associated with gas entrainments in SFRs are: changes in reactivity when the gas reaches the core [11,14,18,[31][32][33], burnout of the fuel pin due to the trap of large bubbles [31], thermal stresses in the reactor structure [32,33], pump cavitation and fluctuations in pump discharge [14,33], reduction in the heat transfer efficiency [11,33], disturbance in the prompt detection of fission products leakage from failed fuel pins [34], disturbance of electromagnetic sensors used for shutdown systems [34], and trouble with acoustic or ultrasonic instrumentation such as boiling noise detector [31,32,34]. The gas entrainment issue is not only applicable to SFRs, but also to Pressurized Water Reactors (PWRs) and Boiling Water Reactors (BWRs). During mid-loop operation, gas may entrain into the Reactor Coolant System (RCS) of PWRs due to a free surface vortex and then be sucked into the Decay Heat Removal System (DHRS), which finally leads to a disturbance of the instrumentation [35]. A total failure of DHRS may occur if the void fraction exceeds 15% [35]. In the BWRs, the gas entrainment due to a free surface vortex may occur at the suction inlet from a condensation chamber/wet-well [36,37].
Computational fluid dynamics (CFD) may help to design a safer process, minimizing the aforementioned risks associated with gas entrainment due to a free-surface vortex. Generally, the previous CFD works available in the literature can be divided into two parts: single-phase and twophase computations. In a single-phase simulation, the deformation of the free surface is not considered and the free surface is defined as a free slip boundary [38,39]. When the mesh resolution is sufficient and the appropriate turbulence model is used, the velocity fields can be calculated using single-phase simulation, as reported by [17,40]. However, to the best of our knowledge, the estimation of the gas entrainment rate has never been performed by a single-phase CFD. In addition, the direct observation of a free-surface vortex from the single-phase simulation is not possible. A post-processing method is required to judge the occurrence and the location of the vortex, e.g., Q criterion [38]: where is the second invariant of velocity gradient tensor, is the vorticity tensor, and is the strain rate tensor. The above equation states that a free surface vortex exists when the strength of rotation is bigger than the local strain rate [38].
In the case of two-phase simulation, usually, the volume of fluid (VOF) model is employed [20,40,41]. The computation is performed in a fixed grid solving only one momentum equation, which is shared by both fluids [42]. Generally, a very fine mesh is required to resolve the interface, e.g., [40]. In the case of bubble entrainment driven by the free-surface vortex, the entrainment process needs to be resolved. This requires a very fine mesh, especially at the tip of the gas core. To reduce the number of computational cells used in the simulation, one may refine the mesh only in the tip region. However, to apply the local refinement, the precise location of the tip should be known before the simulation. This will be difficult to realize in practice since the tip location is transient in nature [43]. Therefore, a fine mesh should be used in a larger predicted region or over the whole computational domain instead of only at the tip region, which leads to very high computational costs. When there is a large computational domain such as a hydropower plant, this becomes a significant disadvantage of the method. An option to overcome this problem is using the Euler-Euler model, wherein the gas and liquid have their own sets of momentum equations. The possible advantage is that a coarser mesh may be used to model this gas entrainment. However, appropriate closure models need to be used, including the gas entrainment model.
The GENTOP (Generalized Two Phase) concept, based on the Euler-Euler model, which has been developed at Helmholtz-Zentrum Dresden-Rossendorf (HZDR), aims to handle various flow conditions where multiscale interfacial structures exist [44]. The potency of this concept has been demonstrated in several types of flows such as impinging jet [44], bubble column [44], a dam-break case with obstacle [45], churn turbulent flow in the vertical pipe [46], and boiling case [47,48]. However, to increase the effectiveness, continuous improvement and validation of various flow conditions are required. For that purpose, the applicability of the GENTOP concept to the flow case of bubble entrainment driven by a free-surface vortex was investigated in this work. A new bubble entrainment model considering the physics of the flow has been developed and implemented in the GENTOP concept in this study.

The GENTOP Concept
The Euler-Euler model was used in this work, whereby the following continuity and momentum equations were solved across the fixed computational mesh [49]:

∇.
(2) In the above equations, the liquid and gas phases are assumed to be adiabatic and incompressible. The subscript j represents the phase, while , , , t, , and p are the volume fraction, density, velocity vector, time, mass source, and pressure, respectively. and denote the interfacial forces and the momentum sources (i.e., due to external body forces), respectively [49].
The flow is represented by a continuous liquid phase, one velocity group for the polydispersed gas phase (dg) and one for the potentially continuous gas (cg) in the frame of the GENTOP concept [44]. The closure models are assigned specifically for each of the gas fields (i.e., dg and cg). For the dispersed gas phase, dg, the closure models based on the baseline model concept [50] are used as given in Table 1. These closure models influence the simulation results obtained in the complex rotating gas-liquid flow case, as discussed, e.g., in [51,52]. To avoid numerical instability problems, the virtual mass force (Equation (17)) was not used in this study. Table 1. Closure models for the polydispersed gas phase based on the baseline model concept [50]. The table is taken from [52].

Force
Formulation Ref. No.
, , , Wall lubrication Turbulent dispersion Virtual mass [61][62][63] (17) 0.5 (18) The detection of the statistically resolved large interface between gas and liquid is done in the GENTOP concept by the free surface function , which is formulated as follows [44]: where ∇ is the volume fraction gradient of the potentially continuous gas and ∇ is its critical value, which is given by [44]: The blending coefficient and were set to 100 and 4, respectively [44]. To support the transition from dg to cg by the agglomeration of cg when it reaches a critical amount, a so-called clustering force is used in the GENTOP concept [44]: where is the coefficient that defines the strength of the force [44]. The coefficient was set to = 0.1, similar to the value used by [44] for the impinging jet case. The blending function is used to define the location where the force should be activated [44]: where is the blending coefficient, set to 20 [44]. The force is blended in when is equal to 0.2 (i.e., , = 0.2) and blended out when it exceeds 0.8 (i.e., , = 0.8) [44]. To account for the surface tension, the following surface tension model of [64] that was implemented by [46] in the GENTOP concept was used: The interfacial area density and drag coefficient for cg in the GENTOP concept are calculated, respectively, as [44]: where is a blending function, formulated as follows [44]: where , is the critical volume fraction of cg, set to 0.3 [44]. The interfacial area density of bubble , and free surface , are calculated as follows [44]: , The drag coefficient of bubble , is calculated using Equation (6), while that of the droplet and free surface are given by [44]: Since the focus of this present work was to investigate the capability of the newly proposed entrainment model (described in the next section) to model the bubble entrainment driven by a freesurface vortex, bubble break-up and coalescence models as well as the complete coalescence model in the GENTOP concept were deactivated in this study. It is also important to note that liquid droplets were not considered in this work since they were not observed or measured in the experiment.

The New Entrainment Model
The pressure balance for a stable gas core is presented in Figure 2. The hydraulic pressure in the liquid phase is higher than in the gas core due to the higher density of the liquid. In the stable condition, the vortex gas core in the rotational flow does not collapse due to the contribution of centrifugal force, which counteracts the pressure from the liquid side. The pressure balance can be formulated as follows: where ∆ is the distance to the liquid level and is the centrifugal force density, which can be calculated based on the forced vortex model as follows: . (31) where is the vorticity magnitude and is the radius of the gas core at the equilibrium state. Substituting Equation (31) to Equation (30) yields: .
(34) Figure 2. Pressure balance for a stable gas core.
The underlying assumption of this entrainment model is that the detachment of the bubble from the gas core depends on the relation between the radius of the gas core related to the turbulent length scale : The entrainment is assumed to be proportional to the ratio between turbulent length scale and the gas core radius : The gas entrainment per unit volume and time then can be formulated as follows: The gas entrainment model is active only when . is a dimensionless coefficient and is the liquid vorticity magnitude. In the GENTOP concept, the entrainment model was then implemented as a source for the continuity equation as follows: This entrainment model controls the conversion of continuous gas cg into dispersed gas dg in the region detected as "interface" in the GENTOP concept, which is defined by . The gas entrainment model formulated in Equation (37) is different from the entrainment model proposed by [65], which has the formulation: In the proposed entrainment model described in Equation (37), the vorticity magnitude is used to take into account the rotation and to identify the location of the gas core tip, which is not included in the entrainment model of [65].

The Turbulence Model
The two-equation k-ω-based shear stress transport (SST) model proposed by [66] and the dispersed zero equation model are usually used for the liquid and gas phase, respectively, in the GENTOP concept [44,45]. However, the eddy viscosity models are insensitive to streamline curvature and system rotation [49] due to their isotropic nature [39]. Therefore, the curvature correction implemented in the SST model by [67] was used for the liquid phase in this work. The correction function is used in the model to control the turbulence production term as follows [49]: where is the empirical function proposed by [68]: in Equation (41) is the scale of curvature correction available in ANSYS CFX, which can be adjusted [49]. The arguments * and ̃ in Equation (43) are given by [49]: * where and represent the strain rate and vorticity tensor, respectively [49]: , 0.09 .
The empirical constants , , and in Equation (43) are set to 1.0, 2.0, and 1.0, respectively [49]. To consider bubble-induced turbulence, additional source terms based on [69,70] were used in the turbulence model.

CFD Set-Up
For the present work, an experimental dataset of bubble-type entrainment of the Ikoma experiment [71,72] was selected. The computational domain for this CFD study was adapted from the work of [72] with a slight modification to improve the boundary condition for the outlet. The computational domain, the mesh, and the defined boundaries are presented in Figure 3. The geometry consists of three major parts: the cylindrical test section, the connection pipe, and the bubble catcher or the buffer tank. In the Ikoma experiment, water enters the cylindrical test section through the tangential inlet and leads to the formation of a free-surface vortex [71]. In the case of bubble-type entrainment, bubbles are detached from the gas core and then entrained into the connection pipe [71]. Finally, the entrained bubbles are separated from the water in the buffer tank [71].
To investigate the influence of mesh size on the gas entrainment rate, two different meshes, namely Mesh A and Mesh B, which consist of 131,712 and 408,728 hexahedral cells, respectively, were used in this work. The meshes are refined in the central region and at the bottom of the cylindrical vessel, at the gas-liquid interface and in the connection pipe (see Figure 3a). The liquid level in the cylindrical test section ℎ and the exit velocity from the vessel are 60 mm and 0.66 m/s, respectively. The definition of this exit velocity is the superficial velocity of water in the connection pipe (see Figure 3b). The water flow rate can be calculated as follows: where D is the diameter of the outlet pipe (D = 8 mm). From the calculation, it is known that the water flow rate is 2 L/min. The inlet velocity can be determined by the following equation: , where is the width of the inlet slit ( = 10 mm). This results in being equal to 0.055 m/s. can be calculated in a similar manner, resulting in a velocity of 0.33 m/s. The boundaries that were applied to the computational domain are shown in Figure 3b. A uniform liquid velocity of 0.055 m/s was defined at the inlet. An opening boundary condition was applied at the top of the cylindrical vessel, while a degassing boundary was used at the top of the buffer tank. Uniform velocity boundary condition of = 0.33 m/s was applied at the outlet. All other boundaries were defined as no-slip wall and free-slip wall for the liquid phase and gas phase, respectively. Resolving the viscous sublayer is not required in the simulations since a wall function assuming a smooth wall is used [73].
In this study four bubble classes, G1-G4, having a diameter from 1 mm to 7 mm, are defined within the dispersed gas velocity group dg, while the last group G5, having a diameter of 9 mm or larger, is classified as continuous velocity group cg (see Table 2). The investigations were performed in several steps. First, the investigation of the influence of the scale of curvature correction was carried out by performing simulations using various . For this investigation, a fixed entrainment coefficient = 2E-4 was used to control the mass transfer from cg to G1. Next, the influence of the entrained bubble size distribution was investigated by performing three separate simulations, namely Sim. A, Sim. B, and Sim. C (see Table 2). In Sim. A, only the smallest bubble class G1 is entrained from the potentially continuous gas. In Sim. B, G1 and G2 with equal entrainment fractions are assigned to be entrained. In Sim. C, all bubble classes are assigned to be entrained with uniform entrainment fraction. The scale of curvature correction = 1 and the entrainment coefficient = 2E-4 were used in this step. The investigation of the influence of the entrainment coefficient was also carried out by performing simulations with two different . The mass transfer from cg to G1 with the scale of curvature correction = 1 was set in the simulations. Finally, the influence of the mesh size was also investigated and presented in the last section.

The Influence of the Turbulence Model
The interfacial deformation and bubble entrainment using various with a fixed value of = 2E-4 are presented in Figure 4a and 4b, respectively. In the case of = 0, the free surface is almost flat. At this condition, almost no bubble entrainment can be observed. Increasing to 0.5 leads to dimple formation. A spot of the dispersed gas fraction is observed, indicating that the entrainment model is activated. However, no bubble entrainment (with equal to or more than 0.02) into the connection pipe can be observed. Increasing to 1 leads to the formation of the gas core. At this condition, the bubble entrainment into the connection pipe is observed, confirming that the entrainment model successfully takes into account the bubble entrainment phenomenon. To have a quantitative description regarding the influence of on the bubble entrainment, the transient plot of continuous, bubble, and total gas entrainment are shown in Figures 5 a-c, respectively. Those plots confirm that almost no gas entrainment is produced other than = 1. The averaged gas entrainment rate shown in Figure 5d shows that the gas entrainment using = 1 is mostly obtained from the entrainment of the dispersed gas phase. This is the expected condition since the bubble type entrainment is observed in the experiment. Figure 5e shows that the gas entrainment rate using = 1, which is the default value in ANSYS CFX, is the most suitable value since it leads to the entrainment rate closest to the experimental value. To reveal the influence of on velocity fields, the velocity vectors represent the tangential projection of the liquid velocity are evaluated on horizontal planes H1, H2, and H3, which are located at an axial distance of 10, 20, and 30 mm from the bottom surface of the cylindrical vessel, respectively. These planes are in the liquid region to avoid the distraction of the velocity fields due to the presence of the gas core. Figure 6 presents the velocity vectors and the liquid vorticity magnitude as the background color. The velocity vector confirms that for all cases the flow is rotated and the center of rotation is not located precisely at the center of the cylindrical vessel. The magnitude of the rotation represented by the vorticity is largely affected by the selected used in the simulation. Increasing leads to a higher vorticity value that contributes to a higher gas entrainment. This can be understood since, according to Equation (37), the gas entrainment rate is proportional to the square of the vorticity magnitude. It can also be observed from the figure that the vorticity magnitude increases as the distance from the bottom of the vessel or the intake pipe decreases.

The Influence of Entrained Bubble Size Distribution
The contour of dispersed gas fractions at t = 6 s plotted on the central vertical plane is presented in Figure 7 for the different simulations listed in Table 2. The results of Sim. A clearly show bubbles being entrained into the connection pipe and then accumulating in the buffer tank (see Figure 7a). Although bubble entrainment can also be observed in Sim. B, it was not as significant as in Sim. A (see Figure 7b). Some of the dispersed gas still collected in the top of the gas core. In Sim. C, the dispersed gas mostly accumulated near the free surface close to the gas core instead of moving downward to the intake pipe (see Figure 7c). This may be because the terminal velocity for larger bubbles is bigger and the downward liquid velocity is not high enough to force bubbles to move into the connection pipe. The fact that some of the bubbles produced in Sim. B and Sim. C did not merge with the potentially continuous gas above the free surface is due to the inactivation of the complete coalescence model. This also shows that the complete coalescence is an important feature in GENTOP to bring the modeling of this flow behavior closer to the real condition, as we know from the literature (e.g., [13]) that there is a condition in which some of the entrained bubbles return to the free surface instead of flowing to the outlet. This inactivation in the simulations elucidates the role of the entrainment model alone.  The transient profile of the potentially continuous gas cg, the dispersed gas dg, and the total gas entrainment for all simulations are presented in Figure 8a-c, respectively. The transient profile of cg entrainment of all simulations shows that cg entrainment is negligible. The transient profile of dg entrainment shows a relatively stable curve after 4 s of simulation in the case of Sim. A and Sim. B. The dg entrainment in Sim. A is always higher than in Sim. B. In contrast, dg entrainment in Sim. C is always negligible. Figure 8d shows the averaged entrainment of cg and dg. This plot confirms that the entrainment of the dispersed gas phase contributes most of the gas entrainment. The total gas entrainment obtained in Sim. A is around 2.6 times that obtained in Sim. B. Sim. A gives the best fit to the experimental value (see Figure 8e). The results show that the entrainment rate that is controlled by the entrainment model is highly influenced by the entrained bubble size distribution. In the present study, the smallest bubble class, G1, is suitable to obtain a higher entrainment.

The Influence of the Entrainment Coefficient
The interfacial deformation and bubble entrainment using two different entrainment coefficients are presented in Figure 9a and 9b, respectively. The isosurface of the potentially continuous gas shows no significant difference in the interfacial deformation between those two cases. In contrast, larger bubble entrainment is observed in the case of = 2E-4. To obtain a quantitative description, the transient profile of potentially continuous gas cg, dispersed gas dg, and total gas entrainment for simulations using two different values of are presented in Figure 10a-c, respectively. The plot of cg entrainment shows almost no entrainment obtained from cg for both simulations. In contrast, dg entrainment shows a significant difference of gas entrainment behavior, where the entrainment obtained for = 2E-4 is around 1.7 times that obtained for = 1E-4. This shows that this coefficient can control the entrainment rate. The average cg and dg entrainment values show that dg entrainment is much higher than cg entrainment for both simulations (see Figure 10d). The result of the simulation using = 2E-4, as shown in Figure 10d, is the same as in Figure 8e, so is in the best agreement with the experiment. The velocity vectors represent the tangential projection of the liquid velocity on horizontal planes H1, H2, and H3, as presented in Figure 11. The background color represents the liquid vorticity magnitude. The influence of on the magnitude of the rotation represented by the vorticity magnitude seems to be minor. For all the evaluation planes, no significant difference can be observed from the simulation when using two different .

The Influence of the Computational Cell Size
The interfacial deformation for different mesh sizes is presented in Figure 12. The deformation and gas core length are significantly influenced by the cell size. For the same (i.e., = 1), the gas core length in the simulation using the finer mesh (i.e., Mesh B) is about 70% larger than the one using Mesh A (see Table 3). The differences in the deformation and gas core length affect the gas entrainment rate. The gas entrainment rate for Mesh B is three times that of Mesh A. As discussed earlier, the coefficient of the curvature correction, has an important influence on the shape of the free-surface vortex. For this reason, this parameter is modified. A similar vortex shape and gas core length are obtained for Mesh A with = 1 and Mesh B with = 0.85. As shown in Table  3, a similar entrainment rate is obtained for this case. This means that the entrainment model itself has the ability to be mesh-independent. The observed mesh dependency mainly results from the different curvature correction required for different mesh sizes.   Table 3 also shows the large dependence of the calculated vorticity magnitude on the mesh size. The dependence is caused by the inability of the curvature correction method to predict the same value of the vorticity magnitude for different cell sizes. As a consequence, the gas entrainment also differs largely when a different cell size is used. The gas entrainment model is expected to provide a mesh-independent solution when the turbulence model is able to estimate the mesh-independent magnitude vorticity.

Conclusions
Turbulence modeling is important for capturing the development of a free-surface vortex and gas entrainment. In the CFD simulation using the SST turbulence model, curvature correction is required to increase the sensitivity on the curvature or rotation. A lower value of leads to only a shallow dimple, while higher provides a better gas core formation. The bubble entrainment increases with the increase in , which can be attributed to the stronger rotation captured by the simulation. The results obtained by using a higher are characterized by the higher vorticity magnitude. The increase of vorticity magnitude closer to the connection pipe is observed for all CFD simulations presented in this paper.
A new entrainment model considering the physics of the phenomena was implemented in the GENTOP concept. In general, a bubble entrainment rate of the same order as in the experiment can be achieved. In the case presented in this work, the CFD simulation where the entrainment model was used to control the entrainment model from the potentially continuous gas into the smallest bubble class group gave the highest gas entrainment and the best fit with the experiment. Using a higher diameter of entrained bubbles led to a smaller gas entrainment. The entrainment coefficient can be adjusted to control the entrainment rate without significantly altering the local velocity fields or the vorticity magnitude. A higher coefficient leads to higher gas entrainment. The gas entrainment rate also significantly affects the gas entrainment rate, which mainly results from the different curvature correction required for different mesh sizes.