Distribution of mean flow and turbulence statistics in plunge pools Free

There is growing consensus regarding the fact that climate change will lead to enhanced extreme flooding in certain areas of the world. This situation must be confronted with spillway capacity and special operational scenarios for large dams. If the capacity of the spillway is insufficient, the dam might be overtopped, thus generating new loading scenarios, and raising questions about potential risky hydrodynamic actions and scour downstream of the dam (Wahl et al. 2008).

Spain is the fourth country in the world according to the number of large dams – it has over 1,200. Fifty percent of these dams were built fifty years ago. In this sense, numerous dams need to be re-evaluated in their safety with respect to potential overflow, in line with what is being done in the USA (FEMA 2014).

When the rectangular jet or nappe flow occurs due to overtopping, the design considerations need to ensure that most energy is dissipated, and that there is minimal to no erosion downstream of the dam. In other words, we need to estimate the hydrodynamic actions on the bottom of the basin where the jet discharges, as a function of the characteristics of the jet (Annandale 2006).

When the jet falls through a long-enough distance, the jet becomes fully developed (see L_b in Figure 1). Castillo et al. (2015) established different equations to calculate the jet energy dissipation in the air and in the water cushion, as a function of the Y/B_j and H/L_b ratios (where Y and H denote the depth of the water cushion at the exit and the total head, respectively, and L_b is the break-up length). During the falling, the energy dissipation is due to the air entrainment into the falling jet and the depth of water upstream of the jet. Energy dissipation in the basin by diffusion effects can only be produced if there is an effective water cushion (Y/B_j> 5.5 for the rectangular jet case (Castillo et al. 2015)).

Within the plunge pool downstream of the impingement point, the flow resembles a flow in a submerged hydraulic jump and a wall jet. However, the situation is complicated here by the air entrainment. Several formulas have been put forward to express the horizontal velocity distribution in the vertical direction. We revisit some of these formulas later in the paper.

There are only a few well-documented references on the numerical simulations of free overflow spillways. Ho & Riddette (2010) analysed the different applications of computational fluid dynamics (CFD) code to hydraulic structures, and identified some limitations that do not allow a comprehensive analysis of the two-phase phenomena. They suggested the following principal lines of future research:

• Validation of the air entrainment along chutes and free-falling jets.

• Extension of models to two-fluid theories to include air entrainment and jet break-up. Validate impact pressures against experimental data.

• Develop parametric studies with turbulence models to identify the level of reliability of the computed pressure.

Castillo et al. (2014, 2015) and Carrillo (2014) have followed the above suggestions and, specifically, they have developed laboratory research on the velocity, air, and pressure fields in the jet and the basin. They undertook experimental observations of jet velocities, and air concentrations with an optical fibre probe, and of pressures in the plunge pool with pressure transducers. Then, they applied ANSYS CFX to simulate overflow nappe impingement jets in a general way, and investigated the different turbulence closures which better represent the data. In those works, the emphasis was put on the turbulence of the falling jet, the pressure distribution near the stagnation point in the water cushion, and the horizontal distance to the stagnation point. A good agreement among numerical simulations and laboratory data was obtained.

In Castillo et al. (2014), the so-called ‘homogeneous’ theoretical model of CFX was employed. It was shown that this model is able to reproduce correctly the jet water velocity, and the averaged pressures in the plunge pool. There is always a challenge in modelling two-phase flows to discern which level of complexity is needed to represent different aspects of the flow (Bombardelli 2004; Bombardelli & Jha 2009; Jha & Bombardelli 2010). One of the objectives of this paper is to determine whether this theoretical model is sophisticated enough to represent velocities in the plunge pool.

Continuing the line of research, this work presents a systematic study which considers specific flows and water cushions in the plunge pool. New laboratory data were obtained and new three-dimensional simulations were specifically performed for this work. ANSYS CFX was again selected due to the variety of turbulence closures available in the code, the previous experience with it and, more importantly, due to the diverse two-phase flow models embedded in the package, which can allow us to expand the research further in the near and long-term future.

From an engineering point of view, knowing the parameters analysed, designers will be able to estimate the scour effects and the stability of the dam with greater certainty.

The experimental facility was constructed at the Hydraulics Laboratory of the Universidad Politécnica de Cartagena, and was described in detail in Castillo & Carrillo (2012, 2013) and Carrillo & Castillo (2014). Here, we revisit its main features. The facility consists of a mobile mechanism which permits variation of the weir height between 1.7 and 3.5 m and flows from 10 to 150 L/s. It has an inlet channel with a length of 4.0 m and width of 0.95 m. The discharge is produced through a sharp-crested weir with a width of 0.85 m and height of 0.37 m.

The plunge pool, in which different water cushions may be simulated, is a 1.3-m high, 1.1-m wide and 3-m long methacrylate box. Turbulent kinetic energy values at the inlet channel were obtained with an acoustic Doppler velocimeter (ADV); mean velocities and air concentrations in different sections of the falling jet were acquired with optical fibre instrumentation; and instantaneous pressure values were measured with piezoresistive transducers located on the basin bottom. In addition, ADV and optical fibre probe were used in the basin to obtain velocity and air concentration profiles, respectively, downstream of the impingement point.

The flow volumetric rate was measured with a V-notch weir and was tested with ADV measurement upstream from the weir. Differences were smaller than 5%.

The setting characteristics of the ADV employed in these tests were determined by considering that the main objective was to characterise the turbulence. The velocity range was selected as ±4.00 m/s (the maximum available in the equipment). With this setting, the ADV was able to measure horizontal velocities up to 5.25 m/s and vertical velocities up to 1.50 m/s. Due to the fact that the Doppler equipment needs to be totally submerged into water, the first 5–6 cm of the water column could not be measured.

We used the turbulent kinetic energy measured 0.50 m upstream of the weir in the experimental facility to specify inlet boundary conditions in the numerical simulations in that location. This distance guarantees that the stream lines are horizontal upstream of the sharp weir crest (0.50 m > 5 h).

The plunge pool was surveyed in cross sections spaced 0.10 m horizontally. Four specific flows were tested (0.023, 0.037, 0.058 and 0.082 m²/s) with different water cushion depths. This covers 24 different configurations generated downstream of a rectangular free-falling jet.

In order to characterize the macro turbulence of the flow in the plunge pool, 5,000 values were recorded in each measured point by using a frequency of 10 Hz (more than 8 minutes of observation). In this way, 2006 points in the symmetrical vertical plane of the basin were obtained. As the flows are highly turbulent, the values obtained with ADV may be affected by spurious signals or ‘spikes’. Each time series must be filtered with the velocity and acceleration thresholds (Wahl 2000). Furthermore, in this particular case, the air may also affect the signal of the ADV. Frizell (2000) studied the air effect measuring concentrations varying from 0 to 3.61%. As the air concentrations increase and bubble sizes increase, correlation values drop dramatically as the acoustic signals used by the probe are absorbed and reflected by the two-phase flow mixture. Matos et al. (2002) also found that air bubbles affect the accuracy of velocity measurements taken with the ADV. However, their experimental results suggest that the ADV can provide reasonable estimates of the velocity for low air concentrations up to 8%.

To measure the air concentration at the falling jet and at the basin, an RBI-instrumentation dual-tip probe optical fibre phase-detection instrument was used. The rise and fall of the probe signal corresponds, respectively, to the arrival and the departure of the gas phase at the tip of the sensor. The thresholding values were set to 1.0 and 2.5 V (Boes & Hager 1998). The void fraction was defined as the ratio of the total time the probe is in the gas (Σt_Gi) to the experiment duration time t.

According to Stutz & Reboud (1997a, 1997b), the equipment enables measurement in water up to 20 m/s flow velocity and the relative uncertainty concerning the void fraction is estimated at about 15% of the measured value. One source of error in the estimation of air presence in the flow is due to the counting statistic of the number of air bubbles in contact with the tips of the probe (Stutz 1996). Therefore, a short duration of the sequence would contribute to an increased inaccuracy of the result. In order to evaluate the minimal measuring duration, André et al. (2005) considered the stabilization of the mean value during the measuring sequence and the quasi-stationary of the air concentration signal as statistically representative for the air concentration. Based on the sensitivity study of the probe behaviour, the authors recommend a sampling sequence of 60 s as a good compromise between accuracy and time consumption.

Boes & Hager (2003) carried out experiments with 4,000 air bubbles and sampling sequences of 30 seconds. The authors considered the accuracy of the air concentration and velocity measurements is related to the variation of the phase N_b, air to water or inversely, rather than to the sampling duration t.

The maximum air concentration is around 12% (at a distance of 21% from the bottom) for the first sections. However, from the section 0.30 m and a distance from the bottom smaller than 70%, the air concentration is below 10%. Concentrations remain high still at the upper portion of the water depth in the basin.

Considering highly turbulent and aerated flow that occurs in the basin of energy dissipation, the Phase-Space Thresholding filter (Goring & Nikora 2002, modified by Castillo 2009) was used. The spikes were replaced on each record by the mean value of the twelve closest points.

The main steps of the method are as follows:

1. Calculate surrogates for the first, ⁠, and second, ⁠, derivatives using a centered differences scheme.

2. Calculate the standard deviations of all three variables, ⁠, ⁠, and ⁠, and thence the expected maxima using the Universal criterion

3. Calculate the rotation angle of the principal axis of versus ⁠, using an expression for the cross correlation. Instead of using the expression by Goring & Nikora (2002):

   

   3a

Castillo (2009) proposed a new relation obtained by a Gauss’ fit:

• 

  3b

• 4. For each pair of variables, calculate the ellipse that has maxima and minima from point 3.

  • • For versus ⁠, the major axis is and, the minor axis is ⁠.

  • • For versus ⁠, the major axis is and, the minor axis is ⁠.

  • • For versus ⁠, the major and minor axes, a and b, respectively, are the solutions of (Goring & Nikora 2002):

    

    4a

    

    4b
    Castillo (2009) proposed the following system of equations instead:

    

    4c

    

    4d

From these plots it can be concluded that the filtering attenuates the standard deviation. Further filtering affects the values of the turbulent kinetic energy.

The theoretical model (1) comes as a result of the addition of the equations of the two phases (Drew & Passman 1999; ANSYS CFX Manual 2015). Further, ⁠. Rigorously speaking, models like this have been found to provide adequate predictions only in relatively-dilute mixtures. Jha & Bombardelli (2010) established that dilute mixtures can extend up to the range 2–5% in the context of sediment-laden flows; something similar could be assumed in bubble flows. For larger concentrations, Jha & Bombardelli (2010) found that the velocity distribution could not be well predicted relatively far from the wall with mixture models. Thus, we expect the ‘homogeneous’ model to be able to represent rather adequately those areas in which air concentrations are not that high (cf. Figure 5).

As said, the code ANSYS CFX has been used, which is based on an element-oriented, finite-volume method. It allows different types of volumes, including tetrahedral and hexahedral volumes. Solution variables are stored at the nodes (mesh vertices). More details are given in the ANSYS CFX Manual (2015).

Finally, the instantaneous values in the numerical simulations show large variations during the establishment of the quasi-steady conditions. Once the quasi-steady conditions are reached, there are small fluctuations around the mean value.

The Standard k-ɛ model is considered to be the traditional turbulence model and it is embedded in the majority of the CFD programs. The Re-Normalisation Group (RNG) k-ɛ model is based on the renormalisation group analysis of the Navier–Stokes equations (Yakhot & Orszag 1986; Yakhot & Smith 1992). The transport equations for turbulent kinetic energy and dissipation rate are similar as those for the standard model, although their respective constants are different.

The k-ω based Shear-Stress Transport (SST) model (Menter 1994) assumes that the eddy viscosity is linked to the turbulence kinetic energy, k, and the turbulent frequency, ω, as ⁠. The SST model takes into account the accuracy of the k-ω model in the near-wall region and the free stream independence of the k-ɛ model in the outer part of the boundary layer. It is considered as a hybrid model (see Rodi et al. 2012).

The free surface model assumes that each control volume contains three possible conditions:

• r_k= 0 if cell is empty (of the k-th fluid);

• r_k= 1 if cell is full (of the k-th fluid);

• 0 < r_k < 1 if cell contains the interface between the fluids.

Tracking of the interface between fluids is accomplished by the solution of the volume fraction equation (see ANSYS CFX Manual 2015).

Like in the simulation of the pressure field on the basin bottom (Castillo et al. 2014), all scenarios were obtained by a calculation time of 60 seconds, using a constant time-step of 0.05 s. The transient statistics were obtained by considering that quasi-steady-state conditions are reached after the first 20 seconds of simulation (Castillo et al. 2014). Hence, statistics were obtained from time-step 400 to time-step 1,200 (800 data in 40 s).

ANSYS CFX has different treatments near the wall. ω-based turbulence models (e.g. SST) use automatic wall treatments which switch between regular wall functions (Pope 2000) and low-Reynolds wall treatment (Menter 1994). Wall functions are used when the wall adjacent vertices are in the log-law layer (y⁺ ≈ 20–200). The low-Reynolds wall treatment is used when the wall adjacent vertices are in the viscous sublayer (ANSYS CFX Manual 2015). Considering the wall treatment used by ANSYS CFX, the mesh sizes close to the solid boundary were smaller than in the rest of the domain. Values of y⁺ were smaller than 40.

For simplicity, only a half part of the model was simulated. The symmetry condition in the longitudinal plane of the plunge pool was used.

The inlet condition considers the volumetric flow rate with a normal direction to the boundary (q = 0.082 m²/s, q = 0.058 m²/s, q = 0.037 m²/s, q = 0.023 m²/s), the turbulent kinetic energy (5.1 × 10^–4 m²/s² for q = 0.082 m²/s, 3.6 × 10^–4 m²/s² for q = 0.058 m²/s, 1.9 × 10^–4 m²/s² for q = 0.037 m²/s, and 1.1 × 10^–4 m²/s² for q = 0.023 m²/s), and the water level height at the upstream deposit (2.313 m for q = 0.082 m²/s, 2.285 m for q = 0.058 m²/s, 2.263 m for q = 0.037 m²/s, 2.237 m for q = 0.023 m²/s).

The outlet condition was considered with flow normal to the boundary and hydrostatic pressure. The water level height at the outlet was modified according to the water cushion depth, Y, in the laboratory device. For all walls of the upper deposit, the weir and the dissipation basin, no slip smooth wall conditions were considered. The roughness of methacrylate was indicated in the walls. In the transverse direction, wall boundary conditions were used.

To judge the convergence of iterations in the numerical solution, we monitored the residuals (Wasewar & Vijay Sarathi 2008). The solution is said to have converged in the iterations if the scaled residuals are smaller than fixed values ranging between 10^–3 and 10^–6. In this work, the residual values were set to 10^–4 for all the variables. Each time-step required around 12 iterations to reach the convergence criterion. With this choice and for 791,354 elements (255,776 in the falling jet), the mean computational time was 7.2 × 10⁵ seconds (≈8.3 days), using a Central Processing Unit (CPU) with sixteen processors (Intel^® Xeon^® E5-2699 v3 @ 2.30 GHz).

In this section, we present some formulas obtained from the literature to interpret the different hydraulic parameters of plunge pools. To that end, we selected expressions for submerged hydraulic jumps and horizontal wall jets.

To characterize the non-dimensional velocity distribution in hydraulic jumps, some authors have proposed different expressions which we include in the Appendix. Using the results obtained by diverse researchers, Wu & Rajaratnam (1996) considered the length scale δ_l for the submerged jumps as a function of the distance to the beginning of the jump. They found that most of the observations for submerged jumps are contained within one standard deviation of the mean value of the wall jet, and only the data points near the end of the jump show an accelerated growth rate.

In general, the measured horizontal velocities, V_x, are greater than the calculated counterparts. The bellies which appear in some profiles seem to be related with the three-dimensional flow features of the flow. This level of agreement is understandable given the observed air concentrations, close to 10%. As said with this level of concentrations of the disperse phase, Jha & Bombardelli (2010) found that homogeneous-type models can only provide approximate values of velocity. In Figure 15(a)–15(c) we can see that the choice of the turbulence model has a significant effect over the result of the horizontal flow profiles. The less accurate results were obtained with the RNG k-ɛ turbulence model. However, the SST and the standard k-ɛ turbulence models produce velocities fairly close to the values measured in the laboratory.

In addition to the mean velocities, the turbulent kinetic energy profiles were also compared (Figure 15(d)). Like in the mean velocity cases, the best results of turbulent kinetic energy profiles were obtained with the SST model. In general, the results from the numerical simulations show the same behaviour as the results obtained in the laboratory. However, the differences are very important close to the stagnation point, where the numerical model may not obtain accurate results due to the relatively important air entrainment into the plunge pool (Figures 3 and 5).

Observing and predicting two-phase flows in hydraulic structures is very complicated due to the rather non-dilute nature of the flow. Under non-dilute conditions, both experiments and simulations cannot be expected to lead to clean comparisons.

In this work, mean velocity and turbulent kinetic energy profiles have been analyzed in a plunge pool located downstream of a rectangular free-falling jet. The energy dissipation in a plunge pool is very high. For the test carried out, the dissipation of the impingement jet energy was between 75 and 95%. This ratio increases when the ratio water cushion/impingement jet thickness decreases.

In general, the CFD simulations provided results fairly close to the values measured in the laboratory, and to the formulas proposed by diverse authors, in spite of having used a simple two-phase flow model. ‘Homogeneous’ model seems to be able to predict rather well areas in which air concentration is not very high. However, in the highly aerated regions rather strong differences appear.

Regarding the modelling of turbulence, three closures were tested. SST turbulence model offered results closer to the laboratory measurements. The RNG k-ɛ model tends to underestimate the turbulent kinetic energy.

It was possible to propose a single mean velocity distribution law for ratios V_x/V_max ≥0.40. For smaller values, there are necessarily diverse distribution laws.

In order to develop this work further, we plan to examine air entrainment in the stilling basin. Comparison of results with diverse CFD codes (open source and commercial ones) against data will be considered.

The first two researchers express their gratitude for the financial aid received from the Ministerio de Economía y Competitividad and the Fondo Europeo de Desarrollo Regional (FEDER), through the Natural Aeration of Dam Overtopping Free Jet Flows and its Diffusion on Dissipation Energy Basins project (BIA2011-28756-C03-02). The first author acknowledges the support of Ministerio de Educación, Cultura y Deporte of España through Estancias de Movilidad de Profesores e Investigadores Senior program (PRX 14/00367), which allowed him to develop a stay as a Visiting Scholar Researcher at the Department of Civil and Environmental Engineering of the University of California, Davis, from April to October 2015.

Expressions of the velocity distribution in wall jets and hydraulic jumps which have been provided by different authors are:

Author .	Formula .
Görtler (1942), cited by Liu et al. (1998)
Rajaratnam (1976)
Ohtsu et al. (1990)
Lin et al. (2012)
where

Author .	Formula .
Görtler (1942), cited by Liu et al. (1998)
Rajaratnam (1976)
Ohtsu et al. (1990)
Lin et al. (2012)
where

View Large