Abstract
The aerator can reduce erosion by mixing a large amount of air into the water in the solid wall area. The effectiveness of erosion reduction is mainly based on air concentration and its bubble size distribution. However, simultaneous simulation of the air concentration and its bubble size distribution in numerical simulations is still a hot and difficult area of research. Aiming at the downstream aerated flow of hydraulic aeration facilities, several numerical models, such as VOF, mixture, Euler, and Population Balance Model (PBM), are compared and verified by experiments. The results show that the CFD-PBM coupled model performs well compared to other conventional multiphase models. It can not only obtain the evolution law of the bubble distribution downstream of the aerator but also accurately simulate the recombination and evolution process of bubble aggregation and breakage. The Sauter mean diameter of the air bubbles in the aerated flow decreases along the way and eventually reaches a stable value. The bubble breakage is the main process in the development of the bubbles. It reveals the aeration law that the small air bubbles are closer to the bottom plate, while the large bubbles float up along the aerated flow, which provides a powerful support for the basic research on the mechanism of aeration and erosion reduction.
1. Introduction
The cavitation erosion induced by high-speed flow is very prominent in hydraulic discharge structures with a high head. With the development of engineering technology, it is an effective method to prevent wall cavitation erosion damage by setting aeration facilities to promote water flow aeration [1, 2]. Both scientific research and engineering practice have proved that a small amount of air entrainment into the water can significantly reduce or even eliminate cavitation erosion in some cases [3]. However, opinions differ as to the internal mechanism of how air bubbles eliminate cavitation erosion. The research on the effect of bubble size on cavitation erosion is still making new progress. In Cheng et al.’s [4] study, PIV technology was used to obtain the flow field information of the gas-liquid two-phase flow in an aeration tank. Brujan et al. [5] found that the interaction between air bubbles and cavitation bubbles was helpful to prevent erosion damage by high-speed photography. Wu et al. [6] proposed that the bubble size has a significant impact on the reduction of cavitation erosion through model experiments research. Small bubbles support to alleviate cavitation erosion, even at the same air concentration. Bai et al. [7] measured and analyzed the development of bubble characteristics downstream of the spillway aeration facility by using a conductivity probe. It was found that the chord length of bubbles decreased sharply in the impact zone. However, the model experimental work is very heavy, and the information obtained is not detailed enough to fully reflect the overall distribution of the bubble size in high-speed aerated flow.
The phenomenon of high-speed flow aeration in the discharge structure belongs to the problem of the gas-liquid two-phase flow. The main characteristics of the flow field are as follows: (1) one phase is continuous and the other is discrete; (2) there are interphase mixing, diffusion, deformation, and relative slip; and (3) there is the exchange of momentum, energy, and mass between phases. The flow aeration is also affected by the boundary condition, so the influence factors are very complex [8]. With the development of computer technology, many new numerical simulation methods have been applied to the study of the two-phase flow. Compared with the model experiments, numerical simulation has the characteristics of short operation periods and less resource consumption. The most important thing is that it can obtain a lot of comprehensive and detailed flow information of the gas-liquid two-phase flow which cannot be measured or difficult to measure in model experiments.
The current mainstream numerical simulation methods for two-phase flows are divided into the Eulerian–Lagrangian model and the Eulerian-Eulerian model. Representative models based on the Lagrangian framework are the DEM and the Direct Simulation Monte Carlo (DSMC). In a study by Peng et al. [9], the CFD-DEM model has been applied to investigate the underlying mechanisms of segregation and mixing of binary systems of solids. Specifically, the simulation reproduced the main features of the gas-solid flow as captured in experiments. Wu et al. [10] investigated the particle distribution in the scrubbing-cooling chamber of the entrained-flow coal gasifier using a three-dimensional Eulerian–Lagrangian model. The collisions between particles were taken into account by means of the Direct Simulation Monte Carlo (DSMC) method based on the hard-sphere model. The results indicate that the axial distribution of the particle number concentration becomes wave-shaped in the pool of the scrubbing-cooling chamber. Big particles are easy to subside, and small particles are easy to suspend. The DSMC method was also applied by Fan et al. [11] to the study of acoustic agglomeration for understanding wave conditions, and the model is shown to be capable of accurately predicting the dynamic acoustic agglomeration process in terms of the detailed evolution of particle size and spatial distribution. Although the form of the equations in the Eulerian–Lagrangian model is simpler, the pursuit of particle positions greatly increases the computational burden. Therefore, in the simulation of high-speed aerated flows at such large scales, with high air concentrations and a large number of grids in hydraulic engineering, the Eulerian-Eulerian model is more commonly used.
Li et al. [12] used the mixture multiphase model and RNG turbulence model to study the air concentration of the bottom plate and sidewalls after the whole cross-section aerator. In Wang and Sun’s [13] study, simulations were performed for upward air-water bubbly flows in a pipe by employing Fluent’s two-fluid model together with an interfacial area transport equation (IATE) model. It is proposed that the lift force is essential to obtain accurate lateral phase distribution. Bak et al. [14] proposed a semitheoretical correlation developed from a steady-state bubble number density transport equation for predicting the distribution of local bubble size using a computational fluid dynamics (CFD) code. However, there are still few studies on bubble size distribution in high-speed aerated flow in hydraulic engineering.
As the traditional Eulerian-Eulerian model is unable to predict statistical information such as bubble collision, aggregation, and breakage in the aerated flow, the CFD-PBM coupled model was born. On the basis of the Eulerian-Eulerian model, many scholars coupled the population balance equation to predict particle size, mass transfer, and other information. The Population Balance Model (PBM) is widely used in petroleum, chemical industry, ship dynamics, and other fields owing to its ability to calculate particles of different sizes in the gas-liquid two-phase flow. Fukuma et al. and Saxena et al. [15, 16] found that the bubble diameter increases with the increase of the superficial gas velocity and finally reaches a maximum value. Li et al. [17] found the same phenomenon and proposed that large bubbles concentrated in the center of the tower, while small bubbles concentrated in the sidewall. The Sauter mean diameter of bubbles increases with the increase of the superficial velocity and decreases from the center of the bubble column to the sidewall.
The present study aims to compare the simulation accuracy of different numerical models for the high-speed aerated flow. The coupling of PBM and CFD methods is used to further predict the bubble distribution in the flow field downstream of the aerator. The mechanism of bubble aggregation and breakage has been considered.
2. Numerical Model and Validation
2.1. Multiphase Model
2.1.1. Euler Model
In the calculation of the gas-liquid two-phase flow, the Euler model takes water and air as the continuous phase filling the calculation area, and the gas phase is composed of spherical bubbles with uniform size, so it is also called the two-fluid model. This model was first proposed by Ishii [18] and then analyzed and discussed by scholars [19]. Each phase in the model has a set of momentum equations and continuity equations. A key aspect of the application of this two-phase flow model is to provide an appropriate closure relationship for the interface exchange term in the equilibrium equation [20].
The continuity equation is as follows:
The momentum equation is as follows:where t is the time, p is the pressure shared by all phases, is the physical density of the phase q, is the velocity of the phase q, and is the qth phase stress-strain tensor.where is the volume fraction, and are the shear and bulk viscosity of the phase q, is an external body force, is a lift force, is a wall lubrication force, is a virtual mass force, and is a turbulent dispersion force.
2.1.2. VOF Model
In the VOF model, the two fluids share a set of momentum equations, and the volume fraction of each fluid in the calculation domain is tracked on each calculation unit [21, 22].
The momentum equation is as follows:where is the velocity component and and are the average density and dynamic viscosity, respectively. The interfacial force source term is equivalent to the pressure increment caused by the interfacial tension at the interface.
The volume fraction is used to describe the fraction of different materials in each cell; varies in [0, 1] and should sum to unity everywhere:
is calculated by the VOF equation:
The density used in momentum equation (4) is calculated from the volume average of all the materials:
2.1.3. Mixture Model
As with the two-fluid model, all phases are treated as a continuous medium penetrating each other in the mixture model. The mixture model solves the momentum equation of mixture and describes the discrete phase by relative velocity [23].
The continuity equation for the mixture is as follows:where is the mass-averaged velocity:and is the mixture density:
The momentum equation for the mixture can be obtained by summing the individual momentum equations for all phases. It can be expressed as follows:where n is the number of the phase, is a body force, and is the viscosity of the mixture:and is the drift velocity for the secondary phase q:
2.1.4. Population Balance Model (PBM)
For the simulation of the gas-liquid two-phase flow after the aerator, the correction of gas-liquid interaction force and turbulent energy has an important impact on the fluid motion. With the increase of flow velocity, the turbulence becomes more severe. The interaction between bubbles and the influence of turbulence vortices make bubbles coalesce and break, and bubbles exist in multiscale in the flow field. Therefore, in addition to the conservation of momentum, mass, and energy, it is necessary to add an equilibrium equation to describe the bubble equation [24]. The expression is as follows:where ni is the bubble number density, ui is the bubble velocity, BB,i represents the source term of the bubble i generated by the breakage of all bubbles larger than the bubble diameter di, DB,i represents the vanishing term caused by the disappearance of the bubble I, BC,i represents the source term of the bubble i generated by the aggregation of all bubbles smaller than the bubble diameter di, and DC,i represents the vanishing term caused by the aggregation of bubbles within or with other groups of bubbles.
The bubble aggregation rate model:where Vi and Vj represent volume of bubbles i and j, is the bubble aggregation rate, is the bubble collision frequency, and is the bubble aggregation efficiency.
The aggregation probability model is as follows:where , is the added mass coefficient, consider 0.5 in this article, and is the Weber number.
According to Prince and Blanch and Luo and Svendsen et al., it is believed that only when the turbulent energy of the turbulent vortex in the gas-liquid two-phase flow is greater than the increment of surface energy after bubble breakage, the bubble breakage rate is defined as the product of bubble breakage probability and collision frequency [25]:where is the probability density of the breakup when the bubble of size d and the turbulent vortex of size collide with each other at the breakage ratio of fv, b(d) is the bubble breakage rate with size d, is the rate density function of the bubble with size d and breakage ratio , is the collision frequency of bubbles and turbulent vortices, is the size of the turbulent vortex, is the average velocity of the turbulent vortex, and is the number density of turbulent vortices with size .
The Luo breakage model [25] is as follows:where is the bubble breakage rate, V is the original bubble volume, is the sub-bubble volume, , , , , and .
A homogeneous discrete method was used for bubble size grouping. In the homogeneous discrete method, the bubble size range of the bubble group is discretized into a finite size interval. The advantage of this method is that the bubble size distribution can be calculated directly. Before the solution, the bubble size distribution can be roughly estimated. When the numerical fluctuation is between 2 and 3 times, the homogeneous discrete method is very effective.
The transport equation of the phase fraction of the discrete element is expressed as follows:
In the homogeneous discrete method, since all discrete elements belong to the secondary phase, all discrete cells call the same phase velocity. Chen et al. [26] used the Algebraic Slip Mixture Model (ASMM) to test this hypothesis. The results show that the calculation results of the ASMM considering different bubble rise velocities are close to those considering a single bubble rise velocity. This suggests that the assumption of using the same local phase velocity for different bubble groups is reasonable. The net source term of all discrete elements due to aggregation and breakage Sbi is zero, which can be expressed as follows:
In order to combine the PBM with the whole hydrodynamics problem, the Sauter mean diameter can be used to represent the bubble diameter of the air phase. For the discrete method, the Sauter mean diameter d32 is defined as follows:where Ni is the number of bubbles in the group I and Li is the bubble size of the group i.
When the PBM is set up, the probability density of the bubble size of each component in the air phase at the inlet of velocity can be given by the lognormal distribution function.
The lognormal distribution for the number density n, as a function of the bubble size L, can be written as follows:where and are, respectively, the location and scale parameters of the distribution and can be written as follows:
2.2. Turbulence Model
2.2.1. Standard Model
Due to the complexity of turbulence, the commonly used turbulence models are simplified based on various assumptions. The assumption of the standard method is that the flow is fully turbulent, and the effects of molecular viscosity are negligible. The turbulent viscosity is computed by combining k and as follows:
(turbulent kinetic energy) equation is given as follows:
(dissipation rate of turbulent kinetic energy) equation is given as follows:where is the turbulent viscosity which can be deduced for the turbulence intensity k and energy dissipation rate , , , , , and .
In this model, the turbulent kinetic energy equation represents the characteristic velocity, and the turbulent energy dissipation equation represents the characteristic length scale. Both the characteristic length and the characteristic velocity are solved by a partial differential equation. The turbulence assumption of isotropy is made for the turbulent flow field, which can better describe the complex flow of the channel flow [27].
2.2.2. Reynolds Stress Model (RSM)
The general form of the Reynolds stress equation is given as follows:
Of the various terms in these exact equations, Cij, DL,ij, Pij, and Fij do not require any modeling. However, DT,ij, Gij,, and need to be modeled to close the equations.where ui and uj are time-average velocity components, is the Reynolds stress component, and .
The production terms due to buoyancy are modeled as follows:where is the turbulent Prandtl number for energy, with a default value of 0.85.
The classical approach to modeling uses the following decomposition:where is the slow pressure-strain term, also known as the return-to-isotropy term, is called the rapid pressure-strain term, and is the wall reflection term.where is an additional dilatation dissipation term. The turbulent Mach number in this term is defined as follows:where a is the speed of sound. The scalar dissipation rate is computed with a model transport equation similar to that used in the standard model.
2.3. Drag Model
Drag force is the most important force of momentum transfer between gas and liquid phases, which mainly represents the blocking effect of the surrounding liquid on bubbles. Due to the problem of spherical particles, the Schiller–Naumann drag coefficient model is used in this paper:where FD is the drag force, CD is the drag coefficient of the bubble, is the gas phase velocity, and ul is the liquid velocity.
The expression of Schiller–Naumann drag coefficient is as follows:where Re is the relative Reynolds number.
2.4. Model Validation
In order to verify the accuracy and reliability of the numerical simulation, the physical model experiments were carried out at the State Key Laboratory of Hydraulics and Mountain River Engineering, Sichuan University, China. The three-dimensional model of the experiments is shown in Figure 1. The width of the spillway in the experimental study section is 15 m, and the gradient of a straight section of the upstream gentle slope is 0.015 and that of a straight section of the downstream steep slope is 0.5. The upstream and downstream are connected by a circular arc with radius R = 100m. At the end of the circular arc, the bottom plate descends to form a step aerator. The water discharges of exp1/exp2/exp3 are 4400 m3/s, 3300 m3/s, and 1900 m3/s, respectively, corresponding to the three water levels of high, medium, and low. The scale of the normal model is 1 : 40 based on the Froude number similarity. In order to facilitate observation, the model is made of transparent plexiglass, and the roughness is about 0.0079–0.0083, which meets the requirements of resistance similarity. Due to the severe turbulence, the steel ruler is used to measure the average water level in the aerated area. The water level data are corrected according to the measurement data of flow, velocity, and air concentration. The air concentration is measured by the CQ6-2005 air concentration instrument, as shown in Figure 2.
The results of the experiments are shown in Table 1. Rollers appeared in the aerated cavity under the three experimental conditions but did not reach the elevation of the vent. The cavity shape under the medium water level of Q = 3300 m3/s is shown in Figure 3. With the increase of the discharge, the Froude number decreases, the flip distance of the lower edge of the water nappe increases, the rollers’ length of the bottom cavity increases, and the net cavity length decreases slightly.
(a)
(b)
3. Results and Discussion
The CFD software Fluent is used as the platform to load the numerical model and boundary conditions and realize the CFD-PBM coupling. The VOF model, mixture model, and Euler model are adopted as the multiphase model. The standard turbulent model and RSM turbulent model are adopted as the turbulent model. The Euler multiphase model, turbulent model, and PBM are used in the CFD-PBM coupling model. The numerical calculation conditions are shown in Table 2. da in the table is the characteristic diameter of the air bubble.
The boundary conditions are as follows:(1)Inlet boundary: the inlet is divided into two parts, which are the velocity-inlet boundary for the water phase and the pressure-inlet boundary for the air phase, respectively(2)Outlet boundary: pressure outlet, utilizing the user-defined function file to control the downstream water depth(3)Wall boundary: no-slip velocity boundary condition; the near-wall regions of the flow were analyzed using the method of standard wall function
The mesh of the numerical model is configured as a structural grid. Sensitivity analyses were carried out for three different grid sizes, and the quality of each group of grids was ensured to be above 0.95. The total number of grids corresponding to Grid1/Grid2/Grid3 was 70311/139802/661700, respectively. The pressure and velocity results of different grid computing are shown in Figures 4(a) and 4(b). From the analysis of the simulation results of the hydrodynamic parameters, it can be seen that grid 1 is too sparse and the accuracy is not enough, while grid 2 and grid 3 give relatively reasonable and consistent results. Pressure on the bottom plate increases sharply at the point of impact of the water nappe and then decreases. However, grid 1 is the densest grid, which takes up more resources and takes longer to compute. To ensure both computational accuracy and efficiency, the grid used in this article is grid 2, which has a more appropriate number of grids. The calculated results for grid 2 were compared with the experimental data and subjected to error analysis, and the results are shown in Figure 4(c). The calculated results are in good agreement with the experimental data, with a maximum error of 1.94% for pressure distribution at X = 69m. Grid 2 has a minimum grid size of 0.2 m and a maximum grid size of 1 m. The grid is locally encrypted to 0.2 m near the aerator to ensure sufficient accuracy to capture the cavity rollers’ phenomenon and the water-air interface, and denser grids are used throughout the bottom plate area to ensure the accuracy of the simulation of the air concentration near the bottom plate. The layout of grid 2 is shown in Figure 4(d).
(a)
(b)
(c)
(d)
In order to simulate bubbles of different sizes, it is necessary to group the bubble sizes and give the initial bubble size distribution in the gas-liquid system at the inlet interface. The bubbles were divided into six bins with an equal diameter ratio of 1 : 2, as shown in Table 3. The proportion of initial bubble groups obeys lognormal distribution, and the Sauter mean diameter d32 = 0.01 m, as shown in Figure 5.
3.1. Length of the Aerated Cavity
Considering that there were certain rollers in the cavity under various conditions, the net length of the cavity L and the depth of the rollers hr are taken as the indicators to measure the cavity shape. In order to avoid the randomness of rollers’ length, the average value within 5 s of numerical results is used to analyze. The water-air interface whose VOF value is equal to 0.3 is extracted. The height difference between the lower edge of the water nappe and the rollers’ interface of the cavity is considered as the junction point. The vertical distance from this point to the aerator is the aerated cavity length, as shown in Figure 6. Four sections X = 35 m/40 m/45 m/50 m on the central axis downstream of the aerator are selected to compare the cavity shape and aerated cavity length of the VOF, mixture, Euler, and CFD-PBM coupled model. The typical sections are shown in Figure 7.
Contour maps of the aerated cavity are shown in Figure 8. The result of the Euler- shows that the aerated cavity shape is narrow, and there is a large amount of rollers with relatively high air concentration in the aerated cavity. The length and depth of the aerated cavity in the VOF model and the Euler- are similar to the actual situation relatively, but the upper edge curve of rollers is similar to a circle arc, which is slightly different from the physical model experiments. The aerated cavity shape of the CFD-PBM is closest to the physical model experiments.
(a)
(b)
(c)
(d)
The variation charts of the air concentration along the water depth are shown in Figure 9, and the value of the air concentration greater than 0.9 is considered as an aerated cavity. The lower edge of the air concentration curve can reflect the shape of the aerated cavity. It can be seen in Figure 9(a) that the cavity shape calculated by the VOF-RSM is slender, which is quite different from the actual situation. Compared with the other models, the lower edge of the air concentration curve calculated by the VOF model changes more sharply, and the bottom air concentration value is far from the measured value. The simulation results of the mixture model in Figure 9(b) show that the air concentration at the rollers of the cavity is higher than the experimental value, and the transition of the air concentration curve of the mixture model is smooth. Thus, it can be seen that the cavity shape is narrow. According to Figures 9(c)–9(e), there is a great difference in the simulation of the cavity shape by choosing different characteristic diameters in the Euler model. And, the performances of each characteristic diameter are also different under different water levels. There is no exact characteristic diameter to perfectly simulate the cavity shape and air concentration. The optimal bubble size with the best simulation accuracy and the closest to the experimental data is different. da = 0.005 m, da = 0.01 m, and da = 0.05 m are the optimal bubble diameters at high, medium, and low water levels. Figure 9(f) shows the simulation results of the CFD-PBM coupled model at three water levels. It can be seen that the model can adapt to various water levels. With the variation of inflow parameters such as water level and discharge, the aerated concentration curves are relatively consistent, which are in good agreement with the experimental values.
(a)
(b)
(c)
(d)
(e)
(f)
The relative error of the aerated cavity length between the simulation results of each model and the experimental value is shown in Figure 10. The formula is as follows:
(a)
(b)
(c)
Obviously, the relative errors of the CFD-PBM coupled model are extremely small at all three water levels, which explain the great advantages of the PBM model for the simulation of the aerated cavity length.
The relative root mean square error (RRMSE) is used to compare the errors of the air concentration near the bottom plate between the several simulated results and the experimental value. The formula is as follows:where n is the total number of all data points, yi is the simulated results of each model, and ui is the model experimental value. The RRMSE values of different models are shown in Figure 11, which show that the CFD-PBM coupled model has a minimal error value of 43.2%, while the VOF model has the largest error value of 83.3%.
3.2. Bottom Air Concentration
Due to a distinct interface between the phases obtained in the VOF model, water and air cannot be evenly mixed. When the jet flows out of the step aerator and impinges the bottom plate, some air will be mixed into the water in the form of an air mass or airbag due to the broken water-air interface. The distribution on the bottom plate is shown in Figure 12(a). However, the mixture and Euler models show a decreasing trend in the bottom air concentration, as shown in Figures 12(b) and 12(c). The air concentration calculated by the CFD-PBM coupled model shows a more symmetrical decay trend, with a high air concentration in the bottom plate at the impact zone of the water nappe, as shown in Figure 12(d).
(a)
(b)
(c)
(d)
Since the variation trend of the aeration concentration is symmetrical on both sides of the axis, the air concentration on the centerline of the bottom plate downstream of the aerator at the medium water level is taken, as shown in Figure 13. The air concentration in the aerated cavity is 1, and the air concentration in clean water is 0. Almost all of the air concentrations in the bottom plate in the VOF model simulations are zero. All air concentrations in the mixture model results are high and tend to decrease slowly. With the decrease of the characteristic diameter selected by the Euler model, the air concentration along the bottom plate increases gradually. It is exciting that the result of the CFD-PBM coupled model is very accurate. Both the air concentration value and reduction trend agreed well with the experimental value.
3.3. Results’ Analysis
For the selection of turbulence models, the standard model has good adaptability and fast calculation speed, while the RSM equation is complex and difficult to converge. Thus, the standard turbulence model is used in the subsequent calculation of the Euler model.
For different multiphase flow models, the VOF model can simulate the aerated cavity shape which is close to the experimental phenomenon, but the simulation accuracy for the air concentration is poor. The simulation results of the mixture and the Euler models are related to the characteristic diameter of bubbles. The results show that the Euler model has better accuracy and adaptability to the drag model. Therefore, six different bubble characteristic diameters are calculated, respectively, by using the Euler model. When da is small, the air concentration along the bottom plate is higher, and the cavity shape is narrow; with the increase of the characteristic diameter, the air concentration at the cavity decreases gradually and the cavity shape expands.
For different da calculated by the Euler model, the air phase mainly distributes near the bottom plate for the bubble with small uniform bubble size, while the bubble with larger uniform bubble size tends to escape upward along the way. The larger the bubble size, the more pronounced is its upward escape, resulting in a lower air concentration in the bottom plate. For different inflow conditions, it is difficult to select an optimal diameter for different Froude numbers. Therefore, although the Euler model can also simulate the real aerated cavity shape, it has higher requirements for the characteristic diameter of bubbles. For different parameters of the flow, the numerical simulation results of the Euler model with uniform bubble size are not accurate and inconsistent with the actual situation.
From the CFD-PBM coupled model simulation results with different bubble size bins, the results under each water level are excellent, which are in good agreement with the experimental value, and the distribution law of the air concentration is quite close to the actual situation. Furthermore, it can be seen from Figure 14 that the CFD-PBM coupled model can obtain more information, including the distribution of the bubble diameter in the water downstream of the aerator. Large bubbles are distributed in the middle and surface of the water phase, while small bubbles are more distributed near the bottom plate. After bubble aggregation and breakage, the distribution of bubble size at the outlet section is shown in Figure 15.
(a)
(b)
The Sauter mean diameter of bubbles at different water levels is close to each other. The Sauter mean diameter d32 of the outlet section is 0.0114 for the high water level, 0.0120 for the medium water level, and 0.0123 for the low water level. The regularities of bubble size distribution are similar at different water levels, as shown in Figure 16.
According to the results of the physical model, it can be clearly observed that the rollers in the cavity exhibited three-dimensional flow characteristics. The central part of the cavity rollers was squeezed to the sides and flowed downstream under gravity from the highest point and forms a swirling flow with downstream rollers. The horseshoe-shaped reflux pattern is shown in Figure 17(a). This phenomenon is also observed in the numerical simulation results. The CFD-PBM coupled model simulation results show that the rollers form a gyration on both sides with the central axis as the boundary, and the streamline disperses at the falling point of the water nappe behind the aerator, as shown in Figures 17(b) and 17(c).
(a)
(b)
(c)
4. Model Application
4.1. Computational Model Layout
The CFD-PBM coupled model is of high value for applications in the gas-liquid two-phase flow in hydraulic engineering, and the accuracy and reliability were validated through the physical model experiments of Bai et al. [7, 28]. The experiments were conducted in a 0.25 m wide rectangular chute with a bottom slope of 14°. The measured section is 5m long, including a smooth convergent nozzle with a length of 1.8 m, a width of 0.25 m, and a height of 0.15 m. The height of the offset hs is 0.045 m. The inlet velocity is 6 m/s and 9 m/s, and the Froude number is 4.9/7.4. The layout of numerical simulations is shown in Figure 18, where the model grid is shown in Figure 18(b), and the minimum grid size is 0.005 m.
(a)
(b)
The bubble size is divided into ten bins with a diameter ratio of 1 : 2, as shown in Table 4. The initial proportion of each bin obeys lognormal distribution, as shown in Figure 19. The Sauter mean diameter d32 = 0.01 m. The distribution of bubble chord length at four monitoring points in the impact zone (X/L = 1.06 and X/L = 1.15) and equilibrium zone (X/L = 1.21 and X/L = 1.64) along the centerline of the bottom plate was measured. The measurement points at the same position were also used in the numerical simulation which is shown in Figure 20. The height of the monitoring points was 3 mm from the floor, and the statistical time was 40 seconds which is the same as the physical model experiments. The total number of bubbles per unit volume is as follows:where is the number density function, where x represents the external coordinates of bubbles, that is, the spatial position of bubbles, and represents the internal coordinates of bubbles, which includes the size, composition, and other variables of bubbles.
4.2. Results’ Analysis
4.2.1. Air Concentration
In the cavity zone, the aeration phenomenon only exists at the bottom and top of the flow, and the air concentration gradually decays to 0 from the lower surface up to the zone of black water. When X/L > 1, the distribution of the air concentration gradually becomes more uniform as the distance increases, which is due to the bubbles entering from the bottom cavity gradually floating upwards in the water to mix evenly. The air concentration increases vertically along the bottom and reaches a peak, after which it gradually decays. A comparison of the section air concentrations in the model experiments and numerical simulations is shown in Figure 21. The numerical simulation results are in good agreement with the model experiment results.
The results of the experimental and numerical simulation for the bottom air concentration Cb at the bottom of the chute (Z = 1.5 mm) on the central axis are shown in Figure 22. The air concentration decays slowly along the axial axis in the equilibrium zone, varying from 0.17 at X/L = 1.2 to 0.08 at X/L = 2.0; the numerical simulation data agrees well with the model experimental value.
Figure 23 shows the water phase in the axial plane. After the impact zone, the air concentration of the bottom plate decays along the way, and the large bubbles float up under the influence of buoyancy, drag force, and lift force. Therefore, the distance from the position of the maximum aerated concentration in each section to the bottom gradually increases, as shown in Figure 24.
4.2.2. Bottom Bubble Diameter Distributions
According to the model experiment results, the size of bubbles with dense distribution is between 1 mm and 10 mm, and the diameter of the conductivity probe used in model experiments is 0.25 mm, so the data collected for bubbles less than 1 mm may be less than the actual situation. The distributions of bottom bubble chord length were skewed with a large proportion of small bubbles. Compared with the impact zone, the distribution range of bubble chord length in the equilibrium zone is reduced. The curve of distribution of the bottom bubble chord sizes was also sharper and narrower along the chute. For bubbles larger than 1 mm in the impact and equilibrium zones, the power law between the bubble number density and the diameter is fitted. In the impact zone, showed a slight increase from −1.2 at X/L = 1.06 to −5/3 at X/L = 1.06. In the equilibrium zone, showed a slight increase from −5/3 at X/L = 1.21 to −2.4 at X/L = 1.64. The bubble size distribution of the four monitoring points obtained by CFD-PBM coupled model numerical simulation is shown in Figure 25. The slope of the fitting curve is −1.436/−1.609/−1.532/−2.116, which is in good agreement with the experimental value. The trend in bubble size distribution is similar to the experimental values, while the number of small bubbles at each point in the numerical model is slightly higher than the experimental values. As the number of bubbles pierced by the probe is counted at each point in the experiment, small bubbles may not be pierced at the tip of the needle due to the tip effect, resulting in missing statistics. But the numerical simulation counts all bubble sizes passing through this grid.
(a)
(b)
(c)
(d)
From the distribution of each bubble bins, the bubble size is larger when the air is introduced into the water in the cavity zone, and the bubbles in bin 0 account for a larger proportion. After the water nappe impacts the bottom plate, the large bubbles break up into smaller bubbles due to turbulent vortex collisions, viscous shear forces, and surface instability of the large bubbles. Therefore, in the equilibrium zone, the bubbles below 10 mm increase sharply, while the bubbles above 30 mm almost disappear. The contour map of bubble diameter distribution is shown in Figure 26. The buoyancy, drag force, lift force, turbulent diffusion force, and virtual mass force are different for bubbles of different sizes so that small and large bubbles will show different kinematic behaviors. The large size bubbles have significantly higher upward velocity than the small size bubbles, while the small bubbles are still distributed near the bottom plate by moving around the wall. The small bubbles of bins 6/7/8/9 smaller than 1 mm are all distributed near the bottom plate behind the equilibrium zone, which is an important guarantee for maintaining a certain air concentration near the bottom plate to achieve the purpose of aeration and erosion reduction. This phenomenon indicates that small bubbles contribute more to the erosion reduction effect in the equilibrium area and the far zone, and the large bubbles will play a smaller role along the way. The contour maps of each bubble component distribution are shown in Figure 27.
(a)
(b)
(c)
(d)
(e)
(f)
Two other points were taken from 3 cm before and after the four monitoring points. The Sauter mean diameter at 12 points along the way decreases, and the slope tends to be gentle. The results of polynomial fitting are as follows: . Finally, the Sauter mean diameter is close to 0.01 m, as shown in Figure 28.
5. Conclusions
The complexity of the high-speed aerated flow is reflected in the changes in bubble size distribution due to bubble breakage and aggregation under severe turbulence, and the experimental experience is not fully sufficient to capture the entire flow field information. In order to investigate the air concentration and bubble size distribution in the high-speed aerated flow downstream of the aerator, several different numerical models including the CFD-PBM coupled model are compared. The bubble size distribution and its evolution along the high-speed aerated flow were focused on. The conclusions are as follows:(1)For the gas-liquid two-phase flow numerical simulation of the high-speed aerated flow, the VOF model can simulate the shape of the aerated cavity. The simulation error for the length of the aerated cavity is small, but the simulation effect for the air concentration is poor. The maximum relative root mean square error can reach 83.3%.(2)In the Euler model, different bubble characteristic diameters have a great influence on the simulation results; the smaller the bubble characteristic diameter, the more serious the air entrainment with higher air concentration. For different Froude numbers, different characteristic diameters should be selected to ensure accuracy.(3)In addition to accurately predicting the air concentration, the CFD-PBM coupled model can also accurately predict the bubble size distribution in the entire flow field, providing more comprehensive information on the flow field than the model experiments. The results of the CFD-PBM coupled model are accurate for different inflow conditions, and the regularities of bubble size distribution are similar at different water levels with a similar Sauter mean diameter.(4)Through physical model experiments observation and CFD-PBM coupled model calculation, it is found that the flow characteristics at the rollers of the aerated cavity are three-dimensional flow. After the water nappe impacts the bottom plate, a pair of cyclonic flows are formed along the central axis of the floor, forming a dynamic balance.(5)Downstream of the aerator, the air concentration near the bottom plate decays along the course after the impact zone, and the position of the maximum air concentration in each section shifts slightly upwards. The reason is that the buoyancy, drag force, lift force, turbulent diffusion force, and virtual mass force vary for different bubble sizes so that different kinematic behaviors occur for small and large bubbles. Bubbles larger than 10 mm float at a significantly higher rate than small bubbles smaller than 10 mm. The larger diameter bubbles will escape upward, while the small bubbles are still distributed near the bottom plate by moving around the wall, which is an important guarantee for maintaining a certain air concentration near the bottom plate to achieve the purpose of aeration and erosion reduction.(6)In the cavity zone, the bubble size is larger when the air is introduced into the water. After the water nappe impacts the bottom plate, the large bubbles break up into smaller bubbles due to turbulent vortex collisions, viscous shear forces, and surface instability of the large bubbles. In the impact zone and the equilibrium zone, the number of bubbles versus bubble diameter varies as a power law of −5/3. Due to bubble breakage, the number of bubbles exceeding 10 mm after the equilibrium zone is greatly reduced, and the Sauter mean diameter gradually decreases along the axial direction and eventually tends to stabilize. The farther away from the impact zone, the larger the proportion of small bubbles, and after the equilibrium zone X/L = 1.64, the Sauter mean diameter tends to be 0.01 m.
The development of the multiphase model is still imperfect. In the current simulation, most of the nucleation, growth, breakage, and aggregation models used in the PBM model are developed in the low-velocity multiphase flow, such as the chemical industry. The parameter selection of each model in the application of the high-speed aerated flow in hydraulic engineering is worth further discussion. Furthermore, the bubbles in the high-speed flow may have different geometry shapes under different inflow conditions, and different geometry shapes may have a great influence on the interphase force. These are the future research directions.
Nomenclature
Symbols| a: | Speed of sound (m/s) |
| b(d): | Bubble breakage rate (s−1) |
| : | Rate density function of bubble with size d and breakage ratio (s−1) |
| BB,i: | Source term of bubble generated by the breakage |
| BC,i: | Source term of bubble generated by the aggregation |
| C: | Air concentration |
| Cb: | Bottom air concentration |
| CD: | Drag coefficient of the bubble |
| da: | Air bubble diameter (m) |
| d32: | Sauter mean diameter (m) |
| DB,i: | Vanishing term caused by the disappearance of bubbles |
| DC,i: | Vanishing term caused by the aggregation of bubbles |
| : | Breakage ratio of the bubble |
| : | External body force (N) |
| FD: | Drag force (N) |
| : | Lift force (N) |
| : | Wall lubrication force (N) |
| : | Virtual mass force (N) |
| : | Turbulent dispersion force (N) |
| : | Acceleration of gravity (m·s−2) |
| h0: | Approach flow depth |
| hr: | Depth of the rollers (m) |
| L: | Net length of the aerated cavity (m) |
| Lexp: | Experimental value of aerated cavity length (m) |
| Mt: | Turbulent Mach number |
| n: | Number |
| ni: | Bubble number density |
| p: | Pressure (Pa) |
| Prt: | Turbulent Prandtl number |
| Re: | Reynolds number |
| Sbi: | Net source term of all discrete elements due to aggregation and breakage |
| t: | Time (s) |
| ui: | Bubble velocity (m/s) |
| : | Velocity (m/s) |
| : | Drift velocity (m/s) |
| V: | Volume (m3) |
| : | Sub-bubble volume (m3) |
| : | Weber number |
| : | Volume fraction |
| : | The power law dead mass coefficient |
| : | Density (kg/m3) |
| : | Dynamic viscosity (Pa·s) |
| : | Turbulent viscosity (Pa·s) |
| : | Bulk viscosity (Pa·s); size of the turbulent vortex (m) |
| : | Energy dissipation rate (m2·s−3) |
| : | Stress-strain tensor |
| : | Slow pressure-strain term |
| : | Rapid pressure-strain term |
| : | Wall-reflection term |
| : | Collision frequency of bubbles (s−1) |
| a: | Air |
| b: | Bubble; bottom |
| exp: | Experiment |
| : | Gas |
| i, j: | Axial component; bubble group |
| l: | Liquid |
| m: | Mixture; mass |
| q: | The qth phase |
| r: | Roller. |
Data Availability
Some or all data, models, and codes generated or used during the study are available from Jiaming Lei upon request ([email protected]).
Conflicts of Interest
The authors declare that they have no conflicts of interest regarding the publication of this paper.
Acknowledgments
This work was supported by the National Science Fund for Distinguished Young Scholars (No. 51625901) and the National Natural Science Foundation of China (No. 51879176).