Modellierung der Turbulenz-Chemie- Interaktion in technischen Brennkammern Von der Fakultät für Luft- und Raumfahrttechnik der Universität Stuttgart zur Erlangung der Würde eines Doktors der Ingenieurwissenschaften (Dr.-Ing.) genehmigte Abhandlung Vorgelegt von Peter Theisen aus Wittlich Hauptberichter: Prof. Dr.-Ing. Manfred Aigner Institut für Verbrennungstechnik Deutsches Zentrum für Luft- und Raumfahrt e.V. Mitberichter: Prof. Dr.-Ing. Bernhard Weigand Institut für Thermodynamik der Luft- und Raumfahrt Universität Stuttgart Tag der mündlichen Prüfung: 10. August 2000 Institut für Verbrennungstechnik des Deutschen Zentrums für Luft- und Raumfahrt e.V. 2000 Inhaltsverzeichnis 1 EINLEITUNG................................................................................................................ 1 2 THEORETISCHE GRUNDLAGEN ........................................................................... 3 2.1 MODELLIERUNG DER CHEMISCHEN REAKTION ......................................................... 3 2.1.1 Grundlagen der Beschreibung chemischer Reaktionen und Reaktionssysteme 3 2.1.2 Vereinfachungen der detaillierten Kinetik ...................................................... 6 2.2 MODELLIERUNG DER STRÖMUNG ........................................................................... 12 2.2.1 Beschreibung der Strömungsprozesse ........................................................... 12 2.2.2 Mischungsmodell........................................................................................... 17 2.2.3 Beschreibung turbulenter Strömungsprozesse............................................... 18 2.2.4 Schließungsansätze........................................................................................ 20 2.2.5 Mehrdimensionalität ...................................................................................... 30 2.3 STRÖMUNGSLÖSER ................................................................................................. 32 3 AUFGABENSTELLUNG........................................................................................... 34 4 ENTWICKLUNG EINES CHEMIE-MODULS ...................................................... 37 4.1 ORGANISATIONSSTRUKTUR DER CHEMIE-ROUTINEN IM CHEMIE-MODUL.............. 37 4.2 OPTIMIERTE DATENSPEICHERUNG.......................................................................... 42 4.3 GLEICHGEWICHTSMODUL ....................................................................................... 47 5 NUMERISCHE AUSWERTUNG VON MISCHUNGSBRUCH-MESSUNGEN.54 5.1 VORSTELLUNG DES ANWENDUNGSFALLES ............................................................. 54 5.2 EXPERIMENTELLE UND NUMERISCHE UNTERSUCHUNG DES STRÖMUNGSFELDES ... 60 5.3 AUSWERTUNG DER ZEITLICH GEMITTELTEN MESSUNGEN....................................... 66 6 NUMERISCHE AUSWERTUNG VON SPEZIES-MESSUNGEN ....................... 76 6.1 EXPERIMENTELLE UNTERSUCHUNG........................................................................ 77 6.2 VERGLEICH DER GEMESSENEN PDF’S MIT BERECHNETEN WAHRSCHEINLICHKEITSDICHTEFUNKTIONEN ..................................................................... 89 7 SIMULATION DES TECFLAM-DRALLBRENNERS ........................................ 104 7.1 BERECHNUNG DER MISCHUNG ............................................................................. 104 7.1.1 Problemstellung ........................................................................................... 104 7.1.2 Vorgehensweise bei der Rechnung.............................................................. 106 7.1.3 Ergebnisse der Strömungssimulation .......................................................... 107 7.2 BERECHNUNG DER WÄRMEFREISETZUNG............................................................. 111 7.2.1 Rechnung mit der Annahme chemischen Gleichgewichts .......................... 111 7.2.2 Rechnung mit Berücksichtigung von Nichtgleichgewichtseffekten ........... 113 8 ZUSAMMENFASSUNG UND AUSBLICK ........................................................... 115 LITERATURVERZEICHNIS ......................................................................................... 118 Nomenklatur Lateinische Zeichen: A Fläche, Oberfläche Ar vorexponentieller Faktor der Arrheniusgleichung amr, bmr stöchiometrische Koeffizienten cp spezifische Wärme bei konstantem Druck cµ , cε Modellkonstanten D Diffusionskoeffizient Er Aktivierungstemperatur der Arrheniusgleichung F beliebige Feldgröße f Mischungsbruch g Mischungsbruch-Fluktuation, spezifische freie Enthalpie h spezifische Enthalpie hf 0 Standard-Bildungsenthalpie I innere Energie J Wärmefluss-Vektor K Wärmeleitkoeffizient k turbulente kinetische Energie kfr, kbr Reaktionsrate, Hin- und Rückreaktion l Länge p Druck P Wahrscheinlichkeit Pr Prandl-Zahl Qr Reaktionswärme ℜ universelle Gaskonstante s spezifische Entropie t Zeit T Temperatur u,v,w Komponenten des Geschwindigkeitsvektors u Geschwindigkeitsvektor V Volumen W Molmasse x Ortsvektor Griechische Zeichen: β Temperatur-Exponent in der Arrheniusgleichung γ Isentropen-Exponent ε Dissipationsrate der turbulenten kinetischen Energie µ Viskosität, Mittelwert der Gauß-Verteilung ρ Dichte ρm Partialdichte der Spezies m σ Spannungstensor, Standardabweichung der Gauß-Verteilung ωr chemische Produktion Ψ beliebige Zufallsgröße Φ skalare Feldgröße ξjm Anzahl von Atomen des Elementes j im Molekül m Indizes, tiefgestellt: b die Rückreaktion betreffend f die Hinreaktion betreffend m m-te Spezies im Gasgemisch r r-te Reaktion des Reaktionsmechanismus Indizes, tiefgestellt: c Chemie, chemisch eq chemisches Gleichgewicht 0 Standardzustand G Gauß aG abgeschnittener Gauß Zusammenfassung Die Modellierung von turbulenten Strömungen mit chemischer Reaktion in praxisnahen Brennkammern stellt nach wie vor hohe Ansprüche an heutige Rechnersysteme. Diese Ar- beit liefert einen Beitrag zur sinnvollen Vereinfachung der bestehenden Modelle und zeigt Wege auf, adäquate Ergebnisse in vertretbaren Rechenzeiten zu erreichen. Hierfür werden Ergebnisse zeitlich aufgelöster Messungen an Hand numerischer Methoden analysiert. Die Messergebnisse liefern die Wahrscheinlichkeiten, mit der die verschiedenen Kombinationen der Einflussparameter auf die chemische Reaktion auftreten. Diese Größen, die normaler- weise vom Strömungslöser bereitgestellt werden, sind hier die Grundlage für den Vergleich verschiedener Modellierungsmethoden für die Bestimmung der mittleren chemischen Quell- terme. Insbesondere werden verschiedene Ansätze für die Modellierung der Wahrschein- lichkeitsdichtefunktion der Einflussparameter auf die chemische Reaktion mit den gemesse- nen Schwankungen verglichen. Aus dem Vergleich der Ergebnisse ergeben sich Kriterien, in welchen Fällen vereinfachte Modelle zur Verringerung der Rechenzeit Anwendung fin- den können, ohne Einbußen in der Qualität der Ergebnisse hinnehmen zu müssen. Um die Vielzahl von möglichen Berechnungsmethoden flexibel einsetzen zu können, wurde eine Programm-Struktur entwickelt, die das modulare Hinzufügen neuer Modelle zulässt. Durch die neue Programm-Struktur kann das entsprechende Chemiemodell während der Berechnung in Abhängigkeit der lokal herrschenden thermodynamischen Zustände frei ge- wählt werden. Zur weiteren Effizienzsteigerung wurde ein n-dimensionales Tabellensystem entwickelt, welches sich nach Angabe der geforderten Genauigkeit in Punkto Zugriffszeit und Speicherplatzbedarf automatisch optimiert. Das entstandene Programmsystem zur Be- reitstellung der mittleren chemischen Quellen kommt an einem praxisnahen Methan- gefeuerten Drallbrenner zur Anwendung. Abstract Numerical simulation of turbulent combustion sets high demands even on todays modern computer systems. This work delivers a contribution to simplification of the existing models and shows ways to reach good results in justifiable computation times. Existing instationary measurements in combustion chambers show the distribution of prob- abilities for different combinations of parameters in the chemical reactions. These parame- ters, that can also be computed by the CFD-Code, are the base for the comparison of differ- ent numerical modells for the calculation of the mean chemical sources. The focus lies on methods for modellierung the probability density function of the parameters for chemical reactions. The comparison gives criteria in which cases simplified modells can be used to reduce computation time without sacrificing quality in the results. In order to be able to choose different chemical models flexible, a new program structure was developed, that allows to add new models quickly. Now, the corresponding chemistry model can freely be selected during the calculation in dependence of the locally thermody- namical state. To increase the performance of the calculation, a multidimensional tabulation system was developed, which automatically optimizes itself according to the required access time and the possible memory allocation. The developed program system for the calculation of the mean chemical sources is used to simulate the flame in a methane-fired swirl-burner. 1 1 Einleitung Seit etwa zweihundert Jahren hat sich das Leben auf der Erde drastisch verändert. Durch die Nutzbarmachung von Energie wurde viele Dinge möglich, die vorher undenkbar gewesen wären. Diese industrielle Revolution war jedoch auch notwendig, weil durch den enormen Zuwachs der Weltbevölkerung unter anderem neue Methoden für Produktion und Transport von Nahrungsmitteln benötigt wurden. Energie ist also ein wesentlicher Baustein für den Fortbestand der heutigen Gesellschafts- formen. Der notwendige immense Energiebedarf wird zu einem wesentlichen Anteil aus der Verbrennung von fossilen Brennstoffen gewonnen. Hierbei entstehen jedoch eine Reihe von Schadstoffen, die in den auftretenden großen Mengen starke Schäden an der Umwelt verursachen. Ferner trägt das zwangsläufig entstehende Kohlendioxid zu einer Erwärmung der Erdatmosphäre bei, deren Auswirkungen auf das Weltklima heute noch nicht vorherzu- sagen sind. In der Vergangenheit wurden daher große Anstrengungen unternommen, alter- native Energiequellen zur Verbrennung fossiler Brennstoffe zu finden. Die einzige wirt- schaftlich sinnvolle Methode in die Energiegewinnung aus Kernenergie, die allerdings auf- grund ihrer nicht kalkulierbaren Risiken wieder im Abbau begriffen ist. Andere Energie- quellen wie beispielsweise Sonnenenergie, Windenergie oder Erdwärme können bisher bei weitem nicht die Energiemengen liefern, die die Industrienationen benötigen. Das bedeutet, dass trotz intensiver Forschungsaktivitäten im Bereich der regenerativen Energielieferanten derzeit noch keine Methode in Sicht ist, die die Verbrennung fossiler Brennstoffe auf Dauer ersetzen könnte. Das bedeutet, dass es notwendig ist, die Technik der Verbrennung so weit zu optimieren, dass bei möglichst hohem Wirkungsgrad möglichst wenig Schadstoffe entstehen. Nun könnte man glauben, dass Forschung und Entwicklung in diesem Bereich bereits seit Jahrzehnten abgeschlossen sind und keine wesentlichen Verbesserungen der Prozesse mehr gefunden werden können. Immerhin handelt es sich bei der Technik der Verbrennung um eine der ältesten Techniken, die sich der Mensch zu eigen gemacht hat. Schaut man sich jedoch die Entwicklungen in diesem Bereich genauer an, so muss man feststellen, dass die meisten Fortschritte auf rein empirischen Untersuchungen beruhten. Vielfach wussten die Entwickler nicht, wie sich die einzelnen Teilprozesse der Verbrennung gegenseitig beein- flussen und verbesserten die Prozesse nur anhand ihrer Erfahrung und durch bloßes Auspro- bieren. Die Komplexität des Problems ist bis heute noch nicht vollständig erfasst. Weitere Entwicklungen zur Verbesserung der Technik der Verbrennung setzen jedoch ein tieferge- hendes Verständnis der ablaufenden Prozesse voraus. So ergibt sich hier ein weites Betäti- gungsfeld für die moderne Forschung, dem sich auch diese Arbeit widmet. Woraus resultieren denn die besonderen Schwierigkeiten bei der Entwicklung neuer Tech- niken für die Verbrennung ? Der Prozess der Verbrennung, wie er in heutigen Brennkammern abläuft, setzt sich aus mehreren Teilprozessen zusammen. Diese Teilprozesse zeichnen sich durch sehr komplexe physikalische Vorgänge aus. Die einfachste Art der Verbrennung ist die chemische Umset- zung eines brennbaren Gases (Methan, Propan, Wasserstoff,...). Schon hier liegt eine sehr schwierig zu beschreibende turbulente Gasströmung zu Grunde. Die Strömung vermischt die brennbaren Gase mit Luft. Dieses Gemisch reagiert dann unter Energiefreisetzung zu den Verbrennungsprodukten. Ein komplexes chemisches System beschreibt diesen Reakti- 2 onsvorgang. Auch die chemischen Systeme, die für jeden Brennstoff anders aussehen, sind heute noch Gegenstand intensiver Forschungen. Weitere Teilprozesse sind beispielsweise die Entstehung von Ruß, die Verdunstung von Kraftstoff-Tropfen, Strahlungswärme- Übergänge und der Wärmeübergang an der Wand der Brennkammer. Da all diese Prozesse gleichzeitig ablaufen und sich gegenseitig stark beeinflussen, ist es schwierig, bei Untersuchungen die zu Grunde liegenden Teilprozesse für bestimmte beo- bachtete Phänomene von der Vielzahl von ablaufenden Prozessen zu extrahieren. Diese Aufgabe kann nur anhand aufwendiger Messmethoden bei gleichzeitigem Rechnereinsatz gelöst werden. Neben der sehr komplexen Problematik stellt sich allerdings ein weiteres Problem: Sowohl die numerischen Methoden, als auch die Messmethoden sind unter dem Gesichtspunkt der Genauigkeit und der räumlichen und zeitlichen Auflösung nicht weit genug fortgeschritten, als dass sie alleine belastbare Aussagen über alle Detailfragen herbeiführen könnten. Die numerischen Methoden müssen eine Vielzahl von Rechenoperationen durchführen, so dass eine exakte Lösung auch mit heutigen Hochleistungsrechnern nicht mehr in akzeptabler Zeit gelöst werden kann. Aufgrund der hohen Temperaturen in Brennkammern können meist nur optische Messmethoden eingesetzt werden, die ihrerseits in realistischen Konfigu- rationen durch Streulicht, Ruß, Dichtegradienten, usw. in ihrer Genauigkeit erheblich beein- trächtigt werden. Bei hoher räumlicher und zeitlicher Auflösung entstehen sowohl in der Rechnung als auch in der Messung erhebliche Datenmengen, die kaum noch auswertbar sind. Trotz aller Schwierigkeiten ist es dennoch möglich, einen tatsächlichen Wissenszuwachs durch intensive Forschungsbemühungen in diesem Bereich zu erreichen. Mehr denn je ist hierfür eine enge Zusammenarbeit zwischen Experiment und numerischer Simulation not- wendig. Der Messung kommt hierbei die wichtige Aufgabe zu, Momentaufnahmen von Zu- ständen zu machen, die aus den vielen komplexen Prozessen resultieren. Sie kann zwar Teilprozesse nur bedingt untersuchen, jedoch ist es hiermit möglich, sich einen Überblick über Strukturen und Niveaus der Zustände zu verschaffen, die dann wichtige Erkenntnisse für die Verbrennungsrechnung liefern. Die numerischen Modelle können an Hand der Messergebnisse validiert werden. Sind daraus die Stärken und Schwächen der Rechnersimulation bekannt, können mit numerischen Experimenten räumliche Detailprozesse untersucht, Einzelprozesse extrahiert und Parameterstudien durchgeführt werden, um die wesentlichen Einflussparameter für die optimierte Verbrennung zu finden, die dann für konstruktive Veränderungen an der Brennkammerauslegung verantwortlich sind. In ieser Arbeit wird der Einfluss der Turbulenz der Gasströmung auf die chemische Reak- tion in reinen Gasverbrennungen untersucht. Zeitlich aufgelöste Messergebnisse werden numerisch ausgewertet und mit Programmen zur Berechnung der chemischen Reaktion ge- koppelt. Ziel der Untersuchung ist es, aus der detaillierten Analyse der tatsächlich herr- schenden Zustände in der Brennkammer Rückschlüsse für eine geeignete numerische Be- schreibung der Turbulenz-Chemie-Interaktion zu erhalten. Besonderer Wert soll hierbei auf die Effizienzsteigerung der Modelle gelegt werden. Maßgabe ist, jeweils das einfachste Mo- dell zu finden, welches genügend genaue Ergebnisse liefert. Nachdem Kapitel 2 die für das Verständnis der Arbeit notwendige Theorie erklärt, wird in Kapitel 3 die Aufgabenstellung näher erläutert werden. 3 2 Theoretische Grundlagen In diesem Kapitel werden die theoretischen Grundlagen erläutert, die für das Verständnis dieser Arbeit notwendig sind. 2.1 Modellierung der chemischen Reaktion 2.1.1 Grundlagen der Beschreibung chemischer Reaktionen und Reaktionssysteme Treffen in einem Gas zwei Moleküle aufeinander, so sind sie unter gewissen Voraussetzun- gen in der Lage, eine chemische Bindung einzugehen und somit ein neues Molekül zu bil- den. Bei der Bildung neuer Moleküle kann Wärme an die Umgebung abgegeben (exotherme Reaktion) oder von der Umgebung abgezogen (endotherme Reaktion) werden. Chemische Reaktionen werden durch Reaktionsgleichungen beschrieben. Hierbei wird festgelegt, wie viele Mole der reagierenden Spezies (Reaktanden) zu wie vielen Molen der Reaktionspro- dukte reagieren. Folgende Gleichung zeigt die allgemeine Form einer Reaktionsgleichung [Libby&Williams,1994]: r ra x b xB B B B B B jœ œ (2.1) x repräsentiert ein Mol der beteiligten Spezies α. amr und bmr sind die sogenannten stöchio- metrischen Koeffizienten der Reaktion r, die folgende Bedingung erfüllen: 0r ra b WB B B B  œ (2.2) Wα bedeutet hier die Molmasse der Spezies α. Diese Arbeit beschäftigt sich insbesondere mit der Interaktion zwischen chemischer Reakti- on und Strömung. Für die Strömung ist primär von Bedeutung, wie sich die Spezieszusam- mensetzung während eines Strömungszeitschritts bedingt durch chemische Reaktion ändert. Diese sogenannte Produktionsrate einer Reaktion setzt sich aus der Vorwärts- und Rück- wärtsgeschwindigkeit der Reaktion zusammen und kann wie folgt beschrieben werden: r ra br fr brk c k cB BB B B B X    (2.3) Der linke Term beschreibt die Vorwärtsreaktion und der rechte Term stellt die Rate der Rückwärtsreaktion dar. Der Produktterm der Konzentrationen beschriebt die Wahrschein- lichkeit, mit der Moleküle aufeinander stoßen und eine Reaktion eingehen können. Die Af- 4 finität der Moleküle zur chemischen Reaktion wird durch die Koeffizienten der Hin- und Rückreaktion kfr und kbr beschrieben. Der Arrhenius-Ansatz stellt den Zusammenhang die- ser Koeffizienten zu Stoffeigenschaften und Temperaturniveau dar: fr Tfrn fr fr E A T ek  ; brTbrnbr br E T ek A  (2.4) Efr und Ebr sind die Aktivierungstemperaturen. Die vorexponentiellen Faktoren und die Ak- tivierungstemperaturen können für eine Vielzahl von Elementarreaktionen in der Literatur nachgeschlagen werden [Frenklach Group,1999]. Hier wird schon eine Grundproblematik der Beschreibung chemischer Reaktionen deutlich: Die vom Strömungslöser benötigte Re- aktionsrate ist durch eine stark nichtlineare Funktion beschrieben. Zum einen hängt die Größe der Arrheniuskoeffizienten sehr stark vom gerade herrschenden Temperaturniveau ab und zum anderen gehen in das Massenwirkungsgesetzt (2.3) die stöchiometrischen Koeffi- zienten potenziell ein. Dadurch erhält man auch hier bei stöchiometrischen Koeffizienten größer eins einen nichtlinearen Zusammenhang zwischen der Reaktionsrate und der gerade herrschenden Zusammensetzung. Bei der Verbrennung von technisch interessanten Brennstoffen, wie Kerosin oder Methan wird die Reaktion nicht nur durch eine Reaktionsgleichung beschrieben, sondern durch ei- nige hundert Elementarreaktionen, die in der Summe einen sogenannten Reaktionsmecha- nismus bilden. Während des Reaktionswegs vom unverbrannten Brennstoff-Luft-Gemisch bis hin zu der Verteilung von Spezies, bei der sich in der Summe keine Konzentrationen mehr ändern, dem sogenannten Gleichgewichtszustand, wird eine Vielzahl von Zuständen durchfahren. Zu jedem Zeitpunkt der chemischen Reaktion kann die Brutto-Bildungsrate jeder Spezies bestimmt werden. Der Satz von HESS [Atkins, 1992] besagt, dass Reaktionsgleichungen wie mathematische Gleichungen behandelt werden können. Somit kann durch Summation aller Elementargleichungen die Brutto-Bildungsrate, die der Strömungslöser benötigt, für jede einzelne Spezies wie folgt bestimmt werden: c r r r r W a bB B B BS X œ  (2.5) Bei der Bildung jeder einzelnen Spezies wird entweder Wärme frei oder aus der Umgebung abgezogen. Die Summe dieser Energien wird Wärmefreisetzung der Reaktion genannt und ist eine der Zielgrößen der Berechnung chemischer Reaktionen, da sie der Grund der Tem- peraturerhöhung des Gases während der Verbrennung ist. Die pro Zeiteinheit gebildete E- nergie wird durch folgende Gleichung beschrieben: 0c r r rf r Q a b hB B BB X  %œœ  (2.6) Hierbei ist (∆hf0)α die Standardbildungsenthalpie der Spezies α, die Tabellenwerken [Stull and Prophet,1971] entnommen werden kann. Im chemischen Gleichgewicht einer Reaktion 5 ist die Vorwärtsgeschwindigkeit gleich der Rückwärtsgeschwindigkeit. Dieser Zustand stellt sich bei jeder Reaktion nach einiger Zeit ein, wenn die Randbedingungen für die Reaktion nicht von außen geändert werden. In vielen Strömungsfällen sind die Bedingungen für die chemische Reaktion so, dass sich der Gleichgewichtszustand sehr schnell einstellt. Bei ent- sprechend langsamen Strömungen ist die Änderung der Randbedingungen so langsam, dass die Annahme des chemischen Gleichgewichtszustandes am Ende eines jeden Strömungs- zeitschrittes gerechtfertigt ist. Man stellt fest, dass jede Änderung der Mischung infolge von Speziestransport durch die Strömung quasi sofort verbrannt ist. Der chemische Gleichge- wichtszustand ist für die Strömungssimulation daher interessant, weil er numerisch sehr einfach bestimmbar ist. In diesem Fall gilt folgender einfacher Zusammenhang [At- kins,1992]: / mr mrb a r cm m m W K TS   (2.7) Die Gleichgewichtskonstante Krc der Reaktion r bestimmt sich aus den Arrheniuskoeffizien- ten wie folgt [Libby&Williams,1994]: frr c br k K k  (2.8) Da die Arrheniuskoeffizienten oft nur für die Vorwärtsrichtung bekannt sind, werden in dieser Arbeit die Gleichgewichtskonstanten über die freie Enthalpie, die sogenannte Gibbs- Enthalpie, bestimmt. In Kapitel 4.3 wird die Bestimmung des Gleichgewichtszustandes nä- her erläutert werden. Oftmals ist die Strömungsgeschwindigkeit jedoch so schnell, dass der Gleichgewichtszu- stand in einem Strömungszeitschritt nicht erreicht werden kann. Daher wird im folgenden etwas weiter auf Reaktionsmechanismen eingegangen. Neben komplexen Brennstoffen (z.B. bei Müllverbrennung) sind auch Reaktionsmechanis- men, die die Verbrennung von gängigen Brennstoffen, wie Methan, Propan oder Kerosin beschreiben, immer noch Gegenstand laufender Forschungsaktivitäten. Insbesondere die Verbrennung auf hohem Druckniveau, wie sie in Verbrennungsmotoren und Gasturbinen- Brennkammern vorkommen, konnte bisher noch nicht in ausreichender Form durch kom- plexe Reaktionsmechanismen beschrieben werden. Aus diesem Grund müssen bekannte Mechanismen, wie beispielsweise die GRI-Mechanismen [Frenklach Group,1999], ständig erweitert werden. 6 Bild 2.1.1 zeigt einen noch als einfach zu bezeichnenden Reaktionsmechanismus, und zwar den der Reaktion von Wasserstoff mit Luft. Reaktionsgleichung Hinreaktion Rückreaktion Af mf Ef Ab mb Eb H+H+M ⇔ H2+M 6.4e17 -1 0 5.5e18 -1 51987 O+O+M ⇔ O2+M 6.0e13 0 -503.3 7.2e18 -1 59340 H+O2 ⇔ OH+O 2.2e14 0 8455 1.5e13 0 0 H2+O ⇔ OH+H 1.8e10 1 4479 3.0e13 0 4429 OH+H2 ⇔ H2O+H 2.2e13 0 2592 8.4e13 0 10116 H2O+M ⇔ H+OH+M 5.2e21 -1.5 59386 2.2e22 -2 0 O+H+M ⇔ OH+M 7.1e18 -1.0 0 8.5e18 -1 50830 H2O+O ⇔ OH+OH 5.8e13 0 9059 6.3e12 0 549 H2+O2 ⇔ OH+OH 1.7e13 0 24232 5.7e11 0 14922 N+O2 ⇔ NO+O 6.4e9 1 3170 2.4e11 -0.5 19200 N+NO ⇔ N2+O 1.6e13 0 0 5.0e13 0 37940 N+OH ⇔ NO+H 4.5e13 0 0 1.7e14 0 24500 Bild 2.1.1 Wasserstoff-Luft-Mechanismus [Mühleck,1995] Die Spezies M repräsentiert die möglichen Stoßpartner, die zur Aktivierung der Elementar- reaktion benötigt werden. In der Praxis wird im allgemeinen die Konzentration aller anderen Spezies zur Bestimmung der Reaktionsgeschwindigkeit als Konzentration von M gesetzt. Diese Reaktionen gehören oftmals auch zu den druckabhängigen Reaktionen. Druckabhän- gig sind die Reaktionen, deren Summe der stöchiometrischen Koeffizienten ungleich Null ist. Vergleicht man die Arrheniuskoeffizienten verschiedener Literaturquellen einer Reaktion miteinander, so stellt man fest, dass diese sich oft unterscheiden. Dies liegt daran, dass ver- schiedene Koeffizienten oftmals ähnliche Verläufe liefern. Ferner sind Reaktionsmechanis- men immer nur im Zusammenhang zu sehen, weil diese meist nur als Ganzes validiert wer- den. Es sollte also vermieden werden, einzelne Reaktionen verschiedener Herkunft zu mi- schen. Die Verläufe der Geschwindigkeitskoeffizienten von Reaktionen über der Temperatur zei- gen, ob eine Reaktion endotherm oder exotherm ist. Steigt der Geschwindigkeitskoeffizient bei hohen Temperaturen stark an, so ist die Reaktion endotherm und bedarf einer hohen Temperatur um Produkte zu erzeugen. Der Umkehrschluss gilt in diesem Fall auch. Es exis- tieren allerdings auch Reaktionen, die über weite Bereiche der Temperatur einen mittleren Geschwindigkeitskoeffizienten besitzen. Hier liegen sowohl die Produkte als auch die Eduk- te in weiten Temperaturbereichen in nennenswerten Konzentrationen vor. 2.1.2 Vereinfachungen der detaillierten Kinetik Durch die Leistungssteigerung der Prozessoren in den letzten Jahren können mehrere hun- dert Reaktionsgleichungen zur Bestimmung des Verlaufes einer chemischen Reaktion mit- einbezogen werden. Hierdurch konnten viele Phänomene, wie beispielsweise die Entstehung von Ruß-Vorläufersubstanzen, genauer erforscht werden. 7 Für die Simulation praxisnaher Brennkammern müssen allerdings eine Vielzahl von Be- rechnungen gleichzeitig durchgeführt werden [ vgl. Kapitel 4 ], so dass diese Probleme auch heute noch die Kapazitäten jedes Großrechners erschöpfen. Daher wurde in der Vergangen- heit nach Lösungen gesucht, die Rechenzeit, die für die Bestimmung der chemischen Quel- len benötigt wird, durch sinnvolle Vereinfachungen zu reduzieren. Im Folgenden sind die wesentlichen Strategien zur Vereinfachung des chemischen Systems dargestellt. a) Chemisches Gleichgewicht Das wohl einfachste heute noch verwendete Verbrennungsmodell ist die Annahme unend- lich schneller chemischer Reaktion. Dieses Verbrennungsmodell berücksichtigt weder den zeitlichen Verlauf der Verbrennung, noch die Streckung der Flammenfront. Der chemische Gleichgewichtszustand hängt lediglich von der atomaren Zusammensetzung und einer ther- modynamischen Größe, wie beispielsweise der Temperatur oder der Enthalpie, ab. Das Gleichgewichtsmodell eignet sich sehr gut für Verbrennungsrechnungen, bei denen nur die Strömung von Interesse ist, weil das Temperaturniveau, welches den wesentlichen An- teil der Rückwirkung der Reaktion auf die Strömung darstellt, mit geringem rechnerischen Aufwand hiermit hinreichend genau bestimmt werden kann. Schadstoffe, wie beispielsweise Stickoxide, werden mit diesem Modell nicht sehr gut vorhergesagt. b) Reduzierte Reaktionsmechanismen Der wohl üblichste Weg, die Anzahl der Einflussparameter zu verringern, ist die Reduktion des Reaktionsmechanismus. Es gibt zwei verschiedene Methoden: 1. Reduktion des Mechanismus durch physikalische Betrachtungen: Hierbei werden die Einflüsse der einzelnen Elementargleichungen auf das Gesamtgeschehen analysiert. Für die Strömungsberechnung interessiert meist nur der gleichgewichtsnahe Be- reich. Hier liegen viele der Elementarreaktionen bereits im Gleichgewichtszustand vor, wo- hingegen andere weiterhin den Verlauf der Bruttoreaktion bestimmen. Durch Sensitivitäts- analysen können einzelne Reaktionsschritte eliminiert werden. Bild 2.1.2 zeigt einen solchen reduzierten Mechanismus für ein Methan-Luft-System. Er ergibt sich aus algebraischen Umformungen von 68 Elementargleichungen mit Hilfe von Annahmen, wie beispielsweise dem stationären Zustand der Radikale, und wird von Paczko, Lefdal und Peters (21. Combustion Symposium, 1986) wie folgt beschrieben: 8 CH4 + 2 H + H2O ⇔ CO + 4 H2 CO + H2O ⇔ CO2 + H2 2 H + M ⇔ H2 + M O2 + 3 H2 ⇔ 2 H + 2 H20 Bild 2.1.2 Reduzierter Methan-Luft-Mechanismus [Paczko et al., 1986] 2. Reduktion des Mechanismus durch mathematische Betrachtungen: Um die unter 1. dargestellte Methode durchführen zu können, bedarf es sehr viel Wissen und Erfahrung in Bezug auf den Ablauf von chemischen Reaktionen. Die aufwendige Ana- lyse der einzelnen Reaktionsgleichungen muss für jeden Reaktionsmechanismus neu durch- geführt werden. Durch die Vielzahl der möglichen Brennstoffe und damit der in Frage kommenden Reaktionsmechanismen, ist es wünschenswert, den Vereinfachungsprozess zu automatisieren. Ein sehr eleganter Weg hierfür stellt die ILDM-Methode (Intrinsic Low- Dimensional Manifolds) dar. Diese formale, mathematisch motivierte Methode wurde von Maas und Pope [Maas & Pope,1992] Anfang der 90’ger Jahre aufgegriffen und zählt seither in vielen Instituten weltweit zu den bevorzugten Methoden zur Simulation von Verbren- nungsprozessen mit Nichtgleichgewichtseffekten. Wie bereits erwähnt, handelt es sich bei einem Reaktionsmechanismus um ein nichtlineares Gleichungssystem gewöhnlicher Differentialgleichungen in der Zeit. Durch eine Eigenwert- analyse dieses Systems können Informationen über die Geschwindigkeiten der einzelnen Reaktionen gefunden werden. Auf Grund der Nichtlinearität muss das Gleichungssystem hierzu durch eine Taylor-Reihenentwicklung linearisiert werden. Die Eigenwerte der entste- henden Jakobi-Matrix geben Aufschluss darüber, wie die Elementarreaktionen an Hand ih- rer Geschwindigkeit geordnet werden können. Diese Information entspricht dem Ergebnis aus der Sensitivitätsanalysen der Vereinfachung des Mechanismus von Hand. 9 Bild 2.1.3, Trajektorien der CH4-Verbrennung im Zustandsraum für verschiedene Startpunkte [Schmitt et al.,1995] Trägt man, wie in Bild 2.1.3 dargestellt, den Reaktionsfortschritt des chemischen Systems auf einer zweidimensionalen Mannigfaltigkeit ( ≡ Projektion auf einen Unterraum) auf, so erkennt man, dass sich unabhängig vom Startpunkt die Zusammensetzung im Zustandsraum letztlich nicht nur in einem Gleichgewichtspunkt endet, sondern dass die Reaktionspfade ab einer für jeden Mechanismus charakteristischen Zeit in einer Linie zusammenlaufen. Die Reaktionspfade ergeben sich aus der Eigenwertanalyse [Warnatz et al., 1997], die zu jedem Zeitpunkt über die Eigenvektoren die Richtung der Reaktion beschreibt. Dieses Zusammen- laufen der Linien führt dazu, dass die Anzahl der erforderlichen unabhängigen reaktionsbe- schreibenden Variablen mit zunehmender Annäherung an den Gleichgewichtszustand ab- nimmt. Bei der Berechnung der ILDM-Mechanismen wird der Zeitpunkt bestimmt, an dem die Konzentrationen aller Spezies nur noch vom Mischungsbruch und einigen Spezieskon- zentrationen abhängt. Diese Spezies werden Fortschrittsvariablen genannt, weil sie ein Maß für den Reaktionsfortschritt darstellen. So kann die Lösung des komplexen chemischen Sys- tems durch die Lösung der ILDM ersetzt werden. Die Anzahl der im Strömungslöser zu lösenden Speziesgleichungen ist nicht mehr gleich der Anzahl der im detaillierten Mecha- nismus enthaltenen Spezies, sondern nur noch gleich der Anzahl der Fortschrittsvariablen der ILDM. Die zwei transportierten Spezies werden auch Fortschrittsvariablen genannt, weil sie ein Maß für den Ausbrand der Reaktion sind. In der Praxis werden Haupt-Spezies, wie 10 H2O oder CO2, als Fortschrittsvariablen bevorzugt, weil diese in weiten Bereichen der Brennkammer auch tatsächlich vorkommen. Durch Hinzunahme von weiteren Fortschrittsvariablen ist es möglich, für die Berechnung von noch schnelleren Strömungen mit kürzeren Zeitskalen die charakteristische Zeit weiter nach vorne zu verlegen. Damit erhöhen sich die Rechenzeiten allerdings drastisch. Der Vorteil der ILDM-Methode ist, dass jeglicher Reaktionsmechanismus durch diese Me- thode schnell untersucht und mit wenigen Variablen automatisch parametrisiert werden kann. Die aufwendige physikalische Analyse der Prozesse, die für jeden Brennstoff neu durchgeführt werden muss, entfällt und wird durch einen formalisierten Prozess ersetzt. Ferner ist der Ansatz nahe an den tatsächlich ablaufenden Prozessen. Anders ist dies bei der folgenden Methode, bei der die Lösung eines eindimensionalen Ersatzsystems für die Be- schreibung der Verbrennung in dreidimensionalen Strömungen herangezogen wird. Als ein Nachteil der Methode ist sicherlich zu sehen, dass die numerische Lösung des aus ihr resul- tierenden Eigenwertproblems oft zu Schwierigkeiten führt. Es ist leider in der Praxis nicht der Fall, dass quasi auf Knopfdruck eine Tabelle beliebiger Brennstoffe erstellt werden kann. c) Flamelets: Bei den reduzierten Mechanismen bestimmen die Konzentrationen einiger weniger Spezies den gerade herrschenden Abstand vom chemischen Gleichgewichtszustand. In der Flamelet- Methode [Flamelets ≡ Flämmchen] wird ein anderer Weg beschritten, diesen Abstand zu modellieren: Die Gleichgewichtsannahme besagt, dass die Strömung so langsam ist, dass der Reaktion genügend Zeit zur Verfügung steht, um vollständig auszureagieren. Werden die Geschwin- digkeiten erhöht, dann wird mehr unverbranntes Gas in das Kontrollvolumen transportiert, als die Reaktion umsetzen kann. Je weiter man die Geschwindigkeiten erhöht, desto weiter entfernt sich der Zustand im Kontrollvolumen vom Gleichgewichtszustand. Diese Eigen- schaft macht man sich in der Flamelet-Methode zu Nutze, um Nichtgleichgewichtseffekte zu beschreiben. Hierbei wird die dreidimensionale turbulente Flamme durch ein Ensemble von laminaren, eindimensionalen Flammen ersetzt. Es wird angenommen, dass sich die turbulente Flamme durch das Ensemble von vielen kleinen, dünnen Flämmchen annähern lässt. Diese Flamelets lassen sich in einem Modellexperiment erzeugen und sowohl numerisch als auch experimentell gut beschreiben. Bild 2.1.4 zeigt eine Gegenstrom-Diffusionsflamme, die als Ersatzsystem für die Beschrei- bung von Flamelets dient. Hiermit können Flammenfronten erzeugt werden, bei denen die Geschwindigkeit der Mischung frei einstellbar ist. Das heißt, dass die Geschwindigkeit der Strömung variiert und so verschiedene Nichtgleichgewichtszustände in Abhängigkeit der Mischung und der Geschwindigkeit erzeugt werden können. 11 Bild 2.1.4 Schematische Darstellung einer Gegenstrom-Diffusionsflamme Brennstoff und Oxidator werden senkrecht aufeinander geblasen, so dass eine ebene, lami- nare Diffusionsflamme entsteht. Numerisch kann dieses Problem eindimensional behandelt werden. Die beschreibenden Gleichungen dieser eindimensionalen Flamme können bei Warnatz, Maas & Dibble [Warnatz et al., 1997] nachgeschlagen werden. Nun werden die Strömungsgeschwindigkeiten variiert und die Ergebnisse der Flamelet-Rechnung in Flame- let-Bibliotheken in Abhängigkeit des Mischungsbruchs und der sogenannten skalaren Dissi- pationsrate χ=2D( ∂f / ∂xi )2 abgelegt. Hierbei stellt D den Diffusionskoeffizienten, f den Mischungsbruch (Definition siehe Kapitel 2.2.2) und xi den Ortsvektor dar. In Flammen ist der Geschwindigkeitsgradient a = ∂u/∂r ein charakteristisches inverses Zeitmaß für die Streckung der Flammen. Obwohl in den Transportgleichungen die skalare Dissipationsrate χ auftritt, wird als Streckungsparameter bevorzugt die Scherrate a herange- zogen, da sie das Verlöschen der Flamme besser parametrisiert. Je größer die Streckung der Flamme, desto stärker ist der Effekt des lokalen Verlöschens [Landenfeld, 1999]. Die Flamelet-Methode eignet sich besonders in Bereichen großer Scherkräfte, wie dies bei- spielsweise bei seitlicher Einblasung von Sekundärluft in einen Primärkanal in Fett-Mager- Brennkammern [vgl. Kapitel 5] der Fall ist. Sie wird jedoch in dieser Arbeit nicht benutzt werden, da die in Kapitel 6 behandelten Raman-Messungen die Spezieskonzetrationen lie- fern. Für die detaillierte Auswertung dieser Messungen ist die ILDM-Methode prädestiniert, da sie als Eingabe-Parameter lediglich den Mischungsbruch und die Spezieskonzentrationen benötigt. Die von der Flamelet-Methode erwarteten Geschwindigkeiten wurden jedoch nicht zeitgleich mit den Spezieskonzetrationen gemessen, so dass diese Methode für die Auswer- tung der zur Verfügung stehenden zeitlich aufgelösten Messungen nicht geeignet ist. 12 2.2 Modellierung der Strömung 2.2.1 Beschreibung der Strömungsprozesse Die Gleichungen der Thermodynamik können die Felder der Dichte, der Bewegung und der Temperatur bestimmen. Dafür müssen Anfangs- und Randbedingungen und der Satz von beschreibenden Differentialgleichungen, denen die Felder genügen, bekannt sein. In diesem Kapitel werden die Erhaltungsgleichungen angeführt, welche notwendig sind, um turbulente Strömungsprozesse zu beschreiben. In der Modellvorstellung kann die Strömung eines kontinuierlichen Mediums als Vielzahl von kleinen Volumina konstanter Masse angesehen werden, die sich im Raum bewegen und dabei verschiedenen Kräften ausgesetzt werden. Bild 2.2.1 zeigt ein sogenanntes materielles Volumen V mit den beeinflussenden Größen. Bild 2.2.1 Beispiele von verändernden Größen an einem materiellen Volumen V Der Strömungslöser, der im Rahmen dieser Arbeit zur Anwendung kommt, löst Bilanzglei- chungen für solche sich im Raum bewegenden materiellen Volumina. Diese Bilanzglei- chungen können durch das Reynoldsche Transporttheorem [Schade,1970 ] in Gleichungen für raumfeste Kontrollvolumina (Eulerform) überführt werden. Da die meisten Strömungs- löser sich der Eulerschen Formulierung bedienen, wird im Folgenden nicht weiter auf die Lagrangesche Formulierung eingegangen. Die Transportgleichungen für materielle Volumi- na können in „compressible-fluid dynamics“ [Thompson, 1972] nachgelesen werden. Bei der Herleitung der Eulerschen Transportgleichungen wird von einem raumfesten Volu- men ausgegangen, welches extensive Größen als Eigenschaften besitzt. Diese extensiven Größen F(t) lassen sich durch die zugehörige Dichte f(r,t)=dF/dV durch folgende Integration bestimmen: 13 ,F t f r t dV 8  ¨ (2.9) r ist der Ortsvektor, t die Zeit und dV ein differenzielles Volumenelement im Volumen Ω. Eine zeitliche Veränderung der extensiven Größen F des Volumens Ω kann drei Ursachen haben: Änderung durch eine Stromdichte Φf . dA durch die Oberfläche ∂Ω (Diffusion, Wärmelei- tung, Reibungskräfte, Konvektion, ...) Änderung durch einen Quellterm qf im Inneren des Volumenelementes, wobei qf die pro Zeit- und Volumeneinheit gebildete Menge an F beschreibt (z.B. durch chemische Reakti- on). Änderung durch Fernwirkung sf von außerhalb in das Innere des Volumenelementes (z.B. durch Gravitation, Strahlungsabsorption, ...) Die gesamte zeitliche Änderung der jeweils betrachteten Größe F F fdV t t8 s ss s¨ (2.10) lässt sich durch Integration des Flusses über die gesamte Oberfläche ∂Ω und Integration der Quellterme über das gesamte Volumenelement Ω berechnen mit f f A f dV dA q dV s dV t8 8 8 s '  s¨ ¨ ¨ ¨ (2.11) Mit Hilfe des Gaußschen Integralsatzes, der eine Beziehung zwischen einem Oberflächenin- tegral und einem Volumenintegral herstellt, ergibt sich folgender Zusammenhang: f f f f dV dV q dV s dV t8 8 8 8 s ‹ ¸'  s¨ ¨ ¨ ¨ (2.12) Dies ist nun die allgemeine Form einer Transportgleichung für extensive Zustandsgrößen eines Volumenelementes dV. In der Literatur wird anstatt der integralen Form meist die differenzielle Form angegeben. Die differenzielle Form der Bilanzgleichung lautet: f f f f q s t s ‹ ¸'  s (2.13) 14 Aus dieser Gleichung lassen sich die Bilanzgleichungen für Impuls, Masse und Energie wie folgt ableiten: a) Massenerhaltung Für den Transport der Gesamtmasse m ergibt sich folgender Zusammenhang: fm = ρ Φm = ρu qm = 0 sm = 0 0u t S Ss ‹ ¸ s (2.14) b) Spezies-Massenerhaltung Für den Transport einer einzelnen Spezies m,i ergibt sich folgender Zusammenhang: fm,i = ρ i Φm,i = ρ i u i + ji qm,i = Mi ωi sm,i = 0 i i i i iu j Mt S S Xs ‹ ¸ ‹ ¸ s (2.15) ji stellt die Diffusionsstromdichte der Spezies dar und wird später näher erläutert. c) Impulserhaltung Für den Transport des Impulses mu ergibt sich folgender Zusammenhang: fmu = ρu Φmu = ρ u u + p qmu = 0 smu = ρg u uu p g t S S Ss ‹ ¸ ‹ ¸ s (2.16) 15 p stellt hierbei den Drucktensor dar, der später näher erläutert wird. g ist die Erdbeschleu- nigung. d) Energieerhaltung Für den Transport der Energie ergibt sich folgender Zusammenhang: fI = ρI ΦI = ρ u I + p u + jq qI = ρ u g sI = qr : rqI uI j p u q ugtS S Ss ‹ ¸ ‹  s (2.17) jq stellt hierbei die Wärmestromdichte dar. Die oben gezeigten Transportgleichungen sind erst dann geschlossen, wenn für den Wärme- fluss jq, den Diffusionsfluss ji und den Drucktensor p geeignete Ansätze gefunden werden. Hier werden folgende empirische Gesetze herangezogen: a) Das Newtonsche Spannungsgesetz An Hand experimenteller Untersuchungen wurde folgender Ansatz für den Drucktensor p gefunden: p pE T (2.18) Dabei ist E der Einheitstensor, der mit p multipliziert den hydrostatischen Anteil des Druck- tensors beschreibt. Für den viskosen Anteil ergibt sich unter Vernachlässigung der Volu- menviskosität folgender Zusammenhang: 2 3 Tu u u ET N   ¯  ‹ ‹  ‹ ¸¡ °¡ °¢ ± (2.19) b) Das Fouriersche Wärmeleitungsgesetz In Verbrennungsprozessen wird die Wärmestromdichte durch folgenden Ansatz beschrie- ben: q i i i j T h jM ‹ œ (2.20) 16 Der erste Term stellt die Wärmeleitung dar. Der zweite Anteil beschreibt den Enthalpiefluss durch Speziesdiffusion. λ ist der Wärmeleitkoeffizienten, hi ist die spezifische Enthalpie der Spezies i und ji ist der Diffusionsfluss der Spezies i. c) Das Ficksche Diffusionsgesetz Der wesentliche Anteil der Diffusion ist in Verbrennungsprozessen die Druckdiffusion. Der Diffusionsfluss der Spezies i wird hierbei analog zum Fickschen Diffusionsgesetz beschrie- ben: M ii i i i w j D x x S ‹ (2.21) DiM ist hierbei der mittlere Diffusionskoeffizient der Spezies i in Bezug zum restlichen Gasgemisch. xi stellt den Molenbruch der Spezies i dar. 17 2.2.2 Mischungsmodell Durch die im Kapitel 2.2.1 dargestellte Vorgehensweise kann die Verteilung der einzelnen Spezies-Konzentrationen an jedem Ort in der Brennkammer für jeden Zeitpunkt berechnet werden. Bei den gängigen Verbrennungsmodellen genügt meistens die Kenntnis einiger weniger Spezies und das Verhältnis, in dem sich eintretende Stoffströme in der aktuellen Zelle vermischt haben, um die chemischen Quellen zu berechnen. Um dieses Verhältnis zu beschreiben, wird der Mischungsbruch f definiert. Dieser leitet sich wie folgt ab [Lib- by&Williams, 1994]: Die Massenkonzentration Zj eines Atoms j in einem Kontrollvolumen ist definiert durch: / m mj j jm m Z W W wYœ (2.22) Wj ist das Molgewicht des Atoms j, Wm stellt das Molgewicht der Spezies m dar und wm ist der Massenanteil der Spezies im Volumen und ξjm ist die Anzahl der Atome j im Molekül m. Unter der Annahme, dass alle Diffusionskoeffizienten gleich groß sind, ist es unerheblich, ob die Elemente rechnerisch in atomarer Form oder in ihrer molekularen Bindung transpor- tiert werden. Multipliziert man nun die Speziesgleichungen mit ρ(Wj/Wm)ζjm und summiert über alle Spezies auf, so ergibt sich für jedes Element j eine gleichlautende Erhaltungsglei- chung. Sie lautet: j j jZ Z D Zt S S Ss   ¯ ‹  ‹ ‹¡ °¢ ±s u (2.23) Der chemische Quellterm, wie er in den Speziesgleichungen vorkommt, tritt in den Glei- chungen des atomaren Transports nicht auf, da auf Grund der Massenerhaltung die Anzahl der Elemente bei chemischen Reaktionen gleich bleibt. Wenn als weitere Voraussetzung zwei unterschiedlich zusammengesetzte, in den Brennraum eintretende Stoffströme ange- nommen werden, kann mit der Normierung von Zj mit den atomaren Konzentrationen am Eintritt Zj1 und Zj2 der Mischungsbruch eingeführt werden: 2 1 2 j j j j Z Z f Z Z   (2.24) Der Mischungsbruch ist eine massenbezogene Größe. Die im Kontrollvolumen vorliegen- den Atome sind das Resultat der Mischung der beiden Eintrittsströme. Der Mischungsbruch sagt nun aus, in welchem Massenverhältnis die beiden Stoffströme gemischt werden müs- sen, um genau die in der Zelle vorliegenden Atomkonzentrationen zu erhalten. Das bedeutet hier auch umgekehrt, dass durch die Kenntnis der Zusammensetzung der eintretenden Stoff- ströme und des Mischungsbruchs die atomare Zusammensetzung der Zelle bekannt ist. 18 Die Transportgleichung für den Mischungsbruch leitet sich aus Gleichung (2.15) wie folgt ab: < >f f D f t S S Ss ‹  ‹ ‹s u (2.25) 2.2.3 Beschreibung turbulenter Strömungsprozesse Betrachtet man die Geschwindigkeiten über einen längeren Zeitraum an einem festen Punkt in einer turbulent durchströmten Brennkammer, so erscheint ein zeitlich stochastischer Ver- lauf. Der Grund hierfür sind turbulente Schwankungsbewegungen, die die Werte der Ge- schwindigkeiten von ihrem zeitlichen Mittel abweichen lassen. Aus diesen Schwankungen der Geschwindigkeit resultieren auch Schwankungen aller anderen Größen, wie Temperatur, Dichte oder Spezieskonzentrationen. Bild 2.2.2 Schwankungen einer fluktuierenden Größe In Bild 2.2.2 ist der zeitliche Verlauf einer fluktuierenden Größe dargestellt. Diesen zeitli- chen Verlauf an Hand der Transportgleichungen, die in Kapitel 2.2.1 beschrieben wurden, numerisch zu berechnen, ist auch mit den heutigen Rechnern nur für sehr einfache Anwen- dungsfälle möglich. Daher wird zunächst folgende Aufspaltung der schwankenden Größen φ(t) vorgenommen [Reynolds, 1885]: (2.26) φ t( ) =φ t( )+ ′ φ t( ) 19 Der Momentanwert φ(t) bestimmt sich aus seinem zeitlichen Mittelwert tG und dem Schwankungswert tGa . Per Definition ist der zeitliche Mittelwert der Schwankungen gleich Null ( 0Ga  ). Diese Vorgehensweise wird Reynolds-Mittelung genannt. Für den Mittelwert des Produktes zweier schwankender Größen gilt: 1 2 1 2 1 2GG G G GGa a ¸ (2.27) Der zweite Term der rechten Seite wird als Korrelation der turbulenten Schwankungsgrößen bezeichnet. Diese Korrelation wird in dieser Arbeit eine zentrale Rolle spielen, da bei der Modellierung dieses Terms untersucht werden muss, ob die Schwankungen beider Größen in einem funktionalen Zusammenhang stehen oder ob sie vollkommen unabhängig vonein- ander sind [ vgl. Kapitel 2.2.4 ]. Bei Verbrennungsprozessen treten typischerweise große Dichteschwankungen auf. Es er- weist sich hier als zweckmäßig, die sogenannte Favre-Mittelung [Favre, 1969] einzuführen, in der bei der Mittelung der fluktuierenden Größe die Dichteschwankungen mit einbezogen werden. Eine Favre-gemittelte fluktuierende Größe G berechnet sich wie folgt: SGG S  (2.28) Auch hier ergibt sich wieder die Aufteilung in Mittelwert und Schwankungswert: t t tG G Gaa  (2.29) Im Gegensatz zu der Reynolds-Mittelung ist bei der Favre-Mittelung 0Gaa v aber 0SGaa  . In der folgenden Gleichung ist das Produkt zweier fluktuierender Größen einmal Favre-gemittelt und einmal Reynolds-gemittelt worden. Die Reynolds-Mittelung liefert im Fall turbulenter Strömung mit Verbrennung mehr Terme als die Favre-Mittelung, so dass hierbei mehr Terme modelliert werden müssen. 1 2 1 2 1 2 1 2 1 2 1 2 2 1SG G SG G SGG SG G SG G G S G G S Ga a a a a a a a    (2.30) Daher wird bei der Simulation turbulenter Brennkammer-Strömungen meist die Favre- Mittelung benutzt. Der Nachteil der Favre-Mittelung ist, dass die modellierten Terme mess- technisch nur schwer zu bestimmen sind, da die Dichteschwankungen gleichzeitig mit er- fasst werden müssen. Die in Kapitel 2.2.1 hergeleiteten Transportgleichungen beschreiben den Verlauf von Mo- mentanwerten. Für turbulente Strömungen können diese Gleichungen durch einsetzen von 20 Gleichung (2.29) in die folgende gemittelte Form umgewandelt werden [Lib- by&Williams,1994]: Massenerhaltung: 0u t S Ss ‹ ¸ s  (2.31) Spezieserhaltung: i i i i i iw uw D w u w MtS S S S Xs aa aa ‹ ¸ ‹ ¸  ‹ s   (2.32) Impulserhaltung: u uu p u u gtS S S Ss aa aa ‹ ¸ ‹ ¸ s   (2.33) Energieerhaltung: rTh p uh h qt t S S SOs s ‹ ¸ ‹ ¸ ‹ s s    (2.34) Eine Transportgleichung für die Varianz einer skalaren Strömungsgröße lässt sich aus den Gleichungen für Konvektion und Diffusion und Modellannahmen für die Produktion und Dissipation dieser Größe ableiten. Für die Varianz g einer konservativen skalaren Größe wie beispielsweise dem Mischungsbruch ergibt sich [Mühleck,1995]: 21 2u Pr g gg g g g C f C g t k S N FS N S   ¯ ¬s ­ž¡ °­ž ‹ ¸  ‹ ‹ ‹ ­¡ °ž ­­žs Ÿ ®¡ °¢ ± (2.35) 2.2.4 Schließungsansätze Mittels der in Kapitel 2.2.3 hergeleiteten Gleichungen können turbulente Strömungen simu- liert werden. Jedoch müssen die Gleichungen zuerst noch vollständig geschlossen werden. 21 Es ist notwendig, sowohl die mittleren Reaktionsraten i iM X , als auch die Korrelationen iu wS aa aa und i ju uS aa aa zu modellieren. Das Geschwindigkeitsmoment zweiter Ordnung i ju uaa aa wird als Reynoldsscher Span- nungstensor bezeichnet. Seine Modellierung ist Gegenstand der sogenannten Turbulenzmo- delle. Das wohl bekannteste Turbulenzmodell, welches auch in dieser Arbeit zum Einsatz kommt, ist das k-ε-Modell nach Launder und Spalding [Launder&Spalding,1972]. Hierbei wird die turbulente kinetische Energie k definiert als: 21 2 iuk S S aa œ (2.36) Die Dissipationsrate ε der turbulenten kinetischen Energie k wird wie folgt definiert [Laun- der&Spalding,1972]: :TT u uF O aa aa ‹ ‹ (2.37) Neben diesen beiden Größen wird die Wirbelviskositätshypothese von Boussinesq [Boussi- nesq,1877] benutzt, die einen linearen Zusammenhang zwischen den durch Turbulenz ent- stehenden Schubspannungen und dem Gradienten der mittleren Geschwindigkeiten an- nimmt. Der Faktor νT ist der turbulente Anteil der Viskosität und wird nach folgender Glei- chung bestimmt: 2 T kcOO S F   (2.38) Die Wirbelviskosität wird in vielen gebräuchlichen Turbulenzmodellen als isotrop ange- nommen, was nach Jones [Jones, 1980] eine geeignete Approximation für viele praktische Anwendungen ist. Die Größen k und ε werden an Hand folgender Transportgleichungen bestimmt: :Tk uk k u u ut S S SO S SFs ‹ ¸ ‹ ¸ ‹  ‹ aa aas      (2.39) 21 2:Tdiv u div c u u u ct k kSF F FS F SO F S Ss  ‹   ‹ aa aas        (2.40) Die Konstanten cν , c1 und c2 werden an Hand experimenteller und rechnerischer Betrach- tungen empirisch bestimmt. 22 Das k-ε-Modell ist das bekannteste und wohl auch das meistverwandte Turbulenzmodell. Dies liegt daran, dass es für viele Anwendungsfälle gute Ergebnisse liefert und die Rechen- zeiten nicht sehr hoch liegen. Bei stark verdrallten Strömungen reicht das k-ε-Modell oft nicht mehr aus. Für diese Art von Strömungen kommt daher in Rahmen dieser Arbeit ein modifiziertes k-ε-Modell (MKE) zum Einsatz [Hirsch, 1995]. Es verwendet gegenüber dem Standard-k-ε-Modell (STKE) die selben Transportgleichungen für die turbulente kinetische Energie k und ihre Dissipationsrate ε. Es führt jedoch eine Drallkorrektur des im STKE verwendeten Boussinesq-Schließungsansatzes ein, welcher eine Proportionalität zwischen dem Reynolds-Spannungs-Tensor und dem Schergeschwindigkeitstensor des mittleren Ge- schwindigkeitsfeldes annimmt. In Drallströmungen erweist sich diese Annahme als eine Schwachstelle, wobei insbesondere die Zentrifugalwirkung auf die turbulente Strömung dadurch schlecht beschrieben wird. Im MKE-Modell wird dieser Einfluss berücksichtigt, indem unter gewissen Annahmen eine Drallkorrektur für das k-ε-Modell aus dem Algebrai- schen Spannungsmodell abgeleitet wird. In der Vergangenheit konnte gezeigt werden, dass die Berechnung stark verdrallter Strö- mungen mit diesem Modell wesentlich besser modelliert werden konnten [Ohl, 1997]. Bei den zu schließenden Termen fehlt an dieser Stelle noch der gemittelte chemische Quell- term i iM X . Der Momentanwert hängt dabei stark nichtlinear von seinen Einflussparametern ab. Einflussparameter sind hierbei alle Spezies und eine thermodynamische Größe, wie bei- spielsweise Temperatur oder Enthalpie. Durch die Nichtlinearität des Momentanwertes der chemischen Quelle folgt: , , , ...r r mT YX X Sv  (2.41) Um den zeitlichen Mittelwert der chemischen Quellen adäquat zu bestimmen, müssten also die Quellterme für alle durch turbulente Schwankungen vorkommenden Kombinationen von Momentanwerten der Einflussgrößen gemittelt werden. Das bedeutet auch, dass bekannt sein muss, wie oft eine bestimmte Kombination von Einflussgrößen vorkommt. Als Maß für diese Anzahl von Ereignissen wird im Folgenden die sogenannte Wahrscheinlichkeitsdich- tefunktion eingeführt: Die Verteilung einer fluktuierenden Größe wird in einer Strömung durch die Wahrschein- lichkeitsdichtefunktion PDF (Propability Density Function) beschrieben. Ihr Funktionswert gibt die Wahrscheinlichkeit an, dass die fluktuierende Größe Φ in ein Intervall [a,b] fällt [Libby&Williams,1994]: < > Prob b a a b P dG G Gb b  ¨ (2.42) Da die Wahrscheinlichkeit, dass die fluktuierende Größe Φ im Intervall [-∞,∞] liegt gleich eins ist, muss die Normierungsbedingung gelten: 1 P dG G d d  ¨ (2.43) 23 Der Mittelwert der Zufallsgröße Φ ergibt sich nun aus der Integration der wahrscheinlich- keitsgewichteten Zufallswerte über den gesamten Wertebereich der Zufallsgröße: P dG G G G d d  ¨ (2.44) Die Varianz ist definiert durch: 2' 2 P dG G G G Gd d  ¨ (2.45) Ist die Wahrscheinlichkeitsdichtefunktionen der Einflussgrößen bekannt, können die aus den Momentanwerten resultierenden chemischen Quellen gemittelt werden. Jede einzelne der berechneten chemischen Quellen wird mit der zugehörigen aus der PDF entnommenen Wahrscheinlichkeit gewichtet und über alle möglichen Zustände wie folgt aufintegriert: 1 2 1 2 1 2 1 2, , , , , n n n ni i x x x x x x P x x x dx dx dxG G ¨ ¨ ¨! ! ! ! (2.46) Wie die Momentanwerte φi(xi) der chemischen Quellen bestimmt werden, wurde in Kapitel 2.1 beschrieben. Im Folgenden soll auf die Ermittlung der Wahrscheinlichkeitsdichtefunkti- on eingegangen werden. Nach Gleichung (2.46) ist es notwendig, die Wahrscheinlichkeit jeder Kombination von Einflussgrößen der Chemie zu kennen. Es ist möglich, eine Transportgleichung für diese mehrdimensionale PDF aufzustellen und an Hand dieser Gleichung die PDF an jeder Stelle im Raum zu ermitteln [Pope, 1976]. Eine genauere Beschreibung der Transportgleichung für die Wahrscheinlichkeitsdichte kann bei Mühleck (1995) nachgelesen werden. Aufgrund der hohen Dimension der transportierten PDF, müssen extrem viele Rechenopera- tionen zur Lösung durchgeführt werden. Auch wenn diese Probleme meist über einen Mon- te-Carlo-Ansatz optimiert werden, ist es bisher noch nicht gelungen, dieses Problem für praktische Anwendungsfälle in angemessener Rechenzeit zu lösen. Daher wird die berechnete PDF durch eine angenommene, analytisch bestimmbare Wahr- scheinlichkeitsdichtefunktion ersetzt. Dies hat den Vorteil, dass pro Bestimmungsparameter der PDF nur eine Transportgleichung gelöst werden muss. In der Praxis werden mehrere verschiedene analytisch bestimmbare PDF’s verwandt. Wel- che dieser verschiedenen Formen eingesetzt wird, hängt vom physikalischen Problem ab. Die Auswirkung der Form der PDF auf die Lösung des Chemiemoduls wird Gegenstand von Kapitel 5 und 6 sein. Obwohl in der Vergangenheit eine Reihe von verschiedenen PDF-Formen angewandt wur- den (Borghi 1988), haben sich folgende drei Formen weitgehend durchgesetzt: 24 a) Doppel/Dreifach-Delta-Funktion Delta-Funktionen haben den großen Vorteil, dass sie nur sehr wenig Rechenzeit in An- spruch nehmen. Auch die Integration über die Delta-PDF benötigt nur wenig Rechenzeit, da pro Dimension der PDF nur für jeden Delta-Peak der chemische Zustand bestimmt werden muss. Die Gleichung der Doppel-Delta-PDF lautet [Serag-el-Din,1977]: 1 2P x a x x b x xE E   (2.47) P(x) hängt von insgesamt vier Parametern ab, und zwar von der Höhe der Delta-Peaks und deren Position. Daher müssen vier Bestimmungsgleichungen gefunden werden, um die PDF eindeutig zu bestimmen. Aus der Normierungsbedingung (2.43) und den Bedingungen für die ersten beiden Momente, Mittelwert (2.44) und Varianz (2.45) und der Annahme gleicher Höhen der Peaks ergibt sich folgender Zusammenhang für P(x): 2 2P x a x x x b x x xE E  ¯   ¯a a    ¡ ° ¡ °¡ ° ¡ °¢ ± ¢ ± (2.48) mit a = ½ , b = ½. Unterschreitet der Wert der Standardabweichung ein Tausendstel des Betrages des Mittel- wertes, so wird lediglich ein Peak mit der Höhe eins an die Stelle des Mittelwertes gesetzt. Die Einflussparameter der Chemie liegen in einem endlichen Gültigkeitsbereich. So bewegt sich der Mischungsbruch beispielsweise zwischen 0 und 1 und die Fortschrittsvariablen zwischen 0 und ihrem Gleichgewichtswert, der wiederum vom momentanen Mischungs- bruch abhängig ist. Es kommt häufig vor, dass einer der beiden äußeren Peaks außerhalb dieses Gültigkeitsbereichs liegt. In diesem Fall wird dieser Peak genau auf den Rand des Bereichs gelegt. Nach Khalil [Khalil,1976] berechnen sich die Parameter a und b in diesem Fall wie folgt: Für 2 0x x a b : 2 2 xa xx x x a  ¬a ­ž ­ ž ­ž ­­žŸ ® 2 xb xx x  a (2.49) 25 Für 2 0x x a p : 2 1 1 1 xa xx x  a  2 2 21 xb x x a a (2.50) Bild 2.2.3 zeigt ein Beispiel für eine typische Doppel-Delta-Funktion für den Mischungs- bruch: Bild 2.2.3 Doppel-Delta-Funktion für den Mischungsbruch f Sind neben der Normierungsbedingung und den Gleichungen für die ersten beiden Momente noch weitere Bestimmungsgleichungen bekannt, so kann die Doppel-Delta-Funktion durch weitere Peaks erweitert werden. Somit können zum Beispiel Größen wie die Schiefe einer PDF mit berücksichtigt werden. Zu diesem Thema sind jedoch bisher keine Literaturstellen bekannt. 26 b) Abgeschnittene Gauß-Verteilung Eine weitere oft verwandte und wohl die bekannteste Form für die Bestimmung einer Wahr- scheinlichkeitsdichte-Funktion ist die Gauß-Verteilung. Viele Phänomene in der Natur las- sen sich durch die Gauß-Annahme sehr gut beschreiben. Die Bestimmungsgleichung hierfür lautet: 2 2 1( ) exp 2 2G P G NG T Q T  ¬ ­ž ­ ž ­ž ­žŸ ® (2.51) Die Gauß-Verteilung hat einen Definitionsbereich von -∞ bis +∞. Da die Einflussparameter der Chemie einen endlichen Gültigkeitsbereich [a,b] besitzen, wird in der Verbrennungs- Simulation eine Sonderform der Gauß-Verteilung benutzt, und zwar die sogenannte abge- schnittene Gauß-Verteilung. Hierbei ist der Bereich der Gauß-Verteilung, der außerhalb des Definitionsbereichs der Einflussparameter liegt, durch Delta-Peaks an den Rändern ersetzt. Damit erweitert sich Gleichung (2.51) wie folgt: 221( ) exp2 2 aaG bP A B G NG E G G E G GT Q T  ¬ ­ž ­   ž ­ž ­žŸ ® (2.52) mit a GA P dG G d  ¨ und G b B P dG G d  ¨ (2.53) Bild 2.2.4 zeigt eine abgeschnittene Gauß-Verteilung für den Mischungsbruch. Diese Wahr- scheinlichkeitsdichtefunktion ist von vier Parametern abhängig: dem Mittelwert µ, der Standardabweichung σ und den beiden Wichtungsparametern der Dirac-Peaks A und B. 27 Bild 2.2.4 Abgeschnittene Gauß-Verteilung für den Mischungsbruch f Zur Bestimmung dieser Parameter sind vier Gleichungen erforderlich. Dies sind neben der Bedingung für Mittelwert und Standardabweichung wiederum die Normierungsbedingung und zusätzlich die Definition des Gaußschen Fehlerintegrals. Das Einsetzen von Gleichung (2.52) in die Momentengleichungen (2.44) und (2.45) führt zunächst zu 1 2 2 0 exp 2 2 fff B df N T Q T  ¬ ­ž ­ ž ­ž ­žŸ ®¨ (2.54) und 1 222 2 2 0 exp 2 2 fff B df f N T Q T  ¬ ­ža ­  ž ­ž ­žŸ ®¨ (2.55) Mit der Substitution fz NT  und a az NT  , b bz NT  (2.56) 28 und dem Gaußschen Fehlerintegral 21 exp 2 2 x sx dsQ d  ¬­ž'   ­ž ­­žŸ ®¨ (2.57) ergibt sich: 221 1 exp exp 2 2 2 a b a b zzf z z TN N Q   ¯ ¬ ¬ ­­ žž¡ °  '  ¸'    ­­ žž ­ ­žŸ ®¡ °Ÿ ®¢ ± (2.58) und 2 2 2 2 2 22 22 2 1 1 exp exp 2 2 2 2 exp exp 2 2 2 a b a b a b a b f z z zzz z zz f T N T N T Q TN Q a   ¸'  ¸'   ¯ ¬ ¬ ­­ žž¡ °    ­­ žž ­ ­žŸ ®¡ °Ÿ ®¢ ±   ¯ ¬ ¬ ­­ žž¡ °    ­­ žž ­ ­žŸ ®¡ °Ÿ ®¢ ± (2.59) Aus diesen beiden transzendenten Gleichungen werden die beiden Parameter µ und σ be- stimmt. Dies erfolgt mittels eines auf Intervalschachtelung basierenden Iterationsschemas (Press, 1989). Diese Berechnung nimmt etwa ein Drittel der gesamten für die Integration über die PDF erforderlichen Rechenzeit in Anspruch. Soll die Integration online, d.h. wäh- rend der Strömungssimulation durchgeführt werden, so ist es sinnvoll, die Abbildung Mit- telwert und Standardabweichung der zu integrierenden Größe zu den Parametern µ und σ der Gauß-Verteilung im Preprozessing zu berechnen und in Tabellen abzulegen. Für die Beträge der Dirac-Peaks ergeben sich aus der Normierungsbedingung aA z ' und 1 bB z ' (2.60) Damit kann die PDF durch die Angabe von Mittelwert und Varianz und der Angabe der Grenzen des Gültigkeitsbereichs eindeutig bestimmt werden. 29 c) β-PDF Eine dritte, sehr oft angewandte Form einer Wahrscheinlichkeitsdichtefunktion ist die β- PDF. Sie ist vom Rechenzeitbedarf her etwas günstiger als die Gauß-Verteilung, da sich ihre beiden Parameter an Hand folgender einfachen Abbildung direkt aus Mittelwert und Vari- anz der Einflussgrößen ergeben: 2 1 1 G GB G G  ¬ ­ž ­ž ¸  ­ž ­ž a ­žŸ ® und 211 1G GC G G  ¬ ­ž ­ž  ¸  ­ž ­ž a ­žŸ ® (2.61) Die Bestimmungsgleichung der β-PDF lautet: 11 1 11 0 1 1 P d CB CB G GG G G G    ¨ (2.62) Diese PDF ist im Bereich 0 ≤ Φ ≤ 1 definiert und kann verschiedene typische Formen an- nehmen, die in Bild 2.2.5 zu erkennen sind. Im Gegensatz zur Gauß-Verteilung kann die β- PDF auch asymmetrische und bimodale Verteilungen abbilden. Allerdings liegen die dann entstehenden Maxima immer am Rand des Definitionsbereichs. Die bimodale Verteilung (α=0, β=0) kommt bei Verbrennungen nie vor und ist daher im Rahmen dieser Arbeit nicht interessant. Numerische Schwierigkeiten ergeben sich bei Eingabe von Mittelwerten, die sich nahe an den Definitionsgrenzen befinden. Hier sind Fallunterscheidungen notwendig, um die Bestimmung der PDF numerisch robust zu halten. 30 Bild 2.2.5 Mögliche Formen der β-PDF [ Warnatz et al.,1997 ] Die β-PDF liefert im Grunde genommen ähnliche Ergebnisse wie die Gauß-PDF. Werden bei der Bestimmung der Gauß-Verteilung die Parameter µ und σ aus Tabellen gelesen, dann hat die β-PDF auch keine rechentechnischen Vorteile mehr durch kürzere Rechenzeiten. Daher wird die β-PDF im weiteren Verlauf der Arbeit nicht mehr weiter in die Betrachtun- gen mit einbezogen. 2.2.5 Mehrdimensionalität Die gerade vorgestellten Formen von Wahrscheinlichkeitsdichtefunktionen sind hier nur als eindimensionale Funktion gezeigt worden. In der Realität liegen jedoch mehrdimensionale Verteilungen vor, da jeder Einflussparameter nicht nur eine eigene Schwankungscharakte- ristik aufweist, sondern diese auch in Abhängigkeit von anderen Parametern ändert. Diese mehrdimensionale Form der PDF erhält man bei der Lösung der Transportgleichung für die komplette PDF [Pope,1976]. Löst man aus Gründen des Rechenzeitbedarfs jedoch nur die Gleichungen für die ersten beiden Momente, so fehlt die Information über mögliche Korre- lationen der einzelnen Parameter untereinander. Hier wird daher die Annahme der statistischen Unabhängigkeit getroffen werden. Statistisch unabhängig bedeutet, dass die Einflussparameter unkorreliert sind. Ob zwei Einflussparame- 31 ter x und y miteinander in Korrelation stehen, zeigt der sogenannte Korrelationskoeffizient rxy [Bronstein & Semendjajew, 1989]: 2 2 i i xy i xy x y i i i i x x y ys r s s x x y y    ¸  ¸  œ œ œ (2.63) Dabei ist sxy die Kovarianz zweier Zufallsgrößen x mit y. sx ist die Varianz von x und sy ist die Varianz von y. Der Korrelationskoeffizient sagt aus, inwieweit sich eine größere Anzahl von Zahlenpaaren durch eine Gerade annähern lässt [Dubbel,1981]. Strebt er gegen plus oder minus eins, so ergibt sich eine Geradengleichung, die den Zusammenhang zwischen den Zahlenpaaren beschreibt. Strebt der Korrelationskoeffizient gegen Null, so sind die Funktionen vollkom- men unkorreliert und es besteht kein Zusammenhang zwischen den Zahlenpaaren. Wenn keine Korrelation zwischen den Wertepaaren der Einflussparameter vorliegt und die Einflussparameter normalverteilt sind, dann sind sie als statistisch unabhängig anzusehen [Bronstein & Semendjajew, 1989]. Im Fall der statistischen Unabhängigkeit kann die mehrdimensionale Wahrscheinlichkeits- dichtefunktion P(ϕ1, ϕ2, ..., ϕn) durch das Produkt einzelner marginaler PDF’s ersetzt wer- den. Definition der statistischen Unabhängigkeit: ( ) ( )i i P PK K (2.64) wobei der Vektor K die statistisch unabhängigen Zufallsgrößen und P(K) deren Erwar- tungswert beinhaltet. Gerade in der Nähe des chemischen Gleichgewichtszustandes stimmt diese Annahme nicht mehr. Die Theorie der reduzierten Mechanismen zeigt, dass die Anzahl der unabhängigen Variablen von Reaktionsmechanismen stetig mit Annäherung an den Gleichgewichtszustand abnimmt, bis schließlich noch die reine Mischungsabhängigkeit verbleibt. Daher ist es not- wendig zu wissen, wie groß der Fehler ist, der aus der Annahme der statistischen Unabhän- gigkeit resultiert. Dies wird in Kapitel 6 näher erläutert werden. Soll eine eventuelle statistische Abhängigkeit der Parameter mit berücksichtigt werden ohne die komplette PDF zu transportieren, so könnte es eine Kompromisslösung geben, die bisher noch nicht näher untersucht wurde. Die Transportgleichungen für die Querkorrelationen der Einflussparameter sind nur unwesentlich komplizierter als die der zweiten Momente. Trans- portiert man diese mit, ergeben sich weitere Bestimmungsgleichungen für die Bestimmung von Mehr-Delta-PDF’s. Jede neue Bestimmungsgleichung macht es möglich, weitere Peaks zu setzen, die dann Phänomene wie die Schiefe einer PDF oder Ähnliches simulieren kön- nen. 32 2.3 Strömungslöser Alle in dieser Arbeit gezeigten Strömungsberechnung wurde an Hand des DLR-eigenen Strömungslösers TRUST durchgeführt. TRUST löst die instationären Navier-Stokes- Gleichungen sowohl für den kompressiblen als auch für den inkompressiblen Bereich, wo- durch der Löser für viele praktische Anwendungen geeignet ist. Die meisten Strömungslöser bedienen sich der Eulerschen Formulierung der Navier-Stokes- Gleichungen. Diese Formulierung hat sich als vorteilhaft für die Beschreibung stationärer Probleme herausgestellt. Da TRUST ursprünglich für die Berechnung instationärer Strö- mungen mit Verbrennung und Spray entwickelt wurde, sind die Transportgleichungen in der Langrangeschen Schreibweise formuliert. Dies hat den Vorteil, dass während der Lösung der Transportgleichungen die schwierig zu modellierenden konvektiven Terme nicht expli- zit auftreten. So erweist sich der Löser auch für die gekoppelte Berechnung vieler einzelner physikalischer Phänomene als sehr robust. Die Robustheit und das Konvergenzverhalten des Lösers ist gerade bei der Berechnung von Strömung mit Verbrennung sehr wichtig, da hier große Dichtegradienten auftreten. TRUST beinhaltet eine ganze Reihe von physikalischen Modellen, wie beispielsweise Strahlungsmodelle, Modelle für Spray-Verdunstung und -Verbrennung, turbulente Gleich- gewichtschemie-Chemie und diverse Turbulenz-Modelle. Jedes dieser einzelnen Modelle liefert einen eigenen Zeitschritt. TRUST wählt nun aus diesen Teilzeitschritten den kleins- ten Zeitschritt aus und bestimmt hiernach eine CFL-Zahl [Zurmühl,1965]. Das damit ver- bundene Stabilitätskriterium bestimmt nun automatisch in jeder Rechenzelle eine lineare Wichtung zwischen expliziter und impliziter Lösung der Transportgleichungen. In den Be- reichen, wo die physikalischen Prozesse einen großen Zeitschritt erlauben wird implizit ge- rechnet. Der implizite Löser bedient sich eines Newten-CG-Verfahrens und einer Quasi- Second-Order Diskretisierung [KIVA, 1998]. Die Netzgenerierung erfolgte mit den bekannten kommerziellen Programmen ICEM-CFD und IDEAS. Die Netze sind unstrukturiert und basieren auf Hexaedern. Hexaeder haben gegenüber Tetraedern den Vorteil, dass in Nähe der Berandung ihre Flächen oft annähernd parallel oder normal zur Strömung liegen. Dies begünstigt das Konvergenzverhalten des Lösers, was sich nicht nur positiv auf die Genauigkeit der Lösung auswirkt, sondern auch eine Reduzierung der Rechenzeiten zur Folge hat. Ferner benötigen Hexaeder deutlich we- niger Elemente zum Auflösen einer Geometrie als Tetraeder. Der Nachteil der Hexaeder liegt in der Netzgenerierung. Für sehr komplexe Geometrien ist es mittlerweile möglich, quasi auf Knopfdruck Tetraedernetze zu generieren. Unstrukturierte Hexaedernetze in kom- plexen dreidimensionalen Geometrien zu erstellen, ist jedoch oft sehr schwierig und bedarf sehr viel Erfahrung des Bearbeiters. Neuere Netzgenerierungsprogramme wie beispielswei- se MARC [MSC] können schon automatische Hexaeder-Netze erzeugen. Bei komplexen Geometrien werden hierbei jedoch sehr viele Zellen erzeugt, die zum Teil sehr ungünstige Formen annehmen können. Dieses Problem wird erst dann gelöst sein, wenn Strömungslö- ser lokale Verfeinerungen mit hängenden Knoten verarbeiten können. Die Netzgenerierung ist nach wie vor eins der zentralen Probleme bei der Simulation von technischen Brenn- kammern. Unzulänglichkeiten des Netzes können mit den besten physikalischen Modellen nicht mehr ausgeglichen werden. Bei allen gezeigten Berechnungen wurden als Randbedingungen am Einlass Geschwindig- keitsrandbedingungen gesetzt. Für die Bestimmung des Druckes im Auslass wurde sich ei- ner eindimensionalen Wellengleichung bedient [Rudy&Strikwerda,1980]. An der Wand wurde ein turbulentes Wandgesetz benutzt [KIVA,1998]. Als Turbulenzmodell stand das 33 modifizierte k-ε-Modell (MKE) [Hirsch,1996] zur Verfügung, welches insbesondere für stark verdrallte Strömungen gute Ergebnisse liefert.