\chapter*{Danksagung}

Hiermit möchte ich allen danken, die zum Gelingen dieser Diplomarbeit beigetragen haben. Im Besonderen gilt
mein Dank folgenden Personen:
\begin{itemize}
\item Professor Wilhelm Kley für die Betreuung dieser Diplomarbeit
\item Professor David Wharam für die Weckung der Neugierde an diesem Thema
\item Juha Ruokolainen und Peter R\aa{}back für ihre Unterstützung bei der Benutzung des Elmer-Softwarepakets
\item Volker Geyer und Daniel Duffner für ihre Lieferung von experimentellen Daten zu Lichtmühlen
\item Christin Flischikowski für das Korrekturlesen dieser Arbeit
\item meinen Eltern Elfriede Freund-Nadler und Harald Nadler für die finanzielle Unterstützung meines Studiums
\end{itemize}

\begin{center}
Vielen Dank
\end{center}

\tableofcontents

\chapter{Einleitung}

\section{Zusammenfassung}

Im Rahmen dieser Diplomarbeit wurden Erweiterungen für einen mit der Finiten"=Elemente"=Methode arbeitenden
Navier"=Stokes"=Löser geschrieben, um die Physik verdünnter Gase simulieren zu können. Es wurden im
wesentlichen zwei Effekte -- die thermische Transpiration nach Maxwell und der Temperatursprung nach
Smoluchowski -- implementiert. Nachdem erfolgreich Testfälle anderer Autoren nachgerechnet wurden, wurde
versucht die Drehmomente, die auf die Flügel einer Lichtmühle wirken, zu simulieren. Die Ergebnisse wurden
mit experimentellen Daten verglichen. Dabei ergaben sich zum Teil gute qualitative Übereinstimmungen mit
bestimmten gemessenen Proportionalitäten aber auch deutliche Unterschiede in den genauen Zahlenwerten.

\section{Was ist thermische Transpiration?}

Thermische Transpiration ist ein durch einen Temperaturgradienten hervorgerufener Gasfluss in Richtung der
höheren Temperatur. Er tritt aber nur in stark verdünnten Gasen auf. Diese Eigenschaft von Gasen wird durch
 die in Gleichung \ref{Kn} definierte Knudsenzahl $\kn$ charakterisiert. $d$ ist dabei eine
"`charakteristische"' Abmessung des Systems -- in einem langen Rohr wäre das der
Durchmesser -- und $\lambda$ die mittlere freie Weglänge der Gasteilchen.
\begin{equation}
\kn = \frac{\lambda}{d}
\label{Kn}
\end{equation}

In diesem Sinn bedeutet "`verdünnt"' also geringe Dichte und damit eine große mittlere freie Weglänge
$\lambda$ oder das Strömen um oder in sehr kleinen Geometrien.
Bei großen $\kn$, das heißt nahe 1, kann die thermische Transpiration zu einem dominierenden Effekt werden,
auch wenn sie sonst kaum eine Rolle spielt. So kann in Kanülen mit Durchmessern in der Größenordnung von
\textmu m eine
Druckdifferenz von mehreren 1000\,Pa erreicht werden auch wenn das Gas eine gewöhnliche Dichte hat. Ein
anderes Beispiel sind Lichtmühlen. Bei den großen mittleren freien Weglängen, wie sie in
den
Kammern von Lichtmühlen vorkommen, werden eben diese durch thermische Transpiration nur mit Tageslicht zur
Bewegung gebracht.



\section{Bedeutung des thermischen Transpirationseffekts}

Das heutige Interesse am Verhalten von Fluiden bei hohen Knudsenzahlen wird vor allem durch die
Fortschritte in der Miniaturisierung erklärt.  Miniaturisierte Geräte, die mechanische und elektrische
Effekte kombinieren -- engl. "`microelectromechanical systems"' (MEMS), sollen in Zukunft für eine Fülle von
Anwendungen nützlich sein. Vorstellbar sind zunächst viele Möglichkeiten, wie mit MEMS die Messungen von
Effekten in Mikro- und Nanostrukuren verbessert werden. So könnten in Zukunft ganze Versuchsaufbauten auf
Mikrochips passen, die Dank der Massenfertigungskapazitäten der Halbleiterindustie sogar relativ günstig sein
könnten. Großes Potential wird auch in der Medizintechnik gesehen. Mit MEMS könnten Körperfunktionen
überwacht,
unterstützt oder ersetzt werden.
Das Verständnis von Effekten bei hohen $\kn$"=Werten ist auch wichtig, selbst wenn man sie nicht bewust
ausnutzen
will. So stoßen zum Beispiel Festplatten heute in Größenskalen vor, bei denen Wissen über diese Effekte nötig
ist.

Mit thermischer Transpiration im Speziellen werden sogenannte Knudsenpumpen gebaut. Durch Anlegen eines
Temperaturgradienten werden Druckdifferenzen erzeugt. Da Knudsenpumpen ohne bewegliche mechanische Teile
auskommen, erwartet man sehr gute Miniaturisierungsmöglichkeiten. Des Weiteren benötigen sie keine
Schmierstoffe oder Arbeitsgase und versprechen daher, das zu verdünnende oder zu verdichtende Gas nicht zu
verunreinigen. Auf Grund ihrer geringen Größen können leicht mehrere in Reihe geschaltet werden um große
Druckgradienten zu erreichen.

\section{Geschichte des thermischen Transpirationseffekts und der Lichtmühlen}

Die Lichtmühle besteht aus einem Rotor, der sich in einem Glaskolben befindet, der mit Luft sehr geringer
 Dichte beziehungsweise Druck gefüllt ist. Abbildung \ref{fig:lightmill} zeigt eine handelsübliche Lichtmühle.
Die vier senkrecht
stehenden, am Rotor befestigten Flügel sind je auf einer Seite silbrig glänzend und auf der anderen schwarz;
meist ist das Schwarz einfach eine Schicht Ruß. Wenn die Lichtmühle beleuchtet wird, dreht sich der Rotor
wobei die silbrigen Seiten in Drehrichtung zeigen.

\begin{figure}
\capstart\centering
\includegraphics[width=0.34\linewidth]{fig/lightmill.jpg}
\caption{Eine Lichtmühle. Das Foto stammt von Wikipedia.}
\label{fig:lightmill}
\end{figure}

Die Lichtmühle wurde um 1870
von William Crookes -- daher auch der alternative Name "`Crookes Radiometer"' -- zum ersten mal ausführlicher
wissenschaftlich untersucht. Ursprünglich wurden sie zur Messung von Lichtintensität verwendet. Es kursieren
bis heute viele falsche Erklärungen zur Funktionsweise der Lichtmühle, obwohl schon bis Ende der 70er
Jahre des 19. Jahrhunderts
praktisch alle falschen Erklärungen widerlegt wurden. Crookes selbst versuchte die Bewegung der Flügel
zunächst mit Lichtimpuls zu erklären. Das ist offensichtlich falsch, da sich die Flügel dann mit den
schwarzen Seiten
voran bewegen müssten, was nicht der Fall ist. Wenn man perfekte Absorption von Licht auf den schwarzen Seiten
und perfekte Reflexion auf den silbrigen Seiten annimmt, muss die silbrige Seite den doppelten Impuls der
schwarzen Seite bekommen und somit müssen sich die schwarzen Seiten voran bewegen.

Eine weitere falsche Erklärung lautet: Das Gas nahe der schwarzen Seite hat eine höhere Temperatur als das
auf der silbrigen. Da der Molekülimpuls mit der Temperatur wächst, treffen die Luftmoleküle mit höherem Impuls
auf die schwarze Seite als auf die silbrige. Es gibt dann einen netto Impuls, der die Flügel mit den
silbrigen Seiten voran bewegt. Der Fehler in dieser Erklärung ist die Tatsache, dass die höhere Temperatur
auch zu
einer niedrigeren Dichte auf der schwarzen Seite führt, weil der Druck auf Grund der offenen Geometrie auf
beiden Flügelseiten gleich ist. Somit haben die Teilchen auf der schwarzen Seite zwar einen höheren Impuls
aber
es sind auch weniger. Genaue Rechnungen zeigen, dass sich diese beiden Effekte gerade wegheben und es kein
netto Impulsübertrag auf den Flügel gibt.

Auch dass die Bewegung durch Ausdampfen von Gasen aus den Flügeloberflächen verursacht wird, ist falsch.

1878 und 1879 erschienen Artikel von Osborne Reynolds und James Clerk Maxwell \cite{maxwell}, die die
Bewegung der
Lichtmühle mit einem ganz neuen Effekt, der thermischen Transpiration, erklären. Diese Idee stammt von
Reynolds: er war der erste, der mit Überlegungen zur kinetischen Gastheorie auf die Existenz der thermischen
Transpiration schloss und auch den Begriff prägte. Maxwells Arbeiten waren deutlich ausführlicher und besser
begründet. Von ihm stammen die Geschwindigkeitsrandbedingungen, die direkt in hydrodynamische Gleichungen
eingesetzt werden können und dort für ein Gleiten des Gases auf der Wand in Richtung eines
Temperaturgradienten
sorgen.

Um 1910 wurden ausführliche Experimente von Knudsen zur Untersuchung der thermischen Transpiration
durchgeführt. Er verwendete unter anderem durch poröse Materialien (Meerschaum) verbundene Kammern. Wurde
eine Kammer erwärmt, kam es zu einem Gasfluss durch den Meerschaum in die wärmere Kammer. Er erwähnte auch die
Möglichkeit, mit diesem Effekt Pumpen ohne mechanische Teile zu bauen.

In den zwanziger Jahren gab es weitere Artikel zur Natur der Lichtmühlenkräfte -- unter anderem von Albert
Einstein und Paul Epstein. In den dreißiger Jahren erschienen Lehrbücher über kinetische Gastheorie von Loeb
und Kennard \cite{kennard}, die die damaligen Erkenntnisse zu Effekten in verdünnten Gasen und
der Kräfte in Lichtmühlen in je einem eigenen Kapitel zusammenfassten. Das Resümee war, dass verschiedene
Effekte bei der Bewegung der Lichtmühle eine Rolle spielen und die thermische Transpiration in jedem Fall
einen wichtigen Teil beiträgt.




\chapter{Theorie}
\label{theorie}

\section{Hydrodynamik bei hohen Knudsenzahlen}

Hohe Knudusenzalen $\kn = \frac{\lambda}{d}$ bedeuten, dass die mittlere freie Weglänge $\lambda$ groß ist im
Vergleich zu einer charakteristischen Systemgröße $d$. In Rohrströmungen wird typischerweise der Durchmesser
des Rohres, bei Strömung durch Quader dessen Höhe als $d$ definiert. Bei offeneren Geometrien gibt es keine
eindeutigen Regeln zu Bestimmung von $d$ in der Literatur. Für die Hydrodynamik bei hohen $\kn$"=Werten müssen eigene
Modelle und Theorien entwickelt werden, da etablierte Theorien wie die Navier-Stokes-Gleichungen, dort
falsche Ergebnissen liefern.

Wenn $\lambda$ groß ist im Vergleich zu $d$ dann wechselwirken die Moleküle vor allem mit den Wänden und
deutlich seltener mit sich selbst als bei kleinen $\kn$"=Werten. Dadurch ist klar, dass Rand- beziehungsweise
Oberflächeneffekte eine vorherrschende Rollen spielen und deren möglichst genaue Modellierung sehr wichtig
ist.
Eine
weitere Konsequenz aus großen mittleren freien Weglängen ist die Nicht"=Gültigkeit eines lokalen
thermodynamischen Gleichgewichts.

Viele thermodynamische beziehungsweise hydrodynamische Gleichungen sind unter der
Annahme hergeleitet worden, dass sich das zu beschreibende System im thermodynamischen Gleichgewicht befindet
oder in einer sehr kleinen Abweichung davon. Im Fall einer Abweichung vom thermodynamischen Gleichgewicht,
soll zumindest ein lokales gelten. Das heißt, dass sich das Fluid, auch wenn makroskopische thermodynamische
Variablen wie Druck oder Temperatur sich global im System unterscheiden, in einem gewissen "`kleinen"'
Untervolumen in der durch die dortigen Werte der thermodynamischen Variablen bestimmten
Gleichgewichtsverteilung befindet. In anderen Worten bedeutet das, dass sich die Moleküle durch
die Maxwell"=Boltzmann"=Verteilung beschreiben lassen.

Wenn $\lambda$ groß ist, befinden sich auch in einem kleinen Volumen des Fluids Moleküle, die weitgehend
ungestört von Gebieten mit unterschiedlichen thermodynamischen Größen gekommen sind. Somit gibt es kein
lokales
Gleichgewicht und eine wichtige Voraussetzung für viele thermo- beziehungsweise hydrodynamische Modelle, wie
zum
Beispiel die Navier"=Stokes"=Gleichungen, sind nicht mehr gültig.

\subsection{Hydrodynamische Regime}

Es gibt zwei verschiedene Herangehensweisen um Modelle für Systeme mit hohen Knudsenzahlen zu erzeugen:
Entweder versucht man Kontinuumsgleichungen anzupassen beziehungsweise neue zu finden und mit entsprechenden
Randbedingungen auszustatten und sie somit in Übereinstimmung mit der Physik bei großen $\kn$"=Werte zu
bringen,
oder man geht von mechanischen Modellen aus und versucht diese mit entsprechenden Näherungen für immer
größere Teilchenzahlen praktikabel zu machen.

In Tabelle \ref{tab:hydroReg} wird eine typische Übersicht der Gültigkeit beziehungsweise Praktikabilität von
physikalischen Modellen bei bestimmten $\kn$ gezeigt, wie sie es zum Beispiel in \cite{karniadakis2005man}
gibt. Man nennt diese Einteilungen auch hydrodynamische Regime. Tabelle \ref{tab:hydroReg} ist nur als eine
sehr grobe Einteilung zu verstehen, die auf bisher gemachten Erfahrungen mit Simulationen und Experimenten
beruht. Im Prinzip ist für jeden Effekt und jede Geometrie eine eigene Untersuchung von passenden Modellen bei
verschiedenen $\kn$"=Werten nötig.

\begin{table}
\capstart\centering
\begin{tabular}{p{0.23\linewidth}|l|p{0.4\linewidth}}
Name & $\kn$ & Übliche Modelle \\
\hline
Klassisches Kontinuums"=Regime & $\kn < 0.001$ & Navier-Stokes-Gleichungen\\
\hline
Slip-Regime & $0.001 < \kn < 0.1$ & N-S-Gl. mit angepassten Randbedingungen z.B. Slip\\
\hline
Übergangsregime & $0.1 < \kn < 10$ & Mögliche Kandidaten u.a.: DSMC, Burnett-Gl, Lattice"=Boltzmann \\
\hline
Freier Molekülfluss & $10 < \kn$ & Mechanische Teilchenmodelle \\
\end{tabular}
\caption{Liste hydrodynamischer Regime mit möglichen Modellen zu darin stattfindenden physikalischen
Vorgängen}
\label{tab:hydroReg}
\end{table}

Eine Möglichkeit zur Herleitung von Kontinuumsgleichung, die besser für hohe Knudsenzahlen
geeignet sind, ist die Chapman"=Enskog"=Methode, die eine Reihe von näherungsweisen Lösungen der
Boltzmann"=Gleichungen erzeugt. In Kapitel \ref{ch-en} wird darauf näher eingegangen.


Ein Beispiel für ein mechanisches Teilchenmodell, das auch noch die relativ großen Teilchenzahlen des
Übergangsregimes auf aktuellen Rechnern handhaben kann, ist die "`Direct Simulated Monte Carlo"'"=Methode
(DSMC). In dieser Methode werden nur Wand"=Teilchen"=Wechselwirkungen deterministisch berechnet. Stöße
zwischen Teilchen werden dagegen stochastisch modelliert, in dem allen Teilchen in einem diskretisierten
Volumenteil $\Delta V$ nach einem Zeitschritt $\Delta t$ neue Geschwindigkeiten $\bm v_i$ zugewiesen werden.
Die neue Geschwindigkeitsverteilung entspricht den im Gittervolumen $\Delta V$ herrschenden makroskopischen
Zustandvariablen. Welches Teilchen welche der neuen $\bm v_i$ bekommt, wird mit einem Zufallsgenerator
bestimmt. Außerdem steht ein simuliertes Teilchen für mehrere echte Moleküle.

Eine wiederum andere Methode hat den Namen Lattice"=Boltzmann. Mit diesem Schema wird die (linearisierte)
Boltzmann"=Gleichung direkt numerisch gelöst, anstatt erst Näherungsgleichungen zu Entwickeln um diese
numerisch zu lösen.


In dieser Diplomarbeit wurden die normalen Navier"=Stokes"=Gleichungen und Wärmeleitgleichung benutzt und mit
Maxwells Thermischen"=Transpirationsrandbeding"=ung und Smoluchowskis Temperatursprungrandbedingung erweitert.
Obwohl dieses Modell  in der Literatur nur für $\kn$ Zahlen bis zur Größenordnung 0.1 empfohlen wird, wurde
in dieser Arbeit versucht, physikalische Vorgänge im Bereich $\kn \approx 1$ zu simulieren.

\subsection{Weitere Effekte bei hohen Knudsenzahlen}

Es gibt eine Vielzahl von weiteren Effekten, die bei kleinen $\kn$ fast keine, bei großen $\kn$ dagegen eine
wichtige Rolle spielen. Einige sind bis jetzt nur theoretisch vorhergesagt worden und sehr schwer in
Experimenten nachzuweisen beziehungsweise von einander zu unterscheiden. Ein Beispiel wäre der "`thermal
stress slip
flow"', der sich theoretisch aus einer Entwicklung der Boltzmann"=Gleichungen bis zur 2. Ordnung ergibt und
für
eine Korrektur der Slip"=Geschwindigkeit abhängig von einem Temepraturgradient senkrecht zur Wand sorgt. Somit
kann dieser Effekt für Gasbewegungen sorgen auch falls auf der Wand eine gleichförmige Temperatur herrscht und
keine thermische Transpiration stattfindet.

\section{Phänomenologische Ableitung der Navier"=Stokes"=Gleichungen}
\label{N-S-Gl}

Man kann die Navier-Stokes-Gleichungen als Verallgemeinerung der Newtonschen Mechanik von Punktmassen $m_i$
auf eine kontinuierliche Dichteverteilung $\rho (\bm r)$ unter der Berücksichtigung der Erhaltungssätze
für Masse, Energie und Impuls sehen. Entsprechend besteht die Navier-Stokes-Gleichung aus drei Teilen, die
jeweils Massen-, Energie- und Impulserhaltung ausdrücken. Die folgende Herleitung der
Impulserhaltungsgleichung basiert auf \cite{huang1964sm}.

Gleichung \ref{newFluid} ist die Bewegungsgleichung entsprechend dem zweiten Newtonschen Gesetz für ein
kleines Fluidelement mit dem Volumen $\d V$ und der Geschwindigkeit $\bm u(\bm r, t)$ und der Masse $m$.
\begin{equation}
m \ddt \bm u = \bm F_\mathrm{ges}
\label{newFluid}
\end{equation}
Wird die Zeitableitung in Gleichung \ref{newFluid} ausgeführt und $m$ mit $\rho (\bm r, t) \d V$ und die
gesamte auf $\bm u$ wirkende Kraft $\bm F_\mathrm{ges}$ durch $( \bm F_\mathrm{ext} + \bm G ) \d V$ ersetzt,
ergibt sich
Gleichung \ref{newFluid2}.
\begin{equation}
\rho \left( \pddt + \bm u \nabla \right) \bm u = \bm F_\mathrm{ext} + \bm G
\label{newFluid2}
\end{equation}
$\bm F_\mathrm{ext}$ ist eine äußere Kraft pro Volumen und $\bm G$ eine durch innere Reibung und statischen
Druck
hervorgerufene Kraft pro Volumen.

Zur Herleitung von $\bm G$ wird ein Fluid"=Würfel mit Kanten parallel zu den kartesischen Achsen angenommen,
auf dessen Seiten die von den Nachbarwürfeln verursachten Kräfte pro \textit{Fläche} $\bm T_i$ wirken. Die
Vorzeichen der $\bm T_i$ hängen von der gewählten Konvention für die Richtung der Flächennormalen ab. 
Auf die beiden Seiten senkrecht zu $x_i$ wirkt dann Ausdruck \ref{Ti}.
\begin{equation}
\bm T_i, \quad -\left( \bm T_i + \frac{\partial \bm T_i}{\partial x_i} \d x_i \right) \quad i \in \{ 1,2,3 \}
\label{Ti}
\end{equation}
Somit ergibt die gesamte, von allen Seiten aus auf den Fluid"=Würfel wirkende Kraft pro Volumen $\bm G$
Gleichung \ref{GdV}.
\begin{equation}
\bm G \, \d V = - \left( \frac{\partial \bm T_1}{\partial x_1} + \frac{\partial \bm T_2}{\partial x_2} +
\frac{\partial \bm T_3}{\partial x_3} \right) \d V
\label{GdV}
\end{equation}

Werden die Komponenten der Vektoren $\bm T_i$ wie in Gleichung \ref{TSi} benannt, kann man die Komponenten von
$\bm G$ wie in Gleichung \ref{Gi} ausdrücken.
\begin{align}
\bm T_1 &= ( S_{1,1}, S_{1,2}, S_{1,3} )\\
\bm T_2 &= ( S_{2,1}, S_{2,2}, S_{2,3} )\\
\bm T_3 &= ( S_{3,1}, S_{3,2}, S_{3,3} )
\label{TSi}
\end{align}
\begin{equation}
G_i = - \sum_{j=1}^3 \frac{\partial S_{i,j}}{\partial x_j} \quad \text{abgekürzt als} \quad \bm G = -\nabla
\mathbf S
\label{Gi}
\end{equation}

Nimmt man an, dass das betrachtete Fluid isotrop ist, muss $S_{1,1} = S_{2,2} = S_{3,3} \equiv P$ gelten. $P$
ist dabei gerade der hydrostatische Druck. Man kann $S_{i,j}$ formal, wie in Gleichung \ref{Pstrich}, in
diagonale und nicht"=diagonale Elemente aufspalten.
\begin{equation}
S_{i,j} = \updelta_{i,j} P + S_{i,j}'
\label{Pstrich}
\end{equation}
Die Matrix $\mathbf S'$ hat offensichtlich die Spur 0. Wenn das betrachtet Fluid im Volumen $\d V$ kein
inneren Drehimpuls hat, folgt dass $\mathbf S$ und damit $\mathbf S'$ symmetrische Matrizen sind. Diese
Annahme macht Sinn, denn das Volumen $\d V$ soll ja beliebig klein sein.

An dieser Stelle wird der Zusammenhang zwischen wirkenden Scherkräften und Deformationsgeschwindigkeiten des
Fluids als empirische Tatsache eingeführt. Eine auf den Fluidwürfel parallel zu einer Koordinate wirkende
Scherspannung $s$ streckt ihn mit der Winkelgeschwindigkeit $\frac{\d \phi}{\d t}$ in ein Parallelepiped.
Es gilt $s = \mu \frac{\d \phi}{\d t}$.
$\mu$ ist dabei die dynamische Viskosität und $\phi$ der Winkel, um den der Würfel in einen Parallelepiped
gekippt wurde. Wenn die angreifende Scherspannung $S_{2,1}'$ ist, dann ergibt sich diese wie in Gleichung
\ref{S21} beziehungsweise eine beliebige Komponente von $\mathbf S'$ wie in Gleichung \ref{Sstrich}.

\begin{equation}
S_{2,1}' = -\mu \left( \frac{\d \phi_1}{\d t} + \frac{\d \phi_2}{\d t} \right)
= -\mu \left( \frac{\partial u_2}{\partial x_1} + \frac{\partial u_1}{\partial x_2} \right)
\label{S21}
\end{equation}
\begin{equation}
S_{i,j}' = -\mu \left( \frac{\partial u_i}{\partial x_j} + \frac{\partial u_j}{\partial x_i} \right)
\label{Sstrich}
\end{equation}

Die Definition von $\mathbf S'$ Gleichung \ref{Sstrich} ist noch unvollständig, weil sie im Gegensatz zur
Konstruktion von $\mathbf S'$ in Gleichung \ref{Pstrich} nicht die Eigenschaft $\Sp \mathbf S' = 0$ hat.
$\mathbf S'$ muss noch um einen Term erweitert werden, damit die Spur wieder 0 ist. Somit ergibt sich
$\mathbf S'$ zu Gleichung \ref{Sstrich2} und damit die vollständige Matrix $\mathbf S$ wie in Gleichung
\ref{S}.
\begin{equation}
S_{i,j}' = -\mu \left( \frac{\partial u_i}{\partial x_j} + \frac{\partial u_j}{\partial x_i} \right)
- \updelta_{i,j} \frac{2}{3} \mu \nabla \bm u 
\label{Sstrich2}
\end{equation}
\begin{equation}
S_{i,j} = -\mu \left( \frac{\partial u_i}{\partial x_j} + \frac{\partial u_j}{\partial x_i} \right)
- \updelta_{i,j} \frac{2}{3} \mu \nabla \bm u + \updelta_{i,j} P
\label{S}
\end{equation}

Jetzt ist Gleichung \ref{newFluid2} vollständig bestimmt und wird zu Gleichung
\ref{N-S-Kompl}. Darin gilt $\bm F_\mathrm{ext} = \rho \bm a_\mathrm{ext}$ mit einer externen Beschleunigung
$\bm a_\mathrm{ext}$, weil $\bm F_\mathrm{ext}$ eine externe Kraft pro Volumen ist.

\begin{equation}
\rho \left( \pddt \bm u + \bm u \nabla \bm u \right)  = - \nabla \mathbf S + \bm F_\mathrm{ext}
\label{N-S-Kompl}
\end{equation}


\section{Herleitung von Transporterscheinungen aus der statistischen Physik}

Wenn ein System nicht im thermodynamischen Gleichgewicht ist beziehungsweise die Verteilungsfunktion $f$ nicht
die
Maxwell"=Verteilung ist, wird das System durch den Transport von Energie, Masse oder Impuls ins Gleichgewicht
kommen. In diesem Abschnitt wird genauer auf die von Nicht"=Gleichgewichtszuständen hervorgerufenen
Transporterscheinungen eingegangen. Es wird gezeigt, dass die Boltzmann"=Gleichung, welche die zeitliche
Entwicklung eines Systems -- und damit dessen Weg ins Gleichgewicht -- auf eine sehr allgemeine
Art beschreiben, durch
eine Reihenentwicklung angenähert werden kann. Diese Entwicklung nennt man Chapman"=Enskog"=Methode und
findet nach Potenzen der mittleren freien Weglänge statt. Die thermische Transpiration in der mathematischen
Formulierung wie sie in dieser Diplomarbeit implementiert wurde, ist ein Transportphänomen erster Ordnung.
Die Herleitungen in diesem Kapitel orientieren sich an der Darstellung in \cite{huang1964sm}.

\subsection{Die Boltzmann-Gleichung}

Es ist eine der wesentlichen Erkenntnisse der statistischen Mechanik, dass sich alle thermodynamischen Größen
eines Systems aus der Kenntnis einer Verteilungsfunktion $f( \bm r, \bm v , t)$ herleiten lassen. $f$ gibt an,
wie viele Teilchen sich in einem bestimmten Zustand $(\bm r, \bm v, t)$ befinden, so dass die gesamte Zahl der
Teilchen $N$ im System zu einem Zeitpunkt $t_1$ gerade $\int_\infty^\infty f( \bm r, \bm v , t_1)\, \d V \d^3
v$
ist. Es ist also nicht nötig die exakten Orte und Geschwindigkeiten aller Teilchen zu kennen, sondern es
genügt zu wissen wieviele Teilchen sich in einem bestimmten Orts- Geschwindgkeitszustand befinden. Weil $f$
eine Wahrscheinlichkeitsdichte ist, kann man den Erwartungswert einer makroskopischen Größe $A$ mit
Gleichung \ref{Aausf} bestimmen.

\begin{equation}
\langle A \rangle = \frac{\int A f \,\d^3 v}{\int f \,\d^3 v}
\label{Aausf}
\end{equation}


Die Boltzmann"=Gleichung \ref{boltzmann} beschreibt, wie sich die Verteilungsfunktion $f$ zeitlich entwickelt.
Die in Gleichung \ref{boltzmann} angegebene Form gilt für ein System mit einer Sorte Teilchen. $\mt$ ist die
Masse dieser Teilchensorte, $\bm F$ eine äußere Kraft pro Teilchen. Der Kollisionsterm $I_\mathrm{coll}$
enthält die
Beschreibung der Teilchenwechselwirkungen.
\begin{equation}
\left( \pddt + \bm v \nabla + \frac{\bm F}{\mt} \nabla_{\bm v} \right) f( \bm r, \bm v , t) =
I_\mathrm{coll}(f)
\label{boltzmann}
\end{equation}

Die Boltzmann"=Gleichung ist für die meisten sinnvollen Kollisionsterme $I_\mathrm{coll}$ eine schwer zu
handhabende Intero"=Differenzialgleichung und wird daher selten zur Lösung realistischer Probleme genutzt. Das
sieht man, wenn man $I_\mathrm{coll}$ für ein einfaches ideales Gas angibt. Dazu sollen folgende Näherungen
gelten:
\begin{itemize}
\item
Es gibt nur Zweierstöße, das heißt es wirken nie 3 Teilchen gleichzeitig miteinander.
\item
Wände werden vernachlässigt.
\item
Der Einfluss einer äußeren Kraft $\bm F$ auf die Teilchenstöße werden vernachlässigt.
\item
Geschwindigkeit und Ort eines Teilchens sind unkorreliert: es gilt molekulares Chaos. Dadurch kann man die
Zahl der Teilchenpaare, die in einem Volumenelement $\d V$ liegen und deren Geschwindigkeiten $\bm v_1$
und $
\bm v_2$ in den Volumenelementen des Geschwindigkeitsraums $\d^3 v_1$ und $\d^3 v_2$ liegen, als $f( \bm r,
\bm
v_1 , t)\, \d V \d^3 v_1 \cdot f( \bm r, \bm v_2 , t)\, \d V \d^3 v_2 $ schreiben.
\end{itemize}
Dieser auf Boltzmann selbst zurückgehende Ansatz führt auf den Kollisionsterm in Gleichung \ref{Icoll}.
$\sigma$ ist darin der differenzielle Wirkungsquerschnitt für den Zweierstoß, der $(\bm v_1, \bm v_2 )$ zu 
$(\bm v_1', \bm v_2')$ macht.


\begin{equation}
I_\mathrm{coll} = \iint \sigma(\varOmega) \left| \bm v_1 - \bm v_2 \right| \left( f( \bm r, \bm v_1' , t)
 f( \bm r, \bm v_2' , t) - f( \bm r, \bm v_1 , t)f( \bm r, \bm v_2 , t)\right) \, \d^3 v_2 \, \d \varOmega
 \label{Icoll}
\end{equation}



\subsubsection{Die Verteilungsfunktion im Gleichgewicht (Maxwell"=Verteilung)}


Die Lösung der Boltzmann"=Gleichung für $t \rightarrow \infty$ ist die Gleichgewichtsverteilung $f_0$. Sie
wird Maxwell"=Verteilung genannt. Da die Verteilungsfunktion für ein abgeschlossenen System im
thermodynamischen Gleichgewicht weder vom Ort noch von der Zeit abhängen kann, vereinfacht sich die
Boltzmann"=Gleichung auf Gleichung \ref{bg0}. Und da ein Integral für ein Integrationsintervall ungleich 0 nur
0 sein kann, wenn der Integrand 0 ist, gilt die noch einfachere Gleichung \ref{bg02}. Bei dieser
Vereinfachung wurde vorausgesetzt, dass $\sigma > 0$ gilt, was für reale Teilchen offensichtlich zutrifft.

\begin{equation}
0 = \iint \sigma(\varOmega) \left| \bm v_1 - \bm v_2 \right| \left( f( \bm v_1') f(\bm v_2') - f(\bm
v_1)f(\bm v_2)\right) \, \d^3 v_2 \, \d \varOmega
\label{bg0}
\end{equation}

\begin{equation}
f( \bm v_1') f(\bm v_2') - f(\bm v_1)f(\bm v_2) = 0
\label{bg02}
\end{equation}

Nimmt man den Logarithmus von \ref{bg02} erhält man Gleichung \ref{herMax1}. Darin sieht man, dass sich $\ln
f_0$ gleich verhält wie die Erhaltungsgrößen vor und nach dem Stoß $(\bm v_1,\bm v_2) \rightarrow
(\bm v_1',\bm v_2')$, denn die Erhaltungsgrößen $\chi$ verhalten sich bei einem Stoß wie $\chi (\bm v_1) +
\chi(\bm v_2) \rightarrow \chi (\bm v_1') + \chi(\bm v_2')$.

\begin{equation}
\ln f_0(\bm v_1) + \ln f_0(\bm v_2) = \ln f_0(\bm v_1') + \ln f_0(\bm v_2')
\label{herMax1}
\end{equation}

Wenn die Erhaltungsgrößen beziehungsweise Stoßinvarianten und $\ln f_0$ sich bei einem Stoß gleich verhalten
und beide
 nur von der Geschwindigkeit $\bm v$ abhängen, muss es eine lineare Funktion geben, die die beiden verknüpft.
Somit hängt $\ln f_0$ durch eine Linearkombination von $\bm v$ (Impulserhaltung), $v^2$ (Energieerhaltung)
und einer Konstante, die wir einfach $\ln C$ (Massenerhaltung) nennen, ab. Ein solcher Ansatz für $\ln f_0$
ist

\begin{equation}
\ln f_0(\bm v) = -A(\bm v - \bm v_0 )^2 + \ln C \Longleftrightarrow f_0(\bm v) = C \e^{-A(\bm v - \bm v_0
)^2}.
\label{herMax2}
\end{equation}

Man kann die Konstanten aus Gleichung \ref{herMax2} bestimmen, indem man mit Gleichung \ref{Aausf}
makroskopische Größen berechnet, deren Wert man für den Gleichgewichtszustand kennt. So weiß man, dass $\bm u
= \left\langle \bm v \right\rangle = 0$ im Gleichgewicht gilt.
Daraus ergibt sich, dass $\bm v_0 = 0$ sein muss. Weitere Integrale, die man löst, wären das zur Bestimmung
von
der Teilchendichte $n$, das $N/V$ ergeben muss, und das zur Bestimmung der mittleren Energie pro Teilchen
$\epsilon$. Sind so
alle
Konstanten bestimmt, ist die Verteilungsfunktion des Gleichgewichts -- also die Maxwell-Verteilung

\begin{equation}
f_0(\bm v) = \frac{N}{V}\left(  \frac{\mt}{2\uppi kT}\right)^{\frac{3}{2}} \exp \left( \frac{m \bm
v^2}{2kT}\right).
\label{maxwellVer}
\end{equation}

\subsection{Transporterscheinungen}

Auf dem Weg ins thermodynamische Gleichgewicht werden in einem Fluid so lange Masse, Impuls und Energie über
Stöße ausgetauscht, bis alle makroskopischen Größen überall gleich sind. Die mittlere freie Weglänge
$\lambda$,
also der Weg, den ein Teilchen durchschnittlich zwischen zwei Stößen zurücklegt, beziehungsweise die mittlere
Stoßzeit
$\tau$, also die Zeit, die durchschnittlich zwischen zwei Stößen vergeht, bestimmen die Geschwindigkeit mit
der ein lokales Ungleichgewicht ausgeglichen wird.

In Gleichung \ref{mfwDef} wird $\lambda$ mit der
Teilchendichte $n$, der Stöße pro Volumen pro Zeit $Z = 4n^2 \sigma_\mathrm{tot} \sqrt{\frac{kT}{\uppi \mt}}$
und der
mittleren Geschwindigkeit $ \bar v$ definiert. $Z$ kann man aus dem
Kollisionsterm $I_\mathrm{coll}$ berechen und $\bar v$ aus der Maxwellverteilung, was zu Gleichung
\ref{mg} führt. $\tau$ ist dann entsprechend
$\tau = \lambda / \bar v$. 
\begin{equation}
\lambda = \frac{n}{2Z} \bar v
\label{mfwDef}
\end{equation}
\begin{equation}
\bar v= \sqrt{\frac{8kT}{\uppi\mt}}
\label{mg}
\end{equation}
Wird nun $\bar v$ und das aus $I_\mathrm{coll}$ errechnete $Z$ in Gleichung \ref{mfwDef} eingesetzt, ergibt
sich Gleichung \ref{mfwDef2}. Der rechte Term in Gleichung \ref{mfwDef2} ergibt sich, wenn man noch $n$ mit
der Hilfe der idealen Gasgleichung ersetzt und $\sigma_\mathrm{tot} = \uppi \sigma^2$ setzt, wobei
$\sigma$ der Teilchendurchmesser im Modell harter Kugeln ist. Streng genommen gilt das so bestimmte $\lambda$
wieder nur für ein System im Gleichgewicht, also wenn $f$ die Maxwell"=Verteilung ist.

\begin{equation}
\lambda = \frac{1}{ \sqrt 2 n \sigma_\mathrm{tot}} = \frac{kT}{\sqrt{2} \uppi \sigma^2 P}
\label{mfwDef2}
\end{equation}

Mit Gleichung für $\tau$ kann man abschätzen wie schnell sich Ungleichheiten von makroskopischen Größen in
der Entfernung von $\approx \lambda$ ausgleichen. Im Fall von $\mathrm
H_2$ bei Normaldruck und Zimmertemperatur ist $\lambda \approx 1 \, \text{\textmu}\mathrm{m}$ und damit wären
Unterschiede in dieser Entfernung nach $\tau \approx 1 \, \mathrm{ns}$ ausgeglichen.

Um Transportphänomene beschreiben zu können, soll nun eine zeitabhängige Lösung $f$ der Boltzmann"=Gleichung
gefunden werden. Weil es in unserem System nur elastische Stöße gibt, muss jede Lösung $f$ bestimmte aus den
Erhaltungssätzen abgeleitete Eigenschaften erfüllen. Wenn $\chi$ eine beliebige Erhaltungsgröße ist, dann gilt
für sie Gleichung \ref{chiErhalt}. Das ist plausibel, da $I_\mathrm{coll}$ die Stoßprozesse beschreibt
und $\chi$ davon unverändert bleiben muss.

\begin{equation}
\int \chi (\bm r, \bm v ) I_\mathrm{coll} \, \d^3 v = 0
\label{chiErhalt}
\end{equation}

Um einen der Boltzmann"=Gleichung zugeordneten Erhaltungssatz zu bekommen, wird Gleichung \ref{boltzmann} mit
$\chi$ multipliziert und über $\bm v$ integriert, wegen Gleichung \ref{chiErhalt} ist die rechte Seite dann
0. Setzt man dann noch vorraus, dass $|\bm v| \rightarrow 0$ gilt und sich Mittelwerte $\left\langle A
\right\rangle$ nach Gleichung \ref{Aausf} berechnen lassen, ergibt sich Gleichung \ref{boltzErhalt}. $n(\bm
r, \bm v)$ ist darin die Teilchenzahldichte.

\begin{equation}
\pddt \left\langle n \chi \right\rangle + \pddxi \left\langle n v_i \chi \right\rangle 
- n \left\langle v_i \pddxi \chi \right\rangle - \frac{n}{\mt} \left\langle F_i \pddvi \chi \right\rangle 
 - \frac{n}{\mt} \left\langle \chi \pddvi F_i \right\rangle =0
\label{boltzErhalt}
\end{equation}

Setzt man nun für $\chi$ Masse, Impuls und Energie wie in Gleichung \ref{chiMIE} ein, führt einige kleinere
Umformungen durch und nimmt die in Gleichung \ref{subst1} bis \ref{subst2} definierten Substitutionen durch,
ergeben sich drei Erhaltungssätze (Gleichung \ref{erhaltS1} bis \ref{erhaltS2}), die, sobald $f$ bekannt ist,
sofort verwendbare hydrodynamische Gleichungen liefern.

\begin{equation}
\chi = \left( \mt, \mt v_i, \frac{1}{2} \mt |\bm v - \bm u(\bm r, t)|^2 \right)
\label{chiMIE} 
\end{equation}

\begin{align}
\label{subst1}
\rho(\bm r,t) &\equiv \mt \int f(\bm r,\bm v ,t)\,\d^3v \quad\text{Massendichte}\\
\bm u(\bm r,t) &\equiv \left\langle \bm v \right\rangle \quad \text{Mittlere Geschwindigkeit} \\
\theta(\bm r,t) &\equiv \frac{1}{3} \mt \left\langle |\bm v - \bm u|^2 \right\rangle \equiv kT
\quad\text{Temperatur}\\
\bm q(\bm r,t) &\equiv \frac{1}{2} \mt \rho \left\langle (\bm u - \bm v) |\bm v - \bm u|^2 \right\rangle
\quad\text{Wärmestrom} \\
P_{i,j} &\equiv \rho \left\langle (v_i - u_i)(v_j - u_j) \right\rangle \quad\text{Drucktensor} \\
\varLambda_{i,j} &\equiv \frac{1}{2} \mt \left( \frac{\partial u_i}{\partial x_j} + \frac{\partial
u_j}{\partial x_i} \right)
\label{subst2}
\end{align}

\begin{align}
\label{erhaltS1}
\pddt \rho + \nabla (\rho \bm u) &= 0 \quad\text{Massenerhaltung}\\
\rho\left( \pddt + \bm u \nabla \right) \bm u &= \frac{\rho}{\mt}\bm F - \nabla \mathbf P \quad
\text{Impulserhaltung} \\
\rho \left( \pddt + \bm u \nabla \right) \theta &= -\frac{2}{3} \nabla \bm q - \frac{2}{3} \sum_{i,j} P_{i,j}
\varLambda_{i,j}
\quad\text{Energieerhaltung}
\label{erhaltS2}
\end{align}

Was nun immer noch fehlt ist die Kenntnis von $f$. In den folgenden Abschnitten werden die resultierenden
hydrodynamischen Gleichungen vorgestellt, die sich aus immer besseren Näherungen für die exakte Lösung
$f$ der Boltzmann"=Gleichung ergeben.

\subsection{Transporterscheinungen in Näherung nullter Ordnung}

Die einfachste näherungsweise Lösung der Boltzmann"=Gleichung für ein System, dass sich nicht im Gleichgewicht
befindet, ist die Gleichgewichtsverteilung $f_0$. Die Maxwell"=Verteilung wurde schon in diesem Kapitel schon
bestimmt und kann deshalb direkt in die Gleichungen \ref{subst1} is \ref{subst2} eingesetzt werden. Dadurch
ergeben sich die Transportgleichungen für ein Nicht"=Gleichgewichtssystem in nullter Näherung zu Gleichung
\ref{euler1} bis \ref{euler2}. Diese nullte Näherung sind die Euler"=Gleichungen.

\begin{align}
\label{euler1}
\pddt \rho + \nabla (\rho \bm u) &= 0 \\
\left( \pddt + \bm u \nabla \right) \bm u + \frac{1}{\rho} \nabla P &= \frac{\bm F}{\mt} \\
\left( \pddt + \bm u \nabla \right) \theta + \frac{k}{c_V \mt} (\nabla \bm u )\theta &=0
\label{euler2}
\end{align}

\subsection{Transporterscheinungen in Näherung erster Ordnung}

Um die hydrodynamische Gleichungen in erster Näherung zu finden, nehmen wir an, die Lösung der
Boltzmann"=Gleichung würde sich nur durch eine geringe Störung $g$ der Verteilungsfunktion des Gleichgewichts
$f_0$ unterscheiden. Es soll also gelten $f = f_0 +g$. Wenn man diesen Ansatz in $\I$, also Gleichung
\ref{Icoll} einsetzt, bei der Berechnung dann quadratische Terme vernachlässigt und auch von den linearen
Termen nur einen berechnet, kommt man auf die Abschätzung für $\I$ in Gleichung \ref{boltzO1}.
\begin{equation}
\left( \pddt + \bm v \nabla + \frac{\bm F}{\mt} \nabla_{\bm v} \right) (f_0 + g) = I_\mathrm{coll}(f_0 + g)
\approx -\frac{g}{\tau}
\label{boltzO1}
\end{equation}
Weil $g \ll f_0$ kann man $g$ auf der rechten Seite von Gleichung \ref{boltzO1} auch gleich wieder
vernachlässigen. Die gesuchte Störung $g$ der Gleichgewichtsverteilung wird dann zu Gleichung \ref{g}.
\begin{equation}
g = -\tau\left( \pddt + \bm v \nabla + \frac{\bm F}{\mt} \nabla_{\bm v} \right) f_0
\label{g}
\end{equation}

Stellt man $f_0$ als eine Funktion von $\theta$, $\rho$ und $\bm u$ dar und formt die Ableitungen in
Gleichung \ref{g}
entsprechend um, erhält man Gleichung $\ref{g2}$. Die darin vorkommenden Abkürzungen $\D$ und $\bm U$ sind in
Gleichung \ref{subst3} definiert.
\begin{equation}
g = -\tau f_0 \left( \frac{1}{\rho} \D \rho + \frac{1}{\theta}\left( \frac{\mt}{2\theta} U^2 - \frac{3}{2}
 \right) \D \theta + \frac{\mt}{\theta} U_j \D u_j - \frac{1}{\theta} \bm F \bm U \right) 
\label{g2}
\end{equation}

\begin{equation}
\D X \equiv \left( \pddt + v_i \pddxi \right) X \quad \bm U \equiv \bm v - \bm u  
\label{subst3}
\end{equation}

Mit Hilfe der Euler"=Gleichungen \ref{euler1} bis \ref{euler2} kann man explizite Ausdrücke für $\D \rho$, $\D
\theta$ und $\D u_j$ finden. Setzt man diese in Gleichung \ref{g2} ein, ergibt sich mit der Substitution $P =
\rho \theta /\mt$ sowie weiteren Umformungen und Vereinfachungen

\begin{equation}
g = -\tau \left(  \frac{1}{\theta}\pddxi \theta U_i \left( \frac{\mt}{2\theta} U^2 -
\frac{5}{2} \right) + \frac{1}{\theta} \varLambda_{i,j} (U_i U_j - \frac{1}{3} \updelta_{i,j} U^2) \right)
f_0.
\label{g3}
\end{equation}

Setzt man $g$ aus \ref{g3} in die Erhaltungsgleichungen \ref{erhaltS1} bis \ref{erhaltS2} ein, erhält man für 
den Wärmestrom $\bm q$ Gleichung \ref{q} und für den Drucktensor $\mathbf P$ Gleichung \ref{P}.

\begin{equation}
\bm q = -K \nabla \theta \quad \text{mit} \quad K = \frac{\mt^2 \tau}{6\theta} \int U^4 \left(
\frac{\mt}{2\theta} U^2 - \frac{5}{2} \right) f_0 \, \d^3 U = \frac{5}{2} \tau \theta n
\label{q}
\end{equation}

\begin{align}
\label{P}
P_{i,j} = \updelta_{i,j} P- \frac{2\mu}{\mt} \left( \varLambda_{i,j} - \frac{\mt}{3} \updelta_{i,j} \nabla 
\bm u \right) \\ \text{mit} \quad \mu = \frac{\mt^2 \tau}{\theta} \int U^2 (\bm v_1) U^2 (\bm v_2) f_0 \,
\d^3 U =  \tau \theta n
\end{align}

Die Größen $\mu$ und $K$ sind die Viskosität und die Wärmleitfähigkeit. Der Quotient der beiden $K/\mu=5/2$
ist ein Konstante, die mit den Freiheitsgraden des Gases zusammenhängt. Falls es sich nicht um ein
ein"=atomiges
Gas handelt, ist der Zusammenhang zwischen den beiden komplizierter. Er ist in Gleichung \ref{wlT}
dargestellt.
Das in Gleichung \ref{wlT} definierte $\kappa$ hat nicht die gleiche Einheit wie $\mu$. Die Einheit
unterscheidet sich um den Faktor $\mathrm{J/(kg\,K)}$, da $\kappa$ in Bezug auf $T$ und nicht in Bezug auf
$\theta = kT$ und außerdem pro Masse definiert wird.

Wenn man den am Anfang diese Kapitels vorgestellten Ausdruck für die mittlere Stoßzeit $\tau$ in die Formeln
für $\mu$ und $K$ aus Gleichung \ref{q} und \ref{P} einsetzt, ergibt sich die Beziehung in Gleichung
\ref{muK}. Darin zeigt sich, dass $\mu$ und $K$ beziehungsweise $\kappa$ nach auf der statistischen
Physik beruhenden Überlegungen von der Temperatur $T$, der Teilchenmasse $\mt$ und vom Teilchendurchmesser
$\sigma$ abhängen.

\begin{equation}
\mu \sim K \sim \frac{\sqrt{\mt kT}}{\sigma^2}
\label{muK}
\end{equation}

Um nun die gesuchten hydrodynamischen Gleichungen in erster Näherung zu erhalten, müssen noch $\bm q$ und
$\mathbf P$ aus Gleichung \ref{q} und \ref{P} in die allgemeinen Erhaltungsgleichungen \ref{erhaltS1} bis
\ref{erhaltS2} eingesetzt werden. Nach einigen Umformungen sowie der Vernachlässigung aller Terme, die
Ableitungen von $\mu$ oder $K$ enthalten oder quadratisch in $\bm u$ sind, ergeben sich Gleichung \ref{Aausf}
bis \ref{Aausf}. Das sind die Navier"=Stokes"=Gleichungen. Darin gilt $\bm F/\mt = \bm a_\mathrm{ext}$
mit einer externen Beschleunigung $\bm a_\mathrm{ext}$, weil $\bm F$ eine externe Kraft pro Teilchen ist.

\begin{align}
\label{NSstat1}
\pddt \rho + \nabla (\rho \bm u) &= 0 \\
\left( \pddt \bm u + \bm u \nabla \bm u \right)  &= - \frac{1}{\rho} \nabla \left( P-\frac{\mu}{3} \nabla
\bm u \right) + \frac{\mu}{\rho} \nabla^2 \bm u + \frac{\bm F}{\mt} \\
c_V \frac{\mt}{k} \left( \pddt + \bm u \nabla \right) \theta +  (\nabla \bm u )\theta &= \frac{K}{\rho}
\nabla^2 \theta
\label{NSstat2}
\end{align}


\subsection{Transporterscheinungen in Näherung höherer Ordnung}
\label{ch-en}

In den letzten beiden Abschnitten wurden die hydrodynamischen Gleichungen als Näherungen der
Boltzmann"=Gleichung in nullter und erster Ordnung hergeleitet. Ein systematisches Verfahren, Näherungen der
Verteilungsfunktion $f$ von immer höherer Ordnung $f_k$ zu erzeugen, ist die Chapman"=Enskog"=Methode. Sobald
ein $f_k$ bekannt ist, ergeben sich dann die Transportgleichungen durch Einsetzten von $f_k$ in allgemeinen
Erhaltungsgleichungen \ref{erhaltS1} bis \ref{erhaltS2}.

Zur Erinnerung noch einmal die Definition der Boltzmann"=Gleichung \ref{boltzmann2} mit dem Kollisionsterm
$\I$ in Gleichung \ref{Icoll2}.
\begin{equation}
\pddt f + \D f = \I(f,f) \quad \text{mit} \quad \D f = \left( \bm v \nabla + \frac{\bm F}{\mt} \nabla_{\bm v}
\right) f
\label{boltzmann2}
\end{equation}
\begin{equation}
\I(f,g) = \iint \sigma(\varOmega) \left| \bm v - \bm v_1 \right| \left( f( \bm v')
g( \bm v_1') - f( \bm v)g( \bm v_1)\right) \, \d^3 v_1 \, \d \varOmega
\label{Icoll2}
\end{equation}
In einer ersten Einschränkung soll $f$ nur implizit über die Teilchendichte $n$, die Fluid"=Geschwindigkeit
$\bm u$ und die Temperatur $\theta$ von der Zeit $t$ abhängen. Formal wird für $f$ ein Reihenansatz wie in
Gleichung \ref{fAnsatz} gewählt, in dem der Beitrag von $f_k$ zu $f$ immer kleiner wird mit steigendem
$k$.

\begin{equation}
\label{fAnsatz}
f = \frac{1}{\zeta} \sum_{k=0}^\infty \zeta^k f_k
\end{equation}

Des Weiteren wird gefordert, dass $n$, $\bm u$ und $\theta$ vollständig durch $f_0$ bestimmt sind
beziehungsweise
Gleichungen \ref{nut1} bis \ref{nut2} gelten. Das hat die wichtige Konsequenz, dass Masse, Energie und
Impuls vollständig in der nullten Näherung enthalten sind und nicht zum Beispiel ein Teil der Teilchen -- also
der Masse -- durch weglassen der höheren Näherungsglieder verschwindet. Diese Einschränkung ist der einfachste
Weg, dass die fundamentalen Erhaltungssätze gelten, egal bei welchem Glied die Reihe abgebrochen wird.


\begin{align}
\label{nut1}
\int f_0 \,\d^3 v &= n \\
\frac{1}{n} \int f_0 \bm u \,\d^3 v &= \bm u \\
\frac{1}{n} \int f_0 \frac{\mt}{3} | \bm v - \bm u |^2 \,\d^3 v &= \theta \\
\int f_k \cdot (1, \bm v, v^2 ) \,\d^3 v &= (0,0,0) \quad \forall k > 0
\label{nut2}
\end{align}

Die Erhaltungssätze \ref{erhaltS1} bis \ref{erhaltS2} werden mit dem Reihenansatz von $f$ zu Gleichung
\ref{erhaltR1} bis \ref{erhaltR2}. Die darin vorkommenden Größen $\mathbf P_k$ und $\bm q_k$ sind in Gleichung
\ref{Pk} und \ref{qk} definiert.

\begin{align}
\label{erhaltR1}
\pddt \rho + \sum_i \pddxi(\rho u_i) &= 0 \\
\rho \left( \pddt u_i + \sum_j u_j \frac{\partial}{\partial x_j} u_i \right) &= \frac{\rho}{\mt} F_i -
\sum_{k=0}^\infty \zeta^k \sum_j \frac{\partial P_{i,j,k}}{\partial x_j} \\
\rho \left( \pddt \theta + \sum_j u_j \frac{\partial \theta}{x_j} \right) &= -\frac{2}{3} \sum_{k=0}^\infty
\zeta^k \left( \sum_i \pddxi q_{i,k} + \sum_{i,j} \varLambda_{i,j} P_{i,j,k} \right) 
\label{erhaltR2}
\end{align}

\begin{align}
\label{Pk}
P_{i,j,k} &= \mt \int f_k (v_i- u_i)(v_j- u_j) \,\d^3 v \\
q_{i,k} &= \frac{\mt^2}{2} \int f_k (v_i- u_i) |\bm v - \bm u|^2 \,\d^3 v
\label{qk}
\end{align}

Als nächstes werden auch die Differenzialoperatoren in der Boltzmann"=Gleichung durch Reihenentwicklungen
dargestellt. Zunächst in Gleichung \ref{DR} der Operator $\D$ aus Gleichung \ref{boltzmann2}.
\begin{equation}
\D f = \frac{1}{\zeta} \sum_{k=0}^\infty \zeta^k \D f_k
\label{DR}
\end{equation}
Weil $f$ nur über $\rho = \mt n$, $\theta$ und $\bm u$ von der Zeit abhängt, kann man $\pddt f$ umschreiben
zu
\begin{equation}
\pddt f = \frac{\partial f}{\partial \rho} \pddt \rho + \frac{\partial f}{\partial u_i} \pddt u_i
+ \frac{\partial f}{\partial \theta} \pddt \theta.
\label{pddtf}
\end{equation}
Die Entwicklung der Ableitungen $\frac{\partial f}{\partial X}$ mit $X \in \{ \rho, \theta, u_i \}$ lautet
\begin{equation}
\frac{\partial f}{\partial X} = \frac{1}{\zeta} \sum_{k=0}^\infty \zeta^k \frac{\partial f_k}{\partial X}.
\label{pdfdXR}
\end{equation}
Es fehlen noch die Entwicklungen der Ableitungen der Art $\pddt X$. Diese werden aus der $k$ten"=Näherung der
Erhaltungssätze bestimmt und der Tatsache, dass diese schon für $k = 0$ erfüllt sein müssen, bestimmt. Sie
werden in Gleichung \ref{pddtR1} bis \ref{pddtR2} definiert.

\begin{align}
\label{pddtR1}
\frac{\partial_0}{\partial t} \rho &= - \sum_i \pddxi(\rho u_i) \\ \frac{\partial_k}{\partial t} \rho &= 0
\quad
\forall k > 0 \\
\frac{\partial_0}{\partial t} u_i &= - \sum_j \frac{\partial u_i}{\partial x_j} + \frac{F_i}{\mt} -
\frac{1}{\rho} \sum_j \frac{\partial P_{i,j,0}}{\partial x_j} \\ \frac{\partial_k}{\partial t} u_i
&= -\frac{1}{\rho} \sum_j \frac{\partial P_{i,j,k}}{\partial x_j} \quad \forall k > 0 \\
\frac{\partial_0}{\partial t} \theta &= - \sum_i u_i \pddxi \theta - \frac{2}{3} \frac{1}{\rho} \sum_i \pddxi
q_{i,0} - \frac{2}{3} \frac{1}{\rho} \sum_{i,j} \varLambda_{i,j} P_{i,j,0} \\
\frac{\partial_k}{\partial t} \theta &= - \frac{2}{3} \frac{1}{\rho} \left( \sum_i \pddxi
q_{i,k} + \sum_{i,j} \varLambda_{i,j} P_{i,j,k} \right) \quad \forall k > 0
\label{pddtR2}
\end{align}

Gleichung \ref{IcollR} ist, analog zu den bisherigen Entwicklungen, der Reihenansatz des Kollisionsterms. Das
$k$-te Glied von $\I$ wird über Gleichung \ref{IcollnR} definiert.
\begin{align}
\I(f,f) = \frac{1}{\zeta^2} \sum_{k,l=0}^\infty \zeta^{k+l} \I(f_k,f_l)
\label{IcollR}
\end{align}
\begin{align}
I_{\mathrm{coll},k}(f_0,f_1,\dots, f_k) = \sum_{r,s \text{ mit } r+s = k} \I(f_r,f_s)
\label{IcollnR}
\end{align}
Wenn $\I$ über Gleichung \ref{IcollnR} ausgedrückt wird, erhält man
\begin{equation}
\I(f,f) = \frac{1}{\zeta^2} \sum_{k=0}^\infty \zeta^{k} I_{\mathrm{coll},k}.
\label{IcollRn}
\end{equation}

Nach all diesen Definitionen kann nun die gesamte Reihenentwicklung der Boltzmann"=Gleichung
geschrieben werden als

\begin{multline}
\frac{1}{\zeta}\left( \left( \frac{\partial_0}{\partial t} + \zeta\frac{\partial_1}{\partial t} + \zeta
\frac{\partial_2}{\partial t} + \dots \right) + \D\right) ( f_0 + \zeta f_1 + \zeta^2 f_2 + \dots ) \\
= \frac{1}{\zeta^2} \left( I_{\mathrm{coll},0}(f_0) + \zeta I_{\mathrm{coll},1}(f_0,f_1) + \zeta^2
I_{\mathrm{coll},2}(f_0,f_1,f_2) + \dots \right).
\label{boltzmannR}
\end{multline}

Die $k$"=te Lösung $f_k$ kann nun gefunden werden indem gefordert wird, dass $\zeta = 1$. Es ergibt sich dann
die
Folge von rekursiv lösbaren Näherungen der Boltzmann"=Gleichung. Die Gleichungen \ref{boltzmannR1} bis
\ref{boltzmannR2} zeigen die ersten drei Glieder dieser Folge.

\begin{align}
\label{boltzmannR1}
0 &= I_{\mathrm{coll},0}(f_0) \\
\frac{\partial_0}{\partial t} f_0 + \D f_0 &=  I_{\mathrm{coll},1}(f_0,f_1) \\
\frac{\partial_0}{\partial t} f_1 + \frac{\partial_1}{\partial t} f_0 + \D f_0 &= 
I_{\mathrm{coll},2}(f_0,f_1,f_2)
\label{boltzmannR2}
\end{align}

Die Lösung von Gleichung \ref{boltzmannR1} ist natürlich die Maxwell"=Verteilung $f_0$ und die Entsprechenden
Transportgleichungen sind die Euler"=Gleichungen. Mit dem bekannten $f_0$ kann die Korrektur erster Ordnung
$f_1$ berechnet werden. Die dazugehörigen Transportgleichungen sind die Navier"=Stokes"=Gleichungen. Mit $f_0$
und $f_1$ kann dann die Korrektur zweiter Ordnung errechnet werden, was auf die Burnett"=Gleichungen führt.

\section{Randbedingungen}

Die gewöhnlichen Geschwindigkeitsrandbedingungen für hydrodynamische Gleichungen sind $\bm u(\varGamma) = \bm
u_\varGamma$, wobei $\varGamma$ der Rand des Definitionsgebietes ist. Das heißt ein Fluid hat an der Wand die
Geschwindigkeit der Wand. Im Fall eine ruhenden Wand
gilt natürlich $\bm u_\varGamma = 0$. Diese Annahme ist in gutem Einklang mit experimentellen Daten.

Es sind zwei Extreme vorstellbar, in denen $\bm u(\varGamma) = 0$ kein gutes Modell der Realität ist:
Der eine Fall tritt ein, wenn die interne Viskosität einer Flüssigkeit viel größer ist als die Reibung mit der
Wand. So kann
man bei einem ein Gel, das aus einer Tube gedrückt wird, erwarten, dass das Gel über seine
Begrenzungen
gleitet. Der andere, wenn ein Gas stark verdünnt ist. Dann kann sich ein Gasteilchen sehr lange parallel zur
Wand bewegen, ohne dass es mit anderen Teilchen, die die Information der nahen Wand vermitteln könnten, in
Kontakt kommt. In dieser Arbeit wird nur auf die letztere Möglichkeit eingegangen.

Knudt und Warburg haben schon 1875 aus Messungen geschlossen, dass es in verdünnten Gasen zu einem Gleiten
über Wände kommt. Die Existenz dieses Gleitens oder "`slip"' wurde experimentell für viele verschiedene
Gas"=Wand"=Kombinationen nachgewiesen. Eine umfangreiche Auflistung der Resultate von Gleit"=Messungen an
verschiedenen Flüssigkeiten bietet \cite{lauga-2007}. Eine Messung der Gleitlänge von Luft über Glas
wurde von \cite{maali:027302} durchgeführt.

Die folgenden Herleitung einer Gleitrandbedingung mit thermischer Transpiration für Gase geht auf Maxwell
\cite{maxwell} zurück, ist aber in der hier stark vereinfachten Form \cite{kennard} entnommen.

\subsubsection{Slip-Randbedingung}

Zunächst wird die Annahme gemacht, dass ein Gasmolekül mit der Wand entweder eine einfache Reflexion
ausführt oder von ihr absorbiert und mit unter einen gleich"=verteilten Winkel reemittiert wird; es muss also
ein
Impulsübertrag stattfinden. Der Anteil der Moleküle, der absorbiert und reemittiert wird, wird mit $\sigma_p$
bezeichnet.

Sei $u_y(x)$ die Geschwindigkeit mit der ein Gas eine Wand parallel zur y"=Achse vorbeiströmt, wobei es zu
einem Gleiten entlang dieser Wand kommt. Es soll nun ein Ausdruck gefunden werden, der den Übertrag der
y"=Komponente des Impulses pro Zeit pro Fläche auf die Wand beschreibt. In einem sich mit $u_y$ mitbewegenden
System wird $\frac{1}{2} \mu \frac{\partial u_y}{\partial x}$ auf die Wand pro Zeit pro Fläche übertragen.
Dazu kommt noch der von dem Gleiten auf der Wand $u_y(\varGamma)$ übertragene Impuls. Dieser ist $\mt
u_y(\varGamma) \cdot \frac{1}{4} n\bar v$, denn $\frac{1}{4} n\bar v$ sind die Einschläge auf eine
Seite einer beliebigen
Fläche im Gas pro Zeit pro Fläche. In Gleichung \ref{slip1} werden diese beiden Anteile addiert und mit
$\sigma_p$
multipliziert, da nur dieser Anteil der Teilchen nach der Übereinkunft im letzten Abschnitt Impuls auf die
Wand übertragen kann.

\begin{equation}
\sigma_p \left( \frac{1}{2}\mu \frac{\partial u_y}{\partial x} + \mt
u_y(\varGamma) \frac{1}{4} n\bar v  \right) = \mu \frac{\partial u_y}{\partial x}
\label{slip1}
\end{equation}

Aufgrund der Impulserhaltung muss die rechte Seite von \ref{slip1} gerade $\mu \frac{\partial u_y}{\partial
x}$ sein, denn die Scherspannung im Gas ist gerade der tangentiale Impuls pro Fläche pro Zeit, der vom Gas in
Richtung Wand transportiert wird. Löst man Gleichung \ref{slip1} nach $u_y(\varGamma)$ auf, ersetzt $\bar v$
mit Gleichung \ref{mg} und $n\mt$ mit $\frac{P}{\Rs T}$ (ideale Gasgleichung), ergibt sich

\begin{equation}
u_y(\varGamma) = \frac{2-\sigma_p}{\sigma_p} \frac{\sqrt{\uppi \Rs T} \mu}{ \sqrt 2 P} \frac{\partial
u_y}{\partial x}.
\label{slip2}
\end{equation}

Kennard benutzt den Ausdruck $\mu = c\rho\bar v \lambda$ für die Viskosität, wobei $c$ ein Koeffizient ist,
der aus einer auf Enskok zurückgehende Reihenentwicklung bestimmt wird und je nach Genauigkeit der Entwicklung
zwischen 0.491 und 0.499 liegt. Wenn $c=5\uppi/32$, dann ist Kennards Modell der Viskosität identisch
mit dem in dieser Diplomarbeit verwendeten und in Gleichung \ref{viskT} dargestelltem Modell.

Wird nun Kennards Viskositätsmodell in Gleichung \ref{slip2} eingesetzt und der dann darin
auftauchende Faktor $2c$ mit 1 genähert, ergibt sich Gleichung \ref{slip}, welche die Randbedingung für
Gleitgeschwindigkeiten nach Maxwell ist.

\begin{equation}
u_y(\varGamma) = \frac{2-\sigma_p}{\sigma_p} \lambda \frac{\partial u_y}{\partial x}
\label{slip}
\end{equation}


\subsection{Thermische Transpiration}

\subsubsection{Qualitative Beschreibung}

Zunächst widerspricht eine Gasströmung von einem kalten in einen warmen Bereich, ohne dass eine
Druckdifferenz vorliegt der Intuition. Die folgende
Abschätzung soll es
plausibel machen, dass es so etwas wie thermische Transpiration geben muss.

Man stelle sich zwei Container verbunden mit einem Kanal vor. In den Containern herrschen unterschiedliche
Temperaturen. Wenn der Kanal sehr dick ist, ist der offensichtliche Gleichgewichtszustand dieses Systems,
dass der Druck überall gleich ist, die unterschiedlichen Temperaturen aber wegen des idealen Gasgesetzes
durch entsprechende Dichten ausgeglichen werden. Es gilt also $T_1/T_2 = \rho_2/\rho_1$. Nun soll der Kanal
enger sein und das Gas so dünn, dass es nicht mehr mit sich selbst im Kanal wechselwirkt. Für die Dichte
$\rho$ in jedem Container gilt $\rho \sim n$ und für die mittleren Teilchengeschwindigkeiten gilt $ v^2 \sim
T$. Für den Massenfluss von Container 1 nach Container 2 gilt: $\sim \mt n_1 v_1$. Und entsprechend für den
gegenläufigen Massenfluss: $\sim \mt n_2 v_2$. Wenn sich ein Gleichgewicht unter der Aufrechterhaltung der
unterschiedlichen Temperaturen eingestellt hat, muss Gleichung \ref{PPTT} für die beiden Massenflüsse gelten.

\begin{equation}
\mt n_1 v_1 = \mt n_2 v_2 \Longleftrightarrow \frac{\mt n_1 v_1}{ \mt n_2 v_2} = 1
\label{PPTT}
\end{equation}

Mit denen im letzten Abschnitt genannten Proportionalitäten kann man Gleichung \ref{PPTT} umformen zu

\begin{equation}
1 = \frac{\mt n_1 v_1}{\mt n_2 v_2} = \frac{\rho_1}{\rho_2} \sqrt{\frac{T_1}{T_2}} = \frac{P_2}{P_1}
\sqrt{\frac{T_1}{T_2}} \Longleftrightarrow \frac{P_1}{P_2} = \sqrt{\frac{T_1}{T_2}}.
\label{PPTT2}
\end{equation}

Gleichung \ref{PPTT2} macht deutlich, dass der Temperaturunterschied in den beiden Containern im
Gleichgewicht zu einem Druckunterschied in eben diesen führen muss. In der Stärke, wie ihn Gleichung
\ref{PPTT2} vorhersagt, ist er allerdings nur für den Fall $\kn \rightarrow \infty$. Bei endlichen
$\kn$"=Werten
ist das Verhältnis der Drücke kleiner. Bei $\kn = 0$ gilt dann $P_1 = P_2$.

\begin{figure}
\capstart\centering
\includegraphics[width=0.55\linewidth]{fig/whatyouneedDE}
\caption{Schematische Darstellung von Teilchen, die mit verschiedenen Geschwindigkeiten (schwarze Pfeile) auf
das Element $a$ einer Wand mit Temperaturgradient (rot) eintreffen. }
\label{fig:whatyouneed}
\end{figure}

Als nächstes folgt eine qualitative Beschreibung der Vorgänge auf Teilchenebene im Kanal. Es soll weiterhin
die Annahme gelten, dass $\lambda$ so groß ist, dass Wechselwirkungen der Teilchen miteinander gegenüber denen
mit der Wand zu vernachlässigen sind. Der Kanal sei parallel zur x"=Achse wie es in Abbildung
\ref{fig:whatyouneed} skizziert ist. Entlang x"=Achse sei ein positiver Temperaturgradient. Betrachtet man die
Teilchen, die auf eine Stelle $a$ auf der Wand eintreffen, so ist ihr Impuls, allein
von der Temperatur, der Stelle von der sie stammen, bestimmt, weil sie eben unterwegs nicht mit anderen
Teilchen zusammengestoßen sind. Wenn ein Teilchen mit der Wand in Wechselwirkung tritt, sei diese
diffus, das heißt das
Teilchen wird von der Wand absorbiert und reemittiert und zwar mit einer Geschwindigkeit, die nur von der
Temperatur an der Stelle $a$ des Wandstoßens abhängt. Des Weiteren soll die Richtung, mit der sie die Wand
verlassen, gleichverteilt sein. Also verlassen die Teilchen die Wand im Mittel genau senkrecht.

Teilchen die von negativer x"=Richtung kommen, tragen weniger Impuls als Teilchen von positiver x"=Richtung
wegen der unterschiedlichen Temperaturen, die in den Regionen herrschen aus denen sie kommen. Im Durchschnitt
verlassen die Teilchen die Stelle $a$ aber mit einer tangentialen Impulskomponente von $\bar p_x = 0$. Es muss
also im Schnitt tangentialer Impuls in negativer x"=Richtung auf Wand übertragen worden sein. Da die Wand
festgehalten wird, sorgt die Impulserhaltung dafür, dass das Gas einen Impuls in positive x"=Richtung erfährt.
Es gibt also, aufgrund des Temperaturgradienten, eine Gasströmung in Richtung der höheren Temperatur.



\subsubsection{Maxwells Randbedingung}


Ein Gas, in dem ein Temperaturgradient in y"=Richtung herrscht, soll durch die Maxwellverteilung $f_0$ und ein
Korrekturterm erster Ordnung $f_1$, der die Abweichung vom thermischen Gleichgewicht enthält, beschrieben
werden. Die Gestalt von $f_0$ und $f_1$ ist
\begin{equation}
f_0 = A\e^{-\beta^2 \bm v^2} \quad \text{und} \quad f_1 = C v_x \left( \frac{5}{2} - \beta^2 \bm v^2 \right)
\e^{-\beta^2 \bm
v^2}.
\label{tt1}
\end{equation}
Wobei $A = \beta^3 \uppi^{-3/2}$ und $\beta^2 = \mt/(2kT)$ ist und $C$ in Gleichung \ref{defC} definiert ist.
$C$ ergibt sich aus längerer Rechnung durch den Vergleich der durch $f_0$ und $f_1$ hervorgerufenen
Kollisionsraten und enthält die Physik des Wärmetransports.
\begin{equation}
C=\frac{3}{4}\frac{\beta^8}{\uppi \sqrt 2} \frac{1}{n^2 \uppi \sigma} \cdot \frac{5}{4}
\frac{\uppi^{\frac{3}{2}}}{\beta^7} \frac{nA}{T} \frac{\partial T}{\partial s}
\label{defC}
\end{equation}

Mit den Integralen in \ref{tt2} ist es möglich zu Prüfen, ob die Korrektur erster Ordnung zu einem netto
Impulsfluss im Gas
kommt. Dabei stellt man fest, dass es zu keiner Bewegung des Gases kommt.

\begin{equation}
\left\langle p_x\right\rangle  \sim \int_\infty^\infty v_x^2 ( f_0 + f_1 ) \,\d^3 v = 0 \quad \left\langle
p_y\right\rangle  \sim \int_\infty^\infty v_x v_y ( f_0 + f_1 ) \,\d^3 v = 0
\label{tt2}
\end{equation}

Ganz anderes sieht die Situation aus, wenn man ein Gas in der Nähe einer Wand betrachtet. Die Wand soll
parallel zur y"=Achse sein und auf ihr ein Temperaturgradient herrschen, der sich auf das Gas übertragen hat.
Das Gas soll wieder durch $f_0 + f_1$ beschrieben werden. Weil wir aber nun alle Teilchen betrachten, die sich
unter jedem möglichen Winkel $\alpha$ wie in Abbildung \ref{fig:ttskizze} auf die Wand zu bewegen, wird $v_x$
in $f_1$ durch $v_x \cos(\alpha) +
v_y \sin(\alpha)$ ersetzt. Ausdruck \ref{tt3} zeigt die Impulskomponente in y"=Richtung
pro Zeit pro Flächenelement, die durch Teilchen, die eine
Geschwindigkeitskomponente in Richtung der Wand ($-v_x$) haben, verursacht wird.

\begin{figure}
\capstart\centering
\includegraphics[width=0.3\linewidth]{fig/ttskizze}
\caption{Teilchen, die unter dem Winkel $\alpha$ auf eine Stelle einer Wand auftreffen, und dabei den
Temperaturgradienten $\frac{\partial T}{\partial s}$ durchqueren, der eine Projektion des eigentlichen
Gradienten auf der Wand parallel zur y-Achse ist. }
\label{fig:ttskizze}
\end{figure}

\begin{equation}
-\mt n\int\limits^\infty_{-\infty} \int\limits^\infty_{-\infty} \int\limits^0_{-\infty} v_x v_y \left(
A+C ( v_x \cos(\alpha) + v_y \sin(\alpha) ) \left( \frac{5}{2} - \beta^2 \bm v^2 \right)   \right)
\e^{-\beta^2 \bm v^2} \,\d v_x \d v_y \d v_z
\label{tt3}
\end{equation}

Wertet man das Integral in Ausdruck \ref{tt3} aus, ergibt sich Ausdruck \ref{tt4}. Setzt man darin $C$
aus Gleichung \ref{defC} ein, nachdem man darin vorher $\sigma$ mit Gleichung \ref{viskT} ersetzt hat, ergibt
sich Ausdruck \ref{tt5}.

\begin{equation}
-\frac{\uppi}{8}\frac{n\mt C}{\beta^6}\sin \alpha
\label{tt4}
\end{equation}

\begin{equation}
-\frac{3}{4\uppi}\frac{\mu}{\beta^2 T}\sqrt\frac{\uppi \mt}{8kT}\frac{\partial T}{\partial s} \sin \alpha
\label{tt5}
\end{equation}

Ausdruck \ref{tt5} ist der Impulsfluss (y"=Komponente des Impulses pro Zeit pro Fläche), der auf die Wand
übertragen wird. Das stetige Entziehen von Impuls des Gases kann nur durch eine makroskopische Bewegung $u_y$
des Gases ausgeglichen werden, da weiter weg von der Wand ja kein netto Impulstransport stattfindet. Die
y"=Komponente des Impulses pro Zeit pro Fläche, die ein Gasfluss in Richtung des Temperaturgradienten
aufnehmen kann, ist $\frac{1}{4} \rho \bar v u_y$; dies muss mit Ausdruck \ref{tt5} gleichgesetzt werden, was
Gleichung \ref{tt6} ergibt.

\begin{equation}
-\frac{1}{4} \rho \bar v u_y=-\frac{3}{4\uppi}\frac{\mu}{\beta^2 T}\sqrt\frac{\uppi \mt}{8kT}\frac{\partial
T}{\partial s} \sin \alpha
\label{tt6}
\end{equation}

Löst man Gleichung \ref{tt6} nach $u_y$ auf, ergibt sich Gleichung \ref{tt7}. $\beta$ und $\bar v$ wurde
durch ihre in diesem Kapitel eingeführten Definitionen ersetzt und $\frac{\partial T}{\partial s}\sin \alpha$
wurden durch $\frac{\partial T}{\partial y}$ ersetzt, weil $s$ der Weg auf ein Element der Wand unter allen
Winkeln $\alpha$ ist.

\begin{equation}
u_y = \frac{3}{4} \frac{\mu}{\rho T} \frac{\partial T}{\partial y}
\label{tt7}
\end{equation}

Gleichung \ref{tt7} ist der thermische Transpirationsterm von Maxwells Randbedingungen. Er muss zum
Gleitterm aus Gleichung \ref{slip} hinzuaddiert werden um Maxwells vollständige Randbedingungen zu erhalten.
Bei beliebigen 2D Geometrien muss der Gleitterm noch um die zweite mögliche, partielle Ableitung erweitert
werden. Gleichung \ref{tt} zeigt Maxwells vollständige Randbedingung in einem tangential"=normalen
Koordinatensystem mit $t$ und $n$ als Koordinaten, wie sie zum Beispiel in \cite{PhysRevE.70.017303} zu finden
ist.

\begin{equation}
u_t(\varGamma) = \frac{2-\sigma_p}{\sigma_p} \lambda 
\left( \frac{\partial u_t}{\partial n} + \frac{\partial u_n}{\partial t} \right) 
+ \frac{3\mu}{4 \rho T} \frac{\partial T}{\partial t}
\label{tt}
\end{equation}


\subsection{Temperatursprung-Randbedingungen}

Sehr ähnliche Überlegungen, wie die bei der Herleitung der Gleitrandbedingung, kann man auch über die
Temeraturrandbedingung anstellen. Es folgt dann, dass ein Gas direkt an eine Wand nicht genau die
Temperatur der Wand annimmt. Entsprechende Randbedingungen für die Wärmeleitgleichung wurden von Smoluchowski
hergeleitet \cite{smoluchowski} und auch Experimente zum Temperatursprung zwischen Wänden und verdünnten
Gasen wurden von ihm durchgeführt. Die hier präsentierte Herleitung richtet sich nach Kennard \cite{kennard}.

Wie in den anderen Herleitungen in diesem Kapitel sollen Teilchen, die eine Geschwindigkeitskomponente in
Richtung einer Wand haben, diese Wand relativ ungestört erreichen, und ihre Energie und ihr Impuls dem Ort
entsprechen, von dem sie herkommen. Die Wechselwirkung wird, wie bei den Gleitrandbedinungen, durch einen
Koeffizienten $\sigma_T$ beschrieben, der aussagt, wie groß der Anteil der Teilchen ist, die einen
Energieaustausch mit der Wand vornehmen; also deren Temperatur annehmen. $1-\sigma_T$ ist entsprechend der
Anteil der Moleküle, die an der Wand nur einen Reflexion vornehmen und keine Energie austauschen. $\sigma_T$
wird "`thermal accomodation coefficient"' genannt, da er den Wärmeenergieaustausch beziehungsweise
Temperaturaustausch
beschreibt. $\sigma_p$ dagegen beschreibt den Austausch des tangentialen Impulses.

Wir betrachten einen Teilchenstrom bestehend aus Teilchen mit einer Geschwindigkeitskomponente in Richtung der
Wand mit der Energie pro Zeit pro Fläche $E_\mathrm i$. Die Teilchen, die von der Wand zurückkommen, ohne mit
ihr Energie ausgetauscht zu haben, besitzen die Energie pro Zeit pro Fläche $E_\mathrm r$ und die Teilchen,
die Energie ausgetauscht haben $E_\mathrm W$. Gleichung \ref{sm0} zeigt den Zusammenhang zwischen den
Energien dieser Teilchenströme.

\begin{equation}
E_\mathrm i - E_\mathrm r = \sigma_T (E_\mathrm i - E_\mathrm W )
\label{sm0}
\end{equation}


Die Gasteilchen mit einer Geschwindigkeitskomponente in Richtung der Wand sollen die Temperatur $T_\mathrm
G$
haben. Sie transportieren auf zwei Arten Energie auf die Wand: Durch die kinetische Energie und durch
Diffusion. Der Diffusionsbeitrag ist $\frac{1}{2}\kappa\frac{\partial T}{\partial n}$. $n$ ist die Koordinate
senkrecht auf der Wand.

Die Translationsenergie pro Masse des Teilchenstroms ist $2\Rs T_\mathrm G$. Diese Zahl ist $\frac{4}{3}$ mal
größer
als die spezifische innere Energie $U_\mathrm G$ eines Gases im Gleichgewicht. Ausdruck \ref{sm1} ist die
Summe aus Translationsenergie und innerer Energie des Gases, das sich in Richtung der Wand
bewegt, während Ausdruck \ref{sm2} die Energie ist, die sich von der Wand, die die Temperatur $T_\mathrm W$
hat, wegbewegt.
\begin{equation}
m( 2RT_\mathrm G + U_\mathrm G )
\label{sm1}
\end{equation}
\begin{equation}
m( 2RT_\mathrm W + U_\mathrm W )
\label{sm2}
\end{equation}
Eine gute Näherung der Differenz der beiden Terme \ref{sm1} und \ref{sm2} ist Ausdruck \ref{sm3}.
\begin{equation}
m\left(  c_V + \frac{1}{2}\Rs \right) (T_\mathrm G - T_\mathrm W )
\label{sm3}
\end{equation}

Wenn man zu Ausdruck \ref{sm3} noch die Diffusion addiert, ergibt sich Gleichung \ref{sm4}: Die Differenz
der einlaufenden Energie pro Zeit pro Fläche und der auslaufenden Energie pro Zeit pro Fläche, die mit der
Wand energetisch gewechselwirkt hat.
\begin{equation}
E_\mathrm i - E_\mathrm W = \frac{1}{2} \kappa \frac{\partial T}{\partial n} + S\left(  c_V + \frac{1}{2}\Rs
\right) (T_\mathrm G - T_\mathrm W )
\label{sm4}
\end{equation}
$S$ ist in Gleichung \ref{sm4} die Masse, die pro Zeit auf einer Fläche ankommt, also $S = \frac{1}{4} \rho
\bar v = P/\sqrt{2\uppi\Rs T}$. Setzt man das in Gleichung \ref{sm4}, zusammen mit $\Rs = (\gamma -1) c_V$
ein, ergibt sich
\begin{equation}
E_\mathrm i - E_\mathrm W = \frac{1}{2} \kappa \frac{\partial T}{\partial n} + \frac{1}{2}(\gamma -1)
\frac{c_V(T_\mathrm G - T_\mathrm W)P}{\sqrt{2\uppi\Rs T}}.
\label{sm5}
\end{equation}

Weil es keine netto Strömung im Gas gibt, muss der gesamte Energietransport durch eine gedachte Fläche
parallel
zur Wand ausschließlich über Diffusion stattfinden. Das heißt, es muss Gleichung \ref{sm6} gelten. $E_\mathrm
i- E_\mathrm r$ ist ja gerade die Differenz aus der gesamten Energie, die sich auf die Wand zu bewegt und der
Energie, die von der Wand zurückkommt ohne mit ihr kinetisch Energie ausgetauscht zu haben. Also muss
$E_\mathrm i - E_\mathrm r$ die Diffusion von Wärme senkrecht nur Wand sein:
\begin{equation}
E_\mathrm i - E_\mathrm r = \kappa \frac{\partial T}{\partial n}
\label{sm6}
\end{equation}

Setzt man Gleichung \ref{sm5} und Gleichung \ref{sm6} in Gleichung \ref{sm0} ein ergibt sich
\begin{equation}
\kappa \frac{\partial T}{\partial n} = \sigma_T\left( \frac{1}{2} \kappa \frac{\partial T}{\partial n} +
\frac{1}{2}(\gamma -1)
\frac{c_V(T_\mathrm G - T_\mathrm W)P}{\sqrt{2\uppi\Rs T}}\right).
\label{sm7}
\end{equation}
Löst man Gleichung \ref{sm7} nach $T_\mathrm G - T_\mathrm W$ auf und benutzt die Definitionen von
$\lambda$ (Gleichung \ref{mfwDef2}), $\mu$ (Gleichung \ref{viskT}) und der Prandtl"=Zahl $\mathit{P\!r} = \mu
c_P /
\kappa$, ergibt sich
\begin{equation}
T_\mathrm G - T_\mathrm W = \frac{2-\sigma_T}{\sigma_T} \frac{2\gamma}{\gamma+1}\frac{\lambda}{\mathit{P\!r}}
\frac{\partial T}{\partial n}.
\label{sm}
\end{equation}

Gleichung \ref{sm} ist Smoluchowskis Temperatursprungrandbedingung. In dieser Form wurde sie in dieser
Diplomarbeit implementiert.



\subsection{Randbedingungen höherer Ordnung}

Auch für Randbedingungen kann man einen allgemeinen Reihenansatz -- ähnlich wie bei der
Chapman-Enskok"=Methode
 -- machen und iterativ immer genauere Randbedingungen bestimmen. Im Bild einer solchen Reihenentwicklung sind
alle in diesem Kapitel vorgestellten Randbedingungen eine Näherung erster Ordnung. Maxwells
Randbedingung in zweiter Ordnung enthält zweite Ableitungen von $T$, $\rho$, $u_n$ und $u_t$, wie man zum
Beispiel \cite{PhysRevE.70.017303} entnehmen kann. Wegen ihrer Herleitung passen diese Randbedingungen
zweiter Ordnung zu den Burnett"=Gleichungen. Weil Randeffekte eine so große Rolle in der Simulation
verdünnter Gase spielen, gibt es auch Vorschläge, Randbedingungen zweiter Ordnung mit
Navier"=Stokes"=Gleichungen zu benutzen, um so numerischen Aufwand zu sparen. Die Vor- und Nachteile dieser
Herangehensweise sind genauso wie andere numerische Modelle für verdünnte Gase noch nicht genau bekannt,
vor allem weil es zu wenig experimentelle Daten gibt.


\chapter{Numerische Methoden und Implementierung}

\section{Überblick}

Für die Implementierung wurde die multi"=physikalische Finite"=Elemente Simulationssoftware Elmer \cite{elmer}
genutzt. Das Programm Gmsh \cite{gmsh} diente zur Gittergenerierung. Mit in C++
selbstgeschriebenen
Nachbearbeitungsprogrammen wurden aus den Simulationsergebnissen weitere physikalische Größen berechnet, die
einen Vergleich mit experimentellen Daten ermöglichen. Um Simulationen über einen großen Parameterbereich zu
iterieren, wurden selbstgeschriebene Python"=Skripte genutzt.

\section{Finite-Elemente-Methode}

Die Finite"=Elemente"=Methode (FEM) ist eine Methode zur Lösung von partiellen Differenzialgleichungen. Sie
geht von der variationellen Formulierung eines Randwertproblems aus. Ein großer Vorteil der FEM ist, dass die
Gitterelemente beliebige Polygone mit unterschiedlicher Orientierung und Größe sein können. Somit lassen
sich auch sehr unregelmäßige Geometrien leicht beschreiben und die Erhöhung der Elementdichte in Bereichen von
großem Interesse ist immer möglich. Im Folgenden wird die FEM anhand der von 
$n$ Variablen abhängigen allgemeinen eliptischen Differenzialgleichung 2. Ordnung mit
Dirichlet"=Randbedingungen erklärt. Dieses Kapitel
basiert auf dem Vorlesungsskript "`Numerisches Lösen von Differenzialgleichungen"' von Prof. Lubich. Es
werden folgende Konventionen für Indizes verwendet:

\begin{itemize}
\item
$N$ ist die Anzahl der gesamten Knotenpunkte; $i$ und $j$ sind Indizes, die von 1 bis $N$ laufen
\item
$R$ ist die Anzahl der Knotenpunkte pro Element; $r$ und $s$ sind Indizes, die von 1 bis $R$ laufen
\item
$E$ ist die Anzahl der Elemente; $e$ ist ein Index, der von 1 bis $E$ läuft.
\end{itemize}


Gesucht ist die Lösung $u$ der Differenzialgleichung \ref{poisson}. Gleichung \ref{varPoisson} ist die
variationelle Formulierung von Gleichung \ref{poisson}. Wenn ein $u$ gefunden wird, so dass Gleichung
\ref{varPoisson} für alle $v$ erfüllt ist, ist $u$ auch die Lösung von Gleichung \ref{poisson}. $\varOmega$
ist das Definitionsgebiet der Differenzialgleichung, $\varGamma$ der Rand des Gebietes und $\bar
\varOmega$ die Vereinigung aus $\varGamma$ und $\varOmega$.
\begin{equation}
  \left\{\begin{alignedat}{2}
  \sum_{i,j}^n a_{i,j} \frac{\partial}{\partial x_i} \frac{\partial}{\partial y_i} u &=f& \quad \text{in}\
&\varOmega\ ,\\
  u &= 0& \quad \text{auf}\ &\varGamma,
  \end{alignedat}\right.
  \quad\text{$\varOmega$ stückweise $C^1$, $f:  \bar\varOmega \to \mathbb R$ stetig.}
\label{poisson}
\end{equation}

\begin{align}
a(u, v) &= l(v)\quad\forall v \in V 
\label{varPoisson} 
\intertext{wobei gilt:}
a(u, v) &= \int_\varOmega \left\{\nabla u^\T \mathcal{A} \nabla v + a_0 uv\right\} \, \d \bm x \\
l(v) &= \int_\varOmega fv\, \d\bm x
\intertext{mit}
\mathcal{A} &= a_{i,j}\big|^n_{i,j=1}
\end{align}
Bei der variationellen Formulierung werden Dirichlet"=Randbedinungen in den Raum $V$ miteingebaut. Deshalb
muss $V$ der Sobolev"=Raum 1. Ordnung sein für dessen Elemente $v$ zusätzlich gilt $v(\varGamma)=0$. 


Die variationelle Formulierung der Differenzialgleichung wird in Gleichung \ref{approxPoisson} mit dem
Galerkin"=Verfahren angenähert.
\begin{equation}
a(u_N,v_N)=l(v_N)\quad\forall v_N\in V_N
\label{approxPoisson}
\end{equation}
Der unendlichdimensionale Funktionenraum $V$ wurde dabei durch einen endlichdimensionalen Raum $V_N$ ersetzt.
Die Dimension $N$ entspricht dabei der Knotenzahl der erforderlichen Diskretisierung von $\varOmega$. Die
Finite"=Elemente"=Methode ist die Anwendung des Galerkin"=Verfahrens mit dem Raum der stückweise stetigen
Polynome als Wahl für $V_N$. Also wird auch die gesuchte Funktion $u$ durch stückweise stetige Polynome
genähert.

Ein finites Element besteht aus Knotenpunkten $q_r$ und auf dem Element definierte Basisfunktionen $\phi_r$,
mit deren Hilfe die gesuchte Funktion $u$ angenähert wird. Ein Element benötigt so viele Knotenpunkte $q_r$
und Basisfunktionen $\phi_r$, wie der Funktionenraum $P$ der Polynome auf dem Element gewünschte Dimensionen
hat. Wählt man Dreiecke als Elemente und lineare Funktionen als Basisfunktionen, so sind drei Knotenpunkte
nötig; bei quadratischen Basisfunktionen auf Dreiecken braucht man sechs Knoten. Für die Basisfunktionen
$\phi_r$ muss gelten $\phi_r (q_s) = \updelta_{r,s}$. Das heißt, sie sind auf allen Knotenpunkten bis auf
einen 0. Damit läßt sich Polynom $p$ eines Elementes darstellen als $p = \sum_{r=1}^R p(q_r)\phi_r$.

Als erster Schritt zur Bestimmung von $u$ muss das Gebiet $\varOmega$ mit den finiten Elementen ausgefüllt
werden. Dabei dürfen sich die Elemente nie schneiden, jede Seite eines Elements muss mit der Seite eines
weiteren Elements oder einem Teil des Randes von $\varOmega$ übereinstimmen. Außerdem müssen Knoten auf
gemeinsamen Seiten übereinstimmen.

Als nächstes muss ein Gleichungssystem aufgestellt werden, das die Galerkin"=Approximation repräsentiert.
Die Galerkin"=Approximation \ref{approxPoisson} ist äquivalent zum Gleichungssystem $\mathbf A \bm\mu = \bm
b$,
wenn dieses wie in Gleichung \ref{steif1} bis \ref{defMu} definiert ist.

\begin{align}
  \mathbf A &=  a(\phi_j, \phi_i) \Big|^N_{i, j = 1}\\
  \label{steif1}
  \bm b &= l(\phi_i)\Big|_{i=1}^N \\
  \mu_i &= u_N(q_i) \quad u_N = \sum_{i=1}^N \mu_i \phi_i
  \label{defMu}
\end{align}

Aus historischen Gründen nennt man $\mathbf A$ die Steifigkeitsmatrix und $\bm b$ den Lastenvektor. Die
Einträge $a_{i,j}$ von $\mathbf A$ in Gleichung \ref{steif2} werden über die Elementmatrizen \ref{steifEle}
berechnet.

\begin{equation}
  a(\phi_i, \phi_j)
  = \int_\varOmega \left\{ \nabla \phi_i^\T \mathcal{A} \nabla \phi_j
  + a_0 \phi_i \phi_j \right\}\, \d \bm x
  \label{steif2}
\end{equation}

\begin{equation}
  a_{rs}^e = \int_{K^e} \left\{ (\nabla \phi_r^e)^\T \mathcal{A} \nabla \phi_s^e
  + a_0 \phi_r^e \phi_s^e\right\} \, \d \bm x
  \label{steifEle}
\end{equation}

Im Allgemeinen sind die Basisfunktionen der Elemente in einem lokalen Elementkoordinatensystem vorhanden und
müssen erst mit $\mathbf F^e$ auf das globale Koordinatensystem transformiert werden. Die Elementmatrix in
globalen Koordinaten aber in Abhängigkeit von lokalen Koordinaten ist in Gleichung \ref{steifEleLok} zu
sehen. Die lokalen Größen sind dabei mit einem \textasciicircum{} markiert. Der Operator $\mathrm D$ in
Gleichung \ref{steifEleLok} leitet jeden Eintrag in jeder Zeile von $\mathbf F$ partiell nach der Position
in der Zeile ab. In Gleichung \ref{steifEleLok} wurde der Index $e$ zwecks einfacherer Lesbarkeit weggelassen.

\begin{equation}
  a_{r,s} = \int_{\hat{K}} \left\{ \nabla \hat{\phi}_r^\T \hat{\mathcal{A}} \nabla \hat{\phi}_s
  + \hat{a}_0 \hat{\phi}_r \hat{\phi}_s\right\}
  \lvert \det \mathrm D \mathbf F \rvert \, \d\hat{\bm x}
\label{steifEleLok}
\end{equation}

Der globale Lastenvektor $\bm b$ \ref{lastEle} wird ebenfalls aus den lokalen Lastenvektoren \ref{lastEleLok}
berechnet. 

\begin{equation}
b_i = \int_\varOmega f \phi_i \, \d \bm x
\label{lastEle}
\end{equation}

\begin{equation}
b_i^e = \int_{\hat{K}} \hat{f} \hat{\phi}_r^e \lvert\det \mathrm D\mathbf F^e \rvert \, \d\hat{\bm x}
\label{lastEleLok}
\end{equation}

Um das LGS zu vervollständigen müssen noch die Randbedingungen $g$ in den Lastenvektor gesetzt werden.
Dazu werden die genäherten Randbedingungen $\tilde u_0$ wie in Gleichung \ref{randApprox} bestimmt und in das
LGS, wie Gleichung \ref{LGSmR} zeigt, eingebaut.

\begin{equation}
\tilde{u}_0 = \sum_{j\in J_0} g(q_j) \phi_j
\label{randApprox}
\end{equation}
\begin{equation}
\mathbf A\bm\mu = \bm b -  \sum_{j\in J_0} g(q_j) a(\phi_j, \phi_i)\Big|_{i=1}^N
\label{LGSmR}
\end{equation}

Als letzter Schritt muss noch das LGS gelöst werden. Dadurch erhält man die gesuchte näherungsweise Lösung
der Differenzialgleichung $u_N$. Da die Zahl der zu bestimmenden Variablen im LGS gleich der Zahl der
Knotenpunkte ist, wird das LGS in den meisten realen Anwendungen zu groß sein um es direkt lösen zu können.
Dass die Steifigkeitsmatrix $\mathbf A$, wegen der Eigenschaft der $\phi_i$ fast überall 0 zu sein, schwach
besetzt
ist, kann von den Lösungsalgorithmen ausgenutzt werden. Übliche iterative Verfahren zur Lösung von LGS, wie
sie bei der Lösung von Differenzialgleichungen auftreten, sind Krylow"=Unterraum"=Verfahren wie konjugierte
Gradienten oder Multigrid"=Verfahren. Falls als finite Elemente gerade Quadrate mit linearen Basisfunktionen
$\phi_r$ gewählt werden, ist die FEM
identisch zum Finite"=Differenzen"=Verfahren.

\subsection{Zusammenfassung der FEM}

\begin{itemize}
\item
Ausfüllen des Gebietes mit finiten Elementen
\item
Bestimmung der Basisfunktionen $\phi$ mit Hilfe der Knoten $q$
\item
Errechnung der nicht 0 Elemente der Steifheitsmatrix $\mathbf A$ mit numerischer Integration
\item
Errechnung der nicht 0 Elemente des Lastenvektors $\bm b$ mit numerischer Integration
\item
setze Randbedingungen in $\bm b$
\item
Lösen von $\mathbf A \bm\mu = \bm b$
\end{itemize}

\subsection{Fehlerquellen}

Bei der Finite-Elemente-Methode gibt es folgende Fehlerquellen:

\begin{itemize}
\item
Galerkin"=Approximation bei der der unendlichdimensionale Funktionenraum $V$ durch den endlichdimensionalen
Raum $V_n$ ersetzt wird
\item
numerische Integration bei Aufstellung der Steifigkeitsmatrix $\mathbf A$ und des Lastenvektors $\bm b$
\item
Approximation des Gebietrandes auf der die Differenzialgleichung definiert ist, falls er kein Polygon
beziehungsweise
Polyeder ist
\item
Lösen des linearen Gleichungssystems, falls kein direktes Verfahren verwendet wird
\item
Rundungsfehler
\end{itemize}


\section{Verwendete Modelle, Implementierung und Parameter}

\subsection{Implementierte Gleichungen}
\label{modelle}
Elmer löst die Navier"=Stokes \ref{NS} zusammen mit der Kontinuitätsgleichung \ref{konti} und die
Wärmeleitgleichung für kompressible Gase \ref{WL}.
\begin{equation}
\rho \left( \frac{ \partial \bm u}{\partial t} + ( \bm u \nabla ) \bm u \right) - \nabla \mathbf S = \rho
\bm a_\mathrm{ext}
\label{NS}
\end{equation}

\begin{equation}
\rho \left( \frac{ \partial \rho}{\partial t} + ( \bm u \nabla ) \rho \right) + \rho (\nabla \bm u) = 0
\label{konti}
\end{equation}

\begin{equation}
\rho c_V \left( \frac{ \partial T}{\partial t} + \bm u \nabla T \right) - \nabla ( \kappa \nabla T )
= -P \nabla \bm u + \phi( \mu \mathbf E )
\label{WL}
\end{equation}
Der Spannungstensor $\mathbf S$ in Gleichung \ref{ST} ist so definiert, dass er ein Newtonsches Fluid
voraussetzt. Der Term $\phi( \mu \mathbf E )$ in Gleichung \ref{WL} beschreibt die Erwärmung durch Reibung.
Äußere Beschleunigungen $\bm a_\mathrm{ext}$ wurden in dieser Diplomarbeit nicht genutzt.

\begin{equation}
\mathbf S = 2\mu \mathbf E - \frac{2}{3} \mu ( \nabla \bm u ) \mathbf I - P \mathbf I \quad
\text{mit }  E_{i,j} =  \frac{1}{2}  \left( \frac{\partial u_i}{\partial x_j} + \frac{\partial
u_j}{\partial x_i} \right)
\label{ST}
\end{equation}

Als Modell für die Kompressibilität dient die ideale Gas"=Gleichung \ref{idGas}.
\begin{equation}
PV = NkT \Longleftrightarrow P = \rho \Rs T
\label{idGas}
\end{equation}
$\Rs$ ist die spezifische Gaskonstante. Sie ergibt sich aus universeller Gaskonstante $R$ und molarer Masse
$M$, wie in Gleichung \ref{Rs} gezeigt wird.
\begin{equation}
\Rs = \frac{R}{M} = \frac{k}{\mt} = c_P - c_V = c_P (1 - \frac{1}{\gamma})
\label{Rs}
\end{equation}

Elmer bietet die Möglichkeit, selbstgeschriebene Fortran-90"=Routinen von der zentralen Konfiguratonsdatei für
die Simulation aufzurufen. Damit wurden zunächst temperaturabhängige Modelle für Viskosität $\mu$ und
Wärmeleitfähigkeit $\kappa$ implementiert. Die benutzten Gleichungen für Viskosität \ref{viskT} und
Wärmeleitfähigkeit \ref{wlT} sind mit Hilfe der statistischen Mechanik hergeleitet und wurden
\cite{hirschfelder} entnommen. Auf Grund ihrer Herleitung gelten sie eigentlich nur für einatomige ideale
Gase.

\begin{equation}
\mu = \frac{5}{16} \frac{k}{\sigma^2} \sqrt{\frac{T}{\Rs \uppi}}
\label{viskT}
\end{equation}


\begin{equation}
\kappa = \frac{15}{4} \Rs \mu \left( \frac{4}{15} \frac{c_V}{\Rs} + \frac{3}{5}\right) 
= \frac{5}{16} \frac{k}{\sigma^2} \sqrt{\frac{T \Rs}{\uppi}} \left(  \frac{1}{\gamma-1} + \frac{9}{4} \right) 
\label{wlT}
\end{equation}

Am wichtigsten ist natürlich die Implementierung von Maxwells thermische Transpirationsrandbedingung,
da sie erst für die Existenz des Effekts in der Simulation sorgt. Gleichung \ref{maxwell} ist Maxwells
Randbedingung für beliebige zweidimensionale Geometrien. In dieser Form ist sie \cite{PhysRevE.70.017303}
entnommen.
\begin{equation}
u_t(\varGamma) = \frac{2-\sigma_p}{\sigma_p} \lambda 
\left( \frac{\partial u_t}{\partial n} + \frac{\partial u_n}{\partial t} \right) 
+ \frac{3\mu}{4 \rho T} \frac{\partial T}{\partial t}
\label{maxwell}
\end{equation}

\begin{equation}
\lambda = \frac{kT}{\sqrt{2} \uppi \sigma^2 P}
\label{mfw}
\end{equation}

In Gleichung \ref{maxwell} ist $\varGamma$ die Grenzlinie zwischen Gas und Wand und $n$ und $t$ sind die
Koordinaten, die normal und tangential auf $\varGamma$ stehen. $\lambda$ ist in Gleichung \ref{mfw} für
ein"=atomige ideale Gase definiert. Maxwells Randbedingung ist eine Neumann"=Randbedingung, da sie Aussagen
über
die Ableitung der Geschwindigkeit auf dem Rand macht. Neumann"=Randbedingung für den Navier"=Stokes"=Löser von
Elmer müssen in Form von Spannungen gegeben werden. Werden in der Konfigurationsdatei Neumann"=Randbedingungen
aktiviert, gilt Gleichung \ref{elmerRand}, wobei der Slip"=Koeffizient $a$ und die externe Spannung $f$ vom
Benutzer in der Konfigurationsdatei oder in Form von selbstgeschriebenen Fortran 90 Routinen gesetzt werden
können.

\begin{equation}
a_i u_i - (\mathbf S \bm{\hat n})_i = f_i \quad i \in \{ x,y,z \} \text{ oder } i \in  \{n,t_1,t_2\}
\label{elmerRand}
\end{equation}

Da $(\mathbf S \bm{\hat n})_i$ die $i$-te Komponente des Produktes aus Spannungstensor \ref{ST} und dem
normierten Normalenvektor ist, gilt für zwei Dimensionen
\begin{equation}
(\mathbf S \bm{\hat n})_i = \mu \left(\frac{\partial u_x}{\partial y}+\frac{\partial u_y}{\partial x} \right).
\label{STi}
\end{equation}

Wenn man Gleichung \ref{maxwell} mit $a_i$ multipliziert, erhält man
\begin{equation}
a_i u_i = a_i \frac{2-\sigma_p}{\sigma_p} \lambda \left(\frac{\partial u_x}{\partial y}+\frac{\partial
u_y}{\partial x} \right)
+ a_i \frac{3\mu}{4 \rho T} \frac{\partial T}{\partial t}.
\label{maxMod}
\end{equation}
Gleichung \ref{maxMod} hat nun die selbe Struktur wie Gleichung \ref{elmerRand} und 
die Summanden können miteinander identifiziert werden. Setzt man die rechte Seite 
von Gleichung \ref{maxMod} und $(\mathbf S \bm{\hat n})_i = \mu\left(\frac{\partial u_x}{\partial y}
+ \frac{\partial u_y}{\partial x} \right)$ 
in Gleichung \ref{elmerRand} ein, erhält man:
\begin{align}
\label{aHerleitung1}
& a_i \frac{2-\sigma_p}{\sigma_p} \lambda \left(\frac{\partial u_x}{\partial y} 
+ \frac{\partial u_y}{\partial x} \right)
+ a_i \frac{3\mu}{4 \rho T} \frac{\partial T}{\partial t}
- \mu\left(\frac{\partial u_x}{\partial y}+\frac{\partial u_y}{\partial x} \right) = f_i \\
\Longleftrightarrow &a_i = \frac{\mu}{\beta \lambda} \quad \text{mit} \quad \beta =
\frac{2-\sigma_p}{\sigma_p}
\label{aHerleitung2}
\end{align}

Setzt man nun in Gleichung \ref{aHerleitung2} $\mu$ wie in Gleichung \ref{viskT} und $\lambda$ wie in
Gleichung \ref{mfw} definiert ein, ergibt sich der Slip"=Koeffizient $a_i$ zu:

\begin{equation}
a_i = \frac{5 \sqrt{2}}{16} \frac{P}{\beta} \sqrt{\frac{\uppi}{\Rs T}}
\label{ai}
\end{equation}

In Gleichung \ref{aHerleitung2} erkennt man ebenfalls, dass die externe Spannung $f_i$, die man in Elmer
setzten kann, gerade $a_i \frac{3\mu}{4 \rho T} \frac{\partial T}{\partial x}$ sein muss. Ersetzt man auch
hier $a_i$, $\mu$ und $\rho$ durch ihre jeweiligen Definitionen ergibt sich für $f_i$ schließlich 


\begin{equation}\label{fi}
\begin{split}
f_i &= a_x \frac{3\mu}{4 \rho T} \frac{\partial T}{\partial t} \\
    &= \frac{\mu}{\beta \lambda} \frac{3\mu}{4 \rho T} \frac{\partial T}{\partial t} \\
    &= \frac{3}{4} \beta^{-1} \frac{\sqrt{2} \uppi \sigma^2 P}{kT} \frac{5^2}{16^2} \frac{\uppi}{\Rs T}
\frac{k^2 T^2}{\uppi^2 \sigma^4} \frac{\Rs T}{P} \frac{1}{T} \frac{\partial T}{\partial t} \\
    &= \frac{3}{4} \frac{5^2}{16^2} \beta^{-1} \frac{\sqrt{2} k}{\sigma^2} \frac{\partial T}{\partial t} \\
    &= \frac{75 \sqrt{2}}{1024} \beta^{-1} \frac{k}{\sigma^2} \frac{\partial T}{\partial t}.
\end{split}
\end{equation}


Des Weiteren wurden Smoluchowskis Temperaturrandbedingungen \ref{TJ} implementiert. Das Vorzeichen auf der
rechten Seite von Gleichung \ref{TJ} ist so gewählt, dass die Randbedingung die von Elmer verwendete
Konvention der Richtung des Normalenvektors $\bm n$ übereinstimmt.
\begin{equation}
T_\mathrm{Gas} - T_\mathrm{Wand} = -b(T,P) \frac{\partial T}{\partial n} \quad 
\text{mit } b(T,P) = \frac{2-\sigma_T}{\sigma_T} \frac{2\gamma}{\gamma+1} \frac{\lambda}{\mathit{P\!r}}
\label{TJ}
\end{equation}
Gleichung \ref{elmerTrand} zeigt die Wärmeflussrandbedingung, die Elmer verwendet. Es herrscht formale
Gleichheit mit Gleichung \ref{TJ}, wenn man $b$ mit $\frac{\kappa}{\alpha}$ und $T_\mathrm{Wand}$ mit 
$T_\mathrm{ext}$ identifiziert. 
\begin{equation}
-\kappa \frac{\partial T}{\partial n} = \alpha ( T - T_\mathrm{ext} )
\label{elmerTrand}
\end{equation}
In Rechnung \ref{alpha} sieht man, wie der Koeffizient $\alpha$ durch eine Fortran 90 Routine gesetzt werden
muss, wenn man Smoluchowskis Randbedingung mit Hilfe von Elmers Wärmeflussrandbedingung implementieren will.
Es wurden dabei das temperaturabhängige Modell für die Wärmeleitfähigkeit $\kappa$ aus Gleichung \ref{wlT} und
die Definition der Prandtlzahl für ideale Gase $\mathit{P\!r} = \frac{4 \gamma}{9 \gamma -5}$ benutzt.
\begin{equation}\label{alpha}
\begin{split}
\alpha = \frac{\kappa}{b} &= \frac{5}{16} \frac{k}{\sigma^2} \sqrt{\frac{T \Rs}{\uppi}} \left( 
\frac{1}{\gamma -1} + \frac{9}{4} \right) \frac{\sigma_T}{2-\sigma_T} \frac{\gamma+1}{2\gamma} \frac{\sqrt{2}
\uppi \sigma^2 P}{kT} \frac{4\gamma}{9\gamma -5} \\
&= \frac{5}{8} \sqrt \frac{2\uppi \Rs}{T} \frac{\sigma_T}{2-\sigma_T} \left(  \frac{1}{\gamma -1} +
\frac{9}{4} \right)
\frac{\gamma +1}{9\gamma -5} P \\
&= \frac{5}{8} \sqrt \frac{2\uppi \Rs}{T} \frac{\sigma_T}{2-\sigma_T} \frac{\gamma+1}{4(\gamma-1)} P \\
&= \frac{5 \sqrt{2\uppi}}{32 \left(  \frac{2}{\sigma_T}-1\right) } \frac{\gamma+1}{\gamma-1}
\sqrt{\frac{\Rs}{T}} P
\end{split}
\end{equation}


\subsection{Elmer}

Elmer ist eine Finite"=Elemente"=Software zur Lösung von partiellen Differenzialgleichungen und kann viele
verschiedene physikalische Gleichungen lösen. Diese können auch miteinander kombiniert werden. Für diese
Diplomarbeit wurden die Löser für die Navier"=Stokes"=Gleichung \ref{NS} mit Kontinuitätsgleichung \ref{konti}
und für die Wärmeleitgleichung \ref{WL} genutzt. Elmer geht beim Numerischen Lösen von Differenzialgleichungen
nach folgendem Schema vor:

\begin{enumerate}
\item
Zeitintegration: Die äuserste Iterationshierarchie ist die Zeitintegration, diese kann über die
Crank"=Nicolson"=Methode
oder eine Backward Differences Formulae (BDF) erfolgen.
\item
Gleichgewichtszustand für aktuellen Zeitschritt erreichen: In dieser Iterationsebene findet das Linearisieren
der wahrscheinlich nicht linearen Zusammenhänge zwischen den gesuchten Variablen, die durch die Kopplung von
mehreren physikalischen Gleichungen entstehen, statt. Des Weiteren wird die Synchronisation zwischen den zu
lösenden gekoppelten Gleichungen durchgeführt. Die Linearisierung wird ähnlich durchgeführt wie im nächsten
Schritt. Falls in dem aus den gekoppelten Gleichungen erzeugte Gleichungssystem $\mathbf A \bm x = \bm b$
 $\mathbf A$ und $\bm b$ von $\bm x$ abhängen, wird es in das Gleichungssystem $ \mathbf A(\bm x_{i-1}) \bm
x_i = \bm b(\bm x_{i-1})$ überführt.

\item
Linearisieren der einzelnen Gleichungen: Falls die mit jedem einzelnen Löser zugeschalteten Gleichung in sich
nicht"=lineare DGLs sein
sollten, werden diese linearisiert. Das heißt, falls im von der Finiten"=Elemente"=Methode erzeugte
Gleichungssystem $\mathbf A \bm x = \bm b$ $\mathbf A$ und $\bm b$ von $\bm x$ abhängen, wird es in das
Gleichungssystem $\mathbf A(\bm x) \bm x = \bm b(\bm x) \rightarrow \mathbf A(\bm x_{i-1}) \bm x_i = \bm b(\bm
x_{i-1})$ überführt, bei dem in jedem einzelnen Iterationsschritt ein Gleichungssystem zu lösen ist, das
linear in $\bm x_i$ ist.

\item
LGS lösen: Die innerste Iterationshierarchie ist das Lösen der linearen Gleichung $\mathbf A \bm x = \bm b$.
Dazu stehen
verschiedene numerische Lösungsmethoden aus 3 verschiedenen Methoden"=Familien zur Verfügung. Implementiert
sind
unter anderem die direkten LGS"=Löser LAPACK und UMFPACK \cite{992206}, die Krylov"=Unterraum"=Methoden
"`Konjungierte Gradienten"' (CG), BiCGStab und GMRES sowie die Multigitter"=Methoden GMG und AMG. Um die
Konvergenz der Krylov"=Unterraum"=Methoden zu verbessern, können Vorkonditionierer wie die unvollständige
LR"=Zerlegung benutzt werden.
\end{enumerate}


Die Zeitintegration kann auch weggelassen werden. Die Gleichgewichtsiteration ist dann die äußerste Hirarchie.
In diesem Modus versucht Elmer eine Lösung der DGL zu finden, die nicht mehr von der Zeit abhängt, also dem
Gleichgewichtszustand entspricht. Es ist dann aber nicht mehr möglich eine Aussage über die Zeit zu machen,
die
das System benötigt hat um in den Gleichgewichtszustand zu kommen.

Um die Chancen auf Konvergenz beim Linearisieren in Schritt 2 oder 3 zu erhöhen, kann ein "`Relaxation
Factor"' $\alpha$ gewählt
werden. Damit werden die $\bm x_i$ eines Iterationsschrittes durch $\bm x_i' =\alpha \bm x_i + (1-\alpha ) \bm
x_{i-1}$ ersetzt.


\subsection{Gewählte Parameter}

Die materialabhängigen Parameter wurden so gewählt, dass sie trockener Luft entsprechen. Die in Kapitel
\ref{modelle} verwendeten Gleichungen sind aber aus den Voraussetzungen für ideale einatomige Gase aus
"`harten Kugeln"' hergeleitet. Luft ist aber ein vor allem aus zwei zwei"=atomigen Molekülsorten ($\mathrm
N_2, \mathrm O_2$) bestehendes Gemisch. Da die beiden wichtigsten Parameter -- der Durchmesser eines Moleküls
im
Hartkugelmodell $\sigma$ und die Massen der Moleküle von Sauerstoff und Stickstoffen -- sehr änlich sind, wird
kein zu großer Fehler gemacht, wenn Durchschnittswerte dieser Parameter in Modelle für einatomige ideale Gase
eingesetzt werden. Das Massenverhältnis ist $m_\mathrm O / m_\mathrm N = 1.14$ und das Verhältnis der
Moleküldurchmesser ist $\sigma_{\mathrm N_2}/\sigma_{\mathrm O_2} = 1.039$. Die Zahlen für den
Moleküldurchmesser des Hartkugelmodells wurden \cite{hirschfelder} entnommen; sie stammen aus
Viskositätmessungen. Wo die größere Anzahl von Freiheitsgraden von Molekülen gegenüber Atomen eine Rolle
spielt,
wie zum Beispiel in der Formel für die Wärmeleitfähigkeit, wurde sie durch entsprechende Werte $\gamma$ und
$c_P$ berücksichtigt.

Für Werte für den tangential momentum accommodation coefficient (TMAC) $\sigma_p$ gibt es nur wenige Quellen
und es wurden keine gefunden, die ihn für über Wände aus Glimmer strömende Luft messen. In einer Doktorarbeit 
\cite{Arkilic97measurementof} wurden aber der TMAC von verschiedenen über Silizium strömende Gase
untersucht, wobei diese immer um 0.8 lag und in einem Artikel \cite{maali:027302} wurde $\sigma_p$ für über
Glas strömende Luft abgeschätzt; das Ergebnis ist ebenfalls 0.8. Es ist zu vermuten, dass der TMAC bei
Lichtmühlen sehr hoch liegt, sonst würden sie ja nicht gut funktionieren, deshalb wurde in dieser
Diplomarbeit $\mathrm{TMAC} = 1$ gewählt. Wie im Theorieteil dargelegt, beschreiben TMAC und TAC das
Verhältnis von diffuser und reflektierender Wechselwirkung zwischen Teilchen und Wand. Daher müssen TMAC und
TAC für die selbe Wand"=Gas"=Kombination sehr ähnlich sein. Aus diesen Gründen wurde der TAC auch auf 1
gesetzt.

In Tabelle \ref{tab:parameter} sind alle benutzten Materialgrößen zusammengefasst. Falls
nicht ausdrücklich bei der Präsentation von Ergebnissen andere Werte genannt werden, wurden diese Werte
verwendet.

\begin{table}
\capstart\centering
\begin{tabular}{p{0.5\linewidth}|l}
Name						& Wert \\
\hline
Adiabatenexponent $\gamma$	& 1.4 	 \\
\hline
Isobare spezifische Wärmekapazität $c_P$	& $1.01\,\mathrm{\frac{kJ}{kg\,K}}$  \\
\hline
Moleküldurchmesser (Hartkugelmodell) $\sigma$	& 0.37\,nm \\
\hline
TMAC $\sigma_p$		& 1.0   \\
\hline
TAC $\sigma_T$		& 1.0 
\end{tabular}
\caption{Werte verwendeter Materialparameter}
\label{tab:parameter}
\end{table}



\section{Verifikation der Implementierung}

\subsection{Test von Maxwells Geschwindigkeitsrandbedingungen}

Um zu prüfen, ob die Implementierung von Maxwells Randbedingung \ref{maxwell} richtige Ergebnisse liefert,
wurde ein von \cite{karniadakis2005man} und \cite{JeCz07} benutzter Testfall gerechnet und die Resultate mit
diesen Quellen verglichen. Der Testfall besteht aus aus zwei Containern, die mit einem Kanal verbunden sind.
Er ist in Abbildung \ref{fig:weightGrid} dargestellt. In Tabelle \ref{tab:Testcase} sind die Werte der
physikalischen Parameter dieses Testfalls aufgelistet. Auf den Wänden des Kanals wurden ein Temperaturwerte
gesetzt, die linear mit $x$ von $T_\mathrm L$ nach $T_\mathrm R$ verlaufen; auf den Kanalwänden herrscht somit
ein linearer Temperaturgradient.

\begin{figure}
\capstart\centering
\includegraphics[width=0.4\linewidth,angle=-90]{fig/weightGrid3}
\caption{Maßstabgetreue Abbildung des Gitters des Verifikationstestfalls. Der Kanal hat ein
Höhen"=Längen"=Verhältnis von 1:10, relativ dazu haben die Container je eine Höhe von 11 Kanalhöhen und eine
Länge von 5 Kanalhöhen}
\label{fig:weightGrid}
\end{figure}

\begin{table}
\capstart\centering
\begin{tabular}{l|l}
Name						& Wert \\
\hline
Referenzdruck $P_0$	& 101325\,Pa	 \\
\hline
Temperatur links $T_\mathrm L$	& 300\,K \\
\hline
Temperatur rechts $T_\mathrm R$	& 400\,K \\
\hline
Moleküldurchmesser (Hardkugelmodell) $\sigma$	& 0.3678\,nm  \\
\end{tabular}
\caption{Verwendete Parameter für den Verifikationstestfall. Weitere benötigte Parameter sind identisch mit
Tabelle \ref{tab:parameter}.}
\label{tab:Testcase}
\end{table}

In \cite{karniadakis2005man} wird eine "`Direct Simulated Monte Carlo"' (DSMC) Methode verwendet und in
\cite{JeCz07} der gleiche Ansatz wie in dieser Diplomarbeit: Verwendung von Maxwells Randbedingung mit
Navier"=Stokes"=Gleichung.

In der Simulation führt die thermische Transpiration zu einem Gasfluss an den Kanalwänden in Richtung des
wärmeren, rechten Containers. Daraus ergibt sich ein höherer Druck im rechten Container, der für einen
Gasfluss in der Mitte des Kanals in Richtung des linken Containers sorgt. Im Gleichgewichtszustand sind die
Massenflüsse in beide Richtungen gleich groß und die Druckdifferenz zwischen den beiden Containern bleibt auf
einem konstanten Wert.

Diese Simulation wurde für drei verschiedene Knudsenzahlen durchgeführt. Als systembeschreibende Knudsenzahl
wurde die Knudsenzahl in der Mitte des Kanals gewählt dort war der Druck gerade der Referezdruck $P_0$ und die
Temperatur $\frac{T_\mathrm R + T_\mathrm L}{2} = 350\,\mathrm K$. Als "`Abmessung des Systems"' $d$ diente
die Kanalhöhe $h$. $\kn$ wurde nur über die Kanalhöhe variiert, der Referenzdruck $P_0$ und die Temperaturen
blieben immer die gleichen. Da immer das gleiche Gitter verwendet wurde, bestimmte die Kanalhöhe $h$ auch die
physikalischen Abmessungen des Gitters.

In Abbildung \ref{fig:weightv} sieht man, wie die thermische Transpiration für einen Gasfluss entlang der
Kanalwände in Richtung des wärmeren Containers sorgt. Die Geschwindigkeiten sind um so höher, je größer $\kn$
ist. Im dort abgebildeten Gleichgewicht sind auch die Rückflüsse mit steigender $\kn$ immer größer, da sie
durch eine höhere Druckdifferenz in den Containern verursacht werden. Die verschiedenen Druckdifferenzen
$\Delta P$, die sich für bestimmte $\kn$ einstellen, sind in Tabelle \ref{tab:deltaP} aufgelistet, und mit
den
Werten einer DSMC Simulation von \cite{karniadakis2005man} verglichen. Während bei niedrigen $\kn$ sehr gute
Übereinstimmung herrscht, werden die Abweichungen für hohe $\kn$ immer größer. Das kann von der in Kapitel
\ref{theorie} besprochenen besseren Eignung von Teilchenmodellen gegenüber Kontinuumsmodellen für 
hohe $\kn$ herrühren.

\begin{table}
\capstart\centering
\begin{tabular}{l|l|l|l|l}
Kanalhöhe $h$ (\textmu m) & $\kn$ & $\Delta P_\mathrm{lit}$ (Pa) nach \cite{karniadakis2005man} & $\Delta P$
(Pa) & $\frac{\Delta P_\mathrm{lit} - \Delta P}{\Delta P_\mathrm{lit}}$ \\
\hline
1.5 & 0.053 & 336.0 & 335.9 & 0.0004\\
\hline
0.65 & 0.122 & 1409.0 & 1352.6 & 0.05\\
\hline
0.22 & 0.36 & 8832.0 & 6425 & 0.27\\
\end{tabular}
\caption{Vergleich der Druckdifferenzen der eigenen Simulation mit den numerischen Ergebnissen von
\cite{karniadakis2005man}. }
\label{tab:deltaP}
\end{table}

\begin{figure}
\capstart\centering
\include{w_y_vx}
\caption{Geschwindigkeitsprofile für die 3 verschiedenen Kanalhöhen $h$ beziehungsweise Knudsenzahlen $\kn$
des
Verifikationstestfall im Gleichgewichtszustand. Es wurde die Parallele zur y"=Achse ausgeschnitten, die in
der Mitte des jeweiligen Kanals liegt.}
\label{fig:weightv}
\end{figure}

Abbildung \ref{fig:weightp} zeigt den Verlauf des Drucks im Gleichgewichtszustand für die drei
unterschiedlichen
$\kn$"=Werte aus Tabelle \ref{tab:Testcase}. Man sieht deutlich die Zunahme der Druckdifferenzen zwischen den
beiden Containern für steigende Knudsenzahlen. Die Druckwerte sind auf deiner Parallelen zur x"=Achse durch
die Mitte des Kanals.

\begin{figure}
\capstart\centering
\include{w_x_p}
\caption{Durch $P_0$ normierter Druckverlauf durch das gesamte Gitter entlang einer durch die Mitte
des Kanals verlaufenden Parallelen zur x"=Achse. Die Simulation hat für alle drei $\kn$ den
Gleichgewichtszustand erreicht.}
\label{fig:weightp}
\end{figure}

\subsection{Test von Smoluchowskis Temperaturrandbedingungen}

Smoluchowskis Temperaturrandbedingungen wurden mit einem anderen Testfall erprobt, da in einem langen Kanal,
in dem die gegenüberliegenden Seiten die gleiche Temperatur haben, $\frac{\partial T}{\partial n}$ fast exakt
0
ist. Es wurde ein Testfall aus \cite{JeCz07} nachgerechnet, der aus einem Rechteck besteht, auf dessen obere
Wand
350\,K und auf allen anderen 3 Wänden 300\,K gesetzt sind. Wären auf den Rändern normale
Dirichlet"=Randbedingungen gesetzt, würde der Temperaturverlauf im Gas von unten nach oben bei genau 300\,K
beginnen und bei genau 350\,K enden.
Mit Smoluchowskis Temperaturrandbedingungen weicht die Temperatur des
Gases an der Wand von der Temperatur der Wand ab. Diese Abweichung wird mit steigender Knudsenzahl $\kn$
größer. Abbildung \ref{fig:testTJ} zeigt dieses Verhalten als Ergebnis des Testfall aus \cite{JeCz07}. Es
wurde der Temperaturverlauf für verschiedene $\kn$"=Werte auf einer Linie durch die Mitte des Rechtecks von
unten nach oben gezeichnet. Der Referenzdruck blieb immer konstant. Somit wurde die Knudeszahl nur über die
Rechteckshöhe variiert. In Tabelle \ref{tab:testTJ} stehen die für diesen Testfall verwendeten Parameter.
Einige Parameter des Temperatursprung Testfalls fehlten in \cite{JeCz07} beziehungsweise waren falsch
angegeben. Das
wurde in einer E"=Mail"=Korrespondenz mit einem der Autoren geklärt.

Das Rechteck hat ein Längen"=Höhenverhältnis von 30:1 für alle verwendeten Höhen $h$. $\kn$ wurde gleich wie
im Testfall für Maxwells
Randbedingung bestimmt: Also die mittlere freie Weglänge $\lambda$ in der Mitte des Rechtecks mit Gleichung
\ref{mfw} errechnet und durch die eingestellte Höhe $h$ des Rechtecks geteilt.

Der Verlauf der Kurve in \ref{fig:testTJ} stimmt für $\kn = 0.1$ sehr gut mit \cite{JeCz07} überein. Für
$\kn =
0.5$ und $\kn = 0.77$ weichen die Schnittpunkte der Kurven mit der $y = 0$ Achse leicht von denen in
\cite{JeCz07} ab. Die relative Abweichung der Schnittpunkte liegt unter 1\%. Ob diese kleinen Abweichungen
tatsächlich durch die unterschiedlichen Implementierungen
hervorgerufen wurden oder es trotz des E"=Mail"=Kontakts mit den Autoren kleine Unterschiede zwischen den
verwendeten Parametern gab, konnte nicht abschließend geklärt werden.


\begin{figure}
\capstart\centering
\include{testTJ}
\caption{Vergleich der durch Temperatursprungrandbedingung erzeugten Temperaturverläufe durch die Mitte des
 Rechtecks für verschiedene Höhen $h$ beziehungsweise $\kn$"=Werte und normalen Dirichlet"=Randbedingungen.}
\label{fig:testTJ}
\end{figure}

\begin{table}
\capstart\centering
\begin{tabular}{l|l}
Name						& Wert \\
\hline
Referenzdruck $P_0$	& 101325\,Pa	 \\
\hline
Temperatur oben	& 350\,K \\
\hline
Temperatur rechts links unten	& 300\,K \\
\hline
Moleküldurchmesser (Hartkugelmodell) $\sigma$	& 0.31623\,nm  \\
\hline
Isobare spezifische Wärmekapazität $c_P$	& $1006\,\mathrm{\frac{J}{kg\,K}}$ \\
\end{tabular}
\caption{Verwendete Parameter für den Temperatursprungtestfall. Weitere benötigte Parameter sind identisch
mit Tabelle \ref{tab:parameter}.}
\label{tab:testTJ}
\end{table}



\chapter{Ergebnisse}

\section{Simulation einer Lichtmühle}
\label{simlicht}
Mit der getesteten Implementierung der Navier"=Stokes"=Gleichungen mit Maxwells und Smoluchowskis
Randbedingungen
wurden die Drehmomente simuliert, die auf eine handelsübliche Lichtmühle wirken. Die Abmessungen
und Parameter der Simulation wurden an die Experimente von \cite{duffner} und \cite{geyer} angepasst, die
zum Teil die gleiche Lichtmühle vermessen haben. Die Form und Ausrichtung der Flügel in den Experimenten
ist gleich wie diejenige in Abbildung \ref{fig:lightmill}. Weil nur 2D"=Simulationen durchgeführt
wurden,
mussten die 3D Abmessungen der echten Lichtmühle übertragen werden. Der quaderförmigen Flügel
der experimentell untersuchten Lichtmühle misst $2.32\,\mathrm{cm} \times 2.32\,\mathrm{cm}\times
0.015\,\mathrm{cm}$. Deshalb hat der simulierte 2D Flügel die Maße $2.3\,\mathrm{cm}\times
0.015\,\mathrm{cm}$, wobei auf je zwei gültige Stellen gerundet wurde. Der simulierte Flügel
liegt in einem abgeschlossenen rechteckigen Bereich, der mit Luft gefüllt ist und dessen senkrechter Abstand
2\,cm zu den Flügelseiten ist. Die simulierte Konfiguration ist in Abbildung \ref{fig:paddle} dargestellt.

\begin{figure}
\capstart\centering
\includegraphics[width=0.8\linewidth]{fig/trptest.png}
\caption{Gitter mit 181466 Elementen des $2.3\,\mathrm{cm}\times 0.015\,\mathrm{cm}$ großen Flügels, der in
einem 2\,cm von jeder Flügelseite entfernten Kasten liegt. In der Nähe der Flügelkanten ist die Elementdichte
besonders hoch.}
\label{fig:paddle}
\end{figure}

Um unphysikalische Spitzen in den Kräften pro Koten $\bm F_i$ auf den Eckknoten zu vermeiden, wurden die vier
Ecken des
Flügels durch Viertelskreise mit einem Radius $r$ von einem sechstel der Flügeldicke $d$ ersetzt.
$d=0.015$\,cm führt deshalb auf Eckkreisradius $r=25\,\text{\textmu}\mathrm m$. Das Aussehen einer
abgerundeten
Flügelkante ist
in Abbildung \ref{fig:paddle2}, die einen stark vergrößerten Ausschnitt aus Abbildung \ref{fig:paddle} zeigt, 
zu sehen. Leider traten diese Spitzen in den $\bm F_i$ weiterhin auf, wenn auch in abgeschwächter Form. In
Abschnitt \ref{kraft} wird darauf genauer eingegangen.

\begin{figure}
\capstart\centering
\includegraphics[width=0.7\linewidth]{fig/trpCloseUp}
\caption{Stark vergrößerter Ausschnitt des Bereiches um die linke Flügelkante von Abbildung
\ref{fig:paddle}}
\label{fig:paddle2}
\end{figure}

Das Flügelinnere wurde nicht mitsimuliert, sondern nur die Strömung um den
Flügel mit allen Modellen aus Kapitel \ref{modelle}. Das heißt, der Flügel bestand in der Simulation nur aus
den auf seinem Rand gesetzten Randbedingungen. Auch die einfallende Lichtstrahlung wurde nicht simuliert,
sondern es wurde angenommen, dass das einstrahlende Licht zu einer höheren Temperatur auf der schwarzen Seiten
führt.
Also wurden auf die beiden langen Seiten des 2D"=Flügels zwei unterschiedliche Temperaturen gesetzt und auf
den beiden kurzen Seiten ein linearer Temperaturverlauf von der kalten zur warmen Seite.

Die Lichtmühle in der Simulation ist statisch, das heißt der Flügel bewegt sich nicht. Die durch Maxwells
Randbedingungen hervorgerufenen Strömungen um die Flügelkanten sorgen aber für eine viskose Kraft $\bm F$ auf
den Flügel, die in der Simulation errechnet wird und aus der sich ein auf den Flügel wirkendes Drehmoment
bestimmen lässt.

Die Kraft $F_x$, die durch eine viskose Strömung in x"=Richtung über eine weit ausgedehnte Fläche $A$, die in
der x"=z"=Ebene liegt, erzeugt wird, ist
\begin{equation}
F_x = \mu \frac{\partial u_x}{\partial y}A.
\label{visF}
\end{equation}
Verallgemeinert man die Vorrausetzungen von Gleichung \ref{visF} auf eine beliebige Strömung $\bm u$, die
über eine beliebige Oberfläche $A$ fließt, ergibt sich
\begin{equation}
\bm F = \int_A \mathbf S\hat{ \bm{n}} \,\d o.
\label{visF3D}
\end{equation}
In Gleichung \ref{visF3D} ist $\hat{ \bm{n}}$ der auf 1 normierte Normalenvektor und $\mathbf S$ der
Spannungstensor aus Gleichung \ref{ST}. Man sieht in Gleichung \ref{visF}, dass die Richtung der Kraft
nicht von der Richtung der Strömung sondern von dem Vorzeichen der Ableitung abhängt. Thermische
Transpiration ist eine Strömung, die auf der Wand am stärksten ist und abnimmt je weiter man sich von ihr
entfernt, deshalb erzeugt sie eine Kraft auf die Wand entgegen ihrer Fließrichtung.

Elmer bietet mehrere Methoden, die durch viskose Strömung auf Ränder wirkende Kräfte pro Knoten
$\bm F_i$, zu bestimmen; eine davon ist die numerische Ausführung des Integrals aus Gleichung \ref{visF3D}.
Mit den $\bm F_i$ ergibt sich
das Drehmoment $\bm D$ über das Kreuzprodukt $\bm D = \sum_i \bm r_i \times \bm F_i $, wobei $i$ über alle
Knoten des Flügelrandes läuft. $\bm r_i$ sind die Abstände von der Drehachse zu den Flügelknoten. Die
Drehachse wurde auch nicht mitsimuliert, sondern ist nur eine Koordinate bei der Berechnung des Drehmoments.
Diese Koordinate wurde so gewählt, dass der Abstand der Drehachse zum Flügelschwerpunkt in Experiment und
Simulation gleich ist. Dieser Abstand ist 2.02\,cm. Somit sollte sich die Drehachse am Ende einer geraden 
Linie mit 0.86\,cm Länge, die senkrecht auf der Mitte der rechten kurzen Seite steht, befinden. Die
0.86\,cm wurden auf 1\,cm gerundet.

Es spielt keine Rolle, dass die Lichtmühle im Experiment relativ zur Drehachse rautenförmig, aber in
der Simulation dagegen  mit einer der beiden langen Seiten parallel zur Drehachse steht, so lange der Abstand
von Schwerpunkt und Drehachse gleich groß ist. Das liegt an der Symetrieeigenschaft des Trägheitstensors eines
Quaders.

Die 2D"=Simulation produziert Strömungen um den Flügel, die einem in die dritte Dimension (z"=Richtung)
unendlich ausgedehnten Flügel entsprechen. Die von Elmer errechneten $\bm F_i$ sind so normiert, dass sie den
Kräften auf einem in z"=Richtung 1\,m großen Teilstück entsprechen. Um die Drehmomente der Simulation mit dem
Experiment vergleichen zu können, wurde das resultierende Drehmoment mit dem Faktor 0.023 multipliziert; also
entsprechend den 2.3\,cm, die der reale Flügel eigentlich hoch ist. Des Weiteren wird das Drehmoment noch
mit zwei multipliziert, da in es in der Simulation nur zwei Kanten mit Temperaturgradient gibt, bei einem
3D"=Flügel dagegen vier. Als letztes wird noch ein Faktor 4 hinzumultipliziert, weil die Lichtmühle im
Experiment vier Flügel hat. Das resultierende Drehmoment $D_\mathrm{res}$ ist also $D_\mathrm{res} = 2\cdot4
l_\mathrm{rel} D_\mathrm{sim}$ mit der relativen Länge des Flügels $l_\mathrm{rel} = 0.023$.




\section{Experimentelle Daten}
\label{exp}
Die Parameter der Lichtmühlensimulation sind so gewählt, dass die Simulation mit Messungen einer
Zulassungsarbeit und einer Diplomarbeit vergleichbar ist. In Herr Geyers Zulassungsarbeit
\cite{geyer} wurden Druck"=Drehfrequenz- sowie Druck"=Drehmoment"=Messungen durchgeführt. In Herr Duffners
Diplomarbeit \cite{duffner} wurde zusätzlich noch die Lichtleistung variiert. Der Artikel \cite{arenas}
lieferte ebenfalls Lichtleistung"=Drehfrequenz"=Daten. Die Drehmomente wurden aus Messungen der Abklingzeiten
der Bewegung nach Ausschalten der Beleuchtung unter Annahme von Stokesscher Reibung ($F_\mathrm{reib} \sim v$)
bestimmt. Aus diesen Messungen lassen sich folgende Zusammenhänge ableiten:

\begin{itemize}
\item Sowohl Drehgeschwindigkeit als auch Drehmoment $D$ wachsen linear mit der von der Beleuchtungsanlage
aufgenommenen Leistung.
\item
Wie in Abbildung \ref{fig:pfdata} zu sehen ist, hat die Drehfrequenz $f$ ein klares Maximum bei einem
bestimmten Druck $P_\mathrm{opt}$. Die Frequenz $f$ fällt sowohl bei gegen 0 strebende $P$ als auch bei
großen $P$ auf 0 ab, wobei der Abfall von $f$ für große $P$ exponentiell zu sein scheint.
\item
$P_\mathrm{opt}$ ist nicht von der Stärke der Beleuchtung abhängig.
\end{itemize}

Da $f$ und $D$ linear mit der Lichtintensität wachsen, muss die geschwindigkeitsabhängige Reibung noch sehr
klein sein. In Abbildung \ref{fig:pfdata} scheint $P_\mathrm{opt}$ für die beiden Kurven unterschiedlich zu
sein. Weitergehende Messungen in \cite{duffner} zeigen aber, dass es zwar Schwankungen bei $P_\mathrm{opt}$
gibt,
sich aber kein Trend bei wachsenden Lichtleistungen ableiten lässt.

\begin{figure}
\capstart\centering
\include{pfdata}
\caption{Experimentelle Daten aus \cite{geyer} und \cite{duffner}: Umdrehungen pro Sekunde $f$ der Lichtmühle
in Abhängigkeit des Drucks $P$ in der Kammer für zwei verschieden starke Beleuchtungen}
\label{fig:pfdata}
\end{figure}

Abbildung \ref{fig:Dvolker} zeigt das Drehmoment, das auf die Lichtmühle wirkt, die in dieser Diplomarbeit
simuliert wurde. Außerdem vergleicht sie die Form der Drehmomentkurve mit der Form der Frequenzkurve. Das ist
wichtig, da die Drehmomentkurve in Abbildung \ref{fig:Dvolker} die einzige ist, die zum Vergleich vorliegt,
und die meisten weiteren Zusammenhänge aus Drehfrequenzkurven gewonnen werden mussten. Eine weitere
Drehmoment"=Druckkurve zeigt Abbildung \ref{fig:Ddaniel}. Die Lichtmühle, die in dieser Messung verwendet
wurde, wurde allerdings nicht simuliert.

\begin{figure}
\capstart\centering
\include{Dvolker}
\caption{Drehmoment-Druckverlauf aus \cite{geyer} verglichen mit der Form des Drehfrequenz"=Druckverlaufs der
Lichtmühle, nach der die simulierte Lichtmühle dieser Diplomarbeit modelliert wurde.}
\label{fig:Dvolker}
\end{figure}

\begin{figure}
\capstart\centering
\include{Ddaniel}
\caption{Drehmoment-Druckverlauf aus \cite{duffner} verglichen mit der Form des Drehfrequenz"=Druckverlaufs
einer Lichtmühle, die \textbf{nicht} simuliert wurde.}
\label{fig:Ddaniel}
\end{figure}

Auch wenn die Drehmoment"=Druckverläufe das Maximum bei etwas höheren Drü"-cken als die Frequenzkurven haben,
wurde davon ausgegangen, dass die Verläufe ähnlich sind. Die Unterschiede zwischen den Drehmomentverläufen
und Frequenzverläufen können auch leicht durch Messfehler hervorgerufen werden oder sie würden vielleicht
besser übereinstimmen, wenn Vergleichswerte für größere Druckbereiche vorliegen würden.


\subsubsection{Abschätzung des Temperaturgradienten}

Leider wurden von \cite{geyer} und \cite{duffner} keine Messungen der Temperatur, die auf Lichtmühlenflügel
herrscht, durchgeführt und es wurden auch keine anderen Quellen für solche Daten gefunden. Aus den
experimentell bestimmten Drehmomenten lässt sich mit Gleichung \ref{fi} aber eine Abschätzung der
Temperaturdifferenz vornehmen. Ersetzt man die Scherspannung $f$ in Gleichung \ref{fi} durch $F/A$ und die
Ableitung der Temperatur durch den Differenzenquotienten $\Delta T/d$ und setzt $\sigma_p = 1$ ergibt sich

\begin{equation}
\frac{F}{A} =\frac{75 \sqrt{2}}{1024} \frac{k}{\sigma^2} \frac{\Delta T}{d}.
\label{absch1}
\end{equation}

Weil $A$ in Gleichung \ref{absch1} die Fläche ist auf der in der Lichtmühle thermische Transpiration
stattfindet, gilt $A = 4\cdot 4 ld$. Denn die "`Kantenfläche"' eines Fügels ist gerade vier mal seine Dicke
$d$ mal seine Länge $l$ und des Weiteren gibt es vier Flügel in der Lichtmühle. Außerdem kann man die Kraft 
auf die Lichtmühlenflügel $F$ mit $F=D/a$ durch das gemessene Drehmoment $D$ und der Abstand von Drehachse zu
Flügelschwerpunkt $a$ ersetzten. Wenn man $A$ und $F$ in Gleichung \ref{absch1} so ersetzt und nach $\Delta T$
auflöst ergibt sich

\begin{equation}
\Delta T = \frac{64}{75\sqrt{2}} \frac{D}{al} \frac{\sigma^2}{k}.
\label{absch}
\end{equation}

Interessant ist, dass die Abschätzung von $\Delta T$ in Gleichung \ref{absch} nicht mehr von der Flügeldicke
$d$ abhängt. Setzt man die in diesem Kapitel genannten Abmessungen in Gleichung \ref{absch} ein -- also
$a=2.02$\,cm und $l=2.3$\,cm -- ergibt sich mit dem maximalen Drehmoment $D=0.25\,\text{\textmu}\mathrm{Nm}$
aus
Abbildung \ref{fig:Dvolker} ein $\Delta T$ von 3.2\,K.



\section{Simulationsresultate}

Es folgen nun die Simulationsergebnisse für die in Kapitel \ref{simlicht} beschriebene Lichtmühle mit den
Parametern
aus Tabelle \ref{tab:parameter}. In der Simulation wurde ein Referenzdruck $P$ und eine initiale Gastemperatur
$T = 293\,\mathrm K$ eingestellt, damit gibt die initiale Dichte im Behälter mit der idealen Gasgleichung
\ref{idGas}. Auf die Außenwände und die "`silbrige Seite"', die in der Simulation immer die untere Seite des
Flügels ist, wurde ebenfalls $T= 293\,\mathrm K$ gesetzt. Die "`schwarze Seite"' hatte eine um 20\,K höhere
Temperatur, also 313\,K. Auf die beiden kurzen Seiten wurde eine mit der y-Koordinate linear von 293\,K
auf 313\,K steigende Temperatur gesetzt.
Da keine Informationen über die tatsächliche Temperaturdifferenz $\Delta T$ auf
Lichtmühlenflügeln vorlag,
wurden einfach $\Delta T$-Werte gewählt, die Drehmomente erzeugen, die zumindest in der selben
Größenordnung liegen wie die experimentellen Daten aus Abbildung \ref{fig:Dvolker}.

Abbildung \ref{fig:paddleTemp} zeigt die Temperaturverteilung, die sich im Gleichgewichtszustand und einem
Referenzdruck von 4\,Pa eingestellt hat, und Abbildung \ref{fig:paddleSpeed} die entsprechende
Geschwindigkeitsverteilung. Das Einzige, das das Gas zur Bewegung bringt, sind Temperaturgradienten. Weil aber
Smoluchowskis Randbedinungen dafür sorgen, dass im Gas auch entlang der langen Seiten Temperaturgradienten
existieren, gibt es auch dort Strömungen. Wie es nach der Theorie der thermischen Transpiration zu erwarten
ist, kommt es zu einer Strömung an den Ecken des Flügels in Richtung der wärmeren, oberen Seite. Weil die
Strömung auf dem Rand des Flügels am stärksten ist und mit zunehmender Entfernung vom Flügel abnimmt, ist die
viskose Kraft nach unten, also entgegen der Richtung der Strömung, gerichtet. Das heißt die simulierte
Lichtmühle würde sich wie eine echte mit der silbrigen Seite voran drehen. Die beschriebene Strömung um die
Flügelkannte herum ist in Abbildung \ref{fig:paddleVelocity} zu sehen.

\begin{figure}
\capstart\centering
\includegraphics[width=0.8\linewidth]{fig/trptemp.png}
\caption{Temperaturverteilung (K) um den simulierten Flügel. Wegen Smoluchowskis Randbedingungen erreicht das
Gas an keiner Stelle die 313\,K, die auf die obere Flügelseite gesetzt sind. Der Referenzdruck ist 4\,Pa und
$\Delta T=20$\,K.}
\label{fig:paddleTemp}
\end{figure}
\begin{figure}
\capstart\centering
\includegraphics[width=0.8\linewidth]{fig/trpspeed.png}
\caption{Betrag der Geschwindigkeiten (m/s) für die gleiche Simulation wie in Abbildung \ref{fig:paddleTemp}}
\label{fig:paddleSpeed}
\end{figure}
\begin{figure}
\capstart\centering
\includegraphics[width=1.0\linewidth]{fig/velvecTJ.png}
\caption{Vergrößerter Ausschnitt der rechten Flügelkante aus Abbildung \ref{fig:paddleSpeed}. Die Pfeile
stehen für die Geschwindigkeitsvektoren. Nur ihre Länge ist proportional zum Betrag der Geschwindigkeit; ihre
Häufigkeit hängt von der Knotendichte ab.}
\label{fig:paddleVelocity}
\end{figure}

Abbildung \ref{fig:v_x} und \ref{fig:v_y} zeigen die x- beziehungsweise y"=Komponente der Gleitgeschwindigkeit
auf dem
Flügelrand. Man kann in ihnen deutlich besser die Werte und Verteilung der Geschwindigkeiten sehen als in
Abbildung \ref{fig:paddleSpeed} oder \ref{fig:paddleVelocity}. Natürlich ist Geschwindigkeitskomponente, die
senkrecht auf dem Rand liegt immer 0, da dies die Randbedingung einer undurchdringlichen Wand ist.

\begin{figure}
\capstart\centering
\include{v_x}
\vspace{-2\baselineskip}
\caption{x"=Komponente der Geschwindigkeit $u_x$ (Punkte) auf der rechten Kante des Flügels (durchgezogene
Linie) für die selbe Simulation wie in Abbildung \ref{fig:paddleTemp}.}
\label{fig:v_x}
\end{figure}

\begin{figure}
\capstart\centering
\include{v_y}
\vspace{-2\baselineskip}
\caption{y"=Komponente der Geschwindigkeit $u_y$ (Punkte) auf der rechten Kante des Flügels (durchgezogene
Linie) für die selbe Simulation wie in Abbildung \ref{fig:paddleTemp}.}
\label{fig:v_y}
\end{figure}

\subsection{Viskose Kräfte auf dem Flügel}
\label{kraft}

Die auf den Flügel wirkende viskosen Kräfte sind das wichtigste Detail der Simulation, da mit ihnen das
Drehmoment, also die Vergleichsgröße mit den Experimenten bestimmt wird. In anfänglichen Testrechnungen gab es
große Probleme mit den simulierten Kräften. Der Flügel wurde dort zunächst mit einem Rechteck ohne
abgerundete Ecken modelliert, was zu sehr starken Ausreißern in den Kraftwerten $\bm F_i$ auf den vier
Eckknoten führte. Die Ausreißer waren auflösungsabhängig und damit sicher unphysikalisch. Weil als Grund für
ihr Auftreten die Nicht-Differenzierbarkeit von Ecken vermutet wurde, wurde versucht durch Abrunden der Ecken
dieses Problem zu Lösen.

Abbildung \ref{fig:sp_x} und \ref{fig:sp_y} zeigen die x- beziehungsweise y"=Komponente der Kraft pro Länge
($F/l$) auf
einen kleinen Ausschnitt der linken, schon abgerundeten Flügelkante, die als durchgezogene Linie ebenfalls zu
sehen ist. Man erkennt klar, dass die größten Kräfte auf die Kante wirken und auf den langen Flügelseiten
abfallen. Wegen des Kreuzproduktes $\bm D= \bm r \times \bm F$ und der Position der Drehachse, hat die
x"=Komponente von $\bm F$ nur Auswirkungen im Promillbereich auf das resultierende Drehmoment $D$.

\begin{figure}
\capstart\centering
\include{sp_x}
\vspace{-2\baselineskip}
\caption{x"=Komponente der Kraft pro Länge $F_x/l$ (Punkte) auf der linken Kante des Flügels
(durchgezogene Linie) in willkürlichen Einheiten für die selbe Simulation wie in Abbildung
\ref{fig:paddleTemp}.}
\label{fig:sp_x}
\end{figure}

\begin{figure}
\capstart\centering
\include{sp_y}
\vspace{-2\baselineskip}
\caption{y"=Komponente der Kraft pro Länge $F_y/l$ (Punkte) auf der linken Kante des Flügels
(durchgezogene Linie) in willkürlichen Einheiten für die selbe Simulation wie in Abbildung
\ref{fig:paddleTemp}.}
\label{fig:sp_y}
\end{figure}

Die negativen $F_y$-Werte in Abbildung \ref{fig:sp_y} sind der wichtigste Beitrag zum Drehmoment und würden
im Fall eines beweglichen Flügels für eine Bewegung in negative y"=Richtung sorgen. Zum Teil werden sie durch
die positiven $F_y$-Werte auf den runden Ecken des Flügels wieder kompensiert.

Man erkennt in Abbildung \ref{fig:sp_x}, dass es dort eine Spitze im $F_x$"=Verlauf gibt, wo das gerade Stück
 der kurzen Kante in die Abrundungen der Ecken übergeht. Analog dazu gibt es auch Spitzen im $F_y$"=Verlauf in
Abbildung \ref{fig:sp_y} an den Stellen, wo die langen Seiten in die Viertelskreise münden. Wie auch schon
die nicht gezeigten Spitzen in den Kraftverläufen auf den Ecken von perfekt rechteckigen Flügeln, wachsen
auch Spitzen auf den abgerundeten Flügeln mit der Auflösung, wenn auch deutlich schwächer. Die stärkere
Ausprägung der Spitzen sieht man deutlich in Abbildung \ref{fig:sp_xH} und \ref{fig:sp_yH}, die den
Kraftverlauf in eine sehr hohe Auflösung zeigen. Abbildung \ref{fig:sp_xH} und \ref{fig:sp_yH} zeigt einen
Flügel mit einer anderen Dicke $d$ als die bisher gezeigten und wurden ohne Temperatursprungrandbedingungen
erzeugt, was die Spitzen deutlicher hervortreten lässt. Die Spitzen wuchsen aber auch in den mit der
bisherigen Konfiguration getätigten Simulationen. In Kapitel \ref{ohneTJ} wird genauer auf Simulationen ohne
Smoluchowskis Randbedingungen eingegangen.

\begin{figure}
\capstart\centering
\include{sp_xH}
\vspace{-2\baselineskip}
\caption{x"=Komponente der Kraft pro Länge $F_x/l$ (Punkte) auf der linken Kante des Flügels
(durchgezogene Linie) in willkürlichen Einheiten bei sehr hoher Elementdichte zur Verdeutlichung der Spitzen
in den $F_x$-Werten.}
\label{fig:sp_xH}
\end{figure}

\begin{figure}
\capstart\centering
\include{sp_yH}
\vspace{-2\baselineskip}
\caption{y"=Komponente der Kraft pro Länge $F_y/l$ (Punkte) auf der linken Kante des Flügels
(durchgezogene Linie) in willkürlichen Einheiten bei sehr hoher Elementdichte zur Verdeutlichung der Spitzen
in den $F_y$-Werten.}
\label{fig:sp_yH}
\end{figure}

Übergänge von Kreisabschnitten in ihre Tangenten sind nur einmal stetig differenzierbar. Weil die
Navier-Stokes-Gleichungen Differenzialgleichungen zweiter Ordnung sind, kann es sein, dass dies der Grund für
das Auftreten der Spitzen an den Übergängen von geraden und kreisförmigen Abschnitten des Flügels ist.
Leider konnte dieser Sachverhalt aus Zeitmangel nicht mehr eindeutig geklärt werden. In jedem Fall
führen die Spitzen in $\bm F$, solange ihre Ursache nicht bekannt ist, zu Unsicherheiten in allen weiteren
Ergebnissen.

Die Abrundung der Ecken brachte dennoch eine Verbesserung gegenüber den ursprünglich verwendeten und hier
nicht gezeigten rechtwinkligen Ecken. So sieht man in Tabelle \ref{tab:konvergenz}, dass das resultierende
Drehmoment $D$ mit steigenden Knotenzahlen $N$ zu konvergieren scheint, auch wenn die Spitzen in den Kräften
mit $N$ größer werden. Diese Konvergenz von $D$ war bei Simulationen, die rechtwinklige Ecken verwendeten
deutlich schlechter. $N$ ist in Tabelle \ref{tab:konvergenz} nicht die gesamte Knotenzahl, sondern die Zahl
der
Knoten auf dem Flügelrand. Die Simulationseinstellungen sind in Tabelle \ref{tab:konvergenz} die gleichen wie
in Abbildung \ref{fig:sp_xH} und \ref{fig:sp_yH}; diese Abbildungen sind mit den Daten der Simulation mit
$N=2840$ erstellt worden.

\begin{table}
\capstart\centering
\begin{tabular}{l|l}

Knoten $N$ & Drehmoment $D$ (Nm)\\
\hline
696 & $1.47025518\cdot 10^{-7}$\\
\hline
1070 & $1.45789741\cdot 10^{-7}$\\
\hline
1342 & $1.45618222\cdot 10^{-7}$\\
\hline
1772 & $1.45569079\cdot 10^{-7}$\\
\hline
2840 & $1.45559889\cdot 10^{-7}$\\
\end{tabular}
\caption{Abhängigkeit des resultierenden Drehmoments $D$ von der Zahl der Knoten auf dem Flügelrand $N$. Man
sieht, dass $D$ relativ konstant bleibt, obwohl die Spitzen in den Kräften auf dem Flügel mit $N$ wachsen.}
\label{tab:konvergenz}
\end{table}

\subsection{Simulierte Druck"=Drehmomentkurven}


Wenn die Simulation für ein $P$ den Gleichgewichtszustand erreicht hatte, wurde das Drehmoment $D$ aus den
viskosen Kräften auf den Flügel berechnet und die Simulation mit einem neuen $P$ gestartet. Die sich daraus
ergebenden Druck"=Drehmomentkurven bieten Vergleichsmöglichkeiten mit den Experimenten.

Abbildung \ref{fig:Dp} und \ref{fig:Dp2} zeigen die Variation des Betrags des Drehmoments $D$ mit dem
Referenzdruck $P$ für zwei verschiedene Druckbereiche. Dabei zeigt sich zunächst, dass die simulierten Kurven
die gleichen qualitativen Eigenschaften wie die gemessenen haben: Für $P\rightarrow 0$ und $P\rightarrow
\infty$ geht $D$ gegen 0. Der Abfall des Drehmoments für große $P$ ist exponentiell und es gibt ein Maximum
des Drehmoments $D_\mathrm{max}$ bei niedrigen Drücken, nämlich bei 4\,Pa.

\begin{figure}
\capstart\centering
\include{tD20K}
\caption{Drehmoment $D$ auf die simulierte Lichtmühle in Abhängigkeit des Drucks $P$ bei $\Delta T = 20$\,K}
\label{fig:Dp}
\end{figure}

\begin{figure}
\capstart\centering
\include{tD20K4000}
\caption{Drehmoment $D$ auf die simulierte Lichtmühle in Abhängigkeit des Drucks $P$ bei $\Delta T = 20$\,K}
\label{fig:Dp2}
\end{figure}
 
Quantitativ gibt es allerdings deutliche Unterschiede zwischen Simulation und Experiment. So ist der Druck,
bei
 dem das maximale Drehmoment auftritt, in der Simulation doppelt so hoch wie im Experiment. Der Abfall von
$D$
bei steigendem $P$ fällt in der Simulation deutlich schwächer aus. Der Vergleich ist zunächst nicht so
einfach, weil, wie man in Abbildung \ref{fig:Dvolker} sehen kann, nur Messungen bis 4\,Pa vorliegen. Aus
Abbildung \ref{fig:Ddaniel}, die an einer anderen Lichtmühle gemessen wurde, kann man aber schließen, dass der
Abfall des Drehmomentes zumindest ähnlich wie der Abfall der Drehfrequenz $f$ ist, und diese zeigen, wie man
in
Abbildung \ref{fig:pfdata} sehen kann, einen sehr viel schnelleren Abfall bei steigendem Druck. Des Weiteren
ist das maximale Drehmoment bei $\Delta T=20$\,K zu klein. Erst bei ca. 40\,K Temperaturdifferenz stimmen die
maximalen Drehmomente überein. Diese Temperatur ist eine ganze Größenordnung höher als das Resultat der
Abschätzung in Gleichung \ref{absch}, welches $\Delta T=3.2$\,K ergab. Damit scheinen die Drehmomente, die die
Simulation produziert, zu niedrig zu sein. Aber ohne das Wissen über die tatsächlichen Temperaturen auf einem
Lichtmühlenflügel ist es nicht möglich, eine klare Aussage über die Güte der Simulation in Bezug auf das
maximale Drehmoment zu treffen.




\subsection{Variation der Flügelform}

Als nächstes wurden Drehmoment"=Druckkurven für verschiedene Flügelformen untersucht. Dazu wurden zwei weitere
Gitter mit je verschiedenen Fügeldicken $d=0.25$\,mm und $d=1$\,mm erzeugt. Der Radius der Viertelskreise in
den vier Ecken des $d=0.25$\,mm Flügels war wie beim $d=0.15$\,mm Flügel ein sechstel seiner Dicke, beim
$d=1$\,mm Flügel dagegen wurde der gleiche Radius $r$ wie beim $d=0.25$\,mm Flügel benutzt -- also
$r=41.\bar 6\,\text{\textmu} \mathrm m$. Die Länge $l$ war bei allen drei Flügeln gleich, also 2.3\,cm.

In Abbildung \ref{fig:verglform} werden die Drehmoment"=Druckkurven für die drei verschiedenen Flügel
verglichen. Man sieht, dass sich unterschiedliche Kurvenverläufe ergeben. Insbesondere wird das maximale
Drehmomente $D_\mathrm{max}$ bei verschiedenen Drücken erreicht und nimmt unterschiedliche Werte an.
Außerdem ist der Abfall der Kurven bei großen $P$ unterschiedlich
stark. Die genauen Werte von $D_\mathrm{max}$ kann man Tabelle \ref{tab:Popt} entnehmen. Es ist
interessant, dass der $d=1$\,mm Flügel zwar das größte maximale Drehmoment hat, $D$ aber so schnell für hohe
$P$ abfällt, dass er bei 100\,Pa kleinere Drehmomente als die anderen beiden Flügel erzeugt.

\begin{figure}
\capstart\centering
\include{verglform}
\caption{Drehmoment-Druckverlauf bei $\Delta T = 20$ für alle drei Flügeldicken $d$. Unterschiedliche $d$
führen zu deutlich unterschiedlichen Kurvenverläufen.}
\label{fig:verglform}
\end{figure}

Zwar haben die gemessenen Drehmoment"=Druckkurven von zwei verschiedenen Lichtmühlen in Abbildung
\ref{fig:Dvolker} und \ref{fig:Ddaniel} auch unterschiedliche Verläufe, man kann aber die Unterschiede
dennoch nicht mit der Simulation vergleichen, da sich die zwei Lichtmühlen nicht nur in der Dicke der 
Flügel sondern auch in der Länge und vor allem auch im Material unterscheiden, was unterschiedliche TMAC- und
TAC"=Werte zur Folge haben kann.



\subsection{Variation der Temperaturdifferenz}


Als nächstes wurden die Temperaturdifferenzen $\Delta T$ variiert, in dem die Temperatur der warmen, oberen
Flügelseite geändert und alle anderen Temperaturen wie zuvor auf $T=293$\,K festgehalten wurden. Wenn
ansonsten alle Parameter gleich bleiben, führen höhere $\Delta T$ zu einer Drehmoment"=Druckkurve, deren
Drehmomente um einen Faktor größer sind, der proportional zur Änderung des $\Delta T$-Wertes ist. Das kann man
in Abbildung \ref{fig:vergldT} gut sehen. Das maximale
Drehmoment wächst also linear mit $\Delta T$ und wird immer beim gleichen Druck $P_\mathrm{opt}$
erreicht. Diese beiden Sachverhalte sind in Tabelle \ref{tab:Popt} und in Abbildung \ref{fig:deltaT-Dmax}
dargestellt.

\begin{figure}
\capstart\centering
\include{trp11-5-10-15}
\caption{Drehmoment-Druckverlauf für drei verschiedene Temperaturdifferenzen $\Delta T$. Die Drehmomente $D$
skalieren linear mit $\Delta T$.}
\label{fig:vergldT}
\end{figure}

\begin{table}
\capstart\centering
\begin{tabular}{l|l|l|l|l|l}
\multicolumn{2}{c|}{$d=0.15$\,mm} & \multicolumn{2}{|c|}{$d=0.25$\,mm} & \multicolumn{2}{|c}{$d=1$\,mm}\\
\hline
$\Delta T$\,(K) & $P_\mathrm{opt}$\,(Pa) & $\Delta T$\,(K) & $P_\mathrm{opt}$\,(Pa) &$\Delta T$\,(K) &
$P_\mathrm{opt}$\,(Pa)\\
\hline
5 & 4 & 5 & 7 & 5 & 6.5\\
\hline
10 & 4 & 10 & 7 & 10 & 6.5\\
\hline
15 & 4 & 15 & 7 & 15 & 6.5\\
\hline
20 & 4 & 20 & 7 & 20 & 6.5\\
\hline
50 & 4 & 50 & 6.5 & 50 & 6.5\\
\hline
100 & 4 & 100 & 6.5 & 100 & 6.5\\
\end{tabular}
\caption{Der Druck $P_\mathrm{opt}$, bei dem das maximale Drehmoment $D_\mathrm{max}$ erreicht wurde, für
simulierten Fügeldicken $d$ und Temperaturdifferenzen $\Delta T$. Man sieht, dass $P_\mathrm{opt}$ von der
Flügelform jedoch nicht von $\Delta T$ abhängt.}
\label{tab:Popt}
\end{table}

\begin{figure}
\capstart\centering
\include{deltaT-Dmax}
\caption{Zusammenhang zwischen gesetzter Temperaturdifferenz $\Delta T$ und dem maximal erreichten
Drehmoment $D_\mathrm{max}$ für alle drei Flügeldicken $d$}
\label{fig:deltaT-Dmax}
\end{figure}

Bei dem $d=0.25$\,mm Flügel scheint $P_\mathrm{opt}$ für große $\Delta T$ in Tabelle
\ref{tab:Popt} kleiner zu werden. Die Änderung von $P_\mathrm{opt}$ ist allerdings sehr klein. Wahrscheinlich
ist es nur ein Resultat des Samplings des Drucks bei der Parameterstudie, das 0.5\,Pa war. Es ist zu
vermuten, dass $P_\mathrm{opt}$ des $d=0.25$\,mm Flügels eigentlich zwischen 7 und 6.5\,Pa liegt und
wahrscheinlich würden feinere Parameterstudien dies bestätigen.

Das linare Wachstum von $D_\mathrm{max}$ mit $\Delta T$ und die Konstanz von $P_\mathrm{opt}$ sind in gutem
Einklang mit den experimentellen Daten.



\subsection{Variation der Akkomodationskoeffizienten}

Weil die Akkomodationskoeffizienten $\sigma_p$ und $\sigma_T$ in den Experimenten nicht bekannt sind und die
Literatur von hohen Werten ausgeht, wurden sie in allen bisherigen Simulationen auf 1 gesetzt. In diesem
Abschnitt wird kurz auf ihren Einfluss auf die simulierten Drehmomente eingegangen. Aufgrund der Theorie
erwartet man
ein Absinken der durch thermische Transpiration verursachten Strömung und damit des wirkenden Drehmoments $D$,
weil der Effekt auf einem Impulsaustausch mit der Wand basiert und $\sigma_p$ ja gerade der Anteil der
Teilchen ist, bei denen ein Impulsaustausch mit der Wand stattfindet. Abbildung \ref{fig:ac} zeigt deutlich,
dass kleinere Akkomodationskoeffizienten zu kleineren $D$ führt. Weiterhin zeigt Abbildung \ref{fig:ac}, dass
es zu einem Abflachen der gesamten Drehmomentkurve kommt und $D_\mathrm{max}$ sich zu größeren Drücken hin
verschiebt. Diese beiden Konsquenzen sind nicht direkt aus der Definition der Akkomodationskoeffizienten
ersichtlich.

\begin{figure}
\capstart\centering
\include{ac}
\caption{Der Einfluss der Akkomodationskoeffizienten $\sigma_p$ und $\sigma_T$ auf den Verlauf der
Druck"=Drehmomentkurve für den $d=0.25$\,mm Flügel bei $\Delta T=20$\,K.}
\label{fig:ac}
\end{figure}

Zur besseren Übersicht enthält Abbildung \ref{fig:ac} nur die Ergebnisse des $d=0.25$\,mm Flügel bei
$\Delta T=20$\,K. Die Auswirkungen von kleinen Akkomodationskoeffizienten waren bei den anderen Flügeln und
Temperaturdifferenzen gleich wie in dem einen gezeigten Beispiel.




\section{Ergebnisse ohne Smoluchowskis Randbedingungen}
\label{ohneTJ}

In diesem Abschnitt werden die Ergebnisse von Simulationen ohne Smoluchowskis Randbedingungen vorgestellt,
also Simulationen, die Dirichlet"=Randbedingungen für die Wärmeleitgleichung benutzen. Das soll
veranschaulichen, welche Effekte in den Ergebnissen von welcher Randbedingung -- also Maxwells und
Smoluchowskis
-- erzeugt werden, denn wegen der Ähnlichkeit ihrer Herleitungen ist anzunehmen, dass die Annahmen auf denen
sie basieren, grundsätzlich gleichzeitig auftreten oder gar nicht.

Zunächst ist natürlich die Temperaturverteilung eine andere. In Abbildung \ref{fig:paddleTemp2} sieht man,
dass im Vergleich zu Abbildung \ref{fig:paddleTemp} das Gas an der Wand nun genau die Temperatur der Wand hat
und
es Temperaturgradienten im Gas nur noch auf den kurzen Kanten und nicht mehr auf den langen Seiten gibt. Auch
liegt jetzt der volle Gradient von 20\,K auf den Kanten, während in Abbildung \ref{fig:paddleTemp} zwischen
keinen zwei Punkten im Gas die vollen 20\,K Differenz vorhanden sind. Das hat zur Folge, dass auch nur noch
auf den beiden kleinen Kanten Strömungen induziert werden. Diese haben allerdings höhere Beträge der
Geschwindigkeit $u$ als die Strömungen mit Smoluchowskis Randbedingung, weil der im Gas
herrschende Temperaturgradient größer ist. Die veränderte Form und Stärke der Strömung ist in Abbildung
\ref{fig:paddleSpeed2} und \ref{fig:paddleVelocity2} zu sehen. Die stärkeren Strömungen führen entsprechend
auf größere viskose Kräfte auf die Flügel, die auch auf die beiden kleinen Kanten beschränkt sind. Das ist in
Abbildung \ref{fig:sp_xH} und \ref{fig:sp_yH} zu sehen, wenn auch für einen Flügel mit $d=0.25$\,mm anstatt
$d=0.15$\,mm wie die anderen Abbildungen in diesem Abschnitt.

\begin{figure}
\capstart\centering
\includegraphics[width=0.8\linewidth]{fig/trpp4T.png}
\caption{Temperaturverteilung (K) um den simulierten Flügel mit Dirichlet"=Randbedingungen für die
Wärmeleitgleichung aber ansonsten die gleichen Bedingungen wie in Abbildung \ref{fig:paddleTemp}.}
\label{fig:paddleTemp2}
\end{figure}
\begin{figure}
\capstart\centering
\includegraphics[width=0.8\linewidth]{fig/trpp4v.png}
\caption{Betrag der Geschwindigkeit (m/s) für die gleiche Simulation wie in Abbildung
\ref{fig:paddleTemp2}. Es ist nicht das gesamte Gitter sondern nur ein Ausschnitt um den Flügel abgebildet.}
\label{fig:paddleSpeed2}
\end{figure}
\begin{figure}
\capstart\centering
\includegraphics[width=1.0\linewidth]{fig/velvec.png}
\caption{Vergrößerter Ausschnitt der linken Flügelkante aus Abbildung \ref{fig:paddleSpeed2}. Die Pfeile
stehen für die Geschwinkdigkeitsvektoren. Nur ihre Länge ist proportional zum Betrag der Geschwindigkeit; ihre
Häufigkeit hängt von der Knotendichte ab.}
\label{fig:paddleVelocity2}
\end{figure}


Nach dem Gesagten ist klar, dass die Drehmomente in einer Simulation ohne Smoluchowskis Randbedingung größer
ausfallen müssen. Abbildung \ref{fig:mitohneTJ} vergleicht die Druck"=Drehmomentkurven von Simulationen mit
und ohne Smoluchowskis Randbedingung und zeigt deutlich die höheren $D$. Weil Smoluchowskis Randbedingung für
große $\lambda$ mehr und mehr in eine "`normale"' Dirichlet"=Randbedingung übergeht, nähern sich die beiden
Kurven für große $P$ immer mehr an und liegen, was in Abbildung \ref{fig:mitohneTJ} nicht mehr zu sehen ist,
ab ca. 4000\,Pa aufeinander.

\begin{figure}
\capstart\centering
\include{mitohneTJ}
\caption{Drehmoment $D$ in Abhängigkeit des Referenzdrucks $P$. Es wird die Kurve aus Abbildung \ref{fig:Dp}
mit einer Simulation verglichen, bei der alle Parameter bis auf das Aktivieren von Smoluchowskis
Randbedingungen
gleich sind, also $d=0.15$\,mm und $\Delta T=20$\,K. Ab ca. 4000\,Pa liegen die beiden Kurven aufeinander.}
\label{fig:mitohneTJ}
\end{figure}

Man sieht in Abbildung \ref{fig:mitohneTJ} deutlich, dass das Drehmoment $D$ für $P\rightarrow 0$ immer
weiter steigt und im Gegensatz zum Experiment nicht auf 0 abfällt. Somit ist klar, dass das "`verschmieren"'
des Temperaturgradienten durch Smoluchowskis Randbedingungen der einzige Grund für das Abfallen von $D$ in
allen anderen gezeigten Druck"=Drehmomentkurven ist. Wie man in Gleichung \ref{maxwell} sehen kann, wird die
thermische Transpiration immer stärker mit abnehmender Dichte $\rho$; daher wird auch die viskose Kraft immer
größer, da sie nur von Geschwindigkeitsgradienten abhängt und nicht vom Massenfluss, der natürlich immer
kleiner wird.

In der Natur gibt es auch ohne den Temperatursprung in verdünnten Gasen einen Mechanismus, der die durch
thermische Transpiration hervorgerufenen Kräfte am Divergieren für $P\rightarrow 0$ hindert. Das in dieser
Diplomarbeit verwendete Modell für die dynamische Viskosität (Gleichung \ref{viskT}) wird für extrem verdünnte
Gase, also für $\kn \gg 1$, ungültig. Die Wechselwirkung zwischen einem Gas und einer Oberfläche wird dann
besser ausgedrückt durch
\begin{equation}
Z = \frac{\sigma_p}{2-\sigma_p} \frac{P}{\sqrt{2\uppi\Rs T}} = \frac{\sigma_p}{2-\sigma_p} \rho
\sqrt{\frac{T\Rs}{2\uppi}}.
\label{molvisk}
\end{equation}
Gleichung \ref{molvisk} definiert die "`Viskosität der freien Moleküle"' $Z$, wie sie zum Beispiel von
\cite{kennard} hergeleitet wird. Die Kraft $F_x$, die durch eine Strömung freier Moleküle in x"=Richtung $u_x$
über eine weit ausgedehnte Fläche $A$, die in der x"=z"=Ebene liegt, erzeugt wird, ist $F_x =Zu_x A$. In
dieser Gleichung sieht man, dass $Z$ die Dimension $\mathrm{Pa\,s/m}$ hat, die verschieden ist von der
Dimension $\mathrm{Pa\,s}$ von $\mu$. Bei der Berechnung von $F_x$ ist kein Geschwindigkeitsgradient mehr
enthalten, weil in dem Regime in dem $Z$ gültig ist, Stöße zwischen Molekülen so selten sind, dass eine
gasinterne Viskosität vernachlässigbar ist beziehungsweise das Konzept einer viskosen Strömung also nicht mehr
anwendbar
ist.
Weil $Z$ für kleine $\rho$ gegen 0 strebt, schwächt das die thermische Transpiration selbst (vergleiche
Gleichung \ref{maxwell}) als auch der viskosen Kraft auf die Flügel bei sehr kleinen $\rho$ stark ab. Da
aber Maxwell selbst bei der
Herleitung der thermische Transpiration von "`normaler"' Viskosität ausging, sind die implementierten Modelle
konsistent.

Der Vergleich von Ergebnissen mit und ohne Smoluchowskis Randbedingung macht noch einmal den großen Einfluss
von Randeffekten deutlich und wirft die Frage auf, ob die implementieren Modelle noch ausreichend sind für den
in
dieser Dipomarbeit untersuchten Druckbereich.

\chapter{Schlussfolgerungen und Ausblick}

\section{Was erreicht wurde}

Im Rahmen dieser Diplomarbeit wurden spezielle Randbedingungen zur Simulation von Effekten in verdünnten
Gasen -- namentlich Maxwells Thermische"=Transpirationsrandbedingung und Smoluchowskis
Temperatursprungrandbedingung -- zur Verwendung mit dem Finite"=Elemente"=Löser Elmer
implementiert. Die Implementierung wurde erfolgreich an Simulationsergebnissen anderer Programme getestet. Sie
ist in ihrer momentanen Form auf beliebig geformte 2D"=Geometrien anwendbar und durch ihren modularen Aufbau
leicht auf 3D"=Geometrien und weitere Effekte erweiterbar.

Des Weiteren wurde die Simulation genutzt um die auf eine Lichtmühle wirkenden Drehmomente zu errechnen. Die
simulierten Druck"=Drehmomentkurven haben zumindest ähnliche Charakteristiken wie experimentell bestimmte.
Zwei Eigenschaften der numerischen Resultate -- die Konstanz des optimalen Drucks der Lichtmühle bei
verschiedenen Leistungen und das lineare Wachstum der gesamten Druck"=Drehmomentkurve mit der eingestrahlten
Leistung -- sind sogar in sehr guter Übereinstimmung mit den Experimenten.


\section{Noch vorhandene Probleme}

Es gibt auch noch viele Unstimmigkeiten bei der numerischen Simulation von Lichtmühlen. Der wichtigste Punkt
sind die in Kaptitel \ref{kraft} aufgezeigten unphysikalischen Spitzen in den Kräften auf den Flügeln.
Solange ihre Ursachen nicht geklärt ist, bleiben Zweifel ob die Implementierung korrekt ist.

Die folgende Liste nennt die vorhandenen Abweichungen zwischen Simulation und Experimente:
\begin{itemize}
\item
Wenn die Abschätzung der nötigen Temperaturdifferenz $\Delta T$ aus Kapitel \ref{exp} eine realistische
Größenordung liefert, sind die von der Simulation errechneten Drehmomente um eine Größenordnung zu klein.
\item
Der Abfall der simulierten Druck"=Drehmomentkurven für große $P$ ist viel geringer als der Abfall der
Drehfrequenzen $f$ im Experiment.
\end{itemize}

Diese beiden Abweichungen könnten auch verknüpft sein: Wenn die Simulation für sehr große $P$ richtige
Ergebnisse liefert, das maximale Drehmomente aber eine Größenordnung höher sein müsste, wäre der Abfall der
Kurve für steigende $P$ automatisch stärker.

\section{Gründe für Abweichungen vom Experiment und Verbesserungsmöglichkeiten}

Es gibt eine Vielzahl von Gründen unterschiedlicher Natur für die Abweichungen zwischen Theorie und
Experiment.

Als erstes sollte die Ursache für die unphysikalischen Spitzen in den Kräften auf den Flügeln gefunden
werden. Ob, wie in Kapitel \ref{exp} erwähnt, die Ursache in der nur einmaligen Differenzierbarkeit von
Übergängen eines Kreisabschnitts in seine Tangente liegt, kann durch ein neues geometrisches Modell für den
Flügel geklärt werden. Modelliert man den Flügel anstatt mit Kreisen und Geraden mit kubischen Splines, die
zweimal stetig differenzierbar sind, wäre diese Frage geklärt. Das ist im Prinzip mit dem für die
Gittergenerierung verwendeten Programm -- Gmsh -- ohne Probleme machbar; konnte aber leider nicht mehr im
Rahmen
dieser Arbeit getestet werden.

Es ist wichtig deutlich zu machen, dass die beiden im letzen Abschnitt genannten Abweichungen zwischen
Simulation und Experiment zum Teil auch durch zu wenige Messdaten versucht sein können. Die Vermutung, dass
die maximalen Drehmomente zu klein sind basiert nur auf einer Abschätzung der Temperaturen auf den
Lichtmühlenflügeln, es
liegen leider keine experimentellen Daten zu diesen Temperaturen vor. Auch in der Literatur wurden keine
gemessenen Flügeltemperaturen gefunden. Der zu schwache Abfall der Druck"=Drehmomentkurve für hohe $P$, wird
vor allem aus dem Vergleich mit Drehfrequenzdaten abgeleitet. Dass diese sich gleich verhalten wie die
Drehmomente ist zwar plausibel, aber aus den Messdaten nicht eindeutig abzulesen. Man sieht in Abbildung
\ref{fig:Dvolker} und \ref{fig:Ddaniel}, dass nur zwei beziehungsweise drei Drehmomentwerte rechts vom
jeweiligen
Maximum stehen. Es ist schwierig mit so wenig Daten auf einen gleichen Verlauf von Drehmoment und Drehfrequenz
zu schließen.

Nun zu den Gründen, warum die bisherige Implementierung wahrscheinlich noch nicht ausreichend ist. Die zwei
wichtigsten sind sicher:
\begin{itemize}
\item
Die Navier"=Stokes"=Gleichungen sind nicht exakt genug in dem untersuchten Regime. Die Knudsenzahl liegt
in der Lichtmühle sicher um 1 und für die sehr kleine Drücke $P<1$\,Pa gilt $\kn > 10$.  Damit sehr eng
verwandt
ist der nächste Punkt:
\item
Es wurden zu wenige Effekte berücksichtigt, die in diesem Regim eine Rolle spielen. Zwei wurden in dieser
Diplomarbeit erwähnt -- die mit $\rho$ kleiner werdende Viskosität und der "`thermal stress slip flow"'.
\end{itemize} 
Diese zwei Punkte sind miteinander verbunden, weil die natürlichen Randbedingungen der einer Erweiterung die
Navier"=Stokes"=Gleichungen, also der Brunett"=Gleichungen, Randbedingungen zweiter Ordnung sind, die Effekte
wie den "`thermal stress slip flow"' enthalten. Natürlich kann man auch versuchen zusätzliche Effekte in die
Navier"=Stokes"=Gleichungen einzubauen, das wäre einfacher und würde weniger Rechnerresourcen benötigen,
könnte aber auch wiederum nicht ausreichend sein. Ein deutlicher Hinweis auf die Unzulänglichkeit der
aktuellen Implementierung bei hohen $\kn$ liefert der Vergleich mit den Ergebnissen einer DSMC"=Simulation 
von \cite{karniadakis2005man} in Tabelle \ref{tab:Testcase}. Darin gibt es schon deutliche Abweichungen
für $\kn = 0.365$. Allerdings ist das natürlich kein Vergleich mit experimentellen Daten und man weiß
nicht, wie gut \cite{karniadakis2005man} seine Simulation mit Experimenten verglichen hat.

Es gibt eine viele weitere Verbesserungsmöglichkeiten der momentanen Implementierung. Die meisten davon wären
weniger aufwändig zu testen, wie die beiden schon genannten. Einige Möglichkeiten sind:
\begin{itemize}
\item 
3D"=Rechnungen anstatt 2D"=Rechnungen verwenden. Weil die Flügelform deutliche Auswirkungen auf den Verlauf
der Druck"=Drehmomentkurven hatte, ist anzunehmen, dass auch dieser Schritt drauf einen Einfluss hat.
\item
Flügelgeometrie besser modellieren. Wenn man sich die Lichtmühlenflügel genau anschaut, sieht man, dass sie
nicht
glatt sind, was eine größere Oberfläche als die in der Simulation benutzen zur Folge hat. Auch das kann die
Druck"=Drehmomentkurven beeinflussen.
\item
Eine vollständige Lichtmühle mit vier Flügeln und rundem Behälter simulieren. Ein solches 2D"=Gitter wurden
schon erstellt und erste Testrechnungen legen nahe, dass der Einfluss auf die Drehmomente sehr klein ist,
allerdings wurden noch keine Parameterstudien durchgeführt.
\end{itemize}

Natürlich kann auch die Tatsache, dass nur statische Flügel und keine sich bewegenden simuliert wurden eine
Rolle spielen. Dynamische Flügel können aber nicht als eine Erweiterung der bisherigen Implementierung
realisiert werden, sondern nur als ein ganz neues Projekt. Es besteht auch immer die Gefahr, dass sich Fehler
in die Implementierung eingeschlichen haben, wie zum Beispiel beim Setzen der Randbedingungen. Es hat sich bei
Testrechnungen gezeigt, dass selbst wenn nur einzelne Knoten -- vier auf dem Ganzen Flügel -- andere
Randbedingungen bekommen, die Druck"=Drehmomentkurven andere Verläufe haben können.

\section{Schlusswort}
Abschließend ist zu sagen, dass der Versuch Lichtmühlen zu modellieren vielleicht etwas zu ambitioniert war,
 weil sie bei recht hohen $\kn$"=Werten operieren. Es wäre sehr interessant die Implementierung mit
"`einfacheren"' Experimenten, wie zum Beispiel Strömungen durch klar definierte Kanülen im Bereich $\kn \leqq
0.1$, zu vergleichen.

Andererseits wurde trotz ausgiebiger Recherche nur eine einzige weitere Arbeit zu numerischen
Simulation von Lichtmühlen gefunden. Allerdings benutzt dieser Artikel \cite{scheisse} sehr unrealistische
Parameter wie eine Tempereraturdifferenz von 600\,K und ein"=dimensionale Flügel ($d=0\,\mathrm m$). In diesem
Sinn stellt diese Diplomarbeit eine Verbesserung auf einem bisher kaum bearbeiteten Gebiet dar.

\chapter{Verbeserungen}

In diesem Abschnitt liste ich Verbesserungen von inhaltlichen Fehlern auf, die mir erst nach der Abgabe
dieser Arbeit aufgefallen sind:
\begin{itemize}
\item In Kapitel \ref{exp} ist die Rede von einem exponentiellen Abfall der Druck"=Drehmomentkurven. Aus der
Theorie von Maxwells Randbedingen geht aber hervor, dass es ein Abfall proportional zu $1/x$ sein sollte
(wurde leider nicht explizit im Theorieteil hergeleitet). Da weder von mir noch von den Autoren von
\cite{geyer} und \cite{duffner} versucht wurde durch mathematische Tests die stärke des Abfalls zu bestimmen, spricht
ad hoc nichts gegen die Annahme, dass er in der Tat proportional zu $1/x$ und nicht zu $\e^{-x}$ ist.
\end{itemize}

%\bibliographystyle{plain}
\bibliographystyle{IEEEtran.bst}
\bibliography{referenzen.bib}