Modellierung reaktiver Strömungen unter Verwendung automatisch reduzierter Reaktionsmechanismen [PDF]

Zitiervorschau

INAUGURAL-DISSERTATION zur Erlangung der Doktorwurde der Naturwissenschaftlich - Mathematischen Gesamtfakultat der Ruprecht - Karls - Universitat Heidelberg

vorgelegt von Diplom-Physiker Dietmar Schmidt geboren in Geislingen an der Steige Tag der mundlichen Prufung: 09.07.1996

Thema

Modellierung reaktiver Stromungen unter Verwendung automatisch reduzierter Reaktionsmechanismen

Gutachter: Prof. Dr. Jurgen Warnatz Prof. Dr. Jurgen Wolfrum

Inhaltsverzeichnis 1 Einleitung 2 Mathematische Modellierung reaktiver Stromungen

5 8

2.1 Einfuhrung : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : 8 2.2 Erhaltungsgleichungen und empirische Gesetze : : : : : : : : : : : : : : : : : : : : : 9

3 Automatische Reduktion von Reaktionsmechanismen

3.1 Einfuhrung : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : 3.2 Dynamik chemischer Reaktionssysteme : : : : : : : : : : : : : : : : : : : : : 3.2.1 Allgemeine chemisch reagierende Systeme : : : : : : : : : : : : : : : 3.2.2 Homogene Reaktionssysteme : : : : : : : : : : : : : : : : : : : : : : 3.2.3 Begrisbildung : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : 3.2.4 Chemische Kinetik am Beispiel der Methanoxidation : : : : : : : : : 3.2.5 Erhaltungsgroen, partielle Gleichgewichte und stationare Zustande 3.2.6 Eigenwertanalysen chemischer Systeme : : : : : : : : : : : : : : : : : 3.3 Intrinsische niedrig-dimensionale Mannigfaltigkeiten : : : : : : : : : : : : : 3.3.1 Einfuhrung : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : 3.3.2 Mathematische Formulierung der ILDM : : : : : : : : : : : : : : : : 3.3.3 Numerische Berechnung niedrig-dimensionaler Mannigfaltigkeiten : : 3.3.4 Kopplung von chemischer Reaktion und Transport : : : : : : : : : :

4 Anwendungen

4.1 Einfuhrung : : : : : : : : : : : : : : : : : : : : : : : : : : : 4.2 Homogene Reaktionssysteme : : : : : : : : : : : : : : : : : 4.2.1 Einfuhrung : : : : : : : : : : : : : : : : : : : : : : : 4.2.2 Erhaltungsgleichungen und Losungsansatz : : : : : : 4.2.3 Geschlossenes homogenes System : : : : : : : : : : : 4.2.4 Oener Reaktor und Zundung : : : : : : : : : : : : 4.2.5 Oener Reaktor, Mischung mit unverbranntem Gas 4.3 Simulation eindimensionaler laminarer Flammen : : : : : : 4.3.1 Einfuhrung : : : : : : : : : : : : : : : : : : : : : : : 4.3.2 Erhaltungsgleichungen und Losungsansatz : : : : : : 4.3.3 Laminare CH4 -Vormisch ammen : : : : : : : : : : : 4.3.4 Laminare C7H16 Vormisch ammen : : : : : : : : : : 4.3.5 Laminare Methanol-Flamme : : : : : : : : : : : : : 4.3.6 Zundung eines Methanol-Luft-Gemisches : : : : : : 4.4 Simulation von Hyperschallstromungen : : : : : : : : : : : : 4.4.1 Einfuhrung : : : : : : : : : : : : : : : : : : : : : : : 4.4.2 Erhaltungsgleichungen und Losungsansatz : : : : : : 4.4.3 Stromung um einen Wiedereintrittskorper : : : : : : 3

: : : : : : : : : : : : : : : : : :

: : : : : : : : : : : : : : : : : :

: : : : : : : : : : : : : : : : : :

: : : : : : : : : : : : : : : : : :

: : : : : : : : : : : : : : : : : :

: : : : : : : : : : : : : : : : : :

: : : : : : : : : : : : : : : : : :

: : : : : : : : : : : : : : : : : :

: : : : : : : : : : : : : : : : : :

: : : : : : : : : : : : : : : : : : : : : : : : : : : : : : :

: : : : : : : : : : : : : : : : : : : : : : : : : : : : : : :

: : : : : : : : : : : : : : : : : : : : : : : : : : : : : : :

: : : : : : : : : : : : : : : : : : : : : : : : : : : : : : :

: : : : : : : : : : : : : : : : : : : : : : : : : : : : : : :

13

13 14 14 14 14 16 25 28 31 31 33 49 56

65

65 67 67 67 68 69 71 72 72 72 74 84 85 88 93 93 93 95

4

5 Zusammenfassung 6 Symbolverzeichnis

INHALTSVERZEICHNIS

99 108

Kapitel 1

Einleitung Chemisch reagierende Stromungen sind charakterisiert durch das komplexe Zusammenspiel von chemischer Reaktion, Konvektion und molekularem Transport. Vor allem auf dem Feld der detaillierten Reaktionsmechanismen, aber auch der Beschreibung von Transportprozessen, sind in den letzten Jahren groe Anstrengungen unternommen worden, um z.B. Verbrennnungsprozesse genauer zu beschreiben. Die mathematische Modellierung hat dabei inzwischen einen Stellenwert eingenommen, der es erlaubt, ein tiefergehendes Verstandnis der komplizierten Zusammenhange zu erlangen. Dadurch ist eine systematische Untersuchung der komplexen Vorgange in typischen Verbrennungssystemen mit Hilfe des Computers moglich. Insbesondere auf dem Gebiet der laminaren Stromungen ist der Kenntnisstand inzwischen so gut, da die auf dem Computer gewonnenen Ergebnisse sehr gut mit Experimenten ubereinstimmen [1{9]. Da aber die verwendeten Reaktionsund Transportmodelle immer komplexer, und damit auch immer realitatsnaher, werden, ist man allerdings immer noch auf sehr vereinfachte Geometrien beschrankt. So sind Simulationen in zwei Raumrichtungen, die detaillierte Reaktionsmechanismen und detaillierten Transport berucksichtigen, auf heutigen Rechnern kaum moglich [3, 10{12]. Geht man noch einen Schritt weiter und untersucht die in der Praxis hau ger vorkommenden turbulenten Stromungen, so stot man schon sehr schnell an die Grenzen des Machbaren. Insbesondere die Kopplung an die chemische Reaktionen bereitet doch erhebliche Schwierigkeiten [13{15]. Der Grund fur die Probleme numerischer Simulationen mit Hilfe detaillierter Reaktionsmechanismen ist, da fur jede chemische Spezies (um so mehr, je detaillierter und besser der Reaktionsmechanismus ist) eine Erhaltungsgleichung gelost werden mu. Es kommt erschwerend hinzu, da die unterschiedlichen Reaktionen mit Zeitskalen ablaufen, die sich uber mehr als 10 Zehnerpotenzen erstrecken. Das erfordert fur die numerische Simulation implizite Losungsmethoden, die sehr Rechenzeit-intensiv sind. Bei technischen Verbrennungsvorgangen liegt ein komplexes Gemisch hoherer Kohlenwasserstoe vor. Dabei erhalt man aber leicht mehr als 100 chemische Spezies mit uber 1000 Reaktionen. Das ist auf realitatsnahen Gittern auf Jahre hinaus rechnerisch nicht zu bewaltigen. Deshalb hat in den letzten Jahren die Entwicklung vereinfachter Kinetiken an Bedeutung gewonnen. Dabei mu naturlich gewahrleistet sein, da die vereinfachte Kinetik die Realitat ausreichend genau beschreibt. Die ersten Arbeiten auf dem Gebiet der Mechanismenreduktion gehen dabei auf Bodenstein [16] und Semenov [17] zuruck. Sie beobachteten, da schnelle Reaktionsprozesse zu quasistationaren Zustanden und partiellen Gleichgewichten fuhrten. Diese Idee wurde vielerorts aufgegrien, um detaillierte Reaktionsmechanismen zu vereinfachen [14, 18{23]. Dabei werden zwei Ziele gleichzeitig verfolgt: zum einen die Reduktion der Anzahl der Spezies fur die man die Erhaltungsgleichungen losen mu (jede Annahme von Quasistationaritat und partiellem Gleichgewicht erniedrigt die Dimension des Gleichungssystems um eins, da die entsprechende Dierentialgleichung in eine algebraische Gleichung uberfuhrt wird), und zum zweiten damit verbunden eine Reduktion der Steifheit des Gleichungssystems. Um die heutigen komplexen Reaktionsmechanismen mit Hil5

6

KAPITEL 1. EINLEITUNG

fe der (globalen) Annahmen der Quasistationaritat und partieller Gleichgewichte zu vereinfachen, ist allerdings ein sehr hoher Kenntnisstand der ablaufenden Prozesse und der chemischen Kinetik notwendig. Auerdem ist der Arbeitsaufwand sehr hoch, da fur jedes neue Problem der detaillierte Reaktionsmechanismus neu untersucht werden mu. So andern sich z.B. bei nicht vorgemischter Verbrennung die Reaktionspfade beim U bergang von einem fetten Bereich in einen mageren. Fur jeden Reaktionspfad mu aber der Reaktionsmechanismus erneut vereinfacht werden. Die auf diesen globalen Annahmen beruhenden vereinfachten Mechanismen sind also nur begrenzt einsetzbar. Dies ist besonders schwerwiegend, wenn man technisch relevante Verbrennungsprozesse, die in der Regel turbulenter Natur sind, beschreiben will. Dort hat man den kompletten Parameterbereich zu beschreiben, d.h. sehr fette und magere Bedingungen in stark unterschiedlichen Temperaturenbereichen. Die meisten Probleme bei der Reduktion der Reaktionsmechanismen lassen sich vermeiden, wenn man den detaillierten Reaktionsmechanismus einer mathematischen Analyse unterzieht. Lam und Goussis haben die Idee der singularen Storungstheorie zur Identi kation partieller Gleichgewichte und quasistationaren Zustanden auf Reaktionsmechanismen angewandt [24,25]. Diese Informationen konnen dann zur Entwicklung reduzierter Reaktionsmechanismen benutzt werden [26]. Das Verfahren kann auch direkt zur Losung der chemischen Reaktionsgleichungen angewendet werden. Damit wird die Losung des Dierentialgleichungssystems sehr ezient gelost, da man durch die mathematische Analyse automatisch die schnellen Zeitskalen identi zieren und entkoppeln kann. Allerdings bietet die Methode keine Moglichkeit, auch die Dimension des Gleichungssystems zu reduzieren. Fur praktische Anwendungen ist es notwendig, die Zahl der Freiheitsgrade zu verringern, da es auf heutigen Rechneranlagen unmoglich ist, bei anwendungsorientierten Geometrien mehrere hundert Erhaltungsgleichungen zu losen. Deshalb ist das Ziel bei der Entwicklung reduzierter Reaktionsmechanismen, sowohl die Steifheit als auch die Zahl der reaktiven Variablen zu verkleinern, ohne einen zu groen Verlust an Genauigkeit in Kauf nehmen zu mussen. In der vorliegenden Arbeit soll nun ein Verfahren vorgestellt und auf komplexe Reaktionsmechanismen angewandt werden, das die obigen Bedingungen erfullt [27{33]. Das Prinzip des Verfahrens beruht auf einer mathematischen Analyse des Reaktionssystems, die Informationen uber die Zeitskalen der Reaktionsprozesse liefert. Dabei erkennt man, da die chemischen Reaktionen unterschiedlich schnell sind (typischerweise von 10,10s bis 1 s). Stromung und Transport laufen in typischen Verbrennungssystemen dagegen mit Zeitskalen im Bereich von 100 s bis 10 ms ab. Zwischen den sehr schnell ablaufenden chemischen Prozessen und den physikalischen Zeitskalen des molekularen Transports ndet keine direkte Kopplung statt. D.h., die sehr schnellen chemische Vorgange fuhren zur Ausbildung von quasistationaren Zustanden und partiellen Gleichgewichten. Fur diese Prozesse lassen sich, basierend auf der mathematischen Analyse, die Gleichgewichtsbedingungen nden, so da die schnellsten Prozesse vom Dierentialgleichungssystem entkoppelt werden konnen. Die einzige Information, die das Verfahren benotigt, sind dabei die Anzahl der Gleichgewichtsbedingungen, oder in anderen Worten, der Grad der Vereinfachung. Die komplette Analyse basiert auf einem mathematischen Algorithmus. Dadurch erreicht man, da fur die Durchfuhrung der Vereinfachung eines Reaktionsmechanismus kein Fall-zu-Fall-Verstandnis mehr notwendig ist. Auerdem erfordert die Methode keine globalen Vorausannahmen, so da ein und derselbe reduzierte Mechanismus uber einen weiten Parametersatz der Stochiometrie und der Temperatur Gultigkeit hat. Das hier angedeutete Verfahren kommt in dieser Arbeit zur Anwendung und wird in Kap. 3.3 eingehend vorgestellt. Zuvor werden in Kap. 2 die mathematischen Modelle reaktiver Stromungen behandelt. Da die Reduktionsmethode auf der Analyse der chemischen Kinetik basiert, wird in Kap. 3.2 auf die allgemeine chemische Kinetik eingegangen. Inbesondere wird auf die Kinetik der Methanoxidation vorgestellt, da diese den Hauptbestandteil der Anwendungen in dieser Arbeit reprasentiert. In diesem Zusammenhang wird auf eine anschauliche Art und Weise gezeigt, da die chemische Kinetik im allgemeinen von einer vieldimensionalen Beschreibung auf ein Beschrei-

7 bung durch wenige Variable projiziert werden kann. Die vereinfachte Darstellung entspricht einer Bewegung in niedrig-dimensionalen Unterraumen, sogenannten intrinsischen niedrig-dimensionalen Mannigfaltigkeiten (ILDM), des Zustandsraumes, der durch die Anzahl aller im Reaktionssystem vorkommenden Spezies charakterisiert ist. Die niedrig-dimensionalen Mannigfaltigkieten, werden dagegen durch wenige Reaktionsfortschrittsvariablen beschrieben. Das mathematische Modell und die Ableitung der Methode der intrinsischen niedrig-dimensionalen Mannigfaltigkeiten wird in Kap. 3.3 vorgestellt. Dabei zeigt sich, da besonderer Wert auf die Kopplung der reduzierten chemischen Dynamik mit der Stromung und molekularem Transport gelegt werden mu. Die Validierung der Methode und Anwendung auf komplexe Reaktionsmechanismen sind in Kap. 4 zusammengestellt. Hier werden vereinfachte Mechanismen der Methan-, Methanol-, und Heptanoxidation und der Dissoziation von Luft auf technisch relevante Simulationsbeispiele angewandt. Diese sind die Simulation perfekt durchmischter Reaktoren und eindimensionale Flammen (sowohl stationare Flammen als auch instationare Flammenausbreitung). Da die Methode nicht nur auf Verbrennungssysteme anwendbar ist, wird bei der Simulation einer Hyperschallumstromung eines Raumgleiters beim Wiedereintritt in die Erdatmosphare gezeigt.

Kapitel 2

Mathematische Modellierung reaktiver Stromungen 2.1 Einfuhrung Die Dynamik chemisch reaktiver Systeme wird durch die mathematische Modellierung von Erhaltungsgleichungen beschrieben. Die zugrunde liegenden Gleichungen fur die Variablen der Gesamtmasse, der verschiedenen Teilchenmassen, des Impulses und der Energie bilden ein gekoppeltes nichtlineares Dierentialgleichungssystem, welches mit numerischen Methoden gelost werden kann. Erhaltungsgleichungen lassen sich makroskopisch (siehe Kap. 2.2) oder mikroskopisch [34] fur die extensiven Variablen herleiten. Diese extensiven Variable sind charakterisiert durch ihre Stomengenabhangigkeit (z.B. Gesamtmasse, Energie und Impuls). Intensive Variablen, wie Temperatur oder Druck, sind dagegen stomengenunabhangig. Sf

n

φf

qf

θ δθ

Abbildung 2.1: A nderung einer extensiven Groe in einem Volumenelement. Ein Volumenelement  sei durch seine Ober ache @ begrenzt. In  andert sich eine extensive Groe F bzw. seine Dichte f (~r; t) durch verschiedene Prozesse (vergleiche Abb. 2.1):  Flu durch die Ober ache @ (Diusion, Warmeleitung, Reibung) ~f ~ndS , wobei ~f = Stromdichte, ~n = Normalenvektor auf der Ober ache @ und dS = dierentielles Ober achenelement bedeutet.  Produktion im Volumenelement (chemische Produktion) qf .  Fernwirkung (Strahlung, Gravitation) Sf . Bilanziert man die verschiedenen Prozesse auf, erhalt man fur die zeitliche A nderung der Dichte von F @f + div ~  = q + S : (2.1) f f f @t 8

2.2. ERHALTUNGSGLEICHUNGEN UND EMPIRISCHE GESETZE

9

2.2 Erhaltungsgleichungen und empirische Gesetze Kontinuitatsgleichung

Die Erhaltungsgleichung der Gesamtmasse wird durch die Kontinuitatsgleichung beschrieben,

@ + div(~v) = 0 : (2.2) @t Die Erhaltungsgroe ist die Dichte  und ~v ist der Geschwindigkeitsvektor des stromenden Mediums. Da chemische Reaktionen weder Gesamtmasse erzeugen noch sie vernichten, gibt es in der Kontinuitatsgleichung keinen Quellterm.

Teilchenerhaltung und chemische Reaktion

Die Dichte  derPMischung setzt sich aus den Partialdichten i aller Spezies i in der Mischung zusammen:  = i i wi, mit wi = Massenbruch der Spezies i. Fur die Partialdichten erhalt man die Erhaltungsgleichung @i + div( ~v ) = M ! : (2.3) i i i i @t

Mi!i bezeichnet den Produktionsterm der Spezies i infolge chemischer Reaktionen mit Mi = molare Masse, !i = chemischer Umsatz der Spezies i und ~vi = Geschwindigkeit der Spezies i. Sie ist durch einen konvektiven Anteil und der Diusionsgeschwindigkeit gegeben: ~vi = ~v + ~viD . Fur eine

Multikomponentenmischung hangt die Diusionsgeschwindigkeit einer Spezies von den Gradienten der anderen Spezieskonzentrationen ab. Da deren Berucksichtigung mit einem groen Rechenaufwand verbunden ist, nahert man dies meistens durch eine Diusion der Spezies i in die Mischung an. Die Diusionsgeschwindigkeit der Spezies i hangt in diesem Fall nur noch vom Gradienten des Molenbruchs der Spezies i ab [35]. Detaillierte Transportmodelle berucksichtigen zudem noch die Thermodiusion der Teilchen. Der Diusionsstrom ergibt sich somit zu

~ji = ~jiD + ~jiT i D rx , DiT rT : = , w x i i T i

(2.4)

Dabei ist wi = Massenbruch, xi = Molenbruch der Spezies i, T = Temperatur. Die Transportkoezienten Di und DiT folgen aus der kinetischen Gastheorie verdunnter Gase [34,36]. Fur weite Bereiche lassen sie sich als Polynome der Temperatur darstellen. Die Teilchenerhaltungsgleichung lat sich zusammenfassen zu

@i + div( ~v) + div~j  = M ! : i i i i @t

(2.5)

Die Reaktionsgeschwindigkeiten der Spezies erhalt man durch Bilanzierung aller chemischen Reaktionen, 0 1

!i =

ncr X @ l=1

(^ail , ail ) kl

ns X

j =1

caj jl A :

Dabei sind ncr = Anzahl der chemischen Reaktionen im System, a^il ; ail = die stochiometrischen Koezienten der Edukte und Produkte, kl = Geschwindigkeitskoezienten des Reaktionssystems der Spezies Xi in der Form ns X i=1

kl

ail Xi !

ns X i=1

^ail Xi

l = 1; : : :; ncr :

10

 KAPITEL 2. MATHEMATISCHE MODELLIERUNG REAKTIVER STROMUNGEN

Der Geschwindigkeitskoezient kl wird durch einen erweitertem Arrheniusansatz ausgedruckt,

kl(T ) = Al T l exp

 Eal  R T ;

Gas Faktor, l =

wobei T = Temperatur, Al = praexponentieller Temperaturexponent, Eal = Aktivierungsenergie der Reaktion l und RGas = allgemeine ideale Gaskonstante sind.

Impulserhaltung Die Impulserhaltungsgleichungen sind gegeben durch die Navier-Stokes-Gleichungen

@~v + div(~v ~v ) + div(p) = ~g : (2.6) @t ~v ist dabei das Geschwindigkeitsfeld und ~g die Erdbeschleunigung (die im folgenden vernachlassigt wird). p ist der Drucktensor, der sich aus dem hydrostatischem Druck pE und dem viskosen Spannungstensor  zusammensetzt (E = Einheitstensor). Empirisch ergibt sich fur den Spannungstensor nach dem Newtonschen Spannungsgesetz [34,36]





 = , grad(~v ) + grad(~v )T , 23 div(~v )E :

(2.7)

Ein Quellterm des Impulses existiert nicht.  ist der Viskositatskoezient, der analog zum Diusionskoezient als Polynom der Temperatur dargestellt wird.

Energieerhaltung Die Erhaltungsgleichung der inneren Energie lautet

@u + div(u~v ) + div ~j  + p : grad(~v ) = q : q r @t ~jq steht fur die Warmestromdichte, p fur den Drucktensor und qr fur einen Strahlungsquellterm.

In typischen Verbrennungssystemen wird allerdings meist eine Erhaltungsgleichung fur die spezi sche Enthalpie h, die durch h = u + p (mit p = Druck) gegeben ist, verwendet. Die Erhaltungsgleichung der spezi schen Enthalpie lautet

@h , @p + div h~v + ~j  + p : grad(~v) , div(p~v ) = q : q r @t @t

(2.8)

Die Warmestromdichte ergibt sich empirisch aus dem Fourierschem Warmeleitungsgesetz [34,36] und dem diusiven Transport der spezi schen Enthalpie. Vernachlassigt man die Druckabhangigkeit, ist sie gegeben durch ns ~jq = , grad(T ) + X hi~ji : (2.9) i=1

 ist der Warmeleitfahigkeitskoezient und hi die spezi sche Enthalpie der Spezies i. Sie ergibt sich durch Integration der spezi schen Warmekapazitat der Spezies i uber die Temperatur

ZT ,
0 hi (T ) = hi + cp;i T 0 dT 0 T0

:

h0i ist die Bildungsenthalpie der Spezies i zur Temperatur T0 (in der Regel gilt T0 = 298 K). cp;i ist die spezi sche Warmekapazitat der Spezies i bei konstantem Druck p. Fur die spezi sche Warmekapazitat cp und Enthalpie h der Mischung gelten cp =

ns X i=1

wi cp;i

und

h=

ns X i=1

wi hi :

2.2. ERHALTUNGSGLEICHUNGEN UND EMPIRISCHE GESETZE

11

Ideales Gasgesetz Um das vorgestellte Dierentialgleichungssystem zu schlieen, benotigt man noch eine Beziehung, die den Druck mit der Temperatur und der Dichte verknupft. In weiten Bereichen der technischen Anwendungen kann die Gasmischung durch das ideale Gasgesetz angenahert werden, da in typischen Verbrennungssystemen hohe Temperaturen bei maigem Druck vorherrschen. Das ideale Gasgesetz lautet p = RGas T (2.10)

M mit M = mittlere molare Masse der Mischung, RGas = Gaskonstante,  = Druck und T = Temperatur.

Primitive Form der Erhaltungsgleichungen

Die aufgestellten Erhaltungsgleichungen fur die Dichte , Partialdichten der Spezies i , des Impulses ~v und der spezi schen Enthalpie h sind jeweils fur die konservativen Variablen formuliert. D.h., die unabhangige Variable entspricht jeweils der Dichte der Erhaltungsgroe. Oft ist es aber von Vorteil, die Dierentialgleichungen fur die primitiven Variablen Dichte , Massenbruche wi, Geschwindigkeiten ~v und Temperatur T zu formulieren. In vielen Fallen ist die Angabe von Anfangsund Randbedingungen dann sehr viel einfacher. Allerdings ist die Anwendbarkeit der primitiven Formulierung beschrankt. Zu ihrer Ableitung mu die Kontinuitatsgleichung uberall erfullt sein. Betrachtet man aber z. B. Hyperschallstromungen, so treten im Rechengebiet Diskontinuitaten auf. In der Stowelle - beim U bergang von U berschall- zu Unterschallgeschwindigkeiten - andern sich Temperatur und Dichte innerhalb der mittleren freien Weglange der Gasmischung nahezu sprunghaft um mehrere Groenordnungen. Eine primitive Formulierung lat hier keine numerische Losung mehr zu (siehe hierzu Kap. 4.4). Da in dieser Arbeit bzw. in vielen technischen Anwendungen aber primitive Erhaltungsgleichungen zur Simulation herangezogen werden, sollen sie hier ohne Ableitung aufgefuhrt werden:  Gesamtmasse (Kontinuitatsgleichung)

@ + div(~v) = 0 : @t

(2.11)

 Speziesmassen (Gleichungen fur die Massenbruche) @wi + 1 div~j  + ~v grad(w ) = 1 M ! : i i @t   i i

(2.12)

 Geschwindigkeitsfeld @~v + (~vgrad) ~v + 1 grad(p) + 1 div( ) = 0 : @t  

(2.13)

 Temperaturgleichung @T , 1 @p @t cp @t ns ns   X X + c1 Mi !ihi + c1 cp;i~ji grad(T ) + c1 div ~jq p i=1 p i=1 p 1 1 + ~v grad(T ) , c ~v grad(p) + c p : grad(~v ) = c1 qr : p

p

p

(2.14)

12

 KAPITEL 2. MATHEMATISCHE MODELLIERUNG REAKTIVER STROMUNGEN

In dieser Arbeit kommen die Verfahren der niten Dierenzen fur die raumliche Diskretisierung [7] zur Anwendung. Die zeitliche Integration erfolgt mit Hilfe des impliziten Extrapolationsverfahrens LIMEX [37, 38] im Fall der Flammenrechnungen. Zur Simulation des Wiedereintrittproblems wurde dazu ein Runge-Kutta-Verfahren verwandt, das die Quellterme durch eine semi-implizite Formulierung berechnet [39,40].

Kapitel 3

Automatische Reduktion von Reaktionsmechanismen 3.1 Einfuhrung In Kap. 2 wurden die Erhaltungsgleichungen vorgestellt, die man numerisch losen mu, wenn man reaktive Stromungen simulieren will. Die mathematische Modellierung ist dabei auf einem Niveau angelangt, das es erlaubt, laminare Verbrennungsprozesse mit detaillierten Reaktionsmechanismen und Transportmodellen zu berechnen [1{9]. Die U bereinstimmung zwischen experimentellen Messungen und den Ergebnissen mathematischer Modellierung ist dabei sehr gut. Allerdings ist man bei der Modellierung immer noch auf stark vereinfachte Geometrien festgelegt, da sich der Rechenzeitaufwand fur realistische Anwendungen nicht vertreten lat. Um die detaillierte Kinetik zu berechnen, mussen fur alle im Gemisch auftretenden Spezies Erhaltungsgleichungen gelost werden. Dies bedeutet zum einen eine sehr groe Anzahl von Gleichungen, die zu modellieren sind. Fur realistische Brennstoe wie z. B. Dodekan sind das mehr als 1000. Ein anderes Problem stellen die unterschiedlichen Zeit- und Langenskalen realer Verbrennungssysteme dar. Realistische Abmessungen technischer Verbrennungsanlagen liegen in der Groenordnung von Millimeter bis Meter. Auf der anderen Seite sind sehr kleine Strukturen wie Flammenfronten im m- bis mm-Bereich aufzulosen. Dies setzt z.B. adaptive Gitter vorraus. Eine weitere (oft viel groere) Schwierigkeit sind die groen Unterschiede in den Zeitskalen. Die Gleichungen der verschiedenen Spezies sind stark an die chemischen Quellterme gekoppelt, welche unterschiedlich schnell fur verschiedene Reaktionen sind. Die NO-Bildung dauert u.U. uber 1 Sekunde, wahrend andere chemische Reaktionen oft nur Nanosekunden benotigen. Dies bedingt eine starke Steifheit des Dierentialgleichungssystems chemisch reaktiver Stromungen. Die numerische Berechnung setzt deshalb implizite Losungsverfahren voraus. Die in dieser Arbeit verwendeten Methoden zur Losung der Dierentialgleichungen (Newton-Verfahren, semi-implizites Extrapolationsverfahren) verlangen eine Bestimmung der Jacobimatrizen, was aber einer der Rechenzeit-intensivsten Schritte der numerischen Simulation [12] ist (wachst mit der Anzahl der Spezies zum Quadrat). Es ist einsichtig, da fur realistische Brennstoe und Geometrien typischer Verbrennungssysteme eine vernunftige Simulation bzw. Parameterstudien nicht mehr moglich sind. Seit einigen Jahren bildet deshalb die "Ruckentwicklung\ von detaillierten Reaktionsmechanismen zu vereinfachten Mechanismen mit nur wenigen Spezies und Bruttoreaktionen einen Schwerpunkt in der Verbrennungsforschung [14,18{25,27,28]. Verschiedene Ansatze werden in Kap. 3.3.1 kurz vorgestellt. Zuvor wird aber in Kap. 3.2 etwas detaillierter auf die allgemeine chemische Dynamik eingegangen (insbesondere auf die Dynamik der Methanoxidation). In Kap. 3.3 wird der mathematische Ansatz der Methode der intrinsischen niedrig-dimensionalen Mannigfaltigkeiten der Reduktionsmethode, die in dieser Arbeit zur Anwendung kommt - ausfuhrlich diskutiert. 13

14

KAPITEL 3. AUTOMATISCHE REDUKTION VON REAKTIONSMECHANISMEN

3.2 Dynamik chemischer Reaktionssysteme 3.2.1 Allgemeine chemisch reagierende Systeme

Die mathematische Modellierung reaktiver Stromungen erfolgt durch das Aufstellen der Erhaltungsgleichungen, die die zeitliche A nderung eines Zustands (charakterisiert durch die Dichte , den Teilchendichten i, der spezi schen Enthalpie h und dem Geschwindigkeitsfeld ~v) beschreiben. Das Dierentialgleichungssystem wird durch algebraische Gleichungen, die die Temperatur T , die Dichte  und den Druck p mit den Zustandsvariablen verknupfen, geschlossen (siehe Kap. 2.2 fur die Erhaltungsgleichungen). Um eine klare Trennung zwischen chemischer Kinetik und physikalischen Prozessen (Stromung, molekularer Transport, etc.), raumlichen Gradienten und dem Geschwindigkeitsfeld zu bekommen, sollen fur die folgenden Ableitungen zunachst homogene Systeme betrachtet werden.

3.2.2 Homogene Reaktionssysteme Das mathematische Modell zur Mechanismenreduktion wird anhand einfacher homogener chemischer Reaktionssysteme abgeleitet. Es sei an dieser Stelle bemerkt, da fur die spateren Anwendungen keine derartige Einschrankungen notwendig sind. Die Methode lat sich auf ortsabhangige chemisch reaktive Stromung anwenden, wie in Kap. 4 noch gezeigt wird. Die Idee der Methode lat sich aber anschaulicher darlegen, wenn das Stromungsfeld von der "puren\ chemischen Dynamik entkoppelt ist. Um eine anschauliche Betrachtung zu erlauben, werden im folgenden zusatzlichen Annahmen gemacht: Das chemische System sei adiabat, isobar, homogen und geschlossen. Das Differentialgleichungssystem fur den Druck (Dichte), die spezi sche Enthalpie und die Speziesmassen vereinfacht sich in diesem Fall zu

@h = 0 @t @p = 0 (3.1) @t @wi = Mi!i i = 1; : : :; ns :  @t Es ist naturlich genauso erlaubt, isotherme (@T=@t = 0) anstatt adiabatischer und isochore (@=@t = 0) anstatt isobarer Systeme zu betrachten.

3.2.3 Begrisbildung

Zustands-, Mischungs-, Reaktionsraum Der Zustand eines chemisch reaktiven Gemisches ist durch seine Geschwindigkeit, die Dichte, die Enthalpie und die lokale chemische Zusammensetzung bestimmt. Da uns hier zunachst nur die chemische Dynamik interessiert, genugt es an dieser Stelle, das Geschwindigkeitsfeld auer Acht zu lassen, da die Anbindung aller Prozesse an die chemischen Reaktionen uber lokale A nderungen von Dichte, Druck und spezi scher Enthalpie erfolgt. Diese Groen, zusammen mit der chemischen Zusammensetzung (z.B. Massen- oder Molenbruche, spezi sche Molzahlen i = wi =Mi ), werden in dieser Arbeit als Zustandsvariablen bezeichnet (Achtung: Unterschiedliche De nition in der Thermodynamik). Das chemische System wird also uber die (2 + ns ) Zustandsvariablen charakterisiert, die sich im Zustandsvektor ~ = (h; p; 1; 2; : : :; ns )T zusammenfassen lassen. Ob man als Zustandsvariable die spezi sche Enthalpie oder die Temperatur bzw. die Dichte anstatt des Druckes wahlt, hangt vom betrachteten System ab (ob adiabat, isotherm, isochor oder isobar). Wahlt man sich das System der Variablen geeignet, reduziert sich der Raum, in dem die chemischen Reaktionen statt nden, u.U. um zwei Dimensionen (z.B. @h=@t = 0 und @p=@t = 0 fur adiabatische

3.2. DYNAMIK CHEMISCHER REAKTIONSSYSTEME

15

und isobare Systeme). D.h., das System ist allein durch die Dynamik der Spezieszusammensetzung ~ = (1 ; 2; : : :; ns ) bestimmt. Die Dynamik beschrankt sich dann auf den ns -dimensionalen Mischungsraum. Ist das chemische System geschlossen, nden darin nur chemische Reaktionen statt. Nach dem Gesetz der Elementerhaltung sind dabei deren molare Massen konstant. Damit beschrankt sich die Dynamik auf den (nr = ns , ne )-dimensionalen Reaktionsraum (ne = Anzahl der Elemente des Systems).

Element- und Reaktionsvektoren Element- und Reaktionvektoren besitzen eine besondere Bedeutung im Mischungs- und Reaktionsraum. Deshalb sollen sie hier kurz vorgestellt werden. ij sind die Anzahl der Atome des Elements i in der Spezies j . Der Elementvektor des Elements i ist dann ~i = (i1; i2; : : :; ins )T . Damit lassen sich die spezi schen Elementmolzahlen direkt aus der chemischen Zusammensetzung berechnen,

ei

=

ns X j =1

ij j = ~Ti ~ :

In geschlossenen homogenen Systemen sind die spezi schen Elementmolzahlen konstant (Elementerhaltung bei der chemischen Reaktion). Die Elementvektoren stehen also immer senkrecht auf den Trajektorien der chemischen Reaktionen im Mischungsraum. Chemische Reaktionen konnen durch Reaktionsvektoren ~i beschrieben werden: Die i-te Reaktion der Spezies Xj mit den stochiometrischen Koezienten aij ; a^ij ist durch

a1iX1 + a2i X2 + : : : + ans Xns ! a^1i X1 + : : : + a^nsi Xns gegeben. Der Reaktionsvektor der i-ten Reaktion wird durch die Dierenz der stochiometrischen Koezienten ji = aji , a^ji gebildet: ~i = (1i ; 2i ; : : :; nsi )T . Fur die Element- und Reaktionsvektoren gilt die Orthogonalitatsbedingung

~i ~Tj = 0

i = 1; : : :; ncr ; j = 1; : : :; ne :

Betrachtet man die Element- und Reaktionsvektoren im Zustandsraum, so ergibt sich fur sie

~^i = (0; 0; 1i; 2i; : : :; nsi )T und ~^j = (0; 0; 1j ; 2j ; : : :; nsj )T . Im folgenden wird allerdings keine

Unterscheidung zwischen den Element- und Reaktionsvektoren im Mischungs- bzw. Zustandsraum gemacht.

Erlaubte Zustande Die erlaubten Zustande im Mischungsraum sind durch die chemisch physikalischen Variablen bestimmt:

 

T ~ >0   p ~ >0    ~ >0  wi  1 Pns w = Pns0 M Pne e e i=1 i i=1 i i = j =1 Mj j = 1 :

16

KAPITEL 3. AUTOMATISCHE REDUKTION VON REAKTIONSMECHANISMEN

3.2.4 Chemische Kinetik am Beispiel der Methanoxidation In diesem Abschnitt soll eingehender die chemische Dynamik, insbesondere der Methanoxidation, erlautert werden und eine anschauliche Einfuhrung in die Idee der Methode der intrinsischen niedrig-dimensionalen Mannigfaltigkeiten gegeben werden. Wie schon angedeutet, wird das chemische Reaktionssystem als geschlossen, adiabat, isobar und homogen betrachtet und durch das Dierentialgleichungssystem (3.2) beschrieben. Die chemischen Reaktionen bewegen einen Zustand ~ = (h; p; ~)T des Systems entlang von Trajektorien im Zustandsraum, die durch ihren Anfangszustand charakterisiert sind. Da in diesem Fall das System homogen, isobar, adiabat und geschlossen ist, sind nur die chemischen Reaktionen fur die Bewegung im Reaktionsraum verantwortlich. Nach unendlich langer Zeit ist das System schlielich im Gleichgewichtspunkt angelangt, in dem alle Trajektorien enden. Fur obige Annahmen ist der Gleichgewichtspunkt eine eindeutige Funktion der spezi schen Enthalpie, des Druckes und der Elementzusammensetzung. Im Reaktionsraum ist der chemische Gleichgewichtswert ein Punkt, im Mischungsraum eine (ne , 1)-dimensionale und im Zustandsraum eine (ne + 1)-dimensionale Mannigfaltigkeit. Dieser Gleichgewichtswert wird fur viele praktische Anwendungsfalle der Simulation reaktiver (turbulenter) Stromungen als chemisches Modell angenommen. Da die meisten realistischen Verbrennungssysteme sehr komplexer Natur sind, ist es nahezu unmoglich solche Systeme vollstandig zu beschreiben, nicht zuletzt wegen der sehr steifen Chemie. Deshalb wird sehr oft die Annahme der Null-Schritt-Chemie "gemischt gleich verbrannt\ benutzt [41], d.h. die Annahme, da das chemische Gleichgewicht im Vergleich zu den physikalischen Zeitskalen sehr schnell erreicht wird. Naturlich ist mit derartigen Annahmen keine Aussage mehr uber die chemische Dynamik moglich, insbesonders wenn die Chemie der geschwindigkeitsbestimmende Faktor im reaktiven Gemisch ist. Das andere Extrem sind die detaillierten Reaktionsmechanismen. Sie beschreiben alle moglichen Zeitskalen der chemischen Dynamik, auch dann, wenn sie in keiner Weise geschwindigkeitsbestimmend sind. Bei Verwendung detaillierter Reaktionsmechanismen wird die zeitliche Entwicklung eines Zustands durch (ns + 2) Zustandsvariablen beschrieben. Die Annahme "gemischt gleich verbrannt\ dagegen beschreibt den Zustand nur durch die spezi sche Enthalpie, den Druck und der Elementzusammensetzung, wobei aber samtliche chemische Dynamik vernachlassigt wird. Wunschenswert sind aber Ansatze, die durch eine geringe Anzahl sogenannter Reaktionsfortschrittvariablen ~ die Dynamik gut beschreiben konnen, auf der anderen Seite aber die groe Dimension des detaillierten Dierentialgleichungssystems reduzieren. Der Zustand wird dann durch weniger Variablen beschrieben ~ red = ~ (h; p; ~; ~)T . An dieser Stelle soll nun anschaulich gezeigt werden, wie die chemische Dynamik im Reaktionsraum Antworten liefern kann, wie ein Reaktionsmechanismus auf geeignete Art und Weise reduziert werden kann. Das soll am konkreten Beispiel der Methanoxidation vorgestellt werden, bevor im nachsten Kapitel die mathematische Identi zierung der intrinsischen niedrig-dimensionalen Mannigfaltigkeiten beschrieben wird. Das chemische Reaktionssystem, das eingehender diskutiert wird, ist ein stochiometrisches Gemisch aus Methan und Luft. Methan ist der niedrigste Kohlenwassersto (CH4 ) und bietet erste Einblicke in die Verbrennung realistischer Brennstoe (wie z.B. Heptan oder Dodekan) [42]. Auerdem ist es mit mehr als 96% im Erdgas enthalten, so da man aus Ergebnissen von CH4 -O2 -N2 Simulationen auch auf Ergebnisse der Erdgasverbrennung schlieen kann. Der zugrunde liegende Reaktionsmechanismus beinhaltet 34 chemische Spezies (ns = 34) und die 4 Elemente (ne = 4) H, O, C, N. Die Spezies sind uber 288 Reaktionen miteinander gekoppelt [42] (es existieren auch schon Mechanismen der Methanoxidation, die z.B. die Stickoxidbildung beinhalten. Dann ist man selbst fur das Methan-Luft-System bei uber 100 Spezies angelangt [43]). Der Mischungsraum ist also 34-dimensional und der Reaktionsraum 30-dimensional. Als konkretes Beispiel wird eine stochiometrische Mischung aus Methan und Luft (9.5% CH4 , 19.0% O2 ,

3.2. DYNAMIK CHEMISCHER REAKTIONSSYSTEME

17

71.5% N2 ) bei Atmospharendruck untersucht, da dieses zum einen experimentell gut zuganglich ist und zum anderen ein durchaus realistisches Reaktionssystem darstellt. Die spezi sche Enthalpie des unverbrannten Gemisches bei einer Temperatur von T = 298 K betragt h = ,2:58 105 J/kg. Die Gleichgewichtsmischung ist fur dieselbe spezi sche Enthalpie ca. 2240 K hei. Fur dieses System wurden zunachst fur verschiedene Anfangsbedingungen einige Simulationen der Reaktionsgleichungen (3.2) mit Hilfe des detaillierten Reaktionsmechanismus durchgefuhrt. Das verwendete numerische Verfahren ist ein implizites Extrapolationsverfahren mit variabler Ordnungs- und Schrittweitensteuerung [37, 38]. Fur diese Berechnungen ergeben sich verschieden Trajektorien im Reaktionsraum fur die verschiedenen Anfangsbedingungen wie sie in Tab. 3.1 zusammengestellt sind. Alle Trajektorien enden fur sehr groe Zeiten im Gleichgewichtspunkt. Da man unmoglich den hoch-dimensionalen Reaktionsraum graphisch darstellen kann, werden verschiedene zweidimensionale Projektionen der Phasenraume gezeigt (siehe Abb. 3.1 und Abb. 3.2). An ihnen lat sich sehr anschaulich die Dynamik des chemischen Systems diskutieren. Tabelle 3.1: Anfangsbedingungen der Trajektorien, wi sind Massenbruche. Symbol    M + N 

  H

wCO

2

.106 .136 .151 .076 .076 .076 .113 .151 .121 .144

wH O wO 2

.118 .078 .108 .093 .108 .077 .077 .093 .118 .121

2

.022 .047 .014 .055 .041 .069 .055 .027 .016 .005

wCO

.029 .0096 .0 .048 .048 .048 .024 .0 .019 .005

wH

2

.0007 .0052 .0017 .0035 .0017 .0052 .0052 .0035 .0007 .0004

T (K) 2080 1831 2184 1716 1864 1567 1733 2040 2144 2267

Betrachtet man die Trajektorien im CO2 -CO-Phasenraum, so bilden sie naherungsweise eine Gerade mit der Steigung von ungefahr -1. Diese beiden Spezies beinhalten fast den gesamten Kohlenstoanteil. Die anderen kohlenstohaltigen Spezies treten nur in Spuren auf. Die Abnahme des CO-Anteils fuhrt also direkt zu einer Zunahme von CO2 . Die Trajektorien in der CO2 -H2 -Projektion verlaufen dagegen im Anfangsstadium fast vertikal. Das bedeutet, da zunachst der Wassersto abreagiert, bevor dann, wie im CO2 -CO-Diagramm zu sehen, CO zu CO2 aufoxidiert wird. Dasselbe Bild ist auch fur die CO2 -O2 -Projektion zu sehen, wobei hier die Trajektorien im Anfangsbereich nicht ganz so steil verlaufen wie beim Wassersto. Das deutet darauf hin, da zunachst der Wassersto verbraucht wird und danach die O2 -Molekule (ohne da CO2 gebildet wird). Erst danach wird das CO zu CO2 aufoxidiert. Dieser Sachverhalt zeigt sich sehr schon im Phasendiagramm der Molekule H2 und O2 . Das CO2 -H2 O -Phasenraumdiagramm zeigt ebenso deutlich, da zuerst das Wasser gebildet wird, bevor CO zu CO2 oxidiert (steiler Anstieg der Trajektorien im Anfangsbereich). Dies wird durch das H2 O -CO-Phasendiagramm unterstrichen. Erst nach anfanglicher H2O -Bildung erfolgt ein merklicher Ruckgang der CO-Konzentration in der Nahe des Gleichgewichtspunktes. Die Radikalbildung sieht man sehr schon, wenn man diese uber eine stabile Spezies auftragt. In Abb. 3.2 sind die Phasendiagramme CO2 -O und CO2 -H zu sehen. Die Trajektorien beginnen immer auf der CO2 -Achse, da als Anfangsbedingung immer eine Mischung der stabilen Spezies CO2 , H2 O , CO, O2 , N2 und H2 angenommen wurde. In allen Phasenraumdiagrammen Radikal uber stabile Spezies ist ein sehr schneller Anstieg der Radikale zu sehen, wahrend sich CO2 oder H2 O nicht andern. Dieser Anstieg charakterisiert den Reaktionsbeginn.

KAPITEL 3. AUTOMATISCHE REDUKTION VON REAKTIONSMECHANISMEN

18

wC O

wH 2

0.05

0.006

0.04 0.004

0.03

0.02 0.002

0.01

0 0.06

0.08

0.1

0.12

0.14

0.16

0 0.06

0.08

0.1

wC O 2

0.12

0.14

0.16

wC O 2

wO 2

wH 2 0.006

0.07

0.004 0.035 0.002

0 0.06

0 0.08

0.1

0.12

0.14

0

0.16

0.035

wC O 2

0.07

wO 2

wC O

wH 2 O

0.05

0.13 0.12

0.04

0.11 0.03 0.1 0.02 0.09 0.01

0.08 0.07 0.06

0.08

0.1

0.12 wC O 2

0.14

0.16

0 0.07

0.08

0.09

0.1

0.11

0.12

0.13

wH 2 O

Abbildung 3.1: Trajektorien im Reaktionsraum fur das stochiometrische CH4 -O2 -N2 -System fur verschiedene Anfangsbedingungen.  kennzeichnet den chemischen Gleichgewichtspunkt.

3.2. DYNAMIK CHEMISCHER REAKTIONSSYSTEME

wO

19

wH

0.005

0.0009

0.004 0.0006

0.003

0.002 0.0003

0.001

0 0.06

0.08

0.1

0.12

0.14

0 0.06

0.16

0.08

0.1

wC O 2

0.12

0.14

0.16

0.003

0.004

0.005

wC O 2

wC H O

wO H

7 10-7

0.006

6 10-7 5 10-7 0.004 4 10-7 3 10-7 0.002

2 10-7 1 10-7

0

0 0

0.001

0.002

0.003

0.004

0

0.005

0.001

0.002

wO

wO

T(K)

T(K)

2200

2200

2000

2000

1800

1800

1600

1600

0.06

0.08

0.1

0.12 wC O 2

0.14

0.16

0.07

0.08

0.09

0.1

0.11

0.12

0.13

wH 2 O

Abbildung 3.2: Trajektorien im Reaktionsraum fur das stochiometrische CH4 -O2 -N2 -System fur verschiedene Anfangsbedingungen.  kennzeichnet den chemischen Gleichgewichtspunkt.

20

KAPITEL 3. AUTOMATISCHE REDUKTION VON REAKTIONSMECHANISMEN

Einen genaueren U berblick uber die chemische Kinetik kann man sich verschaen, wenn man im Phasenraum Radikalkonzentrationen gegeneinander auftragt. Das Ergebnis ist meist ein kompliziertes Diagramm, wobei naturlich bedingt durch die Anfangsbedingungen (keine Radikale) die Trajektorien im Ursprung beginnen. Zuletzt sei noch die Temperatur in Abhangigkeit der stabilen Spezies diskutiert. Ein merklicher Temperaturanstieg ist erst mit der einsetzenden CO2 -Produktion zu sehen. Dieser Sachverhalt ist auf die stark exotherme Bruttoreaktion 2CO + O2 ! 2CO2 zuruckzufuhren. Die Temperatur ist in der Hauptsache eine Funktion der CO2 -Bildung. Weiter oben wurde gesagt, da in praktischen Anwendungen sehr oft die Annahme "gemischt gleich verbrannt\ verwendet wird. Das bedeutet nichts anderes, als da die gesamte chemische Kinetik vernachlassigt wird und der chemische Umsatz nur durch die chemische Gleichgewichtslage beschrieben wird. Der Gleichgewichtspunkt ist in Abb. 3.1 und Abb. 3.2 durch den groen Kreis gekennzeichnet. Das System reagiert bei Verwendung dieser Vereinfachung unendlich schnell, bzw. es wird vorrausgesetzt, da alle physikalischen Prozesse (z.B. Diusion, Warmeleitung oder Turbulenz) sehr viel groere Zeitskalen besitzen als alle chemische Reaktionen. Ist das der Fall, ist die Annahme chemischen Gleichgewichts durchaus erlaubt. Das wird mit Abb. 3.3 und Abb. 3.4 veranschaulicht. Hier sind dieselben Phasenraumdiagramme dargestellt wie in Abb. 3.1 und Abb. 3.2, nur sind die ersten 5 ms der zeitlichen Entwicklung der Trajektorien abgeschnitten worden. Was von der komplexen Struktur ubrigbleibt entspricht in etwa dem chemischen Gleichgewichtspunkt. Besitzen also molekularer Transport und Warmeleitung immer groere Zeitskalen als mehrere Millisekunden, kann man im Falle der Methanoxidation annehmen, da die chemische Kinetik schon ins Gleichgewicht relaxiert ist, bevor physikalischer Transport einsetzt. Dann kann die detaillierte Chemie durch die "Null-dimensionale Dynamik\ gemischt gleich verbrannt ersetzt werden. In typischen Verbrennungssystemen beanspruchen die physikalischen Prozesse einen etwas anderen Bereich an Zeitskalen (vergleiche auch Abb. 3.9). Er liegt in etwa zwischen mehreren Millisekunden und 10,4 bis 10,5 Sekunden. D.h., die obige drastische Vereinfachung der Gleichgewichtschemie ist nicht erlaubt, da die Einstellung des chemischen Gleichgewichts ebenfalls mehrere Millisekunden benotigt. Die chemische Kinetik koppelt somit direkt an Stromung und molekularen Transport. Um diese Kopplung zu beschreiben mu man wissen, an welchem Punkt im Reaktionsraum man sich be ndet, bzw. wie die genaue chemische Zusammensetzung und der chemische Umsatz aussehen. D.h., man mu die komplette chemische Dynamik beschreiben. Sind die Zeitskalen der physikalischen Prozesse groer als die Zeitskalen der chemischen Gleichgewichtseinstellung, so ist das Modell "gemischt gleich verbrannt\ ein durchaus brauchbares chemisches Modell. D.h., auf reale physikalische Zeitskalen ubertragen (mehrere Millisekunden), genugt es chemische Reaktionen die langere Zeit benotigen als z.B. 50s zu berucksichtigen (unter der Vorraussetzung, da man nicht explizit an der sehr schnellen Chemie interessiert ist). Man betrachtet also nicht die komplette chemische Dynamik, sondern nur Prozesse mit Zeitskalen groer als 50s bis zur Einstellung des chemischen Gleichgewichts. Alle schnelleren chemischen Prozesse sind zu diesem Zeitpunkt schon ins lokale Gleichgewicht hineinrelaxiert. Das ist wiederum fur dieselben Phasenraumportraits dargestellt (vergleiche Abb. 3.1 und Abb. 3.2). In Abb. 3.5 und Abb. 3.6 sind im Gegensatz zu den vollstandigen Trajektorien die ersten 50s der chemischen Kinetik abgeschnitten worden. Man erhalt wiederum ein stark vereinfachtes Bild der chemischen Dynamik. Die moglichen chemischen Zustande, die sich nach 50s noch einstellen konnen, liegen alle auf einer Linie im Reaktionsraum. Der Zustand des chemischen Systems lat sich durch eine eindimensionale Mannigfaltigkeit beschreiben, wenn man keine schnelleren Prozesse als 50s au osen will. Sind aber doch noch schnellere Zeitskalen interessant, so wird das System durch eine Flache (entsprechend eine zweidimensionalen Mannigfaltigkeit), danach durch eine dreidimensionale Hyper ache (entsprechend einer dreidimensionalen Mannigfaltigkeit) im Reaktionsraum beschrieben. Fuhrt man das zu immer kleiner werdenden Zeitskalen weiter, ist man irgendwann bei der detaillierten Beschreibung angelangt, die ja alle moglichen Zeitskalen beinhaltet. Fur gegebene spezi sche

3.2. DYNAMIK CHEMISCHER REAKTIONSSYSTEME

wC O

21

wH 2

0.05

0.006

0.04 0.004

0.03

0.02 0.002

0.01

0 0.06

0.08

0.1

0.12

0.14

0.16

0 0.06

0.08

0.1

wC O 2

0.12

0.14

0.16

wC O 2 wH 2

wO 2

0.006

0.07

0.004 0.035 0.002

0 0.06

0 0.08

0.1

0.12

0.14

0

0.16

0.035

wC O 2

0.07

wO 2

wH 2 O

wC O

0.13

0.05

0.12

0.04

0.11 0.03 0.1 0.02 0.09 0.01

0.08 0.07 0.06

0.08

0.1

0.12 wC O 2

0.14

0.16

0 0.07

0.08

0.09

0.1

0.11

0.12

0.13

wH 2 O

Abbildung 3.3: Trajektorien im Reaktionsraum fur das stochiometrische CH4 -O2 -N2 -System fur verschiedene Anfangsbedingungen, wobei hier die ersten 5 ms der chemischen Dynamik abgeschnitten wurden.  kennzeichnet den chemischen Gleichgewichtspunkt.

KAPITEL 3. AUTOMATISCHE REDUKTION VON REAKTIONSMECHANISMEN

22

wO

wO

0.06

0.06

0.04

0.04

0.02

0.02

0 0.06

0.08

0.1

0.12

0.14

0.16

0 0.06

0.08

0.1

wC O 2

0.12

0.14

0.16

0.003

0.004

0.005

wC O 2 wC H O

wO H

7 10-7

0.006

6 10-7 5 10-7 0.004 4 10-7 3 10-7 0.002

2 10-7 1 10-7

0

0 0

0.001

0.002

0.003

0.004

0

0.005

0.001

0.002

wO

wO

T(K)

T(K)

2200

2200

2000

2000

1800

1800

1600

1600

0.06

0.08

0.1

0.12 wC O 2

0.14

0.16

0.07

0.08

0.09

0.1

0.11

0.12

0.13

wH 2 O

Abbildung 3.4: Trajektorien im Reaktionsraum fur das stochiometrische CH4 -O2 -N2 -System fur verschiedene Anfangsbedingungen, wobei hier die ersten 5 ms der chemischen Dynamik abgeschnitten wurden.  kennzeichnet den chemischen Gleichgewichtspunkt.

3.2. DYNAMIK CHEMISCHER REAKTIONSSYSTEME

wC O

23

wH 2

0.05

0.006

0.04 0.004

0.03

0.02 0.002

0.01

0 0.06

0.08

0.1

0.12

0.14

0.16

0 0.06

0.08

0.1

wC O 2

0.12

0.14

0.16

wC O 2 wH 2

wO 2

0.006

0.07

0.004 0.035 0.002

0 0.06

0 0.08

0.1

0.12

0.14

0

0.16

0.035

wC O 2

0.07

wO 2

wC O

wH 2 O

0.05

0.13 0.12

0.04

0.11 0.03 0.1 0.02 0.09 0.01

0.08 0.07 0.06

0.08

0.1

0.12 wC O 2

0.14

0.16

0 0.07

0.08

0.09

0.1

0.11

0.12

0.13

wH 2 O

Abbildung 3.5: Trajektorien im Reaktionsraum fur das stochiometrische CH4 -O2 -N2 -System fur verschiedene Anfangsbedingungen, wobei hier die ersten 50s der chemischen Dynamik abgeschnitten wurden.  kennzeichnet den chemischen Gleichgewichtspunkt.

24

KAPITEL 3. AUTOMATISCHE REDUKTION VON REAKTIONSMECHANISMEN

wO

wH

0.005

0.0009

0.004 0.0006

0.003

0.002 0.0003

0.001

0 0.06

0.08

0.1

0.12

0.14

0 0.06

0.16

0.08

0.1

wC O 2

0.12

0.14

0.16

0.003

0.004

0.005

wC O 2

wC H O

wO H

7 10-7

0.006

6 10-7 5 10-7 0.004 4 10-7 3 10-7 0.002

2 10-7 1 10-7

0

0 0

0.001

0.002

0.003

0.004

0

0.005

0.001

0.002

wO

wO

T(K)

T(K)

2200

2200

2000

2000

1800

1800

1600

1600

0.06

0.08

0.1

0.12 wC O 2

0.14

0.16

0.07

0.08

0.09

0.1

0.11

0.12

0.13

wH 2 O

Abbildung 3.6: Trajektorien im Reaktionsraum fur das stochiometrische CH4 -O2 -N2 -System fur verschiedene Anfangsbedingungen, wobei hier die ersten 50s der chemischen Dynamik abgeschnitten wurden.  kennzeichnet den chemischen Gleichgewichtspunkt.

3.2. DYNAMIK CHEMISCHER REAKTIONSSYSTEME

25

Enthalpie, Druck und Elementzusammensetzung entspricht das einer Beschreibung durch eine nr dimensionale Mannigfaltigkeit. In anderen Worten ausgedruckt: Je nachdem, wie viele Zeitskalen berucksichtigt werden sollen (je kleiner desto detaillierter die Beschreibung), wird die chemische Dynamik durch eine Bewegung auf einer m-dimensionalen Mannigfaltigkeit im Reaktionsraum angenahert. Die Dynamik im viel-dimensionalen Reaktionsraum kann also durch eine Dynamik auf niedrig-dimensionalen Mannigfaltigkeiten ersetzt werden. Analog zum System gemischt gleich verbrannt (der Gleichgewichtspunkt entspricht einer null-dimensionalen Mannigfaltigkeit im Reaktionsraum) konnen die Mannigfaltigkeiten identi ziert werden und zur Stromungssimulation benutzt werden. Die chemische Dynamik lat sich also, ist man nur am Zeitskalenbereich von 50s und groer interessiert, vollstandig durch die spezi sche Enthalpie, den Druck, der Elementzusammensetzung und einer weiteren rektiven Variable beschreiben, die die Bewegung auf der Linie in Abb. 3.5 und Abb. 3.6 beschreibt. Ist diese Variable bekannt, kennt man den genauen Zustand des Systems im Reaktionsraum. Um die Dynamik des reaktiven Gemisches zu beschreiben, mu man nur die zeitliche A nderung der einen reaktiven Variablen kennen. Auf diese Art hat man die chemische Kinetik in zwei Systeme unterteilt. Zum einen, in eine schnelle Relaxationsphase von jedem beliebigen Anfangszustand ausgehend bis zum Erreichen der Trajektorien solcher niedrig-dimensionaler Unterraume. Das zweite System ist bestimmt durch die anschlieende langsame Bewegung in den niedrig-dimensionalen Unterraumen (hier Zeitskalen > 50s), bis schlielich das chemische Gleichgewicht erreicht ist. Ist man in der Lage, diese zwei Phasen der chemischen Dynamik zu trennen, kann die Kinetik vieler realer Verbrennungssysteme durch wenige reaktive Fortschrittsvariablen beschrieben werden, wenn die chemischen Einschwingvorgange nicht interessieren. Der Grund hierfur sind die sehr vielen Reaktionen, die sehr schnell ablaufen und deshalb sehr schnell lokale Gleichgewichtsprozesse im Reaktionsraum erreicht werden. Die Folge ist, da Nettoreaktionsgeschwindigkeiten aus Hin- und Ruckreaktion verschwinden. Daraus folgen partielle Gleichgewichte verschiedener Reaktionen. Auch quasistationare Zustande fur Spezieskonzentrationen konnen daraus entstehen. Man spricht davon, wenn chemische Spezies durch langsame Reaktionen gebildet werden, danach aber sehr schnell wieder verbraucht werden. Fur beide Falle erhalt man algebraische Gleichungen fur einzele Spezies, die dann nicht mehr durch Integration von Dierentialgleichungen aufgelost werden mussen. Fur den gerade angesprochenen Fall der Methanoxidation und interessanten Zeitskalen > 50s (vergleiche Abb. 3.5 und Abb. 3.6) sind dann genau 29 partielle Gleichgewichte und stationare Zustande erreicht. Damit ist das System durch h; p; H; O ; C ; N und einer Reaktionsfortschrittsvariablen beschreibbar, und damit alle anderen im chemischen Reaktionssystem vorkommenden Spezies bestimmt.

3.2.5 Erhaltungsgroen, partielle Gleichgewichte und stationare Zustande

Die Dynamik allgemeiner Reaktionssysteme wird durch die Reaktionsgleichungen fur den (ns + 2)dimensionalen Zustandsvektor beschrieben,

@ ~ = F~  ~  : @t

~ = (h; p; 1; : : :; ns ) beschreibt an jedem Ort und Zeit den Zustand der chemischen Kinetik.

Weiter oben wurde schon angedeutet, da im Zustandsraum bestimmte Vektoren eine Sonderrolle einnehmen. Das sind die Element- und Reaktionsvektoren. Bewegt man sich entlang von Elementvektoren ~i , andert sich die Elementzusammensetzung, was bei den hier betrachteten homogenen Systemen nicht erlaubt ist. Entlang der Reaktionsvektoren ~i andert sich nur die Spezieszusammensetzung durch chemische Reaktionen, nicht aber die Elementzusammensetzung. Die zugrundeliegenden Gleichungen sollen wieder fur adiabatische, isobare und homogene Systeme gelten. D.h., die spezi sche Enthalpie, der Druck und die Elementzusammensetzung sind konstant. Chemische

26

KAPITEL 3. AUTOMATISCHE REDUKTION VON REAKTIONSMECHANISMEN wH 2 O 0.13 schnell

0.12

langsam

schnell 0.11 langsam schnell

0.1 0.09

schnell 0.08

schnell schnell

0.07 0.06

0.08

0.1

0.12

0.14

0.16

wC O 2

Abbildung 3.7: Relaxation von Trajektorien auf eine eindimensionale Mannigfaltigkeit. Reaktionen nden nur im (ns , ne )-dimensionalen Reaktionsraum statt, in dem die Spezieszusammensetzungen, die durch chemische Reaktionen ineinander uberfuhrt werden konnen, beschrieben werden. Der Reaktionsraum wird durch die Reaktionsvektoren aufgespannt. An dieser Stelle wird gezeigt, da sich der Mischungsraum (siehe Kap. 3.2.3) nicht nur durch die Spezieszusammensetzung (z.B. spezi sche Molzahlen, i = wi =Mi ) beschreiben lat. Dazu wird die Basis (die im allgemeinen durch die Eigenvektoren der Jacobimatrix gebildet wird) in eine alternative Basis uberfuhrt, die sich aus ns linear unabhangigen Vektoren aufspannen lat [24]. Die Dynamik des Mischungsraumes, beschrieben durch

@ ~ = F~ ~  = R~r ; @t

(3.2)

wird in eine alternative Basis durch die Transformationsgleichung

B ~ = ~

(3.3)

mit den neuen Variablen ~ uberfuhrt. B ist die invertierbare Transformationsmatrix, R die Matrix der stochiometrischen Koezienten der chemischen Reaktionen und ~r der Vektor der Reaktionsgeschwindigkeiten der Spezies. R ist im Mischungsraum eine (ns  ncr )-dimensionale Matrix mit den Reaktionsvektoren ~^1 ; : : :; ~^ncr als Spaltenvektoren (ncr = Anzahl der Reaktionen),

0 j ^ R=B @ ~ 1 j

1 j j C ~^2 : : : ~^ncr A : j j

Die Erhaltungsgleichungen (3.2) transformieren sich zu



~_ = B ,1 F~ ~ = B,1 R~r : Die Transformationsmatrix sei speziell durch die ne Elementvektoren ~ i und nr = ns , ne Reaktionsvektoren ~i , die den Reaktionsraum aufspannen (jeweils als Spaltenvektoren geschrieben)

3.2. DYNAMIK CHEMISCHER REAKTIONSSYSTEME

27

aufgespannt,

0, BB BB , 1 B =B BB ,, BB @, ,

0 1 j j j j j j B=B @ ~1 ~2 : : : ~ne ~1 ~2 : : : ~nr CA j j j j j j

~ L1 , 1

CC C ~Lne , C C: CC ~1L , C .. C . ,A ~nLr , .. .

B ,1 ist die Inverse von B, mit den Vektoren ~Li und ~iL als Zeilenvektoren, die sich durch die Invertierung von B ergeben. Element- und Reaktionsvektoren stehen senkrecht aufeinander (die

Elementzusammensetzung andert sich nicht durch Reaktion, nur die Spezieszusammensetzung), ~Ti ~j = 0 i = 1; : : :; ne j = 1; : : :; ncr ; oder allgemeiner ~ Ti R = 0 i = 1; : : :; ne : (3.4) Fuhrt man nun die Transformation (3.3) aus, so erhalt man

0, BB BB ~_ = B,1 R~r = B BB ,, BB @, ,

~ L1 , 1

CC 0 C j ~Lne , C C B@ ~^1 CC j ~1L , C .. C . ,A ~nLr , .. .

0 1 BB j CB B : : : ~^n A B B j B B@

j ^~ 2 j

cr

r1 1

CC .. C . C CC : CC A

rn

cr

Mit den obigen Orthogonalitatsbedingungen (3.4) und fi = ~^i R~r erhalt man fur das Dierentialgleichungssystem 0  1 0 0 1 BB ..1 CC BB .. CC B . CC BB . CC @ ~ = B B ne CC = BB 0 CC : (3.5) BB ne+1 CC BB f1 CC @t B B@ ... CA B@ ... CA

ne+nr

fnr

Die Komponenten von ~ entsprechen den Anteilen in den Richtungen der Element- und Reaktionsvektoren. Dabei sind die ersten ne Komponenten i Koordinaten der Elementzusammensetzung und die restlichen nr Komponenten i Koordinaten des Reaktionsraumes,

0  BB ..1 BB . e ~ = B BB nn+1 BB e. @ .. ns

1 0 CC BB CC BB CC = BB CC BB CA B@

1 1 .. C . C C

CC CC : .. C C . A

ne 1

nr

Das Gleichungssystem (3.5) wird dann zu ~_ = 0  oder ~ = konstant ~_ = f~ h; p; ~; ~ :

(3.6) (3.7)

28

KAPITEL 3. AUTOMATISCHE REDUKTION VON REAKTIONSMECHANISMEN

In der Basis der Element- und Reaktionsvektoren ndet also eine Aufspaltung der Variablen statt. Durch die ersten ne algebraischen Gleichungen werden die Reaktionen auf den Reaktionsraum (ein Unterraum des Mischungsraumes) beschrankt. Innerhalb des Reaktionsraumes wird die Dynamik durch die nr Dierentialgleichungen (3.7) beschrieben. Fur das hier vorliegende Reaktionssystem ergeben sich automatisch die Erhaltungsgleichungen fur die Elemente. Zusatzlich konnen partielle Gleichgewichte fur Reaktionen gewahlt werden. D.h., man nimmt an, da Hin- und Ruckreaktionen gleich schnell sind. Fur Konzentrationen bestimmter Spezies ergeben sich dann aus dem Massenwirkungsgesetz der Thermodynamik algebraische Gleichungen als Funktion anderer Spezieskonzentrationen. Jede Annahme eines partiellen Gleichgewichtes reduziert somit die Anzahl der Dierentialgleichungen, die zu losen sind, um Eins. Physikalisch wird ein System, das aus dem Gleichgewicht gestort wird, sehr schnell wieder in den Zustand partieller Gleichgewichte zuruckrelaxieren. Im Fall der vorgestellten alternativen Darstellung in der Basis der Element- und Reaktionsvektoren konnen die partiellen Gleichgewichte sehr einfach spezi ziert werden [24,25, 44{47]. Sie sind einfach dadurch gegeben, da die Komponenten der Reaktionsgeschwindigkeiten in Richtung der betreenden Reaktionsvektoren verschwinden. Quasistationare Spezies sind Zwischenprodukte, die sofort, nachdem sie gebildet wurden, wieder aufgebraucht werden. Die zeitliche A nderung dieser Spezies ist dann gleich Null (Fi = 0). Auch mit dieser Annahme kann die Dimension des betreenden Dierentialgleichungssystems um Eins reduziert werden. Fur praktische Anwendungen liegt meistens eine Kombination aus np partiellen Gleichgewichten und nq Quasistationaritatsannahmen fur Spezies vor. Die chemische Kinetik beschrankt sich also auf einen (m = ns , ne , np , nq )-dimensionalen Unterraum des Reaktionsraumes. Liegt z.B. ein Ein-Schritt-Mechanismus vor, so wird die Dynamik durch die Bewegung auf einer Linie im Reaktionsraum beschrieben. Die Gleichgewichts-Chemie entspricht einem Punkt, also einem Null-SchrittMechanismus, im Reaktionsraum. Anhand der chemischen Dynamik, wie sie hier vorgestellt wurde, lassen sich qualitativ solche Vorhersagen treen und Reaktionsmechanismen vereinfachen. Dazu ist allerdings eine sehr gute Kenntnis der verwendeten Mechanismen und der physikalisch-chemischen Vorgange notig. Auerdem sind die Mechanismen nur sehr global reduziert. Einer A nderung von Reaktionspfaden (z.B. durch Temperaturanderung oder beim U bergang von einem fetten ins magere Verbrennungsregime) wird nur unter erheblichem Aufwand Rechnung getragen.

3.2.6 Eigenwertanalysen chemischer Systeme Wie ein chemisches System auf eine Veranderung bzw. eine Storung reagiert, lat sich mit Hilfe der Eigenwertanalyse vorhersagen. Dabei nimmt man an einem bestimmten Punkt im Zustandsraum eine kleine Storung vor und betrachtet das Verhalten in der Umgebung des Punktes. Das System der partiellen Dierentialgleichungen lat sich in seiner einfachsten Form in vektorieller Darstellung schreiben,

@ ~ = F~  ~  @t

~ (t = 0) = ~0 :

(3.8)

Es soll die Kinetik des chemischen Systems um den Punkt ~0 untersucht werden. Dazu wird an diesem Punkt eine Storung ~ aufgebracht und die Reaktion auf das Dierentialgleichungssystem betrachtet. Ferner wird angenommen, da sich das Gleichungssystem um den Punkt ~0 linearisieren lat:

@ ~ = F~  ~  + F  ~ , ~  + O  ~ , ~ 2 + : : : 0 0 0 ~ @t 0

~ (t = 0) = ~0 :

(3.9)

3.2. DYNAMIK CHEMISCHER REAKTIONSSYSTEME

29

In (3.9) ist F die Jacobimatrix des Dierentialgleichungssystem (3.8) am Punkt ~ 0, die folgendermaen de niert ist 0 @F @F @F 1 2 1

@ 2 @F2 @ 2









@ n @F2 @ n

@Fn @ 1

@Fn @ 2






@Fn @ n

1 1

1

BB @@F F =B BB @ .. @ .

.. .

...

1

CC C: .. C . C A

Der Ausdruck O(( ~ , ~0)2 ) bezeichnet Terme zweiter und hoherer Ordnung. Er wird bei der Linearisierung vernachlassigt, was aufgrund der Kleinheit der aufgepragten Storung gerechtfertigt ist. Die Dynamik des gestorten Systems ~ = ~ + ~ ergibt sich zu

@ ~ = F~  ~  + F  ~ , ~  ~(t = 0) = ~0 + ~(t = 0) = ~0; : 0  0 ~ @t Fur die Dynamik der Storung ~ = ~  , ~ erhalt man @  ~ , ~ = F  ~ , ~  ~(t = 0) , ~ (t = 0) = ~0; , ~0  ~ @t  @~ = F ~ ~(t = 0) = ~0 : ~ @t 0

(3.10)

0

0

Die Losung ergibt sich auf analytischem Weg als ~(t) = exp(F t)~0 : Fur diagonalisierbare Jacobimatrizen gilt

(3.11)

F ~ = UDU ,1 0

mit U = Matrixform der rechten Eigenvektoren (als Spaltenvektoren), U ,1 = Matrixform der linken Eigenvektoren (als Zeilenvektoren) zu den jeweiligen Eigenwerten der Jacobimatrix am Punkt ~ 0. Diese Eigenwerte (der Einfachheit halber seien sie alle nicht entartet) stehen in der Diagonalen von D. Fur (3.11) ergibt sich





~(t) = exp UDU ,1 t ~0 = U exp(Dt)U ,1~0 : Transformiert man das System nun in die Basis der Eigenvektoren, ~^ = U ,1~ ; so erhalt man wegen der Identitat U ,1 U = I U ,1~^ = U ,1 U exp(Dt)U ,1~0 = exp(Dt)U ,1~0 = exp(Dt)~^0 : In Komponenten geschrieben lautet das ^i = exp(it)^0;i i = 1; : : :; n : i sind die Eigenwerte der Jacobimatrix und ^i die Komponenten der Storung in der Basis der Eigenvektoren. Der Betrag der aufgepragten Storung entwickelt sich gema

  j^ij = exp Real t j^0;ij i

(3.12)

30

KAPITEL 3. AUTOMATISCHE REDUKTION VON REAKTIONSMECHANISMEN

mit Real i als Realteil des i-ten Eigenwerts. Die Deutung ist nun folgende: Das chemische System wird in Richtung des i-ten Eigenvektors der Jacobimatrix gestort. Das System reagiert in die entsprechende Richtung. Der Betrag des Realteils des zugehorigen Eigenwerts entscheidet uber die Zeitkonstante, mit der die Storung anwachst oder abgedampft wird. Das Vorzeichen des Eigenwerts erlaubt drei verschiedene Moglichkeiten der Entwicklung:  Real i > 0: Storung wachst exponentiell mit der Zeit.  Real i = 0: Storung bleibt zeitlich konstant.  Real i < 0: Storung wird exponentiell kleiner. Es sei an dieser Stelle noch angemerkt, da die Jacobimatrix fur diese Betrachtungen nicht notwendigerweise diagonalisierbar sein mu [48]. Dasselbe Ergebnis kann mit einer Jordanzerlegung ebenfalls erreichen [48,49]. Diese Analyse gibt ein anschauliches Bild der chemischen Dynamik, auch wenn diese stark vereinfacht dargestellt ist. Aber man sieht, da betragsmaig stark negative Eigenwerte fur eine schnelle Relaxation im Reaktionssystem verantwortlich sind. In Abb. 3.8 ist ein Eigenwertspektrum fur ein Knallgas-System mit 8 Spezies uber eine laminare vorgemischte ache Flamme aufgetragen. Lenkt man einen Punkt im Reaktionsraum aus, geben die Eigenwerte der Jacobimatrix Aussagen Eigenwerte (1/s) -10 -1

-10 4

-10

9

0.002

0.004

0.006

0.008

0.010

r (m)

Abbildung 3.8: Eigenwertspektrum der Knallgas-Reaktionen in einer vorgemischten achen laminaren Flamme. uber die zeitliche Entwicklung der Auslenkung. Fur negative Eigenwerte geht die Auslenkung in den zugehorigen charakteristischen Richtungen im Reaktionsraum (Eigenvektoren zu den jeweiligen Eigenwerten) zuruck. Je groer der Betrag der negativen Eigenwerte, desto schneller ist der Relaxationsproze.

3.3. INTRINSISCHE NIEDRIG-DIMENSIONALE MANNIGFALTIGKEITEN

31

3.3 Intrinsische niedrig-dimensionale Mannigfaltigkeiten 3.3.1 Einfuhrung

In diesem Kapitel soll nun die Methode der intrinsischen niedrig-dimensionalen Mannigfaltigkeiten eingefuhrt werden, um bestehende Reaktionsmechanismen, basierend auf einer mathematischen Analyse, zu vereinfachen. Zuvor aber noch ein paar Worte zu alternativen Vorgehensweisen. Die einfachste Moglichkeit einen vorhandenen Reaktionsmechanismus zu reduzieren ist naturlich die Vernachlassigung unwichtiger Spezies. Durch die Halbierung der Speziesanzahl erhalt man naturlich auch eine Reduzierung der CPU-Zeit bei der numerischen Simulation, die bei impliziten Verfahren ja immerhin mit dem Quadrat der Speziesanzahl eingeht. Durch so eine willkurlichen Eliminierung reaktiver Spezies vom Reaktionssystem konnen in bestimmten Bereichen aber auch unerwunschte Eekte auftreten, wie z.B. ungunstige Zeitschrittweitensteuerung oder eine Zunahme von Iterationsschritten. Auerdem hat man zwar die Anzahl der Spezies verkleinert, die Steifheit des partiellen Dierentialgleichungssystems aber bleibt, weswegen auch weiterhin implizite Methoden benotigt werden. Um die Steifheit der chemischen Kinetik zu reduzieren, mu man den Reaktionsmechanismus schon sehr genau kennen und analysieren, will man ihn mit konventionellen Methoden vereinfachen. Wie in Abb. 3.9 schematisch dargestellt ist, laufen die chemischen Prozesse auf stark Chemische Physikalische 1234 Zeitskalen langsam (NO-Bildung)

1234 1234 1234 1234 1234 1234

100s 10-2s 10-4s

schnell, part. Glgew., stat. Zustand

12345 12345 12345 12345 12345 12345 12345 12345 12345 12345 12345 12345 12345 12345 12345

Zeitskalen von molekularem Transport

10-6s 10-8s 10-10s ▼

entkoppeln

Abbildung 3.9: Schematische Darstellung der Zeitskalen der chemischen und physikalischen Prozess. unterschiedlichen Zeitskalen ab, die sich um mehrere Groenordnungen unterscheiden. Es gibt sehr schnelle Reaktionen im Nanosekundenbereich, wogegen die Schadstobildung oft sehr lang dauert. Die physikalischen Prozesse, wie molekularer Transport, Turbulenz, etc., sind dagegen auf einen sehr engen Bereich in der Groenordnung von Millisekunden beschrankt. Viele chemische Reaktionen sind aber sehr schnell. D.h., zu dem Zeitpunkt, an dem die physikalischen Terme dominant werden, sind die meisten chemischen Reaktionen schon lange im partiellen Gleichgewicht. Es nden zu diesem Zeitpunkt nur noch Bewegungen in ganz bestimmte Richtungen des Zustandsraumes statt, welche durch die langsamen chemischen Reaktionen bestimmt werden. Ein Hauptaugenmerk bei der Reduktion von Reaktionsmechanismen mu also auf die Bestimmung der schnellen chemischen Reaktionen gelegt werden, um diese dann in geeigneter Weise vom Dierentialgleichungssystem zu entkoppeln. Hat man dies erreicht, so reduziert man nicht nur die Anzahl der zu losenden Erhaltungsgleichungen, sondern auch ganz erheblich die Steifheit des Dierentialgleichungssystems. Um dies zu erreichen, wurden in den letzten Jahren verschiedene Ansatze modelliert. Zahlreiche reduzierte Reaktionsmechanismen gehen auf erste Arbeiten von Bo-

32

KAPITEL 3. AUTOMATISCHE REDUKTION VON REAKTIONSMECHANISMEN

denstein und Semenov [16, 17] zuruck. Sie basieren alle auf Quasistationaritatsannahmen und der Annahme partieller Gleichgewichte fur bestimmte Spezies bzw. bestimmter Reaktionen [14,18{23].

Partielle Gleichgewichte Das Reaktionssystem wird insgesamt durch die langsamen, geschwindigkeitsbestimmenden Prozesse bestimmt. Schnelle Reaktionen mit groen Quelltermen konnen dann schon im partiellen Gleichgewicht sein, wenn die Nettorate aus Hin- und Ruckreaktion sehr klein ist. Um diese schnellen chemischen Prozesse zu identi zieren, sucht man also nach Reaktionen, fur die 1 ,  > rrhin < 1 +  ruck

gilt, fur einen kleinen Toleranzwert . rhin; rruck sind die Reaktionsgeschwindigkeiten der Hin- bzw. der Ruckreaktion. Fur solche Reaktionen im partiellen Gleichgewich sind dann nur noch algebraische Gleichungen zu losen, da die an den Reaktionen beteiligten Spezieskonzentrationen uber das Massenwirkungsgesetz zusammenhangen. Der Grenzfall ist durch den chemischen Gleichgewichtspunkt gegeben, wo insgesamt nr Reaktionen im partiellen Gleichgewicht sind. In der Basis der Element- und Reaktionsvektoren, wie sie in Kap. 3.2.5 vorgestellt wurde, heit das, da alle Komponenten der Reaktionsgeschwindigkeiten in Richtung des zu der Reaktion gehorenden Reaktionsvektors verschwindet. Diese Methode der Mechanismenreduktion birgt aber einige Unsicherheitsfaktoren und Fehlerquellen in sich. Mit der Zeit bilden sich in der Flamme an verschiedenen Orten verschiedene Bedingungen aus. So ist es denkbar, da es Bereiche gibt, an denen fette Bedingungen vorherrschen, in benachbarten aber magere Verbrennung statt ndet. Beim U bergang von fetter zu magerer Verbrennung andert sich aber der Reaktionskanal (bei mageren und stochiometrischen Mischungen ndet z.B. der Abbau von Methan uber CH3 weiter zu den oxidierten Spezies statt. Unter fetten Bedingungen verschiebt sich der Reaktionskanal immer weiter hin zur C2 -Chemie). Auch andert sich die Bildung und Rekombination von Zwischenprodukten beim Durchgang uber eine Flammenfront. Hier durchlauft man einen groen Temperaturbereich, in dem fur die verschiedenen Temperaturen ganz unterschiedliche Reaktionen wichtig sein konnen. Die Annahme von globalen partiellen Gleichgewichten ist also fur diese stark unterschiedlichen Parametersatze, wie sie in technischen Verbrennungssystemen vorherrschen, global zu ungenau.

Quasistationaritat von Spezies Eine andere Moglichkeit, einen detaillierten Reaktionsmechanismus zu vereinfachen, ist, schnelle und langsame Zeitskalen als eine Funktion der Spezieskonzentrationen selbst zu betrachten. Die Radikalbildung und der Verbrauch der Zwischenprodukte ist ein sehr schneller Proze. Diese konnen in einer Gruppe von Reaktionen zusammengepat werden. Betrachtet man den Nettoquellterm fur ein Radikal, so ist er in der Regel sehr klein gegenuber dem Produktions- und Verbrauchswert der Elementarreaktionen. Setzt man nun Quasistationaritat bestimmter Spezies vorraus, so ist der Nettoquellterm dieser Spezies gleich Null. Reduzierte Mechanismen, die anhand der Annahmen der Quasistationaritat und partieller Gleichgewichte entwickelt wurden, sind schon weit verbreitet. Deren ausfuhrliche Diskussion und Anwendungen in Verbrennungssystemen konnen in [22,23] nachgelesen werden.

CSP Eine auf einem mathematischen Formalismus aufbauende Methode wurde von Lam und Goussis [24, 25] vorgeschlagen. Sie erlaubt die Entkopplung schneller Zeitskalen vom reaktiven System. Dabei wird die Steifheit des Dierentialgleichungssystems wesentlich reduziert. Die Zeitskalen der

3.3. INTRINSISCHE NIEDRIG-DIMENSIONALE MANNIGFALTIGKEITEN

33

reaktiven Variablen des detaillierten Reaktionsmechanismus werden uber eine mathematische Analyse klassi ziert und nach aufsteigendem Betrag sortiert. Durch den Anwender wird dann ein CutO-Wert fur die Zeitskalen festgelegt. Schnellere Zeitskalen als dieser Cut-O-Wert werden als abgelaufene schnelle Moden klassi ziert, die beteiligten Reaktionen als partielle Gleichgewichte und die Spezies als quasistationare Spezies identi ziert. Diese Prozedur wird lokal im Flammenprogramm durchgefuhrt. Das hat gegenuber den konventionellen Methoden der globalen Gleichgewichtsannahmen den Vorteil, da immer eine lokale Analyse des Reaktionsmechanismus statt ndet. Ein Nachteil ist, da aber das gesamte System der Dierentialgleichungen fur das skalare Feld mitgefuhrt werden mu. Es erfolgt also keine Reduktion der Dimension des Dierentialgleichungssystems.

ILDM Eine Methode, basierend auf der mathematische Analyse des Reaktionsmechanismus, analog der Prozedur der dynamischen Systeme, erlaubt eine Reduktion der Variablen, die die chemische Dynamik beschreiben [27, 28, 50, 30, 31, 33]. Hierbei wird der vieldimensionale Zustandsraum auf einen durch langsame Zeitskalen charakterisierten niedrig-dimensionalen Unterraum projiziert. Der niedrig-dimensionale Unterraum ist dabei allein durch den detaillierten Reaktionsmechanismus bestimmt, wobei die Richtung der Bewegung der chemischen Kinetik in diesem Unterraum nicht mehr unbedingt durch bestimmte Spezies oder Reaktionen beschrieben wird. Der mathematische Formalismus benutzt die Methoden der qualitative Theorie gewohnlicher Dierentialgleichungen. Hierbei wird das System in alle moglichen Richtungen des Zustandsraumes gestort, wobei die Reaktion des chemischen Systems auf die Storung analysiert wird, welche Aussagen uber die Zeitskalen der chemischen Prozesse und charakteristische Richtungen im Zustandsraum der chemischen Dynamik macht. Schnelle Prozesse konnen lokal im Zustandsraum identi ziert und vom Gleichungssystem entkoppelt werden. Somit erhalt man an jedem Punkt im Zustandsraum eine optimale Reduktion der chemischen Kinetik, da die entkoppelten schnellen Prozesse nicht mehr mit bestimmten chemischen Spezies oder Reaktionen assoziiert sind. Im folgenden Kapitel wird ausfuhrlich auf die Methode der intrinsischen niedrig-dimensionalen Mannigfaltigkeiten eingegangen. Die Methode soll sowohl anschaulich als auch mathematisch beschrieben werden. Anhand eines sehr einfachen analytischen Sytems wird die Methode dargelegt. Danach wird im Detail die Implementierung in verschiedene CFD-Anwendungen diskutiert. Hierbei wird der niedrig-dimensionale Unterraum als Funktion von Reaktionsfortschrittsvariablen in tabellarischer Form zur Verfugung gestellt. Die chemische Dynamik wird nur noch als Funktion dieser Reaktionsfortschrittsvariablen beschrieben. Der Zustand auf den niedrig-dimensionalen Unterraumen ist einzig abhangig von diesen Variablen. Die Validierung der Methode der intrinsischen niedrig-dimensionalen Mannigfaltigkeiten als reduzierte Mechanismen erfolgt dann in Kap. 4. Hier werden verschiedene chemische Reaktionssysteme reduziert und an unterschiedlichen technisch relevanten Simulationsbeispielen verwendet. Die Ergebnisse, die man mit der reduzierten Kinetik erhalt, werden dabei direkt mit Ergebnissen verglichen, die man mit detaillierten Reaktionsmechanismen erhalt.

3.3.2 Mathematische Formulierung der ILDM Die mathematische Ableitung zur Beschreibung der Dynamik der niedrig-dimensionalen Mannigfaltigkeiten im Zustandsraum geht von einem linearisierten System aus. Das Reaktionssystem im (n = 2 + ns )-dimensionalem Zustandsraum wird durch das Dierentialgleichungssystem

@ ~ = F~  ~  @t

(3.13)

34

KAPITEL 3. AUTOMATISCHE REDUKTION VON REAKTIONSMECHANISMEN

beschrieben. ~ = (h; p; 1; 2; : : :; ns )T ist der Zustandsvektor und F~ der Vektor der rechten Seiten. An jedem Punkt ~ im Zustandsraum ist die lokale Jacobimatrix F de niert als

0 @F BB @@F F =B BB @ .. @ .

1 CC C: .. C . C A

2 1

@F1 @ 2 @F2 @ 2









@F1 @ n @F2 @ n

@Fn @ 1

@Fn @ 2






@Fn @ n

1 1

.. .

...

Die (n n)-dimensionale Jacobimatrix besitzt n Eigenwerte 1; : : :; n mit n zugehorigen Eigenvektoren ~u1 ; : : :; ~un . Aus diesen Eigenvektoren soll nun eine alternative Basis im Zustandsraum konstruiert werden. Sind alle Eigenwerte i reell und nicht entartet, so sind die zugehorigen Eigenvektoren ~ui alle linear unabhangig voneinander. In diesem Fall bilden die Eigenvektoren eine alternative Basis im Zustandsraum. Fur komplexe oder entartete Eigenwerte lassen sich sogenannte verallgemeinerte Eigenvektoren konstruieren, die ebenfalls alle voneinander linear unabhangig sind und somit ebenfalls eine alternative Basis des Zustandsraumes bilden. Fur die folgenden Betrachtungen sollen die EigenReal Real werte nach kleiner werdenden Realteilen der Eigenwerte sortiert sein: Real 1  2  : : :  n . Real Real Im Im Fur komplexe und entartete Eigenwerte (i = i+1 ) gilt die Reihenfolge i  i+1. Die zugehorigen linken Eigenvektoren der verallgemeinerten Eigenvektoren ~ui seien mit ~uLi bezeichnet. Die Jacobimatrix lat sich mit Hilfe der Basis der verallgemeinerten Eigenvektoren zerlegen, wobei in der Matrix U die Eigenvektoren ~ui (als Spaltenvektoren) und in der inversen U ,1 die linken Eigenvektoren ~uLi (als Zeilenvektoren geschrieben) stehen. Es lat sich zeigen, da eine Jordanzerlegung der Jacobimatrix immer existiert [49],

0 L1 ~u1 B L C F = UJU ,1 = (~u1~u2 : : :~un) J B B@ :~u:2: CCA : u~ Ln

Die Matrix J besitzt dabei Jordansche Normalform, fur die gilt:

 Wenn alle Eigenwerte reell sind und nicht entartet, stehen in der Hauptdiagonalen der Jordanschen Matrix die Eigenwerte Jii = i und sonst uberall Nullen.  i; i+1 sind komplex konjugiert und nicht entartet, dann ist ! ! Real Im J J  i;i i;i +1 i i Ji = J = ,Im Real : i+1;i Ji+1;i+1 i i  Fur reelle Eigenwerte, die m-fach entartet sind, gilt 0 J Ji;i+1 Ji;i+2 : : : Ji;i+m,1 i;i B J J J B i+1;i i+1;i+1 i+1;i+2 : : : Ji+1;i+m,1 B J J J i+2;i i+2;i+1 i+2;i+2 : : : Ji+2;i+m,1 Ji = B B .. .. . . B ... .. @ .. . . Ji+m,1;i Ji+m,1;i+1 Ji+m,1;i+2 : : : Ji+m,1;i+m,1

1 0 CC BB 0i CC = BB 0 CC BB . A @ ..

1 0 i 1 0 i .. .. . . 0 0 0

 Komplexe Eigenwerte, m-fach entartet bilden sich nach obigen Konstruktionsregel.

1

::: 0 ::: 0 C C ::: 0 C CC . . . ... C A

: : : i

3.3. INTRINSISCHE NIEDRIG-DIMENSIONALE MANNIGFALTIGKEITEN

35

Fur die folgenden Betrachtungen soll es ausreichen, reelle, nicht entartete Eigenwerte zu betrachten. Die Jordanmatrix wird in diesem Fall zur Diagonalmatrix 1 0 1 0 : : : 0 BB .. C 0  . C CC : B 2 D=B ... CA B@ ... 0 : : : 0 n Somit ist die Basis der verallgemeinerten Eigenvektoren ~ui mit den linken Eigenvektoren ~uLi konstruiert. Das Dierentialgleichungssystem (3.13) sei nun im Punkt ~0 linearisierbar (Glieder der zweiten Ordnung und hoher werden in folgenden vernachlassigt),

@ ~ = F~  ~  + F  ~ , ~  ~ (t = 0) = ~0 : 0 0 ~ @t Mit der Abkurzung '~ = ~ , ~0 wird das zu @~' = F~  ~  + F '~ '~ (t = 0) = 0 : (3.14) 0 ~ @t Das lat sich in die Basis der verallgemeinerten Eigenvektoren U = (~u1 : : :~un ) transformieren. Mit '~ = U ~ wird (3.14) zu @ U ~ = U @ ~ @t @t~  ~ = F 0 + F ~ U ~ ~(t = 0) = 0 : 0

0

0



Multiplikation von links mit der Matrix der linken Eigenvektoren U ,1 = u~ L1 ; ~uL2 ; : : :; ~uLn (mit UU ,1 = I = Einheitsmatrix)

@ ~ = U ,1 F~  ~  + U ,1 F U ~ 0 ~ @t 0

T

ergibt

~(t = 0) = 0 :

Berucksichtigt man die Zerlegung der Jacobimatrix in der Basis der Eigenvektoren, wird das zu

@ ~ = U ,1 F~  ~  + D~ 0 @t

~(t = 0) = 0 :

Weiter vorne wurde festgelegt, da die Eigenwerte nach groer werdendem Betrag sortiert sein sollen. D.h., 1 ist der grote und n der am negativste Eigenwert. Das Dierentialgleichungssystem lat sich dann in zwei Systeme aufspalten. Das eine korrespondiert mit den langsamen Prozessen (entsprechend den Eigenwerten 1 bis n,nf ), das zweite mit den schnellen, welche durch die Eigenwerte n,nf +1 bis n charakterisiert werden,

@ ~s = U ,1 F~  ~  + D ~ ~s(t = 0) = 0 0 s s s @t @ ~f = U ,1 F~  ~  + D ~ ~f (t = 0) = 0 : f f 0 f @t Die Matrizen Us ; Uf ; Us,1 und Uf,1 entsprechen den Untermatrizen der Basis 0 1 0 1 j j j j j j Us = B @ ~u1 ~u2 : : : ~un,nf CA Uf = B@ ~un,nf +1 ~un,nf +2 : : : ~un CA : j j j j j j

(3.15) (3.16)

36

KAPITEL 3. AUTOMATISCHE REDUKTION VON REAKTIONSMECHANISMEN

0 ~uL BB ~u1L2 Us,1 = B B@ ...

~uLn,nf

1 CC CC A

0 ~uL BB ~unLn,,nnf +1 f +2 Uf,1 = B B@ ... ~uLn

1 CC CC : A

Die Diagonalmatrix ist analog in zwei Anteile Ds und Df unterteilt,

0 BB 1 2 Ds = B @


n,nf 0 BB n,nf +1 n,nf +2 Df = B @

1 CC CA

(n,nf )(n,nf )






n

1 CC CA

: (nf )(nf )

Die Indizes s und f stehen hierbei fur s=langsam (slow) und f=schnell (fast). Ds charakterisiert also die Zeitskalen der (n , nf ) langsamsten Prozesse im langsamen Unterraum Us , Df die der schnellen im Unterraum der schnellen Relaxationsprozesse Uf . Die Basis, die den Raum aufspannt, in dem die Dynamik der langsamen Bewegungen statt ndet, ist gegeben durch die (ersten) verallgemeinerten Eigenvektoren, die zu den langsamen Eigenwerten 1; : : :; n,nf gehoren. Wie in Kap. 3.2 schon anschaulich gezeigt (vergleiche Abb. 3.1 bis Abb. 3.6) relaxieren die Trajektorien der chemischen Dynamik sehr schnell auf niedrig-dimensionale Unterraume, auf denen die weitere Bewegung dann langsam weitergeht, bis das chemische Gleichgewicht erreicht ist. Ist man an den Prozessen der Relaxation nicht interessiert, sondern genugt die Dynamik in einem niedrigdimensionalen Unterraum (sagen wir auf einer (m = n , nf )-dimensionalen Mannigfaltigkeit), so gelten fur die nf schnellsten Prozesse die folgenden Gleichgewichtsbedingungen

 

~uLi F~ ~ = 0

i = n , nf + 1; : : :; n :

(3.17)

Das bedeutet, da die nf schnellsten Prozesse als relaxiert betrachtet werden. Der chemische Quellterm F~ in Richtung der schnellsten Zeitskalen ~un,nf +1 ; : : :; ~un verschwindet. Die Gleichungssysteme (3.15) und (3.16) werden zu

@ ~s = U ,1F~  ~  + D ~ ~s(t = 0) = 0 0 s s s @t @ ~f = D ~ ~f (t = 0) = 0 : f f @t

Mit der Annahme, da die schnellsten Zeitskalen im lokalen Gleichgewicht sind, gilt deshalb

und damit

@ ~f = 0 ; @t ~f (t) = 0 :

Die Dynamik des chemischen Systems beschrankt sich also auf die Dynamik der langsamsten Prozesse, charakterisiert durch die Eigenwerte 1 : : :n,nf . Die schnellen Prozesse sind schon ins lokale Gleichgewicht relaxiert. Findet man im Zustandsraum Punkte fur die die Gleichgewichtsbedingungen (3.17) erfullt sind, hat man die Punkte auf der intrinsischen niedrig-dimensionalen Mannigfaltigkeit identi ziert.

3.3. INTRINSISCHE NIEDRIG-DIMENSIONALE MANNIGFALTIGKEITEN

37

Um die niedrig-dimensionalen Unterraume im Zustandsraum zu nden, genugt es also, eine Eigenwertanalyse der lokalen Jacobimatrizen durchzufuhren. Die Mannigfaltigkeit ist, wenn die obigen Gleichgewichtsbedingungen erfullt sind, automatisch durch die Anzahl der gewunschten Freiheitsgrade n , nf de niert. In realen Systemen treten aber im allgemeinen komplexe und entartete Eigenwerte auf, was eine schlecht konditionierte Jacobimatrix zur Folge hat. Bei der numerischen Realisierung der Mannigfaltigkeiten ergeben sich dadurch Schwierigkeiten, so da die Mannigfaltigkeiten oft nicht identi ziert werden konnen. Es ist dadurch sinnvoll, eine alternative Basis zur Beschreibung der Unterraume heranzuziehen. Hiefur bietet sich die Basis der modi zierten reellen Schur-Vektoren an.

Alternative Basis: Schur-Vektoren In realen Systemen ist die Jordan-Zerlegung infolge nahezu entarteter Eigenwerte im allgemeinen schlecht konditioniert, was die numerische Beschreibung in der Eigenvektor-Basis erschwert. Da bei der numerischen Realisierung der Methode die Schur-Vektoren verwendet werden, soll die Transformation in deren Basis vorgestellt werden. Es lat sich zeigen, da eine Schur-Zerlegung der Jacobimatrix immer existiert. Diese spannt denselben invarianten Unterraum der langsamen Zeitskalen auf, wie die verallgemeinerte Eigenvektorbasis Us . Um die numerischen Schwierigkeiten zu umgehen, die aus der Existenz entarteter und komplexer Eigenwerte entstehen, wird die Beschreibung der Dynamik auf den niedrig-dimensionalen Mannigfaltigkeiten in der Basis der Schur-Vektoren durchgefuhrt. Deshalb soll hier gezeigt werden, da die obige Beschreibung der niedrig-dimensionalen Mannigfaltigkeiten basierend auf den Eigenvektoren aquivalent zur Beschreibung mit den reellen Schur-Vektoren ist. Analog zur Zerlegung in der Basis der Eigenvektoren existiert eine Schurzerlegung der Jacobimatrix, !   T Q QT F Q = QsT F Qs Qf = Z mit

f

0 ! BB z011 zz12 22 Z = Z0s ZZ2 = B . B f @ ..

1

: : : z1n : : : z2n C CC

.. C : . A

...

0 0 : : : znn Wie vorher angedeutet, spannen die Schur-Vektoren ~qi ; i = 1; : : :; n , nf denselben invarianten Raum auf wie die Eigenvektoren ~ui ; i = 1; : : :; n , nf . Der Unterraum der langsamen Zeitskalen (Index s) wird also durch die ersten n , nf Schur-Vektoren beschrieben: ~q1 ; : : :; ~qn,nf . Die SchurVektoren konnen direkt aus der Basis der Eigenvektoren bestimmt werden [48],

 und Daraus folgt

Qs Qf Qs Qf

!

 A B! = Us Uf 0 C

 

^ = A 0

B^ C^

!

^ f,1 : QTf = CU

Us,1 Us,1

!

:

Aus den Gleichgewichtsbedingungen (3.17) fur die schnellen Zeitskalen Uf,1 F = 0 folgt in der Basis der Schur-Vektoren ^ f,1 F = QTf F = 0 : CU (3.18)

38

KAPITEL 3. AUTOMATISCHE REDUKTION VON REAKTIONSMECHANISMEN

Betrachtet man die Prozedur genauer, erkennt man, da die einzige Information, die man benotigt, um die Mannigfaltigkeit zu identi zieren, die Dimension der Mannigfaltigkeit selbst ist. Der detaillierte Reaktionsmechanismus wird analysiert und die Mannigfaltigkeiten, basierend auf einer Eigenwertanalyse, uber die notwendigen Gleichgewichtsbedingungen (3.18) identi ziert.

Eigenschaften In Kap. 3.3.2 wurde die mathematische De nition der intrinsischen niedrig-dimensionalen Mannigfaltigkeiten, basierend auf der Eigenwertanalyse der lokalen Jacobimatrizen im Zustandsraum, ausfuhrlich diskutiert. Daraus folgen bestimmte Eigenschaften der Unterraume, die hier dargestellt werden sollen.

Bewegung hin zur Mannigfaltigkeit: Die Frage, wie sich ein chemisches System verhalt,

wenn es sich in einem Zustand be ndet, der nicht auf der Mannigfaltigkeit liegt, ist von sehr groer Bedeutung. Es wird noch gezeigt, da Prozesse existieren, die einen Zustand dergestalt storen, da dieser von der Mannigfaltigkeit wegbewegt wird. Auch wird man bei der numerische Simulation von Stromungsprozessen nicht immer optimale Anfangs- oder Randbedingungen angeben konnen, so da der Startpunkt einer Berechnung eventuell auerhalb der Mannigfaltigkeit liegt. Da dies keine prinzipiellen Einschrankungen der Methode darstellt, wird im folgenden gezeigt. Das chemische System wird durch die Dierentialgleichungen @ ~ =@t = F~ ( ~ ) beschrieben. Der Punkt ~0 sei ein Punkt auf der Mannigfaltigkeit und soll durch eine kleine Storung aus dem Unterraum fortbewegt werden. Um den Ein u der Storung auf das chemische Reaktionssystem betrachten zu konnen, wird das Dierentialgleichungssystem um ~0 linearisiert,

@ ~ = F~  ~  + F  ~ , ~  0 0 ~ @t

~ (t = 0) = 0 ;

0

mit '~ = ~ , ~ 0

@~' = F~  ~  + F '~ : 0 ~ @t 0

Nun wird das Dierentialgleichungssystem von links mit der transponierten Teilmatrix der SchurVektoren QTf multipliziert, wobei QTf die Untermatrix der transponierten Schur-Vektoren ist, die zu den schnellen Zeitskalen gehoren,

0 ~qT 1 n,nf +1 BB ~qnT,n +2 CC f QTf = B B@ ... CCA :

Das DGL-System wird dann zu

q~nT

@ QT'~  = QTF~  ~  + QT F '~ : 0 ~ f f @t f 0

(3.19)

~0 auf der Mannigfaltigkeit liegt, gelten fur ihn die Gleichgewichtsbedingungen Da der Punkt   QTf F~ ~ 0 = 0 (siehe (3.18)) und das System wird zu

@ QT '~  = QT F '~ : ~ f @t f Unter Verwendung der Schur-Zerlegung QT F Q = Z mit der oberen Dreiecksmatrix Z folgt fur die zeitliche A nderung von '~ @ QT'~  = QT QZQT'~ : f @t f 0

3.3. INTRINSISCHE NIEDRIG-DIMENSIONALE MANNIGFALTIGKEITEN

39

Mit den folgenden Umformungen und m = n , nf

0 ~qT 1 m+1 T C B q ~ QTfQ = B C ( ~q1 ~q2 : : : ~qm ~qm+1 : : : ~qm+nf ) B@ m..+2 C . A q~nT

= ( 0mn Inf n )(nnf )

QTf QZ = ( 0(mn) I(nf n) )(nnf ) =

QTf QZQT

=





0(mnf ) Zf nf nf (

0(mnf ) Zf nf nf (



)

(

(

) )

(nn)

QTsnm ! QTfnnf (nn) (

(nnf )

)

)

(

(nnf )



!

Zs mm Z2 nf m 0(mnf ) Zf nf nf

(

ZfT QTf

= ; erhalt man fur das DGL-System der Storung,

)

)

@ QT '~  = Z T QT'~ : f f @t f Dies lat sich integrieren und man erhalt mit '~ (t = 0) = '~ 0 QTf '~ (t) = eZ tQTf '~ 0 : Aus der De nition der Matrix Z und der Sortierung der Eigenwerte erkennt man, da in Zf die Eigenwerte der schnellen Prozesse stehen, d.h. Real  0. Damit folgt fur die Storung, da sie i f

sich immer sehr schnell und exponentiell der Mannigfaltigkeit nahert. Dies gilt mit einer Zeitskala, die dem langsamsten Eigenwert von Zf entspricht, d.h. mit der Zeitskala, des ersten entkoppelten Eigenwerts s = m+1 . Wird also ein Zustand durch einen physikalischen Proze von der Mannigfaltigkeit heruntergebracht, so relaxiert der gestorte Zustand mit   '~ (t) = QTf F~ ~  Konst: e t '~ (t = 0) auf die Mannigfaltigkeit zuruck. Der Sachverhalt wird wichtig, wenn die Sprache auf die Kopplung der reduzierten Kinetik mit physikalischem Transport (siehe Kap. 3.3.4) kommt. Diusion will namlich den Zustand von der Mannigfaltigkeit fortbewegen. Schnelle Relaxationen fuhren die Storungen aber immer wieder auf die Mannigfaltigkeiten zuruck. s

Bewegung auf der Mannigfaltigkeit: Be ndet man sich auf der Mannigfaltigkeit, wird die

Dynamik nur noch durch die langsamen Prozesse kontrolliert. Die schnellen Relaxationsprozesse beschranken die Bewegung nur auf die Unterraume. Aber bleibt das System, das auf die Mannigfaltigkeit relaxiert ist, auch  ~auf ihr bis das chemische Gleichgewicht erreicht ist? Das ist der Fall, wenn ~  die Anderung von F immer tangential zur Mannigfaltigkeit ist. Mathematisch bedeutet das, da die Mannigfaltigkeit ein invarianter Unterraum des Reaktionsraumes sein mu. Die Bedingung ist fur lineare Mannigfaltigkeiten erfullt. Fur nichtlineare Mannigfaltigkeiten liegt der Fall etwas komplizierter. Es lat sich aber mathematisch zeigen, da eine invariante Mannigfaltigkeit vorliegt, wenn das Verhaltnis der schnellen Zeitskalen zu den langsamen gegen unendlich geht. Fur reale Verbrennungssysteme ist das nie erfullt. Da aber dennoch sehr groe Unterschiede zwischen den langsamen und schnellen Reaktionen existieren (in der Regel mehrere Groenordnungen), ist die Bedingung der invarianten Mannigfaltigkeiten fur praktische Anwendungen hinreichend gut erfullt.

40

KAPITEL 3. AUTOMATISCHE REDUKTION VON REAKTIONSMECHANISMEN

Die Dimensionen der Mannigfaltigkeiten: Aus der De nitionsgleichung der Mannigfaltig  T~ ~

keit Qf F = 0 folgt, da niedrig-dimensionale Mannigfaltigkeiten immer Teilmengen der um eine Dimension erhohten Mannigfaltigkeiten sind. D.h., die nulldimensionale Mannigfaltigkeit ist Teilmenge der eindimensionalen Mannigfaltigkeit, diese wiederrum ist Teilmenge der zweidimensionalen, etc. Das ist sehr schon an den Trajektorien in Abb. 3.1 und Abb. 3.2 zu sehen. Die Trajektorien erreichen zunachst eine Flache (zweidimensionale Mannigfaltigkeit) im Reaktionsraum bevor sie sich auf dieser Flache in eine Linie (eindimensionale Mannigfaltigkeit) bundeln und schlielich im Gleichgewichtspunkt (nulldimensionale Mannigfaltigkeit) enden. Die Gleichgewichts-Chemie und der detaillierte Reaktionsmechanismus stellen also die "Eckpunkte\ der Theorie dar. Die reduzierte Chemie lat einen kontinuierlichen U bergang zwischen diesen Extremfallen zu.

Eindeutigkeit und Existenz der Mannigfaltigkeiten: Meistens wei man aus physikalisch-

chemischen U berlegungen, ob und wo die Mannigfaltigkeiten existieren. Das ist wiederrum sehr anschaulich uber die Trajektorienberechnung (siehe Abb. 3.1 und Abb. 3.2) zu sehen. Da an jedem Punkt im Reaktionsraum die lokale Jacobimatrix existiert und es zu jeder Jacobimatrix eine Schur-Zerlegung gibt, vereinfacht sich die Untersuchung auf Existenz und  ~Eindeutigkeit der ManT ~ nigfaltigkeit. Findet man eine Losung fur das Gleichungssystem Qf F = 0, existiert auch die niedrig-dimensionale Mannigfaltigkeit. Im Rahmen dieser Arbeit soll es dann auch ausreichen, die Mannigfaltigkeiten auf numerischem Weg zu bestimmen, ohne extra deren Existenz nachzuweisen.

Stetigkeit und Dierenzierbarkeit der Mannigfaltigkeiten: Stetigkeit und Dierenzierbar-

keit der invarianten Unterraume ist immer dann erfullt, wenn sich die Eigenwerte der Jacobimatrix nicht kreuzen. D.h., ein schneller Proze bleibt schnell und ein langsamer bleibt langsam. Ist das der Fall, ist die Jacobimatrix beliebig oft dierenzierbar und die Eigenvektoren andern sich kontinuierlich. Schwierigkeiten erhalt man, wenn n,nf = n,nf +1 wird. In diesem Fall andern sich die De nitionsgleichungen der niedrig-dimensionalen Mannigfaltigkeit. Aber auch hier soll es im Rahmen dieser Arbeit genugen, die Stetigkeit und Dierenzierbarkeit wahrend der numerischen Bestimmung (siehe Kap. 3.3.3) der Mannigfaltigkeiten zu uberprufen.

Analytisches Beispiel

Betrachtet werden soll in einem einfachen aber doch genug komplexen Beispiel die Methode der intrinsischen niedrig-dimensionalen Mannigfaltigkeiten. Das Reaktionssystem sei gegeben durch: X

a1 ,

Y

,

@@Ra2 b - Z c

Das Dierentialgleichungssystem ist

@x = , (a + a ) x 1 2 @t @y = a x , by + cz 1 @t @z = a x + by , cz 2 @t

und in vektorielle Schreibweise

@ ~ = F~  ~ @t

x(t = 0) = x0 y (t = 0) = y0 z (t = 0) = z0 ;

3.3. INTRINSISCHE NIEDRIG-DIMENSIONALE MANNIGFALTIGKEITEN

0x1 ~ = @yA

0 , (a1 + a2) x 1 F = @ a1x , by + cz A : z a2x + by , cz 0 ,a1 , a2 0 0 1 J = @ a1 ,b c A a2 b ,c

Die Jacobimatrix ist mit den Eigenwerten den rechten Eigenvektoren

1 = 0

  ~ ~

2 = ,b , c

3 = ,a1 , a2 ;

001 001 0 ,a ,a +b+c 1 a ,b ~u1 = @ cb A ~u2 = @ ,1 A ~u3 = @ aa ,,cb A 1

1

2

2 1 2

1

1

und den linken Eigenvektoren

~uL1 = ( b+b c b+b c b+b c ) , ~uL2 = (b+c)(a,ab,,aa c+b+c) b,+bc ~uL3 = ( ,a ,a a,+b b+c 0 0 ) : 1

1

c
b+c

2 2

2

1

Die Diagonalmatrix erhalt man zu

2

00

D = U ,1JU = @ 0 0

0

,b , c 0

0 0

,a1 , a2

1 A:

Die Transformation in die Eigenvektorbasis U ~ = ~ liefert

0 1 0 ~_ = U ~_ = UJ ~ = UJU ,1~ = D~ = @ , (a1 + a2) 2 A , (b + c) 3

mit der analytische Losung in der Eigenvektorbasis

0 x0 1 ~(t) = @ y0e,(b+c)t A :

z0 e,(a +a )t Rucktransformation in die ursprungliche Basis U ~ = ~ ergibt 1

2

1 0 x(t) 1 0 ,a ,a +b+c z0 e,(a +a )t a ,b ~ (t) = @ y (t) A = B @ cb x0 + y0e,(b+c)t , aa ,,cb z0e,(a +a )t CA : 1

2

1

2

2

z(t)

1 2

1

2

x0 + y0e,(b+c)t + z0e,(a1 +a2 )t

Die Integrationskonstanten ergeben sich aus

0 x0 1 @ y0 A = U ~0 :

z0 Das chemische Gleichgewicht stellt sich fur t ! 1 ein, x(t ! 1) = 0 y (t ! 1) = b +c c

z (t ! 1) = b +b c :

41

42

KAPITEL 3. AUTOMATISCHE REDUKTION VON REAKTIONSMECHANISMEN z(t)

y(t) 0.7

0.7

0.6

0.6

0.5

0.5

0.4

0.4

0.3

0.3

0.2

0.2

0.1

0.1

0

0 0

0.2

0.4

0.6

0.8

1

0

0.2

0.4

0.6

0.8

1

0.8

1

x(t)

x(t) z(t)

x(t), y(t(, z(t)

0.7

1

0.6 0.8

0.5 0.6

0.4 0.3

0.4

0.2 0.2

0.1 0

0

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0

0.2

y(t)

0.4

0.6 t(s)

Abbildung 3.10: Trajektorien des einfachen analytischen Systems. In Abb. 3.10 sind einige Trajektorien in den verschiedenen Projektionen der Variablen fur die Geschwindigkeitskonstanten a1 = 10; a2 = 15; b = 1; c = 2 abgebildet. Zusatzlich ist fur eine Anfangsbedingung (x(0) = 1; y (0) = z (0) = 0) die zeitliche Entwicklung der drei Variablen zu sehen. Die Gleichung der niedrig-dimensionalen Mannigfaltigkeit erhalt man durch Multiplikation der schnellen Untermatrix Uf,1 der linken Eigenvektoren (die zu den schnellen Zeitskalen gehoren) mit der rechten Seite des Ausgangsgleichungssystems. In diesem Beispiel soll nur die dritte Zeitskala entkoppelt werden. Das bedeutet, da a1 + a2  b + c sein soll: 0 ,(a1 + a2)x 1   ~uL3 F~ ~ = ( 1 0 0 ) @ a1x , by + cz A =! 0 ; a2x + by , cz woraus man die Bedingung fur die Mannigfaltigkeit erhalt:

x=0: Mit der Normierungsbedingung (alle Massenbruche addieren sich zu Eins) x + y + z = 1 erhalt man fur die Mannigfaltigkeit die Gleichungen

x = 0 und y = 1 , z : In Abb. 3.10 kann man erkennen, da die verschiedenen Trajektorien genau auf diese Kurve hinrelaxieren. Damit hat man die eindimensionale Mannigfaltigkeit im dreidimensionalen Reaktionsraum der Variablen identi ziert.

3.3. INTRINSISCHE NIEDRIG-DIMENSIONALE MANNIGFALTIGKEITEN

43

Beispiel: Methanoxidation Im obigen Beispiel war die Dynamik absichtlich so einfach gewahlt, so da die Methode der intrinsischen niedrig-dimensionalen Mannigfaltigkeiten anschaulich und analytisch darstellbar ist. Fur reale Verbrennungssysteme ist die Gesamtdynamik aber weitaus komplexer. Die Reaktionsmechanismen beschreiben die Kinetik in einem sehr hoch-dimensionalen Zustandsraum. Die Variablen der Spezieskonzentrationen sind auf sehr komplexe nichtlineare Art und Weise miteinander verknupft. Das hier angefuhrte Beispielsystem der Methanoxidation in Luft beinhaltet fur den detaillierten Reaktionsmechanismus 34 Spezies und 288 Reaktionen. Der Mischungsraum ist also 34-dimensional. Insgesamt kommen 4 verschiedene Elemente vor, so da fur konstante Enthalpie, Druck und Elementzusammensetzung der Reaktionsraum 30-dimensional ist. Fur den Fall der Gleichgewichtschemie (gemischt gleich verbrannt) sind alle 30 verschiedenen Zeitskalen im lokalen Gleichgewicht. Fur 29 Zeitskalen im Gleichgewicht erhalt man eine eindimensionale Mannigfaltigkeit und fur nf = 28 eine zweidimensionale.

Eindimensionale Mannigfaltigkeit: In Abb. 3.11 und Abb. 3.12 ist eine eindimensionale Man-

nigfaltigkeit fur die Methanoxidation abgebildet. Zusatzlich sind verschiedene Trajektorien (vergleiche Kap. 3.2, Abb. 3.1 und Abb. 3.2) mit eingezeichnet. Da die Dimension des Reaktionsraumes fur eine graphische Darstellung zu hoch ist, wurden wieder Projektionen in verschiedene Ebenen des Phasenraumes benutzt. Der Gleichgewichtswert ist wieder durch den Kreis gekennzeichnet und als Teil der eindimensionalen Mannigfaltigkeit sichtbar. Man erkennt sofort, da die Mannigfaltigkeit einen Attraktor fur die Trajektorien darstellt, wie schon in Kap. 3.2 beobachtet wurde und in Abb. 3.11 und Abb. 3.12 fur verschiedene stabile, aber auch reaktive Spezies zu sehen ist. Das ist ein wichtiger Aspekt, wenn man Schadstoe berechnen will. Um z.B. die NO-Bildung in Flammen zu berechnen, ist die O-Konzentration ein wichtiger Faktor. Bei guter Vorhersage fur die O-Atome kann man sehr gute Aussagen uber die NO-Bildungsgeschwindigkeit machen, da diese (Quasistationaritat der N-Atome vorrausgesetzt) nur noch eine Funktion der atomaren Sauerstokonzentration ist. In Abb. 3.13 sind zwei Eigenwerte uber die Reaktionsfortschrittsvariable wCO der eindimensionalen Mannigfaltigkeit abgebildet. Die untere Linie ist der langsamste (betragsmaig kleinste) entkoppelte Eigenwert. Er gibt die Zeitskala ( = 1=) an, mit der die Relaxationsprozesse die Mannigfaltigkeit erreichen. Dieser ist uberall kleiner als ,10+4 s, uber weite Strecken sogar um die ,10+5 s. Dies wird erwartet, wenn man mit den Aussagen aus Kap. 3.2 vergleicht, wo man die ersten 50s der Trajektorien einfach abgeschnitten hatte (siehe Abb. 3.5). Zu diesem Zeitpunkt sind die Trajektorien schon alle auf die eindimensionale Linie hinrelaxiert. Die Gute der Mannigfaltigkeit lat sich aus der Dierenz der abgebildeten Eigenwerte ablesen. Die obere Linie in Abb. 3.13 stellt den Eigenwert dar, der die schnellste Bewegung entlang der Mannigfaltigkeit charakterisiert. Die untere ist dagegen die Zeitskala des Relaxationsprozesses auf die Mannigfaltigkeit, also der langsamste entkoppelte Eigenwert. Je weiter diese beiden Eigenwerte voneinander entfernt sind, desto besser ist die Annahme der Entkopplung der schnellen Prozesse. 2

Zweidimensionale Mannigfaltigkeit: Liegen die oben angesprochenen Eigenwerte zu dicht

beieinander, ist die Entkopplung der schnelleren der beiden Zeitskalen nicht mehr erlaubt. Um die Genauigkeit der Simulation reaktiver Stromungen zu erhohen mu man fur die Beschreibung der chemischen Kinetik einen zusatzlichen Freiheitsgrad betrachten. Man benutzt dann eine zweidimensionale Mannigfaltigkeit, entsprechend einer zweidimensionalen Flache im Reaktionsraum. Der langsamste entkoppelte Eigenwert der eindimensionalen Mannigfaltigkeit ist jetzt kein schneller Proze mehr, sondern die zugehorigen Prozesse werden in die chemische Dynamik mit einbezogen. Am Beispiel der Methanoxidation hat es sich als gunstig erwiesen, die Spezieskonzentrationen von CO2 und H2 O als Reaktionsfortschrittsvariablen zu wahlen. Alle anderen Spezies sowie die Tem-

KAPITEL 3. AUTOMATISCHE REDUKTION VON REAKTIONSMECHANISMEN

44

wC O

wH 2

0.05

0.006

0.04 0.004

0.03

0.02 0.002

0.01

0 0.06

0.08

0.1

0.12

0.14

0.16

0 0.06

0.08

0.1

wC O 2

0.12

0.14

0.16

wC O 2 wH 2

wO 2

0.006

0.07

0.004 0.035 0.002

0 0.06

0 0.08

0.1

0.12

0.14

0

0.16

0.035

wC O 2

0.07

wO 2

wH 2 O

wC O 0.05

0.13 0.12

0.04

0.11 0.03 0.1 0.02 0.09 0.01

0.08 0.07 0.06

0.08

0.1

0.12 wC O 2

0.14

0.16

0 0.07

0.08

0.09

0.1

0.11

0.12

0.13

wH 2 O

Abbildung 3.11: Eindimensionale Mannigfaltigkeit im Reaktionsraum fur das stochiometrische CH4 -O2 -N2 -System. Gestrichelte Linien sind Trajektorien,  kennzeichnet den chemischen Gleichgewichtspunkt.

3.3. INTRINSISCHE NIEDRIG-DIMENSIONALE MANNIGFALTIGKEITEN

wO

45

wH

0.005

0.0009

0.004 0.0006

0.003

0.002 0.0003

0.001

0 0.06

0.08

0.1

0.12

0.14

0 0.06

0.16

0.08

0.1

wC O 2

0.12

0.14

0.16

0.003

0.004

0.005

wC O 2

wC H O

wO H

7 10-7

0.006

6 10-7 5 10-7 0.004 4 10-7 3 10-7 0.002

2 10-7 1 10-7

0

0 0

0.001

0.002

0.003

0.004

0

0.005

0.001

0.002

wO

wO

T(K)

T(K)

2200

2200

2000

2000

1800

1800

1600

1600

0.06

0.08

0.1

0.12 wC O 2

0.14

0.16

0.07

0.08

0.09

0.1

0.11

0.12

0.13

wH 2 O

Abbildung 3.12: Eindimensionale Mannigfaltigkeit im Reaktionsraum fur das stochiometrische CH4 -O2 -N2 -System. Gestrichelte Linien sind Trajektorien,  kennzeichnet den chemischen Gleichgewichtspunkt.

46

KAPITEL 3. AUTOMATISCHE REDUKTION VON REAKTIONSMECHANISMEN λ (1/s) -10

-100

-1000

-104

-105 0.1

0.11

0.12

0.13

0.14

0.15

wC O 2

Abbildung 3.13: Eigenwerte des stochiometrischen Methan-Luft-Gemisches einer eindimensionalen Mannigfaltigkeit uber die Reaktionsfortschrittsvariable wCO . Die untere Linie entspricht dabei dem langsamsten entkoppelten Eigenwert, die obere dem schnellsten der Prozesse, die die Bewegung auf der Mannigfaltigkeit beschreiben. 2

peratur und chemischen Quellterme sind somit eindeutige Funktionen von CO2 und H2 O . Das ist graphisch in dreidimensionalen Schaubildern in Abb. 3.14 und Abb. 3.15 dargestellt. Die gekrummten Flachen stellen die zweidimensionale Mannigfaltigkeit im Reaktionsraum fur eine konstante Enthalpie, Druck und Elementzusammensetzung dar. Aufgetragen sind verschiedene Massenbruche in Abhangigkeit der beiden Reaktionsfortschrittsvariablen. Zusatzlich sind noch einige Trajektorien eingezeichnet, die irgendwo auf der Mannigfaltigkeit beginnen. Man kann sehr schon erkennen, da die Trajektorien, die auf der zweidimensionalen Flache verlaufen, sich auf einer Linie zusammenbundeln und dann schlielich im Gleichgewichtspunkt (Kreis) enden. Auerdem wurden der schnellste nicht entkoppelte Eigenwert und der langsamste entkoppelte Eigenwert uber der Mannigfaltigkeit aufgetragen (siehe Abb. 3.16). Im Bereich groen Abstands produziert die Mannigfaltigkeit gute Ergebnisse. Im kalten Reaktionsgebiet (bei wCO = wH O = 0 hat man die Temperatur von 300 K) ist die chemische Kinetik sehr langsam. Streng genommen reagieren hier Methan und Luft erst nach sehr langer Zeit. Es ist einzusehen, da die Vorraussetzungen zur Berechnung der Mannigfaltigkeiten hier nicht gegeben sind. An den vorgestellten Graphiken lassen sich einige Eigenschaften der niedrig-dimensionalen Mannigfaltigkeiten festmachen. Zum einen sind sie gekrummte nichtlineare Hyper achen im Reaktionsraum. Das ist klar, da die chemische Kinetik extrem nichtlineare Terme in den Erhaltungsgleichungen besitzt. Zweitens erkennt man in den Pro len in Abb. 3.14 und Abb. 3.15 lokale Minima und Maxima. Diese deuten auf Wechsel in den Reaktionskanalen hin, d.h. andere Reaktionen werden bei A nderung der Reaktionsfortschrittsvariablen geschwindigkeitsbestimmend. An solchen Punkten kreuzen sich im allgemeinen die Eigenwerte. Da die Methode der intrinsischen niedrig-dimensionalen Mannigfaltigkeiten auf einer lokalen Analyse im Reaktionsraum basiert, wird solchen Eekten Rechnung getragen, was bei konventionell reduzierten Reaktionsmechanismen ublicherweise nicht der Fall ist, da diese nur einer globalen Betrachtung unterzogen werden. 2

2

3.3. INTRINSISCHE NIEDRIG-DIMENSIONALE MANNIGFALTIGKEITEN

wCO

47

wO 6.0E-3

8.0E-2

4.0E-3

4.0E-2

2.0E-3

0.0E0

0.0E0

-4.0E-2

0.01 0.06

0.01

-8.0E-2

0.01 0.01

0.11

wH2O

0.11

wH2O

0.05

0.06

0.11

wH2O 0.09

0.05 0.09

0.13

0.13

wCO2

wCO2 wH2

wH

1.2E-2

2.0E-3

8.0E-3

1.0E-3

4.0E-3 0.0E0

0.0E0

0.01 0.06

0.01

0.11

0.01 0.06

0.01

wH2O

0.05

0.05

0.09

0.09

0.13

0.13

wCO2

wCO2 wO2

wOH 6.0E-3

2.0E-1

4.0E-3

1.6E-1

2.0E-3 0.0E0 0.01

1.2E-1

0.01

8.0E-2

0.06 0.11

wH2O

4.0E-2 0.0E0

0.05

0.01 0.09

0.01

0.06

0.11

wH2O

0.05 0.09

0.13 0.13

wCO2

wCO2

Abbildung 3.14: Zweidimensionale Mannigfaltigkeit im Reaktionsraum fur das stochiometrische CH4 -O2 -N2 -System. Zusatzliche Linien sind Trajektorien,  kennzeichnet den chemischen Gleichgewichtspunkt.

48

KAPITEL 3. AUTOMATISCHE REDUKTION VON REAKTIONSMECHANISMEN

wH2O2

wHO2

6.0E-5

4.0E-5 3.0E-5

4.0E-5

2.0E-5

2.0E-5

1.0E-5 0.0E0

0.0E0

0.01 0.06

0.01

0.11

wH2O

0.05

0.01 0.06

0.01

0.11

wH2O

0.05

0.09

0.09

0.13

0.13

wCO2

wCO2 wCHO

wCH4 4.0E-2

2.0E-6

2.0E-2 0.0E0

1.0E-6 0.0E0

0.01

wH2O

0.06

0.01

0.11

0.01 0.06

0.01

0.11

0.05

wH2O

0.05

0.09 0.09 0.13 0.13

wCO2 wCO2 wCH2O

T (K)

1.2E-5 8.0E-6 4.0E-6 0.0E0 0.01

0.05

0.09

0.01

2.0

E3

1.5

E3

1.0

E3

5.0

E2 0.01

0.06 0.11

wH2O

0.1

0.0

5 0.0

9 0.1

wCO2

wH2O

1 0.0

0.13

1

6 0.0

3

wCO2

Abbildung 3.15: Zweidimensionale Mannigfaltigkeit im Reaktionsraum fur das stochiometrische CH4 -O2 -N2 -System. Zusatzliche Linien sind Trajektorien,  kennzeichnet den chemischen Gleichgewichtspunkt.

3.3. INTRINSISCHE NIEDRIG-DIMENSIONALE MANNIGFALTIGKEITEN λm (1/s)

49

λm+1 (1/s) 0.0E0

0.0E0

5 -1.0E

5 -1.0E -2.0E

5

-3.0E

5

5 -2.0E 5 -3.0E 5 -4.0E

5 -4.0E -5.0E

5

-6.0E

5

wH2O

0.01

0.02

0.07

0.12

5 -5.0E

wH2O

5 -6.0E 5 -7.0E 0.01

0.02

0.06

0.10

0.05

0.05

0.09

0.09

0.13

0.13

wCO2

wCO2

Abbildung 3.16: Eigenwerte des stochiometrischen Methan-Luft Gemisches einer zweidimensionalen Mannigfaltigkeit uber die Reaktionsfortschrittsvariablen CO2 und H2 O . Das rechte Pro l entspricht dabei dem langsamsten entkoppelten Eigenwert, das linke dem schnellsten der Prozesse, die die Bewegung auf der Mannigfaltigkeit beschreiben.

3.3.3 Numerische Berechnung niedrig-dimensionaler Mannigfaltigkeiten Einfuhrung

In den letzten Kapiteln wurde gezeigt, da es zu einem chemischen Reaktionssystem

@ ~ = F~  ~  @t

niedrig-dimensionale Unterraume gibt, in die die chemische Kinetik hineinrelaxiert. Die niedrigdimensionale Mannigfaltigkeiten werden durch die Gleichgewichtsbedingungen (3.18)

   

QTf ~ F~ ~ = 0 de niert. Eine direkte Verwendung solcher Gleichgewichtsbedingungen in den Programmen der chemisch reaktiven Stromungssimulation ist aber infolge der komplexen Struktur chemischer Reaktionskinetik und der damit verbundenen aufwendigen Berechnung der Eigenwerte und Eigenvektorbasen nicht ratsam. Deshalb werden die intrinsischen niedrig-dimensionalen Mannigfaltigkeiten vorab berechnet und in Form tabellarischer Bibliotheken den jeweiligen Programmen zur Verfugung gestellt. (3.18) ist eine implizite Darstellung der Mannigfaltigkeit (Hyper ache im Zustandsraum). Fur eine vollstandige Beschreibung mu das Gleichungssystem noch vervollstandigt werden, wie im Anschlu beschrieben wird. Danach werden Probleme aufgezeigt, die bei einer numerischen Losung des Gleichungssystem auftreten konnen und deren Konsequenzen.

Parametergleichungen

Um die Gleichgewichtsbedingungen (3.18), die die m-dimensionale Mannigfaltigkeit im n-dimensionalen Zustandsraum beschreiben, vollstandig zu losen, mussen m = n , nf zusatzliche Parametergleichungen eingefuhrt werden. nf ist dabei die Anzahl der schnellsten Zeitskalen, die vom chemischen System entkoppelt werden sollen. Wie die Parametrisierung ausgefuhrt wird, hat fur

50

KAPITEL 3. AUTOMATISCHE REDUKTION VON REAKTIONSMECHANISMEN

die Mannigfaltigkeit selbst keinen Ein u. Die Parametergleichungen sollen mit





C~ ~ ;~ = 0 beschrieben sein, wobei ~ = (1 ; 2; : : :; m ) die ausgesuchten Parameter sind. Die Parametergleichungen konnen in zwei Klassen eingeteilt werden. Die erste beschreibt direkt die Erhaltungsgroen, wie z. B. (siehe vorheriges Kapitel) die spezi sche Enthalpie h, den Druck p und die Elementmolzahlen i fur geschlossene, homogene, isobare, adiabatische Systeme. Die (ne +2) Parametergleichungen lauten in diesem Fall 0 = h , h 0 = p , p   0 = i ~ , i

i = 1; : : :; ne

mit den Parametern h ; p; i . Naturlich konnen auch analoge Gleichungen fur isochore oder isotherme Systeme aufgestellt werden. Insgesamt hat man also das Gleichungssystem





0 = C~ 1 ~ ;~E   0 = C~ 2 ~ ;~C   0 = QTf F~ ~ ; zu losen mit ~E = (h ; p;  ; : : :), C~ 1 = Parametergleichungen der Erhaltungsgroen. C~ 2 bzw. ~C bezeichnen die noch zu bestimmenden Parameter bzw. ihre Gleichungen. Letztere sind frei wahlbar, wenn es gewahrleistet ist, da Existenz und Eindeutigkeit erfullt sind. Der Einfachheit halber nimmt man in der Regel spezi sche Molzahlen bestimmter Spezies, obwohl die Wahl der Parametrisierung keinen Ein u auf die Mannigfaltigkeit selbst hat. Wahlt man spezi sche Molzahlen bestimmter Spezies als zusatzliche Parameter, entspricht die Konstruktion der Mannigfaltikeit einer Projektion auf bestimmte Koordinaten des Zustandsraumes. Fur die Parameter erhalt man z.B. 1

0 = i , i

i = 1; : : :; nc :

Betrachten wir das fur das Methan-Luft System aus dem vorhergehenden Kapitel mit den 34 reaktiven Spezies. Das System wird beschrieben durch

n ne ns nr

Dimension des Zustandsraumes Anzahl der Elemente Anzahl der Spezies Dimension des Reaktionsraumes

= = = =

36 4 34

ns , ne = 30 :

Will man die Dynamik mit nur m = 2 Reaktionsfortschrittsvariablen beschreiben, so benotigt man nf = nr , m = 28 Gleichgewichtsbedingungen. Zur Parametrisierung des Gleichungssystems benotigt man ne + 2 = 6 Parametergleichungen fur die Erhaltungsgroen h; p; H; O ; C ; N und zwei weitere Gleichungen fur die reaktiven Freiheitsgrade. Fur diese sollen in diesem speziellen Fall (wie auch in den Beispielen, die weiter hinten vorgestellt werden) die spezi schen Molzahlen der Spezies CO2 und H2O gewahlt werden. Insgesamt gilt es also, das Gleichungssystem 0 = h , h 0 = p , p 0 = H , 

H

3.3. INTRINSISCHE NIEDRIG-DIMENSIONALE MANNIGFALTIGKEITEN

51

0 = O ,  0 = C ,  0 = N ,  0 = CO , CO 0 = H O , H O 0 = QTf F~ ~ zu losen. Sind Werte fur i bekannt, lassen sich Punkte auf der Mannigfaltigkeit im Zustandsraum ermitteln. Fur verschiedene Werte der dynamischen Freiheitsgrade (in diesem Fall fur CO2 und H2 O ) spannt man ein Gitter uber die m-dimensionale Mannigfaltigkeit. Diese kann dann in tabellarischer Form abgelegt werden. O

C

N

2

2

2

2

Losungsverfahren zur numerischen Identi kation Um das vorgestellte Gleichungssystem aufzulosen, bieten sich verschiedene Methoden an. Hier sollen die angewandte Methoden kurz vorgestellt werden.

Pseudo-Zeititeration:

1. Anfangswert ~0 (am besten eignet sich der Gleichgewichtswert, da dieser immer eine Losung der Mannigfaltigkeit darstellt), der die Parametergleichungen erfullt. 2. Fur i = 0; 1; 2; : : : 3. Berechne die Jacobimatrix F i von (3.3.3) zum Zeitpunkt i und deren reelle Schur-Zerlegung QTf F Q = Z der Jacobimatrix fur das gegeben ~ (i). 4. Lose das Gleichungssystem   0 = C~ ~ (i);~



  ~ (i) QTf @ @t = QTfF~ ~ (i) furt ! 1 :

5. Falls

~ (i) , ~ (i,1)

<  : Konvergenz, stop sonst: Gehe zu Punkt 3 . Als Solverpakete stehen das implizite Extrapolationsverfahren LIMEX [37, 38] und das BDF-Verfahren DASSL [51] zur Verfugung.

Newton-Iteration: Verwendet wurde ein modi ziertes Newton-Verfahren [52,53].

1. Anfangswert ~0 der die Parametergleichungen erfullt. 2. Fur i = 0; 1; 2; : : : 3. Berechne die Jacobimatrix F i von (3.3.3) zum Zeitpunkt i und die reelle Schur-Zerlegung QTf F Q = Z der Jacobimatrix fur den Zustand ~ (i). 4. Lose das System     QTf F ~ (i+1) , ~ (i) = ,F~ ~ (i);~     C i ~ (i+1) , ~ (i) = ,C~ ~ (i);~ : C i = Jacobimatrix der Parametergleichungen. ( )

( )

52

KAPITEL 3. AUTOMATISCHE REDUKTION VON REAKTIONSMECHANISMEN



5. Falls

~ (i) , ~ (i,1)

<  : Konvergenz, stop sonst: Gehe zu Punkt 3 .

Fortsetzungsverfahren: Hier geht man von einem Punkt ~0 (~0) aus, von dem man wei, da

er auf der Mannigfaltigkeit liegt. Dies ist in jedem Fall der Gleichgewichtspunkt, den man immer bestimmen kann. Danach werden sukzessive weitere Punkte gesucht, die die Gleichungen ebenfalls erfullen. Es ist also der Zustand ~ 0(1) auf der Mannigfaltigkeit gesucht, der zu einem neuen Parametersatz ~1 gehort. Dies lat sich als eindimensionales Problem formulieren. Man verfolgt eine Linie auf der Mannigfaltigkeit ausgehend von einem bekannten Punkt. Das zu losende Gleichungssystem ist     0 = QTf ~ F~ ~   0 = C~ ~ ;~ ~ = ~0 + (~1 , ~0) mit als Fortsetzungsparameter. Dies wurde in den Paketen ALCON und PITMAN [54{57] realisiert. Diese Pakete stehen allerdings nur fur 2-dimensionale Mannigfaltigkeiten bereit. Ein groer Vorteil dieser Methoden ist, da man hier Informationen uber kritische Punkte, wie Wendepunkte oder Bifurkationen, sprich uber Existenz und Eindeutigkeit, erhalt.

Methoden zur Berechnung von Eigenvektorbasen: Bei der Berechnung der Mannigfaltig-

keiten mussen standig Zerlegungen der Jacobimatrix vorgenommen werden. Benutzt man die Basis der Schur-Vektoren, so mu die Zerlegung

0 z11 z12 B 0 z22 QT F Q = Z = B B@ ...

: : : z1n 1 : : : z2n C C

. . . ... C A 0 : : : 0 znn bestimmt werden. Q ist darin die orthogonale Matrix der Schur-Vektoren ~qi . In der Diagonalen der oberen Dreiecksmatrix Z stehen im Falle reeller Eigenwerte die Eigenwerte i selbst. Fur den Fall komplex konjugierter Eigenwerte steht an dieser Stelle eine (2  2)-Matrix, bestehend aus Realund Imaginarteil der Eigenwerte. Auerdem sollen die Eigenwerte nach absteigendem Realteil in der Matrix Z angeordnet sein (zii  zi+1;i+1 ). Um die Schur-Zerlegung zu bestimmen, wird eine modi zierte QR-Zerlegung durchgefuhrt [48]. Danach mussen die Schur-Vektoren sortiert werden.

Tabellierung der Mannigfaltigkeit Es wurde gerade beschrieben, wie ein Punkt auf der Mannigfaltigkeit numerisch bestimmt werden kann. Vollzieht man die Prozedur an verschiedenen Punkten der reaktiven Variablen und legt so ein Gitter uber die Mannigfaltigkeit, kann diese in tabellarischer Form abgelegt werden. Sie kann dann als Bibliothek den verschiedensten Programmen zur Simulation reaktiver Stromungen zur Verfugung gestellt werden kann.

Bestimmung des Tabellierungsbereiches: Im Voraus mu geklart werden, in welchem Bereich

physikalisch realisierbare Parameter des Zustandsraumes liegen. Die Massenbruche der Spezies sind nach oben und unten beschrankt: 0  wi  1. Ebenfalls beschrankt sind die physikalischen Parameter fur Temperatur T > 0, Druck p > 0 und Dichte  > 0. Fur die oberen Grenzen der physikalischen Parameter mussen Extremwerte bestimmt werden, die die Parametergleichungen erfullen. Im einzelnen bedeutet das Einschrankungen fur die physikalisch-chemische Variablen.

3.3. INTRINSISCHE NIEDRIG-DIMENSIONALE MANNIGFALTIGKEITEN

53

Erlaubter Bereich der Elementmassenbruche: 0   i i = 1; : : :; ne P n s 1 = j =1 Mje xj : wobei Mje = molare Massen der Elemente j . Ferner sind gegebenenfalls zusatzliche Beschrankungen durch das betrachtete System zu erfullen. Z. B. gilt fur homogene Reaktoren oder Stromungen mit vernachlassigbarem molekularem Transport, da die Elementzusammensetzung konstant ist, da diese sich durch chemische Reaktionen nicht andern kann. Dasselbe gilt fur vorgemischte Flammen unter der Annahme vernachlassigbarer dierentieller Diusion.

Erlaubte Bereiche von Enthalpie, Druck und Dichte:  ~  ~  ~  T

;p

;

>0

ist eine absolute Begrenzung. Oft ist bei Flammenbetrachtungen fur kleine Stromungsgeschwindigkeiten (v