數(shù)值解法在燃燒室出口燃?xì)鉁囟扔?jì)算中的應(yīng)用

韓 冰，王明瑞，李亞娟，馬 征

（中國航發(fā)沈陽發(fā)動(dòng)機(jī)研究所，沈陽110015）

0 引言

隨著航空發(fā)動(dòng)機(jī)推重比的不斷提升，燃燒室出口溫度越來越高，燃?xì)夥治龇ㄗ鳛?種間接測溫方法，在高溫升燃燒室出口溫度場的測試中應(yīng)用越來越廣泛[1]。燃?xì)夥治鰷y溫法是集熱力學(xué)、數(shù)值解技術(shù)和程序設(shè)計(jì)于一體的工程測試技術(shù)。高溫燃燒產(chǎn)物是多成分的復(fù)雜混合物，根據(jù)熱力學(xué)知識(shí)利用燃?xì)庵懈鞒煞纸⑽镔|(zhì)守恒、化學(xué)平衡、壓力和能量守恒方程構(gòu)成非線性方程組，對(duì)其求解即可得到化學(xué)平衡組分和燃?xì)鉁囟萚2]。

目前非線性方程組的求解沒有通用的解析方法，只能通過數(shù)值計(jì)算方法獲得，非線性方程組的解法長期以來都是數(shù)值計(jì)算領(lǐng)域的重要研究內(nèi)容，已經(jīng)產(chǎn)生了許多有效求解理論和方法，牛頓-拉夫森法是求解非線性方程組的經(jīng)典算法[3-4]，突出優(yōu)點(diǎn)是收斂速度快，缺點(diǎn)是為了得到方程組確切解，需要保證給定的初值接近實(shí)際解，否則迭代不易收斂。針對(duì)燃燒室各種工作條件，比如燃燒室中不同的油氣比、工作壓力以及燃料組成等相關(guān)非線性方程組的計(jì)算求解，通用的初值很難確定，只能根據(jù)具體算例反復(fù)嘗試給定[5-7]。在求解非線性方程組時(shí)，如何確定初值和迭代是否收斂以及收斂速度的快慢是燃?xì)鉁囟扔?jì)算所關(guān)注的問題。本文采用雙變量迭代法和赫夫法，不僅易確定初值，而且2種算法的迭代收斂速度都很快，通過VB編程實(shí)現(xiàn)非線性方程組的求解方法，并將算法成功應(yīng)用于燃燒室出口燃?xì)鉁囟葴y試中。

1 非線性方程組的建立

燃燒室燃燒產(chǎn)生的高溫反應(yīng)產(chǎn)物是1種多成分的復(fù)雜混合物，對(duì)于燃燒室中的穩(wěn)態(tài)燃燒過程，通常采用這樣的假設(shè)：由于氣體的停留時(shí)間較長，燃?xì)獬煞值幕瘜W(xué)平衡能夠建立，可以認(rèn)為燃燒室出口燃?xì)馓幱诨瘜W(xué)平衡狀態(tài)。而碳?xì)淙剂虾涂諝獾娜紵a(chǎn)物平衡成分一般由 CO2、CO、H2O、Ar、O2、N2、H2、OH、NO、O、H、N共12種成分組成，根據(jù)這些成分計(jì)算燃燒溫度可滿足計(jì)算需求。

通過列出物質(zhì)守恒、壓力、化學(xué)平衡和能量守恒方程構(gòu)建非線性方程組[8-9]。

1.1 物質(zhì)守恒方程

設(shè)定燃料分子式為CHaObNc，根據(jù)1 kg原始反應(yīng)物與燃燒后生成的1 kg燃?xì)飧髟匚镔|(zhì)的量相等原則列出碳、氫、氧、氮、氬的守恒方程

式中：As=p·Mg，Mg為燃?xì)獾奈镔|(zhì)的量；f和 ω 分別為燃燒室油氣比和進(jìn)口空氣的質(zhì)量含濕量；RO2、SN2、TCO2和DAr分別為O2、N2、CO2和Ar在空氣中的體積分?jǐn)?shù)；Ma、Mf、Mω分別為空氣、燃料和水蒸氣的物質(zhì)的量；pi為燃?xì)庵懈鞒煞址謮海籒i表示1 kg原始反應(yīng)物中各元素的物質(zhì)的量。

1.2 壓力方程

式中：p為燃燒室的工作壓力。

1.3 化學(xué)平衡方程

根據(jù)化學(xué)反應(yīng)式（4）建立以平衡常數(shù)表示的化學(xué)平衡方程

式中：Kpi為按氣體分壓計(jì)算的化學(xué)平衡常數(shù)；pˉi為pi與標(biāo)準(zhǔn)狀態(tài)壓力p0之比，p0=101325 Pa。

1.4 能量守恒方程

碳?xì)淙剂吓c空氣燃燒的通用化學(xué)反應(yīng)式為

式中：n0為參與燃燒空氣的物質(zhì)的量，根據(jù)已知油氣比可通過式（7）求出；ng為燃燒生成燃?xì)獾奈镔|(zhì)的量，在計(jì)算出燃?xì)飧鹘M分的物質(zhì)的量百分?jǐn)?shù)后通過任意1種元素的物質(zhì)守恒方程求出，如可利用碳元素守恒式（8）求出；hH2O為燃燒室進(jìn)口空氣中的物質(zhì)的量含濕量。

式中：a、b與CHaObNc中的下標(biāo)是同一含義，為燃料中氫和氧的原子數(shù)；L0為單位質(zhì)量燃料完全燃燒所需理論空氣質(zhì)量；rCO2和rCO為燃?xì)庵蠧O2和CO的物質(zhì)的量百分?jǐn)?shù)，其中

根據(jù)式（6）建立能量守恒方程

式中：Qf為燃料的低熱值；η為燃?xì)鉄釕B(tài)平衡時(shí)的燃燒效率，根據(jù)燃?xì)庵形慈冀M分，包括離解組分可以計(jì)算其值；d0為定溫燃燒焓差；Hi,T為各組分在溫度T的焓值；T3為燃燒室進(jìn)口空氣溫度；T4為燃?xì)鉁囟取?/p>

式（1）、（3）、（5）、（9）構(gòu)成非線性方程組，其中未知變量包括12種氣體組分分壓pi、As和燃?xì)鉁囟萒4共14個(gè)。

2 數(shù)值解法的應(yīng)用

求解非線性方程組的根本目的是獲得燃?xì)鉁囟龋氲碾p變量迭代法和赫夫法主要用于求解化學(xué)平衡組分，下面具體介紹2種數(shù)值解法在燃?xì)鉁囟扔?jì)算中的應(yīng)用[10-11]。

2.1 雙變量迭代法

雙變量迭代法的基本方程式可以參考文獻(xiàn)[12]，該算法的基本原理是給出pN2和pO2的初值，通過該算法變換后的方程組計(jì)算出其他變量和pN2、pO2新的數(shù)值。經(jīng)過大樣本計(jì)算驗(yàn)證一般取pN2=0.780881p、pO2=0.209495p作為迭代初值，可以保證對(duì)燃燒室各種工作條件下的迭代計(jì)算都有效收斂。針對(duì)貧油燃?xì)猓捎秒p變量松弛法同時(shí)迭代方案計(jì)算更為快捷，迭代公式為

當(dāng)T＞3000 K或者富油時(shí)，采用雙變量同時(shí)迭代一般不易收斂，此時(shí)采用雙變量分別迭代可滿足計(jì)算收斂，即先固定pN2值，對(duì)pO2進(jìn)行迭代，然后再對(duì)pN2進(jìn)行迭代，如此反復(fù)迭代直至滿足精度要求為止。雙變量分別迭代法與同時(shí)迭代法的方程式形式基本相同，不同的是雙變量分別迭代公式采用雙點(diǎn)弦割法，如式（12）所示。雙變量分別迭代方法編程相比于雙變量同時(shí)迭代法要稍微復(fù)雜。

通常給定精度小于10-9p的情況下，雙變量迭代法的迭代次數(shù)在數(shù)次到十幾次。目前針對(duì)主燃燒室貧油燃?xì)獾墓ぷ鳁l件，雙變量同時(shí)迭代方案基本能滿足計(jì)算要求，而且其計(jì)算收斂時(shí)間比雙變量分別迭代要快得多。因此一般在計(jì)算主燃燒室出口溫度時(shí)采用雙變量同時(shí)迭代法基本能滿足計(jì)算收斂。

2.2 赫夫法

赫夫法是將牛頓法各方程式和變量取對(duì)數(shù)形式后求解非線性方程組，這樣處理使得燃?xì)饨M分的初值可以任意給定，此外該方法對(duì)任何原始反應(yīng)物和反應(yīng)產(chǎn)物都是通用的，迭代收斂很好，而且適用于富油燃燒[13-15]。

將式（1）整理成

式中：Δi為各變量的修正量；Ri、fi為過程計(jì)算量，無實(shí)際意義。

根據(jù)式（3）整理成

將式（5）整理成

將式（1）、（3）、（5）以及其中的變量先進(jìn)行對(duì)數(shù)處理后，對(duì)修正量來說處理后的由式（13）、（16）、（17）構(gòu)成的方程組是線性的，而對(duì)方程組的系數(shù)矩陣可以采用高斯主元消去法求解，首先需要給定氣體組分分壓和As的初值，可以選擇按照貧油完全燃燒產(chǎn)物的相應(yīng)數(shù)值給定，這樣可以減少迭代次數(shù)，但是這樣確定初值增加了計(jì)算量；可以任意給定1組初值，但是取值不能為0，一般在計(jì)算中取氮?dú)獾姆謮撼踔禐?.8 p，其余氣體和As的分壓初值為10-5p，經(jīng)過大量抽樣計(jì)算驗(yàn)證，這種初值取值能滿足絕大多數(shù)計(jì)算工況。經(jīng)高斯主元消去法得到各修正量并代入式（20）和（21）得到新的各分壓值，再將新的數(shù)值重新代入系數(shù)矩陣中計(jì)算，直至相鄰2次計(jì)算結(jié)果滿足給定的精度要求為止。在給定精度小于 的情況下，赫夫法的迭代次數(shù)在數(shù)次到二十幾次，該算法的迭代收斂時(shí)間要大于雙變量迭代法的。

2.3 溫度迭代計(jì)算

利用雙變量迭代法或赫夫法求解化學(xué)平衡成分的計(jì)算過程實(shí)質(zhì)是在求解由式（1）、（3）、（5）構(gòu)成的非線性方程組，計(jì)算的前提條件是燃燒室工作壓力、油氣比和燃?xì)鉁囟纫阎怯捎谟?jì)算過程中變量燃?xì)鉁囟仁俏粗模虼嗽陂_始迭代計(jì)算時(shí)，需要先假定1個(gè)燃?xì)鉁囟纫员愕?jì)算。將能量守恒方程式（9）整理為

圖1 燃?xì)鉁囟鹊挠?jì)算流程

式中：fT為溫度T的函數(shù)。

在計(jì)算出燃?xì)馄胶獬煞趾螅瑢⒏鞒煞值奈镔|(zhì)的量百分?jǐn)?shù)代入式（22），如果|fT|在限定的精度范圍內(nèi)，則結(jié)束計(jì)算；如果|fT|的值超出限定值，則采用牛頓法迭代公式（23）計(jì)算出新的溫度，再重新計(jì)算化學(xué)平衡組分，直至|fT|滿足精度要求即可。燃?xì)鉁囟鹊挠?jì)算流程如圖1所示。

3 算例計(jì)算

選用算例對(duì)采用雙變量同時(shí)迭代法、雙變量分別迭代法和赫夫法計(jì)算得到的燃?xì)馄胶饨M分和燃?xì)鉁囟冗M(jìn)行比較分析。設(shè)定燃料的氫碳比為1.923，低熱值為42650 kJ/kg，進(jìn)口溫度為500℃，燃燒室工作壓力為0.5 MPa，按貧油和富油2種情況油氣比分別取0.05453和0.07496。采用前面所述的數(shù)值解法得到的溫度計(jì)算結(jié)果見表1、2。

從表1、2中結(jié)果可得出：

（1）針對(duì)同一工作條件，采用數(shù)值解法得到的化學(xué)平衡組分和燃?xì)鉁囟然疽恢拢?/p>

（2）從計(jì)算時(shí)間上比較，由快到慢分別為雙變量同時(shí)迭代法、雙變量分別迭代法和赫夫法，但是在燃燒室溫度場試驗(yàn)中，測試軟件采集數(shù)據(jù)的頻率是1 s采集1組數(shù)據(jù)，因此任意1種數(shù)值解法在測試軟件中的應(yīng)用都滿足測試需求；

表1 數(shù)值解法的結(jié)果對(duì)比（油氣比0.05453）

表2 數(shù)值解法的結(jié)果對(duì)比（油氣比0.07496）

（3）富油條件下雙變量同時(shí)迭代法無法完成迭代計(jì)算。

在某5頭部扇形燃燒室出口溫度場試驗(yàn)中，將熱電偶與燃?xì)夥治龅臏y量結(jié)果進(jìn)行比較。燃燒室采用燃料的氫碳比為1.923，燃料低熱值為42650 kJ/kg，進(jìn)口溫度為400℃，工作壓力為0.55 MPa，油氣比為0.03。燃燒室中間3個(gè)頭部的熱電偶與燃?xì)夥治鰷y得的出口溫度場結(jié)果對(duì)比如圖2所示，具體數(shù)值見表3，其中燃?xì)夥治鼋Y(jié)果包括前面所述采用數(shù)值解法得到的溫度結(jié)果。

圖2 油氣比0.03時(shí)熱電偶與燃?xì)夥治鰷囟葘?duì)比

表3 熱電偶與燃?xì)夥治鰷囟冉Y(jié)果對(duì)比

從圖2和表3中可見：

（1）熱電偶所測溫度與燃?xì)夥治鲇?jì)算得到的溫度沿著周向變化的趨勢(shì)基本一致，應(yīng)用雙變量迭代法和赫夫法的燃?xì)夥治鼋Y(jié)果相同；

（2）燃?xì)夥治龅臏y溫結(jié)果偏大于熱電偶的測溫結(jié)果，最大相差89℃，最小相差0.4℃，平均相差34.9℃。

4 結(jié)論

本文利用雙變量迭代法和赫夫法迭代求解用于計(jì)算燃?xì)鉁囟鹊姆蔷€性方程組，得出如下結(jié)論：

（1）雙變量迭代法計(jì)算簡單，收斂速度快，初值易確定并能滿足所有條件下的計(jì)算收斂，針對(duì)不同的燃燒情況可以采用雙變量同時(shí)迭代和分別迭代法完成計(jì)算。針對(duì)貧油燃?xì)猓捎秒p變量同時(shí)迭代法收斂更快。

（2）赫夫法是對(duì)牛頓-拉夫森法的改進(jìn)優(yōu)化，變化后的赫夫法各變量初值幾乎可以任意給定，迭代收斂較好并且該方法適用范圍非常廣，缺點(diǎn)是每次迭代都需要重新計(jì)算1次系數(shù)矩陣，計(jì)算量較大。

（3）通過算例驗(yàn)證可知，所有數(shù)值解法得到的計(jì)算結(jié)果基本一致，但從計(jì)算迭代收斂的速度來看，雙變量同時(shí)迭代法計(jì)算最快，赫夫法計(jì)算相對(duì)較慢。

（4）針對(duì)某5頭部扇形燃燒室中間3個(gè)頭部出口溫度的對(duì)比結(jié)果可見，燃?xì)夥治鲇?jì)算得到的溫度大于熱電偶測得的溫度。