\section{Einführung}

\subsection{Warum Gitterverfeinerung?}

Die einfachsten und zugleich am weitesten verbreiteten Umsetzungen von Simulationen mit Lattice Boltzmann Methoden basieren auf uniformen Gittern, in denen Zellen immer den gleichen Abstand zu ihren Nachbarzellen haben.

Die Genauigkeit von Lattice Boltzmann basierenden Simulationen hängt maßgeblich von der Auflösung des verwendeten Gitters ab. Bei Außerachtlassung weiterer wichtiger Faktoren wie dem verwendeten Kollisionsterm und Randkonditionen kann im Allgemeinen davon ausgegangen werden, dass eine feinere Auflösung des Gitters zu besseren Ergebnissen führt.

In praktischen Beispielen können innerhalb eines Modells große Unterschiede in der Strömungskomplexität existieren. So kann es große Gebiete eines Modells geben, die mit einem vergleichsweise groben Gitter gut simuliert werden können, während in anderen Gebieten \textendash{} beispielsweise in komplexen Geometrien und an Rändern \textendash{} ein vielfach feineres Gitter zur adäquaten Behandlung benötigt wird. In uniformen Gittern muss jedoch das gesamte Modell unabhängig der lokalen Situation mit der maximal benötigten Auflösung abgebildet werden.

Da die Anzahl der benötigten Gitterpunkte sich maßgeblich auf den Speicherbedarf und Rechenaufwand auswirkt, ist es wünschenswert die Anzahl der Gitterpunkte zu minimieren. Ein Ansatz, dies zu erreichen, ist die lokale Variation der Gitterauflösung.

\newpage
\section{Grundlagen}

In diesem Kapitel werden wir die, dem weiteren Verlauf dieser Arbeit zugrunde liegende, Lattice Boltzmann Methode in 2D nachvollziehen.

\subsection{Lattice Boltzmann Methode}

Grundlage und Namensgeber von Simulationen mit Lattice Boltzmann Methoden ist die Boltzmann Gleichung. Sie beschreibt das Verhalten von Gasen auf mesoskopischer Ebene als Verteilungsfunktion der Masse von Partikeln in einer Raumregion mit gegebener Geschwindigkeit.

\begin{Definition}[Die Boltzmann-Gleichung]
Sei \(f(x,\xi,t)\) die Verteilungsfunktion der Partikelmasse zu Zeit \(t\) in Ort \(x \in \R^2\) mit Geschwindigkeit \(\xi \in \R^2\), \(\rho\) die Dichte und \(F \in \R^2\) eine etwaige äußere Kraft. Die Boltzmann-Gleichung beschreibt die zeitliche Veränderung der Verteilungsfunktion anhand des totalen Differential \(\Omega(f)\):
\[ \left( \partial_t + \xi \cdot \partial_x + \frac{F}{\rho} \cdot \partial_\xi \right) f = \Omega(f) \left( = \partial_x f \cdot \frac{dx}{dt} + \partial_\xi f \cdot \frac{d\xi}{dt} + \partial_t f \right) \]

Hierbei handelt es sich um eine Advektionsgleichung wobei der Term \(\partial_t f + \xi \cdot \partial_x f\) die Strömung der Partikelverteilung mit Geschwindigkeit \(\xi\) und \(\frac{F}{\rho} \cdot \partial_\xi f\) einwirkende Kräfte darstellt. Der Term \(\Omega(f)\) beschreibt, entsprechend als Kollisionsoperator bezeichnet, die kollisionsbedingte lokale Neuverteilung von \(f\).
\end{Definition}

Zentrale Anforderung an den Kollisionsoperator ist die Impuls- und Masseerhaltung.

Die im Folgenden betrachtete Lattice Boltzmann Methode verwendet die übliche BGK Approximation der Boltzmann-Gleichung ohne äußere Kraft von Bhatnagar, Gross und Krook (siehe \citetitle{krueger17}~\cite[Kap.~3.5.3]{krueger17}).

Grundlegendes Element dieser Approximation ist der BGK Operator
\[\Omega(f) := -\frac{f-f^\text{eq}}{\tau} \Delta t\]
welcher die Partikelverteilung mit Rate \(\tau\) gegen eine Equilibriumsverteilung \(f^\text{eq}\) relaxiert. Ohne Beschränkung der Allgemeinheit setzen wir dabei  im Folgenden \(\Delta t = 1\).

Wenden wir den BGK Operator auf die Boltzmann-Gleichung ohne äußere Kräfte an, erhalten wir die BGK Approximation:

\begin{Definition}[BGK Approximation]
Sei \(\tau \in \R_{\geq 0}\) eine Relaxionszeit und \(f^\text{eq}\) die von der Maxwell-Boltzmann Verteilung gegebene Equilibriumsverteilung.
\[ (\partial_t + \xi \cdot \nabla_x) f = -\frac{1}{\tau} (f(x,\xi,t) - f^\text{eq}(x,\xi,t)) \]
\end{Definition}

An dieser Stelle bemerken wir, dass die BGK Approximation für beliebige \(\xi \in \R^2\) definiert ist. Da wir die LBM auf einem endlichen Rechner umsetzen wollen, müssen wir die Menge der betrachteten Geschwindigkeiten auf eine endliche Menge diskretisieren.

Eine übliche Menge diskreter Geschwindigkeiten in 2D ist \emph{D2Q9} wobei \emph{D2} die Anzahl der Dimensionen und \emph{Q9} die Anzahl der Geschwindigkeiten verschlüsselt.

\begin{Definition}[D2Q9 Modell]
\[ \{\xi_i\}_{i=0}^8 = \left\{ \V{0}{0}, \V{-1}{\phantom{-}1}, \V{-1}{\phantom{-}0}, \V{-1}{-1}, \V{\phantom{-}0}{-1}, \V{\phantom{-}1}{-1}, \V{1}{0}, \V{1}{1}, \V{0}{1} \right\} \]
\end{Definition}

\begin{figure}[t]
\centering
\input{img/d2q9.tikz}
\caption{Umgebung einer Zelle in D2Q9}
\end{figure}

Mithilfe einer solchen endlichen Menge diskreter Geschwindigkeiten lässt sich die BGK Approximation bezüglich der Geschwindigkeit diskretisieren:

\begin{Definition}[BGK Geschwindigkeitsdiskretisierung]
\label{def:disVelBGK}
Seien \(\xi_i\) Vektoren einer Menge mikroskopischer Geschwindigkeiten wie z.B. D2Q9 und \(f_i(x,t) \equiv f(x,\xi_i,t)\). Dann ist
\[ (\partial_t + \xi_i \cdot \nabla_x) f_i(x,t) = -\frac{1}{\tau} (f_i(x,t) - f_i^\text{eq}(x,t)) \]
die Diskretisierung der BGK Approximation entlang der Geschwindigkeiten.
\end{Definition}

Hierbei ist die diskrete Equilibriumsverteilung \(f_i^\text{eq}\) wie folgt definiert:

\begin{Definition}[Diskrete Equilibriumsverteilung]
\label{def:fieq}
Seien \(\rho \in \R_{\geq 0}\) die Dichte, \(u \in \R^2\) die Gesamtgeschwindigkeit, \(\xi_i\) die \(i\)-te diskrete Geschwindigkeitskomponente, \(w_i\) das Gewicht jener Komponente bzgl. des Lattice und \(c_s\) die Lattice-Schallgeschwindigkeit.
\[f_i^\text{eq} = w_i \rho \left( 1 + \frac{u \cdot \xi_i}{c_s^2} + \frac{(u \cdot \xi_i)^2}{2c_s^4} - \frac{u \cdot u}{2c_s^2} \right)\]
\end{Definition}

Die Werte von \(u = u(x,t)\) und \(\rho = \rho(x,t)\) in Ort \(x\) zu Zeit \(t\) ergeben sich dabei aus den \emph{Momenten} der Verteilungsfunktion \(f_i\):

\begin{Definition}[Momente der Verteilungsfunktion]
\label{def:Momente}
\[\rho(x,t) = \sum_{i=0}^{q-1} f_i(x,t) \text{ und } \rho u(x,t) = \sum_{i=0}^{q-1} \xi_i f_i(x,t)\]
\end{Definition}

Für D2Q9 ergeben sich nach \cite[Gl.~3.60 bzw. Tab.~3.3]{krueger17} die Gewichte:
\[w_0 = \frac{4}{9}, \ w_{2,4,6,8} = \frac{1}{9}, \ w_{1,3,5,7} = \frac{1}{36}\]

Weiter folgt zusammen mit der Bedingung \(\sum_{i=1}^{q-1} w_i (\xi_i)_a (\xi_i)_b = c_s^2 \delta_{a,b}\) aus \cite[Gl.~3.60]{krueger17} die Schallgeschwindigkeit \(c_s = \sqrt{1/3}\) des Lattice. Konditionen zur Bestimmung dieser gitterspezifischen Konstanten sind hierbei die Erhaltung von Impuls und Masse sowie die Forderung von \emph{Rotationsisotropie}.

Zur Entwicklung einer \emph{implementierbaren} expliziten BGK Gleichung können wir nun die Geschwindigkeitsdiskretisierung~\ref{def:disVelBGK} integrieren:
\[ f_i(x+\xi_i, t+1) - f_i(x,t) = \int_0^1 \Omega_i(x+\xi_i s,t+s) ds \]

Wobei \(\Omega_i(x,t)\) hier die diskrete Formulierung des BGK Kollisionsoperators darstellt:
\[\Omega_i(x,t) := -\frac{1}{\tau} ( f_i(x,t) - f_i^\text{eq}(x, t) )\]

Da sich die exakte Lösung des Integrals auf der rechten Seite schwierig gestaltet, wird es in der Praxis nur approximiert. Während es dazu vielfältige Ansätze gibt, beschränken wir uns an dieser Stelle auf Anwendung der Trapezregel:
\begin{align*}
f_i(x+\xi_i,t+1) - f_i(x,t) &= \frac{1}{2} \left( \Omega_i(x,t) + \Omega_i(x+\xi_i,t+1) \right) \\
&= -\frac{1}{2\tau} \left( f_i(x+\xi_i,t+1) + f_i(x,t) - f_i^\text{eq}(x+\xi_i,t+1) - f_i^\text{eq}(x,t) \right)
\end{align*}

Zur expliziten Lösung dieser impliziten Gleichung benötigen wir nun nur noch eine geeignete Verschiebung von \(f_i\) und \(\tau\):

\begin{Definition}[Diskrete LBM BGK Gleichung]
\label{def:LBGKeq}
Seien \(\overline{f_i}\) und \(\overline\tau\) definiert:
\begin{align*}
\overline{f_i} &= f_i + \frac{1}{2\tau}(f_i - f_i^\text{eq}) \\
\overline\tau &= \tau + \frac{1}{2}
\end{align*}

Setzen wir diese verschobenen Variablen in das Ergebnis der Trapezregel ein, erhalten \cite[Kap.~A.5 mit \(\Delta t=1\)]{krueger17} wir die die vollständig diskretisierte LBM BGK Gleichung:

\[\overline{f_i}(x+\xi_i,t+1) = \overline{f_i}(x,t) - \frac{1}{\overline\tau} (\overline{f_i}(x,t) - f_i^\text{eq}(x,t))\]
\end{Definition}

Bemerkenswert ist an dieser Stelle, dass die Momente der Verteilungen mit \(\overline{f_i}\) analog zu \ref{def:Momente} berechnet werden können:
\begin{align*}
\sum_{i=0}^{q-1} \overline{f_i} = \sum_{i=0}^{q-1} f_i + \frac{1}{2\tau} \sum_{i=0}^{q-1} (f_i - f_i^\text{eq}) &= \rho \\
\sum_{i=0}^{q-1} \xi_i \overline{f_i} = \sum_{i=0}^{q-1} \xi_i f_i + \frac{1}{2\tau} \sum_{i=0}^{q-1} (f_i - f_i^\text{eq}) &= \rho u
\end{align*}

\newpage
\subsubsection{Implementierung}

Bei der Implementierung der {diskreten LBM BGK Gleichung}~\ref{def:LBGKeq} auf einem Computer ist die Aufteilung in Kollisions- und Strömungsschritt üblich.

\begin{Definition}[Kollisionsschritt]
Annäherung der Verteilungsfunktion an die lokal berechnete Equilibriumsverteilung entsprechend dem BGK Kollisionsoperator.
\[f_i^\text{out}(x,t) = f_i(x,t) - \frac{1}{\tau}(f_i(x,t) - f_i^\text{eq}(x,t))\]
\end{Definition}

\begin{Definition}[Strömungsschritt]
Strömen der neuen Verteilungen auf die benachbarten Zellen entsprechend der jeweiligen diskreten Geschwindigkeit.
\[f_i(x+\xi_i,t+1) = f_i^\text{out}(x,t)\]
\end{Definition}

Bemerkenswert ist hierbei, dass der Kollisionsschritt nur lokale Informationen der jeweiligen Zelle benötigt und sich somit sehr gut zur parallelen Verarbeitung eignet.

\subsubsection{Chapman-Enskog Analyse}

Ziel der beschriebenen Lattice Boltzmann Methode ist die möglichst gute Approximation der schwach-kompressiblen Navier-Stokes Gleichungen auf der Simulationsdomäne.

\begin{Definition}[Schwach-kompressible Navier-Stokes Gleichungen]
Sei \(\rho\) die Dichte, \(u\) die Geschwindigkeit und \(p\) der Druck zu Zeit \(t\) sowie \(\nu\) die kinematische Viskosität und \(\mathrm{S}\) der Verzerrungstensor.
\begin{align*}
\partial_t  \rho + \nabla \cdot (\rho u) &= 0 \\
\partial_t u + (u \cdot \nabla) u &= -\frac{1}{\rho} \nabla p + 2\nu\nabla \cdot (\mathrm{S})
\end{align*}

Dabei sind Druck \(p\), kinetische Viskosität \(\nu\) und Verzerrungstensor \(\mathrm{S}\) definiert als:
\begin{align*}
p &= c_s^2 \rho \\
\nu &= c_s^2 \tau \\
\mathrm{S} &= \frac{1}{2} (\nabla u + (\nabla u)^\top) \numberthis\label{eq:Verzerrungstensor}
\end{align*}
\end{Definition}

Nach \cite[Kap.~4.1]{krueger17} kann die asymptotische Äquivalenz von LBM BGK Gleichung und schwach-kompressiblen Navier-Stokes Gleichungen mit der Entwicklung von Chapman-Enskog gezeigt werden.

\begin{Definition}[Chapman-Enskog Ansatz]
\label{def:ChapmanEnskog}
Der Chapman-Enskog Ansatz besteht in der Annahme, dass die Verteilungsfunktion \(f_i\) als leicht gestörte Equilibriumsverteilung dargestellt werden kann: \[f_i = f_i^\text{eq} + \epsilon f_i^{(1)} + \mathcal{O}(\epsilon^2)\]
Hierbei ist \(\epsilon f_i^{(1)}\) mit \(\epsilon \ll 1\) der Störterm 1. Ordnung. Dieser ist gegeben als:
\[\epsilon f_i^{(1)} = \frac{w_i}{2 c_s^4} \mathrm{Q}_i : \mathrm{\Pi}^{(1)} \numberthis\label{eq:firstOrderPertubation}\]
Wobei der Geschwindigkeitstensor \(\mathrm{Q}_i\) und das Störmoment \(\mathrm\Pi^{(1)}\) nach \cite[Kap.~4.1.3]{krueger17} definiert sind als:
\begin{align*}
\mathrm{Q}_i &= \xi_i \xi_i - c_s^2 I \numberthis\label{eq:velocityTensor} \\
\mathrm\Pi^{(1)} &= \sum_{i=0}^{q-1} \xi_i \xi_i \epsilon f_i^{(1)} = -2 c_s^2 \rho \overline\tau \mathrm{S} \numberthis\label{eq:pertubationMoment}
\end{align*}
\end{Definition}

\begin{Definition}[Nicht-Equilibriumsverteilung]
Die Relaxion der Verteilungsfunktion \(f_i\) gegen die Equilibriumsverteilung \(f_i^\text{eq}\) impliziert eine Dekomposition in Equilibriums- und Nicht-Equilibriumsverteilung: \[f_i := f_i^\text{eq} + f_i^\text{neq}\]
\end{Definition}

Unter Vernachlässigung von Störtermen der Ordnung \(\mathcal{O}(\epsilon^2)\) ergibt sich eine Näherung der Nicht-Equilibriumsverteilung:
\[f_i^\text{neq} \cong \epsilon f_i^{(1)}\]

Diese Darstellung können wir unter Verwendung von Definition~\ref{def:ChapmanEnskog} ausführen als:
\[f_i^\text{neq} \cong -\frac{w_i \rho \overline\tau}{c_s^2} \mathrm{Q}_i : \mathrm{S} \numberthis\label{eq:approxFneq}\]

\newpage
\subsection{Herangehensweisen an Gitterverfeinerung}

Grundsätzlich existieren mit der Multi-Grid bzw. Multi-Domain Herangehenweise zwei verschiedene Ansätze für Gitterverfeinerung in Lattice Boltzmann Methoden \cite[Kap.~3.1]{lagrava12}. Im Westenlichen unterscheiden die Ansätze sich in der Struktur der variabel aufgelösten Teilgitter der Simulationsdomäne. Weitere Unterschiede folgen dann aus dieser grundlegenden Struktur.

\subsubsection{Multi-Grid Ansatz}

Bei Verfahren des Multi-Grid Ansatzes existiert das gröbste Gitter über der gesamten Domäne. Feinere Teilgitter werden über gröberen Gittern platziert und nicht aus deren Verarbeitung ausgeschlossen. Somit existieren über der gesamten \emph{Fläche} feinerer Gitter auch die Knoten gröberer Gitter.

\begin{figure}[h]
\centering
\input{img/multi_grid.tikz}
\caption{Teilgitter in der Multi-Grid Herangehensweise}
\end{figure}

Vorteil dieser Herangehensweise ist es, dass feinere Teilgitter im Verlauf der Simulation ohne aufwendige Restrukturierung verschoben werden können, um beispielsweise komplexere Strömungsstrukturen zu \emph{verfolgen}. Nachteil ist, dass das Einsparpotentiale in Speicher- und Rechenbedarf durch Mehrfachabdeckung von Teilen der Simulationsdomäne mittels mehrerer Gitter nicht voll ausgenutzt werden können.

\subsubsection{Multi-Domain Ansatz}

Kern des Multi-Domain Ansatzes ist es, außerhalb von etwaigen verfahrensbedingten Übergangsbereichen, jede Position in der Simulationsdomäne durch genau ein Teilgitter abzubilden. Konkret werden also bereits durch feinere Gitter abgedeckte Bereiche aus gröberen Teilgittern ausgespart.

\begin{figure}[h]
\centering
\input{img/multi_domain_plain.tikz}
\caption{Teiligitter in der Multi-Domain Herangehensweise}
\end{figure}

Vorteil gegenüber des Multi-Grid Ansatzes ist hier der weiter reduzierte Speicherbedarf sowie erwartete Einsparungen in der benötigten Rechenzeit. Erkauft werden diese Vorteile durch aufwendigere Kopplung der verschiedenen Teilgitter in den Übergangsbereichen.

\newpage
\section{Verfeinerungsmethode}

Wie in Kapitel~\ref{sec:olbRefinementChoice} noch näher begründet werden wird, bieten sich der Multi-Domain Ansatz als Grundlage für Gitterverfeinerung in OpenLB an. Passend zu dieser Wahl sowie der, im Rahmen dieser Arbeit getroffenen, Einschränkung auf zweidimensionale LBM mit BGK-Kollisionsoperator haben Lagrava et al. 2012 in \citetitle{lagrava12}~\cite{lagrava12} ein solches Verfeinerungsverfahren entwickelt. Die Stuktur dieses Verfahrens, mit potenziell austauschbaren Restriktions- und Interpolationsverfahren im zentralen Kopplungsschritt, erscheint dabei sogleich auch als Fundament eines Multi-Domain Gitterverfeinerungsframeworks in OpenLB.

\subsection{Übersicht}

Das Verfahren basiert auf der Multi-Domain Herangehensweise \cite[Kap.~3.1]{lagrava12} , bei welcher die feiner aufgelösten Teilbereiche der Simulationsdomäne aus dem gröber aufgelösten Gitter ausgeschlossen werden und sich nur in Übergangsbereichen überlappen.

\begin{figure}[h]
\centering
\input{img/multi_domain.tikz}
\caption{Multi-Domain Herangehensweise mit Übergangsbereich \cite[vgl.~Abb.~3]{lagrava12}}
\label{fig:MultiDomainOverlap}
\end{figure}

In diesen Übergangsbereichen, welche eine Breite von mindestens einer Einheit des gröberen Zellabstands haben, liegt die Hauptarbeit des Verfeinerungsverfahrens.

Während der Übergang vom feinen zum groben Gitter sich im Wesentlichen auf eine skalierte und gefilterte Restriktion der Verteilungen beschränkt, gestaltet sich der Übergang vom groben zum feinen Gitter aufwendiger, da feine Knoten, für deren Position kein grober Knoten existiert, aus den übrigen Daten interpoliert werden müssen.

Entsprechend liegt der Fokus des von Lagrava et al. entwickelten Algorithmus auf der Auswahl des Interpolationsverfahrens sowie der Skalierung der physikalischen Werte zwischen den unterschiedlich aufgelösten Verteilungen. Um diese Kopplung der verschiedenen Gitterauflösungen theoretisch erfassen zu können, müssen wir zunächst die Gitter selbst konkreter definieren:

\begin{Definition}[Diskretisierung der Gitter]
Wir betrachten zwei Gitter \(\G\) und \(\F\) als Diskretisierung der analytischen Simulationsdomänen \(D_g\) bzw. \(D_f\). Die Domänen seien so gewählt, dass sie gerade die konvexen Hüllen ihrer Diskretisierungsgitter darstellen.
\begin{align*}
\G &\subset D_g \cap \{ x \in \R^2 | \exists i \in \Z^2 : x = \delta x_g \cdot i \} && \text{Gröberes Gitter} \\
\F &\subset D_f \cap \{ x \in \R^2 | \exists i \in \Z^2 : x = \delta x_f \cdot i \} && \text{Feineres Gitter}
\end{align*}
\(\delta x_g = \delta x_g / 2 \in \R_+\) seien die Diskretisierungsauflösungen im Verhältnis \(1:2\).
\end{Definition}

Zur Betrachtung der Gitterkopplung fordern wir, dass sich \(\G\) und \(\F\) um eine grobe Gitterweite \(\delta x_g\) überlappen, vgl. Abbildungen \ref{fig:MultiDomainOverlap} und \ref{fig:OverlapZone}. Die Seitenlängen der konvexen Hüllen \(D_g\) und \(D_f\) sind ganzzahlige Vielfache von \(\delta x_g\) und \(\delta x_f\). Wir benötigen diese in Form der Gitter diskretisierten Mengen, um die Gitterknoten der Übergangsbereiche näher zu klassifizieren:

\begin{Definition}[Gitterknoten der Übergangsbereiche]
\label{def:OverlapGridNodes}
\begin{align*}
\U_g &:= \G \cap \F && \text{Grobe Knoten im Übergangsbereich} \\
\U_f &:= \F \cap (D_f \cap D_g) && \text{Feine Knoten im Übergangsbereich} \\
\U_{g \to f} &:= \partial D_f \cap \U_f && \text{Knoten mit Übertragung von grob nach fein} \\
\U_{f \to g} &:= \partial D_g \cap \U_g && \text{Knoten mit Übertragung von fein nach grob}
\end{align*}
\end{Definition}

Die Übertragungsrichtungen in \(\U_{g \to f}\) und \(\U_{f \to g}\) ergeben sich aus den jeweils fehlenden Verteilungsfunktionen an den Rändern der Gitter. So fehlen beispielsweise zur Kollision der groben Gitterknoten in \(\U_{f \to g}\) Verteilungsfunktionen in Richtung des feinen Gitters, während die feinen Zellen in dieser Menge noch vollständig definiert sind, da sie im Inneren des feinen Gitters liegen.

Mit diesem Argument lässt sich auch die Notwendigkeit eines Übergangsbereiches \(\U_g \cup \U_f\) der Mindestbreite \(\delta x_g\) erklären: Gäbe es diesen nicht, so fehlten an der Grenze zwischen grobem und feinem Gitter Verteilungsfunktionen in beide Richtungen zugleich.

\begin{figure}[h]
\centering
\input{img/overlap_zone.tikz}
\caption{Skizze des Übergangsbereich \cite[vgl.~Abb.~4]{lagrava12}}
\label{fig:OverlapZone}
\end{figure}

Die Aufgabe der Skalierungs-, Restriktions- und Interpolationsschritte des Verfahrens besteht also \emph{nur} darin, diese jeweils fehlenden Verteilungsfunktionen möglichst gut zu rekonstruieren.

\newpage
\subsection{Skalierung}

Während die Skalierung räumlicher Größen durch die Festlegung des Verfahrens auf Übergänge im Verhältnis \(1:2\) definiert ist, eröffnen sich für die zeitliche Skalierung zwei Möglichkeiten: Konvektive oder diffusive Skalierung. Unterschied der beiden Ansätze ist dabei das jeweilige Verhältnis zwischen räumlicher und zeitlicher Auflösung.

\begin{Definition}[Konvektive Skalierung]
Sei \(\delta t > 0\) die zeitliche und \(\delta x > 0\) die räumliche Diskretisierung. Dann gilt bei konvektiver Skalierung das Verhältnis:
\[ \delta t \sim \delta x \]
Es besteht also eine lineare Proportionalität.
\end{Definition}

\begin{Definition}[Diffusive Skalierung]
Sei \(\delta t > 0\) die zeitliche und \(\delta x > 0\) die räumliche Diskretisierung. Dann gilt bei diffusiver Skalierung das Verhältnis:
\[ \delta t \sim \delta x^2 \]
Es besteht hier also eine quadratische Proportionalität. Im Vergleich zu einer konvektiven Skalierung ist die zeitliche Auflösung somit um eine Ordnung feiner.
\end{Definition}

Es ist klar zu erkennen, dass diffusive Skalierung einen deutlich größeren numerischen Aufwand gegenüber der konvektiven Skalierung nach sich zieht. Vorteil der bei diffusiver Skalierung erhöhten Schrittanzahl pro Zeiteinheit sind kleinere Fehler.

Für die Autoren des hier erörterten Gitterverfeinerungsverfahrens überwog jedoch das Argument der numerischen Effizienz, weshalb auch wir hier nun die konvektive Skalierung betrachten wollen. Die Austauschbarkeit des Skalierungsverfahrens sollte jedoch bei der Implementierung in OpenLB beachtet werden, da dieser Aspekt eine weitere prinzipiell flexible Komponente des Verfahrens darstellt.

\bigskip

Aus der Wahl von konvektiver Skalierung ergibt sich zunächst:
\[\frac{\delta t_g}{\delta x_g} = \frac{\delta t_f}{\delta x_f} \land \delta x_f = \frac{\delta x_g}{2} \implies \delta t_f = \frac{\delta t_g}{2} \numberthis\label{eq:gridTime}\]
Auf dem feinen Gitter müssen also doppelt so viele Zeitschritte wie auf dem groben Gitter durchgeführt werden. Geschwindigkeit, Druck und Dichte sind stetig im Gitterübergang. Dies gilt nicht für die kinetische Viskosität \(\nu = c_s^2 \tau\), was wir bei der Herleitung der feinen Relaxionszeit \(\tau_f\) aus \(\tau_g\) beachten müssen.

\begin{Definition}[Physikalische Reynolds-Zahl]
\label{def:PhysicalReynoldsNumber}
Seien \(U\) die charakteristische Geschwindigkeit, \(L\) die charakteristische Länge und \(\nu\) die kinetische Viskosität in physikalischen Einheiten. Dann ist die Reynolds-Zahl definiert als: \[\text{Re} := \frac{U L}{\nu}\]
\end{Definition}

\begin{Definition}[Lattice Reynolds-Zahl]
\label{def:LatticeReynoldsNumber}
Sei \(\# \in \{f, g\}\) Symbol des feinen oder groben Gitters.
Seien \(U_\# := \delta t_\# / \delta x_\# \cdot U\) die charakteristische Geschwindigkeit, \(L_\# := 1 / \delta x_\# \cdot L\) die charakteristische Länge und \(\nu_\# := c_s^2 \tau_\#\) die kinetische Viskosität in Lattice-Einheiten. Dann ist die \emph{Lattice} Reynolds-Zahl des feinen bzw. groben Gitters definiert als: \[ \text{Re}_\# := \frac{U_\# L_\#}{\nu_\#} = \frac{\delta t_\# U L}{(\delta x_\#)^2 \nu_\#} \]
\end{Definition}

Wir erzwingen nun mit \(\text{Re}_g = \text{Re}_f\) die Unabhängigkeit von Reynolds-Zahl und Gitterauflösung. Diese Gleichsetzung ist sinnvoll, da die Reynolds-Zahl gerade die Vergleichbarkeit von Strömungen verschiedener Modellgrößen ermöglicht. Durch Einsetzen von Definition~\ref{def:LatticeReynoldsNumber} erhalten wir eine Verknüpfung der Relaxionszeiten \(\tau_f\) und \(\tau_g\):
\begin{align*}
\text{Re}_g = \text{Re}_f &\iff \frac{\delta t_g U L}{(\delta x_g)^2 \nu_g} = \frac{\delta t_f U L}{(\delta x_f)^2 \nu_f} \\
&\iff \frac{\delta t_g}{(\delta x_g)^2 \nu_g} = \frac{\delta t_f}{(\delta x_f)^2 \nu_f} \\
&\iff \frac{\delta t_g}{(\delta x_g)^2 c_s^2 \tau_g} = \frac{\delta t_f}{(\delta x_f)^2 c_s^2 \tau_f} \\
&\iff \tau_f = \frac{\delta t_f (\delta x_g)^2}{(\delta x_f)^2 \delta t_g} \tau_g \\
&\iff \tau_f = 2 \tau_g \numberthis\label{eq:gridTau}\\
\end{align*}

Für die zur expliziten Lösung der diskreten LBM BGK Gleichung in Definition~\ref{def:LBGKeq} verschobenen Relaxionszeiten ergibt sich somit:
\[\overline{\tau_f} = 2 \overline{\tau_g} - \frac{1}{2} \numberthis\label{eq:gridTauShift}\]

Die Equilibriumsverteilung \(f_i^\text{eq}\) ergibt sich nach Definition~\ref{def:fieq} aus Geschwindigkeit \(u\) und Dichte \(\rho\). Sie sind also explizit unabhängig der Gitterauflösung und, wie erwähnt, stetig im Gitterübergang. Diese Aussage gilt nicht für die Nicht-Equilibriumsverteilung \(f_i^\text{neq}\), da diese nach (\ref{eq:approxFneq}) von dem Geschwindigkeitsgradienten \(\nabla u\) abhängt.

Bezeichnen nun \(f_{f,i}^\text{neq}\) und \(f_{g,i}^\text{neq}\) die gitterspezifischen Nicht-Equilibriumsanteile und \(\mathrm{S}_f\) sowie \(\mathrm{S}_g\) die entsprechenden Verzerrungstensoren. Zur Skalierung von \(f_{f,i}^\text{neq}\) suchen wir ein \(\alpha \in \R\) s.d. gilt: \[f_{f,i}^\text{neq} = \alpha f_{g,i}^\text{neq} \numberthis\label{eq:scaleFneqReq}\]

Mit Hilfe von \ref{eq:approxFneq} lässt sich diese Relation nun nach \(\alpha\) auflösen:
\begin{align*}
f_{f,i}^\text{neq} = \alpha f_{g,i}^\text{neq} &\iff -\frac{w_i \rho \overline{\tau_f}}{c_s^2} \mathrm{Q}_i : \mathrm{S}_f = -\alpha \left( \frac{w_i \rho \overline{\tau_g}}{c_s^2} \mathrm{Q}_i : \mathrm{S}_g \right) \\
&\iff \overline{\tau_f} \mathrm{Q}_i : \mathrm{S}_f = \alpha \overline{\tau_g} \mathrm{Q}_i : \mathrm{S}_g \\
&\iff \frac{\overline{\tau_f}}{\delta t_f} \mathrm{Q}_i : \mathrm{S} = \alpha \frac{\overline{\tau_g}}{\delta t_g} \mathrm{Q}_i : \mathrm{S} \\
&\iff \alpha = \frac{\delta t_g}{\delta t_f} \frac{\overline{\tau_f}}{\overline{\tau_g}}\\
\end{align*}

Auflösen dieses \(\alpha\) in (\ref{eq:scaleFneqReq}) und Einsetzen der Relationen (\ref{eq:gridTime}) sowie (\ref{eq:gridTauShift}) ergibt dann:
\begin{align*}
f_{f,i}^\text{neq} &= \frac{\delta t_g}{\delta t_f} \frac{\overline{\tau_f}}{\overline{\tau_g}} f_{g,i}^\text{neq} \\
&= 2 \frac{2\overline{\tau_g} - \frac{1}{2}}{\overline{\tau_g}} f_{g,i}^\text{neq} \\
&= \left( 4 - \frac{1}{\overline{\tau_g}} \right) f_{g,i}^\text{neq} \numberthis\label{eq:scaleFneq}
\end{align*}

Insgesamt haben wir hiermit die Skalierung der Diskretisierungen in Raum und Zeit, der Relaxionszeit sowie der Nicht-Equilibriumsverteilung zwischen den Gittern \(\F\) und \(\G\) hergeleitet.

\bigskip

Seien \(x_{f \to g} \in \U_{f \to g}\) und \(x_{g \to f} \in \U_{g \to f}\) die Knoten aus dem Übergangsbereich mit Übertragung von fein nach grob bzw. von grob nach fein. Dann gilt:
\begin{align}
f_{f,i}(x_{g \to f}) &= f_i^\text{eq}(\rho(x_{g \to f}), u(x_{g \to f})) + \left(4-\frac{1}{\overline{\tau_g}}\right) f_{g,i}^\text{neq}(x_{g \to f}) \label{eq:basicG2F}\\
f_{g,i}(x_{f \to g}) &= f_i^\text{eq}(\rho(x_{f \to g}), u(x_{f \to g})) + \left(4-\frac{1}{\overline{\tau_g}}\right)^{-1} f_{f,i}^\text{neq}(x_{f \to g}) \label{eq:basicF2G}
\end{align}

Die zusammengesetzten Verteilungsfunktionen von Übergangsknoten des einen Gitters lassen sich also durch Skalierung des Nicht-Equilibriumsanteils der Verteilungsfunktionen des jeweils anderen Gitters rekonstruieren. Leider reicht dies noch nicht zur vollständigen Beschreibung eines Gitterverfeinerungsverfahrens, da nicht für alle feinen Gitterknoten im Übergangsbereich passende grobe Gitterknoten existieren -- vgl. dazu Abbildung~\ref{fig:OverlapZone}. Auch der Übergang von fein nach grob gestaltet sich trotz passenden feinen Knoten potenziell komplizierter, als eine bloße Skalierung, wie wir im nächsten Kapitel zeigen werden.

\newpage
\subsection{Restriktion}

Kraft seiner höheren Auflösung enthält das feine Gitter mehr Informationen als das umgebende grobe Gitter. Der Übergang von fein nach grob stellt also eine Restriktion der Verteilungsfunktionen dar.

Konkret suchen wir nach einer sinnvollen Definition der in \(x_{f \to g} \in \U_{f \to g}\) fehlenden Verteilungsfunktionen \(f_{g,i}\). Eine Solche ergibt sich aus der skalierten Dekomposition \ref{eq:basicF2G} durch Ersetzen der einfachen Nicht-Equilibriumsverteilung \(f_{f,i}^\text{neq}(x_{f \to g})\) mit einer restringierten Variante ebendieser.
\[f_{g,i}(x_{f \to g}) = f_i^\text{eq}(\rho(x_{f \to g}), u(x_{f \to g})) + \alpha^{-1} \rfneq{i}{x_{f \to g}} \numberthis\label{eq:restrictedF2G}\]

Wir bemerken an dieser Stelle, dass nur die Nicht-Equilibriumsverteilung durch die Restriktionsoperation \(\boldsymbol{r}\) eingeschränkt wird, während der Equilibriumsanteil unangetastet bleibt. Dies ist damit zu begründen, dass Dichte und Geschwindigkeit bei der verwendeten konvektiven Skalierung im Gitterübergang stetig bleiben.

Die skalierte Dekomposition \ref{eq:basicF2G} lässt sich in der Schreibweise von \ref{eq:restrictedF2G} formulieren, wenn die Identität als Restriktionsoperation eingesetzt wird: \[\rfneq{i}{x_{f \to g}} := f_{f,i}^\text{neq}(x_{f \to g})\]

Die für unser Verfahren \cite[Kap.~3.3]{lagrava12} beschriebene Restriktion ist der Mittelwert aller umliegenden gerichteten Nicht-Equilibriumsanteilen: \[\rfneq{i}{x_{f \to g}} := \frac{1}{q} \sum_{j=0}^{q-1} f_{f,j}^\text{neq}(x_{f \to g} + \delta x_f \xi_j) \numberthis\label{eq:neqAvgRestrictionF2G}\]

\subsection{Interpolation}

Zunächst ergänzen wir die Gitterteilmengen aus Definition~\ref{def:OverlapGridNodes} um eine Unterscheidung zwischen Knoten, für die ein, der Übertragung von grob nach fein dienlicher, passender grober Knoten existiert und alleinstehenden feinen Knoten.

\begin{Definition}[Gitterknoten mit Übertragung von grob nach fein]
\begin{align*}
\U_{g \to f}^g &:= \U_{g \to f} \cap \U_g && \text{Grobe Knoten mit Übertragung von grob nach fein} \\
\U_{g \to f}^f &:= \U_{g \to f} \setminus \U_g && \text{Alleinstehende feine Knoten mit Übertragung von grob nach fein}
\end{align*}
\end{Definition}

\newpage
\section{Implementierung in OpenLB}

\subsection{Auswahl der Verfeinerungsmethode}\label{sec:olbRefinementChoice}

\subsection{Struktur des Gitterverfeinerungsframework}

\subsection{Umsetzung des Verfahrens von Lagrava et al.}

\newpage
\section{Evaluierung}

\subsection{Wahl der Beispiele}

\subsection{Rohrströmung}