Verbesserung Der Numerischen Simulation Der Mischung Von Triebwerksstrahlen

Universität der Bundeswerhr München

Fakultät für Luft- und Raumfahrttechnik

Institut für Thermodynamik
Verbesserung der numerischen Simulation

der Mischung von Triebwerksstrahlen
Sebastian F. Saegeler
Vollständiger Abdruck der von der Fakultät für Luft- und Raumfahrttechnik der Universität
der Bundeswehr München zur Erlangung des akademischen Grades eines
Doktor-Ingenieures (Dr.-Ing.)
genehmigten Dissertation
Vorsitzender: Prof. Dr.-Ing. R. Förstner

1. Berichterstatter Prof. Dr.-Ing. C. Mundt
2. Berichterstatter Prof. Dr.-Ing. S. Staudacher
Diese Dissertation wurde am 8.10.2013 bei der Universität der Bundeswehr München eingereicht und
durch die Fakultät für Luft- und Raumfahrttechnik am 16.10.2013 angenommen.
Tag der mündlichen Prüfung war der 28.03.2014.
Danksagung
Die vorliegende Arbeit entstand während meiner Tätigkeit als wissenschaftlicher Mitar-
beiter am Institut für Thermodynamik der Universität der Bundeswehr München.
Ich möchte mich ganz herzlich bei Herrn Professor Dr.-Ing. Christian Mundt für die
persönliche Betreuung während dieser Zeit bedanken. Die intensive fachliche und per-
sönliche Unterstützung sowie die zahlreichen anregenden Diskussionen während dieser
Zeit waren sehr hilf- und lehrreich. Dadurch hatte ich stets die dafür nötige Motivation
um an der Fertigstellung dieser Dissertation festzuhalten.
Ebenso möchte ich mich bei Herrn Professor Dr.-Ing. Stephan Staudacher für die Über-
nahme der Zweitprüferschaft aufs herzlichste bedanken.
Herrn Professor Dr.-Ing. Roger Förstner danke ich ganz herzlich für die Übernahme des
Prüfungsvorsitzes.
Gefördert wurde die Arbeit durch das Luftfahrtforschungsprogramm der Bundesrepu-
blik Deutschland sowie durch Rolls-Royce Deutschland Ltd. & Co. KG im Rahmen des
gemeinsamen Projekts OPTITHECK. Vor allem die Zusammenarbeit und der wissen-
schaftliche Austausch mit Rolls-Royce hat im wesentlichen Maße zur Verwirklichung
dieser Arbeit beigetragen. Besonderer Dank gilt Dr.-Ing. Jan Lieser von Rolls-Royce für
dessen Unterstützung.
Ferner bedanke ich mich noch bei meinen Kolleginen und Kollegen für zahlreiche inspi-
rierende und anregende Diskussionen sowie dem involvierten technischen Personal, vor
allem des Rechenzentrums der UniBW.
Eine solche Arbeit wäre ohne den Rückhalt und Antrieb von Familie und Freunden
kaum erfolgreich gewesen. Deshalb Danke ich den Personen, die mir während der Arbeit
an meiner Dissertation nahe standen und mich unterstützt haben. Eine ganz besonders
herzliche Danksagung gilt meiner Freundin Anne. Ihr uneingeschränkter Beistand und
Rückhalt diente mir als größtmögliche Motivation und Antrieb.
ii
Inhaltsverzeichnis
1. Einleitung 1
1.1. Einführung . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1
1.2. Stand der Forschung . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4
1.3. Zielsetzung der Arbeit . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
2. Prinzip und Leistungsbewertung von Schubdüsen mit internem Mischer 9

2.1. Grundprinzip interner Strahlmischung . . . . . . . . . . . . . . . . . . . 9
2.2. Quantitative Leistungsbewertung . . . . . . . . . . . . . . . . . . . . . . 11
2.2.1. Durchflusskoeffizient . . . . . . . . . . . . . . . . . . . . . . . . . 11
2.2.2. Geschwindigkeitskoeffizient . . . . . . . . . . . . . . . . . . . . . . 12
3. Numerische Verfahren 15
3.1. Navier-Stokes Gleichungen . . . . . . . . . . . . . . . . . . . . . . . . . . 15
3.2. Gemittelte Navier-Stokes Gleichungen . . . . . . . . . . . . . . . . . . . . 17
3.3. Turbulenzmodellierung . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19
3.3.1. Allgemeine Schließungsansätze . . . . . . . . . . . . . . . . . . . . 19
3.3.2. Turbulenzmodelle . . . . . . . . . . . . . . . . . . . . . . . . . . . 24
3.4. Strömungsberechnung und Lösungsverfahren . . . . . . . . . . . . . . . . 27
3.4.1. Numerische Flussberechnung . . . . . . . . . . . . . . . . . . . . . 28
3.4.2. Preconditioning . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28
3.5. Randbedingungen . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31
3.5.1. Einlass . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31
3.5.2. Auslass . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33
3.5.3. Fernfeld . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33
3.5.4. Wände . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33
3.5.5. Seitenränder . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34
3.6. Zeitschrittverfahren . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35
3.6.1. Lokales Zeitschrittverfahren . . . . . . . . . . . . . . . . . . . . . 35
3.6.2. Runge-Kutta Zeitschrittverfahren . . . . . . . . . . . . . . . . . . 36
3.7. Gittergenerierung . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37
3.7.1. Geometrieerzeugung . . . . . . . . . . . . . . . . . . . . . . . . . 37
3.7.2. Vernetzung . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38
iv
3.7.3. Fertigstellung . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40
4. Simulation der Düsenströmung mit OpenFOAM 42

4.1. Simulationsgrundlagen . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42
4.1.1. Berechnungsdetails . . . . . . . . . . . . . . . . . . . . . . . . . . 42
4.1.2. Randbedingungen . . . . . . . . . . . . . . . . . . . . . . . . . . . 44
4.1.3. Simulationsablauf . . . . . . . . . . . . . . . . . . . . . . . . . . . 45
4.1.4. Konvergenzkriterien . . . . . . . . . . . . . . . . . . . . . . . . . . 46
4.1.5. Postprocessing . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46
4.2. Simulation der Düsenströmung mit Temperaturkorrekturmodellen . . . . 47
4.2.1. Prinzip der Temperaturkorrektur . . . . . . . . . . . . . . . . . . 47
4.2.2. Ergebnisse . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50
4.2.3. Zusammenfassung - Temperaturkorrektur . . . . . . . . . . . . . . 61
4.3. Simulation der Düsenströmung mit verbessertem Wärmestromterm . . . 63
4.3.1. Verbesserter turbulenter Wärmestromschließungsansatz . . . . . . 63
4.3.2. Ergebnisse . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68
4.3.3. Zusammenfassung - Verbessertes turbulentes Wärmestrommodell 76
4.4. Koppelung von Temperaturkorrektur und verbessertem Wärmestromterm 78
4.4.1. Experimentelle Messungen . . . . . . . . . . . . . . . . . . . . . . 78
4.4.2. Rechenergebnisse mit existierenden Modellen . . . . . . . . . . . . 81
4.4.3. Koppelungsmodell . . . . . . . . . . . . . . . . . . . . . . . . . . 85
4.4.4. Ergebnisse . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 87
4.4.5. Zusammenfassung - gekoppeltes Mischungsmodell . . . . . . . . . 96
5. Parameteruntersuchung der Düsenströmung mit ANSYS FLUENT 98

5.1. Rechenmodell und Randbedingungen . . . . . . . . . . . . . . . . . . . . 98
5.2. Ergebnisse Parameteruntersuchung . . . . . . . . . . . . . . . . . . . . . 101
5.2.1. Gittergröße . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101
5.2.2. Periodische Ränder/Mischertaktung . . . . . . . . . . . . . . . . . 103
5.2.3. Konstante/variable spezifische Wärmekapazität . . . . . . . . . . 105
5.2.4. Turbulenzmodell und dessen Diskretisierung . . . . . . . . . . . . 107
5.3. Zusammenfassung - Parameteruntersuchung . . . . . . . . . . . . . . . . 111
6. Zusammenfassung 112
Literaturverzeichnis 116
A. Kompressibel-turbulenter skalarer Transport 126
B. Modifizierte ω-Gleichung für Temperaturkorrekturmodell 130
v
Nomenklatur
Lateinische Buchstaben
A Fläche
AF an Querschnittfläche am Eintritt des Fanstroms
ACore Querschnittfläche am Eintritt des Kernstroms
AW and Wandoberfläche
Aexit Düsenquerschnittsfläche am Auslass
Athroat Düsenhalsfläche
a Schallgeschwinigkeit
Cjθ konvektiver Transport von u′′j θ
CD Durchflusskoeffizient
CN Konditionszahl
CT Temperaturkorrekturfaktor
CT g Modellkonstante für Temperaturkorrekturfaktor
CV Geschwindigkeitskoeffizient
Csθ Modellkonstante für differentielles Wärmetransportmodell
Cµ Modellkonstante für Turbulenzmodelle
CF L Courant-Friedrichs-Lewy Zahl
cp spezifische Wärmekapazität bei konsantem Druck
cv spezifische Wärmekapazität bei konstantem Volumen
c1θ Modellkonstante für differentielles Wärmetransportmodell
c2θ Modellkonstante für differentielles Wärmetransportmodell
t
Djθ turbulent diffusiver Transport von u′′j θ
v
Djθ viskos diffusiver Transport von u′′j θ
E Totalenergie
vi
e innere Energie
F~ Vektor konvektiver Flüsse
FB,x Bruttoschub in x-Richtung
f~ext Vektor externer Körperkräfte
f~n Ortsvektor zum Mittelpunkt einer Zellfläche
G~ Vektor diffusiver Flüsse
H Totalenthalpie
H(x) Heaviside Sprungfunktion
h Enthalpie
I turbulente Intensität
k turbulente kinetische Energie
l turbulentes Längenmaß
Ma Machzahl
M aT turbulente Machzahl
ṁ Massenstrom
ṁid idealer Massenstrom
~n Auswärts gerichteter normierter Oberflächennormalenvektor
N Anzahl an Zellen im Rechengebiet
P gemittelter statischer Druck
Pr laminare Prandtlzahl
P rT turbulente Prandtlzahl
P~ Ortsvektor zum Mittelpunkt einer Zelle
Pjθ,1 Produktion von u′′j θ infolge von Temperaturgradienten
Pjθ,2 Produktion von u′′j θ infolge von Scherung
p statischer Druck
pt Totaldruck
p∗ kritischer Druck
p∞ Umgebungsdruck
Q reduzierter Massenstrom/Q-Funktion
Qideal idealer reduzierter Massenstrom/Q-Funktion
Q̇ Wärmestrom
~
Q Vektor primitiver Variablen
qj Wärmestromvektor in Indexnotation
vii
qLj laminarer Wärmestromvektor in Indexnotation
qT j turbulenter Wärmestromvektor in Indexnotation
~q Wärmestromvektorvektor
R spezifische Gaskonstante
Re Reynoldszahl
~
R Residuenvektor
S Sutherlandkonstante
Sij Deformationsgeschwindigkeit
T statische Temperatur
Tg normierter Totaltemperaturgradient
T gExp Modellkonstante für normierten Totaltemperaturgradient
Tt Totaltemperatur
t Zeit
tij gemittelter viskoser Stresstensor
∆t Zeitschrittweite
∆tinv reibungsfreie Zeitschrittweite
∆tvis viskose Zeitschrittweite
Ur Referenzgeschwindigkeit
ei
U gemittelten Geschwindigkeit in Indexnotation
~
U Geschwindigkeitsvektor
u, v, w x-, y- und z-Komponente des Geschwindigkeitsvektors
uexit Geschw. der Strömung am Düsenauslass in x-Richtung
uexit,id ideale Geschw. der Strömung am Düsenauslass in x-Richtung
ui Momentangeschwindigkeit in Indexnotation
u′ i flukt. Anteil (Reynolds) der Momentangeschwindigkeit in Indexnota-
tion
u′′ i flukt. Anteil (Favre) der Momentangeschwindigkeit in Indexnotation
V Volumen
Ẇ Arbeit des Fluids im Kontrollvolumen
~
W Vektor konservativer Variablen
x, y, z kartesische Raumkoordinaten
xi Positionsvektor in Tensornotation
∆x charakteristische Zellänge
y+ dimensionsloser Wandabstand
viii
Griechische Buchstaben
ℵ(ui ) Navier-Stokes Operator

α laminare Temperaturleitfähigkeit
αT turbulente Temperaturleitfähigkeit
β Runge-Kutta Koeffizient
δij Kronecker Delta
ǫ turbulente Dissipation
ǫjθ molekulare Dissipation von u′′j θ
Γ Precondition-Matrix
Γnc Precondition-Matrix für nichtkonservatives System
κ Isentropenexponent
λ laminare Wärmeleitfähigkeit
λT turbulente Wärmeleitfähigkeit
µ laminare dynamische Viskosität
µT turbulente dynamische Viskosität
ν laminare kinematische Viskosität
νT turbulente kinematische Viskosität
ω spezifische turbulente Dissipation
Φ∗jθ Pressure-Scrambling von u′′j θ
ρ Dichte
Σ(c) Transportgleichungs-Operator für Skalar c
σ Spannungstensor
σij Spannungstensor in Indexnotation
σx Dichte der Oberflächenkräfte in x-Richtung
τij spezifischer Reynolds Stresstensor
Θ Precondition-Parameter
θ fluktuierender Anteil der statischen Temperatur
ϕ
~ Richtungsvektor für Geschwindigkeit am Einlass
Alle physikalische Größen sind in SI-Einheiten angegeben
ix
Abkürzungen
1stO Genauigkeit erster Ordnung im Raum

2ndO Genauigkeit zweiter Ordnung im Raum
CAD Computer-Aided Design
CNPR Nebenstromdüsendruckverhältnis
DES Detatched Eddy Simulation
DNS Direct Numerical Simulation
EDM Eddy Diffusivity Model
FANS Favre gemittelte Navier-Stokes Gleichungen
ISA International Standard Atmosphere
LES Large Eddy Simulation
OT OPTITHECK
PS Nebenstrom- zu Kernstromtotaldruckverhältnis
RANS Reynolds gemittelte Navier-Stokes Gleichungen
RR Rolls-Royce Düsengeometrie
RSM Reynolds Stress Modell
SFC spezifischer Treibstoffverbrauch
SKE standard k-ǫ Turbulenzmodell
SST Shear Stress Transport
STP Standard for the Exchange of Product Model Data
TC Temperaturkorrekturverfahren von Abdol-Hamid
TC+DSM gekoppeltes Modell aus Temperaturkorrekturverfahren und DSM
TKE turbulente kinetische Energie
TR Kern- zu Nebenstromtotaltemperaturverhältnis
VTK Visualization Toolkit
x
Kurzfassung
In der vorliegenden Arbeit wird das Mischungsverhalten von Abgassystemen mit inte-
griertem Blütenmischer numerisch untersucht. Zum Einsatz kommen dabei Verfahren
und Ansätze, die eine Verbesserung der Beschreibung anisothermer Phänomene berück-
sichtigen. Die Implementierung der Modelle erfolgte in den Open Source Strömungslöser
OpenFOAM, mit welchem auch die Großzahl der Berechnungen durchgeführt wurden.
Bei den bisher gängigen Reynolds bzw. Farve gemittelten Berechnungsverfahren spielt
die Fluktuation von Dichte und Temperatur praktisch keine Rolle. Diese Schwankungs-
größen können aber vor allem bei stark anisothermen Strömungen einen entscheidenden
Einfluss haben und auf die Entstehung und Evolution von Turbulenz und somit auch
auf das Strömungsbild einwirken. Durch ein Temperaturkorrekturverfahren wird die-
ser Einfluss nachmodelliert indem in Bereichen großer Totaltemperaturgradienten die
turbulente Viskosität und somit die Ausmischung in der anisothermen Scherschicht er-
höht wird. Die höhere Mischungsrate zwischen Bypass- und Kernstrom wirkt sich aber
auch auf das Leistungsverhalten der Düse aus. So wurde dargestellt, dass die größere
Ausmischung der beiden Ströme und somit eine homogenere Temperaturverteilung am
Düsenausgang, zu einem Anstieg des Geschwindigkeitskoeffizienten führt. Gleichzeitig
nimmt der Wert des Durchflusskoeffizienten aufgrund von auftretenden Verlusten ab.
Vergleicht man die Leistungsdaten mit experimentell ermittelten Absolutwerten liegen
die Koeffizienten tendenziell etwas unterhalb der Messungen. Allerdings wird durch das
Temperaturkorrekturverfahren der qualitative Verlauf der Kurven besser wiedergegeben.
Des Weiteren wurde der bisher gängige turbulente Prandtlzahlansatz durch ein differen-
tielles Modell zur Berechnung des turbulenten Wärmestroms ersetzt. Mit diesem Ver-
fahren ist eine Bestimmung und Festlegung einer turbulenten Prandtlzahl nicht mehr
nötig. Außerdem können dadurch auftretende physikalische Phänomene besser beschrie-
ben werden. Die Simulationsergebnisse zeigen, dass das differentielle Verfahren lokal
differenzierter mischt als das turbulente Prandtlzahlmodell. So erhöht sich bei diesem
Ansatz folgerichtig die thermische Mischung in stark turbulenten dreidimensionalen Ge-
bieten gegenüber Bereichen im Strömungsfeld, wo eine eher zweidimensionale und ge-
schichtet verlaufende Düsenströmung vorherrscht. Allerdings wurde beobachtet, dass die
Konvergenz dieses Ansatzes ist im Vergleich zum algebraischen Prandtlzahlansatz leicht
schlechter ist und es wurde eine erhöhte Totalenthalpieflussdifferenz zwischen den Ein-
und Auslassflächen festgestellt. Ursache könnte die geringfügig stärkere Aufheizung der
xii
Kurzfassung
Wandgrenzschicht beim differentiellen Modell sein. Der festgestellte leichte Versatz der
Leistungskurven von Geschwindigkeits- und Durchflusskoeffizient gegenüber dem turbu-
lenten Prandtlzahlansatz könnte damit möglicherweise erklärt werden.
In einem nächsten Schritt wurde eine Kombination des Temperaturkorrekturverfahrens
mit der differenziellen Wärmestromberechnung ausgearbeitet. Zur Kalibrierung des Mo-
dells wurde auf gegebene Totaltemperaturmessungen zurückgegriffen. Obwohl mit dem
neuen Modell die thermische Ausmischung deutlich gesteigert wurde, konnte die stark
diffusive Totaltemperaturverteilung der Düsenströmung aus den Messungen nicht repro-
duziert werden. Allerdings zeigt ein Vergleich mit experimentellen Totaldruckmessungen,
dass bereits bei den gängigen Standardmodellen diese Werte deutlich besser mit den Si-
mulationsergebnissen korrelieren. Es muss deshalb gefolgert werden, dass eventuell auch
die Totaltemperaturmessungen in Frage gestellt werden sollten.
Abschließend wurde noch mit dem kommerziellen Strömungslöser FLUENT eine Para-
meteruntersuchung durchgeführt. In dieser wurde zum einen der Einfluss des Rechengit-
ters auf die Ergebnisse untersucht. Während die qualitativen Unterschiede überschaubar
gering sind, lassen sich für feinere Gitter sowohl größere Durchfluss- als auch Geschwin-
digkeitskoeffizienten erzielen. Ferner wurde noch der geometrische Einfluss der Mischer-
taktung untersucht. Mit einer periodischen statt einer symmetrischen Randbedingung
an den Seitenränder des Rechengebiets konnte bei ansonsten identischem Netz eine ge-
änderte Mischergeometrie simuliert werden. Während sich qualitativ dadurch das Strö-
mungsbild deutlich ändert, bleiben die Leistungsbeiwerte der Düse aber unverändert. Da
der verwendete Löser in OpenFOAM nur mit einer konstanten spezifischen Wärmeka-
pazität kompatibel ist, wurde auch dieser Einfluss noch untersucht. Im Ergebnis zeigte
sich, dass für eine von der Temperatur abhängige Wärmekapazität bei einem Totaltem-
peraturverhältnis zwischen Kern- und Nebenstrom größer eins der Durchflusskoeffizient
gegenüber einem konstanten Wert sinkt und der Geschwindigkeitskoeffizient steigt. Ein
nennenswerter qualitativer Unterschied in der Lösung ist jedoch nicht erkennbar. Letzt-
lich wurde noch der Einfluss des Turbulenzmodells und dessen räumliche Diskretisierung
untersucht. Für den Vergleich wurde das k-ǫ und k-ω-SST Modell verwendet. Bei einer
Diskretisierung des Turbulenzmodells erster Ordnung zeigte sich, dass die Ausmischung
des Kern- und Nebenstroms für das k-ǫ Modell am schwächsten ausfällt, wohingegen sie
für das k-ω-SST am stärksten ist. Bei zweiter Ordnung nähern sich beide Turbulenz-
modelle an und das Ergebnis ist fast identisch. Die qualitative Beobachtung bezüglich
der Ausmischung spiegelt sich auch quantitativ in den berechneten Leistungsbeiwerten
wider.
xiii
Abstract
In this thesis, the numerical investigation of mixing hot and cold streams in an exhaust
nozzle system with installed lobed forced mixer is presented. For an improved simulati-
on of such flows, advanced models are used which better incorporate high anisothermal
phenomena. The implementation of the new models was conducted in the open source
flow solver OpenFOAM, which was also used to perform the most of the simulations.
In the commonly used Reynolds or Favre averaged equations, the fluctuation of density
or temperature is not considered. However these fluctuating quantities can play an role
important in high anisothermal flows and may effect the production and evolution of
turbulence and finally impact the flow. Therefore, a temperature correction method was
implemented to take into account such effects. The model increases the turbulent visco-
sity in regions, where high total temperature gradients occur and thus enhance mixing
in anisothermal shear layers. The higher mixing rate between bypass and core flow also
influences the performance of the nozzle. It was shown, that the increased mixing and
the therefore more homogenous temperature distribution at the nozzle exit, leads to
increased velocity coefficients. Simultaneously the discharge coefficient decreases due to
increasing losses. If the numerical results are compared to experimental measurements
it can be noticed, that the coefficients from CFD are slightly below the measured values.
However with the temperature correction method, the qualitative distribution of the
curves is better.
Furthermore, the commonly used turbulent Prandtl number approach hast been sub-
stituted by a differential model to calculate the turbulent heat-flux. With this method,
it is not necessary to define or calculate a turbulent Prandtl number. Moreover, this
approach has potential to better describe the physics than the turbulent Prandtl num-
ber approach does. The simulation results show that the differential model mixes locally
more differentiate than the Prandtl number model. In regions where the flow is turbulent
and three dimensional, the thermal mixing of the differential model is higher compared
to zones where the character of the flow is more layered and two dimensional. However it
was noticed that convergence of the differential approach seems to be slightly worse than
for the algebraic Prandtl number method. The calculated positive total enthalpy flux
difference between fan plus core inlets and nozzle exit of the advanced model is higher.
A reason for this extra amount of energy within the nozzle control volume might be the
slightly stronger heating of the wall shear layer for this model. This might also explain
xiv
Abstract
the little offset of the curves of the performance parameters compared to the algebraic
model.
In a next step, a combination of the temperature correction method and the differential
heat-flux model was elaborated. To calibrate the new model, total temperature measure-
ments were given. Although for the new model the thermal mixing of the streams could
be highly increased, it was not possible do reproduce the strongly diffusive profile of the
total temperature distribution from the measurements. However the comparison of CFD
data to experimental total pressure plots shows clearly better agreement, even for the
commonly used standard models. This leads to the conclusion, that the performed total
temperature measurements might be questioned as well.
Finally, a parametric study has been conducted with the commercial solver FLUENT.
Here the influence of the mesh was studied. While the qualitative influence of the num-
ber of cells within the domain was little, there was more impact on the performance
parameters. For a finer mesh, higher velocity and discharge coefficients can be achieved.
Furthermore the influence of the mixer geometry has been studied. Using a periodic
instead of a symmetric boundary condition at the domain bounds, for the existing mesh
an alternative clocking of the mixer lobes can be simulated without changing the mesh.
It was demonstrated how the change of mixer geometry effects the evolution of the flow.
However there was no essential influence on performance parameters. As the used flow
solver in OpenFOAM is only compatible with a constant specific heat capacity, the effect
of a temperature depending value was tested. Wile the qualitative differences between
both cases is negligible, the performance parameters behave different. For increasing
total temperature ratios between core and fan stream, the velocity coefficient is higher,
but lower for the discharge coefficient if using variable heat capacity. In the end the
influence of the turbulence model and its spatial discretization scheme was studied. For
this comparison, the k-ǫ and k-ω-SST model was used. If discretising both models first
order accurate, the mixing of hot core and cold bypass streams is weakest for the k-ǫ
model, however strongest for the k-ω-SST. If switching to second order accuracy, k-ǫ and
k-ω-SST deliver almost identical results. This behavior is also mirrored in the computed
performance parameters.
xv
1. Einleitung
1.1. Einführung
Im Zuge abnehmender fossiler Ressourcen und der zunehmenden Umweltbelastung durch
konventionelle Verbrennungsmotoren wird die Forderung nach innovativen Antriebskon-
zepten immer größer. Steigende Benzin- oder Kerosinpreise stehen einem stetig wach-
senden Personen- und Warenbeförderungsaufkommen gegenüber. Ebenso fand vor allem
in den großen Industrieländern ein verstärktes ökologisches Bewusstsein Einzug. Beides
lässt den Ruf nach effizienten, emmisionsarmen und nachhaltig betriebenen Antriebsfor-
men lauter werden.
Stark gefordert mit zusätzlichen Entwicklungszielen ist vor allem auch die Luftfahrt-
branche. Während die ersten Turbojettriebwerke hauptsächlich den thermodynamischen
Vorteil des Gasturbinenprozesses gegenüber den Kolbenmotoren demonstrierten, wiesen
diese im Vergleich zu modernen Strahltriebwerken noch eine verhältnismäßig bescheide-
ne Effizienz auf. Durch den Gang der letzten Jahrzehnte hin zu immer größer werdenden
Bypassverhältnissen wurde eine starke Reduzierung des spezifischen Treibstoffverbrauchs
(SFC) realisierbar. Abbildung 1.1 stellt den Zusammenhang zwischen unterschiedlichen
Triebwerksgeneration und dem spezifischen Treibstoffverbrauch über die Jahre dar.
Eine immer gewichtigere Rolle bei der Entwicklung moderner Flugtriebwerke spielt ne-
ben der ökonomischen Bilanz die Lärmemission. Zwar ging die Vergrößerung des By-
passverhältnisses auch mit einer Reduzierung des Strahllärms einher, allerdings muss die
weitere Verringerung der Lärmemission noch weiter vorangetrieben werden. Der stetig
wachsende Flugverkehr, vor allem in den großen und dicht besiedelten Ballungsräumen,
führt zu einer immer größer werdenden Belastung durch den durch Flugtriebwerke verur-
sachten Lärm. So gehen manche Flughafenbetreiber dazu über, ihre Landegebühren nach
Lärmkriterien zu staffeln. Das heißt also, dass die Lärmemission der Fluggeräte direkt an
höhere Kosten für die Airlines gebunden ist. Den steigenden Kostendruck versuchen die
Fluggesellschaften wiederum zum Teil an die Flugzeug- und Triebwerkshersteller weiter-
zugeben, indem sie diese auffordern deren Forschungs- und Entwicklungsanstrengungen
dahingehend zu erhöhen.
Der Forderungskatalog an die Luftfahrtbranche umfasst somit in Zukunft nicht nur eine
1
Einleitung
Abbildung 1.1.: Relativer spezifischer Treibstoffverbrauch unterschiedlicher Triebwerks-

generationen bezogen auf das erste Tubojettriebwerk [1].
weitere Steigerung der Effizienz und der Reduktion des Treibstoffverbrauchs wegen rein
ökonomischen Gründen, sondern neuerdings auch ökologische Aspekte. Unter Letztere
fallen einerseits die Ziele eines geringeren Schadstoffausstoßes oder des ressourcenscho-
nenden Betriebs der Triebwerke zur Erhaltung einer gesunden Umwelt. Andererseits
gehört dazu auch die Lärmentwicklung zu reduzieren, um so die Belastung für Mensch
und Umwelt so gering wie möglich zu halten.
Ein Großteil der Lärmemission eines Flugzeugs als Gesamtsystem geht von den Triebwer-
ken aus [2]. Folglich liegt ein Groß des Optimierungspotentials in den Antriebssystemen.
Allerdings ist die Entwicklung moderner Antriebssysteme für Flugkörper heute bereits
weit voran geschritten. Möglichkeiten zur Optimierung konventioneller Luftstrahltrieb-
werke finden sich daher heute zu großen Teilen in Detaillösungen wieder.
Viele Lösungen zur Lärmminderung beeinflussen unter anderem direkt den Triebwerks-
prozess. Eine Reduzierung des Strahllärms bei gleichzeitigem gleich halten der Leistungs-
daten und ohne den Antriebszyklus zu beeinflussen, war vor allem bei Triebwerken mit
getrenntem Abgassystem jahrelang eine große Herausforderung. Ein geeignetes Konzept
hierfür scheint jedoch mit dem Einsatz von Chevrondüsen gefunden zu sein. Unter Ver-
wendung dieser gezackten Düsenhinterkanten kann der Strahllärm signifikant gesenkt
werden, während die Leistungseinbußen dadurch minimal bleiben [2]. Beim Einsatz von
2
Einleitung
Chevrons wird der Fan- und der Umgebungsstrom oder Kern- und Fanstrom deutlich
besser vermischt als bei konventionellen gerade verlaufenden Düsenhinterkanten [3, 4].
Dadurch werden die Scherschichten zwischen den Luftstrahlen früher aufgebrochen, die
Luft besser vermischt und letztlich der Strahllärm gesenkt [5]. Durch entsprechendes
Takten der Verzahnungsreihen kann der gewonnene Vorteil durch die Chevrons opti-
miert werden [6]. Abbildung 1.2 zeigt ein Triebwerk mit getrenntem Abgassystem und
Chevrondüsen.
Abbildung 1.2.: Triebwerk mit Chevrondüsen.
Bei Zweistromtriebwerken mit mittlerem Bypassverhältniss werden heutzutage haupt-

sächlich interne Mischer eingesetzt. Solche Triebwerke kommen typischerweise bei Ge-
schäftsflugzeugen (z. B. Gulfstream G650), als auch bei Kurz- und Mittelstreckenflug-
zeugen (z. B. Airbus A320, Boeing 717) zum Einsatz. Durch die Verwendung von soge-
nannten Blütenmischern, oder zu englisch Lobed Mixer, wird durch gezieltes ineinander
lenken von heißem Kernstrom und kaltem Fanstrom eine starke Durchmischung beider
Strahlen erzwungen. Dabei ist ein Teil der Mischungsweglänge im Nachlauf des Mischers
vermantelt und befindet sich somit innerhalb des für den Schub verantwortlichen Kon-
trollvolumens des Triebwerks. Dadurch kann neben der Reduzierung des Strahllärms
auch eine Schubsteigerung erzielt werden [7, 8]. Dies wurde bereits früh erkannt und
schon in den 60er Jahren versucht, die Vorteile von interner Mischung umzusetzen [9].
Allerdings muss beachtet werden, dass durch die Strahlmischung auch Mischungsver-
luste entstehen [10]. Durch eine sorgfältige Auslegung des Mischer- und Düsendesigns
3
Einleitung
kann aber die Mischung optimiert und resultierende Druckverluste können minimiert
werden. Ferner führt die zusätzliche Verkleidung bei zunehmender Größe des Triebwerks
zu höher werdenden Gewicht und Widerstand, so dass der Vorteil durch interne Mi-
schung ab einer bestimmten Triebwerksgröße geringer wird und man daher wieder zur
offenen Strahlmischung übergeht. Abbildung 1.3 zeigt ein Beispieltriebwerk mit internem
Mischer.
Abbildung 1.3.: Triebwerk mit internem Mischer aus der BR700-Serie von Rolls-Royce.
1.2. Stand der Forschung

Eine immer tragendere Rolle bei der Entwicklung und Optimierung von Antriebssyste-
men nehmen rechnergestützte Verfahren ein und sind aus der Prozesskette der Entwick-
lung nicht mehr wegzudenken.
Trotz der Verwendung moderner Hochleichstungsrechner stößt man bei der Berechnung
komplexer Systeme jedoch schnell an die Grenzen des Möglichen. Um die physikalisch
limitierte Rechenleistung zu kompensieren, werden Kompromisse bei der Rechengenauig-
keit in Kauf genommen. Die zur Beschreibung der physikalischen Phänomene hinterleg-
4
Einleitung
ten Gleichungen werden daher zum Teil sehr stark vereinfacht und durch Modellannah-
men angenähert. Grundsätzlich lassen sich die Berechnungsmethoden in drei Kategorien
einteilen [11]:
• DNS (Direct Numerical Simulation)
• LES (Large Eddy Simulation)
• RANS (Reynolds Averaged Navier-Stokes Gleichungen)
Bei der DNS werden die Erhaltungsgleichungen räumlich und zeitlich vollständig vom
Rechner gelöst [12, 13]. Dies setzt aber ein hinreichend feines Gitter zur adäquaten Auf-
lösung aller auftretenden Phänomene voraus. Außerdem kann DNS momentan nur bei
sehr kleinen Reynoldszahlen eingesetzt werden [14]. Der Rechenaufwand zur Lösung der
Gleichungssysteme ist für die meisten Anwendungen deshalb deutlich zu hoch.
Bei LES und RANS nimmt man eine gewisse Ungenauigkeit beim Lösen der Gleichungen
in Kauf. Sie erlauben auch eine geringere Auflösung des Gitters. Die aufgrund des grö-
beren Rechennetzes nicht berücksichtigte Physik wird nun durch Modelle angenähert.
Während bei LES die Gleichungen bis zu einer vorgegebenen räumlichen Auflösung
vollständig gelöst werden, werden die Gleichungen bei RANS von vorn herein zeitlich
gemittelt, so dass fluktuierende Größen eliminiert werden. Um diese „verlorenen“ Infor-
mationen bei der RANS durch die Mittelung und deren Einfluss auf die Physik teilweise
wieder zurück zu gewinnen, wird bei den sehr populären Zweigleichungsmodellen der
RANS zumindest die Schwankung der Geschwindigkeitskomponenten bzw. die daraus
resultierende turbulente kinetische Energie mit deren Einflüssen nachmodelliert. Die
LES löst die Wirbel dagegen bis zu einer bestimmten Größe vollständig auf, die kleine-
ren Wirbel werden hingegen modelliert [12, 15].
Aufgrund des hohen Rechenaufwands bei der Verwendung von DNS, findet dieses Verfah-
ren in der industriellen Nutzung kaum Anwendung [13]. Weiter verbreitet sind dagegen
LES und vor allem RANS, wobei auch Ersteres für viele Ingenieursanwendungen wegen
des hohen Rechenaufwands eine eher untergeordnete Rolle spielt. RANS hingegen hat
sich für sehr viele Problemstellungen als ausgezeichneter Kompromiss hinsichtlich Ge-
nauigkeit und rechnerischem Aufwand herausgestellt. Selbst äußerst komplexe Probleme
lassen sich mit diesem Verfahren mit angemessenem Aufwand lösen. Zwischen LES und
RANS ist noch die Detatched Eddy Simulation (DES) Methode angesiedelt. Die DES
ist ein hybrides Verfahren und eine Mischung von RANS und LES. Aber auch die DES
ist für die meisten Ingenieursanwendungen mit zu großem Rechenaufwand verbunden.
Speziell die Ansprüche an einen Strömungslöser zur Berechnung von Flugzeugdüsenströ-
mungen sind sehr hoch und es muss ein breites Anforderungsspektrum abgedeckt werden.
Das Verfahren muss zum einen in der Lage sein, stark kompressible Strömungen lösen
zu können. In den meist transonischen Strömungen können auch lokal Machzahlen von
5
Einleitung
über eins erreicht werden, wo Phänomene wie Verdichtungsstöße auftreten. Zum anderen
können auch bei generell hohen Strömungsgeschwindigkeiten im Rechengebiet vereinzelt
Bereiche mit sehr geringen Machzahlen auftreten. So z. B. in Rezirkulationsgebieten, in
Grenzschichten oder beim Austritt der Düsenströmung in ein ruhendes Fluid.
DNS kommt bei solchen Düsenströmungen auf Grund des hohen Rechenaufwands nicht
zum Einsatz. Auch die Verwendung von LES-Methoden ist eher selten und wenn dann
meist nur zu akademischen Zwecken. Das gängigste Verfahren ist das Lösen von zeit-
lich gemittelten Gleichungen. Da die Düsenströmungen hauptsächlich kompressibel sind,
werden die Navier-Stokes Gleichungen zusätzlich zur zeitlichen Mittelung auch dichte-,
bzw. massengemittelt. Diese sogenannten Favre-gemittelten Navier-Stokes Gleichungen
(FANS) fallen in die Klasse der RANS, wobei RANS häufig auch als Synonym im Bezug
auf die kompressiblen gemittelten Erhaltungsgleichungen benutzt wird, wenn eigentlich
von Favre Mittelung die Rede ist. Die Mittelung, egal welcher Art, führt zu zusätzlichen
Termen in den Navier-Stokes Gleichungen, die mit Hilfe geeigneter Turbulenzmodelle
geschlossen werden müssen. Die häufigste Anwendung finden hier wie bereits erwähnt
die populären Zweigleichungsturbulenzmodelle.
Das Mischungsverhalten von Kern- und Nebenstrom eines Flugtriebwerks spielt für des-
sen Leistungsverhalten sowie Lärmemission eine wesentliche Rolle. So kann eine verbes-
serte Ausmischung des heißen Kernstroms mit dem Nebenstrom zu einer Schuberhöhung
führen sowie die durch den Strahl verursachte Lärmemission reduzieren [7]. Genaue nu-
merische Verfahren sind daher von hoher Bedeutung zum Studium der auftretenden
Phänomene und für eine verlässliche Vorauslegung eines Triebwerks.
Generell lässt sich sagen, dass die Mischungsrate von Triebwerksstrahlen von den RANS
Modellen unterschätzt wird. Bei heißen Triebwerksstrahlen ist dieser Fehler besonders
groß [16, 17]. Dabei ist z. B. die Vorhersage der Länge des so genannten heißen Poten-
tialkerns besonders wichtig für die Lärmbestimmung [18, 19].
Ein Grund für dieses Defizit liegt unter anderem darin, dass die verschiedenen Zwei-
gleichungsmodelle bei RANS fast ausschließlich darauf abzielen, die Wiedergabe des
Strömungsfeldes hinsichtlich der kinetischen Energie zu verbessern. Die bekanntesten
Vertreter hierfür sind die Zweigleichungsturbulenzmodelle k-ω, k-ǫ oder das Hybridmo-
dell k-ω-SST. Sie liefern für kalte Strahlen im niedrigen Machzahlbereich gute Lösungen.
Für Hochtemperaturstrahlen kommen hier jedoch wesentliche physikalische Größen gar
nicht erst zum Tragen. So werden zum Beispiel die Schwankungsgrößen von Tempera-
tur und Dichte quasi nicht berücksichtigt, da diese durch die Favre-Mittelung eliminiert
werden.
Seiner et al. [20] and Thomas et al. [21] zeigen in ihren Untersuchungen jedoch, dass
hohe Totaltemperaturgradienten zu einer schnelleren Vermischung und einer stärkeren
Aufspreizung eines Strahls führen. Ferner zeigen Tam und Ganesan [22] in einer Stabili-
tätsanalyse, wie die Dichtedifferenz zweier Fluide zu einer Kelvin-Helmholtz-Instabilität
6
Einleitung
führt und die Mischungsscherschicht beeinflusst. Abdol-Hamid [18] nennt direkt im Zu-
sammenhang mit der Präsenz hoher Temperatur- bzw. Dichtegradienten das Auftreten
von Turbulenz. Die Fluktuation von Druck und Temperatur führt zwangsläufig zu einer
Dichteschwankung und wirkt sich somit als Quelle von Turbulenz auf das Geschwindig-
keitsfeld aus.
Ein weiterer Grund für die unterschätzte Ausmischung von Heiß- und Kaltstrahlen liegt
an dem zumeist verwendeten einfachen Modell für den turbulenten Wärmestrom. Zur
Schließung des in der Energiegleichung auftretenden turbulenten Wärmestromterms wird
fast ausschließlich das Eddy Diffusivity Modell (EDM) herangezogen. Bei diesem Mo-
dell wird angenommen, dass der turbulente Wärmestrom proportional zum mittleren
Temperaturgradienten ist und setzt die Kenntnis einer turbulente Prandtlzahl P rT vor-
aus. In den meisten Fällen wird diese Größe als konstant angenommen. Jedoch variiert
diese turbulente Prandtlzahl je nach Strömungszustand. Nach Wilcox [13] und Birch et
al. [19] können die Werte im Strömungsgebiet zwischen 0,9 und 0,4 liegen. Ein einzi-
ger konstanter Wert ist demzufolge nicht in der Lage alle Bereiche im Strömungsgebiet
adäquat abzudecken. So gibt es mehrere Versuche unterschiedlicher Autoren einen va-
riablen Prandlzahlansatz zu finden, sowohl als algebraisches Modell [23, 24, 25, 26, 27],
als auch in differentieller Form [28, 29]. Allerdings hat das EDM einen entscheidenden
Nachteil selbst wenn die turbulente Prandtlzahl en jedem Ort des Rechengebiets kor-
rekt bestimmt wird. So ist das Eddy Diffusivity Modell quasi nur eindimensional, d.h.
turbulente Wärme wird nur in eine Richtung transportiert, von hohen zu niedrigen Tem-
peraturbereichen. Experimentelle Messungen [30, 31] haben jedoch ergeben, dass selbst
wenn der Temperaturgradient in Strömungsrichtung viel kleiner ist als normal dazu, der
turbulente Wärmestrom parallel zur Fließrichtung um ein Vielfaches größer sein kann.
In der Summe führen die oben genannten Gründe bei den Standardturbulenzmodellen
und den Standardschließungsansätzen so zu einer zu geringen Ausmischung in anisother-
men Scherschichtströmungen und folglich zu Ungenauigkeiten in den Ergebnissen.
1.3. Zielsetzung der Arbeit

Gegenstand der Arbeit ist die qualitative und quantitative Untersuchung der Strahl-
mischung im Nachlauf eines Blütenmischers. Dabei sollen unterschiedliche numerische
Modelle zur verbesserten Berechnung in solchen anisothermen Strömungen zum Einsatz
kommen.
Zum einen sollen bereits existierende Modelle eingesetzt werden, wie z. B. ein Tempera-
turkorrekturverfahren zur Berücksichtigung von Temperaturfluktuationen. Ebenso soll
eine bessere Formulierung des turbulenten Wärmestromvektors verwendet werden, ba-
sierend auf einem differentiellen Ansatz, welcher keine Definition bzw. Berechnung einer
7
Einleitung
turbulenten Prandtlzahl benötigt. Die Implementierung der Modelle soll in den Open
Source Code OpenFOAM erfolgen. Die unterschiedlichen Ansätze sollen dann anhand
von ihrem Verhalten und ihrer Wirksamkeit auf die Strahlmischungssimulation analy-
siert und bewertet werden. Die Ergebnisse werden neben der qualitativen Auswertung
auch hinsichtlich leistungsrelevanter Parameter wie dem Durchflusskoeffizient CD und
dem Geschwindigkeitskoeffizient CV für verschiedene Betriebszustände untersucht. Die
numerischen Ergebnisse können mit experimentell ermittelten Leistungsdaten verglichen
werden.
Später soll ein neues Modell entwickelt werden, welches die Vorteile des differentiellen
Ansatzes mit dem Temperaturkorrekturverfahren zur Berücksichtigung von Temperatur-
schwankungen kombiniert. Das neue Modell kann mit weiteren experimentellen Messun-
gen verglichen und validiert werden. So stehen hierfür experimentelle Totaltemperatur-
und Totaldruckmessungen zur Verfügung.
In einer Parameteruntersuchung sollen des Weiteren unterschiedliche numerische und
physikalische Einflüsse auf die Strömungsberechnung untersucht werden. Diese Fallstu-
die soll mit dem kommerziellen Strömungslöser ANSYS FLUENT durchgeführt werden.
So soll zum einen eine Gitterstudie realisiert werden, um den Einfluss der Anzahl und
Größe der Gitterzellen auf das Ergebnis zu untersuchen. Ferner soll die Auswirkung un-
terschiedlicher Mischertaktungen auf die Leistungskoeffizienten überprüft sowie die Wahl
und die Diskretisierungsform des Turbulenzmodells analysiert werden. Da der eingesetz-
te Strömungslöser in OpenFOAM nur mit einer konstanten spezifischen Wärmekapazität
kompatibel ist, soll in FLUENT auch der Einfluss eines temperaturabhängigen Werts
untersucht werden.
8
2. Prinzip und Leistungsbewertung
von Schubdüsen mit internem
Mischer
2.1. Grundprinzip interner Strahlmischung

In Triebwerken mit mittlerem Nebenstromverhältnis werden zur Steigerung der Effizienz
heutzutage meist gemischte Abgassysteme eingesetzt. Durch eine effektive Mischung des
heißen Kernstroms mit dem kalten Nebenstrom innerhalb des ummantelten Düsensys-
tems kann die Leistungsfähigkeit des Triebwerks gesteigert werden. Der Einsatz gemisch-
ter Abgassysteme bietet prinzipiell folgende Potentiale, die einzeln oder in Kombination
mehrerer Potentiale ausgenutzt werden können [7, 9, 32, 33, 34]:
• Schuberhöhung
• Reduzierung des Strahllärms
• Erhöhung des Vortriebswirkungsgrads
• Reduzierung des spezifischen Treibstoffverbrauchs
• Reduzierung der Einbaulänge des Abgassystems
• Reduzierung der Infrarotstrahlung
Um die thermische Ausmischung zwischen Heiß- und Kaltstrahl zu steigern, kommen
Zwangsmischer zum Einsatz. Ihre meist mäanderförmige Kontur erinnert beim Blick
von hinten in das Düsensystem an eine Blumenblüte, weshalb sie auch als Blütenmi-
scher bezeichnet werden. Diese Mischer erzeugen in ihrem Nachlauf sowohl eine größe-
re Scherschicht als auch ein komplexes Wirbelsystem, und mischen so die heißen und
kalten Strömungsbereiche in einem größeren Gebiet als es mit molekularen oder tur-
bulenten Transport möglich ist [7]. Ein Zugewinn an Schub lässt sich im Prinzip erst
durch einen Temperaturunterschied zwischen Haupt- und Nebenstrom erzeugen [32].
9
Prinzip und Leistungsbewertung von Schubdüsen mit internem Mischer
Der Grund hierfür liegt in der Divergenz der Linien konstanten Druckes im Enthalpie-
Entropie-Diagramm mit steigender Enthalpie [35]. Durch die Übertragung der höheren
thermischen Energie vom Kernstrom auf den Nebenstrom lässt sich somit eine Stei-
gerung des Schubes erzielen [7]. Bei idealer Mischung und ohne jegliche Verluste ist
durch Ausnutzung dieses thermodynamischen Effekts eine theoretische Treibstoffeinspa-
rung von 5-7% möglich [36]. Unter Berücksichtigung von Mischungsverlusten beträgt das
Einsparpotential noch ca. 3-5% [32, 37, 38]. Gleichzeitig sinkt die Verlustleistung des Ab-
gasstrahls durch die quadratische Abhängigkeit von dessen Geschwindigkeit [39, 40]. Als
Folge steigt dadurch der Vortriebswirkungsgrad des Triebwerks [7]. Real ist es jedoch
schwer sich den genannten Werten anzunähern und die Potentiale voll auszuschöpfen.
Selbst wenn es gelingt die viskosen und thermodynamischen Verluste in Grenzen zu
halten geht das Prinzip gemischter Abgassysteme, und im Speziellen mit dem Einbau
von Blütenmischern, mit einer Gewichtszunahme und gesteigerten Kosten gegenüber
vergleichbaren offenen Systemen einher. Es muss also bei der Auslegung genau darauf
geachtet werden, dass der thermodynamische Nutzen der durch die Mischung der heißen
und kalten Strahlen entsteht in angemessenem Verhältnis zur erhöhten Systemkomple-
xität und den damit verbundenen Ausgaben bleibt.
Ein weiterer großer Vorteil gemischter Abgassysteme ist die dadurch mögliche Lärm-
reduzierung [41, 42, 43]. Für ein Triebwerk mit mittlerem Bypassverhältnis geht der
Strahllärm hauptsächlich vom Heißstrahl und seiner hohen Geschwindigkeit aus. Die In-
tensität des Strahllärms steigt mit der achten Potenz der Strahlgeschwindigkeit an [44].
Bei heißen Strahlen kann dieser Exponent kleiner werden und bis zur sechsten Potenz
sinken. Durch den durch die Mischung erzielten Energietransfer vom heißen, schnel-
len Kernstrom auf den Nebenstrom, wird der Heißstrahl verlangsamt. Der Breitband-
Strahllärm, der durch die turbulenten Geschwindigkeits- und Dichteschwankungen in der
Scherschicht entsteht, sinkt dadurch stark überproportional zu den Geschwindigkeitsän-
derungen [7]. Ferner wird durch die Mischung die Oberfläche des Strahls vergrößert. Dies
bewirkt eine Frequenzverringerung, so dass ein Teil der Lärmfrequenzen in einen Bereich
verschoben wird, der unterhalb der menschlichen Hörschwelle liegt. Zum anderen wer-
den hohe Frequenzen in Bereiche verschoben, in denen sie von der Atmosphäre besser
absorbiert werden können [41]. Durch die reduzierte Lärmemission können z. B. in eini-
gen Ländern die Landegebühren gesenkt und darüber hinaus der durch die Lärmgrenzen
bedingte zulässige maximale Startschub erhöht werden [32].
Die Mischer wurden in den letzten Jahren immer weiter hinsichtlich Ausmischung, Leis-
tungsverhalten und Lärmemission untersucht und optimiert [45, 46, 47, 48, 49, 50, 51]. Zu
erwähnen sind auch die aus zahlreichen Untersuchungen hervorgegangenen Mischerfor-
men wie das „scalloping“ und „scarfing“. Bei Ersterem wird an den Mischerseitenrändern
Material ausgespart, so dass es zu einer früheren und somit verstärkten Mischung von
Heiß- und Kaltstrahl kommt. Diese Matereialaussparung an den Flanken der Mischer-
blüten kann relativ extrem ausfallen, deshalb werden diese Mischer aufgrund ihrer Form
10
auch „Boomerang Scallops“ genannt [46]. Beim scarfing werden die Mischerblüten alter-
nierend schräg nach hinten geschnitten, was ebenfalls die Vermischung verstärken soll.
Gleichzeitig führt die Materialersparnis in beiden Fällen zu einer Gewichtsreduzierung.
Heute werden bis zu einem Nebenstromverhältnis von etwa 6,5 gemischte Abgassysteme
eingesetzt [52], danach überwiegen meist die Vorteile offener Systeme.
2.2. Quantitative Leistungsbewertung

Zur Beurteilung der Leistungsfähigkeit gemischter Abgassysteme werden im Folgenden
Abschnitt Leistungsparameter definiert. Mit ihnen soll die Effizienz des Düsensystems
für unterschiedliche Betriebszustände quantifiziert werden können.
2.2.1. Durchflusskoeffizient
Der Durchflusskoeffizient CD definiert sich aus dem Verhältnis von verlustbehafteten

Massenstrom ṁ aus der Rechnung oder Messung, zum theoretisch ideal erreichbaren
Massenstrom ṁid durch die Düse. Verluste werden z. B. durch die Wandrauigkeit, durch
die geometrische Form der Düse und durch die Strömungsqualität beeinflusst. Bedingt
durch die Verluste ist der reale Massenstrom durch die Düse geringer als bei der idealen,
isentropen Strömung. Die sich daraus ergebende Definition für CD kann z. B. in [53]
gefunden werden:
√
Tt
ṁ 1 pt
ṁ
CD = = . (2.1)
ṁid Athroat Qideal
Athroat ist die Querschnittfläche im Düsenhals und Q ist der reduzierte Massenstrom
oder auch die sogenannte Q-Funktion. Diese berechnet sich mit
√
Tt
ṁ pt
Q= . (2.2)
Athroat
Das ideale Q hingegen ist definiert als

v !
√
κ1 u κ−1
ṁideal pTt t p u 2κ p κ
Qideal = = t 1− , (2.3)
Athroat pt R(κ − 1) pt
11

pt /p ≤ pt /p∗ , p = p∞
mit ,
pt /p > pt /p∗ , p = p∗
wobei pt /p∗ das kritische Druckverhältnis, κ der Isentropenexponent und R die spezifi-
sche Gaskonstante sind. Das kritische Druckverhältnis ergibt sich zu [41]
κ−1
κ
pt κ−1
= 1+ . (2.4)
p∗ 2
Die Fallunterscheidung für das kritische Druckverhältnis wird vorgenommen um den

idealen Massenstrom zu korrigieren, je nachdem ob die Düse gesperrt ist oder nicht.
Diese Definition für CD is nur auf einstrahlige Düsensysteme anwendbar. Wie in [53]
nachzulesen ist, kann der Durchflusskoeffizient aber auch als Verhältnis der effektiven
zur geometrischen Düsenhalsfläche ausgedrückt werden. Die effektive Fläche (Aef f =
CD · A) kann unter Verwendung des gemessenen bzw. berechneten Massenstroms gemäß
Gleichung 2.1 berechnet werden. Im Falle von gemischten Strahlen ist es gängige Praxis,
den Durchflusskoeffizienten mit der Summe der effektiven Flächen beider Ströme zu
bilden:
 √ ! √ ! 
P Tt Tt
Aef f 1 ṁ pt
ṁ pt
CD = =  + . (2.5)
Athroat Athroat Qideal Qideal
f an core
2.2.2. Geschwindigkeitskoeffizient
Der Geschwindigkeitskoeffizient CV stellt per Definition das Verhältnis von tatsächli-

cher Düsenaustrittsgeschwindigkeit u9 zur idealen Düsenaustrittsgeschwindigkeit u9,id
dar [41]:
uexit
CV = . (2.6)
uexit,id.
Dieser Leistungsparameter repräsentiert den Effekt von viskosen Verlusten in der Grenz-
schicht der Düse und ist ein Gradmesser für die Effizienz des Düsensystems [53]. Der vis-
kose Einfluss ist abhängig von der Größe der umströmten Oberfläche, von der Machzahl
und der damit verbundenen Stoß-Grenzschicht-Wechselwirkungen in Wandnähe sowie
von der Reynoldszahl. Allerdings kann der Reynoldszahleinfluss in der Regel vernach-
lässigt werden. Danach bleibt als einziges die geometrische Abhängigkeit, da auch die
12
Machzahl längs des divergenten Düsenteils nur vom Flächenverhältnis Aexit /Athroat ab-
hängt . In Anlehnung an die eben erläuterte Definition soll CV hier allerdings über den
Bruttoschub der Düse berechnet werden. Dafür wird Gleichung (2.6) im Zähler und
Nenner jeweils um den ermittelten Massenstrom erweitert:
ṁ uexit FB,x
CV = = . (2.7)
ṁ uexit,id. ṁ uexit,id
Im Zähler steht nun der erreichte Bruttoschub der Düse FB,x in x-Richtung. Die ideale
Austrittsgeschwindigkeit kann gemäß der Formel von de Saint-Venant und Wantzel [54]
bestimmt werden:
v !
u 1−κ
u 2 κ Tt R pt κ
uexit,id =t 1− . (2.8)
κ−1 p∞
Der Bruttoschub wird mithilfe der Impulsbilanz über das maßgebliche Kontrollvolumen
der Düse ermittelt. Das Kontrollvolumen ist in Abbildung 2.1 schematisch dargestellt.
Die Impulsbilanz in x-Richtung stellt sich wie folgt auf:
Z Z Z
2 2
FB,x = ρu dA + ρu dA + (p − p∞ )dA
AF an ACore AF an
Z Z (2.9)
+ (p − p∞ )dA + σx dA,
ACore AW and
wobei σx die Dichte der Oberflächenkräfte in x-Richtung darstellt, welche viskose Kräfte
und Druckkräfte beinhaltet.
Analog zur Berechnung des Durchflusskoeffizienten für gemischte Abgassysteme in (2.5)
wird nun auch der Geschwindigkeitskoeffizient bestimmt:
FB,x
CV = . (2.10)
(ṁ uexit,id )F an + (ṁ uexit,id )Core
13
AF an
AN acelle Aexit
AM ixer
ACore
ABullet
Abbildung 2.1.: Schematischer Schnitt durch die Düse mit Mischer mit eingezeichnetem
Kontrollvolumen.
14
3. Numerische Verfahren
Für die Berechnung der Düsenströmung wurden zwei unterschiedliche CFD-Programme
verwendet. Zum einen wurde der kommerzielle Strömungslöser FLUENT 14.0 der Firma
ANSYS benutzt sowie der Open Source Code OpenFOAM. Beide Codes verwenden das
Finite-Volumen-Verfahren. Der Schwerpunkt des folgenden Kapitels liegt aber auf den in
OpenFOAM verwendeten und implementierten Verfahren. Für mehr Details bezüglich
FLUENT wird später auf die entsprechende Dokumentation der Software im FLUENT
Theory Guide [55], bzw. im FLUENT User’s Guide [56] verwiesen.
Der verwendete Löser in OpenFOAM basiert auf dem von Oliver Borm geschriebenen
Code [57, 58] für OpenFOAM-1.6-extend. Dieser diente als Basis und wurde für die in
dieser Arbeit durchgeführten Rechnungen zum Teil angepasst und erweitert.
3.1. Navier-Stokes Gleichungen

In der Strömungsmechanik spielen im wesentlichen drei Erhaltungsgleichungen eine wich-
tige Rolle: die Kontinuitätsgleichung, die Impulserhaltungsgleichung und die Energieer-
haltungsgleichung.
Die integrale Form der Kontinuitätsgleichung für ein beliebiges Kontrollvolumen lautet
[12]:
Z I
∂ρ ~
dV + ρ U · ~n dA = 0. (3.1)
V ∂t A
Diese Gleichung repräsentiert die Massenerhaltung. Der erste Term auf der linken Seite
stellt die unmittelbare Änderung der Masse innerhalb des Kontrollvolumens dar. Der
zweite Term gibt den Nettomassenfluss an, der aus dem Kontrollvolumen heraus fließt.
Die Impulserhaltungsgleichung innerhalb eines beliebigen Kontrollvolumens in integraler
Form ist definiert als [12]
Z I I Z
∂ ~
ρU dV + ~ ~
ρU U · ~n dA = − p ~n dA + ρ f~ext dV. (3.2)
V ∂t A A V
15
Numerische Verfahren
Der erste Term auf der linken Seite gibt die momentane Änderungsrate des Impulses
im Kontrollvolumen an, der zweite Term repräsentiert den Netto-Impulsfluss aus dem
Kontrollvolumen heraus. Auf der rechten Seite stehen zuerst der Netto-Druckkraft-Term
und letztlich noch der Term, der die externen Körperkräfte auf das Fluid berücksichtigt.
Zuletzt folgt die Energieerhaltungsgleichung gemäß [12]:
Z I
∂ ~ · ~n dA = Q̇ − Ẇ
(ρE) dV + ρH U (3.3)
V ∂t A
E und H stellen die Totalenergie und die Totalenthalpie dar. Die Terme im Einzelnen,
angefangen auf der linken Seite, sind die augenblickliche Änderungsrate der Totalenergie,
gefolgt vom Nettofluss der Totalenthalpie aus dem Kontrollvolumen. Auf der rechten
Seite stehen die Summen der zeitlichen Änderung des Wärmeflusses sowie die zeitliche
Änderung der Arbeit des Fluids im Kontrollvolumen. Die Totalenergie wird wie folgt
berechnet:
1 →2
E = e + U . (3.4)
2
Die Totalenthalpie H berechnet sich zu
1 →2 p
H = h + U = E + . (3.5)
2 ρ
Der statische Druck p lässt sich über die thermische Zustandsgleichung idealer Gase
bestimmen [59]:
p = ρRT. (3.6)
Der oben aufgeführte Satz an Erhaltungsgleichungen (3.1)-(3.3) lässt sich in kompakter

Notation wie folgt wiedergeben:
Z ~ I
∂W
dV + F~ − G
~ dA = 0. (3.7)
V ∂t A
~ beinhaltet den Satz der sogenannten konservativen Variablen

Der Spaltenvektor W
[60]:
16
 
ρ
~ =  ρU
W ~  (3.8)
ρE
F~ ist der Vektor konvektiver Flüsse:

 
ρ U ~ · ~n
 
~  ~ ~ 
F = ρU U · ~n + p~n  (3.9)
 
(ρE + p) U ~ · ~n
~ Auf den Inhalt dieses

Die restlichen Terme befinden sich im diffusiven Flussvektor G.
Vektors soll an späterer Stelle eingegangen werden.
3.2. Gemittelte Navier-Stokes Gleichungen

Wie bereits in der Einleitung erwähnt, ist für technische Strömungsprobleme die Direkte
Numerische Simulation (DNS) zu aufwändig. Die bei solchen Problemen auftretenden
Strömungen sind zumeist turbulent. Diese Turbulenz charakterisiert sich durch sehr
kleine Zeit- und Längenskalen, welche auch durch moderne Hochleistungscomputer per
DNS-Rechnung nicht vollständig aufzulösen sind. Folglich bedient man sich zumeist eines
statistischen Ansatzes. Eine beliebige Größe im Strömungsfeld wird dann durch einen
Mittel- und einen Schwankungswert beschrieben. Die Erhaltungsgleichungen werden da-
zu statistisch gemittelt. Am häufigsten findet die zeitliche Mittelung Anwendung. Diese
findet zumeist Gebrauch bei inkompressiblen Strömungen und wird auch als sogenann-
te Reynolds-Mittelung bezeichnet. Die Reynolds gemittelten Navier-Stokes Gleichungen
(RANS) führen dann zu zusätzlichen nicht geschlossenen Termen, welche durch geeignete
Turbulenzmodelle geschlossen werden müssen. Bei kompressiblen Strömungen muss zu-
sätzlich die Schwankung der Dichte des Fluids noch mit berücksichtigt werden. Bei einer
rein zeitlichen Mittelung, also der Reynolds Mittelung, würden nochmals weitere nicht
geschlossene Schwankungsterme hinzukommen, was die Komplexität des Schließungsan-
satzes enorm steigern würde [13]. Daher wird für kompressible Strömungen die Favre
Mittelung verwendet. Bei diesem Ansatz werden alle Zustandsgrößen, bis auf die Dichte
und den statischen Druck, zeitlich und dichtegewichtet gemittelt. Dadurch werden die
Schwankungsterme der Dichte eliminiert und man erhält somit die Favre-gemittelten
Navier-Stokes Gleichungen (FANS). Für die Herleitung der gemittelten Navier-Stokes
17
Gleichungen [13] über deren integrale Formulierung [12] sowie die dabei getroffenen An-
nahmen soll auf die entsprechende Literatur verwiesen werden.
Die Favre-gemittelten Navier-Stokes Gleichungen bilden die Grundlage für die in dieser
Arbeit durchgeführten Rechnungen. Die aus [13] bzw. [61] übernommenen Gleichungen
in differentieller Form schreiben sich in Indexschreibweise wie folgt:
∂ρ ∂ e
+ ρUi = 0, (3.10)
∂t ∂xi

ei
∂ ρU ∂ e e ∂P ∂ e
+ ρUj Ui = − + tji − ρu′′j u′′i , (3.11)
∂t ∂xj ∂xi ∂xj

e
∂ ρE ∂ e h e i ∂
+ Uj ρE + P = −qLj − ρu′′j h′′ + tji u′′ i − ρu′′ j k
∂t ∂xj ∂xj (3.12)
∂ he e i
+ Ui tij − ρui uj .
′′ ′′
∂xj
Der mittlere statische Druck berechnet sich analog zu (3.6):
P = ρRTe (3.13)
ρ und P sind die Reynolds-gemittelte Dichte und der Druck. Eine Tilde weist hingegen
auf Favre Mittelung hin. Der instantane Wert einer physikalischen Größe wird aufgeteilt
in einen gemittelten und einen fluktuierenden Anteil. So z. B. bezogen auf die Geschwin-
digkeit:
ei + u′′ i .
ui = U (3.14)
ei ist der gemittelte Geschwindigkeitsanteil, u′′ i ist ihr schwankender Anteil1 .

U
Für ein kalorisch perfektes Gas ergibt sich die Totalenthalpie dann zu
e = cp Te + 1 U
H ei U
ei , (3.15)
2
1
Eine doppelt gestrichene Größe (z. B. u′′i ) steht für den fluktuierenden Anteil der Favre Mittelung,
während eine einfach gestrichene Variable die Schwankungsgröße der Reynolds Mittelung ist (z. B.
u′i ).
18
und die Totalenergie zu
e = cv Te + 1 U
E ei U
ei . (3.16)
2
Der in Gleichung (3.12) auftretende Ausdruck k steht für die turbulente kinetische Ener-
gie (TKE). Sie ist definiert als
1
ρk = ρu′′i u′′i . (3.17)
2
Es sei an dieser Stelle nochmals erwähnt, dass die Favre Mittelung lediglich dazu dient,
Dichteschwankungen zu eliminieren um zusätzliche komplizierte Terme, welche Schlie-
ßung erfordern würden, zu vermeiden. Sie ist eine rein mathematische Vereinfachung,
keine physikalische [13]. Für mehr Einzelheiten bezüglich Mittelung sei auf die entspre-
chende Literatur verwiesen.
3.3. Turbulenzmodellierung
Wie im vorangegangenen Kapitel bereits angeführt, treten bei den gemittelten Navier-
Stokes Gleichungen zusätzliche Terme auf, für deren Schließung ergänzend Turbulenz-
modelle von Nöten sind. Die verwendeten Schließungsansätze sollen in den folgenden
Teilkapiteln kurz erläutert werden.
3.3.1. Allgemeine Schließungsansätze
Viskoser und turbulenter Spannungstensor
Der gängigste Ansatz zur Modellierung des turbulenten Spannungstensors für Null-, Ein-
und Zweigleichungsturbulenzmodelle gründet auf der Boussinesq-Annahme [13]. Die Idee
von Boussinesq basiert auf der Beobachtung, dass der Impulsaustausch in einer turbulen-
ten Strömung durch die Ausmischung durch große energiereiche Wirbel dominiert wird
und stellt einen Zusammenhang zwischen Spannungen und Deformationsgeschwindig-
keit her [61]. Die laminaren und turbulenten Spannungsanteile zusammengefasst, erhält
man
19
!
ek
1 ∂U 2
eij = e
σ tij − ρu′′i u′′j = e
tij + ρτij = 2 (µ + µT ) Sij − δij − ρkδij . (3.18)
3 ∂xk 3
µ ist die laminare dynamische Viskosität bzw. die molekulare Zähigkeit des Fluids. Im
Gegensatz dazu steht die turbulente Viskosität µT in keinem direkten Zusammenhang
mit der molekularen Zähigkeit, obwohl ihre Bezeichnung darauf hindeutet. Sie wird
auch als Austauschgröße oder als Wirbelviskosität bezeichnet [62] und beschreibt kei-
ne physikalische Charakteristik des Fluids. Vielmehr ist sie eine Funktion der lokalen
Strömungsbedingungen. Je nach Turbulenzmodell wird diese Größe unterschiedlich be-
rechnet. Sij ist die Deformationsgeschwindigkeit und k stellt wie besagt die turbulente
kinetische Energie dar. Auf die Berechnung von µT und k wird bei der Besprechung der
Turbulenzmodelle näher eingegangen.
Die laminare Viskosität µ ist eine stoffabängige Größe. Sie kann als konstanter, oder bes-
ser als temperaturabhängiger Wert formuliert werden. Eine sehr gute Näherung ist das
temperaturabhängige Sutherlandgesetz [63]. Es kann mit zwei oder drei Koeffizienten
beschrieben werden [56]. Mit drei Koeffizienten hat das Gesetz folgende Form:
3/2
T Tref + S
µ = µref (3.19)
Tref T +S
Die Größen µref und Tref sind konstante Referenzwerte, S ist die Sutherlandkonstante
mit der Dimension der Temperatur. Die Sutherlandformel aus (3.19) wurde hier als
temperaturabhängiges Polynom approximiert:
kg kg kg
µ = 2, 8474 · 10−6 + 5, 8034 · 10−8 · T − 2, 076 · 10−11 · T 2. (3.20)
ms msK msK 2
Der Verlauf von µ(T ) ist in Abbildung 3.1 grafisch dargestellt.
Molekularer und turbulenter Wärmestrom
Zur Berechnung des laminaren bzw. molekularen Wärmestroms wird zu allermeist das
Fourie’sche Wärmeleitgesetz verwendet [12]. Unter der Annahme eines kalorisch perfek-
ten Gases lautet dieses
∂T
qLj = −λ . (3.21)
∂xj
20
−5
x 10
4.5
3.5
3
µ [kg/ms]
2.5
1.5
1
200 300 400 500 600 700 800 900 1000
T [K]
Abbildung 3.1.: Polynom µ(T ).
Mit Einführung einer laminaren Prandtlzahl P r kann die Temperaturleitfähigkeit λ ab-

hängig von der Viskosität µ und der spezifischen Wärmekapazität cp berechnet werden:
cp µ
λ= . (3.22)
Pr
Die laminare Prandtlzahl ist eine Konstante und hat für Luft in der Regel einen Wert
von P r = 0, 72. Die Berechnung von λ soll hier wieder durch ein geeignetes, von der
Temperatur abhängiges Polynom, erfolgen:
W W W
λ = 3, 38045 · 10−3 + 7, 919012 · 10−5 2
· T − 1, 457318 · 10−8 3
· T 2 . (3.23)
mK mK mK
Das Polynom wurde aus Rechnungen von Rolls-Royce Deutschland übernommen. Der
Verlauf des Polynoms λ(T ) ist grafisch in Abbildung 3.2 wiedergegeben.
Zur Berechnung des turbulenten Wärmestroms qTj wird Gleichung (3.22) analog auf den
turbulenten Fall übertragen. Der turbulente Wärmestrom berechnet sich dann zu
21
0.07
0.06
0.05
λ [W/mK]
0.04
0.03
0.02
0.01
200 300 400 500 600 700 800 900 1000
T [K]
Abbildung 3.2.: Polynom λ(T ).
cp µT ∂ Te
qTj = ρu′′i h′′ = − (3.24)
P rT ∂xj
Hier bedient man sich der bereits oben angesprochenen turbulenten Viskosität, anstatt
der laminaren. Ebenso wird eine turbulente Prandtlzahl definiert. Aus gründen der Ein-
fachheit wird in der Regel eine konstante turbulente Prandtlzahl für das gesamte Rechen-
gebiet von P rT = 0, 9 angenommen. Tatsächlich variiert diese Zahl aber. So ist dieser
Wert in wandnahen Grenzschichten eine gute Näherung. In Freistrahlgebieten hingegen
sind Werte um P rT = 0, 5 angemessener [13]. Auf den Wert der turbulenten Prandtl-
zahl soll gegen später noch einmal ausführlicher eingegangen werden. Der in Gleichung
(3.24) gezeigte Schließungsansatz für den turbulenten Wärmestrom wird fortan als Ed-
dy Diffusivity Model (EDM) bezeichnet. Dieser Ansatz wird ebenfalls zu einem späteren
Zeitpunkt noch einmal genauer diskutiert. Der gesamte Wärmestrom kann dann mit der
laminaren und turbulenten Temperaturleitfähigkeit, α und αT , folgendermaßen ausge-
drückt werden:
22
∂e
h
qj = − (α + αT ) , (3.25)
∂xj
wobei α und αT wie folgt berechnet werden:
λ µT
α + αT = + . (3.26)
c p P rT
Molekulare Diffusion und turbulenter Transport
Übrig bleiben nun die noch nicht geschlossenen Ausdrücke tji u′′ i und ρu′′ j k in der Ener-
gieerhaltungsgleichung (3.12). Dabei handelt es sich um den molekular diffusiven und
turbulenten Transport von turbulenter kinetischer Energie. Der häufig verwendete An-
satz nach Wilcox [13] zur Schließung hierfür ist

µT ∂k
tji u′′ i − ρu′′ jk = µ+ , (3.27)
σk ∂xj
wobei σk ein Proportionalitätsfaktor ist und hier einen Wert von eins hat. Dieser Schlie-
ßungsterm wird von vielen Autoren häufig vernachlässigt [61], kann aber bei hypersoni-
schen Strömungen durchaus eine Rolle spielen [13].
Geschlossene Gleichungen
Mit den oben angeführten Schließungsansätzen (Gl. (3.18), (3.25) und (3.27)) kann nun
~ aus Gleichung (3.7) angegeben werden [58]:
auch der diffusive Flussvektor G
 
0

~ =  ~n ·
σ 
G  (3.28)
~ − ~n · ~q + ~n · µ +
~n · σ · U µT
∇k.
σk
Der geschlossene Satz gemittelter Navier-Stokes Gleichungen in differenzieller Form er-

gibt sich nun wie folgt:
∂ρ ∂ e
+ ρUi = 0, (3.29)
∂t ∂xi
23

ei
∂ ρU ∂ e e ∂P ∂ e
+ ρUj Ui = − + tji + ρτji , (3.30)
∂t ∂xj ∂xi ∂xj

e
∂ ρE ∂ e h e i ∂ µT ∂k
+ Uj ρE + P = −qj + µ +
∂t ∂xj ∂xj σk ∂xj (3.31)
h
∂ e e i
+ Ui tij + ρτij .
∂xj
3.3.2. Turbulenzmodelle
Hauptaufgabe der Turbulenzmodelle ist es, die bisher noch nicht quantifizierten Größen
wie die turbulente kinetische Energie k und die turbulente Viskosität µT zu berechnen.
Neben Eingleichungsturbulenzmodellen, wie das von Spalart und Allmaras [64], oder
dem Reynolds Stress Modell (RSM) [65, 66, 67] findet die Klasse der Zweigleichungs-
modelle bei der Berechnung technischer Strömungen die häufigste Anwendung. Letztere
sind vom Ansatz nicht so aufwändig wie das Reynolds Stress Modell, weil bei ihnen z. B.
isotrope Turbulenz vorausgesetzt wird, liefern aber für ein breites Spektrum an Anwen-
dungen eine hinreichende Genauigkeit bei gleichzeitig angemessenem Rechenaufwand.
Hier sollen die in dieser Arbeit verwendeten Zweigleichungsturbulenzmodelle vorgestellt
werden. Zum einen das k-ǫ Turbulenzmodell, anschließend das k-ω-SST Modell von Men-
ter.
k-ǫ Turbulenzmodell
Das wohl am weitesten verbreitete Turbulenzmodell ist das k-ǫ Modell. Die größten Bei-
träge hierzu lieferten vor allem Jones und Launder [68, 69], sowie Launder und Spalding
[70]. Aufgrund seiner weitreichenden Verbreitung und Anwendung wird das hier ver-
wendete Turbulenzmodell auch häufig als standard k-ǫ Modell (SKE) bezeichnet [71].
Ebenfalls gängige Vertreter der k-ǫ Gruppe sind neben dem SKE z. B. das realizable k-ǫ
Modell [72] und das RNG k-ǫ Modell [73].
Beim k-ǫ Turbulenzmodell wird je eine Transportgleichung für die turbulente kinetische
Energie sowie eine Gleichung für deren Dissipation ǫ gelöst. In differentieller Form lauten
beide:
24

∂ ∂ e ej
∂U ∂ µT ∂k
(ρk) + ρk Ui = ρτij − ρǫ + µ+ , (3.32)
∂t ∂xi ∂xi ∂xj σk ∂xj

∂ ∂ e ǫ ej
∂U ǫ2 ∂ µT ∂ǫ
(ρǫ) + ρǫUi = C1ǫ ρτij − C2ǫ ρ + µ+ . (3.33)
∂t ∂xi k ∂xi k ∂xj σǫ ∂xj
Für Strömungen mit höheren Machzahlen wird Turbulenz von der Kompressibilität der
Strömung durch die sogenannte „Ausdehnungsdissipation“ (engl. dilatation dissipation)
beeinflusst, was für gewöhnlich in den inkompressiblen Gleichungen vernachlässigt wird
[55]. Diese Art der Dissipation führt zu einer Abnahme der Ausdehnung für steigende
Machzahlen. Um diesen Kompressibilitätseffekt zu berücksichtigen wird hier der Ansatz
von Sarkar [74] verwendet, was Einfluss auf die Gleichung für die turbulente kinetische
Energie hat. Gemäß Sarkar ist die Ausdehnungsdissipation proportional zum Quadrat
der turbulenten Machzahl M aT . Gleichung (3.32) wurde daher um die Kompressibili-
tätskorrektur von Sarkar erweitert und lautet nun

∂ ∂ e ej
∂U 2
∂ µT ∂k
(ρk) + ρk Ui = ρτij − ρǫ 1 + M aT + µ+ . (3.34)
∂t ∂xi ∂xi ∂xj σk ∂xj
Die turbulente Machzahl ist eine Funktion der turbulenten kinetischen Energie und der
Schallgeschwindigkeit a:
r
2k
M aT = . (3.35)
a
Letztlich ist noch die turbulente Viskosität zu bestimmen. Die Berechnung von µT erfolgt
über die Verknüpfung der TKE und ihrer Dissipation:
k2
µT = ρCµ , (3.36)
ǫ
wobei Cµ eine Konstante ist.
Die Modellkonstanten für das k-ǫ Turbulenzmodell haben folgende Werte: C1ǫ = 1, 44,
C2ǫ = 1, 92, Cµ = 0, 09, σk = 1, 0 und σǫ = 1, 3.
25
k-ω-SST Turbulenzmodell
Das k-ω-Shear-Stress-Transport (SST) Turbulenzmodell von Menter [75, 76, 77] kombi-
niert das k-ω-Turbulenzmodell von Wilcox [13, 78] mit dem k-ǫ Modell in einem Tur-
bulenzmodell. Dabei versucht das k-ω-SST Modell die positiven Eigenschaften beider
Turbulenzmodelle zu vereinen. In wandnahen Bereichen bietet das k-ω Turbulenzmo-
dell bessere Eigenschaften gegenüber dem k-ǫ Modell. Im Grenzschichtnachlauf sowie in
Freistrahlen hingegen ist das k-ǫ dem k-ω Turbulenzmodell überlegen [61]. Über geeig-
nete Mischungsfunktionen wird beim SST Modell zwischen dem k-ω Modell und dem in
die k-ω Form transformierte k-ǫ Turbulenzmodell hin- und hergeschaltet. Eine weitere
Besonderheit des SST Turbulenzmodells ist auch die modifizierte Berechnung der tur-
bulenten Viskosität. Sie soll die Genauigkeit des Modells in Bereichen erhöhen, wo große
Druckgradienten und druckbedingte Grenzschichtablösungen auftreten.
Die Gleichungen für die turbulente kinetische Energie und die spezifische turbulente
Dissipationsrate ω lauten in differentieller Form

∂ ∂ e ej
∂U ∂ ∂k
(ρk) + ρk Ui = ρτij ∗
− β ρkω + (µ + σk µT ) , (3.37)
∂t ∂xi ∂xi ∂xj ∂xj

∂ ∂ e γ ej
∂U 2 ∂ ∂ω
(ρω) + ρω Ui = ρτij − βρω + (µ + σω µT )
∂t ∂xi νT ∂xi ∂xj ∂xj (3.38)
1 ∂k ∂ω
+ 2 (1 − F1 ) ρσω2 .
ω ∂xi ∂xi
Die Größen γ, β, σk und σω sind Konstanten. Je nachdem, in was für einem Bereich
sich die Strömung befindet, wird zwischen ihrem k-ω und k-ǫ Wert über die Verschnei-
dungsfunktion F1 hin- und hergeschaltet. Ihre Werte für die k-ω Formulierung sind
γ1 = 0, 5532, β1 = 0, 075, σk1 = 0, 85 und σω1 = 0, 5. Für die k-ǫ Formulierung lauten
sie γ2 = 0, 44, β2 = 0, 0828, σk2 = 1, 0 und σω2 = 0, 856. Für die Konstante β ∗ gilt
stets β ∗ = 0, 09. νT ist die kinematische turbulente Viskosität und kann anhand von µT
berechnet werden: νT = µT /ρ. Die dynamische turbulente Viskosität µT für das SST
Modell ist definiert zu
ρa1 k
µT = . (3.39)
max (a1 ω, SF2 )
a1 ist eine Konstante mit dem Wert 0, 31. S wird aus der Deformationsgeschwindigkeit
berechnet, F2 ist wiederum eine Verschneidungsfunktion. Für mehr Details sei wieder
26
auf die entsprechende Literatur [77] verwiesen.

Zur Berücksichtigung auftretender Kompressibilitätseffekte auf die Turbulenz in Berei-
chen hoher Machzahlen wurde auch im k-ω-SST Turbulenzmodell die Kompressibilitäts-
korrektur von Sarkar implementiert. Im Gegensatz zum k-ǫ Turbulenzmodell wirkt sich
die Korrektur auf beide Gleichungen des Turbulenzmodells aus [79]:
∂ ∂ e ej
∂U
(ρk) + ρk Ui = ρτij − β ∗ ρkω 1 + M a2T (1 − F1 )
∂t ∂xi ∂xi
(3.40)
∂ ∂k
+ (µ + σk µT ) ,
∂xj ∂xj
∂ ∂ e γ ej
∂U
(ρω) + ρω Ui = ρτij − βρω 2 + β ∗ ρω 2 M a2T (1 − F1 )
∂t ∂xi νT ∂xi
(3.41)
∂ ∂ω 1 ∂k ∂ω
+ (µ + σω µT ) + 2 (1 − F1 ) ρσω2 .
∂xj ∂xj ω ∂xi ∂xi
Die turbulente Machzahl M aT berechnet sich wieder gemäß Gleichung (3.35).
3.4. Strömungsberechnung und Lösungsverfahren

Historisch gesehen gibt es im Prinzip zwei unterschiedliche Arten die Navier-Stokes Glei-
chungen numerisch zu lösen [80]. Der druckbasierte Ansatz wurde zum Lösen inkom-
pressibler bzw. schwach kompressibler Strömungen mit geringen Machzahlen entwickelt.
Einige dieser Verfahren wurden dann immer weiter entwickelt, um auch kompressible
Strömungen mit mittleren Machzahlen lösen zu können. Alternativ hierzu wurden dich-
tebasierte Lösungsmethoden für transonische, stark kompressible Probleme entwickelt.
Auch diese Verfahren wurden dann später erweitert, um auch nahezu inkompressible
Strömungen mit dichtebasierten Lösern rechnen zu können. Bei beiden Verfahren, sowohl
beim druck- als auch beim dichtebasierten Ansatz, erhält man das Geschwindigkeitsfeld
aus der Impulsgleichung. Die Temperaturverteilung wird aus der Energiegleichung ge-
wonnen. Beim dichtebasierten Löser wird die Dichte aus der Lösung der Massenerhal-
tungsgleichung gewonnen. Mit Hilfe einer thermodynamischen Zustandsgleichung (z. B.
der idealen Gasgleichung) lässt sich dann daraus das Druckfeld bestimmen. Druckbasier-
te Verfahren lösen hingegen eine Druckgleichung. Aus dieser wird dann wiederum mit
der erwähnten Zustandsgleichung die Dichte berechnet.
27
Im Gegensatz zu den Rechnungen mit FLUENT, bei denen ein druckbasierter Löser
verwendet wurde, wurden für die Rechnungen mit OpenFOAM ein dichtebasiertes Ver-
fahren eingesetzt.
3.4.1. Numerische Flussberechnung
Für die Berechnung des konvektiven Flusses F~ kommt ein approximativer Riemannlöser
zum Einsatz. Zu den bekanntesten zählen unter anderem der Rusanov, der HLLE [81]
und HLLC [82, 83], die Reihe der AUSM [84, 85, 86, 87], sowie der Riemannlöser von Roe
[88, 89]. Letzterer wurde hier mit der Entropiekorrektur von Harten [90, 91] verwendet.
Das Verfahren von Roe gehört zur Klasse der Fluss-Differenzenverfahren. Der konvektive
Fluss F~ an der Zellfläche zweier benachbarter Zellen wird über deren Zustandsgrößen
durch Lösung des Riemann-Problems berechnet. Im Prinzip kann das Verfahren von
Roe auch als ein zentrales Differenzenverfahren mit einem zusätzlichen Dissipationsterm
angesehen werden. Aufgrund der hohen Genauigkeit, vor allem auch in Grenzschicht-
strömungen und bei der Auflösung von Stößen, findet die Berechnungsmethode von Roe
besonders häufig Anwendung [61].
Zur Erzielung einer Genauigkeit höherer Ordnung müssen bei diesen Verfahren die Zu-
standsgrößen an den Zellflächen rekonstruiert werden. Dies geschieht durch lineare Re-
konstruktion der Zustandsgrößen vom Zellzentrum aus mittels einer Taylerreihenexpan-
sion. Zur Vermeidung von nichtphysikalischen Überschwingern müssen die Gradienten
noch limitiert werden, sodass lokal keine neuen Minima oder Maxima entstehen. Die
Limitierung kann mithilfe sogenannter Slope-Limiter erfolgen [92, 93].
Da der verwendete Solver auf der Arbeit von Borm [57, 58] aufbaut und dessen Lö-
serarchitektur größtenteils übernommen wurde, soll für mehr Details auf dessen Arbeit
verwiesen werden.
3.4.2. Preconditioning
Zur Erhöhung der Genauigkeit der Lösung und für bessere Stabilität wurde der Rie-
mannlöser von Roe um ein Vorkonditionierungsverfahren bzw. ein sogenanntes Precon-
ditioning erweitert [94].
Prinzipiell sind dichtebasierte Zeitschrittverfahren nicht die erste Wahl, wenn es darum
geht ein Verfahren auszuwählen um inkompressible Strömungen mit sehr geringer Mach-
zahl zu lösen. Eines der Probleme ist, dass inkompressible Systeme in der Regel nicht voll
hyperbolisch sind, was zur Folge hat, dass der Druck nicht über die Zustandsgleichung
aktualisiert werden kann [80]. Außerdem resultiert aus der hohen Diskrepanz zwischen
der vorherrschenden Strömungsgeschwindigkeit und der Schallgeschwindigkeit in solchen
28
Strömungsfeldern eine schlechte Konditionierung des Systems. Die hohen Konditionszah-

len CN , definiert über das Verhältnis des maximalen zum minimalen Eigenwert
~|+a
|U
CN = , (3.42)
~|
|U
erhöhen die Steifigkeit des Gleichungssystems, was den Konvergenzprozess verlangsamt

[61, 80]. Ferner sind dichtebasierte Systeme nicht in der Lage den nötigen Betrag an
künstlicher Dissipation für kleiner werdende Machzahlen korrekt zu bestimmen. Dar-
unter leidet im allgemeinen die Genauigkeit für dichtebasierte Strömungslöser bei Strö-
mungsfeldern mit geringen Machzahlen [61].
Weiss und Smith [80] versprechen mit ihrem Preconditioning eine Methode, welche die
eben aufgeführten Probleme, die mit dichtebasierten Solvern einhergehen, zu lösen. Zum
einen wird bei dem Verfahren durch Induzierung einer zeitabhängigen Druckgleichung
eine künstliche Kompressibilität erzeugt und das System somit hyperbolisch gemacht.
Zum anderen werden die Eigenwerte so skaliert, dass sie alle dieselbe Größenordnung
besitzen. Dies führt dazu, dass große Konditionszahlen verschwinden und sich Werte um
eins annähern. Des Weiteren führt die geschickte Skalierung der Eigenwerte dazu, dass in
Bereichen mit niedrigen Geschwindigkeiten die Dissipation entscheidend erhöht wird.
Im Folgenden soll das Verfahren und die Vorgehensweise kurz vorgestellt werden.
Zur Generierung einer zeitabhängigen Druckgleichung wird das Gleichungssystem für die
konservativen Variablen aus (3.7) in ein System mit primitiven Variablen überführt.
~ Z ~ I
∂W ∂Q
dV + F~ − G
~ dA = 0, (3.43)
~
∂Q V ∂t A
mit
 
p
~ = U
Q ~ . (3.44)
T
Durch Multiplikation einer geeigneten Matrix K̄¯ erhält man dann folgende neue Matrix,
die als Faktor vor der zeitlichen Ableitung der primitiven Variablen steht:
29
 
ρ,p 0 0 0 ρ,T
0 ρ 0 0 0 
~  
¯ ∂W = 0 0  (3.45)
K̄
~ 0 0 ρ .
∂Q 0 0 0 ρ 0 
−1 0 0 0 ρcp
Die Ableitung der Dichte nach dem Druck ρ,p tritt als Faktor für die zeitliche Ableitung
des Drucks auf und bestimmt somit die Ausbreitungsgeschwindigkeit der Schallwellen im
System. Für ein ideales Gas gilt ρ,p = 1/RT = κ/a2 . Wobei für inkompressible Strömun-
gen ρ,p = 0 gilt, was wiederum konsistent ist mit der Annahme, dass sich Druckwellen in
einem solchen Medium mit praktisch infinitisimaler Geschwindigkeit ausbreiten. Ersetzt
man die Ableitung der Dichte nach dem Druck jedoch durch einen Term, der proportional
zur Inversen der lokalen Geschwindigkeit zum Quadrat ist, kann man so die Eigenwer-
te des Systems geschickt beeinflussen. Das System wird folglich vorkonditioniert indem
man Gleichung (3.45) durch die Precondition-Matrix ersetzt:
 
Θ 0 0 0 ρ,T
0 ρ 0 0 0 
 
Γnc =
0 0 ρ 0 0  . (3.46)
0 0 0 ρ 0 
−1 0 0 0 ρcp
Der Precondition-Parameter Θ wird angegeben mit

1 ρ,T 1 1
Θ= 2
− = 2
− . (3.47)
Ur ρcp Ur T cp
Die Referenzgeschwindigkeit Ur wird durch folgende Beziehung berechnet:

~ µ
Ur = min max |U |, , a · 10 −6
,a . (3.48)
ρ∆x
Die Berechnung dieser Größe ist letztlich maßgebend für die Bestimmung und Steue-
rung der Eigenwerte des Systems. ∆x ist eine charakteristische Länge, die für jede Zelle
berechnet wird. Die Berechnung dieser Größe wird in Abschnitt 3.6.1 erläutert.
Das nun vorkonditionierte System erhält man durch erneute Multiplikation des Glei-
¯ −1 :
chungssystems mit K̄
30
Z ~ I
∂Q
Γ dV + F~ − G
~ dA = 0, (3.49)
V ∂t A
mit
¯ −1 Γ .
Γ = K̄ (3.50)
nc
Letztlich erhält man das neue Gleichungssystem für die primitiven Variablen durch Mul-
tiplikation von (3.49) mit der inversen Matrix von Γ:
Z ~ I
∂Q
dV + Γ−1 F~ − G
~ dA = 0. (3.51)
V ∂t A
Bei der Verwendung des Verfahrens von Roe zur Flussberechnung ist darauf zu achten,
dass sich durch die Vorkonditionierung die Eigenwerte des Systems und somit auch die
Einträge in der Roe-Matrix ändern. Die neue Flussberechnung für das vorkonditionierte
System kann der entsprechenden Literatur entnommen werden [61, 80].
3.5. Randbedingungen
Im folgenden Abschnitt sollen die numerischen Randbedingungen kurz erläutert werden.
Zur besseren Übersicht sind die Ränder in Abbildung 3.3 grafisch dargestellt. Bei den
Randbedienungen wurde hauptsächlich auf bereits existierende Implementierungen zu-
rückgegriffen. Für den Fall, dass die Strömung sowohl sub- als auch supersonisch austritt,
wurde eine neue Randbedingung implementiert.
3.5.1. Einlass
Die Bestimmung des statischen Drucks erfolgt an den Einlässen über eine Neumann-
Randbedingung. So wird hier ein Druckgradient von Null vorgegeben. Die statische
Temperatur berechnet sich über eine Isentropenbeziehung. Durch die zuvor definierte
Totaltemperatur und den Totaldruck lässt sich die die Temperatur wie folgt ermitteln:
κ−1
p κ
T = Tt . (3.52)
pt
31
Abbildung 3.3.: Darstellung der Ränder des Rechengebiets.
Der Geschwindigkeitsvektor wird unter Vorgabe der Einlassrichtungsvektors und der

Totaltemperatur bzw. statischen Temperatur bestimmt:
q
~ =ϕ
U ~ 2cp (Tt − T ). (3.53)
Die Randbedingungen für die Turbulenzmodelle berechnen sich zum Teil aus Größen aus
dem Strömungsfeld. Die Bestimmung der turbulenten kinetischen Energie am Eintritt
erfolgt durch
3 ~ 2
k = I 2 U . (3.54)
2
I ist die turbulente Intensität und muss vorgegeben werden.
Die Berechnung der turbulenten Dissipationsrate ǫ erfolgt unter Vorgabe eines zuvor
spezifizierten turbulenten Längenmaßes l:
3
k2
3
ǫ = Cµ .
4
(3.55)
l
Ähnlich wie in Gleichung (3.55) wird auch die spezifische turbulente Dissipationsrate ω
unter Vorgabe des turbulenten Längenmaßes am Einlass bestimmt:
1
k2
ω= 1 . (3.56)
Cµ4 l
32
Cµ hat den bereits bekannten Wert von 0,09.
3.5.2. Auslass
Am Auslass des Fernfelds werden die Randbedingungen für Temperatur, Geschwindig-

keit sowie für die turbulenten Größen als Nullgradienten-Randbedingungen formuliert.
Um eine definierte Strömungsgeschwindigkeit zu erhalten, wird der statische Druck am
Rand als konstant definiert falls die Machzahl an der Austrittsebene M a ≤ 1 beträgt.
Wenn die Machzahl am Austritt größer ist als eins, wird vom Löser automatisch auf eine
Nullgradienten-Randbedingung gewechselt.
3.5.3. Fernfeld
Am oberen Fernfeldrand (siehe Abbildung 3.3) wurde wieder für alle Größen bis auf
den statischen Druck eine Nullgradienten-Randbedingung gewählt. Der statische Druck
wird als konstant definiert. Je nach Strömungszustand ist am Fernfeldrand ein Ein- oder
Ausströmen des Fluids möglich.
3.5.4. Wände
An den Wänden gilt für die Geschwindigkeit die Haftbedingung. Aus diesem Grund
wird die Fluidgeschwindigkeit dort zu Null gesetzt. Für den statischen Druck und die
statische Temperatur sowie die turbulente kinetische Energie wird eine Nullgradienten-
Randbedingung verwendet. Da die verwendeten Rechennetze so konzipiert wurden, dass
sich an den Wänden y + -Werte von ca. 30 ergeben, können für die restlichen Turbul-
nezgrößen sogenannte Wandfunktionen eingesetzt werden. Diese speziellen Funktionen
werden in OpenFOAM für den wandnahen Bereich ab ca. y + > 11 verwendet. Falls
<
einzelne Zellen doch kleinere Werte aufweisen, also y + ∼ 11, wird die turbulente Visko-
sität µT sowie der Produktionsterm des Turbulenzmodells gleich Null gesetzt. Für mehr
Details zur Wandbehandlung soll auf die jeweilige Stelle im Quelltext in OpenFOAM
verwiesen werden.
33
3.5.5. Seitenränder
Um den Rechenaufwand möglichst gering zu halten wurde die Düse nicht als volle 360◦ -
Geometrie gerechnet, sondern als 45◦ -Teilstück ausgeführt (siehe Abbildung 3.4). Die
Seitenränder wurden entweder als symmetrisch oder periodisch definiert.
Abbildung 3.4.: 45◦ -Stück der Düsengeometrie mit Darstellung der Seitenränder. Blick
von vorne in die Düse.
34
3.6. Zeitschrittverfahren
Hier soll das verwendete Zeitschrittverfahren kurz erläutert werden. Zur Konvergenz-
beschleunigung kommt zum einen ein lokales Zeitschrittverfahren zum Einsatz, zum
anderen wird die Runge-Kutta Zeitschrittmethode verwendet.
3.6.1. Lokales Zeitschrittverfahren
Explizite Zeitschrittverfahren sind nur bis zu einer gewissen Zeitschrittgröße stabil. Des-
halb richtet sich bei zeitechten Simulationen der verwendete Zeitschritt in der Regel nach
dem kleinsten Zeitschritt, bei dem die Rechnung gerade noch stabil ist. Für stationäre
Lösungen kann die Konvergenz beschleunigt werden, indem statt eines globalen Schritts
für jede Zelle ein lokaler Zeitschritt berechnet wird. Dieser für jede Zelle individuelle
Zeitschritt ist so zu wählen, dass dieser gerade noch das Stabilitätskriterium erfüllt. Die
Stabilitätsbedingung ist die sogenannte Courant-Friedrichs-Lewy Bedingung, oder kurz
CFL-Bedingung [95, 96].
Zur Berechnung des lokalen Zeitschritts findet man in der Literatur zahlreiche Varian-
ten (z. B. [61, 80, 97, 98]). Für die hier durchgeführte Berechnung wurde der Zeitschritt
analog zum Verfahren von Weiss & Smith [80] ermittelt. Für die Bestimmung des ma-
ximalen Zeitschritts ∆t ist eine reibungsfreie, ∆tinv , und eine viskose Zeitschrittweite,
∆tvis , zu bestimmen. Der reibungsfreie Zeitschritt berechnet sich mit
∆x
∆tinv = . (3.57)
~|+a
|U
Der viskose Zeitschritt ist gegeben durch
(∆x)2 ρ
∆tvis = . (3.58)
µ + µT
Die Größe ∆x in (3.57) und (3.58) ist hier eine charakteristische Zelllänge. Für jede
Zellfläche wird der Abstandsvektor vom Zellmittelpunkt zum Mittelpunkt der Zellfläche
bestimmt. Die minimale Abstandslänge in jeder Zelle dient schließlich als charakteristi-
sche Länge:

∆x = min f~n − P~ · ~nn . (3.59)
35
Die physikalische Zeitschrittweite erfolgt nun über die Bestimmung des minimalen Zeit-
schritts aus Gleichung (3.57) und (3.58) in jeder Zelle:
∆t = CF L · min [∆tinv , ∆tvis ] . (3.60)
Mit der Vorgabe der CF L-Zahl kann die Zeitschrittweite während der Rechnung regu-
liert werden.
Ein lokales Zeitschrittverfahren macht nur Sinn, wenn die Zellen bzw. die berechneten
Zeitschrittweiten innerhalb des Rechengebiets unterschiedlich groß sind. Je größer die
Unterschiede zwischen den einzelnen Zellen, desto größer fällt folglich die Konvergenz-
beschleunigung aus.
3.6.2. Runge-Kutta Zeitschrittverfahren
Eine weitere, sehr häufig verwendete Methode zur Verbesserung der Konvergenz ist das
Runge-Kutta Verfahren [95]. Dieses mehrstufige Verfahren unterteilt den Zeitschritt in
mehrere Unterschritte. Der Ablauf eines m-stufigen Runge-Kutta Verfahrens ist wie
folgt:
~ (0) = Q
Q ~ (i)
~ (1) = Q
Q ~ (0) − β1 ∆t ~ (0)
R
~ (2) = Q
Q ~ (0) − β2 ∆t ~ (1)
R (3.61)
..
.
~ (i+1) = Q
Q ~ (m) = Q
~ (0) − βm ∆t R
~ (m−1) .
~ (i) ist die Lösung zum vorherigen Zeitschritt, Q

Q ~ (i+1) ist die neue Lösung. Der Vorteil
dieser Methode ist, dass abhängig von der Anzahl der Zwischenschritte und der Wahl
der Runge-Kutta Koeffizienten β die Stabilität, Genauigkeit und der maximal mögliche
Zeitschritt erhöht werden können [61]. Zwar sind durch das Runge-Kutta Verfahren
prinzipiell höhere CFL-Zahlen möglich, es sei jedoch noch angemerkt, dass eine größere
Anzahl an Unterstufen auch einen erhöhten Rechenaufwand zur Folge hat.
36
3.7. Gittergenerierung
Einen großen Teil der Arbeit nimmt die Gittergenerierung in Anspruch. Besonders die
komplexe Geometrie des gescarften Mischers erhöht den Arbeitsaufwand und den Schwie-
rigkeitsgrad der Vernetzung. Für die Gittererzeugung wurde die kommerzielle Software
ICEM CFD 14.0 der Firma ANSYS benutzt. Mittels dieser Software ist es möglich das
gesamte Rechengebiet inklusive Mischer blockstrukturiert zu vernetzen.
3.7.1. Geometrieerzeugung
Im Rahmen dieser Arbeit wurden zwei Geometrien vernetzt, welche beide von Rolls-
Royce Deutschland zur Verfügung gestellt wurden. Zum einen kommt die Düsengeome-
trie aus dem Forschungsprojekt OPTITHECK (OT) zum Einsatz, zum anderen wurde
eine reale Rolls-Royce (RR) Düsengeometrie bereitgestellt. Die Geometrien wurden bei-
de als als STP-Datei bereitgestellt. Dieses CAD-File kann in ICEM CFD geladen und
geöffnet werden. Abbildung 3.5 zeigt eine Visualisierung einer solchen STP-Datei, bei der
schon die Aufteilung der einzelnen Geometriesegmente in Parts (farblich hervorgehoben)
in ICEM CFD vorgenommen wurde.
Abbildung 3.5.: Darstellung der einzelnen Parts der OT-Düsengeometrie in ICEM CFD.
Für die Simulation der Düsenströmung mit einem äußerem Co-Flow muss im Vernet-
zungsprogramm auch die Größe des Fernfelds definiert werden. Der Fokus dieser Arbeit
liegt speziell auf der Vermischung innerhalb der Düse und nahe dem Düsenauslass, eben-
so wie auf der Bestimmung der Leistungsdaten. Daher kann das Rechengebiet um die
Düse herum klein gehalten werden um Zellen und somit Rechenzeit einzusparen. Das
Fernfeld wurde so gewählt, dass sich der numerische Auslass je ca. 1,5 Düsendurchmesser
hinter dem Düsenauslass befindet und der obere Fernfeldrand in radialer Richtung in
beiden Düsengeometrien 2,5 Düsendurchmesser beträgt (siehe Abbildung 3.6).
37
Abbildung 3.6.: Darstellung der OT-Geometrie mit Fernfeld im Vernetzungsprogramm.
3.7.2. Vernetzung
Bisher wurde keine hier bekannte strukturierte Vernetzung einer Triebwerksdüse mit in-
ternem komplex gescarften Mischer realisiert. Grund dafür sind die stark mäanderförmi-
gen Flächen des Mischers, die z. B. bei Vernetzungstools wie Gambit zu stark verzerrten
Zellen führen [99]. Mit ICEM CFD ist es aber möglich ein Blocking so zu erstellen, dass
sich das gesamte Rechengebiet blockstrukturiert vernetzen lässt. Möglich macht dies die
Verwendung von sogenannten O-Grids um den Mischer herum. Dadurch lässt sich das
Netz gut um den Mischer ziehen, ohne Zellen schlechter Qualität zu generieren (das
Blocking und das Netz im Bereich des Mischers sind jeweils in Abbildung 3.7 und 3.8
dargestellt). Die Geometrie des Mischers ist lediglich jeweils als Fläche ausgeführt, sie
besitzt also keine physikalische Dicke und damit keine Hinterkante. Diese Gegebenheit
erleichtert die strukturierte Vernetzung in ICEM CFD zusätzlich.
Erwähnenswert ist auch die Behandlung des Netzes nahe der Strahlachse. Die Gitter-
zellen wurden nicht bis zur Achse geführt. Stattdessen wurde in Achsennähe eine Art
Welle „ausgestanzt“. Die Ausführung ist anschaulich in Abbildung 3.9 dargestellt. Da-
durch wird verhindert, dass die Gitterzellen auf der Strahlachse singulär in einem Punkt
zusammenlaufen, was zu Problemen beim Strömungslöser führen kann. Der Durchmesser
dieser Aussparung beträgt ca. 1% des Düsendurchmessers. Da auch der Einlasskanal des
Kernstroms bis zur Achse geht, wurde auch dort die Geometrie entsprechend vernetzt. Es
wird angenommen, dass diese Art der Vernetzung keinen oder nur einen äußerst geringen
Einfluss auf die Strömungssimulation hat. Nachteil dieses Verfahrens ist allerdings, dass
dadurch sehr lange und schmale Stabzellen in Achsennähe entstehen. Dies führt unter
anderem dazu, dass sich dadurch sehr kleine Zeitschritte ergeben (siehe Kapitel 3.6.1),
38
Abbildung 3.7.: Verlauf des Blockings um den Mischer.
Abbildung 3.8.: Visualisierung des Gitters im Bereich des Mischers.
39
was die Konvergenz der Simulation verlangsamt. Eine Alternative zu diesem Ansatz für
folgende Netze wäre stattdessen ein O-Grid an den unteren Teil des Plugs anzuflanschen.
Dadurch würde eine prismenförmige Zelle direkt auf der Strahlachse liegen. So wäre das
Rechengebiet einerseits bis zur Strahlachse vernetzt, andererseits ergeben sich durch die
größeren Zellen entlang der Strahlachse wiederum größere Zeitschritte.
Abbildung 3.9.: Gitternetz auf dem Plug mit Detailansicht in Achsennähe.
Im wandnahen Bereich wurde darauf geachtet, dass die Zellhöhen so gewählt wurden,
dass sich bei den Rechnungen y + -Werte von ca. 30 ergeben. Dadurch können bei der spä-
teren Simulation Wandfunktionen eingesetzt werden, wodurch die Grenzschicht nicht voll
aufgelöst werden muss. D.h. der Abstand der ersten Zelle zur Wand wird so gewählt,
dass die viskose Unterschicht und die Pufferschicht in der Grenzschicht von den Wand-
funktionen überbrückt werden [71]. Dieser Ansatz erhöht deutlich das Potential, Zellen
im Rechengebiet einzusparen.
3.7.3. Fertigstellung
Nachdem das Gitter als blockstrukturiertes Netz in ICEM CFD erzeugt wurde, wird es
in ein unstrukturiertes Netz konvertiert. Dieser Vorgang erfolgt ebenfalls in ICEM CFD.
Der Grund hierfür ist, dass sowohl OpenFOAM als auch FLUENT nur Netze verarbei-
ten können, welche als unstrukturierte Gitter gespeichert werden. Anschließend kann das
40
unstrukturierte Netz in ICEM CFD in eine FLUENT-Mesh-Datei (.msh) umgewandelt

werden.
Bevor das Netz in das OpenFOAM Format konvertiert wird, muss es in diesem Falle
mindestens einmal in FLUENT geöffnet werden. Der Grund liegt in der bereits erwähn-
ten Tatsache, dass der Mischer nur als Flächenelement ausgeführt ist. Beim Öffnen der
.msh-Datei erzeugt FLUENT aus der bisher einteiligen Mischerfläche eine zusätzliche
Shadow-Fläche. Der Mischer besteht nun aus zwei Flächen die quasi identisch sind und
aufeinander liegen. Allerdings wird eine Fläche der Oberseite (z. B. Fanstromseite), die
andere Fläche der Unterseite (z. B. Kernstromseite) zugeordnet. Nach diesem Vorgang
kann das Netz als FLUENT-Case-Datei (.cas) gespeichert werden, um es anschließend
mit der OpenFOAM Anwendung fluent3DMeshToFoam in das OpenFOAM Datenformat
umzuwandeln.
41
4. Simulation der Düsenströmung
mit OpenFOAM
4.1. Simulationsgrundlagen
4.1.1. Berechnungsdetails
Für die Berechnungen wurden zwei Gitternetze eines jeweils 45◦ -Tortenstücks der Düse
erstellt. Das Netz für die OT-Geometrie wurde mit ca. 1,25 Mio. Zellen diskretisiert. Die
Gittergröße der RR-Düse beträgt ca. 1,47 Mio. Zellen. Beide Geometrien wurden mit
einem Faktor von ca. 1/5 gegenüber ihrem 1:1 Modell skaliert. Die Rechennetze sind in
den Abbildungen 4.1 und 4.2 dargestellt.
Die OT-Düse wurde für die Leistungsrechnungen und zur Untersuchung des Einflusses
der unterschiedlichen Modelle auf die Strahlmischung verwendet. Die Geometrie der
RR-Düse wurde erst sehr spät gegen Ende dieser Arbeit von Rolls-Royce Deutschland
bereitgestellt. Sie diente für Vergleichs- und Kalibrierungsrechnungen mit experimentell
ermittelten Totaltemperatur- und Totaldruckprofilen.
Die Berechnung des Strömungsfelds erfolgte stets mit einer Genauigkeit von zweiter Ord-
nung im Raum. Die Turbulenzmodelle wurden mit erster oder zweiter Ordnung diskre-
tisiert. Um mögliche Oszillationen bei der Berechnung zweiter Ordnung vorzubeugen,
wurde für das Strömungsfeld der von Borm implementierte Slope-Limiter von Venka-
takrishnan [93] verwendet. Ferner wurden zusätzlich nochmals alle Gradienten mittels
einem OpenFOAM-eigenen Zell-zu-Zell-Limiter limitiert, um das Auftreten neuer loka-
ler Minima oder Maxima bei zweiter Ordnung zu vermeiden. Zur Gradientenberechnung
wurde das Verfahren von Gauss verwendet [100].
Für die Runge-Kutta Zeitschrittmethode wurde ein dreistufiges Verfahren gewählt. Die
verwendeten Koeffizienten 0, 1481, 0, 4 und 1, 0 aus [61] stellten sich als besonders effek-
tiv bezüglich des Konvergenzverhaltens heraus.
Als Strömungsmedium diente Luft als ideales Gas mit einer konstanten spezifischen
Wärmekapazität von cp = 1007 J/kgK.
42
Simulation der Düsenströmung mit OpenFOAM
Abbildung 4.1.: 45◦ -Gitternetz der OT-Düsengeometrie. Links: Schnitt durch das ge-
samte Rechengebiet. Rechts: Darstellung der Geometrie mit Seiten-
randfläche.
Abbildung 4.2.: 45◦ -Gitternetz der RR-Düsengeometrie. Links: Schnitt durch das ge-
samte Rechengebiet. Rechts: Darstellung der Geometrie mit Seiten-
randfläche.
43
4.1.2. Randbedingungen
Leistungsrechnungen
Für die Leistungsrechnungen wurden unterschiedliche Betriebszuständen simuliert. Die

Arbeitspunkte der Düse definieren sich im Prinzip über die Totaltemperatur und den
Totaldruck am jeweiligen Einlass sowie dem Umgebungszustand im Fernfeld. Um die
Betriebszustände später parametrisch angeben zu können werden folgende Variablen
definiert:
• Nebenstromdüsendruckverhältnis, CNPR=pt,f an /p∞

• Nebenstrom- zu Kernstromtotaldruckverhältnis, PS = pt,f an /pt,core
• Kern- zu Nebenstromtotaltemperaturverhältnis, TR = Tt,core /Tt,f an
Mit diesen Größen und nachfolgender Nomenklatur lassen sich jetzt sämtliche durch-
geführten Betriebszustände beschreiben: <CNPR>-<PS>-<TR>-<Turbulenzmodell>.
Also z. B. 1.6-1.1-2.4-kOmegaSST steht für CNPR=1,6, PS=1,1, TR=2,4, durchgeführt
mit dem k-ω-SST Turbulenzmodell.
Die Totaltemperatur an den Einlässen vom Nebenstrom und dem Fernfeld wurde immer
auf 288,15 K gesetzt. Der Umgebungsdruck hat einen Totalwert von pt,∞ = 101325 Pa.
Die Außenströmung wurde durch Absenken des statischen Drucks auf p∞ = 101148 Pa
so eingestellt, dass sich eine Außenströmung mit ca. M a∞ = 0, 05 ergibt. Mit diesen
Angaben und den Werten der nachfolgenden Tabelle 4.1, welche die simulierten Betrieb-
spunkte beinhaltet, können die restlichen Totalzustände an den Einlässen berechnet
werden.
CNPR PS TR
1.0
1.5
1.6 1.1 2.0
2.4
1.0
1.5
2.6 1.1 2.0
2.4
2.55
Tabelle 4.1.: Simulierte Betriebspunkte.
44
Als Randbedingungen für die Turbulenzmodelle wurde je von einem Turbulenzgrad am

Einlass von Fan- und Kernstrom von 5% ausgegangen. Am Fernfeld von 0, 1%. Für das
turbulente Längenmaß an Fan- und Kerneinlass wurde 5 · 10−5 m gewählt und für das
Fernfeld ein Wert von 0, 0015 m.
Die Seitenränder des Rechengebiets wurden als Symmetrieebenen definiert.
Sideline-Fall
Der Vergleich von experimentellen Totaltemperatur- und Totaldruckmessungen zu Kali-

brierungszwecken wurde für einen sogenannten Sideline-Fall berechnet. Dieser simuliert
einen Überflug des Flugzeugs in geringer Höhe und wird in der Regel für Lärmemissi-
onsmessungen gewählt.
Für die Temperaturmessungen im Labor wurde für diesen Fall kein Co-Flow eingestellt,
was in der Simulation berücksichtigt wurde. Allerdings wurde zugunsten einer besseren
numerischen Konvergenz eine geringe Machzahl für die Umgebungsströmung von ca.
M a∞ = 0, 05 eingestellt.
Die Randbedingungen für die Turbulenzmodelle wurden wieder wie bei den Leistungs-
rechnungen gesetzt. Die Seitenränder wurden als periodisch definiert.
4.1.3. Simulationsablauf
Um die Erzeugung einer Simulationslösung mit OpenFOAM zu beschleunigen, wurde

für jede Rechnung eine Startlösung mit FLUENT erzeugt. Aufgrund der in FLUENT
zahlreich implementierten Methoden (z. B. Löser mit Geschwindigkeits- und Druckkop-
pelung, Multigrid [56]) kann die Konvergenz hin zum stationären Zustand der Simula-
tion gegenüber OpenFOAM deutlich beschleunigt werden. Für die Transformation der
Lösung aus FLUENT stellt OpenFOAM die Anwendung fluentDataToFoam bereit. Al-
lerdings werden mit dieser Applikation nicht alle benötigten Felder übernommen. Um
alle relevanten Größen aus FLUENT in das OpenFOAM-Format zu importieren wurde
deshalb die Anwendung fluentDataToFoam modifiziert bzw. erweitert.
Nach erfolgtem Datenimport aus FLUENT und dem Setzen der Randbedingungen kann
die Simulation gestartet werden. Um ein stabileres Anrechnen mit der Startlösung zu
realisieren, wurden die Rechnungen immer mit kleineren CFL-Zahlen begonnen und die-
se dann langsam gesteigert. Als guter Startwert stellte sich CFL=0,3 heraus. Innerhalb
der nächsten ca. 2000 bis 3000 Iterationen wurde die CFL-Zahl dann langsam auf Ma-
ximalwerte von bis zu 1,5 gesteigert, meistens wurde jedoch mit Werten von CFL=1,3
gerechnet.
45
4.1.4. Konvergenzkriterien
Zur Überprüfung der Rechnung auf Konvergenz wurde auf die absoluten Residuen der
Erhaltungsgleichungen geachtet. Ebenso waren vor allem die Massenströme an den Ein-
lässen von Kern- und Fanstrom für die Konvergenzbewertung relevant. Während bei den
Residuen bei fortgeschrittenem Simulationsablauf praktisch keine Änderungen mehr zu
beobachten waren, dauert es erfahrungsgemäß noch einige Zeit bis sich auch für die
Massenströme an den Inlets von Kern- und Nebenstrom konstante Werte einstellten.
Zahlreiche Tests und Überprüfungen haben ergeben, dass dieser Zustand bei den ver-
wendeten Gittern nach ca. 40000 Iterationen erreicht wird. Längeres Rechnen führt dann
zu keiner Veränderung mehr. Sind die Massenströme konstant, ändert sich weder das
Strömungsfeld, noch ändern sich dann die Leistungskoeffizienten CD und CV bei weite-
rem Rechenfortschritt. Auch die Gesamtmassenstrombilanz über alle Ein- und Auslässe
ist dann ausgeglichen. Die Differenz zwischen ein- und ausströmender Masse hat je nach
Rechnung eine Größenordnung von ca. 10−1 bis 10−4 %.
4.1.5. Postprocessing
Als Postprocessing-Software bzw. zur grafischen Auswertung der Rechnungen wurde das
Open Source Programm ParaView, Version 3.14 und 3.98, verwendet. Dazu wurden die
Daten in das von ParaView lesbare VTK-Format mittels der Applikation foamToVTK
umgewandelt. Für die quantitative Auswertung der Daten, z. B. zur Berechnung der
Leistungsparameter aus Kapitel 2.2, wurden teils eigene Skripte geschrieben, teilweise
konnten Skripte und Anwendung der OpenFOAM Software verwendet werden. Für die
Ausgabe der Felder, welche zur Berechnung der Leistungsparameter nötig waren wurden
hauptsächlich neue Skripte geschrieben oder existierende modifiziert. Die Integration
konnte dann mit einer Applikation von OpenFOAM durchgeführt werden.
46
4.2. Simulation der Düsenströmung mit

Temperaturkorrekturmodellen
4.2.1. Prinzip der Temperaturkorrektur
In der Einleitung wurde bereits auf die Defizite gängiger Turbulenzmodelle bei der Si-
mulation von heißen Abgasstrahlen hingewiesen. Aus diesem Grund gibt es zahlreiche
Modifikationen für die weitverbreiteten Zweigleichungsturbulenzmodelle. Die Änderun-
gen reichen von modifizierten Modellkonstanten [101] bis hin zu größeren Eingriffe in
der Modellierung [74, 102, 103, 104]. Eine kleine Übersicht über die Modelle geben [105]
und vor allem [106].
Korrekturmodelle die direkt auf die Verbesserung hinsichtlich Anisothermie abzielen lie-
fern Abdol-Hamid [18] sowie Tam und Ganesan [22]. Beide nehmen das k-ǫ Turbulenz-
modell als Basis und verstärken die Ausmischung zwischen Heiß- und Kaltstrahl indem
in der Mischungsscherschicht die turbulente Viskosität erhöht wird. Tam und Ganesan
begründen die verstärkte Mischung in diesem Bereich mit einer in der Realität verstärk-
ten Kelvin-Helmholtz Instabilität in der Scherschicht. Einen Beweis für diesen Effekt
erbringen sie anhand einer Stabilitätsanalyse. Das Modell von Tam und Ganesan ist so
aufgebaut, dass an Orten hoher vorherrschender Dichte- und Geschwindigkeitsgradien-
ten zum Standardwert der turbulenten Viskosität (vgl. Gleichung (3.36)) ein zusätzlicher
Anteil addiert wird. Neben geänderten Koeffizienten ändern Tam und Ganesan auch den
Gleichungssatz, in dem die Korrektur von Pope [103] berücksichtigt ist. Abdol-Hamid
hingegen erklärt die von den Standardgleichungsmodellen unterschätze Ausmischung
zwischen heißen und kalten Strahlen mit der Abwesenheit von modellierter Turbulenz
thermodynamischer Größen. Die Eliminierung von Dichteschwankungen durch die Fa-
vre Mittelung vereinfacht zwar das Schließungsproblem und liefert hinreichend genaue
Ergebnisse für viele Anwendungen, vernachlässigt aber den Einfluss von Dichte- bzw.
Temperaturturbulenz bei anisothermen Strömungen. Abdol-Hamid erklärt, dass große
Fluktuationen einer thermodynamischen Variable (Druck oder Temperatur) auch di-
rekten Einfluss auf die Dichte und somit auf die Ausdehnung des Fluids haben. Diese
Dichteschwankungen und daraus folgende Ausdehnungsschwankungen beeinflussen dann
die Entstehung und Ausbreitung von Turbulenz, die in der Realität zu stärkerer Aus-
mischung anisothermer Strömungen führt. Im Gegensatz zu Tam und Ganesan benutzt
Abdol-Hamid das standard k-ǫ Modell und modifiziert lediglich die turbulente Viskosität.
Auch hier wird zur turbulenten Standardviskosität ein zusätzlicher Anteil addiert, der
bei Abdol-Hamid aber abhängig vom Totaltemperaturgradienten ist. Diese Form zur
Modellierung von Temperaturturbulenz deckt sich auch mit der Aussage von Béguier
[107]. Dieser erklärt, dass analog zu den Geschwindigkeitsfluktuationen, welche durch
47
den mittleren Geschwindigkeitsgradienten verursacht werden, Temperaturfluktuationen

durch den mittleren Temperaturgradient generiert werden.
Als favorisiertes Modell für die Implementierung in OpenFOAM wurde aufgrund der
größeren Universaltauglichkeit das Temperaturkorrekturmodell von Abdol-Hamid aus-
gewählt. Wie bereits erwähnt, benutzt Abdol-Hamid das standard k-ǫ Turbulenzmodell,
während Tam und Ganesan auf eine modifizierte Variante von selbigem von Theis und
Tam zurückgreifen [102]. Der dort beinhaltete Vortex-Stretching Term wurde ursprüng-
lich von Pope [103] vorgeschlagen, da dieser Term die Berechnung einfacher axialsymme-
trischer Jets verbessern soll. Für komplexere Jets kann die Berücksichtigung des Vortex-
Stretching Terms allerdings zu größeren Fehlern führen. Sogar selbst bei einfacheren
Geometrien können numerische Fehler auftreten, wenn die Strömung heiß und kompres-
sibel ist [19]. Für komplizierte Geometrien, wie z. B. dem Blütenmischer, scheint dieses
Modell deshalb nicht oder zumindest nur bedingt geeignet. Sadiki [108] gibt zusätzlich
an, dass dieser eingeführte Term zu einer thermodynamischen Inkonsistenz führen kann.
Die Korrektur von Pope mag zwar für manche Fälle zu Verbesserungen der Simulation
von Abgasstrahlen führen, beeinflusst allerdings das Verhalten der Gleichungen in der
wandnahen Grenzschicht. Da angenommen wird, dass die Standardturbulenzmodelle gu-
te Eigenschaften hinsichtlich der Genauigkeit im wandnahen Bereich haben, führt der
Term von Pope vermutlich zu einer Beeinflussung der Grenzschichteigenschaften. Eine
genaue Berechnung der viskosen Wandkräfte ist jedoch in diesem Falle wichtig für die
Bestimmung der Düsenleistungsparameter.
Ein weiteres Argument, welches für die Verwendung der Temperaturkorrektur von Abdol-
Hamid spricht, ist die Abhängigkeit des Modells von der Totaltemperatur. Statt wie Tam
und Ganesan in ihrem Modell einen Zusammenhang der turbulenten Viskosität zum
Dichte- und Geschwindigkeitsgradienten aufzustellen wird hier gewährleistet (adiabate
Wandrandbedingungen vorausgesetzt), dass in der Grenzschicht wieder die für das k-ǫ
Turbulenzmodell definierte Standardviskosität verwendet wird, da Totaltemperaturgra-
dienten in Wandnähe verschwinden. Außerdem ist durch den Ansatz von Abdol-Hamid
das Korrekturverfahren praktisch unabhängig vom Strömungszustand, da die Totaltem-
peratur als Erhaltungsgröße z. B. auch über Stöße oder Expansionen hinweg für ideale
Gase konstant bleibt.
Temperaturkorrektur für k-ǫ Turbulenzmodell
Als Grundlage dient dem Temperaturkorrekturmodell von Abdol-Hamid das k-ǫ Turbu-
lenzmodell aus Kapitel 3.3.2 mit den Gleichungen (3.33) und (3.34). Jedoch wird die
Konstante Cµ neu definiert in Abhängigkeit eines Faktors CT .
48
Cµ = 0, 09 · CT . (4.1)
Der Faktor CT wiederum ist eine hauptsächlich von der Temperatur abhängige Funktion.
Diese beinhaltet den normierten Totaltemperaturgradienten Tg .

Tg3
CT = 1 + . (4.2)
0, 041 + f (M aT )
Obwohl hier Cµ direkt modifiziert wird, geschieht letztlich bei diesem Korrekturverfahren
nichts anderes als die turbulente Viskosität aus Gleichung (3.36) mit dem Faktor CT zu
multiplizieren. Tg ist eine Funktion des lokalen Totaltemperaturgradienten und wird
durch das lokale turbulente Längenmaß sowie die Totaltemperatur normiert:
|∇Tt | k 3/2
Tg = . (4.3)
Tt ǫ
Die Funktion f (M aT ) die im Nenner von (4.2) auftaucht soll Kompressibilitätseffekte

berücksichtigen und die Ausmischung bei hohen turbulenten Machzahlen reduzieren:

f (M aT ) = M a2T − M a2T 0 H (M aT − M aT 0 ) . (4.4)
H(x) ist die Heaviside-Sprungfunktion. Die Konstante M aT 0 beträgt 0,1. Der Maximal-
wert für CT wurde wie in [109] auf 5 limitiert. Die Limitierung hilft vor allem beim
Anrechnen mit dem Modell, hier können sehr hohe Werte für CT erreicht werden und
zu instabilen Rechnungen führen.
Temperaturkorrektur für k-ω-SST Turbulenzmodell
Die Temperaturkorrektur für das k-ǫ Turbulenzmodell soll auch für das k-ω-SST Modell
von Menter (Kapitel 3.3.2) anwendbar gemacht werden. Dies ist deshalb möglich, weil
das zweischichtige SST Turbulenzmodell im wandfernen Bereich in seiner Formulierung
dem k-ǫ Turbulenzmodell entspricht. Im wandnahen Bereich, wo das SST Modell in die
k-ω Form übergeht ist die Temperaturkorrektur nicht aktiv, der Faktor CT beträgt hier
eins. Mit der Beziehung
ǫ = β ∗ kω. (4.5)
49
kann Gleichung (4.3) in die ω-Form transformiert werden:
|∇Tt | k 1/2
Tg = . (4.6)
Tt β ∗ ω
Der Faktor CT aus Gleichung (4.2) berechnet sich nun mit (4.6). Wie schon bei der
Temperaturkorrektur in Verbindung mit dem k-ǫ Turbulenzmodell wurde auch hier CT
auf einen Maximalwert von 5 beschränkt. Die neue turbulente Viskosität für das k-ω-SST
Modell inklusive Faktor CT lautet dann
ρa1 k
µT = C T . (4.7)
max (a1 ω, SF2 )
Diese neue Viskosität für das SST Modell kann jedoch nicht mit der Transportgleichung
(3.41) für ω verwendet werden. Um konsistent mit dem Effekt der Temperaturkorrektur
für das k-ǫ Modell zu sein, muss diese Gleichung angepasst werden. Der Grund dafür
wird beim Blick auf die exakte Transformation der ǫ-Gleichung (3.33) in die ω Schreib-
weise (siehe z. B. [13] und Anhang B) ersichtlich. Für eine konsistente Übertragung der
Temperaturkorrektur auf das SST Modell müssen der Produktionsterm sowie der Cross-
Diffusionsterm in Gleichung (3.41) ebenfalls mit CT multipliziert werden. Andernfalls ist
eine Unterproduktion von ω bzw. Überproduktion der turbulenten kinetischen Energie
in Gebieten die Folge, wo der Faktor CT ungleich eins ist. Die neue Transportgleichung
für ω mit der Temperaturkorrektur von Abdol-Hamid lautet dann
∂ ∂ e γ ∂Uej
(ρω) + ρω Ui = CT ρτij − βρω 2 + β ∗ ρω 2 M a2T (1 − F1 )
∂t ∂xi νT ∂xi
(4.8)
∂ ∂ω CT ∂k ∂ω
+ (µ + σω µT ) + 2 (1 − F1 ) ρσω2 .
∂xj ∂xj ω ∂xi ∂xi
4.2.2. Ergebnisse
Neben einer qualitativen Bewertung der Ergebnisse wurden die Lösungen auch hinsicht-
lich der Leistungsparameter quantitativ eingeordnet. Dazu dienen die in Kapitel 2.2
definierten Leitungsparameter der Düse. Die Leistungsbeiwerte aus den Simulationen
können mit den Ergebnissen experimenteller Messungen von ASE FluiDyne Aerotest
Laboratory verglichen werden. Des Weiteren wurde noch der Totaldruckverlust am Dü-
senausgang ermittelt. Hier soll der durch die temperaturkorrigierten Turbulenzmodelle
50
hervorgerufene Verlust gegenüber ihrem jeweiligen standard-Pendant quantifiziert wer-

den.
Während die Berechnung des Strömungsfelds mit einer Genauigkeit zweiter Ordnung er-
folgte, wurden die Turbulenzmodelle aus Stabilitäts- und Konvergenzgründen mit erster
Ordnung diskretisiert.
Qualitativer Einfluss der Temperaturkorrekturmethode
Konturplotts der Totaltemperaturverteilung im Längsschnitt durch die Düse für alle ge-
testeten Turbulenzmodelle sind in Abbildung 4.3 dargestellt. Dabei ist die Strömung
für den Fall 1.6-1.1-2.4 gezeigt. Die Variante mit Temperaturkorrektur ist jeweils in der
unteren Bildhälfte dargestellt. Vor allem beim SST Modell ist eine verstärkte Ausmi-
schung im Vergleich zum k-ǫ Ansatz zu erkennen. Die unterschiedliche Form des Mischers
jeweils in der oberen und unteren Bildhälfte ist der Darstellung geschuldet. Um die spie-
gelbildliche Abbildung zu realisieren, wurde in ParaView für die Darstellung der unteren
Bildhälfte das Bild um 180◦ rotiert, ebenso wie die Geometrie des Mischers. Daher ist
dieser in der oberen und unteren Bildhälfte jeweils einmal von der einen, das andere Mal
von der anderen Seite zu sehen, was aber nicht weiter irritieren soll.
Besser erkennen lässt sich die Ausmischung bei Schnitten an verschiedenen x/D Positio-
nen in Bild 4.4 (x/D = 0, 0 befindet sich am Düsenausgang). Das Temperaturprofil beim
SST Modell scheint gegenüber dem k-ǫ Modell prinzipiell diffusiver zu sein. Dies setzt
sich auch bei den Modellen mit Temperaturkorrektur fort. Vergleicht man das k-ǫ und
k-ǫ-TC miteinander, fällt hauptsächlich die Verkürzung des heißen Potentialkerns auf.
Beim k-ω-SST-TC Modell ist die Wirkung der Temperaturkorrektur am stärksten. Be-
reits bei x/D = −0, 7, also auf etwa halbem Weg zwischen Mischer- und Düsenaustritt,
sind Temperaturspitzen vor allem in den Mischerwirbeln bereits deutlich abgebaut. Ein
Grund für die stärkere Ausmischung ist die erhöhte Turbulenz im Nachlauf des Mischers.
Bild 4.5 zeigt wiederum den Längsschnitt durch die Düse. Hier ist die turbulente kine-
tische Energie dargestellt. Es ist zu sehen, dass das SST Modell auch im Nachlauf der
Düsengondel mehr Turbulenz erzeugt. Besser zu erkennen ist der Einfluss der Modelle
auf die Turbulenz im Mischernachlauf in Abbildung 4.6 an unterschiedlichen x/D Sta-
tionen. Während das k-ǫ-TC gegenüber dem standard k-ǫ Modell nur unwesentlich mehr
Turbulenz erzeugt ist der Unterschied zwischen dem k-ω und k-ω-TC deutlich größer.
Besonders auffällig ist, dass die hohe Turbulenz direkt im Nachlauf des Mischers beim
SST-TC Modell ab Verlassen der Düse zu einer geringeren Produktion von k führt als
beim standard SST Modell. Die starke Aufspreizung direkt hinter dem Mischer reduziert
offensichtlich die lokalen Geschwindigkeitsgradienten so, dass nachfolgend die Produkti-
on der turbulenten kinetischen Energie sinkt.
Der höhere Einfluss der Totaltemperatur beim SST Modell liegt hauptsächlich an der
51
Abbildung 4.3.: Totaltemperaturverteilung der OT-Düse für den Fall 1.6-1.1-2.4. Oberes
Bild, oberer Teil: k-ǫ, unterer Teil: k-ǫ-TC. Entsprechend unteres Bild,
obere Hälfte: k-ω-SST, untere Hälfte: k-ω-SST-TC.
Formulierung des normierten Temperaturgradienten in (4.3), bzw. (4.6). Hier geht ne-
ben der Totaltemperatur noch das turbulente Längenmaß k 3/2 /ǫ, bzw. k 1/2 /(β ∗ ω) mit
ein. Letzteres ist im Nachlauf des Mischers größer als beim k-ǫ Turbulenzmodell und
52
Abbildung 4.4.: Totaltemperaturverteilung der OT-Düse für den Fall 1.6-1.1-2.4 an un-
terschiedlichen x/D Positionen. Von links nach rechts: k-ǫ-TC, k-ǫ, k-ω-
SST-TC, k-ω-SST. Von oben nach unten: x/D = −1, 4; −0, 7; 0, 0; 0, 7;
1, 4;
53
Abbildung 4.5.: Verteilung der turbulenten kinetischen Energie der OT-Düse für den Fall
1.6-1.1-2.4. Oberes Bild, oberer Teil: k-ǫ, unterer Teil: k-ǫ-TC. Entspre-
chend unteres Bild, obere Hälfte: k-ω-SST, untere Hälfte: k-ω-SST-TC.
führt somit zu einem größeren Temperaturkorrekturfaktor CT . Die Verteilung von CT

beider temperaturkorrigierter Turbulenzmodelle ist in Abbildung 4.7 zu sehen und liegt
im Mittel bei Werten zwischen ungefähr zwei und drei.
54
Abbildung 4.6.: Verteilung der turbulenten kinetischen Energie der OT-Düse für den Fall
1.6-1.1-2.4 an unterschiedlichen x/D Positionen. Von links nach rechts:
k-ǫ-TC, k-ǫ, k-ω-SST-TC, k-ω-SST. Von oben nach unten: x/D = −1, 4;
−0, 7; 0, 0; 0, 7; 1, 4;
55
Abbildung 4.7.: Verteilung des temperaturabhängigen Faktors CT der OT-Düse für den
Fall 1.6-1.1-2.4. Oberer Teil: k-ǫ-TC, untere Hälfte: k-ω-SST-TC.
Quantitativer Einfluss der Temperaturkorrekturmethode
Im folgenden Abschnitt soll die Düsenströmung mittels berechneter Leistungsparameter

bewertet und mit Experimenten verglichen werden. Die experimentell gemessenen Leis-
tungsbeiwerte wurden von Rolls-Royce Deutschland zur Verfügung gestellt. Die Mes-
sungen erfolgten im Channel 11 bei ASE FluiDyne Aerotest Laboratory in Plymouth,
Minnesota in den USA. Für mehr Details zu den Experimenten siehe [110].
In den Abbildungen 4.8 und 4.9 sind die Durchflusskoeffizienten über dem Totaltem-
peraturverhältnis geplottet. Für die temperaturkorrigierten Modelle wurden hier keine
Rechnungen für TR=1,0 durchgeführt, da die Ergebnisse in vorangegangenen Testrech-
nungen immer identisch zu den Standardmodellen waren. Es ist zu erkennen, dass der
Durchflusskoeffizient für höhere Temperaturen des Kernstroms abfällt. Dies liegt an der
Zunahme von Verlusten, die mit der Erhöhung von TR zunehmen. Für eine verlustfreie
Düsenströmung würde der Durchflusskoeffizient eins betragen. Die ansteigenden Verlus-
te bei Erhöhung der thermischen Energie des Kernstroms können z. B. damit begründet
werden, dass der steigende Wärmestrom vom Heiß- zum Kaltstrahl zu einer höheren
Entropieproduktion innerhalb des Kontrollvolumens der Düse führt [111, 112]. Außer-
dem führt die Erhöhung der Temperatur auch gleichzeitig zu einer Beschleunigung des
56
Abbildung 4.8.: CD über TR für CNPR=1,6.
Kernstroms. Somit werden sich dadurch auch die viskosen Verluste erhöhen, die wieder-
um den Massenstrom reduzieren. Die Modelle mit Temperaturkorrektur führen durch
ihre verbesserte Ausmischung für höhere TR zu geringeren Durchflusskoeffizienten als
das jeweilige Standardmodell. Dass das k-ω-SST-TC am stärkesten mischt kann eben-
falls in den Diagrammen abgelesen werden. Es führt zu den größten Verlusten bzw. zum
geringsten CD .
Verglichen mit den Ergebnissen aus den Experimenten liefert die CFD durchweg zu ge-
ringe Werte für den Durchflusskoeffizienten. Allerdings scheint der qualitative Verlauf
der Kurven mit Temperaturkorrektur die Experimente besser zu treffen.
Der Verlauf des Geschwindigkeitskoeffizienten über TR ist in den Abbildungen 4.10 und
4.11 dargestellt. CV als Verhältnis von erreichtem Bruttoschub des gemischten Systems
zum ideal ungemischten Schub kann hier auch als eine Art Mischungseffizienz gesehen
werden. Je effizienter der Mischer funktioniert und der Grad der Ausmischung ist, desto
höher fallen die Werte für CV aus. Die berechneten Werte für CV die mit dem SST Mo-
dell erreicht werden liegen leicht oberhalb der Kurven des k-ǫ Turbulenzmodells. Somit
spiegelt der Verlauf auch hier die prinzipiell höhere Ausmischung in den Konturplots
des SST gegenüber dem k-ǫ Modell wider. Ebenso lässt sich in den Kurven für den Ge-
57
Abbildung 4.9.: CD über TR für CNPR=2,6.
schwindigkeitskoeffizienten die Wirkung der Temperaturkorrektur gut erkennen.

Für die unterkritische Düsenströmung ist wie zuvor bei CD ein Versatz der Kurven aus
der CFD zum Experiment zu sehen. Die Ergebnisse aus der Simulation liegen unterhalb
der FluiDyne Messungen. Für den Fall CNPR=2,6 ist die Übereinstimmung besser. Le-
diglich die Punkte für TR=1,0 weichen stärker ab.
Lenkt man den Blick noch einmal auf die Kaltrechnung (TR=1,0), ist es interessant zu
sehen, dass sowohl das k-ǫ als auch das k-ω-SST Modell nahezu identische Werte lie-
fern. CV und CD sind sowohl für die über- als auch unterkritische Düsenströmung fast
gleich. Die im Strömungsfeld höhere Turbulenz bei den SST Modellen scheint keinen
großen Einfluss bei den Kaltrechnungen zu haben, wohl aber bei den Heißrechnungen.
Es ist deshalb naheliegend zu folgern, dass die Unterschiede bei den Leistungsparametern
hin zu heißeren Rechnungen (TR>1,0) hauptsächlich aus der Energiegleichung rühren.
Nämlich durch den Zusammenhang zwischen Turbulenz (in Form von turbulenter Vis-
kosität) und der Wärmeleitung. So führt wohl zum größten Teil die bessere thermische
Ausmischung der einzelnen Turbulenzmodelle zu der Diskrepanz von CV und CD bei
den Heißrechnungen zwischen SST und SKE.
Die Abweichung der Ergebnisse zwischen den Messungen und CFD kann mehrere Grün-
58
Abbildung 4.10.: CV über TR für CNPR=1,6.
de haben. Zum einen weicht die im Experiment verwendete Geometrie von der aus der
Simulation ab. Aufgrund eines Fertigungsfehlers seitens FluiDyne hat der Mischer im
CFD Modell eine andere Taktung der Mischerblüten als im Experiment. Das Netz für die
von Rolls-Royce neu verteilten Geometrie wurde nicht mehr rechtzeitig fertiggestellt. Al-
lerdings kann diese geometrische Abweichung wohl vernachlässigt werden und hat keinen
Einfluss auf die Leistungsparameter, wie in Kapitel 5 noch gezeigt wird. Eine weitere
Erklärung für das Abweichen der Messungen könnte an der Position der Datenerfas-
sung liegen. Während die Experimente an den Charging-Planes, kurz vor Eintritt in
den Mischer, ausgewertet wurden, wurden diese Ebenen bei der Netzgenerierung nicht
berücksichtigt und die Daten in der CFD an den numerischen Einlässen erfasst. Eine
nachträgliche Auswertung der CFD Ergebnisse an den Charging-Planes wäre nur mit
sehr großem Aufwand möglich gewesen. Ferner wurde im Experiment ein sogenannter
Spacer installiert. Dieser soll die thermische Ausdehnung des experimentellen Aufbaus
berücksichtigen und wurde für ein Totaltemperaturverhältnis von 2,6 optimiert. Somit
sollten die experimentellen Daten mit TR nahe diesem Wert am genauesten sein. Weitere
Gründe für Abweichungen können auch Messungenauigkeiten sein. So wurden die Tem-
59
Abbildung 4.11.: CV über TR für CNPR=2,6.
peraturen für die kalten Fälle über eine Korrekturformel bestimmt [110]. Wie in Kapitel
5 gezeigt wird, führt auch die Annahme einer konstanten spezifischen Wärmekapazität
in der Simulation bei höheren Totaltemperaturverhältnissen zu Abweichungen. Zuletzt
gibt es in der CFD prinzipiell immer Unsicherheiten auf Grund numerischer Einflüsse
wie z. B. Gitterabhängigkeit, Diskretisierungsschemata, etc.. Die Einflüsse dieser Para-
meter werden zum Teil in der Falluntersuchung in Kapitel 5 noch analysiert.
Abschließend soll noch der Totaldruckverlust der Korrekturmodelle beziffert werden.
Gemeint ist die Differenz des massengemittelten Totaldrucks an der Düsenauslassfläche
des jeweiligen temperaturkorrigierten Turbulenzmodells gegenüber der nicht korrigierten
Berechnung. Das Ergebnis ist in Schaubild 4.12 dargestellt. Wie zu sehen ist, führen die
temperaturkorrigierten Modelle zu zunehmendem Totaldruckverlust für steigende TR.
Der Anstieg beim k-ω-SST scheint degressiv zu sein, während er beim k-ǫ Modell pro-
gressiv ist. Ob die Düsenströmung über- oder unterkritisch ist, scheint dabei allerdings
eine vernachlässigbare Rolle zu spielen. Die Verluste können sich zum einen durch die hö-
here Turbulenz und turbulente Viskosität innerhalb des Düsenkontrollvolumens erklären
lassen. Zum anderen steigt bei den Modellen mit Temperaturkorrektur der Wärmetrans-
fer in der Mischungsscherschicht, was zu einer höheren Entropieproduktion führt. Es sei
60
Abbildung 4.12.: Prozentualer Totaldruckverlust am Düsenaustritt der temperaturkor-

rigierten Turbulenzmodelle gegenüber den Standardmodellen.
hier darauf hingewiesen, dass das hier präsentierte Ergebnis hinsichtlich des Totaldruck-
verlusts auf Grund eines Auswertungsfehlers von dem bereits veröffentlichten Resultat
in [113] abweicht.
4.2.3. Zusammenfassung - Temperaturkorrektur
Die Berücksichtigung thermischer Turbulenz führt zu einer stärkeren Ausmischung ani-

sothermer Strömungen. Mit dem Ansatz von Abdol-Hamid konnte gezeigt werden, wel-
chen Einfluss ein solches Modell auf eine Düsenströmung hat. Zum einen werden unter
Berücksichtigung des Temperaturturbulenzmodells der heiße Kern- und der kalte Fan-
strom qualitativ besser gemischt. Zum anderen hat die effektivere Vermischung aber
auch einen quantitativen Einfluss auf die Leistungsparameter. Die größere Mischungs-
rate führt bei Heißrechnungen zu geringeren Durchflusskoeffizienten und zu höheren
Geschwindigkeitskoeffizienten. Es wurde gezeigt, dass das k-ω-SST Turbulenzmodell in
dieser Konfiguration hier prinzipiell stärker mischt als das k-ǫ Modell. Die Berücksich-
tigung von Temperaturturbulenz führt nochmals zu einer Verstärkung dieser Tatsache.
In Kapitel 5 wird allerdings noch gezeigt, dass dieses Verhalten in Verbindung mit der
räumlichen Diskretisierung in Zusammenhang gebracht werden muss.
61
Bei den Kaltrechnungen wurde beobachtet, dass hier beide Turbulenzmodelle fast diesel-
ben Ergebnisse für CV und CD liefern. Dies legt den Schluss nahe, dass die Unterschiede
zwischen SST und k-ǫ Modell hauptsächlich aus der Energiegleichung kommen. Da das
SST Modell im Nachlauf des Mischers mehr Turbulenz erzeugt, verbessert es auch die
thermische Ausmischung.
Die Werte aus den Messungen konnten nur für CV und die überkritische Düsenströmung
gut getroffen werden. Für die restlichen Rechnungen liefert die CFD zu geringe Werte.
Allerdings scheinen die Modelle mit Temperaturkorrektur nun qualitativ besser mit den
FluiDyne Messungen übereinzustimmen. Um näher an den experimentellen Daten zu
liegen müssen außerdem numerische Fehler und Ungenauigkeiten minimiert werden. In
Kapitel 5 wird der Einfluss solcher Parameter zum Teil quantifiziert. Es soll jedoch noch
betont werden, dass auch die experimentellen Daten nicht über alle Zweifel erhaben sind.
So müssen vor allem die Messwerte für die kalten Düsenströmungen kritisch betrachtet
werden.
Bei der Auswertung des massengemittelten Totaldrucks am Düsenauslass wurde gezeigt,
dass die Modelle mit Temperaturkorrektur gegenüber den Standardmodellen zu einem
Totaldruckverlust führen. Für das SKE nimmt dieser mit steigendem TR progressiv zu,
für das SST hingegen degressiv.
Die hier gezeigten Ergebnisse zur Temperaturkorrekturmethode wurden in zwei Konfe-
renzen vorgestellt und veröffentlicht [110, 113].
62
4.3. Simulation der Düsenströmung mit verbessertem

Wärmestromterm
4.3.1. Verbesserter turbulenter Wärmestromschließungsansatz
In Kapitel 4.2 wurde versucht, die Ausmischung anisothermer Scherschichten zu verbes-

sern, indem ein Modell zur Berücksichtigung von Temperaturturbulenz eingesetzt wur-
de. Es konnte gezeigt werden, dass durch das Modell von Abdol-Hamid [18] der Heiß-
und Kaltstrahl besser ausmischt und sich die effektivere Vermischung auch auf die Leis-
tungsparameter der Düse auswirkt. Der Ansatz des Modells von Abdol-Hamid, ebenso
wie der von Tam und Ganesan [22], den Einfluss einer stark anistothermen Scherschicht
auf die Strömung durch eine gezielte Erhöhung der turbulenten Viskosität zu erfassen,
erscheint sinnvoll. Wird nämlich angenommen, dass die in Experimenten [20, 21] be-
obachtete verstärkte Ausmischung gegenüber isothermer Scherschichten auf zusätzliche
Turbulenzeffekte zurückzuführen ist, sollte ein entsprechendes Modell auch Einfluss auf
die Impulsgleichungen haben. Über eine lokal modifizierte turbulente Viskosität kann
dies gewährleistet werden.
Allerdings benutzen beide Modelle zur Schließung des turbulenten Wärmestromvektors
das in Kapitel 3.3.1 bereits beschriebene Eddy Diffusivity Modell (EDM). Diese ver-
gleichsweise sehr einfache Schließungsmethode zur Modellierung des turbulenten Wärme-
stromterms in der Energieerhaltungsgleichung bringt einige Nachteile mit sich. So ist für
diesen Ansatz eine turbulente Prandtlzahl zu bestimmen. Sowohl Abdol-Hamid als auch
Tam und Ganesan nehmen diese als konstant an mit P rT = 0, 9, bzw. 0, 422. Allerdings
ist die Annahme einer konstanten turbulenten Prandtlzahl nicht für das gesamte Strö-
mungsgebiet zutreffend. In [13] wird z. B. ein Wert von P rT = 0, 9 in der Grenzschicht
als geeignet angegeben. In Freistrahlen hingegen sei ein Wert von P rT = 0, 5 angemesse-
ner. Birch et al. [19] dagegen differenzieren den Wert für die turbulente Prandtlzahl im
Freistrahl noch einmal genauer. In axialsymmetrischen Strahlen liegt der Wert bei ca.
0,7, in voll turbulenten Freistrahlen beträgt P rT gemäß Birch nur 0,4. Dembowski et al.
[114] zeigen in einer Studie unter anderem welchen Einfluss unterschiedliche turbulente
Prandtlzahlen auf die Länge des Potentialkerns haben.
Die Verwendung einer konstanten turbulenten Prandtlzahl in einem der Temperaturkor-
rekturmodelle führt also zwangsläufig dazu, dass die turbulente Viskosität nicht überall
korrekt berechnet wird. Im Falle von Abdol-Hamid kann dies in Freistrahlen zu einer
Überschätzung des tatsächlichen µT führen. Bei Tam und Ganesan wird überall dort,
wo in der Realität größere Werte als P rT = 0, 422 vorkommen, die turbulente Viskosität
eventuell unterschätzt.
Im Folgenden sollen nun zwei unterschiedliche Philosophien zur Herangehensweise der
63
Modellierung eines verbesserten turbulenten Wärmestroms vorgestellt und kritisch be-

urteilt werden.
Modellierung einer turbulenten Wärmeleitzahl/Prandtlzahl
Für eine bessere Bestimmung der turbulenten Prandtlzahl bzw. Verbesserung der Wär-
meleitung basierend auf dem EDM existieren mehrere Modelle. Die Ansätze reichen
von algebraischen Modellen zur Vorhersage bis hin zu differentiellen Methoden. Kays
[23] findet eine algebraische Formel für eine variable turbulente Prandtlzahl. Weitere
algebraische Modelle finden sich in [24, 25, 26, 27]. Interessant ist auch der Ansatz
zur differentiellen Bestimmung der turbulenten Prandtlzahl oder der Wärmeleitzahl. So
z.B. die Autoren Nagano und Kim [115], und darauf aufbauend Youseff et al. [116] sowie
Deng et al. [117]. Diese formulieren die turbulente Wärmeleitfähigkeit als Funktion des
Zeitmaßstabs der Geschwindigkeits- und der Temperaturturbulenz in folgender Form:
m !n
k θ2
λT ∝ k , (4.9)
ǫ ǫθ
wobei m und n hier konstante Exponenten sind. Die Größen θ2 und ǫθ sind die gemittel-
te Temperaturvarianz sowie deren Dissipationsrate. Die beiden Größen werden mittels
entsprechender Transportgleichungen bestimmt.
Chidambaram et al. [28] verwendeten diesen Ansatz wahrscheinlich als Erste zur Vor-
hersage der turbulenten Prandtlzahl von Heißstrahlströmungen. Anstatt einer üblichen
konstanten turbulenten Prandtlzahl wird hier ein variables P rT formuliert. Mit λT =
cp µT /P rT und in Anlehnung an Gleichung (4.9) wurde die turbulente Prandtlzahl wie
folgt definiert:
r
Cµ k 2ǫθ
P rT = . (4.10)
Cλ ǫ θ2
Cµ ist die Konstante aus dem k-ǫ Modell, Cλ ist eine weitere Modellkonstante. Auch hier
sind zusätzlich zu einem Turbulenzmodell wie dem k-ǫ Modell die Lösung der Transport-
gleichungen für die Temperaturvarianz und deren Dissipationsrate notwendig. Während
Chidambaram et al. mehr die generelle Anwendbarkeit des Modells an einfachen Geome-
trien testen, versuchen sich Kenzakowski et al. [29] an anwendungsnäheren Testrechnun-
gen. Jedoch zeigt das Modell nur teilweise Erfolg. Für eine singuläre Düse mit M a = 2, 0
scheint das Modell zu funktionieren. Bei Simulationen von Koaxialdüsen, die mit einer
64
Austrittsgeschwindigkeit von M a = 0, 85 getestet wurden konnten keine sinnvollen Er-

gebnisse erzielt werden. In der Scherschicht zwischen Bypass- und Kernstrom scheint das
Modell zu versagen. Weitere Tests des Modells mit der Anwendung auf Düsenströmun-
gen brachten trotz einigen Modellverfeinerungen allerdings weiterhin nur wenig Erfolg
[118, 119, 120].
Modellierung einer turbulenten Wärmestromgleichung
Selbst wenn es gelingen sollte ein Modell zu finden, das die turbulente Prandtlzahl in al-
len Bereichen des Rechengebiets korrekt vorhersagt, bringt das Eddy Diffusivity Modell
immer noch ein unumgängliches Defizit mit sich. Dieses liegt an der möglichen Wär-
meleitrichtung die sich gemäß des EDM aus Gleichung (3.24) ergeben kann. So ist die
Richtung der Wärmeübertragung praktisch vorgegeben durch den Temperaturgradien-
ten. Messungen [30, 31] ergeben jedoch, dass selbst wenn der Temperaturgradient in
Strömungsrichtung viel kleiner ist als normal dazu, der turbulente Wärmestrom parallel
zur Fließrichtung um ein Vielfaches größer sein kann.
Ein potentieller Lösungsweg könnte die direkte Formulierung des turbulenten Wärme-
stromterms über einen differentiellen Ansatz darstellen. Von Launder [66, 121] wird
dieses Verfahren als Diffenrential Second-Moment Closure (DSM) bezeichnet. Die Diffe-
rentialgleichung für den turbulenten Wärmestrom wird über eine Momentenmethode ge-
wonnen. Dabei ist die Vorgehensweise analog zur Herleitung des Reynolds-Stressmodells.
Die Gleichung für das DSM kann gewonnen werden, indem die Impulsgleichung mit der
fluktuierenden Größe der Temperatur θ und die thermische Energiegleichung mit der
fluktuierenden Geschwindigkeitsgröße multipliziert wird. Beide Gleichungen addiert und
gemittelt ergibt die gewünschte Transportgleichung für den turbulenten Transport der
Größe θ. Einen gute Übersicht über diesen Ansatz und unterschiedliche Schließungsme-
thoden bietet z. B. [122].
Dieser Ansatz bietet grundsätzlich großes Potential den turbulenten Wärmestromvek-
tor physikalischer zu modellieren als es durch das Eddy Diffusivity Modell möglich ist.
Zum einen weil keine turbulente Prandtlzahl bestimmt werden muss, zum anderen weil
die Multidimensionalität der Wärmeleitrichtung bei dieser Methode besser berücksich-
tigt wird. Aufgrund der Vorteile und des grundsätzlichen Potentials des DSM soll dieser
Ansatz in dieser Arbeit verwendet werden.
Als Basis für das verwendete Modell dient die Veröffentlichung von Lai und So [123].
Diese geben für den inkompressiblen Transport von u′i θ folgende Gleichung an:
65
" #
∂ e ′ ∂ ∂ θ∂u ′
j u ′
j ∂θ
Uk uj θ = − u′ u′ θ + ν +λ
∂xk ∂xk j k ∂xk ∂xk ∂xk
(4.11)
∂ Te ej
∂U θ ∂p′ ∂u′j ∂θ
−u′j u′k −u′k θ − − (ν + α) ,
∂xk ∂xk ρ ∂xj ∂xk ∂xk
oder in symbolischer Schreibweise
t v
Cjθ = Djθ + Djθ + Pjθ,1 + Pjθ,2 + Φ∗jθ − ǫjθ . (4.12)
Die Bedeutung der Terme in Gleichung (4.12) ist von links nach rechts der konvektive
Transport, der turbulent-diffusive Transport, der viskos-diffusive Transport, Produktion
durch Temperaturgradienten, Produktion durch Scherung, Pressure-Scrambling und die
molekulare Dissipation. Da in den in dieser Arbeit verwendeten Rechennetzen die Grenz-
schicht nicht voll aufgelöst, sondern mit Wandfunktionen überbrückt wird, kann für das
DSM die sogenannte High-Reynolds Version verwendet werden. Gemäß [124] können für
High-Reynolds Strömungen und isotroper Turbulenz der zweite und sechste Term auf der
rechten Seite von Gleichung (4.12) vernachlässigt werden. Die einzigen Terme die es noch
zu modellieren gilt sind somit die turbulente Diffusion und das Pressure-Scrambling.
Die modellierten Terme und Konstanten wurden alle von [123] übernommen. Der tur-
bulent diffusive Transport wurde mit dem Ausdruck
!
k ′ ′ ∂uk θ
′ ∂u′j θ
− u′j u′k θ = cθs uj ul + uk ul
′ ′
(4.13)
ǫ ∂xl ∂xl
modelliert. Dieser Ansatz wird auch in [125] vorgeschlagen. Die Modellkonstante cθs wird
mit 0,11 angegeben. Laut [126] gibt dieses Modell prinzipiell eine gute Übereinstimmung
mit Messwerten.
Für die Modellierung des Pressure-Scrambling Terms wurde nur der High-Reynolds An-
teil gemäß [123] berücksichtigt. So ergibt sich
Φ∗jθ = Φjθ , (4.14)
mit der Unterteilung
Φjθ = Φjθ,1 + Φjθ,2 + Φjθ,w . (4.15)
66
Φjθ,1 repräsentiert den fluktuierenden Anteil, Φjθ,2 die Spannungskomponente und Φjθ,w
ist der Wandreflektionsterm. Φjθ,1 ist modelliert durch
ǫ
Φjθ,1 = −c1θ u′j θ (4.16)
k
und Φjθ,2 durch
ej
∂U
Φjθ,2 = −c2θ Pjθ,2 = c2θ u′k θ , (4.17)
∂xk
mit den Modellkonstanten c1θ = 3, 0 und c2θ = 0, 4. Lai und So erzielten sehr gute Ergeb-
nisse unter Vernachlässigung von Φjθ,w und stellen diesen Term stark in Frage. Sie folgern
aus ihren Ergebnissen, dass der Wandreflektionsterm unnötig ist, zumindest unter Ver-
wendung von Gleichung (4.16) und (4.17) für das Modell des Pressure-Scramblings. Auf
Grund dieser Folgerung und der Unsicherheit bei der Berücksichtigung von Φjθ,w , wurde
dieser Term hier ebenfalls vernachlässigt.
Die Verwendung eines differentiellen Schließungsansatzes für den turbulenten Wärme-
strom, zusammen mit einem Zweigleichungsturbulenzmodell und in Verbindung mit einer
stark kompressible Strömung, ist vermutlich der erste Kombinationsversuch dieser Art.
Folglich gibt es auch in der Literatur keine speziell auf diese Anwendung abgestimmten
Modelle und Koeffizienten, denen man sich bedienen könnte. So variieren z. B. die in
anderen Arbeiten verwendete Koeffizienten für c1θ mit Werten zwischen 2,5 und 9,7 recht
stark [124, 127, 128]. Es wurde daher gerade für diese erste Anstrengung versucht, den
Ansatz für die Formulierung des turbulenten Wärmestroms im Sinne von Allgemeingül-
tigkeit und im Hinblick auf etwaige Modellierungsunsicherheiten eher konservativer zu
gestalten. Auch das Handbook of Fluid Dynamics [124] gibt diese hier gewählte For-
mulierung für das Pressure-Scrambling mit den gewählten Koeffizienten c1θ = 3, 0 und
c2θ = 0, 4 als den wohl gebräuchlichsten Ansatz an.
Mit allen Schließungsansätzen und getroffenen Annahmen lautet die modellierte Trans-
portgleichung inklusive zeitlicher Ableitung jetzt
" !#
∂ ∂ ∂ k ∂u′ θ ∂u′j θ
u′j θ + Uek u′ θ = cθs u′j u′l k + u′k u′l
j
∂t ∂xk ∂xk ǫ ∂xl ∂xl
(4.18)
∂ Te ej
∂U ǫ ej
∂U
− u′j u′k − u′k θ − c1θ u′j θ + c2θ u′k θ .
∂xk ∂xk k ∂xk
67
Zum Berechnen der Düsenströmung muss Gleichung (4.18) allerdings auch für kom-
pressible Strömungen anwendbar sein. D.h. die Gleichung muss in Erhaltungsform für
kompressible Medien umgeschrieben werden. Gleichung (4.18) zusammen mit der An-
nahme ρτij = −ρu′′i u′′j aus (3.18) ergibt für den turbulenten Transport von θ für den
kompressiblen Fall
" !#
∂ ∂ ∂u′′j θ
ρu′′j θ + Uek ρu′′ θ = ∂ cθs k −ρτjl
∂u′′k θ
− ρτkl
j
∂t ∂xk ∂xk ǫ ∂xl ∂xl
(4.19)
∂ Te ej
∂U ǫ ej
∂U
+ ρτjk − ρu′′k θ − c1θ ρu′′j θ + c2θ ρu′′k θ .
∂xk ∂xk k ∂xk
Strenggenommen bedarf es im Gegensatz zur inkompressiblen, Reynolds gemittelten

Gleichung (4.11) für den kompressiblen Fall die Favre gemittelte Transportgleichung
inklusive entsprechender Schließung. Die Herleitung der Favre gemittelten Transport-
gleichung ist im Anhang A dargestellt. Hier wird gezeigt, dass die Gleichung wie sie in
(4.19) steht für kompressible Strömungen Gültigkeit besitzt.
Für den turbulenten Wärmestromvektor wird letztlich folgende Beziehung festgelegt:
qT j = cp ρu′′j θ. (4.20)
4.3.2. Ergebnisse
Wie zuvor beim Temperaturkorrekturverfahren werden die Ergebnisse wieder qualitativ

und mittels berechneter Leistungsparameter aus Kapitel 2.2 quantitativ bewertet. Zum
Vergleich sind bei den Leistungskurven wieder die bereits aus Abschnitt 4.2 bekannten,
experimentell ermittelten Leistungsparameter geplottet.
In diesem Abschnitt wird für die Berücksichtigung kinetischer Turbulenz das k-ǫ Modell
verwendet. Die Berechnung des Strömungs- und des Turbulenzfelds erfolgt mit einer Ge-
nauigkeit von zweiter Ordnung im Raum. Ebenso die Transportgleichung in (4.19), bis
auf den Term für die turbulente Diffusion. Aus Stabilitätsgründen wird hier die Diver-
genz des Tensorfelds mit erster Ordnung diskretisiert. Da adiabate Wände angenommen
werden, wird der Wandwärmestrom zu null gesetzt: u′′j θ|w = (0, 0, 0).
68
Qualitative Bewertung des verbesserten Wärmestrommodells
Die Konturplots der statischen Temperatur sind in Abbildung 4.13 im Längsschnitt

durch die Düsengeometrie dargestellt. Es ist zu erkennen, wie die Temperaturverteilung
für P rT = 0, 5 im mittleren und das differentielle Modell im unteren Bild insgesamt
diffusiver wird als für die Referenzrechnung mit P rT = 0, 9. Bei einem direkten Ver-
Abbildung 4.13.: Verteilung der statischen Temperatur der OT-Düse für den Fall 1.6-1.1-
2.4 im Längsschnitt. Oberes Bild P rT = 0, 9; mittleres Bild P rT = 0, 5;
unteres Bild DSM.
gleich des EDM mit P rT = 0, 5 und dem DSM Modell fällt bei genauerer Betrachtung
auf, dass beide Modelle in verschiedenen Bereichen unterschiedlich stark mischen. So
mischt das DSM Modell z. B. im direkten Nachlauf des Mischers stärker als das EDM.
69
In diesem Bereich wird angenommen, dass die Strömung voll turbulent und dreidimen-
sional ist. Gemäß Birch et al. [19] macht dieses Verhalten des differentiellen Modells hier
Sinn. Birch et al. geben an, dass die thermische Mischung in voll turbulenten Scher-
schichten sehr groß werden kann. Als Zahlenwert wird hier eine turbulente Prandtlzahl
von ca. 0,4 angegeben und liegt somit unterhalb des Werts für das Eddy Diffusivity
Modells. Entlang des heißen Potentialkerns in der Nähe des Plugs liegt die Mischung
des DSM hingegen zwischen den Rechnungen mit P rT = 0, 9 und P rT = 0, 5. Auch
dies erscheint sinnvoll, denn in Achsennähe nahe des Plugs nähert sich der Strahl einer
axialsymmetrischen Form. Die Stromlinien verlaufen mehr geschichtet und haben einen
eher zweidimensionalen Charakter. Der in [19] für axialsymmetrische Strömungsbereiche
angegebene Wert für die turbulente Prandtlzahl liegt bei ca. 0,7 und damit zwischen den
verwendeten Werten für das EDM. Besser zu erkennen ist der Vergleich in Abbildung
4.14, wo die Temperaturkonturen an unterschiedlichen Ebenen entlang der x-Achse ge-
zeigt werden. Im Nachlauf des Plugs mischt das EDM mit P rT = 0, 5 bis in etwa zum
Düsenausgang leicht stärker als das DSM. Danach nimmt die Mischung des differenti-
ellen Modells gegenüber dem algebraischen Modell wieder leicht zu. Gut erkennbar ist
auch der Unterschied zwischen EDM und DSM im Nachlauf des Mischers. Besonders gut
zu erkennen ist dies auf ca. halbem Weg zwischen Mischeraustritt und Düsenausgang bei
x/D=-0,7, wo die Temperaturspitzen beim DSM in den Mischerwirbeln am geringsten
sind.
Eine Verteilung der berechneten Größe u′′j θ ist in Abbildung 4.15 zu sehen. In diesen
Konturplotts wird die Dreidimensionalität des turbulenten Wärmestromvektors deut-
lich (dabei zeigt x in Strömungsrichtung, y in die Bildebene hinein und z nach oben). Es
ist ersichtlich, dass die Werte für u′′1 θ parallel zur Hauptströmungsrichtung große Werte
annehmen können, was einen großen thermischen Transport qT1 in diese Richtung zur
Folge hat. Der Betrag von u′′1 θ scheint in den meisten Gebieten des Strömungsfelds größer
zu sein als u′′2 θ und u′′3 θ, die hauptsächlich normal zur Hauptströmungsrichtung zeigen.
Dadurch, dass die Richtung der Wärmeübertragung beim EDM nur vom Temperatur-
gradienten abhängig ist, wäre die Wärmestromkomponente qT1 hier erwartungsgemäß
quasi null.
Quantitative Bewertung des verbesserten Wärmestrommodells
In Schaubild 4.16 und 4.17 sind die Durchflusskoeffizienten dargestellt. Die Charakte-
ristik der Kurven entspricht der aus Abschnitt 4.2 und zeigt sinkende Durchflusskoeffi-
zienten für steigende Totaltemperaturverhältnisse. Der Vollständigkeit halber sind auch
wieder die von FluiDyne experimentell ermittelten Leistungsbeiwerte aufgetragen. Auf
deren Verläufe und Abweichung zu den Simulationen soll in diesem Abschnitt nicht mehr
eingegangen werden. Es wird hier auf die Diskussion in Kapitel 4.2 verwiesen. Die beiden
70
Abbildung 4.14.: Temperaturprofile der OT-Düse für den Fall 1.6-1.1-2.4. Von links nach
rechts: P rT =0,9, P rT =0,5, DSM. Von oben nach unten: x/D=-1,4, -
0,7, 0,0, 0,7, 1,4 (x/D=0,0 ist Düsenaustrittsfläche).
71
Abbildung 4.15.: Verteilung von u′′j θ der OT-Düse für den Fall 1.6-1.1-2.4 im Längs-
schnitt. Oberes Bild u′′1 θ; mittleres Bild u′′2 θ; unteres Bild u′′3 θ.
Kurven des EDM liefern dieselben Ergebnisse für die Kaltrechnungen (TR=0,0). In Rich-
tung steigender TR divergieren beide Kurven aber auseinander und so liefert das Modell
mit besserer Ausmischung für P rT = 0, 5 gegenüber 0,9 geringere Werte für CD . Für
das DSM fällt allerdings auf, dass bereits bei den Kaltrechnungen die Punkte unterhalb
des algebraischen Modells liegen. Erhöht man die Temperatur am Kernstrom nähern
sich die Ergebnisse denen für P rT = 0, 5 an. Bei der Darstellung des Geschwindigkeits-
koeffizienten in Abbildung 4.18 und 4.19 zeigt sich ähnliches. Die Kurvencharakteristik
ist aus Kapitel 4.2 wieder bekannt. Bei derselben Totaltemperatur am Fan- und Kern-
stromeintritt ist CV quasi unabhängig von der turbulenten Prandtlzahl. Erhöht man
72
Abbildung 4.16.: Durchflusskoeffizient über TR.
hingegen TR, spiegelt sich die bessere Ausmischung in stärker ansteigenden Geschwin-
digkeitskoeffizienten wider. Auch bei den Plots von CV sieht man, dass das DSM bereits
für die Kaltrechnungen höhere Werte liefert. Folgt man der Kurve zu höheren Totaltem-
peraturverhältnissen, liegt der Geschwindigeitskoeffizent am Ende leicht oberhalb des
EDM mit P rT = 0, 5. Bei der Auswertung der Rechnungen fiel auf, dass das einfachere
algebraische Schließungsmodell für den turbulenten Wärmestrom leicht besser konver-
giert als das differentielle Verfahren. Während die Massenstrombilanz über alle Ein-
und Auslässe bei allen Rechnungen auf einem ähnlich niedrigen Niveau liegt, fällt bei
der Totalenthalpieflussbilanz ein kleiner Unterscheid auf. Tabelle 4.2 zeigt die positive
Totalenthalpieflussdifferenz zwischen den Eintrittsflächen von Bypass- und Kernstrom
gegenüber dem Düsenaustritt. Für das Eddy Diffusivity Modell liegt die Differenz für
alle Betriebspunkte auf einem ähnlichen Niveau. Die Totalenthalpie steigt innerhalb der
Düse zwischen Eintritts- und Austrittsflächen etwa gleich stark an. Beim DSM fällt auf,
dass die Differenz bei den Kaltrechnungen größer ist als bei den Heißrechnungen. Für
die hohen Totaltemperaturverhältnisse gleichen sich die Werte für die Totalenthalpie-
flussdifferenz zwischen Ein- und Austritt dem algebraischen Modell in etwa an. Diese
zusätzliche Energie, welche beim differentiellen Modell innerhalb des Düsenkontrollvo-
73
Abbildung 4.17.: Durchflusskoeffizient über TR.
lumens entsteht und vermutlich aus zu hoher numerischer Dissipation resultiert, könnte
die Diskrepanz der Ergebnisse für CD und CV bei den Rechnungen mit geringerem TR
erklären, da die Leistungsparameter erfahrungsgemäß äußerst sensibel sind. Die Unter-
schiede zwischen Heiß- und Kaltrechnungen beim DSM könnten durch leicht stärkere
Überschwinger der Temperaturverteilung in den thermischen Grenzschichten erklärbar
sein. Es zeigt sich, dass sich die Wandgrenzschicht beim differentiellen Ansatz leicht
stärker erhitzt als beim algebraischen Modell. Die Wandtemperaturen für den Fan- und
Kernstrombereich sind exemplarisch für den Fall 1.6-1.1-2.4 in Abbildung 4.20 darge-
stellt. Da bei den Heißrechnungen die Totalenthalpie aufgrund des heißeren Kernstroms
insgesamt auf einem deutlich höheren Niveau ist, fallen möglicherweise die Einflüsse aus
der Grenzschicht weniger ins Gewicht. Lai und So zeigen, dass mit ihrem Modell die Wie-
dergabe der Temperatur nahe der Wand prinzipiell sehr gut ist. Allerdings wird in deren
Arbeit das Modell für eine voll aufgelöste Wandgrenzschicht verwendet. Eventuell hat
der hier verwendete High-Reynolds Ansatz, also mit der Überbrückung des wandnahen
Bereichs in diesem Gebiet, Defizite. Während die Kaltrechnungen für die Leistungsbe-
wertung der betrachteten Düse eher akademische Rechenbeispiele darstellen, sind die
Heißrechnungen um TR=2,5 hingegen am interessantesten. Für diese realitätsnäheren
74
Abbildung 4.18.: Geschwindigkeitskoeffizient über TR.
Temperaturverhältnisse scheint das DSM hingegen verlässlichere Ergebnisse zu liefern.
Rechnung P rT =0,9 P rT =0,5 DSM

1.6-1.1-1.0 0,033% 0,032% 0,150%
1.6-1.1-1.5 0,035% 0,034% 0,148%
1.6-1.1-2.0 0,031% 0,031% 0,094%
1.6-1.1-2.4 0,029% 0,028% 0,043%
2.6-1.1-1.0 0,037% 0,034% 0,159%
2.6-1.1-1.5 0,042% 0,038% 0,157%
2.6-1.1-2.0 0,044% 0,042% 0,117%
2.6-1.1-2.4 0,041% 0,039% 0,081%
2.6-1.1-2.55 0,040% 0,038% 0,065%
Tabelle 4.2.: Totalenthalpieflussdifferenz zwischen Fan- und Kernstromeintritt sowie

Duesenausgang.
75
Abbildung 4.19.: Geschwindigkeitskoeffizient über TR.
4.3.3. Zusammenfassung - Verbessertes turbulentes

Wärmestrommodell
Über die direkte Formulierung des turbulenten Wärmestromterms als Differentialglei-

chung kann das Dilemma der Findung einer geeigneten turbulenten Prandtlzahl umgan-
gen werden. Außerdem kann bei korrekter Formulierung aller Schließungsansätze davon
ausgegangen werden, dass der turbulente Wärmestromvektor besser wiedergegeben wird
als mit dem einfachen Eddy Diffusivity Modell. Während beim EDM nur der Tem-
peraturgradient die Richtung des Wärmestroms bestimmt, geht beim DSM auch der
Geschwindigkeitsgradient in die Berechnung mit ein. Dies bewirkt eine ebenfalls starke
Abhängigkeit des Wärmestromvektors vom Geschwindigkeitsfeld. In den Konturplotts
für die Temperaturverteilung wurde gezeigt, wie das DSM Modell im Vergleich zum
EDM in verschiedenen Regionen im Rechengebiet unterschiedlich stark mischt. Beson-
ders im turbulenten Nachlauf des Mischers, wo die Strömung stark dreidimensional ist,
mischt das differentielle Modell stärker als das EDM mit einer turbulenten Prandtlzahl
von 0,5. Bei den Heißrechnungen um TR=2,5 liegen die Durchflusskoeffizienten für das
76
Abbildung 4.20.: Wandtemperaturverteilung der OT-Düse für den Fall 1.6-1.1-2.4.

Oben: Fanstrombereich. Unten: Kernstrombereich.
Jeweils links: P rT = 0, 9. Jeweils rechts: DSM.
DSM auf etwa demselben Niveau wie für P rT = 0, 5. Die Geschwindigkeitskoeffizienten

sind beim DSM etwas höher. Allerdings fällt auf, dass für geringere TR beim DSM ein
gewisser Versatz zu den Kurven des EDM zu erkennen ist. Erklärbar wäre dieses Verhal-
ten durch die Beobachtung, dass das differentielle Modell etwas schlechter konvergiert
als der algebraische Ansatz bzw. eine größere Totalenthalpieflussdifferenz zwischen Fan-
und Kernstromeinlass und dem Düsenausgang aufweist. Grund hierfür ist möglicherwei-
se die stärkere Erwärmung der Wandgrenzschicht beim DSM.
In Kapitel 4.2 wurde erläutert wie sich Temperaturturbulenz auf die Strömung auswirkt.
Es wurde erörtert, dass sich eine Berücksichtigung von thermischer Turbulenz auch auf
die Impulsgleichungen auswirken sollte. Dies ist bei dem hier untersuchten Modell nicht
der Fall. Der Einfluss beschränkt sich in erster Linie auf die Energiegleichung. Eine Kop-
pelung mit einem Temperaturkorrekturmodell wie in Kapitel 4.2 könnte daher sinnvoll
sein und noch Verbesserungspotential bieten. Diese Idee soll im nachfolgenden Kapitel
verfolgt werden.
Teile der Ergebnisse in diesem Kapitel wurden auf der 21. AIAA Computational Fluid
77
Dynamics Conference in San Diego, USA, veröffentlicht [129].
4.4. Koppelung von Temperaturkorrektur und

verbessertem Wärmestromterm
In diesem Teil sollen die in Abschnitt 4.2 vorgestellte Temperaturkorrekturmethode und
das in 4.3 eingesetzte differentielle Modell für den turbulenten Wärmestrom kombiniert
werden. Zum Vergleich und zur Kalibrierung wurden von Rolls-Royce Deutschland Mess-
daten zur Verfügung gestellt. Es sind sowohl Messdaten für die Totaltemperatur als auch
für den Totaldruck an zwei Ebenen verfügbar. Ziel ist es, das Temperaturkorrekturverfah-
ren und das differentielle Wärmestrommodell so zu koppeln, dass sich eine bestmögliche
Übereinstimmung mit den Totaltemperaturprofilen aus den Messungen ergibt.
4.4.1. Experimentelle Messungen
Die von Rolls-Royce Deutschland zur Verfügung gestellten experimentellen Daten sind
Messungen eines gemischten Rolls-Royce Düsensystems. Die Messungen wurden für den
in 4.1.2 beschriebenen Sideline Fall durchgeführt, und wurden bei QinetiQ in England
vorgenommen.
Die experimentellen Messungen erfolgten an einem Düsenmodell, welches wie das ent-
sprechende CFD-Modell (RR-Düse) einen Maßstab von ca. 1:5 gegenüber der Original-
düse aufweist. Gegenüber der simulierten Geometrie ist die Düse im Experiment leicht
geknickt. Dies kann später anhand der Konturplotts gut erkannt werden. So wandern
das Totaldruckprofil und das Totaltemperaturfeld in den Darstellungen mit zunehmen-
dem Abstand vom Auslass in Richtung der oberen rechten Bildecke. Für die Messungen
wurde je eine Traverse mit 12 Messsensoren verwendet. In radialer Richtung besitzen die
Sensoren je einen Abstand von 10 mm zueinander. Um das komplette Feld zu erfassen
wurde der Messrechen in Intervallen von je 2◦ gedreht. Gemessen wurde an zwei Posi-
tionen in axialer Richtung. Einmal 5 mm hinter dem Düsenausgang, das zweite Profil
wurde mit Abstand von einem Düsendurchmesser hinter der Düsenauslassebene gene-
riert.
Die Messdaten wurden zuerst mit der kommerziellen Software Tecplot 360 2012R1 aus-
gewertet und anschließend mit ParaView visualisiert. Da nichts genaueres bekannt ist
wird davon ausgegangen, dass die Datenfelder jeweils die Profile wie üblich mit Blick in
Richtung in die Düse zeigen.
78
Totaltemperaturmessung
Die Messung des Totaltemperaturfelds ist in ein Datenfeld für den Kernstrom und den
Bypassstrom unterteilt. Beide Felder wurden in Tecplot zu einem Feld zusammengefasst.
Das Netz der Messsensoren ist in Abbildung 4.21 dargestellt. Die Datenpunkte sind für
die jeweiligen Gitterknoten erhältlich.
Abbildung 4.21.: Messgitter für die Totaltemperaturmessung. Das Netz ist unterteilt in
einen Bypass- und Kernstrombereich.
Die achsennahen Sensoren befinden sich bei +5 mm und -5 mm in radialer Richtung. Bei
einer 360◦ Rotation der Traverse wird das Feld in Achsennähe also überlappend bzw.
doppelt gemessen. Für die Auswertung wurden deshalb für die Werte bei +5 und -5 mm
Mittelwerte berechnet. Im zur Verfügung gestellten Datenfeld für x/D=1,0 (x/D=0 ist
Düsenauslassebene) waren manche Temperaturmessungen offensichtlich fehlerhaft und
weisen Werte weit unterhalb von 0◦ C auf. Diese fehlerhaften Messpunkte konnten teil-
weise durch interpolierte Werte umliegender Messwerte ersetzt werden. Die Messung 5
mm hinter dem Düsenauslass ist in Abbildung 4.22 dargestellt. Die Daten für x/D=1,0
sind in Bild 4.23 zu sehen. Weiß eingerahmt sind hier bereits die grafischen Vergleichs-
flächen für die CFD. Diese entsprechen in etwa einem numerisch simulierten Segment
der Düse.
79
Abbildung 4.22.: Experimentelle Totaltemperaturmessung 5 mm hinter dem Düsen-

austritt. Die weiß eingerahmte Fläche entspricht in etwa dem numerisch
berechneten Düsensegment.
Abbildung 4.23.: Experimentelle Totaltemperaturmessung bei x/D=1,0. Die weiß einge-

rahmte Fläche entspricht in etwa dem numerisch berechneten Segment.
80
Totaldruckmessung
Das Messgitter für den Totaldruck ist in Darstellung 4.24 abgebildet. Hier befinden
Abbildung 4.24.: Messgitter für die Totaldruckmessung.
sich die Sensoren nahe der Achse bei -10 mm, 0 mm und +10 mm. Da das volle 360◦
Feld gemessen wurde, überlappen sich die Messpunkte auch hier in der Mitte. Für die
grafische Auswertung wurden daher die Messwerte bei -10 mm ignoriert. Somit reduziert
sich die Anzahl der Datenpunkte in radialer Richtung auf elf Punkte. Ferner weisen die
übermittelten Werte zum Teil unphysikalisch hohe Werte auf. Dies betrifft jeweils den
neunten Messfühler in radialer Richtung (aufsteigend und ausgehend vom Ursprung in
der Achsenmitte). Diese Werte wurden alle ersetzt, indem zwischen dem jeweils weiter
innen liegenden und dem außen liegenden Messpunkt der Traverse interpoliert wurde. Die
sich nun ergebenden Totaldruckprofile sind in den Bildern 4.25 und 4.26 dargestellt.
4.4.2. Rechenergebnisse mit existierenden Modellen
Hier werden die numerischen Lösungen für den Sideline-Fall der RR-Düsenströmung
präsentiert.
Zur Berechnung der kinetischen Turbulenz kommt wieder das k-ǫ Turbulenzmodell zum
Einsatz. Um den turbulenten Wärmestrom zu schließen werden die aus Abschnitt 4.2
und 4.3 bekannten und getesteten Modelle verwendet. Sowohl das Strömungsfeld, als
81
Abbildung 4.25.: pt /p∞ aus Experiment 5 mm hinter dem Düsenaustritt. Die weiß
eingerahmte Fläche entspricht in etwa dem numerisch berechneten
Düsensegment.
Abbildung 4.26.: pt /p∞ aus Experiment bei x/D=1,0. Die weiß eingerahmte Fläche ent-
spricht in etwa dem numerisch berechneten Segment.
82
auch die kinetischen und thermischen Modelle werden vollständig mit zweiter Ordnung
im Raum diskretisiert. Auch sonst werden die bekannten, bereits beschriebenen Ein-
stellungen verwendet. Wie in 4.1.2 erläutert, wird das 45◦ -Segment der RR-Düse mit
periodischen Seitenrandbedingungen gerechnet um der im Experiment getesteten Geo-
metrie zu entsprechen.
Gerechnet wurde mit den Modellen EDM (P rT = 0, 9 und 0, 5), DSM und TC. Die Ver-
teilungen der Totaltemperatur für die Ebene nahe dem Düsenausgang sind in Abbildung
4.27 gegenübergestellt. Die Resultate bei x/D=1,0 sind in Bild 4.28 wiedergegeben. Die
Simulationsergebnisse können jeweils mit den experimentellen Daten in Bild 4.22 und
4.23 verglichen werden.
Abbildung 4.27.: Totaltemperaturverteilung der CFD-Rechnungen für den Sideline Fall

5 mm hinter dem Düsenaustritt. Von links nach rechts: EDM (P rT =
0, 9), EDM (P rT = 0, 5), DSM, TC.
Vergleicht man die Rechnungen untereinander fällt auf, dass die schwächste thermische
Ausmischung für das EDM mit P rT = 0, 9 erzielt wird. Die Ergebnisse für das EDM mit
P rT = 0, 5 und das DSM nahe dem Düsenauslass sehen sehr ähnlich aus. Das Tempera-
turkorrekturmodell von Abdol-Hamid scheint vor allem nahe dem heißen Potentialkern
stark zu mischen. Der Vergleich mit der experimentellen Totaltemperaturmessung of-
fenbart jedoch, dass keine der Simulationen den thermischen Mischungsgrad aus dem
Versuch aufweist. Die im Experiment ermittelte Totaltemperaturkontur scheint deutlich
diffusiver und der Vergleich zeigt, dass in den numerischen Modellen zu geringe ther-
mische Mischung erzielt wird. Deutlich ist der Unterschied auch bei x/D=1,0. Das sich
aus den Messungen ergebende Temperaturbild ist viel stärker vermischt. Temperatur-
spitzen befinden sich hier nur noch im Bereich des heißen Potentialkerns, während in
den Rechnungen auch in den Wirbeln noch hohe Temperaturen herrschen. Abgesehen
vom generell höheren Totaltemperaturbetrag im Vergleich zum Experiment, scheint auch
83
Abbildung 4.28.: Totaltemperaturverteilung der CFD-Rechnungen für den Sideline Fall

bei x/D=1,0. Von links nach rechts: EDM (P rT = 0, 9), EDM (P rT =
0, 5), DSM, TC.
das Strömungsbild nicht ganz dem experimentell gemessenen zu entsprechen. So ist der
Heißstrahl des Wirbelpaars aus den beiden simulierten Mischerblüten bei x/D=1,0 im
Nachlauf des Mischers noch deutlich abgrenzbar zu erkennen. Im Experiment ergibt sich
hier ein anderes Bild.
Die Abweichung der Totaldruckprofile zwischen CFD und Experiment sind hingegen
geringer. Die Konturen sind in den Abbildungen 4.29 und 4.30 gegeben. Vor allem die
Rechnung mit Temperaturkorrekturverfahren, bei dem das Modell auch Einfluss auf die
Impulsgleichungen hat, kommt den Messungen am nächsten.
Abbildung 4.29.: pt /p∞ aus CFD-Rechnungen für den Sideline Fall 5 mm hinter dem
Düsenaustritt. Von links nach rechts: EDM (P rT = 0, 9), EDM (P rT =
0, 5), DSM, TC.
84
Abbildung 4.30.: pt /p∞ aus CFD-Rechnungen für den Sideline Fall bei x/D=1,0. Von
links nach rechts: EDM (P rT = 0, 9), EDM (P rT = 0, 5), DSM, TC.
Aufgrund der festgestellten Diskrepanz zwischen den experimentellen Messungen und

den CFD Ergebnissen soll nun versucht werden ein verbessertes Mischungsmodell zu
entwickeln. Das neue Modell soll vor allem mit Hinblick auf die thermische Ausmischung
besser mit dem Experiment korrelieren.
4.4.3. Koppelungsmodell
Im vorangegangenen Abschnitt 4.4.2 wurden die Unterschiede zwischen den Messungen

und Rechnungen mit bereits existierenden Modellen dargelegt. Mittels eines neuen Mo-
dells zur besseren thermischen Ausmischung sollen die Simulationsergebnisse nun dem
Experiment angenähert werden. Dafür werden das Temperaturkorrekturmodell aus Ab-
schnitt 4.2 und das differentielle Wärmetransportmodell aus 4.3 miteinander gekoppelt
und ein neues Mischungsmodell ausgearbeitet. Die Modellkonstanten des neuen Modells
werden dann gemäß den Totaltemperaturmessungen von QinetiQ kalibriert.
Aufstellen des Modells
Die Formulierung für den turbulenten Wärmestrom wird gemäß Gleichung (4.19) über-
nommen. Lediglich die Koeffizienten sollen hier angepasst werden. Der Ansatz von
Abdol-Hamid für das TC wird für die Koppelung leicht modifiziert und neue Konstanten
werden eingeführt. Die Berechnung des Faktors CT aus Gleichung (4.2) wird geändert.
Als neue Gleichung wird folgende Formulierung verwendet:
85

CT = 1 + CT g · TgT gExp . (4.21)
Im Gegensatz zur bisherigen Definition von CT wird hier auf die Kompressibilitätsfunk-
tion f (M aT ) aus Gleichung (4.4) verzichtet. Grund ist, dass diese Funktion im rele-
vanten Mischungsgebiet keine Rolle spielt. Abbildung 4.31 zeigt, dass die Funktion nur
im direkten Nachlauf der Düsenhinterkante aktiv ist. Der in Gleichung (4.2) im Nenner
Abbildung 4.31.: Darstellung der Heaviside Sprungfunktion aus Gleichung (4.4). Simu-
lation des Sideline Falls mit RR-Düse.
auftretende Faktor 0,041 wird in den Zähler verschoben und durch die neue Modellkon-
stante CT g ersetzt. Ferner wird der Exponent in Gleichung (4.2) ebenfalls durch einen
variablen Parameter T gExp ersetzt. Die übrigen Größen Cµ und Tg berechnen sich wie
in den Gleichungen (4.1) und (4.3) bereits definiert.
Das neue gekoppelte Modell soll fortan als TC+DSM bezeichnet werden.
Variation der Modellkonstanten
Mit Gleichung (4.21) und (4.19) ergeben sich nun insgesamt fünf Variablen zur Kali-
brierung des neuen Modells: T gExp , CT g , Csθ , c1θ und c2θ . Zur Unterscheidung der ein-
zelnen Rechnungen wird für die Sideline Rechnung folgende Nomenklatur festgelegt:
SL-<Modell>-<T gExp >-<CT g >-<c1θ >-<c2θ >-<Csθ >. Alle Einstellungen bzgl. der Si-
mulation sind analog zu den Rechnungen in Abschnitt 4.4.2.
In zahlreichen Simulationen wurden diese Parameter nun systematisch verändert und es
wurde versucht, dem Totaltemperaturprofil aus den experimentellen Messungen näher
86
zu kommen. Bereits bei vorhergehenden Testrechnungen hatte sich herausgestellt, dass

unter den drei Modellkonstanten des DSM der Faktor c1θ den größten Einfluss auf den
thermischen Transport und somit auf die Temperaturverteilung hat. Deshalb wurden
hauptsächlich die Konstanten T gExp , CT g und c1θ variiert. Die Simulationen wurden in
vier Gruppen eingeteilt und entsprechend dem Exponenten T gExp kategorisiert. In jeder
Gruppe wurden dann die Faktoren CT g und c1θ durchvariiert. Für jede Kategorie wurde
jeweils eine Startkonfiguration definiert. Dabei wurde der Faktor CT g je so gewählt, dass
sich für den Zahlenwert des berechneten CT im Strömungsfeld eine ähnliche Größenord-
nung ergibt. Für die Koeffizienten der differentiellen Wärmestromgleichung wurden die
aus Kapitel 4.3 bekannten Standardkoeffizienten gewählt. Das Ergebnis für diese jewei-
lige Startkonfiguration kann in Abbildung 4.32 gesehen werden. Während der Wert für
CT bei allen Rechnungen ähnlich groß ist, kann hier der Einfluss des Exponenten T gExp
gut betrachtet werden. Eine Verringerung dieser Größe führt dazu, dass auch kleinere
Werte für den normierten Totaltemperaturgradienten Tg stärker gewichtet werden. Für
T gExp = 3, 0 sind die Übergänge von niedrigen zu hohen CT schärfer, während die Ver-
teilung für Werte von T gExp nahe eins breiter ausfällt.
Um die thermische Ausmischung entsprechend zu beeinflussen, mussten die Werte für
den Faktor CT g zum Teil deutlich erhöht und bis zum mehr als des achtfachens gegenüber
dem jeweiligen Startwert eingesetzt werden. Noch höhere Werte für CT g führten aller-
dings zu instabilen Lösungen oder sogar zum Absturz der Rechnung. Um die Rechnungen
trotz der hohen Werte für CT g etwas stabiler zu machen, wurde der berechnete Wert für
CT auf maximal 15 beschränkt. Dieses Limit wurde für hohe CT g in den meisten Fällen
bei der auskonvergierten Rechnung fast ausschließlich direkt an der Mischerhinterkante
erreicht, wo hohe Totaltemperaturgradienten auftreten. Die Limitierung ist aber auch
sinnvoll speziell beim Anrechnen eines jeden neuen Falles. Hier können gerade bei den
ersten Iterationen sehr hohe CT im Rechengebiet entstehen, bis die Werte auf Grund
der stärkeren Spreizung des Strahls und der Temperaturverteilung wieder sinken. Über
den Faktor c1θ konnte der turbulente Wärmestrom noch einmal direkt stark beeinflusst
werden, während die Größen c2θ und Csθ wie bereits erwähnt einen eher geringen Einfluss
haben.
4.4.4. Ergebnisse
Es sei an dieser Stelle vorweggenommen, dass für die getesteten Parametereinstellungen

die im Experiment gemessene Totaltemperaturverteilung nicht erreicht werden konnte.
Die für die Verteilung maßgeblichen Parameter CT g und c1θ konnten nicht entsprechend
angepasst werden. Dies liegt zum einen - wie bereits erwähnt - daran, dass Rechnungen
für größe CT g abstürzen, zum anderen führen aber auch zu geringe Werte für c1θ zu
87
Abbildung 4.32.: Berechneter Temperaturkorrekturfaktor CT . Von oben nach unten:

SL-TC+DSM-3,0-24,3-3,0-0,4-0,11;
SL-TC+DSM-2,0-11,0-3,0-0,4-0,11;
SL-TC+DSM-1,5-8,0-3,0-0,4-0,11;
SL-TC+DSM-1,0-5,5-3,0-0,4-0,11
schlechterer Konvergenz mit höheren Residuen. Exemplarisch werden hier zwei Rechen-
ergebnisse mit unterschiedlichem T gExp gezeigt. Sie zeigen die jeweils besten Ergebnisse,
88
welche für T gExp = 3, 0 und 1, 5 möglich waren. Für Letzteres wurde die insgesamt
beste Gesamtübereinstimmung mit den experimentellen Totaltemperaturmessungen al-
ler Kalibrierungsrechnungen erzielt. Die Plots für diese beiden Fälle sind in Abbildung
4.33 dargestellt. Während die diffusive Totaltemperaturverteilung aus dem Experiment
nicht erreicht wurde, konnte dennoch eine deutlich stärkere thermische Ausmischung
gegenüber den Modellen aus Teilkapitel 4.4.2 erreicht werden. So sind z. B. die Tem-
Abbildung 4.33.: Totaltemperaturverteilung von ausgewählten Kalibrierungsrechnung-

en. Oben: 5 mm hinter Düsenauslass. Unten: x/D=1,0.
Links: SL-TC+DSM-3,0-200,0-2,0-0,4-0,11.
Rechts SL-TC+DSM-1,5-40,0-1,5-0,4-0,11.
peraturspitzen in den Mischerwirbeln sowie auch die Temperaturen nahe dem heißen
Potentialkern geringer, aber gegenüber dem Experiment immer noch zu hoch. Abbil-
dung 4.34 zeigt die sich aus den Konfigurationen ergebenden Totaldruckprofile an den
jeweiligen Ebenen.
Trotz der vorgenommenen Anstrengungen scheint es offensichtlich nicht möglich mit dem
gewählten numerischen Verfahren die Totaltemperaturmessungen zu reproduzieren. Ver-
gleicht man die Standardmodelle mit den modifizierten Modellen aus Kapitel 4.2 (TC)
89
Abbildung 4.34.: pt /p∞ von ausgewählten Kalibrierungsrechnungen. Oben: 5 mm hinter

Düsenauslass. Unten: x/D=1,0.
Links: SL-TC+DSM-3,0-200,0-2,0-0,4-0,11.
Rechts SL-TC+DSM-1,5-40,0-1,5-0,4-0,11.
und 4.3 (DSM), führen letztere erwartungsgemäß zu stärkerer Ausmischung zwischen

Heiß- und Kaltstrahl. Das TC hat darüber hinaus sogar noch Einfluss auf die Impulsglei-
chungen. Die Begründung des Einflusses eines Temperatur- oder hier Totaltemperatur-
gradienten auf die Impulsgleichungen wurde in dieser Arbeit bereits mehrmals genannt
und scheint logisch. Bei dem in diesem Kapitel dargestellten neuen Mischungsmodell
(TC+DSM) wurde versucht über die Anpassung der Koeffizienten eine stärkere thermi-
sche Ausmischung und Anpassung an das Experiment zu erzwingen. Dafür wurden die
Koeffizienten zum Teil deutlich gegenüber ihren Standardeinstellungen verändert. Der
aus dem Temperaturkorrekturverfahren von Abdol-Hamid stammende Faktor analog
zu CT g wurde dabei um bis zum über achtfachen gegenüber der Standardkonfiguration
erhöht. Es stellt sich daher die Frage, ob eine weitere Erhöhung für noch stärkere Aus-
mischung überhaupt physikalisch ist. Schließlich geht die Höhe des Faktors CT g direkt
einher mit der Größe des Einflusses der Temperaturturbulenz auf die Strömung. Größere
Werte konnten hier in den Rechnungen allerdings auch nicht getestet werden, ohne dass
90
dies einen Absturz der Rechnung zur Folge hat.

Es sollte auch in Betracht gezogen werden, die Totaltemperaturmessungen zu hinterfra-
gen. Schließlich scheinen die experimentellen Totaldruckprofile prinzipiell besser mit den
Simulationen zu korrelieren. Es wäre z. B. interessant zu wissen, ob ein dichteres Netz an
Messsonden im Strahl zu einer besseren Auflösung und Wiedergabe des Totaltempera-
turfelds führen würde. Zum Zweck besserer Vergleichbarkeit wurde daher exemplarisch
eine CFD Rechnung auf das Messgitter interpoliert. Hierfür wurde aus der Rechnung mit
dem Temperaturkorrekturverfahren (siehe Abbildung 4.27, bzw. 4.29) an der Messebene
nahe dem Düsenauslass für die Totaltemperatur und den Totaldruck eine 360◦ Kontur
erstellt und auf das Netz aus den Messungen übertragen. Das Ergebnis ist in Abbil-
dung 4.35 und 4.36 dargestellt. Gegenüber den ursprünglichen Darstellung in 4.27 und
Abbildung 4.35.: Totaltemperaturprofil aus CFD Simulation mit Temperaturkorrektur-

verfahren 5 mm hinter dem Düsenaustritt. Lösung wurde zur 360◦ Kon-
tur vervollständigt und auf das Messgitter interpoliert.
4.29 fällt auf, dass das Ergebnis nach der Übertragung auf das gröbere Messgitter dif-
fusiver wird. Zwar weicht die Totaltemperaturverteilung im Vergleich zum Experiment
immer noch stark ab, die Totaldruckverteilung passt dadurch allerdings nochmals etwas
besser zu den Messungen. Eine zweite Hypothese, wie es bei den Messungen zu einem
solch diffusiven Temperaturfeld kommen kann wäre, dass eine Temperaturübertragung
innerhalb der Traverse auf die einzelnen Messsonden durch Wärmeleitung stattgefunden
haben könnte. Durch solch einen Vorgang würde wahrscheinlich eine diffusivere Total-
91
Abbildung 4.36.: pt /p∞ aus CFD Simulation mit Temperaturkorrekturverfahren 5 mm

hinter dem Düsenaustritt. Lösung wurde zur 360◦ Kontur vervollstän-
digt und auf das Messgitter interpoliert.
temperaturverteilung mit geringeren Maximalwerten im Heißstrahlbereich gemessen, als

sie im Experiment tatsächlich vorhanden ist.
An dieser Stelle soll ein interessanter Vergleich eingeschoben werden. Ohne dabei bis ins
Detail zu gehen wird die SL-Rechnung mit der Lösung auf einem unstrukturierten Gitter
verglichen. Das unstrukturierte Netz wurde von Rolls-Royce Deutschland erstellt. Dabei
handelt es sich um eine 180◦ Geometrie der RR-Düse mit etwas früher abgeschnittenem
Bullet und einer wie im Experiment geknickten Düse. Das Fernfeld wurde deutlich grö-
ßer ausgeführt als in den hier verwendeten strukturierten Gittern. Bild 4.37 zeigt einen
Längsschnitt durch die Geometrie im düsennahen Bereich. Die Rechnung des Sideline
Falls wurde in FLUENT mit dem gekoppelten druckbasierten Löser durchgeführt. Für
die Berechnung des turbulenten Wärmestrom dient das EDM mit P rT = 0, 85. Wie in
den bisher gezeigten Rechnungen wurde das k-ǫ Turbulenzmodell verwendet, auch al-
le anderen Einstellungen sind mit den Simulationen in OpenFOAM vergleichbar. Das
Ergebnis für die Totaltemperaturprofile in Abbildung 4.38 zeigt, dass das verwendete
unstrukturierte Netz zu einer deutlich diffusiveren Verteilung der führt. Diese Lösung
entspricht insgesamt deutlich besser der im Experiment gemessenen Totaltemperatur-
kontur, wobei jetzt im Bereich der Mischerwirbel die Verteilung noch diffusiver ausfällt
92
Abbildung 4.37.: Längsschnitt durch das unstrukturierte Netz der RR-Düse im Düsen-
nahen Bereich.
Abbildung 4.38.: Lösung für das Totaltemperaturfeld auf dem unstrukturierten Netz der
RR-Düse. Links: 5 mm hinter dem Düsenausgang, rechts bei x/D=1,0.
als im Experiment. In der Nähe des heißen Potentialkerns wird die Ausmischung aus dem
Experiment allerdings immer noch nicht erreicht. Die Profile für den Totaldruck werden
in Bild 4.39 gezeigt. Auch die Totaldruckverteilung zeigt sich gegenüber der Messung
verstreuter.
Mit unstrukturierten Netzen bestehend aus Triangeln, Tetraedern oder Pyramiden ist
es in der Regel schwerer den selben Grad an Genauigkeit zu erzielen wie für ein struktu-
riertes Hexaedernetz. Grob gesagt resultiert die schlechtere Genauigkeit unstrukturierter
93
Abbildung 4.39.: Lösung für pt /p∞ auf dem unstrukturierten Netz der RR-Düse. Links:
5 mm hinter dem Düsenausgang, rechts bei x/D=1,0.
Netze aus dem größeren numerischen Abbruchfehler bedingt durch deren Topologie. Be-
sonders gravierend wird der Nachteil dieser Zellen gegenüber Hexaedern in Bereichen
scharfer Gradienten wie in Wandgrenzschichten oder ausgeprägten Scherschichten [130].
In stark dreidimensionalen Strömungsgebieten können allerdings auch strukturiert ver-
netzte Geometrien an Rechengenauigkeit verlieren, nämlich wenn die Zellflächen (Faces)
der Hexaeder nicht mit der Ausbreitungs-/Strömungsrichtung fluchten. Dann entsteht
wieder ein zusätzlicher Abbruchfehler. Die in Bild 4.38 und 4.39 zu sehende starke Diffusi-
vität der Konturplots kann hier also mit den schlechteren Eigenschaften unstrukturierter
Netze in Verbindung gebracht werden, weil die Scherschicht nicht ausreichend aufgelöst
wird. Die Lösung für die hier verwendeten strukturierten Netze sollte deshalb genauer
sein (gleiche Zellgrößen vorausgesetzt, was hier in etwa der Fall ist), ist aber tatsäch-
lich weiter entfernt vom Experiment. Zumindest für die Verteilung der Totaltemperatur.
Auch dieser eingeschobene Vergleich legt nahe, dass die Temperaturmessungen eventuell
fehlerbehaftet sein könnten.
Um nach diesem Einschub wieder zurück zu den Simulationen in OpenFOAM zu kommen
werden jetzt die Auswirkungen der modifizierten Modelle noch quantitativ eingeordnet
und für die durchgeführten Sideline Rechnungen noch die Leistungsparameter berechnet.
Diese sind in Abbildung 4.40 dargestellt. Beim gekoppelten Modell (TC+DSM) führt die
gesteigerte thermische Ausmischung gegenüber den vier bereits bekannten Mischungs-
modellen zu einem deutlich geringeren CD , aber zu einem höheren CV . Die Ergebnisse
passen damit weiterhin in die aus den Kapiteln 4.2 und 4.3 beobachtete Logik, dass
eine höhere thermische Ausmischung zu einem höheren CV , aber zu einem geringeren
CD führt. Diese Beobachtung ist vor allem deshalb interessant, weil die hohen Werte
für CT g direkt zu mehr turbulenter Viskosität innerhalb des Kontrollvolumens der Düse
führen, wie in Abbildung 4.41 gezeigt, aber der Geschwindigkeitskoeffizient trotzdem
94
stark ansteigt.
Abbildung 4.40.: Durchfluss- und Geschwindigkeitskoeffizient für Sideline Rechnung mit

verschiedenen Modellen.
95
Abbildung 4.41.: Turbulente Viskosität für den Sideline Fall. Von oben nach unten:
x/D=-1,5; -0,75; 0,0 (x/D=0,0 ist Düsenaustrittsfläche). Von links nach
rechts: EDM (P rT = 0, 9), TC, TC+DSM-3,0-200,0-2,0-0,4-0,11,
TC+DSM-1,5-40,0-1,5-0,4-0,11.
4.4.5. Zusammenfassung - gekoppeltes Mischungsmodell
Bei einem Vergleich zwischen experimentellen Messungen sowie Rechenergebnissen und

den in vorherigen Abschnitten bereits gezeigten Modellen (EDM, TC, DSM) wurde eine
größere Diskrepanz in den Totaltemperaturprofilen festgestellt. Gegenüber der CFD wei-
sen die Messungen eine stärker diffusive Temperaturverteilung auf. Aus diesem Grund
sollte in diesem Abschnitt das aus Kapitel 4.2 bekannte Temperaturkorrekturverfahren
von Abdol-Hamid und das in 4.3 angewandte differentielle Wärmestrommodell mitein-
ander gekoppelt und zu einem neuen Mischungsmodell verschmolzen werden. Während
beim DSM nur die Koeffizienten angepasst wurden, wurde die Berechnung des Faktors
CT gegenüber dem TC modifiziert. Das neue Modell sollte dann anhand der experimen-
tellen Messungen, welche von Rolls-Royce Deutschland zur Verfügung gestellt wurden,
96
justiert werden.
Durch zahlreiche Kalibrierungsrechnungen konnte mit dem TC+DSM die thermische
Ausmischung gegenüber bisherigen Rechnungen nochmal deutlich gesteigert werden. Al-
lerdings musste der Faktor CT g dazu sehr stark erhöht werden. Ebenso wird der Faktor
c1θ gegenüber seinem Standardwert aus dem DSM stark verändert. Trotz der großen
Modifikationen konnte aber das stark diffusive Bild der Totaltemperaturverteilung aus
den Messungen numerisch nicht erreicht werden. Größere Werte für CT g führten bei den
Kalibrierungsrechnungen allerdings zu instabilen Lösungen und sogar zur Divergenz der
Rechnung. Die gezeigten Lösungen für das TC+DSM stellen die bestmögliche Überein-
stimmung mit den Experimenten dar.
Während das Temperaturprofil aus den von Rolls-Royce zur Verfügung gestellten Daten
nicht reproduziert werden konnte, war die Übereinstimmung mit den Totaldruckprofilen
zwischen Experiment und Numerik prinzipiell besser. Statt den Grund für die Diskre-
panz der Totaltemperaturprofile ausschließlich in der Numerik zu suchen, könnten daher
auch die Temperaturmessungen in Frage gestellt werden. Zum einen wäre es interessant,
ob ein feineres Netz an Messpunkten ein weniger diffusives Totaltemperaturprofil erge-
ben würde. Zweitens wäre es theoretisch möglich, dass die Messergebnisse aufgrund von
Wärmeleitung innerhalb der Traverse verfälscht wurden. Der eingeschobene Vergleich
mit der Simulation auf einem unstrukturierten Netz zeigt bessere Übereinstimmung mit
der experimentellen Totaltemperaturmessung. Es wird allerdings angenommen, dass die
zusätzliche Diffusivität durch das unstrukturierte Hexaedernetz eingebracht wird und
die Lösung auf dem strukturierten Netz genauer sein sollte. Obwohl ungenauer, passt
also das Totaltemperaturfeld besser, aber das Totaldruckprofil scheint nun zu diffusiv
zu sein. Dies bekräftigt die Vermutung, dass die Totaltemperaturmessung fehlerbehaftet
sein könnte.
Die Auswertungen der Leistungsdaten für den Sideline Fall zeigen jedoch, dass sich die
Ergebnisse für CV und CD mit dem neuen Modell weiterhin an die zuvor beobachtete
Gesetzmäßigkeit halten. Obwohl das neue Modell zum Teil deutlich mehr turbulente
Viskosität im Strömungsfeld erzeugt, steigt der Geschwindigkeitskoeffizient für größere
TR aufgrund der stärkeren thermischen Ausmischung deutlich an, während der Durch-
flusskoeffizient sinkt.
97
5. Parameteruntersuchung der
Düsenströmung mit ANSYS
FLUENT
Neben all den Verbesserungen im Detail besteht grundsätzlich die Unsicherheit welchen
Einfluss bestimmte Parameter und getroffene Annahmen auf die Simulation haben.
In Kapitel 4 wurde gezeigt, dass der charakteristische Verlauf der Leistungsparameter
durch Modelloptimierung im Vergleich zu experimentellen Daten zwar verbessert werden
kann, jedoch ein prinzipieller Versatz zwischen beiden Datenreihen liegt. Eine Diskrepanz
zwischen Experiment und CFD kann grundsätzlich zwar an fehlerhaft erfassten Messda-
ten liegen, häufiger jedoch liegt die Abweichung an den Daten aus der Numerik. Durch
die meist modellierte Physik, wie auch durch die verwendeten numerischen Modellen,
Algorithmen, Diskretisierungsschemata, etc. kommen eine Vielzahl potentieller Fehler-
quellen ins Spiel. Im folgenden Kapitel soll nun versucht werden den Einfluss numerischer
und physikalischer Annahmen für die Berechnung einer kompressiblen Düsenströmung,
im Hinblick auf die Berechnung ihrer Leistungsparameter, zu quantifizieren.
5.1. Rechenmodell und Randbedingungen

Im Gegensatz zu den bisher durchgeführten Rechnungen dient für diese Parameterun-
tersuchung als Strömungslöser der kommerzielle Code FLUENT 14.0 von ANSYS. Als
Lösungsverfahren wurde der druckbasierte Strömungslöser mit Druck- und Geschwin-
digkeitskoppelung verwendet. Die bereits bekannte OT-Düsengeometrie wurde mit 4,3
Mio. Zellen diskretisiert. Für die Gittergenerierung wurde das Netz aus den vorher ge-
zeigten Rechnungen mit 1,25 Mio. Zellen als Basis verwendet und die Zellenzahl in alle
drei Raumrichtungen mit dem Faktor 1,5 expandiert. Der Skalierungsfaktor beträgt wie
bisher ca. 0,2 gegenüber der Originalgröße. Als Turbulenzmodell für die Referenzrech-
nungen dient das in FLUENT implementierte k-ω-SST Modell. Während die Strömungs-
gleichungen mit zweiter Ordnung im Raum diskretisiert wurden, wurde für die Turbulenz
erste Ordnung verwendet. Im Gegensatz zu der bisherigen Annahme einer konstanten
98
Parameteruntersuchung der Düsenströmung mit ANSYS FLUENT
spezifischen Wärmekapazität wird dieses Mal für die Referenzrechnungen ein von der
Temperatur abhängiges Polynom verwendet:
J J J J
cp (T ) = 1035, 887 − T · 0, 256 2
+ T 2 · 0, 625 · 10−3 3
− T 3 · 2, 628 · 10−7 .
kgK kgK kgK kgK 4
(5.1)
Der Verlauf der Funktion ist in Schaubild 5.1 dargestellt.
1150
1100
c [J/kgK]
p
1050
1000
200 300 400 500 600 700 800 900 1000
T [K]
Abbildung 5.1.: Polynom cp (T ).
Sonstige Annahmen und Randbedingungen für den Referenzfall:
• ISA Bedingungen mit pt,∞ = 101325 Pa, T∞ = 288, 15 K

• p∞ = 101148 Pa, M a∞ = 0, 05
• Polynom für Wärmeleitfähigkeit, siehe (3.23), Sutherlandformel für Viskosität
• EDM mit P rT = 0, 85 für turbulenten Wärmestrom
• symmetrische Seitenrandbedingung
99
• Einlassturbulenz: Fan und Core: 5 %; Fernfeld: 0,1 %

• turbulenter Längenmaßstab am Eintritt:
Fan und Core: 5 · 10−5 m; Fernfeld: 0, 0015m
• Gradientenberechnung: Gauss, zellbasiert
Weitere Details zum Löser können dem FLUENT User’s [56] Guide und Theory Guide
[55] entnommen werden.
Tabelle 5.1 gibt eine Übersicht über die simulierten Betriebspunkte und entspricht qua-
si den in Kapitel 4 gerechneten Betriebspunkten für die Leistungsrechnung. Es wurde
aber nicht immer für alle Parameter die volle Matrix in Tabelle 5.1 gerechnet, sondern
manchmal einzelne Punkte ausgewählt. Die Nomenklatur für die Bezeichnung der Be-
triebszustände ist wieder aus Kapitel 4.1.2 zu entnehmen.
CNPR PS TR
1.0
1.5
1.6 1.1 2.0
2.4
1.0
1.5
2.6 1.1 2.0
2.4
Tabelle 5.1.: Simulierte Betriebspunkte der Parameterstudie.
Tabelle 5.2 gibt noch einmal eine Übersicht über die Referenzeinstellungen und zeigt die
variierten Parameter.
100
Parameter Referenz Variation

Gittergröße 4,3 Mio. 1,25 Mio.; 15 Mio.
Seitenrandbedingung
symmetrisch periodisch
(bzw. Mischertaktung)
spez. Wärmekapazität cp = f (T ) cp = 1007 J/kgK
Turbulenzmodell SST 2ndO;
und Diskretisierung SST 1stO SKE 1stO, 2ndO
(mit std. Wandfkt.)
Tabelle 5.2.: Parametervariation: Referenzwerte und variierte Parameter.
5.2. Ergebnisse Parameteruntersuchung
5.2.1. Gittergröße
Für die Gitterstudie wurden drei unterschiedlich feine Netze erzeugt. Als Basis dient
das bereits in den Kapiteln 4.2 und 4.3 verwendete OT-Rechennetz. Von ihm ausgehend
wurde die Anzahl der Zellen jeweils mit einem Expansionsfaktor von 1,5 in alle drei
Raumrichtungen vergrößert. Somit ergeben sich die drei Netze mit je 1,25 Mio., 4,3 Mio.
und 15 Mio. Zellen.
Die Plots der Totaltemperatur sind in Abbildung 5.2 zu sehen. Trotz der großen Dis-
krepanz in der Gitterzahl zwischen dem feinsten und dem gröbsten Gitter sind bei der
qualitativen Betrachtung die Unterschiede bei der Totaltemperaturverteilung überschau-
bar. Bei genauerer Betrachtung ist zu sehen, dass die Verteilung bei 1,25 Mio. Zellen
leicht diffusiver ist als bei den anderen beiden. Die Unterschiede zwischen dem Gitter
mit 4,3 Mio. und 15 Mio. Zellen sind sehr gering.
Deutlicher fällt die Differenz der drei Gitter bei den Leistungsparametern aus, wie in
den Schaubildern 5.3 zu sehen ist. Hier sind die CD und CV Werte über der invertierten
Gitteranzahl geplottet. Sowohl Durchfluss- als auch Geschwindigkeitskoeffizient steigen
für ein feineres Rechengitter an. Die Differenz von CD und CV zum jeweils gröbsten
Gitter ist in Tabelle 5.3 gegeben. Wie aus den Werten ersichtlich ist, scheint es keine
große Rolle zu spielen ob eine Heiß- (TR>1,0) oder Kaltrechnung (TR=1,0) durchge-
führt wird. Die Anstiege der Werte hin zum feineren Netz sind in etwa identisch. D.h.
eine bessere Auflösung des Temperaturgradienten durch das feinere Netz scheint eine
eher untergeordnete Rolle bei der Leistungsrechnung zu spielen.
101
Abbildung 5.2.: Tt -Plot für den Fall 1.6-1.1-2.4. Von links nach rechts: 1,25; 4,3; 15 Mio.
Zellen. Von oben nach unten: x/D=-1,4; -0,7; 0,0; 0,7; 1,4 (x/D=0,0 ist
Düsenaustrittsebene).
102
Abbildung 5.3.: CD und CV über 1/Zellenanzahl.
# Zellen N ∆ CD % ∆ CV % CNPR TR
4,3 Mio 0,13 0,08 1,6 1,0
15,0 Mio 0,18 0,12
4,3 Mio 0,10 0,09 1,6 2,4
15,0 Mio 0,16 0,14
4,3 Mio 0,05 0,06 2,6 1,0
15,0 Mio 0,09 0,08
4,3 Mio 0,06 0,05 2,6 2,4
15,0 Mio 0,10 0,08
Tabelle 5.3.: Differenz von CD and CV verglichen zum kleinsten Netz mit 1,25 Mio Zellen.
5.2.2. Periodische Ränder/Mischertaktung
Durch die Definition einer periodischen Randbedingungen an den Seitenrändern statt

einer Symmetriebedingung lässt sich für die gegebene Geometrie der Einfluss einer geän-
derten Mischertaktung simulieren, ohne tatsächlich die Geometrie bzw. das Rechennetz
zu ändern.
In Bild 5.4 ist deutlich zu erkennen wie sich durch die periodischen Ränder das Wirbel-
bild bzw. die Totaltemperaturverteilung ändert.
Auf die Leistungsdaten scheint die nun geänderte Taktung der Mischerblätter allerdings
keinen Einfluss zu haben. Obwohl das sich nun ergebende Strömungs- bzw. Tempe-
raturfeld ein anderes ist, bleiben die Leistungsparameter dabei nahezu identisch. Die
Graphen sind in Schaubild 5.5 geplottet. In Kapitel 4.2 wurde erwähnt, dass die expe-
rimentelle Bestimmung der Düsenleistungsdaten mit einer anderen Mischertaktung wie
im CFD Modell durchgeführt wurde. Wie hier gezeigt, führt aber eine andere Taktung
nicht zwangsläufig zu geänderten Leistungsdaten. Somit sind die unterschiedlichen Wer-
103
Abbildung 5.4.: Tt -Plot für den Fall 1.6-1.1-2.4. Links: symmetrische Seitenränder,
rechts: periodische Ränder. Von oben nach unten: x/D=-1,4; -0,7; 0,0;
0,7; 1,4 (x/D=0,0 ist Düsenaustrittsebene).
104
te für CD und CV zwischen Experiment und CFD wohl nicht mit der unterschiedlichen
Geometrie zu erklären.
Abbildung 5.5.: Leistungskoeffizienten über TR.
5.2.3. Konstante/variable spezifische Wärmekapazität
Für die Rechnungen in OpenFOAM wurde eine konstante turbulente Wärmekapazität

von cp = 1007 J/kgK verwendet. Daher soll hier der Einfluss einer von der Tempe-
ratur abhängigen Wärmekapazität untersucht werden. FLUENT bietet die Möglichkeit
für cp ein temperaturabhängiges Polynom einzusetzen, was in allen Referenzrechnungen
gemacht wurde. Das Polynom ist in Gleichung (5.1) gegeben. Der Vergleich zeigt den
Unterschied zu dem konstanten Wert von 1007 J/kgK.
Abbildung 5.6 zeigt eine Gegenüberstellung für 1.6-1.1-2.4 und stellt die Totaltempera-
turverteilung dar. Qualitativ sind die Unterschiede zwischen beiden Rechnungen aller-
dings minimal.
Die Leistungsbeiwerte sind jedoch deutlich beeinflusst. Die Kurven divergieren bei höhe-
ren Totaltemperaturverhältnissen auseinander. Die Werte für CD sinken für steigende TR
bei einem temperaturabhängigen cp stärker ab. Der Geschwindigkeitskoeffizient hinge-
gen steigt gleichzeitig stärker an als für ein konstantes cp . Der Grund für die Abweichung
liegt in der Vorgabe der Totaltemperatur an den numerischen Einlässen. Da mit stei-
genden Temperaturen auch die spezifische Wärmekapazität zunimmt (siehe Abbildung
5.1), steigt damit die Enthalpie am Einlass für ein variables cp und führt zu anderen
Leistungswerten. Die Kurven sind in Bild 5.7 dargestellt.
105
Abbildung 5.6.: Totaltemperaturverteilung für 1.6-1.1-2.4. Links: cp = cp (T ); rechts:

cp = konst.. Von oben nach unten: x/D=-1,4; -0,7; 0,0 (x/D=0,0 ist
Düsenaustrittsebene).
Abbildung 5.7.: Leistungsbeiwerte über TR.
106
5.2.4. Turbulenzmodell und dessen Diskretisierung
In Kapitel 4.2 wurde bereits mit unterschiedlichen Turbulenzmodellen gerechnet und

die Ergebnisse verglichen. Dabei wurden in diesem Teil beide Turbulenzmodelle mit ei-
ner Genauigkeit von erster Ordnung im Raum diskretisiert. Hier soll der Vergleich der
Turbulenzmodelle nochmals zusammengefasst und zusätzlich der Einfluss der Diskreti-
sierungsordnung untersucht werden.
Während das k-ω-SST Turbulenzmodell mit räumlicher Diskretisierung erster Ordnung
(1stO) als Referenz dient, wurde der Einfluss der Turbulenz bei zweiter Ordnung (2ndO)
untersucht. Ebenso wurden Rechnungen mit dem standard k-ǫ Modell erster und zweiter
Ordnung durchgeführt. Beim SKE wurde mit Standardwandfunktionen (siehe [55, 56])
gerechnet.
Beim Blick auf die Lösung offenbart sich ein äußerst interessantes Verhalten der Mo-
delle, welches so nicht unbedingt erwartet wurde. Der Plot für die Totaltemperatur ist
in Abbildung 5.8 dargestellt. Am diffusivsten verhält sich offensichtlich das SST Mo-
dell bei erster Ordnung, während das SKE bei erster Ordnung dagegen die geringste
Diffusion der Totaltemperatur aufweist. Bei einer Genauigkeit zweiter Ordnung schei-
nen sich beide Turbulenzmodelle hinsichtlich der thermischen Verteilung anzugleichen.
Die Begründung für diese Totaltemperaturverteilung liefert ein Blick auf die Verteilung
der turbulenten Viskosität in Bild 5.9. Das SST Modell erzeugt bei erster Ordnung die
meiste turbulente Viskosität. Beim SKE Turbulenzmodell ist sie bei dieser Genauigkeit
insgesamt offensichtlich am geringsten. Rechnet man bei beiden Turbulenzmodellen mit
zweiter Ordnung, ergeben sich für die Verteilung der turbulenten Viskosität nahezu iden-
tische Lösungen.
Ähnliches spiegeln auch die Leistungsparameter für die Heißrechnungen wider (siehe
Abbildung 5.10 und 5.11). Bei den Kaltrechnungen sind die Durchflusskoeffizienten fast
unbeeinflusst vom Turbulenzmodell und der Ordnung. Größer ist der Einfluss bei den
Heißrechnungen. Liegen die Werte für CD bei CNPR=1,6 und erster Ordnung noch
weiter auseinander, rücken die Lösungen beider Turbulenzmodelle bei zweiter Ordnung
näher zusammen. Für CNPR=2,6 liefern SKE und SST bei zweiter Ordnung sogar nahe-
zu identische Ergebnisse. Ähnlich verhalten sich die Geschwindigkeitskoeffizienten. Bei
den Kaltrechnungen scheint weder das Turbulenzmodell noch die Ordnung einen Einfluss
zu haben. Bei den Heißrechnungen liefern SST und SKE praktisch identische Werte für
CV , wenn mit zweiter Ordnung gerechnet wird.
Obwohl wie in den Plots in 5.9 gezeigt, die Größenordnung der Verteilung von µT unter-
schiedlich ist, kann aus dem beobachteten Verhalten geschlossen werden, dass die sich im
Strömungsfeld befindliche turbulente Viskosität auf die Leistung der Düse vernachlässig-
baren Einfluss hat. Erst bei den Heißrechnungen kommt die berechnete Wirbelviskosität
entscheidend zum tragen, da diese aufgrund des Schließungsansatzes für den turbulenten
Wärmestrom die Wärmeübertragung vom Heiß- auf den Kaltstrahl maßgeblich beein-
107
Abbildung 5.8.: Tt Verteilung für 1.6-1.1-2.4. Von links nach rechts: SST 1stO, SST
2ndO, SKE 1stO, SKE 2ndO. Von oben nach unten: x/D=-1,4; -0,7;
0,0; 0,7; 1,4 (x/D=0,0 ist Düsenaustrittsebene).
108
Abbildung 5.9.: Verteilung von µT für 1.6-1.1-2.4. Von links nach rechts: SST 1stO, SST
2ndO, SKE 1stO, SKE 2ndO. Von oben nach unten: x/D=-1,4; -0,7; 0,0;
0,7; 1,4 (x/D=0,0 ist Düsenaustrittsebene).
109
Abbildung 5.10.: CD über CNPR.
Abbildung 5.11.: CV über CNPR.
flusst.
Des Weiteren muss aufgrund des hier gezeigten Resultats das Ergebnis aus Kapitel 4.2
differenzierter bewertet werden. Dort wurde festgestellt, dass das SST gegenüber dem
SKE Turbulenzmodell stärker ausmischt, was sich bei den Heißrechnungen auf die Leis-
tungsbeiwerte auswirkt. Der Mischungseffekt hatte sich dann, wie in 4.2 ebenfalls gezeigt,
mit dem Temperaturkorrekturmodell nochmals verstärkt. Allerdings wurden in allen
Rechnungen die Turbulenzmodelle mit erster Ordnung diskretisiert. Mit dem hier gefun-
denen Zusammenhang zwischen Turbulenzmodell und der Ordnung kann man schlussfol-
gern, dass die Ergebnisse aus Kapitel 4.2 bei zweiter Ordnung für beide Turbulenzmodelle
näher zusammen liegen würden. Gleiches gilt für das Temperaturkorrekturmodell.
110
5.3. Zusammenfassung - Parameteruntersuchung

In der Parameteruntersuchung wurde der Einfluss physikalischer als auch numerischer
Größen auf die Rechnungen untersucht. Hier werden die wichtigsten Erkenntnisse zu-
sammengefasst.
Bei den drei unterschiedlich feinen Netzen wurden mit dem feinsten Gitter die höchsten
Werte für CD und CV erzielt. Extrapoliert man die Werte der Leistungsparameter aus
Schaubild 5.3 hin zur Ordinate, scheinen für ein noch feineres Netz als eines mit 15 Mio.
Zellen nochmals leicht höhere Werte für die Leistungskoeffizienten möglich. Qualitativ
bleibt der Einfluss der Gittergröße allerdings überschaubar.
Eine geänderte Mischerblütentaktung wurde durch die Definition einer periodischen
Randbedingung an den Seitenrändern simuliert. Während das Strömungsbild zu sym-
metrischen Rändern deutliche Unterschiede aufweist, sind CD und CV hingegen quasi
unbeeinflusst. Die in Kapitel 4.2 abweichende Mischertaktung zwischen CFD und Ex-
periment sollte folglich keinen Einfluss auf die Leistungsparameter haben. Somit sollten
beide, CFD und Messungen diesbezüglich vergleichbar sein.
Wählt man für die spezifische Wärmekapazität ein temperaturabhängiges Polynom statt
einem konstanten Wert von 1007 J/kgK führt das bei Heißrechnungen zu höheren Er-
gebnissen für CV . CD hingegen sinkt für steigende TR. Die Unterschiede stehen mit
einer höheren Totalenthalpie an den numerischen Einlässen in Zusammenhang, da hier
als Randbedingung die Totaltemperatur vorgegeben wird.
Interessantes gibt der Vergleich des SKE mit dem SST Modell mit unterschiedlicher
Diskretisierungsordnung preis. Wie bereits in Kapitel 4.2 festgestellt, mischt das k-ω-
SST Modell bei erster Ordnung besser als das k-ǫ Turbulenzmodell. Führt man dieselbe
Rechnung bei zweiter Ordnung durch, liefern beide Turbulenzmodelle fast identische
Ergebnisse. D.h. dass die in Kapitel 4.2 durchgeführten Vergleichsrechnungen für das
SST und SKE Modell für größere TR bei zweiter Ordnung prinzipiell näher beieinander
liegen sollten.
111
6. Zusammenfassung
Der Fokus dieser Arbeit lag auf der Simulation und Untersuchung des Mischungspro-
zesses des heißen Kern- und des kälteren Nebenstrom. In einem ersten Schritt wurden
die Defizite von bisher gängigen Berechnungsmodellen bei solchen stark anisothermen
Strömungen lokalisiert. Um die Qualität der Simulationsergebnisse zu verbessern wur-
den dann in den Open Source Strömungslöser OpenFOAM geeignete Modelle imple-
mentiert und mit den Ergebnissen bestehender Verfahren verglichen. Zum einen wurde
ein Temperaturkorrekturmodell verwendet um den Effekt von Temperaturturbulenz zu
berücksichtigen. Des Weiteren wurde ein differentieller Ansatz zur Berechnung des tur-
bulenten Wärmestroms verwendet. Letztlich wurde versucht eine neues Modell zu entwi-
ckeln, dass eine Kombination des Temperaturkorrekturverfahrens und der differentiellen
Wärmestromvorhersage darstellt. Zur Abrundung dieser Arbeit wurde der kommerzielle
Löser FLUENT verwendet um den Einfluss unterschiedlicher numerischer und physi-
kalischer Parameter zu untersuchen. Die Lösungen wurden jeweils qualitativ als auch
quantitativ mittels der Definition von Leistungsparametern bewertet und zum Teil mit
experimentellen Messungen verglichen.
Aufgrund der für kompressible Strömungen durchgeführten Favre Mittelung der Glei-
chungen werden Dichtefluktuationen eliminiert. Bei vielen Ingenieursanwendungen, z. B.
wenn isotherme Strömungen vorherrschen, kann die Vernachlässigung von thermodyna-
mischer Turbulenz eine akzeptable Annahme sein. Bei stark anistothermen Strömungen
wie z. B. bei der Ausmischung von heißen und kalten Abgasstrahlen können diese phy-
sikalischen Effekte allerdings einen deutlichen Einfluss auf das Strömungsbild haben.
In der Realität führen stark anisotherme Scherschichten zu einer höheren Turbulenz.
Übertragen auf das Beispiel der Abgasstrahlmischung haben die höheren Instabilitä-
ten eine verstärkte Ausmischung zur Folge, die aber von den gängigen Modellen in
den Simulationen bisher nicht berücksichtigt werden. Durch die Implementierung ei-
nes Temperaturkorrekturmodells sollen die Effekte von Temperaturturbulenz modelliert
werden. Dieses für das k-ǫ Modell kalibrierte Verfahren wurde zusätzlich für das k-ω-
SST Turbulenzmodell von Menter erweitert. Es konnte gezeigt werden, dass durch das
Temperaturkorrekturverfahren die Ausmischung zwischen dem heißen Kernstrom und
dem kalten Fanstrahl prinzipiell gesteigert wird. Die Folge ist eine homogenere Tem-
peraturverteilung am Düsenausgang. Die bessere Mischung hat aber auch Einfluss auf
112
Zusammenfassung
die Leistungsparameter der Düse. So steigt bedingt durch die bessere Mischungsrate der
Geschwindigkeitskoeffizient an. Der Durchflusskoeffizient sinkt allerdings aufgrund der
höheren Mischungsverluste ab. Bei den durchgeführten Rechnungen zeigte sich, dass das
SST Modell stärker mischt als das SKE Modell. Dieser Effekt verstärkt sich nochmals bei
Einsatz des Temperaturkorrekturverfahrens mit dem entsprechenden Turbulenzmodell.
Gleichwohl wurden in den Rechnungen beide Turbulenzmodelle mit erster Ordnung dis-
kretisiert. In einer später durchgeführten Parameteruntersuchung wurde jedoch gezeigt,
das sich die Ergebnisse beider Turbulenzmodelle bei zweiter Ordnung annähern.
Ein Vergleich der Simulationen mit experimentell ermittelten Leistungswerten weist eine
gewisse Diskrepanz zwischen beiden Daten auf. So liegen CD und CV praktisch immer
unterhalb der experimentellen Werte. Neben numerischen Gründen wie z. B. Netzab-
hängigkeit, Diskretisierungsschemata oder der Modellierung der Physik können auch
geometrische Ursachen oder Messungenauigkeiten für die Unterschiede verantwortlich
sein. Es zeigt sich jedoch, dass der qualitative Verlauf der Rechnungen mit Temperatur-
korrekturverfahren generell besser wiedergegeben wird.
Ein weiteres Defizit liegt in dem in den RANS/FANS Modellen häufig verwendeten ein-
fachen Schließungsansatz für den turbulenten Wärmestrom. Der Nachteil dieses weit ver-
breiteten Eddy Diffusivity Modells liegt zum einen darin, dass eine turbulente Prandtl-
zahl definiert werden muss. Zum anderen, dass dadurch nur Wärmeströme berücksichtigt
werden, welche von Temperaturgradienten verursacht werden. Deshalb wurde der tur-
bulente Wärmestromvektor durch einen differentiellen Ansatz modelliert und in Open-
FOAM implementiert. Vergleichsrechnungen haben gezeigt, dass das DSM Modell diffe-
renzierter mischt als das EDM. So ist die Mischung bei der verwendeten Düsengeome-
trie im stark turbulenten und dreidimensionalen Strömungsgebiet hinter dem Mischer
leicht stärker als beim EDM mit einer turbulenten Prandtlzahlen von 0,5. Im direkten
Nachlauf des Plugs in Achsennähe, wo die Strömung eher zweidimensionalen Charakter
aufweist, ist die Ausmischung dagegen weniger stark. Die berechneten Leistungsbeiwerte
des DSM für die Heißrechnungen weisen gegenüber dem EDM einen gewissen Versatz
auf. So liegen die Geschwindigkeitskoeffizienten leicht über, die Durchflusskoeffizienten
leicht unterhalb der EDM Ergebnisse. Für TR um 2,5 liegen die erzielten Leistungspa-
rameter des DSM allerdings nahezu auf den Resultaten für das EDM mit P rT = 0, 5. Es
wurde beobachtet, dass die Konvergenz des Differentialgleichungsmodells etwas schlech-
ter ist als das algebraische Verfahren. So wurde für das DSM eine grundsätzlich höhere
Totalenthalpieflussdifferenz zwischen den Einlassflächen von Fan- plus Kernstrom so-
wie Düsenauslassfläche ausgemacht. Diese Diskrepanz könnte aus der Wandgrenzschicht
herrühren, die sich für das differentielle Modell leicht stärker aufheizt, als es beim al-
gebraischen Ansatz der Fall ist. Am größten ist die Totalenthalpieflussdifferenz bei den
eher akademischen Kaltrechnungen. Bei den realitätsnäheren Heißrechnungen scheinen
die Ergebnisse des DSM deutlich verlässlicher.
113
Zusammenfassung
Während das Temperaturkorrekturverfahren durch seine Modellierung auch folgerich-

tig die Impulsgleichungen beeinflusst, ist dies beim differentiellen Wärmestrommodell
nicht der Fall. Allerdings setzt die Temperaturkorrektur auf das Eddy Diffusivity Mo-
dell. Durch die zusätzliche Annahme einer konstanten turbulenten Prandtlzahl kann dies
dazu führen, dass die turbulente Viskosität lokal über- oder unterschätzt wird. Deshalb
scheint eine Verwendung des DSM Verfahrens weiterhin attraktiv, weil zum einen keine
turbulente Prandtlzahl benötigt wird und zum anderen dadurch die Richtung des Wär-
mestromvektors physikalischer wiedergegeben werden kann. Eine Kombination beider
Ansätze mit Temperaturkorrektur und DSM scheint deshalb äußerst sinnvoll, da damit
die genannten Nachteile beider Modelle angegangen werden könnten.
Aus diesem Grund wurde das DSM mit dem Temperaturkorrekturverfahren kombiniert
und ein neues Mischungsmodell ausgearbeitet. Das neue Modell sollte dann anhand
von Rolls-Royce Deutschland bereitgestellten experimentellen Daten validiert werden.
In zahlreichen Kalibrierungs- und Validierungsrechnungen wurde versucht das neue Mo-
dell den Totaltemperaturmessungen von Rolls-Royce anzupassen. Trotz der zum Teil
sehr großen Abweichungen der Einstellungen von den beiden Ursprungsmodellen (DSM
und TC) konnte eine Anpassung an das experimentelle Ergebnis allerdings nicht er-
zielt werden. Während das stark diffusive Totaltemperaturprofil aus den Messungen im
Versuch nicht zufriedenstellend simuliert werden konnte, war die Übereinstimmung mit
experimentell ermittelten Totaldruckprofilen insgesamt deutlich besser. Es ist daher na-
heliegend den Grund für die Diskrepanz der Totaltemperaturprofile nicht ausschließlich
in der Numerik zu suchen, sondern ebenfalls die Temperaturmessungen in Frage zu stel-
len. Zum einen wäre es interessant, ob ein feineres Netz an Messsonden ein weniger
diffusives Totaltemperaturprofil ergeben würde. Eine zweite Hypothese für die Abwei-
chung ist, dass die Messergebnisse aufgrund von Wärmeleitung innerhalb der Traverse
verfälscht worden sein könnten. Der eingeschobene Vergleich mit einer Rechnung auf
einem unstrukturierten Netz bekräftigt diese Vermutung nochmals.
Um die in dieser Arbeit durchgeführten Untersuchungen abzurunden und die Unter-
schiede zwischen Numerik und Experiment besser einordnen zu können wurde eine Pa-
rametervariation durchgeführt. Diese soll Aufschluss auf den Einfluss numerischer und
physikalischer Parameter geben.
So wurde z. B. gezeigt, dass ein feineres Netz zu größeren Leistungsbeiwerten führt. Das
in OpenFOAM zur Leistungsrechnung verwendete Rechennetz weist zwar qualitativ sehr
gute Übereinstimmung mit viel feineren Netzen auf, quantitativ wurden für CD und CV
auf besser aufgelöstem Gitter allerdings nochmals höhere Werte erzielt.
Die für die Leistungsrechnungen verwendete Geometrie des Mischers hat eine andere
Taktung als im Experiment. Daher wurde der Einfluss der Mischertaktung auf das Er-
gebnis untersucht, indem im CFD Modell die ansonsten symmetrische Randbedingung
als periodischer Rand ausgeführt wurde. Obwohl sich dadurch das Strömungsbild än-
114
Zusammenfassung
dert, hat diese Maßnahme quasi keinen Einfluss auf die Leistungskoeffizienten.
Des Weiteren wurde der Einfluss einer variablen spezifischen Wärmekapazität untersucht
und ein von der Temperatur abhängiges cp verwendet. Während sich dadurch am Strö-
mungsbild und der Temperaturverteilung zwar praktisch nichts ändert konnte gezeigt
werden, dass eine variable Wärmekapazität für steigende TR zu höheren CV führt, CD
hingegen sinkt.
Interessantes konnte hinsichtlich dem Verhalten bezüglich der Turbulenzmodelle und
deren räumlicher Diskretisierung beobachtet werden. Bei den in OpenFOAM durchge-
führten Leistungsrechnungen mit Temperaturkorrektur wurde die Turbulenz noch mit
erster Ordnung diskretisiert. Bei erster Ordnung mischen SST und SKE unterschied-
lich stark, d. h. das SST mischt deutlich stärker als das SKE. In der Parameterstudie
wurde allerdings gezeigt, dass sich das SST und das SKE Modell beim Umschalten
auf zweite Ordnung unterschiedlich verhalten. So nimmt die Ausmischung bei zweiter
Ordnung beim SST Modell ab, während sie beim SKE zunimmt. Beide Modelle liefern
dann bei zweiter Ordnung fast identische Ergebnisse. Es ist also anzunehmen, dass auch
die Rechenergebnisse mit Temperaturkorrekturverfahren von SST und SKE bei zweiter
Ordnung näher zusammen liegen würden.
Als Resumee kann zusammenfassend festgehalten werden, dass trotz nicht uneinge-
schränkter Übereinstimmung der Simulationsergebnisse mit experimentellen Messun-
gen, der Einfluss unterschiedlicher numerischer Ansätze zur besseren Berücksichtigung
anisothermer Phänomene im Bezug auf Abgastrahlmischung deutlich wurde. So wurde
zunächst die Auswirkung von Temperaturturbulenz berücksichtigt und deren Wirken
auf das Strömungs- und Mischungsverhalten studiert. Ebenso wurde mit einem differen-
tiellen Ansatz zur Berechnung des turbuleten Wärmestoms eine logische differenziertere
lokale thermische Ausmischung erreicht, als dies beim einfachen turbulenten Prandtlzahl-
ansatz der Fall ist. Allerdings muss festgehalten werden, dass der Rechenaufwand bei
einem differentiellen Modell deutlich größer ist als mit einem algebraischen Ansatz. Mit
der abschließenden Parameteruntersuchung konnten noch numerische und physikalische
sowie geometrische Einflüsse untersucht und eingeordnet werden.
115
Literaturverzeichnis
[1] Calmon, J.: From Sir Frank Whittle to the Year 2000 - Whats is new in Propul-
sion? In: Zeitrschrift für Flugwissenschaften und Weltraumforschung 12 (1988), S.
302–312
[2] Huff, D. L.: Noise Reduction Technologies for Turbofan Engines. In: AIAA Paper
214495 (2007)
[3] Bridges, J. ; Wernet, M. P.: Turbulence Measurements of Separate Flow Nozz-
les With Mixing Enhancement Features. In: AIAA Paper 2484 (2002)
[4] Kenzakowski, D. C. ; Shipman, J. ; Dash, S. M.: Turbulence Model Study of
Laboratory Jets with Mixing Enhancements for Noise Reduction. In: AIAA Paper
16136 (2000)
[5] Tide, P. S. ; Srinivasan, K.: Novel Chevron Nozzle Concepts for Jet Noise
Reduction. In: Proceedings of the Institution of Mechanical Engineers 223 (2009),
S. 51–67
[6] Mengele, V. G.: Relative Clocking of Enhancement Mixing Devices for Jet Noise
Benefit. In: AIAA Paper 0996 (2005)
[7] Wegener, M. ; Lieser, J. ; Mundt, Ch.: Lärm- und leistungsoptimierter Strahl-
mischer. In: DGLR Jahrestagung Berlin DGLR-JT99-82 (1999)
[8] Mundt, Ch. ; Wegner, M.: Improvement of Forced Mixing in Turbofan Engines
with Respect to both Noise and Performance. In: 7. CEAS European Propulsion
Forum, Pau (1999)
[9] Frost, T. H.: Practical Bypass Mixing Systems for Fan Jet Aero Engines. In:
The Aeronautical Quarterly (1966), S. 141–159
[10] Mundt, Ch. ; Lieser, J.: Performance Improvement of Propulsion Systems by
Optimization of the Mixing Efficiency and Pressure Loss of Forced Mixers. In: 8.
CEAS European Propulsion Forum, Nottingham (2001)
[11] Shaw, C. T.: Using Computational Fluid Dynamics. Prentice Hall, 1992
116
[12] Wilcox, David C.: Basic Fluid Dynamics. 1. DCW Industries, 1997
[13] Wilcox, David C.: Turbulence Modeling for CFD. 2. DCW Industries, 1998
[14] Desideri, J. A. ; Chetverushkin, B. N. ; Kuznetsov, Y. A. ; Periaux, J.
; Stoufflet, B.: Experimentation, Modelling and Computation in Flow, Turbu-
lence and Combustion. John Wiley & Sons, 1996
[15] Wessling, Pieter: Principles of Computational Fluid Dynamics. 1. Springer
Verlag, 2001
[16] Koch, L. D. ; Bridges, J.: Flowfield Comparisons From Three Navier-Stokes
Solvers for an Axisymmetric Separate Flow Jet. In: AIAA Paper 0672 (2002)
[17] Engblom, William. A. ; Khavaran, A. ; Bridges, J.: Numerical Prediction of
Chevron Nozzle Noise Reduction using WIND-MGBK Methodology. In: AIAA
Paper 2979 (2004)
[18] Abdol-Hamid, Khaled S. ; Pao, S. P. ; Massey, Steven J. ; Elmiligui, Alaa:
Termperature Corrected Turbulence Model for High Temperature Jet Flow. In:
Journal of Fluids Engineering 126 (2004), S. 844–850
[19] Birch, Stanley F. ; Lyubimov, D. A. ; Secundov, A. N. ; Yakubovsky, K. Y.:
Numerical Modeling Requirements for Coaxial and Chevron Nozzle Flows. In:
AIAA Paper 3287 (2003)
[20] Seiner, J. M. ; Ponton, M. K. ; Jansen, B. J. ; Lagen, T. N.: The Effects of
Temperature on Supersonic Jet Noise Emission. In: AIAA Paper 92-02-046 (1992)
[21] Thomas, R. H. ; Kinzie, K. W. ; Pao, S. P.: Computational Analysis of a Pylon-
Chevron Core Nozzle Interaction. In: AIAA Paper 2001-2185 (2001)
[22] Tam, Christopher K. W. ; Ganesan, Anand: Modified k-ǫ Turbulence Model for
Calculation Hot Jet Mean Flows and Noise. In: AIAA Journal 42 (2004), S. 26–34
[23] Kays, W. M.: Turbulent Prandtl Number - Where are We? In: ASME Journal of
Heat Transfer 116 (1994), S. 284–295
[24] Daly, B. J. ; Harlow, F. H.: Transport Equations in Turbulence. In: Physics of
Fluids 13 (1970), S. 2634–2649
[25] Abe, K. ; Kondoh, T. ; Nagano, Y.: A Two-Equation Heat Transfer Model
Reflecting Second-Moment Closures for Wall and Free Turbulend Flows. In: Int.
Journal of Heat Fluid Flow 17 (1996), S. 228–237
[26] So, R. M. C. ; Sommer, T. P.: An Expicit Algebraic Heat-Flux Model for the
Temperature Field. In: Int. J. Heat Fluid Flow 7 (1995), S. 455–465
117
[27] Ronki, M. ; Gatski, T. B.: Predicting Turbulent Convective Heat Transfer in

Fully Developed Duct Flows. In: Int. J. Heat Fluid Flow 22 (1981), S. 381–392
[28] Chidambaram, N. ; Dash, S. M. ; Kenzakowski, D. C.: Scalar Variance Trans-
port in the Trubulence Modeling of Propulsive Jets. In: Journal of Propulsion and
Power 17 (2001), S. 79–84
[29] Kenzakowski, D. C. ; Papp, J. ; Dash, S. M.: Evaluation of Advanced Turbu-
lence Models and Variable Prandl/Schmidt Number Methodology for Propulsive
Flows. In: AIAA Paper 16720 (2000)
[30] Bremhorst, K. ; Bullock, K. J.: Spectral Measurementes of Temperature and
Longitudinal Velocity Fluctuations in Fully Developed Oipe Flow. In: Internationl
Journal of Heat and Mass Transfer 13 (1970), S. 1313–1329
[31] Bremhorst, K. ; Bullock, K. J.: Spectral Measuremente of Turbulent Heat
and Momentum Transfer in Fully Developed Pipe Flow. In: Internationl Journal
of Heat and Mass Transfer 16 (1973), S. 2141–2154
[32] Lösch, S.: Influence of Flow Mixing on the Performance of a Turbofan Engine of
High By-Pass Ratio. In: Zeitschrift für Flugwissenschaften und Weltraumforschung
14 (1990), S. 333–341
[33] Mao, R. ; Yu, S. C. M. ; Zhou, T. ; Chua, L. P.: On the Vorticity Characteristics
of lobe-Forced Mixer at Different Configurations. In: Experiments in Fluids 46
(2009), S. 1049–1066
[34] Sventitskiy, A. ; Mundt, Ch.: Simulation of the Infrared Emission from Nons-
cattering Aircraft Engine Exhaust Plumes Using Statistical Narrowband Models
in Conjunction with Model Parameters Based on the Modern Spectrosopic Data.
In: AIAA Paper 2012-3315 (2012)
[35] Babbitt, R.R. ; Cohn, J. A. ; Fleming, K. J.: Advanced High Bypass Mixed-
Flow Exhaust System Design Study. In: AIAA Paper AIAA-91-2242 (1991)
[36] Pearson, H.: Mixing of Exhaust and By-Pass Flow in a By-Pass Engine. In:
Journal of the Royal Aeronautical Society 66 (1962)
[37] Kerrebrock, J. L.: Aircraft Engines and Gas Turbines. 2. The MIT Press, 1996
[38] Oates, G. C.: Aerothermodynamics of Gas Tubine and Rocket Propulsion - Revi-
sed and Enlarged. 3. AIAA Education Series, 1988
[39] Münzberg, Hans-Georg: Flugantriebe. Springer-Verlag, 1972
118
[40] Staudacher, S.: Luftfahrtantriebe - Vorlesungsskript. 2005. – Universität Stutt-

gart
[41] Bräunling, W. J. G.: Flugzeugtriebwerke. Springer-Verlag, 2001
[42] Mengele, V. G. ; Dalton, W. N.: Lobed mixer Design for Noise Suppression:
Acoustic and Aerodynamic Test Data Analysis. In: NASA/CR-2002-210823 1
(2002)
[43] Mengele, V. G. ; Baker, V. D. ; Dalton, W. N.: Lobed mixer Design for Noise
Suppression: Plume, Aerodynamic and Acoustic Data. In: NASA/CR-2002-210823
2 (2002)
[44] Lighthill, M. J.: On Sound Generated Aerodynamically, General Theory (I) and
Turbulence as a Source of Sound (II). In: Proc. Roy. Soc A221 (1954), S. 564–587
[45] Mengele, V. G.: Anomalous Effefct of Nozzle Length Reduction on Jet Noise of
Forced Mixers. In: AIAA Paper 99-1968 (1999)
[46] Mengele, V. G.: Jet Noise Reduction by Lobed Mixers with Boomerang Scallops.
In: AIAA Paper 99-1923 (1999)
[47] Yu, S. C. M.: Some Flow Characteristics of Lobed Forced Mixers at Higher
Velocity Ratios. In: Experiments in Fluids 20 (1996), S. 234–237
[48] Banzhaf, B.: Untersuchungen zum Betriebsverhalten gemischter Abgassysteme in
Luftstrahltriebwerken im niedrigen Lastbereich. 2009. – Dissertation, Universtität
Stuttgart
[49] Deinert, B.: Mischungsverstärker in konstant gekreuzten Scherschichten und im
Abgasblütenmischer. 2006. – Dissertation, Technische Universität Berlin
[50] Klötzer, A. ; Staudacher, S. ; Weißschuh, M.: Investigations on the Influ-
ence of Altitude Conditions on Mixed Exhaust Systems. In: ASME Turbo Expo
GT-2011-46175 (2011)
[51] Banzhaf, B. ; Staudacher, S.: Untersuchungen zum Betriebsverhalten des
Triebwerksmischers im niedrigen Lastbereich. In: DGLR-2006-095 (2006)
[52] Klötzer, A.: Untersuchungen zum Betriebsverhalten gemischter Abgassysteme in
Luftstrahltriebwerken unter Höhenbedingungen. 2012. – Dissertation, Universität
Stuttgart
[53] Mattingly, J. D. ; Heiser, W. H. ; Pratt, D. T.: Aircraft Engine Design. 2.
American Institute of Aeronautics and Astronautics, 2002
[54] Windemuth, Eberhard: Strömungstechnik. Springer-Verlag, 1984
119
[55] ANSYS FLUENT 14.0 Theory Guide 2011

[56] ANSYS FLUENT 14.0 User’s Guide 2011
[57] Borm, O. ; Jemcov, A. ; Kau, H.-P.: Density Based Navier-Stokes Solver for
Transonic Flows. 2011. – 6th OpenFOAM Workshop, PennState University, USA
[58] Borm, O.: Instationäre numerische Untersuchung der aerodynamischen Rotor-
Stator-Interaktion an einem Radialverdichter. 2012. – Dissertation, TU München
[59] Frohn, A.: Einführung in die Technische Thermodynamik. 3. Konrad Wittwer,
Stuttgart, 1998
[60] Hirsch, C.: Numerical Computation of Internal and External Flows. Volume 1:
Funamentals of Numerical Discretization. John Wiley & Sons, 1988
[61] Blazek, J.: Computational Fluid Dynamics: Principles and Applications. 2.
Elsevier, 2005
[62] Oertel, H. ; Böhle, M. ; Reviol, T: Strömungsmechanik. 6. Vieweg + Teubner,
2011
[63] Sutherland, William: The Viscosity of Gases and Molecular Force. In: Philo-
sophical Magazine 36 (1893), S. 507–531
[64] Spalart, P. ; Allmaras, S.: A One-Equation Turbulence Model for Aerodynamic
Flows. In: AIAA Paper 92-0439 (1992)
[65] Gibson, M. M. ; Launder, B. E.: Ground Effects on Pressure Fluctuations in
the Atmospheric Boundary Layer. In: Journal of Fluid Mechanics 86 (1978), S.
491–511
[66] Launder, B. E.: Second-Moment Closure: Present... and Future? In: Int. Journal
of Heat and Fluid Flow 10(4) (1989), S. 282–300
[67] Launder, B. E. ; Reece, G. J. ; Rodi, W.: Progress in the Development of a
Reynolds-Stress Turbulence Fluxes. In: Journal of Fluid Mechanics 68(3) (1975),
S. 537–566
[68] Jones, W. ; Launder, B. E.: The Prediction of Laminarization with a Two-
Equation Model of Turbulence. In: Int. Journal of Heat and Mass Transfer 15
(1972), S. 301–314
[69] Jones, W. ; Launder, B. E.: The Prediction of Low-Reynolds-Number Pheno-
mena with a Two-Equation Model of Turbulence. In: Int. Journal of Heat and
Mass Transfer 16 (1973), S. 1119–1130
120
[70] Launder, B. E. ; Spalding, D. B.: Lectures in Mathematical Models of Turbu-

lence. Academic Press, London, England, 1972
[71] Versteeg, H. K.: Computational Fluid Dynamics: The Finite Volume Method.
2. Pearson, 2007
[72] Shih, Tsan-Hsing ; Liou, William W. ; Shabbir, Aamir ; Yang, Zhiang ; Zhu,
Jiang: A New k-ǫ Eddy Viscosity Model for High Reynolds Number Turbulent
Flows. In: Computers Fluids 24 (1995), S. 227–238
[73] Yakhot, V. ; Orszag, S. A.: Renormalization Group Analysis of Turbulence I
Basic Theory. In: Journal of Scientific Computing 1(1) (1986), S. 1–51
[74] Sarkar, S. ; Lakshmanan, B.: Application of a Reynolds Stress Turbulence
Model to the Compressible Shear Layer. In: AIAA Journal 16 (1991), S. 743–749
[75] Menter, F. R.: Zonal Two Equation k-omega Turbulence Models for Aerodyna-
mic Flows. In: AIAA Paper 93-2906 (1993)
[76] Menter, F. R.: Two-Equation Eddy-Viscosity Turbulence Models for Engineering
Applications. In: AIAA Journal 8 (1994), S. 269–289
[77] Menter, F. R. ; Kuntz, M. ; Langtry, R.: Ten Years of Industrial Experience
with the SST Turbulence Model. In: Heat and Mass Transfer 4 (2003)
[78] Wilcox, D. C.: Reassessment of the Scale-Determining Equation for Advanced
Turbulence Models. In: AIAA Journal 26 (1988), S. 1299–1310
[79] Suzen, Y. B. ; Hoffmann, K. A.: Investigation of supersonic Jet Exhaust Flow
by One- and Two-Equation Turbulence Models. In: AIAA Paper 98-0322 (1998)
[80] Weiss, Jonathan M. ; Smith, Wayne A.: Preconditioning Applied to Variable
and Constant Density Flows. In: AIAA Journal 33 (1995), S. 2050–2057
[81] Einfeldt, B.: On Godunov-Type Methods for Gas Dynamics. In: SIAM J.
Numer. Anal. 25(2) (1988), S. 294–318
[82] Batten, P. ; Leschziener, M. A. ; Goldberg, U. C.: Average-State Jacobians
and Implicit Methods for Compressible Viscous and Turbulent Flows. In: Journal
of Computational Physics 137 (1997), S. 38–78
[83] Toro, E. F. ; Spruce, M. ; Spears, W.: Restoration of the Contact Surface in
the HLL-Riemann Solver. In: Shock Waves 4 (1992), S. 25–34
[84] Liou, M. S.: On a New Class of Flux Splittings. In: Proc. 13th Int. Conf. on
Numerical Methods in Fluid Dynamics, Rome, Italy (1992)
121
[85] Liou, M. S. ; Steffen, C. J. J.: A New Flux Splitting Scheme. In: Journal of
Computational Physics 107 (1993), S. 23–39
[86] Liou, M. S. ; Steffen, C. J. J.: A Sequal to AUSM: AUSM+. In: Journal of
Computational Physics 129 (1996), S. 364–382
[87] Liou, M. S.: A Sequal to AUSM, Part II: AUSM+up for all Speeds. In: Journal
of Computational Physics 214 (2006), S. 137–170
[88] Roe, P. L.: Approximate Riemann Solvers, Parametric Vectors, and Difference
Schemes. In: Journal of Computational Physics 43 (1981), S. 357–372
[89] Roe, P. L. ; Pike, J.: Efficient Construction and Utilisation of Approximate
Riemann Solutions. North Holland Publishing, 1984
[90] Harten, A. ; Lax, P. D. ; Van Leer, B.: On Upstream Differencing and
Godunov-Type Schemes for Hyperbolic Conservation Laws. In: Soc. Industrial
and Applied Mathematics Rew. 25(1) (1983)
[91] Harten, A. ; Hayman, J. M.: Self Adjusting Grid Methods for One-Dimensional
Hyperbolic Conservation Laws. In: Journal of Computational Physics 50 (1983),
S. 235–269
[92] Barht, Timothy J. ; Jesperen, Dennis C.: The Design and Application of
Upwind Schemes on Unstructured Meshes. In: AIAA Paper AIAA-89-0366 (1989)
[93] Venkatakrishnan, V.: On the Accuracy of Limiters and Convergence to Steady
State Solutions. In: AIAA Paper 93-0880 (1993)
[94] Saegeler, S. ; Mundt, Ch.: Implementation of a Preconditioner for a Density-
Based Navier-Stokes Solver. In: 7th OpenFOAM Workshop, Technische Universität
Darmstadt (2012)
[95] Hirsch, C.: Numerical Computation of Internal and External Flows. Volume 2:
Computational Methods for Inviscid and Viscous Flows. John Wiley & Sons, 1990
[96] Toro, E. F.: Riemann Solvers and Numerical Methods for Fluid Dynamics - A
Practical Introduction. 3. Springer-Verlag, 2009
[97] Arnone, Andrea ; Liou, Meng-Sing ; Povinelli, Louis A.: Multigrid Time-
Accurate Integration of Navier-Stokes Equations. In: NASA Forschungsbereicht
NASA TM 106373 (1993)
[98] Mavriplis, D. J. ; Jameson, A. ; Martinelli, L.: Multigrid Solution of the
Navier-Stokes Equations on Triangular Meshes. In: AIAA Paper 89-0120 (1989)
122
[99] Möller, Christine: Bestimmung von Mischer- und Düsenkoffizenten mittels CFD-
Methoden und Bewertung des Einflusses auf die Modellierung des Betriebsverhal-
tens eines Triebwerkes mit Strahlmischung. 2006. – Diplomarbeit, Universität
Stuttgart
[100] Ferziger, J. H. ; Perić, M.: Computational Mehtods for Fluid Dynamics. 3.
Springer, 2002
[101] Chien, K.-Y.: Predictions of Channel and Boundary Layer Flows with a Low-
Reynolds-Number Turbulence Model. In: AIAA Journal 20 (1982), S. 33–38
[102] Thies, Andrew. T. ; Tam, Christopher K. W.: Computation of Turbulent Axisym-
metric and Nonaxisymmetric Jet Flows Using the k-ǫ Model. In: AIAA Journal
34 (1996), S. 309–316
[103] Pope, S. B.: An Explanation of the Round-Jet/Plane-Jet Anomaly. In: AIAA
Journal 16 (1978), S. 279–281
[104] Engblom, William A. ; Georgiadis, Nicholas J. ; Khavaran, Abbas: Inves-
tigation of Variable-Diffusion Turbulence Model Correction for Round Jets. In:
[105] Georgiadis, Nicholas J. ; Yoder, Dennis A. ; Engblom, William A.: Evaluation
of Modified Two-Equation Turbulence Models for Jet Flow Predictions. In: AIAA
Paper 490 (2006)
[106] Georgiadis, Nicholas J. ; DeBonis, James R.: Navier-Stokes Analysis Methods
for Tubulent Jet Flows with Applications to Aircraft Exhaust Nozzles. In: Progress
in Aerospace Sinces 42 (2006)
[107] Beguier, C. ; Dekeyser, I. ; Launder, B. E.: Ratio of Scalar and Velocity
Dissipation Time Scales in Shear Flow Turbulence. In: Phys. Fluids 21(3) (1978),
S. 307–310
[108] Sadiki, Amsini: Thermodynamik und Turbulenzmodellierung - Habilitations-
schrift. 1998. – Technische Universität Darmstadt
[109] Massey, Steven J. ; Thomas, Russell H. ; Abdol-Hamid, Khaled S. ; Elmili-
gui, Alaa A.: Computational and Experimental Flowfield Analyses of Seperate
Flow Chevron Nozzles and Pylon Interaction. In: AIAA Paper 2003-3212 (2003)
[110] Saegeler, S. ; Lieser, J. ; Mundt, Ch.: Improved Modelling of Vortical Mixing
for the Simulation of Efficient Propulsion Systems. In: 28th International Congress
of the Aeronautical Sciences (2012)
123
[111] DeGroot, Sybren R.: Thermodynamik Irreversibler Prozesse. B.I. Mannheim -

Hochschultaschenbücher, 1960
[112] DeGroot, Sybren R. ; Mazur, Peter: Grundlagen der Thermodynamik Irrever-
sibler Prozesse. B.I. Mannheim - Hochschultaschenbücher, 1969
[113] Saegeler, S. ; Mundt, Ch.: Einfluss der Modellierung von Temperaturturbulenz
auf die Mischung in Triebwerksstrahlen. In: 61. Deutscher Luft- und Raumfahrt-
kongress (2012)
[114] Dembowski, Mary A. ; Georgiadis, Nicholas J.: An Evaluation of Parameters
Influencing Jet Mixing Using the WIND Navier-Stokes Code. In: AIAA Paper
221727 (2002)
[115] Nagano, Y. ; Kim, C.: A Two-Equation Model for Heat Transport in Wall Tur-
bulent Shear Flows. In: Journal of Heat Transfer 110 (1988), S. 583–589
[116] Youssef, M. S. ; Nagano, Y. ; Tagawa, M.: A Two-Equation Heat Transfer
Model for Predicting Turbulent Thermal Fields under Arbitrary Wall Thermal
Conditions. In: Int. Journal of Heat and Mass Transfer 35 (1992), S. 3095–3104
[117] Deng, B. ; Wu, W. ; Xi, S.: A Near-Wall Two-Equation Heat Transfer Model for
Wall Turbulent Flows. In: Int. Journal of Heat and Mass transfer 24 (2000), S.
691–698
[118] Kenzakowski, Donald C.: RANS Modeling Improvements for Jets Using Scalar
Variance Equations. In: AIAA Paper 491 (2006)
[119] Brinckman, Kevin W. ; Kenzakowski, Donald C. ; Dash, Sanford M.: Progress
in Practical Scalar Fluctuation Modeling for High-Speed Aeropropulsive Flows. In:
[120] Brinckman, Kevin W. ; William H. Calhoon, Jr. ; Mattick, Stephen J. ;
Tomes, Jeremy ; Dash, Sanford M.: Scalar Variance Model Validation for High-
Speed Variable Composition Flows. In: AIAA Paper 715 (2006)
[121] Launder, B. E.: On the Computation of Convective Heat Transfer in Complex
Turbulent Flows. In: Journal of Heat Transfer 110 (1988), S. 1112–1128
[122] Launder, B. E.: Topics in Applied Physics - Turbulence. Springer-Verlag, 1976
[123] Lai, Y. G. ; So, R. M. C.: Near-Wall Modeling of Turbulent Heat Closures. In:
Int. J. Heat and Mass Transfer 33(7) (1990), S. 1429–1440
[124] Johnson, R. W.: Handbook of Fluid Dynamics. 1. Springer-Verlag, 1998
124
[125] Launder, B. E.: Scalar Property Transport by Turbulence. In: Rep. No.
HTS/73/26, Dep. Mech. Engng, Imperical College, London (1973)
[126] Launder, B. E. ; Samaraweera, D. S. A.: Application of a Second-Moment
Turbulence Closure to Heat and Mass Transport in Thin Shear Flows - I. Two-
Dimensional Transport. In: Int. J. Heat Mass Transfer 22 (1979), S. 1631–1643
[127] Donaldson, C. duP. ; Sullivan, R. D. ; Rosenbaum, H.: A Theoretical Study
of the Generation of Atmospheric-Clear Air Turbulence. In: AIAA Journal 10,2
(1972), S. 162–170
[128] Wyngaard, J. C. ; Coté, O. R.: The Evolution of a Convective Planetary Boun-
dary Layer: A Higher-Order-Closure Model Study. In: Boundary Layer Meteorol.
7 (1974), S. 289–307
[129] Saegeler, S. ; Mundt, Ch.: Advanced Numerical Simulation of Mixing Hot
Core and Cold Bypass Flow in Modern Propulsion Systems with Internal Lobed
Forced Mixer. In: 21st AIAA Computational Fluid Dynamics Conference, San
Diego, USA (2013)
[130] Kim, Sung-Eun ; Markov, Boris ; Caraeni, Doru: A Mulit-Dimensional Linear
Reconstruction Scheme for Arbitrary Unstructured Grids. In: AIAA Paper 3990
(2003)
125
A. Kompressibel-turbulenter
skalarer Transport
Hier soll die Herleitung der Korrelation u′′i c′′ für den kompressiblen Fall durchgeführt wer-
den. u′′i ist wie üblich die turbulente Geschwindigkeitsgröße. c′′ repräsentiert die Schwan-
kungsgröße eines beliebigen Skalars: c = C̃ + c′′ .
ℵ(ui ) sei der Navier-Stokes Operator und Σ(c) der Transportgleichungs-Operator für
einen beliebigen Skalar:
∂ui ∂ui ∂p ∂ 2 ui
ℵ(ui ) = ρ + ρuk + − µ 2 = 0, (A.1)
∂t ∂xk ∂xi ∂xk
∂c ∂c ∂ 2c
Σ(c) = ρ + ρuk −γ = 0, (A.2)
∂t ∂xk ∂xk 2
mit γ als Diffusitität des Skalars c.

Um die Korrelation für u′′i c′′ zu erhalten wird ℵ(ui ) mit der fluktuierenden Größe c′′ und
Σ(c) mit der fluktuierenden Geschwindigkeit u′′i multipliziert, anschließend addiert und
dann gemittelt:
c′′ ℵ(ui ) + u′′i Σ(c) = 0. (A.3)
Um die Vorgehensweise nachvollziehbarer zu machen, sollen im Folgenden die Terme aus

(A.3) einzeln betrachtet werden. Begonnen wird mit dem instationären Term1 :
1
Die Indizes , t und , k bedeuten eine Differenzierung nach der Zeit bzw. nach dem Ort.
126
Anhang A
∂ui ∂c
ρc′′ + ρu′′i = ρc′′ (Ũi + u′′i ),t + ρu′′i (C̃ + c′′ ),t
∂t ∂t
= ρc′′ u′′i,t + ρu′′i c′′,t
∂ρ (A.4)
= (ρc′′ u′′i ),t − u′′i c′′
∂t
∂(ρui c )
′′ ′′
∂ρ
= − u′′i c′′ .
∂t ∂t
Der konvektive Anteil ergibt folgende Lösung:
∂ui ∂c
ρc′′ + ρu′′i =ρc′′ (Ũk + u′′k )(Ũi + u′′i ),k + ρu′′i (Ũk + u′′k )(C̃ + c′′ ),k
∂xk ∂xk
=ρc′′ (Ũk Ũi,k + Ũk u′′i,k + u′′k Ũi,k + u′′k u′′i,k )
+ ρu′′i (Ũk C̃,k + Ũk c′′,k + u′′k C̃,k + u′′k c′′,k )
=ρc′′ u′′k Ũi,k + ρu′′i u′′k C̃,k + ρc′′ u′′i,k uk + ρc′′,k u′′i uk
=ρc′′ u′′k Ũi,k + ρu′′i u′′k C̃,k + ρ(c′′ u′′i ),k uk
=ρc′′ u′′k Ũi,k + ρu′′i u′′k C̃,k + (ρc′′ u′′i uk ),k − c′′ u′′i (ρuk ),k
=ρc′′ u′′k Ũi,k + ρu′′i u′′k C̃,k + (ρc′′ u′′i u′′k ),k + (ρc′′ u′′i Ũk ),k − c′′ u′′i (ρuk ),k
∂ Ũi ∂ C̃ ∂
=ρc′′ u′′k + ρu′′i u′′k + (ρc′′ u′′i u′′k )
∂xk ∂xk ∂xk
∂
+ (ρc′′ u′′i Ũk ) − c′′ u′′i (ρuk ),k .
∂xk
(A.5)
Der Druckgradiententerm ergibt
∂p
c′′ = c′′ (P + p′ ),i
∂xi
∂P
= c′′ + (c′′ p′ ),i − c′′,i p′
∂xi
(A.6)
∂P ∂ ′′ ′ ∂c′′
= c′′ + (c p ) − p′
∂xi ∂xi ∂xi
∂P ∂p ′
= c′′ + c′′ .
∂xi ∂xi
127
Anhang A
Letztlich fehlt noch die Betrachtung des viskosen Anteils (unter der Annahme konstanter
bzw. geringer Veränderlichkeit der Stoffeigenschaften):
∂ 2 ui 2
′′ ∂ c
−µc′′ − γu i = −µc′′ ui,kk − γu′′i c,kk
∂x2k ∂x2k
= −µc′′ (Ũi + u′′i )i,kk − γu′′i (C̃ + c′′ ),kk
= −µc′′ u′′i,kk − γu′′i c′′,kk (A.7)

= − (µc′′ u′′i,k ),k − µc′′,k u′′i,k − (γu′′i c′′,k ),k − γu′′i,k c′′,k
" #
∂ ∂u ′′
i ∂c ′′ ∂u′′ ∂c′′
=− µc′′ + γu′′i + (µ + γ) i .
∂xk ∂xk ∂xk ∂xk ∂xk
Durch Substitution der Gleichungen (A.4)-(A.7) in (A.3) erhält man somit die Favre
gemittelte Gleichung für den turbulenten Transport einer turbulenten skalaren Größe
c′′ . Um letztlich zur vollständigen Gleichung zu gelangen kann man von der Tatsache
Gebrauch machen, dass die Summe der letzten Terme auf der rechten Seite der Gleichun-
gen (A.4) und (A.5) verschwindet, da ihre Summe proportional zu den beiden Termen
in der Kontinuitätsgleichung ist. Die finale Gleichung lautet dann
∂ ∂ ′′ ′′ ∂ Ũi ∂ C̃ ∂
ρu′′i c′′ + ρc ui Ũk = − ρc′′ u′′k − ρu′′i u′′k − ρc′′ u′′i u′′k
∂t ∂xk ∂xk ∂xk ∂xk
∂P ∂p ′
− c′′ − c′′ (A.8)
∂xi ∂xi
" #
∂ ∂u′′ ∂c′′ ∂u′′ ∂c′′
+ µc′′ i + γu′′i − (µ + γ) i .
Für den turbulenten Transport von Temperaturturbulenz θ analog zu Gleichung (A.8)

erhält man dann folglich:
∂ ∂ ∂ Ũi ∂ T̃ ∂
ρu′′i θ + ρθu′′i Ũk = − ρθu′′k − ρu′′i u′′k − ρθu′′i u′′k
∂t ∂xk ∂xk ∂xk ∂xk
∂P ∂p ′
−θ −θ (A.9)
∂xi ∂xi
" #
∂ ∂u′′i ∂θ ∂u′′ ∂θ
+ µθ + λu′′i − (µ + ρα) i .
128
Anhang A
Beim Vergleich von Gleichung (A.9) mit der inkompressiblen Gleichung (4.11) aus [123]
fällt auf, dass im Prinzip lediglich der durch den Druckgradient zusätzlich entstande-
ne Term −θ∂P/∂xi in beiden Gleichungen abweicht. Alle anderen Terme sind quasi
identisch. Eine Schließung dieses zusätzlichen durch die Korrelation mit dem Druck-
gradienten hervorgerufenen Terms ist nicht möglich, da hierfür kein Schließungsansatz
gefunden wurde.
Ähnliches lässt sich auch bei einem Vergleich der kompressiblen mit der inkompressiblen
Transportgleichung für die turbulente kinetische Energie beobachten [13]. Sowohl die
kompresiblen als auch die inkompressiblen Gleichungen sind praktisch identisch bis auf
Terme die aus einer Korrelation der Geschwindigkeit mit dem Druck resultieren. In [13]
werden hierfür zwar spärlich Schließungsterme mit mehr oder weniger guter Validierung-
grundlage vorgeschlagen, in der Praxis finden diese aber praktisch kaum Anwendung und
werden in der Regel vernachlässigt. Es wird daher aufgrund der Struktur von (A.9) im
Vergleich zur inkompressiblen Gleichung und wegen der Analogie der Problematik bei
den gängigen Turbulenzmodellen angenommen, dass Gleichung (4.11) in konservativer
Form auch voll für den kompressiblen Fall anwendbar ist und diese Druckkorrelation
vernachlässigt wird
129
B. Modifizierte ω-Gleichung für
Temperaturkorrekturmodell
In diesem Kapitel wird die Begründung für den zusätzlichen Temperaturkorrekturfaktor
CT in Gleichung (4.8) gegeben. Dieser ist nötig um die Temperaturkorrektur, welche für
das k-ǫ Turbulenzmodell entwickelt wurde, konsistent mit dem k-ω-SST Modell verwen-
den zu können.
Ausgangspunkt ist die exakte Transformation der ǫ-Gleichung in die ω-Form mittels der
bereits bekannten Beziehung ǫ = β ∗ kω (vgl. Gleichung (4.5)). Auf den Rechenweg der
Transformation soll hier verzichtet werden. Stattdessen wird die in [13] bereits transfor-
mierte, inkompressible Form übernommen:
ej
∂ω e ∂ω ω ∂U 2 ∂ ∂ω
+ Ui =α τij − βω + (ν + σνT )
∂t ∂xi k ∂xi ∂xj ∂xj
(B.1)
(ν + σνT ) ∂k ∂ω ω ∂ ∗ ∂k
+2 + (σ − σ ) νT .
k ∂xi ∂xi k ∂xj ∂xj
α, β, σ und σ ∗ sind einfache Modellkonstanten. Legt man den Fokus auf die Berechnung
von Freistrahlen, kann die molekulare Viskosität ν vernachlässigt werden (bzw. hier soll
sie nur im Cross-Diffusions Term nicht berücksichtigt werden). Geht man weiter davon
aus, dass die Modellkonstanten mit σ = σ ∗ gleich sind, vereinfacht sich obige Gleichung
zu
ej
∂ω e ∂ω ω ∂U ∂ ∂ω νT ∂k ∂ω
+ Ui = α τij 2
− βω + (ν + σνT ) + 2σ . (B.2)
∂t ∂xi k ∂xi ∂xj ∂xj k ∂xi ∂xi
Hält man der nun erhaltenen Formulierung die inkompressible Gleichung von Menter
entgegen (in Anlehnung an Gleichung (3.38), ferner wird nur ihr Anteil fernab fester
Wände betrachtet),
130
Anhang B
ej
∂ω e ∂ω γ ∂U ∂ ∂ω 1 ∂k ∂ω
+ Ui = τij 2
− βω + (ν + σω2 νT ) + 2σω2 , (B.3)
∂t ∂xi νT ∂xi ∂xj ∂xj ω ∂xi ∂xi
fallen die Unterschiede beider Formulierungen auf. Sieht man von den modifizierten Kon-
stanten α und γ, sowie σ und σω2 in beiden Gleichungen ab, wird jedoch die unterschied-
liche Formulierung des ersten und letzten Terms auf der rechten Seite augenscheinlich.
Im Produktionsterm von ω wurde ω/k von Menter offensichtlich durch 1/νT ersetzt.
Ebenfalls wurde der Cross-Diffusions-Term modifiziert. Hier wurde νT /k durch 1/ω aus-
getauscht. Diese Substitution durch Menter ist gültig, sofern man für die Berechnung
der turbulenten Viskosität folgende Formel annimmt
k
νT = . (B.4)
ω
Sie ist ebenfalls gültig für Menters neue turbulente Viskosität (siehe Gleichung (3.39) für
das SST Modell). Diese geht fernab von festen Wänden wieder in Gleichung (B.4) über.
Modifiziert man die in Gleichung (B.4) gegebene Beziehung, z.B. durch einen beliebigen
Faktor C, muss dies für eine korrekte Substitution folglich auch in der Transportgleichung
für ω berücksichtigt werden. Statt Gleichung (B.4) gilt jetzt nämlich
k
νT = C · . (B.5)
ω
Mit Gleichung (B.2), in Verbindung mit der Beziehung (B.5), lautet die neue Transport-
gleichung für ω nun
ej
∂ω e ∂ω C ∂U ∂ ∂ω C ∂k ∂ω
+ Ui = γ τij 2
− βω + (ν + σνT ) + 2σ . (B.6)
∂t ∂xi νT ∂xi ∂xj ∂xj ω ∂xi ∂xi
131

Verbesserung Der Numerischen Simulation Der Mischung Von Triebwerksstrahlen

Hochgeladen von

Copyright:

Verfügbare Formate

Verbesserung Der Numerischen Simulation Der Mischung Von Triebwerksstrahlen

Hochgeladen von

Dokumentinformationen

Originaltitel

Copyright

Verfügbare Formate

Dieses Dokument teilen

Dokument teilen oder einbetten

Freigabeoptionen

Stufen Sie dieses Dokument als nützlich ein?

Sind diese Inhalte unangemessen?

Copyright:

Verfügbare Formate

Verbesserung Der Numerischen Simulation Der Mischung Von Triebwerksstrahlen

Hochgeladen von

Copyright:

Verfügbare Formate

Universität der Bundeswerhr München

Fakultät für Luft- und Raumfahrttechnik

Verbesserung der numerischen Simulation

Vorsitzender: Prof. Dr.-Ing. R. Förstner

2. Prinzip und Leistungsbewertung von Schubdüsen mit internem Mischer 9

4. Simulation der Düsenströmung mit OpenFOAM 42

5. Parameteruntersuchung der Düsenströmung mit ANSYS FLUENT 98

A. Kompressibel-turbulenter skalarer Transport 126

B. Modifizierte ω-Gleichung für Temperaturkorrekturmodell 130

ℵ(ui ) Navier-Stokes Operator

Alle physikalische Größen sind in SI-Einheiten angegeben

1stO Genauigkeit erster Ordnung im Raum

Abbildung 1.1.: Relativer spezifischer Treibstoffverbrauch unterschiedlicher Triebwerks-

Abbildung 1.2.: Triebwerk mit Chevrondüsen.

Bei Zweistromtriebwerken mit mittlerem Bypassverhältniss werden heutzutage haupt-

1.2. Stand der Forschung

1.3. Zielsetzung der Arbeit

2.1. Grundprinzip interner Strahlmischung

2.2. Quantitative Leistungsbewertung

Der Durchflusskoeffizient CD definiert sich aus dem Verhältnis von verlustbehafteten

Das ideale Q hingegen ist definiert als

Die Fallunterscheidung für das kritische Druckverhältnis wird vorgenommen um den

Der Geschwindigkeitskoeffizient CV stellt per Definition das Verhältnis von tatsächli-

3.1. Navier-Stokes Gleichungen

Der oben aufgeführte Satz an Erhaltungsgleichungen (3.1)-(3.3) lässt sich in kompakter

~ beinhaltet den Satz der sogenannten konservativen Variablen

F~ ist der Vektor konvektiver Flüsse:

~ Auf den Inhalt dieses

3.2. Gemittelte Navier-Stokes Gleichungen

Der mittlere statische Druck berechnet sich analog zu (3.6):

ei ist der gemittelte Geschwindigkeitsanteil, u′′ i ist ihr schwankender Anteil1 .

und die Totalenergie zu

3.3.1. Allgemeine Schließungsansätze

Viskoser und turbulenter Spannungstensor

Der Verlauf von µ(T ) ist in Abbildung 3.1 grafisch dargestellt.

Molekularer und turbulenter Wärmestrom

Abbildung 3.1.: Polynom µ(T ).

Mit Einführung einer laminaren Prandtlzahl P r kann die Temperaturleitfähigkeit λ ab-

Abbildung 3.2.: Polynom λ(T ).

wobei α und αT wie folgt berechnet werden:

Molekulare Diffusion und turbulenter Transport

Der geschlossene Satz gemittelter Navier-Stokes Gleichungen in differenzieller Form er-

auf die entsprechende Literatur [77] verwiesen.

Die turbulente Machzahl M aT berechnet sich wieder gemäß Gleichung (3.35).

3.4. Strömungsberechnung und Lösungsverfahren

3.4.1. Numerische Flussberechnung

Strömungsfeldern eine schlechte Konditionierung des Systems. Die hohen Konditionszah-

erhöhen die Steifigkeit des Gleichungssystems, was den Konvergenzprozess verlangsamt

Der Precondition-Parameter Θ wird angegeben mit

Die Referenzgeschwindigkeit Ur wird durch folgende Beziehung berechnet:

Abbildung 3.3.: Darstellung der Ränder des Rechengebiets.

Der Geschwindigkeitsvektor wird unter Vorgabe der Einlassrichtungsvektors und der

Cµ hat den bereits bekannten Wert von 0,09.

Am Auslass des Fernfelds werden die Randbedingungen für Temperatur, Geschwindig-

3.6.1. Lokales Zeitschrittverfahren

Der viskose Zeitschritt ist gegeben durch

∆t = CF L · min [∆tinv , ∆tvis ] . (3.60)

3.6.2. Runge-Kutta Zeitschrittverfahren

~ (i) ist die Lösung zum vorherigen Zeitschritt, Q