Volume of Fluid Computations of Gas Entrainment and Void Fraction for Plunging Liquid Jets to Aerate Wastewater

: Among various mechanisms for enhancing the interfacial area between gases and liquids, a vertical liquid jet striking a still liquid is considered an e ﬀ ective method. This method has vast industrial and environmental applications, where a signiﬁcant application of this method is to aerate industrial e ﬄ uents and wastewater treatment. Despite the huge interest and experimental and numerical e ﬀ orts made by the academic and scientiﬁc community in this topic, there is still a need of further study to realize improved theoretical and computational schemes to narrow the gap between the measured and the computed entrained air. The present study is a numerical attempt to highlight the air being entrained by water jet when it intrudes into a still water surface in a tank by the application of a Volume of Fluid (VOF) scheme. The VOF scheme, along with a piecewise linear interface construction (PLIC) algorithm, is useful to follow the interface of the air and water bubbly plume and thus can provide an estimate of the volume fraction for the gas and the liquid. Dimensionless scaling derived from the Fronde number and Reynolds number along with geometric similarities due to the liquid jet’s length and nozzle diameter have been incorporated to validate the experimental data on air entrainment, penetration and void fraction. The VOF simulations for void fraction and air-water mixing and air jet’s penetration into the water were found more comparable to the measured values than those obtained using empirical and Euler-Euler methods. Although, small overestimates of air entrainment rate compared to the experiments have been found, however, VOF was found e ﬀ ective in reducing the gap between measurements and simulations.


Introduction
The significance of the liquid plunging jet's entrainment of a gas lies in achieving a large interfacial area between gas and liquid, which is useful for a large number of environmental and industrial applications. Fuel consumption optimization of ships is compromised due to the formation of a wake on the water surface associated to air being entrained across the boundary layer [1,2]. Pressurized water reactors (PWRs) face partially filled cold legs associated to emergency core cooling as a result of gas entrainment when the coolant strikes the water's surface [3,4]. Plunging jets have been proved efficient methods to aerate industrial effluents [5][6][7][8]. Other applications include mineral-processing flotation cells [9,10] and oxygenation of cell bioreactors [11]. In another industrial process involving the casting of polymers and glass [12], promoting enhanced interaction between gas and liquid aeration is significant for the efficient formation of products.
Among the several applications of the plunging jet mentioned above, aeration and floatation of industrial effluents and wastewater treatment [8,13] represent a huge application of this method. For aerobic treatment processes, such as the activated sludge process, plunging jet aeration systems provide a simple and inexpensive method of supplying oxygen for wastewater treatment [13,14]. main emphasis is to determine the distribution of air void fraction at various depths of the plume from the water-free surface and its corresponding impact on the values of the air entrainment rate.

Physical Parameters
The schematic illustration in Figure 1 shows a vertical plunging water jet of diameter D j , which intrudes into the still surface of a water tank at a velocity u j . The relevant dimensionless numbers, that describe the physics of the phenomenon are the Reynolds number (Re), Froude number (Fr), and Weber number (We). Reynolds number can be expressed as: where, u j is the jet's velocity when impacting the still surface of the water pool, D j is the jet's diameter at the instance when the jet impacts the still surface of the water, ρ w is the water density, i.e., 1000 kg/m 3 and µ w is the water's viscosity, i.e., 1 × 10 −3 pa.s. Froude number is expressed as: where, g is the gravitational acceleration, i.e., 9.8 m 2 /s. Weber number is written as: where, σ is the surface tension corresponding to the air-water interface, 0.072 Nm −1 .
ChemEngineering 2020, 4, x FOR PEER REVIEW 3 of 19 with PLIC proved effective to characterize the distinction between the gas and liquid phase, which can prove useful in the present case to determine the bubbly plume and void fraction along the plume. However, our main emphasis is to determine the distribution of air void fraction at various depths of the plume from the water-free surface and its corresponding impact on the values of the air entrainment rate.

Physical Parameters
The schematic illustration in Figure 1 shows a vertical plunging water jet of diameter Dj, which intrudes into the still surface of a water tank at a velocity uj. The relevant dimensionless numbers, that describe the physics of the phenomenon are the Reynolds number (Re), Froude number (Fr), and Weber number (We). Reynolds number can be expressed as: where, is the jet's velocity when impacting the still surface of the water pool, Dj is the jet's diameter at the instance when the jet impacts the still surface of the water, is the water density, i.e., 1000 kg/m 3 and is the water's viscosity, i.e., 1 × 10 −3 pa.s. Froude number is expressed as: where, g is the gravitational acceleration, i.e., 9.8 m 2 /s. Weber number is written as: where, is the surface tension corresponding to the air-water interface, 0.072 Nm −1 .

Numerical Scheme
The computational simulation has been conducted utilizing the commercial code Fluent applying VOF to characterize the multiphase flow. The VOF model, developed by Hirt and Nicholas

Numerical Scheme
The computational simulation has been conducted utilizing the commercial code Fluent applying VOF to characterize the multiphase flow. The VOF model, developed by Hirt and Nicholas [27], can track the interface between the two phases, so it can be used in all those studies where an interface between two phases exists. A two-equation, k-ε scheme using standard wall values in the vicinity of the wall layer has been applied to describe the turbulence.
The VOF has been applied here to solve multiphase flows. The VOF method deals with flow equations through averaging them in terms of volume to attain a group of equations, whereas the interface between the gas and liquid can be followed by applying a phase marker (α), better named volume function, which can be expressed as: (i) If α = 1 (i.e., computational cell occupies phase 1, e.g., liquid), (ii) If α = 0 (computational cell occupies phase 2, e.g., gas), (iii) 0 < α < 1 (both phases are present and an interface exists within the computational cell). For those cells containing either phase 1 or 2, the formulation consists of mass conservation equations and momentum conservation equations containing properties of the phase present in the cell. If both phases exist within the cell, then the properties of the air and water mixture are evaluated through a mean voidage of the fluids present within the respective cell, thus ensuring to acquire a single set of equations by averaging the variables, for the present case, these variables are density and viscosity, their averaging is expressed as: Using averaging variables and considering that across the interface, the velocity corresponding to the two phases is uniform, the continuity and momentum equations [28,29] expressed as [30,31] ∂ρ ∂t where, F S symbolizes the continuum surface force (CSF) [30,32]. This force acts across the interface, containing the interfacial and subsequent surrounding cells. F S can be expressed as: where κ represents the curve of the surface due to air or water on it. κ can be computed [31,32] as: where, → n is the unit vector representing the surface perpendicular to the interface, however, at the wall, the surface normal to the computational cell, adjacent to the wall will be considered [30,32,33], → n is expressed as: where θ is the angle between the interface and the wall and → n x and → n y stand for the unit vectors along tangential and normal to the wall, respectively. However, this is to confirm that all simulations that have been obtained here, are at a contact angle of 90 • from the wall.
Since the volume fraction, α, is a property of the fluid, which changes along with the fluid, thus, its evolution is dictated by the simple advection relation, expressed as [34]: Due to the fact the VOF model is more flexible and efficient compared to others used to track liquid-gas interfaces, the VOF has become useful in examining a wide range of gas-liquid problems [21,35,36]. Progression of the interface as a function of time between the gas and liquid is tracked in the computational cell using scheme, piecewise linear interface construction (PLIC). As prescribed by Youngs [37], the method considers a linear form of the interface between the two phases within all those cells where the two phases meet each other. Linear behavior of the interface within a cell ( Figure 2) and distributions of the normal and tangential velocity at the face have been used to compute the advection of fluid (Equation (10)) across each face of the cell. Subsequently, the void fraction in every cell (α) is determined owing to the fluxes, computed in the previous cell. the computational cell using scheme, piecewise linear interface construction (PLIC). As prescribed by Youngs [37], the method considers a linear form of the interface between the two phases within all those cells where the two phases meet each other. Linear behavior of the interface within a cell ( Figure  2) and distributions of the normal and tangential velocity at the face have been used to compute the advection of fluid (equation 10) across each face of the cell. Subsequently, the void fraction in every cell (α) is determined owing to the fluxes, computed in the previous cell.

Turbulence Equations
The liquid leaving the nozzle has a velocity high enough to be in a turbulent regime. The gas entrapped due to the liquid jet interacting the surface of the still liquid pool and the gas-liquid two phase flow thus formed, plunge into the pool of liquid at a mixed velocity in the range of turbulent flow. Thus, it is reasonable to compute the gas-liquid jet flow in the liquid pool by application of equations for turbulent quantities acquired by the liquid and gas. The flow variables (e.g., velocity, pressure, etc.), contain mean and fluctuating components. Thus, incorporating the instantaneous equations (Equations (6a,b)) having decomposed variables and utilizing averaging of Reynolds stresses to the Navier-Stokes equations, along with the kinetic energy, k, of turbulence, and its dissipation rate, ε. After dropping the overbar on the mean variables for convenience, The Reynoldsaveraged equation after ignoring the overbar for mean quantities can be expressed in a general form as: where represents the scalar variable, Γ is the diffusion coefficient and is the source term. The representation of these quantities for different variables are given in Table 1. Also, the radial (r) and axial (z) momentum equations along with k and ε equations and values of the respective constants appearing in these equations are also part of this Table.

Turbulence Equations
The liquid leaving the nozzle has a velocity high enough to be in a turbulent regime. The gas entrapped due to the liquid jet interacting the surface of the still liquid pool and the gas-liquid two phase flow thus formed, plunge into the pool of liquid at a mixed velocity in the range of turbulent flow. Thus, it is reasonable to compute the gas-liquid jet flow in the liquid pool by application of equations for turbulent quantities acquired by the liquid and gas. The flow variables (e.g., velocity, pressure, etc.), contain mean and fluctuating components. Thus, incorporating the instantaneous equations (Equation (6a,b)) having decomposed variables and utilizing averaging of Reynolds stresses to the Navier-Stokes equations, along with the kinetic energy, k, of turbulence, and its dissipation rate, ε. After dropping the overbar on the mean variables for convenience, The Reynolds-averaged equation after ignoring the overbar for mean quantities can be expressed in a general form as: where ϕ represents the scalar variable, Γ ϕ is the diffusion coefficient and S ϕ is the source term. The representation of these quantities for different variables are given in Table 1. Also, the radial (r) and axial (z) momentum equations along with k and ε equations and values of the respective constants appearing in these equations are also part of this Table. µ Transport equations relating to momentum are solved to determine turbulence quantities of the two-phase flow by sharing the variables k and ε as well as the Reynolds stresses (i.e., −ρu'u', −ρu'v', u' and v' are fluctuating velocities along with the axial and transverse flow) by the phases across the flow. However, in the governing equations for turbulent production (k) and its dissipation (ε), presented in Table 1, the consequence of turbulence (i.e., Reynolds stresses) is linked to the average flow through the standard k-ε turbulence model [38]. The standard version of the turbulence model has been chosen in our study here due to the sensible computing time and is useful for high Reynolds number flows. Table 1. Source terms and values of coefficients of Generic Transport Equations [37].

Geometry, Mesh, and Boundary Conditions
The 2d meshed domain of the front section of the plunging liquid jet assembly is shown in Figure 3. The meshed domain is used to compute the liquid jet plunging into the pool of still water with subsequent gas entrainment as well as the development of the bubbly plume. The computational domain consists of a 0.15 m radius, where the height varies between 0.397, 0.42 and 0.47 m, which comprises the liquid pool's depth (0.32 m), nozzle pipe (50 mm), and liquid jet's length, with values of 27, 50 and 100 mm. The liquid jet's lengths adopted here are due to the nozzle sizes, i.e., 5, 6.83, 12.5 and 25 mm, which have been chosen here to follow Chanson's [39] experimental setup. The dimensions of the Chanson's experimental tank are sufficiently large so as not to have any impact on the contact of the water jet in the vicinity of the water-free surface. Here, the dimensions of the domain, including sizes of the water jet and the water pool (i.e., r = 0.15 m and h = 0.47 m) are sizeable enough that the boundaries have no influence, however, the jet's domain is much smaller than the pool dimensions. Due to the preference of 2d geometrical scheme and symmetry boundary conditions, the meshed region is reduced to about 25% of the entire geometry, which is useful to reduce the computing time as well as falling within the limitations of the computing grid storage.
The inward flowrate exiting the nozzle located at the top of the water free surface and the out rate of the geometry is located at the bottom, also, wall and outlet pressure constitute the boundary conditions [40]. Water flowing through the nozzle follows a constant mean velocity profile. At the wall, no slip has been assumed. Further, an open boundary condition is considered to model the top of the domain with gauge pressure 0 Pa. However, the outlet flow from the horizontal section of the tank should follow pressure outlet conditions.
Clustering is applied towards the centerline of the liquid jet and the top water free surface to track the phases owing to the liquid and the gas phases, along with the deformation of the interface between the gas and the liquid since the tiniest gas bubble experienced in the experiment involving plunging water jet utilizing tap water is nearly 0.7 mm [39]. Thus, the cell sizes being retained are of asize smaller than the tiniest bubbles within the plume. The fluid domain can be divided into two regions: the region that has to inhibit the jet impingement as well as the interaction regions, appropriately referred here as region of interest (ROI) between the impinged jet and the surrounding water, the mesh size in the region of interest (ROI) may approximately be within 0.01-0.03 mm, whereas in the region that is external to the main flow region (ROI) the mesh size has been keep a little coarse which is around 0.7-0.85 mm. In this way we have tried to reduce the total number of meshes with the best-attained results for the parameters of our concern (i.e., jet impingement length and the jet diameter) with the optimized meshing. The fluid domain can be divided into two regions: the region that has to inhibit the jet impingement as well as the interaction regions, appropriately referred here as region of interest (ROI) between the impinged jet and the surrounding water, the mesh size in the region of interest (ROI) may approximately be within 0.01-0.03 mm, whereas in the region that is external to the main flow region (ROI) the mesh size has been keep a little coarse which is around 0.7-0.85 mm. In this way we have tried to reduce the total number of meshes with the best-attained results for the parameters of Along the region surrounding the centerline, the cells are quadratic having excellent computational prospects owing to skewness factor and orthogonal quality [41]. At all faces of the computational cell, boundary values for physical and flow properties are specified, which include gas and liquid velocity components, acting pressure (which is atmospheric pressure in the present set-up), the level function, the gas number density, the turbulent kinetic energy, k, and the turbulent dissipation, ε. To model the liquid jet at the top surface of the domain, the vertical share of the jet's velocity is set as the nozzle outlet velocity, whereas, all remaining components of the velocity relating to the liquid and the gas as well as the function owing to the bubble's number density are considered as zero. As stated, before that for pressure, a zero-gradient boundary condition has been set, however, the function associated with the level set has been considered such that the periphery of the liquid's jet coincides with the level function having the zero, whereas, its value elsewhere by application of the level function is chosen to specify its zero value that coincides with that of the periphery of the liquid jet, whereas, its value elsewhere is equivalent to the distance from the level setting as zero. The model variables due to the turbulence, k, and ε, relating to their free-stream values as 9 × 10 −4 /Re and 0.9, respectively [21] have been set, which, have been verified to be significantly small so as to not affect the simulations.
For computations relating to fluid volumes, the mass rate at outlet boundaries has been set as less than 10 −6 and the 2-D derivatives were used in the conservation equations are discretized by the use of the quick scheme, while, geometrical re-construction has been applied to achieve discretization of the void fraction in advection Equation (10). The pressure-velocity coupling has been achieved through the algorithm, pressure-inherent with the splitting of operators, PISO [42], whereas, time discretization is used as the inherent first order function applicable to the conservation equations.
An extensive grid independence test was conducted by varying the mesh size with the same Hex type of mesh elements. The varying mesh size with the total number of meshes along with the its effect on the penetration length and jet diameter has been presented in Table 2. Grid independence tests show that upon varying the density of the meshes from 0.31 million to 0.75 million, the highest values of the jet penetration length to the diameter of the jet has been found to be equal to its highest values i.e., (jet penetration length/jet dia 84.58-99.98/19.55-24.98).
The highest values of the jet penetration length/jet dia has been achieved at 0.57 million meshes. After that computations were obtained for higher mesh sizes, however, the values for both quantities only differed by 1.0-0.01/0-0.01. Thus in order to save computational resources and time it was decided to continue with the 0.57 million meshes with the provision that a higher mesh resolution has to be used in the area where the jet impingement occurred all the way till the end of bubbly jet. As the total number of cases ranged to 9 for each mesh density thus we presented the cumulative results of the grid independence in the form of the ranges for the properties of interest i.e., the jet penetration length and the jet dia. in Table 2.

Numerical Results
To achieve similarity between the computations and the experimental results, non-dimensionized equations have been set to achieve right values for the jet's size and the jet velocity that lead to Reynold's numbers of 12,050-105,600 and Froude numbers ranging from 5 till about 9, which correspond to the experimental liquid jet velocity (1.79 m/s-4.4 m/s) analogous to the empirical nozzle dia (6.83, 12.5 and 25 mm). These values were considered for the purpose of comparison with the simulations relating to the void fraction. For evaluation of simulated outcomes for penetration of bubbly plume, experimental values [43] related to nozzle dia (0.3, 1.3, 2.4 mm) and jet velocity, ranging from 6.79-15.7 m/s, were used to evaluate the accuracy of the simulations. All these details including the geometrical dimensions and fluids properties along with the respective values of Re and Fr are summarized in Table 3.
Simulations have been performed initially iterating for a steady liquid phase alone. A steady-state solution has been achieved after time iterations of approximately 2000 with 0.015 dimensionless time, which is equivalent to real-time of nearly 0.1 ms. The steady conditions have been verified by inspecting no change in the time-averaged quantities by raising the number of iterations. Thus, before moving on to the simulations for the bubbly flow, a steady flow of the liquid phase alone has been set as the initial condition.
The simulation involving bubbly flow has been initiated by setting on the multiphase model. The bubbly flow simulation has been obtained with the time step the same as of the single phase, and bubbles travel vertically down across the computational domain towards the region where the bubbles' buoyancy is sufficient to enable the bubbles to move up. However, for simulating gas entrainment and void fraction profiles along with a depth of water pool at 2D j from the contact surface, here air entrainment and void fraction measurements have been achieved experimentally. The bubbles have taken nearly 600 time periods to cover the computational domain equivalent to 2D j . The time-averaged void fractions have been obtained after attaining a steady value after approximately 3000 time periods. The measurements such as gas entrainment and void fraction, made by [39] involving water jet sinking into the water pool have been used for comparison with the simulations.

Gas Entrainment and Developing Bubbly Regime
Sequential illustrations of the simulations presented in Figure 4 show the formation of a gas plume within the liquid pool. VOF simulation exhibits entrainment of the air by the water jet, that has been driven deeper in the tank. The entrained gas bubbles have been advected by the downward flow to form a plume. The depth of the plume associates with the momentum of the water jet, however, bubbles can only be driven down to a depth where the bubble's buoyancy overcomes the momentum owned by the liquid jet, thus, bubbles rise due to the buoyancy to reach the surface of the pool, where they the gas separates the liquid. Similar computational results were found by [18,21,25]. The volumetric downward gas rate at a given level below the surface is given by: where q g is the volumetric rate flux of the gas. Since the 1980s, empirical predictive modeling equations relating to plunging liquid jets have been proposed [15,44,45], however, simple relations based on geometry conditions have been proposed in these studies and absolutely no influence due to the physical properties of the fluids into the predicting modeling has been emphasized in them. Bin [13] was the first who has realized the impact of fluid properties into his empirical predictive modeling: where Q g is the downward vertical gas flux density or gas (air) entrainment rate (m 3 /s), Q j is the liquid volumetric rate (m 3 /s) and h is the depth of the plume from the impact (m). The above equation has been found agreeable with the data following conditions falling within, The present simulated air entrainment ( values along with Chanson [39] experimentally measured values and Bin's [13] correlations as well as Euler-Euler simulated values [25], are summarized in Table 4. The VOF simulations have revealed a closer agreement with air entrainment values determined from Bin's correlation and experiments compared to the Euler-Euler simulations. This shows that VOF scheme is effective to trace the air-water interface.  Moreover, the gas entrainment in the present situation is modeled using the advection (Equation (10)) that provides sufficient protection from the errors in entrainment caused due to the numerical method, whereas, within the region across the developing bubbly plume, the distribution of the bubbles is not influenced by the interfacial forces as only surface tension has been included in the momentum equation. For improvement in the bubble's distribution within the plume's profile across Moreover, the gas entrainment in the present situation is modeled using the advection (Equation (10)) that provides sufficient protection from the errors in entrainment caused due to the numerical method, whereas, within the region across the developing bubbly plume, the distribution of the bubbles is not influenced by the interfacial forces as only surface tension has been included in the momentum equation. For improvement in the bubble's distribution within the plume's profile across the depth of the water pool, there is a need to improve the modeling by enhancing the interfacial interaction between gas and liquid.

Penetration Depth
After nearly four seconds of physical time, the gas plume reaches a steady penetration depth, h p , below the water free surface. The simulated values for the penetration of bubbly plume can be seen in Figure 5, also recent data for the circular plunging jet from Guyot [43] on the penetration of bubbly plume has been included in this figure. It's interesting to express here the empirical correlation proposed by Bin [13] for plume's penetration (h p ): ChemEngineering 2020, 4, x FOR PEER REVIEW 8 of 19 The relation has been used to determine the depth of the gas plume as per the sizes of the configuration. The relation is applicable for liquid jet nozzle sizes that vary from 3.9 mm to 12 mm, whereas, and for a pool depth (h) about 0.50 m or less, with the error of ±20%. VOF simulated values for ℎ higher than both experimental [39] and predictive values of the above equation. This equation can also be expressed alternatively [18] as: The relation has been used to determine the depth of the gas plume as per the sizes of the configuration. The relation is applicable for liquid jet nozzle sizes that vary from 3.9 mm to 12 mm, whereas, and for a pool depth (h) about 0.50 m or less, with the error of ±20%. VOF simulated values for h p higher than both experimental [39] and predictive values of the above equation. This equation can also be expressed alternatively [18] as: The relation reveals that the bubbly plume is independent of the quantity of air being entrained and it only depends on the jet's velocity and the nozzle size. Table 5 provides a summary of the penetration depth of bubbly plume below the impact surface for a jet velocity of 2.54 m/s and nozzle dia of 5 mm. The VOF computed penetration depth that follows the continuous entrainment regime against the operating condition is 0.3 m, which exceeds both the empirical correlation (Equation (14)) and the experimental values. However, with the Euler-Euler scheme, the plume depth doesn't reach as much as the experiment (Equation (14)).

Void Fraction
The instant void fraction values have been averaged over time to obtain a steady mean void fraction over 3000 iterations, which have been averaged again along the circumference. The mean simulated values of void fraction are compared with averaged values of void fraction experimentally measured by Chanson [39], obtained under similar operating conditions as in the simulations described here. The experimental void fractions determined by Chanson [39] have been measured along a horizontal axis that passes through the jet's center. These measurements have been averaged by including results either sides of the jet to curtail the error owing to the measurements. It is known that the flow attains the status of fully developed at 40D J [46], which is much outside the computational domain being used here. Subsequently, the plume's depth such as 0.8D j , 1.2D j , and 2D j from the impact surface, have been chosen for simulated void fraction profiles across the bubbly plume.
It has been described under the above topics that the bubbles have been entrained through the impact between the liquid jet and the receiving pool. Further from the impact surface, the bubbles move both inside and outside radially, as the two-phase mixture travels deep into the water pool. All this can be seen in a sequence of snaps shown in Figure 4 for the operating conditions involving a liquid jet velocity of 3.5 m/s and a nozzle size of 25 mm.
The simulated peaked void fraction along with the depth (0.8D j , 1.2D j , and 2.0D j ) from the impact surface along with the corresponding experimental measured values of void fraction are presented in Figures 6-8 The simulated peaked void fraction along with the depth (0.8Dj, 1.2Dj, and 2.0Dj) from the impact surface along with the corresponding experimental measured values of void fraction are presented in Figures 6-8, where the predicted values have been averaged in time. The computed void fraction profiles are obtained for liquid jet velocities of 3.5, 4.1, and 4.4 m/s, to compare both with experimentally measured as well as other simulated void faction profiles.       Since the simulated bubbly jet is spreading at depth, 0.8Dj is small, however, the agreement with measurements suggests that the VOF model accurately depicts the location of air entrainment owing to their sources. Further, it can be seen in these figures that the maximum value of the void fraction agrees well with the measured value, which confirms the accuracy of the VOF model in computing the quantity of air that has been entrained.
.      Since the simulated bubbly jet is spreading at depth, 0.8D j is small, however, the agreement with measurements suggests that the VOF model accurately depicts the location of air entrainment owing to their sources. Further, it can be seen in these figures that the maximum value of the void fraction agrees well with the measured value, which confirms the accuracy of the VOF model in computing the quantity of air that has been entrained.
The plots further down to lower depths from 0.8D j are presented in Figures 10 and 11, which also indicate good agreement between the experiments and the presented simulations. This indicates that VOF accurately predicts the spreading of bubbly jet and the transport of bubbles which are entrapped by the liquid jet at the impact surface. Transport of bubbles is powered by the liquid velocity along with a balance between turbulent dispersion and drag forces acting on a bubble. Also, the present simulated results make good agreement with the analytical results derived from the simple diffusion scheme adopted by Chanson [39]. Such evidence has been derived by Drew and Passman [47], who have indicated that in circumstances when the dispersed phase momentum under the resulting dominance impact of drag force and turbulent dispersion, the mass conservation should be modeled by a diffusion equation.
Thus, the simulations are obtained at 0.8D j , 1.2D j , and 2.0D j to enable comparison with measured values of the void fraction by Chanson [39] at the same downstream distance from the impact surface. The comparisons between the time-averaged steady simulated values of void fraction and the measured values are presented in Figures 9-11.

Conclusions
CFD simulations of water jet impacting into a still water pool are performed using the VOF formulation through the commercial code FLUENT, whereas tracing of the interface between air and water is conducted by applying a PLIC algorithm. The simulated results are confirmed by comparing them with Chanson [39] experiments through considering uniformity conditions of experimental three jet sizes along with three jet velocities and three jet diameters. The simulations agreed well with experimental data. Focusing VOF simulations on the developing region of the bubbly plume, it has been inferred that the location of peak void fraction in both experiments and the simulations match. In the experiment, this maximum decays exponentially for deeper locations. However, in the present simulations, it remains almost constant as the two phases move with the same velocity. This fact is inherent to the VOF modeling which resolves only one equation. Hence, no slip velocity between the two phases can be modeled. Slight overestimates of air entrainment rate compared to the experiments of Chanson [39] are noticed for all the three cases, however, VOF simulations came close to the experimental values as compared to the empirical correlation and Euler-Euler scheme. This indicates that the entrainment should be regulated by more parameters than just mere geometric similarities. Hence, small scale-models should be investigated keenly to describe the clear picture of the air entrainment.