**Numerical Study of the Collapse of Multiple Bubbles and the Energy Conversion during Bubble Collapse**

**Jing Zhang 1,2,3, Lingxin Zhang 1,2,3,\* and Jian Deng 1,2,3**


Received: 10 December 2018; Accepted: 29 January 2019; Published: 31 January 2019

**Abstract:** This paper investigates numerically the collapses of both a single cavitation bubble and a cluster consisting of 8 bubbles, concerning mainly on the conversions between different forms of energy. Direct numerical simulation (DNS) with volume of fluid (VOF) method is applied, considering the detailed resolution of the vapor-liquid interfaces. First, for a single bubble near a solid wall, we find that the peak value of the wave energy, or equivalently the energy conversion rate decreases when the distance between the bubble and the wall is reduced. However, for the collapses of multiple bubbles, this relationship between the bubble-wall distance and the conversion rate reverses, implying a distinct physical mechanism. The evolutions of individual bubbles during the collapses of multiple bubbles are examined. We observe that when the bubbles are placed far away from the solid wall, the jetting flows induced by all bubbles point towards the cluster centre, while the focal point shifts towards the solid wall when the cluster is very close to the wall. We note that it is very challenging to consider thermal and acoustic damping mechanisms in the current numerical methods, which might be significant contributions to the energy budget, and we leave it open to the future studies.

**Keywords:** cavitation bubble, multiple bubble collapse, pressure wave energy, energy conversion rate

### **1. Introduction**

As a natural phenomenon, cavitation occurs frequently in a variety of engineering applications. As the local pressure in liquid is lower than its saturated vapor pressure, gas nucleus can develop into cavitation bubbles. When the local pressure is recovered or driven by high pressure waves, the cavitation bubbles collapse, leading to a series of physical scenarios such as micro jetting flows, shock wave emission, heat release etc. In hydraulic systems, cavitation can cause material damage, loss of power, and induce noises. On the other hand, bubble collapse can also be utilized as a tool in the fields of medical treatments and mining, such as drug delivery, ultrasound operation and cavity jet exploits [1]. Numerous experimental and numerical studies have been performed to study the bubble dynamics and the collapses of cavitation bubbles, in order to control or reinforce the effects of cavitation bubble collapse. However, the dynamics of multiple cavitation bubbles and their collapses, particularly the underlying behavior of energy conversion have not been well understood.

Historically, studies on bubble dynamics can date back to 1917, when Reyleigh proposed a theoretical spherical-bubble model [2]. This model was for an ideal spherical bubble, neglecting the effects of viscosity, surface tension, compressibility and mass transfer. It was followed by a series of modifications and perfections [3–8].

In the past few decades, a large amount of works have been undertaken to understand the characteristics of bubble collapse [9–24]. Lauterborn recorded the process of a single bubble collapsing near a rigid wall, in which the shock wave emission was captured by a schlieren system [25–27]. Moreover, the particle image velocimetry (PIV) technique was adopted to monitor the velocity fields by lauterborn & Vogel [28], who captured the micro jetting flows as the bubble collapses asymmetrically. They revealed the relationship between the noises induced by the bubble collapse and the dimensionless distance *γ* (*γ* = *L*/*Rmax*, where *L* is the distance between the bubble centre and the rigid wall, and *Rmax* is the maximum bubble radius). Gonzalez-Avila measured the pressure emission of laser-induced bubble collapse [29], in which the pressure amplitude was found to be up to 1*k* Bar when the bubble was close to the rigid boundary [30]. Merouan [31] measured the temperature variations during the growth and collapse of a single cavitation bubble. More recently, Fortes et al. proposed an analytic model to examine the pressure wave emission during bubble collapse, with the conversion rate of the wave energy found to be strongly influenced by gas content in the bubble [32].

However, not many efforts have been made to study the pressure wave emission from the collapses of multiple bubbles due to the difficulties in experimental setup and the consideration of compressibility of liquid in numerics [33–36]. Generally, a bubble cluster collapses inward, generating a high level impulsive pressure wave, and the symmetries for individual bubbles are lost due to the presence of the rigid wall and their nearby bubbles [37,38]. It has been shown that the shock wave emission focuses inward during the bubble cloud collapse, leading to a very high peak pressure, which can exceed that of the collapse of a single bubble with the equivalent volume [39]. Tiwari et al. [40] simulated the expansion and subsequent collapse of a hemispherical cluster of 50 bubbles adjacent to a plane rigid wall, with the detailed compressible-fluid mechanics of bubble-bubble interactions considered. They found that the peak pressure was associated with the centremost bubble, which caused a corresponding peak pressure on the nearby wall. Despite these advances in the few existing researches, the detailed bubble-scale dynamics of a cluster collapse remain poorly understood.

Here, we investigate the characteristics of the collapses of multiple bubbles by solving directly the Navier-Stokes equations considering the compressiblility of the fluid, with particular concerns on the energy conversion. The volume of fluid method is involved to capture the interfaces between the liquid and the cavitation bubbles. The simulation results are organized into two parts: in the first part, a single bubble is studied as a reference to our further multiple-bubble studies. In the second part, we focus the collapses of multiple bubbles, presenting the unique dynamic behaviors in pressure wave propagation and energy conversion, in contrast to the single bubble situation.

### **2. Numerical Methods**

The numerical work is performed using a finite-volume based code. Assuming that the liquid is compressible and the cavitation bubbles are filled with pure vapour, and considering viscosity and surface tension, the governing equations are formularized as the following

$$\frac{d\rho}{dt} + \rho \nabla \overline{\mathcal{U}} = 0,\tag{1}$$

$$
\rho \frac{\partial \overrightarrow{\boldsymbol{U}}}{\partial t} + \rho \overrightarrow{\boldsymbol{U}} \cdot \nabla \overrightarrow{\boldsymbol{U}} = \rho \overrightarrow{\boldsymbol{g}}' - \nabla p + 2\nabla \cdot \left(\mu \overrightarrow{\boldsymbol{D}}\right) - \frac{2}{3} \nabla \left(\mu \nabla \cdot \overrightarrow{\boldsymbol{U}}\right) + \sigma \kappa \overrightarrow{\boldsymbol{N}}' \tag{2}
$$

$$\frac{d\mathbf{a}}{dt} = \frac{\partial \mathbf{a}}{\partial t} + \overrightarrow{\mathbf{U}} \cdot \nabla \mathbf{a} = \mathbf{0},\tag{3}$$

where *<sup>k</sup>* is surface curvature, *<sup>σ</sup>* is the surface tension coefficient, *<sup>D</sup>* is the strain rate tensor, −→*<sup>N</sup>* is the unit normal vector of the interface, *α* is the volume fraction of the liquid, and *ρ* and *μ* are the density and the viscosity respectively of the mixture fluid, which are obtained by weighting of the volume fractions:

$$
\rho = \alpha \rho\_1 + (1 - \alpha) \rho\_2. \tag{4}
$$

$$
\mu = \kappa \mu\_1 + (1 - \alpha) \mu\_2. \tag{5}
$$

where subscripts 1 and 2 denote liquid phase (water) and gas phase (vapour), respectively. The densities of each phase are calculated by

$$
\rho\_1 = \rho\_{10} + \psi\_1 p\_\prime \tag{6}
$$

$$
\rho\_2 = \rho\_{20} + \psi\_2 \mathbf{p}\_{\prime} \tag{7}
$$

where *ψ*<sup>1</sup> and *ψ*<sup>2</sup> are two constants associated with the respective compressibilities of water and vapour.

In our simulations, the space discretizations are second-order upwind for the convection terms and central differences for the Laplacian terms, respectively. The time discretization is first-order implicit Euler. The pressure-velocity coupling is obtained using the Pressure Implicit Split Operator (PISO) scheme. The preconditioned conjugate gradient (PCG) method is used to treat the pressure equation and the preconditioned biconjugate gradient (PBiCG) method is used for the velocity equations.

The computational domain is a cylinder with a radius of 25 mm, and a height of 50 mm. The grid distributions on some specific sections are shown in Figure 1. The bottom of the cylindrical flow domain is set as a rigid wall, the others are far field boundaries. The initial pressure of the liquid phase is 101,325 Pa, and the pressure of vapor is 3154 Pa. The initial flow is quiescent, and the initial bubbles are spherical, with the radius of 2 mm. The geometrical setup of the multiple-bubble case is shown in Figure 2. The dimensionless wall distance is defined as *γ* = L/*Rmax*. After carefully carrying out self-consistency tests, we find that the current grid with a cell number of 25,581,192, and an initial maximum time step size of *dt* = <sup>1</sup> × <sup>10</sup>−<sup>7</sup> s, are sufficient to assure satisfactory independence of the results with respect to both mesh and time discretizations. We note that the time step is adjustable during the simulations to meet the requirement of local courant number *Co* = 0.35.

(**a**) Top view. (**b**) Side view.

**Figure 1.** (**a**) Top view and (**b**) side view of the grid distributions within the computational domain. Note that the grid in the central area of the cylindrical domain is refined where the bubbles are placed, and only a coarse grid with 10% of the cells is displayed to diagram the grid distributions.

**Figure 2.** Geometrical description for the simulations of multiple bubbles.

To further validate the spatial resolutions, we carry out simulations with different numbers of cells along the diameter of a single bubble. Figure 3 presents the time histories of the dimensionless radius *R*∗, where *R*<sup>∗</sup> = *R*/*Rmax*, and the time is non-dimensionalized according to *T*<sup>∗</sup> = *t*/*tc*. Here, *Rmax* is the maximum or the initial radius of the bubble, and t*c* is the collapsing time predicated from the Reyleigh-plesset equation. It is observed that the result converges at 34 cells, which accords well with the theoretical solution. We show in Figure 4a the evolutions of bubble shapes during collapse simulated based on the current resolution of 34 cells. The numerical results are consistent with the previous experimental observation [41].

**Figure 3.** Time histories of radius for a single bubble resolved by different number of cells distributed along the diameter. The reference curve is the theoretical solution of the Rayleigh-Plesset equation.

**Figure 4.** Variations of the bubble shapes during collapse obtained by (**a**) the present numerical simulations and (**b**) previous experiments [41].

In the following sections, according to [42], the potential energy of a bubble is defined as

$$E\_{\rm pot} = \frac{4}{3} \pi R^3 \Delta p\_\prime \tag{8}$$

where *R* is the bubble radius and Δ*p* is the pressure difference between two sides of the bubble interface, then the maximum potential energy of a cavitation bubble is defined as *Epot*−*max*. The total kinetic energy of the flow domain is the integral:

$$E\_k = \int \frac{1}{2} \rho l I^2 dV.\tag{9}$$

We also define the wave energy emitted by the bubble collapse, following an acoustic approach [28,43], with the following form:

$$E\_{\text{wave}} = \int \frac{\Delta P^2}{\left(\rho c\right)^2} dV\_\prime \tag{10}$$

where *c* is the sound speed in water, *c* = 1500 m/s in the current simulations. The conversion rate of the wave generation is defined as

$$
\eta = \frac{E\_{\text{wave}}}{E\_{\text{pot}-\text{max}}} \tag{11}
$$

### **3. Result**

### *3.1. Collapse of a Single Bubble*

In this section, we study a single cavitation bubble collapsing near a solid wall. To understand the presence of the solid wall on the non-spherical deformation, and consequently the energy conversion, we carry out a series of simulations by varying the distance between the bubble and the solid wall. The dimensionless distance *γ* ranges from 1.5 to 4.0.

In Table 1, we present the simulation results. First, we observe that the energy conversion rate reaches its maximum of 29.29% at *γ* = 12.5, which is not difficult to understand since the *γ* = 12.5 case corresponds to a distance very far from the solid wall, which we believe that the wall effect can be reasonably neglected. The values of maximum kinetic energy for all cases are around 4.0 mJ, with slight variations among different distances. We note that the energies and the conversions at *γ* = 2.5 and *γ* = 4.0 are very close, exhibiting nontrivial bounding effects of the solid wall.

**Table 1.** Energies and their conversion rates for various distances for the collapse of a single bubble (note the different units for the two energies).


In Figure 5a,b, we present the time histories of different forms of energy at *γ* = 12.5 and *γ* = 1.5 respectively, to make a direct comparison between a small and a large distances. In Figure 5a, we observe that the potential energy of the bubble decreases during the process of collapse as the bubble shrinks, while the kinetic energy grows from zero, representing the accelerated process of the bubble collapse. As the bubble collapses to a singular point at 199 μs, the kinetic energy reaches its maximum value, while the potential energy drops to nearly zero value. During this whole period, the total energy is approximately equal to the initial bubble potential energy. In other words, the total mechanical energy remains approximately constant, or is conserved. During this process of bubble shrinkage, the potential energy is lost to the re-entrant jetting towards the bubble centre.

After the bubble vanishes, around 200 μs, a pressure wave is emitted and travels outwards, resulting in an intensive wave energy observed in Figure 5a. However, since we consider only limited forms of energy in the simulations and neglect thermal and acoustic damping mechanisms, which might be significant, the total energy is not conserved anymore. We understand that it is very challenging to quantify in detail any switch-over of these dominant damping mechanisms

in such violent bubble collapses. We leave it open to the future studies. To help understand this process more intuitively, we present in Figure 6 the flow fields for three time instants demonstrating the alternative dominance of different energy forms. In Figure 6a, at the early stage, the bubble potential energy dominates, while it is converted to kinetic energy when the bubble starts to collapse (see Figure 6b), and in Figure 6c the pressure wave is emitted with only the wave energy prominent in the energy budget.

**Figure 5.** Time histories of different forms of energy in the flow domain during the collapse of a single bubble at (**a**) *γ* = 12.5 and (**b**) *γ* = 1.5. For both cases, *Rmax* = 2 mm and ΔP = 101,325 Pa.

For comparison, we present the time histories of different forms of energy at *γ* = 1.5 in Figure 5b, in which case the initial bubble is very close to the solid wall. Apparently, the changing trend of each form of energy is roughly the same with its large distance counterpart (see Figure 5a). They still differ in some aspects: first, the wave energy appears later than the large distance case due to the confinement of the solid wall, which delays the collapse of the bubble. Second, during the whole process, even after the kinetic energy drops, it is always higher than the wave energy, implying a different physical mechanism. We conjecture that the presence of a rigid boundary can lead to asymmetric collapse, which requires more time. In the process of asymmetric collapse, high speed micro jetting flows are induced towards the rigid boundary, the liquid can thus be pushed away along the wall instead of collapsing towards a singular point. Therefore, it is easy to understand that the kinetic energy dominates even after the bubble collapse.

To get a comprehensive understanding of the distance effects, in Figure 7 we present the variations of wave energy with time for various distances corresponding to Table 1. When the bubble is located far away from the rigid boundary, e.g., *γ* = 12.5, the wave energy reaches the maximum value of 0.95 mJ at *t* = 1.85 ms within a very short period of 0.1 ms, and with a very sharp peak. As the distance is reduced, the peak value of wave energy decreases accordingly, and the sharp peak is replaced by a broad distribution of the high wave energy. In specific, the peak value of wave energy for *γ* = 1.5 is around 0.3 mJ, only one-third of that for *γ* = 12.5. It is interesting to find that there are secondary peaks in the wave energy distributions at *γ* = 2.5–4.0. We believe that these secondary peaks are induced by the pressure waves reflected by the solid boundary, which are difficult to be distinguished as we further reduce the distance. The comparison between small and large distances suggests that the near-wall collapses generate less wave energy due to the non-spherical characteristics, but with a longer duration in the computational domain.

**Figure 6.** Flow fields with pressure contours for three time instants at *γ* = 12.5, demonstrating different dominant forms of energy: (**a**) potential energy, (**b**) kinetic energy and (**c**) wave energy.

Figure 8 shows the variations of kinetic energy with time for various distances. Again, their overall trends are the same, as we have discussed in Figure 5, and they are consistent with the wave energy variations shown in Figure 7. There are also some slight differences among various distances. In short, the wall effects on kinetic energy can be concluded as: first, the drop of kinetic energy is delayed as the bubble is placed closer to the wall. Second, the sharp peak is replaced by broader distributions of the kinetic energy as the bubble stays closer to the wall. Both conclusions are consistent with the wave energy variations, as shown in Figure 7.

**Figure 7.** Wave energy variations in the flow fields at different values of *γ*.

**Figure 8.** (**a**) Kinetic energy variations in the flow fields at different values of *γ*, and the enlarged views around the peak regions are shown in (**b**).

### *3.2. Collapses of Multiple Bubbles*

The asymmetric bubble collapse can be brought by placing a solid wall, as we have studied in the last section for a single bubble, it can also be caused by the presence of nearby bubbles. In this section, we study the collapses of multiple bubbles, or more specifically, 8 bubbles, which are placed in a cubic region with two layers, and 4 bubbles on each layer.

In Table 2, we present the peak values of both the kinetic and wave energies, as well as their conversion rates. The value of E*<sup>k</sup>* slightly grows with the distance *γ*, while *η* decreases as *γ* increases. In specific, the energy conversion rate at *γ* = 1.5 is approximately twice of that at *γ* = 11.25.

To examine the dynamic evolutions of individual bubbles, in Figures 9 and 10, we present the instantaneous shapes of the bubbles in three subsequent time instants. As shown in Figures 9a and 10a for the initial bubble configurations, the two cases represent respectively the small (*γ* = 1.5) and large (*γ* = 11.5) distances of the bubbles away from the solid wall. For the small distance, as shown in Figure 9, the wall effect leads to asymmetric collapses of the bubble cluster, with the bubbles in the upper layer collapsing faster than that in the lower layer, due to the solid wall below the lower layer. While at *γ* = 12.5, as shown in Figure 10, the asymmetry is brought only by the presence of nearby

bubbles, therefore, the collapses and the evolutions of bubbles are symmetric with respect to the symmetry plane between the two layers, which is a straightforward demonstration that the solid wall does not affect the bubble collapse markedly when they are spaced sufficiently far away. The two cases differ in the focal point. For the large distance case, the jetting flows from all bubbles point to the centre of the bubble cluster, while for the small distance case, the focal point shifts towards the solid wall.

**Table 2.** Energies and their conversion rates for multiple bubbles.

**Figure 9.** Deformations of multiple bubbles at *γ* = 1.5 for (**a**) the initial shapes (side view), (**b**) T = 80 μs, (**c**) T = 340 μs and (**d**) T = 400 μs.

Figure 11a presents the potential energy variations of multiple bubbles. The drops of potential energy slow down at the small distances, which is apparent, because the collapses of bubbles have been delayed due to the confinement of the solid wall. Comparing to the single bubble, for the multiple bubbles, the solid wall plays a similar role, as the bubbles are very close to the wall. However, the trend of wave energy conversion is quite different with that of the single bubble, which might be raised by the strongly nonlinear interactions between the bubbles. Similar comparisons can be made in the kinetic energy variations, as shown in Figures 8b and 11b.

**Figure 10.** Deformations of multiple bubbles at *γ* = 11.5 for (**a**) the initial shapes (side view), (**b**) T = 80 μs, (**c**) T = 315 μs and (**d**) T = 340 μs.

**Figure 11.** (**a**) Potential energy and (**b**) kinetic energy variations in the flow fields for the collapses of multiple bubbles.

The most remarkable difference between the single bubble and the multiple bubbles lies in the wave energy. Figure 12 presents the wave energy variations in the simulations of multiple bubbles. In contrast to that of the single bubble (see Figure 7), here, the ascending and descending branches of the wave energy variations are nearly symmetric with respect to the peak location. The peak values of wave energy for the small distance cases are higher than that of the large distance cases, resulting in the same behaviour for the energy conversion rate. This relationship between the peak wave energy is opposite to that of the single bubble.

**Figure 12.** Variations of wave energy in the flow fields during the collapsed of multiple bubbles.

### **4. Conclusions**

In the present work, we carry out numerical simulations to study the collapses for both a single bubble and a bubble cluster with 8 bubbles. Different forms of energy are evaluated to understand the underlying physical mechanisms of bubble collapse.

The first part concentrates on a single bubble collapsing near a solid wall. The intuitive physical rule tells that the energy should be conserved, which has indeed been observed during the first stage of collapse. As the bubble shrinks, the potential energy determined by the vapour volume converts into kinetic energy. Further, as the bubble vanishes, a part of kinetic energy starts to convert into other energy forms such as wave energy, or is lost due to the thermal and acoustic damping mechanisms, which have unfortunately not been considered in the current simulations. We report that a rigid boundary can affect the process of bubble collapse by deforming the bubble into non-spherical shapes, delaying the collapse, and consequently decreasing the conversion rate.

For a cluster consisting of 8 bubbles, their collapses are more complicated than the single bubble case, because the individual bubbles can be affected by both the solid wall and the surrounding bubbles, and both can break the symmetry of their geometric configurations. By examining the evolutions of individual bubbles during the collapse, we observe that when the bubbles are far away from the solid wall, the jetting flows from all bubbles point towards the cluster centre, while the focal point shifts towards the solid wall as the cluster is located very close to the wall. Moreover, we find that the collapse differentiates from the single bubble case mainly in the relationship between the bubble-wall distance and the wave energy, or equivalently the energy conversion rate.

**Author Contributions:** Conceptualization, J.Z. and L.Z.; Methodology, L.Z. and J.Z.; Investigation, J.Z.; Data Curation, J.Z.; Writing-Original Draft Preparation, J.Z. and L.Z.; Writing-Review and Editing, J.D.; Visualization, J.Z.; Funding Acquisition, L.Z. and J.D.

**Funding:** This research was funded by the National Natural Science Foundation of China (No. 11772298, No. 11272284) and the State Key Program of National Natural Science of China (No. 11332009).

**Conflicts of Interest:** The authors declare no conflict of interest.

### **References**


c 2018 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

### *Article* **Numerical Analysis on the Hydrodynamic Performance of an Artificially Ventilated Surface-Piercing Propeller**

### **Dongmei Yang \*, Zhen Ren, Zhiqun Guo and Zeyang Gao**

College of Shipbuilding Engineering, Harbin Engineering University, Harbin 150001, China; forfanqiang@126.com (Z.R.); guozhiqun@hrbeu.edu.cn (Z.G.); gaozeyangpc@hotmail.com (Z.G.) **\*** Correspondence: yangdongmei@hrbeu.edu.cn; Tel.: +86-594-8258-8360

Received: 11 September 2018; Accepted: 18 October 2018; Published: 23 October 2018

**Abstract:** When operated under large water immersion, surface piercing propellers are prone to be in heavy load conditions. To improve the hydrodynamic performance of the surface piercing propellers, engineers usually artificially ventilate the blades by equipping a vent pipe in front of the propeller disc. In this paper, the influence of artificial ventilation on the hydrodynamic performance of surface piercing propellers under full immersion conditions was investigated using the Computational Fluid Dynamics (CFD) method. The numerical results suggest that the effect of artificial ventilation on the pressure distribution on the blades decreases along the radial direction. And at low advancing speed, the thrust, torque as well as the efficiency of the propeller are smaller than those without ventilation. However, with the increase of the advancing speed, the efficiency of the propeller rapidly increases and can be greater than the without-ventilation case. The numerical results demonstrates the effectiveness of the artificial ventilation approach for improving the hydrodynamic performance of the surface piercing propellers for high speed planning crafts.

**Keywords:** surface-piercing propeller; artificial ventilation; hydrodynamic performance; numerical simulation

### **1. Introduction**

Surface-piercing propellers (SPPs) are also known as Surface penetrating propellers. The SPP is so named because part of the propeller is above the water surface and the rest is under the water during normal operation. Compared with the conventional propellers, the SPPs mainly have the following three advantages (Ding et al. [1]): (1) the resistance of the appendages, such as the paddle shaft and the shaft bracket, is minimized; (2) propeller diameter is no longer limited by some parameters such as soak depth and stern frame; and (3) the cavitation erosion on the blade surface is substantially reduced. Due to these advantages, the SPPs become an optimal choice for high-speed boats and some shallow draft ships, and thus have really good application prospects.

However, in some cases the SPPs might be operated under the large water immersion, e.g., the boats are at the bow-up status, which makes the SPPs be in heavy load conditions. One of the practical approaches to solve this problem is setting an aeration pipe in front of SPPs to ventilate the blades. This approach has been adopted in some actual boats and achieved significant results. Nonetheless, few attentions have been paid on how the ventilation of blades improves the hydrodynamic characteristics of SPPs. To this end, this paper is dedicated to investigate the influence of the ventilation of blades on the hydrodynamic performance of SPPs using numerical methods.

In some simplified studies, the surface-piercing process of an SPP blade can be deemed as the water-entry process of a 2D profile, which is easy to be understood. Zhao et al. [2] used the non-linear boundary element method to investigate the characteristics of the sprays generated by the water-entry of wedges with different bottom-angles. It was found that the increase of the wetted area directly leads to a great slamming pressure on the wetted surface. And when the slope angle of the bottom is larger than 30◦, there is no typical slamming pressure concentration phenomenon due to the smaller wetted surface. Young et al. [3], in his doctoral dissertation, carried out similar studies on the water-entry of wedges, and analyzed the pressure distribution on the wetted surface of a flat plate under different entry angles.

Yari et al. [4] conducted a preliminary study on a wedge entering the water under different inclination angles. The obtained pressure distribution, free surface deformation, and so on agree well with the experimental data. Yu et al. [5] carried out a numerical study on the water entry problem of a 2D cross-section profile, and found that the cup shape of the trailing edge can enhance the hydrodynamic load on the edge, while reducing the lift-to-drag ratio as well as the efficiency of the propeller.

Ghassemi et al. [6] investigated the hydrodynamics of SPPs named SPP-1 and SPP-2 under full and half immersed conditions using a boundary element method (BEM). The numerical results agree well with experimental ones. In the study, they found that the Weber number has significant influence on the ventilation status of the SPP. Kinnas et al. [7] studied the bubble flow around the hydrofoil and SPPs using a BEM, in which they investigated the hydrodynamic effects of the bubble flow on the hydrofoil and SPPs with various blade profiles. Young et al. [8,9] carried out a series of studies on large-scale SPPs and super-vacuum paddles also by using the BEM, which contributed to the knowledge of the hydrodynamics of these propellers.

With the development of the computer technology, the CFD methods become popular and play the role of benchmark in the numerical studies on surface-piercing propellers. Young et al. [3] predicted the hydrodynamic performance of large-scale surface paddles using a coupled boundary element method–finite element method (BEM–FEM), and found that the numerical results compared well with those obtained by the Reynolds Averaged Navier–Stokes (RANS) solver (Fluent software), in which the Reynolds stress tensor was modeled using the SST form of the *k* − *ω* turbulence model. Shi et al. [10] investigated the wake of surface-piercing propellers, as well as the pulsating pressure on the paddles and the deformation of the free surface after the paddles. Alimirzazadeh et al. [11] employed the OpenFOAM software to study the hydrodynamic performance of SPP-841B surface paddle under different sway angles and immersion depths. In the study, the *k* − *ω* based Shear Stress Transport (SST) model with automatic wall functions (mixed formulation) was employed. The numerical results agree with experimental ones well under different immersion depths, while poorly under different sway angles. They found that increasing the yaw angle would decrease the thrust coefficient and torque coefficient, while the efficiency was improved. In addition, the maximum efficiency point has not been achieved at the zero shaft yaw angle, due to the implementation of SPPs.

Obviously, in the existing works the most attention has been paid to investigating the hydrodynamic performance of SPPs under the natural operation conditions, while very little attention has been paid on the artificial ventilated SPPs. It is believed that the hydrodynamic characteristics of the artificial ventilated SPPs should be significantly different from the without-ventilation case. Thereby, we were motivated to study the SPPs under artificial ventilation conditions. In this paper, a right-handed three-blade paddle was proposed for the investigation, and an aeration pipe was installed in front of the paddle disk for the ventilation purpose. The numerical simulation was carried out based on the commercial CFD software Star-CCM+ [12]. The rotation of the blades was realized by using the overlapped mesh technique.

### **2. SPP Model and Numerical Setup**

The submergence ratio of SPPs is defined as *It* = *h*/*D*, where *h* denotes the immersion depth of the propeller disc, *D* is the diameter of the propeller, as shown in Figure 1.

**Figure 1.** The submergence ratio of surface-piercing propellers (SPPs), which is defined as the ratio of the immersion depth *h* to the propeller diameter *D*.

Generally, the SPPs work on the water-air interface, and the blades rush into and out of the water in a staggered manner. In this process, the air is sucked into the water forming an air cavity, which connects the atmosphere and the paddle. This phenomenon is known as the ventilation of SPPs.

When the submergence ratio is sufficient large, the torque on the paddle will dramatically increase. For the adjustable SPPs, the torque can be reduced by decreasing the immersion depth. For the fixed SPPs, however, it is commonly difficult to adjust the immersion depth. Although the paddle may pump some air into water when it rotates, the torque can still be very large. For the purpose of torque reduction, an aeration pipe is generally installed in front of the paddle to ventilate the blades and thus reduce the torque.

### *2.1. SPP Model*

The SPP proposed for the calculation is the SPP-1 type surface-piercing propeller model given in Table 1 [6], with the rotation direction being right-handed and the profile being S-C type. The diameter of the aeration pipe is 50 mm, and the distance from the outlet of the aeration pipe to the paddle surface is 100 mm. The sketch of the propeller and the vent pipe are shown in Figure 2.

**Table 1.** SPP-1 propeller geometric parameters [6].

**Figure 2.** The sketch of the propeller and the vent pipe.

### *2.2. Control Equation*

Under the full immersion condition, ventilation through the aeration pipe allows the paddle always being in an unsteady air-water two-phase flow field. The VOF (volume of fluid) model [13] was selected for capturing the interface between water and air. Air and water in the flow field are deemed incompressible.

Mass conservation equation can be written as:

$$\frac{\partial p}{\partial t} + \nabla(\rho \stackrel{\rightarrow}{\mu}) = \mathcal{S}\_m \tag{1}$$

The above equation is a general expression of the mass conservation equation, where the term *Sm* can be any user-defined source term added to the continuity term.

The momentum equations are given as:

$$\frac{\partial(\rho\mu\_i)}{\partial t} + \frac{\partial(\rho\mu\_i\mu\_j)}{\partial \mathbf{x}\_j} = -\frac{\partial p}{\partial \mathbf{x}\_i} + \frac{\partial}{\partial \mathbf{x}\_j} \left(\mu \frac{\partial \mu\_i}{\partial \mathbf{x}\_j} - \rho \overline{\mu\_i'\mu\_j'}\right) + \mathcal{S}\_i \tag{2}$$

where *μ<sup>i</sup>* is the velocity component in the xi direction (*i*, *j* = 1, 2, 3) of the Cartesian coordinate system, *p* the fluid pressure, *ρ* the fluid density, *μ* the dynamic viscosity coefficient, *t* the time, and *Si* the volume force.

In this paper, the *k* − *ω* SST [14,15] was employed as the turbulence model, which could take the transport characteristics of turbulent shear forces into account, and make possible to obtain more accurately results of the flow separation in the counter pressure gradient region [16].

### *2.3. Computational Domain and Grid Division*

In order to make the flow field around the propeller be fully developed, and the numerical results be accurate and reliable, the size of the flow domain in the simulation should be large enough. In practice, the entire flow domain is divided into two parts: a cylindrical rotating zone containing the paddle and an outer cylindrical stationary zone. The outer diameter of the stationary zone is 5.0*D*, the inlet is 3.5*D* from the paddle plane, and the outlet is 7.0*D* from the paddle plane. The diameter of the rotation domain is 1.2*D*, and both end surfaces of the cylinder are 0.45*D* from the paddle surface. Figure 3 shows an oblique view of the propeller's computational domain.

**Figure 3.** Oblique view of propeller and computational field.

To simulate the rotation of the paddle in the flow field, the overlapped mesh method was employed. At the overlapping domain between the stationary and the rotational domains, the information of the flow field is exchanged through the overlapping grids. The mesh in the overlapping domain should be the same size as much as possible. We set the polyhedral mesh in the inner rotational domain, while the cutting body mesh in the outer stationary domain. Moreover, local mesh densification was performed around the overlapping domain. Finally, the total number of

grids is 3.5 million, in which the rotation domain is 1.7 million, and the stationary domain is 1.8 million. The min and max size of computational grids are 7.65 × <sup>10</sup>−<sup>4</sup> mm and 11 mm, respectively. The surface mesh is shown in Figure 4.

**Figure 4.** Mesh on the propeller and aeration pipe.

### *2.4. Grid Convergence*

Let *Y*+ be the non-dimensional wall distance defined as *Y*+ ≡ *u*∗*y*/*ν*, where *u*<sup>∗</sup> is the friction velocity at the nearest wall, *y* the distance to the nearest wall, and *ν* the local kinematic viscosity of the fluid. In this sub-section the sensitivity of numerical results with respect to *Y*+ is studied. Generally, the range of dimensionless wall distance *Y*+ ≈ 30 ∼ 100 is acceptable for blade surfaces. Table 2 lists the error of the thrust coefficient *Kt* and torque coefficient 10*Kq* compared with the test value with different *Y*+ value. It can be seen that three errors reach minimum when the *Y*+ is equal to 60. Therefore, we set *Y*+ = 60 in the meshing setup.


**Table 2.** Results for the *Y*+ for selection.

Grid convergence study is important in the mesh generation process. Here four grid numbers, 1.77 million (G1), 2.5 million (G2), 3.5 million (G3), and 5 million (G4) are selected for the convergence study. The grid number and simulation results are given in Table 3. From these results, it is clearly noticed that with the growth of grid number, the results get closer to the experimental ones, i.e., the grid convergence can be guaranteed in these numbers. However, the greater the grid number, the lower computational efficiency for the numerical simulation. It is found that the accuracy was only slightly improved in the G4 case at the expense of a lot of calculation time. Hence, when determining the grid number, one should make a compromise between accuracy and efficiency. In this paper, the fine mesh of G3 was adopted in the present simulations, which achieves both computational efficiency and accuracy.


**Table 3.** Results from the grid convergence.

### *2.5. Boundary and Initial Conditions*

Due to the air-suck phenomenon, the flow around the surface-piercing paddle is usually an air-water two-phase unsteady one. The density of water and air in the simulation is 997.56 kg/m<sup>3</sup> and 1.18 kg/m3, respectively. The dynamic viscosity of water and air is 8.89 × <sup>10</sup>−<sup>4</sup> Pa·s and 1.86 × <sup>10</sup>−<sup>5</sup> Pa·s, respectively. The VOF approach was employed to track the air-water interface. The principle of this approach is to determine the interface by counting the volume ratio function of the fluid in each cell rather than the motion of the particle on the free surface.

The fluid density can be expressed as:

$$m\frac{dv}{dt} = \sum a\_m \rho\_m \tag{3}$$

where *α<sup>m</sup>* is the volume fraction of the *m*th type of fluid in each cell, which satisfies the following condition:

$$\sum\_{m=1}^{n} a\_m = 1\tag{4}$$

In this paper, velocity inlet was divided into air velocity inlet and water velocity inlet. Both inlets have the same velocity vector. The volume fraction of the overall inlet was set as the mixture of air and water fractions, where the volume fraction of air at the air velocity inlet was set to 1, and the water volume fraction was set to 0. The outlet is set as a pressure outlet and the standard atmospheric pressure is made as the reference pressure. The outer wall of the stationary zone is set as a plane of symmetry, and the surfaces of the surface-piercing paddle and the vent pipe are set to be non-slip and non-penetrable. The outer boundary of the rotation domain is set to overlapped grids.

The second-order schemes were applied for spatial discretization and linear interpolation, as well as the time integration. The time step is 6.94 × <sup>10</sup>−<sup>5</sup> s, during which the propeller rotates 1◦ under the given rotational speed. As suggested by Blocken and Gualtieri [17], this time step makes the condition CFL ≤ 1 satisfied. The numerical calculation was performed on a PC with an Intel Xeon CPU X5690 (6 cores, 3.46 GHz), and one case simulation approximately costs 40 h.

### *2.6. Validation of the Numerical Setup*

To verify the aforementioned numerical setup, the hydrodynamic performance of a propeller in the open water was investigated and compared with experimental results from Ghassemi et al. [6]. The main parameters of the SPP-1 propeller are listed in Table 1. The hydrodynamic characteristics are plotted using solid lines (see Figure 5). Let *T*, *Q*, *n*, and *VA* be the thrust, torque, rotational speed, and forward speed of the propeller, respectively. The relevant parameters for the propeller can be nondimensionalized as follows.

Speed coefficient:

$$J = \frac{V\_A}{nD} \tag{5}$$

Thrust coefficient:

$$K\_l = \frac{T}{\rho n^2 D^4} \tag{6}$$

Torque coefficient:

$$m\frac{dv}{dt} = f\tag{7}$$

Promote efficiency:

$$
\eta = \frac{K\_t}{K\_q} \frac{I}{2\pi} \tag{8}
$$

In the numerical and experimental setup, the submergence ratio for the propeller is *It* = 1/3, and the rotational speed is *n* = 2400 rpm. The speed coefficient *J* ranges from 0.85 to 1.45, which was obtained by varying the magnitude of the incoming flow velocity *VA*.

By using the CFD solver, we obtained the thrust coefficient *Kt*, the torque coefficient 10*Kq* of the propeller, and the open water efficiency *η*. 10*Kq* is used due to the fact that *Kq* is one order smaller than *Kt* and *η*. Figure 5 compares the CFD results with the experimental ones given in Table 1 [6].

**Figure 5.** Comparison of CFD and experimental results at *It* = 1/3.

As shown in Figure 5, the thrust coefficient *Kt* obtained by the CFD is slightly larger than the experimental one when *J* = 0.85–1.30, and then the discrepancy gradually decreases with the growth of *J*. The torque coefficient 10*Kq* has a similar trend but the discrepancy between CFD and experimental results is slightly larger than that in the *Kt* case. In contrast, the open water efficiency η obtained by the CFD is always slightly smaller than the experimental one when 0.85 < *J* < 1.30, but the error rapidly grows with *J* when 1.30 < *J* < 1.45 and reaches its maximum 12% at *J* = 1.45. This study demonstrates that the numerical setup in the CFD solver is effective and can provide the required precision for the hydrodynamic simulation of the propeller with the existence of air-water interface.

### **3. Numerical Results of the Ventilated SPP**

### *3.1. Calculation of Hydrodynamic Coefficient of Wigley Ship*

Figure 6 depicts the wake of the ventilated SPP at *J* = 1.15, from which it can be seen that the majority of the air stream flow bypasses the hub, while the rest of the flow is swung away by the rotating blades. Figure 6 displays a snapshot of the air volume fraction (green color) flowing around the rotating propeller. One can find that the trailing edges are enclosed by air bubbles due to the relatively low pressure around them. The air bubbles on the trailing edges generally stretch from the root to the tip of the blades, though on the leading edges the bubbles may only concentrate around the root domain of the blades. In the wake, three helical air bubbles are forming and rotating about the centric air stream. In the near propeller F field, the helical air bubbles and the centric air stream are connected to each other. With the flow moving downstream, the helical air bubbles gradually separate from the centric air stream, and then both helical air bubbles and centric air stream keep breaking into smaller bubbles until all of them disappear in the far field. Through setting up the interface and defining two-phase flow displayable distribution, it also shows the air and water distribution on the two cross sections *x*/*D* = −0.5 and *x*/*D* = −1.0 after the paddle plane. The red and blue colors refer to air and water, respectively. It can be seen that the air-content ratio gradually decreases along the downstream direction of the wake.

**Figure 6.** The wake after the ventilated propeller at *J* = 1.15 when *x*/*D* = −1.0 (left section) and *x*/*D* = −0.5 (right section).

### *3.2. Effect of Ventilation on the Hydrodynamic Performance of the SPP*

Table 4 compares the numerical results for the ventilated SPP with the unventilated one under full submergence condition.


**Table 4.** Numerical results for the SPP under full submergence conditions.

From Table 4, it can be seen that the thrust coefficient *Kt* and torque coefficient 10*Kq* of the propeller decrease after ventilation. With the growth of the speed coefficient *J* (from 0.85 to 1.30), the influence of ventilation on the thrust coefficient *Kt* gradually diminishes, while the effect on the torque coefficient 10*Kq* is reinforced. On the other hand, the efficiency *η* of the ventilated propeller is slightly less than the unventilated one under low speed coefficient *J*. However, with the growth of the speed coefficient *J*, the efficiency *η* of the ventilated propeller quickly increases as compared to the unventilated one.

As one sees, the thrust *T* and the torque *Q* decrease after ventilation. This is mainly due to the fact that the air bubbles attached to the surfaces of the blades after ventilation. However, with the growth of the speed factor, the effect of the ventilation on the thrust coefficient *Kt* is getting weaker, while the effect on the torque coefficient 10*Kq* becomes stronger, which makes the efficiency *η* of the ventilated propeller greater than the unventilated one.

### *3.3. Effect of Ventilation on the Pressure Distribution*

Figure 7 portrays the pressure distribution cloud on the blades of the unventilated SPP at *J* = 0.85. Figure 7a,b show the pressure distribution on the pressure surface and the suction surface, respectively. The arrows in the figures indicate the rotational direction of the propeller. From Figure 7a, one notes that on the pressure surface the high pressure is mainly concentrated on the region near the leading edge of the blades. Obviously, the closer the leading edge, the higher the pressure. While the closer the trailing edge, the lower the pressure. On the thick trailing edge, the pressure may even be negative. In contrast, as shown in Figure 7b, the lower pressure mainly locates on the region near the leading edge of the blades.

**Figure 7.** The pressure distribution on the blades at *J* = 0.85 under unventilated condition. (**a**) The pressure distribution on the pressure surface; (**b**) The pressure distribution on the suction surface.

As shown in Figure 7b, there exists a pressure jump near the leading edge, which is caused by the geometric characteristics of the SPP, i.e., the relative large pitch of the SPP induces a vortex formed after the leading edge, as shown in Figure 8. Figure 8a,b depict a non-dimensional velocity distribution cloud at the 0.5*R* section and a non-dimensional tangential velocity vector near the leading edge, respectively. The flow velocity is nondimensionalized with respect to the forward speed of the propeller *VA*. From Figure 8a,b, one can find that there exists a high flow velocity zone on the suction surface near the leading edge. The local high velocity produces a sudden pressure drop on the suction surface. Moreover, it can be seen from Figure 8a that around the root of the thick trailing edge, the flow velocity is also very high, which makes the pressure in this zone be low and thus the air bubbles can adhere to this zone.

**Figure 8.** Non-dimensional velocity distribution cloud at 0.5*R* section. (**a**) The non-dimensional velocity cloud; (**b**) The non-dimensional velocity vector cloud near the leading edge.

Figure 9 shows the pressure distribution on the surfaces of the ventilated SPP. Figure 9a,b are the pressure distribution on the pressure surface and the suction surface, respectively. Comparing Figure 9a with Figure 7a, it can be found that the pressure on the trailing edge of the pressure surface increases after ventilating due to the fact that air bubbles adhere to it. On the other hand, comparing Figure 9b with Figure 7b, it can be found that the area of the high pressure zone on the suction surface is significantly increased. This is because of the adsorption of air on the suction surface after ventilation. The increase in the area of the high pressure zone on the suction surface results in a reduction in the thrust and torque of the SPP.

Figure 10 plots the pressure distribution at *J* = 0.85 on various profiles (0.24*R*, 0.50*R* and 0.70*R*) of the ventilated/unventilated SPP. *Cp* is the pressure coefficient. Blue and red lines refer to the pressure distribution profiles of the unventilated and ventilated SPP, respectively. From Figure 10a one can see that on the profile near the root (0.24*R*), when there is no ventilation, the difference between the pressure on the pressure surface and the one on the suction surface is very large, especially near the leading edge. After ventilating, the pressure on the suction surface of the blade significantly increases and evenly distributes along the entire chord direction. The pressure almost remains unchanged because the air bubbles completely cover the root region of the blades, and in the vicinity of the leading edge (*x*/*C* ≤ 0.1), the pressure on the pressure surface and the suction surface are almost the same. This is because both the pressure surface and suction surface near the leading edge are covered by the air bubble. The same phenomenon can be observed around the trailing edge.

**Figure 9.** *Cont.*

**Figure 9.** The pressure distribution on the blades at *J* = 0.85 under ventilated condition. (**a**) The pressure distribution on the pressure surface; (**b**) The pressure distribution on the suction surface.

(**b**)

ƺ ƶ

**Figure 10.** *Cont.*

ƺ ƶ

ƺ ƶ

**Figure 10.** The pressure distribution profiles at *J* = 0.85. (**a**) The pressure distribution profile at 0.24*R*; (**b**) The pressure distribution profile at 0.50*R*; (**c**) The pressure distribution profile at 0.70*R*.

As shown in Figure 10b, on the profile 0.50*R*, the pressures on the pressure surface of the ventilated and unventilated SPP are almost the same. And only on the suction surface the pressure from the ventilated SPP is slightly bigger than the unventilated one due to the ventilation.

On the profile closer to the tip of the blade (0.70*R*), see Figure 10c, the pressure distribution on the pressure surface of the ventilated and unventilated SPP are completely the same. Only in the vicinity of the guide edge were the pressures on the suction surface of the ventilated and unventilated SPP slightly different, which is due to the non-uniformity of the flow field.

Comparing Figure 10a–c, one observes that the effect of ventilation on the pressure distribution of the SPP is limited to a zone in the vicinity of the blade root. In this zone, the ventilation has significant effect on the pressure. And along the outward radius direction, the effect of the ventilation gradually diminishes.

### **4. Conclusions**

In this paper, the effect of artificial ventilation on the hydrodynamic performance of a surface-piercing propeller (SPP) was investigated using the CFD technique. The present work is different from previous works as follows. In previous works, such as Ghassemi et al. [6,7], Young et al. [9–11], and Alimirzazadeh et al. [13], the SPPs work under half immersed conditions with the propeller naturally ventilated. In such conditions, all blades experience periodic water-entry and water-exit processes. While in the present work, the SPP works under full immersion conditions with artificial ventilation through a vent pipe in front of the propeller disc. In such condition, all blades are partially surrounded by air bubbles, in which the phenomenon as well as the hydrodynamic performance of the SPP are significantly different from those given in literatures. The main conclusions are summarized as follows:


the thrust coefficient is smaller than on the torque coefficient, which makes the efficiency (*η*) of the ventilated SPP, be significantly greater than that of unventilated SPP. The numerical results demonstrate the effectiveness of the ventilation approach for improving the hydrodynamic performance of the SPPs for high speed planning crafts.

**Author Contributions:** D.Y. performed the numerical calculations; Z.R. and Z.G. (Zhiqun Guo) completed the geometric modeling and mesh generation; Z.G. (Zhiqun Guo) and Z.G. (Zeyang Gao) analyzed the post-processing data; D.Y. wrote the paper.

**Funding:** This project is supported by the National Natural Science Foundation of China (Grant No. 51509053, No. 51579056, and No. 51579051).

**Conflicts of Interest:** The authors declare no conflict of interest. The founding sponsors had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, and in the decision to publish the results.

### **References**


© 2018 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

### *Article* **Study of Cavitation Bubble Collapse near a Wall by the Modified Lattice Boltzmann Method**

### **Yunfei Mao, Yong Peng \* and Jianmin Zhang**

State Key Laboratory of Hydraulics and Mountain River Engineering, Sichuan University, Chengdu 610065, China; maoyunfeigy@163.com (Y.M.); zhangjianmin@scu.edu.cn (J.Z.)

**\*** Correspondence: pengyongscu@foxmail.com; Tel.: +86-187-8023-3156

Received: 2 September 2018; Accepted: 3 October 2018; Published: 12 October 2018

**Abstract:** In this paper, an improved lattice Boltzmann Shan-Chen model coupled with Carnahan-Starling equation of state (C-S EOS) and the exact differential method (EDM) force scheme is used to simulate the cavitation bubble collapse in the near-wall region. First, the collapse of a single cavitation bubble in the near-wall region was simulated; the results were in good agreement with the physical experiment and the stability of the model was verified. Then the simulated model was used to simulate the collapse of two cavitation bubbles in the near-wall region. The main connection between the two cavitation bubble centre lines and the wall surface had a 45◦ angle and parallel and the evolution law of cavitation bubbles in the near-wall region is obtained. Finally, the effects of a single cavitation bubble and double cavitation bubble on the wall surface in the near-wall region are compared, which can be used to study the method to reduce the influence of cavitation on solid materials in practical engineering. The cavitation bubble collapse process under a two-dimensional pressure field is visualized, and the flow field is used to describe the morphological changes of cavitation bubble collapse in the near-wall region. The improved lattice Boltzmann Method (LBM) Shan-Chen model has many advantages in simulating cavitation problems, and will provide a reference for further simulations.

**Keywords:** collapse near a wall; double cavitation bubble; tilt distribution cavitation; parallel cavitation; pseudopotential lattice Boltzmann model

### **1. Introduction**

When a liquid is heated at a constant temperature or depressurized by static or dynamic methods at a constant temperature, steam bubbles or vapour-filled cavitation bubbles appear and develop over time. Cavitating water in a low-pressure zone involves a large amount of vapour forming a two-phase flow, and when the water flows through a region with a higher pressure downstream, the cavitation bubbles collapse under the effect of the pressure or temperature. In the collapse event, the cavitation phenomenon includes the nascent development and collapse of cavitation bubbles, which is an unsteady compressible multiphase turbulent flow phenomenon involving mass transfer between the gas and liquid phases. During the collapse of the cavitation bubble, high pressure of up to thousands of atmospheres is generated. When the collapse of the cavitation bubble occurs within a certain distance from a solid sidewall, the sidewall is subjected to continuous impact, causing fracture or fatigue damage of the material, which leads to erosion and cavitation. Corrosion is associated with a series of problems such as mechanical efficiency reduction and equipment damage in devices. Therefore, the formation and collapse of cavitation bubbles near a wall has been a focus of cavitation research.

Cavitation phenomena were first observed experimentally using a flow field. Kling and Hammitt [1] and Lauterborn [2,3] used high-speed photography techniques to study spark- and laser-induced cavitation bubble collapse processes. Lauterborn [3] studied the collapse of cavitation bubbles near a wall. Some scholars have also conducted in-depth physical research on cavitation [4–13]. However, although the bubble collapse process can be observed through this test method, data such as pressure field and velocity field are not available. With the development of numerical simulation technology, more methods have emerged to simulate the bubble collapse process, such as the interface tracking method, volume-of-fluid (VOF) method, and level set method. LBM is another new and efficient numerical simulation technology that has been developed as a digital model. Relative to traditional computational fluid dynamics (CFD), LBM has substantial advantages. Sukop and Or [14] first simulated the expansion and collapse of cavitation bubbles using the lattice Boltzmann Shan-Chen model; Chen et al. [15] used the lattice Boltzmann method to simulate two-dimensional cavitation within static and shear flows. The Rayleigh-Plesset equation results were compared with their simulation results. Shan et al. [16] and Zhou et al. [17] used LBM to simulate the growth and collapse of a single bubble under a two-dimensional pressure field and the sound field at a plane rigid wall. Kucera and Blake [18] used the mirror image method to simulate cavitation erosion of a plane rigid interface upon single bubble collapse near the interface and simulated bubbles near a rigid interface at different angles but did not consider the influence of cavitation on the flow field. Li et al. [19] used shadow photography and light deflection to investigate the motion of a single cavitation bubble core and the collapse time under different cone angles. Zhang et al. [20] investigated the three-dimensional (3D) cavitation bubble phenomenon at a low liquid pressure and successfully reproduced bubble growth in low-pressure water. Mishra et al. [21] investigated the coupling between the hydrodynamics of a collapsing cavity and the resulting solute chemical species introduced by cavitation based on the Shan-Chen multiphase model. Other scholars have carried out some numerical simulations on cavitation [22–33].

Previous simulations of cavitation have focused on the relationship between bubble collapse and the position of a single cavitation bubble relative to the boundary or the collapse of multiple cavitation bubbles, but no simulation of the density field or pressure field has been performed. In this paper, based on the Shan-Chen lattice Boltzmann method, a C-S state equation and an exact difference method that can accurately derive the external force term are coupled. The non-equilibrium extrapolation format and pressure boundary are used to simulate the collapse of a bubble in the near-wall region. First, the collapse evolution law of a single cavitation bubble is obtained, which is consistent with the physical experiment. On this basis, the collapse evolution law of two cavitation bubbles in the near-wall region is studied, and the collapse evolution law and flow field changes of two cavitation bubbles at different pressures are obtained.

### **2. Basic Principle of the LBM**

LBM is considered to be one of the most effective methods to solve multiphase flow problems. Commonly used two-dimensional models are D2Q7, D2Q9, etc., where D is the dimension of space and Q is the number of discrete velocities. In this paper, the two-dimensional case is simulated, and in order to ensure the calculation accuracy, nine discrete velocities are used, so the D2Q9 model is used for simulation. The specific arrangement is shown in Figure 1.

**Figure 1.** D2Q9 model.

*Water* **2018**, *10*, 1439

In the Boltzmann model, discrete particle distribution function *fi* is used to replace the fluid particle distribution function:

$$f\_i(\mathbf{x} + c\mathbf{c}\_i \Delta t, t + \Delta t) = f(\mathbf{x}, t) - \frac{1}{\mathbf{r}} [f\_i(\mathbf{x}, t) - f\_i^{eq}(\mathbf{x}, y)] + F\_i(\mathbf{x}, t) \tag{1}$$

where *fi*(*x*, *t*) is a single particle density distribution function, *f eq <sup>i</sup>* (*x*, *t*) is the equilibrium particle distribution function, *τ* is the relaxation time, and the kinematic viscosity is *ν* = *c*<sup>2</sup> *<sup>s</sup>*(*τ* − 0.5)Δ*t*. The discrete speed can be defined as

$$
\begin{bmatrix}
\varepsilon\_1, \varepsilon\_2, \varepsilon\_3, \varepsilon\_4, \varepsilon\_5, \varepsilon\_6, \varepsilon\_7, \varepsilon\_8, \varepsilon\_9
\end{bmatrix} = \mathbf{c} \begin{bmatrix}
0 & 1 & 0 & -1 & 0 & -1 & -1 & 1 \\
0 & 0 & 1 & 0 & -1 & 1 & -1 & -1
\end{bmatrix}' \tag{2}
$$

where *c* = Δ*x*/Δ*t* is the grid velocity and Δ*x* and Δ*t* are the grid step and time step, respectively. The lattice sound velocity is *cs* = *c*/ <sup>√</sup>3*.* In the D2Q9 model, the equilibrium distribution function can be expressed as

$$f\_i^{cq}(\mathbf{x}, t) = \omega\_i \rho [1 + \frac{c\_i u}{c\_s^2} + \frac{\left(c\_i u\right)^2}{2c\_s^4} - \frac{u^2}{2c\_s^2}],\tag{3}$$

where *u* is the fluid velocity and *ω<sup>i</sup>* is the weighting factor of the equilibrium distribution function, where *ω*<sup>1</sup> = 4/9, *ω<sup>i</sup>* = 1/9(*i* = 2, 3, 4, 5), *ω<sup>i</sup>* = 1/9(*i* = 6, 7, 8, 9). The discrete velocity weight coefficient is related to the discrete strategy. In the discretized velocity space, the macroscopic velocity and density of the fluid can be expressed as

$$\rho = \sum\_{i} f\_i(\mathbf{x}, t) \tag{4}$$

$$
\rho \mu = \sum\_{i} f\_i(x, t) e\_i. \tag{5}
$$

The single-component multiphase flow interaction force *Fi*(*x*, *t*) is modelled as follows:

$$F(X,t) = -G\psi(X,t)\sum\_{i=2}^{9} \xi\psi(X + c\_i\Delta t, t)e\_{i\star} \tag{6}$$

where *G* is the interaction strength and *ξ<sup>i</sup>* = 1/9(*i* = 2, 3, 4, 5), *ξ<sup>i</sup>* = 1/36(*i* = 6, 7, 8, 9). According to Yuan and Schaefer (2006) [34], the interaction potential can be expressed as *<sup>ψ</sup>* = 2(*<sup>p</sup>* − *<sup>ρ</sup>cs*)/*Gc*<sup>2</sup> *s*, where *cs* <sup>=</sup> <sup>√</sup>*RT* is the lattice sound velocity. Through this interaction potential, different equations of state can be applied.

The C-S state equation can be expressed as follows:

$$p = \rho RT \frac{1 + b\rho/4 + \left(b\rho/4\right)^2 - \left(b\rho/4\right)^3}{\left(1 - b\rho/4\right)^3} - a\rho^2,\tag{7}$$

where *a* = 0.4963*R*2*T*<sup>2</sup> *<sup>c</sup>* /*pc* and *b* = 0.18727*RTc*/*pc*. The critical parameters can be expressed as follows:

$$T\_c = \frac{0.3773a}{bR} \tag{8}$$

$$p\_c = \frac{0.0706}{b^2} \tag{9}$$

The forces of the fluid act all over the objects placed in it, and how to integrate these forces into the LBM is very important. An EDM with second-order accuracy is suitable. The Boltzmann equation is discretized in the velocity space, and the term in accordance with the LB equation is derived by the EDM. The external force term in Equation (1) can be expressed as

$$F\_i = f\_i^{eq}(\rho, u + F\Delta t/\rho) - f\_i^{eq}(\rho, u),\tag{10}$$

where *F* is the total interaction force between the fluid and the solid.

### **3. Physical Model**

In order to obtain the evolution of cavitation bubbles in the near-wall region, the physical models are shown in Figure 2.

**Figure 2.** Physical model. (**a**) (single bubble model); (**b**) (double bubbles model) (R0—bubble initial radius; bi—distance; Pv—vapour pressure in bubble; P∞—ambient pressure).

The figure above shows the initial layout of the simulation. Here, (a) is the physical model of the study of the evolution of a single cavitation bubble in the near-wall region. *R*<sup>0</sup> is the radius of the cavitation bubble, *b* is the distance from the centre of the cavitation bubble, *Pv* is the pressure in the cavitation bubble, and *P*∞ is the pressure outside the cavitation bubble. The left and right boundaries are infinite areas, the upper boundary is the pressure inlet, and the lower boundary is the rigid wall. (b) is the physical model for the study of the evolution of the double cavitation bubble in the near-wall region. *b*<sup>1</sup> is the distance from the centre of the left cavitation bubble to the rigid wall, *b*<sup>2</sup> is the distance from the centre of the right cavitation bubble to the rigid wall, and *b*<sup>3</sup> is the horizontal distance of the centre of the two cavitation bubbles. The other settings are the same as those of the physical model of a single cavitation bubble in the study of the collapse of the near-wall region.

### **4. Simulation Content and Parameter Initialization Settings**

In this paper, the numerical simulation of cavitation collapse in the near-wall region mainly includes three parts: First, the evolution law of the collapse of a single cavitation bubble in the near-wall region is obtained by numerical simulation, and the bubble deformation is compared with the results of classical experiments to verify the model. Second, the pressure field of the bubble is simulated and theoretically study the mechanism of the erosion of the bubble on the solid sidewall. Third, the previously validated model is used to study the evolution of the double cavitation bubble in the near-wall region. In this paper, if there is no special explanation, the unit in the text always uses the grid unit, the length unit is lu, the time unit is ts, the mass unit is mu, and the temperature unit is tu; thus, the density unit is mu·lu−<sup>3</sup> , and the pressure unit is mu·lu−<sup>1</sup> ts<sup>−</sup>2. The simulated initial layout is the same as the physical model. In the simulation, *T*/*Tc* = 0.689 is used to simulate and the equilibrium pressure *<sup>p</sup>* = 0.0028 mu·lu−<sup>1</sup> ts−<sup>2</sup> is gotten by the equal area rule is. In the C-S state of equation, *a* = 1,

*b* = 4 and *R* =1[34] are adopted. The initial temperature is set to a specific temperature, the velocity is zero, and the density field is initialized as follows:

$$\rho(\mathbf{x}, \mathbf{y}) = \left(\rho\_{\text{liquid}} + \rho\_{\text{gas}}\right) / 2 + \left(\rho\_{\text{liquid}} - \rho\_{\text{gas}}\right) / 2 \cdot \left[\tanh\left(\left(2\left(\sqrt{\left(\mathbf{x} - \mathbf{x}\_1\right)^2 + \left(\mathbf{y} - \mathbf{y}\_1\right)^2} - R\_0\right)/W\right)\right), \tag{11}$$

where x1 and y1 is the location of the middle of the bubble at the initial moment, the hyperbolic tangent function *tanh* = (*e<sup>x</sup>* − *<sup>e</sup>*−*x*)/(*e<sup>x</sup>* + *<sup>e</sup>*−*x*) and the phase interface width is *<sup>W</sup>* = 4. Since the cavitation bubble radius is approximately ten times the width of the phase interface, it can obtain better numerical stability of the model, but it is essential to ensure that the bubble collapse is not affected by other boundaries and that the calculation cost is reduced. Therefore, our simulation calculation area is 401 × 401.

Through the C-S state equation, the P-V curve of the gas-liquid isotherm curve can be obtained, and then the Maxwell construction should be used, which can be stated as

$$\int\_{V\_{\rm m,l}}^{V\_{\rm m,g}} P dV\_{\rm m} = p\_0 (V\_{\rm m,g} - V\_{\rm m,l}) \, \tag{12}$$

where *P* is the pressure in the EOS and *p*<sup>0</sup> is a constant pressure. When the equation is established, *p*<sup>0</sup> is the equilibrium pressure, and *Vm*,*g* and *Vm*,*<sup>l</sup>* are the physical quantities that characterize the equilibrium gas pressure and the equilibrium liquid pressure, respectively. In addition, the coexistence densities *ρ<sup>v</sup>* and *ρ<sup>l</sup>* of gas and liquid, respectively, can also be determined by phase separation simulation with a slight random disturbance of the initial density. In the calculation process, *ρ<sup>l</sup>* need to be slightly adjusted to ensure *ρ<sup>l</sup>* has the same density as the pressure boundary so that an additional pressure difference between the inside and the outside of the bubble is obtained after the fluid balance in the entire calculation domain.

The collapse of cavitation bubbles under the influence of a single wall surface is studied, so infinite areas are needed on both the left and right sides to ensure the computing domain is unaffected by the boundary fluctuations. Under the premise of the minimum calculation area, the unbalanced extrapolation format is well suited for our needs. Therefore, the non-equilibrium extrapolation format is used for the left and right borders. The bottom boundary uses a standard bounce-back format. In addition, a pressure boundary condition is applied at the inlet, and the liquid pressure in the calculation zone is equal to the pressure boundary pressure.

### **5. Study of the Evolution of a Single Cavitation Bubble**

Detailed experimental data have been obtained for the evolution of cavitation bubble collapse at different distances from the wall [35]. However, traditional experimental data and numerical simulation have great limitations. For example: Plesset and Chapman [36] made six assumptions, including negligible surface tension, constant vapour pressure and ambient pressure, an incompressible liquid, non-viscous flow, and no permanent gas, which are challenging to satisfy in experiments and LBM simulations, resulting in inapplicability to specific practical problems. In this paper, using the improved Shan-Chen model, numerical simulation data of cavitation bubble collapse at different locations are obtained at a specific temperature. This part of the study mainly consists of two parts: the first part verifies the simulation results by comparison with physical experiments, and the second part obtains the evolution law of the collapse of a single cavitation bubble in detail. The following is a comparative analysis of density field images of LBM numerical simulations and experimental images of cavitation bubbles collapsing at two different positions with the same radius and pressure. Through comparison with the experimental data of Philipp [35], the LBM calculation results which are shown in Figures 3 and 4 are found to be consistent with the qualitative analysis of the experimental data which are shown in Figures 5 and 6.

Figure 3 shows the change of density field during the collapse of cavitation bubbles. In this simulation, a dimensionless quantity *λ* = *b*/*R*<sup>0</sup> was introduced, which is the amount that characterizes the distance from the centre of the bubble to the wall. The figures below show cases with *λ* = 1.6,

*R*<sup>0</sup> = 80 lu and *λ* = 2.5, *R*<sup>0</sup> = 80 lu. The cavitation bubble is initially a circle. The bubble size and the thickness of the gas-liquid boundary layer are controlled by Equation (11). Since the pressure difference exists inside and outside the bubble, the bubble is deformed by extrusion. Due to the influence of the bottom rigid sidewall, the longitudinal flow is blocked, and a negative pressure forms under the bubble to induce longitudinal expansion of the bubble. Due to the shrinkage and deformation of the bubbles, the volume is decreasing, and the surrounding liquid fills the space created by the bubbles, resulting in a decrease in the density and pressure around them. Then, the pressure of the upper pressure boundary is first transmitted to the upper surface of the bubble, and a high-pressure zone is formed in the upper part of the bubble that acts together with the low-pressure zone of the cavitation bubble to form a depression at the upper portion (*t* = 470). Due to the rebound effect of the liquid and the relatively high speed of movement of the upper portion of the bubble, a relatively large conical high-pressure region is formed in the upper portion, which is crucial in the subsequent deformation. Over time, the depression continues to expand (*t* = 530), and with the influence of the surrounding high pressure, the bubble gradually shrinks, assuming a crescent shape (*t* = 570). When the sag causes the upper surface of the bubble to touch the lower surface, a large pressure difference directly breaks down the empty bubble, forming a micro shock wave, which has a destructive effect on the wall surface (*t* = 588). At the same time, a complex sound field is generated, which causes additional damage to the rigid wall. By comparison, the morphology of cavitation bubble collapse differs for different dimensionless parameters *λ* and Δ*P* is confirmed. In our simulation, a crescent-shaped bubble is formed at *λ* = 1.6 and is then broken down to form two bubbles, which generates a micro-shock; however, at *λ* = 2.5, the bubble is squashed directly. A crescent-shaped bubble is formed, but the bubble does not break in the middle and finally collapses in the form of a small bubble. The calculation results are consistent with the results of Philipp [35]. When the bubble is too far from the wall surface, the sidewall has a small retarding effect on the bubble, and a strong negative pressure is unlikely to form under the bubble, so when the bubble collapses, the bubbles are gradually crushed and collapsed, and the impact of the formed pressure on the wall surface is also alleviated.

**Figure 3.** Density field of the LBM simulation: λ = 2.5, *R*<sup>0</sup> = 80 lu.

**Figure 4.** Density field of the LBM simulation: λ = 1.6, *R*<sup>0</sup> = 80 lu.

**Figure 5.** Experimental images: <sup>λ</sup> = 2.5, *<sup>R</sup>*<sup>0</sup> = 1.45 <sup>×</sup>10−<sup>3</sup> m.

**Figure 6.** Experimental images: <sup>λ</sup> = 1.6, *<sup>R</sup>*<sup>0</sup> = 1.45 <sup>×</sup>10−<sup>3</sup> *<sup>m</sup>*.

### **6. Study of the Evolution of a Double Cavitation Bubble**

In the previous section, the collapse evolution of a single cavitation bubble in the near-wall region was discussed, which was found to be in good agreement with the physical experiment, thus verifying the stability of the numerical simulation model. However, in actual engineering, cavitation bubbles do not appear alone. Cavitation clouds typically form in cavitation-concentrated areas, and hundreds of cavitation bubbles interact and collapse under the extra pressure; thus, the erosion of the wall is more complex. This research area merits further study and is instructive for the possibility of reducing cavitation damage. In our study, using the previously validated mathematical model, the collapse evolution law by simulating the collapse of two cavitation bubbles in the near-wall region under pressure induction was obtained. Taking two of these cases as examples, the evolution of the two cavitation bubbles in the near-wall region under pressure induction was analysed. For convenience of description, the left and right bubbles are designated left bubble (LB) and right bubble (RB), respectively. The spatial direction is set as shown in Figure 1. For example, the upper right is the e6 direction. Similar to the previous simulation, the dimensionless quantities *λ*1, *λ*2, and *λ*<sup>3</sup> is introduced, where *λ*<sup>1</sup> = *b*1/*R*0, which is the amount that characterizes the distance of the LB centre from the wall; *λ*<sup>2</sup> = *b*2/*R*0, which is the amount that characterizes the distance of the LB centre from the wall surface; and *λ*<sup>3</sup> = *b*3/*R*0, which is the amount that represents the horizontal distance between the LB and RB centres.

### *6.1. Case 1: Numerical Simulation of Tilt Distribution Cavitation with Two Cavitation Bubbles*

When *λ*<sup>1</sup> = 1.2, *λ*<sup>2</sup> = 2.7, and *λ*<sup>3</sup> = 1.5, the line connecting the centres of the two cavitation bubbles is at an angle of 45◦ to the wall surface. The density field and the pressure field coupling with the velocity field of the cavitation bubble collapse process are shown in Figures 7 and 8, respectively.

**Figure 7.** Density field of cavitation bubble collapse.

**Figure 8.** Pressure field and velocity field of cavitation bubble collapse.

The cavitation bubble distribution at the initial moment is as shown by *t* = 1 ts, and the specific parameter settings are the same as the previous settings. The collapse process of each cavitation bubble was analysed. When the simulation proceeds to *t* = 200 ts, the change of LB is similar to the collapse of a single cavitation bubble. Due to the blockage of the wall surface, a low-pressure zone is formed in the bottom region, and LB exhibits an elongated deformation state in the e5 direction. The e6 direction of the LB and the e8 direction of the RB are attracted to each other, and the opposite direction is deformed. If the gas density and the liquid density are in equilibrium at this time, the two cavitation bubbles would attract each other and eventually merge into a single bubble. However, since the liquid has been pressurized during the simulated initialization, the additional generated pressure prevents the emergence of the two bubbles. As the simulation progresses, the deformation is further aggravated, and the bubbles decrease under the action of additional pressure. Between them, the change of LB is the most easily detected. When the time reaches *t* = 350 ts, the LB is elongated and deformed due to the low pressure resulting from the blockage of the wall surface; the e6 direction of RB appears the same as the change of a single cavitation bubble in the near-wall region, and a depression occurs. The change becomes more apparent at 400 ts. Around the shrinking bubble, the liquid fills the newly available space, and a slightly lower pressure appears around the cavitation bubble. The transition can be observed in the velocity field. Around the bubble, the velocity direction tends to support cavitation. Note that in the figure of the pressure field and the velocity field, for the sake of convenience, only the velocity of the liquid is shown and we do not plot the velocity inside the bubble. The LB exhibits an elliptical shape that exhibits an inclination in the e6–e8 direction under the action of the sidewall, the attraction of RB and the external extreme pressure. The RB is also elliptical due to the attraction of the two bubbles, but the deformation is not as strong as with the bubble on the left. At this time, relative with the collapse rule of a single cavitation bubble in the near-wall region, LB is equivalent

to rigidity avoidance for RB, so a low-pressure zone is generated between the two bubbles, and high pressure is generated by the upper pressure boundary. A pressure difference is generated in the e6–e8 direction of the RB caused the deformation to continue to produce. Inspection of the velocity field reveals that the maximum velocity occurs in the same direction of the RB. At *t* = 450 ts, the cavitation bubble is crescent-shaped. The bubble is strongly compressed from the upper part of the depression until the entire bubble is completely collapsed from the upper part of the bubble, which occurs because the flow field is not completely symmetrical on both sides of the e6–e8 direction. The upper part is first crushed by the pressure transferred from the pressure boundary. The other reason is that the most attractive part of LB is located at the nearest position of the two bubbles. The velocity field indicates that the maximum velocity direction has been deflected, and RB begins to collapse from top to bottom under the pressure difference, producing a huge jet and complex sound field. The huge pressure generated by the collapse acts in conjunction with other factors to promote the continued collapse of LB. At *t* = 550 ts, the pressure field and velocity field images indicate that the flow field changes due to the collapse pressure of RB, and the micro-jet generated by RB in the e6 direction collides with the downward flowing liquid to generate a high voltage and noise, and energy begins to be consumed; furthermore, the micro-jet in the e7 direction overlaps with the original flow field, and the flow state becomes complex. As the simulation proceeds, the flow field exhibits vortices, and LB begins to collapse. First, the cavitation bubble e6 direction begins to shrink under the influence of the pressure of RB, which causes a protrusion in the upper portion of the cavitation bubble. At this time, the energy generated by the RB collapse is insufficient to continue to compress LB. Under the joint action of the upper pressure boundary pressure and the high pressure generated by the RB collapse, a high-pressure region (*t* = 550 ts) is formed at the upper convex portion, and collapse from the upper portion of the cavitation bubble is induced. As RB is destroyed, the e9 direction finally collapses, and a jet from the e7 to the e9 direction is generated in the flow field, so that the entire flow field exhibits vortices. Finally, LB collapses under the combined action of various factors, and the generated pressure impacts the wall surface to avoid impact. The pressure field diagram of *t* = 700 ts indicates that a high pressure is generated on the wall surface. Note that after RB collapse, LB collapse occurs under the influence of various factors, and the mechanism of action is relatively complex. The LB high-pressure zone rotates around the cavitation bubble. This visualization is a powerful way to illustrate the complex vortices and other phenomena in the flow field, involving other physical quantities. The parameters of this study cannot describe the more detailed characteristics. The specific mechanism needs further study.

### *6.2. Case 2: Numerical Simulation of Two Parallel Cavitation Bubbles*

When *λ*<sup>1</sup> = 1.2, *λ*<sup>2</sup> = 1.2, and *λ*<sup>3</sup> = 2.4, that is, the line connecting the centres of the two cavitation bubbles is parallel to the wall surface, and the density field of the cavitation bubble collapse process is shown in Figure 9 as well as the pressure field and the velocity field shown in Figure 10.

**Figure 9.** Density field of cavitation bubble collapse.

**Figure 10.** Pressure field and velocity field of cavitation bubble collapse (800 ts is the pressure field and streamline diagram).

The figure shows the change in the two cavitation bubbles in the near-wall region and the initial boundary and density settings are the same as described above. Since the changes of the two bubbles are basically the same, the LB is used as an example for analysis. At *t* = 100 ts, low pressure is created at the bottom of the cavitation bubbles due to the retardation of the rigid wall. At the same time, the two cavitation bubbles are attracted to each other. The cavitation bubbles that are not attracted to the rigid wall and the cavitation bubbles begin to shrink under the action of the additional pressure, and the new space that is generated from the reduced area was filled by the surrounding liquid. Through the analysis of the pressure field, the conclusion that the pressure decreases in the e9 direction of LB can be obtained. The bottom region of the symmetry axis is defined as the pressure change region, which is the region where the pressure change is most obvious during the cavitation process, except for the pressure region generated by the collapse of the cavitation bubble itself. In the simulation, under the condition of the specific additional pressure and the relative position of the cavitation bubble, a new cavitation bubble is generated in the pressure change zone, but it will collapse quickly. Although the new cavitation bubbles are produced for a short period of time, the effect on the flow field is very important and will be explained in detail in the analysis of the maximum wall pressure behind. When the simulation is carried out to *t* = 200 ts, for LB, the high pressure generated by the collapse of the new cavitation bubble in the pressure change zone continues to act on the e9 direction of the cavitation

bubble, showing a slightly lifted shape. As the simulation progresses to *t* = 400 ts, the deformation continues; the deformation of the cavitation bubbles is already evident, and the cavitation bubbles are inclined at an inclination angle of 45◦, and the cavitation bubbles become crescent-shaped. This is mainly due to the interaction between the cavitation bubbles, the interaction between the cavitation bubbles and the wall surface, and the collapse of the new cavitation bubbles in the pressure change zone. In addition, the velocity field indicates that the velocity of the LB in the e7 direction is the largest. At *t* = 600 ts, the LB appears asymmetric on the axis in the e7–e9 direction. The upper part is strongly influenced by the pressure boundary and the lower part is strongly influenced by the rigid wall. Therefore, the position where the pressure difference between the upper side and the lower side of the LB is the largest is located at the centre, so that the cavitation bubbles collapse from the upper portion until they completely collapse. For each bubble, a complex eddy current phenomenon occurs in the flow field because the collapse does not occur instantaneously but spreads from the top to the adjacent region. The 800 ts streamline diagram clearly shows the eddy currents in the flow.

### **7. Maximum Wall Pressure**

The maximum wall pressure of a single cavitation bubble and two parallel distributed cavitation bubbles collapsed in the near-wall region was identified in this study. The initial conditions are set as follows: For the case of a single cavitation bubble, the value of *λ* (from 1.05–2.0) with a series of gradients is used to simulate the maximum wall pressure resulting from single cavitation bubble collapse under different initial conditions. Since the maximum pressure of a single cavitation bubble on the wall is only likely to occur below the cavitation bubble (e5 direction), the pressure of the wall below the cavitation bubble is recorded to find the maximum value. For the case of a single cavitation bubble, a lambda value (from 1.05–2.0) with a series of gradients is used to simulate the maximum wall pressure resulting from a single cavitation bubble collapse under different initial conditions. Since the maximum pressure of a single cavitation bubble on the wall is only likely to occur below the cavitation bubble (e5 direction), the pressure of the wall below the cavitation bubble is recorded to find the maximum value. For the case of two cavitation bubbles, *λ*<sup>3</sup> = 2.4 is used and, likewise, the value of λ (from 1.05 to 2.0) with a series of gradients was used to simulate the maximum wall pressure resulting from the collapse of two cavitation bubbles under different initial conditions, where *λ* = *λ*<sup>1</sup> = *λ*2. Since the cavitation bubble collapse process is complicated, the maximum wall pressure is not generated at the same position for the collapse process under different initial conditions. The wall pressure below the cavitation bubble (e5 direction) and the pressure change zone is recorded to find the maximum value The comparison of the maximum wall pressure under different initial conditions is shown in Figure 11.

**Figure 11.** Comparison of the maximum wall pressure under different initial conditions.

Some conclusions are given by analysis:


In summary, the maximum wall pressure generated by the collapse of the two cavitation bubbles in the near-wall region is smaller than the single cavitation bubble at *λ* < 1.7; the closer the bubble is to the wall, the greater the influence on the wall surface. When *λ* > 1.7, the distance between the cavitation bubble pair and the wall surface has little effect on the maximum wall pressure. This is a point worthy of further study with different initial conditions. It can be used as a reference for reducing cavitation damage in practical engineering.

### **8. Conclusions**

Based on the improved lattice Boltzmann Shan-Chen model coupled with C-S EOS, the collapse phenomenon of cavitation bubbles in the near-wall region is studied. In this paper, the collapse of a single bubble in the near-wall region is simulated and compared with the physical experiment; the results are in good agreement. Then, using the verified mathematical model, the simulation of the collapse of two cavitation bubbles in the near-wall region under different pressure-induced conditions is carried out, and the collapse evolution law is obtained. And analysed the effect of wall pressure on different initial conditions. By analysing in detail the density field changes and pressure field changes during bubble collapse, the following conclusions can be drawn:

1. For the case where a single cavitation bubble collapses in the near-wall region, when λ = 1.6, a crescent-shaped bubble is formed that is broken down to form two bubbles, and a micro-shock is generated; when *λ* = 2.5, the bubbles are crushed to form crescent-shaped bubbles, but the

bubbles do not break in the middle but rather ultimately collapse in the form of a small bubble. The shape of the cavitation bubble is related to the distance of the cavitation bubble from the rigid wall.


The results indicate that the improved LBM Shan-Chen model has many advantages in simulating cavitation problems, providing a reference for further simulation.

All the abbreviations are explained in the Appendix A (Table A1).

**Author Contributions:** Conceptualization, Y.P. and J.Z.; Methodology, Y.P. and J.Z.; Software, Y.M.; Investigation, Y.M.; Data Curation, Y.M.; Writing-Original Draft Preparation, Y.P. and J.Z.; Writing-Review & Editing, Y.M.; Visualization, Y.M.; Funding Acquisition, Y.P.

**Funding:** This research was funded by the National Natural Science Foundation of China (51579166) and the National Key Research and Development Program of China (2016YFC0401705).

**Conflicts of Interest:** The authors declare no conflict of interest.

### **Appendix A**


**Table A1.** All variables with definitions.


**Table A1.** *Cont.*

### **References**


© 2018 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

### *Article* **Three-Dimensional Aerators: Characteristics of the Air Bubbles**

### **Shuai Li 1,2, Jianmin Zhang 3,\*, Xiaoqing Chen 1,2 and Jiangang Chen 1,2**


Received: 26 August 2018; Accepted: 29 September 2018; Published: 12 October 2018

**Abstract:** Three-dimensional aerators are often used in hydraulic structures to prevent cavitation damage via enhanced air entrainment. However, the mechanisms of aeration and bubble dispersion along the developing shear flow region on such aerators remain unclear. A double-tip conductivity probe is employed in present experimental study to investigate the air concentration, bubble count rate, and bubble size downstream of a three-dimensional aerator involving various approach-flow features and geometric parameters. The results show that the cross-sectional distribution of the air bubble frequency is in accordance with the Gaussian distribution, and the relationship between the air concentration and bubble frequency obeys a quasi-parabolic law. The air bubble frequency reaches an apex at an air concentration (*C*) of approximately 50% and decreases to zero as *C* = 0% and *C* = 100%. The relative location of the air-bubble frequency apex is 0.210, 0.326 and 0.283 times the thickness of the layers at the upper, lower and side nappes, respectively. The air bubble chord length decreases gradually from the air water interface to the core area. The air concentration increases exponentially with the bubble chord length. The air bubble frequency distributions can be fit well using a "modified" gamma distribution function.

**Keywords:** three-dimensional aerator; air concentration; air bubble frequency; air bubble chord length

### **1. Introduction**

Cavitation erosion caused by high velocity flows is a common phenomenon in spillways or chutes with high-head dams. A commonly adopted counter-measure to reduce or prevent cavitation damage is aeration [1–6]. Three-dimensional aerators are used to enhance air entertainment into water and to prevent cavitation erosion on spillways. Full interfacial aeration is commonly observed at the free surface of the upper, lower and side nappe along the free jet (Figure 1). The effect of air entrainment is based on the jet disintegration process due to turbulence and secondary interactions with the surrounding atmosphere.

**Figure 1.** Air water flow in three-dimensional aerators: (**a**) Three Gorges Project; (**b**) Xiang Jia-Ba Project.

Air entrainment at the interfacial area is a fundamental parameter in numerical models for simulating two-phase flows. Likewise, it is also an essential consideration when designing optimum aerators. During the last few decades, numerous studies have been conducted to investigate the air concentration and free-surface aeration on an aerator. A series of experiments were conducted to study the free-surface aeration and diffusion characteristics of the bottom aeration devices [7–11]. Results indicate that the quantities and geometric scales of air bubbles are important parameters in characterizing the air concentration in the interfacial area. Chanson [12] described the flow structure of the developing aerated region in a flat chute and found that the maximum mean air content in different cross-sections was 12%. Toombes [13] showed that the maximum bubble count rates typically correspond to air concentration between 40% and 60%. The distribution of velocity indices at the different nappes, which stays at the downstream of an expanding chute aerator, were found to be different [14] and a solution was also found to improve both the bottom and lateral cavities' length [15]. The longitudinal position of the upper trajectory apex was found to be further downstream than that of the lower trajectory apex [16]. Kramer and Hager [17] found that both the entrainment rate and dimensionless bubble count rate increased with jet length. The mean air concentration is affected by aerators, but the air concentration at the bottom is limited [18]. Furthermore, Toombes and Chanson [19] studied flows past a backward-facing step, where the void fraction distribution, bubble count rate, local air distribution and water chord size were indicated. The results showed that the relationship between bubble count rate and void fraction can be defined by a quasi-parabolic function, while probability density functions of local chord size exhibited a quasi-log-normal shape. Some researchers [20–24] also established analytical equations to predict the air concentration distributions at the aerators' downstream. Numerical models such as the lattice Boltzmann method were found to successfully simulate two-phase flows and free-surface flows [25–27]. These models may be used to complement physical models and experimental measurements after careful validation and verification of model parameters such as the void fraction distribution, velocity, and air-bubble chord sizes. Chanson [28] reviewed the recent progress in turbulent free-surface flows and the mechanics of aerated flows and highlighted that physical modelling tests are efficient methods to validate phenomenological, theoretical and numerical models.

Although numerous studies regarding air–water flows have been conducted over the years, only a few researchers, however, have systematically investigated the air bubble count rate and bubble size of high velocity air–water flows downstream of the pressure outlet joined by a three-dimensional aerator (sudden vertical drop and lateral enlargement) in the tunnel. There is also a lack of knowledge regarding the influence of the bubble count rate and bubble size on the air concentration in two-phase flows, which is necessary for improving turbulence models. Thus, further insight is needed regarding the details of air entrainment in high-speed flows and the related core area features. In this study, a series of experimental investigations were systematically conducted to further the knowledge of

air water flows among the upper, lower and side nappes downstream of a three-dimensional aerator (Figure 2). The experimental data and equations obtained from this study provide novel information and a quantitative description of the aeration flow structure at the high-speed interfacial area.

**Figure 2.** Schematic with the relevant parameters: (**a**) side view and (**b**) plan view.

### **2. Experimental Setup**

Experiments were performed at the State Key Laboratory of Hydraulic and Mountain River Engineering, Sichuan University, China. The experimental setup (Figure 3) includes an upstream reservoir, a pressure section, a sudden vertical drop and lateral enlargement aerator, a free flow section, a tail-water section, and an underground reservoir. Water is discharged (can be up to 0.5 m3/s) by a continuous and stable water supply system (2.5 m wide, 3.5 m long, and 6.3 m high). The components of the pressure inlet were assembled with a smooth steel plate. The pressure section is a rectangular pipe (0.25 m wide, 0.15 high and 2.0 m long) with a variable inclination angle varied with the downstream chute. The sudden fall-expansion aerator (height *h* × width *b*) was installed at the end of the pressure section. The free-flow section was fabricated from polymethyl methacrylate (PMMA) with a roughness height of 0.008 mm. The upper, lower, and side nappes were exposed to atmospheric pressure. The constant water levels corresponding to a particular head were maintained in the upstream reservoir.

Experiments were carried out for mean velocities (*V*0) of 6 m/s, 7 m/s, 8 m/s, and 9 m/s; chute slopes of 0%, 10%, and 25%; and aerator sizes (*h*, *b*) of 0.025m, 0.045m, and 0.065 m. The test cases are summarized in Table 1 (0.9 < *q*<sup>w</sup> < 1.35, 0.89 < *R*<sup>e</sup> < 1.34). It should be noted that some combinations (series 7–10, and 13 in Table 1) pertain to cases involving three-dimensional (vertical drop and lateral enlargement) aerators. The other cases involve only either vertical drop aerators or lateral enlargement aerators. Upstream flows were supercritical (4.95 < *Fr* < 7.42) for all of the investigated flow conditions. Downstream of the aerator, a free jet, with the use of a ventilated air cavity below and/or beside it, was formed.

**Figure 3.** Sketch of the experimental model and test equipment.


**Table 1.** Summary of the operating conditions for the experiments.

The air concentration, bubble count rate, and chord length were recorded using a CQY—Z8a measurement instrument [14,29,30] with a double-tip conductivity probe (Figure 4). Different voltage indices between air and water were measured by the platinum tip. The probe consisted of two identical tips, including an external stainless steel electrode with 0.7-mm-diameter and an internal concentric platinum electrode with 0.1-mm-diameter. The tips are aligned in the flow direction and the distance between the two tips is 12.89 mm. Both tips were connected to an electronics device with a response time less than 10 μs. The vertical translation of the probe was dominated by a fixed device with an accuracy of 1 mm. The probe measurements were taken at regular intervals of 5 mm along the vertical direction from the chute bottom to the free surface (or side nappe surface) and 100 mm along the flow direction. At each location, signals were recorded at a scan rate (*f*) of 100 kHz per channel for a scan period (*T*) of 10 s. The scan period was based on the findings of Toombes [13] which indicated that a scanning period of 10 s was long enough to provide a reasonable representation of the flow characteristics while maintaining realistic time and data storage constraints.

Figure 4 presents a typical aeration structure detected by a probe at a fixed location along the aeration flow. Assuming each segment is either air or water, the air concentration can be presented as the probability of any discrete element being air. The air concentration is computed as the encountering air's probability at the leading tip of the probe.

**Figure 4.** Sketch of the double-tip conductive probe.

The number of data (*N*) can be obtained by multiplying the frequency (*f*) of the measurement instrument with the sampling time *T*, i.e., *N* = *f* × *T*. Those data were classified according to two categories, air and water signals. The air concentration (*C*) was computed as the encountering air's probability at the leading tip of the probe, which can be expressed as [26]

$$C = \frac{\sum\_{i=1}^{N} R\_i}{N} \tag{1}$$

where *Ri* represents the ratio of bubbles measured all throughout the whole measurement time. When the probe tip is completely immersed in air bubbles *Ri* = 1, otherwise *Ri* = 0.

The length of bubble chord is defined as the straight distance between the two intersections of the interface [12]. Note that an air bubble is defined as a volume of air (i.e., air entity), which can be detected by the leading tip of the probe between two continuous air–water interface events. The bubble chord length is calculated by

$$d\_i = \frac{v \times n\_i}{f} \tag{2}$$

where *ni* is the number of detected bubbles recorded; and *v* is the bubble velocity equal to the local mean aeration velocity (i.e., no slip between the air and water phases).

### **3. Results**

The air bubble frequency and chord length are two important parameters used to characterize the air concentration of aerated flows. We first describe air bubble frequency among the upper, lower and side air–water mixed layers in the experimental flows with various conditions. In this section, the relationships between air concentration with air bubble frequency, and with the relative location at which the maximum air bubble frequency are also discussed. Next, we identify the effects of initial flow velocity on air bubble chord length. In addition, the relationships between the air concentration and air bubble chord length, and between air bubble count rate and air bubble size are also evaluated.

### *3.1. Air Bubble Frequency*

The air bubble frequency characterizes the flow fragmentation, which is proportional to the specific area of the interface between two phases. The air bubble frequency, or air bubble count rate, is a function of bubbles' shape and size, surface tension, and shear forces of the fluid. Predicting how air concentration affects the bubble frequency is a complex problem [19]. A simplified analogy would

be to consider on aeration flows past a fixed probe, as a series of discrete one-dimensional air and water elements (Figure 4). Here, the air bubble frequency (*Fa*) is defined as the number of air-structures (*Na*) per unit time (*t*) detected by the leading tip of the probe sensor (i.e., *Fa* = *Na*/*t*).

Typical air bubble frequency distributions at each cross-section along the upper, lower, and side nappes are presented in Figure 5. The air bubble frequency distributions exhibit unitary self-similarity along the thickness of the mixed layer and the horizontal distance. The trend of the cross-sectional distribution tends to initially increase and then decrease from the aeration interface to the inside of the water following a Gaussian distribution that is defined by the following functions:

$$\frac{F\_a}{\left(F\_a\right)\_{\text{max}}} = \begin{cases} \frac{1}{\sqrt{2\pi}\sigma\_0} e^{-\left(\frac{Z-Z\_0}{Z\_{00}} - \mu\_0\right)^2} & \text{the upper nappe} \\\\ \frac{1}{\sqrt{2\pi}\sigma\_0} e^{-\left(\frac{Z-Z\_0}{Z\_2-Z\_{00}} - \mu\_0\right)^2} & \text{the lower nappe} \\\\ \frac{1}{\sqrt{2\pi}\sigma\_0} e^{-\left(\frac{Y-Y\_{00}}{Z\_2-Y\_{00}} - \mu\_0\right)^2} & \text{the side nappe} \end{cases} \tag{3}$$

where, *Fa*/(*Fa*)max is the dimensionless air bubble frequency and (*Fa*)max is the maximum bubble frequency in the cross-section; The characteristic value of *μ*<sup>0</sup> reflects the depth corresponding to the air bubble frequency; σ<sup>0</sup> represents the degree of deviation between the relative depth and its mean value.

**Figure 5.** *Cont*.

**Figure 5.** Air bubble frequency distributions compared with the theoretical distribution obtained from Equation (3). Upper nappe (**a**) (*h*, *b*) = (4.5 cm, 4.5 cm), *θ* = 25%. (**b**) (*h*, *b*) = (6.5 cm, 6.5 cm), *θ* = 10%; lower nappe (**c**) (*h*, *b*) = (4.5 cm, 4.5 cm), *θ* = 25%. (**d**) (*h*, *b*) = (6.5 cm, 6.5 cm), *θ* = 10%; side nappe (**e**) (*h*, *b*) = (4.5 cm, 4.5 cm), *θ* = 25%. (**f**) (*h*, *b*) = (6.5 cm, 6.5 cm), *θ* = 10%.

For all experimental results, the air bubble frequency distributions correspond well with the data reported by Toombes and Chanson [19], as illustrated in Figure 5. Figure 6 presents a plot of experimental data, *Fa*/(*Fa*)max, with calculated non-dimensional air bubble frequency using Equation (3). The calculated results at upper, lower, and lateral nappes were in qualitative agreement with experimental data, although there was some scatter. It is noted that the characteristic values of *μ*<sup>0</sup> and *σ*<sup>0</sup> differed among the upper, lower, and side nappes (Table 2). For the upper nappe, the air bubbles do not easily diffuse and instead transport to the inside of the water due to the buoyancy. This result in the air bubbles being drawn closer to the free surface (*μ*<sup>0</sup> = 0.210). For the side nappe, the air bubbles can easily diffuse and transport to the interior of the water (*μ*<sup>0</sup> = 0.283) because of the effect of large transverse turbulent diffusion [15]. For the lower nappe, the air bubbles are dragged by the buoyant force to the interior of the water (*μ*<sup>0</sup> = 0.326).

**Figure 6.** Comparison of Equation (3) with the experimental data. (**a**) upper nappe; (**b**) lower nappe; (**c**) side nappe.

**Table 2.** Values of *μ*<sup>0</sup> and σ0.


Air concentration distribution is just the external manifestation of air–water flows, whereas the internal factors include the count rate and size of the air bubbles. So, it is crucial to understand the features of the air bubble frequency and the air bubble chord length.

The air bubble frequency distributions at various positions along the jet obtained in this study and those reported by Chanson [12] and Toombes and Chanson [19] are shown in Figure 7. It is evident that the air bubble frequency initially increases and then decreases with the increase of air concentrations in the upper, lower, and side nappes. At each cross-section, the air bubble frequency profiles reach to an apex which corresponds to air concentrations of approximately 50%. The profiles tend to zero at extremely low and extremely high air concentrations. Overall, the distributions of the dimensionless air bubble frequency can be well fitted by the parabolic function:

$$\frac{F\_a}{(F\_a)\_{\text{max}}} = 4\mathbb{C}(1-\mathbb{C})\tag{4}$$

Qualitatively, the calculated results of Equation (4) correspond with the results of Chanson [12] and Toombes and Chanson [19] (Figure 7). It is worth noting that the data [12] slightly deviate from the fully developed supercritical flows, suggesting that Equation (4) has broad applicability.

**Figure 7.** *Cont*.

**Figure 7.** Dimensionless air bubble frequency as a function of the air concentration fitted with Equation (4). Upper nappe (**a**) (*h*, *b*) = (4.5 cm, 4.5 cm), *θ* = 25%. (**b**) (*h*, *b*) = (6.5 cm, 6.5 cm), *θ* = 10%; lower nappe (**c**) (*h*, *b*) = (4.5 cm, 4.5 cm), *θ* = 25%. (**d**) (*h*, *b*) = (6.5 cm, 6.5 cm), *θ* = 10%; side nappe (**e**) (*h*, *b*) = (4.5 cm, 4.5 cm), *θ* = 25%. (**f**) (*h*, *b*) = (6.5 cm, 6.5 cm), *θ* = 10%.

The majority of the data fall within the ±20% error lines as shown in Figure 8. The calculated non-dimensional air bubble frequency using Equation (4) are in reasonably good agreement with the experimental data at upper, lower, and lateral nappes.

**Figure 8.** Comparison of experimental data with calculated values of dimensionless air bubble frequency. (**a**) upper nappe; (**b**) lower nappe; (**c**) side nappe.

The relationships between air concentration and the relative distance, which corresponds to the maximum air bubble frequency are presented in Figure 9 where a large fluctuation of air concentration (*C* = 30% ∼ 60%) in the initial regime of aeration can be observed. The fluctuation may be due to the aeration layer located close to the pressure outlet, which may be thin and unstable. The air–water layer becomes stable with the development of the aeration further downstream. The air concentration corresponding to the maximum air bubble frequency gradually moves toward to 50% line. The result is consistent with the findings of Figure 7. Actually, the maximum air bubble frequency does not always coincide with *C* = 50% accurately. Factors such as the average size and length scales of discrete air and water elements may affect the local air concentration and flow conditions in air–water flows [13,31].

The relative location corresponding to the maximum air bubble frequency is shown in Figure 10. It can be seen that the relative location of the maximum air bubble frequency is 0.210, 0.326 and 0.283 times the thicknesses of the air–water layers in the upper, lower and side nappes, respectively. The results are consistent with the results of Figure 5. For the upper nappe, the air can easily be

drawn into the water because the air–water external interface directly interacts with the atmosphere. However, the air bubbles cannot easily diffuse and transport to the inside of the jet because of the effect of buoyancy. For the lower nappe, the air bubbles can easily diffuse and move into the interior of the jet due to the local air rotation and eddy currents, which are driven by the air–water interfacial turbulence coupled with the positive effect of buoyancy.

**Figure 9.** Air concentration corresponding to the maximum air bubble frequency. Upper nappe (**a**) (*h*, *b*) = (4.5 cm, 4.5 cm). (**b**) *θ* = 10%; lower nappe (**c**) (*h*, *b*) = (4.5 cm, 4.5 cm). (**d**) *θ* = 10%; side nappe (**e**) (*h*, *b*) = (4.5 cm, 4.5 cm). (**f**) *θ* = 10%.

**Figure 10.** Relative location corresponding to the maximum dimensionless air bubble frequency. Upper nappe (**a**) (*h*, *b*) = (4.5 cm, 4.5 cm). (**b**) *θ* = 10%; lower nappe (**c**) (*h*, *b*) = (4.5 cm, 4.5 cm). (**d**) *θ* = 10%; side nappe (**e**) (*h*, *b*) = (4.5 cm, 4.5 cm). (**f**) *θ* = 10%.

### *3.2. Air Bubble Chord Length Distributions*

Air bubble size is another parameter reflecting the characteristics of air–water flows. It is difficult to measure the size of air bubbles since their shapes vary greatly, are complex and highly changeable. This value can only be indirectly determined from the bubble chord length which is obtained using the double-tip conductivity probe. Here, the chord length (*chab*)mean is adopted to assess the size of the air bubbles, which is defined as follows:

$$(ch\_{ab})\_{\text{mean}} = \frac{\sum\_{j=1}^{N} (ch)\_j n\_j}{\sum\_{j=1}^{N} n\_j} \tag{5}$$

where *nj* is the count of air bubbles in the bubble chord length interval Δ(*chab*)*<sup>j</sup>* (Δ(*chab*)*<sup>j</sup>* = 0.1 mm in the experiments); and (*chab*)*<sup>j</sup>* is the average value of the chord length interval (e.g., the probability of chord length from 1.0 to 1.1 mm is represented by the label 1.05 mm).

The air bubble chord lengths (*chab*)mean are presented in Figure 11 at various positions in the upper, lower and side nappes. It is note that the air bubble chord length distributions obtained at the upper, lower, and side nappes have similar shapes. The air bubble chord length in the vicinity of the air–water interface is largest and then gradually decreases toward the inside of the water. A broad spectrum of air bubble chord lengths, i.e., from less than 0.1 mm to greater than 100 mm, is observed at each location and cross-section. At high air concentration (i.e., *C* ≥ 90%), chord length of many air bubbles can reach up to 100 mm or even more. These large values may be large air packets and air volumes surrounding the water structure (e.g., droplets). The reason is that the regime with high air concentration is always close to the air–water interfaces, which become uneven under the influence of turbulent forces. Air in the vicinity of the free surface that are trapped by the generated waves were treated as air bubbles by the conductivity probe. In addition, the shape of the air bubble is not completely spherical, in most cases, it is oval or banding, causing it to be measured to be much larger than it actually is. Indeed, it is nearly impossible to distinguish between air-bubbles and an un-enclosed bubble structure. The results indicate that the large air bubbles are constantly sheared and tore as they move towards the interior of the water. This results in air bubbles located deeper inside the body of water to be much smaller in size. Figure 9 also suggests that for *C* < 90%, bubble chord lengths will most likely be no more than 20 mm and that bubbles close to the pressure outlet will have small sizes.

A positive correlation is observed between the air concentration and air bubble chord length among the upper, lower and side nappes (Figure 12). The distribution can be defined according to the function:

$$(ch\_{ab})\_{\text{mean}} = \frac{1}{a\_1 + a\_2 \ast \mathbb{C}^{a\_3}}\tag{6}$$

where *a*1, *a*<sup>2</sup> and *a*<sup>3</sup> are empirical coefficients that may be related to the interaction. For the present experimental data, *a*<sup>1</sup> = 0.212, *a*<sup>2</sup> = −0.006 and *a*<sup>3</sup> = 0.767 were used. The coefficients of determination were set to be 0.85. 0.93 and 0.91 for the upper, lower and side nappes respectively.

Figure 13 presents a plot of experimental data, (*chab*)mean, with calculated results using Equation (6). Since most of the values of (*chab*)mean calculated from Equation (6) fall around the line of perfect agreement in an area within ±20% error lines, although some discrepancy exists. The results exhibit agreement with the theoretical distribution curve of Equation (6) for all flow conditions at various positions.

As the relationships between the air concentration and air bubble frequency, as well as the air bubble chord length have been discussed above, the relationship between the air bubble count rate and bubble size is also considered (Figure 14). It is evident that the correlations are clearly non-symmetric, unimodal and positively skewed and may therefore be represented by a probability density function defined by the Gamma distribution. According to mathematical statistics, the gamma distribution is a continuous probability function that can be written as follows:

$$G a(\mathbf{x}) = \frac{1}{\Gamma(\alpha)} \beta^{-\alpha} \mathbf{x}^{\alpha - 1} e^{-\frac{\mathbf{x}}{\beta}}, \mathbf{x} > 0 \tag{7}$$

where *α* is the shape parameter, *β* is the scale parameter, and Γ(*x*) = <sup>∞</sup> <sup>0</sup> *<sup>t</sup>x*−1*e*−*<sup>t</sup> dt* is the gamma Function, and Γ(*x* + 1) = *x*Γ(*x*).

**Figure 11.** Air bubble chord length distributions. Upper nappe (**a**) (*h*, *b*) = (4.5 cm, 4.5 cm), *θ* = 0%. (**b**) (*h*, *b*) = (2.5 cm, 2.5 cm), *θ* = 10%; lower nappe (**c**) (*h*, *b*) = (4.5 cm, 4.5 cm), *θ* = 0%. (**d**) (*h*, *b*) = (2.5 cm, 2.5 cm), *θ* = 10%; side nappe (**e**) (*h*, *b*) = (4.5 cm, 4.5 cm), *θ* = 0%. (**f**) (*h*, *b*) = (2.5 cm, 2.5 cm), *θ* = 10%.

**Figure 12.** Measured air bubble chord length as a function of the air concentration. The solid lines are determined from Equation (6). Upper nappe (**a**) (*h*, *b*) = (4.5 cm, 4.5 cm), *θ* = 0%. (**b**) (*h*, *b*) = (2.5 cm, 2.5 cm), *θ* = 10%; lower nappe (**c**) (*h*, *b*) = (4.5 cm, 4.5 cm), *θ* = 0%. (**d**) (*h*, *b*) = (2.5 cm, 2.5 cm), *θ* = 10%; side nappe (**e**) (*h*, *b*) = (4.5 cm, 4.5 cm), *θ* = 0%. (**f**) (*h*, *b*) = (2.5 cm, 2.5 cm), *θ* = 10%.

**Figure 13.** Comparison of experimental data with calculated values of air bubble chord length. (**a**) upper nappe; (**b**) lower nappe; (**c**) side nappe.

**Figure 14.** *Cont*.

**Figure 14.** Dimensionless bubble frequency as a function of the bubble chord length compared with the theoretical curve derived from Equation (9). Upper nappe (**a**) (*h*, *b*) = (4.5 cm, 4.5 cm), *θ* = 10%. (**b**) (*h*, *b*) = (6.5 cm, 6.5 cm), *θ* = 10%; lower nappe (**c**) (*h*, *b*) = (4.5 cm, 4.5 cm), *θ* = 10%. (**d**) (*h*, *b*) = (6.5 cm, 6.5 cm), *θ* = 10%; side nappe (**e**) (*h*, *b*) = (4.5 cm, 4.5 cm), *θ* = 10%. (**f**) (*h*, *b*) = (6.5 cm, 6.5 cm), *θ* = 10%.

Equation (7) is a classic two-parameter gamma distribution. Considering the complexity of the air bubble frequency and chord length in aeration flows, the five-parameters form of the generalized gamma distribution (also known as the generalized gamma function with four parameters after translation) is used in this study. The probability density function is written as follows:

$$Ga(\mathbf{x}) = \frac{1}{\Gamma(a\_1)} a\_2^{\eta\_1 a\_1} (\mathbf{x} - a\_3)^{\eta\_2 a\_1} e^{-a\_4 \ast (\mathbf{x} - a\_3)^{a\_5}} \qquad \mathbf{x} > \mathbf{0} \tag{8}$$

where *η*<sup>1</sup> and *η*<sup>2</sup> are the coefficients; *a*1, *a*2, *a*3, *a*4, and *a*<sup>5</sup> are the five parameters, of which *a*<sup>1</sup> and *a*<sup>5</sup> denote the shape parameters, *a*<sup>2</sup> and *a*<sup>4</sup> denote the scale parameters, and *a*<sup>3</sup> is the threshold parameter; and Γ(·) is the gamma function.

Based on the experimental data and the theoretical distribution curve of Equation (8), the relationship between the air bubble frequency and air bubble chord length can be fitted as follows:

$$\frac{F\_a}{(F\_a)\_{\text{max}}} = 115 [(ch\_{ab})\_{\text{mean}} - 1.5]^{3.8} \times e^{-7.65 [(ch\_{ab})\_{\text{mean}} - 1.5]^{0.25}} \tag{9}$$

Examples of a "modified" five-parameters gamma distribution that were used to fit the air bubble chord length distributions (Figure 14). The data is reasonably well represented by a "modified" gamma function distribution, although there is significant scatter.

Additionally, Figure 15 presents a plot of experimental data, *Fa*/(*Fa*)max, with calculated non-dimensional air bubble frequency using Equation (9). Since most of the values of *Fa*/(*Fa*)max calculated from Equation (9) fall around the line of perfect agreement in an area within ±20% error lines, suggesting that Equation (9) can be considered adequate in calculating the air bubbles chord length distribution. The "modified" gamma function is skewed to the left with a peak that lies in the small chord size values. The mode was found to be between 0.5 and 30 mm.

**Figure 15.** Comparison of experimental data with calculated values of dimensionless bubble frequency. (**a**) upper nappe; (**b**) lower nappe; (**c**) side nappe.

### **4. Discussion**

Air concentration is a macro parameter that is used to describe air–water flows. However, it is still insufficient in providing insight that would lead to a deeper understanding of the said phenomena. The process of aeration defines the general developing process of air moving into a body of water, i.e., the air bubbles. As such, the micro characteristics of air bubbles such as its count rate and size are crucial parameters in the study of aeration flows.

A double-tip conductivity probe used in present study to detect the air–water interfaces at different fixed points downstream of a three-dimensional (vertical drop and lateral enlargement) aerator. Some simple expressions were developed to simulate the distributions of air bubble frequency and bubble chord sizes. The relationships between the air concentration, the air bubble frequency and the air bubble chord length were also established. The obtained air concentration distributions satisfy the analytical model of Chanson [12], Equation (4) for air bubble frequency distributions as observed in supercritical open channel flows and backward-facing step flows (Toombes and Chanson [19]). A "modified" gamma probability density function was found to provide a good fit to the air bubble chord length distributions measured from three-dimensional aerator flows. This is different from the log-normal probability density functions that Chanson (1997) used for supercritical open channel flows. The distributions of the air bubble chord length showed a positive correlation with the void fraction.

Physical modelling may provide some information on the flow motion if a suitable dynamic similarity is selected. Air–water flows depend on Froude, Reynolds, and Weber numbers. Strictly speaking, dynamic similarity of air concentration and air bubble dissipation for free-surface aeration on aerator flows on a reduced-scale model is not possible (Chanson [32]) because it is impossible to satisfy simultaneously Froude and Reynolds similarities. Hence, the experimental results cannot be directly extrapolated to prototypes unless working at full scale. Nevertheless, with the assumption of Froude similarity, the Reynolds number (Re) should be larger than 1 × 105 to minimize the scale effects (Chanson [32]; Heller [33]). This limit is considered in the present study (Re ~0.89–1.34 × 106). Actually, results from scale experiments using Froude similarity are on the safe side for engineers with respect to air concentrations due to aeration tends to be underestimated in the smaller experiments.

Due to the complexity of bubble distribution in the air–water interfaces, the present results only serve as a preliminary investigation of the air bubble properties in the air–water mixtures downstream of a three-dimensional aerator. Hence, more detailed research is needed to study the air bubbles in aeration flows.

### **5. Conclusions**

A number of experiments were performed under various configurations and approach flow velocities for a 3D aerator. The air concentration, air bubble frequency and bubble chord length were

recorded in the developing shear layer of the upper, lower and side nappes downstream. The results are summarized as follows:


**Author Contributions:** S.L. and J.Z. proposed the research ideas and methods of the manuscript and writing, X.C. put forward the revise suggestion to the paper. J.C. is responsible for data collection and creating the figures.

**Funding:** This research was funded by the Chinese Scholarship Council (CSC), grant number 201804910306; the CAS "Light of West China" Program, grant number Y6R2220220; the National Science Fund for Distinguished Young Scholars, grant number 51625901; and the Natural Science Foundation of China, grant number 51579165.

**Conflicts of Interest:** The authors declare no conflicts of interest.

### **Notations**

The following symbols are used in this paper:



### **References**


© 2018 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

### *Article* **Experimental Study on the Air Concentration Distribution of Aerated Jet Flows in a Plunge Pool**

### **Weilin Xu \*, Chunqi Chen and Wangru Wei**

State Key Laboratory of Hydraulics and Mountain River Engineering, Sichuan University, Chengdu 610065, China; ccq-cedric@163.com (C.C.); wangru\_wei@hotmail.com (W.W.)

**\*** Correspondence: xuwl@scu.edu.cn

Received: 20 October 2018; Accepted: 30 November 2018; Published: 4 December 2018

**Abstract:** There is a lack of knowledge on the air concentration distribution in plunge pools affected by aerated jets. A set of physical experiments was performed on vertical submerged aerated jet flows impinging a plunge pool. The air concentration distribution in the plunge pool was analyzed under different inflow air concentrations, flow velocities, and discharge rate conditions. The experimental results show that the air concentration distribution follows a power-law along the jet axis, and it is independent of the initial flow conditions. A new hypothetical analysis model was proposed for air diffusion in the plunge pool, that is, the air concentration distribution in the plunge pool is superposed by the lateral diffusion of three stages of the aerated jet motion. A set of formulas was proposed to predict the air concentration distribution in the plunge pool, the results of which showed good agreement with the experimental data.

**Keywords:** submerged jets; aerated flow; air concentration; plunge pool

### **1. Introduction**

The impingement of a jet flow on a floor is a common issue that is important for many engineering applications, such as the flood discharging of hydraulic structures, and heat and mass transfer in industrial operations. Due to the presence of the floor being located at a distance from the nozzle exit, the diffusion of jets in a plunge pool is different from that of submerged free jets. Based on how the jet impacts the floor, the flow field can be divided into the free jet region, the impingement region, and the free shear region. The jet flow velocity conforms to the distribution law of submerged free jets in the free jet region. In the impingement region, the jet flow velocity decreases rapidly and reaches zero at the stagnation point [1,2].

Aerated and non-aerated impinging jet flows have been studied by many associated investigators [3–5]. Ervine summarized the characteristics of free turbulent jets, and compared the jet diffusion in the plunge pool with submerged and impinging jets [6]. Ervine collected data for both mean and fluctuating components of the pressure field on a plunge pool floor subjected to jet impingement, and compared the data for circular jets with those for wide rectangular nappes and rectangular slot jets [7]. Castillo analyzed the degree of break-up of a rectangular jet before entering the basin, and its relationship to the pressure fluctuation on the plunge pool floor [8]. Chanson discussed the characteristics and similitude of the air–water flow caused by impinging jets [9]. Wei analyzed the influence of aeration and initial water thickness on the axial velocity attenuation of jet flows using experiments and numerical calculations [10,11]. The research concludes that increasing the air concentration of jet flows and decreasing the initial jet thickness are effective ways to improve the axial velocity attenuation of jet flows. Some researchers focused on the pressure on the floor of the plunge pool by an aerated and non-aerated jet flow [12–15]. The results showed that the decreasing effect of air concentration on the time-averaged impact pressure is minimal, if any. However, the fluctuating pressure rises as the air concentration increases. Duarte compared the jet impingement on the center and side of a block embedded on the floor using plunging and submerged jets [16,17]. Others have inquired into the scouring of the plunge pool floor due to jets [18–20]. In addition to the velocity and pressure fields in the plunge pool, the process of air bubble entrainment at the impinging point of the water surface has also been studied extensively [21,22]. From the interaction of the jet flow and the plunge pool surface, an air sheet forms around the impinging point before air bubbles are entrained into the water surface [23]. Several correlative and predictive models on the air entrainment ratio have been reported, due to plunging jets under different working conditions [24–26]. Some researchers used numerical simulations to study aerated jet flows. Brouillio performed a series of two-dimensional numerical experiments on a translating impacting jet, and studied the dynamics of air entrainment by the jets impacting on the surface of the water [27]. Samanta studied the effect of the permeable wall on secondary cross-stream flow by performing direct numerical simulations of the fully developed turbulent flow through a porous square duct [28]. In model tests, the jet velocity profiles and pressure fluctuations on the plunge floor have been studied extensively, with detailed results available. Numerical simulations can essentially reflect characteristics such as the flow field and air entrainment process on the water surface. However, for mass and heat transfer characteristics at the air–water interface under complex conditions, a large number of systematic tests are still needed.

Few researchers have studied the air concentration distribution in the entire plunge pool, including the free shear region. Dong described the concentration diffusion of an aerated jet in the plunge pool with different plunging angles [29]. The air concentration along the jet axis gradually decreased, following a hyperbolic curve. Chanson conducted separate experiments with a vertical supported jet and a horizontal hydraulic jump [30]. The results showed that the air concentration in both cases exhibits a Gaussian distribution, and the longitudinal maximum air content declines exponentially. Khaled used the Particle Image Velocimetry (PIV) technique to characterize the flow field of both aerated and non-aerated jet flows beneath the water surface [31]. However, analyses and predictions of air concentration distribution in the entire plunge pool, especially where affected by the swirling jet, are still missing.

In the present paper, measurements and analyses of the air concentration distribution are presented for a plunge pool. A new analysis model is proposed to describe the air diffusion, with jet swirling in the entire plunge pool. A set of formulas has been built up to predict the air concentration distribution, and the results are believed to be useful for understanding and predicting the air concentration distribution in both the plunge pool and downstream.

### **2. Experimental Methods**

Experiments were carried out at the State Key Laboratory of Hydraulics and Mountain River Engineering (Sichuan University). An experimental facility was set up, as shown in Figure 1. The jet was aerated before being plunged vertically into the still water in a rectangular glass flume of 4.00 m × 0.25 m × 0.80 m. The airflow was provided using an air compressor with a power rating of 7.5 kW at a pressure of 0.8 MPa. The inflow jet was controlled using valves and measured with a TDS-100H ultrasonic flowmeter, of which the measurement range was 0.8–128.0 L/s, and the relative error of measurement was ±1%. The air was controlled using valves, and measured with a SLDLB-Y150D vortex flowmeter, of which the measurement range was 2.0–60.0 L/s, and the relative error of measurement was ±1%. To form a stable and well-mixed air-water flow, the air was forced through micro-perforated plates with a diameter of 1 mm, and mixed into the water flow in the air entrainment region. The water in the flume flowed out at both ends, where the water level was maintained at *H* = 0.55 m. The incident angle of the jet flow was 90 degrees. The rectangular nozzle, 0.20 m in width and 0.02 m–0.05 m at the opening, was installed slightly below the water surface to ensure a precise, stable initial air concentration.

**Figure 1.** Experimental facility.

The coordinate origin of an *x*-0-*y* coordinate system was located on the flume floor, the horizontal *x*-axis faced the downstream direction of the plunge pool, and the vertical *y*-axis coincided with the jet axis. In total, there were 264 measurable points for the air concentration in the plunge. A total of 22 rows were arranged along the x-axis direction. First, six rows were set at an interval of 5 cm (0 cm ≤ *x* ≤ 25 cm), and the following 16 rows were set at an interval of 10 cm (30 cm ≤ *x* ≤ 180 cm). In total, 12 rows of measuring points were arranged along the *y*-axis direction, with the first row at *y* = 1 cm and the last row at *y* = 54 cm (1 cm from the flume floor). The remaining 10 rows of measuring points along the *y*-axis were set between *y* = 5 cm and *y* = 50 cm at intervals of 5 cm. The air concentration was tested using a CQ6-2004 resistance-type air concentration meter (China Institute of Water Resources and Hydropower Research) with a resolution of 1%. The measurement principle of the resistance-type air concentration meter was to determine the air concentration by detecting the clear water resistance and the aerated water resistance between the two electrodes, and averaging the air concentration values in the selected integration time. Accurate flow measurements can be an important factor in obtaining adequate conclusions [32]. Different integration times for the air concentration measurements were tested beforehand, as shown in Figure 2. It can be seen that the mean air concentration data fall in the measurement accuracy zone if the integration time of the air concentration meter exceeds 60 seconds. Therefore, the sampling period for each measurement point was set to 60 seconds to minimize random error.

**Figure 2.** The uncertainty of air concentration measurements in the plunge pool.

The experiments were performed under different conditions, including the jet flow thickness *d*<sup>0</sup> (2 cm, 3 cm, 4 cm, 5 cm), the initial air concentration *C*<sup>0</sup> (10%, 15%, 30%, 40%, 50%, 60%), the clear water flow rate *Q*w (14.0 L/s–35.2 L/s), and the air flow rate *Q*a (2.1 L/s–42.0 L/s), as listed in Table 1.


**Table 1.** Experimental conditions.

In the experiment, the aerated jet was vertically injected and outflowed simultaneously at both ends of the flume, which was slightly wider than the nozzle exit. Jets could develop on the bottom plate of the flume under all experimental conditions. Thus, the development of the aerated jet in the plunge pool could be simplified as a two-dimensional flow, which is symmetrical to the jet axis (*y*-axis). For typical working conditions (*d*<sup>0</sup> = 2 cm, *C*<sup>0</sup> = 30%, *Q*<sup>w</sup> = 16.8 L/s), an experimental photograph and a contour map of the air concentration distribution on one side of the jet axis are shown in Figure 3. The majority of the air bubbles flow with the jet flow in the free jet region, before being blocked by the flume floor and deflected by 90 degrees into the wall jet region on both sides. Afterwards, the air bubbles flow gradually towards the water surface, with some eventually escaping into the surrounding air, whereas the rest are carried into the swirling region.

(**a**) Photograph (**b**) Air concentration contour map

**Figure 3.** Air concentration distribution for *d*<sup>0</sup> = 2 cm, *C*<sup>0</sup> = 30%, *Q*<sup>w</sup> = 16.8 L/s.

### **3. Results and Discussion**

### *3.1. Air Concentration Distribution in the Free Jet Region*

Along the jet axis, the air concentration distribution should be closely related to the jet velocity distribution. The jet velocity decreases sharply as the aerated jet approaches the plunge pool floor. Due to the pressure gradient, the jet carrying air bubbles drifts toward both sides of the plunge pool, resulting in a steep decrease in air concentration in the impingement region. The region higher than one-fifth of the water depth is considered to be the free jet region, where the attenuation of air concentration was studied in the test.

Considering the influencing factors, the initial water velocity *V*w, the initial air concentration *C*0, and the clear water flow rate *Q*<sup>w</sup> were adopted, and a series of experiments were performed to test each effect on the air concentration attenuation, as shown in Figure 4. Under a different initial water velocity (3 m/s < *V*w < 8 m/s) and a clear water flow rate (21.1 L/s < *Q*w < 35.2 L/s), the attenuation of *C*m/*C*<sup>0</sup> varies in a similar pattern along the jet axis. The reason for the results is that the velocity attenuation along the jet axis changes little with different flow conditions, including the initial jet velocity and the clear water flow rate [33].

**Figure 4.** Influence of *V*w and *Q*w on axial air concentration distribution.

By theoretically analyzing the concentration flux of the aerated jet based on the self-similarity of the velocity and air concentration distribution, the power-law relation was found between the axial air concentration and the flow distance in the free jet region [29]. The relationship between the axial air concentration *C*<sup>m</sup> and the distance from the nozzle *x'* = *H* − *y*, normalized by the initial air concentration *C*<sup>0</sup> and the jet flow thickness *d*<sup>0</sup> in the plunge pool, respectively, is written as:

$$\frac{\mathcal{C}\_{\rm m}}{\mathcal{C}\_{0}} = k\_{0} \times \left(\frac{\mathbf{x}'}{d\_{0}}\right)^{-0.5} = k\_{0} \times \sqrt{\frac{d\_{0}}{H - y}}\tag{1}$$

where *k*<sup>0</sup> is a coefficient, suggested as approximate 1.2 in this study.

The calculated results of the axial air concentration distribution in the free jet region using Equation (1) was compared with the experimental data, as shown in Figure 5a. Independent of the initial jet velocity and the clear water flow rate, the axial air concentration sharply decreased to about 50% of *<sup>C</sup>*<sup>0</sup> in a short distance ( *<sup>H</sup>*−*<sup>y</sup> <sup>d</sup>*<sup>0</sup> ≈ 5) from the nozzle exit. Due to the influence of the pressure gradient near the plunge floor, further from the nozzle exit, the air concentration distribution was more dispersed in the free jet region. It can be derived from the experimental results that the air concentration distribution along lateral cross-sections in the free jet region conformed to the Gaussian distribution law, as shown in Figure 5b. In the free jet region, no significant influence of different initial flow conditions was observed on the cross-sectional air concentration distribution, which tended to be self-similar. Downstream of the free jet region, the cross-sectional air concentration distribution was affected by the swirling jet in the plunge pool.

**Figure 5.** Air concentration distribution in the free jet region.

For a two-dimensional constant jet along the *x*' direction in an *x*'–0–*y*' coordinate system, the continuity equation of the air concentration diffusion is:

$$\frac{\partial(\text{CLI})}{\partial \mathbf{x'}} + \frac{\partial(\text{CV})}{\partial \mathbf{y'}} = \frac{\partial}{\partial \mathbf{y'}} \left( D\_y \frac{\partial \mathbf{C}}{\partial \mathbf{y'}} \right) \tag{2}$$

where *C* is the air concentration, *U* is the velocity in the *x*' direction, *V* is the velocity in the *y*' direction, and *Dy* is the turbulent diffusion coefficient.

The longitudinal velocity *U* can be approximated using the average velocity *u*. The transverse velocity can be negligible, i.e., *V* ≈ 0.

According to Prandtl's hypothesis of free turbulence theory [34], *Dy* can be approximated as:

$$D\_y = k\overline{\pi}b$$

where *k* is the coefficient, and *b* is the thickness of the jet, which is assumed to expand linearly and can be written as:

$$b = mx' = \tan \theta \times \mathbf{x'} \tag{4}$$

where *θ* is the diffusion angle.

Previous research suggests that with respect to the conventional diffusion angle (about eight degrees), the diffusion angle of aerated water increases due to the influence of bubbles [10]. Twelve degrees is taken as the diffusion angel in this study.

Due to hardly any influence of initial flow conditions on both the axial attenuation and the lateral distribution of air concentration in the free jet region, it is assumed that the air concentration distribution in the free jet region is self-similar, resembling that in the velocity distribution, which can be written as:

$$\frac{\mathbb{C}}{\mathbb{C}\_{\text{m}}} = f\_1(\eta), \quad \eta = \frac{y'}{b} \tag{5}$$

Combining Equations (2), (3) and (5), *<sup>C</sup> <sup>C</sup>*<sup>m</sup> can be expressed as:

$$f\_1 = \frac{\mathcal{C}}{\mathcal{C}\_{\rm m}} = \mathbf{A} \times \exp\left[-\frac{m}{2k} \times \eta^2\right] \tag{6}$$

where A is the integral constant and determined by boundary conditions. When *η* = *<sup>y</sup> <sup>b</sup>* = 0, *C* = *C*m, A = 1, Equation (6) can be rewritten as:

$$\frac{\mathcal{C}}{\mathcal{C}\_m} = \exp\left[-\frac{m}{2k} \times \left(\frac{y'}{b}\right)^2\right] \tag{7}$$

The concentration half-width is defined as:

$$\left. b\_{1/2} = y' \right|\_{c=0.5C\_m} \tag{8}$$

It is assumed that the concentration half-width extends linearly along the path, thus:

$$b\_{1/2} = \frac{1}{2}d\_0 + \tan\theta \times x'\tag{9}$$

Combining Equations (7) and (8):

$$
\ln\left(\frac{1}{2}\right) = -\frac{m}{2k} \times \left(\frac{b\_{\frac{1}{2}}}{b}\right)^2 = -0.693\tag{10}
$$

Considering:

$$-\frac{m}{2k} \times \left(\frac{y'}{b}\right)^2 = -\frac{m}{2k} \times \left(\frac{b\_{\frac{1}{2}}}{b}\right)^2 \times \left(\frac{y'}{b\_{\frac{1}{2}}}\right)^2 = -0.693 \times \left(\frac{y'}{b\_{1/2}}\right)^2\tag{11}$$

Combining Equations (7) and (11), the air concentration distribution on the flow cross-section can be calculated by:

$$\frac{\text{C}}{\text{C}\_{\text{m}}} = \exp\left[-0.693 \times \left(\frac{y'}{b\_{1/2}}\right)^2\right] \tag{12}$$

where *C* is the transverse air concentration diffusion, *C*m is the corresponding air concentration along the jet axis, and *y*' is the transverse distance from the jet axis.

The calculated result of Equation (12) was compared with the experimental results, as shown in Figure 5b. In the free jet region, the calculated results of the lateral diffusion were essentially consistent with the experimental data. Further away from the axis, due to the influence of the swirling jet, the experimental value of the air concentration was greater than the calculated result obtained using Equation (13). Therefore, the next section will focus on the establishment of an analysis model for air concentration diffusion in the plunge pool.

In the *x*–0–*y* coordinate system in this study, the air concentration distribution in the free jet region can be calculated by:

$$\begin{aligned} \mathcal{C} &= \mathcal{C}\_0 \times 1.2 \sqrt{\frac{d\_0}{x}} \times \exp[-0.693 \left(\frac{y'}{\frac{b\_1}{2}}\right)^2] \\ &= \mathcal{C}\_0 \times 1.2 \sqrt{\frac{d\_0}{H - y}} \times \exp[-0.693 \left(\frac{x}{b\_{1/2}}\right)^2] \end{aligned} \tag{13}$$

### *3.2. Air Concentration Distribution in the Swirling Region of the Plunge Pool*

For one side of the plunge pool, the experimental photographs of the submerged jet motion and air concentration distribution are shown in Figure 6. After being plunged into the still water, the aerated jet moves along the vertical axis, and the aerated wall jet forms after a 90-degree deflection of the jet flow near the flume floor. While the wall jet moves downstream, it is gradually developed toward the free surface under the combined effects of turbulence and buoyancy. Due to the existence of the downstream board, part of the aerated flow moves upstream along the water surface and forms a swirling region in the plunge pool. At the downstream end of the swirling region, the direction of

the aerated flow movement is almost vertically upward. While moving with the water flow, the air bubbles expand continuously and laterally into the swirling region under the effects of turbulence.

(**a**) *d*0 = 2 cm, *C*0 = 15%, *Q*w *=* 16.1 L/s

(**b**) *d*0 = 3 cm, *C*0 = 30%, *Q*w = 21.1 L/s

(**c**) *d*0 = 2 cm, *C*0 = 60%, *Q*w = 14.5 L/s

**Figure 6.** Air concentration distribution in the plunge pool.

According to Melo [13], the flow patterns in the plunge pool impinged by vertical jets are composed of three distinct regions, as shown in Figure 7a. In the free jet region, the vertical jet develops by shearing with the surrounding water body without the influence of the plunge floor. The impact of the jet against the floor builds up the pressure gradient and deflects the jet parallel to the plunge floor into the wall jet region. Due to the limited area of the plunge pool, the horizontal wall jet develops upwards and flows into the swirling region.

**Figure 7.** Generalized calculation of air concentration in the swirling region.

When calculating the concentration field in the swirling region, the motion and diffusion of the air bubbles in the plunge pool are generalized into the following three stages, as shown in Figure 7b:

Stage 1: The aerated jet moves along the jet axis until it reaches the floor of the plunge pool, categorized as the free jet region in Figure 7a. The attenuation and lateral diffusion of the air concentration occurs along the path, which is consistent with the distribution law in the free jet region mentioned above. At this stage, the air concentration that diffuses laterally into the swirling region is taken as *C*1;

Stage 2: The aerated flow is deflected by 90 degrees at the bottom of the plunge pool, and moves along the axis near the bottom plate, categorized as the wall jet region in Figure 7a. Regardless of the local and instantaneous air concentration loss, the initial air concentration at this stage is assumed to be the same as the air concentration at the end of the axis in the previous stage. The air concentration of the aerated jet continues to decrease along the path, until the jet reaches the end section of the swirling region. At this stage, the air concentration that diffuses vertically into the swirling region is taken as *C*2;

Stage 3: The aerated flow moves upstream along the free surface until it reaches the original jet axis at stage 1, categorized as the swirling region in Figure 7a. The initial air concentration value in this stage is obtained by approximating the axis of the aerated jet in the swirling region to a rectangle with a length of *L*' and a width of *H*. The air concentration that diffuses vertically into the swirling region is taken as *C*3.

In the plunge pool, excluding the jet axis, the air concentration at a certain point can be considered as the superposition of *C*1, *C*2, and *C*3. When calculating the axial air concentration attenuation and the lateral diffusion in the three stages mentioned using Equation (13), separately, the jet flow distance *x*' and the transverse distance from the axis *y*' can be expressed as:

$$\text{stage 1}: \ x\_1' = H - y, \quad y\_1' = \ge \tag{14}$$

$$\text{stage 2: } \ x\_2' = H + \text{x, } \ y\_2' = y \tag{15}$$

$$\text{stage 3}: \ x\_3' = 2H \, + \ 2L' - x, \quad y\_3' = H - y \tag{16}$$

In order to locate the end section in stage 2 of the air concentration diffusion, the swirling region length *L*' is defined as the distance between the jet axis and the vertical cross-section, where the mean air concentration equals 5% of *C*0.

$$L' = \mathfrak{x}|\_{c\_{\text{mean}} = 0.05c\_0} \tag{17}$$

where *C*mean is the mean air concentration.

Under current experimental conditions, the relationship between the swirling region length and the initial water flow velocity is shown in Figure 8, which can be correlated with the linear function as:

$$\frac{L\prime}{d\_0} = 8 \times \frac{V\_{\rm w}^2}{gH} + 11, \quad 0.6 < \frac{V\_{\rm w}^2}{gH} < 6\tag{18}$$

where *V*<sup>w</sup> is the water velocity, and *g* is the acceleration of gravity. The correlation coefficient is 0.80. It should be noted that the length of the swirling region was deliberately defined using the statistical mean air concentration of the vertical cross-section. The swirling region was actually determined by the flow field. In the downstream end of the plunge pool, the decrease in flow velocity led to changes in the following behaviors of the air bubbles. This resulted in the deviation of measurements in the experiment.

**Figure 8.** Distribution of swirling region length *L*'.

In stage 2 of the air concentration diffusion, the position of the maximum air concentration values along the vertical cross-sections in the wall jet can be estimated as:

$$y\_{\rm ms} = \frac{1}{3} \times b$$

where *b* is the thickness of the jet.

According to Beltaos [35], the relation between the maximum velocity of the wall jet *u*m1 and the initial velocity *u*<sup>0</sup> is:

$$\frac{u\_{\rm m1}}{u\_0} \sqrt{\frac{H}{d\_0}} = 2.77 \left\{ 1 - \exp \left[ -38.5 \left( \frac{\chi}{H} \right)^2 \right] \right\}^{\frac{1}{2}} \tag{20}$$

Considering the symmetrical aerated flow on both sides of the plunge pool, and the relation:

$$\frac{1}{2}u\_{\text{m1}} \times b = \frac{1}{2}u\_0 \times d\_0 \tag{21}$$

Combining Equations (19)–(21), *y*ms can be expressed as:

$$y\_{\rm ms} = \frac{1}{3} \times \frac{\sqrt{H \times d\_0}}{2.77 \left\{ 1 - \exp \left[ -38.5 (\text{x}/H)^2 \right] \right\}^{1/2}} \tag{22}$$

The air concentration can be linearly solved on the vertical cross-sections between the axis of the wall jet and the plunge pool floor, where the air concentration value is zero.

*Water* **2018**, *10*, 1779

According to the above analysis, at any point (*x*, *y*) in the swirling region of the plunge pool, the calculation of the lateral air concentration diffusion *C*<sup>1</sup> from the jet axis can be written as:

$$\mathcal{C}\_1 = \mathcal{C}\_0 \times 1.2 \sqrt{\frac{d\_0}{H - y}} \times \exp[-0.693 \left(\frac{\text{x}}{b\_{1/2}}\right)^2] \tag{23}$$

The calculation of the vertical air concentration diffusion *C*<sup>2</sup> from near the plunge pool floor can be written as:

$$\mathcal{C}\_2 = \mathcal{C}\_0 \times 1.2 \sqrt{\frac{d\_0}{H + \text{x}}} \times \exp[-0.693 \left(\frac{y - y\_{\text{ms}}}{b\_{1/2}}\right)^2] \tag{24}$$

The calculation of the vertical air concentration diffusion *C*<sup>3</sup> from the free surface can be written as:

$$\mathcal{C}\_3 = \mathcal{C}\_0 \times 1.2 \sqrt{\frac{d\_0}{2H + 2L' - \mathbf{x}}} \times \exp\left[-0.693 \left(\frac{H - y}{b\_{1/2}}\right)^2\right] \tag{25}$$

Thus, the air concentration value at the point (x, y) can be obtained by:

$$\mathcal{C}\Big|\_{(x,y)} = \mathcal{C}\_1 + \mathcal{C}\_2 + \mathcal{C}\_3 \tag{26}$$

with 12 degrees of diffusion angles used in calculating *C*<sup>1</sup> and *C*2, and zero degrees in calculating *C*3.

Combing Equations (23)–(26), it is possible to calculate the air concentration at any point in the swirling region of the plunge pool.

### *3.3. Comparison of Calculated and Experimental Results*

The calculated and experimental results of the air concentration distribution along the vertical cross-sections in the plunge pool are compared in Figure 9. When the initial air concentration was small, such as *C*<sup>0</sup> < 40%, the calculated results agreed well with the experimental measurements. When the initial air concentration was large, such as *C*<sup>0</sup> > 40%, the calculated results could essentially reflect the variation of the air concentration distribution in the region close to the jet axis, such as *x*/*d*<sup>0</sup> < 15–20. In the downstream region, on the other hand, such as *x*/*d*<sup>0</sup> > 15–20, the calculated results were somewhat different from the experimental measurements.

The reason for the deviation is that the bubble motion and diffusion process in the hypothetical analysis model is different from that in the actual jet flow. The influence of the bubble floating process by buoyancy on the actual jet flow was not considered, leading to a decrease in the adaptability of the idealized hypothetical analysis in this paper. However, the absolute value of the air concentration in the downstream region of the plunge pool was relatively small. For example, for *d*<sup>0</sup> = 2 cm, *C*<sup>0</sup> = 60%, *Q*<sup>w</sup> = 14.5 L/s, the absolute value of the air concentration was essentially less than 0.05 in the downstream region, such as *x*/*d*<sup>0</sup> > 10. Consequently, the present analysis model can essentially predict the air concentration distribution of aerated jet flows in the swirling region of the plunge pool.

**Figure 9.** Comparison of calculated and experimental results.

### **4. Conclusions**

In this paper, a series of experimental studies on the air concentration distribution of vertical submerged aerated jet flows in a plunge pool were carried out. Based on the distribution law of the air concentration attenuation along the axial and lateral directions, an analysis model of the air concentration distribution in the plunge pool was proposed. The conclusions were as follows.

In the free jet region, the axial air concentration attenuation followed a power-law distribution, and the air concentration distribution along the lateral cross-sections conformed to the Gaussian distribution law. Both the axial attenuation and the lateral distribution of air concentration in the free jet region were almost unchanged for different initial water velocities, air concentrations, and flow rate conditions. Based on the observation and analysis of the flow patterns in the swirling region in the plunge pool, the air concentration was assumed to be a superposition of lateral diffusion in three stages of the jet flow. The position of the maximum air concentration values along the vertical cross-sections in the wall jet and the length of the swirling region, as well as their relationship with the initial flow conditions, were introduced. A set of formulas was proposed for the prediction of the air concentration distribution in the plunge pool. The present analysis model, the results of which have been compared with experimental data, can essentially predict the air concentration distribution of aerated jet flows in the swirling region of a plunge pool.

In water conservancy projects, the jet flows are aerated, and they form an air concentration distribution inside the plunge pool, which has a significant effect on the downstream energy dissipation and the total dissolved gas. Meanwhile, the increase in aeration results in changes to the mass and heat transfer process at the air–water interface. This research improves our understanding of the air dissolution process downstream practical applications, and provides a new model to predict the air concentration distribution in plunge pools caused by aerated jets. The results of this study reveal the complexity of the air concentration distribution in a plunge pool caused by aerated jets. In order to further improve the prediction accuracy for the air concentration distribution in plunge pools, coupling analyses of bubble floating by buoyancy and jet motion in the air diffusion process should be the focus of future studies.

**Author Contributions:** The paper is the product of the joint efforts of the authors who worked together on experimental model tests. W.X. conceived the experiments, having a scientific background in applied hydraulics. C.C. and W.W. conducted the experimental investigations and analyzed the data under the supervision of W.X. C.C. wrote the manuscript with the co-authors.

**Funding:** Resources to cover the Article Processing Charge were provided by the National Natural Science Foundation of China (Grant No. 51609162).

**Conflicts of Interest:** The authors declare no conflict of interest.

### **References**


© 2018 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

### *Article* **Spatial Distribution Characteristics of Rainfall for Two-Jet Collisions in Air**

### **Hao Yuan, Weilin Xu \*, Rui Li, Yanzhang Feng and Yafeng Hao**

State Key Laboratory of Hydraulics and Mountain River Engineering, Sichuan University, Chengdu 610065, China; 18996152721@163.com (H.Y.); li.rui.sai@gmail.com (R.L.); 13835141175@163.com (Y.F.); 13258387216@163.com (Y.H.)

**\*** Correspondence: xuwl@scu.edu.cn

Received: 12 September 2018; Accepted: 4 November 2018; Published: 7 November 2018

**Abstract:** Many researchers have studied the energy dissipation characteristics of two-jet collisions in air, but few have studied the related spatial rainfall distribution characteristics. In this paper, in combination with a model experiment and theoretical study, the spatial distributions of rainfall intensity of two-jet collisions, with different collision angles and flow ratios, are systematically studied. The experimental results indicated that a larger collision angle corresponds to a larger rainfall intensity distribution. The dimensionless maximum rainfall intensity sharply decreased with the flow ratio, while the maximum rainfall intensity slightly increased when the flow ratio was greater than 1.0. A theoretical equation to compute the location of maximum rainfall intensity is presented. The range of rainfall intensity distribution sharply increased with the flow ratio. When the flow ratio was greater than 1.0, the range of longitudinal distribution slightly increased, whereas the lateral distribution remained unchanged or slowly decreased. A formula to calculate the boundary lines of the *x*-axis is proposed.

**Keywords:** rainfall intensity distribution; two water jets; collision in air; Gaussian distribution; trajectory line

### **1. Introduction**

At present, the high dam projects that have been built, are under construction, or are planned in China have the common features of large drops, a narrow river valley, and a large flood discharge flow. The flood discharge energy dissipation mode of the dam body has been mostly adopted by two-jet collisions and downstream water reservoirs, such as Ertan (Figure 1), Xiaowan, Xiluodu, Goupitan, and the Jinping grade I project. Two-jet collisions in air change the motion track of each of the two jets, significantly disperse the water flow, effectively alleviate the dynamic pressure impact of the water flow on the bottom plate of the water cushion reservoir, and have a remarkable energy dissipation effect [1–7]. However, compared with the traditional energy dissipation method, two-jet collisions in air appear to be a serious problem with regard to flood discharge atomization in the existing engineering examples, which has caused great damage to the normal operation, traffic safety, surrounding environment, and even the stability of the downstream bank slopes. It is necessary to study the spatial distribution characteristics of rainfall intensity in two-jet collisions in air, in order to improve the prediction accuracy of the discharge atomization rain strength and influence range, as well as make a relevant protection classification and protection design scheme in the design planning [8–10].

**Figure 1.** Simultaneous discharge of the surface and the deep hole of the Ertan hydropower station in 2010.

Many former researchers who have studied two-jet collisions in air have focused on the energy dissipation effect of collisions. Xiong [11] introduced the main factors that affect energy dissipation during an in-air collision between the surface and the deep hole. Guo [12] theoretically analyzed the effect of the collision angle on energy dissipation. Using the momentum integral equation of fluid mechanics, Liu [13] derived the related formula of energy dissipation for the collision of two water jets in air in detail, and introduced the concept and calculation method for the energy dissipation rate. Diao [14] studied the relationship between the collision angle, the flow ratio, and the collision energy dissipation effect through model testing, and proposed that the main function of jet collisions is to disperse the water flow. Sun [15] analyzed hydraulic characteristics, such as the energy loss of upand down-collisions, three-dimensional (3-D) diffusion and leakage collision of the surface, and deep orifice slugs. The 3-D collision velocity, collision efficiency, and collision energy loss were derived using the momentum equation and variation law of the air diffusion width of the water tongue, and the optimal hydraulic conditions for collision were obtained. Sun [16] applied the turbulence jet theory and flow momentum equation, in order to introduce the calculation formula of the collision velocity vector and collision energy dissipation efficiency of the combined water tongue when it collides with the left and right sides of the air.

At present, the research on flood discharge atomization is mainly about ski-jump atomization. Liang [17] proposed a calculation model of atomized water flow, and obtained various calculation formulas and methods in the field of atomized water flow influence. Liu [18,19] studies the diffusion law, velocity distribution, and energy loss of a water jet. According to the prototype observation data of the flood discharge atomization and some experimental data of the model, based on the comprehensive analysis of various factors that affect the atomization range, Li [20] discusses the rough estimation method of the range of the energy dissipation atomization precipitation area. Liu [21–24] introduces several key problems with regard to the shape and the numerical simulation of atomized flows, such as the jet calculation of the water tongue, collision between the water tongue and the water surface, and calculation of the fog source quantity. By collecting, inducing, and summarizing the partial flood discharge atomization of engineering prototype observation data, Sun [25] found the longitudinal boundary flood discharge atomization average discharge flow, and a good relationship between the water flow velocity and the water entry angle of the water tongue. Based on the dimensional analysis, a method was established to estimate the rainfall intensity, flood discharge atomization, and longitudinal boundary experience relationship. Liu [26] proposed an atomization forecast model based on artificial neural networks. The Monte Carlo method was applied to add the effect of

environmental, wind, and topographic factors to the atomization mathematical model of overpass discharge by Lian [27]. Sun [28] revised the resistance coefficient of particles, using the research results of raindrop movement, and obtained the resistance coefficient of the splashing water droplet movement. Then the effect of the diameter of a splashing water droplet on the splashing length is preliminarily discussed, which deepens the understanding of the splashing movement. Using the generalized model test in 2013 by Wang [29]—under different hydraulic conditions, to select a tongue that fell into the water downstream of an atomized water source area of rainfall intensity—the flood discharge atomization fog source area of the plane distribution features was measured and analyzed, to determine the reasons for the formation of different areas around the atomization source and the plane distribution of rainfall intensity. Lian [30] proposes a method to predict the spray atomization of spillage by combining a physical model and theoretical analysis. Zhang [31] experimentally studied the motion track, distribution range, water point shape, and velocity of single-strand jets under different air dosages. Although their methods are accurate in predicting the atomization of a ski-jump flow, they are not good at predicting the atomization of two-jet collisions, which add a collision zone and change the trajectory of two jets. In rocket propulsion, extensive experiments—mostly based on physical model investigations—have been conducted on two-jet collisions, including the liquid membrane [32–35], droplet size, and velocity distribution [36–41].

Most domestic and foreign scholars have focused on the atomization of flood discharge by single-strand pick flows. However, compared with the single-strand pick flow, two-jet collisions in air increase the area of collision atomization, where both range and intensity of the atomization increase. For example, Ertan applied this method to flood discharge and cause landslide by rainfall. Therefore, in this paper, through theoretical analysis and a physical model experiment, we studied the distribution characteristics of rainfall intensity after the collision of two jets, the trajectory lines of the water tongue, and the range of the *x*-axis after collision for different flow ratios and collision angles.

### **2. Model Design and Experiment**

The experiment was conducted in the State Key Laboratory of Hydrodynamics and Mountain River Engineering in Sichuan University. Figure 2 shows the schematic diagram of the model test. The water flux was supplied by a large rectangular tank (width of 4.0 m, length of 4.0 m, and height of 5.0 m). When laying the rainfall receiving platform, we considered the symmetry of the surface and deep hole, and only one side was arranged. Water was fed into the iron box through the pump behind the iron box, and the water level of the iron box was maintained through the valve in order to achieve a stable deep and surface flow.

**Figure 2.** *Cont*.

**Figure 2.** (**A**) Experiment layout: side view; (**B**) experiment layout: front view.

In order to be consistent with the prototype which adopts two-jet collisions in air, this experiment used a constant surface hole angle *θ*<sup>1</sup> (a pitch-down angle of −30◦) and adjusted the deep hole pick angle *θ*<sup>2</sup> (0◦~45◦) to change the collision angle *β*. The width of the surface deep hole was identical, the height of the deep hole could be adjusted, and the range was 0~5 cm. The maximum single-width flux of the model table hole was 0.1236 m2/s, the maximum single-width flux of the deep hole was 0.1431 m2/s, and the flux was measured by a triangular thin-walled weir (±1%). Four horizontal planes, under the collision points of 36 cm, 56 cm, 76 cm, and 96 cm, were set to measure the rainfall intensity. The platform was supported by a slide rail system, and the accuracy in the vertical platform position was less than 0.5 cm. Table 1 lists 18 working conditions in the experiment, including almost 30,000 measurement points. To obtain the four different sections, we first used image processing to find the collision point, and then adjusted the height of the measuring platform. Detailed visual observations of the collision points were conducted using a Canon EOS80D DSLR camera (Cannon, Chengdu, China).

Rainwater was collected in the horizontal plane (*xy* plane in Figure 3) at the different elevations below the collision point and 0.15 m above the faceplate. The rainfall collection platform included the movable base, support, and rainfall collection panel, which was made of an organic glass panel and a PVC pipe. In each plane, up to 121 PVC pipes (each with a diameter of 2.0 cm) were used simultaneously to collect rain within a period varying from 60 s to 1000 s, depending on the rain intensity in that region, which was found to be sufficient to ensure stable results. Zhang and Zhu [31] indicated that the measurement error could be ignored with different pipe diameters. In their test, the bottles (25 mL with a bottle-neck diameter of 1.52 cm) were tested to give the same (3% difference) result of rain intensity distribution as those of much larger plastic bottles (250 mL with a bottle-neck diameter of 3.26 cm). To avoid the impact of water splashing on the back of the faceplate, the faceplate was hollow, except for the holes and supporting frames. The lower end of the short pipe, whose range was 15~3000 mL, used a measuring cylinder to collect rainfall. Rainfall intensity varies from point to point, so different cylinders were used. To obtain the whole region's rainfall intensity, the platform was first placed in one region to collect rainwater, then the entire platform was moved 1.2 m in the longitudinal direction to the next region. To ensure measurement accuracy, two measurements were made in the same working conditions, and the average was taken.


**Table 1.** Working conditions in experiments.

*q*m is the flow discharge per unit width of deep hole. *q*c is the flow discharge per unit width of surface hole. *f* is the ratio of deep hole flow to surface hole flow, that is *f* = *q*m/*q*c. *Q* is the integral of measurement value.

**Figure 3.** Rainfall collection platform in 2017.

In each test, the total rain flux in the horizontal plane was integrated and compared to the total water flux of the surface and the deep holes. The total water/rain flux *Q* was calculated from *<sup>Q</sup>* <sup>=</sup> <sup>+</sup><sup>∞</sup> −∞ +∞ −∞ *I*d*x*d*y*, where *I* was the measured rain intensity. Comparisons of the integrated values with the water flux at the surface and the deep holes were listed in Table 1. For all tests, the conservation ratio of rain flux on average was 101.4 ± 7.0%. In runs T2-6 and T3-6, the loss of rain flux was slight larger (±9.3%), which may be related to the fact that the rain had a wider range in the test. Hence, the above comparison showed the reliability of the measurements.

### **3. Experimental Results and Discussion**

### *3.1. Spatial Distribution Characteristics of Rainfall Intensity*

Figures 4–6 show the distribution of rainfall intensity for several flow ratios, collision angles, and heights of the rainfall intensity collection platform from the collision point. The rainfall intensity boundary defined 5% of the maximum rainfall intensity under each working condition [31].

In Figure 4, the flow ratios vary, and the height and collision angles from the collision point to the measurement platform are constant. Note that on the section with the same vertical distance from the collision point, the rainfall intensity range showed an increasing trend on the *y*-axis when the flow ratio increased. On the *x*-axis, the rainfall intensity range increased first, and subsequently decreased slowly. For *f* = 0.3, the rainfall intensity range was concentrated in the range 0 cm ≤ *x* ≤ 22 cm; for *f* = 0.9, the rainfall intensity range was concentrated in the range of 10 cm ≤ *x* ≤ 100 cm; while when *f* = 3.2, the rainfall intensity range was concentrated in the range of 40 cm ≤ *x* ≤ 115 cm. Under all of the working conditions, the maximum rainfall intensity was on the *x*-axis.

**Figure 4.** *Cont*.

**Figure 4.** Rainfall intensity distribution at different flow ratios on the *x*-*o*-*y* measuring platform. (**a**) *θ*<sup>1</sup> = −30◦, *θ*<sup>2</sup> = 30◦, *β* = 77◦, *z* = 56 cm, and *f* = 0.3; (**b**) *θ*<sup>1</sup> = −30◦, *θ*<sup>2</sup> = 30◦, *β* = 77◦, *z* = 56 cm, and *f* = 0.6; (**c**) *θ*<sup>1</sup> = −30◦, *θ*<sup>2</sup> = 30◦, *β* = 77◦, *z* = 56 cm, and *f* = 0.9; (**d**) *θ*<sup>1</sup> = −30◦, *θ*<sup>2</sup> = 30◦, *β* = 77◦, *z* = 56 cm, and *f* = 1.3; (**e**) *θ*<sup>1</sup> = −30◦, *θ*<sup>2</sup> = 30◦, *β* = 77◦, *z* = 56 cm, and *f* = 2.5; (**f**) *θ*<sup>1</sup> = −30◦, *θ*<sup>2</sup> = 30◦, *β* = 77◦, *z* = 56 cm, and *f* = 3.2.

**Figure 5.** *Cont*.

**Figure 5.** Rainfall intensity distribution at different collision angles on the *x*-*o*-*y* measuring platform. (**a**) *θ*<sup>1</sup> = −30◦, *θ*<sup>2</sup> = 0◦, *β* = 48◦, *z* = 56 cm, and *f* = 1.3; (**b**) *θ*<sup>1</sup> = −30◦, *θ*<sup>2</sup> = 30◦, *β* = 77◦, *z* = 56 cm, and *f* = 1.3; (**c**) *θ*<sup>1</sup> = −30◦, *θ*<sup>2</sup> = 45◦, *β* = 90◦, *z* = 56 cm, and *f* = 1.3.

In Figure 5, the flow ratios and height from the collision point to the measurement platform are constant, and the collision angles vary. Note that with the increase in the collision angle, the rainfall intensity range on the *x*- and *y*-axes increased—for example, when *β* = 48◦, *L*<sup>x</sup> = 70 cm, and *L*<sup>y</sup> = 20 cm; when *β* = 77◦, *L*<sup>x</sup> = 108 cm, and *L*<sup>y</sup> = 68 m; when *β* = 90◦, *L*<sup>x</sup> = 163 cm, and *L*<sup>y</sup> = 80 cm. However, the maximum rainfall intensity decreased with the increased collision angle—for example, when *<sup>β</sup>* = 48◦, *<sup>I</sup>*max = 5.78 × <sup>10</sup><sup>5</sup> mm/h; and for *<sup>β</sup>* = 90◦, *<sup>I</sup>*max = 5.78 × <sup>10</sup><sup>5</sup> mm/h.

In Figure 6, the flow ratio and collision angle were constant, and the collision point varied from the height of the measurement platform. Note that (1) the range of rainfall intensity on the *x*- and *y*-axes gradually increases along the *z* direction, and (2) the maximum rainfall intensity value continues decreasing. At *<sup>z</sup>* = 56 cm, the maximum value is 5.74 × <sup>10</sup><sup>5</sup> mm/h; at *<sup>z</sup>* = 76 cm, the maximum value is 5.32 × <sup>10</sup><sup>5</sup> mm/h; and at *<sup>z</sup>* = 96 cm, the maximum value is 4.22 × 105 mm/h. The maximum value decreases from 5.74 × <sup>10</sup><sup>5</sup> mm/h to 4.22 × <sup>10</sup><sup>5</sup> mm/h, which indicates that the maximum rainfall was continuously spreading during the falling process.

**Figure 6.** *Cont*.

**Figure 6.** Rainfall intensity distribution at different section heights on the *x*-*o*-*y* measuring platform. (**a**) *θ*<sup>1</sup> = −30◦, *θ*<sup>2</sup> = 30◦, *β* = 77◦, *z* = 56 cm, and *f* = 0.3; (**b**) *θ*<sup>1</sup> = −30◦, *θ*<sup>2</sup> = 30◦, *β* = 77◦, *z* = 76 cm, and *f* = 0.3; (**c**) *θ*<sup>1</sup> = −30◦, *θ*<sup>2</sup> = 30◦, *β* = 77◦, *z* = 96 cm, and *f* = 0.3.

The maximum rainfall of each section (*q*max) is shown in Figure 7. Note that (1) the maximum rainfall intensity decreases with an increase in distance of *z*, (2) *q*max sharply decreases for *f* < 1.4 and slightly increases for *f* > 1.4, and (3) the results indicated that *q*max decreases with an increase of *β*. In Figure 7a,c, for *<sup>f</sup>* = 1.3 and *<sup>z</sup>* = 96 cm, *<sup>q</sup>*max/*Q*<sup>t</sup> × <sup>10</sup><sup>3</sup> = 1.6 at *<sup>β</sup>* = 48◦, while *<sup>q</sup>*max/*Q*<sup>t</sup> × <sup>10</sup><sup>3</sup> decreased to 0.45 at *β* = 90◦, where *Q*<sup>t</sup> is the sum flow discharge of surface hole and deep hole outlets.

**Figure 7.** The maximum rainfall varies with the collision angles and flow ratios at different sections. (**a**) *θ*<sup>1</sup> = −30◦, *θ*<sup>2</sup> = 0◦, and *β* = 48◦; (**b**) *θ*<sup>1</sup> = −30◦, *θ*<sup>2</sup> = 30◦, and *β* = 77◦; (**c**) *θ*<sup>1</sup> = −30◦, *θ*<sup>2</sup> = 45◦, and *β* = 90◦.

### *3.2. Rainfall Intensity Distribution on the x-Axis*

The rainfall intensity distributions on the *x*-axis are shown in Figure 8. Note that the rainfall intensity data along the *x*-axis was best correlated with the Gaussian distribution. In addition, the maximum rainfall intensity value gradually decreases with an increase along the *z*-axis. For example, in Figure 8a, the height of the section increases from 36 cm to 96 cm, and the maximum rainfall intensity value decreases from 2.48 × <sup>10</sup><sup>5</sup> mm/h to 1.43 × 105 mm/h. With an increase in section height, the rainfall intensity distribution along the *x*-axis flattens, and the range on the *x*-axis is bigger.

**Figure 8.** Rainfall intensity distributions along the *x*-axis and comparison with Equations (1) and (2). (**a**) *θ*<sup>1</sup> = −30◦, *θ*<sup>2</sup> = 30◦, *β* = 77◦, and *f* = 0.9; (**b**) *θ*<sup>1</sup> = −30◦, *θ*<sup>2</sup> = 45◦, *β* = 90◦, and *f* = 3.2.

For two-jet collisions in air on the *x*- and *z*-axis, the rainfall intensity distribution along the *x*-axis after the collision can be expressed by the following formulas (Figure 8):

$$\log I = \lg I\_{\text{max}} \times \varepsilon^{\left(-a \times \left(\frac{x - x\_{\text{max}}}{z}\right)^2\right)} \tag{1}$$

$$\mathfrak{a} = 1.7e^{(-2.1 \times f)} + 2\cos\beta\tag{2}$$

where *I* is the rainfall intensity and *I*max is the maximum rainfall intensity, *α* is a coefficient.

### *3.3. Trajectory Line after the Collision on the x-Axis*

As shown in Figure 9, based on previous research on collision flow [9], the surface water tongue and deep water tongue completely collide in air; then, two water tongues combine and shoot in one direction. *V*<sup>1</sup> is the surface hole water tongue velocity, *θ*<sup>1</sup> is the angle of depression, and *q*<sup>1</sup> is the discharge per unit width; in addition, *V*<sup>1</sup> is the deep hole water tongue velocity, *θ*<sup>2</sup> is the pick angle, and *q*<sup>2</sup> is the discharge per unit width. Two water jets meet at the unit of *M*. At the confluence, the flow velocity of the upper edge of the water tongue at the confluence is *V*1*M*; the angle between the upper edge of the water tongue and the horizontal direction is *β*1; the flow velocity of the inner edge of the water tongue is *V*2*M*; the angle between the inner edge of the water tongue and the horizontal direction is *β*2; the flow velocity of the mixed water tongue after the phase confluence is *VM*, and the angle between the water tongue and the horizontal direction is *βM*.

The spatial position of collision point *M* can be obtained by calculating the trajectories of the inner edge of the surface-hole jet and outer edge of the deep-hole jet. The trajectories of the deep and surface water tongues were obtained from Equation (3).

$$z(\mathbf{x}) = z\_0 + \tan(a)\mathbf{x} - \frac{\mathbf{g}\mathbf{x}^2}{2V\_0^2\cos^2(a)}\tag{3}$$

where *V*<sup>0</sup> is the initial speed, *g* is gravity acceleration.

**Figure 9.** Schematic diagram of the collision between the surface-hole jet and the deep-hole jet in air.

The coordinate systems *x*10*z*<sup>1</sup> and *x*20*z*2, in Figure 9 have the following geometric relationship:

$$
\Delta \mathbf{x}\_1 = \mathbf{x}\_2 + \Delta \mathbf{x}\_0,\\
\mathbf{z}\_1 = \mathbf{z}\_2 + \Delta \mathbf{z}\_0 \tag{4}
$$

Using the projectile principle, the trajectory equation of the surface water jet was:

$$z\_1 = \frac{\mathcal{G}}{2V\_1^2 \cos^2 \theta\_1} \mathbf{x}\_1^2 + \mathbf{x}\_1 t \mathbf{g} \theta\_1 = A\_1 \mathbf{x}\_1^2 + \mathbf{x}\_1 t \mathbf{g} \theta\_1 \tag{5}$$

where *V*<sup>1</sup> is the surface hole jet initial speed, *x*<sup>1</sup> is the transverse distance, *z*<sup>1</sup> is the vertical distance, *<sup>A</sup>*<sup>1</sup> = *<sup>g</sup>* 2*V*<sup>2</sup> <sup>1</sup> cos2 *<sup>θ</sup>*<sup>1</sup> .

The energy velocities of *V*1*<sup>M</sup>* and *V*2*<sup>M</sup>* before the impact of the surface water jet and deep water jet became:

$$V\_{1M} = \varphi\_{\mathbf{a}} \sqrt{2\mathbf{g}(z\_{1M} + \frac{V\_1^2}{2\mathbf{g}})} \tag{6}$$

$$V\_{2M} = \wp\_{\mathfrak{a}} \sqrt{2\lg(z\_{2M} + \frac{V\_2^2}{2\mathfrak{g}})} \tag{7}$$

where *ϕ<sup>a</sup>* is the velocity coefficient of air resistance (*ϕ<sup>a</sup>* ≈ 0.95), *z*1M is the vertical distance from the surface hole to control volume, *z*2M is the vertical distance from the deep hole to control volume. The control volume was taken around point *M*, as shown in Figure 9, and the momentum integral equation and continuous equation of fluid mechanics were used.

$$\iiint\_{M} \frac{\partial \rho}{\partial t} \bar{\mathbf{v}} d\mathbf{V} + \iiint\_{M} \rho \bar{\mathbf{v}} (\bar{\mathbf{v}} \cdot \bar{\mathbf{n}}) d\mathbf{A} = \Sigma \,\, \overline{F} \tag{8}$$

$$\iiint\_{M} \frac{\partial \rho}{\partial t} d\mathbf{V} + \iiint\_{M} \rho (\vec{\boldsymbol{\nu}} \cdot \vec{\boldsymbol{n}}) d\mathbf{A} = 0 \tag{9}$$

where *<sup>ρ</sup>* is the density of water, and <sup>→</sup> *F* is the external force that acts on the control volume. In the absence of air resistance, *β<sup>M</sup>* and *VM* could be obtained with

$$
\tan \beta\_M = \frac{V\_{1M} q\_1 \sin \beta\_1 - V\_{2M} q\_2 \sin \beta\_2}{V\_{1M} q\_1 \cos \beta\_1 + V\_{2M} q\_2 \cos \beta\_2} \tag{10}
$$

$$V\_M = \frac{V\_{1M}q\_1 \cos \beta\_1 + V\_{2M}q\_2 \cos \beta\_2}{q\_M \cos \beta\_M} \tag{11}$$

with the assumption that *V*1*<sup>M</sup>* = *V*2*M*. Thus, Equations (10) and (11) become

$$\tan \beta\_M = \frac{q\_1 \sin \beta\_1 - q\_2 \sin \beta\_2}{q\_1 \cos \beta\_1 + q\_2 \cos \beta\_2} \tag{12}$$

$$V\_M = \frac{V\_{1M}(q\_1 \cos \beta\_1 + q\_2 \cos \beta\_2)}{q\_M \cos \beta\_M} \tag{13}$$

The comparisons of the measured locations of *I*max with Equations (12) and (13) are shown in Figure 10, which indicates good agreement.

**Figure 10.** Comparison of experimental and calculated values of the maximum rainfall intensity points. (**a**) *θ*<sup>1</sup> = −30◦, *θ*<sup>2</sup> = 0◦, and *β* = 48◦; (**b**) *θ*<sup>1</sup> = −30◦, *θ*<sup>2</sup> = 30◦, and *β* = 77◦; (**c**) *θ*<sup>1</sup> = −30◦, *θ*<sup>2</sup> = 45◦, and *β* = 90◦.

Based on the above results, we assume that the velocity magnitude was identical in all directions after the collision, but that the directions vary. Therefore, the velocities in all directions were *VM*. Downstream of the collision, the trajectory boundaries (5% × *I*max) are shown in Figure 11. Hence,

using Equation (3), *α*<sup>1</sup> and *α*<sup>2</sup> can be obtained—*α*<sup>1</sup> is the angle between the *z*-axis and the inner edge line, and *α*<sup>2</sup> is the angle between the *x*-axis and the outer edge line.

**Figure 11.** Comparison of the experimental and calculated values of *L*X1 and *L*X2. *L*X1 is the distance of inner edge water tongue in x axis, *L*X2 is the distance of outer edge water tongue in x axis. (**a1**) *θ*<sup>1</sup> = −30◦, *θ*<sup>2</sup> = 0◦, and *β* = 48◦; (**a2**) *θ*<sup>1</sup> = −30◦, *θ*<sup>2</sup> = 0◦, and *β* = 48◦; (**b1**) *θ*<sup>1</sup> = −30◦, *θ*<sup>2</sup> = 30◦, and *β* = 77◦; (**b2**) *θ*<sup>1</sup> = −30◦ and *θ*<sup>2</sup> = 30◦; (**c1**) *θ*<sup>1</sup> = −30◦, *θ*<sup>2</sup> = 45◦, and *β* = 90◦; (**c2**) *θ*<sup>1</sup> = −30◦, *θ*<sup>2</sup> = 45◦, and *β* = 90◦.

The calculated values of *α*<sup>1</sup> and *α*<sup>2</sup> are shown in Figure 12. Note that (1) with the identical collision angle *β*, *α*<sup>1</sup> increased, and with the increase in flow ratio, *α*<sup>2</sup> decreased; and (2) with the identical flow ratio, *α*<sup>1</sup> increased with the increase in collision angle *β*, and *α*<sup>2</sup> decreased with the increase in collision angle *β*. It can also be seen from Figure 12 that *α*<sup>1</sup> approximately shows exponential growth with the flow ratio, and *α*<sup>2</sup> shows logarithmic growth with the flow ratio. Data analysis indicates that the values of *α*<sup>1</sup> and *α*<sup>2</sup> are obtained as

$$w\_1 = 16e^{(f\sin^2\beta)}\tag{14}$$

$$
\alpha\_2 = \beta \times (-0.182) \times \ln(f) - \beta \times \sin(\beta - 80^\circ) \tag{15}
$$

*Water* **2018**, *10*, 1600

**Figure 12.** Comparison of the values obtained from Equations (14) and (15) with the present data. (**A**) Inner angle of *α*1, and (**B**) outer angle of *α*2.

The computed and measured values correlated to *R*<sup>2</sup> = 0.97 (Figure 12A) and 0.96 (Figure 12B). However, up until now, there was a paucity of information on the rainfall characteristics of the two-jet collisions. This present study provides a new method for understanding the rainfall characteristics of the two-jet collision, which will be useful to extend to other models or even prototypes for further research.

### **4. Conclusions**

In this paper, the flow characteristics of two-jet collisions with different collision angles and flow ratios were investigated systematically by a series of model tests. The present investigation illustrated the development of the spatial distribution of rainfall intensity. Observations indicated that the maximum rainfall intensity sharply decreased with an increase in flow ratio, while the maximum rainfall intensity slightly increased, as the flow ratio was greater than 1.0.

On the *x*-axis, the distribution of rainfall intensity followed a Gaussian distribution. The range of rainfall intensity tended to reach the maximum for the flow ratio = 1.0. A theoretical equation to compute the locations of maximum rainfall intensity was presented, which is in good agreement with the model tests. The formulas to calculate the boundary lines of the *x*-axis were proposed. The present analysis was based upon experiments performed in a model with 48◦ < *β* < 90◦, 0.3 < *f* < 3.2, and 36 cm < *z* < 96 cm. The rainfall intensity distributions on the *y*-axis will be further studied.

**Author Contributions:** Conceptualization, H.Y. and W.X.; Methodology, W.X.; Software, Validation, H.Y.; W.X. and R.L.; Formal Analysis, H.Y.; Experiment, H.Y., Y.F., Y.H.; Writing-Original Draft Preparation, H.Y.; Writing-Review & Editing, H.Y.; Supervision, W.X.; Project Administration, W.X.; Funding Acquisition, W.X.

**Funding:** This work was supported by the National Key Research and Development Program of China (Grant No. 2016YFC0401707) and National Natural Science Foundation of China (Grant 51609162 and Grant 51409181).

**Conflicts of Interest:** The authors declare no conflicts of interest.

### **References**


© 2018 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

### *Article* **Development of Bubble Characteristics on Chute Spillway Bottom**

### **Ruidi Bai 1, Chang Liu 1, Bingyang Feng 1, Shanjun Liu <sup>1</sup> and Faxing Zhang 1,\***

State Key Laboratory of Hydraulics and Mountain River Engineering, Sichuan University, Chengdu 610065, China; bairuidiscu@163.com (R.B.); sculiuchangvip@163.com (C.L.); 18625801588@163.com (B.F.); drliushanjun@vip.sina.com (S.L.)

**\*** Correspondence: zhfx@scu.edu.cn; Tel.: +86-159-2803-8306

Received: 30 July 2018; Accepted: 23 August 2018; Published: 24 August 2018

**Abstract:** Chute aerators introduce a large air discharge through air supply ducts to prevent cavitation erosion on spillways. There is not much information on the microcosmic air bubble characteristics near the chute bottom. This study was focused on examining the bottom air-water flow properties by performing a series of model tests that eliminated the upper aeration and illustrated the potential for bubble variation processes on the chute bottom. In comparison with the strong air detrainment in the impact zone, the bottom air bubble frequency decreased slightly. Observations showed that range of probability of the bubble chord length tended to decrease sharply in the impact zone and by a lesser extent in the equilibrium zone. A distinct mechanism to control the bubble size distribution, depending on bubble diameter, was proposed. For bubbles larger than about 1–2 mm, the bubble size distribution followed a—5/3 power-law scaling with diameter. Using the relationship between the local dissipation rate and bubble size, the bottom dissipation rate was found to increase along the chute bottom, and the corresponding Hinze scale showed a good agreement with the observations.

**Keywords:** chute aerator; spillway bottom; air concentration; air bubble frequency; air bubble chord length

### **1. Introduction**

Chute aerators, regarded as a cost-effective method, have been widely used owing to their proven success as a countermeasure in the Grand Coulee dam since 1960 [1]. Previous studies have provided a significant amount of information on the macroscopic air-water flow properties [2–13], such as the air concentration, air entrainment and detrainment. Rutschmann and Hager [2], combining both model and prototype data, proposed two approaches to calculate the air entrainment coefficient *β* = *qa*/*q*, where *qa* = air discharge and *q* = water discharge. The air entrainment coefficient was also investigated by, among others, Shi [14], Pinto [15,16], and Low [17]. Chanson [3] investigated the characteristics of chute aerator flow on the distributions of various bubble parameters (void fraction, velocity, air discharge, turbulent diffusivity, etc.) for different inflow Froude numbers. He also provided a better understanding of the air entrainment processes and found a strong detrainment process occurred in the impact zone. Kramer et al. [5] divided the chute aerator into five zones: inflow region; air detrainment region; minimum air concentration region; air entrainment region; and uniform mixture flow region, and systematically investigated the air detrainment gradients in the far-field of chute aerators. Pfister and Hager [6–8] measured the air concentration and divided the chute aerator into three zones: jet zone, reattachment and spray zone, and far–field zone. The authors also discussed the air entrainment and proposed two comprehensive methods to calculate the average and bottom air concentration development using deflectors and offsets along the chute. Avoiding the upper aeration effect, Bai et al. [9] divided the chute aerator into three zones and found that an

intensive air detrainment occurred in the impact zone and the detrained air escaped into the cavity zone. The previous experimental investigations are summarized in Table 1.


**Table 1.** Experimental investigations of Chute Spillways.

Many studies concerning the minimum or critical air concentration to avoid cavitation erosion have been conducted [18,19]. Peterka [18] found that when the average air concentration was between 0.01 and 0.06, no cavitation erosion was observed. Up to now, there is no reliable design guideline for the distance required between two aerators to produce sufficient bottom air concentration [5]. Based on prototype observations on Russian dams, Semenkovand Lantyaev [20] provided an average bottom air concentration decay of 0.40% to 0.80% per chute meter.

However, air concentration governs the comprehensive performance of the bubbles, including their chord length and frequency form. In particular, Robinson [21], Xu et al. [22], and Brujan [23] adopted advanced techniques to investigate the mechanism of erosion damage and found that the interaction between the air bubbles and cavitation bubbles was the mechanism that assisted in the prevention of erosion damage. These studies indicate that the bubble size and frequency are significant parameters for the mechanism of erosion damage. However, negligible research has been conducted on the spillway bottom air microcosmic characteristics such as the bottom bubble frequency and chord size. Such investigations are difficult because of the complex nature of the flow and strong impaction of the bottom. In the present study, the development of the bottom bubble characteristics was measured and analyzed.

### **2. Hydraulics Model**

The experiments were conducted in a rectangular chute that was 0.25-m wide with a variable bottom slope 5.71◦ ≤ *α* ≤ 18.2◦ (Figure 1). The measured section was 5-m long and consisted of a smooth convergent nozzle with a width of 0.25 m and height of 0.15 m. The approach flow depth was *h*<sup>0</sup> = 0.15 m and the emergence angle was 0◦ ≤ *θ*<sup>0</sup> ≤ 14.1◦. The water discharge was supplied by a large tank (3.1-m wide, 3.2-m long, and 9.2-m high) and was measured using a rectangular sharp-crested weir with an accuracy of approximately 1%. In order to gain sufficient velocity, the nozzle was located in the bottom of the tank. A range of approach flow velocities was set from 4.0 m/s to 9.0 m/s. Observations indicated that the approach of the flow was stable. Hence, the Froude number F0 (*V*0/(*gh*0) 0.5) was within the range of 3.3 to 7.4 (Table 2), where *g* was the gravity. *x* was the streamwise coordinate along the chute bottom, *hs* was the offset height, *L* was the cavity length, and *Lm* and *LD* were the specific lengths [9].

(**b**)

**Figure 1.** Definition sketch with relevant parameters (**a**) and side view of chute aerator flow with high speed camera (**b**).

**Table 2.** Experimental conditions of chute aerator flow.


All measurements were performed on the channel centre line. The air-water flow properties were measured with a phase-detection needle probe. The working principle of a conductivity probe was based on the difference in the voltage indices at the platinum tip between the air and water phases. The response time of this sensor was less than 10 μs. The signals from the conductivity probe were recorded at a scan rate of 100 kHz per for a 40 s scan period. The accuracy of the air concentration was Δ*C/C* = 3%, and the bubble frequency was Δ*f*/*f* = 1%, where *C* was the air concentration and *f* was the bubble frequency. The length of the impact zone was very short. A few studies in the impact zone have been performed with a space of 0.03 m along the chute bottom.

### **3. Basic Chute Aerator Flow Properties**

### *Air Concentration and Bubble Frequency Profiles*

A typical photo of the offset-aerator flow is shown in Figure 2, where four zones were introduced. The typical distributions of the air concentration and bubble frequency in the different zones are depicted in Figure 2 as a function of the dimensionless distance perpendicular to the bottom, i.e., *z*/*h*0. Owing to a sufficient depth of the approach flow, the unaerated water prevented the upper aeration effect.

**Figure 2.** Distributions of air concentration (**a**) and bubble frequency (**b**) along the chute (S5: *V*<sup>0</sup> = 9 m/s, F0 = 7.4).

In the cavity zone (0 < *x*/*L* < 1), air entrainment occurs in the form of air bubbles entrapped into the lower flow and the air concentration decreased from the lower surface to the vertical direction. At *x*/*L* = 1, the lower jet impacted on the chute bottom. Further downstream (*x*/*L* > 1), the distributions of air concentration were different from the cavity zone. Figure 2a shows that the air concentration distributions downstream of the cavity zone, exhibiting an air concentration peak in the turbulent shear region, whereas the air concentration increases from the chute bottom to its maximum value, and then decreases in the vertical direction.

Note that: (1) for *x* > *L*, the maximum value air concentration *Cm* decreased sharply in the impact zone, while in the equilibrium zone, the peak air concentration *Cm* decreased slightly. (2) Experimental results indicated the existence of a maximum value *fm* in bubble frequency (Figure 2b), but its location differed from the locus of the maximum void fraction as noticed by Bai et al. [10]. In the cavity, the location of *fm* approached *C* = 0.50 with the development of the low jet aeration. The location of *fm* was approximately the same as that of the section maximum air concentration *Cm* at a section in the impact zone, but lower than that in the equilibrium zone.

For more than 25 experiments conducted in this model, only one set of experiments is shown herein to illustrate the basic chute aerator flow property [9]. However, equations presented hereafter are based on all experiments.

### **4. Bottom Air Bubble Characteristics**

The entrained air from the chute aerator aims to prevent erosion damage. Recent investigations have indicated that the bubble size and frequency are the significant factors preventing erosion damage [22,23]. The air bubbles generated in the cavity zone travelled through the impact zone to continue downstream. For the chute aerator flows, a large amount of air entrained in the cavity zone was only partially transported to the downstream flow. Under the effect of rollers, large bubbles broke into small bubbles before being transported downstream zone simultaneously in the impact zone.

### *4.1. Bottom Air Concentration and Bubble Frequency*

Bottom air concentration *Cb* was defined experimentally using the values measured closest to the model chute bottom (*z* = 1–1.5 mm). The experimental data of bottom air concentration *Cb* is presented in Figure 3. Note the following: (1) the decay rate of *Cb* was different in the impact (*L* < *x* < *Lm*) and equilibrium zones (*Lm* ≤ *x* < *LD*); *Cb* decreased sharply in the impact zone and decreased slightly in the equilibrium zone. In the impact zone, bottom air concentration *Cb* = 0.64 at *x*/*L* = 1.06 and it decreased to 0.18 at *x*/*L* = 1.16. In contrast, the bottom air concentration decreased slightly in the equilibrium zone, with concentration *Cb* = 0.145 at *x*/*L* = 1.31 decreasing to 0.064 at *x*/*L* = 1.91. (2) With the increase in F0, the values of *Cb* increase in the present experiments, in agreement with the results of Pfister and Hager [7]. (3) The comparison of the measured values of *Cb* with the values from the formula proposed by Pfister and Hager [7] is shown in Figure 3, which displays a similar trend downstream of the chute aerator. However, the present values are generally larger than those of Pfister and Hager [7]. At a sufficient water depth, the detrained air bubbles cannot escape into the atmosphere through the upper surface; it can only escape into the cavity zone, which may cause the level of *Cb* to be larger than that of previous studies. In addition, for a shorter length of the impact zone, the fewer measurement sections may cause an insufficient sharp decrease. (4) It can be concluded that the air transportation regularity is different in the impact and equilibrium zones.

(**a**)

**Figure 3.** *Cont*.

**Figure 3.** Bottom air concentration along the chute and compared with Pfister [7] (**a**): S2; (**b**): S5.

From the experimental data analysis, bottom air concentration *Cb* can be expressed by the following equation (Figure 3):

$$\frac{\mathcal{C}\_b}{(\mathcal{C}\_b)\_L} = \left(\frac{\mathcal{X}}{L}\right)^{0.985F\_0 - 17.5} \text{ for } L < \text{ x } < \text{ } L\_{m\prime}\text{ R}^2 = 0.89 \tag{1}$$

$$\frac{\mathcal{C}\_b}{(\mathcal{C}\_b)\_{L^m}} = \left(\frac{\mathbf{x}}{L\_m}\right)^{0.071F\_0 - 1.73} \text{ for } \mathbf{x} > L\_{m\_\prime} \ R^2 = 0.97\tag{2}$$

where *R* is the correlation coefficient, (*Cb*)*<sup>L</sup>* is the bottom air concentration at *x* = *L*, (*Cb*)*Lm* is the bottom air concentration at *x* = *Lm*.

The experimental data of bottom bubble frequency *fbh*0/*V*<sup>0</sup> along the chute are displayed in Figure 4. There are a few reports on the chute bottom bubble frequency of the lower jet downstream of the chute aerator. Note the following: (1) *fbh*0/*V*<sup>0</sup> decreased in both the impact and equilibrium zones, whereas decay rate of *fbh*0/*V*<sup>0</sup> differed slightly. For F0 = 7.4, *fbh*0/*V*<sup>0</sup> = 11.2 at *x*/*L* = 1.06, and it decreased to 10.88 at *x*/*L* = 1.17 in the impact zone (Figure 4d). In contrast, the bottom bubble frequency decreased slightly more in the equilibrium zone, with *fbh*0/*V*<sup>0</sup> = 9.44 at *x*/*L* = 1.31 decreasing to 6.22 at *x*/*L* = 1.91. (2) In contrast, *fbh*0/*V*<sup>0</sup> decreased slightly in the impact zone, and there was a large variation in air bubble concentration *Cb*. (3) Furthermore, *fbh*0/*V*<sup>0</sup> increased with increasing F0.

**Figure 4.** *Cont*.

**Figure 4.** Bottom air bubble frequency along the chute (**a**): S2; (**b**): S3; (**c**): S4; (**d**): S5.

Overall, the present data could be described by

$$\frac{f\_b h\_0}{V\_0} = 2.81 F\_0 \times \exp\left(-\frac{\text{x}}{L} F\_0^{-0.27}\right) \text{ for x} > L, \text{ } R^2 = 0.94 \tag{3}$$

where *R* was the correlation coefficient.

### *4.2. Bottom Bubble Chord Length Distributions*

The chord length distribution of the air bubbles shows microscopic variations passing a measurement point in the stream-wise direction. For all the flow conditions, the experimental results exhibited a broad range of the chord length, varying from less than 0.1 mm to more than 30.0 mm. In this study, the chord length interval was 0.1 mm. The chord lengths are not the bubble diameters, but are the characteristic stream-wise bubble sizes.

Garrett et al. [24] and Deane et al. [25] indicated that the bubble chord length distribution *Na* depended on the dissipation rate *ε* and bubble diameter *dab*, and the scaling law followed as

$$N\_{\rm a} \propto A \varepsilon^{1/3} d\_{ab}^{-5/3} \tag{4}$$

where A was a constant, *Na* = *Nt*/*t*, *Nt* was the number of air bubbles detected in the scan period *t* (*t* = 40 s).

The distributions of the bottom bubble chord length in the impact and equilibrium zones are depicted in Figure 5. In the impact zone, the air concentration was within the range of 0.10 < *C* < 0.90, such that the intermediate region neither corresponded to the high-air concentration region nor to the bubbly flow region. In the equilibrium zone, the maximum air concentration was less than 0.3, so that the zone corresponded to the bubbly flow region. Note the following: (1) the distributions of bottom bubble chord length were skewed with a large proportion of small bubbles. Compared with the impact region, the range of bubble chord length distributions decreased. In the impact zone, the largest probability of bubble chord lengths was concentrated in the range 0.1 ≤ *dab* ≤ 30 mm, while in the equilibrium zone, the highest probability of bubble chord lengths was concentrated in the range 0.1 ≤ *dab* ≤ 10 mm. (2) The curve of distribution of the bottom bubble chord sizes was also sharper and narrower along the chute. The power-law scaling of bubble density on diameter is about *β =* −5/3 for larger than about 1–2 mm bubbles in the impact and equilibrium zones. (3) In the impact zone, *β* showed a slight increase from −1.2 at *x*/*L* = 1.06 to −5/3 at *x*/*L* = 1.15 due to the air detrainment. In the equilibrium zone, *β* showed a slightly increase from –5/3 at *x*/*L* = 1.21 to −2.4 at *x*/*L* = 1.64 due to the buoyancy effect.

**Figure 5.** Bottom bubble chord length distributions in the impact zone (**a**) and equilibrium zone (**b**) (S5, *V*<sup>0</sup> = 8 m/s, F0 =6.6).

### *4.3. Bottom Turbulent Dissipation Rate*

Xu et al. [22,23] and Brujan et al. proved that the size of air bubbles near the bottom was a significant parameter for the collapse process of cavitation bubbles. There was little information to investigate the air bubble variations on the chute bottom. The equilibrium zone belonged to the bubbly flow region and the dynamic behavior of the flow was considered to be a steady state air-water two-phase flow.

Hinze [26] indicated that any fluid gas bubbles or droplets were likely to fragment if the pressure forced on its surface exceeds the restoring forces associated with surface tension.

$$
\rho \frac{\mu^2}{\sigma} d\_\varepsilon = \mathcal{W}\_\mathbb{C} \tag{5}
$$

where W*c* = the critical Weber numbers, *u* = the turbulent velocity field on the scale of the bubble. If the fluctuating velocity field was assumed to be described by Kolmogorov's inertial subrange, so that *u*<sup>2</sup> = 2*ε*2/3*d*ab2/3, then only bubbles with diameters larger than the Hinze scale would fragment [17–19].

$$d\_c = \left(\frac{\mathcal{W}\_{\rm c} \sigma}{2\rho}\right)^{3/5} \varepsilon^{-2/5} \tag{6}$$

where *dc* = the Hinze scale, *ε* = turbulent dissipation rate, which was estimated rather than measured directly. Recent experiments suggest that Wc lied in the range 0.585–4.7, and we have used the value Wc = 1.0 (Table 3). For this to exceed a critical value, Wc required that *d* > *dc*.


**Table 3.** Critical Weber numbers for the splitting of air bubbles in water flows.

If the Hinze scale was known, the turbulent energy dissipation rate *ε* on the chute bottom could be estimated from Equation (7). Hinze [26] used the definition that 95% of the air is contained in bubbles with a diameter less than *dc* and the calculated turbulent energy dissipation rate *ε* for Deane's [18] data was much smaller. Garrett et al. [24] used the results of Martínez-Bazán et al. [30,31] and established a one-to-one relationship between the local dissipation rate and bubble size.

$$\varepsilon = \left(1/2\right)^{5/2} \left(\frac{\sigma}{\rho}\right)^{3/2} \int\_0^\infty \left(\frac{d\_{ab}}{2}\right)^{1/2} N\_d \mathbf{d}d\_{ab} / \int\_0^\infty \left(\frac{d\_{ab}}{2}\right)^3 N\_d \mathbf{d}d\_{ab} \tag{7}$$

The bottom dissipation rate *ε* calculated from Equation (7) is shown in Figure 6 and tends to increase in the measured distance. From Equation (6), a change in dissipation rate *ε* from 0.10 to 0.75 m<sup>2</sup> s−<sup>3</sup> corresponded to a change in Hinze scale from 5.5 to 2.4 mm. Thus, the estimate of 2.4 mm for the Hinze scale was in good agreement with the observed change in power-law scaling of the bubble distribution at around 1–2 mm (Figure 5b). Because bubbles smaller than the Hinze scale are stabilized by surface tension forces, a reasonable estimate for this scale is the diameter of the smallest bubbles observed to fragment. In the present study, the bubble chord lengths were smaller than the bubble diameters, due to the measurement technique utilizing conductivity probes. A challenge remains to observe the bubble diameter in a model facility, or even a prototype facility.

**Figure 6.** Calculated values of turbulent energy dissipation rate ε on the chute bottom.

### **5. Conclusions**

This study presented new information on the bottom air bubble characteristics in the downstream of chute offset-aerators, which were measured with an accurate measurement technique utilizing conductivity probes. For eliminating the upper aeration, the bottom concentration, bubble frequency, bubble chord length, and bottom dissipation rate were analyzed along the chute longitudinal direction. The present study complemented the existing studies on the bubble properties of chute aerator flows.

The bottom air concentration was found to decrease sharply in the impact zone and to decrease slightly in the equilibrium zone. The bottom air bubble frequency was found to decrease slightly along the chute bottom. A formula to predict the bottom bubble frequency in the impact and equilibrium zones was proposed.

Observations showed that the range of probability of bottom bubble chord lengths tended to decrease sharply in the impact zone and to a lesser extent in the equilibrium zone. The bubbles larger than the Hinze scale were subject to fragmentation and showed −5/3 power-law scaling with diameter. Using the relationship between the local dissipation rate and bubble size [24], the bottom dissipation rate was found to increase along the chute bottom. The Hinze scale estimated from Equation (6) showed a good agreement with the observed change in power-law scaling of the bubble distribution.

**Author Contributions:** Conceptualization, R.B., F.Z. and S.L.; Methodology, R.B. and F.Z.; Software, C.L.; Validation, R.B.; Formal Analysis, R.B. and F.Z.; Investigation, R.B.; Resources, S.L.; D.C., B.F.; Writing-Original Draft Preparation, R.B.; Writing-Review & Editing, R.B.; Visualization, B.F.; Supervision, F.Z.; Project Administration, R.B., F.Z. and S.L.; Funding Acquisition, R.B., F.Z. and S.L.

**Funding:** This work was supported by the National Natural Science Foundation of China (Grant No. 51709293 and Grant No. 51679157), the Fundamental Research Funds for the Central Universities (Grant No. 20826041A4305), and the National Key Research and Development Program (Grant No. 2016YFC0401707).

**Conflicts of Interest:** The authors declare no conflict of interest.

### **References**


© 2018 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

### *Article* **Experimental Study on the Impact Characteristics of Cavitation Bubble Collapse on a Wall**

### **Jing Luo, Weilin Xu \*, Jun Deng, Yanwei Zhai and Qi Zhang**

State Key Laboratory of Hydraulics and Mountain River Engineering, Sichuan University, Chengdu 610065, China; luojing@scu.edu.cn (J.L.); djhao2002@scu.edu.cn (J.D.); zaiyanwei0422@163.com (Y.Z.); scu\_zq@163.com (Q.Z.)

**\*** Correspondence: xuwl@scu.edu.cn; Tel.: +86-028-8540-1301

Received: 20 August 2018; Accepted: 14 September 2018; Published: 15 September 2018

**Abstract:** As a hydrodynamic phenomenon, cavitation is a main concern in many industries such as water conservancy, the chemical industry and medical care. There are many studies on the generation, development and collapse of cavitation bubbles, but there are few studies on the variation of the cyclic impact strength on walls from the collapse of cavitation bubbles. In this paper, a high-speed dynamic acquisition and analysis system and a pressure measuring system are combined to study the impact of a cavitation bubble generated near a wall for various distances between the cavitation bubble and the wall. The results show that (1) with the discriminating criteria of the impact pressure borne by the wall, the critical conditions for the generation of a micro-jet in the collapse process of the cavitation bubbles are obtained, and therefore collapses of cavitation bubbles near the wall are mainly divided into primary impact area collapses, secondary impact area collapses and slow release area collapses; (2) it can be seen from the impact strength of the cavitation bubble collapse on the wall surface that the impact of cavitation bubbles on the wall surface during the first collapse decreases as *γ* (the dimensionless distance between the cavitation bubble and the wall) increases, but the impact of the second collapse on the wall surface increases first and then decreases sharply. When *γ* is less than 1.33, the impact on the wall surface is mainly from the first collapse. When *γ* is between 1.33 and 2.37, the impact on the wall surface is mainly from the second collapse. These conclusions have potential theoretical value for the utilization or prevention and control technologies for cavitation erosion.

**Keywords:** cavitation bubble; high-speed photography; impact pressure; micro-jet

### **1. Introduction**

When the collapse of cavitation bubbles occurs near a wall, the wall will bear the cyclic impact from the collapse of the cavitation bubble. The characteristic of the interaction between the cavitation bubble and the wall has become the focus of many technical fields, such as water conservancy, shipbuilding, the chemical industry and many other industries.

Theoretical study about cavitation bubbles can be tracked back to the last century [1–3]. The study of the cavitation field mainly focuses on the damage and destruction of a solid wall surface upon cavitation bubble collapse [4,5].

Shaw et al. [6] used a laser to generate a cavitation bubble and combined this with Schlieren photography to conduct a study on the impact performance of the impact wave and micro-jet generated during the first collapse of a cavitation bubble on a wall. Adopting high-speed photographic equipment with a camera speed of up to 100 million frames per second to study the evolution properties of cavitation bubbles around a wall, Lindau and Lauterborn [7] discovered the annular form during the collapse of a cavitation bubble, the formation of a contrajet and the shock wave. Liu et al. [8] used the method of laser-inducing a cavitation bubble to conduct a study on the impacts of a micro-jet caused by laser-induced cavitation bubble collapse on a metallic copper surface, and found that the liquid-jet was the main damage mechanism in cavitation erosion. Xu et al. [9] studied the cavitation bubble dynamic characteristics on an aluminum surface by using the technology of a sensitive fiber-optic sensor based on the optical beam deflection principle. Ren et al. [10] used a hydrophone to conduct a study on the strength mechanism of impact from cavitation collapse within the range of the relative bubble–wall dimensionless distance *γ* (*γ* = *h/Rmax*, *Rmax* is the maximal radius of the cavitation bubble, *h* is the distance between the center of the cavitation bubble and the wall) which is between 0.5 and 2.5. The above studies on the collapse of cavitation bubbles either focus on the first collapse or adopt the indirect method of hydrophone measurement. In addition, from the research on the impact of the micro-jet or high-velocity droplet on a wall with the use of high-speed photography, only qualitative rules can be obtained. However, measuring the microscopical area of the wall under the impact of the micro-jet or high-velocity droplet with the pressure test system provides a new quantitative method for solving millimeter-sized impact-related issues [11].

In recent years, some new research directions on cavitation have been proposed, such as extracorporeal shock-wave lithotripsy [12] and ultrasonic cleaning [13,14]. In addition, studies on the interaction between cavitation bubbles and cavitation bubbles, cavitation bubbles and air bubbles, cavitation bubbles and particles, and cavitation bubbles and ice blocks have also been reported. Fong et al. [15] adopted the underwater discharge technique to induce multiple cavitation bubbles so as to study the coupling effect between two cavitation bubbles generated at different instants. Pain et al. [16] studied the jet flow in an air bubble induced by a nearby cavitation bubble and discovered that the speed of such a jet flow can be up to 250 m/s. Luo et al. [17,18] used the method of high-voltage discharge technology to generate a cavitation bubble to study the interaction between a cavitation bubble and an air bubble, and the interaction between two cavitation bubbles and a rigid wall. Goh et al. [19] studied the interactions between cavitation bubbles horizontally placed below the underwater slab and air bubbles attached to it and discovered that the ratio between the oscillation time of the air bubbles and the oscillation time of the cavitation bubbles is the important parameter impacting the micro-jet caused by cavitation bubble collapse. The interaction between a cavitation bubble and a particle was researched by adopting a high-voltage discharge technology to induce a cavitation bubble [20] and low-voltage underwater discharge technology to induce a cavitation bubble [21]. The interaction between ice and cavitation bubbles was researched by Cui et al. [22], and the direction of the micro-jet and the propagation of shock waves were captured.

Due to the very small volume and short life cycle of cavitation bubbles during the evolution, great difficulties have been imposed on related studies. Therefore, numerical simulation studies have been performed by various researchers. Coupled boundary element and finite element methods have been performed on bubble–structure interactions [23–25]. The boundary element method is used for research on the toroidal configuration and jet impact characteristics in the later stage of the collapse of cavitation bubbles in the inviscid and irrotational field [26–28] numerically simulated the late rebound phenomena of a cavitation bubble by the method of vortex lines. Instead of using a vortex surface, Klaseboer et al. [29] simulated the ring stage of a cavitation bubble with a vortex ring, and profoundly understood the dynamic characteristics of a cavitation bubble and the interaction with walls for different distances between the bubble and the wall, which was in agreement with the experimental results.

For the past few decades, many researchers have tried to find out the impulses associated with cavitation bubble collapse. In most existing studies, these results are only qualitative due to the non-uniform spatial sensitivity of the transducer. Therefore, in this paper, cavitation bubbles are induced by adopting a low voltage spark-discharge system, the cavitation bubble evolution and corresponding wall pressure change processes are measured simultaneously by the piezoresistive pressure sensor. The impact strength of the first collapse and rebound–regeneration collapse of the cavitation bubble on the wall surface is mainly analyzed with the high-speed photographic image and pressure testing data of the wall surface, and the influence of the dimensionless distance of the bubble wall on the impact strength of the wall surface is revealed.

### **2. Experimental Setups**

### *2.1. Cavitation Bubble Induction Device*

In order to get the impact characteristics of the cavitation bubble collapse on the wall, the experimental setup requires a cavitation bubble induction device, a high-speed dynamic collection and analysis system and a transient pressure testing system, as shown in Figure 1.

**Figure 1.** Schematic diagram of experimental facilities.

There are three ways to induce cavitation bubbles in quiescent water: by focusing a laser, by a high-frequency ultrasound wave or by an underwater discharge. The laser technique to induce a cavitation bubble has many advantages, such as its higher accuracy in controlling the focal spot, but the focusing method of the laser is rather complex and costly [30]. In the current work, the method of low-voltage underwater discharge was adopted. The method can accurately control the position and size of the cavitation bubble. To minimize possible interference, the diameter of the electrodes used is 0.15 mm, and the bubbles created have a maximum radius of 30 times the diameter of the electrodes. The copper electrodes (length 50 mm) are fixed in the trestles. The trestles are placed in a water tank (40 cm × 30 cm × 40 cm), which is filled with twice deionized water (at a constant temperature of 22 ◦C). To create the bubble, two capacitors are fully charged to the specified voltage (the adjustable range is 0 V~100 V) through a 1 KΩ resistor, then switched to the discharge circuit, and the capacitor is discharged through the pair of crossed electrodes.

The life cycle of cavitation bubbles is very short; in order to study the impact process of the collapse of the cavitation bubble on the wall surface, a high-speed dynamic collection and analysis system and a transient pressure testing system must be adopted for synchronized recording. The high-speed dynamic collection and analysis system is composed of a high-speed camera, lens assembly, lighting equipment, etc., and the transient pressure testing system is composed of the resistance-type pressure sensor and the data acquisition instrument. These two systems are respectively connected to a computer through a synchronous trigger device. The Fastcam SA-Z high-speed camera (maximum acquisition rate: 1,000,000 fps, Photron Inc., Tokyo, Japan,) is used. The radius is of millimeter-size when the cavitation bubble generated after the discharge between electrodes in the water body evolves to its maximum volume. Thus, in order to obtain the clear outline of the cavitation bubbles, a micro-lens and LED lamp (150 W) must be adopted in the acquisition system and synchronized with the high-speed photography system. A framing rate of 150,000 frames per second is used to capture the bubble dynamics in this paper, with a shutter speed of 5.06 μs. The width of the frame is 256 pixels.

### *2.2. Measurement of Cavitation Bubbles During the Evolution Process and the Impact Process*

Considering the size of the cavitation bubbles (the radius of cavitation bubbles produced from the system in the experiment ranges from 5.00 mm to 12.00 mm) and the time for the impact of the collapse of the bubble on the wall surface being very short, the size of the sensor in the transient pressure testing system must be small enough, and the sampling frequency and accuracy must be high enough; in addition, by observing the high-speed photographic images, the time for the impact of the micro-jet on the wall due to the collapse of the cavitation bubbles induced by the experimental system lasted about 100 μs. The minimum rise time was at least 10 μs. The pressure sensor used in the experimental system has a high-frequency bandwidth (200 kHz). Therefore, with the existing sensor technology, the piezoresistive pressure sensor used in this paper can meet the requirements for the measurement of the impact of the collapse of cavitation bubbles on the wall surface. The radius of the sensing area of the sensor adopted in the test is 1.50 mm (the maximum measuring range is 40 MPa). The sampling frequency is 2 MHz, and the accuracy is 0.5% of the maximum measuring range.

The method of image gray processing is used to obtain the maximum radius of the cavitation bubble. In this paper, a circle radius with the same number of the pixels is regarded as the equivalent radius *Req* of the cavitation bubble (*Rmax* is *Req* when the cavitation bubble radius reaches maximum), and the centroid of the cavitation bubble in the high-speed photographic image is regarded as the center of the cavitation bubble. The dimensionless distance *γ* is the ratio of the distance *h* between the centers of the cavitation bubble and pressure sensor and the maximum value of the maximal equivalent radius *Rmax* of the cavitation bubble. Figure 2 shows the data processing process: therein, the scatter diagram shows the relationship between the first cycle of the cavitation bubble time and the corresponding *Rmax*, and it can be seen that this data processing method provides good reproducibility.

**Figure 2.** Processing method of characteristic parameters. (**a**) Cavitation bubble photograph; (**b**) characteristic parameters.

### **3. Results and Discussion**

### *3.1. Influence of Distance Between the Cavitation Bubble and Wall on the Development of Micro-Jetting*

A selected set of high-speed photographic images of the cavitation bubbles in the collapse stage for different *γ* (dimensionless distances between the cavitation bubble and the wall) are shown in Figure 3. Therein, the wall surface is located at the bottom of the image. In the three sets of tests, the dimensionless distances *γ*s are 2.91, 1.91 and 0.83, respectively, the maximum values *Rmax* of the

equivalent radius *Req* are 9.40 mm, 9.31 mm and 9.47 mm respectively, and the distances *h* from the center of the cavitation bubbles to the wall's surface when the cavitation bubbles expand to their maximum radii are 27.35 mm, 17.76 mm and 7.86 mm, respectively.

**Figure 3.** High-speed photographic images of the cavitation bubbles in the contraction collapse stage under the conditions of different distances *γ* (frame-rate: 150,000 fps; exposure time: 5.06 μs; frame width: 32.00 mm).

When the cavitation bubble expands to the maximum radius in the testing, as shown in Figure 3 (a1,b1,c1), the cavitation bubbles begin the contraction stage through an action of the external water body. With the gradual contraction of the cavitation bubble, due to the existence of the wall under the cavitation bubble, the surrounding water begins to fill the surrounding space released by the cavitation bubbles in their contraction process, and the upper and lower surfaces of the cavitation bubble begin to form asymmetrical shapes gradually. The wall surface under the cavitation bubble blocks the filling of water, resulting in a slow contraction speed of the lower surface of the cavitation bubble, while the upper surface belongs to the unbounded domain and contracts quickly. Thus, this asymmetric contraction gradually begins to form from the upper surface (as shown in a9, b9 and c11 in Figure 3). The development of the asymmetric collapse of the surface has been accompanied by the shrinkage of cavitation bubbles to the minimum volume (as shown in a11, b11 and c13 in Figure 3). The cavitation bubbles will rebound–regenerate after they shrink to the minimum volume, and as shown in Figure 3, the surface of the rebound cavitation bubble is not very smooth (as shown in a12 and b12 in Figure 3), but just as in the first cycle, the cavitation bubble will again undergo the expansion–contraction process. In the second expansion–contraction process, the surface is not smooth enough, but this non-smoothness is not enough to change the collapse direction of cavitation bubbles. In Figure 3b, the rebound cavitation bubbles move quickly toward the wall surface in the expansion–collapse process and impact onto the wall in the second collapse.

From the above three tests, the following can be concluded: when the dimensionless distance *γ* is small, the micro-jet directly impacting the wall surface can be formed during the first asymmetric contraction of the cavitation bubble; as the dimensionless distance *γ* has been increasing continuously, the cavitation bubbles will move to the wall surface quickly and will eventually impact on the wall surface. With the further increase of *γ*, cavitation bubbles do not impact on the wall surface through two or three times of contraction. Therefore, the collapse of cavitation bubbles near the wall surface can be divided into main impact area collapse, secondary impact area collapse and slow release area collapse as per the parameters of the distance between the cavitation bubble and wall.

### *3.2. Impact of Cavitation Bubble Collapse on Wall*

The collapse position of cavitation bubbles near the wall surface can be obtained directly through high-speed photography for the development of the micro-jet under the above different *γ* conditions, but the impact strength of cavitation bubbles on the wall surface during the first collapse, whether the impact can reach the wall surface during the second collapse, and the impact strength of the second collapse on the wall surface cannot be ascertained. With these questions, the high-speed dynamic collection and analysis system and transient pressure testing system are combined in this part, obtaining the impact process of the first and second cavitation bubble collapses on the wall surface (that is, the main impact area and secondary impact area). Figure 4 shows the impact process of the two groups of cavitation bubble collapses on the wall surface when the dimensionless distances *γ* are 0.91 and 1.79 respectively. The wall surface is located at the bottom of the image, and the pressure sensor is placed inside of the wall surface.

The maximum radius of the cavitation bubble in Figure 4a is 9.98 mm, and when the cavitation bubble reaches the maximum volume, the distance from the center to the wall surface is 9.08 mm. When the electric spark discharges into the water during the induction of the cavitation bubble, the pressure value appears to fluctuate slightly, and the pressure returns to a normal condition after discharge. The subsequently generated cavitation bubbles will begin the expansion–contraction–collapse–rebound stage. In the cavitation bubble development process, the pressure on the wall surface is stable; when the cavitation bubble contracts to the minimum volume (the first expansion-contraction cycle is 2.91 ms), the pressure on the wall surface sharply increases, abd from the contents of the previous section it can be seen that because of the short dimensionless distance *γ*, the collapse generated due to asymmetrical cavitation bubble contraction directly forms the micro-jet impacting the wall surface (shown as the image in Figure 4a), and the maximum peak

value of impact pressure of the first cavitation bubble collapse on the wall surface is 19.37 MPa. After that, the peak pressure gradually decreases when the cavitation bubble is in the rebound–regeneration stage, the pressure on the wall surface is relatively stable for the whole rebound regeneration stage, and the pressure borne on the wall surface greatly increases again when the cavitation bubble contracts to the minimum volume again. The time period of the whole regeneration rebound is 1.23 ms; the wall surface bears the collapse impact of the rebound cavitation bubble again, and by this time, the maximum pressure on the wall surface is 6.50 MPa. It can be seen from the whole impact process of the cavitation bubbles on the wall surfaces in the condition of *γ* = 0.91 that the strength of the two impacts of cavitation bubbles on the wall surfaces will gradually reduce from 19.37 MPa in the first impact to 6.50 MPa in the second impact.

**Figure 4.** The impulsive process generated by cavitation bubble collapse for (**a**) *γ* = 0.91; (**b**) *γ* = 1.79.

The maximum radius of the cavitation bubble in Figure 4b is 11.02 mm, and the distance from the center to the wall surface is 19.73 mm when the cavitation bubble reaches the maximum radius. The first peak and the second peak are the impact on the wall surface during the first collapse and rebound collapse of the cavitation bubble, with the peak pressure values of 6.38 MPa (2.64 ms) and 21.89 MPa (4.01 ms), respectively. The minimum distance to the wall surface when the cavitation bubble shrinks to the minimum volume for the first time is 13.81 mm, and at this moment, the dimensionless distance from the center of the cavitation bubble to the wall surface is 5.06. The cavitation bubble will quickly move to the wall surface in the rebound evolution process; when the cavitation bubble contracts to the minimum volume once more, the cavitation bubble is close to the wall surface, which is as shown in Figure 4b.

It can be clearly seen from Figure 4 that when the dimensionless distance *γ* is small, the impact of the first collapse of the cavitation bubble on the wall surface is greater than the impact during rebound collapse; when *γ* can meet a certain condition and the cavitation bubble completes the impact of the first collapse on the wall surface, the rebound cavitation bubble will quickly move to the wall surface, and the impact pressure greater than that of the first collapse will act on the wall surface once more.

Shaw et al. [6] used laser-induced cavitation bubbles and a pressure (voltage) transducer to obtain the waveform of the impact of the first collapse of cavitation bubbles on the wall when *γ* is between 0.56 and 1.5. The waveform of the impact of the first collapse of cavitation bubbles on the wall reported in the literature is very similar to that in Figure 4a. The similarities between the two are as follows: (1) there are three peaks in the pressure rise. Except the maximum peak, the other two are distributed on both sides of the maximum pressure value; (2) the time interval between the two peaks before the pressure in the experiment increases to the maximum value is 8 μs, and the time interval during the pressure drop is 3 μs; (3) the time of the pressure rise is slightly greater than that of the pressure drop. The differences between the two are as follows: (1) the time span between several peaks is different. The time span in the experiment in this paper is greater than the corresponding value in the literature; (2) the time span of the pressure rise and that of the pressure drop differ. The time of the pressure rise in this experiment is about 13 μs, and that of the pressure drop about 6 μs, which are both greater than that in the literature. The above phenomenon is mainly caused by the influence of the size of the cavitation bubbles. The maximum radius of the corresponding cavitation bubble is 1.17 mm when *γ* is 0.89 in the literature, and 9.98 mm when *γ* is 0.91 in this paper.

The impact of the jet on the wall will cause a water hammer effect. The water hammer pressure keeps a linear relation with the jet velocity and is calculated as follows [3]:

$$P\_{\rm wh} = \frac{\rho\_{\rm w} \mathfrak{c}\_{\rm w} \rho\_{\rm s} \mathfrak{c}\_{\rm s}}{\rho\_{\rm w} \mathfrak{c}\_{\rm w} + \rho\_{\rm s} \mathfrak{c}\_{\rm s}} \mathfrak{v}\_{\rm l} \tag{1}$$

where *ρ<sup>w</sup>* and *cw* are the density and sound velocity of the jet medium, respectively; *ρ<sup>s</sup>* and *cs* are the density and sound velocity of the solid wall, respectively; and *vt* is the velocity of the jet as it impacts the wall. The wall used in this paper is made of plexiglass, with a density and sound velocity of 2700 kg/m<sup>3</sup> and 2692 m/s, respectively. The density and sound velocity of the jet medium are 1000 kg/m<sup>3</sup> and 1435 m/s, respectively.

Figure 5 is a high-speed photographic image of the collapse of cavitation bubbles in Figure 4b, where *v* is the micro-jet velocity. The frame-rate used in the experiment was 150,000 frames per second, and the pictures in Figure 5 were acquired every 15 frames based on the experimental rate. The micro-jet was formed during the first collapse of the cavitation bubble when the velocity of the micro-jet was approximately 74.96 m/s. When the cavitation bubble entered the rebound stage, a micro-jet was again formed with impact on the wall. The micro-jet developed gradually and eventually imposed an impact on the wall, with its speed decreasing to approximately 13.00 m/s when arriving at the wall. According to the formula for the water hammer pressure, the impact pressure of the micro-jet formed from the rebound cavitation bubbles against the wall is estimated to be 21.36 MPa, and the pressure peak measured by the pressure sensor is 21.89 MPa, both of which are very close. It can be seen that the impact of the micro-jet formed from the rebound cavitation bubbles against the wall mainly comes from the impact of the micro-jet. For the impact of the micro-jet from the first collapse

on the wall when the distance between the bubble and the wall is small, due to the limited speed of the high-speed camera used in this experiment, no clear shock wave can be obtained. Therefore, it is currently impossible to tell whether the main reason for the pressure peak of the impact of the cavitation bubble on the wall from the first collapse is the shock wave or the micro-jet.

**Figure 5.** Process of the impact of micro-jet of rebound cavitation bubbles on a wall (frame-rate: 150,000 fps; exposure time: 5.06 μs; *Rmax* = 11.02 mm, *γ* = 1.79).

### *3.3. Influence of the Distance Between the Cavitation Bubble and Wall on the Impact Strength of the Wall Surface*

In the whole process of cavitation bubble evolution, when the dimensionless distance *γ* meets a certain condition, the second collapse impact strength on the wall surface is greater than the first collapse impact strength. Thus, in this section, the dimensionless distance *γ* and the impact strength on the wall surface are systematically studied to gain the relationship of the peak impact pressure on the wall surface at cavitation bubble collapse stages under different *γ* conditions.

In order to obtain the evolution process and impact process of the cavitation bubble with different dimensionless distances *γ* on the wall surface, the distance from the discharge electrode to the wall surface (with pressure sensor buried inside) will be slightly adjusted in the testing, and the characteristic radius *Rmax* of the cavitation bubble will be changed under each *h* condition. For the peak value of the impact strength of the cavitation bubble on the wall surface, the maximum pressure of the first impact and second impact during the impact of the cavitation bubble on the wall surface is selected in this section to reflect the impact strength of the cavitation bubble on the wall surface.

Figure 6 shows the relationship between the dimensionless distance *γ* and the wall peak pressure; therein, the horizontal axis indicates the dimensionless distance *γ* and the vertical axis indicates the corresponding pressure peaks on the wall surface during the first and second collapses of the cavitation bubble. It can be seen from Figure 6 that for the first collapse of the cavitation bubble, the relationship between the dimensionless distance *γ* and pressure peak on the wall surface shows an exponential-type distribution overall: with the gradual decrease of *γ*, the maximum value of the impact on the wall surface rapidly increases. Then, the peak value will increase first and then rapidly decrease as a whole during the impact of the second collapse of the cavitation bubble on the wall surface. When the distance between the cavitation bubble and the wall increases, the peak pressure on the wall can be divided into the following three parts: when *γ* is less than 1.33, the impact of the first collapse of the cavitation bubble on the wall surface is greater than the impact of the second collapse on the wall surface; when *γ* between the cavitation bubble and wall is greater than 1.33 and less than 2.37, the impact strength on the wall surface caused by the second collapse of the cavitation bubble will be greater than the impact strength of the first collapse on the wall surface; when *γ* is greater than 2.37, the impact of the first collapse of the cavitation bubble on the wall surface is greater than the impact of the second collapse on the wall surface, and the maximum of the first impact peak is less than 4 MPa.

**Figure 6.** The relationship between the dimensionless distance (γ) and the wall peak pressure.

It can be seen from the above analysis that for the collapse of cavitation bubbles near the wall surface, when the dimensionless distance *γ* is less than 2.37, the impact of the first and second collapses of the cavitation bubble on the wall surface is huge; when *γ* is greater than 1.33 and less than 2.37, the impact on the wall surface caused due to the second collapse of the cavitation bubble is greater than the impact on the wall surface caused due to the first collapse. The reason why the above impact characteristics are seen can be analyzed from the following two aspects: on the one hand, when the bubble is very close to the wall surface, during the first bubble-asymmetric shrinkage, a micro-jet with a direct impact on the wall is formed. When the cavitation bubble is a little further away from the wall surface, the impact of the wall on the first shrinkage form of the cavitation bubble is weakened, and a micro-jet is formed during the shrinkage of the cavitation bubble. The micro-jet towards the wall surface drives the cavitation bubble, and finally, during the second shrinkage, the micro-jet is formed. This micro-jet can impose an impact on the wall. Although the energy of the cavitation bubble itself is attenuated after the first cycle of evolution, the impact of the micro-jet from the second shrinkage on the wall is still considerably large. As the distance between the cavitation bubble and the wall surface further increases, the bubble during the second shrinkage obviously moves toward the wall surface, but the resulting micro-jet from the second collapse is not sufficient to impose an impact on the wall, and then the bubble enters the next expansion and contraction stage. On the other hand, the cavitation bubble has a certain energy after the birth. The energy is dissipated due to the viscosity of the water body and other factors after the first and second expansion. After the two-fold expansion and shrinkage, the energy is significantly reduced. As a result, less obvious movement toward the wall is observed and the impact on the wall is reduced greatly in the subsequent evolution of the cavitation bubble.

Philipp and Lauterborn [3] studied the deterioration of the wall under different bubble–wall conditions when the cavitation bubble's *Rmax* = 1.45 mm with the laser-induced cavitation technology. There were over 100 groups of samples. By analyzing the deterioration and bubble–wall distance *γ*, it is found that when *γ* ≥ 2.2, no deterioration was found on the solid wall; when 2.2 ≥ *γ* ≥ 1.5, the area of deterioration increased with the decrease of *γ*. By comparing the deterioration in the literature and the impact of the cavitation bubble on the wall obtained in this paper, it can be found that the critical values are basically consistent. Since a low-pressure discharge-induced cavitation bubble is used here, it is more difficult to achieve cavitation under the condition of smaller values of *γ*. Furthermore, the cavitation bubble was larger with the low-voltage underwater discharge technology used in this paper, and the size in the literature is smaller, which may cause the slight difference in the

critical value between the paper and the literature. This difference, however, does not affect the zoning characteristics of the impact of the cavitation bubble on the wall.

### **4. Conclusions**

Through the combination of the high-speed dynamic collection and analysis system and transient pressure testing system, the impact characteristics of the collapse of the low-voltage spark-discharge cavitation bubble on the wall surface are studied in this paper, drawing the following important conclusions:

(1) The dimensionless distance between the cavitation bubble and wall (*γ*) has an important influence on the formation of the micro-jet in the late evolution stage of the cavitation bubble. When the dimensionless distance *γ* is less than 1.33, a micro-jet directly impacting the wall surface can be formed during the first contraction of the cavitation bubble; subsequently, the collapse of the main impact area will occur. When *γ* is within 1.33~2.37, the second collapse of the cavitation bubble after rebound will cause a greater impact on the wall surface, then the collapse of the secondary impact area will occur. When *γ* is greater than 2.37, the influence of wall on the development of a micro-jet is not obvious, that is, collapse in the slow release area occurs.

(2) Concerning the peak value of the impact pressure of the cavitation bubble collapse on the wall surface, the impact of cavitation bubbles on the wall surface during the main impact gradually decreases with the increase of the distance between the cavitation bubble and the wall, and the impact of the cavitation bubble on the wall surface during the secondary impact increases first and then decreases. More specifically, when the dimensionless distance *γ* is less than 1.33, the impact on the wall surface due to cavitation bubble collapse is mainly from the first collapse, and with increasing *γ*, the impact of the second collapse on the wall surface increases. When *γ* is within 1.33~2.37, the impact of the cavitation bubble on the wall surface during the first collapse will be less than the impact strength of the re-collapse after rebound on the wall surface with increase of *γ*. When *γ* is greater than 2.37, the impact of the first and second collapses of the cavitation bubble on the wall surface will decrease considerably.

The above conclusions have potential theoretical value for the anti-erosion design of hydraulic engineering, the erosion resistance of blades in hydraulic machinery and the better exertion of the cavitation effect in ultrasonic cleaning technology.

**Author Contributions:** Conceptualization, W.X. and J.L.; Methodology, J.D. and J.L.; Software, J.L.; Validation, W.X. and J.D.; Formal Analysis, J.L.; Investigation, J.L.; Resources, J.L.; Data Curation, Y.Z. and Q.Z.; Writing-Original Draft Preparation, J.L.; Writing-Review & Editing, J.D.; Visualization, J.L.; Supervision, W.X.; Project Administration, W.X.; Funding Acquisition, J.D.

**Funding:** This research was funded by the National Key R&D Program of China (Grant No.2016YFC0401901) and the National Natural Science Foundation of China (Grant No.51409180).

**Conflicts of Interest:** The authors declare no conflict of interest.

### **References**


© 2018 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

### *Article* **Numerical and Experimental Comparative Study on the Flow-Induced Vibration of a Plane Gate**

**Chunying Shen 1,2, Wei Wang 1,\*, Shihua He <sup>2</sup> and Yimin Xu <sup>2</sup>**


Received: 13 September 2018; Accepted: 29 October 2018; Published: 31 October 2018

**Abstract:** A numerical method is applied here to simulate the unstable flow and the vibration of a plane gate. A combination of the large eddy simulation (LES) method and the volume of fluid (VOF) model is used to predict the three-dimensional flow field in the vicinity of a plane gate with submerged discharge. The water surface profile, the streamline diagrams, the distribution of turbulent kinetic energy, the power spectrum density curve of the fluctuating pressure coefficient at typical points underneath the gate, and the complete vortex distribution around the gate are obtained by LES-VOF numerical calculation. The vibration parameters of the gate are calculated by the fluid-structure coupling interface transferring the hydrodynamic load. A simultaneous sampling experiment is performed to verify the validity of the algorithm. The calculated results are then compared with experimental data. The difference between the two is acceptable and the conclusions are consistent. In addition, the influence of the vortex in the slot on the flow field and the vibration of the gate are investigated. It is feasible to replace the experiment with the fluid-structure coupling computational method, which is useful for studying the flow-induced vibration mechanism of plane gates.

**Keywords:** LES-VOF method; flow-induced vibration; plane gate; numerical analysis

### **1. Introduction**

With the wide application of a high head gate in hydraulic engineering, the vibration problem of the plane gate is increasingly prominent. Essentially, the flow-induced vibration of the plane gate is a complex fluid-structure interaction phenomenon. The time varying hydrodynamic load that is caused by the unstable flow around the gate is the main excitation source. The coexistence of the vortices at the gate's bottom edge, in the gate slot and in the gate downstream, further complicates the flow structure of the sluice flow. However, it is difficult to capture these vortices simultaneously and to catch the comprehensive information of the vortex-induced vibration by using an experimental approach. Researchers, such as Hardwick [1], Kolkman [2], Jongeling [3], and Ishii [4], have been making efforts to explore the mechanism of gate vibration. Hardwick conducted a model test investigation on the plane gate vibration. It is believed that this nonlinear vertical vibration response is caused by the resonance between the fluid shear force acting at the bottom edge of the gate and the gate. Kolkman et al. believed that the flow fluctuation and the variation of the pressure due to the inertia effect of sluice flow might be the main reason for the vertical vibration of the gate. Jongeling proved that the streamwise vibration of the gate is induced by the instability of the bottom edge shear layer. Ishii considered that the vortices at the bottom edge of the gate and the interaction with the structure lead to the streamwise vibration of the gate. These research results are mostly limited to the qualitative description of the physical process of the gate vibration, and tend to separate the vertical vibration and the streamwise flow vibration. Therefore, it is especially important to analyze the vibration of the gate combining the vertical and the streamwise vibration. Simultaneously, these researches mainly focus on investigating the vibration mechanism that resulted from the vortex at the bottom edge of the gate.

With the rapid development of computational fluid dynamics (CFD) and computer technology, the accuracy of numerical calculation of flow field and the complexity of solving problems have been greatly improved, and it has been possible to numerically simulate the sluice flow in some extent. The method adopted is mainly based on the finite volume scheme. The SIMPLE series algorithm is used in the calculation of pressure field and velocity field. When dealing with turbulence problems, the turbulence model is employed to solve the Reynolds time-averaged Navier-Stokes equation. Pani et al. [5] used a finite element technique to investigate the fluid-structure interaction effect on the hydrodynamic pressure. Erdbrink et al. [6] also assessed the local velocities and pressures of the two-dimensional field by using the finite element method. Liu et al. [7] adopted the two-dimensional (2-D) mathematical model of Reynolds-averaged equations to study the turbulent flow behind a sluice gate. Kostecki [8–10] predicted the two-dimensional flow field in the vicinity of an underflow vertical lift gate by using a combinative numerical model of the vortex method and the boundary element method. Kazemzadeh [11] applied the smoothed fixed grid finite element method to simulate free surface flow in gated tunnels. In recent years, some novel computational methods and models have been developed. For example, Zhu et al. [12,13] proposed a new slip model, which was used in the study of lid-driven cavity flows in both continuum and transition flow regimes. However, there is still a lack of theories and methods to analyze effectively and deal with the complex vibration phenomenon of plane gates.

The development of parallel algorithms and supercomputers has strongly promoted the progress of direct numerical simulation (DNS) and large eddy simulation (LES) methods. Analyzing the pressure pulsation of turbulent flow using LES method has received more and more attention. However, the application of CFD technology, especially simulating pressure pulsation, in the sluice flow is still preliminary.

Different from the arc gate, the side walls of a plane gate have the gate slots that are used for the lifting movement of gate. In the case of the submerged sluice flow, the random movement of multiple vortices behind the gate, in the gate slot and at the gate's bottom edge will aggravate the pulsation of velocity and the pressure. When water flow crosses the slot, such as that which exists in the hydraulic plane gate, the vortices in the slot will interact with the separated flow near the downstream corner of the slot due to local boundary mutation in the slot area, which can easily lead to flow cavitation and even endanger the safe operation of hydraulic projects [14]. The vortex-induced vibration of a blunt body or other solid boundary caused by vortex shedding has been studied most [15–18], while studies of the gate vibration induced by the vortex in the gate slot has been reported rarely.

Based on the theory of fluid-structure interaction, the present investigation studies the characteristics of the flow-induced vibration of the hydraulic plane gate with a submerged discharge. The three-dimensional numerical simulation method with a large eddy simulation turbulence model can obtain a complete flow structure, including a funnel vortex in the slot. The vibration parameters of the gate are acquired by the fluid-structure coupling interface transferring the hydrodynamic load. The fluid-structure coupling method is validated by the experimental evidence.

### **2. Computational Modeling**

### *2.1. Turbulence Model*

The large eddy simulation (LES) method uses instantaneous flow control equations to directly simulate the large scale vortex in a turbulent flow field. After processing by using the filtering function, the instantaneous Navier-Stokes equation and continuity equation of control fluid flow become:

$$\frac{\partial}{\partial t}(\rho \overline{\mathbf{u}}\_i) + \frac{\partial}{\partial \mathbf{x}\_j}(\rho \overline{\mathbf{u}}\_i \overline{\mathbf{u}}\_j) = -\frac{\partial \overline{p}}{\partial \mathbf{x}\_i} + \rho g\_i + \frac{\partial}{\partial \mathbf{x}\_j} \left[ \mu \left( \frac{\partial \overline{u}\_i}{\partial \mathbf{x}\_j} + \frac{\partial \overline{u}\_j}{\partial \mathbf{x}\_i} \right) - \rho \mathbf{r}\_{ij} \right] \tag{1}$$

$$\frac{\partial \rho}{\partial t} + \frac{\partial}{\partial x\_i} (\rho \overline{u}\_i) = 0 \tag{2}$$

where *i*, *j* = 1, 2, 3; *ρ* is the specific mass; *μ* is the dynamic viscosity; *gi* is the mass force; and, *p* and *ui* are the filtered pressure and velocity components, respectively. *<sup>τ</sup>ij* <sup>=</sup> \_\_\_\_\_ *uiuj* − *uiuj* are the components of the sub grid scale (SGS) stress tensor representing the effect of the small scale motion on the large scale motion. According to the basic SGS model, Smagorinsky [19] assumes the SGS stress in the following form:

$$
\pi\_{\rm ij} - \frac{1}{3} \pi\_{kk} \delta\_{\rm ij} = -2\mu\_t \overline{R}\_{\rm ij} \tag{3}
$$

The Smagorinsky-Lilly model is used in this paper. In the Smagorinsky-Lilly [20] model, the eddy-viscosity is modeled by:

$$
\mu\_t = (\mathbb{C}\_s \Delta) \left| \overline{\mathbb{R}} \right| \tag{4}
$$

where *R* = - 2*RijRij*, *Rij* = <sup>1</sup> 2 *∂ui ∂xj* + *∂uj ∂xi* , *Cs* is the Smagorinsky constant; and, Δ is the local grid scale, which is computed according to the volume of the computational cell using Δ = *V*1/3.

It is more difficult to simulate the free surface accurately, and the volume of fluid (VOF) method is more advantageous to track the complex free surface in the multiphase flow models. Due to the downstream flow directly in contact with the atmosphere, the VOF model is used to track the air-water surface in the present investigation.

The tracking of the interface between the phases is accomplished by solving the continuity equation of the volume fraction of one or more phases. The volume fraction of water is a function of time and space, so the VOF model must be solved by transient. Assuming that *α<sup>w</sup>* denotes the volume fraction for water phase, *α<sup>w</sup>* = 0 indicates the exclusion of the water phase in the region, *α<sup>w</sup>* = 1 indicates that the region is full of water phase, and 0 < *α<sup>w</sup>* < 1 indicates the interface for water phase and other phases.

The tracking equations of the interface are as follows:

$$\frac{\partial \alpha\_w}{\partial t} + u\_i \frac{\partial \alpha\_w}{\partial x\_i} = 0 \tag{5}$$

### *2.2. Structural Model*

The gate structure is assumed to consist of linear elastic small deformation material, which is made of organic glass. Under the Lagrangian system, structural vibration control equations with a small deformation assumption are as follows:

$$\rho^s \frac{\partial^2 \mathbf{u}^s}{\partial t^2} = \text{div}\left(J\sigma^s F^{-\mathbf{T}}\right) + \rho^s \mathbf{g} \tag{6}$$

$$
\sigma^s = f^{-1} F(\lambda^s(\text{tr}E)I + 2\mu^s E)F^T \tag{7}
$$

The *F*, *E*, and *J* are defined as:

$$F = I + \nabla \boldsymbol{u}^s \tag{8}$$

$$E = \frac{1}{2} \left( F^T F - I \right) \tag{9}$$

$$J = \det F \tag{10}$$

where *ρ<sup>s</sup>* is the density of structure, *u*<sup>s</sup> is the vibration displacement, *σ*<sup>s</sup> is Cauchy's stress tensor of the gate structure, *<sup>λ</sup><sup>s</sup>* = *<sup>υ</sup>Es*/((<sup>1</sup> + *<sup>υ</sup>*)(<sup>1</sup> − <sup>2</sup>*υ*)) and *<sup>μ</sup><sup>s</sup>* = *<sup>E</sup>s*/(2(<sup>1</sup> + *<sup>υ</sup>*)) are the constants of the gate material, *E<sup>s</sup>* is material elastic modulus, and *υ* is Poisson's ratio.

For complete fluid and structure coupling, the geometry coordination condition and the force equilibrium condition on the fluid-structure coupling interface should be satisfied.

The geometry coordination condition on the fluid-structure coupling interface *<sup>S</sup>f s* = <sup>Ω</sup>*<sup>f</sup>* ∩ <sup>Ω</sup>*<sup>s</sup>* is:

$$
\mu\_i^f = \mu\_i^s \text{ on } S^{fs} \tag{11}
$$

The force equilibrium condition on the fluid-structure coupling surface *<sup>S</sup>f s* = <sup>Ω</sup>*<sup>f</sup>* ∩ <sup>Ω</sup>*<sup>s</sup>* is:

$$
\sigma\_{\rm ij}^f n\_{\rm j}^{fs} = \sigma\_{\rm ij}^s n\_{\rm j}^{fs} \quad \text{on} \quad S^{fs} \tag{12}
$$

where *nf s <sup>j</sup>* is the normal unit vector on the fluid-structure interface, <sup>Ω</sup>*<sup>f</sup>* is the fluid region, and <sup>Ω</sup>*<sup>s</sup>* is the structure region. Due to the water flow effect, the gate structure contacting with the fluid vibrates, which leads to flow field changes. The changing flow field acts on the gate structure again, which forms fluid-structure interaction problems.

### *2.3. Meshing and Numerical Approach*

The direction of flow is the *x* direction, the depth direction of flow is the *y* direction, and the direction across the flume is the *z* direction. Submerged discharge occurs downstream of the gate. The flow inlet is the speed boundary, and the flow outlet is the pressure boundary converted by the head. The top of the downstream flow is the air boundary, and the relative pressure is zero. The flume bottom and the sidewalls are no-slip solid boundary conditions, and the given normal velocity is zero. The method of the wall function is used for processing the near-wall viscous sublayer.

Due to the symmetry of the gate and the flume, half the model is selected to be meshed to reduce the calculating work. In this paper, the cases of the calculation are shown in Table 1, where *e* is the opening of the gate and *H* is the upstream water head. The underflow ratio *U*, is used to access the characteristics of the gate vibration caused by the submerged flow. Comprehensively considering the upstream and the downstream water head, underflow ratio *U* is defined, as follows:

$$
\Delta U = \frac{h\_t - h\_c''}{H\_0 - h\_c''} \tag{13}
$$

where *ht* is the downstream water depth, *hc* is the contraction section depth, *h <sup>c</sup>* is the conjugate water depth of the contraction section, and *H*<sup>0</sup> is the total head.

**Table 1.** A list of the cases of calculation.


For *U* = 0.059, the half grid of fluid with tetrahedral mesh has 252,758 units and a total of 49,112 nodes (Figure 1). The time step length of unsteady flow is 0.001 s, and the maximum number of iterations of each step is 20 times. At the same time, the numerical simulation of turbulent flow under the unslotted condition is carried out to calculate the vortex effect of the gate slot. The half grid of the plane gate is shown in Figure 2a, and the size of the plane gate is shown in Figure 2b.

The ANSYS Worckbench17.0 is used to solve the system of dependent variables. After constructing a three-dimensional fluid-structure coupling model, the flow field in the vicinity of the lift gate with a submerged discharge condition is simulated by the finite volume method. A combination of the k − double-equation turbulence model and the VOF model of the multiphase flow is used to calculate the flow field before 10.5 s. With the result as the initial conditions of the Smagorinsky-Lilly subgrid model and the LES method after 10.5 s, the VOF model is still adopted to determine of the water-air interface. Velocity pressure coupling by the PISO algorithm and the momentum equation in discrete form by the second-order windward scheme are used to obtain the hydrodynamic pressure distribution that acts on the gate surface. Under the action of the surface pressure load, the vibration response of the gate structure is analyzed by the finite element method. To analyze the calculation results visually, the magnitude and the materials of the calculation model are consistent with the synchronous experiment model. The density, elastic modulus, and the Poisson's ratio of the gate material are 1400 kg/m3, 3 × 109 N/m2, and 0.4, respectively. Constraints are imposed on the top and sides of the gate. The width-to-depth ratio of the gate slot is 1.78.

**Figure 1.** The grid of fluid.

**Figure 2.** The plane gate: (**a**) the grid; and, (**b**) the size (unit: mm).

### **3. The Results of Numerical Simulation**

To verify the validity of the numerical simulation method, the numerical results are compared with the experimental data. The flow field behind the plane gate with submerged discharge, the hydrodynamic pressure process at typical points underneath the gate and the vibration parameters of the gate are acquired by the combination experiment [21] of three-dimensional particle image velocimetry (3D-PIV) made by TSI Company, the digital pressure sensor with high precision and the multichannel vibrating data acquisition system made by Lance Measurement Technologies Company, respectively.

The PIV test uses the Insight3G software to start the laser, and simultaneously captures the flow field particle images with two PIV-specific cross-frame CCD cameras. The autocorrelation or cross-correlation principle is used to extract the image of flow characteristics. Finally vector diagram of transient flow velocity in the measured range is obtained by professional post-processing software, such as Tecplot and Matlab. The maximum emission frequency of the laser is set to 14.5 Hz. In this test, a certain concentration of SiO2 is selected as the tracer particle, and the particle size is 10–15 μm. The three-axis accelerometer is directly attached to the upstream surface of the gate. The accuracy of the pressure sensor and accelerometer both are ±0.1%FS, and the sampling interval of the pulsating pressure and the vibration displacement are 0.01 s and 0.001 s, respectively. Three instruments are tested simultaneously.

The size of the organic glass flume is 4880 mm × 100 mm × 190 mm, and the length behind the plane gate is approximately 3000 mm. The calculations and the experimental cases are consistent with each other. The general arrangement of the experiment is shown in Figure 3.

**Figure 3.** General arrangement of the experiment.

### *3.1. Flow Field Distribution*

The calculation length after the lift gate is 920 mm and width B of the flume is 100 mm. The computation results of the time-averaged surface profiles with the whole flow passage of the modeling after the gate in *U* = 0.059 are shown in Figure 4. Blue represents the water phase and red represents the gas phase.

**Figure 4.** The color nephogram of the time-averaged flow profile.

Figure 5 shows the time-averaged flow profile comparison behind the gate of section *z* = 0 mm for the calculation and experiment of *U* = 0.059. The maximum and minimum of relative difference of the height between the calculation and the test values of the flow profile are 4.8% and 0.16%, respectively, which shows that the VOF method can simulate the water-air interface well.

**Figure 5.** The time-averaged flow profile.

The calculational time-averaged streamlines on sections *z* = 0 mm and *z* = 35 mm in different cases are presented in Figure 6. According to (a), (b), and (c) in Figure 6, there is only a vortex on section *z* = 35 mm in the time-averaged flow field behind the plane gate in the *x* direction, but there are two vortices on section *z* = 0 mm that are within the 300 mm range. At different sections in the same case, the vertical distances of the vortex center from the *x*-axis are approximately the same. With the *e/H* increasing, the vertical distance of the vortex center from the *x*-axis increases gradually.

**Figure 6.** The calculational time-averaged streamlines.

The size of the PIV calibration target is 200 mm × 200 mm. Due to the size limitation, it is impossible to test the vortices in the slot and at the bottom edge of the gate by PIV, so the computational results of only the vortices behind the gate are compared carefully with those of the three-dimensional PIV test. Both the simulated and the experimental velocity vector diagrams are achieved at 20.6 s.

Figure 7 shows the instantaneous streamlines behind the gate on the same section *z* = 0 mm for the calculation and the experiment. The numerical simulation can make up for the drawback of losing the velocity vector by 3D-PIV positioning. When compared with the experimental results, the vortex center, which is closest to the gate obtained by the simulation, is lower in the *y* direction under *U* = 0.059, but this distance is basically same for the other two cases. On the whole, the numerical simulation still obtained the most ideal vortex information.

**Figure 7.** The instantaneous streamlines (*z* = 0 mm).

It can be seen from the instantaneous streamlines in Figure 7a–c that there are the same number of vortexes on section *z* = 0 mm between the computation and the experiment. Under *U* = 0.059, the numerical simulation captures the phenomenon of the vortex detachment from the gate's bottom edge, but there is no such phenomenon when *U* = 0.355 or *U* = 0.730. Therefore, when the *e/H* and *U* reach a certain extent, the vortex shedding phenomenon can be generated from the bottom edge of the gate.

Normalizing the turbulent kinetic energy by the square of the mean flow velocity of the gate downstream section, Figure 8 presents the computational normalized turbulent kinetic energy distribution on section *z* = 0 mm and section *z* = 35 mm. Turbulent kinetic energy under the gate on section z = 0 mm is higher than the value on section *z* = 35 mm. The different magnitude of turbulent kinetic energy is one of the significant factors affecting the gate vibration extremum. The experiment proved that the extremum of vibration displacement on section *z* = 0 mm is larger than that on section *z* = 35 mm.

**Figure 8.** Normalized turbulent kinetic energy distribution.

### *3.2. Hydrodynamic Pressure*

The measuring points of hydrodynamic pressure just below the gate are point 1 (0, 0, 0 mm) and point 2 (0, 0, 35 mm). The digital pressure sensor with high precision is used to test the hydrodynamic pressure with a sampling frequency of 0.01 s, while the monitoring points are also set up in the numerical calculation.

Figure 9 shows the time history curves of the fluctuating pressure coefficient of points 1 and 2 in 10 s. *Cp* = (*p* − *p*)/0.5*ρv2*, where *Cp* is the fluctuating pressure coefficient, p is the instantaneous pressure, *p* is the average pressure, *ρ* is the water flow density, and *v* is the average velocity of the exit section. Under the same case, the fluctuating pressure of point 1 is stronger than that of point 2. The time to reach the extreme of the two points is out of synchronization.

Figure 10 shows the normalized power spectrum density curve of the pressure for point 1 in different cases. *St* is the Strouhal number and *St* = f · L/*v*, Where *f* is the frequency of the pressure pulsating, *L* is the characteristic length and taken as the gate opening, and *v* is the mean velocity of section just below the gate. With the *e/H* or *U* increasing, the energy of pulsating pressure of point 1 decreases. The computational Strouhal numbers for dominant frequencies of the pressure are 0.332 × <sup>10</sup>−3, 2.217 × <sup>10</sup>−3, 1.545 × <sup>10</sup>−<sup>3</sup> with the corresponding cases of *<sup>U</sup>* = 0.059, 0.355, 0.730. Accordingly, the experimental Strouhal numbers for dominant frequencies of the pressure are 0.337 × <sup>10</sup><sup>−</sup>3, 2.242 × <sup>10</sup><sup>−</sup>3, 1.552 × <sup>10</sup>−<sup>3</sup> in corresponding cases. The Strouhal number difference for main frequency of the pressure between the computational and the experimental results is very small.

**Figure 9.** The time history curve of the pressure pulsating coefficient for different point.

**Figure 10.** The power spectrum density (PSD) curve of the pressure for point 1 in different cases.

In order to illustrate the influence of the gate slot on the flow field and the gate vibration, both calculation and experiment are carried out on the filling of the slots underneath the gate (unslotted). Figure 11 shows the time history curve of the slotted and unslotted pressure pulsating coefficients in *U* = 0.059. It can be found that the maximum values of the slotted pressure pulsation coefficients for both points 1 and 2 were larger than those of the without a slot. The characteristic values of the pressure pulsating coefficients (slotted and unslotted) for both the experiment and for the calculation in *U* = 0.059 are presented in Figure 12. When comparing the value of the computational and the synchronous experimental results, the relative difference is within 10%, which indicates that the numerical simulation results meet the basic requirements.

**Figure 11.** The time history curve of the pressure pulsating coefficients (slotted and unslotted).

**Figure 12.** The comparison of characteristic value of the pressure pulsating coefficient.

In the experiment, the dominant frequency of the pressure at the measuring points underneath the gate are only 0.025~1.245 Hz, and the dominant frequency of the pressure at the measuring points behind the gate are only 0.183~14.63 Hz. The natural frequencies of the gate vibration by the modal analysis are 247.43 Hz, 529.85 Hz, 1086.6 Hz, 1280.6 Hz, 1605.1 Hz, and 1963.1 Hz, respectively. The natural frequency of the gate vibration is far from the dominant frequency of the flow pressure, so it would not have produced resonance.

Figure 13a presents the computational pressure distribution and streamlines in the gate slot for different level cross sections (*y*/*e* = 0.067, *y*/*e* = 0.2, *y*/*e* = 0.333, *y*/*e* = 0.467, *y*/*e* = 0.6, *y*/*e* = 0.733) from the bottom to the top in *y* direction for *U* = 0.059. It can be seen from the superimposed vortex diagram that the low-pressure zone of the vortex gradually enlarges from the bottom to the top, which is shaped like a funnel and it shifted to the downstream. The pressure is lowest in the vortex center where it is also most prone to generate cavitation in the gate slot. Due to the water impacting the slot from the left to the right, the pressure is highest at the right of the slot. The closer to the bottom the place is, the higher the pressure is. Under *U* = 0.355 (Figure 13b) and *U* = 0.730 (Figure 13c), the range of the low pressure and the high pressure zones in the gate slot is much smaller than that in *U* = 0.059. The funnel vortex of *U* = 0.355 still exists, but the center of the vortex is shifted to the upstream from the bottom to the top, and the center pressure of the bottom vortex is relatively large. For *U* = 0.730, the centerline of the vortex is similar to the spiral line, and the funnel vortex is not obvious.

**Figure 13.** Pressure (pa) distribution and streamlines at 22.0 s in the gate slot.

The maximum vibration displacement of the gate with a submerged condition in the *x* direction is not due to the resonance, but it is due to the mixed vortex-induced process, including the vortices behind the gate, in the gate slot and at the gate's bottom edge. It can be seen from the above analysis that the numerical simulation method in this paper can obtain not only the pressure fluctuation distribution, but also a more complete flow structure, such as multiple vortices of the sluice flow, and it forms the excitation force of the plane gate vibration.

### *3.3. Vibration of the Gate*

Define the non-dimensional coefficient *Kd* = (*d* − *d*)/*σd*, where *d* is vibration displacement, *d* is the average vibration displacement, and *σ<sup>d</sup>* is the vibration standard deviation. Point A, which is located on the central axis of the upstream surface of the gate and 10 mm away from the bottom edge, is used to study the response of the gate vibration.

The normalized power spectrum density curve of vibration displacement for point A in *U* = 0.059, which is obtained by the computation and the experiment, is shown in Figure 14. In the experiment, the acceleration sensor and the multichannel data acquisition instrument were used to monitor the acceleration of the gate vibration in the *x* direction and *y* direction. The vibration displacement is achieved by the second integral of the acceleration. The computational Strouhal number for dominant frequency of vibration displacement in *x* direction is 2.582 with the corresponding experiment of 2.660. In addition, the computational Strouhal number for dominant frequency of vibration displacement in *y* direction is 3.724 with the corresponding experiment of 3.768. The Strouhal number for main frequency of the vibration displacement between the computational and the experimental results is closer.

**Figure 14.** The power spectrum density (PSD) curve of vibration displacement for point A.

Figure 15 plots the time-history curve of the vibration displacement coefficient *Kd* in the *x* direction and in the *y* direction of point A, with and without slots, by using the numerical simulation. It can be seen that the maximum value of the vibration displacement coefficient of the plane gate with a slot is larger than that of the plane gate without a slot in both the *x* and *y* directions. When there is no slot, the maximum value of the *x* direction vibration displacement coefficient appears at 12.7 s and the minimum value appears at 21.8 s. The maximum value in the *y* direction appears at 14.7 s and the minimum value appears at 17.6 s. The moment when the vibration displacement coefficient of the *x* direction and the *y* direction appears is out of sync.

**Figure 15.** The time history curve of the vibration displacement coefficient *Kd* of point A.

### **4. Conclusions**

In this study, numerical simulations were performed to compute the flow field and the water surface profile for a plane gate with a submerged discharge by the combination of the LES and the VOF methods. Overall, the numbers and locations of the vortex are consistent with the synchronous experimental results. When compared with the experimental values, the relative difference of both the standard deviation pressure and the maximum pressure of the calculation is within 10%, which verifies the validity of the numerical simulation methods that were adopted in the present work. Simultaneously, the maximum vibration displacement coefficient was acquired from a calculation of the vibration displacement process of Point A, which was very close to that of the experiment. So, the fluid-structure coupling computational method used in certain conditions can replace the complex experiment and make up for the lack of experiment.

The turbulent kinetic energy under the gate on section *z* = 0 mm was higher than the value on section *z* = 35 mm. The change trend of the turbulent kinetic energy was consistent with the amplitude of the gate vibration displacement, which shows that the turbulent kinetic energy under the gate plays an important role in the amplitude of the gate vibration. The maximum value of the vibration displacement coefficient of the plane gate with a slot was larger than that of the plane gate without a slot in both *x* and *y* directions. Thus, the vortex in the gate slot aggravates the vibration of the gate. Turbulence pressure pulsation and multiple vortices in the flow field are the main causes of gate vibration. The mixed vortex-induced process includes the vortices behind the gate, in the gate slot, and at the gate's bottom edge.

**Author Contributions:** C.S. and W.W. conceived and designed the model; C.S. and S.H. set up and debugged the programmes; C.S. analysed the results and wrote the paper; C.S. and Y.X. provided editorial improvements to the paper.

**Acknowledgments:** Financial support provided by the National Key Research and Development Program (Grant No. 2016YFC0401603) and the National Natural Science Foundation of China (Grant No. 51369013) are gratefully acknowledged.

**Conflicts of Interest:** The authors declare no conflict of interest.

### **References**


© 2018 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

### *Article* **Experimental Optimization of Gate-Opening Modes to Minimize Near-Field Vibrations in Hydropower Stations**

### **Yong Peng 1, Jianmin Zhang 1,\*, Weilin Xu <sup>1</sup> and Matteo Rubinato <sup>2</sup>**


Received: 26 July 2018; Accepted: 10 October 2018; Published: 12 October 2018

**Abstract:** Multi-Horizontal-Submerged Jets are successfully applied to dissipate energy within a large-scale hydropower station. However, notable near-field vibrations are generated when releasing high discharges through the gates, which is generally typical in a flooding case scenario. Under these conditions, the magnitude of the vibrations varies when applying different gate-opening modes. To investigate and find optimized gate-opening modes to reduce the near-field vibration, multiple combinations were tested by varying gate-opening modes and hydraulic conditions. For each of the tests conducted, fluctuating pressures acting on side-walls and bottoms of a stilling basin were measured. The collected datasets were used to determine the maximum and minimum fluctuating pressure values associated with the correspondent gate-opening mode and a detailed comparison between each of the gate-opening modes was completed. The paper presents the quantitative analysis of the discharge ratio's effect on fluctuating pressures. It also investigates the influence of different gate-opening modes by including side to middle spillways and upper to lower spillways configurations. The flow pattern evolutions triggered by each different gate-opening mode are discussed and optimal configurations that minimize near-field vibrations at high discharges are recommended to support both the design of new systems and assessment of the performance of existing ones.

**Keywords:** multi-horizontal-submerged jets; energy dissipation; near-field vibration; fluctuating pressure; gate-opening modes

### **1. Introduction**

Since 2000, over 300 hydropower projects have been completed or are under construction in China [1]. Flood discharge and energy dissipation are crucial to maintaining the safety of these hydropower installations. Furthermore, lessons must be drawn from the past dam disasters, e.g., Vajont Dam [2,3] and Banqiao Dam [4]. Decades ago, hydropower stations were not constructed in proximity to residential areas, and as a result, research at that time mainly focused on the stability and the functionality of the hydraulic structures. However, in recent years, due to increasing urbanization [5], some large hydropower stations had to be located in the proximity of residential areas. Moreover, due to the climatic scenario [6], it is likely that in the future there will be numerous and more frequent intense rainfall events, and consequently large volumes of water may need to be released from dam reservoirs. Therefore, it is essential to provide solutions to a need for increasing the stability of these hydraulic structures and surrounding areas, minimizing environmental impacts due to extreme flow discharges [7].

The purpose of this work is to optimize a recently developed novel technique for energy dissipation, Multi-Horizontal-Submerged Jets (simplified as MHSJ). MHSJ have been studied and improved by many researchers in recent years [8–17]. It has been demonstrated that MSHJ possesses three main advantages when compared with the hydraulic jump, which is the traditional energy dissipater: (1) lower atomization; (2) a higher rate of energy dissipation; (3) more flexibility of gate-opening modes [12–14]. Thus, it seems feasible to implement MHSJ to tackle the problem of energy dissipation within large hydropower stations that are typically characterized by higher water heads and larger flowrates. The adoption of MHSJ was tested by Deng et al. [12] and successfully applied as energy dissipation in large-scale hydropower stations [14].

However, despite positive progress on solving energy dissipation and erosion problems, there is still a lack of studies on the considerable flow-induced vibrations that occur in this kind of structure and propagate towards the near-field when releasing the greater flowrates. These near-field vibrations can be considered a serious threat for the stability of the structure, and hence this phenomenon needs to be further investigated. Yin and Zhang [18] and Yin et al. [19] studied the relationship between the discharge flow and the induced vibrations and simulated the results with a finite element method. Yin et al. [19] found that at higher flowrates during flooding situations, the flow is complex and unsteady, and this is the main cause of near-field vibrations. Therefore, if the vibrations are caused by releasing higher flowrates, optimizing gate-opening modes could reduce this effect and near-field vibrations could decrease correspondingly. To date, this issue has not been investigated in detail and to address this gap, this paper presents an experimental study to explore the relationship between complex flow scenarios and gate-opening modes in order to recommend an optimal solution to minimize the formation of vibrations.

### **2. Methods**

The experimental facility, shown in Figure 1, was constructed based on a large-scale hydropower station (Xiangjiaba hydropower station. It sits on the Jinsha River, a tributary of the Yangtze River in Yunnan Province and Sichuan Province, Southwest China). The experimental facility follows Froude similarity, i.e., the ratio between inertial and gravitational forces remains constant, with the scale of 1:80. The experimental model is composed of the upper and lower reaches of the river, a dam, mid-level channels, spillways and two stilling basins with identical layout (the difference between two stilling basins is the their location on the dam, being parallels). The single stilling basin includes 5 mid-level channels and 6 spillways. Both of mid-level channels and spillways were arranged alternately and division piers were located between them.

In the model, the outlets of spillways are horizontal with symmetrical contraction of 1 m and width of 8 m. Also, the mid-level channels' outlets are horizontal, without any contraction, and their width is 6 m. The division piers are 3 m wide. The offset of mid-level channels, spillways, and division piers from the bottom of the stilling basin are 8 m, 16 m, 26 m, respectively. The test section (L = 15 m, W = 2 m and H = 1 m) is made of plexiglass. The modelled stilling basin is 3 m long and a gate is adopted to control tail-water level downstream of the model. To measure the fluctuating pressure within the system, 135 measurement taps for pressure transducers are located in the right-side wall and bottom of the stilling basin, as well as in the right-side wall and bottom of spillway (1# spillway in Figure 1) and mid-level channel (1# mid-level channel in Figure 1). For each test, the fluctuating pressure was measured synchronously by sensors for pressure transducers at each measurement point. The accuracy of the pressure sensors utilized in this study is 0.1%, with a 0–50 kPa working range and a frequency response of 100 Hz.

**Figure 1.** Sketch of experimental arrangement (single stilling basin). (**a**) 3D view of gate session; (**b**) plan view; (**c**) side view.

Multiple experimental flooding scenarios were completed within the experimental facility, varying flowrate and changing gate-opening modes. The area under investigation includes the: (i) stilling basin bottom (BB); (ii) stilling basin side-wall (BS); and (iii) piers, spillways and mid-level channels (PSO). Discharges considered varied gradually between 800 and 24,500 m3/s. The downstream water elevation was measured within the model for each test.

### **3. Results and Discussion**

### *3.1. Relation between Fluctuating Pressure and Flood Discharge*

After calculating the root-mean-square (RMS) of the fluctuating pressure measured for each hydraulic condition tested, the arithmetic mean of RMS values obtained for every part of the BB, BS and PSO were calculated individually and these final values are used in the following sections.

Figure 2 shows (a) the maximum and minimum RMS of fluctuating pressure on BB against the flood discharge scenario and (b) the gate-opening modes considered. In Figure 2a, the best fitting line and its correspondent equation are displayed. The configuration of gate-opening modes has a significant impact on the magnitude of the fluctuating load acting on the BB during the release of flows in the range of 2000 to 12,000 m3/s. This is a clear indicator that the fluctuating load can be reduced effectively by optimizing gate-opening modes. As a consequence, the near-field vibrations induced by releasing higher flow associated with flooding scenario can also be reduced. However, this cannot be confirmed when the flood discharge increases to 15,000 m3/s into a single stilling basin (as shown in Figure 2a) due to the limited adjustment of the gate-opening modes available within this experimental facility. Additionally, Figure 2b highlights that the use of side spillways seems favorable to reduce fluctuating pressure on the BB.

**Figure 2.** (**a**) Maximum and minimum RMS of fluctuating pressure on BB and (**b**) corresponding gate opening modes with flood discharge. (Black indicates closed and white indicates open).

Figure 3 displays the maximum and minimum RMS of fluctuating pressure on the BS during flood discharge scenarios, similar to Figure 2. Outcomes displayed in Figure 3 are clear indicators that it is beneficial to reduce pressure values in the BS by reducing the discharge associated with the side spillways as they are directly connected with BS, and hence cause a higher impact. For example, the maximum is more than 7 times higher than the minimum on BS, with *Q* of 12,000 m3/s.

**Figure 3.** (**a**) Maximum and minimum RMS of fluctuating pressure on BS and (**b**) corresponding gate-opening modes with flood discharge. (Black indicates closed and white indicates open).

Figure 4 indicates maximum and minimum RMS of fluctuating pressure on the PSO and corresponding gate opening modes. Although fluctuating pressure on the PSO is larger than those on the BB and BS, the main sources of near-field vibration are the BB and BS. The effect of PSO is relatively small according to Yin and Zhang [18] and Yin et al. [19].

**Figure 4.** *Cont*.

**Figure 4.** (**a**) Maximum and minimum RMS of fluctuating pressure on PSO and (**b**) corresponding gate-opening modes with flood discharge. (Black indicates closed and white indicates open).

### *3.2. Effect of Flow Discharge Ratio and Configuration Adopted to Release the Flow on Fluctuating Pressure*

To quantitatively study the effect of different flowrates, more experiments were carried out. Specifically, three kinds of gate-opening modes were tested at higher flowrates typical of flooding scenarios: (i) mid-level channels separately, (ii) spillways separately, and (iii) combining both previous cases (i) and (ii). Based on measured fluctuating pressure and flow pattern evolution of these three different configurations, recommended gate-opening modes are proposed for various flood discharge rates.

The following notations are used to simplify the understanding of outcomes described in the upcoming sections: *Q* is total flood discharge for single stilling basin (m3/s); *Q*<sup>u</sup> and *Q*<sup>d</sup> are the total discharge of spillways and mid-level channels, respectively (m3/s); *Q*ds and *Q*dm are the total discharge of two side mid-level channels (1# and 5# as shown in Figure 1) and three middle ones (2–4#), respectively (m3/s); *Q*us and *Q*um are total flood discharge of two side spillways (1# and 6#) (m3/s) and four middle ones (2–5#) (m3/s).

Figure 5 illustrates the effect of the ratio *Q*ds/*Q*dm on fluctuating pressure on the BB and BS. It indicates that fluctuating pressure tends to increase on the BS and to reduce on the BB as *Q*ds increases slightly if the mid-level channels are used separately. This result suggests that mid-level channels should be opened uniformly to avoid this disruption.

**Figure 5.** Effect of the ratio *Q*ds/*Q*dm on RMS of fluctuating pressure on BB and BS.

Figure 6 shows the effect of the ratio *Q*us/*Q*um on fluctuating pressure on BB and BS. By slightly increasing *Q*us, if spillways are used separately, Figure 6 shows that the corresponding fluctuating pressure generated increases on the BS and reduces on the BB. This effect is consistent if applied to the

mid-level channels configuration. Therefore, the same as the mid-level channels, spillways should be opened uniformly to reduce the magnitude of pressure on the BB.

**Figure 6.** Effect of the ratio *Q*us/*Q*um on RMS of fluctuating pressure on BB and BS.

Figure 7 displays the effect of the ratio *Q*u/*Q*<sup>d</sup> on fluctuating pressure on the PSO, BB, and BS. When combination of the spillways and mid-level channels is adopted, fluctuating pressure increases on the BS and reduces on the BB and PSO as *Q*<sup>u</sup> increases maintaining constant *Q*. Additionally, by applying this configuration, fluctuating load varies more rapidly on the PSO, and this effect is clearer as *Q* increases. When the ratio *Q*u/*Q*<sup>d</sup> is between 2 and 3, values of fluctuating pressure in all of three areas of study are relatively small.

**Figure 7.** Effect of the ratio *Q*u/*Q*<sup>d</sup> on RMS of fluctuating pressure on BB, BS and PSO.

*3.3. Effect of Flow Discharge Ratio and Configuration Adopted to Release the Flow on Flow Patterns in the Stilling Basin*

In this section, the flow patterns in the stilling basin generated by the different hydraulic conditions and gate-opening configurations tested are discussed. As observed by previous studies [10], there is a correlation between the flow fluctuations and the flow patterns in the stilling basin. This explains why there is a significant difference, considering the same discharge flowrate, between fluctuating pressure on the PSO, on BB, and on BS. Figure 8 displays some examples of flow patterns observed in the stilling basin when the mid-level channels were used separately for diverse discharge flowrates. This is an indicator that various jets tend to quickly merge together and are likely to swing in the stilling

basin if two or three of the middle mid-level channels are opened symmetrically. The formation of these flow patterns generates an increase in fluctuating pressure on the side-walls of the stilling basin. This phenomenon is mainly due to the outlet' width of each mid-level channel that is much smaller than the width of the stilling basin. When replicating this hydraulic condition of using mid-level channels separately, a backflow characterized by an "S" shape emerges in the stilling basin. The flow stability is relative good in the case of combination of 1#, 3#, 5# orifices (Figure 8c), or 1#, 2#, 4#, 5# orifices (Figure 8d), but the flow is not stable when all 1–5# orifices (Figure 8f) are fully open. This result suggests that mid-level channels can only be used individually for small flowrates during flooding events, and mid-level channels should only be opened symmetrically.

**Figure 8.** Flow pattern evolution in the stilling basin with opening of mid-level channels separately for different flow rates. (**a**) Opening partly of 2# and 4# *Q* = 2000 m3/s; (**b**) Opening partly of 2#, 3# and 4# *Q* = 4000 m3/s; (**c**) Opening fully of 1#, 3# and 5# *Q* = 5300 m3/s; (**d**) Opening fully of 1#, 2#, 4# and 5# *Q* = 7000 m3/s; (**e**) Opening partly of 1–5#, *Q* = 7000 m3/s; (**f**) Opening fully of 1–5#, *Q* = 8874 m3/s.

(**e**) (**f**)

Figure 9 displays flow patterns that tend to develop when spillways are opened separately. It shows that the flow pattern is stable when the spillways are uniformly opened and separately regardless of the gate-opening height. Moreover, flow patterns are still stable if the middle spillways

are opened symmetrically by closing side spillways. Finally, Figure 10 shows flow patterns generated when combining discharging orifices and spillways. In this case, flows discharged and flow in the stilling basin fully mix, and the flow patterns produced confirm the stability of the flow in the stilling basin. This stability continues until the outlet section of the stilling basin.

**Figure 9.** Flow pattern evolution in the stilling basin with opening of spillways separately for different flowrates. (**a**) Opening fully of 2–5#, *Q* = 4000 m3/s; (**b**) Opening fully of 1#, 3#, 4#, 6#, *Q* = 8634 m3/s; (**c**) Opening partly of 1–6#, *Q* = 2000 m3/s; (**d**) Opening fully of 1–6#, *Q* = 5868 m3/s.

(**c**) (**d**)

**Figure 10.** Flow pattern evolution in the stilling basin for combining mid-level channels with spillways for different flowrates. (**a**) Opening fully of 1–6# spillways and partly of 1–5# mid-level channels *Q* = 17,800 m3/s, (**b**) Opening fully of 1–6# spillways and 1–5# mid-level channels *Q* = 24,200 m3/s.

To summarize, flow patterns are stable in the stilling basin when combinations of mid-level channels and spillways are applied. The lack of clear and repeatable flow patterns in the stilling basin generated under different scenarios when mid-level channels are opened separately suggest that such a scenario should be avoided. Despite this, there is a consistent good agreement between the flow patterns and the measured fluctuating pressure values acting on the overflow walls.

### **4. Conclusions**

According to the experimental results, the relationship between gate-opening modes and fluctuating pressure on the overflowing walls was explored and optimal gate-opening modes were recommended. Based on the above study, the following conclusions can be stated:


**Author Contributions:** J.Z. and W.X. conceived and designed the experiments; Y.P. performed the experiments; J.Z. and Y.P. analyzed the data; Y.P. and M.R. wrote the paper.

**Funding:** This research was funded by the National Natural Science Foundation of China (51579166) and the National Key Research and Development Program of China (2016YFC0401705).

**Conflicts of Interest:** The authors declare no conflict of interest.

### **Abbreviations**


### **References**


© 2018 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

### *Article* **River Bathymetry Model Based on Floodplain Topography**

### **Ludek Bures 1,\*, Petra Sychova 1, Petr Maca 1, Radek Roub <sup>1</sup> and Stepan Marval <sup>2</sup>**


Received: 16 April 2019; Accepted: 17 June 2019; Published: 20 June 2019

**Abstract:** An appropriate digital elevation model (DEM) is required for purposes of hydrodynamic modelling of floods. Such a DEM describes a river's bathymetry (bed topography) as well as its surrounding area. Extensive measurements for creating accurate bathymetry are time-consuming and expensive. Mathematical modelling can provide an alternative way for representing river bathymetry. This study explores new possibilities in mathematical depiction of river bathymetry. A new bathymetric model (Bathy-supp) is proposed, and the model's ability to represent actual bathymetry is assessed. Three statistical methods for the determination of model parameters were evaluated. The best results were achieved by the random forest (RF) method. A two-dimensional (2D) hydrodynamic model was used to evaluate the influence of the Bathy-supp model on the hydrodynamic modelling results. Also presented is a comparison of the proposed model with another state-of-the-art bathymetric model. The study was carried out on a reach of the Otava River in the Czech Republic. The results show that the proposed model's ability to represent river bathymetry exceeds that of his current competitor. Use of the bathymetric model may have a significant impact on improving the hydrodynamic model results.

**Keywords:** DEM; hydrodynamic modelling; river bathymetry; floods; Bathy-supp

### **1. Introduction**

Knowledge of terrain morphology is crucial for the hydrodynamic modelling of floods. The accuracy and applicability of hydrodynamic models is driven by the nature, availability, and accuracy of source topographic data [1–3].

The digital elevation models (DEMs) are required as a main input for hydrodynamic modelling. Topographic mapping is conventionally conducted by ground surveying. The main advantage of such a method is its high accuracy. Among the major limitations of measured data acquisition are its high cost and time-consuming data collection. Therefore, ground mapping is increasingly being replaced by remote-sensing methods. Radar and laser altimetry are among the most commonly used remote-sensing techniques [4–6].

A description of the river channel and its surrounding area is necessary to create a DEM for purposes of hydrodynamic modelling. The DEM must have precise vertical accuracy and spatial resolution [1,7]. DEMs obtained from satellites are commonly used on a global scale, but the spatial resolution of these models often does not meet the requirements of precise hydrodynamic modelling. The most commonly cited DEMs used for hydrodynamic modelling are the Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER), and the Shuttle Radar Topography Mission (SRTM) [8–11].

Aerial laser scanning (ALS) can be another source of input data for a DEM intended for hydrodynamic modelling. This data source can produce DEMs with high spatial resolution [4]. The method is based on Light Detection and Ranging (LiDAR) technology. ALS methods usually use an infrared laser beam, which is unable to penetrate the water surface or scan the river bed [12]. This problem can be solved by using dual LiDAR (DiAL) technology, which uses a combination of two laser beams [13–15]. However, in some cases, this technology fails because the green laser beam (used for scanning water areas) is not able to penetrate an aquatic environment characterized by high turbidity [16–18].

The most commonly used method for creating a DEM with correct bathymetry is to merge topographic (e.g., ALS) data with another source of bathymetric data. Often used for this purpose are acoustic Doppler current profilers (ADCP), Sound Navigation and Ranging (SONAR) techniques, theodolites, and Global Positioning System (GPS) stations [3,19,20]. The main advantage of these methods is their high accuracy of the acquired data. The main disadvantage lies in the cost of the data's acquisition.

Interpolation approaches have been developed in an attempt to reduce the volume of measured data needed to represent river bathymetry [21,22]. Some interpolation techniques for the creation of bathymetry use cross-sections with various spacings as input data [21,23,24]. Nevertheless, measured data are still required.

Another option for representing river bathymetry can be to use mathematical modelling methods. In works of Dutch and Wang [25] and James [26], the cross-section shape of fluvial sediment deposits was estimated using analytical curves. Moramarco et al. [27] have used entropy-based methods for estimation of the cross-section. Roub et al. [28] introduced a bathymetric model, which estimates the river channel on the basis of hydraulic parameters. The river channel is thereby schematized into a trapezoidal shape. The flow area of this channel is determined using Chezy's equation.

This paper explores new possibilities in mathematical representation of river bathymetry. A new model of river bathymetry based on topographic data describing the area surrounding the river channel is proposed. In this work, the model's ability to represent actual bathymetry is evaluated. A two-dimensional (2D) hydrodynamic model is used to evaluate the influence of the bathymetric model on the results of hydrodynamic modelling. Also presented is a comparison of the proposed model with another state-of-the-art model, the Cross-section Solver ToolBox (CroSolver) [28].

### **2. Materials and Methods**

The proposed bathymetry model is based on analytical curves. The curves are bent into the shapes of the cross-sections. The most precise values of the model parameters are needed in order to obtain the best description of the river bathymetry.

### *2.1. Input Data for the Bathy-Supp Model*

Three types of input data are needed for a river bathymetric model. The first is a DEM, describing the floodplain topography. Also extracted from the DEM are the terrain characteristics, the altitudes for height adjustment of the new cross-sections, and a definition for the aquifer area of the river channel. The second type of input data is the design flow rate. The design flow rate is the flow rate at the time of elevation data acquisition. The last, but not least, type of input data is Manning's roughness coefficient for the river channel.

### *2.2. Construction of the Bathy-Supp Model*

The proposed bathymetric model is constructed in four main steps:

1. The user defines the number and location of the new cross-sections (location of cross-section endpoints) from which a new bathymetry model will be composed.


### 2.2.1. Location of the New Cross-Sections

After displaying the input DEM, it is possible to identify the position of the river channel itself as a no-data region, which thereby divides the model into two (or more) parts. In this no-data region, the user defines the number of new cross-sections and their locations. The distance between the first and last cross-section defines that part of the river for which new bathymetry will be estimated. Each cross-section is defined by two endpoints. The distance between these endpoints defines the width of the channel. The endpoints have coordinates X and Y. A DEM point with the lowest altitude is searched in the circular space around each endpoint. The altitude of this lowest point is used as the Z coordinate of the endpoint. The radius of the circular space is called the lowest point search radius. An example of the cross-section location can be seen in Figure 1.

**Figure 1.** Defining the position of the new cross-sections. (width (W), endpoint (E), lowest point search radius (R)).

### 2.2.2. Explaining Model Parameters

Usually, the best-fit model parameters (m1, m2) are unknown. To establish them (without using an inverse problem-solving method), it is necessary to find the relationship between the search parameters and other explanatory variables. Spatial terrain characteristics can be used as possible explanatory variables [21,26]. The overall curvature, planar curvature, profile curvature, overall slope, slope in the x-direction, and slope in the y-direction were, in the present case, selected for predicting individual model parameters. The characteristics were determined as described by Zevenbergen and Thorne [29]. The source of terrain characteristics was the DEM input. Terrain characteristics were calculated for the all raster cells located behind the bank lines around the river. The terrain characteristics of the nearest raster cell were used to estimate the parameters of a given endpoint. As explanatory variables, the left bank terrain characteristics for model parameter m1 and the right bank terrain characteristics for parameter *m*<sup>2</sup> were used.

Statistical learning methods were used for finding the dependence between model variables and terrain characteristics. Those statistical learning methods studied were multiple linear regression, extended linear regression, and random forest (RF).

A simple linear model (LM) is used for finding a linear relationship between a response and its predictor, and, in the case of multiple linear regressions, this relationship is based upon more than a single predictor. The LM is described in Equation (1):

$$y = \beta\_0 + \sum \beta\_i \mathbf{x}\_i + \varepsilon\_i, \ \varepsilon\_i \sim \mathcal{N}(0, \sigma^2), \ i = 1, \dots, M,\tag{1}$$

where β*x* are the linear parameters, *xi* are predictor variables, ε*<sup>i</sup>* is the error term, and *M* is the number of predictor variables. A mixed selection procedure, as described by Gareth et al. [30], was adopted for choosing the optimal number of variables.

An extended linear model with no random effects (GLS) was also used. This method extends linear regression with an ability to fit models with heteroscedastic and correlated within-group errors, but with no random effects [31]. The extended formula of the LM is described in Equation (2):

$$y = \beta\_0 + \sum \beta\_i x\_i + \varepsilon\_i, \ \varepsilon\_i \sim \mathcal{N}(0, \sigma^2 A^i), \ i = 1, \dots, M,\tag{2}$$

where *Ai* are positive-definite matrices composed using variance and covariance matrices, β*<sup>x</sup>* are the linear parameters, *xi* are predictor variables, ε*<sup>i</sup>* is the error term, and *M* is the number of predictor variables. Again, the mixed selection procedure, as described by Gareth et al. [30], was adopted for choosing the optimal number of variables.

Random forest (RF) is a combined machine learning method for classification and regression. This method is based on an ensemble of a regression tree (RT) algorithm. RT deals with tree structure by dividing the dataset into homogenous groups. That division is driven by some classification criterion, such as minimizing the variance of a given set of variables. In the case of RF, a dataset is divided into multiple sub-datasets by a bootstrap aggregating algorithm. For each sub-dataset, an RF of its own is constructed. This creates a group of random trees, termed RF. For each predictor variable, a measure of variable importance can be determined. Based on variable importance values, it is possible to decide which variables have significant impact for the response and which can be omitted [32].

All statistical analyses in this work were performed using statistical software R. The extended linear model with no random effects was applied using the package nlme [33], and the package randomForest [34] was used for random forests.

The coefficient of determination (*R*2) was used to evaluate the reliability of model parameters. As a second quality assessment, vertical differences between models based on the best model parameters and a model based on estimated parameters were calculated. For this comparison, a similar approach is used in Section 2.1.

### 2.2.3. Cross-Section Construction and Transformation

Once the model parameters are estimated, new cross-sections can be constructed. A tested theoretical cross-section model Equation (3) is able to estimate the natural river cross-section on the basis of estimated parameters [35]. The studied theoretical model of the river cross-section is explained as:

$$z(d) = m\_1 m\_2 d^{m\_1 - 1} (1 - d^{m\_1})^{m\_2 - 1},\tag{3}$$

where *z(d)* is the depth of water at a distance d from the left endpoint of the cross-section, while m1 and m2 are theoretical model parameters that are unique for each river cross-section. An example of the estimated shapes of the cross-section is shown in Figure 2.

**Figure 2.** Example showing ability of the proposed bathymetric model to schematize a cross-section.

Due to the mathematical nature of the proposed model, a new cross-section width of 0–1 [35] must be selected in the first step. Note that the value 0 represents the left side of the cross-section. Additional points are inserted between the points 0 and 1. The user decides upon the number of points to insert. Once new stationing is defined, the depth value for each cross-section point is computed by applying Equation (3). New cross-sections produced by the proposed bathymetric model are in normative state (width and flow area are equal to (1)). Width transformation is simple. Stationing of the new cross-section is multiplied by the distance (W) between its endpoints (Figure 1). For the flow area transformation, an adequate flow area must be identified. This adequate flow area defines the flow area of the river channel required for the transfer of the design flow. The adequate flow area is determined using Chezy's equation on the basis of Manning's roughness coefficient, water surface slope, and design flow rate. This adequate flow area is compared with the area of the new cross-section. If the adequate flow area is smaller than the cross-sectional area, then the depths of the cross-section points are multiplied by the area multiplication parameter. This step is repeated until the adequate flow area is equal to or less than the area of the cross-section. Manning's coefficient, design flow rate, and area multiplication parameter are the input parameters. The water longitudinal surface slope and water surface elevation of each cross-section are extracted from the DEM. A similar approach had been used in the work of Roub et al. [28].

The final transformation step is to add the new cross-section into the coordinate system used. The XY coordinates of the first and last point of the cross-section correspond to the coordinates of the endpoints E1 and E2, for which the cross-section has been created. All internal points are placed on the line connecting the endpoints. Therefore, the coordinates of the internal points can be calculated based on the coordinates of the endpoints and its station value. The method of calculating the coordinates of internal points may vary depending on the coordinate system used. The lower of the altitudes of both endpoints is determined as the water level of this cross-section. The altitude (Z value) of each point is obtained by subtracting the water level and its depth. All cross-section points have coordinates X, Y, and Z and stationing after transformation.

### 2.2.4. River Bed Reconstruction

To create a three dimensional (3D) bathymetric model composed of isolated linear structures, such as cross-sections, spatial interpolation between these cross-sections must be used. Many different bed reconstruction algorithms are available in the literature [23,24,36,37]. In this contribution, we adopted the approach of Caviedes–Voullième [23]. This bed reconstruction algorithm is based upon cubic Hermite splines (CHS). Spline trajectory connects two depth points in two consecutive cross-sections, and it is driven by their normal vectors. The spatial interpolation (X, Y) of new bathymetric points is performed using CHS, and linear interpolation is performed for the height (Z). For a more detailed description, see the work of Caviedes–Voullième [23].

### *2.3. Model Suitability*

Global optimization methods are used for estimating the best model parameters. The global optimization schemes are based on heuristics inspired by natural processes. Methods of differential evolution [38] are used for determining the solutions of inverse problems related to the parameters of a mathematical model of river bed surfaces. Differential evolution is a population-based stochastic optimization search algorithm that iteratively estimates a candidate solution with regard to a given measure of quality [38]. The analyzed objective function is a least squares method, determining the differences between the depths of the real cross-sections and the modelled cross-sections. We employed the *best1bin* algorithm in this work [38]. Analysis of the model's suitability was made in the C++ programming environment.

### Model Suitability Evaluation

For evaluating the model's suitability, a vertical differences comparison between the model cross-sections (with the best model parameters) and the measured cross-sections (see Section 2.5.2) was made. The root mean square error (RMSE) and the mean absolute error (MAE) were calculated for this purpose. The equations for these evaluation criteria are as follows:

$$RMSE = \sqrt{\frac{1}{N} \sum\_{i=1}^{N} \left( Elev\_{MOD} - Elev\_{REF} \right)^2},\tag{4}$$

$$MAE = \frac{1}{N} \sum\_{i=1}^{N} |Elev\_{MOD} - Elev\_{REF}| \tag{5}$$

where *ElevMOD* represents the elevation (m) obtained from each model cross-section and *ElevREF* represents the equivalent reference point obtained from the measured cross-section (see Section 2.5.2). *N* represents the number of the cross-section points.

### *2.4. Case Study Area*

A part of the Otava River in the Czech Republic was chosen as an area of interest (Figure 3). The full length of the Otava River is 111.7 km with the basin area of 3840 km2. The studied river reach is located into the lower part of the river (between 22.83 and 24.58 river km) and is 1.75 km long. The average depth of the river reach fluctuates around 1 m. The average bankfull width is between 22.8 and 52.7 m. The average annual flow rate is 23.4 m3/s. The average water level is 354.84 m above sea level. The flow rates for N–year floods in the Otava River reach are shown in Table 1. Qa.a. denotes the average annual discharge.

**Table 1.** The flow rates for N-year floods in the Otava River reach [39].


**Figure 3.** The study area: Otava River reach.

### *2.5. Ground and Bathymetry Data*

This section of the paper describes the topographic source data and bathymetry data source used for identification of the model parameters and creating the DEMs.

### 2.5.1. Aerial Laser Scanning Data

The technology used for ALS data collection in this study was aerial laser scanning with one infrared laser beam. These ALS data were provided by the Czech Office for Surveying, Mapping, and Cadastre. Provided data were in the form of the list of elevation points. The point density per square meter varies depending on the slope of the terrain described. The point density starts with one point (flat area) and ends with dozens of points (steep area). Vertical data accuracy is +/−0.18 m in open areas and +/−0.3 m in forest areas. [40].

### 2.5.2. Acoustic Doppler Current Profiler Data

The data obtained by using an ADCP measurement device have the depth range 0.2–80 m and accuracy of +/−1% of the depth. The ADCP technology was successfully applied in recent works [18,41,42]. The data were sampled as a set of single cross-sections. The distance between these cross-sections was approximately 5 m. The distance between points within the cross-sections was less than 0.5 m. A total of 375 cross-sections were measured. These data were used for testing the suitability of the bathymetric model as well as for determining the best model parameters.

### 2.5.3. Compared DEMs

A total of four DEMs, based upon different sources of bathymetric data, were compared in this practical demonstration. The first DEM (ALS) was created only from the ALS data. It contains no bathymetric data and represents models of that sort produced by remote-sensing techniques.

The second DEM (CRO) was composed using ALS data, which describe the floodplain, and data from the Cross-section Solver ToolBox (CroSolver) bathymetric model. CroSolver is a software tool for improving an existing DEM in the riverbed area. The newly formed channel is schematized as trapezoidal or rectangular. The Chezy equation is used to determine the channel flow area for the design flow rate. The depth of the new channel is determined on the base of the size of the flow area, channel width, and selected channel schematization [28]. Table 2 provides the settings overview for the CroSolver model.


**Table 2.** The Cross-section Solver ToolBox (CroSolver) software specific settings.

The third DEM (BAT) was composed using ALS data, which describes the floodplain and the Bathy-supp model proposed in this paper. The cross-section distance for the bathymetry creation was about 100 m. Table 3 provides the settings overview for the proposed model.



The last DEM (ADP) is a model composed using ALS data and ADCP data. This model is used as a reference model because it is based on measured data (i.e., on data with assured accuracy). Table 4 provides background information on the DEMs being compared.


**Table 4.** Overview of compared digital elevation models (DEMs).

### 2.5.4. DEM evaluation

For evaluation of the DEM quality, the vertical differences between compared DEMs (ALS, CRO, BAT), and the reference DEM (ADP) were assessed. This comparison was made for all corresponding cells in the river channel area. Comparison of the models in the floodplain is not relevant. Identical data were used for its description. RMSE Equation (4) and MAE Equation (5) values were employed for evaluation.

### *2.6. Hydrodynamic Modelling*

All hydrodynamic simulations were performed using a 2D HEC-RAS model [43]. HEC-RAS 2D allows computations for steady and unsteady flow conditions. The hydrodynamic results are determined on the base of solving of 2D Saint–Venant equations. Many authors successfully employed HEC-RAS 2D for hydrodynamic simulations [44,45].

In this study, hydrodynamic simulations for four different DEMs were made (see Section 2.5.3). For each DEM, four simulations were made with chosen N-year flow rates. Manning's roughness coefficients were set separately for each land cover category. Values in the range 0.035–0.10 s/m1/<sup>3</sup> were determined for inundations and the value for the main channel was 0.031 s/m1/3. Selection of the Manning's values was verified by a calibration-verification process. The validation of model parameters was based on a flood event from December 2012, where, for discharges of 143 m3/s, the recorded water level (hydrological station located in the model river reach) was equal to 356.43 m. The normal

depth was used as the downstream boundary condition. The results of the basic hydrodynamic model (topography source was ADP DEM) for the given discharge provided a difference of 2 cm in the water level. That verified the accuracy of the model setup. All simulations began from a minimal start discharge, which then grew until eventually becoming steady [18,46,47]. Table 1 shows the selected N-year floods that were used as boundary conditions.

Water surface elevation (WSE) and Inundation areas (IAs) were evaluated to determine the influence of channel bathymetry representation. WSE was evaluated by comparing the differences between vertical measurements for the raster-based WSE (ALS, CRO, BAT) and those of the reference WSE (ADP). The RMSE (Equation (4) and MAE (Equation (5) values were used in the evaluation.

For IAs, the following compliance criterion was used:

$$IA\_{dif} = \frac{|IA\_{DEM} - IA\_{REF}|}{IA\_{REF}} \times 100\tag{6}$$

where *IAdif* represents the difference in inundation areas (%), *IADEM* represents the inundation area (km2) of compared models (ALS, CRO, BAT), and *IAREF* represents the inundation area of the ADP model (km2).

### **3. Results and Discussion**

### *3.1. Bathymetric Model Suitability*

The ability of the proposed bathymetric model to schematize the measured cross-section was evaluated. The best model parameters were used for this purpose. The parameter range m1 ranged from 1.0202 to 1.8219, and for the *m*<sup>2</sup> parameter from 1.0929 to 2.0376. Overall, 375 measured cross-sections were evaluated. The mean RMSE and MAE values were 0.16 m and 0.11 m, respectively. The variability of the RMSE and MAE values is shown in Figure 4.

**Figure 4.** Comparing ability of the bathymetric model (with ideal parameters) to represent a measured cross-section. Shown are variances of root mean square error (RMSE) and mean absolute error (MAE) values.

### *3.2. Explaining Bathymetric Model Parameters*

High-quality coefficient estimation is a prerequisite for correct bathymetric model creation. Therefore, various methods for its estimation were evaluated. For LM and GLS methods, overall curvature, planar curvature, overall slope, slope in the x direction, and slope in the y direction were considered as predictor variables. For the RF method, overall curvature, planar curvature, profile curvature, overall slope, slope in the x direction, and slope in the y direction were considered. LM and GLS methods provided similarly poor results. In contrast, the RF method produced a good result. Table 5 provides an overview of determination coefficients for the compared techniques. In view of these findings, the RT method was adopted for creating the bathymetric model.

**Table 5.** Coefficients of determination for the best-fit model parameters in comparing the estimation techniques linear model (LM), the extended linear model with no random effects (GLS), and the random forest (RF).


Differences in parameter estimation quality may be due to nonlinear relationships in the data structure. More detailed analysis of the suitable regression structures and evaluation of its results are planned in follow-up research.

### *3.3. DEM Comparison*

This comparison was made for cross-sections extracted from compared DEMs. The smallest divergence from the reference model (ADP) was achieved by the BAT model. The CRO model is unable to closely simulate the depth variability across the cross-section. This is due to the trapezoidal schematization of the channel. [28]. The ALS model almost completely neglects the riverbed area, which is specific for the sampling method used [4,12,40,48]. A graphical comparison is shown in Figure 5.

**Figure 5.** Visual comparison: Random cross-section derived from evaluated digital elevation models.

Table 6 describes the RMSE and MAE values achieved for the river channel (raster comparison). Both evaluated errors for the ALS model are >1. Note that the mean water depth in the channel at the time of acquiring ALS data was around 1 m. This suggests that the ALS model neglected almost the entire flow area of the river channel [16,20,28,42]. The RMSE and MAE values for the CRO model were 0.46 and 0.36 m, respectively. The smallest error values, RMSE 0.30 and MAE 0.23 m, were achieved by the BAT model. The raster cell difference variability is shown in Figure 6. Greater difference variability was achieved by the ALS model, and the smallest variability was achieved by the BAT model.

**Table 6.** RMSE and MAE errors for the compared DEMs (channel comparison).


**Figure 6.** Variance of raster cell differences among compared digital elevation models.

Channel difference maps are shown in Figure 7. The biggest differences can be seen for the ALS model. The largest difference values are located along the thalweg (or stream centerline) location. Indeed, description of the river bottom was unrealistic. The error map presented for the CRO model identifies some areas with high cell difference values. These areas are located near embankment areas. The differences would probably be less significant if the cross-section were actually to have a trapezoidal shape (reflecting technical modification of the channel). For a natural river, the differences may be more significant. The BAT difference map provides the best results. Even here, however, it is possible to identify areas with a poor match. In this case, there was an excessive recess of one cross-section in the BAT model. The water surface slope derived from inundation data was slightly underestimated. This may be due to the accuracy of the data used to describe the floodplain [40].

**Figure 7.** Error maps comparing evaluated digital elevation models (ALS, CRO, BAT) and the reference model (ADP).

### *3.4. Thalweg Comparison*

The ALS model provided a poor match with the ADP model. Thalweg derived from the ALS model was shifted toward the water surface, where it fluctuated. Thus, it is evident that almost the entire flow area was neglected [16,20]. The CRO model had higher deviations relative to the BAT model. These deviations may occur in places where the channel narrows or expands locally. This is due to the fact that the local changes in flow velocity are not reflected in the CroSolver software [28]. The graphical comparison is shown in Figure 8, and the results of the numerical comparison of thalwegs are shown in Table 7.

**Figure 8.** Visual comparison of thalwegs.

**Table 7.** RMSE and MAE values for the compared DEMs (thalweg comparison).


### *3.5. Water Surface Elevations Comparison*

Figure 9 presents the variability in WSE errors between the compared DEMs (ALS, CRO, BAT) and the ADP DEM. The greatest deviations in WSE were seen in the ALS model, which manifested significant overestimation of WSE at all modelled flow rates. At flow Qa.a., the RMSE was more than 1.2 m, although the WSE variability was comparable to that for other models. For models CRO and BAT, the medians of the differences, as well as their variability, decreased with the increasing flow rate. The BAT model provided the smallest median values and the smallest variance for all rated flows. Table 8 presents RMSE and MAE values for WSE comparison.



**Figure 9.** Variance of raster cell differences in water surface elevation between the compared digital elevation models (ALS, CRO, BAT) and the reference model (ADP).

### *3.6. Inundation Areas Comparison*

Table 9 presents a comparison of inundation areas. The ALS model in all cases overestimated inundation vis-à-vis the reference ADP model. This overestimation was >50% in the case of Qa.a. These differences were diminishing with increasing flow, but nevertheless were >20% for the Q100 flow rate. Similar results had been presented in the works of Bures et al. [18] and Roub et al. [42]. In both those cases, software-modified DEMs (CRO, BAT) provided better results than the ALS model. The CRO model underestimated the ADP model values by as much as 4.5%. The BAT model underestimated them by as much as 3.8%. The fact that CRO and BAT consistently underestimate the results may reflect their similar model settings. The flow area for model cross-section transformation was set the same for the two models.


**Table 9.** Inundation area (IA) differences for compared DEMs.

Overall, the ALS model showed the poorest results. That is because the technology used is unable to capture terrain under water. It was a LiDAR technology that, in our case, used an infrared laser beam [40]. Similar poor results could be expected from other remote-sensing methods [8–11]. The error produced by these models depends on the size of the neglected flow area of the river channel [18,42]. Although differential absorption LiDAR (DiAL) technologies may constitute exceptions to this general case [14,15], even using DiAL technology does not ensure any increase in accuracy. That is particularly the case when a high-turbidity river channel is measured [16–18,20].

The ALS model used in this study is a basic DEM, representative of the group of other remote-sensing models.

In the BAT model, the longitudinal water surface slope and water surface elevations were determined for all constructed cross-sections. Both of these parameters are derived from DEM inputs. The uncertainty in determining these parameters depends upon the accuracy of the data from which they are derived. It is possible, however, to assume that the development of remote-sensing techniques will reduce the uncertainty in determining water surface slope and elevations. Uncertainty in determining roughness parameters is the same for both models (CRO, BAT) and depends upon user experience. Further research in the field of model parameter estimation may bring improvement in this area. More appropriate regression relationships can be found, as can other explanatory variables. These variables might, for example, be the curvature of the flow centerline [21,38], mutual tilt of consecutive cross-sections [23], or ratio of the bank line distances.

In comparing the quality of the DEMs, we can see that the BAT model produced better results than the CRO model. The extent of the differences between the CRO and BAT models will probably depend upon the nature of the channel described. With respect to its adaptability, it can be assumed that the BAT model will produce better results in the case of natural channels. If the river channel has been technically adjusted (into trapezoidal shape), it can be assumed that the BAT model will produce similar results to a CRO model.

In the present bathymetric model, the interpolation mechanism introduced by Caviedes– Voullième [23] was used. Several interpolation methods can alternatively be used for this step [24,37,38]. Which of these methods will yield more reliable results remains an open question.

In the presented study, the importance of river bathymetry was evaluated mainly from the perspective of flood modeling. However, the proposed bathymetric model can find its place in other scientific disciplines, such as understanding the morphological changes of the river segment or understanding the hydrodynamic behavior of the river.

### **4. Conclusions**

In this study, the new theoretical bathymetric model Bathy-supp was presented. The model is based on the analytical curves, which schematizes the river cross-sections. The shape of the analytical curves is driven by the floodplain topography. In the case study, the practical usability of the model was evaluated.

The vertical differences in the model's ability to represent the cross-sections were 0.16 m for RMSE and 0.11 m for MAE. For estimation of the model parameters, the three regression methods were used. The best parameter estimates were provided by the RF method. However, the regression relationships between model parameters and topographical characteristics need to be further investigated.

The study's results also show that, in the DEM quality assessment, the BAT model provides the best results in the river channel and thalweg representation. When assessing the WSE and IA, the BAT model provided better results than the CRO model, especially for high design flows rates.

Both DEMs created by merging ALS data with mathematical bathymetric models (CRO and BAT) provided significantly better results. The extent of this improvement was directly proportional to the size of the flow area neglected by using the ALS data itself.

The newly proposed Bathy-supp model was also compared with the state-of-the-art CroSolver model. The results showed that the Bathy-supp model describes the river bathymetry more accurately than the CroSolver model. This is because the analytical curves used in the Bathy-supp tool are more able to simulate the individual cross-sections of the channel than the trapezoidal shapes used by the CroSolver tool. Another advantage of the Bathy-supp model is in determining the flow area for each single cross-section from which the model is composed. This allows the model to take into account local narrows or expansions of the river channel.

The merging of bathymetric models with other types of remote-sensing data can significantly improve the DEMs, thereby enhancing their practical usefulness. Hence, the methods for mathematical schematizing of riverbeds should be further studied and improved.

The Bathy-supp model as an appropriate mechanism to improve DEM quality when other bathymetry measurements that are not available can be used. However, it will never be able to fully replace large bathymetric measurements.

**Author Contributions:** Conceptualization, L.B. and P.M.; methodology, L.B., P.M., and R.R.; validation, L.B., P.S., and R.R.; investigation, P.S.; data curation, L.B.; writing—original draft preparation, L.B.; writing—review and editing, P.S., S.M.; visualization, L.B. and S.M.; supervision, P.M. and R.R.; funding acquisition, R.R.

**Funding:** This research was supported from the European Union by the Operational Programme Prague–Growth Pole of the Czech Republic, project No. CZ.07.1.02/0.0/0.0/17\_049/0000842, Tools for effective and safe management of rainwater in Prague city–RainPRAGUE and the Internal Grant Agency (IGA) of the Faculty of Environmental Sciences (CULS) (IGA/20164233).

**Conflicts of Interest:** The authors declare no conflict of interest.

### **References**


© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

### *Article* **Application of Python Scripting Techniques for Control and Automation of HEC-RAS Simulations**

### **Tomasz Dysarz**

Department of Hydraulic and Sanitary Engineering, Poznan University of Life Sciences, ul. Wojska Polskiego 28, 60-637 Poznan, Poland; dysarz@up.poznan.pl or tdysarz@gmail.com; Tel.: +48-061-846-6586; Fax: +48-061-848-7726

Received: 1 August 2018; Accepted: 30 September 2018; Published: 2 October 2018

**Abstract:** The purpose of the paper was to present selected techniques for the control of river flow and sediment transport computations with the programming language Python. The base software for modeling of river processes was the well-known and widely used HEC-RAS. The concepts were tested on two models created for a single reach of the Warta river located in the central part of Poland. The ideas described were illustrated with three examples. The first was a basic simulation of a steady flow run from the Python script. The second example presented automatic calibration of model roughness coefficients with Nelder-Mead simplex from the SciPy module. In the third example, the sediment transport was controlled by Python script. Sediment samples were accessed and changed in the sediment data file stored in XML format. The results of the sediment simulation were read from HDF5 files. The presented techniques showed good effectiveness of this approach. The paper compared the developed techniques with other, earlier approaches to control of HEC-RAS computations. Possible further developments were also discussed.

**Keywords:** river flow modeling; sediment transport simulation; automation of flow modeling; HEC-RAS controller; python scripting

### **1. Introduction**

River flow modeling is a very popular approach in science and education, as well as in small, medium-size, and large technical projects in the areas of water management [1,2], river regulation [3,4], sediment transport [5], flood protection [6], and many others [7,8]. There are several professionally prepared and widely used commercial and non-commercial models for river flow, e.g., HEC-RAS [9], MIKE 11 [10], Delft Sobek [11], BASEMENT by ETH [12], and SRH1D [13]. The first of the mentioned models, HEC-RAS, is very specific. The acronym HEC means Hydrologic Engineering Center, and it is the name of the software developer. The second element, RAS, means River Analysis System. It defines the software application area which is focused on modeling of flow and transport processes in rivers, floodplains, and reservoirs. All of the modeled elements are solved with highly accurate numerical methods [9]. Additionally, HEC-RAS is freeware with a professionally developed graphical user interface. It makes application of the software easier, and it is considered to be the reason that HEC-RAS is so popular. Today a strong focus on modeling of uncertainty related to river processes is observed. It is a particularly crucial problem in the areas of model calibration [14,15], sediment transport [16,17], and flood hazard mapping [18,19]. The effective management of such computational tasks is not possible without automation of simulation and integration with other tools, e.g., GIS software. The examples analyzed in the paper present crucial automation techniques for the mentioned applications.

The main purpose of this paper is to present selected techniques for controlling HEC-RAS computations with the programming language Python. Control of river flow simulations with external scripts opens several new opportunities in the fields of flood hazard assessment, water management, and river hydraulics. It also enables application of more sophisticated methods in these fields, such as sensitivity analysis, automatic model calibration, and probabilistic flood mapping. The control of the HEC-RAS computations is possible due to the fact that this package is installed with the Component Object Module (COM). This COM library of methods is called HECRASController. The COM library enables outer access to the computational elements of HEC-RAS. The discussion presented here is focused on three specific examples. The ideas presented are not limited to the basic running of flow model simulation. Some concepts are aimed at integration of HEC-RAS with specific tools available in Python modules, e.g., SciPy, while others extend the capabilities of HECRASController and open access to more sophisticated data formats, e.g., XML (eXtensible Markup Language) and HDF5 (Hierarchical Data Format 5).

The HECRASController is a collection of Visual Basic subroutines and functions within the HEC-RAS package. It may be effectively used with VBA (Visual Basic for Applications) in Excel [20,21]. The main capabilities of this collection include running HEC-RAS and particular HEC-RAS editors, running computations plans, reading results of flow simulation, etc. One very interesting feature is the possibility of manipulation of roughness coefficients. It is a very important problem studied by many researchers, e.g., Yang et al. [22] and Lacasta et al. [23]. This capability is implemented in AHYDRA software [24]. The AHYDRA supports calibration of roughness coefficients and boundary conditions for steady flow conditions. The software cooperates with the initially configured HEC-RAS model, enabling automatic change of settings.

The capabilities of HECRASController are still limited. Even the main author promoting the use of this programming tool, Goodell [20,21], advises manipulation of the HEC-RAS data files in ASCII format. Another problem is the use of Visual Basic/VBA. Although it is a very handy programming language, its use is rather limited and still narrowing, e.g., VBA was replaced in ArcGIS by Python and its extension ArcPy. Hence, other languages are also implemented to handle HECRASController, e.g., MATLAB [25].

Python [26] is one of the most popular programming languages used today worldwide. It is easily applied in practical cases, as well as in highly scientific problems [27,28]. Its popularity follows from two factors: (1) relative simplicity and (2) a broad range of application areas. Considering the specific problem considered here, Python has several advantages and unique features. They include very easy access to text files in standard form, e.g., [26]; several specific libraries included, e.g., [29–33]; opportunities to handle the components of the Windows operation system, e.g., [34,35]; access to sophisticated data formats, e.g., [36,37]; and linking with a number of Windows application, e.g., [38,39].

It is no surprise that the first approach to control of HEC-RAS with Python was developed. The PyRAS module prepared by Peña-Castellanos [40] is one such approach. Unfortunately, PyRAS only transforms objects and functions of HECRASController available in VBA to Python. Hence, the limitations of HECRASController still apply. Another disadvantage is that the module was prepared for older versions of HEC-RAS, namely 4.1 and 5.0.0. It has not been developed further to date. Despite these facts, the codes of PyRAS could be an inspiration for independent development. Another open-source application linking Python and HEC-RAS is the example included in the package PyFloods [41]. Unfortunately, it is the only example for the previous version of HEC-RAS with some important elements hard-coded.

No general approach to linking Python and HEC-RAS has been proposed as yet, and it seems to us that it should be done on the basis of available Python modules for management of COM libraries. Another requirement is to develop the potential of HECRASController with additional Python modules, e.g., optimization of model parameters, access to specific data in XML files, and reading and handling of results in HDF5 files. These problems are addressed in the present paper.

The next section, Methods, includes general description of all important elements in the presented manuscript. Theoretical aspects of open channel flow, as well as Python scripting issues related to the applied codes, are discussed. In the section Results and Discussion, examples illustrating the invented concepts are described in detail. At the beginning of this section, the HEC-RAS models used are presented. In the final section, the summary, concluding remarks, and further developments are discussed.

### **2. Methods**

### *2.1. HEC-RAS*

HEC-RAS is a well-known hydrodynamic model for rivers and water reservoirs. The concepts applied in the package are well described by Brunner [9]. The HEC-RAS is applied for simulation of flow and transport processes in river networks, including floodplains and reservoirs. The modeled flow conditions include steady and unsteady longitudinal flow. The first is based on the simple equation of energy balance. The form of this equation implemented in HEC-RAS is shown below:

$$z\_1 + h\_1 + a\_1 \frac{\mu\_1^2}{2g} = z\_2 + h\_2 + a\_2 \frac{\mu\_2^2}{2g} + h\_{c\_2} \tag{1}$$

where *z* is bottom elevation, *h* is the depth, and *u* is mean velocity in the channel cross-section. *α* is called the St. Venant coefficient and plays the role of correction factor, including the effect of velocity profile non-uniformity. *g* is well known acceleration of gravity. The Equation (1) is written for gradually-varying flow, when the assumption of hydrostatic pressure distribution may be suitable. The subscripts 1 and 2 denote two different cross-sections in the same channel reach. The basic assumption is that the cross-section number 1 is located upstream of the cross-section number 2. The first three terms of both sides represent the potential energy of the stream, work of pressure forces, and kinematic energy of the stream. The last element of the right side, *he*, describes friction losses due to bed and banks influence on flowing water. It also includes effects of channel contraction and extension. If the depth is known in one cross-section, on this basis, the depth may be also determined in the second cross-section. To calculate the distribution of the depth along the channel, its value must be known in one cross-section a priori. In practical cases, this cross-section is the inlet or outlet boundary of the channel reach. Hence, the condition is frequently called "boundary condition", but it is rather hydraulic jargon than strict mathematical terms. The discharge *Q* is a very important parameter of the Equation (1). The kinematic energy terms, as well as friction losses, depend on the magnitude of the flow. The influence of floodplains is included in the calculation of the St. Venant coefficients and calculation of weighted distance between cross-sections. The last element is important for the determination of the total losses *he*. More detailed description of this equation and its implementation may be found in Reference [9].

For description of the second model, the unsteady flow, the numerical solution of the St. Venant system of equations for 1D flow is implemented. This system consists of two partial differential equations shown below.

$$\frac{\partial A}{\partial t} + \frac{\partial (\phi \mathcal{Q})}{\partial \mathbf{x}\_{\mathcal{E}}} + \frac{\partial [(1 - \phi)\mathcal{Q}]}{\partial \mathbf{x}\_{f}} = 0 \tag{2}$$

$$\frac{\partial \mathcal{Q}}{\partial t} + \frac{\partial}{\partial \mathbf{x}\_{\varepsilon}} \left( \phi^2 \frac{\mathcal{Q}^2}{A\_{\varepsilon}} \right) + \frac{\partial}{\partial \mathbf{x}\_{f}} \left[ (1 - \phi)^2 \frac{\mathcal{Q}^2}{A\_{f}} \right] + \mathcal{g}A\_{\varepsilon} \left( \frac{\partial H}{\partial \mathbf{x}\_{\varepsilon}} + \mathcal{S}\_{f\varepsilon} \right) + \mathcal{g}A\_{f} \left( \frac{\partial H}{\partial \mathbf{x}\_{f}} + \mathcal{S}\_{ff} \right) = 0 \tag{3}$$

The Equation (2) describes mass balance in the open channel flow. The second equation represents the momentum balance. It is written for gradually-varying flow conditions in the channel. There are two independent variables: *t*—time and *x*—distance. The distance is measured not only along the channel, but also along the floodplains. Hence, there are denotations *xc* and *xf* for channel and floodplain flow paths, respectively. It is important that the pair of distance (*xc*, *xf*) is unique for each cross-section. *Q* is still the total discharge in the cross-section and *g* is acceleration of gravity. *H* is the water surface elevation, which is the sum of the bottom elevation *z* and the depth *h* shown in the previously discussed Equation (1). It is assumed that the transversal slope of the water surface is zero in the single cross-section. Hence, *H* is constant in the cross-section, but varies along the channel. Denotation *A* without subscript is the total cross-section area in a single cross-section. However, *Ac* and *Af* are the area of channel and floodplain parts of the cross-section, respectively. The hydraulic slopes describing the friction losses along the channel and along the floodplains are denoted as *Sfc* and *Sff*. The HEC-RAS specific element is coefficient *φ*. It describes the distribution of the discharge between the parts of the cross-section. Hence, *φQ* is the flow along the channel and (1 − *φ*)*Q* is the flow along the floodplains. The coefficient *φ* depends on the hydraulic parameter called conveyance, which is the function of the water surface elevation *H* and other elements of the cross-section, such as shape and roughness. In the Equation (2), the first term represents the local storage of water, whilst two other terms describe the net inflow of water along the channel and floodplains. In Equation (3), the first three terms represent the inertia forces, the terms including gradient of *H* describe the gravity and pressure forces, and the terms with hydraulic slopes *Sfc* and *Sff* represent the friction forces. The solution of the system (2)–(3) are two functions of time *t* and distance pair (*xc*, *xf*). These are discharge *Q* and water surface elevation *H*. To solve the system (2)–(3), the initial and boundary conditions have to be specified. The basic forms of the initial condition are known distributions of *Q* and *H* along the modeled reach. In each inlet and outlet cross-section, there should be one boundary condition imposed in gradually-varying flow. These conditions may have a basic form of discharge or water stage hydrographs. They can also be imposed as other more complex relationships of discharge and water stage, e.g., normal depth, rating curve, etc. The Preissmann scheme is used for the approximation. More detailed description of the St. Venant equations used in HEC-RAS and their implementation may be found in Reference [9].

The number of hydro-structures, such as bridges, culverts, etc., is also available for incorporation with the prepared water system. Both these mathematical models of open channel flow are well described in many books [42–45], as well as software manuals mentioned above [9–13]. In HEC-RAS, there is also a module enabling so-called quasi-unsteady flow simulation, i.e., simplified flow simulation in unsteady conditions on the basis of the energy equation. In the latest versions of the package, a 2D flow module is also available. The HEC-RAS modules for simulation of transport processes include different transport of solutes, heat, and sediments with deposition and erosion. Heat and solute transport is based on numerical solution of convection–dispersion equations with source terms. Models of this kind are described in some books on mathematical modeling of river processes, e.g., References [9,46,47]. The QUICKEST scheme with the ULTIMATE limiter is applied as an approximation method [48]. The sediment transport is modeled with the standard Exner equation, enabling simulation of different fractions of sediments. The Exner equation is shown below:

$$(1 - p)B\frac{\partial z}{\partial t} + \frac{\partial Q\_s}{\partial x} = 0 \tag{4}$$

where *x* is the distance along the channel and *t* is time. *z* represents bottom elevation, *B* is the width of this bottom part, where the sediments are deposited or removed from. *p* is porosity of the bed layer. *Qs* is volumetric intensity of sediment transport. This element is described by subjectively chosen empirical formulae, e.g., Meyer-Peter and Müller (MPM). It depends on the hydraulic parameters and the sediment characteristics. The first term represents the local change of the sediment deposited. The second term is the net inflow of sediments. The solution of Equation (4) is the change of bottom elevation *z* in time *t,* and along the channel *x*. The initial and one boundary condition have to be imposed to solve the problem. The initial condition is simple initial bottom elevation profile along the modeled channel. The boundary condition is imposed in the inlet cross-section. It may have the form of bottom elevation changes, sediment inflow intensity varying in time, or equilibrium sediment load in the stream. In HEC-RAS, the applied solution method is the basic upwind scheme [9,45]. More details may be found in Reference [9].

The HEC-RAS package also includes several useful tools for data preparation and results processing. These tools include the module for GIS data processing called RAS Mapper [49].

### *2.2. HECRASController*

HECRASController is a part of the HEC-RAS application programming interface (API). It is a collection of programming tools: Classes, functions and subroutines. Access to the HEC-RAS elements is possible because this package is compiled as the Component Object Module (COM). Any program able to read COM DLL (Dynamic-Link Library) may be used to control HEC-RAS computations. The most convenient approach is to use Visual Basic in any form, i.e., as a separate suite for developers or linked with another application, e.g., MS Excel.

The most detailed description of the above-mentioned collection of programming tools was presented by Goodell [20,21]. Appendix A of his book [20] is a very good source of information about syntax, and usage of HECRASController functions and subroutines. The examples presented there, illustrate the great potential of these concepts. VBA in MS Excel is very powerful in computing, as well as in data analysis [50]. Linking HECRASController with VBA in Excel opens new opportunities for handling simulation data and results, as presented by Goodell [20,21]. The controller makes it possible to open HEC-RAS projects, run simulations, read results, and store them in other specific formats. It is also possible to open HEC-RAS editing tools for manual configuration of input data. The controller also enables direct manipulation of project data, but its capabilities in this area are rather limited. There are three procedures for changing roughness coefficients, and one for setting the area of the element called storage area.

Although this concept is not new and seems to be very useful, there are still areas of HEC-RAS computations in which controller functions are not able to help. The basic HEC-RAS data, such as geometry and flow boundary conditions, are stored in ASCII files. They may be accessed by manually written code. Goodell [20] has advised such an approach. However, some HEC-RAS data are stored in formats other than pure text, e.g., sediment data in XML format. Access to such data requires application of additional modules. The results of more complex simulations are also beyond HECRASController access, e.g., 1D sediment and 2D flow simulation results are stored only in HDF files.

In addition, the opportunities afforded by the controller depend on the application areas of VBA. This powerful language is widely used. However, there are also very important areas where VBA has been withdrawn, e.g., scripting in ArcGIS. There are also areas in which VBA has never been present, e.g., scripting in QGIS. Relying on VBA for automation of HEC-RAS simulation could limit possible applications.

### *2.3. Python Scripting and Applied Modules*

Python is one of the most popular programming languages today, and it is still under development. In this paper, version 2.7.12 was used. Its usefulness and popularity in many areas have been reported in References [51–53]. This scripting language is relatively simple for the beginner, but it is also very powerful if applied by an experienced coder. A description of this language may be found in many books or on internet websites, e.g., Downey et al. [54], Python Software Foundation [26]. Considering the problem considered here, Python has several specific advantages:


(5) There are Python modules integrating this language with geoprocessing software, such as ArcGIS [38] and QGIS [39].

VBA is a compiled language, which makes it faster in comparison with any scripting language in stand-alone applications. However, in the analyzed cases, the programming language is used to run an external application, HEC-RAS. Hence, the efficiency of the external application determines the efficiency of the whole concept, which reduces this advantage of VBA. The above list of Python features is the list of undoubted advantages of this scripting language over standard VBA. In some cases, the application of Python is as easy as VBA, e.g., access to COM libraries or XML data. Other elements, such as libraries for numerical computations, make Python implementation more effective. There are also features that are unique to Python, e.g., access to HDF5 data, and linking with Fortran and C. The undeniable advantage of Python is also clearly visible in the area of geoprocessing. This area is important for any hydrological or hydraulic modeling of flow and transport processes in natural water systems. As mentioned above, the former applications of VBA in older versions of ArcGIS are now replaced by modern scripting languages, such as Python [38], JavaScript [55], and R [56]. Python is the primary scripting language used in commercial Esri software, such as ArcGIS 10.x and ArcGIS Pro. This is also the only scripting language implemented in QGIS, which is the open source and cross-platform equivalent of Esri's packages.

Several Python modules are used in this research, including: (1) Pywin32, (2) NumPy, (3) SciPy with optimization, (4) ElementTree XML API, and (5) HDF 5 for Python. The first of them, the Pywin32 module, is used to access COM objects and COM servers in Python code [34,35]. Brief descriptions of available Pywin32 packages may be found in the PythonCOM Documentation Index [35]. The part of the module win32com.client is used in this paper. To access COM objects in the presented scripts, the Dispatch function is applied.

NumPy is a well-known module including objects and tools for scientific computing [20,57]. The functions and objects from the NumPy module work faster in numerical applications than standard Python functions and variables. The reason is that the main elements of NumPy were translated from Fortran and implemented in Python. This mechanism of method transfer may be continued by the user with the F2PY package [32] or ctypes module [33]. In the present research, only the array constructors available in NumPy were applied. These are array and empty methods. The string, integer, and float element arrays were used. However, extension of some presented examples, e.g., Examples 2 and 3, enables application of one's own methods and opens the way for more sophisticated techniques.

Another important package used here was SciPy [29,31]. The optimization part of SciPy was implemented in Example 2. The scipy.optimization module contains several optimization methods and root finding algorithms [31]. In this study, the most basic and well-known Nelder-Mead downhill simplex algorithm was used. A detailed description of it was provided by Press et al. [58].

The sediment simulation data in HEC-RAS are stored in XML format. It is a text-based markup language, derived from the Standard Generalized Markup Language. A detailed description of XML format and usage may be found in Reference [59]. To access data stored in an XML file, the ElementTree module is applied as described in References [35,60]. This package includes objects of the Element type. They are container objects designed to store hierarchical data structures, such as XML. The Element type is a cross between a list and a dictionary in Python. The module imported is xml.etree.ElementTree.

The results of the sediment simulation are stored in HDF5 format. Large sets of numerical data are easily accessed in HDF5 through the well-described hierarchy. A description of this format may be found in HDF Group [37]. To access these data the h5py module is applied [61].

### *2.4. General Remarks on Analyzed Examples*

Three examples were tested. The first was a code of a basic simulation composed in a way that enabled comparisons with the examples presented by Goodell [20,21]. The comparison revealed the most important differences between VBA and Python, such as initialization of HECRASController, access to controller methods, handling of data and results.

The second example was the calibration of the steady flow model. Although there is a procedure for the calibration of roughness coefficients in HEC-RAS, it works only with the unsteady flow model. However, calibration of the steady flow model is useful in many cases. Additionally, the approach presented in this example could also be applied in sensitivity analysis of the model performance or probabilistic flood risk mapping. The Nelder-Mead simplex algorithm was applied to search for the proper roughness coefficients. This is routine from the optimization part of the SciPy module, namely scipy.optimize. It works for similar purposes as the AHYDRA application [24], but the optimization algorithm is applied instead of manual alteration of roughness coefficients. The code is prepared in such a way that the applied Nelder-Mead method may be replaced with any optimization algorithm, e.g., the Genetic Algorithm as proposed by Leon and Goodell [25] with MATLAB.

The objective of the third example was to automatically change sediment samples in sediment data and run a sediment simulation. Such an approach has not been found by the author in the literature to date. The sediment sample data describe the percentage of sediment fractions in a sample. In HEC-RAS, they are called gradations. This function plays the role of a model parameter. It is necessary for the calculation of the total sediment transport, as well as transport of particular fractions. However, in practical applications it should be a piece of material sampled from the river bed. Collection of sediment samples is not very common. The measurements are local and extrapolated to a region, e.g., the cross-section. Hence, such measurements are highly uncertain, which is an additional reason for presentation of Example 3. This example is very simple, but may be extended for more sophisticated cases. It requires handling of XML and HDF5 data with special Python modules. Although, the control of sediment simulation is a more complex problem, the example is focused on change of the sediment samples, because these data are accessible in the most difficult way. If the sediment fractions may be changed by the Python script, all other elements of the sediment simulation, e.g., sediment transport formula, parameters of chosen formula, method for calculation fall velocity, etc., may be changed too. Hence, the adjustment of sediment sample in the XML file is a good illustrative example of sediment data management in general.

There are two main application areas for the mechanism presented in Example 3. The first is obviously the uncertainty or sensitivity analysis of sediment transport simulation. Considering the complexity of such computations, the analysis of sensitivity is too difficult to perform in the standard manual way. However, it is necessary if the sediment transport models are going to be used for forecast purposes, which sometimes happens today. The sediment, as well as flow data, are too uncertain to neglect this problem. The second application area is related to historical data and their processing. Sometimes data collected in the past include the full longitudinal profiles, as well as cross-sections, for different moments in time. The present author had at his disposal such data for Ner river collected for 1983, 2003, and 2007 [62]. These data were collected by the company BIPROWODMEL, for the analysis of regulation requirements. They illustrate the effects of the sedimentation process along the analyzed channel. Unfortunately, the company taking measurements did not collect sediment samples at the same time. Hence, the modeling of sediment transport in the Ner river is possible only if the model is calibrated with respect to these data. There is one important aspect which should be carefully considered. The calibration should also consider other factors, such as the sediment transport formula. However, the differences in the role and structure of the factors influencing the sediment transport simulation lead to the design of separate computational processes for each factor. One such process should be searching for an optimal sediment sample or samples, with other factors kept as constants. The third example presents the powerful capabilities of Python in this area.

An example with change of data stored in ASCII files of HEC-RAS is not presented here. Such examples have been discussed by other authors [20,25] with VBA and MATLAB. In the author's personal opinion, the handling of the ASCII files with Python is easier than with the other languages mentioned above. Hence, such an example would not bring anything new. The third and more difficult example, with XML and HDF files, may be treated as a help or guide in dealing with ASCII files. If Python enables reading and writing of information stored in such complex formats as XML and

HDF5, and it is quite easy to access the data in the ASCII files, it means that the presented scripting mechanisms could be applied to manipulate all data and results of any HEC-RAS model. This is also a huge advantage of the presented approach.

### *2.5. Applied Rules of Coding*

There are important differences between VBA and Python, which have to be discussed before the particular examples are presented. To explain them, the well-known examples given by Goodell [20,21] are quoted and compared with the devised Python code.

The first difference is declaration of variables. In Python, the variables are created dynamically during running of the program. Their type is recognized as the values that are assigned to them. The assignment operator may change this type, any time in the script. There are exceptions to these rules, e.g., special variables such as NumPy arrays. Their shape and type have to be declared before the variable is used. In VBA, dynamic creation is also available. However, it is more convenient to switch off this mechanism with the command Option Explicit at the beginning of the module. Then the variables have to be declared statically with the Dim statement.

To use the HECRASController object in the VBA code, it has to be declared as a variable [20,21]. In Python, the object is created in a similar way by calling the Dispatch function from the module win32com.client. Then all the functions and subroutines of the HEC-RAS controller are available in Python. The access is similar to the VBA calling of the HECRASController. The scheme is shown in Scripts 1–4.

In Script 1. HECobject is the HECRASController variable in the project, arguments is the list of input parameters and return\_values is the list of the values we would like to get from the subroutine. The HECobject has to initialized with the Dispatch function of the win32com.client module. Although it seems to be quite simple, there is one technical problem. The subprograms in VBA and Python have different structures and different mechanisms of data exchange with the main body of the program. There are two types of arguments used by VBA subroutines. These are the "called-by-value" and "called-by-references" arguments. Exchange of information between the subprogram and the main program based on "by-value" and "by-reference" mechanisms are well known. The arguments of the first type are used only as inputs. The second calling mechanism may work as input and output. When the function construction is used, another method for output is available, i.e., the returned value of the function. Although these mechanisms are very common and are applied in many programming languages, e.g., VBA, Pascal, Fortran, C/C++, etc., the subprograms used in Python are constructed in a different way. In general, the Python subprograms are functions. The parameters set in the function head are only inputs. The output is exchanged with the outer code only by returned values of the function. It is important to indicate that the function may return a list of values and each of them can be of different type.

```
import win32com.client 
HECobject = win32com.client.Dispatch("RAS503.HECRASController")