10  Transportprozesse

Lässt man in zellulären Automaten aus Kapitel 9 auch kontinuierliche Werte zu, und betrachtet die einzelnen Zellen als Diskretisierung eines kontinuierlichen Bereichs (z.B. eines Sees), so kann man mit der “Automatenstruktur” auch Strömungsprozesse beschreiben. Von besonderer Bedeutung sind hierbei die Advektion (von lat. advehi = heranbewegen; bezeichnet allgemein das Heranführen von Dingen) und die Diffusion. Mit Advektion bezeichnet man z.B. den Transport von einer in Wasser gelösten Substanz bzw. das durch eine Strömung transportierte Wasser.

Beispiel: Man nimmt einen durchsichtigen Gartenschlauch und injiziert an einer Stelle etwas Tinte in den Schlauch. Macht man dieses Experiment in der Realität, wird man feststellen, das die Bande aus Tinte mit der Zeit breiter und heller wird. Dies liegt daran, dass die Tintenteilchen im realen System das Bestreben haben, die Konzentration im Schlauch auszugleichen. Dieses Bestreben nennt man Diffusion. Dreht man nun das Wasser auf, so strömt die Tinte entsprechend der angelegten Geschwindigkeit durch den Schlauch. Dieses bezeichnet man als Advektion.

10.1 Advektion

Die Beschreibung der Advektion (Strömung), die durch die Strömungsgeschwindigkeit \(v\) bestimmt ist, lässt sich im eindimensionalen Fall wie folgt beschreiben: Angenommen, wir betrachten ein kleines Volumenelement \(V= \Delta x \cdot \Delta y \cdot \Delta z\) des Schlauchs (Abbildung 10.1), in das mit der Geschwindigkeit \(v\) (in m/s) Wasser strömt. Das Wasser sei inkompressibel und die Konzentration bzw. Dichte der Tinte sei durch die noch unbekannte Funktion \(c(t,x)\) (in mg/m\(^3\)) beschrieben. Vereinfachend nehmen wir an, dass in y- und z-Richtung nichts passiert. Die Strömung verlaufe nur in x-Richtung.

Abbildung 10.1: Das Volumenelement \(V= \Delta x \cdot \Delta y \cdot \Delta z\), in das in x-Richtung mit der Geschwindigkeit \(v\) Wasser strömt..

Die Masse an Tinte, die nun in einem Zeitschritt \(\Delta t\) durch die Strömung mit der Strömungsgeschwindigkeit \(v\) in das Volumen \(\Delta x \cdot \Delta y \cdot \Delta z\) gelangt, ist durch

\[ m_{in} = c(t,x-\Delta x) \cdot v(t,x-\Delta x) \cdot \Delta t \cdot \Delta y \cdot \Delta z \] wobei \(t\) die aktuelle Zeit und \(x-\Delta x\) der linke Rand des Volumenelements ist. Ist die Flüssigkeit inkompressibel, fließt dasselbe Volumen raus. Die Masse, die herau sfließt, ist nun durch \[ m_{out} = c(t,x) \cdot v(t,x) \cdot \Delta t \cdot \Delta y \cdot \Delta z \] gegeben Die Massendifferenz beträgt also \[ \Delta m = \left [ c(t,x) \cdot v(t,x) -c(t,x-\Delta x) \cdot v(t,x-\Delta x) \right ] \cdot \Delta t \cdot \Delta y \cdot \Delta z \]

Sie ist positiv, wenn Masse aus dem Volumen heraus fließt. Die zeitliche Änderung der Masse im Volumen \(V= \Delta x \cdot \Delta y \cdot \Delta z\) innerhalb eines Zeitschritts entspricht also gerade dem negativen von diesem Wert, da die Masse in V abnimmt, wenn etwas heraus fließt:

\[ m(t+\Delta t,x)-m(t,x) = - \Delta m \]

Durch Erweiterung der rechten Seite mit \(\frac{\Delta x}{\Delta x}\) erhält man schließlich

\[ m(t+\Delta t,x)-m(t,x) =-\frac{c(t,x)\cdot v(t,x) -c(t,x-\Delta x)\cdot v(t,x-\Delta x)}{\Delta x} \cdot V \cdot \Delta t \]

Division durch \(V \cdot \Delta t\) ergibt für konstantes \(V\) und mit \(c=\frac{m}{V}\)

\[ \frac{c(t+\Delta t,x)-c(t,x)}{\Delta t} = - \frac{c(t,x)\cdot v(t,x) -c(t,x-\Delta x)\cdot v(t,x-\Delta x)}{\Delta x} \tag{10.1}\]

Die Grenzübergänge \(\Delta t \rightarrow 0\) und \(\Delta x \rightarrow 0\) führen nun zu

\[ \frac{ \partial c(t,x)}{\partial t} = - \frac{ \partial (c(t,x) \cdot v(t,x))}{\partial x} \]

Mit der Produktregel erhält man

\[ \frac{ \partial c(t,x)}{\partial t} = - \left [ c(t,x)\cdot \frac{ \partial v(t,x)}{\partial x} +v(t,x) \cdot \frac{ \partial c(t,x)}{\partial x}\right ] \]

Für konstante Geschwindigkeit \(v\) gilt daher

\[ \frac{ \partial c(t,x)}{\partial t} = - v \cdot \frac{ \partial c(t,x)}{\partial x} \tag{10.2}\]

Gleichung 10.2 heißt Advektionsgleichung (für konstante Geschwindigkeiten). Es handelt sich um eine partielle Differentialgleichung, da die Ableitungen nach \(t\) und nach \(x\) auftauchen.

10.1.1 Diskretisierung der Advektionsgleichung

Um nun Gleichung 10.2 numerisch zu lösen, muss die Gleichung diskretisiert werden. Der Schlauch wird durch ein eindimensionales Gitter simuliert. Hierzu wird die die x-Achse in Abschnitte der Länge \(\Delta x\) unterteilt. Jeder Abschnitt beschreibt ein kurzes Stückchen des Schlauches (Gitterzelle), in dem die Substanz als homogen verteilt angenommen wird. Es wird nun die zeitliche Änderung der Konzentration in der \(i\)-ten Gitterzelle bestimmt. Hierzu werden zuerst die Gradienten zu benachbarten Gitterzellen diskretisiert, indem man Differenzenquotienten betrachtet:

\[ \frac{c(t+\Delta t ,x)-c(t ,x)}{\Delta t} + v \cdot \frac{c(t,x) - c(t,x-\Delta x) }{\Delta x} =0 \tag{10.3}\]

Dies entspricht wieder Gleichung 10.1, was aber nicht zwingend ist. Es gibt andere Diskretisierungen. In Gleichung 10.3 ist die räumliche Ableitung nicht symmetrisch sondern stromaufwärts verschoben. Man nennt dies daher auch Upstream-Schema. Solch ein Schema macht Sinn, wenn \(v\) nicht die Richtung wechselt (\(v>0\)). Es ist höchstens dann stabil, wenn

\[ v \cdot \frac{\Delta t}{\Delta x} \le 1 \qquad \text{\textbf{Courant-Friedrichs-Lewy-Kriterium (CFL)}} \] gilt. Die Größe \(\vert v \cdot \frac{\Delta t}{\Delta x} \vert\) heißt Courant-Zahl. Sie beschreibt, welche Wegstrecke (\(v \cdot \Delta t\)) in einem Zeitintervall \(\Delta t\) zurückgelegt wird. Diese Strecke darf die räumliche Diskretisierung \(\Delta x\) nicht überschreiten.

Abbildung 10.2: Diskretisierung des Schlauchs in Gitterzellen der Länge \(\Delta x\), in die in x-Richtung mit der Geschwindigkeit \(v\) Wasser strömt.

Betrachtet man nun Gitterzellen mit Abstand \(\Delta x\), erhält man folgendes Rechenschema. \[c(t+\Delta t ,i)= c(t ,i) + v \cdot \frac{c(t,i-1) - c(t,i) }{\Delta x} \cdot \Delta t \tag{10.4}\]

Diese Art der Diskretisierung nennt man auch Methode der finiten Differenzen. Das Problem bei der Diskretisierung der Advektion ist, dass \(\Delta x\) und \(\Delta t\) genau zueinander “passen” müssen. Dies kann man aber höchstens zufällig bei konstanter Geschwindigkeit \(v\) erreichen. Ansonsten “verschmiert” das Signal, die Bande aus Tinte wird immer breiter. Dieses sieht ähnlich wie eine Diffusion aus. Daher wird dieses Phänomen, das sich aufgrund der Diskretisierung ergibt, auch numerische Diffusion genannt. Man kann zeigen, dass es sich bei dem Fehler, den man durch die Diskretisierung macht, tatsächlich um einen Diffusionsterm handelt.

10.2 Diffusion

Die Diffusion kann man wie folgt formalisieren: Angenommen wir betrachten den Schlauch als eindimensional (x-Richtung) und geben die Tinte als Konzentration \(c\) an. Dann ist \(c(t,x)\) die Konzentration zur Zeit \(t\) am Ort \(x\).

Nach dem 1. Fickschen Gesetz ist die Teilchenstromdichte (der Fluss) \(J\) proportional zum Konzentrationsgradienten entgegen der Diffusionsrichtung:

\[J = -D \frac{\partial\ c}{\partial x}\]

Hierbei ist \(D\) die Proportionalitätskonstante, die Diffusionskonstante. Das Minuszeichen bewirkt, dass der Fluss entgegen dem Konzentrationsgradienten wirkt, also von höherer zu niedrigerer Konzentration. Ist die Konzentration c z.B. in \(\frac{mol}{m^3}\) angegeben und die Zeit in s, so ist die Einheit des Flusses J \(\frac{mol}{m^{2}s}\). Die Einheit der Diffusionskonstanten ist \(\frac{m^2}{s}\).

Die zeitliche Änderung der Konzentration ergibt sich aus der Kontinuitätsgleichung, die besagt, dass keine Masse verloren geht oder entsteht, bzw. aus der Advektionsgleichung 10.2, bei der der Fluss \(c\cdot v\) nun durch \(J\) ersetzt wird:

\[\frac{\partial c}{\partial t} +\frac{\partial J}{\partial x} =0\]

\[\frac{\partial c}{\partial t} = -\frac{\partial J}{\partial x} = - \frac{\partial }{\partial x} \left [ -D \cdot \frac{\partial\; c}{\partial x} \right ]\]

Für eine Diffusionskonstante unabhängig von x ergibt sich

\[\frac{\partial c}{\partial t} =D \cdot \frac{\partial^2 c}{\partial x ^2}\]

Um die Diffusionsgleichung für ein System zu lösen, muss man eine Anfangsbedingung vorgeben. In unserem Fall die Konzentration an jeder Stelle \(x\) zur Zeit \(t=0\) , also \(c(0,x)\) für alle \(x\).

Bisher haben wir noch keine Angabe über die Länge des Schlauches gemacht. Wir sind von einem unendlich langen Schlauch ausgegangen. Will man die Verteilung in einem endlichen Schlauchstück beschreiben, so muss man auch Randbedingungen definieren. Hier gibt es zwei Möglichkeiten: Entweder man setzt die Konzentration am Anfang (x=0) und am Ende des Schlauchstücks (x=L) fest, also

\[c(t,0) \text{ und } c(t,L) \text{ f.a } t \ge 0 \qquad \textbf{Dirichlet-Problem}\]

oder man gibt die Ableitungen am Rand vor

\[\frac{\partial c}{\partial t}c(t,0) \text{ und }\frac{\partial c}{\partial t}c(t,L) \text{ f.a } t \ge 0 \qquad \textbf{Neumann-Problem}\]

Der Unterschied ist wie folgt: Sind die Konzentrationen am Rand vorgegeben, so bestimmen diese letztlich nach einer gewissen Zeit die Konzentration im Schlauchstück. Vom Rand her diffundiert Substanz (Tinte) in den Schlauch hinein. Hält man den Gradienten am Rand fest, z.B. setzt in gleich null, dann bedeutet dies, dass keine Substanz aus dem Schlauch oder in den Schlauch diffundieren kann. Die Menge im Schlauch bleibt konstant und verteilt sich mit der Zeit gleichmäßig.

Im allgemeinen ist es sehr schwierig, Anfangs- und Randwertprobleme analytisch zu lösen. Daher ist es meistens sinnvoll solche Probleme numerisch zu behandeln. Hierzu gibt es hervorragende Verfahren. Das Grundprinzip soll im folgenden mit einem sehr elementaren Verfahren erläutert werden:

Diskretisierung der Diffusionsgleichung

Der Schlauch wird wieder durch ein eindimensionales Gitter simuliert. Hierzu wird die die x-Achse in Abschnitte der Länge \(\Delta x\) unterteilt. Jeder Abschnitt beschreibt ein kurzes Stückchen des Schlauches (Gitterzelle), n dem die Substanz als homogen verteilt angenommen wird. Es wird nun die zeitliche Änderung der Konzentration in der \(i\)-ten Gitterzelle bestimmt. Hierzu werden zuerst die Gradienten zu den benachbarten Gitterzellen diskretisiert, indem man die jeweiligen Differenzenquotienten betrachtet:

\[\begin{aligned} \frac{\partial c(t,x)}{\partial t} &\approx& D \cdot \frac{\frac{\partial c(t,x+\Delta x)}{\partial x} - \frac{\partial c(t,x)}{\partial x} }{\Delta x} \\ &&\\ &\approx& D \cdot \frac{\frac{ c(t,x+\Delta x) - c(t,x)}{\Delta x} - \frac{c(t,x) - c(x-\Delta x)}{\Delta x} }{\Delta x}\\ &&\\ &=& D \cdot \frac{c(t,x+\Delta x) - 2 \cdot c(t,x) + c(x-\Delta x)}{\Delta x ^2}\\ \end{aligned}\]

Die zeitliche Diskretisierung führt zu

\[ c(t+\Delta t ,x)= c(t ,x) + D \cdot \frac{c(t,x+\Delta x) - 2 \cdot c(t,x) + c(t,x-\Delta x)}{\Delta x ^2} \cdot \Delta t \tag{10.5}\]

Betrachtet man nun Gitterzellen mit Abstand \(\Delta x\) so erhält man folgende Formel:

\[c(t+\Delta t ,i)= c(t ,i) + D \cdot \frac{c(t,i+1) - 2 \cdot c(t,i) + c(i-1)}{\Delta x ^2} \cdot \Delta t\]

Diese Art der Diskretisierung nennt man forward in time und centered in space (FTCS). Man muss nun \(\Delta x\) und \(\Delta t\) so wählen, dass tatsächlich ein Konzentrationsausgleich geschaffen wird, und dass das diskretisierte System nicht oszilliert. Dies ist für

\[D< \frac{1}{2}\cdot \frac{\Delta x^2}{\Delta t}\]

der Fall. Dies versteht man leicht, wenn man berücksichtigt, dass pro Zeitschritt aus einer Zelle nicht mehr (Masse) herausfließen darf als drin ist. Dies bedeutet, dass

\[\frac{2 D \cdot c(t,x) }{\Delta x ^2} \cdot \Delta t < c(t ,x)\]

gelten muss, da in Gleichung 10.5 die Konzentration in den Nachbarzellen schlimmstenfalls null ist.

10.2.1 Zweidimensionale Advektion und Diffusion 

Eine sehr einfache Möglichkeit, die Transportprozesse im zweidimensionalen Fall zu beschreiben, funktioniert ähnlich wie bei den Automaten. Man betrachtet jetzt nicht nur die Nachbarzellen rechts und links, sondern auch oben und unten. Man definiert also für jede Zelle vier Austauschflüsse. Nehmen wir vereinfachend an, dass die räumliche Diskretisierung in beide Koordinatenrichtung gleich \(\Delta x=\Delta y\), so gilt für die Verfahren erster Ordnung für die Advektion mit den Geschwindigkeiten \(v_x >0\) und \(v_y>0\) in die Koordinatenrichtungen

\[\begin{aligned} c(t+\Delta t ,i,j)&= c(t ,i,j) \\[1ex]\nonumber &+ v_x \cdot \frac{c(t,i-1,j) - c(t,i,j) }{\Delta x} \cdot \Delta t \\[1ex]\nonumber &+v_y \cdot \frac{c(t,i,j-1) - c(t,i,j) }{\Delta x} \cdot \Delta t\\\nonumber \end{aligned}\]

Das CFL-Kriterium muss nun für beide Raumrichtungen und die Raumdiagonale erfüllt sein. Man kann zeigen, dass dies für

\[ (v_x+v_y) \cdot\frac{\Delta t}{\Delta x} \le 1 \]

gilt. Man sollte aber nicht “=1” ausreizen, da die Numerik der verwendeten Programmiersprache (Zahldarstellung) zu falschen Ergebnissen führen kann.

Für die Diffusion ergibt sich das Schema

\[\begin{aligned} c(t+\Delta t ,i,j) &= c(t ,i,j) \\[1ex]\nonumber &+ D \cdot \frac{c(t,i+1,j) - 2 \cdot c(t,i,j) + c(t, i-1,j) }{\Delta x ^2} \cdot \Delta t\\[1ex]\nonumber &+ D \cdot \frac{c(t,i,j+1) - 2 \cdot c(t,i,j) + c(t,i,j-1)}{\Delta x ^2} \cdot \Delta t\\\nonumber \end{aligned}\]

mit der Stabilitätsbedingung

\[D< \frac{1}{4}\cdot \frac{\Delta x^2}{\Delta t}\]

10.2.2 Numerische Diffusion 

Die angegebene Diskretisierung der Advektion (Gleichung 10.1) lässt sich auch herleiten, indem man die Konzentration als Taylorpolynom (siehe Abschnitt A.12) in \(t\) und \(x\) darstellt. Hierbei muss wiederum die Asymmetrie der räumlichen Ableitung berücksichtigt werden.

\[\begin{aligned} c(t+\Delta t, x) &= c(t,x)+ \frac{\partial c(t,x)}{\partial t}\cdot \Delta t + \frac{1}{2} \frac{\partial^2 c(t,x)}{\partial t^2} \cdot \Delta t^2 + \mathcal{O}(\Delta t^3) \\[2ex] c(t, x-\Delta x) &= c(t,x)- \frac{\partial c(t,x)}{\partial x}\cdot \Delta x + \frac{1}{2} \frac{\partial^2 c(t,x)}{\partial x^2} \cdot \Delta x^2 + \mathcal{O}(\Delta x^3) \end{aligned}\]

Bei der hier dargestellten Entwicklung bis zum quadratischen Term sind die Fehler \(\mathcal{O}\) von dritter Ordnung und werden in der nachfolgenden Rechnung vernachlässigt. Löst man die Gleichungen nach \(\frac{\partial c(t,x)}{\partial t}\) bzw. \(\frac{\partial c(t,x)}{\partial x}\) auf und setzt das Ergebnis in Gleichung 10.2 ein, so erhält man

\[\begin{aligned} \nonumber &\frac{c(t+\Delta t, x) - c(t,x)}{\Delta t} - \frac{1}{2} \frac{\partial^2 c(t,x)}{\partial t^2} \cdot \Delta t \\[2ex] + v \cdot &\left ( \frac{c(t, x) - c(t,x- \Delta x)}{\Delta x} + \frac{1}{2} \frac{\partial^2 c(t,x)}{\partial x^2} \cdot \Delta x \right ) =0 \end{aligned} \tag{10.6}\]

Die zweite zeitliche Ableitung wird man wieder los, indem man Gleichung 10.2 berücksichtigt und annimmt, dass \(v\) konstant ist:

\[\begin{aligned} \nonumber \frac{\partial c(t,x)}{\partial t} &=- v \cdot \frac{\partial c(t,x)}{\partial x}\\[2ex]\nonumber &\text{also}\\[2ex]\nonumber \frac{\partial^2 c(t,x)}{\partial t^2}&=-v \cdot \frac{\partial}{\partial t} \left [\frac{\partial c(t,x)}{\partial x} \right ] = -v \cdot \frac{\partial}{\partial x} \left [\frac{\partial c(t,x)}{\partial t} \right ] = -v \cdot \frac{\partial}{\partial x} \left [-v \cdot \frac{\partial c(t,x)}{\partial x} \right ]\\[2ex]\nonumber & = v^2 \cdot \frac{\partial c(t,x)}{\partial x^2}\\ \end{aligned}\]

Setzt man dies nun in Gleichung 10.6 ein und fasst zusammen, so erhält man

\[\begin{aligned} &\frac{c(t+\Delta t, x) - c(t,x)}{\Delta t} + v \cdot \frac {c(t, x) - c(t,x- \Delta x)}{\Delta x} \\[2ex]\nonumber = &\frac{1}{2} \cdot v^2 \cdot \frac{\partial^2 c(t,x)}{\partial x^2} \cdot \Delta t - \frac{1}{2} \cdot v \cdot \frac{\partial^2 c(t,x)}{\partial x^2} \Delta x \\[2ex]\nonumber = &-\frac{1}{2} \cdot v \cdot \left ( \Delta x - v \cdot \Delta t\right )\cdot \frac{\partial^2 c(t,x)}{\partial x^2}\\[2ex]\nonumber \end{aligned} \tag{10.7}\]

Dieses Schema ergibt nun bis auf die höheren Terme die richtige Lösung. Vergleicht man es mit dem Upstream-Schema fällt der Term auf der rechten Seite auf. Beim Upstream-Schema wird also zusätzlich das Negative dieses Terms berechnet:

\[\begin{aligned} \frac{1}{2} \cdot v \cdot \left ( \Delta x - v \cdot \Delta t\right )\cdot \frac{\partial^2 c(t,x)}{\partial x^2}\\[2ex]\nonumber \end{aligned} \tag{10.8}\]

Hierbei handelt es sich um einen Diffusionsterm mit der Diffusionskonstanten \(D= \frac{1}{2} \cdot v \cdot \left (\Delta x- v \cdot \Delta t \right ) > 0\). Dies ist die sogenannte numerische Diffusion. Aus der Tatsache, dass eine Diffusionskonstante positiv sein muss, folgt sofort dass das Schema nur für \(v>0\) angewendet werden darf. Dreht man die Stromrichtung um, muss man das numerische Schema anpassen. Für negatives \(D\) fliegt uns das ganze sonst um die Ohren, da keine Massenerhaltung vorliegt. Weiterhin folgt für \(v>0\) sofort, dass \(\Delta x -v \cdot \Delta t> 0\) gelten muss und man erhält die Stabilitätsbedingung 10.4.

Durch die Vernachlässigung der Terme zweiter Ordnung beim Upstream-Schema entsteht ein Fehler, der sich wie eine Diffusion auswirkt. Man kann das Verfahren verbessern, indem man den Term auf der rechten Seite von Gleichung 10.7 berücksichtigt, indem man ihn wie in Abschnitt 10.2 diskretisiert. Das Verfahren, das man so erhält, heißt Lax-Wendroff-Verfahren:

\[\begin{aligned} c(t+\Delta t, x)&= c(t,x) +a\cdot \left ( c(t,x- \Delta x)-c(t,x)\right )\\[2ex]\nonumber &-\frac{a-a^2}{2}\cdot \left ( c(t,x+ \Delta x) + c(t,x- \Delta x)-2\cdot c(t,x) \right ) \end{aligned} \tag{10.9}\]

mit der Courantzahl \(a=\frac{v\cdot\Delta t}{\Delta x}\).