# Simulation of Cavitation Water Flows

Authors: Piroz Zamankhan
Journal: Mathematical Problems in Engineering
 (2015)

---

## Abstract

The air-water mixture from an artificially aerated spillway flowing down to a canyon may cause serious erosion and damage to both the spillway surface and the environment. The location of an aerator, its geometry, and the aeration flow rate are important factors in the design of an environmentally friendly high-energy spillway. In this work, an analysis of the problem based on physical and computational fluid dynamics (CFD) modeling is presented. The numerical modeling used was a large eddy simulation technique (LES) combined with a discrete element method. Three-dimensional simulations of a spillway were performed on a graphics processing unit (GPU). The result of this analysis in the form of design suggestions may help diminishing the hazards associated with cavitation.


---

## Body

## 1. Introduction

In spillway engineering there are numerous challenging obstacles. One of the most determining factors is the geological geometry in which the dam and its spillway have to be built in. The geometry, such as length and slope of the spillway, can range from short and steep to long and flat or vice versa. Generally, the longer and steeper is the spillway the higher is the gradient of the energy head of the flowing water. In some cases the velocity of the water might exceed 30 m/s and cavitation might occur. Used as an illustrative example, the high-energy spillway of the Kárahnjúkar Dam which supports the Hálsson Reservoir in eastern Iceland displays numerous problems in spillway engineering. An aerial perspective of this spillway can be seen in Figure1.Figure 1
Aerial view of Kárahnjúkar spillway with bottom aerator, eastern Iceland.When a fluid changes state from liquid to vapor, its volume increases by orders of magnitude and cavities (or bubbles) can form. Indeed, cavities, which are filled with water vapor, are formed by boiling at the local pressure which equals the water vapor pressure. If the cavity is filled with gases other than water vapor, then the process is named gaseous cavitation [1].There is a technical difference between boiling and cavitation. On the one hand, boiling is defined as the process of phase change from the liquid state to the vapor state by changing the temperature at constant pressure. On the other hand, cavitation is the process of phase change from the liquid state to the vapor state by changing the local pressure at constant temperature.Cavitation is a consequence of the reduction of the pressure to a critical value, due to a flowing liquid or in an acoustical field. Spillways are subject to cavitation. Inertial cavitation is the process in which voids or bubbles rapidly collapse, forming shock waves. This process is marked by intense noise. The shock waves formed by cavitation can significantly damage the spillway face [1]. Cavitation damage on the spillway face is a complex process [2]. Richer concrete mixes are used to increase the resistance to cavitation damage and erosion [3].Falvey [1] suggested that impurities and microscopic air bubbles in the water are necessary to initiate cavitation. As can be seen from Figure 1, the color of water accumulated behind the dam indicates that it may contain various types of impurities including suspended fine solid particles. The density of the water accumulated behind the dam was determined from a sample using a hydrometer. The relative density of the sample with respect to water has been found to be approximately 1.08.The melting water which transports suspended clay particles flows into the side channel and then follows an approximately 414 m long chute with a transition bend and a bottom aerator.The most important task is to reduce the energy head of the flowing water until it enters the river. In this example the jet should not hit the opposite canyon wall where it increases erosion. An artificial aeration device “cushions” the water; the widening at the chute’s end reduces the flow velocity and increases the impact area in the plunge pool which in turn decreases rock scouring. Additionally, seven baffles and lateral wedges at the main throw lip abet the jet disintegration. In general, every high-energy spillway has to be analyzed for cavitation risk.The water in spillways contains air bubbles and impurities, such as suspended particles, in a range of sizes. As discussed later these air bubbles and impurities induce cavitation. Vaporization is the most important factor in bubble growth. At a critical combination of flow velocity, water pressure, and the vapor pressure of the flowing water, cavitation starts. The cavitation number is used to define this starting point. The equation for this parameter is derived from the Bernoulli equation for a steady flow between two points. In dimensionless terms, the comparable equation results in a pressure coefficient,Cp. The value of this parameter is a constant at every point until the minimum pressure at a certain location is greater than the vapor pressure of water. The pressure at a certain location will not decrease any further once the vapor pressure is reached. Taking everything into account, such as temperature, clarity of the water, and safety reasons, the pressure coefficient is set to a minimum value and then called cavitation number σ. The cavitation number is a dimensionless number which expresses the relationship between the difference of a local absolute pressure from the vapor pressure and the kinetic energy per volume. It may be used to characterize the potential of the flow to cavitate.It has been suggested that no cavitation protection for a spillway would be needed if the cavitation number is larger than 1.8. If the cavitation number is in the range of 0.12–0.17, then the spillway should be protected by additional aeration grooves. In Figure2 the cavitation indices for different discharges have been calculated.Figure 2
Cavitation numbers at different discharges on the Kárahnjúkar spillway [4]. The location of the outlet is shown in Figure 1. Here, blue diamonds: 200 m3/s; red squares: 400 m3/s; green triangle: 800 m3/s; blue squares: 950 m3/s; yellow square with green star: 1350 m3/s; and orange dot: 2250 m3/s.As a precaution the cavitation zone has been set toσ=0.25 because of the lower atmospheric pressure, low water temperature, and overall safety reasons. Thus, as a precaution, the aerator was placed 125 m before the end of the chute. Nevertheless, incipient white water far upstream from the aerator allows the assumption that in such a long spillway fully self-aerated flow may have already been developed. Therefore, the purpose of this aerator is mainly to reduce the impact energy in the plunge pool.The main idea is to develop a fully 3D simulation of an aerated spillway in the future. But to achieve that the physics of bubble growth due to cavitation and aeration and the formation and distribution of the bubbles have to be understood and simulations have to be validated. In this work it is shown that LES combined with a discrete element method on GPUs yields promising results. Cavitation is defined as the formation of the vapor phase and the subsequent immediate implosion of small liquid-free zones in a liquid, called voids or bubbles [1].The first section deals with forced cavitation and the dynamics of cavitation bubbles. An experiment was carried out in which cavitation was forced to develop with clear water in glass and galvanized steel pipes. These findings were compared to computational simulations. The mathematical model and simulation results of the cavitation bubble dynamics will be described.In the second part of the paper the flow behavior of an aerated spillway is computed. In this section, the reliability of the predicted cavitation zone as shown in Figure2 is examined.


## 2. Experimental Methods and Observations

Cavitation occurs when the fluid pressure is lower than the local vapor pressure [5]. In the cavities, the vapor phase replaces the liquid phase. The surrounding liquid then experiences evaporative cooling. For water these thermal effects can have a significant influence on the dynamics of bubbles in the cavitation process. However, for simplicity, in this paper the isothermal problem will be investigated which retains a number of interesting qualitative features of the cavitation process.For water flow over a uniform roughness in a spillway, incipient cavitation will occur randomly over the spillway surface. If the cavitation number is lowered below incipient conditions, then cavitation will occur in sheets whose locations are influenced by the magnitude and spectrum of turbulent fluctuation within the boundary layer [1].To enhance the predictive reliability of computational models for cavitation, which are still considered to be immature, further validation based on experimental data is required. Figure3 shows a simple system in which sheet cavitation may be observed.Figure 3
The experimental setup for investigating cavitation. (a) Schematic of the experimental setup. Inset (on the right-hand side): an image of the aquarium stone used in the cavitation test. Inset (lower): the minimum pressure measured was 4.55 kPa (absolute). (b) An image of sheet cavitation at the dischargeQ = 12.5 L/s.
(a)
(b)Surface roughness, cracks, and offsets may initiate cavitation. The local pressure on the step’s surfaces or directly at the step in a high velocity spillway can reach the vapor pressure. This may result in sheet cavitation similar to that shown in Figure3. As stated earlier, cavitation can cause severe damage to the spillway concrete. The extent of cavitation damage is a function of the local cavitation number in the spillway chute and the duration of the flow.In brief, experiments were conducted in a test setup as illustrated in Figure3(a). The tests were performed in a horizontal galvanized steel pipe with an outer diameter of D=0.0254 m and the length of Lp = 0.3 m. Two very similar aquarium stones with a length of l= 0.0215 m and a height of h=0.01 m were attached to the bottom of the pipe in the entrance section. The stones produced a sudden offset whose leading edge is approximately an ellipse. The dimensions of one of the stones are depicted in the inset of Figure 3(a).Clear water was used in the experiments; its physical properties are given in Table1. The discharge Q=12.5 L/s was provided with a heavy-duty high-pressure wash-down pump with a pressure range of 2.5 bar to 140 bars (maximum pressure). A high-pressure hose was used to connect the test tube to the high-pressure pump.Table 1
Physical properties.
Material
Dynamic viscosity (Pa s)
Density (kg/m3)
Vaporization pressure (mmHg) [18]
Water  (liquid)
0.001
996.4
log
10
⁡
P
=
8.07131  -1730.63/T-39.724
Water (vapor)
1.34
×
1
0
-
5
0.5542
Air
1.79
×
1
0
-
5
1.175The average pressure at the pipe inlet was approximately 2.5 bars (absolute) and the minimum measured pressure was 4.55 kPa (absolute). Figure3(b) shows the monometer connected to the pressure measuring port (PMP) as shown in the second inset of Figure 3(a). The port was located Lm = 0.048 m downstream from the inlet of the test tube. Here, Lc represents the length of the sheet. The low pressure measured at the sampling port indicates that cavitation may have occurred in the test tube.The above-mentioned test with the same discharge ofQ = 12.5 L/s and the same average inlet pressure of 2.5 bar (absolute) was repeated in a cylindrical borosilicate glass pipe with a diameter of D=0.026 m. The pipe wall thickness was 0.001 m. A camera was used to take images of sheet cavitation from the underside of the test tube.Figure3(b) shows an image of sheet cavitation with an average inlet channel pressure of approximately 2.5 bars (absolute). As can be seen from Figure 3(b), the sheet appears to be attached to the stone with a length of Lc>5D. The sheet consists of a number of fuzzy white clouds. A flash photograph reveals that a cloud consists of individual small bubbles [1].In the following section, simulation results will be presented, giving a more detailed insight and showing an ongoing process of cavity formation and deformation.It is difficult to obtain accurate measurements using the simple setup as shown in Figure3. The stones attached to the bottom of the pipes may not be stable over long periods of time.Additional experimental studies in cavitation tunnels are required to supply a more reliable and convincing data basis for the validation of computational approaches to model cavitation.


## 3. Mathematical Model and Simulation Results

If sheet cavitation, as described in the experiments in the preceding section, occurs in the spillway chute, it can produce extensive damage to its concrete linings [1]. Note that cavitation damage in tunnels and conduits is less likely to result in dam failure.Experimental observations such as those shown in Figure3(b) are indispensable for investigating cavitation problems. However, numerical simulation can provide possibilities to study cavity formation and deformation in greater detail than what is affordable by experiments.


### 3.1. Bubble Dynamics in a Quiescent Fluid

The importance of microscopic bubbles as cavitation nuclei has been known for a long time [1]. Consider a single bubble, which contains the vapor of the liquid, in equilibrium in anincompressible, Newtonian fluid that is at rest at infinity. In this case, the bubble radius RE must satisfy the following condition [6]:(1)pv-p∞(0)RE3-2γRE2+G=0.The bubble radius at the limit point, where (pv-p∞(0))=(32γ3/27G)1/2, can be expressed in terms of material parameters as [5] (2)Rcritical=3G2γ1/2.If RE≤Rcritical, then the bubble isstableagainst infinitesimal changes in its radius. The bubble becomes unstable if RE>Rcritical. Note that the bubble may contain a “contaminant” gas which dissolves very slowly compared to the time scales associated with changes in its size.Apparently small bubbles whose radii are smaller than the critical radius behave like rigid spheres in accelerated flows [7].


### 3.2. Bubble Dynamics in Accelerated Flows

For the present study the interaction of a small bubble whose size is characterized by the critical bubble radius with its neighboring bubbles is assumed to be subjected to the Lennard-Jones (LJ) condition [8]. Figure 5(a) shows the Lennard-Jones 12-6 pair potential. The force between two bubbles with diameter dB located at ri and rj is given as(3)Fij=48εrij2dBrij12-12dBrij6rij.A simplified momentum equation for theith spherical bubble may be given as(4)π6ρvdB3dvBidt︸momentum transfer=π6ρwdB2CDReB,AcvBi-vwvBi-vw︸viscous+π12ρwdB3DvwDt-dvBidt︸added mass+32ρwdB2πvw1/2︸Basset×∫-∞tDvwDt-dvBidtdτt-τ1/2︸histroy term+π6ρwdB3DvwDt︸pressure term-π6dB3ρv-ρwg︸buoyancy+∑j=1NFij︸bubble-bubbleinteraction.In the current effort the aim is to combine the Lagrangian (bubble-based) model (4) and large eddy simulation (LES) in order to achieve more accurate simulations of cavitation water flows. The filtered continuity momentum for an isothermal 3D flow of water may be given as(5)∂ρ¯∂t+∇·ρ¯vw~=-Sv,(6)∂ρ¯vw~∂t+∇·ρ¯vw~vw~=-1ρw∇ρ¯p¯+1ρw∇·ρ¯σ~+∇·τ+fbw¯+ρ¯g.The viscous stress tensor in (6) is defined as(7)σ=-23μw∇·vwI+μw∇vw+∇vwT.In (6), τ may be given as(8)τ=-ρ¯vw~vw~~-vw~~vw~~-ρ¯vw′′vw~~+vw~vw′′~-ρ¯vw′′vw′′~.The SGS stress tensor, τ, is required to close the equations for the large-scale fields on a grid small enough (but much larger than the Kolmogorov scale) to provide reasonable resolutions. In contrast to the filtered single-phase equations, a conceptual restriction arising in the present approach is that the filter width, Δ, should strictly be larger than the length scale characteristic of the small bubbles. Thus, an appropriate choice of Δ should provide a sufficiently large-scale resolution while not violating the aforementioned restriction.Following Sauer and Schnerr [9], the mass flow rate through the surface of a bubble with the radius of R, Sv, may be given as(9)Sv=4πnBR21+43πnBR3-2ρwρvρR˙.Here, the changes in the bubble radius may be calculated using the famous Rayleigh-Plesset equation [6](10)pv+GT/R3-p∞ρw-1ρw2γR+4μwRR˙=RR¨+32R˙2.The vapor transport equation may be given as(11)∂∂tϕvρv+∇·ϕvρvVB=Sv.It is known that the presence of bubbles contributes to the processes of energy removal from the resolved scales of the liquid phase. This two-way coupling effect may be modeled by superposing bubble-induced SGS energy dissipation to that induced by shear. A tentative first attempt at closure, the momentum equation (6), may be closed with an SGS model for, τ, given as [10] (12)τmod⁡=-ρ¯vw~vw~~-vw~~vw~~+ρ¯2CsΔ2S~∶S~1/2+CsΔvslip~×S~-13tr⁡S~I.Here, the value of the backscatter parameter, Cb0, is set to Cb0=0.2. The model for τ has been validated for dense gas-particle flows [11–13].During each time step of the simulation, the forces acting on each particle should be calculated. This requires knowledge of the local values of the fluid velocity components at the position of the bubbles. These variables are only known at spatial grid points in the computational domain. A tricubic interpolation has been used for calculating the fluid velocity components at the centers of the bubbles. The bubbles are propagated using the generalized Verlet algorithm, whose parallelization for the GPU is described in detail in [10].


### 3.3. Simulation Results

#### 3.3.1. Air Bubbles

Figures4(b) and 4(c) show a computed large bubble in an aerated tank [7]. The bubble is a cluster of small bubbles with a size of the critical air bubble radius. Figure 4(d) depicts the denoised version of the bubble using multiscale image denoising based on a multiscale representation of the images [14]. The denoised version of the bubble-like structure in Figure 4(b) is very similar to the small bubbles observed in an aerated tank [7].Figure 4
The interaction of small bubbles in an aerated tank. (a) The Lennard-Jones 12-6 pair potential. (b) Computed bubble-like structure in an aerated tank at the Reynolds number of 29500. (c) A different view of (b). (d) The denoised version of (c). The length of the bubble is approximately 1 cm and its height is 0.65 cm. (e) Configuration of a three bubble-like structure in an aerated tank. (f) The bubbles in (e) coalesce to form a larger bubble after 0.3 sec. (g) The denoised version of (f).Figure 5
A CFD model for a stone in Figure3. (a) Taking 2D images of the stone by moving camera’s position at regular intervals. (b) Automatically aligning the images with the others. (c) Reconstruction of a 3D model from 2D image.
(a)
(b)
(c)Figures4(e) and 4(f) are snapshots of bubbles separated by t=0.3 sec. As can be seen from Figure 4(f), three bubbles in Figure 4(e) coalesce to form one large bubble. The denoised version depicted in Figure 4(g) is very similar to the bubbles with complex geometries shown in [7].To further assess the quality of the model described in the preceding section, the cavitation water flow in Figure3 is simulated in the following section.


#### 3.3.2. 3D Model Using 2D Images

Apparently, it is not possible to reconstruct a 3D model from a single image of the stone as shown in the inset of Figure5(a). To get around this limitation, an optical processing algorithm can be developed that employs multiple photographs taken at every 15 degrees, from the above and from the side to create a 3D model of the stone. Figure 5(a) illustrates the positions of the camera for taking about 40 photographs of the stone. The images acquired are the input of the Autodesk 123D Catch 3D scanning software [15].Figure5(b) shows the basics of the stitching process for the accurate creation of a 3D model of the stone. Figure 5(c) represents a medium quality 3D mesh of the stone which is used to simulate sheet cavitation in the following section.


#### 3.3.3. Sheet Cavitation

Aspects of the Simulation. As mentioned earlier, turbulence modeling has a critical role in cavitation prediction. Cavitation flow in the liquid-vapor region is locally compressible [9, 16]. To capture the shedding dynamics and the unsteadiness of cavitation the modified form of LES as detailed in the preceding section would be required. Note that (5) and (6) look like the equations of motion of a fluid with variable density ρ=ρw1-ϕv.In this section, the mathematical model described in the preceding section will be used to analyze sheet cavitation in a pipe flow as shown in Figure3. The discharge was set to Q=12.5 L/s. The tube inner diameter was D=0.024 m and its length was L=0.3 m.The length and the height of the stones in the simulations are exactly the same as those used in the experiments. The model of the stone is shown in Figure5(c). It is likely that cavitation as shown in Figure 3(b) starts at minute cracks on the otherwise smooth surface of the stones. The aforementioned system may be characterized as belonging to the singular roughness category [1].In brief, the Lagrangian (bubble-based) and the Eulerian (grid-based) methods are used to simulate the cavitation singular roughness in a tube flow. Figure6(a) illustrates the model used in the simulations. Here, the dimensions are Lp = 0.3 m, D=0.024 cm, and hs = 0.004 m. Figure 6(b) depicts a top view of the entrance region of the pipe in which the positions of the stones are clearly shown. The total length of the pipe covered by stones is l2 = 0.037 m, and the pipe’s length covered by the first stone is l1 = 0.0185 m. The length and height of the stones are exactly the same as those shown in Figure 3(a). Figure 6(c) shows the grid which consists of more than 5×106 tetrahedral meshes.Figure 6
A CFD model for the tube in Figure3. (a) Schematic of the pipe used in sheet cavitation simulations with some nomenclatures. (b) A top view of the entrance region of the pipe. (c) The grid used in the cavitation simulations. Insets: the grid is magnified and replotted.
(a)
(b)
(c)The filtered Navier-Stokes equations (5) and (6) are discretized with the finite-volume method on a staggered Cartesian grid. Convective and diffusive fluxes are approximated with central differences of second-order accuracy and time advancement is achieved by a second-order, explicit Adams-Bashforth scheme. The equation for coupling the pressure to the velocity field is solved iteratively using the SIP method [17]. The subgrid-scale stresses appearing in the filtered Navier-Stokes equations are computed using (12). No-slip boundary conditions are employed on the surface of the stone and the tube walls. The physical properties used in the simulations are listed in Table 1. The temperature was assumed to be constant and it was set to T=302 K. Note that the vaporization pressure [18] at that temperature is 4500 Pa (absolute). The acceleration due to gravity is g=9.81 m/s2.The Lagrangian (bubble-based) simulation was performed to calculate the motion of bubbles in the liquid using (4). In brief, equivalent spherical bubbles were determined from irregular vapor structure by marking the cells with ϕw<0.95. Equation (4) was integrated using the generalized Verlet algorithm. The time step used in the simulation was Δt = 5×10-7 s. The linear interpolation routines were used to communicate the information from grid nodes to particle positions and vice versa. The simulations were performed by using GPU computing [10].Results. Obviously, a single parameter such as the cavitation number cannot describe many of the complexities of cavitation. Falvey [1] suggested that for flow past stones as shown in Figure 7(a), cavitation will not occur if the cavitation number is greater than about 1.8. Figure 7(a) shows that cavitation bubbles form within the flow at the discharge Q=12.5 L/s. This indicates that the cavitation number was below 1.8. Figure 7(b) represents a view from below of an instantaneous configuration of bubbles in the cavitation water flow. In this case, a sheet consists of a large number of small bubbles which are extensively developed downstream of the step formed by the stones in the pipe. Here, cavitation is formed by turbulence in the shear zone which is produced by the sudden change of flow direction at the stone face. Figure 7(b) also shows large nonspherical cavities with patterns similar to those shown in Figure 3(b).Figure 7
Computer experiments at a cavitation number ofσh=0.9 in a duct. (a) A top view of the instantaneous configuration of bubbles in a cavitating water flow of water through a pipe with inner diameter of D=0.024 m at the discharge of Q=12.5 L/s and at the cavitation number of σh=0.9. (b) A view from below of the bubbles in the cavitating flow of water through the pipe. Here, the large cavities consisting of small individual bubbles are illustrated. (c) Computed contours of the averaged vapor volume fraction on thexz-plane for the cavitating water flow in (a). Here, a color code as illustrated below the stones is used to indicate the values of the averaged vapor volume fraction on the surfaces of the stones. (d)–(g) Variations of ϕ-v as a function of lw at the sampling ports 1 through 16, respectively.
(a)
(b)
(c)
(d)
(e)
(f)
(g)Following Falvey [1], it may be estimated that the system in Figure 7(a) operates at a cavitation number of σh=0.9. For this case, the reference location is immediately upstream of the first stone and at its maximum height. The computed length of the sheet is quite comparable with that shown in Figure 3(b). It would be expected that for even lower values of the cavitation number, the clouds form one long supercavitating pocket.Figure7(c) depicts the computed instantaneous contours of the averaged vapor volume fraction in thexz-plane. In this figure, the locations at which cavitation occurred on the surfaces of the stones are presented by the use of color coding. As can be seen from Figure 7(c), the process starts with the occurrence of swarm of bubbles in the small regions of the first stone at sampling port 3 located at x=0.015 m. However, no evidence of cavitation can be found at sampling port 4 located at x=0.022 m in a narrow cleft between the stones where the local pressure is much higher than that of the vapor pressure. This figure also shows that a large number of small bubbles are extensively produced on the trailing edge of the second stone close to sampling port 8 located at x=0.045 m. Sampling ports 1 through 16 are all located in the xz-plane. The distances of these ports from the inlet of the pipe are x=5×10-5, 0.0077, 0.015, 0.022, 0.025, 0.029, 0.035, 0.045, 0.055, 0.065, 0.09, 0.12, 0.2, 0.25, 0.29, and 0.298 m, respectively.Figure7(d) illustrates the averaged vapor volume fraction, ϕ-v, as a function of the vertical distance from the pipe wall, lw, at sampling ports 1 through 4. Here, the squares, circles, diamonds, and left triangles represent the computed vapor volume fraction at ports 1, 2, 3, and 4, respectively. This figure indicates the occurrence of a swarm of bubbles at sampling port 3.The squares, circles, diamonds, and left triangles in Figure7(e) represent variations of ϕ-v as a function of lw at sampling ports 5, 6, 7, and 8, respectively. This figure indicates that cavitation reoccurred at sampling port 7 and became extensive at port 8. These complexities cannot be described with a single parameter such as cavitation number.Figure7(f) shows variations of ϕ-v as a function of lw at sampling ports 9, 10, 11, and 12 using the squares, circles, diamonds, and left triangles, respectively. As can be seen from Figure 7(f), the vapor bubbles roll up into a larger volume and become cloudy as they are transported further downstream and leave the entrance region of the pipe.Figure7(g) represents the final stage of the process which is an ongoing process of bubble formation and deformation. This figure depicts variations of ϕ-v as a function of lw at sampling ports 13, 14, 15, and 16 using the squares, circles, diamonds, and left triangles, respectively.In brief, the appearance of visible cavitation in the flowing water in the pipe as shown in Figure7 was preceded by the occurrence of small bubbles in the small area on the surface of the stones. This observation highlights the importance of bubbles as cavitation nuclei which has been known for a long time [6].Figure7(d) indicated that a swarm of small bubbles occurred at sampling port 3. However, no evidence of cavitation was found at sampling port 4 which was located further downstream from port 3. To address this issue, the pressure field is illustrated in thexz-plane in Figure 8(a). Here, the dimensionless pressure Π* is defined as Ln(p)/Ln(pv), where p represents the local pressure. Figure 8(b) depicts variations of the dimensionless pressure p* (defined as p/pv) as a function of lw at sampling ports 1, 2, 3, and 4 using the squares, circles, diamonds, and left triangles, respectively. This figure indicates that the local pressure in the narrow cleft between the stones is approximately 50 times higher than that of the vapor pressure.Figure 8
Cavitation in a tube. (a) The instantaneous computed contours of the dimensionless pressureΠ* on thexz-plane bubbles in a cavitation flow of water at the cavitation number of σh=0.9. The locations of the sampling ports are the same as those in Figure 7(c). (b)–(e) Variations of the dimensionless pressure p* as a function of lw. (f) The instantaneous computed velocity field of a cavitating water flow through the pipe with inner diameter of D=0.024 m at the discharge of Q=12.5 L/s.
(a)
(b)
(c)
(d)
(e)
(f)Figure8(c) illustrates variations of the dimensionless pressure, p*, as a function of lw at sampling ports 5, 6, 7, and 8 using squares, circles, diamonds, and left triangles, respectively. As can be seen from this figure, the local pressure at port 7 on the surface of the second stone reached the vapor pressure. Consequently, the small bubbles were formed on the trailing edge of the second stone as shown in Figure 7(e). The appearance of visible cavitation as shown in Figure 7(a) was preceded by the occurrence of these bubbles in the small area around sampling port 5 on the surface of the second stone. The computed low pressure at port 8 (which is the same as that measured at the PMP in Figure 3(a)) provides a favorable situation for a sheet, which consists of a large number of small bubbles, to be developed and transported downstream as depicted in Figures 7(a) and 7(b). The low pressure condition extends from sampling port 9 to sampling port 14 as shown in Figures 8(d) and 8(e). The squares, circles, diamonds, and left triangles in Figure 8(d) represent variations of p* as a function of lw at sampling ports 9, 10, 11, and 12, respectively.In the aforementioned region the vapor bubbles roll up into a larger volume and become less extensive as they are transported further downstream. The volume fraction of the vapor bubble is lower at the upper part of the pipe between ports 9 and 14 due to the fact that the pressure slightly increases withlw.The squares, circles, diamonds, and left triangles in Figure8(e) represent variations of p* as a function of lw at sampling ports 13, 14, 15, and 16, respectively. This figure indicates that the pressure increased from ports 14 to 16 and reached the atmospheric pressure at the pipe outlet at which the vapor bubbles vanished.The area-weighted average pressure at the inlet of the pipe was found to be 257 kPa. The results show agreement between simulation and experiment. The model is therefore deemed to be sufficiently flexible to capture a number of the interesting qualitative features of the cavitation process.Figure8(f) shows the computed velocity vector filed in thexz-plane. This figure indicates a complex shearing flow around sampling port 8 at which cavitation was formed. As stated earlier, the system in Figure 7(a) operates at a cavitation number of σh=0.9. Note that the computed pressure at sampling port 1 immediately upstream of the first stone and at its maximum height is approximately 300 kPa. In this case, the computed average water velocity is 25 m/s. The area-weighted average pressure at the inlet of the pipe was found to be 257 kPa.In the following section, the model is employed to check the reliability of the predicted cavitation zone as shown in Figure2.


#### 3.3.4. Cavitation in a Spillway Chute

The previous section described the complexity of the occurrence of cavitation. As stated earlier, Figure2 indicates that the cavitation number should be quite low at approximately 135 m from the end of the chute of Kárahnjúkar spillway. Apparently, the designers of the Kárahnjúkar spillway decided to use an aerator at that location (as shown in Figure 1) in order to protect the spillway surface from cavitation.Recall that Semenkov and Saranchev [19] suggested that the beginning of the spillway is the cavitation-hazardous section. Particularly for large discharges, natural air entrainment is absent in more than 60 m from the inlet [19]. The bottom turbulent layer is shielded above by the undisturbed potential flow core, which reduces diffusion air entrainment [19]. Consequently, on the entire length of the cavitation-hazardous section of the spillway passage the volume air concentration in the boundary layer with thickness 0.2–0.25 m would be less than 7-8%. In this case, artificial air entrainment is needed to protect the surface of the spillway [20].In this section, the usefulness of using artificial air at 125 m from the end of the chute of the Kárahnjúkar spillway is investigated using the mathematical model as detailed in the preceding section. To this end, a simplified model of Kárahnjúkar spillway with a hypothetical hole with diameter of 0.02 m is used in the simulation. Figure9(a) shows the hypothetical hole which is essential to initiate cavitation at the cavitation number of 0.25. The simplified model of Kárahnjúkar spillway is 414 m long with a slope of approximately α = 11°.Figure 9
Cavitation in a spillway chute. (a) Schematic of a hypothetical hole with a diameterdh=0.02 on the surface of a model of Kárahnjúkar spillway at a distance from the inlet of Lh=300 m. The total length of the spillway is LT=414 m and its slope is α = 11°. Inset: the computed contours of velocity in the spillway. The distances of the sampling ports 1, 2, 3, and 4 from the inlet are 0, 107, 300, and 414 m, respectively. (b) The averaged air volume fraction as a function of lw at the sampling ports 1, 2, 3, and 4. Here, the squares, circles, diamonds, and left triangles represent the air volume fraction at the sampling ports 1 through 4, respectively. (c) The computed velocity vector field around the hole at the sampling port 3. Here, a color code as illustrated on the left-hand side of the hole is used to indicate the values of the local pressure.
(a)
(b)
(c)The inlet boundary condition atx=0 was set to V=20 m/s. The height of the water at the inlet was 3 m. Furthermore, it was assumed that the top surface of the model is open to the atmosphere. The outlet boundary condition was set to p=patm. Periodic boundary conditions were used in the z-direction. The no-slip boundary condition was applied on the walls of the model. To capture the effects of viscous boundary layers, a wall-function was used to specify velocity at the forcing points as [21](13)uu*=y+y+≤5-1.15576+1.47869y+-0.0527848y+2+0.000655444y+35<y+≤30ln⁡y+Κ+5.1y+>30,where y+ represents dimensionless, sublayer-scaled distance, and Κ is the von Kármán constant (typically the value 0.41 is used).The inset of Figure9(a) represents the computed velocity contours in the model of the Kárahnjúkar spillway. The computed bubble volume fraction as a function of vertical distance from the wall lw at four sampling ports is illustrated in Figure 9(b). In this case, the bubble is an air bubble which may contain vapor. As can be seen from Figure 9(a), the hypothetical hole is located at 135 m from the outlet of the model. Figure 9(b) indicates that the thickness of the air-water mixture increases with x. Here, the squares, circles, diamonds, and left triangles represent the bubble volume fraction at sampling ports 1, 2, 3, and 4, respectively. As can be seen from Figure 9(b), the bubble volume fraction at sampling port 2 is lower than 7%. This finding highlights the need of artificial aerators between sampling ports 1 and 2 in order to protect the surface of the spillway.Figure9(b) indicates that the bubble volume fraction at the hypothetical hole (i.e., at sampling port 3) is approximately 38% which is in accordance with Semenkov and Saranchev [19] and very much above the “cavitation safe” value of 7% of entrained air. In this case, the vapor volume fraction is quite low.Figure9(c) illustrates the computed velocity vector field. A color code as illustrated on the left hand side of the hole in Figure 9(c) is used to indicate the values of the local pressure. This figure indicates that the water flows mainly over the hole. A small part of the water stream hits the hole’s face and reticulates inside the hole. As can be seen from Figure 9(c), the pressure reaches a maximum of 2500 kPa and diverts the flow to an upward direction over a little edge and in a downward direction which initiates the water swirl in the hole. The local pressure in the center of the hole is slightly lower than 5 kPa. This finding indicates that cavitation might occur at 125 m from the end of the chute of the Kárahnjúkar spillway. Indeed, cavitation will occur for the flow past the sudden into-the-flow offset introduced by the hypothetical hole as noted in Figure 9(c). This figure indicates a low-pressure region and thus potential cavitation might also occur shortly after passing the hole’s edge. However, the air volume fraction at the hole as shown in Figure 9(b) is approximately 38%. This finding suggests that the use of an aerator at port 3 in order to protect the surface of the spillway might not be needed. Adding extra air at sampling port 3 would produce a very thick air-water mist flow at the outlet of the spillway. This will be further analyzed in the following section.Figure9(b) also indicates that in the absence of an artificial aerator the computed thickness of the air-water mixture at the outlet of the spillway is approximately 8 m. In addition, the average bubble volume fraction is higher than 50%.


#### 3.3.5. Artificial Aeration

The designers of the Kárahnjúkar spillway were challenged in their task due to the lack of literature on air entrainment for spillways. As mentioned in the preceding section, an aerator must be introduced to prevent cavitation in the Kárahnjúkar spillway at high discharges. The air volume fraction in the spillway downstream from sampling port 2 in Figure9(b) is well above the 7% minimum requirement of its protection against cavitation. Hence, artificial aeration would only be needed between sampling ports 1 and 2 illustrated in Figure 9(a).In this section, the mathematical model described in the preceding section will be used to analyze the Kárahnjúkar spillway with a bottom aerator. Figure10 illustrates a CFD model which can be imagined as a 1 : 15 Froude scale of the Kárahnjúkar bottom aerator. Figure 1 shows an aerial view of the Kárahnjúkar spillway and its bottom aerator located in eastern Iceland.Figure 10
A CFD model for the Kárahnjúkar bottom aerator. (a) Perspective view of the simulation model and selected nomenclature. (b) The front view of the spillway with bottom aerator. Here,δ = 0.02 m, L = 2.42 m, Ha = 0.44 m, Hw = 0.11 m, L1 = 0.11 m, L2 = 2.15 m, h0 = 0.022, Hout = 0.55 m, θ1 = 20°, and θ2=30°.
(a)
(b)By implementing periodic boundary conditions in they-direction, a small part of the system that is far from the vertical side walls was simulated. The model flow rate was set to vw0 = 20 m/s at the inlet with a size of Hw = 0.11 m, corresponding to Q= 2250 m3/s.As can be seen from Figure11, the water stream is deflected by a deflector in order to ventilate the nappe to the atmosphere. The connection to the atmosphere is made via a channel under the ski jump of the aerator. Note that the air entraining capacity of the aerator depends on the take-off distance from the free flow. The air pressure is atmospheric at the inlets and the outlet.Figure 11
Computed contours of the velocity field for the bottom aerator of Kárahnjúkar spillway. Inset: a pair of vortices in the ventilation zone.Figure12(a) and its insets illustrate the computed contours and vector velocity field of a more optimized design for the aerator of the Kárahnjúkar spillway. Here, vortices can be seen to be created at the ventilation gate and to detach periodically from its other side. A pair of vortices is magnified and replotted in the inset. The air enters the ventilation zone with a normalized velocity of v*=3.81. Here, the velocity is normalized with the water velocity at inlet vw0. Figure 12(b) represents an instantaneous configuration of bubbles in the spillway with an optimized bottom aerator. The high stream of air tends to move the water stream up. In this case, the air concentration in the ventilated zone appears to be more than that needed to prevent cavitation. Figure 1 shows a highly aerated water stream which is likely to serve as a source for an intense canyon erosion.Figure 12
The computed flow field for the aerator of Kárahnjúkar spillway. (a) Computed contours of velocity field of a more optimized design for the aerator of Kárahnjúkar spillway. Inset: vortical structures in the ventilation zone. (b) Instantaneous configuration of bubbles in a spillway with ski jump aerator. Here,Lnb = 0.62 m. (c) Air concentration distribution versus z in the aeration region. Here, right triangles, squares, left triangles, and diamonds represent the air concentration at the sampling ports S1, S2, S3, and S4, respectively.
(a)
(b)
(c)When water flows over a spillway, vortex sheet instability at the air-water interface known as Kelvin-Helmholtz instability can occur. The instability will be in the form of waves being generated on the water surface which can initiate natural air entrainment.The air entrainment process is strengthened when a growing boundary layer in the water flow reaches the free surface. In this case, the entrained air from the free surface can reach the spillway surface to minimize the danger of cavitation. As can be seen from Figure12, the details of air entrainment through the upper air-water interface of the jet can be predicted by the model. In addition, the details of air entrainment through the lower interface are very similar to those presented in [10]. Here, Lnb represents the length at which the air entrainment through the lower interface reaches that through the upper interface. Figure 12(c) depicts the variations of air concentration as a function of z at four different sampling ports as shown in (b).The region at which the natural air entrainment effectively eliminates the risk of damage due to cavitation is called the fully aerated region. The artificial aeration process discussed above would not be needed in the fully aerated region.Figure13(a) illustrates the computed contours of the velocity field for a hypothetical case in which the discharge of the Kárahnjúkar spillway is Q= 3150 m3/s. Figure 13(b) indicates that the air enters the ventilation zone with a normalized velocity of v*=5.6. The vortex shedding of the aeration gate is quite intense. Therefore, the air-water mixture entering the canyon could cause extensive erosion. Recall that the canyon shown in Figure 1 suffered a disastrous erosion in 2009.Figure 13
The computed flow field for a hypothetical case with a discharge ofQ=3150 m3/s. (a) Computed contours of velocity field. (b) Computed vector plot of vertical structures in the ventilation zone. (c) and (d) Variations of the normalized velocity v* as a function of z. Squares, circles, diamonds, left triangles, gradients, and deltas represent the normalized velocity at the sampling ports S1, S6, S5, S2, S3, and S4, respectively.
(a)
(b)
(c)
(d)


## 3.1. Bubble Dynamics in a Quiescent Fluid

The importance of microscopic bubbles as cavitation nuclei has been known for a long time [1]. Consider a single bubble, which contains the vapor of the liquid, in equilibrium in anincompressible, Newtonian fluid that is at rest at infinity. In this case, the bubble radius RE must satisfy the following condition [6]:(1)pv-p∞(0)RE3-2γRE2+G=0.The bubble radius at the limit point, where (pv-p∞(0))=(32γ3/27G)1/2, can be expressed in terms of material parameters as [5] (2)Rcritical=3G2γ1/2.If RE≤Rcritical, then the bubble isstableagainst infinitesimal changes in its radius. The bubble becomes unstable if RE>Rcritical. Note that the bubble may contain a “contaminant” gas which dissolves very slowly compared to the time scales associated with changes in its size.Apparently small bubbles whose radii are smaller than the critical radius behave like rigid spheres in accelerated flows [7].


## 3.2. Bubble Dynamics in Accelerated Flows

For the present study the interaction of a small bubble whose size is characterized by the critical bubble radius with its neighboring bubbles is assumed to be subjected to the Lennard-Jones (LJ) condition [8]. Figure 5(a) shows the Lennard-Jones 12-6 pair potential. The force between two bubbles with diameter dB located at ri and rj is given as(3)Fij=48εrij2dBrij12-12dBrij6rij.A simplified momentum equation for theith spherical bubble may be given as(4)π6ρvdB3dvBidt︸momentum transfer=π6ρwdB2CDReB,AcvBi-vwvBi-vw︸viscous+π12ρwdB3DvwDt-dvBidt︸added mass+32ρwdB2πvw1/2︸Basset×∫-∞tDvwDt-dvBidtdτt-τ1/2︸histroy term+π6ρwdB3DvwDt︸pressure term-π6dB3ρv-ρwg︸buoyancy+∑j=1NFij︸bubble-bubbleinteraction.In the current effort the aim is to combine the Lagrangian (bubble-based) model (4) and large eddy simulation (LES) in order to achieve more accurate simulations of cavitation water flows. The filtered continuity momentum for an isothermal 3D flow of water may be given as(5)∂ρ¯∂t+∇·ρ¯vw~=-Sv,(6)∂ρ¯vw~∂t+∇·ρ¯vw~vw~=-1ρw∇ρ¯p¯+1ρw∇·ρ¯σ~+∇·τ+fbw¯+ρ¯g.The viscous stress tensor in (6) is defined as(7)σ=-23μw∇·vwI+μw∇vw+∇vwT.In (6), τ may be given as(8)τ=-ρ¯vw~vw~~-vw~~vw~~-ρ¯vw′′vw~~+vw~vw′′~-ρ¯vw′′vw′′~.The SGS stress tensor, τ, is required to close the equations for the large-scale fields on a grid small enough (but much larger than the Kolmogorov scale) to provide reasonable resolutions. In contrast to the filtered single-phase equations, a conceptual restriction arising in the present approach is that the filter width, Δ, should strictly be larger than the length scale characteristic of the small bubbles. Thus, an appropriate choice of Δ should provide a sufficiently large-scale resolution while not violating the aforementioned restriction.Following Sauer and Schnerr [9], the mass flow rate through the surface of a bubble with the radius of R, Sv, may be given as(9)Sv=4πnBR21+43πnBR3-2ρwρvρR˙.Here, the changes in the bubble radius may be calculated using the famous Rayleigh-Plesset equation [6](10)pv+GT/R3-p∞ρw-1ρw2γR+4μwRR˙=RR¨+32R˙2.The vapor transport equation may be given as(11)∂∂tϕvρv+∇·ϕvρvVB=Sv.It is known that the presence of bubbles contributes to the processes of energy removal from the resolved scales of the liquid phase. This two-way coupling effect may be modeled by superposing bubble-induced SGS energy dissipation to that induced by shear. A tentative first attempt at closure, the momentum equation (6), may be closed with an SGS model for, τ, given as [10] (12)τmod⁡=-ρ¯vw~vw~~-vw~~vw~~+ρ¯2CsΔ2S~∶S~1/2+CsΔvslip~×S~-13tr⁡S~I.Here, the value of the backscatter parameter, Cb0, is set to Cb0=0.2. The model for τ has been validated for dense gas-particle flows [11–13].During each time step of the simulation, the forces acting on each particle should be calculated. This requires knowledge of the local values of the fluid velocity components at the position of the bubbles. These variables are only known at spatial grid points in the computational domain. A tricubic interpolation has been used for calculating the fluid velocity components at the centers of the bubbles. The bubbles are propagated using the generalized Verlet algorithm, whose parallelization for the GPU is described in detail in [10].


## 3.3. Simulation Results

### 3.3.1. Air Bubbles

Figures4(b) and 4(c) show a computed large bubble in an aerated tank [7]. The bubble is a cluster of small bubbles with a size of the critical air bubble radius. Figure 4(d) depicts the denoised version of the bubble using multiscale image denoising based on a multiscale representation of the images [14]. The denoised version of the bubble-like structure in Figure 4(b) is very similar to the small bubbles observed in an aerated tank [7].Figure 4
The interaction of small bubbles in an aerated tank. (a) The Lennard-Jones 12-6 pair potential. (b) Computed bubble-like structure in an aerated tank at the Reynolds number of 29500. (c) A different view of (b). (d) The denoised version of (c). The length of the bubble is approximately 1 cm and its height is 0.65 cm. (e) Configuration of a three bubble-like structure in an aerated tank. (f) The bubbles in (e) coalesce to form a larger bubble after 0.3 sec. (g) The denoised version of (f).Figure 5
A CFD model for a stone in Figure3. (a) Taking 2D images of the stone by moving camera’s position at regular intervals. (b) Automatically aligning the images with the others. (c) Reconstruction of a 3D model from 2D image.
(a)
(b)
(c)Figures4(e) and 4(f) are snapshots of bubbles separated by t=0.3 sec. As can be seen from Figure 4(f), three bubbles in Figure 4(e) coalesce to form one large bubble. The denoised version depicted in Figure 4(g) is very similar to the bubbles with complex geometries shown in [7].To further assess the quality of the model described in the preceding section, the cavitation water flow in Figure3 is simulated in the following section.


### 3.3.2. 3D Model Using 2D Images

Apparently, it is not possible to reconstruct a 3D model from a single image of the stone as shown in the inset of Figure5(a). To get around this limitation, an optical processing algorithm can be developed that employs multiple photographs taken at every 15 degrees, from the above and from the side to create a 3D model of the stone. Figure 5(a) illustrates the positions of the camera for taking about 40 photographs of the stone. The images acquired are the input of the Autodesk 123D Catch 3D scanning software [15].Figure5(b) shows the basics of the stitching process for the accurate creation of a 3D model of the stone. Figure 5(c) represents a medium quality 3D mesh of the stone which is used to simulate sheet cavitation in the following section.


### 3.3.3. Sheet Cavitation

Aspects of the Simulation. As mentioned earlier, turbulence modeling has a critical role in cavitation prediction. Cavitation flow in the liquid-vapor region is locally compressible [9, 16]. To capture the shedding dynamics and the unsteadiness of cavitation the modified form of LES as detailed in the preceding section would be required. Note that (5) and (6) look like the equations of motion of a fluid with variable density ρ=ρw1-ϕv.In this section, the mathematical model described in the preceding section will be used to analyze sheet cavitation in a pipe flow as shown in Figure3. The discharge was set to Q=12.5 L/s. The tube inner diameter was D=0.024 m and its length was L=0.3 m.The length and the height of the stones in the simulations are exactly the same as those used in the experiments. The model of the stone is shown in Figure5(c). It is likely that cavitation as shown in Figure 3(b) starts at minute cracks on the otherwise smooth surface of the stones. The aforementioned system may be characterized as belonging to the singular roughness category [1].In brief, the Lagrangian (bubble-based) and the Eulerian (grid-based) methods are used to simulate the cavitation singular roughness in a tube flow. Figure6(a) illustrates the model used in the simulations. Here, the dimensions are Lp = 0.3 m, D=0.024 cm, and hs = 0.004 m. Figure 6(b) depicts a top view of the entrance region of the pipe in which the positions of the stones are clearly shown. The total length of the pipe covered by stones is l2 = 0.037 m, and the pipe’s length covered by the first stone is l1 = 0.0185 m. The length and height of the stones are exactly the same as those shown in Figure 3(a). Figure 6(c) shows the grid which consists of more than 5×106 tetrahedral meshes.Figure 6
A CFD model for the tube in Figure3. (a) Schematic of the pipe used in sheet cavitation simulations with some nomenclatures. (b) A top view of the entrance region of the pipe. (c) The grid used in the cavitation simulations. Insets: the grid is magnified and replotted.
(a)
(b)
(c)The filtered Navier-Stokes equations (5) and (6) are discretized with the finite-volume method on a staggered Cartesian grid. Convective and diffusive fluxes are approximated with central differences of second-order accuracy and time advancement is achieved by a second-order, explicit Adams-Bashforth scheme. The equation for coupling the pressure to the velocity field is solved iteratively using the SIP method [17]. The subgrid-scale stresses appearing in the filtered Navier-Stokes equations are computed using (12). No-slip boundary conditions are employed on the surface of the stone and the tube walls. The physical properties used in the simulations are listed in Table 1. The temperature was assumed to be constant and it was set to T=302 K. Note that the vaporization pressure [18] at that temperature is 4500 Pa (absolute). The acceleration due to gravity is g=9.81 m/s2.The Lagrangian (bubble-based) simulation was performed to calculate the motion of bubbles in the liquid using (4). In brief, equivalent spherical bubbles were determined from irregular vapor structure by marking the cells with ϕw<0.95. Equation (4) was integrated using the generalized Verlet algorithm. The time step used in the simulation was Δt = 5×10-7 s. The linear interpolation routines were used to communicate the information from grid nodes to particle positions and vice versa. The simulations were performed by using GPU computing [10].Results. Obviously, a single parameter such as the cavitation number cannot describe many of the complexities of cavitation. Falvey [1] suggested that for flow past stones as shown in Figure 7(a), cavitation will not occur if the cavitation number is greater than about 1.8. Figure 7(a) shows that cavitation bubbles form within the flow at the discharge Q=12.5 L/s. This indicates that the cavitation number was below 1.8. Figure 7(b) represents a view from below of an instantaneous configuration of bubbles in the cavitation water flow. In this case, a sheet consists of a large number of small bubbles which are extensively developed downstream of the step formed by the stones in the pipe. Here, cavitation is formed by turbulence in the shear zone which is produced by the sudden change of flow direction at the stone face. Figure 7(b) also shows large nonspherical cavities with patterns similar to those shown in Figure 3(b).Figure 7
Computer experiments at a cavitation number ofσh=0.9 in a duct. (a) A top view of the instantaneous configuration of bubbles in a cavitating water flow of water through a pipe with inner diameter of D=0.024 m at the discharge of Q=12.5 L/s and at the cavitation number of σh=0.9. (b) A view from below of the bubbles in the cavitating flow of water through the pipe. Here, the large cavities consisting of small individual bubbles are illustrated. (c) Computed contours of the averaged vapor volume fraction on thexz-plane for the cavitating water flow in (a). Here, a color code as illustrated below the stones is used to indicate the values of the averaged vapor volume fraction on the surfaces of the stones. (d)–(g) Variations of ϕ-v as a function of lw at the sampling ports 1 through 16, respectively.
(a)
(b)
(c)
(d)
(e)
(f)
(g)Following Falvey [1], it may be estimated that the system in Figure 7(a) operates at a cavitation number of σh=0.9. For this case, the reference location is immediately upstream of the first stone and at its maximum height. The computed length of the sheet is quite comparable with that shown in Figure 3(b). It would be expected that for even lower values of the cavitation number, the clouds form one long supercavitating pocket.Figure7(c) depicts the computed instantaneous contours of the averaged vapor volume fraction in thexz-plane. In this figure, the locations at which cavitation occurred on the surfaces of the stones are presented by the use of color coding. As can be seen from Figure 7(c), the process starts with the occurrence of swarm of bubbles in the small regions of the first stone at sampling port 3 located at x=0.015 m. However, no evidence of cavitation can be found at sampling port 4 located at x=0.022 m in a narrow cleft between the stones where the local pressure is much higher than that of the vapor pressure. This figure also shows that a large number of small bubbles are extensively produced on the trailing edge of the second stone close to sampling port 8 located at x=0.045 m. Sampling ports 1 through 16 are all located in the xz-plane. The distances of these ports from the inlet of the pipe are x=5×10-5, 0.0077, 0.015, 0.022, 0.025, 0.029, 0.035, 0.045, 0.055, 0.065, 0.09, 0.12, 0.2, 0.25, 0.29, and 0.298 m, respectively.Figure7(d) illustrates the averaged vapor volume fraction, ϕ-v, as a function of the vertical distance from the pipe wall, lw, at sampling ports 1 through 4. Here, the squares, circles, diamonds, and left triangles represent the computed vapor volume fraction at ports 1, 2, 3, and 4, respectively. This figure indicates the occurrence of a swarm of bubbles at sampling port 3.The squares, circles, diamonds, and left triangles in Figure7(e) represent variations of ϕ-v as a function of lw at sampling ports 5, 6, 7, and 8, respectively. This figure indicates that cavitation reoccurred at sampling port 7 and became extensive at port 8. These complexities cannot be described with a single parameter such as cavitation number.Figure7(f) shows variations of ϕ-v as a function of lw at sampling ports 9, 10, 11, and 12 using the squares, circles, diamonds, and left triangles, respectively. As can be seen from Figure 7(f), the vapor bubbles roll up into a larger volume and become cloudy as they are transported further downstream and leave the entrance region of the pipe.Figure7(g) represents the final stage of the process which is an ongoing process of bubble formation and deformation. This figure depicts variations of ϕ-v as a function of lw at sampling ports 13, 14, 15, and 16 using the squares, circles, diamonds, and left triangles, respectively.In brief, the appearance of visible cavitation in the flowing water in the pipe as shown in Figure7 was preceded by the occurrence of small bubbles in the small area on the surface of the stones. This observation highlights the importance of bubbles as cavitation nuclei which has been known for a long time [6].Figure7(d) indicated that a swarm of small bubbles occurred at sampling port 3. However, no evidence of cavitation was found at sampling port 4 which was located further downstream from port 3. To address this issue, the pressure field is illustrated in thexz-plane in Figure 8(a). Here, the dimensionless pressure Π* is defined as Ln(p)/Ln(pv), where p represents the local pressure. Figure 8(b) depicts variations of the dimensionless pressure p* (defined as p/pv) as a function of lw at sampling ports 1, 2, 3, and 4 using the squares, circles, diamonds, and left triangles, respectively. This figure indicates that the local pressure in the narrow cleft between the stones is approximately 50 times higher than that of the vapor pressure.Figure 8
Cavitation in a tube. (a) The instantaneous computed contours of the dimensionless pressureΠ* on thexz-plane bubbles in a cavitation flow of water at the cavitation number of σh=0.9. The locations of the sampling ports are the same as those in Figure 7(c). (b)–(e) Variations of the dimensionless pressure p* as a function of lw. (f) The instantaneous computed velocity field of a cavitating water flow through the pipe with inner diameter of D=0.024 m at the discharge of Q=12.5 L/s.
(a)
(b)
(c)
(d)
(e)
(f)Figure8(c) illustrates variations of the dimensionless pressure, p*, as a function of lw at sampling ports 5, 6, 7, and 8 using squares, circles, diamonds, and left triangles, respectively. As can be seen from this figure, the local pressure at port 7 on the surface of the second stone reached the vapor pressure. Consequently, the small bubbles were formed on the trailing edge of the second stone as shown in Figure 7(e). The appearance of visible cavitation as shown in Figure 7(a) was preceded by the occurrence of these bubbles in the small area around sampling port 5 on the surface of the second stone. The computed low pressure at port 8 (which is the same as that measured at the PMP in Figure 3(a)) provides a favorable situation for a sheet, which consists of a large number of small bubbles, to be developed and transported downstream as depicted in Figures 7(a) and 7(b). The low pressure condition extends from sampling port 9 to sampling port 14 as shown in Figures 8(d) and 8(e). The squares, circles, diamonds, and left triangles in Figure 8(d) represent variations of p* as a function of lw at sampling ports 9, 10, 11, and 12, respectively.In the aforementioned region the vapor bubbles roll up into a larger volume and become less extensive as they are transported further downstream. The volume fraction of the vapor bubble is lower at the upper part of the pipe between ports 9 and 14 due to the fact that the pressure slightly increases withlw.The squares, circles, diamonds, and left triangles in Figure8(e) represent variations of p* as a function of lw at sampling ports 13, 14, 15, and 16, respectively. This figure indicates that the pressure increased from ports 14 to 16 and reached the atmospheric pressure at the pipe outlet at which the vapor bubbles vanished.The area-weighted average pressure at the inlet of the pipe was found to be 257 kPa. The results show agreement between simulation and experiment. The model is therefore deemed to be sufficiently flexible to capture a number of the interesting qualitative features of the cavitation process.Figure8(f) shows the computed velocity vector filed in thexz-plane. This figure indicates a complex shearing flow around sampling port 8 at which cavitation was formed. As stated earlier, the system in Figure 7(a) operates at a cavitation number of σh=0.9. Note that the computed pressure at sampling port 1 immediately upstream of the first stone and at its maximum height is approximately 300 kPa. In this case, the computed average water velocity is 25 m/s. The area-weighted average pressure at the inlet of the pipe was found to be 257 kPa.In the following section, the model is employed to check the reliability of the predicted cavitation zone as shown in Figure2.


### 3.3.4. Cavitation in a Spillway Chute

The previous section described the complexity of the occurrence of cavitation. As stated earlier, Figure2 indicates that the cavitation number should be quite low at approximately 135 m from the end of the chute of Kárahnjúkar spillway. Apparently, the designers of the Kárahnjúkar spillway decided to use an aerator at that location (as shown in Figure 1) in order to protect the spillway surface from cavitation.Recall that Semenkov and Saranchev [19] suggested that the beginning of the spillway is the cavitation-hazardous section. Particularly for large discharges, natural air entrainment is absent in more than 60 m from the inlet [19]. The bottom turbulent layer is shielded above by the undisturbed potential flow core, which reduces diffusion air entrainment [19]. Consequently, on the entire length of the cavitation-hazardous section of the spillway passage the volume air concentration in the boundary layer with thickness 0.2–0.25 m would be less than 7-8%. In this case, artificial air entrainment is needed to protect the surface of the spillway [20].In this section, the usefulness of using artificial air at 125 m from the end of the chute of the Kárahnjúkar spillway is investigated using the mathematical model as detailed in the preceding section. To this end, a simplified model of Kárahnjúkar spillway with a hypothetical hole with diameter of 0.02 m is used in the simulation. Figure9(a) shows the hypothetical hole which is essential to initiate cavitation at the cavitation number of 0.25. The simplified model of Kárahnjúkar spillway is 414 m long with a slope of approximately α = 11°.Figure 9
Cavitation in a spillway chute. (a) Schematic of a hypothetical hole with a diameterdh=0.02 on the surface of a model of Kárahnjúkar spillway at a distance from the inlet of Lh=300 m. The total length of the spillway is LT=414 m and its slope is α = 11°. Inset: the computed contours of velocity in the spillway. The distances of the sampling ports 1, 2, 3, and 4 from the inlet are 0, 107, 300, and 414 m, respectively. (b) The averaged air volume fraction as a function of lw at the sampling ports 1, 2, 3, and 4. Here, the squares, circles, diamonds, and left triangles represent the air volume fraction at the sampling ports 1 through 4, respectively. (c) The computed velocity vector field around the hole at the sampling port 3. Here, a color code as illustrated on the left-hand side of the hole is used to indicate the values of the local pressure.
(a)
(b)
(c)The inlet boundary condition atx=0 was set to V=20 m/s. The height of the water at the inlet was 3 m. Furthermore, it was assumed that the top surface of the model is open to the atmosphere. The outlet boundary condition was set to p=patm. Periodic boundary conditions were used in the z-direction. The no-slip boundary condition was applied on the walls of the model. To capture the effects of viscous boundary layers, a wall-function was used to specify velocity at the forcing points as [21](13)uu*=y+y+≤5-1.15576+1.47869y+-0.0527848y+2+0.000655444y+35<y+≤30ln⁡y+Κ+5.1y+>30,where y+ represents dimensionless, sublayer-scaled distance, and Κ is the von Kármán constant (typically the value 0.41 is used).The inset of Figure9(a) represents the computed velocity contours in the model of the Kárahnjúkar spillway. The computed bubble volume fraction as a function of vertical distance from the wall lw at four sampling ports is illustrated in Figure 9(b). In this case, the bubble is an air bubble which may contain vapor. As can be seen from Figure 9(a), the hypothetical hole is located at 135 m from the outlet of the model. Figure 9(b) indicates that the thickness of the air-water mixture increases with x. Here, the squares, circles, diamonds, and left triangles represent the bubble volume fraction at sampling ports 1, 2, 3, and 4, respectively. As can be seen from Figure 9(b), the bubble volume fraction at sampling port 2 is lower than 7%. This finding highlights the need of artificial aerators between sampling ports 1 and 2 in order to protect the surface of the spillway.Figure9(b) indicates that the bubble volume fraction at the hypothetical hole (i.e., at sampling port 3) is approximately 38% which is in accordance with Semenkov and Saranchev [19] and very much above the “cavitation safe” value of 7% of entrained air. In this case, the vapor volume fraction is quite low.Figure9(c) illustrates the computed velocity vector field. A color code as illustrated on the left hand side of the hole in Figure 9(c) is used to indicate the values of the local pressure. This figure indicates that the water flows mainly over the hole. A small part of the water stream hits the hole’s face and reticulates inside the hole. As can be seen from Figure 9(c), the pressure reaches a maximum of 2500 kPa and diverts the flow to an upward direction over a little edge and in a downward direction which initiates the water swirl in the hole. The local pressure in the center of the hole is slightly lower than 5 kPa. This finding indicates that cavitation might occur at 125 m from the end of the chute of the Kárahnjúkar spillway. Indeed, cavitation will occur for the flow past the sudden into-the-flow offset introduced by the hypothetical hole as noted in Figure 9(c). This figure indicates a low-pressure region and thus potential cavitation might also occur shortly after passing the hole’s edge. However, the air volume fraction at the hole as shown in Figure 9(b) is approximately 38%. This finding suggests that the use of an aerator at port 3 in order to protect the surface of the spillway might not be needed. Adding extra air at sampling port 3 would produce a very thick air-water mist flow at the outlet of the spillway. This will be further analyzed in the following section.Figure9(b) also indicates that in the absence of an artificial aerator the computed thickness of the air-water mixture at the outlet of the spillway is approximately 8 m. In addition, the average bubble volume fraction is higher than 50%.


### 3.3.5. Artificial Aeration

The designers of the Kárahnjúkar spillway were challenged in their task due to the lack of literature on air entrainment for spillways. As mentioned in the preceding section, an aerator must be introduced to prevent cavitation in the Kárahnjúkar spillway at high discharges. The air volume fraction in the spillway downstream from sampling port 2 in Figure9(b) is well above the 7% minimum requirement of its protection against cavitation. Hence, artificial aeration would only be needed between sampling ports 1 and 2 illustrated in Figure 9(a).In this section, the mathematical model described in the preceding section will be used to analyze the Kárahnjúkar spillway with a bottom aerator. Figure10 illustrates a CFD model which can be imagined as a 1 : 15 Froude scale of the Kárahnjúkar bottom aerator. Figure 1 shows an aerial view of the Kárahnjúkar spillway and its bottom aerator located in eastern Iceland.Figure 10
A CFD model for the Kárahnjúkar bottom aerator. (a) Perspective view of the simulation model and selected nomenclature. (b) The front view of the spillway with bottom aerator. Here,δ = 0.02 m, L = 2.42 m, Ha = 0.44 m, Hw = 0.11 m, L1 = 0.11 m, L2 = 2.15 m, h0 = 0.022, Hout = 0.55 m, θ1 = 20°, and θ2=30°.
(a)
(b)By implementing periodic boundary conditions in they-direction, a small part of the system that is far from the vertical side walls was simulated. The model flow rate was set to vw0 = 20 m/s at the inlet with a size of Hw = 0.11 m, corresponding to Q= 2250 m3/s.As can be seen from Figure11, the water stream is deflected by a deflector in order to ventilate the nappe to the atmosphere. The connection to the atmosphere is made via a channel under the ski jump of the aerator. Note that the air entraining capacity of the aerator depends on the take-off distance from the free flow. The air pressure is atmospheric at the inlets and the outlet.Figure 11
Computed contours of the velocity field for the bottom aerator of Kárahnjúkar spillway. Inset: a pair of vortices in the ventilation zone.Figure12(a) and its insets illustrate the computed contours and vector velocity field of a more optimized design for the aerator of the Kárahnjúkar spillway. Here, vortices can be seen to be created at the ventilation gate and to detach periodically from its other side. A pair of vortices is magnified and replotted in the inset. The air enters the ventilation zone with a normalized velocity of v*=3.81. Here, the velocity is normalized with the water velocity at inlet vw0. Figure 12(b) represents an instantaneous configuration of bubbles in the spillway with an optimized bottom aerator. The high stream of air tends to move the water stream up. In this case, the air concentration in the ventilated zone appears to be more than that needed to prevent cavitation. Figure 1 shows a highly aerated water stream which is likely to serve as a source for an intense canyon erosion.Figure 12
The computed flow field for the aerator of Kárahnjúkar spillway. (a) Computed contours of velocity field of a more optimized design for the aerator of Kárahnjúkar spillway. Inset: vortical structures in the ventilation zone. (b) Instantaneous configuration of bubbles in a spillway with ski jump aerator. Here,Lnb = 0.62 m. (c) Air concentration distribution versus z in the aeration region. Here, right triangles, squares, left triangles, and diamonds represent the air concentration at the sampling ports S1, S2, S3, and S4, respectively.
(a)
(b)
(c)When water flows over a spillway, vortex sheet instability at the air-water interface known as Kelvin-Helmholtz instability can occur. The instability will be in the form of waves being generated on the water surface which can initiate natural air entrainment.The air entrainment process is strengthened when a growing boundary layer in the water flow reaches the free surface. In this case, the entrained air from the free surface can reach the spillway surface to minimize the danger of cavitation. As can be seen from Figure12, the details of air entrainment through the upper air-water interface of the jet can be predicted by the model. In addition, the details of air entrainment through the lower interface are very similar to those presented in [10]. Here, Lnb represents the length at which the air entrainment through the lower interface reaches that through the upper interface. Figure 12(c) depicts the variations of air concentration as a function of z at four different sampling ports as shown in (b).The region at which the natural air entrainment effectively eliminates the risk of damage due to cavitation is called the fully aerated region. The artificial aeration process discussed above would not be needed in the fully aerated region.Figure13(a) illustrates the computed contours of the velocity field for a hypothetical case in which the discharge of the Kárahnjúkar spillway is Q= 3150 m3/s. Figure 13(b) indicates that the air enters the ventilation zone with a normalized velocity of v*=5.6. The vortex shedding of the aeration gate is quite intense. Therefore, the air-water mixture entering the canyon could cause extensive erosion. Recall that the canyon shown in Figure 1 suffered a disastrous erosion in 2009.Figure 13
The computed flow field for a hypothetical case with a discharge ofQ=3150 m3/s. (a) Computed contours of velocity field. (b) Computed vector plot of vertical structures in the ventilation zone. (c) and (d) Variations of the normalized velocity v* as a function of z. Squares, circles, diamonds, left triangles, gradients, and deltas represent the normalized velocity at the sampling ports S1, S6, S5, S2, S3, and S4, respectively.
(a)
(b)
(c)
(d)


## 3.3.1. Air Bubbles

Figures4(b) and 4(c) show a computed large bubble in an aerated tank [7]. The bubble is a cluster of small bubbles with a size of the critical air bubble radius. Figure 4(d) depicts the denoised version of the bubble using multiscale image denoising based on a multiscale representation of the images [14]. The denoised version of the bubble-like structure in Figure 4(b) is very similar to the small bubbles observed in an aerated tank [7].Figure 4
The interaction of small bubbles in an aerated tank. (a) The Lennard-Jones 12-6 pair potential. (b) Computed bubble-like structure in an aerated tank at the Reynolds number of 29500. (c) A different view of (b). (d) The denoised version of (c). The length of the bubble is approximately 1 cm and its height is 0.65 cm. (e) Configuration of a three bubble-like structure in an aerated tank. (f) The bubbles in (e) coalesce to form a larger bubble after 0.3 sec. (g) The denoised version of (f).Figure 5
A CFD model for a stone in Figure3. (a) Taking 2D images of the stone by moving camera’s position at regular intervals. (b) Automatically aligning the images with the others. (c) Reconstruction of a 3D model from 2D image.
(a)
(b)
(c)Figures4(e) and 4(f) are snapshots of bubbles separated by t=0.3 sec. As can be seen from Figure 4(f), three bubbles in Figure 4(e) coalesce to form one large bubble. The denoised version depicted in Figure 4(g) is very similar to the bubbles with complex geometries shown in [7].To further assess the quality of the model described in the preceding section, the cavitation water flow in Figure3 is simulated in the following section.


## 3.3.2. 3D Model Using 2D Images

Apparently, it is not possible to reconstruct a 3D model from a single image of the stone as shown in the inset of Figure5(a). To get around this limitation, an optical processing algorithm can be developed that employs multiple photographs taken at every 15 degrees, from the above and from the side to create a 3D model of the stone. Figure 5(a) illustrates the positions of the camera for taking about 40 photographs of the stone. The images acquired are the input of the Autodesk 123D Catch 3D scanning software [15].Figure5(b) shows the basics of the stitching process for the accurate creation of a 3D model of the stone. Figure 5(c) represents a medium quality 3D mesh of the stone which is used to simulate sheet cavitation in the following section.


## 3.3.3. Sheet Cavitation

Aspects of the Simulation. As mentioned earlier, turbulence modeling has a critical role in cavitation prediction. Cavitation flow in the liquid-vapor region is locally compressible [9, 16]. To capture the shedding dynamics and the unsteadiness of cavitation the modified form of LES as detailed in the preceding section would be required. Note that (5) and (6) look like the equations of motion of a fluid with variable density ρ=ρw1-ϕv.In this section, the mathematical model described in the preceding section will be used to analyze sheet cavitation in a pipe flow as shown in Figure3. The discharge was set to Q=12.5 L/s. The tube inner diameter was D=0.024 m and its length was L=0.3 m.The length and the height of the stones in the simulations are exactly the same as those used in the experiments. The model of the stone is shown in Figure5(c). It is likely that cavitation as shown in Figure 3(b) starts at minute cracks on the otherwise smooth surface of the stones. The aforementioned system may be characterized as belonging to the singular roughness category [1].In brief, the Lagrangian (bubble-based) and the Eulerian (grid-based) methods are used to simulate the cavitation singular roughness in a tube flow. Figure6(a) illustrates the model used in the simulations. Here, the dimensions are Lp = 0.3 m, D=0.024 cm, and hs = 0.004 m. Figure 6(b) depicts a top view of the entrance region of the pipe in which the positions of the stones are clearly shown. The total length of the pipe covered by stones is l2 = 0.037 m, and the pipe’s length covered by the first stone is l1 = 0.0185 m. The length and height of the stones are exactly the same as those shown in Figure 3(a). Figure 6(c) shows the grid which consists of more than 5×106 tetrahedral meshes.Figure 6
A CFD model for the tube in Figure3. (a) Schematic of the pipe used in sheet cavitation simulations with some nomenclatures. (b) A top view of the entrance region of the pipe. (c) The grid used in the cavitation simulations. Insets: the grid is magnified and replotted.
(a)
(b)
(c)The filtered Navier-Stokes equations (5) and (6) are discretized with the finite-volume method on a staggered Cartesian grid. Convective and diffusive fluxes are approximated with central differences of second-order accuracy and time advancement is achieved by a second-order, explicit Adams-Bashforth scheme. The equation for coupling the pressure to the velocity field is solved iteratively using the SIP method [17]. The subgrid-scale stresses appearing in the filtered Navier-Stokes equations are computed using (12). No-slip boundary conditions are employed on the surface of the stone and the tube walls. The physical properties used in the simulations are listed in Table 1. The temperature was assumed to be constant and it was set to T=302 K. Note that the vaporization pressure [18] at that temperature is 4500 Pa (absolute). The acceleration due to gravity is g=9.81 m/s2.The Lagrangian (bubble-based) simulation was performed to calculate the motion of bubbles in the liquid using (4). In brief, equivalent spherical bubbles were determined from irregular vapor structure by marking the cells with ϕw<0.95. Equation (4) was integrated using the generalized Verlet algorithm. The time step used in the simulation was Δt = 5×10-7 s. The linear interpolation routines were used to communicate the information from grid nodes to particle positions and vice versa. The simulations were performed by using GPU computing [10].Results. Obviously, a single parameter such as the cavitation number cannot describe many of the complexities of cavitation. Falvey [1] suggested that for flow past stones as shown in Figure 7(a), cavitation will not occur if the cavitation number is greater than about 1.8. Figure 7(a) shows that cavitation bubbles form within the flow at the discharge Q=12.5 L/s. This indicates that the cavitation number was below 1.8. Figure 7(b) represents a view from below of an instantaneous configuration of bubbles in the cavitation water flow. In this case, a sheet consists of a large number of small bubbles which are extensively developed downstream of the step formed by the stones in the pipe. Here, cavitation is formed by turbulence in the shear zone which is produced by the sudden change of flow direction at the stone face. Figure 7(b) also shows large nonspherical cavities with patterns similar to those shown in Figure 3(b).Figure 7
Computer experiments at a cavitation number ofσh=0.9 in a duct. (a) A top view of the instantaneous configuration of bubbles in a cavitating water flow of water through a pipe with inner diameter of D=0.024 m at the discharge of Q=12.5 L/s and at the cavitation number of σh=0.9. (b) A view from below of the bubbles in the cavitating flow of water through the pipe. Here, the large cavities consisting of small individual bubbles are illustrated. (c) Computed contours of the averaged vapor volume fraction on thexz-plane for the cavitating water flow in (a). Here, a color code as illustrated below the stones is used to indicate the values of the averaged vapor volume fraction on the surfaces of the stones. (d)–(g) Variations of ϕ-v as a function of lw at the sampling ports 1 through 16, respectively.
(a)
(b)
(c)
(d)
(e)
(f)
(g)Following Falvey [1], it may be estimated that the system in Figure 7(a) operates at a cavitation number of σh=0.9. For this case, the reference location is immediately upstream of the first stone and at its maximum height. The computed length of the sheet is quite comparable with that shown in Figure 3(b). It would be expected that for even lower values of the cavitation number, the clouds form one long supercavitating pocket.Figure7(c) depicts the computed instantaneous contours of the averaged vapor volume fraction in thexz-plane. In this figure, the locations at which cavitation occurred on the surfaces of the stones are presented by the use of color coding. As can be seen from Figure 7(c), the process starts with the occurrence of swarm of bubbles in the small regions of the first stone at sampling port 3 located at x=0.015 m. However, no evidence of cavitation can be found at sampling port 4 located at x=0.022 m in a narrow cleft between the stones where the local pressure is much higher than that of the vapor pressure. This figure also shows that a large number of small bubbles are extensively produced on the trailing edge of the second stone close to sampling port 8 located at x=0.045 m. Sampling ports 1 through 16 are all located in the xz-plane. The distances of these ports from the inlet of the pipe are x=5×10-5, 0.0077, 0.015, 0.022, 0.025, 0.029, 0.035, 0.045, 0.055, 0.065, 0.09, 0.12, 0.2, 0.25, 0.29, and 0.298 m, respectively.Figure7(d) illustrates the averaged vapor volume fraction, ϕ-v, as a function of the vertical distance from the pipe wall, lw, at sampling ports 1 through 4. Here, the squares, circles, diamonds, and left triangles represent the computed vapor volume fraction at ports 1, 2, 3, and 4, respectively. This figure indicates the occurrence of a swarm of bubbles at sampling port 3.The squares, circles, diamonds, and left triangles in Figure7(e) represent variations of ϕ-v as a function of lw at sampling ports 5, 6, 7, and 8, respectively. This figure indicates that cavitation reoccurred at sampling port 7 and became extensive at port 8. These complexities cannot be described with a single parameter such as cavitation number.Figure7(f) shows variations of ϕ-v as a function of lw at sampling ports 9, 10, 11, and 12 using the squares, circles, diamonds, and left triangles, respectively. As can be seen from Figure 7(f), the vapor bubbles roll up into a larger volume and become cloudy as they are transported further downstream and leave the entrance region of the pipe.Figure7(g) represents the final stage of the process which is an ongoing process of bubble formation and deformation. This figure depicts variations of ϕ-v as a function of lw at sampling ports 13, 14, 15, and 16 using the squares, circles, diamonds, and left triangles, respectively.In brief, the appearance of visible cavitation in the flowing water in the pipe as shown in Figure7 was preceded by the occurrence of small bubbles in the small area on the surface of the stones. This observation highlights the importance of bubbles as cavitation nuclei which has been known for a long time [6].Figure7(d) indicated that a swarm of small bubbles occurred at sampling port 3. However, no evidence of cavitation was found at sampling port 4 which was located further downstream from port 3. To address this issue, the pressure field is illustrated in thexz-plane in Figure 8(a). Here, the dimensionless pressure Π* is defined as Ln(p)/Ln(pv), where p represents the local pressure. Figure 8(b) depicts variations of the dimensionless pressure p* (defined as p/pv) as a function of lw at sampling ports 1, 2, 3, and 4 using the squares, circles, diamonds, and left triangles, respectively. This figure indicates that the local pressure in the narrow cleft between the stones is approximately 50 times higher than that of the vapor pressure.Figure 8
Cavitation in a tube. (a) The instantaneous computed contours of the dimensionless pressureΠ* on thexz-plane bubbles in a cavitation flow of water at the cavitation number of σh=0.9. The locations of the sampling ports are the same as those in Figure 7(c). (b)–(e) Variations of the dimensionless pressure p* as a function of lw. (f) The instantaneous computed velocity field of a cavitating water flow through the pipe with inner diameter of D=0.024 m at the discharge of Q=12.5 L/s.
(a)
(b)
(c)
(d)
(e)
(f)Figure8(c) illustrates variations of the dimensionless pressure, p*, as a function of lw at sampling ports 5, 6, 7, and 8 using squares, circles, diamonds, and left triangles, respectively. As can be seen from this figure, the local pressure at port 7 on the surface of the second stone reached the vapor pressure. Consequently, the small bubbles were formed on the trailing edge of the second stone as shown in Figure 7(e). The appearance of visible cavitation as shown in Figure 7(a) was preceded by the occurrence of these bubbles in the small area around sampling port 5 on the surface of the second stone. The computed low pressure at port 8 (which is the same as that measured at the PMP in Figure 3(a)) provides a favorable situation for a sheet, which consists of a large number of small bubbles, to be developed and transported downstream as depicted in Figures 7(a) and 7(b). The low pressure condition extends from sampling port 9 to sampling port 14 as shown in Figures 8(d) and 8(e). The squares, circles, diamonds, and left triangles in Figure 8(d) represent variations of p* as a function of lw at sampling ports 9, 10, 11, and 12, respectively.In the aforementioned region the vapor bubbles roll up into a larger volume and become less extensive as they are transported further downstream. The volume fraction of the vapor bubble is lower at the upper part of the pipe between ports 9 and 14 due to the fact that the pressure slightly increases withlw.The squares, circles, diamonds, and left triangles in Figure8(e) represent variations of p* as a function of lw at sampling ports 13, 14, 15, and 16, respectively. This figure indicates that the pressure increased from ports 14 to 16 and reached the atmospheric pressure at the pipe outlet at which the vapor bubbles vanished.The area-weighted average pressure at the inlet of the pipe was found to be 257 kPa. The results show agreement between simulation and experiment. The model is therefore deemed to be sufficiently flexible to capture a number of the interesting qualitative features of the cavitation process.Figure8(f) shows the computed velocity vector filed in thexz-plane. This figure indicates a complex shearing flow around sampling port 8 at which cavitation was formed. As stated earlier, the system in Figure 7(a) operates at a cavitation number of σh=0.9. Note that the computed pressure at sampling port 1 immediately upstream of the first stone and at its maximum height is approximately 300 kPa. In this case, the computed average water velocity is 25 m/s. The area-weighted average pressure at the inlet of the pipe was found to be 257 kPa.In the following section, the model is employed to check the reliability of the predicted cavitation zone as shown in Figure2.


## 3.3.4. Cavitation in a Spillway Chute

The previous section described the complexity of the occurrence of cavitation. As stated earlier, Figure2 indicates that the cavitation number should be quite low at approximately 135 m from the end of the chute of Kárahnjúkar spillway. Apparently, the designers of the Kárahnjúkar spillway decided to use an aerator at that location (as shown in Figure 1) in order to protect the spillway surface from cavitation.Recall that Semenkov and Saranchev [19] suggested that the beginning of the spillway is the cavitation-hazardous section. Particularly for large discharges, natural air entrainment is absent in more than 60 m from the inlet [19]. The bottom turbulent layer is shielded above by the undisturbed potential flow core, which reduces diffusion air entrainment [19]. Consequently, on the entire length of the cavitation-hazardous section of the spillway passage the volume air concentration in the boundary layer with thickness 0.2–0.25 m would be less than 7-8%. In this case, artificial air entrainment is needed to protect the surface of the spillway [20].In this section, the usefulness of using artificial air at 125 m from the end of the chute of the Kárahnjúkar spillway is investigated using the mathematical model as detailed in the preceding section. To this end, a simplified model of Kárahnjúkar spillway with a hypothetical hole with diameter of 0.02 m is used in the simulation. Figure9(a) shows the hypothetical hole which is essential to initiate cavitation at the cavitation number of 0.25. The simplified model of Kárahnjúkar spillway is 414 m long with a slope of approximately α = 11°.Figure 9
Cavitation in a spillway chute. (a) Schematic of a hypothetical hole with a diameterdh=0.02 on the surface of a model of Kárahnjúkar spillway at a distance from the inlet of Lh=300 m. The total length of the spillway is LT=414 m and its slope is α = 11°. Inset: the computed contours of velocity in the spillway. The distances of the sampling ports 1, 2, 3, and 4 from the inlet are 0, 107, 300, and 414 m, respectively. (b) The averaged air volume fraction as a function of lw at the sampling ports 1, 2, 3, and 4. Here, the squares, circles, diamonds, and left triangles represent the air volume fraction at the sampling ports 1 through 4, respectively. (c) The computed velocity vector field around the hole at the sampling port 3. Here, a color code as illustrated on the left-hand side of the hole is used to indicate the values of the local pressure.
(a)
(b)
(c)The inlet boundary condition atx=0 was set to V=20 m/s. The height of the water at the inlet was 3 m. Furthermore, it was assumed that the top surface of the model is open to the atmosphere. The outlet boundary condition was set to p=patm. Periodic boundary conditions were used in the z-direction. The no-slip boundary condition was applied on the walls of the model. To capture the effects of viscous boundary layers, a wall-function was used to specify velocity at the forcing points as [21](13)uu*=y+y+≤5-1.15576+1.47869y+-0.0527848y+2+0.000655444y+35<y+≤30ln⁡y+Κ+5.1y+>30,where y+ represents dimensionless, sublayer-scaled distance, and Κ is the von Kármán constant (typically the value 0.41 is used).The inset of Figure9(a) represents the computed velocity contours in the model of the Kárahnjúkar spillway. The computed bubble volume fraction as a function of vertical distance from the wall lw at four sampling ports is illustrated in Figure 9(b). In this case, the bubble is an air bubble which may contain vapor. As can be seen from Figure 9(a), the hypothetical hole is located at 135 m from the outlet of the model. Figure 9(b) indicates that the thickness of the air-water mixture increases with x. Here, the squares, circles, diamonds, and left triangles represent the bubble volume fraction at sampling ports 1, 2, 3, and 4, respectively. As can be seen from Figure 9(b), the bubble volume fraction at sampling port 2 is lower than 7%. This finding highlights the need of artificial aerators between sampling ports 1 and 2 in order to protect the surface of the spillway.Figure9(b) indicates that the bubble volume fraction at the hypothetical hole (i.e., at sampling port 3) is approximately 38% which is in accordance with Semenkov and Saranchev [19] and very much above the “cavitation safe” value of 7% of entrained air. In this case, the vapor volume fraction is quite low.Figure9(c) illustrates the computed velocity vector field. A color code as illustrated on the left hand side of the hole in Figure 9(c) is used to indicate the values of the local pressure. This figure indicates that the water flows mainly over the hole. A small part of the water stream hits the hole’s face and reticulates inside the hole. As can be seen from Figure 9(c), the pressure reaches a maximum of 2500 kPa and diverts the flow to an upward direction over a little edge and in a downward direction which initiates the water swirl in the hole. The local pressure in the center of the hole is slightly lower than 5 kPa. This finding indicates that cavitation might occur at 125 m from the end of the chute of the Kárahnjúkar spillway. Indeed, cavitation will occur for the flow past the sudden into-the-flow offset introduced by the hypothetical hole as noted in Figure 9(c). This figure indicates a low-pressure region and thus potential cavitation might also occur shortly after passing the hole’s edge. However, the air volume fraction at the hole as shown in Figure 9(b) is approximately 38%. This finding suggests that the use of an aerator at port 3 in order to protect the surface of the spillway might not be needed. Adding extra air at sampling port 3 would produce a very thick air-water mist flow at the outlet of the spillway. This will be further analyzed in the following section.Figure9(b) also indicates that in the absence of an artificial aerator the computed thickness of the air-water mixture at the outlet of the spillway is approximately 8 m. In addition, the average bubble volume fraction is higher than 50%.


## 3.3.5. Artificial Aeration

The designers of the Kárahnjúkar spillway were challenged in their task due to the lack of literature on air entrainment for spillways. As mentioned in the preceding section, an aerator must be introduced to prevent cavitation in the Kárahnjúkar spillway at high discharges. The air volume fraction in the spillway downstream from sampling port 2 in Figure9(b) is well above the 7% minimum requirement of its protection against cavitation. Hence, artificial aeration would only be needed between sampling ports 1 and 2 illustrated in Figure 9(a).In this section, the mathematical model described in the preceding section will be used to analyze the Kárahnjúkar spillway with a bottom aerator. Figure10 illustrates a CFD model which can be imagined as a 1 : 15 Froude scale of the Kárahnjúkar bottom aerator. Figure 1 shows an aerial view of the Kárahnjúkar spillway and its bottom aerator located in eastern Iceland.Figure 10
A CFD model for the Kárahnjúkar bottom aerator. (a) Perspective view of the simulation model and selected nomenclature. (b) The front view of the spillway with bottom aerator. Here,δ = 0.02 m, L = 2.42 m, Ha = 0.44 m, Hw = 0.11 m, L1 = 0.11 m, L2 = 2.15 m, h0 = 0.022, Hout = 0.55 m, θ1 = 20°, and θ2=30°.
(a)
(b)By implementing periodic boundary conditions in they-direction, a small part of the system that is far from the vertical side walls was simulated. The model flow rate was set to vw0 = 20 m/s at the inlet with a size of Hw = 0.11 m, corresponding to Q= 2250 m3/s.As can be seen from Figure11, the water stream is deflected by a deflector in order to ventilate the nappe to the atmosphere. The connection to the atmosphere is made via a channel under the ski jump of the aerator. Note that the air entraining capacity of the aerator depends on the take-off distance from the free flow. The air pressure is atmospheric at the inlets and the outlet.Figure 11
Computed contours of the velocity field for the bottom aerator of Kárahnjúkar spillway. Inset: a pair of vortices in the ventilation zone.Figure12(a) and its insets illustrate the computed contours and vector velocity field of a more optimized design for the aerator of the Kárahnjúkar spillway. Here, vortices can be seen to be created at the ventilation gate and to detach periodically from its other side. A pair of vortices is magnified and replotted in the inset. The air enters the ventilation zone with a normalized velocity of v*=3.81. Here, the velocity is normalized with the water velocity at inlet vw0. Figure 12(b) represents an instantaneous configuration of bubbles in the spillway with an optimized bottom aerator. The high stream of air tends to move the water stream up. In this case, the air concentration in the ventilated zone appears to be more than that needed to prevent cavitation. Figure 1 shows a highly aerated water stream which is likely to serve as a source for an intense canyon erosion.Figure 12
The computed flow field for the aerator of Kárahnjúkar spillway. (a) Computed contours of velocity field of a more optimized design for the aerator of Kárahnjúkar spillway. Inset: vortical structures in the ventilation zone. (b) Instantaneous configuration of bubbles in a spillway with ski jump aerator. Here,Lnb = 0.62 m. (c) Air concentration distribution versus z in the aeration region. Here, right triangles, squares, left triangles, and diamonds represent the air concentration at the sampling ports S1, S2, S3, and S4, respectively.
(a)
(b)
(c)When water flows over a spillway, vortex sheet instability at the air-water interface known as Kelvin-Helmholtz instability can occur. The instability will be in the form of waves being generated on the water surface which can initiate natural air entrainment.The air entrainment process is strengthened when a growing boundary layer in the water flow reaches the free surface. In this case, the entrained air from the free surface can reach the spillway surface to minimize the danger of cavitation. As can be seen from Figure12, the details of air entrainment through the upper air-water interface of the jet can be predicted by the model. In addition, the details of air entrainment through the lower interface are very similar to those presented in [10]. Here, Lnb represents the length at which the air entrainment through the lower interface reaches that through the upper interface. Figure 12(c) depicts the variations of air concentration as a function of z at four different sampling ports as shown in (b).The region at which the natural air entrainment effectively eliminates the risk of damage due to cavitation is called the fully aerated region. The artificial aeration process discussed above would not be needed in the fully aerated region.Figure13(a) illustrates the computed contours of the velocity field for a hypothetical case in which the discharge of the Kárahnjúkar spillway is Q= 3150 m3/s. Figure 13(b) indicates that the air enters the ventilation zone with a normalized velocity of v*=5.6. The vortex shedding of the aeration gate is quite intense. Therefore, the air-water mixture entering the canyon could cause extensive erosion. Recall that the canyon shown in Figure 1 suffered a disastrous erosion in 2009.Figure 13
The computed flow field for a hypothetical case with a discharge ofQ=3150 m3/s. (a) Computed contours of velocity field. (b) Computed vector plot of vertical structures in the ventilation zone. (c) and (d) Variations of the normalized velocity v* as a function of z. Squares, circles, diamonds, left triangles, gradients, and deltas represent the normalized velocity at the sampling ports S1, S6, S5, S2, S3, and S4, respectively.
(a)
(b)
(c)
(d)


## 4. Conclusions

A previous study revealed that the cavitation number of the Kárahnjúkar spillway is below 0.25 at a discharge of 2250 m3/s. Hence, the spillway was protected by additional aeration grooves located some 300 m from its inlet. In this work, an analysis of the problem based on physical and computational fluid dynamics (CFD) modeling is presented. The numerical modeling is based on the large eddy simulation technique (LES) combined with the discrete element method.A three-dimensional simulation of a cavitation tube flow at cavitation number of ca. 0.9 was performed to highlight the importance of artificial aeration for protecting a spillway. More simulations were performed using a CFD model which is a 1 : 15 Froude scale of the Kárahnjúkar bottom aerator in order to gain insights into how to diminish the hazards associated with cavitation and minimize canyon erosion. It is recommended to use the artificial aerators only at the distances less than 100 m from the inlet in order to minimize canyon erosion.The main finding was that the design of the aerator of the Kárahnjúkar spillway should be improved. It is highly unlikely that the artificial aerators would be needed at the distances of 300 m from the inlet of the spillway.



---

*Source: 872573-2015-06-23.xml*