Theorie und Numerik von Erhaltungsgleichungen

Dieses Buch steht im Regal Mathematik.

Erhaltungssätze

Bearbeiten

Erhaltungssätze sind physikalische Gesetze, die aussagen, dass eine bestimmte Größe, die Erhaltungsgröße, in einem abgeschlossenen System konstant ist, also erhalten bleibt. Sie spielen eine zentrale Rolle in vielen Bereichen der Physik. Neben starker experimenteller Absicherung lassen sie sich nach dem Noether-Theorem mit bestimmten grundlegenden Symmetrien in Verbindung bringen und sind somit auch theoretisch stark fundiert. Etwa gehört zur Energieerhaltung die Translationsinvarianz physikalischer Aussagen in der Zeit. Beispiele sind

  1. Energieerhaltung
  2. Impulserhaltung
  3. Drehimpulserhaltung
  4. Erhaltung von elektrischer Ladung

Kein echter Erhaltungssatz ist die Massenerhaltung. Diese gilt nur bei nichtrelativistischen Geschwindigkeiten; bei hohen Geschwindigkeiten muss die spezielle Relativitätstheorie in Betracht gezogen werden. Die berühmte Formel

 

wird in der Physik als Teil der Energieerhaltung betrachtet: Masse ist dann eine spezielle Form der Erhaltungsgröße „totale Energie“. Ebenso ergeben sich bei zerfallenden Atomkernen nach dem Zerfall kleinere Massensummen der Tochterteilchen.

Keine Erhaltungsgrößen sind beispielsweise Druck, Temperatur oder elektrische Spannung.

Erhaltungsgleichungen

Bearbeiten

Die mathematische Formulierung eines Erhaltungssatzes wird als Erhaltungsgleichung bezeichnet. Hier werden wir uns mit Strömungsmechanik beschäftigen. Eine der großen wissenschaftlichen Leistungen des 19. Jahrhunderts war die Erkenntnis, dass sich Strömungen durch Kopplung der Erhaltungsgleichungen für Masse, Energie und Impuls, zusammen mit einer Zustandsgleichung, komplett beschreiben lassen. Ausgangspunkt ist dabei die integrale Form. Wir betrachten als Beispiel die Massenerhaltung. Gegeben sei ein Volumen  . Die darin enthaltene Masse   ist dann gegeben durch das Integral der Dichte   über  :

 

Aussage des Erhaltungssatzes der Masse ist, dass Masse nicht zerstört oder erschaffen wird. Änderungen der Masse M im Volumen   mit der Zeit sind also Resultat von Massentransport über den Rand von  . Dies lässt sich formulieren als

 

wobei   die Massenflussdichte beschreibt. Diese lässt sich präzisieren als  , wobei   die Geschwindigkeit eines kleinen Massenvolumens beschreibt. Damit ergibt sich

 

Zu beachten ist hierbei, dass dabei außer der Integrierbarkeit von   und   keine Voraussetzungen an die Funktionen   und   oder an das Volumen   gestellt wurden. Eine zweite wichtige Form einer Erhaltungsgleichung ist die differentielle Form. Dazu wird die letzte Gleichung mittels des Integralsatzes von Gauß umgeformt:

 

Ist das Volumen   nicht zeitabhängig, lässt sich Zeitdifferentiation und Integration vertauschen. Da die Aussage für beliebige Volumina gilt, ergibt sich die differentielle Form der Massenerhaltung

 

Allgemeiner haben Erhaltungsgleichungen die Form

 

oder, wenn zusätzlich Volumenkräfte modelliert werden müssen,

 

Letztere Form wird im Englischen auch als Balance law bezeichnet (deutsch: Bilanzgleichung).

Eine dritte wichtige Form der Erhaltungsgleichung ist die schwache Form. Dazu multiplizieren wir die differentielle Form (hier für den eindimensionalen Fall) mit glatten Testfunktionen  , also aus dem Raum der stetig differenzierbaren Funktionen mit kompaktem Träger. Integration liefert

 

Partielle Integration ergibt aufgrund des kompakten Trägers von  

 

Eine Funktion   wird dann schwache Lösung der Erhaltungsgleichung genannt, wenn die obige Gleichung für alle   erfüllt ist.

Beispiele für Erhaltungsgleichungen

Bearbeiten

Das wichtigste Beispiel der Strömungsmechanik ist das System der Navier-Stokes-Gleichungen

 
 
 

Als Vereinfachung ergeben sich bei Vernachlässigung der Terme 2. Ordnung die Euler-Gleichungen, bei denen die Flüsse   nicht von   abhängen, also ein System erster Ordnung:

 
 
 

Analysis von Erhaltungsgleichungen

Bearbeiten

Erhaltungsgleichungen in differentieller Form sind partielle Differentialgleichungen (PDEs) und lassen sich somit mit deren Maschinerie analysieren. Dazu ist es notwendig, verschiedene Typen von PDEs auseinanderzuhalten. Zunächst bezeichnet man als Ordnung einer PDE den höchsten Grad der auftretenden Ableitungen. Man spricht von einem System, wenn   ein Vektor ist, ansonsten von einer skalaren Gleichung. Sind   und   linear heißt die Gleichung linear, ansonsten nichtlinear. Die Navier-Stokes-Gleichungen sind also ein nichtlineares System zweiter Ordnung, die Euler-Gleichung ein nichtlineares System erster Ordnung. Eine Gleichung heißt instationär, wenn sie eine Zeitabhängigkeit enthält und stationär, wenn das nicht der Fall ist. Sie heißt mehrdimensional, wenn es mehrere Ortsvariablen gibt.

Schließlich teilt man die Gleichungen noch in hyperbolisch, elliptisch, parabolisch und gemischte Formen ein. Grob kann man sagen, dass diese sich durch die Arten der erlaubten Randbedingungen unterscheiden.

  • Elliptische Gleichungen erlauben nur Randbedingungen, keine Anfangsbedingungen. Damit gibt es für elliptische Gleichungen keine Variable, die man einer Zeit zuordnen könnte; Sie entsprechen stationären Zuständen und Störungen der Nebenbedingungen machen sich sofort in der gesamten Lösung bemerkbar.
  • Parabolische Gleichungen erlauben Anfangs- und Randbedingungen. Sie haben eine Zeitabhängigkeit, Störungen der Nebenbedingungen können sich sofort in der kompletten Lösung bemerkbar machen. Ferner haben parabolische Gleichungen eine Glättungseigenschaft: Sind die Anfangsdaten zum Zeitpunkt Null endlich oft stetig differenzierbar, ist es möglich, dass die Lösung für Zeiten größer Null unendlich oft stetig differenzierbar ist.
  • Hyperbolische Gleichungen erlauben ebenfalls Anfangs- und Randbedingungen. Sie haben eine Zeitabhängigkeit, Störungen der Nebenbedingungen breiten sich mit endlicher Geschwindigkeit in ausgezeichnete Richtungen aus. Im Gegensatz zu parabolischen Gleichungen, die Anfangsdaten glätten, können hier auch bei beliebig glatten Anfangsdaten unstetige Lösungen auftreten.

Wir werden hier vor allem hyperbolische Erhaltungsgleichungen betrachten und nebenbei noch parabolische. Aus Gründen der Einfachheit untersuchen wir zunächst ausschließlich eindimensionale Probleme. Dieses Vorgehensweise ist deswegen sinnvoll, weil ein numerisches Verfahren, das schon im eindimensionalen nicht funktioniert, in mehreren Dimensionen erst recht nicht funktioniert. Zunächst betrachten wir die skalare nichtlineare hyperbolische Erhaltungsgleichung mit  :

  auf dem Gebiet   und
  auf  

Entropielösungen

Bearbeiten

Nichtlineare hyperbolische Erhaltungsgleichungen lassen in der Regel mehrere Lösungen zu. Da Erhaltungsgleichungen physikalische Gesetze beschreiben sollen, ist es sinnvoll, aus diesen mit einer Zusatzbedingung jene Lösungen herauszuziehen, die physikalisch sinnvoll sind. Hierbei hat sich die Entropiebedingung als nützlich erwiesen. Nach dem zweiten Hauptsatz der Thermodynamik ist die Gesamtentropie eines Systems nichtsinkend. Ferner steigt die Entropie eines Gases, wenn man sich über einen physikalischen Schock bewegt. Dies führt zur Entropiebedingung

 

  ist die sogenannte Entropiefunktion und   der Entropiefluss. An dieses Entropie-Entropiefluss-Paar stellt man die Bedingung, dass   glatt und konvex ist,   glatt und ferner die Beziehung

 

gilt. Dann nennt man eine Lösung, die die obige Ungleichung erfüllt, Entropielösung der Erhaltungsgleichung. Analog kann man eine schwache Form herleiten:

 

Eine schwache Lösung, die die obige Gleichung für alle zulässigen Entropiepaare erfüllt, heißt schwache Entropielösung der Erhaltungsgleichung. Diese Lösung ist dann sogar eindeutig und hängt in der  -Norm stetig von den Anfangsdaten ab!

Die Lösung der verschwindenden Viskosität

Bearbeiten

Eine alternative Herangehensweise ist die Betrachtung der Gleichung

 

für  . Diese Gleichung ist parabolisch und hat eine eindeutige Lösung  . Im Grenzwert   ergibt sich die Lösung der verschwindenen Viskosität (vanishing viscosity solution). Diese erfüllt die Entropiebedingung.

Erste Aussagen zur Numerik

Bearbeiten

Dazu sei ein Gitter aus Gitterzellen   gegeben mit   und diskrete Zeitpunkte  . Die Schrittweiten beschreiben wir mit   und  . Mit   bezeichnen wir die numerische Approximation an die exakte Lösung   auf  .

Die zentrale physikalische Eigenschaft die Erhaltungsgleichungen modellieren, ist offensichtlich die Erhaltung. Das numerische Verfahren muss dies berücksichtigen, wenn Hoffnung bestehen soll, eine physikalisch sinnvolle Approximation zu berechnen. Wir werden dies später im Satz von Lax-Wendroff präzisieren, aber zunächst die Definition:

Definition: Konservatives Verfahren

Bearbeiten

Ein numerisches Verfahren zur Lösung der Erhaltungsgleichung   heißt konservativ, wenn es die Form

 

hat. Die Funktion   wird auch als numerische Flussfunktion bezeichnet. Hintergrund ist, dass die numerische Näherung an die Erhaltungsgröße gegeben ist durch

 

Sind die Flüsse an den Randpunkten also physikalisch sinnvoll, erhält das numerische Verfahren die Erhaltungsgröße. Entscheidend für ein sinnvolles numerisches Verfahren ist die Wahl der Flüsse. Die Definition ist offensichtlich auch erfüllt, wenn wir die Nullfunktion als Fluss wählen, eine gute numerische Approximation ist davon aber nicht zu erwarten. Um dem abzuhelfen, braucht es eine Eigenschaft, die verlangt, dass unsere numerische Flussfunktion   den physikalischen Fluss   approximiert.

Definition: Konsistentes Verfahren

Bearbeiten

Eine Flussfunktion   mit den Argumenten  ,   bis   heißt konsistent, wenn die beiden folgenden Eigenschaften erfüllt sind:

  1.  
  2.   ist Lipschitz-stetig in allen Komponenten.

Ausgestattet mit diesen beiden Eigenschaften lässt sich bereits ein wichtiges Resultat beweisen.

Satz von Lax-Wendroff

Bearbeiten

Da nichtlineare Erhaltungsgleichungen in der Regel mehrere schwache Lösungen zulassen, ist die Frage der Konvergenz nicht einfach zu beantworten. Es kann sogar sein, dass ein numerisches Verfahren konvergiert, aber nicht gegen eine schwache Lösung der Erhaltungsgleichung. Mit diesem Problem beschäftigt sich der Satz von Lax-Wendroff.

Satz von Lax-Wendroff (Lax und Wendroff, 1960)

Gegeben sei eine Folge von Gittern mit festen Maschenweiten   und   mit   für  . Gegeben sei ein konsistentes und konservatives numerisches Verfahren, die entsprechende numerische Lösung auf Gitter   sei  . Sei ferner die totale Variation von   beschränkt in dem Sinne, dass zu jedem   ein   existiert, so dass

  für alle  

und das Supremum über alle Unterteilungen   der reellen Achse genommen wird. Schliesslich konvergiere   gegen eine Funktion   in folgendem Sinne: Auf jeder beschränkten Menge   gelte

 , für  

Dann ist   eine schwache Lösung der Erhaltungsgleichung.

Beweis

Zu zeigen ist, dass   die schwache Form der Erhaltungsgleichung erfüllt:

  für alle  

Wir betrachten zunächst das feste Gitter   und vernachlässigen zur Vereinfachung der Notation die zugehörigen Indizes. Die numerische Lösung   ist für alle  ,   gegeben durch die konservative Vorschrift

 

Dies multiplizieren wir mit einer beliebigen Testfunktion   und erhalten mit der Bezeichnung   für die Einschränkung von   auf die Zelle   und das Zeitintervall  :

 

Summieren wir die obige Formel über alle   und   auf, ergibt sich

 

Diese Summen schreiben wir nun so um, dass die Differenzen nun von Testfunktionen genommen werden. Hierbei ist zu beachten, dass die Summen nur formal unendlich sind, da die Testfunktionen einen kompakten Träger haben und damit für   und für   Null sind. Was wir also tun ist, Summen der Form   umzuschreiben als  . Damit ergibt sich

 

und nach weiteren Umformungen

 

Diese Gleichung gilt auf jedem Gitter  . Wir betrachten nun den Grenzwert  . Definieren wir Anfangsdaten beispielsweise derart, dass  , so konvergiert die rechte Seite aufgrund der Glattheit von   und der Konvergenz von   gegen   gegen  . Der erste Term in der Summe auf der linken Seite konvergiert ebenfalls wegen der Glattheit von   und der Konvergenz von   gegen  .

Der zweite Term in der Klammer enthält auch die Flussfunktion und ist deswegen etwas weniger glatt und schwieriger. Die Lipschitz-Stetigkeit von  , gemeinsam mit der beschränkten Variation von   ist allerdings ausreichend, um die Konvergenz von   gegen   fast überall zu garantieren. Damit ergibt sich dann auch für diesen Term Konvergenz gegen das gewünschte Integral   und damit der Beweis des Satzes.

Entropiebedingung

Bearbeiten

Das Problem des Satzes von Lax-Wendroff ist, dass die schwache Lösung der Erhaltungsgleichung nicht eindeutig sein muss und er macht keine Aussagen darüber, zu welcher dieser Lösungen man konvergiert. Noch genauer kann es sein, dass unterschiedliche Folgen von Gittern zu unterschiedlichen schwachen Lösungen konvergieren. Ebensowenig sagt er etwas über den Abschneidefehler aus, mit dem wir in Anbetracht von endlich großen Maschenweiten konfrontiert sind oder liefert überhaupt eine Konvergenzaussage. Für das erste Problem ist die Lösung die so genannte Entropiebedingung. Im kontinuierlichen Fall stellt diese eine eindeutige Lösung bereit und auch im diskreten Fall lässt sich eine Entropiebedingung angeben. Gegeben sei also ein Entropiepaar  ,   und dazu ein diskretes Analogon  , das zu   konsistent ist (im bekannten Sinne). Dann ist die diskrete Entropiebedingung:

 

Erfüllt die numerische Flussfunktion diese Bedingung lässt sich der Beweis für den Satz von Lax-Wendroff in analoger Weise durchführen und zeigen, dass die Grenzfunktion dann die schwache Entropielösung ist.

Stabilität

Bearbeiten

Wir haben nun Bedingungen an ein numerisches Verfahren an der Hand, die im Falle der Konvergenz die Konvergenz gegen eine physikalisch sinnvolle Lösung garantieren. Was fehlt, ist ein Konvergenzbeweis und dafür ist Stabilität des Verfahrens notwendig. Genauer liefert Konsistenz, dass der lokale Abschneidefehler   bei verschwindender Zeitschrittweite gegen Null geht. Was zur Konvergenz noch fehlt ist dann eine Eigenschaft die besagt, dass sich die bereits gemachten Fehler nicht aufschaukeln. Dies bedeutet, dass der durch das numerische Verfahren produzierte Fehler bei gegebenen Schrittweiten   und   für   in der entsprechenden Norm beschränkt bleibt.

Zur numerischen Lösung   zum Zeitpunkt   gehöre der Fehler  . Das numerische Verfahren   liefert dann eine Näherung   zum Zeitpunkt   mit Fehler  . Also haben wir

 

Der zweite Term ist der lokale Abschneidefehler und Stabilität ist die Eigenschaft, die den ersten Term kontrolliert. Wir betrachten zunächst Fall einer linearen Gleichung und dafür ein lineares numerisches Verfahren. Dann ist

 

und es ist ausreichend die Wirkung des numerischen Verfahrens auf den Fehler   zu betrachen. Bleibt der Fehler beschränkt, lässt sich Konvergenz zeigen (Äquivalenzsatz von Lax).

Von-Neumann-Stabilitätsanalyse

Bearbeiten

Wir betrachten zunächst Stabilität in der  -Norm. Die theoretische Rechtfertigung dieser Form der Stabilitätsanalyse basiert auf linearen räumlich periodischen Problemen, für die sich dann sehr starke Aussagen treffen lassen. Allerdings lässt sich die Technik auch sehr gut auf andere Probleme anwenden und liefert sehr gute praktische Ergebnisse, womit sie das Standardanalyseverfahren für Stabilität von numerischen Verfahren für zeitabhängige PDEs darstellt. Gegeben sei auf dem Intervall   die Gleichung

  mit  .

mit periodischen Randbedingungen und die auf ganz   periodisch fortgesetzte Lösung  . Das numerische Verfahren definiert dann eine Evolution des Fehlers in der Zeit und lässt sich über eine Amplifikationsmatrix darstellen, deren Dimension allerdings mit kleiner werdender Maschenweite ansteigt. Die Idee ist nun, den periodischen Fehler zum Zeitpunkt   im Diskretisierungspunkt   in eine Fourier-Reihe

 

zu entwickelt. Dies diagonalisiert die Amplifikationsmatrix und wir können dann die Entwicklung der einzelnen Fourierkoeffizienten in der Zeit analysieren. Nach der parsevalschen Ungleichung gilt  . Die  -Stabilitätsbedingung reduziert sich dann darauf, dass das numerischer Verfahren genau dann stabil ist, wenn der Spektralradius der Amplifikationsmatrix   betragsmäßig kleiner gleich eins ist.

Beispiel

Der einfachste Fall ist die lineare Advektionsgleichung

  mit  

In der Zeit nehmen wir das explizite Euler-Verfahren   gekoppelt mit zentralen Differenzen auf einem äquidistanten Gitter im Raum. Der zweite Term wird also mittels

 

approximiert. Insgesamt ergibt sich

 

was auch die Entwicklung der Fehler definiert und auch jedes einzelnen Terms der Fourier-Reihenentwicklung. Betrachten wir den j-ten Summanden, so ergibt Einsetzen in die obige Formel und Division durch   mit  :

 

Die Amplifikationsmatrix ist nun gegeben durch

 

Das Verfahren ist  -stabil, falls   für alle  , was hier nicht der Fall ist, da   Damit ist das Verfahren unabhängig von der Wahl der Schrittweiten instabil. Dieses Verhalten beobachteten die Mitarbeiter des Manhattan-Projekts, was von Neumann zur Entwicklung der Stabilitätsanalyse führte. Wird im Raum das Upwind-Verfahren

 

benutzt, so ergibt sich die Courant-Friedrichs-Lewy-Bedingung

 ,

also bedingte Stabilität. Im nichtlinearen Fall oder im Falle variabler Koeffizienten kann das Verfahren durch Linearisierung und Einfrieren der Koeffizienten angewandt werden. Allerdings liefert die Analyse im Allgemeinen Fall nur noch eine notwendige Bedingung für Stabilität, in Spezialfällen auch eine hinreichende. Unter Verwendung eines leicht anderen nichtäquivalenten Stabilitätsbegriffes wird manchmal eine andere Bedingung an Stabilität gestellt:

 .

Ein allgemeines Verfahren zur vollständigen Stabilitätsanalyse nichtlinearer Gleichungen ist nicht bekannt.

TV-Stabilität

Bearbeiten

Im nichtlinearen Fall braucht man in der Regel nichtlineare numerische Verfahren und dann funktioniert die obige Fehlerdarstellung nicht. Ist das numerische Verfahren eine Kontraktion, also gilt

 

so ergibt sich immerhin

 

womit sich Stabilität zeigen lässt. Diese Eigenschaft gilt insbesondere für die später erwähnten monotonen Verfahren, die aber maximal erster Ordnung sind. Entsprechend brauchen wir einen besseren Stabilitätsbegriff, nämlich die TV-Stabilität, bei der es um Konvergenz in der  -Norm geht. Hintergrund ist, dass Mengen der Form

  und  

in   kompakt sind.

Definition TV-stabil

Wir nennen ein numerisches Verfahren nun TV-stabil, falls die numerischen Lösungen für alle   in einer Menge der obigen Form sind mit M und R fix.

Hierzu gilt folgender Satz:

Satz

Gegeben sei ein konservatives numerisches Verfahren mit einem Lipschitz-stetigen numerischen Fluss und für die Anfangsdaten   existiere ein   und ein  , so dass für alle n und alle   mit   und   gelte:

 

Dann ist das numerische Verfahren TV-stabil.

Beweis

Siehe LeVeque, Seite 250.

Hiermit können wir nun ein Konvergenzresultat beweisen:

Satz

Die numerische Lösung U sei von einem TV-stabilen, konservativen und konsistenten numerischen Verfahren mit fester Zeitschrittweite   generiert. Dann konvergiert die Lösung für   in der  -Norm gegen eine schwache Lösung der Erhaltungsgleichung.

Beweis

Wir nehmen an, das die Aussage falsch wäre, es gäbe also eine Nullfolge  , so dass die zugehörige Folge von numerischen Lösungen nicht gegen eine schwache Lösung konvergiert. Da das Verfahren TV-stabil ist, liegt die Folge   der numerischen Lösungen in einer in   kompakten Menge und entsprechend existiert eine  -konvergente Teilfolge, die die Voraussetzungen des Satzes von Lax-Wendroff erfüllt. Diese konvergiert also gegen eine schwache Lösung der Erhaltungsgleichung.

Monotone Verfahren

Bearbeiten

Eine Eigenschaft der hier betrachteten Erhaltungsgleichungen ist, dass der Lösungsoperator invers monoton ist: gilt für zwei Anfangsdaten  , dann gilt für die entsprechenden Lösungen  . Ein konservatives Verfahren, das die entsprechende diskrete Eigenschaft hat, nennt man monoton. Genauer heisst das Verfahren monoton, falls für die numerische Lösung aus

  für alle   folgt, dass
  für alle  

gilt. Dies ist äquivalent mit

  für alle  ,

wobei   mittels   definiert ist.

Monotone Verfahren haben den grossen Vorteil, dass sie umfassende Konvergenzaussagen erlauben. Harten, Hyman und Lax bewiesen 1976, dass monotone Verfahren im Falle der Konvergenz gegen die schwache Entropielösung konvergieren. 1980 bewiesen dann Crandall und Majda folgendes Resultat:

Satz

Die durch ein monotones Verfahren generierte Lösung konvergiert in   gegen die schwache Entropielösung der Erhaltungsgleichung für   gegen Null mit   fix. Dieses Resultat gilt auch im mehrdimensionalen.

Kröner und Ohlberger bewiesen 2000 eine grundlegende Fehlerabschätzung für monotone Verfahren. Um 2000 wurde von verschiedenen Autoren folgendes mehrdimensionale Resultat erarbeitet, was es auf der einen Seite tatsächlich erlaubt, a priori Fehlerschätzer zu erarbeiten, aber auf der anderen Seite das große Problem monotoner Verfahren vor Augen führt.

Satz

Sei   in   und   eine schwache Entropielösung der Erhaltungsgleichung. Sei   eine durch ein monotones numerisches Verfahren mit dem expliziten Eulerverfahren als Zeitintegration generierte numerische Lösung, mit einem in bestimmter Weise beschränkten Zeitschritt. Sei   eine bestimmte Teilmenge von  . Dann existiert eine Konstante  , so dass

 

Im eindimensionalen Fall gilt die Aussage mit  , diese Konvergenzrate ist optimal.

Numerisch beobachtet man häufig in O(h)-Verhalten, nichtsdestotrotz ist diese Konvergenzrate zu niedrig. Der Grund für die niedrige Rate liegt in der Eigenschaft der Monotonie. Diese ist sehr weitgehend, was die starken Konvergenzaussagen erst erlaubt, allerdings auch eine extrem starke Einschränkung an die Verfahren darstellt. Ein monotones Verfahren hat nebenbei noch die folgenden Eigenschaften:

  1. Die  -Norm der numerischen Lösung ist mit der Zeit nichtsteigend ( -kontraktiv).
  2. Damit gilt auch: Die Totale Variation ist nicht steigend (TVD).
  3. Und daraus folgt: Monotonie der Lösung bleibt erhalten (Monotonieerhaltend).

Diese Eigenschaften bilden eine Hierarchie von Verfahrensklassen, die echte Teilmengen bilden. Es liegt also auf der Hand, die Anforderungen an das Verfahren abzuschwächen, um eine höhere Konvergenzordnung zu erzielen.

Monotonieerhaltende Verfahren

Bearbeiten

Die Lösung   einer skalaren hyperbolischen Erhaltungsgleichung hat die Eigenschaft, dass aus Monotonie der Anfangsdaten Monotonie der Funktionen   für alle Zeiten folgt. Das diskrete Analogon davon sind die monotonieerhaltenden Verfahren, die bei monotonen Anfangsdaten   monotone Näherungslösungen   für alle   liefern.

Godunow bewies in den 1950er Jahren, dass lineare Einschrittverfahren zweiter Ordnung für die lineare Advektionsgleichung nur in uninteressanten Spezialfällen monotonieerhaltend sein können. Später wurde die Aussage von anderen auch für lineare Mehrschrittverfahren bewiesen. Wir zeigen den Beweis von Godunow.

Jedes lineare Einschrittverfahren lässt sich in der Form

 

schreiben (das Verfahren ist eine lineare Abbildung von   auf  . Zunächst ein Hilfssatz

Satz

Das obige Verfahren ist genau dann monotonieerhaltend, wenn

  für alle  .
Beweis

ObdA sei   und   nichtfallend.

Hinrichtung: Sei das Verfahren monotonieerhaltend, aber es existiere ein p mit   Wir wählen nun Anfangsdaten:

 ,   und  ,  

Dann ist   und damit   nicht nichtfallend, was einen Widerspruch ergibt.

Rückrichtung: Seien alle  . Dann ist für alle j:

 

also ist das Verfahren monotonieerhaltend.

Satz von Godunow

Bearbeiten

Wir können nun den Satz von Godunow beweisen:

Lineare Einschrittverfahren zweiter Ordnung für die lineare Advektionsgleichung können nicht monotonieerhaltend sein, es sei denn  

Beweis

Sei   und die entsprechenden numerischen Anfangsdaten  . Die exakte Lösung ist ein Polynom zweiten Grades und ergibt sich durch Advektion der Anfangsdaten:

 

In Godunows Sinne bedeutet, dass ein Verfahren zweiter Ordnung diese Lösung exakt reproduziert, es muss also gelten:

 

Nehmen wir an, das Verfahren sei Monotonieerhaltend, also   Wegen   folgt

  für alle  .

Für   und größer Null wählen wir nun   mit  . Dann gilt

 

und damit ein Widerspruch.

Ein etwas stärkerer Beweis stammt von Roe aus dem Jahr 1986, siehe Wesseling, Seite 346 für eine Skizze. Der Satz gilt ebenfalls für lineare Mehrschrittverfahren. Die Aussage ist also: Lineare Verfahren können nicht zweite Ordnung erreichen und um hohe Ordnung zu erreichen sind nichtlineare Verfahren unabdingbar. Hier haben sich insbesondere die TVD-Verfahren als nützlich erwiesen.

TVD-Verfahren

Bearbeiten

Ein numerisches Verfahren hat die TVD-Eigenschaft (von Total Variation Diminishing), wenn für alle  

 

gilt. Diese Eigenschaft gilt analog für die kontinuierliche Lösung und so ist es sinnvoll, dies auch vom numerischen Verfahren zu verlangen. Ist die totale Variation zum Anfangszeitpunkt beschränkt, so ist sie das bei einem TVD-Verfahren auch für alle späteren Zeitpunkte. Damit gilt nach den bereits gebrachten Aussagen zur TV-Stabilität: Konservative Verfahren mit Lipschitz-stetigem Fluss die die TVD-Eigenschaft haben sind TV-stabil und damit konvergent gegen eine schwache Lösung der Erhaltungsgleichung.

Die TVD-Eigenschaft stellt eine Bedingung an die kombinierte Raum-Zeit-Diskretisierung dar. Hierzu gibt es nützliche Sätze von Harten:

Satz

Gegeben sei das explizite Drei-Punkt-Verfahren

 

wobei   und   von   abhängen. Wenn die Koeffizienten für alle   die Bedingung  ,   und   erfüllen, ist das Verfahren TVD.

Beweis

Wir betrachten

 

Nach Voraussetzung ist dann

 

Summation über   ergibt genau die totale Variation zum Zeitpunkt  :

 

und damit die Aussage des Satzes.

Hierbei ist anzumerken, dass die Bedingung   eine Art CFL-Bedingung ist, damit gibt es also keinen Wiederspruch zur Konvergenz des Verfahrens. Für implizite Verfahren ergibt sich die Aussage

Satz

Gegeben sei das implizite Drei-Punkt-Verfahren

 

wobei   und   von   abhängen, während   und   von   abhängen. Wenn die Koeffizienten die Bedingung   und   erfüllen, ist das Verfahren TVD.

Der Beweis teilt das Verfahren in seinen expliziten und den impliziten Anteil auf und nutzt für den expliziten Anteil die vorherige Aussage. Die zusätzlichen Bedingungen sind also dazu da, die totale Variation des impliziten Anteils ausreichend zu beschränken.

Beide Sätze lassen sich auf breitere Diskretisierungssterne erweitern. Ferner kann der explizite Satz genutzt werden, um Bedingungen an TVD-Runge-Kutta-Verfahren zu ermitteln und solche zu konstruieren. Dazu später mehr.

Wir wissen aufgrund des Satzes von Godunow, dass lineare monotonieerhaltende Verfahren nicht von zweiter Ordnung sein können, und zwar global. Entsprechend sucht man den Ausweg über nichtlineare Verfahren. Leider stellt sich heraus, dass auch diese Grenzen haben: Osher bewies 1984, dass TVD-Verfahren an glatten Extrema zu erster Ordnung degenerieren und Goodman und LeVeque bewiesen 1985, dass ein TVD-Verfahren in mehreren Dimensionen maximal erste Ordnung haben kann. Nichtsdestotrotz sind TVD-Verfahren, die höhere Ordnung über Limiter oder über ENO-Konstruktionen erzielen, sehr nützliche Hilfsmittel.

Jameson-Schmidt-Turkel-Verfahren (JST-scheme), Teil 1

Bearbeiten

Neben den erwähnten Limitern und ENO-Verfahren besteht eine andere Möglichkeit darin, die instabile zentrale Diskretisierung als Ausgangspunkt zu nehmen und über nichtlineare Zusatzterme zu stabilisieren. Jameson, Schmidt und Turkel entwickelten hierzu das nach ihnen benannte JST-Verfahren, zum Einsatz in nichtlinearen Mehrgitterverfahren, da Limiter in diesem Kontext die Konvergenz verlangsamen können. Betrachten wir das semidiskrete Verfahren

 

mit  . Der Term künstlicher Dissipation wird definiert mittels

 

Der Koeffizient   wird dabei so konstruiert, dass er in Regionen mit glatter Lösung von der Ordnung   ist, während   dort O(1) ist, so dass beide Terme dort von der Ordnung   sind, womit sich in glatten Regionen eine niedrige Diffusion ergibt. In der Nähe eines Schocks sollte   sein und damit das Verfahren erster Ordnung.

Das originale JST-Verfahren war nicht TVD, funktionierte aber gut in der Praxis, seitdem wurde es erheblich verbessert. Zunächst ein Einschub.

Positive Verfahren/LED-Verfahren

Bearbeiten

Wie bereits erwähnt sind TVD-Verfahren in mehreren Dimensionen maximal erster Ordnung, es ist also zu überlegen, ob man die TVD-Eigenschaft noch einmal relaxiert oder eine zumindest im eindimensionalen äquivalente Eigenschaft sucht. Dies führt auf positive Verfahren (Spekreijse 1987), beziehungsweise Local extremum diminishing schemes (LED, Jameson). Sei also das semidiskrete Verfahren

 

gegeben. Befindet sich im   ein lokales Maximum, so ist damit  . Dann gilt   und damit, dass ein lokales Maximum nicht weiter wächst, wenn alle  . Dies gilt offensichtlich auch im mehrdimensionalen, für unstrukturierte Gitter und für lokale Minima. Ein Verfahren mit der ersten Eigenschaft heisst LED, ein Verfahren mit der zweiten positiv. Manchmal redet man statt LED auch von einem diskreten Maximumsprinzip. In 1-D folgt aus der LED-Eigenschaft die TVD-Eigenschaft (Tadmor 1988), nicht aber im mehrdimensionalen. Spekrijse nutzte positive Verfahren, um mehrdimensionale MUSCL-Verfahren zweiter Ordnung zu konstruieren, ein Konvergenzbeweis im mehrdimensionalen steht aber noch aus.

Für das JST-Verfahren gilt folgender Satz:

Satz

Es gelte immer dann, wenn   oder   ein lokales Extremum ist, für die Koeffizienten des JST-Verfahrens:

  und  

wobei der Koeffizient   gegeben sei durch

  für   und
  für  .

Dann ist das Verfahren LED.

Beweis

Sei   ein Maximum und damit die Koeffizienten   und   Null. Dann ergibt sich das semidiskrete Verfahren zu

 

und damit, dass alle Koeffizienten positiv sind.

Damit wissen wir, wie die Koeffizienten   und   zu konstruieren sind.

JST-Verfahren, Teil 2

Bearbeiten

Das JST-Verfahren wurde im Laufe der Jahrzehnte immer wieder angepasst, die aktuelle Version ist folgende, bei der die Koeffizienten definiert werden als:

  und  

An einem Extremum sollte   also Eins sein. Dazu sei

  und  

mit

 

  ist dabei eine positive ganze Zahl. Wenn   und   entgegengesetzte Vorzeichen haben (entspricht einem Extremum), ist  , ansonsten kleiner.

Systeme von Gleichungen

Bearbeiten

Bisher haben wir ausschliesslich skalare Gleichungen betrachtet, manchmal in mehreren Dimensionen. Nun sind die meisten praxisrelevanten Gleichungen Systeme, die Frage ist also, welche Aussagen noch übrig bleiben. Zunächst gilt der Satz von Lax-Wendroff weiterhin. Schon schwieriger wird es bei der Entropiebedingung, diese muss man bei Systemen auf so genannte genuinely nonlinear systems einschränken, was in etwa einer Konvexitätsbedingung an die Flussfunktion entspricht. Damit ist noch keine Konvergenz des numerischen Verfahrens gezeigt, was der Situation im kontinuierlichen entspricht. Ausser für Spezialfälle (starke Einschränkungen an die Anfangsdaten oder das System) gibt es keine Aussagen zu Lösungen für nichtlineare Systeme.

Für Konvergenz des numerischen Verfahren brauchen wir Stabilität und wir wollen die Probleme an der totalen Variation erläutern. Für lineare Systeme ist es möglich die totale Variation so zu definieren, dass das diagonalisierte System betrachtet wird und dann die totale Variation die Summe der totalen Variationen der einzelnen Lösungskomponenten. In jeder einzelnen Komponente ist die skalare Theorie anwendbar und damit auch die gesamte totale Variation nichtsteigend.

Für nichtlineare Systeme funktioniert das nicht mehr, da sich die Diagonalisierungsmatrizen mit der Lösung ändern. In der Tat kann die totale Variation schon für die kontinuierliche Lösung ansteigend sein. Ein Beispiel sind zwei kollidierende Schockwellen bei den Euler-Gleichungen, die nach der Kollision zwei auseinanderlaufende Schockwellen unter Erhöhung der totalen Variation erzeugen. Gerade bei den Euler-Gleichungen in mehreren Dimensionen tritt noch das Problem auf, dass ein Eigenwert mehrfach vorkommt, bei der Erweiterung der MHD-Gleichung ist dann sogar die schwache Entropielösung des Riemann-Problems nicht mehr eindeutig. Einer der wenigen Konvergenzbeweise für nichtlineare Systeme ist Glimm's Random Choice Method.

Zeitintegration

Bearbeiten

Wir betrachten ab jetzt die semidiskrete Gleichung

 

wobei   ein Gittervektor ist und   aus der räumlichen Diskretisierung stammt. Es handelt sich also um ein nichtlineares System gewöhnlicher Differentialgleichungen, die mit einer der vielen Verfahren für Anfangswertprobleme gelöst werden können.

TVD-Runge-Kutta-Verfahren

Bearbeiten

Eine Frage ist, unter welchen Bedingungen die Zeitdiskretisierung die TVD-Eigenschaft der räumlichen Diskretisierung erhält bzw. wann das Gesamtverfahren TVD ist. Wir hatten bereits Bedingungen an Verfahren bei Nutzung des expliziten Euler-Verfahrens hergeleitet. Diese Technik kann genutzt werden, um zu zeigen, dass das Heun-Verfahren die TVD-Eigenschaft erhält. Dieses ist ein Runge-Kutta-Verfahren 2. Ordnung der Form:

  1.  
  2.  

Dies lässt sich als Folge zweier Schritte des expliziten Euler-Verfahrens schreiben:

  1.  
  2.  
  3.  

Wenn nun das Verfahren, dass durch   TVD ist (etwa die Bedingungen des Satzes von Harten für das explizite Euler-Verfahren erfüllt), sind die einzelnen Schritte TVD und es gilt

 

und damit

 

das Gesamtverfahren ist dann also TVD und hat dasselbe Stabilitätsgebiet. Dieses ist sogar bei dritter Ordnung möglich (Shu und Osher 1988):

  1.  
  2.  
  3.  

Strong stability preserving methods

Bearbeiten

Allgemeiner heißt ein Verfahren strong stability preserving (SSP), wenn es genau wie das explizite Euler-Verfahren in der  -Norm stabil ist, also

 

unter einer CFL-artigen Stabilitätsbedindung. Gottlieb, Shu und Tadmor bewiesen 2001, dass in der Klasse der SSP-Runge-Kutta-Verfahren  -ter Ordnung mit   Stufen keine im Vergleich zum expliziten Euler relaxierte Stabilitätsbedingung möglich ist. Das Verfahren von Heun und das genannte Verfahren 3. Ordnung sind in dieser Hinsicht also optimal.

Stationäre Rechnungen, explizite Zeitintegration

Bearbeiten

Bei stationären Rechnungen ist es nicht von Belang, wie die Zwischenergebnisse aussehen sondern nur, ob die numerische Näherung an die stationäre Lösung sinnvoll ist, entsprechend kann dann auf die TVD-Eigenschaft der Zeitintegration verzichtet werden. Eine zweite Idee ist es, die konvektiven und die diffusiven Terme mit unterschiedlichen Verfahren zu behandeln. Zunächst betrachten wir explizite Zeitintegration. Die Idee ist es, die Koeffizienten der Runge-Kutta-Verfahren so zu wählen, dass das Stabilitätsgebiet maximal wird. Jameson schlägt folgendes Verfahren vor. Die rechte Seite wird aufgespalten in den konvektiven und den diffusiven Anteil

 .

Dann wird das Runge-Kutta-Verfahren definiert mittels

 
 , für  
 

Hierbei ist  ,   und ferner   sowie  

Die Koeffizienten   für den konvektiven Anteil werden nun so gewählt, dass das Stabilitätsgebiet entlang der imaginären Achse vergrößert wird. Für ein 5-stufiges Schema schlägt Jameson folgendes vor:

 

Die Koeffizienten   für den diffusiven Anteil werden so gewählt, dass das Stabilitätsgebiet entlang der negativen reellen Achse vergrößert wird. Für ein 5-stufiges Schema schlägt Jameson folgendes vor:

 .

Im Vergleich zum klassischen Runge-Kutta-Verfahren ist das Stabilitätsgebiet auf der negativen reellen Achse etwa dreimal so groß.

Implizite Verfahren

Bearbeiten

Als zweite Klasse wollen wir hier noch implizite Verfahren erwähnen, die dann zum Einsatz kommen können, wenn explizite Verfahren aufgrund ihres Stabilitätsgebiets einen zu kleinen Zeitschritt haben. In der numerischen Strömungsmechanik sind die betrachteten Gleichungen in der Regel sehr steif, interessant sind also solche Verfahren, ein möglichst großes Stabilitätsgebiet haben. Das ist für A-stabile Verfahren erfüllt und so sind BDF-2, die implizite Mittelpunktsregel und das Crank-Nicholson-Verfahren sehr beliebt. Diese Verfahren sind alle zweiter Ordnung, was nach der Zweite Dahlquist-Barriere auch das Optimum für A-stabile lineare Mehrschritt-Verfahren ist. Bei Verfahren höherer Ordnung gibt es im wesentlichen zwei Ansätze. Entweder werden BDF-Verfahren höherer Ordnung genommen, in der Hoffnung, dass die Exlusion der imaginären Achse nicht so schlimm ist. Sowohl BFD-3 als auch BDF-4 wurden erfolgreich eingesetzt.

Die zweite Möglichkeit ist, implizite Runge-Kutta-Verfahren zu nutzen. Diese können grundsätzlich A-stabil sein, allerdings ist es dann erforderlich, ein nichtlineares Gleichungssystem in jeder Stufe zu lösen. Bijl et al (2002) haben so genannte ESDIRK-Verfahren genutzt und gezeigt, dass diese Verfahren kompetitiv sind und besser als die vergleichbaren BDF-Verfahren sein können.

Gemein ist allen Verfahren, dass nichtlineare Gleichungssysteme gelöst werden müssen.

Nichtlineare Gleichungssystemlöser

Bearbeiten

Zur Lösung von Systemen nichtlinearer Gleichungssysteme wurden zahlreiche Gleichungssystemlöser entwickelt. Wir werden hier die Löser betrachten, die zur numerischen Lösung der Euler- und Navier-Stokes-Gleichungen eingesetzt werden.

Newton-Verfahren

Bearbeiten

Das Newton-Verfahren (auch Newton-Raphson-Verfahren, aber auch Simpson trug wesentlich zur Entwicklung bei) löst Nullstellengleichungen   und beruht auf einer Linearisierung der Funktion  . Die Iterationsvorschrift ist gegeben durch

 

wobei   die Jacobi-Matrix von   an der Stelle   ist. Da die Inverse in der Regel nicht zur Verfügung steht, wird diese Iteration in zwei Schritten durchgeführt, nämlich der Lösung eines linearen Gleichungssystems und dem Lösungsupdate:

  1. Löse  
  2. Berechne  

Das Newton-Verfahren konvergiert lokal quadratisch. Die konkrete Form des Jacobi-Matrix hängt insbesondere vom Zeitintegrationsverfahren ab und wie dieses formuliert ist. Betrachten wir das  -Verfahren ( ):

 

Nach einer Finite-Volumen-Diskretisierung erhalten wir für das komplette System auf dem gesamten Rechengitter

 

wobei   die Diagonalmatrix bezeichnet, bei der auf der Diagonalen die entsprechenden Volumina zu finden sind. Dies lässt sich auf verschiedene Weisen in eine Nullstellengleichung für den Vektor der Unbekannten   umformen:

 

Die Struktur der Jacobi-Matrix ist (auch für andere Zeitintegrationsverfahren) immer die, dass sie aus der Summe der Ableitung des Vektors   und der Ableitung der numerischen Flussfunktion besteht, wobei die Faktoren   und   nur bei einem der Terme auftauchen. Es ergibt sich dann für das lineare Gleichungssystem in diesem Fall

 

Hierbei hat man noch die Auswahl, in welchem Variablensatz man das Newton-Verfahren durchführen will (konservative oder primitive Variablen etwa). Bei konservativen Variablen ist   eine Diagonalmatrix, dafür sind die Ableitungen der numerischen Flussfunktion in primitiven Variablen häufig etwas einfacher.

Die obige Matrix ist aufgrund der Eigenschaften von   unsymmetrisch, indefinit, keine M-Matrix, für hinreichend große   schlecht konditioniert. Sie ist eine Blockmatrix, wobei die Größe der Blöcke der Anzahl der Variablen des Systems der Strömungsgleichungen entspricht (je nach Raumdimension 3-5).

Die Aussage der lokal quadratischen Konvergenz gilt für den Fall, wo die Ableitung exakt ausgerechnet wurde und das lineare Gleichungssystem in jedem Schritt exakt berechnet wird. Wird auf eins oder beides verzichtet, so heißt das Verfahren inexakt und verliert man die quadratische Konvergenz, dafür kann die benötigte Rechenzeit drastisch reduziert werden. Erstes Beispiel ist das vereinfachte Newton-Verfahren, bei dem die Ableitung nur im ersten Schritt exakt berechnet wird. Dieses Verfahren kann immerhin superlineare konvergieren. Bei inexakten Newton-Verfahren ist eine quantitative Konvergenzaussage schwierig, grob gilt, je schlechter approximiert wird, desto schlechter ist das Konvergenzverhalten. Wichtigster Faktor ist die Geschwindigkeit, mit der die linearen Gleichungssysteme auf befriedigende Genauigkeit gelöst werden können.

Newton-Krylow-Verfahren

Bearbeiten

Die Variante, bei der das lineare Gleichungssystem mittels eine Krylow-Unterraum-Verfahrens gelöst werden, nennt man Newton-Krylow-Verfahren. Aufgrund der Eigenschaften der Systemmatrix kommen nicht alle Verfahren in Frage, sondern nur solche für allgemeine Systeme, insbesondere GMRES und BiCGSTAB. Nach praktischer Erfahrung ist BiCGSTAB hier die bessere Wahl (Meister, Vömel 2001).

Matrixfreie Verfahren

Bearbeiten

Krylow-Unterraum-Verfahren haben die Eigenschaft, dass sie nie die Matrix   explizit benötigen, sondern immer nur Matrix-Vektor-Produkte. Ist die Matrix, wie in diesem Fall eine Jacobi-Matrix, repräsentiert das Matrix-Vektor-Produkt eine Richtungsableitung und kann damit durch Funktionsauswertungen approximiert werden:

 

für  . Das sich ergebende Verfahren braucht somit deutlich weniger Speicher und kann schneller sein als ein Verfahren bei dem die Matrix-Vektor-Produkte ausgerechnet werden. Als grundlegendes Krylow-Unterraum-Verfahren ist hier GMRES die beste Wahl. Dies liegt daran, dass hier letztlich eine Störung der Matrix vorliegt und damit der aufgespannte Krylow-Unterraum   ebenfalls gestört wird. GMRES ist dabei aufgrund seiner Minimierungseigenschaft robuster als etwa BiCGSTAB gegen solche Störungen: In jedem Schritt wird die Norm   über den aufgespannten Krylow-Unterraum minimiert, GMRES liefert also immer in jedem Schritt etwas sinnvolles.

In unserem konkreten Fall ergibt sich:

 

und damit

 

Pro Matrix-Vektor-Multiplikation muss also eine weitere Auswertung der rechten Seite erfolgen, da   während eines Newton-Schritts konstant ist.

Vorkonditionierung

Bearbeiten

Entscheidender als die Wahl des Krylow-Unterraum-Verfahrens ist die Wahl der Vorkonditionierung. Unterschieden wird zwischen Links- und Rechtsvorkonditionierung (und beidseitiger):

 

entspricht Linksvorkonditionierung und

 ,  

Rechtsvorkonditionierung. Es geht darum, einen Vorkonditionierer zu finden, der

  1. In der Anwendung wenig kostet.
  2. Die inverse von A gut approximiert.

Jedes lineare iterative Verfahren ist als Vorkonditionierer geeignet, da eine Anwendung eines iterativen Verfahrens einer Approximation von   entspricht. Für Strömungsprobleme haben sich hier die unvollständige LU-Zerlegung (ILU), bei der die Matrix A selbst approximiert wird und das symmetrische Gauß-Seidel-Verfahren (SGS,  ) etabliert. SGS ist im supersonischen Bereich extrem gut, da dann die Matrix fast tridiagonal wird und SGS damit ein exakter Löser wäre. Im Bereich kleiner Machzahlen ist dies nicht der Fall und die Performance sinkt. Ferner hat sich herausgestellt, dass Rechtsvorkonditionierung effektiver ist als Linksvorkonditionierung.

Rechtsvorkonditioniertes GMRES

GMRES basiert darauf, mit Hilfe der Arnoldi-Prozedur eine orthogonale Basis des Krylow-Unterraums zu berechnen und dann über das Minimierungsproblem die Lösung auszurechnen. GMRES rechnet im Krylow-Unterraum   und wendet bei der Berechnung der Lösung einmalig den Vorkonditionier nochmal an:

Gegeben  , initialisiere   mit Nullen.

Führe dann den Arnoldi-Prozess durch:

  1. Berechne  . If  , then END.  .  .
  2. For  
    Berechne  .
    Berechne  
    For   do  
     
    Berechne   und  
    Berechne die Norm der Näherungslösung über  . Ist diese kleine genug, STOP.
  3. Definiere  .

Berechne anschliessend die Näherungslösung:

  mit   wie oben. Wenn vorgesehen: RESTART, sonst ENDE.

Im letzten Schritt wird die Lösung als Linearkombination der vorkonditionierten Arnoldi-Basisvektoren   berechnet.

FGMRES

Eine weitere Variante sind Vorkonditionierer, die sich mit jedem Schritt ändern, dazu gehören insbesondere nichtlineare Mehrgitterverfahren. Dazu ist es nötig, das Krylow-Unterraum-Verfahren etwas anzupassen: Oben haben wir die Vektoren   nach ihrer Berechnung nicht weiter abspeichern müssen. Es reichte, die   und den Vorkonditionierer abzuspeichern. Dies ändert sich bei varieblen Vorkonditionierern, bei denen   gilt. Eine Möglichkeit so vorzugehen bietet FGMRES (flexible GMRES) (Saad 1993).

Gegeben  , initialisiere   mit Nullen.

Führe dann den Arnoldi-Prozess durch:

  1. Berechne  . If  , then END.  .  .
  2. For  
    Berechne  .
    Berechne  
    For   do  
     
    Berechne   und  
    Berechne die Norm der Näherungslösung über  . Ist diese kleine genug, STOP.
  3. Definiere  .

Berechne anschliessend die Näherungslösung:

  mit   wie oben. Wenn vorgesehen: RESTART, sonst ENDE.

Der Zusatzaufwand ist hier das Abspeichern der Vektoren   (die  ) müssen eh gespeichert werden). Dafür ist es nun möglich, flexible Vorkondonditionierer einzusetzen. Eine andere Möglichkeit ist GMRES-R von van der Vorst und Vuik, 1994.

Mehrgitter-Verfahren

Bearbeiten

Die wesentliche Alternative zu Newton-Verfahren sind Mehrgitter-Verfahren. Die schnellsten Verfahren zur Berechnung der stationären Euler-Gleichungen stammen aus dieser Klasse. Die Grundidee ist, eine Hierarchie von gröber werdenden Gittern aufzubauen und diese zu nutzen, um eine Approximation des unbekannten Fehlers einer Anfangsnäherung der Lösung zu berechnen.

Lineare Mehrgitter-Verfahren

Bearbeiten

Zunächst sei auf dem feinsten Gitter ein lineares Gleichungssystem   mit regulärer Matrix   gegeben. Eine Iteration des Mehrgitter-Verfahrens sieht dann folgendermaßen aus:

 
if  ,   (Löse exakt auf gröbstem Gitter)
else
  (Vorglättung)
  (Restriktion)
 
Für     (Berechnung der Grobgitterkorrektur)
  (Prolongation der Grobgitterkorrektur)
  (Nachglättung)
end if
End

Was ein sinnvoller Glätter   ist, hängt von der zu lösenden partiellen Differentialgleichung ab. Für viele sind Gauß-Seidel- oder Jacobi-Relaxation geeignet. Da niederfrequente Fehleranteile auf feinen Gittern hochfrequenten Fehleranteilen auf gröberen Gittern entsprechen, ergibt sich mit der Residuumsgleichung

 

ein Problem mit ähnlicher Struktur wie das Ursprungsproblem, allerdings mit kleinerer Dimension. Dies wird rekursiv bis zum gröbsten Gitter wiederholt (  entspricht einem V-Zyklus,   einem W-Zyklus), wo das Gleichungssystem direkt gelöst werden kann.

 
V-Zyklus
 
W-Zyklus

Der berechnete Fehler wird dann sukzessive mittels Prolongation P auf die feineren Gitter rückprolongiert und die ursprüngliche Näherung der Lösung korrigiert. Dies funktioniert bei einem linearen Problem   mit Lösung  , da dann der Fehler   der Näherungslösung   über die Residuumsgleichung   berechnet werden kann. Es folgt noch Nachglättung.

Für viele elliptische und parabolische Gleichungen gibt es umfassende Konvergenztheorie zu linearen Mehrgitterverfahren. Dies ist bei den Euler- und Navier-Stokes-Gleichungen nicht der Fall, sie finden aber manchmal als Vorkonditionierer oder als Löser in Newton-Verfahren Verwendung.

Full Approximation Scheme

Bearbeiten

Auf ein nichtlineares Problem   lässt sich das obige Vorgehen nicht direkt übertragen, da die Residuumsgleichung   gar nicht lösbar sein muss. Deswegen löst man da auf dem groben Gitter statt dessen  , was im linearen Fall äquivalent wäre. Es ergibt sich dann

 
if  ,  
else
 
 
Wähle   und  
 
Für    
 
 
end if
End

  beschreibt dabei eine Approximation an die Lösung und   wird so klein gewählt, dass das entsprechende nichtlineare Gleichungssystem lösbar ist.   entspricht dem Full Approximation Scheme (FAS) von Achi Brandt (1977). In Hackbusch (81, 82) und seinem Buch Multigrid Methods and Applications (85) findet sich eine ausführlichere Beschreibung von Auswahlmöglichkeiten.

Nichtlineares Mehrgitterverfahren von Jameson & al

Bearbeiten

Eine weitere Variante wurde von Jameson und Turkel im Laufe von Jahrzehnten für die stationären Euler-Gleichungen entwickelt. Dabei wird auf Nachglättung verzichtet ( ) und als Vorglätter wird ein Schritt mit dem oben erwähnten speziellen Runge-Kutta-Verfahrens gemacht, das nämlich neben einem grossen Stabilitätsgebiet noch Glättungseigenschaften hat. Da die Restriktion auf ein gröberes Gitter führt, ist es dort möglich, den Zeitschritt deutlich zu erhöhen und so sehr rasche Konvergenz zum stationären Zustand zu erzielen. Dazu wird das Runge-Kutta-Verfahren umformuliert als

 
 
 

Dabei ist  .

Als Zyklus wird ein 4- oder 5-Gitter-W-Zyklus gewählt. Die Konvergenz wird dabei durch weitere Maßnahmen noch beschleunigt. Zunächst indem nicht in jeder Zelle mit demselben Zeitschritt gerechnet wird, sondern mit einer global gültigen CFL-Zahl, mittels derer dann der lokale Zeitschritt ausgerechnet wird (local time stepping). Allgemeiner haben Lee, Turkel, Roe und Van Leer nichtlineare Vorkonditionierer entwickelt, mit denen die Zeitableitung multipliziert wird.

Als zweites wird Residuenglättung genutzt (Jameson, Baker 1984), was den Glättungseffekt erhöht. Dazu werden in jeder Stufe des Runge-Kutta-Verfahrens die Flüsse   vor der Aufdatierung von   durch   ersetzt, wobei   durch die folgende Gleichung mit   klein gegeben ist:

 

Diese Gleichung lässt sich leicht lösen und erhält den stationären Zustand, ist aber für die Berechnung instationärer Probleme komplett ungeeignet.

Es ist möglich, mit diesem Verfahren stationäre Lösungen der Euler-Gleichungen für Tragflächenumströmungen mit drei bis fünf W-Zyklen zu berechnen (Jameson, Caughey 2001). Die stationären Navier-Stokes-Gleichungen sind schwieriger zu lösen, je nach Problem sind 30 bis 100 Schritte notwendig.

Dual Timestepping

Bearbeiten

Das Problem bei dem oben beschriebenen Mehrgitter-Verfahren ist, dass ein Großteil der Beschleunigungstechniken nur für stationäre Gleichungen funktioniert. Um die Technik auf instationäre Probleme zu übertragen wird das Dual Time Stepping Verfahren genutzt (Jameson 91). Gegeben sei dazu die Semidiskrete instationäre Gleichung

 

Wird diese mittels eines impliziten Verfahrens integriert, ergibt sich ein nichtlineares Gleichungssystem

 

wobei   von der konkreten Wahl der Zeitintegration abhängt. Wird dazu eine zusätzliche Zeitableitung addiert in Dualzeit, ergibt sich das modifizierte Problem

 

Die Lösung des eigentlichen Problems ist dann die Lösung des letzten Problems wenn die duale Zeitableitung verschwindet, also der stationäre Zustand in Dualzeit erreicht ist. Dafür wird das oben erwähnte Mehrgitterverfahren angewandt. Aufgrund der zusätzlichen Terme wird dabei das Spektrum des Operators   etwas verschoben, was die Konvergenz negativ beeinflusst. Im Fall der RANS-Gleichungen sind 50 bis 250 Iterationen notwendig.

Konvergenztheorie im hyperbolischen

Bearbeiten

Eine Konvergenztheorie im hyperbolischen Fall ist für keines der Verfahren bekannt. Ein wesentliches Problem sind Schocks. Durch die Kopplung auf dem groben Gittern werden Informationen von beiden Seiten des Schocks theoretisch nach Prolongation über das komplette Gitter verteilt. Dies erschwert die theoretische Analyse bisher entscheidend.

Fluid-Struktur-Interaktion

Bearbeiten

Als letztes betrachten wir ein Thema, bei dem hyperbolische Erhaltungsgleichungen nur eine Nebenrolle spielen, nämlich gekoppelte Systeme von partiellen Differentialgleichungen, bei denen eine die Strömung eines Fluids beschreibt und eine zweite eine Struktur. Das meistuntersuchte Problem ist dabei die Kombination, bei der Fluid und Struktur über Kräfte interagieren, handelt es sich um Luft, spricht man von Aeroelastizität. Das Fluid verformt also durch mechanische Kräfte das Fluid, was wiederum Kräfte in der Struktur bewirkt. Diese können sich in nichtlinearer Weise aufschaukeln. Beispiele sind

Eine andere Form ist die thermomechanische Kopplung, bei der an der Struktur aufgrund von Hitze plastische Verformungen zu befürchten sind oder hervorgerufen werden sollen. Beispiele sind

  • der Wiedereintritt von Raumfahrzeugen in die Atmosphäre, bei denen die erhitzte Atmosphäre das Gefährt erhitzt,
  • thermische Belastungen in Motoren und Triebwerken oder
  • Abkühlungsprozesse mittels Gasen bei Stahlumformung.

Modellierung

Bearbeiten

Gegeben sind zwei partielle Differentialgleichungen auf verschiedenen Gebieten   und  . Dazu seien Kopplungsbedingungen gegeben für die Menge  . Konkreter haben wir ein System für das Fluid

 

und eines für die Struktur

 

Die Variablensätze 1 und 2 sind dabei den verschiedenen Gebieten zugeordnet und die Kopplungsbedingung ist in den einzelnen Abbildungen enthalten. Zum Beispiel könnte f den Navier-Stokes- oder Euler-Gleichungen entsprechen und g der Wärmeleitungsgleichung in einer Struktur. Kopplungsbedingung wäre dann, dass auf dem gemeinsamen Rand die Temperatur T und der Wärmefluss q übereinstimmen. Dies ist ein erstes nützliches Modellproblem vor der echten thermomechanischen Kopplung.

Kopplungsverfahren

Bearbeiten

Grundsätzlich wird zwischen zwei Lösungsansätzen unterschieden: monolithischen, bei denen das gekoppelte System als ein Gleichungssystem gelöst wird, und partitionierten, bei denen bestehende Lösungsverfahren zur Lösung der Einzelsysteme eingesetzt werden. Wesentlicher Vorteil der letzteren Ansätze ist, dass bestehende Löser für die jeweiligen Probleme verwendet werden können, was den Entwicklungsaufwand drastisch reduziert. Für gekoppelte Probleme existiert dagegen in der Regel nichts, so dass die Implementierung eines monolithischen Verfahrens sehr hohen Aufwand bedeuten kann. Wir werden uns hier auf partitionierte Verfahren beschränken.

Die Einzelsysteme seien also mit dem jeweiligen Verfahren im Raum diskretisiert, wobei wir der Einfachheit halber die Notation für   und   beibehalten:

 
 
Lose Kopplung

Einfachster Fall ist nun die lose oder schwache Kopplung. Gegeben sei für die jeweiligen Systeme ein Zeitintegrationsverfahren  , was aus der Lösung   zum Zeitpunkt   die Lösung   zum Zeitpunkt   generiert. Die erste Möglichkeit ist dann folgender Algorithmus, der parallel ausgeführt werden kann:

 
 

Eine üblichere Variante, die nur sequentiell funktioniert, ist:

 
 

Beide Kopplungsverfahren sind maximal erster Ordnung in der Zeit und explixit, so dass sich erhebliche Zeitschritteinschränkungen ergeben können. Beide Varianten kann man innerhalb eines Zeitschritts iterieren, um die Genauigkeit zu erhöhen.

Einschub Thermische Kopplung

Wir betrachten nun konkret die Kopplung der Navier-Stokes-Gleichungen mit der Wärmeleitungsgleichung. Im diskreten Fall ist es nicht möglich, in Struktur und Fluid sowohl Temperatur als auch Wärmefluss am Rand vorzuschreiben, da dies für keine der Gleichungen sinnvolle Randbedingungen ergibt. Es gibt also die Möglichkeit, dem Fluid und der Struktur Temperatur oder Wärmefluss vorzuschreiben.

Die eine Möglichkeit bedeutet: Die Wärmeleitungsgleichung hat als Randbedingung eine Temperaturvorgabe und gibt dem Fluid einen Wärmefluss, während die Strömung als Randbedingung einen Wärmefluss erhält und eine Randtemperatur an die Struktur zurückgibt. Laut Giles 1997 ist dieses Verfahren instabil.

Deswegen erhält die Strömung am Rand Temperaturdaten und die Struktur einen Wärmefluss.

Starke Kopplung

Schwache Kopplung hat verschiedene Probleme: Geringe Ordnung, möglicherweise starke Zeitschrittbeschränkungen und schließlich erfüllt das Verfahren keine zusätzlichen Bedingungen wie Energieerhaltung. Eine Lösung ist eine implizite Herangehensweise, wobei dann die Frage ist, wie ein partitioniertes Lösungsverfahren aussieht. Folgende Formulierung bietet sich an (siehe beispielsweise Mathies, Steindorf 2002):

 
 

Ein Newton-Schritt für dieses System wäre dann:

 

Nun liegen die entsprechenden partiellen Ableitungen nicht vor. Diese können allerdings entweder durch eine matrixfreie Technik in einem Krylov-Unterraum-Verfahren wie oben bereitgestellt werden oder durch Newton-Schritte im lokalen Verfahren.

Literatur

Bearbeiten

Als Literatur ist zu empfehlen:

  • C. Hirsch: Numerical computation of internal and external flows, 2 Bände, 1988, Wiley & Sons, New York. Dieses Buch beschreibt umfassend den Stand der Wissenschaft zum damaligen Zeitpunkt.
  • Randall J. LeVeque: Numerical Methods for Conservation Laws: Lectures in Mathematics, 1992, Birkhäuser Verlag. Dies ist ein recht kurzes Buch, das aus einer Vorlesung LeVeques an der ETH Zürich entstanden ist und stellt die vermutlich beste Einführung in das Thema dar.
  • Randall J. LeVeque: Finite Volume Methods for Hyperbolic Problems, 2002, Cambridge Texts in Applied Mathematics, ISBN 0-521-00924-3. Zehn Jahre später das ganze nochmal, aber in ausführlich.
  • P. Wesseling: Principles of Computational Fluid Dynamics, 2000, Springer Verlag. Wesseling geht sehr gewissenhaft vor und behandelt Strömungsmechanik und Verfahren für inkompressible und kompressible Strömungen.

Zusammenfassung des Projekts

Bearbeiten

  „Theorie und Numerik von Erhaltungsgleichungen“ ist nach Einschätzung seiner Autoren zu 90 % fertig

  • Zielgruppe: Studierende der Mathematik, Physik oder des Maschinenbaus ab dem 5. Semester mit grundlegenden Kenntnissen in Numerik.
  • Lernziele: Grundlegende theoretische Aussagen zu Erhaltungungsgleichungen und der zugehörigen Numerik mit Schwerpunkt nichtlinearer Gleichungen der Strömungsmechanik. Darüberhinaus grundlegendes Wissen zu Finite-Volumen-Verfahren und nichtlinearen Lösungsverfahren.
  • Buchpatenschaft / Ansprechperson: Benutzer:P. Birken
  • Sind Co-Autoren gegenwärtig erwünscht? Nein, aber die Korrektur von Tippfehlern oder Anmerkungen zu Fehlern auf der Diskussionsseite ist erwünscht.

Diese Seite hat mittlerweile eine Größe erreicht, die es als geeignet erscheinen lässt, sie in mehrere einzelne Seiten zu zerlegen. In welche Teile diese Seite zerlegt werden könnte, kann auf der Diskussionsseite besprochen werden. Wie man es macht, steht im Wikibooks-Lehrbuch im Abschnitt Buch in Kapitel untergliedern.