Tudományos Diákköri Dolgozat
SZABÓ BOTOND
Arrhenius-paraméterek becslése közvetett és közvetlen mérések alapján
Turányi Tamás. Zsély István Gyula Kémiai Intézet
Eötvös Loránd Tudományegyetem Természettudományi Kar Budapest, 2009
Tartalomjegyzék 1. Bevezetés................................................................................................................................ 3 2. Irodalmi áttekintés.................................................................................................................. 5 2.1. Reakciókinetikai paraméterek meghatározása mérési adatok alapján ............................ 5 2.2. Az Arrhenius-paraméterek bizonytalansága ................................................................... 6 2.3. A Monte Carlo latin hiperkocka mintavételezési eljárás................................................. 8 2.4. A rangkorreláció.............................................................................................................. 9 2.5. A felhasznált programok ............................................................................................... 10 3. Az Arrhenius-paraméterek eloszlásának meghatározására .................................................. 12 3.1. Az algoritmus ................................................................................................................ 12 3.2. A modellrendszer kiválasztása ...................................................................................... 14 3.3. Eredmények a H+O2→O+OH és a H+O2+M→HO2+M reakcióra............................... 15 3.4. Továbbfejlesztési lehetőségek....................................................................................... 18 4. Összefoglalás........................................................................................................................ 20 Irodalom ................................................................................................................................... 22
2
1. Bevezetés A reakciókinetikában gyakran használnak nagy reakciómechanizmusokat [1]. Egy közepes méretű mechanizmus esetén mintegy 50 anyagfajta van jelen és a reakciók száma nagyjából 500. Kémiai folyamatok vagy kísérletek szimulációja közönséges vagy parciális differenciálegyenletek megoldása segítségével történhet. A reakciómechanizmusokban szereplő paraméterek kapcsolódhatnak anyagfajtákhoz (például képződési entalpia, entrópia, hőkapacitás, diffúziós együttható), és reakciólépésekhez (például A, n, E Arrheniusparaméterek, esetenként a reakció nyomásfüggését leíró paraméterek). A szimulációk során fel kell használni a rendszer fizikai modelljét leíró paramétereket is (például hőmérséklet, nyomás, tartózkodási idő). A reakciómechanizmusok paramétereit sokféle módszerrel határozhatják meg. A dolgozatban leírt számítások során egyes reakciólépések Arrheniusparamétereit fogjuk becsülni. A k sebességi együttható hőmérsékletfüggését legáltalánosabban a kiterjesztett Arrhenius-egyenlet adja meg: k (T ) = AT n exp(− E / RT ) ,
(1.1)
ahol A, n, E az Arrhenius-paraméterek, R a gázállandó (8,314 J K-1 mol-1) és T a hőmérséklet (K). A sebességi együtthatók kísérleti vagy elméleti úton is meghatározhatók. Egy-egy reakció sebességi együtthatójának meghatározására gyakran több független mérést illetve számítást végeznek, és ezek eredménye általában többé-kevésbé összhangban van egymással [2]. Ugyanakkor az egyes sebességi együtthatók még mindig jelentősen bizonytalanok [3] és ez a szimulációk során számított eredményeket jelentősen bizonytalanná teszi [4-9]. Turányi és munkatársai megmutatták [4-9], hogy a sebességi együtthatók bizonytalansága okozza legnagyobb részt a reakciókinetikai szimulációk eredményének bizonytalanságát. Az egyes reakciólépések sebességi együtthatójának meghatározására irányuló kísérleteket nevezhetjük közvetlen méréseknek. A reakciómechanizmusokat reakciólépésekből állítják össze, majd úgynevezett közvetett kísérletekkel szemben tesztelik. Ezeknél a kísérleteknél olyanok a körülmények, hogy megfelelő szimulációs eredmény eléréséhez több reakció paraméterének kell megfelelőnek lennie, így ezek a kísérletek a teljes mechanizmus ellenőrzésére alkalmasak. Ilyen közvetett kísérletek például a gyulladási idő mérések, a lángsebesség mérések, vagy a koncentrációprofilok meghatározása csőreaktorban, lángokban. Természetesen a vizsgált kísérleti körülményektől függ, hogy milyen anyagok és mely reakciók paramétereinek pontos ismerete szükséges a kísérletek pontos szimulációjához. Az égéskémia területén az egyik legegyszerűbb rendszer a hidrogén−levegő elegy égése, ezért ezt választottuk modellrendszernek. A hidrogén−levegő elegy égését leíró, a közelmúltban publikált reakciómechanizmusok [10-12] sok mérést elfogadható pontossággal leírnak, de egyik sem képes az összes mérési eredményt hibahatáron belül reprodukálni. Meglepő módon számos reakció esetén jelentősen eltérő értékeket használnak a sebességi együtthatókra. Ez nem egyedi jelenség, különböző kutatócsoportok sokszor jelentősen eltérő rekaciómechanizmusokat közöltek ugyanarra a kinetikai rendszerre [13]. A fentiek alapján látható, hogy a mai napig hiányzik egy olyan reakciómechanizmus, amely akár a hidrogén viszonylag egyszerű égését kellő pontossággal képes lenne leírni minden kísérletet figyelembe véve. A dolgozat célja, hogy egy egyszerű, néhány reakcióból (így kevés paraméterből) álló modellrendszeren kifejlesszünk egy olyan módszert, ami képes az összes (közvetett és közvetlen) mérés figyelembe vételével egyes reakciólépések Arrhenius-paramétereinek meghatározására. Célunk továbbá a kapott paraméterek kovariancia struktúrájának,
3
eloszlásának meghatározása. A paraméterek kovariancia bizonytalanságának megbecslésében lehet segítségünkre.
4
mátrixa
a
szimulációk
2. Irodalmi áttekintés 2.1. Reakciókinetikai paraméterek meghatározása mérési adatok alapján A nagy, részletes reakciómechanizmusok fejlesztésének alapvető módszere, hogy felhasználják az egyes reakciókinetikai paraméterek közvetlen mérésekkel meghatározott értékét. A közvetlenül mért sebességi együtthatók meghatározásának jellemző hibája viszonylag magas, például 30%-os teljes hiba (szisztematikus és statisztikai hiba együtt) jellemzőnek tekinthető. A közvetlen mérések során kapott reakciókinetikai paramétereket a mechanizmus készítésénél általában már nem változtatják meg, és a reaktorokban kapott közvetett mérési adatokat (pl. gyulladási idő, koncentráció−idő görbék) csak ellenőrzésre veszik figyelembe. Egy másfajta mechanizmusfejlesztési módszer kidolgozását a Gas Research Institute (GRI) támogatta. A GRI megbízása alapján Michael Frenklach és munkatársai egy olyan módszert dolgoztak ki [14-16], mely szerint nem csak a sebességi együttható (k) közvetlen méréseken alapuló meghatározását veszi figyelembe a reakciómechanizmus paramétereinek megválasztásánál, hanem a közvetett mérési adatokat (pl. gyulladási idő, koncentráció−idő görbék) is. Modellrendszernek a metán égését választották. Összegyűjtöttek csaknem száz megbízhatónak tartott kísérletet, amelyek a metán égését különböző szempontokból vizsgálták. Ezen kísérletek eredményének pontos reprodukciója volt az optimalizációs eljárásuk célja. Érzékenységanalízissel kiválasztották azokat a reakciókat, amelyek jelentősen befolyásolták a eredményeket. Ezen reakciók sebességi együtthatójához közvetlen mérések alapján meghatároztak egy-egy [kmin, kmax] intervallumot úgy, hogy a sebességi együttható várható értéke bizonyosan beleessen ezekbe az intervallumokba. Ezenkívül más megkötést nem tettek fel k-ra a közvetlen mérések alapján. A paraméterbecslés során ezen reakciók preexponenciális paraméterét (A Arrhenius-paraméter) illesztették a többi paraméter rögzített értéke mellett. Az így kapott reakciómechanizmus az úgynevezett GRI-mechanizmus [15], mely továbbfejlesztett változatai [14] és [16] a mai napig nagy népszerűségnek örvendenek az égéskémikusok körében. A módszer egyik hiányossága, hogy a közvetlen mérési adatokat csak a tartomány meghatározására használták, így nem ritkán előfordult, hogy az optimalizált érték az elfogadási intervallum határára esett. Ez arra utal, hogy a paraméter értéke nem optimális, vagy az illesztési intervallum volt helytelenül megválasztva. Kémiai szempontból a legnagyobb probléma, hogy az A Arrhenius-paraméter illesztése rögzített n és E Arrheniusparaméterek esetén nem értelmes. Ekkor ugyanis a sebességi együttható értékét csak eltoljuk, annak hőmérsékletfüggését nem változtatjuk meg. Az összes Arrhenius-paraméter illesztése sokkal célravezetőbb, kémiailag értelmesebb megoldás lenne.
5
2.2. Az Arrhenius-paraméterek bizonytalansága Ebben az alfejezetben Nagy Tibor PhD értekezésének [17] a dolgozat témájával egybevágó részét foglaltuk össze. Induljunk ki az (1.1) egyenletben megadott kiterjesztett Arrheniusegyenletből és vegyük annak a logaritmusát: ln k (T ) = ln A + n ⋅ ln T − E R ⋅ T −1 .
(2.1)
Vezessük be az α = ln A , κ = ln k és ε = E/R jelöléseket. Az α, n, és ε paraméterekre a (2.2) egyenlet lineáris:
κ (T ) = α + n ⋅ ln T − ε ⋅ T −1 .
(2.2)
A κ(T) hőmérsékletfüggő valószínűségi változó, így a (2.2) egyenletből meghatározható α, n, ε származtatott Arrhenius-paméterek is valószínűségi változók lesznek. Miután az Arrhenius-paraméterek és így transzformáltjaik is egy adott hőmérséklet tartományon fizikai állandók, így az α, n, ε valószínűségi változók együttes sűrűségfüggvénye hőmérsékletfüggetlen. Ebből következően várható értékeik ( α , n , ε ), szórásnégyzeteik ( σ α2 , σ n2 , σ ε2 ) és korrelációik ( rαn , rαε , rεn ) is hőmérsékletfüggetlenek. A korrelációs együtthatókra és szórásokra teljesül, hogy
0 ≤ σα ,σ n ,σ ε (2.3)
− 1 ≤ rαn , rαε , rnε ≤ +1 Vezessük be a p = (α, n, ε ) jelölést a belőlük képzett vektorra, ekkor a Σ p kovariancia mátrixuk a következőképpen határozható meg:
T
Σ p = (p − p )(p − p )
σ α2 = rαn σ α σ n rαε σ α σ ε
rαn σ α σ n
σ n2 rnε σ n σ ε
rαε σ α σ ε rnε σ n σ ε . σ ε2
(2.4)
A kovariancia mátrix pozitív szemi-definit, így a korrelációs együtthatókra teljesül, hogy 0 ≤ 1 − rα2n − rαε2 − rn2ε + 2 rαn rαε rnε .
(2.5)
Tegyük fel, hogy a κ(T) hőmérsékletfüggő valószínűségi változó eloszlása ismert a [T1 ,T2 ] hőmérsékletintervallumon. Adott T∈ [T1 ,T2 ] hőmérsékleten jelöljük κ (T ) -vel a várható értékét és σ κ2 (T ) -vel a szórásnégyzetét. Így a (2.2)-es egyenletből az alábbi összefüggés teljesül a sebességi együttható logaritmusának adott hőmérsékleten vett várható értékére és a származtatott Arrhenius-paraméterek várható értékére: κ (T ) = α + n ⋅ ln T − ε ⋅ T −1 .
6
(2.6)
A (2.2) és (2.6) egyenletekből az alábbi kapcsolat írható le sebességi együttható és a származtatott Arrhenius-paraméterek kovarianciái között: 2
σ κ2 (T ) = (κ (T ) − κ (T ))2 = ((α + n ⋅ ln T − ε ⋅ T −1 ) − (α + n ⋅ ln T − ε ⋅ T −1 )) = = σ α2 + σ ε2T −2 + σ n2 ln 2 T − 2rαε σ ασ ε T −1 − 2rεnσ ε σ nT −1 ln T + 2rαnσ ασ n ln T
(2.7)
A (2.7) képletből látható, hogy a sebességi együttható szórásnégyzete pontosan akkor hőmérsékletfüggetlen, ha σ α2 -n kívül a jobb oldalon minden paraméter értéke nulla, azaz egyedül az α paraméter bizonytalan és a többi értéke pontosan ismert. A következőkben megmutatjuk, hogyan lehet az Arrhenius-paraméterek eloszlását meghatározni a különböző hőmérsékleten vett sebességi együtthatók eloszlásából. Vezessük be a κi = κ(Ti) jelölést. Így az alábbi lineáris összefüggés írható fel a különböző hőmérsékleten vett sebességi együtthatók és az transzformált Arrhenius-paraméterek között:
1 ln T1 1 ln T2 1 ln T3
− T1−1 α κ1 − T2−1 n = κ 2 − T3−1 ε κ 3
(2.8)
Vezessük be továbbá az alábbi jelöléseket: 1 ln T1 Θ = 1 ln T2 1 ln T3
κ 1 κ = κ 2 . κ 3
α p = n , ε
− T1−1 − T2−1 , − T3−1
(2.9)
Ezekkel a jelölésekkel a (2.8) egyenlet Θp = κ alakba írható fel. Miután a Θ mátrix nem szinguláris, az átparaméterezett eloszlások sűrűségfüggvényeinek kapcsolatából a következő összefüggés írható fel az fκi(κ) sebességi együttható logaritmusának Ti hőmérsékleten vett sűrűségfüggvénye és az fp (p) átparaméterezett Arrhenius-paraméterek sűrűségfüggvény között: f κ (κ )d 3κ = f κ (Θp ) det Θ d 3 p ⇒ f p (p ) = f κ (Θp ) det Θ .
(2.10)
A kapcsolat természetesen a másik irányba is meghatározható:
(
)
(
)
(
)
f p (p )d 3 p = f p Θ −1κ d 3 p = f p Θ −1κ det Θ −1 d 3κ ⇒ f κ (κ ) = f p Θ −1κ det Θ −1 .
(2.11)
A csonkolt normális eloszlás természetes feltevés adott hőmérsékleten a sebességi együttható logaritmusának eloszlására. Általában m σ κ (T ) -nél szokták levágni a normális eloszlást, ahol m=2 vagy 3. Tehát a T∈ [T1 ,T2 ] tartományban a κ(T) változó sűrűségfüggvénye a [ κ (T ) -m σ κ (T ) , κ (T ) + m σ κ (T ) ] intervallumon g1 (κ ; T ) =
κ (T ) − κ (T ) 2 , exp− 2σ κ2 (T ) 2π σ κ (T )
(
c1
7
)
(2.12)
ahol a c1>1 állandó a levágás miatt szereplő szorzótényező. A fent megadott intervallumon kívül nulla a sűrűségfüggvény értéke. A κ(T) és p homogén lineáris kapcsolatából következően a származtatott Arrheniusparaméterek ( p = (α, n, ε ) ) eloszlására is feltehető, hogy normális, ha a κ(T) eloszlása normális. Ennek oka, hogy normális eloszlások lineáris kombinációja is normális. Ezek alapján a származtatott Arrhenius-paraméterek levágás nélküli sűrűségfüggvénye: f (p ) =
1 N /2
(2π )
1 T exp − (p − p ) Σp−1 (p − p ) . detΣp 2
(2.13)
A mi esetünkben azonban κ(T) csonkolt normális eloszlású volt, de miután a csonkolás m σ κ (T ) -nél történik, ahol m ≥ 2 , amely miatt κ(T) értékeknek legalább 95%-a megbízhatósági intervallumba esik, így az eloszlása jó közelítéssel normális, ezért a csonkolás csak kis mértékben változtatja meg a származtatott Arrhenius-paraméterek kovarianciáját és normáltságát. Megállapítható, hogy az (α, n, ε ) paraméterek normálisából a megbízhatósági tartományon kívül eső κ(T) értékeket levágva származtatott eloszlás jó közelítéssel megadja a κ(T)-re feltételezett csonkolt normális eloszlást. Az égéstudományi adatbázisokban szereplő összefüggések a sebességi együttható logaritmusának szórása és bizonytalanságok között a következő: f (T ) =
3 σ κ (T ) . ln 10
(2.14)
Vezessük be az egyenlet jobb oldalára az Mσκ(T) jelölést. A bizonytalanságot legalább hat hőmérsékleten ismernünk kell, hogy a származtatott Arrhenius-paraméterek kovariancia mátrixát ki tudjuk számítani. A kiszámításra a következő egyenletet használhatjuk: f (T ) = M σ α2 + σ ε2T −2 + σ n2 ln 2 T − 2rαε σ α σ ε T −1 − 2rεnσ ε σ n T −1 ln T + 2rαnσ α σ n ln T .
(2.15)
Természetesen vigyáznunk kell még arra, hogy (2.3) és (2.5) teljesüljön a szórásokra és korrelációs együtthatókra.
2.3. A Monte Carlo latin hiperkocka mintavételezési eljárás Egy adott valószínűségi eloszlás alapján minta generálására rengeteg módszer lehetséges (például Monte Carlo eljárás [18] vagy Markov-lánc Monte Carlo módszer [19]. A kapott minta általában a nagy valószínűséggel bekövetkező értékekre fog koncentrálódni, azaz kevés minta esik majd a határokra. Számunkra szintén a nagyobb valószínűségű tartományok az érdekesek, de nem szeretnénk elveszteni az eloszlás kis valószínűségű, határon lévő értékeiből származó információt sem. Tegyük fel, hogy egy d dimenziós, független koordinátájú eloszlásból szeretnénk mintavételezni. Legyen az i-edik koordináta sűrűségfüggvénye fi és osszuk fel a sűrűségfüggvény tartóját n azonos valószínűségű intervallumra fi szerint. Így nd d-dimenziós téglatestet kapunk (nem feltétlenül azonos méretűeket), melyek valószínűsége megegyezik. A latin hiperkocka módszerben ezek közül véletlenszerűen választunk ki n darabot a következő eljárás szerint. Tekintsük minden koordináta esetén az n azonos valószínűségű intervallum egy-egy permutációját. Az egyes koordinátákra kapott permutációk alapján állítsuk össze a d-dimenziós téglatesteinket úgy, hogy az i-edik téglatestet az egyes koordinátákra kiválasztott intervallumpermutációk i-edik tagjai határozzák meg. Az így 8
kapott n téglatestből ezután az eredeti eloszlásnak a téglatestre megszorított változata szerint választunk ki egy-egy pontot. Ezt a mintavételezési eljárást nevezzük Monte Carlo latin hiperkocka módszernek [1]. Az így kapott adatsor az eloszlásból számolt mintának tekinthető, és megvan az az előnye, hogy kevésbé koncentrálódnak az adatok, mint hagyományos mintavételezés esetén. Az alábbi két ábrán szemléltetjük a különbséget a hagyományos és a latin hiperkocka módszerrel végzett mintavételezés között. A teljesen véletlen mintavételezésnél egyes tartományok üresen maradhatnak (nagy körök az 1. ábrán), más tartományokban pedig egymáshoz közel kerülnek a pontok (kis körök). Ilyen probléma kevésbé jelentkezik a latin hiperkocka mintavételnél (l. 2. ábra).
1. ábra
2. ábra
100 elemű minta véletlen mintavételezéssel
100 elemű minta latin hiperkocka mintavételezéssel
2.4. A rangkorreláció A rangkorrelációt két sorozat együttváltozásának mérésére szokás alkalmazni. A rangkorreláció mértékét a rangkorrelációs együttható írja le, mely −1 és +1 közötti értékeket vehet fel. Amennyiben az egyik sorozat ugyanazon helyeken vesz fel nagyobb illetve kisebb értékeket, mint a másik, akkor a rangkorrelációs együttható pozitív értéket vesz fel. Minél erősebb ez a kapcsolat, annál közelebb lesz +1-hez az együttható értéke. Fordított esetben, ha az egyik sorozat ott vesz fel kisebb értékeket, ahol a másik nagyobbakat és fordítva, akkor a rangkorrelációs együttható negatív lesz és minél erősebb ez a kapcsolat annál közelebb lesz −1-hez az együttható értéke. A rangkorrelációt széles körben alkalmazzák eloszlások összehasonlítására. Ilyenkor a sűrűségfüggvény értékeit szokták összehasonlítani véletlenszerűen választott pontokban. A rangkorrelációt a minta megfelelő megválasztásával annyiban lehet manipulálni, hogy az eloszlások általunk fontosnak tartott tartományaira helyezzen nagyobb hangsúlyt, például az egyik eloszlás nagyobb valószínűségű tartományaira. A legnépszerűbb rangkorrelációs együtthatók a Spearman-féle ρ [20] és Kendall-féle τ [21]. Induljunk ki a sűrűségfüggvények az xi (i=1,…,n) pontokban felvett fi és gi értékeiből és jelöljük yi-vel és zi–vel a függvényértékek növekvő sorozatában betöltött sorszámukat. Legyen továbbá di az yi- zi sorszámok különbsége. Ekkor a Spearman-féle ρ rangkorrelációs együttható értéke a
9
ρ = 1−
6∑ d i n( n 2 − 1)
(2.16)
képlettel számolható ki. A Kendall-féle τ képletének megadásához először képezzük az összes lehetséges (fi ,gi) és (fj ,gj) párt, ahol i<>j. Egy ilyen párt megfelelőnek nevezünk, ha sgn( f i − g i ) = sgn( f j − g j ) ,
(2.17)
sgn( f i − g i ) = − sgn( f j − g j ) ,
(2.18)
és eltérőnek nevezünk, ha
ahol a sgn x függvény negatív értékeknél −1-t, pozitív értékeknél +1-t és nullában nullát vesz fel. Jelöljük ezután nc-vel a megfelelő, és nd-vel az eltérő párok számát a sorozatban. Ekkor a Kendall-féle τ képlete az alábbi lesz: τ =2
nc − nd . n( n − 1)
(2.19)
A képletben szereplő n(n-1)/2-es osztó az összes lehetséges párosítások számát adja meg. A számlálóban a megfelelő és eltérő párok számának különbsége szerepel, így ha a párok többsége megfelelő, akkor a hányados közel lesz egyhez, míg ha a párok többsége eltérő, akkor a τ értéke −1-hez fog közelíteni. A dolgozat során mi végül a Spearman-féle ρ rangkorrelációs együtthatóval dolgoztunk, de a számításokat a Kendall-féle τ-val elvégezve is hasonló eredményt kaptunk.
2.5. A felhasznált programok A MATLAB programcsomagot a FORTRAN nyelvből fejlesztették ki eredetileg azzal a céllal, hogy a tömbműveletek indexelési hibáiból származó hosszú hibakeresést minimalizálják. Az elmúlt évtizedek során a program nagyon sokoldalúvá vált, és a mérnöki területen ismerete alapvető. Tudományos kutatás során is kitűnően használható, mivel magasszintű programnyelvként nagyon könnyű megtanulni a használatát, ugyanakkor a mátrixműveletek igen egyszerűen programozhatók vele, és a programfejlesztés nagyon gyors benne. Meglévő FORTRAN vagy C nyelvű programokat is egyszerűen integrálni lehet MATLAB programokba. Ezen előnyös tulajdonságai miatt választottuk ezt a környezetet az algoritmus fejlesztésére. Az általunk kidolgozott algoritmus helyességének ellenőrzésére felhasználtunk egy szélsőérték kereső eljárást is. Az irodalomban fellelhető sokféle szélsőérték kereső módszer közül mi a MATLAB programcsomagba beépített fminsearch eljárást választottuk ki. Az fminsearch szubrutin szimplex kereső módszert használ [22]. Az n-dimenziós tér szimplexei azok a halmazok, melyeket n+1 adott térbeli pont feszít ki. Például a kétdimenziós tér szimplexei a háromszögek, míg háromdimenzióban a tetraéderek alkotják a szimplexek halmazát. Az fminsearch eljárás során kiindulunk egy szimplexből és minden egyes lépésben egy új pontot generálunk a szimplex közelében, vagy magán a szimplexen belül és kiszámítjuk a függvény értékét az így kapott pontban. A függvényértéket összehasonlítjuk a szimplexet kifeszítő pontokban számolt értékekkel, és amennyiben kisebb értéket ad valamelyik feszítő pontban számolt értéknél, úgy az új pontot kicseréljük a generált pontra. A gyakorlat azt mutatja, hogy a lépések többségében történik pontcsere. Az algoritmus addig tart, míg a szimplex mérete megfelelően kicsi lesz. Az eljárás tartalmaz egy maximális 10
iteráció számot megadó paramétert, így kisebb futásidő alatt is kaphatunk közelítően jó eredményt. Az eljárás egy fontos tulajdonsága, hogy az esetek többségében nem okoz neki gondot, ha a függvényünk nem folytonos (főleg ha a megoldás egy környezetében az). A mi esetünkben az eloszlás határán lévő pontok környezetében a kikötések és levágások miatt a függvényünk nem feltétlenül lesz folytonos, de mint azt korábban láttuk, ez nem befolyásolja az algoritmust [23]. A latin hiperkocka mintavételezéses Monte Carlo módszerrel készült mintát Zádor Judit részben módosított FORTRAN programjával [8] állítottuk elő. A program bemenő paraméterei a változók száma, a generálandó minták száma, a paraméterek várható értéke, szórása, valamint a vágási határok voltak. A program független paraméterek mintavételezésére szolgál. Az általunk vizsgált paraméterek részlegesen korreláltak, így ezt a módszert közvetlenül nem tudtuk alkalmazni. Azonban többdimenziós korrelált normális eloszlás lineáris transzformációjával nem korrelált normális eloszlást tudunk kapni. Legyen ugyanis X~N(µ,Σ) d-dimenziós normális eloszlású valószínűségi változó, ahol a µ ddimenziós vektor a várható értéket és a Σ dxd-es szimmetrikus pozitív definit mátrix a kovariancia mátrixot jelöli. Legyen továbbá A qxd-es mátrix, ekkor Y=AX~N(Aµ,AΣAT) qdimenziós normális eloszlású valószínűségi változó [24]. Válasszuk meg az A mátrixot úgy, hogy az AΣAT mátrix diagonális legyen. Miután a Σ szimmetrikus mátrix, a normált sajátvektoraiból alkotott A unitér mátrixra AΣAT diagonális lesz. Az AΣAT transzformációval kapott mátrix főátlójában a Σ mátrix sajátértékei szerepelnek [25]. Tehát megfelelő lineáris transzformáció segítségével független koordinátájú többdimenziós normális eloszlás nyerhető (Y), melyre már alkalmazható a mintagenerálásra használt algoritmus. A kapott mintát pedig visszatranszformálva az eredeti normális eloszlás szerinti mintát kapunk. A visszatranszformálásra az AT mátrixot használhatjuk, mert definíció szerint az A mátrix unitér volt, így ATA=I (ahol I az egységmátrixot jelöli) és ebből következően ATY= ATAX=X~N(µ,Σ) az eredeti többdimenziós normális eloszlást adja vissza. Az égéskémia területén az első átfogó szimulációs környezet a CHEMIN programcsomag volt. Ennek második verzióját, a CHEMKIN-II csomagot [26] a mai napig elterjedten használják a tudományos kutatásban. Elterjedtségének oka, hogy – ellentétben a CHEMKIN program későbbi változataival – a forráskódja rendelkezésre áll, így működése teljesen átlátható és a program szükség szerint módosítható. Miután a felhasználok is hozzáfértek a forráskódhoz, a program csaknem összes hibáját kijavították az évek folyamán. A folyamatos hibajavításoknak köszönhetőn a programcsomag nagyon megbízható [1]. A MATLAB képes FORTRAN programokat futtatni úgy, hogy egy saját fordító segítségével úgynevezett .mex fájllá alakítja és ennek tartalmát közvetlenül, a belső MATLAB függvényekhez hasonlóan hívja. A kinetikai egyenletek megoldását végző SENKIN programot ilyen módszerrel használtuk. A bemenő paraméterek a hőmérséklet és a reakciókhoz tartozó paraméterek voltak. Ez a program eredményül a gyulladási időt szolgáltatta.
11
3. Az Arrhenius-paraméterek eloszlásának meghatározására 3.1. Az algoritmus Ebben az alfejezetben ismertetjük egy tetszőleges modell paramétereinek meghatározására kidolgozott módszerünket. Célunk egy olyan módszer kidolgozása volt, amely a közvetett és közvetlen mérési eredményeket egyaránt figyelembe veszi a paraméterek meghatározásánál. Ezt úgy szeretnénk elérni, hogy minimalizáljuk a paraméterek felhasználásával kapott szimulációs eredmények eltérését a közvetett mérési eredményektől, és a paraméterekkel számolt sebességi együtthatók eltérését a közvetlen méréssel kapott értékektől. A következőkben a módszer egy általános leírása következik, amelyet a 3.3. szakaszban alkalmazunk majd az általunk kiválasztott modellre. Tekintsünk egy olyan modellt, amelyben r reakció szerepel, és a közvetett kísérletekből s mérési adatsort ismerünk. Minden egyes reakciólépéshez tartozik legalább egy, a sebességi együtthatót közvetlenül meghatározó méréssorozat. Tegyük fel, hogy az i-edik közvetlen méréssorozat ni (i=1,..,r) és a j-edik közvetett méréssorozat mj (j=1,…,s) különbözőző hőmérsékleten történő mérést tartalmaz. Jelöljük p-vel a meghatározandó paraméterekből készített vektort. Legyen kim (Tji) az i-edik reakcióhoz tartozó Tji (j=1,….ni) hőmérsékleten kapott sebességi együttható mérési eredménye (az m betű utal arra, hogy mérési nem pedig számolt eredményről van szó). Továbbá jelöljük ki(Tji)-vel az i-edik reakcióhoz tartozó Tji (j=1,….ni) hőmérsékleten adott Arrhenius-paraméterekkel (adott esetben beleértve a nyomásfüggő reakciók paraméterezésére használt egyéb paramétereket is) számolt sebességi együttható értékét a p paramétervektor esetén. Vezessük be az yim, j és yi,j jelöléseket, melyek rendre az i-edik közvetlen kísérlet (j=1,….mi) mérési illetve szimulációs eredményeit jelentik, ahol a szimulációkor a p paraméter vektort használtuk. Első lépésben tekintsük az i-edik reakcióra a mért és a számolt ki értékek különbségét minden egyes Tji (j=1,…, ni) hőmérsékleten és osszuk le az így kapott különbség négyzetét a mérési eredmény négyzetével. Így egy normalizált, relatív eltérését kapjuk a mért és számolt értékeknek. Vegyük az így kapott tagok összegét és osszuk le a tagok számával, azaz ni-vel. Az ni-vel történő leosztás segítségével elérjük, hogy a különböző számú mérést tartalmazó méréssorozatok esetén kapott mérőszámok összehasonlíthatók legyenek. Ezután szorozzuk be az így kapott értéket egy ci konstans szorzóval, mely a kísérlet relatív megbízhatóságát fejezi ki. Az i-edik közvetlen mérésekre vonatkozó rész így: ci ni
ni
∑
(k (T ) − k i j
i
(T ji )
m i i 2 j
2
)
(3.1)
kim (T )
j =1
Végül az egyes reakciókra a fenti módon kapott értékek összegét véve állítsuk elő a sebességi együttható közvetlen mérésére vonatkozó célfüggvényt: r
ci ∑ i =1 ni
ni
∑ j =1
(k (T i
i j
) − k im (T ji )
k im (T ji ) 2
2
)
(3.2)
A közvetlen mérésekhez hasonló módon járjunk el a közvetett mérések esetében is. Első lépésben vegyük az i-edik közvetett kísérlethez tartozó mért és szimulált értékek különbségeinek négyzeteit és osszuk le az így kapott értékeket rendre a mérési eredmények négyzeteivel. Így a szimulált és mért értékek relatív négyzetes eltérését kaptuk. Ezután 12
vegyük normalizált négyzetes eltérések átlagát és szorozzuk be egy ki konstanssal, mely az egyes méréssorozatok relatív megbízhatóságát fejezi ki. Így az i-edik közvetett kísérletre kapott érték: ki mi
mi
∑
(y
i, j
− yim, j yim, j
j =1
2
)
2
.
(3.3)
Vegyük a közvetett kísérletekhez tartozó értékek összegét, mely a közvetett mérésekre vonatkozó célfüggvényt alkotja: s
ki ∑ i =1 mi
mi
∑
(y
i, j
j =1
− yim, j yim, j
2
2
)
.
(3.4)
A közvetlen és közvetett mérésekre vonatkozó együttes célfüggvény létrehozásakor lehetőségünk nyílik új súlyfaktorok bevezetésére (s1 és s2), melyek azt fejezik ki, hogy a közvetett illetve közvetlen mérési eredményeinket milyen arányban szeretnénk figyelembe venni. Ez a két súlyfaktor beépíthető a ci és ki konstansokba is. A következő összeg adja meg az algoritmusunk során felhasznált célfüggvényt: r
c f cél ( p) = s1 ∑ i i =1 ni
ni
∑ j =1
(k (T i
i j
) − k im (T ji )
k im (T ji ) 2
2
)
s
k + s2 ∑ i i =1 mi
mi
∑ j =1
(y
i, j
− yim, j yim, j
2
2
)
(3.5)
Célunk egy olyan paraméterkészlet megtalálása, mely a közvetett és közvetlen mérésekre vonatkozó együttes célfüggvény értékét minimalizálja. Továbbá szeretnénk egy olyan eloszlást hozzárendelni a p paramétervektorhoz, melyben nagy súlyt kapnak az alacsony célfüggvény értékű paramétervektorok, míg a nagy célfüggvény értékkel rendelkezőek kis valószínűségűek lesznek. A probléma megoldásához használhatunk relatívan gyors többdimenziós lokális szélsőérték kereső eljárást, ami könnyen rossz megoldást adhat, vagy alkalmazhatunk globális szélsőérték kereső eljárást, amely bonyolultabb célfüggvény esetén nagyon lassúak lehetnek. Továbbá szélsőérték kereső eljárással nem kapunk információt a paraméterek kovariancia struktúrájáról. A 2.2. alfejezetben leírtak alapján egy reakciólépés származtatott Arrheniusparamétereinek vektora többdimenziós normális eloszlású lesz. Természetesen adódik a feltételezés, hogy a kiválasztott reakciólépések Arrhenius-paramétereinek együttes eloszlása is normálishoz közeli lesz, ezért kiindulhatunk abból a feltételezésből, hogy a p paramétervektor normális eloszlású. Első lépésben mintát generálunk a (2.13)-ban megadott képlet segítségével Monte Carlo latin hiperkocka mintavételezéssel. Itt az egyes reakciólépéseket még függetlenek tekintjük egymástól. Minden reakciólépés esetében külön generálunk egy Arrhenius-paraméter vektort, melyeket majd összeillesztünk. A mintát érdemes megfelelően nagynak, több ezereleműnek választani a gyorsabb konvergencia érdekében. A következő lépésben súlyozni fogjuk a mintát a célfüggvény értékük alapján. Legyen fi=fcél(pi) a pi paramétervekorral számolt célfüggvény érték és ennek segítségével rendeljük a pi paramétervektorhoz a gi=exp(-γ fi) súlyt, ahol γ>0 egy optimalizálandó paraméter. Az ilyen alakú súlyozást az irodalomban Gibbs-mértéknek nevezik és elterjedten használják fizikai kísérletek leírásánál [27]. Az így kapott gi súlyok minden pi paramétervektorra pozitívak lesznek és kis célfüggvény érték esetén nagy, míg nagy célfüggvény érték esetén kis értéket vesznek fel. Ezután kiszámítjuk a súlyozott minta empirikus várható értékét és kovariancia mátrixát az alábbi képletek segítségével:
13
n
µˆ = ∑ g i i =1
n Σˆ = ∑ g i i =1
−1 n
∑g
−1 n
∑g (p i
i
pi ,
(3.6)
i =1
i
− µˆ )( pi − µˆ ) T .
(3.7)
i =1
Miután a gi súlyok a γ>0 paraméter függvényei voltak, így a µˆ empirikus várható érték és Σˆ empirikus kovariancia mátrix is a γ paraméter függvényei lesznek. Miután feltettük, hogy a paraméterek eloszlása többdimenziós normális, így a kapott várható érték és kovariancia mátrix egyértelműen meghatározza az eloszlást. Célunk azon γ paramétert meghatározni, melyhez tartozó normális eloszlás rangkorrelációban leginkább megközelíti a célfüggvény inverzét. Nagy rangkorreláció esetén teljesül, hogy a nagy célfüggvény értékű paramétervektorok kis súlyt, míg a kis célfüggvény értékű paramétervektorok nagy súlyt kapnak. Természetesen minket a nagyobb valószínűségű paraméterek viselkedése jobban érdekel, mint a kevésbé valószínűeké, de ez a módszerünkbe automatikusan bele van építve azzal, hogy generált mintával dolgozunk. Így a nagyobb sűrűségű paramétertartományokat több mintaelem képviseli, mint a kisebb sűrűségű paramétertartományokat. Amennyiben a módszer megkívánja egy további fi-től fordítottan arányosan függő súlyfaktort is bevezethetünk a rangkorrelációs képletünkbe, ezzel is nagyobb hangsúlyt fektetve a kis célfüggvény értékű paramétervektorok viselkedésének pontosabb leírására. A példában a 2.4 alfejezetben megadott Spearman-féle ρ-val számoltunk, de a Kendall-féle τ is alkalmazható lenne. A γ paraméter optimális értékének meghatározása után az így kapott normális eloszlásból újra generáljunk ezres nagyságrendű mintát. A korábbi mintánkat is érdemes megtartani a nagyobb információ mennyiség miatt, a súlyozás hatására a kevésbé valószínű paraméterek úgyse befolyásolják nagyon az eloszlás paraméterét. Az új mintaelemekre is számoljuk ki a célfüggvény értékeket és a γ paraméter függvényében számoljuk ki az empirikus várható érték és kovariancia mátrixot. Ezután γ-ban optimalizálva a maximális rangkorreláció szerint egy újabb többdimenziós normális eloszlást kaptunk. A módszert iteráljuk, amíg az empirikus paraméterek megfelelően nem konvergálnak.
3.2. A modellrendszer kiválasztása Az égéskémia területén az egyik legegyszerűbb rendszer a hidrogén−levegő elegy égése, ezért ezt választottuk modellrendszernek. A szénhidrogének égése során is kulcsfontosságú a hidrogén égési mechanizmusa, így annak pontos ismerete alapfeltétele a szénhidrogének reakcióinak kvantitatív leírásához. A közelmúltban számos kutatócsoport publikált hidrogén égési mechanizmust [10-12]. Ezek mindegyike számos kísérleten ellenőrzi a javasolt reakciómechanizmust, de egyik sem képes az összes mérési eredményt hibahatáron belül reprodukálni. A három idézett mechanizmus közül kiemelkedik a H. Curran által publikált, ugyanis Curran és kutatócsoportja az elmúlt években erre a mechanizmusra építve egészen 5 szénatomszám nagyságig közölt szánhidrogének égését leíró reakciómechanizmusokat, amelyek pontosan ugyanazt a hidrogén-részmechanizmust tartalmazzák. A dolgozatunkban felhasználtuk Sedyo Inez TDK dolgozatának eredményeit [28], melyben 45 kísérleti adatsor kiértékelésével számolta ki az egyes reakciólépésekhez tartozó Arrhenius-paraméter értékeket. Az egyes reakciólépésekhez tartozó f(T) értékeket több hőmérsékletre meghatározta 14
és a (2.7) és (2.14) képletek segítségével, figyelembe véve a (2.3) és (2.5) feltételeket, kiszámította a reakciólépésekhez tartozó kovariancia mátrixokat. A módszer fejlesztése során igyekeztünk olyan kis modellt választani, amely mérete folytán nem akadályozza az algoritmus fejlesztését, de már elég nagy ahhoz, hogy a paraméterillesztés eredményét kémiai szempontból ellenőrizni tudjuk. Olyan közvetett kísérleteket kerestünk, amelyeknél a szimulációs eredmények néhány reakcióra a többinél sokkal érzékenyebbek, így a paraméterek változtatása biztosan a számított értékek jelentős változását okozza. Ez alapján választottunk ki az irodalomból Petersen és munkatársai [29] valamint Slack és munkatársai [30] gyulladási idő méréseit, amelyeknél a számított gyulladási idő minden esetben a következő két reakcióra volt a legérzékenyebb: H+O2 → O+OH ,
(R1)
H+O2(+M) → HO2(+M) .
(R2)
Pontosabban az érzékeny paraméterek az (R1) reakció Arrhenius-paraméterei és az (R2) nyomásfüggő reakció alacsony nyomású határának Arrhenius-paraméterei. Sedyo Inez dolgozatában az (R1) reakcióra ajánlott Arrhenius-paraméterek: A1= 1,01*1015 dm3 s-1 mol-1, n1= −0,270, E1/R= 8063 K, míg az (R2) reakcióra ajánlott Arrhenius-paraméter értékek: A2= 1,73*1019 dm6 s-1 mol-2, n2= −1,23, E2/R= −39,6 K voltak. Az irodalomban a harmadik test paraméter értékét általában 0.67-nek veszik, de [28]-ből kiderül, hogy ez az érték a kísérleti adatok alapján inkább 0,47 és 0,51 közötti érték. A dolgozatban mi az m=0,5-ös értékkel számolunk. Sedyo Inez az alábbi korrelációs mátrixokat határozta meg az Arrhenius-paraméterekhez. Az 1. táblázat mutatja a p = (α, n, ε ) paramétereknek az (R1), a 2. táblázat az (R2) reakciókra vonatkozó kovariancia mátrixait.
8,644
-1,128
1171,47
10,71
-1,531
324,16
-1,128
0,1477
-155,2
-1,531
0,22
-45,89
324,16
-45,89
10000,0
1171,47
-155,2 169826,4 1. táblázat
(R1) reakcióhoz tartozó származtatott Arrhenius-paraméterek kovariancia mátrixa
2.táblázat (R2) reakciótartozó származtatott Arrhenius-paraméterek kovariancia mátrixa
3.3. Eredmények a H+O2→O+OH és a H+O2+M→HO2+M reakcióra Ebben a fejezetben a 3.1. szakaszban ismertetett algoritmusunk működését fogjuk az előző szakaszban megadott modellre bemutatni. Jelen esetben két reakciónk ((R1) és (R2)) és két közvetett mérési adatsorunk van (Petersen és Slack gyulladási idő mérései). A közvetlen mérési adatsorunkat a [29] és [30] cikkek megfelelő táblázataiból, míg a közvetett mérési adatsorunkat a [31] és [32] cikkek táblázataiból nyertük. Mindkét reakcióhoz három
15
Arrhenius-paraméter tartozik, így a feladatunk egy hatdimenziós eloszlás közelítése hatdimenziós normális eloszlással. Tegyük fel, hogy mind a közvetett, mind a közvetlen mérések ugyanolyan mértékben megbízhatóak, ezért a c1, c2, k1, k2, s1 és s2 súlyokat mind egységnyinek választottuk. A 3. táblázat tartalmazza az algoritmusunkkal kapott normális eloszlás várható értékét. H + O2 → O + OH reakció H + O2 (+M) → HO2 (+M) reakció ln (A/cm3mol-1s-1) n E/R ln (A/cm3mol-1s-1) n E/R 38,61 -0,771 8478 K 50,79 -2,15 -70,7 K 3. táblázat Ezentúl optimalizált paramétervektorként fogunk a többdimenziós normális eloszlás várható értékére hivatkozni, miután normális eloszlás esetén a várható érték egyben a legnagyobb sűrűségfüggvény értékkel rendelkező vektor az eloszlás tartójában. A kapott paraméterekkel számított gyulladási idő szimulációs eredményeket a 4.a és 4.b ábra, míg az Arrhenius-paraméterekből számolt sebességi együtthatókat a 3.a. és 3.b. ábra mutatja.
28
-1 -1
ln(k / cm mol s )
29
3
27 26 25 24 0.0005 0.0006 0.0007 0.0008 0.0009 0.0010 -1
T /K
-1
3.a ábra
3.b ábra
Az (R1) reakcióra vonatkozó sebességi együtthatójának ajánlott paraméterekre számolt értéke: — optimalizált paraméterekre számolt értéke: … közvetlen mérései: o [29], ∆ [30]
Az (R1) reakcióra vonatkozó sebességi együtthatójának ajánlott paraméterekre számolt értéke: — optimalizált paraméterekre számolt értéke: - közvetlen mérései (generált értékek): o
16
4.a ábra
4.b ábra
A Petersen-kísérlet gyulladási idejének ajánlott paraméterekre szimulált értéke: — optimalizált paraméterekre szimulált értéke: - közvetett mérési eredményei: o
A Slack-kísérlet gyulladási idejének ajánlott paraméterekre szimulált értéke: — optimalizált paraméterekre szimulált értéke: - közvetett mérései eredményei: o
Látható, hogy az optimalizált paraméterekkel számított gyulladási idő szimulációk eredménye lényegesen jobban megközelítik a reaktormérési eredményeket, mint az ajánlott paraméterekkel számított gyulladási idő szimulációk eredménye. Az (R1) reakcióra kapott közvetlen mérési eredményeink így jobban eltérhetnek az optimalizált paraméterekkel számított sebességi együtthatóktól, mint az ajánlott Arrhenius-paraméterekre kapottaktól, de ez a romlás az alkalmazott c1=c2=k1=k2=1 súlyok mellett kisebb, mint a közvetett mérések szimulálásánál és az (R2) reakciónál tapasztalt javulás. Az eljárás során kapott paramétervektor optimális voltát érdemes szélsőérték kereső eljárással leellenőrizni. Az optimalizált paramétervektorból kiindulva, a 2.5. alfejezetben leírt, fminsearch MatLab szélsőérték kereső eljárást futtattuk le és az így kapott legkisebb célfüggvény érték kevesebb, mint 1%-kal volt kisebb, mint az optimalizált paraméterekkel számolt célfüggvény érték. Így elmondható, hogy jó közelítéssel lokális minimumot talált az eljárásunk. Ahhoz, hogy a jó közelítéssel globális optimum tulajdonságát is belássuk az optimalizált paramétervektornak egy több pontból elindított szélsőérték kereső eljárást alkalmaztunk. Első lépésben Monte Carlo latin hiperkocka mintavételezéssel tízezres nagyságrendű mintát generáltunk, majd a legkisebb célfüggvény értéket adó paramétervektorokból elindítottuk az fminsearch szélsőérték kereső eljárást. Az így kapott legkisebb célfüggvény értékű vektort tekintettük globális minimum helynek, míg a célfüggvény értékét globális minimumnak. Ez ugyan csak egy közelítő eljárás, de az algoritmust többször megismételve hasonló nagyságú célfüggvény értéket adó vektorokat kaptunk eredményül, melyből levonható a következtetés, hogy nagy valószínűséggel a minimális célfüggvény értékhez közeli értékeket kaptunk az eljárás során. Az így előállított globális minimum érték elfogadható, 5%-nál kisebb mértékben tért el az optimalizált paraméterekkel számolt célfüggvény értéknél. Az optimális érték kiszámításán kívül érdemes a célfüggvényt magát is ábrázolni, hogy több információhoz jussunk a paraméterek viselkedéséről. Jelen esetben a célfüggvényünk egy hatdimenziós leképezés, ezért teljes ábrázolása lehetetlen, de kettesével tudjuk ábrázolni a paramétereket, ha a többi négy paramétert az optimalizált értéken lerögzítjük. Így bármely két paraméter feltételes korrelációjáról képet alkothatunk. Az ábrázolás során a paraméterekhez tartozó ajánlott érték középpontú és 3σ fél-sávszélességű intervallumokban vizsgáltuk a paraméterek értékét. Az intervallumon kívül eső paraméterek vizsgálatától eltekintünk, hiszen
17
a hozzájuk tartozó ajánlott valószínűségi érték elhanyagolható. Számunkra a kis célfüggvény értékkel rendelkező paramétervektorok az érdekesek, ezért érdemes a célfüggvény helyett valamilyen, vele fordítottan arányos függvényt ábrázolni. Miután a súlyokat exp(-γ fcél(pi)) alakban adtuk meg, így kézenfekvőnek tűnik az exp(- fcél(p)) függvényt ábrázolni.
6.a ábra
6.b ábra
α1- n 2 paraméterekre megszorított exp(- fcél(p))
α1- α 2 paraméterekre megszorított exp(- fcél(p))
Az ábrákon látható, hogy az exp(- fcél(p)) függvény kétdimenziós megszorítottjai egy-egy nyerget határoznak meg. A paraméterek közötti korreláció határozza meg a nyereg irányát az ábrákon. A 3.1. alfejezetben leírt módszerünk nem csak a normális eloszlás várható értékét, azaz egy optimalizált paramétervektort határoz meg, hanem a paraméterek egy kovariancia szerkezetét is. A rangkorreláció akkor lesz minimális a célfüggvény inverze és a normális eloszlás között, ha a paraméterek közötti korreláció közel azonos a két függvényben. Az eljárásunkkal az alábbi korrelációs mátrixot kaptuk az R1 és R2 reakció illesztett Arrhenius-paramétereire: lnA1 n1 E1/R lnA2 n2 E2/R lnA1 1.000 -1.000 0.999 0.156 -0.147 -0.165 n1 -1.000 1.000 -1.000 -0.158 0.150 0.168 E1/R 0.999 -1.000 1.000 0.169 -0.161 -0.177 lnA2 0.156 -0.158 0.169 1.000 -1.000 -0.983 n2 -0.147 0.150 -0.161 -1.000 1.000 0.983 E2/R -0.165 0.168 -0.177 -0.983 0.983 1.000 4. táblázat
3.4. Továbbfejlesztési lehetőségek Az algoritmusunk megad egy olyan paramétervektort, melynek meghatározásánál egyszerre vettük figyelembe a közvetett és közvetlen méréssel kapott kísérleti adatokat, és ezáltal az optimalizált paramétervektor jobban leírja a modellünket, mint a kizárólag közvetlen mérési eredményekből számított paraméterértékek. Egy olyan korrelációmátrixot is kapunk, mely jól leírja egy, a célfüggvénnyel fordított arányban álló függvény paramétereinek korrelációszerkezetét. Folyamatban van egy olyan módszer kidolgozása, melyben a mérések bizonytalanságát is figyelembe vesszük a paraméterek szórásainak meghatározásakor. Ezzel a 18
módszerrel a meglévő korrelációmátrixot felhasználva egy pontosabb kovarianciamátrixot kaphatunk a paraméterekre, mint amit jelenlegi módszerünk biztosít. A 3.2. alfejezetben megadott modell ugyan viszonylag jó leírást ad a hidrogén−levegő elegy égésére, de több kísérleti adattal és reakciólépés figyelembe vételével lényegesen pontosabb leírást adhatunk a jelenségről. Az algoritmusban a paraméterek számától függetlenül csak egy egydimenziós szélsőérték kereső eljárást futtatunk (γ-ban optimalizálunk), ezért nem kell attól tartanunk, hogy több paraméter esetén a módszer nagyon lelassulna, vagy esetleg le sem futna. Emiatt újabb szimulációkkal, kísérleti adatokkal és paraméterekkel szeretnénk kibővíteni a modellünket a hidrogén-levegő elegy égésének minél jobb leírásához. A továbbiakban szeretnénk a módszert bonyolultabb rendszerekre is alkalmazni. A következő lépésben alacsonyabb szénatomszámú szénhidrogének (például metán vagy etán) égését is szeretnénk modellezni. Egy másik elképzelés szerint a célfüggvény egy módosított alakját lehetne használni a modell paramétereinek meghatározására. A célfüggvény felépítése hasonló lenne, de most nem a paraméterekkel számolt sebességi együtthatónak és szimulált gyulladási időnek a mérési eredményektől való átlagos euklideszi távolságát adná meg, hanem azt, hogy hány mérési eredmény esik hibaszázalékon belül a számolt értékekre. A kirívó mérési adatok így nem tudnák elvinni a célfüggvény értékét, mert sok mérési adat esetén nagyon kis súly esne a nagyon rossz mérési eredményekre. A módszer jó működéséhez nagyon sok mérési adatra van szükség, ezért kevés mérési adat esetén az euklideszi távolságos módszert célszerű alkalmazni, míg sok kísérleti adat esetén a fentebb leírt módszer lehet a célravezetőbb. A módszert leíró algoritmust már programoztuk Matlab-ban és a tesztelése folyamatban van.
19
4. Összefoglalás A kémiai reakciók általában sok reakciólépésből tevődnek össze. Minden egyes reakciólépéshez tartozik egy sebességi együttható (k), amelynek segítségével, a tömeghatáskinetika feltételezésével kiszámítható a reakciólépés sebessége. Ez a sebességi együttható az esetek többségében hőmérsékletfüggő és esetleg nyomásfüggő. A reakciómechanizmusok leírásához szükségünk van az egyes reakciólépésekhez tartozó sebességi együtthatók (k) meghatározására. Ezek mérése történhet közvetlenül a sebességi együtthatók meghatározására szolgáló kísérletekkel, vagy olyan kísérletekkel, amelyek eredménye csak az egész reakciómechanizmussal értelmezhető. A dolgozatban szereplő mérések közül ilyen utóbbi típusú volt a gyulladási idő mérése. Az ilyen méréseket a továbbiakban közvetett méréseknek fogjuk nevezni. Miután részben a reakciólépésekhez tartozó sebességi együtthatók határozzák meg a reakciómechanizmussal kapott szimulációs eredményeket, ezért a gyulladási idő mérése az egyes sebességi együtthatókra is tartalmaz információt. A reakciók sebességi együtthatójának meghatározására gyakran csak a közvetlenül a sebességi együtthatóra vonatkozó méréseket veszik figyelembe és figyelmen kívül hagyják a közvetett mérési eredményeket. Ezzel a módszerrel sok értékes adat elveszik, pedig azok segítséget adhatnának a sebességi együtthatók pontosabb meghatározására. Egy másik megközelítés szerint [14] ezzel ellentétben a közvetlen méréseket csak a becsült paraméterek elfogadási intervallumának meghatározására használják fel és csak a közvetett mérésekből határozzák meg a sebességi együtthatókat. Ennél a módszernél a közvetlen mérésekkel kapott adatokat hagyják figyelmen kívül és ezáltal szintén sok hasznos adatot veszítenek el. Az irodalomban nem találtunk olyan módszert, amely a közvetett és közvetlen mérési eredményeket is egyaránt figyelembe venné a reakciók sebességi együtthatóinak meghatározása során. Az általunk kifejlesztett módszerben a sebességi együtthatót meghatározó Arrhenius-paramétereket és egyéb, a vizsgált reakciókhoz tartozó paramétereket becsülhetjük a közvetett és közvetlen mérések segítségével. Eljárásunk során a paramétervektor olyan sűrűségfüggvényét keressük, mely kisebb értéket vesz fel nagyobb célfüggvény érték esetén és nagyobb értéket vesz fel kisebb célfüggvény érték esetén. Módszerünkben többdimenziós normális eloszlással közelítjük a kiválasztott reakciók Arrhenius-paramétereinek együttes eloszlását. A probléma megoldására egy iteratív algoritmust adtunk, mely minden egyes iterációs lépésben a legnagyobb rangkorrelációs együtthatójú, és így a célfüggvényt fordított arányban legjobban leíró normális eloszlást választja ki. Az algoritmus egy olyan eloszlást adott eredményül, melynek várható értéke és így normális eloszlásról révén szó legnagyobb súlyú paramétervektora közel minimalizálja a célfüggvény értékét. Ennek ellenőrzéséhez egy szélsőérték kereső eljárást indítottunk az optimalizált paramétervektorból. Az eljárásunk továbbá egy korrelációs mátrixot is szolgáltatott a paraméterekre, amely segítségével képet kaphatunk a kiválasztott reakciólépések Arrhenius-paramétereinek kapcsolatáról. A paraméterek viselkedésének jobb megértéséhez továbbá egy a célfüggvénnyel fordítottan arányos függvény kétdimenziós megszorításait is vizsgáltuk, amelyekből a paraméterek korrelációjára kapott értékeket tudjuk leellenőrizni. A program futásidejének csökkentésére, illetve a minél homogénebb mintavételezéshez a Monte Carlo latin hiperkocka mintavételezési eljárást alkalmaztuk. Az így a kapott minta
20
jobban szóródik, mint hagyományos mintavételezés esetén, és ezáltal a célfüggvényértékek feltérképezése is hatásosabb. A módszert alkalmaztuk a hidrogén gyulladása két kritikus reakciójának, a H + O2 → O + OH H + O2 (+M) → HO2 (+M)
(R1) (R2)
reakciók Arrhenius-paramétereinek meghatározására. Ehhez felhasználtuk az irodalomból Petersen és munkatársai [29] valamint Slack és munkatársai [30] gyulladási idő méréseit, továbbá Joe V. Michael és munkatársai [31-32] sebességi együtthatójának meghatározását. Pontosabban az (R2) reakció esetén a közvetlen mérési adatok a nagynyomású határértékre vonatkoztak és a közvetett mérések körülményeinél végzett szimulációs eredmények is az (R2) reakció nagynyomású határérték sebességi együtthatójának Arrhenius-paramétereire voltak érzékenyek. A módszerrel a következő reakciókinetikai paramétereket határoztuk meg: H + O2 → O + OH reakció ln (A/cm3mol-1s-1) n E/R 38,61 -0,771 8478 K
H + O2 (+M) → HO2 (+M) reakció ln (A/cm3mol-1s-1) n E/R 50,79 -2,15 -70,7 K
Az algoritmusunk segítségével az alábbi korrelációmátrixot számítottuk ki a kiválasztott (R1) és (R2) reakciólépések származtatott Arrhenius-paramétereire: lnA1 n1 E1/R lnA2 n2 E2/R lnA1 1.000 -1.000 0.999 0.156 -0.147 -0.165 n1 -1.000 1.000 -1.000 -0.158 0.150 0.168 E1/R 0.999 -1.000 1.000 0.169 -0.161 -0.177 lnA2 0.156 -0.158 0.169 1.000 -1.000 -0.983 n2 -0.147 0.150 -0.161 -1.000 1.000 0.983 E2/R -0.165 0.168 -0.177 -0.983 0.983 1.000
21
Irodalom [1] [2] [3]
[4]
[5] [6] [7] [8] [9]
[10] [11]
[12] [13]
[14]
Turányi Tamás Reakciómechanizmusok vizsgálata Akadémiai Kiadó, Budapest (2010) http://kinetics.nist.gov/kinetics/index.jsp (utolsó elérés: 2010. március 9) Miller, J.A., M.J. Pilling, J. Troe Unravelling combustion mechanisms through a quantitative understanding of elementary reactions. Proceedings of the Combustion Institute, 30, 43-88 (2005) T. Turányi, L. Zalotai, S. Dóbé, T. Bérces Effect of the uncertainty of kinetic and thermodynamic data on methane flame simulation results Phys.Chem.Chem.Phys., 4, 2568-2578 (2002) Zsély, I.G., J. Zádor, T. Turányi Uncertainty analysis of NO production during methane combustion International Journal of Chemical Kinetics, 40, 754-768 (2008) Zádor, J., I.G. Zsély, and T. Turányi Local and global uncertainty analysis of complex chemical kinetic systems Reliability Engineering and System Safety, 91, 1232-1240 (2006) Zsély, I.G., J. Zádor, T. Turányi Uncertainty analysis of updated hydrogen and carbon monoxide oxidation mechanisms Proceedings of the Combustion Institute, 30, 1273-1281 (2005) Zádor, J., I.G. Zsély, T. Turányi, M. Ratto, S. Tarantola, A. Saltelli Local and global uncertainty analyses of a methane flame model Journal of Physical Chemistry A, 109: p. 9795-9807 (2005) Zádor, J., T. Turányi, K. Wirtz, M. J. Pilling Quantitative assessment of uncertainties for a model of tropospheric ethene oxidation using the European Photoreactor (EUPHORE) Atmospheric Environment, 39, 2805-2817 (2005) Li, J., Zhao, Z.; Kazakov, A.; Dryer, F. L. An updated comprehensive kinetic model of hydrogen combustion International Journal of Chemical Kinetics, 36, 566-576 (2004) O Conaire, M., J.M. Simmie, H.J. Curran Experiments used to validate kinetic mechanisms: an appraisal of N2 as a bathgas and interpreting selected experiments Proceedings of the European Combustion Meeting (2005) Konnov, A.A., Remaining uncertainties in the kinetic mechanism of hydrogen combustion Combustion and Flame, 152, 507-528 (2008) Hughes, K.J., T. Turányi, A. Clague, M.J.Pilling Development and testing of a comprehensive chemical mechanism for the oxidation of methane Int.J.Chem.Kinet., 33, 513-538 (2001) Smith, G.P. , D. M. Golden, M. Frenklach, N. W. Moriarty, B. Eiteneer, M. Goldenberg, C. T. Bowman, R. K. Hanson, S. Song, W. C. Gardiner, Jr., V. V. Lissianski, Z. Qin GRI-Mech 3.0. Elérési hely: http://www.me.berkeley.edu/gri_mech/.
22
(utolsó elérés: 2010. március 9) [15] Frenklach, M. H. Wang, C.-L. Yu, M. Goldenberg, C.T. Bowman, R.K. Hanson, D.F. Davidson, E.J. Chang, G.P. Smith, D.M. Golden, W.C. Gardiner and V. Lissianski GRI-Mech 1.2. Elérési hely: http://www.me.berkeley.edu/gri_mech/. (utolsó elérés: 2010. március 9) [16] Bowman, C.T. , R.K. Hanson, D.F. Davidson, W.C. Gardiner, Jr., V. Lissianski, G.P. Smith, D.M. Golden, M. Frenklach and M. Goldenberg GRI-Mech 2.11. Elérési hely: http://www.me.berkeley.edu/gri_mech/. (utolsó elérés: 2010. március 9) [17] Nagy Tibor Reakciókinetikai modellek bizonytalanságanalízise és redukciója Ph.D. értekezés, ELTE (2009) [18] http://mathworld.wolfram.com/MonteCarloMethod.html (utolsó elérés: 2010. március 9) [19] http://omega.albany.edu:8008/cdocs/ (utolsó elérés: 2010. március 9) [20] http://mathworld.wolfram.com/SpearmanRankCorrelationCoefficient.html (utolsó elérés: 2010. március 9) [21] Hervé Abdi The Kendall rank correlation coefficient Elérési hely: http://www.utdallas.edu/~herve/Abdi-KendallCorrelation2007-pretty.pdf (2007) (utolsó elérés: 2010. március 9) [22] Lagarias, J.C., J. A. Reeds, M. H. Wright, P. E. Wright Convergence properties of the nelder-mead simplex method in low dimensions SIAM Journal of Optimization, 9,, 112-147(1998) [23] MATLAB Help [24] http://en.wikipedia.org/wiki/Multivariate_normal_distribution (utolsó elérés: 2010. március 9) [25] http://en.wikipedia.org/wiki/Eigendecomposition_(matrix) (utolsó elérés: 2010. március 9) [26] Kee, R.J., F.M. Rupley, J.A.Miller CHEMKIN-II: A FORTRAN Chemical kinetics package for the analysis of gas-phase chemical kinetics Sandia National Laboratories (1989) [27] http://en.wikipedia.org/wiki/Gibbs_measure (utolsó elérés: 2010. március 9) [28] Sedyo, I. Néhány gázkinetikai reakció paraméterei bizonytalanságának hőmérsékletfüggése ELTE Kémiai Intézet, TDK dolgozat (2009) [29] Petersen, E. L., Davidson, D. F., Röhrig, M. Hanson, R. K. Shock Waves 941-946 (1996) [30] Slack, M. W. Rate coefficient for H + O2 + M = HO2 + M evaluated from shock tube measurements of induction times Combust. Flame 28, 241 (1977) [31] Pirraglia, A N, Michael, J V, Sutherland, J W, Klemm, R B
23
A flash photolysis-shock tube kinetic study of the H atom reaction with O2: H + O2 = OH + O (962K =< T =< 1705K) and H + O2 + Ar -> HO2 + Ar (746K =< T =< 987K) Journal of Physical Chemistry, 93, 282-291 (1989) [32] Kuan S. Shin and J. V. Michael Rate constants for the reactions H + O2 + OH + 0 and D + O2+ OD + O over the temperature range 1085-2278 K by the laser photolysis-shock tube technique Journal of Chemical Physics, 95, 262-273 (1991)
24
