利用PP法簡(jiǎn)化液體姿軌控發(fā)動(dòng)機(jī)化學(xué)反應(yīng)機(jī)理

王大銳，程圣清，張 楠

（北京航天動(dòng)力研究所，北京100076）

0 引言

近年來(lái)，耦合詳細(xì)化學(xué)反應(yīng)的液體火箭發(fā)動(dòng)機(jī)燃燒過(guò)程研究成為熱點(diǎn)，掌握主要推進(jìn)劑組合的化學(xué)反應(yīng)機(jī)理是開(kāi)展該研究的基礎(chǔ)。目前，甲基肼（MMH）/四氧化二氮（NTO）,液氫（LH2）/液氧（LO2）和煤油（Kerosene）/液氧（LO2）是3種主要的航天推進(jìn)劑組合。與后兩種推進(jìn)劑組合相比，MMH/NTO屬于自燃推進(jìn)劑，其詳細(xì)機(jī)理很難獲得。此前，對(duì)于此類(lèi)推進(jìn)劑組合火箭發(fā)動(dòng)機(jī)燃燒CFD計(jì)算通常采用總包反應(yīng)代替，如O.Knab將四步總包反應(yīng)加入到發(fā)動(dòng)機(jī)燃燒模型中[1]，對(duì)400 N液體姿軌控火箭發(fā)動(dòng)機(jī)燃燒過(guò)程進(jìn)行了數(shù)值模擬。研究后發(fā)現(xiàn)，總包反應(yīng)對(duì)于中間反應(yīng)忽略過(guò)大，數(shù)值計(jì)算結(jié)果與試驗(yàn)測(cè)量結(jié)果出入較大，因此需要構(gòu)建詳細(xì)機(jī)理對(duì)發(fā)動(dòng)機(jī)燃燒過(guò)程進(jìn)行模擬，其中Laurent構(gòu)建了一套含有82種組分和403步基元反應(yīng)的詳細(xì)機(jī)理（機(jī)理未公布）[2-3]，該機(jī)理通過(guò)與Seamans的爆炸理論在點(diǎn)火延遲時(shí)間上計(jì)算結(jié)果的比較[4]，驗(yàn)證了該機(jī)理的正確性，但該機(jī)理過(guò)于復(fù)雜，并不適合耦合在湍流場(chǎng)中進(jìn)行發(fā)動(dòng)機(jī)燃燒CFD計(jì)算。因此，需要將詳細(xì)機(jī)理進(jìn)行適當(dāng)簡(jiǎn)化，以滿(mǎn)足發(fā)動(dòng)機(jī)燃燒CFD計(jì)算要求。

目前對(duì)于復(fù)雜燃燒機(jī)理的簡(jiǎn)化大體可以分為2大類(lèi)，一類(lèi)是以去除影響較小的基元反應(yīng)來(lái)進(jìn)行簡(jiǎn)化，如敏感性分析法SA[5-6]，準(zhǔn)穩(wěn)態(tài)假定QSSA[7]，主成分分析法PCA[8-9]等；另一類(lèi)則是以去除不重要的組分來(lái)進(jìn)行簡(jiǎn)化，如直接關(guān)系圖法DRG[10]以及多代通量法PFA[11-12]等。以上各簡(jiǎn)化方法各有優(yōu)勢(shì)，但也都有局限性。

本文提出將PCA和PFA方法聯(lián)合在一起使用，對(duì)43種組分和201步基元反應(yīng)的甲基肼/四氧化二氮反應(yīng)機(jī)理進(jìn)行簡(jiǎn)化研究，并將簡(jiǎn)化機(jī)理與國(guó)外文獻(xiàn)以及詳細(xì)機(jī)理在典型工況下點(diǎn)火延遲時(shí)間、平衡溫度以及主要平衡產(chǎn)物計(jì)算結(jié)果進(jìn)行對(duì)比。最后，對(duì)比不同簡(jiǎn)化方法對(duì)MMH和NTO，CH4和O2以及H2和O2三種推進(jìn)劑組合化學(xué)反應(yīng)機(jī)理的簡(jiǎn)化結(jié)果，說(shuō)明PP法簡(jiǎn)化效果更好，適用范圍更廣。

1 主成分聯(lián)合多代通量機(jī)理簡(jiǎn)化方法

“第一步簡(jiǎn)化”先對(duì)含有m種組分和p步基元反應(yīng)機(jī)理的速度敏感性系數(shù)矩陣，其形式如下所示：

矩陣F由下式計(jì)算得到：

其中

Rj=kjrj(c)

式中：fi為組分i的生成速率；kj為第j步反應(yīng)的反應(yīng)速率常數(shù)；υi,j為基元反應(yīng)j中組分i的化學(xué)計(jì)量系數(shù)；Rj為第j步反應(yīng)的生成物的生成速率；rj(c)為第j步反應(yīng)中反應(yīng)物濃度的乘積；Fi,j為第j步反應(yīng)中組分i的生成或者消耗速率與組分i的凈生成率的比值。再將FTF對(duì)角化，得到矩陣

式中λ1,λ2,...,λp為FTF的特征值，且λ1>λ2>...>λp。根據(jù)主軸定理，有：

式中：λi為主成分Ψi中各基元反應(yīng)的重要性，即λi大所對(duì)應(yīng)的特征向量中，元素越大其代表的基元反應(yīng)越重要。此時(shí)，給定簡(jiǎn)化機(jī)理基元反應(yīng)個(gè)數(shù)p1，并要求第一步簡(jiǎn)化后機(jī)理的著火延遲時(shí)間和平衡溫度以及重要組分百分比與原機(jī)理計(jì)算結(jié)果在一定誤差之內(nèi)，具體表達(dá)式如下：

式中：PAB為組分B相對(duì)于A的生成相關(guān)系數(shù)；PA和CA分別為組分A的生成通量和消耗通量。類(lèi)似地，二代和三代相關(guān)系數(shù)表達(dá)式如下：

式中Mi1和Mi2為中間組分。

從文獻(xiàn) [12]上可以知道，理論上代數(shù)越多，簡(jiǎn)化結(jié)果越好，但通過(guò)對(duì)GRI-Mech3.0機(jī)理簡(jiǎn)化結(jié)果來(lái)看，四代計(jì)算精度提高并不明顯，所以本文仍采用三代通量計(jì)算來(lái)進(jìn)行簡(jiǎn)化。最后定義rAB為組分B對(duì)組分A的一代、二代及三代相關(guān)系數(shù)相加之和，表達(dá)式如下：

判斷組分B對(duì)A的重要程度時(shí)，將rAB與設(shè)定的閾值ε作比較，如果大于該閾值，則認(rèn)為該組分對(duì)組分A重要，予以保留，否則在簡(jiǎn)化機(jī)理時(shí)應(yīng)予以剔除。當(dāng)全部對(duì)A重要的其他組分被找到后，再以這些組分替換A，按照上面的方法繼續(xù)尋找對(duì)該組分重要的其他組分，直到最后一次尋找重要組分時(shí)沒(méi)有找到新的重要組分為止，將此時(shí)得到的簡(jiǎn)化機(jī)理2帶入式 (4)等驗(yàn)證，驗(yàn)證后的簡(jiǎn)化機(jī)理作為PP方法的最終簡(jiǎn)化結(jié)果。

2 簡(jiǎn)化結(jié)果與分析

本文詳細(xì)反應(yīng)機(jī)理MN-Detailed包含43種組分201步基元反應(yīng)，動(dòng)力學(xué)參數(shù)參考文獻(xiàn) [13]，熱力學(xué)參數(shù)參考Burcat熱力學(xué)庫(kù)[14]，燃燒數(shù)值求解由Chemkin軟件完成，簡(jiǎn)化機(jī)理由自編程序PP-Beta計(jì)算得出。

2.1 簡(jiǎn)化機(jī)理與文獻(xiàn)計(jì)算結(jié)果對(duì)比

首先對(duì)比詳細(xì)機(jī)理和簡(jiǎn)化機(jī)理與其他文獻(xiàn)計(jì)算結(jié)果。采用Chemkin封閉全混均質(zhì)反應(yīng)器對(duì)甲基肼/四氧化二氮自燃著火過(guò)程進(jìn)行模擬，模擬反應(yīng)工況為初始溫度為T(mén)=298 K，初始?jí)毫閜=242×104Pa，定容燃燒，得到不同機(jī)理著火溫度曲線如下圖1所示。

圖1 不同機(jī)理著火過(guò)程溫度對(duì)比Fig.1 Comparison of temperature in combustion process of different mechanisms

圖1中實(shí)線為PP法簡(jiǎn)化結(jié)果MN-Reduced機(jī)理著火曲線，該機(jī)理著火延遲時(shí)間為3.33 ms；虛線為MN-Detailed機(jī)理著火曲線，著火延遲時(shí)間為3.5 ms；點(diǎn)劃線為Catorie利用量子化學(xué)得到的Catorie-Detailed機(jī)理[2]，著火延遲時(shí)間為3.82 ms；Seamans利用爆炸理論計(jì)算結(jié)果為 3.6 ms[4]，MN-Detailed機(jī)理和MN-reduced機(jī)理著火延遲時(shí)間與量子化學(xué)詳細(xì)機(jī)理和爆炸理論計(jì)算的結(jié)果符合較好。

為進(jìn)一步驗(yàn)證簡(jiǎn)化結(jié)果的正確性，對(duì)比了不同壓力、不同混合比、不同初溫下MN-Reduced機(jī)理與Seamans根據(jù)爆炸理論著火延遲時(shí)間計(jì)算的結(jié)果，如圖2所示其中離散點(diǎn)為MN-Reduced機(jī)理計(jì)算結(jié)果。在溫度和混合比一定情況下，著火時(shí)間隨著壓力成負(fù)指數(shù)下降，而MN-Reduced機(jī)理計(jì)算壓力對(duì)著火時(shí)間影響規(guī)律與Seamans理論計(jì)算結(jié)果也符合較好。

圖2 不同壓力下著火延遲時(shí)間對(duì)比Fig.2 Comparison of ignition delay time at different pressure

2.2 發(fā)動(dòng)機(jī)特定燃燒工況下PP法簡(jiǎn)化結(jié)果與原機(jī)理計(jì)算結(jié)果比較

要將甲基肼/四氧化二氮簡(jiǎn)化機(jī)理運(yùn)用到液體火箭發(fā)動(dòng)機(jī)燃燒CFD中，還需要比較簡(jiǎn)化機(jī)理與原詳細(xì)機(jī)理在發(fā)動(dòng)機(jī)燃燒典型工況下的計(jì)算結(jié)果。液體火箭發(fā)動(dòng)機(jī)燃燒是定壓燃燒過(guò)程，所以著火過(guò)程采用封閉全混均質(zhì)反應(yīng)器，求解器采用定壓求解。表1中Φ為氧化劑和推進(jìn)劑質(zhì)量比。

表1 不同壓力下推進(jìn)劑燃燒初始條件Tab.1 Initial conditions of propellant combination combustion at different pressure

表1和圖3是對(duì)初始溫度一定、混合比一定時(shí)，不同壓力下簡(jiǎn)化機(jī)理與詳細(xì)機(jī)理對(duì)著火延遲時(shí)間和平衡溫度計(jì)算結(jié)果的對(duì)比。可以看出，簡(jiǎn)化機(jī)理與詳細(xì)機(jī)理在不同壓力下著火過(guò)程均符合較好，著火延遲相差均在3×10-5s以?xún)?nèi)，平衡溫度相差0.1 K以?xún)?nèi)。

圖3 不同壓力下簡(jiǎn)化機(jī)理和詳細(xì)機(jī)理著火過(guò)程比較Fig.3 Combustion process comparison of simplified mechanism and detailed mechanism at different pressure

表2和圖4是對(duì)初始溫度一定、燃燒壓力一定，不同混合比下簡(jiǎn)化機(jī)理與詳細(xì)機(jī)理對(duì)著火延遲時(shí)間和平衡溫度計(jì)算結(jié)果的對(duì)比?？梢钥闯?，簡(jiǎn)化機(jī)理與詳細(xì)機(jī)理在不同混合比下著火過(guò)程均符合較好，著火延遲相差均在2×10-5s以?xún)?nèi)，平衡溫度相差0.1 K以?xún)?nèi)。

表2 不同氧燃比下推進(jìn)劑燃燒初始條件Tab.2 Initial conditions of propellant combination combustion at different mixing ratio

圖4 不同氧燃比下簡(jiǎn)化機(jī)理與詳細(xì)機(jī)理燃燒過(guò)程比較Fig.4 Combustion process comparison of reduced mechanism and detailed mechanism at different Φ

圖5 不同壓力下簡(jiǎn)化機(jī)理與詳細(xì)機(jī)理H2O摩爾百分?jǐn)?shù)計(jì)算結(jié)果對(duì)比Fig.5 Comparison of H2O mole fraction calculation results of simplified mechanism and detailed mechanism at different pressure

圖6 不同壓力下簡(jiǎn)化機(jī)理與詳細(xì)機(jī)理N2摩爾百分?jǐn)?shù)計(jì)算結(jié)果對(duì)比Fig.6 Comparison of N2mole fraction calculation results of simplified mechanism and detailed mechanism at different pressure

圖5和圖6是質(zhì)量混合比1.65，初始溫度298 K，不同壓力下主要平衡產(chǎn)物H2O和N2簡(jiǎn)化機(jī)理和詳細(xì)機(jī)理計(jì)算結(jié)果的對(duì)比?？梢钥闯觯麄€(gè)著火過(guò)程中簡(jiǎn)化機(jī)理計(jì)算結(jié)果與詳細(xì)機(jī)理計(jì)算結(jié)果均符合較好。

綜上，利用PP法得到的簡(jiǎn)化機(jī)理在著火延遲時(shí)間、平衡溫度以及平衡組分摩爾百分?jǐn)?shù)都與詳細(xì)機(jī)理計(jì)算結(jié)果保持很好的一致性，這為耦合MMH/NTO化學(xué)反應(yīng)的發(fā)動(dòng)機(jī)燃燒CFD研究奠定了良好的基礎(chǔ)。

2.3 PP法與主成分分析和多代通量方法的比較

主成分法 (PCA)和多代通量分析法 (PFA)都是比較常用的化學(xué)反應(yīng)機(jī)理簡(jiǎn)化方法，但各自在使用時(shí)均有一定的局限性，本文采用PCA，PFA 及 PP法對(duì)MMH/NTO，H2/O2[14]及CH4/O2[15]推進(jìn)劑化學(xué)反應(yīng)詳細(xì)機(jī)理進(jìn)行簡(jiǎn)化，對(duì)比簡(jiǎn)化結(jié)果如表3所示。

表3中，MMH/NTO詳細(xì)機(jī)理在質(zhì)量混合比1.6～1.8，壓力 1.6～2.0 MPa，初始溫度 298 K 范圍內(nèi)簡(jiǎn)化；H2/O2詳細(xì)機(jī)理在摩爾混合比1～4，壓力0.1～1 MPa，著火溫度為1 000 K范圍內(nèi)簡(jiǎn)化；CH4/O2詳細(xì)機(jī)理在摩爾混合比1～2，壓力0.1～1 MPa,著火溫度為1 300～1 500 K范圍內(nèi)簡(jiǎn)化。從表3中可看出：一方面PCA簡(jiǎn)化機(jī)理結(jié)果基元反應(yīng)數(shù)少于PFA簡(jiǎn)化結(jié)果，但組分種類(lèi)要明顯多于PFA，這是因?yàn)镻CA通過(guò)去除非重要基元反應(yīng)達(dá)到簡(jiǎn)化詳細(xì)機(jī)理目的，而PFA則是通過(guò)去除非重要組分種類(lèi)來(lái)簡(jiǎn)化詳細(xì)機(jī)理，PP法將它們的特點(diǎn)結(jié)合起來(lái)，既能夠去除非重要基元反應(yīng)又能夠去除非重要組分，因而得到更為簡(jiǎn)單的簡(jiǎn)化機(jī)理；另一方面，從H2/O2和CH4/O2簡(jiǎn)化結(jié)果來(lái)看，PCA、PFA存在一定局限性，只能對(duì)一部分推進(jìn)劑反應(yīng)機(jī)理進(jìn)行簡(jiǎn)化，而PP法可以簡(jiǎn)化機(jī)理的范圍比前兩者都大，具有發(fā)展成一種通用簡(jiǎn)化方法的潛力。

表3 不同簡(jiǎn)化方法對(duì)3種推進(jìn)劑組合簡(jiǎn)化結(jié)果對(duì)比Tab.3 Comparison of simplified results of three combined propellants using different simplification methods

3 結(jié)論

本文運(yùn)用PP方法對(duì)MMH/NTO推進(jìn)劑組合43種組分201步基元詳細(xì)機(jī)理進(jìn)行了簡(jiǎn)化，簡(jiǎn)化機(jī)理與國(guó)外文獻(xiàn)以及詳細(xì)機(jī)理在多種典型工況下著火延遲時(shí)間、平衡溫度以及平衡組分摩爾百分?jǐn)?shù)計(jì)算結(jié)果都符合較好，說(shuō)明簡(jiǎn)化機(jī)理可以有效代替詳細(xì)機(jī)理進(jìn)行燃燒計(jì)算，并為進(jìn)行耦合化學(xué)反應(yīng)的液體姿軌控火箭發(fā)動(dòng)機(jī)燃燒CFD研究奠定了良好基礎(chǔ)；PP方法相比于PCA和PFA方法具備可同時(shí)減少基元反應(yīng)個(gè)數(shù)和組分種類(lèi)的特點(diǎn)，簡(jiǎn)化效果好，簡(jiǎn)化機(jī)理適用性也更廣。

[1]KNAB O,PRECLIK D,ESTUBLIER D.Flow field prediction within liquid film cooled combustion chambers of storable bipropellant rocket engines,AIAA 98-3370[R].USA：AIAA,1998.

[2]CATOIRE L,CHAUMEIX N,PAILLARD C.Chemical kinetic model for monomethylhydrazine/nitrogen tetroxide gas-phase combustion and hypergolic ignition[J].Journal of Propulsion and Power,2004,20(1)：87-92.

[3]CATOIRE L,SWIHART M T.Thermochemistry of species produced from monomethylhydrazine in propulsion andspace-relatedapplications[J].Journal ofPropulsion and Power,2002,18(6)：1242-1252.

[4]SEAMANS T F,VANPEE M,AGOSTA V D.Development of a fundamental Model of hypergolic ignition in space-ambient engines[J].AIAA Journal,1967,l5(9)：1616-1624.

[5]TURANYI T,BERCES T,VAIDA S.Reaction rate analysis of complex kinetic system[J].International Journal of Chemical Kinetics,1989,21(2)：83-99.

[6]喬瑜,徐明厚.基于敏感性分析的甲烷反應(yīng)機(jī)理優(yōu)化簡(jiǎn)化[J].華中科技大學(xué)學(xué)報(bào)(自然科學(xué)版),2007,35(5)：85-87.

[7]PETERS N,WILLIAMS F A.The asymptotic structure of stoichiometric methane-airflame[J].Combustion and flame,1987,68(2)：185-207.

[8]LU T F,CHUNG K L.A directed relation graph method for mechanism reduction[J].Proceedings of the Combustion Institute,2005,30：1333-1341.

[9]VAJDA S,VALKO P,TURANYI T.Principle component analysis of kinetic model[J].International Journal of Chem-Kinet,1985,17：55-81.

[10]董剛,邱榕,蔣勇,張和平.利用主成分分析法簡(jiǎn)化甲烷/空氣層流預(yù)混火焰的反應(yīng)機(jī)理[J].火災(zāi)科學(xué),2004,13(3)：158-162.

[11]SUN W T,CHEN Z,GOU X L,et al.A path flux analysis method for the reduction of detailed chemical kinetic mechanisms[J].Combustion and Flame,2010,157(7)：1298-1307.

[12]茍小龍,王衛(wèi),桂瑩,施萬(wàn)玲.一種多代路徑通量分析化學(xué)機(jī)理簡(jiǎn)化方法[J].推進(jìn)技術(shù),2012,33(3)：412-417.

[13]ANDERSON W R.,MCQUAID M J,MICHAEL J,et al.A detailed,finite-rate,chemical kinetics mechanism for monomethylhydrazine-red fuming nitric acid systems,ARL-TR88-502010[R].[S.l.]：ARL-TR,1988.

[14]MATSUO A,FUJII K,FUJIWARA T.Flow features of shock-induced combustion around projectile traveling at hypervelocities[J].AIAA Journal,1995,33(6)：1056-1063.

[15]YANG B,POPE S B.An investigation of the accuracy of manifold methods and splitting schemes in the computational implementation of combustion chemistry[J].Combustion and Flame,1998,112：16-32.

[16]聶萬(wàn)勝,豐松江.液體姿軌控火箭發(fā)動(dòng)機(jī)燃燒動(dòng)力學(xué)模型與數(shù)值模擬[M].北京：國(guó)防工業(yè)出版社,2011.