Finite-Volumen-Verfahren

Das Finite-Volumen-Verfahren ist ein numerisches Verfahren zur Diskretisierung von Erhaltungsgleichungen, also von speziellen, partiellen Differentialgleichungen, denen ein Erhaltungssatz zugrunde liegt.

Bei korrekter Anwendung des Finite-Volumen-Verfahrens werden die Erhaltungseigenschaften dieser Gleichungen bewahrt, man spricht deswegen auch von einem konservativen Diskretisierungs-Verfahren. Vor allem aus diesem Grund hat sich das Finite-Volumen-Verfahren in der numerischen Strömungsmechanik durchgesetzt. Es wird aber auch in der numerischen Strukturmechanik und Elektrotechnik angewendet.

Das Berechnungsgebiet wird bei diesem Verfahren durch finite Volumen diskretisiert, die eine beliebige polygonale oder polyedrische Gestalt aufweisen können. Deswegen können auch komplizierte Geometrien einfach vernetzt werden. Approximiert werden die Funktionswerte der gesuchten Größen in den Mittelpunkten dieser finiten Volumen.

Herleitung

Ein Erhaltungssatz ist durch eine Gleichung der Art

${\frac {\partial }{\partial t}}u(x,t)+\nabla \cdot f(u(x,t))=g(u(x,t))$

auf einem Gebiet $\Omega \subset {\mathbb {R}}^{d}\times [0,\infty )$ gegeben, wobei $\nabla$ den Nabla-Operator bezeichnet, welcher hier für die Divergenz steht. Der gewöhnliche Fall, der hier betrachtet wird, ist g(u)=0 . Die Herleitung für Gleichungen mit weiteren Termen erfolgt analog. Zunächst wird das Gebiet in eine endliche (finite) Zahl an Volumenelementen (vgl. Gitterzellen) zerlegt. In jeder dieser Zellen gilt der Erhaltungssatz. Erfüllt jede der Zellen die Bedingungen des gaußschen Integralsatzes, etwa Lipschitz-Stetigkeit des Randes, so ergibt die Integration über eine Zelle $\Omega _{i}$ und Umwandlung des Integrals der Divergenz in ein Oberflächenintegral:

$\int _{{\Omega _{i}}}{\frac {\partial }{\partial t}}u\,{\mathrm d}\Omega +\int _{{\partial \Omega _{i}}}f(u)\cdot n\,{\mathrm d}S=0$ .

Die Veränderung einer erhaltenen Größe (z. B. der Energie) in einer Zelle kann also nur durch Ab- oder Hinzufließen (in diesem Fall von Energie) über den Rand der Zelle passieren. In jeder Zelle berechnet man nun den Mittelwert der Erhaltungsgröße in dieser Zelle $u_{i}={\frac {1}{|\Omega _{i}|}}\int _{{\Omega _{i}}}u\,{\mathrm d}\Omega$ und erhält im Falle, dass sich die Zelle mit der Zeit nicht verändert, eine Gleichung, welche die Veränderung der Mittelwerte in den Zellen mit der Zeit beschreibt:

${\frac {\partial }{\partial t}}u_{i}=-{\frac {1}{{|\Omega _{i}|}}}\int _{{\partial \Omega _{i}}}f(u)\cdot n\,{\mathrm d}S$ .

In numerischen Verfahren wählt man üblicherweise polygonal berandete Zellen, so dass sich das Integral über den Rand als Summe von Oberflächenintegralen über einfache Gebilde (im zweidimensionalen Fall gerade Kanten), darstellen lässt.

Lösung der Gleichung

Zur Berechnung der Oberflächenintegrale wird im Regelfall eine Gauß-Quadratur zweiter Ordnung genommen. Nach Mittelung der Werte in den einzelnen Zellen ergibt sich das Problem, dass die numerische Lösung entlang der Gitterkanten unstetig ist. Allerdings lässt sich die Situation an der Kante als Riemann-Problem auffassen. Die Verwendung eines approximativen Riemann-Lösers erlaubt dann die Berechnung der Flüsse. Hierbei wird Konsistenz des Riemann-Lösers verlangt, was in diesem Fall zum einen die Stetigkeit oder sogar Lipschitz-Stetigkeit bedeutet, sowie die Bedingung, dass er bei identischen Daten aus beiden Zellen den physikalischen Fluss liefert.

Diese liefert das dann zu erstellende System von gewöhnlichen Differentialgleichungen nur, wenn noch eine Entropie-Bedingung hinzugenommen wird. Denn die rein mathematische Betrachtung der Unstetigkeit am Zellenrand erlaubt neben der für das Riemann-Problem korrekten Lösung vermöge eines Verdichtungsstosses auch den unphysikalischen Verdünnungsstoss. Die Entropiebedingung aber schließt den Verdünnungsstoss aus. Das Riemann-Problem wird dann mittels numerischer Methoden für gewöhnliche Differentialgleichungen (unter Beachtung der Entropiebedingung) näherungsweise (z.B. Osher) oder iterativ-exakt (Godunov) gelöst.

1D Beispiel

Es wird die folgende Konvektionsgleichung betrachtet:

${\frac {\partial u}{\partial t}}=-v{\frac {\partial u}{\partial x}}\,,\qquad t\geq 0\,,\qquad x\in \left[x_{\mathrm {min} },x_{\mathrm {max} }\right]$

Dabei beschreibt v die Konvektionsgeschwindigkeit. Diese partielle Differentialgleichung soll mit Hilfe des Finite-Volumen-Verfahrens entlang der Ortskoordinate diskretisiert werden. Dazu wird zunächst die Ortskoordinate in n diskrete Abschnitte zerlegt

$x_{\mathrm {min} }=x_{\frac {1}{2}}<x_{\frac {3}{2}}<\dots <x_{n-{\frac {1}{2}}}<x_{n+{\frac {1}{2}}}=x_{\mathrm {max} }$

Die Kontrollvolumina $\Omega _{i}$ , deren Mittelpunkte $x_{i}$ und die Größen der Kontrollvolumina $\Delta x_{i}$ sind definiert durch

$\Omega _{i}=\left[x_{i-{\frac {1}{2}}},x_{i+{\frac {1}{2}}}\right]$

$x_{i}={\frac {x_{i-{\frac {1}{2}}}+x_{i+{\frac {1}{2}}}}{2}}$

$\Delta x_{i}=x_{i+{\frac {1}{2}}}-x_{i-{\frac {1}{2}}}$

Hier beziehen sich ganzzahlige Indizes auf den Mittelpunkt eines Kontrollvolumens und nicht-ganzzahlige Indizes beziehen sich auf den Rand eines Kontrollvolumens.

Nun wird jeder Term der zu diskretisierenden Differentialgleichung über ein beliebiges inneres Kontrollvolumen integriert. Für den Akkumulationsterm folgt

$\int _{x_{i-{\frac {1}{2}}}}^{x_{i+{\frac {1}{2}}}}{\frac {\partial u}{\partial t}}\,\mathrm {d} x={\frac {\mathrm {d} }{\mathrm {d} t}}\int _{x_{i-{\frac {1}{2}}}}^{x_{i+{\frac {1}{2}}}}u\,\mathrm {d} x={\frac {\mathrm {d} u_{i}}{\mathrm {d} t}}\Delta x_{i}$

Dabei darf nach der Leibnizregel für Parameterintegrale die Integration und Differentiation vertauscht werden, sofern das Kontrollvolumen zeitlich unveränderlich ist. Anschließend wird der diskretisierte Funktionswert $u_{i}$ mit Hilfe des Mittelwertsatzes der Integralrechnung als integraler Mittelwert des tatsächlichen Funktionsverlaufes im Kontrollvolumen approximiert.

Anschließend folgt für den konvektiven Term mit Hilfes des Gaußschen Integralsatzes

$-v\int _{x_{i-{\frac {1}{2}}}}^{x_{i+{\frac {1}{2}}}}{\frac {\partial u}{\partial x}}\,\mathrm {d} x=-v\left(\left.u\right|_{x_{i+{\frac {1}{2}}}}-\left.u\right|_{x_{i-{\frac {1}{2}}}}\right)$

Hier müssen die Funktionswerte auf den Rändern der Kontrollvolumina $\left.u\right|_{x_{i+{\frac {1}{2}}}}$ als Funktion der bekannten Funktionswerte $u_{i}$ approximiert werden. Eine einfache Methode dazu ist das sog. UPWIND-Verfahren 1. Ordnung, das besagt, dass der Funktionswert auf dem Rand durch den nächsten bekannten stromaufwärtsliegenden Funktionswert approximiert wird, d.h.

$\left.u\right|_{x_{i+{\frac {1}{2}}}}={\begin{cases}u_{i}&:v>0\\u_{i+1}&:v<0\end{cases}}$

Werden die beiden Terme wieder zusammengesetzt, erhält man den Satz gewöhnlicher Differentialgleichungen

${\frac {\mathrm {d} u_{i}}{\mathrm {d} t}}=-{\frac {v}{\Delta x_{i}}}\left(u_{i}-u_{i-1}\right)$

der mithilfe eines beliebigen Verfahrens für gewöhnliche Differentialgleichungen gelöst werden kann.

Verfahren Höherer Ordnung

Das bisher beschriebene Verfahren ist durch die Mittelung der Werte in jeder Zelle nur erster Ordnung. Höhere Ordnung wird dadurch erreicht, dass Polynome höherer Ordnung in den Zellen angesetzt werden. D.h. es wird eine Verteilung (konstant, linear, parabolisch usw.) angenommen, die den Integralwert erhält.

Die zentrale Schwierigkeit hierbei ist, dass Verdichtungsstösse bzw. Schocks in der Lösung zu Oszillationen führen können. Zur Vermeidung dessen werden Total-Variation-Diminishing-Verfahren (TVD-Verfahren) eingesetzt, die die totale Variation nicht erhöhen und so keine neuen Extrema zulassen (Da Polynome Unstetigkeiten i.a. mit Überschwingern interpolieren). Die wichtigsten Klassen von Verfahren sind hier die Flux-Limiter-Verfahren und die ENO-Verfahren (bzw. WENO).

Konvergenztheorie

Finite-Volumen-Verfahren lassen sich für elliptische Gleichungen als spezielle Finite-Elemente-Verfahren auffassen, bei denen man stückweise konstante bzw. stückweise lineare Ansatzfunktionen wählt, die auf den Zellen und nicht auf den Gitterpunkten leben. Dies erlaubt mit Hilfe der dortigen umfassenden Theorie eine Konvergenzanalyse.

Für parabolische oder hyperbolische Gleichungen wie die Euler- oder Navier-Stokes-Gleichungen ist die mathematische Konvergenztheorie allerdings weniger weit fortgeschritten.

Hyperbolische Erhaltungsgleichungen

Bei hyperbolischen Problemen treten insbesondere Stöße auf, die die Analyse erheblich erschweren. Eine weitere Schwierigkeit hier besteht darin, dass die Lösung der Gleichungen im Regelfall nicht eindeutig ist. Der Satz von Lax-Wendroff liefert, dass ein Finite-Volumen-Verfahren bei Konvergenz tatsächlich gegen eine schwache Lösung der Gleichung konvergiert. Entropiebedingungen bzw. numerische Viskosität werden dann genutzt, um zu zeigen, dass dies tatsächlich die physikalisch sinnvolle ist.

Die Konvergenz einer numerischen Approximation $u_{\Delta t}$ zu einer exakten Lösung ist mittels des globalen Fehlers $E_{\Delta t}(x,t):=u_{\Delta t}(x,t)-u(x,t)$ definiert:

$\|E_{\Delta t}(\cdot ,t)\|\longrightarrow 0\quad {\text{für}}\;\Delta t\rightarrow 0,\quad t\in \mathbb {R} ^{+}.$

Eine andere Aussage, die für alle Finite-Volumen-Verfahren gilt, ist die notwendige CFL-Bedingung, dass der numerische Abhängigkeitsbereich den tatsächlichen Abhängigkeitsbereich enthalten muss. Andernfalls ist das Verfahren instabil.

Insbesondere für mehrdimensionale Gleichungen ist die Konvergenztheorie schwierig. Im Eindimensionalen gibt es auch für Verfahren höherer Ordnung Resultate, die darauf beruhen, dass der Raum der Funktionen mit beschränkter Variation kompakte Mengen in L₁ liefert.

Konsistenz

Ein numerisches Verfahren $H_{\Delta t}$ ist eine Vorschrift zur Konstruktion der numerischen Approximation für den nächsten Zeitschritt:

$u_{j}^{n+1}=H_{\Delta t}(U^{n},j),$

wobei der Zeitindex, der Ortsindex und $U^{n}$ die approximative Lösung zum vorigen Zeitpunkt aller Ortsindizes darstellt:

$t=t_{n},\quad x=x_{j},\quad U^{n}=\left(u_{1}^{n},u_{2}^{n},\cdots \right).$

Das numerische Verfahren $H_{\Delta t}$ kann auch kontinuierlich definiert sein:

$u(x,t+\Delta t)=H_{\Delta t}\left(u(\cdot ,t),x\right).$

Der lokalen Abschneidefehler $L_{\Delta t}$ des Verfahrens $H_{\Delta t}$ ist definiert als:

$L_{\Delta t}(x,t):={\frac {1}{\Delta t}}\left(u(x,t+\Delta t)-H_{\Delta t}\left(u(\cdot ,t),x\right)\right),$

wobei als exakte Lösung zur Zeit angenommen wird. Zur Definition der Konsistenz (Numerik) des Verfahrens $H_{\Delta t}$ nutzt man nun den lokalen Abschneidefehler $L_{\Delta t}$ :

$\|L_{\Delta t}(\cdot ,t)\|\longrightarrow 0\quad {\text{mit}}\;\Delta t\rightarrow 0,\quad t\in \mathbb {R} ^{+}.$

Das Verfahren $H_{\Delta t}$ hat Konsistenzordnung $n\in \mathbb {N}$ , falls es ein C>0 gibt, sodass:

$\|L_{\Delta t}(\cdot ,t)\|\leq C(\Delta t)^{n}\quad {\text{mit}}\;\Delta t\rightarrow 0,\quad t\in \mathbb {R} ^{+}$

gibt.

Geschichte

Die grundlegenden theoretischen und praktischen Ideen wurden ab den 1950er Jahren für die Raumfahrt entwickelt, insbesondere von dem Russen Godunow und dem Ungarn Lax. Die ersten Finite-Volumen-Verfahren stammen von Richard S. Varga (1962) und Preissmann (1961). Der Terminus Finite-Volumen-Verfahren wurden dann unabhängig voneinander 1971 von McDonald und 1972 von MacCormack und Paullay für die Lösung der zeitabhängigen zweidimensionalen Euler-Gleichungen eingeführt.

Die Idee der approximativen Riemann-Löser tauchte erstmals in den 1980ern auf, als Roe, Osher, van Leer und andere ebenfalls unabhängig voneinander solche Verfahren vorstellten.

Basierend auf einem Artikel in:

Wikipedia.de