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

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

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

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

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

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

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

0 引言

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

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

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

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

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

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

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

1 渦輪試驗(yàn)臺(tái)裝置

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

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

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

2 加熱器的平衡流模型

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

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

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

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

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

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

這2點(diǎn)假設(shè)的不同會(huì)帶來新模型中溫度及壓強(qiáng)算法的不同。

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]。

上述方程為常微分方程的初值問題，沿時(shí)間推進(jìn)求解時(shí)步驟如下：

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

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

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

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

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

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

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

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

本文采用Fortran2008語言編程，將試驗(yàn)臺(tái)系統(tǒng)所涉及的主要模塊面向?qū)ο蠡⒉捎眯碌钠胶饬髂Ｐ蛯?duì)加熱器進(jìn)行改進(jìn)，摻混點(diǎn)火區(qū)和流動(dòng)燃燒區(qū)采用基于總能量（包含化學(xué)能）守恒的平衡流模型并考慮傳熱，此方案稱為新方案。

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

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

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

圖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ù)渦輪試驗(yàn)臺(tái)氣路系統(tǒng)運(yùn)行工況，在30 s，120 s及260 s分別截取3個(gè)測(cè)點(diǎn)的總溫?cái)?shù)據(jù)，并與試驗(yàn)數(shù)據(jù)進(jìn)行對(duì)比，如表1和表2所示。

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

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

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

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

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

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

4 結(jié)論

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

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

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

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

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

5）后續(xù)考慮對(duì)加熱器后續(xù)管道進(jìn)行平衡流模型改造，以期能夠?qū)崿F(xiàn)當(dāng)?shù)鼐W(wǎng)格組分及物性參數(shù)的緩變效應(yīng)，更好地描述渦輪試驗(yàn)臺(tái)氣路系統(tǒng)中溫度沿流動(dò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]晏政,劉澤軍,程玉強(qiáng),等.航天器推進(jìn)系統(tǒng)模塊化建模方法[J].國防科技大學(xué)學(xué)報(bào),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ā)動(dòng)機(jī)組態(tài)建模技術(shù)研究[J].航空動(dòng)力學(xué)報(bào),2010,25(10)：2372-2378.

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

[9]高玉閃,杜正剛,金平,等.熱沉式燃燒室熱載分析[J].航空動(dòng)力學(xué)報(bào),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]陳陽,唐振宇,蔡國飆,等.渦輪試驗(yàn)臺(tái)燃燒裝置建模與仿真(Ⅰ)基于混合物分?jǐn)?shù)的分區(qū)模型[J].推進(jìn)技術(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]呂翔,何國強(qiáng),劉佩進(jìn),等.RBCC發(fā)動(dòng)機(jī)性能分析模型改進(jìn)方法研究[J].固體火箭技術(shù),2010,33(4)：387-390．

[15]汪球,趙偉,滕宏輝,等.高焓激波風(fēng)洞噴管流場(chǎng)非平衡特性研究[J].空氣動(dòng)力學(xué)學(xué)報(bào),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]陳陽,高芳,張振鵬,等.準(zhǔn)一維可壓縮瞬變管流的有限體積模型(Ⅱ)管壁溫度場(chǎng)的有限體積模型[J].航空動(dòng)力學(xué)報(bào),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)域?yàn)橐后w推進(jìn)系統(tǒng)動(dòng)力學(xué)與仿真