Untitled
Untitled
Untitled
von
Stefan Lipp
Dissertation, Karlsruher Institut für Technologie
Fakultät für Maschinenbau
Tag der mündlichen Prüfung: 23. März 2011
Impressum
ISBN 978-3-86644-678-6
Numerische Simulation turbulenter
reaktiver Strömungen mit einem
hybriden CFD/transported PDF Modell
Die vorliegenden Arbeit entstand während meiner Tätigkeit als akademischer Mitarbeiter
am Institut für Technische Thermodynamik (ITT) des Karlsruher Instituts für Technologie
(KIT). Zahlreiche Personen haben zum Gelingen dieser Arbeit beigetragen, denen ich allen
an dieser Stelle herzlich danken möchte. Wenn ich im folgenden einige wenige namentlich
erwähne mag doch niemand, der sich nicht in der Aufzählung findet und mich trotzdem in
welcher Form auch immer unterstützt und begleitet hat, sich geringgeschätzt fühlen. Ich
bin für jede auch noch so kleine Unterstüzung sehr dankbar. Es liegt aber leider in der
Natur der Sache, dass einzelne zu erwähnen immer bedeutet andere, die auch einen Beitrag
geleistet haben, unerwähnt zu lassen.
Als erstes sei meinem wissenschaftlichen Betreuer Prof. Dr. rer. nat. habil. Ulrich Maas für
die Möglichkeit zur Durchführung der Arbeit an seinem Institut gedankt. Die interssan-
te Aufgabenstellung, seine persönliche Unterstützung und die vielen interessanten wissen-
schaftlichen Diskusionen, die ich mit ihm führen durfte, haben entscheidend zum Gelingen
der Arbeit beigetragen.
Zusätzlich danke ich Prof. Dr.-Ing. habil. Andreas Class für die freundliche Übernahme des
Korreferats und die hilfreichen Hinweise zur Arbeit.
Allen akutellen und ehemaligen Kollegen am ITT danke ich für die freundschaftliche Ar-
beitsatmosphäre und die große Hilfsbereitschaft bei allen Problemen. Besonders erwähnen
möchte ich meinen langjährigen Zimmerkollegen Florian Zieker mit dem ich viele fachlich
und menschlich hoch interessante Diskusionen zu vielfältigen Aspekten der technischen Ver-
brennung und des Alltagslebens führen durfte. Ich danke ihm für das viele, dass ich dabei
fachlich und menschlich von ihm lernen konnte. Für die Hilfe und Unterstützung besonders
in der Anfangszeit als akademischer Mitarbeiter danke ich meinen ehemaligen Kollegen
Rainer Stauch und Alexander Schubert sehr herzlich. Fachlich hat die Arbeit insbesondere
in der zweiten Hälfte wesentlich von den Diskusionen und Anregungen von Gerd Steinhil-
ber profitiert. Hierfür sei ihm neben mancher anderer Unterstützung beim Thema PDF
Modellierung ein besonders großer Dank ausgesprochen.
Mark Scherr danke ich für die stets rasche und fundierte Hilfe bei allen LATEX Problemen,
die bei der Erstellung der schriftlichen Ausarbeitung dieser Arbeit aufgetreten sind. Dar-
überhinaus stand er mir gemeinsam mit Christian Hofrath auch immer wieder bei Unklar-
heiten, Fragen und Problemen zu der faszinierend Welt der Thermodynamik zur Verfügung.
Beiden sei dafür an dieser Stelle besonders herzlich gedankt.
Ebenfalls danke ich Peter Lammers (ehemals Höchstleistungsrechenzentrum Stuttgart) sehr
für seine Unterstützung bei der Portierung des Simulationscodes auf die NEC SX-8 am
HLRS und die sehr interessanten Diskussionen und Hinweise zum Thema Höchstleistungs-
rechnen.
I
II
Franco Magagnato vom Fachgebiet Strömungsmaschinen des KIT danke ich für die freund-
liche Überlassung des Strömungslösers SPARC und Unterstützung bei der Einarbeitung in
diesen Simulationscode.
Allen Studien- und Diplomarbeitern, die ebenfalls einen Teil zum Gelingen dieser Arbeit
beigetragen haben, sei auch herzlich gedankt.
Die Arbeit wurde finanziell durch die DFG im Rahmen des Paketprojekts „Flammenbe-
schleunigung in Wirbelröhren“ und den Sonderforschungsbereich 606 im Rahmen ver-
schiedener Teilprojekte gefördert, sowie durch das Höchstleistungsrechenzentrum Stuttgart
(HLRS) im Rahmen des vom Bundesminstrerium für Bildung und Forschung geförderten
Vorhabens „Teraflop Workbench“ durch die Bereitstellung von Rechenzeit unterstützt.
Abschließend sei noch ein besonderer Dank an meine Eltern ausgesprochen, die mich in
allem meinem Tun stets vorbehaltlos und liebevoll unterstützt haben und mich viel Gutes
und Nützliches für das Leben gelehrt haben.
Abbildungsverzeichnis VII
Tabellenverzeichnis XI
Abstract XIII
1 Einleitung 1
1.1 Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2
1.2 Phänomenologie . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
1.3 Zielsetzung . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
2 Turbulente Flammen 7
2.1 Phänomenologie turbulenter Strömungen . . . . . . . . . . . . . . . . . . . . 7
2.2 Charakterisierung turbulenter vorgemischter Flammen . . . . . . . . . . . . 10
2.3 Charakterisierung turbulenter nicht-vorgemischter Flammen . . . . . . . . . 12
III
IV INHALTSVERZEICHNIS
3.6.2 Verbundwahrscheinlichkeitsdichtefunktion . . . . . . . . . . . . . . . 34
3.6.3 PDF Transportgleichung . . . . . . . . . . . . . . . . . . . . . . . . . 35
3.6.4 Numerische Lösung der PDF Transportgleichung . . . . . . . . . . . 37
3.6.5 Modellierung der ungeschlossenen Terme . . . . . . . . . . . . . . . . 39
3.6.6 Zusammenfassung der Gleichungen auf Partikelebene . . . . . . . . . 50
3.7 Beschreibung der Reaktionskinetik . . . . . . . . . . . . . . . . . . . . . . . 51
3.7.1 Bruttoreaktionsmechanismen . . . . . . . . . . . . . . . . . . . . . . 52
3.7.2 Detaillierte Elementarreaktionsmechanismen . . . . . . . . . . . . . . 53
3.7.3 Reduzierte reaktionskinetische Modelle . . . . . . . . . . . . . . . . . 55
Literaturverzeichnis 143
Abbildungsverzeichnis
VII
VIII ABBILDUNGSVERZEICHNIS
XI
Abstract
The mathematical modelling of turbulent flames is a challenging task due to the intense
coupling between turbulent transport processes and chemical kinetics. In this thesis a
mathematical model focused upon the turbulence-chemistry interaction is presented. This
model is intended to be generic in terms of being applicable to different combustion regimes
(e.g. premixed and non-premixed flames, low and high turbulence levels, steady and non-
steady flows, flames determined either by chemical or physical time scales) and predictive.
In the scientific community as well as in industrial applications an increasing demand
for reliable predictive models can be observed. Such models need to treat combustion
chemistry without the use of oversimplified kinetic models. This is vitally important for
the improvement of industrial combustion devices like gas turbines or internal combustion
engines towards an improved fuel efficiency and new low emission goals. The standard
models for non-reacting flows (RANS/LES) cannot satisfactorily tackle the problem of
the strong non-linearity of the chemical source term, which is often modeled by making
use of oversimplified closure assumptions. Additionally these models often suffer from a
poor modelling of the turbulence-chemistry interaction. Since probability density functions
(PDF) allow a detailed modelling of this effect they are used in the model presented in this
work. Showing a high capability for modelling turbulent combustion processes by treating
convection and finite-rate chemistry exact only the effect of molecular mixing has to be
modelled.
In the presented model is a hybrid CFD/transported PDF scheme is derived. In order
to get the mean velocity and pressure fields together with the turbulent time scale the
compressible Navier-Stokes equations are solved by a finite volume solver. This solver is
coupled with a PDF solver to optain the thermokinetic state of the fluid. Chemical kinetics
are taken into account by using automatically reduced chemical reaction mechanisms. The
detailled mechanisms are reduced with the recently developed REDIM method which allows
a sound description of the reaction progress using only a very small number of variables
(two in this case) and already accounts for the effect of molecular transport in state space.
The model was validated with a broad experimental data base taken from literature. First
of all a statically stationary non-premixed and a statistically stationary premixed flame were
investigated. Subsequently the model was used to simulate a flame flashback phenomenon
in a swirl stabilized lean premixed flame.
XIII
XIV ABSTRACT
For both the premixed and the non-premixed flame the results of the model validation sim-
ulations show a good agreement of the experimental data and the results of the simulations.
Spacial profiles of velocity, mixture fraction, temperature and different species depending
on the available experimental data base are studied. The global operation parameters of
both flames cover a wide range of Reynolds and Damköhler numbers and thus a wide range
of turbulence intensities and ratios of physical and chemical time scales.
Based upon the validation calculations the model was used for analysis and prediction of
flame flashback in lean premixed combustion devices. The device investigated was a sim-
plified model of a real combustion chamber used within an industrial gas turbine. The
flame is stabilized by a swirling main flow and a backward facing step at the intersection
of the premixing duct and the combustion chamber. The occuring flashback mode differs
notably from many modes already broadly discused in the literature. It is named com-
bustion induced vortex breakdown (CIVB). Such a flashback is initiated by enriching the
mixture. Detailed analysis of a single flashback event reveals that enriching first leads to
constriction of a small backflow bubble from the main backflow area. At further enrichment
this bubble travels upstream transporting the flame into the premixing duct. As driving
mechanism the negative azimutal vorticity which is produced at the tip of the backflow
bubble was identified. The negative azimutal vorticity induces a momentum in negative
axial direction on the backflow bubble. This term is determined by the cross product of
the density and the pressure gradient. Parameterical studies should clarify the influence of
two global operation parameters namely the premixing temperature and the total mass flux
upon the flashback behaviour of the flame. Two general trends where observed. An increase
in the premixing temperature shifts the flashback point towards leaner mixtures due to the
increased gradient at the flame front by the increasing preheat temperature. Larger inlet
mass fluxes shift the flashback point towards richer mixtures. The higher inlet mass flux
leads to a higher momentum in positive axial direction which has to be outbalanced by
the momentum in negative axial direction created by the negative azimutal vorticity a the
flame tip.
Kapitel 1
Einleitung
Die Entdeckung und Nutzbarmachung des Feuers in der Geschichte der Menschheit reicht
nach heute gesicherten wissenschaftlichen Erkenntnissen schon nahezu eine Million Jah-
re zurück [51] und kann sogar bei den Vorfahren des modernen Menschen nachgewiesen
werden [52]. Der Einsatz der Verbrennung als Technologie findet sich schon früh in nahe-
zu allen Kulturen weltweit wieder. Über die Jahrhunderte und Jahrtausende hinweg hat
sich jedoch die Art des Einsatzes stark verändert. In den frühen Kulturen dienten die offe-
nen Feuerstellen hauptsächlich der Aufbereitung von Nahrungsmitteln und des Heizens von
Wohnbereichen oder Lagerstätten. Heute sind die Einsatzgebiete weit vielfältiger. Verbren-
nungsprozesse in Kolbenmaschinen oder Flugtriebwerken dienen der Fortbewegung oder
dem fördern und pumpen von Fluiden, Verbrennungsprozesse in Kraftwerken werden zur
Erzeugung von elektrischem Strom eingesetzt, dezentrale Haushaltsbrenner sorgen für die
Versorgung von Haushalten mit Warmwasser und Wärme und in großtechnischen Feue-
rungsanlagen werden unter anderem Haushalts- und Industrieabfälle thermisch behandelt.
Diese Aufzählung ließe sich noch nahezu beliebig ergänzen.
Problematisch am Einsatz von Verbrennungsprozessen sind im wesentlichen drei Dinge.
Zum einen verursachen Verbrennungsprozesse häufig eine ungewollte Emission von Schad-
stoffen in die Atmosphäre. So kann es durch unvollständige Verbrennung zur Emission
von auf den Menschen toxisch wirkendem Kohlenmonoxid oder zum Ausstoß von toxischen
oder cancerogenen unverbrannten Kohlenwasserstoffen kommen. Zusätzlich können durch
Verunreinigungen im Brennstoff weitere Schadstoffe entstehen. Als Beispiel seien hier die
durch Schwefeleinlagerungen in Kohle oder Benzin entstehenden Emissionen von Schwe-
feloxiden genannt. Weiterhin sind die Resourcen der eingesetzten fossilen Brennstoffe wie
Kohle, Erdöl oder Erdgas endlich und fordern, dass technische Verbrennungsprozesse immer
weiter hinsichtlich ihrer Effizienz hin optimiert werden müssen. Das letzte Problem stellen
die bei der Verbrennung fossiler Brennstoffe unvermeidlichen Kohlendioxidemissionen dar.
Aktuelle Studien und wissenschaftliche Arbeiten legen nahe, dass die von Menschen verur-
sachten Kohlendioxidemissionen Auswirkungen auf das globale Klima der Erde haben [60].
1
2 KAPITEL 1: EINLEITUNG
Die Auswirkungen der Klimaänderung können in weiten Teilen eine sehr negative Wirkung
auf die Lebensbedingungen und Lebensräume vieler Menschen weltweit haben.
Vor diesem Hintergrund kommt der Forschung im Bereich technischer Verbrennung in Wis-
senschaft und Technik heute und auch in absehbarer Zukunft eine wichtige Rolle zu. Weitere
Optimierung bereits etablierter Verbrennungssysteme fordert einen immer tiefgehenderen
Einblick in die in diesen Systemen ablaufenden sehr komplexen physikalischen und chemi-
schen Prozesse.
1.1 Motivation
Um einen tiefgehenden Einblick in die physikalischen und chemischen Prozesse während
einer Verbrennung erreichen zu können, sind detaillierte Untersuchungen und Analysen not-
wendig. Hierbei ist sowohl an die Analyse von mittleren oder instantanen Geschwindigkeits-
oder Temperaturfeldern sowie der Konzentration verschiedener Spezies, als auch an die Un-
tersuchung verschiedener äußerer Einfluss- und Betriebsparameter zu denken.
Die numerische Simulation kann dabei eine wertvolle Ergänzung zu experimentellen Un-
tersuchungen darstellen. Sie ermöglicht oftmals einen tieferen Einblick in den Charakter
physikalischer und chemischer Prozesse und deren Kopplung als dies bei experimenellen
Untersuchungen aufgrund von messtechnischen Schwierigkeiten erreichbar ist. Man denke
nur an die Problematik des Designs eines nicht inversiven Experiments oder die Schwierig-
keiten bei der Schaffung einer optischen Zugänglichkeit des zu untersuchenden Objekts. Die
numerische Simulation ermöglicht eine detaillierte Untersuchung der Interaktion zwischen
chemischer Kinetik und physikalischen Transportprozessen. Diese Untersuchungen sind un-
erlässlich für das theoretische Verständnis von Transportprozessen und somit notwendig
zur Neu- und Weiterentwicklung bestehender Modelle und technischer Systeme. Zusätzlich
lassen sich im Rahmen von Simulationen äußere Einflussfaktoren ohne übermäßigen Mehr-
aufwand gezielt und separat variieren. Sensitivitätsanalysen erlauben darüberhinaus eine
detaillierte Analyse des äußeren Einflussfaktoren. Neben den diskutierten Vorteilen bietet
die numerische Simualtion noch eine Reihe weiterer Vorteile [86]. Nichts desto trotz stellt
die numerische Simulation keinesfalls einen Ersatz für experimentelle Untersuchungen dar,
sondern sie kann vielmehr als eine sinnvolle Ergänzung dieser betrachtet werden. So erfolgt
beispielsweise die Verifizierung mathematischer Modelle und numerischer Simulationen im-
mer durch den Vergleich der Resultate mit experimentellen Ergebnissen.
In vielen technischen Anwendungen besteht ein großer Bedarf nach zuverlässigen prädik-
tiven Modellen turbulenter Flammen [56, 73]. Beispielhaft sei hier nur die Entwicklung
moderner schadstoffarmer Brennerkonzepte für Gasturbinen [8] und aktuelle Forschungsar-
beiten für verbrauchsarme Verbrennungsmotoren [1, 7, 142] und optimierte Großfeuerungs-
anlagen [147] angeführt. Gekennzeichnet sind alle diese Anwendungsfälle durch die komplexe
Wechselwirkung von Turbulenz und chemischer Kinetik, die Maßgeblich das Erscheinungs-
1.2 PHÄNOMENOLOGIE 3
1.2 Phänomenologie
Die Vorgänge in turbulenten Flammen sind gekennzeichnet durch eine starke Interaktion
von Strömung, Turbulenz und Chemie. Sie spielen sich zudem auf einem weiten Bereich
von Zeit- und Längenskalen ab. So kommt es auf den großen Längenskalen zur Produkti-
on von Turbulenz. Der hierfür relevante Längenbereich liegt im Bereich charakteristischer
makroskopischer Abmessungen von einigen Millimetern bis in extremen Fällen hin zu eini-
gen Metern. Die Dissipation findet auf der Kolomogorovschen Längenskala statt. Es gibt
sehr schnelle chemische Reaktionen wie zum Beispiel Radikalreaktionen und sehr langsame
Reaktion wie beispielsweise die Bildung von Stickoxiden. Hinzukommt, dass die Zeit- und
Längenskalen wie in Abbildung 1.1 zu erkennen zusätzlich noch in einem realen Verbren-
nungsprozess stark überlappen.
Die Standardverfahren für nicht-reagierende Strömungen (RANS, LES) können das Pro-
blem der starken Nichtlinearität des chemischen Quellterms nicht zufriedenstellen behan-
deln und die Qualität der Ergebnisse leidet oftmals unter einer mangelhaften Modellierung
der Turbulenz-Chemie Interaktion, wohingegen zum Beispiel durch die Anwendung von
Wahrscheinlichkeitsdichtefunktionsmethoden (PDF) eine detaillierte Modellierung dieses
4 KAPITEL 1: EINLEITUNG
10 0 s 10 0 s
−2
10 s 10 −2s
Strömung
−4 −4
10 s Transport 10 s
10 −6s 10 −6s
−8
10 s 10 −8s
Abbildung 1.1: Charakteristische Zeitskalen der physikalischen und chemischen Prozesse sowie
typische Zündverzugszeiten und motorische Zeitskalen [159]
Effekts möglich ist. PDF Methoden zeigen eine hohe Leistungsfähigkeit bei der Modellie-
rung turbulenter Flammen, da sie sowohl den Einfluß von Konvektion als auch chemische
Reaktionen exakt behandeln können [124, 126] und nur der Einfluß von molekularer Mi-
schung modelliert werden muss [98, 99, 132].
Die in der Literatur zu findenden PDF Modelle lassen sich grob in zwei Klassen unterglie-
dern. Einige Autoren nutzen ein „stand alone“-PDF Verfahren, bei welchem alle hydrodyna-
mischen und thermokinetischen Eigenschaften des Fluids aus einer Verbundwahrscheinlich-
keitsdichtefunktion (JPDF) berechnet werden (z.B. [82, 82, 133, 135, 153]). Andere Autoren
wiederum bestimmen nur einen bestimmten Teil des hydrodynamischen und thermokineti-
schen Zustandsvektors aus einer JPDF und berechnen die übrigen Eigenschaften mit einem
gewöhnlichen CFD Löser für die Navier-Stokes Gleichungen (z.B. [10, 63, 104, 167]).
Grundsätzlich ist die Form der PDF in einem turbulenten Verbrennungsprozess a priori un-
bekannt. Ein allgemeingültiger Ansatz zur Bestimmung der PDF besteht in der Lösung einer
Transportgleichung für die PDF selbst. Auf dieses Verfahren wird auch in der vorliegenden
Arbeit zurückgegriffen. Vereinfachende Anätze mit einer a priori Formannahme (z.B. eine
Betafunktion) für die PDF selbst und dem anschließenden ausschließlichen Bestimmen der
Momente der Verteilungsfunktion kommen wegen ihrer nicht allgemeinen Gültigkeit nicht
zum Einsatz. Die Transportgleichung für die PDF kann aus den Navier-Stokes Gleichungen
abgeleitet werden [124]. Wie bereits oben erwähnt tritt sowohl der chemische Quellterm
als auch der Term der die Volumenkräfte und den Einfluss des mittleren Druckgradienten
beschreibt in geschlossener Form auf. Lediglich der Term der den Effekt des fluktuieren-
den Druckgradienten und des molekularen Mischens beschreibt muss modelliert werden.
Numerisch unterscheidet sich die Behandlung der PDF Transportgleichung deutlich von
den Navier-Stokes Gleichungen. Im Gegensatz von dem System partieller Differentialglei-
chungen welches die Navier-Stokes Gleichungen darstellen, ist die Transportgleichung der
PDF eine multidimensionale skalare Erhaltungsgleichung. Sie besitzt im allgemeinen 7 + ns
unabhänige Variablen. Dies sind drei Dimenisionen im Ortsraum, drei Dimensionen im
1.3 ZIELSETZUNG 5
Geschwindigkeitsraum, die Zeit und die Zahl der Skalare ns die notwendig sind um den
thermochemischen Zustand des Fluids zu beschreiben. Wegen ihrer hohen Dimensionalität
ist es (zumindest in den allermeisten Fällen) nicht praktikabel die Gleichung mit Finite Vo-
lumen oder Finite Differenzen Verfahren zu lösen. Aus diesem Grund werden Monte Carlo
Verfahren verwendet, die in der Numerik dort weite Verbreitung finden wo Probleme hoher
Dimensionalität zu lösen sind. Ihr wesentlicher Vorteil ist, dass der numerische Aufwand
nur linear mit der Anzahl an Dimensionen zunimmt.
Das Lösungsverfahren nutzt eine spezielle Eigenschaft der PDF. Die kontinuierliche PDF
kann durch eine diskrete Verteilungsfunktion angenähert werden. Die diskrete Verteilungs-
funktion setzt sich aus verschiedenen möglichen Einzelrealisierungen der Strömung zusam-
men, deren zeitliche und räumliche Entwicklung durch stochatische Prozesse beschrie-
ben wird. Jede dieser Einzelrealisierungen stellt ein sogenanntes stochatisches Partikel
dar. Durch ein Ensemble solcher stochastischer Partikel kann die PDF rekonstruiert wer-
den [122]. Die Transportgleichung der PDF wird in diesem Zusammenhang in ein System
gewöhnlicher (stochastischer) Differentialgleichungen transformiert. Die verfügbaren Mon-
te Carlo Verfahren verursachen ein starkes Rauschen beim Lösen dieses Differentialglei-
chungssystems, welches zu Stabilitätsproblemen bei der Bestimmung des mittleren Druck-
gradienten führen kann. Eine Möglichkeit diese Stabilitätsprobleme zu beseitigen ist die
Partikelmethode mit einem konventionellen Finite Volumen Löser zu koppeln um den mitt-
leren Druckgradient aus den Navier-Stokes Gleichungen zu berechnen. Diese sogenannten
hybriden PDF/CFD Verfahren wurde bereits von verschiedenen Autoren für eine Vielzahl
verschiedener Flammentypen eingesetzt (z.B. [62, 103, 131, 151]).
In der vorliegenden Arbeit wird ebenfalls ein hybrides Lösungsschema benutzt. Die kom-
pressiblen Navier-Stokes Gleichungen werden mit einem Finite Volumen Löser gelöst um die
mittleren Geschwindigkeiten und den mittleren Druckgradienten zu erhalten. Dieser Löser
wird mit einem PDF Verfahren gekoppelt um den thermokinetischen Zustand des Fluids
bestimmen zu können. Der Einfluss der chemischen Kinetik wird durch automatisch re-
duzierte detaillierte Reaktionsmechanismen erfasst. Der detaillierte Reaktionsmechanismus
wird mit dem kürzlich entwickelten REDIM Verfahren reduziert [24]. Dieses erlaubt eine
zuverlässige Beschreibung der chemischen Kinetik bereits mit einer sehr geringen Zahl von
Parametern. Im vorliegenden Fall sind dies zwei Parameter. REDIM berücksichtigt dabei
bereits den Einfluss von molekularem Transport im Zustandsraum.
1.3 Zielsetzung
Zielsetzung dieser Arbeit ist die Entwicklung und Validierung eines Gesamtmodells zur Si-
mulation turbulenter Verbrennungsprozesse. In der Literatur sind eine große Zahl solcher
Modelle zu finden. Wesentliches Augenmerk soll deshalb auf hochwertigen Modellen für
die Wechselwirkung von Chemie und Turbulenz sowie auf der zuverlässigen Beschreibung
6 KAPITEL 1: EINLEITUNG
Turbulente Flammen
• Sie beinhalten eine große Menge an Wirbelstärke. Die Streckung von Wirbeln ist ein
wesentlicher Mechanismus durch den zusätzliche Turbulenz produziert wird.
7
8 KAPITEL 2: TURBULENTE FLAMMEN
• Turbulenz erhöht das Maß in dem passive oder reaktive Skalare miteinander vermischt
werden. Dies geschieht dadurch, dass Fluidpakete mit unterschiedlicher Konzentra-
tion der Skalare miteinander in Kontakt gebracht werden. Die eigentliche Mischung
geschied zwar durch Diffusion, nichts desto trotz wird diese Eigenschaft turbulenter
Strömungen dennoch oftmals als „diffusiv“ bezeichnet.
• Durch Erhöhung der Mischung bringt die Turbulenz Fluidpakete in Kontakt miteinan-
der, die einen stark unterschiedlichen Impuls haben. Die dadurch bedingte Reduktion
der durch die Viskosität des Fluids verursachten Geschwindigkeitsgradienten veringert
die kinetische Energie der Strömung. Mit anderen Worten ausgedrückt ist die Strö-
mung „dissipativ“. Die kinetische Energie wird hierbei irreversibel in innere Energie
des Fluids umgewandelt.
Alle diese Phänomene sind wichtig. Die durch Turbulenz verursachten Effekte können in
einem bestimmten Anwendungsfall gewollt oder ungewollt sein. Eine intensive Mischung
zweier Skalare kann gewünscht sein bei chemischen Reaktion oder in Strömungen bei wel-
chen ein großer Wärmeübergang erzielt werden soll. Auf der anderen Seite führt eine ver-
stärkte Mischung des Impulses zu zusätzlichen Reibungskräften und somit ist die Leistung,
die zum Beispiel notwendig ist, um ein Fluid zu pumpen höher oder der Widerstand eines
umströmten Körpers größer.
In einer turbulenten Strömung ist das Strömungsfeld durch starke räumliche und zeitliche
Fluktuationen der hydrodynamischen Größen gekennzeichnet. Ob eine Strömung unter ge-
gebenen Randbedingungen turbulent sein wird läßt sich anhand einer charakteristischen
dimensionslosen Kennzahl, der Reynoldszahl abschätzen. Übersteigt sie einen gewissen
Schwellenwert ist die Strömung turbulent. Die absolute Größe des Schwellenwerts hängt
zusätzlich noch vom betrachteten Strömungsproblem (Grenzschichtströmung, Freistrahl,
Innenströmung, etc.) ab [136]. Berechnet wird sie aus dem Verhältnis von Trägheits- und
Reibungskräften. Die Reynoldszahl
U ·L
Re = (2.1)
ν
hängt von einer charakteristischen Geschwindigkeit U , einer für das betrachtete Strömungs-
problem charakteristischen Länge L und der kinematischen Viskosität des Fluids ab. Als
charakteristischen Länge dient zum Beispiel für eine Rohrströmung der Rohrdurchmessen
2.1 PHÄNOMENOLOGIE TURBULENTER STRÖMUNGEN 9
und für einen Strömung über eine ebene Platte die Lauflänge der Grenzschicht [136]. Zum
Umschlag von laminarer und turbulenter Strömung kommt es, wenn die Trägheitskräfte
deutlich größer als die Reibungskräfte sind. Die viskosen Kräfte reichen dann nicht mehr
dazu aus, die destabilisierenden Trägheitskräfte zu kommpensieren und es kommt zu einer
ungeordneten, stochatisch schwankenden Strömung in welcher sich Wirbel bilden. Der gan-
ze Prozess wird durch kleinste Störungen initiert, die in der Natur stets vorhanden sind.
Die Wirbel bilden sich durch die senkrecht zur Hauptströmung entstehenden Geschwindig-
keiten. Der Stofftransport und Impulstransport innerhalb des Fluids wird so intensiviert.
Die dreidimensionalen Wirbel werden sowohl durch Scherung der Strömung als auch durch
benachbarte Wirbel ähnlicher oder größerer Abmessungen verformt. Hierdurch verkleinern
sich ihre charakteristischen Abmessungen und sie zerfallen letztendlich. Die Energie wird
aufgrund der Drehimpulserhaltung von den gößeren auf die kleinern Wirbel weitergeben.
Die Umdrehungsgeschwindigkeit der einzelnen Wirbel nimmt dabei immer mehr und mehr
zu. Ab einen bestimmten Punkt wird der Einfluß der Viskosität bemerkbar. Die kleinsten
Wirbel beginnen zu dissipieren und geben ihre Energie in Form von Wärme an das Fluid
ab. Die turbulente kinetische Energie wird folglich auf den großen Längenskalen produziert
und auf den kleinen Längenskalen dissipiert [12, 46].
Die Weitergabe der Turbulenzenergie von großen Wirbeln an immer kleiner werdende Wir-
bel veranschaulicht die in Abbildung 2.1 dargestellte Turbulenzkaskade. In diesem Dia-
gramm ist das Energiespektrum E(κ, , ν) über der Wellenzahl κ = 2πl (Wirbel pro Län-
geneinheit) aufgetragen. Sie zeigt wie die turbulente kinetische Energie über Wirbel ver-
6E(κ, , νl )
Wellenzahl
κ = 2πl
-
Produktions- energietragende Inertial-/Trägheits- Dissipationsbereich
bereich Wirbel lt bereich Kolmogorov-Maß lη
Abbildung 2.1: Normierte, spektrale Energieverteilung E(κ, , νl ) als Funktion der Wellenzahl
κ (nach [75])
10 KAPITEL 2: TURBULENTE FLAMMEN
schiedener Größe aufgeteilt ist. Während der Produktion von Turbulenz steigt die Energie
der Wirbel immer weiter an, da sich aus kleinern Wirbelstrukturen große Wirbel bilden.
Diese energietragenden Wirbel geben über den Kaskadenprozess ihre Energie an immer
kleinere Wirbel weiter. Innerhalb dieses Inertialbereichs findet ausschließlich Transfer von
Turbulenzenergie von großen auf kleinere Strukturen statt. Bei sehr kleinen Wirbeln be-
ginnt schließlich der Dissipationsbereich. Hier lösen sich die sehr kleinen Wirbel auf und
geben ihre kinetische Energie in Form von Wärme an das Fluid ab dessen innere Energie
dadurch zunimmt [78, 130, 150].
Abbildung 2.2: Regime einer turbulenten vorgemischten Verbrennung nach Borghi [16]
2.2 CHARAKTERISIERUNG TURBULENTER VORGEMISCHTER FLAMMEN 11
und beschreibt das Verhältnis aus einem turbulenten Zeitmaß τt für die Strömung und
einem charakteristischen Zeitmaß τC für die chemischen Reaktionen [159]. Die turbulente
Karlowitzzahl Kat bestimmt sich aus dem Quotienten eines chemischen Zeitmaßes τC und
dem Kolmogorovzeitmaß τη zu
τC τC 1
Kat = ≈ Ret2 (2.4)
τη τt
wobei das Kolmogorovzeitmaß mit der Umdrehungsdauer der kleinsten in der Strömung
vorhandenen Wirbel korrespondiert.
Mit Hilfe dieser dimensionslosen Kennzahlen lassen sich beliebige Flammen in das Borghi-
diagramm einordnen. Die einzelnen Regime unterscheiden sich durch die Feinstruktur der
Flammenfront. Die Zuordung der einzelnen Regime zu den jeweiligen Zahlenwerten von
Dat und Kat kann unmittelbar aus der Abbildung 2.2 entnommen werden. Die einzelnen
Regime sind dabei durch die folgenden Phänomene gekennzeichnet.
• Laminare, leicht gewellte Flammenfronten (u < uL )
Durch die gegenüber der Flammengeschwindigkeit kleinen Geschwindigkeitsfluktua-
tionen und den im Vergleich zu der Flammenfrontdicke größeren Turbulenzelementen
wird die Flammenfront nur wenig beeinflußt. Sie ist gewellt und wirkt im zeitlichen
Mittel verdickt („corrugated flamelets“).
• Stark gewellte Flammen mit Inselbildung (u > uL und Kat < 1)
Verursacht durch die zunehmende Schwankungsgeschwindigkeit kommt es zu einer
Auffaltung der Flammenfront. Die Flammenfront bleibt aber gegenüber den Wirbe-
labmessungen dünn. Durch die verstärkte Krümmung der Flammenfront kann es hier
zu lokalem Verlöschen und Wiederzünden kommen. Wodurch sich Inseln von Frisch-
gemisch und Rauchgas abschnüren können („wrinkled flamelets with pockets“).
kleinen Abmessung in diese eindringen und dort den diffusiven Austausch verstärken.
Dies führt zu einer Verdickung der Flammenfront. Bis hin zu Dat ≥ 1 können immer
größere Wirbel in die Flammenfront eindringen, wodurch es lokal zu Verlöschungen
bis hin zum Zerreißen der Flammenfront kommen kann („distributed reaction zones“).
Angemerkt sei, dass sich der Zustand in einer realen vorgemischten Flamme in der Regel
nicht durch einen Punkt im Borghidiagramm beschreiben läßt, da die Längen- und Zeit-
skalen stark variieren können. Somit findet man in realen Flammen meist auch ein gewisses
Spektrum an verschiedenen turbulenten Karlowitz- und Damköhlerzahlen vor. Der daraus
im Borghidiagramm entstehende Bereich kann dabei durchaus auch verschiedene Regime
beinhalten, wodurch es zu mehreren der oben erwähnten Effekte innerhalb derselben Flam-
me kommen kann.
Abbildung 2.3: Regime einer turbulenten nicht-vorgemischten Verbrennung nach Borghi [17]
dem Diffusionskoeffizienten D. Unter der Annahme, dass nur eine chemische Zeitskala exi-
stiert und das Verhältnis der Temperaturleitfähigkeit a und des Diffusionskoeffizienten D
gleich eins ist (Lewiszahl Le = Da = 1), ergibt sich die laminare Flammendicke zu
δL = D · τC (2.5)
Die laminare Flammengeschwindigkeit bestimmt mit welcher Geschwindigkeit sich eine un-
gestreckte Flammenfront durch ein vorgemischtes Brennstoff/Luft Gemisch bewegt. Auch
sie ist proportional zu dem Diffusionskoeffizienten und der chemischen Zeitskala und ergibt
sich mit der erwähnten Vereinfachungen zu
D
uL = (2.6)
τC
Die turbulente Reynoldszahl wird analog zu Gleichung 2.2 berechnet. Laminare Flammen
befinden sich im Diagramm 2.3 unterhalb der Linie Ret = 1. Eine turbulente Flamme wird
durch drei Zeitskalen in die verschiedenen Flammenregime eingeordnet. Im einzelnen sind
dies: die Zeitskala chemischer Reaktion τC , die kolomogorovsche Zeitskala τK , durch die
Umdrehungszeit der kleinsten Wirbel repräsentiert wird und die turbulenten Zeitskala τt ,
die die Umdrehungsdauer der größten Wirbel darstellt.
Die einzelnen Regime sind durch verschiedene charakteristische Phänomene gekennzeichnet,
die im folgenden kurz dargestellt werden.
14 KAPITEL 2: TURBULENTE FLAMMEN
15
16 KAPITEL 3: NUMERISCHE SIMULATION REAKTIVER STRÖMUNGEN
Sg
r
n
qg
uur
¶W Fg
W
Rand ∂Ω begrenzt. Eine beliebige Zustandgröße läßt sich zum Zeitpunkt t als das Volu-
menintegral
G (t) = g (x, t) (3.1)
Ω
über die zugehörige Dichteverteilung g(x, t) mit x als Ortsvektor beschreiben. Eine Ände-
rung der extensiven Größe G(x, t) kann in einer reaktiven Strömung durch drei Prozesse
erfolgen (n bezeichnet den Normalenvektor des Randes ∂Ω des Volumenelementes Ω; dV
und dS jeweils ein Volumen- bzw. ein Flächenelement):
g · n · dS durch die Oberfläche (z.B. durch Diffusion, Wärmeleitung
• durch einen Fluß φ
oder Reibung),
• durch Produktion q̇g im Inneren des Volumenelementes (z.B. durch chemische Reak-
tion),
Die Berücksichtigung all dieser Prozesse führt für die zeitliche Änderung der Größe G zu
folgendem integralen Ausdruck
∂g
· dV + φg · n · dS = q̇g · dV + sg · dV (3.2)
∂t
Ω ∂Ω Ω Ω
der mit Hilfe des Gauss’schen Integralsatzes [20] in differentieller Form dargestellt werden
kann.
∂g
+ div φg = q̇g + sg (3.3)
∂t
Mit dieser allgemeinen Formulierung lassen sich die speziellen Erhaltungsgleichungen für
Masse, Impuls, Energie und Spezieserhaltung mathematisch ableiten. Auf eine detaillierte
Herleitung soll an dieser Stelle verzichtet werden. Sie findet sich an verschiedenen Stellen
in der Literatur. Genannt seien hier [11, 86, 109].
3.1 NAVIER-STOKES GLEICHUNGEN 17
Massenerhaltung
Speziesmassenerhaltung
Impulserhaltungsgleichung
In der Impulserhaltungsgleichung
∂ (ρ · v )
+ div (ρ · v ⊗ v ) + div p̄¯ = ρg (3.6)
∂t
gibt es ebenfalls einen konvektiven Term div (ρ · v ⊗ v ). Mit ⊗ wird das dyadische Produkt
zweier Vektoren bezeichnet [20]. Hinzu kommt der durch Druck- und Reibungskräfte ver-
¯ wobei p̄¯ den Drucktensor bezeichnet. Im Unterschied zu den bisher
ursachte Anteil div p̄,
betrachteten Erhaltungsgleichungen tritt in der Impulserhaltungsgleichung die Fernwirkung
ρg auf, welche den Einfluß von Gravitation oder einer anderen äußeren Beschleunigung be-
schreibt.
Energieerhaltungsgleichung
Einfluß von Fernwirkungen wie zum Beispiel Strahlung wurde in obiger Gleichung bereits
vernachlässigt. Strahlungseffekte spielen bei den in dieser Arbeit untersuchten Flammen
nur eine untergeordnete Rolle. Dies liegt an der sehr geringen Rußbildung, den moderaten
Temperaturen magerer Flammen und dem Mangel äußerer Strahlungsquellen oder starkt
Strahlung absorbierender Flächen. Die Energieerhaltungsgleichung kann somit quellterm-
frei formuliert werden. Die Gesamtenergie setzt sich aus der inneren Energie, der kinetischen
Energie und der potentiellen Energie zusammen. Der Beitrag der potentiellen Energie zur
Gesamtenergie ist für die untersuchten Fälle sehr gering und kann deshalb in guter Nähe-
rung vernachlässigt werden. Die als Gesamtenergie verwendetet Größe ist also nur noch die
Summe aus innerer und kinetischer Energie des Fluids.
1 1
ρe = ρu + ρ |v |2 + ρg = ρu + ρ |v |2 (3.8)
2
2
=0
Zustandsgleichung
ihren Mittelwert ḡ und eine Schwankung g gegenüber diesem Mittelwert aufzuteilen. Dieses
Mittelungsverfahren wird Reynoldsmittelung genannt.
g (x, t) = g (x, t) + g (x, t) (3.10)
Häufig wird als Mittelwert ein zeitlicher Mittelwert benutzt. Dies ist aber nicht zwingend
notwendig. Völlig analog könnte auch ein Ensemblemittelwert verwendet werden. Hier soll
nur die Verwendung eines zeitlichen Mittelwerts näher beschrieben werden. Der zeitliche
Mittelwert ist als
1
Δt
definiert. Das Mittelungsinterval Δt muss so groß gewählt werden, dass der Mittelwert
zeitunabhängig wird. Der zeitliche Mittelwert einer Schwankungsgröße wird dann zu Null.
g = 0 (3.12)
Die Anwendung dieses Mittelungsverfahrens auf die abhängigen Variablen der Navier-Stokes
Gleichungen führt zu den sogenannten Reynoldsgleichungen. Dabei entstehen neben den
Mittelwerten der abhänigen Variablen neue ungeschlossene Kreuzkorrelationsterme sowohl
in der Impulserhaltungsgleichung als auch in der Energieerhaltungsgleichung. Diese Terme
(in der Impulserhaltungsgleichung in der Form ui uj ) beinhalten den Einfluss der Turbulenz
auf die mittleren Größen des Strömungsfeldes und müssen modelliert werden. Details hierzu
finden sich in Kapitel 3.3.
Eine für Verbrennungsprozesse typische Eigenschaft ist das Auftreten großer Dichteschwan-
kungen. Dementsprechend müßte die Dichte ebenfalls in einen Mittelwert ρ und eine
Schwankungsgröße ρ zerlegt werden. Er erweist sich als zweckmäßig (siehe unten) einen
weiteren Mittelwert einzuführen, die Favre-Mittelung (dichtegewichteter Mittelwert, [41]).
Dieser ist für eine beliebige Strömungsgröße gegeben als
ρg
g̃ = bzw. ρ g̃ = ρ g (3.13)
ρ
Analog zu 3.10 läßt sich eine Größe wieder in Mittelwert und Schwankungsgröße bezüglich
des Favremittels aufspalten
g (x, t) = g̃ (x, t) + g (x, t) (3.14)
Mittelwert der Schwankung bei Reynolds-Mittelung gilt g = 0. Eine Gleichung für die
Umrechung des Reynolds-Mittelwertes einer Variablen in ihren Favre-Mittelwert abzuleiten,
ist durch einsetzen von Gleichung 3.13 in Gleichung 3.10 einfach möglich. Es folgt
ρg (ρ + ρ ) (g + g ) ρ g + ρ g + ρ g + ρ g
g̃ = = =
ρ ρ ρ
und daraus
ρ g
g̃ = g + . (3.16)
g
Hierzu muss jedoch die Korrelation der Schwankungen der Größe g und der Dichte bekannt
sein. Wesentlicher Vorteil der Verwendung der Favare-Mittelung ist, dass sie häufig eine
kompaktere Schreibweise gestattet. Dies sei am folgenden Beispiel des dichtegewichteten
Mittelwerts zweier Größen u und v veranschaulicht.
Während die Reynoldsmittelung für den Ausdruck ρ u v liefert
ρ u v = (ρ + ρ ) + (u + u ) + (v + v )
= ρ u v + ρ u v + ρ u v + ρ u v + ρ u v + ρ u v + ρ u v + ρ u v
= ρ u v + ρ u v + u ρ v + v ρ u + ρ u v (3.17)
jeweils für die Masse, die Speziesmasse1 , den Impuls und die Gesamtenergie des Systems.
Zum Schließen des Gleichungssystems wird noch eine thermische Zustandsgleichung benö-
tigt. Hierzu wird die gemittelte Zustandsgleichung eines idealen Gases in der Form
T
p=ρ·R· (3.23)
M
verwendet.
In den Erhaltungsgleichungen treten verschiedene ungeschlossene Terme auf. In der Ener-
gieerhaltungsgleichung ist dies der Term ρuj e und der Term uj p. Auf die Darstellung
möglicher Modellierungsansätze für diese Terme wird hier verzichtet, da bei dem in dieser
Arbeit verwendeten Gesamtmodell auf eine Lösung der Energiegleichung verzichtet werden
kann. Die Temperatur des Fluids wird direkt aus der Lösung der PDF Transportgleichung
ermittelt (näheres hierzu in Abschnitt 3.6.3 und Abschnitt 4).
In der Impulserhaltungsgleichung ist die Kreuzkorrelation ρui uj ungeschlossen. Sie stellt,
wie bereits oben erwähnt, den Einfluss der turbulenten Fluktuationen auf den mittleren
Geschwindigkeitsvektor dar. Dieser Term wird auch als turbulenter Scheinschubspannungs-
tensor bezeichnet. Zur Lösung des Gleichungssystems muß er modelliert werden. Hierzu
finden sich vielfältige Ansätze in der Literatur. Einige wenige werden in Abschnitt 3.3 nä-
her erläutert.
Die Spezieserhaltungsgleichung beinhaltet drei ungeschlossene Terme
Das wesentliche Problem bei der Mittelung des chemischen Quellterms ist, dass dieser
stark nichtlinear von der Temperatur und der Zusammensetzung abhängt. Mathematisch
ausgedrückt heißt das:
˜
Mα ω̇α = Sα (T, φ) = Sα (T̃ , φ) . (3.24)
μ t = ρ · u∗ · l ∗ (3.26)
Somit muss ein Modell gefunden werden, welches eine Bestimmung von u∗ und l∗ für die
lokalen Bedingungen in einer Strömung ermöglicht. Hierfür existieren in der Literatur ver-
schiedenste Modellierungsansätze, die sich nach der Zahl der zusätzlich zu lösenden Diffe-
rentialgleichungen in
• Eingleichungsmodelle und
2
δij = 1 für i = j und δij = 0 für i = j
3.3 STATISTISCHE MODELLIERUNG DER TURBULENZ 23
• Zweigleichungsmodelle
gruppieren lassen. So wird bei algebraischen Modellen (als typischer Vertreter wird das
Mischungswegmodell vorgestellt) l∗ auf Basis der Geometrie der Strömung spezifiziert. Bei
Ein- und Zweigleichungsmodellen werden zusätzliche Transportgleichungen für bestimmte
Turbulenzgrößen gelöst, aus denen sich das Produkt von u∗ und l∗ berechnen lässt. Ein in der
Literatur oft verwendeter Vertreter dieser Modellgruppen ist das k− Modell [64,65] bei ihm
werden u∗ und l∗ durch Lösung zweier zusätzlicher Transportgleichungen für die turbulente
kinetische Energie k und ihre Dissipationsrate bestimmt. Für alle drei oben genannten
Modellgruppen soll im folgenden ein typischer Vertreter näher beschrieben werden.
Zusätzlich existieren höherwertige Turbulenzmodelle, welche nicht auf der Bousinesq-
Annahme beruhen. Beispiele sind Reynoldsspannungsmodelle [72] und Large-Eddy Simula-
tionen [46]. Auch hierzu finden sich im folgenden einige kurze grundlegende Bemerkungen
(siehe Abschnitt 3.3.4).
3.3.1 Nullgleichungsmodelle
Nullgleichungsturbulenzmodelle (oder algebraische Turbulenzmodelle) haben den wesent-
lichen Vorteil einfach in bestehende numerische Verfahren integrierbar zu sein und wenig
rechenaufwand zu verursachen, da nur einfache algebraische Gleichungen und keine gewöhn-
lichen oder partiellen Differentialgleichungen gelöst werden müssen. Andererseits hängt die
turbulente Viskosität nur von lokalen (mittleren) Geschwindigkeitsprofilen ab. Ein Einfluss
des Zustands stromauf oder stromab der betrachteten Stelle wird nicht berücksichtigt, was
die Anwendbarkeit solcher Modelle stark einschränkt.
Ein bekannter Vertreter dieser Gruppe von Turbulenzmodellen ist das auf dem Prandlschen
Mischungswegansatz [130] beruhende Modell von Baldwin und Lomax [6]. Es wurde sehr
erfolgreich zur Simulation von Tragflügelumströmungen eingesetzt und soll grob beschrieben
werden. Für weitere Details sei auf die angegebene Literatur verwiesen.
Der Mischungswegansatz nach Prandtl bestimmt die turbulente Viskosität aus dem Ge-
schwindikeitsgradienten und einer noch zu bestimmenden Mischungslänge lm (x), die eine
Funktion des Ortes ist.
dũ
μt = ρ · lm ·
2
(3.27)
dy
Die Mischungsweglänge ist dabei definiert als diejenige Weglänge, die ein Fluidelement zu-
rücklegen kann bis es sich mit seiner Umgebung vollständig vermischt hat und somit gleich-
sam seine Identität verloren geht [130]. Im Baldwin-Lomax Modell wird die Mischungsweg-
länge mit der Gleichung
y+
lm = κ · y · 1 − exp − + (3.28)
A
24 KAPITEL 3: NUMERISCHE SIMULATION REAKTIVER STRÖMUNGEN
bestimmt. Darin sind κ eine Konstante mit dem Wert 0, 4, y der Wandabstand, y + ein
dimensionsloser Wandabstand und A+ eine Konstante mit dem Wert 26. Details zur Her-
leitung und physikalischen Begründung dieser Gleichung finden sich in [6].
In der allgemeinen Form ist das Mischungswegmodell auf alle Arten turbulenter Strömungen
anwendbar. Der größte Nachteil ist aber die Unvollständigkeit des Modells. Die Mischungs-
wegslänge lm muss bekannt sein und hängt wesentlich von der Geometrie der betrachteten
Strömung ab. Eine Abschätzung der Mischungsweglänge in einer komplexen Strömung ist
meist sehr schwierig und mit großen Fehleren behaftet wodurch die Qualität der Ergebnisse
meist sehr schecht sein wird. Anders sieht es aus bei technisch relevanten Strömungen, die
messtechnisch gut untersucht sind, so dass angepasste Korrelation für die Mischungsweglän-
ge existieren. Als Beispiel sei hier nochmals die turbulente Umströmung eines Tragflügels
oder andere Arten von klassischen Grenzschichtströmung genannt, bei welchen das Modell
sehr gute Resultate liefert.
3.3.2 Eingleichungsmodelle
Bei Eingleichungsmodellen wird die turbulente Viskosität durch Lösen einer zusätzli-
chen Transportgleichung bestimmt. Ein sehr verbreitetes Modell stammt von Spalart und
Allmaras [139], das eine zusätzliche Transportgleichung für die an die turbulente Viskosi-
tät angelehnte Hilfsgröße ν̃ einführt. Außer in Wandnähe stimmt ν̃ mit der turbulenten
Viskosität νt überein.
⎡ 2 ⎤
∂ (ρν̃) ∂ (ρν̃ui ) 1 ⎣ ∂ ∂ ν̃ ∂ ν̃ ⎦
+ = (μ + ρν̃) + Cb2 ρ
∂t ∂xi σν̃ ∂xj ∂xj ∂xj
2
ν̃
− Cω1 ρfω + Cb1 ρS̃ ν̃ (3.29)
d
Aus dieser Gleichung kann dann die turbulente Viskosität nach
μt = ρ · ν̃ · fν1 (3.30)
berechnet werden. Darin beschreibt fν1 eine viskose Dämpfung in der Nähe einer Wand und
kann durch Einführen eines dimensionslosen Parameters χ = ν̃ν bestimmt werden.
χ3
fν1 = 3
(3.31)
χ3 + Cν1
Die Modellkonstante Cν1 hat dabei den Wert 7, 1. Auf eine weitere Darstellung des Modells
und der zum Lösen der Gleichung 3.29 außdem notwendigen Konstanten und Modellpara-
meter wird hier verzichtet. Sie findet sich in [139].
3.3 STATISTISCHE MODELLIERUNG DER TURBULENZ 25
3.3.3 Zweigleichungsmodelle
Das meist verbreitete und eingesetzte Zweigleichungsmodell ist das Standard k − Modell
[64, 65]. Es löst eine Transportgleichung für die mittlere turbulente kinetische Energie
1 2
k̃ = · u + v 2 + w2 , (3.32)
2
die aus den Impulserhaltungsgleichungen abgeleitet werden kann und eine Transportglei-
chung für die mittlere Dissipationsrate ˜ der turbulenten kinetischen Energie. ˜ ist definiert
durch den Ausdruck
∂ui ∂ui
˜ = ν · (3.33)
∂xk ∂xk
mit ν = μρ als der laminaren kinematischen Viskosität. Die Transportgleichung selbst wird
dabei empirisch abgeleitet. Die beiden Gleichungen lauten
∂(ρk̃)
+ div(ρv k̃) − div(μt gradk̃) = Gk − ρ˜ (3.34)
∂t
∂(ρ˜) ˜
+ div(ρv ˜) − div(μt grad˜) = (C1 Gk − C2 ρ˜) . (3.35)
∂t k̃
Die turbulente Viskosität wird dann aus
k̃ 2
μt = Cμ · ρ · (3.36)
˜
berechnet. Der Quellterm der turbulenten kinetischen Energie ist durch die Gleichung
∂ ũi ∂ ũj
G k = Cμ ρ μ t + (3.37)
∂xj ∂xi
gegeben. C1 , C2 und Cμ sind empirische Modellkonstanten. Ihre Werte für das Standard-
modell [64] sind Tabelle 3.1 zu entnehmen.
Cμ 0,09
C1 1,44
C2 1,92
Auch in der vorliegenden Arbeit wurde ein Zweigleichungsmodell zur Modellierung der Tur-
bulenz verwendet. Zum Einsatz kam dabei nicht das Standard k − Modell sondern eine
Variante davon. Die verwendete Variante stammt von Speziale et al. [166] und löst zusätzlich
zur Transportgleichung für die turbulente kinetische Energie eine Transportgleichung für ein
turbulentes Zeitmaß τ . Das Modell berücksichtigt zusätzlich noch den Einfluss von Wänden
26 KAPITEL 3: NUMERISCHE SIMULATION REAKTIVER STRÖMUNGEN
auf die turbulente Viskosität. Das direkte Lösen einer Transportgleichung für ein turbulen-
tes Zeitmaß gleicht einen wesentlichen Nachteil des k − Modells aus. Für Epsilon ist die
Formulierung einer physikalisch sinnvollen Randbedingung an einer Wand sehr schwierig.
Unmittelbar ableiten läßt sich nur eine Randbedinung, die die zweite Ableitung der turbu-
lenten kinetischen Energie an der Wand beinhaltet. Die Anwendung dieser Randbedingung
kann zu einer hohen numerischen Steifigkeit des Problems führen [114]. Die turbulente Zeits-
kala hingegen geht an einer Wand gegen Null und somit kann für die Transportgleichung
einfach eine physikalisch begründetet Randbedingung abgeleitet werden [114, 166].
Die turbulente Viskosität ist bei diesem Modell durch die Gleichung
μt = Cμ ρ fμ k̃ τ̃ (3.38)
gegeben. Der Faktor fμ beschreibt den Einfluss einer Wand und berechnet sich aus
3, 45 y+
fμ = 1+ √ · tanh . (3.39)
Ret 70
Dabei ist l die mittlere freie Weglänge. Die verwendete Transportgleichung für die turbulente
kinetische Energie hat dieselbe Form wie beim k − Modell und wird durch die Gleichung
3.34 bestimmt. Die Transportgleichung für die turbulente Zeitskala lautet
∂(ρτ̃ ) ∂(ρτ̃ ) τ̃ ∂ ũi ∂ νt ∂ τ̃
+ ũi = (C
1 − 1) τij + (C
2 f2 − 1) + ν+
∂t ∂xi k̃ ∂xj ∂xi στ 2 ∂xi
2 νt ∂ τ̃ ∂ τ̃ 2 νt ∂ k̃ ∂ τ̃
− ν+ + ν+ (3.41)
τ̃ στ 2 ∂xi ∂xi k̃ στ 1 ∂xi ∂xi
Bei dem Faktor f2 handelt es sich um einen Dämpfungsterm, der den Einfluss einer Wand
modelliert und folgendermaßen bestimmt wird
2
2 Re2 y+
f2 = 1 − exp − t 1 − exp − . (3.42)
9 36 5
Die Werte der empirischen Modellkonstanten finden sich in Tabelle 3.2. Werden die Werte
der Modellkonstanten στ 1 , στ 2 und σk zu Null gesetzt, so reduziert sich das k − τ Modell
auf das k − Modell wobei dann direkt in Gleichung 3.41 die turbulente Zeitskala durch
τ = k
ersetzt werden kann.
3.3 STATISTISCHE MODELLIERUNG DER TURBULENZ 27
Cμ 0,09
C
1 1,44
C
2 1,83
στ 1 1,36
στ 2 1,36
σk 1,36
Reynoldsspannungsmodelle
2 ∂u
Rij = δij p k (3.47)
3 ∂xk
Die Gleichungen zeigen auf, dass bei RSM zum einen stets ein Modell für die Dreierkorre-
lationen der Geschwindigkeitsfluktuationen im Term Tkij und im Term Rij für die Kreuz-
korrelation aus Druckfluktuationen und des Gradienten der Geschwindigkeitsfluktuationen
gefunden werden muss. Detaillierte Betrachtungen zu Modellierungsansätzen dieser Terme
finden sich verschiedentlich in der Literatur (z.B. [129]).
Large-Eddy Simulation
Numerisch deutlich aufwendiger als der Einsatz von Reynolds-Spannungs-Modellen ist die
Anwendung der Large-Eddy Simulation (LES) [46, 129]. Sie kann sowohl vom Rechenauf-
wand her als auch vom Modellierungsgrad zwischen RSM und der direkten numerischen
Simulation (DNS) angesiedelt werden. Die der LES zu Grunde liegende Idee ist die einer
turbulenten Wirbelkaskade [146], die schematisch in Abbildung 3.2 dargestellt ist.
Die in einer turbulenten Strömung existierenden Wirbel teilen sich in unterschiedliche Grö-
ßen ein. Es gibt große energietragende Wirbelstrukturen, die durch die speziellen Strö-
mungsverhältnisse des betrachteten Problems bestimmt werden und den Bereich der kleinen
dissipierenden Wirbel, der für alle Strömungsverhältnisse selbstähnlich ist. Die Wirbelkas-
kade stellt den Ablauf des kontinuierlichen Zerfalls großer Wirbel in mehrere kleinskalige
Wirbel dar. Die großskaligen Wirbel können dabei durch unterschiedliche Prozesse entste-
hen. So können viskose Scherkräfte zwischen zwei Fluiden mit verschiedenen Geschwindig-
keiten, ein Querschnittssprung in einer Rohrströmung oder der Nachlaufeffekt eines geome-
trischen Hindernisses (z.B. Staukörper) in einer Strömung zur Bildung großskaliger Wirbel
Abbildung 3.2: Spektrum der turbulenten kinetischen Energie (Wirbelkaskade) [146, 159]
3.4 SCHLIEßUNGSPROBLEM DES CHEMISCHEN QUELLTERMS 29
führen. Die Größe der gebildeten Wirbel korreliert dabei mit den Längenskalen der sie
erzeugenden makroskopischen Strukturen. Eine Betrachtung der einzelnen Längenskalen
zeigt, daß hin zu kleinen Skalen die Wirbelstruktur mehr und mehr geometrisch unabhänig
wird. Es wird gemeinhin angenommen, dass ab einer bestimmten kleinen Skala die Turbu-
lenz als isotrop angesehen werden kann. Dies deutet auf die Existenz einer mittleren Skala
hin unterhalb der die Strömung als isotrop angesehen werden kann, und oberhalb der die
Strömung geometrieabhänig und somit auch problemabhänig ist. Dies ist die wesentliche
der LES zu Grunde liegende physikalische Motivation.
Im Unterschied zu den Mittelungsmethoden (RANS, URANS) werden bei der LES die hy-
drodynamischen Größen nicht in Mittelwert und Schwankungsgröße aufgeteilt, sondern es
wird eine Filterfunktion definiert, durch die die hydrodynamischen Größen in einen auflös-
baren Anteil (Grobstruktur, grid scale) und einen nicht auflösbaren Anteil (Feinstruktur,
sub-grid scale) aufgeteilt wird. Geht die Filterweite gegen Null (Δ → 0) dann geht die LES
in eine DNS über. Als Vorteil gegenüber RSM werden in der LES Instabilitäten auf großen
Skalen nicht gedämpft, sondern direkt repräsentiert. Bei Umströmungen von Staukörpern
oder Strömungen mit Verbrennung können diese Effekte signifikant sein. Die Grobstruk-
tur beinhaltet inhomogene und anisotrope Turbulenzstrukturen sowie energiereiche Wirbel.
Diese können direkt berechnet werden. Kleinere Strukturen, wie zum Beispiel energiearme
Wirbel, können nicht direkt aufgelöst werden. Modelliert werden muss somit der Einfluss
der Feinstruktur auf die Grobstruktur. Die Filterweite Δ sollte sorgsam gewählt werden da-
mit die gewünschte Größe der Struktur aufgelöst werden kann. Mit keinerem Δ wird aber
auch das Rechengitter engmaschiger und die Rechenzeit länger. LES Methoden unterschei-
den sich zum einen in der Art der verwendeten Filterfunktion sowie in den Modellen zur
Beschreibung der Interaktion großer und kleiner Skalen. Eine Übersicht über verschiedene
Modellansätze findet sich unter anderem in [129].
beschreiben. Sie beruhen im wesentlichen auf der stark vereinfachenden Annahme, dass die
Reaktionsrate durch die turbulente Dissipation bestimmt wird („mixed is burnt“). Analog
zum Abfall der Turbulenz wird eine Rate postuliert, die den Zerfall von Bereichen unver-
branntem Gases in Bereiche verbranntem Gases beschreibt. Diese Bruchstücke haben dann
ausreichen Kontakt mit bereits verbranntem Fluid und eine ausreichend hohe Temperatur,
so daß sie reagieren können. Für die Reaktionsrate eines beliebigen Brennstoffs F ergibt
sich somit
M w˙F = −ρ CF wF2 . (3.48)
k
Hierbei ist CF eine empirische Konstante der Größenordnung 1. Problematisch wird dieses
einfache Modell unter anderem in der Nähe von Wänden, wo sowohl k̃ als auch ˜ gegen Null
gehen und ihr Quotient somit nicht mehr einfach bestimmt ist. Deshalb sind häufig noch
zusätzliche Modellannahmen notwendig um zu Verhindern, dass der chemische Quellterm
nahe an einer Wand physikalisch nicht sinnvolle Werte annimmt [97].
zurückgegriffen werden. Die Gleichungen sollen für eine beliebige Größe φ formuliert sein
und es wird angenommen, dass sowohl das Geschwindigkeitsfeld als auch alle Eigenschaften
des Fluids bekannt sind.
Eine Disketisierung mittels Finite Volumen ist vorteilhaft, da die Navier-Stokes Gleichungen
integralen Charakter haben [42]. Das Lösungsgebiet wird durch die Diskretisierung in eine
endliche Zahl von Kontrolvolumen aufgeteilt. Der Einfachheit halber werden hier zur Er-
klärung zweidimensionale kartesische Gitter verwendet. Die Methodik läßt sich jedoch auch
für beliebige Gittertypen adaptieren. Für gewöhnlich werden die Kontrollvolumina durch
3.5 NUMERISCHE LÖSUNGSVERFAHREN FÜR DIE NAVIER-STOKES GLEICHUNGEN 31
b b b b
b b b b
b b b b
b b b b
b b b b
ein geeignetes Gitter festgelegt und die Zellmittelpunkte als Stützpunkte bei der Lösung der
Gleichungen verwendet (Abbildung 3.3). Die integrale Form der Erhaltungsgleichung muss
sowohl für jedes Kontrollvolumen als auch für das gesamte Lösungsgebiet erfüllt sein. Um
algebraische Ausdrücke für jedes Volumenelement zu erhalten müssen nun die Oberflächen-
und Volumenintegrale numerisch approximiert werden.
wobei f für die Normalkomponente des konvektiven (ρφu · n) oder des diffusiven (Γgradφ ·
n) Flusses steht. Die notwendigen Ableitungen können zum Beispiel über einen zentralen
Differenzenoperator bestimmt werden.
∂f fW − fE
= + O Δx2 (3.51)
∂x P 2 Δx
Als einfachste Approximation der Integrale wird dann beispielsweise die Mittelpunktsregel
verwendet. Exemplarische dargestellt für die Seite e lautet sie
Fe = Se f dS = fe Se ≈ fe Se . (3.52)
Se steht hierbei für die Seitenfläche des Kontrollvolumens, fe errechnet sich aus zwei be-
nachbarten Knotenpunkten beispielsweise nach
fE + fP
fe = . (3.53)
2
32 KAPITEL 3: NUMERISCHE SIMULATION REAKTIVER STRÖMUNGEN
Δx -
Nr
6
n 6
W
r w Pr e - Er Δy
s ?
?
y Sr
6
-x
Durch die Approximation der Volumen- und Oberflächenintegrale können die Erhaltungs-
gleichungen numerisch gelöst werden.
3.6 Wahrscheinlichkeitsdichtefunktion
Ansätze zur Beschreibung chemischer Reaktionen in einer turbulenten Strömung, die Fluk-
tuationen der Konzentration berücksichtigen, lassen sich anhand der zugrunde liegenden
Annahme über die Verteilungsfunktion der Fluktuationen der beteiligten Spezies charak-
terisieren. So beruht zum Beispiel, dass in Abschnitt 3.4 diskutierte Eddy Breakup Modell
für Vormischflammen auf der Annahme, dass die Reaktionen so schnell ablaufen, dass nur
verbranntes oder unverbranntes Gemisch vorliegen kann [140]. Die dazu gehörige Vertei-
lungsfunktion besteht folglich aus zwei δ-Funktionen. Werden hingegen Wahrscheinlich-
keitsdichtefunktionen (PDF) verwendet, so wird im Gegensatz dazu direkt versucht die
PDF zu berechnen. Klassifiziert werden können PDF Methoden grob in zwei Gruppen:
solche die a priori eine bestimmte Form der PDF annehmen (z.B. β-Funktion) und dann
nur Transportgleichungen für einige wenige Momente der Verteilungsfunktion lösen und
solche die die PDF direkt aus ihrer Transportgleichung bestimmen. Der zweite Ansatz ist
3.6 WAHRSCHEINLICHKEITSDICHTEFUNKTION 33
deutlich allgemeingültiger, da die Form der PDF im allgemeinen bei einem Verbrennungs-
prozess nicht als bekannt angesehen werden kann und zusätzlich sich bei der Annahme
einer bestimmten Form für die PDF das Problem ergibt, dass zur Darstellung von Ver-
bundwahrscheinlichkeitsdichtefunktionen Produktansätze verwendet werden müssen. Dies
ist bei statistisch nicht unabhänigen Größen wie zum Beispiel verschiedenen Spezieskon-
zentrationen eigentlich nicht zulässig. Aus diesen Gründen wird in dieser Arbeit nur der
zweite Ansatz verwendet. Er soll im folgenden näher erklärt werden.
∂F (V, t)
f (V ; t) ≡ (3.56)
∂V
vollständig beschrieben werden, wobei die Verteilungsfunktion die Wahrscheinlichkeit be-
schreibt mit der die Zufallsvariable U zu einem bestimmten Zeitpunkt kleiner als die Pro-
benvariable V ist.
In einer turbulenten Strömung ist der Geschwindigkeitsvektor U (x, t) ein zeitabhäniges
Zufallsvektorfeld. Die dazugehörige Verteilungsfunktion an einem Punkt in Ort und Zeit
lautet dementsprechend
F (V ; x, t) = Prob [Ui (x, t) < Vi ] (3.57)
mit i = 1, 2, 3. Ihre gemeinsame Wahrscheinlichkeitsdichtefunktion (joint PDF) ergibt sich
zu
∂ 3 F (V ; x, t)
f (V ; x, t) = . (3.58)
∂V1 ∂V2 ∂V3
An jedem Punkt in Ort und Zeit beschreibt diese PDF den Zufallsgeschwindigkeitsvektor
vollständig, aber beinhaltet ebenfalls keine gemeinsame Information über zwei oder mehr
34 KAPITEL 3: NUMERISCHE SIMULATION REAKTIVER STRÖMUNGEN
Wobei die zweite Zeile eine im folgenden verwendete Kurzschreibweise des Terms
∞ ∞ ∞
(· · · ) dV1 dV2 dV3 (3.61)
−∞ −∞ −∞
darstellen soll.
Die PDF besitzt zwei wichtige Eigenschaften. Zum einen ist sie immer positiv
und zum anderen ergibt eine Integration über der gesamten Zustandsraum, auf dem sie
definiert ist, die Normierungseigenschaft
f (V ; x, t) dV = 1 . (3.63)
3.6.2 Verbundwahrscheinlichkeitsdichtefunktion
Zur Beschreibung einer turbulenten Strömung mit überlagerter chemischer Reaktion, ist
als Erweiterung des bisher dargestellten, die Betrachtung der Verbundwahrscheinlichkeits-
dichtefunktion der Geschwindigkeit und des Zustandsvektors notwendig. Die Verbundwahr-
scheinlichkeitsverteilungsfunktion beschreibt die Wahrscheinlichkeit mit der an einer Stelle
in Ort und Zeit der Geschwindigkeitsvektor kleiner als der Probenvektor und der Zustands-
vektor ebenfalls kleiner als der dazugehörige Probenvektor ist.
Die JPDF ergibt somit die Wahrscheinlichkeit mit welcher der Zustandsvektor im Interval
ψ ≤ φ < ψ + dψ und der Geschwindigkeitsvektor im Interval V ≤ U < V + dV ist.
Zur später beschriebenen Lösung der PDF Transportgleichung werden bedingte Wahr-
scheinlichkeiten beziehungsweise bedingte Erwartungswerte benötigt. Die bedingte Wahr-
scheinlichkeit beschreibt wie groß die Wahrscheinlichkeit des Eintritts eines Ereignisses B
ist, sofern zuvor das Ereignis A eingetreten ist.
Prob(AB)
Prob(B | A) = (3.68)
Prob(B)
Dargestellt wird sie durch einen senkrechten Strich. Das hinter dem Strich stehende Ereignis
ist das vorausgesetzte. Prob(AB) ist die Wahrscheinlichkeit des gemeinsamen Auftretens
von A und B.
Analog lassen sich bedingte Wahrscheinlichkeitsdichtefunktionen definieren, z.B. die Wahr-
scheinlichkeit einer bestimmten Zusammensetzung unter der Voraussetzung, daß eine be-
stimmte Geschwindigkeit vorliegt.
fU φ (ψ, V ; x, t)
fφ (ψ | V = U ; x, t) = (3.69)
fU (V ; x, t)
Dies stellt exemplarisch den Erwartungswert der Größe Q für eine bestimmte Geschwindig-
keit V = U dar. Hierdurch kann wie später gezeigt wird der mittlere chemische Quellterm
unter der Annahme, daß seine PDF bekannt ist, exakt bestimmt werden. Die stark nicht-
lineare Abhänigkeit des chemischen Quellterms von der Temperatur kann somit detailliert
berücksichtigt werden.
werden. Eine detaillierte Herleitung findet sich unter anderem in [124,129] oder im Anhang
dieser Arbeit.
∂f ∂f ∂
p ∂f
ρ(ψ) + ρ(ψ)Vj + ρ(ψ)gj − +
∂t
∂xj
∂xj
∂Vj
1 2 3
∂
[ρ(ψ)Sα (ψ)fU φ (V, φ; x, t)] =
∂ψα
4
∂ ∂τij ∂p
− + fU φ (V, φ; x, t) +
∂Vj ∂xi ∂xi
5
∂ ∂Jiα
− | V, ψ fU φ (V, φ; x, t) (3.71)
∂ψα ∂xi
2. Änderung der PDF durch konvektiven Transport aufgrund des stochastischen Ge-
schwindigkeitsfeldes
4. chemischer Quellterm
Alle Terme auf der linken Seite der Gleichung, einschließlich des chemischen Quellterms,
erscheinen in geschlossener Form und bedürfen keiner weiteren Modellierung. Lediglich die
beiden Terme auf der rechten Seite der Gleichung müssen modelliert werden. Jedoch ge-
stalltet sich die Modellierung dieser beiden Terme, wie Abschnitt 3.6.5 ausgeführt, deutlich
einfacher als das Modellieren des (ungeschlossenen) mittleren chemischen Quellterms, wie
es bei vielen anderen Methoden zur Simulation turbulenter reaktiver Strömungen nötig ist.
Allerdings handelt es sich bei der PDF Transportgleichung um eine multidimensionale
Transportgleichung. Die Zahl der unabhänigen Variablen beträgt im dreidimensionalen Fall
7 (jeweils drei Dimensionen im Ortsraum und im Geschwindigkeitsraum, sowie die Zeit) plus
die Zahl der zur Beschreibung des thermochemischen Zustands notwendigen Variablen nS .
Wegen ihrer hohen Dimensionalität wird die PDF Transportgleichung im allgemeinen nicht
mit einem Finite Volumen oder Finite Differenzen Verfahren gelöst, da der Rechenauf-
wand zu groß wäre. Generell steigt bei solchen Verfahren der Rechenaufwand exponentiell
3.6 WAHRSCHEINLICHKEITSDICHTEFUNKTION 37
mit der Dimension der betrachteten Gleichung an. Deshalb werden zur Lösung der PDF
Transportgleichung Monte Carlo Verfahren eingesetzt. Diese finden weite Anwendung in
verschiedenen Bereichen der rechnergestützten Lösung physikalischer Probleme und haben
den wesentlichen Vorteil, dass bei ihnen der Rechenaufwand nur linear mit der Zahl der
betrachteten Dimensionen ansteigt. Die Details des auf numerischen Partikeln basierenden
Lösungsverfahrens finden sich im folgenden Abschnitt.
Zur Erklärung sei zunächst eine Massendichtefunktion (mass density function, MDF) ein-
geführt. Das tatsächlich verwendetet Lösungsverfahren für die PDF Transportgleichung
basiert auf dieser MDF.
Sei
x; t) ≡ ρf (V , ψ;
F(V , ψ, x, t) (3.73)
die erwähnte MDF, dann sind ihre beiden wichtigsten Eigenschaften
x; t) dV dψ
F(V , ψ, =
ρ (3.74)
x; t) dV dψ
F(V , ψ, dx = M (3.75)
wobei M die Gesamtmasse des Systems darstellt. Der massengewichtete Mittelwert einer
ergibt sich zu:
beliebigen Funktion Q von V und ψ
ρ Q̃ = F(V , ψ,
Q(V , ψ) x; t) dV dψ
(3.76)
nicht erfüllt. Die Größe F/M ist jedoch eine PDF, deren physikalische Interpretation wie
folgt ist:
38 KAPITEL 3: NUMERISCHE SIMULATION REAKTIVER STRÖMUNGEN
• Das gesamte betrachtete System wird auf einen Punkt mit der Gesamtmasse M re-
duziert, dessen Ort und Zustand durch eine vektorielle Zufallsvariable beschrieben
werden kann.
dx die Wahrscheinlichkeit, dass sich der Orts-, Geschwindigkeits-
• Dann ist F dV dψ
und Zustandsvektor jeweils im Interval [x, x + dx], V , V + dV bzw. ψ, ψ
+ dψ
befinden.
Zur Ableitung des formalen Zusammenhangs zwischen MDF und Zufallsgröße des Monte-
Carlo Verfahrens soll lediglich eine Realisierung der Strömung betrachtet werden. Das ge-
samte System reduziert sich somit auf ein sogenanntes stochastisches Partikel dessen mo-
mentaner Zustand als
F∗ (V , ψ, ∗ − V ) δ(φ
x; t) = M δ(U ∗ − ψ)
δ(x∗ − x) (3.78)
dargestellt werden kann. Das Partikel wird einem stochastischen Prozess unterworfen. Der
Erwartungswert des Partikelzustands ist dann die PDF F/M . Für die MDF kann analog
zu Gleichung 3.71 eine Transportgleichung formuliert werden. Sie findet sich unter anderem
in [124, 162]. Hierbei muss um vom stochastischen Prozess zur MDF zu kommen stets die
mathematische Operation der Erwartungswertbildung vorgenommen werden. In der kon-
kreten Umsetzung wird die Erwartungswertbildung durch Bildung des Mittelwertes ersetzt.
Um dabei zu belastbaren Ergebnissen mit einem möglichst geringen stochatischen Fehler
zu kommen ist eine gewisse Mindestanzahl an stochastischen Elementen bzw. stochatischen
Partikeln notwendig [165]. Das verwendete Ensemble nP der stochatischen Partikel, ergibt
eine Approximation FN der MDF F.
nP
FN = mP ∗(i) − V ) δ(φ
δ(U ∗(i) − ψ)
δ(x∗(i) − x) (3.79)
n=1
Wobei
M
mP = (3.80)
nP
die durch ein Partikel repräsentierte Masse ist. Da gilt
F =
FN (3.81)
handelt es sich bei FN um eine diskrete Approximation der kontinuierlichen MDF F durch
eine Summe von δ-Funktionen, die als ein Ensemble von stochatischen Partikeln inter-
pretiert werden können. Womit die Gültigkeit von Gleichung 3.72 gezeigt wurde und ihre
Bedeutung als diskete Approximation der PDF, die Abbildung 3.5 darstellt, einsichtig wird.
Gemäß den Anfangsbedingungen, die durch die PDF zum Zeitpunkt t = 0 gegeben sind,
wird im physikalischen Raum und im Zustandsraum ein Ensemble stochastischer Partikel
generiert, welches den Gleichungen 3.74 und 3.76 genügt. Die Partikel werden dann sto-
chastischen Prozessen unterworfen. Die kontinuierliche PDF läßt sich im einfachsten Fall
3.6 WAHRSCHEINLICHKEITSDICHTEFUNKTION 39
Abbildung 3.5: Wahrscheinlichkeitsdichtefunktion als Summe von δ-Funktionen (links) und die
dazugehörige Verteilungsfunktion (rechts) [9]
Abbildung 3.6: Schematische Darstellung eines allgemeinen Rechnengitters und der stochasti-
schen Partikel
lichsten Modellierungsansätze, die in der Literatur zu finden sind, sollen im folgenden dis-
kutiert werden.
Molekulares Mischen
Der Term
∂ ∂J α
− i | V, ψ fU φ (V, φ; x, t) (3.82)
∂ψα ∂xi
beschreibt den Einfluss von molekularem Mischen (Diffusion und dadurch eventuell be-
dingter Wärmeleitung) auf die Form der PDF und kommt in ungeschlossener Form vor, da
die PDF nur statistische Informationen über einen Punkt, nicht aber über die bedingten
Wahrscheinlichkeiten der lokalen Gradienten enthält.
Es existieren eine Vielzahl verschiedener Modellierungsansätze in der Literatur. Eine ver-
gleichende Studie verschiedener Modelle für nicht-vorgemischte Flammen aus jüngster Zeit
findet sich unter anderem in [99,144]. Aus Modellierungssicht ist es sehr entscheidend für das
molekulare Mischen (häufig auch Mikromischen genannt) allgemeingültige und möglichst
exakte Ansätze zu finden. Die zu modelliernden Prozesse spielen sich auf den kleinsten Ska-
len ab und verursachen enorme Schwierigkeiten bei der Beschreibung [126]. Von Fox [43],
3.6 WAHRSCHEINLICHKEITSDICHTEFUNKTION 41
Pope [129] und Subramaniam und Pope [145] werden folgende Kriterien genannt, die ein
optimales Mischungsmodell erfüllen sollte.
(iii) In homogener Turbulenz soll die gemeinsame PDF inerter Skalare auf eine Gaussfunk-
tion relaxieren
(iv) Die Realisierbarkeit muss Berücksichtigung finden, so sind z.B. die Molenbrüche auf
das Interval [0; 1] beschränkt
(vii) Die reale Abhängigkeit von skalaren Längenskalen muss dargestellt werden
(viii) Eine korrekte Abhänigkeit von Reynolds-, Schmidt- und Damköhler-Zahl ist notwen-
dig
Die ersten Modellierungsansätze für das molekulare Mischen gehen auf die Arbeiten von [32]
zurück. Das dort vorgestellte Coalescence-Dispersion (CD) Modell soll hier nicht näher
dargestellt werden. Detailliert vorgestellt werden im folgenden das IEM (interaction by
exchange with mean) Modell, das Curl Modell sowie das in dieser Arbeit verwendete mo-
difizierte Curl Modell von Janicka et al. [61]. Abschließend werden dann noch einige Be-
merkungen zu aktuell entwickelten Mischungsmodellen gemacht. Die Erläuterungen sollen
sich dabei auf die Betrachtungsebene der stochatischen Prozesse (stochatischen Partikeln)
beschränken.
IEM Modell
Ausgangspunkt des IEM Modells [34, 35, 155] ist die Beobachtung, dass die Zusammen-
setzungsfluktuationen immer in Richtung des Mittelwertes relaxieren. Die Geschwindigkeit
dieses Prozesses wird durch den empirisch bekannten Varianzabbau festgelegt. Auf Parti-
kelebene bedeutet dies
11 ∗
Δφ∗α = − (φ −
φα ) Δt (3.83)
2 τφ α
mit
1
τφ = τt . (3.84)
Cφ
42 KAPITEL 3: NUMERISCHE SIMULATION REAKTIVER STRÖMUNGEN
φ2
6
u
Z
Z
Z
Z u
Z
~
Z
)
φ2 u
Z
}
Z
Z
Z
Z
Zu
u
-
φ1 φ1
Abbildung 3.7: Turbulente Vermischung eines Partikelensembles nach dem IEM Modell
φα ist dabei der Massenbruch der Spezies α, Cφ eine Modellkonstante mit der die Mischungs-
geschwindigkeit des Modells an empirische Daten angepasst wird und τt die turbulente
Zeitskala. Abbildung 3.7 stellt für ein System mit zwei Skalare φ1 und φ2 das Verhalten
des IEM Modells dar. Jeder Punkt stellt ein stochatisches Partikel dar. Während eines
Zeitschritts Δt relaxieren alle Partikel mit einer individuell verschiedenen Geschwindigkeit
hin zu dem gemeinsamen Erwartungswert des Partikelensembles. Die Änderung der Par-
tikelzusammensetzung durch die modellierte Relaxation sollen die Pfeile in der Abbildung
darstellen.
Dieses Modell erfüllt die Anforderungen (i), (ii), (iv) und (v), jedoch kann es nicht die
Relaxation zu einer Gaussfunktion darstellen (iii) und ebenfalls keine Lokalität des Mi-
schungsprozesses im Zustandraum gewährleisten (vi). Der wesentliche Nachteil des Modells
ist jedoch, dass es die Form der PDF erhält, was eine falsche Beschreibung der realen
physikalischen Bedingungen darstellt. Dies kann bei homogener Turbulenz zu sehr schlech-
ten Ergebnissen führen, wenn keine richtige Anfangsbedingung für die PDF vorgegeben
wird. Jedoch bleibt festzuhalten, dass dieser Nachteil bei inhomogener Turbulenz deutlich
weniger ins Gewicht fällt, da hier durch Markomischungsvorgänge ebenfalls die Form der
lokalen PDF der Skalare beeinflusst wird. Deshalb und wegen seiner Einfachheit und Effizi-
enz ist das IEM Modell immer noch das am meisten verwendete Mischungsmodell [105,106].
Curl Modell
Beim Curl Modell handelt es sich um ein stochastisches Partikelmodell. In seiner Grundform
wurde es von Curl veröffentlicht [32]. Weitere Ausführungen dazu finden sich auch in [123].
3.6 WAHRSCHEINLICHKEITSDICHTEFUNKTION 43
Der Algorithmus läuft wie folgt ab: Für einen Zeitschritt Δt wird aus der Gesamtzahl NP
des Partikelensembles in einer Zelle des Rechengebiets eine Anzahl zu mischender Partikel
Nmix gemäß
1
Nmix = β NP Δt (3.85)
τφ
mit
1
τφ = τt (3.86)
Cφ
ausgewählt. Aus den so gewählten Partikeln werden nun zufällig Paare Q und P gebildet
und miteinander vermischt.
1
φneu
Q = (φQ + φP ) (3.87)
2
1
φneu
P = (φQ + φP ) (3.88)
2
Die Konstante β wird dabei so festgelegt, dass die Varianzabnahme mit der Zeit korrekt
wiedergegeben werden kann. Allerdings hat das Modell in dieser Form den wesentlichen
Nachteil, dass es nur diskrete Mischungszustände wiedergeben kann. Leicht einsichtig wird
dies am Beispiel der Vermischung eines passiven Skalars der als Anfangsverteilung zwei
δ-Funktionen bei φ = 0 und φ = 1 besitzt. Beim Mischen entstehen nun neue diskrete
Zustände. Nach dem ersten Mischungsschritt existieren somit die Zustände φ = 0, φ = 12
und φ = 1. Beim zweiten Mischungsschritt kommen die Zustands φ = 14 und φ = 34 hinzu.
Setzt man diese Reihe beliebig lange fort, so ist leicht einsichtig, dass alle dabei generierte
Zustände Vielfache von n−2 , mit n als Anzahl der Mischungsschritte, sind. Offensichtlich
entspricht dies nicht der Realität, da ein Skalar jeden Wert beliebigen zwischen 0 und 1
einnehmen kann.
Dieser Nachteil wird durch eine in [61] dargestellte Modellerweiterung behoben. Dabei wer-
den die Partikelpaare nicht mehr in einem festen Verhältnis miteinander vermischt, sondern
mit einem Mischungsgrad α.
1
φneu
Q = (1 − α) φQ + (φQ + φP ) (3.89)
2
1
φneu
P = (1 − α) φP + (φQ + φP ) (3.90)
2
Der Mischungsgrad α ist eine auf dem Interval [0..1] gleichverteilte Zufallsvariable. Die An-
zahl der zu mischenden Partikel (Gleichung 3.85) sowie die turbulente Mischungszeitskala
(Gleichung 3.86) bleiben dabei jeweils unverändert zum Curl Modell. Lediglich die Modell-
konstante β muss an diese Modifikation angepasst werden und hat in der Literatur den
Wert 3 [61]. Abbildung 3.8 stellt einen mit dem modifizierten Curl Modell beschriebenen
44 KAPITEL 3: NUMERISCHE SIMULATION REAKTIVER STRÖMUNGEN
φ2
6 u u
u
u Q
uφ
u
u neu u
φQ
neu u
φ P
u
P u
φ
u
-
φ1
Abbildung 3.8: Abbildung der Mischung in einem vereinfachten Zustandsraum mit dem modi-
fizierten Curl Modell
Weitere Modellierungsansätze
Ohne Anspruch auf Vollständigkeit sollen im folgenden kurz einige weitere Modellierungs-
ansätze für molekulare Mischungsprozesse dargestellt werden.
3.6 WAHRSCHEINLICHKEITSDICHTEFUNKTION 45
mit
⎡ ⎛
⎞⎤
11 ⎣ φ2
Aφ = 1 + K0 ⎝1 − 2 ⎠⎦ (3.92)
2 τφ φmax
1 (φ∗ −
φα )2 2
Bφ = 1− α 2 φ (3.93)
τφ φmax
Durch den Term Aφ wird die Relaxation in Richtung des Mittelwertes beschrieben,
Term Bφ stellt einen stochastischen Diffusionsterm dar. Die Modellkonstante K0 hat
einen Wert von 0, 7 und der Term dW b ist ein binomiales Wienerinkrement, welches
√
durch den Term ξ dt modelliert wird. ξ ist dabei eine binomialverteilte Zufallsvaria-
ble. Das Modell erfüllt neben den Bedingungen (i), (ii), (iv) und (v) noch Bedingung
(iii), die eine Relaxation der Verteilung hin zu einer Gaussverteilung in homogener
Turbulenz fordert. Allerdings erfüllt es ebenfalls nicht die Anforderung der Lokali-
tät im Zustandsraum. Problematisch ist ebenfalls die relative hohe Rechenzeit des
Modells bedingt durch die aufwändige Generierung binomialverteilter Zufallszahlen
sowie seine schwere Erweiterbarkeit auf Systeme mit mehreren Skalaren. Eine Stu-
die zur Leistungsfähigkeit eines solchen Mischungsmodells im Vergleich zu anderen
Mischungsmodellen findet sich unter anderem in [112].
• EMST Modell
Ein Modell welches zusätzlich noch die Lokalitätsbedingung im Zustandsraum erfüllt
ist das von Subramanian und Pope [145] vorgestellte EMST Modell. Es basiert auf
der mathematischen Konstruktion von „Euclidian minium spanning trees“ im Zu-
standsraum. Durch einen EMST ist die unmittelbare Umgebung eines Partikels im
Zustandsraum bestimmt. Die Mischung eines Partikels findet dann längs einer durch
den lokalen EMST des Partikels festgelegten Weges statt.
• PSP Modell
Das von Meyer und Jenny [100] vorgestellte PSP Modell arbeitet mit parametrisierten
Skalarprofilen (PSP). Die grundlegende Idee hinter diesem Modell ist die Parame-
trisierung eines eindimensionalen Skalarprofils. Die Mischung zweier Partikel erfolgt
46 KAPITEL 3: NUMERISCHE SIMULATION REAKTIVER STRÖMUNGEN
dann längs eines so konstruierten Skalarprofils. Weitere Details und eine genauere
Beschreibung findet sich in der angegebenen Literaturstelle.
u2
6
u
@
@
@
@
@
@
Re
@
e
I
@
@
@
@
@
@
@u
-
u1
Abbildung 3.9: Effekt der Mischung zweier stochatischer Partikel im Geschwindigkeitsraum (hier
vereinfacht nur die u1 , u2 Ebene dargestellt)
Der isotrope Teil (− 23 δij ) der Gleichung sorgt für eine Abnahme der Normalspannungen
( u1 u1 , u2 u2 und u3 u3 ), wohingegen der anisotrope Teil Rjk dazu dient Energie zwi-
schen den einzelnen Reynoldsspannungen umzuverteilen. Die Spur von Rjk ist Null. Die
beobachtete Tendenz der Reynoldsspannungen isotrop zu werden, läßt sich somit durch
diesen Term modellieren. Rotta [134] schlägt hierfür die Gleichung
uj uk − 23 k δjk
Rjk = −C2 (3.99)
τ
vor. Hierin ist C2 eine empirische Modellkonstante. Ein gutes Modell für die turbulenten
Druckfluktuationen und viskosen Scherkräfte muss somit sicherstellen, daß die PDF der
Geschwindigkeit zu einer Gaussfunktion relaxiert, sich der Mittelwert der Verteilung nicht
ändert und die zweiten Momente der Verteilung sich (in etwa) wie in Gleichung 3.99 ver-
halten. Dies kann durch eine Kombination des stochastischen Mischungsmodells und eines
stochastischen Reorientierungsmodells erreicht werden.
Das stochastische Mischungsmodell funktioniert anlaog zum Mischungsmodell für die ska-
lare Dissipationsrate. Mit der Wahrscheinlichkeit Δt
τu
= Δt
τ
Cu N mischen zwei Partikel mit-
einander. Diese Partikelpaare werden zufällig (ohne zurücklegen) aus dem Partikelensemble
ausgewählt. Cu ist ein Modellparameter. Wie Abbildung 3.9 veranschaulicht wird die Ge-
schwindigkeit jeden Partikels durch den Mittelwert beider Partikel ersetzt.
U ∗(P ) (t + Δt) = 1 U
∗(Q) (t + Δt) = U ∗(Q) (t) + U
∗(P ) (t) (3.100)
2
48 KAPITEL 3: NUMERISCHE SIMULATION REAKTIVER STRÖMUNGEN
u2
6
u e
@
@
@
@
@
@
@
R
@
@
I
@
@
@
@
@
@
e @u
-
u1
Abbildung 3.10: Effekt der Mischung und stochatischen Reorientierung zweier stochatischer
Partikel im Geschwindigkeitsraum (hier vereinfacht nur die u1 , u2 Ebene dargestellt)
mit ΔU P Q (t) = U
∗(P ) (t)−U
∗(Q) (t) als dem Geschwindigkeitsunterschied der beiden Partikel.
ξ ist ein zufällig orientierterter Einheitsvektor. Dieser in Abbildung 3.10 veranschaulichte
Prozess ändert weder den Mittelwert noch den Geschwindigkeitsunterschied zweier Partikel
und beeinflusst somit weder den Impuls noch die turbulente kinetische Energie. In [124]
wird dargestellt, dass dieses Modell den gleichen Effekt auf die Reynoldsspannungen wie
das durch Gleichung 3.99 beschriebene Modell hat.
3.6 WAHRSCHEINLICHKEITSDICHTEFUNKTION 49
Ebenso kann der Einfluss von Druckfluktuationen und viskosen Scherkräften auf die PDF
auch durch einen einen Langevin Ansatz [29, 71] der Form
∂
∂ ∂τij ∂p 1 ∂ 2f
− + |V , φ f = −Gij Vj − Ũj f + B (3.104)
∂Vj ∂xi ∂xi ∂Vj 2 ∂Vj ∂Vi
modelliert werden. Im Unterschied zum stochatischen Mischungs- und Reorientierungsmo-
dell beruht dieser Modellansatz auf Partikelebene nicht auf der Interaktion von zwei zufällig
ausgewählten Partikeln, sondern leitet eine Gleichung für die zeitliche Entwicklung der Ge-
schwindigkeit eines Partikels ab. Die Bestimmung der Koeffizienten Gij und B in Gleichung
3.104 erfolgt aus grundlegenden Hypothesen der Turbulenztheorie [101].
Auch durch ein Langevinmodell soll der Abbau der Geschwindigkeitfluktuationen, was
gleichbedeutend mit dem Abklingen der Turbulenz und der Abnahme der Varianz der Ge-
schwindigkeitsverteilung ist, richtig beschrieben werden können. Die turbulente kinetische
Energie soll deshalb ebenfalls nach der Beziehung
∂k
= − (3.105)
∂t
abklingen. Unter Annahme isotroper Turbulenz läßt sich zeigen [124] das dies durch den
Term
1 3
Gij = − + C0 δij (3.106)
2 4 k
gewährleistet wird.
Der Koeffizient B wird so gewählt, dass die zeitliche Kovarianz zweier Fluidpartikel mit den
Vorhersagen der Turbulenztheorie übereinstimmt [101]. Für die zeitliche Kovarianz zweier
Fluidpartikel gilt demnach
Δδt Ui (t)Δδt Uj (t) = C0 δt δij (3.107)
(t) = U
Δδt U (t + δt) − U
(t) (3.108)
wobei C0 eine universelle Konstante mit dem Wert 2, 1 ist und δt ein Zeitintervall darstellt,
dass wesentlich größer als die Kolomogorovzeitskala τk und wesentlich kleiner als die tur-
bulente Zeitskala τt ist. Die Kolomogorovzeitskala kann als Umdrehungszeit der kleinsten
Wirbel einer turbulenten Strömung interpretiert werden und als turbulente Zeitskala soll
das integrale turbulente Zeitmaß der Strömung herangezogen werden.
τk δt τt (3.109)
Einsetzten des Langevin Ansatzes in die PDF Transportgleichung liefert für die zeitliche
Kovarianz die Beziehung
Δδt Ui (t)Δδt Uj (t) = B δt δij (3.110)
50 KAPITEL 3: NUMERISCHE SIMULATION REAKTIVER STRÖMUNGEN
Im Unterschied zum vereinfachten Langevin Modell (SLM) ist der Term Gij beim verallge-
meinerten Langevin Modell (GLM) kein Skalar sondern ein Tensor. Er ergibt sich als
α1 δij + α2 δij ∂
Uk
Gij = + Hijkl (3.113)
τ ∂xl
worin
Hijkl = β1 δij δkl + β2 δij δjl + β3 δil δjk
+ γ1 δij δkl + γ2 δij δjl + γ3 δil δjk
+ γ4 δij δkl + γ5 δij δjl + γ6 δil δjk (3.114)
die elf Koeffizienten αm , βm und γm müssen noch bestimmt werden. Die Herleitung und die
wesentlichen Überlegungen hierzu finden sich in [124]. Ihre Zahlenwerte sind zusätzlich im
Anhang B angegeben.
gestellt. Wie bereits erwähnt wird die Lösung der multidimensionalen PDF Transportglei-
chung ersetzt durch die Lösung einer Reihe von gewöhnlichen (stochatischen) Differential-
gleichungen für die stochatischen Monte Carlo Partikel. Die Gleichungen lassen sich als Be-
wegungsgleichungen jeweils für den physikalischen Raum, den Geschwindigkeitsraum und
den Zustandsraum interpretieren
• Bewegung im Geschwindigkeitsraum
Zur Beschreibung der Bewegung eines Partikels in Geschwindigkeitsraum wird ein
vereinfachtes Langevin Model (Simplified Langevin Model, SLM [124]) verwendet.
$ % dt
∂
p 1 3 C0 k
dUj∗ = − dt − + C0 Uj∗ −
Uj + dWj (3.116)
∂xj 2 4 τ τ
Die Gleichung ist hier allgemein für die Komponente Uj des Geschwindigkeitsvektors
in Richtung des Vektors xj formuliert. C0 ist eine Modellkonstante, die in SLM den
Wert 2, 1 hat. Bei τ handelt es sich um ein turbulentes Zeitmaß, k ist die turbulente ki-
netische Energie und Wj die j-Komponenten eines vektorwertigen Wienerinkrements.
• Bewegung im Zustandsraum
Die Bewegung im Zustandsraum wird durch die Summe zwei Terme beschrieben: M
und S.
dφ
= M + S (3.117)
dt
Dabei ist M der Einfluss von molekularer Mischung. Dieser muss modelliert werden.
Hierzu wird ein modifiziertes Curl-Modell verwendet [61]. Der chemische Quellterm
S tritt in geschlossener Form auf und bedarf keiner weiteren Modellierung.
3.7.1 Bruttoreaktionsmechanismen
Ein einfaches Beispiel für eine Bruttoreaktion, ist die Oxidation von Methan mit Sauerstoff
nach der Gleichung
CH4 + 2 O2 −→ CO2 + 2 H2 O
Diese Reaktion kann auf molekularer Ebene nicht ablaufen. Eine detaillierte Betrachtung
zeigt, dass hierfür bei einem Stoß der drei Moleküle eine Vielzahl von Bindungen gebro-
chen und neu gebildet werden müssten. Die Wahrscheinlichkeit hierfür ist verschwindent
gering. Tatsächlich werden eine ganze Reihe von Zwischenprodukten wie O, OH, CH und
CO gebildet. Zusätzlich ist die mathematische Beschreibung des Reaktionsgesetzes einer sol-
chen Reaktion meist sehr schwierig. Ihre Reaktionsordung ist in der Regel nicht ganzzahlig
und kann zusätzlich eine Funktion der Zeit und der Reaktionsbedingungen sein [87, 159].
Nichts desto trotz werden einfache Reaktionsmechanismen basierend auf Bruttoreaktion,
auch globalen Reaktionsstufen genannt, formuliert. Ein Beispiel sei hier ein Reaktionsme-
chanismus, der die Verbrennung von Methan mit Luft durch vier globale Reaktionsschritte
beschreibt [116, 119].
CH4 + 2 H + H2 O −→ CO + 4 H2
CO + H2 O −→ CO2 + H2
H+H+M −→ H2 + M
3.7 BESCHREIBUNG DER REAKTIONSKINETIK 53
O2 + 3 H2 −→ 2 H + 2 H2 O
Alles diese Reaktionen können ebenfalls nicht auf molekularer Ebene ablaufen. Aller-
dings kann hiermit schon die Entstehung zwischen Zwischenprodukten dargestellt werden,
die auch je nach Reaktionsbedingungen als Endprodukte der Verbrennung übrig bleiben
können.
OH + H2 −→ H2 O + H
(i) Unimolekulare Reaktionen beschreiben die Umlagerung oder Dissoziation eines Moleküls.
A −→ Produkte
Unimolekulare Reaktionen haben ein Zeitgesetz erster Ordung, dass heißt verdoppelt
sich die Anfangskonzentration dann verdoppelt sich ebenfalls die Reaktionsrate.
(ii) Bimolekulare Reaktionen sind die am häufigsten vorkommenden Reaktionen. Sie laufen
nach der Reaktionsgleichung
A + B −→ Produkte
ab und haben immer ein Zeitgesetz zweiter Ordnung, was bedeutet daß eine Verdopp-
lung der Anfangskonzentration jedes Reaktionspartners zu einer Vervierfachung der
Reaktionsrate führt.
54 KAPITEL 3: NUMERISCHE SIMULATION REAKTIVER STRÖMUNGEN
(iii) Trimolekulare Reaktionen folgen einem Zeitgesetz dritter Ordnung. Gewöhnlich han-
delt es sich bei ihnen um Rekombinationsreaktionen.
A + B + C −→ Produkte
gegeben ist, dann kann das Zeitgesetz für die Bildung der Spezies i in Reaktion r durch
den Ausdruck
S
∂ci (P ) (E) & νrs (E)
= kr νri − νri cs (3.119)
∂t chem,r s=1
(E) (P )
dargestellt werden. Dabei sind νrs und νrs die stöchiometrischen Koeffizienten der Edukte
und Produkte und cs bezeichnet die Konzentration der S verschiedene Spezies s.
Charakteristisch für chemische Reaktionen ist, dass ihre Geschwindigkeitskoeffizienten stark
nichtlinear von der Temperatur abhängen. Nach Arrhenius [3] kann diese Abhängigkeit mit
einer relativ einfachen exponnential Beziehung ausgedrückt werden.
E
k = A · exp − a (3.120)
RT
Neuere Arbeiten zeigen eine Temperaturabhänigkeit des präexponnentiellen Faktors A wes-
halb das Arrheniusgesetz häufig in modifizierter Form verwendet wird.
b Ea
k = A T · exp − (3.121)
RT
Hierin ist jeweils Ea die Aktivierungsenergie der Reaktion, T die Temperatur und R die
universelle Gaskonstante. Die Aktivierungsenergie korrespondiert mit der beim Start der
Reaktion zu überwindenden Energiebariere. Ihr Zahlenwert hängt von der Bindungsenergie
der beteiligten Moleküle ab und kann je nach dem wieviele und welche Art von Bindungen
bei der Reaktion aufgebrochen und neu gebildet werden müssen deutlich unterschiedliche
Größen annehmen.
Zusätzlich zur Temperatur hängt die Reaktionsgeschwindigkeit noch vom Druck ab. Dies ist
leicht einsichtig es bei bei höherem Druck zu mehr Stößen pro Zeiteinheit der Moleküle un-
tereinander kommt. Die Berücksichtigung dieses Effekts auf die Reaktionsgeschwindigkeit
soll hier nicht näher beschrieben werden. Der Einfluss des Drucks auf die Reaktionsge-
schwindigkeit ist aber deutlich schwächer als die Temperaturabhängigkeit [5, 79].
3.7 BESCHREIBUNG DER REAKTIONSKINETIK 55
• auf einer von Hand durchgeführten detaillierten Analyse des chemischen Reaktions-
mechanismus beruhen, wie das chemical lumping [2, 115],
• auf einer detaillierten Simulation des vollständigen Zustandsraumes beruhen und eine
(problemangepasste) möglichest optimale Parameterisierung des Systems mit wenigen
Freiheitsgraden ansteben. Als Beispiel sei die FGM (Flamelet generated Manifold)
Methode [110, 111, 117] genannt,
• auf einer effizienten Tabellierungs- und Lookupstrategie für den detailierten berechne-
ten chemischen Quellterm beruhen und damit die notwendige Simulationszeit deutlich
reduzieren können, als Vertreter diese Gruppe kann das ISAT (In-situ adaptive tabu-
lation) Verfahren [81, 128] angeführt werden
• und Methoden, die auf einer Zeitskalenanalyse des detaillierten chemischen Reaktions-
systems basieren. Anzuführen ist hier die QSSA (Quasi-Stationarität) [118] Annahme,
das CSP (Computational singular pertubation) Verfahren [54], MIM (Method of in-
tegral/invariant Manifolds) Methode [49], ILDM Konzept [89, 90] und die in dieser
Arbeit verwendete REDIM Methode [22, 24].
Gewählt wurde letzteres Verfahren, da es eine sehr zuverlässige Beschreibung der Reakti-
onskinetik bereits mit sehr wenigen Parametern erlaubt und zusätzlich zu dem Einfluss von
chemischen Prozessen noch den Einfluss von molekularem Transport auf den Zustandsraum
berücksichtigt. Verwendet wurden lediglich zwei Parameter. Aber damit lassen sich sowohl
Majoritäten als auch Minoritätenspezies mit hoher Genauigkeit vorhersagen. Ebenso gleicht
sie einige Defizite des in vorherigenen Arbeiten [9, 104, 162] gemeinsam im Kontext eines
hybriden CFD/transported PDF Modells verwendeten ILDM Verfahrens aus.
56 KAPITEL 3: NUMERISCHE SIMULATION REAKTIVER STRÖMUNGEN
Das ILDM Verfahren liefert eine gute Approximation der Dynamik des detaillierten Sy-
stems für homogene Systeme. In diesem Grenzfall relaxieren die schnellen chemischen Pro-
zesse auf eine niedrigdimensionale Mannigfaltigkeit, die die Dynamik des langsamen Teils
des Systems beinhaltet und seine Wechselwirkung mit Konvektion und Diffusion. Dieser
Sachverhalt ist in Abbildung 3.11 dargestellt. Sie zeigt eine Darstellung der Dynamik eines
reaktiven Systems in einer Projektion auf einen zweidimensionalen Zustandsraum, beschrie-
ben durch den Massenbruch von CO2 und H2 O. Ausgehend von einem beliebigen, durch
die Anfangspunkte der roten Pfeile dargestellenten, Punkt bewegt sich das System durch
die schnellen Prozesse auf die grün dargestellte niedrigdimensionale Mannigfaltigkeit. Von
dort ausgehend bringt die Dynamik der langsamen Prozesse das System in das chemische
Gleichgewicht (blauer Punkt). Jedoch gibt es Probleme mit der Existenz der ILDM und
ihrer attraktiven Eigenschaften in der gesamten für die Anwendung in einer realen turbu-
lenten Flamme interessierenden Domäne des Zustandsraums. Hierfür gibt es verschiedene
Gründe. So wird zum einen im Niedertemperaturbereich der chemische Quellterm vernach-
lässigbar klein und zu anderen gibt es einen Übergangsbereich intermediärer Zeitskalen in
welchem chemische und physikalische Zeitskalen überlappen und das System weg von der
niedrigdimensionalen Mannigfaligkeit gestört wird. Anderes ausgedrückt kann eine Man-
nigfaltigkeit, die nur durch Betrachtung des Reaktionsterms abgeleitet wird, selbst wenn
sie überall definiert sein sollte, die Interaktion zwischen Transport- und Reaktionsprozessen
in allgemeiner Art und Weise nur unzureichend beschreiben. Teilweise lassen sich wie von
Bykov und Maas dargestellt [23, 25] diese Probleme durch eine geschickte Erweiterung der
ILDM in den restlichen Zustandsraum beheben. Hierbei wird der interessierende Zustands-
raum in drei Bereiche aufgeteilt. Im ersten Bereich ist die chemische Kinetik bestimmend
für das Verhalten des Systems, im zweiten Bereich sind Reaktion, Konvektion und Diffu-
3.7 BESCHREIBUNG DER REAKTIONSKINETIK 57
sion stark miteinander gekoppelt und im dritten Bereich (unendlich) langsamer Chemie
wird das System ausschließlich durch Konvektions- und Diffusionsprozesse bestimmt. Alle
drei Bereiche werden durch verschiedene niedrigdimensionale Mannigfaltigkeiten beschrie-
ben. Das wesentliche Manko dieser Betrachtungsweise ist, dass angenommen wird, dass die
mittlere Domäne asymtotisch in die Grenze zwischen erster und dritte Domäne schrumpft.
Was nur im Falle einer starken nichtlinearen Abhänigkeit des Quellterms von den System-
parametern gegeben ist. Diese Einschränkung entfällt durch den in folgenden beschriebenen
REDIM („Reaction-Diffusion Manifold“) Ansatz.
Die dem REDIM Formalismus zu Grunde liegende Idee, ist das Konzept der invarianten
Mannigfaltigkeiten im Lösungsraum einer allgemein formulierten Konvektions-Reaktions-
Diffusionsgleichung. In symbolischer Vektorschreibweise ist das typische System partieller
Differentialgleichungen (PDEs) zur Beschreibung einer reaktiven Strömung durch folgende
Gleichung gegeben
∂Ψ 1
= F (Ψ) − v · grad Ψ − div (D · grad Ψ) ≡ Φ(Ψ) (3.122)
∂t ρ
worin v das Geschwindigkeitsfeld, ρ die Dichte und D die (n × n)-dimensionale Matrix der
Transportkoeffizienten [48] darstellt. Der Zustandsvektor Ψ ist der (n = ns +2)-dimensionale
w1 ws
Vektor Ψ = (h, p, M 1
,··· , M s
) mit der Enthalpie h, dem Druck p, den ns Speziesmassenbrü-
chen w1 , · · · , ws und den molaren Massen M1 , · · · , Ms . F (Ψ) ist der n-dimensionale Vektor
des thermochemischen Quellterms und t bezeichnet die Zeit. Jede andere Formulierung des
Zustandsvektors wäre in diesem Zusammenhang ebenfalls anwendbar.
Angenommen die Lösung des Systems im Zustandsraum ist nahe oder Teil einer durch eine
explizite Funktion Ψ(θ) bestimmte ms -dimensionaler Mannigfaltigkeit. Diese Annahme ist
valide, da typische Verbrennungssysteme von stark verschiedene Zeitskalen geprägt sind,
was zu einer Zerlegung des Systems in langsame und schnelle Prozesse führt. Mathematisch
wird diese Mannigfaltigkeit beschrieben als
M = Ψ : Ψ = Ψ(θ), Ψ : ms → n . (3.123)
Ψ+Φ stellt die Moore-Penrose pseudoinverse Matrix von ΨΦ dar [50]. Berechnet wird sie für
eine reguläre Matrix Ψ+Φ · ΨΦ durch
−1
Φ = ΨΦ · ΨΦ
Ψ+ +
· ΨTΦ (3.126)
und existiert immer, wenn die Spalten einer Matrixdarstellung von ΨΦ linear unabhänig
sind. Diese Bedingung kann durch eine geeignete Wahl der lokalen Koordinaten stets er-
füllt werden. Von Bykov und Maas [24] wird eine Umformulierung von Gleichung 3.125
in ein System partieller Differentialgleichungen für Ψ = Ψ(Φ) zur Lösung der Gleichung
vorgeschlagen,
∂Ψ(θ)
= I − ΨΦ Ψ+ Φ · Φ (Ψ(θ)) (3.127)
∂t
so dass eine stationäre Lösungsmenge Ψ(θ, ∞) des obrigen partiellen Differentialgleichungs-
systems die gewünschte Mannigfaltigkeit ergibt. Hierzu wird das System (Gleichung 3.127)
ausgehend von einer geratenen Startlösung solange integriert bis sich eine stationäre Lösung
ergibt. Die Dynamik des Systems sich vollständig als eine Bewegung innerhalb der der aus
diese Weise bestimmten Mannigfaltigkeit darstellen.
Durch Umschreiben des Diffusionsterms und mit der Annahme gleicher Diffusivität für alle
Spezies (D = d · I) kann das Gleichungssystem zur Bestimmung der Mannigfaltigkeit in
vereinfachter Form geschrieben werden
∂Ψ(θ) 1
= I − ΨΦ (θ)Ψ+
Φ (θ) · F (Ψ) − d Ψθθ ◦ grad θ ◦ grad θ . (3.128)
∂t ρ
Die Gleichung verdeutlich unter welchen Umständen sich mit Reduktionsmethoden, die den
Einfluß von Transportprozessen auf das reduzierten Modell vernachlässigen, sinnvolle Er-
gebnisse erzielen lassen. Dies ist der Fall in allen Bereich des Zustandsraums, in welchen der
chemische Quellterm F (Ψ) deutlich größer als der Transportterm ρ1 d Ψθθ ◦ grad θ ◦ grad θ
und somit bestimmend für die Dynamik des Systems ist. Dort wird die Form und Dimen-
sion der durch die Lösungsmenge definierten Mannigfaltigkeit im wesentlichen durch F (Ψ)
bestimmt. Allerdings existieren prinzipiell immer auch Bereiche, in denen chemischer Quell-
term und Transportterm in die gleiche Größenordnung fallen und somit der Transportterm
wesentlich Form und Dimension der Mannigfaltigkeit mitbestimmt. Für solche Bereiche
führt die zusätzliche Berücksichtigung des Transportterms zu einer deutlichen Verbesser-
ung der Resultate [24].
Durch numerisches Lösen dieser Gleichung kann die REDIM Mannigfaltigkeit bestimmt
werden. Sie ergibt sich, wie bereits erwähnt, als stationäre Lösung von Gleichung 3.128. Als
3.7 BESCHREIBUNG DER REAKTIONSKINETIK 59
4
6E 05
5E 05
3
4E 05
H2O
HCO
2
3E 05
2E 05
1
1E 05
0 0
0 2 4 0 2 4
CO2
CO2
Abbildung 3.12: Simulation einer laminaren 1D Methan/Luft Flamme (schwarze Punkte: detail-
lierte Lösung, blau: Gradientenschätzung aus Flameletrechnung grün: Gerechnet mit im Vergleich
zu blau um Faktor 10 größeren Gradienten, rot: Gerechnet mit im Vergleich zu blau um Faktor
10 kleineren Gradienten) [24]
Anfangsbedingung dient dabei die erweiterte ILDM [25] was zu einer schnellern Konver-
genz und stabileren numerischen Lösung führt. Für das Lösungsschema wird als Parameter
lediglich eine Abschätzung für die typischen räumlichen Gradienten von θ benötigt. Wie
in [24] gezeigt wird, ist die die REDIM Methode wenig sensitv auf die tasächliche Wahl der
Gradienten. Eine geeignete Grandientenabschätzung kann zum Beispiel aus einer Flame-
letrechnung oder aus DNS Rechnungen turbulenter Flammen gewonnen werden.
In wie weit die REDIM zur reduzierten Beschreibung der Dynamik eines reaktiven Systems
eingesetzt werden kann, sollen die folgenden Resulate für eine eindimensionale laminare Ma-
than/Luft Flamme veranschaulichen. Verglichen werden in Abbildung 3.12 die Ergebnisse
einer mit einen detaillierten Mechanismus durchgeführten Simulationsrechnung (schwarze
Punkte) mit den mit einem durch die REDIM Methode reduzierten Mechanismus durch-
geführten Simulationsrechnungen (farbige Linien). Die Ergebnisse zeigen jeweils einen Teil
des Zustandsraums. Es ist zu kennen, daß sowohl für Majoritätenspezies wie hier darge-
stellt H2 O und CO2 als auch für Minoritätenspezies wie das HCO Radikal sich sehr gute
Übereinstimmungen mit den detaillierten Rechnungen ergeben. Insbesondere die gute Be-
schreibung der Radikale, die für die Dynamik von Verbrennungsprozessen sehr entscheidend
sein können, ist ein wesentlicher Vorteil des REDIM Verfahrens. Zusätzlich erkennt man
noch, daß die Methodik sehr insensitiv auf eine veränderte Abschätzung des Gradienten rea-
giert. Wie oben dargestellt ist diese Größe die a priori einzig notwendige Information über
das betrachtete System. Die lediglich geringe Abhänigkeit von der Gradientenabschätzung
unterstreicht die Robustheit und weite Einsatzmöglichkeit der verwendeten Methodik.
Aus den reduzierten Reaktionsmechanismen werden, zur Implementierung der Methodik
in Simulationen für reaktive Strömungen vorintegrierte, Lookuptabellen erstellt. Details zu
60 KAPITEL 3: NUMERISCHE SIMULATION REAKTIVER STRÖMUNGEN
Erstellung und Inhalt dieser Tabellen und ihrer Einbindung in das Gesamtverfahren finden
sich in Abschnitt 4.4.
Kapitel 4
Nach der Darstellung der einzelnen Modellteile in den vorherigen Abschnitten, stellt der
folgende Abschnitt das entwickelte Gesamtmodell zusammenfassend dar. Eine schematische
Darstellung zeigt zunächst auf wie die einzelnen Teile miteinander interagieren und welche
Daten sie miteinander austauschen. Zusätzlich werden die zur Verbesserung der numeri-
schen Stabilität und Beschleunigung der Konvergenz verwendeten Verfahren beschrieben.
Hierbei kommt dem Vermindern des statistischen Rauschens des Monte Carlo Verfahrens
eine wesentliche Bedeutung zu. Realisiert wird dies durch die geschickte Kombination von
Zeitmittelungs- und Iterationsmittelungsmethoden. Ebenfalls muss die Konsistenz zwischen
den Teilmodellen sichergestellt sein. Dies gilt insbesondere für doppelt auftauchende Größen
wie die Geschwindigkeit oder die Dichte. Die Konsistenz wird durch spezielle Korrektur-
alogorithmen gewährleistet. Abschließend zeigt dieser Abschnitt wie die automatisch redu-
zierten chemischen Kinetiken effizient tabelliert werden können und somit eine geeignete
Implementierung in das Gesamtmodell möglich ist.
61
62 KAPITEL 4: DARSTELLUNG DES GESAMTMODELLS
∂ p̄i
ũ, ṽ, w̃, ∂xi
, k̃, τ̃
φ(t), Δt
?
-
CFD PDF
REDIM
Navier-Stokes Gleichungen PDF Transportgleichung
Tabellenzugriff
Finite-Volumen Verfahren Monte Carlo Verfahren
6 φ(t + Δt)
T
M
,φ
die spezifische Molzahl von CO2 als Reaktionsfortschrittsvariable verwendet. Der zeitlich
Fortschritt der Reaktion wird durch Zugriff auf eine im Vorfeld erstellte Tabelle bestimmt.
Abbildung 4.1 zeigt ein Schema des Gesamtmodells. Die Simulation beginnt mit einem CFD
Teilschritt. Dabei erfolgt das Lösen der Navier-Stokes Gleichungen auf einem bockstruktu-
rierten Gitter mittels eines Finite-Volumen Verfahrens. Als Zwischenergebnis werden der
favregemittelte Geschwindigkeitsvektor, der mittlere Druckgradient und die turbulente ki-
netische Energie sowie das turbulente Zeitmaß für den aktuellen Zeitschritt an den PDF Teil
übergeben. Dort erfolgt unter Verwendung dieser Größen die Lösung der PDF Transport-
gleichung mittels eines auf stochatischen Partikeln basierenden Monte Carlo Verfahrens.
Der Fortschritt der chemischen Reaktion wird dabei durch einen Tabellenzugriff bestimmt.
Die Tabellen sind im Vorfeld der eigentlichen Simulation erstellt worden und bereits vor-
integriert, so dass der Reaktionsfortschritt ohne weitere Rechenoperationen direkt aus der
Tabelle ausgelesen werden kann. Der Reaktionsfortschritt wird nur für die reduzierten Zu-
standsvariablen berechnet. Der vollständige Zustandsvektor für den neuen Zeitschritt kann
ebenfalls unmittelbar aus der Tabelle ausgelesen werden. Zuruckgegeben an den CFD Teil
wird der favregemittelte Quotienten aus Temperatur und molarer Masse gemeinsam mit
dem vollständigen Zustandsvektor. Notwendig zur Kopplung des PDF Teils an den CFD
Teil ist nur der Quotienten aus Temperatur und molarer Masse. Dieser wird im CFD Teil in
die thermische Zustandsgleichung des idealen Gases eingesetzt und legt nach der Gleichung
T
p = ρ·R· (4.1)
M
den Druck fest. Der CFD Löser bestimmt das Dichte- und Geschwindigkeitsfeld unter Be-
rücksichtigung der Randbedingungen mit dieser neuen Temperatur so, dass die Impuls-
und Massenerhaltung erfüllt sind. Das benötigte Favremittel für den Quotienten wird im
PDF Teil bereits direkt berechnet. Kompressiblitätseffekte werden vernachlässig, da nur
4.2 NUMERISCHE STABILITÄT UND MITTELUNG 63
Strömungen mit relativ kleiner Machzahl untersucht wurden. Deshalb kann auf das Lösen
der Energieerhaltungsgleichung im CFD Teil verzichtet werden.
Das beschriebene iterative Verfahren wird für eine stationäre Simulation so lange wiederholt
bis die globale Konvergenz erreicht wird und im instationären Fall so lange wiederholt bis
die maximale Zahl an Zeitschritten erreicht ist.
FV Löser -
PDF Löser -
6 (feste Anzahl Iterationen) (feste Anzahl Iterationen)
6 imax
'
an+1 = 1
imax
ai
i=1
an+1 = β · an+1 + (1 − β) · an
Mittelung über die gesamten Iterationen (nmax )
Iterationen lag hier typischerweise bei etwa 100. Diese Zahl ist im Schema mit imax bezeich-
net. Die Größe a wird zur Reduzierung des oben beschriebenen stochastischen Rauschens
des partikelbasierten Lösungsverfahrens über alle PDF Iterationen nach der Beziehung
1 i
max
an+1 = ai (4.3)
imax i=1
gemittelt. Der neue Wert an+1 kann sich durch den im PDF Teil berücksichtigten Einfluss
von molekularer Mischung und chemischer Reaktion sehr stark von dem alten Wert an
unterscheiden. Würde dieser Wert unmittelbar an den FV Löser zurückgegeben könnten
dadurch starke numerisch verursachte Druckwellen entstehen, die zum Absturz des FV
Lösers führen würden. Um dies zu vermeiden und einen möglichst glatten Verlauf von a
zu erhalten, wird a zusätzlich noch mit einem gleitenden Durchschnitt über die einzelnen
Aufrufe des PDF Lösers gemittelt. Der neue Wert an+1 berechnet sich somit zu
Der Mittelwert des Partikelensembles in einer Zelle des Rechengebiets wird so korregiert,
dass er stets der mittleren Geschwindigkeit der CFD Lösung entspricht. Das aus der Im-
pulsgleichung durch ein Finite Volumen Verfahren bestimmte mittlere Geschwindigkeitsfeld
wird nicht durch stochatisches Rauschen des Lösungsverfahrens beeinflusst. Deshalb wird
die mittlere Geschwindigkeit des Partikelensembles auf diesen Wert korrigiert. Dabei wird
die Verwendung von Favre-gemittelten Größen in der Finite Volumen Methode entsprechend
berücksichtigt. Der erste Schritt hierbei ist die Bestimmung der Schwankungsgeschwindig-
keit jedes einzelnen Partikels.
ρi ui
ui ,alt = ui − (4.5)
ρi
Die neue (korrigierte) Partikelgeschwindigkeit berechnet sich aus
ukorr
i = ũF V + ui ,alt · A . (4.6)
66 KAPITEL 4: DARSTELLUNG DES GESAMTMODELLS
Nach dieser Korrektur stimmen die Mittelwerte der Finite Volumen und der Partikellösung
überein. Der Varianzkorrekturfaktor A wird so bestimmt, dass auch die Varianzen beider
Teilmodelle gleich sind (siehe unten).
Die im FV und PDF Teil verwendeten Submodelle, die die Varianz des Geschwindigkeits-
feldes abbilden, sind a priori inkonsistent. Eine Korrektur der Varianz des Geschwindig-
keitsfelds ist somit notwendig. Die Inkonsistenz der beiden Modelle beruht auf den unter-
schiedlichen Annahmen für den Turbulenzschluss. Es läßt sich zeigen, dass das im PDF Teil
verwendete SLM (Simplified Langevin Modell) auf der FV Ebene einem Turbulenzschluss
zweiter Ordnung mittes des Rotta-Modells entsprechen würde [127]. Benutzt wird im FV
Teil aber zum Turbulenzschluss lediglich ein Zweigleichungsmodell. Deshalb kann eine voll-
ständige Konsistenz der Varianz des Geschwindigkeitsfeldes nicht hergestellt werden1 . Aus
diesem Grund kann die Varianz des Partikelfeldes nur so korrigiert werden, dass die turbu-
lente kinetische Energie des Partikelfeldes und der Finite Volumen Lösung übereinstimmen.
Unter Berücksichtigung der Verwendung favregemittelter Größen soll der Korrekturalgo-
rithmus also dafür sorgen, dass die Bedingung
k̃P = k̃FV (4.7)
für jede Zelle des Rechengebiets erfüllt ist. Unter der Annahme isotroper Turbulenz wird
daraus die Bedingung
ui ,alt = uF V . (4.8)
Vor der Korrektur werden die beiden Terme (die beiden Seiten der Gleichung) im allge-
meinen Fall ungleich sein. Die Varianz des Partikelfeldes wird also mit einem noch zu
bestimmenden Korrekturfaktor A multipliziert, so dass die Gleichung erfüllt ist.
ui ,alt · A = uF V (4.9)
Die Schwankungsgröße uF V kann nicht unmittelbar aus den Ergebnissen der Finite Volumen
Rechnung bestimmt werden. Lediglich der Ausdruck
2
ρu2 F V = ρk̃F V (4.10)
3
kann aus der dichtegemittelten turbulenten kinetischen Energie ebenfalls wieder unter der
Annahme isotroper Turbulenz berechnet werden. Der zu ρu2 F V analoge Term ist die dich-
tegemittelte Varianz des Partikelgeschwindigkeitsfelds. Sie kann durch die Gleichung
' 2
ρuP
2
wi ρ i u i − ρP
ρu = ' (4.11)
P wi
1
Dieselbe Aussage gilt auch für die Reynoldsspannungen, auch hier kann keine vollständige Konsistenz
hergestellt werden, da das Zweigleichungsmodell lediglich die turbulente kinetische Energie und das turbu-
lente Zeitmaß, nicht aber direkt die einzelnen Reynoldsspannungen selbst, liefert.
4.3 KONSISTENZ UND KORREKTURALGORITHMEN 67
bestimmt werden. Der Term wi stellt dabei das numerische Gewicht eines Partikels dar. Bei
gleichem numerischen Gewicht aller Partikel wird wi für alle Partikel zu eins und der Term
zur Berechnung der Varianz des Partikelgeschwindigkeitsfeldes vereinfacht sich zu
2
2
ρuP
ρu = ρi ui − . (4.12)
P
ρP
Wird eine unverändert bleibenden Dichteverteilung in der jeweiligen Zelle des Rechengebiets
angenommen, so kann aus den oben abgeleiteten Termen der Korrekturfaktor A durch die
Beziehung
u V ρu2 F V 2
ρ k̃F V
A = F,alt = 2 = 3
(4.13)
ui
ρu P
ρu2 P
bestimmt werden. Als korrigierte Varianz des Partikelfeldes ergibt sich somit
ui ,neu = ui ,alt · A . (4.14)
führt. Die Korrekturgeschwindigkeit wird in dem von [103] vorgeschlagenen Verfahren aus
der Gleichung
∂φ ∂ Q̃ ∂Q
Uic =− − aU0 L ζ + (1 − ζ) (4.17)
∂xi ∂xi ∂xi
⎧
ρF V
⎨ 1 falls ρT A
≥ ζ
ζ=⎩ ρF V
(4.18)
0 falls ρT A
< ζ
berechnet. Hierin sind U0 und L eine charakteristische Geschwindigkeit und eine Längens-
kala, bei a handelt es sich um einen dimensionslose Parameter. Die hierfür verwendenden
Werte finden sich alle in Tabelle 4.1. φ stellt ein Korrekturpotential dar, dessen zeitliche
Entwicklung durch eine Differntialgleichung beschrieben wird.
∂φ
= bU02 Q (4.19)
∂t
In dieser Gleichung ist b ein dimensionsloser Parameter. Sein Wert findet sich ebenfalls in
Tabelle 4.1. Q̃(x, t) ist der in Ort und Zeit geglättete Dichteunterschied Q. Er wird ebenfalls
durch Lösen einer Differentialgleichung bestimmt. Sie lautet
∂ Q̃ U ∂ 2 Q̃
0
= − Q̃ − Q c + f U0 L (4.20)
∂t L ∂xi ∂xi
mit c und f als zusätzlichen dimensionslosen Parametern. Bei dieser Gleichung handelt es
sich um einen Zeitmittelungsoperator mit einem zusätzlichen Diffusionsterm. Auf diese Art
und Weise ist Q̃ das sowohl in der Zeit als auch im Ort geglättete Q. Beim Term ζ in
Gleichung 4.17 handelt es sich um eine Umschaltfunktion, die Q̃ durch Q ersetzt, so dass
der Algorithmus rasch reagieren kann sofern die mittlere Dichte des Finite Volumen Lösers
kleiner ist als die mittlere Dichte des Partikelfeldes.
Erreicht φ einen stationären Zustand, so gilt für den zeitlichen Mittelwert von Gleichung
4.19
QT A = 0 (4.21)
und somit auch
ρF V =
ρT A . (4.22)
Die Konsistenzbedingung wird folglich im zeitlichen Mittel exakt erfüllt sobald die statisti-
sche Stationarität erreicht wurde.
Wie bereits erwähnt finden sich in Tabelle 4.1 die in der Literatur verwendeten Werte für die
einzelnen in dem Korrekturverfahren vorkommenden Parameter. Für die Parameter kf , kb ,
NTc A , (CF L)P und ζ werden die in [103] vorgeschlagenen Werte 3, 8, 20, 0.4 und 0.25 ver-
wendet. Δx ist eine charakteristische Gitterweite und |U |max die maximale Geschwindigkeit
in der Lösungsdomäne.
4.4 BESTIMMUNG DES REAKTIONSFORTSCHRITTS 69
Δx
L π
U0 |U |max
1 1
c π (CF L)P NTc A
f kf c
b kb f = kb kf c
a 1 + cb2
ψ 0 = ψ 0 (θ) (4.23)
• die REDIM ergibt sich dann als stationäre Lösung von Gleichung 4.24 zu
Erzeugung einer
Erzeugung m-dimensionalen
einer Menge an Punkten Flameletrechungen
Anfangslösung
Interpolation auf ein
m-dimensionales Gitter
?
Berechnung der niedrigdimensionalen
Mannigfaltigkeit
?
Wahl einer geeigneten
Reaktionsfortschrittsvariablen
?
Tabelliere Reaktionsraten und andere
Informationen
ψ = ψ(θ) (4.26)
grad θ(θ) = ψθ+ (θ) grad ψ(θ) (4.27)
Die generalisierten Koordinaten ermöglichen später einen sehr effizienten Zugriff auf die
Tabelle. Als Koordinaten können zum Beispiel die Gitterindizes verwendet werden.
Die so bestimmte Mannigfaltigkeit bildet die Startlösung für die anschließenden REDIM
Iterationen. Hierbei wird Gleichung 4.24 mittels eines geeigneten numerischen Verfahrens
gelöst. Das numerische Verfahren muss große Systeme steifer Differentialgleichungen effektiv
lösen können. Die niedrigdimensionale REDIM Mannigfaltigkeit ergibt sich als stationäre
Lösung der Bestimmungsgleichung 4.24. Abbildung 4.5 stellt Startlösung und Endlösung
der Mannigfaltigkeit dar.
Im Anschluss muss eine geeignete Reaktionsfortschrittsvariable zur Parametrisierung der
REDIM ausgewählt werden. In vorliegenden Fall wird die spezifische Molzahl von CO2 als
Reaktionsfortschrittsvariable benutzt und die spezifische Molzahl von N2 zur Beschreibung
des lokalen Mischungszusandes verwendet. Um einen effizienten Zugriff auf die Tabellen zu
ermöglichen und gleichzeitig den Speicherplatzbedarf der Tabellen gering zu halten, wird die
REDIM mit einem problemangepassten Gitter diskretisiert. Abbildung 4.6(a) stellt einen
Ausschitt des Tabellierungsgitters dar. 2-Var steht hierbei für die spezifische Molzahl von
72 KAPITEL 4: DARSTELLUNG DES GESAMTMODELLS
CO2 und 1-Var für die spezifische Molzahl von N2 . Wie zu erkennen ist, wird das Gitter
im Bereich der Stöchiometrie, wo eine starke Dynamik der Reaktion zu erwarten ist, lokal
verfeinert. Abgelegt auf den Gitterpunkten ist, wie Abbildung 4.6(b) veranschaulicht, die
Reaktionsrate der Fortschrittsvariable zusammen mit dem hier nicht dargestellten vollstän-
digen Zustandsvektor.
Zur einfacheren Implementierung und Verbesserung der numerischen Genauigkeit bei dem
Zugriff auf den Reaktionsfortschritt wird die Tabelle der Reaktionsrate anschließend noch
über die Zeit integiert. Abgelegt werden dann die (zeitlich) neuen Werte zu verschiedenen
Zeitschritten (Zeitinkrementen). Abbildung 4.7 stellt eines dieser Zeitinkremente dar. Er-
kennbar ist unteranderem die unterschiedliche Geschwindigkeit der chemischen Reaktion
Abbildung 4.7: Werte der Reaktionsfortschrittsvariablen für einen definierten Zeitschritt (Zei-
tinkrement) von 128 μs
4.4 BESTIMMUNG DES REAKTIONSFORTSCHRITTS 73
Das entwickelte Simulationsmodell soll anhand von zwei verschiedenen stationären Flam-
mentypen validiert und anschließend zur Simulation instationärer Verbrennungsprozesse
eingesetzt werden. Um die breite Anwendbarkeit des Modells zu zeigen, wurden als Testfäl-
le sowohl eine vorgemischte als auch eine nicht-vorgemischte Flamme verwendet. In diesem
Kapitel werden zunächst die Ergebnisse der nicht-vorgemischten Flamme vorgestellt, die
Resultate für die vorgemischte Flamme finden sich in Kapitel 6 und die Simulationen der
instationären Verbrennungsprozesse in Kapitel 7.
Als Testfall einer nicht-vorgemischten Flamme wurde eine durch Masri et al. [92–95] experi-
mentell untersuchte Methan/Luft Flamme verwendet. Die Flamme wird durch einen Pilot-
brenner stabilisiert ist statistisch stationär. Ausgewählt wurde dieser Testfall aus zwei Grün-
den: zum einen ist in der Literatur eine detaillierte Datenbasis zur Validierung zu finden
und zum anderen zeigen numerische Untersuchungen anderer Arbeitsgruppen (z.B. [96]),
dass sich dieser Fall gut als Testfall für hybride CFD/transported PDF Modelle eignet.
Ein Vergleich mit diesen Arbeiten ermöglicht es zusätzlich die Güte des in dieser Arbeit
verwendeten Simulationsmodells bewerten zu können.
75
76 KAPITEL 5: ERGEBNISSE FÜR EINE STATIONÄRE NICHT-VORGEMISCHTE FLAMME
verwendet. Außen um die Brennstoffzufuhr herum ist zentrisch ein Pilotbrenner angeord-
net. Die Pilotflamme besteht aus einem nahezu vollständig abreagierten Ethin, Wasser-
stoff und Luft Gemisch. Das C zu H Verhältnis des Pilotbrenners entspricht dem eines
stöchiometrischen Methan/Luft Gemischs. Zentrisch um den Pilotbrenner herum sitzt die
Verbrennungsluftzufuhr.
Für den Brenner können verschiedene globale Betriebszustände durch Variation der Aus-
trittsgeschwindigkeit des Brennstoffs bei gleichbleibender Austrittsgeschwindigkeit des Pi-
lotbrenners und der Verbrennungsluft eingestellt werden. Die hierdurch entstehenden ver-
schiedenen Betriebszustände werden als Flamme K, L und M bezeichnet. Die genauen
Betriebsparameter lassen sich Tabelle 5.1 entnehmen.
Mit steigender Austrittsgeschwindigkeit des Brennstoff steigt von Flamme K über L bis
M die Turbulenzintensität an. Durch die erhöhte Turbulenz wird der Mischungsvorgang
zwischen den drei Teilströmen intensiviert. Das Verhältnis der turbulenten Zeitskala und
der chemischen Zeitskala beschrieben durch die Damköhlerzahl
τt
Da = (5.1)
τC
wird kleiner.
5.1 BESCHREIBUNG DES UNTERSUCHTEN BRENNERS 77
Der Testfall erlaubt somit eine Aussage darüber ob das Simulationsmodell über einen wei-
teren Bereich verschiedener Reynoldszahlen (Turbulenzintensität) und Damköhlerzahlen
anwendbar ist und zuverlässige Resultate liefert.
5.1.2 Gittergenerierung
Für die später dargestellten Simulationsergebnisse wurde ein zweidimensionales achsensy-
metrisches Rechengitter verwendet. Das Gitter besteht aus etwa 10000 Zellen. Abbildung
5.2 beinhaltet eine schematische Darstellung des Rechengitters und zeigt seine Lage in Re-
0,5 m
Luft
Pilot
Brennstoff
0,6 m
lation zum Brenner auf. Gezeichnet ist hier der Übersichtlichkeit halber nur jede vierte
Gitterlinie. Um Randeffekte weitestgehend ausschließen zu können ersteckt sich das Gitter
in x-Richtung über 0, 6 m und in y-Richtung über 0, 5 m. Das Gitter wurde in Bereichen, in
denen starke Gradienten zu erwarten sind lokal verfeinert.
5.1.3 Randbedingungen
Für das Rechengitter müssen an vier Seiten Randbedingungen spezifiziert werden. Am Ein-
lass wird eine Dirichletrandbedingung bestehend aus radialen Profilen für die Geschwindig-
keit, die Turbulenzgrößen und die Skalare verwendet. Abbildung 5.3 zeigt die verwendeten
Profile für die Axialgeschwindigkeit und die turbulente kinetische Energie. Die Profile sind
aus der Literatur [96] entnommen und stammen aus Messungen. Die Gradienten an den
Übergangen des Brenngasstroms und des Piloten sowie des Piloten und des Luftstroms sind
in analoger Weise zu den Siumlationen von [96] aus numerischen Gründen leicht abgeflacht
78 KAPITEL 5: ERGEBNISSE FÜR EINE STATIONÄRE NICHT-VORGEMISCHTE FLAMME
0.04 0.04
K-Flamme K-Flamme
L-Flamme L-Flamme
M-Flamme M-Flamme
0.035 0.035
0.03 0.03
0.025 0.025
y/m
y/m
0.02 0.02
0.015 0.015
0.01 0.01
0.005 0.005
0 0
0 10 20 30 40 50 60 70 0 5 10 15
u / m/s k / m²/s²
worden. Die erwähnte Literaturstelle zeigt auf, dass dies keinen signifikanten Einfluss auf
die Güte des Simulationsergebnisses hat. Die Differenz des gemessenen und des tätsächlich
verwendeten Profils ist sehr gering. Bei der Geschwindigkeit beispielsweise macht der Un-
terschied nur wenige Prozent aus. Für den Verlauf des turbulenten Zeitmaßes τ sind leider
keine Messungen vorhanden. Hierfür wird auf eine in der Literatur von Nau [104] für diesel-
be Flamme verwendet Abschätzung zurückgegriffen. Die Dissipationsrate der turbulenten
kinetischen Engergie kann danach durch
3 3
Cμ4 · k 2
= (5.2)
0, 25 · dchar
abgeschätzt werden und daraus kann dann nach der Formel
k
τ= (5.3)
das turbulente Zeitmaß errechnet werden. Hierin ist Cμ eine Modellkonstante mit dem Wert
0, 09 (siehe Abschnitt 3.3.3), k die turbulente kinetische Energie und dchar ein charakteri-
stischer Durchmesser. Als charateristischer Durchmesser dchar dient für den Brenngasstrom
sein Durchmesser, für den Piloten sein hydraulischer Durchmesser und für die Verbren-
nungsluft der Durchmesser der Auslassöffnung. Die Wahl der richtigen Randbedingung für
das turbulente Zeitmaß erwies sich in den Simulationen als kritisch. Wie die Ergebnisse in
den folgenden Abschnitten zeigen, wären hier weitere Informationen aus den experimen-
tellen Untersuchungen nötig gewesen. Deshalb zeigen die Ergebisse den deutlichen Einfluss
dieser Größe auf die Resultate in Form einer Parameterstudie auf.
An der Ober- und Unterseite des Rechengebiets wird eine Symmetrierandbedingung gesetzt,
was bedeutet, dass der Gradient aller Größen zu Null gesetzt wird. Am Auslass wird lediglich
5.2 ERGEBNISSE FÜR DIE K-FLAMME 79
der statische Druck vorgegeben und für alle anderen Größen ebenfall ein Gradient von Null
normal zur Austrittsfläche als Randbedingung vorgegeben.
typ [143] kommen zu dem Schluss, daß es sich bei Cφ nicht um eine feste Modellkonstante
handelt, sondern dass sie je nach betrachtetem Problem stark unterschiedliche Werte an-
nehmen kann. Aus diesem Grund wurde für die im folgenden gezeigten Parameterstudien
der Wert von Cφ variiert. Für die K-Flamme stellte sich bei einem Wert von 15 für Cφ die
beste Übereinstimmung mit den gemessenen Werten ein. Dieser Wert ist in guter Über-
einstimmung mit dem in [143] genannten Wert von 10 für die beste Reproduktion der
experimentellen Daten.
Abbildung 5.4 stellt die Simualtionsergebnisse der K-Flamme dar. In der linken Spalte sind
die Werte der ersten Messebene bei x = 0, 144 m und in der rechten Spalte die Werte
der zweiten Messebene bei x = 0, 36 m dargestellt. Die erste Zeile zeigt den Verlauf der
Axialgeschwindigkeit längs der y-Koordinate senkrecht zur Symmetrieachse des Brenners.
In der ersten Ebene stimmen die Simulationsdaten und die gemessenen Profile nahezu
vollständig überein. Lediglich in der zweiten weiter stromab liegenden Messebene ist nahe
an der Symmetrieline eine geringe Abweichung festzustellen. Generell wird das mittlere
Strömungsfeld durch die Simulation sehr gut wiedergegeben.
Der experimentelle Verlauf der Axialgeschwindigkeitsfluktuation kann unmittelbar aus den
Messdaten [96] entnommen werden. Als Vergleichsgröße dient die turbulente Schwankungs-
geschwindigkeit aus den Simulationsergebnissen. Sie wird aus der turbulenten kinetischen
Energie unter der Annahme isotroper Turbulenz berechnet. In der ersten Messebene stim-
men die beiden Verläufe noch sehr gut überein, in der zweiten Ebene zeigt sich bei
y/m < 0, 02 eine deutliche Abweichung. Der qualitative Verlauf wird jedoch richtig wie-
dergegeben. Allerdings sagen die Simulationsdaten eine um etwa 25% bis 30% zu hohe
Geschwindigkeitsfluktuation vorraus.
Die Profile für u und u sind jeweils für ein Cφ von 2, 1 dargestellt. Der Einfluss einer Var-
iation von Cφ auf diese beiden Größen ist sehr gering, weshalb die anderen Ergebnisse hier
nicht dargestellt sind. Unter Einbeziehung der unten gezeigten Temperaturverläufe, die sich
bei verschiedenen Werten für Cφ deutlich unterscheiden, ist dieser Befund zunächst über-
raschend. Zumindest für die Axialgeschwindigkeit würde man aufgrund der Kopplung über
die Dichte an die Temperatur einen deutlichen Einfluss erwarten. Eine genaue Analyse der
Simulationsergebnisse zeigt allerdings, dass das Axialgeschwindigkeitsfeld tatsächlich nahe-
zu unverändert bleibt, sogar bei einer Rechnung für eine nichtreagierde Strömung zeigt sich
ein zu den Messungen und reagierenden Simulationen nahezu identischer Verlauf. Allerdings
unterscheiden sich die Radialgeschwindigkeiten sehr massiv. Die Profile der Radialgeschwin-
digkeit sind wegen fehlender Messdaten nicht dargestellt. Es zeigt sich aber, dass die bei
Änderung der Temperatur offenkundige Änderung der axialen Impulsströme durch eine
radiale Ausgleichsströmung kompensiert wird, so dass die Profile der Axialgeschwindigkeit
nahezu unverändert bleiben. Dasselbe gilt auch bei den in den folgenden beiden Abschnitten
dieses Kapitels dargestellten Simulationsregebnissen für die L- und M-Flamme.
5.2 ERGEBNISSE FÜR DIE K-FLAMME 81
0.04 0.04
0.03 0.03
y/m
y/m
0.02 0.02
0.01 0.01
0 0
0 10 20 30 0 10 20 30
u / (m/s) u / (m/s)
0.04 0.04
0.03 0.03
y/m
y/m
0.02 0.02
0.01 0.01
0 0
0 1 2 3 4 5 6 0 1 2 3 4 5 6
u’ / (m/s) u’ / (m/s)
(b) Schwankungsgröße u (x = 0, 144 m) (e) Schwankungsgröße u (x = 0, 36 m)
0.04 0.04
c_phi = 2.1 c_phi = 2.1
c_phi = 5.1 c_phi = 5.1
c_phi = 10.0 c_phi = 10.0
c_phi = 15.0 c_phi = 15.0
Messung Messung
0.03 0.03
y/m
y/m
0.02 0.02
0.01 0.01
0 0
500 1000 1500 2000 500 1000 1500 2000
T/K T/K
Abbildung 5.4: Vergleich der gemessenen (Punkte) und simulierten (Linie) Profile für die
Flamme K
82 KAPITEL 5: ERGEBNISSE FÜR EINE STATIONÄRE NICHT-VORGEMISCHTE FLAMME
Gut zu erkennen ist er Einfluss der Variation des turbulenten Zeitmaßes bei den in der
untersten Zeile in Abbildung 5.4 dargestellten Temperaturverläufen. Die (turbulente) me-
chanische Zeitskala τt wird aus dem Turbulenzmodell des CFD Teils berechnet. Je größer
der Wert für Cφ wird, desto kleiner wird die (turbulente) chemische Zeitskala τφ im Ver-
gleich zur mechanischen Zeitskala. Im Mischungsmodell sorgt eine kleine chemische Zeits-
kala für mehr Mischungsvorgänge (siehe Gleichung 3.85). Beim Standardwert von 2, 1 wird
in der ersten Messebene oberhalb des Temperaturmaximums in y-Richtung eine zu ho-
he Temperatur vorrausgesagt und unterhalb des Temperaturmaximums eine zu niedrige
Temperatur. Mit steigendem Cφ nimmt durch vermehrte Mischungsprozesse der numeri-
schen Partikel die Temperatur oberhalb des Maximums weiter ab bis bei einem Cφ von 15
schließlich das Profil der Messung nahezugetroffen wird. Unterhalb des Maximums führt
die Steigerung von Cφ zu einer Zunahme der Temperatur und somit wird dort ebenfalls bei
Cφ = 15 der Verlauf der Messdaten sehr gut wiedergegeben. Unterhalb eines Punktes bei
etwa y/m = 0, 005 an dem sich die Linien der Simulationsdaten schneiden zeigt sich ein
umgekehrter Effekt: mit steigendem Cφ nimmt die Temperatur ab. Grund für beide Phäno-
mene sind die vermehrten Mischungsereignisse zwischen den Partikeln. Im Bereich oberhalb
des Schnittpunktes sind noch sehr viele heiße Partikel aus dem Pilotstrom enthalten. Ein
Mehr an Mischungsprozessen führt zu einer Steigerung der Temperatur. Im unteren Bereich
sieht es anders aus. Hier sind viele (kalte) Partikel aus dem Brennstoffstrom vorhanden. Die
Mehrzahl an Mischungsereignissen führt hier zu einem Verlöschen oder zumindest starkem
Abkühlen, der durch Turbulenz aus dem Pilotstrom eingemischten heißen Partikel. In der
zweiten Messebene zeigt sich oberhalb des Temperaturmaximums nur ein geringer Einfluß
der Konstante des Mischungsmodells. Lediglich bei Cφ = 15 ist die Temperatur etwas zu
gering. Bei den anderen Werten sind die Verläufe nahezu gleich. Schwach ist ebenfalls der
Trend einer abnehmenden Temperatur mit steigendem Cφ zu beobachten. Unterhalb des
Maximums nimmt die Temperatur mit steigendem Cφ wieder ab. Grund hierfür sind die
oben beschriebenen häufiger auftretenden Mischungsprozesse. Der Temperaturverlauf wird
ebenfalls durch ein Cφ von 15 über den gesamten Bereich gesehen am besten wiedergegeben.
Als Fazit kann festgehalten werden, dass sich Axialgeschwindigkeit und turbulente Schwan-
kungsgeschwindigkeit gut reproduzieren lassen. Lediglich bei dem Temperaturverlauf wirkt
sich die Unsicherheit hinsichtlich des turbulenten Zeitmaßes deutlich aus. Hier wären für
eine abschließende Beurteilung detailliertere Informationen aus den Messungen notwen-
dig gewesen. Der Wert für Cφ , bei dem die beste Übereinstimmung mit den gemessenen
Temperaturverläufen erzielt werden konnte, weicht zwar deutlich von dem in der Litera-
tur angegebenen Standardwert ab, deckt sich aber mit verschiedenen anderen Arbeiten an
ähnlichen Flammentypen [9, 143].
5.3 ERGEBNISSE FÜR DIE L-FLAMME 83
0.04 0.04
0.03 0.03
y/m
y/m
0.02 0.02
0.01 0.01
0 0
10 15 20 25 30 35 40 45 10 15 20 25 30 35 40 45
u / (m/s) u / (m/s)
0.04 0.04
0.03 0.03
y/m
y/m
0.02 0.02
0.01 0.01
0 0
0 2 4 0 2 4
u’
’ / (m/s) u’
’’ / (m/s)
(b) Schwankungsgröße u (x = 0, 144 m) (e) Schwankungsgröße u (x = 0, 36 m)
y/m
0.02 0.02
0.01 0.01
0 0
500 1000 1500 500 1000 1500
T/K T/K
Abbildung 5.5: Vergleich der gemessenen (Punkte) und simulierten (Linie) Profile für die
Flamme L
5.4 ERGEBNISSE FÜR DIE M-FLAMME 85
Zusammenfassend zeigt sich für die L-Flamme eine sehr gute Übereinstimmung der Simu-
lation mit den Messdaten für die mittlere Geschwindigkeit und die Geschwindigkeitsfluk-
tuation. Das Temperaturprofil kann durch die sich dort stark auswirkende experimentelle
Unsicherheit bei der Bestimmung einer zuverlässigen Randbedingung für das turbulente
Zeitmaß nur mit einer gewissen Unsicherheit berechnet werden.
0.04 0.04
0.03 0.03
y/m
y/m
0.02 0.02
0.01 0.01
0 0
10 15 20 25 30 35 40 45 50 10 15 20 25 30 35 40 45
u / (m/s) u / (m/s)
0.04 0.04
0.03 0.03
y/m
y/m
0.02 0.02
0.01 0.01
0 0
0 2 4 0 2 4
u’
’ / (m/s) u’
’’ / (m/s)
(b) Schwankungsgröße u (x = 0, 144 m) (e) Schwankungsgröße u (x = 0, 36 m)
y/m
0.02 0.02
0.01 0.01
0 0
500 1000 1500 500 1000 1500
T/K T/K
Abbildung 5.6: Vergleich der gemessenen (Punkte) und simulierten (Linie) Profile für die
Flamme M
5.4 ERGEBNISSE FÜR DIE M-FLAMME 87
Verwendung nur einer Fortschrittsvariablen bei der Modellierung der chemischen Reakti-
on eine Rolle spielen. Nahe an der Symmetrieachse liegt die höchste Turbulenzintensität
vor, die dazu noch deutlich über der maximalen Turbulenzintensität der K- und L-Flamme
liegt. Der Zustand der Flamme sollte hier also nahe an der Streckgenze liegen. Durch die
turbulenten Fluktuationen kann es in dieser Zone somit sehr wahrscheinlich zu lokalem
partiellen oder vollständigen Verlöschen der Flamme kommen. Um diese Zustände in der
Chemiemodellierung besser darstellen zu können müsste auf eine REDIM Tabelle zugegrif-
fen werden, die zwei oder mehr chemische Fortschrittsvariablen beinhaltet. Möglicherweise
könnten dadurch die gemessenen Temperaturfelder der M-Flamme deutlich besser durch
die Simulation reproduziert werden.
Zusammenfassend kann festgehalten werden, daß es auch bei der Flamme M sehr gut mög-
lich ist den Strömungszustand und die turbulenten Geschwindigkeitsflukationen mit dem
verwendeten Simulationsmodell darzustellen. Das Temperaturfeld läßt sich nur unzurei-
chend wiedergeben, da eine Information aus den Experimenten über den Verlauf des tur-
bulenten Zeitmaßes leider nicht vorhanden ist. Die Abweichungen im Temperaturfeld sind
bei der M-Flamme deutlicher als bei der K- und L-Flamme. Dies wird maßgeblich durch
die deutlich größere Turbulenzintensität verursacht.
Kapitel 6
89
90 KAPITEL 6: ERGEBNISSE FÜR EINE STATIONÄRE VORGEMISCHTE FLAMME
Hauptflamme stabilisiert. Die Luftzahl von Pilot und Hauptflamme ist gleich. Das Ex-
periment wurde mit einem stark verbreiterten Pilotbrenner durchgeführt um ein lokales
verlöschen und abheben der Flamme zu vermeiden. In den Experimenten wird zwischen
drei verschiedenen Flammentypen (Flamme F1, F2 und F3) unterschieden, die jeweils ei-
ne unterschiedliche Austrittsgeschwindigkeit des Hauptstroms haben. Es variiert somit der
Turbulenzgrad und das Verhältnis von physikalischem und chemischem Zeitmaß der ein-
zelnen Flammen untereinander. In Tabelle 6.1 sind die globalen Betriebsparameter der
verschiedenen Flammentypen zusammen gefasst.
6.1.2 Gittergenerierung
Für die folgend gezeigten Simulationsergebnisse wurde ein zweidimensionales achsensyme-
trisches Rechengitter verwendet. Das Gitter besteht insgesamt aus 3604 Zellen und ist in
Abbildung 6.2 in Relation zur Brennergeometrie dargestellt. Die Abbildung zeigt dabei
nicht das gesamte Rechengebiet sondern stellt nur den Ausschnitt in der unmittelbaren Nä-
he des Brenners dar. Um Randeffekte auf die Resultate minimieren zu können ersteckt sich
das Gitter in beiden Raumrichtungen über eine Länge von etwa 100 mal des Durchmessers
des Hauptgasstroms. Dargestellt ist hier nur ein kleiner Ausschnitt des Gitters. Lokale Ver-
feinerungen sollen die Gitterauflösung in Bereichen in denen starke räumliche Gradienten
6.1 BESCHREIBUNG DER UNTERSUCHTEN FLAMME 91
Pilot
Hauptgas
zu erwarten sind verbessern. Zusätzlich wurde der hinter der Austrittsebene des Brenners
liegende Bereich ebenfalls vernetzt um ein mitreißen der Umgebungsluft durch den ausströ-
menden Jet („Entrainment“) berücksichtigen zu können. Die später dargestellten Simula-
tionsrechnungen zeigten, dass dieser Effekt recht groß ist. Der Brenner saugt deutlich Luft
aus der umgebenden Gasströmung an. In den gezeigten Simulationen war die Vernetzung
dieses Bereichs notwendig um die experimentellen Daten validieren zu können.
6.1.3 Randbedingungen
Die zur Simulation notwendigen Randbedingungen konnten aus den in der Literatur veröf-
fentlichten Messdaten entnommen werden [26]. Vermessen wurden in der Austrittsebene des
Brenners in radialer Richtung von der Symmetrieachse bis zum Ende des Hauptgasstroms
die Axialgeschwindigkeit u (Abb. 6.3(a)) und die Geschwindigkeitsfluktuationen u und v .
u und v sind dabei die Geschwindigkeitsfluktuationen in axialer und radialer Richtung des
Brenners. Die nicht vermessene Geschwindigkeitsfluktuation w in Umfangsrichtung wurde
für die in den Simulationen verwendeten Randbedinungen zu Null gesetzt. Die turbulente
kinetische Energie (Abb. 6.3(b)) kann somit über die Beziehung
1 2
k= u +v2 (6.1)
2
berechnet werden. Um das turbulente Zeitmaß
k
τt = (6.2)
ε
92 KAPITEL 6: ERGEBNISSE FÜR EINE STATIONÄRE VORGEMISCHTE FLAMME
Flamme 1 Flamme 1
0.03 Flamme 2 0.03 Flamme 2
Flamme 3 Flamme 3
0.025 0.025
0.02 0.02
y/m
y/m
0.015 0.015
0.01 0.01
0.005 0.005
0 0
20 40 60 80 10 20 30 40 50
u / (m/s) k / (m²/s²)
am Rand zu erhalten, muss zunächst die Dissipationsrate der turbulenten kinetischen Ener-
gie ε nach der Gleichung
3
u2
ε= (6.3)
llat
berechnet werden. Das laterale integrale Längenmaß llat ist aus dem Experiment [26] für
die Mitte des Brenneraustritts bekannt (llat = 2, 4 mm) und wurde als konstant über den
Brennerradius angenommen.
Um die Einlassbedingungen des Pilotstroms richtig in der Simulation abbilden zu können
sind einige Vorüberlegungen notwendig. Die Austrittsebene des Piloten besteht aus einer
Lochscheibe mit 1175 Löchern mit einem Durchmesser von jeweils 1 mm. Experimentell
bestimmt wurde die Austrittsgeschwindigkeit aus einer Lochscheibe im Falle einer kalten
Strömung mit einem Wert von 0, 84 ms . In der Simulation wird als Randbedingung ein Block-
profil für die Axialgeschwindigkeit verwendet. Die verwendete Geschwindigkeit errechnet
sich aus dem in der Literatur angegebenen Versperrungsverhältnis und dem Quotienten
aus der Dichte der kalten und heißen Strömung. Somit ergibt sich ein errechneter Wert
von 4, 75 ms . Näheres hierzu findet sich in [160]. Die turbulente kinetische Energie und das
turbulente Zeitmaß wurden analog zu den Messdaten jeweils auf einen sehr kleinen Wert
gesetzt.
0.2 TA
/K
2100
2000
1900
1800
1700
0.15 1600
1500
1400
1300
1200
y/m
1100
1000
0.1 900
800
700
600
500 A B C D E
400
0.05
0
0 0.05 0.1 0.15
x/m
Tabelle 6.2: Axiale Position der Messebenen für Flamme F1, F2 und F3
94 KAPITEL 6: ERGEBNISSE FÜR EINE STATIONÄRE VORGEMISCHTE FLAMME
0.15
u
u / (m/s)
80
75
70
65
60
55
0.1 50
45
40
35
30
y/m
25
20
15
0.05 10
5
0
-5
-10
0 0.05 0.1
x/m
0.15
TT/K
2086
1974
1863
1751
1639
1527
0.1 1416
1304
1192
1080
969
y/m
857
745
633
0.05 522
410
0 0.05 0.1
x/m
Abbildung 6.6: Farbverlauf der der Temperatur für F1 als Vergleich von Experiment (unten)
und Simulation (oben)
Für die in diesem Abschnitt dargestellten Ergebnisse, als auch für die in den beiden folgen-
den Abschnitte gezeigten Ergebnisse für Flamme F2 und F3, wurde für das Mischungsmo-
dell jeweils mit dem Standardwert für Cφ von 2, 1 gerechnet. Eine Variation dieses Wertes
zur richtigen Wiedergabe der Temperatur- und Speziesverläufe in der Simulation war im
Gegensatz zu den Rechnungen der K-, L- und M-Flamme im vorherigen Kapitel hier nicht
notwendig. Der Grund dafür sind die deutlich detaillierteren Messdaten, die zur Bestim-
mung der Randbedingungen für die Simulationen herangezogen werden können. So ist das
turbulente Längenmaß am Auslass des Brenners gemessen worden und kann zur Bestim-
mung einer validen Randbedingung (siehe Abschnitt 6.1.3) für die im Turbulenzmodell
gelöste Gleichung für das turbulente Zeitmaß herangezogen werden. Die bei der K-, L- und
M-Flamme durch die Abschätzung des turbulenten Zeitmaßes verursachte Unsicherheit be-
steht für die Flammen F1, F2 und F3 somit nicht.
Ein Vergleich der gemessenen und berechneten radialen Profile in den fünf Messebenen
soll eine detaillierte Validierung der Simulationsergebnisse ermöglichen. Aus der Vielzahl
an vorhandenen Daten wurden zur Validierung des Strömungsfeldes die Messdaten für die
Axialgeschwindigkeit und die turbulente kinetische Energie in den Ebenen A bis C ausge-
wählt. Zur Validierung der skalaren Größen sollen exemplarisch die Temperatur und die
spezifische Molzahl von CH4 herangezogen werden. Dargestellt sind hier die Ergebnisse
in den Ebenen C bis E. Die anderen nicht dargestellten Ebenen zeigen jeweils qualitativ
ähnliche Ergebnisse [160]. Die Temperatur soll eine Orientierungsgröße für den Fortschritt
der chemischen Reaktion darstellen und zusätzlich den Einfluss der Einmischung des hei-
ßen Abgases aus dem Pilotbrenner aufzeigen. Die spezifische Molzahl von CH4 liefert eine
Information über den Abbau des aus dem Hauptgasstroms stammenden Brennstoffs.
96 KAPITEL 6: ERGEBNISSE FÜR EINE STATIONÄRE VORGEMISCHTE FLAMME
Der Vergleich zwischen gemessenen und berechneten Profilen ist in Abbildung 6.7 für die
Axialgeschwindigkeit und die turbulente kinetische Energie gezeigt. Aufgetragen ist jeweils
der Radius über der betrachteten Größe. Die einzelnen Zeilen stellen eine der Messebe-
nen A bis C dar. In der linken Spalte findet sich die Axialgeschwindigkeit in der rechten
jeweils die turbulente kinetische Energie. Die Symbole stellen die Messergebnisse dar, die
Linien die Resultate der Simulationen. In allen drei Messebenen stimmen die Profile der
Axialgeschwindigkeit von Messung und Simulation sehr gut überein. Sowohl der radiale
Verlauf als auch der absolute Wert werden gut getroffen. Bei der turbulenten kinetischen
Energie sind die Ergebnisse ähnlich gut, allerdings ergeben sich in allen drei Ebenen leichte
Abweichungen der Simulationsergebnisse. In der ersten Ebene nahe am Brenner wird die
radiale Höhe des Scherschicht leicht unterschätzt und somit liegt auch das Maximum der
turbulenten kinetischen Energie bei einem zu niedrigen y−Wert. Ab der mittleren Ebene
wird die Höhe wieder richtig berechnet. In Ebene B kann der Verlauf von k nahezu optimal
wiedergegeben werden. Lediglich unmittelbar an der Symmetrieachse ist der absolute Wert
zu klein. In der letzten dargestellten Ebene wird das Maximum der turbulenten kineti-
schen Energie um etwa 10 % überschätzt. Ebenso wird hier im Bereich der Symmetrieachse
ein zu niedriger Wert von k berechnet. Diskrepanzen der simulierten Werte im Bereich
der Symmetrieachse werden allem Anschein nach durch die in der Simulation notwendige
Randbedingung verursacht. Sie versuracht einen zur Symmetrieachse senkrechten Gran-
dienten von k. Zusammenfassend kann festgehalten werden, dass das durch die Simulation
berechnete Strömungsfeld sehr gut die experimentell gemessenen Werte wiedergibt.
Der Vergleich der Simulationsergebnisse für Temperatur und spezifische Molzahl von CH4
mit den dazugehörigen experimentellen Daten ist in Abbildung 6.8 dargestellt. In der lin-
ken Spalte der Abbildung ist jeweils der Abstand von der Symmetrieachse des Brenners
über der Temperatur aufgetragen, die rechte Spalte zeigt den Abstand über der spezifi-
schen Molzahl von CH4 . Die einzelnen Zeilen stellen die Messebenen C bis E dar. Das
Temperaturprofil kann in allen drei Ebenen sehr gut wiedergegeben werden. Lediglich in
der am weitesten vom Brenner entfernten Messebene wird der Maximalwert der Tempe-
ratur um knapp 10% unterschätzt. An der Symmetrieachse zeigt sich durch den Einfluss
der Randbedingung ebenfalls eine stärkere Abweichung zwischen den gemessenen und si-
mulierten Werten. Der prinzipielle Verlauf der spezifischen Molzahl von Methan wird wie
die rechte Spalte von Abbildung 6.8 zeigt in allen Messebenen qualitativ richtig wieder-
gegeben. Die an die experimentellen Datenpunkte gezeichneten Fehlerbalken geben die in
der Literatur angegebene Messunsicherheit von 10% bis 15% wieder [26]. Die experimen-
tellen Werte der spezifischen Molzahl basieren auf Ramanmessungen. Es ergeben sich in
allen Ebenen zumindest geringe Abweichung zwischen den experimentellen und den simu-
lierten Werten. Betrachtet man diese Abweichung, so liegen die simulierten Werte für alle
Messebenen zumindest nahezu innerhalb der angegebenen Messunsicherheit. Zudem muss
vermutet werden, dass die Messunsicherheit in Bereichen kleiner spezifischer Molzahlen,
6.2 ERGEBNISSE FÜR FLAMME F1 97
0.03 0.03
0.02 0.02
y/m
y/m
0.01 0.01
0 0
0 20 40 60 80 0 20 40 60 80 100 120 140
u / (m/s) k / (m²/s²)
0.02 0.02
y/m
y/m
0.01 0.01
0 0
0 20 40 60 80 0 20 40 60 80 100 120 140
u / (m/s) k / (m²/s²)
0.02 0.02
y/m
y/m
0.01 0.01
0 0
0 20 40 60 80 0 20 40 60 80 100 120 140
u / (m/s) k / (m²/s²)
Abbildung 6.7: Vergleich der gemessenen (Punkte) und simulierten (Linie) Axialgeschwindigkeit
und turbulenten kinetischen Energie in den Messebenen A bis C in des Brenners
98 KAPITEL 6: ERGEBNISSE FÜR EINE STATIONÄRE VORGEMISCHTE FLAMME
0.03 0.03
0.02 0.02
y/m
y/m
0.01 0.01
0 0
400 600 800 1000 1200 1400 1600 1800 0.5 1 1.5 2 2.5 3 3.5
T/K n_CH4 / (mol/kg)
0.02 0.02
y/m
y/m
0.01 0.01
0 0
400 600 800 1000 1200 1400 1600 1800 0.5 1 1.5 2 2.5 3 3.5
T/K n_CH4 / (mol/kg)
0.02 0.02
y/m
y/m
0.01 0.01
0 0
400 600 800 1000 1200 1400 1600 1800 0.5 1 1.5 2 2.5 3 3.5
T/K n_CH4 / (mol/kg)
Abbildung 6.8: Vergleich des gemessenen (Punkte) und simulierten (Linie) Temperaturprofils
und der spezifischen Molzahl des Brennstoffs in den Messebenen C bis E des Brenners
6.3 ERGEBNISSE FÜR FLAMME F2 99
wo die Diskrepanz der experimentellen und simulierten Werte in der Tendenz stärker ist,
durch die dort vorhandene geringe Signalstärke in der Messungen tendenziell eher größer
sein sollte. In Bereichen nahe der Strahlachse sind keine Messwerte vorhanden, so dass über
die Güte der Simulation in diesem Bereich keine unmittelbare Aussage möglich ist.
Zusammenfassend bleibt festzuhalten, dass die Strömungsgrößen sehr gut mit den experi-
mentellen Daten übereinstimmen. Somit liegt der Schluss nahe, dass die verwendeten CFD
Teilmodelle das untersuchte Strömungsproblem ausreichend gut beschreiben können. Die
gute Übereinstimmung bei Temperaturfeld und anderen Skalarfeldern zeigt zum einen, daß
die Kopplung zwischen CFD und PDF Teil des Modells auch für vorgemischte Flammen
gut funktioniert und sowohl die reaktionskinetischen Modelle als auch die Modelle zur Be-
schreibung des molekularen Transports auf Seiten des PDF Lösers für vorgemischte Systeme
gut funktionieren.
Um diese Aussagen zu untermauern wurden zusätzlich Simulationsrechungen für anderen
globale Betriebsparameter des Brenners (Flamme F1 und F2) durchgeführt. Die Resultate
hiervon sind in den folgenden beiden Abschnitten zu finden.
0.15
u / (m/s)
u
60
55
50
45
40
35
0.1 30
25
20
15
10
y/m
5
0
-5
0.05
0 0.05 0.1
x/m
0.15
T /TK
2100
2000
1900
1800
1700
1600
0.1 1500
1400
1300
1200
1100
y/m
1000
900
800
0.05 700
600
500
400
0 0.05 0.1
x/m
Abbildung 6.10: Farbverlauf der der Temperatur für F2 als Vergleich von Experiment (unten)
und Simulation (oben)
Ein ähnliches Bild zeigt sich für den Temperaturverlauf. Ein Vergleich der Simulationser-
gebnisse mit den experimentellen Daten ist in Abbildung 6.10 zu sehen. Ebenfalls stimmen
hier die Breite des Strahls sowie der qualitative Temperaturverlauf gut miteinander über-
ein. Der Beginn der Hauptflamme liegt aufgrund der geringeren Strömungsgeschwindigkeit
etwas niedriger als bei Flamme F1 und fällt etwa in die Mitte des durch Messungen ab-
gedeckten Bereichs. Der Flammenbeginn liegt in der Simulation in axialer Richtung etwas
zu hoch. Die Übereinstimmung mit der Lage im Experiment ist jedoch gut. Die absolu-
6.3 ERGEBNISSE FÜR FLAMME F2 101
ten berechneten Temperaturen stimmen, wie die Farbkodierung zeigt, ebenfalls über den
gesamten messtechnisch untersuchten Bereich, sehr gut mit den Messungen überein. Die
Einlasstemperatur des Pilotbrenners wird zur Berücksichtigung von Wärmeverlusten am
Brennerkopf herabgesetzt. Dieses Vorgehen deckt sich mit den in der Literatur veröffent-
lichten Ergebnissen anderer Arbeitsgruppen [143]. Auch dort ließen sich die experimentellen
Werte nur unter Berücksichtigung des Wärmeverlustes und der dadurch verursachten ge-
ringeren Eintrittstemperatur des Piloten reproduzieren.
Qualitativ gesehen ist die Güte der Simulationergebnisse für diesen Betriebszustand der
Flamme vergleichbar hoch wie für Flamme F1. Zu einer detaillierten Validierung der simu-
lierten Werte werden im folgenden die radialen Profile aus den Messungen in verschiedenen
Ebenen im Freistahlgebiet der Flamme herangezogen. Abbildung 6.11 zeigt diesen Vergleich.
In der linken Spalte finden sich die Geschwindigkeitsprofile von Messung und Simulation, in
der rechten Spalte die Profillinien der turbulenten kinetischen Energie. Die einzelnen Zeilen
stellen jeweils eine Messebene dar. Während das Geschwindigkeitsprofil in allen Ebenen
sehr gut berechnet werden kann, zeigt sich bei der turbulenten kinetischen Energie eine
gewisse Diskrepanz zwischen Simulation und Experiment. Besonders in der ersten Ebene
weichen die Werte stark von einander ab. In der Simulation wird zunächst die Produkti-
on der turbulenten kinetischen Energie deutlich überschätzt weshalb ein deutlich zu hoher
Wert herauskommt. Die turbulente kinetische Energie dissipiert danach aber rasch wieder,
so dass ab der zweiten Messebene B die beiden Profile wieder gut zu einander passen.
Auch in Ebene C sind lediglich geringe Abweichungen zu beobachten. Qualitativ wird der
Profilverlauf in allen Ebenen durch die Simulation korrekt wiedergegeben.
Aus den verschiedenen als Messdaten vorhandenen Skalarprofilen wurden die radialen Pro-
file der Temperatur und der spezifischen Molzahl von CH4 zusätzlich im Detail mit den
Simulationsergebnissen verglichen. Dargestellt ist dieser Vergleich in Abbildung 6.12. In
der linken Spalte sind die Temperaturprofile und in der rechten Spalte das Profil der spe-
zifischen Molzahl von CH4 dargestellt. Die Zeilen zeigen jeweils die Messebenen B bis D.
Im Bereich nahe der Symetrielinie bis etwa y = 0, 01 m ist die Übereinstimmung von ge-
messenem und berechnetem Temperaturverlauf sehr gut. Lediglich in Bereich des Maxi-
mums der Temperatur gibt es in allen Ebenen eine geringe Abweichung. Grund hierfür ist
voraussichtlich, dass es sich, wie oben dargestellt, als schwierig erweist exakte Randbedin-
gungen für den Pilotbrenner aus der Literatur abzuleiten (siehe Abschnitt 6.1). Dies fällt
bei Flamme F2 mehr ins Gewicht als bei Flamme F1 da die Ausströmgeschwindigkeit des
Hauptstroms geringer ist und somit auch der Enthalpie- und Impulsunterschied zwischen
beiden Strömen geringer wird. Erst weiter stromab im Rechengebiet wird also der Einfluss
der wahrscheinlich zumindest bedingt unphysikalischen Randbedingung abgeklungen sein.
Dadurch wirken sich die vorhandenen Unsicherheiten der Randbedingungen stärker auf die
Simulationsergebnisse aus. Bei den auf der rechten Seite in Abbildung 6.12 aufgetragenen
CH4 Profilen sind den Messwerten Fehlerbalken hinzugefügt worden. Die Breite der Fehler-
102 KAPITEL 6: ERGEBNISSE FÜR EINE STATIONÄRE VORGEMISCHTE FLAMME
0.03 0.03
0.02 0.02
y/m
y/m
0.01 0.01
0 0
0 20 40 60 80 0 20 40 60 80 100
u / (m/s) k / (m²/s²)
0.02 0.02
y/m
y/m
0.01 0.01
0 0
0 20 40 60 80 0 20 40 60 80 100
u / (m/s) k / (m²/s²)
0.02 0.02
y/m
y/m
0.01 0.01
0 0
0 20 40 60 80 0 20 40 60 80 100
u / (m/s) k / (m²/s²)
Abbildung 6.11: Vergleich der gemessenen (Punkte) und simulierten (Linie) Axialgeschwindig-
keit und turbulenten kinetischen Energie in den Messebenen A bis C der Flamme F2
6.3 ERGEBNISSE FÜR FLAMME F2 103
0.03 0.03
0.02 0.02
y/m
y/m
0.01 0.01
0 0
400 600 800 1000 1200 1400 1600 1800 0 0.5 1 1.5 2 2.5 3 3.5
T/K n_CH4 / (mol/kg)
0.02 0.02
y/m
y/m
0.01 0.01
0 0
400 600 800 1000 1200 1400 1600 1800 0 0.5 1 1.5 2 2.5 3 3.5
T/K n_CH4 / (mol/kg)
0.02 0.02
y/m
y/m
0.01 0.01
0 0
400 600 800 1000 1200 1400 1600 1800 0 0.5 1 1.5 2 2.5 3 3.5
T/K n_CH4 / (mol/kg)
Abbildung 6.12: Vergleich des gemessenen (Punkte) und simulierten (Linie) Temperaturpro-
fils und der spezifischen Molzahl des Brennstoffs in den Messebenen C bis E des Brenners im
Betriebspunkt F2
104 KAPITEL 6: ERGEBNISSE FÜR EINE STATIONÄRE VORGEMISCHTE FLAMME
balken entspricht der in der Literatur [26] angegebenen Messunsicherheit von 15 %. Auch
für Flamme F2 sind keine Messdaten nahe der Symmetrieachse vorhanden. Die berechneten
Profile passen qualitativ in allen drei Ebenen gut zu den experimentellen Daten. Im Bereich
mittlere Abstände von der Symmetrieachse liegen die Simulationsergebnisse innerhalb der
Messunsicherheit des Experiments. Lediglich in Bereichen etwas weiter weg von der Sym-
metrielinie ist ein geringer Unterschied zu den experimentellen Werten zu sehen. Ob in
diesem Bereich kleiner Absolutwerte von CH4 und somit auch geringerer Signalstärken im
Experiment immer noch derselbe prozentuale Fehler auftritt wird in der Literatur nicht
diskutiert, bleibt jedoch zumindest fraglich. Generell bleibt es somit schwierig die Güte der
Simulationsergebnisse in diesem Bereich der Flamme zu beurteilten. Es bleibt festzuhalten,
daß die Rechungen auch bei einer detaillierten Betrachtung der Profillinien eine sehr gute
Übereinstimmung mit den experimentellen Daten zeigen.
Das abschließende Fazit der Validierungsrechnungen für den Betriebszustand F2 fällt ähn-
lich positiv aus wie für Flamme F1. Die qualitativen Verläufe sowohl der Temperatur, der
spezifischen Molzahl von Methan, der Axialgeschwindigkeit als auch der turbulenten ki-
netischen Energie ließen sich sehr gut validieren. Die quantitativen Profile der genannten
Größen passen ebenfalls über einen weiten Bereich der Flamme sehr gut zu den experi-
mentellen Befunden. Abweichungen die ausserhalb der Messunsicherheit liegen sind bis zu
einem gewissen Grad auf die unvermeidbaren Unsicherheiten bei der Festlegung der Rand-
bedingungen zurückführen. Für die dargestellten Resulate konnte in allen Teilmodellen mit
den aus der Literatur bekannten Standardparametern gerechnet werden. Insbesondere gilt
dies für das Mischungsmodell im PDF Teil als auch für das Turbulenzmodell im CFD Teil,
für welche keine Anpassungen der Modellparameter vorgenommen worden sind.
0.15
u / (m/s)
u
35
32
29
26
23
20
0.1 17
13
10
7
4
y/m
1
-2
-5
0.05
0 0.05 0.1
x/m
0.15
T /TK
2100
2000
1900
1800
1700
1600
0.1 1500
1400
1300
1200
1100
y/m
1000
900
800
0.05 700
600
500
400
0 0.05 0.1
x/m
Abbildung 6.14: Farbverlauf der Temperatur für F3 als Vergleich von Experiment (unten) und
Simulation (oben)
läßt sich auch das Temperaturfeld der Flamme in der Simulation darstellen. In Abbildung
6.14 ist der Vergleich der Simulationsergebnisse und der experimentellen Daten gezeigt.
Beide Farbverläufe stimmen sehr gut miteinander überein. Der Beginn der Hauptflamme
wird richtig vorhergesagt. Er sitzt auf Grund der geringeren Austrittsgeschwindigkeit im
Vergleich zu F1 und F2 näher an der Austrittsebene des Brenners bei etwa x = 0, 05 m.
Um die Simulationsdaten im Detail besser mit den experimentellen Werten vergleichen
zu können sind im folgenden wieder Profillinien verschiedener Größen aus Simulation und
106 KAPITEL 6: ERGEBNISSE FÜR EINE STATIONÄRE VORGEMISCHTE FLAMME
Messung gegeneinander aufgetragen. Die Diagramme zeigen jeweils Ebenen die parallel zur
Austrittsebenen des Brenners liegen. In Abbildung 6.15 ist ein Vergleich der Profile von
Axialgeschwindigkeit und Temperatur in verschiedenen Höhen über dem Brenneraustritt
aufgetragen. In der linken Spalte ist die Geschwindigkeit und in der rechten Spalte die
Temperatur abgebildet. Die einzelnen Zeilen stellen jeweils eine Messebene oberhalb des
Brenners dar. Der angegebene Abstand entspricht dem Abstand zur Austrittsfläche. Aufge-
tragen in den Diagrammen ist der Abstand von der Symmetrielinie über der untersuchten
Größe. In den Ebenen A und B passen die simulierten Geschwindigkeitsprofile sehr gut zu
den Messergebnissen. Es sind lediglich miminale Abweichungen zu beobachten. In der letz-
ten Ebene C gibt es bei Radien größer als 0, 01 m eine geringe Diskrepanz zwischen dem aus
der Simulation erhaltenen und dem gemessenen Geschwindigkeitsverlauf. Der rechts darge-
stellte Temperaturverlauf passt in den Ebenen B und C nahe zu exakt zu den gemessenen
Profilen. Geringe Abweichungen zeigen sich lediglich in einem leicht zu großen Maximalwert
in Ebene C und in einer etwas zu großen Temperatur bei kleinen Radien in Ebene B. In
Messebene A hingegen ist die Abweichung zwischen Experiment und Simulation deutlich
größer. Zwar passt der Profilverlauf qualitativ sehr gut zu den Messergebnissen, quantitativ
ist jedoch doch eine deutliche Abweichung festzustellen. Hier wirken sich die vorhandenen
Unsicherheiten bei der Bestimmung der experimentellen Randbedingungen und ihren Ad-
aption für die Simulationen deutlich aus. So spielt möglicherweise der bereits im vorherigen
Abschnitt angeführte Wärmeverlust am Brennerkopf eine wesentliche Rolle.
Die nächste Abbildung soll verdeutlichen, dass zusätzlich zu den bisher dargestellten Größen
mit dem in dieser Arbeit entwickelten Simulationsmodell noch die Verteilung aller weiterer
Spezies berechnet werden kann. Prinzipiell lassen sich Werte für alle im detaillierten Reakti-
onsmechnismus vorkommende Spezies bestimmen. Ausgewählt wurden ein paar wenige für
die zusätzlich Messdaten vorhanden sind. Abbildung 6.16 stellt die Profile der spezifischen
Molzahl von CH4 , O2 und H2 O dar. Die linke Spalte zeigt radiale Profile in Ebene B, die
rechte Spalte radiale Profile in Ebene C. Die einzelnen Zeilen zeigen jeweils eine der Größen.
Aufgetragen ist in jeder Einzelabbildung der radiale Abstand von der Symmetrielinie über
der spezifischen Molzahl der jeweiligen Größe. Die Ergebnisse der Simulationsrechnungen
sind als Linie eingezeichnet und die Punkte stellen die Messdaten dar. Die eingezeichneten
Fehlerbalken bei den experimentellen Werten entsprechen der in der Literatur angeführten
experimentellen Unsicherheit von 15%. Betrachtet man zunächst das Profil von CH4 so
liegt der berechnete Verlauf bei nahezu allen Messpunkten innerhalb der experimentellen
Unsicherheit. Lediglich kleine Unterschiede sind zu beobachten. Das Methanprofil ergibt
sich als Überlagerung des Effekts der Lufteinmischung von außen und des Abbaus von
Methan durch die Reaktion. Ganz ähnlich sieht dies bei dem in der zweiten Zeile darge-
stellten Sauerstoffprofil aus. Es folgt ebenfalls aus dem Verbrauch von Sauerstoff durch die
chemische Reaktion und der Einmischung aus der Umgebungsluft. Die simulierten Werte
liegen ebenfalls für nahezu alle Punkte innerhalb des experimentell bestimmten Bereichs.
6.4 ERGEBNISSE FÜR FLAMME F3 107
0.03 0.03
0.02 0.02
y/m
y/m
0.01 0.01
0 0
0 20 40 60 80 400 600 800 1000 1200 1400 1600 1800 2000
u / (m/s) T/K
0.02 0.02
y/m
y/m
0.01 0.01
0 0
0 20 40 60 80 400 600 800 1000 1200 1400 1600 1800 2000
u / (m/s) T/K
0.02 0.02
y/m
y/m
0.01 0.01
0 0
0 20 40 60 80 400 600 800 1000 1200 1400 1600 1800 2000
u / (m/s) T/K
Axialgeschwindigkeit Temperatur
Abbildung 6.15: Vergleich der gemessenen (Punkte) und simulierten (Linie) Axialgeschwindig-
keit und der Temperatur in den Messebenen A bis C der Flamme F3
108 KAPITEL 6: ERGEBNISSE FÜR EINE STATIONÄRE VORGEMISCHTE FLAMME
0.03 0.03
0.02 0.02
y/m
y/m
0.01 0.01
0 0
0.5 1 1.5 2 2.5 3 3.5 0 0.5 1 1.5 2 2.5 3 3.5
n_CH4 / (mol/kg) n_CH4 / (mol/kg)
0.03 0.03
0.02 0.02
y/m
y/m
0.01 0.01
0 0
1 2 3 4 5 6 7 1 2 3 4 5 6 7
n_O2 / (mol/kg) n_O2 / (mol/kg)
0.03 0.03
0.02 0.02
y/m
y/m
0.01 0.01
0 0
1 2 3 4 5 6 7 8 1 2 3 4 5 6 7 8
n_H2O / (mol/kg) n_H2O / (mol/kg)
Ebene B Ebene C
Abbildung 6.16: Vergleich der gemessenen (Punkte) und simulierten (Linie) spezifischen Molzahl
verschiedener Skalare in den Messebenen B und C für Flamme F3
6.4 ERGEBNISSE FÜR FLAMME F3 109
Jedoch zeigen sich analog zum Verlauf von Methan auch hier kleine Unterschiede zwischen
Simulation und Experiment. Die in der letzten Zeile dargestellte spezifische Molzahl von
Wasser ist eine Kenngröße für den Reaktionsfortschritt. Er wird in beiden dargestellten
Ebenen richtig berechnet. Das durch Simulation bestimmte H2 O Profile ist generell etwas
niedriger als der gemessene Verlauf. Allerdings liegen die berechneten Werte noch innerhalb
der experimentellen Unsicherheit.
Als Fazit bleibt festzuhalten, dass sich auch die Flamme F3 sowohl qualitativ als auch quan-
titaiv sehr gut durch das in dieser Arbeit entwickelte Simulationsmodell beschreiben läßt.
Dies gilt grundsätzlich für alle Höhen über dem Brenneraustritt und auch für nahezu alle
betrachteten Größen. Das Modell eignet sich also prinzipiell zur Berechnung vorgemischter
Flammen und konnte mit den gezeigten Fällen F1 bis F3 über einen weiten Reynolds- und
Damköhlerzahlbereich validiert werden.
Kapitel 7
Simulation des
verbrennungsinduzierten
Wirbelaufplatzen
Die Validierungsrechungen in den beiden vorherigen Kapiteln zeigen, dass es mit dem in
dieser Arbeit entwickelten Modell grundsätzlich möglich ist statistisch stationäre turbulente
Flammen zu simulieren. Dies gilt sowohl für vorgemischte als auch für nicht-vorgemischte
Systeme. Der Strömungszustand der untersuchten Flammen umfasst dabei einen weiten Be-
reich unterschiedlicher Turbulenzgrade und Verhältnisse von chemischen und physikalischen
Zeitskalen. Die Übereinstimmung mit den experimentellen Daten war in allen Fällen sehr
gut. Nun soll gezeigt werden, dass sich das entwickelte Modell prinzipiell auch zur Beschrei-
bung instationärer Prozesse in turbulenten Flammen eingesetzt werden kann. Als Beispiel
für einen instationären Prozess wurden Flammenrückschläge durch verbrennungsinduziertes
Wirbelaufplatzen [44, 70] untersucht. Dieses Phänomen kann in vorgemischten Drallflam-
men auftreten, die nahe an der mageren Löschgrenze betrieben werden. Technisch relevant
sind solche Flammen insbesondere bei modernen schadstoffarmen Brennerkonzepten für
Gasturbinen wie zum Beispiel LPP1 Brennern. Das verbrennungsinduzierte Wirbelaufplat-
zen führt dabei dazu, dass die Flamme bei Schwankung der Betriebsparameter von ihren
stabilen Betriebspunkt in der Brennkammer sehr rasch und unaufhaltsam in die Vormisch-
strecke des Brenners wandert. Dort kann die Flammeneinwirkung zu massiven Schäden bis
hin zur vollständigen Zerstörung des technischen Systems führen.
Das folgende Kapitel fasst zunächst einige Grundlagen zur Phänomenologie des verbren-
nungsinduzierten Wirbelaufplatzens zusammen und stellt das untersuchte Brennersystem
dar. Dabei handelt es sich um ein eng an reale Brennersysteme angelehntes Modellsystem.
Die beiden folgenden Abschnitte des Kapitels zeigen dann eine detailierte Untersuchung
eines exemplarisch ausgewählten Rückschlagspunktes und eine Parameterstudie über den
1
Lean Premixed Prevaporized
111
112 KAPITEL 7: SIMULATION DES VERBRENNUNGSINDUZIERTEN WIRBELAUFPLATZEN
• Flammenrückschlag in Grenzschichten
In Grenzschichten kann die Strömungsgeschwindigkeit so gering werden, dass es zu
einem Stromaufpropagieren der Flamme entlang zum Beispiel der Grenzschicht einer
Wand kommen kann. Dieser Effekt ist durch Flammenlöschen in Folge von Quenching
begrenzt. Die Arbeiten von Lewis und von Elbe [76] formulieren ein Kriterium beste-
hend aus dem Geschwindigkeitsgradient an einer Wand und dem minimalen Abstand
(Quenchabstand) von einer Wand bei dem eine Flamme noch existieren kann und
treffen somit Vorraussagen ob es zum Flammenrückschlag in einer Wandgrenzschicht
kommen kann. Dieser Effekt tritt vorwiegend bei drallfreien niedrig turbulenten Strö-
mungen [45] oder lagsamen katalytischen Verbrennungen auf [121].
deutlich größer als bei einer laminaren Flamme ist. Eine größere Flammenfronto-
berfläche führt vereinfacht gesagt zu einem schnellern Umsatz des Brennstoffs pro
Zeit- und Volumeneinheit und somit zu einer höheren Brenngeschwindigkeit. Beob-
achtbar ist dieser Effekt zum Beispiel in hochturbulenten Drallflammen. Dort kann
in bestimmten Bereichen die turbulente Brenngeschwindigkeit größer als die mittle-
re Ausströmgeschwindigkeit des Brenners werden, was zu einem Flammenrückschlag
zum Beispiel auf der Brennerachse führen kann [53].
Für einen sicheren und störungsfreien Betrieb muss in technischen Anwendungen dafür
gesorgt werden, dass es nicht zu Flammenrückschlagen kommen kann. Bei manchen mo-
deren Brennerkonzepten für schadstoffarme Verbrennung (z.B. LPP, Lean Premixed Pre-
vaporized) kann das Problem auftreten, dass die dort verwendeten mageren vorgemischten
Drallflammen Flammenrückschläge zeigen, die nicht einer der oben genannten Kategorien
zugeordnet werden können. Die Flammen werden aerodynamisch durch einen definiert er-
zwungen Wirbelzerfall der Strömung am Brennkammereintritt stabilisiert. Bei leichter An-
fettung solch einer Flamme ausgehend von einem (stabilen) mageren Betriebspunkt kann
es zu einem Rückschlagen der Flamme in die Vormischstrecke kommen. Veranschaulicht ist
114 KAPITEL 7: SIMULATION DES VERBRENNUNGSINDUZIERTEN WIRBELAUFPLATZEN
dies in Abbildung 7.1 anhand eines vereinfachten Modells eines solchen Brenners. Das Bren-
nersystem besteht im wesentlichen aus drei Teilen: einem Plenum, einer Vormischstrecke
und der Brennkammer. Die Hauptströmungsrichtung ist in der Abbildung von links nach
rechts. Dem Plenum wird ein perfekt vorgemischtes Brennstoff/Luft Gemisch zugeführt,
welches anschließend über einen Drallrezeuger in die Vormischstrecke hinein strömt. Am
Ende der Vormischstrecke kommt es durch den Querschnittsprung und die ausreichend ho-
he Drallzahl zu einem Aufweiten der Hauptströmung. Ist, wie hier dargestellt, die Drallzahl
ausreichend hoch, kommt es zu einem Wirbelaufplatzen, das heißt durch die Zentrifugal-
kraft wird das Fluid nach außen hin weg von der Brennerachse gedrückt, die Dichte sinkt
lokal auf der Brennerachse ab und es entsteht eine entgegen der Hauptströmung gerichtete
Ausgleichsströmung. In dieser Rückströmzone (geschlossenen Rückströmblase) sammelt sich
heißes Abgas welches kontiuierlich dem aus der Vormischstrecke ausströmenden Frischgas
zugemischt wird und somit für ein stabiles Brennen der Flamme sorgt. Wird die Luftzahl
der Flamme nun ausgehend von diesem stabilen Betriebspunkt aus abgesenkt (Anfettung)
so kommt es bei Unterschreiten einer kritischen Luftzahl zum Rückschlagen der Flamme
in die Vormischstrecke. Dabei handelt es sich um einen inhärent transienten Prozess. Die
Flamme bleibt nicht an einer Stelle in der Vormischstrecke stehen sondern schlägt bis an den
Drallerzeuger zurück und stabilisiert sich erst dort wieder durch einen anderen Mechanis-
mus. Der Rückschlag erfolgt unmittelbar auf der Brennerachse mit einer Geschwindigkeit,
die um etwa eine Größenordung oberhalb der lokalen turbulenten Brenngeschwindigkeit
liegt, was einen Rückschlag in Folge turbulenter Flammenausbreitung in der Kernströmung
ausschließt. Ebenfalls zeigt sich in experimentellen Untersuchungen, dass keine akustischen
Instabilitäten zu beobachten sind [44, 70]. Das Wirbelaufplatzen der isothermen Strömung
und der Strömung mit magerer Flamme findet in der Brennkammer statt. Erst beim Anfet-
ten kommt es zu der Bewegung der Flamme in die Vormischstrecke hinein. Somit handelt
es sich auch nicht um den oben beschriebenen Flammenrückschlag durch (isothermen) Wir-
belzerfall. Als Abgrenzung davon wird dieser Rückschlagmodus als verbrennungsinduzierter
Wirbelzerfall (CIVB2 ) bezeichnet [44, 70].
desto größer wird auch der negative Impuls auf das Rückströmgebiet. Ist der Dichtegradient
hinreichend groß, oder gleichbedeutend die Flamme hinreichend fett, überwiegt der indu-
zierte Impulsstrom den Impuls der Hauptströmung und es kommt zum Stomaufpropagieren
der Flamme.
Die Hypothese stüzt sich auf Untersuchungsergebnisse an isothermen Drallströmungen, die
das entstehen des Wirbelzerfalls und des damit einhergehenden Staupunkts untersuchen
[21, 33, 40, 55, 74, 83]. Parameterstudien sollen den Einfluss globaler Betreibsparameter (wie
z.B. der Vorwärmtemperatur des Brennstoff/Luft Gemisches und der thermische Leistung
des Brenners) auf das kritische Luftverhältnis, bei welchen es zu einem Flammenrückschlag
durch verbrennungsinduzierten Wirbelzerfall kommt, zeigen.
Die numerische Simulation des verbrennungsinduzierten Wirbelzerfalls stellt einen sehr her-
ausfordernden Testfall für das in dieser Arbeit entwickelte Modell dar. So muss das mittlere
Strömungsfeld einer Drallflamme berechnet werden, die komplexe Interaktion von Chemie
und Turbulenz, die wichtig für das auftreten von CIVB ist, detailliert modelliert werden
und die chemische Kinetik muss auch nahe der mageren Löschgrenze noch gut beschrieben
werden.
Abbildung 7.3: Lage der Messebenen für die Axial- und Tangentialgeschwindigkeit in der Brenn-
kammer (A: x = 29 mm; B: x = 54 mm; C: x = 79 mm; D: x = 104 mm; E: x = 129 mm)
ratur der Luft von 298 K verwendet worden. Die Abbildung 7.3 zeigt als Farbverlauf die
Axialgeschwindigkeit in der Brennkammer. Simuliert wurde lediglich der Teil oberhalb der
Symmetrieachse, der untere Teil wurde zur besseren Übersichtlichkeit durch Spiegeln der
Resultate hinzugefügt. Die Simulation wurde für den nicht-reagierenden Fall durchgeführt.
Rein qualitativ ist die große geschlossene Rückströmzone zu beobachten, die an einem Stau-
punkt auf der Symmetrieachse am Ende der Vormischstrecke beginnt. Im reagierenden Fall
wird sich durch das in dieser Zone befindliche heiße, verbrannte Abgas die Flamme stabi-
lisieren können. In den fünf mit den Buchstaben A bis E gekennzeichneten Ebenen sind
Messdaten für die mittlere Umfangs- und Axialgeschwindigkeit vorhanden [44,70]. Ein Ver-
gleich von gemessenen und simulierten Geschwindigkeitsprofilen ist für die verschiedenen
Ebenen in Abbildung 7.4 und 7.5 dargestellt. In allen Abbildungen ist die Tangentialge-
schwindigkeit oder Axialgeschwindigkeit über dem Radius aufgetragen. In der linken Spal-
te befindet sich jeweils das Diagramm der Axialgeschwindigkeit auf der rechten Seite das
zugehörige Profile der Umfangsgeschwindigkeit. Die unterschiedlichen Zeilen in den bei-
den Abbildungen stehen für die verschiedenen Ebenen. Die roten Punkte markieren die
Messwerte, die blauen Linien stellen die Resultate der Simulationsrechnung dar. Für die
Axialgeschwindigkeit zeigt sich über alle Ebenen hinweg eine sehr gute Übereinstimmung
der Ergebnisse der Simulation mit den experimentellen Daten. Lediglich nahe der Symme-
trieachse (bei y = 0 m) ist eine Abweichung zu erkennen, die mit zunehmender Lauflänge
in der Vormischstrecke zunimmt. Die simulierte Geschwindigkeit ist dort etwas zu gering.
Eine mögliche Ursache hierfür könnte das verwendetet Zweigleichungsmodell sein. Zwei-
gleichungsmodelle führen in Drallströmungen oft zu einer Überschätzung der turbulenten
kinetischen Energie was zu einem abflachen der lokalen Geschwindigkeitsgradienten führt.
Die Tangentialkomponente wird ebenfalls sehr gut durch die Simulation wiedergegeben. Je-
doch kann die leichte Überhöhung der Tangentialgeschwindigkeit in den ersten vier Mess-
ebenen im Bereich der Symmetrieachse nicht richtig wiedergegeben werden. Die dort zu
erkennenden zu starke Dissipation der Umfangsgeschwidigkeit wird vermutlich ebenfalls
durch das verwendete Zweigleichungsturbulenzmodell verursacht. Diese Modelle verhalten
sich wie oben angedeutet und vielfach in der Literatur beschrieben bei Drallströmungen,
die in der Regel einen komplexen dreidimensionalen Turbulenzzustand mit nicht verschwin-
denen Normalspannungen haben, meist so, dass sie eine zu starke Abnahme der Umfangs-
geschwindigkeit vorraussagen. Dies ist hier ebenfalls zu beobachten. Dieses Defizit muss
aber, da aus rechenzeitgründen keine höherwertigen Turbulenzmodelle zum Einsatz kom-
men sollen, akzeptiert werden. Allerdings ist er für die Voraussage und Untersuchung des
Rückschlagverhaltens des Brenners wahrscheinlich von untergeordneter Bedeutung. Die für
den Rückschlag verantwortlichen wirbeldynamischen Prozesse finden erst deutlich hinter
der letzten dargestellten Messebene im Übergangsbereich zwischen Vormischstrecke und
Brennkammer statt. Wie die letzte vorhandene Messebene E vermuten lässt wird dort das
Profil der Umfangsgeschwindigkeit richtig wiedergegeben. In der letzten Messebene ist auch
7.4 DETAILIERTE ANALYSE EINES FLAMMENRÜCKSCHLAGS 119
0.04 0.04
0.03 0.03
y/m
y/m
0.02 0.02
0.01 0.01
0 0
0 10 20 30 40 50 0 10 20 30
u / (m/s) w / (m/s)
0.03 0.03
y/m
y/m
0.02 0.02
0.01 0.01
0 0
0 10 20 30 40 50 0 10 20 30
u / (m/s) w / (m/s)
0.03 0.03
y/m
y/m
0.02 0.02
0.01 0.01
0 0
0 10 20 30 40 50 0 10 20 30
u / (m/s) w / (m/s)
Axialgeschwindigkeit Tangentialgeschwindigkeit
Abbildung 7.4: Vergleich der gemessenen (rote Punkte) und simulierten (blaue Linie) Axial-
und Tangentialgeschwindigkeitsprofile in den Messebenen A bis C in der Vormischstrecke
120 KAPITEL 7: SIMULATION DES VERBRENNUNGSINDUZIERTEN WIRBELAUFPLATZEN
0.04
0.03
0.03
0.02
y/m
y/m
0.02
0.01 0.01
0 0
10 20 30 40 50 0 10 20 30
u / (m/s) w / (m/s)
0.04 0.04
0.03 0.03
y/m
y/m
0.02 0.02
0.01 0.01
0 0
0 10 20 30 40 50 0 10 20 30
u / (m/s) w / (m/s)
Axialgeschwindigkeit Tangentialgeschwindigkeit
Abbildung 7.5: Vergleich der gemessenen (rote Punkte) und simulierten (blaue Linie) Axial-
und Tangentialgeschwindigkeitsprofile in den Messebenen D und E in der Vormischstrecke
in der Messung die Überhöhung der Umfangsgeschwindigkeit dissipiert und die experimen-
tellen Daten und die Simulationsergebnisse zeigen eine sehr gute Übereinstimmung.
Das isotherme Strömungsfeld in der Brennkammer läßt sich somit wie gezeigt insgesamt
gut in der Simulation darstellen. Darauf aufbauend wurden die im folgenden gezeigten
Simulationen für den reagierenden Fall durchgeführt. Hierfür sind keine Messungen der Ge-
schwindigkeitsfelder vorhanden, die zur Validierung der Simulationsergebnisse herangezogen
werden könnten. Es sollte aber davon ausgegangen werden können, dass die Geschwindig-
keitsfelder auch im reagierenden Fall richtig wiedergegeben werden. Aufbauend auf den
isothermen Rechnungen soll zunächst für feste Werte der globalen Betriebsparameter ein
7.4 DETAILIERTE ANALYSE EINES FLAMMENRÜCKSCHLAGS 121
0.2 0.2
T T/ K
1900
1800
u
u / (m/s) 1700
0.1 0.1
1600
y/m
y/m
55
50 1500
44 1400
39 1300
33 1200
28 1100
0 0 1000
22
17 900
11 800
6 700
0 600
-5 500
-0.1 -0.1
0 0.1 0.2 0.3 0.4 0.5 0 0.1 0.2 0.3 0.4 0.5
x/m x/m
(a) λ = 1, 45 (d) λ = 1, 45
0.2 0.2
T T/ K
1900
1800
u
u / (m/s) 1700
0.1 0.1
1600
y/m
y/m
55
50 1500
44 1400
39 1300
33 1200
28 1100
0 0 1000
22
17 900
11 800
6 700
0 600
-5 500
-0.1 -0.1
0 0.1 0.2 0.3 0.4 0.5 0 0.1 0.2 0.3 0.4 0.5
x/m x/m
(b) λ = 1, 40 (e) λ = 1, 40
0.2 0.2
T T/ K
1900
1800
u
u / (m/s) 1700
0.1 0.1
1600
y/m
y/m
55
50 1500
44 1400
39 1300
33 1200
28 1100
0 0 1000
22
17 900
11 800
6 700
0 600
-5 500
-0.1 -0.1
0 0.1 0.2 0.3 0.4 0.5 0 0.1 0.2 0.3 0.4 0.5
x/m x/m
(c) λ = 1, 30 (f) λ = 1, 30
Die kleine Rückstromblase beginnt stromaufzuwandern und zieht die Flamme hinter sich
in die Vormischstrecke hinein. Durch die als Einlassrandbedingung fest vorgegebenen Ge-
schwindigkeitsprofile kommt die Flamme in der Simulation am Anfang der Vormischstrecke
zum stehen. Im Experiment schlägt die Flamme vollständig bis in das Plenum des Bren-
ners zurück. Vergleicht man die numerisch bestimmte kritische Luftzahl von λ = 1, 30 so
stimmt diese hervorragend mit dem im Experiment beobachteten Rückschlagsgrenze von
λkrit = 1, 25 − 1, 30 überein. Wie später in Abschnitt 7.5 noch gezeigt werden wird konnten
die experimentellen Rückschlaggrenzen nicht nur für diesen einen Betriebspunkt sondern
auch bei Variation der globalen Betriebsparameter reproduziert werden.
Um verstehen und erklären zu können wie und warum es zu einem verbrennungsinduzierten
Rückschlagen (CIVB) einer Drallflamme kommen kann sollen zunächst einmal die in der Li-
teratur zu findenen Erklärungen wie es in einer isothermen Strömung zum Aufplatzen eines
Wirbel kommt kurz betrachtet werden. Die Erklärung hierfür stützt sich im wesentlichen
auf die von Darmofal [33] durchgeführten Arbeiten, der die Rolle der wirbeldynamischen
Effekte auf das Wirbelaufplatzen systematisch untersucht.
Hinlänglich bekannt ist, dass es bei einer isothermen Wirbelströmung zum Beispiel in einer
Röhre zum Aufplatzen der Wirbelströmung oberhalb einer kritischen Drallzahl kommen
kann [55]. Gemeinhin wird als kritische Drallzahl ein Wert von 0, 5 angegeben. Verstärkt
wird die Neigung der Strömung zum Aufplatzen durch das vorhandensein von Querschnitts-
sprüngen. Eine gute Übersicht über die wirbeldynamischen Effekte finden sich zum Beispiel
in [40,55,74]. Aufbauend auf den Arbeiten von Lopez et al. [21,83] entwickelt Darmofal [33]
aus analytischen Betrachtungen und numerischen Simulationen die Theorie, dass es zu ei-
nem Wirbelzerfall und dem daraus resultierenenden Staupunkt und Rückströmgebiet in
einer Wirbelströmung stets in Bereichen großer negativer azimutaler Wirbelsträrke kommt.
Der Wirbelstärkevektor selbst ist definiert als die Rotation des Geschwindigkeitsvektors.
⎛ ∂w ⎞
∂y
− ∂v
∂z
⎜ ∂u ⎟
ω = ∇ × v = ⎜ ∂w ⎟
⎝ ∂z − ∂x ⎠ (7.1)
∂v ∂u
∂x
− ∂y
Die azimutale Wirbelstärke ist hierin die dritte Komponente des Wirbelstärkevektors.
∂v ∂u
ωz = − (7.2)
∂x ∂y
Abbildung 7.7(a) und 7.7(b) sollen zeigen, dass es auch hier im Bereich großer azimu-
taler Wirbelstärke zum Aufplatzen der Strömung kommt. Dargestellt ist der Verlauf der
azimutalen Wirbelstärke sowohl in einem stabilen Betriebszustand als auch während des
Flammenrückschlags. Zusätzlich ist zur Orientierung in beiden Abbildungen die Isolinie
der Axialgeschwindigkeit Null eingetragen. Der stabile Betriebszustand und der Zustand
während des Flammenrückschlags zeigen einen wesentlichen Unterschied bei der Lage des
124 KAPITEL 7: SIMULATION DES VERBRENNUNGSINDUZIERTEN WIRBELAUFPLATZEN
0.25 0.25
vorticity_z vorticity_z
20 20
17 17
15 15
0.2 12 0.2 12
9 9
6 6
4 4
1 1
0.15 -2 0.15 -2
-5 -5
-7 -7
y/m
y/m
-10 -10
0.1 0.1
0.05 0.05
0 0
0.1 0.2 0.3 0.1 0.2 0.3
x/m x/m
Gebiets negativer azimutaler Wirbelstärke. Im stabilen Betrieb ist die azimutale Wirbel-
stärke unmittelbar hinter dem Staupunkt auf der Symetrieachse und im vorderen Teil der
Rückströmzone negativ. Dies führt zum ortsfesten Aufplatzen der Stömung und somit zu
einem stabilen Betrieb des Brenners. Wird nun durch Anfetten der Flamme ein Rückschlag
initiert so ändern sich die Verhältnisse an der Spitze der Blase deutlich. Abbildung 7.7(b)
stellt den Verlauf der azimutalen Wirbelstärke während des Flammenrückschlags dar. Die
Zone negativer azimutaler Wirbelstärke liegt nun nicht mehr innerhab der Rückstömbla-
se, sondern stromauf des momentanen Staupunkts. Dort wird es lokal zum neuerlichen
Aufplatzen der Strömung kommen. Die Rückströmblase rutsch so ein Stück weiter in die
Vormischstrecke hinen und transportiert die Flamme mit sich. Dieser Mechnismus scheint
maßgeblich für den Transport der Flamme in die Vormischstrecke hinein und somit das
Rückschlagen der Flamme zu sein.
Zu untersuchen bleibt nun noch wie es durch Anfetten der Flamme zur Entstehung der
großen negativen azimutalen Wirbelstärke an der Blasenspitze kommen kann. Grundsätzlich
können aus der Wirbeltransportgleichung
∂ω 1
= (u · ∇) ω
+ u grad ω −ω
(∇ · u) + (∇ρ + ∇p) (7.3)
∂t
ρ2
Steckung Dilatation
baroklines Drehmoment
drei potentielle Quellterme3 für die Wirbelstärke identifiziert werden. Die Quellterme wer-
den durch Streckung, Dialtation und barokline Einflüsse auf die Wirbelfäden verursacht.
Die Arbeiten von Ashurst [4] über die dem Flammentransport in Wirbelröhren zu Grunde
3
die Betrachtungen sind hier für gravitationsfreie und reibungsfreie Strömungen, wodurch die in Folge
dieser Effekte verursachten Quellterme nicht berücksichtigt werden
7.4 DETAILIERTE ANALYSE EINES FLAMMENRÜCKSCHLAGS 125
0.25 0.25
s_exp_wz s_dilat_wz
5.0000E+04 2.5000E+00
0.2 4.0000E+04 0.2 2.0000E+00
3.0000E+04 1.5000E+00
2.0000E+04 1.0000E+00
1.0000E+04 5.0000E-01
0.0000E+00 0.0000E+00
0.15 -1.0000E+04 0.15 -5.0000E-01
-2.0000E+04 -1.0000E+00
-3.0000E+04
y/m
y/m
0.1 0.1
0.05 0.05
0 0
0.1 0.2 0.3 0.1 0.2 0.3
x/m x/m
0.25
s_baro_wz
1.0000E+06
0.2 -4.1667E+05
-1.8333E+06
-3.2500E+06
-4.6667E+06
-6.0833E+06
0.15 -7.5000E+06
-8.9167E+06
-1.0333E+07
y/m
-1.1750E+07
-1.3167E+07
-1.4583E+07
0.1 -1.6000E+07
0.05
0
0.1 0.2 0.3
x/m
Abbildung 7.8: Quellterme der Wirbeltransportgleichung für die azimutale Komponente wäh-
rend des Flammenrückschlags
tors dargestellt. Die exakten Gleichungen zur Berechung der Terme finden sich in Anhang
D. Dieser durch das Kreuzprodukt von lokalem Dichte- und Druckgradienten bestimmte
Term ist im Bereich der Blasenspitze an höchsten und wird durch anfetten der Flamme
und das damit einhergehende absinken der Dichte zusätzlich vergrößert. Somit erklärt sich
warum es durch anfetten zum einem Flammenrückschlag kommen kann. Der lokale Dichte-
gradient an der Blasenspitze wird immer größer, führt zu einem größeren Quellterm durch
das barokline Drehmoment und einer immer größer werdenden negative azimutalen Wirbel-
strärke an der Flammenspitze. Irgendwann wird der dadurch entstehende negative axiale
Impuls auf die Rückströmblase so groß, dass es zu einem Transport der Blase und somit
auch der Flamme stromauf in die Vormischstrecke hinein kommt.
Numerischen Untersuchungen anderer Arbeitsgruppen an derselben Geometrie [67] zeigen
die gleichen Resultate und untermauern ebenfalls die durch [44, 45, 70] aufgestellte Theorie
zur Entstehung eines Flammenrückschlags durch einen verbrennungsinduzierten Wirbelzer-
fall aus dem durch anfettung der Flamme steigendem baroklinen Drehmoment.
Abbildung 7.9: Vergleich numerisch und experimentell bestimmter Rückschlaggrenzen für ver-
schiedene globale Betriebsparameter
Luftzahl kann durch den bei höheren Durchsätzen gößeren positiven Geschwindigkeiten in
axialer Richtung erklärt werden. Dieser zusätzliche Impuls in positiver axialer Richtung
wirkt dem durch das barokline Drehmoment verursachten negativen Impuls entgegen. Ein
größerer Dichtegradient und somit eine fettere Mischung ist notwendig um die Rückström-
blase in die Vormischstrecke zu schieben und den Flammenrückschlag einzuleiten. Da es
bei höherer Vorwärmtemperatur schon bei geringeren Luftzahlen zu einem CIVB kommt
erklärt sich in ähnlicher Weise. Durch die höhere Vorwärmtemperatur steigt die Flammen-
temperatur an und der lokale Dichtegradient an der Flammenfront wird vergößert. Ein
größerer Dichtegradient führt zu einem größeren baroklinen Quellterm in der Wirbeltrans-
portgleichung und somit zu einem erhöhten negativen azimuthalen Wirbelstärke. Dieses
führt wiederrum zu einem gößeren negativen Impuls auf die Rückstromblase und läßt die
Flamme bereits bei magereren Luftverhältnissen zurückschlagen.
Kapitel 8
In dieser Arbeit wurde ein Modell zur Simulation stationärer und instationärer turbulenter
Strömungen mit chemischen Reaktionen entwickelt. Bei dem Modell handelt es sich um
ein hybrides CFD/transported PDF Modell. Es besteht im wesentlichen aus zwei Teilen:
einem Finite Volumen Löser für die Navier Stokes Gleichungen und einem auf stochati-
schen Partikeln basierenden Monte Carlo Löser für die Transportgleichung der gebunde-
nen Wahrscheinlichkeitsdichtefunktion von Geschwindigkeit und Zusammensetzung. Der
Schwerpunkt des Modells liegt auf der detaillierten Modellierung der kompelexen Interakti-
on von Turbulenz und chemischer Reaktion durch Lösen der PDF Transportgleichung. Die
chemische Kinetik wird durch automatisch reduzierte Reaktionsmechanismen beschrieben.
Als Reduktionsmethode wird das kürzlich entwickelte REDIM Verfahren eingesetzt [24].
Es erlaubt eine zuverlässige Beschreibung der Dynamik der chemischen Kinetik bereits mit
sehr wenigen Parametern und berücksichtig zusätzlich den Einfluss der Diffusion im Zu-
standsraum. Das in dieser Arbeit entwickelte PDF Modell basiert auf den Arbeiten von
Nau [104,105] und Bender [9,10]. Jedoch wurden einige substanzielle Änderungen und Wei-
terentwicklungen vorgenommen. So wurde ein neuer Strömungslöser an das bestehende PDF
Modell gekoppelt. Der Strömungslöser bietet den Vorteil verschiedender Erweiterungsmög-
lichkeiten, genannt sei hier nur die prinzipiell mögliche Verwendeung von LES auf der CFD
Seite. Der Strömungslöser arbeitet mit blockstrukturierten Gittern. Das PDF Modell wurde
um diese Funktionaliät erweitert, um auf diese Art und Weise auch komplexere Geometrien
abbilden zu können. Zusätzlich wurde noch eine Gleichung für die Umfangsgeschwindigkeit
im PDF Teil hinzugeführt um so axialsymmetrische drallbehaftete Strömungen berechnen
zu können. Außerdem wurden die Konvergenz des Gesamtvefahrens durch Untersuchung ge-
eigneter Mittelungsmethoden und die Verbesserung der Konsistenz der beiden Modellteile
durch umfangreiche neue Implementierungen weiter verbessert.
Das Modell sollte an einer möglichst vielfältigen Datenbasis an vorgemischten und nicht-
vorgemischten Flammen validiert werden. Anschließend wurde es zur Berechnung von Flam-
129
130 KAPITEL 8: ZUSAMMENFASSUNG UND AUSBLICK
dem untersuchten Beispielsystem handelt es sich um ein von Fritz [44] und Kröner [70] unter-
suchtes vereinfachtes Modell einer realen Gasturbinenbrennkammer. Die Flamme wird dort
durch Verdrallung der Hauptströmung und das sich dadurch im Bereich eines Querschnitts-
sprungs am Übergang der Vormischstrecke in die Brennkammer bildende Rückströmgebiet
stabilisiert. Der auftretende Rückschlagsmechanismus unterscheidet sich wesentlich von vie-
len bereits eingängig in der Literatur diskutierten Ursachen für Flammenrückschäge. Er
wird als verbrennungsindzierter Wirbelzerfall (CIVB) bezeichnet ist in der Literatur unter
anderem in [44, 67, 70] untersucht worden. Initiiert wird solch ein Rückschlag durch Anfet-
tung der Flamme ausgehend von einem stabilen (mageren) Betriebspunkt. Eine detaillierte
Analyse eines einzelnen Rückschlagereignisses zeigt, dass es durch anfetten zunächst zum
Abschnüren einer kleinen Rückstromblase von der Hauptrückstromzone kommt, die bei wei-
terer Anfettung in die Vormischstrecke hineinwandert und die Flamme mittransportiert.
Als treibender Mechanismus hierfür konnte die im Bereich der Spitze der Rückströmzo-
ne massiv produzierte negative azimutale Wirbelstärke identifiziert werden. Diese führt
zu einem starken Impuls in negativer x-Richtung auf die Rückströmblase. Als wesentliche
Quelle der negativen azimutalen Wirbelstärke zeigt sich in detaillierten Analysen der Si-
mualtionsergebnisse das barokline Drehmoment. Dieser Term wird durch das Kreuzprodukt
des lokalen Dichtegradienten und des lokalen Druckgradienten bestimmt. Diese Ergebnisse
sind in guter Übereinstimmung mit den Resultaten anderer Gruppen [158]. Zusätzlich wur-
den anhand experimenteller Untersuchungen Parameterstudien zum Einfluss der globalen
Betriebsparameter Einlassmassenstrom in die Brennkammer und Vorwärmtemperatur des
Brennstoff/Luft Gemischs durchgeführt. Ziel dieser Parameterstudie war die numerische Va-
lidierung des experimentell ermittelten Stabilitätsbereiches des Brenners. Die numerischen
Ergbnisse geben hierbei über den gesamten experimentell untersuchten Parameterbereich
die gemessenen kritschen Luftzahl, bei welcher es zum Flammenrückschlag kommt, sehr
gut wieder.
Die erfolgreichen Validierungsrechnungen für recht verschiedenartige Flammentypen und
die sehr erfolgreiche Anwendung des Gesamtmodells zur Simulation des verbrennungsindu-
zierten Wirbelzerfalls in mager vorgemischten Brennersystemen lassen an eine Reihe von
zukünftigen Ideen, Anwendungen und Erweiterungen des entwickelten Modells denken. So
kann zukünftig recht einfach eine Erweiterung des bisher auf zweidimensionale achsensym-
metrische Probleme beschränkten Modells auf drei Raumdimensionen durchgeführt werden.
Dadurch ließen sich auch komplexe reale Brennergeometrien simulieren. Der verwendetet
Strömungslöser kann bereits mit drei Dimensionen in Raum und Geschwindigkeit arbeiten,
lediglich der PDF Löser muß noch entsprechend erweitert werden. Die 3D Erweiterung wird
dann zukünftig auch den Einsatz von LES auf der Strömungslöserseite ermöglichen. Pas-
send dazu kann das PDF Verfahren so modifiziert werden, dass es die gefilterte Transport-
gleichung der Wahrscheinlichkeitsdichtefunktion (FDF) löst. Ein so entstehendes hybrides
LES/FDF Modell könnte dann sowohl auf der hydrodynamischen Seite durch den Einsatz
132 KAPITEL 8: ZUSAMMENFASSUNG UND AUSBLICK
von LES als auch bei der Modellierung der Chemie-Turbulenz-Interaktion auf hochwertige
Modelle zurückgreifen. Die durchgeführte Erweiterung auf instationäre Strömungsphäno-
mene lassen auch zukünftige Anwendungen des Modells im Bereich der Simulation von
selbstzündenden Freistrahlen zu. Hierbei ist es ebenfalls interessant die Möglichkeiten und
Grenzen der durch das REDIM Verfahren automatisch reduzierten Reaktionsmechanismen
zur Beschriebung der komplexen Kinetik bei Zündprozessen untersuchen zu können.
Anhang A
Die dargestellte Herleitung orientiert sich eng an den in [88] gemachten Ausführungen.
Weitere Details und Informationen zur Herleitung einer PDF Transportgleichung finden
sich unter anderem in [43, 124, 129].
Ausgangspunkt der Betrachtungen sind die Erhaltungsgleichungen für das Geschwindig-
keitsfeld (Impulserhaltungsgleichung) und das Skalarfeld (Spezieserhaltungsgleichung) in
ungemittelter Form.
• Impulserhaltungsgleichung
∂
ρ + U grad U = − div
Π
− grad p+ g = A (A.1)
∂t
• Spezieserhaltungsgleichung
∂
ρ + U grad φ = − divj + ρS = B (A.2)
∂t
Transport Reaktion
Diese Gleichungen konnen durch Einführen der substanziellen Ableitung einer beliebigen
Größe η
Dη ∂η
= + U grad η (A.3)
Dt ∂t
und dem Zusammenfassen der rechten Seite einfacher geschrieben werden als
DU
ρ = A (A.4)
Dt
Dφ
ρ = B . (A.5)
Dt
133
134 KAPITEL A: HERLEITUNG DER PDF TRANSPORTGLEICHUNG
Nun sei genau eine Realisierung einer Strömung betrachtet, bei welcher das Geschwindig-
keitsfeld und das Skalarfeld gegeben sind durch U und φ.
f ∗ (V, ψ; x, t) = δ (U1 (x, t) − V1 ) · · · δ (U3 (x, t) − V3 ) · δ (φ1 (x, t) − ψ1 ) · · · δ (φn (x, t) − ψn )
= δ (U (x, t) − V ) · δ (φ(x, t) − ψ) (A.6)
Die Änderung der PDF läßt sich durch ein vollständiges Differential angeben.
Df ∗ ∂f ∗ DU1 ∂f ∗ DU3 ∂f ∗ Dφ1 ∂f ∗ Dφn
ρ = ρ + ··· + ρ +ρ + ··· + ρ
Dt ∂U1 Dt ∂U3 Dt ∂φ1 Dt ∂φn Dt
∂f ∗ DU ∂f ∗ Dφ
= ρ +ρ
∂U Dt ∂φ Dt
∂δ(U − V )δ(φ − ψ) DU ∂δ(U − V )δ(φ − ψ) Dφ
= + (A.7)
∂U Dt ∂φ Dt
Nach Anwenden der allgemeinen Rechenregel
∂ ∂
δ(x − y) = − δ(x − y) (A.8)
∂x ∂y
folgt somit
Df ∗ ∂f ∗ ∂f ∗
ρ =− A− B . (A.9)
Dt ∂V ∂ψ
Bis zu diesem Punkt wurde immer noch eine deterministische Gleichung, die exakt eine
Realisierung der Strömung darstellt, betrachtet. Die Strömung besteht aber tatsächlich aus
einer Vielzahl verschiedener Realisierungen. Deshalb ist die Bildung eines Erwartungswertes
sinnvoll und notwendig.
+∞
∗
f = f ∗ (V, ψ) f (V̂ , ψ̂) dV̂ dψ̂
−∞
+∞
Anschaulich betrachtet stellt die PDF somit den Erwartungswert der PDFs der Einzel-
realisierungen dar. Wird nun der Erwartungswert von Gleichung A.9 gebildet, so ergibt
sich
Df ∗ ∂f ∗ ∂f ∗
ρ = − A− B . (A.11)
Dt ∂V ∂ψ
Der Erwartungswert der linken und rechten Seite dieser Gleichung kann durch Einsetz-
ten der in Gleichung A.10 gezeigten Rechenregel durch Integration über den gesamten
Zustandsraum bestimmt werden. Es folgt dann der Ausdruck
+∞
Df ∗ (V, V̂ , ψ, ψ̂)
ρ(ψ) f (V̂ , ψ̂) dV̂ dψ̂ =
Dt
−∞
⎡ ⎤
+∞
∂f ∗ (V, V̂ , ψ, ψ̂) ∂f ∗ (V, V̂ , ψ, ψ̂) ⎦
f (V̂ , ψ̂) ⎣− A− B dV̂ dψ̂ . (A.12)
∂V ∂ψ
−∞
135
D
+∞
umgeformt werden und vereinfacht sich durch anwenden von Gleichung A.10 weiter zu
D D
ρ(ψ) f (V̂ , ψ̂) f ∗ (V, V̂ , ψ, ψ̂) = ρ(ψ)
f ∗ . (A.14)
Dt Dt
Die rechte Seite von Gleichung A.12 vereinfacht sich durch dieselben Umformungen in
analoger Weise zu
+∞
+∞
∂f ∗ (V, V̂ , ψ, ψ̂) ∂f ∗ (V, V̂ , ψ, ψ̂)
RS = − f (V̂ , ψ̂) A dV̂ dψ̂ − f (V̂ , ψ̂) B dV̂ dψ̂
∂V ∂ψ
−∞ −∞
∂ ∗
+∞ +∞
∂ ∗
= − f (V, V̂ , ψ, ψ̂) f (V̂ , ψ̂) A dV̂ dψ̂ − f (V, V̂ , ψ, ψ̂) f (V̂ , ψ̂) B dV̂ dψ̂
∂V ∂ψ
−∞ −∞
∂ ∂
= −
f ∗ A −
f ∗ B . (A.15)
∂V ∂ψ
Der Erwartungswerte von Gleichung A.9 wird folglich zu
D ∗ ∂ ∂
ρ(ψ)
f = −
f ∗ A −
f ∗ B . (A.16)
Dt ∂V ∂ψ
Ungeschickt ist an dieser Stelle, dass noch die Einzel-PDF auftritt. Sie kann allerdings
durch den Zusammenhang
f ∗ = f (ψ) in geeigneter Weise ersetzt werden.
D ∂ ∂
ρ(ψ) f =−
f ∗ A −
f ∗ B (A.17)
Dt ∂V ∂ψ
f ∗ A(U, φ, grad U, grad φ) =
δ(U − V ) δ(φ − ψ) A(U, φ, grad U, grad φ) (A.18)
Beide Terme
f ∗ A und
f ∗ B stellen also bedingte Erwartungswerte dar. Allgemein kann
ein bedingter Erwartungswert berechnet werden durch
+∞
+∞
Q = Q(U, φ, · · · )|U = V̂ , φ = ψ̂ fU φ (V̂ , ψ̂) dV̂ dψ̂ . (A.19)
−∞ −∞
+∞
+∞
= δ(U − V )δ(φ − ψ)A(U, φ, grad U, grad φ)|U = V̂ , φ = ψ̂ fU φ (V̂ , ψ̂) dV̂ dψ̂
−∞ −∞
+∞
+∞
= δ(V̂ − V )δ(ψ̂ − ψ)A(U, φ, grad U, grad φ)|U = V̂ , φ = ψ̂ fU φ (V̂ , ψ̂) dV̂ dψ̂
−∞ −∞
=
A(U, φ, grad U, grad φ)|U = V, φ = ψ fU φ (V, ψ) . (A.20)
Df (V, ψ) ∂
ρ(ψ) = − fU φ (V, ψ)
A(U, φ, grad U, grad φ)|U = V, φ = ψ
Dt ∂V
∂
− fU φ (V, ψ)
B(U, φ, grad U, grad φ)|U = V, φ = ψ (A.21)
∂ψ
Die beiden hierbei auftauchenden bedingten Erwartungswerte sind ungeschlossene Terme.
Mögliche Modellierungsansätze für diese Terme wurden in Abschnitt 3.6.5 dargestellt.
Anhang B
In Tabelle B.1 sind die Koeffiziente des verallgemeinerten Langevin Modells dargestellt. Sie
werden aus experimentellen Daten gewonnen. Weitere Details hierzu findet sich in [58,124].
Die Modellkonstante γ ∗ ergibt sich aus den angegebenen Konstanten durch die Beziehung
γ ∗ = γ2 + γ 3 + γ 4 + γ 6 . (B.1)
α1 − 1
2
+ 34 C0 − α2 bij bji − τ γ ∗ bkl bil ∂U
∂xl
k
α2 3,7
β1 − 15
4
β2 5
β3 − 15
γ1 -1,28
γ2 3,01
γ3 -2,18
γ4 0
γ5 4,29
γ6 -3,09
137
Anhang C
0.03
y/m
0.02
0.01
0 20 40
u / (m/s)
139
140 KAPITEL C: RANDBEDINGUNGEN FÜR DIE FLAMMENRÜCKSCHLAGSRECHUNGEN
0.03
y/m
0.02
0.01
0 5 10 15 20 25
w / (m/s)
0.03
y/m
0.02
0.01
0 5 10 15
k / (m²/s²)
0.03
y/m
0.02
0.01
141
142 KAPITEL D: BERECHNUNG DER QUELLTERME DER WIRBELTRANSPORTGLEICHUNG
[2] Ahmed, S.S. ; Mauss, F. ; G., Moreac ; T., Zeuch: A comprehensive and compact n-
heptane oxidation model derived using chemical lumping. Physical chemistry chemical
physics 9 (2007), Nr. 9, S. 1107–1126
[3] Arrhenius, S.: Über die Reaktionsgeschwindigkeit bei der Inversion von Rohrzucker
durch Säuren. Zeitschrift Physikalische Chemie 4 (1889), S. 226
[4] Ashusrt, WM. T.: Flame Propagation Along a Vortex: the Baroclinic Push. Com-
bustion, Science and Technology 112 (1996), S. 175–185
[6] Baldwin, B.S. ; Lomax, H.: Thin Layer Approximation and Algebraic Model for
Separated Turbulent Flows. AIAA Journal 78 (1978), S. 257
[8] Bauer, H.J.: New Low Emission Strategies and Combustion Designs for Civil Ae-
roengine Applications. Progress in Computational Fluid Dynamics 4 (2004), Nr. 3-5,
S. 130–142
[9] Bender, R.: Modellierung der Kopplung von chemischer Reaktion und turbulen-
ter Mischung bei turbulenten Vormischflammen, Universität Stuttgart, Fakultät für
Energietechnik, Dissertation, 2003
143
144 LITERATURVERZEICHNIS
[11] Bird, R.B. ; Steward, W.E. ; Lightfoot, E.N.: Transport phenomena. New York,
Chichester, Brisbane, Toronto, Singapore : John Wiley & Sons, 1960
[12] Bird, R.B. ; Stewart, W.E. ; Lightfoot, E.N.: Transport phenomena. 2. ed..
Aufl. New York : Wiley, 2002. – ISBN 0–471–41077–2 ; 0–471–36474–6
[14] Blouch, J.D. ; Chen, J.Y. ; Law, C.K.: A joint scalar PDF study of nonpremixed
hydrogen ignition. Combustion and Flame 135 (2003), S. 209–225
[16] Borghi, R.: On the Structure and Morphology of Turbulent Premixed Flames.
Recent Advance in the Aerospace Sciences (1984)
[17] Borghi, R.: Turbulent Combustion Modelling. Progress in Energy and Combustion
Science 14 (1988), S. 245–292
[18] Boussinesq, J.: Theorie de lâĂŹEcoulement Tourbillant. Mem. Presentes par Divers
Savants Acad. Sci. Inst. Fr. 23 (1877), S. 46–50
[20] Bronstein, I.N. ; Semendjajew, K.A.: Taschenbuch der Mathematik, 25. 1991
[21] Brown, G.L. ; Lopez, J.M.: Axisymmetric vortex breakdown part 2: physical me-
chanisms. Journal of Fluid Mechanics 221 (1990), S. 553–576
[22] Bykov, V. ; Gol’dshtein, V. ; Maas, U.: Global Quasi Linearization (GQL) for the
automatic reduction of chemical kinetics. In: Proceedings of the European Combustion
Meeting 2007, 2007
[24] Bykov, V. ; Maas, U.: The extension of the ILDM concept to reaction-diffusion
manifolds. Combustion Theory and Modelling 11 (2007), Nr. 6, S. 839–862
LITERATURVERZEICHNIS 145
[25] Bykov, V. ; Maas, U.: Extension of the ILDM method to the domain of slow
chemistry. In: Proceedings of the Combustion Institute Bd. 31, 2007, S. 465–472
[28] Choi, K.S. ; Lumley, J.L.: The return to isotropy of homogeneous turbulence.
Journal of Fluid Mechanics 436 (2001), S. 59–84
[29] Coffey, W.T. ; Kalmykov, Y.P. ; Waldron, J.T.: The Langevin Equation: with
applications in Physics, Chemistry and Electrical Engineering. 1996 (World Scientific
Series in Contemporary Chemical Physics - Vol. 10)
[30] Comte-Bellot, G. ; Corrsin, S.: The use of a contraction to improve the isotropy
of grid-generated turbulence. Journal of Fluid Mechanics 25 (1966), S. 657–682
[31] Comte-Bellot, G. ; Corrsin, S.: Simple Eulerian time correlation of full and
narrow-band velocity signals in grid-generated ’isotropic’ turbulence. Journal of Fluid
Mechanics 48 (1971), S. 273–337
[32] Curl, R.L.: Dispersed Phase Mixing: 1. Theory and Effects in Simple Reactors.
A.I.Ch.E. Journal 9 (1963), S. 175–181
[33] Darmofal, D.L.: The Role of Vorticity Dynamics in Vortex Breakdown. AIAA
Journal (1993)
[34] Dopazo, C.: Relaxation of Initial Probability Density Functions in the Turbulent
Convection of Scalar Flieds. Physics of Fluids 22 (1979), S. 20–30
[35] Dopazo, C.: Recent developments in pdf methods. In: Libby, P.A. (Hrsg.) ; Wil-
liams, F.A. (Hrsg.): Turbulent Reacting Flows. Academic Press, New York, 1994, S.
375–474
[36] Drain, L. E. ; Monzingo, T. W. R. A. & Miller A. R. A. & Miller Ṙ. A. & Miller M.
R. A. & Miller (Hrsg.): The laser Doppler techniques. 1980
[37] Dreeben, T.D. ; Pope, S.B.: PDF/Monte Carlo Simulation of near-wall turbulent
flows. Journal of Fluid Mechanics 357 (1997), S. 141–166
146 LITERATURVERZEICHNIS
[38] Eckbreth, A.C.: Laser Diagnostics for Combustion Temperature and Species, 2nd
edition. Gordon & Breach, 1996
[39] Ehrly, M.: Numerische Simulation einer nicht-vorgemischten Flamme mittles eines
hybriden CFD/transported PDF Verfahrens, Institut für technische Thermodynamik,
Universität Karlruhe, Diplomarbeit, 2009
[40] Escudier, M.: Vortex breakdown: Observations and explanations. Progress in Ae-
rospace Sciences 25 (1978), S. 189–229
[41] Favre, A.J.: The equations of compressible turbulent gases. USAF Contract AF61
(052)-772, AD 622097 (1965)
[42] Ferziger, J.H. ; Peric, M.: Computational Methods for Fluid Dynamics. 2. Aufl.
Springer Verlag, 1997
[43] Fox, R.O.: Computational Models for Turbulent Reacting Flows. Cambridge : Cam-
bridge University Press, 2003
[48] Giovangigli, V.: Convergent iterative methods for multicomponent diffusion. Im-
pact Comput. Sci. Eng 3 (1991), S. 244–276
[49] Goldstein, V. ; Sobolev, V.: Integral manifolds in chemical kinetics and combusti-
on. In: Gindikin, S.G. (Hrsg.): Singularity Theory and Some Problems of Functional
Analysis. Providence, RI : American Mathematical Society, 1992, S. 73–92
[50] Golub, G.H. ; Loan, C.F. van: Matrix Compuations. Baltimore : The Hopkins
University Press, 1989
LITERATURVERZEICHNIS 147
[52] Gowlett, J.A.J.: The early settlement of northern Europe: Fire history in the
context of climate change and the social brain. Comptes Rendus - Palevol 5 (2006),
Nr. 1-2, S. 299–310
[53] Guin, C.: Characterization of autoignition and flashback in premixed injection sy-
stems. In: AVT Symposium on Gas Turbine and Engine Combustion, Emissions and
Alternative Fuels TP1998-227 ONERA, 1998
[54] Hadjinicolaou, M. ; Goussis, D.A.: Asymptotic solution of stiff PDEs with the
CSP method: The reaction diffusion equation. SIAM journal of Scientific Computing
20 (1999), S. 781–910
[55] Hall, M.G.: Vortex breakdown. Annual Review of Fluid Mechanics 4 (1972), S.
195–218
[56] Hawkes, E.R. ; Sankaran, R. ; Sutherland, J.C. ; Chen, J.H.: Direct numerical
simulation of turbulent combustion: fundamental insights towards predictive models.
Journal of Physics: Conference Series 16 (2005), S. 65–79
[57] Haworth, D.C.: Progress in probability density function methods for turbulent
reacting flows. Progress in Energy and Combustion Science 36 (2010), Nr. 2, S. 168
– 259
[58] Haworth, D.C. ; Pope, S.B.: A Generalized Langevin Model for Turbulent Flows.
Physics of Fluids 29 (1986), S. 387–405
[60] IPCC: Climate Change 2007: Synthesis Report / Intergovernmental Panel on Climate
Change. 2007. – Forschungsbericht
[61] Janicka, J. ; Kolbe, W. ; Kollmann, W.: Closure of the Transport Equation of the
Probability Density Function of Turbulent Scalar Flieds. Journal of Non-Equilibrium
Thermodynamics 4 (1979), S. 47–66
[64] Jones, W.P. ; Launder, B.E.: The prediction of laminarization with a two-equation
model of turbulence. International Journal of Heat and Mass Transfer 15 (1972), S.
301–314
[65] Jones, W.P. ; Whitelaw, J.H.: Modelling and measurement in turbulent combusti-
on. In: Twenth Symposium (International) on Combustion The Combustion Institute,
1985, S. 233–249
[70] Kröner, M.: Einfluss lokaler Löschvorgänge auf den Flammenrückschlag durch ver-
brennungsinduziertes Wirbelaufplatzen, Technische Universität München, Fakultät für
Maschinenwesen, Dissertation, 2003
[71] Langevin, P.: On the Theory of Brownian Motion. C. R. Acad. Sci. (Paris) 146
(1908), S. 530–533
[72] Launder, B.E. ; Reece, G.J. ; Rodi, W.: Pogress in the development of a Reynolds-
stress turbulence closure. Journal of Fluid Mechanics 68 (1975), Nr. 3, S. 537–566
[73] Law, C.K.: Combustion at a crossroads: Status and prospects. In: Proceedings of the
Combustion Institute Bd. 31, 2007, S. 1–29
[74] Leibovich, S.: The structure of vortex breakdown. Annual Review of Fluid Mecha-
nics 10 (1978), S. 221–246
LITERATURVERZEICHNIS 149
[75] Leuckel, W. ; Zarzalis, N.: Theorie turbulenter StrÃűmungen ohne und mit Ãij-
berlagerter Verbrennung, Skriptum zur Vorlesung. Engler-Bunte-Institut Bereich Ver-
brennungstechnik, Uni Karlsruhe, 2003
[76] Lewis, B. ; Elbe, G. von: Combustion, Flames and Explosions of Gases. New York,
Academic, 1987
[77] Libby, P.A. ; Williams, F.A.: Turbulent Reacting Flows. Academic Press, 1994
[78] Libby, Paul A.: Introduction to turbulence. Washington, DC [u.a.] : Taylor and
Francis, 1996 (Combustion)
[79] Lindemann, F.A.: Discussion on „The radiation theory of chemical action “. Tran-
sactions of the Faraday Society 17 (1922), S. 598–606
[80] Lipp, S. ; Magagnato, F. ; Maas, U.: A hybrid transported PDF/CFD model for
turbulent flames using REDIM. In: Proceedings of the European Combustion Meeting,
Vienna, Austria, 2009
[81] Liu, B.J.D. ; Pope, S.B: The performance of in situ adaptive tabulation in com-
putations of turbulent flames. Combustion Theory and Modelling 9 (2005), Nr. 4, S.
549–568
[82] Liu, K. ; Pope, S.B. ; Caughey, D.A.: Calculations of Bluff-Body Stabilized Flames
Using a Joint Probability Density Function Model with Detailed Chemistry. Combu-
stion and Flame 141 (2005), S. 89–117
[83] Lopez, J.M.: Axisymmetric vortex breakdown part 1: confined swirling flow. Journal
of Fluid Mechanics 221 (1990), S. 533–552
[84] Lovett, J.A. ; Abuaf, N.: Emissions and stability characteristics of flameholders
for lean premixed combustion. In: Proc. International Gas Turbine and Aeroengine
Congress, 1992
[85] Lumley, J.: Computational modeling of turbulent flows. Advances in Applied Me-
chanics 18 (1978), S. 124–174
[87] Maas, U.: Chemie der Verbrennung. Vorlesungsmanuskript, ITT, Universität Karls-
ruhe (TH), 2007
150 LITERATURVERZEICHNIS
[91] marzougui, H. ; Khlifi, H. ; Lili, T.: The extension of the Launder, Reece and
Rodi model on compressible shear flow. The European physical journal B, Condensed
matter physics 45 (2005), Nr. 1, S. 147–154
[92] Masri, A.R. ; Bilger, R.W.: Conditional Probability Density Functions Measured
in Turbulent Nonpremixed Flames of Methane Near Extinction. Combustion and
Flame 74 (1988), S. 267–284
[93] Masri, A.R. ; Bilger, R.W.: Turbulent non-premixed flames of hydrocarbon fuels
near extinction: Mean structure from probe measurements. In: Twenty-First Sym-
posium (International) on Combustion The Combustion Institute, Pittsburgh, PA,
1988, S. 1511–1520
[94] Masri, A.R. ; Bilger, R.W.: Turbulent Nonpremixed Flames of Methane Near
Extinction: Probability Density Functions. Combustion and Flame 73 (1988), S. 261–
285
[95] Masri, A.R. ; Bilger, R.W. ; Dibble, R.W.: Turbulent Nonpremixed Flames of
Methane Near Extinction: Mean Stucture from Raman Measurements. Combustion
and Flame 71 (1988), S. 245–266
[96] Masri, A.R. ; Pope, S.B.: PDF Calculations of Piloted Turbulent Nonpremixed
Flames of Methane. Combustion and Flame 81 (1990), S. 13–29
[98] Merci, B. ; Roekaerts, D. ; Naud, B.: Study of the perfomance of three micro-
mixing models in transported scalar PDF simulations of a piloted jet diffusion flame
(“Delft Flame III”). Combustion and Flame 144 (2006), S. 476–493
[100] Meyer, D.W. ; Jenny, P.: A mixing model for turbulent flows based on paramete-
rized scalar profiles. Physics of Fluids 18 (2006), S. 035105–1–15
[101] Monin, A. ; Yaglom, A.: Statistical Fluid Dynamics. MIT Press, Cambridge, 1971
[103] Muradoglu, M. ; Pope, S.B. ; Caughey, D.A.: The Hybid Method for the PDF
Equations of Turbulent Reactive Flows: Consistency Conditions and Correction Al-
gorithms. Journal of Computational Physics 172 (2001), S. 841–878
[104] Nau, M.: Berechnung turbulenter Diffusionsflammen mit Hilfe eines Verfahrens zur
Bestimmung der Wahrscheinlichkeitsdichtefunktion und automatisch reduzierter Re-
aktionsmechanismen, Universität Stuttgart, Fakultät für Energietechnik, Dissertati-
on, 1997
[107] Nooren, P.A. ; Wouters, H.A. ; Peeters, T.W.J. ; Roekaerts, D. ; Maas, U.:
Monte Carlo PDF modelling of a turbulent natural-gas diffusion flame. Combustion
Theory and Modelling 1 (1997), S. 79–96
[109] Oertel, H.: Prandtl-Führer durch die Strömungslehre. Grundlagen und Phänomene.
Braunschweig : Vieweg Verlag, 2002
[110] Oijen J.A. van ; Goey L.P.H. de: Modelling of premixed laminar flames using
flamelet- generated manifolds. Compustion Science and Technology 161 (2000), S.
113–137
[111] Oijen J.A. van ; Goey L.P.H. de: Modelling of premixed counterflow flames using
the flamelet-generated manifold method. Combustion Theory and Modelling 6 (2002),
S. 463–478
152 LITERATURVERZEICHNIS
[112] Orbegoso, E.M. ; Silva, L.F.F. da: Study of stochastic mixing models for combu-
stion in turbulent flows. Proceedings of the Combustion Institute 32 (2009), Nr. 1, S.
1595 – 1603
[114] Patel, V.C. ; Rodi, W. ; Scheuerer, G.: Turbulence Models for Near-Wall and
Low Reynolds Number Flows: A Review. AIAA Journal 23 (1985), S. 1308–1319
[118] Peters, N. ; Roggs, B.: Reduced Kinetic Mechanisms for Applications in Combu-
stion Systems. In: Lecture Notes in Physics. Berlin : Springer-Verlag, 1992
[120] Pitz, R.W. ; Daily, J.W.: Cmbustion in a turbulent mixing layer frmed at a rearward
facing step. AIAA Journal 21 (1983), S. 1565–1570
[122] Pope, S.B: A Monte Carlo Method for PDF Equations of Turbulent Reactive Flow.
Combustion, Science and Technology 25 (1981), S. 159–174
[123] Pope, S.B: An Improved Turbulent Mixing Model. Combustion, Science and Tech-
nology 28 (1982), S. 131–135
[124] Pope, S.B.: PDF Methods for Turbulent Reactive Flows. Progress in Energy Com-
bustion Science 11 (1985), S. 119–192
[125] Pope, S.B.: Computations of Turbulent Combustion: Progress and Challenges. In:
Twenty-third Symposium (International) on Combustion The Combustion Institute,
Pittsburgh, PA, 1990
LITERATURVERZEICHNIS 153
[126] Pope, S.B: Lagrangian PDF Methods for Turbulent Flows. Annual Review of Fluid
Mechanics 26 (1994), S. 23–63
[127] Pope, S.B: On the relationship between stochastic Lagrangian models of turbulence
and second moment closures. Physic of Fluids A 6 (1994), S. 973
[134] Rotta, J.C.: Statistische Theorie nichthomogener Turbulenz. Zeitschrift der Physik
129 (1951), S. 547–572
[139] Spalart, P.R. ; Allmaras, S.R.: A One-Equation Turbulence Model for Aerody-
namic Flows. AIAA Journal (1992), Nr. AIAA-92-0439
154 LITERATURVERZEICHNIS
[140] Spalding, D.B.: Mixing and Chemical Reaction in Steady Confined Turbulent Fla-
mes. In: Thirteenth Symposium (International) on Combustion The Combustion In-
stitute, 1971, S. 649–657
[143] Stöllinger, M. ; Heinz, S.: Evaluation of scalar mixing and time scale models in
PDF simulations of a turbulent premixed flame. Combustion and Flame 157 (2010),
S. 1671–1685
[144] Subramaniam, S. ; Pope, S.B: Comparison of Mixing Model Performance for Non-
premixed Turbulent Reactive Flow. Combustion and Flame 117 (1999), Nr. 4, S.
732–754
[145] Subramaniam, S. ; Pope, S.B: A Mixing Model for Turbulent Reactive Flows based
on Euclidean Minimum Spanning Trees. Combustion and Flame 115 (1999), Nr. 4,
S. 487–514
[146] Szablewzki, W.: Über das Energiespektrum isotroper Turbulenz. Analen der Physik
33 (1976), Nr. 4, S. 300–308
[147] Szegö, G.G. ; Dally, B.B. ; Nathan, G.J.: Scaling of NOx emissions from a
laboratory-scale mild combustion furnace. Combustion and Flame 154 (2008), Nr.
1-2, S. 281 – 295
[150] Tennekes, H. ; Lumley, J.L.: A first course in turbulence. 1. Aufl. MIT Press,
1972. – 300 S.
[151] Tolpadi, A.K. ; Hu, I.Z. ; Correa, S.M. ; Burrus, D.L.: Coupled Lagrangian
Monte Carlo PDF-CFD Computation of Gas Turbine Combustor Flowfields with
LITERATURVERZEICHNIS 155
Finite-Rate Chemistry. Journal of Engineering for Gas Turbines and Power 119
(1997), S. 519–526
[152] Valino, L. ; Dopazo, C.: A binomial langevin model for turbulent mixing. Physics
Fluids A 3 (1991), S. 3034–3037
[153] Van Slooten, P.R. ; Pope, S.B: Application of PDF Modeling to Swirling and
Nonswirling Turbulent Jets. Flow Turbulence and Combustion 62 (1999), Nr. 4, S.
295–334
[160] Wassmer, D.: Validierung eines hybriden CFD/transported PDF Verfahrens mittels
einer turbulenten vorgemischten Methan/Luft Flamme, Institut für technische Ther-
modynamik, Universität Karlruhe, Studienarbeit, 2009
[163] Wörz, B.: Validierung eines hybriden CFD/transported PDF Verfahrens mittels
einer turbulenten nicht-vorgemischten Methan/Luft Flamme, Institut für technische
Thermodynamik, Universität Karlruhe, Studienarbeit, 2009
[165] Xu, J. ; Pope, S.B: Assessment of Numerical Accuracy of PDF/Monte Carlo Methods
for Turbulent Reacting Flows. Journal of Computational Phsics 152 (1999), Nr. 1, S.
192–230
[166] Zhang, H.S. ; So, R.M.C. ; Speziale, C.G. ; Lai, Y.G.: A near-wall two-equation
model for compressible turbulent flows. In: Aerospace Siences Meeting and Exhibit,
30th, Reno, NV AIAA, 1992, S. 23
[167] Zhang, Y.Z. ; Haworth, D.C.: A general mass consistency algorithm for hybrid
particle/finite-volume PDF methods. Journal of Computational Physics 194 (2004),
S. 156–193
Lebenslauf
Personalien
Geburtsdatum: 17.09.1978
Geburtsort: Heidelberg
Werdegang