scieee Science in your language
[en] (orig)
Effiziente Simulation komplexer Str¨
omungen auf
semi-strukturierten Gittern
vorgelegt von
M. Sc.
Jianping Yan
aus Zhejiang, V. R. China
Von der Fakult¨
at V – Verkehrs- und Maschinensysteme
der Technischen Universit¨
at Berlin
zur Erlangung des akademischen Grades
Doktor der Ingenieurwissenschaften
– Dr.-Ing. –
genehmigte Dissertation
Berlin 2003
D83
Effiziente Simulation komplexer Str¨
omungen auf
semi-strukturierten Gittern
vorgelegt von
M. Sc.
Jianping Yan
aus Zhejiang, V. R. China
Von der Fakult¨
at V – Verkehrs- und Maschinensysteme
der Technischen Universit¨
at Berlin
zur Erlangung des akademischen Grades
Doktor der Ingenieurwissenschaften
– Dr.-Ing. –
genehmigte Dissertation
Promotionsausschuss:
Vorsitzender: Prof. Dr. rer.nat. Dr.-Ing. habil. A. Dillmann
Berichter: Prof. Dr.-Ing. H. Eickhoff
Berichter: Prof. Dr.-Ing. F. Thiele
Tag der wissenschaftlichen Aussprache: 22. Januar 2003
Berlin 2003
D83
i
Vorwort
Die vorliegende Arbeit entstand w¨
ahrend meiner T¨
atigkeit als studentischer und wis-
senschaftlicher Mitarbeiter am Hermann-F¨
ottinger-Institut f¨
ur Str¨
omungsmechanik im
Fachgebiet Numerische Methoden der Thermo- und Fluiddynamik (CFD) von Herrn Prof.
Dr.-Ing. Frank Thiele.
F¨
ur seine langj¨
ahrige Unterst¨
utzung und F¨
orderung, die diese Arbeit erm¨
oglichte, m¨
ochte
ich Herrn Prof. Dr.-Ing. F. Thiele herzlich danken. Besonderer Dank gilt auch Herrn Prof.
Dr.-Ing H. Eickhoff f¨
ur die freundliche ¨
Ubernahme der Begutachtung und Herrn Prof. Dr.
rer.nat. Dr.-Ing. A. Dillmann f¨
ur die ¨
Ubernahme des Vorsitzes im Promotionsausschuss.
Weiterhin danke ich allen, die zum Gelingen dieser Arbeit beigetragen haben, im beson-
deren meinem Kollegen Markus Schatz f¨
ur die Durchsicht der Arbeit und die wertvollen
Verbesserungsvorschl¨
age.
F¨
ur ihre Unterst¨
utzung m¨
ochte ich meiner Frau Xianghong Li herzlich danken. Meinen
T¨
ochtern Julia Ziyi und Eveline Ziling ist diese Arbeit gewidmet.
Berlin, im Januar 2003 Jianping Yan
ii
Kurzfassung
In dieser Arbeit werden Techniken zur biharmonischen Gittergenerierung und effizienten
Gl¨
attung entwickelt und auf komplexe, praktisch relevante zwei- und dreidimensionale
Konfigurationen angewendet. Validierungen belegen, dass die durch biharmonische Ver-
fahren generierten Netze von hoher Qualit¨
at sind und dass Mehrgitterverfahren nicht nur
eine deutliche Konvergenzbeschleunigung, sondern auch eine verbesserte Stabilit¨
at ge-
gen¨
uber dem Eingitterverfahren erzielen.
Zur Beschleunigung der Str¨
omungssimulation mit dem SIMPLE-Algorithmus wird in der
vorliegenden Arbeit ein modifiziertes Full Multigrid (FMG) Verfahren mit V-Zyklus ent-
wickelt, bei dem die Startl¨
osung inklusive der Turbulenzvariablen nicht vom feinen Gitter
¨
ubertragen, sondern vom vorherigen Zyklus ¨
ubernommen werden. Die Validierung erfolgt
an ausgew¨
ahlten zwei- und dreidimensionalen Testf¨
allen.
Als wichtige technische Anwendungsf¨
alle werden die Simulation einer Stufenk¨
orper-Dif-
fusionsflamme und die instation¨
are Umstr¨
omung eines generischen KFZ-Außenspiegels
vorgestellt. F¨
ur die Simulation der Diffusionsflamme wird das Flamelet-Modell implemen-
tiert und Ergebnisse f¨
ur verschiedene Turbulenzmodelle erzeugt. Beim KFZ-Außenspiegel
wird der DES Ansatz im Vergleich zu einer URANS Simulation eingef¨
uhrt. Bei beiden
Testf¨
allen wird deutlich, dass die Wiedergabe der physikalischen Vorg¨
ange stark vom
Turbulenzmodellierungsansatz abh¨
angig ist.
Abstract
In this work, techniques for the biharmonic generation and efficient smoothing of grids
are developed and applied to complex, practically relevant two- and three-dimensional
configurations. The validation proves the grid generated by biharmonic methods to be
good quality, and that the multigrid method results in not only a clear convergence acce-
leration, but also an improved stability compared to the single grid method.
For the acceleration of flow simulations with the SIMPLE-algorithm, a modified Full Mul-
tigrid (FMG)-method with V-cycle has been developed, in which the starting solution,
including the turbulence quantities, is not taken from the finer grid but from the previous
cycle. Some two- and three-dimensional cases are validated.
As important technical applications the simulations of a bluff-body flame and the unsteady
flow around a generic car wing mirror are presented. For the diffusion flame simulation
the Flamelet-model is implemented and the results for different turbulence models are
obtained, while in the case of the mirror as an alternative to RANS, the DES formulation
is applied. Both test cases indicate that the reproduction of the physical precesses is
strongly dependent on the turbulence model used.
Inhaltsverzeichnis
Seite
1 Einleitung 1
2 Grundlagen der Str¨
omungssimulation 9
2.1 Gleichungen zur Gittergenerierung . . . . . . . . . . . . . . . . . . . . . . 9
2.1.1 Biharmonische Gleichungen . . . . . . . . . . . . . . . . . . . . . . 9
2.1.2 Gleichungen der Gittergenerierung . . . . . . . . . . . . . . . . . . 10
2.2 Str¨
omungsmechanische Bilanzgleichungen . . . . . . . . . . . . . . . . . . . 12
2.3 Modellierung turbulenter Str¨
omungen..................... 13
2.3.1 Lineare Wirbelz¨
ahigkeitsmodelle . . . . . . . . . . . . . . . . . . . . 14
2.3.2 Explizite Algebraische Spannungsmodelle . . . . . . . . . . . . . . . 18
2.3.3 Detached Eddy Simulation . . . . . . . . . . . . . . . . . . . . . . . 20
2.4 Verbrennungsmodelle.............................. 20
2.4.1 Flamelet-Modell............................. 22
2.4.2 Flamelet-Bibliothek und Kopplungsmodell . . . . . . . . . . . . . . 24
3 Numerische Methoden 25
3.1 Gittergenerierung................................ 25
3.1.1 Berechnung der Kontrollfunktionen . . . . . . . . . . . . . . . . . . 25
3.1.2 Erzeugung einer vorgegebenen Gitterschicht . . . . . . . . . . . . . 27
3.1.3 Approximation der Differentialgleichungen . . . . . . . . . . . . . . 28
3.1.4 Effektive Behandlung der Randbedingungen . . . . . . . . . . . . . 29
3.1.5 Globalgl¨
attung ¨
uber Blockgrenzen . . . . . . . . . . . . . . . . . . . 30
3.1.6 Datenstruktur und Datenzugriff . . . . . . . . . . . . . . . . . . . . 31
3.2 Str¨
omungssimulation .............................. 35
3.2.1 Finite Approximation der Bilanzgleichungen . . . . . . . . . . . . . 35
3.2.2 Modifizierte Behandlung der Wandrandbedingungen . . . . . . . . . 38
4 Grundlage des Mehrgitterverfahrens 41
4.1 Grundprinzip .................................. 41
4.2 Gittervergr¨
oberung und Gittertransfer . . . . . . . . . . . . . . . . . . . . 43
4.3 Mehrgitterzyklen ................................ 44
4.4 FMG-Methode ................................. 45
iii
iv Inhaltsverzeichnis
5 Gittergenerierung und Gl¨
attung mit Mehrgitterverfahren 47
5.1 Komponenten von Mehrgitterverfahren . . . . . . . . . . . . . . . . . . . . 47
5.2 Rechenablauf des Mehrgitterverfahrens . . . . . . . . . . . . . . . . . . . . 48
5.3 Anwendung des Gittergl¨
attungsverfahrens . . . . . . . . . . . . . . . . . . 51
5.3.1 2D ω-TypGeometrie.......................... 51
5.3.2 2D Fl¨
ugelkonfiguration......................... 55
5.3.3 3D Fl¨
ugelkonfiguration......................... 57
5.3.4 Fl¨
ugel-Rumpf Konfiguration . . . . . . . . . . . . . . . . . . . . . . 60
5.3.5 Fl¨
ugel-Rumpf-Leitwerk Konfiguration . . . . . . . . . . . . . . . . . 62
6 Mehrgitterverfahren zur L¨
osung der Navier-Stokes Gleichungen 65
6.1 Komponenten von Mehrgitterfahren . . . . . . . . . . . . . . . . . . . . . . 65
6.2 Modifizierte Full-Multigrid Methode . . . . . . . . . . . . . . . . . . . . . . 66
6.3 Grobgittergleichungen ............................. 67
6.4 Maßnahmen zur Turbulenzsimulation . . . . . . . . . . . . . . . . . . . . . 70
7 Validierung des Simulationsverfahrens 71
7.1 2D laminare Str¨
omung in rechteckiger Nische . . . . . . . . . . . . . . . . . 72
7.2 2D laminare Str¨
omung in verzerrter Nische . . . . . . . . . . . . . . . . . . 73
7.3 2D laminare Str¨
omung ¨
uber eine zur¨
uckspringende Stufe . . . . . . . . . . 76
7.4 2D turbulente Str¨
omung ¨
uber ein H¨
ugelmodell................ 79
7.5 3D laminare Nischenstr¨
omung......................... 84
7.6 3D Str¨
omung durch einen Rohrkr¨
ummer ................... 86
8 Simulation komplexer Str¨
omungen 95
8.1 CH4/H2Stufenk¨
orper-Diffusionsflamme . . . . . . . . . . . . . . . . . . . . 95
8.2 Generischer KFZ-Außenspiegel . . . . . . . . . . . . . . . . . . . . . . . . . 103
9 Zusammenfassung 109
Literaturverzeichnis ................................111
Kapitel 1
Einleitung
Komplexe Str¨
omungsph¨
anomene treten in vielen ingenieurtechnischen Anwendungen auf.
Zu ihrer Beschreibung werden neben den traditionellen Werkzeugen wie theoretischen
Methoden und experimentellen Untersuchungen auch numerische Str¨
omungssimulationen
verwendet. Die numerische Str¨
omungssimulationen ist dabei mittlerweile als ein integra-
ler Bestandteil sowohl von str¨
omungsphysikalischen Grundlagenuntersuchungen, wie z.B.
turbulenter Str¨
omungen mit chemischen Reaktionen, als auch von Auslegungsprozessen
f¨
ur praxisrelevante Konfigurationen, wie z.B. beim Entwurf von Verkehrs- und Transport-
flugzeugen zu sehen. Die heutigen Berechnungsmethoden erm¨
oglichen in der Vorentwurfs-
phase eine Reduzierung umfangreicher, zeitaufwendiger und kostspieliger experimenteller
Untersuchungen zur gezielten ¨
Uberpr¨
ufung und Verfeinerung des Entwurfs.
Ingenieurtechnische Problemstellungen sind meistens sowohl mit komplexen dreidimen-
sionalen Geometrien als auch mit komplexen turbulenten Austauschmechanismen ver-
bunden. Eine genaue numerische Simulation verlangt in jedem Fall den Einsatz von sehr
feinen Rechengittern guter Qualit¨
at und h¨
oherwertiger Modelle zur Beschreibung der
str¨
omungsphysikalischen Transportph¨
anomene. Dabei entstehen zwangsl¨
aufig große Sy-
steme gekoppelter diskreter Gleichungen. Die effiziente L¨
osung derartiger Gleichungssy-
steme, die durch große Datenmengen und aufwendige Rechenoperationen gekennzeichnet
ist, l¨
asst sich nur auf Hochleistungsrechnern unter Verwendung von optimierten nume-
rischen Algorithmen und Konvergenzbeschleunigungstechniken realisieren, wobei die nu-
merischen L¨
osungsverfahren an die Rechnerarchitektur angepasst werden m¨
ussen.
Dementsprechend wird die Effizienz und G¨
ute von numerischen Verfahren zur Str¨
omungs-
simulation im wesentlichen von den Teilbereichen Gittergenerierung, numerisches L¨
osungs-
verfahren und Modellbildung beeinflusst. Diese Teilbereiche k¨
onnen nicht unabh¨
angig von-
einander betrachtet werden, sondern stehen in Wechselwirkung miteinander. So k¨
onnen
Defizite z.B. in der Gittergenerierung auch durch sehr gute numerische L¨
osungsverfah-
ren nicht grundlegend kompensiert werden. Genauso erfordern Turbulenzmodelle, die den
wandnahen Str¨
omungsbereich erfassen, ein entsprechend angepasstes Grenzschichtgitter.
Ausschlaggebend f¨
ur die Effizienz eines Verfahrens zur Str¨
omungssimulation sind dement-
sprechend nicht nur die verwendeten numerischen Einzeltechniken, sondern ebenso die
Ber¨
ucksichtigung der Anforderungen aus der Modellierung und der Str¨
omungskonfigura-
1
2Kapitel 1. Einleitung
tion durch Gittergenerierung und L¨
osungsalgorithmus.
Gittergenerierung
Die numerische L¨
osung der Navier-Stokes Gleichungen erfordert insbesondere in der Grenz-
schicht eine hohe Gitteraufl¨
osung. Automatisch generierbare unstrukturierte Gitter ge-
n¨
ugen in vielerlei Hinsicht nicht den Genauigkeitsanforderungen, die von numerischen
Str¨
omungssimulationsverfahren gestellt werden. Gegen einen Einsatz von unstrukturierten
Gittern spricht in erster Linie die mangelhafte Approximationsgenauigkeit der turbulenten
Austauschmechanismen im k¨
orpernahen Bereich (Thompson 1996), weshalb technische
Anwendungen in der Aerodynamik derzeit meist auf k¨
orperangepassten strukturierten
Mehrblockgittern basieren. Zur Generierung derartiger Gitter verwendet man im wesent-
lichen algebraische Methoden, die im Gebietsinneren auf einer Interpolation zwischen den
Randpunkten basieren und Methoden, die elliptische oder hyperbolische Differentialglei-
chungssysteme l¨
osen. Die Erzeugung geeigneter strukturierter Rechennetze f¨
ur komplexe
Konfiguration ist somit weiterhin eine zeitaufwendige und kostspielige Aufgabe.
Aufgrund des vergleichsweise geringen Rechenaufwands eignen sich die algebraischen Me-
thoden in besonderem Maße zur interaktiven Netzgenerierung, die bei der rechnerischen
Erfassung komplexer dreidimensionaler Geometrien eine wichtige Rolle spielt. Der gr¨
oßte
Nachteil der algebraischen Netzgenerierung liegt jedoch darin, dass ihr keine Gl¨
attungs-
eigenschaften anhaften, denn Diskontinuit¨
aten im Verlauf der Gitterlinien, die durch die
Berandung aufgepr¨
agt werden, pflanzen sich unged¨
ampft ins Gebietsinnere fort. Im Ver-
gleich zu elliptischen erlauben hyperbolische Methoden eine bessere Steuerung der Ortho-
gonalit¨
at und Maschenweite. F¨
ur Geometrien mit vorgegebener Form im Außenbereich
sind sie dagegen aus mathematischen Gr¨
unden weniger gut geeignet (Weatherill 1990).
Die Motivation zur Benutzung elliptischer Gleichungssysteme liegt in deren Gl¨
attungsei-
genschaft. Viele elliptische Verfahren bauen auf den Laplace-Gleichungen auf. Der Nachteil
der auf diese Weise erzeugten Gitter ist deren Tendenz, Gitterlinien ¨
uber konvexen Beran-
dungen zu konzentrieren und ¨
uber konkaven R¨
andern auseinanderzuziehen. Daher werden
Kontrollfunktionen in die Laplace-Gleichungen eingef¨
uhrt und somit Poisson-Gleichungen
gel¨
ost. Hier besteht das Hauptproblem darin, auf welche Weise die Kontrollfunktionen auf
dem Rand zu bestimmen sind und wie sie in das Gebietsinnere interpoliert werden k¨
onnen.
Dazu wurden in der Vergangenheit unterschiedliche Verfahren entwickelt. Sie k¨
onnen bei-
spielsweise durch Projektion aus einem algebraischen Netz bestimmt werden (Thompson
1987). Mit solchen Kontrollfunktionen ergeben sich glatte Netze, die im wesentlichen die
Eigenschaften des algebraischen Netzes aufweisen. Daneben besteht eine zweite Methode
darin, die Kontrollfunktionen aus der Punkteverteilung auf der Berandung zu berechnen
und anschließend ins Gebietsinnere zu interpolieren (Steger und Sorenson 1979; Thomas u.
a. 1980; Thompson 1987; White 1990). Die Verteilung der Kontrollfunktionen im gesamten
Rechengebiet ist dabei einzig und allein durch die Punkteverteilung auf den Randfl¨
achen
bestimmt. Zwar kann sie dadurch auf einfache Weise direkt ermittelt werden, jedoch hat
man wenig M¨
oglichkeiten der Einflussnahme auf das Netz bez¨
uglich Maschenweite und
Orthogonalit¨
at am Rand und im Gebietsinneren k¨
onnen keine glatten Kontrollfunktionen
3
garantiert werden, was die L¨
osungsstabilit¨
at stark negativ beeinflusst.
Im Vergleich zum elliptischen Verfahren scheint der biharmonische Operator ∆∆φ= 0
das pr¨
adestinierte Verfahren zu sein (Sparis 1985). Er kann auf einfache Weise in eine
Poisson-Gleichung und eine Laplace-Gleichung zerlegt werden, auf deren Grundlage schon
jeweils sehr gute Gitteroperatoren entwickelt wurden. Wegen der Differentialgleichung
vierter Ordnung stehen f¨
ur die Gittergenerierung nicht nur Dirichlet- sondern zus¨
atz-
lich auch Neumann-Randbedingungen zur Verf¨
ugung. Durch die zweite Randbedingung
gibt es mehrere M¨
oglichkeiten, die Eigenschaften des Gitters zu steuern. Es k¨
onnen so-
wohl die Schnittwinkel der Gitterebenen mit den R¨
andern als auch die Maschenweiten
am Rand vorgegeben werden. Besonders hervorzuheben ist, dass die Kontrollfunktionen
einer Differentialgleichung gen¨
ugen und sie deshalb automatisch glatt sind. Im allgemei-
nen ist dies eine entscheidende Voraussetzung daf¨
ur, dass numerische Verfahren auch f¨
ur
komplexe Konfigurationen zu robusten Algorithmen f¨
uhren. Andererseits ist die L¨
osung
empfindlich gegen¨
uber Randbedingungen. Dies f¨
uhrt zu schlechter Stabilit¨
at und niedri-
ger Konvergenzrate des Verfahrens, weshalb effektivere Verfahren ben¨
otigt werden.
Aufgrund des elliptischen Charakters liefert die L¨
osung der oben erw¨
ahnten Differenti-
algleichungssysteme zur Gittergl¨
attung mittels eines direkten Verfahrens, z.B. der Gauß-
Elimination, in einem Schritt die geschlossene L¨
osung f¨
ur jeden Punkt. Bei großen Punkt-
anzahlen ist diese Vorgehensweise jedoch sehr rechenzeitintensiv. Effizienter ist dagegen
ein iteratives Verfahren, denn die Konvergenzgeschwindigkeit kann dabei durch Mehrgit-
terverfahren deutlich beschleunigt werden. Das grundlegende Prinzip des Mehrgitterver-
fahrens der Gittergenerierung wurde schon von St¨
uben und Linden (1985) vorgestellt und
von Spitaleri (1990; 1991) und Yan u. a. (1998) verwendet.
Turbulenzmodellierung
Die in der industriellen Praxis vorkommenden Str¨
omungen sind oftmals turbulent und
durch instation¨
are, dreidimensionale und stochastische Schwankungsbewegungen gekenn-
zeichnet. Die Beschreibung der turbulenten Ph¨
anomene kann in verschiedener Weise er-
folgen. Die Direkte Numerische Simulation (DNS) l¨
ost die exakten Erhaltungsgleichungen
f¨
ur die Masse, den Impuls und die Energie ohne Turbulenzmodellierung. Hierbei m¨
ussen
auch die kleinsten Skalen der turbulenten Schwankungsbewegung r¨
aumlich und zeitlich
aufgel¨
ost werden, w¨
ahrend bei der Grobstruktur- oder Large Eddy Simulation (LES) in
einer dreidimensionalen, instation¨
aren Rechnung nur alle ¨
uber einer bestimmten Gr¨
oße
liegenden Wirbel direkt numerisch behandelt werden. Kleinere, dissipative Skalen, die in
ihrer Struktur weniger anisotrop sind als die großen, energietragenden Wirbel, werden
herausgefiltert und durch sogenannte Feinstrukturmodelle (subgrid-scale model) erfasst.
Einen neueren ¨
Uberblick ¨
uber die LES-Verfahren ist bei Piomelli und Chasnov (1999) zu
finden. Probleme bereitet die LES nicht nur wegen des großen Bedarfs an Rechenleistung,
sondern auch wegen der Formulierung der zeitabh¨
angigen Randbedingungen sowie der
Behandlung des wichtigen wandnahen Bereichs, wenn die Skalen kleiner als die Filterwei-
te werden. Da bei r¨
aumlich variabler Filterweite zus¨
atzliche Terme in den Gleichungen
entstehen, die diese verkomplizieren, wurde die LES bisher haupts¨
achlich in Verbindung
4Kapitel 1. Einleitung
mit Wandfunktionen angewendet. Generell l¨
asst sich sagen, dass DNS und LES f¨
ur die
meisten praktisch relevanten Problemf¨
alle noch einen zu hohen numerischen Aufwand be-
deuten.
In technischen Anwendungen sind zumeist nur die statistischen Eigenschaften einer turbu-
lenten Str¨
omung von Bedeutung. Numerische Verfahren zur Berechnung praxisrelevanter
turbulenter Str¨
omung basieren daher fast ausnahmslos auf einer statistischen Betrach-
tung der Str¨
omung. Ausgangspunkt der Str¨
omungssimulation mit Turbulenzmodellen bil-
den die Reynolds-gemittelten Navier-Stokes Gleichungen (RANS). Neben den gemittelten
Geschwindigkeiten treten in diesen Gleichungen Korrelationsterme der Geschwindigkeits-
schwankungen von zweiter Ordnung auf, die sogenannten Reynoldsspannungen, die mit
Hilfe eines Turbulenzmodells auf gemittelte Impulsgr¨
oßen und eventuell auftretende Mo-
dellparameter zur¨
uckgef¨
uhrt werden. Reynoldsspannungsmodelle (RSM) l¨
osen f¨
ur jede
Komponente des Reynoldsspannungstensors eine eigene Transportgleichung und sind da-
mit aufgrund der h¨
oheren Modellierungsebene und der besseren Beschreibung der physi-
kalischen Vorg¨
ange, wie beispielsweise der Turbulenzanisotropie, den nachstehend behan-
delten Wirbelviskosit¨
atsmodellen ¨
uberlegen. Diese RSM Vorgehensweise wird aufgrund
des damit verbundenen Rechnenzeit- und Speicherbedarfs sowie numerischer Probleme in
praxisnahen Anwendungen nur vereinzelt eingesetzt.
Lineare, isotrope Wirbelviskosit¨
atsmodelle (Eddy Viscosity Model, EVM) st¨
utzen sich
auf einen nach Boussinesq (1887) benannten isotropen Zusammenhang zwischen den
Reynoldsspannungen und den Scherraten. Die bekanntesten Varianten davon sind z.B.
die Modelle von Baldwin und Lomax (1978) und Spalart und Allmaras (1992), das k−²
Modell (Jones und Launder 1972) und das k−ωModell (Wilcox 1988). Sie bilden die Stan-
dardmodelle f¨
ur derzeitige ingenieurtechnische Anwendungen. Aufgrund der beschr¨
ankten
M¨
oglichkeit, Anisotropiezust¨
ande der Reynoldsspannungen oder Stromlinienkr¨
ummungs-
einfluss wiederzugeben, liefern lineare Wirbelviskosit¨
atsmodelle bei komplexer dreidimen-
sionaler Str¨
omung mit stark gekr¨
ummten Stromlinien jedoch keine befriedigenden Ergeb-
nisse.
Zur Verbesserung der str¨
omungsphysikalischen Modellbildung werden deshalb nichtlinea-
re Erweiterungen der Zwei-Parameter-Wirbelz¨
ahigkeitsmodelle verfolgt. Diese Ans¨
atze
st¨
utzen sich auf einen im Unterschied zu isotropen EVM nichtlinearen Zusammenhang
zwischen den Reynoldsspannungen und den Geschwindigkeitsgradienten. Die nichtlinea-
ren Wirbelz¨
ahigkeitsbeziehungen k¨
onnen auf unterschiedliche Weise konstruiert werden.
Ein solcher Ansatz wird bei der Herleitung von Expliziten Algebraischen Spannungsmo-
dellen (EASM) verfolgt. Die Entwicklung und Analyse anisotroper Wirbelz¨
ahigkeitsbe-
ziehungen auf der Basis expliziter algebraischer Spannungsmodelle k¨
onnen den Arbeiten
von Rung u. a. (1999) und Rung (2000) entnommen werden.
Bei instation¨
aren Str¨
omungen stellt sich das Problem der Trennung der zeitlich variie-
renden mittleren Str¨
omung und der turbulenten Schwankungen. Trotz dieser nicht genau
definierten Trennung zwischen Turbulenz und mittlerer Str¨
omung werden Simulation auf
Basis von RANS h¨
aufig mit Erfolg bei instation¨
aren Problemen in der Form von Unsteady
5
RANS (URANS) eingesetzt. Die LES verlangt dagegen keine Aufteilung in Mittelwert und
Turbulenzanteil, sondern unterscheidet zwischen simulierten großen Schwankungen und
nicht-simulierten kleinen Schwankungen. Die LES ist dadurch f¨
ur instation¨
are Probleme
besser geeignet. Da der Anteil der nicht-simulierten Schwankungsbewegungen m¨
oglichst
klein gehalten werden sollte, ist allerdings der Aufwand gegen¨
uber RANS wesentlich h¨
oher.
Dies f¨
uhrte zur Entwicklung von kombinierten Verfahren, die die Vorteile von RANS und
LES vereinen sollen. Die Detached Eddy Simulation (DES) nutzt in Wandn¨
ahe ein Turbu-
lenzmodell zur Wiedergabe der kleinen Schwankungsbewegungen, simuliert gr¨
oßere dreidi-
mensionale Schwankungen im Fernfeld aber direkt wie eine LES, so dass sich Großr¨
aumige
Wirbelabl¨
osungen mit der Qualit¨
at einer LES erfassen lassen, wobei Gitteranforderungen
und Rechenzeiten mit denen einer RANS vergleichbar sind. Die Theorie und Anwendun-
gen von DES findet man in den Arbeiten von Spalart (1999; 2001), Shur u. a. (1999) und
Strelets (2000).
Verbrennungsmodellierung
Der Verbrennung selbst von relativ kleinen Kohlenwasserstoffen liegen sehr umfangreiche
Reaktionsmechanismen zugrunde. Zum Teil sind mehrere tausend Elementarreaktionen
(z.B. bei der Selbstz¨
undung von Cetan) am Gesamtgeschehen beteiligt. Das Wechselspiel
dieser Elementarreaktionen beeinflusst den gesamten Verbrennungsvorgang. Unabh¨
angig
von den spezifischen Eigenschaften des Brennstoffes weisen alle Reaktionsmechanismen
Eigenschaften auf, die f¨
ur Verbrennungsprozesse charakteristisch sind. So ergibt sich z.B.,
dass selbst bei großen Reaktionsmechanismen nur wenige Elementarreaktionen die Ge-
samtgeschwindigkeit beeinflussen. Auf dieser Grundlage wird in der numerischen Ver-
brennungssimulation oftmals eine Vereinfachung der Reaktionsmechanismen verwendet.
Dies ist von besonderem Interesse, da die Simulation von realen Systemen (dreidimen-
sional, turbulent wie z.B. in Motoren oder Feuerungen) sonst nur mit einem nicht zu
bew¨
altigenden Rechenzeitaufwand durchgef¨
uhrt werden k¨
onnten. Trotz der Vereinfachung
enthalten die Reaktionsmechanismen noch zweistellige Anzahlen von Spezies, die sich oft
durch die differentiellen Gleichungen l¨
osen lassen. Außerdem besteht f¨
ur die turbulente
Verbrennung noch das Schließungsproblem der Reaktionsgeschwindigkeit. Die statistische
Behandlung mit Hilfe von Wahrscheinlichkeitsdichtefunktionen (PDF) bietet hier einen
m¨
oglichen Ansatz. Allerdings ist die numerische Simulation der turbulenten Verbrennung
mit PDF-Transportgleichungen (Pope 1981; Pope 1985) oder Conditional Moment Closure
(CMC) (Bilger 1993) zu zeitaufwendig und deshalb nicht f¨
ur ingenieurtechnische Probleme
geeignet. F¨
ur ingenieurtechnische Anwendungen wird nicht nur Effizienz sondern auch die
F¨
ahigkeit erfordert, die Verteilung einzelner Spezies zu bestimmen. Wegen der Effizienz
und Einfachheit kommt deswegen in dieser Arbeit das Flamelet-Modell zur Anwendung.
6Kapitel 1. Einleitung
Numerische Effizienz
Mehrgitterverfahren
Die numerische Simulation von laminaren Str¨
omungen auf Basis der Navier-Stokes Glei-
chungen erfordert bereits eine große Anzahl von Gitterpunkten, um komplexe Str¨
omungs-
ph¨
anomene aufzul¨
osen und eine vom Gitter unabh¨
angige L¨
osung zu erhalten. F¨
ur turbu-
lente Str¨
omungen werden vor allem bei hoher Reynoldszahl noch feinere Gitter ben¨
otigt,
um die turbulente Scherschicht aufzul¨
osen. Die finite Diskretisierung der differentiellen
Gleichungen (RANS) f¨
uhrt zu einem System von algebraischen Gleichungen, dessen Gr¨
oße
jedoch von der Anzahl der Gitterpunkte abh¨
angt. Solche Gleichungssysteme k¨
onnen nu-
merisch effizient nur mit iterativen Verfahren gel¨
ost werden. Die am weitesten verbrei-
tete Methode zur iterativen L¨
osung der Navier-Stokes Gleichungen f¨
ur inkompressible
Str¨
omungen ist der SIMPLE Algorithmus (Patankar 1980). Wie bei den meisten der
existierenden L¨
osungsalgorithmen verlangsamt sich die Konvergenzgeschwindigkeit des
SIMPLE-Algorithmus mit zunehmender Gitterdichte.
Viele fortgeschrittene L¨
oser wie z.B. die Strongly Implicit Procedure (Stone 1968) bie-
ten einige Vorteile, sind aber nicht immer gleichm¨
aßig leistungsf¨
ahig. Von den bekannten
Beschleunigungstechniken haben sich in den letzten Jahren Mehrgitterverfahren als sehr
effiziente und allgemein anwendbare Verfahren herauskristallisiert, die die Konvergenzge-
schwindigkeit unabh¨
angig von der Gr¨
oße des zu l¨
osenden algebraischen Systems verbessern
k¨
onnen. Mehrgitterverfahren eignen sich deshalb besonders f¨
ur komplizierte Str¨
omungen,
wo sehr feine Diskretisierungen ben¨
otigt werden.
Seit den achtziger Jahren werden Mehrgitterverfahren in vielen Simulationsprogrammen
f¨
ur laminare Str¨
omungen eingesetzt und in letzter Zeit gewinnen sie auch bei der Berech-
nung turbulenter Str¨
omungen immer mehr an Bedeutung. Untersuchungen haben gezeigt,
dass Mehrgitterverfahren nicht nur f¨
ur die Berechnung einfacher laminarer Str¨
omungen
sondern auch f¨
ur Str¨
omungen mit starker W¨
arme¨
ubertragung, z.B. Auftriebsstr¨
omun-
gen in Hohlr¨
aumen (Demirdˇzi´c u. a. 1992), Kristallwachstumsprozesse (Leister und Peri´c
1992) und Str¨
omungen mit Strahlungsw¨
arme¨
ubertragung (Demirdˇzi´c u. a. 1996) effizi-
ent sind. Bei laminaren Str¨
omungen werden abh¨
angig von der Problemstellung und den
benutzten Gitterebenen Beschleunigungen um bis zu zwei Gr¨
oßenordnungen erzielt. Die
Anwendung des Mehrgitterverfahrens zur Berechnung turbulenter Str¨
omungen (Bai und
Fucks 1992; Lien und Leschziner 1991; Lien u. a. 1994a; Mavripilis 1995; Zheng u. a. 1997;
Huurdeman 1999) ergibt hingegen viel niedrigere Beschleunigungsfaktoren.
Mehrgitterverfahren f¨
ur die Navier-Stokes Gleichungen basieren meist auf einem Full Mul-
tigrid (FMG)–Full Approximation Storage (FAS) Algorithmus, wobei die Variablen und
Residuen vom feinen Gitter auf das n¨
achste gr¨
obere restringiert werden. Die Initialmas-
senfl¨
usse ¨
uber die Kontrollvolumenfl¨
achen auf dem groben Gitter ergeben sich durch Sum-
mation der entsprechenden zwei oder vier Feingittermassenfl¨
usse. Bei dieser Technik er-
gibt sich das Problem, dass die Massenfl¨
usse, die restringierten Geschwindigkeiten und
die Turbulenzvariablen auf unterschiedlichen Gitterebenen nicht zusammenpassen. Um
7
dieses Problem zu ¨
uberwinden, wird in der vorliegenden Arbeit ein modifizierter Mehrgit-
teralgorithmus entwickelt (Yan and Thiele 1998; Yan 1999), wobei nur die Residuen von
einem feinen Gitter auf die n¨
achst gr¨
obere Gitterstufe transferiert werden. Auf diese Weise
l¨
asst sich das Problem vermeiden und die Mehrgitterstrategie sowie die Programmstruk-
tur vereinfachen sich. Das Verhalten des modifizierten Algorithmus wird in Kombination
mit Konvektionsschmeta h¨
oherer Ordnung und auf verschiedenen Gittern untersucht.
Parallele Effizienz
Die schnelle Entwicklung der Computertechnologie ¨
außert sich nicht nur in der Steige-
rung der Rechengeschwindigkeit und Speicherkapazit¨
at des einzelnen Rechners. Sie ¨
außert
sich auch in der erh¨
ohten Daten¨
ubertragungsgeschwindigkeit und in der Verbesserung
der Netzwerktopologie f¨
ur Mehrprezessormaschinen. Dementsprechend haben sich meh-
rere Parallelisierungstrategien in der numerischen Str¨
omungssimulation entwickelt. Eine
davon wird als raumbezogene Parallelisierung bezeichnet und stellt heutzutage die Stan-
dardtechnik der Parallelisierung dar. Hierbei wird das Rechengebiet in etwa gleichgroße
zusammenh¨
angende Gebiete unterteilt. Jeder Prozessor bearbeitet eine dieser Gruppen,
womit die l¨
osbare Problemgr¨
oße nur durch die gesamte Speicherkapazit¨
at aller beteilig-
ten Prozessoren begrenzt ist. Die Datenmenge, die bei dieser Strategie zur ¨
Ubertragung
anf¨
allt, ist proportional zur Anzahl der bei der Gebietszerlegung entstandenen Trenn-
fl¨
achen. Bei einem dreidimensionalen Problem ist das eine relative große Datenmenge, die
sich ohne deutliche Einbußen in der Gesamtleistung nur durch eine effiziente Datenkom-
munikationstechnik verarbeiten l¨
asst.
Zielsetzung
Entsprechend den dargestellten Problemen sowohl bei der Gittergenerierung als auch bei
der numerischen Simulation f¨
ur komplexe Str¨
omungskonfigurationen sollen das Mehrgit-
terverfahren und ein modifiziertes Mehrgitterverfahren zur effizienten L¨
osung der diffe-
rentiellen bzw. Bilanzgleichungen entwickelt werden. Durch das Mehrgitterverfahren soll
die biharmonische Gittergenerierung bzw. Gittergl¨
attung beschleunigt werden, so dass
Gitter f¨
ur komplexe Konfigurationen mit diesem entwickelten Verfahren gegl¨
attet werden
k¨
onnen. Zur Behandlung von Mehrblockstrukturen soll die Methode dar¨
uber hinaus auf
die Gl¨
attung der Gitterlinien ¨
uber Blockgrenzen erweitert werden.
Zur Simulationsbeschleunigung f¨
ur die Anwendung auf komplexe Str¨
omungskonfiguratio-
nen wird ein Standard Mehrgitterverfahren unter Verwendung eines FMG-Algorithmus
weiter entwickelt. Der modifizierte Algorithmus soll durch die Simulation von laminaren
und turbulenten Str¨
omungen zun¨
achst grundlegend untersucht und hinsichtlich der Mehr-
gitterbeschleunigung validiert werden.
Zum Abschluss sollen als wichtige technische Anwendungsf¨
alle die Simulation einer Stu-
fenk¨
orper-Diffusionsflamme und die instation¨
are Umstr¨
omung eines generischen KFZ-
Außenspiegels vorgestellt werden. Bei beiden Testf¨
allen soll verdeutlicht werden, dass
8Kapitel 1. Einleitung
die Wiedergabe der physikalischen Vorg¨
ange stark vom Turbulenzmodellierungsansatz
abh¨
angig ist.
Gliederung
Das zweite Kapitel ist den Grundlagen und somit den Gleichungen der Gittergenerierung
und der Str¨
omungsmechanik gewidmet, wonach im Kapitel drei die numerischen Me-
thoden zur Gittergenerierung und Str¨
omungssimulation inklusiver der Behandlung der
R¨
ander beschrieben werden. Dabei wird auf ein wichtiges Problem der Gittergenerierung
eingegangen, n¨
amlich die Bestimmung der Kontrollfunktionen am Rand. Hier werden die
Gleichungen der Kontrollfunktionen am Rand durch Einf¨
uhrung der gew¨
unschten und ei-
ner virtuellen Gitterschicht hergeleitet, mit denen die Kontrollfunktionen iterativ berech-
net werden. Es wird ebenfalls das Prinzip beschrieben, wie Gitterlinien ¨
uber Blockgrenzen
gegl¨
attet werden k¨
onnen.
Ein Schwerpunkt der Arbeit ist die Entwicklung und Anwendung von Mehrgitterverfahren.
Im vierten Kapitel werden deshalb die f¨
ur die Entwicklung und Anwendungen relevanten
Grundlagen des Mehrgitterverfahrens er¨
ortert. Die Validierung der Gittergenerierung und
Gl¨
attung bezieht sich dabei vorwiegend auf praxisbezogene Konfigurationen, welche einen
vertieften Einblick in das Potential der Gittergenerierung erlauben. Die diesbez¨
uglichen
Anwendungen sind begleitend in den laufenden Text von Kapitel f¨
unf integriert. Hierzu
z¨
ahlen unter anderem die Gittergenerierung und Gittergl¨
attung von Ω-Typ Geometrien,
zwei- und dreidimensionale Fl¨
ugelkonfigurationen, eine Fl¨
ugel-Rumpf Konfiguration so-
wie eine Fl¨
ugel-Rumpf-Leitwerk Konfiguration.
In Kapitel sechs wird zun¨
achst kurz das Standard FMG-FAS dargestellt und darauf auf-
bauend ein modifizierter FMG-Algorithmus zur Anwendung in der Str¨
omungssimulation
entwickelt. Kapitel sieben beinhaltet die Verifizierung der Leistungsf¨
ahigkeit dieses Mehr-
gitterverfahrens an ausgew¨
ahlten akademischen Beispielen, deren Ergebnisse bekannt sind
und welche einen Geschwindigkeitsvergleich erlauben. Im Mittelpunkt stehen hier die la-
minare Str¨
omung in zweidimensionalen und dreidimensionalen Nischen, eine laminare
Str¨
omung ¨
uber eine zur¨
uckspringende Stufe, eine turbulente Str¨
omung ¨
uber ein zweidi-
mensionales H¨
ugelmodell, sowie Rohrkr¨
ummerstr¨
omungen mit starkem Sekund¨
arwirbel.
In Kapitel acht wird die Simulation zweier technische Anwendungsf¨
alle dargestellt und im
Rahmen der abschließenden Bemerkungen, Kapitel neun, werden die Schwerpunkte der
Arbeit zusammengefasst.
Kapitel 2
Grundlagen der
Str¨
omungssimulation
Im diesem Kapitel werden die der Simulation zugrunde liegenden Gleichungen zur Gitter-
generierung und die str¨
omungsmechanischen Bilanzgleichungen eingef¨
uhrt. Anschließend
werden die verwendeten Turbulenzmodelle und das Verbrennungsmodell in seinem we-
sentlichen Teilen beschrieben.
2.1 Gleichungen zur Gittergenerierung
F¨
ur die numerischen Str¨
omungssimulationen werden die Differentialgleichungen, die so-
wohl der Gittererzeugung als auch der Beschreibung der letztendlich zu berechnenden
physikalischen Vorg¨
ange zugrunde liegen, h¨
aufig in allgemein krummlinigen Koordina-
ten dargestellt. Die Bestimmung eines k¨
orperorientierten krummlinigen Koordinatensy-
stems reduziert sich schließlich auf das Problem, Transformationsbeziehungen xi(ξj) zu
finden, die die kartesischen Koordinaten xider Netzpunkte im physikalischen Raum in
Abh¨
angigkeit von den krummlinigen Koordinaten ξjdarstellen. Dazu verwendet man im
wesentlichen algebraische Methoden und solche, die auf der L¨
osung elliptischer oder hyper-
bolischer Differentialgleichungen basieren. Ein ¨
Ubersicht zu den gebr¨
auchlichen Verfahren
findet sich bei Thompson u. a. (1985). Die Verwendung von biharmonischen Gleichungen
zur Gittergenerierung wurde von St¨
uben u. a. (1986) vorgeschlagen und z.B. von Sparis
(1985) sowie Sparis und Karkanis (1992) eingesetzt. In der vorliegenden Arbeit konzen-
triert sich die Gittergenerierung nur auf das biharmonische Verfahren, da es Netze mit
guter Qualit¨
at erzeugen kann.
2.1.1 Biharmonische Gleichungen
Die biharmonischen Gleichungen wurden zuerst von Bell u. a. (1982) und Shubin u. a.
(1982) zur Gittergenerierung vorgeschlagen, wobei sie im transformierten Bereich f¨
ur die
physikalischen Koordinaten gel¨
ost wurden:
44x= 0 .(2.1)
9
10 Kapitel 2. Grundlagen der Str¨
omungssimulation
An dieser Stelle sei jedoch bemerkt, dass es kein Extremwert-Prinzip f¨
ur das Gleichungs-
system gibt. Infolgedessen kann die L¨
osung m¨
oglicherweise nicht vom transformierten zum
physikalischen Bereich zur¨
uck abgebildet werden. Deshalb wird die L¨
osung der transfor-
mierten biharmonischen Gleichungen im physikalischen Raum durchgef¨
uhrt:
44ξ= 0 (2.2)
mit den zwei vorgegebenen Randbedingungen
ξ=f , (2.3)
∂ξ
∂n =g . (2.4)
Dabei sind ξdie krummlinigen Koordinaten. Die biharmonischen Gleichungen k¨
onnen
grunds¨
atzlich in zwei gekoppelte Gleichungen vom Poisson-Typ aufgespalten werden:
4ξ=P , (2.5)
4P= 0 ,(2.6)
wobei Pdie Kontrollfunktionen sind.
2.1.2 Gleichungen der Gittergenerierung
Die direkte Anwendung der Gleichungen (2.5) und (2.6) kann bei komplexen Konfigura-
tionen große Werte der Kontrollfunktionen am Rand produzieren, die die Stabilit¨
at der
numerischen L¨
osung ung¨
unstig beeinflussen. Die Gr¨
oße der Kontrollfunktionen l¨
asst sich
durch die folgenden Gleichungen reduzieren:
4ξl=gllPl,(2.7)
4Pl= 0 .(2.8)
Die modifizierten Gleichungen k¨
onnen im transformierten Bereich wie folgt geschrieben
werden:
n
X
i,j=1
gijxl,ij +
n
X
i=1
giiPixl,i = 0 ,(2.9)
n
X
i,j=1
gijPl,ij +
n
X
i=1
giiPiPl,i = 0 .(2.10)
Hierbei sind gij die kontravarianten Metrikkoeffizienten. Die Gleichungen wurden auch
von Findling und Hermann (1991) benutzt, wobei zus¨
atzliche Parameter - die sogenann-
ten D¨
ampfungsfunktionen hij - eingef¨
uhrt werden, so dass in den Gleichungen (2.9 - 2.10)
2.1. Gleichungen zur Gittergenerierung 11
gij durch ˜gij =hijgij ersetzt wurden.
Durch die Multiplikation mit der Jacobi-Determinaten lassen sich die obigen Gleichungen
(2.9 - 2.10) wie folgt umformen:
F¨
ur den zweidimensionalen Fall:
α11x,ξξ + 2α12x,ξη +α22x,ηη +α11P1x,ξ +α22P2x,η = 0 ,(2.11)
α11P,ξξ + 2α12P,ξη +α22P,ηη +α11P1P,ξ +α22P2P,η = 0 ,(2.12)
mit
α11 =x,η ·x,η , α12 =−x,ξ ·x,η , α22 =x,ξ ·x,ξ .(2.13)
F¨
ur den dreidimensionalen Fall:
α11x,ξξ +α22x,ηη +α33x,ζζ + 2α12x,ξη + 2α13x,ξζ + 2α23x,ηζ (2.14)
+α11P1x,ξ +α22P2x,η +α33P3x,ζ = 0 ,
α11P,ξξ +α22P,ηη +α33P,ζζ + 2α12P,ξη + 2α13P,ξζ + 2α23P,ηζ (2.15)
+α11P1P,ξ +α22P2P,η +α33P3P,ζ = 0 ,
mit
α11 =g22g33 −g23g23 , α22 =g11g33 −g13g13 ,
α12 =g13g23 −g12g33 , α23 =g13g12 −g11g23 ,(2.16)
α13 =g12g23 −g13g22 , α33 =g11g22 −g12g12 .
Die kovarianten Metrikkoeffizienten gij sind von xabh¨
angig und werden bestimmt durch:
g11 =x,ξ ·x,ξ , g12 =x,ξ ·x,η , g13 =x,ξ ·x,ζ ,
g22 =x,η ·x,η , g23 =x,η ·x,ζ, g33 =x,ζ ·x,ζ .
(2.17)
Die Gleichungen (2.11) und (2.14) sind die als Inhomogeneous Thompson-Thames-Mastin
(ITTM)-Formeln bekannt und sind in der elliptischen Gittergenerierung weit verbreitet.
In der jetzigen Methode werden nicht nur die physikalischen Koordinaten xsondern auch
die Kontrollfunktionen Pdurch die Differentialgleichungen (2.12) und (2.15) bestimmt.
12 Kapitel 2. Grundlagen der Str¨
omungssimulation
2.2 Str¨
omungsmechanische Bilanzgleichungen
Die Str¨
omung eines Fluids wird mathematisch durch die differentiellen Bilanzgleichungen
in Zeit tund Raum x f¨
ur Masse, Impuls und Energie beschrieben. In ihrer konservativen
Form lauten sie:
Massenbilanz (Kontinuit¨
atsgleichung):
∂ρ
∂t +5·(ρu) = 0 .(2.18)
Impulsbilanz (Navier-Stokes Gleichungen f¨
ur den Geschwindigkeitsvektor u ):
∂(ρu)
∂t +5·(ρuu −π) = Su.(2.19)
Energiebilanz (Formulierung f¨
ur die Enthalpie h):
∂(ρh)
∂t +5·(ρuh +q) = Sh.(2.20)
F¨
ur Newtonsche Fluide lautet der Spannungstensor:
π=−pI + 2µS −2
3µ(5·u)I , (2.21)
mit dem Einheitstensor Iund dem Deformationsgeschwindigkeitstensor
S=1
2£5u+5Tu¤,(2.22)
wobei pden Druck und µdie dynamische Viskosit¨
at bezeichnen. Der Term Suist die
Volumenkraftdichte und beinhaltet z.B. die Wirkung der Schwerkraft.
F¨
ur die W¨
armestromdichte qverwendet man den Fourierschen Ansatz:
q=−λ5T. (2.23)
Die W¨
armeleitf¨
ahigkeit λist eine vom thermodynamischen Zustand abh¨
angige Stoffgr¨
oße.
F¨
ur thermisch und kalorisch ideale Stoffe mit der spezifischen Gaskonstante R, dem kon-
stanten Isentropenexponent γund der konstanten isobaren W¨
armekapazit¨
at cpgelten
dann die Zustandsgleichungen:
p=ρRT und h=cpTmit cp=γR
γ−1.
Mit der obigen Beziehung l¨
asst sich die W¨
armestromdichte durch
q=−µ
Pr 5h(2.24)
2.3. Modellierung turbulenter Str¨
omungen 13
darstellen, wobei Pr =µcp
λdie Prandtlzahl ist.
Die dynamische Viskosit¨
at µkann bei Gasen z.B. durch die halbempirische Formel von
Sutherland (1893)
µ=µref
Tref +Ts
T+TsµT
Tref 3
2
(2.25)
dargestellt werden. Darin ist µref die Viskosit¨
at f¨
ur die Referenztemperatur Tref und Ts
eine stoffabh¨
angige Konstante. F¨
ur Luft findet man (Jischa 1982):
Ts= 110 K, Tref = 280 K , µref = 1.75 ·10−5kg/ms . (2.26)
Der Quellterm Shf¨
ur die Enthalpie lautet:
Sh= (π+pI)··S+∂p
∂t +u·5p+wh.(2.27)
Hierin ist whdie volumenspezifische W¨
armequelle, z.B. durch W¨
armesstrahlung.
Bei gen¨
ugend kleiner Machzahl ( Ma ¿1 ), dem Verh¨
altnis der Str¨
omungs- zur Schall-
geschwindigkeit, k¨
onnen die Erhaltungsgleichungen f¨
ur inkompressible Str¨
omungen ver-
einfacht werden. Ist außerdem das Str¨
omungsproblem von der Temperatur unabh¨
angig,
dann entf¨
allt die Energiegleichung, und man erh¨
alt ein vereinfachtes System partieller Dif-
ferentialgleichungen zur Berechnung laminarer und turbulenter Str¨
omungen newtonscher
Fluide unter der Annahme konstanter Dichte ρund Viskosit¨
at µ.
2.3 Modellierung turbulenter Str¨
omungen
Die meisten technisch relevanten Str¨
omungen sind turbulent, wobei sehr hohe Reynolds-
zahlen auftreten k¨
onnen. Die Anwendung von DNS und LES, bei denen die physikalisch
relevanten Skalen ausreichend aufgel¨
ost werden m¨
ussen, beschr¨
ankt sich wegen des hohen
Rechenaufwands gegenw¨
artig noch auf kleinere bis mittlere Reynoldszahlen und eignet
sich damit nicht zur Berechnung praxisrelevanter Str¨
omungen. Deshalb werden bei tech-
nisch relevanten Anwendungen die zuf¨
alligen turbulenten Schwankungen der verschiede-
nen Str¨
omungsgr¨
oßen meist statistisch behandelt. Dazu bedient man sich der schon 1895
von Reynolds (Reynolds 1895) vorgeschlagenen Mittelungstechnik, bei der man zuerst die
Momentangr¨
oßen in Mittelwert φund Schwankungswert φ0teilt:
φ=φ+φ0.(2.28)
Bei der zeitlichen Mittelung der Bilanzgleichung verschwinden die in den Schwankungs-
gr¨
oßen linearen Terme. Die Anwendung dieses Prinzips auf die Impulsbilanzgleichungen
(Navier-Stokes Gleichungen) f¨
uhrt auf die sogenannten Reynoldsschen Gleichungen, hier
in ihrer inkompressiblen Form gegeben:
∂(ρ u)
∂t +5·³ρu u +ρu0u0−π´=Su.(2.29)
14 Kapitel 2. Grundlagen der Str¨
omungssimulation
Die daraus zu berechnenden Variablen der Reynoldsschen Gleichungen sind nicht mehr
die mit stochastischen Schwankungen behafteten physikalischen Geschwindigkeiten u, son-
dern ihre Mittelwerte u.
Die Gleichungen (2.29) gelten in der Form nur f¨
ur inkompressible Str¨
omungen, aber eine
formal identische Darstellung kann man mit der dichtegewichteten Mittelung (Favre 1965)
˜
φ=ρφ
ρ(2.30)
f¨
ur kompressible Str¨
omungen gewinnen. Dazu werden in die Reynoldsschen Gleichungen
anstelle der Ensemble-gemittelten Geschwindigkeiten udie dichtegewichteten Geschwin-
digkeiten ˜ueingesetzt, und die Dichte ρwird durch ihren Ensemble-Mittelwert ρersetzt.
Hierbei ¨
andert sich auch die Definition der Schwankungsgeschwindigkeiten
u0=u−˜u . (2.31)
Wegen der Nichtlinearit¨
at der Navier-Stokes Gleichungen treten bei der Mittelung zus¨
atz-
liche Unbekannte in Form des Reynoldsspannungstensors u0u0auf. Es handelt sich hier um
einen symmetrischen Tensor mit sechs Unbekannten, die durch Schließungsannahmen im
Rahmen der Turbulenzmodellierung auf die gemittelten Gr¨
oßen zur¨
uckgef¨
uhrt werden. Bis
heute ist eine große Anzahl von Turbulenzmodellen unterschiedlicher Komplexit¨
at und f¨
ur
verschiedene Anwendungsbereiche entwickelt worden. Abgesehen von den Reynoldsspan-
nungsmodellen, bei denen die Komponenten des Reynoldsspannungstensors mit Hilfe von
Transportgleichungen gel¨
ost werden, wird ¨
ublicherweise der Ansatz von Boussinesq (1877)
angewandt, wie im folgenden Kapitel dargestellt. Eine Zusammenfassung ¨
uber gebr¨
auchli-
che Turbulenzmodelle findet sich z.B. bei Wilcox (1993), im Rahmen dieser Arbeit erfolgt
lediglich eine kurze Zusammenstellung der verwendeten Turbulenzmodellierungsans¨
atze,
wobei zwischen linearen Wirbelz¨
ahigkeitsmodellen, Expliziten Algebraischen Spannungs-
modellen (EASM) und der Detached Eddy Simulation (DES) unterschieden wird.
2.3.1 Lineare Wirbelz¨
ahigkeitsmodelle
Analog zur mittleren freien Wegl¨
ange in einen Gas wird bei der Boussinesq Approximation
eine mittlere Mischungsl¨
ange lmund eine charakteristische Geschwindigkeit der turbulen-
ten Fluktuationsbewegung νmeingef¨
uhrt, um die Wirbelviskosit¨
at µtzu bestimmen:
µt∝ρνmlm.(2.32)
Der Reynoldsspannungstensor l¨
asst sich somit durch
τt=−ρu0u0=−2
3ρkI + 2µtS−2
3µt(5·u)I(2.33)
approximieren, wobei sich die turbulente kinetische Energie
k=1
2u0·u0(2.34)
2.3. Modellierung turbulenter Str¨
omungen 15
proportional zum Quadrat der charakteristischen Fluktuationsgeschwindigkeit k∝ν2
m
verh¨
alt. Die Modellierung reduziert sich demnach auf die Bestimmung der Gr¨
oßen lmund
νm, was z.B. in Form von ein oder zwei zus¨
atzlichen Transportgleichungen vorgenommen
werden kann.
SALSA–Modell
Bekannte Vertreter des Konzepts von Eingleichungs-Modellen sind die Ans¨
atze basierend
auf einer Transportgleichung f¨
ur die turbulente Viskosit¨
at oder eine vergleichbaren Gr¨
oße.
Hier m¨
ussen geeignete Annahmen f¨
ur die Bestimmung der Mischungsl¨
ange gefunden wer-
den. Zu erw¨
ahnen sind die Modelle von Baldwin und Lomax (1978) und von Spalart
und Allmaras (1992). Die hier durchgef¨
uhrten Simulationen verwenden eine Weiterent-
wicklung in Form des Strain Adaptive Linear Spalart-Allmaras (SALSA) Modells nach
Rung u. a. (2003). Ausgangspunkt ist das Modell von Spalart und Allmaras (1992) in der
von Edwards (1996) angegebenen Modifikation der Wandd¨
ampfungsterme. Die wesent-
lichen Bestandteile dieses Turbulenzmodells sind zum einen das Wirbelz¨
ahigkeitsprinzip
f¨
ur inkompressible Medien
ρuiuj=2
3ρkδij −2µtSij mit µt=fνρνt,(2.35)
sowie zum anderen eine semi-empirische Transportgleichung zur Bestimmung der Wir-
belz¨
ahigkeit νt:
∂ρνt
∂t +5·(ρuνt) = 5··µµ+µt
Prνt5νt¸+Pν+Dν−Yν.(2.36)
Hierin bezeichnen Pνden Produktionsterm, Dνden nicht-konservativen Diffusionsterm
und Yνden Vernichtungsterm. Dvund Yvwerden wie folgt definiert:
Dν=ρ5νt·5νt
Cb2
Prνt
, Yν=ρfw·Cb1
κ2+1 + Cb2
Prνt¸ν2
t
l2
n
.(2.37)
Die Darstellung des Produktionsterms lautet:
Pν=ρCb1νtS ,
Cb1= 0.1355√τ, τ = min(1.25,max(γ, 0.75)), γ = max(α1, α2),
α1= max h0,1−tanh ³νt
68´i0.65 , α2=µ1.01 νt
κ2l2
nS0.65
.(2.38)
Hierin kennzeichnet κdie von-K´arm´an-Konstante, lnden Wandnormalenabstand und S
die Invariante des Deformationsgeschwindigkeitstensors S:
S=q2S··S .
16 Kapitel 2. Grundlagen der Str¨
omungssimulation
Die angef¨
uhrten D¨
ampfungsfunktionen ergeben sich aus
fν=ν3
t
C3
ν1+ν3
t
, fw=g·1 + C6
w3
g6+C6
w3¸1/6
,
ν+
t=νt
ν, g =r£1 + Cw2(r5−1)¤, r = 1.313tanh µνt
Sκ2l2
n.(2.39)
Die verwendeten Konstanten lauten:
Cb2= 0.622, Cν1= 7.1, Cw2= 0.3, Cw3= 2.0, κ = 0.41, Prνt=2
3.(2.40)
k−²Modell
Den gr¨
oßten Anteil an Untersuchungen und Entwicklungen in der Turbulenzmodellierung
erfuhren in der Vergangenheit die Zweigleichungs-Turbulenzmodelle. Diese Modelle basie-
ren nicht nur auf der Berechnung der turbulenten kinetischen Energie k, sondern zus¨
atzlich
auf der Berechnung eines turbulenten L¨
angenmaßes oder einer davon abgeleiteten Gr¨
oße.
Da nach der Boussinesq-Approximation die turbulenten Austauschmechanismen von zwei
Gr¨
oßen (νmund lm) abh¨
angen, bildet ein Zweigleichungsmodell das einfachste Komplett-
system, um diesen Vorgang zu erfassen.
F¨
ur die Gleichung zur Bestimmung der turbulenten L¨
angen- Skala gibt es mehrere Ans¨
atze.
Abgesehen von Modellen f¨
ur die Dissipationsrate ²(Jones und Launder 1972) oder die
spezifische Dissipationsrate ω(Wilcox 1988) gibt es noch Modelle, die auf der turbulenten
L¨
angenskala lmoder der turbulenten Dissipationszeit τbasieren.
Die bekanntesten Vertreter der Zweigleichungsmodelle sind die sogenannten k−²und k−ω
Modelle, die jeweils aus zwei Transportgleichungen und einer Gleichung f¨
ur die turbulente
Viskosit¨
at µtbestehen. In der vorliegenden Arbeit wird das k−²Modell untersucht. Dies
kann in folgender Form geschrieben werden:
∂ρk
∂t +5·(ρuk) = 5··µµ+µt
σk5k¸+Pk−ρ² , (2.41)
∂ρ²
∂t +5·(ρu²) = 5··µµ+µt
σ²5²¸+C²1
²
kPk−C²2ρ²2
k,(2.42)
µt=Cµρk2
².(2.43)
Die verwendeten Modellkoeffizienten lassen sich aus theoretischen ¨
Uberlegungen und ex-
perimentellen Ergebnissen gewinnen. Launder und Spalding (1974) geben f¨
ur das k−²
Modell folgende Konstanten an:
σk= 1.0, σ²= 1.3, Cµ= 0.09, C²1= 1.44, C²2= 1.92 .(2.44)
2.3. Modellierung turbulenter Str¨
omungen 17
LL k−²Modell
Da das k−²Modell in seiner Originalform f¨
ur hohe Reynoldszahlen eine unrealistisch
hohe Turbulenzproduktion in der viskosen Unterschicht und den sich daran anschließend
¨
Ubergangsbereich bewirkt, wird das Modellverhalten mit Hilfe von geeigneten D¨
amp-
fungsfunktionen in diesem Bereich verbessert. Dazu wurden verschiedene Modelle ent-
wickelt wie z.B. Launder und Sharma (1974), Lam und Bremhost (1981), Chien (1982)
sowie Yang und Shih (1993). Als Variablen zur Steuerung der D¨
ampfungsfunktionen, die
hier mit fµ,f²1und f²2bezeichnet sind, werden die beiden turbulenten Reynoldszahlen
Rtund Rkoder der normierte Wandabstand Y+
Rt=k2
²ν , Rk=√kln
ν, Y +=lnpτw/ρ
ν(2.45)
verwendet. Der Einfluss der D¨
ampfungsfunktionen und ihre physikalische Bedeutung fin-
det sich in den Arbeiten von Patel u. a. (1985), Rodi und Mansour (1993) sowie Rousseau
u. a. (1997) .
Das von Lien und Leschziner (1993) entwickelte k−²Modell empfiehlt sich aufgrund seiner
hohen numerischen Stabilit¨
at und seiner g¨
unstigen wandnahen L¨
angenmaßeigenschaften.
Das Modell verfolgt eine ¨
ahnliche Philosophie wie die sogenannten Zweischichtenmodelle
(Rodi 1991), ohne jedoch explizit zwischen zwei verschiedenen Bereichen mit unterschied-
lichen Modellierungstechniken zu differenzieren. Die Transportgleichungen f¨
ur kund ²
lauten in allgemeiner Form:
∂ρk
∂t +5·(ρuk) = 5··µµ+µt
σk5k¸+Pk−ρ² , (2.46)
∂ρ²
∂t +5·(ρu²) = 5··µµ+µt
σ²5²¸+f²1C²1
²
kPk−f²2C²2ρ²2
k,(2.47)
µt=fµCµρk2
².(2.48)
Die Koeffizienten und D¨
ampfungsfunktionen sind wie folgt definiert:
σk= 1.0, σ²= 1.3, C²1= 1.44, C²2= 1.92, Cµ= 0.09,
f²1= 1 + P∗
Pk
, f²2= 1 −0.3e−R2
t, fµ=1−e−αµRk
1−e−α²Rk,(2.49)
P∗=ρf²2C²2k3/2
C²1l²
e−αdR2
k, l²=κC−3/4
µln(1 −e−α²Rk).
αd= 0.00222, α²= 0.263, αµ= 0.016 .
Im Vergleich zu den algebraischen Turbulenzmodellen und den Eingleichungs-Turbulenz-
modellen ist der Anwendungsbereich der Zweigleichungsmodelle als deutlich gr¨
oßer anzu-
sehen.
18 Kapitel 2. Grundlagen der Str¨
omungssimulation
LEA k−²Modell
Das zuvor beschriebene k−²Modell hat einen konstanten Anisotropieparameter Cµ=
0.09, der aus einer einschr¨
ankenden Annahme folgt. Im folgenden Kapitel 2.3.2 wird kurz
auf die Expliziten Algebraische Spannungsmodelle eingegangen, die auf Grund ihrer allge-
meing¨
ultigeren Herleitung einen variablen Anisotropieparameter anwenden. Beim EASM
wird der Reynoldsspannungstensor durch Gleichung (2.50) gegeben, und zur Berechnung
der Wirbelviskosit¨
at µtwird Cµdurch Gleichung (2.52) bestimmt. Verwendet man dieses
Berechnungsvorschrift ebenfalls f¨
ur das Cµdes k−²Modells, so beschr¨
ankt man sich nur
auf den linearen Term des Zusammenhanges zwischen Spannungen und Scherratentensor
des EASM. Es ergibt sich ein Modellierungsansatz, der als Lineares Explizites Algebrai-
sches (LEA) k−²Modell bezeichnet wird, der ¨
ahnliche Eigenschaften wie das EASM
besitzt, wie in Kapitel 8.1 gezeigt wird.
2.3.2 Explizite Algebraische Spannungsmodelle
Da alle bisher beschriebenen Modelle auf der Boussinesq Approximation basieren, unter-
liegt dabei die Reynoldsspannung einer linearen Proportionalit¨
at zum Deformationsge-
schwindigkeitstensor S:
−ρu0u0= 2µtS−2
3(ρk +µt5·u)I .
Diese lineare Beziehung zwischen den Reynoldsspannungen und den Scherraten ist aber
nicht allgemeing¨
ultig. Sie ist zum Beispiel f¨
ur Str¨
omungen, in denen extreme ¨
Ande-
rungen des Hauptspannungstensors Sauftreten, oder bei starker Stromlinienkr¨
ummung
nicht mehr erf¨
ullt und f¨
uhrt zu ungleichen Reynolds-Normalspannungen u0
iu0
i. Beispiele
hierzu sind drallbehaftete und dreidimensionale Str¨
omungen. Komplexe dreidimensionale
Str¨
omungen sind durch stark anisotrope Spannungszust¨
ande gekennzeichnet. Sowohl f¨
ur
k−²als auch f¨
ur k−ωModelle gibt es viele Korrekturterme, die f¨
ur einen oder mehrere
dieser Str¨
omungstypen verbesserte Simulationsergebnisse erm¨
oglichen.
Eine allgemeinere Bestimmung des Reynoldsspannungstensors besteht in Reynoldsspan-
nungsmodellen (Launder u. a. 1975), die die Anisotropie des Spannungstensors oder
Kr¨
ummungseffekte ohne Korrektur wiedergeben k¨
onnen, da f¨
ur jede Komponente eine
Transportgleichung gel¨
ost wird. Eine im Vergleich zu den Zweigleichungsmodellen noch
gr¨
oßere Zahl an Modellfunktionen und -konstanten erfordert jedoch noch intensive Unter-
suchungen, um in Richtung eines ”universellen” Reynoldsspannungsmodells zu gelangen.
Einen guten ¨
Uberblick ¨
uber die physikalischen Aspekte der einzelnen Terme findet man
in der Arbeit von Launder und Shima (1989) und Beispiele f¨
ur Anwendungen und Mo-
dellentwicklungen in den Arbeiten von Gibson und Younis (1986), Lai (1995), Launder
und Shima (1989) sowie Obi u. a. (1991).
Aus den impliziten Transportgleichungen der Reynoldsspannungen kann unter der Annah-
me homogener Turbulenzstruktur, also bei Vernachl¨
assigung der substantiellen Ableitung
2.3. Modellierung turbulenter Str¨
omungen 19
der Reynoldsspannungen durch darstellungstheoretische ¨
Uberlegungen ein explizites alge-
braisches Spannungsmodell in Gestalt eines nichtlinearen EVM formuliert werden. Dieses
kann zwar ebenfalls keine Kr¨
ummungseffekte wiedergeben, vermag aber Anisotropieefek-
te besser darzustellen. Der große Rechenaufwand impliziter RSM l¨
asst sich somit durch
die Verwendung vom EASM erheblich reduzieren, wobei die G¨
ute der Modellierung weit-
gehend erhalten bleibt. In dieser Arbeit wird exemplarisch das ”Realizable Quadratic
Eddy Viscosity Model” (Rung u. a. 1999) verwendet, das ein EASM darstellt, welches
zus¨
atzlich noch die Realizability-Bedingungen zu Gew¨
ahrleistung physikalischer Plausibi-
lit¨
at erf¨
ullt. Zur linearen Abh¨
angigkeit der Reynoldsspannungen von Skommen hier noch
Terme h¨
oherer Ordnung hinzu:
−ρu0u0= 2µtS−2
3(ρk +µt5·u)I+ 2µtB , (2.50)
mit
B=β2
k
²T2−β3
k
²T3,
wobei T2 und T3 quadratische Terme des Deformationsgeschwindigkeitstensors Sund des
Drehtensors Ω sind:
T2 = S·Ω−Ω·S, T3 = S·S−1
3(S··S)I . (2.51)
Der Anisotropieparameter Cµwird nicht konstant gehalten sondern in Abh¨
angigkeit von
Geschwindigkeitsgradienten bestimmt:
Cµ=3β1
3−2η2+ 6ξ2,(2.52)
mit
η2= (β3S)2/8, ξ2= (β2Ω)2/2.
Dabei werden die Modellkoeffizienten wie folgt angegeben:
β1= (2/3−C2/2)/g, β2= (1 −C4/2)/g, β3= (2 −C3)/g, g =C1+C5−1,
C1= 2.6, C2= 0.45, C3= 0.37, C4= 0.5, C5= 1.98 .
Die quadratische Abh¨
angigkeit der Reynoldsspannungen von den Scherraten ist in Ver-
bindung mit dem k−²Modell implementiert, wobei hierzu der Koeffizient C²2= 1.83
statt C²2= 1.92 eingesetzt wird.
Der bei der numerischen Umsetzung anders als die linearen Anteile explizit zu behandeln-
de, nichtlineare Anteil der Reynoldsspannungen wird im folgenden mit
τtn = 2µtB
bezeichnet.
20 Kapitel 2. Grundlagen der Str¨
omungssimulation
2.3.3 Detached Eddy Simulation
DES stellt in gewisser Weise einen Mittelweg zwischen LES und RANS dar. Die großska-
ligen Turbulenzstrukturen, deren L¨
angenmaß gr¨
oßer als das gew¨
ahlte numerische Gitter
ist, wird mit den Gitterweiten wie bei der LES berechnet, und die kleinskaligen Struktu-
ren, f¨
ur die das Gitter zu grob ist, werden durch RANS simuliert.
Die in der vorliegenden Arbeit verwendete DES basiert auf dem SALSA Modell, wobei das
L¨
angenmaß auf den Wandnormalenabstand lnzur¨
uckgef¨
uhrt wird. Die DES-Formulierung
erfordert im Wesentlichen das Ersetzten von lndurch einen neu definierten Abstand lDES,
der wie folgt berechnet wird (Spalart 1999):
lDES = min(ln, CDES4),4= max(4X, 4Y, 4Z),(2.53)
wobei 4X,4Yund 4Zdie Gitterweiten darstellen und die empirische Konstante CDES
zu 0.65 gesetzt ist. F¨
ur kleine Wandabst¨
ande, ln< CDES∆ (in der Grenzschicht), wird das
urspr¨
ungliche SALSA Modell verwendet. In der freien Str¨
omung ( ln> CDES∆)n¨
ahert
sich lDES der lokalen Gitterweite 4an. Damit funktioniert das SALSA Modell in dem
Bereich wie ein Feinstrukturmodell mit einer Mischungsl¨
ange, die gleich der Gitterweite
ist. Die Konstante ist dabei auf CDES = 0.65 gesetzt.
2.4 Verbrennungsmodelle
Turbulente Str¨
omungen in Kombination mit chemischen Reaktionen spielen eine wichti-
ge Rolle bei vielen technischen Verbrennungsprozessen. Bei der numerischen Simulation
reagierender Str¨
omungen werden zus¨
atzlich zu den Navier-Stokes Gleichungen die Glei-
chungen der Massenbr¨
uche Yider einzelnen Spezies iund oft der Temperatur Tgel¨
ost:
∂ρYi
∂t +5·(ρuYi) = 5·(ρDi5Yi) + ˙wi,(2.54)
∂ρT
∂t +5·(ρuT) = 5·(ρDT5T)
+1
cpÃN
X
i=1
cpiDi5Yi·5T−
N
X
i=1
hi˙wi+∂p
∂t + ˙qR!.(2.55)
Dabei ist N die Gesamtzahl aller auftretenden chemischen Komponenten, ˙wiist die Mas-
senproduktion jeder einzelnen Komponente infolge chemischer Reaktionen und hiist ihre
Enthalpie. cpbezeichnet die spezifische W¨
armekapazit¨
at der Gasmischung bei konstantem
Druck und ˙qRden Strahlungsverlust.
In der Theorie der Diffusionsflammen l¨
asst sich mit Hilfe der Elementmassenbr¨
uche Ziein
Mischungsbruch Z durch Z=(Zi−Zi2)/(Zi1−Zi2) definieren (die Indizes 1 und 2 bezeichnen
den Strom von Brennstoff bzw. Oxidationsmittel). In dieser Arbeit werden gleiche Diffu-
sionskoffizienten Df¨
ur alle chemischen Komponenten angenommen. Setzt man zus¨
atzlich
2.4. Verbrennungsmodelle 21
die Gleichheit von Stoff- und W¨
armefdiffusion voraus, wird das Diffusionsverhalten aller
chemischer Komponenten durch nur einen einzigen universalen Diffusionskoeffizienten D
beschrieben
D=Di=DT=λ
ρcp
=µ
SC
,1≤i≤N,(2.56)
und so sind die Lewiszahlen Leif¨
ur alle chemischen Komponenten identisch eins:
Lei=Le =λ
CpρD = 1 ,1≤i≤N.(2.57)
Mit diesem Ansatz gilt das folgende Erhaltungsgesetz:
∂ρZ
∂t +∇·(ρuZ) = ∇·(ρD∇Z),(2.58)
mit den Randbedingungen Z=0 im Oxidationsmittel und Z=1 im Brennstoffstrom.
Die Favre-gemittelten Gleichungen f¨
ur turbulente reagierende Str¨
omungen lauten:
∂ρe
Yi
∂t +5·(ρeue
Yi) = 5·(ρD 5e
Yi)−5·(ρg
u00Y00
i) + e˙wi,(2.59)
∂ρ e
T
∂t +5·(ρeue
T) = 5·(ρD 5e
T)−5·(ρg
u00T00)
+1
cpÃn
X
i=1
cpiD5e
Yi·5e
T−
N
X
i=1
hie˙w+∂p
∂t + ˙qR!.(2.60)
Hier beschreiben die Terme g
u00Y00
iund g
u00
jT00 den turbulenten Transport. Sie sind unbe-
kannt und m¨
ussen ¨
ahnlich wie die Geschwindigkeitskorrektur in den Reynoldsgleichungen
modelliert werden. Dazu wird oft der Gradientenansatz verwendet:
−ρg
u00
jY00
i=µt
Sct
∂˜
Yi
∂xj
,−ρg
u00
jT=µt
Sct
∂˜
T
∂xj
,(2.61)
wobei Sctdie turbulente Schmidtzahl ist.
Einer L¨
osung der gemittelten Erhaltungsgleichungen steht jetzt nur noch die Bestimmung
der mittleren Massenproduktion ˙wiim Wege. Diese ist eine hoch nichtlineare Funktion
der Massenbr¨
uche und der Temperatur. Durch die Nichtlinearit¨
at der zugrunde liegenden
physikalisch-chemischen Prozesse ist die Bestimmung von e˙wischwierig. Die statistische
Behandlung mit Hilfe von Wahrscheinlichkeitsdichtefunktionen (PDF) bietet hier einen
m¨
oglichen Ansatz.
Aus Sicherheitsgr¨
unden werden bei der industriellen Feuerung und in Brennern ¨
uberwie-
gend Diffusionsflammen eingesetzt. F¨
ur die numerische Simulation von Verbrennungspro-
zessen mit diesem Flammentyp sind zur Zeit das Verfahren der PDF-Transportgleichung
(Pope 1981; Pope 1985; Pope 1990) und das Flamelet-Modell (Peters 2000; Warnatz 1993)
weit verbreitet. Letzteres wird in dieser Arbeit verwendet.
22 Kapitel 2. Grundlagen der Str¨
omungssimulation
2.4.1 Flamelet-Modell
Das Flamelet-Konzept basiert auf der Annahme, dass eine turbulente Flamme aus vie-
len gewinkelten laminaren Fl¨
ammchen besteht. Deshalb m¨
ussen die chemischen Zeit- und
L¨
angenskalen nicht vom Str¨
omungsl¨
oser aufgel¨
ost werden (Peters 2000). Wegen seiner
Einfachheit ist dieses Konzept weit verbreitet und wird in Verbindung mit hochentwickel-
ten Turbulenzmodellen in den vorliegenden Untersuchungen verwendet.
Die Verwendung des Flamelet-Modells basiert auf der Annahme einer β-f¨
ormigen PDF,
deren Form durch den Favre-gemittelten Mischungsbruch e
Zund seine Varianz g
Z002wie
folgt bestimmt wird:
e
P(Z) = Zα−1(1 −Z)β−1Γ(α+β)
Γ(α)Γ(β).(2.62)
Dabei sind α,βund γ
α=e
Zγ , β = (1 −e
Z)γ , γ =e
Z(1 −e
Z)
g
Z002,(2.63)
und die Gamma-Funktion
Γ(x) = Z∞
0
tx−1e−tdt . (2.64)
F¨
ur diese zwei Variablen e
Zund g
Z002werden folgende Transportgleichungen gel¨
ost:
∂ρ e
Z
∂t +∇·(ρeue
Z) = ∇·(µ
Sc +µt
Sct
)∇e
Z , (2.65)
∂ρg
Z002
∂t +∇·(ρeug
Z002) = ∇·(µ
Sc +µt
Sct
)∇g
Z002+ 2 µt
Sct
(∇e
Z)2−ρeχ , (2.66)
mit
eχ=cχ
²
kg
Z002.(2.67)
Dabei werden die Schmidtzahl Sc und die turbulente Schmidtzahl Sctzu 0.7 bzw. 1.0
gesetzt. Der Modellparameter cχwird zu 2.0 gesetzt (Peters 2000).
Die Flamelet-Gleichungen k¨
onnen nach (Peters 2000) durch die Annahme der Existenz ei-
ner d¨
unnen, von der Turbulenz ungest¨
orten Reaktionszone mit einer Crocco-Transformation
in ein lokales, auf das Diffusionsflamelet bezogenes Koordinatensystem abgeleitet werden.
F¨
ur den Massenbruch Yiund die Temperatur Tlauten sie wie folgt:
ρ∂Yi
∂t −ρχ
2
∂2Yi
∂Z2−˙wi= 0 ,(2.68)
2.4. Verbrennungsmodelle 23
ρ∂T
∂t −ρχ
2
∂2T
∂Z2−ρχ
2cpÃN
X
i=1
cpi
∂Yi
∂Z +∂cp
∂Z !∂T
∂Z +1
cpÃN
X
i=1
hi˙wi−˙qR!= 0 .(2.69)
Weil die Lewiszahlen identisch Eins sind, treten sie hier in diesen Gleichungen nicht auf.
Diese Gleichungen sind eindimensional und scheinbar außer in der Zeit tnur noch in einer
einzigen weiteren unabh¨
angigen Variablen, dem Massenbruch Z, formuliert. Aber wegen
der Zustandgleichung p=ρRT ist der Druck pauch eine zus¨
atzliche Variable f¨
ur die
Flamelet-Gleichungen. Die verwendete Modellierung von ˙qRf¨
ur Wasserstoffflammen in
der d¨
unnen Reaktionszone wird von Pitsch u. a. (2000) gegeben:
˙qR= 2σ(T4−T4
u)(pCO2αp,CO2+pH2Oαp,H2O).(2.70)
Dabei bezeichnet σdie Stefan-Boltzmannkonstante, Tudie umgebende Temperatur, pi
den Partialdruck und αp,i den Plankschen Absorbtionskoeffizient der Spezies i.
In den Flamelet-Gleichungen wird die skalare Dissipationsrate χdurch folgenden Gradienten-
Diffusion-Ansatz definiert:
χ= 2D(∇Z)2.(2.71)
Sie stellt den Einfluss des Mischungsfeldes und somit der turbulenten Str¨
omung auf die
laminare Flammenstruktur dar. Sie hat die Einheit s−1und kann als Inverse einer cha-
rakteristischen Diffusionszeit interpretiert werden.
Die Abh¨
angigkeit der zur L¨
osung der Flamelet-Gleichungen ben¨
otigten skalaren Dissipa-
tionsrate vom Mischungsbruch wird modelliert als (Peters 2000)
χ=χst
f(Z)
f(Zst)(st : am st¨ochiometrischen Mischungsbruch Zst),(2.72)
mit
f(Z) = e−2[erfc−1(2Z)]2.(2.73)
Dabei ist erfc−1die Inverse der komplement¨
aren Fehlerfunktion. Gleichung (2.72) zeigt,
dass χnur vom Mischungsbruch Zsowie dem Wert am st¨
ochiometrischen Punkt χst
abh¨
angig ist. Da die beiden Gr¨
oßen χst und Z statistisch unabh¨
angig sind, kann die
turbulente Mittelung der skalaren Dissipationsrate wie folgt geschrieben werden:
eχ=Zχst
χst e
P(χst)dχst ZZ
f(Z)
f(Zst)e
P(Z)dZ . (2.74)
Dabei bezeichnet e
Pdie Favre-gemittelte PDF. In dieser Gleichung definiert die erste
Integration die gemittelte skalare Dissipationsrate am st¨
ochiometrischen Punkt als
< χst >=Zχst
χst e
P(χst)dχst .(2.75)
24 Kapitel 2. Grundlagen der Str¨
omungssimulation
Unter Verwendung von Gleichung (2.67) lautet die am st¨
ochiometrischen Mischungsbruch
gemittelte skalare Dissipationsrate:
< χst >=cχ²
kg
Z002f(Zst)
R1
0f(Z)e
P(Z)dZ .(2.76)
Mit Annahme von χst =< χst >kann das zur L¨
osung der Flamelet-Gleichungen ben¨
otigte
χwie folgt ermittelt werden:
eχZ=χ=cχ²
kg
Z002f(Z)
R1
0f(Z)e
P(Z)DZ =eχf(Z)
R1
0f(Z)e
P(Z)DZ .(2.77)
Dann erh¨
alt man f¨
ur die turbulente reagierende Str¨
omung:
ρ∂Yi
∂t −ρeχZ
2
∂2Yi
∂Z2−˙wi= 0 ,(2.78)
ρ∂T
∂t −ρeχZ
2
∂2T
∂Z2−ρeχZ
2cpÃN
X
i=1
cpi
∂Yi
∂Z +∂cp
∂Z !∂T
∂Z +1
cpÃN
X
i=1
hi˙wi−˙qR!= 0 ,(2.79)
mit der bedingt gemittelten skalaren Dissipationsrate eχZ.
Es sollte beachtet werden, dass zur Schließung der Flamelet-Gleichungen die Annahme
der bedingt gemittelten skalaren Dissipationsrate benutzt wird. Daraus folgt, dass es sich
auch bei der L¨
osung um bedingt gemittelte Gr¨
oßen handelt.
2.4.2 Flamelet-Bibliothek und Kopplungsmodell
Die Flamelet-Gleichungen besitzen die zwei Parameter eχZund p. Durch die L¨
osung der
Gleichungen mit unterschiedlichen Werten wird eine sogenannte Flamelet-Bibliothek er-
zeugt. Darin werden zu jedem Parameterpaar fχZBIB und pBIB die Massenbruchverteilung
Yi(Z) und die Temperaturverteilung T(Z) gespeichert. F¨
ur die numerische Simulation ei-
ner turbulenten reagierenden Str¨
omung gen¨
ugt es nun, aus den bekannten Feldern k,²,e
Z
und g
Z002die Parameter fχZund pzu bestimmen. Mit zwei ermittelten Variablen eχZund
pkann aus der Flamelet-Bibliothek die entsprechende Massenbruch- und Temperaturver-
teilung im Mischungsbruchraum Z(0 – 1) ausgelesen werden.
Liegt die Massenbruchverteilung Yi(Z) der chemischen Komponente iund die Tempe-
raturverteilung T(Z) im Mischungsbruchraum Zals Ergebnis der Flamelet-Gleichungen
vor, so erh¨
alt man die Werte ρ,e
Yiund e
Tim physikalischen Raum durch Integration mit
der PDF ¨
uber den Mischungsbruch:
1
ρ=Z1
0e
P(Z)
ρ(Z)dZ mit ρ(Z) = p
RT(Z)PN
i=1
Yi(Z)
Mi
,
e
Yi=Z1
0
Yi(Z)e
P(Z)dZ , e
T=Z1
0
T(Z)e
P(Z)dZ . (2.80)
Dabei bezeichnet Rdie allgemeine Gaskonstante (R=8,314 J/mol K) und Midie molare
Masse des Stoffes i.
Kapitel 3
Numerische Methoden
Dieses Kapitel befasst sich mit den numerischen Methoden zur Gittergenerierung und zur
Str¨
omungssimulation. Das Verfahren zur Gittergenerierung wird in seinen wesentlichen
Komponenten beschrieben, wobei speziell auf die Weiterentwicklung hinsichtlich der Be-
rechnung der Kontrollfunktion eingegangen wird. W¨
ahrend das numerische Verfahren zur
Str¨
omungssimulation auf der Grundlage der Bilanzgleichungen nur in einer ¨
Ubersicht dar-
gestellt ist, wird die Implementierung von Wandrandbedingungen f¨
ur stark schiefwinklige
Gitterlinien weiterentwickelt.
3.1 Gittergenerierung
In diesem Abschnitt wird zuerst die Berechnung der Kontrollfunktionen am Rand be-
schrieben und danach die Diskretisierung der Differentialgleichungen vorgestellt. Das hier
verwendete Verfahren von Linden (1984) zeichnet sich durch eine erh¨
ohte L¨
osungsstabi-
lit¨
at hinsichtlich der Randbedingungen sowie durch eine erh¨
ohte Konvergenzgeschwindig-
keit aus. Anschließend wird die Gl¨
attung ¨
uber Blockgrenzen beschrieben.
3.1.1 Berechnung der Kontrollfunktionen
Die biharmonischen Gleichungen zeichnen sich durch zwei Randbedingungen aus. Die erste
Randbedingung bezieht sich auf die Gitterpunktverteilung am Rand. Als zweite Randbe-
dingungen k¨
onnen die ersten Ableitungen x,ξ,x,η und x,ζ am Rand vorgegeben werden,
d.h. die Schnittwinkel und Maschenweiten. W¨
ahrend sich die Randbedingung f¨
ur die Glei-
chungen (2.11) und (2.15) direkt aus den Koordinaten der Randgitterschicht ergibt, kann
die Randbedingung f¨
ur die Gleichungen (2.12) und (2.16) mit der oben genannten zweiten
Randbedingung gewonnen werden. Um die Kontrollfunktionen Pam Rand herzuleiten,
werden hier, wie in Abbildung 3.1 gezeigt, drei Gitterschichten XB+1,XBund XB−1ein-
gef¨
uhrt. Davon ist XBdie Randgitterschicht, XB−1ist eine virtuelle Gitterschicht außer-
halb des Rechengebiets und XB+1 ist die gew¨
unschte erste Gitterschicht. Bei vorgegebener
erster Ableitung am Rand hat die virtuelle Gitterschicht immer den gleichen Abstand δx
zwischen den realen Punkten xB+1 und den gew¨
unschten Punkten XB+1. Die Schnitt-
winkel zum Rand werden bei der elliptischen Gl¨
attung deshalb weitgehend festgehalten.
Diese Vorgehensweise l¨
asst sich am Beispiel der Berechnung der Kontrollfunktionen f¨
ur
25
26 Kapitel 3. Numerische Methoden
den zweidimensionalen Fall demonstrieren.
P
N
S
WE





























Β+1
Β
Β−1
i
i−1
η
ξ
i+1
δx
δx
Gitterpunkt am Rand
gew¨
unschter Gitterpunkt
realer Gitterpunkt
virtueller Gitterpunkt
Abbildung 3.1: Schematische Darstellung der Behandlung der Randbedingung im zweidi-
mensionalen Fall
Mit den Beziehungen
δx =xB+1 −XB+1 , xB−1=XB−1+δx (3.1)
werden die Ableitungen der Koordinaten am Rand mit zentralen Differenzenschemata
berechnet:
xP
,ξ =(xE−xW)
2h=XP
,ξ , xP
,η =(xN−xS)
2h=XP
,η .(3.2)
Die Maschenweite h der kontravarianten Koordinaten wird ¨
ublicherweise zu eins gesetzt.
Zur vereinfachten Darstellung wird im folgenden der Hochindex Pweggelassen. Nach
Gleichung (3.2) ¨
andert sich die Metrikkoeffizienten
α11 =X,η ·X,η , α12 =−X,ξ ·X,η , α22 =X,ξ ·X,ξ (3.3)
nicht. Die zweiten Ableitungen x,ξξ und x,ξη bleiben ebenfalls konstant:
x,ξξ =X,ξξ , x,ξη =X,ξη .(3.4)
Nur x,ηη erf¨
ahrt eine ¨
Anderung infolge der Abweichung δx der randn¨
achsten Gitterschicht:
x,ηη =X,ηη + 2δx
h2.(3.5)
Damit ergibt sich aus Gleichung (2.11) die folgende Beziehung:
α11X,ξξ + 2α12X,ξη +α22X,ηη + 2α22 δx
h2+α11PB
1X,ξ +α22PB
2X,η = 0 .(3.6)
Weil die Metrikkoeffizienten jedoch bekannt sind, h¨
angen die Kontrollfunktionen dann
nur von δx ab und k¨
onnen am Rand iterativ bestimmt werden.
F¨
ur den dreidimensionalen Fall kann die Randbehandlung in analoger Weise f¨
ur jede
Fl¨
ache fortgesetzt werden. Es ergibt sich eine ¨
ahnliche Beziehung f¨
ur die Kontrollfunktio-
nen am Rand.
3.1. Gittergenerierung 27
3.1.2 Erzeugung einer vorgegebenen Gitterschicht
Eine vorgegebene Gitterschicht dient zur Berechnung der Kontrollfunktion am Rand,
wobei in der Praxis eine sinnvolle Wandabstandsverteilung der randn¨
achsten Gitterlinie
¨
uber alle W¨
ande vorzugeben ist. In vielen F¨
allen erf¨
ullt die einfache Vorgabe eines Ab-
stands nicht die Anforderungen, da die Simulation turbulenter Str¨
omungen angepasste
Wandabst¨
ande ben¨
otigt. Unter der Voraussetzung, dass die Abstandsverteilung eines al-
gebraischen Netzes bereits gut ist, kann diese aus dem Startgitter berechnet werden. Das
Verfahren ist f¨
ur den dreidimensionalen Fall in Abbildung 3.2 skizziert, wobei B, W, E, S,
N die Gitterpunkte auf der Wand sind und P der wandn¨
achste Gitterpunkt ist, w¨
ahrend
XB+1 den gew¨
unschten Gitterpunkt kennzeichnet.
B+1
X
a
b
c
n
n
d
d
W
S
N
E
B
P
Abbildung 3.2: Schematische Darstellung der Erzeugung einer vorgebenden Gitterschicht
mit dem orthogonalen Wandabstand des Ausgangsgitters
Ausgehend von den vier Abst¨
anden d1,d2,d3und d4des Punktes P auf den vier Ebenen
BSE, BEN, BNW, BWS kann der kleinste Wert davon als der gesuchte Abstand definiert
werden, d.h. dn=min(d1, d2, d3, d4). Alternativ l¨
asst sich auch der Mittelwert der vier
Abst¨
ande dn= (d1+d2+d3+d4)/4 verwenden. Zur Berechnung des Normalenvektors ~n
wird das Kreuzprodukt benutzt:
~n1=~a ×~
b, ~n2=~
b×~c, ~n3=~c ×~
d, ~n4=~
d×~a, ~n =1
4
4
X
i=1
~ni
|~ni|.
Mit dem Abstand dnund dem berechneten Normalenvektor ~n ergibt sich der gew¨
unschte
Gitterpunkt XB+1 zu:
XB+1 =XB+dn~n . (3.7)
28 Kapitel 3. Numerische Methoden
3.1.3 Approximation der Differentialgleichungen
Die Gleichungen (2.11-2.12) und (2.13-2.14) werden ¨
ublicherweise nach dem Finite-
Differenzenschema approximiert und mit dem iterativ arbeitenden SOR-Verfahren
(Successive-Over-Relaxation) gel¨
ost. Die Konvergenz wird maßgeblich von der Struktur
der zugeordneten Koeffizientenmatrix und damit sowohl von den Koeffizienten der Diffe-
rentialgleichung als auch von der Form der Differenzenapproximation beeinflusst. Neben
der Benutzung zentraler Differenzenquotienten f¨
ur die ersten partiellen Ableitungen wer-
den alternativ auch ”Upwind”-Approximationen verwendet, wodurch sich Konvergenzver-
halten und L¨
osungsqualit¨
at beeinflussen lassen.
Die zentrale Differenzenformel hat folgende Form:
φ,ξ =(φi+1,j,k −φi−1,j,k)
2h.(3.8)
In Abh¨
angigkeit vom Vorzeichen der Kontrollfunktion Pimit positiven αii lautet die vor-
oder r¨
uckw¨
artige Upwind-Approximation:
φ,ξ =φi,j,k −φi−1,j,k
hwenn Pi<0 ;
φ,ξ =φi+1,j,k −φi,j,k
hwenn Pi>0.
(3.9)
Auf der feinsten Gitterebene ist die Maschenweite hgleich eins und wird beim ¨
Ubergang
zur n¨
achst feineren Gitterebene mit zwei multipliziert.
F¨
ur die in den Differentialgleichungen f¨
ur xund Pvorkommenden zweiten sowie gemisch-
ten Ableitungen werden zentrale Differenzenquotienten verwendet, z.B.:
φ,ξξ =(φi+1,j,k −2φi,j,k +φi−1,j,k)
h2,(3.10)
φ,ξη =(φi+1,j+1,k −φi+1,j−1,k −φi−1,j+1,k +φi−1,j−1,k)
4h2.(3.11)
Die so entstandenen Differenzengleichungen lassen sich in der folgenden allgemeinen Form
darstellen:
−apφp+X
nb
anbφnb +Sφ= 0 .(3.12)
Hierbei steht φf¨
ur eine der Variablen xbzw. P. Der Index nb l¨
auft ¨
uber alle Nachbar-
punkte und apund anb sind die Koeffizienten. Der Term Sφbeinhaltet die diskrete Form
der gemischten Ableitung. Alle Koeffizienten sind abh¨
angig von den Metrikkoeffizienten
und damit indirekt auch von der L¨
osung selbst. F¨
ur den zweidimensionalen Fall lautet
die Gleichung (3.12):
−apφi,j +aeφi+1,j +awφi−1,j +anφi,j+1 +asφi,j−1+Sφ= 0 .(3.13)
3.1. Gittergenerierung 29
Hierbei bezeichnen die Indizes w,s,eund ndie West-, S¨
ud-, Ost- und Nordseite des Punk-
tes P. Die Koeffizienten ap,ae,aw,asund anbzw. der Term der gemischten Ableitungen
Sφwerden wie folgt bestimmt:
ae=α11 + 0.5α11h[(1 −Cupw)P1+Cupw(P1+|P1|)] ,
aw=α11 −0.5α11h[(1 −Cupw)P1+Cupw(P1−|P1|)] ,
an=α22 + 0.5α22h[(1 −Cupw)P2+Cupw(P2+|P2|)] ,
as=α22 −0.5α22h[(1 −Cupw)P2+Cupw(P2−|P2|)] ,
ap=ae+aw+an+as,
Sφ=−.5α12(φi+1,j+1 −φi−1,j+1 −φi+1,j−1+φi−1,j−1).
(3.14)
Bei der finiten Approximation wird die urspr¨
ungliche Approximationsgleichung mit h2
multipliziert. P1und P2sind die Kontrollfunktionen des Gitterpunktes P. Der Koeffizient
f¨
ur das Approximationsschema Cupw kann einen beliebigen Wert zwischen 0 und 1 anneh-
men. Setzt man Cupw = 1, wird die UPWIND-Approximation verwendet, mit Cupw = 0
wird die zentrale Differenzenformel benutzt. F¨
ur 0 < Cupw <1 ergibt sich ein gemischtes
Schema.
Gleichung (3.12) l¨
asst sich f¨
ur alle Punkte des Rechengebietes angeben, wobei das resul-
tierende Gleichungssystem in Matrizenschreibweise AX =Blautet. Hier bezeichnet Adie
zugeordnete Koeffizientenmatrix und Bdie rechte Seite. Xist der L¨
osungsvektor, der die
gesuchten kartesischen Koordinaten bzw. Kontrollfunktionen enth¨
alt. Das algebraische
Gleichungssystem wird ¨
ublicherweise iterativ gel¨
ost. Dazu wird hier das SOR-Verfahren
betrachtet, das durch folgenden Algorithmus beschrieben wird:
apφn+1
p=ωÃq
X
nb=1
anbφn+1
nb +
N
X
nb=q+1
anbφn
nb +S(φ)!+ (1 −ω)apφn
p,(3.15)
wobei das Superskript ndie Iterationsstufe des SOR-Verfahrens bezeichnet und ωder
Relaxationsfaktor ist. F¨
ur ω= 1 erh¨
alt man das Gauß-Seidel-Verfahren.
3.1.4 Effektive Behandlung der Randbedingungen
Die L¨
osung der biharmonischen Gleichungen ist empfindlich gegen¨
uber den verwendeten
Randbedingungen, was zu einer schlechten Stabilit¨
at und langsamen Konvergenz f¨
uhren
kann. Deswegen wurden spezielle Verfahren mit der Intention entwickelt, die aus den
Randbedingungen resultierenden Einschr¨
ankungen zu ¨
uberwinden. So haben Brandt und
Dym (1994) z.B. vorgeschlagen, einige zus¨
atzliche Iterationen im randnahen Bereich hin-
zuzuf¨
ugen. Bei dem Verfahren von Linden (1984) wird eine gleichzeitige Relaxation f¨
ur die
Randbedingung mit einem benachbarten inneren Punkt verlangt. Beide Verfahren haben
sich als effektive Techniken zur Behandlung der Randbedingung erwiesen.
Im Verfahren von Linden werden die Kontrollfunktionen Pder Randgitterschicht mit den
Koordinaten xund den Kontrollfunktionen Pder randn¨
achsten Gitterpunkte zu einem
algebraischen Gleichungssystem zusammengefasst, das im zweidimensionalen Fall wie folgt
lautet:
30 Kapitel 3. Numerische Methoden
a11 0a13 a14 0 0
0a22 a23 a24 0 0
0 0 a33 0a35 0
0 0 0 a44 0a46
a51 0 0 0 a55 a56
0a62 0 0 a65 a66
xB+1
1
xB+1
2
PB+1
1
PB+1
2
PB
1
PB
2
=
d1
d2
d3
d4
d5
d6
(3.16)
Dabei werden nur PB
i,xB+1
iund PB+1
ials Unbekannte betrachtet, w¨
ahrend die ¨
ubrigen
Abh¨
angigkeiten von xund Pauf die rechte Seite geschrieben und explizit behandelt wer-
den. Wird das Gleichungssystem mit Hilfe des Gaußschen Eliminationsverfahrens gel¨
ost,
muss wegen der starken ¨
Anderung des Wertes bei den ersten Iterationen unterrelaxiert
werden:
x=xold +ω1(x∗−xold), P =Pold +ω2(P∗−Pold).
Die Ermittlung der Koordinaten xerfolgt schließlich in einem iterativen Prozess:
(1) Startnetz aus algebraischem Verfahren;
(2) Start-Kontrollfunktionen aus algebraischem Netz;
(3) Berechnung der Kontrollfunktionen PB
iam Rand mit Kontrollfunktionen Piund
Koordinaten xider ersten randn¨
achsten Gitterpunkte mit dem Verfahren von Lin-
den;
(4) Gl¨
attung des inneren Gebiets nach Gleichung (3.15);
(5) Iteration von (3) bis (4) bis zur Konvergenz.
Prinzipiell kann ein Netz mit dem in den vorigen Abschnitten beschriebenen Verfahren
blockweise erzeugt und gegl¨
attet werden. Im allgemeinen ist die Gl¨
attung sehr zeitaufwen-
dig, besonders f¨
ur komplexe dreidimensionale Konfigurationen mit großer Gitterpunkt-
zahl. Die Konvergenzgeschwindigkeit kann jedoch durch ein Mehrgitterverfahren erh¨
oht
werden, das im Kapitel 4 dargestellt ist.
3.1.5 Globalgl¨
attung ¨
uber Blockgrenzen
Um komplexe K¨
orper wie z.B. komplette Flugzeuge mit einem strukturierten Gitter ver-
netzen zu k¨
onnen, ben¨
otigt man eine Aufteilung der K¨
orperoberfl¨
ache und des umgeben-
den Raumes in Bl¨
ocke, die untereinander auch unstrukturiert angeordnet werden k¨
onnen.
Im Inneren jedes einzelnen Blockes wird ein strukturiertes Netz erzeugt, wobei die mei-
sten Gittergeneratoren algebraische Verfahren verwenden. Zur Str¨
omungssimulation ist es
vorteilhaft, diese algebraischen Netze im Hinblick auf eine optimale L¨
osungsqualit¨
at zu
gl¨
atten. Im allgemeinen ist der ¨
Ubergang der Gitterlinien ¨
uber Blockgrenzen nicht glatt.
Die Gitterpunkte auf den Blockgrenzfl¨
achen zwischen benachbarten Bl¨
ocken sollen da-
her mit den Gitterpunkten im Blockinneren zusammen gegl¨
attet werden. Dazu muss der
Gl¨
attungsalgorithmus entsprechend erweitert werden, wobei die Indizes der benachbarten
Gitterpunkte in den angrenzenden Bl¨
ocken ben¨
otigt werden.
3.1. Gittergenerierung 31
3.1.6 Datenstruktur und Datenzugriff
Die Gitterpunkte sind innerhalb der einzelnen Gitterbl¨
ocke strukturiert angeordnet. F¨
ur
die unstrukturierte Blockverbindung ist jedoch eine ¨
uber die strukturierte Speicherung
hinausgehende Datenstruktur und ein auf Indexzeigern basierender Datenzugriff erforder-
lich. Der strukturierte Anteil der Daten wird nach Rechenkoordinaten (I, J, K) f¨
ur jeden
Gitterblock M auf einer Gitterstufe L strukturiert gespeichert. Die Indizierung eines Git-
terpunktes ist in Tabelle 3.1 zusammengefasst.
Laufindex Bedeutung Obergrenze
L Gitterstufe NG Anzahl der Gitterstufen
M Gitterblock NB Anzahl der Gitterbl¨
ocke
I Rechenkoordinate NI Dimension in I-Richtung
J Rechenkoordinate NJ Dimension in J-Richtung
K Rechenkoordinate NK Dimension in K-Richtung
Tabelle 3.1: Indizierung des strukturierten Anteils des Datenfeldes im Gittergenerator
Die f¨
unfdimensionale Indizierung (I, J, K, L, M) wird entsprechend der Speicheradresse im
Gitterprogramm einem eindimensionalen Laufindex IJK zugeordnet. Die Gitterpunkte auf
den Blockgrenzfl¨
achen, die dieselben Koordinaten aber unterschiedliche Indizes aufweisen,
werden als Gitterpunktpaar definiert. Als Beispiel ist eine K-Fl¨
ache im zweidimensionalen
Fall mit spezieller Orientierung in Abbildung (3.3) dargestellt. Der Gitterpunkt F mit den
Indizes (I1, J1, K1, N1, M1) und der Gitterpunkt G mit den Indizes (I2, J2, K2, N2, M2)
bilden dabei ein Gitterpunktpaar.
Wenn man die Gitterpunkte F und G wie die Gitterpunkte A, E, I im Gebietsinneren
gl¨
atten m¨
ochte, ben¨
otigt man die Indizes der umgebenden Punkte A, B, C, D, E, H, I, J,
K, L. Sofern sie bekannt sind, k¨
onnen die ersten, zweiten und die gemischten Ableitungen
berechnet werden. Setzt man die Abteilungen, die Koordinaten und Kontrollfunktionen
der umgebenden Gitterpunkte mit bekannten Indizes in die Gleichungen (2.11-2.12) ein,
erh¨
alt man Gleichung (3.12). Mit Gl¨
attung nach (3.15) ergeben sich neue Koordinaten
xF. Da die Gitterpunkte G und F dieselbe Position haben, setzt man xG=xF.
Zur Identifikation der Gitterpunkte auf den benachbarten Blockgrenzfl¨
achen wird ein
f¨
unfdimensionales Feld Paarindex(2, I1, I2, N, M) eingef¨
uhrt. Hierin ist N die Nummer
der Seite im Block M. I1 und I2 sind die lokalen Koordinaten auf der Blockfl¨
ache. Die
Bedeutung der Seitennummer N sowie die Definition der lokalen Koordinaten I1 und I2
sind in der folgenden Tabelle 3.2 dargestellt.
Der Index (I1, I2, N, M) bezeichnet den Gitterpunkt mit den lokalen Koordinaten I1 und
I2 auf der Blockfl¨
ache N von Block M (z.B. Gitterpunkt F in Abbildung 3.3). Im Paar-
index(1, I1, I2, N, M) wird dann der Index des Gitterpunkts G auf der benachbarten
Blockfl¨
ache gespeichert, w¨
ahrend im Paarindex(2, I1, I2, N, M) der Index des Gitter-
32 Kapitel 3. Numerische Methoden




Abbildung 3.3: Gitterpunkte auf den benachbarten Blockgrenzfl¨
achen
N I1 I2 Dimension von I1 Dimension von I2
1 J K NJ(NG, M) NK(NG, M)
2 J K NJ(NG, M) NK(NG, M)
3 K I NK(NG, M) NI(NG, M)
4 K I NK(NG, M) NI(NG, M)
5 I J NI(NG, M) NJ(NG, M)
6 I J NI(NG, M) NJ(NG, M)
Tabelle 3.2: Abh¨
angigkeit der lokalen Koordinaten I1 und I2 von Fl¨
ache N
punkts H, also des zu G benachbarten Gitterpunkts, gespeichert ist.
Aus Abbildung (3.3) geht die Berechnung der Indizes hervor. Die Orientierung der Lau-
findizes IJKF und IJKG f¨
ur die Gitterpunkte F und G k¨
onnen aus den jeweiligen Indizes
berechnet werden. Dann wird der Laufindex IJKH des umgebenen Gitterpunktes H wie
folgt berechnet:
3.1. Gittergenerierung 33
N2 = 1 →IJKH = IJKG + 1
N2 = 2 →IJKH = IJKG −1
N2 = 3 →IJKH = IJKG + NI(NG, M2)
N2 = 4 →IJKH = IJKG −NI(NG, M2)
N2 = 5 →IJKH = IJKG + NI(NG, M2)*NJ(NG, M2)
N2 = 6 →IJKH = IJKG −NI(NG, M2)*NJ(NG, M2)
Sofern IJKG und IJKH bekannt sind, werden sie im Datenfeld Paarindex gespeichert:
Paarindex(1, I1, I2, N1, M1) = IJKG, Paarindex(2, I1, I2, N1, M1) = IJKH.
Hierbei sind I1 und I2 bez¨
uglich N1 des Blocks M1 bekannt. Der Gitterpunkt F bildet ein
Gitterpunktpaar mit G. Gitterpunkt E ist der benachbarte Gitterpunkt von G, dessen
Index durch
N1 = 1 →IJKE = IJKF + 1
N1 = 2 →IJKE = IJKF −1
N1 = 3 →IJKE = IJKF + NI(NG, M1)
N1 = 4 →IJKE = IJKF −NI(NG, M1)
N1 = 5 →IJKE = IJKF + NI(NG, M1)*NJ(NG, M1)
N1 = 6 →IJKE = IJKF −NI(NG, M1)*NJ(NG, M1)
bestimmt wird. Die beiden Laufindizes IJKF und IJKE werden auch im Datenfeld Paarindex
gespeichert:
Paarindex(1, I1, I2, N2, M2) = IJKF, Paarindex(2, I1, I2, N2, M2) = IJKE .
Nicht alle Gitterpunktpaare geh¨
oren zu Gitterpunkten auf den Blockgrenzfl¨
achen, die mit
benachbarten Gitterpunkten zusammen gegl¨
attet werden sollen. Wie in Abbildung (3.4)
zu sehen ist, sind C und D in Abbildung (3.4a) bzw. B und F in Abbildung (3.4b) zwar
Gitterpunktpaare, ihre Position ist aber festgelegt.
Die Gl¨
attung ¨
uber Blockgrenzen wird nur auf dem feinsten Gitter durchgef¨
uhrt. Weil
f¨
ur jede Blockfl¨
ache (N, M) die zur Gl¨
attung ¨
uber Blockgrenzen ben¨
otigten Indizes der
benachbarten Gitterpunkte schon bekannt sind, kann die Gl¨
attung auf den Blockgrenzen
wie im Inneren behandelt werden. Sofern sich xund Pan einem Punkt ¨
andern, m¨
ussen
sie dem anderen Punkt des Gitterpunktpaars ¨
ubergeben werden. Damit ist sichergestellt,
dass beide Gitterpunkte dieselben Koordinaten besitzen.
34 Kapitel 3. Numerische Methoden




















!
!"
"
#
#$
$
%%
%%&
&
'
'(
(
)
)*
*
+
+,
,
--
--.
.
/
/0
0
11
11
22
22
3
34
4
5
56
6
7
78
8
99
99:
:
;;
;;<
<
=
=>
>
??
??
@@
@@
A
AB
B
C
CD
D
EE
EEF
F
G
GH
H
Block M1
Block M2
A
B
C
D
EF
G H
I J
K L
MN
Gitterpunktpaar
Gitterpunkt im Block M1
Gitterpunkt im Block M2
(a)










Block M2
Block M1
H
D
C
G
B
F
E
A
(b)
Abbildung 3.4: Schematische Darstellung der Gitterpunktpaare
3.2. Str¨
omungssimulation 35
3.2 Str¨
omungssimulation
Zur Str¨
omungssimulation werden die von Rung u. a. (1999) und Xue (1998) entwickelten
zweidimensionalen und dreidimensionalen Finite-Volumen Navier-Stokes Verfahren einge-
setzt. Beide benutzen einen Algorithmus, der auf der streng konservativen Formulierung
der Grundgleichungen f¨
ur primitive Variablen in Bezug auf allgemein-krummlinige Koor-
dinaten basiert. Zur effizienten Erfassung komplexer Geometrien wird das Integrationsge-
biet durch ein blockstrukturiertes Rechengitter diskretisiert, dessen Verbindungsfl¨
achen
unstrukturiert sein k¨
onnen.
Die in den Integro-Differentialgleichungen auftretenden Integrale werden durch die von
zweiter Ordnung genaue Mittelpunktregel angen¨
ahert. Zur Approximation konvektiver
Terme werden die Variablenwerte in den Fl¨
achenschwerpunkten durch stromauf gerich-
tete, monotone TVD-Ans¨
atze (Harten 1983) interpoliert. Diese sind von maximal drit-
ter Genauigkeitsordnung und erm¨
oglichen eine oszillationsfreie Darstellung der f¨
ur ho-
he Reynoldszahlen dominanten konvektiven Prozesse bei gleichzeitiger Beachtung ihrer
physikalischen Orientierung. Aus Stabilit¨
atsgr¨
unden st¨
utzt sich die Implementierung der
TVD-Approximation auf die deferred-correction Technik (Lien und Leschziner 1994c). Die
Untersuchungen von Lilek und Peri´c (1995) belegen, dass die hybride Approximations-
technik einen optimalen Kompromiss zwischen Effizienz und Genauigkeit darstellt.
Das Verfahren benutzt eine voll implizite zeitliche Integration zweiter Ordnung. Die Be-
stimmung eines Druckfeldes, das vertr¨
aglich mit der Kontinuit¨
atsgleichung ist, basiert auf
einem kompressiblen Druckkorrekturalgorithmus (Demirdzi´c u. a. 1993). Alle Str¨
omungs-
gr¨
oßen sind in den Zellzentren gespeichert. Die enge Kopplung von Druck, Geschwindigkeit
und Reynoldsspannungen wird durch eine verallgemeinerte Rhie & Chow Interpolation
(Rhie und Chow 1983) auf Grundlage des apparent-pressure / apparent-viscosity Prinzips
(Obi u. a. 1991) sichergestellt.
Die Flamelet-Bibliotheken werden durch die in dem CFD-Code implementierte modifizier-
te Version des Programms RIF (Representative Interactive Flamelets) bestimmt, das vom
Institut f¨
ur Technische Mechanik, RWTH Aachen zur Verf¨
ugung gestellt wurde (Peters
2000).
3.2.1 Finite Approximation der Bilanzgleichungen
Die Finite-Volumen-Approximation wird anhand einer allgemeinen Bilanzgleichung f¨
ur
eine generische Gr¨
oße φdurchgef¨
uhrt:
∂ ρ φ
∂ t +∇·(ρ u φ)−∇·(Γφ∇φ) = Sφ.(3.17)
Hierin kennzeichnet Γφden Diffusionskoeffizienten f¨
ur die Gr¨
oße φ. Die Quellterme und
diejenigen Terme der Bilanzgleichung, die nicht durch den Gradientenfluß −∇ · Γφ∇φ
dargestellt werden k¨
onnen, sind in Sφzusammengefasst.
36 Kapitel 3. Numerische Methoden
W
E
S
N
T
B
we
b
t
s
n
P
δΑe
Abbildung 3.5: Bilanzvolumen und Kontrollvolumenfl¨
achen eines strukturierten Gitters,
w: west; e: east; s: south; n: north, b: bottom, t: top
F¨
ur ein Kontrollvolumen, wie in Abbildung 3.5 skizziert, f¨
uhrt die numerische Approxi-
mation auf ein algebraisches Gleichungssystem f¨
ur die Variable φP:
AP
PφP−X
f
ANf φNf =SP.(3.18)
Hierin bezeichnen hochgestellte Indizes den Ort, an dem die Gr¨
oße ausgewertet werden soll
und der untere Index der Koeffizienten markiert die Fl¨
ache, an der die Bilanz ausgef¨
uhrt
werden soll. Die Koeffizientenmatrix enth¨
alt die Hauptdiagonale APund die Nebendia-
gonalen −ANf . Da die Variablen nur an den Volumenzentren vorhanden sind, muss das
Diskretisierungsschema die Werte auf den Fl¨
achen der Kontrollvolumen mit den Werten
in den Volumenzentren in Verbindung setzen. In SP
φsind neben dem urspr¨
unglichen Quell-
term alle Anteile aus Diffusion und Konvektion enthalten, die explizit behandelt werden.
Auf das Diskretisierungsschema wird nachfolgend genauer eingegangen, so dass hier der
Index Nf, siehe Tabelle 3.3, die Beziehung zwischen dem Wert auf der Oberfl¨
ache und
dem benachbarten Volumenzentrum veranschaulicht.
f w e s n b t
Nf W E S N B T
Tabelle 3.3: Zuordnung der Fl¨
achenwerte zu Volumenzentren
3.2. Str¨
omungssimulation 37
Die Berechnung der Koeffizienten ergibt sich aus der Integration der Transportgleichung
(3.17) unter Ber¨
ucksichtigung des Integralsatzes von Gauß:
ZδV
∂ ρφ
∂ t dV +ZδA
(ρ u φ −Γ∇φ)·dA =Sφ,(3.19)
wobei δA die Oberfl¨
ache des Volumens bezeichnet und dA das aus dem Volumen heraus-
zeigende Fl¨
achenelement der Oberfl¨
ache ist. F¨
uhrt man die Integration f¨
ur ein Kontrollvo-
lumen mit dem Mittelwertsatz der Integralrechnung aus, wird nur ein Punkt im Volumen
bzw. auf den Fl¨
achen ben¨
otigt:
µδV ∂ ρφ
∂ t P
+X
f
δAf·(ρ u φ −Γ∇φ)f=SP
φ.(3.20)
Die Summation muss ¨
uber alle Teiloberfl¨
achen des Kontrollvolumens um Perfolgen. δV
gibt den Volumeninhalt des Kontrollvolumens an und δA ist der nach außen gerichtete
Fl¨
achenvektor der Teiloberfl¨
ache f. Mit der diskreten Form der Massenbilanz, d.h. φ= 1,
q= 0 und Γ = 0:
µδV ∂ ρ
∂ t P
+X
f
Ff
m= 0 ,(3.21)
in der der Massenfluß durch Fm=δA ·ρu definiert ist, lautet die diskrete Transportglei-
chung φwie folgt:
µδV ρ ∂ φ
∂ t PX
f
FfφP+X
f
Ffφf−X
f
[δA ·(Γ ∇φ)]f=SP
φ.(3.22)
Der in der Impulsbilanz enthaltene Reynoldssche Spannungsterm h¨
angt in starkem Maße
von den Geschwindigkeitsgradienten ab. Die lineare Abh¨
angigkeit zwischen Reynolds-
spannungen und Geschwindigkeitsgradienten werden nach der Mischungsweghypothese
mit dem molekularen Spannungsterm zusammengefasst:
πg=−pgI+µgS−2
3µg(∇·u)I+τtn , µg=µ+µt, pg=p+2
3ρk ,
wobei τtn den nichtlinearen Anteil der Reynoldsspannungen aus dem EASM bezeichnet.
Bei der Umsetzung der Reynoldsschen Gleichungen wird der Term µg∇uimplizit behan-
delt. Alle anderen Terme werden explizit bestimmt und kommen auf die rechte Seite der
Gleichung. Bei der Bilanzierung eines Kontrollvolumens lautet der explizite Teil:
−δA ·pg+2
3µg(∇·u)¸f
−δA ·³µ∇Tu+τtn´f.
F¨
ur die zu l¨
osenden str¨
omungsmechanischen Bilanzgleichungen m¨
ussen auf den Berandun-
gen des Rechengebiets Randbedingungen vorgegeben werden. Da es sich bei dem Finite-
Volumen-Verfahren um die L¨
osung von Integro-Differentialgleichungen handelt, wird die
Fl¨
achenintegration dieser Gleichungen an Randfl¨
achen durchgef¨
uhrt und an den rand-
n¨
achsten Bilanzvolumen teils implizit, teils explizit ber¨
ucksichtigt.
38 Kapitel 3. Numerische Methoden
3.2.2 Modifizierte Behandlung der Wandrandbedingungen
Bei stark schiefwinkligem Gitterlinienverlauf kann die Implementierung der Wandrandbe-
dingungen die L¨
osungsqualit¨
at des Gesamtverfahrens erheblich beeinflussen. Im Rahmen
dieser Arbeit wird deshalb ein Konzept weiter entwickelt, das bei stark schiefwinkligem
Gitterlinienverlauf verwendet werden kann.
R
d
sP
A= A nδ
δ
δ
Pdn
(a)
δΑ=δΑ
δS
PRn
d
n
P
d
(b)
Abbildung 3.6: Veranschaulichung eines Kontrollvolumens mit physikalischer R¨
andfl¨
ache
R
Um eine genaue Approximation an der Wand zu erreichen, wird zun¨
achst die Geschwin-
digkeit an einem Punkt P0berechnet, der sich im Abstand lnsenkrecht ¨
uber dem Zentrum
der Randfl¨
ache R befindet, wie Abbildung 3.6a zeigt. Der Abstandsvektor d0ist die Pro-
jektion des Abstandsvektor din Fl¨
achennormalenrichtung:
d0=lnn , ln=d·n .
Der Abstandsvektor vom Finite-Volumen-Zentrum zum Punkt P0wird mit δs bezeichnet
und ist Null bei rechtwinkligen Gittern:
δs =d−d0.
Die Geschwindigkeit am Punkt P0wird mit explizit berechneten Geschwindigkeitsgradi-
enten bestimmt:
uP0=uP+δs ·(5u)P.(3.23)
3.2. Str¨
omungssimulation 39
An der Wand gilt als physikalische Randbedingung der Impulsgleichungen die Wandhaft-
bedingung u= 0. Damit l¨
asst sich die von der Wand auf das Finite-Volumen wirkende
Fl¨
achenkraft durch
δA ·τg=−pδA +µgδA∂us
∂n ≈ −pRδA −µgδA
ln
uP0
s(3.24)
berechnen. Mit uswird die Komponente des Geschwindigkeitsvektors parallel zur Wand
bezeichnet.
Diese Behandlung ergibt einen relativ großen Fehler bei stark schiefwinkligem Gitter-
linienverlauf an der Wand. Wie in Abbildung 3.6b zu erkennen ist, kann der Punkt P’
weit außerhalb des Volumens liegen. G¨
unstig ist es, hier im Rahmen eines Kompromisses
zwischen Genauigkeit und Stabilit¨
at eine modifizierte Variante zu nutzen. Bei solch stark
schiefwinkligem Gitterlinienverlauf an der Wand wird die Fl¨
achenkraft direkt mit der
Geschwindigkeit des FV-Zentrums durch
δA ·τg=−pδA +µgδA∂us
∂n ≈ −pRδA −µgδA
ln
uP
s(3.25)
gebildet. Er ergibt sich aus der Differenz zwischen der Gesamtgeschwindigkeit uund dem
zur Wand senkrecht stehenden Teil des Geschwindigkeitsvektors un:
un=unn= (u·n)n . (3.26)
Der wandparallele Anteil der Geschwindigkeit am Punkt Pist:
uP
s=uP−(uP·n)n . (3.27)
Die Fl¨
achenkraftkomponente der Impulsgleichung f¨
ur die Geschwindigkeitskomponente
uαist dann mit Cw=δAµg/lngegeben:
−pRδAα−Cw(1 −n2
α)uP
α
| {z }
Term1
−pRδAα+Cwnα(bβuP
β+nγuP
γ)
| {z }
Term2
.(3.28)
F¨
ur die zyklisch laufenden Indizes α,βund γgilt keine Summationskonvention. Term1
mit uP
αwird implizit und Term2 explizit behandelt. Bei einer laminaren Str¨
omung ist
die Gesamtviskosit¨
at µggleich der dynamischen Viskosit¨
at µdes Fluids, w¨
ahrend bei
turbulenter Str¨
omung folgende Fallunterscheidung gilt:
(a) Das wandn¨
achste Finite-Volumen-Zentrum Pbefindet sich außerhalb der z¨
ahen Un-
terschicht. Hier wird das logarithmische Wandgesetz (High-Re Hypothese)
u
uτ
=1
κlogeY++B , Y +=lnuτ
µ, κ = 0.41, B = 5.0,(3.29)
herangezogen und daraus eine Beziehung zwischen der Gesamtviskosit¨
at µgund der
turbulenten kinetischen Energie khergeleitet:
µg=kρlnµτ
loge(βY +)mit uτ=C1/4
µ√k , β = exp(κB).(3.30)
Diese Beziehung gilt f¨
ur Y+>11.
40 Kapitel 3. Numerische Methoden
(b) Das wandn¨
achste Finite-Volumen-Zentrum Pbefindet sich in der z¨
ahen Unterschicht
mit Y+<11:
µg=µ . (3.31)
An der Wand sind die turbulente kinetische Energie kund ihr Gradient senkrecht zur
Wand Null, so dass die Konvektion und die Fl¨
achenintegration ¨
uber den Diffusionsterm
keinen Beitrag liefern. Bei der turbulenten Dissipationsrate ²bzw. der spezifischen Dis-
sipationsrate ωwird f¨
ur das wandn¨
achste (wn) Finite-Volumen keine Bilanzgleichung
gel¨
ost. Vielmehr werden dort die Werte von ²wn in Abh¨
angigkeit vom Wandnormallenab-
stand lnund der lokalen turbulenten kinetischen Energie kwn im Finiten-Volumen direkt
vorgegeben:
²wn =k3/2
wn
κC−3/4
µln
.(3.32)
Diese Werte lassen sich durch
²wn =2νkwn
l2
n
(3.33)
bestimmen, wenn sich die wandn¨
achste Gitterzelle innerhalb der z¨
ahen Unterschicht be-
findet.
Beide o.a. Randbedingungstechniken sind nicht universell einsetzbar, weil sie bestimm-
te Str¨
omungsverh¨
altnisse in den wandnahen Gitterschichten voraussetzen. Insbesondere
muss bei der Gittergenerierung die Dicke der sich ausbildenden semi-viskosen Str¨
omungs-
bereiche betrachtet werden, worin der entscheidende Nachteil beider Methoden besteht.
Die Ausdehnung der Wandturbulenz h¨
angt nicht nur von dem zu simulierenden Str¨
omungs-
feld ab, sondern variiert erfahrungsgem¨
aß auch stark entlang von K¨
orperkonturen. Nahe-
zu diskontinuierliche Variationen k¨
onnen z.B. durch Stoß/Grenzschicht Wechselwirkungen
auftreten. Eine lokale Verletzung des G¨
ultigkeitsbereichs der verwendeten Randbedingung
kann mit erheblichen globalen Nachteilen f¨
ur die Stabilit¨
at und Qualit¨
at der Simulation
verbunden sein. Ferner ist mit starken Einschr¨
ankungen bei der Verwendung von Mehrgit-
tertechniken, deren Gitter teilweise stark unterschiedliche Aufl¨
osungsverm¨
ogen besitzen,
zu rechnen. Eine adaptive Randbedingung f¨
ur die gemeinsame Behandlung von Low- und
High-Re Bereichen wurde von Rung (1999) entwickelt und wird hier bei der Str¨
omungs-
simulation verwendet.
Kapitel 4
Grundlage des Mehrgitterverfahrens
Als eine allgemeine Methode f¨
ur die Beschleunigung iterativer L¨
osungsverfahren wurden
Mehrgitterverfahren schon auf eine ganze Reihe praktisch relevanter Probleme angewandt,
so etwa die elliptische Gittergenerierung sowie die L¨
osung der Navier-Stokes Gleichungen
oder der Euler Gleichungen. In Rahmen dieser Arbeit konzentriert sich die Darstellung
der Methode auf die Anwendung in der Gittergenerierung und Str¨
omungssimulation. Das
Prinzip des Mehrgitterverfahrens wird grob skizziert. F¨
ur eine ausf¨
uhrlichere Darstellung
sei auf Hackbusch (1985) und Wesseling (1992) verwiesen. Die Mehrgitterkomponenten
werden in den folgenden Kapiteln beschrieben.
4.1 Grundprinzip
Die finite Approximation einer beliebigen linearen oder nichtlinearen Differentialgleichung
auf einem Gitter, gekennzeichnet mit dem hochgestellten Index h, f¨
uhrt zu einem Glei-
chungssystem vom Typ
L(φh) = fh,(4.1)
wobei Leinen Operator darstellt, der alle funktionalen Zusammenh¨
ange der Unbekannten
φhrepr¨
asentiert, und fhist die von φhunabh¨
angige rechte Seite des Gleichungssystems.
Bei einem nichtlinearen Verfahren ist Lvon φabh¨
angig. Dieses System wird normalerweise
iterativ gel¨
ost. Ausgehend von einer Startsch¨
atzung erh¨
alt man nach einigen Iteration mit
einem Iterationsverfahren eine N¨
aherungsl¨
osung ˜
φh, die das Gleichungssystem (4.1) bis
auf ein Residuum Rherf¨
ullt:
L(˜
φh) = fh+Rh.(4.2)
Bez¨
uglich der exakten L¨
osung φhist die N¨
aherungsl¨
osung ˜
φhmit einen Fehler δφhbehaftet,
so dass die folgende Gleichung erf¨
ullt ist:
L(˜
φh+δφh) = fh.(4.3)
Subtraktion dieser Gleichung von Gleichung (4.2) liefert die Feingitterfehlergleichung:
L(˜
φh+δφh) = L(˜
φh)−Rh.(4.4)
41
42 Kapitel 4. Grundlage des Mehrgitterverfahrens
Wenn Lein linearer Operator ist, vereinfacht sich die Bestimmungsgleichung f¨
ur δφh:
L(δφh) = −Rh.(4.5)
Es kann nun festgestellt werden, dass die meisten iterativen L¨
osungsverfahren in der Lage
sind, hochfrequente Fehler, d.h. eine Verteilung des Residuums im Gitter, die sich von
einer Gitterzelle zur n¨
achsten stark ¨
andert, sehr schnell zu reduzieren. Es stellt sich dann
eine relativ glatte niederfrequente Fehlerverteilung ein, die von den gebr¨
auchlichen ite-
rativen L¨
osern w¨
ahrend des Iterationsprozesses nur langsam reduziert wird. Man spricht
deshalb auch von den Gl¨
attungseigenschaften eines iterativen L¨
osers bzw. eines Gl¨
atters.
Die Mehrgittermethode macht sich diese Eigenschaft iterativer L¨
oser zunutze. Bei ihr wer-
den L¨
osungen des Gleichungssystems auf unterschiedlich feinen Gittern bestimmt. Man
berechnet zun¨
achst eine N¨
aherung der L¨
osung auf dem feinsten Gitter. Es stellt sich
schnell eine glatte Fehlerverteilung ein. Auf einem gr¨
oberen Gitter wird dann die Fehler-
verteilung dieser L¨
osung berechnet. Dies ist aus zweierlei Gr¨
unden auf dem groben Gitter
einfacher. Zum einen hat eine niederfrequente Fehlerverteilung auf einem feinen Gitter in
einem groben Gitter bereits eine h¨
ohere Fehlerverteilungsfrequenz. Zum anderen bestehen
die korrespondierenden Gleichungssysteme auf dem gr¨
oberen Gitter aus sehr viel weniger
Unbekannten als auf dem feinen Gitter.
Der Mehrgitter-Berechnungsalgorithmus setzt sich aus folgenden Schnitten zusammen:
(1) Vorgl¨
attung: Auf dem Feingitter wird eine N¨
aherungsl¨
osung ˜
φhdurch ν1Iteratio-
nen erzeugt. Sie hat bez¨
uglich der exakten L¨
osung φhdes Gleichungssystems einen
Fehler:
δφh=φh−˜
φh.(4.6)
Diese Fehlerverteilung kann auf dem Grobgitter, das mit dem hochgestellten Index
Hgekennzeichnet ist, berechnet werden. Dazu wird das lokale Residuum ben¨
otigt,
das folgendermaßen definiert ist:
Rh=Lh(˜
φh)−fh.(4.7)
(2) Berechnung der Grobgitterfehlerverteilung: Es gibt zwei Standardverfahren: FAS
und CS (Correction Scheme). Beim nichtlinearen FAS-Verfahren werden die N¨
ahe-
rungsl¨
osung ˜
φhund das Residuum Rhdes Feingitters auf das Grobgitter Hrestrin-
giert. Im folgenden bezeichnet [IH
h] den Restriktionsoperator f¨
ur die Gittervariablen,
[ˆ
IH
h] bezeichnet den Restriktionsoperator f¨
ur das Residuum:
φH= [IH
h]˜
φh, RH= [ˆ
IH
h]Rh.(4.8)
Die Grobgitterl¨
osung ˜
φHwird durch
LH(˜
φH) = LH(φH)−[ˆ
IH
h]Rh(4.9)
4.2. Gittervergr¨
oberung und Gittertransfer 43
berechnet, wobei φHund ˜
φHdie Startl¨
osung bzw. die N¨
ahrungsl¨
osung auf dem
Grobgitter sind. Daraus wird der Grobgitterfehler bestimmt:
δφH=˜
φH−φH.(4.10)
Beim linearen CS-Verfahren kann der Grobgitterfehler direkt berechnet werden
durch:
LH(δφH) = −[ˆ
IH
h]Rh.(4.11)
(3) Prolongation der Fehlerverteilung: Nachdem die Fehlerverteilung auf dem Grob-
gitter bestimmt wurde, kann sie auf das feine Gitter transferiert und damit die
Feingitterl¨
osung korrigiert werden:
φh=˜
φh+ [Ih
H]δφH.(4.12)
Hier bezeichnet [Ih
H] den Prolongationsoperator vom groben auf das feine Gitter.
(4) Nachgl¨
attung: Die so korrigierte L¨
osung φhauf dem feinen Gitter erf¨
ullt im all-
gemeinen das Gleichungssystem noch nicht, da durch die Prolongation neue Fehler
eingebracht werden. Das Residuum ist allerdings hochfrequent und kann daher leicht
durch einige zus¨
atzliche Iterationszyklen reduziert werden.
Es sollte erw¨
ahnt werden, dass der Fehler auch oft als Korrektur bezeichnet wird. Im
folgenden wird deshalb ebenfalls der Begriff Korrektur verwendet.
4.2 Gittervergr¨
oberung und Gittertransfer
In dem hier vorgestellten Ansatz ben¨
otigt man ein Verfahren zur Gittervergr¨
oberung, das
f¨
ur die Gittergl¨
attung und die Str¨
omungssimulation in unterschiedlicher Weise gebildet
wird. Die bisherigen Betrachtungen sind auch weitgehend unabh¨
angig von der verwen-
deten Diskretisierungsmethode. Mehrgitterverfahren lassen sich f¨
ur Finite-Differenzen,
Finite-Volumen und Finite-Elemente-Verfahren in nahezu analoger Weise definieren, wo-
bei jedoch insbesondere bei der Interpolation und Restriktion auf die jeweilige Besonder-
heit der Diskretisierung R¨
ucksicht genommen werden muss. Das Verfahren zur Gitter-
vergr¨
oberung und die Gittertransferoperationen werden im folgenden beispielhaft f¨
ur die
Gittergenerierung und die Str¨
omungssimulation basierend auf Differenzen bzw. Finite-
Volumen beschrieben.
44 Kapitel 4. Grundlage des Mehrgitterverfahrens
4.3 Mehrgitterzyklen
Da die L¨
osung auf dem Grobgitter normalerweise ebenfalls iterativ bestimmt werden muss,
wird das oben beschriebene Verfahren rekursiv angewendet, d.h. die L¨
osung auf dem Grob-
gitter kann berechnet werden, indem zun¨
achst eine N¨
aherungsl¨
osung bestimmt und diese
nach dem gleichen Verfahren mit dem Fehler eines noch gr¨
oberen Gitters verbessert wird.
Wenn auf jeder Gitterebene das Zweigitterverfahren nur einmal angewendet wird, erzeugt
es einen sogenannten V–Zyklus. Bei zweimaliger rekursiver Anwendung des Zweigitter-
verfahrens wird das entstandene Mehrgitterverfahren W–Zyklus genannt. In Abbildung
4.1 sind V– und W–Zyklus schematisch f¨
ur ein Mehrgitterverfahren mit drei Gitterebe-
nen dargestellt. Im allgemeinen werden jedoch V-Zyklen verwendet, da sie bez¨
uglich der
Stabilit¨
at ausreichend sind.
3
ν
ν
ν
ν
ν
1
12
2
0
ν ν
ν
ν
νν
ν1
1
0
2
0
2
2
grob
fein
1
2
Vor/Nach-Glättung
Glättung auf dem gröbsten Gitter
W-Zyklus
V-Zyklus
Prolongation
Restriktion
Gitterebene
Abbildung 4.1: Schematischer Ablauf eines V–Zyklus und W– Zyklus
In Abbildung 4.1 symbolisiert ein abw¨
arts gerichteter Pfeil den ¨
Ubergang vom Fein– zum
Grobgitter und beinhaltet die Berechnung der Residuen Rhauf dem Feingitter, die Re-
striktionen und die Berechnung der rechten Seite fHder Grobgittergleichung. Das Symbol
%verk¨
orpert den ¨
Ubergang vom Grob– zum Feingitter und beinhaltet die Berechnung
des Fehlers auf dem Grobgitter, seine Prolongation und die L¨
osungskorrektur auf dem
Feingitter.
Die Konvergenz des Verfahrens wird verbessert, indem man mehrere Zyklen hinterein-
ander durchf¨
uhrt, d.h. der Zyklus aus Berechnung der Feingitterl¨
osung (Vorgl¨
attung),
Restriktion, Berechnung der Grobgitterkorrektur, Prolongation und Nachgl¨
attung wird
mehrfach wiederholt, bis das Residuum auf dem feinsten Gitter eine vorgegebene Grenze
unterschreitet.
4.4. FMG-Methode 45
4.4 FMG-Methode
Zur L¨
osung der Navier-Stokes Gleichungen wird oft die sogenannte FMG-Methode ver-
wendet. Hier wird bei der Rechnung nicht vom feinsten Gitter ausgegangen und dann
auf das Grobgitter restringiert, sondern die L¨
osung wird zuerst auf dem gr¨
obsten Gitter
ermittelt. Durch Prolongation der Grobgitterl¨
osung ergibt sich dann eine m¨
oglichst gute
Startl¨
osung f¨
ur das n¨
achst feinere Gitter. Dort bestimmt man eine N¨
aherungsl¨
osung und
¨
ubertr¨
agt diese gemeinsam mit dem Residuum zur¨
uck auf das grobe Gitter. Das Ergeb-
nis einer Grobgitterkorrektur wird nun auf das feine Gitter prolongiert, um die Feingit-
terl¨
osung zu korrigieren. Dieser Algorithmus ist rekursiv. Man beginnt auf dem gr¨
obsten
Gitter und berechnet dann sukzessive mit FAS oder CS die L¨
osung auf mehreren Gitter-
ebenen bis die L¨
osung auf dem feinsten Gitter gefunden wird. F¨
ur ein FMG-Verfahren
mit vier Gitterebenen und V-Zyklen ergibt sich das in Abbildung 4.2 gezeigte Flussdia-
gramm. Da die FMG-Methode algorithmisch leicht implementiert werden kann, wird sie
h¨
aufig eingesetzt.
Konvergierte Lösung Glättung auf dem gröbsten Gitter
Vor/Nach−Glättung
Abbildung 4.2: Schematischer Ablauf eines Full Multigrid-Verfahrens
Bei der Zyklenwahl muss noch festgelegt werden, wie viel Vor- und Nachgl¨
attungen durch-
gef¨
uhrt werden sollen.
Startl¨
osung
Beim FAS-Algorithmus wird die Restriktion der N¨
aherungsl¨
osung ˜
φhauf dem feinen Gitter
als Startl¨
osung φHauf dem groben Gitter benutzt. Im Prinzip gibt es mehrere M¨
oglich-
keiten, die Startl¨
osung φHzu w¨
ahlen. Als eine Alternative kann z.B. die N¨
aherungsl¨
osung
des letzten Zyklus auf dem Grobgitter als Startl¨
osung φHverwendet und der rechtseitige
Term L(φH) in Gleichung (4.9) daraus berechnet werden. Dieses Konzept f¨
uhrt zu einem
modifizierten Mehrgitterverfahren f¨
ur die Navier-Stokes Gleichungen, welches im Kapitel
6 n¨
aher beschrieben wird.
46 Kapitel 4. Grundlage des Mehrgitterverfahrens
Kapitel 5
Gittergenerierung und Gl¨
attung mit
Mehrgitterverfahren
Dieses Kapitel befasst sich mit der Gittergenerierung und der Gl¨
attung mit dem Mehr-
gitterverfahren. Das Prinzip des Mehrgitterverfahrens in der Gittergenerierung wurden
schon von St¨
uben und Linden (1984) angegeben. Obwohl in der Vergangenheit Mehrgit-
terverfahren erfolgreich zur L¨
osung elliptischer Gleichungen angewendet wurden, sind nur
wenige Anwendungen in der Gittergl¨
attung zu finden. Diese betreffen meistens einfache
Konfigurationen (Jain 1986; Spitaleri 1990; Spitaleri 1991). Hier wurde nur der nicht-
lineare FAS-Algorithmus verwendet, f¨
ur den linearen CS-Algorithmus findet man keine
Anwendung. Im Rahmen einer Fortf¨
uhrung des aktuellen Kenntnisstandes wird in dieser
Arbeit nicht nur FAS sondern auch CS zur Gittergl¨
attung verwendet. Die Gl¨
attung ba-
siert auf dem Verfahren mittels biharmonischen Gleichungen, das hier entwickelt wurde.
Das Verhalten der beiden Algorithmen wird mit V– und W–Zyklus untersucht.
In diesem Kapitel werden zun¨
achst die Komponenten des Mehrgitterverfahrens in der
Gittergenerierung vorgestellt und dann der Rechenablauf im Detail beschrieben. Zur Va-
lidierung werden einige komplexe zwei- und dreidimensionale Konfigurationen ausgew¨
ahlt.
Damit wird in dieser Arbeit erstmals ein Konzept zur differentiellen Gl¨
attung f¨
ur praktisch
relevanten Konfigurationen angewendet.
5.1 Komponenten von Mehrgitterverfahren
Wie in Abbildung 4.1 gezeigt ist, bestehen beim Einsatz von Mehrgitterverfahren eine
Reihe von M¨
oglichkeiten: es muss festgelegt werden, welches Gl¨
attungsverfahren und wel-
che Interpolationen und Restriktionen verwendet werden. Desweiteren muss die Art der
Gittervergr¨
oberung gew¨
ahlt werden und die Anzahl der Vor- und Nachgl¨
attungen ist zu
bestimmen. Weiterhin muss entschieden werden, wann die Gittervariablen auf das n¨
achst
gr¨
obere oder n¨
achst feinere Gitter interpoliert werden. Wie in Kapitel 4 vorgestellt wur-
de, wird das in der Gittergenerierung ¨
ublicherweise verwendete SOR-Verfahren als das
Gl¨
attungsverfahren benutzt.
47
48 Kapitel 5. Gittergenerierung und Gl¨
attung mit Mehrgitterverfahren
Gittervergr¨
oberung: Weil die Differentialgleichungen im krummlinigen Koordinaten-
system diskretisiert werden, kann die gleichf¨
ormige Gittervergr¨
oberung mit einer Verdop-
pelung der Gittermaschenweite, d.h. H= 2h, durchgef¨
uhrt werden. Hierbei ist hdie
Maschenweite des n¨
achst feineren Gitters und es gilt h= 1 auf dem feinsten Gitter. Die
Grobgitterpunkte fallen mit den Feingitterpunkten zusammen.
Restriktion: Bei der Restriktion werden die Residuen und beim FAS-Verfahren auch
die Feingitterl¨
osung (Koordinaten und Kontrollfunktionen) von dem Fein- auf das Grob-
gitter mit Hilfe des Operators [IH
h]¨
ubergeben. Die einfachste Form der Restriktion ist die
sogenannte Injektion, wobei sowohl das Residuum als auch die L¨
osung direkt vom kor-
respondierenden Feingitterpunkt ¨
ubernommen werden. Restriktionsvorschriften h¨
oherer
Ordnung beziehen weitere Umgebungspunkte dieser Feingitterpunkte mit ein. In der vor-
liegenden Arbeit wird das Full-Weighting Schema verwendet, da diese Vorgehensweise von
der Stabilit¨
at her gesehen als die sicherste gilt. Das Schema wird im zweidimensionalen
Fall wie folgt definiert:
[ˆ
IH
h] = 1
16
121
242
121
.
Prolongation: Die Prolongation dient zum Transfer der berechneten Korrektur von dem
Grobgitter auf das Feingitter. Die bilineare Interpolationsformel gilt f¨
ur zweidimensionale,
die trilineare Interpolationsformel f¨
ur dreidimensionale F¨
alle. Im zweidimensionalen Fall
berechnet sich f¨
ur den Feingitterpunkt der interpolierte Wert auf den Grobgitterpunkt
nach:
δφh
2i+l,2j+m=1
4(δφH
i,j +δφH
i+l,j +δφH
i,j+m+δφH
i+l,j+m),(5.1)
mit l= 0,1 und m= 0,1.
5.2 Rechenablauf des Mehrgitterverfahrens
Das Mehrgitterverfahren wird nur zur Berechnung der Koordinaten verwendet. In der
vorliegenden Arbeit werden das CS- und das FAS-Verfahren untersucht. Der prinzipielle
Ablauf des Mehrgitterverfahrens wurde bereits in Kapitel 4 dargestellt. Als Beispiel wer-
den hier das verwendete Gl¨
attungsverfahren, der ¨
Ubergang vom Fein– zum Grobgitter
und vom Grob– zum Feingitter zusammenfassend beschrieben.
Das diskretisierte linearisierte Gleichungssystem auf einer Gitterebene mit Maschenweite
hkann wie folgt geschrieben werden:
−ah
pφh
p+X
nb
ah
nbφh
nb +Sh
φ=fh.(5.2)
Hierbei stellt φeine der Variablen xoder Pauf dem feinsten Gitter, die Variable x
beim FAS-Verfahren oder die Gitterkorrektur δx beim CS-Verfahren auf dem Grobgitter,
5.2. Rechenablauf des Mehrgitterverfahrens 49
dar. Die rechte Seite fφist auf dem feinsten Gitter Null und auf dem groben Gitter von
Null verschieden. Auf dem groben Gitter bleibt sie konstant f¨
ur die Grobgitterkorrek-
turrechnung. Sie wird sp¨
ater genau angegeben. Zur L¨
osung der Gleichung (5.2) wird das
SOR-Verfahren verwendet:
an
pφn+1
p= (1 −ω)an
p)φn
p+ωÃq
X
nb=1
anbφn+1
nb +
N
X
nb=q+1
anbφn
nb +Sh
φ−fh!.(5.3)
Auf dem feinsten Gitter wird zuerst die im Kapitel 4 beschriebene Behandlung der Rand-
bedingungen f¨
ur alle Randgitterpunkte und ihre randn¨
achsten Nachbarn angewendet und
danach die Gl¨
attung f¨
ur die inneren Gitterpunkte durchgef¨
uhrt. Die Metrikkoeffizienten,
die Koeffizienten ah
nb und Sh
φin Gleichung (5.2) werden nach jeder Iteration aktualisiert.
Auf dem Grobgitter ist die Gl¨
attung nur im Gebietsinneren durchzuf¨
uhren. Zuerst wer-
den die Koordinaten und Kontrollfunktionen auf dem groben Gitter vom feinen Gitter
restringiert. Weil das Mehrgitterverfahren nur zur Berechnungen der Koordinaten ver-
wendet wird, bleiben die Kontrollfunktionen auf allen Grobgittern konstant. Hier wird
das DCA-Verfahren (Discretization Coarse Grid Approximation) (Hackbusch 1985) ver-
wendet, d.h. die Metrikkoeffizienten, die Koeffizienten apund anb in Gleichung (5.2) auf
dem Grobgitter werden direkt aus den restringierten Koordinaten und den Kontrollfunk-
tionen berechnet. Beim FAS–Verfahren sind die in den Koeffizienten ap,anb bzw. Sφ
vorkommenden Metrikkoeffizienten von der L¨
osung abh¨
angig und werden deswegen ak-
tualisiert, sobald die neuen Koordinaten berechnet werden. Beim CS-Verfahren werden sie
f¨
ur einen V–oder W–Zyklus festgehalten. Das heißt, dass die Metrikkoeffizienten immer
dann ge¨
andert werden, wenn ein neuer Zyklus beginnt.
Beim ¨
Ubergang vom Feingitter mit der Maschenweite hzum Grobgitter mit der Maschen-
weite Hsind folgende Rechenschritte durchzuf¨
uhren:
•Berechnung der Residuen Rhauf dem Feingitter h: Durch ν1-malige Vorgl¨
attung
wird eine N¨
aherungsl¨
osung ˜
φermittelt, die die Approximationsgleichung bis auf das
Residuum Rh
φerf¨
ullt,
Rh
φ=−ap˜
φh
p+X
nb
anb ˜
φh
nb +Sh
φ−fh
φ.(5.4)
Hierin sind die Koeffizienten ah
p,ah
nb und Sh
φaus den N¨
aherungsl¨
osungen xhund den
Kontrollfunktionen Phberechnet.
•Restriktion der Koordinaten und der Kontrollfunktionen durch
xH= [IH
h]˜xh, PH= [IH
h]Ph(5.5)
und der Residuen Rhdurch
RH
φ= [ˆ
IH
h]∗Rh
φ.(5.6)
50 Kapitel 5. Gittergenerierung und Gl¨
attung mit Mehrgitterverfahren
mit [ˆ
IH
h]∗= (H/h)2[ˆ
IH
h] = 4[ˆ
IH
h].
Die Koordinaten ˜xhund die Kontrollfunktionen Phwerden auf das Grobgitter
Hinterpoliert. Die Modifikation des Restriktionsoperators ist deshalb notwendig,
dass bei der finiten Approximation die urspr¨
unglichen Gleichungen mit h2und der
Jacobi–Determinante multipliziert werden. W¨
ahrend die Jacobi–Determinante auf
Grob– und Feingitter etwa gleich bleibt, vergr¨
oßert sich hbeim ¨
Ubergang vom Fein–
zum Grobgitter immer um den Faktor 2.
xHund PHwerden beim CS-Verfahren nur zur Berechnung der Metrikkoeffizi-
enten und der Koeffizienten anb verwendet. Beim FAS-Verfahren wird xHals die
Startl¨
osung der Koordinaten auf dem Grobgitter benutzt. Beim CS-Verfahren wird
die Startl¨
osung zu Null gesetzt.
•Berechnung der rechten Seite der Grobgittergleichungen fH: Vor der Gl¨
attung auf
dem Grobgitter wird zuerst die rechte Seite fHberechnet. Beim FAS-Verfahren ist
fHwie folgt bestimmt:
fH
φ=−apφH
p+X
nb
anbφH
nb +SH
φ−RH
φ.(5.7)
Hierin werden die Koeffizienten ap,anb und SH
φaus der Startl¨
osung xHund PH
berechnet. Beim CS-Verfahren wird sie durch
fH
φ=−RH
φ.(5.8)
ersetzt. F¨
ur die Berechnung der Grobgitterkorrektur bleibt sie konstant.
Der ¨
Ubergang vom Grobgitter Hzum Feingitter hist wie folgt durchzuf¨
uhren:
•Berechnung der Grobgitterkorrektur δφH: Nach ν2-maliger Iteration auf dem Grob-
gitter Hwird eine N¨
aherungsl¨
osung ˜
φHermittelt. Die Korrektur wird beim FAS–
Verfahren durch
δφH=˜
φH−φH(5.9)
und beim CS-Verfahren durch
δφH=˜
φH(5.10)
berechnet.
5.3. Anwendung des Gittergl¨
attungsverfahrens 51
•Prolongation der Korrektur und Korrektur der Feingitterl¨
osung: Nach der Berech-
nung der Grobgitterkorrektur wird sie auf das Feingitter interpoliert, um die N¨
ahe-
rungsl¨
osung des Feingitters zu korrigieren:
δφh= [Ih
H]δφH, φh=˜
φh+λδφh.(5.11)
Um den Mangel an Gl¨
attung des Grobgitterfehlers auszugleichen wird hier zur Sta-
bilisierung ein Unterrelaxationsfaktor λ(0 < λ ≤1) f¨
ur den Fehler verwendet. Die
verbesserte L¨
osung wird als neue Startl¨
osung zur Nachgl¨
attung auf dem Feingitter
benutzt.
5.3 Anwendung des Gittergl¨
attungsverfahrens
In diesem Abschnitt werden die Anwendungen des oben beschriebenen Verfahren zur Git-
tergenerierung mit dem Mehrgitterverfahren und der effektiven Behandlung der Randbe-
dingungen demonstriert. Die mittels transfiniter Interpolation generierten Gitter dienen
als Startl¨
osung. Die zweidimensionalen Testf¨
alle wurden auf einem Rechner vom Typ
SiliconGraphics Indy mit einem 100 MHZ IP22 Prozessor durchgef¨
uhrt. F¨
ur die dreidi-
mensionalen Testf¨
alle wurde ein Rechner vom Typ SiliconGraphics mit einem 175 MHZ
IP28 Prozessor benutzt. Die CPU-Zeit wird in [sec] berechnet.
5.3.1 2D ω-Typ Geometrie
Die Geometrie des ersten zweidimensionalen Testfalls ist in Abbildung 5.1 dargestellt.
Wegen der Spitzkehren am unteren Rand ist es schwierig, mit einem Standardverfahren
ein g¨
unstiges Gitter zu erzeugen. Wenn man die Geometrie in angemessene Bl¨
ocke auf-
teilt, kann man sicherlich ein gutes Netz erzeugen. In den meisten F¨
allen m¨
ochte man
jedoch Gitter mit einfacher Topologie erzeugen, d.h. mit einem Block die ganze Geome-
trie erfassen. Abbildung 5.1a zeigt ein algebraisch generiertes Netz. Es ist zu erkennen,
dass die Gitterlinien sich ¨
uberschneiden. Dieses Netz wird als Startl¨
osung verwendet.
Beim Mehrgitterverfahren werden drei Gitterebenen mit 73 ×41, 37 ×21 und 19 ×11
Gitterpunkten verwendet. Die Vorgl¨
attungs- und Nachgl¨
attungszahl werden auf 2 bzw. 1
fixiert. Abbildung 5.1 (b) und (c) zeigt das gegl¨
attete Gitter und ein Detail in der N¨
ahe
der Spitzkehren.
In den Tabellen 5.1 - 5.2 werden f¨
ur die Gl¨
attung mit den Unterrelaxationsfaktoren ωf¨
ur
das SOR-Verfahren und λf¨
ur die Korrektur sowie dem Zyklustyp die Rechenzeit beim
Eingitter (EG)- und Mehrgitterverfahren gegen¨
ubergestellt. Mit Hilfe des Mehrgitterver-
fahrens l¨
asst sich eine ¨
uber 10-fache Beschleunigung gegen¨
uber dem Eingitterverfahren
erreichen. Die Anwendung des Unterrelaxationsfaktors λist wichtig, um die Stabilit¨
at zu
verbessern und eine konvergente L¨
osung zu erhalten.
52 Kapitel 5. Gittergenerierung und Gl¨
attung mit Mehrgitterverfahren
GRIDPLLX
GRID BLOCKS 1
4/03/97 13:46:32
X
Y
-1.51 -0.71 0.09 0.89 1.69
-0.45
0.35
1.15
1.95
GRIDPLLX
GRID BLOCKS 1
3/14/97 16:32:31 X
Y
-1.52 -0.72 0.08 0.88 1.68
-0.34
0.46
1.26
2.06
GRIDPLLX
GRID BLOCKS 1
3/14/97 16:32:31 X
Y
-0.72 -0.52 -0.32 -0.12 0.08
-0.14
0.06
0.26
0.46
b)
c)
a)
Abbildung 5.1: Gitterdarstellung f¨
ur 2D ω-Typ Geometrie: 1. a) Startl¨
osung; b) Globale
Darstellung nach der Gl¨
attung; c) Detail in der N¨
ahe der Spitzkehren
Gitterebenen 73 ×41, 37 ×21, 19 ×11
Relaxationsfaktoren ω= 0.6, λ = 0.7ω= 0.8, λ = 0.6
Zyklus-Typ EG V W EG V W
CPU-Zeit 32.435 2.315 2.633 22.803 1.944 2.203
Beschleunigungsfaktoren - 14.011 12.319 - 11.730 10.351
Tabelle 5.1: Mehrgitterparameter und Beschleunigungsfaktoren des CS-Verfahrens f¨
ur 2D
ω-Typ Geometrie
5.3. Anwendung des Gittergl¨
attungsverfahrens 53
Gitterebenen 73 ×41, 37 ×21, 19 ×11
Relaxationsfaktoren ω= 0.6, λ = 0.7ω= 0.8, λ = 0.6
Zyklus-Typ EG V W EG V W
CPU-Zeit 32.435 2.688 2.616 22.803 2.211 2.061
Beschleunigungsfaktoren - 12.067 12.399 - 10.313 11.064
Tabelle 5.2: Mehrgitterparameter und Beschleunigungsfaktoren des FAS-Verfahrens f¨
ur
2D ω-Typ Geometrie
Tabelle 5.1 zeigt den Leistungsvergleich des CS-Verfahrens. Wie aus dieser Tabelle her-
vorgeht, ist die Rechenzeit bei der Nutzung des Eingitterverfahrens und des Mehrgitter-
verfahrens mit ω= 0.8 deutlich geringer als mit ω= 0.6. Mit ω= 0.6 und λ= 0.7
erh¨
alt man jedoch h¨
ohere Beschleunigungsfaktoren, wobei ein Beschleunigungsfaktor von
14 beim V-Zyklus gegen¨
uber 12 beim W-Zyklus vorliegt. Mit ω= 0.8 und λ= 0.6 l¨
asst
sich ein Beschleunigungsfaktor von 11.7 beim V-Zyklus gegen¨
uber 10.4 beim W-Zyklus
erreichen. Beim CS-Verfahren bringt ein W-Zyklus keine Konvergenzbeschleunigung ge-
gen¨
uber dem V-Zyklus.
Tabelle 5.2 stellt das Leistungsverhalten des FAS-Verfahrens dar. Beim FAS-Verfahren
zeigt der W-Zyklus ein besseres Konvergenzsverhalten mit Beschleunigungsfaktoren von
12.4 bei ω= 0.6, λ= 0.7 und von 11.0 bei ω= 0.8, λ= 0.6 gegen¨
uber dem V-Zyklus mit
Beschleunigungsfaktoren von 12.0 bei ω= 0.6, λ= 0.7 und von 10.3 bei ω= 0.8, λ= 0.6.
Es ist interessant, dass in diesem Testfall das CS-Verfahren mit V-Zyklus effektiver ar-
beitet als mit W-Zyklus, w¨
ahrend das FAS-Verfahren besser mit W-Zyklus zu benutzen ist.
Abbildung 5.2 zeigt den Residuenverlauf ¨
uber der Zyklenanzahl f¨
ur das CS- bzw. FAS-
Verfahren. Wie hier zu sehen ist, bringt das Mehrgitterverfahren nicht nur Konvergenz-
beschleunigung sondern auch gr¨
oßere Stabilit¨
at gegen¨
uber dem Eingitterverfahren. Bis zu
100 Zyklen schwankt das lokale maximale Residuum sehr stark und sinkt danach sehr
langsam beim Eingitterverfahren, w¨
ahrend es beim Mehrgitterverfahren mit linearer Ge-
schwindigkeit abnimmt.
54 Kapitel 5. Gittergenerierung und Gl¨
attung mit Mehrgitterverfahren
0 20 40 60 80 100
Zyklen
10−5
10−4
10−3
10−2
max. Residuum
Eingitter
CS−V
CS−W
0 20 40 60 80 100
Zyklen
10−5
10−4
10−3
10−2
max. Residuum
Eingitter
FAS−V
FAS−W
Abbildung 5.2: Residuenverl¨
aufe ¨
uber der Zyklenanzahl f¨
ur 2D ω-Typ Geometrie
5.3. Anwendung des Gittergl¨
attungsverfahrens 55
5.3.2 2D Fl¨
ugelkonfiguration
Abbildung 5.3a zeigt ein Gitter f¨
ur eine zweidimensionale Fl¨
ugelkonfiguration, die eine
typische Geometrie in der Aerodynamik repr¨
asentiert. Dieses Gitter kann erfolgreich durch
verschiedene Verfahren wie z.B. algebraische oder elliptische Verfahren erzeugt werden.
GRIDPLLX
GRID BLOCKS 1
4/03/97 12:08:49
X
Y
-0.57 -0.07 0.43 0.93 1.43
-0.88
-0.38
0.12
0.62
GRIDPLLX
GRID BLOCKS 1
4/03/97 12:01:29 X
Y
-0.116 -0.053 0.009 0.072 0.134
-0.039
0.024
0.086
0.149
GRIDPLLX
GRID BLOCKS 1
4/03/97 12:15:50 X
Y
0.763 0.888 1.013 1.138 1.263
-0.353
-0.228
-0.103
0.022
a)
b)
c)
Abbildung 5.3: Gitterdarstellung um eine 2D Fl¨
ugelkonfiguration: a) Globale Darstellung
nach der Gl¨
attung; b), c) Detail in der N¨
ahe der Nase und der Hinterkante
Die CS- und FAS-Verfahren werden jeweils mit V- und W-Zyklus untersucht. Beim Mehr-
gitterverfahren werden vier Gitterebenen mit 321 ×65, 161 ×33, 81 ×17 und 41 ×9
Gitterpunkten verwendet. Die globale Darstellung des generierten Netzes und Details in
der N¨
ahe der Nase und der Hinterkante werden in Abbildung 5.3 dargestellt. Das Gitter
weist gute Orthogonalit¨
at und den gew¨
unschten Wandabstand in der ersten Zelle auf.
In Tabelle 5.3 werden die Mehrgitterparameter, CPU-Zeit und Beschleunigungsfaktoren
gegen¨
ubergestellt. Beim Eingitterverfahren werden f¨
ur die Gl¨
attung ¨
uber 1000s ben¨
otigt,
56 Kapitel 5. Gittergenerierung und Gl¨
attung mit Mehrgitterverfahren
Gitterebenen 321 ×65, 161 ×33, 81 ×17, 41 ×9
Relaxationsfaktoren ω= 0.6, λ = 1.0ω= 0.8, λ = 0.8
Zyklus-Typ EG V W EG V W
CPU-Zeit 1695.611 46.927 50.079 1070.104 55.826 48.545
Beschleunigungsfaktoren - 36.133 33.858 - 19.169 22.044
Tabelle 5.3: Mehrgitterparameter und Beschleunigungsfaktoren des CS-Verfahrens f¨
ur
eine 2D Fl¨
ugelkonfiguration
w¨
ahrend beim Mehrgitterverfahren weniger als 60s erforderlich sind, bis das maximale
Residuum auf dem feinsten Gitter unter ein vorgegebenes Konvergenzkriterium von 10−5
gefallen ist. In Abh¨
angigkeit von den Relaxationsfaktoren und dem Zyklus-Typ werden
Beschleunigungsfaktoren von 19.2 bis 36.1 beim CS-Verfahren erreicht.
Wie in Tabelle 5.4 zu sehen, ist die Leistungsf¨
ahigkeit des FAS-Verfahrens etwas niedri-
ger als die des CS-Verfahrens. Es werden nur Beschleunigungsfaktoren von 18.6 bis 28.4
erreicht.
Gitterebenen 321 ×65, 161 ×33, 81 ×17, 41 ×9
Relaxationsfaktoren ω= 0.6, λ = 1.0ω= 0.8, λ = 0.9
Zyklus-Typ EG V W EG V W
CPU-Zeit 1661.802 65.292 59.780 1069.481 55.747 57.561
Beschleunigungsfaktoren - 25.970 28.364 - 19.196 18.591
Tabelle 5.4: Mehrgitterparameter und Beschleunigungsfaktoren des FAS-Verfahrens f¨
ur
eine 2D Fl¨
ugelkonfiguration
Abbildung 5.4 zeigt die Residuenverl¨
aufe ¨
uber der Zyklenanzahl f¨
ur CS- bzw. FAS-Verfahren
mit ω= 0.8. Es ist zu erkennen, dass das Eingitterverfahren stabil ist, aber das maximale
Residuum sehr langsam sinkt. Beim Mehrgitterverfahren nimmt hingegen das maximale
Residuum mit einer fast linearen Geschwindigkeit ab.
5.3. Anwendung des Gittergl¨
attungsverfahrens 57
0 40 80 120 160
Zyklen
10−4
10−3
10−2
max. Residuum
Eingitter
CS−V
CS−W
0 40 80 120 160
Zyklen
10−4
10−3
10−2
max. Residuum
Eingitter
FAS−V
FAS−W
Abbildung 5.4: Residuenverl¨
aufe ¨
uber der Zyklenanzahl f¨
ur eine 2D Fl¨
ugelkonfiguration
5.3.3 3D Fl¨
ugelkonfiguration
Der numerischen Simulation in der industriellen Konfigurationsaerodynamik sind heute
noch enge Grenzen gesetzt. Eines der gr¨
oßten Probleme stellt die Erzeugung geeigneter
Rechennetze f¨
ur komplexe Konfigurationen dar, da Ergebnisse der Str¨
omungssimulation
wesentlich von der Gitterqualit¨
at abh¨
angen. Durch hochentwickelte elliptische Gittergene-
rierungsverfahren wie z.B. das biharmonische Verfahren kann die Gitterqualit¨
at deutlich
verbessert werden.
58 Kapitel 5. Gittergenerierung und Gl¨
attung mit Mehrgitterverfahren
Das im Rahmen dieser Arbeit vorgestellte biharmonische Verfahren zur Gittergenerierung
verspricht deshalb die Anwendungsm¨
oglichkeit von Simulation zu erweitern. Besondere
Herausforderung bei der Generierung solcher Gitter sind die Erzeugungen wandorthogo-
naler Gitterlinien und die Gl¨
attung des inneren Gebiets.
Zur Untersuchung der Leistungsf¨
ahigkeit des Mehrgitterverfahrens f¨
ur dreidimensionale
Fl¨
ugelkonfigurationen wird ein Testfall mit einem Netz f¨
ur die Simulation mit Euler-
Gleichungen verwendet. Beim Mehrgitterverfahren werden drei Gitterebenen jeweils mit
121 ×25 ×65, 61 ×13 ×33 und 31 ×7×17 Gitterpunkten gebildet. Das gegl¨
attete Netz
ist in Abbildung 5.5 dargestellt. Orthogonalit¨
at in Wandn¨
ahe und der Wandabstand sind
wie gew¨
unscht.
In Tabelle 5.5 werden wiederum die Mehrgitterparameter, CPU-Zeit und Beschleunigungs-
faktoren gegen¨
ubergestellt. Bei der Gl¨
attung mit dem Mehrgitterverfahren wird λ= 1.0
gesetzt. Bei diesem Testfall ist das FAS effizienter als das CS-Verfahren. Es wird ein Be-
schleunigungsfaktor von 13.0 beim FAS-Verfahren mit dem W-Zyklus erreicht, gegen¨
uber
11.1 beim CS-Verfahren mit dem V-Zyklus. Die Residuenverl¨
aufe in Abbildung 5.6 zeigen,
dass das FAS-Verfahren stabiler als das CS-Verfahren ist und das maximale Residuum
beim FAS-Verfahren schneller als beim CS-Verfahren sinkt.
Multi-Block Plot
3D XYZV MB= 1 I= 34 plane
9 Feb 97 17:14:59
0.321 0.352 0.383 0.415 0.446
-0.129
-0.098
-0.067
-0.035
Abbildung 5.5: Gitterdarstellung einer 3D Fl¨
ugelkonfiguration
5.3. Anwendung des Gittergl¨
attungsverfahrens 59
Gitterebenen 121 ×25 ×65, 61 ×13 ×33, 31 ×7×17
EG CS FAS
Zyklus-Typ V W V W
CPU-Zeit 15781 1427 1500 1222 1214
Beschleunigungsfaktoren - 11.10 10.52 12.91 13.00
Tabelle 5.5: Mehrgitterparameter und Beschleunigungsfaktoren des Mehrgitterverfahrens
mit ω= 0.8 und λ= 1.0 f¨
ur eine 3D Fl¨
ugelkonfiguration
0 40 80 120 160
Zyklen
10−8
10−7
10−6
max. Residuum
Eingitter
CS − V
CS − W
0 40 80 120 160
Zyklen
10−8
10−7
10−6
max Residuum
Eingitter
FAS − V
FAS − W
Abbildung 5.6: Residuenverl¨
aufe ¨
uber der Zyklenanzahl f¨
ur eine 3D Fl¨
ugelkonfiguration
60 Kapitel 5. Gittergenerierung und Gl¨
attung mit Mehrgitterverfahren
5.3.4 Fl¨
ugel-Rumpf Konfiguration
Abbildung 5.7 zeigt einen Block mit 161 ×17 ×65 Gitterpunkten von einem gegl¨
attete
Netz um eine Fl¨
ugel-Rumpf-Konfiguration. Dieses Gitter ist f¨
ur Navier-Stokes Rechnung
geeignet. Das generierte Netz ist im ganzen Gebiet glatt und an allen W¨
ande orthogonal.
Tabelle 5.6 zeigt die Leistungsf¨
ahigkeit des Mehrgitterverfahrens mit drei Gitterebenen
(161 ×17 ×65, 81 ×9×33 und 41 ×5×17 Gitterpunkte). In diesem Testfall erreicht das
Mehrgitterverfahren nur etwa f¨
unf fache Beschleunigung gegen¨
uber dem Eingitterverfah-
ren. Die Leistungsf¨
ahigkeit ist relativ gering im Vergleich zu zweidimensionalen Testf¨
allen.
Der Grund daf¨
ur ist die geringe Gitterausdehnung in j-Richtung von 17 Punkten im fein-
sten und nur 5 Punkten im gr¨
obsten Gitter, durch die keine vollst¨
andige Erfassung der
Komplexit¨
at m¨
oglich ist. Abbildung 5.8 zeigt die Residuenverl¨
aufe ¨
uber der Zyklenanzahl.
Die Anwendung des W-Zyklus ben¨
otigt nur wenige Iterationen auf dem feinsten Gitter,
daf¨
ur jedoch mehr CPU-Zeit.
Plotted by jianping@effendi on Mon Mar 2 11:37:39 MET 1998, using the ICEM CFD LEO Grid Visualizer 3.3.2.
X
Y
Z
Abbildung 5.7: Gitterdarstellung einer Fl¨
ugel-Rumpf-Konfiguration
5.3. Anwendung des Gittergl¨
attungsverfahrens 61
Gitterebenen 161 ×17 ×65, 161 ×17 ×65, 161 ×17 ×65
EG CS FAS
Zyklus-Typ V W V W
CPU-Zeit 16122 3053 4010 5234 5374
Beschleunigungsfaktoren - 5.28 4.02 3.08 3.20
Tabelle 5.6: Mehrgitterparameter und Beschleunigungsfaktoren des Mehrgitterverfahrens
mit ω= 0.6 und λ= 0.6 f¨
ur eine Fl¨
ugel-Rumpf Konfiguration
0 200 400 600 800
Zyklen
10−2
10−1
100
101
102
max. Residuum
Eingitter
CS−V
CS−W
0 200 400 600 800
Zyklen
10−2
10−1
100
101
102
max. Residuum
single grid
3 level, FAS−V
3 level, FAS−W
Abbildung 5.8: Residuenverl¨
aufe ¨
uber der Zyklenanzahl f¨
ur eine Fl¨
ugel-Rumpf Konfigu-
ration
62 Kapitel 5. Gittergenerierung und Gl¨
attung mit Mehrgitterverfahren
5.3.5 Fl¨
ugel-Rumpf-Leitwerk Konfiguration
In einem letzten Testfall wird eine generische Fl¨
ugel-Rumpf-Leitwerk Konfiguration (H8Y)
behandelt. Abbildung 5.9 zeigt die Gitterdarstellung der Fl¨
ugel-Rumpf-Leitwerk Konfi-
guration mit 289 ×89 ×89 Gitterpunkten nach der Gl¨
attung. Zur Gl¨
attung wird nur
das CS-Verfahren mit vier Gitterebenen verwendet. Hier gibt es keinen Vergleich der Lei-
stungsf¨
ahigkeit des Mehrgitterverfahrens gegen¨
uber dem Eingitterverfahren.
Abbildung 5.9: Gitterdarstellung einer generischen Fl¨
ugel-Rumpf-Leitwerk Konfiguration
Abbildung 5.10 zeigt einen k¨
orpernahen Schnitt des Gitters vor der Gl¨
attung, sowie nach
der Gl¨
attung. Wie in Abbildung 5.10(a) zu erkennen ist, sind die Gitterlinien vor der Git-
tergl¨
attung stark verschert, dies kann zu erh¨
ohter numerischer Diffusion in die Simulation
f¨
uhren. Abbildung 5.10(b) zeigt den selben Schnitt des Gitters nach der Gl¨
attung. Die
Verscherung der Gitterlinien hat deutlich abgenommen.
5.3. Anwendung des Gittergl¨
attungsverfahrens 63
X
Y
Z
X
Y
Z
(a) (b)
Abbildung 5.10: Gitterdetail einer Fl¨
ugel-Rumpf-Leitwerk Konfiguration: (a) vor der
Gl¨
attung; (b) nach der Gl¨
attung
Die oben dargestellten Ergebnisse zeigen, dass das jetzige Verfahren eine effektive Ge-
staltung des Gitters im wandnahen Bereich erm¨
oglicht. Sowohl der Schnittwinkel zu den
R¨
andern als auch die Maschenweite an den R¨
andern k¨
onnen festgelegt werden. Weil die
Kontrollfunktionen im Gebietsinneren durch die elliptische differentielle Beziehung be-
stimmt werden, haben sie im gesamten Gitterbereich einen glatten Verlauf. Das Mehrgit-
terverfahren kann die Konvergenzgeschwindigkeit deutlich beschleunigen.
Das Mehrgitterverfahren wird nur f¨
ur die Berechnung der Koordinaten eingesetzt. Eine
gleichzeitige Anwendung f¨
ur die Kontrollfunktionen kann der Erfahrung nach die Konver-
genzbeschleunigung nicht verbessern und beim FAS-Verfahren wird manchmal ¨
uberhaupt
keine Konvergenz erzielt.
Bemerkung
F¨
ur die praktische Gittergl¨
attung gibt es kein allgemeing¨
ultiges Kriterium zur Kover-
genzbeurteilung. Das Gitter wird meist einige Gl¨
attungszyklen unterworfen und dann
64 Kapitel 5. Gittergenerierung und Gl¨
attung mit Mehrgitterverfahren
manuell bez¨
uglich seiner Qualit¨
at bewertet. Bei unzul¨
anglicher G¨
ute wird die Gl¨
attung
entsprechend fortgef¨
uhrt. Dieses Vorgehen macht die Vorgabe eines Konvergenzkriteriums
unn¨
otig. Wegen des akademischen Ansatzes dieser Arbeit wurde hier aber ein Konvergenz-
kriterium vorgegeben.
Kapitel 6
Mehrgitterverfahren zur L¨
osung der
Navier-Stokes Gleichungen
Zur Beschleunigung der numerischen L¨
osung der Navier-Stokes Gleichungen eignen sich
besonders Mehrgitterverfahren, wobei meist FMG-FAS Anwendung findet. Bei bisherigem
Umsetzung des Verfahrens ergibt sich das Problem, dass die Massenfl¨
usse, restringierten
Geschwindigkeiten und die Turbulenzvariablen auf unterschiedlichen Gitterebenen nicht
zusammenpassen. Um dieses Problem zu ¨
uberwinden, wird hier ein modifizierter Mehrgit-
teralgorithmus entwickelt, der auf einem FMG-Algorithmus mit V-Zyklus ohne Restrikti-
onsprozedur f¨
ur die Gittervariablen mit Ausnahme der Residuen basiert. Diese Eigenschaft
f¨
uhrt sowohl zur Vereinfachung der Mehrgitterstrategie als auch der Programmstruktur.
Im folgenden wird zun¨
achst eine Einf¨
uhrung in die Komponenten von Mehrgitterverfah-
ren zur L¨
osung der Navier-Stokes Gleichungen gegeben, und dann werden der FMG-FAS-
sowie der modifizierter FMG-Algorithmus vorgestellt. Anschließend werden die Grobgit-
tergleichungen f¨
ur die Navier-Stokes Gleichungen und die Behandlung der Druckkorrek-
tur mit dem modifizierten Algorithmus im Detail beschrieben. Die Unterschiede dieses
modifizierten Algorithmus zum Standardverfahren liegen sowohl in der Behandlung der
Druckkorrektur auf dem Grobgitter als auch in der Behandlung anderer Variablen, wie
z.B. der Geschwindigkeitskomponenten.
6.1 Komponenten von Mehrgitterfahren
Gittervergr¨
oberung: Zur Generierung von Gittern verschiedener Ebenen wird folgen-
de Vorgehensweise praktiziert: Ausgehend von dem der Geometrie und den Str¨
omungs-
verh¨
altnissen angepassten feinsten Gitter werden rekursiv gr¨
obere Gitter gebildet, wo-
bei nur jede zweite Gitterlinie in jeder Richtung des Feingitters zur Bildung des n¨
achst
gr¨
oberen Gitters verwendet wird. Ein Grobgitterbilanzvolumen umfasst damit jeweils vier
Feingitterbilanzvolumen in zweidimensionalen F¨
allen und acht Feingitterbilanzvolumen in
dreidimensionalen F¨
allen. Eine Konsequenz davon ist, dass die Residuen auf dem feinen
Gitter in einfacher Weise auf das Grobgitter restringiert werden k¨
onnen.
65
66 Kapitel 6. Mehrgitterverfahren zur L¨
osung der Navier-Stokes Gleichungen
Restriktion: Beim FV-Verfahren wird die Restriktion wie folgt berechnet: Das Resi-
duum einer Grobgitterzelle wird durch Summierung der Residuen der Feingitterzellen,
die dem fusionierten Grobgittervolumen entsprechen, bestimmt. Werden Gr¨
oßen, der zu
einer Grobgitterzelle geh¨
orenden Feingitterzellen mit hgekennzeichnet, dann lautet die
Restriktion des Residuums:
RH=XRh.(6.1)
Die f¨
ur das FAS notwendigen Variablenwerte auf dem Grobgitter werden durch die vo-
lumengewichtete Mittelung aus den Werten von der zu einem Grobgitter geh¨
orenden
Feingitterzellen berechnet:
φH= [IH
h]φh=PδV hφh
PδV h(6.2)
Beim Standardverfahren ist nicht nur die Startl¨
osung auf dem Grobgitter erforderlich,
sondern es ist ebenfalls notwendig, die Massenfl¨
usse an den Kontrollvolumenfl¨
achen zu be-
stimmen. Dabei ist es besonders wichtig, den konservativen Charakter des FV-Verfahrens
auch auf dem Grobgitter zu erhalten. Dies erreicht man dadurch, dass der Massenfluß
an einer Kontrollvolumenfl¨
ache auf dem Grobgitter als Summe der Massenfl¨
usse an den
entsprechenden Kontrollvolumenfl¨
achen auf dem n¨
achste feineren Gitter gebildet wird:
˙mH=X˙mh.(6.3)
Prolongation: Die Prolongation dient der ¨
Ubertragung der Grobgitterkorrekturen auf
das Feingitter. Beim FV-Verfahren auf strukturierten Gittern werden ¨
ublicherweise die
Korrekturen f¨
ur Feingitterst¨
utzstellen aus den Werten der umliegenden Grobgitterst¨
utz-
stellen linear interpoliert. Die Prolongation der Grobgitterkorrektur auf das feine Gitter
erfolgt hier bei zweidimensionalen Gittern mit einer bilinearen Interpolation und bei drei-
dimensionalen Gittern mit einem auf Gradienten der Variablen basierten Verfahren. Diese
Vorgehensweise wurde von Xue (1998) beschrieben.
6.2 Modifizierte Full-Multigrid Methode
Die FMG-Methode wird oft zur Berechnung der Navier-Stokes Gleichung benutzt. Ab-
bildung 4.2 zeigt ein Flussdiagramm mit drei Gitterebenen. Dabei ist die Startl¨
osung
des Grobgitters vom n¨
achst feineren Gitter restringiert. Alternativ dazu k¨
onnen auch
Startwerte φHvom vorherigen Zyklus ¨
ubernommen werden. Nur die Residuen auf dem
feinsten Gitter werden auf das Grobgitter ¨
ubertragen. Dies f¨
uhrt zu einer modifizierten
FMG-Methode mit V-Zyklus, wie in Abbildung 6.1 dargestellt. ν1und ν2bezeichnen die
Vorgl¨
attungs- bzw. Nachgl¨
attungszahl und ν0die Zahl der Relaxationsiteration auf dem
gr¨
obsten Gitter.
Im modifizierten FMG-Algorithmus wird die Grobgittergleichung ganz genau wie Glei-
chung (4.9) konstruiert. Sobald die N¨
aherungsl¨
osung berechnet ist, wird die Grobgitter-
korrektur durch Gleichung (4.10) bestimmt. Wenn die Residuen auf dem feinen Gitter
6.3. Grobgittergleichungen 67
.
Prolongation
Konvergierte Lösung
Vor−/Nachglättung
Restriktion
Übernahme der Variablen
vom vorherigen Zyklus
Mod. V−Zyklus
u,v,w,p ...
Modifiziertes FMG
Glättung auf dem gröbsten Gitter
Abbildung 6.1: Schematische Darstellung des modifizierten Full Multigrid-Verfahrens
verschwinden, ist Gleichung (4.9) identisch erf¨
ullt und es wird keine Korrektur zur Fein-
gitterl¨
osung produziert, da ˜
φHdurch φHinitialisiert wird. Zuerst wird die L¨
osung auf dem
gr¨
obsten Gitter erreicht. Da die L¨
osung auf dem Grobgitter nur wenige Kontrollvolumen
ben¨
otigt, sind sowohl die Zahl der erforderlichen Iterationen als auch die Rechenzeit klein.
Die konvergierte L¨
osung wird dann zum n¨
achst feineren Gitter extrapoliert und dient als
Startl¨
osung f¨
ur den Zweigitteralgorithmus. Nachdem die L¨
osung auf diesem Gitter erreicht
ist, wird sie entsprechend auf das folgende feinere Gitter ¨
ubertragen. Dies wird solange
fortgesetzt, bis das feinste Gitter erreicht ist.
In Abbildung 6.1 symbolisiert der Pfeil &die Residuen¨
ubertragung vom feinen Gitter
auf das folgende gr¨
obere Gitter. Der Pfeil %entspricht den Prozessen beim ¨
Ubergang
vom groben Gitter zum folgenden feineren Gitter, die die Korrekturberechnung auf dem
groben Gitter, die Prolongation der Grobgitterkorrektur und die Korrektur der vorherge-
henden N¨
aherungsl¨
osung auf dem feinen Gitter umfassen. Durch den gestrichelten Pfeil
werden die Prozesse gekennzeichnet, bei denen Variablen vom vorhergehenden Zyklus di-
rekt ¨
ubernommen werden.
Der einzige Unterschied zwischen dem Standard-FMG und dem modifizierten FMG-
Algorithmus besteht in der verwendeten Startl¨
osung. Im modifizierten FMG-Algorithmus
wird kein Restriktionsoperator IH
hf¨
ur die Gittervariable ben¨
otigt. Die Restriktion des
Residuums kann einfach durch Summierung der entsprechenden Residuen auf dem feinen
Gitter bestimmt werden.
6.3 Grobgittergleichungen
Bei der L¨
osung der Navier-Stokes Gleichungen ist zu beachten, dass der konvektive Term
in der Impulsgleichung nichtlinear von der Geschwindigkeit abh¨
angig ist, der Druckgra-
dient hingegen linear vom Druck. Zur ¨
ubersichtlicheren Darstellung der Diskretisierung
68 Kapitel 6. Mehrgitterverfahren zur L¨
osung der Navier-Stokes Gleichungen
werden die Impulsgleichung (3.18) f¨
ur die Geschwindigkeitskomponente uiwie folgt um-
geschrieben:
L{ui}+Di(p) + S{ui}=f{ui},(6.4)
hierin repr¨
asentiert L{ui}den impliziten Anteil der Impulsbilanz, Di(p) die Abh¨
angigkeit
vom Druck und S{ui}den expliziten Anteil in Abh¨
angigkeit von der Geschwindigkeit.
Die rechte Seite f{ui}stellt den von den L¨
osungsvariablen unabh¨
angigen Anteil dar, der
auf dem feinsten Gitter Null und auf dem Grobgitter von Null verschieden wird.
Nach der Durchf¨
uhrung einiger SIMPLE-Iterationszyklen auf dem feinen Gitter wird eine
N¨
aherungsl¨
osung der Geschwindigkeit und des Druckes ermittelt, die die Impulsgleichung
bis auf ein Residuum ˜
Rh{˜uh
i}erf¨
ullt,
˜
Lh{˜uh
i}+˜
Dh
i(˜ph) + ˜
Sh{˜uh
i}=fh{ui}+˜
Rh{˜uh
i}.(6.5)
Hierbei bezeichnet ˜uh
idie N¨
aherungsl¨
osung der Geschwindigkeit und ˜phdie des Druckes.
W¨
ahrend der Iteration werden Korrekturen der Geschwindigkeiten und des Druckes ein-
gef¨
uhrt,
uh
i= ˜uh
i+δuh
i, ph= ˜ph+δph,(6.6)
um das Residuum ˜
Rh{ui}verschwinden zu lassen,
Lh{uh
i}+Dh
i(ph) + Sh{uh
i}=fh{ui}.(6.7)
Subtraktion der Gleichung (6.5) von Gleichung (6.7) liefert eine Beziehung, die als Basis
f¨
ur die Mehrgitterkopplung dient:
Lh{uh
i}+Dh
i(ph) + Sh{uh
i}=˜
Lh{˜uh
i}+˜
Sh{˜uh
i}+˜
Dh
i{˜ph}− ˜
Rh{˜uh
i}
| {z }
konstant
.(6.8)
F¨
ur eine skalare Variable φ, wie z.B. die turbulente kinetische Energie koder die turbulente
Dissipation ²erh¨
alt man folgende ¨
ahnliche Gleichung:
Lh{φh}+S{φh}=˜
Lh{˜
φh}+˜
S{˜
φh}− ˜
Rh{˜
φh}
| {z }
konstant
.(6.9)
Durch Verwendung dieser Kopplungsgleichungen k¨
onnen die Korrekturen auf dem Grob-
gitter berechnet werden. Daf¨
ur k¨
onnen FAS- oder CS-Verfahren verwendet werden. Der
FMG-FAS Verfahren und der modifizierte Mehrgitteralgorithmus werden hier verwendet.
Auf dem Grobgitter mit dem Index Hlautet die Gleichung der Geschwindigkeit ui:
˜
LH{˜uH
i}+˜
DH
i(˜pH) + ˜
SH{˜uH
i}=fH{ui}
=LH{uH
i}+DH
i{pH}+SH{uH
i}−[IH
h]˜
Rh{˜uh
i}
| {z }
konstant
,(6.10)
6.3. Grobgittergleichungen 69
und die f¨
ur skalare Variablen φ:
˜
LH{˜
φH}+˜
SH{˜
φH}=fH{φ}=LH{φH}+SH{φH}−[IH
h]˜
Rh{˜
φh}
| {z }
konstant
.(6.11)
In diesen Gleichungen sind die Variablen uH
i,pHund φHdie Startwerte der Geschwin-
digkeitenskomponente, des Druckes und der Skalare, w¨
ahrend die ˜uH
i, ˜pHund ˜
φHdie
N¨
aherungsl¨
osung auf dem Grobgitter bezeichnen. Im modifizierten Algorithmus werden
uH
i,pHund φHdirekt vom letzten Zyklus ¨
ubernommen. Auf der rechten Seite stehen
die Terme [IH
h]˜
Rh{˜uh
i}und [IH
h]˜
Rh{˜
φh}die restringierten Residuen, die auf dem feinen
Gitter berechnet werden. Die andere Terme werden basierend auf der Startl¨
osung f¨
ur das
Grobgitter berechnet. Die gesamte rechte Seite ist konstant auf dem Grobgitter H.
Zur Erf¨
ullung der Kontinuit¨
atsgleichung wird neben den Impulsgleichungen auch eine
Druckkorrekturgleichung auf dem Grobgitter gel¨
ost. Die Behandlung der Druckkorrek-
turgleichung ben¨
otigt besondere Aufmerksamkeit innerhalb des Mehrgitteralgorithmus.
Die Gleichung f¨
ur die Druckkorrektur im SIMPLE-Algorithmus kann in einer allgemeinen
Form folgendermaßen
Qm=L{p0}+Q∗
m=f{p0},(6.12)
geschrieben werden. Hierbei ist die rechte Seite f{p0}auf dem feinsten Gitter Null und auf
dem Grobgitter von Null verschieden. Q∗
mund Qmsind die Massenflussdefekte, die sich
aus den Massenfl¨
ussen vor ( ˙m∗) und nach ( ˙m) der Verbesserung durch die Druckkorrektur
ergeben:
Q∗
m= ˙m∗
e+ ˙m∗
w+ ˙m∗
n+ ˙m∗
s, Qm= ˙me+ ˙mw+ ˙mn+ ˙ms.(6.13)
Nach der Durchf¨
uhrung von einigen SIMPLE-Iterationszyklen auf dem feinen Gitter ist
die L¨
osung noch nicht konvergiert, Gleichung (6.12) f¨
uhrt zum Residuum Rh{p0}
Rh{p0}=Qh
m−fh{p0}.(6.14)
Die Gleichung f¨
ur die Druckkorrektur auf dem Grobgitter lautet wie folgt
QH
m=L{p0H}+Q∗H
m=fH{p0}=QH
m−[IH
h]Rh{p0h}
| {z }
konstant
.(6.15)
Hierbei sind QH
mdie initialen Massenflussdefekte, die vom letzten Zyklus ¨
ubernommen
werden. Weil Gleichung (6.15) mit einer von Null verschiedenen rechten Seite fH{p0}
auf dem Grobgitter gel¨
ost wird, soll der globale Massenflussdefekt durch die Summe von
fH{p0}auf dem ganzen Gebiet PfH{p0}korrigiert werden.
Mit den Initialwerten uH
i,pHund φHwerden ν1SIMPLE-Relaxationszyklen durchgef¨
uhrt.
Als Ergebnis ergibt sich die Grobgitterl¨
osung inklusive dem Grobgitterdruck. Aus dieser
k¨
onnen die Korrekturen
δuH
i= ˜uH
i−uH
i, δpH= ˜pH−pH, δφH=˜
φH−φH(6.16)
70 Kapitel 6. Mehrgitterverfahren zur L¨
osung der Navier-Stokes Gleichungen
bestimmt, auf das feine Gitter prolongiert
δuh
i= [Ih
H]δuH
i, δph= [Ih
H]δpH, δφh= [Ih
H]δφH.(6.17)
und dort zur N¨
aherungsl¨
osung addiert werden
uh
i= ˜uh
i+λuδuh
i, ph= ˜ph+λpδph, φh=˜
φh+λφδφh.(6.18)
Um den Mangel an Gl¨
attung der Grobgitterkorrektur auszugleichen wird hier zur Stabi-
lisierung ein Unterrelaxationsfaktor f¨
ur die Korrektur von λ(0 < λ ≤1) verwendet. Die
verbesserte L¨
osung wird als Startl¨
osung auf dem feinen Gitter hf¨
ur die Nachgl¨
attung
oder f¨
ur den n¨
achsten Zyklus weiter benutzt.
Die Grobgitterkorrektur des Druckes δpHwird direkt durch Gleichung (6.16) bestimmt.
Diese Vorgehensweise unterscheidet sich vom Verfahren aus der Arbeit (Lien und Lesch-
ziner 1994a, Orth 1991), wobei die sogenannte Korrektur der Druckkorrektur berechnet
wird. Zudem wird die Massenflusskorrektur aus der Geschwindigkeitskorrektur δuh
ibe-
rechnet und zum vorher bestimmten Massenfluß addiert.
6.4 Maßnahmen zur Turbulenzsimulation
Aufgrund des hochgradig nichtlinearen Charakters der Transportgleichungen ist eine ge-
sonderte Behandlung der turbulenten Scheinz¨
ahigkeit und der Quellterme bei der An-
wendung vom Mehrgitterverfahren erforderlich. In der Arbeit von Yan (1999) ist eine
Vorgehensweise, die am Beispiel von zwei- und dreidimensionalen turbulenten Str¨
omun-
gen beschrieben wird, angegeben. Sie weist folgende charakteristische Eigenschaften auf:
•Festhalten der turbulenten Scheinz¨
ahigkeit auf dem Grobgitter: Nachdem die kon-
vergierte L¨
osung bekommen wird, bleibt die dazu entsprechende turbulente Schein-
z¨
ahigkeit auf dem Grobgitter bei weiterer Mehrgitterprozedur unver¨
andert.
•Die rechte Seite fH{φ}der Grobgittergleichung (6.11) wird wie folgt behandelt:
Wenn fH{φ}>0, wird sie auf der rechten Seite ber¨
ucksichtigt, ansonsten wird sie
als ”Pseudodissipation” zur St¨
arkung der Diagonaldominanz in der Form fH{φ}
φHzum
zentralen Koeffizienten geschlagen.
•Um zu vermeiden, dass die Turbulenzgr¨
oßen kund ²nach der Prolongation negativ
werden, wird f¨
ur sie folgende Formulierung verwendet:
φh=|˜
φh+λφ[Ih
H]δφH|.
Nach den bisher gemachten Erfahrungen bringt diese Behandlung eine gute Stabilit¨
at.
Die Randbedingungen werden auf allen Gitterstufen genauso wie in jedem einzelnen Git-
teralgorithmus behandelt. Die Neumannbedingung, die f¨
ur die Druckkorrektur an allen
R¨
andern auf dem feinsten Gitter verwendet werden, sind auch auf die Druckkorrektur auf
dem gr¨
oberen Gitter anwendbar.
Kapitel 7
Validierung des
Simulationsverfahrens
Das oben beschriebene Mehrgitterverfahren sowie die modifizierte Variante f¨
ur die L¨
osung
der Navier-Stokes Gleichungen wird nun zur Simulation von Str¨
omung und W¨
arme¨
uber-
gang in komplexen Konfigurationen verwendet. Zur zuverl¨
assigen Bewertung des Al-
gorithmus wird dieser im Zusammenhang mit Konvektionsschemata h¨
oherer Ordnung,
QUICK (Quadratic Upstream Interpolation for Convection Kinematics) (Leonard 1979)
und TVD-typ MUSCL (Monotonic Upstream-centred Scheme for Conservation Laws)
(Van Leer 1977) getestet. Der Gittereinfluss wird durch die Verwendung mehrerer ex-
emplarischer Gitter mit sehr speziellen Eigenschaften gepr¨
uft. Verwendet werden karte-
sische, ungleichm¨
aßige, stark verzerrte und nicht-orthogonale Gittervarianten. Schließlich
wird das Mehrgitterverfahren in Kombination mit zwei verschiedenen Turbulenzmodellen
eingesetzt, dem k−²Modell mit Wandfunktionen sowie dem LL k−²Modell von Lien
und Leschziner (1993). Untersucht werden: eine zweidimensionale Nischenstr¨
omung, die
Str¨
omung ¨
uber eine r¨
uckw¨
artsgerichtete Stufe, bzw. einen zweidimensionalen Modellh¨
ugel
und ferner eine dreidimensionale Rohrkr¨
ummerstr¨
omung mit starkem Sekund¨
arwirbel.
Bei allen Simulationen wird die L¨
osung als konvergiert betrachtet, wenn die Summe der
absoluten Residuen von Impuls-, Kontinuit¨
ats- und Turbulenztransportgleichung bezogen
auf dem Gesamtmassenfluss auf dem feinsten Gitter unter ein Konvergenzkriterium von
10−3gefallen ist. Im Fall der Nischenstr¨
omung ergibt sich der Gesamtmassenfluss dabei
aus dem Produkt von Geschwindigkeit und Breite der bewegten Wand.
Zur Verbesserung der L¨
osungsstabilit¨
at werden alle Gleichungen relaxiert. Der Unterre-
laxationsfaktor betr¨
agt 0.7 f¨
ur die Geschwindigkeit, 0.3 f¨
ur den Druck und 0.7 f¨
ur die
Turbulenzvariablen. Innerhalb des Mehrgitteralgorithmus wird der Unterrelaxationsfak-
tor der Grobgitterkorrektur λf¨
ur alle Gittervariablen im zweidimensionalen Fall auf 1.0
und im dreidimensionalen Fall auf 0.7 gesetzt.
Alle zweidimensionalen Simulationen werden auf einer Silicon Graphics Workstation durch-
gef¨
uhrt, die dreidimensionalen auf einer Cray-T3E. Die Rechengeschwindigkeit wird an-
hand der ben¨
otigten CPU-Zeit [sec] angegeben. Als Referenz dient eine Simulation auf
71
72 Kapitel 7. Validierung des Simulationsverfahrens
Basis des Eingitterverfahrens, auf die alle weiteren Rechnungen bezogen werden.
In diesem Kapitel sollen die numerischen Eigenschaften und die Leistungsf¨
ahigkeit des
Mehrgitterverfahrens untersucht werden. Deshalb werden nicht die Str¨
omungsph¨
anomene
im Detail analysiert, sondern die Auswertung konzentriert sich vielmehr auf den Residu-
enverlauf und die erforderlichen Rechenzeiten.
7.1 2D laminare Str¨
omung in rechteckiger Nische
Der erste Testfall zur Bewertung der Mehrgittermethode ist die durch eine bewegte Wand
getriebene Nischenstr¨
omung mit quadratischem Querschnitt. Die Leistung des Standard-
verfahrens und der modifizierten Methode k¨
onnen daran verglichen werden. Die Simula-
tion wird f¨
ur verschiedene Reynoldszahlen von 100 bis 1000, gebildet mit der Geschwin-
digkeit der bewegten Wand und der Nischenbreite, sowie mit den Konvektionsschemata
MUSCL und QUICK durchgef¨
uhrt. Es kommen Gitter mit bis zu f¨
unf Gitterebenen von
16×16 bis 256×256 Gitterzellen und ungleichm¨
aßiger Gitterverteilung zum Einsatz. Die
Vorgl¨
attungszahl ν1und die Nachgl¨
attungszahl ν2sind auf 3 bzw. 2 fixiert, die Relaxati-
onszahl auf dem gr¨
obsten Gitter ν0betr¨
agt 8.
In Tabelle 7.1 ist die ben¨
otigte CPU-Zeit bei Verwendung des Eingitterverfahrens dar-
gestellt. Es zeigt sich, dass beim QUICK-Schema weniger CPU-Zeit ben¨
otigt wird als
bei MUSCL, um eine konvergierte L¨
osung zu erhalten. F¨
ur die Simulation des Eingit-
terverfahrens bei Re=100 ben¨
otigt das QUICK-Schema ungef¨
ahr 10% weniger CPU-Zeit
als das MUSCL-Schema. Bei h¨
oherer Reynoldszahl und auf feineren Gittern f¨
allt dieser
Unterschied unter 5%. In den folgenden Simulationen werden deshalb beide Schemata
verwendet.
Wie weiterhin zu erkennen ist, ben¨
otigt das Eingitterverfahren auf groben Gittern mit
Zunahme der Reynoldszahl mehr CPU-Zeit w¨
ahrend auf feineren Gittern die gegenteilige
Tendenz beobachtet wird. Dieses Verhalten ist in ¨
Ubereinstimmung mit den Ergebnissen
von Orth (1991), Lien und Leschziner (1994a) sowie Lilik u. a. (1997).
Die in den Tabellen 7.2-7.5 dargestellten Ergebnisse der Anwendung der FMG-Methode
zeigen ein erhebliches Anwachsen des Beschleunigungsfaktors mit wachsender Gitter-
punktszahl. Im Allgemeinen vermindert sich die Effektivit¨
at, wenn die Reynoldszahl erh¨
oht
wird. Dieses Ph¨
anomen ist in Einklang mit den Beobachtungen von Vanka (1986b) sowie
Lien und Leschziner (1994a). Dabei ist darauf hinzuweisen, dass der modifizierte FMG-
Algorithmus ein hohes Potential f¨
ur feine Gitter mit 128×128 und 256 ×256 Gitterzellen
besitzt. Der Beschleunigungsfaktor ist ungef¨
ahr 40 % h¨
oher als beim Standard FMG, das
in der Arbeit von Lien und Leschziner (1994a) auch eine ¨
ahnliche Beschleunigung erreicht.
Die Gr¨
unde f¨
ur die hier h¨
ohere Effektivit¨
at liegen m¨
oglicherweise in der unterschiedlichen
Programmstruktur. Zu einem ¨
ahnlichen Schluss kann man auf Grund der Ergebnisse auf
dem feinsten Gitter von 256×256 Gitterzellen und drei Gitterstufen kommen. Die Anwen-
dung des Mehrgitteralgorithmus mit vier oder f¨
unf Gitterebenen ergibt einen sehr hohen
7.2. 2D laminare Str¨
omung in verzerrter Nische 73
MUSCL QUICK
Gitter Re=100 Re=400 Re=1000 Re=100 Re=400 Re=1000
16 ×16 0.118 0.169 0.195 0.107 0.159 0.179
32 ×32 1.241 1.396 1.437 1.077 1.196 1.348
64 ×64 21.928 20.708 19.723 19.856 18.718 18.264
128 ×128 598.289 520.074 468.733 567.453 492.232 444.678
256 ×256 10994.886 8454.691 7581.088 9689.093 8204.583 7573.642
Tabelle 7.1: Ben¨
otigte CPU-Zeit des Eingitterverfahrens in Abh¨
angigkeit von Konvekti-
onsschema und Reynoldszahl auf verschiedenen Gittern f¨
ur die rechteckige
Nischenstr¨
omung
Beschleunigungsfaktor. Abbildung 7.1 zeigt Gitter und Stromlinien f¨
ur die laminare Ni-
schenstr¨
omung bei Re = 400.
mod. FMG FMG-FAS
Gitter Re 100 400 1000 100 400 1000
CPU-Zeit 0.118 0.169 0.195 0.107 0.159 0.179
16 ×16 Be. Fakt. 1.000 1.000 1.000 1.000 1.000 1.000
CPU-Zeit 0.367 0.540 0.733 0.381 0.562 0.842
32 ×32 Be. Fakt. 2.934 2.215 1.839 2.826 2.128 1.601
CPU-Zeit 1.626 1.698 1.957 1.774 2.083 2.766
64 ×64 Be. Fakt. 12.212 11.024 9.333 11.192 8.986 6.603
CPU-Zeit 9.333 7.480 7.423 11.539 9.922 10.848
128 ×128 Be. Fakt. 60.801 65.806 59.905 49.177 49.610 40.992
CPU-Zeit 42.375 30.480 26.710 53.804 43.778 44.792
256 ×256 Be. Fakt. 228.651 288.864 283.551 180.081 187.400 169.085
Tabelle 7.2: CPU-Zeit und Beschleunigungsfaktor in Abh¨
angigkeit von Mehrgitteralgo-
rithmen und Gitterebenen mit QUICK Schema f¨
ur die rechteckige Nischen-
str¨
omung, gr¨
obstes Gitter 16 ×16
In Tabelle 7.6 ist zu erkennen, dass sich die durch MUSCL und QUICK ermittelten
Geschwindigkeiten an der vertikalen Mittelachse in ¨
Ubereinstimmung mit den Ergebnissen
von Vanka (1986b) befinden, die mit einem HYBRID Schema berechnet wurden. Es wird
deutlich, dass die mit Hilfe Konvektionsschemata h¨
oherer Ordnung errechneten Daten
bereits f¨
ur weniger Gitterzellen eine gitterunabh¨
angige L¨
osung liefern.
7.2 2D laminare Str¨
omung in verzerrter Nische
Viele technische Anwendungen erfordern stark verzerrte Gitter, durch die die Effizienz der
L¨
osungsprozedur beeinflusst wird. Deshalb wird der Effekt eines schr¨
agen Gitters auf den
74 Kapitel 7. Validierung des Simulationsverfahrens
mod. FMG FMG-FAS
Gitter Re 100 400 1000 100 400 1000
CPU-Zeit 1.077 1.196 1.348 1.077 1.196 1.348
32 ×32 Be. Fakt. 1.000 1.000 1.000 1.000 1.000 1.000
CPU-Zeit 2.675 2.627 3.204 3.350 2.909 3.323
64 ×64 Be. Fakt. 7.422 7.125 5.700 5.927 6.435 5.496
CPU-Zeit 10.538 9.031 8.150 12.834 11.739 10.338
128 ×128 Be. Fakt. 53.615 54.505 54.562 44.024 41.931 43.014
CPU-Zeit 52.657 41.035 33.371 62.891 58.344 47.435
256 ×256 Be. Fakt. 184.004 214.563 226.953 154.062 140.624 159.664
Tabelle 7.3: CPU-Zeit und Beschleunigungsfaktor in Abh¨
angigkeit von Mehrgitteralgo-
rithmen und Gitterebenen mit QUICK Schema f¨
ur die rechteckige Nischen-
str¨
omung, gr¨
obstes Gitter 32 ×32
mod. FMG FMG-FAS
Gitter Re 100 400 1000 100 400 1000
CPU-Zeit 19.856 18.718 18.264 19.856 18.718 18.264
64 ×64 Be. Fakt. 1.000 1.000 1.000 1.000 1.000 1.000
CPU-Zeit 37.843 31.827 29.811 41.563 32.683 30.892
128 ×128 Be. Fakt. 14.930 15.466 14.917 13.594 15.060 14.395
CPU-Zeit 83.477 67.190 58.692 117.661 106.853 82.871
256 ×256 Be. Fakt. 116.069 131.040 129.040 82.348 76.786 91.215
Tabelle 7.4: CPU-Zeit und Beschleunigungsfaktor in Abh¨
angigkeit von Mehrgitteralgo-
rithmen und Gitterebenen mit QUICK Schema f¨
ur die rechteckige Nischen-
str¨
omung, gr¨
obstes Gitter 64 ×64
modifizierten Mehrgitteralgorithmus anhand einer um 27◦verscherten Nische untersucht.
In Abbildung 7.2a sind die verwendeten vier Mehrgitterstufen von 16 ×16 bis 128 ×128
dargestellt. Die Berechnung f¨
ur Re=1000 wird anstelle von QUICK und MUSCL mit Hilfe
des flux-blending Schema mit einem Faktor von 0.5 durchgef¨
uhrt.
In Tabelle 7.7 wird die CPU-Zeit des Eingitterverfahrens aufgelistet. F¨
ur diesen Testfall
ben¨
otigt das QUICK-Schema wieder ungef¨
ahr 10% weniger CPU-Zeit als das MUSCL-
Schema. Tabelle 7.8 zeigt CPU-Zeit und den Beschleunigungsfaktor des Mehrgitterverfah-
rens in Abh¨
angigkeit von Reynoldszahl, Konvektionsschema und Mehrgitteralgorithmus.
Zum Vergleich ist der von Lien und Leschziner (1994a) f¨
ur dieselbe Str¨
omung ermittel-
te Beschleunigungsfaktor angegeben. Auf den groben Gittern von 32 ×32 sowie 64 ×64
ben¨
otigen beide FMG-Methoden ungef¨
ahr gleiche CPU-Zeit. Auf dem feinen Gitter von
128 ×128 ist das modifizierte FMG deutlich effizienter (ungef¨
ahr 40%) als das Standard
FMG-FAS. Die Resultate zeigen, dass das modifizierte FMG ein hohes Potential auf dem
7.2. 2D laminare Str¨
omung in verzerrter Nische 75
Re mod. FMG FMG-FAS Lien (1994a)
Gitter 100 400 1000 100 400 1000 100 400 1000
16 ×16 1.000 1.000 1.000 1.000 1.000 1.000 1.0 0.9 0.8
32 ×32 2.934 2.215 1.839 2.826 2.128 1.601 3.9 2.5 1.8
64 ×64 12.212 11.024 9.333 11.192 8.986 6.603 19.9 10.6 6.2
128 ×128 60.801 65.806 59.905 49.177 49.610 40.992 43.7 26.5 16.8
Tabelle 7.5: Beschleunigungsfaktor in Abh¨
angigkeit von Mehrgitteralgorithmen und Git-
terebenen mit QUICK Schema f¨
ur die rechteckige Nischenstr¨
omung
(a) (b)
Abbildung 7.1: (a) Geometrie und vier Gitterebenen; (b) Stromlinien f¨
ur die rechteckige
Nischenstr¨
omung bei Re=400
(a) (b)
Abbildung 7.2: (a) Geometrie und vier Mehrgitterebenen; (b) Stromlinien f¨
ur die
Str¨
omung in einer verzerrten Nische bei Re=400
feinen Gitter besitzt. Dieses Verhalten entspricht dem f¨
ur die quadratische Nische und
orthogonalem Gitter ermittelten. Die Gitterschr¨
age hat offenbar keinen großen Einfluss
76 Kapitel 7. Validierung des Simulationsverfahrens
mod. FMG (MUSCL) mod. FMG (QUICK) HYBRID
Gitter Re=100 Re=1000 Re=100 Re=1000 Gitter Re=100 Re=1000
16 ×16 -0.193 -0.250 -0.191 -0.260 40 ×40 -0.202 -0.258
32 ×32 -0.206 -0.322 -0.206 -0.321 80 ×80 -0.209 -0.338
64 ×64 -0.211 -0.363 -0.211 -0.367 160 ×160 -0.212 -0.381
128 ×128 -0.212 -0.381 -0.212 -0.379 320 ×320 -0.213 -0.387
256 ×256 -0.213 -0.382 -0.213 -0.380 320 ×320 -0.213 -0.387
Tabelle 7.6: Vergleich der maximal auftretenden R¨
uckstr¨
omungsgeschwindigkeit auf der
vertikalen Nischenmittelachse f¨
ur unterschiedliche Schemata
MUSCL QUICK FLUX-B.
Gitter Re = 100 Re = 400 Re=100 Re=400 Re=1000
16 ×16 0.130 0.249 0.117 0.234 0.177
32 ×32 1.429 2.104 1.325 1.864 1.677
64 ×64 26.793 32.463 24.719 29.733 31.212
128 ×128 794.761 880.566 761.931 846.387 949.273
Tabelle 7.7: Eingitterverfahren: CPU-Zeit in Abh¨
angigkeit von Reynoldszahl und Kon-
vektionsschema f¨
ur die Str¨
omung in einer verzerrten Nische
auf die Effizienz des modifizierten Mehrgitteralgorithmus.
Im Vergleich zur Untersuchung von Lien und Leschziner (1994a) ist der Beschleunigungs-
faktor des jetzigen FMG-FAS f¨
ur das vierstufige Gitter bis zu 60% h¨
oher. Auf dem groben
Gitter ist er hingegen geringer. Abbildung 7.2b zeigt Stromlinien f¨
ur Re=400, die mit
MUSCL berechnet wurden.
7.3 2D laminare Str¨
omung ¨
uber eine zur¨
uckspringen-
de Stufe
Die Str¨
omung ¨
uber eine zur¨
uckspringende Stufe ist schematisch in Abbildung 7.3 darge-
stellt. Hinter der Stufe bildet sich eine große Rezirkulationszone aus. In der gew¨
ahlten
Konfiguration betr¨
agt die Stufenh¨
ohe Hgerade die H¨
alfte der Kanalh¨
ohe. Das Rechen-
gebiet hat eine L¨
ange von 10 Stufenh¨
ohen, um am Ausgang eine Neumannsbedingung
verwenden zu k¨
onnen. Am Eintritt wird ein parabolisches Geschwindigkeitsprofil vorgege-
ben. Die Reynoldszahl, gebildet mit der durchschnittlichen Geschwindigkeit und der Ein-
gangsh¨
ohe, betr¨
agt Re = 100. Um Vergleichbarkeit mit der Untersuchung von Orth (1991)
sicherzustellen, werden dieselben Relaxationsfaktoren, n¨
amlich 0.7 f¨
ur die Geschwindig-
keit und 0.2 f¨
ur die Druckkorrektur benutzt.
Tabelle 7.9 zeigt das Verhalten des Mehrgitterverfahrens auf dem gr¨
obsten Gitter von
7.3. 2D laminare Str¨
omung ¨
uber eine zur¨
uckspringende Stufe 77
mod. FMG (QUICK) FMG-FAS
QUICK Flux-B. QUICK Flux-B. Lien (1994a)
Gitter Re 100 400 1000 100 400 1000 100 400
CPU-Zeit 0.130 0.249 0.177 0.117 0.234 0.177 – –
16 ×16 Be. Fakt. 1.000 1.000 1.000 1.000 1.000 1.000 – –
CPU-Zeit 0.404 0.851 0.547 0.445 0.908 0.571 – –
32 ×32 Be. Fakt. 3.280 2.190 3.066 2.978 2.053 2.937 3.3 2.1
CPU-Zeit 1.827 2.639 1.888 2.102 3.398 2.766 – –
64 ×64 Be. Fakt. 13.530 11.267 16.532 11.760 8.750 11,284 15.5 11.7
CPU-Zeit 12.425 12.000 15.886 15.436 16.804 23.358 – –
128 ×128 Be. Fakt. 61.322 70.532 59.755 49.361 50.368 40.640 27.2 31.3
Tabelle 7.8: CPU-Zeit und Beschleunigungsfaktor f¨
ur die Simulation einer verzerrten Ni-
schenstr¨
omung mit verschiedenen Mehrgitterprozeduren, Gitterebenen und
Konvektionsschemata, gr¨
obstes Gitter 16 ×16
(a)
(b)
Abbildung 7.3: (a) Gitterverteilung auf dem gr¨
obsten Gitter; (b) Stromlinien f¨
ur die la-
minare Str¨
omung ¨
uber eine zur¨
uckspringende Stufe bei Re=100
40 ×20 Zellen. Die Simulation mit dem QUICK-Schema best¨
atigt die Ergebnisse der Ni-
schenstr¨
omungen. Dargestellt ist der Beschleunigungsfaktor des modifizierten FMG, der
bis zu 80% h¨
oher ist als der des FMG-FAS auf dem 160×80 Gitter. F¨
ur das feinste Gitter
lassen sich mit dem modifizierten FMG und vier Gitterebenen ein Beschleunigungsfaktor
von 167 und eine CPU-Zeit von 41 Sekunden erzielen.
78 Kapitel 7. Validierung des Simulationsverfahrens
CPU-Zeit Beschleunigungsfaktor
Gitter EG m. FMG FMG-FAS m. FMG FMG-FAS Orth (1991)
40 ×20 0.900 0.900 0.900 1.000 1.000 —
80 ×40 12.863 2.541 4.052 5.062 3.174 —
160 ×80 360.272 9.201 19.715 38.923 18.274 9.850
320 ×160 6840.364 40.793 — 167.684 — —
Tabelle 7.9: CPU-Zeit und Beschleunigungsfaktor in Abh¨
angigkeit von Mehrgitterproze-
dur und Konvektionsschema f¨
ur die Str¨
omung ¨
uber eine zur¨
uckspringende
Stufe bei Re=100, gr¨
obstes Gitter 40 ×20
Hier sollte erw¨
ahnt werden, dass sich mit dem FMG-FAS und vier Gitterebenen keine
konvergente L¨
osung erzielen l¨
asst. Ein Grund daf¨
ur liegt darin, dass die vom feineren Git-
ter restringierte Geschwindigkeit nicht zum parabolischen Profil am Eintritt passt. Dies
beeintr¨
achtigt m¨
oglicherweise die Stabilit¨
at der Standardprozedur des FMG-FAS. Außer-
dem gibt es einen großen Unterschied zwischen dem Beschleunigungsfaktor des FMG-FAS
und dem von Orth (1991). Die Gr¨
unde liegen m¨
oglicherweise in der Programmstruktur.
0 20 40 60 80 100 120
Zyklen
10-3
10-2
10-1
100
101
102
max. Residuum
MUSCL, Eingitter
MUSCL, mod. FMG, 3 Ebenen
MUSCL, FMG-FAS, 3 Ebenen
Abbildung 7.4: Residuenverlauf ¨
uber der Zyklenanzahl f¨
ur die Str¨
omung ¨
uber eine zur¨
uck-
springende Stufe mit 160 ×80 Gitterzellen bei Re=100
Tabelle 7.10 zeigt das Verhalten des modifizierten FMG auf dem gr¨
obsten Gitter mit
7.4. 2D turbulente Str¨
omung ¨
uber ein H¨
ugelmodell 79
MUSCL QUICK
CPU-Zeit CPU-Zeit
Gitter SG mod. FMG Be. Fakt. SG mod. FMG Be. Fakt.
80 ×40 14.620 14.620 1.000 12.683 12.683 1.000
160 ×80 368.645 35.857 10.281 360.272 34.533 10.433
320 ×160 7102.170 71.409 99.457 6840.364 71.242 96.016
Tabelle 7.10: CPU-Zeit und Beschleunigungsfaktor in Abh¨
angigkeit von Mehrgitterproze-
dur f¨
ur die Str¨
omung ¨
uber eine zur¨
uckspringende Stufe bei Re=100, gr¨
obstes
Gitter 80 ×40
80 ×40 Zellen. Die Rechnung wurde mit MUSCL bzw. QUICK durchgef¨
uhrt. Es gibt
leichte Unterschiede zwischen den Beschleunigungsfaktoren. Das QUICK-Schema ben¨
otigt
ein bisschen weniger Rechenzeit, w¨
ahrend MUSCL einen h¨
oheren Beschleunigungsfaktor
erreicht. Der in Abbildung 7.4 dargestellte Konvergenzablauf unterst¨
utzt die Beobachtung,
dass die modifizierte FMG-Prozedur ohne Restriktion der Geschwindigkeit eine stabile
und effiziente Konvergenz zeigt.
7.4 2D turbulente Str¨
omung ¨
uber ein H¨
ugelmodell
Das hier gew¨
ahlte Beispiel einer inkompressiblen, zweidimensionalen, turbulenten Innen-
str¨
omung ¨
uber ein H¨
ugelmodell wird vielfach zur Bewertung von Turbulenzmodellen ver-
wendet und soll hier Bewertung des Mehrgitterverfahrens verwendet werden. Die unter-
suchte Konfiguration ist der experimentellen Studie von Almeida u.a. (1993) entnommen.
Die mit der H¨
ugelh¨
ohe und der mittleren Einstr¨
omgeschwindigkeit gebildete Reynoldszahl
betr¨
agt Re=6000. Als Einstr¨
ombedingung wird eine voll entwickelte turbulente Kanal-
str¨
omung angenommen. Die Str¨
omung l¨
ost im Bereich der H¨
ugelspitze ab und es bildet
sich stromabw¨
arts eine Rezirkulationszone aus, die sich im Experiment bis x = 135 mm
erstreckt.
F¨
ur diese Konfiguration werden numerische Simulationen auf unterschiedlichen Gittern
unter Verwendung des k−²Modells mit Wandfunktion sowie dem LL k−²Modell
durchgef¨
uhrt. Das Mehrgitterverfahren basiert auf f¨
unf Stufen, von 25 ×10 bis 400 ×160
Gitterzellen ungleichm¨
aßiger Verteilung. F¨
ur die beiden Turbulenzmodelle werden zwei
unterschiedliche Gittersysteme verwendet. Die Parameter des Mehrgitterverfahrens wer-
den wie folgt eingestellt: Vorgl¨
attungszahl ν1= 4, Nachgl¨
attungszahl= 2, Iterationszahl
f¨
ur das gr¨
obste Gitter ν0= 8. Die Anwendung von FMG-FAS in Kombination mit dem
Standard k−²Modell liefert keine konvergente L¨
osung, wenn das gr¨
obste Gitter nur
25 ×10 Gitterzellen besitzt. Der dennoch ausgewertete Beschleunigungsfaktor f¨
ur den
modifizierten FMG-Algorithmus betr¨
agt 36, wie aus Tabelle 7.11 zu entnehmen ist.
Die Tabellen 7.12–7.14 zeigen das Verhalten des modifizierten FMG bzw. des Standard
80 Kapitel 7. Validierung des Simulationsverfahrens
Gitter SG mod. FMG Be. Faktor
25 ×10 0.377 0.377 1.000
50 ×20 3.171 2.273 1.395
100 ×40 30.563 13.610 2.246
200 ×80 772.776 72.007 10.732
400 ×160 12468.871 328.169 36.067
Tabelle 7.11: CPU-Zeit und Beschleunigungsfaktor des modifizierten FMG f¨
ur die turbu-
lente Str¨
omung ¨
uber ein 2D H¨
ugelmodell mit dem k−²Modell, gr¨
obstes
Gitter 25 ×10
FMG-FAS mit dem k−²Modell auf verschiedenen Gitterebenen. Das modifizierte FMG
erreicht auf vier Gitterebenen einen Beschleunigungsfaktor von 38.4, w¨
ahrend dieser beim
Standard FMG-FAS 12.1 betr¨
agt.
CPU-Zeit Be. Fakt.
Gitter SG mod. FMG FMG-FAS mod. FMG FMG-FAS
50 ×20 3.171 3.171 3.171 1.000 1.000
100 ×40 30.563 15.145 23.648 1.939 1.292
200 ×80 772.776 78.272 164.372 9.873 4.704
400 ×160 12468.871 324.583 1029.491 38.415 12.112
Tabelle 7.12: CPU-Zeit und Beschleunigungsfaktor verschiedenstufiger Mehrgitterproze-
duren f¨
ur die Str¨
omung ¨
uber ein H¨
ugelmodell mit dem k−²Modell, gr¨
obstes
Gitter 50 ×20
CPU-Zeit Be. Fakt.
Gitter SG mod. FMG FMG-FAS mod. FMG FMG-FAS
100 ×40 30.563 30.563 30.563 1.000 1.000
200 ×80 772.718 157.580 278.800 4.904 2.772
400 ×160 12468.871 599.381 1580.110 20.803 7.891
Tabelle 7.13: CPU-Zeit und Beschleunigungsfaktor verschiedenstufiger Mehrgitterproze-
duren f¨
ur die Str¨
omung ¨
uber ein H¨
ugelmodell mit dem k−²Modell, gr¨
obstes
Gitter 100 ×40
7.4. 2D turbulente Str¨
omung ¨
uber ein H¨
ugelmodell 81
CPU-Zeit Be. Fakt.
Gitter SG mod. FMG FMG-FAS mod. FMG FMG-FAS
200 ×80 772.718 772.718 772.718 1.000 1.000
400 ×160 12468.871 2233.317 3656.589 5.573 3.410
Tabelle 7.14: CPU-Zeit und Beschleunigungsfaktor zweistufiger Mehrgitterprozeduren f¨
ur
die Str¨
omung ¨
uber ein H¨
ugelmodell mit dem k−²Modell
Die Tabellen 7.15–7.17 zeigen die Rechenzeiten der beiden Mehrgitterverfahren im Ver-
gleich in Kombination mit dem dem LL k−²Modell. Wie aus den Ergebnissen zu erkennen
ist, hat der modifizierte Algorithmus ein besseres Verhalten als das Standardverfahren,
allerdings ist die Konvergenzrate nicht gleichbleibend wie bei den laminaren F¨
allen.
CPU-Zeit Be. Fakt.
Gitter SG mod. FMG FAS-FMG mod. FMG FAS-FMG
50 ×20 4.525 4.525 4.525 1.000 1.000
100 ×40 47.462 14.893 19.646 3.186 2.415
200 ×80 923.487 85.894 110.457 10.751 8.361
400 ×160 14361.002 396.117 774.637 36.254 18.539
Tabelle 7.15: CPU-Zeit und Beschleunigungsfaktor verschiedenstufiger Mehrgitterproze-
duren f¨
ur die Str¨
omung ¨
uber ein H¨
ugelmodell mit dem LL k−²Modell,
gr¨
obstes Gitter 50 ×20
CPU-Zeit Be. Fakt.
Gitter SG mod. FMG FAS-FMG mod. FMG FAS-FMG
100 ×40 44.590 44.590 44.590 1.000 1.000
200 ×80 923.487 158.329 280.230 5.833 3.295
400 ×160 14361.022 992.602 1255.192 14.468 11.441
Tabelle 7.16: CPU-Zeit und Beschleunigungsfaktor verschiedenstufiger Mehrgitterproze-
duren f¨
ur die Str¨
omung ¨
uber ein H¨
ugelmodell mit dem LL k−²Modell,
gr¨
obstes Gitter 100 ×40
82 Kapitel 7. Validierung des Simulationsverfahrens
CPU-Zeit Be. Fakt.
Gitter SG mod. FMG FAS-FMG mod. FMG FAS-FMG
200 ×80 923.487 923.487 923.487 1.000 1.000
400 ×160 14361.022 2459.859 4014.104 5.838 3.578
Tabelle 7.17: CPU-Zeit und Beschleunigungsfaktor zweistufiger Mehrgitterprozeduren f¨
ur
die Str¨
omung ¨
uber ein H¨
ugelmodell mit dem LL k−²Modell
0 20 40 60 80 100
Zyklen
10-3
10-2
10-1
100
101
max. Residuum
MUSCL, Eingitter
MUSCL, mod. FMG, 4 Ebenen
MUSCL, FMG-FAS, 4 Ebenen
Abbildung 7.5: Residuenverlauf ¨
uber der Zyklenzahl f¨
ur die Str¨
omung ¨
uber ein H¨
ugel-
modell mit dem LL k−²Modell
Abbildung 7.5 zeigt den Konvergenzablauf der Rechnung mit dem LL k−²Modell f¨
ur
Eingitterverfahren bzw. Mehrgitterverfahren. Es wird deutlich, dass das Maximalresidu-
um bei der modifizierten Methode schneller reduziert wird. Dieses wird m¨
oglicherweise
dadurch verursacht, dass das Flussfeld in diesem Fall nicht zu den restringierten Turbu-
lenzgr¨
oßen auf dem groben Gitter passt. Die turbulente Viskosit¨
at µt=ρCµk2/² auf dem
groben Gitter passt ebenfalls nicht zur restringierten turbulenten Energie kHbzw. der
Dissipationsrate ²H.
Stromlinien f¨
ur die Simulation mit dem k−²Modell sowie die horizontale Str¨
omungsge-
schwindigkeit an einem Schnitt von x=0 mm im Vergleich zu den experimentellen Daten
werden in Abbildung 7.6 dargestellt.
7.4. 2D turbulente Str¨
omung ¨
uber ein H¨
ugelmodell 83
x
y
-80. -40. 0. 40. 80. 120. 160. 200. 240. 280. 320.
0.
34.
68.
102.
136.
170.
-1.0 0.0 1.0 2.0 3.0 4.0
u
0
40
80
120
160
y
Experiment
Simulation
(a)
x
y
-80. -40. 0. 40. 80. 120. 160. 200. 240. 280. 320.
0.
34.
68.
102.
136.
170.
-1.0 -0.5 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0
u
0
20
40
60
80
100
120
140
160
y
Experiment
Simulation
(b)
x
y
-80. -40. 0. 40. 80. 120. 160. 200. 240. 280. 320.
0.
34.
68.
102.
136.
170.
-1.0 -0.5 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0
0
20
40
60
80
100
120
140
160
Experiment
Simulation
(c)
x
y
-80. -40. 0. 40. 80. 120. 160. 200. 240. 280. 320.
0.
34.
68.
102.
136.
170.
-1.0 -0.5 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0
0
20
40
60
80
100
120
140
160
Experiment
Simulation
(d)
Abbildung 7.6: Stromlinien f¨
ur die Str¨
omung ¨
uber einen H¨
ugel, als Ergebnisse von
Simulationen auf verschiedenen Gittern mit: (a) 25 ×10; (b) 50 ×20; (c)
100 ×40; (d) 200 ×80 Gitterzellen
84 Kapitel 7. Validierung des Simulationsverfahrens
7.5 3D laminare Nischenstr¨
omung
Um die modifizierte Mehrgittermethode f¨
ur dreidimensionale Simulationen bewerten zu
k¨
onnen, wird im folgenden die Str¨
omung in einem w¨
urfelf¨
ormigen Hohlraum (siehe Abbil-
dung 7.7a) untersucht. F¨
ur diese Konfiguration sind umfangreiche numerische Vergleichs-
untersuchungen, wie z.B. von Vanka (1986a), Orth (1991), Lien und Leschziner (1994a)
sowie Lilek u. a. (1997) verf¨
ugbar. F¨
ur die Simulation der Str¨
omung bei drei betrachteten
Reynoldszahlen (Re=100, 400, 1000, gebildet mit der Geschwindigkeit der oberen Wand
und der W¨
urfelkantenl¨
ange) werden zwei verschiedene vierstufige, kartesische Gittersy-
steme mit 323, 163, 83, 43sowie 403, 203, 103und 53Gitterzellen verwendet. Zur Wand
hin sind alle Gitter verdichtet. Konvektive Fl¨
usse werden mit Hilfe des QUICK-Schemas
erfasst.
u
z
y
x
(a) (b) y-z Ebene bei x=0.5
(c) x-y Ebene bei z=0.5 (d) x-z Ebene bei y=0.5
Abbildung 7.7: (a) Geometrie; (b)-(d) Str¨
omungsvektoren an unterschiedlichen Schnitt-
ebenen f¨
ur 3D laminare Nischenstr¨
omung
Abbildung 7.8 zeigt den logarithmisch aufgetragenen Residuenverlauf f¨
ur die Ein- und
Mehrgitterberechnung mit 403Volumenzellen bei Re=400. Bei der Eingitterberechnung
wird das Residuum zu Beginn des Iterationsablaufs relativ schnell reduziert, danach sinkt
die Konvergenzgeschwindigkeit jedoch deutlich. Insgesamt werden 83 Zyklen ben¨
otigt, um
das festgeschriebene Konvergenzkriterium (² < 10−3) zu erreichen. Im Fall der Mehrgit-
terberechnung bleibt die Konvergenzgeschwindigkeit konstant hoch.
7.5. 3D laminare Nischenstr¨
omung 85
0 10 20 30 40 50 60 70 80 90 100
Zyklen
10−3
10−2
10−1
100
max. Residuum
QUICK, Eingitter
QUICK, mod. FMG
Abbildung 7.8: Residuenverlauf ¨
uber der Zyklenzahl f¨
ur 3D laminare Nischenstr¨
omung
Rechenzeit und Beschleunigungsfaktoren des modifizierten FMG sowie der Untersuchun-
gen von Lien und Leschziner (1994a) bzw. Lilek u. a. (1997) sind in den Tabellen 7.18 –
7.19 dargestellt. F¨
ur Re=100 l¨
asst sich mit dem modifizierten FMG ein Beschleunigungs-
faktor von 13.05 f¨
ur Gitter mit 323Volumenzellen bzw. von 20.595 f¨
ur Gitter mit 403
Volumenzellen erreichen. Diese liegen deutlich ¨
uber den Vergleichswerte in den Arbeiten
von Lien und Leschziner (1994a) und Lilek u. a. (1997).
SG mod. FMG
Gitter Re=100 Re=400 Re=1000 Re=100 Re=400 Re=1000
3231217.889 1080.142 990.745 93.322 95.450 113.866
4033847.424 3335.171 2903.639 186.813 190.468 190.213
Tabelle 7.18: Erforderliche CPU-Zeit f¨
ur Eingitterverfahren und modifiziertes FMG auf
vier Gitterebenen f¨
ur 3D laminare Nischenstr¨
omung
mod. FMG Lien (1994a) Lilek (1997)
Gitter Re=100 Re=400 Re=1000 Re=100 Re=100 Re=1000
32313.050 11.316 8.701 6.36 6.67 2.55
40320.595 17.510 15.265 10.08 – –
Tabelle 7.19: Vergleich der Beschleunigungsfaktoren f¨
ur 3D laminare Nischenstr¨
omung
Abbildungen 7.7a, b, c zeigen Str¨
omungsvektoren bei Re=1000 in drei ausgesuchten
Schnittebenen. In der Symmetrieebene (y=0.5) wird die Struktur durch einen großen
Wirbel in der Mitte und kleinere Wirbel in den unteren Ecken charakterisiert. Das kom-
plexe sekund¨
are Wirbelsystem ist an der Zirkulationszone in den Schnittebenen x=0.5
und z=0.5 deutlich zu erkennen.
86 Kapitel 7. Validierung des Simulationsverfahrens
7.6 3D Str¨
omung durch einen Rohrkr¨
ummer
In diesem Abschnitt soll das modifizierte FMG Verfahren am Beispiel der Str¨
omung durch
einen dreidimensionalen Rohrkr¨
ummer mit starkem Sekund¨
arwirbel untersucht werden.
F¨
ur diese Str¨
omungskonfiguration existieren experimentelle Ergebnisse als Validierungs-
basis sowohl f¨
ur laminare als auch f¨
ur turbulente Str¨
omungen. Hiermit wird die Leistung
des Mehrgitterverfahrens in dreidimensionalen Str¨
omungen mit starkem Sekund¨
arwirbel
demonstriert. Abbildung 7.9 zeigt die Rohrgeometrie des Experiments von Taylor u. a.
(1982). Um einen 90-Grad-Kr¨
ummer mit einem Mittelradius von 92 mm gruppieren sich
die aufw¨
arts und abw¨
arts gerichteten Tangenten mit 0.3 mbzw. 2.0 mL¨
ange. Die Dimen-
sion des Querschnitts im Kr¨
ummer (40 mm ×40 mm) entspricht dem der Tangenten.
Die Volumengeschwindigeit Ubbetragt 1.98 cm/s f¨
ur die laminare bzw. 1.00 m/s f¨
ur die
turbulente Str¨
omung, entsprechend einer Reynoldszahl von 790 bzw. 40000, gebildet mit
Einstr¨
omgeschwindigkeit Ubund den hydraulischen Durchmesser D. Das Geschwindig-
keitsprofil am Rohreingang wird entsprechend dem Experiment als Blockprofil angenom-
men.
i
ro
o
y
r
00.1
0.3
0.5
0.9
0.7
*
1
X = y/D
X = x/D
*
*
r
z
o i o
*
x
2 m
0.3 m
0
0
θ=0
θ=90
z
r =(r-r )/(r - r )
z = 2z/D
H
H
i
o
r = 0.112m
r = 0.072m
Abbildung 7.9: Geometrie des quadratischen Rohrkr¨
ummers, Taylor u. a. (1982)
Das Gitter mit 128 ×32 ×32 Zellen wird in axialer Richtung im Bereich der Kr¨
ummung
verdichtet. Der Zulauf und der senkrecht stehende Ausfluss sind in Hauptstr¨
omungsrich-
tung relativ grob aufgel¨
ost. In axialer Richtung ist es in 8 Bl¨
ocke eingeteilt, woraus drei
Gitterebenen mit 16 ×32 ×32, 8 ×16 ×16 und 4 ×8×8 Zellen gebildet werden.
7.6. 3D Str¨
omung durch einen Rohrkr¨
ummer 87
Laminare Str¨
omung
Zun¨
achst werden Simulationen einer laminaren Durchstr¨
omung des Rohrkr¨
ummers bei
Re=790 durchgef¨
uhrt, wobei zum Vergleich ein Eingitterverfahren sowie das modifizierte
FMG zum Einsatz kommen. Die Behandlung der Konvektion erfolgt mit einem TVD-
MUSCL Schema. Tabelle 7.20 zeigt das Verhalten des Mehrgitterverfahrens mit drei Git-
terebenen.
CPU-Zeit Be. Fakt.
Gitter SG mod. FMG mod. FMG
128 ×32 ×32 10972.07 2100.07 5.22
Tabelle 7.20: CPU-Zeit und Beschleunigungsfaktor f¨
ur die Berechnung der laminaren
Rohrkr¨
ummerstr¨
omung bei Re=790
0 20 40 60 80 100 120
Zyklen
10-4
10-3
10-2
10-1
100
101
102
max. Residuum
MUSCL, Eingitter
MUSCL, mod. FMG
Abbildung 7.10: Residuenverlauf ¨
uber der Zyklenzahl f¨
ur laminare Rohrkr¨
ummer-
str¨
omung bei Re=790
Abbildung 7.10 zeigt den logarithmisch aufgetragenen Residuenverlauf f¨
ur die Ein- und
Mehrgitterberechnungen. Beim Mehrgitterverfahren wird das Residuum zu Beginn sehr
schnell reduziert, danach sinkt die Konvergenzgeschwindigkeit und das Residuum f¨
allt
linear. Trotzdem wird eine Leistungssteigerung um den Faktor 5.22 erreicht.
In Abbildungen 7.11–7.12 sind die auf dem feinsten Gitter ermittelten axialen Geschwin-
digkeitsprofile im Vergleich zu den Messungen dargestellt. Die numerischen Ergebnisse
stimmen sehr gut mit den experimentellen ¨
uberein.
88 Kapitel 7. Validierung des Simulationsverfahrens
0.00.51.01.50.00.51.01.50.00.51.01.50.00.51.01.50.00.51.01.52.0
U/Ub
1.0
0.8
0.6
0.4
0.2
0.0
0.8
0.6
0.4
0.2
z*
Simulation (mod. FMG)
r*=.1 r*=.3 r*=.5 r*=.7 r*=.9
0.0
0.8
0.6
0.4
0.2
0.0
Experiment Taylor et al
XH=-.25
θ=300
θ=600
Abbildung 7.11: Axiale Geschwindigkeitsprofile U/Uban den Rohrquerschnitten XH=
−0.25, θ= 30◦und θ= 60◦f¨
ur laminare Rohrkr¨
ummerstr¨
omung bei
Re=790
Abbildung 7.13 zeigt die Sekund¨
arstr¨
omung in verschiedenen Rohrquerschnitten. Deut-
lich ist zu erkennen, dass in der am Anfang weitgehend zweidimensionalen Str¨
omung
durch die Rohrkr¨
ummung starke axiale Wirbel erzeugt werden. Das Maximum der Axi-
algeschwindigkeit wandert sehr schnell von der Rohrmitte zur Außenseite des Kr¨
ummers.
Die Sekund¨
arstr¨
omung dagegen bewegt sich entlang der Rohrseitenw¨
ande in Richtung
der Innenwand und entlang der Symmetrieebene wieder zur¨
uck zur Außenwand. Diese
großr¨
aumige axiale Wirbelstruktur teilt sich im hinteren Teil des Kr¨
ummers in zwei Wir-
bel und setzt sich in den senkrecht stehenden Rohrteil hinein fort.
7.6. 3D Str¨
omung durch einen Rohrkr¨
ummer 89
0.00.51.01.50.00.51.01.50.00.51.01.50.00.51.01.50.00.51.01.52.0
U/Ub
1.0
0.8
0.6
0.4
0.2
0.0
0.8
0.6
0.4
0.2
z*
Simulation (mod. FMG)
r*=.1 r*=.3 r*=.5 r*=.7 r*=.9
0.0
0.8
0.6
0.4
0.2
0.0
Experiment Taylor et al
XH=2.5
XH=.25
θ=77.50
Abbildung 7.12: Axiale Geschwindigkeitsprofile U/Uban den Rohrquerschnitten θ=
77.5◦,XH= 0.25 und XH= 2.5 f¨
ur laminare Rohrkr¨
ummerstr¨
omung
bei Re=790
90 Kapitel 7. Validierung des Simulationsverfahrens
(a) (b)
(c) (d)
(e) (f)
Abbildung 7.13: Richtungsvektoren der Sekund¨
arstr¨
omung in den Rohrquerschnitten (a)
XH=−0.25, (b) θ= 30◦, (c) θ= 60◦, (d) θ= 77.5◦, (e) XH= 0.25 und
(f) XH= 2.5 f¨
ur laminare Rohrkr¨
ummerstr¨
omung
7.6. 3D Str¨
omung durch einen Rohrkr¨
ummer 91
Turbulente Str¨
omung
Die turbulente Str¨
omung durch den oben beschriebenen Rohrkr¨
ummer bei einer Reynolds-
zahl von Re=40000 wird mit dem LL k−²Modell berechnet. Die Vorgl¨
attungszahl ν1
betr¨
agt 2 und die Nachgl¨
attungszahl ν2wird auf 1 gesetzt. Auf dem gr¨
obsten Gitter
ist die Gl¨
attungszahl ν0= 8. Tabelle 7.21 zeigt CPU-Zeit und Beschleunigungsfaktor
f¨
ur eine Simulation auf drei Gitterebenen und einem Konvergenzkriterium von 10−5. Bei
dieser komplizierten turbulenten Str¨
omung l¨
asst sich ein Beschleunigungsfaktor von 6.6
erreichen.
CPU-Zeit Be. Fakt.
Gitter SG mod. FMG mod. FMG
128 ×32 ×32 24444.58 3690.31 6.62
Tabelle 7.21: CPU-Zeit und Beschleunigungsfaktor f¨
ur die Simulation einer turbulenten
Rohrkr¨
ummerstr¨
omung mit dem LL k−²Modell bei Re=40000
Abbildung 7.14 zeigt den Residuenverlauf f¨
ur die Berechnung der turbulenten Str¨
omung
bei Re=40000 mit dem LL k−²Modell und TVD-MUSCL Schema. Wie in Abbildungen
7.15 – 7.16 zu erkennen ist die ¨
Ubereinstimmung der Geschwindigkeitsprofile zwischen
Rechnung und den Messungen sehr gut.
0 50 100 150 200 250 300
Zyklen
10-4
10-3
10-2
10-1
100
101
102
max Residual
MUSCL, Eingitter
MUSCL, mod. FMG
Abbildung 7.14: Residuenverlauf ¨
uber der Zyklenzahl f¨
ur turbulente Rohrkr¨
ummer-
str¨
omung mit dem LL k−²Modell und MUSCL bei Re=40000
In Abbildung 7.17 sind Richtungsvektoren der Sekund¨
arstr¨
omung in verschiedenen Rohr-
querschnitten dargestellt. Die Sekund¨
arstr¨
omung ist etwas weniger stark ausgepr¨
agt als
92 Kapitel 7. Validierung des Simulationsverfahrens
bei der laminaren Str¨
omung. Die Position der maximalen Axialgeschwindigkeit liegt am
Anfang des Kr¨
ummers nahe der Innenwand. Sie wandert innerhalb des Kr¨
ummungsbe-
reichs langsam nach außen. Am Ausgang des Kr¨
ummungsteils liegt sie etwa in der Mitte
des Rohres.
0.00.51.01.50.00.51.01.50.00.51.01.50.00.51.01.50.00.51.01.52.0
U/Ub
1.0
0.8
0.6
0.4
0.2
0.0
0.8
0.6
0.4
0.2
z*
Simulation k-ε (mod. FMG)
r*=.1 r*=.3 r*=.5 r*=.7 r*=.9
0.0
0.8
0.6
0.4
0.2
0.0
Experiment Taylor et al
XH=-.25
θ=300
θ=600
Abbildung 7.15: Axiale Geschwindigkeitsprofile U/Uban den Rohrquerschnitten XH=
−0.25, θ= 30◦und θ= 60◦f¨
ur die turbulente Rohrkr¨
ummerstr¨
omung
bei Re=40000
Das Wachstum der turbulenten Grenzschicht in dem Zulauf ist langsamer als bei lami-
narer Str¨
omung. Deshalb ist bei XH=−0.25 die Grenzschicht d¨
unner, wobei sich die
Kernstr¨
omung im Vergleich zum laminaren Fall weiter in Richtung zur inneren Wand-
oberfl¨
ache bewegt. Sie bleibt an dieser Stelle, ohne Einfluss auf die Str¨
omung geringen
Impulses die sich bis θ= 60◦auf dieser Seite akkumuliert. Danach ist die schnelle Ausbil-
dung einer Region niedriger Hauptstr¨
omungskomponente an der inneren Wand bis zum
Ausgang zu beobachten. Dies korrespondiert mit einer Verlagerung der Kernstr¨
omung zur
¨
außeren Wand.
7.6. 3D Str¨
omung durch einen Rohrkr¨
ummer 93
0.00.51.01.50.00.51.01.50.00.51.01.50.00.51.01.50.00.51.01.52.0
U/Ub
0.8
0.6
0.4
0.2
0.0
1.0
0.8
0.6
0.4
0.2
0.0
z*
Computation k-ε (mod. FMG)
r*=.1 r*=.3 r*=.5 r*=.7 r*=.9
0.8
0.6
0.4
0.2
0.0
Experiment Taylor et al
XH=2.5
XH=.25
θ=77.50
Abbildung 7.16: Axiale Geschwindigkeitsprofile U/Ubin den Rohrquerschnitten θ= 77.5◦,
XH= 0.25 und XH= 2.5 f¨
ur die turbulente Rohrkr¨
ummerstr¨
omung bei
Re=40000
94 Kapitel 7. Validierung des Simulationsverfahrens
(a) (b)
(c) (d)
(e) (f)
Abbildung 7.17: Richtungsvektoren der Sekund¨
arstr¨
omung an den Rohrquerschnitten (a)
XH=−0.25, (b) θ= 30◦, (c) θ= 60◦, (d) θ= 77.5◦, (e) XH= 0.25 und
(f) XH= 2.5 f¨
ur die turbulente Rohrkr¨
ummerstr¨
omung bei Re=40000
Kapitel 8
Simulation komplexer Str¨
omungen
In diesem Kapitel wird das Simulationsverfahren beispielhaft auf zwei industriell relevante
Anwendungsf¨
alle angewendet, wobei zun¨
achst Ergebnisse f¨
ur eine Stufenk¨
orperflamme
unter Verwendung des Flamelet–Modells mit verschiedenen Turbulenzmodellen gezeigt
werden. Weiterhin wird die Simulation der instation¨
aren Umstr¨
omung eines generischen
KFZ-Außenspiegels mit einem Eingleichungsmodell bzw. einer DES durchgef¨
uhrt.
8.1 CH4/H2Stufenk¨
orper-Diffusionsflamme
Modelliert wird ein generischer Stufenk¨
orper-Brenner, der große ¨
Ahnlichkeit zu prakti-
schen Brennkammern hat, wie sie in vielen industriellen Bereichen verwendet werden.
Der Brenner wird in einem Koaxialstrom zentriert und besteht im allgemeinen aus ei-
nem kreisf¨
ormigen Stufenk¨
orper mit einer D¨
use f¨
ur den Brennstoff in der Mitte. Ein
kompliziertes Str¨
omungssystem entsteht hinter dem K¨
orper, wobei sich eine bzw. zwei
R¨
uckstr¨
omungszonen ausbilden. Diese m¨
ussen gen¨
ugend heiße Gase zur Flammenstabi-
lisierung produzieren. Die Konfiguration dient oft als Modellproblem, weil sie viele Ei-
genschaften mit praktischen Brennkammern gemeinsam hat und dabei verh¨
altnism¨
aßig
einfache und gut definierte Randbedingungen aufweist. Sie eignet sich besonders gut f¨
ur
die Studie von turbulent-chemischen Interaktion und ist dem international Workshop
on Measurement and Computation of Turbulent Nonpremixed Flames (TNF Workshop,
www.ca.sandia.gov/tdf/Workshop.html) entnommen. Experimentelle Untersuchung wur-
de an der Universit¨
at Sydney/Sandia durchgef¨
uhrt.
Abbildung 8.1 gibt Auskunft ¨
uber die Simulationsparameter. Die Anordnung besteht aus
einem Brennstoffstrahl mit einem Durchmesser von 3.6 mm, der von einem Zylinder mit ei-
nem Durchmesser von D=50 mm umgeben ist. Außen str¨
omt Koaxialluft um den Zylinder.
Die Messwerte f¨
ur Vergleich sind unter der Adresse http://www.mech.eng.usyd.edu.au/
research/energy zum Download verf¨
ugbar. Die f¨
ur den Vergleich herangezogenen experi-
mentellen Daten stammen vom Fall B4F3A. In diesem Fall besteht der Brennstoff volu-
metrisch aus 50% CH4und 50% H2. Der st¨
ochiometrische Mischungsbruch ist Zst = 0.050.
Die mittlere axiale Geschwindigkeit des Brennstoffstrahls betragt 118 m/s und die mittlere
axiale Geschwindigkeit der Koaxialluft ist 40 m/s.
95
96 Kapitel 8. Simulation komplexer Str¨
omungen


r
(0,0)
1D 6D
x
1.8 mm
20 D = 1 m
Ausgang
25 mm (D=50 mm)
Symmetrie
Brennstoff (50 % CH4, 50% H2)
Druckrandbedingung
Brennstoff
Luft
Abbildung 8.1: Brennerkonfiguration und numerische Randbedingungen f¨
ur Simulation
der Stufenk¨
orper-Diffusionsflamme
Die Flamelet-Bibliothek h¨
angt nun nicht nur von den bedingt gemittelten skalaren Dissi-
pationsraten eχZsondern auch vom Druck ab. Mit der Betrachtung des Strahlungsverlusts
werden die Flamelet-Gleichungen in eχZ-Raum (0.01 – 60.0) und in Druckraum (0.995
– 1.015 bar) gel¨
ost. Die turbulente kinetische Energie kund ihre Dissipationsrate ²am
Einstr¨
omrand werden ¨
uber konstante Verteilungen spezifiziert, die einer Intensit¨
at von
Tu= 3% und νt= 10νentsprechen.
Abbildung 8.2(a) zeigt das mit dem EASM ermittelte Str¨
omungsfeld der axialsymme-
trischen Konfiguration. Hinter dem Stufenk¨
orper gibt es eine R¨
uckstr¨
omungszone. Dies
produziert gen¨
ugend heiße Gase, um die Flamme zum Brenner hin zu stabilisieren. In Ab-
bildung 8.2 (b) sind die Messpositionen skizziert, entlang derer die numerischen Resultate
mit den experimentellen Daten verglichen werden.
Zuerst wird der Einfluss des Turbulenzmodells, bzw. Anisotropieparameters untersucht.
Dazu werden die numerisch ermittelten Profile der axialen Geschwindigkeit, der Tempera-
tur, der Massenbr¨
uche der Hauptprodukte H2O und CO2mit den experimentellen Daten
an der Position x/D=1.3 in Abbildung 8.3 verglichen. Hier bezeichnen die durchgezogenen
Linien Ergebnisse des Standard k−²Modells, die gestrichenen Linien Ergebnisse des LEA
k−²Modells, und die Punkt-gestrichenen Linien Resultate des EASMs. Im allgemeinen
ergeben das EASM und das LEA k−²Modell ¨
ahnliche Ergebnisse. An der Innenseite
f¨
ur r/R <0.15 ergibt das LEA k−²Modell sogar bessere Geschwindigkeitsverteilungen
als das EASM. Im Gebiet von 0.45 <r/R <0.9 kann eine deutliche Abweichung beim
Standard k−²Modell von den ermittelten Ergebnissen des LEA k−²Modells und des
EASMs sowie von den experimentellen Daten beobachtet werden. Die Temperaturvertei-
lung bliebt im Fall des Standard k−²Modells weit unter den experimentellen Daten
zur¨
uck. Auch alle anderen Modelle k¨
onnen den experimentellen Maximalwert der Tempe-
8.1. CH4/H2Stufenk¨
orper-Diffusionsflamme 97
(a)
x/D=.1 x/D=.6
x/D=.9
x/D=1.3
x/D=1.8
x
r
r/R=1.2
(b)
Abbildung 8.2: (a) Stromlinien hinter dem Stufenk¨
orper, Ergebnis einer Simulation mit
EASM; (b) Messpositionen f¨
ur die Stufenk¨
orper-Diffusionsflamme (Syd-
ney/Sandia)
ratur von T≈1800K nicht richtig reproduzieren. Die Position des Maximums wird jedoch
durch das EASM bzw. das LEA k−²Modell richtig vorhergesagt, w¨
ahrend die durch das
k−²Modell ermittelte Position des Maximums (ca. r/R=0.5) deutlich von der experi-
mentellen Position (ca. r/R=0.7) abweicht. Die Massenbr¨
uche von H2O werden durch das
LEA k−²Modell und das EASM befriedigend erfasst. Es gibt wieder eine relativ große
Abweichung beim k−²Modell. Die ¨
Ubereinstimmung der ermittelten Massenbr¨
uche von
CO2mit den experimentellen Daten ist unbefriedigend. Eine große Abweichung tritt nicht
nur beim Maximum, sondern auch in der N¨
ahe der Achse auf. Es wird offensichtlich, dass
das Turbulenzmodell einen dominanten Einfluss hat, und dass, wenn die Vorhersage der
Str¨
omungsgeschwindigkeit bereits defizit¨
ar ist, auch die Ergebnisse f¨
ur die Temperatur
und die Massenbr¨
uche nicht mit den experimentellen Daten ¨
ubereinstimmen.
Das Ph¨
anomen, dass es kaum Unterschiede zwischen den L¨
osungen des LEA k−²Mo-
dells und des EASMs gibt, kann durch theoretische Analyse begr¨
undet werden. Da die
Str¨
omung nicht verdrallt ist, lautet die U-Impulsgleichung in Zylinderkoordinaten wie
98 Kapitel 8. Simulation komplexer Str¨
omungen
0.0 0.3 0.6 0.9 1.2
r/R
0
0.02
0.04
0.06
0.08
0.1
0.12
0.14
0.16
Massenbruch H2O
Exp. Sydney/Sandia
k-ε
LEA k-ε
EASM
-10
0
10
20
30
40
50
60
U [m/s]
0.0 0.3 0.6 0.9 1.2
r/R
0
0.01
0.02
0.03
0.04
0.05
0.06
0.07
Massenbruch CO2
300
600
900
1200
1500
1800
2100
T [K]
x/D=1.3
x/D=1.3
x/D=1.3
x/D=1.3
Abbildung 8.3: Einfluss des Anisotropieparameters bei den k−²Modellen: Ermittelte
Geschwindigkeit, Temperatur, Massenbr¨
uche von H2O und CO2sowie ex-
perimentelle Daten f¨
ur die Stufenk¨
orper-Diffusionsflamme
folgt:
∂(rρeueu)
∂x +∂(rρeveu)
∂r =∂(rp)
∂x +∂
∂x µrµ∂eu
∂x+∂
∂r µµr∂eu
∂x−∂rρu0u0
∂x −∂(rρu0v0)
∂r
In diesem Fall dominiert die Scherung, die f¨
ur großen radialen Gradient der Axialgeschwin-
digkeit, d.h. ∂
e
u
∂r >> ∂
e
u
∂x ,∂
e
v
∂x und ∂
e
v
∂r , im Folgenden abgesch¨
atzt werden soll. Die Kompo-
nenten des Deformationsgeschwindigkeitstensors bzw. des Drehtensors ergeben sich n¨
ahe-
rungsweise wie folgt:
S=
0S12 0
S12 0 0
0 0 0
, W =
0W12 0
−W12 0 0
0 0 0
,(8.1)
wobei
S12 =1
2
∂eu
∂r , W12 =1
2
∂eu
∂r .
Daraus k¨
onnen kann Bij wie folgt berechnet werden:
Bij =k
²[β2(SikWkj +SjkWki)−β3(SikSkj −SlkSklδij)] ,
8.1. CH4/H2Stufenk¨
orper-Diffusionsflamme 99
und dann kann die Scherspannung wie folgt bestimmt werden:
−ρu0v0∼2µtS12
|{z}
LTerm
+ 2µtB12
|{z}
NLTerm
.
Durch theoretische Analyse l¨
asst sich zeigen, dass im Vergleich zum linearen Term (LTerm)
der nichtlinearen Term (NLTerm) viel kleiner ist, d.h. die Spannung wird sehr stark mit
dem linearen Term gekoppelt:
−ρu0v0∼2µtS12 ∼2µt
∂eu
∂r .
Weil der Term ∂(rρu0u0)
∂x viel kleiner als ∂(rρu0v0)
∂r ist, wird der Einfluss der Turbulenz auf die
Geschwindigkeit stark durch den linearen Term bestimmt. Beim Flamelet-Modell werden
die Temperatur und Massenbr¨
uche durch die bedingt gemittelte skalare Dissipationsra-
te eχZund die β-f¨
ormige PDF bestimmt, die durch die Turbulenzgr¨
oße kund ²sowie
den Mischungsbruch e
Zund seine Varianz g
Z002definiert werden, welche wiederum durch
die Transportgleichungen gel¨
ost werden. Deshalb sind die Resultate f¨
ur Geschwindigkeit,
Temperatur und Massenbr¨
uche f¨
ur das LEA k−²Modell und das EASM ¨
ahnlich. Die
Abweichung zwischen den Ergebnissen der k−²und LEA k−²Modelle wird durch
die verschiedenen Werte des Anisotropieparameters Cµverursacht. Der variable Aniso-
tropieparameter des LEA k−²Modells bzw. des EASMs wird durch eine Funktion der
Geschwindigkeitsgradienten berechnet. Es besitzt entscheidenden Einfluss bei der Simu-
lation scherdominierter Str¨
omungen.
0.0 0.3 0.6 0.9 1.2
r/R
0.00
0.10
0.20
0.30
0.40
0.50
Exp. Sydney/Sandia
k-ε
LEA k-ε
EASM
0.0 0.3 0.6 0.9 1.2
r/R
0.00
0.20
0.40
0.60
0.80
mixture fraction
x/D=0.9 x/D=1.3
Abbildung 8.4: Ermittelte Mischungsbr¨
uche und experimentelle Daten an den axialen Po-
sitionen x/D=0.9 und 1.3 f¨
ur Stufenk¨
orper-Diffusionsflamme
Im Folgenden werden die Mischungsbr¨
uche, die axialen Geschwindigkeiten, die Tempera-
turen und die Massenbr¨
uche von H2O, CO2sowie H2mit den experimentellen Daten an
zwei unterschiedlichen, ausgew¨
ahlten axialen Positionen verglichen. Werden Vergleiche an
anderen Positionen durchgef¨
uhrt, so best¨
atigen sie lediglich die im Folgenden beschrie-
benen Beobachtungen. Auf der Symmetrieachse gibt es keinen Vergleich, da dort keine
100 Kapitel 8. Simulation komplexer Str¨
omungen
0.0 0.3 0.6 0.9 1.2
r/R
-20
0
20
40
60
80
0.0 0.3 0.6 0.9 1.2
r/R
-30
0
30
60
90
120
U [m/s]
Exp. Sydney/Sandia
k-ε
LEA k-ε
EASM
x/D=1.3x/D=0.9
Abbildung 8.5: Ermittelte Geschwindigkeiten und experimentelle Daten an den axialen
Positionen x/D=0.9 und 1.3 f¨
ur Stufenk¨
orper-Diffusionsflamme
experimentellen Daten vorhanden sind.
Abbildung 8.4 zeigt die ermittelten Profile der Mischungsbr¨
uche an den axialen Positionen
x/D=0.9 und 1.3. Nah dem Strahlausgang f¨
ur x/D ≤0.9 sind die numerischen Resultate
aller Turbulenzmodelle relativ ¨
ahnlich und die ¨
Ubereinstimmung mit den experimentellen
Daten ist gut. Stromabw¨
arts f¨
ur x/D >0.9 wird ein st¨
arkerer Einfluss der Turbulenzmo-
delle beobachtet. An der Außenseite r/R >0.8 ergeben alle Turbulenzmodelle eine gute
¨
Ubereinstimmung mit den Messwerten, da in diesem Gebiet die Luftstr¨
omung dominie-
rend ist und keine Verbrennung stattfindet.
0.0 0.3 0.6 0.9 1.2
r/R
400
800
1200
1600
2000
0.0 0.3 0.6 0.9 1.2
r/R
400
800
1200
1600
2000
T[K]
Exp. Sydney/Sandia
k-ε
LEA k-ε
EASM
x/D=0.9 x/D=1.3
Abbildung 8.6: Ermittelte Temperatur und experimentelle Daten an den axialen Positio-
nen x/D=0.9 und 1.3 f¨
ur Stufenk¨
orper-Diffusionsflamme
Abbildung 8.5 zeigt den Vergleich der ermittelten axialen Geschwindigkeitsprofile mit den
Messwerten. Auf der Achse wird die Geschwindigkeitsverteilung mit allen Turbulenzmo-
dellen gut wiedergegeben. Bei x/D=0.9 sind die numerischen Ergebnisse nahezu gleich und
stimmen gut mit den experimentellen Profilen ¨
uberein. Auf der Achse liefern alle Turbu-
8.1. CH4/H2Stufenk¨
orper-Diffusionsflamme 101
lenzmodelle besonders das LEA k−²Modell h¨
ohere Geschwindigkeiten als die Messungen.
Stromabw¨
arts f¨
ur x/D=1.3 reproduzieren das LEA k−²Modell und das EASM mit den
Messwerten vergleichbare Geschwindigkeitsprofile, w¨
ahrend es im Gebiet von 0.45 <r/R
<0.8 eine deutliche Abweichung f¨
ur das k−²Modell gibt.
0.0 0.3 0.6 0.9 1.2
r/R
0.00
0.04
0.08
0.12
0.16
0.20
0.0 0.3 0.6 0.9 1.2
r/R
0.00
0.04
0.08
0.12
0.16
0.20
mass fraction H2O
Exp. Sydney/Sandia
k-ε
LEA k-ε
EASM
x/D=0.9 x/D=1.3
Abbildung 8.7: Ermittelte Massenbr¨
uche von H2O und experimentelle Daten an den axia-
len Positionen x/D=0.9 und 1.3 f¨
ur Stufenk¨
orper-Diffusionsflamme
Die Temperaturprofile werden in Abbildung 8.6 dargestellt. Im allgemeinen stimmen die
Temperaturprofile f¨
ur das LEA k−²Modell und das EASM mit den Messwerten ¨
uberein,
w¨
ahrend das k−²Modell relativ große Abweichungen im Gebiet von 0.45 <r/R <0.8 er-
gibt. Stromabw¨
arts ist die Maximaltemperatur f¨
ur die Profile bei x/D=0.9 und 1.3 kleiner
als die experimentellen Daten, liegen aber f¨
ur das LEA k−²Modell und das EASM an
den richtigen Stellen. Die Positionen des Temperaturmaximums (ca. r/R = 0.75) werden
noch nicht richtig durch das k−²Modell repr¨
asentiert. In den numerischen Simulationen
ergeben alle Turbulenzmodelle niedrigere H¨
ochstwerte. M¨
oglicherweise ist die Ungenauig-
keit des verwendeten Strahlungsmodells eine Ursache daf¨
ur, weil dder Strahlungsverlust
hier wahrscheinlich zu hoch ist.
Abbildung 8.7 zeigt die ermittelten radialen Massenbruchprofile des Hauptprodukts H2O
und die experimentellen Daten. Das LEA k−²Modell und das EASM ergeben eine hin-
reichende ¨
Ubereinstimmung mit den Messwerten, w¨
ahrend sich beim k−²Modell relative
große Abweichungen von den Messwerten im Gebiet von 0.45 <r/R <0.8 ergeben. Bei
x/D=0.9 sind die numerischen Ergebnisse des LEA k−²Modells und des EASMs etwas
h¨
oher als die Messwerte, dagegen sind die Resultate des k−²Modells etwas niedriger als
die experimentellen Daten. An der Außenseite r/R >1.0 werden die Massenbr¨
uche durch
alle Turbulenzmodelle untersch¨
atzt. Im allgemeinen k¨
onnen jedoch gute ¨
Ubereinstimmung
mit den experimentellen Daten durch das LEA k−²Modell und das EASM erzielt werden.
In Abbildung 8.8 werden die Massenbruchprofile des Hauptprodukts CO2dargestellt. Die
CO2-Verteilung ist nicht so gut wie die bei H2O, denn die Massenbr¨
uche von CO2werden
deutlich untersch¨
atzt. An der Außenseite f¨
ur r/R >0.8 werden die Profile gut reprodu-
ziert, w¨
ahrend sie an der Innenseite r/R <0.4 deutlich unter den experimentellen Daten
102 Kapitel 8. Simulation komplexer Str¨
omungen
0.0 0.3 0.6 0.9 1.2
r/R
0.00
0.02
0.04
0.06
0.08
0.0 0.3 0.6 0.9 1.2
r/R
0.00
0.02
0.04
0.06
0.08
0.10
mass fraction CO2
Exp. Sydney/Sandia
k-ε
LEA k-ε
EASM
x/D=0.9 x/D=1.3
Abbildung 8.8: Ermittelte Massenbr¨
uche von CO2und experimentelle Daten an den axia-
len Positionen x/D=0.9 und 1.3 f¨
ur Stufenk¨
orper-Diffusionsflamme
0.0 0.3 0.6 0.9 1.2
r/R
0.00
0.10
0.20
0.0 0.3 0.6 0.9 1.2
r/R
0.00
0.10
0.20
mass fraction O2
Exp. Sydney/Sandia
k-ε
LEA k-ε
EASM
x/D=0.9 x/D=1.3
Abbildung 8.9: Ermittelte Massenbr¨
uche von O2und experimentelle Daten an den axialen
Positionen x/D=0.9 und 1.3 f¨
ur Stufenk¨
orper-Diffusionsflamme
liegen.
Der Vergleich der Massenbruchprofile von O2mit den experimentellen Daten ist in Ab-
bildung 8.9 gezeigt. Auf der Mittellinie ist die Strahlstr¨
omung dominierend und es gibt
keinen Sauerstoff. Alle Modelle liefern dort eine richtige Verteilung. An der Außenseite
f¨
ur r/R >1.0 ist die Luftstr¨
omung dominierend und es wird ein Massenbruch wie in Luft
erzielt. Im Bereich von 0.3 <r/R <0.8 dagegen spielen die Turbulenzmodelle eine große
Rolle. Große Abweichungen k¨
onnen f¨
ur das k−²Modell beobachtet werden, w¨
ahrend
durch das LEA k−²Modell und das EASM gute ¨
Ubereinstimmungen mit den Messdaten
erzielt werden.
Abschließend erfolgt der Vergleich der Massenbruchprofile einer Komponente vom Brenn-
stoff H2mit den Messwerten in Abbildung 8.10. An allen Positionen sind die Massenbr¨
uche
H2an der Innenseite r/R = 0.3 ¨
ubersch¨
atzt. Bei x/D=0.9 ergeben alle Modelle sehr ¨
ahn-
liche Profile, w¨
ahrend bei x/D=1.3 ein Einfluss des Turbulenzmodells erkennbar ist. Das
8.2. Generischer KFZ-Außenspiegel 103
0.0 0.3 0.6 0.9 1.2
r/R
0.00
0.02
0.04
0.06
0.08
0.0 0.3 0.6 0.9 1.2
r/R
0.00
0.02
0.04
0.06
0.08
0.10
mass fraction H2
Exp. Sydney/Sandia
k-ε
LEA k-ε
EASM
x/D=0.9 x/D=1.3
Abbildung 8.10: Ermittelte Massenbr¨
uche von H2und experimentelle Daten an den axia-
len Positionen x/D=0.9 und 1.3 f¨
ur Stufenk¨
orperflamme
LEA k−²Modell und das EASM ergeben abermals ¨
ahnliche L¨
osungen, w¨
ahrend beim
k−²Modell die Abweichung zum Experiment gr¨
oßer ist.
Zusammenfassend kann festgestellt werden, dass das implementierte Flamelet–Modell zu-
friedenstellende Ergebnisse liefert, aber unter der Voraussetzung, dass das verwendete
Turbulenzmodell die Str¨
omung mit einer ausreichenden Qualit¨
at vorhersagt.
8.2 Generischer KFZ-Außenspiegel
Aufgrund der in den letzten Jahren erzielten Fortschritte bei der Reduzierung von Motor-
und Reifenger¨
ausch sind Windger¨
ausche im Fahrzeug heutzutage bereits ab 100 km/h
eine wesentliche L¨
armquelle. Dabei steht neben der Absenkung des Gesamtpegels vor
allem die Beeinflussung des subjektiven Ger¨
auschcharakters im Vordergrund. Aus dem
Spektrum herausragende Anteile, wie sie beispielsweise infolge von Str¨
omungsabl¨
osungen
bei der Umstr¨
omung des Außenspiegels entstehen, werden von den Fahrzeuginsassen als
besonders st¨
orend empfunden. F¨
ur die Vorhersage des aerodynamisch induzierten L¨
arms
wird vom Simulationsverfahren die Erfassung der Schallquelle und Schallausbreitung ver-
langt. Wegen des erheblichen numerischen Aufwands der gleichzeitigen Ber¨
ucksichtigung
von Schallentstehung und Schallausbreitung werden gegenw¨
artig oft entkoppelte Berech-
nungen durchgef¨
uhrt, in denen zuerst das Str¨
omungsfeld simuliert und daran anschließend
das Schallfeld mit Hilfe akustischer Analogien (Ffowcs Williams und Hawkings 1969) be-
rechnet wird. An das simulierte Str¨
omungsfeld wird dabei die Bedingung gestellt, die f¨
ur
die Schallentstehung wesentlichen Mechanismen, insbesondere die Druckschwankungen
wiederzugeben. An Hand erster Ergebnisse f¨
ur das Simulationsverfahren soll diese Eigen-
schaft gezeigt werden.
Abbildung 8.11 stellt die vereinfachte Modellgeometrie eines Außenspiegels. Messungen
hierzu wurden im Aeroakustik-Fahrzeugwindkanal am FKFS an der Universit¨
at Stuttgart
durchgef¨
uhrt. Eine Beschreibung dieser Messungen findet man in den Arbeiten von Sie-
104 Kapitel 8. Simulation komplexer Str¨
omungen
-0.5
0
0.5
1
X
-0.5
0
0.5
Y
9.64E+04 9.66E+04 9.69E+04 9.71E+04 9.73E+04 9.76E+04 9.78E+04 9.81E+04
Anströmungsrichtung
Platte
Modellspiegel 1.6
-0.8
-0.9
0.8
Abbildung 8.11: Prinzipskizze des Außenspiegelmodells, farblich ist die numerisch ermit-
telte Druckverteilung auf Platte und Spiegel dargestellt
gert u. a. (1999) und H¨
old u. a. (1999). Ziel des ausgew¨
ahlten Falls ist die Untersuchung
unterschiedlicher Turbulenzmodellellierungsans¨
atze im Hinblick auf eine nachfolgende Be-
rechnung des abgestrahlten L¨
arms. Dazu werden eine URANS mit dem SALSA Modell
und eine DES auf Basis des SALSA Modells untersucht. Alle Simulationen nutzten das
gleiche Gitter und die gleiche Zeitschrittweite, wodurch der Aufwand nahezu identisch ist
und der Einfluss des Ansatzes isoliert werden kann. Als Randbedingung wird der Ein-
str¨
omrand 1.2 m vor dem Beginn der Platte definiert, und die seitliche Begrenzung des
Gitters ragt auf jeder Seite um 0.7 m ¨
uber die Grundplatte hinaus. Der Ausflussrand
befindet sich 0.6 m hinter dem Plattenende. Die Unterseite ohne Grundplatte, die obere
Seite des Rechengebiets und die seitlichen R¨
ander werden als Symmetrieebenen definiert.
Das von ADtrans zur Verf¨
ugung gestellte Gitter wurde auch in den Simulation von Daim-
lerChrysler mit NSFLEX verwendet (H¨
old u. a. 1999). Es besteht aus ca. 1.2 ·106Punkten
und ist in 39 Bl¨
ocken unterteilt (siehe Abbildung 8.12), die dann von mehreren Prozes-
soren bearbeitet werden. Die Simulation werden auf 28 Knoten des institusteigenen PC
Clusters durchgef¨
uhrt. Die daraus resultierenden Oberfl¨
achendruckverteilungen dienen als
Basis f¨
ur die Berechnung des abgestrahlten Ger¨
ausches.
8.2. Generischer KFZ-Außenspiegel 105
(3D) 17 May 2001 WALL SURFACE
-0.5 0
X
Y
X
Z
(3D) 17 May 2001 WALL SURFACE
Abbildung 8.12: Ausschnitt des Gitters von DaimlerChrysler
Alle Simulationen werden mit einer Zeitschrittweite von 10−4sdurchgef¨
uhrt. Der Ein-
str¨
omrand wird mit einem konstanten Blockprofil mit u= 38.9 m/s (=140 km/h) und
νt= 1.44 ·10−4m2/s belegt.
Eine detaillierte Beurteilung des instation¨
aren Druckfelds erlauben Zeitschriebe an ein-
zelnen Messpunkten. In Abbildungen 8.14 und 8.15 sind die Simulationsergebnisse f¨
ur
die Positionen 121 und 123 (siehe Abbildung 8.13) aufgetragen. Deutlich sind die Un-
terschiede zwischen den Modellen zu erkennen. Das SALSA Modell ergibt ein wesentlich
regelm¨
aßigeres und schw¨
acheres Drucksignal, w¨
ahrend die DES ein stochastisches Signal
liefert, das dem Experiment eher gleicht. Im Hinblick auf eine nachfolgende Berechnung
der Schallabstrahlung des Spiegels sind die in den Drucksignalen enthalten Frequenzantei-
le von Interesse, weshalb eine Discrete Fourier Transform (DFT) Analyse der Zeitschriebe
vorzunehmen ist. An der Position 121 zeigt die DES deutlich gr¨
oßere Amplituden ¨
uber
den gesamten Frequenzbereich, Abbildung 8.16. Die ¨
Ubereinstimmung f¨
ur Frequenzen un-
terhalb 100 Hz ist gut, danach fallen alle Ergebnisse leicht ab. Die RANS ergibt wesentlich
kleinere Amplituden, die ¨
uber den gesamten Frequenzbereich nicht mit dem Experiment
¨
ubereinstimmen. An der Position 123 folgen beide Modelle dem Experiment ¨
uber den
gesamten Frequenzbereich. Bei 2 kHz ist die gr¨
oßte Abweichung ca. 6 dB.
Die beiden Messpunkte 121 und 123 repr¨
asentieren unterschiedliche Str¨
omungszust¨
ande.
106 Kapitel 8. Simulation komplexer Str¨
omungen
b121 123
b
R=0.1 m 120 122
119
6
Körper
b =0.1 m
Anströmung
Abbildung 8.13: Messpositionen auf der Grundplatte
W¨
ahrend der Punkt 121 im Abl¨
osungsgebiet liegt, befindet sich der Punkt 123 außerhalb
oder am Rand des Abl¨
osungsgebiets. Dies begr¨
undet das unterschiedliche Verhalten der
beiden Ans¨
atze zur Behandlung der Turbulenz an den unterschiedlichen Punkten, denn die
DES entfaltet ihre St¨
arken erst in Gebieten massiver Abl¨
osung, so dass die Vorhersagen
der DES im Punkt 121 deutlich n¨
aher am Experiment sind. Das kann den Abbildungen
8.16 und 8.17 entnommen werden. Da der Rechenaufwand auf identischen Gittern nahezu
unver¨
andert bleibt, besitzt die DES somit ein hohes Anwendungspotential zur Simulation
der Schallquellen in komplexen Geometrien.
Anzumerken ist, dass die DFT Analyse der Drucksignale sehr stark vom Simulationsver-
lauf anh¨
angt. Durch die im Experiment beobachtete Variation des dominanten Hiebtones
im Bereich von 20 - 40 Hz wird dem Signal eine Schwebung ¨
uberlagert, die eine lange Simu-
lationsdauer erforderlich macht. Die DFT Analyse selbst reagiert stark auf unvollst¨
andige
Perioden im Zeitschrieb. Bei der ¨
Uberlagerung verschiedener Frequenzen im Drucksignal
kann aber nicht ausgeschlossen werden, dass Perioden einzelner Anteile unvollst¨
andig
sind. Daher kann durch eine Auswahl der Eingangsl¨
ange des Signals ein g¨
anzlich anderes
Ergebnis erzielt werden. W¨
ahrend der Simulation wurde fortlaufend das Druckspektrum
betrachtet und festgestellt, dass es innerhalb der Simulation stark variiert.
8.2. Generischer KFZ-Außenspiegel 107
0 0.02 0.04 0.06 0.08 0.1
t[s]
-1200
-1000
-800
-600
-400
-200
0
200
400
600
800
1000
Ampl. P[Pa]
SALSA
DES
Abbildung 8.14: Druckverlauf an der Position 121
0 0.02 0.04 0.06 0.08 0.1
t[s]
-600
-500
-400
-300
-200
-100
0
100
200
300
400
Ampl. P[Pa]
SALSA
DES
Abbildung 8.15: Druckverlauf an der Position 123
108 Kapitel 8. Simulation komplexer Str¨
omungen
101102103
f [Hz]
10-2
10-1
100
101
102
Ampl. P. [Pa]
Messung
SALSA
DES
Abbildung 8.16: Vergleich der Spektren zwischen Rechnung und Messung im Punkt 121
101102103
f [Hz]
60
70
80
90
100
110
120
130
140
150
Ampl. P [DB]
Messung
SALSA
DES
Abbildung 8.17: Vergleich der Spektren zwischen Rechnung und Messung im Punkt 123
Kapitel 9
Zusammenfassung
Diese Arbeit befaßt sich mit der Entwicklung der biharmonischen Gittergenerierung mit
Beschleunigung durch Mehrgitterverfahren, sowie der Entwicklung eines Mehrgitterver-
fahrens zur L¨
osung der Navier-Stokes Gleichungen.
Bei der Gittergenerierung wird in der vorliegenden Arbeit die biharmonische Gleichung
zur Gl¨
attung verwendet, wobei die Kontrollfunktionen mittels der zweiten Randbedingung
der Differentialgleichungen vierter Ordnung bestimmt werden, die durch die Einf¨
uhrung
einer virtuellen Gitterschicht iterativ bestimmt werden. Zur Validierung werden einige
zweidimensionale bzw. dreidimensionale praxisrelevante Rechengitter gew¨
ahlt, an denen
gezeigt wird, dass die Gitterqualit¨
at durch das entwickelte Verfahren im Vergleich zu
herk¨
ommlichen Methoden deutlich verbessert wird. Das Verfahren besitzt folglich großes
Anwendungspotential in der strukturierten Gittergenerierung. Die Stabilit¨
at und Konver-
genzgeschwindigkeit werden durch die effektive Behandlung der Randbedingungen und
das Mehrgitterverfahren verbessert. Das Verhalten des Mehrgitterverfahrens auf Basis
von CS/FAS-Verfahren ist jeweils mit dem V- bzw. W-Zyklus untersucht worden.
Zur Beschleunigung der L¨
osung der Navier-Stokes Gleichungen wird das Mehrgitterverfah-
ren in Kombination mit dem SIMPLE-Algorithmus verwendet. Dazu wird ein FMG-FAS
Verfahren benutzt, bei dem die Startl¨
osung inklusive Massenfluss des Grobgitters vom
feinen Gitter restringiert wird. Dies f¨
uhrt zu dem Problem, dass Geschwindigkeiten und
Massenfl¨
usse nicht zueinander passen. Bei der Simulation turbulenter Str¨
omungen wer-
den zus¨
atzlich die Turbulenzvariablen vom feinen Gitter restringiert, wodurch ebenfalls
die Beziehungen zwischen der Wirbelviskosit¨
at und der turbulenten Energie bzw. ihrer
Dissipationsrate nicht mehr richtig erf¨
ullt werden. In dieser Arbeit wurde daf¨
ur ein modi-
fiziertes FMG Verfahren mit V-Zyklus entwickelt, bei dem die Startl¨
osung inklusive der
Turbulenzvariablen nicht vom feinen Gitter ¨
ubertragen, sondern vom vorherigen Zyklus
¨
ubernommen wird. Der modifizierte Algorithmus hat gegen¨
uber dem Standard FMG-FAS
den Vorteil, dass er numerisch stabil ist und einfach implementiert werden kann. Er wur-
de anhand einer Reihe von zwei- und dreidimensionalen Testf¨
allen validiert, wobei f¨
ur
laminare Str¨
omungen Effizienzsteigerungen des Gesamtverfahrens um Faktoren von 6 bis
180 erreichbar sind. Besonders f¨
ur feine Gitter ist der entwickelte Algorithmus deutlich
effizienter als das Standard FMG-FAS. F¨
ur turbulente Str¨
omungen werden ebenfalls sehr
109
110 Kapitel 9. Zusammenfassung
bemerkenswerte Beschleunigungsfaktoren von 6 bis 36 erreicht.
Als Anwendungsf¨
alle werden eine Stufenk¨
orper-Diffusionsflamme und die instation¨
are
Umstr¨
omung eines generischen KFZ-Außenspiegels simuliert, die neben den schon ge-
zeigten Einfl¨
ussen von Numerik und Gitteraufl¨
osung noch den Einfluss der Turbulenzmo-
dellierung dokumentiert.
Literaturverzeichnis
[1] Almeida, G.P., Durao, D.F.G. und Heitor, M.V. (1993). Wakes flow behind two-
dimensional model hills. Experimental Thermal and Fluid Science,7, 87-101.
[2] Bai, X.S. und Fucks, L. (1992). Fast multigrid method for 3-D turbulent incom-
pressible flows. Int. J. Num. Meth. Heat Fluid Flow,2, 127-137.
[3] Baldwin, B.S. und Lomax, H. (1978). Thin layer approximation and algebraic
model for separated turbulent flows. AIAA paper 78-0257.
[4] Bell, J.B., Shubin, G.R. und Stephen, A.B. (1982). A segmentation approach to
grid generation using Biharmonics. J. Comp. Phys. 47, 463-472.
[5] Bilger, R.W. (1993). Conditional moment closure for turbulent reacting flows.
Phys. Fluids A5, 436-444.
[6] Boussinesq, J. (1877). Th´eorie de l’´
Ecoulement tourbillant. Pr´sent´es par Divers
Savantes Acad. Sci. Inst. Fr. 23, 46-50.
[7] Brandt, A. und Dym, J. (1994). Effective boundary treatment for the biharmo-
nic dirichlet problem. Proc. Seventh Copper Mountain Conference on Multigrid
Methods.
[8] Chien, K.Y. (1982). Predictions of channel and boundary-layer flows with a Low-
Reynolds-Number turbulence model. AIAA J. 20, 33-38.
[9] Demirdˇzi´c, I., Lilek, ˇ
Z. und Peri´c, M. (1992). Fluid flow and heat transfer tests
problems for non-orthogonal grids: bench-mark solutions. Int. J. Numer. Meth.
Fluids,15, 329-354.
[10] Demirdˇzi´c, I., Lilek, ˇ
Z. und Peri´c, M. (1993). A colloacted finite volume method
for predicting flows at all speeds. Int. J. Numer. Meth. Fluids,16, 1029-1050.
[11] Demirdˇzi´c, I. und Peri´c, M. (1996). Numerical study of grey-body surface radiation
coupled with fluid flow for general geometries using a finite volume multigrid solver.
Int. J. Heat Fluid Flow,6, 3-18.
[12] Edwards, J.R. (1996). Comparison of eddy viscosity-transport turbulence models
for three dimensional, shock-separated flow fields. AIAA J. 34, 756-763.
111
112 LITERATURVERZEICHNIS
[13] Favr´e, A. (1965). Equations des gaz turbulents compressibles. Jounal de
Mech´eanique,4, 361-390.
[14] Findling, A. und Hermann, U. (1991). Development of an efficient and robust
solver for elliptic grid Generation, In Numerical grid generation in computational
fields simulation and related fields, (Eds. Arcilla, A.S., H¨
auser, J., Eisemann P.R.
and Thompson J.F.), pp. 781-792 , Proceedings of the 3rd International Grid
Conference, Barcelona, June.
[15] Gibson M.M. und Younis, B.A. (1986). Calculation of swirling jets with a Reynolds
stress closure. Phys. Fluids,29, 38-48.
[16] Hackbusch, W. (1985). Multigrid methods and applications. Springer–Verlag New
York Heidelberg Berlin.
[17] Harten, A. (1983). High resolution schemes for hyperlic conservation laws. J.
Comp. Phys. 49, 357-393.
[18] H¨
old, R., Brenneis, A., Eberle, A., Schwarz, V. und Siegert, R. (1999). Numerical
simulation of aeroacoustic sound generated by generic bodies placed on a plate:
Part I - Prediction of aeroacoustic sources. AIAA paper 99-1896.
[19] Huurdeman, R. (1999). Numerische Simulation inkompressibler turbulenter
Str¨
omungen mit Mehrgitterverfahren auf unstrukturierten Gittern. Dissertation,
Universit¨
at Stuttgart, Germany.
[20] Jain, R.K. (1986). Generation of body-fitted grids around airfoils using multi-
grid method. In Numerical grid generation in computaional fluid dynamics (Eds.
H¨
auser, J. and Taylor, C.).
[21] Jischa, M. (1982). Konvektiver Impuls-, W¨
arme- und Stoffaustausch. In Grundla-
gen der Ingenieurwissenschaften, Vieweg, Braunschweig/Wiesbaden.
[22] Jones, W.P. und Launder, B.E. (1972). Prediction of laminarization with a two-
equation turbulence model. Int. J. Heat and Mass Transfer,15, 301-314.
[23] Lai, Y.G. (1995). Computational method of second-mpment turbulence closures
in complex geometries. AIAA J. 33, 1425-1432.
[24] Lam, C.K.G. und Bremhorst, K. (1981). A modified form of the k−²model for
predicting wall turbulence. J. Fluids Engineering,103, 456-460.
[25] Launder, B.E. und Spalding, D.B. (1974). The numerical computation of turbulent
flows. Comput. Meth. Appl. Mech. Engrg. 3, 29-289.
[26] Launder, B.E. und Sharma, B.I. (1974). Application of the energy dissipation
model of turbulence to the calculation of flow near a spinning disc. Letters in Heat
and Mass Transfer,1, 131-138.
LITERATURVERZEICHNIS 113
[27] Launder, B.E., Reece, G. und Rodi, W. (1975). Progress in the development of a
Reynolds stress turbulence closure. J. Fluid Mech. 68, 537-566.
[28] Launder, B.E. und Shima, N. (1989). Second-moment closure for the near-wall
sublayer: Development and application. AIAA J. 27, 1319-1325.
[29] Leister, H.J. und Peri´c, M. (1992). Numerical simulation of a 3D Czochralski melt
flow by a finite volume multigrid algorithm. J. Crystal Growth,123, 567-574.
[30] Leonard, B.P. (1979). A stable and accurate convective modelling procedure based
on quadratic upstream interpolation. Comput. Meth. Appl. Eng. 19, 59-98.
[31] Lien, F.S. und Leschziner, M.A. (1991). Multigrid convergence acceleration for
complex flow including turbulence. Int. Series of Numerical Mathematics,98, 277-
288.
[32] Lien, F.S. und Leschziner, M.A. (1993). Computational modelling of 3D turbulent
flow in S-diffuser and transition ducts. In Engineering Turbulence Modelling and
Experiments 2, pp. 217-228, Elsevier, Amsterdam.
[33] Lien, F.S. und Leschziner, M.A. (1994a). Multigrid acceleration for recirculating
laminar and turbulent flows computed with a non-orthogonal col-located finite-
volume scheme. Comput. Meth. Appl. Mech. Engrg. 118, 351-371.
[34] Lien, F.S. und Leschziner, M.A. (1994b). Assessment of turbulence-transport mo-
dels including non-linear RNG eddy-viscosity formulation and second-moment
closure for flow over a backward-facing step. Compt. Fluids,23, 983-1004.
[35] Lien, F.S. und Leschziner, M.A. (1994c). A general non-orthogonal collocated finite
volume algorithm for turbulent flow at all speeds incorporating second-moment
turbulent-transport closure. Comput. Meth. Appl. Mech. Engrg. 114, 123-167.
[36] Lilek, ˇ
Z. und Peri´c, M. (1995). A fourth-order finite-volume method with collocated
variable arrangement. Comput. Fluids,24, 239-252.
[37] Lilek, ˇ
Z. Muzaferija, S. und Peri´c, M. (1997). Efficiency and accuracy aspects of
a full-multigrid SIMPLE algorithmus for three-dimensional flows. Numerical Heat
Transfer, Part B,31, 23-42.
[38] Linden, J. (1984). A multigrid method for solving the biharmonic equation on
rectangular domains. In Advanced in Multigrid Methods, pp. 64-76, Proceedings of
the Conference held in Oberwolfach, Dec. 8-13. F. Vieweg & Sohn, Braunschweig,
Wiesbaden.
[39] Mavripilis, D.J. (1995). Three-dimensional multigrid Reynold-averaged Navier-
Stokes solver for unstructured meshes, AIAA J. 33, 445-453.
[40] Obi, S., Peri´c, M. und Scheuerer, G. (1991). Second-moment calculation procedure
for turbulent flows with collocated variable arrangement. AIAA J. 29, 585-590.
114 LITERATURVERZEICHNIS
[41] Orth, A. (1991). Mehrgittermethod zur Berechnung inkompressibler, station¨
arer
Str¨
omungen mit krummlinigen Berandungen. Dissertation, Karlsruhe, Germany.
[42] Patankar, S.V. (1980). Numerical heat transfer and fluid flow. Hemisphere, Wa-
shington D. C.
[43] Patel, V.C., Rodi, W. und Scheuerer, G. (1985). Turbulence models for near-wall
and low Reynolds numer flows: a review. AIAA J. 23, 1308-1319.
[44] Peters, N. (2000). Turbulent Combustion. Combridge University Press.
[45] Pitsch, H., Riesmeier, E. und Peters, N. (2000). Unsteady flamelet modelling of
soot formation in turbulent diffusion flames. Combust. Sci. and Tech. 158, 389-
406.
[46] Piomelli, U. und Chasnov, J.R. (1996). Large-eddy simulation, theory and app-
lications. In Turbulence and Transition Modelling (Eds. Hallb¨
ack, M. et al.), pp.
269-336. Kluwer Academic Publishers. Dordrecht Boston London.
[47] Pope, S.B. (1981). A Monte Carlo method for the PDF equations of turbulent
reactive flow’ Combust. Sci. and Tech. 25, 159-174.
[48] Pope, S.B. (1985). PDF Methods for turbulent reactive flows. Prog. Energy Com-
bust. Sci. 11, 119-192.
[49] Pope, S.B. (1990). Computations of turbulent combustion: Progress and Challen-
ges. In Twenty-Third Symposium (International) on Combustion, pp. 591-612. The
Combustion Institue, Pittsburgh.
[50] Rhie, C.M. und Chow, W.L. (1983). Numerical Study of the Turbulent Flow Past
an Airfoil with Trailing Edge Separation. AIAA J. 21, 1325-1332.
[51] Rodi, W. (1991). Two-layer models of turbulence. AIAA Paper 91-0216.
[52] Rodi, W. und Mansour, N.N. (1993). Low-Reynolds number k−²modelling with
the aid of direct numerical simulation data. J. Fluid Mech. 250, 509-524.
[53] Rung, T., Thiele, F. und Fu, S. (1999). On the realizabililty of nonlinear stree-
strain relationships for Reynolds stress closures. Flow, Turbulence and Combusti-
on,60, 333-359.
[54] Rung, T. (1999). Formulierung universeller wandrandbedingung f¨
ur transport-
gleichungsturbulenzmodelle. Institusbericht, 2-99, Hermann-F¨
ottinger-Institut f¨
ur
Str¨
omungsmechanik, Technische Universit¨
at Berlin, Germany.
[55] Rung, T. (2000). Entwicklung anisotroper Wirbelz¨
ahigkeitsbezeiehung mit Hilfe
von Projektionstechniken. Dissertation, Technische Universit¨
at Berlin, Germany.
[56] Rung, T., Bunge, U., Schatz, M. und Thiele, F. (2003). Restatement of the Spalart-
Allmaras Eddy-Viscosity Model in a Strain-Adaptive Formulation. To be published
in AIAA J.
LITERATURVERZEICHNIS 115
[57] Rousseau, A.N., Albright, L.D. und Torrance, K.E. (1997). A short comparison
of damping functions of standard low-Reynolds-number k−²models. J. Fluids
Engineering,119, 460-463.
[58] Shur, M., Spalart, P.R., Strelets, M. und Travin, A. (1999). Detached-Eddy Simu-
lation of an airfoil at high angle of attack. In Turbulent Shear Flows (Eds. Rodi,
W. and Laurence, D.) pp. 669-678. Elsevier Science Ltd.
[59] Siegert, R., Schwarz, V. und Reichenberger, J. (1999). Numerical simulation of
aeroacoustic sound generated by generic bodies placed on a plate: Part II – Pre-
diction of radiated sound pressure. AIAA paper 99-1895.
[60] Steger, J.L. und Sorenson, I.N. (1979). Automatic mesh point clustering near a
boundary in grid generation with elliptic partial differential equations. J. Comp.
Phys. 33, 405-410.
[61] Spalart, P.R. und Allmaras, S.R. (1992). A one-equation turbulence model for
aerodynamic flows. AIAA Paper 92-439.
[62] Spalart, P.R. (1999). Strategies for turbulence modelling and simulation. Enginee-
ring Turbulence Modelling and Experiments 4, 3-17, Elsevier, Amsterdam.
[63] Spalart, P.R. (2001). Young-person’s guid to Detached-Eddy Simulation grids.
NASA/CR-2001-211032.
[64] Sparis, P.D. (1985). A method for generating boundary-orthogonal curviliear coor-
nidate systems using the biharmonic equation. J. Comp. Phys. 61, 445-462.
[65] Sparis, P.D. und Karkanis, A. (1992). Bounday-orthogonal biharmonic grids via
preconditioned gradient methods. AIAA J. 30, 671-678.
[66] Spitaleri, R.M. (1990). A multigrid method for grid generation with line-spacing
control. Applied Mathematics and Computation,40, 125-134.
[67] Spitaleri, R.M. (1991). Multigrid grid generation in scalar and parallel environ-
ments. In Numerical Grid Generation in Computational Fluid Dynamics and Re-
lated Fields (ed. Hauser, A.S, Eisemann P.R. und Thompson, J.F.), pp. 805- .
Proceedings of the 3rd International Grid Conference, Barcelona, June.
[68] Stone, H.L. (1968). Iterative solution of implicit approximations of multidimen-
sional partial differential equations. SIAM J. Numer. Anal. 5, 530-558.
[69] Strelets, M., Travin, A., Shur, M. und Spalart, P.R. (2000). Detched-eddy simula-
tions past a circular cylinder. Flow, Turbulence and Combustion,63, 293-313.
[70] St¨
uben, K. und Linden, J. (1986). Multigrid methods: an overview with emphasis
on grid generation processes. In Numerical Grid Generation in Computational
Fluid Dynamics (Eds. H¨
auser, J. and Taylor, C.), pp. 483-509. Proceeding of the
International Conference, Landshut, West Germany, July 14-19.
116 LITERATURVERZEICHNIS
[71] Shubin, G.R., Stephen, A.B. und Bell, J.B. (1982). Three dimensional grid ge-
neration using biharmonic. In Numerical Grid Generation, pp. 761-774. Elsevier
Sc.
[72] Sutherland, D.M. (1893). The viscosity of gases and molecular forces. Phil. Mag.
Ser. 36(5), 507-530.
[73] Taylor, A.M., Whitelaw J.H. und Yianneskis, M. (1982). Curved Ducts With
Strong Secondary Motion: Velocity Measurements of Developing Laminar and Tur-
bulent Flow. J. Fluids Engineering,104, 350-359.
[74] Thomas, P.D. und Middlecoef, J.F. (1980). Direct control of the grid point distri-
bution in meshes generated by elliptic equations. AIAA J. 18, 652-656.
[75] Thompson, J.F., Warsi, Z.U.A. und Mastin, C.W. (1985). Numerical grid gene-
ration: foundations and applications. North-Holland, New York ·Amsterdam ·
Oxford.
[76] Thompson, J.F. (1987). A General Three-Dimensional Elliptic Grid Generation
System on A Composite Block Structure. Comput. Meth. Appl. Mech. Engrg. 64,
377-411.
[77] Thompson, J.F. (1996). A reflection on grid generation in the 90s: trends, needs,
and influences. Proc. of 5th Int. Grid Conference: Numerical Grid Generation in
Computational Fluid Dynamics, pp. 1029-1110.
[78] Vanka, S.P. (1986a). A calculation procedure for three-dimensional steady recir-
culating flows using multigrid methods. Comput. Meth. Appl. Mech. Engrg. 55,
321-338.
[79] Vanka, S.P. (1986b). Block-implicit multigrid solution of Navier-Stokes equations
in primitive variables. J. Comput. Phys. 65, 138-158.
[80] Warnatz, J. und Mass, U. (1993). Technische Verbrennung. Springer-Verlag, 1993.
[81] Weatherill, N.P. (1990). Numerical grid generation. Lecture series 1990-06, Von
Karman institut for fluid dynamics.
[82] White, J.A. (1990). Elliptic Grid Generation with Orthogonality and Spacing Con-
trol On An Arbitrary Number of Boundaries. AIAA paper 90-1568
[83] Ffowcs Williams, J.E. und Hawkings, D.L. (1969). Sound generated by turbulence
and surfaces in arbitrary motion. Philosophical Transaction of the Royal Societym.
264, 321-342.
[84] Wesseling, P. (1992). An introduction to multigrid methods. John Wiley & Sons,
Chichester ·New York ·Brisbane ·Toronto ·Singapore.
[85] Wilcox, D.C. (1988). Reassessment of the scale determining equation for advanced
turbulence models. AIAA J. 26, 1299-1310.
LITERATURVERZEICHNIS 117
[86] Wilcox, D.C. (1993). Turbulence modelling for CFD. DCW Industries, Inc., La
Canada, CA.
[87] Xue, L. (1998). Entwicklung eines effizienten parallelen L¨
osungsalgorithms zur
dreidimensionalen Simulation komplexer turbulenter Str¨
omungen. Dissertation,
Technische Universit¨
at Berlin, Germany.
[88] Yang, Z. und Shih, T.H. (1993). New time scale based k−²model for near-wall
turbulence. AIAA J. 31, 1191-1198.
[89] Yan, J., Xue, L. und Thiele, F. (1998). Performance of multigrid method in bi-
harmonic grid generation. In Numerical Grid Generation in Computational Fields
Simulation and Related Fields (ed. Cross, M., Soni, B.K., Thompson, J.F., Hau-
ser, J. und Eiseman, P.R.), pp. 303-312. Proceedings of the 6th International Grid
Conference, University of Greenwich, July.
[90] Yan, J. und Thiele, F. (1998). Performance and accuracy of a modified full multi-
grid algorithm for fluid flow and heat transfer. Numerical Heat Transfer, Part B,
34, 323-338.
[91] Yan, J. (1999). A modified full multigrid algorithm for the Navier-Stokes equations.
Institusbericht, Hermann-F¨
ottinger-Institut f¨
ur Str¨
omungsmechanik, Technische
Universit¨
at Berlin, Germany.
[92] Zheng, X., Liao, C., Liu, C., Sung C.H. und Huang, T.T. (1997a). Multigrid com-
putation of incompressible flows using two-equation turbulence models: part I –
numerical method. J. Fluids Enginnering,119, 893-899.
[93] Zheng, X., Liao, C., Liu, C., Sung, C.H. und Huang, T.T. (1997b). Multigrid
computation of incompressible flows using two-equation turbulence models: part
II – applications. J. Fluids Enginnering,119, 900-905.