25 Zeitschrittverfahren

Sei $\psi: A \subseteq \mathbb{R}^3 \times \mathbb{R} \to \mathbb{R}^n$ der atmosphärische Zustandsvektor. Alle fundamentalen physikalischen Gesetze sind erster Ordnung in der Zeit, daher kann man schreiben

\[ \begin{align} \frac{d\psi}{dt} = F\left(\psi\right), \end{align} \]

wobei $F$ die Physik enthält. Diskretisiert man dies auf ein zeitliches Gitter $n\Delta t$ mit $n \in \mathbb{Z}$ und interpoliert zwischen zwei Zeitschritten $n$, $n + l$, erhält man

\[ \begin{align} \psi_{n + l} - \psi_{n} = \int_{n\Delta t}^{\left(n + l\right)\Delta t}F\left(\psi\right)dt. \end{align} \]

Die Anzahl der dabei involvierten Zeitschritte ist $l + 1$. Man kann nun das Integral auf der rechten Seite auf vielfälitge Art und Weise diskretisieren. Das explizite Euler-Verfahren ergibt sich mit $l = 1$ über die Näherung

\[ \begin{align} \int_{n\Delta t}^{\left(n + 1\right)\Delta t}F\left(\psi\right)dt \approx F\left(\psi_{n}\right)\Delta t \end{align} \]

mit $\psi_{n} \coloneqq \psi\left(n\Delta t\right)$. Das implizite Euler-Verfahren hingegen erhält man über

\[ \begin{align} \int_{n\Delta t}^{\left(n + 1\right)\Delta t}F\left(\psi\right)dt \approx F\left(\psi_{n + 1}\right)\Delta t \end{align} \]

In diesem Kapitel soll verständlich gemacht werden, wie sich die zeitliche Diskretisierung einer linearen Differenzialgleichung auf die Dispersionsrelation und die Stabilität einer Wellenlösung auswirken kann. Es wird die lineare Advektionsgleichung

\[ \begin{align} \frac{\partial f}{\partial t} + U\frac{\partial f}{\partial x} = 0\tag{25.5}\label{eq:test_lineare_instab} \end{align} \]

betrachtet mit $U\not = 0$ als homogener Advektionsgeschwindigkeit und $f$ als irgendeiner Teilcheneigenschaft. Die analytische Lösung lautet in der komplexen Schreibweise

\[ \begin{align} f = \exp\left(i\left(kx - \omega t\right)\right), \end{align} \]

setzt man dies ein folgt $-i\omega + Uik = 0$, also als Dispersionsrelation

\[ \begin{align} \omega = Uk. \end{align} \]

Phasen- und Gruppengeschwindigkeit sind also gleich $U$. Nun werden Ort und Zeit diskretisiert, man schreibt $f_l^{(m)} = f\left(l\Delta x, m\Delta t\right)$ und setzt an

\[ \begin{align} f_{l}^{(m)} = e^{Cm\Delta t}e^{i\left(kl\Delta x - \omega m\Delta t\right)}\tag{25.8}\label{eq:ansatz_time_stepping_1d_analysis} \end{align} \]

mit $f_{l}^{(0)} = e^{ikl\Delta x}$ als Anfangszustand. Hierbei sind $C$, $k$ und $\omega$ reelle Zahlen. $k$ kann hierbei aus Symmetriegründen keinen Imaginärteil haben. Es sollen $C$ und die Dispersion $\omega\left(k\right)$ bestimmt werden.

25.1 Rein explizite Verfahren

25.1.1 Explizites Euler-Verfahren

Zunächst wird folgende Diskretisierung untersucht:

\[ \begin{align} \frac{f_l^{(m + 1)} - f_l^{(m)}}{\Delta t} + U\frac{f_{l + 1}^{(m)} - f_{l - 1}^{(m)}}{2\Delta x} = 0 \end{align} \]

Dies bezeichnet man als explizites Euler-Verfahren. Einsetzen ergibt

\[ \begin{align} \frac{e^{C\Delta t}e^{-i\omega\Delta t} - 1}{\Delta t} + U \frac{e^{ik\Delta x} - e^{-ik\Delta x}}{2\Delta x} = \frac{e^{C\Delta t}\cos\left(\omega\Delta t\right) - ie^{C\Delta t}\sin\left(\omega\Delta t\right) - 1}{\Delta t} + \frac{iU\sin\left( k\Delta x\right)}{\Delta x} = 0. \end{align} \]

Der Realteil ist

\[ \begin{align} e^{C\Delta t}\cos\left(\omega\Delta t\right) = 1, \end{align} \]

hieraus folgt i. A. $C> 0$, es liegt Instabilität vor. Hierbei handelt es sich um lineare Instabilität, da Glg. (25.5) linear in $f$ ist. Der Imaginärteil ist

\[ \begin{align} e^{C\Delta t}\sin\left(\omega\Delta t\right) = \frac{U\Delta t\sin\left(k\Delta x\right)}{\Delta x}. \end{align} \]

Hieraus folgt

\[ \begin{align} \tan\left(\omega\Delta t\right) = \frac{U\Delta t\sin\left(k\Delta x\right)}{\Delta x}, \end{align} \]

also die Dispersion

\[ \begin{align} \omega\left(k\right) = \frac{1}{\Delta t}\arctan\left(\frac{U\Delta t\sin\left(k\Delta x\right)}{\Delta x}\right). \end{align} \]

$\lambda$ sei die Wellenlänge, also $k = \frac{2\pi}{\lambda}$. Für $\lambda\gg \Delta x$ gilt $\sin\left(k\Delta x\right) \approx k\Delta x$, also $\omega \approx \frac{1}{\Delta t}\arctan\left(Uk\Delta t\right)$. Im Fall $U\Delta t\ll \lambda$, also $Uk\Delta t\ll 1$, ist $\omega \approx \frac{1}{\Delta t}U\Delta t k = Uk$ wie im kontinuierlichen Fall. Lange Wellen werden bei hinreichend kleinem Zeitschritt also gut aufgelöst. Bei der kürzesten aufgelösten Welle $\lambda = 2\Delta x$ jedoch folgt $\omega = \frac{1}{\Delta t}\arctan\left(\frac{U\Delta t\sin\pi}{\Delta x}\right) = 0$, die Welle ist also stationär. Abb. 25.1 zeigt die Dispersionsrelation für den Fall $U = 10$ m/s. Man erkennt, dass Wellen mit $k\leq \frac{1}{\Delta x}$, also $\lambda\geq 2\pi\Delta x$ anständig simuliert werden.

25.1.1.1 First-order-upwind-Verfahren

Im Fall $U>0$ bezeichnet man die Diskretisierung

\[ \begin{align} \frac{f_l^{(m + 1)} - f_l^{(m)}}{\Delta t} + U\frac{f_{l}^{(m)} - f_{l - 1}^{(m)}}{\Delta x} = 0 \end{align} \]

als First-order-upwind-Verfahren. Dies führt auf

\[ \begin{align} \frac{e^{C\Delta t}e^{-i\omega\Delta t} - 1}{\Delta t} + U \frac{1 - e^{-ik\Delta x}}{\Delta x} = \frac{e^{C\Delta t}\cos\left(\omega\Delta t\right) - ie^{C\Delta t}\sin\left(\omega\Delta t\right) - 1}{\Delta t} + U\frac{1 - \cos\left( k\Delta x\right) + i\sin\left(k\Delta x\right)}{\Delta x} = 0. \end{align} \]

Der Imaginärteil hiervon ist

\[ \begin{align} &\frac{e^{C\Delta t}\sin\left(\omega\Delta t\right)}{\Delta t} = U\frac{\sin\left(k\Delta x\right)}{\Delta x}\nonumber\\ &\Leftrightarrow e^{C\Delta t}\sin\left(\omega\Delta t\right) = U\frac{\Delta t\sin\left(k\Delta x\right)}{\Delta x}\nonumber\\ &\Leftrightarrow e^{C\Delta t} = U\frac{\Delta t\sin\left(k\Delta x\right)}{\Delta x\sin\left(\omega\Delta t\right)}. \end{align} \]

Die Diskretisierung ist genau dann stabil, wenn $C\leq 0$ ist, also für

\[ \begin{align} U\frac{\Delta t\sin\left(k\Delta x\right)}{\Delta x\sin\left(\omega\Delta t\right)}\leq 1. \end{align} \]

Wegen

\[ \begin{align} c = \frac{\omega}{k}\leq\frac{\Delta x}{\Delta t}\Leftrightarrow\omega\Delta t\leq k\Delta x \end{align} \]

folgt weiter die Abschätzung

Diese Gleichung bezeichnet man als Courant-Friedrichs-Lewy-Kriterium (CFL-Kriterium). Man definiert die Courant-Zahl oder auch CFL-Zahl $c$ durch

\[ \begin{align} c\coloneqq\frac{U\Delta t}{\Delta x}. \end{align} \]

Das CFL-Kriterium lautet sprachlich

25.1.2 Leapfrog-Verfahren

Verwendet man einen zentralen zeitlichen Differenzenquotienten, schreibt also

\[ \begin{align} \frac{f_l^{(m + 2)} - f_l^{(m)}}{2\Delta t} + U\frac{f_{l + 1}^{(m + 1)} - f_{l - 1}^{(m + 1)}}{2\Delta x} = 0, \end{align} \]

so ergibt dies eingesetzt

\[ \begin{align} \frac{e^{2C\Delta t}e^{-2i\omega\Delta t} - 1}{2\Delta t} + U \frac{e^{ik\Delta x}e^{C\Delta t}e^{-i\omega\Delta t} - e^{-ik\Delta x}e^{C\Delta t}e^{-i\omega\Delta t}}{2\Delta x} &= 0 \nonumber\\ \Leftrightarrow \frac{e^{C\Delta t}e^{-i\omega\Delta t} - e^{-C\Delta t}e^{i\omega\Delta t}}{2\Delta t} + U \frac{e^{ik\Delta x} - e^{-ik\Delta x}}{2\Delta x} &= 0. \end{align} \]

Der Realteil ergibt

\[ \begin{align} e^{C\Delta t}\cos\left(\omega\Delta t\right) = e^{-C\Delta t}\cos\left(\omega\Delta t\right). \end{align} \]

Das bedeutet $C = 0$, die Lösung ist also stabil. Der Imaginärteil wird zu

\[ \begin{align} \frac{\sin\left(\omega\Delta t\right)}{2\Delta t}\left(e^{-C\Delta t} + e^{C\Delta t}\right) = \frac{\sin\left(\omega\Delta t\right)}{\Delta t} = U\frac{\sin\left(k\Delta x\right)}{\Delta x}. \end{align} \]

Daraus folgt

\[ \begin{align} \sin\left(\omega\Delta t\right) = U\frac{\Delta t}{\Delta x}\sin\left(k\Delta x\right)\Rightarrow\omega = \frac{1}{\Delta t}\arcsin\left[\frac{U\Delta t}{\Delta x}\sin\left(k\Delta x\right)\right], \end{align} \]

also eine andere Dispersionsrelation als im expliziten Fall. Allgemein scheint die Stabilität der Lösung aus dem Realteil ersichtlich zu sein, während die Dispersion aus dem Imaginärteil zu folgen scheint. Schemata sind weiterhin entweder stabil ($C \leq 0$), labil ($C > 0$) oder bedingt stabil (stabil für hinreichend kleine $\Delta t$). Dies wird vom Zeitschrittverfahren beeinflusst. Der Optimalfall ist $C = 0$, er ergibt sich für das Leapfrogverfahren, daher ist dieses Verfahren das geeignetste.

25.1.3 Allgemeine Multischritt-Verfahren

25.1.4 Runge-Kutta-Verfahren

Sei $f: \mathbb{R}^l\times\mathbb{R} \to \mathbb{R}^m$ mit $l, m \geq 1$ eine unendlich oft stetig differenzierbare Funktion, welche das Anfangswertproblem

\[ \begin{align} \frac{df}{dt} = F\left[f\left(t\right)\right], & {} & f_0 = f_0 \end{align} \]

erfüllt. Sogenannte Runge-Kutta-Verfahren basieren darauf, mehrere Werte für $df/dt$ zu berechnen und so zu so überlagern, dass möglich viele Glieder der Taylor-Reihe abgedeckt werden.

25.1.4.1 Linearer Fall

Die Taylor-Entwicklung von $f$ lautet

\[ \begin{align} f\left(t\right) = \sum_{k = 0}^\infty\frac{f^{(k)}\left(0\right)}{k!}t^k. \end{align} \]

Die Funktion $F$ macht aus einer Funktion $f$ deren zeitliche Ableitung, somit folgt

\[ \begin{align} f\left(t\right) = \sum_{k = 0}^\infty\frac{F^k\left(0\right)}{k!}t^k \end{align} \]

mit $F^0 = f$. Unter der Annahme, dass $F$ weiterhin linear ist, kann man dies zu

\[ \begin{align} f\left(t\right) &= \sum_{k = 0}^n\frac{F^k\left(0\right)}{k!}t^k + \sum_{k = n + 1}^\infty\frac{F^k\left(0\right)}{k!}t^k\nonumber\\ &= f_0 + tF\left\{f_0 + \frac{t}{2}F\left[f_0 + \dots + \frac{t}{n}F\left(f_0\right)\right]\right\} + O\left(t^{n + 1}\right)\tag{25.30}\label{eq:rk_lin} \end{align} \]

umformulieren. Es werden also die ersten $n$ Ableitungen von $f$ durch dieses Verfahren approximiert, was man als lineares Runge-Kutta-Verfahren dritter Ordnung bezeichnet.

25.1.4.2 Nichtlinearer Fall

Definiert man

\[ \begin{align} f_1 \coloneqq \frac{t}{2}F\left[f_0 + \dots + \frac{t}{n}F\left(f_0\right)\right], \end{align} \]

kann man Glg. (25.30) in der Form

\[ \begin{align} f\left(t\right) &= f_0 + tF\left(f_0 + f_1\right) + O\left(t^{n + 1}\right) \end{align} \]

notieren. Ist die Funktion $F$, welche die Zeitentwicklung des Systems festlegt, nichtlinear, gilt nicht mehr

\[ \begin{align} tF\left(f_0 + f_1\right) = tF\left(f_0\right) + tF\left(f_1\right), \end{align} \]

sondern

\[ \begin{align} tF\left(f_0 + f_1\right) = tF\left(f_0\right) + tF'\left(f_0\right)f_1 + tO\left(f_1^2\right). \end{align} \]

modifizieren. Somit gilt in diesem Fall

\[ \begin{align} f\left(t\right) &= f_0 + tF\left(f_0 + f_1\right) + O\left(t^3\right). \end{align} \]

Im nichtlinearen Fall ist eine Diskretisierung in Form von Glg. (25.30) also immer zweiter Ordnung.

25.2 Verfahren mit impliziten Anteilen

25.2.1 Implizites Euler-Verfahren

Verwendet man einen linksseitigen Differenzenquotienten für die Zeitableitung (dies bezeichnet man als implizites Euler-Verfahren), schreibt also

\[ \begin{align} \frac{f_l^{(m + 1)} - f_l^{(m)}}{\Delta t} + U\frac{f_{l + 1}^{(m + 1)} - f_{l - 1}^{(m + 1)}}{2\Delta x} = 0, \end{align} \]

so ergibt dies eingesetzt

\[ \begin{align} & \frac{e^{C\Delta t}e^{-i\omega\Delta t} - 1}{\Delta t} + U \frac{e^{ik\Delta x}e^{C\Delta t}e^{-i\omega\Delta t} - e^{-ik\Delta x}e^{C\Delta t}e^{-i\omega\Delta t}}{2\Delta x} = 0 \nonumber\\ &\Leftrightarrow\frac{1 - e^{-C\Delta t}e^{i\omega\Delta t}}{\Delta t} + U \frac{e^{ik\Delta x} - e^{-ik\Delta x}}{2\Delta x} = 0. \end{align} \]

Der Realteil ergibt

\[ \begin{align} 1 = e^{-C\Delta t}\cos\left(\omega\Delta t\right). \end{align} \]

Das bedeutet $C<0$, die Lösung ist also immer stabil, sie konvergiert jedoch gegen Null. Der Imaginärteil wird zu

\[ \begin{align} \frac{e^{-C\Delta t}\sin\left(\omega\Delta t\right)}{\Delta t} = U\frac{\sin\left(k\Delta x\right)}{\Delta x}. \end{align} \]

Daraus folgt

\[ \begin{align} \tan\left(\omega\Delta t\right) = U\frac{\Delta t}{\Delta x}\sin\left(k\Delta x\right), \end{align} \]

also die gleiche Dispersionsrelation wie beim expliziten Euler-Verfahren.

25.2.2 Crank-Nicolson-Verfahren

Das Crank-Nicolson-Verfahren ist eine Kombination aus explizitem und implizitem Euler-Verfahren:

\[ \begin{align} \frac{f_l^{(m + 1)} - f_l^{(m)}}{\Delta t} + U\frac{f_{l + 1}^{(m)} - f_{l - 1}^{(m)} + f_{l + 1}^{(m + 1)} - f_{l - 1}^{(m + 1)}}{4\Delta x} = 0 \end{align} \]

Setzt man hier Glg. (25.8) ein, erhält man

\[ \begin{align} & \frac{e^{C\Delta t}e^{-i\omega\Delta t} - 1}{\Delta t} + U \frac{e^{ik\Delta x} - e^{-ik\Delta x} + e^{ik\Delta x}e^{C\Delta t}e^{-i\omega\Delta t} - e^{-ik\Delta x}e^{C\Delta t}e^{-i\omega\Delta t}}{4\Delta x} = 0\nonumber\\ &\Leftrightarrow\frac{e^{C\Delta t}\cos\left(\omega\Delta t\right) - ie^{C\Delta t}\sin\left(\omega\Delta t\right) - 1}{\Delta t} + U \frac{2i\sin\left(\omega\Delta t\right) + 2i\sin\left(k\Delta x\right)e^{C\Delta t}e^{-i\omega\Delta t}}{4\Delta x} = 0\nonumber\\ &\Leftrightarrow\frac{e^{C\Delta t}\cos\left(\omega\Delta t\right) - ie^{C\Delta t}\sin\left(\omega\Delta t\right) - 1}{\Delta t} + U \frac{2i\sin\left(\omega\Delta t\right) + 2ie^{C\Delta t}\sin\left(k\Delta x\right)\cos\left(\omega\Delta t\right) + 2e^{C\Delta t}\sin\left(k\Delta x\right)\sin\left(\omega\Delta t\right)}{4\Delta x} = 0. \end{align} \]

Der Realteil ergibt

\[ \begin{align} &\Leftrightarrow\frac{e^{C\Delta t}\cos\left(\omega\Delta t\right) - 1}{\Delta t} + U\frac{2e^{C\Delta t}\sin\left(k\Delta x\right)\sin\left(\omega\Delta t\right)}{4\Delta x} = 0\nonumber\\ &\Leftrightarrow\frac{\cos\left(\omega\Delta t\right) - e^{-C\Delta t}}{\Delta t} + U\frac{\sin\left(k\Delta x\right)\sin\left(\omega\Delta t\right)}{2\Delta x} = 0\nonumber. \end{align} \]

Der Imaginärteil ergibt

\[ \begin{align} &\Leftrightarrow-\frac{ie^{C\Delta t}\sin\left(\omega\Delta t\right)}{\Delta t} + U\frac{2i\sin\left(\omega\Delta t\right) + 2ie^{C\Delta t}\sin\left(k\Delta x\right)\cos\left(\omega\Delta t\right)}{4\Delta x} = 0\nonumber\\ &\Leftrightarrow-\frac{e^{C\Delta t}}{\Delta t} + U\frac{1 + e^{C\Delta t}\sin\left(k\Delta x\right)\cot\left(\omega\Delta t\right)}{2\Delta x} = 0\nonumber. \end{align} \]

Das bedeutet $C < 0$, die Lösung ist also stabil, sie verschwindet jedoch irgendwann. Der Imaginärteil wird zu

\[ \begin{align} \frac{e^{-C\Delta t}\sin\left(\omega\Delta t\right)}{\Delta t} = U\frac{\sin\left(k\Delta x\right)}{\Delta x}. \end{align} \]

Daraus folgt

\[ \begin{align} \tan\left(\omega\Delta t\right) = U\frac{\Delta t}{\Delta x}\sin\left(k\Delta x\right), \end{align} \]

also die gleiche Dispersionsrelation wie beim expliziten Euler-Verfahren.

25.3 Mischung von Zeitschrittverfahren

Man ist nicht darauf beschränkt, eine differenzialgleichung entweder nur explizit oder nur implizit zu lösen. Sei $\newtilde{f}$ das diskretisierte Analogon eines kontinuierlichen Feldes $f$ und seien $\mathbf{r}_i$ die Gitterpunkte. Ein Zwei-Schritt-Verfahren für eine differenzialgleichung erster Ordnung in der Zeit für $f$ lässt sich in der Form

\[ \begin{align} \newtilde{f}\left(\mathbf{r}_i, \left(n + 1\right)\Delta t\right) = \newtilde{f}\left(\mathbf{r}_i, n\Delta t\right) + \Delta t\sum_{j = 1}^N\newtilde{F}_j\left[\newtilde{f}\left(\left(n + m_j\right)\Delta t\right)\right] \end{align} \]

notieren. Hierbei seien $\newtilde{F}_j$ für $j = 1, \dots, N$ die auf $\newtilde{f}$ wirkenden diskretisierten Forcings. Dabei soll gelten

\[ \begin{align} m_j = \begin{cases} 0,\text{ Term $j$ wird explizit behandelt,}\\ 1,\text{ Term $j$ wird implizit behandelt.} \end{cases} \end{align} \]

Es wurde stillschweigend davon ausgegangen, dass $f$ das einzige auftretende Feld ist. Verallgemeinernd kann man $f$ als Tupel von Feldern ansehen, also als Zustandsvektor. Weiterhin kann man $0 \leq m_j \leq 1$ zulassen, hierbei ist $m_j$ der implizite Anteil. Das Crank-Nicolson-Verfahren bedeutet $m_j = \frac{1}{2}$.

25.4 Aspekte partieller Differenzialgleichungen

25.4.1 Forward-backward-Verfahren

25.5 Energetische Selbstkonsistenz im idealen Gas

25.5.1 Advektion von Impuls

25.5.1.1 Advektion kinetischer Energie

Da $\mathbf{v}\cdot\left[\left(\nabla\times\mathbf{v}\right)\times\mathbf{v}\right] = 0$ gilt, ist der verallgemeinerte Coriolis-Term energetisch unerheblich. Für eine selbstkonsistente zeitliche Diskretisierung der Advektion von kinetischer Energie verwendet man daher das Gleichungssystem

\[ \begin{align} \frac{\partial\mathbf{v}}{\partial t} = -\nabla\frac{\mathbf{v}^2}{2}, & {} & \frac{\partial\rho}{\partial t} = -\nabla\cdot\left(\rho\mathbf{v}\right). \end{align} \]

Im kontinuierlichen Fall gilt

\[ \begin{align} \rho\mathbf{v}\cdot\frac{\partial\mathbf{v}}{\partial t} &= -\rho\mathbf{v}\cdot\nabla\frac{\mathbf{v}^2}{2} = -\nabla\cdot\left(\rho\mathbf{v}\frac{\mathbf{v}^2}{2}\right) + \frac{\mathbf{v}^2}{2}\nabla\cdot\left(\rho\mathbf{v}\right) = -\nabla\cdot\left(\rho\mathbf{v}\frac{\mathbf{v}^2}{2}\right) - \frac{\mathbf{v}^2}{2}\frac{\partial\rho}{\partial t}\nonumber\\ \Rightarrow \rho\mathbf{v}\cdot\frac{\partial\mathbf{v}}{\partial t} + \frac{\mathbf{v}^2}{2}\frac{\partial\rho}{\partial t} &= \frac{\partial\left(\rho\mathbf{v}^2/2\right)}{\partial t} = -\nabla\cdot\left(\rho\mathbf{v}\frac{\mathbf{v}^2}{2}\right) \end{align} \]

für die Auswirkung der Flussdichte von kinetischer Energie auf die lokale kinetische Energiedichte. Diese Umformungsschritte müssen sich auf die Diskretisierung übertragen lassen. Man definiert nun $\delta_0, \delta_1, \delta_2, \delta_3, \delta_4, \delta_5 \in \left\{0, 1\right\}$, um die zeitliche Diskretisierung der kinetischen Energie und der Massenflussdichte in der Form

\[ \begin{align} \frac{\mathbf{v}^2}{2} = \frac{\mathbf{v}_{\delta_0}\cdot\mathbf{v}_{\delta_1}}{2}, & {} & \rho\mathbf{v} = \frac{\rho_{\delta_2}\mathbf{v}_{\delta_3} + \rho_{\delta_4}\mathbf{v}_{\delta_5}}{2} \end{align} \]

zu notieren. Damit erhält man zunächst als Gleichungssystem

\[ \begin{align} \frac{\mathbf{v}_1 - \mathbf{v}_0}{\Delta t} = -\nabla\frac{\mathbf{v}^2}{2}, & {} & \frac{\rho_1 - \rho_0}{\Delta t} = -\nabla\cdot\left(\rho\mathbf{v}\right). \end{align} \]

Hieraus folgt

\[ \begin{align} &\frac{\rho_{\delta_2}\mathbf{v}_{\delta_3} + \rho_{\delta_4}\mathbf{v}_{\delta_5}}{2}\cdot\frac{\mathbf{v}_1 - \mathbf{v}_0}{\Delta t} = -\frac{\rho_{\delta_2}\mathbf{v}_{\delta_3} + \rho_{\delta_4}\mathbf{v}_{\delta_5}}{2}\cdot\nabla\frac{\mathbf{v}^2}{2} = -\frac{\rho_{\delta_2}\mathbf{v}_{\delta_3} + \rho_{\delta_4}\mathbf{v}_{\delta_5}}{2}\cdot\nabla \frac{\mathbf{v}_{\delta_0}\cdot\mathbf{v}_{\delta_1}}{2}\nonumber\\ &= -\nabla\cdot\left(\frac{\rho_{\delta_2}\mathbf{v}_{\delta_3} + \rho_{\delta_4}\mathbf{v}_{\delta_5}}{2}\frac{\mathbf{v}_{\delta_0}\cdot\mathbf{v}_{\delta_1}}{2}\right) + \frac{\mathbf{v}_{\delta_0}\cdot\mathbf{v}_{\delta_1}}{2}\nabla\cdot\frac{\rho_{\delta_2}\mathbf{v}_{\delta_3} + \rho_{\delta_4}\mathbf{v}_{\delta_5}}{2}\nonumber\\ &= -\nabla\cdot\left(\frac{\rho_{\delta_2}\mathbf{v}_{\delta_3} + \rho_{\delta_4}\mathbf{v}_{\delta_5}}{2}\frac{\mathbf{v}_{\delta_0}\cdot\mathbf{v}_{\delta_1}}{2}\right) - \frac{\mathbf{v}_{\delta_0}\cdot\mathbf{v}_{\delta_1}}{2}\frac{\rho_1 - \rho_0}{\Delta t}\nonumber\\ &\Leftrightarrow\frac{\rho_{\delta_2}\mathbf{v}_{\delta_3} + \rho_{\delta_4}\mathbf{v}_{\delta_5}}{2}\cdot\frac{\mathbf{v}_1 - \mathbf{v}_0}{\Delta t} + \frac{\mathbf{v}_{\delta_0}\cdot\mathbf{v}_{\delta_1}}{2}\frac{\rho_1 - \rho_0}{\Delta t} = -\nabla\cdot\left(\frac{\rho_{\delta_2}\mathbf{v}_{\delta_3} + \rho_{\delta_4}\mathbf{v}_{\delta_5}}{2}\frac{\mathbf{v}_{\delta_0}\cdot\mathbf{v}_{\delta_1}}{2}\right)\nonumber\\ &\Leftrightarrow \frac{\rho_{\delta_2}\mathbf{v}_{\delta_3}\cdot\mathbf{v}_1 + \rho_{\delta_4}\mathbf{v}_{\delta_5}\cdot\mathbf{v}_1 - \rho_{\delta_2}\mathbf{v}_{\delta_3}\cdot\mathbf{v}_0 - \rho_{\delta_4}\mathbf{v}_{\delta_5}\cdot\mathbf{v}_0 + \rho_1\mathbf{v}_{\delta_0}\cdot\mathbf{v}_{\delta_1} - \rho_0\mathbf{v}_{\delta_0}\cdot\mathbf{v}_{\delta_1}}{2\Delta t} = -\nabla\cdot\left(\frac{\rho_{\delta_2}\mathbf{v}_{\delta_3} + \rho_{\delta_4}\mathbf{v}_{\delta_5}}{2}\frac{\mathbf{v}_{\delta_0}\cdot\mathbf{v}_{\delta_1}}{2}\right)\nonumber\\ &\Leftrightarrow \frac{\rho_{\delta_2}\mathbf{v}_{\delta_3}\cdot\mathbf{v}_1 + \rho_{\delta_4}\mathbf{v}_{\delta_5}\cdot\mathbf{v}_1 + \rho_1\mathbf{v}_{\delta_0}\cdot\mathbf{v}_{\delta_1} - \rho_{\delta_2}\mathbf{v}_{\delta_3}\cdot\mathbf{v}_0 - \rho_{\delta_4}\mathbf{v}_{\delta_5}\cdot\mathbf{v}_0 - \rho_0\mathbf{v}_{\delta_0}\cdot\mathbf{v}_{\delta_1}}{2\Delta t} = -\nabla\cdot\left(\frac{\rho_{\delta_2}\mathbf{v}_{\delta_3} + \rho_{\delta_4}\mathbf{v}_{\delta_5}}{2}\frac{\mathbf{v}_{\delta_0}\cdot\mathbf{v}_{\delta_1}}{2}\right) \end{align} \]

Es muss gelten

\[ \begin{align} \rho_{\delta_2}\mathbf{v}_{\delta_3}\cdot\mathbf{v}_1 + \rho_{\delta_4}\mathbf{v}_{\delta_5}\cdot\mathbf{v}_1 + \rho_1\mathbf{v}_{\delta_0}\cdot\mathbf{v}_{\delta_1} - \rho_{\delta_2}\mathbf{v}_{\delta_3}\cdot\mathbf{v}_0 - \rho_{\delta_4}\mathbf{v}_{\delta_5}\cdot\mathbf{v}_0 - \rho_0\mathbf{v}_{\delta_0}\cdot\mathbf{v}_{\delta_1} = \rho_{1}\mathbf{v}_{1}\cdot\mathbf{v}_1 - \rho_0\mathbf{v}_{0}\cdot\mathbf{v}_{0}. \end{align} \]

Setzt man $\delta_0 = \delta_1 = 0$ ein, folgt

\[ \begin{align} \rho_{\delta_2}\mathbf{v}_{\delta_3}\cdot\mathbf{v}_1 + \rho_{\delta_4}\mathbf{v}_{\delta_5}\cdot\mathbf{v}_1 + \rho_1\mathbf{v}_{0}\cdot\mathbf{v}_{0} - \rho_{\delta_2}\mathbf{v}_{\delta_3}\cdot\mathbf{v}_0 - \rho_{\delta_4}\mathbf{v}_{\delta_5}\cdot\mathbf{v}_0 - \rho_0\mathbf{v}_{0}\cdot\mathbf{v}_{0} = \rho_{1}\mathbf{v}_{1}\cdot\mathbf{v}_1 - \rho_0\mathbf{v}_{0}\cdot\mathbf{v}_{0}. \end{align} \]

Mit $\delta_2 = \delta_3 = 1$ folgt

\[ \begin{align} \rho_{1}\mathbf{v}_{1}\cdot\mathbf{v}_1 + \rho_{\delta_4}\mathbf{v}_{\delta_5}\cdot\mathbf{v}_1 + \rho_1\mathbf{v}_{0}\cdot\mathbf{v}_{0} - \rho_{1}\mathbf{v}_{1}\cdot\mathbf{v}_0 - \rho_{\delta_4}\mathbf{v}_{\delta_5}\cdot\mathbf{v}_0 - \rho_0\mathbf{v}_{0}\cdot\mathbf{v}_{0} = \rho_{1}\mathbf{v}_{1}\cdot\mathbf{v}_1 - \rho_0\mathbf{v}_{0}\cdot\mathbf{v}_{0}. \end{align} \]

Hieraus wiederum folgt $\delta_4 = 1$, $\delta_5 = 0$. Die Paare $\left(\delta_2, \delta_3\right)$, $\left(\delta_4, \delta_5\right)$ sind vertauschbar, was aber keine neue Möglichkeit darstellt, da die Massenflussdichte sich dabei nicht ändert.

Setzt man $\delta_0 = \delta_1 = 1$ ein, folgt

\[ \begin{align} \rho_{\delta_2}\mathbf{v}_{\delta_3}\cdot\mathbf{v}_1 + \rho_{\delta_4}\mathbf{v}_{\delta_5}\cdot\mathbf{v}_1 + \rho_1\mathbf{v}_{1}\cdot\mathbf{v}_{1} - \rho_{\delta_2}\mathbf{v}_{\delta_3}\cdot\mathbf{v}_0 - \rho_{\delta_4}\mathbf{v}_{\delta_5}\cdot\mathbf{v}_0 - \rho_0\mathbf{v}_{1}\cdot\mathbf{v}_{1} = \rho_{1}\mathbf{v}_{1}\cdot\mathbf{v}_1 - \rho_0\mathbf{v}_{0}\cdot\mathbf{v}_{0}. \end{align} \]

Mit $\delta_2 = 0$, $\delta_3 = 1$ folgt

\[ \begin{align} \rho_{0}\mathbf{v}_{1}\cdot\mathbf{v}_1 + \rho_{\delta_4}\mathbf{v}_{\delta_5}\cdot\mathbf{v}_1 + \rho_1\mathbf{v}_{1}\cdot\mathbf{v}_{1} - \rho_{0}\mathbf{v}_{1}\cdot\mathbf{v}_0 - \rho_{\delta_4}\mathbf{v}_{\delta_5}\cdot\mathbf{v}_0 - \rho_0\mathbf{v}_{1}\cdot\mathbf{v}_{1} = \rho_{1}\mathbf{v}_{1}\cdot\mathbf{v}_1 - \rho_0\mathbf{v}_{0}\cdot\mathbf{v}_{0}. \end{align} \]

Hieraus wiederum folgt $\delta_4 = 0$, $\delta_5 = 0$.

Setzt man $\delta_0 = 0$, $\delta_1 = 1$ ein, folgt

\[ \begin{align} \rho_{\delta_2}\mathbf{v}_{\delta_3}\cdot\mathbf{v}_1 + \rho_{\delta_4}\mathbf{v}_{\delta_5}\cdot\mathbf{v}_1 + \rho_1\mathbf{v}_{0}\cdot\mathbf{v}_{1} - \rho_{\delta_2}\mathbf{v}_{\delta_3}\cdot\mathbf{v}_0 - \rho_{\delta_4}\mathbf{v}_{\delta_5}\cdot\mathbf{v}_0 - \rho_0\mathbf{v}_{0}\cdot\mathbf{v}_{1} = \rho_{1}\mathbf{v}_{1}\cdot\mathbf{v}_1 - \rho_0\mathbf{v}_{0}\cdot\mathbf{v}_{0}. \end{align} \]

Mit $\delta_2 = \delta_3 = 0$ folgt

\[ \begin{align} \rho_{0}\mathbf{v}_{0}\cdot\mathbf{v}_1 + \rho_{\delta_4}\mathbf{v}_{\delta_5}\cdot\mathbf{v}_1 + \rho_1\mathbf{v}_{0}\cdot\mathbf{v}_{1} - \rho_{0}\mathbf{v}_{0}\cdot\mathbf{v}_0 - \rho_{\delta_4}\mathbf{v}_{\delta_5}\cdot\mathbf{v}_0 - \rho_0\mathbf{v}_{0}\cdot\mathbf{v}_{1} = \rho_{1}\mathbf{v}_{1}\cdot\mathbf{v}_1 - \rho_0\mathbf{v}_{0}\cdot\mathbf{v}_{0}. \end{align} \]

Hieraus wiederum folgt $\delta_4 = 1$, $\delta_5 = 1$.

Zusammenfassend ergeben sich folgende drei Möglichkeiten:

\[ \begin{align} \frac{\mathbf{v}^2}{2} = \frac{\mathbf{v}_0\cdot\mathbf{v}_0}{2}, & \rho\mathbf{v} = \frac{\rho_1\mathbf{v}_0 + \rho_1\mathbf{v}_1}{2}\\ \frac{\mathbf{v}^2}{2} = \frac{\mathbf{v}_1\cdot\mathbf{v}_1}{2}, & \rho\mathbf{v} = \frac{\rho_0\mathbf{v}_0 + \rho_0\mathbf{v}_1}{2}\\ \frac{\mathbf{v}^2}{2} = \frac{\mathbf{v}_0\cdot\mathbf{v}_1}{2}, & \rho\mathbf{v} = \frac{\rho_0\mathbf{v}_0 + \rho_1\mathbf{v}_1}{2} \end{align} \]

25.5.1.2 Generalisierte Coriolis-Terme

Für die Zeitentwicklung der kinetischen Energie gilt

\[ \begin{align} \frac{\mathbf{v}_1^2 - \mathbf{v}_0^2}{2\Delta t} &= \frac{\left(\mathbf{v}_1 + \mathbf{v}_0\right)\cdot\left(\mathbf{v}_1 - \mathbf{v}_0\right)}{2\Delta t} = \mathbf{v}_{1/2}\cdot\frac{\mathbf{v}_1 - \mathbf{v}_0}{\Delta t} \end{align} \]

mit

\[ \begin{align} \mathbf{v}_{1/2} \coloneqq \frac{\mathbf{v}_0 + \mathbf{v}_1}{2}. \end{align} \]

Soll der verallgemeinerte Coriolisterm nicht zur Evolution der kinetischen Energie beitragen, muss dieser also mit dem Geschwindigkeitsfeld $\mathbf{v}_{1/2}$ berechnet werden. Die Überlegungen in diesem Abschnitt (25.5.1) wurden erstmals von Gassmann und Herzog in [23] vorgebracht.

25.5.1.3 Zeitentwicklung der inneren Energie

Für die Zeitentwicklung der inneren Energiedichte $\newtilde{I}$ kann man notieren

\[ \begin{align} \frac{\Delta\newtilde{I}}{\Delta t} &= \frac{\newtilde{I}_1 - \newtilde{I}_0}{\Delta t} = c^{(v)}\rho_0\frac{T_1 - T_0}{\Delta t} + c^{(v)}T_1\frac{\rho_1 - \rho_0}{\Delta t}\tag{25.64}\label{eq:deriv_temp_evol_model_1}\\ &= \frac{\newtilde{I}_1 - \newtilde{I}_0}{\Delta t} = c^{(v)}\rho_1\frac{T_1 - T_0}{\Delta t} + c^{(v)}T_0\frac{\rho_1 - \rho_0}{\Delta t}.\tag{25.65}\label{eq:deriv_temp_evol_model_2} \end{align} \]

Zunächst geht man von Glg. (25.64) aus. Mittels der Kontinuitätsgleichung erhält man

\[ \begin{align} \md{\rho} = -\rho_0\nabla\cdot\mathbf{v}. \end{align} \]

Aus Glg. (9.14) folgt

\[ \begin{align} \frac{T_1 - T_0}{\Delta t} = -\frac{R_d}{c^{(v)}}T_0\nabla\cdot\mathbf{v}. \end{align} \]

Beide Male wurde dabei davon ausgegangen, dass die skalaren Variablen vom alten Zeitschritt verwendet werden, wie es z. B. bei einem Forward-backward-Verfahren der Fall ist. Setzt man diese beiden Feststellungen in Glg. (25.64) ein, erhält man

\[ \begin{align} \frac{\newtilde{I}_1 - \newtilde{I}_0}{\Delta t} &= -c^{(v)}\rho_0\frac{R_d}{c^{(v)}}T_0\nabla\cdot\mathbf{v} - c^{(v)}T_1\rho_0\nabla\cdot\mathbf{v} = -c^{(v)}\rho_0\nabla\cdot\mathbf{v}\left(\frac{R_d}{c^{(v)}}T_0 + T_1\right)\nonumber\\ &= -c^{(p)}\rho_0\nabla\cdot\mathbf{v}\left(\frac{R_d}{c^{(p)}}T_0 + \frac{c^{(v)}}{c^{(p)}}T_1\right). \end{align} \]

Hieraus folgert man, dass die Temperatur in der Gewichtung

\[ \begin{align} T^\star \coloneqq \frac{R_d}{c^{(p)}}T_0 + \frac{c^{(v)}}{c^{(p)}}T_1\tag{25.69}\label{eq:c_v_c_p} \end{align} \]

in die diskretisierte Evolution der inneren Energiedichte eingeht.

Geht man jedoch von Glg. (25.65) aus, führt dies mit

\[ \begin{align} \md{\rho} = -\rho_1\nabla\cdot\mathbf{v} \end{align} \]

und

\[ \begin{align} \frac{T_1 - T_0}{\Delta t} = -\frac{R_d}{c^{(v)}}T_1\nabla\cdot\mathbf{v} \end{align} \]

auf

\[ \begin{align} \frac{\newtilde{I}_1 - \newtilde{I}_0}{\Delta t} &= -c^{(p)}\rho_1\nabla\cdot\mathbf{v}\left(\frac{c^{(v)}}{c^{(p)}}T_0 + \frac{R_d}{c^{(p)}}T_1\right),\\ T^\star &\coloneqq \frac{c^{(v)}}{c^{(p)}}T_0 + \frac{R_d}{c^{(p)}}T_1.\tag{25.73}\label{eq:r_d_c_p} \end{align} \]

Dies sind die beiden Möglichkeiten für energetisch selbstkonsistente zeitliche Diskretisierungen der inneren Energie.