Hydraulic jump in open channels occurs when flow transits from supercritical to subcritical. The nondimensional Froude number (Fr) which is the ratio of inertial to gravitational force determines if the flow is supercritical (Fr > 1) or subcritical (Fr < 1). Flow in this transitional region is turbulent and accompanied by air entrainment coupled with vortex development. In civil engineering applications, a hydraulic jump is created to dissipate the energy in the flow. Any standard textbook in fluids can provide relevant background theory on hydraulic jumps.
Over the last three decades, numerical modeling of open channel flows with hydraulic jump has drawn the attention of many researchers. The complexity of the solved flow equations ranged from Bernoulli’s energy equation to three dimensional Navier-Stokes equations. Hydraulic models typically solve either the energy equation or shallow water equations. General-purpose hydraulic models are relatively simple to use and can broadly serve the purpose for many applications. Popular one-dimensional hydraulic models include HEC-RAS, WSPG, DHM, MIKE 11, TELEMAC and SWMM among others. For surface water applications where the flow is predominantly in one direction, these models have been shown to provide a reasonably accurate simulation. These models can be applied for testing various “what if” flow scenarios and are not constrained by the required CPU time.
The availability of increased computation resources has given rise to Computational Fluid Dynamics (CFD) software, and these algorithms are providing opportunities for researchers to capture the physics of flow at microscales. CFD tools focus on solving the three dimensional Navier-Stokes equations across varying spatial and time scales. The CFD models are providing new windows of opportunity to better understand the flow in complex situations in engineering and science disciplines. Application of CFD models needs significant amounts of dedicated computational resources. However, since they can better model the flow, the involved computational costs are a fraction of the prototype physical models, which will continue to motivate the practioner audience to use CFD models. The code in each of these models solves a system of equations based upon conservation of mass, energy, and momentum, typically using either finite difference, finite volume or finite element numerical techniques. Popular CFD models include OpenFOAM, FLOW-3D, and TUFLOW.
In this work, we compare the performance characteristics of the two popular hydraulic models HEC-RAS and WSPG with the CFD OpenFOAM, for a flow situation with a hydraulic jump. The outline of this paper is as follows. In the literature review, the focus is on publications that used OpenFOAM for modeling hydraulic jump and other civil engineering applications. All three models have been briefly described. In the review of the CFD OpenFOAM, the focus was laid on its focus on turbulence modeling components. The target test problem and the associated boundary conditions are next detailed. The CFD results of Bayon et al.  have been used as reference data for validating the performance of the selected hydraulic models. In the results section, the depth profiles from the three models are analyzed and result from sensitivity analysis on channel bottom roughness coefficient illustrated.
2. Literature Review
Bayon et al.  simulated the salient characteristics in a hydraulic jump (Froude number = 6.1) by solving the three-dimensional equations using OpenFOAM. The turbulence was modeled using three Reynolds-averaged Navier-Stokes (RANS) models: Standard k-epsilon, RNG k-epsilon, and SST k-omega. Various jump characteristics like sequent depths, efficiency, roller length, free surface profile, were compared with published studies. Bayon et al.  used OpenFOAM to model the characteristics of a hydraulic jump and compared the model results with their experimental data. While using experimental data across a rectangular channel in a laboratory setting has been a common practice, they conducted experiments on a prototype channel to simulate diversion of flow due to the construction of high-speed rail infrastructure. Their prototype channel was configured to simulate a combination of curved transition, stilling basin, weir, and a stabilization reach. The air-water interface was defined using a Eulerian volume method. Using their experimental data as a benchmark, they compared the surface profile among different turbulence models at the target stilling basin.
At the transition of supercritical and subcritical flows, because of the turbulence, pockets of air are captured in water, and they move in recirculatory motion. Simultaneous bubble breakup and their mergers occur in the turbulent shear section of the recirculating region. Large air bubbles, because of velocity gradients, can experience multiple breakups. If the bubble is big, buoyancy lifts it to the surface, while smaller bubbles remain in the lower portions. The volume of the air bubbles changes continuously. Witt et al.  modeled the two and three-dimensional void fraction and bubble size in a transient hydraulic jump with a Froude number of 4.82. They used a realizable k-epsilon turbulence model. Their fluid dynamics video depicts the air entrainment characteristics and bubble behavior within the hydraulic jump. All these characteristics were also present in their OpenFOAM CFD model.
Bayon et al.  assessed the accuracy of the solution resulting from their OpenFOAM and FLOW-3D models with experimental and published data. They compared the 3D swirling turbulent flow details of the jump with an incoming Froude number of 6.0. Their numerical results show that while OpenFOAM can better reproduce the structure of the jump, the strength of FLOW-3D model lies in its ability to better model the interaction between the supercritical and subcritical flow along with the derived variables. Both the models could nor reproduce the physics of flow in the roller region, where flow swirling occurs. Castillo et al.  modeled free and submerged jumps (Froude numbers of 2.69 and 1.80) downstream of the sluice gate using ANSYS CFX, FLOW-3D, and OpenFOAM by comparing the model results with their experimental data. The chosen turbulence closure model, k-ε was common for the three codes. In ANSYS CFX and OpenFOAM, the air and water were solved as a homogeneous model, and in FLOW-3D, the free surface was modeled with a simplified Volume of Fluid technique. The results indicate that the reliability of the models depends on the flow variable of interest and the location of the chosen cross-section within the jump. All three models could not reproduce the experimental observations in the recirculation area of the jump. The authors concluded that the k-ε turbulence model is not appropriate for solving flows with jumps, and they recommended using other turbulence models like k-ω.
Romagnoli et al.  modeled the hydraulic jump created by a sluice gate by solving the RANS equations with interFOAM solver. Their numerical values for the streamwise velocity and turbulent energy are close to the experimental data. Martins et al.  simulated both using OpenFOAM and experimentally, hydraulic jump in a gully located in a rectangular channel. Their focus was to analyze the structural part of the gully, as the standard drainage models use simplified models. The results provide a complete three-dimensional insight into the hydraulic behavior of the flow inside the gully. Lopes et al.  investigated the ability of OpenFOAM to reproducing drainage efficiency of a continuous transverse gully with a grate, across a combination of 40 varying flow rates and slopes. The validation of the numerical results was accomplished by running the model to experimental real-scale data sets. The close agreement of the results showed the promise of the CFD model to act as an alternative to costly experimental investigations.
Egea  experimentally and numerically investigated submerged type hydraulic jumps and the associated turbulence coupled with air entrainment. The RANS equations were closed using the k-e turbulence model using OpenFOAM. The InterFoam numerical solver was used for solving the equations. This solver uses a single set of Navier-Stokes equations for the two fluids independently for water and air, and additional equations required to simulate free-surface where the velocity is shared by both the phases. After initializing the required flow and model parameters, the solver first calculated the time step is calculated based on Courant condition, and the adaptive time control technique implemented. At any time step, the fraction equation is first solved to guarantee the smoothness of the air-water interface. The mass flux and the density of the cells are obtained for the given time step. With this information, the momentum equation is solved to yield the velocity profiles. The solution is then advanced to the next time step until the desired time level is reached. The CFD model results indicate that the downstream bed velocity peak, in terms of magnitude, at specific locations, is well predicted. However, the model could not accurately simulate the spatial dissipation of the bed velocity peak and surface velocity components.
Other investigators who have used OpenFOAM for hydraulic applications include the works of Teuber et al.  who modeled flow over a two-dimensional ground sill. They tested different RANS and LES simulations and compared the numerical results with analytical and experimental data. Kramer et al.  studied the fundamental physics of a two-dimensional, flat plate at low Froude numbers by performing nonlinear numerical simulations. The results were compared with a linearized potential-flow model in order to determine the effects of nonlinearity and wave breaking. Beg et al.  used OpenFOAM to study the effect of manhole surcharge on storm drain manhole head loss coefficients and manhole-gully discharge coefficients. The CFD model results were compared with discharge and water depth data at the manhole and the ADV velocity data at the gully, measured from the physical model. They used the interFoam solver to track the free surface or interface location between two fluids. The RANS equations were solved, and turbulence phenomena were simulated using the Standard k-ε model and RNG k-ε model to replicate the turbulence condition at the gully and manhole respectively. Schulze and Thorenz  reviewed the strengths and weaknesses of OpenFOAM for hydraulic applications. They reviewed the functionalities and capabilities of the software for hydraulic engineering applications, including a short description of the meshing process, the numerics of the solver as well as a short overview of the applicability and the limitations of the interFoam solver. They underscored the importance of grid quality and chosen discretization, as crucial parameters that affect the accuracy of the results. Small grid size can enable capturing minor bubbles or droplets. In the absence of appropriate mesh size, the software cannot model air entrainment or bubble transport and detrainment, which are encountered in some hydraulic applications. Higuera et al.  validated the ability of OpenFOAM to simulate multiple physical processes. These processes include three-dimensional dynamic pressure induced by a solitary wave on a vertical structure, transient wave group, rip current on a barred beach and run up on a conical island.
3. Numerical Models
The salient characteristics of the three numerical models considered in this study are listed below.
Open Field Operation And Manipulation (OpenFOAM) is a free, open-source software for CFD. It has an extensive range of features to solve anything from complex fluid flows involving turbulence and multi-phase  . Its versatile C++ toolbox for the Linux operating system enables developing customized, efficient numerical solvers and pre-/post-processing utilities for all kinds of CFD applications by solving the Navier-Stokes equations. OpenFOAM uses a cell-centered Finite Volume Method (FVM) to solve the partial differential equations of continuum mechanics and fluid flow. In this approach, the equations are integrated over each of the control volumes (cells) on the mesh, and volume integrals that contain a divergence term are converted to surface integrals using Gauss’s theorem. The surface integrals can then be evaluated by summing the contributions from each of the cell faces. This approach to the solution of the equations requires a method to extrapolate the velocity stored at the centroid of the cell to the value of the velocity at the face of the cell. The time integration can be done through Backward Euler, Steady-state solver, Crank-Nicholson. The available gradient, divergence, Laplacian, and interpolation schemes are the second-order central difference, Fourth-order central difference, First order upwind and First/second-order upwind. The available solvers, options in specifying the boundary conditions, mesh generation tools, flow visualization software, and extensive documentation is making OpenFOAM popular among engineering and sciences modeling community. It has a large user base across most areas of engineering and science, from both commercial and academic organizations. For multiphase flows, the accuracy of the solution is dependent on the free surface reconstruction algorithm, and for complex flows, this continues to be a difficult task.
Modeling Turbulence Using OpenFOAM
Flow in a hydraulic jump is turbulent. In turbulent flows, the field properties in the vicinity of the jump are random functions of space and time. A feature of the flows is the presence of small-scale, high-frequency random fluctuation, which is superimposed on the main flow that has a primary flow axis. Although the magnitude of these fluctuations is small, they tend to have a major impact on some of the jump characteristics. In applications, where analyzing the turbulent characteristics of the hydraulic jump are essential, application of three dimensional CFD models is required. OpenFOAM provides a variety of turbulence model options from Reynolds-Averaged Navier-Stokes (RANS) to Large Eddy Simulation (LES) and Direct Numerical Simulation (DNS).
Reynolds averaged Navier-Stokes simulation (RANS) also known as Reynolds averaged simulations have been widely used for solving the time-averaged flow equations. The assumption behind the RANS equations is that the time-dependent turbulent velocity fluctuations in Navier-Stokes equations can be separated from the mean flow velocity. The new set of unknowns called the Reynolds stresses are functions of the velocity fluctuations. Solving these require using a turbulence model. The standard k-ε  , k-ω  and the k-ω Shear Stress Transport (SST) model  have been widely used. The equations in the k-ε model can be written as 
where k is turbulent kinetic energy,
is the dissipation rate of k, t is time,
is the coordinate in the
is dynamic viscosity,
is turbulent dynamic viscosity,
is the production of turbulent kinetic energy,
is the buoyancy effect,
is the dilatation effect, and
are the moduli of mean rate-of-strain tensor.
The above limitation in RANS equations has given rise to the Large Eddy Simulation (LES) method as an alternative to the RANS equations. In LES, the turbulent scales under a specific filter length scale are modeled (RANS approach), while for the larger ones Navier Stokes equations are resolved (DNS approach). The computational costs associated with LES are higher when compared to RANS models. Some of the available LES turbulence models in OpenFOAM include Smagorinsky, kOmegaSSTDES, WALE, DeardorffDiffStress and dynamicKEqn.
In Direct numerical simulation (DNS), the full Navier-Stokes equations are solved by resolving the whole range of spatial and temporal scales of turbulence. Although DNS appears to be the preferred approach, the limiting factor in using this is the expensive computing cost. This is due to the small mesh size that is required for capturing turbulence which occurs at varying spatial scales coupled with the use of higher-order accurate numerical techniques. The dnsFoamsolver facilitates DNS simulations in OpenFOAM.
HEC-River Analysis System (RAS) facilitates one-dimensional steady flow, one and two-dimensional unsteady flow, sediment transport/mobile bed computations, and water temperature/water quality modeling  . The one-dimensional river analysis components are steady flow water surface profile computations and unsteady flow simulation. In steady-state mode, the energy equation between successive cross-sections is solved. Energy losses are evaluated using friction and contraction/expansion coefficients. The momentum equation may be used in situations where the water surface profile is rapidly varied, as in the present case. HEC-RAS has been widely used across various hydraulic simulations, and its results often act as a benchmark for other models.
Water Surface Pressure Gradient (WSPG) model is perhaps the first numerical model that was developed by the Los Angeles County Department of Public Works. It solves the Bernoulli energy equation between any two cross-sections, using the standard step method  . The program computes uniform and non-uniform steady flow water surface profiles. As part of the solution, it can compute the jump characteristics.
As mentioned earlier, we have used the results of Bayon et al.  as the benchmark data for validating the hydraulic models. While we refer readers (for a detailed understanding of the equations, modeling approach, parameters, results) to their work, some salient aspects of their effort that will suffice for our current model comparison effort are summarized below. They developed a three-dimensional computational fluid dynamics model (using OpenFOAM) for analyzing hydraulic jumps in a horizontal smooth rectangular prismatic channel. The Reynolds-averaged Navier-Stokes (RANS) equations were solved. The tested turbulence models were the Standard k-ε, RNG k-ε, and SST k-ω. Sensitivity tests by varying the mesh size, turbulence parameters, and boundary condition location were conducted. Four different mesh sizes ranging from 7:00-8:75 mm, were used. These mesh sizes resulted in total cells ranging across 3.47 to 6.51 million in the computational domain. Since they assumed the channel to be smooth, their model did not consider the channel bottom roughness coefficient. The variables that they studied include sequent depths, efficiency, roller length, free surface profile, and the turbulence model accuracy.
The boundary conditions used in the model are stated below
· OpenFOAM: Bayon et al.  validated their CFD model for a hypothetical channel (6 × 0.5 × 0.75 m) at a flow rate of 0.177 m3/s. At the upstream end, the authors used a flow depth of 0.07 m, and a velocity profile was imposed using a Dirichlet boundary condition. The pressure profile was hydrostatic. The constants in the RANS turbulence modeled were assigned a low value, and a short stretch of the channel was added so that the flow is well developed before the jump forms. At the downstream end, instead of using the subcritical flow depth of 0.553 m, a velocity profile is imposed, so that hydrostatic pressure profile develops. As long as the mass is conserved in the system, this downstream boundary approach will translate to the required subcritical flow depth. The channel bottom was assumed to be smooth. A no-slip condition is imposed at the walls, and roughness is not considered. At the top of the channel surface, an atmospheric boundary condition is imposed which allows fluids to enter and leave the channel. The density and the kinematic viscosity are ρ = 1000 kg/m3 and 10−6 m/s2.
· HEC-RAS: At the upstream and downstream ends of the domain, the water surface elevations of 0.07 m and 0.553 m were specified. The flow in the channel was 0.177 m3/s.
· WSPG: At the downstream end (system outlet), a flow depth of 0.553 m was assigned. At the upstream end (system headwork), a depth of 0.07 m was specified. The flow in the channel was specified as 0.177 m3/s.
Figure 1 is a plot of the stationary depth profile along the length of the channel for the three subject models. In the hydraulic models, the channel length was divided into 100 cross-sections, and a manning’s roughness coefficient of 0.013 was used. While for both the HEC-RAS and WSPG models, the jump forms across two adjacent nodes, the result of the CFD model is more in agreement with the reported experimental data, where the jump forms over an elongated reach. Although all the three turbulence models performed well, RNG k-ε solution was identified as more accurate, and hence we chose it for validating the hydraulic model's output. As shown by Bayon et al.  , the CFD model can capture additional details of the flow that occur along the other two dimensions. These include bubble breakup and coalescence, fluid mixing, free-surface turbulent interactions, surface wave formation and breaking processes. Although the considered hydraulic models cannot capture these details, the location of the jump which is an important variable in design computations is reasonably predicted. The sensitivity of the end solution to the roughness coefficient for the RAS and WSPG modela is shown in Figure 2 and Figure 3.
Figure 1. Predicted stationary depth profiles for subject hydraulic jump with Froude number of 6.125.
Figure 2. Effect of bottom roughness value on HEC-RAS predicted depth profile.
Figure 3. Effect of bottom roughness value on WSPG predicted depth profile.
Computational methods have evolved significantly over the last decade with computer capabilities greatly increasing thus enabling the solution of massive matrix systems to be economically solved. As a result, the class of differential equations becoming commonplace for use in the typical analysis have also greatly evolved from the solution of the classic Bernoulli’s energy equation to now the full Navier-Stokes equations. With this evolution of technology, it is important to view the new technology modeling outcomes with respect to the prior more traditional modeling approach outcomes. Such a comparison between modeling technology levels is provided for the situation of a relatively high Froude number flow with a hydraulic jump. Other such comparisons involving other related topics are important to be examined in order to provide continuity between modeling and advances in technology.
The article compared the predicted steady-state flow profile of a hydraulic jump. The solution from two one-dimensional hydraulic models was compared with a published benchmark outcome, produced from a three dimensional CFD model. The CFD model solved the three-dimensional RANS equations using computational model OpenFOAM. The hydraulic models solved the standard Bernoulli’s energy equation. Based on the outcomes from all three models, it can be concluded that 1) the hydraulic model depth profiles are similar to the CFD outcome 2) the hydraulic models fail to adequately predict the length of the jump and 3) the solution from hydraulic models is sensitive to channel roughness value.