Numerical Simulation of Hydraulic Jumps. Part 2: Recent Results and Future Outlook

Abstract

During the past two decades, hydraulic jumps have been investigated using Computational Fluid Dynamics (CFD). The second part of this two-part study is devoted to the state-of-the-art of the numerical simulation of the hydraulic jump. First, the most widely-used CFD approaches, namely the Reynolds-Averaged Navier–Stokes (RANS), the Large Eddy Simulation (LES), the Direct Numerical Simulation (DNS), the hybrid RANS-LES method Detached Eddy Simulation (DES), as well as the Smoothed Particle Hydrodynamics (SPH), are introduced pointing out their main characteristics also in the context of the best practices for CFD modeling of environmental flows. Second, the literature on numerical simulations of the hydraulic jump is presented and discussed. It was observed that the RANS modeling approach is able to provide accurate results for the mean flow variables, while high-fidelity methods, such as LES and DES, can properly reproduce turbulence quantities of the hydraulic jump. Although computationally very expensive, the first DNS on the hydraulic jump led to important findings about the structure of the hydraulic jump and scale effects. Similarly, application of the Lagrangian meshless SPH method provided interesting results, notwithstanding the lower research activity. At the end, despite the promising results still available, it is expected that with the increase in the computational capabilities, the RANS-based numerical studies of the hydraulic jump will approach the prototype scale problems, which are of great relevance for hydraulic engineers, while the application at this scale of the most advanced tools, such as LES and DNS, is still beyond expectations for the foreseeable future. Knowledge of the uncertainty associated with RANS modeling may allow the careful design of new hydraulic structures through the available CFD tools.

1. Introduction

Numerical modeling has been a revolution in most civil engineering disciplines. With the boom of available computer power and the continuous algorithm improvement, classic Partial Differential Equations (PDE) have been numerically addressed, allowing accurate case-dependent solutions for a wide range of problems. However, fluid dynamics equations remain obscure to both researchers and practitioners, and highly accurate numerical solutions have been obtained uniquely for very simple geometries. Computational Fluid Dynamics (CFD) has benefited from simplifications, which mainly came in the form of turbulence models, hence allowing the study of large, complex problems. Nonetheless, the focus is on the accuracy for complex hydraulic engineering challenges. For such problems, modeler needs to be aware of the associated uncertainties and take them into consideration in the design stage. The accuracy of CFD for such specific problems can deviate from typical values obtained in simpler flows, reported in the classic literature (e.g., [1,2]).

Common hydraulic problems usually involve shear regions, alone or interacting with each other, free surfaces, and under certain conditions, air entrainment. All these phenomena are observed in hydraulic jumps (Figure 1). Traditionally, hydraulic jumps have been widely studied experimentally, as depicted in Part 1 of this study. As proposed in Part 1 [3], hydraulic jumps arise as an adequate workbench for numerical modeling performance with both mean and turbulent variables largely reported in the literature. While numerical modeling allows easy access to the complete flow field, deviations from reality may occur (note the aeration differences observable in Figure 1). Systematic reporting of numerical modeling accuracy may allow the study of more complex problems with a certain level of confidence (e.g., stilling basins; see Figure 2).

Part 2 of this study presents several three-dimensional numerical modeling techniques and their achievements to date. In Section 2, focus is put on giving the basic concepts, while the adventurous reader is addressed to the relevant experts’ reviews and books for each of the techniques [8,9,10,11,12,13]. Some comments on best practices are presented, as well. Section 3 brings to light a large variety of studies conducted during the past twenty years, disclosing aspects related to the expected accuracy of the employed numerical modeling methods. Conclusions are presented in Section 4, followed by Section 5, which includes a brief future outlook.

2. Numerical Methods for Hydraulic Jump Research

2.1. Overview

Incompressible two-phase flows are governed by the equations of conservation of mass and momentum, defined for both air and water phases. These laws are represented through the Navier–Stokes (N-S) equations, which, in their original form, encompass all known internal and external effects of the motion of a fluid. Further details can be found in general fluid mechanics books such as Pope [8], among others. Any form of numerical modeling of hydraulic jumps was achieved by means of the N-S equations or any of their simplifications. N-S equations are inherently nonlinear, time-dependent, and three-dimensional PDE. Three-dimensionality is essential in these equations, because vortex stretching is not possible in a two-dimensional space.

The numerical modeling of turbulent flows involves considerable difficulties, related to the accurate computation of the vortex dynamics, the determination of the instantaneous position and configuration of the free surface, turbulence effects, and air entrainment [13,14]. The hydraulic jump has been studied using mathematical models based on the N-S equations, free surface tracking methods such as the Volume Of Fluid (VOF) method [15], or the level-set method and different turbulence models (see, e.g., [16,17,18]).

The VOF method uses a function $f(x,y,z,t)$ to assign the free surface. The function f, representing a fraction of fluid, is obtained from the following equation:

$$\frac{\partial f}{\partial t}+U_j\frac{\partial f}{\partial x_j}=0$$

(1)

U being the convection velocity.

A unit value of f corresponds to a cell full of fluid, while a zero value indicates that the cell contains no fluid. Cells with f values between zero and one contain a free surface. The function f can be solved with different methods, and different schemes have been implemented to improve the accuracy of the free surface reconstruction. More complicated free surface tracking schemes are also possible; see for instance the methods used in the study of Mortazavi et al. [19], Fuster et al. [20], or the extensive review of Scardovelli and Zaleski [21]. It must be noted that the free surface is a movable interface where both density and viscosity change abruptly, being prone to significant numerical diffusion. Accurate estimation of surface tension effects directly depends on proper representation of free surface curvature.

2.2. One- and Two-Dimensional Modeling

One- and two-dimensional modeling was attempted for hydraulic jumps, obtaining great insight into the flow structure. Non-hydrostatic approaches are noteworthy, providing considerable understanding of the free surface profile, velocity fields, and even characteristic frequencies. Some one- and two-dimensional modeling includes Madsen and Svendsen [22], Valiani [23], Madsen et al. [24], Castro-Orgaz and Hager [25], and Richard and Gavrilyuk [26]. For the sake of briefness, the interested reader is addressed to the extensive book of Castro-Orgaz and Hager [27] for a thorough description on non-hydrostatic methods.

2.3. Three-Dimensional Modeling

2.3.1. RANS Approach

The most commonly-applied approach to simulate a turbulent flow is that based on the time averaging, herein Reynolds-averaging, of the N-S equations, leading to the Reynolds Averaged N-S (RANS) equations. In this statistical approach, the instantaneous values of velocity and pressure are depicted as the sum of a mean and fluctuating component [28]. Any instantaneous velocity $u_i$ at a given coordinate $\mathbf{x}$ and time t can be written as:

$$u_i(\mathbf{x},t)=\overline{u_i}(\mathbf{x},t)+u_i^{\prime }(\mathbf{x},t)$$

(2)

$\overline{u_i}$ and $u_i^{\prime }$ being the mean and fluctuating velocity, respectively.

Note that the mean velocity is still meant to depend on the time t, which is possible when the time averaging is defined as [2]:

$$\overline{u_i}(\mathbf{x},t)=\frac{1}{T}{\int }_t^{t+T}{u}_i(\mathbf{x},t)\mathrm{d}t,T_1\ll T\ll T_2$$

(3)

$T_1$ is a characteristic time-scale of the turbulent fluctuations and $T_2$ a characteristic time-scale of the unsteadiness. While $T\ll T_2$, some flow unsteadiness can be reproduced (e.g., jump toe fluctuations, as shown by Bayon et al. [5]). Hence, we are assuming that a difference between both time-scales exists and is of several orders of magnitude. Consequently, the time derivative of the RANS equations does not necessarily vanish, and transient flows can be also modeled through a time-averaged procedure:

$$\rho \left(\frac{\partial \overline{u_i}}{\partial t}+\overline{u_j}\frac{\partial \overline{u_i}}{\partial x_j}\right)=-\frac{\partial \overline{p}}{\partial x_i}+\frac{\partial }{\partial x_i}\left(\mu \frac{\partial \overline{u_i}}{\partial x_i}+R_{ij}\right)$$

(4)

$\mu$ is the dynamic viscosity, $\rho$ the density, and $\overline{p}$ the mean pressure. Due to the non-linearity of the N-S equations, averaging leads to unknown correlations between scalar quantities known as Reynolds stresses, which are represented in the previous formula with the Reynolds stress tensor $R_{ij}=-\rho \overline{u_i^{\prime }{u}_j^{\prime }}$ [28]. The Reynolds stresses are six additional unknowns introduced by the averaging procedure, which, nonetheless, do not introduce any new equation. This is the well-known closure problem of the RANS equations [2,8]. How the Reynolds stresses are modeled has led to a wide, eclectic range of turbulence models.

The simplest RANS models are the eddy viscosity models (based on the Boussinesq hypothesis). Following the Boussinesq hypothesis, Reynolds stresses are modeled using an eddy (or turbulent) viscosity ${\nu }_t$ (equivalently, ${\mu }_t$), which, ultimately, means that the effect of turbulence is to act on the mean flow as an increased viscosity [29]. Hence, apparent stresses $R_{ij}$ are nulled, and an increased viscosity is contemplated instead. Based on dimensional analysis, ${\nu }_t$ can be determined from a turbulent time-scale (or velocity scale) and a turbulent length-scale. Note that two parameters are needed to close the system properly through an eddy viscosity model. The most commonly-used variables are listed in

Table 1. Variables used in the most common turbulence models to calculate ${\mu }_t$.

Dimensional Scale	Unit	Equation
Turbulent kinetic energy	L${}^2$/T${}^2$	$k=\overline{u_i^{\prime }{u}_i^{\prime }}/2$
Turbulence dissipation rate	L${}^2$/T${}^3$	$ϵ=\nu \overline{\frac{\partial u_i^{\prime }}{\partial x_j}\frac{\partial u_i^{\prime }}{\partial x_j}}$
Specific dissipation rate	1/T	$\omega =ϵ/k$

.

Different turbulence models estimate ${\nu }_t$ based on different turbulent variables. Some examples are as follows:

• The Spalart–Allmaras model [30] solves a transport equation for an auxiliary turbulent viscosity ${\nu }_t\sim \tilde{\nu }$.
• The standard $k-ϵ$ [31,32] and the RNG $k-ϵ$ [33,34] solve transport equations for k and another for $ϵ\to {\nu }_t\sim k^2/ϵ$. The former is the most widely-used engineering turbulence model for industrial applications because of its robustness and reasonable accuracy for a wide range of flows. It is a two-equation transport model for k and $ϵ$, their coefficients being empirically derived. The RNG version adds a term for the $ϵ$ equation, which is known to be responsible for differences in its performance [8]. It is oftentimes mentioned that accuracy for rapidly-swirling flows is improved, despite fewer validations having been conducted [1].
• Standard $k-\omega$ (2008) and SST $k-\omega$ [35] solve transport equations for k and $\omega \to {\nu }_t\sim k/\omega$. The $k-\omega$ model is a two-equation transport model based on solving equations for k and $\omega$. The original version was known to be excessively sensitive to free stream quantities, which is an undesired feature [36], while the revised version of Wilcox [37] attempts to solve this issue. The SST $k-\omega$ is a variant of the former. The SST $k-\omega$ model uses a blending function to transition gradually from the standard $k-\omega$ model near the wall to a high Reynolds number version of the $k-ϵ$ model in the outer portion of the boundary layer [35].

In addition to the eddy viscosity models, the Reynolds stress tensor can be modeled (based on some assumptions) by formulating an equation for each of the terms. The large dimension of the resulting equations and the small gain in accuracy has led to a reduced use in comparison to eddy viscosity models [38].

2.3.2. LES and DES Approach

Large Eddy Simulation (LES) has been developed considering that the large scale is strongly energetic, whereas the small scale, which is more dissipative, is universal and represents only a small part of the turbulent spectrum. Hence, the large scale of turbulence is resolved, while the small scale is modeled [39]. The main difference between RANS and LES approaches is that the N-S equations are, in the former, time-averaged and, in the latter, space filtered. Notwithstanding that, both sets of equations take a similar form because a stress tensor is created by the time-averaging and filtering processes. In the LES approach, the stress tensor resulting from low-pass filtering is:

$$R_{ij}^{SGS}=\rho \left(\overline{u_i{u}_j}-\overline{u_i}\overline{u_j}\right)$$

(5)

These sub-grid scale stresses, differently than in RANS models, attempt to approximate only the part of the turbulence that is unresolved, i.e., below the size of the low-pass filter window. It is noteworthy that while some unsteadiness can be reproduced by RANS models, their computation does not change (in theory) with the cell size, while in LES, $R_{ij}^{SGS}\to 0$ and, consequently, converges to a DNS. In a similar manner to RANS models, the stresses appearing in the flow equations can be modeled in different ways [39].

Some hybrid LES-RANS approaches were also proposed in the literature. They have the general goal of combining the advantages of LES and RANS, minimizing their shortcomings. For the Detached Eddy Simulation (DES) approach [40], the near-wall region is simulated using the RANS method, while away from the wall, LES is used [39].

2.3.3. DNS Approach

DNS numerically addresses the N-S and continuity equations with physically consistent accuracy in space and time, in such a way as to resolve all the essential turbulent scales. If the mesh is fine enough, the time step is short enough, and the numerical scheme is designed to minimize the dispersion and dissipation errors, one obtains an accurate three-dimensional time-dependent solution of the governing equations, in which the only errors are those introduced by the residual approximations incorporated in the numerical scheme and in the number-representation technology of the computing machine [41]. The contributions of DNS to turbulence research in the last few decades have been impressive [42]. As one can see, the biggest problem regarding DNS is their overwhelming requirement for computer power in a sense of both the processor’s performance and the size of the memory for storing intermediate results.

2.4. Best Practices for CFD Modeling of Environmental Flows

The application of CFD methods to the environmental fluid mechanics area generally involves many different issues related to the existing uncertainty of geometry and boundary conditions, drag coefficients, driving forces, and the interactions among different processes and inputs [43]. This is particularly important in the numerical simulation of the hydraulic jump, which involves fluctuating boundaries, as well as a multiphase flow. These uncertainties involve almost every aspect of the modeling process, and it may therefore be very difficult to assess model performance, as these uncertainties also apply to the experimental model (see Part 1). This implies that CFD application to environmental flows may have research priorities very different from other CFD applications, bearing in mind that the closer they are to the left side of

Table 2. Turbulence modeling, indicator for resolved/modeled turbulence (R: Resolved; M: Modeled), accuracy, and computational cost of RANS, LES, Detached Eddy Simulation (DES), and DNS.

	RANS	DES	LES	DNS
Modeling approach	Reynolds-averaged	RANS-LES coupling	Sub-grid scale	Not necessary
Large scales/inertial Subrange/dissipation scale	M/M/M	R/M/M	R/R/M	R/R/R
Expected accuracy	Mean variables	Allow inertial’ modeling	Allow inertial’ modeling	Up to the Kolmogorov scale
Computational cost	Low	High	High	Extremely high

, the more necessary proper validation will be.

Roache [44] discussed the quantification of uncertainty in CFD methods, pointing out the difference between verification, which means solving the equations right, and validation, which is solving the right equations, concluding that a code cannot be validated, but only a calculation or a range of calculations with a code can be validated. Later, Roache [45] noted that calibration, which is the adjustment or tuning of free parameters in a model to fit the model output with experimental data, which is a sometimes necessary component of model development, cannot be considered as validation, unless when the previously-calibrated model predictions are evaluated using a dataset not used in the tuning.

Some criteria for the application of CFD codes to open channel flows were proposed by Knight et al. [46], while Lane et al. [47] expanded the guidelines specified by ASME [48] for the specific case of open channel flows. They demonstrated that the ASME criteria may not be sufficient for many practical fluvial applications, but can provide a minimum framework. In particular, Lane et al. [47] argued that the ASME criterion on model validation should be replaced by a more comprehensive framework based on sensitivity analysis, including mesh sensitivity study, benchmarking, and flow visualization. For mesh sensitivity analysis, methods such as that of Celik et al. [49], recommended by the ASME, should be preferred. Celik et al.’s work [49] was based on the Richardson extrapolation and allows computation of the numerical uncertainty of a given variable. A common error in numerical studies is the misjudgment in detecting the relevant variable to be subject to the numerical uncertainty analysis. Preferably, the numerical uncertainty analysis should be performed on all the results presented or on variables with an expected similar level of uncertainty. A study should not assess the numerical uncertainty of the water depth and claim that turbulence fluctuations would present a similar level of mesh dependence. If the aim of the study is to reproduce aeration, air concentrations or turbulence quantities should be the main variables of the sensitivity analysis.

Finally, as suggested by Blocken and Gualtieri [43], it is advisable that future best practice guidelines for CFD applications to environmental flows include a larger application of benchmarking because this provides relevant information for reliable comparison among different CFD techniques and approaches. Therefore, specific benchmarks should be developed in the different contexts and practical applications of CFD methods for environmental flows.

3. Hydraulic Jump Numerical Modeling Studies

3.1. General Comments

The study of aerated flows in hydraulic engineering is still challenging in terms of either experimental investigation or numerical modeling considering the large number of relevant equations, parameters, and their complexity (see Part 1). Thus, the hydraulics of two-phase flows can greatly benefit from the insights provided by numerical simulations, and the last few decades have seen the development of powerful numerical capabilities with direct applications into gas-liquid flows [11].

Depending on the flow field description, two approaches can be distinguished: Lagrangian methods and Eulerian methods. The first approach follows a fluid particle at each point. The latter treats the fluid properties as a function of time and space. Both approaches have been used in fluid mechanics to describe hydraulic jump flows. Among the meshless Lagrangian techniques, one of the most versatile techniques is Smoothed Particle Hydrodynamics (SPH), which solves flow equations for a set of moving particles (with a certain mass). Use of a kernel function to interpolate flow variables becomes critical, but the determination of the free surface is natural, unlike grid-based techniques (Eulerian methods), which need an additional method to track the free surface (i.e., VOF by Hirt and Nichols [15]).

Within Eulerian methods, the RANS approach is the most widely used, followed by LES. Hybrid RANS-LES methods, such as DES, are growing in popularity. For the sake of briefness, the aforementioned methods are not extensively described in this study, but the reader is addressed to previous works, as referenced in

Table 3. Milestone references on different numerical modeling approaches.

Flow Field Description	Approach	References
Lagrangian	Smoothed Particle Hydrodynamics	Gingold and Monhagan [50]
		Monaghan [51]
		Violeau [52]
		Violeau and Rogers [53]
Eulerian	Reynolds-Averaged Navier–Stokes	Reynolds [54]
		Rodi [28]
		Spalart [36]
		Wilcox [2]
	Detached Eddy Simulation Large Eddy Simulation	Spalart [55]
		Labourasse et al. [56]
		Rodi et al. [39]
	Direct Numerical Simulation	Moin and Mahesh [42]
		Prosperetti and Tryggvason [11]
		Alfonsi [41]

and

Table 4. Most widely-used turbulence models in RANS applications and original contributions.

Turbulence Model	References
Standard $k-ϵ$	Jones and Launder [31]
	Launder and Sharma [32]
RNG $k-ϵ$	Yakhot and Orszag [33]
	Yakhot et al. [34]
$k-\omega$	Wilcox [2]
	Wilcox [37]
SST $k-\omega$	Menter [35]

.

Table 5. Numerical simulation studies of hydraulic jumps and inlet Froude number (${\mathrm{F}}_1$).

Flow Field Description	References	Year	Numerical Approach	Turbulence Model	Resolution Dependence Study	${{F}}_1$
Lagrangian Methods	Lopez et al. [57]	2010	SPH	$k-ϵ$	Meshless method	3.25, 3.41, 3.58, 3.61, 3.64, 4.12, 4.45, 4.88, 5.13, 5.24, 6.62, 7.06, 7.16
	De Padova et al. [58]	2013	XSPH	mixing length and $k-ϵ$	Meshless method	3.90, 8.30
	De Padova et al. [59]	2018	XSPH	mixing length and $k-ϵ$	Meshless method	2.58, 2.66, 2.70, 3.89, 4.99, 5.11
Eulerian Methods	Chippada et al. [60]	1994	RANS	STD $k-ϵ$	−	2.00, 4.00
	Zhao et al. [61]	2004	RANS	STD $k-ϵ$, $k-l$	−	1.46
	Gonzalez and Bombardelli [62]	2005	RANS	STD $k-ϵ$	√	2.00, 2.50, 3.32
	Carvalho et al. [14]	2008	RANS	RNG $k-ϵ$	−	6.00
	Abbaspour et al. [63]	2009	RANS	STD $k-ϵ$ RNG $k-ϵ$	−	4.00, 4.70, 5.00, 5.70, 5.80, 6.10, 7.00, 7.20, 8.00
	Ma et al. [64]	2011	RANS	SST $k-\omega$	−	1.98
	Ebrahimi et al. [65]	2013	RANS	STD $k-ϵ$	−	3.00, 3.30, 3.60, 3.60, 5.00, 5.70, 6.70, 8.00
	Bayon-Barrachina and [66]	2015	RANS RANS	STD $k-ϵ$ RNG $k-ϵ$ SST $k-\omega$	√	6.10
	Witt et al. [67]	2015	RANS	realizable $k-ϵ$	√	2.43, 3.65, 4.82
	Bayon et al. [5]	2016	RANS	RNG $k-ϵ$	√	6.50
	Witt et al. [68]	2018	RANS	realizable $k-ϵ$	√	2.43, 3.65, 4.82
	Harada and Li [69]	2018	RANS	$k-ϵ$, $k-\omega$	√	5.80
	Valero et al. [6]	2018	RANS	RNG $k-ϵ$, $k-ϵ$	√	3.12, 3.88, 4.20, 6.17, 6.37, 6.47, 8.27, 9.52
	Ma et al. [64]	2011	DES		√	1.98
	Jesudhas et al. [70]	2018	DES		√	8.5
	Gonzalez and Bombardelli [62]	2005	LES		√	2.00, 2.50, 3.32
	Lubin et al. [71]	2009	LES		−	5.09
	Mortazavi et al. [19]	2016	DNS	Full solution of all the scales	√	2.00

lists several numerical CFD studies dealing with hydraulic jumps. They are presented chronologically and are classified according to the flow field description and the type of numerical approach used. RANS-based methods prevail due to their lower computational cost, but, in turn, turbulence model dependence is stronger (we mention: standard $k-ϵ$, termed in Table 5 as STD $k-ϵ$, RNG $k-ϵ$, and SST $k-\omega$).

3.2. Reynolds-Averaged Navier–Stokes Approach

The early study of Chippada et al. [60] studied two Froude numbers (${\mathrm{F}}_1$ = 2, 4) using the standard standard $k-ϵ$ model and 2D simulations, together with a Eulerian-Lagrangian flow description helping the determination of the moving free surface.

Zhao and Misra [61] focused on the mean flow motions of a hydraulic jump using the numerical model RIPPLE [72]. Two turbulence models were used: the standard $k-ϵ$ model and the $k-l$ model of Zhao et al. [61]. Zhao and Misra focused on weak hydraulic jumps, their main interest being on bores. Model results were compared to the LDV laboratory measurements of Bakunin [73] and Svendsen et al. [74] with an overall adequate agreement, despite underestimation of downstream streamwise velocities.

Gonzalez and Bombardelli [62] simulated the entire two-phase flow of low Froude number hydraulic jumps (mean flow, turbulence, and air entrainment). Air features were incorporated through the free surface. Gonzalez and Bombardelli [62] compared 2D and 3D RANS and LES models to the experimental observations of Liu et al. [75], including both mean flow and turbulence, and some qualitative analysis of the air fractions.

Carvalho et al. [14] carried out a 2D RANS simulation of a hydraulic jump with a higher Froude (${\mathrm{F}}_1$ = 6). The jump was contained in a stilling basin downstream of a spillway, and the comparison was conducted against their own experimental data and Hager [76] empirical relations. The numerical model combined high-order schemes [77] for the advection terms, a refined VOF algorithm [78] for free surface tracking, the Fractional Area-Volume Obstacle Representation (FAVOR) method for representing internal obstacles (as defined by Hirt and Sicilian [79]), and the RNG $k-ϵ$ turbulence model. The numerical models of Carvalho et al. [14] were not able to properly represent either the free surface or the pressures. Additionally, shear stresses were underestimated. Nevertheless, the mean velocities were satisfactorily predicted, except in the zone where vortices develop. The maximum values of the velocity and pressure in the flow domain remained below those measured in the experimental installation, but their locations were well predicted.

Abbaspour et al. [63] used 2D numerical models to study the hydraulic jump over a rough bed using both standard and RNG $k-ϵ$ models, together with a VOF method to determine the free surface location. Characteristics of the jump, such as water surface profile, length of the jump, velocity profile in different sections, and bed shear stress were evaluated for a range of Froude numbers 4 < ${\mathrm{F}}_1$ < 8 through a total of 12 simulations. The authors found a close agreement with both the experimental values obtained by Ead and Rajaratnam [80] in terms of mean relative error: about 1%–8.6% for free surface profiles and about 1%–12% for the length of the jump. The bed shear stresses were computed using the measured depths up- and down-stream of the jump. Differences between total shear stress in experimental and numerical models remained small. Abbaspour et al. [63] concluded that the best agreement was obtained using the RNG model.

Void fraction predictions in a hydraulic jump were the main goal of Ma et al. [64], who used a sub-grid air entrainment model coupled with RANS and DES models to study an ${\mathrm{F}}_1$ = 1.98 hydraulic jump. Sub-grid models arise from the hypothesis that the numerical model alone will not be able to reproduce aeration. Thus, some additional “physics” are input; likewise, a turbulence model aims to represent the small scale fluctuations. Free surface was modeled using a level-set method and turbulent viscosity evaluated using the SST model. Results were compared to the experimental model of Murzyn et al. [81], concluding that air transported in the shear layer was better reproduced than in the upper region due to RANS modeling limitations to capture free surface fluctuations.

Ebrahimi et al. [65] carried out 2D simulations (together with the standard $k-ϵ$ and VOF) of hydraulic jumps on a rough bed for a wide range of Froude numbers (3 < ${\mathrm{F}}_1$ < 7) across a total of 16 simulations. The authors found a close agreement with the experimental measurements from Elsebaie and Shabayek [82], with mean relative errors around 0–4.4% for water surface profiles and 0–6.7% for the length of the jump. The model also provided good estimates of the total shear force.

Bayon-Barrachina and Jimenéz [66] used the open source OpenFOAM code with standard $k-ϵ$, RNG $k-ϵ$, and SST $k-\omega$ models to study a 3D classic hydraulic jump with ${\mathrm{F}}_1$ = 6.1. Results were compared to the experimental relations of Hager [76] and Wu and Rajaratnam [83], among others, obtaining high accuracies. The authors found that the RNG $k-ϵ$ model was performing better, followed by the SST $k-\omega$ and the $k-ϵ$; although in all cases, errors were lower than 4%.

Witt et al. [67] simulated the air entrainment characteristics of three Froude number hydraulic jumps using the realizable $k-ϵ$ model, with VOF treatment for the free surface through 2D and 3D simulations. Velocity profiles, void fraction profiles, and Sauter mean diameter were compared to the experimental data from Murzyn et al. [81], Liu et al. [75], and Lin et al. [84], among others. Velocity profiles showed good agreement with the experimental observations for jumps. Void fractions were well predicted with differences rarely over 10%. Witt et al. [67] observed that at least eight cells per bubble diameter are necessary to reproduce the associated aeration. Furthermore, 3D simulation improved void fraction predictions in the shear layer, and consequently, the average free surface elevation.

A systematic performance analysis of the OpenFOAM and FLOW-3D^® codes was presented by Bayon et al. [5], investigating a steady hydraulic jump of a Froude number of 6.5. The RNG $k-ϵ$ model and the VOF method were used. Performance assessment was conducted by comparing several hydraulic jump variables with the literature and their own experimental data [76,85,86,87,88,89,90]. Numerical uncertainty was assessed following Celik et al. [49], finding that FLOW-3D^® presented a faster convergence with coarser meshes, whereas OpenFOAM had a more stable solution with refinement. All mean variables studied (not including aeration) obtained accuracies above 90%, except for the roller length, which remained at around 80%, and the negative velocities at the upper jump region. An autocorrelation analysis was conducted, finding good computation of jump toe frequencies, thus showing that RANS models can also reproduce the unsteadiness typical of hydraulic jumps.

Witt et al. [68] extended the analysis of the work of Witt et al. [67], studying some air-water mean and turbulent flow features, showing that satisfactory prediction of turbulence quantities can also be achieved through RANS modeling.

Harada and Li [69] modeled a 2D hydraulic jump of ${\mathrm{F}}_1$ = 5.8 using the standard $k-ϵ$ model and the VOF method. Air concentrations were compared to the experimental model of Kucukali and Chanson [91], obtaining errors below 10% roller length against Hager [76].

Valero et al. [6] studied the flow field (RNG $k-ϵ$ model and the VOF method) inside a United States Bureau of Reclamation (USBR) Type III stilling basin and compared numerical results, contemplating eight Froude numbers (in the range of 3.12–9.52), to the data of Frizell and Svoboda [7]. Good agreement was obtained for sequent depth relations and sweep-off point. The numerical model allowed insight into the distribution of forces on the different elements of the basin, thus explaining the experimental observation of Frizell and Svoboda [7] on the capacity of the Type III to hold inside jumps even for very low tailwater levels.

Figure 3 shows the accuracy of previous RANS models with standard $k-ϵ$, RNG $k-ϵ$, and the SST $k-\omega$ turbulence models. In Figure 3, the results are also classified by Froude number; following the jump typology of Bradley and Peterka [92] (see Part 1). The RNG $k-ϵ$ model comparatively provided the best performances for all four jump parameters, despite is also showing the largest scatter. The standard $k-ϵ$ model provided more accurate results in terms of the roller length and the free surface.

3.3. Detached Eddy Simulation and Large Eddy Simulation

Ma et al. [64], besides RANS simulations, used also the DES approach to predict air entrainment in a hydraulic jump of ${\mathrm{F}}_1$ = 1.98. By comparing the numerical results with the experimental data of Murzyn et al. [81], Ma et al. [64] found that the averaged DES results predicted the observed void fraction profiles in both the lower shear layer and the upper roller region, while the RANS approach missed the latter. This is due to the better capabilities of DES to capture strong fluctuations. Nonetheless, Bayon et al. [5] showed that the RANS model can reproduce satisfactorily big-scale fluctuations inside hydraulic jumps.

Recently, Jesudhas et al. [70] studied an ${\mathrm{F}}_1$ = 8.5 hydraulic jump using a DES approach together with a VOF and a high-resolution capturing scheme. Not only were the numerical results compared to past experimental literature (among others, [84,89,93,94])—with excellent agreement—but also new insights into the turbulence structure were presented, as for instance how the developing shear layer leads to intense free surface deformations. The study of Jesudhas et al. [70] proved that DES may suffice to study turbulent features.

Concerning the LES part of the study of Gonzalez and Bombardelli [62], satisfactory agreement was found in terms of mean flow, turbulence, and air entrainment. However, LES outcomes showed similar patterns in the time-averaged flow field as those obtained with the $k-ϵ$ model.

Large eddy simulation of the hydraulic jump was also executed by Lubin et al. [71]. The numerical results were then compared with experimental data collected on a physical model. Important handicaps were noticed as a result of the strong interface deformations, break-ups, and high shear levels. The flow did not tend to a stationary state because of some numerical diffusion at the free surface.

3.4. Direct Numerical Simulation

Mortazavi et al. [19] recently presented the first Direct Numerical Simulation (DNS) of a stationary turbulent hydraulic jump with ${\mathrm{F}}_1$ = 2 and a low Reynolds number (11,000) with no boundary layer development at the inlet. A non-dissipative geometric VOF method was used to track the strongly-distorted interface, while level-set equations were also solved in order to calculate the interface curvature for the surface tension forces’ computation. They compared the numerical results with the experimental data by Murzyn et al. [95]. It was shown that the Reynolds number did not influence the DNS statistics. The error of the mean velocity, fluctuations, and Reynolds stresses computed was within ±10% accuracy. Reasonably good agreement was observed with experiments for void fraction measurements and interface length-scale measurements. Mortazavi et al. [19] found that bubble entrainment presented a periodic nature, bubbles being generated in patches with a specific frequency associated with the roll-up frequency of the roller at the toe of the jump. This highlights the relevance of experimental studies focused on hydraulic jump frequencies to better understand the air entrainment process. The study of Mortazavi et al. [19] also shed light on the energy transfer processes, yielding better understanding of the energy dissipating properties of the jump. Additionally, Mortazavi et al. [19] analyzed how the effects of scale affect the entrainment of air by modeling the same hydraulic jump with the same flow parameters, except the Reynolds number (5000 and 11,000). Despite the limited Reynolds number range, this type of information can be priceless for the study of scale effects.

3.5. Smoothed Particle Hydrodynamics

Lopez et al. [57] investigated the capabilities of Smoothed Particle Hydrodynamics (SPH) to reproduce a mobile hydraulic jump for different Froude numbers. Numerical simulations were compared to their own physical model. Results showed good agreement for ${\mathrm{F}}_1$ < 5. Applying the standard $k-ϵ$ model, a considerable improvement in the model simulating hydraulic jumps with ${\mathrm{F}}_1$ > 5 was also obtained. The main disadvantage of this approach is the doubling of the computational time. Finally, it was found that SPH provides correct estimates of the average pressures at the boundaries, but exhibits large dispersion for instantaneous water depths.

De Padova et al. [58,59] carried out different studies on the numerical modeling of the hydraulic jump using the SPH method and a weakly-compressible XSPHscheme, which includes a mixing-length turbulence model and a two-equation turbulence model to represent turbulent stresses. First, they reproduced the formation of different undular jumps [58]. Instantaneous and time-averaged flow variables were compared with acoustic Doppler velocity measurements for ${\mathrm{F}}_1$ equal to 3.9 and 8.3. The SPH model slightly overestimated the measurements of the velocity, in the zone of the lateral shock wave and in the vortex zone near the lateral wall. The predicted free surface elevations and velocity profiles showed a satisfactory agreement with measurements and most of the particular features of the flow, such as the trapezoidal shape of the wave front and the flow separations at the toe of the oblique shock wave along the side walls. The numerical results were compared with their own experimental investigations. The comparison proved that, when simulating hydraulic jumps without a surface roller, a simple turbulence model based on the mixing-length hypothesis suffices to yield results that, in terms of water depths and average velocity predictions, were like those obtained using the standard $k-ϵ$ turbulence model.

De Padova et al. [59] investigated the oscillating characteristics and cyclic mechanisms in hydraulic jumps. Both the mixing-length turbulence model and the standard $k-ϵ$ model were able to predict some features of the flow, but comparison with experimental results showed that the mixing-length results were closer to the experimental data [96] than those from the $k-ϵ$ in, for instance, predicting the spectra of the surface elevations upstream and downstream of the jump.

4. Conclusions

Hydraulic jumps have been traditionally studied mainly by means of experimental modeling, with some limited contributions from analytical studies (see Part 1 of this study [3]). In the past twenty years, numerical research activity has experienced a large burst in terms of quantity and quality. Two different types of studies can be found: those that aim to validate the numerical tool for a set of flow conditions (mainly RANS modeling techniques) and those that aim to shed light on some unresolved flow features (DES, LES, and DNS).

Surprisingly, RANS modeling approaches have yielded accurate results, with accuracies over 90% for mean flow variables, including air concentrations in some cases. Besides, some fluctuations have been observed, contrary to the early RANS studies’ conclusions. It must be noted, however, that aeration has been studied mainly in terms of air concentrations. Two different solutions can present the same air concentration profile with a completely different distribution of bubble sizes. Future validations in terms of air entrainment should address this key issue. Nonetheless, the number of available validations is limited (Figure 3) with very sparse approaches used, and hence, it cannot be stated that a single model has stood out over the rest. From an engineering point of view, any of the two-equation transport models presented could be used for design purposes as long as the uncertainty reported by previous studies is kept in mind on the final results. In this regard, it must be noted that the recirculation length should hold a larger uncertainty than the variables reported in Figure 3.

High-fidelity approaches, such as DES and LES, have shown that turbulence quantities can be properly reproduced with these models. It is noteworthy that the first DNS study on hydraulic jumps has been released [19], leading to important conclusions in terms of hydraulic jump structure and scale effects, despite the low Froude and Reynolds numbers.

5. Future Outlook

With recent years, numerical studies have become more physically sound and follow more strict analysis procedures. Upcoming model validations should include turbulent variables, and larger, community-based benchmarking may help to cover the gaps of Figure 3. Once having properly established the deficiencies of low fidelity models, they should be compensated with improved turbulence and sub-grid air entrainment models.

Results of the high fidelity models can compensate handicaps of multiphase flow instrumentation. Yet, it is still to be exploited how the information of higher accuracy models (and more computationally expensive, in turn comparable to experimental modeling) is used to formulate more accurate—physically-based—sub-grid and turbulence models to be used in more affordable techniques (e.g., within RANS models or improved DES models).

Solving the shortcomings of less demanding models could anticipate a new era of massive 3D numerical modeling of prototype-scale problems (with scales reaching hundreds of meters and velocities a few dozens of meters per second); which is of great relevance for practitioners, for which LES (and obviously DNS) could remain out of their hardware capabilities for some decades.