Analysis of Two Phase Flows on Stepped Spillways

Self-areation is a very relevant phenomenon that occurs in flows on stepped chutes, and is one of the reasons of building such structures. In this sense, the main purposes of stepped chute structures is to increase the energy dissipation of the flows by internal friction and to protect the bed of spillways through the presence of the air absorbed from the atmosphere. It is a “lucky coincidence” that, while the steps significantly increase the energy dissipation rate, they also increase the air uptake, reducing the risks of cavitation along the bed of the channel (Peterka, 1953; Frizell & Melford, 1991). An additional consequence of the use of stepped chutes is the reduction of the size of the energy dissipation basins at the toe of the spillways, implying in reduction of the total costs of the spillway structures. Physically, it is observed that the superficial air uptake begins at the cross section where the turbulent boundary layer attains the water surface, as described for example by Chanson (1997), and predictions of the location where water depth and boundary layer thickness coincide, and of the amount of air absorbed by the water, are important for the design of spillways. As a consequence of the abovementioned characteristics, stepped chutes are extensively used in hydraulics and environmental engineering. The air-water two-phase flow along the stepped channels is generally classified accordingly to its regime, as Skimming flow, Nappe flow and Transition flow, and the change from one regime to the other is still matter of discussion in the literature. In this chapter, a procedure is followed, aiming to adjust numerical tools in order to adequately represent the air-water interface along stepped chutes. In this case, it implies in finding conditions that allow the “breaking” of the surface, so that air can be incorporated by the water. Such a tool is interesting for the prediction of gas transfer in stepped spillways, and it is expected to be used by designers that need to quantify the so called selfaerated flows. As mentioned, the air is incorporated into the water through the surface, and knowing the distribution of the entrained air and the flow conditions, the gas transfer along the different positions in the cross section of the flow can be estimated, until it attains the bottom of the chute. In addition, having a good numerical code, it may be possible to obtain the pressure distribution along the structure and to localize the potential cavitation regions, allowing to “test numerically” possible solutions to enhance the air content and to reduce


Introduction
Self-areation is a very relevant phenomenon that occurs in flows on stepped chutes, and is one of the reasons of building such structures.In this sense, the main purposes of stepped chute structures is to increase the energy dissipation of the flows by internal friction and to protect the bed of spillways through the presence of the air absorbed from the atmosphere.It is a "lucky coincidence" that, while the steps significantly increase the energy dissipation rate, they also increase the air uptake, reducing the risks of cavitation along the bed of the channel (Peterka, 1953;Frizell & Melford, 1991).An additional consequence of the use of stepped chutes is the reduction of the size of the energy dissipation basins at the toe of the spillways, implying in reduction of the total costs of the spillway structures.Physically, it is observed that the superficial air uptake begins at the cross section where the turbulent boundary layer attains the water surface, as described for example by Chanson (1997), and predictions of the location where water depth and boundary layer thickness coincide, and of the amount of air absorbed by the water, are important for the design of spillways.As a consequence of the abovementioned characteristics, stepped chutes are extensively used in hydraulics and environmental engineering.The air-water two-phase flow along the stepped channels is generally classified accordingly to its regime, as Skimming flow, Nappe flow and Transition flow, and the change from one regime to the other is still matter of discussion in the literature.In this chapter, a procedure is followed, aiming to adjust numerical tools in order to adequately represent the air-water interface along stepped chutes.In this case, it implies in finding conditions that allow the "breaking" of the surface, so that air can be incorporated by the water.Such a tool is interesting for the prediction of gas transfer in stepped spillways, and it is expected to be used by designers that need to quantify the so called selfaerated flows.As mentioned, the air is incorporated into the water through the surface, and knowing the distribution of the entrained air and the flow conditions, the gas transfer along the different positions in the cross section of the flow can be estimated, until it attains the bottom of the chute.In addition, having a good numerical code, it may be possible to obtain the pressure distribution along the structure and to localize the potential cavitation regions, allowing to "test numerically" possible solutions to enhance the air content and to reduce cavitation damages (Gomes, 2006;Olinger, 2001).The obtained results presented in this chapter are still approximate, but already show that the problem is treatable using numerical tools.The main problems, as shown in this chapter, are related to the behavior of the free surface.This study considers the skimmimg flow regime, which generally occurs for large specific water discharges (Chanson, 2000).It is expected that the bed profile of the inlet structure affects the position of the inception point of the aeration (Boes & Hager, 2003a;Chanson, 2001).This study presents results of computational analyses for two inlet conditions of a stepped chute with an ogee crest.The aeration of the flows was investigated for an entry profile containing 1) steps with the same size, and 2)"transition steps" (steps with sizes increasing between the smooth bed and the stepped bed).In the last few decades, the computational fluid dynamics, CFD, has imposed itself as one of the best tools for the prediction of flow fields with accuracy and detail (Versteeg & Malalasekera, 1995).But even considering its success, complex flows still need many efforts to be completely understood.The description of the physical processes of a two phase flow involves a broad variety of computational methodologies to predict the quantities of the constituent components and the behavior of the interface.In the present study, the numerical procedures for the stepped channel flow simulations were conducted using the open softwares Salome and OpenFoam® (Open Source Field Operation and Manipulation), for mesh generation and CFD, respectively.In a general way, it may be commented that the implementation of CFD modeling is still limited, due, in part, to the relatively high costs of the commercial softwares.OpenFoam® is an open CFD software that helps to fill this gap.For the purposes of this chapter, the particularly relevant characteristic of the mentioned software is its ability to simulate flows interacting with hydraulic structures.The skimming flows along stepped chutes are inherently turbulent (Lima, 2003, Chamani, 2000), and their "turbulent characteristics" were investigated here using the "k-ε" model of turbulence, and adequate wall functions for the wall boundaries.The VOF method (Volume Of Fluid), was used for the interface-capturing procedures in the two-fluid model.Further, the PISO method (Pressure-Implicit with Splitting of Operators) was used as the pressurevelocity coupling scheme.In the calculations performed for this chapter, the rasInterFoam solver was chosen because, traditionally, the flow along spillways is considered incompressible (its velocity is much lower than the speed of the sound in water).The main purposes of this chapter were to show the implementation of the mentioned open softwares for the problem of aerated flows over stepped chutes, and to compare some calculated hydrodynamics characteristics with the literature, like the pressure distribution on the steps.The literature addresses the position of the inception point by empirical correlations (see, for example, Chanson (2006)).In the present study, numerical models are applied to evaluate the interface behavior, air concentrations along the flow, and to investigate the dependence between the inception point position and the bed profile characteristics.First results are shown, which show that these calculations are possible, and that more studies are necessary to improve the results.

Structure design
The parameters to be adjusted for stepped spillways during the design phase, such as step height, flow regime, and the desired air flow rate will of course depend on the application proposed for the steps.In general the stepped design takes into account the higher energy dissipation efficiency of the stepped chutes, when compared with smooth chutes, because it implies in reduction of costs of the dissipation basin (Arantes, 2007).Fig. 1 shows structures of spillways with different geometries, which illustrate that the purposes of the structure define the kind of spillway to be used.The spillway profiles used in this chapter are based on the WES profile (Waterways Experiment Station).The WES standard spillway profile is shown in Fig. 2. It has a basic slope of 0.78H: 1V, and the curved shape upstream from the crest is a composition of two circumferences.For the sizes considered in this study, Fig. 2 shows that the tangency point between the two circumferences occurs at (x= 0.61m, y=0.48m).Transition steps were proposed by García & Mateos (1995) after the crest.The steps grow continuously from a relatively small height to the constant step height.Fig. 2. WES standard spillway profile with the transition steps proposed by García & Mateos (1995).
The crest inlet profile strongly depends on the structure discharge operation capacities.And the flow over the structure is analyzed based in the structure inlet design to achieve a minor impact on the structure, to reduce cavitation risks and to optimize stilling basin.With the pressure diagrams, it is possible to compute an operation curve to optimize the flow over the structure.In this matter, all parameters as inflow conditions, reservoir volume, outflow discharge and the maximum discharge capacity are enrolled in the optimization process.As already mentioned, basically two distinct flow regimes occur on stepped spillways, as a function of the discharge and the step geometry (Povh, 2000).In the nappe flow, the steps act as a series of falls and the water plunges from one step to the other.Nappe flows are representative of low discharges capacities and large steps.On the other hand, small steps and large discharges are addressed by the skimming flow regime.This flow regime is characterized by a main stream that skims over the steps, which are usually assumed as forming the so called "pseudo-bottom".In the cavities formed by the steps and the pseudobottom, recirculation vortices are generated by the movement of the main flow.Perhaps the simplest form to quantify the transition from nappe to skimming flow is to express it through the ratio between the critical flow depth Yc and the step height Sh.Rajaratnam (1990) suggested the occurrence of the skimming flows for Yc/Sh > 0.8.This theme was also a matter of study of Stephenson (1991), who introduced a Drop term, D= q 2 /gSh 3 to distinguish between both regimes: for nappe flow, D< 0.6, and for skimming flow, D > 0.6.Chanson (2006) proposed limits for both flow regimes.He suggested, as a limit for the nappe flow, the approximated equation = .− .
, while the limit for the skimming flow was given by the following approximation: = .− . .This short presentation shows that this transition is still a matter of studies, and that a more general numerical tool is desirable to overcome the difficulties of defining a priori the flow regime along a stepped chute, (Chinnarasri & Wongwises, 2004;Chanson, 2002).Fig. 3 shows the main characteristics of the two flow regimes mentioned here.

Mathematical aspects 3.1 Free surface flow
The most difficult part of the interface simulation procedures is perhaps to obtain a realistic free surface flow solution.In a free surface flow, special numerical techniques are required to keep the position of the interface between the two phases.There are many free surface techniques available in the literature, which involve different levels of difficulties and several procedures to obtain a solution.In general, the numerical methods are under constant modifications, in order to improve their results and to avoid, as best as they can, nonphysical representations.Because the computational tool itself is under constant improvement, increasing both, the storage and the calculation speed capacities, this situation of constant improvement of the numerical methods is understood as a "characteristic" of this methodology of study.In this sense, the prediction of the behavior of interfaces is one of the problems that is being "constantly improved", so that different "solutions" can be found in the literature.Some of the procedures devoted "to the capture of the interface" introduce an extra term that accounts for the "interface compression", which acts just in the thin interface region between the phases.Some improvements related to the stability and efficiency calculations are described, for example, in Rusche (2002).In the present study a time step was adjusted to impose a maximum Courant number, and the prediction of the movement of the interface was viewed as a consequence of drag forces and mass forces acting between the phases.More details of the numerical procedures followed in this chapter are presented in the Section 4.

Self-aerated flow
The phenomenon of air entrainment and bubble formation is initiated by the entrapment of air volumes at the water surface, which are then "closed into bubbles", or really entrained into the flow.At the upstream end of usual spillways, the flow is smooth and no air entrainment occurs, so that it is called briefly as "black water".The flow turns into the so called "white water", or two-phase flow, only after a distance has been traversed, which may involve several steps, (see Fig. 4).The inception point of aeration is generally defined as the position where the boundary layer formed at the bed of the chute attains the water surface, (Carvalho, 1997;Boes & Hager, 2003b).The distance travelled by the water until attaining the inception point is commonly named as "black water length".Downstream, after a distance in which the air is transported from the surface of the flow to the bottom of the chute (named "transition length" by Schulz &Simões, 2011, andSimões et al., 2011), the flow attains a "uniform regime".Fig. 1b illustrates this uniform flow far from the inception point (in this example, it is impossible to attest for the uniformity of the velocity profiles, but it is possible to verify that the white-water global characteristics are maintained).The air uptake is generally described as a consequence of the turbulent movement at the water surface.The black water length (position of the aeration inception point) is a matter of continuing studies.Tombes & Chanson (2005), for example, furnished the predictive Eq. 1.
. * . ( It is based on the distance L i with origin at the crest of the spillway, the flow depth Y i , = (which represents the step depth per unit width), where q w is the discharge per unit width, g is the gravity acceleration, is the angle between the bed of the chute and the horizontal, S h is the step height and F * is the dimensionless discharge.Other expressions for the location of the inception point have been proposed by several authors.See, for example, Chanson (1994).The air concentration distribution, downstream of the inception point, may be obtained, for example, using a diffusion model (Arantes et al., 2010), as proposed by Chanson (2000), which leads to Eq. 4.

= − ℎ
´ C is the void fraction, tanh is the notation for hyperbolic tangent, y it is a transverse coordinate with origin at the pseudo bottom, D′ is a dimensionless turbulent diffusivity, K′ is an integration constant, Y 90 is the normal distance to the pseudo-bottom where C is equal to 90%.D′ e K′ are functions of the depth-averaged air concentration, C. As can be seen, such equations involve constants which must be obtained from measurements.
In this study, the transition between the black water and the white water is of the major interest.That is, the simulation of the transition of the smooth surface to the turbulent multiphase surface flow, and the calculation of the void fraction distribution along the flow are important for the present purposes.A sketch of the mentioned region is seen in Fig. 6.Fig. 6.Sketch of superficial disturbances and formation of drops and bubbles around the interface.

Two phase flow
A two phase flow is essentially composed by two continuing fluids at different phases, which form a dispersed phase in some "superposition region".According to the volume fraction of the dispersed phase, the prediction of multi-phase physical processes may differ substantially.According to Rusche (2002), the CFD methodologies for dispersed flows have been focused on low volume fraction so far.Processes which operate with large volumes of the dispersed fraction present additional complexities to predict momentum transfer between the phases and turbulence.
In the self-aerated flow along a stepped spillway, the flow behaves like an air-water jet that becomes highly turbulent after the air entrainment took place, where the flow can be topologically classified as dispersed.To computationally represent the localization of the interface with the dispersed phase (composed by drops and bubbles) is an extremely difficult task.Considering the traditional procedures, the optimization of the design of stepped structures crucially depends on the measurement of the inception point and the void fraction distribution (Tozzi, 1992).Both are mean values (mean position and mean distribution), obtained from long term observations.Energy aspects are obtained from mean depths and mean velocities (Christodoulou, 1999;Peterka, 1984).But, as such measurements are generally made in reduced models, the scaling to prototype dimensions may introduce deviations.When considering the cavitation risks, the scaling up questions may be more critical (Olinger & Brighetti, 2004).The main numerical simulation advantages rely, in principle, on the lower time consumption and lower costs, in comparison to those of the experimental measurements (Chatila & Tabbara, 2004).However, when simulating the flow, it is necessary to generate first a stable surface, continuous and unbroken, and then allow its disruption, generating drops, bubbles, and a highly distorting interface.Further, this disruption must happen in a "mean position" that coincides with the observation, and it must be possible to obtain void fraction profiles that allow concluding about cavitation risks.As can be seen, many numerical problems are involved in this objective.

Numerical aspects 4.1 Equation for the movement of fluids
The two phase flow along a stepped spillway is modeled by the averaged Navier Stokes equations and the averaged mass conservation equation, complemented by an equation to address fluid deformations and stresses.In this chapter the heat and mass transfers, as well as the phase changes, were not considered.
The flows are inherently turbulent and their characteristics were investigated here using the k-ε turbulence model and adequate wall functions for the wall boundaries.
The averaged Navier Stokes equations neglect small scale fluctuations from the two phase model.However, if the two phase flow has small particles in the dispersed phase, it has to be taken into account in the analysis, in order to achieve an accurate prediction of the void fraction.
Eq. 5 represents the Navier-Stokes equations for incompressible, viscous fluids (it is represented here in vectorial form, thus for the usual three components).It is not discretized in this chapter, because the rules for discretization may be found in many basic texts.But it is understood that it is important to show that the analysis considers the classic concepts in fluid motion.
is the velocity, is the density, is pressure, is time, represent the shear forces and represent the body forces.An outline of the boundary conditions applied to solve the equations is presented in Section 3.4.

Closure problem
As known, the averaging procedures applied to the Navier Stokes equations introduce additional terms in the transport equations, that involve correlations of the fluctuating components.These terms require new equations, which is known as the problem of "turbulence closure".Based on the Boussinesq hypothesis, it is possible to express the turbulent stresses and the turbulent fluxes as proportional to the gradient of the mean velocities and mean concentrations or temperatures, respectively.For the two phase flow it is still required to examine the effects of the dispersed phase on the turbulent quantities.There is a very wide spectrum of important length and time scales in such situations.These scales are associated with the microscopic physics of the dispersed phase in addition to the large structures of turbulence.The complexity of such flows may still be hardly increased if considering high compressibility and the simultaneous resolution of the large scale motion and the flow around all the individual dispersed particles.

Two phase methods
The physical representation of the air inception point is an application which considers the primary complex phenomena of the breakup of a liquid jet.When considering the simulation of the jet surface, it is necessary to track it.In general, the tracking methodologies are classified into different categories.We have, for example the Volume Tracking Methodology where the method maintains the interface position, the fluids are marked by the volume fraction and conserves the volume.The volume is also conserved in the moving mesh method, but the mesh is fitted to follow the fluid interface (Rusche, 2002).In this study an Euler-Euler methodology is applied, in which each phase is addressed as a continuum and both phases are represented by the introduction of the phase fractions in the conservation equations.An "interface probability" is considered, and closure methods are adopted to account for the terms that involve transfer of momentum between the continuous and the dispersed phases.In general, the numerical models are able to predict the mean movement of the free surface, but they fail to predict the details of the interface (which are important, for example, to incorporate air).In this study the disruption of the interface was imposed, in order to verify if it is possible to generate realistic interface behaviors and to obtain mean values of the relevant parameters.To attain these objectives, the conservation equations were discretized using the finite volume method, and the PISO (Pressure-Implicit with Splitting of Operators) algorithm was adopted as the pressure-velocity coupling scheme.Considering the Volume of Fluid Method, VOF, the fluids are marked by the volume fraction to represent the interface and it is based on convective schemes.The volume fraction are bounded between the values 0 and 1 (values that correspond to the two limiting phases).
As mentioned, in this study the VOF method was used.The volume of the fluid 1 in each element is denoted by V 1 , while the volume of fluid 2 in the same element is denoted by V 2 .Defining α = V 1 /V, where V is the volume of the cell or element, it implies that V 2 =1-.In this chapter, if the cell is completely filled with water, α =1, and if the cell is completely filled with the void phase, α =0.As usual, mass conservation equation (relevant for the mentioned volume considerations) is given by: Where, is the density, is time, is velocity.
A sharp interface can be achieved in the solver activating the term to interface compression.Eq. 7 illustrates the mass conservation equation with the additional compression term.
Where is a velocity field suitable to compress the interface.Literature examples show that the mathematical model for two phase flows used by interFoam (OpenFoam® two phase flow solver), allowed to obtain appropriated solutions when using the mentioned interface capturing methodology.For example, when simulating the movement of bubbles in bubble columns, two types of bubble trajectories were obtained: a helical trajectory, for bubbles larger than 2mm, and a zigzag trajectory, for smaller bubbles.Rusche (2002) mentions the agreement of the terminal velocity of the air bubbles with literature empirical correlations, which are also based on the bubbles diameters, with diameters between 1 and 5 mm.Although in the present analysis the problem of isolated bubbles is not considered, the mentioned agreement is a positive conclusion that points to the use of this method.

Boundary conditions
In the traditional CFD methodologies, the wall boundary condition is highly depending on the mesh size.The no-slip condition can be applied when the size near the wall is very fine.
On the other hand the slip condition is used when the near wall mesh is very coarse.Most of the times, the use of wall functions are appropriate and it imposes a source term at the boundary faces (Versteeg & Malalasekera, 1995).
For the inlet boundary, two conditions were adopted here.1) An initial condition similar to a dambreak problem, with a column of water having a predefined finite height above the weir crest.2) A constant water discharge having a uniform velocity profile.The water surface elevation upstream of the spillway crest was not specified as a boundary condition because this height is part of the numerical solution.The constant discharge (or constant flow rate) was imposed by a "down entrance" of water into the domain.
The dam break problem has been studied by theoretical, experimental, and numerical analysis in hydraulic engineering due to flow propagation along rivers and channels (Chanson & Aoki, 2001).However, in this study we were interested in the shape of the flow generated by an abrupt break of a dam gate, which then flows over the stepped spillway.Note that, if a constant water height would be defined upstream and far from the weir crest, only the transient related to the growing of the depth would be observed.So, a "water column" was imposed, and the growing and decreasing of the water depth along the spillway was observed.The initial boundary conditions and a subsequent moment of the flow can be visualized in Fig. 7a   At the outlet boundary, an extrapolation of the velocities was applied.It was applied locating the outlet of the flow far from the flow region of main interest.The local phase fraction varies accordingly to the diagram of Fig. 7, and it was observed that it tends to reach a more uniform characteristic at the structure toe (this was better observed for the constant inflow condition).When using the k-ε model, the turbulent kinetic energy and the energy dissipation rate must be imposed at the inlet boundaries.The actual value of these two variables is not easy to estimate.A too high turbulence level is not desirable since it would take too much time to dissipate.In this study, as the flow has a "visual laminar" behavior at the inlet, this condition allowed some simplifications.The closing equations originated from k-ε model are described by Eqs.8,9,10,11 and 12.Although known, the above equations are presented here to stress that the present chapter considers the classical ad hoc approximations for turbulent flows.The k-ε constant models are described at Table 1 in Section 6.1.2.

Simulation tools
As mentioned the open softwares Salome and OpenFoam® were used here for mesh generation and CFD, respectively.A Table with the software versions applied for simulations in this study is described at the Appendix II, as well as those respective websites for download.

Mesh generation
To simulate the flow over the domain, a structured mesh was generated at Salome software that is produced by OpenCascade.Some difficulties arise for the mesh generation in classical spillways, which are associated with the shape of the crest of the inlet structure (used to provide a sub-critical inflow condition).The structured mesh generation was a choice to have a mesh with hexahedral elements, which are highly recommended in the literature for treatment of free surface problems.There are many free softwares for mesh generation available for download, therefore many of them don´t have the ability to generate a hexahedral mesh.In this way, it must also be mentioned that some of the aroused difficulties may be related to the limitations of the specific software.In this case, a structured mesh was generated, based on the software characteristics for hexahedral algorithms.
Limiting the y+ value (One of the parameters that indicate mesh refinement), it is possible to reach more accurate results.However, the computational costs may significantly increase.Some general characteristics are mentioned here, as the case of adopting a too coarse grid, which leads to the situation that the results obtained are rather independent of the turbulence model, because the numerical diffusion dominates over the turbulent diffusion.
As an auxiliary tool to the mesh algorithm the domain were partitioned with many horizontals and vertical plans to direct the algorithm in the specific regions as small steps and crest.In this way the mesh is extremely refined in these regions, as shown in Fig. 8.The mesh used to represent the domain in this study has approximately 3,5 million of hexahedral elements.It is a tri-dimension mesh.However the domain in the z-direction is simulated for 1m of length and the mesh for this direction has a reduced number of eight columns.
After generating the mesh, and creating the faces to apply the boundary conditions at Salome software, it can be easily exported to OpenFoam® through the ".*UNV" file format.
At OpenFoam® all mesh properties can be checked and at the boundary file at the "\case\constant\polymesh\" path the boundary condition properties for any face can be visualized and edited.

Simulation and flow details 6.1 Setup
In this study, OpenFoam-1.7.1 were used to conduct the numerical results (as described in Appendix II).All the numerical analyses were processed considering that the chute is composed by a prismatic channel with rectangular cross section and having a stepped bed.
The flow calculations involved discharges varying from 0.5m 3 /s to 20m 3 /s.The discharge condition directly influences the flow regime over the structure.This range was used to test the performance of the program and the adopted procedures.
The inlet conditions for the gas phase fraction and the liquid velocity were taken directly from inlet water discharge measurements found in the literature.By The maximum specific discharge applied for a good hydraulic performance is 25-30 m 3 /s m (Boes & Minor, 2000).The interface momentum transfer (Rusche, 2002) is used to account for the lift force and predict the phase fraction.
The surface tension coefficient between the two phases is set with the value of 0.07N/m.The gravitational acceleration had the value of 9.81 m/s 2 .The outlet velocity was fixed with a zero gradient far from the interesting flow region and all walls were treated with wall boundaries lawyers.The water properties were considered at the temperature of 298K.
The time precision was automatically adjusted by the solver.Initially it was set to have the value of 0.0001s.Therefore with the solver feature to automatically adjust the time step, it suffered changes trough the time to attempt the maximum Courant Number specified as a precision acceptable to the numerical simulation.
In the phase fraction field the solver allow to set a number of sub-cycles in which the phase fraction equation is solved without doing an extremely reduction in the time step and hardly increasing the cost with time precision.
In the PISO algorithm, the number of correction for the pressure was set to three.To the initial conditions group of sets, all patches defined as faces in the process of mess generation to represent the described domain in the Salome software had a value assigned for volume fraction, pressure and velocity respecting the interaction between them to represent the initial physical properties of the domain described in the Section 4.4.The divergence terms in the velocity equations and in phase fraction equations were discretized using a central difference scheme.While for the Laplacian terms a Gauss upwind scheme was defined.For temporal discretization a Crank-Nicolson method was used.

Discretization schemes
An accurate numerical scheme is very important to obtain a solution.In order to check the influence of the numerical scheme on the solution, two different configurations for the numerical discretization have been used.With the present model, the two phase flows are better predicted with a higher order scheme.The first order scheme is not able to correctly predict the flow around the spillway crest without producing very large oscillations, as shown in Fig. 9. On the other hand, the second order schemes produce a large and stable recirculation bubble zone between the steps (while the first order did not generate such recirculation zone).Both schemes were used, because the main characteristics of the flow were still being adjusted, but it is necessary to use a more adequate scheme to reproduce the interfaces between air and water (free flow and between the steps).Due to the large size of the domain, and the related computational costs, only some more critical regions were refined.The numerical mesh was refined, for example, around the ogee profile, where the oscillations were observed.The different discretization schemes and mesh refinements were also used to check the sensibility of the model.It was observed, as expected, that not refined discretization lead to unreal interface waves, which increase significantly.As pressure condition in the flow, the hydrostatic pressure gradient was used.Further, most of the tests in this study were run using an upwind scheme.Depending on the boundary conditions applied, the flow can be always supercritical, or it can achieve a supercritical condition inside the domain.In this case, the relevant boundary condition is the discharge applied at the inlet, which directly influences the flow regime.The distance covered by the subcritical flow, for the adopted inlet conditions, allowed investigating the surface characteristics of the flow along the hydraulic structures for both the proposed inlet profiles, as shown in Fig. 5a and Fig. 5b, and Fig. 10a and Fig. 10b, respectively.Although in all situations a region with "unbroken liquid" is observed, followed by the broken liquid, the calculated surface did not maintain a smooth characteristic, although with smaller instabilities in relation to Fig. 9.As the breaking of the interface was an objective of this study, the numerical schemes were not chosen to guarantee the unconditional stability of the surface.However, as mentioned, they produced instabilities from the very beginning of the accelerated region (spillway entrance), and not only after attaining the position where the boundary layer coincides with the surface of the liquid.
The algorithm PISO used in this study is based on the assumption that the momentum discretization may be safely kept through a series of pressure correctors.It needs small timesteps.As a consequence, the PISO algorithm is also sensitive to the mesh quality.Aiming to guarantee the convergence, the Courant number, Co, must keep a low value (it is function of the time step, the magnitude of the velocity through the element, and the element size) (Versteeg & Malalasekera, 1995).In order to check the influence of the mesh size in the solution, a more refined mesh can be tested.

Turbulence properties
As mentioned, the k- averaged turbulence model was used.As expected in the RANS computations, the simulations converged to steady solutions.Table 1 shows the values of the constants usually adopted for the k- model (see Eq. 7 through 11).

Pressure diagram distribution
The pressure distributions along the steps are important to study the risk of cavitation in stepped chutes (Franc & Michel, 2004).It is not only the position of the lower pressure that is important, but also its value and frequency.In particular, the results calculated in the present study show no qualitative differences between the measured literature data, when considering the position and the mean pressure distribution analysis (Amador, 2005).As mentioned, the frequency of the low values is important, which points to the need of instantaneous values, and not only the mean profiles.
The negative pressure plotted in Fig. 11 were also predicted by Chen et al. (2002) for a WES structure with the transition steps.The pressure distributions over the steps in a stepped chute depend on the position of the faces of the steps.In the vertical face a separation flow region occurs and the pressure achieve negative values.Over the horizontal face the pressure distribution is influenced by the recirculation in the cavity, which results in a positive value for the mean pressure.The last step of Fig. 11a, amplified in Fig. 11b, shows this difference between the behavior of the pressure along the horizontal and vertical surfaces.

Velocity distribution
The velocity of an aerated flow is expected to be higher than the velocity of an unaerated flow, because the entrained air reduces the wall friction (Steven & Gulliver, 2007).Fig. 12 shows a velocity diagram obtained in the present simulations.The values obtained, though not directly compared to the results of other sources Cain & Wood (1981), are of the same order.Amador (2005) described results of experiments of velocity fields over steps located in the developing flow region, were the growth of the boundary layer was analyzed (upstream of the inception point).The study shows the presence of large size eddies at the corner of the step faces, and recirculation areas.Recirculation and air concentration distributions were also investigated by Matos et al. (2001).The authors mentioned differences for the flows between the step edges and in the steps cavities, describing it and mentioning the effects of turbulence intensity.
In this study, a calculated velocity field is shown in Fig. 12, where the region of maximum velocity, for example, is easy to be observed.As can be seen, the mentioned recirculation is adequately reproduced, which corresponds to the negative horizontal velocity that can happen between the steps.Fig. 12. Velocity distribution along the structure.

Void fraction distribution
As described before, the attention of this chapter is more concentrated in the phenomena related with the air transfer between the two phases interface.Fig. 13 shows an instantaneous of the position of the interface and of the void distribution, obtained following the present procedures.In this case, the example is for constant steps.The imposed flow rate was 15 m 3 /s m.As can be seen, the surface is highly distorted, reproducing the conditions found in real flows.
The deformations characterize entrapped air, as defined by Wilhelm & Gulliver (2005), and not entrained air, but they also allow to calculate a mean "void concentration", a procedure followed to check the methodology which is being tested here.As showed by Lima et al. (2008), measured concentration profiles may involve a relative high percentage of entrapped air, which may induce to mistakes in the calculations of the air flow rate introduced by the water flow.In this sense, the present methodology shows that this "entrapped air concentration profile" may be calculated to correct the evaluation of the entrained air.In  other words, the air flow discharge is obtained.This quantification is relevant, allowing to check experimental results and predictive equations obtained from simplified models, conducting (in a more complete stage of these studies) to a better understanding of the air entrainment phenomena on stepped chutes.In addition, concentration fields like that shown in Fig. 15a will allow to correlate pressure and void fraction distributions, and to recognize regions where the air content is not enough to reduce cavitation damage (that is, the the regions under risk of cavitation).From Peterka (1953), the air concentration that is sufficient to avoid cavitation damage is about 5 to 8%.The inlet of air increases the mixture compressibility and the flow becomes able to absorb the impact of collapsing vaporized bubbles.Fig. 15a also shows the concentration regions defined by Chanson (1997), considered here to verify the possibility of obtaining practical information of such simulations.Accordingly to Chanson (1997), three concentration ranges can be considered in the concentration analyses, as shown in Table 2.

Continuing studies
The results presented in this chapter show that the proposal of studying the breaking of the air-water interface and the formation of a dispersed phase is a challenging issue, and is still not definitively quantified.It was possible to show that instabilities can be produced at the interface in a way to maintain the main flow stable, but these instabilities did not follow a realistic behavior.The aim is to break the surface in the position of the inception point of aeration, but the instabilities were formed in the very beginning of the accelerated region of the flow, at the spillway entrance.
Bubbles and drops were "formed" in the numerical simulations, as shown in Fig. 16.However, the sizes of both depend of the numerical mesh used, and, in this sense, the reproduction of these geometrical characteristics is still limited by the mesh.In practice, smaller drops and bubbles are formed.The numerical solution were processed with the discharge of 20 m 3 /s m.Depending on the conditions of the flow, the abovementioned characteristics are damped, or not observed.This is the case shown in Fig. 17, where the instabilities were minimized and no drops or bubbles were formed along the spillway.As a conclusion, the hydrodynamic characteristics still need to be better reproduced in the simulations, in order to guarantee the formation of the dispersed phase.The numerical solution were processed with the discharge of 15 m 3 /s m.
Fig. 17.Velocity distribution along the structure.In this case, no drops and bubbles are formed.
Numerical simulations over stepped spillways which take into account the stable and unstable regions of the interface, disrupting it and forming the dispersed phase, have not been yet presented in the literature.However high quality measurements of different characteristics of such flows, including cavitation predictions, are found in the literature, and can serve as guide for the mentioned simulations.

Conclusions
Skimming flow over stepped chutes were simulated to verify the possibility of obtaining results of air entrainment in stepped spillways.The stable/unstable behavior of the airwater interface was discussed, and the difficulties to obtain a realistic breaking of the interface and the formation of the dispersed phase were exposed.Two situations related to the geometry of the spillways were considered: constant step sizes and varying step sizes.Also two situations related to the flow conditions were considered: constant inlet water discharge, and a "sluice-gate-brake", that is, a kind of "dam-break" situation flowing over the stepped spillway.
Results of pressure fields, velocity fields, void fraction profiles, along the simulated spillways were furnished, showing behaviors similar to those observed and simulated in the literature.For the void fraction profiles, also a "entrapped air" profile was furnished, showing that the methodology allows to confirm experimental studies about entrapped air found in the literature.Unrealistic behavior of the surface instabilities were observed.They were formed in the accelerated region of the flow over the spillway, and for some flow conditions they produced large disturbances, a negative characteristic that still need to be solved.Drops and bubbles could be formed in the flow over the spillway, but their size was shown to be limited by the mesh characteristics.The two fluid model adopted neglects the effects of large scales and uses an average momentum transfer between the phases.Although the surface deformation could be reproduced, the interface tracking calculation showed to be difficult to apply when considering the formation of a dispersed phase.The open codes "OpenFoam®" and "Salome" were used, showing that such codes are capable to furnish useful results for engineers and researchers in hydraulics and fluid mechanics.

APPENDIX I: Software information
The OpenFoam® software is a continuum mechanics simulation tool, but it is also an objectoriented numerical code, and a toolkit for programming.It is written in C++ programming language and can be freely download at the website https://www.opencfd.co.uk.The CFD library is extensive and incorporates several types of solvers.The software also contains tools for mesh conversion from many mesh generation softwares, including commercial and free codes.The open source code has the advantage of the "transparency", that is, it allows following the methodology adopted for the calculus, and the assumptions made to obtain a solution (ever present in any code).The results also can be converted to be visualized and post processed in a number of different softwares.A version of the Paraview software is integrated in OpenFoam® as its standard post processing tool.

APPENDIX III: Numerical procedures in CFD software
The numerical simulations using the OpenFoam® software are basically performed by file editing.The software is structured using three basic directories.The constant subdirectory contains all properties folders for the fluids settings and mesh characteristics.Some physical constants.For InterFoam solver this directory contains three files designed as: transportProperties, enviromentalProperties and dynamicMeshDict and a specific folder for mesh properties named as polyMesh.Into the transportProperties all information about the properties of the two phases can be edited.The dynamicMeshDict file allows to enable the mesh movements.
The polyMesh constant sub directory keeps all information about the mesh.
In the OpenFoam® the system sub-directory contains exactly five files: ControlDict, DecomposeParDict, fvSchemes and setFieldsDict.The contolDict file allows users set the simulation process information like starting, end and time step.User also can choose the methods to save the simulation data.In the fvSchemes system file, the descritization schemes can be setted.The fvSolution file is responsible for all solver settings like the number of iterations.
In the free surface problems, the initial condition for the volume fraction of the two fluids also needs to be defined.The setFieldsDict allows setting different density values for each part of the domain.The decomposeParDict file describes the decomposition for multiprocessing simulation.

Fig. 4 .
Fig. 4. Self-aeration: inception point and the distinct boundary layers.Fig.5 shows two sketches of the position of the inception point, for the two geometries considered here (only large steps in figure5a, and transition steps in figure5b).As can be seen, air is captured by the water only after the boundary layer has attained the surface.

Fig. 5 .
Fig. 5. Boundary layer growth and the differences of the inlet profiles used in this study.a) With same size for all steps, b) Initial steps with smaller size.
and Fig 7b.The second moment was taken close to the end of the flow of the phenomenon.(Physically, it corresponds to the time of 16s).

Fig. 7 .
Fig. 7. Phase fraction diagram.a) Initial water column (As an initial condition for phase fraction, InterFoam solver requires that both phases exist into the domain, at least into some volume cells), b) Water discharge following the structure slope.

Fig. 9 .
Fig. 9.The influence of the first order discretization scheme along the first part of the surface for the flow over the spillway with varying step sizes.

Fig. 10 .
Fig. 10.Inlet profile: a) Initial structure steps with smaller size.b) With same size for all structure steps.

Fig. 11 .
Fig. 11.Negative Pressure.a) Along stepped chute, b) In a step along the structure.

Fig. 13 .
Fig. 13.Skimming flow regime.Result for constant steps at the entrance.

Fig. 14 .
Fig. 14.Volumetric fraction along the flow.The position 0.0 corresponds to the bottom, and the position 1.0 corresponds to the surface.a) Histogram of void fraction where alpha1_average represents a statistic temporal alpha average, b) presents first results of air entrainment at any location of the spillway, allowing the prediction of the air concentration profiles (void fraction profiles) for all calculation domains.

Fig. 15
Fig.15is an example of shown in Fig.14b, obtained for the position indicated in Fig.15b.With the distribution of entrained air, the gas transfer along the flow can be estimated, or, in Classification C < 0.3 to 0.4 Clear water comprising air bubbles C > 0.6 to 0.7 Air flow comprising water droplets C > 0.7 Two-phase flow with equal air and water contents Table 2. Main regions of concentration profiles.(Chanson, 1997) As already mentioned, two examples are considered: The opening of a dam gate, a transient situation, and a steady-state simulation that allows evaluating changes of the void ratio along the self-aerated spillway flow.Air concentration profiles C(y) are taken in several cross sections at steps downstream from the inception point.

Fig. 15 .
Fig. 15.Phase fraction a) Contours of isosurfaces, b) Distribution of void fraction for the 4 th step.

Fig. 16 .
Fig. 16.General characteristics observed in a simulated flow: 1) Unrealistic superficial instabilities in the accelerated region of the flow; 2) Formation of drops and bubbles, as expected; 3) Cavities between the steps.

Table 3 .
Some useful links with open software information are addressed in Table3with the respective version of software applied in this study.Open software.