江 凡，陳 陽，樊雪松，蔡國飆

（北京航空航天大學(xué)宇航學(xué)院，北京100191）

一種新穎的基于總能量守恒的化學(xué)平衡流算法

江 凡，陳 陽，樊雪松，蔡國飆

（北京航空航天大學(xué)宇航學(xué)院，北京100191）

在推進劑燃燒的建模上，傳統(tǒng)的熱力計算方法一般基于總焓守恒求解定壓絕熱燃燒溫度和平衡組分，不能考慮壁面?zhèn)鳠?；在燃氣流動的建模上，通常采用的凍結(jié)流模型認為本地的組分及熱物理性質(zhì)與燃燒室瞬時一致，忽略了這些參數(shù)因來流氣體與本網(wǎng)格滯留氣體摻混帶來的隨時間的緩變效應(yīng)。提出了一種新穎的可以考慮壁面?zhèn)鳠岬幕诳偰芰渴睾愕幕瘜W(xué)平衡流計算方法，運用Fortran2008語言，采用面向?qū)ο缶幊谭椒ń⒘嘶瘜W(xué)平衡流燃氣發(fā)生器管道的模塊化仿真模型，并將該模型應(yīng)用到一個包含42個組件的渦輪試驗臺氣路系統(tǒng)的建模與仿真中。與早期模型仿真結(jié)果及試驗數(shù)據(jù)的對比發(fā)現(xiàn)，新模型的仿真結(jié)果有一定改進，更加接近試驗數(shù)據(jù)。

熱力計算；化學(xué)平衡流；總能量守恒；渦輪試驗臺系統(tǒng)；數(shù)值仿真

0 引言

系統(tǒng)仿真[1-7]是研究發(fā)動機工作特性的重要方法，目前NASA的推進系統(tǒng)數(shù)值仿真NPSS（Numerical Propulsion System Simulation）項目代表著當(dāng)今發(fā)動機系統(tǒng)仿真的最高水平[1-2]。Amesim[3]，GFSSP[4]及Matlab/Simulink[5-6]等系統(tǒng)仿真軟件在發(fā)動機系統(tǒng)仿真領(lǐng)域應(yīng)用較多，為了在系統(tǒng)仿真的準確性和運算速度之間做到合理的平衡，通常采用一維甚至零維的仿真模型。

在發(fā)動機燃燒室及噴管內(nèi)流場仿真計算中，熱力計算[8-11]是一種重要方法，流動計算中凍結(jié)流模型[12-14]，化學(xué)非平衡流模型[15-17]應(yīng)用較多。

傳統(tǒng)的熱力計算方法一般采用總焓守恒求解，不能考慮傳熱。例如高玉閃使用布林克萊法對氫/氧推進劑進行熱力計算，確定各組分的摩爾數(shù)和燃燒產(chǎn)物溫度[9]。Kayadelen采用化學(xué)平衡常數(shù)法研究CαHβOγNδ類型燃料水噴注燃燒現(xiàn)象，并基于總焓守恒求解絕熱燃燒溫度[10]。Kady在二維的N-S方程中，采用化學(xué)平衡常數(shù)法基于總焓守恒求解燃燒室溫度[11]。

呂翔針對以H2O2(L)/JP-10(L)為推進劑的RBCC發(fā)動機，建立了火箭燃燒室熱力計算模型，并采用凍結(jié)流假設(shè)進行噴管熱力計算[14]。汪球采用5組分17步反應(yīng)機理，研究高焓激波風(fēng)洞中典型狀態(tài)下氣流的熱化學(xué)非平衡流動特性[15]。Wei采用Ar元素的相關(guān)反應(yīng)，研究Ar等離子在噴管非平衡流動中的膨脹過程[16]。

凍結(jié)流模型認為組分及物性不變，當(dāng)?shù)鼐W(wǎng)格的組分和熱物理性質(zhì)參數(shù)與上游瞬時一致，忽略了這些參數(shù)因來流氣體與本網(wǎng)格滯留氣體摻混帶來的隨時間的緩變效應(yīng)；化學(xué)非平衡流從化學(xué)反應(yīng)動力學(xué)出發(fā)，考慮反應(yīng)機理，雖然更加精細，但是會帶來剛性問題，在系統(tǒng)級仿真中使計算效率大大降低，并且其準確性取決于所采用反應(yīng)機理的詳細程度和適用范圍，例如高溫機理無法模擬低溫燃燒。相對于上述2種模型，平衡流[18]模型燃燒產(chǎn)物的組分摩爾分數(shù)及熱物性參數(shù)由當(dāng)?shù)鼐W(wǎng)格參數(shù)求解，且不會明顯降低計算效率。

運用模塊化建模思想[7,12]，文獻 [13]針對渦輪試驗臺氣路系統(tǒng)建立了包含42個組件的數(shù)值仿真模型，加熱器及其出口段（燃燒裝置，以下簡稱加熱器）燃燒區(qū)采用壁面絕熱的基于總焓守恒的熱力計算方法，流動區(qū)采用凍結(jié)流模型。加熱器總溫的動態(tài)仿真結(jié)果與試驗數(shù)據(jù)一致，但加熱器后續(xù)管道中的管流總溫明顯高于試驗數(shù)據(jù)。

針對以往研究中的不足，本文提出了一種基于總能量守恒的熱力計算方法，并從平衡流物理化學(xué)實質(zhì)出發(fā)，進一步提出了一套新穎的基于總能量守恒的化學(xué)平衡流計算方法，采用該方法改進文獻 [12-13]中的加熱器模型，并通過與試驗數(shù)據(jù)的對比，驗證新方法的有效性。

1 渦輪試驗臺裝置

圖1所示為由42個組件組成的渦輪試驗臺氣路系統(tǒng)的仿真模型，模塊化建模[13]時劃分為2個流體源 (FS1-2)，20個氣體管道 (GP1-20)，5個氣體容積 (GVol1-5)，1個氣體摻混器(GMA1)，9個氣體閥門 (GV1-9)，1個氣體噴注器 (GI1)，1個氣體穩(wěn)定器（GSA1），1個化學(xué)平衡流燃氣發(fā)生器管道 (CEFGGP1)，1個氣渦輪(GTurbo1)和1個轉(zhuǎn)子 (Rotor1)。Tt13/Tt14（2個傳感器在同一橫截面不同位置測量總溫，后面類似）、Tt51/Tt52及Tt33/Tt36分別為加熱器下游、放氣旁路管道GP16及主路測量段管道GP17總溫測點。

加熱器（CEFGGP1）、管道GP16、GP17的長度（單位：m）、外徑、壁厚（單位：mm）如圖1所示，其他元件結(jié)構(gòu)參數(shù)參見文獻 [13]，管道類元件每100 mm劃分一個網(wǎng)格，氣體閥門與氣體容積元件劃分為2個標準網(wǎng)格，總長200 mm。FS1和FS2分別為入口流體源和出口流體源。其中，管道GP15，閥門GV7及管道GP16組成放氣旁路，管道GP17，閥門GV8及管道GP18組成渦輪主路工作段為渦輪輸送高溫高壓燃氣，最終帶動轉(zhuǎn)子做功。

試驗和算例設(shè)置、調(diào)節(jié)時序情況參見文獻[13]，仿真時給定一個初始流場和溫度場分布，待冷流穩(wěn)定后，10 s時刻在加熱器中（CEFGGP1）點火，整個仿真過程持續(xù)260 s，重點關(guān)注燃燒裝置在20～260 s的仿真結(jié)果。FS1～2兩個流體源邊界條件采用試驗測量數(shù)據(jù)，即燃油質(zhì)量流量Qmf、大氣壓強p0的時間-參數(shù)點列形式的試驗測量數(shù)據(jù)。整個系統(tǒng)工作狀況隨閥門開度及流體源邊界條件而變化。本文針對其中的加熱器（CEFGGP1）開展研究工作。

2 加熱器的平衡流模型

平衡流模型假設(shè)認為在流動過程中，物質(zhì)化學(xué)反應(yīng)速率無限快，燃燒物質(zhì)時刻處于一種組分和能量的化學(xué)平衡狀態(tài)。

2.1 網(wǎng)格劃分和基本假設(shè)

加熱器型號為BK-1型單管燃燒室，分別以空氣、RP-3煤油為氧化劑和燃料。如圖2所示為根據(jù)加熱器建立的化學(xué)平衡流燃氣發(fā)生器管道模塊的網(wǎng)格圖，將加熱器劃分為摻混點火區(qū)與流動燃燒區(qū)。氧化劑和燃料進入摻混點火區(qū)發(fā)生化學(xué)反應(yīng)，然后燃燒物質(zhì)流入到流動燃燒區(qū)。

相較于文獻 [12-13]中的加熱器模型，基于總能量守恒的平衡流模型有以下不同的假設(shè)：

1）摻混點火區(qū)采用總能量微分方程代替壓強微分方程，燃燒過程不再是絕熱的，可以考慮壁面?zhèn)鳠幔?/p>

2）流動燃燒區(qū)組分和物性不再凍結(jié)不變，而是隨當(dāng)?shù)鼐W(wǎng)格的混合比等狀態(tài)參數(shù)而變化。

這2點假設(shè)的不同會帶來新模型中溫度及壓強算法的不同。

2.2 基于總能量守恒的平衡流模型

式中：nmolis,i為i網(wǎng)格1 kg物質(zhì)中組分is的摩爾數(shù)；ns為組分總數(shù)；his,i為1摩爾第is種組分的焓值（J/mol），是溫度的單值函數(shù)，可由NASA七系數(shù)或九系數(shù)擬合公式得到[18]。

上述方程為常微分方程的初值問題，沿時間推進求解時步驟如下：

1）流動燃燒區(qū)常微分方程變量Ev,i，ρi，fmix,及Wj由式 (5)～式 (8)沿時間積分求解，uj由式(9)得到，這些參數(shù)對于本時層的代數(shù)方程求解來說成為定值。對于代數(shù)方程變量來說：hi是溫度Ti和組分摩爾數(shù)nmolis,i的函數(shù)，而需要熱力計算得到的nmolis,i；Ri又是fmix,i，pi的函數(shù)，因此需聯(lián)立式 (12)～式 (15)迭代求解溫度Ti。摻混點火區(qū)Ev,0，ρ0及fmix,0由式 (1)、式 (3)及式 (4)沿時間積分求解，u0由上游組件決定，其他步驟相同。

2）不論冷流、熱流，均在式 (14)中采用二分法求解i網(wǎng)格的溫度Ti，二者不同之處在于焓值hi和組分數(shù)nmolis,i的計算。

熱流時選取試算溫度Tmiddle，在平衡組分初值庫中，由當(dāng)前時層混合比rof=1/fmix,i-1預(yù)估Tmiddle下平衡組分初值nmolg,i，由式 (12)和式 (13)得到試算壓強pmiddle，采用基于化學(xué)平衡常數(shù)法或最小吉布斯自由能法的熱力計算方法求解 Tmiddle和pmiddle下平衡組分值，由式 (15)計算燃燒混合氣的焓hi，并采用式 (13)更新Ri，最終根據(jù)式(14)采用二分法迭代求解溫度Ti。

冷流時，忽略煤油對冷流的貢獻，只考慮空氣流動。視空氣為一定比例的Ar、O2及N2三種組分的混合物，則nmolis,i和nmolg,i為定值，因此hi是溫度的單值函數(shù)，聯(lián)立式 (13)和式 (15)，最終根據(jù)式 (14)采用二分法迭代求解溫度Ti。

3）得到網(wǎng)格溫度Ti后，由式 (12)求壓強pi。

4）熱流時由熱力計算得到最終的燃燒產(chǎn)物組分數(shù)及熱物性參數(shù)；冷流時，各物質(zhì)組分即為空氣組分，物性參數(shù)由上游相連組件參數(shù)賦值。

3 模型應(yīng)用效果

文獻 [13]采用Fortran95語言和面向過程的編程方法，建立了渦輪試驗臺氣路系統(tǒng)的模塊化仿真模型。其中加熱器模型燃燒區(qū)采用壁面絕熱的基于總焓守恒的熱力計算模型，流動區(qū)采用基于總能量（不包含化學(xué)能）守恒的凍結(jié)流模型，傳熱模型采用文獻 [19]在圓柱坐標系下建立的一種計算管壁傳熱的軸對稱二維有限體積模型。此方案稱為舊方案。

本文采用Fortran2008語言編程，將試驗臺系統(tǒng)所涉及的主要模塊面向?qū)ο蠡?，并采用新的平衡流模型對加熱器進行改進，摻混點火區(qū)和流動燃燒區(qū)采用基于總能量（包含化學(xué)能）守恒的平衡流模型并考慮傳熱，此方案稱為新方案。

加熱器和氣體管道壁面?zhèn)鳠岵捎幂S對稱二維傳熱模型，氣體閥門、噴注器、氣體穩(wěn)定器及氣渦輪為徑向一維傳熱模型，氣體容積為零維傳熱模型，氣體摻混器為絕熱模型。

在編程過程中還進行了其他許多改進以優(yōu)化代碼，例如改進了熱力計算過程中確定溫度的二分法方案，基本思想是采用上一時層的溫度確定本時層溫度二分法求解初值，在收斂精度提高到0.01 K的情況下反而減小了二分法迭代步數(shù)。數(shù)值試算表明，面向?qū)ο蠡投址ǜ倪M使得計算速度提高了約20%。

如圖3～圖5所示為采用化學(xué)平衡常數(shù)法的2種方案仿真結(jié)果與3個測點處試驗數(shù)據(jù)的對比。

圖3中nT為轉(zhuǎn)子轉(zhuǎn)速，GGP-Tt(n/2)、CEFGGPTt(n/2)分別為舊方案和新方案加熱器下游中間網(wǎng)格處總溫。圖4中GP16-Tt(n/2)、GP16-Tt(n-1)分別為管道GP16中間網(wǎng)格與末端網(wǎng)格總溫。圖5中以此類推。

根據(jù)渦輪試驗臺氣路系統(tǒng)運行工況，在30 s，120 s及260 s分別截取3個測點的總溫數(shù)據(jù)，并與試驗數(shù)據(jù)進行對比，如表1和表2所示。

從表1可以看出，新方案的加熱器下游總溫相對于舊方案有一定下降，在30 s，120 s及260 s時，分別下降12.94 K，12.31 K及12.4 K。盡管如此，從圖3可以看出，2種方案的總溫仿真曲線整體上都落在了加熱器下游同一測點2個傳感器Tt13和Tt14測得的總溫曲線的范圍內(nèi)，說明新舊方案在加熱器下游總溫計算上都與試驗相符。

加熱器總溫下降一方面是摻混點火區(qū)采用式(1)求解總能量，在計算空氣、煤油焓值時，采用加熱器入口溫度計算；而舊方案中，統(tǒng)一采用煤油初溫Tinit計算，Tinit=282 K是煤油實測溫度的平均值?？梢姴捎眉訜崞魅肟跍囟雀咏咏鼘嶋H，空氣入口溫度在仿真涉及的260 s過程中約為270 K，相對于Tinit小12 K，引起仿真得到的加熱器總溫下降約10 K。另一方面，摻混點火區(qū)考慮壁面?zhèn)鳠?，?dǎo)致燃氣總溫又下降2～3 K。最終加熱器總溫相較原始方案約有12～13 K的下降。

表1 試驗數(shù)據(jù)與新舊方案總溫仿真結(jié)果Tab.1 Test data and total temperature simulation results of new and old cases

表2 新舊方案總溫仿真結(jié)果與試驗數(shù)據(jù)的對比結(jié)果Tab.2 Comparison of total temperature simulation results of new and old cases with test data

加熱器總溫下降引起了后續(xù)管道總溫的下降，如表1所示，在這3個時刻點，新方案放氣旁路管道GP16總溫相對于舊方案分別有8.06 K，7.8 K及9.41 K下降；主路測量段管道GP17總溫則分別有2.49 K，3.15 K及8.91 K下降。

從圖4可以看出，與試驗數(shù)據(jù)總溫Tt51和Tt52對比，新方案中放氣旁路總溫仿真結(jié)果更加接近試驗曲線。如圖5顯示與試驗數(shù)據(jù)總溫Tt33和Tt36對比，新方案中主路測量段總溫仿真結(jié)果更加接近試驗曲線；只是在62.9 s時刻GV8閥門開啟后，在約70～150 s時間段內(nèi)，仿真結(jié)果總溫突然上升，試驗總溫緩慢上升，二者趨勢仍不一致。分析原因為：鑒于當(dāng)前算例中后續(xù)組件仍然采用凍結(jié)流模型，組分及熱物理性質(zhì)參數(shù)根據(jù)加熱器末端網(wǎng)格計算結(jié)果直接賦值，導(dǎo)致后續(xù)組件網(wǎng)格中的這些參數(shù)是突變的，而實際情況應(yīng)該是物性隨時間緩變的。后續(xù)針對該差異，進一步細化加熱器后續(xù)管道模型，更好地描述氣路系統(tǒng)中溫度沿流動方向的分層現(xiàn)象。

4 結(jié)論

通過上述分析與比較，得出以下結(jié)論：

1）新方案加熱器摻混點火區(qū)分別采用空氣、煤油入口端溫度代替原始方案中的煤油初溫計算空氣和煤油焓值，更加接近物理化學(xué)實際，且導(dǎo)致加熱器下游總溫下降約10 K。

2）加熱器摻混點火區(qū)采用基于總能量守恒的平衡流模型，相對于傳統(tǒng)基于總焓守恒的方法，物理意義更加接近實際。其中，考慮壁面?zhèn)鳠釋?dǎo)致加熱器下游總溫計算結(jié)果較原始方案降低2～3 K。

3）加熱器下游總溫降低，導(dǎo)致后續(xù)組件中放氣旁路總溫、主路測量段總溫仿真曲線也整體降低，更加接近試驗曲線。

4）通過大量的文獻調(diào)研，尚未發(fā)現(xiàn)考慮壁面?zhèn)鳠岬幕诎瘜W(xué)能總能量守恒的類似方法，因此本文提出的方法具有一定的新穎性，在燃燒室及燃氣流動裝置建模與仿真領(lǐng)域具有一定的工程應(yīng)用價值。

5）后續(xù)考慮對加熱器后續(xù)管道進行平衡流模型改造，以期能夠?qū)崿F(xiàn)當(dāng)?shù)鼐W(wǎng)格組分及物性參數(shù)的緩變效應(yīng)，更好地描述渦輪試驗臺氣路系統(tǒng)中溫度沿流動方向的分層現(xiàn)象。

[1]HEATH C M.Characterization of swirl-venturi lean direct injection designs for aviation gas turbine combustion[J]. Journal of propulsion and power,2014,30(5)：1334-1356.

[2]CHAPMAN J W.A process for the creation of T-MATS propulsion system models from NPSS data：AIAA 2014-3931[R].USA：AIAA,2014.

[3]晏政,劉澤軍,程玉強,等.航天器推進系統(tǒng)模塊化建模方法[J].國防科技大學(xué)學(xué)報,2012,34(4)：28-32.

[4]BANDYOPADHYAY A,MAJUMDAR A.Network flow simulation of fluid transients in rocket propulsion systems [J].Journalofpropulsionandpower,2014,30(6)：1646-1653.

[5]OSIPOV V V,DAIGLE M J,MURATOV C B,et al.Dynamical model of rocket propellant loadingwith liquid hydrogen[J].Journal of spacecraft and rockets,2011,48(6)：987-998.

[6]FREDERICK D K.A new method for constructing fast models of jet engines in Simulink：AIAA 2009-5419[R]. USA：AIAA,2009.

[7]陶金偉,黃金泉,周文祥.航空發(fā)動機組態(tài)建模技術(shù)研究[J].航空動力學(xué)報,2010,25(10)：2372-2378.

[8]楊亞晶,何茂剛,徐厚達.水沖壓發(fā)動機的熱力性能分析[J].推進技術(shù),2009,30(2)：240-245.

[9]高玉閃,杜正剛,金平,等.熱沉式燃燒室熱載分析[J].航空動力學(xué)報,2010,25(11)：2619-2624.

[10]KAYADELEN H K,UST Y.Prediction of equilibrium products and thermodynamic properties in H2O injected combustion for CαHβOγNδ type fuels[J].Fuel,2013, 113：389-401.

[11]AL KADY M A，OWIS F M.Numerical simulation of nozzle flow with chemical equilibrium[J].Journal of engineering and applied sciences,2013,8(5)：11-18.

[12]陳陽,唐振宇,蔡國飆,等.渦輪試驗臺燃燒裝置建模與仿真(Ⅰ)基于混合物分數(shù)的分區(qū)模型[J].推進技術(shù), 2012,33(6)：981-986.

[13]CHEN Y,CAI G B,ZHANG Z Z,et al．Multi-field coupling dynamic modeling and simulation of turbine test rig gas system[J].Simulation modelling practice and theory, 2014,44(4)：95-118.

[14]呂翔,何國強,劉佩進,等.RBCC發(fā)動機性能分析模型改進方法研究[J].固體火箭技術(shù),2010,33(4)：387-390．

[15]汪球,趙偉,滕宏輝,等.高焓激波風(fēng)洞噴管流場非平衡特性研究[J].空氣動力學(xué)學(xué)報,2015,33(1)：66-71．

[16]WEI F Z,WANG H X,MURPHY A B,et al.Numerical modelling of the nonequilibrium expansion process of argon plasma flow through a nozzle[J].Journal of physics D applied physics,2013,46(50)：505205.

[17]SOLSVIK J,HAUG-WARBERG T,JAKOBSEN H A. Implementation of chemical reaction equilibrium by Gibbs and Helmholtz energies in tubular reactor models：application to the steam-methane reforming process[J]. Chemical engineering science,2016,140：261-278.

[18]GORDON S,MCBRIDE B J.Computer program for calculation of complex chemical equilibrium compositions and applications：NASA RP-1311[R].USA：NASA,1996.

[19]陳陽,高芳,張振鵬,等.準一維可壓縮瞬變管流的有限體積模型(Ⅱ)管壁溫度場的有限體積模型[J].航空動力學(xué)報,2008,23(2)：317-322.

（編輯：馬 杰）

A novel chemical equilibrium flow algorithm based on total energy conservation

JIANG Fan,CHEN Yang,FAN Xuesong,CAI Guobiao
（School of Astronautics,Beihang University,Beijing 100191,China）

As for the modeling of propellant combustion,the traditional thermodynamic algorithms generally use the total enthalpy conservation to solve constant-pressure adiabatic combustion temperature and equilibrium compositions,which can not consider the wall heat transfer. As for the modeling of combustion-gas flow,the frozen-flow model is usually employed with the assumption that local compositions and the thermophysical properties are instantaneous consistent with those in the combustion chamber,which ignores the slowly-varying effect of these parameters caused by the mixing process of the incoming flow and the residual gas in local grid.In view of the abovementioned weaknesses,a novel chemical equilibrium flow algorithm based on total energy conservation is proposed,which can consider the wall heat transfer.By using Fortran 2008 language and object-oriented programming method,a modular simulation model of chemical equilibrium flow gas generator pipe is established,and is applied to the modeling and simulation of a turbine test rig gassystem with 42 components.In comparison with the earlier simulation results and experimental data, the simulation results of the new model make a certain improvement and are closer to the experimental data.

thermodynamic calculation;chemical equilibrium flow;total energy conservation; turbine test rig system;numerical simulation

V430-34

A

1672-9374(2017)01-0024-08

2016-09-26；

2016-11-27

國家自然科學(xué)基金（11101023）

江凡（1991—），男，碩士，研究領(lǐng)域為液體推進系統(tǒng)動力學(xué)與仿真