Simulation of Hydrodynamics and
Transport Processes in Surface
Water Systems in Berlin
vorgelegt von
AYMAN JOURIEH, Ing.
aus Tartous, Syrien
von der Fakultät VI Planen | Bauen | Umwelt der Technischen Universität Berlin zur Erlangung des akademischen Grades
Doktor der Ingenieurwissenschaften Dr.-Ing.
genehmigte Dissertation Promotionsausschuss:
Vorsitzender: Prof. Dr.-Ing. Matthias Barjenbruch, TU Berlin Berichter: Prof. Dr.-Ing. Reinhard Hinkelmann, TU Berlin
Prof. Dr. rer. nat. Gunnar Nützmann, IGB Berlin Tag der wissenschaftlichen Aussprache: 17.12.2013
On the one hand, urban water systems are stressed by numerous loads, for example contaminations from treated wastewater, industry, runoff of streets and on the other hand their functioning is highly desirable to maintain the
urban environment. The impact of the contaminations can be predicted
with numerical models which help to understand and identify the relevant processes and which can be applied to design engineering measurements. This doctoral thesis deals with two- and three-dimensional numerical simu-lations of hydrodynamics and conservative transport processes in two urban surface water systems of Berlin, a section of the river Spree and the Un-terhavel water system, both characterized by slow flow velocities. For this purpose, the TELEMAC modelling system was used and 2D as well as 3D simulations were carried out considering various conditions: low, mean and high discharge, tracer injection points and diffusivity. Finally, the impact of mean and strong winds was considered for the Unterhavel.
The section of the river Spree has a very simple geometry and is stressed by combined sewer overflow which is idealized as conservative tracer. The project area under investigation here and within the research project SPREE-2011 is located in Berlin, Friedrichshain-Kreuzberg. This study aims to show the impacts of installed storage tanks on the hydraulics and water quality. The numerical results show that the impacts of the tanks on the hydraulics are small, and the transport is strongly advection-dominated except for low flow conditions. The results indicate that a 3D model is only necessary in the direct area surrounding the injection point.
The geometry in the Unterhavel water system is highly complex and therefore
it is a challenge for stable 3D simulations. This study aims to improve
understanding 2D and 3D flow and transport processes in the Unterhavel and to determine dominated processes and parameters. The results show that the Wannsee is strongly influenced by the injection at the Teltow channel but only slightly influenced by the injection at Pichelssee. In addition, a tracer injected in the north at Pichelssee mainly flows from North to South through the system. Moreover, the transport is mainly advection-dominated,
large parts of the domain, except in special stagnation areas. Finally, in the case of strong wind, strong 3D flow and transport effects occur with different flow directions in a profile at the surface (following the wind direction) and opposite flow direction at the bottom as well as complex horizontal and vertical circulations.
Auf der einen Seite sind urbane Wassersysteme von zahlreichen Belastun-gen betroffen: gereinigtes Abwasser, Industrie, Zuflüsse aus versiegelten Flä-chen. Auf der anderen Seite ist deren Funktionieren unabdingbar, um un-sere Umwelt nachhaltig zu schützen. Die Auswirkungen der Kontaminatio-nen könKontaminatio-nen mit Hilfe von numerischen Modellen vorhergesagt werden. Diese Modelle tragen zum Verständnis der Gewässersysteme und deren relevanten Prozesse bei und können zur Entscheidung zu ingenieurtechnischen Maßnah-men herangezogen werden.
Die vorliegende Dissertation beschäftigt sich mit zwei- und dreidimensio-nalen Strömungs- und Transportmodellierungen konservativer Stoffe in den beiden Berliner Gewässern Spree und Unterhavel. Bei diesen beiden Ge-wässerabschnitten handelt es sich um sehr langsam fließende Gewässer. Für die 2D und 3D Berechnungen der Strömungs- und Transportprozesse de das Programmsystem TELEMAC eingesetzt. In den Simulationen wur-den verschiewur-dene hydrodynamische sowie Transportbedingungen untersucht: Niedrig-, Mittel- und Hochwasserabfluss, Einleitungsstellen der Kontamina-tionen und Diffusivität. Zusätzlich wurden die Auswirkungen mittlerer und starker Windereignisse in der Unterhavel simuliert.
Der Flussabschnitt der Spree, der eine einfache Geometrie hat, wird durch Überläufe aus der Mischwasserkanalisation belastet. Das gereinigte Abwas-ser wurde in dieAbwas-ser Arbeit als konAbwas-servativer Tracer idealisiert. Das Untersu-chungsgebiet befindet sich in Berlin, Friedrichshain-Kreuzberg. Diese Arbei-ten waren in das interdisziplinäre BMBF-Verbundprojekt SPREE2011 ein-gebunden. Ziel dieser Untersuchungen war, die Auswirkungen von Mischwas-serspeichern auf die Hydrodynamik und Wasserqualität zu analysieren. Die Ergebnisse zeigen, dass der Einfluss der Mischwasserspeicher auf die Hydro-dynamik vernachlässigbar ist. Die Transportprozesse werden von der Advek-tion dominiert, außer bei Niedrigwasserabfluss. Die Ergebnisse zeigen, dass ein 3D-Modell nur für die unmittelbare Umgebung des Einleitungspunkts notwendig ist.
Kom-Strömungs- und die Transportprozessen in der Unterhavel zu verbessern sowie dominante Prozesse und Parameter herauszuarbeiten. Die Ergebnis-se zeigen, dass der WannErgebnis-see stark von Einleitungen in den Teltow-Kanal, aber nur wenig von Einleitungen in den Pichelssee beeinflusst wird. Dar-über hinaus fließt ein Tracer, der im Norden des Pichelssees eingeleitet wird, hauptsächlich von Nord nach Süd. Weiterhin ist der Transport überwie-gend akvektionsdominiert, außer in Ruhewasserzonen, wo Diffusion an Be-deutung gewinnt. Wenn kein Wind vorhanden ist, sind die 3D Profile von Strömungsgeschwindigkeit und Transport in weiten Bereichen des Untersu-chungsgebiets parabolisch, außer in Ruhewasserzonen. Bei Starkwindereig-nissen treten starke dreidimensionale Strömungs- und Transporteffekte auf mit unterschiedlichen Strömungsrichtungen an der Gewässeroberfläche (in Windrichtung) und am Gewässerboden sowie komplexe horizontale und ver-tikale Zirkulationsfelder.
Die vorliegende Arbeit ist im Rahmen eines Promotionsstipendiums im Fach-gebiet Wasserwirtschaft und Hydrosystemmodellierung, Institut für Bauin-genieurwesen der Technischen Universität Berlin entstanden.
Mein besonderer Dank gilt meinem Doktorvater, Herrn Prof. Dr.-Ing. Rein-hard Hinkelmann für die Möglichkeit dieses Dissertationsthema zu bearbei-ten. Er hat die Entstehung dieser Arbeit mit viel Interesse, Vertrauen und Geduld unterstützt. Ohne seine ständige Unterstützung wäre diese Arbeit nicht möglich gewesen.
Weiterhin möchte ich Herrn Prof. Dr. rer. nat. Gunnar Nützmann (IGB Berlin) für die freundliche Übernahme des Zweitgutachtens danken. Seine fachlichen Kenntnisse haben mir wichtige Anregungen geliefert.
Auch möchte ich mich bei Herrn Prof. Dr.-Ing. Matthias Barjenbruch für die Übernahme des Vorsitzes des Promotionskomitees bedanken.
Bei allen aktuellen und ehemaligen Kollegen des Fachgebiets Wasserwirt-schaft und Hydrosystemmodellierung möchte ich mich für ihre Hilfsbereit-schaft in wissenHilfsbereit-schaftlichen und technischen Fragen, sowie für die freundliche Arbeitsatmosphäre bedanken.
Meinen Eltern und Geschwistern, die mir das Studium des Bauingenieur-wesens ermöglicht haben, fühle ich mich zu besonderem Dank verpflichtet. Meine Eltern haben mich durch Rat und Tat und vor allem unzählbare Ge-bete unterstützt.
Schließlich gilt mein ganz besonderer Dank meiner Frau Shirin und meiner Tochter Sarah für ihr unbeschreibliches Verständnis und ihre Geduld insbe-sondere an den zahlreichen Tagen, die ich im Institut verbracht habe. Berlin, Dezember 2013
Nomenclature XI
List of Figures XIII
1 Introduction 1
1.1 Numerical simulation of hydro- and environmental systems . 1
1.2 Why urban surface water systems? . . . 3
1.3 Why Berlin? . . . 6
1.4 Why multi-dimensional modeling? . . . 10
1.5 Study objectives . . . 11
1.6 Organization of this work . . . 11
2 Numerical Modeling of Surface Water Systems 13 2.1 Model choice . . . 13 2.1.1 Model testing . . . 14 2.1.2 Model dimensions . . . 15 2.1.3 Modeling systems . . . 18 2.2 Grid Generation . . . 19 2.2.1 Introduction . . . 19 2.2.2 Delaunay-Voronoi Method (DVM) . . . 21
2.2.3 Advancing Front Method (AFM) . . . 21
2.2.4 Grid quality . . . 22
2.3 Discretization Methods . . . 23
2.3.1 Introduction . . . 23
2.3.2 Time discretisation . . . 24
2.3.4 Initial and boundary conditions . . . 32
3 Modeling System 35 3.1 TELEMAC modeling system . . . 35
3.1.1 Introduction . . . 35
3.1.2 TELEMAC-2D . . . 39
3.1.3 TELEMAC-3D . . . 41
3.1.4 Pre- and postprocessing . . . 45
3.2 Grid generator JANET . . . 46
4 Multi-dimensional Numerical Simulation of Hydrodynamics and Transport Processes in the Spree 47 4.1 Project area and computational domains . . . 47
4.2 Computational domain A . . . 52
4.2.1 2D simulations of hydrodynamics for mean discharge condition . . . 54
4.3 Computational domain B . . . 61
4.3.1 2D simulations of hydrodynamics and transport for mean discharge condition . . . 63
4.3.2 2D simulations of transport for variable discharge con-ditions . . . 68
4.3.3 Investigation of the diffusivity . . . 72
4.3.4 3D simulation of hydrodynamics and transport . . . . 76
5 Multi-dimensional Numerical Simulation of Hydrodynamics and Transport in the Unterhavel 83 5.1 Project area and computational domain . . . 83
5.2 Grid generator JANET and pre-processing stage . . . 86
5.3 2D simulation of hydrodynamics . . . 92
5.3.1 Investigation of mean discharge condition . . . 92
5.3.2 Investigation of high discharge condition . . . 105
5.4 2D simulation of transport processes . . . 109
5.4.1 Investigation of mean discharge condition . . . 109
5.4.3 Investigation of high discharge condition . . . 127
5.5 3D simulation of hydrodynamics and transport processes . . . 132
5.5.1 3D grid . . . 132
5.5.2 Investigation of mean discharge condition . . . 138
5.6 Multi-dimensional simulation of wind-driven flow and trans-port processes . . . 148
5.6.1 Conceptual approach and wind data . . . 148
5.6.2 Results of 2D simulations . . . 152
5.6.3 Results of 3D simulations . . . 159
6 Conclusions and Outlook 171 6.1 Conclusions . . . 171
6.2 Outlook . . . 174
Abbreviations
1D one-dimensional
2D two-dimensional
3D three-dimensional
AFM advanced-front method
a.s.l. above sea level
BC boundary condition
BCGSTAB Biconjugate Gradient Stabilized method
CFD computational fluid dynamics
CFD Courant-friedrichs lewy
CG Conjugate Gradient
DVM Delaunay-Voroni method
fig. figure
FDM Finite Difference Method
FEM Finite Element Method
FVM Finite Volume Method
GMRES Generalised Minimum RESidual
MPI message passing interface
PCG Preconditioned Conjugate Gradient
sec. section
SUPG Streamline Upwind Petrov-Galerkin
Terms with latin letters
−
→u velocity vector [m/s]
u, v velocity components [m/s]
T tracer concentration [g/l]
g gravitational acceleration [m/s2]
zs free surface elevation [m]
x, y, z horizontal and vertical space coordinates [m]
Sx, Sy source or sink terms (wind, Coriolis force,
bottom friction, a source / sink of momentum)
[kg/m3s]
ST source or sink of tracer [kg/m3s]
υt turbulent viscosity [m2/s]
υT turbulent diffusivity [m2/s]
fx, fy source or sink terms, denoting bottom friction,
wind and Coriolis force
[m/s2] λ Taylor.friction coefficient [-] cD wind-stress coefficient [-] ρa density of air [kg/m3] ρ0 reference density [kg/m3] 4ρ density difference [kg/m3] p pressure [Pa]
vax, vay horizontal components of wind velocity [m/s]
υt turbulent viscosity tensor [m2/s]
υT turbulent diffusivity tensor [m2/s]
υh horizontal turbulent viscosity [m2/s]
υv vertical turbulent viscosity [m2/s]
υT h horizontal turbulent diffusivity [m2/s]
υT v vertical turbulent diffusivity [m2/s]
vax, vay components of wind velocity [m/s]
ρair
ρw ratio of air and water densities [-]
1.1 Overview of urban water systems (NÜTZMANN et al., 2013) 5
1.2 Overview of the Germany waterways (top) and Elbe and Oder
waterways (bottom) (Bundeswasserstraßen, 2013) . . . 8
1.3 Overview of the surface-water system, flow directions and
wastewater treatment plants (WWTP) in Berlin (modified)
(SenGUV, 2006) . . . 9
2.1 Model dimensions and types . . . 17
2.2 Explicit and (fully) implicit methods, after (HELMIG, 2000;
HINKELMANN, 2005) . . . 27
2.3 Discretisation methods (after HINKELMANN (2005)) . . . . 32
3.1 TELEMAC system . . . 38
4.1 Project area and the computational domains . . . 48
4.2 Natural and simplified cross section in domain A . . . 50
4.3 Project area and available vertical profiles (values in [m]) . . 51
4.4 Bottom elevation and 2D grids; top: scenario A, middle:
sce-nario B and bottom: scesce-nario C . . . 53
4.5 Computational domain A and boundary conditions . . . 54
4.6 Free surface for three scenarios after 2 hours, top: scenario A,
middle: scenario B and bottom: scenario C . . . 57
4.7 Scalar velocity for three scenarios after 2 hours, top: scenario
A, middle: scenario B and bottom: scenario C . . . 58
4.8 Temporal variation of scalar velocity after 2 hours at points
4.9 Spatial variation of scalar velocity at sections X-X and Y-Y,
after 2 hours, for three scenarios A, B and C . . . 60
4.10 Computational domain B and boundary conditions . . . 61
4.11 Bottom elevation and 2D grid . . . 62
4.12 Flow field for mean discharge condition and tracer
concentra-tion after 600 seconds . . . 65
4.13 2D tracer concentration after three time steps for mean
dis-charge condition . . . 66
4.14 Temporal variation of tracer concentrations at points X, Y
and Z within 10 hours for mean discharge condition . . . 67
4.15 2D tracer concentration after 5 hours for low, mean and high
discharge conditions . . . 70
4.16 Temporal variation of tracer concentration at point X for low,
mean and high discharge conditions . . . 71
4.17 2D tracer transport after three time steps, left: νT = 0.01 m2/s,
right: νT = 0.0001 m2/s . . . . 74
4.18 Left: spatial variation of tracer concentration at section B-B
[left middle: νT = 0.01 m2/s, left bottom: νT = 0.0001 m2/s].
Right: temporal variation of tracer concentration at point A:
[right middle: νT = 0.01 m2/s, right bottom: νT = 0.0001 m2/s]
. . . 75
4.19 Top: 3D grid at the section B-B, bottom: the considered
sec-tions . . . 79
4.20 Flow field at section Y-Y for scenario A . . . 80
4.21 Distribution of tracer concentration in vertical section Y-Y for
three time steps . . . 81
4.22 Flow field and distribution of tracer concentration in vertical section Z-Z (shown in fig. 4.19) for scenario B for four time
steps . . . 82
5.1 Study area, computational domain and model concepts . . . 85
5.2 DTM and preprocessing stage . . . 87
5.4 Bottom elevation and details of the 2D grid . . . 91
5.5 Computational domain with boundary conditions for the mean
discharge conditions . . . 93
5.6 Comparison of scalar velocity for three friction coefficients at
a section A-A (fig. 5.11) . . . 95
5.7 Spatial variation of water depth for the mean discharge
condi-tion after 3 days; left Jungfernsee, middle whole domain, right
top Schwanenwerder island, right bottom Wannsee . . . 98
5.8 Temporal variation of water depth and free surface at
consid-ered points X, Y and Z for mean discharge condition . . . . 99
5.9 Velocity field for the mean discharge condition in the whole
domain and several subdomains . . . 100
5.10 Temporal variation of scalar velocity (left) and water depth (right) at the considered points X, Y and Z in the first 3 days
for mean discharge condition . . . 101
5.11 Spatial variation of water depth and scalar velocity at different
sections after 3 days for mean discharge condition . . . 102
5.12 Spatial variation of flowrate (discharge per meter width) at different sections in the first 3 days for mean discharge
condi-tion . . . 103
5.13 Correlations between scalar velocity and water depth for mean
discharge condition . . . 104
5.14 Computational domain and boundary conditions for high
dis-charge condition . . . 106
5.15 Flow field for high discharge condition . . . 107 5.16 Temporal variation of scalar velocity (left) and water depth
(right) at three considered points X, Y and Z in the first three days for mean and high discharge condition . . . 108 5.17 Computational domain and boundary conditions for scenarios
A, B and C . . . 112 5.18 2D tracer transport after 4 days for mean discharge condition
5.19 2D tracer transport after 30 days for mean discharge condition
for scenarios A, B and C . . . 114
5.20 2D tracer transport after 120 days for mean discharge condi-tion for scenarios A, B and C . . . 115 5.21 2D tracer transport at Wannsee after 4 days for mean
dis-charge condition for scenarios A, B and C . . . 116 5.22 2D tracer transport at Wannsee after 30 days for mean
dis-charge condition for scenarios A, B and C . . . 117 5.23 2D tracer transport at Wannsee after 120 days for mean
dis-charge condition for scenarios A, B and C . . . 118 5.24 Temporal variation of tracer concentration at three
consid-ered points for three scenarios A, B and C for mean discharge
condition . . . 119
5.25 Path lines of two particles after 120 days . . . 120
5.26 2D tracer transport after 5 days; left: νT = 0.0001 m2/s,
mid-dle: νT = 0.001 m2/s, right: νT = 0.01 m2/s . . . 123
5.27 2D tracer transport after 10 days; left: νT = 0.0001 m2/s,
middle: νT = 0.001 m2/s, right: νT = 0.01 m2/s . . . 124
5.28 2D tracer transport after 20 days; left: νT = 0.0001 m2/s,
middle: νT = 0.001 m2/s, right: νT = 0.01 m2/s . . . 125
5.29 Temporal variation of tracer concentration at three considered
points for the cases: νT = 0.0001 m2/s, νT = 0.001 m2/s and
νT = 0.01 m2/s . . . 126
5.30 2D tracer transport after 2 days (top) and 10 days (bottom)
for mean and high discharge conditions . . . 129
5.31 2D tracer transport after 30 days (top) and 120 days (bottom)
for mean and high discharge conditions . . . 130
5.32 Temporal variation of tracer concentrations at three
consid-ered points for mean and high discharge conditions . . . 131
5.33 Vertical discretization using A): z-method and B): σv-method; elements (left) and nodes (right) . . . 135
5.34 Vertical discretization using (A): z-method; (B): σv-method: planes depend on free surface and bottom, and (C) σv-method:
planes depend only on bottom . . . 136
5.35 Vertical distribution of tracer concentration using the σv-method
(top) and the z-method (bottom) . . . 137
5.36 CPU time for parallel computing . . . 142 5.37 Vertical velocity profile at section H-H simulated with
TELEMAC-3D (middle); theoretical profile for vertically averaged flow (right) . . . 144 5.38 3D flow field at three horizontal planes: A) free surface; B)
mid-depth and C) bottom . . . 145 5.39 3D flow field and vertical distribution of tracer concentrations
at different time steps at section A-A . . . 146
5.40 Temporal variation of velocity and tracer concentrations for
different time steps and at three different water depths . . . 147
5.41 Wind velocities and directions at Wannsee: period 2001-2009 (top), period January-December 2007 (middle), January 17-19, 2007 (bottom, with mean condition (A) and extreme
con-dition (B)) . . . 151
5.42 Wind field (left), 2D spatial variation of free surface (top) and flow field (bottom) at Wannsee with south wind (middle) and without wind (right), after 5 hours . . . 153 5.43 Wind field (left), 2D flow field at Wannsee with south wind
(top) and and without wind (bottom), after 5 hours . . . 154
5.44 Wind field (left), 2D tracer concentration and flow field at Wannsee with south wind (top) and without wind (bottom),
after 5 hours . . . 155
5.45 Wind field (left), 2D spatial variation of free surface (middle top: with south wind, middle bottom: without wind, and right top), scalar velocity (right middle) and tracer concentration (right bottom) at Wannsee , with and without southern wind
5.46 Path lines of two particles after 30 days, top: with south wind,
bottom: with west wind. . . 158
5.47 Theoretical wind-induced velocity profile (after LUI and PEREZ (1971); SHANAHAN et al. (1981)) . . . 159 5.48 Velocity field: (A: horizontal circulation at section T-T, B:
vertical section at section P-P), induced by west wind after 5
hours . . . 160
5.49 Left: Vertical distribution of tracer and velocity (top: without wind, bottom: with west wind), right: Spatial distribution of tracer concentration and velocity at 0.5 m under free surface, at section C-C (top: without wind, bottom: with west wind)
after 5 hours west wind . . . 163
5.50 Vertical distribution of tracer and velocity (m/s) for different time steps at the section S-S (see fig. 5.48); left: without wind, right: with west wind, after 5 hours west wind . . . 164 5.51 Left: Vertical distribution of tracer and velocity (top: without
wind, bottom: with west wind), right: Spatial distribution of tracer concentration and velocity at 0.5 m under free surface at section A-A [shown in figure 5.51, middle]: (top: without
wind, bottom: with west wind) after 5 hour west wind . . . 165
5.52 Left: Vertical distribution of tracer and velocity (top: without wind, bottom: with west wind), right: Spatial distribution of tracer concentration and velocity at 0.5 m under free surface, at section B-B (top: without wind, bottom: with west wind)
after 5 hour west wind . . . 166
5.53 Left: Vertical distribution of tracer and velocity (top: without wind, bottom: with west wind, right: Spatial distribution of tracer concentration and velocity at 0.5 m under free surface at section C-C (top: without wind, bottom: with west wind)
after 5 hour west wind . . . 167
5.54 West wind condition and flow field at three horizontal levels at Wannsee: A) free surface, B) mid-depth, C) bottom, after 5 hour . . . 168
5.55 South wind condition and the flow field at three horizontal levels at Wannsee: A) free surface, B) mid-depth, C) bottom,
This chapter gives a short overview on numerical modeling of hydro- and environmental systems, including the explanation of the importance of urban water systems. The surface water systems in Berlin are briefly introduced. Finally, the need for multi-dimensional modeling, objectives and the organi-zation of this work are also addressed.
1.1 Numerical simulation of hydro- and environmental
systems
“Modeling is a little like art in the words of Pablo Picasso. It is never com-pletely realistic; it is never the truth. But it contains enough of the truth, hopefully, and enough realism to gain understanding about environment sys-tems. Schooner, 1996” (ZHEN-GANG, 2008).
There has been an increasing development of numerical models as they be-come more prevalent and have found their way to simulate the hydrodynam-ics and water quality of water systems.
Every surface water system has its special conditions yet many systems have similar hydraulic and environmental problems: flood and droughts, con-tamination, sedimentation, eutrophication. These problems, in turn, can adversely affect water quantity- quality and maintenance of environmental systems.
The Navier-Stockes equations present the fluid dynamics which is used to describe the transport of sediments, salinity, other constituents and interac-tion processes between constituents (FERZIGER and PERICA, 1996). In
this work the shallow water equations which are special case of the Navier-Stockes equations are chosen as governing equations. For a simple geometry these equations can be solved analytically to describe a system. However, numerical methods are necessary for complex geometries of natural water systems with varying boundary conditions such as tides, floods and winds.
Generally, measured data alone is rarely sufficient to make decisions on water quantity and quality management plans, especially when the water systems are large, complex and boundary conditions are changeable. Due to time, cost and technical constraints, field measurements are often limited to certain small areas and for certain periods of time.
Models enable the selection of a better, more physically based choice among alternatives and the evaluation of the more effective method in solving the water quantity or quality problem. Models also serve to identify data gaps in characterizing water systems. Furthermore, they can be used to analyze the impact of different management alternatives and to select the one that results in the least adverse impact on the environment. If there is a lack in understanding a complex water system, a numerical model can help to get insight, and therefore assist in the estimation of important hydraulic and environmental parameters which in turn enables the better estimation of suitable measures and decisions (ZIELKE et al., 1999).
Consequently, numerical models, supported by numerous pre- and post-processing tools, help to gain a better understanding of the hydrodynam-ical behaviour of water systems, to provide information for the engineering studies, to compute/compare different scenarios and then to assess the best variant, which helps decision-making (based on simulation) and enables to respect the environmental systems. Finally, it should be kept in mind, that a good numerical simulation of the hydrodynamics plays a vital role in a successful modeling application, especially if water quality problems are in-vestigated.
1.2 Why urban surface water systems?
Although urban water contributes only a very small part to the whole water cycle, it is very important for maintaining our environment. These water systems do not only consist of rivers, creeks, canals, (artificial) lakes and groundwater that may serve the water supply in an urban area, but also (combined) the sewer systems, water treatment plants and water distribu-tion systems that transport the water to where the demands are located. Urban water systems are generally contaminated by treated waste water discharges, private households, industry, shipping, recreation activities and run-off from streets. Furthermore, these different water systems are in more or less strong interaction with combined sewer overflows in rivers or bank filtration from lakes (see fig.1.1). Moreover because of the continuous move-ment of populations to urban areas for a better standard of living, these urban water systems are under increasing threat of contamination (SINGH,
1995; KOPMANN and MARKOFSKY, 2000). The impact of the above
mentioned conditions can be predicted by use of mathematical models. From both an hydrodynamic and environmental point of view, modeling of water- waste water flows, water quality, water treatment and water distri-bution is a challenging task, not only because of its hydraulic complexity, but also because of stochastic inputs, changing boundary conditions and demands on the systems, for example: climate and demographic changes (e.g. floods, droughts), variations in water use, new pollutants (e.g. drugs, micro-pollutants), new water technologies (e.g. water treatment), changes in managing urban water (e.g. more decentral infiltration, more storage and use). In addition, the ecological condition of a waterway is an example of a complex system of related variables in which a change in any condition may cause a reaction in many system processes (SCHEUERMANN, 2008). The ecosystem of urban waters is an interacting system that includes hydro-dynamic (e.g. water depth and flow velocity), chemical (e.g. solids, dissolved oxygen, nutrients), and biological (e.g. fish, benthos, plants) characteristics. These processes (physical, chemical, and biological) also vary both in time
(temporal variations may have long-term (years) and short-term (minutes) time scales) and in space (spatial variations largely depend, for example on the complexity of topography, structure and external discharges) (KOP-MANN and MARKOFSKY, 2000).
The improvements and techniques for the management and protection of surface and groundwater resources must be based on the European Water Framework Directive (EWFD) which has the overarching goal of achieving a good ecological state for all water systems. As many urban water systems have been strongly altered in the past and are strongly degraded, the classifi-cation heavily modified water body (HMWB) has been introduced for which a good ecological potential should be achieved instead of a good ecological state, being less restrictive. Nevertheless reaching the aims of the EWFD is a challenging task for many cities.
Taking the aforementioned conditions into consideration, it can be concluded that the simulation of water quantities and qualities in urban water systems provides necessary information for the predictive measures to improve water quantity and quality, managing and sustainable using urban water with a long-term perspective to protect the water resources and the environment.
1.3 Why Berlin?
Berlin is located between the two water systems of the Elbe and Oder rivers (see fig. 1.2). The centre of the city is located on the Spree, which runs from east to west across the city. The river Dahme joins the Spree at Köpenick, in the eastern area of the city Berlin. At Spandau the Spree joins the Unter-havel, which flows north to south along the city’s western boundary (see fig. 1.3). Further important creeks are Panke, Wuhle and Erpe (STEEG, 2006; SenGUV, 2006; ZIEGLER, 2001; WIESE, 2007; MÖLLER, 2008;
LINDEN-SCHMIDT and FRÖHLICH, 2000).
Berlin waters are strongly backwater regulated due to many channels. The most important channels in Berlin are the Teltow channel, the Landwehr channel and the Spandau shipping channel mostly running from east to west between the Spree and Unterhavel rivers, which are the most important rivers in Berlin. The Teltow channel runs from the Dahme south of Köpenick through the southern suburbs of Berlin to an arm of the river Unterhavel east of Potsdam whereas the Spandau shipping channel runs from the Spree near the main station to the river Havel north of Spandau. The Landwehr channel - a shorter channel – runs from the Spree between Treptow and Kreuzberg and joins the Spree again in Charlottenburg (STEEG, 2006; SenGUV, 2006; ZIEGLER, 2001; WIESE, 2007; LINDENSCHMIDT and FRÖHLICH, 2000;
MÖLLER, 2008).
Berlin is characterised by moderate precipitation (∼ 600 mm/year). Berlin waters are denoted as “Gewässerreich, aber wasserarm”, meaning that there are many water bodies, but having overall low discharges.
In Berlin, the lakes are shallow and the rivers are slow-running, carrying little water in summer, yet, large amounts of polluted storm water flowing through the drainage systems may overflow and reach the water systems when it rains. This makes Berlin’s water bodies ecologically very sensitive systems. The biggest lakes are Tegeler See and Müggelsee. Further, more than 90% of the drinking water for Berlin is obtained from groundwater. In the suburbs of the city, two-thirds of the water extracted actually consist of bank filtrate or
artificially recharged surface water. Additionally, it is estimated that about 40% of groundwater reserves are threatened by contamination, mainly by some old polluted sites and areas not yet connected to the sewage system (STEEG, 2006; SenGUV, 2006; ZIEGLER, 2001; LINDENSCHMIDT and FRÖHLICH, 2000; WIESE, 2007; MÖLLER, 2008).
Finally, it should be kept in mind, that the quality of both the surface and groundwater is very important for the drinking water supply. Further, the main water quality problems are eutrophication, the associated growth of algae, and damage to the oxygen balance caused by eight major sewage plants as well as heavy metals, phenol, hydrocarbons, and ammonium. Moreover, improvements have been made in the last years in Berlin reducing phosphate and nitrogen spreading from wastewater treatment plants in order to obtain the clean-up goal of water quality class II (slightly eutrophic). However, the aim of EWFD is hardly to achieve for many urban systems (STEEG, 2006; SenGUV, 2006; WIESE, 2007; MÖLLER, 2008).
After this short introduction to Berlin´s water systems, it can be pointed out that a sustainable water resource management is urgently necessary for both the quantity and quality of the water to fit the demands of today and the future.
Berlin
Germany waterways
Figure 1.2: Overview of the Germany waterways (top) and Elbe and Oder waterways (bottom) (Bundeswasserstraßen, 2013)
WWTP flow
dire
cti
o
n
input WWTP discharge Pressure
pipe com putational d om ain Figure 1.3: Ov erview of the surface-w ater system, flo w directi on s and w astew ater treatmen t p lan ts (WWTP) in Berlin (mo dified) (S e n GUV, 2006)
1.4 Why multi-dimensional modeling?
Multi-dimensional models are: two-dimensional (vertically averaged here), three-dimensional based on hydrostatic or non-hydrostatic pressure assump-tion or fully three-dimensional. These type of models have been developed to compute: hydrodynamics, transport of mass, heat and water quality (ZIELKE, 2000).
The overarching goals of multi-dimensional models are i) better understand-ing of physical, chemical and biological processes in rivers, estuaries, lakes, and coastal waters, ii) to define intrinsic relationships, for example between
pollutant loadings and water quality. Furthermore, application of
multi-dimensional models detects the gaps which induce their further development to represent the water systems more and more realistically.
The available hydrodynamic and water quality models range from simpli-fied one-dimensional, steady-state models to complex instationary
three-dimensional models of hydrodynamics, transport and water quality. On
the other hand along the development of realistic three-dimensional mod-els computational requirements have changed from supercomputers, to high-end workstations, and nowadays to standard PCs (HINKELMANN, 2005; ZIELKE et al., 1999).
The primary advantage of using multi-dimensional numerical models is their capability of giving details about the flow and transport processes in two or three directions. It is straightforward that a three-dimensional modeling re-quires much more CPU time compared to a two-dimensional one. Therefore, a three-dimensional model should only be chosen if required, as in the case of density and wind-induced flow. In cases where only the horizontal variability of the flow is important, for example flood modeling, a two-dimensional ap-proach is generally sufficient. Furthermore, the modeling costs are generally small compared to the cost of implementing a measure.
Well-known multi-dimensional modeling systems are, for example: TELEM-AC (HERVOUET, 2007; EDF-BRD, 2012), Deltares (DELFT-HYDRAULICS, 1992) and Mike family of DHI (DHI, 1993). These modeling systems include
numerical simulators as well as tools for data analysis, grid generation and postprocessing.
1.5 Study objectives
The main objective of this thesis is to study two- and three-dimensional hy-drodynamic and transport processes in two surface water systems in Berlin: 1. a section of the Spree river being influenced by combined sewer overflow
and
2. a part of Unterhavel river in Berlin being influenced by treated waste water.
These two water systems are investigated for various conditions, such as low, mean, high discharge conditions and for different wind intensities as well as for distinct contaminants. The main purpose of this work consists of understanding the processes, identifying important and less important effects, considering the linkage between flow and transport, showing impacts of technical water systems to natural ones and providing a basis for the estimation of the impacts of possible future measures.
1.6 Organization of this work
This work is divided into five chapters:
Chapter 1 (i.e this chapter) presents an introduction and the motivation
of this work. The need for numerical simulation of hydro- and
environ-mental systems is explained and the importance of urban water systems is emphasized. The surface water system in Berlin is briefly introduced. Fur-thermore, multi-dimensional modeling as well as the objectives of this work are addressed here.
Chapter 2 starts with the method to choose a model through which model testing, model dimensions and modeling systems are addressed. Aspects of
grid generation are shortly dealt with and discretisation methods for time and space (FDM, FEM and FVM) are briefly introduced together with their advantages and disadvantages.
Chapter 3 briefly describes the chosen computational tools. The TELEMAC modeling system, i.e TELEMAC-2D and TELEMAC-3D and its pre- and postprocessor, and further computational tools are briefly introduced. Chapter 4 presents a natural urban water system having a simple geometry which is stressed by combined sewer overflow. The computational domain is a section of river Spree in Berlin. The project SPREE-2011 is briefly intro-duced. This chapter explains the 2D and 3D model of hydrodynamics and transport in detail starting from the grid generation step, setting boundary conditions, carrying out the simulations and visualizing the results. Various conditions are also investigated: low, mean and high discharge, contamina-tions.
Chapter 5 presents a natural urban water system with complex geometry, which is stressed by treated waste water. The computational domain is a part of the Unterhavel water system in Berlin. This chapter explains the 2D and 3D model of hydrodynamics and transport in detail starting from the grid generation step, setting boundary conditions, carrying out the simulations and visualizing the results. Various conditions are investigated: low, mean, and high discharge as well as contaminations. Finally, the effect of wind is also investigated.
In Chapter 6, the content of the work is summarized, conclusions are drawn and an outlook on future research is given.
Water Systems
This chapter starts with choosing a model whereby model testing, model di-mensions and modeling systems are addressed. Aspects of grid generation are shortly dealt with. Finally, discretisation methods for time and space (FDM, FEM and FVM) are briefly introduced together with their advantages and disadvantages.
2.1 Model choice
Numerical models should produce results, which provide information to un-derstand the behavior of the considered water systems. The quantification of the model goodness and the model error is very important. Model goodness means the accuracy of the results, furthermore, the effort which is necessary to set-up the model (ZIELKE, 2000).
The selection of an existing model to be used in any project depends in part on the processes that will be modeled, the data available and the data required by the model.
Finally, it is recommended to set up models that: 1- are easy to understand
2- are compatible with available data
3- provide information about input data, processing, visualizing and analysing the results
4- are easily calibrated and validated (if possible) 5- provide the information needed (ZIELKE, 2000).
2.1.1 Model testing
The model testing is carried out in three steps: verification, calibration and validation:
Verification: The verification is done by comparisons with analytical
solu-tions, which are only available for simple systems and single processes, by plausibility tests, for example, checking the global mass conservation, and by ensuring that programming errors have been removed. In other words, the term verification means proving that the numerical results are correct (HINKELMANN, 2005).
Calibration: During the calibration, the numerical results are compared
with experimental or available field data. Calibration is done for a certain time period. Calibration parameters are the friction coefficient (very im-portant) and the turbulent viscosity (generally minor imim-portant) in a free surface flow simulation and the turbulent diffusivity in a transport simulation in free surface flow. During the calibration these parameters are adapted in physically rational range to produce a sufficient agreement of the simulated output to the field measurements. The calibration is often performed by hand on the basis of engineering experience. The calibration offers a first indication of the model´s quality (HINKELMANN, 2005; ZIELKE, 2000).
Validation: In the validation step, the model is then used to simulate an
independent period of time for which experimental or field data under dif-ferent conditions are available. The same parameters as in the calibration step are chosen and the results of the validation case are compared with the available data. If these results are compatible, it can be concluded that the model is valid. If not, the parameters are modified, and the calibration and
validation steps are repeated. This is done iteratively until the results and field observations are close enough. Generally, it is desirable to validate a model with field data. (HINKELMANN, 2005; FERZIGER and PERICA, 1996; ZIELKE, 2000).
2.1.2 Model dimensions
Based on the number of dimensions, in which hydrodynamic and transport processes occur, the numerical models can be classified as follows:
One-dimensional model: One-dimensional models are the simplest ones
(fig. 2.1). An example of this type is the computation of flow processes in a simple channel. The channel is divided by stream cross sections, and for each time step and the model computes the water depth and cross-sectional averaged velocity at each cross section.
The two main advantages of these models are easy implementation and fast computation whereas the main disadvantage lies in the fact that they can not describe details of horizontal and vertical velocity distributions, such as strat-ification or circulations (ZIELKE et al., 1999; KOPMANN and MARKOF-SKY, 2000).
Two-dimensional model: Two-dimensional models go a step further by
con-sidering the flow directions of the water systems in both horizontal directions. 2D models offer a better comprehension of the physical processes, like advec-tion and diffusion etc. Triangular or quadrilateral elements can be chosen. The triangles are more flexible to adapt to complex boundaries and inner structures (fig. 2.1).
In comparison with one-dimesional models these models require more field data, higher capacity of computation and they are more time consuming to set up and run. These models are not able to compute vertical circula-tions resulting from wind influence or stratification (FRANKE et al., 1987; LYNCH, 1986).
2.5 dimensional model: For 3D-simulations, the mesh is often unstruc-tured in the horizontal plane, and mapping it several times over the vertical direction, it is structured in the vertical plane. This can be called 2.5 dimen-sional, whereby nearly all three-dimensional shallow water models are 2.5 D, however often the term 2.5 D is not used (ZIELKE et al., 1999; WANG, 1988)
Three-dimensional model: In the three-dimensional models, as in 2.5
di-mensional models, not only the horizontal dimensions are modeled, but also
the vertical one. The vertical dimension is generally modeled through a
number of layers. These vertical layers can be equally distributed along the vertical direction or a refinement close to the bottom and the water surface is possible. The latter is generally recommended to better approximate the processes close to the vertical boundaries as the number of layers can be in the range of single digit to small double digit.
Pentahedra or hexahedra can be chosen as elements, while the pentahedras are more flexible to adapt to complex boundaries and complex inner struc-tures (fig. 2.1). Three-dimensional models can be used to predict water level and 3D velocity fields, as they occur in cases of stratification and vertical circulation, for example caused by wind (KOPMANN and MARKOFSKY, 2000; PUTZAR et al., 2010).
Generally, arbitrary three-dimensional elements such as tetrahedron are not used for the three-dimensional shallow water flow.
line traingle quadriteral hexahedron tetrahedron pentahedron
3D model:
2D model:
1D model:
2.1.3 Modeling systems
Following the developments in the computer technology – from microcomput-ers and workstations to supercomputmicrocomput-ers – in the last decades, the extension of existing as well as the development of new models has taken place. We can classify these modeling systems as follows:
• the processes considered: flow, transport of components, sediments, heat and further processes (i.e., water quality, waves, etc).
• the dimensions considered: 1D, 2D, 3D.
• the discretisation methods chosen: FDM, FEM, FVM (see 2.3.3) or others.
A large number of modeling systems are available to simulate the flow and transport processes in surface water. Some are developed for research ob-jectives, whereby others are developed for commercial purposes. Some well-known modeling systems are listed below:
• TELEMAC: It has been developed by National Hydraulics and En-vironmental Laboratory (Laboratoire national d’hydraulique et envi-ronnement de Electricité de France - LNHE-EDF) (EDF-BRD, 2002). This code can model two- and three-dimensional shallow water flow, mass, sediment transport and heat as well as groundwater flow, waves and water quality. It is mainly based on the finite-element method (FEM). Since 2010, it is freely available. As there is experience at the chair of water Chair of Water Resources Management and Model-ing of Hydrosystems TU Berlin in the development and application
(e.g. (HINKELMANN and ZIELKE, 2000; HINKELMANN, 2005;
JOURIEH and HINKELMANN, 2012; OMAR, 2009)), this model has been chosen for this work, a more detailed description is given in section 3.
• MIKE family: These models have been developed and are distributed by DHI (Danish Hydraulic Institute). They model 1D, 2D and 3D shallow water flow, mass, sediment and heat transport as well water
quality and waves and groundwater flow. The models are mainly based on Finite-Difference method (FDM). The code is commercial (DHI, 1993).
• Delft family: These models have been developed and are distributed by Deltares (former Delft Hydraulics), to perform simulations of hy-draulics, mass, sediment and heat transport as well as water quality, waves, groundwater in rivers, lakes, estuaries, bays, coastal areas, seas and other water systems. This code is a commercial code (DELFT-HYDRAULICS, 1992).
• There are further well-known one- and multidimensional models for shallow water flow and transport, for example: HEC-RAS (BEAVERS, 1994), HMS (SIMONS et al., 2012), Hydro-As-2D (Hydrotec, 2012), UnTRIM (CASULLI and LANG, 2002), QUAL2E (RECKHOW, 2003). Generally, shallow water models use the hydrostatic pressure assumption. For selected cases, for example around hydraulic structures (pillars, groynes) a non-hydrostatic extension is necessary, which is included in some of the
models, e.g. TELEMAC. Newest models for the near-field of hydraulic
structures are two-phase (water-gas) Navier-Stokes solvers, e.g. OpenFOAM (OpenFOAM, 2012).
2.2 Grid Generation
2.2.1 Introduction
The first step in numerical modeling is to discretize the computational do-main in one, two or three dimensions, representing the system through a one-, two- or three-dimensional grid. The grid dimensions and the represen-tation of the model quantities affect the quantity of data and time required to set it up.
Based on the following characteristics, the grids can be classified in several categories:
Shape: In 2D grids can be triangular or quadrilateral and in 3D according to the shape of the elements, for example: pentahedral or hexahedral (see fig. 2.1)
Orthogonality: If the angle between crossing grid lines is always 90 degree,
the grid is orthogonal, otherwise non-orthogonal.
Dynamics: Grids can be dynamic, for example in 3D when it adapts to the
moving free water surface or to moving bottom levels in case of erosion and deposition.
Structure: The grids can be structured or unstructured.
Structured grids generally consist of rectangular or quadratic elements in 2D
and bricks or cubes in 3D. In the code direct addressing is possible which
is very fast. This type of grid is not suitable for complex geometry i.e.
complex inner structure and boundaries. However, they require less memory
compared to unstructured grids. The finite difference method (FDM) is
based on structured grids.
For complex- and natural hydro- and environmental systems (as the Unter-havel river in this study, see fig. 5.4), unstructured grids should be chosen. These grids are capable to adapt to complex geometries i.e. complex inner structure and boundaries. In the code addressing is indirect and a little more memory and run time is required compared to structured grids. The FEM and FVM are based on unstructured grids. However also applicable on structured grids (HINKELMANN, 2005).
In this work only unstructured grids have been used.
In the following the two most common methods for generating unstructured grids are briefly introduced: Delaunay-Voronoi Method (DVM) and Advanc-ing Front Method (AFM). They also have been applied in this work.
2.2.2 Delaunay-Voronoi Method (DVM)
If a set of points is considered, a Delaunay triangulation for this set is a triangulation such that there is no point in the plane inside the circumcircle of any triangle. The Delaunay triangulation was developed by Boris Delaunay in 1934 (DELAUNAY, 1934; TANIGUCHI et al., 1992). A preparatory grid is generated based on the available boundary nodes as well as possible further nodes in the domain. The grid must be further processed by, for example the Advancing Front Method (AFM) or other techniques (see sec. 3.2) The properties of Delaunay triangulations can be summarized as follows (GEORGE and BOROUCHAKI, 1998; SHEWCHUK, 1996):
• The minimum angle of all triangles is maximized. • Delaunay triangulations avoid skinny triangles.
• The circumcircle of a triangle contains no other input points. • Every edge is shared by at most two triangles.
• The triangulation must satisfy the on-circle criterion that no point of the considered set points is interior to the circumcircle of any triangle.
2.2.3 Advancing Front Method (AFM)
The basic idea of the Advancing Front Method is dividing the boundaries of the grid into edges in 2D and triangular faces in 3D. The triangles or tetrahedra are generated one-by-one, starting from boundary edges or faces and proceeding forward towards the center of the domain (MÖLLER and HANSBO, 1995; SELLERHOF, 2002; ZIELKE, 2000).
The grid is constructed by progressively adding elements starting at the boundaries. This results in a propagation of a front which is the border (internal boundary) between the meshed and the unmeshed region. The difficulty of this method lies in the merging of the advancing front. A new triangular/tetrahedral element is added by inserting one new point. The
location of this point is the important step and is determined by the following criteria (LO, 1985; HO-LE, 1988; MÖLLER and HANSBO, 1995):
• Quality of the resulting element
• Desired mesh spacing given by the control space
• Neighbourhood constraints like other parts of the boundary / front • The point must be inside the domain.
2.2.4 Grid quality
The quality of results is highly depending on the quality of the grid. There-fore, the quality of the grid has a great influence on: solution accuracy, solver convergence and CPU time required.
The goodness of the grid quality is crucial for a successful simulation, e.g, results of a simulation based on a too coarse mesh cannot be trusted. The ratio between the smallest and largest element edge should not differ by too many orders of magnitude. In surface water systems, the length scales on the horizontal planes (kilometers) are generally much larger than in the vertical direction (meters/ tens of meters). Moreover, it is desirable that
elements are equilateral and angles close to 60o for triangles and 90o for
quadrilaterals (SCHINDLER, 2001; RUPPERT, 1994; SELLERHOF, 2002). Finally, the property common to all classes of numerical methods is that the local accuracy of the approximation is controlled by the grid resolution. Consequently, grid refinement is usually the most obvious way to improve the accuracy of a numerical model. A mesh convergent solution is defined as a solution where the results do not change when the grid is further refined in space and time (HINKELMANN, 2005; MITCHELL and VAVASIS, 1992).
2.3 Discretization Methods
2.3.1 Introduction
Discretization methods transform the partial differential equations to alge-braic equations and approximate the solution function at discrete points of the domain. The properties of these discretization methods are as follows: * Consistency:
A discretization method must be consistent, i.e. the discretization error δ tends to zero for vanishing temporal and spatial discretization sizes:
δ = lim
∆t,∆x→0= |W −W | = 0 (2.1)
Here W is the partial differential equation andW is the discretized equation.
* Stability:
Discretization methods must be stable, i.e. numerical solutions should not diverge. Stabilization techniques add artificial diffusion to achieve the sta-bility. However, they generally lead to numerical diffusion/dispersion. * Convergency:
Numerical solutions should be as close as possible to the exact solutions. This means, that the total error E should vanish for space and time steps tending to zero:
E = lim
∆t,∆x→0|g − g| = 0 (2.2)
Here, g andg denote the solution of the partial differential equation and the
solution of discretized equation, respectively.
Stability and consistency imply convergence to the correct solution. In other words, if a method is stable and consistent, convergence is mathematically proven.
* Monotonicity:
Monotonicity means that no over- or undershooting occurs, i.e. a nodal value should lie between the values of its neighbours. Low order schemes can more easily to be monotonous, but also less precise and more diffusive. Higher order schemes can more precise and less diffusive, but more easily tend to oscillations.
* Conservation:
Numerical methods should ensure the conservation of physical quantities such as mass, momentum, and energy. Some methods are locally conserva-tive, others (only) globally or even not necessary conservative.
For more information about the properties of the discretization methods, see (HINKELMANN, 2005).
2.3.2 Time discretisation
The time domain is subdivided into a number of constant or variable time
steps 4tn where the solution functions are determined. For the time to0, the
initial conditions for all primary variables must be given:
tn= to+ n4tn (2.3)
Here, tn denotes the time after the nth time step and n is the number of
time steps (HINKELMANN, 2005). Often a constant time step is chosen. The time-dependent equations can be written in generalized form in the following way:
∂e
∂t = Ae (2.4)
In this equation, e stands for a scalar or vector entity, and the operator A contains the spatial derivatives as well as the other terms (HINKELMANN,
2005). The first derivative of e is descretised using forward differencing FD (see sec. 2.3.3, eq. 2.12).
The numerical schemes can be classified in two primary categories depending on the approach used to discretise the equations in time:
The first one is explicit, also called Forward Euler Method, where on the right-hand side only values from the previous time step are considered:.
en+1− en
4tn = Ae
n (2.5)
Here, we can compute the unknowns on the new time level n+1 by consid-ering only the known values from the previous time level n. The order of consistency is only O (4t). It is important to mention that explicit meth-ods generally have restriction concerning the time step for stability reasons. Furthermore, the computational effort for one time step is low.
The fully implicit or Backward Euler Method is obtained, if we determine e on the right-hand side on the new time level n+1,
en+1− en
4tn = Ae
n+1 (2.6)
In this case, the unknowns on the new time level n+1 depend on each other (fig. 2.2) leading to a system of equations to be solved. Here, the computa-tional effort for one time step is high and there is no limit to time step size for stability reasons. Further, the order of consistency is also only O (4t). If we compute e between the previous time level n and the new time level n+1, the Crank-Nicholson Method is obtained with the Crank-Nicholson factor θ with 0 < θ < 1:
en+1− en
4tn = θAe
The computational effort per time step is high and the method is also im-plicit. For 0.5 < θ < 1, stability is ensured; for 0 < θ < 0.5 this is not
the case. If θ equals to 0.5, the order of consistency is O (4t2), otherwise
only O (4t). The value of θ is chosen between 0.55 and 0.65, especially on unstructured grids, where values for θ closer to 0.5 may lead to stability problems (HINKELMANN, 2005).
Finally, it is pointed out, that explicit schemes have to fulfill stability criteria Courant-Friedrichs-Lewy or CFL condition must be ensured, i.e. for shallow
water flow the Courant number Cr must be smaller than one:
Cr = (c + v)
4x/4t ≤ 1 ⇐⇒ 4t ≤
4t
v + c (2.8)
v denotes the flow velocity and c =√gh the wave velocity with g and h being
the gravitational acceleration and the water depth respectively, and 4x is a characteristic element length.
For shallow water transport, the Courant number also must be smaller than one.
Cr = v
4x/4t≤ 1 (2.9)
For shallow water flow and transport, the Neumann criterion also must be fulfilled, i.e. the Neumann number Ne must be smaller than 0.5:
N e = νt.4t 4x2 ≤ 0.5 ⇐⇒ 4t ≤ 4x 2 νt (2.10) N e = νT.4t 4x2 ≤ 0.5 ⇐⇒ 4t ≤ 4x 2 νT (2.11)
Here νtand νT denote the turbulent viscosity and diffusivity, respectively.
Moreover, implicit schemes are not subject to CFL limitations and they are stable regardless of the time-step size which should be limited for accuracy reasons. The accuracy of implicit schemes can be higher as the computational
effort when compared with explicit schemes.CHAPTER 3. EFFICIENT NUMERICAL METHODS 58
h fully implicit x k k+1 t x h h h h explicit y i-1,j i,j k k k k k i,j+1 k+1 i,j i+1,j i,j-1 h y
Figure 3.3: Explicit and (fully) implicit methods, after HELMIG (2000 [99])
As this method is also implicit, the computational effort per time step is high. For 0:5θ1,
stability is ensured; for 0<θ<0:5 this is not the case. Ifθequals 0:5, the order of consistency
is O(∆t
2
), otherwise only O(∆t)(see sec. 3.1.3). For many practical applications, especially on
unstructured grids,θis chosen between 0:55 and 0:65, because values forθcloser to 0:5 may
lead to stability problems.
If further terms of the Taylor series (see sec. 3.1.3) are taken into account, the accuracy is
increased. If the term with∆t2is introduced, the Lax-Wendroff Method, which is explicit and
has an order of consistency O(∆t
2 ), is obtained: en+1 =(1+∆t n A+ 1 2(∆t n ) 2 A2)e n , en+1 en ∆tn =Ae n + 1 2∆t n A2en (3.12)
If the terms with up to∆t4are taken into account, an explicit Runge-Kutta Method , which has
an order of consistency O(∆t
4
), is determined:
1 1 1
Figure 2.2: Explicit and (fully) implicit methods, after (HELMIG, 2000; HINKELMANN, 2005)
2.3.3 Space discretisation
Discretization can be accomplished under three different methods: Finite-Difference Methods FDM, Finite-Element Methods FEM and Finite-Volume Methods FVM:
Finite-Difference Methods: FDM is the oldest method for numerical
18th century. Richardson presented a paper on the first FDM solution in 1910, at the Royal Society of London. The first applications of FDM in CFD were from Courant, Friedrichs and Lewy in 1928 and from Lax and Wendroff in 1960 (CHUNG, 2010). Finite-Difference Methods depend on Taylor series expansions to define the value taken by a variable (e.g. h, u, v) at a given point, as a function of the values at neighbouring points and of local deriva-tives of increasing orders. These Taylor series are then combined to produce approximate expressions for the derivatives involved in the shallow water or transport equations, as a function of a finite number of neighbouring point values (ZIELKE et al., 1999).
The computational domain in FDM is discretized by cells, these cells can be either rectangular or quadrilateral. The unknown variables are defined in the nodes, which are placed either at the centers of the elements or at the intersection points of element boundaries (see fig. 2.3, left), (HELMIG, 2000; HINKELMANN, 2005).
In the FDM the differential quotients are substituted by difference quotients. There are three different ways to obtain the first derivative of a function e, these ways are called: forward (FD), backward differencing (BD) and central differencing (CD): F D : ∂ei,j ∂x = ei+1,j − ei,j 4x − 4x 2 ∂2ei,j ∂x2 − 4x2 6 ∂3ei,j ∂x3 − 4x3 24 ∂4ei,j ∂x4 − O 4x5 (2.12) ⇐⇒ F D : ∂ei,j ∂x = ei+1,j− ei,j 4x + O (4x) (2.13) BD : ∂ei,j ∂x = ei,j− ei−1,j 4x + O (4x) (2.14) CD : ∂ei,j ∂x = ei+1,j − ei−1,j 24x + O 4x2 (2.15)
The second derivative is then:
∂2e
i,j
∂x2 =
ei+1,j− 2ei,j+ei,,−1,j
4x2 + O
4x2 (2.16)
We note, that the FD and BD is of first order or O(4x), while the CD and
the second derivative is of second order or O(4x2) (HINKELMANN, 2005).
FDM has the following characteristics:
1. FDM is easy to formulate and comparatively easy to translate into a computer program, it is also well suitable to show principle numerical effects.
2. It is also the easiest method to use for simple geometries (structured, rectangular mesh). This restriction to simple geometries is a signifi-cant disadvantage in regard of complex boundaries and complex inner structures.
3. FDM does not necessarily guarantee mass or momentum conservation. 4. Harmonic averaging is possible.
5. There are often problems with boundary conditions in 2D and 3D at in-terfaces with jumps in coefficients (more important in subsurface flow). 6. Due to the disadvantages it is hardly chosen today for space and time discretization. However, in the context of semi-discrete FEM and FVM it is widely used for the time discretisation.
Finite-Element Methods: In the FEM the domain is broken into a set of
finite elements that are in the majority of cases unstructured. In 2D triangles or quadrilaterals can be chosen, while in 3D, tatraherda or hexahedra are often used. The unknowns are defined at the nodes, at the center of elements or the center of the edges (SCHINDLER, 2001; JANKOWISKY, 1999). The first applications of FEM in CFD were undertaken by Zienkiewicz and Cheung 1965 (CHUNG, 2010). The principle of the Finite-Element Method
is based on two primary methods: an interpolation method and a variational method (weighted residuals), these two methods deal with the transforma-tion from the continuous domain to a discrete domain. In the interpolatransforma-tion method the unknowns are discretized, while the equation is discretized in the variational method.
The general idea is to replace an unknown function u, which belongs to an in-finite dimensional space, by an approximation defined on a in-finite dimensional space.
The special property of FEM is that the equations should be multiplied with a weighting function and then they are integrated over the domain. Moreover, the use of the interpolation functions for the approximation of the solution over each element provides the continuity of the solution within element boundaries. (SCHINDLER, 2001; HINKELMANN, 2005)
FEM has the following characteristics:
1. Concerning the mathematical derivation, FEM is a complex method. 2. FEM is well suitable for problems with complex boundaries and
com-plex inner structures.
3. FEM guarantees global mass and momentum conservation which is a consequence of the formulation, because the latter states that the residual should vanish in the whole computational domain.
4. Harmonic averaging is not possible in this method, only arithmetic averaging.
5. Using the shape functions it is possible to compute the unknowns in each point of the domain.
6. In recent years it has lost importance in the context of flow simulation (as FVM is superior in same aspects)
Finite-Volume Methods: In the FVM the space is divided into so-called
as in the FEM are chosen. The shallow water flow and transport equations (in conservative form) are integrated over each control volume to produce equations in which volume integral terms of fluxes through the control
vol-ume boundaries occur. At the centroid of each finite volume there is a
computational node at which the variable values are to be calculated. Flux values across a given boundary are used for both control volumes separated by the boundary, resulting in theoretically perfect local mass and momen-tum conservation of the approach, i.e. a flux into a finite volume through a boundary is always equal to a flux out of a neighbouring finite volume through the same boundary. In other words, the basic idea of the FVM consists of integrating the differential equations in their conservative form over all the control volumes. As result, one obtains algebraic equations for each finite volume, in which a number of neighbour nodal values appear (HINKELMANN, 2005; ZIELKE et al., 1999).
FVM are widely used in modeling of shallow water flow with increasing tendency in the last years. This can be explained by their advantages in terms of local conservation, geometric flexibility, good approximation of sharp front and harmonic averaging. (HINKELMANN, 2005; SCHANKAT, 2009).
FVM has the following characteristics:
1. Concerning the mathematical derivation, FVM is more complex than FDM and less complex than FEM.
2. FVM is well suitable to solve problems with complex boundaries and complex inner structure.
3. FVM guarantees local (in the control volume) and global conservation.
4. Harmonic averaging is possible.
boundary of element
boundary of domain element
FEM, FVM
Figure 2.3: Discretisation methods (after HINKELMANN (2005))
2.3.4 Initial and boundary conditions
Initial conditions: The initial conditions describe the state of the model at
the start of simulation, they are defined as:
e (x, y, z, t = 0) = e0(x, y, z) (2.17)
with e0 being the initial values of the entity (e.g. water level)
(HINKEL-MANN, 2005).
Initial conditions must be set at each node of the computational grid for each primary variable, here: water level, velocity components and concentration.
Boundary conditions: Several boundary conditions can be applied to
im-pose the flow or a tracer to enter and exit computational domains. They describe the interaction of the computational domain with its surrounding. The boundary conditions are specified as one of three types:
- Dirichlet-boundary condition: If this condition is chosen, the value of a
