Numerical analysis of hydraulic jumps using OpenFOAM

ANN

artificial neural network

SPH

smoothed particle hydrodynamics

CFD

computational fluid dynamics

DNS

direct numerical simulation

LES

large eddy simulation

RANS

Reynolds-averaged Navier–Stokes

RNG

re-normalization group

SST

shear stress transport

FVM

finite volume method

VOF

volume of fluid

freestream velocity

compression velocity

shear velocity

pressure

time

inlet flow rate

Cartesian reference system component

body forces

fluid density

water depth

hydraulic head

energy drop

hydraulic jump efficiency

dimensionless water depth

turbulent kinetic energy

dissipation rate of turbulent kinetic energy

specific dissipation rate of turbulent kinetic energy

dynamic viscosity

turbulent eddy dynamic viscosity

kinematic viscosity

turbulent eddy kinematic viscosity

production of turbulent kinetic energy

effect of buoyancy

effect of dilatation

⁠,

modulus of mean rate-of-strain tensor

⁠, ⁠,

model parameters

⁠,

model parameters

acceleration of gravity

mesh element size

model parameter

water fraction in mesh element

Froude number

water-height-based Reynolds number

Mach number

Courant number

channel width

position of hydraulic jump toe

position of roller end

position of hydraulic jump end

dimensionless longitudinal coordinate

length of hydraulic jump

length of roller

generic flow property of fluid “i”

ratio of sequent depths

dimensionless wall distance

dimensionless velocity

error in variable i results

coefficient of determination

superindex relative to classical hydraulic jump

subindex relative to supercritical flow

subindex relative to subcritical flow

Characterizing and analyzing this phenomenon is of paramount importance from both the technical and the environmental point of view. Hager (1992) and Chanson (2013) performed extensive reviews on the attempts to study this phenomenon throughout history. Some of these works focused on a theoretical comprehension of the characteristic features of the classical hydraulic jump. The so-called classical hydraulic jump is defined by Hager (1992) as the hydraulic jump that occurs in smooth horizontal prismatic rectangular channels. Resch & Leutheusser (1972) performed a thorough study on air entrapment and energy dissipation processes depending on the inlet flow characteristics. Gualtieri & Chanson (2007) extended this analysis to a wider range of inlet flow conditions.

Most of the studies performed up until now have focused on the analysis of easily-measurable external macroscopic variables using an experimental approach. This can be explained partially by the difficulty in measuring certain variables using non-intrusive acoustic and optical methods in highly aerated flows (Murzyn et al. 2005). However, since the 1970s, coinciding with the emergence of computational fluid dynamics (CFD), more and more studies on the hydraulic jump are conducted by means of numerical methods. In this regard, computational techniques brought a brand new approach to water engineering modeling. For example, some hard-to-model phenomena, such as heat transfer (Thomas et al. 1990) or coupled biological processes (Muttil & Chau 2006), could be for the first time implemented thanks to numerical methods. This implied a whole paradigm shift and so the literature on this topic is vast: e.g., Ma et al. (2011), among others, modeled hydraulic jumps using different CFD techniques; Caisley et al. (1999) managed to reproduce accurately a hydraulic jump in a canoe chute using FLOW-3D; and Bombardelli et al. (2011) used this commercial software to successfully model a stepped spillway following a similar approach to that proposed here. Other approaches different from CFD have also been used in all kinds of water engineering applications, e.g., De Padova et al. (2013) successfully reproduced a hydraulic jump using smoothed particle hydrodynamics (SPH) techniques.

However, as Murzyn & Chanson (2009a) state, mathematical models still have problems in reproducing the physics of certain hydraulic phenomena, although they can contribute to their better comprehension. As Romagnoli et al. (2009) remark, an entire comprehension of the hydraulic jump internal flow features and turbulence structures has not been achieved so far. For Murzyn & Chanson (2009a), the main features of hydraulic jumps that have not been fully understood are the following: fluid mixing, bubble break-up and coalescence, free surface turbulent interactions, and wave formation and breaking processes.

Other authors have used a CFD approach to analyze hydraulic jumps in terms of shear stresses, potential erosion on stream boundaries and other more practical applications (Chanson 2000). Liu & Garcia (2008) published a model combining the CFD code OpenFOAM and the Exner Equation to model erosion and sedimentation processes in hydraulic structures using mesh deformation.

Nevertheless, despite the increasing number of publications in this area, in the numerical modeling of hydraulic structures and water engineering applications, deterministic models (e.g., CFD) are overwhelmed in number by their statistical counterparts, also known as ‘black box’ models (Chau et al. 2005). Thus, several authors used artificial neural networks (ANNs) to successfully predict the scour that occurs at hydraulic structures, such as bridge piers (Toth & Brandimarte 2011) or culvert outlets (Liriano & Day 2001). Taormina et al. (2012) and Cheng et al. (2005) used ANN to successfully predict aquifer discharge processes. Farhoudi et al. (2010) used fuzzy logic methods to analyze the scour downstream of stilling basins. Other authors used one-dimensional and two-dimensional approaches to reproduce the flow in similar geometries (Dewals et al. 2004). However, in hydraulic engineering, flows are generally strongly three-dimensional (Ahmed & Rajaratnam 1997). Therefore, the use of a fully three-dimensional deterministic model, such as proposed here, allows its application to a wider range of cases. Comparisons among different numerical methods to model hydrological and hydraulic phenomena can be found in the literature (Chen & Chau 2006; Wu et al. 2009).

The main goal of this work is to propose a fully three-dimensional CFD-based method to model classical hydraulic jumps using the open source platform OpenFOAM. The use of freely available open source codes allows a continuous community-based improvement of the model and avoids having to pay for costly software licenses. In this sense, other models can be found, also reproducing hydraulic jumps using OpenFOAM (Romagnoli et al. 2009; Witt et al. 2013). However, different outlet boundary conditions are herein presented. This change constitutes an asset as it allows the model boundaries to be brought closer to the hydraulic jump, which involves a significant saving of computational resources.

A case study is also conducted, where several variables of interest, such as hydraulic jump efficiency or roller length, are computed and compared to previous analytical and experimental studies. The sensitivity of the model to certain parameters, such as mesh element size, turbulence model used or boundary conditions, is assessed.

As discussed in later sections, given the accuracy of the results that have been achieved, this model is fully applicable to more complex geometries where hydraulic jumps have to be investigated, such as dam spillways, stilling basins and river rapids to name but a few.

In this model of the hydraulic jump, the geometry to discretize is rather simple: the domain consists of a prismatic rectangular channel. For this discretization, two big categories of approaches are normally used, namely: unstructured and structured meshing.

Unstructured meshes are generally better suited for a selective refinement, so preventing the over-refinement of regions where no large gradients of property are expected (Kim & Boysan 1999). Besides, these meshes fit better into complex geometries, show less closure issues and their arbitrary topology makes automatizing the meshing process easier (Biswas & Strawn 1998). Nevertheless, none of these advantages applies to the case under study, as the geometry is extremely simple and no mesh refinement is required.

Some authors state that mesh non-orthogonality does not affect results as long as the skewness of its elements is kept low enough (Huang & Prosperetti 1994). Nevertheless, structured meshes tend to be more accurate than their unstructured counterparts ceteris paribus (Biswas & Strawn 1998). Besides, structured meshing algorithms are generally more straightforward to implement and faster to execute. According to Keyes et al. (2000) structured meshing algorithms present a more regular access to memory, which significantly reduces its latency. Also, as discussed below, in multiphase flows, topologically orthogonal meshes with their axis aligned with the fluid interface tend to show fewer numerical problems. For all these reasons, a static structured rectangular hexahedral mesh is considered the best choice.

In some cases, meshes can be slightly refined in the vicinity of solid boundaries for accurately resolving the flow features in boundary layers, where larger property gradients occur. This may result in the formation of highly skewed elements, although this is not a real issue as long as orthogonality between the mesh axes and solid boundaries is ensured (Hirsch 2007). However, in this case, experience demonstrates that the mesh element size necessary to capture the freestream flow features is not smaller than that necessary to resolve the boundary layer features. As a consequence, cubic mesh elements of uniform size are used throughout the entire domain (see Figure 1). The optimum mesh element size is highly case specific, so it is determined by means of a mesh sensitivity analysis.

A wealth of algorithms has been developed to approximate numerically the Navier–Stokes Equations during the last decades. Nevertheless, none of them constitutes a perfect solution as their performance is highly case specific. Indeed, this is a topic extensively discussed in the literature (Jang et al. 1986; Barton 1998). It is important to remark that the algorithm performances are generally assessed in terms of computation requirements and stability as, eventually, all algorithms should converge to a similar solution. The most widely used algorithm to execute stationary simulations and normally the default option in all CFD codes is the SIMPLE algorithm (Patankar & Spalding 1972). Several improvements to its original implementation, such as SIMPLER or SIMPLEC, have been made since the model was developed. One of the most used variations is the PISO algorithm (Issa 1985). However, problems may arise when dealing with multiphase flows where phase changes are abrupt or the density difference is large (Brennan 2001). To overcome this issue, an algorithm was developed combining the best features of SIMPLE and PISO; the so-called PIMPLE (OpenFOAM 2011). This algorithm merges the outer-correction tools of SIMPLE with the inner-corrector loop of PISO in order to achieve a more robust and generalizable pressure-velocity coupling (Rodrigues et al. 2011).

Hence, the PIMPLE algorithm is used here as a good compromise between computation requirements and stability. This algorithm is implemented in OpenFOAM, a freely available open source platform constituted by all sorts of C ++ applications and libraries to solve all kinds of continuum mechanics problems (Weller et al. 1998). This code uses a tensorial approach and object-oriented programing techniques following the widely known finite volume method (FVM), first used by McDonald (1971). An in-depth explanation of the algorithm implementation can be found in the PhD theses of Jasak (1996), Ubbink (1997) and Rusche (2002).

The coexistence and interaction of several fluids and the way that the interface among them is defined is of paramount importance in numerical modeling of multiphase flows. Complex algorithms must be developed to model this phenomenon, whose stability and accuracy have a strong influence on the model final results (Hyman 1984). Surface tracking methods fall into two families of approaches, namely: the surface methods and the volume methods. On the one hand, surface methods explicitly define the free interface either using a Lagrangian approach, i.e., tracking a set of surface marker particles (Daly 1969), or using a Eulerian approach, i.e., defining functions that determine the free surface position (Osher & Sethian 1988). These methods present topology issues when dealing with highly deformed flows and breaking surfaces. For this reason, they are not considered appropriate to model hydraulic jumps.

In the present model, an interface compression algorithm is implemented in order to overcome the aforementioned issues. This method adds an extra term in the left-hand side of Equation (4): ⁠, where is a compression velocity with normal direction to the fluid interface. Multiplying by ensures that it will only affect those regions where the flow fraction is close to (Rusche 2002).

The aeration of a water flow modifies its volume, depth, density, and compressibility (Carvalho 2002), thus affecting the momentum transfer. It also reduces the scour risk by cavitation (Bung & Schlenkhoff 2010) and the shear stresses on the channel boundaries (Chanson 1994). Therefore, this is a phenomenon of paramount importance in highly aerated flows such as hydraulic jumps, bores, breaking waves, etc. Unfortunately, surface tracking methods per se cannot reproduce phenomena smaller than the mesh element size, such as bubbles and droplets, or the entrapment of large amounts of air (Toge 2012). To overcome this issue, additional air entrainment models are implemented. In low-aerated flows, a Eulerian–Lagrangian approach is possible, where the Navier–Stokes Equations are resolved and air is treated as a set of discrete particles. With larger air fractions, this approach is no longer possible and an entirely Eulerian method is necessary. Eulerian–Eulerian approaches yield better results than their Eulerian–Lagrangian counterparts in the latter case. Despite the fact that they require longer computation times, in entirely Eulerian approaches, buoyancy, drag and lift forces are taken into account. For this reason, in the case of the model proposed here, a Eulerian–Eulerian approach is implemented.

Turbulence features can either be resolved down to their lowest scales (direct numerical simulation (DNS)), if the mesh is accordingly fine, or modeled under a wealth of different approaches. Despite the application of DNS models to multiphase flows having been reported (Borue et al. 1995), in most cases turbulence features are partial or completely modeled in common engineering applications.

The RNG model formulation differs from that of the Standard ⁠, essentially, in the values of the aforementioned parameters. These changes seem to improve the model results to such an extent that, according to Bradshaw (1996), the RNG is the most used model in hydraulic applications.

Several authors claim that models are not suitable to model large adverse-pressure gradient flows (Menter 1993; Wilcox 1998). To overcome this issue, models were first introduced by Wilcox (1998). Their implementation is significantly different from that of ⁠, as the dissipation rate of turbulence kinetic energy is not modeled. Instead, a transport equation for its relative value is implemented. Among them, the SST (Menter 1993) proved to perform better than the Standard and the baseline (BSL) ⁠.

The suitability of one model or another is highly case specific and differences from using one model or another are normally remarkable. Hence, in order to determine which model performs best at a reasonable computational cost, a sensitivity analysis is conducted. To do so, simulations are run using the three RANS models discussed above ceteris paribus.

The boundary conditions imposed to force the hydraulic jump to occur consist of a supercritical flow inlet, a subcritical flow outlet, smooth bottom and side walls and an upper open air patch (see Figure 2). At the inlet, in order to fulfill the desired Froude number, a water depth and a potential velocity profile are imposed using a Dirichlet boundary condition. The pressure value is defined as a null von Neumann boundary condition, so forcing a hydrostatic profile. As regards the inlet variables of the RANS model, i.e., k, ⁠, and ⁠, they cannot be directly estimated from measurements. Instead, they are set to an arbitrary low value and a short initial stretch of channel is added in order for the flow to develop while approaching the hydraulic jump.

As regards the outlet, the subcritical water height that forces the hydraulic jump to occur within the simulated domain (⁠⁠) has to be imposed. This variable has to be obtained by iteratively testing values until the resulting hydraulic jump remains stable within the domain. Normally, a subcritical water height and a hydrostatic profile should be imposed at the outlet by means of a Dirichlet water level boundary condition. This, combined with a null Von Neumann boundary condition for velocity, would allow the flow to leave the domain freely. However, the imposition of a subcritical outlet by means of this approach in OpenFOAM appears to cause stability issues.

Indeed, to the knowledge of the authors, all cases of hydraulic jump simulations using OpenFOAM reported in the literature, such as Romagnoli et al. (2009) or Witt et al. (2013), have had to bypass this issue. To do so, they added an additional stretch of channel with an obstacle on the streambed, such as a step, a gate or a ramp, followed by a conventional free outlet.

In the present model, this problem is overcome by imposing a velocity profile at the outlet and so letting the hydrostatic profile develop, as is done at the inlet boundary condition. Assuming mass conservation, this approach univocally produces a given water height. This avoids having to model the aforementioned extra stretch of channel. As this implies bringing the boundary conditions closer to the phenomenon under study, comparative simulations are run in order to assess the model sensitivity to the boundary condition type.

A no-slip condition is imposed at the walls and roughness is not considered (Hager 1992). An atmospheric boundary condition is imposed at the top of the channel to allow fluids to enter and leave the channel. This is achieved by imposing a null Von Neumann condition to all variables except for pressure, which is set to zero (atmospheric pressure). Figure 2 summarizes the model's boundary conditions and some of its most relevant variables.

In terms of accuracy, the best choice would be to use a low Reynolds number model with no wall function at all. However, this implies refining the mesh to such an extent that the computational cost may become unaffordable. There is a vast literature on improvements to the original implementation of wall functions, such as Johnson & Launder (1982), but most of the solutions proposed have not been adopted by most CFD codes. This is due to the fact that, despite these approaches being valid from a theoretical point of view, many of them may cause stability issues (Blocken et al. 2007).

In this research, a high Reynolds number wall function for RANS models and smooth solid surfaces is implemented. The boundary layer in a case of these characteristics is likely to be slightly skewed (Taylor 1959). Nevertheless, as the flow mainstream direction is completely longitudinal, a bi-dimensional wall function is used for the sake of simplicity.

As regards the discretization schemes used to make the CFD model partial differential equations numerically approximable, a good choice always must be a good compromise between accuracy and stability. In spatial discretization, upwind models are generally preferred to downwind approaches as the latter tend to show severe stability issues. Compared to central differencing schemes, upwind approaches are slower, but also less diffusive and so more accurate. The problem is that, when abrupt property gradients occur, the latter schemes may require limiters in order to prevent spurious oscillations (Blazek 2005). Once limited, upwind schemes, such as Van Leer's (1982), are very appealing to discretize abruptly varied properties. The only drawback is that, when limited, these schemes become first-order accurate.

In the present model, the fluid fraction divergence terms are discretized using a limited Van Leer approach due to its abruptly variable nature. In the case of the RANS model variables, divergence terms are discretized using an unlimited upwind approach as they are less prone to cause stability issues. The velocity divergence terms are discretized using a central differencing scheme in order to avoid possible instabilities as well as to reduce computation times. Also, all gradient and interpolation terms of the model are discretized using this approach. Gaussian standard FVM is used to interpolate the variable values from cell centers to their faces. To save computational resources, as this mesh is strictly Cartesian, no orthogonality corrector is applied.

As regards the discretization of time derivatives, explicit schemes tend to be computationally lighter than their implicit homologues. However, they are also more unstable, especially, in simulations with skewed meshes (Blazek 2005) or when solving RANS equations (Lafon & Yee 1992). Therefore, implicit time discretization schemes are preferred in this model. This implies slightly longer computation times and eventual accuracy problems due to phenomena such as wave damping (Casulli & Cattani 1994), but also higher stability. Hence, a first-order accurate bounded implicit Euler scheme is used to discretize time derivative terms.

The time step length is variable throughout the simulation resolution process. Its value is automatically updated after every time step in order to ensure that the maximum Courant number never overcomes a threshold of ⁠, ensuring convergence and stability of simulations.

The variables analyzed and compared to previous studies are now discussed. The reference values to which model results are compared are denoted by a superscript asterisk (*).

A case study is conducted for validation purposes. The simulated case consists of a prismatic rectangular channel of dimensions (length, width, and height). The inlet flow is and the supercritical depth is Hence, the inlet velocity is ⁠. The subcritical depth is obtained following the procedure discussed above in this section. The density and the kinematic viscosity are and ⁠.

The approaching Froude and Reynolds numbers are and ⁠. A case of between and is considered optimum for model validation as it corresponds to the middle of the range of values recommended by the US Bureau of Reclamation (Peterka 1984). This approaching Froude number is assumed to be representative of the behavior of all stabilized hydraulic jumps, as described by Hager (1992).

As discussed above in this section, a mesh, turbulence and boundary condition model sensitivity analysis is conducted. Each of the turbulence models mentioned above (Standard ⁠, RNG ⁠, and SST ⁠) is tested in four different sized meshes. The mesh element sizes assessed are ⁠, ⁠, ⁠, and ⁠, which means meshes of ⁠, ⁠, ⁠, and million cells, respectively. To fulfill the wall treatment function hypothesis, it is checked in all cases that the coordinate mostly remains in the range of values between and ⁠.

A visual analysis of the model results leads to the conclusion that a stabilized hydraulic jump is reached (see Figure 3). All the characteristic features of this kind of jump described by Hager (1992) can be observed, namely: compact and stable appearance, low wave generation, gradual bubble deaeration, vortex formation within the roller, no flow separation in the entering jet, etc. Figure 3 shows how, downstream of the hydraulic jump where bubble deaeration occurs, hydrostatic pressure and velocity profiles quickly reappear. Also, the deaeration of large bubbles can be observed throughout the region where streamlines cut the water free surface. Downstream of that, despite waves and small bubbles not disappearing completely, the characteristics of developed flows can be observed again.

Figure 3 shows the wide span of bubble sizes that occur in the turbulent shear and the recirculation region of hydraulic jumps. Chanson (1994) found experimentally that the range of bubble sizes in hydraulic jumps can extend over several orders of magnitude. The average bubble size rapidly decreases longitudinally. This is due to the fact that large bubbles cannot stay long in the recirculation region as shear stresses break them up and buoyancy forces tend to expel them (Babb & Aus 1981). Small bubbles are not deaerated so quickly. Indeed, they can be dragged by advection forces throughout long distances until buoyancy finally expels them.

As discussed above, a mesh, turbulence, and boundary condition model sensitivity analysis is conducted in order to determine the best combination to achieve accurate results at an affordable computational cost.

As regards the outlet boundary condition used, Figure 4 shows examples of hydraulic jumps simulated using both approaches. A closer comparison between them shows no significant effect on the model outcome accuracy. Although an instant comparison, such as that in Figure 4, shows differences in water level profile and flow aeration, these differences completely disappear when results are averaged. No undesirable effects, such as wave formation, occur despite the new approach implying bringing the boundary conditions significantly closer to the phenomenon under study. The domain reduction achieves computation times up to 30% shorter in some cases.

The three tested turbulence models show small influence on the sequent depth estimations. The most accurate model is the RNG ⁠, followed by the SST and the Standard ⁠, although all errors are below ⁠. The inflexion point in the accuracy of the all models can be clearly observed at a mesh size of ⁠, thus the RNG model with size mesh is the most accurate approach.

As regards the estimation of roller lengths, the SST model appears not to be able to capture accurately this variable. The Standard model shows a reasonable accuracy (all errors are below ⁠) and low sensitivity to mesh size, which is an asseet. However, RNG is even more accurate and shows a perfect monotonically decreasing trend in errors, although the model is also highly sensitive to mesh size variations. The RNG model with size mesh appears to be the most accurate approach in the roller length prediction.

The prediction of the hydraulic jump efficiency achieves the highest accuracy values, with the error of all models below ⁠. The Standard with size mesh is the most accurate but, as observed in Figure 5, the sensitivity of this variable to the model parameters is extremely low.

All further discussion on the results and model validation is exclusively conducted using the results of the RNG model with size mesh. This is because this approach achieved the most accurate results in the most sensitive variable (i.e., roller length) while being reasonably accurate in the prediction of the less sensitive variables.

The mean subcritical water depth obtained from the CFD model is ⁠, which leads to a mean ratio of sequent depths of ⁠. Following Equation (10), the reference value for this variable in classical hydraulic jumps is ⁠. This means that the model yields a value approximately lower than that obtained using the expression proposed by Hager & Bremen (1989).

Regarding the mean roller length, the model yields a mean value of ⁠, being the value computed using Equation (11). This implies an underestimation of only ⁠. The accuracy in the prediction of this variable is crucial, as the largest shear stresses on the streambed generally occur within this stretch. In stilling basin design cases, this can be very helpful in order to determine the region of the structure that must be protected against scour.

Following Equation (13), the mean efficiency of the hydraulic jump is ⁠. This agrees with the result obtained from Equation (14) (Hager & Sinniger 1985), which estimates an efficiency of ⁠, with only a error. The good agreement in the computation of these three variables compared to other works is depicted in Figure 6.

The comparison of water surfaces to previous studies (Bakhmeteff & Matzke 1936) proves the consistency of the model presented herein. Figure 7 shows the mean water levels computed using the three different turbulence models. The most accurate RANS model is the Standard (⁠⁠), followed by the RNG (⁠⁠) and the SST (⁠⁠). It can be deduced from the above that turbulence models exert very little effect on the water free surface definition in average terms.

Nevertheless, an instant observation of the evolution of this variable shows that SST models produce a more unstable and bursting water surface with high bubble and spray production. Both models produce smoother surfaces, even though the Standard is also the model that yields a more uniform and less turbulent free surface.

Table 1 shows the accuracy of the model according to the variable predicted using the most accurate approach (RNG turbulence model with size mesh). It can be deduced from such accurate results that the model proposed herein is validated and so can be applied to real-life design cases.

A three-dimensional CFD model for transient multiphase incompressible flow is developed in order to predict the behavior of classical hydraulic jumps. After analyzing the effects on the results of several model parameters, such as the mesh refinement degree, the turbulence model or the boundary conditions, a stable hydraulic jump and accurate results are obtained. The model is built using exclusively free open source code, which implies avoiding expensive software licenses. Also, a problem found in other cases where OpenFOAM is used to model hydraulic jumps is addressed and solved. This problem involves the outlet boundary condition, where an additional stretch of channel with an obstacle attached to its bottom has to be added in order to force the subcritical flow to occur. Using the approach proposed in this paper, this additional stretch of channel can be removed by modifying the outlet boundary condition, with the subsequent saving of computational resources (up to 30% in some cases) with no effects on the model accuracy whatsoever.

A case study of approaching Froude number is simulated and the results are compared to previous studies of similar characteristics, such as that of Hager (1992) and Wu & Rajaratnam (1996). The roller length appears to be the most sensitive variable to model parameters (the SST is not even able to capture this magnitude). Some hydraulic jump variables are better reproduced in comparison with other authors' results, such as the sequent depth (error of only ⁠), whereas others show lower accuracies, e.g., the hydraulic efficiency (error of ⁠). The water free surface is accurately reproduced by all turbulence models in average terms, the Standard being the most accurate approach. An instant observation of the results shows that the SST model surface looks more turbulent than its counterparts. In any case, the accuracy of all of the variables analyzed is above in all cases so the model can be considered validated.

In the light of the results, the model is ready to be applied to real-life design cases, such as dam stilling basins, stepped spillways, river rapids, meandering channels, etc. As discussed above, the most accurate turbulence model in this kind of application is the RNG ⁠, although very fine meshes are necessary to ensure good performance and this model proved to be slightly slower than the Standard ⁠. The latter turbulence model could be a better choice in cases where low computational requirements are preferred without compromising accuracy excessively. The Standard also proved capable of reproducing the average water free surface slightly better.

As for future work, the model is currently being used in similar applications, both theoretical, such as triangular, circular and radial hydraulic jumps, and real-life cases. Also, the air entrainment and concentration distribution in hydraulic jumps is being studied using this model and compared with experimental data.

This research was conducted thanks to the funding provided by the VALi + D R&D Program of the Generalitat Valenciana (Spain). It would not have been possible without the contribution of Daniel Valero and Beatriz Nacher of the Hydraulics Laboratory of the School of Civil Engineering (Universitat Politècnica de València).