An improved air entrainment model for stepped spillways
AAn improved air entrainment model for stepped spillways
Silje Kreken Almeland a, ∗ , Timofey Mukha b , Rickard E. Bensow b a Norwegian University of Science and Technology, Department of Civil and Environmental Engineering, NO-7491Trondheim, Norway b Chalmers University of Technology, Department of Mechanics and Maritime Sciences, H¨orsalsv¨agen 7A, SE-412 96Gothenburg, Sweden
Abstract
Numerical modelling of flow in stepped spillways is considered, focusing on a highly economical approachcombining interface capturing with explicit modelling of air entrainment. Simulations are performed onspillways at four different Froude numbers, with flow parameters selected to match available experimentaldata. First, experiments using the model developed by Lopes et al. (Int. J. Nonlin. Sci. Num., 2017) areconducted. An extensive simulation campaign is used to carefully evaluate the predictive accuracy of themodel, the influence of various model parameters, and sensitivity to grid resolution. Results reveal that, atleast for the case of stepped spillways, the number of parameters governing the model can be reduced. Acrucial identified deficiency of the model is its sensitivity to grid resolution. To improve the performance ofthe model in this respect, modifications are proposed for the interface detection algorithm and the trans-port equation for the volume fraction of entrained air. Simulations using the improved model formulationdemonstrate better agreement with reference data for all considered flow conditions. A parameter-free cri-terion for predicting the inception point of air entrainment is also tested. Unfortunately, the accuracy ofthe considered conventional turbulence models proved to be insufficient for the criterion to work reliably.
Keywords:
Air entrainment modelling, numerical modelling, CFD, OpenFOAM, self-aeration, steppedspillway
1. Introduction
Along with a renewed interest in stepped spillways as a flood overflow structure and energy dissipator inhydraulic engineering, attempts at gaining a better physical description of spillway flows have also intensified.A process that is especially challenging to study by means of both physical and numerical experiments, isthe self-aeration of the spillway. Yet, since large quantities of entrained air lead to higher flow depths, releaseenergy, and reduce the potential for damage caused by cavitation, accurate prediction of aeration is crucialfor spillway design. In this work, mathematical modelling and simulation of air entertainment is in focus.To put the present contribution into context, a brief review of the physics of air entrainement in spillwaysis given below, followed by an overview of past attempts of accounting for them in a numerical setting.Generally, air entrainment is driven by turbulent motion and occurs when the turbulent forces at the freesurface overcome the stabilizing effects of surface tension and buoyancy [10]. Applied to spillways, it hassince the early work of Straub and Anderson [26] been widely accepted that the onset of self-aeration takesplace when the turbulent boundary layer, developed from the crest, reaches the free surface. This locationis commonly referred to as the ‘inception point’. Several contributions consider the onset of the aeration indetail [6; 33; 28], and empirical relations exist for the distance to the inception point from the spillway crest ∗ Corresponding author
Email addresses: [email protected] (Silje Kreken Almeland), [email protected] (Timofey Mukha), [email protected] (Rickard E. Bensow)
Preprint submitted to Applied Mathematical Modelling September 11, 2020 a r X i v : . [ phy s i c s . f l u - dyn ] S e p
2; 18; 4]. Boes and Hager [2] proposed a computable definition of the inception point as the location wherethe pseudo-bottom air concentration is 0 . R (cid:13) and has been used in publications on stepped spillways [27; 9]. Lopes et al. [16] incorporated the entrainedair flux estimator developed in [17] into a VoF framework by introducing a separate transport equation forthe entrained air. The solver is implemented in OpenFOAM R (cid:13) , and results from stepped spillway simulationsare presented in the article.In this work, we present further developments of the model by Lopes et al. [16]. To motivate the need forimprovements, results from a simulation campaign, in which the model in its current formulation is appliedto spillway flow at four different Froude numbers, are presented. The simulations constituting the campaignvary in the employed grid resolutions and parameters of the model. Results from the second campaign,in which the model is modified as proposed here, are then presented, demonstrating better accuracy androbustness. The article is supplemented by a dataset containing all the simulation results. The reminder of this article is structured as follows. Section 2 describes the computational methods usedin the performed computations. In Section 3, the setup of the stepped spillway simulations is discussed.Section 4 presents the air entrainment model developed by Lopes et al. [16]. Section 5 contains resultsfrom the simulation campaign in which the model of Lopes et al. [16] is used. Improvements to the airentrainment model are then proposed in Section 6, and corresponding simulation results are provided inSection 7. Concluding remarks are given in Section 8.
2. Computational fluid dynamics methods
The flow was simulated using the Volume of Fluid (VoF) multiphase modelling technique [14], in which asingle set of governing equations is solved for all phases and the location of the interface is identified based onthe values of the cell volume fraction of the liquid phase, α l . Both fluids are considered incompressible and ∂ρu∂t + ∇ · ( ρu ⊗ u ) = −∇ p ρgh − gx ∇ ρ + ∇ (cid:0) µ (cid:0) ∇ u + ( ∇ u ) T (cid:1) − ρu (cid:48) ⊗ u (cid:48) (cid:1) + f s (1) ∇ · u = 0 . (2)Here the overbar denotes the Reynolds average, ρ is the density, µ is the dynamic viscosity, u is the velocityvector, p ρgh = p − ρg · x is the dynamic pressure, and f s is the surface tension force. The latter is approx-imated using the Continuous Force Model, see [3] and also [24] for a detailed discussion in the context ofOpenFOAM R (cid:13) . The term ρu (cid:48) ⊗ u (cid:48) represents the Reynolds stresses, which are to be approximated by theturbulence model.An algebraic approach to account for the evolution of α l is adopted, with the associated transportequation originally formulated as ∂α l ∂t + ∇ · ( uα l ) + ∇ · ( u c (1 − α l ) α l ) = 0 , (3)in OpenFOAM R (cid:13) ’s VoF solvers. The third term in the equation is artificial and is meant to introduceadditional compression of the interface to ensure its sharpness. However, here the formulation of this termis modified in order to accommodate it into the air entertainment modelling framework. The details areprovided in Section 4. The definition of u c is nevertheless not altered: It is aligned with the interface normal,and its magnitude is computed as C α | ¯ u | , where C α is an adjustable constant, here set to 1.Given α , the material properties of the fluids are readily obtained as ρ = α l ρ l + (1 − α l ) ρ air , µ = α l µ l + (1 − α l ) µ air . (4)The indices l and air are used to refer to the water and air, respectively.What remains to be discussed is the choice of turbulence model, which for the case of the stepped spillwayis far from trivial. In principle, the model should be able to properly account for the interaction between theturbulent and multiphase structures in order to provide accurate prediction in the aerated region of the flow.None of the closures that have found widespread use were developed with this goal in mind. Nevertheless,it is common for conventional two-equation models to be used for aerated flows. In [16], which is the workthis article largely builds upon, the k - ω SST model [19] is used for stepped spillway simulations. In [17],it is employed in simulations of a plunging jet and a hydraulic jump. For the latter, many studies also usethe k - (cid:15) model and its variations, a comprehensive list of references can be found in Table 5 in [29]. Qianet al. [23] found the realisable k - (cid:15) model to be favourable for stepped spillway flow. Based on the above, weconsider both the k - ω SST and the realisable k - (cid:15) model [25] and test which of them leads to better predictiveaccuracy.It is pointed out in [11] that in the implementation of the above (and also other) turbulence models inOpenFOAM R (cid:13) the viscous diffusion terms are not treated consistently in regions with a non-zero densitygradient. Since in the simulations of the spillway presented here a density gradient is present across thewhole aerated part of the flow, this issue can have a significant effect on the results. The authors of [11]also provide alternative implementations, in which the inconsistency is resolved. Here, we test using boththe default and the improved implementations. The computations are performed using OpenFOAM R (cid:13) version 5, provided by the OpenFOAM Founda-tion. This CFD tool is based on cell-centered finite volume discretization, which is de facto the industrystandard. Two custom solvers are used in the study, implementing the air entraiment modelling presentedin Sections 4 and 6. Both of them represent modifications of the solver interFoam , which is distributedwith OpenFOAM R (cid:13) . This solver implements the VoF methodology discussed in Section 2.1. The governingequations are solved in segregated manner using a variant of the PISO algorithm [15].3 crucial component of the numerical setup is the selection of the spatial interpolation and time integra-tion schemes. Generally, linear interpolation can be used in space except when considering convective fluxes.In the momentum equation, the latter are interpolated using the limitedLinearV scheme, which is a TVDscheme based on the Sweby limiter. The limiter is computed based on the direction of most rapidly changinggradient and then applied to all three velocity components. This improves stability but at a certain expensein terms of accuracy. For the convection of α in equation (3), a TVD scheme using the SuperBee limiteris used. The van Leer limiter was also considered, but SuperBee led to better results on coarser meshesdue to being more compressive. Unfortunately, in a multi-dimensional setting, using a TVD scheme doesnot guarantee that the values of α will be bounded between 0 and 1. Therefore, OpenFOAM R (cid:13) utilizes anadditional flux limiting technique, referred to as MULES. It is based on the Flux Corrected Transport theorydeveloped Zalasak [32], more details can be found in [8]. The convective fluxes in the turbulence equationsare discretized using the second order upwind scheme, called linearUpwind . This scheme is unbounded,but no significant effect of parasitic oscillations was observed even on coarse grids. Finally, as discussed inSection 4.2 below, the air entrainment model adds an advection-diffusion equation for the flow variable α g ,meant to indicate the distribution of the volume fraction of entrained air , to the system. Here, a TVDscheme using the van Leer limiter is employed.The first-order implicit Euler scheme was used for time-stepping. The choice is not of particular impor-tance here because the flow eventually arrives to an essentially steady state. Nevertheless, a CFL number ≤
3. Simulation cases
This section presents the setup of the stepped spillway simulations used to evaluate the performance ofthe entrainment modelling. In order to have a reference with respect to which the accuracy of simulationresults can be analysed, the flow and spillway parameters are selected to match those in the experimentsof Bung [4]. These were performed on four different spillways combining two selections for the angle ( θ =18.4,26.6 ◦ ) with two for the step height ( s = 0 . , .
06 m). For each spillway, measurements were made for threeflow discharge values ( q = 0 . , . , .
11 m s − ). The parameters θ , s , and q can be used to construct thestep Froude number, which can be considered the main controlling parameter of the flow [18; 5], F s = q (cid:112) g sin θK . (5)Here K = s cos θ is the step induced macro-roughness and g is the acceleration due to gravity. The experi-ments of Bung cover twelve different Froude numbers in the range, 2 . ≤ F s ≤
13. For the simulations fourvalues fairly evenly distributed across this range have been selected: 2 .
7, 4 .
6, 8 .
3, and 13 .
0. The values of θ , s , and q in the four simulation cases are provided in Table 1. This table also provides the values of someauxiliary geometrical parameters, the definition of which can be found in Figure 1. The figure also showsthe origin and orientation of the employed Cartesian coordinate system. Table 1: Setup for the different simulation cases. The number of cells are given in 10 . F s (-) θ ( ◦ ) s (m) q (m s − ) L x (m) L s (m) h win (m) n cells,G n cells,G n cells,G n cells,G igure 1: Sketch of the geometry of the simulation case, identifying the geometric parameters, and also the employed boundaryconditions. on dense grids, no air entrainment is resolved by the VoF. This is shown by the interFoam computationson grid G4. By contrast, in a 3D setting, it is more likely that some interface perturbations eventuallystart getting captured. Investigating the performance of the entrainment modelling in such a scenario is ofinterest, but lies out of scope of the current work.The same boundary conditions were used for all cases, with the exception of the discharge q prescribedat the water inlet. The height of the inlet was adjusted to ensure sub-critical inflow conditions. A zerogradient condition was used for the pressure, while Dirichlet conditions were applied to k , (cid:15) , and ω . Thevalues were set assuming 5% turbulent intensity and 10% of the critical height as the turbulent length scale.For the outlet, a zero gradient condition was prescribed for velocity, pressure, α l and α g . For k , (cid:15) , and ω the OpenFOAM inletOutlet boundary condition was applied. It acts as a zero gradient condition in caseof outflow, but for backflow a homogeneous Dirichlet condition is applied instead.No slip conditions were used at the walls, with a zero gradient condition set for α l and α g . Theturbulent quantities were estimated by regular wall laws, in OpenFOAM named as kqRWallFunction , epsilonWallFunction , omegaWallFunction and nutkWallFunction , for k , (cid:15) , ω , and ν t respectively.For the top boundary the total pressure was fixed, and a pressureInletOutletVelocity condition wasapplied for the velocity. Similar to inletOutlet , this imposes a zero gradient for outflow, whereas forbackflow, it assigns a velocity based on the flux in the patch normal direction. The inletOutlet boundarycondition was used for α l , α g , k , (cid:15) , and ω .The material properties of the fluids were set to correspond to air and water. The values are providedin Table 2.The computational grids were constructed using Pointwise R (cid:13) , and consist of square cells with the excep-tion of a small strip close the top boundary, where unstructured meshing was necessary to account for theslope of the geometry. Four grids with increasing cell density, denoted G1, G2, G3, and G4, were constructedfor each of the four spillways. In each consecutive grid the edge length of the square cells is halved. On thecoarsest grid G1, the edge length is 5 mm, which corresponds to what was used in the simulations by Lopeset al. [16]. This can be related to the the critical height of the spillway flow, defined as h c = (cid:0) q /g (cid:1) / .5 able 2: Simulation parameters. Property Value
Liquid density, ρ Gas density, ρ Liquid kinematic viscosity, ν · − m / sGas kinematic viscosity, ν . · − m / sSurface tension coefficient, σ h c is discretised by either 15 or 21 cells. The numbers forthe G4 grid are, respectively, 126 and 171. The densities are not adjusted to remain equal with respectto h c across all flow conditions, because experiments showed that the relevant parameter for entrainmentmodelling is the resolution of the interface. The number of cells in each mesh is given in Table 1.In conclusion, additional characteristic scales of spillway flow are defined. These will be used for non-dimensionalising the results. At a given x , the height h is defined as the y -coordinate of the point where α air = 0 .
9. The velocity u is defined as the x -component of the mean velocity vector at y = h . Similarscales can be defined with respect to other α air values, e.g. h .
4. Air entrainment modelling
This section presents the air entrainment model developed by Lopes et al. [16]. One can split the modelinto three components: an estimator for the flux of entrained air, a transport equation for the volumefraction of entrained air, and a coupling procedure between the model and the VoF framework. Sections 4.1,4.2, and 4.3 each focus on one of these components. Additionally, for the stepped spillway, estimating thelocation of the inception point is necessary and this constitutes an additional component of the model, whichis treated in Section 4.4.
A key component of the model is the estimation of the quantity of entrained air carried passed someimaginary surface located below the interface. The estimate was proposed by Ma et al. [17]: q = a · Pos ( ∇ ( u · n ) · n ) , (6)where a is a length scale associated with the roughness of the interface due to turbulence,Pos( x ) = (cid:26) x, x ≥ , x < , and n is the interface normal defined as n = ∇ α l / ( |∇ α l | + ε ) . (7)Here ε is a small number added for numerical stability.It is assumed that entrainment is confined to a layer of thickness φ ent close to the surface. Therefore, inorder to obtain a volumetric air entrainment rate, q can be divided by φ ent . Note, however, that (6) is bydefinition not restricted to being non-zero only in the vicinity of the interface. Theoretically, entrainment canbe incorrectly predicted in regions where it should not take place. For this reason, in [16], q is additionallymultiplied by some function δ fs , which is non-zero only close to the interface. The final form of the volumetricair entrainment rate estimate is S g = aφ ent Pos ( ∇ ( u · n ) · n ) δ fs . (8)6t remains to define how a , φ ent , and δ fs are computed. The common approach for a is to equate it tothe turbulent length scale as predicted by the RANS model. The value of φ ent should be related to somecharacteristic length scale of the problem.Within the VoF framework, the α l -field stands out as the natural choice as a basis for the developmentof an interface indicator function such as δ fs . Typically, the interface is defined as the isosurface α l = 0 . α , which can be expected to reach its maximum close to theboundary between the continuous air region and the air-water mixture. H¨ansch et al. [12] used the gradientof α and a function based on tanh to find the interface as part of their air entrainment model. This functionwas adopted by Lopes et al. [16] and reads as δ fs = 12 tanh (cid:2) β ∆ x ( |∇ α l | − |∇ α l | cr ) (cid:3) + 0 . . (9)Here |∇ α l | cr is a constant representing the critical value of the gradient that is expected to be reached in theinterface cells. Its estimate can be computed based on the size of the grid cell, ∆ x : |∇ α l | cr = 1 / (4∆ x ). Theparameter β can be used to control the extent of the interface region with respect to the chosen |∇ α l | cr , andthus provides an opportunity to broaden or restrict the number of cells in which the source term is active. α g -equation The source term (8) is introduced into an additional equation for the modelled volume fraction ofentrained air, α g : ∂α g ∂t + ∇ · ( u g α g ) + ∇ · ( ν t ∇ α g ) = S g . (10)Here ν t denotes the turbulent viscosity. The velocity of the entrained air, u g , is either set to be equal to u or alternatively modified according to [7]: u g = u + u r , (11)where the correction velocity u r is calculated based on a bubble radius according to u r = − r . b g, if 0 < r b ≤ × − m − . g, if 7 × − < r b ≤ . × − m − . r . b g if r b > . × − m . (12)The inclusion of the diffusion term in (10) is considered optional.Additionally, Lopes et al. [16] argue that to properly account for the break-up of bubbles at the freesurface, α g should be set to zero when α air exceeds a certain threshold value, referred to as the BBA. Thesuggested value to use is 0 . α g and its relation to α l are somewhat elusive. In [16], the authors discussthe possibility of using the entrainment model without backward coupling to the VoF solver. In this case,the situation is clear: α g shows the modelled distribution of the volume fraction of entrained air, whichcannot be captured by the VoF. However, when the coupling is two-way (particulars presented below), theidea is that the entrained air should be captured in the α l field, and α g is essentially reduced to a buffer-fieldused to propagate the effect of S g onto α l . Here, we are interested in applying the model in a two-way coupling regime, meaning that the model’spredictions should be propagated into the distribution of α l . The premise is that the VoF simulation byitself does not resolve any entrainment, and therefore all of it is accounted for by a subgrid model based onthe α g equation (10). The overall idea is that α l should be reduced in regions where α g is large, and in amanner that does not disrupt mass conservation. 7ere this is done through a modification of the artificial compression term introduced into the α l -equation (3): ∇ · ( u c (1 − α l ) α l ) . (13)The term (1 − α l ) α l = α air α l is originally meant to serve as an indicator for cells constituting the interface,in which the compression is to be applied. The key observation is that multiplying u rj by some negativenumber instead would lead to interface expansion and thus a region occupied by a mixture. The goal isthen to correlate α g with the change in sign in the term in front of u rj . The most obvious way to do that isexchange α air α l for ( α air − α g ) α l . Note that since (13) is a transport term, mass conservation is guaranteed.The modified α l -equation then reads ∂α l ∂t + ∇ · ( uα l ) + ∇ · ( u c ( α air − α g ) α l ) = 0 . (14)Under the definition above, the model is active only when α g > α air , which is reasonable. It is also worthmentioning that otherwise (13) recovers its original compressive function. This occurs even in the regionsoccupied by a mixture, which can be called into question. As part of the work on improving the model,some experiments have been conducted in which positive values of α air − α g where cut to 0, however theexhibited results were inaccurate, and introducing such a discontinuity is probably best avoided. As discussed in the introduction, surface aeration initiates when the turbulence perturbations exceedthe stabilizing forces of surface tension and buoyancy at the free surface. In the model of Lopes et al. [16]no attempt is made to explicitly compute the force balance. Instead, two model parameters, k c and u c areintroduced, where the subscript c stands for critical. The inception is considered to occur when k > k c and u · n > u c and u · g > u c . (15)Appropriate values for k c and u c are extremely difficult to predict a priori, since the selection clearlydepends not only on the flow conditions (see Section 5.1), but also on the turbulence model and its predictionof k . Careful calibration with respect to the selected model is therefore necessary. In [16], the authors never-theless suggest u c = 0 . k c = 0 . / s , referring to previous experimental results. Unfortunately,how these values relate to the characteristic length and velocity scales of the flow is not discussed.
5. Stepped spillway simulations with the original model
This section demonstrates results from simulations performed using the model of Lopes et al. [16] de-scribed in the previous section. In the original source the model was reported to reproduce experimentaldata on a stepped spillway with F s =2.7. However, the simulations were performed on relatively coarse grids,with the grid sensitivity study performed using only the baseline VoF solver. Furthermore, initial testingwithin this work indicated that its effect was reduced upon grid refinement, which motivates the analysesherein.The implementation of the corresponding solver, called airInterFoam was kindly provided by P. Lopesvia personal communication. Below, we abbreviate airInterFoam to AIF. In the following sections, thesolver is evaluated in terms of sensitivity to flow conditions (Section 5.1), grid resolution (Section 5.2), andseveral parameters of the entrainment model (Section 5.3). F s Here, results from four AIF simulations varying in the prescribed Froude number of the spillway floware presented. To highlight the effect of the entrainment model, data obtained with the baseline VoF solver interFoam (abbreviated IF henceforth) are also provided. The employed numerical parameters are basedon [16], where good results for the F s =2.7 case (simulated in 3D) are presented. In particular, the G1 grid is8mployed, k - ω SST is used for turbulence modelling, the diffusion term is omitted in the α g equation, and k c = 0 . /s , u c = 0, u g = u l , BBA = 0 .
1, and φ ent = 0 . h c .Figure 2 shows the obtained values of α air in the uniform flow region. Good accuracy is achieved for F s =2.7 and 4 .
6, but for the two higher F s the model fails to predict the reduced penetration of air into thecorners of the steps. As a result, in terms of magnitude, the errors in the IF and AIF simulations are similar,although in the case of IF the diffusion of the interface is a purely numerical effect. It is also interesting tonote that for F s =2.7 the accuracy is on par with the 3D simulations using similar model settings conductedin [16]. air [-] y / h [ - ] AIFIFexpr (a) F s =2.7 air [-] y / h [ - ] (b) F s =4.6 air [-] y / h [ - ] (c) F s =8.28 air [-] y / h [ - ] (d) F s =13 Figure 2: Vertical void fraction profiles for uniform flow conditions. Spillway flows with different Froude numbers at the coarsestgrid G1. AIF simualtions compared to IF simulations and experimental results by Bung [4].
The evolution of the surface elevation, measured as h , is shown in Figure 3. The elevation’s value inuniform conditions is well-predicted for all Froude numbers. However, the location of the inception pointsare not captured as consistently. The difference in the obtained values with respect to the experimentaldata of Bung is provided in each plot of the figure: ∆ n i stands for the difference in terms the step number,and ∆ L i in terms of x . It should be noted that when comparing across different F s , using ∆ L i is moreappropriate, since for the two higher Froude numbers, the length of the step is halved relatively to the lowerFroude number cases. The obtained incetion point locations for F s =2.7 and F s =4.6 are reasonably accurate,but, unfortunately, at higher F s the disagreement with the experiment becomes larger. Furthermore, thepredicted inception point for F s =8.3 is further downstream as compared to that for F s =13, whereas theexperimental data exhibits the opposite trend.Figure 3 additionally shows the experimental values of h w , which is the equivalent clear water depth,i.e. the surface elevation that should be predicted by IF. In the obtained results, IF somewhat over-predicts h w , the reason being the coarseness of the grid.The predicted profiles of the streamwise velocity are shown in Figure 4. Remarkably, no effect of airentrainment modelling is visible, and accurate profiles can be predicted with IF. This result was reproducedin all the simulations in this paper, and, for that reason, velocity profiles are not further presented ordiscussed.The principle conclusion from the obtained results is that the settings used in [16] for the F s =2.7 casefail to provide consistently accurate results as the Froude number becomes larger. This indicates that someparameter values of the model, for example k c , should be made a function of F s . To explore the grid sensitivity of AIF, the F s =2.7-case was run on all four grids G1-G4. The resulting α air and h profiles are shown in Figure 5. Clearly, the behaviour obtained on the coarse grid G1 issubstantially changed as the grid gets denser. With increasing resolution, less air is distributed towards the9 n edge y [ m ] L i = 0.27 m n i = 2 AIFIF h expr h w , expr x [m] (a) F s =2.7 n edge y [ m ] L i = 0.19 m n i = 1 x [m] (b) F s =4.6 n edge y [ m ] L i = 0.57 m n i = 6 x [m] (c) F s =8.28 n edge y [ m ] L i = -0.57 m n i = -6 x [m] (d) F s =13 Figure 3: Surface elevation plots, h . Spillway flows with different Froude numbers simulated by AIF on the coarsest gridG1, compared to IF simulations and physical model results by Bung [4]. The difference from the experimentally measuredinception point is annotated in meters (∆ L i ) and in steps (∆ n i ). .0 0.5 1.0 u x / u x , 90 [-] y / h [ - ] AIFIFexpr (a) F s =2.7 u x / u x , 90 [-] y / h [ - ] (b) F s =4.6 u x / u x , 90 [-] y / h [ - ] (c) F s =8.28 u x / u x , 90 [-] y / h [ - ] (d) F s =13 Figure 4: Vertical profiles for the streamwise velocity for uniform flow conditions. AIF and IF simulation results compared toexperimental results by Bung [4]. pseudo-bottom, and for the densest grid only a tiny air layer is found close to the surface. The profiles, both h and α air , approach the corresponding solutions obtained with IF. Obviously, this is caused by the factthat less numerical diffusion contribute to the transport of α g as the grid is refined, but inspection showsthat this is also caused by the shrinkage of the area in which the source term S g is non-zero. This, in turn,is controlled by δ fs , which makes the definition of this function a contributor to the grid sensitivity. A moreelaborate discussion follows in Section 6. On the other hand, with respect to h w , the IF solution consistentlyimproves with grid refinement. On the G4 grid the interface is perfectly sharp and the h w profile is verywell-matched. u r , bubble breakup criterion, and α g -diffusion Here we explore the effects of the entrainment model parameters that could arguably be considerednon-essential or optional. First, the impact of the air bubble drift velocity model, as given in eq. (12), isinvestigated. This is followed by an analysis of the bubble breakup criteria, BBA. Finally, the model istested in terms of activation of the diffusion term in the α g -equation (10). The rest of the numerical setupis similar to that used in the Froude number sensitivity study, see Section 5.1.Here we restrict the analysis to α air profiles in the uniform flow region. For the u r model, three valuesof the bubble radius r b are considered, along with setting u r to 0. The result is shown in Figure 6a.Clearly, including u r leads to reduced air entrainment, and this effect becomes larger when the input bubbleradius is increased. This is expected, since u r is, by definition, directed upwards. Since the intensity of airentrainment is already heavily dependent on the formulation of δ fs , having an additional controller in termsof u r introduces unnecessary complication. At least in the case of spillway flow, setting u r = 0 can thereforebe recommended.Three values of the BBA are considered. The first is 0.1, which is recommended by Lopes et al. [16], andthe other two are 0.05 and 0. Recall that the chosen value refers to the minimal volume fraction of waterfor which α g is allowed to be non-zero. The results from the three simulations are shown in Figure 6b. Asexpected, slightly more air is entrained close to the h when the bba value is reduced. Note that comparedto the experimental data, the air fractions predicted by AIF close to the interface (from y/h ≈ .
8) aretoo low. Thus, deactivating the BBA-criterion produces a small improvement in the accuracy meaning thatthis parameter can be safely removed from the formulation of the model.Finally, the effect of the diffusion term in (10) is investigated. Two simulations, with and without thediffusion term included, have been conducted, see Figure 6c. The effect of the diffusion appears to becompletely negligible. Consequently, it is possible to remove it from the model.11 n edge y [ m ] AIF G AIF G AIF G AIF G IF G IF G IF G IF G h expr h w , expr x [m] (a) Surface elevation plot, h -surface. air [-] y / h [ - ] (b) Step 7 air [-] y / h [ - ] (c) Step 11 air [-] y / h [ - ] (d) Step 15 air [-] y / h [ - ] (e) Step 19 Figure 5: AIF simulations for F s =2.7 at different grids (G1-G4) compared to physical model results by Bung [4]. Figure 5ashows the surface elevation, Figures 5b-5e void fraction profiles for different steps. air [-] y / h [ - ] u r = 0 r b =10 r b =10 r b =10 IFexpr (a) air [-] y / h [ - ] bba=0.1bba=0.05bba=0expr (b) air [-] y / h [ - ] No diffusionWith diffusionexpr (c)
Figure 6: Sensitivity of α air profiles to the slip velocity model, bubble breakup criterion, and α g -diffusion.
12n summary, u r , BBA, and α g -diffusion can be excluded from the model, significantly simplifying itsformulation. β As shown in Section 5.2, at high grid resolutions the effects of the model becomes negligible, or evendeactivated. The most important controller of the breadth of the region (in terms of ∇ α l ) where entrainmentis introduced is the form of δ fs , see Eq. (9). In particular, the parameter β determines the breadth of thetanh function, meaning that with smaller β the region of the model activation becomes larger. In principle,it may thus be possible to maintain an appropriate level of entrainment at high grid resolutions by adjusting β accordingly. However, for this to be possible in practice, the necessary change to β should be easy topredict a priori.Multiple simulations across different Froude numbers and grid resolutions have been conducted in anattempt to determine whether clear guidelines for setting β could be established. Unfortunately, theseefforts were fruitless, and the results of using a given β value change significantly depending on the flow andparameters of the simulation.As an illustrative example, Figure 7 shows the α air profiles produced using β = 10 and 25 in simulationsat different Froude numbers on the G2 grid, and compares it to the corresponding profile produced by AIF,where a β -value of 100 is used as default. For F s =2.7 the sensitivity to β is rather small, but for F s =4.6 asignificant increase in aeration occurs. For larger Froude numbers a lower β mainly leads to more air beingpresent in the corner of the steps. While the change in the α air profiles appears to be rather small, uniformflow conditions are not reached, and its effect continues to grow further downstream. Similar sporadicbehaviour was observed with respect to other simulation parameters. For example, contrary to what isobserved when using the G2 grid, using G3 instead leads to the F s =2.7 case becoming very sensitive β . air [-] y / h [ - ] = 10= 25AIFexpr (a) F s =2.7 air [-] y / h [ - ] = 10= 25AIFexpr (b) F s =4.6 air [-] y / h [ - ] = 10= 25AIFexpr (c) F s =8.28 air [-] y / h [ - ] = 10= 25AIFexpr (d) F s =13 Figure 7: Sensitivity on β for AIF evaluated on grid G2 for different step Froude numbers. Vertical α air profiles are shown.The data is extracted from what should be the start of the uniform flow region according to the experimental data [4].
6. Proposed model developments
In the previous sections, two issues with the entrainment model proposed by Lopes et al. [16] have beenidentified. Perhaps the most critical one is the successive deactivation of the model upon grid refinement.The other one is the difficulty in prescribing the k c value in order to get a good prediction of the inceptionpoint. In this section, improvements to the model are proposed aiming at alleviating these problems. First,an alternative formulation for the δ fs function is introduced in Section 6.1. Afterwards, amplification ofthe diffusion term in the α g is argued for in Section 6.2. Finally, a different way of predicting the inceptionpoint is discussed in Section 6.3. 13 .1. Free surface detection Before discussing the formulation of δ fs , it is necessary to consider what kind of restriction of S g in spaceis needed in the case of spillway flow. A typical distribution of an unrestricted S g is shown in the left plotof Figure 8. The first thick yellow line is located right below the interface and represents the region wherethe entrainment can be expected to take place. However, a discontinuous strip of non-zero values is alsoobserved close to the psuedobottom along with more or less randomly distributed points of activation inthe corners of the steps. Physically, no entrainment can occur in these regions, and the δ fs function shouldfilter them out. Note that the spatial separation between the correct and non-physical regions of sourceterm activation is not large, which explains why defining δ fs in a universal way that fits all flow conditionsand numerical settings is not trivial.To arrive to a better formulation for δ fs , it is important to clearly understand the deficiencies of theoriginal definition, see Eq. (9). With ∇ α l,cr set to 1 / (4∆ x ), the distribution of δ fs over ∇ α l depends on twoquantities: β and ∆ x . It is instructive to see how δ fs changes shape when the values of these parametersare changed. In the left plot of Figure 9, δ fs ( ∇ α l ) is shown for ∆ x values corresponding to grids G1-G3,and the two values of β considered in the sensitivity study in Section 5.4. The tanh function defining thetransition region of δ fs from 0 to 1 is centred at ∇ α l,cr , shown in the figure with black vertical lines. As thegrid is refined, this location is shifted to the right, and for larger values of β , the tanh only spans a limitedrange of high ∇ α l values. On the other hand, for smaller β , the tanh becomes so wide that δ fs remainsnon-zero everywhere.This behaviour of δ fs should be related to how ∇ α l is typically distributed across y , see the right plotin Figure 8. Note that the high values of ∇ α l are always restricted to to a relatively thin region close to theinterface. Consequently, for a large β , for example β = 100 as proposed in [16], δ fs will be restricted to anincreasingly smaller region in space when the grid gets refined. This explains why the diminishing effect ofthe entrainment with grid refinement observed in Section 5.2.As mentioned above, it is easy to make the δ fs function less restrictive by lowering β . However, theissue here is that the β and the effective cut-off value in terms of ∇ α l are not intuitively related. Thisis problematic given that the margin of error is quite small, as discussed above in relation to the typicaldistribution of S g . l [m ] y / h [ - ] G3G2G1
Figure 8: Left: Typical distribution pattern of an unrestricted S g . Right: Typical profiles of the gradient of α l . Alternative options for the expression of δ fs were explored, including a parabola based function, and astep-shaped function defined purely based on the distance from a defined interface. For both alternatives,the idea was to set the value of the critical gradient of α l relatively tight, to capture the upper peak in thegradient plot in Figure 8, and then expand its prevalence away from the these locations according to anappropriate function or logic. In the distance based alternative, δ fs was set to 0 or 1 for a particular cell,depending on its distance from the defined interface. If this distance was less then the interface thickness, φ ent , δ fs was set to 1, and otherwise its value was set to 0. However, this led to step-shaped profiles of α l ,14
100 200 300 400 500 l f s G1, = 25G2, = 25G3, = 25G1, = 10G2, = 10G3, = 10 0 50 100 150 200 250 l Figure 9: Surface indicator functions δ fs . where too much air was entrained below the interface, and indicated a need to apply some functionality toreduce the effects of the source term as the gradient of α l is reduced below its critical value.Acknowledging the above, a parabola based function was considered. The possibility to define ∇ α cr as a top point of the function, and to define a cut-off value for the gradient as an intersection point, wasviewed as beneficial features of this function in the current setting. The latter leading to the possibility ofavoiding the long tail in the tanh-function, with the corresponding uncontrollable potential for generationof none-zero values of S g within the steps. The parabola-based δ fs formulation is proposed as δ fs ( ∇ α l ) = (cid:40) Pos (cid:16) − d (cid:0) |∇ α l | − |∇ α l,cr | (cid:1) + 1 (cid:17) if ∇ α l < |∇ α l,cr | . (16)Here, d refers to the distance from the vertex of the parabola to its focus, which can be computed as d = 0 .
25 ( |∇ α l,cr | − |∇ α l,cut | ) , where ∇ α l,cut is an input parameter explicitly defining the lowest ∇ α l for which the source term mayassume non-zero values. The proposed δ fs , computed for grids G1-G3 and two different values of ∇ α l,cut ,is shown in the right plot of Figure 9. The non-zero values of the function are always fixed to the interval[ ∇ α l,cut , ∇ α l,cr ], which expands upon grid refinement. Unfortunately, due to this expansion, even with thisnew δ fs , the region of non-zero S g values shrinks as the grid is refined. However, the process is slowed, since δ fs is left to be non-zero at lower ∇ α l .During initial testing, it was observed that the parabola-based δ fs was very effective at filtering out thesporadic source term activation in the corner of the steps, while preserving the region where entrainmentis expected. However, the secondary strip of non-zero S g values close to the psuedobottom (see left plot inFigure 8) would sometimes still be left unfiltered. This is likely related to the secondary peak in ∇ α l . Torectify this, the parabolic surface indicator function in Eq. (16) is combined with a distance based approachlike the one outlined above.In details, the values of δ fs computed according to Eq. 16 are additionally manipulated as follows. First,the cells in which δ fs ≥ . S g should be activated. For the remaining cells, the distance to the nearest cellswith δ fs ≥ . φ ent , the value of δ fs in thecell is set to 0. Overall, this combined formulation gave improvements compared to the original formulationbased on tanh. It is noted that even without the distance cut-off modification, improved results with theparabolic δ fs could be achieved, and the δ fs ≥ . ∇ α l,cr and ∇ α l,cut grid independent parameters has briefly been tested, but aban-doned due to the sensitivity of the results to the chosen values. In general, our extensive efforts andexperimenting with δ fs formulations and other simulation settings showed that constructing a robust andaccurate δ fs is very challenging. The sensitivity of the results to any changes in the simulation parametersor flow conditions tends to be very strong. Nevertheless, as shown in the Section 7, the proposed δ fs doesrepresent an improvement with respect to prior art. By definition, the employed entrainment model is meant to account for aeration occurring close to thefree surface, within some layer of thickness φ ent . However, experimental results clearly show that in the caseof the stepped spillway, air penetrates all the way down to the surface of the steps, see e.g. the experimentalprofiles in Figure 2. The physical mechanism through which this occurs is described by Pfister and Hager[22]. Inspection of video recordings from their experiments reveals the occasional generation of air troughsthat extend from the surface into the bulk flow. These troughs penetrate deep enough to hit the step edges,and when they do, the air is distributed into the steps.Capturing this intrinsically transient process in a steady state model is not straightforward. Here, weconsider using a somewhat ad-hoc approach, taking advantage of the diffusion term in the α g equation (10), ∇ · ( ν t ∇ α g ). Recall that in Section 5.3 it was shown that the effect of this term on the solution is essentiallynegligible. However, the effect can be easily amplified by pre-multiplying it with some constant C t . Theincreased diffusion of α g will then lead to air being redistributed more evenly across y , and consequentlyresult in stronger aeration closer to the pseudobottom.The difficulty lies in the choice of the value of C t since there is no clear physical analogy between themodelled phenomenon and diffusion. In light of this, it was attempted to search for a suitable value throughexperimentation, to see whether one leading to improved results across all the considered Froude numberscould be found. As a result, C t = 150 was selected. Even if AIF predicts the aeration onset correctly when appropriate values for k c and u c are supplied,the inception point estimation using this method completely depends on user input. As shown in Section5.1, the appropriate critical value k c depends on the flow conditions. An alternative approach, found in thework of Hirt [13], is to directly consider the balance between the energy of turbulent motion and that ofgravity and surface tension. Defining P t = ρk, (17) P d = ρ | g | a + σa , (18)the source term S g is activated only in cells where P t > P d . This way the inception point prediction requiresno user input. However, it completely relies on the correct prediction of k and a by the turbulence model.
7. Simulations with the improved model
This section is dedicated to evaluating the effects of the model improvements proposed above. To thatend, a new solver incorporating these changes has been implemented. Reflecting the focus on stepped spillwaysimulations, the solver is called spillwayFlow , which is abbreviated to SPF below. The robustness of themodel with respect to grid resolution is evaluated in Section 7.1. Results from application of the model tospillway flow at all four considered Froude numbers are presented Section 7.2. The proposed criterion forinception point location is tested separately, in Section 7.3. https://github.com/siljekre/spillwayFlow .1. Grid sensitivity Here the new model is put to the same grid sensitivity analysis as presented in Section 5.2 for AIF. Tosimplify the analysis, the new inception point prediction approach is not employed, and the original criterionbased on k c is used instead. Simulation results for the F s =2.7 case obtained on grids G1-G4 are shown inFigure 10. Here, C t = 0 and the diffusion term in the α g equation is thus inactive. As anticipated, the resultsstill depend on the grid, and the general trend is convergence towards the IF solution. Note that in the α air profiles, the reduction of aeration manifests itself predominately at y/h (cid:47) .
6. Closer to h results remainacceptable even on the G4 grid. Related to the absence of air in the lower parts, the predictions of h itselfare more sensitive, and unfortunately on the finer grids the accuracy is poor. Nevertheless, compared to theoriginal AIF results (see Figure 5) the robustness of the model is improved. n edge y [ m ] SPF G IF G IF G IF G IF G SPF G SPF G SPF G h expr h w , expr x [m] (a) Surface elevation plot, h -surface. air [-] y / h [ - ] (b) Step 7 air [-] y / h [ - ] (c) Step 11 air [-] y / h [ - ] (d) Step 15 air [-] y / h [ - ] (e) Step 19 Figure 10: SPF simulated with no α g -diffusion ( C t = 0) and F s =2.7 at different grids (G1-G4) compared to physical modelresults by Bung [4]. Figure 10a shows the surface elevation, Figures 10b-10e, the void fraction profiles at different steps. Figure 11 shows the results obtained with C t = 150. As expected, a comparatively more even distributionof α air across y is achieved, in particular for the simulations on denser grids. Furthermore, a very clearimprovement in the robustness of the model is evident, with much more similar results obtained on all fourgrids. Now, the results obtained with the new model are presented for all four considered values of the Froudenumber. The simulations were run on the densest mesh, G4. As in the section above, the new criterion forlocating the inception point is not used here, and instead k c is adjusted for each case in order to match the17 n edge y [ m ] IF G IF G IF G IF G SPF G SPF G SPF G SPF G h expr h w , expr x [m] (a) Surface elevation plot, h -surface. air [-] y / h [ - ] (b) Step 7 air [-] y / h [ - ] (c) Step 11 air [-] y / h [ - ] (d) Step 15 air [-] y / h [ - ] (e) Step 19 Figure 11: SPF simulated with α g -diffusion ( C t = 150) and F s =2.7 at different grids (G1-G4) compared to physical modelresults by Bung [4]. Figure 11a shows the surface elevation, Figures 11b-11e, the void fraction profiles at different steps in thedeveloping region. C t = 0 and C t = 150are shown in all the figures. For comparison, they also include results from AIF simulations on the G1 grid,and also from IF simulations on the G4 grid.The α air profiles are discussed first, see Figure 12. Qualitatively, the same behaviour with respect to C t is observed for all F s . With C t = 0 the distribution of α air across y is close to step-wise, with decentagreement with experimental data for y/h (cid:39) .
7. When C t = 150, the profiles are smoothed out, whichgenerally increases the accuracy. The exception is the F s =8.28 case, for which the air volume fraction at low y/h becomes excessive. Yet even for this case the agreement is better than what could be achieved withAIF. air [-] y / h [ - ] (a) F s =2.7 air [-] y / h [ - ] (b) F s =4.6 air [-] y / h [ - ] (c) F s =8.28 air [-] y / h [ - ] (d) F s =13 Figure 12: Vertical void fraction profiles for uniform flow conditions. Spillway flows at different Froude numbers. SPFsimualtions (G4) compared to IF (G4), AIF (G1) and experimental results by Bung [4]. All plots are showing profiles at stepedges.
The surface elevation plots are shown in Figure 13. Overall, C t = 150 leads to better results, which agreewell with the experimental data in the uniform flow region. An interesting exception is the F s =4.6, for which C t barely has influence on h in the uniform flow region, whereas in the developing region C t = 0 leads tovery good agreement with the experiment. However, it is unlikely that this is explained by any fundamentalproperty of the flow or the model. Compared to AIF, the accuracy of the new model is generally on par.AIF curves are marginally closer to experimental data for the two lower F s , and the other way around forthe two higher F s . In this section, the source term activation criterion presented in Section 6.3 is tested. Recall that theinception point location is determined from the flow, and depends heavily on the employed turbulence model.Here we show results from simulations using four models: the k - ω SST and realisable k - (cid:15) from the standardOpenFOAM library, and their respective counterparts in the library by Fan and Anglart [11], in which thedensity gradient is properly accounted for in the transport equations. The latter are referred to by adding varRho to the name of the model. Simulations using standard turbulence modelling and the k c -criteriapurposed by Lopes et al. [16] are added as reference. The simulations are performed on the G3 grid.The results are summarized in Table 3, and the corresponding distributions of α air are shown in Fig-ure 14. For the two models from the standard library, the inception of entrainment is triggered immediately(downstream the crest) for all the considered Froude numbers. When the varRho variants are employedinstead, the inception point shifts downstream. This reflects the fact that these model predict significantlylower values of k . For k - ω SST the agreement with experimental data is nevertheless poor, but for therealisable k - (cid:15) model the results are more promising.In the results produced using standard turbulence and k c = 0 . / s , the inception point is predictedin relatively good agreement with the experimental results for all cases but the F s =13 case, where the19 n edge y [ m ] AIF G IF G SPF G , C t = 150SPF G , C t = 0 h expr x [m] (a) F s =2.7 n edge y [ m ] x [m] (b) F s =4.6 n edge y [ m ] x [m] (c) F s =8.28 n edge y [ m ] x [m] (d) F s =13 Figure 13: Surface elevation plots, h . Spillway flows at different Froude numbers. SPF simualtions (G4) compared to IF(G4), AIF (G1) and experimental results by Bung [4]. F s =8.3 and F s =13. Compared to the k c =0 . / s criterion, the automatic activation criterion performs on par using the realisable k - (cid:15) turbulencemodel and the variable density turbulence framework.Overall, none of the tested models perform well enough to use the proposed source term activationcriterion. The improved performance of the varRho models shows the importance of properly accountingfor the density gradient in the transport equations for the modelled flow quantities. Improving turbulencemodelling accuracy near interfaces is an active area of research. The results presented here warrant a deeperinvestigation of what models are appropriate for the stepped-spillway flow. Table 3: Inception points found using the source term activation criteria given inEq. (17)-(18) for SPF simulations on grid G3 using different models for turbulencemodelling, compared to physical model results by Bung [4]. F s Grid Turbulence L i,sim L i,expr ∆ L i ∆ n i k - ω SST, k c = 0 . k - (cid:15) , k c = 0 . k - ω SST 0.00 0.67 -0.67 -5.02.7 G3 realisable k - (cid:15) varRho / k - ω SST 1.34 0.67 0.67 5.02.7 G3 varRho /realisable k - (cid:15) k - ω SST, k c = 0 . k - (cid:15) , k c = 0 . k - ω SST 0.00 0.95 -0.95 -5.04.6 G3 realisable k - (cid:15) varRho / k - ω SST 2.00 0.95 1.05 5.54.6 G3 varRho /realisable k - (cid:15) k - ω SST, k c = 0 . k - (cid:15) , k c = 0 . k - ω SST 0.00 1.14 -1.14 -12.08.3 G3 realisable k - (cid:15) varRho / k - ω SST 1.70 1.14 0.56 -5.98.3 G3 varRho /realisable k - (cid:15) k - ω SST, k c = 0 . k - (cid:15) , k c = 0 . k - ω SST 0.00 2.00 -2.00 -21.013 G3 realisable k - (cid:15) varRho / k - ω SST 2.28 2.00 0.28 2.913 G3 varRho /realisable k - (cid:15)
8. Conclusions
This study presents developments in numerical modelling of self-aeration in stepped spillways. Themodel of Lopes et al. [16] is taken as baseline, and a large simulation campaign is conducted in order toexplore its properties: robustness with respect to flow conditions, grid resolution, as well as sensitivity tothe model parameters. The simulations were performed for spillway geometries and inflow discharge valuesused in the experiments of Bung [4], and cover four step Froude numbers in the range from 2 . a) F s =2.7 – realisable k - (cid:15) (b) F s =2.7 – k - ω SST(c) F s =4.6 – realisable k - (cid:15) (d) F s =4.6 – k - ω SST(e) F s =8.3 – realisable k - (cid:15) (f) F s =8.3 – k - ω SST(g) F s =13 – realisable k - (cid:15) (h) F s =13 – k - ω SST
Figure 14: α l -fields (starting at the psuedo-bottom) illustrating the inception point locations predicted using Eq. (17)-(18)for the different F s cases on grid G3. The middle sub-figures show simulations using standard turbulence modelling, whilst inthe lower sub-figures the variable density formulation is used. The top sub-figures refers to results using standard turbulencemodelling and k c = 0 . / s , as suggested in [16]. The vertical lines indicate the experimental inception points. δ fs , which is the function used for limiting the activation region of the volumetric air entrainment sourceterm, see Eq. (9).The region of non-zero values of δ fs shrinks as the grid gets refined, and, in the limit,the source term is set to zero in the whole domain, regardless of flow conditions. To address this issue, anew formulation for δ fs is proposed, combining a parabolic profile with distance-based cut-off. Simulationsreveal that while fundamentally the results still depend on the grid resolution in the same manner, underthe new definition of δ fs the robustness of the model is improved.As an additional modification, amplifying the diffusion term in the α g -equation (10) is proposed in orderto account for the propagation of entrained air into corners of the steps. Results reveal that this leads toa significant improvement in predictive accuracy, the new model performing better than the original [16]across the whole considered range of F s numbers. Furthermore, the robustness of the model with respect togrid resolution improves significantly as well. It should be acknowledged that the selection of the value ofthe diffusion coefficient, C t = 150, is currently not physically motivated and can, therefore, be called intoquestion. Nevertheless, we believe that the possibility to use the same value across different flow conditionsand the clearly demonstrated advantages in terms of the performance of the model are sufficiently strongarguments in favour of adopting the proposed modification. Finding a more rigorous connection between C t and the characteristic scales of the flow remains as a line of future work.Finally, an algorithm for automatic estimation of the inception point is tested. The criterion for theinception point is based on energy balance, as proposed by [13]. The performance of the algorithm areheavily dependant on the underlying turbulence model. Unfortunately, for the four considered models thepredictions were not reliable, highlighting the need for a more careful investigation of what turbulencemodelling is appropriate for self-aerating multiphase flows. Acknowledgments
This work was supported by grant number P38284-2 from the Swedish Energy Agency. The simulationswere performed on resources provided by Chalmers Centre for Computational Science and Engineering(C3SE), the NTNU IDUN/EPIC computing cluster, and UNINETT Sigma2 – the National Infrastructurefor High Performance Computing and Data Storage in Norway. The authors are thankful to Daniel Bungfor sharing the data files from his physical experiments and to Pedro Lopes for sharing the source code ofthe airInterFoam solver. 23 eferences [1] Boes, R.M., Hager, W.H., 2003a. Hydraulic design of stepped spillways. Journal of Hydraulic Engineering 129, 671–679.URL: https://doi.org/10.1061/(asce)0733-9429(2003)129:9(671) , doi: .[2] Boes, R.M., Hager, W.H., 2003b. Two-phase flow characteristics of stepped spillways. Journal of HydraulicEngineering 129, 661–670. URL: https://doi.org/10.1061/(asce)0733-9429(2003)129:9(661) , doi: .[3] Brackbill, J., Kothe, D., Zemach, C., 1992. A continuum method for modeling surface tension. Journal of ComputationalPhysics 100, 335 – 354. URL: , doi: http://dx.doi.org/10.1016/0021-9991(92)90240-Y .[4] Bung, D.B., 2011. Developing flow in skimming flow regime on embankment stepped spillways. Journal of HydraulicResearch 49, 639–648. URL: https://doi.org/10.1080/00221686.2011.584372 , doi: .[5] Chanson, H., 1993. Stepped spillway flows and air entrainment. Canadian journal of civil engineering 20, 422–435. URL: https://doi.org/10.1139/l93-057 , doi: .[6] Cheng, X., Gulliver, J., Zhu, D., 2014. Application of displacement height and surface roughness length to determinationboundary layer development length over stepped spillway. Water 6, 3888–3912. URL: https://doi.org/10.3390/w6123888 ,doi: .[7] Clift, R., Grace, J., Weber, M., Weber, M., 1978. Bubbles, Drops, and Particles. Academic Press.[8] Dami´an, S.M., 2013. An extended mixture model for the simultaneous treatment of short and long scale interfaces. Ph.D.thesis. National University of the Littoral.[9] Dong, Wang, Vetsch, Boes, Tan, 2019. Numerical simulation of air–water two-phase flow on stepped spillways behindX-shaped flaring gate piers under very high unit discharge. Water 11, 1956. URL: https://doi.org/10.3390/w11101956 ,doi: .[10] Ervine, D., Falvey, H., 1987. Behaviour of turbulent water jets in the atmosphere and in plunge pools. Proceedings of theInstitution of Civil Engineers 83, 295–314. URL: https://doi.org/10.1680/iicep.1987.353 , doi: .[11] Fan, W., Anglart, H., 2020. varRhoTurbVOF: A new set of volume of fluid solvers for turbulent isothermal multiphase flowsin OpenFOAM. Computer Physics Communications 247, 106876. URL: https://doi.org/10.1016/j.cpc.2019.106876 ,doi: .[12] H¨ansch, S., Lucas, D., Krepper, E., H¨ohne, T., 2012. A multi-field two-fluid concept for transitions between differentscales of interfacial structures. International Journal of Multiphase Flow 47, 171–182. URL: https://doi.org/10.1016/j.ijmultiphaseflow.2012.07.007 , doi: .[13] Hirt, C., 2003. Modeling turbulent entrainment of air at a free surface. Flow Science, Inc .[14] Hirt, C., Nichols, B., 1981. Volume of fluid (VOF) method for the dynamics of free boundaries. Journal of ComputationalPhysics 39, 201 – 225. URL: http://doi.org/10.1016/0021-9991(81)90145-5 , doi: .[15] Issa, R.I., 1986. Solution of the implicitly discretised fluid flow equations by operator-splitting. Journal of ComputationalPhysics 62, 40–65. URL: https://doi.org/10.1016/0021-9991(86)90099-9 .[16] Lopes, P., Leandro, J., Carvalho, R.F., 2017. Self-aeration modelling using a sub-grid volume-of-fluid model. InternationalJournal of Nonlinear Sciences and Numerical Simulation 18. URL: https://doi.org/10.1515/ijnsns-2017-0015 , doi: .[17] Ma, J., Oberai, A.A., Drew, D.A., Lahey, R.T., Hyman, M.C., 2011. A comprehensive sub-grid air entrainment modelfor RaNS modeling of free-surface bubbly flows. The Journal of Computational Multiphase Flows 3, 41–56. URL: https://doi.org/10.1260/1757-482x.3.1.41 , doi: .[18] Matos, J., Pinheiro, A.N., Quintela, A.C., Frizell, K.H., 2001. On the role of stepped overlays to increase spillway capacityof embankment dams, in: Proceedings ICOLD European SymposiumDams in a European Context.[19] Menter, F., 1993. Zonal two equation k- ω turbulence models for aerodynamic flows, in: 23rd Fluid Dynamics, Plasmady-namics, and Lasers Conference, American Institute of Aeronautics and Astronautics. URL: https://doi.org/10.2514/6.1993-2906 , doi: .[20] Moraga, F., Carrica, P., Drew, D., Lahey, R., 2008. A sub-grid air entrainment model for breaking bow waves and navalsurface ships. Computers & Fluids 37, 281–298. URL: https://doi.org/10.1016/j.compfluid.2007.06.003 , doi: .[21] Mortazavi, M., Chenadec, V.L., Moin, P., Mani, A., 2016. Direct numerical simulation of a turbulent hydraulic jump:turbulence statistics and air entrainment. Journal of Fluid Mechanics 797, 60–94. URL: https://doi.org/10.1017/jfm.2016.230 , doi: .[22] Pfister, M., Hager, W.H., 2011. Self-entrainment of air on stepped spillways. International Journal of Multiphase Flow 37,99–107. URL: https://doi.org/10.1016/j.ijmultiphaseflow.2010.10.007 , doi: .[23] Qian, Z., Hu, X., Huai, W., Amador, A., 2009. Numerical simulation and analysis of water flow over stepped spillways.Science in China Series E: Technological Sciences 52, 1958–1965. URL: https://doi.org/10.1007/s11431-009-0127-z ,doi: .[24] Rusche, H., 2002. Computational fluid dynamics of dispersed two-phase flows at high phase fractions. Ph.D. thesis.Imperial College of Science, Technology and Medicine. URL: http://hdl.handle.net/10044/1/8110 .[25] Shih, T.H., Liou, W.W., Shabbir, A., Yang, Z., Zhu, J., 1995. A new k - (cid:15) eddy viscosity model for high reynolds numberturbulent flows. Computers & Fluids 24, 227–238. URL: https://doi.org/10.1016/0045-7930(94)00032-t , doi: .
26] Straub, L.G., Anderson, A.G., 1958. Experiments on self-aerated flow in open channels. J.Hydraul.Div. 84, 1–35.[27] Valero, D., Bung, D., 2015. Hybrid investigation of air transport processes in moderately sloped stepped spillway flows,in: E-proceedings of the 36th IAHR World Congress 28 June - 3 July, 2015, The Hague, the Netherlands, pp. 1 – 10.[28] Valero, D., Bung, D.B., 2016. Development of the interfacial air layer in the non-aerated region of high-velocity spillwayflows. Instabilities growth, entrapped air and influence on the self-aeration onset. International Journal of MultiphaseFlow 84, 66–74. URL: https://doi.org/10.1016/j.ijmultiphaseflow.2016.04.012 , doi: .[29] Viti, N., Valero, D., Gualtieri, C., 2018. Numerical simulation of hydraulic jumps. Part 2: Recent results and futureoutlook. Water 11, 1–18. URL: https://doi.org/10.3390/w11010028 , doi: .[30] Wardle, K.E., Weller, H.G., 2013. Hybrid multiphase CFD solver for coupled dispersed/segregated flows in liquid-liquidextraction. International Journal of Chemical Engineering 2013. URL: https://doi.org/10.1155/2013/128936 , doi: .[31] Yan, K., Che, D., 2010. A coupled model for simulation of the gas–liquid two-phase flow with complex flow patterns.International Journal of Multiphase Flow 36, 333 – 348. URL: https://doi.org/10.1016/j.ijmultiphaseflow.2009.11.007 , doi: http://dx.doi.org/10.1016/j.ijmultiphaseflow.2009.11.007 .[32] Zalesak, S.T., 1979. Fully multidimensional flux-corrected transport algorithms for fluids. Journal of ComputationalPhysics 31, 335–362. URL: https://doi.org/10.1016/0021-9991(79)90051-2 , doi: .[33] Zhang, G., Chanson, H., 2016. Hydraulics of the developing flow region of stepped spillways. I: Physical modeling andboundary layer development. Journal of Hydraulic Engineering 142, 04016015. URL: https://doi.org/10.1061/(asce)hy.1943-7900.0001138 , doi: ..