diffusion-konvektion.tex

\documentclass[11pt, ngerman,colorback,accentcolor=tud2d]{tudreport}
\usepackage[utf8]{inputenc}
\usepackage[ngerman]{babel}
\usepackage{amsmath}

\usepackage{microtype}

\usepackage{pgfplots}
\pgfplotsset{compat=1.3}

\begin{document}
%\tableofcontents
\chapter{1D-Diffusionsgleichung mit FVM}

\begin{equation} \label{base-formula}
\alpha \frac{\partial^2 \phi}{\partial x^2} = f
\end{equation}

Bei Integralbildung über den Bereich eines Kontrollvolumens von $x_w$ bis $x_e$
und der Anwendung des Gaußschen Integralsatzes kann das wie folgt umgeschrieben werden.

\begin{align*}
\int_{x_w}^{x_e}\alpha \frac{\partial^2 \phi}{\partial x^2} dx
&= \alpha \frac{\partial \phi}{\partial x}\Big|_{x_w}^{x_e}\\
&= \alpha \frac{\partial \phi}{\partial x}\Big|_{x_e}
 - \alpha \frac{\partial \phi}{\partial x}\Big|_{x_w}\\
&= \alpha (\phi'_e - \phi'_w)
\end{align*}

Nun müssen die auftretenden Ableitungen durch Variablenwerte in den Mittelpunkten
der Kontrollvolumen ausgedrückt werden.

\begin{equation*}
\int_{x_w}^{x_e}f dx \approx f_P (x_e - x_w)
\end{equation*}

\begin{equation*}
\phi'_e \approx \frac{\phi_E - \phi_P}{x_E - x_P}
\end{equation*}

\begin{equation*}
\phi'_w \approx \frac{\phi_P - \phi_W}{x_P - x_W}
\end{equation*}

Eingesetzt in Gleichung~\ref{base-formula} ergibt dies:

\begin{equation}
\alpha \frac{\phi_E - \phi_P}{x_E - x_P} - \alpha \frac{\phi_P - \phi_W}{x_P - x_W}
= f_P (x_e - x_w)
\end{equation}

Diese Gleichung kann in folgende Beziehung umgeformt werden.

\newcommand{\diffEP}{(x_E-x_P)}
\newcommand{\diffPW}{(x_P-x_W)}
\newcommand{\diffew}{(x_e-x_w)}

\begin{equation}
a_P \phi_P + a_E \phi_E + a_W \phi_W = b_P
\end{equation}

Die Koeffizienten $a_i$ und $b_P$ nehmen hierbei die folgenden Werte an.

\begin{align*}
  a_P &= -\frac{\alpha}{\diffEP \diffew} - \frac{\alpha}{\diffPW \diffew}\\
a_E &= \frac{\alpha}{\diffEP \diffew}\\
a_W &= \frac{\alpha}{\diffPW \diffew}\\
b_P &= f_P
\end{align*}

Für äquidistante Gitter ergeben sich vereinfachte Terme für die Koeffizienten $a$.

\begin{align*}
  a_W &= \frac{\alpha}{\Delta x^2} = a\\
  a_E &= \frac{\alpha}{\Delta x^2} = a\\
  a_P &= - \frac{2 \alpha}{\Delta x^2} = -2a
\end{align*}

Betrachtet man nun alle Kontrollvolumen, kann mit diesen Beziehungen ein 
Gleichungssystem der Form $A\phi = b$ aufgestellt werden.

\section{Randbedingungen}
\label{sec:Randbedingungen}

Es wird angenommen, das der Randwert $\phi_w= \phi^0$ gegeben ist. Damit kann die
Ableitung $\frac{\partial \phi}{\partial x}$ wie folgt approximiert werden:

\begin{equation*}
\left(\frac{\partial \phi}{\partial x}\right)_W = \frac{\phi_P - \phi_W}{x_P - x_W} =
  \frac{\phi_P - \phi^0}{x_P - x_w}
\end{equation*}

Die Koeffizienten für ein Randkontrollvolumen ergeben sich damit zu:

\begin{align*}
  a_W &= 0\\
  a_P &= -\frac{\alpha}{\diffEP \diffew} - \frac{\alpha}{(x_P-x_w) \diffew}\\
  a_E &= \frac{\alpha}{\diffEP \diffew}\\
  b_P &= f_P - \frac{\alpha}{(x_P -x_w)\diffew}\phi^0
\end{align*}

beziehungsweise für den äquidistanten Fall:

\begin{align*}
  a_W &= 0\\
  a_P &=  -\frac{3 \alpha}{\Delta x^2} = -3a\\
  a_E &= \frac{\alpha}{\Delta x^2} = a\\
  b_P &= f_P - \frac{\alpha}{\Delta x \cdot \frac{1}{2} \Delta x}\phi^0 = f_P - 
  \frac{2\alpha}{\Delta x^2} \phi^0 = f_P - 2a \phi^0
\end{align*}


\begin{figure}
\begin{tikzpicture}
    \begin{axis}[
        height=7cm,
        width=0.9*\textwidth,
        grid=major,
        xlabel=$x$,
        ylabel=$\phi$,
    ]
    \addplot[tud2d, mark=*] file {data/solution_20_aequi.txt};
    \end{axis}
\end{tikzpicture}
\caption{Lösung des Diffusionsproblems mit Quellterm aus Gleichung \ref{eq:quellterm}}
\end{figure}

\begin{figure}
\begin{tikzpicture}
    \begin{axis}[
        height=7cm,
        width=0.9*\textwidth,
        grid=major,
        xlabel=$x$,
        ylabel=$\phi - \phi_{exakt}$,
%scaled y ticks = false,
    ]
    \addplot[tud2d, mark=*] file {data/error_20_aequi.txt};
    \end{axis}
\end{tikzpicture}
\caption{Fehlerverteilung der Lösung}
\end{figure}

\begin{figure}
\begin{tikzpicture}
    \begin{axis}[
        height=7cm,
        width=0.9*\textwidth,
        grid=major,
        xlabel=$x$,
        ylabel=$A\phi_{exakt}-b$,
%scaled y ticks = false,
    ]
    \addplot[tud2d, mark=*] file {data/residuum_20_aequi.txt};
    \end{axis}
\end{tikzpicture}
\caption{Residuum der Lösung}
\end{figure}


\chapter{1D Konvektionsgleichung mit FVM}
\label{chap:1D Konvektionsgleichung mit FVM}

Hier soll nun ein Konvektionsproblem nach Gleichung~\ref{eq:konvektion} betrachtet
werden. Hierbei wird $v$ als gegeben angenommen, da sonst ein iteratives
Lösungsverfahren nötig wäre.

\begin{equation}
  \rho v \frac{\partial \phi}{\partial x} = f
  \label{eq:konvektion}
\end{equation}

Gleichung~\ref{eq:konvektion} wird nun nach der finiten Volumen Methode über ein
Kontrollvolumen $V$ integriert. Danach wendet man den Satz von Gauß an und stellt
anschließend das Oberflächenintegral in die Summe der Integrale der Einzelflächen
überführt.

\begin{align}
  \int_V \rho v \frac{\partial \phi}{\partial x} dV &= \int_V f dV \nonumber \\
  \int_S \rho v \phi n_x dS &=\int_V f dV \nonumber\\
  \sum_c \int_{S_c} \rho v \phi n_{xc} dS_c &=\int_V f dV \label{eq:base_conv_ready}
\end{align}

Das Integral $\int_V f dV$ kann mit der Mittelpunktsregel über den Funtionswert
im Zentrum des Kontrollvolumens approximiert werden.

\begin{equation}
  \int_V f dV \approx f_P (x_e-x_w) \label{eq:source_term}
\end{equation}

Die Oberflächenintegrale müssen im eindimensionalen Fall für die Ost- und Westseite
ausgewertet werden.

\begin{align}
  \sum_c \int_{S_c} \rho v \phi n_{xc} dS_c &=
  \int_{S_e} \rho v \phi n_{xe} dS_e +\nonumber
  \int_{S_w} \rho v \phi n_{xw} dS_w\\
  &= \rho v \phi_e - \rho v \phi_w \label{eq:surface_east_west}
\end{align}

Die Gleichungen enthalten nun die Variablenwerte von $\phi$ an den Rändern der
Kontrollvolumen. Da hier jedoch von einer zellenorientierten Anordnung ausgegangen
wird, müssen diese im nächsten Schritt durch Werte in den Mittelpunkten des
Kontrollvolumens ausgedrückt werden.

\paragraph{Upwind-Verfahren}
Bei Upwind-Verfahren (auch Upwind Differencing Scheme, UDS) wird der Verlauf von $\phi$ durch eine Treppenfunktion approximiert, die
von der Geschwindigkeit $v$ abhängt. Es soll hier nur erwähnt, aber im Beispiel
nicht weiter berücksichtigt werden.

\begin{equation*}
\phi_e=\left\{\begin{array}{cl} \phi_P, & \mbox{falls }v>0\\
\phi_E, & \mbox{falls } v<0\end{array}\right.
\end{equation*}

\paragraph{Zentraldifferenzen-Verfahren}
Das Zentraldifferenzen-Verfahren (auch Central Differencing Scheme, CDS) wird der
Wert $\phi_e$ durch eine lineare Interpolation aus $\phi_P$ und $\phi_E$ gewonnen.

\begin{equation*}
  \phi_e \approx \phi_E\gamma_e + \phi_P (1-\gamma_e)
\end{equation*}

Der Interpolationsfaktor $\gamma_e$ ist hierbei gegeben als:

\begin{equation}
  \gamma_e = \frac{x_e-x_P}{x_E-x_P}
  \label{eq:cds_faktor}
\end{equation}

\paragraph{``Flux-Blending''-Verfahren}

Beim Flux-Blending versucht man Vor- und Nachteile von Verfahren höherer und niedrigerer
Ordnung zu verbinden. So hat das UDS-Verfahren einen Interpolationsfehler erster Ordnung,
während das CDS-Verfahren zu numerischen Oszillationen neigt.

Beim Flux-Blending werden nun beide Verfahren gewichtet zusammengeführt. Der Gewichtungsfaktor
wird hierbei $\beta$ genannt. Gilt $\beta = 1$, so ergibt sich das CDS-Verfahren, für
$\beta = 0$ das UDS-Verfahren.

\begin{equation}
  \phi_e \approx (1-\beta)\phi_e^{UDS} + \beta \phi_e^{CDS}
  \label{eq:flux_blending}
\end{equation}



\paragraph{Weitere Verfahren}
Neben UDS und CDS gibt es weitere Verfahren zur Approximation der Randvariablen durch
die Variablen im Zentrum eines Kontrollvolumens.

So wird beim QUICK-verfahren (Quadratic Upwind Interpolation for Convective
Kinematics) ein quadratisches Polynom, welches auf die Flussrichtung abgestimmt ist,
zur Approximation genutzt.

\section{Gesamtgleichungssystem UDS}

Wird die Treppenfunktion für das UDS-Verfahren in Gleichung~\eqref{eq:surface_east_west}
eingesetzt erbit sich folgende Gleichung:

\begin{equation}
  \rho v (\phi_P - \phi_W) = f_P(x_e-x_w)
\end{equation}

In Standardform gebracht ergeben sich damit folgende Koeffizienten für allgemeine
und äquidistante Gitter.

\begin{tabular}{cll}
  $a_W$ & $-\frac{\rho v}{x_e-x_w}$ & $-\frac{\rho v}{\Delta x}$ \\
  $a_P$ & $\frac{\rho v}{x_e-x_w}$ & $\frac{\rho v}{\Delta x}$ \\
  $a_E$ & $0$ & $0$ \\
  $a_W$ & $f_P$ & $f_P$ \\
\end{tabular}

\section{Gesamtgleichungssystem CDS}

Werden nun die Gleichungen \eqref{eq:surface_east_west} und \eqref{eq:source_term}
in Gleichung~\eqref{eq:base_conv_ready} eingesetzt und das Zentraldifferenzen-Verfahren
zur Approximation der Randvariablen genutzt, so ergibt sich:

\begin{equation*}
  \rho v \left[{\phi_E \frac{x_e-x_P}{x_E-x_P} + \phi_P \frac{x_E-x_e}{x_E-x_P}
  -\phi_W \frac{x_P-x_w}{x_P-x_W} - \phi_P \frac{x_w-x_W}{x_P-x_W}
  }\right]
  = f_P(x_e-x_w)
\end{equation*}

Diese Gleichung kann nun in die Form $a_P \phi_P + a_E \phi_E + a_W \phi_W = b_P$
gebracht werden. Mit den Koeffizienten $a_i$ und $b_i$ aller Kontrollvolumen kann
dann ein Gleichungssystem der Form $A\phi = b$ aufgestellt werden.

Die Koeffizienten nehmen bei allgemeinen, nicht äquidistanten Gittern folgende Werte
an:
\begin{align*}
  a_E &= \frac{\rho v}{x_e-x_w}\ \frac{x_e-x_P}{x_E-x_P}\\
  a_W &=-\frac{\rho v}{x_e-x_w}\ \frac{x_P-x_w}{x_P-x_W}\\
  a_P &= \frac{\rho v}{x_e-x_w} \left({\frac{x_E-x_e}{x_E-x_P}
-\frac{x_w-x_W}{x_P-x_W}}\right)\\
  b_P &= f_p
\end{align*}


Bei äquidistanten Gittern ergeben sich folgende Koeffizienten:

\begin{align*}
  a_E &= \frac{\rho v}{2\Delta x}\\
  a_W &=-\frac{\rho v}{2\Delta x}\\
  a_P &= \frac{\rho v}{\Delta x} \left({\frac{\Delta x}{2}
-\frac{\Delta x}{2}}\right) = 0\\
  b_P &= f_p
\end{align*}


\section{Randbedingungen}

Als Randbedingungen des Problems sollen $\phi_e = \phi^0$ und $\phi_w = \phi^1$ für die jeweiligen
Ränder gegeben sein. Bei Betrachtung von Gleichung~\eqref{eq:surface_east_west}
fällt auf, das für die Randkontrollvolumen nur noch eine Zentraldifferenz eingesetzt
werden muss. Es ergibt sich für den östlichen Rand die folgende Formel und die
Koeffizienten:

\begin{equation*}
  \rho v \left[{
  \phi^0  -\phi_W \frac{x_P-x_w}{x_P-x_W} - \phi_P \frac{x_w-x_W}{x_P-x_W}
  }\right]
  = f_P(x_e-x_w)
\end{equation*}

\begin{align*}
  a_E &= 0\\
  a_W &= -\frac{\rho v}{x_e-x_w}\ \frac{x_P-x_w}{x_P-x_W}\\
  a_P &= -\frac{\rho v}{x_e-x_w}\ \frac{x_w-x_W}{x_P-x_W}\\
  b_P &= f_P - \frac{\rho v}{x_e-x_w} \phi^0
\end{align*}

Beziehungsweise für äquidistante Gitter:

\begin{align*}
  a_E &= 0\\
  a_W &=-\frac{\rho v}{2\Delta x}\\
  a_P &=-\frac{\rho v}{2 \Delta x}\\
  b_P &= f_p- \frac{\rho v}{\Delta x}\phi^0
\end{align*}






\section{Gesamtgleichungssystem ``Flux-Blending''}

Werden die Vorschriften von UDS- und CDS-Verfahren in Gleichung~\eqref{eq:flux_blending}
eingesetzt sowie der Interpolationsfaktor~$\gamma_e$ aus Gleichung~\eqref{eq:cds_faktor} genutzt, so ergibt sich:

\begin{align}
  \phi_e &= (1-\beta)\phi_e^{UDS} + \beta \phi_e^{CDS} \nonumber\\
         &= (1-\beta)\phi_P + \beta \left({\phi_E\gamma_e + \phi_P (1-\gamma_e)}\right) \nonumber\\
         &= (1-\beta \gamma) \phi_P + \beta \gamma \phi_E\\
         &= \left({1-\beta \frac{x_e-x_P}{x_E-x_P}}\right) \phi_P + \beta \frac{x_e-x_P}{x_E-x_P} \phi_E
\end{align}

Für $\phi_w$ ergibt sich auf gleiche Weise:

\begin{align}
  \phi_w &= (1-\beta)\phi_w^{UDS} + \beta \phi_w^{CDS} \nonumber\\
         &= (1-\beta)\phi_W + \beta \left({\phi_W \gamma_w + \phi_P (1-\gamma_w)}\right) \nonumber\\
         &= (1-\beta +\beta \gamma_w)\phi_W + \beta (1-\gamma_w)\phi_P\\
         &= \left({1-\beta + \beta \frac{x_P-x_w}{x_P-x_W}}\right) \phi_W + \beta \frac{x_w-x_W}{x_P-x_W} \phi_P
\end{align}

Werden nun die Gleichungen \eqref{eq:surface_east_west} und \eqref{eq:source_term}
in Gleichung~\eqref{eq:base_conv_ready} eingesetzt sowie $\phi_e$ und $\phi_w$ über
die oben hergeleiteten Approximationen ersetzt, ergibt sich folgende Gleichung.

\begin{equation*}
  \rho v \left[{
    \left({1-\beta \frac{x_e-x_P}{x_E-x_P}}\right) \phi_P + \beta \frac{x_e-x_P}{x_E-x_P} \phi_E
    -\left({1-\beta + \beta \frac{x_P-x_w}{x_P-x_W}}\right) \phi_W - \beta \frac{x_w-x_W}{x_P-x_W} \phi_P
  }\right]
  = f_P(x_e-x_w)
\end{equation*}

Diese Gleichung kann nun in die Form $a_P \phi_P + a_E \phi_E + a_W \phi_W = b_P$
gebracht werden. Mit den Koeffizienten $a_i$ und $b_i$ aller Kontrollvolumen kann
dann ein Gleichungssystem der Form $A\phi = b$ aufgestellt werden.

Die Koeffizienten nehmen bei allgemeinen, nicht äquidistanten Gittern folgende Werte
an:

\begin{align*}
  a_E &= \frac{\rho v \beta}{x_e-x_w}\ \frac{x_e-x_P}{x_E-x_P}\\
  a_W &=-\frac{\rho v}{x_e-x_w}\ \left({1-\beta + \beta \frac{x_P-x_w}{x_P-x_W}}\right)\\
  a_P &= \frac{\rho v}{x_e-x_w} \left({
  \left({1-\beta \frac{x_e-x_P}{x_E-x_P}}\right)- \beta \frac{x_w-x_W}{x_P-x_W}
  }\right)\\
  b_P &= f_p
\end{align*}

Bei äquidistanten Gittern ergeben sich folgende Koeffizienten:

\begin{align*}
  a_E &= \frac{\rho v \beta}{2\Delta x}\\
  a_W &=-\frac{\rho v}{\Delta x}\left({1-\frac{\beta}{2} }\right)\\
  a_P &= \frac{\rho v}{\Delta x}(1-\beta)\\
  b_P &= f_p
\end{align*}

Hier ist auch eine gute Überprüfung der Ergebnisse möglich. Wie oben beschrieben
müssen sich für $\beta = 1$ für das ``Flux-Blending''-Verfahren die gleichen
Koeffizienten wie beim CDS-Verfahren ergeben. Einsetzen und ein Vergleich bestätigen dies.

\subsection{Randbedingungen}

\paragraph{Östlicher Rand}

\begin{align*}
  a_E &= 0\\
  a_W &=-\frac{\rho v}{x_e-x_w}\ \left({1-\beta + \beta \frac{x_P-x_w}{x_P-x_W}}\right)\\
  a_P &=-\frac{\rho v}{x_e-x_w} \left({\beta \frac{x_w-x_W}{x_P-x_W}}\right)\\
  b_P &= f_p-\frac{\rho v}{x_e-x_w} \phi^0
\end{align*}


\paragraph{Westlicher Rand}

\begin{align*}
  a_E &= \frac{\rho v \beta}{x_e-x_w}\ \frac{x_e-x_P}{x_E-x_P}\\
  a_W &=0\\
  a_P &= \frac{\rho v}{x_e-x_w} \left({1-\beta \frac{x_e-x_P}{x_E-x_P}}\right)\\
  b_P &= f_p+\frac{\rho v}{x_e-x_w} \phi^0
\end{align*}


\end{document}