close

Вход

Забыли?

вход по аккаунту

?

4156.Numerische Methoden der Stroemungsmechanik 001 .pdf

код для вставкиСкачать
1
Numerische Methoden der Strömungsmechanik
Dr.-Ing. habil. Dipl.-Phys. Andreas Malcherek1
Vorlesungsskript, Version 5.6
1
Bundesanstalt für Wasserbau, Außenstelle Küste, Wedeler Landstr. 157, 22559 Hamburg, Tel.: 040 / 81908346, email: malcherek@hamburg.baw.de
2
Vorwort
Die Anfänge dieses Skriptes gehen zurück auf meine Teilnahme an einem NATO-Workshop
im Jahr 1993 über ’Computer Modeling of Free-Surface and Pressurized Flows’, dessen Ergebnisse in einem gleichnamigen Buch veröffentlicht wurden. Hier gab mir Prof. Zielke die
Möglichkeit, einen Beitrag über numerische Verfahren zur Advektion zu machen. Im Rahmen dieser Arbeit lernte ich vor allem, daß es unmöglich ist, eine beliebige Funktion exakt
auf einem diskreten Gitter zu verschieben. Dies bedeutet nichts anderes, als daß es nie und
nimmer möglich ist, Strömungen in einem Computermodell naturgleich zu modellieren. Diese
Erkenntnis nahm mir mein Mißtrauen vor dieser Technik und begeisterte mich für die Aufgabe, Natur im Computermodell so gut wie möglich zu simulieren, ohne dabei die Ehrfurcht vor
derselben zu verlieren zu müssen.
Ein Lehrauftrag über ’Numerische Methoden für Strömungen, Stoff- und Wärmetransport’ im Wintersemester 1995/96 am Institut für Strömungsmechanik und elektronisches
Rechnen im Bauwesen der Universität Hannover ließ die Versionen 1.0ff entstehen. Da ich in
diesem Semester außerdem noch eine Vorlesungsreihe über ’Turbulenz- und Stofftransportmodelle für Fließgewässer’ am Institut für Geodynamik der Universität Bonn halten durfte,
bestand gleich zweierlei Anlaß neben der reinen Numerik auch das darzustellen, was zu modellieren ist. Dabei hatte ich Physik und Numerik aus didaktischen Gründen vermischt, so daß
weder das eine noch das andere langweilig werden sollte, gleichzeitig aber der Schwierigkeitsgrad kontinuierlich steigt. Version 1 des Skriptes wuchs im Laufe dreier Vorlesungszyklen auf
über 230 Seiten.
Dieses Skript wurde in der Version 2 im wesentlichen in seinem hydrodynamischen Teil zur
Buchform erweitert. Mit der Version 3 wurde die 400-Seiten-Schwelle überschritten und die
’Hydrodynamik der Fließgewässer’ von der hier vorliegenden Schrift ’Numerischen Methoden der Hydrodynamik’ vollständig inhaltlich aber nicht örtlich getrennt. Version 4 diente der
Verbesserung des hydrodynamischen Teils und reifte schließlich zur Habilitationsschrift.
Mit der Version 5 liegen die numerischen Methoden nun als eigene Schrift vor. Damit sollte den Studenten meines wieder aufgenommenen Lehrauftrages in Hannover zunächst einmal
ein eigenständiges Skript geboten werden, welches den Stoff der Hydromechanik auf das wesentliche beschränkt. Ferner habe ich mit diesem Lehrauftrag zusammen mit den Studenten
Java-Applikationen entwickelt, mit denen die Numerik spielerisch erprobt werden kann.
Die vorliegende Version 5.6 ist durch Anregungen beim Klausurseminar des Graduiertekollegs 615 der Universität Hannover entstanden. Das Galerkinverfahren wird nun am einfacheren Beispiel der Poissongleichung erklärt, zeitabhängigen FE-Methoden ist ein eigenständiges
Kapitel gewidmet.
Inhaltsverzeichnis
1 Einführung
9
2 Zeitschrittverfahren
2.1 Die Diskretisierung der Zeit . . . . . . . . . . .
2.2 Die allgemeine Evolutionsgleichung . . . . . . .
2.3 Einschrittverfahren . . . . . . . . . . . . . . . .
2.3.1 Das Eulerverfahren . . . . . . . . . . . .
2.3.2 Das implizite Verfahren . . . . . . . . .
2.3.3 Das Crank-Nicolson-Verfahren . . . . . .
2.3.4 Taylorverfahren . . . . . . . . . . . . . .
2.3.5 Runge-Kutta-Verfahren . . . . . . . . . .
2.4 Konsistenz . . . . . . . . . . . . . . . . . . . . .
2.5 Mehrschrittverfahren . . . . . . . . . . . . . . .
2.5.1 Numerische Integration . . . . . . . . . .
2.5.2 Leap-Frog-Verfahren . . . . . . . . . . .
2.6 Zeitliche Mittlung der Grundgleichungen . . . .
2.7 Das Entity-Relationship-Modell der Zeit . . . . .
2.7.1 Das Entity-Relationship-Modell . . . . .
2.7.2 Das ER-Modell der Zeit . . . . . . . . .
2.8 Die logistische Differentialgleichung als Beispiel
2.9 Zusammenfassung . . . . . . . . . . . . . . . .
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
13
14
14
15
16
16
16
17
17
18
19
19
20
20
21
21
22
23
27
3 Stabilität
3.1 Stetige Abhängigkeit von den Eingangsdaten . . . . . . . .
3.1.1 Stetige Abhängigkeit beim kontinuierlichen Problem
3.1.2 Stetige Abhängigkeit beim diskretisierten Problem .
3.2 Der Verfahrensoperator . . . . . . . . . . . . . . . . . . . .
3.2.1 Matrixnormen . . . . . . . . . . . . . . . . . . . .
3.2.2 Spektralradius und Stabilität . . . . . . . . . . . . .
3.3 Allgemeine Stabilitätskriterien . . . . . . . . . . . . . . . .
3.3.1 Stabilitätskriterien für explizite Einschrittverfahren .
3.3.2 Stabilitätskriterien für implizite Einschrittverfahren .
3.4 Iterative und Prädiktor-Korrektor-Verfahren . . . . . . . . .
3.5 Operator-Splitting-Verfahren . . . . . . . . . . . . . . . . .
3.6 Konvergenz . . . . . . . . . . . . . . . . . . . . . . . . . .
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
29
30
30
32
33
34
35
35
36
36
37
37
39
3
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
INHALTSVERZEICHNIS
4
4 Die Methode der Finiten Differenzen
4.1 Einfache Differenzenquotienten . . . . . . . . . . .
4.1.1 Differenzenverfahren in Matrixschreibweise .
4.1.2 Die stationäre Transportgleichung . . . . . .
4.2 Konsistenz . . . . . . . . . . . . . . . . . . . . . . .
4.3 Stabilität . . . . . . . . . . . . . . . . . . . . . . . .
4.4 Verallgemeinerungen . . . . . . . . . . . . . . . . .
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
5 Die Advektionsgleichung
5.1 Theorie der Advektionsgleichung . . . . . . . . . . . . . . . . .
5.1.1 Die Advektionsgleichung als Evolutionsgleichung . . .
5.1.2 Differentialgleichungssysteme 1. Ordnung . . . . . . .
5.1.3 Bahnlinien als Charakteristiken der Advektionsgleichung
5.1.4 Das Anfangswertproblem . . . . . . . . . . . . . . . .
5.1.5 Das Randanfangswertproblem . . . . . . . . . . . . . .
5.2 FD-Verfahren für die Advektionsgleichung . . . . . . . . . . .
5.2.1 Explizite Verfahren . . . . . . . . . . . . . . . . . . . .
5.2.2 Implizite Verfahren . . . . . . . . . . . . . . . . . . . .
5.2.3 Prädiktor-Korrektor-Verfahren . . . . . . . . . . . . . .
5.3 Numerische Diffusion . . . . . . . . . . . . . . . . . . . . . . .
5.4 Numerische Dispersion . . . . . . . . . . . . . . . . . . . . . .
6 Interpolation auf Finiten Elementen
6.1 Lagrangesche Interpolationspolynome . . . . . .
6.2 Finite Elemente . . . . . . . . . . . . . . . . . .
6.3 Eindimensionale Elemente . . . . . . . . . . . .
6.3.1 Lineare Interpolation . . . . . . . . . . .
6.3.2 Der Interpolationsfehler . . . . . . . . .
6.3.3 Quadratische Interpolation . . . . . . . .
6.3.4 Positivität und andere Extremalprinzipien
6.4 Zweidimensionale Elemente . . . . . . . . . . .
6.4.1 Dreieck mit linearem Ansatz . . . . . . .
6.4.2 Viereck mit bilinearem Ansatz . . . . . .
6.4.3 Viereck mit biquadratischem Ansatz . . .
6.4.4 Sechseck mit quadratischem Ansatz . . .
6.5 Interpolation im Dreidimensionalen . . . . . . .
6.5.1 Ein Prisma mit linearem Ansatz . . . . .
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
7 Die Diffusionsgleichung
7.1 Die Diffusionsgleichung als Evolutionsgleichung . . . . . . .
7.2 Parabolische semilineare Differentialgleichungen . . . . . . .
7.3 Das Anfangswertproblem (1D) . . . . . . . . . . . . . . . . .
7.4 FD-Verfahren für die Diffusionsgleichung . . . . . . . . . . .
7.4.1 Das Forward-Time-Centered-Space-Verfahren (FTCS)
7.4.2 Das Randanfangswertproblem . . . . . . . . . . . . .
7.4.3 Richardson- und DuFort-Frankel-Verfahren . . . . . .
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
41
41
43
44
45
45
48
.
.
.
.
.
.
.
.
.
.
.
.
51
51
51
52
52
52
53
53
54
58
60
60
61
.
.
.
.
.
.
.
.
.
.
.
.
.
.
63
63
64
65
65
66
68
70
70
70
73
73
74
76
76
.
.
.
.
.
.
.
79
79
79
82
83
83
85
86
INHALTSVERZEICHNIS
7.5
5
7.4.4 Das implizite Verfahren . . . . . . . . . . . . . . . . . . . . . . . .
7.4.5 Das Crank-Nicolson-Verfahren . . . . . . . . . . . . . . . . . . . . .
Bewertung . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
8 Lagrange-Verfahren
8.1 Das Lagrange-Verfahren für die Advektionsgleichung
8.2 Numerische Dämpfung . . . . . . . . . . . . . . . .
8.3 Monotonie . . . . . . . . . . . . . . . . . . . . . . .
8.4 Beispiele . . . . . . . . . . . . . . . . . . . . . . .
8.4.1 Stufenfunktion . . . . . . . . . . . . . . . .
8.4.2 Hunte . . . . . . . . . . . . . . . . . . . . .
8.4.3 Jade . . . . . . . . . . . . . . . . . . . . . .
9 Die Transportgleichung
9.1 Das Anfangswertproblem (1D) . . . . . . . . . . . .
9.2 Das Randanfangswertproblem . . . . . . . . . . . .
9.3 FD-Verfahren für die Transportgleichung . . . . . .
9.3.1 Die stossartige und punktförmige Einleitung
9.3.2 Explizite Verfahren . . . . . . . . . . . . . .
9.3.3 Implizite Verfahren . . . . . . . . . . . . . .
9.4 Lagrange-Verfahren für die Transportgleichung . . .
9.5 Bewertung . . . . . . . . . . . . . . . . . . . . . . .
10 Die Burgersgleichung
10.1 Der Ansatz von E. Hopf . . . . . . . . . . . .
10.2 Der Ansatz von J.D. Cole . . . . . . . . . . .
10.3 FD-Verfahren für die Burgersgleichung . . .
10.4 Lagrange-Verfahren für die Burgersgleichung
10.5 Harmonische Anfangsbedingungen . . . . . .
10.6 Zusammenfassung . . . . . . . . . . . . . .
.
.
.
.
.
.
.
.
.
.
.
.
11 Finite-Volumen-Verfahren
11.1 Die Divergenzform der Grundgleichungen . . . .
11.2 Diskretisierung . . . . . . . . . . . . . . . . . .
11.2.1 Die Behandlung der Zeitableitung . . . .
11.2.2 Die Behandlung des Quellterms . . . . .
11.2.3 Die Behandlung des Flußterms . . . . . .
11.3 Die Diffusionsgleichung auf äquidistantem Gitter
11.4 Die Diskretisierung advektiver Terme . . . . . .
11.4.1 Zentrale Verfahren . . . . . . . . . . . .
11.4.2 Upstream-Verfahren . . . . . . . . . . .
11.4.3 QUICK- und QUICKEST-Verfahren . . .
11.5 Globale und lokale Massenerhaltung . . . . . . .
11.6 Gestaffelte Gitter . . . . . . . . . . . . . . . . .
11.7 Finite Volumen auf Dreiecksnetzen . . . . . . . .
11.7.1 Der äußere Normalenvektor im Dreieck .
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
86
87
87
.
.
.
.
.
.
.
89
90
91
93
95
95
95
95
.
.
.
.
.
.
.
.
99
99
100
100
100
101
103
107
109
.
.
.
.
.
.
111
112
113
114
116
116
119
.
.
.
.
.
.
.
.
.
.
.
.
.
.
121
121
121
125
125
126
126
127
127
128
128
130
130
132
135
INHALTSVERZEICHNIS
6
11.7.2 Upstreaming auf Dreiecksnetzen . . . . . . . . . . . . . . . . . . . . 135
11.8 Finite Differenzen und Finite Volumen . . . . . . . . . . . . . . . . . . . . . 135
12 Die Saint-Venantschen Gleichungen
12.1 Hyperbolizität der Saint-Venant-Gleichungen . . . . . . . . . . . .
12.2 Charakteristiken der 1D-Gleichungen . . . . . . . . . . . . . . . .
12.2.1 Das inverse Charakteristikenverfahren . . . . . . . . . . . .
12.2.2 Randbedingungen . . . . . . . . . . . . . . . . . . . . . .
12.3 Charakteristiken der 2D-Gleichungen . . . . . . . . . . . . . . . .
12.3.1 Die stationären tiefengemittelten Gleichungen . . . . . . . .
12.3.2 Inverse Bicharakteristikenverfahren für die 2D-Gleichungen
12.4 Das semiimplizite Verfahren von Casulli . . . . . . . . . . . . . . .
12.5 Bewertung . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
137
139
139
140
142
144
145
147
148
152
13 Galerkinverfahren
13.1 Das Standard-Galerkinverfahren . . . . . .
13.2 Die Poissongleichung in der Hydrodynamik
13.2.1 Potentialströmungen . . . . . . . .
13.2.2 Druck-Korrektur-Verfahren . . . . .
13.3 Lösung der 1D-Poissongleichung . . . . . .
13.4 Integrationsverfahren . . . . . . . . . . . .
13.4.1 Die Integraltransformationsformel .
13.4.2 Analytische Integration . . . . . . .
13.4.3 Numerische Integration . . . . . . .
13.5 Die Poissongleichung auf Dreiecken . . . .
13.6 Zusammenfassung . . . . . . . . . . . . .
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
153
153
154
154
155
156
158
158
159
159
160
164
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
14 Zeitabhängige Finite-Elemente-Methoden
14.1 Lösung der 1D-Transportgleichung . . . . . . . . . . . . . . . . . .
14.2 Die Diskretisierung der Zeit durch Finite Elemente . . . . . . . . .
14.3 Schwache Lagrangeverfahren . . . . . . . . . . . . . . . . . . . . .
14.4 Upstream-Strategien: Petrov-Galerkin-Verfahren . . . . . . . . . .
14.4.1 Das Verfahren nach Christie . . . . . . . . . . . . . . . . .
14.4.2 Das Streamline Upwind/Petrov-Galerkin (SU/PG) Verfahren
14.5 Zusammenfassung . . . . . . . . . . . . . . . . . . . . . . . . . .
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
165
165
168
169
171
171
172
173
15 Finite Elemente und Funktionalanalysis
15.1 Funktionenräume . . . . . . . . . . . . . . . . . .
15.2 Hilberträume . . . . . . . . . . . . . . . . . . . .
15.2.1 Der Hilbertraum
. . . . . . . . . . .
15.2.2 Schwache Ableitungen und Sobolevräume .
15.2.3 Operatoren in Hilberträumen . . . . . . . .
15.2.4 Eigenwerte linearer Operatoren . . . . . .
15.3 Dualräume . . . . . . . . . . . . . . . . . . . . . .
15.4 Bilinearformen . . . . . . . . . . . . . . . . . . .
15.5 Variationsaufgaben . . . . . . . . . . . . . . . . .
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
175
175
176
177
178
179
180
181
181
184
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
INHALTSVERZEICHNIS
15.5.1 Die Poissongleichung . . . . . .
15.5.2 Die Helmholtzgleichung . . . . .
15.5.3 Die stationäre Transportgleichung
15.6 Gemischte Finite Elemente . . . . . . . .
15.6.1 Das Stokes-Problem . . . . . . .
15.6.2 Die Babuska-Brezzi-Bedingung .
15.7 Konvergenz . . . . . . . . . . . . . . . .
7
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
184
184
185
186
186
187
187
16 Gleichungslöser
189
16.1 Diskretisierung und Systemmatrix . . . . . . . . . . . . . . . . . . . . . . . 189
16.2 Direkte Gleichungslöser . . . . . . . . . . . . . . . . . . . . . . . . . . . . 189
16.3 Iterative Gleichungslöser . . . . . . . . . . . . . . . . . . . . . . . . . . . . 190
8
INHALTSVERZEICHNIS
Kapitel 1
Einführung
Diese Schrift beschäftigt sich mit dem Repertoire numerischer Methoden, die man für den
Entwurf eines hydrodynamisch-numerischen Modells für Oberflächengewässer verwenden
kann. Bevor wir uns diesen im einzelnen zuwenden, fragen wir uns, was man unter einem
hydrodynamisch-numerischen Modell - in Kurzform auch HN-Modell - versteht.
Ein hydrodynamisch-numerisches Modell ist zunächst erst einmal ein numerisches Modell.
Unter einem solchen versteht man ein Verfahren, welches in der Regel große Mengen von
Zahlen verarbeitet. Der Prozess der Verarbeitung von Zahlen kann durch Flußdiagramme
dargestellt werden, welche in einer endlichen Abfolge von Schritten ein Ergebnis liefern. Solche Flußdiagramme sollten immer eine Klasse von Problemen repräsentieren, von denen das
Einzelproblem durch eine bestimmte Konfiguration der Eingangsdaten bestimmt wird. Zu jedem speziellen Problem sollte das Flußdiagramm eine eindeutige Menge von Ausgangsdaten
liefern. Hat dieses das Flußdiagramm die genannten Eigenschaften, spricht man von einem
Algorithmus.
Ein hydrodynamisch-numerisches Modell sollte die guten Eigenschaften eines Algorithmus
besitzen. Die von ihm geliefert numerische Lösung besteht wieder aus einer sehr großen Menge von Zahlen bzw. Daten, deren Struktur der der hydrodynamischen Probleme angepaßt ist.
Hydrodynamische Probleme beziehen sich immer auf ein (natürliches oder künstliches) Modellgebiet in Raum und Zeit, also etwa auf die Elbe zwischen Stromkilometer 400 und 500
in der Zeit vom 1. bis zum 30 Juni 1990. Dazu wird das Modellgebiet mit einem Netz von
Punkten überdeckt, die man als Knoten bezeichnet. An jedem dieser Knoten liefert das numerische Modell Ergebnisse für die gesuchten hydrodynamischen Größen. Das Raumgitter kann
dabei eine drei-, zwei- oder eindimensionale Struktur aufweisen. Bei zweidimensionalen Gittern wird dabei die über die Tiefe, bei eindimensionalen Gittern die über den Fließquerschnitt
integrierte Strömung modelliert.
Der Simulationszeitraum wird ebenfalls durch einzelne Zeitpunkte diskretisiert, zu denen die
Simulationsergebnisse gewonnen werden. So entstehen insgesamt vier-, drei- oder zweidimensionale Datenstrukturen als Ergebnis eines hydrodynamisch-numerischen Modells. Abbildung
1.2 veranschaulicht für ein zweidimensionales Raumgitter die entstehende dreidimensionale
Datenstruktur.
Wir werden uns die in einem HN-Modell benötigten numerischen Methoden konstruktiv d.h.
Schritt für Schritt erarbeiten. Dazu benötigt der Leser kaum Kenntnisse aus dem Bereich der
9
KAPITEL 1. EINFÜHRUNG
10
Eingangsdaten
Algorithmus
Flußdiagramm
Allgemeinheit
Eindeutigkeit
Endlichkeit
Ausgangsdaten
Abbildung 1.1: Die Struktur eines numerischen Modells
t = 2 D t
t = D t
t = 0
ie t
b
e
sg
nu g
ö m
S tr
Abbildung 1.2: Dreidimensionale Datenstruktur für eine tiefenintegrierte Strömung
11
Hydrodynamik der Oberflächengewässer. Methodisch funktioniert das folgendermaßen: Wir
werden uns einfache Differentialgleichungen vornehmen und studieren, was diese simulieren
können. Dann werden wir Verfahren zur numerischen Behandlung derselben konstruieren. Im
Laufe der Schrift werden diese Modelle aus Differentialgleichungen und zugehörigen numerischen Verfahren sich der tatsächlichen Hydromechanik immer weiter annäherten aber auch
komplexer werden.
Unser Weg beginnt bei den gewöhnlichen Differentialgleichungen erster Ordnung. Mit denen kann man dynamische d.h. zeitabhängige aber räumlich homogene Prozesse beschreiben.
Durch die Advektionsgleichung bekommen wir im folgenden ein Werkzeug in die Hand, die
Bewegung von Stoffen oder Eigenschaften mit einer Strömung zu beschreiben. Sie ist von ihrer Struktur die einfachste partielle Differentialgleichung, wir werden aber an ihr schon die
Grenzen der numerischen Verfahren leidvoll studieren können.
Dann werden wir uns dem Prozess der Diffusion zuwenden. Koppeln wir diesen mit der Advektion, dann sind wir schon in der Lage, die verschiedensten Transportvorgänge in Fließgewässern naturnah zu simulieren. So transportiert die Strömung aber auch ihren eigenen Impuls und damit irgendwie sich selbst, wodurch eine Rückkopplung entsteht. Dieses Phänomen
kann man im einfachsten Fall durch die Burgersgleichung beschreiben, die uns in das faszinierende Gebiet der nichtlinearen Dynamik und ihre numerische Simulation einführen wird.
Mit den Saint-Venant-Gleichungen treffen wir dann auf ein System von partiellen Differentialgleichungen, wodurch weitere Anforderungen an die Qualität der numerischen Verfahren
gestellt werden. Das Modell ist aber nun schon in der Lage, tiefengemittelte Strömungen in
Fließgewässern zu simulieren. Letztlich wenden wir uns den Reynolds- und Navier-StokesGleichungen zu, die dann die vollständige Hydrodynamik der Fließgewässer beschreiben.
12
KAPITEL 1. EINFÜHRUNG
Kapitel 2
Zeitschrittverfahren
Die Ergründung des Wesens der Zeit ist eines der faszinierendsten Abenteuer des menschlichen Geistes.
Zunächst sind da die Fragen nach ihrem Anfang, präziser, ob sie überhaupt einen solchen
hatte und wenn, wann er war. Der Kirchenlehrer Augustinus (354 - 430) formulierte diese
Frage theologisch1 : Was machte Gott, bevor er Himmel und Erde machte ? Ich gebe nicht das
zur Antwort, was jemand gesagt haben soll, der mit einem Witzwort der schwierigen Frage
ausweichen wollte und sagte: ’Er baute eine Hölle für die Leute, die zu hohe Dinge erforschen
wollen.’ Das will ich nicht zur Antwort geben, Einsehen und Lachen sind zwei verschiedene
Dinge. ... Ich hingegen sage, du, unser Gott, seist der Sch öpfer aller Kreatur, und wenn man
unter dem Wort ’Himmel und Erde’ alle Kreatur versteht, sage ich mit Zuversicht: Bevor Gott
Himmel und Erde machte, machte er nichts. Denn wenn er etwas machte, was machte er,
wenn nicht ein Geschöpf ? Augustinus fragt sich weiter, warum Gott vor der Schöpfung eine
unbestimmt lange Zeit nichts tat, denn wenn er ’Himmel und Erde’ schon im Sinn hatte, hätte
er ja auch früher anfangen können. Dieser absurde Gedanke läßt nur einen Schluß zu: Gott
hatte zunächst bzw. mit der Schöpfung die Zeit erschaffen. Denn eben diese Zeit hattest du
gemacht, und es konnte keine Zeit vorübergehen, bevor du die Zeiten gemacht hattest. Die
moderne Physik sieht das nicht anders, auch sie geht davon aus, daß es vor der Entstehung des
Universums Zeit so nicht gegegeben, bzw. ihr Begriff keine Gültigkeit hat. In vollkommener
Analogie zu Augustinus gibt es keinen Grund, warum der Urknall nicht auch schon hätte früher
stattfinden können, wenn es so etwas wie die Zeit schon vorher gegeben hätte.
Nachdem er dieses so eindeutig geklärt hat, macht sich Augustinus daran, den Ablauf der
Zeit zu ergründen. Er stellt zunächt einmal fest, daß es Vergangenheit und Zukunft nicht gibt,
Vergangenes existiert nur in der gegenwärtigen Erinnerung und die Zukunft nur in der gegenwärtigen Erwartung: Denn meine Kindheit, die nicht mehr ist, liegt in einer vergangenen
Zeit, die nicht mehr ist; aber ihr Bild, das ich heraufhole, wenn ich von ihr erz ähle, sehe ich im
gegenwärtigen Augenblick, weil es noch in meinem Gedächtnis ist. ... Sagt man vom Zuk ünftigen, es werde gesehen, dann wird nicht dieses Zukünftige selbst gesehen, das noch nicht ist,
sondern wohl dessen Ursachen und Zeichen, die schon sind. ... Aus ihnen erfaßt der Geist das
Zukünftige und kann es dann vorhersagen. ... Es gibt daf ür unzählige Beispiele; ich nehme
nur eines davon: Ich sehe die Morgenröte und sage voraus, die Sonne werde aufgehen. Was
ich sehe ist gegenwärtig; was ich voraussage, ist zuk ünftig. Somit existiert allerhöchstens die
1
Die folgenden Zitate stammen aus Augustinus, Bekenntnisse, Reclam Universal-Bibliothek Nr. 2792, Stuttgart 1989.
13
KAPITEL 2. ZEITSCHRITTVERFAHREN
14
Gegenwart. Doch wie lang ist diese ? Augustinus seziert den gegenwärtigen Tag und findet
darin vergangene und zukünftige Stunden. Zwar geht er iterativ zu kleineren Zeitintervallen,
kommt aber nicht zu einer Grenzwertbildung gegen Null, wodurch auch die Gegenwart als
existent bezweifelt würde. Er läßt sich nicht beirren und bleibt dabei, daß Zeit aus gewissen
Zeitspannen besteht. Tatsächlich schreitet die Zeit in Zeitsprüngen bzw. Zeitquanten voran,
deren Länge man als Plancksche Fundamentalzeit bezeichnet. Sie ist
wobei die universelle Gravitationskonstante, das Plancksche Wirkungsquantum und die
Lichtgeschwindigkeit ist. Da eine Sekunde somit aus sehr vielen Zeitquanten besteht, entsteht
der Eindruck, als fließe die Zeit kontinuierlich dahin.
Für unser weiteres Modellieren des Phänomens Zeit ist die Tatsache entscheidend, daß sich die
Zeit grundsätzlich nicht in Zeitpunkte, sondern nur in Zeitintervalle bzw. Zeitschritte zerlegen
läßt, die aber praktischerweise nicht so kurz wie die Plancksche Fundamentalzeit sein müssen.
2.1 Die Diskretisierung der Zeit
In den Grundgleichungen der Hydrodynamik tauchen zeitliche Ableitungen auf, diese Gleichungen beschreiben also, wie sich entsprechende Größen mit der Zeit entwickleln. Um dieses Geschehen numerisch zu behandeln, gehen wir genau so vor, wie es die Ontologie der Zeit
verlangt: Wir legen einen Anfangszeitpunkt fest und simulieren den Lauf der Zeit durch das
fortwährende Hinzufügen von Zeitschritten , wodurch eine Folge von Zeitpunkten entsteht:
(2.1)
Dies ist die Diskretisierung der Zeit.
Wir gehen davon aus, daß zur Zeit alle physikalischen Größen in Form von Anfangsbedingungen bekannt sind. Gesucht sind Rechenvorschriften, die die Lösungswerte zum Zeitschritt
und dann sukzessive für alle weiteren Zeitpunkte , ,..., ergeben.
Nach der Diskretisierung der Zeit wollen wir nun den Typ der Differentialgleichung spezifizieren, den wir in diesem Kapitel numerisch untersuchen wollen. An diesen stellen wir die
Anforderungen, daß er möglichst allgemein ist, damit wir alle zeitabhängigen Differentialgleichungen der Hydrodynamik erfassen und daß er möglichst einfach ist, damit wir es hier am
Anfang noch nicht zu schwer haben. Diese Bedingungen erfüllt die allgemeine Evolutionsgleichung.
2.2 Die allgemeine Evolutionsgleichung
Wir betrachten eine zeitabhängige Differentialgleichung der Form
(2.2)
die auch als allgemeine Evolutionsgleichung bezeichnet wird. In dieser Darstellung enthält der
Operator nur noch Ortsableitungen, so ergibt sich die x-Navier-Stokes-Gleichung durch
2.3. EINSCHRITTVERFAHREN
15
und
Später werden wir die Ortsableitungen ebenfalls diskretisieren; wir werden sehen, daß der
Operator dann zu einer Matrix degeneriert.
Eine formale Lösung der homogenen ( ) zeitabhängigen Differentialgleichung (2.2) ist
durch
(2.3)
gegeben. Formal ist die Lösung deshalb, weil wir zuerst klären müssen, wie man eine Potenz
mit einem Operator als Exponenten behandelt. Naheliegend ist hier die Bildung der Taylorreihe:
(2.4)
Um alle Formalitäten zu erledigen, bliebe noch die Konvergenzfrage für die Reihe zu untersuchen.
Auch die Lösung der inhomogenen Evolutionsgleichung kann formal angegeben werden; sie
lautet:
¼
Schließlich wollen wir noch eine ähnliche Abschätzung für nichlineare Gleichungen vom Typ
der Navier-Stokes-Gleichungen gewinnen. Dazu schreiben wir die Evolutionsgleichung in der
Form [?]
wobei der Operator linear und nichtlinear ist. Indem die Nichtlinearität als Inhomogenität
behandelt wird, ergibt sich die formale Lösung
¼
¼
2.3 Einschrittverfahren
Zur Bestimmung der Lösung an einem neuen Zeitpunkt wird die Zeitableitung in der
Form
KAPITEL 2. ZEITSCHRITTVERFAHREN
16
approximiert werden.
Die zeitlich semidiskretisierte Gleichung wird zu
und es verbleibt die Frage, zu welchem Zeitpunkt wir und bei inhomogenen Gleichungen
auswerten sollen.
2.3.1 Das Eulerverfahren
Bei diesem Verfahren wird zum bekannten Zeitschritt ausgewertet:
(2.5)
Eine Differentialgleichung ist nun nicht mehr zu lösen; das Problem hat sich auf die Auswertung eines reinen Differentialausdruckes reduziert. Das Eulerverfahren benötigt somit einen
sehr geringen Rechenaufwand pro Zeitschritt, allerdings hat es bei zu hohen Zeitschritten Probleme mit der Stabilität, worauf später noch eingegangen wird.
Wir werden später sehen, daß der Fehler des Eulerverfahrens proportional zum Zeitschritt gegen Null geht.
2.3.2 Das implizite Verfahren
Bei impliziten Verfahren wird zum unbekannten Zeitschritt ausgewertet:
(2.6)
Da nun Differentialausdrücke der unbekannten Funktion verbleiben, ist weiterhin eine
partielle Differentialgleichung zu lösen.
Das implizite Verfahren ist somit wesentlich rechenintensiver als das Eulerverfahren. Sein Vorteil liegt jedoch in den wesentlich besseren Stabilitätseigenschaften, wodurch erheblich größere Zeitschritte gewählt werden können und der Nachteil der Rechenintensität somit wieder
kompensiert wird.
2.3.3 Das Crank-Nicolson-Verfahren
Das Crank-Nicolson-Verfahren besteht aus einer Wichtung des bekannten Zeitschrittes und
des unbekannten Zeitschrittes mit einem Wichtungsparameter :
(2.7)
Genauso wie beim impliziten Verfahren verbleibt weiterhin die Lösung einer partiellen Differentialgleichung zum Zeitschritt . Für geht der Verfahrensfehler proportional zu
2.3. EINSCHRITTVERFAHREN
17
zu Null. Ein kleinerer Zeitschritt führt im Vergleich zum Eulerverfahren hier zu günstigeren Ergebnissen.
Das Crank-Nicolson-Verfahren bietet die Möglichkeit kontinuierlich zwischen den Zeitschritten zu wichten. Es wird daher in den meisten Codes verwendet, wenn nicht das implizite
Verfahren aus Stabilitätsgründen zwingend erforderlich ist.
In Verallgemeinerung dessen lassen sich sogar verschiedene Crank-Nicolson-Faktoren für verschiedene Terme der Differentialgleichungen einführen, stellt diese sich z.B. dar als
so wäre das entsprechende Crank-Nicolson-Verfahren durch
gegeben.
2.3.4 Taylorverfahren
Die formale Lösung der homogenen Evolutionsgleichung kann auch zur Konstruktion numerischer Verfahren verwendet werden. Dazu schreiben wir die formale Lösung für den Zeitschritt
als
und erhalten bei der Berücksichtigung der ersten beiden Gleider der Taylorreihe das Eulerverfahren. Nimmt man noch den dritten Term hinzu
so ergibt sich das Verfahren
(2.8)
welches auch Lax-Wendroff-Verfahren [17] genannt wird.
Vollkommen analog zu dieser Herleitung ist die Berücksichtigung des lokalen Verfahrensfehlers mit umgekehrtem Vorzeichen in der zu diskretisierenden Ausgangsgleichung. Diese
Vorgehensweise wird als Defektapproximation ([9]) bezeichnet.
2.3.5 Runge-Kutta-Verfahren
Entwickelt man die formale Lösung bis zur vierten Ordnung, so entsteht das Verfahren
welches man in der sehr einfachen Form
(2.9)
KAPITEL 2. ZEITSCHRITTVERFAHREN
18
mit
programmieren kann. Das Verfahren wird als Runge-Kutta-Verfahren bezeichnet.
2.4 Konsistenz
Die exakte Lösung der Differentialgleichung ist i.a. keine Lösung der diskreten Gleichung.
Setzt man die exakte Lösung dennoch in letztere ein, so erhält man ein Maß für das Abweichen
des diskreten Verfahrens von der zu lösenden Differentialgleichung, welches man als lokalen
Verfahrensfehler bezeichnet.
Für das allgemeine Einschrittverfahren setzt man den lokalen Verfahrensfehler in der
Form
an, wobei mit die exakte Lösung bezeichnet wird. Offensichtlich wird der Fehler genau dann
Null, wenn die exakte Lösung auch die diskretisierte Gleichung erfüllt.
Wird durch die Taylorreihe
dargestellt, so erhält der lokale Diskretisierungsfehler nach kurzer Rechnung die Form
und da und da 2.5. MEHRSCHRITTVERFAHREN
19
wobei die Linearität von vorausgesetzt werden muß. Betrachtet man nur den führenden Term
der Reihe, so
Das Verhalten der führenden Fehlerterme gibt Anlaß zu folgender Definition:
Def: Die Konsistenzordnung eines Verfahrens ist die jeweils minimale Potenz, mit denen
die Schrittweiten im lokalen Verfahrensfehler auftauchen. Tauchen Fehlerterme auf, die keine
Schrittweite enthalten, heißt das Verfahren nicht-konsistent mit dem gestellten Problem.
Die Konsistenzordnung gibt an, wieviel besser ein Verfahren wird, wenn man die jeweiligen Schrittweiten veringert. Halbieren wir in unserem Beispiel den Zeitschritt, so reduzieren
sich die zeitabhängigen Anteile des Verfahrensfehlers auf die Hälfte. Wählen wir den CrankNicolson-Faktor zu , so führt eine Halbierung des Zeitschrittes sogar zu einer Reduktion des Verfahrensfehlers auf ein Viertel. Die Konsistenzordnung sagt jedoch (noch) nichts
darüber aus, ob die diskrete Lösung gegen die exakte Lösung des Problems konvergiert.
Man ist demzufolge bemüht, Verfahren mit hoher Konsistenzordnung zu konstruieren. So
haben das Eulerverfahren und das implizite Verfahren die Konsistenzordnung 1, das LaxWendroff- und das Crank-Nicolson-Verfahren sind Verfahren 2. Ordnung und das RungeKutta-Verfahren besitzt sogar die Konsistenzordnung 4.
2.5 Mehrschrittverfahren
Bisher wurde zur Berechnung der Strömung zum Zeitpunkt nur der direkt vorangehende Zeitschritt verwendet. Mehrschrittverfahren verwenden auch noch weitere Zeitschritte.
Hierdurch kann man höhere zeitliche Konsistenzordnungen erreichen.
2.5.1 Numerische Integration
In diesem Abschnitt sei nur kurz an die drei grundsätzlichsten numerischen Integrationsverfahren erinnert. Die Mittelpunktsregel verwendet den Funktionswert auf der Intervallmitte
(2.10)
die Trapezregel das Mittel der Funktionswerte an den Integrationsgrenzen
(2.11)
KAPITEL 2. ZEITSCHRITTVERFAHREN
20
und die Simpsonregel ein gewichtetes Mittel aus den Funktionswerten an den Intervallgrenzen
und in der Intervallmitte:
(2.12)
2.5.2 Leap-Frog-Verfahren
Zur Herleitung eines Verfahrens, welches auch den Zeitpunkt auswertet, integrieren wir
die zeitabhängige Differentialgleichung (2.2) zwischen und :
Auf den ersten Term wenden wir den Hauptsatz der Differential- und Integralrechnung, auf
den zweiten die Mittelpunkteregel an:
bzw.
(2.13)
Das Verfahren wird in der Literatur als Leap-Frog-Verfahren bezeichnet, es ist ein explizites
Verfahren und benötigt zwei Zeitebenen zur Berechnung der gesuchten Werte auf einer neuen
Zeitebene. Es besitzt die Konsistenzordnung .
2.6 Zeitliche Mittlung der Grundgleichungen
Ein einfache Anwendung des Hauptsatzes der Differential- und Integralrechnung liefert:
(2.14)
Diese Gleichung gewährt einen tiefen Einblick in das, was wir einer Differentialgleichung
antun, wenn wir sie (in der Zeit) diskretisieren. Man berechnet eigentlich die Ableitung der
Mittlung der Funktion zwischen den Zeitschritten und , die durch
(2.15)
(2.16)
gegeben ist. In diesem Sinne ist die Diskretisierung
2.7. DAS ENTITY-RELATIONSHIP-MODELL DER ZEIT
K ra ft
21
W e rt
E in h e it
O rt
Z e it
Abbildung 2.1: Die Attribute der physikalischen Größe Kraft.
exakt. Daher liegt es nahe, die hydrodynamischen Grundgleichungen über die Zeitintervalle zu mitteln und dann die zeitliche Diskretisierung ohne Diskretisierungsfehler durchzuführen. Nach der zeitlichen Diskretisierung der Navier-Stokes-Gleichungen löst man eigentlich die Reynoldsgleichungen und man hat sich an dieser Stelle Gedanken über die durch die
Numerik erzeugten Reynoldsspannungen zu machen.
2.7 Das Entity-Relationship-Modell der Zeit
Im letzten Abschnitt dieses Kapitels wollen wir uns der Modellierung der Zeit aus der Sicht
der Informatik widmen.
2.7.1 Das Entity-Relationship-Modell
Der erste Schritt zur Modellierung von irgendetwas aus der Realwelt ist der Entwurf eines
konzeptionellen Modells. In der Literatur werden dabei eine Reihe von Architekturen für
den Entwurf des konzeptionenellen Modells vorgeschlagen. Grundlegend ist dabei der EntityRelationship-Ansatz (ER-Ansatz), der 1976 von Chen [7] vorgeschlagen und in der Folge
vielfach erweitert wurde. Die dahintersteckende Methodik wird in der modularen als auch
objektorientierten Softwarekonzeption und dem Entwurf von Datenbanken angewendet.
Die Modellbildung erfolgt dadurch, die Dinge der zu modellierenden Welt, die sogenannten
Entitäten in Entitätstypen zu ordnen. Ein Entitätstyp beschreibt also eine Menge gleichwertiger
materieller oder nichtmaterieller Dinge der realen Welt, die für den Problembereich relevant
sind.
Die Anwendung dieser Betrachtungsweise auf die Physik liefert die physikalischen Größen
als Entitätstypen. Geschwindigkeit und Kraft sind also Entitätstypen, während man als Instanz des Entitätstyps Kraft ansehen würde. In der graphischen Darstellung werden Entitätstypen durch Rechtecke symbolisiert.
Ein Entitätstyp wird durch eine Menge von Attributen charakterisiert. Eine Instanz oder Ausprägung des Entitätstyps erhält man durch die Zuordnung eines Wertes zu jedem Attribut
des Entitätstyp. Die Attribute werden im ER-Diagramm durch Striche am repräsentierenden
Rechteck, manchmal aber auch durch mit ihm verbundene Kreise dargestellt. Die Attribute
einer physikalischen Größe sind zunächst einmal ihr Zahlenwert und ihre Einheit, bei raumzeithaften Problemen kommen noch der Ort und der Zeitpunkt hinzu. Durch die Angabe dieser
vier Attribute wird eine Instanz der physikalischen Größe definiert.
Zwischen den Entitäten der realen Welt bestehen Beziehungen bzw. Relationen, zwischen den
Entitätstypen Beziehungstypen. Diese werden im ER-Diagramm durch Rauten dargestellt. Die
Beziehungstypen zwischen den physikalischen Größen sind die physikalischen Gesetze. Das
KAPITEL 2. ZEITSCHRITTVERFAHREN
22
K ra ft
2 .
N e w to n s c h e s
G e s e tz
M a sse
B e s c h le u n ig u n g
Abbildung 2.2: Das Entity-Relationship-Modell für das Newtonsche Gesetz.
Z e itz o n e
Z e it
1
b e s te h t
a u s
*
Z e itin te rv a ll
1
w ird
b e g re n z t
d u rc h
K a le n d e r
Ja h r
M o n a t
T a g
S tu n d e
Z e itp u n k t
2
S e k u n d e
M in u te
Abbildung 2.3: Das Entity-Relationship-Modell für die Zeit.
Newtonsche Gesetz stellt eine Beziehung zwischen den Entitäten Kraft , Masse ! und Beschleunigung dar, dessen Entity-Relationship-Diagramm in Abbildung 2.2 zu bestaunen ist.
Zwei Entitäten können nur mittels einer Beziehung verbunden werden. Gibt es zwischen zwei
Entitäten keine Beziehung, dann bleiben diese im ER-Diagramm unverbunden. Anders ist dies
bei Beziehungen, zwischen ihnen braucht keine Entität zu stehen. In der Welt der Physik kann
man gekoppelte Gesetze bzw. gekoppelte Gleichungssysteme durch Verbindungen von Relations ohne dazwischenliegende Entity modellieren.
Auch Beziehungstypen können durch Attribute näher charakterisiert werden. Man bezeichnet
sie dann als assoziative Entitätstypen und stellt ihre Attribute an einen die Beziehung substantivierenden Entitätstyp, der graphisch durch eine gerichtete Kante vom Beziehungs- zum
Entitätssymbol dargestellt wird.
2.7.2 Das ER-Modell der Zeit
Die Ontologie der Zeit geht von einem Zeitursprung und dem darauffolgenden Verlauf der Zeit
in Zeitschritten bzw. Zeitintervallen aus.
Erst der Anfang und das Ende eines Zeitintervalls werden durch Zeitpunkte markiert. Zusammenfassend ergibt sich das in Abbildung 2.3 dargestellte ER-Modell für die Zeit.
Zur Messung eines Zeitintervalles hat man seit menschengedenk die Rotation der Erde sowie
die Revolution der Erde um die Sonne als periodische Zeitmaße und die damit verbundenen
archetypischen Erfahrungen von Tag und Nacht und der Jahreszeiten verwendet. Die international gültige Definition des Zeitmaßes basiert heute auf der Spezifikation der Sekunde als
Vielfaches eines periodischen Vorganges aus dem Reich der Atomphysik.
Schwieriger als die Angabe eines Zeitintervalles ist die exakte Spezifikation von Zeitpunkten.
Dies liegt daran, daß es Zeitpunkte physikalisch nicht gibt, sie müssen also durch die Angabe
2.8. DIE LOGISTISCHE DIFFERENTIALGLEICHUNG ALS BEISPIEL
23
eines Zeitintervalles relativ zu einem Bezugsereignis spezifiziert werden. Das Jahr Null, sowie
eine willkürliche Unterteilung des Jahres in Monate und Wochen wird durch die Angabe eines Kalenders spezifiziert. Im westlichen Kulturkreis hat sich der Gregorianische Kalender so
etabliert, daß die Angabe desselben zumeist entfällt.
Wichtiger als die Angabe des Kalenders ist die Spezifikation der Zeitzone und gegebenenfalls
einer Sommerzeit, wodurch die Tageszeit an einem bestimmten geographischen Ort festgelegt
wird. Solche Angaben sind nicht nur bei Softwaresystemen aus dem Flugwesen, sondern auch
in der ozeano- und geographischen Simulation unentbehrlich.
In der Programmiersprache Java wird der Entitätstyp Zeitpunkt durch die Klasse java.util.Date
dargestellt. Die Genauigkeit reicht dabei bis zur Sekunde. Der Entitätstyp Zeitintervall gehört
nicht mehr zum Umfang der Standardsprache und soll daher durch eine Klasse TimeStep realisiert werden. Diese enthält die sechs privaten int-Komponenten years, months, days, hours,
minutes, seconds. Neben dem leeren Konstruktor sollte man Objekte der Klasse TimeStep
durch die Angabe der sechs Werte, sowie durch zwei Datums konstruieren können. Die Methode incrementTime liefert für ein gegebenes Objekt der Klasse Date ein um den Zeitschritt
inkrementiertes Objekt der Klasse Date.
2.8 Die logistische Differentialgleichung als Beispiel
Die logistische Differentialgleichung
" #
" # (2.17)
besteht auf der rechten Seite aus zwei Summanden: Für # wächst die Lösung in der
Zeit proportional zu sich selbst, wir haben es in diesem Fall mit der Differentialgleichung
des exponentiellen Wachstums zu tun. Der zweite Term begrenzt dieses Wachstum jedoch,
er wächst überproportional zu und bekommt gegenüber dem Wachstumsterm umso mehr
Gewicht, desto größer ist.
Die Differentialgleichung beschreibt also das logistische Wachstum, ein Wachstum also, welches durch irgendeinen begrenzenden Faktor, wie etwa die Verknappung von Ressourcen, gebremst wird.
Nach Heuser [15] kann man mit der logistischen Differentialgleichung so verschiedenartige
Dinge wie das Wachsen der Kriegslust in einem Land, die Ausbreitung von Gerüchten, Populationen und Sonnenblumen sowie die Gewichtszunahme bei Ratten modellieren. Aus dem
Bereich der Fließgewässerphysik sei das Wachstum von Riffeln und Dünen als Beispiel genannt, welches durch diese Differentialgleichung beschrieben wird.
Um eine numerische Lösung dieser Differentialgleichung brauchen wir uns nicht zu bemühen,
da sie eine analytische Lösung besitzt. Diese lautet
#
" # wobei die Anfangspopulation zum Zeitpunkt angibt.
Für strebt die Lösung für jeden Anfangswert gegen , im Laufe der Zeit
erreicht jede logistisch begrenzte Population einen Endwert, der umso größer ist, desto größer
die Wachstumsrate " und umso kleiner ist, desto kleiner die Sterberate # ist.
KAPITEL 2. ZEITSCHRITTVERFAHREN
24
Ist dieser Grenzwert erst einmal erreicht, ändert sich die Lösung nicht mehr, d.h. sie sollte und
ist zudem die Lösung der stationären Differentialgleichung
" #
Die Lösung der logistischen Differentialgleichung bezeichnet man als asymptotisch stabil, da
sie für jeden Anfangswert gegen die stationäre Lösung konvergiert. Diesen stationären Wert
kann man naheliegend auch als Attraktor bezeichnen.
Die logistische Differentialgleichung bietet sich als Validierungsbeispiel für Zeitschrittverfahren an, da wir ihre analytische Lösung kennen und mit dieser immer vergleichen können.
Starten wir also mit dem Eulerverfahren, es lautet für die logistische Differentialgleichung:
" #
Diese Form offenbart sofort ein erstes Problem. Beginnen wir mit dem Startwert , so
wird jeder Folgewert Null bleiben. Das Verfahren wird also nur dann gegen den richtigen asymptotischen Wert konvergieren, wenn man mit Startwerten ungleich Null beginnt.
Dies kann praktisch große Probleme aufwerfen, wenn man es mit Phänomenen zu tun hat,
die aus dem Nichts entstehen. Schon allein deshalb sollte das explizite Verfahren nicht zur
numerischen Lösung der logistischen Differentialgleichung herangezogen werden.
Wir wollen dennoch ein Lösungsprogramm schreiben, welches die Wachstums- und Sterberaten " und # , den Startwert und den Zeitschritt aus den Textfeldern einer Benutzeroberfläche übernimmt und die Ergebnisse graphisch darstellt. Diese kann man in Abbildung 2.4
bewundern. Sowohl numerische als auch analytische Lösung konvergieren gegen den richtigen
asymptotischen Wert, die numerische Lösung wächst dabei jedoch erst wesentlich später als
die analytische.
Bei einer Erhöhung des Zeitschrittes auf 2 s schwingt die numerische Lösung um die Asymptote, bei wird ihr Verhalten chaotisch und bei jedem Zeitschritt $ wächst die
numerische Lösung undarstellbar ins Unendliche. Das Verhalten des expliziten Verfahrens ist
also noch keineswegs befriedigend.
Wenden wir uns daher dem impliziten Verfahren zu. Nach kurzem Auflösen nach hat es
die Darstellung:
" # " # # Im Gegensatz zum Eulerverfahren ist das implizite Verfahren für beliebige Zeitschritte stabil
und liefert keine Oszillationen. Es konvergiert zudem auch für den Anfangswert gegen
die richtige Lösung. Es steigt allerdings wesentlich schneller als die analytische Lösung an
und ist daher im Anfangsbereich ebenfalls nicht brauchbar.
Als Verfahren der Konsistenzordnung zwei hatten wir das Lax-Wendroff-Verfahren und der
Ordnung vier das Runge-Kutta-Verfahren kennengelernt, die entsprechenden Schemata sollte man zur Übung selbst entwickeln. Die höhere Ordnung dieser Verfahren zahlt sich in der
besseren Reproduktion des Anstieges aus. Beide Verfahren neigen aber sehr schnell zu oszillierendem, chaotischem und instabilen Verhalten.
Wenn das explizite Verfahren instabil wird und das implizite Verfahren stabil bleibt, dann
sollten wir für den quadratischen Term eine Mischung aus explizitem und impliziten Verfahren
in der Form
2.8. DIE LOGISTISCHE DIFFERENTIALGLEICHUNG ALS BEISPIEL
25
Abbildung 2.4: Vergleich der mit dem Eulerverfahren gewonnenen numerischen (rot) und der
analytischen (grün) Lösung der logistischen Differentialgleichung
Abbildung 2.5: Vergleich der mit dem Eulerverfahren gewonnenen numerischen (rot) und der
analytischen (grün) Lösung der logistischen Differentialgleichung
26
KAPITEL 2. ZEITSCHRITTVERFAHREN
Abbildung 2.6: Vergleich der mit dem Eulerverfahren gewonnenen numerischen (rot) und der
analytischen (grün) Lösung der logistischen Differentialgleichung
Abbildung 2.7: Vergleich der mit dem impliziten Verfahren gewonnenen numerischen (rot)
und der analytischen (grün) Lösung der logistischen Differentialgleichung
2.9. ZUSAMMENFASSUNG
27
Abbildung 2.8: Vergleich der mit dem semiimpliziten Verfahren gewonnenen numerischen
(rot) und der analytischen (grün) Lösung der logistischen Differentialgleichung
" # probieren. Das Verfahren kann man als semiimplizit bezeichnen und es ist genau wie das implizite Verfahren stabil. Aus Abbildung 2.8 ist ersichtlich, daß die Lösung wie beim expliziten
Verfahren der analytischen hinterherhinkt.
In der Hydromechanik kommen nichtlineare quadratische Terme in den tiefen- und querschnittsgemittelten Impulsgleichungen als Sohlschubspannungs- bzw. Reibungsterme vor.
Diese werden in der Regel semiimplizit diskretisiert, weil hiermit der geringste Rechenaufwand bei garantierter Stabilität verbunden ist. Allerdings sollte man im Kopf behalten, daß das
numerische Verfahren die dämpfende Wirkung dieser Terme überschätzen kann, was wir hier
am verzögerten Anstieg der Lösung der logistischen Differentialgleichung gesehen haben.
2.9 Zusammenfassung
Wir haben in diesem Kapitel verschiedene Verfahren zur numerischen Behandlung der Zeitableitung kennengelernt. Als erstes Qualitätskriterium zur Beurteilung dieser Verfahren haben
wir eine möglichst hohe Konsistenzordnung gefordert. Unsere Untersuchungen haben ergeben,
daß das Euler- und das implizite Verfahren 1. Ordnung, Crank-Nicolson-, Lax-Wendroff- und
Leap-Frog-Verfahren 2. Ordnung und das Runge-Kutta-Verfahren sogar 4. Ordnung konsistent
sind.
28
KAPITEL 2. ZEITSCHRITTVERFAHREN
Abbildung 2.9: Vergleich der mit dem Lax-Wendroff-Verfahren gewonnenen numerischen
(rot) und der analytischen (grün) Lösung der logistischen Differentialgleichung
Abbildung 2.10: Vergleich der mit dem Runge-Kutta-Verfahren gewonnenen numerischen
(rot) und der analytischen (grün) Lösung der logistischen Differentialgleichung
Kapitel 3
Stabilität
In der nichtlinearen Dynamik bezeichnet man Systeme dann als stabil, wenn sich für zwei
sehr nahe beieinanderliegenden Anfangszuständen und die Zustände zu
späteren Zeiten höchstens stetig voneinander entfernen, d.h. es gibt zwei sehr kleine Zahlen Æ
und % mit
Æ % Æ Dieses Stabilitätskriterium sieht dem berüchtigten % Æ -Kriterium für Stetigkeit sehr ähnlich. Bei der Stetigkeit steht die Wahl des % im Vordergrund, d.h. zum Zeitpunkt sind die
Unterschiede zwischen und beliebig klein, wenn nur die Anfangswerte hinreichend nahe
aneinander liegen. Der Stabilitätsbegriff weicht diese Anforderung ein wenig auf, da er ledig fordert, daß zwei im Anfang sehr ähnliche Funktionen sich zu späteren Zeitpunkten nur
allmählich unterschiedlich verhalten.
Turbulente Strömungen sind in diesem Sinne stabil, da die turbulenten Fluktuationen in einem
gewissen Intervall um einen wohldefinierten Mittelwert schwanken.
Ein System heißt asymptotisch stabil, wenn für zwei sehr nahe beieinanderliegenden Anfangszuständen sich spätere Zustände irgendwann einmal nicht mehr unterscheiden, d.h. es
gibt irgendein sehr kleines Æ , so daß
Æ für
Wir werden in diesem Kapitel das Phänomen Stabilität bei dynamischen numerischen Verfahren im Rahmen der Operatorentheorie beschreiben und allgemeine Stabilitätskriterien
für explizite und implizite Verfahren aufstellen. Da konsistente numerische Verfahren nicht
notwendig stabil sein müssen, werden wir als weiteres Gütekriterium die Stabilität von Zeitschrittverfahren untersuchen und in den Begriff der Konvergenz einmünden lassen.
Es sei an dieser Stelle betont, daß der Begriff der Stabilität nichts mit dem des Chaos zu tun hat
oder etwa dessen Gegenteil ist. Chaotische Systeme sind in der Regel stabil. Ein Zeichen des
Chaos ist die Sensitivität auf Anfangslösungen, weisen diese nur winzige Unterschiede auf,
dann verhalten sich die Lösungen chaotischer Systeme zu späteren Zeitpunkten vollkommen
anders. Dies ist in stabilen Systemen keinesfalls ausgeschlossen, die Lösungen entfernen sich
hier lediglich nicht abrupt voneinander.
29
KAPITEL 3. STABILITÄT
30
u ,v
u
v
e
d
t
Abbildung 3.1: Stabilität turbulenter Strömungen. Die Anfangswerte unterscheiden sich um
den Wert Æ . Unabhängig davon bleiben die Fluktuationen im Intervall %.
3.1 Stetige Abhängigkeit von den Eingangsdaten
Die Stabilität eines numerischen Verfahrens ist eng mit der Frage verbunden, ob die Eigenschaft der stetigen Abhängigkeit von den Eingangsdaten beim Übergang vom kontinuierlichen
auf das diskrete Problem erhalten bleibt.
3.1.1 Stetige Abhängigkeit beim kontinuierlichen Problem
Eine der wichtigsten Eigenschaften der Lösungen der allgemeinen Evolutionsgleichung ist ihre stetige Abhängigkeit von den Eingangsdaten. Umgangssprachlich bedeutet dies, daß die
Lösungen in ihrem zeitlichen Verhalten keine abrupten Änderungen erfahren, sie also der Anfangslösung umso ähnlicher sind, desto näher man sich am Anfangszeitpunkt befindet.
Setzen wir zur Abkürzung
& (3.1)
so folgt für die Lösungen der homogenen Evolutionsgleichung die Abschätzung
& Zu jedem beliebigen Zeitpunkt ist die Norm der Lösung durch den Faktor & mit der Norm
der Anfangslösung verbunden. Dieser Faktor ist eine stetige Funktion der Zeit , ändert sich
also ebenfalls niemals abrupt. Somit kann sich die Lösung im Laufe der Zeit zwar beliebig
weit von der Anfangslösung entfernen und irgendwann vollkommen anders aussehen, diese
Änderung geschieht aber stetig und nicht aus heiterem Himmel.
Wir wollen nun einen wichtigen Spezialfall untersuchen, der in der Strömungsmechanik sehr
oft vorkommt, obwohl seine Voraussetzungen doch sehr einschränkend zu sein scheinen.
3.1. STETIGE ABHÄNGIGKEIT VON DEN EINGANGSDATEN
31
Satz: Falls der Operator linear ist und nur Eigenwerte ' mit nichtnegativem Realteil besitzt
und die zugehörigen Eigenvektoren eine Orthonormalbasis des Lösungsraumes bilden, gilt
& bzw
Da die Norm bisher nicht spezifiziert wurde, gilt der Satz für jede Operatornorm. Nehmen wir
hier z.B. die Maximumnorm auf einem Intervall
dann besagt der Satz, daß die Extremwerte der Lösung im Laufe der Zeit nie größer werden,
die Lösungsfunktion also immer in den anfänglichen Schranken bleibt.
Beweis: Sei ( die Orthonormalbasis des L ösungsraumes. Die Anfangswerte lassen sich
dann als
(
darstellen. Somit folgt
... und da (
(
( (
( (
( (
(
(
*)
... und da für (
( (
( (
' ( ...
... und da , weil ' , denn die Zeit ist größer Null und die Eigenwerte
waren laut Voraussetzung gr ößer Null, folgt
KAPITEL 3. STABILITÄT
32
(
bzw.
womit die Behauptung folgt.
Hat der Operator der homogenen Evolutionsgleichung nur positive Eigenwerte, so ist die
Wirkung der Evolutionsgleichung in gewisser Form dämpfend, denn die Norm der Lösung
nimmt mit der Zeit immer mehr ab.
Wir betrachten nun Abschätzungen für die inhomogene Evolutionsgleichung. Mit der
Abkürzung
&
¼
¼
(3.2)
gilt für die Lösung der inhomogenen Evolutionsgleichung die Abschätzung:
& &
¼
Und schließlich folgt für die nichtlineare Form:
& &
&
3.1.2 Stetige Abhängigkeit beim diskretisierten Problem
Die stetige Abhängigkeit der kontinuierlichen Lösung von den Eingangsdaten hängt bei der
allgemeinen Evolutionsgleichung mit der Erfüllung der Abschätzung
& &
(3.3)
zusammen. Somit können wir von der zeitlich diskretisierten Lösung
& & verlangen. Diese Bedingung garantiert zudem die Stabilität des Zeitschrittverfahrens.
(3.4)
3.2. DER VERFAHRENSOPERATOR
33
3.2 Der Verfahrensoperator
Einschrittverfahren lassen sich mit der Operatorenschreibweise auf die folgende Form bringen:
+ (
(3.5)
Dabei wird + auch Verfahrensoperator genannt.
Schreibt man das Crank-Nicolson-Verfahren im Rahmen der Operatorentheorie als
) ) dann sieht man leicht, daß dem Eulerverfahren ( )
dem impliziten Verfahren (
+
) )
+
) und dem Crank-Nicolson-Verfahren
+
) )
als Verfahrensoperatoren zugeordnet sind.
Geht man sukzessive von den Anfangsbedingungen bis zum Zeitschritt
Verfahren die Form
+ , so erhält das
+ (
(3.6)
Ein Vergleich mit der Stetigkeitsbedingung liefert die Stabilitätsbedingung
+ &
Ist &
,
(3.7)
so ist diese Bedingung erfüllt, wenn
+ , (3.8)
ist. Falls & und + $ , dann ist die Anzahl der Zeitschritte begrenzt, ab einem
gewissen Zeitschritt wäre das Verfahren instabil. Stabilität bedeutet dann, daß die Wirkung
des Verfahrensoperators nicht unbeschränkt ausartet, sondern die Lösung in einem gewissen
Sinne kontrahierend ist.
Um im folgenden Stabilitätsanalysen durchzuführen, benötigt man noch zwei wesentliche
Sätze. Der erste ist der Form nach schon aus der Algebra bekannt.
Satz: Neumannsche Reihe. Sei ein beschränkter linearer Operator mit Identitätsoperator. Dann existiert der folgende Operator
) , und ) der
(3.9)
KAPITEL 3. STABILITÄT
34
Der nächste Satz ist ein mächtiges Werkzeug für die Analyse numerischer Verfahren.
Satz von Kellogg: Sei
für alle des Hilbertraumes und - . Dann gilt
) - ) -
(3.10)
Sei die dargestellte Theorie auf das Crank-Nicolson-Verfahren
angewendet. Die Gleichung läßt sich zu
)
) umstellen, und definiert man den Verfahrensoperator
+
)
)
so ergibt sich die Stabilität des Crank-Nicolson-Verfahrens sofort aus Kelloggs Satz, wenn
nur der Operator die dortigen Voraussetzungen erfüllt. Die Konsistenzanalyse zeigt, daß das
Verfahren zweiter Ordnung in der Zeit ist. Noch allgemeiner kann man sogar folgenden Satz
beweisen:
Satz: Das Crank-Nicolson-Verfahren ist für beliebige Zeitschritte stabil, wenn und positiv-semidefinit ist.
3.2.1 Matrixnormen
Da in den folgenden Kapiteln auch der Raum diskretisiert wird, stellen die Lösungen und
selbst Vektoren dar. Für lineare Verfahren kann der Verfahrensoperator + somit mit einer
Matrix identifiziert werden. Wir benötigen zur Untersuchung der Stabilität Rechenvorschriften
zur Bildung von Normen von Matrizen. Analog zu den Normen von Vektoren definieren wir
dererlei drei:
Die Spaltensummennorm
+ die Quadratsummen- oder Schurnorm
+
und die Zeilensummennorm
+ + + + (3.11)
(3.12)
(3.13)
3.3. ALLGEMEINE STABILITÄTSKRITERIEN
35
3.2.2 Spektralradius und Stabilität
Der Spektralradius eines Operators bzw. einer Matrix ist durch den Grenzwert
+ + (3.14)
gegeben. Er ist unabhängig von der gewählten Norm und ist die untere Schranke aller Operatornormen, denn es gilt:
+ + (3.15)
Satz: Besitzt der lineare Operator + N Eigenwerte ' und N linear unabhängige Eigenvektoren
d.h.
+ ' so gilt:
+ ' d.h. der Spekralradius ist der Betrag des betragsmäßig größten Eigenwertes.
Da die linear unabhängigen Eigenvektoren eine Basis des Lösungsraumes bilden, läßt sich
die Lösung zum Zeitschritt als Linearkombination
.
darstellen. Für die Lösung gilt somit
+ '
. Damit die Lösung zum Zeitschritt kleiner als die zum Zeitschritt ist, muß der
betragsmäßig größte Eigenwert ' von + kleiner als eins sein. Damit haben wir folgenden
Satz gezeigt:
Satz: Ein numerisches Verfahren, dessen Verfahrensoperator linear ist und dessen Eigenvektoren eine Basis des Lösungsraumes bilden, ist dann stabil, wenn der Spektralradius + des
zugehörigen linearen Verfahrensoperators + kleiner als eins ist.
Hiermit haben wir das Stabilitätsproblem auf die Berechnung von Eigenwerten des Verfahrensoperators, d.h. von Matrizen zurückgeführt.
3.3 Allgemeine Stabilit ätskriterien
Wir wollen nun sehr allgemeine Stabilitätskriterien entwickeln, die wir erst nach der
Einführung der konkreten numerischen Verfahren bei der Stabilitätsanalyse anwenden werden.
KAPITEL 3. STABILITÄT
36
Sie geben aber schon jetzt Hinweise für die Konstruktion stabiler numerischer Zeitschrittverfahren. Insbesondere werden wir feststellen, daß implizite Anteile die Stabilität mit zunehmendem Zeitschritt förden und explizite Anteile die Instabilität begünstigen.
3.3.1 Stabilitätskriterien für explizite Einschrittverfahren
Der Weg, wie das Stabilitätskriterium für das Crank-Nicolson-Verfahren gewonnen wurde,
läßt sich verallgemeinern. Dazu betrachte man ein explizites Einschrittverfahren der Form
(
wobei und Operatoren sind, die das Verfahren repräsentieren. Man sieht leicht, daß der
Verfahrensoperator durch
+
gegeben ist. Seine Norm kann als
+ abgeschätzt werden, so daß + , für
(3.16)
gilt. Offensichtlich wächst mit steigendem Zeitschritt die Gefahr, das Stabilitätskriterium
zu verletzen.
3.3.2 Stabilitätskriterien für implizite Einschrittverfahren
Wir wenden uns nun einem allgemeinen impliziten Verfahren in der Form
(
zu. Der Verfahrensoperator ist durch
+
gegeben. Seine Norm ist
) + ) ) so daß man nach kurzer Rechnung zeigen kann, daß + , für
(3.17)
gilt. Offensichtlich sinkt mit steigendem Zeitschritt die Gefahr, das Stabilitätskriterium zu
verletzen, implizite Verfahren sind daher also wesentlich stabiler als explizite.
3.4. ITERATIVE UND PRÄDIKTOR-KORREKTOR-VERFAHREN
37
3.4 Iterative und Prädiktor-Korrektor-Verfahren
Die Verwendung von impliziten Anteilen in Zeitschrittverfahren erhöht zum einen in der
Kombination mit expliziten Anteilen die Genauigkeit und zum anderen die Stabilität des
Verfahrens. Diese Vorteile werden mit höheren Rechenkosten bezahlt, da bei linearen Gleichungen ein lineares Gleichungssystem und bei nichtlinearen Gleichungen ein nichtlineares
Gleichungssystem zu lösen ist. Iterative und Prädiktor-Korrektor-Verfahren versuchen den Rechenaufwand zu reduzieren, indem quasi-implizite Lösungen bestimmt werden. Dazu wird ein
Zeitschritt mehrfach berechnet, wobei das Ergebnis jeweils zur Berechnung von wiederverwendet wird.
Betrachten wir dazu das explizite Einschrittverfahren
bzw.
) Das Ergebnis wird zur Berechnung von wiederverwendet:
bzw.
) ) Wird das Spiel unendlich oft wiederholt, so konvergiert der Verfahrensoperator gegen
) d.h. gegen den impliziten Verfahrensoperator. Der Satz über die Neumannsche Reihe fordert,
daß dazu die Bedingung
, erfüllt sein muß. Die Stabilität des Verfahrens ist durch die Stabilität des expliziten Einschrittverfahrens begrenzt, also sehr niedrig. Es liegt daher nahe, in den einzelnen Iterationsschritten
unterschiedliche Verfahren zu verwenden, womit wir bei der Klasse der Prädiktor-KorrektorVerfahren angelangt sind. Dabei wird im ersten Iterationsschritt – dem Prädiktor – ein stabiles
Verfahren angewendet und im zweiten Iterationschritt ein Verfahren, welches die Genauigkeit
steigert.
3.5 Operator-Splitting-Verfahren
Die Operator-Splitting-Technik zerlegt die urprünglichen Differentialgleichungen in verschiedene Anteile, die dann nacheinander mit den mathematischen Eigenschaften dieser Anteile
entsprechend angepaßten Verfahren gelöst werden.
KAPITEL 3. STABILITÄT
38
Dazu sei die Differentialgleichung in der Form
darstellbar. kann z.B. die advektiven und die diffusiven Terme in einer Transportglei-
chung repräsentieren. Im ersten Schritt wird die Gleichung
gelöst. Das Zwischenergebnis wird als Startwert für den zweiten Schritt verwendet:
Die Addition der beiden Teilschrittgleichungen führt offensichtlich zu einer konsistenten Diskretisierung des Ausgangsproblems. In den einzelnen Teilschritten können verschiedene Integrationsverfahren für die Zeit angewendet werden. Zur Untersuchung der numerischen Eigenschaften des Operator-Splittings wenden wir in beiden Teilschritten ein Crank-NicolsonVerfahren an:
Im Rahmen der Operatorentheorie lassen sich diese Gleichungen umschreiben zu
) ) ) ) und das Gesamtverfahren kann durch
) ) ) ) dargestellt werden, wobei der Verfahrensoperator durch
+
) ) )
)
gegeben ist. Offensichtlich setzt sich der Verfahrensoperator aus dem Hintereinanderausführen
der einzelnen Crank-Nicolson-Verfahren zusammen. Es ist somit leicht einzusehen, daß sich
der Gesamtverfahrensoperator beim Operator-Splitting aus den Verfahrensoperatoren der
Einzelschritte zusammensetzt. Mit der Beziehung (15.5) erhält man somit den wichtigen
Satz: Das gesamte Operator-Splitting-Verfahren ist dann stabil, wenn jeder Einzelschritt stabil
ist.
3.6. KONVERGENZ
39
Sei nun der lokale Verfahrensfehler in der Zeit betrachtet. Wir fragen, ob auch das
Operator-Splitting die Konsistenzordnung 2 aufrecht erhalten kann. Dazu wird der Verfahrensoperator + für unter Anwendung der Neumannschen Reihe bis zur ersten Ordnung
entwickelt:
+
) ) ) ) ) Das weitere Vorgehen hängt nun entscheidend davon ab, ob die Operatoren
tauschbar sind. Denn nur dann läßt sich
und ver-
schreiben und somit der Verfahrensoperator des Lax-Wendroff-Verfahrens:
+
) Damit erhält man das folgende
Resultat: Das Crank-Nicolson Operator-Splitting ist stabil, wenn die Teilschritte stabil sind
und besitzt die Konsistenzordnung 2, wenn die Teilschritte vertauschbar sind.
Der Vorteil des Operator-Splitting-Verfahrens liegt darin, daß in jedem Schritt genau auf die
Gleichung angepaßte Verfahren angewendet werden können, der Nachteil ist die maximale
Konsistenzordnung von .
3.6 Konvergenz
Die wichtigste Anforderung an ein numerisches Verfahren ist die Konvergenz gegen die exakte
Lösung. Für soll sich die numerische Lösung an die exakte Lösung annähern.
Betrachten wir hierzu die Differenz zwischen numerischer und exakter Lösung / :
/ Da die diskrete Lösung die diskretisierte Gleichung löst, gilt
/ /
/
/ bzw. in Operatorschreibweise
) / ) / und somit:
/ + / Wendet man diesen Zusammenhang vom Ausgangszeitpunkt induktiv bis an, so ergibt
sich die Fehlerdarstellung:
/ + / + KAPITEL 3. STABILITÄT
40
Für die Norm des Fehlers folgt:
/
,
+ + + / + / / / / Zum Zeitschritt geht das Residuum mit der Ordnung des Verfahrensfehlers zu Null,
wenn die Anfangswerte exakt sind. Wir haben somit das Theorem von Lax bewiesen:
Theorem von Lax: Ist ein Zeitschrittverfahren konsistent mit der mit der Ordnung und stabil, so konvergiert die numerische Lösung bei kleiner werdenden Zeitschritt mit der Ordnung
gegen die exakte Lösung.
Kapitel 4
Die Methode der Finiten Differenzen
Die Methode der Finiten Differenzen (FD) ist die älteste Methode zur numerischen Lösung
von partiellen Differentialgleichungen.
Auf einem Gebiet sei eine partielle Differentialgleichung der allgemeinen Form
(4.1)
zu lösen. stellt einen Differentialoperator dar, ist die gesuchte Funktion und enthält die
Inhomogenitäten.
Bei der Methode der Finiten Differenzen wollen wir drei Diskretisierungsschritte unterscheiden.
Im ersten Schritt wird das Lösungsgebiet mit einem rechtwinkligen nicht notwendig äquidistanten Gitter überdeckt (Abb. 4.1). Dabei werden Innenknoten (volle Kreise) von Randknoten
(leere Kreise) unterschieden. Das Lösungsgebiet wird also durch eine Menge von Punkten
0 1 &
ersetzt, die man auch Knoten nennt.
Im zweiten Schritt wird die gesuchte Lösungsfunktion diskretisiert, d.h. sie wird durch den
Vektor ihrer Werte auf den Knoten ersetzt. Den zu der kontinuierlichen Funktion gehörigen Vektor von Knotenwerten bezeichnet man auch als Gitterfunktion.
Ebenso verfährt man mit der rechten Seite , auch sie wird mit Hilfe einer Knotenfunktion auf die Werte an den Knoten reduziert.
Da aus der Lösungsfunktion und der Inhomogenität nun Vektoren geworden sind, muß aus dem
Differentialoperator nun ein algebraischer Operator werden, der bei linearen Problemen
durch eine Matrix dargestellt werden kann.
Mit diesem dritten Schritt wollen wir uns im folgenden beschäftigen.
4.1 Einfache Differenzenquotienten
Zur Herleitung des algebraischen Operators werden alle Ableitungen durch Differenzenquotienten der Knotenwerte ersetzt. Um einen Differenzenquotienten die erste Ableitung in
x-Richtung zu erhalten, entwickeln wir die Funktion am Knoten 0 1 in eine Taylorreihe
41
KAPITEL 4. DIE METHODE DER FINITEN DIFFERENZEN
42
0 0 0
1
1
1
Abbildung 4.1: Diskretisierung eines 2D-Gebietes für Finite Differenzen
0 1 erhalten wir
Für den Entwicklungsknoten
Löst man die erste Taylorreihe nach auf, so erhält man als Approximation der Ableitung
(4.2)
(4.3)
welche man als vordere Differenzenapproximation bezeichnet. Andererseits erhält man aus
der zweiten Taylorreihe das Schema
welches man als rückwärtige Differenzenapproximation bezeichnet.
Zieht man die zweite von der ersten Taylorreihe ab, so erhält man in analoger Weise die zentrale Differenzenapproximation
(4.4)
Addiert man schließlich beide Reihen, so ergibt sich eine Differenzenapproximation für die
zweite Ableitung in der Form
(4.5)
4.1. EINFACHE DIFFERENZENQUOTIENTEN
43
die bei einer äquidistanten Diskretisierung in
(4.6)
übergeht. Zur Herleitung von Differenzenverfahren für Ableitungen dritter Ordnung benötigt
man zusätzlich noch die Tayloreihen für und .
In derselben Weise lassen sich die Ableitungen in andere Richtungen approximieren. Für den
zweidimensionalen Fall kann die ersten Ableitungen in y-Richtung am Knoten 0 1 als vordere Differenz
oder als zentrale Differenz
dargestellt werden. Für die zweite Ableitung in y-Richtung ergibt sich
4.1.1 Differenzenverfahren in Matrixschreibweise
Um ein lineares Gleichungssystem zu erhalten, muß die mehrdimensionale Indizierung der
Gitterknoten auf eine eindimensionale Indizierung abgebildet werden, anders
gesagt, die Knoten müssen durchnumeriert werden. Für Knoten, die außerhalb des Lösungsgebietes liegen, gibt es prinzipiell zwei verschiedene Behandlungsmöglichkeiten:
’Trockene’ d.h. außerhalb des Berechnungsgebietes liegende Knoten werden nicht mitindiziert und fallen damit aus der Rechnung heraus.
’Trockene’ Knoten werden maskiert, d.h. es wird ein INTEGER-Feld der Länge aller
Gitterknoten angelegt, welches eine Null für trockene und eine Eins für Berechnungsknoten enthält. Das Ergebnis wird an geeigneten Stellen mit dieser Maske multipliziert.
Für die Numerierung der übrigen Knoten haben sich drei Varianten durchgesetzt:
Die Knoten werden reihenweise (oder spaltenweise) d.h. in Richtung der Koordinatenachsen durchnumeriert.
Die Knoten nach Schrägzeilen d.h. in Diagonalenrichtung durchnumeriert.
Die Knoten werden schachbrettartig numeriert. Dazu teilt man das Gitter wie ein
Schachbrett in schwarze und weiße Knoten ein und numeriert zuerst die schwarzen Knoten reihenweise und dann die weißen entsprechend oder aber umgekehrt.
KAPITEL 4. DIE METHODE DER FINITEN DIFFERENZEN
44
Die Numerierung der Gitterknoten ist ein besonders diffiziles Problem, für das der Programmierer neben einem glücklichen Händchen viel Erfahrung haben muß, denn sie beeinflußt die
Performance des Programmes entscheidend. An dieser Stelle können keine Pauschalaussagen
gemacht werden, da eine optimale Wahl nicht nur von der zu lösenden partiellen Differentialgleichung, sondern auch vom Gleichungslöser, vom Aufbau des gesamten Softwarepaketes
und vom Computertyp (skalare, Vektor- oder Parallelmaschine) abhängig ist.
Nachdem der Lösungsvektor und die rechte Seite des Gleichungssystems durch Umnumerieren der Knoten des Gitters erzeugt wurden, bleibt der Aufbau der Systemmatrix . Hierfür
sei ein eindimensionales Problem betrachtet. Man sieht leicht, daß dann die Matrixdarstellung
für die erste Ableitung, die durch zentrale Differenzen approximiert ist, durch
..
..
.
..
..
.
..
.
..
.
..
.
.
.
gegeben ist. Für die zweite Ableitung gilt
..
..
.
.
..
..
..
.
.
.
..
..
.
.
Alle Leser sollten zur Übung die Matrixdarstellung der vorderen und rückwärtigen Differenzen
konstruieren.
4.1.2 Die stationäre Transportgleichung
Die Lösung von Differentialgleichungen mit Hilfe von Finiten Differenzen soll nun an einem einfachen Beispiel konkretisiert werden. Dazu betrachten wir die stationäre Transportgleichung
Zur numerischen Lösung der Gleichung stehen auf den ersten Blick (mindestens) drei Verfahren zur Verfügung, da wir für die erste Ableitung zentrale Differenzen
vordere Differenzen
und rückwärtige Differenzen
4.2. KONSISTENZ
45
kennengelernt hatten. Alle drei Verfahren lassen sich durch die Einführung des sogenannten
Upwind-Parameters 2 verallgemeinern:
2
2
(4.7)
Für 2 ergeben sich rückwärtige Differenzen etc.
Um die Qual der Wahl zu erleichtern, wollen wir die Kriterien Konsistenz und Stabilität untersuchen.
4.2 Konsistenz
Als erstes Kriterium für die Güte der FD-Verfahren untersuchen wir wieder den lokalen Verfahrensfehler, den wir erhalten, wenn wir die exakte Lösung in das Differenzenverfahren einsetzen.
Konkretisieren wir dies für die stationäre Transportgleichung. Seien mit , und ausnahmsweise und nur hier die exakte Lösung an den entsprechenden Knoten bezeichnet. Wir
entwickeln sie in Taylorreihen
und setzen diese in das Verfahren (4.7) ein:
2
Zur Bestimmung des lokalen Verfahrensfehlers d.h. dem Unterschied zwischen Differentialund Differenzengleichung mit exakter Lösung ziehen wir die zu lösenden Differentialgleichung ab:
2
Für 2 bzw. 2 besitzt das Verfahren die Konsistenzordnung 1, für 2 sogar die
Konsistenzordnung 2. Die zentralen Differenzen wären somit den vorderen und rückwärtigen
Differenzen vorzuziehen.
4.3 Stabilität
Im Rahmen der Zeitschrittverfahren hatten wir gelernt, daß Instabiltäten dann auftreten, wenn
das numerische Verfahren Lösungen produziert, die die Schranken des Wertebereiches der
exakten Lösung verletzen.
KAPITEL 4. DIE METHODE DER FINITEN DIFFERENZEN
46
Solche Schranken existieren auch auf der Ortsebene: Der Wertebereich elliptischer Differentialgleichungen ist durch die Randwerte (nach oben oder unten) beschränkt. Verletzen die
Lösungen eines numerischen Verfahrens diese Schranken, so treten wieder Instabilitäten auf,
die allerdings nicht notwendig zum Abbruch des Computerprogramms führen müssen.
Diese Instabilitäten sind in Form von Oszillationen bei Transportgleichungen zu beobachten, wenn die Advektion gegenüber der Diffusion so sehr überwiegt, daß die oberstromseitige
Randbedingung bei der gewählten Diskretisierung nicht mehr sinnvoll realisiert werden kann.
Eine charakteristische Größe, die das Verhältnis von Advektion zu Diffusion beschreibt, ist die
sogenannte Pecletzahl (auch manchmal Gitterreynoldszahl genannt), die durch
3 (4.8)
gegeben ist. In einem vereinfachten Beispiel [10], [23],[24],[25] wollen wir sowohl die Oszillationen kennenlernen, als auch ein Kriterium für ihr Verschwinden herleiten. Dazu betrachten
wir die stationäre Transportgleichung
Die Gleichung ist elliptisch; alle Lösungen der stationären Transportgleichung nehmen ihre
Extrema auf den Rändern des Lösungsgebietes an.
Werden die Randbedingungen
gewählt, dann ergibt sich eine analytische Lösung als
Zur numerischen Lösung der Gleichung verwenden wir zentrale Differenzen, wir erhalten das
Schema
mit , 1
und
Indem wir die Verfahrensgleichug mit multiplizieren und durch dividieren, erhält man
die dimensionslose Form
3
3
In Abb. 4.2 sehen wir den Vergleich zwischen numerischer und analytischer Lösung für eine
Pecletzahl von 6.67. Der Vergleich zeigt, daß
die analytische Lösung monoton ist und die numerische Lösung oszillierend,
4.3. STABILITÄT
47
1 .0 0 0
.
.
P e = 6 .6 7
.
0 .7 5 0
.
.
0 .5 0 0
.
.
.
0 .2 5 0
.
-
.
.
0 .0 0 0
-
.
.
.
.
.
0 .0 0 0
.
.
.
.
.
.
.
.
.
.
.
.
0 .7 5 0
0 .5 0 0
0 .2 5 0
1 .0 0 0
n u m e r is c h e L ö s u n g
a n a ly tis c h e L ö s u n g
Abbildung 4.2: Numerische (durchgezogen) und analytische (gestrichelt) Lösung für zentrale
Differenzen, , , .
die analytische Lösung nicht-negativ ist und die numerische Lösung auch negative Werte
annimmt. Das Verfahren ist somit nicht monoton.
Numerische und analytische Lösung haben somit ein qualitativ vollkommen verschiedenes
Verhalten.
In der Regel kennt man aber eine analytische Lösung nicht und da die numerische Lösung
nicht generell instabil ist, bleibt dem Anwender eines Codes die Beantwortung der kniffeligen
Frage, ob die Oszillationen tatsächlich physikalischer Natur sind.
Zur Analyse des Problems schreiben wir das Verfahren inklusive Randbedingungen wieder in
Matrixschreibweise:
..
..
..
.
..
.
..
.
..
.
..
.
.
.
..
.
..
.
..
.
..
.
..
.
Beim Auftreten von Oszillationen wird die numerische Lösung auch negativ. Es läßt sich nun
zeigen, daß dies nicht der Fall ist, wenn die Matrix diagonal dominant ist, d.h. wenn
$ 3
3
Diese Bedingung ist genau dann erfüllt, wenn für die Pecletzahl die Bedingung
3 (4.9)
erfüllt ist.
Wenden wir auf denselben Testfall Upstream-Differenzen in der Form
KAPITEL 4. DIE METHODE DER FINITEN DIFFERENZEN
48
1 .0 0 0
.
.
-
P e = 6 .6 7
0 .5 0 0
-
.
0 .7 5 0
.
.
.
.
.
0 .2 5 0
.
-
.
.
0 .0 0 0
-
.
.
.
.
0 .0 0 0
.
.
.
.
.
.
.
.
.
.
.
.
0 .7 5 0
0 .5 0 0
0 .2 5 0
.
1 .0 0 0
n u m e r is c h e L ö s u n g
a n a ly tis c h e L ö s u n g
Abbildung 4.3: Numerische (durchgezogen) und analytische (gestrichelt) Lösung für
Upstream-Differenzen, , , .
an, so verschwinden die Oszillationen, wie man in Abb. 4.3 erkennt.
Dies ist insofern verwunderlich, da das zentrale Verfahren mit der Konsistenzordnung wesentlich genauer ist, als das Upstream-Verfahren ( ). Um dieses Verhalten zu verstehen, verallgemeinern wir zentrale und Upstream-Differenzen durch das Schema
2
2
Für 2-Werte von 1, 0 und 1/2 ergeben sich in der Reihenfolge das Upstream, das Downstream
sowie das zentrale Verfahren. Schreiben wir die Gleichung nochmals ein wenig in der Form
3 2 um, so zeigt sich, daß dem Verfahren ein künstlicher Diffusionsterm der Form 3 2 zugefügt wurde, der nur im zentralen Fall nicht wirkt. Je größer also die Pecletzahl ist, d.h. umso
kleiner die echte Diffusion gegenüber der Advektion ist, desto mehr künstliche Diffusion wird
der Lösung beim Upstream-Verfahren zugefügt. Hierdurch treten zwar keine Oszillationen auf,
die Lösung wird allerdings durch den künstlichen Diffusionsterm stark gedämpft.
4.4 Verallgemeinerungen
Zur Lösung mehrdimensionaler Probleme werden die Grundverfahren durch die entsprechenden zusätzlich auftretenden Terme erweitert.
Wegen der hohen Rechenleistungen, die allein durch die Anzahl der Knoten entsteht, werden in
3D-Simulationen vielfach explizite Verfahren verwendet, wobei hier vor allem das Upstreamund Lax-Wendroff-Verfahren für die Advektion gekoppelt mit dem FTCS-Verfahren für die
Diffusion zu nennen sind. Bei der impliziten Lösung werden Crank-Nicolson-Verfahren verwendet.
Bei den bisher diskutierten Verfahren wurde vorausgesetzt, daß die Gitterabstände äquidistant
( 4, 4) sind. Ersetzt man die konstanten Gitterabstände durch räumlich
4.4. VERALLGEMEINERUNGEN
49
veränderliche in der Form , so lassen sich alle Verfahren auf nicht-äquidistante
Gitter verallgemeinern. Es ist so möglich, interessante Simulationsgebiete höher aufzulösen.
Allerdings breitet sich die höhere Auflösung in Querrichtung über das gesamte Simulationsgebiet fort.
50
KAPITEL 4. DIE METHODE DER FINITEN DIFFERENZEN
Kapitel 5
Die Advektionsgleichung
Wir beginnen mit der Advektionsgleichung, die Bestandteil jeder Transportgleichung ist:
(5.1)
Die Variablen dieser Gleichung beschreiben ein dreidimensionales Geschwindigkeitsfeld 5 , von dem wir der Einfachheit halber annehmen, daß es homogen und zeitlich konstant ist. Mit diesem Geschwindigkeitsfeld wird eine Konzentrationsverteilung bewegt.
5.1 Theorie der Advektionsgleichung
Zunächst wenden wir uns den theoretischen Aspekten dieser Gleichung zu. Man bestätige,
daß es sich um eine lineare Differentialgleichung handelt. Es wird gleich gezeigt, wie sich die
Advektionsgleichung auf die Grundform einer Evolutionsgleichung bringen läßt und wir werden einsehen, daß die Wirkung der Advektionsgleichung in der Verschiebung einer beliebigen
Anfangsfunktion besteht.
5.1.1 Die Advektionsgleichung als Evolutionsgleichung
Betrachtet man die Advektionsgleichung als Evolutionsgleichung (2.2), so ist der Operator durch
gegeben. Er besitzt die Eigenfunktionen , denn es gilt
0 womit man die Eigenwerte als 0 erkennt. Ohne es zu beweisen, sei ge-
sagt, daß die komplexen Exponentialfunktionen eine Orthogonalbasis des Lösungsraumes mit
kontinuierlichem Spektrum bilden. Da die Eigenwerte einen nichtnegativen Realteil (dieser ist
nämlich Null) haben, hat die Advektionsgleichung bzgl. jeder Norm eine dämpfende Wirkung
auf entsprechende Anfangswerte.
51
KAPITEL 5. DIE ADVEKTIONSGLEICHUNG
52
5.1.2 Differentialgleichungssysteme 1. Ordnung
Def.: Ein allgemeines quasilineares System partieller Differentialgleichungen 1. Ordnung der
Form
0 (5.2)
heißt hyperbolisch, falls die Matrizen N reelle Eigenwerte besitzen. Für jede Raumrichtung ist also eine Matrix zu untersuchen. Sind dagegen komplexe Eigenwerte vorhanden, nennt man das System elliptisch. Eine Differentialgleichung 1. Ordnung mit reellen Koeffizienten ist somit immer hyperbolisch.
Somit gehört die Advektionsgleichung ebenfalls zu den hyperbolischen Differentialgleichungen.
5.1.3 Bahnlinien als Charakteristiken der Advektionsgleichung
Die Besonderheit und Gemeinsamkeit hyperbolischer Differentialgleichungen sind ihre Charakteristiken. Dies sind Kurven oder Flächen in Raum und Zeit, auf denen die hyperbolische
partielle Differentialgleichung auf eine gewöhnliche reduziert werden kann.
Die lineare Advektionsgleichung läßt sich auf einer Bahnlinie, die durch
gegeben ist, als
(5.3)
schreiben. Letztere ist eine gewöhnliche Differentialgleichung. Somit sind die Bahnlinien der
Strömung die Charakteristiken der linearen Advektionsgleichung. Und diese besagt nichts anderes, als daß sich die physikalische Größe längs der Bahnlinien nicht ändert.
5.1.4 Das Anfangswertproblem
Als Anfangswertproblem bezeichnet man für eine zeitabhängige Differentialgleichung das
Problem, diese für vorgegebene Anfangsbedingungen auf einem beliebig großem Gebiet zu
lösen.
Die lineare Advektionsgleichung besitzt für das Anfangswertproblem
die allgemeine eindeutige Lösung
(5.4)
5.2. FD-VERFAHREN FÜR DIE ADVEKTIONSGLEICHUNG
53
Abbildung 5.1: Lösung des Anfangswertproblems für die lineare Advektionsgleichung.
Dies bedeutet, daß die lineare Advektionsgleichung die Anfangslösung in der Zeit um den
Vektor
verschiebt (Abb. 5.1). Die Form der Lösung wird dabei selbst nicht verändert. Damit haben wir
eine zweite Darstellung der dämpfenden Wirkung der linearen Advektionsgleichung auf die
Anfangslösung; ihre Extremwerte werden in keinem Fall vergrößert, sondern bleiben konstant.
Das Problem ist sachgemäß gestellt, kleine Änderungen in der Geschwindigkeit als auch an
der Anfangsbedingung bewirken nur kleine Änderungen der Lösung.
5.1.5 Das Randanfangswertproblem
Ist das Lösungsgebiet begrenzt, so muß neben der Anfangsbedingung am Einstromrand eine
Randbedingung in Form des Funktionswertes zu jedem Zeitpunkt vorgegeben werden. Am
Ausstromrand werden keine Randbedingungen benötigt.
Die physikalische Interpretation dieses Sachverhaltes ist einfach: Da die lineare Advektionsgleichung nur eine Verschiebung der Lösungsfunktion bewirkt, muß man vorgeben, was zu
jedem Zeitpunkt in das Lösungsgebiet hineingeschoben werden soll, das Hinausschieben der
Lösung aus dem Gebiet geschieht dann von selbst.
5.2 FD-Verfahren für die Advektionsgleichung
Die Lösung von zeitabhängigen Differentialgleichungen mit der Methode der Finiten Differenzen soll nun am Beispiel der eindimensionalen linearen Advektionsgleichung
konkretisiert werden. Als Testbeispiel untersuchen wir die Advektion einer normalisierten
Gaussverteilung mit der Standardabweichung - , die zur Anfangszeit die Form
KAPITEL 5. DIE ADVEKTIONSGLEICHUNG
54
haben soll. Die Lösung zu späteren Zeitpunkten bekommt man durch ein Verschieben der
Funktion um den Betrag auf der x-Achse:
Das Verhalten der jeweiligen numerischen Lösung wird dabei wieder mit dieser analytischen
Lösung verglichen, wofür man sich eine kleine Java-Applikation schreiben kann.
5.2.1 Explizite Verfahren
Das instabile FTCS-Verfahren
Unter Anwendung des Eulerverfahrens für die Zeitableitung das sogenannte Forward-TimeCentered-Space-Verfahren (FTCS)
(5.5)
wobei die Ortsableitung durch eine zentrale Differenz zum bekannten Zeitschritt approximiert wurde. Das Verfahren ist explizit, da man den Konzentrationswert am Knoten 0 zum
neuen Zeitschritt direkt als
bestimmen kann, wobei zum Zeitschritt bekannt sind.
Der Vergleich zwischen der analytischen und der numerischen Lösung in Abbildung 5.2 zeigt
schon nach 20 Zeitschritten starke Schwingungen, die den Maximalwert der analytischen
Lösung um das 80fache überragen, das numerische Verfahren wird also instabil.
Daher sollten wir untersuchen, wann das Verfahren stabil ist und schreiben es in der Form
wobei die Einheitsmatrix und
..
..
.
.
..
..
..
.
.
.
..
..
ist. Wegen
.
.
, kann das Stabilitätskriterium
nur für negative Zeitschritte erfüllt werden. Das Verfahren ist somit absolut instabil, wir sollten
es zur Lösung der linearen Advektionsgleichung nicht verwenden.
Nichtsdestotrotz besitzt es die Konsistenzordnung .
5.2. FD-VERFAHREN FÜR DIE ADVEKTIONSGLEICHUNG
55
Abbildung 5.2: Vergleich der mit dem instabilen FTCS-Verfahren gewonnenen numerischen
(rot) und der analytischen (grün) Lösung der Advektionsgleichung
Das Upstream-Verfahren
Beim Upstream-Verfahren wird die Ortsableitung mit Hilfe der stromaufwärts gelegenen Differenz gebildet, für ist dies:
(5.6)
(5.7)
Ist , so muß das Upstream-Verfahren als
gebildet werden. Die führenden Terme des Verfahrensfehlers sind
somit ist in beiden Fällen die Konsistenzordnung .
Die Randwerte werden beim Upstream-Verfahren automatisch richtig berücksichtigt: Ist die
Geschwindigkeit positiv, so muß am ersten d.h. Einstromknoten der Wert vorgegeben
werden. Alle weiteren Knoten werden dann mit dem Upstream-Verfahren berechnet.
Die Stabilität der Upsteam-Verfahrens ergibt sich aus dem allgemeinen Stabilitätskriterium,
wobei die Einheitsmatrix und KAPITEL 5. DIE ADVEKTIONSGLEICHUNG
56
Abbildung 5.3: Vergleich der mit dem Upstream-Verfahren gewonnenen numerischen (rot)
und der analytischen (grün) Lösung der Advektionsgleichung
Damit ergibt sich
und
..
..
..
.
..
.
.
.
..
..
.
.
woraus folgt, daß bei der Anwendung des Upstream-Verfahrens das sogenannte Courantkriterium
*/ (5.8)
erfüllt sein muß. Die soeben eingeführte dimensionslose Zahl */ bezeichnet man als Courantzahl.
Wir wollen das Verfahren an unserem Beispiel testen und nehmen die Einstellungen so vor,
daß die Courantzahl 0.5 ist. Abbildung 5.3 illustriert die Ergebnisse. Im Vergleich zur analytischen ist die numerische Lösung wesentlich breiter geworden und der Maximalwert ist
5.2. FD-VERFAHREN FÜR DIE ADVEKTIONSGLEICHUNG
57
Abbildung 5.4: Vergleich der mit dem Lax-Wendroff-Verfahren gewonnenen numerischen
(rot) und der analytischen (grün) Lösung der Advektionsgleichung
stark reduziert. Diesen Effekt bezeichnet man als numerische Diffusion, wir werden ihn später
analysieren.
Das Lax-Wendroff-Verfahren
Das Lax-Wendoff-Verfahren ist für die lineare Advektionsgleichung
Werden die Ortableitungen durch zentrale Differenzen diskretisiert, so ergibt sich
(5.9)
Das Lax-Wendroff-Verfahren für die lineare Advektionsgleichung besitzt die Konsistenzordnung und das Stabilitätskriterium */ .
Das Lax-Wendroff-Verfahren ist am ausstromseitigen Rand nicht anwendbar. Ist dieser der
letzte Knoten des Gebietes , so würde formal ein Wert für den Knoten benötigt, der
jedoch außerhalb des Lösungsgebiets liegt. Daher muß man am ausstromseitigen Rand auf das
Upstream-Verfahren ausweichen. Dies gilt grundsätzlich auch für alle anderen Verfahren, die
am Rand Knoten verwenden, die außerhalb des Gebietes liegen.
Der Test des Verfahrens (Abbildung 5.4) wird wieder mit */ durchgeführt. Er bestätigt
die Stabilität des Verfahrens, dank der höheren Genauigkeit des Verfahrens wachsen die Werte nicht ins Unendliche, obwohl zentrale Differenzen verwendet wurden. Es zeigt sich aber
wieder eine gewisse numerische Diffusion, das Peakmaximum ist reduziert und aufgeweitet.
Ferner bilden sich im Kielwasser des Primärpeaks Wellen, wodurch die numerische Lösung
KAPITEL 5. DIE ADVEKTIONSGLEICHUNG
58
auch negative Werte annimmt. Dies kann gravierende Folgen haben, wenn die Advektionsgleichung Teil eines simulierten Systems ist, wobei andere Prozesse schlimmstenfalls von der
Quadratwurzel der Konztentration abhängen. Dann würde dieser Prozess zum Abbruch des
Programms führen, wenn die Konzentration einen negativen Wert annehmen würde.
Wir bezeichnen die vom Lax-Wendroff-Verfahren nicht erfüllte Qualitätsanforderung als Positivität, sie verlangt, daß eine numerische Lösung positiv bleibt, wenn die reale Lösung immer
positiv ist und wenn Anfangs- und Randwerte ebenfalls positiv sind.
Das Leap-Frog-Verfahren
(5.10)
weist den Abbruchfehler
auf und besitzt somit die Konsistenzordnung und ist für */
stabil.
Das Leap-Frog-Verfahren benötigt die Abspeicherung der beiden Zeitschritte und und
ist dabei nicht genauer als das Lax-Wendroff-Verfahren. Weiterhin wird für die Berechnung
des ersten Zeitschrittes ein anderes Verfahren benötigt.
Das diffusive oder Lax-Verfahren
(5.11)
besitzt die Konsistenzordnung und ist ebenfalls für */ stabil. Wieder müssen
zwei Zeitschritte im Speicher vorgehalten werden.
5.2.2 Implizite Verfahren
Prinzipiell hat man auch bei impliziten Verfahren für die Advektionsterme die Auswahl zwischen Upwind- und zentralen Differenzen. Die Entscheidung fällt hier jedoch nicht schwer.
Wir erinnern uns, daß wir uns dem Upstreamverfahren einzig und allein aus Stabilitätsgründen zugewendet haben, es ist nämlich den zentralen Differnzen in der Genauigkeit um
eine Größenordnung unterlegen. Da implizite Verfahren aber nur zur Stabilisierung verwendet
werden, brauchen wir nicht eine Fliege mit zwei Klappen erschlagen: Wenn schon implizite
stabile Verfahren, dann auch zentrale Differenzen zweiter Konsistenzordnung für die Advektion.
Das voll implizite Verfahren
Setzen wir ein Verfahren in der Form
(5.12)
5.2. FD-VERFAHREN FÜR DIE ADVEKTIONSGLEICHUNG
59
an, so werden die Knotenwerte am neuen Zeitschritt untereinander gekoppelt, wodurch man
nicht mehr direkt nach einem gesuchten Knotenwert auflösen kann. Das Verfahren heißt deshalb implizit, denn es nötigt zur Lösung eines Gleichungssystems für alle Knoten.
Es läßt sich ohne Berücksichtigung von Randbedingungen in der Form
...
..
.
...
..
..
.
.
..
..
..
.
.
.
..
..
.
.
darstellen. Da die Matrix
..
..
.
.
..
..
..
.
.
.
... .. .. .. + ..
..
.
.
positiv-definit ist, , ist das Verfahren nach dem modifizierten Satz von Kellogs
unabhängig vom Zeitschritt stabil. Desweiteren ist es positiv, da
+ Der Abbruchfehler zum Zeitschritt ist:
Es besitzt also die Konsistenzordnung .
Das Crank-Nicolson-Verfahren
besitzt für die Konsistenzordnung und ist für stabil.
(5.13)
In Simulationsmodelle implementiert man gewöhnlicherweise das Crank-Nicolson-Verfahren
als implizites Verfahren, da über den Crank-Nicolson-Faktor gleichmäßig zwischen dem
instabilen expliziten und dem voll impliziten Verfahren variiert werden kann.
Das Preissmann-Verfahren
Das Preissmann-Verfahren verwendet zwei nicht unbedingt benachbarte Knoten
besitzt die Konsistenzordnung und ist für $ stabil.
(5.14)
0 und 1 . Es
KAPITEL 5. DIE ADVEKTIONSGLEICHUNG
60
5.2.3 Prädiktor-Korrektor-Verfahren
Prädiktor-Korrektor-Verfahren berechnen die Lösung in zwei Schritten: In einem ersten
Prädiktor-Schritt wird eine stabile Anfangslösung erzeugt, die in einem Korrektor-Schritt dann
bezüglich der Konsistenzordnung verbessert wird.
Das MacCormack-Verfahren
Das MacCormack-Verfahren [18] verwendet als Prädiktor das Upstream-Verfahren:
(5.15)
Besteht das Upstream-Verfahren aus Rückwärtsdifferenzen, so werden im Korrektor-Schritt
Vorwärtsdifferenzen (und umgekehrt) verwendet:
(5.16)
Aus beiden Schritten wird die neue Lösung gewichtet:
(5.17)
Das Verfahren besitzt die Konsistenzordnung und ist für */ stabil. Für die
lineare Advektionsgleichung ist es mit dem Lax-Wendroff-Verfahren identisch. Erst die Übertragung der Grundidee auf andere (nichtlineare) Gleichungen führt zu neuen Verfahren. Das
MacCormack-Verfahren ist hervorragend geeignet, Diskontinuitäten in der Lösung darzustellen.
Das Beam-Warming-Verfahren
Als Prädiktor kann wieder ein beliebiges explizites stabiles Verfahren verwendet werden. Die
so erhaltenen Werte seien wieder mit bezeichnet. Der Korrektor ist durch
6
(5.18)
gegeben. Die Zeitableitung wichtet die Zeitebenen und , die Ortableitung entsteht durch
Wichtung des Prädiktorergebnisses und des Zeitschrittes . Das Verfahren ist in Abhängigkeit
vom Prädiktor stabil, die Konsistenzordnung ist durch gegeben.
2 2 6
5.3 Numerische Diffusion
Im Verfahrensfehler taucht oftmals ein Term auf, der die zweite Ableitung in Ortsrichtung
enthält. Dies bedeutet, daß man eigentlich eine Gleichung löst, die einen Term umgekehrten
Vorzeichens enthält, der physikalisch eine Diffusion beschreibt.
Aber auch andere Terme mit zweiten Ableitungen nach der Zeit des Verfahrensfehlers sind
Quellen numerischer Diffusion, da sie sich durch die Gleichung
5.4. NUMERISCHE DISPERSION
61
in Diffusionsterme unformen lassen. Untersuchen wir daraufhin unsere Verfahrensfehler nochmals.
Numerische Diffusion des Upstream-Verfahrens
Sie ist größer Null, wenn die Courantzahl kleiner als eins ist. Für */ ist das Verfahren diffusionsfrei. Für Courantzahlen größer als eins ist der Diffusionsbeiwert negativ. Die Tatsache,
daß Gleichungen mit negativem Diffusionsbeiwert keine Lösungen besitzen, ist ein weiteres
Anzeichen für die Instabilität des Verfahrens, wenn das Courantkriterium nicht erfüllt ist.
Numerische Diffusion des impliziten Verfahrens
5.4 Numerische Dispersion
Für die Advekionsgleichung läßt sich neben der Konsistenz und Stabilität ein weiteres Gütekriterium formulieren. Dieses soll die Qualität des numerischen Verfahrens beschreiben, eine Anfangslösung auch tatsächlich mit der geforderten Advektionsgeschwindigkeit fortzubewegen. Wir werden sehen, daß dies leider nicht der Fall ist. Nehmen wir Wellen als Anfangslösungen, so werden wir zeigen, daß die numerische Advektionsgeschwindigkeit sogar
von der Wellenzahl abhängig ist. Da der Physiker das Phänomen, daß die Ausbreitungsgeschwindigkeit von Wellen von ihrer Frequenz oder Wellenzahl abhängig ist, als Dispersion
bezeichnet, wollen wir das entsprechende numerische Phänomen als numerische Dispersion
bezeichnen.
Da die komplexe Exponentialfunktion der Form
eine Lösung der Advektionsgleichung ist, und man nach dem Satz über komplexe Fourierreihen jede (normale) Funktion durch eine Reihe von komplexen Exponentialfunktionen darstellen kann, wollen wir schauen, ob eine analoge Funktion auch die mit zentralen Differenzen
diskretisierte Gleichung
löst. Den Ortsindex bezeichnen wir hier mit 1 , um ihn von der imaginären Einheit 0 zu unterscheiden. Im Gegensatz zur exakten Lösung führen wir die Advektionsgeschwindigkeit des
numerischen Verfahrens ein und hoffen natürlich, daß diese gleich oder wenigstens ähnlich der
exakten Advektionsgeschwindigkeit ist. Wir wagen also den Ansatz:
£
Einsetzen dieser Probelösung in die diskrete Gleichung
KAPITEL 5. DIE ADVEKTIONSGLEICHUNG
62
0 £
£
£
und Division durch £
0 liefert mit Hilfe der Identität Das Ergebnis ist entmutigend: Für 7 ist die numerische Advektionsgeschwindigkeit
Null. Da heißt dies, daß Wellen der Länge ' nicht transportiert werden,
diese Geister des numerischen Verfahrens werden manchmal als -Wellen bezeichnet. Die
numerische Dispersion der zentralen Differenzen ist also erheblich.
Führt man die gleiche Analyse für das Upstream-Verfahren durch, so lautet die Dispersionsbeziehung
0
0
Die Dispersionseigenschaften der Upstreamdifferenzen sind wesentlich besser als die der zentralen Differenzen, da die numerische Advektionsgeschwindigkeit wegen der Division durch
für im Verhältnis zur Gitterauflösung lange Wellen gegen die exakte Advektionsgeschwindigkeit konvergiert.
Kapitel 6
Interpolation auf Finiten Elementen
Als Interpolationaufgabe bezeichnet die Konstruktion von Funktionen, die durch eine gegebene Menge von Funktionswerten laufen soll. Der Interpolation liegen zwei duale Ideen zugrunde:
Man rekonstruiere eine Funktion aus möglichst wenigen Funktionswerten.
Man ersetze eine Funktion ohne Informationsverlust durch möglichst wenige Funktionswerte (Bsp.: zwei Punkte legen eine Gerade, drei Punkte einen Kreisbogen fest).
In beiden Fällen benötigt man Zusatzinformationen über die Funktion, etwa daß es sich um
eine lineare Funktion, einen Kreisbogen oder ein Polynom gegeben Grades handelt. Ohne
diese Information ist das Problem nicht sinnvoll gestellt, denn es gibt z.B. beliebig viele stetige
Funktionen, die durch endlich viele (also z.B. zwei) Punkte laufen. Eine stetige Funktion ist
durch die Angabe von endlich vielen Funktionswerten nicht eindeutig festgelegt.
Auch bei der numerischen Lösung von partiellen Differentialgleichungen kommt der Interpolation eine grundlegende Aufgabe zu: Wir ersetzen das Lösungsgebiet durch eine endliche
Anzahl von Punkten. Die Lösung der Differentialgleichung wird als Interpolationsfunktion
auf diesen Punkten gesucht. Damit wird einerseits das Lösungsgebiet durch eine endliche Anzahl von Punkten ersetzt und andererseits die Menge der zulässigen Lösungsfunktionen eingeschränkt, indem man nur noch die entsprechenden Interpolationsfunktionen als mögliche
Lösungen in Betracht zieht. Diesen Prozeß bezeichnet man als Diskretisierung, sie hat einen
geometrischen und einen analytischen Aspekt.
6.1 Lagrangesche Interpolationspolynome
Wir wollen die Lösung der partiellen Differentialgleichungen in Form eines Lagrangeschen
Interpolationspolynoms suchen. Diese werden durch die Werte auf diskreten Punkten
aufgespannt, die man auch als Knoten bezeichnet. Es gilt also:
0 (6.1)
Zwischen den Knoten wird für die Lösungsfunktion der Ansatz
(
63
(6.2)
64
KAPITEL 6. INTERPOLATION AUF FINITEN ELEMENTEN
gewählt. Die sogenannten Lagrangeschen Ansatzfunktionen ( sollen somit die Eigenschaft
( für
für
1 0
10
erfüllen, wodurch die Koeffizienten der Linearkombination die gesuchten Funktionswerte auf
den Knoten sind. Für eindimensionale Probleme lassen sich die Lagrangeschen Ansatzfunktionen einfach angeben, sie sind:
( 0
(6.3)
Die Lagrangeschen Interpolationspolynome sind normiert
(
(6.4)
d.h. die Summe aller Ansatzfunktionen ist an jedem Punkt des Raumes eins.
Lagrangesche Interpolationspolynome lassen sich auch für mehrdimensionale Probleme konstruieren, allerdings ist die Anzahl der Knoten im Gegensatz zum eindimensionalen Fall nicht
beliebig. Im n-dimensionalen Raum benötigt man zum vollständigen Aufspannen eines Polynomraumes der Ordnung k
0
Knoten.
Leider sind die Approximationseigenschaften der Lagrangepolynome sehr schlecht: Geht man
davon aus, daß die Funktionswerte auf den Knoten den Bereich der möglichen physikalischen Werte der Funktion abdecken, so sollte der Wertebereich der interpolierten Funktion
möglichst ähnlich sein: Stellt eine Strömungsgeschwindigkeit im Bereich kleiner 1 m/s
dar, dann sollte zwischen den Knoten nicht Werte von 100 m/s interpolieren. Genau das
passiert allerdings und die Tendenz zu diesem Verhalten ist umso größer, desto höher der Grad
des Interpolationspolynoms ist.
Selbst im asymptotischen Verhalten bei immer feiner werdender Belegung mit Knoten kann
es passieren, daß die Folge der Interpolationsfunktionen nicht gegen die zu interpolierende
Funktion konvergiert. Dies besagt der Satz von Faber.
6.2 Finite Elemente
Um das Strömungsgebiet einerseits durch hinreichend viele Knoten zu diskretisieren und andererseits den Grad der Ansatzfunktionen nicht zu hoch zu schrauben, kann man das Gebiet
in einzelne sich nicht überschneidende Teilgebiet zerlegen. Diesen Prozeß nennt man Partitionierung oder manchmal auch Triangulierung, weil man zweidimensionale Gebiete oftmals mit
Dreicksnetzen partitioniert.
6.3. EINDIMENSIONALE ELEMENTE
65
Die zur Partitionierung des Strömungsgebietes gewählten Elemente werden durch vier Eigenschaften gekennzeichnet:
Form der Elemente: Bei eindimensionalen Problemen sind nur Linienelemente vorstellbar, bei ebenen Problemen werden in der Regel Dreiecke oder Vierecke gewählt, bei
räumlichen Problemen eignen sich Prismenelemente oder Hexaeder.
Die Anzahl der Knoten: Die Knoten eines Elementes sind die Stützstellen für die Ansatzfunktion. Sie liegen meistens in den Ecken der finiten Elemente, es ist aber auch
möglich, weitere Knoten auf den Kanten oder im Inneren des Elementes zuzulassen.
Die Wahl der Knotenvariablen: Auf den Knoten werden entweder die Strömungsvariablen (Lagrangeknoten) oder eine ihrer Ableitungen (Hermiteknoten) gesucht.
Die Ordnung des Interpolationspolynoms: Die Interpolationsfunktion ist in der Regel
ein mehrdimensionales Polynom der Ordnung 0, 1 oder 2. Dabei muß die Anzahl der
Knoten pro Element mit der Anzahl der Koeffizienten der Interpolationsfunktion übereinstimmen.
Diese Eigenschaften zeigen, daß sich für die Konstruktion von Finiten Elementen eine Vielfalt
von Möglichkeiten bietet.
6.3 Eindimensionale Elemente
Eindimensionale Elemente sind einfach Strecken (auf der x-Achse) zwischen zwei Knoten
und , ein Element kann sich also über mehrere Knoten erstrecken. Jedes so definierte
Element läßt sich auf ein Einheitselement zwischen 8 und 8 durch
8
transformieren. Diese Transformation auf Einheitselemente ist eine wichtige Technik, die viel
Rechenarbeit erspart.
6.3.1 Lineare Interpolation
Ein Linienelement mit linearem Interpolationsansatz besitzt zwei Knoten auf seinen Rändern.
Auf dem Einheitselement ist der lineare Ansatz durch
8
8
8
8
gegeben. An jedem Punkt 8 des Elements gilt offensichtlich 8 8 . Eigentlich
suchen wir jedoch die Ansatzfunktionen auf einem beliebigen Linienelement der Nummer ,
welches ohne Einschränkung der Allgemeinheit aus den Knoten und bestehen soll. Mit
der eingeführten Einheitstransformation ergeben sich die Ansatzfunktionen als
( ( KAPITEL 6. INTERPOLATION AUF FINITEN ELEMENTEN
66
Für die lineare Interpolation einer Strömungsgröße auf dem Element ergibt sich somit:
Durch Differenzieren dieser Funktion erhalten wir zudem eine numerische Approximation für
die Ableitung auf dem Element :
In gleicher Weise ergibt sich für die Approximation der zweiten Ableitung Null. Wir folgern
daher messerscharf, daß sich zweite Ableitungen nur schlecht mit linearen Ansatzfunktionen
für die Ausgangsfunktion approximieren lassen.
Insgesamt läßt sich jedem Knoten global die Ansatzfunktion (siehe Abb. 6.2)
( für
für
für
und
(6.5)
zuordnen. Offensichtlich erfüllt diese Ansatzfunktion die Lagrangeschen Interpolationseigenschaften am Knoten eins und auf allen übrigen Null zu werden.
6.3.2 Der Interpolationsfehler
Wir kommen nun zu der wichtigen Frage, wie gut die Funktion approximiert. Dazu
müssen wir natürlich irgend etwas über die Funktion wissen bzw. annehmen. Wir nehmen
also an, daß der betragsmäßig größte Wert der zweiten Ableitung durch begrenzt ist:
Die maximale Gitterweite sei:
Zur Berechnung des Approximationsfehlers hat sich folgende Vorgehensweise bewährt: Bei
einer Interpolation n-ter Ordnung berechnet man zuerst den Fehler in der n-ten Ableitung.
Hier betrachten wir also die Abweichung in der ersten Ableitung für :
wird durch eine Taylorreihe mit hinreichend hohem Restglied ersetzt:
8 mit 8
. Damit folgt
6.3. EINDIMENSIONALE ELEMENTE
67
0
1
8
Abbildung 6.1: Allgemeines Linienelement und Einheitselement
(
1
(
1
Abbildung 6.2: Die einem Knoten sowie die einem Element zugeordneten Ansatzfunktionen
1
Abbildung 6.3: Die interpolierte Funktion .
KAPITEL 6. INTERPOLATION AUF FINITEN ELEMENTEN
68
8 8 Aus der Abschätzung
folgt schließlich
(6.6)
Diese Abschätzung ist wertvoll und wird später noch benötigt. Wir gewinnen aus ihr durch
Integration den gesuchten Fehler für die Funktion. Soll der Fehler am Punkt bestimmt werden, dann können wir so wählen, daß .
8 8 8 8 Damit folgt für den Interpolationsfehler die Abschätzung:
(6.7)
Bezeichnet man mit & die maximale Steigung von und mit 9 die maximale Variation
von auf einem Intervall , so findet man zusätzlich die Abschätzungen
9 (6.8)
&
(6.9)
6.3.3 Quadratische Interpolation
Zur Konstruktion quadratischer Interpolationspolynome benötigen wir drei Knoten. Das Einheitselement ist in Abb. 6.4 dargestellt.
Die darüber liegenden Ansatzfunktionen sind leicht nachprüfbar durch
8
8
8
8
8
8 8 8 8 6.3. EINDIMENSIONALE ELEMENTE
69
0
1/2
1
8
Abbildung 6.4: Einheitselement für quadratische Interpolationsfunktion
(
Abbildung 6.5: Quadratische Ansatzfunktionen auf einem Patch
gegeben, die Summe der drei Funktionen liefert an jedem Punkt den Wert eins.
Auf einem beliebigen Element mit den äquidistanten Knoten , und ergeben sich
die Ansatzfunktionen als
( ( ( Hier taucht der Ort des mittleren Knotens nicht mehr auf, weil die Knoten äquidistant
sind. Die approximierte Funktion wird auf diesem Element dargestellt als
Für die erste Ableitung erhalten wir die Approximation
KAPITEL 6. INTERPOLATION AUF FINITEN ELEMENTEN
70
und für die zweite Ableitung gilt:
6.3.4 Positivität und andere Extremalprinzipien
Die Lösungen elliptischer und parabolischer Differentialgleichungen unterliegen verschiedenen Extremalprinzipen, die den Wertebereich der Lösungsfunktion in irgendeiner Form beschränken. So sind die Lösungen der Transportgleichung bei entsprechenden Anfangs- und
Randbedingungen immer positiv. Das sollen sie auch sein, da sie auf der physikalischen Ebene Konzentrationsverteilungen darstellen sollen. Produziert ein diskretes Verfahren negative
Lösungen, so kann die zu beschreibende Physik erheblich durcheinander geraten. Wir könnten dann z.B. mit einer Reaktionskinetik negativer Konzentrationen konfrontiert werden. Das
wollen wir verhindern.
Dazu müssen wir untersuchen, welche Anforderungen die bei der Diskretisierung verwendeten
Ansatzfunktionen erfüllen müssen, damit solche Extremalprinzipien nicht verletzt werden. Wir
fassen das Ergebnis gleich in den folgenden
Satz: Seien die Lagrangeschen Ansatzfunktionen positiv und normiert. Dann werden die Extrema der interpolierten Funktion auf Knoten angenommen.
Bew.: Sei der kleinste Knotenwert. Dann gilt:
(
(
(
Für das Maximum versuche man diesen einfachen Beweis selbst.
Beide Bedingungen werden nur von den linearen Ansatzfunktionen erfüllt. Alle Ansatzfunktionen höherer Ordnung besitzen negative Bereiche. So würde man für die Knotenwertbelegung und bei quadratischen Ansatzfunktionen zwischen den Knoten
0 und 0 negative Werte interpolieren. Diese Situation tritt insbesondere an den Rändern
einer begrenzten Konzentrationswolke auf. Wir behalten diese Schwierigkeit im Hinterkopf.
6.4 Zweidimensionale Elemente
Aus der Vorgehensweise bei den eindimensionalen Elementen wurde offensichtlich, daß zur
Charakterisierung eines Elementtyps die Angabe der Ansatzfunktionen auf dem Einheitselement sowie die Einheitstransformation ausreicht. Daher werden hinfort nur diese angegeben.
6.4.1 Dreieck mit linearem Ansatz
In Abb. 6.6 sind die Bezeichnungen im allgemeinen und das zugehörige Einheitsdreieck dargestellt. Auf dem Einheitsdreieck sind die Lagrangeschen Ansatzfunktionen durch
6.4. ZWEIDIMENSIONALE ELEMENTE
71
:
1 1
8
Abbildung 6.6: Allgemeines Dreieck und Einheitsdreieck
Abbildung 6.7: Die zu einem Knoten gehörige Ansatzfunktion
8 : 8 : 8 : 8:
8
:
gegeben. Über die lineare Transformation
8 : 8 : 8 :
8 :
können diese auf jedes beliebige Dreieck abgebildet werden.
Die Umkehrung dieser Transformation bildet jedes Wertepaar
eck auf einen Punkt im Einheitsdreieck ab:
in einem beliebigen Drei-
KAPITEL 6. INTERPOLATION AUF FINITEN ELEMENTEN
72
8 : Setzt man diese in die Lagrangeschen Ansatzfunktionen auf dem Einheitsdreieck ein, dann
erhält man nach kurzweiliger Umformung die Lagrangeschen Ansatzfunktionen auf beliebigen
Dreiecken in der Form
( (6.10)
mit den Koeffizienten
;
; ; ;
;
; ; ;
;
; ;
;
(6.11)
und der Abkürzung:
; Die obige Darstellung der Gewichte der Ansatzfunktionen kann und sollte man mit dem Lagrangeschen Kriterium verifizieren. Tatsächlich wird die Funktion ( für eins und für
bzw. Null, entsprechendes kann man auch für die anderen beiden Ansatzfunktionen bestätigen. Die Gleichung
( ( ( liefert schließlich den interpolierten Wert an einem beliebigen Ort innerhalb des Dreieckes.
Ferner gilt für die ersten Ableitungen
und
6.4. ZWEIDIMENSIONALE ELEMENTE
73
:
1 1
8
Abbildung 6.8: Allgemeines lineares Viereck und Einheitsviereck
6.4.2 Viereck mit bilinearem Ansatz
Jedes allgemeine Viereck läßt sich mit Hilfe der Transformation
8 : 8 : 8 : 8:
8 : 8:
auf das in Abb. 6.8 dargestellte Einheitsviereck abbilden.
Dort sind die Ansatzfunktionen durch
8 : 8 : 8 : 8 : 8
:
8
:
8:
8 :
gegeben. Diese Ansatzfunktionen entstehen aus der Multiplikation der eindimensionalen linearen Ansatzfunktionen für beide Koordinatenrichtungen. Bilineare Ansatzfunktionen sind
somit keine linearen sondern quadratische Funktionen.
6.4.3 Viereck mit biquadratischem Ansatz
Die Knotentopologie für ein Viereck mit biquadratischem Ansatz ist in Abb. 6.9 dargestellt.
Die Ansatzfunktionen werden durch die Kombination der quadratischen Ansatzfunktionen des
Linienelementes gewonnen und sind:
KAPITEL 6. INTERPOLATION AUF FINITEN ELEMENTEN
74
7
5
6
4
8 9
3
2
1
Abbildung 6.9: Viereck mit quadratischem Ansatz
8 : 8 : 8 : 8 : 8 : 8 : 8 : 8 : 8 : 8
8
:
:
8 8 : : 8 8 : : 8 8 : : 8 8 : : 8 8 : : 8 8 : : 8 8 : : 8 8 : : 6.4.4 Sechseck mit quadratischem Ansatz
Für die vollständige quadratische Interpolation im zweidimensionalen Raum
benötigt man insgesamt sechs Knoten, um die sechs Gewichte eindeutig zu bestimmen. Die
quadratische Interpolation im Zweidimensionalen kann also entweder auf Sechsecken oder
wesentlich allgemeiner auf jeweils vier zusammenhängenden Dreiecken durchgeführt werden.
Wir wollen die Interpolation wieder als Lagrangeinterpolation
(
generieren. Die Lagrangeschen Ansatzfunktionen ( erfüllen die Eigenschaft
( und sind quadratische Polynome der Form
für
für
1 0
10
6.4. ZWEIDIMENSIONALE ELEMENTE
75
Abbildung 6.10: Quadratische Interpolation auf einem aus Dreiecken zusammengesetzten
Sechseck. Der Interpolationspunkt liegt im mittleren Dreieck.
( Die Interpolationsgewichte der sechs Ansatzfunktionen berechnet man aus den Lagrangebedingungen an den sechs Knoten:
Æ
Æ
Æ
Æ
Æ
Æ
Die Lösung des äquivalenten Gleichungsystems
Æ
Æ
Æ
Æ
Æ
Æ
kann man mit Hilfe eines direkten numerischen Verfahrens bestimmen.
In der Praxis wird man kaum Gitter mit sechseckigen Elementen verwenden. Hier ist die quadratische Interpolation auf Sechsecken deshalb wichtig, weil man grundsätzlich vier sich aneinanderschmiegende Dreiecke zu einem Sechseck zusammensetzen kann. Fällt also der In-
KAPITEL 6. INTERPOLATION AUF FINITEN ELEMENTEN
76
terpolationspunkt in ein bestimmtes Dreieck, dann nimmt man sich die an den Kanten angrenzenden Dreicke zu Hilfe, um sich ein Interpolationssechseck zu konstruieren.
Die 36 Lagrangeschen Gewichte werden pro Dreieck in der Initialisierungsphase der Anwendung bestimmt und können dann in der Arbeitsphase jeweils direkt in den Interpolationen
angewendet werden.
6.5 Interpolation im Dreidimensionalen
Die Möglichkeiten, dreidimensionale Elemente zu bilden, sind so vielfältig, daß wir uns verschreckt nur eines näher anschauen wollen. Das Lagrangesche Element der linearen Interpolation ist das Tetraeder, welches von vier Dreiecken begrenzt wird. Allerdings gibt es im Tetraedergitter keine sich über die gesamte Vertikale erstreckenden Kanten, so daß gerade die die
Strömung treibenden Gravitationskräfte schlecht dargestellt werden können. Dem Autoren der
Vorlesung ist in der Literatur kein einziges Modell für Strömungen in Oberflächengewässern
bekannt, welches mit Tetraedern arbeitet.
Als Alternative bieten sich Hexaeder, die von sechs Vierecken begrenzt werden sowie Prismen mit dreieckiger Grundfläche. Werden die Ecken des Hexaeders mit Knoten belegt, so
entsteht ein Element mit trilinearen Ansatzfunktionen, triquadratische Ansatzfunktionen enstehen durch zusätzliche Knoten auf den Kantenmitten und den Flächenschwerpunkten. Das
triquadratische Hexaederelement hat somit 27 Knoten.
6.5.1 Ein Prisma mit linearem Ansatz
Die Grundflächen des Prismas seien jeweils Dreiecke und die Säulen stehen vertikal. Ein solches Prisma ist in Abb. 6.11 mit dem zugehörigen Einheitsprisma dargestellt. Man beachte,
daß zur Abspeicherung eines solchen Prismas nur 12 Koordinatenangaben erforderlich sind.
Die Ansatzfunktionen im Einheitsprisma
8 : < 8 : < 8 : < 8 : < 8 : < 8 : < 8 :
<
<
:
<
8 :
<
8
<
8
<
:
können durch die Transformation
8 : 8 :
8 : 8 :
<
<
8 : 8 : 8 : 8 : 6.5. INTERPOLATION IM DREIDIMENSIONALEN
77
<
4
5
8
6
1
3
2
Abbildung 6.11: Geometrie und Bezeichnungen im Real- und Einheitsprisma.
auf jedes beliebige derartige Prisma abgebildet werden.
:
78
KAPITEL 6. INTERPOLATION AUF FINITEN ELEMENTEN
Kapitel 7
Die Diffusionsgleichung
Wir werden uns nun den sogenannten parabolischen Differentialgleichungen zuwenden. Ihr
Standardbeispiel ist die Diffusionsgleichung:
(7.1)
Neben ihrer eigenständigen Bedeutung zut Simulation des physikalischen Prozesses der Diffusion ist sie Teil anderer Differentialgleichungen, so z.B. der Navier-Stokes-Gleichungen zur
Beschreibung von Bewegungen Newtonscher Fluide oder der Transportgleichungen.
7.1 Die Diffusionsgleichung als Evolutionsgleichung
Betrachtet man die Diffusiongleichung als Evolutionsgleichung (2.2), so ist der Operator
durch
gegeben. Der Operator besitzt die Funktionen 0 und 0 als Eigenfunktionen, denn es gilt z.B.
0 0 0 womit man die Eigenwerte als erkennt. Ohne es zu beweisen, sei gesagt,
daß diese trigonometrischen Funktionen eine Orthogonalbasis des Lösungsraumes mit kontinuierlichem Spektrum bilden. Da die Eigenwerte alle positiv sind, hat die Diffusionsgleichung
bzgl. jeder Norm eine dämpfende Wirkung auf entsprechende Anfangswerte. Wir wissen also
nun schon, daß die Extremwerte der Anfangslösung zu späteren Zeitpunkten nie überschritten
werden.
7.2 Parabolische semilineare Differentialgleichungen
Das besondere Verhalten der Evolutionsgleichung beim Vorliegen nicht-negativer Eigenwerte
wird mit dem Begriff der parabolischen Differentialgleichung verallgemeinert. Hier wird aber
79
KAPITEL 7. DIE DIFFUSIONSGLEICHUNG
80
nur vorausgesetzt, daß alle Eigenwerte des Operators das gleiche Vorzeichen besitzen.
Wir starten mit der allgemeinen Form einer semilinearen Evolutionsgleichung zweiter Ordnung:
(7.2)
Die Gleichung kann also nichtlinear sein, da die Koeffizienten der ersten Ableitungen von
abhängig sein dürfen. Wenn wie hier nur die ersten Ableitungen nichtlinear sind, bezeichnet
man die Gleichung in der Literatur auch als semilinear, was ihre Nähe zu linearen Differentialgleichungen unterstreichen soll. Der Autor würde an dieser Stelle jedoch die Bezeichnung
seminichtlinear vorziehen.
Def.: Die Lösung der Differentialgleichung heißt parabolisch in , falls alle Eigenwerte
der Matrix reell sind und das gleiche Vorzeichen besitzen.
Parabolische Differentialgleichungen machen eine subtile Unterscheidung von dem, was
man unter ’vorwärts in der Zeit’ und ’rückwärts in der Zeit’ bewegen versteht. Dies sei am
Maximumprinzip für die Diffusionsgleichung präzisiert.
Satz: Minimum-Maximum-Prinzip. Genügt eine Funktion einer homogenen parabolischen Gleichung mit , und sind alle Eigenwerte der Matrix negativ, so nimmt
sie ihre positiven Maxima und negativen Minima für jedes + $ auf dem als parabolischen
Grenze bezeichneten Rand des Gebietes )= + an.
Im folgenden Widerspruchsbeweis wird eine Funktion konstruiert, die fast dieselben
Eigenschaften wie die Lösungsfunktion hat. An ihr ist die Widerspruchsführung allerdings
viel einfacher als an der Lösung. Ein klassisch schönes Beweisverfahren.
Bew.: Es sei
& und
! Da das Lösungsgebiet seinen Rand mit einschließt, gilt mit Sicherheit
& !
Wir nehmen nun fälschlicherweise an, daß sein Maximum im Inneren des L ösungsgebietes
an der Stelle annimmt. Somit gilt
&
Wir betrachten nun die Funktion
Sie hat folgende Eigenschaften:
Am Punkt
& ! +
nimmt sie denselben Wert wie an:
&
7.2. PARABOLISCHE SEMILINEARE DIFFERENTIALGLEICHUNGEN
Auf dem Rand ist sie kleiner als & , denn da & !
& !
!
81
+ gilt
& !
,&
Damit hat genauso wie sein Maximum in , jedoch nicht notwendig am Punkt
sondern an einem Punkt . Ein Maximum liegt dann vor, wenn
¼¼
¼¼
¼¼
sowie
Hess
¼¼
¼¼
¼¼
¼¼
¼¼
,
Sind alle Eigenwerte von negativ, folgt
und somit
¼¼
¼¼
¼¼
¼¼
¼¼
¼¼
Aus der Definitionsgleichung f ür folgt aber:
¼¼
¼¼
& ! +
¼¼
¼¼
¼¼
¼¼
¼¼
¼¼
¼¼
¼¼
¼¼
¼¼
¼¼
¼¼
& !
+
Widerspruch. Also nimmt sein Maximum auf dem parabolischen Rand des Gebietes an.
Es ist also entscheidend, daß die Lösung einer parabolischen Differentialgleichung für
$ gesucht wird. Das Problem ist für , nicht korrekt gestellt. Darum beschreiben
parabolische Differentialgleichungen irreversible Prozesse, wie sie etwa mit der Dissipation
von Energie oder Impuls verbunden sind.
Um die Parabolizität der Diffusionsgleichung zu erkennen, können wir die Zeit mit der Variablen , den Ort mit der Variablen identifizieren. Für die Koeffizienten in Gl. (7.2) gilt
somit
KAPITEL 7. DIE DIFFUSIONSGLEICHUNG
82
Die Eigenwerte einer diagonalen Matrix sind die Diagonalelemente, daher sind alle Eigenwerte von gleichen Vorzeichens, nämlich negativ.
7.3 Das Anfangswertproblem (1D)
Ist das Lösungsgebiet genügend groß, so besitzt die Diffusionsgleichung für das Anfangswertproblem
die allgemeine Lösung
7
exp
8
8 8
(7.3)
die man auch als Poissonformel bezeichnet. Zum Studium des physikalischen Verhaltens der
Diffusionsgleichung wollen wir jedoch eine Lösung betrachten, die mit Hilfe der Fourieranalyse gewonnen wird.
Aufgabe: Zeige, daß die Funktionen
und
allgemeine Lösungen der Diffusionsgleichung sind.
Die Aufgabe besagt, daß eine Anfangswelle bzw. Mode der Wellenlänge in der Zeit um
den Betrag gedämpft wird. Umso kleiner die Wellenlänge d.h. umso größer die Wellenzahl ist, desto schneller wird die Mode gedämpft.
Für beliebige Anfangsbedingungen läßt sich das Anfangswertproblem mit Hilfe des Satzes
von Fourier lösen:
Satz von Fourier: Sei eine reellwertige stückweise stetige Funktion mit stückweise stetiger
Ableitung. Weiter sei quadratintegrabel d.h.
, Dann läßt sich darstellen als
wobei
7.4. FD-VERFAHREN FÜR DIE DIFFUSIONSGLEICHUNG
83
Abbildung 7.1: Asymptotisches Verhalten einer Lösung der Diffusionsgleichung.
7
7
gilt.
Wenden wir den Satz von Fourier auf die Anfangsbedingungen an, so erhalten wir eindeutig
bestimmte und . Die Lösung des Anfangswertproblems ergibt sich dann zu:
(7.4)
Für wird nur die Mode nicht gedämpft, d.h. im asymptotischen Verhalten strebt
jede Lösung des Anfangswertproblems gegen den Mittelwert der Anfangsbedingung.
7.4 FD-Verfahren für die Diffusionsgleichung
Betrachten wir nun Finite-Differenzen-Verfahren für die Diffusionsgleichung
(7.5)
Damit repräsentiert die Geschwindigkeiten in den Navier-Stokes Gleichungen (mit Diffusion
ist hier somit die Viskosität gemeint) als auch Konzentrationen in den Transportgleichungen.
7.4.1 Das Forward-Time-Centered-Space-Verfahren (FTCS)
(7.6)
besitzt die Konsistenzordnung . Zur Untersuchung der Stabilität überführen wir
das Verfahren in Matrixschreibweise (ohne Berücksichtigung von Randbedingungen)
KAPITEL 7. DIE DIFFUSIONSGLEICHUNG
84
...
..
..
..
..
..
..
.
..
.
..
... .. .. .. .
..
.
..
.
.
.
Es verbleibt die Aufgabe, die Eigenwerte der Verfahrensmatrix zu bestimmen. Eine Matrix der
Form
..
..
..
.
..
.
..
.
..
.
..
.
.
.
besitzt die Eigenwerte
7
' wenn und reelle Zahlen gleichen Vorzeichens sind.
Aufgabe: Zeige mit Hilfe der Additionstheoreme für die trigonometrischen Funktionen, daß
die Matrix
..
..
..
.
..
.
..
.
.
.
..
..
.
.
die Eigenwerte
' 7
besitzt.
Zur Anwendung der Aufgabe können wir die Systemmatrix in der Form
+
)
darstellen, wobei ) die Einheitsmatrix ist. Wir erinnern uns daran, daß für die Berechnung der
Eigenwerte einer Verkettung bzw. Funktion von Matrizen der folgende Satz gilt: Wenn + die
Eigenwerte ' besitzt, dann besitzt + die Eigenwerte ' . Weiter gilt: Die Eigenwerte
einer Diagonalmatrix sind die Diagonalelemente.
Mit diesen beiden Sätzen bestimmen sich die Eigenwerte der Verfahrensmatrix als
7.4. FD-VERFAHREN FÜR DIE DIFFUSIONSGLEICHUNG
' 85
7
Die Stabilitätsbedingung lautet
7
7
die für
und somit für
(7.7)
erfüllt ist. Die neu eingeführte Stabilitätskennzahl bezeichnet man als Neumannzahl, auch
wenn dieser Name in der Literatur recht selten zu finden ist. Die Neumannzahl ist das Äquvalent zur Courantzahl der Advektionsgleichung, beide quantifizieren die Stabilität des expliziten
Verfahren.
In Abhängigkeit von der Diffusionskonstanten , der Gitterweite und dem Zeitschritt kann das explizite Verfahren für die Diffusionsgleichung also instabil werden. In diesem Falle
hilft nur eine Verkleinerung des Zeitschrittes.
7.4.2 Das Randanfangswertproblem
Ist das Lösungsgebiet begrenzt, so muß neben der Anfangsbedingung am Einstromrand eine
zu jedem Zeitpunkt an beiden
Randbedingung in Form der Ableitung des Funktionswertes Rändern vorgegeben werden. Eine Randbedingung in Form einer Ableitung nennt man auch
Neumannsche Randbedingung.
Es soll nun kurz darauf eingegangen werden, wie eine Neumansche Randbedingung in das
Verfahren eingebaut werden kann. Hierzu schreibt man die zweite Ortsableitung in der Form
Hier sehen wir die vordere und rückwärtige Differenz der Funktion erscheinen. Ist am Knoten
somit die Neumannsche Randbedingung in der Form
vorgegeben, so ergibt sich das Schema
> > In analoger Form versuchen wir, die Randbedingungen auch am anderen Rand und in andere
Schemata einzubauen.
KAPITEL 7. DIE DIFFUSIONSGLEICHUNG
86
7.4.3 Richardson- und DuFort-Frankel-Verfahren
Als Zeitschrittverfahren der zweiten Ordnung haben wir das Leap-Frog-Verfahren schätzen
gelernt, welches in der Zeitableitung den aktuellen Zeitschritt überspringt und diesen zur
Berechnung der sonstigen Terme, hier also des Diffusionstermes verwendet:
Das so entstandene Richardson-Verfahren besitzt die Konsistenzordnung und ist
leider generell instabil.
Im allgemeinen hilft es, Stabilitätsprobleme zu überwinden, wenn man mehr Rechengewicht auf den unbekannten Zeitschritt legt. So besteht die Möglichkeit, im RichardsonVerfahren den Term in zu zerlegen, womit das DuFort-Frankel-Verfahren
(7.8)
entsteht. Zwar ist es nun stabil, aber es ist inkonsistent. In unseren Versuchen, ein Verfahren 2.
Ordnung für die Diffusionsgleichung zu konstruieren, sind wir bisher gescheitert.
7.4.4 Das implizite Verfahren
(7.9)
besitzt die Konsistenzordnung . Unter Nichtbeachtung der Randbedingungen ent-
spricht dem impliziten Verfahren die Lösung des Gleichungssystems
..
..
.
.
..
..
..
.
..
.
.
Zur Stabiltätsanalyse reicht es, formal
...
..
..
..
...
..
..
..
..
.
.
..
..
.
.
..
..
..
.
..
.
..
.
.
.
... .. .. .. ... .. .. .. zu schreiben und man liest die Eigenwerte der Verfahrensmatrix leicht als
' ab. Da der Nenner immer größer als eins ist, ist das Verfahren unabhängig vom Zeitschritt
stabil.
7.5. BEWERTUNG
87
7.4.5 Das Crank-Nicolson-Verfahren
(7.10)
besitzt für die Konsistenzordnung und ist für stabil. Das Crank-
Nicolson-Verfahren verallgemeinert wieder das explizite (FTCS) und implizite Verfahren und
wird daher i.a. in Programmsysteme implementiert.
7.5 Bewertung
Für die Diffusionsgleichung ist die Auswahl eines Verfahrens sicherlich einfach. Sie wird aus
Konsistenz- und Stabilitätsgründen sicherlich auf das Crank-Nicolson-Verfahren fallen. In vielen Anwendungsfällen kann man aber auch eine Verletzung des Neumannkriteriums auschließen, bei denen dann auch das explizite Verfahren hinreichend genau ist.
88
KAPITEL 7. DIE DIFFUSIONSGLEICHUNG
Kapitel 8
Lagrange-Verfahren
Bis jetzt hatten wir in der Konstruktion der numerischen Verfahren nur die Eulersche Sichtweise eingenommen: Wir haben uns auf die Knoten des Gitters gesetzt und die zeitliche Entwick ausgewertet.
lung der Strömungsgrößen berechnet, d.h. wir haben die partielle Ableitung Nun wollen wir Verfahren kennenlernen, die sich an der Lagrangeschen Sichtweise der Hydrodynamik orientieren, d.h. wir versuchen, die zeitliche Änderung der Strömungsgrößen entlang
der Bahnlinien
(8.1)
zu analysieren. Daß sich diese der Strömung optimal angepaßte Sicht numerisch auszahlt,
sehen wir schon daran, daß sich die Anzahl der Terme in den Grundgleichungen erheblich
reduziert. So schreiben sich die dreidimensionalen Navier-Stokes-Gleichungen mit Hilfe der
totalen oder Lagrangeschen Ableitung als
9 (
9 (
>
89
(8.2)
KAPITEL 8. LAGRANGE-VERFAHREN
90
und die tiefenintegrierten Gleichungen von Saint-Venant:
! >
!
"
* 9 (
*
# " " "
!
>
*
!
"
* 9 (
# " " "
(8.3)
Zur Diskretisierung der Lagrangeschen Ableitung müssen der Anfangs- oder Basispunkt einer
Bahnlinie zur Zeit und ihr Endpunkt zur Zeit bekannt sein. Dazu sind die Differentialgleichungen der Bahnlinien
(8.4)
zu integrieren. Dabei gibt es prinzipiell zwei verschiedene Vorgehensweisen:
Die Anfangspunkte der Bahnlinien werden auf die Knoten des Gitters zum Zeitpunkt gelegt. Die Differentialgleichungen der Bahnlinien werden in Zeitrichtung vorwärts bis
zum Zeitpunkt integriert.
Die Anfangspunkte der Bahnlinien werden auf die Knoten des Gitters zum Zeitschritt
gelegt. Die Differentialgleichungen der Bahnlinien werden rückwärts in Zeitrichtung integriert; man bezeichnet diese Verfahren manchmal auch als inverse LagrangeVerfahren und ihnen werden wir im folgenden unsere ganze Aufmerksamkeit widmen.
8.1 Das Lagrange-Verfahren für die Advektionsgleichung
Die Änderung einer Größe auf einer Bahnlinie ist durch
; ; gegeben. Die lineare Advektionsgleichung
besagt also nichts anderes, als daß sich auf einer Bahnlinie, die durch
gegeben ist, nicht ändert. Folgender Algorithmus bietet sich als Lösungsverfahren für die Advektionsgleichung an:
8.2. NUMERISCHE DÄMPFUNG
91
Bahnlinie
Abbildung 8.1: Bahnlinie bei konstanter Geschwindigkeit.
Im ersten Schritt berechnen wir die Bahnlinie eines Teilchens, welches zum Zeitschritt am Knoten ankommt. Dieses Teilchen ist zum Zeitpunkt am Ort gestartet, welchen
wir als Basispunkt der Bahnlinie bezeichnen (Abb. 8.1).
Zur Berechnung des Basispunktes integrieren wir die Gleichung der Bahnlinie über den Zeitschritt:
womit man
und somit
(8.5)
erhält. Hat man zu jedem Knoten seinen Basispunkt berechnet, so brauchen wir nur noch den
Funktionswert auf den Basispunkten zu interpolieren, da der Basispunkt nur in den seltensten Fällen einen Gitterknoten trifft. Dieser interpolierte Wert entspricht der gesuchten Lösung
, da auf der Bahnlinie konstant ist.
Liegt der Basispunkt zwischen den Gitterknoten und , so ergibt eine lineare Interpolation
(8.6)
Entsprechend lassen sich Verfahren mit quadratischer und Interpolationen höherer Ordnung
konstruieren. Mit dieser Interpolation wollen wir uns weiter beschäftigen.
8.2 Numerische Dämpfung
Fällt die Basis der Bahnlinie direkt mit einem Knoten zusammen, ist keine Interpolation erforderlich. Damit sind Lagrangeverfahren für ganzzahlige Vielfache der Courantzahl frei von
KAPITEL 8. LAGRANGE-VERFAHREN
92
numerischer Dämpfung. Da dies jedoch (mit schärfer werdender Brille) selten der Fall ist,
führt die erforderliche Interpolation mit Hilfe der Nachbarknoten zu einer mehr oder weniger
hohen numerischen Dämpfung.
Wir führen die residuelle Courantzahl 2 ein, die den gebrochenen Rest der Courantzahl
bezeichne, es gilt:
2
Fällt die Basis der Bahnlinie auf einen Knoten, dann ist die residuelle Courantzahl also Null.
Da die lineare Advektionsgleichung eine monochromatische Welle
lediglich in der Form
verschiebt, muß ein ideales numerisches Verfahren für den Dämpfungsfaktor
'
die Eigenschaft ' erfüllen. J.D. McCalpin ([20]) hat diese Dämpfungsfaktoren für die
numerische Lösung der eindimensionalen linearen Advektionsgleichung mit Hilfe von Lagrangeverfahren untersucht. Das Subskript gibt dabei die jeweilige Ordnung der LagrangeInterpolation an der Basis der Charakteristik an:
' ' ' ' 2 2 2 2 2 2 2 2 2
2 2 2 2 Dabei steht für . Die Tatsache, daß die Dämpfungsfaktoren kleiner als Eins
sind, zeugt von der Zeitschrittstabilität der Lagrangeverfahren.
Für lange Wellen ( ,, lassen sich die Dämpfungsfaktoren mit
' '
' ' ) und kleine residuelle Courantzahlen
2 2 2 2 approximieren. Da das Wellenzahlverhältnis in die quadratische als auch kubische Interpolation mit der vierten Potenz eingehen, sind die Dämpfungseigenschaften der beiden Verfahren ähnlich.
Die zugehörigen numerischen Diffusivitäten können aus
8.3. MONOTONIE
93
' berechnet werden. Für lange Wellen folgen die numerischen Diffusivitäten:
22 2 2 2 2 2 2 2 2 2 Diese Gleichungen enthüllen Altbekanntes aus der Trickkiste des Modellierers:
Die numerische Dämpfung nimmt mit größerem Zeitschritt ab; oftmals der bequemste
Weg, das Problem zu umschiffen.
Dumm gelaufen, wenn der Zeitschritt zur Vermeidung numerischer Dämpfung so groß
gewählt werden mußte, daß die zu untersuchenden Prozesse nicht mehr aufgelöst werden
konnten. Dann hilft nur noch eine feinere räumliche Diskretisierung, d.h. die elendige
Arbeit einer neuen Netzgenerierung.
Die numerischen Diffusivitäten sind – bis auf den linearen Fall – vom Wellenzahlverhältnis abhängig, d.h. dispersiv. Die Dämpfung nimmt somit dramatisch mit der Wellenlänge ab; lediglich die lineare Interpolation läßt sich von langen Wellenlängen nicht
beeindrucken.
Abb. 8.2 zeigt die relativen d.h. auf bezogen numerischen Diffusivitäten für .
Damit ist die Bedingung der großen Wellenlänge zwar erheblich verletzt, wir können jedoch
den qualitativen Verlauf der numerischen Diffusivitäten für höhere als lineare Interpolationen
noch erkennen, da diese schon für verfahren sind somit hochgradig diffusiv.
unter fallen. Lineare Interpolations-
8.3 Monotonie
Die Lösungen partieller Differentialgleichungen unterliegen verschiedenen Extremalprinzipien, die den Wertebereich der Lösungsfunktion in irgendeiner Form beschränken. So sind die
Lösungen der Transportgleichung bei entsprechenden Anfangs- und Randbedingungen immer
positiv.
Ein Verfahren, welches die zu einer Differentialgleichung gehörigen Extremalprinzipien nicht
verletzt, heißt monoton.
Da die Advektionsgleichung nur als Verschiebungsoperator auf gegebene Anfangsbedingungen wirkt, bleiben die Funtkionswerte zu jedem späteren Zeitpunkt in den Schranken der
Funktionswerte der Anfangsbedingung. Mathematisch formuliert heißt dies, daß die Lösung
an einem Knoten zum Zeitschritt durch die Maxima und Minima der Werte an den
Interpolationsknoten zum Zeitschritt begrenzt sein sollte. Jene indizieren wir mit 1 :
KAPITEL 8. LAGRANGE-VERFAHREN
94
0 ,1 4
0 ,1 2
R e la tiv e D iffu s iv itä t
0 ,1
0 ,0 8
lin e a r
q u a d r a tis c h
k u b is c h
v ie r te O r d n u n g
0 ,0 6
0 ,0 4
0 ,0 2
0
0
0 ,1
0 ,2
0 ,3
0 ,4
0 ,5
0 ,6
0 ,7
0 ,8
0 ,9
1
R e s id u e lle C o u r a n tz a h l
Abbildung 8.2: Numerische Diffusion der Lagrangeverfahren mit verschiedenen Interpolationsschemen.
Lagrangeverfahren für die Advektionsgleichung sind dann monoton, wenn die Interpolation an
der Basis der Bahnlinie linear ist. Interpolationen höherer Ordnung sind grundsätzlich nichtmonoton, was sich in Form von Oszillationen in der Lösung äußert.
Bermejo und Staniforth [2] konstruieren ein quasimonotones Verfahren (besser: monotones
Quasiverfahren) dadurch, daß sogenannte Overshoots $ und Undershoots , durch ihre ’erlaubten’ Grenzen und ersetzt werden.
Priestley schlägt 1992 ein allgemeineres aber rechenaufwendigeres monotones und konservatives Verfahren vor, welches auf der Konvexkombination der Lösungen von monotonem
d.h.i.d.R. linearem und nichtmonotonen d.h. Verfahren höherer Ordnung basiert. Indiziert man
die Lösung des Verfahrens höherer Ordnung mit H und die des linearen Verfahrens mit L, so
setzt man die Lösung am Knoten als
2 $ 2 "
2 an. Der Parameter 2 sollte möglicht nahe an Eins unter Wahrung der Monotonie
2 $ 2 " gewählt werden. Der dabei gewonnene Freiheitsgrad (Wahl von 2 unter den genannten Restriktionen) kann schließlich so gewählt werden, daß die Masse im Gebiet erhalten bleibt.
8.4. BEISPIELE
95
8.4 Beispiele
Die folgenden Beispiele zeigen die mit dem Verfahren TRIM2D gewonnenen Ergebnisse. Ohne die Einzelfälle ausführlich zu diskutieren, werden jeweils das Lagrangeverfahren mit linearer und quadratischer Interpolation gegenübergestellt. Für den Transport wurde das beschriebene monotone Verfahren installiert, der Impuls wurde nicht monoton korrigiert.
8.4.1 Stufenfunktion
Abb. 2 zeigt die Ergebnisse für eine Durchbruchkurve einer Stufenfunktion. Obwohl sie durch
das An- und Ausschalten einer konstanten Randbedingung einfach erzeugt werden kann, ist
sie numerisch eine absolut harte Nuss: Die Fouriertransformierte einer Einheitsstufe am Ort
in der Darstellung
7
ist
Damit erstreckt sich ihr Frequenzspektrum periodisch über den gesamten Bereich der reellen Zahlen. Die hochfrequenten Anteile führen dementsprechend zu einer merklichen numerischen Diffusion in allen Verfahren, wobei diese im quadratischen Verfahren allerdings erheblich reduziert ist.
8.4.2 Hunte
Den dramatischen Einfluß der verschiedenen Verfahren auf die Tidewasserstände in der Hunte
am Pegel Oldenburg zeigt Abb. 8.4. Durch die engen Querschnitte bilden sich in der Hunte ausgeprägte Geschwindigkeitsprofile über die gesamte Breite aus. Die resultierenden Geschwindigkeitsgradienten lassen die Viskosität zu einem empfindlichen Parameter werden, der
die Dissipation der Tideenergie empfindlich beeinflußt. Die numerische Diffusion ist bei einer
örtlichen Diskretisierung von ! und einem Zeitschritt von mit etwa
! bei der linearen Interpolation sehr hoch. Im quadratischen Verfahren ist die numerische Diffusion für alle lunaren Partialtidewellen vernachlässigbar (Tabelle 8.1).
8.4.3 Jade
Bei ausgedehnteren Gebieten sind die horizontalen Geschwindigkeitsgradienten entsprechend
geringer, wodurch die physikalische und numerische Viskosität weniger Einfluß auf die Wasserstände nimmt. Dadurch sind die Änderungen beim Wechsel des Interpolationsschemas am
Pegel Voslapp im Modell der Jade entsprechend gering.
KAPITEL 8. LAGRANGE-VERFAHREN
96
Abbildung 8.3: Durchbruchkurven für die Stufenfunktion bei linearer und quadratischer Interpolation.
Partialtide Frequenz [Hz]
M2
M4
! M6
! M8
Wellenzahl [! ]
[! ] [! ]
1.25
1.25
!! 1.25
1.25
Tabelle 8.1: Numerische Viskositäten der Lagrangeverfahren erster und zweiter Ordnung für
die lunaren Flachwassertiden. Die Wellenzahlen wurden für eine Wassertiefe von 10 m, die
Gitterweite wurde als 10 m, der Zeitschritt als 10 s und die residuelle Courantzahl mit 0.5
angesetzt.
8.4. BEISPIELE
97
Abbildung 8.4: TRIM2D-Modell der Hunte: Gerechnete Pegel Oldenburg und Differenzen für
die Interpolationsverfahren.
98
KAPITEL 8. LAGRANGE-VERFAHREN
Abbildung 8.5: TRIM2D-Modell der Jade: Gerechneter Pegel Voslapp und Differenzen für die
Interpolationsverfahren.
Kapitel 9
Die Transportgleichung
Die Transportgleichung verbindet die Advektionsgleichung mit der Diffusionsgleichung:
(9.1)
Sie ist eine parabolische Differentialgleichung 2. Ordnung, die in erweiterter Form die verschiedensten Transportprozesse in der Strömungsmechanik beschreibt.
9.1 Das Anfangswertproblem (1D)
Durch die Variablentransformation
geht die Transportgleichung in die Diffusionsgleichung über:
Damit ist jede Lösung der Diffusionsgleichung unter der entsprechenden Transformation auch
Lösung der Transportgleichung. Insbesondere liefert die Poissonformel die folgende Lösung
des Anfangswertproblems für die Transportgleichung:
7
exp
8 8 8
Als weitere Lösung der Diffusionsgleichung hatten wir die Fourieranalyse der Anfangsbedingung kennengelernt. Durch die obige Transformation ergibt sich eine allgemeine Lösung der
Transportgleichung in der Form
(9.2)
Sie beschreibt die Verschiebung des durch die Fourierkoeffizienten dargestellten Anfangszustandes im Strömungsfeld unter gleichzeitiger Glättung durch die Diffusion.
99
KAPITEL 9. DIE TRANSPORTGLEICHUNG
100
Dieses gutmütige dämpfende Verhalten kann man auch bei einer Interpretation der Transportgleichung als Evolutionsgleichung herausanalysieren. Der Operator nach Gl. (2.2) ist:
5
grad div grad
Er besitzt die Eigenfunktionen . Die zugehörigen Eigenwerte sind 0 . Da der Realteil der Eigenwerte positiv ist, hat die Transportgleichung
dämpfenden Charakter: Die Extremwerte einer Anfangslösung werden auch zu späteren Zeitpunkten nicht überboten.
9.2 Das Randanfangswertproblem
Ist das Lösungsgebiet begrenzt, so muß am Einstromrand der Funktionswert selbst als Dirichletrandbedingung oder als Fluß, der in der Form
> (9.3)
vorgegeben wird. Dies entspricht dann einer gemischten Randbedingung, da hier sowohl der
Funktionswert als auch seine Ableitung auftauchen.
Am Ausstromrand muß die Ableitung der Funktion vorgegeben werden.
9.3 FD-Verfahren für die Transportgleichung
Wir entwickeln im folgenden Finite-Differenzen-Verfahren für die eindimensionale Transportgleichung:
(9.4)
Doch bevor wir die einzelnen Verfahren analysieren wollen, soll ein wichtiger Anwendungsfall
der Transportgleichung vorgestellt werden, der zudem eine analytische Lösung besitzt.
9.3.1 Die stossartige und punktförmige Einleitung
Wir haben ein Gewässer im Auge, dessen Strömungsfeld angenähert homogen und stationär
sei. Zum Zeitpunkt werde am dem punktförmigen Ort eine trockene Masse ! der
Dichte eingeleitet. Die Ökologie interessiert sich allerdings nicht primär für die eingeleitete
Masse, sondern eher für ihre auf das Trägerfluid Wasser bezogene Konzentration . Diese ist
zum Einleitungszeitpunkt aber unendlich und somit numerisch recht schwierig zu behandeln.
Dennoch wird diese Einleitung durch die Diffusion über ihre anfängliche Beschränktheit auf
einen Punkt auf immer größere Räume ausbreiten und durch die Advektion vom Einleitungsort
fortgetragen werden. Somit wird die zeitliche Entwicklung der Konzentrationsverteilung durch
eine Transportgleichung beschrieben werden können und man bestätige, daß das Problem die
analytische Lösung
9.3. FD-VERFAHREN FÜR DIE TRANSPORTGLEICHUNG
!
101
hat und offensichtlich divergiert auch hier die auf das Trägerfluid bezogene Konzentration .
9.3.2 Explizite Verfahren
Jede unserer bisherigen Grundlgleichungen war eine numerische Kennzahl zugeordnet. So
haben wir für die Advektionsgleichung die Courantzahl als das Verhältnis von diskretisiertem advektiven zum Zeitableitungsterm und für die Diffusionsgleichung die Neumannzahl als
das Verhältnis von diskretisiertem Diffusions- zum Zeitableitungsterm kennengelernt. Bei der
Transportgleichung sollte auch das Verhältnis von diskretisierten advektiven zu diffusiven Termen eine Rolle spielen, welches wir als Pecletzahl
3 (9.5)
bei der stationären Transportgleichung schon kennengelernt haben, an dieser Stelle aber nochmals in Erinnerung rufen wollen.
Das FTCS-Verfahren
Wir rekapitulieren, daß das FTCS-Verfahren für die Advektionsgleichung grundsätzlich instabil ist und ein explizites Verfahren für die Diffusionsgleichung bedingt stabil sein kann. Es
besteht daher eine gewisse Hoffnung, daß das FTCS-Verfahren
(9.6)
für die Transportgleichung stabil ist. Offensichtlich besitzt es die Konsistenzordnung
, aber nur wenn die Bedingung
3 */
eingehalten wird. Eingehende Analysen zeigen, daß das FTCS-Verfahren für
*/
stabil ist. Dieses Kriterium ist allerdings mit dem Pecletkriterium
3 unter einen Hut zu bringen.
In Abbildung 9.1 sind analytische und numerische Lösung der stoßartigen Einleitung für eine
entsprechende stabile Parameterkombination dargestellt. Hier bemerken wir, daß die numerische Lösung steiler als die analytische ist, also insgesamt weniger Diffusion aufweist. Diesen
Effekt hatten wir bisher nur in der umgekehrten Form kennengelernt, wir können ihn also als
negative numerische Diffusion bezeichnen.
KAPITEL 9. DIE TRANSPORTGLEICHUNG
102
Abbildung 9.1: Vergleich der mit dem FTCS-Verfahren gewonnenen numerischen (rot) und
der analytischen (grün) Lösung der Transportgleichung
Das explizite Upstream-Verfahren
(9.7)
wird für verwendet (ansonsten vordere Differenzen für die Advektion) und besitzt die
Konsistenzordnung falls
3 */
gilt. Für das explizite Verfahren gilt das Stabilitätskriterium
*/ Numerische Oszillationen treten beim Upwind-Verfahren nicht auf.
Unser Testbeispiel wird in Abbildung 9.2 präsentiert. Die numerische Lösung ist in der von der
Advektionsgleichung schon gewohnten Form gegenüber der analytischen leicht verschmiert,
weist also eine positive numerische Diffusion auf.
Das Lax-Wendroff-Verfahren
(9.8)
wird aus dem Lax-Wendroff-Verfahren der Advektionsgleichung gebildet, indem man hier
einfach den diffusiven Term mitberücksichtigt. Die Konsistenzordnung ist allerdings mit
in Ortsrichtung geringer, das Verfahren ist für
9.3. FD-VERFAHREN FÜR DIE TRANSPORTGLEICHUNG
103
Abbildung 9.2: Vergleich der mit dem Upstreamverfahren gewonnenen numerischen (rot) und
der analytischen (grün) Lösung der Transportgleichung
*/
stabil. Damit keine numerischen Oszillationen auftreten, muß die Pecletzahl wieder kleiner als
zwei sein.
Die Validierung des Verfahrens am Beispiel der stoßartigen Einleitung (Abbildung 9.3) zeigt
eine exzellente Übereinstimmung von analytischer und numerischer L’osung, so daß man sich
hier nur noch über die Stabilitätsrestriktionen zu beschweren braucht.
Das Leap-Frog-Verfahren
hat die Konsistenzordnung und ist für
*/ (9.9)
stabil. Damit keine numerischen Oszillationen auftreten, muß die Pecletzahl wieder kleiner als
zwei sein.
9.3.3 Implizite Verfahren
Wie bei der Advektionsgleichung werden wir auch hier nur zentrale Differenzen für den advektiven Term verwenden, da die Stabilität schon durch die rechenintensive implizite Zeitdiskretisierung erlangt werden soll.
KAPITEL 9. DIE TRANSPORTGLEICHUNG
104
Abbildung 9.3: Vergleich der mit dem Lax-Wendroff-Verfahren gewonnenen numerischen
(rot) und der analytischen (grün) Lösung der Transportgleichung
Transport einer steilen Front
Wir wollen die Qualität der impliziten Verfahren anahnd einer analytischen Lösung für den
Transport einer steilen Front diskutieren. Die Front advektiere von in Richtung der
positiven x-Achse. Damit sind die Anfangs- und Randbedingungen durch
gegeben. Dieser Testfall ist extrem schwierig, da für niedrige Diffusivitäten die Advektion einer unstetigen Stufe verbleibt.
Das Problem besitzt die analytische Lösung [22]
erfc
"
#
"
" erfc
#
lediglich
(9.10)
Im folgenden sollen einige implizite numerische Verfahren auf einem Gitter der Schrittweite
, der Geschwindigkeit , der Courant-Zahl */ und der Peclet-Zahl
3 mit der analytischen Lösung der steilen Front verglichen werden.
Das Crank-Nicolson-Verfahren
(9.11)
9.3. FD-VERFAHREN FÜR DIE TRANSPORTGLEICHUNG
105
1 .0 0
0 .7 5
0 .5 0
0 .2 5
0 .0 0
0 .0 0
0 .5 0
0 .2 5
1 .0 0
0 .7 5
Abbildung 9.4: Simulation einer steilen Front mit dem impliziten Verfahren mit zentralen Differenzen
1 .0 0
0 .7 5
0 .5 0
0 .2 5
0 .0 0
0 .0 0
0 .5 0
0 .2 5
0 .7 5
1 .0 0
Abbildung 9.5: Simulation einer steilen Front mit dem zentralen Crank-Nicolson-Verfahren
besitzt für die Konsistenzordnung und für gilt .
Das gesamte Verfahren ist für stabil, es treten numerische Oszillationen auf, wenn die
Pecletzahl größer als zwei ist. Dies zeigt Abbildung 9.4 am Beispiel der steilen Front für das
implizite zentrale Verfahren, d.h. der Crank-Nicolson-Faktor ist eins. Die Front ist weniger
verschmiert als bei der Anwendung der Upstream-Differenzen, die numerische Diffusion ist
aber immer unakzeptabel hoch. Es entstehen leichte Oszillationen, wenn die Front den rechten
Rand erreicht.
Die Wahl des Crank-Nicolson-Faktors 0.5 führt zu einer befriedigenderen Approximation
der Steilheit der Front (Abbildung 9.5), was auf die zweite Ordnung des Verfahrens zurückzuführen ist. Es treten allerdings wegen der Verletzung des Pecletkriteriums sehr starke Oszillationen auf.
KAPITEL 9. DIE TRANSPORTGLEICHUNG
106
1 .0 0
0 .7 5
0 .5 0
0 .2 5
0 .0 0
0 .5 0
0 .2 5
0 .0 0
0 .7 5
1 .0 0
Abbildung 9.6: Simulation einer steilen Front mit Defekt-Approximation
Verfahren mit Defekt-Approximation
Um das Verfahren zu verbessern, wird oftmals die Defekt-Approximation [9] vorgeschlagen.
Die Grundidee dabei ist, den Abbruchfehler eines Verfahrens mit umgekehrten Vorzeichen in
der Differentialgleichung zu berücksichtigen und so den Verfahrensfehler weiter zu reduzieren.
Dabei entstehen naturgemäß höhere Ableitungen, die durch einfache Drei-Knoten-Differenzen
nicht mehr wiedergegeben werden können. Wendet man die Defekt-Approximation auf das
gestellte Problem an, so wird dieses zu
Wird hierauf das Crank-Nicolson-Verfahren mit Upwind-Differenzen angewendet, so erhält
man das Schema:
$ %
(9.12)
Die Ergebnisse für unsere steile Front (Abb. 9.6) sind oszillationsfrei, beinhalten aber wieder
eine unbefriedigend hohe numerische Diffusion.
Klassische Finite-Differenzen-Verfahren sind also unbrauchbar, das gestellte Problem befriedigend zu lösen, da sie entweder starke Oszillationen oder eine hohe numerische Diffusion
aufweisen.
9.4. LAGRANGE-VERFAHREN FÜR DIE TRANSPORTGLEICHUNG
107
9.4 Lagrange-Verfahren für die Transportgleichung
Auch für die Transportgleichung läßt sich ein Lagrange-Verfahren konstruieren. Im Unterschied zur Advektionsgleichung ist zusätzlich zur Lagrangeschen Ableitung der Diffusionsterm zu behandeln:
und integrieren sie entlang der Bahnlinie zwischen den Zeitpunkten und :
Alles weitere hängt nun davon ab, wie man das Integral auf der rechten Seite approximiert,
wenden wir also die Trapezregel
an.
$
%
Die zweite Ableitung zum Zeitschritt am Knoten 0 ist:
Wir benötigen noch die zweite Ableitung an der Basis, die sich durch lineare Interpolation aus
den zweiten Ableitungen an den Nachbarknoten ergibt. Ist wieder der Knoten links und rechts der Basis der Bahnlinie, so erhält man eine Diskretisierung in der Form
Somit ergibt sich ein implizites Verfahren der Form
KAPITEL 9. DIE TRANSPORTGLEICHUNG
108
Abbildung 9.7: Bahnlinie verläßt Einstromrand.
Ein weniger rechenintensives explizites Verfahren läßt sich dadurch konstruieren, daß man nur
die Basiswerte zur Darstellung des Diffusionstermes verwendet:
(9.13)
So wie bei der linearen Advektionsgleichung bestimmen wir im ersten Schritt die Basis der
von jedem Knoten ausgehenden Bahnlinien. Im zweiten Schritt interpoliert man den Wert der
Funktion aus den bekannten Funktionswerten. Für den Funktionswert an der Basis ergab
sich mittels linearer Interpolation
falls der Basisknoten zwischen den Knoten und liegt.
Etwas kompliziert wird die Situation, wenn eine Bahnlinie den Einstromrand verläßt (Abb.
9.7). Dazu sei angenommen, daß der Ausgangsknoten im Abstand ? vom Rand liegt. Wir
führen einen reduzierten Zeitschritt in der Form
? */
also
?
*/
ein und interpolieren die Randbedingung auf den Durchstoßpunkt der Bahnlinie:
An einem solchen Knoten rechnen wir dann mit reduziertem Zeitschritt weiter, so wird das
explizite Verfahren hier zu
9.5. BEWERTUNG
109
Abbildung 9.8: Vergleich der mit dem Lagrange-Verfahren gewonnenen numerischen (rot) und
der analytischen (grün) Lösung der Transportgleichung
Die Lagrange-Verfahren zeichnen sich wieder durch ihre optimalen Stabilitätseigenschaften
aus, es ist für
stabil. Numerische Oszillationen treten nicht auf.
Die Qualität des Verfahrens sei durch einen Blick auf Abbildung 9.8 demonstriert. Die besseren Stabilitätseigenschaften lassen hier eine Erhöhung des Zeitschrittes zu.
9.5 Bewertung
Wenn die teilsweise sehr restriktiven Stabilit”atskriterien in speziellen Anwendungen kein
Problem sind, dann liefern unter den expliziten Verfahren das Lax-Wendroff und das LeapFrog-Verfahren sehr gute Ergebnisse. Bei den implizten Verfahren sticht das Crank-NicolsonVerfahren durch seine Konsistenz heraus, versagt aber bei hohen Pecletzahlen. Diese treten in
der Praxis aber nur bei sehr steilen Fronten auf. Die Lagrangeverfahren bestechen bezüglich
ihrer Stabilität, die nur durch die Diffusion begrenzt ist und durch ihre geringe numerische
Diffusion, wenn man ein Interpolationsverfahren zweiter Ordnung anwendet. Ihr Nachteil ist
allerdings die mangelnde Massenerhaltung.
110
KAPITEL 9. DIE TRANSPORTGLEICHUNG
Kapitel 10
Die Burgersgleichung
Die physikalische Größe Impuls wird ebenfalls mit der Strömung transportiert. Wir können
also eine Impulsgleichung der eindimensionalen Strömungsmechanik gewinnen, indem wir
die Transportgleichung
für die Impulsdichte betrachten. Im Falle einer konstanten Fluiddichte erhält man:
(10.1)
Eine leicht andere aber mathematisch äquivalente Form dieser Gleichung wurde zum ersten
Mal von J.M. Burgers 1948 [4] untersucht und zur Illustration der Turbulenztheorie verwendet.
Ihm zu Ehren bezeichnet man obige Gleichung als Burgersgleichung.
Die Burgersgleichung ist noch nicht die Impulsgleichung der eindimensionalen Strömungsmechanik, da die zeitliche Änderung des Impulses durch Krafteinwirkung nicht enthalten ist.
Die rechte Seite der Gleichung beschreibt die Impulsdiffusion, Impuls besitzt also auch die
Tendenz, sich gleichmäßig über das Fluid zu verteilen. Den Koeffizienten bezeichnet man
als Viskosität.
Im Gegensatz zu den bisher behandelten Gleichungen ist diese Gleichung und somit auch die
Navier-Stokes-Gleichungen nichtlinear. Dabei bezeichnet man eine Gleichung als linear, wenn
jede Linearkombination von Lösungen wieder eine Lösung ergibt. Bei der Konstruktion von
Fourierlösungen für die Diffusionsgleichung hatten wir diesen Sachverhalt stillschweigend
verwendet. Sind und zwei verschiedene Lösungen der Burgersgleichung, d.h. es gilt
und
so löst die Summe der beiden Lösungen eben nicht die Burgersgleichung. Der Leser
überzeuge sich davon.
Der Sachverhalt der Nichtlinearität hat tragische Konsequenzen, was schon bei der Lösung des
Anfangswertproblems deutlich wird. Durch die Variablentransformation
111
KAPITEL 10. DIE BURGERSGLEICHUNG
112
geht die Burgersgleichung nicht mehr in die Diffusionsgleichung über.
10.1 Der Ansatz von E. Hopf
Nachdem Burgers die Gleichung an das Licht der mathematischen Öffentlichkeit geholt hatte, war schnell der Ehrgeiz erweckt, doch eine Transformation zur Diffusionsgleichung und
somit explizite Lösungswege zu finden. Diese Meisterleistung wurde von E. Hopf 1950 [?]
veröffentlicht.
Man schreibe dazu die Burgersgleichung in der Form
und führe die Transformation
(
Ê
%
(10.2)
durch. Die Inverse ist
(
(
( und für die benötigten Ableitungen folgt
(
$ & %
( (
(
(
(
(
(
womit sich
ergibt. Die Burgersgleichung erhält so die Gestalt
$ & %
&
(
$ & %
(
(
welche sich durch Integration bzgl. (wobei die Integrationskonstante vernachlässigt werde)
in die Diffusionsgleichung verwandelt:
(
(
Die Anwendung der Poissonformel liefert
10.2. DER ANSATZ VON J.D. COLE
( 113
7
exp
8
( 8 8
Nach Ausführung aller Rücktransformationen erhält man schließlich die explizite Lösung des
Anfangswertproblems für die Burgersgleichung:
'
Ê
%
Ê
%
8
8
Auch dieses Gleichungsmonster hilft nicht so recht weiter, allerdings kann man mit ihm zeigen, daß eine Lösung des Anfangswertproblems existiert und eindeutig ist, falls
stetig ist
und
8 8 für
Die Lösung selbst und ihre erste und zweite Ableitung sind dann sogar stetig. Beide Bedingungen sagen etwas über die notwendigen Bedingungen der Anfangsbedingung, sie soll stetig
sein und die Gesamtgeschwindigkeit - was so etwas wie das Integral über die Geschwindigkeit
ist - soll in den beiden Halbräumen und verschwinden. Beispiele für solche Anfangsbedingungen wären etwa die trigonometrischen Funktionen, da ihr Integral Null ist, nicht
aber die konstante Funktion 4, da diese im Integral eben nicht verschwindet.
Somit suchen wir bei einem numerischen Lösungsverfahren für die Burgersgleichung nicht
nach etwas, was nicht existiert.
10.2 Der Ansatz von J.D. Cole
Ein Jahr später veröffentlichte J.D. Cole folgenden Lösungsweg. Im ersten Schritt wird die
Transformation
(
durchgeführt. Eine ähnliche Transformation ist dem Leser vielleicht aus der Theorie der Potentialströmungen bekannt, dort geht man davon aus, daß sich sogar Geschwindigkeitsvektoren
durch den Gradienten einer skalaren Funktion darstellen lassen. Die Burgersgleichung erhält
nun die Gestalt
(
( ( (
Wieder wird (ohne Integrationskonstante) nach integriert
KAPITEL 10. DIE BURGERSGLEICHUNG
114
( (
(
und eine Lösung der Form
( > versucht. Dies führt auf die gewöhnliche Differentialgleichung
> >
welche sich über die Separation der Variablen in die Gleichungen
>
¼¼
> überführen läßt. Für lassen sich die beiden Gleichungen recht einfach lösen und man
erhält nach Ausführung aller Rücktransformationen eine spezielle Lösung der Burgersgleichung in der Form
Für (10.3)
ergibt sich insbesondere die Lösung
(10.4)
Somit haben wir nicht das Anfangswertproblem gelöst, aber zumindest eine partikuläre
Lösung der Burgersgleichung gefunden. Leider ist diese jedoch recht unschön und sogar unphysikalisch, da sie auf dem Lösungsgebiet immer Pole besitzt, an denen die Funktion unendlich wird und somit nicht numerisch auswertbar ist. Weitere 35 exakte Lösungen der Burgersgleichung haben Benton und Platzman [1] zusammengestellt.
10.3 FD-Verfahren für die Burgersgleichung
Das Lösungsgebiet liege ohne Einschränkung der Allgemeinheit zwischen
werde äquidistant durch die Knoten 0 diskretisiert.
und und
Das Upstream-Verfahren
Im Prinzip lassen sich alle Verfahren für die Transportgleichung auch auf die Burgersgleichung
anwenden, wenn man die Advektionsgeschwindigkeit zum Zeitschritt ansetzt. So erhält man
das explizite Upstream-Verfahren für die inneren Knoten als
für
für
,
(10.5)
10.3. FD-VERFAHREN FÜR DIE BURGERSGLEICHUNG
115
Dabei entstehen natürlich dynamische Stabilitätskriterien, da die Lösung selbst in die Courantzahl eingeht.
Das Crank-Nicolson-Verfahren
In derselben Weise lassen sich auch implizite Verfahren leicht übertragen, wie etwa das CrankNicolson-Verfahren mit zentralen Differenzen
(10.6)
welches wieder die Lösung eines linearen Gleichungssystems erfordert. Dies wird dann anders, wenn man eine korrekte Wichtung beiden Zeitschritten neuem und alten Zeitschritt ansetzt, d.h. auch die Advektionsgeschwindigkeit zum unbekannten Zeitschritt auswertet:
(10.7)
Hier entsteht ein nichtlineares Gleichungssystem, welches mit sehr rechenaufwendigen Newtonverfahren gelöst werden muß. Um dies zu vermeiden, werden iterative oder PrädiktorKorrektor-Strategien angewendet.
Das MacCormack-Verfahren
Das MacCormack-Verfahren verwendet als Prädiktor das Upstream-Verfahren:
(10.8)
Im Korrektor-Schritt wird die entgegengesetzte Differenzenapproximation verwendet
(10.9)
und schließlich wird die Lösung der beiden Schritte gewichtet:
Das Verfahren ist für */ stabil.
(10.10)
KAPITEL 10. DIE BURGERSGLEICHUNG
116
10.4 Lagrange-Verfahren für die Burgersgleichung
Genauso wie für die Transportgleichung als ihr lineares Analogon lassen sich LagrangeVerfahren auch für die Burgersgleichung konstruieren. In der Lagrangeschen Schreibweise
lautet die Gleichung:
Sie besagt, daß sich die Geschwindigkeit entlang einer Bahnlinie
nur durch die Wirkung der Viskosität ändert. Wir wollen nun den Basispunkt der Bahnlinie
bestimmen und integrierten diese Gleichung zwischen den Zeitschritten und :
Im Unterschied zur Transportgleichung ist auf der rechten Seite nun keine Konstante mehr.
Wenden wir also die Trapezregel an, so wird die Bestimmungsgleichung für den Basispunkt:
Dumm nur, daß man weder noch den Wert an der Basis kennt, weil wir ihren Ort ja erst
bestimmen wollen. Aus diesem Dilemma entzieht man sich in der Regel, indem man behauptet, was natürlich umso weniger richtig ist, desto größer der Zeitschritt ist.
Somit berechnet sich der Basispunkt in erster Näherung als
Nun kann durch Interpolation bestimmt und die Burgersgleichung in der Lagrangeschen
Form in vollkommener Analogie zur Transportgleichung gelöst werden.
10.5 Harmonische Anfangsbedingungen
Wir wollen nun das Verhalten der Lösungen der Burgersgleichung studieren, wenn die Anfangsbedingungen eine einfache harmonische Funktion der Form
7
'
ist. Damit diese Funktion ordentlich durch die Diskretisierung aufgelöst wird, sollte der Ortsschritt wesentlich kleiner als die Wellenlänge ' sein, also etwa
, '
Am Einstromrand wollen wir periodische Randbedingungen in der Form
10.5. HARMONISCHE ANFANGSBEDINGUNGEN
117
Abbildung 10.1: Lösung der Burgersgleichung mit dem Upstream-Verfahren für harmonische
Anfangsbedingungen.
9 vorgeben, am Ausstromrand setzen wir die Burgersgleichung ohne Diffusion in der UpwindForm an:
Beginnen wir mit dem Upstream-Verfahren. Das Ergebnis ist in gewohnter Form in Abbildung
10.1 gezeigt. Die anfängliche Sinuswelle hat sich sägezahnförmig aufgesteilt, so daß die Geschwindigkeit an einigen Stellen sprunghaft von einem positiven auf einen negativen Wert
übergeht. Darauffolgend steigt die Geschwindigkeit wieder stetig an.
Diese Lösung ist nicht mit der in Abbildung 10.2 dargestellten zu vergleichen, da dort für das
FTCS-Verfahren die Viskosität auf ! erhöht wurde, um das Stabilitätskriterium
*/
zu erfüllen. Es bestätigt sich auch in den anderen Verfahren, daß die Stabilitätskriterien für
die Transportgleichung sich direkt auf die Burgersgleichung übertragen lassen, wenn man nur
die entsprechende dynamische Courantzahl mit der aktuellen Strömungsgeschwindigkeit bzw.
deren Obergrenzen ansetzt.
Ein Verfahren zweiter Ordnung, welches dieselben Stabilitätskriterien wie das Upstreamverfahren erfüllen muß, ist das MacCormack-Verfahren (siehe Abbildung 10.3). Die größere Genauigkeit zeichnet sich im Vergleich zum Upstream-Verfahren durch noch steilere Sprünge
aus.
118
KAPITEL 10. DIE BURGERSGLEICHUNG
Abbildung 10.2: Lösung der Burgersgleichung mit dem FTCS-Verfahren für harmonische Anfangsbedingungen.
Abbildung 10.3: Lösung der Burgersgleichung mit dem MacCormack-Verfahren für harmonische Anfangsbedingungen.
10.6. ZUSAMMENFASSUNG
119
Abbildung 10.4: Lösung der Burgersgleichung mit dem Lagrange-Verfahren für harmonische
Anfangsbedingungen.
Schließlich sei noch das Lagrangeverfahren getestet, auch hier sieht die Lösung den anderen
sehr ähnlich.
10.6 Zusammenfassung
Alle numerischen Verfahren für die Transportgleichung lassen sich direkt auf die Burgersgleichung übertragen. Dies gilt auch für die Stabilitätskriterien und für die Konsistenzordnung. Für die Burgersgleichung haben wir keine analytischen Lösungen gefunden, die sich
für eine Verifikation der numerischen Verfahren eignen. Die numerischen Lösungen zeigen eine Verselbständigung des Inneren des Lösungsgebietes, hier werden (nicht nur harmonische)
Anfangsbedingungen in Sägezahnfunktionen umgeformt. Die Randbedingungen haben kaum
Einfluß auf die Lösung im Inneren.
Interpretiert man die Sägezahnlösungen der Burgersgleichung als Geschwindigkeiten, so produziert sie recht unrealistische Strömungsfelder. An den Sprungpunkten treffen hohe positive
direkt auf hohe negative Geschwindigkeiten. Wäre mit diesen Strömungen ein Materialtransport verbunden, so würde sich an den Sprungstellen sehr viel Material akkumulieren, während
es an den Nullstellen der Lösung langsam verschwindet, da an diesen in beide Richtung wegtransportiert wird.
Ein Fluid, welches durch die Burgersgleichung beschrieben wird, neigt somit zur Verklumpung.
120
KAPITEL 10. DIE BURGERSGLEICHUNG
Kapitel 11
Finite-Volumen-Verfahren
Finite-Volumen-Verfahren basieren auf der Bilanzierung von Flüssen auf Kontrollvolumina.
Damit ist die Finite-Volumen-Methode ein Verfahren, welches sich direkt aus den kontinuumsmechanischen Eigenschaften einer Strömung ergibt.
Finite-Volumen-Verfahren sind eine Verallgemeinerung der Finiten-Differenzen-Verfahren:
Auf einfachen Geometrien werden sich altbekannte Schemata ergeben. Die Denkweise der
Methode der Finiten Volumen d.h. die Bilanzierung von Flüssen auf Kontrollvolumina,
ermöglicht jedoch auch die systematische Ableitung neuer Schemata.
11.1 Die Divergenzform der Grundgleichungen
Alle Grundgleichungen der Hydrodynamik lassen sich in der sogenannten Divergenzform
5
@
div # (11.1)
darstellen. Tun wir dies für den vollständigen Satz der dreidimensionalen inkompressiblen
Navier-Stokes-Gleichungen im Schwerefeld einschließlich Kontinuitäts- und Transportgleichung, so ergibt sich die in Tabelle 11.1 dargestellte Termbelegung.
Auch für die tiefenintegrierten 2D- und die eindimensionalen Gleichungen von Oberflächenströmungen lassen sich Divergenzformen finden, sie sind in Tabelle 11.2 bzw. 11.3 dargestellt.
11.2 Diskretisierung
Zur Lösung der Grundgleichungen wird das Strömungsgebiet vollständig mit einzelnen Kontrollvolumen überdeckt. Jedes Kontrollvolumen wird durch einen Knoten 0 in seinem Inneren repräsentiert.
Die Grundgleichungen werden in ihrer Divergenzform über die Kontrollvolumina integriert:
div #5 121
@ KAPITEL 11. FINITE-VOLUMEN-VERFAHREN
122
5
@
#
x-Impuls
5 grad 9 (
y-Impuls
9 (
z-Impuls
5 grad Kontinuität
Transport
5 grad 5
5 grad >
0
0
Tabelle 11.1: Die dreidimensionalen Grundgleichungen der Hydrodynamik in der Divergenzform.
11.2. DISKRETISIERUNG
123
5
x-Impuls
y-Impuls
@
#
5 grad 5 grad !
>
Transport
*
!
"
* 9 (
# " " "
>
Kontinuität
!
*
!
"
*# " " " 9 (
5
0
5 grad 0
Tabelle 11.2: Die tiefenintegrierten zweidimensionalen Gleichungen der Gewässerströmungen
in der Divergenzform.
5
Impuls
Kontinuität
Transport
@
#
!
>
"
* *
# " "
Tabelle 11.3: Die eindimensionalen Gleichungen von Saint-Venant in der Divergenzform.
KAPITEL 11. FINITE-VOLUMEN-VERFAHREN
124
1
n
1
w
e
s
1
0
0
0
Abbildung 11.1: Zweidimensionale äquidistante Finite-Volumen-Geometrie.
11.2. DISKRETISIERUNG
125
Der Divergenzterm wird unter Anwendung des Gaußschen Integralsatzes integriert. Wir erinnern uns, daß der Gaußsche Integralsatz das Volumen- auf ein Randintegral unter Wegfallen
der Divergenz überführt:
&
5 A5 #
@ (11.2)
In diesem Schritt liegt das Wesen der Methode der Finiten Volumen. Die Gleichung besagt,
daß die zeitliche Änderung der Größe im Volumen durch die Bilanz der Flüsse #5 über
den Rand (dies ist umgangssprachlich nichts anderes, als die Integration über den Rand) und
durch Quellen bzw. Senken @ im Volumen stattfindet.
Durch die Integration des Divergenztermes fallen in den Grundgleichungen Ableitungen weg,
insbesondere werden die zweiten Ableitungen in den Diffusionstermen zu ersten Ableitungen. Somit kann man die Lösungen der Grundgleichungen in der klassischen Form nur in der
Menge der zweimal differenzierbaren Funktionen suchen. In der Integralform reicht einfache
Differenzierbarkeit, die Lösungsmenge der Funktionen ist damit viel größer und die Integralform somit viel allgemeiner gültig.
Rein technisch verbleibt nur noch die numerische Berechnung der drei Integralterme, die wir
nun einzeln angehen wollen.
11.2.1 Die Behandlung der Zeitableitung
Ist B das Volumen des Kontrollvolumens , so läßt sich die Integration nach Vertauschen von
Differentiation und Integration und Anwendung eines Eulerverfahrens näherungsweise als
B B (11.3)
berechnen. Dabei wurde angenommen, daß sich das Kontrollvolumen zeitlich nicht ändert,
aber die Methode läßt sich leicht auch auf diesen Fall verallgemeinern.
Die sonstigen Integralterme können nun explizit, implizit oder mit Hilfe des Crank-NicolsonVerfahrens berechnet werden.
Wir werden uns hinfort mit der Zeitableitung nicht mehr beschäftigen, was sich im Weglassen
der Indizes und manifestiert, da schon alles wesentliche hierzu vorher gesagt wurde.
11.2.2 Die Behandlung des Quellterms
Auch beim Quellterm wird die Integration über die finiten Volumina durch ein Produkt mit
denselben ersetzt:
@ B @
(11.4)
An dieser Stelle lassen sich auch numerische Integrationsverfahren höherer Ordnung anwenden; dazu benötigt man allerdings mehr Knoten im Kontrollvolumen.
KAPITEL 11. FINITE-VOLUMEN-VERFAHREN
126
11.2.3 Die Behandlung des Flußterms
Die Behandlung des Flußterms besteht aus einer Integration d.h. Aufsummation der Flüsse
über alle Ränder. Dies sei am Beispiel der in Abb. 11.1 dargestellten Finite-VolumenGeometrie vorgeführt. Die vier Ränder sind in der Abbildung mit n(orth), s(outh), e(ast) und
w(est) bezeichnet. Der Fluß der Größe habe die Komponenten
5 #
#
#
Bei der Integration ist zu beachten, daß positive Flüsse auf der e- und der n-Seite aus dem
Gebiet heraus und w- und s-Seite in das Gebiet hineingehen. Entsprechend ergibt sich:
'
5 A5
#
# ( # # ) # # ( # # ) # Das Integral der Flüsse wird also ganz heuristisch berechnet: Man bilanziert, was in x-Richtung
in das Kontrollvolumen fließt und ihm auf der anderen Seite entfleucht, mal Kantenlänge,
entsprechendes für die y-Richtung.
Es verbleibt nun die Aufgabe, die Flüsse auf den jeweiligen Rändern aus den an den Knoten
bekannten Werten zu interpolieren. Dies werden wir in den folgenden Beispielen konkretisieren.
11.3 Die Diffusionsgleichung auf äquidistantem Gitter
Auf dem in Abb. 11.1 dargestellten zweidimensionalen FV-Gitter soll nun als erstes Beispiel
die Diffusionsgleichung
diskretisiert werden. Integrieren wir die Gleichung über ein zweidimensionales finites Volumen , so erhält die Diffusionsgleichung die Form
5 A5 #
Das Volumen eines Elementes ist dann sein Flächeninhalt
Zeitterm
Der diffusive Fluß ist durch
5
# und wir erhalten für den
11.4. DIE DISKRETISIERUNG ADVEKTIVER TERME
127
gegeben. Auf den Rändern ist es naheliegend, die Ableitungen durch Differenzen der Nachbarknoten zu ersetzen, also für die e-Seite:
Somit ergibt sich für den Flußterm
&
5 A5 #
und insgesamt als Verfahren für die Diffusionsgleichung
wobei man nach Division durch erhält schließlich das Schema
(11.5)
erhält, welches uns schon hinlänglich bekannt ist. Es gelten daher dieselben Aussagen über
Konsistenz, wie sie bei der Darstellung der FD-Methode hergeleitet wurden. Es sei hier nochmals darauf hingewiesen, daß obigen Schema Zeitindizes fehlen. Alle LeserInnen werden in
der Lage sein, an dieser Stelle die verschiedensten Zeitschrittverfahren anzusetzen.
11.4 Die Diskretisierung advektiver Terme
Zur Lösung der Transportgleichung
bleibt die Berechnung des advektiven Flusses
5
#
%* auf den Rändern der Finiten Volumen. Je nach Strategie werden wir auch hier auf alte Bekannte
treffen. Damit es nicht zu langweilig wird, werden wir auch ein neues Verfahren vorstellen,
welches die Möglichkeiten der Finiten-Volumen-Methode illustriert.
11.4.1 Zentrale Verfahren
Bei zentralen Verfahren wird der Fluß auf einem Rand durch Mittlung bzw. Interpolation aus
den benachbarten Knoten. Somit erhält man z.B. für den Fluß in x-Richtung auf dem e-Rand:
%* #
KAPITEL 11. FINITE-VOLUMEN-VERFAHREN
128
und entsprechendes für alle anderen Seiten. Für den advektiven Flußterm ergibt sich insgesamt:
&
5 A5 #
%*
(11.6)
Fügen wir diesen Term zur diskretisierten Diffusionsgleichung, teilen durch so erhalten
wir ein Verfahren, welches den zentralen Differenzen sehr ähnlich ist:
11.4.2 Upstream-Verfahren
Da die Anwendung zentraler Differenzen auf die advektiven Terme zu sehr ungünstigen Stabilitätseigenschaften führt, wollen wir nun im Rahmen der Methode der Finiten Volumen eine
Upstream-Strategie einführen: Dazu ersetzen wir den Fluß auf den Rändern durch den stromaufwärts gelegenen Knoten; für die e(ast)-Seite erhalten wir so für $ %* #
und falls , , dann ist wohl jedem/jeder klar, daß:
%* #
Nach bewährtem Rechengang ergibt sich nun das Verfahren
für $ .
11.4.3 QUICK- und QUICKEST-Verfahren
Das einfache Upstream-Verfahren ist wegen seiner geringen Konsistenzordnung mit einer hohen numerischen Diffusion behaftet. Daher versuchte man Verfahren zu konstruieren, die mit
einer geringeren numerischen Diffusion behaftet sind. Leonard schlug daher 1979 vor, den
Fluß auf dem Rand eines Kontrollvolumens mit Hilfe eines quadratischen Polynoms zu interpolieren. Das Verfahren heißt daher ’Quadratic Upstream Interpolation for Convective Kinetics’ (QUICK), woraus allerdings nicht gefolgert werden kann, daß es besonders schnell
ist.
Betrachten wir die Situation in Abb. 11.2. Bei positiver Geschwindigkeit soll für die Berechnung des Flusses auf dem Rand des Kontrollvolumens die Knoten 0 , 0 und 0
verwendet werden, d.h. es wird mehr Information aus der Upstream-Richtung verwendet. An
diesen Knoten sind die Flüsse # , # und # bekannt, wir suchen allerdings den Fluß auf
dem virtuellen Randknoten 0 , der zwischen 0 und 0 liegt.
11.4. DIE DISKRETISIERUNG ADVEKTIVER TERME
u
0 0
129
0
0
x
Abbildung 11.2: Situation beim QUICK-Verfahren
Dazu wird im ersten Schritt ein Interpolationspolynom durch die Knoten
gelegt. Nach der Lagrangeschen Interpolationsformel ist dieses durch
#
0 , 0 und 0
# #
#
gegeben. Der Leser, der dieses äußerst wichtige Verfahren noch nicht kennt, schaue sich die
Konstruktion genau an und bestätige die Interpolationseigenschaften
#
#
#
# #
#
(11.7)
Für den Spezialfall eines äquidistanten Gitters und dem Rand der Kontrollvolumina auf
der Mitte zwischen zwei Knoten ergibt sich für den Fluß auf dem Rand:
+,-
# # #
#
# # #
# #
(11.8)
Das QUICK-Verfahren wird in verschiedenen Variationen im Bereich der FV-Methode häufig
verwendet. Sein wesentlicher Nachteil ist die Instabilität für die reine Advektionsgleichung,
wenn diese explizit in der Zeit gelöst wird. Leonard konstruierte deshalb ein zweites Verfahren,
welches er auf den Namen QUICKEST (Quadratic Upstream Interpolation for Convective
Kinetics with Estimated Stream Terms) taufte. Der Fluss # wird dabei auf unserer Kante nach
+,-./
# +,-
# */
#
# */
#
# #
(11.9)
mit der Courantzahl
*/
Das Verfahren ist im Gegensatz zum QUICK-Verfahren auch für die reine Advektionsgleichung stabil, jedoch muß auch hier das Courantzahlkriterium eingehalten werden. In der Praxis
der Finiten Volumen wird daher letzteres dem ersteren vorgezogen, weil man bei der Modellierung der turbulenten Viskositäten und Diffusivitäten oftmals nicht ausschließen kann, daß
diese auch mal Null werden.
KAPITEL 11. FINITE-VOLUMEN-VERFAHREN
130
11.5 Globale und lokale Massenerhaltung
Global ist die Massenerhaltung eines Schemas garantiert, wenn der Fluß, der eine Randfläche
eines bestimmten Kontrollvolumens verläßt, auch tatsächlich in das benachbarte Kontrollvolumen einströmt.
Finite-Volumen-Verfahren sind immer dann global massenerhaltend, wenn die Flüsse auf den
jeweiligen Rändern als eigenständige Größen berechnet werden.
Die hier vorgestellten Verfahren (zentrales, upstream und QUICK) sind alle global massenerhaltend.
Im Gegensatz zur globalen Massenerhaltung spricht man von lokaler Massenerhaltung, wenn
auf jedem Knoten die Kontinuitätsgleichung erfüllt ist.
11.6 Gestaffelte Gitter
Die numerische Behandlung von Strömungen mit Hilfe finiter Volumina hat uns gezeigt, daß
Flußterme immer auf dem Rand der Kontrollvolumina berechnet werden sollten. Da die physikalischen Größen bisher nur auf Knoten im Zentrum der Volumina definiert wurden, war man
zu Interpolationen gezwungen.
Eine Möglichkeit zur Vermeidung dieser Interpolationen ist die Verwendung gestaffelter Gitter, bei denen
für jede physikalische Größe eigene Knoten definiert
und für jede Grundgleichung eigene Kontrollvolumina definiert werden.
Zur Lösung etwa der Kontinuitätsgleichung der zweidimensionalen Gleichungen von SaintVenant
betrachten wir die Variablenanordnung in Abb. 11.3. Die Variable wurde im Zentrum des
finiten Volumens, die Variable an der e- und w-Seite und die Variable an der n- und s-Seite
angeordnet, da sie jeweils nur dort zu den Flüssen beitragen.
Integriert man über das angegebene Kontrollvolumen, so erhält man das Schema
Dabei wurden die Knoten auf den Rändern mit den -Indizes bezeichnet. Der Wasserstand
wurde auf den Rändern jeweils zentral interpoliert.
Aufgabe: Zeichnen Sie bei derselben Anordnung der Knoten ein Kontrollvolumen für die
Lösung der x-Impulsgleichung und konstruieren Sie das FV-Schema.
11.6. GESTAFFELTE GITTER
131
1
x
x
x
x
1
x
x
x
x
1
x
x
x
x
0
0
0
Abbildung 11.3: Zweidimensionales gestaffeltes Gitter. Ausgefüllte Kreise: -Knoten, Kreuze: -Knoten, leere Kreise: -Knoten. Die dargestellte Aufteilung der Kontrollvolumina eignet
sich zur Lösung der Kontinuitätsgleichung.
KAPITEL 11. FINITE-VOLUMEN-VERFAHREN
132
#
#
# Abbildung 11.4: Geometrie für Finite Volumen auf Dreiecken
11.7 Finite Volumen auf Dreiecksnetzen
Die Methode der Finiten Volumen läßt sich auch auf Dreiecksnetze verallgemeinern. Dazu
betrachten wir das in Abb. 11.4 dargestellte Dreieck als finites Volumen , auf dem wir die
allgemeine Bilanzgleichung
5
div # @
lösen wollen.
Eine glückliche Fügung der Geometrie läßt alle Mittelsenkrechten in einem Dreieck sich in
einem Punkt schneiden. Wir wählen diesen Mittelpunkt des das Dreiecke umschreibenden
Kreises als Stützstelle für die Größen und @ , die im Dreieck 0 die Werte und @ haben. Auf
den Mittelpunkten der Kanten sei die Projektion des Flusses auf den Aussennormalenvektor
5 5 # bekannt. Wir bezeichnen diese Normalenflüsse auf den drei Kanten des Dreieckes 0
#
mit # , # und # . Jeder dieser Flüsse soll den Mittelwert auf der entsprechenden Kante
repräsentieren.
Nun wird die allgemeine Bilanzgleichung über die Dreiecksfläche integriert und die Divergenz mit Hilfe des Gausschen Integralsatzes auf ein Randintegral reduziert
# ? # ? # ? @
5 nicht von seiner Bezugswobei ? die Länge der k-ten Kante im Dreieck 0 ist. Ist der Fluss #
größe abhängig, dann kann eine explizite Diskretisierung für die Zeitableitung gewählt wer
den. Dies ist aber eigentlich nie der Fall, so daß eine implizite Diskretisierung oder ein CrankNicolson-Verfahren vorgezogen werden sollte. Der Übersichtlichkeit halber beschränken wir
uns auf die Darstellung für die implizite Diskretisierung der Zeitableitung:
# ? @
11.7. FINITE VOLUMEN AUF DREIECKSNETZEN
133
Ist der Fluss # linear von der Größe abhängig, also # , dann ist die implizite Diskre
tisierung einfach #
. Oftmals ist der Fluss # allerdings nichtlinear von der Größe
abhängig, also # . Dann kann man den Term durch die Zeitableitung linearisieren,
also
# setzen. In beiden Fällen benötigt man die Größe auch auf den Kanten des Dreieckes, was
durch den Index 0 schon angedeutet wurde, womit die k-te Kante des i-ten Dreieckes bezeichnet wird.
Der Fluss #
auf der k-ten Kante des i-ten Dreieckes muss in irgendeiner Form aus den
Werten der Größe aus den benachbarten beiden Dreiecken interpoliert werden. Wir wollen
hier nicht spezifizieren, um welche beiden Dreiecke es sich da handelt und schreiben diese
Interpolation in der allgemeinen Form:
# #
1 So werden alle zentralen Knoten der ; Dreiecke im Gitter berücksichtigt, wobei bei einer
linearen Interpolation nur zwei Gewichte 1 ungleich Null sind. An dieser Stelle können
aber auch höhere Interpolationsverfahren berücksichtigt werden. Wir setzen nun den Fluß in
die diskretisierte Bilanzgleichung ein und bekommen ihre Enddarstellung:
#
1 ? @
Um sie zu lösen, wollen wir sie in die Normalform eines linearen Gleichungssystems bringen. Die Komponenten der Matrix sind ganz offensichtlich
Æ 1 ? wobei Æ das Kroneckersymbol ist und die Komponenten der rechten Seite haben die Form
Schließlich kann das so entstandene Gleichungsystem mit einem Gleichungslöser gelöst werden.
Es sei ausblickend bemerkt, daß dieses Verfahren auf Elementformen mit beliebig vielen
Ecken ausdehnbar ist. Die einzige Bedingung, die jedes Element erfüllen muß, daß die Mittelsenkrechten eines Elementes sich immer in einem Punkt schneiden. Ein aus solchen Basiselementen aufgebautes Gitter nennt man orthogonal strukturiertes Gitter, ein aus Dreiecken und
Vierecken aufgebautes Teilgitter ist in Abbildung 11.5 dargestellt.
134
KAPITEL 11. FINITE-VOLUMEN-VERFAHREN
Abbildung 11.5: Ein Teil eines orthogonal strukturierten Gitters.
11.8. FINITE DIFFERENZEN UND FINITE VOLUMEN
135
11.7.1 Der äußere Normalenvektor im Dreieck
Zur Bestimmung des Flusses auf den Kanten des Dreieckes benötigen wir die dortigen Außennormalenvektoren. Ist das Dreieck durch die drei Punkte , und bestimmt, dann kann man schnell den einfachen Algorithmus
herleiten, der die Komponenten eines Normalenvektors zur Kante
Dieser Vektor weist aus dem Dreieck hinaus, falls
berechnet.
, ansonsten ist sein Vorzeichen in beiden Komponenten umzukehren.
11.7.2 Upstreaming auf Dreiecksnetzen
Auch auf Dreiecksgittern kann die explizite Diskretisierung der advektiven Terme zu Instabilitäten führen, wenn man kein Upstreaming in das Verfahren einbaut. Sei also der Fluss #, den
wir auf den Kanten bestimmen wollen, von der Größe abhängig, die nur im Mittelpunkt des
Dreieckes bekannt ist, # # . Entsprechend dem Upstream-Verfahren bei strukturierten
Gittern berechnen wir den Fluss auf einer Kante 1 im Dreieck 0 dann mit dem Wert , wenn
der Flußvektor aus dem Dreieck weist. Im umgekehrten Fall wird der Wert aus dem anderen,
an derKante gelegenem Dreieck ? verwendet.
Wie das Analogon auf strukturierten Gittern besitzt dieses Verfahren eine gewisse numerische
Diffusion. Upstream-Verfahren höherer Ordnung, wie z.B. das QUICK-Verfahren, sind auf
Dreiecksnetzen bisher nicht verwirklicht worden.
11.8 Finite Differenzen und Finite Volumen
Finite Differenzen und Finite Volumen haben topologisch gesehen vollkommen unterschiedliche Ansätze: Während erstere die Gleichungen auf einzelne Punkte reduziert, geht zweitere
von einer Bilanzierung auf endlichen Volumen aus. Die Methode der Finiten Volumen ist daher
dem Strömungsgeschehen wesentlich näher als die Methode der Finiten Differenzen.
Die Methode der Finiten Volumen ist vom geometrischen Standpunkt eine Verallgemeinerung der FD-Methode, da sie prinzipiell auch Verfahren für beliebige Gittergeometrien zur
Verfügung stellt. Damit birgt sie das Potential in sich, die erheblich rechenaufwendigere und
instabilere Finite-Elemente-Methode in der numerischen Hydrodynamik abzulösen. Hier besteht allerdings noch Entwicklungsbedarf.
Die Unterscheidung von Finiten Differenzen und Finiten Volumen ist in der numerischen
Strömungsmechanik heute fließend. In vielen Programmsystemen werden daher FD/FVKombinationen verwendet: Auf einfachen z.T. äquidistanten Geometrien wird gestaffelt gerechnet, wobei die Verfahren flußorientiert entwickelt wurden.
136
KAPITEL 11. FINITE-VOLUMEN-VERFAHREN
Kapitel 12
Die Saint-Venantschen Gleichungen
Werden Strömungen durch die Burgersgleichung modelliert, so werden alternierend Ansammlungen und Senken der Fluidmasse produziert. Diese Verhalten wird dadurch möglich, weil
es in der Burgersgleichung keinen Mechanismus gibt, der auf diese Kapriolen der Fluidmasse
reagieren kann und das Fluid in eine gleichmäßigere Verteilung zurückdrängt. Stellt man sich
hierzu eine Strömung in einem Kanal vor, so könnte die rückstellende Kraft in der Form
modelliert werden, wobei die Wassertiefe ist. Steigt also der Wasserstand an einem Ort in
x-Richtung steil an, so wird das Wasser in entgegengesetzter Richtung wieder zurückgedrängt.
In der Hydromechanik der Oberflächengewässer wird gezeigt, daß die fehlende Proportionalitätskonstante die Gravitationsbeschleunigung > sein muß und somit die Gleichung
>
eine naturnahe Strömung in einem kanalartigen Gerinne reproduzieren sollte.
Bleibt nur die Frage nach dem Wasserspiegelgefälle also nach dem zeitlichen Verhalten der
Wasserspiegellinie offen. Für sie benötigen wir eine zweite Bestimmungsgleichung, um das Problem zu schließen. Dabei lassen wir uns von dem Bild leiten, daß Verformungen der Wasseroberfläche konservativ advektiert werden, also durch die Gleichung
beschrieben werden. Die Saint-Venant-Gleichungen stellen eine solche Kanalströmung unter
Vernachlässigung des viskosen Terms dar, sie lauten also:
>
Mathematisch lassen sie sich als nichtlineares partielles Differentialgleichungssystem erster
Ordnung charakterisieren. Jede Gleichung enthält neben der Advektion einen Kopplungsterm
137
KAPITEL 12. DIE SAINT-VENANTSCHEN GLEICHUNGEN
138
mit der jeweiligen anderen Gleichung, im System können wir insgesamt drei nichtlineare Terme identifizieren.
Im zweidimensionalen Fall sind die Saint-Venant-Gleichungen durch
>
>
gegeben. Sie stellen nur einen Teil der tiefen- und breitenintegrierten Hydrodynamik der Fließgewässer dar, da sie in sehr eng umgrenzten Spezialfällen gültig sind. Um die folgenden Entwicklungen allgemeiner zu halten, führen wir unbestimmte Quellterme , und * ein, die
alle vernachlässigten Prozesse enthalten sollen. Somit ist die allgemeinere Form der eindimensionalen Saint-Venant-Gleichungen
>
(12.1)
und die ihr zweidimensionales Analogon lautet:
>
>
*
(12.2)
Eine hinreichend genaue physikalische Beschreibung der eindimensionalen Fließgewässerströmungen ergibt sich für die Belegung:
$
%
> !
"
* *
# " "
Dabei sind der Gewässerquerschnitt, dessen Breite, ! die Lage der Sohle, * der Chezywert der Sohlschubspannung, " die Luftdichte, " die Windgeschwindigkeit und * # der
Windschubkoeffizient. Diese Aufzählung läßt erahnen, wieviele physikalische Effekte in dieser Erweiterung berücksichtigt werden können.
12.1. HYPERBOLIZITÄT DER SAINT-VENANT-GLEICHUNGEN
139
12.1 Hyperbolizität der Saint-Venant-Gleichungen
Nach der Definitionsgleichung (5.2) lassen sich Differentialgleichungssysteme erster Ordnung
in hyperbolische und elliptische einteilen. Um zu entscheiden, welcher Fall bei den SaintVenant-Gleichungen vorliegt, identifizieren wir in der zitierten Gleichung mit , mit und mit .
Die Matrix ist im eindimensionalen Fall dann
> Ihre Eigenwerte sind > und damit grundsätzlich reell, womit die Saint-VenantGleichungen hyperbolisch sind. Matrizen mit reellen Eigenwerten lassen sich grundsätzlich
diagonalisieren, d.h. es gibt ein bestimmtes Koordinatensystem, in welchem die Matrix diagonal ist. Ist die Matrix diagonal, dann taucht in jeder der beiden Differentialgleichungen nur
noch eine Ableitung auf, die partiellen Differentialgleichungen entarten in diesem speziellen
Koordinatensystem also zu gewöhnlichen. Diese Eigenschaft macht sich die Charakteristikenverfahren zunutze, die wir in den nächsten Abschnitten entwickeln wollen.
Im zweidimensionalen Fall hat man die Matrizen und zu untersuchen, sie sind:
>
und
Ihre Eigenwerte sind und > bzw. und chungssystem auch hier hyperbolisch.
>
>, also reell und somit ist das Glei-
12.2 Charakteristiken der 1D-Gleichungen
Wir wollen nun untersuchen, wie wir die eindimensionalen Saint-Venant-Gleichungen in
gewöhnliche Differentialgleichungen umformen können. Dies sollte aufgrund ihrer Hyperbolizität möglich sein.
Ersetzt man in beiden Gleichungen die Wassertiefe durch die sogenannte Wellengeschwindigkeit
> (12.3)
multipliziert die erste Gleichung mit > und teilt die Kontinuitätsgleichung durch , so ergibt
sich
>
Beide Gleichungen werden nun jeweils addiert und subtrahiert:
0 0 140
KAPITEL 12. DIE SAINT-VENANTSCHEN GLEICHUNGEN
0
0
0
Abbildung 12.1: Charakteristiken der 1D-Gleichungen von Saint-Venant
Auf der linken Seite beider Gleichung steht der Form nach eine totale Ableitung der Variablen
bzw. nach der Zeit auf Kurven, die durch
>
(12.4)
>
(12.5)
gegeben sind. Wir nennen sie charakteristische Kurven oder einfach Charakteristiken * und *
. Auf den Charakteristiken reduzieren sich die SaintVenantschen Gleichungen mit Hilfe der totalen Differentiale zu:
>
>
auf
*
auf
*
Mit der ersten gewöhnlichen Differentialgleichung können wir die zusammengesetzte Funktion , mit der zweiten bestimmen. Die geschickte Kombination dieser Lösungen
liefert uns dann die gesuchten Funktionen und bzw. .
12.2.1 Das inverse Charakteristikenverfahren
Zur Lösung der Saint-Venantschen Gleichungen betrachte man den Verlauf der Charakteristiken, die von einem Knoten zum unbekannten Zeitschritt in Richtung des bekannten
Zeitschrittes ausgehen in Abb. 12.1.
Die Basispunkte der Charakteristiken und zur Zeit ergeben sich durch die Lösung
ihrer Differentialgleichungen 12.4 und 12.5. Dazu integrieren wir diese über die Zeit:
> > 12.2. CHARAKTERISTIKEN DER 1D-GLEICHUNGEN
141
Auf die linke Seite wird jeweils der Hauptsatz der Differential- und Integralrechnung, auf die
rechte Seite die Trapezregel angewendet:
!
> !
> !
> (12.6)
> (12.7)
!
In derselben Weise werden die Saint-Venantschen Gleichungen auf den Charakteristiken über
die Zeit integriert:
> > auf
*
auf
*
Zur Abkürzung setzen wir
) * ) * > >
>
> >
> (12.8)
(12.9)
und man erhält:
) * ) * auf
auf
*
*
Addition und Subtraktion der beiden Gleichungen und ein wenig Umformen liefern schließlich
die Lösungen:
>
)
* ) * ) * ) * (12.10)
(12.11)
Das Verfahren ist nun eigentlich vollständig, nur daß schon in der Berechnung der Basispunkte
die gesuchte Lösung als auch die Werte von und an der Basis implizit erscheinen. Wir
gehen daher in jedem Zeitschritt iterativ vor:
1. Schritt: Setze für alle Knoten , .
, , und 2. Schritt: Berechne zu jedem Knoten die Basisknoten und .
,
KAPITEL 12. DIE SAINT-VENANTSCHEN GLEICHUNGEN
142
*
Abbildung 12.2: Zum Einbau von Randbedingungen
3. Schritt: Interpoliere die Basiswerte , , und .
4. Schritt: Berechne die Integrale )
5. Schritt: Berechne und für jeden Knoten.
* und ) * für alle Knoten.
6. Schritt: Zur Fehlerabschätzung berechne nochmals die Basisknoten und vergleiche
mit den alten: Ist die maximale Differenz der beiden zu groß, gehe zurück zu Schritt 2.
Dieses Verfahren ist für jeden Zeitschritt zu wiederholen. Die Konsistenzordnung des Charakteristikenverfahrens ergibt sich dann in der Zeit aus dem Exaktheitsgrad der Zeitintegration
und im Raum aus der Ordnung des verwendeten Interpolationsschemas. Das Verfahren ist unabhängig vom Zeitschritt stabil.
12.2.2 Randbedingungen
An den Randknoten stehen zur Anwendung des Charakteristikenverfahrens jeweils nur eine,
d.h. die in das Gebiet laufende Charakteristik zur Verfügung. Wir betrachten die Situation am
linken Rand am Knoten in Abb. 12.2. Es steht hier als einlaufende Charakteristik * zur
Verfügung. Der Basispunkt berechnet sich als
Auf der Charakteristik gilt:
!
> !
> > !
>
) * (12.12)
(12.13)
Die Gleichung ist explizit lösbar, wenn wir entweder oder kennen. Damit muß entweder die Wassertiefe oder die Geschwindigkeit als Randbedingung vorgegeben werden.
Für einen -Rand ( vorgegeben) folgt dann:
!
> !
>
) * (12.14)
Der rechte Rand erfährt im Charakteristikenverfahren die entsprechende Sonderbehandlung.
12.2. CHARAKTERISTIKEN DER 1D-GLEICHUNGEN
143
*
*
* Abbildung 12.3: Eine linkslaufende Charakteristik verläßt das Strömungsgebiet.
Ein weiterer Sonderfall tritt ein, wenn eine von einem Innenknoten ausgehende Charakteristik
das Gebiet verläßt. Zur Analyse der Situation betrachte man Abb. 12.3. Die Idee besteht darin, die linkslaufende Charakteristik * nur bis zum Rand zu verfolgen. Da von jedem Punkt
in Raum und Zeit zwei Charakteristiken und somit auch vom Randschnittpunkt eine rechtslaufende Charakteristik * ausgeht, verfolgen wir diese weiter bis zur Basis am Zeitschritt
.
Im Programm stellen wir also fest, daß unser Basiswert außerhalb des Strömungsgebietes
ist, d.h. also kleiner als ist. Wir berechnen dann sofort den reduzierten Zeitschritt als
(12.15)
Für einen -Rand (die entsprechende Prozedur für einen -Rand sei dem Leser zur Übung
überlassen) interpolieren wir die Randbedingung an diesem Durchstoßzeitpunkt:
(12.16)
Der Basispunkt der Charakteristik ergibt sich nun aus der Integration der Differentialgleichung
>
zwischen und . Man erhält
! > > Schreiben wir nun die Bedingungen auf den drei Charakteristiken hin:
) * ) * ) * auf
auf
auf
*
*
*
(12.17)
KAPITEL 12. DIE SAINT-VENANTSCHEN GLEICHUNGEN
144
Für die drei Unbekannten , und stehen drei Gleichungen zur Verfügung. Durch
die Auflösung erhält man:
) * )
(12.18)
* ) * (12.19)
) * ) * (12.20)
>
Die expliziten Lösungen für und unterscheiden sich formal nicht vom Normalfall,
bei dem die Basispunkte im Gebiet liegen. Es ändert sich lediglich die Basispunktberechnung
für die auslaufende Charakteristik.
12.3 Charakteristiken der 2D-Gleichungen
Die zweidimensionalen Saint-Venant-Gleichungen lassen sich in Matrixform als
>
>
darstellen.
Das Ziel ist es wieder, die Charakteristiken, d.h. die Kurven in Raum und Zeit zu suchen,
auf denen die drei partiellen Differentialgleichungen zu gewöhnlichen entarten: Diese Kurven
werden Charakteristiken genannt.
Dies kann dadurch geschehen, daß wir die partiellen Ableitungen in der Zeit durch totale
Ableitungen auf den Charakteristiken ersetzen. Dazu versuchen wir, die Matrizen zu diagonalisieren, wobei wir die Eigenwerte der ersten Matrix formal mit %
% und die der zweiten Matrix
%
mit % bezeichnen. Die diagonalisierten Gleichungen nehmen dann die Gestalt
>
Die Eigenwerte der Matrizen bestimmen sich aus
12.3. CHARAKTERISTIKEN DER 2D-GLEICHUNGEN
det >
det >
145
>
>
Folglich lassen sich die tiefengemittelten Gleichungen auf der Charakteristik
und
(12.21)
>
und
>
(12.22)
>
und
>
(12.23)
die der Bahnlinie entspricht, und den Charakteristiken
sowie
die man manchmal auch als Bicharakteristiken bezeichnet, auf gewöhnliche Differentialgleichungen mit totalen Ableitungen reduzieren.
Da die gesamte Herleitung unabh ängig vom gewählten Koordinatensystem sein muß, gelten
die Aussagen auf dem gesamten sogenannten bicharakteristischen Kegel der Form
>
(12.24)
Für einen gegebenen Raum-Zeit-Punkt formen die diesen Punkt schneidenden Kurven
den sogennannten bicharakteristischen Kegel.
Die physikalische Interpretation des bicharakteristischen Kegels ist direkt mit der Hyperbolizität des Problems verbunden: Der von einem beliebigen Raum-Zeit-Punkt in die
Zukunft laufende Kegel umrandet das Gebiet, welches durch die St örung von diesem RaumZeit-Punkt beeinflußt wird (Abb. 12.4) und wird somit Einflußgebiet genannt.
Umgekehrt wird der Punkt aus der Vergangenheit nur von den Gebieten beeinflußt, die innerhalb des r ückwärtig verlängerten bicharakteristischen Kegels liegen. Das so
definierte Gebiet wird daher Abhängigkeitsgebiet genannt.
12.3.1 Die stationären tiefengemittelten Gleichungen
Im stationären Fall, bei dem die Zeitableitungen aus den Gleichungen verschwinden, ändern
sich die Verhältnisse drastisch. Die tiefengemittelten Gleichungen nehmen in der Matrixschreibweise die Form
KAPITEL 12. DIE SAINT-VENANTSCHEN GLEICHUNGEN
146
D o m a in o f in flu e n c e
o f ( x i, y j ) a t t n + 1
n + 1
j
y
n + 1
t
y
t
t
n
j
n
n - 1
y
n - 1
j
x
i
n - 1
D o m a in o f d e p e n d e n c e
o f ( x i, y j ) a t t n - 1
X
Abbildung 12.4: Der bicharakteristische Kegel im Fall einer konstanten Geschwindigkeit 5 .
>
>
Da die Zeit nicht mehr als Variable zur Verfügung steht, wählen wir ohne Beschränkung der
Allgemeinheit die y-Richtung als Parameter für die zu bestimmenden Charakteristiken. Auf
ihnen sind die totalen Ableitungen einer Funktion dann
Löst man nach der lokalen Ableitung auf und ersetzt diese in den Gleichungen
>
>
Der Term mit den verbleibenden lokalen Ableitungen entfällt genau dann, wenn seine Determinante Null ist:
12.3. CHARAKTERISTIKEN DER 2D-GLEICHUNGEN
det 147
>
was auf den Charakteristiken
und
der Fall ist [26]. Im Gegensatz zum hyperbolischen instationären Fall sind die Gleichungen
bei schießenden Strömungen (
und bei strömenden Strömungen (
$
) hyperbolisch
, ) elliptisch.
Dieser Sachverhalt ist bei der numerischen Behandlung der stationären Gleichungen zu
berücksichtigen.
12.3.2 Inverse Bicharakteristikenverfahren für die 2D-Gleichungen
Es sei nun ein numerisches Verfahren vorgestellt, welches die Navier-Stokes-Gleichungen mit
Hilfe der Bicharakteristiken löst. Dazu ist zu untersuchen, wie sich die Lösung auf den Stromlinien als auch den Bicharakteristiken verhält. Die dieses Verhalten darstellenden Gleichungen
heißen Kompabilitätsbedingungen [?].
Indem die Terme der Kontinuitätsgleichung umgestellt werden zu
(12.25)
liest sich die Kompabilitätsbedingung auf einer Stromlinie ab als
(12.26)
Der bicharakteristische Kegel wird durch die Einführung des Winkels ( parametrisiert:
> (
(12.27)
> (
(12.28)
Somit erhält man für die totale Ableitung der Wassertiefe
Kegel mit der Kontinuitätsgleichung
auf einem bicharakteristischen
KAPITEL 12. DIE SAINT-VENANTSCHEN GLEICHUNGEN
148
> ( > (
> ( > ( (12.29)
Und in der selben Art erhält man für die Impulsgleichungen
> ( > (
>
> ( > ( (12.30)
und
> ( > (
> ( > ( (12.31)
>
mit > , (12.30) mit > ( und (12.31) mit > ( multipliziert und
Werden (12.29)
schließlich addiert, so erhält man für die Kompabilitätsbedingung auf den bicharakteristischen
Kurven
>
> ( > ( (
>
( ( (
>
(12.32)
( ( Zur numerischen Lösung kann ein inverses tetrahedrales Netz [?] benutzt werden: Dazu betrachte man drei bicharakteristische Kurven (Abb. 12.5), die durch die willkürlichen Winkel
( 7 7 gegeben seien und sich von einem Knoten ( ) in der Zeit
zurück bewegen (Abb. 17). Wieder werden die Basispunkte der drei Kurven zur Zeit mit
den Gleichungen (12.27) und (12.28) bestimmt und die Werte von , und an diesen Punkten interpoliert. Aus den Kompabilitätsbedingungen (12.32) für alle drei Bicharakteristiken
können nun die gesuchten Werte zum Zeitpunkt durch das Ersetzen der totalen Ableitungen durch Vorwärtsdifferenzen bestimmt werden.
In einem inversen pentahedralen Netz [19] werden vier bicharakteristische Kurven und die
Stromlinie verwendet, so daß auch die partiellen Ableitungen in (12.26) ohne Interpolation
bestimmt werden können.
12.4 Das semiimplizite Verfahren von Casulli
Das semiimplizite Verfahren von Casulli [5] ist ein Lösungsverfahren für Differentialgleichungssysteme, die die Kontinuitätsgleichung enthalten. Diese enthält als Lösungsfunktionen
12.4. DAS SEMIIMPLIZITE VERFAHREN VON CASULLI
(x
(x
(x
(x
n
i- 1
,y
n
j+ 1
)
(x
(x
i
n
n
1
,y
b
n
i- 1
,y
,y
1
n
j
(x
(x
)
j+ 1
)
b
b
n
i+ 1
i
n
)
,y
)
b
,y
3
n
3
j
)
,y
n
j+ 1
i
n +1
)
(x
(x
(x
(x
n
i- 1
,y
n
j- 1
)
n +1
j
,y
149
b
n
i+ 1
n
,y
b
,y
2
i
n
,y
n
)
j
)
2
j-1
(x
n
i+ 1
,y
n
i- 1
)
)
)
Abbildung 12.5: Das inverse tetrahedrale Netzwerk auf einem rechteckigen Gitter.
die Wassertiefe und die Geschwindigkeit , wobei durch letztere eine Kopplung mit der Impulsgleichung entsteht. Ziel ist es, die Kontinuitätsgleichung implizit in und explizit in zu
diskretisieren, wodurch die Gleichungen entkoppelt werden. Aus Stabilitätsgründen sind jedoch in der Kontinuitätsgleichung auch implizite Anteile von zu berücksichtigen. Sie haben
daher die allgemeine zeitdiskrete Form:
Die Operatoren und können dabei durchaus von und abhängen. Problem ist ,
da es an dieser Stelle noch nicht bekannt ist. Es wird aus einem zeitdiskreten Verfahren für die
Kontinuitätsgleichung gewonnen, welches die folgende Form haben muß:
Setzt man dieses Verfahren in die Kontinuitätsgleichung ein
dann bleibt eine implizite Gleichung zurück, in der nur noch zum neuen Zeitschritt gesucht
wird. Nach der Lösung dieser Gleichung läßt sich nun auch mit Hilfe der vorangegangenen
Gleichung bestimmen.
Das Verfahren funktioniert also nur, wenn man die Geschwindigkeitsterme in der Impulsgleichung explizit und stabil diskretisieren kann.
In den von Casulli entwickelten Verfahren wird der semiimplizite Algorithmus auf einem
gestaffelten Finite-Volumen-Gitter angewendet. In den Zentren der Volumina ist dabei jeweils
KAPITEL 12. DIE SAINT-VENANTSCHEN GLEICHUNGEN
150
die Wassertiefe und auf den Rändern die Strömungsgeschwindigkeit abgespeichert, letztere werden durch halbzahlige Indizes gekennzeichnet.
Die Integration der Kontinuitätsgleichung in der konservativen Form
über das Volumen liefert
Dabei ist die Stabilität des Verfahrens garantiert, wenn lediglich im Flußterm implizit
gewählt wird. Auf den Volumenrändern berechnet sich der Wasserstand als Mittelwert der
angrenzenden Volumina:
Die so diskretisierte Kontinuitätsgleichung ist noch nicht lösbar, da uns die Geschwindigkeit
zum neuen Zeitschritt fehlt. Jene gewinnen wir aus der Impulsgleichung, auf die ein
Lagrangeverfahren angewendet wird. Dies ist absolut stabil, in der zeitdiskretisierten Form
schreibt es sich als:
>
Die Geschwindigkeit ist an Knoten mit Halbindizes abgespeichert und die Ableitung der
Wassertiefe wird durch eine zentrale Differenz diskretisiert:
>
(12.33)
Die in der Kontinuitätsgleichung auftauchenden Geschwindigkeiten und werden
schließlich substituiert:
> (12.34)
Der Algorithmus des Gesamtverfahrens besteht aus drei Schritten:
1. Schritt: Berechne die Werte der Geschwindigkeiten an den Basen der Bahnlinien
.
2. Schritt: Löse das Gleichungssystem (12.34) für die Wassertiefe .
3. Schritt: Löse die Impulsgleichungen (12.33) für die Geschwindigkeit .
Den ersten Schritt hatten wir schon eingehend in Kapitel 8 untersucht. Wir wollen nun
die Eigenschaften des Gleichungssystems (12.34) untersuchen. Laufen wir die vorhandenen
Lösungsknoten von 0 bis
durch, so sehen wir, daß die Gleichung für den ersten und
12.4. DAS SEMIIMPLIZITE VERFAHREN VON CASULLI
151
letzten Knoten nicht gültig sein kann, hier fehlen uns Randbedingungen, die Werte und
müssen also vorher bekannt sein. Das Gleichungssystem (12.34) schreibt sich nun in
Matrizenform als:
..
..
.
..
.
.
...
...
Im Inneren der Matrix sind Diagonal- und Extradiagonalelemente durch Gleichung (12.34)
vorgegeben, sie lauten:
> > Die Matrixelemente der rechten Seite sind:
Offensichtlich ist das Innere der Matrix symmetrisch, und die Randbedingungen zerstören
diese Symmetrie. Wir wollen diese also aus dem zu lösenden Gleichungssystem herausnehmen
und reduzieren dieses auf die Knoten zwischen 0 und :
..
.
..
.
..
.
...
...
..
.
..
.
Die Matrix ist nun symmetrisch und positiv definit. Sie zählt daher zu dem denkbar gutmütigsten Matrizentyp und kann mit nahezu allen Gleichunglösern gelöst werden.
Das Verfahren ist unabhängig von Zeitschritt stabil. Der Beweis hierzu ist etwas schreibintensiv und soll daher nur angedeutet werden. Er verwendet die in Abschnitt 3.3.2 dargestellte
Beweistechnik. Die linke Seite der Stabilitätsbedingung
besteht aus dem Betragsquadrat des Vektors und ist daher immer positiv. Die rechte Seite
ist immer negativ. Dies liegt am negativen Vorzeichen der , die die Beziehung erfüllen. Die auf der Diagonalen der Matrix erscheinenden mit negativem Vorzeichen
heben sich jeweils mit dem positiven, aber ansonsten identischen Term der nächsten Zeile auf.
152
KAPITEL 12. DIE SAINT-VENANTSCHEN GLEICHUNGEN
12.5 Bewertung
Für die Saint-Venantschen Gleichungen existieren noch eine Vielzahl von weiteren Lösungsverfahren, die Tan Weiyan in seinem 1992 erschienenen Buch Shallow Water Hydrodynamics
[27] zusammengestellt hat. Die Mehrheit dieser Lösungsverfahren krankt an der Stabilität.
Sowohl die Charakteristikenverfahren als auch das semiimplizite Verfahren von Casulli sind
stabil, letzteres ist in dem zitierten Buch noch nicht beschrieben. Die Implementation des Charakteristikenverfahrens für den zweidimensionalen Fall ist allerdings extrem umständlich und
es ist auch nur schwer erweiterbar für die bisher nicht berücksichtigten Prozesse in tiefenintegrierten Strömungen. Da dies beim Casulliverfahren nicht der Fall ist, soll es als Standardverfahren zur Lösung der Saint-Venantschen Gleichungen empfohlen werden.
Kapitel 13
Galerkinverfahren
Während die Methode der Finiten Differenzen die Differentialoperatoren diskretisiert, wird
bei der Methode der Finiten Elemente der Raum der Lösungsfunktionen diskretisiert. Die Differentialoperatoren bleiben dabei unangetastet.
Der Gewinn dieses Vorgehens gegenüber den Finiten Differenzen besteht im besseren geometrischen Auflösungsvermögen der Methode der Finiten Elemente. Mit ihr lassen sich Differentialgleichungen im Zweidimensionalen auch auf Dreiecksgittern oder im Dreidimensionalen
auf Tetraedern lösen. Dadurch können komplexe Ränder besser aufgelöst und lokale Verfeinerungen des Netzes berücksichtigt werden.
Da aber nun auch Finite Volumen auf unstrukturierten Geometrien (d.h. insbesondere Dreiecksgittern im zweidimensionalen) angewendet wurden, wie sie vorher nur den Finiten Elementen vorbehalten waren, ist ein Wettkampf zwischen diesen beiden Methoden eingeläutet
worden, bei dem aller Voraussicht nach die Finiten Elemente in der reinen Hydromechanik
den Kürzeren ziehen werden. Sie sind nicht nur rechenintensiv und recht instabil, sondern
auch nicht konservativ.
Ein wesentlicher Vorteil der FE-Verfahren ist die Produktion von sogenannten schwachen
Lösungen, die auch mit unstetigen Materialparametern fertig werden, so wie sie an den Kontaktflächen zwischen zwei Medien auftreten. Wenn nur die Wassersäule eines Gewässers modelliert wird, hat man nicht mit unstetigen Materialeigenschaften zu kämpfen. Dies ist in der
Kontinuumsmechanik anders, daher wird die FE-Methode dort auch intensiv genutzt.
In diesem Kapitel werden wir das Galerkinverfahren als Kochrezept zur Finite Elemente Methode kennenlernen.
13.1 Das Standard-Galerkinverfahren
Wir wollen die Lösung der recht abstrakten Gleichung
(13.1)
suchen. sei ein linearer Differentialoperator. Beim Standard-Galerkin-Verfahren wird im
ersten Schritt der unendlichdimensionale Raum der Lösungsfunktionen durch einen endlichdimensionalen Teilraum ersetzt, welcher durch die Basisfunktionen
( (
153
KAPITEL 13. GALERKINVERFAHREN
154
aufgespannt sein soll.
Die gesuchte Lösung (genauer: die Projektion der Lösung) läßt sich dann als
(
darstellen und in die Differentialgleichung einsetzen:
(
Um hieraus ein algebraisches Gleichungssystem zu gewinnen, wird die Gleichung der Reihe
nach mit den Ansatzfunktionen ( , 1 multipliziert und über das Lösungsgebiet
integriert:
(( ( 1 (13.2)
In ihrer Funktion als Multiplikatoren bezeichnet man die Ansatzfunktionen an dieser Stelle
auch als Wichtungsfunktionen. Sie sind der Träger des Integrals, denn nur dort wo die Wichtungsfunktionen nicht Null sind, wird das Integral jeweils ausgewertet.
Es ist somit ein lineares Gleichungssystem mit der Matrix ( ( und der rechten Seite
( entstanden. Gesucht sind die Unbekannten
, ... , , aus denen sich dann die
vollständige Lösung konstruieren läßt.
13.2 Die Poissongleichung in der Hydrodynamik
Eine der wichtigsten partiellen Differentialgleichungen ist die Poissongleichung, sie hat die
Form:
Sie taucht in der Hydromechanik an verschiedenen Stellen auf, wir wollen hier nur zwei Beispiele herausnehmen, die Potentialströmungen, sowie die Druck-Korrektur-Verfahren.
13.2.1 Potentialströmungen
Die Massenerhaltung in inkompressiblen Fluiden wird durch die Kontinuitätsgleichung
beschrieben. In der Theorie der Potentialströmungen geht man davon aus, daß sich das Vektorfeld der Geschwindigkeiten rotationsfrei ist, d.h.:
rot 5 13.2. DIE POISSONGLEICHUNG IN DER HYDRODYNAMIK
155
Dann läßt sich zeigen, daß die Geschwindigkeit 5 durch eine skalarwertige Funktion ( in der
Form
5 grad (
(13.3)
ersetzt werden kann. Setzen wir diese in die Kontinuitätsgleichung ein, dann bekommt man
die Laplacegleichung
div grad ( ( ( (
(13.4)
die der homogene Spezialfall der Poissongleichung ist.
13.2.2 Druck-Korrektur-Verfahren
Eine Poissongleichung für den Druck wird bei Druck-Korrektur-Verfahren gelöst, welche die
Navier-Stokes-Gleichungen unter Anwendung des Operator-Splittings entkoppeln.
Dazu wird der Druck zum unbekannten Zeitschritt in der Form
geschrieben. Dabei werden die Navier-Stokes-Gleichungen in zwei Gleichungen
5
5
5 grad 5 5 5 5
grad 5
grad zerlegt. Die Addition beider Gleichungen ergibt die ursprünglichen Navier-StokesGleichungen, wobei der Druckterm zum unbekannten Zeitschritt berechnet wird.
Der Trick des Verfahrens besteht darin, daß man für einen bekannten Druck einsetzt. Dieser
kann entweder
der hydrodynamische Druck zum Zeitschritt oder der hydrostatische Druck [6] sein.
In beiden Fällen verschwindet der Druck als gesuchte Variable aus den Navier-StokesGleichungen.
Von der zweiten Teilschrittgleichung wird die Divergenz gebildet. Da man div 5 erreichen will, verbleibt die Lösung der Poissongleichung
div 5 (13.5)
Im Gegensatz zur exakten Druck-Poissongleichung taucht in ihr nur der Druck als gesuchte
Variable auf, die rechte Seite ist nach der Lösung der ersten Schrittes bekannt.
Druck-Korrektur-Verfahren existieren in den verschiedensten Variationen für stationäre als
auch instationäre Probleme. Zu nennen sind u.a. die Familie der SIMPLE-Algorithmen [21].
KAPITEL 13. GALERKINVERFAHREN
156
13.3 Lösung der 1D-Poissongleichung
Als einfaches Beispiel für die Anwendung des Galerkinverfahrens wollen wir nun die Poissongleichung im Eindimensionalen
numerisch lösen. Der Leser studiere diesen Abschnitt sehr genau, nur dann wird ihm die Anwendung des Standard-Galerkin-Verfahrens in Fleisch und Blute übergehen.
Im ersten Schritt werden die Funktion und durch Lagrangesche Ansatzfunktionen ( auf
den Knoten 0 interpoliert, d.h.
(
(
Man erhält:
( (
Offensichtlich sind nun die Ableitungen auf die Ansatzfunktionen übertragen worden.
Im zweiten Schritt multiplizieren wir mit den Wichtungsfunktionen, die im StandardGalerkinverfahren die Ansatzfunktionen ( selbst sind:
(
(
((
1 Im dritten Schritt integrieren wir über das gesamte Lösungsgebiet :
(
( ( ( 1 Man sieht hier, daß zur Lösung der Poissongleichung mindestens Ansatzfunktionen zweiter
Ordnung vonnöten sind, da ansonsten die linke Seite verschwinden würde. Allgemein folgern
wir, das die Ordnung einer Differentialgleichung (d.h. die maximale Ableitungsordnung) auch
die Ordnung der Ansatzfunktion bestimmt. Hier läßt sich allerdings als Trick die partielle
Integration anwenden, der die Anwendung linearer Ansatzfunktionen ermöglicht. Für die linke
Seite gilt
Somit erhalten wir:
( (
( ( ( (
13.3. LÖSUNG DER 1D-POISSONGLEICHUNG
( ( Wir bestimmen zuerst den Term
157
( (
( ( ( ( 1 ( ( (
für
1
für
1
Durch die partielle Integration werden also Neumannsche Randbedingungen am ersten und
letzten Knoten automatisch in das Verfahren eingebunden. Für alle Innenknoten reduzieren
sich die Gleichungen zu:
( (
( ( 1 Dieses lineare Gleichungssystem wird übersichtlicher, wenn man zur Abkürzung die Koeffizienten
! ( ( (13.6)
( (
(13.7)
der Massenmatrix & und die Koeffizienten
der Steifigkeits- oder Diffusionsmatrix ; einführt. Letztere hat ihre Bezeichnung in der Vorwegnahme kommender Dinge. Man erhält das sehr einfach aussehende Gleichungssystem
; &
Im vierten Schritt müssen die Koeffizientenmatrizen berechnet werden, was im wesentlichen
auf die Auswertung verschiedener Integrale über Kombinationen der Ansatzfunktionen hinausläuft.
Die Massenmatrix ist symmetrisch, d.h. es gilt ! ! . Desweiteren sind nur die Koeffizienten von Null verschieden, bei denen beide Ansatzfunktionen einen nichtverschwindenden
Überlapp haben. Es verbleibt somit, ! , ! und ! zu berechnen. Für diese folgt nach
langer ermüdlicher Rechnung
KAPITEL 13. GALERKINVERFAHREN
158
! und
!
! Die Steifigkeitsmatrix ist ebenfalls symmetrisch und tridiagonal, mit
( für
für
für
und
ergibt sich sofort
Insgesamt ergibt sich somit an Innenknoten 1
das Verfahren
Somit nimmt auch das FE-Verfahren auf einem äquidistanten Gitter die äußerliche Form eines
FD-Verfahrens an. Die Inhomogenität taucht allerdings in der Form einer räumlichen Wichtung auf. Das so entstandene lineare Gleichungssystem ist nun mit einem geeigneten Löser zu
lösen.
13.4 Integrationsverfahren
Wir wollen nun von dem Spezialfall einer äquidistanten Diskretisierung absehen und zudem auch mehrdimensionale Probleme betrachten. Um das Gleichungssystem der FEDiskretisierung aufzubauen, sind im wesentlichen verschiedene Integrale aus Kombinationen
der Ansatz- und Wichtungsfunktionen und deren Ableitungen auszuwerten.
13.4.1 Die Integraltransformationsformel
Dabei ist auch hier die Transformation auf ein Referenzelement eine wichtige Rechentechnik.
Ist das Element, über welches integriert werden soll und das Referenzelement, so gilt
der Integraltransformationssatz:
wobei C1 die Determinante der Jakobimatrix
8 : < C1 (13.8)
13.4. INTEGRATIONSVERFAHREN
C1
159
2'
3
'
2
3
'
2
3
der Transformation ist.
13.4.2 Analytische Integration
Optimal sind analytische Integrationsbeziehungen, die bei einfachen Elementgeometrien etwa
im 1D-Fall, bei 2D-Dreiecken oder 3D-Tetraedern existieren. So hilft bei der Integration des
Einheitslinienelementes die Formel
8
8 8 (13.9)
Auch auf dem 2D-Einheitsdreieck läßt sich jedes Polynom mittels
8 : (13.10)
integrieren und für das Einheitsquadrat gilt
8 : (13.11)
wobei hier jedoch zu beachten ist, daß bei allgemeinen Vierecken die Jakobideterminante nicht
konstant ist.
13.4.3 Numerische Integration
Bei komplexeren Elementen muß man auf numerische Integrationsverfahren ausweichen. Als
numerische Integrations- oder Quadraturformeln bezeichnet man Ausdrücke der Art (im 2D):
(13.12)
Die Punkte bezeichnet man als Stützstellen und die Vorfaktoren als Gewichte der
Quadraturformel. Der wichtigste Begriff der numerischen Integrationstheorie ist der des Exaktheitsgrades, welcher die Ordnung des Polynoms ist, für die die Quadraturformel exakt ist.
So besitzt eine Integrationsformel den Exaktheitsgrad 2, wenn sie für alle Polynome der Form
den exakten Integralwert liefert. Somit ist zuerst zu entscheiden, wie genau die Integrationsformel sein muß; im Beispiel der Transportgleichung treten Produkte der Ansatzfunktionen in
der Massenmatrix auf. Sind die Ansatzfunktionen linear, dann sollte der Exaktheitsgrad der
Integrationsformel zwei sein.
KAPITEL 13. GALERKINVERFAHREN
160
Über den Begriff des Exaktheitsgrades wird auch die Konstrukionsweise von Quadraturformeln deutlich: Man bestimmt das exakte Integral eines Polynoms gegebener Ordnung und
hieraus die Gewichte bzw. Stützstellen.
Als Beispiel sei eine Integrationsformel des Exaktheitsgrads 1 für das Einheitsdreieck gesucht.
Das Polynom ersten Grades ist
und für das Integral über das Einheitsdreieck gilt exakt
woraus man leicht die Integrationsformel
abliest. Integrationsformeln, die für einen gegebenen Exaktheitsgrad mit der geringsten Anzahl
an Stützstellen auskommen, nennt man Gauß-Quadraturformeln. Diese sind für verschiedene
Elementtypen und für verschiedene Exaktheitsgrade in der Literatur tabelliert.
Das Durchführen der Integration und das Aufstellen des algebraischen Gleichungssystems
ist eine Ursache des wesentlich höheren Rechenaufwandes der FE gegenüber der FD oder
FV-Methode, bei denen dieser Rechenschritt entfällt. Desweiteren ist der Rechenaufwand der
FE-Methode deswegen wesentlich höher, weil die entstehenden Matrizen bei ungünstiger Gitternummerierung oftmals voll besetzt sind, während FDM und FVM Bandmatrizen hervorbringen.
13.5 Die Poissongleichung auf Dreiecken
Die bisherige Darstellungen zum Galerkinverfahren haben weniger mit Finiten Elementen als
vielmehr mit der Herleitung eines neuen Finite Differenzen Verfahrens zu tun. Sie haben aber
gezeigt, wie man mit dem Galerkinverfahren zu brauchbaren Diskretisierungen kommt. Wir
wollen uns nun den eigentlichen Finiten Elementen, die im zweidimensionalen oftmals Dreiecke sind, zuwenden.
Wir gehen davon aus, daß das Lösungsgebiet durch + Dreiecke trianguliert ist, die von
insgesamt Knoten aufgespannt werden. Ein Beispiel für eine ganz einfache Triangulierung
ist in Abbildung 13.1 dargestellt. Die Topologie eines solchen Dreiecksgitters bildet man im
sogenannten Elementverzeichnis ab, welches jedem Element einer gewissen Schlüsselnummer
die zugehörigen Knoten zuordnet. Das Knotenverzeichnis enthält dagegen die Knotennummer
als Schlüssel und ordnet diesen ihre - und -Koordinaten zu. Hilfreich ist es jedoch auch, hier
die angrenzenden Elemente zu verzeichnen.
Bei der Entwicklung der Diskretisierung gibt es nun zwei vollkommen analoge Betrachtungsweisen. Wir hatten im vorhergehenden Abschnitt die knotenbezogene Betrachtungsweise
angewendet, nun wollen wir es mit der elementbezogenen Betrachtungsweise versuchen.
Unser Testbeispiel sei der Einfachheit halber die Poissongleichung im zweidimensionalen:
13.5. DIE POISSONGLEICHUNG AUF DREIECKEN
161
E4
E2
E1
E3
E6
E5
!
Abbildung 13.1: Eine einfache Beispieltriangulation.
Elementschlüssel Knoten 1
1
1
2
1
3
1
4
2
5
1
6
1
Knoten 2 Knoten 3
7
4
4
2
2
3
6
3
3
5
5
7
Knotenschlüssel x-Koodinate y-Koordinate Anzahl der Patchelemente Elemente
5
1,2,3,5,6
1
3
2,3,4
2
3
3
3,4,5
2
1,2
4
2
5,6
5
6
1
4
2
1,6
7
Tabelle 13.1: Element- und Knotenverzeichnis zu der einfachen Beispieltriangulation.
KAPITEL 13. GALERKINVERFAHREN
162
Wir erinnern uns daran, daß wir das allgemeine Interpolationsproblem dadurch gelöst haben,
daß das Gesamtgebiet in Finite Elemente zerlegt wurde und die Interpolation nur auf einem
solchen Teilelement durchgeführt wurde (siehe Kapitel 6). Genauso bauen wir nun die Ansatzfunktion aus den einzelnen Dreiecken zusammen, setzen also
falls (x,y) in Dreieck k
wobei die Funktion im Dreieck die dortige Interpolationsfunktion und ansonsten
Null ist. Auf dem Dreieck ist die Ansatzfunktion durch die entsprechende lineare Interpola-
tionsfunktion
( ( ( definiert, wobei die Koeffizienten durch die Gleichungen (6.11) und (6.10) gegeben sind. Die
Gewichte , und beziehen sich dabei auf die drei Knoten des -ten Dreieckes.
Wir setzen zunächst die Ansatzfunktionen für die gesuchte Lösung und die Inhomogenität in die zu lösende Differentialgleichung und bekommen
div grad ( (
In jedem Dreieck bzw. Finitem Element gibt es nur drei Unbekannte, also müssen wir durch
Multiplikation mit den drei Ansatzfunktionen des Dreieckes drei Gleichungen erzeugen:
( div grad ( Die Integration über das Dreieck führt zu:
( div grad ( ( ( ( ( Spätestens an dieser Stelle sollten wir uns wieder Gedanken darüber machen, daß die Divergenz des Gradienten, also zweite Ableitungen der linearen Ansatzfunktionen Null sind. Somit
verschwindet die ganze rechte Seite. Glücklicherweise steht uns aber auch hier das Mittel der
partiellen Integration in Form der 1. Greenschen Formel
( div grad ( ( grad ( A5 grad ( grad ( zur Verfügung, so das alle Ansatzfunktionen höchstens einmal abgeleitet werden. Unser nun
schon algebraisches Gleichungssystem bekommt die Gestalt:
(
)
( grad ( A5 grad ( grad ( *
+ ( ( 13.5. DIE POISSONGLEICHUNG AUF DREIECKEN
163
Dieses recht komplexe mathematische Objekt wollen wir detaillierter betrachten. Die linke
Seite enthält wieder die Massenmatrix & , bloß daß diese nun auf ein Dreieck beschränkt
ist. Sie wird mit dem auf das Dreieck beschränkte Inhomogenität multipliziert:
(
)
( grad ( A5 *
+
grad ( grad ( & Die linke Seite besteht im hinteren Teil aus einer Matrix-Vektor-Multiplikation mit der
Diffusions- oder Steifigkeitsmatrix ; :
( grad ( A5 ; & Ist das Dreieck durch die Knoten , und aufgespannt, dann bekommt
man nach kurzweiliger Rechnung für die Elemente der Massenmatrix
& ( ( und für die Elemente der Diffusionsmatrix
;
( (
;
;
;
( (
; Die linke Seite beginnt mit einem Integral über den Rand des Dreiecks . Der Randnorma5 steht dabei senkrecht auf dem Rand und weist aus dem Dreieck hinaus. Physilenvektor A
kalisch beschreibt dieser Term den Austausch zwischen benachbarten Dreiecken, ist also für
die Kommunikation zwischen letzteren verantwortlich. Er ist also extrem wichtig und damit
wir an dieser Stelle keine Fehler machen, entledigen wir uns seiner Berechnung durch Summation über alle Elemente. Dabei stossen an jeder Innenkante zwei Elemente zusammen, die
denselben Kommunikationsterm mit umgekehrten Vorzeichen haben, weil die Randnormalenvektoren antiparallel zueinander orientiert sind. Der Randintegralterm bleibt damit nur noch
an den Aussenkanten des Gesamtgebietes stehen:
( grad ( A5 ; & KAPITEL 13. GALERKINVERFAHREN
164
Diese Summation über alle Elemente ist algorithmisch nicht trivial und soll hier nur skizziert
werden. Zunächst erweitern wir die 3x3-Gleichungssysteme zu NxN-Gleichungssystemen, in
dem die zu einem Element gehörigen Knoten an die Stellen ihrer Globalnummerierung eingefügt werden und der Rest der Matrizen mit Nullen aufgefüllt werden. Das MassenmatrixVektor-Produkt für Element 1 aus Abbildung 13.1 hatte dann die folgende Gestalt:
& ; D Diese Gleichungen können dann im zweiten Schritt summiert werden, da sie nun diesselbe
Gestalt haben. Diesen Prozess bezeichnet man als Kompilation der Systemmatrix.
Schließlich muss das so erhaltene algebraische Gleichungssystem gelöst werden.
13.6 Zusammenfassung
Das Standard-Galerkinverfahren zur Lösung einer Differentialgleichung besteht aus vier Teilschritten. Zunächst wird die gesuchte Funktion durch eine Ansatz- oder Testfunktion ersetzt,
die in der Regel aus einer Linearkombination von Interpolationspolynomen besteht. Um soviele Gleichungen wie Unbekannte zu erzeugen, wird im zweiten Schritt die Ausgangsgleichung mit jeder Interpolationsfunktion multipliziert. Dieses Differentialgleichungssystem besteht nun aus genauso vielen Gleichungen wie Unbekannten und wird nun durch Integration
über den Lösungsraum in ein algebraisches Gleichungssystem überführt. Der letzte Schritt
beinhaltet die Lösung dieses Systems.
Kapitel 14
Zeitabhängige Finite-Elemente-Methoden
Die Probleme der Strömungsmechanik sind zumeist zeitabhängig, wir hatten bei der Vorstellung des Galerkinverfahrens bisher nur mit stationären Problemen zu kämpfen. Um die
Zeitabhängigkeit in das Lösungsverfahren zu bekommen gibt es prinzipiell zwei Möglichkeiten. Man wählt eine beliebiges Zeitschrittverfahren, so wie wir diese in Kapitel 2 kennengelernt haben. Die Lösungsfunktion liegt nun diskret auf den in das Verfahren eingehenden
Zeitschritten vor und wird jeweils durch ortsabhängige Ansatzfunktionen approximiert. Dieses Verfahren drückt sich um eine Finite-Elemente-Darstellung der Zeitabhängigkeit, es soll
im folgenden Abschnitt anhand der eindimensionalen Transportgleichung demonstriert werden.
Die zweite Klasse von FE-Verfahren für zeitabhängige Probleme besteht in der Verwendung
orts- und zeitabhängiger Ansatzfunktionen.
14.1 Lösung der 1D-Transportgleichung
Als einfaches Beispiel ziehen wir mal wieder die eindimensionale Transportgleichung
heran.
Die Zeitableitung diskretisieren wir mit einem expliziten Eulerverfahren, wobei es dem Leser
überlassen bleibt, die Methodik auf andere Zeitintegrationsschemata auszudehnen.
Im ersten Schritt werden die Funktionen und nun durch Lagrangesche Ansatzfunktionen ( auf den Knoten 0 interpoliert, d.h.
(
(
Man erhält:
165
KAPITEL 14. ZEITABHÄNGIGE FINITE-ELEMENTE-METHODEN
166
( ( ( Offensichtlich sind nun alle Ableitungen auf die Ansatzfunktionen übertragen worden.
Im zweiten Schritt multiplizieren wir mit den Wichtungsfunktionen, die im StandardGalerkinverfahren die Ansatzfunktionen ( selbst sind:
( ( ( (
(
(
1 Im dritten Schritt integrieren wir über das gesamte Lösungsgebiet :
( ( (
( (
( Genau wie bei der Poissongleichung wird der Diffusionsterm auch hier partiell integriert, damit man mit linearen Ansatzfunktionen arbeiten kann:
( (
(
( ( ( Somit erhalten wir:
(
( ( ( ( ( (( Wir bestimmen zuerst den Term
( ( (
( ( für
1
für
1
Durch die partielle Integration werden also Neumannsche Randbedingungen am ersten und
letzten Knoten automatisch in das Verfahren eingebunden. Für alle Innenknoten reduzieren
sich die Gleichungen zu:
( ( (
( ( (
14.1. LÖSUNG DER 1D-TRANSPORTGLEICHUNG
167
Dieses lineare Gleichungssystem wird übersichtlicher, wenn man zur Abkürzung die Massenmatrix nach Gleichung (13.6) und die Steifigkeitsmatrix nach Gleichung (13.7) verwendet.
Zusätzlich sei hier die Advektionsmatrix mit den Koeffizienten
(
( (14.1)
eingeführt. Man erhält das sehr einfach aussehende Gleichungssystem
& & ;
Im vierten Schritt müssen die Koeffizientenmatrizen berechnet werden, was im wesentlichen
auf die Auswertung verschiedener Integrale über Kombinationen der Ansatzfunktionen hinausläuft.
Die Massenmatrix ist symmetrisch, d.h. es gilt ! ! . Desweiteren sind nur die Koeffizienten von Null verschieden, bei denen beide Ansatzfunktionen einen nichtverschwindenden
Überlapp haben. Es verbleibt somit, ! , ! und ! zu berechnen. Für diese folgt nach
langer ermüdlicher Rechnung
! und
!
! Die Diffusionsmatrix ist ebenfalls symmetrisch und tridiagonal, mit
( für
für
für
und
ergibt sich sofort
Die Advektionsmatrix ist nicht symmetrisch, aber ebenfalls tridiagonal:
Insgesamt ergibt sich somit an Innenknoten 1
" #
das Verfahren
168
KAPITEL 14. ZEITABHÄNGIGE FINITE-ELEMENTE-METHODEN
Somit nimmt auch das FE-Verfahren auf einem äquidistanten Gitter die äußerliche Form eines
FD-Verfahrens an. Ein wesentlicher Unterschied liegt jedoch in der Darstellung der Zeitableitung, die auf beiden Zeitebenen in Form einer räumlichen Wichtung auftaucht. Damit ist das
Verfahren letztlich nicht explizit, obwohl alle Terme auf der bekannten Zeitebene diskretisiert wurden. Wir erhalten somit beim reinen Standard-Galerkin-Verfahren immer implizite
Schemata und müssen folglich ein Gleichungssystem lösen.
Mass-lumping
Mit Mass-lumping-Techniken wird die Massenmatrix diagonalisiert, so daß auch explizite Verfahren konstruiert werden können. Die einfachste Technik ist das Zeilensummenverfahren: Dabei ersetzt man in der Massenmatrix die Diagonalterme durch die Summe der Terme in einer
Zeile. Die Massenmatrix wird so zu
! Æ und es ergibt sich das altbekannte zentrale explizite Verfahren
Mass-lumping kann auch partiell angewendet werden, indem zwischen der ursprünglichen und
der diagonalisierten Matrix kontinuierlich gewichtet wird.
Mass-lumping ist ein Trick zur Konstruktion expliziter Verfahren, welches die Minimaleigenschaft des Residuums zerstört. Man sollte sich damit abfinden, daß bei der Anwendung der
FEM immer ein implizites Gleichungssystem zu lösen ist.
14.2 Die Diskretisierung der Zeit durch Finite Elemente
Die mit der Methode der Finiten Elemente konsistente Diskretisierung der Zeit erfolgt durch
die Einführung zeitabhängiger Ansatz- und Wichtungsfunktionen. Um das Wesentliche dieser
Vorgehensweise herauszuarbeiten, sei wieder die allgemeine homogene Evolutinsgleichung
betrachtet, wobei wir uns um die Ortsabhängigkeit der Lösung und des Operators nicht
scheren wollen.
Die gesuchte Lösung wird im ersten Schritt als Lagrangesche Interpolationsfunktion zwischen
den Zeitschritten und angeschrieben:
( (
Die Ortsabhängigkeit des Problems verbirgt sich nun in den Funktionen und . Wir
wollen zwischen den beiden Zeitschritten lineare Ansatzfunktionen verwenden:
( und
( 14.3. SCHWACHE LAGRANGEVERFAHREN
169
Man bezeichnet eine solche Formulierung dann als konsistent mit der Diskretisierung des
Raumes, wenn beide Ansatzfunktionen gleich geartet, also z.B. linear sind. Setzen wir diese Ansatzfunktion also in die Evolutionsgleichung ein:
( ( Wieder soll mit den Ansatzfunktionen auch gewichtet werden. Es entstünden so zwei Gleichungen, wobei allerdings nur der Wert unbekannt ist. Daher brauchen wir nur mit einer
Funktion zu wichten, wir entscheiden uns zunächst für ( und integrieren gleich zeitlich
zwischen und :
( ( (
((
Wieder können die Integrale analytisch ausgewertet werden. Das Ergebnis führt auf die zeitdiskretisierte Gleichung:
Die mit linearen Ansatzfunktionen konsistente Diskretisierung der Zeit schlägt also einen
Crank-Nicolson-Faktor von vor.
Genauso hätte man nun aber auch die Ansatzfunktion ( als Wichtungsfunktion verwenden
können. In diesem Fall lautete die zeitdiskretisierte Evolutionsgleichung:
Läßt man nun in Verallgemeinerung dessen beliebige Linearkombinationen der Ansatzfunktionen ( und ( als Wichtungsfunktionen zu, dann ist zu folgern, daß der Crank-NicolsonFaktor im Intervall [1/3, 2/3] liegen sollte. Eine rein explizite oder implizite Diskretisierung
der Zeit und nachfolgende Finite-Elemente-Verfahren für die Raumabhängigkeit führen also
zu Inkonsistenzen.
14.3 Schwache Lagrangeverfahren
Wir wollen uns nun den advektiven Termen bzw. der Lagrangeschen Zeitableitung zuwenden.
Zurückblickend hatten sich zu deren Lösung die Lagrangeverfahren in einem gewissen Sinne
als ideal erwiesen, sie waren stabil, weil sie sich stark an der Physik der Advektion orientierten.
Die bisher untersuchten Lagrangeverfahren hatten allerdings zwei Schwachstellen. Zum einen
führt die lineare Interpolation an der Basis der Charakteristiken zu einer gewissen numerischen
Diffusion, die jedoch bei der Verwendung von Interpolationsverfahren zweiter Ordnung weitestgehend unterdrückt werden konnte. Schwierigkeiten bereiten solche Verfahren allerdings
immer noch auf unstrukturierten Gittern. Zum anderen waren die Verfahren nicht massenerhaltend, da die Advektionsgleichung selbst nicht konservativ ist.
Das nun vorgestellte schwache Lagrangeverfahren hat beide Schwächen nicht. Die Konservativität ist dadurch gewährleistet, daß das Verfahren auf der konservativen Form der Advektionsgleichung
KAPITEL 14. ZEITABHÄNGIGE FINITE-ELEMENTE-METHODEN
170
div
5 aufbaut und somit massenerhaltend ist [12], [13].
Die schwache Formulierung wird durch die Multiplikation mit einer orts- und zeitabhängigen
Ansatzfunktion E und entsprechender Integration konstruiert:
div
5 E Die partielle Integration über die Zeit und Anwendung des Greenschen Integralsatzes über den
Ort liefert
E E
E 5gradE E5A
Die Ansatzfunktionen E werden nun so gewählt, daß
E
5
grad E (14.2)
Dies ist unproblematisch, wenn man die nur ortsabhängigen Lagrangeschen Ansatzfunktionen
( zum Zeitschritt ansetzt und die Funktionen zum Zeitschritt als Lösung von Gl. (14.2)
mit dem starken Charakteristikenverfahren bestimmt. Es verbleibt dann das Gleichungssystem
E E E5A
(14.3)
zu lösen. Dies ist jedoch i.a. in höchstem Maße rechenaufwendig, da die advektierten Ansatzfunktionen keine Lagrangefunktionen sind und die Auswertung der Integrale nur numerisch
erfolgen kann.
Für das schwache Lagrangeverfahren läßt sich im 1D-Fall ein Schema unter der Voraussetzung
einer konstanten Courant-Zahl kleiner Eins bestimmen [11]. Unter Einbezug des Diffusionsschrittes ergibt sich dieses als
, -
Das Ergebnis für den Transport einer steilen Front in Abb. 14.1 ist überragend. Es zeigt eine
verschwindende numerische Diffusion und keine mit der Pecletzahl verbundenen Oszillationen.
14.4. UPSTREAM-STRATEGIEN: PETROV-GALERKIN-VERFAHREN
171
1 .0 0
0 .7 5
0 .5 0
0 .2 5
0 .0 0
0 .0 0
0 .2 5
0 .5 0
0 .7 5
1 .0 0
Abbildung 14.1: Simulation einer steilen Front mit dem schwachen Lagrangeverfahren, blau:
analytische, rot: numerische Lösung.
14.4 Upstream-Strategien: Petrov-Galerkin-Verfahren
Da die schwachen Lagrangeverfahren bisher noch nicht auf unregelmäßigen Geometrien implementiert und getestet wurden, weil der Aufwand hierfür sehr groß ist, muß auf andere Verfahren ausgewichen werden, um das Upstream-Prinzip bei der Behandlung der advektiven
Terme nicht zu verletzen.
Der Upstream-Effekt wird durch die Verwendung modifizierter Wichtungsfunktionen erzielt.
Diese sind so geformt, daß in das Schema letztlich mehr Informationen aus der Gegenstromrichtung eingeht.
Ganz allgemein bezeichnet man Verfahren, bei denen die Wichtungsfunktionen von den
Ansatzfunktionen ( verschieden sind als Petrov-Galerkin-Verfahren. Gl. (13.2) wird dann zu:
$ $
%
(
%
(14.4)
14.4.1 Das Verfahren nach Christie
Eine gegen die Strömungsrichtung verzerrte Wichtungsfunktion ist z.B. für den 1D-Fall durch
[8]:
( 2
( 2 für
für
gegeben, wobei sich für 2 die lineare Ansatzfunktion ergibt. Ansatz- und Wichtungsfunktion nach Christie sind in Abb. 14.2 dargestellt.
172
KAPITEL 14. ZEITABHÄNGIGE FINITE-ELEMENTE-METHODEN
14.4.2 Das Streamline Upwind/Petrov-Galerkin (SU/PG) Verfahren
Ein auf beliebige mehrdimensionale Elemente verallgemeinerungsfähiges Upwind-Verfahren
ist das Streamline Upwind/Petrov-Galerkin (SU/PG) Verfahren von Brooks und Hughes [3].
Hier werden die Wichtungsfunktionen nur in der jeweiligen Strömungsrichtung durch den
Ansatz
Ü ( Ü 2
5
grad ( Ü
5
(14.5)
verändert, wobei 2 der Upstream-Parameter ist. Für den linearen 1D-Fall sind die entsprechenden Wichtungsfunktionen in Abb. 14.2 dargestellt.
Zur Funktionsweise des SU/PG-Verfahrens vergleichen wir die Anwendung des StandardGalerkin-Verfahrens auf den advektiven Term
5 grad ( mit der Anwendung des SU/PG-Verfahrens:
5 grad ( 2
5
grad grad ( 5
dann erscheint ein zusätzlicher Term, der die Form eines Diffusionstermes hat. Das SU/PGVerfahren stabilisiert die Advektion also durch die Zufügung einer künstlichen numerischen
Diffusion in Strömungsrichtung. Die numerische Diffusion ist im eindimensionalen Fall durch
2
gegeben.
Für die Wahl des Upwind-Parameters 2 existiert bisher keine geschlossene Theorie. Sie bestimmt jedoch entscheidend die Stabilität des Verfahrens als auch die numerische Diffusion.
Brooks und Hughes [3] selbst schlagen
2
mit einer numerischen Diffusion von
*/
vor, Hervouet und Moulin [14] erhalten mit Hilfe der Fourier-Analyse ein Maximum an Stabilität für
2 */
mit einer numerischen Diffusion von
*/
*/ Offensichtlich ist das Verfahren nach Brooks und Hughes für Courantzahlen größer eins und
das Verfahren nach Hervouet für Courantzahlen kleiner eins weniger diffusiv.
14.5. ZUSAMMENFASSUNG
173
i-1
i
i+1
i-1
i
i+1
Abbildung 14.2: Upwind-Wichtungsfunktionen nach Christie et al. und SU/PG
14.5 Zusammenfassung
In der Praxis lassen sich zeitabhängige Probleme bei der Verwendung von räumlich linearen
Ansatzfunktionen durch eine Vorschaltung des Crank-Nicolson-Verfahrens zur Zeitdiskretisierung lösen. Der Crank-Nicolson-Faktor sollte dann nicht kleiner als 1/3 und nicht größer als
2/3 sein.
Bei der Advektion zeigt das schwache Lagrangeverfahren exzellente Eigenschaften bei eindimensionalen Problemen. Die Ausarbeitung und Übertragung des Algorithmus auf mehrdimensionale Probleme ist daher ein lohnendes Forschungsfeld.
174
KAPITEL 14. ZEITABHÄNGIGE FINITE-ELEMENTE-METHODEN
Kapitel 15
Finite Elemente und Funktionalanalysis
Nun gehts ans Eingemachte: Um zu verstehen, warum die FEM funktioniert, müssen wir tief
in den Hexenkessel der Funktionalanalysis eintauchen.
15.1 Funktionenräume
Die Funktionalanalysis behandelt Funktionen als eine Form von Vektoren. Grundlage hierfür
ist die Tatsache, daß Funktionen mit gewissen Eigenschaften Vektorräume bilden. So ist z.B.
jede Linearkombination von stetigen Funktionen wieder eine stetige Funktion, den so erzeugten Vektorraum der stetigen Funktionen auf einem Gebiet bezeichnet man mit * .
Wir wollen durch den Index 0 die Menge aller Funktionen bezeichnen, die auf dem Rand
eines Gebietes Null sind. * bezeichnet also die auf stetigen Funktionen, die
auf dem Gebietsrand Null sind. Da beliebige Linearkombinationen solcher Funktionen diese
Eigenschaften wieder erfüllen, bilden solche Funktionen ebenfalls einen Vektorraum.
Auch Funktionenräume werden durch Basen aufgespannt. So besagt der Satz von Fourier
nichts anderes, als daß die trigonometrischen Polynome
eine Basis des Vektorraumes * 7 bildet, da sich jeder Vektor (d.h. zwischen und 7
stetige Funktion) durch eine Reihe
d.h. durch eine unendliche Linearkombination von trigonometrischen Polynomen, darstellen
läßt. Und der Satz von Taylor liefert uns eine Basis für den Raum der analytischen d.h. unendlich oft stetig differenzierbaren Funktionen * )=; sie lassen sich nämlich als Reihe
darstellen. Die Polynome
175
176
KAPITEL 15. FINITE ELEMENTE UND FUNKTIONALANALYSIS
bilden somit eine Basis des Raumes der analytischen Funktionen.
Aber schon an dieser Stelle offenbart sich der wesentliche Unterschied zu den ’bekannten’
Vektorräumen der linearen Algebra: Funktionenräume sind unendlichdimensional, d.h. ihre
Basen bestehen aus unendlich vielen Funktionen.
Und hier offenbart sich die Notwendigkeit, numerische Analysis zu treiben. Ohne die Definition des Begriffes Algorithmus heranzuziehen, wird niemand die Zeit haben, einem unendlich
lange währenden Rechenprozeß zu folgen.
Damit werden wir geradezu gezwungen, die mathematische Sünde1 zu begehen, unendlichdimensionale Funktionenräume auf endliche zu reduzieren. Lasset uns konkret werden.
Bezüglich einer beliebigen Basis
( ( ( des Lösungsraumes läßt sich die gesuchte Lösung als Linearkombination
(
(15.1)
darstellen. Der Prozeß der Diskretisierung besteht darin, eine Teilbasis mit endlich vielen sogenannten Ansatzfunktionen auszuwählen und die Funktion durch die approximierte Funktion
(
(15.2)
zu ersetzen. Liest man diese Gleichung von links nach rechts, so drängen sich drei Fragen auf:
Wie gut ist diese Approximation der Lösung ?
Wie bestimmt man die Koeffizienten der Linearkombination und somit die approximierte Lösung ?
Wie wählt man den durch die Teilbasis repräsentierten endlichdimensionalen Funktionenraum sinnvoll aus ?
Bevor wir mit der Bearbeitung der dritten Frage beginnen, vertiefen wir unsere Kenntnisse
über einige spezielle Funktionenräume.
15.2 Hilberträume
Ein Hilbertraum ist eine algebraische Struktur; andere sind die Körper, Ringe oder Vektorräume. Algebraische Strukturen bestehen immer aus bestimmten Mengen, auf denen gewisse Rechenregeln erlaubt sind. Die Mathematiker beschäftigen sich mit solchen sehr allgemeinen und daher sehr abstrakten algebraischen Strukturen aus folgendem Grund: Sie untersuchen, welche Behauptungen und Sätze für diese allgemeine Struktur gelten und zeigen dann
1
Es gibt tatsächlich eine theologische Definition von Sünde als Abstand zum Transzendenten. Das Ziel der
Informatik, immer endliche Algorithmen zu konstruieren, ist somit die reine Sünde.
15.2. HILBERTRÄUME
177
im speziellen, daß gewisse Objekte der Mathematik diese algebraische Struktur aufweisen.
Damit erspart man sich eine Menge Beweisarbeit.
Wir werden diese Schlußweise gleich verstehen. Die algebraische Struktur, in der die Methode
der Finiten Elemente –und der größte Teil der Physik überhaupt– spielt, ist der Hilbertraum.
Def.: Ein Hilbertraum F besteht aus einem vollständigen Vektorraum mit einem Skalarprodukt.
Der Begriff der Vollständigkeit besagt formal, daß der Grenzwert jeder Cauchyfolge Element
des entsprechenden Raumes ist. So wichtig dieser Begriff auch ist – wir werden ihn hinfort
nicht benötigen.
Unter dem Skalarprodukt zweier Vektoren und versteht man eine (durchaus sehr
beliebige) Rechenvorschrift die den beiden Vektoren eine reelle Zahl zuordnet. Sie muß
allerdings folgende Eigenschaften erfüllen:
(SP1):
2 6 2 6 (SP2):
2 6 2 6 (SP3):
In einem Hilbertraum wird der geometrische Begriff des Senkrechten verallgemeinert:
Def.: Zwei Vektoren
gilt.
und heißen orthogonal oder senkrecht zueinander, wenn Wir wollen nun, wie schon angekündigt, die Denkweise der Algebra an folgendem Satz üben:
Satz: Ein Vektor ist genau dann der Nullvektor, wenn er zu jedem anderen Vektor des
Hilbertraumes orthogonal ist. Anders ausgedrückt: Der Nullvektor steht senkrecht auf dem
gesamten Hilbertraum.
15.2.1 Der Hilbertraum ¾ ´ªµ
Wir betrachten die Menge der Funktionen
)=
und
./
0
, Man nennt solche Funktionen quadratintegrabel. Sie bilden genauso wie die stetigen oder nfach differenzierbaren Funktionen einen Vektorraum, denn jede Linearkombination solcher
Funktionen ist wieder quadratintegrabel.
mit dem Skalarprodukt
Zusätzlich ist
> >
(15.3)
KAPITEL 15. FINITE ELEMENTE UND FUNKTIONALANALYSIS
178
ein Hilbertraum. Der Leser bestätige, daß das so definierte Skalarprodukt, welches ja eigentlich
ein Integral ist, auch tatsächlich die Eigenschaften eines Skalarprodukts erfüllt.
Mit dem eingeführten Skalarprodukt können nun Funktionen zueinander orthogonal sein. Und
ohne einen aufwendigen Beweis zu führen, sehen wir sofort, daß die Funktion genau dann
.
überall Null ist, wenn > für alle > 15.2.2 Schwache Ableitungen und Sobolevräume
Zur Definition von Sobolevräumen benötigen wir Multiindizes. Diese sind Indexvektoren der
Form
2
2 ... 2
wobei
2 2 Die partiellen Ableitungen einer Funktion bis zur Ordnung 2 schreibt man dann als
4 4
4 Die partielle Integration verwandelt das Produkt aus einer Funktion eins und einer abgeleiteten
Funktion zwei in die Differenz aus einem Randwertterm und dem Integral über das Produkt
der Funktion zwei und der Ableitung der Funktion eins. Ist eine der Funktionen auf dem Rand
des Integrationsgebietes Null, dann verschwindet dabei der Randwertterm. Dieser Satz gilt
auch für beliebige partielle Ableitungen:
4 4 4 4
4 4 4 wenn
*
Diese Eigenschaft kann man umgekehrt dafür verwenden, eine besondere Form der partiellen
Ableitung zu definieren, die sogenannte schwache Ableitung.
Def.: Die Funktion ; 4 heißt schwache Ableitung von , wenn sie die Bedingung
4
4 4
4 für alle * erfüllt.
Warum diese Verrenkungen ? Schwache Formulierungen nutzen ganz allgemein eine Eigenschaft des Integrationsbegriffes aus, das Verhalten der zu integrierenden Funktion auf Nullmengen zu übersehen. Solche Nullmengen sind z.B. einzelne Punkte im Eindimensionalen,
Linien im zweidimensionalen oder Flächen im Dreidimensionalen. Ist eine Funktion oder eine
ihrer Ableitungen auf einer solchen Nullmenge unstetig, dann ist sie nach der der klassischen
15.2. HILBERTRÄUME
179
Definition der Ableitung dort nicht differenzierbar, die schwache Ableitung exisitiert dahingegen schon.
Als Sobolevraum F bezeichnet man den Hilbertraum
F
1
;4 2 2
Die Sobolevräume F bestehen also aus Funktionen , von denen alle partiellen Ableitungen bis zur Ordnung quadratintegrabel sind. Wir sehen, daß diese sehr theoretisch aussehenden Gebilde lediglich Schreibarbeit ersparen: Hier wird einmal wieder heißer gekocht als
gegessen.
Das Skalarprodukt dieser Hilberträume ist durch
$ Jedem ist natürlich klar, daß
$ 4
;4 ist und zur Übung schreibe man einmal $ aus. Der immer noch nicht geschockte Leser
sieht sofort, daß
$ " Schließlich wollen wir mit F die Einschränkung des Raumes
bezeichnen, die auf dem Rand Null sind.
F
auf Funktionen
15.2.3 Operatoren in Hilberträumen
Die Operatorentheorie versucht, sehr komplexe Operatoren (etwa Differentialoperatoren) wie
die Systemmatrizen algebraischer Gleichungssysteme zu behandeln. Die Ergebnisse dieser
Theorie sehen daher oftmals sehr einfach und wohlbekannt aus, obwohl die Grundlagen der
Theorie doch sehr komplex sind.
Ein Operator ist eine Abbildung in einem Hilbertraum. In Hilberträumen, die Funktionen
enthalten, können diese etwa Differentialoperatoren sein. Man kann nun eine Norm für die
Operatoren definieren, indem man
(15.4)
setzt. Die Norm eines Operators gibt also so etwas wie die maximale Auswirkung des Operators auf den zugrundeliegenden Hilbertraum an. Ist die Norm nicht unendlich, so heißt der
Operator beschränkt. Für die Norm zweier hintereinandergeschalteter Operatoren und gilt die wichtige Beziehung
(15.5)
KAPITEL 15. FINITE ELEMENTE UND FUNKTIONALANALYSIS
180
Der nächste Satz ist ein mächtiges Werkzeug für die Analyse numerischer Verfahren.
Modifizierter Satz von Kellogg: Sei raumes und . Dann gilt für beliebiges ) )
für alle eines reellen Hilbert (15.6)
Beweis: Es gilt
) )
)
) ) ) Die Ungleichung ist erfüllt, wenn der Zähler kleiner als der Nenner ist, also nach kurzer
Rechnung:
bzw.
was der Fall ist, wenn positiv semidefinit und ist.
(15.7)
15.2.4 Eigenwerte linearer Operatoren
Analog dem Vorgehen in der Linearen Algebra versucht man auch für lineare Operatoren
Eigenwerte und Eigenvektoren zu finden, respektive die Gleichung
' (15.8)
zu lösen. Auf jeden durch einen Eigenvektor aufgespannten Unterraum wirkt der Operator
als Streckung um den Faktor ' .
Ideal ist die Eigenwerttheorie dann, wenn die Eigenvektoren ein Orthonormalsystem des
Lösungsraumes bilden. Jede potentielle Lösung läßt sich dann als Linearkombination
. (15.9)
. (15.10)
mit
darstellen. Der Operator nimmt dann die Form
.
(15.11)
15.3. DUALRÄUME
181
an. Dies ist insbesondere für sogenannte kompakt-symmetrische Operatoren möglich.
15.3 Dualräume
Keine lineare Abbildung ist so einfach gebaut, wie die, die einem Vektor eine reelle Zahl
zuordnet. Daher hat man diesen Abbildungen einen eigenen Namen gegeben, man bezeichet
sie als lineare Funktionale.
Als Dualraum B B )= bezeichnet man die Menge aller stetigen linearen Funktionale
von einem Vektorraum B in die reellen Zahlen. Da Linearkombinationen von linearen Funktionalen wieder lineare Funktionale sind, ist B ein Vektorraum.
Mit Hilfe der Operatornorm können wir den Dualraum zu einem normierten d.h. Banachraum
machen. Für ? B setzen wir:
?
5 ¼
? (15.12)
5
5 wird als Dualnorm bezeichnet.
Stellen wir uns die zunächst äußerst dumm erscheinende Frage nach der Anzahl der stetigen
linearen Funktionale auf einem Vektorraum. Antwort gibt der
¼
Darstellungssatz von Riesz: Zu jedem stetigen ?
? und
B existiert genau ein B , so daß
?
5 ¼
5
(15.13)
Somit hat der Dualraum ebensoviel Elemente wie der zugehörige Vektorraum. Desweiteren
haben wir eine einfache Konstruktionsvorschrift; alle stetigen linearen Funktionale lassen sich
mit Hilfe von Skalarprodukten konstruieren.
Den Dualraum zum Sobolevraum F bezeichnet man mit F oder leider etwas
schlampig als F . Auf F lautet die Dualnorm:
(
" (15.14)
15.4 Bilinearformen
Wenn wir eine Gleichung mit linearem Differentialoperator mit Hilfe des Galerkinverfahrens diskretisieren, bekommen wir am Ende Ausdrücke der Form:
Diese sind allerdings nicht notwendig symmetrisch und auch nicht positiv. Die von einem
Skalarprodukt eingeforderten Eigenschaften sind also zu streng, um diesen Ausdruck zu beschreiben.
182
KAPITEL 15. FINITE ELEMENTE UND FUNKTIONALANALYSIS
Genauso wie das Skalarprodukt ordnen Bilinearformen zwei Vektoren eine reelle Zahl zu,
diese muß jedoch nicht positiv und die Abbildung muß nicht symmetrisch sein:
Def: Sei ein Vektorraum B , die Vektoren und der Skalar 2 gegeben.
B B )= heißt Bilinearform, falls sie in beiden Komponenten linear ist. Es gilt
also:
2 2 2 2 2 2 Das bekannteste Beispiel für eine Bilinearform sind die quadratischen Formen
einer Matrix und
bei denen der Vektor zuerst mit
dann mit einem Vektor skalar
multipliziert wird. Man kann somit jeder linearen Abbildung dargestellt durch eine Matrix
eine Bilinearform durch das Skalarprodukt zuordnen.
Def: heißt stetig-beschränkt, falls ein */ $ existiert mit
*/ Für die quadratischen Formen ergibt sich mit Hilfe der Cauchy-Schwarzeschen Ungleichung
der Wert */ .
Wir wollen nun obigen Sachverhalt umkehren und untersuchen, ob jeder Bilinearform in
eindeutiger Weise ein linearer Operator zugeordnet ist. Dies funktioniert nur für die
stetig-beschränkten Bilinearformen:
Satz: Einer stetig-beschränkten Bilinearform ist in eineindeutiger Weise ein stetig-beschränkter linearer Operator zugeordnet mit
und für die Operatornorm gilt:
*/
Einer der wichtigsten Ergebnisse der Theorie der Bilinearformen und deren zugeordneten
Operatoren ist eine exakte Bedingung für die Lösbarkeit der Gleichung , d.h. für die
Existenz der Inversen :
Satz: Der inverse Operator existiert genau dann, wenn die Bedingungen
B B
$ B B
$ %
%
15.4. BILINEARFORMEN
183
oder die sogenannten Babuska-Bedingungen
$ B B %
B $ gelten.
Wir benötigen noch folgende wichtige
Def: Eine Bilinearform heißt V-elliptisch, falls ein * .
*. $ existiert mit:
für alle
B
Ganz offensichtlich ist *. , falls die Bilinearform ein Skalarprodukt und die Norm aus diesem kanonich durch erzeugt wird. Der Begriff der VElliptizität wird also erst für unsymmetrische Bilinearformen interessant. Stetig-beschränkte
V-elliptische Bilinearformen sind somit nach oben und in einem gewissen Sinne auch nach
unten beschränkt.
Der Hauptsatz für die Methode der Finiten Elemente ist der folgende
Satz: Die Bilinearform sei stetig und V-elliptisch. Dann hat das Variationsproblem
suche B mit genau eine Lösung der Form für alle B
in B
5
und es gilt die Abschätzung
*.
5
¼
(15.15)
Damit liegt das Standard-Galerkin-Verfahren offen vor uns. Wir brauchen lediglich den
Lösungsraum B durch einen endlichen Teilraum B ersetzen. Besitzt dieser die Basis ( 0 , so läßt sich die Lösung als
(
darstellen. Die Gleichung ist genau dann für alle erfüllt, wenn sie für
jeden Basisvektor ( in der Form ( ( erfüllt ist, denn schließlich läßt sich auch
jedes als Linearkombination der Basis darstellen.
Unter anderen Umständen läßt sich aber auch ein Minimierungsproblem definieren:
Satz: Die Bilinearform
blem
sei symmetrisch und V-elliptisch. Dann hat das Minimierungspro-
suche B , so daß minimal ist
ein eindeutiges Minimum, welches erfüllt.
KAPITEL 15. FINITE ELEMENTE UND FUNKTIONALANALYSIS
184
15.5 Variationsaufgaben
Wir wollen nun in einer kurzen Zusammenstellung untersuchen, für welche Probleme eine Variationsaufgabe existiert. Dabei versuchen wir uns wieder an die Navier-Stokes-Gleichungen
heranzuarbeiten. Wir werden allerdings nur bis zur Transportgleichung gelangen, da wir noch
keine nichtlinearen Probleme untersucht haben.
15.5.1 Die Poissongleichung
Die Poissongleichung
(15.16)
taucht in der Hydrodynamik als Druck-Poisson-Gleichung auf. Sie wird mit einer Wichtungsfunktion multipliziert und über das Gesamtgebiet integriert:
Setzt man voraus, daß auf dem Rand des Gebietes
ersten Greenschen Formel
Null ist, führt die Anwendung der
grad grad grad grad Mit der Bilinearform
grad 5%
(15.17)
zu
%
grad grad wurde das ursprüngliche Differentialgleichungsproblem in ein Variationsproblem überführt.
Ohne es zu beweisen, glauben wir, daß die Bilinearform beschränkt-stetig, V-elliptisch und
symmetrisch ist.
15.5.2 Die Helmholtzgleichung
ist das Resultat der zeitlichen Diskretisierung der Diffusionsgleichung:
(15.18)
Wieder wird mit einer Wichtungsfunktion multipliziert, die auf dem Rand des Lösungsgebietes Null ist. Dann wird integriert und die erste Greensche Formel angewendet. Es verbleibt
das Problem
15.5. VARIATIONSAUFGABEN
185
grad grad %
zu lösen. Die Bilinearform
grad grad ist stetig-beschränkt, symmetrisch und da
grad $ F
grad -elliptisch.
15.5.3 Die stationäre Transportgleichung
sei in der Form
5 grad (15.19)
mit der Advektionsgeschwindigkeit 5 und der Diffusivität untersucht. Der schon bekannte
Weg führt auf die Bilinearform
5 grad grad grad die stetig-beschränkt aber nicht mehr symmetrisch ist. Damit läßt sich für die stationäre Transportgleichung kein Minimierungsproblem formulieren. Da
5 grad und wenn F
,
5 grad %
folgt
und somit
5 grad 5 grad 5 grad grad grad $ Somit ist nur noch F -elliptisch, d.h. die Lösungsfunktion muß auf dem Rand Null sein,
was einer homogenen Dirichletrandbedingung entspricht. Der nach der zeitlichen Diskretisierung verbleibende Rest der Transportgleichung scheint sich ein gegen die Anwendung der
FE-Methode zu wehren.
KAPITEL 15. FINITE ELEMENTE UND FUNKTIONALANALYSIS
186
15.6 Gemischte Finite Elemente
Auf einem Dreiecksgitter lassen sich mit den drei Eckknoten Ansatzfunktionen erster Ordnung und wenn man z.B. die Kantenmitten noch hinzunimmt, Ansatzfunktionen zweiter Ordnung definieren. Als gemische FE-Methoden bezeichnet man solche, bei denen verschiedene
Variablen durch Ansatzfunktionen verschiedener Ordnung dargestellt werden. Das klassische
Beispiel hierzu ist
15.6.1 Das Stokes-Problem
Es beschreibt eine stationäre advektionslose inkrompressible Strömung:
grad 5 5
(15.20)
div 5 Ein entsprechendes Problem tritt aber auch nach der Anwendung lagrangescher Verfahren auf,
womit das Stokes-Problem nicht nur akademischen Charakter hat.
Um so allgemein wie möglich zu bleiben, multiplizieren wir die Impulsgleichungen mit
5 * , die Kontinuitätsgleichung mit @ , integrieren über , wobei Druck
und viskose Terme partiell integriert werden:
div 5 grad 5 grad 5 5
5 @ div 5 Man beachte, daß nun die vektoriellen Impulsgleichungen zu einer skalaren Gleichung zusammengeschrumpft sind. Führen wir die Skalarprodukte
grad 5 grad 5 div 5 ein, so bekommt das Stokes-Problem in der Variationsformulierung die Form:
5 5 5 5 5
@ 5 Die Diskretisierung des Stokes-Problems erfolgt wieder dadurch, daß wir Geschwindigkeit
und Druck in unterschiedlichen endlichen Räumen B und G suchen wollen.
15.7. KONVERGENZ
187
15.6.2 Die Babuska-Brezzi-Bedingung
In Analogie zur Abschätzung (15.15) kann man folgendes Stabilitätskriterium aufstellen:
* Man kann weiterhin zeigen, daß dieses Stabilitätskriterium dann erfüllt ist, wenn die sogenannte Babuska-Brezzi-Bedingung
@div ( 5
@ $ (15.21)
erfüllt ist. Bei dieser Abschätzung betrachten wir zuerst den wichtigen hinteren Teil, er besagt,
daß der vordere Teil größer als Null sein muß, d.h. er darf nicht Null werden. Daher soll
uns vom vorderen Teil der Abschätzung nur der Zähler interessieren. Das darin befindliche
Integral soll also im Supremum echt größer als Null sein. Dies ist eigentlich sehr einfach, man
braucht also nur den Raum B genügend groß wählen, dann wird sich für jedes @ G schon
irgendein Kandidat B finden, der diese Bedingung erfüllt. Ziel ist es aber, den Raum B
in Abstimmung mit dem Raum G möglichst klein zu wählen.
Normalerweiser erreicht man dies, indem man für die Geschwindigkeitsapproximation Ansatzfunktionen höherer Ordnung als für den Druck wählt. So führt die Wahl quadratischer
Ansatzfunktionen für die Geschwindigkeit und linearer Ansatzfunktionen für den Druck zu
stabilen Verfahren. Werden Ansatzfunktionen gleicher Ordnung gewählt, so ist der Interpolationsraum für den Druck zu reichhaltig, die Lösung fängt an zu oszillieren.
15.7 Konvergenz
Die Konvergenztheorie von FE-Verfahren benötigt Methoden der Funktionalanalysis und Approximationstheorie und daher sollen auch hier nur die wesentlichen Ergebnisse vorgestellt
werden.
Die Konsistenzordnung in der Zeit ist bei der Darstellung der Zeitableitung durch Eulerverfahren 1 und bem Crank-Nicolson-Verfahren 2. Der räumliche Diskretisierungsfehler ist bei
elliptischen (z.B. stationäre Probleme) und parabolischen Problemen durch:
fehler (15.22)
wobei ein Maß für die Kantenlänge ist und die Dreiecke nicht zu spitz sind. Eine sehr gute
Einführung in die Konvergenztheorie der FEM ist [16].
188
KAPITEL 15. FINITE ELEMENTE UND FUNKTIONALANALYSIS
Kapitel 16
Gleichungslöser
Die bei impliziten Verfahren entstehenden Gleichungssysteme sind
linear, da die nichtlinearen advektiven Terme entweder explizit oder mit Hilfe des
Newton-Raphson-Verfahrens linearisiert werden und
schwach besetzt, da lediglich die bei der Diskretisierung der partiellen Ableitungen verwendeten nächsten Knoten nichtverschwindende Matrixelemente erzeugen.
Bei den linearen Gleichungslösern unterscheidet man direkte und iterative Verfahren. Direkte
Verfahren lösen das Gleichungssystem durch eine sukzessive Elimination der gesuchten Variablen, iterative Gleichungslöser bestimmen die Lösung als Grenzwert einer Näherungsfolge.
16.1 Diskretisierung und Systemmatrix
Wesentlich zur Effizienz der Gleichungslöser trägt bei, wie die nichtlinearen advektiven Terme diskretisiert werden. Da bei diesen fast immer upwind-Strategien zu berücksichtigen sind,
entstehen hier unsymmetrische Terme. Dies ist der Grund dafür, daß hier oft explizite Verfahren mit den entsprechenden Stabilitätsrestriktionen verwendet werden, damit die advektiven
Terme nicht in das Gleichungssystem mit einbezogen werden müssen. Die Diskretisierung der
diffusiven Terme erzeugt dagegen immer symmetrische und positiv-definite Matrizen.
Die Numerierung der Knoten entscheidet dann über das spezielle ’Aussehen der Systemmatrizen. Bei zweidimensionalen FD-Verfahren erzeugt eine Numerierung in Richtung der Koordinatenachsen eine Tridiagonalmatrix (Bandweite 2), während eine Diskretisierung in Diagonalen zu den Koordinatenachsen eine Blocktridiagonalmatrix erzeugt.
Auch bei der Anwendung von FE-Methoden versucht man, den Rechenaufwand durch intelligente Numerierung zu minimieren. Hier sind vor allem die Numerierungsalgorithmen von
Rosen und Cuthill-McKee zu nennen, die Matrizen minimaler Bandbreite zu erzeugen versuchen.
16.2 Direkte Gleichungsl öser
Wegen der Größe der zu lösenden Gleichungssysteme kommen direkte Verfahren in mehrdimensionalen hydrodynamisch-numerischen Modellen nur selten zur Anwendung.
189
KAPITEL 16. GLEICHUNGSLÖSER
190
Gleichungslöser
Voraussetzungen
Anzahl der Operationen
Gauß/LU-Zerlegung
Pivotisierbarkeit
- für Bandmatrix
Pivotisierbarkeit
Cholesky-Zerlegung symmetrisch, positiv-definit
- für Bandmatrix symmetrisch, positiv-definit
Tabelle 16.1: Direkte Gleichungslöser
Neben der Frage, welche Voraussetzungen die Matrix zur Anwendung des Gleichungslösers
erfüllen muß, ist man speziell an seiner Lösungsgeschwindigkeit interessiert. Diese wird über
die Anzahl der erforderlichen wesentlichen Operationen d.h. Multiplikationen und Divisionen
quantifiziert. In Tafel 16.1 sind dafür die Anwendungsvoraussetzungen und der führende Term
der wesentlichen Operationen in Abhängigkeit von der Größe des Gleichungssystems für das
Gaußsche Eliminationsverfahren und die Cholesky-Zerlegung dargestellt.
16.3 Iterative Gleichungsl öser
Iterative Gleichungslöser führen die Lösung eines Gleichungssystems auf eine Fixpunktform
+ (16.1)
mit einer Iterationsmatrix + zurück. Das Verfahren konvergiert dann, wenn der Spektralradius
der Iterationsmatrix d.h. der Betrag des größten Eigenvektors kleiner eins ist. Die Konvergenzgeschwindigkeit ist dann umso schneller, desto kleiner der Spektralradius ist.
Die Effektivität eines iterativen Gleichungslösers kann deshalb nicht so einfach bestimmt werden, wie es für direkte Gleichungslöser der Fall war, denn sie setzt sich aus der Konvergenzgeschwindigkeit und dem Rechenaufwand pro Iteration zusammen. So kann ein sehr rechenaufwendiges Verfahren dennoch sehr schnell sein, wenn die Konvergenzgeschwindigkeit entsprechend hoch ist.
Als iterative Gleichungslöser werden vielfach die der Klassen der SOR, ADI und CGVerfahren verwendet, wobei Gesamtschritt- und Einzelschrittverfahren wegen ihrer langsamen
Konvergenz nur selten Anwendung finden.
SOR-Verfahren (Successive Overrelaxation) sind Verallgemeinerungen des Einzelschrittverfahrens, bei denen das Ergebnis des letzteren über einen Relaxationsparameter zwischen den
Iterationen gewichtet wird. Die Konvergenz des SOR-Verfahrens ist dann gesichert, wenn die
16.3. ITERATIVE GLEICHUNGSLÖSER
191
Systemmatrix positiv definit ist. Die Konvergenzgeschwindigkeit hängt von der Wahl des Relaxationsparameters ab, für den im Fall einer sogenannten Diagonal-Block-Tridiaginalmatrix
ein optimaler Wert bestimmt werden kann.
ADI-Verfahren (Alternating-Direction-Implicit) lösen das Gleichungssystem, indem die Systemmatrix in ihre verschiedenen Anteile der Ortsableitungen zerlegt wird. Bei einem zweidimensionalen Problem wird also abwechselnd das Gleichungssystem, welches aus der Diskretisierung der x-Ableitungen und den y-Ableitungen entsteht, gelöst. Dabei wird die Reihenfolge in jedem Iterationsschritt gewechselt. ADI-Verfahren sind mit Operator-Splitting Methoden
verwandt, diese verallgemeinern die Grundidee des ADI-Verfahrens auf die zu lösenden Differentialgleichungen und zerlegen diese etwa auch nach physikalischen Anteilen. ADI-Verfahren
eigenen sich außerdem besonders für die Verarbeitung auf Vektor- und Parallelrechnern.
Konjugierte-Gradienten-Verfahren (CG-Conjugate Gradients) führen die Lösung eines Gleichungssystems auf die Minimierung einer quadratischen Form zurück. Sie benötigen pro Iterationsschritt einen relativ hohen Rechenaufwand, dafür ist ihre Konvergenzgeschwindigkeit
allerdings sehr hoch.
Schließlich seien noch Mehrgitterverfahren erwähnt, die ein Gleichungssystem auf verschiedenen sukzessive grober werdenden Gittern löst. Hierdurch kann ein optimales Konvergenzverhalten erzielt werden, da hochfrequente Fehleranteile auf den feinen und niederfrequente
auf den groben Gittern eliminiert werden. Allerdings haben Mehrgitterverfahren noch keinen
Einzug in kommerzielle hydrodynamisch-numerische Programmpakete genommen.
192
KAPITEL 16. GLEICHUNGSLÖSER
Literaturverzeichnis
[1] Benton, E.R. and Platzman, G.W. A Table of Solutions of the One-Dimensional Burgers
Equation. Quaterly of Applied Mathematics, pages 195–221, July 1972. 10.2
[2] Bermejo, R. and Staniforth, A. The Conversion of Semi-Lagrangian Advection Schemes
to Quasi-Monotone Schemes. Monthly Weather Rev., 120:2622–2632, 1992. 8.3
[3] Brooks, A.N. and Hughes, T.J.R. Streamline Upwind/Petrov-Galerkin Formulations for
Convection Dominated Flows with Particular Emphasis on the Incompressible NavierStokes Equations. Comp. Meth. Appl. Mech. Engng., 32:199–259, 1982. 14.4.2, 14.4.2
[4] Burgers, J.M. A Mathematical Model Illustrating the Theory of Turbulence. In von
Mises, R. and von Karman, T., editors, Advances in Applied Mechanics, pages 171 –
199, New York, 1948. Academic Press. 10
[5] Casulli, V. Semi-implicit finite difference methods for the two-dimensional shallowwater equations. J. Comp. Phys., 86:56–74, 1990. 12.4
[6] Casulli, V. and Stelling, G.S. Simulation of Three-Dimensional, Non-Hydrostatic FreeSurface Flows for Estuaries and Coastal Seas. In Spaulding, M.L. and Cheng, R.T., editors, Estuarine and Coastal Modeling, Proceedings of the 4th International Conference,
1996. 13.2.2
[7] Chen, P.P. Der Entity-Relationship-Ansatz zum logischen Systementwurf. Spektrum Akademischer Verlag, 1991. 2.7.1
[8] Christie, I., Griffiths, D.F., Mitchell, A.R., and Zienkiewicz, O.C. Finite Element Methods for Second Order Differential Equations with Significant First Derivatives. Int. J.
num. Meth. Engng., 10:1389–1396, 1976. 14.4.1
[9] Gärtner, S. Zur diskreten Approximation kontinuumsmechanischer Bilanzgleichungen.
Bericht Nr. 24, Institut für Strömungsmechanik und Elektron. Rechnen im Bauwesen der
Universität Hannover, Universität Hannover, Hannover, 1987. 2.3.4, 9.3.3
[10] Gresho, P.M. and Lee, R.L. Don’t Supress the Wiggles - They’re Telling you Something.
Computers and Fluids, 9:223–253, 1981. 4.3
[11] Hasbani, Y., Livne, E., and Bercovier, M. Finite Elements and Characteristics Applied to
Advection-Diffusion Equations. Computers and Fluids, 11 (2):71–83, 1983. 14.3
193
194
LITERATURVERZEICHNIS
[12] Hervouet, J.M. Application de la méthode des caractéristiques en formulation faible à
la résolution des équations d’advection bidimensionelles sur les maillages grilles. Note
technique E41/84.11, Électricité de France, Direction des Études et Recherches, Laboratoire National d’Hydraulique, Chatou, 1984. 14.3
[13] Hervouet, J.M. Remarques sur la convection faible et son application à l’hydraulique.
Note technique HE41/86.09, Électricité de France, Direction des Études et Recherches,
Laboratoire National d’Hydraulique, Chatou, 1986. 14.3
[14] Hervouet, J.M. and Moulin, C. Nouveaux Schemas de Convection dans Telemac-2D v
2.3: Apport de la Methode SUPG. Note technique HE43/93-27, Électricité de France,
Direction des Études et Recherches, Laboratoire National d’Hydraulique, Chatou, 1993.
14.4.2
[15] Heuser, H. Gewöhnliche Differentialgleichungen. B.G. Teubner, Stuttgart, 1989. 2.8
[16] Johnson, C. Numerical Solution of Partial Differential Equations by the Finite Element
Method. Cambridge University Press, Cambridge, N.Y., 1987. (Exzellente Darstellung
der (Fehler-) Theorie der FE-Methode). 15.7
[17] Lax, P. and Wendroff, B. Systems of Conservation Laws. Comm. Pure Appl. Math.,
13:217–237, 1960. 2.3.4
[18] MacCormack, R.W. Numerical Solution of a Shock Wave with a Laminar Boundary.
Lecture Notes in Physics, 8:151–163, 1971. 5.2.3
[19] Matsoukis, P.-F.C. Tidal Models Using Method of Characteristics. J. Waterway, Port,
Coastal and Ocean Eng., 118 (3):233–248, 1992. 12.3.2
[20] McCalpin, J.D. A Quantitative Analysis of the Dissipation Inherent in Semi-Lagrangian
Advection. Monthly Weather Rev., 116:2330–2336, 1988. 8.2
[21] Patankar, S.V. Numerical Heat Transfer and Fluid Flow. Hemisphere Publishing, New
York, Philadelphia, London, 1980. (Der Schöpfer des SIMPLE-Algorithmus). 13.2.2
[22] Ramakrishnan, C.V. An Upwind Finite Element Scheme for the Unsteady Convective
Diffusive Transport Equation. Appl. Math. Modelling, 3:280–284, 1979. 9.3.3
[23] Roache, P.J. Computational Fluid Dynamics. Hermosa Publishers, Albuquerque, New
Mexico, 1972. (Etwas veraltet, aber immer noch ein Klassiker). 4.3
[24] Spalding, D.B. A Novel Finite Difference Formulation for Differential Expressions Involving Both First and Second Derivatives. Int.J.Num.Meth.Engng., 4:551–559, 1972.
4.3
[25] Thomasset, F. Implementation of Finite Element Methods for Navier-Stokes Equations.
Springer-Verlag, New York, Heidelberg, Berlin, 1981. (Eine didaktische Einführung,
daher kein Überblick). 4.3
LITERATURVERZEICHNIS
195
[26] Urban, C. Ein Finite-Element-Verfahren mit linearen Ansätzen für stationäre zweidimensionale Strömungen. Bericht Nr. 22, Institut für Strömungsmechanik und Elektron.
Rechnen im Bauwesen der Universität Hannover, Hannover, 1986. 12.3.1
[27] Weiyan, T. Shallow Water Hydrodynamics. Number 55 in Elsevier Oceanography Series.
Elsevier, Amsterdam, 1992. (FD-lastig, jedoch ausgezeichneter Überblick zur Theorie
und Lösung der tiefengemittelten Flachwassergleichungen). 12.5
Документ
Категория
Без категории
Просмотров
41
Размер файла
1 706 Кб
Теги
der, methodes, numerische, 001, pdf, stroemungsmechanik, 4156
1/--страниц
Пожаловаться на содержимое документа