Modellierung und Simulation einer nicht-vorgemischten ... · Dar uber hinaus geht mein Dank an...

Bachelor-Thesis:zur Erlangung des akademischen Grades Bachelor of Science (B.Sc.)

Modellierung und Simulation einer nicht-vorgemischtenGleichstrom-Wasserstoff-Verbrennung

Vorgelegt von: Jan Gabriel

Computational Science and Engineering

16. Juni 2015

Gutachter:

Prof. Dr. Dirk Lebiedz

M.Sc. Pascal Heiter

Betreuer:

M.Sc. Pascal Heiter

M.Sc. Hagen Muller

Danksagung

An erster Stelle mochte ich Herrn Prof. Dr. Dirk Lebiedz meinen Dank aussprechen,

als Absolvent fur die Moglichkeit, als interdisziplinarer Student eine ebenso interdis-

ziplinare Bachelor-Thesis schreiben zu durfen, und ebenso als einer seiner Studenten

fur das Aufzeigen neuer Wege und eines alternativen Studiengedankens.

Außerdem mochte ich Hagen Muller danken, fur die zahlreichen Ratschlage, ob per

eMail oder personlich, die mir meinen Weg durch thematisches Neuland wesentlich

erleichtert haben.

Daruber hinaus geht mein Dank an Pascal Heiter, sowohl fur die Muhe als mein

Zweitgutachter und Betreuer zu fungieren, als auch fur lehrreiche Momente in der

Universitat und gesellige Momente in der Freizeit.

Zuletzt mochte ich meiner Familie und meinen Freunden danken, fur jeden Mo-

ment, jedes Gesprach und die mentale Unterstutzung!

I

Zusammenfassung

Die vorliegende Bachelor-Thesis behandelt die Modellierung und Simulation einer

nicht-vorgemischten Wasserstoff-Verbrennung. In einem durch Symmetrieeigenschaf-

ten vereinfachten Stromungsfeld werden die turbuleten Stromungsverhaltnisse un-

tersucht und mit einem passenden Turbulenzmodell modelliert. Die Simulation wird

in der rein Terminal-basierten Open-Source-Software OpenFoam [1] durchgefuhrt,

in welcher die Flamme auf einem Finite-Volumen-Gitter, dessen dritte Dimensi-

on aus einem Element besteht, berechnet wird. Mit Hilfe eines neuen numerischen

OpenFoam-Losers und der Software Cantera konnen die vielfaltigen chemischen Zu-

sammenhange von der Stromung entkoppelt modelliert und anschließend mit dersel-

ben wieder zusammengefuhrt werden. Durch die Anwendung einer Parametervariati-

on des Turbulenzmodells, einer Netz-Konvergenzanalyse und durch die Aufpragung

von funktional angepassten Randbedingungen sind die Verlaufe der Simulations-

ergebnisse noch wesentlich verbessert worden. Die Validierung erfolgt durch den

Vergleich mit experimentellen Daten aus den aktuellen Forschungsergebnissen zu

nicht-vorgemischten Wasserstoff-Flammen von Robert Barlow [3].

Abstract

This present bachelor thesis deals with the modeling and simulation of non-premixed

coflow hydrogen combustion. The radial symmetric flow field is reduced to two di-

mensions and is calculated with the reynolds-averaged Navier-Stokes approach. The

simulation takes place in a terminal-based open source environment named Open-

Foam, which uses a finite volume method. A new solver for OpenFoam is introduced

to split the chemical calculations apart from the Navier-Stokes solution, to reunite

them again after solving. Parameter variation of the k-ε-model, mesh convergence

analysis methods and functional boundary adapting leads to a quite optimized si-

mulation, that can be validated by comparison to current experiments performed

by Robert Barlow [3] in the Sandia National Laboratories, Livermore, CA.

II

Inhaltsverzeichnis Inhaltsverzeichnis

Inhaltsverzeichnis

1. Einleitung 3

2. Grundlagen 5

2.1. Thermodynamik . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5

2.1.1. Hauptsatze der Thermodynamik . . . . . . . . . . . . . . . . . 5

2.1.2. Helmholtz- und Gibbs-Energie . . . . . . . . . . . . . . . . . . 7

2.2. Chemische Kinetik reaktiver Stromung . . . . . . . . . . . . . . . . . 8

2.2.1. Reaktions-Mechanismus . . . . . . . . . . . . . . . . . . . . . 9

3. Modellierung turbulenter Flammen 11

3.1. Die Flamelet Methode . . . . . . . . . . . . . . . . . . . . . . . . . . 11

3.2. Reynolds-Averaged-Navier-Stokes . . . . . . . . . . . . . . . . . . . . 14

3.3. Turbulenzmodelle . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16

4. Simulation einer Flamme 19

4.1. OpenFoam . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19

4.2. Die Gleichstrom-H2-Flamme . . . . . . . . . . . . . . . . . . . . . . . 22

4.2.1. Geometrie . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22

4.2.2. Randbedingungen . . . . . . . . . . . . . . . . . . . . . . . . . 24

4.2.3. Rechner-Architektur . . . . . . . . . . . . . . . . . . . . . . . 31

4.2.4. Netz-Konvergenzanalyse . . . . . . . . . . . . . . . . . . . . . 32

4.2.5. Zeitliche Konvergenz . . . . . . . . . . . . . . . . . . . . . . . 36

5. Auswertung 37

5.1. Mischungsbruch . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37

5.2. Massenbruch H2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39

5.3. Massenbruch O2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40

5.4. Variation der Modellkonstanten . . . . . . . . . . . . . . . . . . . . . 41

6. Fazit 43

Abbildungsverzeichnis 46

1

Inhaltsverzeichnis Inhaltsverzeichnis

A. Matlab Skript zur Berechnung der Geschwindigkeitsverteilung 48

B. Matlab Skript zur Berechnung der k-epsilon-Verteilung 50

C. Matlab Skript zur Daten-Evaluation 54

D. Netz-Konvergenz 58

E. Auswertung: H2O und Temperatur 60

2

Einleitung

1. Einleitung

Die Simulation von reaktiven Stromungen dient in erster Linie dazu, das Verstandnis

fur thermodynamische Zusammenhange selbst in komplexeren Umgebungen zu er-

weitern. In Ottomotoren beispielsweise werden die Flammen des entzundeten Treib-

stoffs simuliert um verschiedene Faktoren, wie Schadstoffausstoß oder Effizienz der

Energieumwandlung zu optimieren. Selbst nicht-beobachtbare Phanomene in der

Realitat konnen in Simulationen untersucht und ausgewertet werden.

In dieser Bachelorarbeit soll eine Grundlage fur weitere Simulationen geschaffen wer-

den, in denen die Losungen eines kompletten chemischen Systems einem reduzierten

System gegenubergestellt werden sollen. Eine reduzierte Wasserstoff-Verbrennung

ist zudem nur der erste Schritt nach dieser Bachelorarbeit, welche auf einem un-

reduzierten chemischen System basiert. Als nachstes werden daraufhin reduzier-

te Methan(CH4)-Verbrennungen und letztendlich reduzierte Kerosin-Verbrennungen

folgen.

Abbildung 1: DLR-Flamme (m.) mit Darstellung durch Rayleigh-Scattering (l.) undPLIF (Planar Laser Induced Fluorescence, r.). [14]

3

Einleitung

Nomenklatur

U Innere Energie [J]

W Arbeit [J]

Q Warme [J]

H Enthalpie [J]

p Druck [Pa]

V Volumen [m3]

S Entropie [J/K]

A Freie Energie/ Helmholtz-Energie [J]

G Freie Entropie/ Gibbs Energie [J/K]

Ea Aktivierungsenergie [J]

t Zeit [s]

ci Konzentration der Spezies i [-]

νri Stochiometrischer Koeffizient der Reaktion r und Spezies i [-]

Z Mischungsbruch [-]

M Molare Masse [kg/mol]

Yk Massenbruch der Spezies k [kg/kg]

Yel,i elementarer Massenbruch des Elements i [kg/kg]

T Temperatur [K]

R Universelle Gaskonstante R = 8.314 [J/molK]

h Spezifische Enthalpie [J/kg]

hk Spezifische Enthalpie der Spezies k [J/kg]

χst Stochiometrische skalare Dissipaionsrate [1/s]

ρ Dichte [kg/m3]

ωk Spezies-Reaktionsquellterm der chemischen Kinetik [kg/m3s]

cp Spezifische Warmekapazitat bei konstantem Druck [J/kgK]

νT Turbulenter Austauschkoeffizient [-]

v, u Geschwindigkeit [m/s]

r, R Radius [mm]

k Turbulente kinetische Energie [m2/s2]

ε Turbulente Dissipationsrate [m2/s3]

4

Grundlagen

2. Grundlagen

In diesem Kapitel werden die theoretischen Grundlagen fur die Simulation reaktiver

Stromung behandelt. Es werden sowohl thermodynamische, als auch reaktionskine-

tische Sachverhalte eingefuhrt

2.1. Thermodynamik

Die Grundbausteine der Thermodynamik sind in erster Linie die drei Hauptsatze

der Thermodynamik und die Definition der Helmholtz- und Gibbs-Energie.

2.1.1. Hauptsatze der Thermodynamik

1.”In einem isolierten System ist die Summe aller Energien konstant.“

Die Anderung der inneren Energie eines solchen Systems ist definiert als die

Summe der darauf ausgeubten Arbeit und der in das System uberfuhrten

Warme:

dU = δW + δQ. (2.1)

In dieser Formel steht das Differential dU fur das totale Differential der inneren

Energie, welches sich aus den partiellen Differentialen δW und δQ zusammen-

setzt. Aus ihr folgt die Tatsache, dass Energie weder erzeugt noch zerstort

werden kann, sondern nur in eine andere Art der Energie umgewandelt wer-

den kann.

Betrachtet man ein offenes System, nehmen außer der Warme und der Arbeit

noch die Massenstrome in Form von Verschiebearbeit Einfluss auf das System.

Zur Bilanzierung wird hierbei die Zustandsvariable Enthalpie verwendet und

der erste Hauptsatz lautet dann

H = U + pV, (2.2)

bzw.

dH = dU + pdV + V dp. (2.3)

5

2.1 Thermodynamik Grundlagen

2.”Thermische Energie ist nicht in beliebigem Maße in andere Energiearten um-

wandelbar.“

Dieser Satz impliziert zweierlei Aussagen: Zum einen die vorgeschriebene Rich-

tung der Warmeanderung. Das bedeutet, dass ein Prozess, in welchem Warme

ausschließlich von einem kalten auf einen warmen Korper transferiert wird,

nicht existiert. Zum anderen wird die Aussage des ersten Hauptsatzes der

Thermodynamik eingeschrankt, da Energie zwar in Form von Arbeit in Warme

vollstandig umgewandelt werden kann aber die umgekehrte Richtung hingegen

so nie vollstandig realisiert wird. An dieser Stelle bedarf es einer Einfuhrung

einer weiteren Zustandsvariable, der Entropie. Die Entropie kann auf verschie-

dene Art und Weise beschrieben werden, wie beispielsweise als”ein Maß fur

die Zahl der zuganglichen, energetisch gleichwertigen Mikrozustande [..., als]

ein Maß fur das Unwissen uber die mikroskopische Struktur “[7], oder mit

den Worten von Warnetz et. al.:”[E]ntropy is seen as a measure for mulecu-

lar chaos“[6]. Die hier relevante Eigenschaft der Entropie ist die Beziehung zu

reversiblen und irreversiblen Prozessen der Energieumwandlung. Die Entropie

wird niemals kleiner, sie kann nur entstehen, wenn beispielsweise mechani-

sche Energie durch Reibung in thermische Energie umgewandelt wird. Durch

irreversible Prozesse steigt die Anzahl an energetisch gleichwertigen Mikro-

zustanden und somit die Entropie selbst. Es gilt allgemein

dS =∂Qrev

T, dS >

∂Qirrev

T, (2.4)

bzw. im geschlossenen System (∂Q = 0)

(dS)rev = 0 , (dS)irrev > 0. (2.5)

3.”Am absoluten Nullpunkt verschwinden die Entropien aller in einem inneren

Gleichgewichtszustand befindlichen reinen Stoffe.“

Als innerer Gleichgewichtszustand wird anschaulich eine ideale Kristallgitter-

Anordnung der Atome bezeichnet. Am absoluten Nullpunkt hat ein solches

Gitter keine Moglichkeit mehr zu schwingen, (ir-)reversible Prozesse konnen

nicht mehr ablaufen, und damit kann es keine Anderung der Entropie mehr

geben. Da der Grenzwert der Temperatur nur asymptotisch erreicht wird, ist es

6

2.1 Thermodynamik Grundlagen

auch nicht moglich die entsprechende Entropie in endlicher Zeit zu erreichen.

limT→0

S = 0. (2.6)

Um (dennoch) einen Referenzpunkt fur Berechnungen zu haben, setzte Max

Planck 1912 den absoluten Nullpunkt der Entropie auf Null [8].

2.1.2. Helmholtz- und Gibbs-Energie

Uber das thermodynamische Gleichgewicht hinaus interessiert zudem das chemi-

sche Gleichgewicht fur praktische Anwendungen. Um die Thermodynamik adaquat

beschreiben zu konnen, mussen zwei weitere extensive Zustandsgroßen eingefuhrt

werden: die Helmholtz- und die Gibbs- Energie.

Die Helmholtz-Energie ist eine Kombination der ersten beiden Hauptsatze der Ther-

modynamik:

dU + pdV − TdS 6 0, (2.7)

im chemischen Gleichgewicht:

dU + pdV − TdS = 0, (2.8)

(dU)V,S = 0. (2.9)

Der Index (V, S) bedeutet ausgeschrieben V = const., S = const.. Aufgrund der

schlechten Messbarkeit der Entropie in der Praxis ist die Forderung S = const.

nur schwer zu erfullen. Daher ersetzt man mithilfe der Kettenregel der Differenti-

alrechnung TdS durch d(TS) − SdT , wodurch sich Gleichung 2.8 umformen lasst

zu

d(U − TS) + pdV + SdT = 0, (2.10)

d(A)V,T = 0, (2.11)

mit A = U − TS, der sogenannten freie Energie oder Helmholtz-Energie. Sie druckt

die Energie aus, die benotigt wird um ein System in ein chemisches Gleichgewicht zu

bringen. Die Gibbs-Energie kann durch eine weitere Anwendung der Kettenregel der

Differentialrechnung auf Gleichung 2.10 mit pdV = d(pV )− V dp formuliert werden

7

2.2 Chemische Kinetik reaktiver Stromung Grundlagen

als

d(U − TS + pV )− V dp+ SdT = 0, (2.12)

mitG = U−TS+pV , der sogenannten freien Enthalpie oder Gibbs-Energie. Sie dient

als Hilfsmittel zur Darstellung des reaktiven Potentials einer Reaktionsgleichung und

deren Energieanderung entlang einer bestimmten Reaktionsrichtung.

Im Fall

dG < 0, (2.13)

lauft die Reaktion direktional exergon, also unter Energieausschuttung ab. Wenn

dG = 0, (2.14)

steht der Reaktionsprozess im energetischen Gleichgewicht, was in Gleichung 2.12

der Fall ist. Die Reaktion verlauft direktional endergon, falls

dG > 0. (2.15)

Sie lauft also unter Energiezufuhr ab.

2.2. Chemische Kinetik reaktiver Stromung

In der chemischen Kinetik werden aus Reaktionsgleichungen mathematische For-

mulierungen abgeleitet und als Differentialgleichungssystem dargestellt. Dies ge-

schieht indem die Anderung der Konzentration einer Spezies in Abhangigkeit zu

den involvierten Spezies gesetzt wird. Die Geschwindigkeit einer solchen Konzentra-

tionsanderung wird durch die Ratengleichung nach Arrhenius beschrieben

k = A′ · exp

(− EaRT

), (2.16)

wobei A′ ein praexponentieller Faktor ist, welcher je nach Molekularitat in einer an-

deren Art und Weise von der Temperatur abhangt. Ea ist die Aktivierungsenergie,

R die universelle Gaskonstante und T die absolute Temperatur in [K].

Chemische Reaktionsgleichungen werden in Elementarreaktionen und Globale Re-

aktionen unterschieden, wobei globale Reaktionen aus einer Vielfalt von Elementar-

8


reaktionen bestehen. Die globale Reaktion

2H2 +O2k−→ 2H2O, (2.17)

besteht beispielsweise aus 38 Elementarreaktionen. Diese sind molekular exakt und

konnen somit mittels Arrhenius-Kinetik als Differentialgleichungssystem formuliert

werden.

Im Allgemeinen lauft eine Reaktion nicht ausschließlich in eine Richtung ab sondern

sie verhalt sich reversibel und pendelt sich dann in einem Gleichgewichtspunkt ein.

A+B + C + ...k1k−1

D + E + F + ... (2.18)

Annahme: Gleichgewicht zwischen Hin- und Ruckrichtung:

d[A]

dt

(h)

=d[A]

dt

(r)

, (2.19)

⇔ k1 · [A]a[B]b[C]c · ... = k−1 · [D]d[E]e[F ]f · ... , (2.20)

⇔ k1k−1

=[A]a[B]b[C]c · ...[D]d[E]e[F ]f · ...

= Kc = exp(−∆RA0

RT). (2.21)

Die sogenannte Gleichgewichtskonstante Kc beinhaltet sowohl Informationen uber

die Beziehung der Konzentrationen der Reaktion im Gleichgewichtspunkt, als auch

uber die resultierende Differenz der freien Helmholtz-Energie ∆RA0.

2.2.1. Reaktions-Mechanismus

Der Mechanismus zur mathematischen Formulierung einer chemischen Reaktion lau-

tet wie folgt: (∂ci∂t

)chem,r

= kr

(ν(p)ri − ν

(e)ri

) S∏s=1

cν(e)rss . (2.22)

Dieser beschreibt die Anderung der Konzentration c einer Spezies i in der entspre-

chenden Reaktion r. Die Arrhenius-Konstante kr gibt in einer reversiblen Reaktion

stets die Geschwindigkeit einer Richtung an und wird aus Gleichung (2.16) errech-

net. Die Arrhenius-Konstante der entsprechend anderen Richtung einer reversiblen

Reaktion muss jedoch uber die Gleichgewichtskonstante Kc aus Gleichung (2.21)

9


bestimmt werden. Die Koeffizienten ν(p)ri und ν

(e)ri bilanzieren die stochiometrische

Anderung der Spezies i wobei sich ν(e)rs im Exponent von cs auf alle Spezies der

Edukt-Seite der Reaktion bezieht.

In einem System von Elementarreaktionen kommen viele Spezies mehrfach vor. Da-

her wird nun uber alle Reaktionen r welche die Spezies i enthalten aufsummiert um

das System vollstandig zu beschreiben:

(dcidt

)chem

=R∑r=1

kr

(ν(p)ri − ν

(e)ri

) S∏s=1

cν(e)rss . (2.23)

Beispiel, aus [6]:

Dieses Reaktionsgleichungssystem zeigt die Verbrennungsreaktion von Wasserstoff

(H2) mit Chlor(Cl2):

Cl2 +Mk1k4Cl + Cl +M,

Cl +H2k2−→ HCl +H,

H + Cl2k3−→ HCl + Cl.

Das Differentialgleichungssystem (2.24) - (2.28) ist anhand des oben eingefuhrten

Mechanismus erstellt worden, wobei der Reaktant M eine besondere Rolle spielt:

In einer unimolekularen Reaktion werden in der Natur stochiometrisch unbeteiligte

Reaktanten hinzugezogen um beispielsweise bei der Spaltung von Cl2 zu Cl+Cl die

freigewordene Energie aufzunehmen. Wird diese Energie nicht aufgenommen kommt

es mit dem extremen Energieuberschuss instantan zu einer Ruckreaktion.

d[Cl]

dt= 2k1[Cl2][M ]− k2[Cl][H2] + k3[H][Cl2]− 2k4[Cl]2[M ], (2.24)

d[Cl2]

dt= −k1[Cl2][M ]− k3[H][Cl2] + k4[Cl]2[M ], (2.25)

d[H]

dt= k2[Cl][H2]− k3[H][Cl2], (2.26)

d[H2]

dt= −k2[Cl][H2], (2.27)

d[HCl]

dt= k2[Cl][H2] + k3[H][Cl2]. (2.28)

10

Modellierung turbulenter Flammen

3. Modellierung turbulenter Flammen

Die Modellierung von Turbulenzen wird in der aktuellen CFD-Simulation mit eini-

gen Methoden durchgefuhrt, wie beispielsweise RANS (Reynolds-Averaged Navier-

Stokes, siehe Kapitel 3.2) oder LES (Large-Eddy Simulation). Eine direkte nume-

rische Simulation wird besonders bei uberkritischen Reynolds-Zahlen sehr rechen-

aufwendig und selbst Hochleistungs-Cluster gelangen bei einer entsprechend feinen

Vernetzung an ihre Grenzen. Die Modellierung turbulenter Flammen wird dahinge-

hend noch verkompliziert, dass auf jedem Knoten simultan die komplette chemische

Kinetik berechnet werden muss. Eine Moglichkeit das Problem des Rechenaufwands

zu losen wird im Folgenden erlautert. Die Methode zur Modellierung tubulenter

Flammen mit dem Flamelet-Ansatz und die Implementierung eines numerischen

Losers zur Berechnung dieser Modellierung in der Opensource-Software OpenFoam

beruht auf der Arbeit von Muller, Ferraro und Pfitzner [2].

3.1. Die Flamelet Methode

Die Quintessenz der Flamelet-Methode ist die Entkopplung der Stromungsberech-

nung mit RANS/LES von der Losung des Differentialgleichungssystems der chemi-

schen Kinetik. Hierbei wird die turbulente Flamme in laminare”Flammchen“ (Fla-

melets) unterteilt, d.h. auf jedem Knoten wird einzeln eine laminare Gegenstrom-

Diffusionsflamme (vgl. Abbildung 2) berechnet, welche die Turbulenz in dem jewei-

ligen Knoten approximiert.

Abbildung 2: Laminare Gegenstrom-Diffusionsflamme, [2].

11

3.1 Die Flamelet Methode Modellierung turbulenter Flammen

Jedes Flamelet wird in einer Tabelle in der Flamelet-Bibliothek gespeichert. Um

diese Tabellen zu erstellen mussen zunachst die Flamelet-Gleichungen

ρ∂Yk∂t− 1

2ρχst

∂2Yk∂Z2

− ωk = 0, (3.1)

ρ∂T

∂t− 1

2ρχst

(∂2T

∂Z2+

1

cp

∂cp∂Z

∂T

∂Z

)+

1

cp

∑hkωk = 0, (3.2)

gelost werden. Diese liefern die Funktionen T (Z, χst), Yk(Z, χst) welche jeweils von

dem Mischungsbruch Z und der stochiometrischen skalaren Dissipationsrate χst

abhangen. Anschließend werden T und Yk mit einer Wahrscheinlichkeits-Dichte-

Funktion (PDF) uber den kompletten Raum integriert um die Massenbruche der

einzelnen Spezies und die Temperatur zu Mitteln, mit

Yk =

∞∫0

1∫0

Yk (Z, χst)P (Z, χst) dZdχst, (3.3)

und

T =

∞∫0

1∫0

T (Z, χst)P (Z, χst) dZdχst. (3.4)

Die Wahrscheinlichkeitsdichte P (Z, χst) kann nach dem Satz der stochastischen

Unabhangigkeit nach de Moivre als Produkt marginaler Dichtefunktionen

P (Z, χst) = P (Z) · P (χst) , (3.5)

geschrieben werden.

Hierdurch konnen beide Variablen auf verschiedene Weise verteilt werden. Die Ver-

teilung der stochiometrischen skalaren Dissipationsrate wird durch eine Dirac-Funk-

tion P (χst) = δ(χst − χst) beschrieben, wahrend der Mischungsbruch β -verteilt

ist:

P (Z) = Zα−1(1− Z)β−1Γ(α + β)

Γ(α)Γ(β), (3.6)

wobei Γ(x) die Gamma-Funktion ist, mit geeigneten formdefinierenden Parametern

α und β.

Aus Kompatibilitatsgrunden mit der verwendeten Software OpenFoam (siehe Ka-

12

3.1 Die Flamelet Methode Modellierung turbulenter Flammen

pitel 4.1) wird aus der Temperatur T die spezifische Enthalpie h berechnet, da aus

dieser samtliche thermodynamischen Großen abgeleitet werden konnen. Die resultie-

renden Gleichungen h(Z, χst, Z′′2) und Yk(Z, χst, Z

′′2), welche jeweils vom gemittel-

ten Mischungsbruch, dessen Varianz und der gemittleten skalaren Dissipationsrate

abhangen, werden in der mehrdimensionalen Flamelet-Library zwischengespeichert,

wo h und Yk zur Laufzeit des Solvers ausgelesen und gegebenenfalls interpoliert

werden. Fur die Stromungsberechnung wurde die spezielle OpenFoam Applikation

flameletFoam entwickelt, welche die Navier-Stokes-Gleichungen unabhangig von der

Chemie berechnet, anschließend anhand von Parametern die Flamelet-Library aus-

liest und jedem Knoten die jeweiligen thermodynamischen Großen zuordnet. Dieses

Prinzip wird in Abbilung 3 nochmals verdeutlicht.

Abbildung 3: Funktionalitat der flameletFoam-Applikation, [2].

13

3.2 Reynolds-Averaged-Navier-Stokes Modellierung turbulenter Flammen

3.2. Reynolds-Averaged-Navier-Stokes

Das Reynolds-Averaged-Navier-Stokes (RANS) Modell beschreibt die Navier-Stokes-

Gleichungen einer turbulenten Stromung auf der Basis gemittelter Zeit- und Dich-

tekomponenten (Favre-Average). Bei der Zeitmittlung wird eine Variable q in ihren

zeitlichen Mittelwert q und ihre Fluktuation q′ aufgeteilt. Hierbei ist wichtig zu

beachten, dass der Mittelwert der Fluktuationen Null sein muss.

q = q + q′, q′ = 0. (3.7)

Die Dichtemittlung (Favre-Average) wird bei großeren Dichteschwankungen, wie sie

in Verbrennungsprozessen vorkommen, verwendet. Sie verhalt sich ahnlich wie die

zeitliche Mittlung, bis auf den Unterschied, dass der Mittelwert mit der Dichte ρ

gewichtet wird

q =ρq

ρ, (3.8)

und die Variable q erneut als ihr (gewichteter) Mittelwert q und ihre Fluktuation q′′

ausgedruckt werden kann.

q = q + q′′, q′′ = 0. (3.9)

Die direkte numerische Simulation wird durch gemittelte Werte dahingehend ver-

bessert, dass die Menge der Daten auf ein Minimum reduziert wird wohingegen sich

die Aussagekraft der Daten relativ gesehen verbessert. Im Gegensatz zur direkten

Simulation bzw. einer LES-Simulation (Large-Eddy-Simulation) werden hier als Re-

sultat der Mittlung die Wirbelstrukturen nicht mitberechnet (vgl. Abbildung 4).

Um diese Gleichungen numerisch losen zu konnen, werden in einigen verbreiteten

Turbulenzmodellen die fluktuierenden Terme der Form ρ~v′′q′′ als turbulenter Trans-

port interpretiert und mit dem Gradienten des Mittelwerts jeder Große angenahert

ρ~v′′q′′ = −ρνT grad qi. (3.10)

Verschiedene Modelle des turbulenten Transports werden im folgenden Kapitel 3.3

durch Variation des turbulenten Austauschkoeffizienten νT eingefuhrt und erlautert.

14

3.2 Reynolds-Averaged-Navier-Stokes Modellierung turbulenter Flammen

Abbildung 4: Schematischer Vergleich einer einer direkten Simulation (a) und einerRANS-Simulation (b), [9].

Die RANS-Gleichungen (nach Warnetz et. al. [6]) setzen sich aus den folgenden

Gesetzen zusammen:

• Massenerhaltung/Kontinuitat

∂ρ

∂t+ div

(ρ~v)

= 0, (3.11)

• Massenerhaltung der einzelnen Spezies

∂(ρYi

)∂t

+ div(ρ~vYi

)+ div

(−ρνT grad Yi

)= Miωi, (3.12)

• Impulserhaltung

∂ (ρvi)

∂t+ div

(ρ~v ⊗ ~v

)+ div

(−ρνT grad ~v

)= ρ~g, (3.13)

• Energieerhaltung

∂(ρh)

∂t− ∂p

∂t+ div

(ρ~vh)

+ div(−ρνT grad h

)= qr, (3.14)

• Elementare Massenerhaltung

∂(ρYel,i

)∂t

+ div(ρ~vYel,i

)+ div

(−ρνT grad Yel,i

)= 0. (3.15)

15

3.3 Turbulenzmodelle Modellierung turbulenter Flammen

Die Gleichung (3.15) beschreibt zusatzlich zu den vier RANS-Erhaltungsgleichungen

die Massenerhaltung der beteiligten Elemente. Der Elementare Massenbruch wird

hier definiert als

Yel,i =∑

µijYi, (3.16)

mit dem Verhaltnis µij zwischen dem Element i und der Spezies j. Aufgrund der

problematischen Handhabung der Reaktionsquellterme Miωi in Gleichung 3.12 kann

auf die Element-Erhaltung ∑µijMiωi = 0, (3.17)

zuruckgegriffen werden. Diese entspricht der ungemittelten rechten Seite der Spezies-

Massenerhaltungsgleichung, welche mittels Multiplikation der gesamten Gleichung

mit µij und Aufsummation uber alle beteiligten Elemente die Gleichung 3.15 ergibt.

Diese muss naturlich ebenso zeitlich gemittelt werden.

3.3. Turbulenzmodelle

• Null-dimensionale Transportgleichungen

In dieser Art der Modellierung wird der turbulente Transport und der tur-

bulente Austauschkoeffizient durch mathematische Terme substituiert ohne

eine zusatzliche Differentialgleichung aufzusetzen. Beispielsweise mit der For-

mel der Mischungslange von Prandtl (1925) [6]. Hierbei wird der Transport

modelliert als

ρ~v′′q′′ = −ρ`2∣∣∣∣∂v∂z

∣∣∣∣ ∂q∂z und νT = `2∂v

∂z. (3.18)

Am Beispiel in Abbilung 5 ist ersichtlich, dass der turbulete Transport einer

Abbildung 5: Beispiel einer Scherstromung mit turbulenter Grenzschicht,[6].

16


Scherstromung von der Anderung der Geschwindigkeit in Scherrichtung (z-

Richtung) und einer problemspezifischen Lange ` , welche von der Scherschicht

δ abhangt. Diese Modellierung ist jedoch heutzutage uberholt.

• Eindimensionale Transportgleichungen

Der turbulente Transport wird, im Gegensatz zur null-dimensionalen Modellie-

rung, mit einer Differentialgleichung beschrieben. Der turbulente Austausch-

koeffizient νT berechnet sich mit der turbulenten kinetischen Energie

k =1

2

ρ∑v′2i

ρ, (3.19)

aus der Formel nach Prandtl (1945), [6]

νT = `√k, (3.20)

wobei sich k als die Losung der Differentialgleichung ergibt und v′2i die Diago-

nalelemente des Reynold’schen Spannungstensors bezeichnet (vgl. Gleichung (4.7)).

Die Mischungslange ` wird analog zur null-dimensionalen Modellierung errech-

net. Diese Modellierung ist heutzutage ebenso uberholt.

• Zweidimensionale Transportgleichungen

Die wohl popularste Transport-Modellierung ist das k-ε-Turbulenzmodell als

Spezialfall zweidimensionaler Transportgleichungen, aufgrund der zusatzlichen

Berechenbarkeit von turbulenten Auswirkungen, wie der Konvektion und der

Diffusion [16].

Die Formulierung des turbulenten Austauschkoeffizienten lautet

νt = Cµk

ε, (3.21)

aus der allgemeinen Formel

νT ∝ z1n k

12−m

n , mit n = −1 und m =3

2, (3.22)

wobei im Gegensatz zu anderen Turbulenzmodellen wie beispielsweise dem

k-ω-Turbulenzmodell die Varieble z der gemittelten turbulenten Dissipations-

17


rate ε entspricht. Das zugehorige Differentialgleichungssystem, welches zusatz-

lich zu den RANS-Gleichungen gelost wird lautet

∂(ρk)

∂t+ div(ρ~vk)− div(ρνT grad k) = Gk − ρε, (3.23)

∂(ρε)

∂t+ div(ρ~vε)− div(ρνT grad ε) = (C1Gk − C2ρε)

ε

k. (3.24)

Die Gleichung (3.24) basiert auf reiner Empirie, ebenso wie die Modellkon-

stanten Cµ = 0.09, C1 = 1.44, C2 = 1.92 und der Konstanten Gk welche sich

aus dem Term

Gk = −ρ~v′′ ⊗ ~v′′ : grad ~v, (3.25)

ergibt. Auf deren Problemabhangigkeit - speziell der Konstanten C1 - wird in

Kapitel 5.4 noch weiter eingegangen und die Anfangswerte des Differentail-

gleichungssystems werden in Kapitel 4.2.2 thematisiert.

18

Simulation einer Flamme

4. Simulation einer Flamme

4.1. OpenFoam

Die quelloffene Software OpenFoam (Open Source Field Operation and Manipulati-

on, [1]) bietet eine Vielfalt an Anwendungen zur Losung laminarer und turbulenter

Stromung, sowohl im kommerziellen als auch in akademischen Bereich. Sie beinhal-

tet einige verbreitete Turbulenzmodelle, wie beispielsweise RAS (Reynolds-Averaged

Simulation), LES (Large Eddy Simulation), DES (Detached Eddy Simulation) oder

DNS (Direct Numerical Simulation) zur numerischen Stromungssimulation, speziell

auch fur reaktive Stromung. Der strukturelle Aufbau eines OpenFoam-Cases ist in

Abbildung 6 ersichtlich. Der system-Ordner beinhaltet samtliche Informationen zur

Abbildung 6: Ordner-Struktur eines OpenFoam-Cases, [15].

numerischen Berechnung und Einstellungen fur die Architektur auf der gerechnet

werden soll. In den Dateien fvSchemes und fvSolution sind alle Angaben zu numeri-

schen Losern gespeichert und in der controlDict-Datei werden, wie im nachfolgenden

Code-Beispiel gezeigt, die verwendete Applikation, die Zeitspanne und einige weitere

Parameter angegeben.

19

4.1 OpenFoam Simulation einer Flamme

1 /*--------------------------------*- C++ -*----------------------------------*\

2 | ========= | |

3 | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |

4 | \\ / O peration | Version: 2.3.1 |

5 | \\ / A nd | Web: www.OpenFOAM.org |

6 | \\/ M anipulation | |

7 \*---------------------------------------------------------------------------*/

8 FoamFile

9 {

10 Version 2.0;

11 Format ascii;

12 class dictionary;

13 location "system ";

14 object controlDict;

15 }

16

17 application flameletFoam;

18

19 startFrom startTime;

20

21 startTime 0;

22

23 stopAt endTime;

24

25 endTime 0.65;

26

27 deltaT 2.5e-7; // initial time step

28

29 writeControl adjustableRunTime;

30

31 writeInterval 1e-2;

32

33 purgeWrite 0;

34

35 writeFormat ascii;

36

37 writePrecision 6;

38

39 writeCompression off;

40

41 timeFormat general;

42

43 timePrecision 6;

44

45 runTimeModifiable true;

46

47 adjustTimeStep yes;

48

49 maxCo 0.3;

20

4.1 OpenFoam Simulation einer Flamme

Hierbei ist wichtig zu erwahnen, dass der Parameter maxCo den Wert 0.3 nicht

uberschreiten darf. Die sogenannte Courant-Zahl Co steht in einer Stromungs-Simu-

lation fur das Verhaltnis zwischen einem Zeitschritt ∂t pro Zellengroße ∂x, multipli-

ziert mit der Mittleren Geschwindigkeit in der jeweiligen Zelle des Gitters, also

Co =∂t|U |∂x

. (4.1)

Ist adjustT imeStep = yes, wird der Zeitschritt ∂t automatisch angepasst um die

Courant-Zahl stets unter 0.3 zu halten. Wird die Courant-Zahl zu groß (kritisch bei

Co > 1), wird in einem Zeitschritt eine ganze Zelle im Netz ubersprungen, was ein

feines Netz unnutz macht.

Im system-Ordner befindet sich in diesem Fall noch eine decomposeParDict-Datei,

mit welcher die Moglichkeit besteht, einen OpenFoam-Case parallel via shared me-

mory oder auch via distributed memory zu berechnen.

Im constant-Ordner des OpenFoam-Cases befindet sich der Ordner polyMesh und

zusatzlich Dateien zur Spezifizierung des Turbulenz-Modells, der Art der Verbren-

nung, der Auswahl der chemischen Tabellen und der thermophysikalischen Großen.

Die Dateien im polyMesh-Ordner sind ein Produkt der Verarbeitung der Datei block-

MeshDict durch die OpenFoam-Applikation blockMesh. Sie beinhalten alle relevan-

ten Informationen zum Aufbau und der Art des Rechengitters.

Will man nun einen OpenFoam-Case berechnen, bedient man sich aus Aufwands-

grunden eines make-files. Das Intervall, in welchem die Losungen in Ordnern gespei-

1 # make -file: Allrun

2

3 . $WM_PROJECT_DIR/bin/tools/RunFunctions #Source tutorial run functions

4

5 runApplication blockMesh #create Mesh

6

7 runApplication canteraToFoam #integrate tables

8

9 runApplication decomposePar #decompose case

10

11 runParallel flameletFoam 12 #start solver

12

13 runApplication reconstructPar #reconstruct case

14

15 rm -rf processor*

16

17 # calculate species mass fractions:

18 runApplication flameletFoamPost -latestTime

21

4.2 Die Gleichstrom-H2-Flamme Simulation einer Flamme

chert werden sollen, wird in der controlDict-Datei via writeInterval angegeben und

alle Simulationsergebnisse landen in ihren jeweiligen time directories. Die Ordner

werden nach ihrer jeweiligen Zeit benannt, also

namei = 0 + i ∗ writeInterval, mit i = 1, . . . ,endTime

writeInterval. (4.2)

Zur Berechnung kann in derselben Datei eine von uber 80 OpenFoam-Applikationen

ausgewahlt werden. Genugt diese Auswahl nicht, besteht durch den einheitlichen

Aufbau des Quellcodes der verschiedenen Stromungs-Losern den OpenFoam-Nutzern

die Moglichkeit eigene Programme zu schreiben und sie in OpenFoam einzubinden.

In dieser Bachelorarbeit wird die Applikation flameletFoam von Muller, Ferraro und

Pfitzner [10] verwendet um eine Gleichstrom-Wasserstoff-Flamme zu simulieren. Das

zugrunde liegende Prinzip dieser Applikation wurde in Kapitel 3.1 erlautert.

4.2. Die Gleichstrom-H2-Flamme

In Anlehnung an die Experimente der Sandia National Laboratories Livermore,

CA wird im Folgenden der Aufbau und die Simulation einer nicht-vorgemischten

Gleichstrom-Wasserstoff-Verbrennung aufgezeigt.

4.2.1. Geometrie

Der Simulationsaufbau besteht wie in den Experimenten an der Einfluss-Seite aus

dem Querschnitt einer Wasserstoff-Duse mit 3,75 mm Durchmesser mit einer exklu-

siven Wandstarke von 0,55 mm. Die Duse wird von einem runden Stromungsfeld

mit Durchmesser 300 mm umgeben durch welches Luft einstromt. Die Lange dieses

zylindrischen Gebiets und somit die Hohe der simulierten Flamme betragt 800 mm.

An der Ausfluss-Seite wurde eine durchlassige Wand definiert, welche die Stromung

nicht am Austreten hindert.

Der Raum in dem die Simulation ablauft wurde dahingehend vereinfacht, dass in

dem radialsymmetrischen Stromungsfeld der Flamme lediglich ein”Kuchenstuck“

mit Innenwinkel α ≈ 3◦ betrachtet wird (siehe Abbildung 7). Um eine zweidimen-

sionale reaktive Stromung berechnen zu konnen, besitzt dieser Ausschnitt die Dicke

eines einzigen Elements. Aufgrund dieser Vereinfachung kann die Dusenwand an der

Einfluss-Seite nicht rund dargestellt werden, sondern wird linear zwischen den an

22


den Seitenwanden befindlichen Knoten interpoliert. Wie in Abbildung 7 ersichtlich,

wurde die Einfluss-Ebene in drei Teile unterteilt:

1. Zuluft im Gleichstrom (Air)

2. Dusenwand (Tube)

3. Wasserstoff Einlass (Hydrogen)

Abbildung 7: Kuchenstuck des Stromungsfeldes, 20-fach in z-Richtung hochskaliert.

Der Raum, der an die Zuluft-Ebene anschließt ist ein Zusammenschluss (patch) aus

zwei Unterraumen, da so das Gebiet der Grenzschicht zwischen Luft und Flamme in

y-Richtung feiner vernetzt werden kann, ohne unnotigerweise den kompletten Raum

zu verfeinern.

23


4.2.2. Randbedingungen

In OpenFoam werden Randbedingungen in der time directory”0 “(vgl. Kapitel 4.1)

definiert. Fur jede der Variablen alphat, chi, epsilon, k, mut, p, T , U , varZ und

Z wird eine Datei mit den Werten an den jeweiligen Flachen (faces) Inlet, Outlet,

Air, Hydrogen, etc. angelegt. Fur jede Flache muss ein Wert angegeben werden,

welcher auf alle Elemente der Flache einzeln aufgepragt wird. In diesem Fall wird

fur die Geschwindigkeit U , die turbulente kinetische Energie k und die turbulente

Dissipationsrate ε nicht nur ein Wert, sondern eine Liste von Werte angegeben um

den Anfangswerten einen nicht-konstanten funktionalen Verlauf zu geben. Das ist

notwendig um die Reibung einer Rohrstromung, wie sie vor dem Wasserstoff-Einlass

in der Wasserstoff-Duse vorkommt abzubilden. Der resultierende funktionale Verlauf

wird in Abbildung 8 gezeigt.

Abbildung 8: Reibungsbedingte Geschwindigkeits-Verteilung einer Rohrstromung[11].

24


Geschwindigkeit u0

Die Formel fur die reibungsbehaftete Geschwindigkeit an der Wasserstoff-Duse lautet

u(r) =(

1− r

R

) 1nUmax. (4.3)

Nach der 1/7-Regel von Nikuradse (1932) [12] ergibt sich somit fur n = 7 eine echt

turbulente Rohrstromung. Um die Geschwindigkeitsverteilung den Messdaten der

Randbedingungen von Barlow et. al. [3] etwas Naher zu bringen, ist ein Wert von

n = 12 passender. Die Verteilung wurde der Messung mit 20% Helium-Verdunnung

entnommen, da fur den 0%-Fall keine Randmessdaten existieren.

Die Geschwindigkeit ist in der Literatur als Mittelwert Uavg = 296ms

angegeben. Um

die Maximalgeschwindigkeit der Formel (4.3) zu berechnen wurde die Gleichung (4.5)

implementiert und mit einem m ≈ 106 ausgewertet.

Uavg =

∑mi=0

(1− r(i)

R

) 1nUmax

m, mit r(i) = i · R

m(4.4)

⇔ Umax = limm→∞

Uavgm∑mi=0

(1− r(i)

R

) 1n

. (4.5)

Mit dem berechneten Umax und der Formel (4.3) kann nun mit Matlab eine Ge-

schwindigkeitsverteilung auf den Elementen der Hydrogen-Flache im OpenFoam-

Format erstellt werden (siehe Matlab-Skript im Anhang A). Hier mit numElements =

25. Zur Veranschaulichung zeigt Abbildung 9 diesen fur OpenFoam definierten Ge-

schwindigkeitsverlauf auf 25 Zellen.

25


1 type fixedValue;

2 value nonuniform List <vector >

3 25

4 (

5 ( 320.127453 0 0 )

6 ( 319.017667 0 0 )

7 ( 317.863716 0 0 )

8 ( 316.661758 0 0 )

9 ( 315.407417 0 0 )

10 ( 314.095680 0 0 )

11 ( 312.720768 0 0 )

12 ( 311.275962 0 0 )

13 ( 309.753395 0 0 )

14 ( 308.143767 0 0 )

15 ( 306.435973 0 0 )

16 ( 304.616596 0 0 )

17 ( 302.669213 0 0 )

18 ( 300.573415 0 0 )

19 ( 298.303389 0 0 )

20 ( 295.825799 0 0 )

21 ( 293.096516 0 0 )

22 ( 290.055331 0 0 )

23 ( 286.616943 0 0 )

24 ( 282.654543 0 0 )

25 ( 277.967138 0 0 )

26 ( 272.206244 0 0 )

27 ( 264.679769 0 0 )

28 ( 253.649115 0 0 )

29 ( 231.458560 0 0 )

30 );

Abbildung 9: Reibungsbedingte Geschwindigkeits-Verteilung uber die Zellen, ∗ be-zeichnet die Knoten.

26


Turbulente kinetische Energie k0

Die Formel fur die turbulente kinetische Energie lautet

k =1

2

(u′1u

′1 + u′2u

′2 + u′3u

′3

), (4.6)

wobei uii, i = 1,2,3, die Diagonalelemente des Reynold’schen Spannungstensors und

somit die drei Hauptspannungen bezeichnen:

− τij = ρ

u′1u′1 u′1u

′2 u′1u

′3

u′2u′1 u′2u

′2 u′2u

′3

u′3u′1 u′3u

′2 u′3u

′3

. (4.7)

Da fur den Fall ohne Helium-Verdunnung keine Messdaten existieren, werden wie-

der die Daten mit 20%-Verdunnung verwendet und via Matlab ausgewertet (siehe

Matlab-Skript im Anhang B).

Die Diskretisierung der Messdaten ist, wie in Abbildung 10 verdeutlicht, mit 50

Messdaten uber den kompletten Durchmesser zu gering gewahlt, um sie auf eine

Zellenanzahl numElements > 25 aufzupragen. Daher werden die Daten via Inter-

polation auf ≈ 500 Datenpunkte erweitert. Anschließend mussen die Daten zuge-

Abbildung 10: Messdaten der turbulen-ten kinetischen Energie.

Abbildung 11: Verfeinerte Messdaten.

schnitten werden, um ausschließlich Werte innerhalb des Durchmessers der Duse zu

betrachten. Der Durchmesser wurde hier auf d = 2, 04mm gesetzt, da die Messstelle

27


sich offenbar nicht direkt an der Duse, sondern etwas dahinter befand. Somit ver-

bleiben noch ≈ 400 Datenpunkte (vgl. Abbildung 11). Die Grenzen der relevanten

Daten werden in den Abbildungen als rote Linien gekennzeichnet. Die Vereinfachung

der Simulationsumgebung basiert auf gewissen Symmetrie-Eigenschaften, welche ei-

ne symmetrische Halfte der Randbedingungen benotigt. Um Symmetrie zu schaffen,

werden die verfeinerten Messdaten k ∈ Rn symmetrisch zur y-Achse gemittelt,

ki =|ki + kn−i+1|

2, mit i = 1, . . . ,

n

2, (4.8)

und anschließend halbiert, wie es in Abbildung 12 und 13 ersichtlich ist.

Abbildung 12: Gemittelte Messdaten. Abbildung 13: Halbierte Messdaten.

Dieser symmetrische Verlauf der turbulenten kinetischen Energie kann dann wie

im vorherigen Kapitel als funktionaler Verlauf uber die Zellen in der Hydrogen-

Flache im OpenFoam-Format erstellt werden (siehe Matlab-Skript im Anhang B).

In Abbildung 14 ist der Verlauf der turbulenten kinetischen Energie uber die Zellen

aufgezeigt. Die Verteilung deckt sich mit der zuvor eingefuhrten Rohrstromung aus

Abbildung 8. Die Normalspannungen sind reibungsbedingt in Wandnahe wesentlich

großer als in der Rohrmitte. Aus Formel (4.6) folgt, dass ebenso die turbulente kine-

tische Energie in Wandnahe wesentlich großer sein muss, da die Normalspannungen

quadratisch in die Gleichung eingehen. Der Einfluss von Reibung auf die turbulente

Dissipationsrate wird im Folgenden erlautert.

28


1 type fixedValue;

2 value nonuniform List <scalar >

3 25

4 (

5 557.580764

6 543.032708

7 541.837083

8 552.179792

9 551.551042

10 534.833854

11 538.982604

12 524.174444

13 498.732847

14 499.842257

15 475.502743

16 475.871771

17 465.603889

18 445.749861

19 440.988958

20 441.800139

21 458.372153

22 469.240278

23 477.565764

24 510.699062

25 563.164459

26 728.768216

27 986.141213

28 1187.031665

29 1161.620548

30 );

Abbildung 14: Reibungsbedingte Verteilung der turbulenten kinetischen Energieuber 25 Zellen.

29


Turbulente Dissipationsrate ε0

Mit der turbulenten Dissipationsrate bezeichnet man die Geschwindigkeit in welcher

turbulente kinetische Energie in innere Energie eines Systems umgesetzt wird, also

die Energieanderung pro Zeit. Sie geht direkt aus der turbulenten kinetischen Energie

hervor

ε = C34µk

32

`, (4.9)

mit dem Standard-Koeffizienten des k-ε-Modells Cµ = 0, 09.

Das turbulente Langenmaß ` wird in Anlehnung an Zimmermann [12] auf 30% des

hydraulischen Durchmessers gesetzt. Der funktionale Verlauf der turbulenten Dis-

sipationsrate aus Abbildung 15, kann dann mit Formel (4.9) wie zuvor auf den

Mittelpunkt jeder Zelle der Hydrogen-Flache aufgepragt werden:

1 type fixedValue;

2 value nonuniform List <scalar >

3 25

4 (

5 1923052.024530

6 1848282.464998

7 1842181.628383

8 1895178.507018

9 1891942.457187

10 1806582.063616

11 1827643.483394

12 1752843.429348

13 1626789.341142

14 1632220.442816

15 1514463.809202

16 1516227.166959

17 1467419.406178

18 1374567.857973

19 1352604.827066

20 1356338.631936

21 1433364.676469

22 1484643.852893

23 1524330.492975

24 1685687.050456

25 1952010.808505

26 2873515.175124

27 4523119.104202

28 5973416.756979

29 5782635.104072

30 );

30


Abbildung 15: Reibungsbedingte Verteilung der turbulenten Dissipationsrate.

4.2.3. Rechner-Architektur

Die Simulationen wurden auf dem High-Performance-Computing(HPC)-Cluster na-

mens Pacioli des Ulmer Zentrums fur Wissenschaftliches Rechnen (UZWR) durch-

gefuhrt. Der Cluster lauft unter dem Linux Derivat”Scientific Linux“und besitzt

insgesamt 39 Compute-Nodes und einen dedizierten Storage-Node. Aus Kompa-

tibilitatsgrunden des aktuell eingerichteten”parallel environment“mit OpenFoam,

konnte nur via shared memory auf einem der Knoten mit maximal 12 Prozessen

gerechnet werden. Die Knoten basieren auf Intels Sandy Bridge EP (Xeon E5-2630)

und verfugen uber jeweils 128 Gigabyte RAM. Im folgenden Kapitel wird erlautert,

wie auf dieser Rechner-Architektur die Rechenzeiten der Simulation dem Ergebnis

gegenubergestellt wird.

31


4.2.4. Netz-Konvergenzanalyse

Einen Kompromiss zwischen Rechenzeit und Aussagekraft zu finden ist in den meis-

ten Simulationen schwierig. Oftmals interessieren, besonders in der Industrie, nur

komplexe Fragestellungen. Einfache Probleme konnen analytisch gelost werden. Mit

Hilfe einer Netz-Konvergenzanalyse wird die Netzdichte des Stromungsfelds der Ge-

nauigkeit der Losung gegenubergestellt um das optimale Gitter zu finden. Die Re-

chenzeiten und die zugehorigen Knoten-/Elementzahlen werden in Abbildung 16

verglichen.

Netz # Elemente # Knoten RechenzeitenNetz01 7000 14000 1 h 06 minNetz02 28000 57000 7 h 59 minNetz03 60000 120000 31 h 25 minNetz04 113000 230000 110 h 12 min

Abbildung 16: Element-/Knotenanzahlen und zugehorige Rechenzeit.

Am Beispiel des Mischungsbruchs Z an den axialen Stellen x/L = 1/4 und

x/L = 3/8 kann man sehr gut erkennen, dass die Netz-Konvergenzanalyse in diesem

Fall keine verwertbare Losung liefert. Die Mittelwerte der Messdaten werden als ∗dargestellt.

(i) x/L = 1/4 (ii) x/L = 3/8

Abbildung 17: Z: Netz-Konvergenz, – Netz01, – Netz02, – Netz03, – Netz04.

32


Die Verteilung (ii) wurde die Annahme bestatigen, dass mit zunehmender Netz-

dichte die Losung sich zunehmend verbessert. Besonders im Bereich 0-20 mm zeigt

sich durch die zunehmende Krummung eine wesentlich bessere Annaherung an die

Messdaten (Netz04 - blau) im Vergleich zu den Ergebnissen geringer Netzdichten

(Netz01 - cyan, Netz02 - grun). In der axial direkt benachbarten Verteilung (i) tritt

jedoch Gegenteiliges auf.

(iii) x/L = 3/8 (iv) x/L = 1/2

Abbildung 18: O2: Netz-Konvergenz, – Netz01, – Netz02, – Netz03, – Netz04.

Entgegen der Erwartungen verlauft die Losung des feinsten Netzes mit geringe-

rer Krummung im Bereich 0-20 mm als die Losung der einfachen Netze und bewegt

sich zudem fast außerhalb der Standardabweichung der Messungen. In vielen Vertei-

lungen der verschiedenen chemischen Spezies ist teilweise zu beobachten, dass sich,

wie im Verlauf des Sauerstoff-Massenbruchs in Abbildung 18, die Losung nur unwe-

sentlich andert. Aufgrund dieser Beobachtungen ware die Wahl des Netzes nahezu

irrelevant. In Anhang D werden noch weitere Beispiele gezeigt, welche diese Aussage

untermauern.

Unter Berucksichtigung der Rechenzeit wird in den folgenden Kapiteln das Netz03

mit ca. 60000 Elementen und ca. 120000 Knoten zur Auswertung verwendet. Die

Standardabweichungen werden zwar von allen Netzen sehr gut getroffen, jedoch bil-

det Netz03 – in Relation gesehen – die Mittelwerte der Messdaten am besten ab.

33


Vernetzung

Wie in Abbildung 19 zu sehen ist, wird sowohl der gesamte Bereich des Kerns

der Flamme als auch die Vernetzung um die Dusenwand so angepasst, dass die

Zellen an den relevanten Stellen am kleinsten sind. In OpenFoam wird das Netz

Abbildung 19: Netz03, gesamt (o.) und Ausschnitt an der Dusenwand (u.).

durch Definition der Randpunkt-Koordinaten, Verknupfung dieser Randpunkte zu

”blocks“und Unterteilung dieser

”blocks“in allen drei Dimensionen, erstellt. Im fol-

genden Code-Beispiel wird die Definition des verwendeten Netzes und die Begren-

zung des Stromungsfeldes durch”patches“beschrieben. Die Zeile 25 beispielswei-

se bezeichnet das Hexaederformige (hex) Stromungsfeld, welches an die Hydrogen-

Flache anschließt. Es wird durch acht Rand-Koordinaten definiert und ist 300-fach

in x-Richtung, 25-fach in y-Richtung und 1-fach in z-Richtung unterteilt. Die Me-

thode grading(4 0.5 1) vergroßert die Differenz zwischen dem ersten und letzten

Elements jeder Dimension um ihren jeweiligen Faktor und interpoliert die Element-

großen linear um eine gewisse Verfeinerung an den relevanten Stellen zu erhalten.

Bei Werten kleiner Null wird um deren Kehrwert in die andere Richtung verfeinert.

34


1 convertToMeters 1e-3;

2 vertices

3 (

4 (0 0 0)

5 (800 0 0)

6 (800 1.8747 -0.03272)

7 (800 2.4247 -0.04232)

8 (800 9.3735 -0.16359)

9 (800 149.977 -2.61786)

10 (0 149.977 -2.61786)

11 (0 9.3735 -0.16359)

12 (0 2.4247 -0.04232)

13 (0 1.8747 -0.03272)

14 (800 1.8747 0.03272)

15 (800 2.4247 0.04232)

16 (800 9.3735 0.16359)

17 (800 149.977 2.61786)

18 (0 149.977 2.61786)

19 (0 9.3735 0.16359)

20 (0 2.4247 0.04232)

21 (0 1.8747 0.03272)

22 );

23 blocks

24 (

25 hex (0 1 2 9 0 1 10 17) (300 25 1) simpleGrading (4 0.5 1)

26 hex (9 2 3 8 17 10 11 16) (300 15 1) simpleGrading (4 1 1)



29 );

30 edges

31 (

32 );

33 boundary

34 (

35 Air

36 {

37 type patch;

38 faces

39 (

40 (15 7 6 14)

41 (16 8 7 15)

42 );

43 }

44 Hydrogen

45 {

46 type patch;

47 faces

48 (

49 (0 9 17 0)

50 );

51 }

52 Tube

53 {

54 type wall;

55 faces

56 (

57 (17 9 8 16)

58 );

59 }

60 ...

61 );

62 mergePatchPairs

63 (

64 );

35


4.2.5. Zeitliche Konvergenz

Neben der Konvergenz des Netzes ist zusatzlich wichtig, dass jegliche Turbulen-

zen der transienten Stromungssimulation sich in eine stationare Losung einpen-

deln. Schon geringe Bewegungen, wie in Abbildung 20 verdeutlicht, behindern einen

adaquaten Vergleich der Simulationsergebnisse mit den experimentellen Daten. Hier-

zu wurde die Endzeit der Simulation um bis zu 0,2 Sekunden verlangert und daruber

hinaus wurde eine Interpolation gewahlt welche die Werte stromaufwarts mit ein-

bezieht (upwind).

Abbildung 20: Flackern der Flamme in der transienten Berechnung, eingefarbt nachder Temperaturverteilung.

36

Auswertung

5. Auswertung

Die experimentellen Daten der Thermodynamik entstammen den Messungen von

Robert Barlow [13], einem Mitglied der Combustion Forschungseinrichtung der San-

dia National Laboratories Livermore, CA, wobei die Messdaten der Geschwindigkei-

ten an der ETH Zurich ermittelt wurden. Die gemeinsame, hier verwendete Veroffentlichung

2.0, [3] ist auf dem Stand von 2003 (Januar). Es wurden drei verschiedene Messun-

gen durchgefuhrt: 0% , 20% und 40% Helium Verdunnung, wobei die Simulation

der Wasserstoff-Flamme dem unverdunnten Fall entspricht. Pro Messung wurden

an jeder der sieben axialen Messstellen (18, 14, 38, 12, 58, 34, 11) ca. 500 Messungen durch-

gefuhrt und in Tabellen gespeichert. Diese Tabellen sind entsprechend nach der

Helium Verdunnung, der axialen und der radialen Position benannt. Beispielswei-

se he0x1409.sct enthalt Daten mit unverdunntem Wasserstoff (he0) an der axialen

Stelle x/L = 1/4 (x14) im radialen Abstand 9 mm (09). Um den radialen Verlauf an

jeder axialen Stelle auswerten zu konnen, muss jeweils das radiale Mittel bzw. die

radiale Standardabweichung berechnet werden. Die Mittelwerte sind im folgenden

als ∗ und die Standardabweichungen als —— dargestellt.

5.1. Mischungsbruch

In den Messdaten wird der Mischungsbruch nach Bilger definiert als

Fblgr =0, 5(YH − Y2H)/MH − (YO − Y2O)/MO

0, 5(Y1H − Y2H)/MH − (Y1O − Y2O)/MO

, (5.1)

wobei YH, YO die Elementar-Massenbruche der H- bzw. O-Atome bezeichnet. Die

Indizes 1 und 2 ordnen die Massenbruche dem Wasserstoff-Strahl und enstprechend

dem Luft-Strom zu. MH = 1, 008 und MO = 16 sind die molaren Massen der Stoffe.

In Abbildung 21 sind die radialen Verlaufe des Mischungsbruchs an den vier aussa-

gekraftigsten axialen Stellen aufgezeigt.

Die Verteilung (i), unmittelbar stromabwarts der Duse, trifft die Standardabwei-

chungen nicht zufriedenstellend. Der Radius der Flamme ist hier ausschlaggebend,

was direkt auf die Streuung an der Wasserstoff-Duse zuruckzufuhren ist. Diese Pro-

blematik wird in Kapitel 5.4 ausfuhrlicher behandelt. In der Verteilung (ii) sowie

in (iii) tritt gutes Krummungsverhalten auf, es kann jedoch besonders in (iii) noch

37

5.1 Mischungsbruch Auswertung

durch ein feineres Netz lokal verbessert werden (siehe Kapitel 4.2.4). An der axialen

Stelle x/L = 1/1 ist die Flamme ebenfalls noch zu breit, was jedoch durch eine wei-

tere Verlangerung der Rechenzeit nicht verbessern lasst, da die Flamme sich nicht

weiter einpendelt.

(i) x/L = 1/8 (ii) x/L = 3/8

(iii) x/L = 5/8 (iv) x/L = 1/1

Abbildung 21: – Ergebnis Netz03, Radialer Verlauf des Mischungsbruchs an einigenaxialen Stellen.

38

5.2 Massenbruch H2 Auswertung

5.2. Massenbruch H2

Das Verhalten des Massenbruchs ahnelt dem des Mischungsbruches sehr. Der Verlauf

nahe der Duse trifft die Standardabweichung ebenfalls in den Bereichen r ∈ [0, 3] und

r ∈ [13, 18] nicht, jedoch zeigen die Verteilungen (ii) und (iii) besonders ab r > 20

sehr gutes Annaherungsverhaltenan die Mittelwerte. An der axialen Stelle x/L = 1/1

ist die H2-Konzentration aufgrund der stromaufwarts stattgefundenen Reaktion sehr

gering, wodurch im Experiment Messungenauigkeiten auftreten konnen.

(i) x/L = 1/8 (ii) x/L = 3/8

(iii) x/L = 5/8 (iv) x/L = 1/1

Abbildung 22: – Ergebnis Netz03, Radialer Verlauf des Wasserstoff-Massenbruchsan einigen axialen Stellen.

39

5.3 Massenbruch O2 Auswertung

5.3. Massenbruch O2

Der Massenbruch des Sauerstoffs bildet eine gewisse Symmetrie zu den Verteilungen

des Wasserstoffs. Die Stelle x/L = 1/8 wird von der Simulation, wie in den vorheri-

gen Verlaufen, sehr schlecht abgebildet. Der experimentelle Verlauf (ii) wird von der

Simulation ausgezeichnet angenahert, jedoch ist die Streuung im Bereich r ∈ [0, 10]

derart gering, dass die Wahrscheinlichkeit die Mittelwerte sehr gut abzubilden daher

sehr hoch ist. An der Stelle (iii) mangelt es wie in der Wasserstoff-Verteilung eben-

(i) x/L = 1/8 (ii) x/L = 3/8

(iii) x/L = 5/8 (iv) x/L = 1/1

Abbildung 23: – Ergebnis Netz03, Radialer Verlauf des Sauerstoff-Massenbruchs aneinigen axialen Stellen.

falls am Krummungsverhalten in der Halfte r < 26, jedoch wird die zweite Halfte

sehr gut abgebildet. Die Sauerstoff-Verteilung an der Ausstrom-Grenze befindet sich

40

5.4 Variation der Modellkonstanten Auswertung

gut im Bereich der Standardabweichung. Es tritt ein ahnliches, jedoch gespiegeltes

Steigungsproblem auf wie im Verlauf des Mischungsbruchs. Weitere Simulationser-

gebnisse fur den Massenbruch YH2Ound die Temperatur sind in Anhang E anghangt.

5.4. Variation der Modellkonstanten

Die Standard-Parameter des k-ε-Modells aus Kapitel 3.3 werden in OpenFoam au-

tomatisch initialisiert und im log-file der Berechnung wie folgt ausgegeben:

1 kEpsilonCoeffs

2 { Cmu 0.09;

3 C1 1.44;

4 C2 1.92;

5 C3 -0.33;

6 sigmak 1;

7 sigmaEps 1.3;

8 Prt 1; }

Nach Untersuchungen von Zimmermann [12] beeinflusst der Parameter C1”die

Streuungsrate einer zylindrischen Duse, auch bekannt als ’round-jet/plane-jet an-

omaly’“. Das bedeutet, eine Variation des Parameters C1 kann den Durchmesser

der Flamme direkt am Duseneingang – und somit auch an allen axial nachfolgenden

Stellen – schmalern oder erweitern. Abbildung 24 zeigt eine wesentliche Verbesserung

des radialen Verlaufs der verschiedenen Großen an der axialen Stelle x/L = 1/8.

41

5.4 Variation der Modellkonstanten Auswertung

(i) Mischungsbruch Z (ii) Wasserstoff YH2

(iii) Wasser YH2O(iv) Temperatur T

Abbildung 24: Variation der k − ε-Modellkonstanten,– C1 = 1,44 , – C1 = 1,5.

42

Fazit

6. Fazit

Die Ergebnisse zeigen, dass die openFoam Applikation flameletFoam [10] die Losung

der nicht-vorgemischten Wasserstoff-Flamme gut abbildet, da sie sich – bis auf we-

nige Ausnahmen – immer innerhalb der Standardabweichung der Messdaten von

Barlow [3] befindet. Diese wenigen Ausnahmen konnen zudem von der Messunge-

nauigkeit der experimentellen Daten, wie in Abbildung 22 gezeigt, beeinflusst sein.

Die Verwendung verschiedener Randbedingungen hat eine gewisse Anfangswert-

Sensibilitat aufgezeigt und besonders deren funktionale Darstellung uber die ein-

zelnen Zellen hat die Losung erheblich verbessert.

Durch die Anwendung einer Netz-Konvergenzanalyse ergab sich ein optimales Netz

von ca. 60000 Elementen und ca. 120000 Knoten, wobei dieses Ergebnis sehr von der

Rechner-Architektur abhangt auf welcher gerechnet wird. Die Simulation mit Open-

Foam uber ein”parallel environment“mit distributed memory konnte aufgrund des

beschrankten Zeitraums dieser Arbeit nicht eingerichtet werden, wodurch die Anzahl

der Prozessoren auf 12 beschrankt blieb. Fur auf diese Bachelorarbeit aufbauende Si-

mulationen ist es durchaus sinnvoll, mit einem variierten Parameter des k-ε-Modells

zu rechnen, da auch das erhebliche Verbesserungen, besonders im dusennahen Be-

reich, gezeigt hat. Außerdem konnen Tabellen eines reduzierten Wasserstoff-Differen-

tialsystems in die Ordner des OpenFoam-Cases eingebettet werden um die zugrunde

liegende Modellreduktion zu Validieren und letztendlich den ersten Schritt in Rich-

tung Kerosin-Simulation auf Basis eines reduzierten chemischen Systems zu leisten.

43

Literatur Literatur

Literatur

[1] The OpenFOAM Foundation http://www.openfoam.org/ c©2011-2015.

[2] H.Muller, F.Ferraro und M.Pfitzner, Implementation of a Steady Laminar Fla-

melet Model for non-premixed combustion in LES and RANS simulations. 8th

International OpenFoam Workshop. Jeju,Korea, 2013.

[3] R. S. Barlow, Sandia H2/He Flame Data - Release 2.0, http://www.sandia.

gov/TNF/DataArch/H2HeData.html, Sandia National Laboratories (2003).

[4] Barlow, R. S. and Carter, C. D., Combust. Flame 97:261-280 (1994).

[5] Barlow, R. S. and Carter, C. D., Combust. Flame 104:288-299 (1996).

[6] J. Warnatz, U. Maas and R.W. Dibble, Combustion, 4th Edition, Springer-

Verlag (2006).

[7] http://de.wikipedia.org/wiki/Entropie , 06.04.2015, 16:28 Uhr MESZ.

[8] Yvonne Kristen, Thermodynamik,

https://www.uni-ulm.de/fileadmin/website_uni_ulm/nawi.inst.251/

Didactics/thermodynamik/INHALT/HS3.HTM , 07.04.2015, 11:48 Uhr MESZ.

[9] Schwarze, R., CFD-Modellierung: Grundlagen und Anwendungen bei

Stromungsprozessen, Springer Verlag, Berlin, Heidelberg, 2013, ISBN 978-3-

642-24378-3.

[10] H.Muller, F.Ferraro und M.Pfitzner, flameletFoam-Applikation,

https://openfoamwiki.net/index.php/Extend-bazaar/solvers/

combustion/flameletFoam, 24.05.2015, 17:19 Uhr MESZ.

[11] Peter Junglas, Stromungslehre 2 Skript, PHWT Vechta/Diepholz/Ol-

denburg, http://www.peter-junglas.de/fh/vorlesungen/skripte/

stroemungslehre2.pdf 01.06.2015, 19:16 Uhr MESZ.

[12] Ilona Zimmermann, Dissertation: Modeling and Numerical Simulation of Par-

tially Premixed Flames (2009), Institut fur Thermodynamik, Universitat der

Bundeswehr Munchen.

44

http://www.openfoam.org/

http://www.sandia.gov/TNF/DataArch/H2HeData.html

http://www.sandia.gov/TNF/DataArch/H2HeData.html

http://de.wikipedia.org/wiki/Entropie

https://www.uni-ulm.de/fileadmin/website_uni_ulm/nawi.inst.251/Didactics/thermodynamik/INHALT/HS3.HTM

https://www.uni-ulm.de/fileadmin/website_uni_ulm/nawi.inst.251/Didactics/thermodynamik/INHALT/HS3.HTM

https://openfoamwiki.net/index.php/Extend-bazaar/solvers/combustion/flameletFoam

https://openfoamwiki.net/index.php/Extend-bazaar/solvers/combustion/flameletFoam

http://www.peter-junglas.de/fh/vorlesungen/skripte/stroemungslehre2.pdf

http://www.peter-junglas.de/fh/vorlesungen/skripte/stroemungslehre2.pdf

Literatur Literatur

[13] Robert Barlow Short Bio

http://crf.sandia.gov/combustion-research-facility/

working-with-the-crf/crf-staff-2/barlow/, 13.06.2015 13:13 Uhr

MESZ.

[14] Rayleigh and PLIF image, http://www.sandia.gov/TNF/DataArch/

DLRflames.html, 16.06.2015, 00:26 MESZ.

[15] OpenFOAM Foundation: OpenFOAM User Guide. Version 2.3.0, 2014.

[16] k-ε-Modell, http://www.cfd-online.com/Wiki/K-epsilon_models,

16.06.2015, 00:26 MESZ.

45

http://crf.sandia.gov/combustion-research-facility/working-with-the-crf/crf-staff-2/barlow/

http://crf.sandia.gov/combustion-research-facility/working-with-the-crf/crf-staff-2/barlow/

http://www.sandia.gov/TNF/DataArch/DLRflames.html

http://www.sandia.gov/TNF/DataArch/DLRflames.html

http://www.cfd-online.com/Wiki/K-epsilon_models

Abbildungsverzeichnis Abbildungsverzeichnis

Abbildungsverzeichnis

1. DLR-Flamme (m.) mit Darstellung durch Rayleigh-Scattering (l.)

und PLIF (Planar Laser Induced Fluorescence, r.). [14] . . . . . . . . 3

2. Laminare Gegenstrom-Diffusionsflamme, [2]. . . . . . . . . . . . . . . 11

3. Funktionalitat der flameletFoam-Applikation, [2]. . . . . . . . . . . . 13

4. Schematischer Vergleich einer einer direkten Simulation (a) und einer

RANS-Simulation (b), [9]. . . . . . . . . . . . . . . . . . . . . . . . . 15

5. Beispiel einer Scherstromung mit turbulenter Grenzschicht,[6]. . . . . 16

6. Ordner-Struktur eines OpenFoam-Cases, [15]. . . . . . . . . . . . . . 19

7. Kuchenstuck des Stromungsfeldes, 20-fach in z-Richtung hochskaliert. 23

8. Reibungsbedingte Geschwindigkeits-Verteilung einer Rohrstromung

[11]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24

9. Reibungsbedingte Geschwindigkeits-Verteilung uber die Zellen, ∗ be-

zeichnet die Knoten. . . . . . . . . . . . . . . . . . . . . . . . . . . . 26

10. Messdaten der turbulenten kinetischen Energie. . . . . . . . . . . . . 27

11. Verfeinerte Messdaten. . . . . . . . . . . . . . . . . . . . . . . . . . . 27

12. Gemittelte Messdaten. . . . . . . . . . . . . . . . . . . . . . . . . . . 28

13. Halbierte Messdaten. . . . . . . . . . . . . . . . . . . . . . . . . . . . 28

14. Reibungsbedingte Verteilung der turbulenten kinetischen Energie uber

25 Zellen. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29

15. Reibungsbedingte Verteilung der turbulenten Dissipationsrate. . . . . 31

16. Element-/Knotenanzahlen und zugehorige Rechenzeit. . . . . . . . . . 32

17. Z: Netz-Konvergenz, – Netz01, – Netz02, – Netz03, – Netz04. . . . . 32

18. O2: Netz-Konvergenz, – Netz01, – Netz02, – Netz03, – Netz04. . . . . 33

19. Netz03, gesamt (o.) und Ausschnitt an der Dusenwand (u.). . . . . . 34

20. Flackern der Flamme in der transienten Berechnung, eingefarbt nach

der Temperaturverteilung. . . . . . . . . . . . . . . . . . . . . . . . . 36

21. – Ergebnis Netz03, Radialer Verlauf des Mischungsbruchs an einigen

axialen Stellen. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38

22. – Ergebnis Netz03, Radialer Verlauf des Wasserstoff-Massenbruchs

an einigen axialen Stellen. . . . . . . . . . . . . . . . . . . . . . . . . 39

23. – Ergebnis Netz03, Radialer Verlauf des Sauerstoff-Massenbruchs an

einigen axialen Stellen. . . . . . . . . . . . . . . . . . . . . . . . . . . 40

46

Abbildungsverzeichnis Abbildungsverzeichnis

24. Variation der k − ε-Modellkonstanten,– C1 = 1,44 , – C1 = 1,5. . . . 42

25. Netz-Konvergenz, – Netz01, – Netz02, – Netz03, – Netz04. . . . . . 58

26. Netz-Konvergenz, – Netz01, – Netz02, – Netz03, – Netz04. . . . . . 59

27. – Ergebnis Netz03, Radialer Verlauf des Wasser-Massenbruchs an ei-

nigen axialen Stellen. . . . . . . . . . . . . . . . . . . . . . . . . . . . 60

28. – Ergebnis Netz03, Radialer Temperaturverlauf an einigen axialen

Stellen. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61

47

Matlab Skript zur Berechnung der Geschwindigkeitsverteilung

Anhang

A. Matlab Skript zur Berechnung der

Geschwindigkeitsverteilung

1 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

2 % Bachelor Thesis: Matlab Script %

3 % Jan Gabriel %

4 % Institute for numerical mathematics %

5 % Ulm University %

6 % file: velocity.m %

7 % Description: %

8 % nonlinear inlet boundary velocity %

9 % due to fricional inlet tube %

10 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

11 close all

12 clc

13 clear

14 h = figure;

15 hold on

16 grid

17

18 nofelements = 10; % grid dependent value

19 Umean = 296; %m/s

20 powerLaw = 1/12;

21

22 nofnodes = nofelements +1;

23 Umax = calc_Umax(Umean ,powerLaw);

24 y = linspace (0 ,1.875, nofnodes);

25

26 % velocity inlet distribution

27 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

28 for j = 1:nofnodes -1

29 z(j) = y(j) + (y(j+1)-y(j))/2;

30 end

31 f = @(z) Umax * (1 - (abs(z)/1.875)).^( powerLaw); % 1/7 or 1/12 ?

32 x = f(z);

33 a = [zeros(nofnodes -1,1) x’];

34 y_scale = linspace (0,1, nofnodes);

35 for j = 1:nofnodes -1

36 z_scale(j) = y_scale(j) + (y_scale(j+1)-y_scale(j))/2;

37 end

38 b = [z_scale ’ z_scale ’];

39 for i = 1:nofnodes -1

40 plot(a(i,:),b(i,:),’Marker ’,’>’);

41 %plot(a(i,:) ,-b(i,:) ,’Marker ’,’>’);

42 end

48

Matlab Skript zur Berechnung der Geschwindigkeitsverteilung

43 a = zeros(nofnodes ,1);

44 b = y_scale;

45 for i = 1: nofnodes

46 plot(a(i),b(i),’*r’);

47 %plot(a(i),-b(:) ,’*r’);

48 end

49 title([’inlet tube flow , power law: 1/’ num2str (1/ powerLaw) ’, nofelements = ’

num2str(nofelements)])

50 xlabel(’velocity ’)

51 ylabel(’r/R’)

52

53 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

54 % save image

55 cd(’output/boundary ’)

56 print(h,’-dpng’,’pipefriction_u.png’)

57

58 % write data in files

59 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

60 A = [x’ zeros(nofnodes -1,1) zeros(nofnodes -1,1)];

61

62 fileID = fopen(’inlet_velo.txt’,’w’);

63 fprintf(fileID ,’\t type \t \t fixedValue ;\n’);

64 fprintf(fileID ,’\t value \t \t nonuniform List <vector > \n’);

65 fprintf(fileID ,’\t %d \n’, nofnodes -1);

66 fprintf(fileID ,’\t (\n’);

67 formatSpec = ’\t ( %f \t %d \t %d )\n’;

68 fprintf(fileID ,formatSpec , A’);

69 fprintf(fileID ,’\t );’);

70 fclose(fileID);

71 cd(’../..’)

1 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

2 % Bachelor Thesis: Matlab Function %

3 % Jan Gabriel %



6 % file: calc_Umax.m %

7 % Description: %

8 % calculating Umax %

9 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

10 function Umax = calc_Umax(Umean ,powerLaw)

11 nofnodes = 1e6; % grid dependent value

12 y = linspace (0 ,1.875, nofnodes);

13 f = @(y) (1 - (abs(y)/1.875)).^( powerLaw); % 1/7 or 1/12 ?

14 x = f(y);

15 sum = 0;

16 for i = 1: nofnodes

17 sum = sum + x(i);

18 end

19 Umax = Umean*nofnodes / sum; end

49

Matlab Skript zur Berechnung der k-epsilon-Verteilung

B. Matlab Skript zur Berechnung der

k-epsilon-Verteilung

1 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


3 % Jan Gabriel %



6 % file: k_eps.m %

7 % Description: %

8 % nonlinear k-epsilon boundary distri -%

9 % bution due to fricional inlet tube %

10 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

11

12 clc

13 close all

14 clear

15 h = figure;

16 hold on

17 format long

18 grid

19

20 nofelements = 10; % 200% noefelements = 0 !!!

21

22 % Import Data

23 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

24 cd(’~/ Schreibtisch/Bachelor/Matlab/s0020 ’)

25 load s0020.dat

26 cd(’..’)

27 u = s0020 (:,5);

28 varu = s0020 (:,6);

29 varv = s0020 (:,8);

30 radial = s0020 (:,3);

31

32 % Calculate turblulent kinetic energy

33 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

34 k = 0.5.*( varu+varv);

35

36 % refining/interpolating

37 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

38 i = 2;

39 k_fine = linspace(k(1),k(2) ,10);

40 while 1

41 if i == 6 || i == 7

42 k_fine = [k_fine linspace(k(i),k(i+1) ,20)];

43 else

44 k_fine = [k_fine linspace(k(i),k(i+1) ,10)];

45 end

46 if radial(i+1) == 2.5600

50


47 k_fine = [k_fine k(i+1)];

48 break;

49 end

50 i=i+1;

51 end

52 % cut off unnecessary data

53 k_fine = k_fine (51:end -52);

54 radial_fine = -2.04:0.01:2.04;

55 % symmetric interpolation

56 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

57 k_side = zeros ((( length(k_fine) -1)/2)+1,1);

58 for j = 1:(( length(k_fine) -1)/2)

59 k_side(j) = abs(k_fine(j)+k_fine(end -j+1))/2;

60 end

61 k_side(end)= k_fine ((( length(k_fine) -1)/2)+1);

62 radial_side = -2.04:0.01:0;

63 % Scaling

64 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

65 k_scale = k_side (6:end); % 200 datapoints left

66 for i = 1: length(k_scale)/(200/ nofelements)

67 for j=2:(200/ nofelements)

68 k_scale ((i-1) *(200/ nofelements) +1) = k_scale ((i-1) *(200/ nofelements) +1)

+ k_scale ((i-1) *(200/ nofelements)+j);

69 end

70 k_scale(i) = k_scale ((i-1) *(200/ nofelements) +1) /(200/ nofelements);

71 end

72 k_scale = k_scale (1: nofelements); % cut off old data

73 radial_scale = linspace (0,1, nofelements +1);

74 upper = radial_scale (2) /2;

75 lower = 1 - upper;

76 step = upper * 2;

77 radial_scale = -(-lower:step:-upper);

78

79 % Calculating turbulent dissipation rate

80 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

81

82 % Lengthscale (30 % hydraulic diameter)

83 L = 0.03 * 3.75; % [mm]

84

85 epsilon = (0.09^(3/4) .* k_scale .^(3/2))/(L*1e-3); % [m^2/s^2]

86

87 % Plots

88 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

89 cd(’output/boundary ’)

90 plot(radial ,k,’o--b’)

91 title(’turbulent kinetic energy , measured ’)

92 plot ([2.04 2.04] ,[0 2000] ,’r’)

93 plot ([ -2.04 -2.04] ,[0 2000],’r’)

94 xlabel(’radius ’)

95 ylabel(’k’)

51


96 print(h,’-dpng’,’turbulentEnergy_measured.png’)

97

98 h = figure;

99 grid

100 hold on

101 plot(radial_fine ,k_fine ,’o--b’)

102 plot ([2.04 2.04] ,[0 2000] ,’r’)

103 plot ([ -2.04 -2.04] ,[0 2000],’r’)

104 title(’turbulent kinetic energy , fine’)


106 ylabel(’k’)

107 print(h,’-dpng’,’turbulentEnergy_fine.png’)

108

109 h = figure;

110 grid

111 hold on

112 plot(radial_side ,k_side ,’o--b’)

113 title(’turbulent kinetic energy , symmetric ’)


115 ylabel(’k’)

116 plot ([2.04 2.04] ,[0 2000] ,’r’)

117 plot ([ -2.04 -2.04] ,[0 2000],’r’)

118 hold on

119 plot(-radial_side ,k_side ,’o--b’)

120 print(h,’-dpng’,’turbulentEnergy_symmetric.png’)

121

122 h = figure;

123 grid

124 hold on

125 plot(k_side ,-radial_side ,’o--b’)

126 title(’turbulent kinetic energy , halved ’)

127 plot ([0 2000] ,[2.04 2.04] ,’r’)

128 xlabel(’k’)

129 ylabel(’radius ’)

130 print(h,’-dpng’,’turbulentEnergy_halved.png’)

131

132 h = figure;

133 grid

134 hold on

135 a = [zeros(nofelements ,1) k_scale ];

136 b = [radial_scale ’ radial_scale ’];

137 for i = 1: nofelements

138 plot(a(i,:),b(i,:),’Marker ’,’>’);

139 % plot(a(i,:) ,-b(i,:) ,’Marker ’,’>’);

140 end

141 title(’turbulent kinetic energy , scaled ’)

142 meshx = zeros(1, nofelements +1);

143 meshy = linspace(-1,0, nofelements +1);

144 plot(meshx ,-meshy ,’*r’)

145 %plot(meshx , meshy ,’*r’)

52


146 xlabel(’k’)


148 % save image

149 print(h,’-dpng’,’pipefriction_k.png’)

150

151

152 h = figure;

153 grid

154 hold on

155 a = [zeros(nofelements ,1) epsilon ];

156 b = [radial_scale ’ radial_scale ’];

157 for i = 1: nofelements

158 plot(a(i,:),b(i,:),’Marker ’,’>’);

159 %plot(a(i,:) ,-b(i,:) ,’Marker ’,’>’);

160 end

161 title(’turbulent dissipation rate , scaled ’)

162 plot(meshx ,-meshy ,’*r’)

163 %plot(meshx , meshy ,’*r’)

164 xlabel(’epsilon ’)


166 % save image

167 print(h,’-dpng’,’pipefriction_eps.png’)

168

169

170 k_scale = flipud(k_scale);

171 epsilon = flipud(epsilon);

172

173 % Write data in files

174 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

175 fileID = fopen(’inlet_k.txt’,’w’);


177 fprintf(fileID ,’\t value \t \t nonuniform List <scalar > \n’);

178 fprintf(fileID ,’\t %d \n’, nofelements);


180 formatSpec = ’\t %f\n’;

181 fprintf(fileID ,formatSpec , k_scale ’);


183 fclose(fileID);

184

185 fileID = fopen(’inlet_epsilon.txt’,’w’);


187 fprintf(fileID ,’\t value \t \t nonuniform List <scalar > \n’);

188 fprintf(fileID ,’\t %d \n’, nofelements);


190 formatSpec = ’\t %f\n’;

191 fprintf(fileID ,formatSpec , epsilon ’);


193 fclose(fileID);

194 cd(’../..’)

53

Matlab Skript zur Daten-Evaluation

C. Matlab Skript zur Daten-Evaluation

1 % usage:

2 % auto_eval( axial position , Mesh , evaluated Value , color)

3 % positions: 18, 14, 38, 12, 58, 34, 11 (without solidus "/", 1/8 = 18)

4 % Mesh options: Mesh01 , Mesh02 , Mesh03 , Mesh04

5 % Values: Z(mixture Fraction), H2, O2, H2O , T(emperature)

6 % colors: g(reen), b(lue), y(ellow), m(agenta), c(yan), (blac)k

7 % printString = ’on ’: prints image of solution in output/plots

8 % ’off ’: ’figure ’ is needed between auto_eval (...)

9

10 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


12 % Jan Gabriel %



15 % file: auto_eval.m %

16 % Description: %

17 % evaluate experimental data and %

18 % simulation data automatically %

19 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

20 function auto_eval(positionString , meshString , valueString , color , printString)

21 if strcmp(printString , ’on’) == 1

22 fig = figure;

23 end

24

25 hold on

26

27 cd(’~/ Schreibtisch/Bachelor/Matlab/he0sct ’)

28

29

30 heights = zeros (58,1);

31 heights (18) = 18;

32 heights (14) = 30;

33 heights (38) = 48;

34 heights (12) = 56;

35 heights (58) = 52;

36 heights (34) = 55;

37 heights (11) = 50;

38

39 switch valueString

40 case ’Z’

41 low = 1;

42 high = 6;

43 titleString = ’mixture fraction ’;

44 case ’H2’

45 low = 33;

46 high = 40;

47 titleString = ’hydrogen ’;

48 case ’O2’

54


49 low = 15;

50 high = 22;

51 titleString = ’oxygen ’;

52 case ’H2O’

53 low = 42;

54 high = 49;

55 titleString = ’water’;

56 case ’T’

57 low = 9;

58 high = 12;

59 titleString = ’temperature ’;

60 otherwise

61 warning(’Unexpected value. Use Z, H2, O2 , H2O or T instead.’);

62 end

63

64 A = dir;

65 a = size(dir);

66 g = a(1,1);

67

68 % Import file -names for evaluation

69 % ***********************************************************

70 for i =1:g

71 B=A(i).name;

72 size(B’);

73 for k = 1:size(B’);

74 a(1,k) = B(1,k);

75 end

76 end

77 % Create file -name -string

78 % ***********************************************************

79 for h = 1 : (g-2)

80 if h ==1

81 meinfile = strcat(A(h+2).name);

82 else

83 meinfile = [meinfile ,strcat(A(h+2).name)];

84 end

85 end

86 % search for char ’h’ in file -names (he0xXXXX.sct)

87 % ***********************************************************

88 strf = strfind(meinfile ,’h’);

89 size(strf);

90 length = size(meinfile);

91 df = length (1,2);

92 as = 1;

93

94 for i =1:(g-3)

95 for h =strf(i):(strf(i)+11)

96 diffn = h - strf(i)+1;

97 b(i,diffn) = meinfile(h);

98 end

55


99 end

100 sortrows(b);

101

102 % extracting raw -data

103 % ***********************************************************

104 y = 0;

105 e = 0;

106 x = 0;

107 for i = 1:(g-4)

108 str = b(i,:);

109 if(str (5:6)== positionString)

110

111 fid = fopen(str);

112

113 tline = fgetl(fid);





118

119 value = double(str2double(tline(low:high)));

120 while ischar(tline)


122 if (tline ~= -1)

123 if str2double(tline(low:high)) < 1e-7

124 value = [value ; 0];

125 else

126 value = [value ; str2double(tline(low:high))];

127 end

128 end

129 end

130 fclose(fid);

131

132 % calculating mean values and standard deviations

133 % ***********************************************************

134 y = [y mean(value)];

135 e = [e std(value) ];

136 x = [x int16(str2double(str (7:8)))];

137 end

138 end

139 y = y(2:end);

140 e = e(2:end);

141 x = x(2:end);

142

143 errorbar(x,y,e,’r*’) % matlab -errorplot

144

145 cd([’../ OpenFoamData/’ meshString ])

146

147 warning(’off’,’MATLAB:codetools:ModifiedVarnames ’)

148 Table = readtable ([ positionString ’.csv’]);

56


149 warning(’on’,’MATLAB:codetools:ModifiedVarnames ’)

150

151 % variety of values to be evaluated (T, YH2 ,YO2 ,YH2O ,...)

152 % ***********************************************************

153 value = Table {1:end ,{ valueString }};

154

155 z = linspace(0,heights(str2num(positionString)),size(value ,1));

156 plot(z,value ,color)

157 title([’Radial ’ titleString ’ distribution at l/L = ’ positionString (1) ’/’

positionString (2)])

158 xlabel(’radius [mm]’)

159 if ( strcmp(titleString , ’oxygen ’) == 1 || strcmp(titleString , ’hydrogen ’) == 1 ||

strcmp(titleString , ’water ’) == 1 )

160 ylabel ([ titleString ’ mass fraction ’])

161 else

162 ylabel(titleString)

163 end

164 hold off

165 % saving image of plot

166 % ***********************************************************

167 if strcmp(printString , ’on’) == 1

168 cd([’../../ output/plots/’ meshString ])

169 print(fig ,’-dpng’,[valueString meshString ’at’ positionString ’.png’])

170 cd(’../../.. ’)

171 else

172 cd(’../..’)

173 end

174 end

57

Netz-Konvergenz

D. Netz-Konvergenz

In Abbildung 25 und 26 sind an jeweils zwei axialen Stellen die Verteilungen der in

Kapitel 4.2.4 nicht erwahnten Variablen gezeigt. Diese untermauern die These, dass

sich das Netz durch weitere Verfeinerung den Mittelwerten der Losung und dem

Krummungsverhalten deren Verlaufe nicht weiter annahert.

(i) H2 : x/L = 1/4 (ii) H2 : x/L = 5/8

(i) T : x/L = 3/8 (ii) T : x/L = 5/8

Abbildung 25: Netz-Konvergenz, – Netz01, – Netz02, – Netz03, – Netz04.

58

Netz-Konvergenz

(i) H2O : x/L = 1/4 (ii) H2O : x/L = 1/1

Abbildung 26: Netz-Konvergenz, – Netz01, – Netz02, – Netz03, – Netz04.

59

Auswertung: H2O und Temperatur

E. Auswertung: H2O und Temperatur

Massenbruch H2O

(i) x/L = 1/8 (ii) x/L = 3/8

(iii) x/L = 5/8 (iv) x/L = 1/1

Abbildung 27: – Ergebnis Netz03, Radialer Verlauf des Wasser-Massenbruchs aneinigen axialen Stellen.

60


Temperatur

(i) x/L = 1/8 (ii) x/L = 3/8

(iii) x/L = 5/8 (iv) x/L = 1/1

Abbildung 28: – Ergebnis Netz03, Radialer Temperaturverlauf an einigen axialenStellen.

61


Eidesstattliche Erklarung

Ich, Jan Gabriel, Matrikel-Nr. 766929, versichere hiermit, dass ich meine Bachelor-

arbeit mit dem Thema

Modellierung und Simulation einer nicht-vorgemischten

Gleichstrom-Wasserstoff-Verbrennung

selbststandig verfasst und keine anderen als die angegebenen Quellen und Hilfsmittel

benutzt habe, wobei ich alle wortlichen und sinngemaßen Zitate als solche gekenn-

zeichnet habe. Die Arbeit wurde bisher keiner anderen Prufungsbehorde vorgelegt

und auch nicht veroffentlicht.

Mir ist bekannt, dass ich meine Bachelorarbeit zusammen mit dieser Erklarung

fristgemaß nach Vergabe des Themas in zweifacher Ausfertigung und in gebundener

Form beim Studiensekretariat der Universitat Ulm einzureichen habe, zuzuglich des

in § 16c Abs. 9 Rahmenordnung erforderlichen elektronischen Exemplars fur die

Studierendenakte.

Ulm, den

Jan Gabriel

62

Modellierung und Simulation einer nicht-vorgemischten ... · Dar uber hinaus geht mein Dank an...

Documents

Transcript of Modellierung und Simulation einer nicht-vorgemischten ... · Dar uber hinaus geht mein Dank an...