NUMERICAL SOLUTION OF VISCOUS 2DV FREE SURFACE FLOWS: FLOW OVER SPILLWAY CRESTS

 

MARTIN F. BÜRGISSER

 

Peter Rutschmann

 

Laboratory of Hydraulics, Hydrology and Glaciology (VAW)

Swiss Federal Institute of Technology Zürich (ETHZ)

CH-8092 Zürich, Switzerland

Tel: +41-1-632 5103, Fax: +41-1-632 1192, e-mail: buergisser@vaw.baum.ethz.ch

 

 

ABSTRACT

An iterative solution to the incompressible two dimensional vertical (2DV) steady Navier-Stokes Equations (NSE) is used to solve for the a priori unknown free surface, the velocity and the pressure fields. The entire waterbody is covered by a unstructured finite element grid which is locally refined. The solution method is tested for different discharges on two standard spillway geometries. The results agree with available experimental data.

 

Keywords: Numerical modelling, Navier-Stokes equations, free surface flow, Finite Element Method, design of hydraulic structures, flow over weirs, spillway geometries.

 

INTRODUCTION

Free surface flows are encountered in hydraulic engineering problems including water jets, weirs and around gates. For the numerical modelling of a steady flow, the area of calculation and therefore the position of the free surface have to be known. To circumvent the difficulties attributed to this additional non-linearity, numerical simulations are often made with the Laplace equation, a subset of the incompressible NSE by neglecting viscosity and vorticity, and solving for the stream function. The profile of the free boundary becomes the other independent quantitiy subject to variation. Such numerical results using the finite element method are available (e.g. Ikegawa & Washizu 1973, Diersch et al. 1977, Dao-Yang & Man-Ling 1979 and Li et al. 1989). The shape of the strongly curved spillway crest was in all cases discretized only with few gridpoints.

In the present work, a vertical plane of the waterbody was modeled using the incompressible 2DV-NSE. The equations were solved for the velocities and pressure fields simultaneously. The position of the free surface was initially assumed and a locally refined unstructured grid then generated. At the free surface the zero pressure boundary condition was assumed. The calculated velocities on the free surface may have normal components which were used for the correction of the boundary shape. A steady solution was attained when the calculated normal velocity components vanished along the entire free surface. Due to node adaption, boundary elements deform significantly and overlap with other elements. To avoid inappro­priate elements, inner nodes were also corrected, either by spring algorithm or by generating a completely new grid.

This approach for calculating free surface flows is general and results for velocity components and pressure for each gridpoint are readily accessible. A benefit of the unstructured grid is that the calculation area can be complex, and a locally refined grid can be used in areas of interest. Flows with more than one deformable surface, such as water jets (Bürgisser 1998a) or bed erosion can be handled simultaneously. The approach can be easily extended to solve for transient flow situations, where the grid follows the motion in an Lagrangian way and a coupling with turbulence modelling is possible with minor changes (Bürgisser 1998b).

 

GOVERNING EQUATIONS

The governing equations for free surface water flows are the isothermal, incompressible viscous NSE. For engineering purposes, turbulent fluctuations are removed by time averaging (Reynolds equations) and the eddy viscosity concept was applied (ASCE 1988). The present results were performed with a constant eddy viscosity for the entire calculation area. The steady state momentum equations for the component directions i=1 and 2 are

, (Eq. 1)

which are solved simultaneously with the continuity equation for the velocity components u1 and u2 and pressure p. The external body force is the acceleration due to gravity g in the vertical direction z of the Cartesian coordinate system. The effective eddy viscosity includes the molecular and the turbulent viscosities, i.e. the flow may be rotational (Bürgisser 1998b).

The Finite Element Method was used for the algebraic representation of the governing equations, which is appropriate for unstructured, locally refined grids. To overcome oscillation problems, which are typical for Eq. (1) when convection terms are dominant, the Galerkin/least-squares upwind technique was used (Hughes et al. 1989). This method is stable for high Reynolds numbers and allows solution of the basic variables with equal order (linear) shape functions. To increase numerical stability and to model a part of turbulent effects, was set constant to 0.001m²/s for the present calculations.

 

FREE SURFACE ADAPTION

The iterative approach for the free surface adaption can be divided into two steps. In a first step the spatial domain is discretized assuming an initial position of the free surface. At the air-water interface, physically-correct boundary conditions include atmospheric pressure and tangential velocities. Because the shape of the surface is initially a rough approximation, only one of the two boundary conditions was incorporated (relative atmospheric pressure p=0), leaving the free surface unconfined. In the second step, the governing equations were solved for pressure and velocity. If the position of the assumed free surface was incorrect and the calculated velocities were not parallel to the free surface, its position was vertically moved and repredicted using the normal velocity components for correction (Fig. 1a). This process was repeated until convergence, that is when both boundary conditions were satisfied.

The interior nodes can be moved as well, but it is often easier to deform only the boundary elements. When the mesh is too distorted, remeshing is necessary and the velocities are interpolated from the original grid to maintain accurate values for the linearization of the momentum equations.

 

TEST CASE: FLOW OVER SPILLWAY CRESTS

The iterative solution method was tested for two different spillway geometries. The first corresponds to the WES three-arc curve upstream quadrant and a power function in the downstream quadrant (USCE 1970), whereas the second shape has a continuous change of curvature over the entire crest (Hager 1987). The numerical results were compared with experiments of Hager (1991) performed in a 0.5m wide and a 0.7m deep channel. These geometries were used for the numerical setup; the first spillway had a design head of H_D=0.1m and a downstream chute slope of 45°, whereas the second spillway had a design head of H_D=0.2m and a spillway slope of 30°. The input data and some results are summarized for the test cases in Table 1. The calculated surface profile of the second spillway compares well with experiments (Vischer & Hager 1998).

At the solid boundary, a slip velocity condition (v_n=0) was applied and the horizontal approach velocity was calculated from discharge Q for H/H_D=1 (Fig. 1b). This boundary condition was corrected after each change of the inflow height h₀. The outflow boundary remained normal to the bottom slope for all iterations. The free surface was initially approximated and then corrected within 100 iteration steps.

The largest gradients for all variables occur at the beginning of the upstream quadrant (Fig. 2a-c). An extremely fine grid was generated along the crest to model the different geometries as accurately as possible (Fig. 2d). The surface elevation and the crest velocity (Fig. 3a,b) are well predicted. The calculated underpressure (Fig. 3c) is larger than measured. The largest values result for the largest curvature of the crest in the upstream quadrant with only few data available. In this region the largest risk of cavitation appears and the shape of the crest could be improved.

 

CONCLUSIONS

Flexible solution methods for steady free surface flows problems are presented. By solving for velocity and pressure fields, the results can be readily compared with data from physical model tests. The close agreement between the model predictions and the laboratory results suggests that the present approach is worth further development. Local grid refinement allows investigations of hydraulic details, which are difficult to observe experimentally.

 

REFERENCES

ASCE, Task Committee on Turbulence Models in Hydraulic Computations (1988). "Turbulence Modeling of Surface Water Flow and Transport: Part I." J. Hydr. Engrg., 114(9), 970-991.

Bürgisser, M. F. (1998a). "Hydrodynamics of 2DV Free Surface Flows." Hydro­informatics' 98, 49-53.

Bürgisser, M. F. (1998b). "Numerische Simulation der freien Wasseroberfläche bei Ingenieurbauten," Mitteilung, Versuchsanstalt für Wasserbau, Hydrologie und Glaziologie (VAW), ETH Zürich (in German).

Dao-Yang, D., and Man-Ling, L. (1979). "Mathematical Model of Flow over a Spillway Dam." Teizième Congrès des Grands-Barrages, Q 50, R 55, 959-976.

Diersch, H.-J., Schirmer, A., and Busch, K.-F. (1977). "Analysis of Flows with Initially Unknown Discharge." J. Hydr. Div., ASCE, 103(HY3), 213-232.

Hager, W. H. (1987). "Continuous Crest Profiles for Standard Spillway." J. Hydr. Engrg., 113(11), 1453-1457.

Hager, W. H. (1991). "Experiments on Standard Spillway Flow." Proc. Instn. Civ. Engrs, 91(2), 399-416.

Hughes, T. J. R., Franca, L. P., and Hulbert, G. M. (1989). "A New Finite Element Formulation for Computational Fluid Dynamics: VIII. The Galerkin/Least-Squares Method for Advective-Diffusive Equations." Comp. Meths. Appl. Mech. Engrg., 73, 173-189.

Ikegawa, M., and Washizu, K. (1973). "Finite Element Method Applied to Analysis of Flow Over a Spillway Crest." Int. J. Num. Meth. Engrg., 6, 179-189.

Li, W., Xie, Q., and Chen, C. J. (1989). "Finite Analytic Solution of Flow over Spill­ways." J. Engrng. Mech., 115(12), 2635-2648.

USCE (1970). "Hydraulic Design Criteria.", US Corps of Engineers, Army Waterways Experiment Station, Vicksburg MI.

Vischer, D. L., and Hager, W. H. (1998). Dam Hydraulics, John Wiley & Sons: Chichester, New York.

 

	H/HD [-]	CD [-]	Q [m3/s]	H [m]	u0 [m/s]	CD [-]
	input values			calculated values
USCE profile	1	0.495	0.0694	0.0988	0.0869	0.501
HD=0.1 [m]	2	0.545	0.2159	0.1897	0.2427	0.576
Hager profile	1	0.495	0.1961	0.1925	0.2197	0.514
HD=0.2 [m]	2	0.535	0.5995	0.3805	0.5549	0.543

 

Table 1:

Input values for the selected testcases. The input value for C_D is from Hager (1991). The head on the spillway H and the horizontal inflow velocity u₀ are numerical results. The resulting coefficient of discharge is calculated as C_D=Q/[2g(H+(u₀²/2g))³]^1/2.

 

 

Fig. 1 Definition sketches: a) Free surface nodes, b) Initial shape for free surface calculation.

 

 

Fig. 2: Variable fields for Hager profile with H/H_D=1: a) horizontal, b) vertical velocity components; c) pressure and d) calculation grid with refinement at spillway crest.

 

 

Fig. 3: Normalized figures: a) free surface profile, b) crest velocity, c) bottom pressure. Symbols are for observations, lines for computations.