Numerical Simulation Study on Distribution of Bubble in Flow near Aerator Based on CFD-PBM Coupled Model in Tunnel

The aerator can reduce erosion by mixing a large amount of air into the water in the solid wall area. The eﬀectiveness 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 diﬃcult area of research. Aiming at the downstream aerated ﬂow of hydraulic aeration facilities, several numerical models, such as VOF, mixture, Euler, and Population Balance Model (PBM), are compared and veriﬁed 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 ﬂow 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 ﬂoat up along the aerated ﬂow, which provides a powerful support for the basic research on the mechanism of aeration and erosion reduction.


Introduction
e 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. e 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 highspeed 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. e phenomenon of high-speed flow aeration in the discharge structure belongs to the problem of the gas-liquid two-phase flow. e 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. e 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. e 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.
e 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 scrubbingcooling chamber of the entrained-flow coal gasifier using a three-dimensional Eulerian-Lagrangian model. e collisions between particles were taken into account by means of the Direct Simulation Monte Carlo (DSMC) method based on the hard-sphere model. e 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. e 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. erefore, in the simulation of highspeed 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. e 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. e 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. e present study aims to compare the simulation accuracy of different numerical models for the high-speed aerated flow. e coupling of PBM and CFD methods is used to further predict the bubble distribution in the flow field downstream of the aerator. e mechanism of bubble aggregation and breakage has been considered. In the calculation of the gas-liquid twophase 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. is 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]. e continuity equation is as follows:

Numerical Model and Validation
e momentum equation is as follows: where t is the time, p is the pressure shared by all phases, ρ q is the physical density of the phase q, v → q is the velocity of the phase q, and τ q is the qth phase stress-strain tensor.
where α q is the volume fraction, μ q and λ q are the shear and bulk viscosity of the phase q, F → q is an external body force, is a virtual mass force, and F → td,q is a turbulent dispersion force.

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]. e momentum equation is as follows: where v → is the velocity component and ρ and μ are the average density and dynamic viscosity, respectively. e interfacial force source term F → is equivalent to the pressure increment caused by the interfacial tension at the interface. e volume fraction α q is used to describe the fraction of different materials in each cell; α q varies in [0, 1] and should sum to unity everywhere: α q is calculated by the VOF equation: e 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. e mixture model solves the momentum equation of mixture and describes the discrete phase by relative velocity [23]. e continuity equation for the mixture is as follows: where v m �→ is the mass-averaged velocity: and ρ m is the mixture density: e 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, F → is a body force, and μ m is the viscosity of the mixture: and v → dr,q is the drift velocity for the secondary phase q:

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.
e interaction between bubbles and the influence of turbulence vortices make bubbles coalesce and break, and bubbles exist in multiscale in the flow field. erefore, in addition to the conservation of momentum, mass, and energy, it is necessary to add an equilibrium equation to describe the bubble equation [24]. e expression is as follows: where n i is the bubble number density, u i is the bubble velocity, B B,i represents the source term of the bubble i generated by the breakage of all bubbles larger than the bubble diameter d i , D B,i represents the vanishing term caused by the disappearance of the bubble I, B C,i represents the source term of the bubble i generated by the aggregation of all bubbles smaller than the bubble diameter d i , and D C,i Mathematical Problems in Engineering represents the vanishing term caused by the aggregation of bubbles within or with other groups of bubbles. e bubble aggregation rate model: where V i and V j represent volume of bubbles i and j, is the bubble collision frequency, and P ag (V i , V j ) is the bubble aggregation efficiency. e aggregation probability model is as follows: where ζ ij � d i /d j , c is the added mass coefficient, consider 0.5 in this article, and We ij 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 P b (f v |d, λ) 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 f v , b(d) is the bubble breakage rate with size d, b(f v , d) is the rate density function of the bubble with size d and breakage ratio f v , ϖ λ (d) is the collision frequency of bubbles and turbulent vortices, λ is the size of the turbulent vortex, u λ is the average velocity of the turbulent vortex, and n λ is the number density of turbulent vortices with size λ. e Luo breakage model [25] is as follows: where Ω br (V, V ′ ) is the bubble breakage rate, V is the original bubble volume, V ′ is the sub-bubble volume, 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. e 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. e 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. e results show that the calculation results of the ASMM considering different bubble rise velocities are close to those considering a single bubble rise velocity.
is suggests that the assumption of using the same local phase velocity for different bubble groups is reasonable. e net source term of all discrete elements due to aggregation and breakage S bi 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 d 32 is defined as follows: where N i is the number of bubbles in the group I and L i 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. e 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:

Turbulence Model
Due to the complexity of turbulence, the commonly used turbulence models are simplified based on various assumptions. e assumption of the standard k − ε method is that the flow is fully turbulent, and the effects of molecular viscosity are negligible. e turbulent viscosity μ t 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 μ t is the turbulent viscosity which can be deduced for the turbulence intensity k and energy dissipation rate ε, 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. e turbulence assumption of isotropy is made for the turbulent flow field, which can better describe the complex flow of the channel flow [27].

Reynolds Stress Model (RSM).
e general form of the Reynolds stress equation is given as follows: Of the various terms in these exact equations, C ij , D L,ij , P ij , and F ij do not require any modeling. However, D T,ij , G ij ,ϕ ij , and ε ij need to be modeled to close the equations.
where u i and u j are time-average velocity components, u i ′ u j ′ is the Reynolds stress component, and σ k � 0.82. e production terms due to buoyancy are modeled as follows: where Pr t is the turbulent Prandtl number for energy, with a default value of 0.85. e classical approach to modeling ϕ ij uses the following decomposition: where ϕ ij,1 is the slow pressure-strain term, also known as the return-to-isotropy term, ϕ ij,2 is called the rapid pressurestrain term, and ϕ ij,w is the wall reflection term.
where Y M � 2ρεM 2 t is an additional dilatation dissipation term. e turbulent Mach number in this term is defined as follows:

Mathematical Problems in Engineering 5
where a is the speed of sound. e scalar dissipation rate ε is computed with a model transport equation similar to that used in the standard k − ε model.

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 F D is the drag force, C D is the drag coefficient of the bubble, u g is the gas phase velocity, and u l is the liquid velocity. e expression of Schiller-Naumann drag coefficient is as follows: where Re is the relative Reynolds number.

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. e three-dimensional model of the experiments is shown in Figure 1. e 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. e 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. e water discharges of exp1/exp2/exp3 are 4400 m 3 /s, 3300 m 3 /s, and 1900 m 3 /s, respectively, corresponding to the three water levels of high, medium, and low. e 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. e water level data are corrected according to the measurement data of flow, velocity, and air concentration. e air concentration is measured by the CQ6-2005 air concentration instrument, as shown in Figure 2.
e 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. e cavity shape under the medium water level of Q � 3300 m 3 /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.

Results and Discussion
e CFD software Fluent is used as the platform to load the numerical model and boundary conditions and realize the CFD-PBM coupling. e VOF model, mixture model, and Euler model are adopted as the multiphase model. e standard k − εturbulent model and RSM turbulent model are adopted as the turbulent model. e Euler multiphase model, k − εturbulent model, and PBM are used in the CFD-PBM coupling model. e numerical calculation conditions are shown in Table 2. d a in the table is the characteristic diameter of the air bubble. e 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 userdefined 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 e 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. e total number of grids  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   Mathematical Problems in Engineering 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. e 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). e 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. e 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. e layout of grid 2 is shown in Figure 4(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. e bubbles were divided into six bins with an equal diameter ratio of 1 : 2, as shown in Table 3. e proportion of initial bubble groups obeys lognormal distribution, and the Sauter mean diameter d 32 � 0.01 m, as shown in Figure 5.    8 Mathematical Problems in Engineering the net length of the cavity L and the depth of the rollers h r 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. e water-air interface whose VOF value is equal to 0.3 is extracted. e height difference between the lower edge of the water nappe and the rollers' interface of the cavity Δh ≤ 0.5 m is considered as the junction point. e vertical distance from this point to the aerator is the aerated cavity length, as shown in Figure 6. e 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. e 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. e 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. us, 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. ere is no exact characteristic diameter to perfectly simulate the cavity shape and air concentration. e optimal bubble size with the best simulation accuracy and the closest to the experimental data is different. d a � 0.005 m, d a � 0.01 m, and d a � 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. e relative error of the aerated cavity length between the simulation results of each model and the experimental value is shown in Figure 10. e formula is as follows: 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.  Mathematical Problems in Engineering e 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. e formula is as follows: where n is the total number of all data points, y i is the simulated results of each model, and u i is the model experimental value. e 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%.

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. e 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). e 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).
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. e 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.

Results' Analysis.
For the selection of turbulence models, the standard k − ε model has good adaptability and fast calculation speed, while the RSM equation is complex and difficult to converge. us, the standard k − ε 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. e simulation results of the mixture and the Euler models are related to the characteristic diameter of bubbles. e results show that the Euler model has better accuracy and adaptability to the drag model.   erefore, six different bubble characteristic diameters are calculated, respectively, by using the Euler model. When d a 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 d a 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. e 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. erefore, 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.
e Sauter mean diameter of bubbles at different water levels is close to each other. e Sauter mean diameter d 32 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. e 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 threedimensional flow characteristics. e 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. e horseshoeshaped reflux pattern is shown in Figure 17(a). is phenomenon is also observed in the numerical simulation results. e 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

Computational Model Layout.
e CFD-PBM coupled model is of high value for applications in the gas-liquid twophase flow in hydraulic engineering, and the accuracy and reliability were validated through the physical model experiments of Bai et al. [7,28]. e experiments were conducted in a 0.25 m wide rectangular chute with a bottom slope of 14°. e 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. e height of the offset h s is 0.045 m. e inlet velocity is 6 m/s and 9 m/s, and the Froude number is 4.9/7.4. e 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. e bubble size is divided into ten bins with a diameter ratio of 1 : 2, as shown in Table 4. e initial proportion of each bin obeys lognormal distribution, as shown in Figure 19. e Sauter mean diameter d 32 � 0.01 m. e 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. e measurement points at the same position were also used in the numerical simulation which is shown in Figure 20. e 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. e total number of bubbles per unit volume is as follows: where n(x, ϕ, t) 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.

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. e 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.  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. erefore, the distance from the position of the maximum aerated concentration in each section to the bottom gradually increases, as shown in Figure 24.

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. e 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. e 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. e bubble size distribution of the four monitoring points obtained by CFD-PBM coupled model numerical simulation is shown in Figure 25. e slope of the fitting curve is − 1.436/ − 1.609/− 1.532/− 2.116, which is in good agreement with the experimental value. e 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.
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. erefore, in the equilibrium zone, the bubbles below 10 mm increase sharply, while the bubbles above 30 mm almost disappear. e contour map of bubble diameter distribution is shown     in Figure 26. e 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. e 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. e 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. is 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. e contour maps of each bubble component distribution are shown in Figure 27.
Two other points were taken from 3 cm before and after the four monitoring points. e Sauter mean diameter at 12 points along the way decreases, and the slope tends to be gentle. e results of polynomial fitting are as follows: y � − 0.01355x 3 + 0.06779x 2 − 0.011265x + 0.07258. Finally, the Sauter mean diameter is close to 0.01 m, as shown in Figure 28.

Conclusions
e 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. e bubble size distribution and its evolution along the high-speed aerated flow were focused on. e 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. e simulation error for the length of the aerated cavity is small, but the simulation effect for the air concentration is poor. e 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. e 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) rough 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. e 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. e 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. e 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.
e 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. e 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. ese are the future research directions.

Symbols a:
Speed of sound (m/s) b(d): Bubble breakage rate (s − 1 ) b(f v , d): Rate density function of bubble with size d and breakage ratio f v (s − 1 ) B B,i : Source term of bubble generated by the breakage B C,i : Source term of bubble generated by the aggregation C: Air concentration C b : Bottom air concentration C D : Drag coefficient of the bubble Greek letters α: Volume fraction β: e power law dead mass coefficient ρ: Density (kg/m 3 ) μ: Dynamic viscosity (Pa·s) μ t : Turbulent viscosity (Pa·s) λ: Bulk viscosity (Pa·s); size of the turbulent vortex (m) ε: Energy dissipation rate (m 2 ·s − 3 ) τ: Stress-strain tensor ϕ ij,1 : Slow pressure-strain term ϕ ij,2 : Rapid pressure-strain term ϕ ij,w : Wall-reflection term ϖ λ (d): Collision frequency of bubbles (s − 1 ) e 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 (jiamingleiscu@163.com).

Conflicts of Interest
e authors declare that they have no conflicts of interest regarding the publication of this paper.