·
·
(1)
Larry Weber1, Heqing Huang2 and Yong Lai3
1Assistant Professor, Civil & Environmental Engineering, Iowa Institute of Hydraulic Research, the University of Iowa, 300 S. Riverside Dr., Iowa City, IA 52242-1585
319-335-5597 (ph); 319-335-5238 (fax); E-mail: larry-weber@uiowa.edu
2Research Assistant, IIHR, E-mail: heqingh@engineering.uiowa.edu
3Research
Engineer, IIHR, E-mail: yong-gen-lai@uiowa.edu
Abstract:
It is found in 1960s that excessive gas
supersaturation happened downstream of dam spillway will lead to lethal bubble
disease for precious fish species. EPA in 1986 has set the total gas
supersaturation standard as 110% for better environment not only for man, but
for fish also. In USA today, there is an urgent need to investigating ways to
bring down gas supersaturation. This numerical model is developed to meet the
need and cutting experimental cost. Based on the recent advances in the fields
of air bubble entrainment, numerical gas saturation, supersaturation and flow
simulations, the method for building 2D and 3D models are discussed in the
paper. Preliminary simulation results shows the proposed comprehensive numerical
model can well simulate the gas supersaturation downstream of dam spillway, as
well as suggest possible methods of solving problem by parameter studies.
Keywords: gas supersaturation, bubble disease, water quality, numerical simulation
Generally, a spillway is built in accompany to a dam to protect it from excessive flooding. It has been found since 1960s that spillway flow will result in gas, especially nitrogen gas supersaturation in water, and this gas supersaturation causes a fatal gas bubble disease for fish and endangers their existence (Ebe, 1969, 1975; Weitkamp and Katz, 1980). EPA (1986) has established water quality standards for dissolved gases at 110% of saturation. Much research has been done since then in the US to reduce the gas supersaturation around dam stilling basin region to meet the standard.
Deflectors have been designed and tested both in laboratory experiments and in the field to alter the flow dynamics and lower the gas concentration. (Johnson and Dawley, 1974; Geldert, 1998). It has been proved an effective method for decreasing the gas supersaturation and performance curve method has been developed for its design (Weber, Haug and Jeske, 2000).
To acquire highly reliable prediction from models, numerous researches have been done about various basic respects in related fields. Chanson (1996) has summarized investigations about air entrapment in a number of typical flow conditions. What he concludes that for most natural flows the flow velocity profiles under air entrapment is still continuous and are almost of the same shape as those of no air entrapment, entrapped air concentration profiles are also continuous and show sign of convective-diffusion control process may be said is the experimental basis for our numerical simulation.
Moog and Jirka (1999) have researched the air-water gas transfer in uniform and non-uniform open channel flows. The way they suggested in their paper about how to decide gas transfer coefficient may be quite helpful for numerical models’ parameter decision.
Gulliver etc. (1997) have done work for measurement of Effective Saturation Concentration (ESC) of gas in the plunge pools at hydraulic structures. Hibbs and Gulliver (1998) have given formulas for prediction of ESC at spillway plunge pools. Basically in their method, first the maximum depth of bubble penetration in unconfined pool is calculated by the spillway geometry and inlet flow conditions. Then the effective bubble depth is determined by considering the confinement of stilling basin constraint and bubble swarm shape to account for the static pressure increase. Finally an average ESC in the stilling basin is predicted by incorporating the transfer efficiency also. Geldert etc. (1998) pointed out the four predominant factors resulting in supersaturation in the stilling basin: turbulent mixing, air-water surface increase due to entrapped air bubbles, residence time of bubbles and increase in hydrostatic pressure. Formulas are also given in their paper for predicting the ESC by incorporating the air-water surface, bubble-water gas transfer mechanisms.
Physical models are generally made to keep the Froude number as a constant for open channel flows. It is difficult to maintain Reynolds and Weber number constant at the same time for accurate scale study. This is a limitation in physical model studies of bubble flows (Ettema, 2000).
Further, Orlins and Gulliver (2000) formulated a primitive 2-D numerical model using the hydrodynamic information such as flow velocity, turbulent fluctuating velocities, and bubble distribution data collected from the hydraulic model study (Mannheim and Weber, 1998) to predict downstream gas supersaturation.
Their 2-D numerical model requires inputs of both hydrodynamic data and information about the distribution of air bubbles in the flow, along with the depth-averaged concentration of TDG upstream of the spillway.
We feel that the following limitations to Orlins and Gulliver’s numerical model can be improved:
l Requiring large quantities of laboratory experiment data input such as flow and turbulent quantities.
l Applicable only to simple geometry because of its the single block and orthogonal character.
l Two dimensional, not accounting for transverse changes.
l Boundary conditions are not clearly discussed.
l Model source parameter β1 and surface transfer parameter β2 leaves much space for human manipulation.
The objective of the our research is try to solve those problems by building a integrate numerical model combining numerical flow simulation and numerical supersaturation simulation at there dimensional level with programming based on multi-block non-orthogonal grid, and lessen the adjustable parameters by associating transfer parameters with relative flow quantities from numerical flow simulation.
The flow around the spillway is essentially turbulent flow. For turbulent flow, theoretically the best numerical approach should be DNS or LES. However the fine grid needed for this engineering purpose is far above the capability of present computer (Ferziger and Peric, 1999). There are numerous reports that RANS with k-e or k-w model can simulate the average open channel flow well (Messelhe and Odgaard, 1997; Messelhe and Sotiroplos, 2000; etc.). Therefore it is to be used in the proposed comprehensive numerical simulation.
Upon accomplishment of the proposed numerical model, it shall be possible to predict the highest supersaturation, the dimension of high concentration area, the evolution of high concentration plume when the spillway gates are opening or being closed (unsteady simulation). It shall provide valuable reference for new dam spillway or deflector design and cut experimental cost.
Here we summarize the main assumptions used in the comprehensive model based on the experimental results by Chanson (1996), so that the capabilities and limitations of the model can be fully understood.
(1) High velocity air-water flows behave as a homogeneous mixture.
(2) Mean air-water velocity is mostly the same as non-aerated, mono-phase flows. The air entrainment has no direct effect upon the mean velocity distributions. This suggests that we may simulate the two physical processes separately, and then combine them for complete gas concentration information.
(3) Air content and mean air-water velocity distributions are smooth, continuous and differentiable functions. Though there exist air bubbles, the continuum assumption for Navier-Stokes equations remains intact. It is possible to get good averaged flow simulation from RANS.
(4) The distribution of dissolved gas concentration is continuous and can be represented by a advection-diffusion theory.
The generic conservation equation in integral form is:
·
·
(1)
where W = space considered, S = the surface of W, q = the source of the species concentration C, and D = the diffusion coefficient. This equation is the basis of the finite volume numerical method. Expressed in Cartesian coordinates in differential form:
(2)
where C = dissolved gas concentration; Ceq = equilibrium concentration accounted already for hydrostatic pressure at the bubble-water interface ; U,V, W = flow velocity in the three coordinate directions respectively, Dx, Dy, Dz = turbulent diffusivity of gas in water in the three coordinate directions respectively, K1 = the mass transfer coefficients across the internal bubble-water interface, a = specific surface area.
When steady state concentration is of interest, implicit Euler method is going to be used with large time step to get steady solution. When unsteadiness is going to be investigated, the implicit, second order accuracy three-time level method is going to be used.
The velocity components U, V, W in the three convective terms are to be obtained through k-ε model RANS simulation with free surface.
In 2D simulation for which the results is presented in this paper, the diffusion coefficient is to be determined from the RMS vertical velocity fluctuation as suggested by Nezu and Nakagawa (1993). The RMS vertical velocity fluctuation n’ is estimated from the turbulent kinetic energy k given by RANS k-ε simulation as:
(3)
For 3D simulation, Diffusion coefficients Dx, Dy, and Dz are to be determined from the vorticity acquired from RANS simulation to account for the heterogeneity of diffusion coefficient as equations (4) shows. It will be discussed in detail in out next paper for presenting 3D simulation results.
(4)
The finite volume method for non-orthogonal grid presented in the newly published book “Computational Methods for Fluid Dynamics” (Ferziger and Peric, 1999) is used to discretize the conservation equation (1). It is summarized here briefly for completeness.
For convective fluxes:
·
(5)
where for example, at a
3-D control volume’s east face, its mass flux me should be
expressed as
(6)
where Si means the grid cell east face’s projection area in the plane perpendicular to i axis. Ci should be the value at the cell face, which can be obtained by linear interpolation.
For diffusive fluxes, Gaussian
theorem is used first, the gradient of C is obtained by
(7)
which demonstrates that the gradient of C with respect to x at the CV center can be obtained by summing the products of C with the x-components of the surface vectors at all faces of the CV and dividing the sum by the CV volume, then a deferred correction method is adopted to eliminate the oscillatory solution:
For example at the CV e face
(8)
This indicates that the angle between the cell face normal n and the line x connecting the cell centers on either side is important.
After discretization, equation (1) shall be expressed as a linear equation system for all grid points, then SIP solver is employed to find the final solution.
There are five type boundary conditions to be considered:
(1) Inlet: constant Dirichlet type. It may be a constant profile from dada available.
(2)
Outlet: If the Peclet Number (
) is sufficiently large, no boundary condition is
needed at an outlet boundary (Patankar, 1980)
(3)
Symmetry: Neumann type as
, where n is the outlet normal direction. This is mainly used in flow simulation
for free surface.
(4) Constant concentration Wall: are to be used for free surface. We think it is more reasonable to assume ample exchange of dissolved gas with air at the surface where the pressure is atmospheric according to Henry’s law. If it is considered a sink around the free surface as Orlins and Gulliver’s model does, all the cell adjacent to the surface is considering leaking supersaturation. It is not physically true though both can give similar results. Therefore we delete it from the source term in equation (2) of Orlins and Gulliver’s model and include it here as boundary condition.
(5) Concentration adiabatic wall: are to be applied to riverbed.
Figure 1 shows the 2-D flow simulation uv vectors and turbulent kinetic energy contours for 3 cases: spillway with no deflector, low deflector and high deflector. It can be seen that the average flow field given by RANS with k-e model simulation catches most of the main features revealed by experiment (c.f. Fig.2, 3 of Orlins and Gulliver’s paper, 2000).
Vertically uniform bubble distribution and an exponential decay with distance used by Orlins and Gulliver (2000) is also used:
(9)
Pbub = probability of bubbles occupying a specific location; P0 = initial probability of a bubble at a given elevation as flow enters the stilling basin; w = bubble rise velocity; x = horizontal distance from the entrance of the stilling basin; and q = specific discharge. Bubble rising velocity is assumed to be at constant 0.25 m/s. The air-bubble water contacting area decreases as bubbles rising and merging can also be reflected through define x to be the vertical distance to spillway jet surface
The saturation concentration Cse in equation (2) for turbulent jet area is approximated by the method suggested in Hibbs and Gulliver’s paper (1997). For the transfer efficiency coefficient needed, it may be treated as a calibration parameter. Here 0.4 is taken from model calibrations.
With the above parameter setting, the time evolution numerical supersaturation simulation results are show in Figure 2 for high-position deflector case. The concentration shows in the figure is the percent above upstream inflow saturation. Figure 3 shows the comparison of Orlins and Gulliver’s Model and ours with field data. Good consistency shows the feasibility of this combined numerical simulation.
For no deflector case, it shows the highest supersaturation averaged about 120% and extended far downstream of the spillway. For low-position deflector case, since the deflector has decreased the jet penetration depth significantly, and air bubbles will rising to the surface faster. Besides, higher diffusivity in the region will also promote dissolved gas return to air through the surface. Therefore, spillway with low-position deflector gives a low gas supersaturation average about 110%. For high-position deflector case, spillway inflow is a surface jet. It shows the lowest supersaturation degree downstream.
Model parameter study may reveal some interesting characters as the following example shows.
Figure 4 show the 3 simulations for case 1 without deflector on spillway. Only assumed uniform gas diffusivity D is different, all other conditions are the same. It can be seen that when the diffusion coefficient is big (D=10 m2/s), not only the supersaturation degree is lower in the stilling basin region, but also the area of supersaturation is largely limited in the small area around the stilling basin. This suggests possible ways of containing the supersaturation by making the dissolved gas diffusivity big and uniform. As the diffusion coefficient of saturated gases grows smaller, the high concentration area grows bigger and extends further downstream of the spillway. When it is very small as D=0.1 m2/s case show, the contours of supersaturating even shows the influence of flow pattern.
Even with the 2-D model, it has shown encouraging results. Next step, it will be extended to three dimensional model with multi-block non-orthogonal grid taken into account all main river features: uneven river bed, ten dam powerhouse exits, and twelve spillway gates and their possible deflectors, so that any combination of the real-world flow of the dam and its gas supersaturation can be simulated. Figure 5 shows CFD uv contours and flow velocity vectors for the downstream reach of Wanapum Dam of interest. Using methods discussed above, especially equation (4) for diffusion coefficient in different direction, and by defining the spillway jet dimension which is the principal source of gas supersaturation, the supersaturation contours for the river reach can be predicted for any cross section for better spillway and deflector design. Those will be presented in our next paper soon.
References
Chanson, H. (1996). Air bubble entrainment in free-surface turbulent shear flows. Academic Press, UK.
Ettema, R. (2000) Hydraulic Modeling—Concepts and Practice. ASCE Manuals and Reports on Engineering Practice No. 97.
Ferziger, J. H. and Peric, M. 1999. Computational Methods for Fluid Dynamics. Springer.
Geldert, D.A., Gulliver, J.S. and Wilhelms, S.C. (1998) “Modeling Dissolved gas supersaturation below spillway plunge pools.” J. Hydraulic Engrg, vol. 124, No.5, 513-521.
Gulliver, J.S., Hibbs, D.E. and McDonald, J.P., (1997) “Measurement of Effective Saturation Concentration for Gas Transfer” J. Hydraulic Engrg, vol. 123, No.2, 86-97.
Hibbs, D.E. and Gulliver, J.S. (1997) “Prediction of Effective saturation concentration at spillway plunge pools”, J. Hydraulic Engrg, vol. 123, No.11, 940-949.
Mannheim, C. and Weber, L.J. (1998). Dissolved Gas Supersaturation Downstream of a Spillway. I: Physical Model. Submitted for review to J. Hydraulic Research.
Messelhe,E.A. & Odgaard,A.J. (1998). 3D Numerical Flow Model for Fish Diversion Studies at Wanapum Dam. J. Hydraulic Engineering. Vol.124, No.12, 1203-1214.
Meselhe, E.A. and Sotiropoulos, F. (2000). Three-dimensional numerical model for open-channels with free-surface variations. J. Hydraulic Engineering. Vol.38, 2000, Vol.125, No.1. 115-121.
Moog, D.B. and Jirka, G.H. (1999) “Air-water gas transfer in uniform channel flow” J. Hydraulic Engineering. Vol.125, No.1. 3-10.
Moog, D.B. and Jirka, G.H. (1999) “Stream reaeration in nonuniform flow: Macroroughness enhancement” J. Hydraulic Engineering. Vol.125, No.1. 11-16.
Nezu, I.,
and Nakagawa, H. (1993). Turbulence in open-channel flows, Balkema,
Rotterdam, the Netherlands. pp347-359.
Orlins, J.J. and Gulliver, J.S. (2000). Dissolved Gas Supersaturation Downstream of a Spillway. II: Computational Model. J. Hydraulic Research, Vol.38, No.2, 151-159.
Fig. 1 2-D flow simulation vector and u velocity contour plot.
Fig. 2 The comparison of filed data with Orlins-Gulliver model and ours.
Fig. 3 Supersaturation contours for no deflect case with three different uniform diffusion coefficients.
Fig. 4 The uv contour and flow velocity vector plot at near free surface from 3-D RANS k-ε simulation for Wanapum Dam downstream reach.