26 Auswahl einer zukunftsfähigen horizontalen Diskretisierung

Die Vielfalt an numerischen Methoden für die Diskretisierung der herrschenden Gleichungen ist groß und der Aufwand, der damit verbunden ist, eine solche Diskretisierung bis hin zur operationellen Anwendbarkeit zu implementieren, ebenfalls. Daher ist es wichtig, eine gute Entscheidung über die Grundstruktur des in Teil VII zu formulierenden dynamischen Kerns zu treffen, bevor man mit der weiteren Entwicklung und Programmierung fortfährt.

26.1 Aus dem CFL-Kriterium folgende Beschränkungen

Das CFL-Kriterium lautet

\[ \begin{align} \Delta t \leq \frac{\Delta r_\text{min}}{c_\text{max}}, \end{align} \]

hierbei sind $\Delta r_\text{min}$ der minimale Gitterpunktabstand und $c_\text{max}$ die maximale Phasengeschwindigkeit des Gleichungssystems. Daher ist es für die Effiezienz des Modells bei feiner werdender Auflösung nützlich, wenn das Verhältnis

\[ \begin{align} r \coloneqq \frac{\Delta r_\text{max}}{\Delta r_\text{min}} \end{align} \]

von maximalem zu minimalem Gitterpunktabstand bei feiner werdender Auflösung bei einem Wert nicht viel größer als Eins beschränkt ist. Ein solches Gitter bezeichnet man als quasi-uniform. Auf dem Länge-Breite-Gitter mit Winkelauflösung $\phi$ gilt

\[ \begin{align} r = \frac{\phi}{\phi\sin\left(\phi\right)} = \frac{1}{\sin\left(\phi\right)}. \end{align} \]

Das Länge-Breite-Gitter muss an dieser Stelle also ausgeschlossen werden.

26.1.1 Platonische Körper

Es ist also notwendig, dass das in Teil VII zu entwickelnde Atmosphärenmodell auf einem quasi-uniformen Gitter formuliert wird. Von tatsächlicher Uniformität spricht man, wenn jeweils zwei Flächen paarweise kongruent sind. Die Körper, welche durch diese Gitter aufgespannt werden, bezeichnet man als platonische Körper. Die Anzahl der Ecken eines solchen Körpers bezeichne man mit $V$, die Anzahl der Kanten mit $E$ und die Anzahl der Flächen mit $F$. Nun geht man davon aus, dass alle Flächen $v$ Ecken und somit auch Kanten haben. Dann gelten

\[ \begin{align} V &= \frac{Fv}{f},\tag{26.4}\label{eq:platon_0}\\ E &= \frac{Fv}{2},\tag{26.5}\label{eq:platon_1}\\ 2E &= fV,\tag{26.6}\label{eq:platon_2}\\ \end{align} \]

mit $f$ als Anzahl der Kanten bzw. Flächen, welche sich an einer Ecke treffen. Glg. (26.6) kann man auch erhalten, indem man die Glg.en (26.4) und (26.5) durcheinander dividiert, Glg. (26.6) enthält daher keine neue Information. Dies gilt soweit noch für alle Gitter, welche nur aus v-Ecken aufgebaut sind. Nun muss man noch eine Gleichung aufstellen, welche das Gitter auf Uniformität festlegt.

(Tetraeder, Würfel, Oktaeder, Dodekaeder, Ikosaeder)

Es existieren also nur fünf uniforme Gitter auf der Kugel. Diese entstehen aus den platonischen Körpern, indem man ihre Ecken vom Mittelpunkt auf die Kugeloberfläche projiziert. All diese Gitter sind zu grob für ein Atmosphärenmodell und können daher nicht direkt verwendet werden. Darüber hinus folgt aus dieser Tatsache, dass alle verwendbaren Gitter irgendwo Punkte, Linien oder Flächen haben, an denen die lokale Gitterstruktur von der allgemeinen Struktur abweicht. An diesen Stellen sind die Fehler der Diskretisierung ebenfalls anders, was zu einer Sichtbarkeit des Gitters in der Lösung führen kann oder Wellen erzeugen kann, die als Lärm durch die Lösung propagieren. Dies bezeichnet man als grid imprinting.

26.1.2 Grid imprinting

Quasi-uniforme Gitter können beispielsweise erzeugt werden, indem man von einem der platonischen Körper ausgeht und anschließend die Flächen durch sukzessive Bisektionen unterteilt. Als Ausgangspunkt hierfür wählt man man üblicherweise den Ikosaeder, da dieser die meisten Ecken hat und daher das homogenste Gitter generiert. Durch Bisketionen der dreieckigen Flächen entstehen weitere Dreicke, die ein sogenanntes Voronoi-Gitter bilden. Die Voronoi-Zentren dieser Zellen erhält man als Schnittpunkte der Mittelsenkrechten auf den Kanten. Verbindet man diese Punkte, erhält man zwölf Fünnfecke (an den Ecken des Ikosaeder) und eine von der Anzahl der Bisektionen abhängige Anzahl von Sechsecken.

26.2 Aus Rotationssymmetrie folgende Beschränkungen

Die herrschenden Gleichungen sind lokal rotationssymmetrisch bezüglich $\mathbf{k}$. Das bedeutet, dass bei Rotation um diese Achse um jeden Winkel $\phi \in \mathbb{R}$ die Gleichungen wieder in sich selbst überführt werden. Dies möglichst gut auch für das Gitter zu erreichen ist daher erstrebenswert. Damit ist gemeint, dass nach Rotation um einen möglichst kleinen Winkel $\phi$ das Gitter wieder in sich selbst überführt werden soll. $2\pi/\phi$ ist die Zähligkeit dieser Drehachse. Diese sollte möglichst hoch sein. Dies ist für sphärische nicht-uniforme Gitter nicht zu erreichen. Jedoch ist es sinnvoll, zu versuchen, sich der Isotropie (Rotationssymmetrie) möglichst anzunäheren. Dies impliziert zum Beispiel, dass die Standardabweichung der Kantenlängen der Zellen in Relation zur mittleren Länge klein sein sollte.

26.3 Ausschluss spektraler Formulierungen

Entwickelt man horizontale Felder wie in Abschn. 24.4 beschrieben nach Kugelflächenfunktionen, ist es aus zwei Gründen notwendig, Transformationen zwischen Spektral- und Realgitter durchzuführen:

1. Felder mit steilen Gradienten, wie zum Beispiel Feuchtegrößen, müssen im Realraum belassen werden, da sonst zu viele Gibbs'sche Phänomene auftreten würden.
2. Für das Herausschreiben der Daten muss auf ein Realgitter tranformiert werden.

Die Transformation zwischen Spektral- und Realgitter ist die in Abschn. besprochene Legendre-Transformation. Diese skaliert $\propto N^2$, wohingegen der numerische Aufwand der horizontalen Numerik eines Gitterpunktmodells $\propto N$ skaliert, hierbei ist $N$ die Anzahl der horizontalen Gitterpunkte. Bei höher werdender Auflösung werden Spektralmodelle also zunehmend ineffizient und werden daher an dieser Stelle ausgeschlossen.

26.4 Arakawa-Gitter

Diskretisiert man ein zweidimensionales Vektorfeld auf ein ebenes Gitter, so müssen keineswegs alle Variablen an den gleichen Punkten ausgewertet werden. Die unterschiedlichen Größen können an unterschiedlichen Stellen der Polygone definiert werden. Fünf Arten, dies zu tun, sind die Arakawa-Gitter [4], die in Abb. 26.1 dargestellt sind.

Hier werden zunächst nur das A-Gitter und das C-Gitter eindimensional untersucht. Hierzu geht man vom System der linearisierten Flachwassergleichungen Glg.en (13.173) - (13.174) ohne Corioliskraft aus:

\[ \begin{align} \frac{\partial h}{\partial t} &= -H\frac{\partial u}{\partial x}\tag{26.8}\label{eq:arakawa_deriv_0}\\ \frac{\partial u}{\partial t} &= -g\frac{\partial h}{\partial x}\tag{26.9}\label{eq:arakawa_deriv_1} \end{align} \]

Als Ansatz macht man hier

\[ \begin{align} h = \newhat{h}\exp\left(ikx - i\omega t\right) & {} & u = \newhat{u}\exp\left(ikx - i\omega t\right) \end{align} \]

mit $\newhat{h}, \newhat{u} \in \mathbb{C}$. Dies führt auf

\[ \begin{align} -i\omega\newhat{h} = -Hik\newhat{u} \Leftrightarrow \omega\newhat{h} - Hk\newhat{u} = 0, & {} & -i\omega\newhat{u} = -gik\newhat{h} \Leftrightarrow \omega\newhat{u} - gk\newhat{h} = 0. \end{align} \]

In Matrixschreibweise wird dies zu

\[ \begin{align} \left(\begin{array}{cc} \omega & -Hk\\ -gk & \omega \end{array}\right)\left(\begin{array}{c}\newhat{h}\\ \newhat{u}\end{array}\right) = \left(\begin{array}{c} 0\\ 0\end{array}\right) \end{align} \]

Nichttriviale Lösungen existieren, falls die Determinante der Matrix verschwindet:

\[ \begin{align} \omega^2 - gHk^2 = 0 \Rightarrow \omega = k\sqrt{gH} \end{align} \]

Hieraus folgt für die Geschwindigkeiten

\[ \begin{align} c_\text{ph} = \frac{\omega}{k} = \sqrt{gH}, & {} & c_\text{gr} = \frac{\partial\omega}{\partial k} = \sqrt{gH} = c_\text{ph}. \end{align} \]

Dies ist bereits die aus Glg. (16.77) bekannte Dispersionsrelation der Flachwasserwellen.

In der nun folgenden Herleitung der numerischen Versionen dieser Gleichungen wird von zeitlich kontinuierlichen und räumlich diskretisierten Gleichungen ausgegangen, um die Resultate nicht durch ein gewähltes Zeitschrittverfahren zu beeinflussen.

26.4.1 A-Gitter

Beim A-Gitter werden die Oberflächenauslenkung $h$ und die Windgeschwindigkeit $u$ an denselben Orten $j\Delta x$ mit $j\in\mathbb{N}$ definiert. Man notiert hier als Ansatz

\[ \begin{align} h_j = \newhat{h}\exp\left(ikj\Delta x - i\omega t\right), & {} & u_j = \newhat{u}\exp\left(ikj\Delta x - i\omega t\right). \end{align} \]

Die Gleichungen (26.8) - (26.9) diskretisiert man wie folgt:

\[ \begin{align} \frac{\partial h_j}{\partial t} = -H\frac{u_{j + 1} - u_{j - 1}}{2\Delta x} & {} & \frac{\partial u_j}{\partial t} = -g\frac{h_{j + 1} - h_{j - 1}}{2\Delta x} \end{align} \]

Setzt man hier den Ansatz ein, erhält man

\[ \begin{align} -i\omega\newhat{h} &= -\frac{H\newhat{u}}{2\Delta x}\left[\exp\left(ik\Delta x\right) - \exp\left(-ik\Delta x\right)\right]\nonumber\\ \Leftrightarrow -i\omega\newhat{h} &= -\frac{H\newhat{u}}{2\Delta x}\left[i\sin\left(k\Delta x\right) + i\sin\left(k\Delta x\right)\right]\nonumber\\ \Leftrightarrow -i\omega\newhat{h} &= -\frac{H\newhat{u}i}{\Delta x}\sin\left(k\Delta x\right)\nonumber\\ \Leftrightarrow \omega\newhat{h} - \frac{H\newhat{u}}{\Delta x}\sin\left(k\Delta x\right) &= 0\tag{26.17}\label{eq:arakawa_deriv_2} \end{align} \]

sowie

\[ \begin{align} -i\omega\newhat{u} &= -\frac{g\newhat{h}}{2\Delta x}\left[\exp\left(ik\Delta x\right) - \exp\left(-ik\Delta x\right)\right]\nonumber\\ \Leftrightarrow -i\omega\newhat{u} &= -\frac{g\newhat{h}}{2\Delta x}2i\sin\left(k\Delta x\right) = -\frac{g\newhat{h}}{\Delta x}i\sin\left(k\Delta x\right)\nonumber\\ \Leftrightarrow \omega\newhat{u} - \frac{g\newhat{h}}{\Delta x}\sin\left(k\Delta x\right) &= 0.\tag{26.18}\label{eq:arakawa_deriv_3} \end{align} \]

Notiert man Glg.en (26.17) - (26.18) in Matrixschreibweise, erhält man

\[ \begin{align} \left(\begin{array}{cc} \omega & -\frac{H}{\Delta x}\sin\left(k\Delta x\right)\\ -\frac{g}{\Delta x}\sin\left(k\Delta x\right) & \omega \end{array}\right)\left(\begin{array}{c}\newhat{h}\\ \newhat{u}\end{array}\right) = \left(\begin{array}{c} 0\\ 0\end{array}\right). \end{align} \]

Nichttriviale Lösungen existieren, falls die Determinante der Matrix verschwindet, also für

\[ \begin{align} \omega^2 - gH\frac{\sin^2\left(k\Delta x\right)}{\left(\Delta x\right)^2} = 0, & {} & \Rightarrow\omega = \sqrt{gH}\frac{\sin\left(k\Delta x\right)}{\Delta x}. \end{align} \]

Hieraus folgt für die Phasengeschwindigkeit

\[ \begin{align} c_\text{ph} = c_\text{ph}^{\text{(real)}}\frac{\sin\left(k\Delta x\right)}{k\Delta x}\tag{26.21}\label{eq:arakawa_deriv_6} \end{align} \]

und für die Gruppengeschwindigkeit

\[ \begin{align} c_\text{gr} = \frac{\partial\omega}{\partial k} = c_\text{gr}^{\text{(real)}}\cos\left(k\Delta x\right).\tag{26.22}\label{eq:arakawa_deriv_7} \end{align} \]

Führt man die dimensionslose Geschwindigkeit

sowie die dimensionslose Wellenzahl

ein, so lassen sich die Glg.en (26.21) - (26.22) in der Form

notieren.

26.4.2 C-Gitter

Beim C-Gitter werden die Windgeschwindigkeit $u$ an Orten $j\Delta x$ und die Oberflächenauslenkung $h$ and Orten $\left(j + \frac{1}{2}\right)\Delta x$ mit $j\in\mathbb{N}$ definiert. Man notiert hier als Ansatz

\[ \begin{align} h_{j + \frac{1}{2}} = \newhat{h}\exp\left(ikj\Delta x - i\omega t\right)\exp\left(\frac{ik\Delta x}{2}\right), & {} & u_j = \newhat{u}\exp\left(ikj\Delta x - i\omega t\right). \end{align} \]

Die Gleichungen (26.8) - (26.9) diskretisiert man wie folgt:

\[ \begin{align} \frac{\partial h_{j + \frac{1}{2}}}{\partial t} = -H\frac{u_{j + 1} - u_{j}}{2\Delta x}, & {} & \frac{\partial u_j}{\partial t} = -g\frac{h_{j + \frac{1}{2}} - h_{j - \frac{1}{2}}}{2\Delta x} \end{align} \]

Setzt man hier den Ansatz ein, erhält man

\[ \begin{align} -i\omega\newhat{h}\exp\left(\frac{ik\Delta x}{2}\right) &= -\frac{H\newhat{u}}{\Delta x}\left[\exp\left(ik\Delta x\right) - 1\right]\nonumber\\ \Leftrightarrow i\omega\newhat{h}\exp\left(\frac{ik\Delta x}{2}\right) - \frac{H\newhat{u}}{\Delta x}\left[\exp\left(ik\Delta x\right) - 1\right] &= 0\tag{26.29}\label{eq:arakawa_deriv_4} \end{align} \]

sowie

\[ \begin{align} -i\omega\newhat{u} &= -\frac{g\newhat{h}}{\Delta x}\exp\left(\frac{ik\Delta x}{2}\right)\left[1 - \exp\left(-ik\Delta x\right)\right]\nonumber\\ \Leftrightarrow i\omega\newhat{u} - \frac{g\newhat{h}}{\Delta x}\exp\left(\frac{ik\Delta x}{2}\right)\left[1 - \exp\left(-ik\Delta x\right)\right] &= 0.\tag{26.30}\label{eq:arakawa_deriv_5} \end{align} \]

Notiert man Glg.en (26.29) - (26.30) in Matrixschreibweise, erhält man

\[ \begin{align} \left(\begin{array}{cc} i\omega\exp\left(\frac{ik\Delta x}{2}\right) & -\frac{H}{\Delta x}\left[\exp\left(ik\Delta x\right) - 1\right]\\ -\frac{g}{\Delta x}\left[1 - \exp\left(-ik\Delta x\right)\right] & i\omega \end{array}\right)\left(\begin{array}{c}\newhat{h}\\ \newhat{u}\end{array}\right) = \left(\begin{array}{c} 0\\ 0\end{array}\right). \end{align} \]

Nichttriviale Lösungen existieren, falls die Determinante der Matrix verschwindet, also für

\[ \begin{align} -\omega^2 - gH\frac{\left[\exp\left(ik\Delta x\right) - 1\right]\left[1 - \exp\left(-ik\Delta x\right)\right]}{\left(\Delta x\right)^2} &= 0\nonumber\\ \Rightarrow \omega^2 &= gH\frac{\left[1 - \exp\left(ik\Delta x\right)\right]\left[1 - \exp\left(-ik\Delta x\right)\right]}{\left(\Delta x\right)^2}\nonumber\\ \Rightarrow \omega^2 &= gH\frac{2 - 2\cos\left(k\Delta x\right)}{\left(\Delta x\right)^2} = gH2\frac{1 - \cos\left(k\Delta x\right)}{\left(\Delta x\right)^2}. \end{align} \]

Hieraus folgt für die Phasengeschwindigkeit

\[ \begin{align} c_\text{ph} = c_\text{ph}^{\text{(real)}}\sqrt{2\frac{1 - \cos\left(k\Delta x\right)}{\left(k\Delta x\right)^2}}\tag{26.33}\label{eq:arakawa_deriv_12} \end{align} \]

und für die Gruppengeschwindigkeit

\[ \begin{align} c_\text{gr} = \frac{\partial\omega}{\partial k} = \frac{1}{2\omega}\frac{\partial\omega^2}{\partial k} = gH2\frac{\sin\left(k\Delta x\right)}{2\omega\Delta x} = gH\frac{\sin\left(k\Delta x\right)}{\omega\Delta x} = c_\text{gr}^{\text{(real)}}\frac{\sin\left(k\Delta x\right)}{\sqrt{2 - 2\cos\left(k\Delta x\right)}}.\tag{26.34}\label{eq:arakawa_deriv_13} \end{align} \]

In Termen der dimensionslosen Größen Glg.en (26.23) - (26.24) kann man die Glg.en (26.33) - (26.34) in der Form

notieren. Abb. 26.2 zeigt die Dispersionsrelationen des A-Gitters sowie des C-Gitters.

26.5 Skalare Erhaltungseigenschaften auf dem C-Gitter

Ist $q$ eine Erhaltungsgröße, also eine Größe, für die $\md{q} = 0$ gilt, so kann man für $\newtilde{q} \coloneqq \rho q$ eine Kontinuitätsgleichung herleiten:

\[ \begin{align} \md{q} = \frac{\partial q}{\partial t} + \mathbf{v}\cdot\nabla q &= 0,\\ \frac{\partial\rho}{\partial t} + \mathbf{v}\cdot\nabla \rho + \rho\nabla\cdot\mathbf{v} &= 0,\\ \Rightarrow \frac{\partial\newtilde{q}}{\partial t} + \nabla\cdot \mathbf{j}_{\newtilde{q}} &= 0, \end{align} \]

wobei die Größe $\mathbf{j}_{\newtilde{q}} \coloneqq \newtilde{q}\mathbf{v}$ den Fluss der Größe $q$ bezeichnet. Definiert man

\[ \begin{align} Q \coloneqq \int_A\newtilde{q}d^3r \end{align} \]

als das globale Integral der Größe $q$, so ist $Q$ unter kinematischen Randbedingungen erhalten:

\[ \begin{align} \frac{dQ}{dt} = -\int_A\nabla\cdot \mathbf{j}_{\newtilde{q}}d^3r = -\int_{\partial A}\mathbf{j}_{\newtilde{q}}\cdot d\mathbf{n} \end{align} \]

Zerlegt man $A$ nun in $N \geq 1$ Gitterboxen $A = A_1 \cup A_2 \cup \dots \cup A_N$, so gilt

\[ \begin{align} Q = \sum_{i = 1}^NQ_i \end{align} \]

mit

\[ \begin{align} Q_i\coloneqq\int_{A_i}\newtilde{q}d^3r. \end{align} \]

Somit gilt

\[ \begin{align} \frac{dQ}{dt} = \sum_{i = 1}^N\frac{dQ_i}{dt} = -\sum_{i = 1}^N\int_{\partial A_i}\mathbf{j}_{\newtilde{q}}\cdot d\mathbf{n} = -\sum_{i = 1}^N\sum_{F \in F_i}\mathbf{j}_{\newtilde{q}}\cdot\mathbf{n}_F, \end{align} \]

wobei $F_i$ die Menge aller Seitenflächen der Gitterbox $i$ ist, und $\mathbf{n}_F$ der auf der Fläche $F$ senkrecht stehende Normalenvektor (nach außen zeigend). Dabei wird davon ausgegangen, dass $\partial A_i$ jeweils aus einer Anzahl glatter Flächen besteht, über die einzeln zu summieren ist, da sie durch Kanten getrennt sind. Da jede Fläche $F$ in der Doppelsumme zweimal auftritt, jedoch der Vektor $\mathbf{n}$ dabei in zwei entgegengesetzte Richtungen zeigt, ergibt sich die Summe zu Null, sodass auf einem C-Gitter auch nach der Diskretisierung

gilt. Flächen, die Teil von $\partial A$ sind, stellen dabei trivialerweise kein Problem da, da durch sie unter kinematischen Randbedingungen per Definition überhaupt kein Massenfluss stattfindet. Auch unter der verallgemeinerten Voraussetzung, dass $\mathbf{j}_{\newtilde{q}}$ nicht nur aus $\newtilde{q}\mathbf{v}$, sondern aus weiteren, insbesondere diffusiven Flüssen besteht, gilt dies noch. Angemerkt werden sollte außerdem, dass selbst bei fehlerhafter Berechnung von Gittergeometrien $Q$ eine Erhaltungsgröße bleibt.

26.6 Gradientenfelder auf dem C-Gitter

Seien $\psi$ ein Skalarfeld $\nabla\psi$ der Gradient dieses Skalarfeldes, dies ist ein Vektorfeld. Auf dem C-Gitter wird die Rotation eines Vektorfeldes $\mathbf{v}$ durch Anwendung des Stokes'schen Satzes auf ein Polygon mit Fläche $A$ und Kantenlängen $l_i$ ermittelt. Es gilt also

\[ \begin{align} \mathbf{k}\cdot\nabla\times\mathbf{v} = \frac{1}{A}\sum_{i}\gamma_il_iv_i, \end{align} \]

wobei $v_i$ der Wert von $\mathbf{v}$ an der Kante $i$ ist und $\gamma_i = 1$ ist, falls der Einheitsvektor der Kante in mathmetisch positiver Zirkulationsrichtung um das Polygon zeigt, andernfalls ist $\gamma_i = -1$. Der Einfachheit halber nimmt man o. B. d. A. an, dass $\gamma_i = 1$ für alle $i$ ist, sodass man

\[ \begin{align} \mathbf{k}\cdot\nabla\times\mathbf{v} = \frac{1}{A}\sum_{i}l_iv_i \end{align} \]

notieren kann. Nun setzt man

\[ \begin{align} \mathbf{v} = \nabla\psi \end{align} \]

an, was auf

\[ \begin{align} \mathbf{k}\cdot\nabla\times\mathbf{v} = \frac{1}{A}\sum_{i}l_i\left(\nabla\psi\right)_i\tag{26.48}\label{eq:c-grid_rot_kons_deriv_0} \end{align} \]

führt. Es gilt

\[ \begin{align} \left(\nabla\psi\right)_i = \frac{\psi_{\text{to}\left(i\right)} - \psi_{\text{from}\left(i\right)}}{l_i}. \end{align} \]

Setzt man dies in Glg. (26.48) ein, erhält man

\[ \begin{align} \mathbf{k}\cdot\nabla\times\psi = \frac{1}{A}\sum_{i}l_i\frac{\psi_{\text{to}\left(i\right)} - \psi_{\text{from}\left(i\right)}}{l_i} = \frac{1}{A}\sum_{i}\psi_{\text{to}\left(i\right)} - \psi_{\text{from}\left(i\right)}.\tag{26.50}\label{eq:c-grid_rot_kons_deriv_1} \end{align} \]

Sortiert man die Kanten in mathematisch positiver Zirkulationsrichtung um das betrachtete Polygon, gilt

\[ \begin{align} \psi_{\text{to}\left(i\right)} = \psi_{\text{from}\left(i + 1\right)}. \end{align} \]

Jeder Wert von $\psi$ tritt in Glg. (26.50) also zweimal auf, einmal mit positivem und einmal mit negativem Vorzeichen. Dies gilt unabhängig von der betrachteten Raumrichtung. Somit gilt auf einem C-Gitter

26.7 Skalarprodukt auf dem C-Gitter

Zunächst wird in diesem Abschnitt von einem rechteckigen C-Gitter ausgegangen. Während auf solch einem Gitter die Diskretisierungen von Gradient und Divergenz eindeutig sind, ist dies beim Skalarprodukt nicht der Fall. Als Forderung an die Diskretisierung des Skalarproduktes verwendet man, dass die Identität

\[ \begin{align} \nabla\cdot\left(\psi\mathbf{v}\right) = \mathbf{v}\cdot\nabla\psi + \psi\nabla\cdot\mathbf{v} \end{align} \]

sich auf die Diskretisierung übertragen soll.

26.7.1 Eindimensionaler Fall

Die Notation wird aus Abschn. 26.4.2 mit der Ersetzung $h \to \psi$ übernommen. Die Mittelung von Zellen auf Kanten ist einfach das arithmetische Mittel

\[ \begin{align} \newtilde{\psi}_{j} = \frac{\psi_{j - 1/2} + \psi_{j + 1/2}}{2}, \end{align} \]

da die Kante $j$ sich mittig zwischen den zwei Zellen $j - 1/2$ und $j + 1/2$ befindet. Man rechnet nun

\[ \begin{align} \nabla\cdot\left(\psi u\right) &= \frac{\newtilde{\psi}_{j + 1}u_{j + 1} - \newtilde{\psi}_{j}u_{j}}{\Delta x} = \frac{\left(\psi_{j + 3/2} + \psi_{j + 1/2}\right)u_{j + 1} - \left(\psi_{j + 1/2} + \psi_{j - 1/2}\right)u_{j}}{2\Delta x}\nonumber\\ &= \frac{\psi_{j + 3/2}u_{j + 1} - \psi_{j - 1/2}u_{j} + \psi_{j + 1/2}\left(u_{j + 1} - u_{j}\right)}{2\Delta x}\nonumber\\ &= \frac{\psi_{j + 3/2}u_{j + 1} - \psi_{j - 1/2}u_{j} - \psi_{j + 1/2}\left(u_{j + 1} - u_{j}\right) + 2\psi_{j + 1/2}\left(u_{j + 1} - u_{j}\right)}{2\Delta x}\nonumber\\ &= \frac{\left(\psi_{j + 3/2} - \psi_{j + 1/2}\right)u_{j + 1} - \left(\psi_{j - 1/2} - \psi_{j + 1/2}\right)u_{j} + 2\psi_{j + 1/2}\left(u_{j + 1} - u_{j}\right)}{2\Delta x}\nonumber\\ &= \frac{1}{2}\left[\frac{\psi_{j + 3/2} - \psi_{j + 1/2}}{\Delta x}u_{j + 1} + \frac{\psi_{j + 1/2} - \psi_{j - 1/2}}{\Delta x}u_{j}\right]+ \psi_{j + 1/2}\frac{u_{j + 1} - u_{j}}{\Delta x}\nonumber\\ &\stackrel{!}{=} \psi\nabla\cdot\left(u\right) + u\frac{\partial\psi}{\partial x}. \end{align} \]

Die Produktregel überträgt sich also, falls man das Skalarprodukt $uv$ zweier Vektorfelder $u$, $v$ durch

\[ \begin{align} \left(uv\right)_{j + 1/2} = \frac{u_jv_j + u_{j + 1}v_{j + 1}}{2} \end{align} \]

definiert.

26.7.2 Dreidimensionaler rechteckiger Fall

Summiert man die Herleitung im vorherigen Abschnitt über alle drei Raumrichtungen, erhält man die analoge Aussage für den dreidimensionalen Fall mit einem Skalarprodukt

\[ \begin{align} \mathbf{u}\cdot\mathbf{v} = \sum_{e\in c}\frac{u_ev_e}{2}, \end{align} \]

wobei $c$ die Menge aller Seitenflächen $e$ der betrachteten Box ist.

26.7.3 Zweidimensionales orthogonales Gitter

Sei $A_c$ der Flächeninhalt der Zelle $c$, $e$ eine Kante der Zelle, $o\left(e\right)$ die an die Kante $e$ grenzende Nachbarzelle von $c$, $d_e$ eine duale Kantenlänge und $l_e$ eine primäre Kantenlänge. Man notiert zunächst

\[ \begin{align} \nabla\cdot\left(\psi\mathbf{v}\right) &= \frac{1}{A_c}\sum_{e\in c}l_e\frac{\psi_{o(e)} + \psi_c}{2}u_e,\\ \psi\nabla\cdot\mathbf{v} &= \frac{\psi_c}{A_c}\sum_{e\in c}l_eu_e. \end{align} \]

Man fordert nun

\[ \begin{align} \mathbf{v}\cdot\nabla\psi &\stackrel{!}{=} \nabla\cdot\left(\psi\mathbf{v}\right) - \psi\nabla\cdot\mathbf{v}\nonumber\\ \Leftrightarrow\mathbf{v}\cdot\nabla\psi &= \frac{1}{A_c}\sum_{e\in c}l_e\frac{\psi_{o(e)} + \psi_c}{2}u_e - \frac{\psi_c}{A_c}\sum_{e\in c}l_eu_e = \frac{1}{A_c}\sum_{e\in c}l_e\frac{\psi_{o(e)} - \psi_c}{2}u_e = \frac{1}{A_c}\sum_{e\in c}\frac{\psi_{o(e)} - \psi_c}{d_e}u_e\frac{l_ed_e}{2}\nonumber\\ &= \sum_{e\in c}\frac{l_ed_e}{2A_c}\frac{\psi_{o(e)} - \psi_c}{d_e}u_e. \end{align} \]

Hieraus leitet man

\[ \begin{align} \left(\mathbf{u}\cdot\mathbf{v}\right)_c = \sum_{e\in c}\frac{l_ed_e}{2A_c}u_ev_e \end{align} \]

als Diskretisierung für das Skalarprodukt ab. Im ebenen Fall ergibt sich für die Summe der Gewichtungsfaktoren

\[ \begin{align} \sum_{e\in c}\frac{l_ed_e}{2A_c} = 2. \end{align} \]

Im gekrümmten Fall, wie z. B. einer Kugeloberfläche, ist dies jedoch im Allgemeinen nicht mehr der Fall.

26.7.4 Dreidimensionales orthogonales Gitter (shallow)

Nun wird die Herleitung im vorherigen Abschnitt auf ein dreidimensionales Gitter verallgemeinert, dessen Schichtdickte $\Delta z_k$ nur vom vertikalen Index $k$ abhängt. Vertikale Vektorkomponenten werden mit $w$ bezeichnet, der Index $u$ bezeichnet den vektoriellen oder skalaren Gitterpunkt über der betrachteten Gitterbox, der Index $l$ den jeweiligen Gitterpunkt darunter. Dies führt auf

\[ \begin{align} \nabla\cdot\left(\psi\mathbf{v}\right) &= \frac{1}{A_c\Delta z_k}\left(\sum_{e\in c}l_e\Delta z_k\frac{\psi_{o(e)} + \psi_c}{2}u_e + A_c\frac{\psi_u + \psi_c}{2}w_u - A_c\frac{\psi_l + \psi_c}{2}w_l\right),\tag{26.63}\label{eq:div_c-grid_shallow_no_terrain}\\ \psi\nabla\cdot\mathbf{v} &= \frac{\psi_c}{A_c\Delta z_k}\left(\sum_{e\in c}l_e\Delta z_ku_e + A_cw_u - A_cw_l\right). \end{align} \]

Man fordert nun wieder

\[ \begin{align} \mathbf{v}\cdot\nabla\psi &\stackrel{!}{=} \nabla\cdot\left(\psi\mathbf{v}\right) - \psi\nabla\cdot\mathbf{v}\nonumber\\ \Leftrightarrow\mathbf{v}\cdot\nabla\psi &= \frac{1}{A_c\Delta z_k}\left(\sum_{e\in c}l_e\Delta z_k\frac{\psi_{o(e)} + \psi_c}{2}u_e + A_c\frac{\psi_u + \psi_c}{2}w_u - A_c\frac{\psi_l + \psi_c}{2}w_l - \psi_c\sum_{e\in c}l_e\Delta z_ku_e - \psi_cA_cw_u + \psi_cA_cw_l\right)\nonumber\\ \Leftrightarrow\mathbf{v}\cdot\nabla\psi &= \frac{1}{A_c\Delta z_k}\left(\sum_{e\in c}l_e\Delta z_k\frac{\psi_{o(e)} - \psi_c}{2}u_e + A_c\frac{\psi_u - \psi_c}{2}w_u + A_c\frac{\psi_c - \psi_l}{2}w_l\right)\nonumber\\ \Leftrightarrow\mathbf{v}\cdot\nabla\psi &= \sum_{e\in c}\frac{l_ed_e}{2A_c}\frac{\psi_{o(e)} - \psi_c}{d_e}u_e + \frac{\Delta z_{k - 1/2}}{2\Delta z_k}\frac{\psi_u - \psi_c}{\Delta z_{k - 1/2}}w_u + \frac{\Delta z_{k + 1/2}}{2\Delta z_k}\frac{\psi_c - \psi_l}{\Delta z_{k + 1/2}}w_u. \end{align} \]

Hieraus leitet man

als Diskretisierung für das Skalarprodukt ab in der shallow atmosphere. Im ebenen Fall ergibt sich für die Summe der Gewichtungsfaktoren

\[ \begin{align} \sum_{e\in c}\frac{l_ed_e}{2A_c} + \frac{\Delta z_{k - 1/2} + \Delta z_{k + 1/2}}{2\Delta z_k} = 3, \end{align} \]

wobei

\[ \begin{align} \Delta z_k &= z_{k - 1/2} - z_{k + 1/2} = \frac{z_{k - 1} + z_k}{2} - \frac{z_{k} + z_{k + 1}}{2} = \frac{z_{k - 1} + z_k - z_{k} - z_{k + 1}}{2}\nonumber\\ &= \frac{z_{k - 1} - z_{k} + z_k - z_{k + 1}}{2} = \frac{\Delta z_{k - 1/2} + \Delta z_{k + 1/2}}{2} \end{align} \]

eingesetzt wurde. Im gekrümmten Fall, wie z. B. einer Kugeloberfläche, ist dies jedoch im Allgemeinen nicht mehr der Fall.

26.7.5 Dreidimensionales orthogonales Gitter (deep)

Nun wird die Herleitung im vorherigen Abschnitt auf die deep atmosphere verallgemeinert, d. h. die shallow-Atmosphere-Approximation wird aufgehoben. Man geht wieder davon aus, dass die Schichtdickte $\Delta z_k$ nur vom vertikalen Index $k$ abhängt. Mit der gleichen Indexnotation wie im vorherigen Abschnitt erhält man

\[ \begin{align} \nabla\cdot\left(\psi\mathbf{v}\right) &= \frac{1}{V_{c, k}}\left(\sum_{e\in c}A_{e, k}\frac{\psi_{o(e)} + \psi_c}{2}u_{e, k} + A_{c, k + 1/2}\frac{\psi_{c, k} + \psi_{c, k - 1}}{2}w_{c, k + 1/2} - A_{c, k + 1/2}\frac{\psi_{c, k - 1} + \psi_{c, k}}{2}w_{c, k + 1/2}\right),\tag{26.69}\label{eq:div_c-grid_deep_no_terrain}\\ \psi\nabla\cdot\mathbf{v} &= \frac{\psi_c}{V_{c, k}}\left(\sum_{e\in c}A_{e, k}u_{e, k} + A_{c, k + 1/2}w_{c, k + 1/2} - A_{c, k + 1/2}w_{c, k + 1/2}\right). \end{align} \]

Die Zellen-Grundfläche $A_{c, k + 1/2}$ ist nun auch vom vertikalen Index abhängig. $A_{e, k}$ ist die vertikale Fläche, deren Teilmenge die Kante $\left(e, k\right)$ ist. Man fordert nun wieder

\[ \begin{align} \mathbf{v}\cdot\nabla\psi &\stackrel{!}{=} \nabla\cdot\left(\psi\mathbf{v}\right) - \psi\nabla\cdot\mathbf{v}\nonumber\\ \Leftrightarrow\mathbf{v}\cdot\nabla\psi &= \frac{1}{V_{c, k}}\Bigg(\sum_{e\in c}A_{e, k}\frac{\psi_{o(e), k} + \psi_{c, k}}{2}u_{e, k} + A_{c, k + 1/2}\frac{\psi_{c, k} + \psi_{c, k - 1}}{2}w_{c, k + 1/2} - A_{c, k + 1/2}\frac{\psi_{c, k - 1} + \psi_{c, k}}{2}w_{c, k + 1/2}\nonumber\\ &- \sum_{e\in c}A_{e, k}\psi_{c, k}u_{e, k} - A_{c, k + 1/2}\psi_{c, k}w_{c, k + 1/2} + A_{c, k + 1/2}\psi_{c, k}w_{c, k + 1/2}\Bigg)\nonumber\\ \Leftrightarrow\mathbf{v}\cdot\nabla\psi &= \frac{1}{V_{c, k}}\left(\sum_{e\in c}A_{e, k}\frac{\psi_{o(e), k} - \psi_{c, k}}{2}u_{e, k} + A_{c, k + 1/2}\frac{\psi_{c, k - 1} - \psi_{c, k}}{2}w_{c, k + 1/2} + A_{c, k + 1/2}\frac{\psi_{c, k} - \psi_{c, k + 1}}{2}w_{c, k + 1/2}\right)\nonumber\\ \Leftrightarrow\mathbf{v}\cdot\nabla\psi &= \sum_{e\in c}\frac{A_{e, k}d_{e, k}}{2V_{c, k}}\frac{\psi_{o(e), k} - \psi_{c, k}}{d_{e, k}}u_{e, k} + \frac{A_{c, k + 1/2}\Delta z_{k - 1/2}}{2V_{c, k}}\frac{\psi_{c, k - 1} - \psi_{c, k}}{\Delta z_{k - 1/2}}w_{c, k + 1/2}\nonumber\\ &+ \frac{A_{c, k + 1/2}\Delta z_{k + 1/2}}{2V_{c, k}}\frac{\psi_{c, k} - \psi_{c, k + 1}}{\Delta z_{k + 1/2}}w_{c, k + 1/2}. \end{align} \]

Hieraus leitet man

als Diskretisierung für das Skalarprodukt ab.

26.8 Festlegung auf das C-Gitter

Das C-Gitter hat gegenüber den anderen Gittern zusammenfassend drei Vorteile:

1. Es hat vorteilhafte Dispersionsrelationen.
2. Erhaltung skalarer Variablen ist leicht zu gewährleisten.
3. Gradientenfelder sind rotationsfrei.
4. Die Produktregel $\nabla\cdot\left(\psi\mathbf{v}\right) = \mathbf{v}\cdot\nabla\psi + \psi\nabla\cdot\mathbf{v}$ reproduziert sich bei passender Diskretisierung des Skalarproduktes.

Daher wird an dieser Stelle die Entscheidung getroffen, ein C-Gitter für den in Teil VII zu entwickelnden dynamischen Kern zu verwenden.

26.8.1 Folgerungen

Die Vektoren befinden sich bei einem C-Gitter auf den Kanten und schneiden diese orthogonal. Verlängert man die Vektorpfeile in beide Richtungen, landet man an den skalaren Datenpunkten. Das so erhaltene Gitter ist das sogenannte duale Gitter. Gitter, die ein duales Gitter haben, bezeichnet man als orthogonal. Dies ist nicht immer der Fall. Um die positiven Eigenschaften des C-Gitters nutzen zu können, legt man sich an dieser Stelle auf ein orthogonales Gitter für den in Teil VII zu formulierenden dynamischen Kern fest.

26.9 Das Problem des Verhältnisses der Freiheitsgrade

Im kontinuierlichen Fall gelten für $N \geq 1$ prognostische Variablen $\psi_i$ mit $1 \leq i \leq N$ ein System aus $N$ prognostischen Gleichungen. Linearisiert man diese, entsteht ein lineares gekoppeltes partielles Differenzialgleichungssystem. Nimmt man keine externen Forcings auf, ist dieses homogen. Setzt man für jede der prognostischen Variablen eine ebene Welle $\psi_i = \newhat{\psi}_i\exp\left(i\mathbf{k}\cdot\mathbf{r} - i\omega t\right)$ mit $\newhat{\psi}_i\in\mathbb{C}$ ein, entsteht homogenes lineares Gleichungssystem für die komplexen Amplituden $\newhat{\psi}_i$. Dies ist ein Eigenwertproblem für die Frequenz $\omega$. Die Eigenwerte findet man über Nullsetzen des charakteristischen Polynoms (der Determinante) dieses Gleichungssystems. Das charakteristische Polynom hat nach dem Hauptsatz der Algebra $N$ komplexe, allerdings nicht notwendigerweise verschiedene Nullstellen. Jede der verschiedenen Nullstellen entspricht einem Zweig der Dispersionsrelation.

Bei der Diskretisierung ist dies insbesondere beim Staggering von Vektorfeldern wichtig. Im Falle der linearisierten shallow water equations hat man drei Gleichungen für die drei Variablen $\left(u, v, h\right)$, man hat also „ doppelt so viele Pfeile wie Punkte“. Hat ein Gitter zu wenig Pfeile, so kann die Dispersionsrelation nicht mehr die nötige Anzahl an Moden beinhalten. Hat man zu viele Pfeile, so könnte die Dispersionsrelation zusätzliche, unphysikalische Zweige beinhalten. Dies kann jedoch durch Anwendung spezieller Mittelungsoperatoren in den prognostischen Gleichungen verhindert werden.

26.9.1 Folgerungen

Auf dem dreieckigen Gitter hat man pro skalarem Gitterpunkt 1,5 Vektorpunkte, wohingegen das korrekte Verhältnis 2 lautet. Solche Gitter werden daher ausgeschlossen. Viereckige Gitter sind in dieser Hinsicht optimal. Bei Gittern mit mehr als vier Ecken muss man die Überspezifikation durch $m$ algebraische Bedingungen an Vektorfelder (man könnte auch sagen: diagnostische Gleichungen oder Zustandsgleichungen) pro skalarem Gitterpunkt eliminieren. Bei $N$ Ecken gilt

\[ \begin{align} m = \frac{N - 4}{2}, \end{align} \]

da zwei Zellen sich jeweils eine Kante teilen. Bei ungeraden $N$ liegt $m$ in der Mitte zwischen zwei natürlichen Zahlen, sodass in diesem Fall eine diagnostische Vektorgleichung für zwei Zellen gelten muss, was die Isotropie des Gitters zerstört. Daher ist es notwendig, dass $N$ gerade ist. Weiterhin möchte man möglichst nur eine algebraische Zusatzbedingung an Vektorfelder pro Zelle aufstellen, daher ist die erste Wahl für $N$ 4 und die zweite Wahl 6.

Die in diesem Abschnitt besprochene Problematik wird im nächsten Kapiel für das hexagonale Gitter näher ausformuliert.

26.10 Ausschluss viereckiger Gitter

Das einfachste viereckige Gitter ist das bereits in Abschn. 26.1 ausgeschlossene Länge-Breite-Gitter. Man stellt fest, dass alle bekannten viereckigen Gitter mindestens eines der folgenden grundlegenden Probleme haben, die bereits als Ausschlusskriterium identifiziert wurden:

• Orthogonalität wird verletzt (so z. B. beim cubed-sphere-Gitter)
• Quasi-Uniformität wird verletzt (so z. B. beim Länge-Breite-Gitter).
• Isotropie wird grob verletzt (so z. B. beim kite grid).
• Das Gitter hat Kanten oder Punkte, an denen die Struktur stark von der Gitterstruktur in anderen Regionen abweicht (so z. B. beim reduzierten Länge-Breite-Gitter).

Daher werden viereckige Gitter an dieser Stelle ausgeschlossen.

26.11 Festlegung auf das sechseckige Ikosaedergitter

Das einzige nach diesem Ausschlussverfahren noch übriggebliebene Gitter ist das sechseckige Ikosaedergitter.