Mathematical Problems in Engineering
Volume 2015 (2015), Article ID 872573, 16 pages
http://dx.doi.org/10.1155/2015/872573

Research Article

Simulation of Cavitation Water Flows

Piroz Zamankhan1,2

¹Institute of Geothermal Energy, Dedan Kimathi University of Technology, 101000 Nyeri, Kenya
²Chemical Engineering Department, School of Engineering, Collage of Agriculture, Engineering and Science, University of KwaZulu-Natal, 4041 Durban, South Africa

Received 20 July 2014; Revised 4 September 2014; Accepted 13 September 2014

Academic Editor: Junuthula N. Reddy

Copyright © 2015 Piroz Zamankhan. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

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.

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 Figure 1.

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, . 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 Figure 2 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 m³/s; red squares: 400 m³/s; green triangle: 800 m³/s; blue squares: 950 m³/s; yellow square with green star: 1350 m³/s; and orange dot: 2250 m³/s.

As a precaution the cavitation zone has been set to 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 Figure 2 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. Figure 3 shows a simple system in which sheet cavitation may be observed.

(a)

(b)

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 discharge Q = 12.5 L/s.

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 Figure 3. 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 Figure 3(a). The tests were performed in a horizontal galvanized steel pipe with an outer diameter of  m and the length of = 0.3 m. Two very similar aquarium stones with a length of 0.0215 m and a height of  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 Table 1. The discharge  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/m³)
Vaporization pressure (mmHg) [18]

Water (liquid)
0.001
996.4

Water (vapor)

0.5542

Air

1.175

The average pressure at the pipe inlet was approximately 2.5 bars (absolute) and the minimum measured pressure was 4.55 kPa (absolute). Figure 3(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 = 0.048 m downstream from the inlet of the test tube. Here, 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 of Q = 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  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.

Figure 3(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 . 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 Figure 3. 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 Figure 3(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 must satisfy the following condition [6]:

The bubble radius at the limit point, where , can be expressed in terms of material parameters as [5]

If , then the bubble is stable against infinitesimal changes in its radius. The bubble becomes unstable if . 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 located at and is given as

A simplified momentum equation for the th spherical bubble may be given as

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

The viscous stress tensor in (6) is defined as

In (6), may be given as

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 , , may be given as

Here, the changes in the bubble radius may be calculated using the famous Rayleigh-Plesset equation [6]

The vapor transport equation may be given as

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]

Here, the value of the backscatter parameter, , is set to . 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

Figures 4(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).

(a)

(b)

(c)

Figure 5: A CFD model for a stone in Figure 3. (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.

Figures 4(e) and 4(f) are snapshots of bubbles separated by  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 Figure 3 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 Figure 5(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].

Figure 5(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 .

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 Figure 3. The discharge was set to  L/s. The tube inner diameter was  m and its length was  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 Figure 5(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. Figure 6(a) illustrates the model used in the simulations. Here, the dimensions are = 0.3 m,  cm, and = 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 = 0.037 m, and the pipe’s length covered by the first stone is = 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 tetrahedral meshes.

(a)

(b)

(c)

Figure 6: A CFD model for the tube in Figure 3. (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.

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  K. Note that the vaporization pressure [18] at that temperature is 4500 Pa (absolute). The acceleration due to gravity is  m/s².

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 . Equation (4) was integrated using the generalized Verlet algorithm. The time step used in the simulation was =  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  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).

(a)

(b)

(c)

(d)

(e)

(f)

(g)

Figure 7: Computer experiments at a cavitation number of 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  m at the discharge of  L/s and at the cavitation number of . (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 the xz-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 as a function of at the sampling ports 1 through 16, respectively.

Following Falvey [1], it may be estimated that the system in Figure 7(a) operates at a cavitation number of . 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.

Figure 7(c) depicts the computed instantaneous contours of the averaged vapor volume fraction in the xz-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  m. However, no evidence of cavitation can be found at sampling port 4 located at  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  m. Sampling ports 1 through 16 are all located in the -plane. The distances of these ports from the inlet of the pipe are , 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.

Figure 7(d) illustrates the averaged vapor volume fraction, , as a function of the vertical distance from the pipe wall, , 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 Figure 7(e) represent variations of as a function of 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.

Figure 7(f) shows variations of as a function of 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.

Figure 7(g) represents the final stage of the process which is an ongoing process of bubble formation and deformation. This figure depicts variations of as a function of 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 Figure 7 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].

Figure 7(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 the xz-plane in Figure 8(a). Here, the dimensionless pressure is defined as , where represents the local pressure. Figure 8(b) depicts variations of the dimensionless pressure (defined as ) as a function of 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.

(a)

(b)

(c)

(d)

(e)

(f)

Figure 8: Cavitation in a tube. (a) The instantaneous computed contours of the dimensionless pressure on the xz-plane bubbles in a cavitation flow of water at the cavitation number of . The locations of the sampling ports are the same as those in Figure 7(c). (b)–(e) Variations of the dimensionless pressure as a function of . (f) The instantaneous computed velocity field of a cavitating water flow through the pipe with inner diameter of  m at the discharge of  L/s.

Figure 8(c) illustrates variations of the dimensionless pressure, , as a function of 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 as a function of 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 with .

The squares, circles, diamonds, and left triangles in Figure 8(e) represent variations of as a function of 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.

Figure 8(f) shows the computed velocity vector filed in the xz-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 . 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 Figure 2.

3.3.4. Cavitation in a Spillway Chute

The previous section described the complexity of the occurrence of cavitation. As stated earlier, Figure 2 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. Figure 9(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°.

(a)

(b)

(c)

Figure 9: Cavitation in a spillway chute. (a) Schematic of a hypothetical hole with a diameter on the surface of a model of Kárahnjúkar spillway at a distance from the inlet of  m. The total length of the spillway is  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 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.

The inlet boundary condition at was set to  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 . Periodic boundary conditions were used in the -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]

where represents dimensionless, sublayer-scaled distance, and is the von Kármán constant (typically the value 0.41 is used).

The inset of Figure 9(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 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 . 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.

Figure 9(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.

Figure 9(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.

Figure 9(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 Figure 9(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. Figure 10 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.

(a)

(b)

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, = 2.42 m, = 0.44 m, = 0.11 m, = 0.11 m, = 2.15 m, = 0.022, = 0.55 m, = 20°, and = 30°.

By implementing periodic boundary conditions in the -direction, a small part of the system that is far from the vertical side walls was simulated. The model flow rate was set to = 20 m/s at the inlet with a size of = 0.11 m, corresponding to 2250 m³/s.

As can be seen from Figure 11, 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.

Figure 12(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 . Here, the velocity is normalized with the water velocity at inlet . 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.

(a)

(b)

(c)

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, = 0.62 m. (c) Air concentration distribution versus in the aeration region. Here, right triangles, squares, left triangles, and diamonds represent the air concentration at the sampling ports S₁, S₂, S₃, and S₄, respectively.

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 Figure 12, 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, 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 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.

Figure 13(a) illustrates the computed contours of the velocity field for a hypothetical case in which the discharge of the Kárahnjúkar spillway is 3150 m³/s. Figure 13(b) indicates that the air enters the ventilation zone with a normalized velocity of . 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.

(a)

(b)

(c)

(d)

Figure 13: The computed flow field for a hypothetical case with a discharge of  m³/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 as a function of . Squares, circles, diamonds, left triangles, gradients, and deltas represent the normalized velocity at the sampling ports S₁, S₆, S₅, S₂, S₃, and S₄, respectively.

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 m³/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.

Nomenclature

: | Archimedes number
: | Backscatter parameter
: | Drag coefficient of the bubble, which is a function of the Reynolds number of the bubble
: | Pressure coefficient
: | Smagorinsky constant, given as
: | Pipe diameter
: | Bubble diameter
: | Hole diameter
: | Force between two bubbles
: | Forces exerted by the liquid on the bubbles per unit volume
: | Height of aquarium stone
: | Step height
: | Gas constant
: | Gravitational acceleration
: | Second-order identity tensor
: | Kolmogorov constant
: | Length of the sheet
: | Distance of the hole from the inlet
: | Axial distance of the pressure measuring port for the pipe inlet
: | Pipe length
: | Total length of spillway chute
: | Length of aquarium stone
: | Length of the stone in the axial direction
: | Total length of the region covered by stones
: | Vertical distance from the pipe wall
: | Bubble number density
: | Dimensionless pressure defined as
: | Vapor pressure at ambient temperature
: | Constant in (1)
: | Ambient pressure in the fluid at a large distance from the bubble
: | Filtered pressure
: | Discharge
: | Bubble critical radius
: | Location of the th bubble
: | Distance between the centers of two bubbles
: | Bubble equilibrium radius
: | Bubble Reynolds number defined as
: | Ratio of the mesh spacing and the Kolmogorov length scale, so that
: | Bubble radius
: | Bubble interface velocity
: | Resolved rate of strain tensor
: | Source term in (5) to account for the mass transfer between liquid and vapor phase
: | Temperature
: | Time
: | Mixture velocity
: | Velocity of the th bubble with a diameter
: | Water velocity in absence of the bubble
: | Subgrid-scale part of based on Favre filtering
: | Density-weighted Favre filter velocity defined as
: | Slip velocity
: | Normalized velocity.

Greek Symbols

: | Chute angle for chute slope
: | Filter width
: | Constant in (3)
: | Surface tension
: | Dynamic viscosity of water
: | Dimensionless pressure defined as
: | Filtered density defined as
: | Bubble volume fraction
: | Vapor volume fraction defined as
: | Air volume fraction
: | Water volume fraction
: | Air density
: | vapor density
: | Water density
: | Cavitation number
: | Cavitation number
: | Subgrid scale (SGS) stress tensor.

Conflict of Interests

The author declares that there is no conflict of interests regarding the publication of this paper.

References

H. T. Falvey, Cavitation in Chutes and Spillways, U.S. Department of the Interior, Bureau of Reclamation, 1990.

W. Lee and J. A. Hoopes, “Prediction of cavitation damage for spillways,” Journal of Hydraulic Engineering, vol. 122, no. 9, pp. 481–488, 1996.

P. Novak, A. I. B. Moffat, C. Nalluri, and R. Narayanan, Hydraulic Structures, Spon Press, London, UK, 2004.

Landsvirkjun, Kárahnjúkar HEP Iceland, Physical Model Investigation on the Kárahnjúkar Dam Spillway, VAW, Laboratory of Hydraulics and Glaciology of the Swiss Federal Institute of Technology Zurich, 2006, http://www.vaw.ethz.ch/people/wb/archive/wb_251_karahnjukar_dam_spillway.

D. D. Joseph, “Cavitation in a flowing liquid,” Physical Review E, vol. 51, Article ID R1649(R), 1995.

L. G. Leal, Advanced Transport Phenomena, Cambridge University Press, Cambridge, UK, 2007.

P. Zamankhan, “Large eddy simulation and PIV experiments of air-water mixing tanks,” Communications in Nonlinear Science and Numerical Simulation, vol. 15, no. 6, pp. 1511–1525, 2010.

M. P. Allen and D. J. Tildesley, Computer Simulations of Liquids, Oxford University Press, Oxford, UK, 1987.

G. H. Sauer and J. Schnerr, “Physical and numerical modeling of unsteady cavitation dynamics,” in Proceedings of the 4th International Conference on Multiphase Flow, New Orleans, La, USA, 2001.

P. Zamankhan, “Solid structures in a highly agitated bed of granular materials,” Applied Mathematical Modelling, vol. 36, no. 1, pp. 414–429, 2012.

P. Zamankhan, “Formation patterns at the air-grain interfaces in spinning granular films at high rotation rates,” Communications in Nonlinear Science and Numerical Simulation, vol. 17, no. 11, pp. 4224–4240, 2012.

P. Zamankhan, “Simulation of self-segregation of a low density spherical particle in a bubbling bed,” Chemical Engineering Journal, vol. 181-182, pp. 842–845, 2012.

P. Zamankhan, “Sinking and recirculation of large intruders in vertically vibrated granular beds,” Advanced Powder Technology, vol. 24, no. 6, pp. 1070–1085, 2013.

C. Chaux, L. Duval, A. Benazza-Benyahia, and J.-C. Pesquet, “A nonlinear Stein-based estimator for multichannel image denoising,” IEEE Transactions on Signal Processing, vol. 56, no. 8, pp. 3855–3870, 2008.

http://labs.autodesk.com/technologies/123d/

O. Coutier-Delgosha, R. Fortes-Patella, and J. L. Reboud, “Evaluation of the turbulence model influence on the numerical simulations of unsteady cavitation,” Journal of Fluids Engineering, vol. 125, no. 1, pp. 38–45, 2003.

H. L. Stone, “Iterative solution of implicit approximations of multidimensional partial differential equations,” SIAM Journal on Numerical Analysis, vol. 5, pp. 530–558, 1968.

C. Antoine, “Tensions des vapeurs, nouvelle relation entre les tensions et les températures,” Comptes Rendus des Séances de l'Académie des Sciences, vol. 107, p. 681, 1888.

V. M. Semenkov and V. O. Saranchev, “Computer program for hydraulic calculation of structures for controlling aeration of the flow on spillways,” Hydrotechnical Construction, vol. 30, no. 5, pp. 234–242, 1996.

R. S. Gal'perin, A. G. Oskolkov, V. M. Semenkov, and G. N. Tsedrov, Cavitation on Hydraulich Structures, Énergiya, Moscow, Russia, 1977, (Russian).

P. Zamankhan, “Rollers in low-head dams—challenges and solutions,” Communications in Nonlinear Science and Numerical Simulation, vol. 17, no. 12, pp. 5273–5285, 2012.