CFD-Berechnungen mit uberlappenden Gittern¨ in Bezug auf ...

Universitat Ulm

Fakultat fur Mathematik und Wirtschaftswissenschaften

CFD-Berechnungen

mit uberlappenden Gittern

in Bezug auf den

Voith-Schneider-Propeller

Diplomarbeit

zur Erlangung des akademischen Grades Diplom-Wirtschaftsmathematiker

an der Fakultat fur Mathematik und Wirtschaftswissenschaften der

Universitat Ulm

vorgelegt von

Andreas Helbrich

am 12. November 2008

Gutachter

Prof. Dr. Karsten Urban

Prof. Dr. Stefan Funken

Inhaltsverzeichnis

Inhaltsverzeichnis i

Abbildungsverzeichnis vi

Tabellenverzeichnis x

Abkurzungsverzeichnis xi

1 Einleitung 1

2 Der Voith-Schneider-Propeller 4

2.1 Funktionsweise . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4

2.2 Bedeutung . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6

2.3 Problemstellung . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6

3 Numerische Stromungsmechanik 8

3.1 Modellgleichungen . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9

3.1.1 Massenerhaltung . . . . . . . . . . . . . . . . . . . . . . . . . 10

3.1.2 Impulserhaltung . . . . . . . . . . . . . . . . . . . . . . . . . . 11

3.1.3 Universelle Transportgleichung . . . . . . . . . . . . . . . . . . 12

3.1.4 Randbedingungen und Anfangswerte . . . . . . . . . . . . . . 12

3.2 Finite Differenzen . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13

3.2.1 Begriffe der Konvergenztheorie . . . . . . . . . . . . . . . . . . 15

i

INHALTSVERZEICHNIS ii

3.3 Finite Volumen . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16

3.4 Finite Elemente . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17

3.4.1 Variationsformulierung . . . . . . . . . . . . . . . . . . . . . . 18

3.4.2 Ritz-Galerkin-Verfahren . . . . . . . . . . . . . . . . . . . . . 20

3.4.3 Zusammengesetzte, stuckweise lineare Basisfunktionen . . . . 21

3.5 FV-Diskretisierung . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22

3.5.1 Approximation von Integralen . . . . . . . . . . . . . . . . . . 23

3.5.2 Konvektion . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24

3.5.3 Diffusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27

3.5.4 Gradientenberechung im Zellmittelpunkt . . . . . . . . . . . . 27

3.5.5 Quellterme . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28

3.5.6 Behandlung der Druckgleichung . . . . . . . . . . . . . . . . . 28

3.5.7 Zeitintegration . . . . . . . . . . . . . . . . . . . . . . . . . . 29

3.5.8 Resultierende algebraische Gleichungen . . . . . . . . . . . . . 30

3.5.9 Druckberechnung mit dem SIMPLE-Algorithmus . . . . . . . 31

3.5.10 Randbedingungen . . . . . . . . . . . . . . . . . . . . . . . . . 33

3.6 Turbulenzmodellierung . . . . . . . . . . . . . . . . . . . . . . . . . . 33

3.6.1 RANSE . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34

3.6.2 k- ε -Turbulenzmodell . . . . . . . . . . . . . . . . . . . . . . . 35

3.7 Losen des Gleichungssystems . . . . . . . . . . . . . . . . . . . . . . . 36

3.7.1 Sequentielles Losungsverfahren . . . . . . . . . . . . . . . . . 37

3.7.2 Unterrelaxation . . . . . . . . . . . . . . . . . . . . . . . . . . 38

3.7.3 LGS-Loser . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38

3.7.4 Losungsalgorithmus . . . . . . . . . . . . . . . . . . . . . . . . 40

INHALTSVERZEICHNIS iii

4 Gebietszerlegungsverfahren 41

4.1 Zerlegung in Teilprobleme . . . . . . . . . . . . . . . . . . . . . . . . 42

4.2 Uberlappende Verfahren . . . . . . . . . . . . . . . . . . . . . . . . . 45

4.2.1 Das alternierende Schwarz-Verfahren . . . . . . . . . . . . . . 45

4.2.2 Additives Schwarz-Verfahren . . . . . . . . . . . . . . . . . . . 48

4.3 Gebietszerlegung ohne Uberlappung . . . . . . . . . . . . . . . . . . . 48

5 Das Overlapping Grid Verfahren fur Finite Volumen mit bewegten

Gittern 51

5.1 Einleitung . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51

5.2 Overlapping Grid Methodik . . . . . . . . . . . . . . . . . . . . . . . 53

5.2.1 Zellfunktionen . . . . . . . . . . . . . . . . . . . . . . . . . . . 54

5.2.2 Bestimmung inaktiver Zellen . . . . . . . . . . . . . . . . . . . 56

5.2.3 Spenderzellensuche . . . . . . . . . . . . . . . . . . . . . . . . 57

5.2.4 Interpolationsschema . . . . . . . . . . . . . . . . . . . . . . . 58

5.3 Losungsverfahren . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60

5.3.1 Behandlung inaktiver Zellen . . . . . . . . . . . . . . . . . . . 60

5.3.2 Starke Koppelung der Teillosungen . . . . . . . . . . . . . . . 60

5.3.3 Korrektur der Massenerhaltung . . . . . . . . . . . . . . . . . 62

5.3.4 Zellwerte an Overlap-Grenzflachen . . . . . . . . . . . . . . . 63

5.4 Modellierung bewegter Gitter . . . . . . . . . . . . . . . . . . . . . . 64

5.4.1 ALE-Systemgleichungen . . . . . . . . . . . . . . . . . . . . . 64

5.4.2 Geometrieerhaltungsgesetz . . . . . . . . . . . . . . . . . . . . 65

5.4.3 Anderungen in der Diskretisierung . . . . . . . . . . . . . . . 66

5.5 Aktuelle OLG-Software . . . . . . . . . . . . . . . . . . . . . . . . . . 68

5.5.1 Tau (DLR) . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68

5.5.2 Chimera Grid Tools (NASA, Air Force) . . . . . . . . . . . . . 68

INHALTSVERZEICHNIS iv

5.5.3 Overture (LLNL) . . . . . . . . . . . . . . . . . . . . . . . . . 68

5.5.4 FEFLO (NRL) . . . . . . . . . . . . . . . . . . . . . . . . . . 69

5.5.5 Comet (CD-adapco) . . . . . . . . . . . . . . . . . . . . . . . 69

6 Anleitung fur Overlap-Rechnungen mit Comet 70

6.1 OLG-Programmbefehle . . . . . . . . . . . . . . . . . . . . . . . . . . 70

6.1.1 Praprozessor-Befehle . . . . . . . . . . . . . . . . . . . . . . . 71

6.1.2 Postprozessor-Befehle . . . . . . . . . . . . . . . . . . . . . . . 74

6.1.3 Usercoding-Befehle . . . . . . . . . . . . . . . . . . . . . . . . 75

6.2 Richtlinien fur die Gitterstruktur im Overlap-Bereich . . . . . . . . . 75

6.3 Parallelisierung . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78

6.4 Simulation des VSP mit Overlapping Moving Grid . . . . . . . . . . . 79

6.4.1 Anderungen der Vernetzung . . . . . . . . . . . . . . . . . . . 79

6.4.2 Anderungen im Usercoding . . . . . . . . . . . . . . . . . . . . 80

6.4.3 Simulationsrechnungen . . . . . . . . . . . . . . . . . . . . . . 80

6.5 Modellierung von Anbauteilen . . . . . . . . . . . . . . . . . . . . . . 80

6.5.1 Streben als Aussparung im Hintergrundgitter . . . . . . . . . 81

6.5.2 Streben als Overlap-Korper . . . . . . . . . . . . . . . . . . . 81

6.5.3 Schutzplatte als Aussparung im Hintergrundgitter . . . . . . . 82

6.5.4 Schutzplatte als Overlap-Korper . . . . . . . . . . . . . . . . . 83

6.5.5 Overlap-Aussparungs-Kombinationen . . . . . . . . . . . . . . 83

6.5.6 Overlap-Overlap-Kombinationen . . . . . . . . . . . . . . . . . 84

6.5.7 Aussparungs-Kombinationen . . . . . . . . . . . . . . . . . . . 84

6.6 Problembehandlung . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85

6.6.1 Bekannte Probleme und Losungen . . . . . . . . . . . . . . . . 86

INHALTSVERZEICHNIS v

7 Ergebnisse der Simulationsrechnungen 88

7.1 Modellvalidierung . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 88

7.1.1 Flugelumstromung mit ursprunglicher Vernetzung . . . . . . . 91

7.1.2 Flugelumstromung mit modifizierter Vernetzung . . . . . . . . 92

7.1.3 VSP Einflugler mit ursprunglicher Vernetzung . . . . . . . . . 93

7.1.4 VSP Einflugler mit modifizierter Vernetzung . . . . . . . . . . 94

7.1.5 VSP Einflugler mit Vernetzung fur SP-Rechnungen . . . . . . 94

7.1.6 VSP Funfflugler . . . . . . . . . . . . . . . . . . . . . . . . . . 95

7.1.7 Laufzeiten VSP-Rechnungen (Vernetzung fur SP) . . . . . . . 96

7.2 Rechnungen mit Anbauteilen . . . . . . . . . . . . . . . . . . . . . . . 96

7.2.1 VSP mit Streben . . . . . . . . . . . . . . . . . . . . . . . . . 97

7.2.2 VSP mit Schutzplatte . . . . . . . . . . . . . . . . . . . . . . 97

7.2.3 VSP mit Streben und Schutzplatte . . . . . . . . . . . . . . . 97

8 Fazit und Ausblick 99

A Abbildungen zur Vernetzung 103

B Graphen der Krafte fur VSP ohne Anbauteile 118

C Graphen der Krafte fur VSP mit Anbauteilen 137

Literaturverzeichnis 146

Abbildungsverzeichnis

2.1 Bewegung eines Flugels (links), Flugelstellung fur Schuberzeugung

(rechts), aus Bartels und Jurgens (2006) . . . . . . . . . . . . . . . . 5

2.2 Schutzplatte mit Streben am Rumpf, aus Bartels und Jurgens (2006) 7

3.1 Stetige, stuckweise lineare Shape-Funktion, aus Bollhofer und Mehr-

mann (2004) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21

3.2 Diskretes Kontrollvolumen, aus Hadzic (2005) . . . . . . . . . . . . . 23

3.3 Lineare Interpolation in der Mitte der KV-Seite, aus Hadzic (2005) . 25

4.1 Diskretisierter Losungsraum bestehend aus zwei Teilgebieten, aus Ma-

thew (2008) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41

4.2 Losungsraum bestehend aus zwei uberlappenden Teilgebieten, aus Ma-

thew (2008) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42

5.1 OLG-System aus Hintergrundgitter und zwei Overlap-Objekten, aus

Hadzic (2005) mit Ubersetzungen des Autors . . . . . . . . . . . . . . 53

5.2 Uberlappungsgebiet mit Bezeichnungen, aus Hadzic (2005) mit Uber-

setzungen des Autors . . . . . . . . . . . . . . . . . . . . . . . . . . . 54

5.3 Aktive Teilbereiche bei Mehrfachuberlappungen, aus Hadzic (2005)

mit Ubersetzungen des Autors . . . . . . . . . . . . . . . . . . . . . . 56

5.4 Neighbor-to-neighbor Suchverfahren, aus Hadzic (2005) mit Uberset-

zungen des Autors . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57

5.5 Massenfluss an der Overlap-Grenzflache, aus Hadzic (2005) mit Uber-

setzungen des Autors . . . . . . . . . . . . . . . . . . . . . . . . . . . 62

vi

ABBILDUNGSVERZEICHNIS vii

5.6 Gradienten-Approximation an Overlap-Grenzflache, aus Hadzic (2005)

mit Ubersetzungen des Autors . . . . . . . . . . . . . . . . . . . . . . 64

5.7 Bewegung eines KV uber drei Zeitschritte mit jeweils uberstrichenem

Volumen, aus Hadzic (2005) . . . . . . . . . . . . . . . . . . . . . . . 66

A.1 Flugel mit ersten Vernetzungsschichten (links), Druckverteilung auf

Flugelblatt (rechts) . . . . . . . . . . . . . . . . . . . . . . . . . . . . 103

A.2 Ursprunglicher Flugelzylinder (links), 166800 Zellen, reduzierter Flugel-

zylinder (rechts), 130000 Zellen . . . . . . . . . . . . . . . . . . . . . 104

A.3 Reduziertes Flugelgitter, Nahaufnahme . . . . . . . . . . . . . . . . . 104

A.4 Schnitt Hintergrundgitter Sliding/Econ (oben), 423400 Zellen, Hinter-

grundgitter Overlap (unten), 174528 Zellen, mit jeweiligem Flugelgitter105

A.5 Overlap-Zellfunktionen (IOVC, vgl. Tabelle 6.1), Schnitt Hintergrund-

gitter (oben), vollstandiges Netz (unten) . . . . . . . . . . . . . . . . 106

A.6 VSP Einflugler: Sliding-Fall (oben), 309882 Zellen Hintergrundgitter,

Overlap-Fall (unten), 435520 Hintergrundzellen mit jeweiligem Flugel 107

A.7 Fur Overlap modifizierter Flugel (links), 102800 Zellen, und entspre-

chender Zylinder fur Sliding-Referenzfalle (rechts), 136100 Zellen . . . 108

A.8 Modifiziertes Flugelgitter, Nahaufnahme Oberseite . . . . . . . . . . . 108

A.9 Unterseite ursprungliches Flugelgitter (oben) und modifiziertes Flugel-

gitter (unten) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 109

A.10 Gestauchtes Flugelgitter fur Rechnungen mit Schutzplatte, 74200 Zel-

len, 8 Zellschichten zur Flugelspitze . . . . . . . . . . . . . . . . . . . 110

A.11 Funfflugler Sliding (765630 Zellen, davon 230340 Hintergrund) (oben),

und Overlap (1210772 Zellen, davon 839572 Hintergrund) (unten) . . 111

A.12 Parallelisierung: Erfolgreiche Partitionierung auf 2 Knoten (links),

Probleme bei 3 Knoten (rechts) . . . . . . . . . . . . . . . . . . . . . 112

A.13 Vereinfachte Streben: Overlapkorper (o.l.), IOVC Hintergrund (o.r.),

Strebe als Aussparung (u.l.), Kontakt Overlap-Rand Aussparung fi-

nale Strebe (u.r.) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 112

ABBILDUNGSVERZEICHNIS viii

A.14 Vereinfachte Schutzplatte: Druckverteilung (o.l.), IOVC Hintergrund

auf Hohe Flugelspitze (o.r.), Gitter fur Overlap-Modellierung (u.l.),

und Aussparung (u.r.) . . . . . . . . . . . . . . . . . . . . . . . . . . 113

A.15 Vereinfachte Schutzplatte mit Strebe: Gemeinsame Aussparung (o.l.),

aktive Zellen fur Strebe und SP als eigene Overlap-Korper (o.r.), ak-

tive Zellen fur Strebe als Overlap, SP als Aussparung (u.) . . . . . . . 114

A.16 Finale Schutzplatte: IOVC-Schnitt Hintergrund mit SP-Netz in Block-

form (o.), IOVC-Verteilung fur verfeinertes SP-Gitter auf Hohe Flugel-

spitze mit Ungleichmaßigkeiten (m. und u.) . . . . . . . . . . . . . . . 115

A.17 VSP Einflugler mit finaler SP und 3 Streben: Oberflachennetz (o.)

und Druckverteilung (m. und u.) . . . . . . . . . . . . . . . . . . . . 116

A.18 VSP Funfflugler mit finaler SP und 3 Streben: Druckverteilung (o.)

mit Unregelmaßigkeit (u.) . . . . . . . . . . . . . . . . . . . . . . . . 117

B.1 Stationare Flugelumstromung mit 12 m/s, Anstellwinkel 4 1/2 . . . 119

B.2 Stationare Flugelumstromung mit 12 m/s, Anstellwinkel 4 2/2 . . . 120

B.3 Rel. Abweichung Stat. Umstromung mit 12 m/s, Anstellwinkel 4 . . 121

B.4 Umstromung mod. Fl.Gitter mit 6,477 m/s, Anstellwinkel 4 1/2 . . 122

B.5 Umstromung mod. Fl.Gitter mit 6,477 m/s, Anstellwinkel 4 2/2 . . 123

B.6 Rel. Abweichung Umstromung mod. Gitter mit 6,477 m/s, Winkel 4 124

B.7 VSP (ein ursprungliches Flugelgitter) 1/2 . . . . . . . . . . . . . . . . 125

B.8 VSP (ein ursprungliches Flugelgitter) 2/2 . . . . . . . . . . . . . . . . 126

B.9 VSP (ein modifiziertes Flugelgitter) 1/2 . . . . . . . . . . . . . . . . 127

B.10 VSP (ein modifiziertes Flugelgitter) 2/2 . . . . . . . . . . . . . . . . 128

B.11 VSP (ein Flugelgitter fur SP-Rechnungen) 1/2 . . . . . . . . . . . . . 129

B.12 VSP (ein Flugelgitter fur SP-Rechnungen) 2/2 . . . . . . . . . . . . . 130

B.13 Relative Abweichung Krafte VSP Einflugler . . . . . . . . . . . . . . 131

B.14 VSP Funfflugler, Krafte auf 1. Flugel 1/2 . . . . . . . . . . . . . . . . 132

B.15 VSP Funfflugler, Krafte auf 1. Flugel 2/2 . . . . . . . . . . . . . . . . 133

ABBILDUNGSVERZEICHNIS ix

B.16 Relative Abweichung Funfflugler, Krafte 1. Flugel . . . . . . . . . . . 134

B.17 VSP Funfflugler, summierte Krafte . . . . . . . . . . . . . . . . . . . 135

B.18 Relative Abweichung Funfflugler, summierte Krafte . . . . . . . . . . 136

C.1 VSP Einflugler mit Schutzplatte . . . . . . . . . . . . . . . . . . . . . 138

C.2 Krafte auf vereinfachte Schutzplatte, VSP Einflugler . . . . . . . . . . 139

C.3 VSP Ein- und Funfflugler mit SP und 3 Streben, Krafte 1. Flugel . . 140

C.4 VSP Funfflugler mit SP und 3 Streben, summierte Krafte . . . . . . . 141

C.5 Krafte auf Streben bei VSP Einflugler . . . . . . . . . . . . . . . . . 142

C.6 Krafte auf Schutzplatte bei VSP Einflugler . . . . . . . . . . . . . . . 143

C.7 Krafte auf Streben bei VSP Funfflugler . . . . . . . . . . . . . . . . . 144

C.8 Krafte auf Schutzplatte bei VSP Funfflugler (eingeschrankter Werte-

bereich) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 145

Tabellenverzeichnis

3.1 Zeitintegrationskoeffizienten in Gleichung (3.64) . . . . . . . . . . . . 31

6.1 Farbtabelle der Zellfunktionen . . . . . . . . . . . . . . . . . . . . . . 74

6.2 Fehlermeldungen OLG . . . . . . . . . . . . . . . . . . . . . . . . . . 85

7.1 Mittlere Krafte einer stationare Flugelumstromung mit 12 m/s, An-

stellwinkel 4 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91

7.2 Mittlere Krafte einer stationaren Flugelumstromung mit 6,477 m/s,

Anstellwinkel 4 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92

7.3 Mittlere Krafte VSP (ein Flugel, ursprungliche Vernetzung) . . . . . 93

7.4 Mittlere Krafte VSP (ein Flugel, modifizierte Vernetzung) . . . . . . 94

7.5 Mittlere Krafte VSP (ein Flugel fur SP-Rechnungen) . . . . . . . . . 94

7.6 Mittlere Krafte VSP (funf Flugel fur SP-Rechnungen) . . . . . . . . . 95

7.7 Laufzeiten VSP fur SP-Rechnungen, 1500 Zeitschritte . . . . . . . . . 96

x

Abkurzungsverzeichnis

CFD Numerische Stromungsmechanik,

engl. Computational Fluid Dynamics

CG Konjugierte Gradienten,

engl. conjugate gradients

DNS direkte numerische Simulation

FD Finite Differenzen

FEM Finite Elemente Methode

FV Finite Volumen

FWK Flugelwinkelkurve

IC Unvollstandige Cholesky-Zerlegung,

engl. incomplete Cholesky, incomplete LU

KV Kontrollvolumen

MPI Message Passing Interface

OLG Gitter mit Uberlappung,

engl. Overlapping Grid

PDE partielle Differentialgleichung,

engl. partial differential equation

RANSE Reynolds-gemittelte Navier-Stokes Gleichungen,

engl. Reynolds-averaged-Navier-Stokes Equations

SCL Geometrieerhaltungsgesetz

engl. Space Conservation Law

SP Schutzplatte

VSP Voith-Schneider-Propeller

VWT Voith Wassertrecker

xi

Kapitel 1

Einleitung

Die Arbeit eines Ingenieurs ist mittlerweile ohne Computerunterstutzung kaum noch

vorstellbar. Seien es Statikberechnungen von Bauwerken, Crashtest-Simulationen in

der Fahrzeugentwicklung oder Stromungsoptimierungen von Flugzeugen und Schiffen

- immer spielen Simulationsrechnungen an Hochleistungscomputern eine unverzicht-

bare Rolle fur die Losung dieser Probleme.

Die mathematisch-physikalischen Modelle, die die Grundlage dieser Rechnungen bil-

den, werden in oft jahrelanger Zusammenarbeit von Mathematikern, Physikern,

Informatikern und Ingenieuren entwickelt und immer weiter verbessert. Auch mit

standig wachsender Rechenleistung der verfugbaren Computer ist es immer noch

nicht moglich, alle Eigenschaften eines technischen Problems exakt wiedergeben zu

konnen - die Modellbildung mit der heutigen Technik erfordert immer eine genaue

Uberlegung, was moglich ist, aber auch welche Detailtreue notig ist, um mit ver-

tretbaren Kosten und Zeitaufwand zu einer vernunftigen Losung fur das vorliegende

Problem zu kommen. Hierbei ist die interdisziplinare Arbeit von Spezialisten ver-

schiedenster Fachrichtungen entscheidend, um die Grenzen des technisch Machbaren

zu erweitern.

Um letztendlich teure Experimente weitgehend durch Simulationsrechnungen zu er-

setzen, oder sogar Probleme angehen zu konnen, die sich aufgrund ihrer Dimensi-

on oder Beschaffenheit nicht experimentell untersuchen lassen, benotigt man Fach-

leute aller beteiligten Teilgebiete: Die Naturwissenschaftler und Physiker mit ihrer

Modellumgebung, um die Phanomene zu beschreiben, Mathematiker, die effizien-

te analytische oder numerische Losungsverfahren fur die Modellgleichungen finden,

Informatiker, um die speicher- und rechenintensiven Losungsalgorithmen geschickt

im Computer zu implementieren und die Daten zu visualisieren, sowie am Ende die

1

KAPITEL 1. EINLEITUNG 2

Ingenieure, um die Ergebnisse der Rechnungen auszuwerten und zu interpretieren.

Die Ubergange zwischen den Teilbereichen sind oft so fließend, dass sich keine klare

Trennung mehr vornehmen lasst.

Die Arbeit von Mathematikern spielt fur viele dieser Teilbereiche eine wichtige Rol-

le: Neben der mathematischen Formulierung, generellen Aussagen uber Losbarkeit

und Eindeutigkeit eines Problems und der Entwicklung von Verfahren zur Losung

ist die mathematische Fehleranalyse der neuen Algorithmen von großer Wichtig-

keit. Fehlerschranken fur Modellvereinfachungen sowie Aussagen uber Stabilitat und

Robustheit der Losungsverfahren bestimmen letztendlich die Wahl einer geeigneten

Modellierung und deren algorithmische Umsetzung. Daruber hinaus ermoglicht eine

formale Analyse die Bewertung der Qualitat von Simulationsergebnissen.

Ein wichtiger Schritt im Prozess der Modellentwicklung von der ersten akademi-

schen Idee zur erfolgreichen Anwendung im Tagesgeschaft eines Ingenieurburos ist

die Verifikation - die Untersuchung, inwiefern sich die Ergebnisse der neuen Model-

lierung mit denen bereits bestehender, anerkannter Verfahren decken, soweit sich

die Anwendbarkeit uberschneidet. Die vorliegende Arbeit ist eine derartige Analyse,

angesiedelt im Bereich der numerischen Stromungsmechanik (CFD, kurz fur engl.

Computational Fluid Dynamics):

Es handelt sich um die Untersuchung der Eignung des Overlapping Grid Verfahrens

fur CFD-Berechnungen im Umfeld der Weiterentwicklung des Voith-Schneider-Pro-

pellers, eines Schiffsantriebs der Firma Voith Turbo Schneider Propulsion. Sie ist

beschrankt auf eine detaillierte Untersuchung der Implementierung des Verfahrens

im Softwarepaket Comet der Firma CD-adapco.

Overlapping Grid steht hierbei fur ein Gebietszerlegungsverfahren, also eine Technik,

das Losungsgebiet von Differentialgleichungen in Stucke aufzuteilen, die sich stellen-

weise uberlappen. Das bringt einerseits die Moglichkeit, die Losungsberechnung auf

den Teilgebieten zu parallelisieren, andererseits den hier viel entscheidenderen Vor-

teil, die Teile getrennt vernetzen zu konnen. Dies macht eine optimale Anpassung des

Netzes an das Problem viel leichter als bei einem zusammenhangenden Gesamtnetz.

Die Teillosungen werden uber Interpolation zwischen Vorder- und Hintergrundgit-

ter zur Gesamtlosung zusammengefugt. Daruber hinaus ermoglicht dieser Ansatz

komplexe Bewegungen der Gitter relativ zueinander, ein entscheidender Vorteil ge-

genuber herkommlichen Verfahren, bei denen die Modellierung großraumiger Bewe-

gungen sehr schwierig werden kann. Nachteile sind u.a. der zusatzliche Aufwand, der

durch den Abgleich der Teillosungen entsteht. Das Verfahren hat seine Tauglichkeit

in einigen Bereichen bereits uber akademische Beispiele hinaus unter Beweis gestellt.

KAPITEL 1. EINLEITUNG 3

Ob es in der vorliegenden Implementierung ausgereift genug ist, die stromungsme-

chanischen Vorgange beim Betrieb des VSP hinreichend zu modellieren, soll diese

Diplomarbeit zeigen. Sie gliedert sich in folgende Teile:

Kapitel 2 beschreibt kurz die Bedeutung und Funktionsweise des Voith-Schneider-

Propellers, sowie die Problemstellung, die die Einfuhrung und Anwendung des neuen

Verfahrens notwendig macht.

In Kapitel 3 werden die Grundlagen der numerischen Stromungsmechanik eingefuhrt

und ein Uberblick uber die aktuell verwendeten Verfahren gegeben.

Kapitel 4 widmet sich den Gebietszerlegungsverfahren. Es wird der mathematische

Hintergrund der Verfahren mit uberlappenden Gittern eingefuhrt sowie die wichtigs-

ten Vertreter der Verfahren mit nicht uberlappenden Gittern genannt.

Das fur diese Arbeit wichtigste Verfahren wird in Kapitel 5 genauer erlautert: Das

Overlapping Grid Verfahren fur Finite Volumen mit bewegten Gittern, das die Basis

fur die im weiteren Verlauf durchgefuhrten Simulationen bildet.

Im folgenden Kapitel 6 werden die Comet-Programmbefehle fur das OLG-Verfahren

eingefuhrt und Richtlinien fur den Aufbau von Rechnungen mit Overlapping Grid

gegeben. Erkenntnisse aus durchgefuhrten Simulationen werden zusammengefasst

und die aufgetretenen Probleme diskutiert.

Kapitel 7 prasentiert die messbaren Resultate der im Rahmen dieser Arbeit durch-

gefuhrten numerischen Simulationen: Zur Modellvalidierung werden Vergleiche von

Overlapping Grid Rechnungen mit den Ergebnissen bisher verwendeter Verfahren

durchgefuhrt. Es werden Krafte von stationaren Flugelumstromungen, dem VSP mit

einem und funf Flugeln, sowie Rechnungen mit Schutzplatte und Streben betrachtet.

Zum Abschluss fasst Kapitel 8 die Ergebnisse zusammen und gibt einen Ausblick auf

zu erwartende Entwicklungen des Overlapping Grid Verfahrens in der numerischen

Stromungsmechanik.

Kapitel 2

Der Voith-Schneider-Propeller

Der Voith-Schneider-Propeller (VSP), der im Mittelpunkt der Untersuchungen dieser

Arbeit steht, ist ein einzigartiger Schiffsantrieb, der auf eine Erfindung Ernst Schnei-

ders von 1925 zuruckgeht. Die Firma Voith griff damals seine Idee auf, entwickelte

sie zusammen mit dem Erfinder bis zur Serienreife, und hat das Prinzip bis heu-

te stetig verbessert. Seit einigen Jahren spielt die numerische Stromungsmechanik

hierbei eine große Rolle. Erfolgreiche Arbeiten in diesem Bereich, wie die erhebliche

Wirkungsgradsteigerung des VSP durch Jurgens u. a. (2007), zeigen, dass auch eine

Idee, an der Ingenieure schon seit uber 80 Jahren arbeiten, noch nicht ausgereizt sein

muss.

Die folgenden Abschnitte geben einen kurzen Uberblick uber die Funktionsweise des

Antriebs, einen Einblick in die Bedeutung des VSP fur die moderne Schifffahrt, sowie

die Problemstellung, die im Rahmen dieser Diplomarbeit untersucht werden soll.

2.1 Funktionsweise

Anders als bei herkommlichen Schiffschrauben, bei denen sich mehrere schrag ge-

stellte Flugel um eine waagrechte Achse drehen und so Schub erzeugen, rotieren

beim VSP senkrecht im Wasser stehende Flugel in einer Kreisbahn um eine vertikale

Drehachse.

Die Flugel fuhren zusatzlich zu ihrer Rotation auf dem so genannten Flugelkreis eine

nach einem bestimmten Gesetz schwingende Bewegung aus, je nach Position auf der

Umkreisbewegung wird der Flugel entsprechend angestellt um Schub zu erzeugen. So

4

KAPITEL 2. DER VOITH-SCHNEIDER-PROPELLER 5

Abbildung 2.1: Bewegung eines Flugels (links), Flugelstellung fur Schuberzeugung

(rechts), aus Bartels und Jurgens (2006)

sind die Flugel beim Halbkreis des Umlaufs, der in Fahrtrichtung liegt, nach außen

gedreht, auf dem anderen nach innen, vergleiche Abb. 2.1.

Uber die Große dieser Winkelausschlage lasst sich der Betrag des Schubs kontrollie-

ren, anders als bei herkommlichen Schraubenpropellerantrieben, bei denen man die

Drehzahl variieren muss, um unterschiedlichen Schub zu erzeugen. Das Gesetz, dem

die Schwenkbewegung folgt, heißt Flugelwinkelkurve (FWK). Ihre Optimierung war

auch der Ansatz der Arbeit von Jurgens u. a. (2007). Eine ausfuhrliche Beschreibung

des Funktionsprinzips nebst Entwicklungsgeschichte findet sich in Jurgens und Fork

(2002).

Die weiteren Moglichkeiten, die sich aus dem Prinzip der Steuerung der Flugelwinkel

ergeben, offenbaren sich allerdings erst auf den zweiten Blick: Es lasst sich nicht nur

Schub nach ”vorne” erzeugen, sondern stufenlos in jede beliebige Richtung. Hervorra-

gende Manovrierfahigkeit wird unter Verwendung eines zweiten VSP in einem Schiff

erreicht, man kann traversieren, d.h. seitwarts fahren, oder auf der Stelle wenden,

und kann im Gegensatz zu rudergelenkten Schiffen auch bei sehr niedrigen Geschwin-

digkeiten exakt manovrieren. Einen Uberblick der Anwendung in der Praxis findet

sich im nachsten Abschnitt:


2.2 Bedeutung

Der VSP wird vor allem dort eingesetzt, wo Schiffe besondere Anforderungen an

Sicherheit und Manovrierfahigkeit erfullen mussen. Heute sind vor allem Voith Was-

sertrecker, Doppelendfahren wie die Staten Island Ferries in New York City, und auch

Minenraumboote, Plattformversorgungsschiffe, Feuerloschboote und Schwimmkrane

mit VSPs verschiedener Große und Bauart ausgestattet.

Besonders die Bedeutung des Voith Wassertreckers (VWT), ein Schlepper mit VSP-

Antrieb, dessen Entwicklung 1952 begann, ist nicht deutlich genug hervorzuheben:

Als Eskortschiff fur Schiffe mit gefahrlicher Ladung wie zum Beispiel Oltanker er-

reicht der VWT fur einen sehr großen Geschwindigkeitsbereich hochste Assistenz-

krafte und hat so die Sicherheit von Schleppoperationen in vielen Teilen der Welt

drastisch erhoht. Vor allem durch die neuen Steuermoglichkeiten konnen einige bis-

her sehr gefahrliche Situationen im Schleppbetrieb vermieden werden, und so sorgt

der VWT dafur, dass große Frachtschiffe, die bei niedrigen Geschwindigkeiten alleine

nur unzureichend manovrieren konnen, mit ihrer Ladung sicher im Hafen ankommen.

Detailliertere Informationen uber die Anwendung des VSP finden sich in Bartels und

Jurgens (2006) sowie Jurgens und Fork (2002).

2.3 Problemstellung

Obwohl der VSP aufgrund seiner Bauart ein sehr robuster und langlebiger Antrieb

ist, muss er in der Praxis durch eine Schutzplatte und Streben vor Beschadigungen

durch Baumstamme o.a. geschutzt werden. Diese Platte, die durch stromungsmecha-

nische Wechselwirkungen auch zur Schuberzeugung beitragt, wird mit den Streben

unmittelbar unter den Flugelenden am Schiffsrumpf befestigt. Sie befindet sich also

in direkter Nahe zu den bewegten Teilen, vgl. Abb. 2.2.

Fur CFD-Rechnungen, wie die Firma Voith sie im Rahmen der Weiterentwicklung

des VSP durchfuhrt, stellt genau das ein Problem dar: In den bisherigen Stromungs-

mechanikrechnungen erweist sich besonders die Modellierung der Streben als nicht

durchfuhrbar: Die Realisierung der Flugelbewegung erfordert eine bestimmte Git-

terstruktur, die sich nicht mit einer sinnvollen Vernetzung der Streben vereinbaren

lasst. Deshalb wurde fur die Rechnungen auf Streben verzichtet. Fur eine Wirkungs-

gradsteigerung des Antriebs ware es allerdings wichtig, die optimale Positionierung

der Streben und der Platte bestimmen zu konnen, die durch ihre Nahe zu den Schub


Abbildung 2.2: Schutzplatte mit Streben am Rumpf, aus Bartels und Jurgens (2006)

erzeugenden Teilen aufgrund gegenseitiger Beeinflussung einen deutlichen Einfluss

auf den Gesamtwirkungsgrad haben.

Das Overlapping Grid Verfahren soll theoretisch die flexible Kombination von beweg-

ten Flugeln in direkter Umgebung von festen Anbauteilen wie Streben und Schutz-

platte beherrschen. Ob es auch praktisch dazu in der Lage ist, soll diese Diplomarbeit

zeigen.

Kapitel 3

Numerische Stromungsmechanik

Die Eigenschaften der Bewegung von Flussigkeiten und Gasen, Stromung von Flui-

den genannt, lassen sich sehr gut durch partielle Differentialgleichungen beschreiben.

Diese konnen allerdings, abgesehen von wenigen Spezialfallen, nicht analytisch gelost

werden. Numerische Naherungslosungen schaffen hier Abhilfe: Die Anwendung einer

Diskretisierungsmethode liefert ein System algebraischer Gleichungen, das die Diffe-

rentialgleichung approximiert, welches dann mit Hilfe eines Computers gelost werden

kann.

Grundprinzip dieser Approximation ist die Berechnung von Losungen auf einer Dis-

kretisierung von Raum und Zeit. Die Wahl dieser Diskretisierung und der Berech-

nungsmethode beeinflusst sowohl die Qualitat der Losung, als auch den numerischen

Aufwand, um sie zu bestimmen. So fuhrt eine feinere Diskretisierung zu einer ge-

naueren Auflosung der Stromungsprozesse, allerdings erkauft man sich diese hohere

Detailtreue durch steigende Rechenzeiten, mit denen man unter Umstanden an die

Kapazitatsgrenzen heutiger Computer stoßt.

Einen Ausweg bieten gezielte Modellvereinfachungen, d.h. man verzichtet auf fur die

Losung unwesentliche Details, indem man weniger bedeutsame physikalische Phano-

mene vernachlassigt oder zusammenfasst. Auf diese Weise gewonnene Modelle der

Realitat mussen allerdings durch Experimente validiert werden. Durch geschickte

Modellierung ist mittlerweile die naherungsweise Losung einer sehr großen Anzahl

von Probleme moglich geworden.

Die Vielseitigkeit der modernen numerischen Stromungsmechanik naher zu beschrei-

ben geht weit uber den Rahmen dieser Arbeit hinaus. Denn obwohl alle Fluidbewe-

gungen im Prinzip den gleichen Gesetzen unterworfen sind, unterscheiden sich die

8

KAPITEL 3. NUMERISCHE STROMUNGSMECHANIK 9

Details der Modellierung einer Luftstromung mit Uberschallgeschwindigkeit uber ein

Flugelprofil wie bei Benek u. a. (1986) signifikant von der Simulation der Interakti-

on von Zellen in Blutgefaßen durch Chung u. a. (2007), bedingt vor allem durch

die unterschiedliche Dichte und Viskositat (Zahigkeit), aber auch andere physikali-

sche Eigenschaften der betrachteten Fluide. Griebel u. a. (1995) zeigen weitere un-

terschiedliche Modellansatze aus verschiedenen Bereichen der Stromungsmechanik,

Meakin (1995) gibt einige Anwendungsbeispiele aus Luft- und Raumfahrt.

Das fur diese Arbeit relevante Teilgebiet der Stromungsmechanik sind inkompres-

sible, newtonsche Flussigkeiten, also Fluide, deren Druck und Viskositat konstant

sind. Die mathematischen Grundlagen finden sich in Ferziger und Peric (2008), sie

werden in den folgenden Abschnitten kurz erklart und zusammengefasst: Erst wer-

den die dem physikalischen Modell zu Grunde liegenden Gleichungen eingefuhrt,

danach wird ein Uberblick uber einige in der Praxis verwendete Diskretisierungsver-

fahren gegeben und die Systemgleichungen mit dem Finite-Volumen-Verfahren (FV)

diskretisiert. Nach einem Abschnitt uber die Modellierung turbulenter Stromungen

wird zum Abschluss des Kapitels auf die Verfahren zur Losung des resultierenden

Gleichungssystems eingegangen.

3.1 Modellgleichungen

Die stromungsmechanischen Basisgleichungen basieren auf zwei Erhaltungsgesetzen,

die fur die Bewegung von Fluiden gelten: Massen- und Impulserhaltung.

Die fur Fluidstromungen geeignetste Form der Untersuchung ist die Betrachtung von

Kontrollvolumen (KV), d.h. meist ortsfesten Raumen mit konstanten Abmessungen,

fur die man die Zu- und Abflusse durch die Begrenzungsflachen sowie innere Senken

und Quellen bilanziert. Diese Methode, auch Eulerscher Ansatz genannt, umgeht

die bei Fluiden schwierige Verfolgung und Betrachtung einer festen Kontrollmasse

im Untersuchungsgebiet, dem in der Festkorperdynamik ublichen Lagrange’schen

Prinzip.

Das Erhaltungsgesetz einer Eigenschaft koppelt die Anderungsrate der betrachteten

Quantitat in einem Kontrollvolumen an außere und innere Einwirkungen: Die Varia-

tion einer Große wird durch das Netto aller Zu- und Abflusse durch die Oberflache

des KV aufgrund von Konvektion oder Diffusion, sowie durch Senken oder Quellen

im Inneren des KV bestimmt.


Im Folgenden werden die Erhaltungssatze fur ortsfeste Kontrollvolumen in koordi-

natenfreier Integral-Formulierung sowie die kartesische Koordinatenform der Impuls-

gleichung eingefuhrt, vgl. Hadzic (2005). Auf die formale Herleitung dieser Gleichun-

gen aus dem Reynolds-Transporttheorem wird an dieser Stelle verzichtet.

3.1.1 Massenerhaltung

Die Massenerhaltungsgleichung (auch: Kontinuitatsgleichung) besagt, dass die Ande-

rungsrate der Masse innerhalb eines Kontrollvolumens V gleich dem Nettofluss

durch die Volumenoberflache S ist. Die Gleichung lautet:

∂

∂t

∫

V

ρ dV +

∮

S

ρv · n dS = 0 (3.1)

Fur Ω ⊂ R3 und R

+ = t ∈ R, t > 0 bezeichnet ρ : Ω × R+ → R die Dichte

und v = (u1 , u2 , u3 )T : Ω × R+ → R

3 das Geschwindigkeitsfeld des Fluids, wobei

ui die Komponente in Koordinatenrichtung i bezeichnet. n : S × R+ → R

3 steht

fur das Feld der Normalen auf S . Die Dichte inkompressibler Flussigkeiten ist defi-

nitionsgemaß konstant und der erste Term in Gleichung (3.1) ist somit Null sofern

das Volumen konstant bleibt. Die Modellierung bewegter Gitter in Kapitel 5 kann je

nach verwendeter numerischer Methode zu einer Volumenanderung fuhren, was eine

Betrachtung des Terms notwendig macht.

Eine Anwendung des Gaußschen Integralsatz transformiert das Oberflachenintegral

in ein Volumenintegral:

∂

∂t

∫

V

ρ dV +

∫

V

div(ρv)dV = 0 (3.2)

Die Divergenz eines Vektorfeldes F: Rn → Rn ist fur ein kartesisches Koordinaten-

system wie folgt definiert:

div F = div (F1, . . . ,Fn) =n∑

i=1

∂

∂xi

Fi (3.3)

Mit dem Ubergang auf ein infinitisemal kleines Kontrollvolumen erhalt man die

koordinatenfreie Divergenzform der Massenerhaltungsgleichung:

∂ρ

∂t+ div(ρv) = 0 (3.4)


3.1.2 Impulserhaltung

Der Impulserhaltungssatz besagt, dass der Gesamtimpuls in einem abgeschlossenen

System konstant ist. Gibt es Zu- und Abflusse an Impuls, so sind sie gleich der Dichte

der Impulsquellen im Inneren des KV. Die Gleichung kann wie folgt geschrieben

werden:

∂

∂t

∫

V

ρv dV +

∮

S

ρvv · n dS =

∮

S

σ · n dS +

∫

V

ρ fb dV (3.5)

Hier steht σ fur den Spannungstensor der Oberflachenkrafte (Druck, Zug, Scher-

spannungen) und fb fur den Vektor der Korper- oder Volumenkrafte (Gravitations-,

Zentrifugal- und Corioliskrafte), die auf das Fluid wirken. Fur newtonsche, inkom-

pressible Flussigkeiten ist die Matrix des Spannungstensors wie folgt gegeben:

σ = µ [gradv + (gradv)T] − pI , (3.6)

wobei p den Druck und µ die dynamische Viskositat des Fluids, sowie I den Ein-

heitstensor bezeichnet.

Um Gleichung (3.6) losen zu konnen, muss der Tensor bezuglich eines bestimmten

Koordinatensystems aufgelost werden. Kartesische Koordinaten liefern eine einfa-

che Form der Impulsgleichung, sie ergibt sich aus Gleichung (3.5) und (3.6) durch

Skalarprodukt mit dem i. kartesischen Einheitsvektor ii :

∂

∂t

∫

V

ρui dV +

∮

S

ρuiv · n dS =

∮

S

µ grad ui · n dS +

∮

S

µ [iTi · (gradv)T] · n dS

−

∮

S

p iTi · n dS +

∫

V

ρ fbidV (3.7)

ui bezeichnet die i. Komponente des Geschwindigkeitsvektors v und fbidie i.

Komponente des Vektors der Korperkrafte.

Wie eben lasst sich Gleichung (3.5) durch Anwendung des Gaußschen Integralsatz

und Grenzubergang zu einem unendlich kleinen Volumenelement in die Form einer

koordinatenfreien Divergenzgleichung uberfuhren:

∂(ρv)

∂t+ div(ρvv) = div σ + ρ fb (3.8)

Fur die Divergenz des Spannungstensors gilt div σ = −grad p+µ∆v . Dies ergibt sich

aus (3.6) und div v = 0 , der Massenerhaltungsgleichung fur inkompressible Fluide.


Der Laplace-Operator ∆u =∑3

i=1 uxixiist hierbei komponentenweise zu verstehen.

Diese Eigenschaft und wiederholte Anwendung der Produktregel liefert nach einiger

Rechnung div(ρvv) = ρ(v · gradv) . Somit ergibt sich folgende Divergenzform der

Impulserhaltungsgleichung fur inkompressible, newtonsche Fluide:

ρ

(∂v

∂t+ v · gradv

)

= −grad p + µ∆v + ρ fb (3.9)

Im Falle kompressibler Stromungen ist daruber hinaus eine separate Energieerhal-

tungsgleichung notwendig, die sowohl thermische als auch mechanische Energie be-

handelt. Fur die hier betrachteten inkompressiblen, isothermen Stromungen ist nur

kinetische Energie von Bedeutung. Die Gleichung der Energie lasst sich aus dem

Skalarprodukt der Impulsgleichung mit dem Geschwindigkeitsvektor ableiten, ist al-

so kein selbstandiges Erhaltungsgesetz.

3.1.3 Universelle Transportgleichung

Die eben eingefuhrten Erhaltungsgleichungen, bekannt unter dem Namen inkom-

pressible Navier-Stokes-Gleichungen beschreiben Stromungen eines inkompressiblen,

newtonschen Fluids sehr genau.

Fur eine einfache mathematische Formulierung bietet sich die Einfuhrung einer uni-

versellen Gleichung an, die den Transport einer generischen, skalaren Variable φ

beschreibt:

∂

∂t

∫

V

ρφ dV +

∮

S

ρφv ·n dS =

∮

S

Γφgradφ ·n dS +

∮

S

qφS ·n dS +

∫

V

qφV dV (3.10)

φ bezeichnet die transportierte Variable, Γφ ihren Diffusions-Koeffizienten. qφS

steht fur Zu- und Abflusse durch die Oberflache und qφV fur Quellen und Senken.

Die in diesem Kapitel eingefuhrten Erhaltungsgleichungen (3.1) und (3.5) lassen sich

alle in Form von (3.10) schreiben. Dies vereinfacht in den kommenden Abschnitten

die numerische Diskretisierung, da nur die universelle Form der Transportgleichung

behandelt werden muss.

3.1.4 Randbedingungen und Anfangswerte

Zur Losung der Modellgleichungen benotigt man passende Randbedingungen und

Anfangswerte, damit das Problem mathematisch korrekt gestellt ist, d.h. eine ein-

deutige Losung existiert, die stetig von ihren Startwerten abhangt.


Die parabolische Natur der unstetigen Gleichungen macht einen vollstandigen Satz

Anfangswerte notwendig: Fur alle abhangigen Variablen muss ein Startwert im ge-

samten Berechnungsraum V vorgegeben werden:

φ(r, t0 ) = φ0 (r), r ∈ V (3.11)

Der bezuglich des Raumes elliptische Charakter der Gleichungen macht vollstandige

Randbedingungen an den Gebietsgrenzen SB notwendig. Die verschiedenen physi-

kalischen Randbedingungen, um Einlass, Auslass, Wand, Symmetrieebene, usw. fur

stromungsmechanische Probleme modellieren zu konnen, lassen sich mathematisch

in folgende zwei Klassen einteilen:

• Dirichlet-Randbedingungen, die den Wert einer abhangigen Variable auf einem

Teil des Randes SDB des Losungsgebietes angeben, d.h.

φ(rB , t) = fj(t), rB ∈ SDB . (3.12)

• Neumann-Randbedingungen, die den Gradienten einer abhangigen Variable auf

einem Teil des Randes SNB spezifizieren:

gradφ(rB , t) = fj(t), rB ∈ SNB . (3.13)

3.2 Finite Differenzen

Ansatz der Finite-Differenzen-Methode (FD) ist die Differentialform der zu losen-

den Gleichung. Das Losungsgebiet wird mit einem Gitter uberdeckt, anhand dessen

die Differentialgleichung nach dem Prinzip der numerischen Differentiation approxi-

miert wird: Die partiellen Ableitungen in jedem Gitterpunkt bzw. Knoten werden

ersetzt durch einen Approximationswert, der aus den Funktionswerten an umliegen-

den Gitterknoten gewonnen wird. Ergebnis ist fur jeden Gitterpunkt eine algebraische

Gleichung, die als Unbekannte die Werte des Knotens selbst, sowie die einiger Nach-

barn enthalt. Man kann eine Taylor-Reihenentwicklung oder einen Polynomansatz

benutzen, um Approximationen der Ableitungen an den Gitterpunkten zu erhal-

ten. Es ergibt sich ein System aus Differenzengleichungen, das unter Verwendung

der bekannten Anfangs- und Randwerte mit den Techniken gelost werden kann, die

in Abschnitt 3.7 behandelt werden. Uber Ferziger und Peric (2008) hinaus geben

Morton und Mayers (2005) eine umfassende Einfuhrung in die FD-Methodik.


Ein Vorteil dieser Methoden ist, dass sie bei strukturierten Gittern sehr einfach und

effektiv sind. Auch lassen sich durch die Wahl geeigneter Approximationsschemata

besonders bei regularen Gittern einfach Verfahren hoherer Ordnung konstruieren.

Bekannte Vertreter sind Einschrittverfahren wie Runge-Kutta, sowie die expliziten

Adams-Bashforth- oder impliziten Adams-Moulton-Methoden. Diese verwenden dis-

krete Werte mehrerer Zeitschritte, gehoren also zur Klasse der Mehrschrittverfahren.

Konvergenztheoretische Analysen auf regelmaßigen Gittern entsprechen haufig einer

Taylor-Restgliedbetrachtung. Die Struktur einiger Modellprobleme lasst auch eine

genauere Untersuchung mit Methoden der Fourier-Analysis zu.

Daruber hinaus ermoglicht das diskrete Maximumsprinzip einige Existenzaussagen.

Ein gravierender Nachteil dieses Diskretisierungsverfahrens ist, dass die Behandlung

komplexer Geometrien mit regelmaßigen Gittern sehr schwierig ist. Außerdem sind

die Verfahren ohne besondere Maßnahmen nicht konservativ, d.h es kann zur Ver-

letzung der Erhaltungsgleichungen kommen. Daruber kann fur viele FD-Verfahren

die Konvergenz gegen eine Losung nur durch starke Einschrankungen bezuglich der

Schrittweite garantiert werden.

Warmeleitung mit FD

Hier soll die Methodik am Beispiel der eindimensionalen Warmeleitungsgleichung

mit homogenen Dirichlet-Randbedingungen naher erlautert werden:

ut = uxx fur t > 0 und 0 < x < 1

u(0, t) = u(1, t) = 0 fur t > 0

u(x, 0) = u0(x) fur 0 ≤ x ≤ 1

(3.14)

Sei nun t0 = 0, . . . , tN eine gleichmaßige Diskretisierung der Zeit: tn = n∆t , sowie

x0 = 0, . . . , xJ = 1 eine gleichmaßige Diskretisierung des Raumes: xj = j∆x .

U(xj , tn) = Unj bezeichne die numerische Approximation der Differentialgleichung

an der Stelle xj zum Zeitpunkt tn .

Dann ist die Kombination aus Vorwarts-Differenzenquotient in der Zeit und zentra-

lem Differenzenquotient der zweiten Ortsableitung eine explizite Approximation der

Warmeleitungsgleichung:

Un+1j − Un

j

∆t=Un

j+1 − 2Unj + Un

j−1

(∆x)2(3.15)


Der Diskretisierungsfehler, der sich aus dem Restglied der Taylor-Entwicklung ergibt,

nimmt folgende Form an:

T (x, t) =1

2utt(x, η)∆t−

1

12uxxxx(ξ, t)(∆x)

2 (3.16)

Hier ist ξ ∈ (x − ∆x, x + ∆x), η ∈ (t, t + ∆t) . Das Verfahren ist also von ers-

ter Ordnung, sofern die Ableitungen beschrankt und genugend glatt sind und eine

Stabilitatsbedingung erfullt ist:

µ =∆t

(∆x)2≤

1

2(3.17)

Diese notwendige und hinreichende Bedingung ergibt sich aus der Formulierung der

Losung als Fourierreihe, vgl. Morton und Mayers (2005), Kapitel 2.7.

Die folgende Crank-Nicholson-Methode ist ein implizites Schema, das auf einer zen-

tralen Zeitdifferenz zum Zeitpunkt tn+1/2 und zentralem Differenzenquotient der

Ortsableitung basiert:

Un+1j − Un

j

∆t=

(

Un+1j+1 − 2Un+1

j + Un+1j−1

(∆x)2+Un

j+1 − 2Unj + Un

j−1

(∆x)2

)

(3.18)

Das Verfahren ist von zweiter Ordnung, wie eine Betrachtung der Taylor-Reihe des

Diskretisierungsfehlers zeigt:

Tn+1/2j = −

1

12

[(∆x)2uxxxx + (∆t)2uttt

]n+1/2

j+ · · · (3.19)

Daruber hinaus lasst sich bedingungslose Stabilitat und Konvergenz zeigen, vgl. Mor-

ton und Mayers (2005), Kapitel 2.10.

Im Unterschied zu expliziten Verfahren kann die Losung im impliziten Fall nicht

mehr direkt von den bekannten Randwerten ausgehend berechnet werden. Fur jeden

Zeitschritt muss ein dunn besetztes Gleichungssystem gelost werden, vgl. Abschnitt

3.7. Dafur entfallt die starke Restriktion der Koppelung der Ortsschrittweite an den

Zeitschritt.

3.2.1 Begriffe der Konvergenztheorie

An dieser Stelle sollen kurz einige nicht nur fur FD-Verfahren wichtige Begriffe der

Konvergenztheorie definiert werden:

Sei Lh eine FD-Approximation eines elliptischen Operators L mit Schrittweite h >

0 . Fur u ∈ C(Ω) sei Lhu die Auswertung des Differenzenoperators auf dem diskre-

tisierten Losungsraum Ωh .


Definition 3.1. I) Ein Differenzenverfahren Lh heißt konsistent mit der Gleichung

Lu = f wenn

Lu− Lhuh→0+

−→ 0 fur alle u ∈ C2(Ω)

II) Ein Verfahren hat Konsistenzordnung m, wenn

Lu− Lhu = O(hm) auf Ωh fur h→ 0 und fur alle u ∈ Cm+2(Ω)

III) Ein Differenzenverfahren heißt stabil, falls fur ein C > 0

‖L−1h ‖ ≤ C <∞ ∀h > 0

Konsistenz entspricht einer Fehlerabhangigkeit von der Schrittweite, sofern diese hin-

reichend klein gewahlt ist.

Stabilitat bedeutet, dass eine kleine Storung der Eingangswerte nur eine kleine

Storung der Losung nach sich zieht.

Aus Konsistenz und Stabilitat folgt Konvergenz, wie im jeweiligen Fall gezeigt werden

kann.

3.3 Finite Volumen

Der Finite-Volumen-Methode liegt die am Anfang dieses Kapitels eingefuhrte Inte-

gralform der Erhaltungsgleichung zu Grunde. Das Losungsgebiet wird in eine endliche

Anzahl sich nicht uberlappender Kontrollvolumina aufgeteilt, ein KV wird hierbei

von einer endlichen Anzahl Zellflachen begrenzt. In der Mitte eines jeden KV liegt der

Knoten, in dem die Variablenwerte berechnet werden. Werte auf der KV-Oberflache

werden aus den KV-Zentren durch Interpolation berechnet, Oberflachen- und Volu-

menintegrale werden mittels geeigneter Quadraturformeln approximiert. Man erhalt

so fur jedes KV eine algebraische Gleichung, die die Variablenwerte aus dem eigenen

Rechenknoten, sowie die von Nachbar-KVs enthalten. Die Losung des resultierenden

Gleichungssystems wird in Abschnitt 3.7 behandelt.

Ein Vorteil der FV-Methode ist die Freiheit in der Vernetzung, man kann jeden

Gittertyp verwenden und so auch komplexe Geometrien behandeln. Daruber hinaus

ist das Verfahren per Definition konservativ, solange man die Berechnung diffusiver

und konvektiver Flusse zwischen benachbarten KV uber ein geeignetes Schema ge-

meinsam durchfuhrt. Dann ist der Abfluss des einen KV durch die gemeinsame Seite


gleich dem Zufluss des anderen, und die Erhaltungsgleichungen sind erfullt. Die FV-

Methode ist außerdem anschaulich zu verstehen: Alle zu approximierenden Terme

haben physikalische Bedeutung.

Hauptnachteil der FV-Methodik ist, dass es im Vergleich zu FD sehr schwierig ist,

Verfahren von hoherer Ordnung als zwei zu entwickeln, besonders im Dreidimen-

sionalen. Da in einem FV-Verfahren numerische Interpolation, Differentiation und

Integration direkt ineinandergreifen, fuhrt die Konstruktion eines Verfahrens hoher-

er Ordnung zu komplexen Abhangigkeiten in den algebraischen Gleichungen, welche

die Losung des Gesamtsystems sehr aufwendig machen. Daruber hinaus existieren

aufgrund der Komplexitat der notigen Analysen nur sehr wenige Resultate der forma-

len mathematischen Konvergenztheorie. Die Validierung neuer Verfahren geschieht

haufig anhand von empirischen Untersuchungen.

3.4 Finite Elemente

Auch beim Finite-Elemente-Verfahren (FE) wird das Losungsgebiet in einen Satz fi-

niter Elemente aufgeteilt, die im Allgemeinen unstrukturiert sind. Anders als bei FD

oder FV, die die Differentialgleichung approximieren und Losungswerte an diskreten

Gitterpunkten liefern, sucht man beim FE-Ansatz aus einem (endlichdimensionalen)

Raum erlaubter Funktionen die Funktion, welche die Losung der Differentialglei-

chung am besten approximiert. Hierzu lost man ein Interpolationsproblem: Die ge-

suchte Losung u soll durch einfache Basisfunktionen des Ansatzraumes interpoliert

werden, die außer auf dem betrachteten Teilelement des Losungsgebietes uberall Null

sind. Aus einer Variationsformulierung der Erhaltungsgleichung ergibt sich ein Glei-

chungssystem fur die Interpolationskoeffizienten, die so genannte Systemmatrix, die

alle Funktionswerte an den Knoten korrekt interpoliert und gewisse Kontinuitatsbe-

dingungen an den Ubergangsknoten zwischen benachbarten Elementen erfullt.

Vorteil der FE-Diskretisierung ist die Moglichkeit, beliebige Geometrien zu behan-

deln. Gitterverfeinerungen sind leicht realisierbar, indem man die Elemente weiter

unterteilt. Mit den Mitteln der Funktionalanalysis ist die mathematische Analyse

des Verfahrens sehr gut moglich, fur bestimmte Typen von Gleichungen ist sie nach-

weislich die beste Methode.

Der Hauptnachteil, wie bei allen Verfahren mit unstrukturierten Gittern, ist die kom-

plizierte Struktur des auftretenden Gleichungssystems, die es schwierig macht, effizi-

ente Losungsmethoden zu finden. Daruber hinaus mussen die Erhaltungsprinzipien


der Stromungsprozesse uber zusatzliche Nebenbedingungen modelliert werden.

Eine umfassende Behandlung der Finiten Elemente findet sich z.B. bei Braess (2003).

Das Verfahren wird hier zwar nicht verwendet, jedoch kommt das zu Grunde liegen-

de Prinzip der Variationsformulierung einer PDE bei der Konvergenzanalyse der

Schwarz-Verfahren in Kapitel 4 zum Einsatz. Es soll deshalb an dieser Stelle ein-

gefuhrt werden.

3.4.1 Variationsformulierung

Die Losung einer PDE lasst sich auch uber ein verwandtes Problem der Variations-

rechnung bestimmen. Die Theorie dieser Variationsformulierung wird am Beispiel

folgender elliptischen Differentialgleichung zweiter Ordnung veranschaulicht:

Lu ≡ −div (a(x ) grad u) + c(x )u = f in Ω

u = 0 auf ∂Ω(3.20)

Hierbei sei Ω ⊂ Rd , und es gelte

0 < a0 ≤ a(x) ∀x ∈ Ω

c(x) ≥ 0 sei glatt, und f(x) ∈ L2(Ω) .

Durch Multiplikation mit einer Testfunktion v(x) , die die gleichen analytischen Ei-

genschaften wie die gesuchte Losung u hat, erhalt man eine neue Gleichung, die uber

das Losungsgebiet integriert wird. Uber partielle Integration lassen sich die hoheren

Ableitungen eliminieren und man erhalt die so genannte schwache Formulierung der

Differentialgleichung (3.20):

Bestimme u ∈ H10(Ω) als Losung von

A(u, v) = F (v) ∀v ∈ H10 (Ω) mit

A(u, v) ≡∫

Ω(a(x) grad u grad v + c(x )u v) dx

F (v) ≡∫

Ωfv dx

(3.21)

Wie allgemein ublich bezeichnet H10 (Ω) den Sobolev-Raum

H10 (Ω) ≡ v ∈ H1(Ω) : v = 0 auf ∂Ω (3.22)

und H1(Ω) den Hilbert-Raum

H1(Ω) ≡ v ∈ L2(Ω) : ‖v‖21,Ω <∞ mit Norm

‖v‖21,Ω ≡

∫

Ω(v2 + |grad v |2) dx

(3.23)


Definition 3.2. Sei H ein Hilbertraum. Eine Bilinearform A : H × H → R heißt

stetig, wenn fur ein C > 0 gilt:

|A(u, v)| ≤ C‖u‖ · ‖v‖ ∀u, v ∈ H (3.24)

Eine stetige Bilinearform A heißt koerziv oder elliptisch, wenn fur ein α > 0 gilt:

A(v, v) ≥ α‖v‖2 ∀v ∈ H (3.25)

Theorem 3.1. (Existenzsatz)

Sei L ein gleichmaßig elliptischer Differentialoperator 2. Ordnung. Dann hat Problem

(3.20) stets eine schwache Losung in H10 (Ω) . Diese ist Minimum des Variationspro-

blems

1

2A(v, v) − F (v) → min

in H10 (Ω) .

Beweis: Satz von Lax-Milgram, siehe Braess (2003), Kapitel II.2

Somit ist die Verbindung zwischen dem Variationsproblem und PDE hergestellt. Fur

PDE mit weiteren Termen ergibt sich aus der schwachen Formulierung - ahnlich wie

bei Lagrange-Optimierung unter Nebenbedingungen - eine Variationsformulierung in

der Gestalt eines Sattelpunktproblems.

Dies soll am Beispiel der Stokes-Gleichung demonstriert werden, die den dreidimen-

sionalen Fluss einer inkompressiblen, zahen Flussigkeit beschreibt:

∆u + grad p = −f in Ω

div u = 0 in Ω

u = u0 auf ∂Ω

(3.26)

u : Ω → R3 bezeichnet das Geschwindigkeitsfeld, p : Ω → R den Druck. Fur die

Sattelpunktformulierung setzt man:

X := H10 (Ω)3, M := L2,0(Ω) :=

q ∈ L2(Ω) :

∫

Ωq dx = 0

A(u,v) :=∫

Ωgradu gradv dx

B(v, q) :=∫

Ωdiv v · q dx

(3.27)

Hierbei ist gradu gradv :=∑

i ,j∂ui

∂xj

∂vi

∂xj


Fur homogene Dirichlet-Randbedingungen u0 = 0 ergibt sich folgendes Sattelpunkt-

problem:

Bestimme (u, p) ∈ X × M , so dass

A(u,v) + B(v, p) = F (v) fur v ∈ X

B(u, q) = 0 fur q ∈ M(3.28)

Auch in diesem Fall heißt ein so bestimmtes Losungspaar (u, p) schwache Losung des

Stokes-Problems. Siehe Braess (2003), Kapitel III.6 fur Aquivalenz und Losbarkeit

der Probleme.

3.4.2 Ritz-Galerkin-Verfahren

Fur die numerische Losung einer PDE geht man im Finite Elemente Ansatz folgen-

dermaßen vor: Man bestimmt die Losung des zu Grunde liegenden Variationspro-

blems nicht auf dem unendlichdimensionalen Funktionenraum H10 (Ω) , sondern auf

einem endlichdimensionalen Unterraum, dem Ansatzraum Sh .

uh ist Losung der Variationsaufgabe

1

2A(v, v) − F (v) → min|Sh

(3.29)

wenn gilt

A(uh, v) = F (v) ∀v ∈ Sh (3.30)

Fur eine Basis ψ1, ψ2, . . . , ψN von Sh ist dies aquivalent zu

A(uh, ψi) = F (ψi) fur i = 1, . . . , N (3.31)

Die eindeutige Darstellung der gesuchten Losung bezuglich der Basis von Sh

uh =

N∑

k=1

zkψk (3.32)

fuhrt zu einem Gleichungssystem in den Koeffizienten zk

N∑

k=1

A(ψk, ψi)zk = F (ψi) fur i = 1, . . . , N (3.33)

In Matrix-Vektorform ergibt sich fur Aik := A(ψk, ψi) und bi := F (ψi)

Az = b (3.34)


A heißt Steifigkeitsmatrix oder Systemmatrix des Problems, sie ist dunn besetzt und

kann mit den Losungsverfahren in Abschnitt 3.7.3 behandelt werden.

Es gibt verschiedene Ansatze, das Problem (3.29) zu losen.

• Rayleigh-Ritz-Verfahren zielen darauf ab, das Variationsproblem in Sh direkt

zu minimieren.

• Galerkin-Verfahren verwenden die schwache Formulierung basierend auf (3.31)

fur nicht notwendigerweise symmetrische Probleme. Im symmetrischen, positiv

definiten Fall (wie er sich fur einen elliptischen Operator A ergibt) spricht man

von Ritz-Galerkin-Verfahren.

• Petrov-Galerkin-Verfahren verwendet einen von Sh verschiedenen Raum von

Testfunktionen, was fur Probleme mit Singularitaten von Vorteil sein kann.

Es soll hier auch nicht weiter auf die Verfahren im Detail eingegangen werden, zum

Abschluss der FEM-Theorie wird noch kurz ein Beispiel fur einen moglichen Funk-

tionenraum Sh gegeben:

Abbildung 3.1: Stetige, stuckweise lineare Shape-Funktion, aus Bollhofer und Mehr-

mann (2004)

3.4.3 Zusammengesetzte, stuckweise lineare Basisfunktionen

Ausgehend von einer Diskretisierung von Ω ⊂ R2 mit kongruenten Dreiecken wahle

den Ansatzraum Sh wie folgt:

Sh := v ∈ C(Ω) : v linear in jedem Dreieck und v = 0 auf ∂Ω (3.35)


v ∈ Sh hat in jedem Dreieck die Form v(x, y) = a + bx + cy und ist durch die

Funktionswerte an den Gitterpunkten eindeutig bestimmt. Zu Knoten P1, . . . , Pm

auf Ω definiere die nodale Basis ψ1, . . . , ψN von Sh durch

ψi(Pj) = δij (3.36)

δij ist das Kronecker-Delta. Die so erhaltenen Basisfunktionen ψi heißen Hutfunk-

tionen oder Shape-Funktionen, sie sind nur in direkter Umgebung des Punktes Pi

von null verschieden, vgl. Abb. 3.1.

Die Theorie lasst sich in gleicher Form auf dreidimensionale Falle erweitern. Die

linearen Shape-Funktionen nehmen auf einem Tetraeder folgende Form an:

ψi(x, y, z) = a+ bx+ cy + dz (3.37)

Wie eben gilt ψi(Pj) = δij .

Bei Verwendung eines Funktionenraumes Sh auf Basis solcher stuckweise linearer

Shape-Funktionen zur Interpolation der gesuchten Losung u spricht man von linea-

rer Interpolation. Andere Basen des Ansatzraums (stuckweise quadratische Funk-

tionen, Polynome, Wavelets,..) sind auch denkbar, allerdings ist der rechnerische

Mehraufwand bei Shape-Funktionen hoherer Ordnung signifikant. Deshalb wird im

Dreidimensionalen haufig der lineare Ansatz gewahlt.

3.5 FV-Diskretisierung

Fur viele stromungsmechanische Untersuchungen hat sich die Finite-Volumen-Me-

thode unter anderem wegen ihrer Erhaltungseigenschaften etabliert. Hier soll nun

das diskrete FV-Modell eingefuhrt werden, das dieser Arbeit zu Grunde liegt: Es

wird das von Hadzic (2005) beschriebene Verfahren ubernommen, das auf einer Dis-

kretisierung des Raumes zweiter Ordnung mit impliziter Zeitschrittintegration fur

instationare Probleme beruht. Die Berechnung der gesuchten Variablen erfolgt se-

quentiell, indem die Gleichungen linearisiert und dann nacheinander mit Konjugierte-

Gradienten-Verfahren gelost werden. Fur die Koppelung von Druck und Geschwin-

digkeit und die Druckberechnung wird der SIMPLE Algorithmus benutzt.

Alle Erhaltungsgleichungen außer der Kontinuitatsgleichung, die getrennt behandelt

wird, nehmen folgende Form an, wenn sie auf ein diskretes Kontrollvolumen ∆V wie


Abbildung 3.2: Diskretes Kontrollvolumen, aus Hadzic (2005)

in Abbildung 3.2 angewandt werden:

∂

∂t

∫

∆V

ρφ dV

︸︷︷︸

Anderungsrate

+N∑

j=1

∫

Sj

ρφv · n dS

︸︷︷︸

Konvektion

=N∑

j=1

∫

Sj

Γφgradφ · n dS

︸︷︷︸

Diffusion

+

N∑

j=1

∫

Sj

qφS · n dS

︸︷︷︸

Oberflachenquellen

+

∫

∆V

qφV dV

︸︷︷︸

Volumenquellen

(3.38)

Die Oberflachenintegrale aus (3.10) stellen sich als Summe uber die einzelnen Zell-

flachen dar, die das KV begrenzen. Die unterschiedliche numerische Berechnung der

Terme in Gleichung (3.38) (Anderungsrate, Konvektion, Diffusion und Quellen) wer-

den im Folgenden beschrieben.

3.5.1 Approximation von Integralen

Fur die Integration uber ein KV wird die Mittelpunktsregel verwendet, das einfachste

Verfahren zweiter Ordnung:

ΨP0=

∫

∆V

ψ dV = ψ∆V ≈ ψP0∆V (3.39)

ψP0bezeichnet den Variablenwert im Zellmittelpunkt P0 , der mit dem Volumen

multipliziert wird.


Auch fur Oberflachenintegrale verwendet man die Mittelpunktsregel, der Fluss Fj

der transportierten Variable durch die Oberflache berechnet sich nach folgender For-

mel:

Fj =

∫

Sj

f · n dS ≈ fj · sj (3.40)

f steht fur den Flussvektor, fj fur den Wert von f im Schwerpunkt der Seiten-

flache Sj . Durch das Skalarprodukt mit sj , dem Oberflachenvektor der Seite j,

ergibt sich der Fluss durch Sj . Der Wert von fj steht nicht direkt zur Verfugung, er

muss aus den Knotenwerten durch Interpolation mindestens zweiter Ordnung berech-

net werden, sonst verliert die Mittelpunkt-Quadraturformel ihre Genauigkeit zweiter

Ordnung.

Der Flachenvektor einer KV-Seite berechnet sich wie folgt: Man zerlegt die Seite

so in Dreiecke, dass sie einen gemeinsamen Eckpunkt haben. Sei r1 der Ortsvektor

dieses Punktes. Dann gilt fur den Oberflachenvektor eines Dreiecks:

sD =1

2

Nv∑

i=3

[(ri−1 − r1) × (ri − r1)] (3.41)

Nv bezeichnet die Anzahl der Eckpunkte der Seite. Den Oberflachenvektor der Seite

sj erhalt man durch Aufsummieren der Flachenvektoren der einzelnen Dreiecke. sj

steht normal zur Seitenflache, der Betrag des Vektors entspricht dem Flacheninhalt

Sj der Seite.

Mit dem Gauss-Theorem lasst sich das Zellvolumen ∆V approximieren:

∆V =

∫

V

dV =

∫

V

∇ · (x i) dV =

∮

S

x i · n dS ≈∑

j

xj sxj (3.42)

sxj bezeichnet die x-Komponente des Oberflachenvektors, es gilt

sj = S jn = sxj i + sy

j j + szjk . (3.43)

3.5.2 Konvektion

Der konvektive Fluss durch die Zellflache j wird mit Mittelpunktregel wie folgt

approximiert:

F cj =

∫

Sj

ρφv · n dS ≈ φjmj (3.44)


φj ist der Wert der Variable im Schwerpunkt der Oberflache, mj bezeichnet den

Massenfluss durch die Oberflache j , er wird mittels Picard-Iteration aus dem vorhe-

rigen Zeitschritt approximiert, vgl. Ferziger und Peric (2008). Fur die Interpolation

von φj aus den Zellzentren stehen folgende Methoden zur Verfugung:

Lineare Interpolation

Eines der einfachsten Verfahren zweiter Ordnung beruht auf linearer Interpolation.

Abbildung 3.3: Lineare Interpolation in der Mitte der KV-Seite, aus Hadzic (2005)

Der Wert des Punktes j′ , der zwischen benachbarten Rechenknoten Pj und P0

liegt, berechnet sie wie folgt (vgl. Abb. 3.3):

φZj′ = φPj

λj + φP0(1 − λj) (3.45)

Der Interpolationsfaktor λj berechnet sich uber

λj =(rj − rP0

) · dj

dj · dj

. (3.46)

dj = rj−r0 ist der Vektor, der P0 mit Pj verbindet, rk bezeichnet den Ortsvektor

von k . Im Punkt j′ ist dies eine Approximation zweiter Ordnung. Falls j und j′

allerdings weit auseinander liegen, muss die Formel wie folgt angepasst werden, um

auch in j die gewunschte Interpolationsgute zu erhalten:

φZj = φ′

j + (gradφ)j ′(rj − rj′) (3.47)

Der Gradient in j′ wird durch Interpolation nach (3.45) aus den Gradienten in

den Zellzentren gewonnen. Die hier beschriebene lineare Interpolation wird auch


Zentraldifferenzenmethode genannt, da die Form mit der FD-Approximation der

ersten Ableitung ubereinstimmt.

In der Stromungsmechanik gibt es bestimmte Umstande, unter denen diese Methode

zu physikalisch inkorrekten, oszillierenden Losungen fuhrt. Man benotigt ein stabi-

leres Verfahren.

Upwind-Interpolation

Ein weitverbreitetes Verfahren, das beschrankte, oszillationsfreie Losungen liefert, ist

die Upwind-Interpolation. Der Wert der Variable auf der Oberflache wird mit dem

nachsten, stromaufwarts gelegenen Rechenknotenwert belegt:

φUj =

φP0: mj ≥ 0

φPj: mj ≤ 0

(3.48)

Das Verfahren ist allerdings nur von erster Ordnung, und es wurde gezeigt dass das

Verfahren numerisch diffusiv ist - Diskretisierungsfehler konnen zu unphysikalischer

Diffusion fuhren. Man benotigt ein sehr feines Gitter um Losungen von akzeptabler

Genauigkeit zu erhalten. Deshalb verwendet man in der Praxis haufig eine Mischung

beider Verfahren:

Gemischte Verfahren

Um die unzureichende Stabilitat des Zentraldifferenzverfahrens zu verbessern, kann

man der Interpolation einen gewissen Anteil des Upwind-Verfahrens ”beimischen”.

Der Wert in der Zelloberflache wird dann wie folgt berechnet:

φj = γ φZj + (1 − γ)φU

j (3.49)

Das Mischverhaltnis γ (engl. blending factor) nimmt einen Wert zwischen 0 und 1

an. Je nach Struktur des Problems muss eine andere Einstellung gewahlt werden um

einen Kompromiss aus Stabilitat und Genauigkeit zu erreichen.

Der Ansatz der verzogerten Korrektur liefert folgende Form der Gleichung:

φj = φUj + γ(φZ

j − φUj )alt (3.50)

Hier kommt nur der unbedingt stabile Teil aus dem Upwind-Verfahren fur die ak-

tuelle Iteration zu tragen, die Differenz zwischen den beiden Verfahren wird unter

Verwendung der Werte aus der vorherigen Iteration bei der Losung des Gleichungs-

systems explizit behandelt, vgl. Ferziger und Peric (2008).


3.5.3 Diffusion

Der diffusive Fluss FDj durch eine Zellflache j wird uber die Mittelpunktregel fol-

gendermaßen berechnet:

FDj =

∫

Sj

Γφgradφ · n dS ≈ Γφj (gradφ)∗j · sj = Γφj

(∂φ

∂n

)

j

Sj (3.51)

Γφ steht fur den Diffusionskoeffizienten, (gradφ)∗j fur eine Approximation des Gra-

dienten im Mittelpunkt der Seitenflache, Sj fur den Inhalt des Oberflachenstuckes.

Den Gradienten in Normalenrichtung aus Ableitungen im Zentrum der Zellen zu

interpolieren birgt die Gefahr von Oszillationen, vgl. Ferziger und Peric (2008), wes-

halb Hadzic (2005) folgende, von Muzaferija (1994) eingefuhrte Form verwendet, die

die Probleme verhindert und die Genauigkeit zweiter Ordnung wahrt:

(∂φ

∂n

)

j

≈φPj

− φP0

|dj|− (gradφ)altj ·

(dj

|dj|−

sj

|sj|

)

(3.52)

(gradφ)altj wird durch Interpolation nach Gleichung (3.45) bestimmt. Der erste Term

der rechten Seite ist eine Zentraldifferenzen-Approximation der Ableitung in direkter

Richtung zwischen den Zellmittelpunkten, der zweite Term korrigiert diese bezuglich

der Normalenrichtung sj zur Zelloberflache.

3.5.4 Gradientenberechung im Zellmittelpunkt

Gradienten im Zellmittelpunkt konnen mit Mittelpunktregel und dem Gaußschen

Integralsatz wie folgt approximiert werden:

∫

V

gradφ dV =

∮

S

φ · n dS ⇒ (gradφ)P0≈

∑

j φjsj

∆VP0

(3.53)

Die Approximation ist zweiter Ordnung und kann auf beliebige Zellformen angewandt

werden.

Eine andere Moglichkeit einer Gradientenapproximation zweiter Ordnung ist die Me-

thode der kleinsten Quadrate: (gradφ)P0kann berechnet werden, indem man eine

lineare Ansatzfunktion an einen Satz Nachbarpunkte von P0 anpasst. Hierzu lost

man die folgende Gleichungen:

dj · (gradφ)P0= φPj

− φP0(j = 1 , . . . ,Nj ) (3.54)


Hierbei ist dj der Verbindungsvektor von P0 zu Pj. Das uberbestimmte Gleichungs-

system wird mit der Methode der kleinsten Quadrate gelost, es ergibt sich:

(grad φ)P0= D−1

∑

j

dTj (φPj

− φP0) (3.55)

Die Matrix D berechnet sich folgendermaßen: D =∑

j

dTj dj . Sie ist symmetrisch,

und ihre Eintrage hangen nur von der geometrischen Gitterstruktur ab, folglich muss

sie fur ein gegebenes Gitter nur einmal berechnet werden.

3.5.5 Quellterme

Integration der Quelldichte qφV einer Variable φ uber das Kontrollvolumen ∆VP0

wird mit (3.39) durchgefuhrt:

QφV =

∫

V

qφV dV ≈ (qφV )P0∆VP0

(3.56)

Genauso wird eine Oberflachenquelle als Integral der Quelle qφS uber die KV-Ober-

flache SP0unter Verwendung von Gleichung (3.40) approximiert:

QφS =

∮

S

qφS · n dS =

Nj∑

j=1

∮

Sj

qφS · n dS ≈

Nj∑

j=1

qφSj· sj (3.57)

3.5.6 Behandlung der Druckgleichung

Der Druckterm aus der Impulsgleichung kann auf zwei Arten behandelt werden.

Zum einen kann der Druck berechnet werden, indem man die Druckkrafte auf allen

Zellflachen berechnet und aufsummiert. Die andere Moglichkeit ist, das Oberflachen-

integral mit dem Gaußschen Integralsatz in ein Volumenintegral zu uberfuhren, das

den Druckgradienten enthalt:

∮

S

pI · n dS =

∫

V

grad p dV (3.58)

Die zweite Methode wird hier ubernommen, sie ist konservativ, solange der Druck-

gradient mit Gleichung (3.53) berechnet wird.


3.5.7 Zeitintegration

Fur instationare Stromungen muss die Integration von Gleichung (3.38) auch bezuglich

der Zeit erfolgen. Hierzu wird sie folgendermaßen umgeformt:

∂Ψ

∂t= f(t, φ) mit Ψ =

∫

V

ρφ dV ≈ ρφP0∆VP0

und φ = φ(r, t) (3.59)

Es gibt zwei Klassen von Zeitschrittverfahren: explizite und implizite. Explizite Me-

thoden sind leicht zu implementieren, da keine linearen Gleichungssysteme gelost

werden mussen, und sie erlauben eine beliebige Genauigkeit in der Zeit. Allerdings

mussen aus Stabilitatsgrunden die Zeitschritte oft viel kleiner gewahlt werden, als es

fur die Genauigkeit notig ware, was sie fur die praktische Anwendung oft ineffizient

macht. Abhilfe schaffen implizite Verfahren, sie sind unbedingt stabil und erlauben

beliebige Zeitschrittweiten, die nur durch die gewunschte Genauigkeit beschrankt

sind, allerdings muss man algebraische Gleichungssysteme losen, die nichtlinear sein

konnen. Hier werden zwei (voll-)implizite Methoden betrachtet: Die Euler-Methode,

ein Verfahren erster Ordnung, und die Drei-Ebenen-Methode zweiter Ordnung.

Das Euler-Schema erhalt man, indem man Gleichung (3.59) uber das Zeitintervall

∆t = tn − tn−1 integriert:

1

∆t

tn∫

tn−1

(∂Ψ

∂t

)

dt ≈Ψn

P0− Ψn−1

P0

∆t= f n (3.60)

Das Euler-Schema ist sehr einfach und unbedingt stabil, es hat Konsistenz- und Kon-

vergenzordnung 1. Allerdings kann seine Struktur aufgrund von numerischer Diffu-

sion zu einer starken Dampfung der Eigenschaften instationarer Stromungen fuhren.

Um die Genauigkeit zu erhohen ist in manchen Situationen ein Verfahren hoherer

Ordnung notig. Hadzic (2005) verwendet die Drei-Ebenen-Methode, die Variablen-

werte von drei Zeitschritten erfordert und die Integration uber ∆t gemittelt um

t durchfuhrt, also von t − ∆t/2 bis t + ∆t/2 . Die bei Ferziger und Peric (2008)

beschriebene quadratische Ruckwartsintegration in der Zeit ergibt die Approximati-

onsformel zweiter Ordnung in tn :

1

∆t

tn+1/2∫

tn−1/2

(∂Ψ

∂t

)

dt ≈3Ψn

P0− 4Ψn−1

P0+ Ψn−2

P0

2∆t= f n (3.61)

Ein wichtiger Parameter fur die Wahl der Schrittweite ∆t ist die so genannte Cou-

rant-Zahl, die beschreibt um wie viele Zellen sich eine betrachtete Große pro Zeit-


schritt maximal fortbewegt:

Co =u∆t

∆x(3.62)

Fur viele Anwendungen fordert man Co ≤ 1 .

3.5.8 Resultierende algebraische Gleichungen

Summiert man alle Terme der Gleichung (3.10) wie beschrieben auf, so erhalt man

fur jedes KV eine algebraische Gleichung, die den Wert der abhangigen Variable φ

im KV-Zentrum mit den Werten der Nachbarvariablen in Verbindung setzt:

aP0φφP0−

Nj∑

j=1

ajφPj= bP0

(3.63)

j lauft uber die Indices der Nachbarknoten, die in den impliziten Integralapproxi-

mationen vorkommen, bP0umfasst die Quellterme, Teile der Zeitableitung und die

Teile der konvektiven und diffusiven Terme, die nach der Methode der verzogerten

Korrektur explizit behandelt werden, sowie die Oberflachenquellterme. Die Werte

der Koeffizienten ergeben sich aus den oben eingefuhrten Approximationsformeln:

aj = Γφ|sj|

|dj|− min (mj , 0)

aP0φ = atP0

+

Nj∑

j=1

aj (3.64)

bP0= QφV +QφS + btP0

− γ mjφZj − [max (mj , 0)φP0

− min (mj , 0)φPj]

+ Γφ(gradφ)j ·

(

sj −|sj|

|dj|dj

)

+∑

B

aBφB

atP0

und btP0entsprechen dem transienten Teil des verwendeten Zeitintegrationssche-

mas, die Werte findet man in Tabelle (3.1). Der Upwind-Anteil der Integralapproxi-

mation ergibt sich aus folgender Beziehung:

mjφAj = max (mj , 0)φP0

+ min (mj , 0)φPj(3.65)

Die Summe∑

B

aBφB lauft uber die Oberflachen, die den Rand des Berechnungsge-

bietes beschreiben, φB bezeichnet den Variablenwert am Rand. Des Weiteren werden

die Koeffizienten und Quellterme in Gleichung (3.64) im Einklang mit dem Picard-

Iterationsschema mit den Werten des vorangehenden Zeitschritts berechnet.


Tabelle 3.1: Zeitintegrationskoeffizienten in Gleichung (3.64)

Koeffizient implizites Euler Drei-Ebenen-Methode

atP0

ρ∆VP0

∆t

3ρ∆VP0

2∆t

btP0

ρ∆VP0φn−1

P0

∆t

4ρ∆VP0φn−1

P0− ρ∆VP0

φn−2P0

2∆t

3.5.9 Druckberechnung mit dem SIMPLE-Algorithmus

Der Berechnung des Druckes kommt bei inkompressiblen Flussen eine besondere

Bedeutung zu. Er kann nicht wie die abhangigen Variablen direkt aus einer der

Systemgleichung berechnet werden. Hadzic (2005) wahlt den SIMPLE-Algorithmus,

ein iteratives Verfahren, das erst die mit Druckwerten aus einer vorherigen Itera-

tion linearisierten Impulsgleichungen lost. Das Massenungleichgewicht, das hierbei

in der Kontinuitatsgleichung entsteht, wird verwendet, um eine Druckkorrektur zu

berechnen, mit der die resultierenden Geschwindigkeiten die Kontinuitatsgleichung

erfullen. Da durch die korrigierten Geschwindigkeiten die Impulsgleichungen verletzt

sind, wird der Prozess iterativ wiederholt bis beide Gleichungen ausreichend genau

erfullt sind.

Geschwindigkeit an der Zelloberflache

Folgende Form der Geschwindigkeitsinterpolation an der KV-Oberflache unterdruckt

mogliche Oszillationen:

v∗

j = vj −

(∆VP0

avP0

)

j

[pPj

− pP0

|dj|−

grad p · dj

|dj|

]|dj|sj

dj · sj

(3.66)

Hier ist vj die lineare Interpolation nach Gleichung (3.45), der Rest ein Korrek-

turterm dritter Ordnung, vgl. Ferziger und Peric (2008) fur Ursprung und nahere

Beschreibung.


Pradiktor-Phase

Mit gegebenem KV-Oberflachengeschwindigkeiten lasst sich der Massenfluss uber

m∗

j = ρv∗

j · sj berechnen, und die in der Regel verletzte Massenerhaltungsgleichung

liefert:

Nj∑

j=1

m∗

j = ∆m (3.67)

∆m bezeichnet das Massenungleichgewicht, das nach erfolgter Korrektur Null sein

sollte. Der SIMPLE-Algorithmus liefert die notwendige Flusskorrektur:

m′

j = ρv′

j · sj = −ρ|sj|

(∆VP0

avP0

)

j

(∂p ′

∂n

)

j

≈ −ρ|sj|

(∆VP0

avP0

)

j

p ′

Pj− p ′

P0

dj · sj

|sj| (3.68)

Hierbei muss folgende Bedingung an die korrigierten Massenflusse erfullt sein:

Nj∑

j=1

m′

j + ∆m = 0 (3.69)

Es ergibt sich eine Druckkorrekturgleichung folgender Form:

ap′Pp′P0

−

Nj∑

j=1

ap′jp′Pj

= bp′ (3.70)

Fur die Koeffizienten gilt:

ap′j= ρ

sj · sj

dj · sj

(∆VP0

avP0

)

j

ap′

P=

Nj∑

j=1

ap′

jbp′ = −∆m = −

Nj∑

j=1

m∗

j (3.71)

Korrektor-Phase

Nach der Losung von Gleichung (3.70) korrigiert man den Massenfluss:

mj = m∗

j + ap′j(p′Pj

− p′P0) (3.72)

Dann berechnet man ein neues Geschwindigkeitsfeld uber

vP0= vpred

P0+ v′

P0= vpred

P0−

∆VP0

avP0

(grad p ′)P0(3.73)

Anschließend aktualisiert man die Druckwerte:

pP0= ppred

P0+ αpp

′

P0(3.74)

Der Unterrelaxationsfaktor αp wird fur stationare Probleme mit Werten von 0,1−0,3

und fur instationare Probleme mit Werten bis zu 1 belegt.


3.5.10 Randbedingungen

Fur den Abschluss des algebraischen Gleichungssystems, das aus der Diskretisierung

hervorgegangen ist, benotigt man noch Bedingungen an den Zellflachen, die das

Berechnungsgebiet begrenzen, sie ergeben sich wie folgt:

Am Einstromrand (engl. inlet) sind die Variablenwerte in der Regel bekannt, wes-

halb Dirichlet-Bedingungen verwendet werden. Konvektive Flusse konnen so direkt

berechnet werden, fur diffusive Flusse benotigt man einseitige FD-Approximationen.

Fur den Ausstromrand (engl. outlet) sind die Variablenwerte nicht bekannt, in der

Regel werden sie extrapoliert. Einzig die Geschwindigkeitswerte bedurfen besonderer

Behandlung, da die Massenerhaltung gewahrt bleiben muss. Dazu werden die extra-

polierten Werte mit dem Verhaltnis mI/mO der gesamt einstromenden zur gesamt

ausstromenden Masse multipliziert.

In einer Symmetrieebene muss die Geschwindigkeitskomponente normal zur Ebene

null sein, vBn(rB , t) = 0 , wie auch die konvektiven Flusse aller Skalargroßen. Dies

erfordert teils Dirichlet- und teils Neumann-Randbedingungen.

Fur eine Wand mit Haftbedingungen werden die Dirichlet-Geschwindigkeitswerte

an der Wand direkt auf das Fluid ubertragen, die ubrigen skalare Großen konnen

entweder durch Dirichlet oder Neumann beschrieben werden.

Die Randbedingungen der Druckkorrekturgleichung (3.70) mussen gesondert behan-

delt werden: In Elementen am Rand, in denen direkt Massenfluss beschrieben wird,

muss die Korrektur null sein, realisiert durch Neumann-Bedingungen. Ist hingegen

eine Druckrandbedingung gegeben, so wird die Nullkorrektur am Rand durch Dirich-

let-Bedingungen umgesetzt.

3.6 Turbulenzmodellierung

Das eingefuhrte Modell der Navier-Stokes-Gleichungen modelliert die Stromungen

eines inkompressiblen Fluids exakt. Allerdings stoßt man bei der numerischen Be-

rechnung auf folgende Problematik: Um alle Details einer Stromung, besonders die

Bildung kleinster Wirbel erfassen zu konnen benotigt man eine sehr feine Diskretisie-

rung des Losungsraumes. Diese so genannte direkte numerische Simulation (DNS) ist

sehr rechenaufwendig. Da fur die meisten stromungsmechanischen Problemen nicht

die Details der kleinsten Wirbel, sondern nur die quantitativen Eigenschaften dieser


turbulenten Stromungen von Bedeutung sind, bietet sich eine Modellvereinfachung

an.

3.6.1 RANSE

Die bekannteste Vereinfachung sind die Reynolds-gemittelten Navier-Stokes-Glei-

chungen (RANSE, engl. kurz fur Reynolds-averaged Navier-Stokes Equations). Die

Reynolds-Mittelung basiert darauf, dass statistisch stationare Stromungen als Sum-

me aus einem zeitgemittelten Wert und Schwankungen um diesen Wert dargestellt

werden konnen, vgl. Ferziger und Peric (2008):

φ(xi, t) = φ(xi) + φ′(xi, t) (3.75)

wobei

φ(xi) = limT→∞

1

T

T∫

0

φ(xi, t)dt (3.76)

t ist die Zeit, T das Mittelungsintervall, das hinreichend groß gewahlt werden muss,

um Schwankungen sinnvoll auszugleichen, so dass fur die mittlere Schwankung φ′ = 0

gilt.

Fur nicht statistisch stationare Stromungen muss eine andere Form der Reynolds-

Mittelung gewahlt werden:

φ(xi) = limN→∞

1

N

N∑

n=1

φn(xi, t) (3.77)

N bezeichnet die Anzahl der Ensemblemitglieder, die groß genug sein muss, um die

Fluktuationseffekte zu beseitigen. Diese Ensemblemittelung kann auf jede Stromungs-

art angewendet werden.

Bei der Mittelung nichtlinearer Ausdrucke entstehen Zusatzterme, bei quadratischen

Termen fuhrt das zu Produkten der Mittelwerte und der Kovarianz:

uiφ = (ui + u′i)(φ+ φ′) = uiφ+ u′iφ′ (3.78)

Es ergibt sich folgende Form der gemittelten Kontinuitats- und Impulsgleichungen fur

inkompressible Stromungen ohne Massenkrafte in Tensornotation und kartesischen


Koordinaten:

∂(ρui)

∂xi

= 0 (3.79)

∂(ρui)

∂t+

∂

∂xj

(ρuiuj + ρu′iu

′

j

)= −

∂p

∂xi+∂τ ij

∂xj(3.80)

τ ij bezeichnet den mittleren, viskosen Spannungstensor:

τ ij = µ

(∂ui

∂xj

+∂uj

∂xi

)

(3.81)

µ steht fur die molekulare Viskositat des Fluids.

Fur die gemittelte Gleichung einer skalaren Große ergibt sich:

∂(ρui)

∂t+

∂

∂xj

(ρujφ+ ρu′jφ

′

)=

∂

∂xj

(

Γ∂φ

∂xj

)

(3.82)

Besondere Beachtung gelte den Reynolds-Spannungen ρu′iu′

j und dem turbulenten

Skalarfluss ρu′iφ′ . Diese Terme konnen nicht durch die gemittelten Einzelgroßen

ausgedruckt werden. Man spricht von einem nicht geschlossenen Problem - es wer-

den zusatzliche Gleichungen benotigt, um die unbekannten Terme zu berechnen.

Man verwendet verschiedene Turbulenzmodelle um das RANSE-System zu schließen

- Approximationen der Unbekannten durch Funktionen der gemittelten Großen und

empirischer Parameter.

3.6.2 k- ε -Turbulenzmodell

Den Simulationen in dieser Arbeit liegt das k−ε -Turbulenzmodell zu Grunde, das

am weitesten verbreitete Wirbelviskositatsmodell.

Turbulenz wird bei Wirbelviskositatsverfahren uber eine (lokal) erhohte Viskositat

des Fluids in Bereichen nicht laminarer Stromung modelliert. Hierzu fuhrt man tur-

bulente Diffusivitats- und Viskositatskoeffizienen ein. Dies sind nichtlineare Funktio-

nen der Flussparameter (Diffusivitat, Geschwindigkeit, Viskositat,...), die im Stromungs-

raum um mehrere Großenordnungen variieren konnen und so den unterschiedlichen

Auspragungen der Turbulenz in der Stromung Rechnung tragen.

Fur die Reynolds-Spannungen gilt:

−ρu′iu′

j = µt

(∂ui

∂xj+∂uj

∂xi

)

−2

3ρδijk (3.83)


Fur den turbulenten Skalarfluss gilt entsprechend:

−ρu′jφ′ = Γt

∂φ

∂xj(3.84)

k steht fur die turbulente kinetische Energie, µt fur die Wirbelviskositat, Γt fur

die Wirbeldiffusivitat, sie mussen fur den Abschluss des RANSE-Systems uber die

Modellapproximationen berechnet werden.

Im hier verwendeten k- ε -Modell gilt:

k =1

2u′iu

′

i =1

2

(u′xu

′

x + u′yu′

y + u′zu′

z

)(3.85)

Die Wirbelviskositat berechnet sich unter Verwendung der dimensionslosen Modell-

konstante Cµ = 0.09 uber folgende Beziehung:

µt = ρCµk2

ε(3.86)

ε bezeichnet die Dissipationsrate.

Die Berechnung des RANSE-Systems erfolgt genau wie die normalen Navier-Stokes-

Gleichungen, allerdings unter Verwendung der effektiven Viskositat µeff = µ + µt ,

fur die zusatzlich der Transport der beiden Skalargroßen k und ε berechnet werden

muss. Diese Transportgleichungen, die nach dem Schema der ubrigen Skalargleichun-

gen diskretisiert werden, finden sich in Ferziger und Peric (2008), Kapitel 9.4, oder

Comet User Manual (2001), Kapitel 3.5.

Die auf diese Weise gewonnene Vereinfachung der Navier-Stokes-Gleichungen ist fur

viele Problemstellungen mehr als ausreichend. Allerdings kann diese Turbulenzmo-

dellierung nie alle Eigenschaften einer nicht laminaren Stromung erfassen, sie sollte

daher stets als Approximation gesehen werden. Eine ausfuhrliche Beschreibung des

Modells sowie anderer in der Praxis verwendeter Verfahren findet sich in Ferziger

und Peric (2008).

3.7 Losen des Gleichungssystems

Die in diesem Kapitel beschriebene Diskretisierung fuhrt zu algebraischen Gleichun-

gen der Form (3.63) fur jede Zelle und jede Variable. Die Gleichungen sind nichtli-

near und gekoppelt, d.h die Koeffizienten hangen von den Variablen ab und in den

Gleichungen steht mehr als eine abhangige Variable. Hadzic (2005) wahlt folgendes

Verfahren der sequentiellen Losung:


3.7.1 Sequentielles Losungsverfahren

Das Losungsverfahren beruht darauf, die gesuchten Variablen iterativ und nach-

einander zu berechnen. Hierfur wird das System mit bekannten Werten aus einer

vorherigen Iteration linearisiert und kann wie folgt geschrieben werden:

AφΦ = bφ (3.87)

Aφ bezeichnet die Koeffizientenmatrix mit N × N Eintragen fur die N Zellen. Φ

ist der Vektor der Variablenwerte φ und bφ der Vektor der Quellterme. Die Matrix

Aφ hat zwei wichtige Eigenschaften:

• Sie ist dunn besetzt, nur Nj + 1 Eintrage pro Zeile sind ungleich Null, mit Nj

der Anzahl in Gleichung (3.63) verwendeter Nachbarpunkte von Pj .

• Sie ist diagonaldominant, d.h. fur alle Zellen gilt:

|aP0φ| ≥∑

j

|ajφ| (3.88)

Daher bieten sich iterative Loser fur das lineare Gleichungssystem (LGS) an, die

die Dunnbesetztheit der Matrix bewahren. Eine exakte Losung ist nicht unbedingt

notig, da durch die Linearisierung ohnehin nur Naherungen der Koeffizienten berech-

net werden. Man iteriert solange, bis ein Abbruchkriterium erfullt ist (Reduktion der

Residuen um eine gewisse Großenordnung, Erreichen einer maximalen Zahl von Ite-

rationen). Iterationen fur die Losung des linearisierten Systems werden auch innere

Iterationen genannt.

Nach der Losung der linearisierten Systeme fur alle Variablen werden die Koeffizi-

enten des Gesamtsystems aktualisiert und der gesamte Vorgang wiederholt, bis ein

Konvergenzkriterium erfullt ist oder die maximale Anzahl dieser außeren Iterationen

erreicht ist. Das Konvergenzkriterium ist erfullt, wenn die entsprechend normierte

Summe der Residuen eine gewisse Schranke unterschreitet. Hadzic (2005) wahlt fol-

gende Abbruchbedingung:

Rkφ =

N∑

i=1

|rkφ|

Mφ=

N∑

i=1

|bP0− aP0φφP0

+Nj∑

j=1

ajφPj|

Nj∑

j=1

|aP0φφP0|

≤ ε (3.89)

Hier steht rkφ fur das Residuum im Punkt P0 in Iteration k , ε ist ein Bruchteil

von R0φ . In der Regel reduziert man das Residuum um drei Großenordnungen, also

bis Rkφ < 0,001R0

φ , es kann aber auch ein strengeres Abbruchkriterium notig sein

um eine konvergierte Losung zu erhalten.


3.7.2 Unterrelaxation

Eine verlangsamte Korrektur der Variablenwerte, Unterrelaxation, tragt zur Stabi-

litat des Losungsverfahrens des gekoppelten Gleichungssystems bei. In der außeren

Iteration wird nur ein Bruchteil der Anderung aus den inneren Iterationen uber-

nommen, starke Schwankungen wie sie besonders zu Anfang auftreten werden so

gedampft. Man setzt:

φmP0

= φm−1P0

+ αφ

(φneu

P0− φm−1

P0

)(3.90)

Der Unterrelaxationsfaktor αφ nimmt Werte zwischen 0 und 1 an. m ist der Index

der außeren Iterationen, φneuP0

ergibt sich wie folgt aus der Losung von Gleichung

(3.63):

φneuP0

=

bP0−

Nj∑

j=1

ajφmPj

aP0

(3.91)

3.7.3 LGS-Loser

Es gibt zwei Klassen von LGS-Losern, direkte und iterative. Direkte Loser bestimm-

ten in einer finiten Anzahl von Rechenoperationen eine exakte Losung, iterative Ver-

fahren beginnen mit einer Startlosung (Nullvektor, Losung einer vorherigen Iterati-

on), die sie in kleinen Schritten verbessern. Direkte Verfahren sind wegen der Große

der auftretenden Gleichungssysteme fur die meisten Stromungsmechanikprobleme

ungeeignet. Außerdem basieren die meisten Verfahren zum Losen von gekoppelten,

nichtlinearen Gleichungen auf linearisierten Approximationen, was die recheninten-

sive exakte Losung der LGS aus den inneren Iterationen unnotig macht. Aus der

Klasse der iterativen Verfahren wahlt Hadzic (2005) das Bi-CGStab-Verfahren mit

unvollstandiger Cholesky-Vorkonditionierung, das hier kurz ins Feld der LGS-Loser

eingeordnet werden soll.

Vorkonditionierer

Fur iterative Losungsverfahren großer, dunnbesetzer Systeme ist die Verwendung ei-

nes Vorkonditionierers ublich. Man bestimmt eine Vorkonditionierungsmatrix P um

die numerischen Eigenschaften der Systemmatrix zu verbessern. Das abgewandelte

System

P−1Ax = P−1b (3.92)


benotigt eine stark reduzierte Anzahl von Losungsiterationen, da es in der Regel

eine geringe Kondition und eine geschicktere Eigenwertverteilung hat. Ubliche Vor-

konditionierer sind einige Schritte eines iterativen LGS-Loser wie Jacobi- oder Gauss-

Seidel, Alternativen sind Zerlegungsverfahren wie Cholesky, das die Matrix in eine

untere und obere Teilmatrix zerlegt. Vergleiche Quarteroni u. a. (2002) fur einen

Uberblick der Verfahren.

Die meisten dieser Verfahren sind fur symmetrische, positiv definite Matrizen ausge-

legt, was auf ein Gleichungssystem zur Diskretisierung eines Stromungsproblems auf

einem unregelmaßigen Gitter in der Regel nicht zutrifft. Hierfur werden modifizier-

te Verfahren wie die unvollstandige Cholesky-Faktorisierung (engl. ILU (incomplete

LU), oder IC (incomplete Cholesky)) verwendet:

(LU)−1Ax = (LU )−1b , (3.93)

mit LU ≈ A eine Naherung der Systemmatrix, bestehend aus einer oberen und

unteren Dreiecksmatrix U und L. Die ILU-Zerlegung, die die Besetzungsstruktur

des Ausgangssystems bewahrt, existiert nicht fur alle invertierbaren Matrizen, fur

die hier auftretenden diagonaldominanten Systeme wurde die Existenz allerdings

gezeigt, vgl. Manteuffel (1980).

CG-Verfahren

Das Verfahren der konjugierten Gradienten (engl. CG, conjugate gradient) aus der

Klasse der Krylow-Unterraum-Verfahren ist ein effizientes Verfahren zum Losen großer,

dunnbesetzter Gleichungssysteme mit symmetrischer und positiv definiter System-

matrix. Es basiert auf der Bestimmung der Losung in einem Unterraum, der in

jeder Iteration um eine Dimension erweitert wird und theoretisch in m Schritten die

endgultige Losung eines m-dimensionalen Problems liefert. In der Praxis fuhren z.B.

Rundungsfehler dazu, dass das Verfahren in m Schritten nicht vollstandig konver-

giert. Es sollte daher als iteratives Verfahren zur Bestimmung von Naherungslosun-

gen gesehen werden. Gerade hierbei fuhrt wie bei allen Krylow-Methoden die Ver-

wendung eines Vorkonditionierers wie IC zu einer deutlich erhohten Konvergenzge-

schwindigkeit. Dieses ICCG-Verfahren und einige andere werden in Quarteroni u. a.

(2002) naher beschrieben.

Die CG-Verfahren sind nur fur symmetrische, positiv definite Matrizen sicher kon-

vergent. Zur Behandlung von unsymmetrischen Fallen gibt es allerdings Erweite-

rungen wie die Bi-konjugierte Gradientenmethode (BiCG) und ihre Stabilisierung,


BiCGStab, die auf van der Vorst (1992) zuruckgeht und in Kombination mit IC-

Vorkonditionierung fur die Rechnungen dieser Arbeit verwendet wird.

3.7.4 Losungsalgorithmus

Die Schritte zur Losung eines stromungsmechanischen Problems lassen sich wie folgt

zusammenfassen:

1. Alle Felder mit Anfangswerten belegen

2. Die linearisierten Impulsgleichungen losen um Werte fur die Geschwindigkeit

zu erhalten

3. Die Druckkorrekturgleichung losen und p′ berechnen

4. Mit Gleichung (3.73) die Geschwindigkeiten, mit Gleichung (3.72) den Masse-

fluss und mit Gleichung (3.74) die Druckwerte korrigieren

5. Die Gleichungen der ubrigen Skalargroßen des Modells losen (turbulente kine-

tische Energie und Dissipation, Verteilung anderer Stoffe,...)

6. Wenn das Konvergenzkriterium nicht erfullt ist zuruck zu Schritt 2.

Fur instationare Stromungen muss ein Zeitschritt festgelegt werden, der von Iteration

zu Iteration variieren kann. Die Simulationen in dieser Arbeit werden allerdings mit

konstanten Zeitschritten durchgefuhrt.

Kapitel 4

Gebietszerlegungsverfahren

Dieses Kapitel dient der Einfuhrung mathematischer Grundprinzipien von Gebiets-

zerlegungsverfahren. Dies sind Losungsverfahren fur PDE die darauf basieren, das

Losungsgebiet zu zerlegen, die Losungen der Probleme auf den kleineren Teilgebie-

ten zu bestimmen, und die Teile zu einer Gesamtlosung zusammenzufugen.

Abbildung 4.1: Diskretisierter Losungsraum bestehend aus zwei Teilgebieten, aus

Mathew (2008)

Diese Technik, die auf Schwarz (1870) zuruckgeht, hat durch die Entwicklung von

Multiprozessorcomputern und Rechnerclustern in den letzten Jahrzehnten stark an

Bedeutung gewonnen. Mit wachsender Komplexitat und Große der Modelle wissen-

schaftlicher Anwendungen wird die Frage nach Parallelisierbarkeit der Losungsbe-

rechnung immer wichtiger, sie ist Kerngebiet aktueller Forschung. Die Entwicklungs-

arbeit erster moderner Gebietszerlegungsverfahren lasst sich anhand von Glowinski

(1988) nachvollziehen, Toselli und Widlund (2005) oder Quarteroni und Valli (1999)

geben einen umfassenden mathematischen Hintergrund der klassischen und heutigen

Methodik. Eine umfassende Behandlung aktueller Gebietszerlegungsverfahren findet

sich in Mathew (2008).

41

KAPITEL 4. GEBIETSZERLEGUNGSVERFAHREN 42

Ein haufiges Nebenprodukt der Entwicklung solcher Verfahren sind Vorkonditionie-

rer, welche die Losung großer linearer Gleichungssysteme vereinfachen. Hierauf soll

in dieser Arbeit jedoch nicht weiter eingegangen werden.

In den folgenden Abschnitten wird das bekannte Schwarz-Verfahren und die dar-

aus hervorgegangenen uberlappenden Gebietszerlegungsverfahren beschrieben, be-

schrankt auf eine Formulierung fur zwei Teilgebiete. Die Betrachtung von Gebiets-

zerlegungsverfahren ohne Uberlappung am Ende dieses Kapitels wird sehr kurz ge-

halten, da sie fur diese Arbeit keine große Rolle spielen.

4.1 Zerlegung in Teilprobleme

Definition 4.1. Zwei offene Teilgebiete Ωi ⊂ Ω fur i = 1, 2 nennt man nichtuber-

lappende Zerlegung von Ω , falls folgende Bedingung erfullt ist:

Ω1 ∪ Ω2 = Ω

Ω1 ∩ Ω2 = Ø

Den Rand der Teilgebiete bezeichnet man mit ∂Ωi , den inneren Teilrand mit B(i) ≡

∂Ωi ∩ Ω , den außeren Teilrand mit B[i] ≡ ∂Ωi ∩ ∂Ω und die Ubergangsstelle mit

B ≡ ∂Ω1 ∩ ∂Ω2 .

Abbildung 4.2: Losungsraum bestehend aus zwei uberlappenden Teilgebieten, aus

Mathew (2008)


Definition 4.2. Zwei offene Teilgebiete Ω∗

i ⊂ Ω fur i = 1, 2 nennt man uberlap-

pende Zerlegung von Ω , falls folgende Bedingung erfullt ist:

Ω∗

1 ∪ Ω∗

2 = Ω

Den Rand der Teilgebiete bezeichnet man mit ∂Ω∗

i , den inneren Teilrand mit B(i) ≡

∂Ω∗

i ∩ Ω und den außeren Teilrand mit B[i] ≡ ∂Ω∗

i ∩ ∂Ω , vgl. Abb. 4.2.

Die Zerlegung wird in der praktischen Anwendung durch die Geometrie von Ω be-

stimmt, oder aufgrund von Regularitatseigenschaften der Losung u gewahlt, falls

solche bekannt sind.

Die Theorie der Gebietszerlegung wird am Beispiel einer elliptischen Differentialglei-

chung zweiter Ordnung veranschaulicht:

Lu ≡ −div (a(x ) grad u) + b(x ) · grad u + c(x )u = f in Ω

u = 0 auf ∂Ω(4.1)

Hierbei sei Ω ⊂ Rd , und es gelte

0 < a0 ≤ a(x) , ∀x ∈ Ω

b(x ) und c(x) ≥ 0 seien glatt, und f(x) ∈ L2(Ω) .

Definition 4.3. Eine Zerlegung der Eins uber den uberlappenden Teilraumen Ω∗

1

und Ω∗

2 besteht aus glatten Funktionen χ1(x) und χ2(x) , die folgendes erfullen:

χi(x) ≥ 0 in Ω∗

i

χi(x) = 0 in Ω \ Ω∗

i

χ1(x) + χ2(x) = 1 in Ω

Fur eine nichtuberlappende Zerlegung Ω1 und Ω2 von Ω ergibt sich

χi(x) ≥ 0 in Ωi

χi(x) = 0 in Ω \ Ωi

χ1(x) + χ2(x) = 1 in Ω

mit Unstetigkeiten der χi am Ubergang B = ∂Ω1 ∩ ∂Ω2 .


Um die Losung u des Gesamtproblems 4.1 auf ein gekoppeltes System von Teillosun-

gen wi(x) auf Ωi (bzw. Ω∗

i ) zuruckfuhren zu konnen, mussen zwei Bedingungen

erfullt sein:

• Konsistenz : ui(x) , die Einschrankung von u auf ein Teilgebiet muss auch das

entsprechende Teilproblem losen, d.h. (u1(x), u2(x)) muss das gekoppelte Sys-

tem losen. Das sichert wi(x) = ui(x) fur i = 1, 2 , Konsistenz der Teillosungen

mit dem Gesamtproblem.

• Korrekte Problemstellung : Die Losung (w1(x), w2(x)) des gekoppelten Pro-

blems muss existieren, eindeutig sein und stetig von den Eingangsdaten abhangen.

Stabilitat und Eindeutigkeit sind wichtig fur eine numerische Approximation

der Losung.

Dann lasst sich die Losung u des Gesamtproblems mit einer geeigneten Zerlegung

der Eins wie folgt ausdrucken:

u(x) = χ1(x)w1(x) + χ2(x)w2(x) . (4.2)

Anstatt die Losung des Gesamtproblems zu bestimmen kann man zwei kleinere Teil-

probleme mit zusatzlichen Nebenbedingungen losen, die folgende Form annehmen:

Lwi = fi in Ωi (bzw. Ω∗

i )

Ti(wi, γ) = gi auf B(i)

wi = 0 auf B[i]

fur i = 1, 2 (4.3)

Ti(wi, γ) bezeichnet einen Randbedingungsoperator, der entsprechende Dirichlet-,

Neumann- oder Robin- Randbedingungen setzt um die Losungen zu koppeln und

Konsistenz und Stabilitat zu wahren. Je nach Problemformulierung nimmt er eine

andere Form an.

Auch die Nebenbedingungen, die die Teillosungen stetig koppeln, konnen unter-

schiedliche Formen annehmen:

Algebraische Gleichungen:

ui − uj = 0 auf ∂Ωi ∩ ∂Ωj im nicht uberlappenden,

ui − uj = 0 auf ∂Ω∗

i ∩ ∂Ω∗

j im uberlappenden Fall,


oder Nebenbedingungen an die Flusse zwischen den Gebietsteilen, formuliert uber

Differentialgleichungen:

ni · (a(x ) grad ui) + nj · (a(x ) grad uj ) = 0 auf ∂Ωi ∩ ∂Ωj

im nicht uberlappenden, sowie

ni · (a(x ) grad ui) − ni · (a(x ) grad uj ) = 0 auf ∂Ω∗

i ∩ ∂Ω∗

j

im uberlappenden Fall.

ni bezeichnet hier den Normalen-Einheitsvektor auf ∂Ωi .

Mathew (2008) fuhrt eine Bedingung ein, um die korrekte Problemstellung der ge-

koppelten Formulierung zu sichern und somit die Losungsberechnung mit stabilen

numerischen Verfahren zu ermoglichen. Fur eine Konstante C > 0 unabhangig von

den Eingangsdaten gelte die Schranke

(‖w1‖ + ‖w2‖) ≤ C(‖|f1‖| + ‖|f2‖| + ‖|g1‖| + ‖|g2‖|) (4.4)

mit ‖ · ‖ und ‖| · ‖| geeigneten Normen fur Losung und Eingangsdaten.

4.2 Uberlappende Verfahren

Die Klasse der uberlappenden Gebietszerlegungsverfahren geht auf die Arbeit des

deutschen Mathematikers H. A. Schwarz zuruck. In Schwarz (1870) entwickelte er

ein iteratives Losungsverfahren fur die Laplace-Gleichung auf einem unregelmaßigen

Losungsgebiet, alternierendes Schwarz-Verfahren genannt. Diese Methode lasst sich

auf eine große Klasse elliptischer Differentialgleichungen erweitern und fuhrte zur

Entwicklung vieler Divide and Conquer-Verfahren.

4.2.1 Das alternierende Schwarz-Verfahren

Motivation

Sei u(x) die Losung des Ausgangsproblems (4.1). Definiere wi(x) := u(x) auf Ω∗

i

fur 1 ≤ i ≤ 2 . Dann ist Lwi = f auf Ω∗

i und w1(x) = w2(x) auf Ω∗

1 ∩ Ω∗

2 , und es

gilt:

Lw1 = f in Ω∗

1

w1 = w2 auf B(1)

w1 = 0 auf B[1]

und

Lw2 = f in Ω∗

2

w2 = w1 auf B(2)

w2 = 0 auf B[2]

(4.5)


Falls dieses System fur w1(x) und w2(x) korrekt gestellt ist, d.h. eine eindeutige

Losung existiert, die stetig von ihren Startwerten abhangt, so kann durch u(x) durch

die Losung der beiden Teilprobleme bestimmt werden.

Es gilt folgender Eindeutigkeitssatz:

Theorem 4.1. Sei c(x) ≥ 0 und div b(x ) ≤ 0 . Sei u(x) eine genugend glatte

Losung von (4.1), w1(x) und w2(x) genugend glatte Losungen des Systems (4.5).

Dann gilt:

u(x) =

w1(x) in Ω∗

1

w2(x) in Ω∗

2

(4.6)

Beweis: Maximumprinzip, siehe Mathew (2008), Kapitel 1.2.

Zur iterativen Bestimmung einer Losung u geht man wie folgt vor:

Algorithmus

1. Sei v(0) die Startapproximation der Losung.

2. Bestimme w(k+1)1 als Losung von

Lw(k+1)1 = f1 in Ω∗

1

w(k+1)1 = v(k) auf B(1)

w(k+1)1 = g auf B[1]

(4.7)

3. Aktualisiere die Losungapproximation v(k+1/2) :

v(k+1/2) ≡

w(k+1)1 in Ω∗

1

v(k) in Ω \ Ω∗

1

(4.8)

4. Bestimme w(k+1)2 als Losung von

Lw(k+1)2 = f2 in Ω∗

2

w(k+1)2 = v(k+1/2) auf B(2)

w(k+1)2 = g auf B[2]

(4.9)

5. Aktualisiere die Losungsapproximation v(k+1) :

v(k+1) ≡

w(k+1)1 in Ω∗

2

v(k+1/2) in Ω \ Ω∗

2

(4.10)



Bedingungen fur die Stabilitat von Problem (4.5) ergeben sich aus aus der Betrach-

tung des Problems mit gestorten Eingangsdaten:

Lw1 = f1 in Ω∗

1

w1 = w2 + r1 auf B(1)

w1 = 0 auf B[1]

und

Lw2 = f2 in Ω∗

2

w2 = w1 + r2 auf B(2)

w2 = 0 auf B[2]

(4.11)

Ist (4.11) eindeutig losbar und gilt eine Bedingung der Form

(‖|w1‖| + ‖|w2‖|) ≤ C(‖f1‖ + ‖f2‖ + ‖r1‖ + ‖r2‖) (4.12)

fur geeignete Normen, dann ist (4.5) korrekt gestellt, vgl. Mathew (2008), Kapitel 15

fur Maximum-Norm-Konvergenz.

Fur elliptische Differentialgleichungen der Form (4.1) lasst sich allgemeine Konver-

genz der Losungsapproximation v(k) gegen u zeigen, vgl. Quarteroni und Valli

(1999), Kapitel 4.6 und 5 oder Mathew (2008), Kapitel 15. Unter gewissen Ein-

schrankungen (Große der Uberlappung, Durchmesser der Teillosungsgebiete) ist die

Konvergenz geometrisch.

Formale konvergenztheoretische Resultate fur Schwarz-Verfahren auf Basis einer FV-

Diskretisierung sind dem Autor nicht bekannt. Einige einfache Konvergenzergebnisse

lassen sich anhand einer FD-Diskretisierungen mit dem Maximumsprinzip beweisen.

Die umfassendsten Ergebnisse bezuglich der Konvergenz von Schwarz-Verfahren ent-

stammen den funktionalanalytischen Mitteln der FE-Methoden. Die hierfur benotig-

te Variationsformulierung des Problems soll der Vollstandigkeit halber genannt wer-

den, obwohl im Rahmen dieser Arbeit nicht weiter auf die Theorie eingegangen wird:

Variationsformulierung der Schwarz-Methode

Die Formulierung der Schwarz-Methode als Problem der Variarionsrechnung geht

zuruck auf P. L. Lions, siehe Glowinski (1988), S.1-42. Hier wird sie am Beispiel der

elliptischen Differentialgleichung (3.20) eingefuhrt:

Sei V := H10 (Ω), V 0

i := H10 (Ω∗

i ) fur i = 1, 2 und

V ∗

i := v ∈ V : v = 0 in Ω \ Ω∗

i (4.13)


Sei u0 ∈ V eine Startapproximation der Losung. Bestimme fur alle k ≥ 0

wk1 ∈ V 0

1 : A1(wk1 , v1) = F1(v1) −A1(u

k, v1) ∀v1 ∈ V 01

uk+1/2 = uk + wk1

wk2 ∈ V 0

2 : A2(wk2 , v2) = F2(v2) −A2(u

k+1/2, v2) ∀v1 ∈ V 02

uk+1 = uk+1/2 + wk2

(4.14)

wki bezeichnet die Fortsetzung von wk

i durch 0 auf Ω \ Ω∗

i ,

Ai und Fi die entsprechenden Einschrankungen A|Ω∗

iund F |Ω∗

i.

4.2.2 Additives Schwarz-Verfahren

Das alternierende Schwarz-Verfahren erfordert eine sequentielle Berechnung der Teillosun-

gen. Eine fur Parallelisierung besser geeignete Variante des Algorithmus, genannt

additives Schwarz-Verfahren ergibt sich wie folgt:

1. Seien w(0)1 , w

(0)2 die Startapproximation der lokalen Losungen von (4.1).

2. Bestimme fur i = 1, 2 parallel w(k+1)i als Losung von

Lw(k+1)i = f in Ω∗

i

w(k+1)i = χ1(x)w

(k)1 (x) + χ2(x)w

(k)2 (x) auf B(i)

w(k+1)1 = 0 auf B[i]

(4.15)


Falls c(x) ≥ c0 > 0 und die Gebiete genugend weit uberlappen, konvergieren die

Approximationen v(k) , definiert durch

v(k) ≡ χ1(x)w(k)1 (c) + χ2(x)w

(k)2 (x) (4.16)

geometrisch gegen die Losung u , siehe Mathew (2008), Kapitel 15.

4.3 Gebietszerlegung ohne Uberlappung

Neben der Schwarz-Methode gibt es eine andere große Klasse von Verfahren zur Auf-

teilung des Losungsgebietes: Nicht uberlappende Verfahren. Sie verwenden Formen

der Teillosungskoppelung, die ohne gemeinsamen Uberlappungsbereich auskommen.

Der Vollstandigkeit halber werden hier die wichtigsten Vertreter genannt, obwohl sie

fur den weiteren Verlauf dieser Arbeit keine Rolle spielen.


• Steklov-Poincare-Verfahren, auch Schur-Komplement-Methode: Hier wird ne-

ben der Forderung einer stetigen Gesamtlosung auch ein stetiger Fluss zwischen

den Gebieten als Nebenbedingung der Teilprobleme direkt umgesetzt: Fur Glei-

chung (4.1) fuhrt man folgende Nebenbedinungen am Ubergang

B = ∂Ω1 ∩ ∂Ω2 der Teilgebiete ein:

w1 = w2 auf B

n1 · (a gradw1 − bw1 ) = n1 · (a gradw2 − bw2 ) auf B(4.17)

ni ist der Normalenvektor auf ∂Ωi .

Zusatzlich wird der Steklov-Poincare-Operator eingefuhrt:

S(g, f1, f2) ≡ n1 · (a gradw1 − bw1 ) − n1 · (a gradw2 − bw2 ) (4.18)

Hierbei ist g(·) ein Satz glatter Dirichlet-Randbedinungen, wi fur i = 1, 2 sind

Losungen von

Lwi = fi in Ωi

wi = 0 auf B[i]

wi = g auf B

(4.19)

Findet sich eine Funktion g(·) , fur die

S(g, f1, f2) = 0 (4.20)

gilt, so erfullen die lokalen Losungen wi die Koppelungs-Nebenbedingungen

(4.17). Die Suche einer Losung u von (4.1) lasst sich auf diese Weise auf die

Losung des Steklov-Poincare-Problems (4.20) reduzieren. Ist diese Losung ge-

funden, dient sie als Randbedingung fur die Probleme auf den Teilgebieten, die

dann parallel und unabhangig voneinander gelost werden konnen.

Bei eine diskrete Systemmatrix entspricht diese Reduktion des Systems einer

blockweisen Gauss-Elimination der Unbekannten im Inneren der Teilgebiete.

Fur die iterativen Losung dieses Schur-Komplement-Systems sei auf Mathew

(2008), Kapitel 3 verwiesen. Dort werden auch gebrauchliche Vorkonditionierer

wie Neumann-Neumann oder BDD (Balancing Domain Decomposition) behan-

delt.

• Lagrange-Verfahren: Eine Formulierung als Optimierungsproblem mit Lagran-

ge-Nebenbedingungen zur Koppelung der Teillosungen kann bei Problemen


zum Einsatz kommen, die das Optimalitatsprinzip erfullen. Dies ist fur eine

Gleichung vom Typ (4.1) nicht unbedingt gegeben, eine Einschrankung auf

b(x ) = 0 und c(x ) ≥ 0 ist erforderlich. Solche Probleme nennt man selbst-

adjungiert und koerziv. Eine Lagrange-Sattelpunktformulierung eines solchen

Problems siehe (3.27).

Verfahren, die Lagrange-Optimierung verwenden sind beispielsweise die

FETI (Finite Element Tearing and Interconnection)-Methode der parallelen,

diskreten Optimierung unter Nebenbedingungen oder die Mortar-Methode zur

Losung von elliptischen Differentialgleichungen auf nicht uberlappenden Git-

tern bei unterschiedlicher Diskretisierung der Teilraume. Eine Beschreibung

dieser Methoden geht uber den Rahmen dieser Arbeit hinaus, man findet sie

in Toselli und Widlund (2005) oder Mathew (2008). Dort finden sich auch

Mischformen, die Teilaspekte von Schur-Komplementmethoden und Lagrange-

Verfahren vereinen. Genannt seien FETI-DP (Finite Element Tearing and In-

terconnection Dual-Primal) oder BDDC (Balancing Domain Decomposition by

Constraints).

• Kontrollmethode der kleinsten Quadrate: Dieser Ansatz aus der Kontrolltheorie

bestimmt auf jedem Teilgebiet eine unbekannte Funktion, welche die PDE unter

Berucksichtigung unbekannter Randbedinungen, den so genannten Kontrollda-

ten lost. Diese Kontrolldaten mussen so gewahlt werden, dass die Teillosungen

zu einer globalen Losung gekoppelt werden. Anwendbar ist dieses Prinzip auf

nicht selbstadjungierte, elliptische Differentialgleichungen vom Typ (4.1) so-

wohl auf uberlappenden als auch nicht uberlappenden Gebietszerlegungen.

Die mathematischen Formulierung ist ein Minimierungsproblem der Differenz

der unbekannten Funktionen auf dem gemeinsamen Bereich bezuglich einer

quadratischen Norm unter der Nebenbedingung, dass die Unbekannten die el-

liptische Gleichung auf ihrem Teilgebiet losen. Neben einer Behandlung als

Sattelpunktproblem ist auch eine Formulierung als unbedingtes Minimierungs-

problem moglich, vgl. Mathew (2008), Kapitel 6.

Kapitel 5

Das Overlapping Grid Verfahren

fur Finite Volumen mit bewegten

Gittern

In diesem Kapitel wird die Methode uberlappender Gitter eingefuhrt, die im weiteren

Verlauf dieser Arbeit bezuglich ihrer Anwendbarkeit auf komplexe Stromungsmecha-

nikprobleme untersucht werden soll. Nach einer allgemeine Beschreibung wird die

mathematische Umsetzung des Overlapping Grid Verfahrens (OLG) vorgestellt, das

Hadzic (2005) fur das in Kapitel 3 eingefuhrte Finite-Volumen-Verfahren entwickelt

hat. Zum Abschluss folgt die Modellierung von Gitterbewegungen sowie ein kurzer

Uberblick aktueller Implementierungen des OLG-Verfahrens.

5.1 Einleitung

Die Entwicklung moderner Verfahren zur Losung von Differentialgleichungen mit

uberlappenden Gittern (auch engl. Overlapping Grid oder Chimera-Methoden), geht

auf verschiedene Autoren zuruck. Zu nennen sind die Arbeiten von Atta (1981),

Henshaw (1985), Chesshire und Henshaw (1990), Dougherty (1985) und Benek u. a.

(1986), in denen Verfahren fur unterschiedliche Anwendungsbereiche entwickelt wer-

den.

Das Funktionsprinzip ist identisch: Das diskrete Losungsgebiet besteht aus einer

Anzahl unabhangiger Teilgitter, die sich beliebig uberlappen. Die Gesamtlosung wird

durch Koppelung dieser Teillosungen gewonnen. Diese getrennte Vernetzung bringt

51

KAPITEL 5. DAS OVERLAPPING GRID VERFAHREN 52

folgende Vorteile, die am Beispiel der spater durchgefuhrten Rechnungen erlautert

werden:

• Optimale Teilvernetzung : Die Teilkomponenten des Systems konnen getrennt

voneinander optimal vernetzt werden. So kann ein Hintergrundgitter aus nume-

risch vorteilhaften, regelmaßigen Hexaedern verwendet werden, wahrend man

die Flugelprofile mit einem an die Geometrie angepassten Netz versieht, das

wichtige Details der Stromung auflosen kann. Das in der Praxis große Pro-

blem der Vernetzung eines regelmaßigen Ubergangs zwischen den beiden Teilen

entfallt, die Koppelung wird vom Algorithmus ubernommen.

• Flexible Gitterbewegungen: Ein großer Vorteil des OLG-Verfahrens ist die Moglich-

keit, Bewegungen ohne Veranderung des Gitters durchfuhren zu konnen. Bei

herkommlichen Verfahren gibt es folgende Moglichkeiten: Man kann durch De-

formation von Zellen, also eine Anpassung des Netzes, Bewegungen erzeugen.

Da die Netzstruktur erhalten bleiben muss um eine stabile Losung zu liefern,

lassen sich auf diese Weise nur kleine Bewegungen realisieren. Großraumige

Bewegungen erfordern eine aufwendige Neugenerierung des Netzes zwischen

den Zeitschritten (Re-Meshing), oder man verwendet die Technik der Sliding

Interfaces, wie sie der kommerzielle CFD-Loser Comet von CD-adapco anbie-

tet: Zwischen stationaren und bewegten Teilen des Netzes wird eine Sliding-

Flache definiert, an der in jedem Zeitschritt die aktuellen Nachbarzellen fur

die Losungsberechnung bestimmt werden. So konnen zumindest regelmaßige

Bewegungen auf vordefinierten Bahnverlaufen realisiert werden. Beispielswei-

se simuliert man die Rotation des VSP uber die Drehung eines zylindrischen

Radkorpers, der sich in einer entsprechenden Aussparung des Hintergrundgit-

ters befindet, wahrend die Schwenkbewegungen der Flugel durch eine zylin-

drische Form der Flugelnetze erfolgt, die sich in Aussparungen im Radkorper

bewegen, siehe Abb. A.11. Mit dieser Technik lassen sich zwar relativ komple-

xe Bewegungen umsetzen, allerdings macht die Form der verwendeten Gitter

die Erweiterung um Strukturen wie Streben unmoglich. Details finden sich in

Kapitel 6. Der Overlapping Grid Algorithmus ermoglicht Bewegungen ohne

Strukturen wie den Radkorper, was sowohl die Gittergenerierung vereinfacht,

als auch Streben in Flugelnahe ermoglicht.

• Interaktion mehrerer Korper : Des Weiteren ermoglicht der Overlapping Grid

Ansatz die beliebige Bewegung mehrerer Objekte relativ zueinander, eine mit

herkommlichen Methoden sehr komplexe bis nicht realisierbare Aufgabe.


Abbildung 5.1: OLG-System aus Hintergrundgitter und zwei Overlap-Objekten, aus

Hadzic (2005) mit Ubersetzungen des Autors

5.2 Overlapping Grid Methodik

Zur Losung eines Stromungsproblems auf einem System von uberlappenden Gittern

sind folgende Schritte notwendig:

1. Definition der Teilgebiete: Ein System uberlappender Gitter besteht in der Re-

gel aus einem großeren, Hintergrundgitter genannten Teil, und kleineren Git-

tern, die Objekte oder Bereiche einschließen, die fur die Rechnung von beson-

derer Bedeutung sind. Sie werden im Folgenden mit Overlap-Gitter bezeichnet,

vgl. Abb. 5.1. Die Grenzen dieser Gebiete sind einerseits durch die Geometrie

der umstromten Korper gegeben. Diese werden durch die ublichen, im Ab-

schnitt 3.5.10 definierten Randbedingungen modelliert. Andererseits werden die

Außenrander der Overlap-Gebiete durch einen neuen Typ Rand definiert, dem

Overlap-Rand, der keine physikalischen Eigenschaften besitzt, sondern dem Al-

gorithmus als Ausgangspunkt fur die Koppelung der Teilgebiete dient.

2. Bestimmung des aktiven Losungsgebietes : Ausgehend vom Rand der Overlap-

Gebiete muss das Gebiet bestimmt werden, das zur Berechnung der Losung

beitragt. So werden die Zellen des Hintergrundgitters, die vom Overlap-Gitter

uberdeckt werden, großenteils von der Losungsberechnung ausgeklammert. Der

Teil der uberdeckten Zellen, der in der Gegend des Overlap-Randes liegt, bleibt

fur die Koppelung der Teilgebiete aktiv.

3. Teilbereich-Randbedingungen setzen: Die Teilbereiche werden durch Interpola-

tion der Stromungswerte zwischen den Gittern gekoppelt: Die Randbedingun-

gen der Overlapgitter ergeben sich aus Zellwerten des Hintergrundgitters in der


Nahe des Overlap-Randes. Genauso werden Werte des Overlapgitters fur die

Randwerte der aktiven Teile des Hintergrundgitters verwendet.

4. Losung der Teilprobleme: Die Teilprobleme sind so mit einem kompletten Satz

Anfangs- und Randbedingungen versehen und konnen theoretisch unabhangig

voneinander gelost werden. Untersuchungen in Abschnitt 5.3.2 zeigen jedoch,

dass eine gemeinsame Losung vorteilhafter ist.

5. Koppelung der Teillosungen: Zusatzlich werden einige Schritte durchgefuhrt,

um die Koppelung der Teillosungen zu verstarken. Interpolation zwischen den

Gittern alleine reicht nicht aus, da z.B. die Massenerhaltung verletzt ist.

Abbildung 5.2: Uberlappungsgebiet mit Bezeichnungen, aus Hadzic (2005) mit Uber-

setzungen des Autors

5.2.1 Zellfunktionen

Eine KV-Zelle nimmt fur die Losungsberechnung eine der folgenden Funktionen an,

vgl. Abb. 5.2:

Aktive Zellen

Die Zelle ist normaler Bestandteil des Losungsgebiets, entsprechend einer nicht uber-

lappenden Rechnung.

Inaktive Zellen

Der Großteil der Zellen, die durch ein anderes Gitter uberdeckt werden, spielt fur die

Losungsberechnung keine Rolle und kann ausgeklammert werden. Bei Rechnungen


ohne Gitterbewegungen konnen sie geloscht werden, fur die Modellierung bewegter

Gitter, bei der sich die uberlappenden Bereiche verandern, ist eine zeitweise Deak-

tivierung solcher Zellen geschickter. Dieser Ansatz wird auch fur unbewegte Gitter

beibehalten, da die betroffenen Bereiche in der Regel nur einen kleinen Teil des ge-

samten Losungsgebietes ausmachen. Ein Algorithmus bestimmt automatisch die zu

deaktivierenden Zellen. Ein gewisser Randbereich des uberdeckten Gebiets bleibt fur

die Koppelung der Teillosungen aktiv.

Interpolationszellen

Die Zellen am Rand des Overlap-Gebietes spielen eine besondere Rolle, hier werden

die Teilgebiete durch Interpolation zwischen den Gittern gekoppelt. Die Zellschicht

des Overlap-Gitters, fur die eine Overlap-Randbedingung gesetzt ist, wird zu einer

Schicht von Interpolationszellen Pi . Diese beziehen ihre Werte uber ein Interpolati-

onsschema aus Zellen Sk des Hintergrundgitters:

φPi=

NS∑

k=1

αkφSk(5.1)

Die Interpolationsgewichte αk sowie die Anzahl der an der Interpolation beteiligten

Zellen NS sind vom verwendeten Verfahren abhangig. Am Ubergang von aktiven

zu inaktiven Zellen des Hintergrundgitters wird eine identische Koppelung durch-

gefuhrt: Die erste Schicht aktiver Zellen interpoliert ihre Werte vom Overlap-Gitter.

Die an der Interpolation beteiligten Zellen des anderen Gitters heißen Spenderzellen

(engl. donor cells). Sie werden uber einen Suchalgorithmus ausgewahlt, der wie das

verwendete Interpolationsschema in Kurze beschrieben wird.

Ubergangszellen

Fur bewegte Gitter kann es zu Situationen kommen, in denen es Sinn macht, die

Werte gewisser Zellen zu berechnen, obwohl sie nicht fur die Losung des aktuellen

Zeitschrittes benotigt werden. Diese zunachst inaktiven Zellen, fur die man aktuelle

Stromungswerte bestimmt, die fur spatere Zeitschritte wichtig sein konnen, nennt

man Ubergangszellen. In der Regel kann man davon ausgehen, dass diese Bereiche

eher klein sind und der Rechenaufwand fur diese zusatzlichen Zellen daher nicht von

Bedeutung ist.

Die Funktion einer Zelle kann sich - insbesondere bei bewegten Gittern - im Verlauf

der Rechnung andern.


5.2.2 Bestimmung inaktiver Zellen

Ein einfacher Algorithmus zur Bestimmung der Zellen, die deaktiviert werden konnen,

funktioniert wie folgt:

Ausgehend von einem Overlap-Rand bleiben alle Zellen außerhalb dieses Bereiches

aktiv. Die Zellen im Inneren werden uberpruft, ob sie im Overlap-Bereich liegen,

dem Bereich, der benotigt wird, um eine gute Koppelung der Losungen zu erreichen.

Die Breite δ0 dieses Bereiches ist abhangig von der verwendeten Diskretisierung, er

entspricht in der Praxis bis zu funf kompletten Zellschichten. Zellen, deren Abstand

zum Overlap-Rand großer als δ0 ist, werden deaktiviert. Hierzu bildet man das

Skalarprodukt aus n , dem nach außen gerichteten Normalenvektor in dem Punkt

P0 des Overlap-Randes, der dem untersuchten Gitterpunkt P am nachsten ist, und

dem Vektor r , der P0 und P verbindet. Es ist negativ, falls P im Inneren liegt.

Falls dann δ = |n · r| > δ0 , so wird P deaktiviert.

Bereiche, die sich außerhalb des Losungsgebietes befinden, werden ebenso deaktiviert.

Hierbei kommen nicht nur Gebietsrander zum Tragen, sondern auch Rander von

Korpern im Overlap-Gitter, vgl. Abbildung 5.3. Die Uberlappung mehrerer Overlap-

Korper wird bei Hadzic (2005) uber eine Hierarchie der Overlap-Gitter geregelt,

die man bei der Generierung angibt. Overlapping-Gitter, die sich bewegen, werden

gesondert behandelt. Da die korrekte Bestimmung einer Hierarchie hier erschwert

ist, werden Zellen dieser Korpergitter nicht deaktiviert, sondern als Ubergangszellen

mitberechnet.

Abbildung 5.3: Aktive Teilbereiche bei Mehrfachuberlappungen, aus Hadzic (2005)

mit Ubersetzungen des Autors


”... if a valid composite grid is created for some ordering of component grids then a

valid composite grid should result for all orderings of component grids (although we

have not attempted to prove this)”, so Chesshire und Henshaw (1990), S.16 uber ih-

ren Koppelungsalgorithmus. Sie halten die Anordnung der Teilgitter zueinander fur

nebensachlich, wenn die Struktur der verwendeten Gitter uberhaupt eine sinnvolle

Koppelung zu einem Gesamtnetz zulasst. Lohner u. a. (2001) schalten bei Mehr-

fachuberlappungen diejenigen Zellen mit dem kleinsten Volumen aktiv, da sie die

Stromungseigenschaften am besten auflosen konnen.

5.2.3 Spenderzellensuche

Hadzic (2005) bestimmt die Spenderzellen uber den neighbor-to-neighbor -Suchalgorithmus

von Lohner (1995). Zu einer Zelle P wird die Zelle S des anderen Gitters gesucht,

die P, den Schwerpunkt von P enthalt. Von einer Startzelle A wird iterativ die

Nachbarzelle bestimmt, die in Richtung von P liegt, bis S erreicht ist. Die ubrigen

Spenderzellen des gewahlten Interpolationsschemas sind aus der direkten Umgebung

von S und erfordern keine komplexe Suche.

Abbildung 5.4: Neighbor-to-neighbor Suchverfahren, aus Hadzic (2005) mit Uberset-

zungen des Autors

Es wird uber alle Seitenflachen der Startzelle das Skalarprodukt aus ni , dem Nor-

malenvektor der Seite und pi gebildet, dem Vektor der den Schwerpunkt der Seite

mit P verbindet, vgl. Abb. 5.4 b). Der kleinste Winkel zwischen beiden Vektoren

bestimmt die Startzelle des nachsten Schrittes. Sobald die Skalarprodukte uber alle


Seiten negativ werden ist die Spenderzelle erreicht. Nicht konvexe Zellformen, wie sie

in unstrukturierten Gittern auch auftreten konnen, benotigen ein anderes Abbruch-

kriterium.

Im Verlauf des Verfahrens kann es abhangig von der Wahl des Startpunktes am

Rand des Gebietes zu Problemen kommen, wie aus Abb. 5.4 a) hervorgeht: Startet

der Algorithmus in B , so stoppt er in C , einer Zelle, die am Gebietsrand liegt

und daher keine Nachbarzelle in Suchrichtung besitzt. Lohner (1995) nennt folgende

Methoden, dieses Problem zu beheben:

• Einen neuen Startpunkt aus den Randzellen wahlen, falls diese bekannt sind.

Hadzic (2005) bestimmt diejenige Randzelle, welche am nachsten an P liegt,

als neue Startzelle. Vgl. Abb. 5.4 a) Zelle D .

• Zuerst das Innere der Gebiete interpolieren und die Punkte in der Nahe der

Rander spater behandeln.

• Den Algorithmus in Nachbarn der Startzelle neu starten.

• Brute Force, falls alle andern Methoden gescheitert sind.

Fur bewegte Gitter treten diese Probleme besonders im ersten Zeitschritt auf, spater

liegen bessere Startwerte fur die Suche vor, z.B. die Position der Spenderzelle im

vorherigen Zeitschritt.

5.2.4 Interpolationsschema

Das verwendete Diskretisierungsverfahren bestimmt die Genauigkeit der Teillosun-

gen. Beim Zusammenfugen der Teile zu einer Gesamtlosung sollte diese Genauigkeit

erhalten bleiben. Dazu tragt nicht nur das gewahlte Interpolationsverfahren, sondern

auch dessen Ordnung bei.

Verfahren hoherer Ordnung benotigen einen großeren Overlap-Bereich, da mehr Spen-

derzellen fur die Interpolation verwendet werden. Aber auch bei Verfahren niederer

Ordnung kann ein Overlap, der großer ist als das Verfahren vorgibt, einen positiven

Effekt haben.

Untersuchungen bezuglich der Interpolation zwischen Gittern kommen zu folgenden

Ergebnissen:


• Henshaw (1985) betrachtet ein elliptisches Problem auf einem eindimensiona-

len, uberlappenden Losungsgebiet. Die Teillosungen sind mit einem Verfahren

zweiter Ordnung diskretisiert. Es wird gezeigt, dass eine Gesamtlosung mit

dieser Konvergenzordnung durch ein Interpolationsverfahren dritter Ordnung

erreicht werden kann. Hierbei gilt δ = dh , mit jeder Verfeinerung der Git-

terweite h verkleinert sich auch der Uberlappungsbereich δ . Gilt δ ≡ d , die

Uberlappung ist großer als eine Konstante, so genugt lineare Interpolation, ein

Verfahren zweiter Ordnung.

• Chesshire und Henshaw (1990) beweisen eine Erweiterung dieser Aussage auf

eindimensionale PDEs hoherer Ordnung und zeigen empirisch am Beispiel der

Poisson-Gleichung Ubertragbarkeit der Ergebnisse ins Zweidimensionale.

• Mastin und McConnaughey (1984) untersuchen fur lineare und bilineare Inter-

polationsverfahren eine Vergroßerung des Overlap-Bereiches von einem auf zwei

Zellschichten. Sie beobachten eine deutliche Verbesserung der Konvergenzge-

schwindigkeit, wenn die Zellgroßen beider Gitter weitgehend ubereinstimmen.

Daruber hinaus ist fur ein dort untersuchtes Laplace-Problem das bilineare

Verfahren einer Interpolation durch Taylorentwicklung uberlegen.

• Benek u. a. (1986) verwenden trilineare Interpolation fur ihr dreidimensionales

Chimera-Verfahren. Sie schlagen einen Overlap von vier bis funf Zellen vor. Es

treten Schwierigkeiten beim Ubergang von Schocks durch die Overlap-Grenzen

auf. Deren Ursache kann nicht eindeutig geklart werden, sie nennen die oben

genannten kritische Parameter (Interpolationsverfahren, Große des Overlapge-

bietes, relative Zellgroßen der Teilgitter).

Hadzic (2005) verwendet lineare Interpolation um die Teilgitter zu koppeln (vgl. li-

neare FEM-Basisfunktionen, Kapitel 3.4.3). Er zeigt anhand einiger Beispiele, dass

sich fur das betrachtete Verfahren Konvergenz zweiter Ordnung einstellt. Dies ent-

spricht den Erwartungen fur das hier verwendete Diskretisierungsverfahren der Ord-

nung zwei.

Allerdings stellt sich folgendes Problem: Durch die Interpolation sind an der Overlap-

Grenze die Erhaltungsgleichungen in der Regel verletzt. Die Massenerhaltung ist

fur die Druckkorrektur im vorgestellten Verfahren aber von großer Bedeutung. Eine

Losung wird in Abschnitt 5.3.3 vorgestellt.


5.3 Losungsverfahren

In diesem Abschnitt wird die Losung der Navier-Stokes-Gleichungen auf einem Sys-

tem uberlappender Gebiete behandelt, die durch Interpolation gekoppelt sind. Dazu

wird das in Kapitel 3 entwickelte Finite Volumen Verfahren wie folgt angepasst:

5.3.1 Behandlung inaktiver Zellen

In den diskreten Systemgleichungen (3.63) werden die Koeffizienten von inaktiven

Zellen wie folgt gesetzt:

aP0Φ := 1

aj := 0

bP0:= ΦP0

(5.2)

Somit andern sich die Variablenwerte in einer Zelle nicht, solange sie inaktiv ist.

Die Abbruchbedingung der außeren Iterationen (3.89) muss entsprechend angepasst

werden, dass die inaktiven Zellen keinen Einfluss auf das Maß der Konvergenz haben:

Dies wird durch ein Array ib realisiert, in dem der Status der Zellen aller Gitter

gespeichert ist:

ib =

1 fur Spenderzellen

0 fur inaktive und Interpolationszellen(5.3)

Das modifizierte Abbruchkriterium ergibt sich wie folgt:

Rkφ =

N∑

i=1

|bP0− aP0φφP0

ib +Nj∑

j=1

ajφPjib|

Nj∑

j=1

|aP0φφP0ib|

(5.4)

5.3.2 Starke Koppelung der Teillosungen

Ein Losungsalgorithmus fur ein System uberlappender Gitter, der auf klassischen

Schwarz-Verfahren basiert wie sie in Kapitel 4 eingefuhrt wurden, kann in seiner

Konvergenzgeschwindigkeit signifikant verbessert werden. Im Folgenden wird eine

Modifikationen des Verfahrens vorgestellt, die zu einer verbesserten, starken Koppe-

lung der Teillosungen fuhrt.


Abwandlung des alternierenden Schwarz-Algorithmus

Erste Overlapping Grid Verfahren zur numerischen Losung von PDEs beruhen auf

dem alternierenden Schwarz-Gebietszerlegungsverfahren, siehe Starius (1977). Dort

wird Konvergenz des Verfahrens fur eine generelle Klasse elliptischer Probleme auf

einem uberlappenden Losungsgebiet gezeigt. Da die auftretenden Gleichungssyste-

me aufgrund ihrer Große ohnehin iterativ behandelt werden mussen, fuhrt folgende

Modifikation des Verfahrens (4.2.1) zu einer deutlichen Verbesserung der Konver-

genzgeschwindigkeit:

Der klassische Algorithmus lost ein Teilproblem erst exakt und aktualisiert dann

die Randwerte der ubrigen Probleme. Außere LGS-Losungsiterationen verwenden

somit haufig veraltete Randwerte. Man spricht von einer schwachen Koppelung der

Teillosungen. Eine Alternative besteht darin, nach jeder außeren Iteration die Rand-

werte der ubrigen Teilgebiete zu aktualisieren. Diese beiden Schritte - eine außere

Iteration und Randwertinterpolation - werden nacheinander auf allen Teilgebieten

durchgefuhrt. Im Gegensatz zur ursprunglichen Form des Verfahrens verwendet man

auf diese Weise immer die aktuellsten Knotenwerte zur Losung der Teilprobleme.

Simultanes Losungsverfahren

Henshaw (1985) realisiert dies, indem er die diskretisierten Systemgleichungen al-

ler Teilprobleme und deren Interpolationsgleichungen fur die Randwerte zu einem

gemeinsamen Gleichungssystem zusammenfugt.

Auch Hadzic (2005) verwendet ein solches stark gekoppeltes, simultanes Losungs-

verfahren. Das linearisierte Gleichungssystem (3.87) jedes Teilproblems wird in ein

gemeinsames System geschrieben.

Interpolationszellen PI sollen ihre Werte nicht von ihren Zellnachbarn sondern von

den Spenderzellen der anderen Teilprobleme beziehen. Dazu wird das gemeinsame

Gleichungssystem wie folgt modifiziert: Die aus der Diskretisierung des Teilproblems

stammende algebraische Gleichung fur PI (3.63) wird durch die entsprechende In-

terpolationsgleichung (5.1) ersetzt:

aPIφφPI−

Nj∑

j=1

ajφPj= bPI

⇒ φPI−

NS∑

k=1

αkφSk= 0 (5.5)


5.3.3 Korrektur der Massenerhaltung

Beim Ubergang zwischen den Gittern sind die Erhaltungsgleichungen in der Regel

verletzt. Diese Fehler spielen allerdings keine große Rolle, wie u.a. Untersuchungen

von Meakin (1994) zeigen. Ihre Großenordnung entspricht der von anderen Diskre-

tisierungsfehlern und kann fur die Impulsgleichung und die Transportgleichungen

skalarer Großen vernachlassigt werden.

Die hier gewahlte Behandlung des Druckes erfordert allerdings eine Korrektur der

Massenerhaltungsgleichung. Ist die globale Massenerhaltung verletzt, so verliert das

System der Druckkorrekturgleichungen (3.70) seine Konsistenz.

Hierzu betrachtet man den Massenfluss durch die Flachen, die Interpolationszellen

und die ubrigen aktiven Zellen auf einem Teilgebiet trennen. Diese Overlap-Grenz-

flachen finden sich auf jedem Teil des uberlappenden Systems, sowohl im Over-

lapgitter als auch im Hintergrundgitter. Der Massenfluss durch eine solche Flache

berechnet sich uber folgende Beziehung, vgl. Abb. 5.5:

mi = ρv∗

i · si (5.6)

Abbildung 5.5: Massenfluss an der Overlap-Grenzflache, aus Hadzic (2005) mit Uber-

setzungen des Autors

Fur die Berechnung der Geschwindigkeit v∗

i uber Gleichung (3.66) wird der Wert

der Interpolationszelle verwendet, der im allgemeinen nicht konservativ ist. Es ergibt


sich folgendes Massenungleichgewicht:

∆Molg =∑

mi (5.7)

Hadzic (2005) korrigiert die Massenflusse (5.6) folgendermaßen:

mkori = mi − βi∆Molg (5.8)

βi berechnet sich als entsprechender Anteil am Gesamtfluss durch die Overlap-

Grenzflache Molg =∑

|mi| :

βi =|mi|

Molg

(5.9)

Gleichung (5.7) ist somit null - das fur die Druckkorrektur kritische Massenungleich-

gewicht ist behoben. Die hier gewahlte Form der Massenkorrektur garantiert al-

lerdings nur globale Massenerhaltung. Die Ungleichgewichte auf den Teilsystemen

spielen fur die Losung allerdings keine Rolle, da sie hier ohnehin gemeinsam gelost

werden. Numerische Ergebnisse in Hadzic (2005), Kapitel 6.2 zeigen daruber hinaus,

dass mit Erreichen einer konvergierten Gesamtlosung die Massenungleichgewichte

auf den Teilgebieten nahezu verschwunden sind.

Eine alternative Behandlungsmoglichkeit ist, das Massenungleichgewicht ∆Molg an-

teilig auf die rechten Seiten der Systemgleichungen zu verteilen:

bkorp′ = bp′ + βs∆Molg (5.10)

Sie schneidet in der Praxis allerdings schlechter ab und wird deshalb nicht weiter

verfolgt.

5.3.4 Zellwerte an Overlap-Grenzflachen

Das vorgestellte Verfahren verwendet eine einzige Schicht Interpolationszellen um die

Werte zwischen den Teilgebieten zu koppeln. Die in Kapitel 3 vorgestellten Approxi-

mationsverfahren fur Flusse benotigen z.B. zur Gradientenapproximation zusatzliche

Zellen. Im Umfeld der Grenzflachen stehen diese aufgrund von deaktivierten Zellen

nicht immer direkt zur Verfugung, vgl. Abb. 5.6.

Eine weitere Schicht Interpolationszellen im Bereich der Grenzflachen wurde das Pro-

blem der inaktiven Zellen zwar losen, geht aber mit einem unnotig hohen Berech-

nungsaufwand einher. Alternativ konnte man die benotigten Werte aus dem Inneren


Abbildung 5.6: Gradienten-Approximation an Overlap-Grenzflache, aus Hadzic

(2005) mit Ubersetzungen des Autors

des Gebietes extrapolieren. Hadzic (2005) erzielt mit einem anderen Weg bessere

Ergebnisse:

Fehlende Zellwerte und Gradienten, die fur eine Approximation benotigt werden

und die nicht direkt zur Verfugung stehen, werden nach dem selben Schema wie fur

Interpolationszellen von den ubrigen Systemgittern interpoliert.

5.4 Modellierung bewegter Gitter

In diesem Abschnitt wird das in Kapitel 3 eingefuhrte FV-Verfahren um bewegliche

Zellen erweitert. Es wird das am weitesten verbreitete Verfahren zur gemeinsamen

Behandlung von stationaren und bewegte Gittern vorgestellt: Die Arbitrary Lagran-

gian-Eulerian, kurz ALE-Methode.

5.4.1 ALE-Systemgleichungen

Das ALE-Verfahren ist das Standardmodell zur Realisierung bewegter Gitter. Die

Erhaltungsgleichungen fur Masse (3.1), Impuls (3.7) und skalare Großen (3.10) wer-

den hierzu fur ein bewegtes Bezugssystem formuliert. Sie nehmen folgende Form

an:

d

dt

∫

V

ρ dV +

∮

S

ρ(v − vs) · n dS = 0 (5.11)


d

dt

∫

V

ρui dV +

∮

S

ρui(v − vs) · n dS =

∮

S

µ grad ui · n dS

+

∮

S

µ [iTi · (gradv)T] · n dS

−

∮

S

p iTi · n dS +

∫

V

ρ fbidV (5.12)

d

dt

∫

V

ρφ dV +

∮

S

ρφ(v − vs) · n dS =

∮

S

Γφgradφ · n dS

+

∮

S

qφS · n dS +

∫

V

qφV dV (5.13)

vs ist die Geschwindigkeit, mit der sich die KV-Oberflache S bewegt. Sie wird

auch Gittergeschwindigkeit genannt. Ist sie null so ergibt sich die aus Kapitel 3 be-

kannte eulersche Form der Gleichungen. Ist die Fluidgeschwindigkeit v gleich der

Oberflachengeschwindigkeit vs , so erhalt man die Lagrange’sche Formulierung. Die

Gleichungen fur bewegte KV unterscheiden sich in zwei Aspekten vom stationaren

Fall: Sowohl Zeitableitung als auch konvektiver Fluss mussen bezuglich der relativen

Fluidgeschwindigkeit v−vs bestimmt werden. Beim ortsfesten KV bezeichnet ∂/∂t

die lokale Anderungsrate in einem festen Punkt, berechnet bezuglich der Geschwin-

digkeit v . Hier steht d/dt fur die Anderungsrate in einem sich bewegenden KV, die

sich bezuglich der relativen Geschwindigeit ergibt.

5.4.2 Geometrieerhaltungsgesetz

Die Oberflachengeschwindigkeiten zur Bestimmung der Flusse werden durch Appro-

ximation berechnet. Hierbei werden bekannte Positionen der Zellflachen zu mehreren

Zeitpunkten verwendet. Dadurch kann die Erhaltung der Masse und anderer Großen

verletzt werden, was zu einer Storung der numerischen Losung fuhrt. Die Bildung

dieser zusatzlichen Massequellen oder Senken wird verhindert, wenn folgendes Geo-

metrieerhaltungsgesetz (engl. space conservation law (SCL), geometric conservation

law) erfullt ist:

d

dt

∫

V

dV −

∮

S

vs · n dS = 0 (5.14)


Eine Diskretisierung des SCL nach dem impliziten Euler-Schema fuhrt zu folgender

Gleichung:

∆V nP0

− ∆V n−1P0

∆t=

Nj∑

j=1

∮

Sj

vs · n dS =

Nj∑

j=1

δV nj

∆t(5.15)

∆V nP0

und ∆V n−1P0

sind die Volumina des betrachteten KV P0 zum Zeitpunkt tn

bzw. tn−1 . δV nj bezeichnet das Volumen, das eine Seite j des KV uberstreicht,

wenn es sich in seine neue Position bewegt, vgl. Abb. 5.7. Die Volumenanderung von

P0 entspricht der Summe der von allen Seiten des KV uberstrichenen Volumina:

∆V nP0

− ∆V n−1P0

=

Nj∑

j=1

δV nj (5.16)

Abbildung 5.7: Bewegung eines KV uber drei Zeitschritte mit jeweils uberstrichenem

Volumen, aus Hadzic (2005)

Eine entsprechende SCL-Diskretisierung zweiter Ordnung mit der impliziten Drei-

Ebenen-Methode nimmt folgende Form an:

3∆V nP0

− 4∆V n−1P0

+ ∆V n−2P0

2∆t=

Nj∑

j=1

∮

Sj

vs · n dS =

Nj∑

j=1

3 δV nj − δV n−1

j

2∆t(5.17)

5.4.3 Anderungen in der Diskretisierung

Die Diskretisierung der Systemgleichungen ist weitgehend identisch zu der in Kapitel

3 eingefuhrten Form. Der Term der Anderungsrate muss allerdings in allgemeinerer


Form behandelt werden, um den Volumenanderungen der KV Rechnung zu tragen,

die bei Deformation des Gitters auftreten konnen. Mit dem Euler-Verfahren und der

Drei-Ebenen-Methode ergeben sich folgende Diskretisierungen der Anderungsrate:

d

dt

∫

V

ρφ dV ≈(ρφ∆VP0

)n − (ρφ∆VP0)n−1

∆t(5.18)

d

dt

∫

V

ρφ dV ≈3(ρφ∆VP0

)n − 4(ρφ∆VP0)n−1 + (ρφ∆VP0

)n−2

2∆t(5.19)

Fur unbewegte Gitter oder Bewegungen, bei denen keine Zellen deformiert werden,

reduzieren sich die Gleichungen (5.18) und (5.19) auf die in 3.5.7 eingefuhrte Form.

Auch die Massenerhaltungsgleichung muss in leicht abgewandelter Form behandelt

werden, dazu betrachte man folgende Form von Gleichung (5.11):∮

S

ρv · n dS =

∮

S

ρvs · n dS −d

dt

∫

V

ρ dV (5.20)

Die linke Seite der Gleichung entspricht der Kontinuitatsgleichung des ortsfesten

Falles. Die rechte Seite behandelt die Anteile der Gitterbewegung. Damit die Mas-

senerhaltung gewahrt bleibt, mussen sich die Terme der rechten Seite gegenseitig

ausloschen. Hadzic (2005) wahlt eine Diskretisierung nach dem Schema des oben

behandelten Geometrieerhaltungsgesetzes, sie erfullt diese Bedingung.

Fur die Berechnung geht man wie folgt vor: Der erste Term der rechten Seite von

Gleichung (5.20) wird gemeinsam mit dem Massenfluss berechnet. Fur den Fluss

durch Seite j ergibt sich somit:

mj = ρ(vj − vs) · sj = ρvj · sj − ρδVj

δt(5.21)

Der zweite Term der rechten Seite von Gleichung (5.20) wird gemeinsam mit den

Quelltermen behandelt:

b′P0= bstatP0

+ bsclP0

(5.22)

Die Beitrage des ortsfesten Teils berechnen sich wie in (3.63), fur die Beitrage des

Geometrieerhaltungsgesetzes bsclP0gilt:

bsclP0=

−ρ∆V n

P0− ∆V n−1

P0

∆tbei Euler-Diskretisierung

−ρ3∆V n

P0− 4∆V n−1

P0+ ∆V n−2

P0

2∆tbei Drei-Ebenen-Methode

(5.23)


5.5 Aktuelle OLG-Software

In diesem Abschnitt wird kurz auf CFD-Softwarepakete mit Overlapping Moving

Grid-Funktionalitat eingegangen.

5.5.1 Tau (DLR)

Die deutsche Gesellschaft fur Luft- und Raumfahrt (DLR) verwendet Tau, einen

dreidimensionalen Finite-Volumen-Loser fur RANSE auf unstrukturierten, uberlap-

penden Netzen. Madrane u. a. (2002) haben hierfur ein Overlapping Moving Grid

Verfahren basierend auf Benek u. a. (1986) entwickelt. Eine Weiterentwicklung von

Madrane u. a. (2004) ermoglicht volle Parallelisierung von Rechnungen uber das MPI-

Protokoll.

5.5.2 Chimera Grid Tools (NASA, Air Force)

Aus den Chimera-Verfahren von Benek u. a. (1986), die im Auftrag der NASA ent-

wickelt wurden, ist das Chimera Grid Tool Paket hervorgegangen. Hierbei erzeugt

das Preprocessing-Modul PEGSUS aus einer Eingabe von uberlappenden Teilgittern

ein Gesamtnetz mit Koppelungsinformationen, das dann mit einem Stromungsloser

(OVERFLOW-D) behandelt werden kann. Anwendungen und Ergebnisse siehe Mea-

kin und Wissink (1999). Eine Analyse der MPI-Parallelisierung findet sich bei Djo-

mehri und Biswas (2003).

5.5.3 Overture (LLNL)

Zeitgleich wurde auf Basis der Arbeit von Chesshire und Henshaw (1990) am kalifor-

nischen Lawrence Livermore National Laboratory (LLNL) das Softwarepaket Over-

ture entwickelt, siehe Brown u. a. (1999). Hierbei handelt es sich um eine Sammlung

von C++ Bibliotheken zur Losung von PDEs auf uberlappenden und nichtuberlap-

penden, strukturierten Gittern mit FV- und FD-Methoden. Die Bibliotheken un-

terstutzen u.a. MPI-Parallelisierung und adaptive Gitterverfeinerung, aktuelle Er-

gebnisse siehe bei Henshaw (2008).


5.5.4 FEFLO (NRL)

Lohner u. a. (2002) haben am Naval Research Laboratory in Washington einen par-

allelen FE-Loser fur kompressible und inkompressible Navier Stokes Gleichungen mit

uberlappenden, unstrukturierten Gittern entwickelt.

5.5.5 Comet (CD-adapco)

Die Implementierung uberlappender Gitter nach Hadzic (2005) im kommerzielle FV-

Loser Comet der Firma CD-adapco soll in den folgenden Abschnitten dieser Arbeit

untersucht werden.

Kapitel 6

Anleitung fur Overlap-Rechnungen

mit Comet

Dieses Kapitel beschreibt die Anwendung des Overlapping Grid Verfahrens fur CFD-

Simulationen in Comet, einem kommerziellen Softwarepaket der Firma CD-adapco,

das auf C und Fortran77 basiert.

Die Beschreibung gliedert sich folgendermaßen: Nach einer Einfuhrung der benotig-

ten Programmbefehle folgen generelle Richtlinien fur den Aufbau und die Struk-

turierung von Overlap-Regionen. Besonderheiten im Falle bewegter Gitter werden

genannt, und die Modellierung von Anbauteilen und hierbei auftretende Schwie-

rigkeiten beschrieben. Ein weiterer Abschnitt widmet sich der Parallelisierung von

Overlap-Rechnungen, die in der vorliegenden Version 2.360.Z des Pakets noch auf

zwei Knoten beschrankt ist. Zum Abschluss des Kapitels wird ein Uberblick uber

bekannte Fehlermeldungen des OLG-Moduls gegeben. Neben den Situationen, in de-

nen sie auftreten, werden die Ursachen der Probleme genannt - sofern sie bekannt

sind. Daruber hinaus werden mogliche Losungsstrategien diskutiert.

6.1 OLG-Programmbefehle

Das Programmpaket besteht aus dem Pra- und Postprozessor Cometpp mit dem Si-

mulationsmodelle erstellt, sowie Ergebnisse des Stromungslosers Comet visualisiert

werden konnen. Uber das Usercoding, Fortran77-Programmcode, der von Comet zwi-

schen den Zeititerationen abgearbeitet wird, hat der Benutzer die Moglichkeit, Ein-

fluss auf die Rechnungen zu nehmen. Beispielsweise wird die Bewegung der Gitter

uber solche Usercoding-Routinen durchgefuhrt.

70

KAPITEL 6. ANLEITUNG FUR OVERLAP-RECHNUNGEN MIT COMET 71

An dieser Stelle sollen die Befehle der OLG-Funktionalitat vorgestellt werden, geglie-

dert in Pra-, Postprozessor- und usercoding-Anweisungen. Das Zeichen < · > steht

hierbei fur einen zu setzenden Parameter, der aus einer Liste der Form ( · | · | . . . )

gewahlt werden kann.

6.1.1 Praprozessor-Befehle

Zelltyp

Die Zellen eines Overlap-Korpers mussen in einer eigenen, zusammenhangenden Sub-

domain liegen, die uber den Befehl subdomain definiert und aktiviert wird:

subdomain < Index >< Name >

Dann mussen die Overlap-Zellen der Subdomain uber den Befehl ctype zugeordnet

werden:

ctype < Index >< Farbe >< Domain >< Subdomain >< Gruppe >< Name >

Ublicherweise behalt das Hintergrundgitter Subdomain 1, die ubrigen Overlap-Korper

werden ab Index 2 aufsteigend nummeriert.

Zu beachten ist hierbei:

• Comet unterstutzt momentan nicht mehr als 9 Subdomains.

• Regionen sind auf eine Subdomain beschrankt. Die Verwendung der selben Re-

gionsnummer auf mehreren Subdomains fuhrt zu einem Fehler bei Ausfuhrung

des Befehls geom.

Overlap-Randbedingung

Den Overlap-Rand eines Korpers (vgl. Kapitel 5.2) definiert man fur eine Region

folgendermaßen:

bdef < Regionsnummer > over

< Suchschema >< Interpolationsordnung >< Flux >< Layer >

Folgende Parameter sind zu setzen:

• Fur das Suchschema steht (linear < Grad > | cons | tight) zur Auswahl.

tight entspricht dem in Abschnitt 5.2.3 eingefuhrten Verfahren zur Spenderzel-

lenbestimmung.


• Die Interpolationsordnung kann (0|1) gesetzt werden.

0 enspricht direkter, 1 linearer Interpolation.

• Der Parameter Flux ist aus (0|1|2) wahlbar. Er bezeichnet die Form der Fluss-

korrektur um die Massenungleichgewichte zu beseitigen, vgl. Abschnitt 5.3.3.

0 deaktiviert die Funktion, 1 entspricht der Variante der Korrektur der rechten

Seite der Systemgleichungen, 2 entspricht der ausfuhrlich beschriebenen Form.

• Die Layer-Option steht standardmaßig auf 1 und wurde nicht weiter untersucht.

Gute Ergebnisse wurden mit den Parametern zur Auswahl des jeweiligen Verfahrens

hochster Ordnung erzielt:

bdef RegNr over

tight 1 2

Zu Beachten ist daruber hinaus folgendes:

• Die Overlap-Regionen sollten erst gesetzt werden, wenn alle Modifikationen

am Hintergrundgitter abgeschlossen sind. Das nachtragliche Ausschneiden von

Streben aus dem Hintergrundgitter hat zu einer nicht anlaufenden Rechnung

gefuhrt. Die Ursache ist jedoch unklar.

LGS-Loser

Der Loser fur nicht symmetrische Gleichungssysteme, wie sie im OLG-Fall auftreten,

wird uber den Befehl solver gewahlt:

solver < SymOption >< NichtSymOption >< Vorkonditionierer >

• Fur symmetrische Matrizen stehen (AMG | CG | CGSTAB) zur Verfugung:

Adaptives Multigrid, Konjugierte Gradienten und das Stabilisierte Bi-CG-Ver-

fahren, vgl. Abschnitt 3.7.3.

• Fur nicht symmetrische Matrizen ist (CGSTAB) die einzige Wahlmoglichkeit.

• Mogliche Vorkonditionierer sind (IC | J) : Incomplete Cholesky und Jacobi.

Uber iter konfiguriert man den Loser:

iter < MaxIt >< ResidTol >< MinIt >


• MaxIt bestimmt die maximale Anzahl außerer Iterationen pro Zeitschritt.

• ResidTol ist die Abbruchschranke fur Konvergenz.

• MinIt setzt eine minimale Anzahl durchzufuhrender außerer Iterationen.

Transiente Rechnung

Eine zeittransiente Berechnung hat sich als deutlich robuster fur OLG-Rechnungen

erwiesen als die zeitstetige Alternative, wie sie fur Falle ohne bewegte Gitter moglich

ware. Die implizite Drei-Ebenen-Methode wurde deshalb sowohl fur stationare als

auch instationare Gitter verwendet. Die Methode wird uber den Befehl time gesetzt:

time < Problemtyp >< Zeitschritte >< Schrittweite >< Schema >

• Der erste Option wahlt aus (steady | transient | pseudotrans) den Typ des

Problems.

• Uber Zeitschritte wahlt man die Anzahl der außeren Zeititerationen.

• Die Schrittweite bestimmt die Lange eines solchen Zeitschrittes.

• An Zeitintegrationsschemata stehen (ieuler | ittl) , Implizites Euler und die

implizite Drei-Ebenen-Methode zur Auswahl.

Wie in Abschnitt 3.7.3 behandelt wird der CGSTAB-Algorithmus mit IC-Vorkondi-

tionierung und der Drei-Ebenen-Methode verwendet. Ublicherweise setzt man funf

außere Iterationen pro Zeitschritt, zusammen mit einer hohen Abbruchschranke. Die

Anzahl innerer Iterationen sind von Comet fest vorgegeben und konnen nicht mo-

difiziert werden (Impuls 200, Masse 30, Turbulenz 50). Fur die Konfiguration des

LGS-Losers ergeben sich folgende Befehlsaufrufe (hier fur 1800 Zeitschritte, wobei

die Schrittweite der Zeit entspricht, die der VSP fur eine Drehung um 1 benotigt):

%Solver-Konfiguration

solver amg cgstab ic

iter 5 1e-13

time tran 1800 1.2503126e-3 ittl


6.1.2 Postprozessor-Befehle

Visualisierung der Overlap-Zellfunktionen

Zu Beginn eines Zeitschrittes wird die Funktion der Zellen im Overlapbereich be-

stimmt, vgl. Abschnitt 5.2.1. Mit folgenden Befehlen kann sie farblich dargestellt

werden:

cset news . . .

% auf Konturplot umstellen:

popo all add cont

% Zellfunktionsdaten laden:

stlo 4 iovc

plex

Siehe Tabelle 6.1 fur eine Aufstellung der auftretenden Farbindices und Funktionen

und Abb. A.5 fur ein Beispiel.

Tabelle 6.1: Farbtabelle der Zellfunktionen

Index Funktion

-1 inaktive Zelle, d.h. kein Teil der Stromungslosung

0 aktive Zelle außerhalb des Overlap-Bereiches

d.h. normaler Bestandteil der Stromungslosung

1 Spenderzelle - von diesen Zellen interpoliert das Overlapgitter

2 Zwischenschicht - aktive Zelle innerhalb des Overlap-

Bereiches ohne Sonderfunktion

3 Interpolationszelle des Hintergrundgitters

4 Interpolationszelle eines Overlap-Korpers

5 Ubergangszelle - ein zusatzlich berechneter Wert

z.B. fur Overlap trifft Overlap

Im Falle des Absturzes einer Stromungsrechnung wahrend eines bestimmten Zeit-

schrittes konnen die iovc-Funktionen der Zellen trotzdem uberpruft werden: Durch

eine Anderung der außeren Iterationen OutIt auf 0 in der prob Datei bestimmt Comet

in den folgenden Zeitschritten nur die Overlap-Funktionen der Zellen.


Aktive Zellen der Stromungslosung

Uber den Befehl cset wahlt man die bei der Stromungslosung aktiven Zellen aus:

cset active

6.1.3 Usercoding-Befehle

Debugging-Ausgaben

Ein Zusatz im Usercoding fuhrt zu erweiterten Kontrollausgaben fur OLG-Rechnun-

gen. Die Datei useinp.f im usercoding-Ordner muss folgende Zeile enthalten:

ioverdbg = 3

Es empfiehlt sich, die Bildschirmausgabe in eine Datei umzuleiten.

Fur eine Liste der Fehlermeldungen siehe Tabelle 6.2 in Abschnitt 6.6, hier werden

auch Ursache und Losung der Probleme behandelt.

Auf manchen Systemen werden im Verlauf der Rechnung die Dateien fort.74 und

fort.75 angelegt. Sie konnen uber 100 MB groß werden. Der Inhalt ist nicht von

Bedeutung, sie konnen im Anschluss an die Rechnung geloscht werden.

Gitter-Koppelung

Die Erweiterung der Datei useinp.f im usercoding-Ordner um die Zeile

lhierarch=.true.

soll eine Verbesserung der Gitter-Koppelung bei Mehrfachuberlappungen bewirken.

Die Koppelungsprobleme bei hier durchgefuhrten Simulationen konnten durch diesen

Zusatz allerdings nicht gelost werden.

6.2 Richtlinien fur die Gitterstruktur im Overlap-

Bereich

Die Untersuchungen, die im Rahmen dieser Arbeit durchgefuhrt wurden, kommen

zu folgenden Ergebnissen bezuglich Aufbau und Strukturierung eines uberlappenden

Gittersystems:


Anzahl uberlappender Zellschichten

Der Overlap-Korper und das Hintergrundgitter sollten sich in mindestens vier Zell-

schichten uberlappen. So viele Schichten werden benotigt, um zu verhindern, dass

eine Gitterzelle gleichzeitig Interpolationszelle und Spenderzelle ist. Werden nicht ge-

nug Zellen gefunden um diese Trennung zu erfullen, so wird auf Verfahren niedrigerer

Ordnung fur die Stromungsberechnung zuruckgegriffen.

Zellgroßenverhaltnis

Vom Overlap-Rand aus einwarts gehend sollten drei bis vier Zellschichten auf beiden

Gittern eine ahnliche Kantenlange aufweisen. Große Differenzen der Zellgroßen von

einem Gitter zum andern fuhren zu einer schlechteren Koppelung der Losung oder

sogar zum Absturz der Rechnung.

Wenn ein Overlap-Rand auf die Nahtstelle einer Zellverfeinerung trifft (econ refi nach

einem Aufruf von cref cset) kommt es haufig zu Absturzen. Die Verfeinerung eines

Gebietes, das den gesamten Bewegungsbereich des Overlap-Korpers einschließt, ist

eine sicherere Alternative.

Randstruktur und Bewegung

Bei bewegten Gittern mit einer sehr feinen Struktur im Overlap-Randbereich kann

es zwischen zwei Zeitschritten - relativ zur Zellgroße gesehen - zu einer großen Bewe-

gung kommen. Bei einer zu feinen Auflosung des Overlap-Randes ist der Algorithmus

nicht mehr in der Lage, die Teilgebiete zu koppeln.

In diesem Falle muss der Zeitschritt, und damit die Bewegung, verkleinert werden.

Alternativ kann ein Overlap-Gitter verwendet werden, das neben der notigen Feinheit

der Vernetzung des Korpers auch die Vergroberung zu einer gleichmaßigen Außen-

struktur aufweist, wie sie fur eine gute Koppelung notig ist.

Kritische Randsituationen

Bei Overlap-Rechnungen kommt es zu folgenden Konstellationen:

• Overlap-Korper trifft auf Gebietsaußenrand : Bei den hier durchgefuhrten Rech-

nungen hatte nur die Oberkante des Overlap-Korpers Kontakt zu den Grenzen


des Hintergrundgebietes. Ein durchgangiger Rand des Korpers vom Typ over

fuhrte hier zu Problemen. Diese Grenzflache zum Außenrand mit der selben

Randbedingung wie dem Rest des Außernrandes zu versehen (slip wall) loste

das Problem. Die ubrigen Randflachen des Overlap-Korpers als over zu defi-

nieren fuhrte zum erwunschten Ergebnis.

• Overlap-Korper trifft auf Loch im Hintergrundgitter : Wenn ein Overlap-Rand

großflachig auf einen Hohlraum trifft, so kommt es sehr haufig zu einem Ab-

sturz der Rechnung. Unter gewissen Umstanden ist das Uberstreichen eines

Hohlraumes moglich, einige Falle werden spater vorgestellt. Kritische Parame-

ter sind eine Kombination aus Zellgroßenverhaltnissen, Anzahlen von Overlap-

Schichten und Große der Aussparung im Gitter. Erfolge in diesem Bereich ba-

sieren aber eher auf try and error als auf sicheren Erkenntnissen.

Um Absturze sicher zu verhindern darf der Rand des Overlap-Korpers nicht

in einen Hohlraum hineinreichen, sondern hochstens bis in die Randschicht der

Aussparung.

• Overlap-Korper trifft auf Overlap-Korper : Fur die Koppelung zweier oder meh-

rerer Overlap-Korper auf einem Hintergrundgitter gelten neben den ublichen

Regeln einige Besonderheiten. Kontakt von over gesetzten Randern ist kein

Problem, der Algorithmus bestimmt dann das aktive Losungsgebiet als Kom-

bination von Hintergrund- und Overlap-Gitterzellen.

Bei einem Hohlraum in einem Overlap-Korper erhoht sich die Anzahl der noti-

gen Zellschichten zwischen Aussparung und dem Overlap-Rand eines anderen

Korpers. Anders als im eben behandelten Fall einer Aussparung im Hinter-

grundgitter mussen drei komplette Zellschichten Abstand gewahrt bleiben, um

Nebeneffekte zu verhindern: Ein kleinerer Abstand fuhrt nicht unweigerlich zu

einem Absturz, aber haufig werden unerwunschte Zellen eines anderen Gitters

im Hohlraum aktiviert, die dann falschlicherweise zur Stromungslosung beitra-

gen.

Econ-Klebeflachen

Drei Situationen sind bei der Verbindung von Teilgittern mit dem Befehl econ zu

unterscheiden:

An Klebeflachen zwischen Regionen (econ regi), die am Ubergang ahnliche Zellgroßen

aufweisen, kommt es zu keinen Schwierigkeiten beim Kontakt mit Overlap-Korpern.


Beruhrungen des Overlap-Randes mit einer Verfeinerungs-Klebeflache (econ refi)

sind zu vermeiden: aufgrund der Zellgroßenunterschiede kommt es hier bei Rech-

nungen haufig zu Absturzen.

Dynamische Sliding-Flachen (econ slid), an denen fur jeden Zeitschritt die aktuellen

Partnerzellen der Teilgebiete gekoppelt werden, lassen sich nicht mit einem Overlap-

Rand kombinieren. Die Rechnung sturzt bei Kontakt sofort ab.

6.3 Parallelisierung

Eine Parallelisierung auf mehr als zwei Knoten ist in Comet momentan fur Rech-

nungen mit OLG nicht moglich (Stand Version 2.360.Z). Diese Einschrankung gilt

hierbei nicht nur fur die Zerlegung des Bereiches, der den Overlap-Korper enthalt,

sondern fur das gesamte Losungsgebiet.

Fur eine Verteilung der Zellen auf zwei Knoten hat man folgende Moglichkeiten:

• Das Overlapgebiet wird vollstandig auf einen der beiden Knoten gelegt. Bei

bewegten Gittern muss das Gebiet den kompletten Bereich enthalten, den der

Overlap-Korper im Verlauf der Rechnung uberstreicht. Dies ist vergleichbar

mit dem Prinzip, dass ein Sliding-Korper in einer parallelisierten Rechnung

komplett auf einem Knoten liegen muss.

• Eine Aufteilung des Overlap-Bereiches kann, anders als im Sliding-Fall, nach

gewissen Regeln erfolgen: Solange die Bewegung der Overlap-Korper in einer

Ebene erfolgt, kann das Gebiet bezuglich dieser Ebene getrennt werden. Hierzu

mussen alle Zellen an der Trennflache auf eine einheitliche Koordinate gebracht

werden (vmod vset), vgl. Abb. A.12.

• Eine Kombination beider Verfahren ist auch moglich.

Automatische Overlap-Partitionierung wird nicht unterstutzt. Die Verteilung auf die

Knoten erfolgt durch manuelle Partitionierung nach entsprechender Zellauswahl:

%Partitionierung

cset news ...

% 1. Teilpartition

ndef manu 1 1 cset

% 2. Teilpartition


cset inve

ndef manu 1 2 cset

Eine moglichst gleichmaßige Verteilung der Zellen auf die beiden Knoten liefert die

besten Ergebnisse.

6.4 Simulation des VSP mit Overlapping Moving

Grid

Ausgehend vom Sliding-Interfaces-Fall mussen fur den Umstieg auf eine Simulation

des VSP mit Overlapping Moving Grids folgende Anderungen durchgefuhrt werden:

6.4.1 Anderungen der Vernetzung

• Jeder Flugel wird als eigener Overlap-Korper behandelt.

• Der Wingtip, der bisher den Flugelzylinder uber ein Sliding-Interface nach

unten abgeschlossen hat, wird fest mit dem Flugelgitter verbunden.

• Die strenge Zylinderform des Flugelgitters muss nicht mehr beibehalten wer-

den, fur die OLG-Rechnungen kann der Flugelzylinder ohne Qualitatsverlust

auf ein Oval reduziert werden.

• Die Radkorper-Struktur mit den Sliding Interfaces wird nicht langer benotigt.

An ihre Stelle tritt ein neues, durchgangiges Hintergrundgitter. Die Feinheit

des Gitters im Bereich, in dem sich die Flugel befinden, ist fur eine optimale

Koppelung an die Zellgroßen der Overlap-Gitter angepasst.

Die feine Gitterstruktur an der Flugelhinterkante, bedingt durch die Zylinderbauwei-

se, ist aufgrund der Schwierigkeiten bei der Koppelung unterschiedlich großer Zellen

nicht optimal fur OLG-Rechnungen. Auch beim Bodenstuck, das den Flugel nach

unten abschließt, muss der Ubergang vom feinen Flugelprofil zu einer gleichmaßigen

Zellgroße am Rand noch verbessert werden.

Fur die hier durchgefuhrten Simulationen wird eine vergroberte Form des Flugelzy-

linders verwendet, bei der die Zellgroßen im Randbereich nicht zu stark variieren. Er

wurde mit gewonnenen Erkenntnissen im Verlauf der Arbeit an die jeweiligen Bedurf-

nisse angepasst. Fur optimale Ergebnisse ist jedoch ein Neudesign des Flugelgitters

bezuglich dieser Aspekte notwendig.


6.4.2 Anderungen im Usercoding

Das Usercoding, das die Bewegung der Flugel nach der FWK steuert, bleibt weitge-

hend erhalten. Zwei kleine Anderungen mussen durchgefuhrt werden:

• Die Aufrufe in der Datei umovgr.f zur Drehung des Radkorpers sind uber-

flussig.

• Die Randbedinungen der Wingtip-Flachen mussen nicht mehr uber die Datei

userbc.f gesetzt werden.

6.4.3 Simulationsrechnungen

Im Rahmen des Validierungsprozesses des OLG-Verfahrens, der in Abschnitt 7.1

beschrieben ist, wurden folgende Simuationen durchgefuhrt:

• Stationare Flugelumstromungen: Die Krafte, die bei Anstromung von vorne

auf einen unbewegten Flugel wirken, wurden mit Sliding Interfaces und einem

Overlapping Grid Modell berechnet. Fur unterschiedliche Anstellwinkel wurden

die Kraftkomponenten in X-Richtung, Y-Richtung sowie das Drehmoment um

die Z-Achse verglichen.

• Einflugler : Ein VSP mit einem Flugel wurde als OLG- und Sliding-Fall berech-

net. Es wurden Schub, Querkraft und Wirkungsgrad verglichen.

• Funfflugler : Fur einen VSP mit funf Flugeln wurden neben einer Untersuchung

von Schub, Querkraft und Wirkungsgrad eines einzelnen Flugels auch ein Ver-

gleich der uber alle Flugel summierten Krafte durchgefuhrt, die sich aus beiden

Varianten der Modellierung ergaben.

Bei diesen Fallen wurde durchweg eine gute Deckung der Ergebnisse erzielt. Im wei-

teren Verlauf dieser Arbeit wurden durch die Hinzunahme von Anbauteilen deutlich

komplexere Situationen untersucht.

6.5 Modellierung von Anbauteilen

Das OLG-Verfahren soll die Modellierung von Anbauteilen wie Schutzplatte und

Streben im direkten Umfeld des VSP ermoglichen. Grundsatzlich stehen hierfur zwei


Techniken zur Verfugung: Die Behandlung eines Anbauteils als eigenen Overlap-

Korper oder die Modellierung durch eine Aussparung im Hintergrundgitter. Fur an-

einander grenzende Teile wie Streben und Schutzplatte stellt sich auch die Frage der

Kombinierbarkeit dieser Verfahren.

Hier sollen nun die Erkenntnisse aus Versuchen in diesem Bereich vorgestellt werden,

die im Rahmen dieser Arbeit durchgefuhrt wurden.

6.5.1 Streben als Aussparung im Hintergrundgitter

Fur die ersten Versuche wurde eine vereinfachte, quadratische Strebe mit einer

Grundflache von 2x2 Zellen aus dem Hintergrundgitter geschnitten, vgl. Abb. A.15

(u.l.). Es wurden verschiedene Abstande und mehrere Verfeinerungsstufen des Lochran-

des mit folgenden Ergebnissen getestet:

• Wird eine zu grob vernetzte Aussparung mit dem Flugelgitter uberstrichen,

fuhrt dies zu einem Absturz der Rechnung.

• Eine Verfeinerung des Randes und Umfeldes der Strebe (cref) ermoglicht ein

Eindringen der Strebe in den Randbereich des Flugelgitters.

Weitere Untersuchungen mit einem großeren, ovalformigen Strebenkorper fuhrten zu

folgendem Resultat:

• Großflachiger Kontakt des Flugelgitter-Overlap-Randes mit der Aussparung

fuhrt trotz passender Zellgroßenverhaltnisse zu einem Absturz der Rechnung,

vgl. Abb. A.13 (u.r.).

6.5.2 Streben als Overlap-Korper

Fur das erste Overlap-Gitter wurde ein Block des Hintergrundgitters dupliziert und

aus dessen Mitte wie im oben untersuchten Fall eine quadratische Strebe mit einer

Grundflache von 2x2 Zellen ausgeschnitten, vgl. Abb. A.15 (o.l.). Es wurden mehrere

Blockgroßen und Verfeinerungsstufen des Gitters mit folgenden Ergebnissen getestet:

• Kontakt des Overlap-Randes der Strebe mit dem Flugelprofil fuhrt zum sofor-

tigen Absturz der Rechnung.


• Das Zellgroßenverhaltnis ist ein sehr kritischer Parameter fur die Stabilitat

einer Rechnung mit mehreren uberlappenden Overlap-Korpern auf einem Hin-

tergrundgitter.

• Wenn ein Overlap-Gebiet bis zum Rand verfeinert wird, fuhrt dies haufig zu

einem Absturz der Rechnung. Bei Verfeinerungen, die erst im Inneren des Over-

lap-Gitters beginnen (zwei Zellschichten Abstand zum Rand), treten diese Pro-

bleme nicht auf.

• Am oberen Ende der Strebe kommt es zu einer ungewollten Aktivierung von

Zellen im ausgesparten Bereich, wenn Flugelgitter und Strebengitter bis in die-

sen Hohlraum hinein uberlappen. Eine Verfeinerung des Strebengitters redu-

ziert zwar die Dicke dieser Schichten aus Hintergrund- und Flugelgitterzellen,

beseitigt das Phanomen aber nicht. Die Ursache konnte nicht endgultig geklart

werden.

6.5.3 Schutzplatte als Aussparung im Hintergrundgitter

Eine erste Schutzplatte wurde als kreisformige Aussparung im Hintergrundgitter ge-

neriert. Ein Abstand zwischen Flugelspitze und Schutzplatte von 1% des Flugelkreis-

durchmessers sollte fur Simulationsrechnungen realisiert werden. Dieses Ziel wurde

unter Gewinnung folgender Erkenntnisse erreicht:

• Falls der Flugel-Overlapkorper in die Aussparung hineinragt, kommt es zu

einem Absturz der Rechnung. Das Flugelgitter ist so zu kurzen, dass es keinen

Kontakt zum Hohlraum hat.

• Gleichzeitig mussen die Zellen unterhalb der Flugelspitze und das Hintergrund-

gitter entsprechend verfeinert werden, um eine hinreichende Anzahl uberlap-

pender Zellen der Gitter zu gewahrleisten.

Fur das Hintergrundgitter wurde hierfur eine cref-Verfeinerung in Z-Richtung so

durchgefuhrt, dass acht Zellschichten zwischen Oberkante der Schutzplatte und Hohe

der Flugelspitze liegen.

Das Flugelgitter wurde so verkurzt, dass der Flugel durch acht Zellschichten nach

unten abgeschlossen wird. Der Boden des Flugelkorpers reicht hochstens bis in die

Zellschicht des Hintergrundgitters, welche die Schutzplatte nach oben begrenzt, vgl.

Abb. A.14.


Weitere Versuche wurden mit der gegebenen Vernetzung einer realistischen Schutz-

platte durchgefuhrt, die ins Hintergrundgitter eingefugt wurde (econ regi):

• Dieses Schutzplattengitter wies eine grobe und ungleichmaßige Vernetzung

der Oberflache auf. Eine sinnvolle Koppelung mit dem Flugelgitter war nicht

moglich. Eine Verfeinerung mit cref im Bereich der Bewegung des Flugelgitters

war notwendig, um ein gutes Ergebnis zu erzielen.

• Eine geeignete Wahl des Verfeinerungsbereiches der Schutzplatte gestaltete sich

allerdings schwierig: Eine Verfeinerung der gesamten Region ergab Bereiche

mit zu kleinen Zellen fur die ubliche Bewegung von 1 pro Zeitschritt. Eine

kleinere Wahl des Gebietes fuhrte anfangs zu Problemen aufgrund der Zell-

großensprunge, die an econ refi-Klebeflachen im Bewegungsbereich des Flugels

auftraten.

6.5.4 Schutzplatte als Overlap-Korper

Ein Block verfeinerter Hintergrundzellen ist die Basis einer Schutzplatte im OLG-

Verfahren. Die vereinfachte Schutzplatte wird nicht aus dem Hintergrundgitter, son-

dern aus den duplizierten Zellen ausgeschnitten. Die Analyse dieses Falles ergab:

• Eine Ausdehnung eines Overlap-Korpers bis in Aussparungen im Nachbargitter

ist zu vermeiden. Es kommt zu Problemen bei der Bestimmung der aktiven

Zellen.

• Der Abstand des Overlap-Randes zur Schutzplatten-Aussparung muss deutlich

vergroßert werden: Im nicht uberlappenden Fall der Schutzplatte konnte der

Flugelkorper bis in die Randschicht der Aussparung ragen. Hier ist ein Ab-

stand von drei kompletten Zellschichten notwendig, damit dem Algorithmus

eine Koppelung der drei Gitter gelingt. Eine weitere Verfeinerung des Hinter-

grundgitters ist im Ubergangsbereich erforderlich, vgl. Abb. A.14.

6.5.5 Overlap-Aussparungs-Kombinationen

Aus der Kombination beider Verfahren ergab sich Folgendes:

Strebe als Overlap-Korper, Schutzplatte als Aussparung im Hintergrund :


• Eine stationare, vereinfachte Strebe kann in den Hohlraum hineinragen, ohne

dass es zu einem Absturz kommt.

• Es kommt zu unerwunschten Zellschichten in der Strebenaussparung am Uber-

gang. Die Modellierung einer durchgangigen Aussparung zwischen Strebe und

Schutzplatte ist fur Mischfalle der Modellierung nicht gelungen.

• Bei Kontakt mit dem Flugelgitter kommt es nicht nur an der Oberkante, son-

dern auch an der Unterkante der Strebe zu einigen Schichten unerwunscht

aktiver Zellen im Hohlraum, vgl. Abb. A.15 (u.).

Schutzplatte als Overlap-Korper, Strebe als Aussparung im Hintergrund :

• Es ist nicht gelungen, eine solche Verbindung herzustellen. Die ausgeschnitte-

nen Zellen am Ubergang von Strebe zu Schutzplatte verursachen einen Absturz

beim Versuch, das fur die Losung aktive Gebiet zu bestimmen.

6.5.6 Overlap-Overlap-Kombinationen

• Bei der Modellierung durch separate Overlap-Gebiete bilden sich am Ubergang

zwischen Strebe und Schutzplatte mehrere Zwischenschichten aus aktiven Hin-

tergrundzellen, vgl. Abb. A.15 (o.r.). Dieses Phanomen tritt unabhangig von

der Flugelposition auf.

• Ragt ein Strebenteil bis in die Aussparung der Schutzplatte, so bleiben diese

Zellen fur die Rechnung aktiv.

• Bei einer gemeinsamen Vernetzung von Schutzplatte und Strebe gelingt ein

sauberer Ubergang zwischen den beiden Teilgebieten. Allerdings kommt es bei

Kontakt von Flugelgitter und Aussparung wieder zu aktiven Hintergrundzellen

im Strebenkorper.

6.5.7 Aussparungs-Kombinationen

Die Kombination der vereinfachten Streben und der Schutzplatte als Aussparungen

im Hintergrundgitter siehe in Abb. A.15 (o.l.). Fur diesen Fall ergaben sich allerdings

keine neuen Erkenntnisse uber die Einzelfalle hinaus.

Daruber hinaus wurden Kombinationen von realistisch geformten Schutzplatten- und

Strebengittern getestet:


• Die komplexeren Gitter werden in entsprechende Aussparungen des Hinter-

grundgitters geklebt (econ regi), vgl. Abb. A.16.

• Klebeflachen am Ubergang (econ regi) ermoglichen eine Kombination von Teil-

netzen.

Die Ergebnisse sind vielversprechend. Eine luckenlose Modellierung des Ubergangs

zwischen Schutzplatte und Strebe war aufgrund der Gitterstruktur der Schutzplatte

in den Simulationen nicht moglich, vgl. Abb. A.17. Auf das notwendige Neudesign

der Teilnetze wurde im Rahmen dieser Arbeit aus Zeitgrunden verzichtet.

6.6 Problembehandlung

Tabelle 6.2 fasst Fehlermeldungen und den jeweiligen Bereich des Problems zusam-

men. In diesem Abschnitt soll auf die Situation ihres Auftretens und bekannte Losun-

gen eingegangen werden.

Tabelle 6.2: Fehlermeldungen OLG

Meldung Problembereich

segmentation fault Suchalgorithmus:

Absturz (Ursache unbekannt)

dirsearch <= 0 ; Suchalgorithmus:

searchinter no couple subdomain Bestimmung Spenderzellen

something wrong in boundary; Zellfunktionen:

something wrong in close Konflikte in der erfolgten Zuordnung

not enough active neighbors; Gescheiterte Spenderzellensuche:

inactive donor cells; Zu große Bewegung pro Zeitschritt resultiert

switch from inactive to . . . in ungenugender Zahl aktiver Spenderzellen

NaN-Werte Numerik-Versagen


6.6.1 Bekannte Probleme und Losungen

Suchalgorithmus

Probleme des Suchalgorithmus zur Bestimmung von Spenderzellen fuhren zu Feh-

lermeldungen der Funktion direction search. Angezeigt wird eine Meldung dirsearch

oder dirse <= 0 . Dies passiert haufig im ersten Zeitschritt, wenn eine erste Kop-

pelung der Gitter hergestellt werden muss. Eine Verbesserung der Gitterstruktur

bezuglich Anzahl uberlappender Zellschichten und besserer Anpassung der Zell-

großen lost das Problem in der Regel. Haufig gelingt dem Algorithmus trotz anfang-

licher Fehlermeldungen eine gultige Koppelung der Teilgitter.

In seltenen Fallen kommt es zu einem vollstandigen Versagen des Suchalgorithmus.

Dies resultiert in einem Absturz mit Meldung segmentation fault, eine Ursache konnte

nicht gefunden werden.

Die genaue Ursache der Fehlermeldung searchinter no couple subdomain konnte nicht

geklart werden.

Im Falle einer Rechnung, die ohne eine Fehlermeldung nicht uber den Anfang der

Koppelungsprozedur hinauskam, brachte eine Parallelisierung auf zwei Knoten den

gewunschten Erfolg. Die Ursache blieb jedoch unklar.

Zellfunktions-Konflikte

Meldungen mit Inhalt something wrong in . . . weisen auf Konflikte in der Zuord-

nung der Zellfunktionen bezuglich OLG hin: Aufgrund der Ungleichmaßigkeit der

Gitter oder der Komplexitat der Situation schafft es der Algorithmus nicht, das ak-

tive Gebiet korrekt zu bestimmen. Es gelingt nicht, ein abgeschlossenes Gebiet des

Hintergrundgitters auszublenden, oder einen geschlossenen Overlap-Korper aktiv zu

schalten.

Ist kein sinnvolles Gitter fur die Rechnung aktiv, so kommt es meist innerhalb des

nachsten Zeitschrittes zu einem Versagen der Numerik.

Haufig gelingt dem Algorithmus allerdings trotz dieser Meldung die Bestimmung ei-

ner gultigen Koppelung. Zur Uberprufung der Gitterkombination bleibt neben der

optischen Analyse der IOVC-Verteilung im entsprechenden Gebiet die erweiterte

Fehlerausgabe: Unter der Meldung invalid jumps wird fur jeden Zeitschritt die An-

zahl der Zellen ausgegeben, fur die es zu unerwartet großen Unterschieden in den


Stromungswerten gekommen ist. Einzelne Zellen sind hierbei kein großes Problem,

eine großere Anzahl deutet jedoch auf eine Unregelmaßigkeit hin, die von verfalschten

Ergebnissen bis zu einem Absturz der Rechnung fuhren kann.

Inaktive Spenderzellen

Falls die Bewegung pro Zeitschritt zu groß ist, bzw. die Zellen fur die Bewegung rela-

tiv zu klein, kommt es zu folgender Situation: Es gelingt dem Suchalgorithmus zwar,

geeignete Spenderzellen fur die Koppelung der Gitter zu finden, diese sind jedoch Teil

des inaktiven Losungsraumes. Da in diesem Fall keine aktuellen Stromungswerte fur

diese Zelle vorliegen, konnen die Gebiete nicht sinnvoll gekoppelt werden. Die Feh-

lermeldungen lauten not enough active neighbors, inactive donor cells sowie switch

from inactive to donor/active without interpolation stage. Einen Ausweg bietet ein

groberes, eventuell auch gleichmaßigeres Gitter, oder eine Reduzierung der Bewegung

zwischen den Zeitschritten.

Numerik-Versagen

NaN-Stromungswerte in Zellen, die nicht von einem ungultigen aktiven Gitter herruhren,

konnen wie bei allen bewegten Gittern haufig durch Veranderungen des Unterrela-

xationsparameter oder des Blending-Faktors des Impulsgleichungslosers vemieden

werden. Standardmaßig bei jeweils 0.7 fuhrt ein Blendingfaktor von 0 sowie ein Un-

terrelaxationsparameter von 0.5 fur die ersten 50 bis 200 Zeitschritte oft zu einer

hinreichenden Stabilisierung der Rechnung.

Kapitel 7

Ergebnisse der

Simulationsrechnungen

Kapitel 6 beschreibt die allgemeinen Erkenntnisse, die aus den Simulationen mit Co-

met im Rahmen dieser Diplomarbeit gewonnen wurden. In diesem Kapitel werden

nun die messbaren Ergebnisse dieser Rechnungen vorgestellt. Die Resultate sind in

zwei Abschnitte untergliedert: Am Anfang stehen Versuche zur Validierung der neu-

en Simulationstechnik: Es werden Vergleiche von OLG-Simulationsergebnissen mit

Daten aus herkommlichen Berechnungsverfahren durchgefuhrt. Die Ergebnisse der

Versuche mit Anbauteilen finden sich im zweiten Abschnitt des Kapitels.

7.1 Modellvalidierung

Jenseits der in Kapitel 4 genannten Konvergenztheorie fur FD- und FE-Verfahren ist

dem Autor keine Quelle einer formalen mathematischen Analyse der Schwarz-Ver-

fahren fur Finite Volumen bekannt. Das schließt auch deren Anwendungen ein, wie

das hier untersuchte Overlapping Grid Verfahren. Es gibt jedoch andere Methoden

zur Validierung.

Hadzic (2005) fuhrt in der Entwicklung des OLG-Verfahrens folgende Schritte zur

Validierung durch:

• Vergleich mit Referenzergebnissen aus Versuchen und numerischen Verfahren

ohne Overlapping Grid

• Vergleich der Ergebnisse bei verschieden großen Uberlappungsbereichen

88

KAPITEL 7. ERGEBNISSE DER SIMULATIONSRECHNUNGEN 89

• Systematische Gitterverfeinerungen in allen Versuchen fur Aussagen uber die

Konvergenzordnung des Verfahrens

Die Ergebnisse der OLG-Versuche weisen durchweg eine sehr gute Deckung zu den

Vergleichsfallen auf, vgl. Hadzic (2005), Kapitel 6. Die Geschwindigkeit der Konver-

genz gegen eine gitterunabhangige Losung ist von zweiter Ordnung. Dies entspricht

den Erwartungen an ein Verfahren, welches auf Approximationen zweiter Ordnung

basiert.

Diese Arbeit soll daruber hinaus die Tauglichkeit des OLG-Verfahrens fur VSP-

Stromungsberechnungen untersuchen. Hierzu wurden Vergleiche der Overlapping

Grid Technik mit den herkommlichen Verfahren (econ-Klebeflachen fur stationare

Gitter, sliding interfaces fur bewegte Rechnungen) in verschiedenen Formen durch-

gefuhrt:

• Stationare Flugelumstromungen: Ein Flugelblatt steht senkrecht im Wasser

und wird von vorne angestromt. Die Krafte, die auf den Flugel wirken, werden

berechnet und ausgegeben. Fur einen festen Anstellwinkel des Flugels betrach-

tet man Fx, Fy und Mz : Die Kraftkomponente in Stromungsrichtung X, die

Querkraft in Y-Richtung, sowie das Moment an der Drehachse des Flugels Z.

Neben einer visuellen Analyse der Graphen wird die relative mittlere Abwei-

chung zum Referenzfall ermittelt. Hierzu vergleicht man die uber die letzten

100 Zeitschritte gemittelten Kraftwerte der Rechnungen, die einen stationaren

Endzustand erreicht haben.

F =1

100

100∑

k=1

F (tN−k+1) (7.1)

Daruber hinaus wird die Entwicklung der relativen Abweichung der Krafte vom

gemittelten Endzustand des Referenzfalles uber die Zeit untersucht:

Frel(t) =F (t) − F ref

F ref

(7.2)

Auf diese Weise sind Vergleiche des Konvergenzverhaltens der Verfahren moglich.

• VSP mit einem Flugel : Aus dem von Voith standardmaßig verwendeten Userco-

ding erhalt man fur jeden Zeitschritt Schub (Kraft in X-Richtung), Querkraft

(in Y-Richtung) und Wirkungsgrad des Propellers. Aufgrund der Bewegung

stellt sich kein stationarer Zustand im eigentlichen Sinn ein. Die Regelmaßig-

keit der Bewegung macht jedoch den Vergleich kompletter Perioden moglich.


Die Rahmenbedingungen der Simulationen sind jeweils identisch: Ein VSP hat

einen Flugelkreisdurchmesser von 1600 mm und eine Flugelblattlange von 1200

mm. Er wird von vorne mit einer Geschwindigkeit von 6,477 m/s angestromt.

Der Radkorper rotiert mit einer Drehzahl von 2,2216666 Hz, bei einem kon-

stanten Zeitschritt von 0,00125031 s entspricht das einer Drehung von 1 pro

Zeitschritt. Es wurden mindestens funf komplette Umdrehungen gerechnet, das

entspricht 1800 Zeitschritten. Auf jeden Zeitschritt kommen funf außere Itera-

tionen des Systemlosers. Der Autor folgt hierbei den bei Voith fur solche Falle

ublichen Vorgaben.

Auch hier werden die mittleren Kraftwerte verglichen: Die letzten 720 Zeit-

schritte sind eine gute Basis fur diese Mittelung, denn hier hat sich in der

Regel ein periodischer Zustand der Rechnung eingestellt.

F =1

720

720∑

k=1

F (tN−k+1) (7.3)

Zusatzlich betrachtet man die Entwicklung der relativen Abweichung der Krafte

vom Referenzfall mit Sliding Interfaces. Die Abweichung wird mit dem Maxi-

malwert des Sliding-Falles der letzten 3 Perioden normiert:

Frel(t) =F (t) − Fslid(t)

Fmaxslid

(7.4)

• VSP mit zwei und mehr Flugeln: Neben den Daten des ersten Flugels stehen

fur Mehrflugler auch die uber alle Flugel summierten Krafte fur Vergleiche zur

Verfugung.

Fur jede dieser Simulationen wurden neben den Kraften weitere Punkte analysiert:

• Plausibilitat des Druckbildes

• Plausibilitat der Zellfunktionen im Overlap-Bereich

• Erweiterte Fehlermeldungen

Aufgrund der Komplexitat der Falle und des hohen Rechenaufwandes der Simulatio-

nen wurde im Rahmen dieser Arbeit auf systematische Gitterverfeinerungen verzich-

tet. Eine solche Analyse hatte eine Neuvernetzung einiger Modellteile im Hinblick auf

eine gleichmaßigere Außenschicht notwendig gemacht, dies wurde aus Zeitgrunden

nicht mehr durchgefuhrt. Die Ergebnisse der Simulationen sind folglich nicht als


gitterunabhangige Losungen zu sehen. Allerdings wurde bereits fur die gegebenen

Gitter eine sehr gute Ubereinstimmung zwischen den Overlap-Rechnungen und dem

jeweiligen Referenzfall erzielt.

Die folgenden Abschnitte fassen die Ergebnisse der durchgefuhrten Simulationen zu-

sammen.

7.1.1 Flugelumstromung mit ursprunglicher Vernetzung

Tabelle 7.1: Mittlere Krafte einer stationare Flugelumstromung mit 12 m/s, Anstell-

winkel 4

Modell F x(N) F y(N) M z(Nm) ∆F x(%) ∆F y(%) ∆M z(%)

Slid 606,46 -9130,20 371,72 - - -

Econ 606,46 -9130,19 371,72 0 0 0

Over 604,91 -9054,25 369,06 -0,3 -0,8 -0,7

OverRed 602,49 -9056,34 370,30 -0,7 -0,8 -0,4

Ein stationarer Flugel mit Flugelblattlange 1200 mm wird mit 12 m/s in X-Richtung

angestromt. Die Krafte an dem mit 4 angestellten Flugelblatt werden mit vier Me-

thoden berechnet und miteinander verglichen: Ein mit Slid bezeichneter Fall mit Sli-

ding Interfaces und der mit Econ betitelte Fall einer festen Verbindung der Teilgitter

(Econ Regi) dienen als Referenzfalle nach herkommlichen Verfahren. Als Overlap-

Korper dienen der komplette Flugelzylinder und eine verschlankte Version. Abbil-

dungen A.1 bis A.4 zeigen die verwendete Vernetzung.

Berechnet wurden 2000 Zeitschritte von jeweils 0,001 Sekunden bei einer außeren

Iteration pro Zeitschritt. Graphen der berechneten Krafte siehe Abb. B.1.

Tabelle 7.1 zeigt die Ergebnisse des Vergleichs der gemittelten Krafte. Den relativen

Kraftevergleich siehe Abb. B.3.

Econ und Sliding-Fall weisen nach etwa 100 Zeitschritten identische Krafteverlaufe

auf und werden als gemeinsamer Referenzfall betrachtet. Die Overlap-Falle nahern

sich innerhalb von 40 Zeitschritten den Kraften der Referenzfalle an und verlaufen

weitgehend gemeinsam. Nach etwa 1000 Zeitschritten ist der stationare Zustand er-

reicht. Die Abweichung der Falle ist kleiner 1% fur alle Kraftkomponenten und bleibt

konstant. Der verbleibende Unterschied zwischen Sliding und Overlap ist auf unter-

schiedliche Feinheiten der verwendeten Hintergrundgitter zuruckzufuhren. Mit einer


Differenz von weniger als 0,4% zwischen den Overlap-Gittern gilt die Moglichkeit der

Verschlankung des Overlap-Korpers als verifiziert.

7.1.2 Flugelumstromung mit modifizierter Vernetzung

Tabelle 7.2: Mittlere Krafte einer stationaren Flugelumstromung mit 6,477 m/s, An-

stellwinkel 4

Modell F x(N) F y(N) M z(Nm) ∆F x(%) ∆F y(%) ∆M z(%)

Slid 203,061 -2706,14 -59,3429 - - -

Over 201,853 -2716,2 -57,9988 -0,6 0,4 -2,3

OverRef 202,058 -2717,75 -58,1963 -0,1 0,4 -1,9

OverAlt 185,647 -2620,78 -40,9884 -8,7 -3,2 -31

Over3600 202,607 -2713,17 -58,5369 -0,2 0,3 -1,4

Das Overlap-Flugelgitter wurde bestmoglich modifiziert, um eine gleichmaßigere Zell-

große am Overlap-Rand zu gewahrleisten. Hierzu wurde das Netz an der Hinterkante

des Flugels vergrobert sowie die Vernetzung der Unterseite des Flugels geandert, um

dort eine gunstigere Zellstruktur fur die Koppelung zu erhalten. Die Abbildungen A.7

bis A.9 zeigen die Veranderungen in der Vernetzung, die auch den Flugelzylinder fur

die Referenz-Sliding-Rechnungen betreffen.

Die Graphen in Abb. B.4 bis B.6 zeigen die Krafte, die auf einen stationaren Flugel

bei einer Anstromgeschwindigkeit vom 6,477 m/s und 4 Anstellwinkel wirken. Es

wurden 2000 Zeitschritte von jeweils 0,001 Sekunden bei einer außeren Iteration pro

Zeitschritt berechnet. Tabelle 7.2 zeigt einen Vergleich der mittleren Krafte. Der

Fall des modifizierten Overlap-Gitters ist mit Over bezeichnet. OverRef ist eine fur

stationare Rechnungen mogliche verbesserte Anpassung des Hintergrundgitters an

das Flugelgitter. Hierzu wurden die noch vorhandenen Zellgroßenunterschiede an der

Hinterkante des Flugels durch eine weitere cref-Verfeinerung des Hintergrundgitters

reduziert. Zum Vergleich sind die Referenzwerte Slid der Sliding-Rechnung, sowie

die Werte des ursprunglichen, verschlankten Overlap-Flugels OverAlt angegeben.

Die verbesserte Zellgroßenanpassung im Fall OverRef zeigt anfangs eine deutlich

schnellere Konvergenz gegen den Referenzwert als die ubrigen Falle. Dies ist ein

Indiz dafur, dass gleichmaßige Zellgroßen am Ubergang von großer Wichtigkeit fur

die Konvergenzgeschwindigkeit sind.


Das ursprungliche Gitter OverAlt weist große Zellgroßendifferenzen zum Hinter-

grundgitter und eine feinere Vernetzung des Flugelprofils auf. Die Krafte auf den

Flugel konvergieren fur dieses Gitter nicht gegen den neuen Referenzfall, was darauf

hindeutet, dass die gewahlte Vernetzung des Sliding-Falles fur eine gitterunabhangige

Losung zu grob ist. Aus oben genannten Grunden wird jedoch auf eine Untersuchung

von feineren Flugelvernetzungen verzichtet.

Die verbleibende Abweichung der ubrigen Overlap-Falle zum Sliding-Fall von bis zu

2% sind auf einen nach 2000 Zeitschritten noch nicht endgultig erreichten stationaren

Zustand zuruckzufuhren. Die Feinheit der verwendeten Gitter verzogert die Konver-

genz. Der Vergleichswert Over3600 nach 3600 Zeitschritten weist deutlich kleinere

Fehler auf.

7.1.3 VSP Einflugler mit ursprunglicher Vernetzung

Tabelle 7.3: Mittlere Krafte VSP (ein Flugel, ursprungliche Vernetzung)

Modell F x(N) F y(N) W.Grad ∆F x(%) ∆F y(%) ∆W.Grad (%)

Slid -16265,40 1153 0,414814 - - -

Over -16043,7 1087,39 0,416207 -1,4 -5,7 0,33

Basierend auf dem ursprunglichen, verschlankten Sliding-Flugelzylinder (vgl. Abb.

A.2) wurde ein VSP mit einem Flugel im OLG-Verfahren simuliert. Es wurden 1800

Zeitschritte bei 1 Drehung und funf außeren Iterationen pro Zeitschritt berechnet.

Abbildung A.6 zeigt die verwendeten Hintergrundnetze fur den Overlap- und den

auf Sliding Interfaces basierenden Referenzfall. Graphen der berechneten Krafte sie-

he in Abb. B.7 bis B.8. Trotz einer weitgehend guten Ubereinstimmung zeigen sich

hier noch Unterschiede, die im relativen Vergleich in Abb. B.13(Graph OverAlt)

deutlicher zu sehen sind als in den gemittelten Kraften in Tabelle 7.3. Die relati-

ven Abweichungen von bis zu 4% gehen auf Schwierigkeiten bei der Koppelung der

sehr fein vernetzten Flugelhinter- und Flugelunterkante mit dem Hintergrundgitter

zuruck, wie weitere Experimente zeigen.


Tabelle 7.4: Mittlere Krafte VSP (ein Flugel, modifizierte Vernetzung)


Slid -16265,40 1153 0,414814 - - -

Over -16186,9 1228,11 0,412221 0,5 6,5 -0,6

7.1.4 VSP Einflugler mit modifizierter Vernetzung

Fur den VSP mit einem modifizierten Flugel gelingt eine deutlich bessere Uberein-

stimmung mit dem Sliding Interfaces Referenzfall. Dies geht allerdings nicht aus den

gemittelten Kraften in Tabelle 7.4 hervor. Eine Betrachtung des Krafteverlaufs in

den Abbildungen B.9 und B.10 und des relativen Fehlers in Abbildung B.13 zeigt

jedoch die bessere Deckung der beiden Falle. Die maximale Abweichung kann von bis

zu 4% im ursprunglichen Fall auf etwa 1% reduziert werden. Wie ublich wurden 1800

Zeitschritte mit funf außeren Iterationen bei 1 Drehung pro Zeitschritt berechnet.

7.1.5 VSP Einflugler mit Vernetzung fur SP-Rechnungen

Tabelle 7.5: Mittlere Krafte VSP (ein Flugel fur SP-Rechnungen)


Slid -16265,40 1153 0,414814 - - -

OverSP -16540,20 1222,99 0,419111 1,7 6,1 1,03

OverOpt -16186,9 1228,11 0,412221 0,5 6,5 -0,6

Fur die Rechnungen mit Schutzplatte musste das Flugelgitter an der Unterkante

deutlich gestaucht werden, Abbildung A.10 zeigt die Modifikationen. Das Hinter-

grundgitter wurde in der Feinheit der Schichten entsprechend angepasst. Die Gra-

phen in Abb. B.11 und B.12 zeigen die Krafte auf einen VSP-Einflugler mit gekurz-

tem Overlap-Flugelgitter sowie den Sliding Interfaces Referenzfall (ohne Modifikati-

on bezuglich Schutzplattenrechnung). Wie immer wurden 1800 Zeitschritte mit funf

außeren Iterationen bei 1 Drehung pro Zeitschritt berechnet.

Es dauert fast 360 Zeitschritte, bis eine Ubereinstimmung mit dem Referenzfall ein-

tritt. Auch der relative Fehler ist im Vergleich zum optimalen, nicht gekurzten Gitter


großer geworden, wie Abbildung B.13 und Tabelle 7.5 zeigen. Allerdings ist die De-

ckung noch deutlich besser als im nicht modifizierten Fall, die Abweichung betragt

meist weniger als 2%.

7.1.6 VSP Funfflugler

Tabelle 7.6: Mittlere Krafte VSP (funf Flugel fur SP-Rechnungen)


Slid1 -10614 -883,566 0,314003 - - -

Over1 -10741 -782,618 0,314361 1,2 -11,4 0,11

OverOpt1 -10640,7 -766,472 0,313851 0,25 -13,25 -0,05

Slid5 -53073,7 -4412,43 0,625686 - - -

Over5 -53683 -3882,86 0,625319 1,1 -12,0 0,06

OverOpt5 -53177,2 -3807,81 0,625228 0,20 -13,7 -0,07

Mit dem gekurzten Flugelgitter wurde ein VSP-Funfflugler gerechnet und mit einem

Sliding Interface Referenzfall verglichen. Einen Uberblick des Aufbaus gibt Abb.

A.11. Zusatzlich wurde ein Funfflugler mit den optimalen, ungekurzten Flugelgittern

simuliert. Wie immer wurden 1800 Zeitschritte zu je funf außeren Iterationen bei 1

Drehung pro Zeitschritt berechnet. Abb. B.14 bis B.18 zeigen die Krafte, Tabelle

7.6 die gemittelten Kraftwerte. Es wurden jeweils die Krafte fur den ersten der funf

Flugel sowie die aufsummierten Werte untersucht.

Nach etwa 360 Zeitschritten tritt Ubereinstimmung mit dem Referenzfall auf, die

Deckung der Kurven ist weitgehend gut. Der Unterschied der beiden Overlap-Falle

ist minimal, weshalb abgesehen von den relativen Fehlerplots in Abb. B.16 und B.18

auf Darstellung des Falles verzichtet wurde.

Die großen Abweichungen der mittleren Fehler der Querkraft sind bezuglich der

Großenordnung der Krafte zu sehen. Der relative Fehler der Kraft uber einen Flugel

liegt unter 4%. Es wird vermutet, dass die großen Abweichungen durch eine feinere

Grenzschichtauflosung behoben werden konnen, die eine genauere Bestimmung der

Richtung der Kraftkomponenten zulasst. Dass die Deckung der Wirkungsgrade mit

einem Fehler von kleiner 1% sehr gut ist scheint diese Theorie zu unterstutzen.


7.1.7 Laufzeiten VSP-Rechnungen (Vernetzung fur SP)

Tabelle 7.7: Laufzeiten VSP fur SP-Rechnungen, 1500 Zeitschritte

Modell Zellen ges. Flugelzellen Node1 Node2 CPU-Zeit Gesamtzeit

Over 913812 74240 100% - 1,212E+5 1,218E+5

OverPar 913812 74240 66% 34% 8,783E+4 8,792E+4

Slid 385062 75180 100% - 3,786E+4 4,272E+4

Over5 1210772 371200 100% - 1,906E+5 1,910E+5

Slid5 765630 375900 100% - 8,373E+4 8,448E+4

Auf formale Analysen des Geschwindigkeitsgewinnes durch Parallelisierung wurde

in dieser Arbeit wegen der herrschenden Beschrankung auf zwei Knoten verzichtet.

Auch Vergleiche der Laufzeit von Overlap und Sliding Interface Fallen erschienen

dem Autor aufgrund der durch die Modellierung bedingten, großen Unterschiede der

Zellenzahl nicht sinnvoll.

Aus Grunden der Vollstandigkeit seien an dieser Stelle in Tabelle 7.7 einige Simula-

tionslaufzeiten angegeben. Neben Ein- und Funfflugler im Sliding und Overlap Fall

wurde die Rechnung des VSP Einfluglers parallelisiert: Hierzu wurde eine Zellkno-

ten-Schicht in Vorder- und Hintergrundgitter auf einheitliche Hohe gebracht (vmod

vset), vgl. Abschnitt 6.3 und Abb. A.12. Die Gittermodifikation ist hierbei mini-

mal, da eine Schicht mit sehr geringen Differenzen der Knotenhohe in Vorder- und

Hintergrund ausgewahlt wurde, die Ergebnisse zeigen keinerlei Abweichung.

7.2 Rechnungen mit Anbauteilen

Als zweiter Kernpunkt der Arbeit sollte die Moglichkeiten ausgelotet werden, Anbau-

teile wie Streben und Schutzplatte in die Modellierung aufzunehmen. Die Ergebnisse

aus erfolgreich durchgefuhrten Simulationen sollen an dieser Stelle zusammengefasst

werden. Aus genannten Grunden gibt es in diesem Bereich allerdings nur wenige

Referenzdaten fur formale Kraftevergleiche.


7.2.1 VSP mit Streben

Abbildung A.13 zeigt die verwendete Vernetzung fur ersten Versuche mit Streben, die

allerdings nicht zu messbaren Ergebnissen fuhrten, da zu Anfang nur kurze Flugel-

bewegungen an einer Strebe vorbei berechnet wurden.

7.2.2 VSP mit Schutzplatte

Abbildung A.14 zeigt die Vernetzung der vereinfachten Schutzplatte, die sowohl als

Overlap-Korper wie auch als Aussparung im Hintergrundgitter mit einem Einflugler

simuliert wurde. Der Abstand der Flugelspitze zur Platte betragt jeweils 1% des

Flugelkreisdurchmessers. Die Krafte der Falle sind in Abbildung C.1 gegeben. Zum

Vergleich ist ein auf Sliding Interfaces basierter Referenzfall mit Schutzplatte ange-

geben.

Die Deckung der Flugelkrafte der beiden simulierten Falle ist weitgehend gut, der

SP-Overlap-Fall zeigt einen etwas weniger glatten Kurvenverlauf. Die Abweichung

zu den Sliding-Werten ist durch die unterschiedliche Gitterstruktur des Falles zu

erklaren, der nicht auf dem Standard-Sliding-Referenzgitter basiert.

Die Krafte, die im Verlauf der Rechnung auf die Schutzplatte wirken, finden sich

in Abbildung C.2. Die Modellierung auf Basis der Aussparung im Hintergrundgitter

weist hier einen deutlich glatteren Krafteverlauf auf als der SP-Overlap-Fall. Der

Sliding-Fall basiert auf einer anderen Schutzplatte, was die deutliche Abweichung

erklart, er ist nur zum Vergleich der Periodizitat der Krafte angegeben.

Durch Modifikation der gegebenen Gitter ist es nicht gelungen, den Krafteverlauf

des SP-Overlap-Falles zu glatten. Im weiteren Verlauf der Arbeit wurde auf die Mo-

dellierung der Schutzplatte als Overlap-Korper verzichtet.

7.2.3 VSP mit Streben und Schutzplatte

Aus der Kombination von vereinfachter Strebe und Schutzplatte wie in Abb. A.15

ergaben sich keine neuen Vergleichswerte. Die Verbindung der beiden Falle war aller-

dings ein wichtiger Zwischenschritt zu einer realistisch vernetzten Schutzplatte mit

Streben wie in Abb. A.17.

Die hierfur verwendete Schutzplattenvernetzung basiert auf einem Gitter mit sehr

ungleichmaßiger Oberflache, wie Abb. A.16 verdeutlicht. Dies fuhrt zu Problemen


bei der Koppelung, wie das Bild der Zellfunktionszuordnung besonders fur den Funf-

flugler zeigt.

Fur den Einflugler machen sich diese Probleme nur schwach bemerkbar, wie der je-

weilige Krafteverlauf zeigt: Der in Abb. C.3 mit OneFull bezeichnete Einflugler mit

kompletten Anbauteilen zeigt einen glatten Verlauf in Schub und Querkraft, ahnlich

der Kurve des Sliding-Einfluglers ohne Streben. Im Verlauf der Schutzplattenkrafte in

Abb. C.6 weisen allerdings deutlich großere Unsauberkeiten als in den Vergleichsfalle

auf die Probleme der Koppelung hin. Die Strebenkrafte in Abb. C.5 weisen eine gute

Periodizitat ohne unerwartete Sprunge auf. Der Verlauf tragt auch den unterschied-

lichen Anstellwinkeln des Flugels beim Passieren der jeweiligen Strebe Rechnung.

Auch fur die Erweiterung auf funf Flugel steht nur ein Sliding-Referenzfall mit leicht

abgewandeltem Gitter zu Verfugung. Abweichungen zu diesem Fall, der daruber hin-

aus keine Streben enthalt, sind unter diesem Aspekt zu sehen. Der Verlauf der Krafte

des ersten Flugels, siehe Abb. C.3, weist kleine Schwankungen auf. Uber funf Flugel

aufsummiert verstarken sie sich deutlich, wie Abb. C.4 zeigt. Der zu Vergleichszwe-

cken phasenverschoben aufsummierte Wert des ersten Flugels zeigt, dass es sich nicht

um ein Problem der Krafteberechnung des Usercodings handelt.

Verlauf und Großenordnung der Krafte ahneln dem Funfflugler ohne Schutzplatte,

vgl. Abb. B.17: Diese Falle weisen Abweichungen in Schub und Querkraft zum Re-

ferenzfall auf, die sich im Wirkungsgrad wieder ausgleichen, so dass hier eine gute

Deckung erzielt wird. Fur den Wirkungsgrad steht fur SP-Rechnungen kein Referenz-

fall zu Verfugung, zum Vergleich ist der Wert des Sliding-Falles ohne SP angegeben.

Beide Kurven weisen einen sehr glatten Verlauf auf. Wie bereits im Fall ohne An-

bauteile erwahnt, wird vermutet, dass eine feinere Profilvernetzung die Abweichung

in den X- und Y-Komponenten der Kraft deutlich reduzieren kann. Großere Aus-

schlage an einigen Stellen im Krafteverlauf, die nicht durch Interaktion mit Streben

zu erklaren sind, sind auf Probleme bei der Koppelung zuruckzufuhren. Sie sollten

durch eine gleichmaßigere SP-Vernetzung leicht zu beheben sein.

Die Krafte auf die Streben in Abb. C.7 sind plausibel, sie zeigen abgesehen von eini-

gen Spitzen ein periodisches Bild. Die Krafte auf die Schutzplatte hingegen verlaufen

sehr ungleichmaßig, die Probleme der Koppelung sind hier am offensichtlichsten, vgl.

Abb. C.8. Die Auswirkungen sind teilweise auch in der Druckverteilung zu sehen, wie

Abb. A.18 zeigt.

Kapitel 8

Fazit und Ausblick

CFD-Berechnungen mit uberlappenden Gittern in Bezug auf den Voith-Schneider-

Propeller : Ziel dieser Diplomarbeit war es, ein Overlapping Grid Verfahren fur Pro-

bleme der Stromungsmechanik zu validieren, wie sie im Rahmen der Weiterentwick-

lung des VSP auftreten. Hierfur wurde eine bestehende Implementierung des Verfah-

rens in Comet untersucht, einem Finite-Volumen-CFD-Loser der Firma CD-adapco.

Hier soll zum Abschluss eine kurze Zusammenfassung der Arbeitsergebnisse gegeben

werden.

Die Erkenntnisse, die sich im Rahmen dieser Arbeit aus dem Studium der theoreti-

schen Grundlagen ergaben, lassen sich zu folgenden Punkten zusammenfassen:

Der CFD-Loser basiert auf einer RANSE-Diskretisierung durch ein FV-Verfahren

2. Ordnung mit vollimpliziter Drei-Ebenen-Zeitdiskretisierung. Es wird das k − ε -

Turbulenzmodell zum Abschluss der Systemgleichungen und der ALE-Ansatz zur

Modellierung bewegter Gitter verwendet.

Das Overlapping Grid Verfahren ist ein Algorithmus, der die Losung der Teilgitter

durch lineare Interpolation miteinander koppelt. Die hierdurch auftretenden Mas-

senungleichgewichte werden durch den SIMPLE-Algorithmus korrigiert. Ein Such-

verfahren ubernimmt die automatische Bestimmung der jeweiligen Interpolations-

partner zwischen den Gittern.

Das OLG-Verfahren ist kein Gebietszerlegungsverfahren im eigentlichen Sinne:

Wahrend es prinzipiell eng mit der Klasse der Schwarz-Verfahren verwandt ist, ver-

zichtet man zu Gunsten einer starken Koppelung der Teillosungen auf die Moglichkeit

der parallelisierten Losungsberechnung auf den Teilgebieten. Das gekoppelte System

99

KAPITEL 8. FAZIT UND AUSBLICK 100

wird in einer gemeinsamen Matrix iterativ durch ein CG-Verfahren fur nicht sym-

metrische Matrizen mit IC-Vorkonditionierung gelost.

Die fur Schwarz-Verfahren in FD- und FEM-Diskretisierung existierenden konver-

genztheoretischen Ergebnisse lassen sich nicht auf FV-Verfahren ubertragen. Die

Modellvalidierung muss auf andere Weise geschehen: In solchen Fallen werden ubli-

cherweise direkte Vergleiche von Simulationsergebnissen in Kombination mit syste-

matischen Gitterverfeinerungen durchgefuhrt.

Die praktischen Ergebnisse aus der Anwendung von Comet lauten wie folgt:

Die Modellierung der Bewegung des VSP konnte in das neue Verfahren ubertragen

werden. Daruber hinaus wurden viele Erkenntnisse bezuglich der Realisierung von

Anbauteilen im direkten Umfeld des Antriebs gewonnen. Sie wurden erfolgreich in

der Simulation eines VSP Funffluglers mit Schutzplatte und drei Streben umgesetzt.

Die Erfahrungen aus Versuchen, die im Rahmen dieser Arbeit durchgefuhrt wurden,

sind in einem ausfuhrlichen Katalog von Richtlinien fur OLG-Simulationen zusam-

mengefasst.

Die zur Modellvalidierung durchgefuhrten Vergleiche umfassen die Analyse von sta-

tionaren Flugelumstromungen und VSP mit einem und funf Flugeln. Neben einer

direkten Analyse des Krafteverlaufs wurden gemittelte Krafte verglichen, sowie die

relative Abweichung zu geeigneten Referenzfallen bestimmt. Da bis auf wenige Aus-

nahmen Abweichungen von unter 2% bei gleicher Konvergenzgeschwindigkeit erzielt

wurden, kann man im Hinblick auf die verwendeten Gitter von einer sehr guten

Ubereinstimmung der Ergebnisse sprechen. Die Flugelvernetzung, die ursprunglich

aus Sliding Interface Fallen stammt, wurde bestmoglich an die OLG-Kriterien ange-

passt, um eine gute Koppelung zu gewahrleisten. Durch eine im Hinblick auf OLG

optimierte, nach außen gleichmaßigere Neuvernetzung ließe sich die Abweichung zu

den Sliding-Referenzfallen noch weiter reduzieren.

Eine uberraschende Erkenntnis war die durch Comet eingeschrankte Parallelisier-

barkeit: In der vorliegenden Version 2.360.Z konnten Falle, die ein Overlap-Gebiet

enthalten, auf maximal zwei Knoten parallelisiert werden. Die Moglichkeit einer deut-

lich erweiterten Parallelisierbarkeit eines FV-OLG-Verfahrens haben Madrane u. a.

(2004) am Beispiel des DLR-Tau-Codes nachgewiesen.

Folgende Punkte, die im Rahmen dieser Arbeit untersucht wurden, konnten nicht

abschließend geklart werden:


Beim Vergleich zwischen Sliding- und Overlap-Modellierung eines VSP kommt es in

den Querkraften besonders fur Funfflugler zu großeren Abweichungen, wahrend die

Deckung des Wirkungsgrades mit einem relativen Fehler von kleiner 1% sehr gut

bleibt. Es wird vermutet, dass sich diese Abweichung durch eine feinere Vernetzung

des Flugelprofils reduzieren lasst.

Aufgrund der Komplexitat der verwendeten Falle und der Struktur der verwende-

ten Gitter konnten keine systematischen Verfeinerungen zum Erreichen einer git-

terunabhangig Losung durchgefuhrt werden. Hierfur ware eine Neuvernetzung der

Flugelgitter bezuglich einer gleichmaßigeren Außenschicht notwendig gewesen.

Bei der Simulation des VSP mit einer realistischen Schutzplatte kam es aufgrund

der Ungleichmaßigkeiten der Gitterstruktur zu Schwierigkeiten der Koppelung. Eine

Neuvernetzung der SP mit gleichmaßigerer Außenschicht sollte diese Probleme jedoch

beheben. Im Fall einer vereinfachten, regelmaßig vernetzten Schutzplatte traten diese

Probleme nicht auf.

Daruber hinaus war es aus Zeitgrunden nicht mehr moglich, einen luckenlosen Uber-

gang der Streben in die Schutzplatte zu modellieren. Es verblieben einige Zellschich-

ten im Zwischenraum, deren Elimination eine Neuvernetzung des entsprechenden

Bereichs erfordert hatten.

Zusammenfassend lasst sich sagen, dass sich der Betrieb des VSP sehr gut mit dem

OLG-Verfahren simulieren lasst. Die Moglichkeiten der Modellierung von Anbautei-

len im Umfeld des Antriebs werden durch die neue Methode deutlich erweitert: Die

Kombinationen von Streben und VSP ist nun moglich.

Die Validierung des Verfahrens konnte im Rahmen dieser Arbeit allerdings noch nicht

abgeschlossen werden: Fur die in dieser Arbeit verwendeten Gitter sind gute Ergeb-

nisse erzielt worden. Damit die Methode allerdings als endgultig verifiziert gelten

kann, mussen weitere Untersuchungen mit feiner vernetzten Gittern durchgefuhrt

werden. Hierzu konnen die umfassenden Erkenntnisse, die im Rahmen dieser Arbeit

uber das OLG-Verfahren gewonnen und dokumentiert wurden, einen wichtigen Bei-

trag leisten.

Eine großes Problem fur die Anwendung des Verfahrens ist die starke Einschrankung

der Parallelisierbarkeit von Overlap-Fallen, die von Comet auf zwei Knoten be-

schrankt wird.


Ausblick

Das OLG-Verfahren hat sich zu einer vielseitigen Methode zur Behandlung von

komplexen Bewegungen in Prozessen der numerischen Stromungsmechanik entwi-

ckelt, wie auch die Ergebnisse dieser Arbeit zeigen. Forschung zur weiteren Verbes-

serung dieser Methodik wird momentan in unterschiedlichen Richungen betrieben:

Ein Punkt ist die adaptive, dynamische Gitterverfeinerung im Umfeld der Overlap-

Korper: Die Modellqualitat kann so deutlich verbessert werden, ohne die Zellenzahl

unnotig zu erhohen.

Daruber hinaus ist die Parallelisierbarkeit von OLG-Rechnungen im Hinblick auf

die Großenentwicklung der Simulationen von entscheidender Wichtigkeit. Henshaw

(2008) prasentiert neueste Ergebnisse aus beiden Bereichen.

Anhang A

Abbildungen zur Vernetzung

Abbildung A.1: Flugel mit ersten Vernetzungsschichten (links), Druckverteilung auf

Flugelblatt (rechts)

103

ANHANG A. ABBILDUNGEN ZUR VERNETZUNG 104

Abbildung A.2: Ursprunglicher Flugelzylinder (links), 166800 Zellen, reduzierter

Flugelzylinder (rechts), 130000 Zellen

Abbildung A.3: Reduziertes Flugelgitter, Nahaufnahme


Abbildung A.4: Schnitt Hintergrundgitter Sliding/Econ (oben), 423400 Zellen, Hin-

tergrundgitter Overlap (unten), 174528 Zellen, mit jeweiligem Flugelgitter


Abbildung A.5: Overlap-Zellfunktionen (IOVC, vgl. Tabelle 6.1), Schnitt Hinter-

grundgitter (oben), vollstandiges Netz (unten)


Abbildung A.6: VSP Einflugler: Sliding-Fall (oben), 309882 Zellen Hintergrundgitter,

Overlap-Fall (unten), 435520 Hintergrundzellen mit jeweiligem Flugel


Abbildung A.7: Fur Overlap modifizierter Flugel (links), 102800 Zellen, und entspre-

chender Zylinder fur Sliding-Referenzfalle (rechts), 136100 Zellen

Abbildung A.8: Modifiziertes Flugelgitter, Nahaufnahme Oberseite


Abbildung A.9: Unterseite ursprungliches Flugelgitter (oben) und modifiziertes

Flugelgitter (unten)


Abbildung A.10: Gestauchtes Flugelgitter fur Rechnungen mit Schutzplatte, 74200

Zellen, 8 Zellschichten zur Flugelspitze


Abbildung A.11: Funfflugler Sliding (765630 Zellen, davon 230340 Hintergrund)

(oben), und Overlap (1210772 Zellen, davon 839572 Hintergrund) (unten)


Abbildung A.12: Parallelisierung: Erfolgreiche Partitionierung auf 2 Knoten (links),

Probleme bei 3 Knoten (rechts)

Abbildung A.13: Vereinfachte Streben: Overlapkorper (o.l.), IOVC Hintergrund

(o.r.), Strebe als Aussparung (u.l.), Kontakt Overlap-Rand Aussparung finale Strebe

(u.r.)


Abbildung A.14: Vereinfachte Schutzplatte: Druckverteilung (o.l.), IOVC Hinter-

grund auf Hohe Flugelspitze (o.r.), Gitter fur Overlap-Modellierung (u.l.), und Aus-

sparung (u.r.)


Abbildung A.15: Vereinfachte Schutzplatte mit Strebe: Gemeinsame Aussparung

(o.l.), aktive Zellen fur Strebe und SP als eigene Overlap-Korper (o.r.), aktive Zellen

fur Strebe als Overlap, SP als Aussparung (u.)


Abbildung A.16: Finale Schutzplatte: IOVC-Schnitt Hintergrund mit SP-Netz in

Blockform (o.), IOVC-Verteilung fur verfeinertes SP-Gitter auf Hohe Flugelspitze

mit Ungleichmaßigkeiten (m. und u.)


Abbildung A.17: VSP Einflugler mit finaler SP und 3 Streben: Oberflachennetz (o.)

und Druckverteilung (m. und u.)


Abbildung A.18: VSP Funfflugler mit finaler SP und 3 Streben: Druckverteilung (o.)

mit Unregelmaßigkeit (u.)

Anhang B

Graphen der Krafte fur VSP ohne

Anbauteile

118

ANHANG B. GRAPHEN DER KRAFTE FUR VSP OHNE ANBAUTEILE 119

Abbildung B.1: Stationare Flugelumstromung mit 12 m/s, Anstellwinkel 4 1/2

-2000

-1500

-1000

-500

0

500

1000

1500

2000

20 40 60 80 100 120 140 160 180 200

Fx

bei 4

Gra

d, 1

2 m

/s

Zeitschritt

EconSlid

OverOverRed

-9000

-8000

-7000

-6000

-5000

-4000

-3000

20 40 60 80 100 120 140 160 180 200

Fy

bei 4

Gra

d, 1

2 m

/s

Zeitschritt

EconSlid

OverOverRed

-300

-200

-100

0

100

200

300

400

500

600

20 40 60 80 100 120 140 160 180 200

Mz

bei 4

Gra

d, 1

2 m

/s

Zeitschritt

EconSlid

OverOverRed


Abbildung B.2: Stationare Flugelumstromung mit 12 m/s, Anstellwinkel 4 2/2

530

540

550

560

570

580

590

600

610

200 400 600 800 1000 1200 1400 1600 1800 2000

Fx

bei 4

Gra

d, 1

2 m

/s

Zeitschritt

EconSlid

OverOverRed

-9200

-9100

-9000

-8900

-8800

-8700

-8600

-8500

-8400

-8300

-8200

200 400 600 800 1000 1200 1400 1600 1800 2000

Fy

bei 4

Gra

d, 1

2 m

/s

Zeitschritt

EconSlid

OverOverRed

354

356

358

360

362

364

366

368

370

372

374

200 400 600 800 1000 1200 1400 1600 1800 2000

Mz

bei 4

Gra

d, 1

2 m

/s

Zeitschritt

EconSlid

OverOverRed


Abbildung B.3: Rel. Abweichung Stat. Umstromung mit 12 m/s, Anstellwinkel 4

-0.12

-0.1

-0.08

-0.06

-0.04

-0.02

0

0.02

600 800 1000 1200 1400 1600 1800 2000

Rel

. Feh

ler

Fx

bei 4

Gra

d, 1

2 m

/s

Zeitschritt

SlidEconOver

OverRed

-0.02

-0.015

-0.01

-0.005

0

0.005

600 800 1000 1200 1400 1600 1800 2000

Rel

. Feh

ler

Fy

bei 4

Gra

d, 1

2 m

/s

Zeitschritt

SlidEconOver

OverRed

-0.014

-0.012

-0.01

-0.008

-0.006

-0.004

-0.002

0

0.002

600 800 1000 1200 1400 1600 1800 2000

Rel

. Feh

ler

Mz

bei 4

Gra

d, 1

2 m

/s

Zeitschritt

SlidEconOver

OverRed


Abbildung B.4: Umstromung mod. Fl.Gitter mit 6,477 m/s, Anstellwinkel 4 1/2

140

160

180

200

220

240

260

280

300

320

340

20 40 60 80 100 120 140 160 180 200

Fx

bei 4

Gra

d, 6

,477

m/s

Zeitschritt

SlidOverAlt

OverOverRef

-3500

-3000

-2500

-2000

-1500

-1000

20 40 60 80 100 120 140 160 180 200

Fy

bei 4

Gra

d, 6

,477

m/s

Zeitschritt

SlidOverAlt

OverOverRef

-140

-120

-100

-80

-60

-40

-20

0

20 40 60 80 100 120 140 160 180 200

Mz

bei 4

Gra

d, 6

,477

m/s

Zeitschritt

SlidOverAlt

OverOverRef


Abbildung B.5: Umstromung mod. Fl.Gitter mit 6,477 m/s, Anstellwinkel 4 2/2

165

170

175

180

185

190

195

200

205

200 400 600 800 1000 1200 1400 1600 1800 2000

Fx

bei 4

Gra

d, 6

,477

m/s

Zeitschritt

SlidOverAlt

OverOverRef

-2750

-2700

-2650

-2600

-2550

-2500

-2450

-2400

200 400 600 800 1000 1200 1400 1600 1800 2000

Fy

bei 4

Gra

d, 6

,477

m/s

Zeitschritt

SlidOverAlt

OverOverRef

-60

-55

-50

-45

-40

-35

-30

-25

200 400 600 800 1000 1200 1400 1600 1800 2000

Mz

bei 4

Gra

d, 6

,477

m/s

Zeitschritt

SlidOverAlt

OverOverRef


Abbildung B.6: Rel. Abweichung Umstromung mod. Gitter mit 6,477 m/s, Winkel

4

-0.12

-0.1

-0.08

-0.06

-0.04

-0.02

0

0.02

500 1000 1500 2000 2500

Rel

. Feh

ler

Fx

bei 4

Gra

d, 6

,477

m/s

Zeitschritt

SlidOver

OverRef

-0.01

-0.005

0

0.005

0.01

0.015

500 1000 1500 2000 2500

Rel

. Feh

ler

Fy

bei 4

Gra

d, 6

,477

m/s

Zeitschritt

SlidOver

OverRef

-0.3

-0.25

-0.2

-0.15

-0.1

-0.05

0

0.05

500 1000 1500 2000 2500

Rel

. Feh

ler

Mz

bei 4

Gra

d, 6

,477

m/s

Zeitschritt

SlidOver

OverRef


Abbildung B.7: VSP (ein ursprungliches Flugelgitter) 1/2

-25000

-20000

-15000

-10000

-5000

0

5000

10000

20 40 60 80 100 120 140 160 180 200

Sch

ub

Zeitschritt

SlidOver

-5000

0

5000

10000

15000

20000

25000

20 40 60 80 100 120 140 160 180 200

Que

rkra

ft

Zeitschritt

SlidOver

-1

-0.8

-0.6

-0.4

-0.2

0

0.2

0.4

0.6

0.8

20 40 60 80 100 120 140 160 180 200

Wirk

ungs

grad

Zeitschritt

SlidOver


Abbildung B.8: VSP (ein ursprungliches Flugelgitter) 2/2

-60000

-50000

-40000

-30000

-20000

-10000

0

10000

200 400 600 800 1000 1200 1400 1600 1800

Sch

ub

Zeitschritt

SlidOver

-40000

-30000

-20000

-10000

0

10000

20000

30000

200 400 600 800 1000 1200 1400 1600 1800

Que

rkra

ft

Zeitschritt

SlidOver

-1

-0.8

-0.6

-0.4

-0.2

0

0.2

0.4

0.6

0.8

1

200 400 600 800 1000 1200 1400 1600 1800

Wirk

ungs

grad

Zeitschritt

SlidOver


Abbildung B.9: VSP (ein modifiziertes Flugelgitter) 1/2

-25000

-20000

-15000

-10000

-5000

0

5000

10000

20 40 60 80 100 120 140 160 180 200

Sch

ub

Zeitschritt

SlidOver

-5000

0

5000

10000

15000

20000

25000

20 40 60 80 100 120 140 160 180 200

Que

rkra

ft

Zeitschritt

SlidOver

-1

-0.8

-0.6

-0.4

-0.2

0

0.2

0.4

0.6

0.8

20 40 60 80 100 120 140 160 180 200

Wirk

ungs

grad

Zeitschritt

SlidOver


Abbildung B.10: VSP (ein modifiziertes Flugelgitter) 2/2

-60000

-50000

-40000

-30000

-20000

-10000

0

10000

200 400 600 800 1000 1200 1400 1600 1800

Sch

ub

Zeitschritt

SlidOver

-40000

-30000

-20000

-10000

0

10000

20000

30000

200 400 600 800 1000 1200 1400 1600 1800

Que

rkra

ft

Zeitschritt

SlidOver

-1

-0.8

-0.6

-0.4

-0.2

0

0.2

0.4

0.6

0.8

1

200 400 600 800 1000 1200 1400 1600 1800

Wirk

ungs

grad

Zeitschritt

SlidOver


Abbildung B.11: VSP (ein Flugelgitter fur SP-Rechnungen) 1/2

-25000

-20000

-15000

-10000

-5000

0

5000

20 40 60 80 100 120 140 160 180 200

Sch

ub

Zeitschritt

SlidOver

-5000

0

5000

10000

15000

20000

25000

30000

20 40 60 80 100 120 140 160 180 200

Que

rkra

ft

Zeitschritt

SlidOver

-1

-0.8

-0.6

-0.4

-0.2

0

0.2

0.4

0.6

0.8

20 40 60 80 100 120 140 160 180 200

Wirk

ungs

grad

Zeitschritt

SlidOver


Abbildung B.12: VSP (ein Flugelgitter fur SP-Rechnungen) 2/2

-60000

-50000

-40000

-30000

-20000

-10000

0

10000

200 400 600 800 1000 1200 1400 1600 1800

Sch

ub

Zeitschritt

SlidOver

-40000

-30000

-20000

-10000

0

10000

20000

30000

200 400 600 800 1000 1200 1400 1600 1800

Que

rkra

ft

Zeitschritt

SlidOver

-1

-0.8

-0.6

-0.4

-0.2

0

0.2

0.4

0.6

0.8

1

200 400 600 800 1000 1200 1400 1600 1800

Wirk

ungs

grad

Zeitschritt

SlidOver


Abbildung B.13: Relative Abweichung Krafte VSP Einflugler

-0.025

-0.02

-0.015

-0.01

-0.005

0

0.005

0.01

0.015

0.02

0.025

0.03

600 800 1000 1200 1400 1600 1800

Rel

. Abw

eich

ung

Sch

ub

Zeitschritt

OverAltOverOptOverSP

-0.04

-0.03

-0.02

-0.01

0

0.01

0.02

0.03

0.04

600 800 1000 1200 1400 1600 1800

Rel

. Abw

eich

ung

Que

rkra

ft

Zeitschritt


-0.04

-0.02

0

0.02

0.04

600 800 1000 1200 1400 1600 1800

Rel

. Abw

eich

ung

Wirk

ungs

grad

Zeitschritt



Abbildung B.14: VSP Funfflugler, Krafte auf 1. Flugel 1/2

-40000

-30000

-20000

-10000

0

10000

20000

100 200 300 400 500 600 700

Sch

ub

Zeitschritt

SlidOver

-25000

-20000

-15000

-10000

-5000

0

5000

10000

15000

20000

25000

100 200 300 400 500 600 700

Que

rkra

ft

Zeitschritt

SlidOver

-1

-0.5

0

0.5

1

100 200 300 400 500 600 700

Wirk

ungs

grad

Zeitschritt

SlidOver


Abbildung B.15: VSP Funfflugler, Krafte auf 1. Flugel 2/2

-40000

-30000

-20000

-10000

0

10000

20000

800 1000 1200 1400 1600 1800

Sch

ub

Zeitschritt

SlidOver

-25000

-20000

-15000

-10000

-5000

0

5000

10000

15000

20000

25000

800 1000 1200 1400 1600 1800

Que

rkra

ft

Zeitschritt

SlidOver

-1

-0.8

-0.6

-0.4

-0.2

0

0.2

0.4

0.6

0.8

800 1000 1200 1400 1600 1800

Wirk

ungs

grad

Zeitschritt

SlidOver


Abbildung B.16: Relative Abweichung Funfflugler, Krafte 1. Flugel

-0.03

-0.02

-0.01

0

0.01

0.02

0.03

600 800 1000 1200 1400 1600 1800

Rel

. Feh

ler

Sch

ub

Zeitschritt

OverOptOverSP

-0.05

-0.04

-0.03

-0.02

-0.01

0

0.01

0.02

0.03

0.04

0.05

600 800 1000 1200 1400 1600 1800

Rel

. Feh

ler

Que

rkra

ft

Zeitschritt

OverOptOverSP

-0.04

-0.02

0

0.02

0.04

600 800 1000 1200 1400 1600 1800

Rel

. Feh

ler

Wirk

ungs

grad

Zeitschritt

OverOptOverSP


Abbildung B.17: VSP Funfflugler, summierte Krafte

-56000

-55000

-54000

-53000

-52000

-51000

-50000

800 1000 1200 1400 1600 1800

Sch

ub

Zeitschritt

SlidOver

-8000

-7000

-6000

-5000

-4000

-3000

-2000

-1000

0

800 1000 1200 1400 1600 1800

Que

rkra

ft

Zeitschritt

SlidOver

0.58

0.59

0.6

0.61

0.62

0.63

0.64

0.65

0.66

0.67

0.68

800 1000 1200 1400 1600 1800

Wirk

ungs

grad

Zeitschritt

SlidOver


Abbildung B.18: Relative Abweichung Funfflugler, summierte Krafte

-0.02

-0.015

-0.01

-0.005

0

0.005

0.01

0.015

0.02

0.025

0.03

1200 1300 1400 1500 1600 1700 1800

Rel

. Feh

ler

Sch

ub

Zeitschritt

OverOptOverSP

-0.25

-0.2

-0.15

-0.1

-0.05

0

0.05

0.1

1200 1300 1400 1500 1600 1700 1800

Rel

. Feh

ler

Que

rkra

ft

Zeitschritt

OverOptOverSP

-0.03

-0.02

-0.01

0

0.01

0.02

0.03

1200 1300 1400 1500 1600 1700 1800

Rel

. Feh

ler

Wirk

ungs

grad

Zeitschritt

OverOptOverSP

Anhang C

Graphen der Krafte fur VSP mit

Anbauteilen

137

ANHANG C. GRAPHEN DER KRAFTE FUR VSP MIT ANBAUTEILEN 138

Abbildung C.1: VSP Einflugler mit Schutzplatte

-60000

-50000

-40000

-30000

-20000

-10000

0

10000

800 1000 1200 1400 1600 1800

Sch

ub

Zeitschritt

SlidSPHoleSPOver

-40000

-30000

-20000

-10000

0

10000

20000

30000

800 1000 1200 1400 1600 1800

Que

rkra

ft

Zeitschritt

SlidSPHoleSPOver

-1

-0.8

-0.6

-0.4

-0.2

0

0.2

0.4

0.6

0.8

1

800 1000 1200 1400 1600 1800

Wirk

ungs

grad

Zeitschritt

SPHoleSPOver


Abbildung C.2: Krafte auf vereinfachte Schutzplatte, VSP Einflugler

1000

2000

3000

4000

5000

6000

7000

8000

9000

800 1000 1200 1400 1600 1800

Fx

Zeitschritt

SlidingSPOverSPHole

-3000

-2500

-2000

-1500

-1000

-500

0

500

1000

1500

2000

800 1000 1200 1400 1600 1800

Fy

Zeitschritt

SlidingSPOverSPHole


Abbildung C.3: VSP Ein- und Funfflugler mit SP und 3 Streben, Krafte 1. Flugel

-60000

-50000

-40000

-30000

-20000

-10000

0

10000

20000

30000

800 1000 1200 1400 1600 1800

Sch

ub

Zeitschritt

FullOneFullOneSlid

Slid

-40000

-30000

-20000

-10000

0

10000

20000

30000

40000

800 1000 1200 1400 1600 1800

Que

rkra

ft

Zeitschritt

FullOneFullOneSlid

Slid

-1

-0.5

0

0.5

1

800 1000 1200 1400 1600 1800

Wirk

ungs

grad

Zeitschritt

FullOneFull


Abbildung C.4: VSP Funfflugler mit SP und 3 Streben, summierte Krafte

-54000

-52000

-50000

-48000

-46000

-44000

-42000

-40000

1200 1300 1400 1500 1600 1700

Sch

ub

Zeitschritt

Full5xOne

Slid

-18000

-16000

-14000

-12000

-10000

-8000

-6000

-4000

-2000

1200 1300 1400 1500 1600 1700

Que

rkra

ft

Zeitschritt

Full5xOne

Slid

0.54

0.56

0.58

0.6

0.62

0.64

0.66

0.68

0.7

1200 1300 1400 1500 1600 1700

Wirk

ungs

grad

Zeitschritt

FullSlidWingonly


Abbildung C.5: Krafte auf Streben bei VSP Einflugler

-500

0

500

1000

1500

2000

2500

800 1000 1200 1400 1600 1800

Fx

Zeitschritt

OneS1OneS2OneS3

-4000

-3000

-2000

-1000

0

1000

2000

3000

4000

5000

800 1000 1200 1400 1600 1800

Fy

Zeitschritt

OneS1OneS2OneS3

-4000

-3000

-2000

-1000

0

1000

2000

3000

4000

5000

800 1000 1200 1400 1600 1800

Mz

Zeitschritt

OneS1OneS2OneS3


Abbildung C.6: Krafte auf Schutzplatte bei VSP Einflugler

0

1000

2000

3000

4000

5000

6000

7000

8000

800 1000 1200 1400 1600 1800

Fx

Zeitschritt

SlidingOneFullSPHole

-1500

-1000

-500

0

500

1000

800 1000 1200 1400 1600 1800

Fy

Zeitschritt

SlidingOneFullSPHole


Abbildung C.7: Krafte auf Streben bei VSP Funfflugler

-1500

-1000

-500

0

500

1000

1500

2000

2500

3000

3500

4000

1200 1300 1400 1500 1600 1700 1800

Fx

Zeitschritt

OneS1OneS2OneS3

-6000

-4000

-2000

0

2000

4000

6000

8000

1200 1300 1400 1500 1600 1700 1800

Fy

Zeitschritt

OneS1OneS2OneS3

-1500

-1000

-500

0

500

1000

1500

2000

1200 1300 1400 1500 1600 1700 1800

Mz

Zeitschritt

OneS1OneS2OneS3


Abbildung C.8: Krafte auf Schutzplatte bei VSP Funfflugler (eingeschrankter Wer-

tebereich)

200

400

600

800

1000

1200

1400

1600

1800

2000

800 1000 1200 1400 1600 1800

Fx

Zeitschritt

SlidingFull

-1000

-800

-600

-400

-200

0

200

800 1000 1200 1400 1600 1800

Fy

Zeitschritt

SlidingFull

Literaturverzeichnis

[Atta 1981] Atta, E. H.: Component adaptive grid interfacing. 1981. – AIAA

Paper Nr. 81-0382

[Bartels und Jurgens 2006] Bartels, J.-E. ; Jurgens, D.: The Voith Schneider

Propeller : current applications and new developments. Heidenheim : Voith Turbo

Marine GmbH & Co. KG, 2006

[Benek u. a. 1986] Benek, J. A. ; Steger, J. L. ; Dougherty, F. C. ; Buning,

P. G.: Chimera: A grid-embedding technique. 1986. – AD-A167466; AEDC-TR-

85-64; NAS 1.15:89246; NASA-TM-89246

[Bollhofer und Mehrmann 2004] Bollhofer, M. ; Mehrmann, V.: Numerische

Mathematik : eine projektorientierte Einfuhrung fur Ingenieure, Mathematiker und

Naturwissenschaftler. Wiesbaden : Vieweg, 2004. – ISBN 9783528032203

[Braess 2003] Braess, D.: Finite Elemente, Theorie, schnelle Loser und Anwen-

dungen in der Elastizitatstheorie. Berlin Heidelberg [u.a.] : Springer, 2003. – ISBN

3-540-00122-0

[Brown u. a. 1999] Brown, D. L. ; Henshaw, W. D. ; Quinlan, D. J.: OVER-

TURE: Object-oriented tools for overset grid applications. 1999. – AIAA Paper Nr.

99-3130, Centre for Applied Scientific Computing, Lawrence Livermore National

Laboratory, Livermore, California

[Chesshire und Henshaw 1990] Chesshire, G. ; Henshaw, W. D.: Composite

overlapping meshes for the solution of partial differential equations. In: Journal

Computational Physics 90 (1990), S. 1–64

[Chung u. a. 2007] Chung, B. ; Johnson, P.C. ; Popel, A.S.: Application of

Chimera grid to modelling cell motion and aggregation in a narrow tube. In:

International journal for numerical methods in fluids ISSN 0271-2091 53 (2007),

Nr. 1, S. 105–128

146

LITERATURVERZEICHNIS 147

[Comet User Manual 2001] Comet User Manual: Version 2.0. 2001. – ICCM

Institute of Computational Continuum Mechanics GmbH, Hamburg

[Djomehri und Biswas 2003] Djomehri, M. J. ; Biswas, R.: Performance en-

hancement strategies for multi-block overset grid cfd applications. 2003. – NAS

Technical Report Nr. 03-011

[Dougherty 1985] Dougherty, F. C.: Development of a Chimera Grid Scheme

with Applications to Unsteady Problems. 1985. – Dissertation, Stanford University

[Ferziger und Peric 2008] Ferziger, J. H. ; Peric, M.: Numerische Stromungs-

mechanik. Berlin Heidelberg [u.a.] : Springer, 2008. – ISBN 9783540675860

[Glowinski 1988] Glowinski, R.: First International Symposium on Domain De-

composition Methods for Partial Differential Equations, Ecole Nationale des Ponts

et Chaussees, Paris, France, January 7 - 9, 1987. Philadelphia, Pa. : Soc. for In-

dustrial and Applied Mathematics, 1988. – ed. by Roland Glowinski. – ISBN

0898712203

[Griebel u. a. 1995] Griebel, M. ; Dornseifer, T. ; Neunhoeffer, T.: Nume-

rische Simulation in der Stromungsmechanik : eine praxisorientierte Einfuhrung.

Braunschweig u.a. : Vieweg, 1995. – ISBN 3528067616

[Hadzic 2005] Hadzic, H.: Development and application of a finite volume method

for the computation of flows around moving bodies on unstructured, overlapping

grids. 2005. – Dissertation, Techn. Univ., Hamburg-Harburg, Institut fur Fluid-

dynamik und Schiffstheorie

[Henshaw 1985] Henshaw, W. D.: Pt.1: The numerical solution of hyperbolic

systems of conservation laws Pt.2: Composite overlapping grid techniques. 1985. –

Dissertation, California Inst.of Technology

[Henshaw 2008] Henshaw, W. D.: Parallel computation of three-dimensional flows

using overlapping grids with adaptive mesh refinement. In: Journal of Computa-

tional Physics 227 (2008), S. 7469–7502

[Jurgens und Fork 2002] Jurgens, B. ; Fork, W.: Faszination Voith-Schneider-

Propeller : Geschichte und Technik. Hamburg : Koehler, 2002. – ISBN 3782208544

[Jurgens u. a. 2007] Jurgens, D. ; Palm, M. ; Singer, S. ; Urban, K.: Numerical

Optimization of the Voith-Schneider-Propeller. In: ZAMM 87 (2007), Nr. 10,

S. 698–710


[Lohner 1995] Lohner, R.: Robust, vectorized search algorithms for interpolation

on unstructured grids. In: Journal of Computational Physics 118 (1995), S. 380–

387

[Lohner u. a. 2001] Lohner, R. ; Sharov, D. ; Lou, H. ; Ramamurti, R.: Over-

lapping unstructured grids. 2001. – AIAA Paper Nr. 01-0439

[Lohner u. a. 2002] Lohner, R. ; Yang, C. ; Cebral, J. ; u.a.: Advances in

FEFLO. 2002. – AIAA Paper Nr. 02-1024

[Madrane u. a. 2002] Madrane, A. ; Heinrich, R. ; Gerhold, T.: Implemen-

tation of the chimera method in the unstructured hybrid DLR finite volume Tau-

Code. 2002. – 6th Overset Composite Grid and Solution Technology Symposium

Ft. Walton Beach, Florida, USA October 8-10, 2002

[Madrane u. a. 2004] Madrane, A. ; Raichle, A. ; Stuermer, A.: Parallel

implementation of a dynamic overset unstructured grid approach. 2004. – Euro-

pean Congress on Computational Methods in Applied Sciences and Engineering,

ECONMAS 2004

[Manteuffel 1980] Manteuffel, T.: An Incomplete Factorization Technique for

Positive Definite Linear Systems. In: Mathematics of Computation 150(34) (1980),

S. 473–497

[Mastin und McConnaughey 1984] Mastin, C. W. ; McConnaughey, H. V.:

Computational Problems on Composite Grids. 1984. – AIAA Paper Nr. 84-1611

[Mathew 2008] Mathew, T. P. A.: Domain Decomposition Methods for the Nu-

merical Solution of Partial Differential Equations. Berlin u.a. : Springer, 2008. –

ISBN 9783540772057

[Meakin 1994] Meakin, R. L.: On the spatial and temporal accuracy of overset

grid methods for moving body problems. 1994. – AIAA Paper Nr. 94-1925

[Meakin 1995] Meakin, R. L.: The Chimera Method of Simulation for Unsteady

Three-Dimensional Viscous Flow. In: Computational Fluid Dynamics Review 1995,

ed. M. Hafez, K. Oshima (1995), S. 70–86

[Meakin und Wissink 1999] Meakin, R. L. ; Wissink, A. M.: Unsteady aerody-

namic simulation of static and moving bodies using scalable computers. 1999. –

AIAA Paper Nr. 99-3302


[Morton und Mayers 2005] Morton, K. W. ; Mayers, D. F.: Numerical solution

of partial differential equations : an introduction. Cambridge u.a. : Cambridge

Univ. Press, 2005. – ISBN 0521607930

[Muzaferija 1994] Muzaferija, S.: Adaptive finite volume method for flow pre-

diction using unstructured meshes and multigrid approach. 1994. – Dissertation,

University of London

[Quarteroni u. a. 2002] Quarteroni, A. ; Sacco, R. ; Saleri, F.: Numerische

Mathematik. Bd. 1. Berlin u.a. : Springer, 2002. – ISBN 3540678786

[Quarteroni und Valli 1999] Quarteroni, A. ; Valli, A.: Domain decomposition

methods for partial differential equations. Oxford u.a. : Clarendon Press, 1999. –

ISBN 0198501781

[Schwarz 1870] Schwarz, H. A.: Uber einen Grenzubergang durch Alternirendes

Verfahren. In: Vierteljahresschrift der Naturforschenden Gesellschaft in Zurich 15

(1870), S. 272–286

[Starius 1977] Starius, G.: Composite mesh difference method for elliptic and

boundary value problems. In: Numerical Mathematics 28 (1977), S. 243–258

[Toselli und Widlund 2005] Toselli, A. ; Widlund, O. B.: Domain Decompo-

sition Methods - Algorithms and Theory. Bd. 34. Berlin u.a. : Springer, 2005. –

Springer Series in Computational Mathematics. – ISBN 3540206965

[van der Vorst 1992] van der Vorst, H.: Bi-CGSTAB: a fast and smoothly

converging variant of Bi-CG for the solution of non-symmetric linear systems. In:

SIAM Journal on Scientific and Statistical Computing 12 (1992), S. 631–644

Ehrenwortliche Erklarung

Ich erklare hiermit ehrenwortlich, dass ich die vorliegende Arbeit selbststandig ange-

fertigt habe; die aus fremden Quellen direkt oder indirekt ubernommenen Gedanken

sind als solche kenntlich gemacht. Die Arbeit wurde bisher keiner anderen Prufungs-

behorde vorgelegt und auch noch nicht veroffentlicht.

Ich bin mir bewusst, dass eine unwahre Erklarung rechtliche Folgen haben wird.

Ulm, den 12.November 2008

(Unterschrift)

CFD-Berechnungen mit uberlappenden Gittern¨ in Bezug auf ...

Documents

Transcript of CFD-Berechnungen mit uberlappenden Gittern¨ in Bezug auf ...