folgt: 8.2.4 Horizontale Impulsdiffusion und hinauf: 8.2 Verfahren nach CASULLI vorher: 8.2.2 Konvektion


8.2.3 3D-Erweiterung

Der von CASULLI und CHENG [20] 1992 vorgenommenen Erweiterung für 3D-Strömungen in flachen Oberflächengewässern liegt die hydrostatische Druckverteilung zugrunde. Diese basiert auf der Annahme, dass der vertikale Druckverlauf allein aus dem Gewicht des Wassers herrührt. Um damit die vertikale Komponente der Impuls-Gl. (6.7) erfüllen zu können, muss angenommen werden, dass die vertikalen Geschwindigkeiten und ihre Änderungen (Beschleunigungen) hinreichend klein sind. Die hydrostatische Druckverteilung hat zur Folge, dass sich bei bekannter Lage des Wasserspiegels der Druck an einer beliebigen Stelle im Wasservolumen aus der Tiefe unter der Wasseroberfläche ergibt: \(p=p_a + \rho \cdot g \left( z_s-z \right) \) mit \(p_a\) Luftdruck an der Wasseroberfläche. Damit kann die Berechnung der Lage des Wasserspiegels als Kern des Verfahrens beibehalten werden.

Die Art der Diskretisierung und die hier verwendeten Bezeichnungen zeigt Bild 15. In das Fluidvolumen werden M horizontale Levels hineingelegt. Auf jedem Level werden die horizontalen Geschwindigkeitskomponenten und die Konzentrationen einschließlich der Dichte diskretisiert. Eine Schicht endet in vertikaler Richtung entweder auf der Mitte zum darüber bzw. darunter liegenden Level oder am Wasserspiegel bzw. Boden. Die vertikalen Geschwindigkeitskomponenten und die Turbulenzgrößen einschließlich der Wirbelviskosität sind auf den Schichträndern diskretisiert. Aus einem Knoten, in der 2D-Berechnung ein Punkt in der Ebene, wird nun eine vertikale Linie. Ebenso wird mit dem Punkt auf der Kantenmitte verfahren. Die diskretisierte Impulsbilanz in Gl. (8.1) umfasst noch die gesamte Wassertiefe und operiert mit einer tiefengemittelten Geschwindigkeit. In der 3D-Berechnung wird für jede der M Schichten eine eigene Bilanz der horizontalen Impulskomponenten wie folgt aufgestellt:

\begin{displaymath} \frac{v_{p.k.j}-v_{p.k-1.j.u}}{\Delta t}= \end{displaymath}


\begin{displaymath} - g \frac{z_{s.k.e}-z_{s.k.a}}{\Delta s} +\frac{ \nu_{j+1}\f... ...z_{j}+\Delta z_{j-1}\right)/2} }{\Delta z_{j}} + f_{h.j} + b_j \end{displaymath} (8.5)

mit
v Geschwindigkeit,
\( \Delta t \) Zeitschrittweite,
g Fallbeschleunigung,
z geodätische Höhe,
\( \Delta s \) Kantenlänge, Ortsschrittweite,
\( \Delta z \) Schichtdicke,
\(\nu\) Wirbelviskosität an der Unterseite der indizierten Schicht,
f Reibung, f friction,
b Einfluss der variablen Dichte, b buoyancy,
Indizes  
p, n Geschwindigkeitskomponente parallel bzw. normal
  zur Kantenrichtung,
k, k-1 aktueller Zeitschritt, vorangegangener Zeitschritt,
j Nummerierung der jeweiligen Schicht 1...M,
u Größe am Ursprung der Bahnlinie,
s Wasserspiegelhöhe, s surface,
a, e Knoten am Anfang bzw. am Ende der Kante und
h horizontal.


Die explizite Diskretisierung der horizontalen Reibung \(f_h\) erläutert Abschnitt 8.2.4. Dichtegetriebene Strömungen (\(b_j\)) werden in dieser Arbeit nicht untersucht.

Es werden nun die Geschwindigkeiten \( v_{p.k.j} \) in allen Schichten an einer Kante zu einem Vektor zusammengefasst. Die Impulsbilanzen in allen Schichten an einer Kante können nach dem Vektor der Geschwindigkeiten aufgelöst werden. Dazu muss eine Matrix invertiert werden, deren Nebendiagonalglieder aus der vertikalen Reibung resultieren, d. h. die horizontalen Geschwindigkeiten in verschiedenen Schichten sind lediglich über die Reibung miteinander verknüpft. Die vertikale Reibung wird also implizit diskretisiert. Bei der Schicht, die mit dem Boden in Berührung kommt, entfällt die Reibung mit der unter ihr liegenden Schicht, anstatt dessen wird die Sohlreibung angesetzt [20]. Entsprechend lässt sich an der Wasseroberfläche eine vom Wind induzierte Schubspannung anbringen.

Nachdem die Impulsbilanzen an allen Kanten nach den Vektoren der Geschwindigkeiten aufgelöst worden sind, können diese in die Massenbilanz der Zelle eingesetzt werden. Das führt dann wieder zu dem folgenden Gleichungssystem für den Wasserspiegel [20]:
\begin{displaymath} A \frac{z_{s.k}-z_{s.k-1}}{\Delta t} + \sum_{i=1}^{N} \sum_{j=1}^{M} v_{p.k.i.j} \cdot \Delta z_{i.j} \cdot d_{i} =0 \end{displaymath} (8.6)

mit
A Zellfläche, umrandet von \(d_1 \dots d_N \),
N Anzahl der mit dem Knoten verbunden Kanten,
  der Zellrandabschnitte,
M Anzahl der horizontalen Schichten und
d Länge des Zellrandabschnitts (s. Bild 13).
Indizes  
i Nummerierung der Kante, des Zellrandabschitts 1...N.
j Nummerierung der jeweiligen Schicht 1...M,


Bei Schichten, in die der Wasserspiegel oder der Boden hineinragt, verringert sich das jeweilige \( \Delta z_{j} \). Wenn die Schicht komplett über dem Wasserspiegel oder unter dem Boden liegt, wird das jeweilige \( \Delta z_{j} \) Null gesetzt. Damit werden ,,trockene`` Schichten deaktiviert. Der Watt-Algorithmus (s. Abschnitte 8.2.1 und 9.2.2) wird dadurch in der 3D-Berechnung fortgesetzt und funktioniert ebenso unproblematisch.

Zur Berechnung der vertikalen Geschwindigkeitskomponente \(v_z\) wird ausschließlich die Massenerhaltung Gl. (6.5) verwendet. Dazu wird eine Zelle in Schichten unterteilt. Die vertikale Geschwindigkeitskomponente ist dort diskretisiert, wo die Knoten-Vertikale die Ober- und Unterseite der Schicht durchstößt. Die Vertikalgeschwindigkeit auf der Oberseite berechnet sich dann wie folgt:
\begin{displaymath} v_{z.j+1} = v_{z.j} - \frac{1}{A}\sum_{i=0}^{N} v_{p.i.j} \cdot d_i \cdot \Delta z_{j} \end{displaymath} (8.7)

mit
\(v_{z.j+1}\) Vertikalgeschwindigkeit auf der Oberseite von Schicht j,
\(v_{z.j}\) Vertikalgeschwindigkeit auf der Unterseite von Schicht j und
\(v_{p.i.j}\) horizontaler Abfluss aus der betrachteten Zelle in Schicht j
  über Zellrandabschnitt i.


Die Vertikalgeschwindigkeit unmittelbar auf dem Boden ergibt sich aus den Horizontalgeschwindigkeiten und der Annahme, dass die Geschwindigkeit parallel zum Boden verläuft. Bei waagerechtem Boden ist die Vertikalgeschwindigkeit demnach Null.

Die vertikale Geschwindigkeitskomponente ist im hier implementierten Verfahren mit keinerlei Trägheit verknüpft. Bei entsprechenden Änderungen der horizontalen Geschwindigkeitskomponenten können vertikale Bewegungen ohne Verzögerung in Gang gesetzt werden. Dies lässt sich als Ursache für Instabilitäten angeben, die in den Testrechnungen in den Abschnitten 9.3.3 und 9.3.4 auftreten. In beiden Beispielen werden Instabilitäten beobachtet, wenn das Problem 3D-diskretisiert wird. Die starken vertikalen Geschwindigkeitsgradienten führen in diesen Beispielen zu nicht mehr vernachlässigbaren vertikalen Geschwindigkeiten. Erst eine Vergleichmäßigung der Randbedingungen erlaubt konvergente Berechnungsergebnisse.

folgt: 8.2.4 Horizontale Impulsdiffusion und hinauf: 8.2 Verfahren nach CASULLI vorher: 8.2.2 Konvektion

Jens WYRWA * 2003-11-05