Katja Rettemeier1, Olaf Bergen2, Andreas Van Linn1 and JÜrgen Köngeter1
1Institute of Hydraulic Engineering and Water Resources Management,
University of Technology Aachen, Germany
2IMS Ingenieurgesellschaft mbH, Hamburg, Germany
Mies van der Rohe Str. , 52056 Aachen, Germany,
Phone: ++49-(0)241-805263,
Fax: ++49-(0)241-8888-348, E-mail: rettemeier@iww.rwth-aachen.de.
Abstract: The paper shows the application of a finite element model to simulate the flow and turbulence characteristics of very complex natural geometry. The flow phenomena in natural geometry such as lakes and reservoirs are very complex and in most cases fully three-dimensional. This is due to the irregular boundaries, the large volume of the water body and the variability of the external forces on the water body. Numerical models must at least be able to simulate the large scale processes in correspondence of the influence of wind forces, momentum induced by in- and outflow and the coriolis force. Another problem is the large spectrum of turbulence existing in large natural geometry. Standard statistic models based on the Reynolds-equations are widely used, but it is very difficult to model the whole spectrum correctly which depends on the individual reservoir. Additionally strong anisotropic flow characteristics demand an anisotropic turbulence model. Direct Numerical Simulation would require more computer capacities than available thus a very promising approach is the Large-Eddy-Simulation technique. This paper presents the numerical method which is based on a Finite Element model in conjunction with a mixed model LES technique. The results of the validation study on flow phenomena in an idealized lake are discussed as well as the turbulence characteristics in large irregular water bodies.
Keywords: CFD, LES, finite element method, reservoir flow, turbulence
Lakes and reservoirs are sometimes subject to contamination. The most common case is the inflow of nutrients that often lead to a strong algae growth and disturb the ecological balance of natural flow systems. Furthermore accidental contamination may occur and challenges any kind of alert systems in environmental hydraulics. Numerical models can be a very helpful tool to develop a better understanding of limnological processes, and so minimize negative impacts by improved flow system management. Furthermore they can be used to evaluate and optimize the efficiency of the sometimes very high investment costs especially in the field of reservoir design.
The flow phenomena in lakes and reservoirs are very complex and in most cases fully three-dimensional. Although they have been studied for several decades (Imberger & Patterson, 1990) only nowadays with increasing computer capacities and improved measurement techniques a more detailed analysis is possible. A trustworthy numerical model must be able to simulate the influences of density dependent flow, wind forces, acting at the water surface, momentum induced by in- and outflow and the coriolis force due to Earth rotation. To fulfill these objectives the numerical model must be able to simulate the large scale processes as well as the entire spectrum of turbulence existing in large natural geometry in full detail. A very promising approach is the Large-Eddy-Simulation technique. The numerical model applied is based on the finite element method and contains a three-dimensional ‘Large-Eddy’ approach. To account for the anisotropy of the flow an anisotropic mixed model is applied in conjunction with a moderate anisotropic Smagorinsky-model.
The first step for the validation of the numerical model was the separate study of the influencing factors, thus the model could be validated with simple test cases before applying it to the complex structures of lakes and reservoirs (Rettemeier et al., 2000). In the sense of a complete validation of the LES model the application to a real lake would be desirable. This is however only reasonable if an extensive array of flow and turbulence data is present. However the available data, which permits a quantitative comparison between the natural and the simulated current is not sufficient so far. Additionally a minimization of the topological influence enables a better understanding of the influence of the main acting forces and the topography on the flow pattern. Thus an idealized lake was studied. Special attention is drawn towards the turbulence characteristics evaluating energy transport and spectral analysis which agreed well with theoretical investigations. A detailed analysis of the flow phenomena in a lake as well as the turbulence characteristics will be presented.
The complexity of the flow phenomena in large-scale geometry is due to the very irregular boundaries inducing high velocity gradients and a large spectrum of turbulent flow. Therefor the model must be capable of representing the whole eddy spectrum and simulating three-dimensional flow in the irregular boundaries. A three-dimensional “Large-Eddy” approach fulfills these requirements.
The “Large Eddy” finite element model PASTIS-3D is based on the spatially filtered continuity and Navier-Stokes equations (eq.1 - eq.2):
and
(1)
with the subgrid scale stresses
(2)
Here
represents the velocity components
in x,y,z-direction, P is the kinematic pressure and
is the kinematic viscosity.
includes all internal forces. The
overbar denotes a spatially filtered variable.
All subgrid scale stresses are modeled with an anisotropic mixed model / scale similarity model as discussed by Aldama (1993) and used by Cottet & Wray (1997). To account for all subgrid scale stresses, an anisotropic mixed model in conjunction with an anisotropic Smagorinsky-model (4) was used (D = filterwidth, cs Smagorinsky-constant).
(3)
with
Finite hexahedron elements with arbitrary shape are used for the spatial discretization. Therefor all integration on the element level are obtained by full 2x2x2 Gauss point integration. The approximation functions are trilinear for the velocities and constant for the pressure respectively (FORKEL, 1995). For the time discretization the semi-implicit trapezoid scheme is used accounting for no numerical diffusion and time-accurate solutions are obtained by applying the very efficient ‘projection 2’ algorithm which was developed and analyzed by GRESHO (1990) and GRESHO & Chan (1990). For more details of the numerical approach see DANIELS (1992), FORKEL (1995) and BERGEN (1999).
In the sense of a complete validation of the LES model the application to a real lake would be desirable. This is however only reasonable if an extensive array of measurement data is present. However the available data, which permits a quantitative comparison between the natural and the simulated current for LES is not sufficient so far. Additionally a minimization of the topological influence enables a better understanding of the influence of the main acting forces and the topography on the flow pattern. Thus an idealized lake was modeled. The selected lake geometry forms a regular octagon, which approximates a circle. Figure 1 presents the discritization and the main dimensions. The horizontal grid spacing corresponds to 0.5 of the local flow depth so the full spectrum of turbulence can be modeled (Bergen, 1999).

Fig. 1 Discritization of the Lake Like Geometry
The initial condition are set to zero for the whole domain. A constant wind force corresponding to 5 m/s wind velocity acts on the water surface. On the west border of the lake on the symmetry axis an inflow boundary condition of constant velocity (0.1 m/s east) is applied over 12.0 m times 0.5 m thus the total discharge equals 0.6 m³/s. The outflow with the same dimension is located on the east border of the lake. The 2 Point Method (Schweim et al., 2000) is applied to the bottom boundary condition of the lake. The wall distance is set to 5 % of the local depth. The simulation period extents over three days with 25 s time steps. For the evaluation of the flow phenomena and the turbulence characteristics 50,000 s are considered. The LES model parameter cs equals 0.16, while the filter width corresponds to the local grid spacing (Bergen, 1999).
Figure 2 presents the simulated averaged velocity field at the water surface (average time 50,000 s). The contour plot represents the x-velocities while the direction of the flow is clarified by vectors. The influence of the Coriolis force, which causes the clear right diversion of the incoming jet as well as the diversion of the surface drift in east-southeast direction, are particularly remarkable. Due to the limited spatial expansion of the lake the surface drift with respect to the wind direction remains limited to an angle of » 25°. Thus it is clearly smaller than the Ekman diversion of 45°, determined on highly simplified assumptions. This corresponds with numerous field observations, which prove clearly reduced surface drift likewise.

Fig. 2 Simulated averaged velocity field – surface
High flow rates can be observed in particular in the northern, eastern and southern bank areas. Here the accelerated current is constricted by both the reduced flow depth and the lateral area bound. The distinct recirculation zone, which is located right before the outflow is particularly characteristic. Its formation can be described clearly on the basis of the mainstream paths in the lake, represented in figure 3. The mainstreams at the water surface are rerouted at the lateral bounds in the eastern third of the lake and funneled toward the outflow. Since both currents face each other, they collide after the plunging below the recirculation zone. This process is linked with high vertical velocities and strong turbulence. This can be observed clearly at the water surface by the interrupted wind drift. The currents do not collide frontally because they both have a velocity component in western direction. Thus local velocities against the wind direction occur at the surface. The main back flow follows the west northwest direction.

Fig. 3 Main flow paths in the water body of a lake
Figure 4 provides an
impression of the resolved, turbulent eddy scales. A snapshot of the transient
flow field is plotted after the simulation duration of two days. The averaged
velocities have been subtracted. The horizontal cross section at the middle of
the flow depth represents the velocity vectors over a contour plot of the
vorticity of the horizontal velocity field. The directly resolved eddy scales of
the turbulent flow are clearly visible in the lake. These eddies have the
highest energy and they develop within the area of the recirculation zone. They
are transported with the back flow from the recirculation zone towards the
northwest bank. By the process of vortex stretching the eddies pass a part of
their energy to small eddy scales. This process is connected with a reduction of
the vorticity of the large eddies. This can be detected clearly by the route of
transportation evident in figure 4.

Fig. 4 Horizontal snapshot of the fluctuation and vorticity of the velocity field
Thus the water body is divided into a strongly and weakly turbulent area by a virtual diagonal. It is to be noted that an estimation of the turbulence with the Reynolds number would supply an exactly opposite arrangement of the turbulence intensities. Therefor the consideration of the Reynolds number for the estimation of the turbulence degree of lake currents must be analyzed at least critically.
Besides the evaluation of the results with respect
to energy transport spectral analysis provides additional validation data of
the developed LES approach.. Figure 5 (left) shows the three-dimensional frequency
spectra for a FE node in the center of the lake. For the maximum lake dimension
of 500 m and the average, convective velocity at the regarded FE node of
m/s a characteristic frequency
of
Hz can be detected with Taylor’s hypothesis on frozen turbulence.
This frequency marks the area of the largest, geometry induced eddies and the
energy input into the system. The energy cascade has the characteristic
gradient, which changes over
- into a
gradient. It is not to be assumed
that the dissipative scales are actually described by the
gradient. On the contrary it shows
the effect of the triangle filter, which takes off energy from scales, which
are larger than the filter width. This is in good agreement with the investigations
of Piomelli et al. (1987) and is
a general characteristic of non-perfect filter functions. The largest represented
frequency of
Hz is explicitly determined by
the spatial resolution.

Fig. 5 Frequency spectra in a lake like geometry
Figure 5 (right) shows the modeled frequency
spectra for a FE node in the bank area. The average speed at this node is equal
to
m/s and causes the slight shift
of the maximum frequency of
to
Hz. It can be noticed that the
spectra is dominated by a
gradient and thus in agreement
with the theoretical and measured spectral characteristics in the bank area.
The simulated results are in good agreement with theoretical distributions (FJORTOFT,
1953 and KRAICHNAN, 1967) as well as measurements in Lake Ontario (PALMER, 1973),
where a
was measured in shallow water (9
m depth) and
in deep water (22 m depth).
The studied idealized lake geometry proved the complexity of the flow phenomena in large-scale geometry. Even though the topography had moderate irregular boundaries high velocity gradients were induced and the simulation results proved the three-dimensionality of the flow. The analysis of the energy transport and spectra proved the large spectrum of turbulence. In conclusion the presented simulation results demonstrate the quality of the numerical model and the LES approach. The model is capable of representing the whole eddy spectrum and simulating three-dimensional flow in the irregular boundaries.
For the future use the application to natural geometry is desired. For the validation and calibration of numerical models measurements have to be carried out. The Acoustic Doppler technique is widely used in natural waters and turbulent phenomena can be evaluated. A measurement campaign is carried out in a German reservoir. The main flow phenomena as well as the evaluated turbulent structures will be compared with simulation results. The presented results prove the validity of PASTIS-3D for the study of flow and turbulence phenomena in lakes and reservoirs. Thus a detailed evaluation of flow and turbulence phenomena in natural waters is possible.
Aldama, A. A. (1993): Leonard and Cross-Term Approximations in the Anisotropically Filtered Equations of Motion, in: Large Eddy Simulation of Complex Engineering and Geophysical Flows, Eds. Galperin, Orszag, Cambridge University Press.
Bergen, OLAF (1999): Die Large-Eddy Simulation von Strömungen in natürlichen Seen und Talsperrenspeichern mit der Finite Elemente Methode, Dissertation, Mitteilungen des Instituts für Wasserbau und Wasserwirtschaft, 119 RWTH, Aachen.
Comte-Bellot (1963): Turbulent flow between two parallel walls. Univ. Grenoble, Thesis, Aero. Res. Council no.31609.
Daniels, H. (1992): PASTIS-3D, Finite-Element Projection Algorithm Solver for Transient Incompressible Flow Simulations, UCRL-MA-111833, Lawrence Livermore National Laboratory, California, USA.
Ekman (1905): On the Influence of the Earth’s Rotation on Ocean Currents. Roy. Swedish Acad. Of Sciences, Stockholm – Upsala.
Fjortoft, R. (1953): On the Changes in the Spectral Distribution of Kinetic Energy for Two-Dimensional Non-Divergent Flow; in: Tellus, Vol. 5, pp. 225-230.
Forkel, C. (1995): Die Grobstruktursimulation turbulenter Strömungs- und Transportprozesse in komplexen Geometrien, Dissertation, Mitteilungen des Instituts für Wasserbau und Wasserwirtschaft, 102, RWTH Aachen.
Gresho, P.M. (1990): On the Theory of Semi-Implicit Projection Methods for Viscous Incompressible Flow and Its Implementation Via a Finite Element Method that also Introduces a Nearly-Consistent Mass Matrix, Part 1: Theory, International Journal for Numerical Methods in Fluids, Vol. 11, J. Wiley & Sons, Inc., pp. 587-620.
Gresho, P.M.; CHAN, S.T. (1990): On the Theory of Semi-Implicit Projection Methods for Viscous Incompressible Flow and Its Implementation Via a Finite Element Method that also Introduces a Nearly-Consistent Mass Matrix, Part 2: Implementation, International Journal for Numerical Methods in Fluids, Vol. 11, J. Wiley & Sons, Inc., pp. 621-660.
Imberger, J.; PATTERSON, J.C. (1990): Physical Limnology, Advances in Applied Mechanics (Ed. T. Wu), 27, pp. 303-475.
Kraichnan, R. H. (1967): Inertial Ranges in Two-Dimensional Turbulence; in: Physics of Fluids, Vol. 10, pp. 1417-1423.
Palmer, M. D. (1973): Some Kinetic Energy Spectra in a Nearshore Region of Lake Ontario; in: Journal of Geophysical Research, Vol. 28, pp. 3585-3595.
Piomelli, U.; Ferziger, J. H.; Moin, P. (1987): Models for Large Eddy Simulations of Turbulent Channel Flows Including Transpiration; Thermosciences Division, Dep. Mech. Eng., Stanf. Univ., CA, Report TF-32.
Rettemeier, K.; van Linn, A.; Köngeter, J. (2000): Numerical Model Investigations to Identify Flow Phenomena and Turbulence in Lakes and Reservoirs, Proceedings of the 4th International Conference on Hydro-Science and -Engineering (ICHE): Seoul, Republic of Korea, September 26-29, 2000 / Ed. by Y.N. Yoon [et al.].
Schweim, C.; Zhou, J.; Prochnow, J.V.; Köngeter, J. (2000): Large-Eddy Simulation of a Lid-Driven Rotating Annular Flume Flow. 4th International Conference on Hydroinformatics, 23-27 July, Iowa.