硝酸羥胺基單組元發(fā)動機(jī)起動過程數(shù)值模擬

孫得川，金東洙，于澤游

(大連理工大學(xué)航空航天學(xué)院，工業(yè)裝備結(jié)構(gòu)分析國家重點(diǎn)實(shí)驗(yàn)室， 遼寧 大連 116024)

單組元液體火箭發(fā)動機(jī)系統(tǒng)簡單、可靠性高，廣泛應(yīng)用在空間飛行器上[1]。我國現(xiàn)有型號的單組元發(fā)動機(jī)大多采用無水肼作為推進(jìn)劑，它不僅昂貴、易燃、易爆，而且有劇毒，不利于人員操作，因此新型無毒單組元發(fā)動機(jī)是當(dāng)前單組元火箭發(fā)動機(jī)的主要研制方向之一。無毒單組元推進(jìn)劑包括過氧化氫、氧化亞氮、硝酸羥胺(HAN)基、二硝酰胺銨(ADN)基、和硝酸烷基酯類(如硝酸異丙酯，IPN)的推進(jìn)劑等[2-7]。目前對無毒推進(jìn)劑單組元發(fā)動機(jī)的研究主要集中在推進(jìn)劑的催化分解研究和推進(jìn)劑性能實(shí)驗(yàn)研究方面。例如，Amrousse[7-8]、Katsumi[9]、Chang Hwan Hwang[10]等人對AND和HAN基推進(jìn)劑的熱分解和燃燒進(jìn)行了研究；陳永康[11]、潘玉竹[12]、任曉光[13]分別對HAN基液體推進(jìn)劑的熱分解動力學(xué)、高壓燃燒特性和催化分解進(jìn)行了研究。對于無毒單組元發(fā)動機(jī)工作過程的研究相對較少，其中主要是通過熱試車對無毒單組元發(fā)動機(jī)的比沖、起動特性等做一些宏觀的測試[14,15,16]，缺乏對單組元發(fā)動機(jī)工作的物理過程的理解。單組元發(fā)動機(jī)的仿真重點(diǎn)在于催化室內(nèi)的物理和化學(xué)過程。雖然針對催化床內(nèi)的流動和傳熱的仿真研究不少[17-21]，但是這些研究主要是針對一般工業(yè)中的催化床，無法直接將其方法和結(jié)論應(yīng)用于單組元火箭發(fā)動機(jī)。因此，本文通過對某一HAN基推進(jìn)劑單組元發(fā)動機(jī)進(jìn)行起動過程和穩(wěn)態(tài)工作過程的仿真研究，旨在探討適用于單組元發(fā)動機(jī)的仿真方法，并加強(qiáng)對單組元發(fā)動機(jī)內(nèi)部流動與傳熱過程的理解。

1 計(jì)算模型

發(fā)動機(jī)的催化床中堆積了數(shù)量眾多的催化劑顆粒，直接對顆粒間的流動和傳熱進(jìn)行仿真很困難，所以本文將催化床簡化為多孔介質(zhì)進(jìn)行計(jì)算，采用如下假設(shè)：

1) 前、后催化床均簡化為多孔介質(zhì)區(qū)域，忽略催化床前、后的多孔隔板；

2) 不考慮推進(jìn)劑的分解過程，認(rèn)為推進(jìn)劑在進(jìn)入催化床后完全分解；

3) 推進(jìn)劑從噴注器均勻流入催化床；

4) 只考慮催化床內(nèi)氣相的流動及氣相和固相間的傳熱；

5) 催化床的物理性質(zhì)各向同性；

6) 發(fā)動機(jī)外壁面只考慮輻射換熱。

本文應(yīng)用FLUENT軟件進(jìn)行計(jì)算。仿真求解的氣相方程為不帶化學(xué)反應(yīng)的多組分N-S方程。氣體成分為推進(jìn)劑催化反應(yīng)后的凍結(jié)組分，其性質(zhì)(如理論溫度、分子量、粘性系數(shù)等物性參數(shù))通過化學(xué)平衡計(jì)算給出。催化床內(nèi)的流動主要依據(jù)多孔介質(zhì)流動的Ergun方程[22]。多孔介質(zhì)內(nèi)的傳熱采用non-equilibrium模型，即氣相和固相的溫度不同，兩者之間存在換熱。催化室前床和后床的孔隙率、慣性阻力系數(shù)、粘性阻力系數(shù)和對流換熱系數(shù)根據(jù)催化床參數(shù)進(jìn)行計(jì)算。

計(jì)算采用基于壓力的分離算法；壓力修正采用SIMPLE方法；其他物理量的空間差分采用二階精度的迎風(fēng)格式。湍流模型采用標(biāo)準(zhǔn)k-ε模型。對于起動過程，時(shí)間推進(jìn)采用非定常的隱式二階精度格式進(jìn)行求解(雙時(shí)間法)。

1) 流體控制方程

在流體計(jì)算中，連續(xù)方程、能量方程和動量方程表示成統(tǒng)一的形式：

(1)

式中：φ為通用變量，Γ為廣義擴(kuò)散系數(shù)，S為廣義源項(xiàng)。

2) 多孔介質(zhì)的動量方程

多孔介質(zhì)對流體的作用通過動量源項(xiàng)在流體動量方程中體現(xiàn)，該源項(xiàng)包括粘性損失項(xiàng)和慣性損失項(xiàng)。第i個(gè)動量方程中的源項(xiàng)Si表示為：

(2)

方程右端第一項(xiàng)是粘性損失項(xiàng)，第二項(xiàng)是慣性損失項(xiàng)，D和C是給定的矩陣。對于均勻的多孔介質(zhì)，

(3)

(4)

式中：Dp為催化劑顆粒的直徑，ε為催化床的孔隙率。

3) 多孔介質(zhì)能量方程

多孔介質(zhì)模型對能量方程的影響體現(xiàn)在對時(shí)間導(dǎo)數(shù)和對流項(xiàng)的修正上，時(shí)間導(dǎo)數(shù)中，引入了固體區(qū)域?qū)Χ嗫捉橘|(zhì)的熱慣性效應(yīng)，能量方程如下：

(5)

4) 多孔介質(zhì)中流固耦合傳熱

多孔介質(zhì)中，氣相與固體的交界面處應(yīng)該滿足能量守恒，即

(6)

式中：λ為固體的導(dǎo)熱系數(shù)，h為流體的對流換熱系數(shù)(根據(jù)文獻(xiàn)[23]計(jì)算)，T為溫度，n為壁面的法向方向，下標(biāo)f和s分別表示流體和固體。

2 計(jì)算條件

2.1 發(fā)動機(jī)參數(shù)

計(jì)算對象為某HAN基推進(jìn)劑單組元試驗(yàn)發(fā)動機(jī)。試驗(yàn)測得該發(fā)動機(jī)的穩(wěn)態(tài)工作壓強(qiáng)(絕對壓力)Pc=0.72 MPa。

發(fā)動機(jī)推力室結(jié)構(gòu)如圖1所示。其中催化室由前床和后床組成；殼體材料為高溫合金GH1131。圖中T1、T2、T3點(diǎn)為試車時(shí)發(fā)動機(jī)外壁面的測溫點(diǎn)。

催化劑顆粒材料為含鋁氧化物。表1給出了催化劑顆粒的相關(guān)參數(shù)，從表1中計(jì)算出前床的孔隙率為ε=0.3，后床的孔隙率為ε=0.375。

2.2 計(jì)算網(wǎng)格及邊界條件

圖2所示為對發(fā)動機(jī)簡化后的二維軸對稱計(jì)算網(wǎng)格，整個(gè)計(jì)算域網(wǎng)格總數(shù)為18 758個(gè)單元，其中流場計(jì)算域的網(wǎng)格為16 790個(gè)單元，固體計(jì)算域的網(wǎng)格數(shù)1 968個(gè)單元。

表1 前、后床顆粒的物性參數(shù)

推進(jìn)劑入口邊界條件采用流量入口。根據(jù)試驗(yàn)數(shù)據(jù)，擬合得到流量與時(shí)間之間的關(guān)系式：

qm=100.97t-6 605.90t2+189 982.62t3-

1 491 000t4+2 650 870t5-0.31 g/s

時(shí)間t的值為0～0.154 s，大于0.154 s后發(fā)動機(jī)的流量為30 g/s。根據(jù)上式編寫UDF作為流量進(jìn)口的邊界條件。出口邊界條件為壓力出口，設(shè)置為試驗(yàn)時(shí)的環(huán)境壓強(qiáng)值(100 Pa)。發(fā)動機(jī)的外壁面為輻射邊界條件，輻射系數(shù)取0.5，環(huán)境溫度假定為300 K。

在時(shí)間推進(jìn)的雙時(shí)間法中，時(shí)間步長和時(shí)間步長內(nèi)迭代的次數(shù)影響計(jì)算穩(wěn)定性和準(zhǔn)確性；經(jīng)過試算驗(yàn)證，本文中時(shí)間步長取為1×10-6s，內(nèi)迭代次數(shù)取為30次，保證每個(gè)時(shí)間步長內(nèi)殘差基本收斂到穩(wěn)定值。

2.3 熱力計(jì)算

因?yàn)镠AN推進(jìn)劑的催化反應(yīng)十分的復(fù)雜，故本文計(jì)算中假設(shè)推進(jìn)劑分解后的氣體為化學(xué)平衡組分，采用熱力計(jì)算的方法對燃?xì)獬煞趾桶l(fā)動機(jī)理論性能進(jìn)行計(jì)算。

熱力計(jì)算時(shí)，推進(jìn)劑的配方依據(jù)發(fā)動機(jī)性能和參考文獻(xiàn)[24]，假定為63.2%硝酸羥胺(HAN)，20.0%三乙醇胺硝酸鹽(TEAN),16.8%水(H2O)。

表2列出了根據(jù)化學(xué)平衡計(jì)算結(jié)果得到的燃?xì)饨M分?jǐn)?shù)據(jù)(略去摩爾含量小于1%的組分)：

表2 燃?xì)獬煞?/p>

通過熱力計(jì)算得到了燃?xì)馄骄肿恿繛?9.382，氣體常數(shù)為428.979 J/(mol·K)，動力黏性系數(shù)為5.397×10-5Pa·s，導(dǎo)熱系數(shù)為0.209 1 W/(m·K)；燃?xì)鉁囟葹? 587 K，發(fā)動機(jī)的理論比沖為228.1 s。

3 計(jì)算結(jié)果討論

本文計(jì)算了發(fā)動機(jī)的起動過程和穩(wěn)態(tài)工作過程，并將計(jì)算結(jié)果與試驗(yàn)測量結(jié)果進(jìn)行了對比。在發(fā)動機(jī)試車前，在發(fā)動機(jī)外壁的不同位置布置3個(gè)熱電偶用于測量外壁面溫度變化(熱電偶的位置如圖1所示)，同時(shí)在電磁閥門上游及催化床下游布置壓力傳感器用于測量壓力變化。

3.1 起動過程

在HAN基推進(jìn)劑發(fā)動機(jī)試車時(shí)，需要對催化床進(jìn)行預(yù)熱，以使推進(jìn)劑能夠順利分解，催化床的有效預(yù)熱溫度為120 ℃。本文計(jì)算中，將計(jì)算區(qū)域的初始溫度統(tǒng)一設(shè)置為120 ℃。

圖3反映了起動過程中推力室對稱面上的溫度變化，圖4為t=2 s時(shí)其內(nèi)部燃?xì)夂痛呋驳妮S向溫度曲線。

觀察這兩圖可知，催化床的升溫過程包括兩個(gè)方面：一是催化床整體的溫度在高溫?zé)峤鈿怏w的傳熱作用下逐漸升高；二是高溫從催化床上游逐漸向下游傳遞。從溫度分布上看，催化床上游溫度高于下游，另外催化床內(nèi)部靠近軸線位置的溫度也高于壁面附近溫度，即在催化床內(nèi)不僅存在軸向溫度梯度，也存在徑向溫度梯度。此外，由于推力室殼體向外界的輻射傳熱，使得殼體溫度顯著低于催化室溫度。

需要指出，因?yàn)楸疚牡挠?jì)算模型沒有考慮液體推進(jìn)劑分解成為高溫燃?xì)獾倪^程，所以計(jì)算得到的溫度在噴注器處呈現(xiàn)最大值，這與實(shí)際情況并不相符。實(shí)際狀態(tài)應(yīng)該是噴注器處的溫度等于或略高于液體推進(jìn)劑溫度，而后在下游液相分解反應(yīng)時(shí)迅速升溫。

圖5給出了某次試驗(yàn)的發(fā)動機(jī)外壁溫度變化，所示實(shí)驗(yàn)值為圖1中T1處的溫度值。從試驗(yàn)曲線看到，因?yàn)闊崃總鬟f到壁面處需要一定時(shí)間，所以0～1 s的溫度變化很??；當(dāng)壁面溫度開始上升時(shí)，其升溫速率基本呈線性。本文計(jì)算了外壁面輻射系數(shù)為不同值時(shí)的外壁面升溫曲線，對比試驗(yàn)數(shù)據(jù)可知當(dāng)熱輻射系數(shù)設(shè)置為0.4時(shí)，計(jì)算得到的升溫速率和試驗(yàn)值接近。

圖6給出了起動過程中推力室的溫度變化過程(測點(diǎn)位于催化床下游的燃燒室內(nèi)，距內(nèi)壁面0.5 mm)。從圖6可知：(1)起動時(shí)燃?xì)饬鬟^催化床的時(shí)間大約為0.094 s；(2)高溫燃?xì)庠诖呋矁?nèi)流動時(shí)向催化床顆粒傳熱，流出催化床時(shí)溫度較低，起動初期小于200 ℃，而后隨著時(shí)間推移逐漸升高，說明催化床的溫度也在升高，燃?xì)馀c催化床之間的換熱量逐漸減小。

圖7給出了起動過程中推力室的壓力Pc變化過程(測點(diǎn)位于催化床下游，距內(nèi)壁面0.5 mm，與圖5的溫度測點(diǎn)相同)?？梢姡?jì)算得到的升壓過程曲線與試驗(yàn)曲線較為吻合，但計(jì)算的穩(wěn)態(tài)值略高于試驗(yàn)值。另外，升壓過程存在比較明顯的拐點(diǎn)(圖7中示出了拐點(diǎn)坐標(biāo))，在拐點(diǎn)前壓力上升非?？欤髣t呈現(xiàn)緩慢爬升。升壓過程的這種表現(xiàn)具有單組元發(fā)動機(jī)的典型特點(diǎn)，與參考文獻(xiàn)[25]描述的肼發(fā)動機(jī)基本一致：即以壓力曲線的拐點(diǎn)為界可將升壓過程分為兩個(gè)階段；前面升壓快的部分對應(yīng)著燃?xì)獬涮钔屏κ业倪^程，燃?xì)庠谶@個(gè)階段中迅速充滿推力室并使噴管達(dá)到壅塞流動狀態(tài)，但是因?yàn)槿細(xì)庀虼呋矀鳠幔谷細(xì)鉄崮軗p失較大，導(dǎo)致拐點(diǎn)壓力并不高；而后面升壓慢的部分則對應(yīng)著催化床的升溫過程，在這個(gè)階段中隨著燃?xì)鈱Υ呋驳膫鳠?，使催化床溫度逐漸升高，進(jìn)而燃?xì)獾臒釗p失逐漸減小，壓力進(jìn)一步升高。因?yàn)榇呋w粒的熱容較大，所以催化床溫度上升較慢，升壓過程的第二階段時(shí)間較長。此外，在升壓曲線中沒有起動峰值出現(xiàn)，說明催化床的催化效能較好，沒有產(chǎn)生推進(jìn)劑積存現(xiàn)象。

從計(jì)算結(jié)果來看，升壓曲線拐點(diǎn)的出現(xiàn)比試驗(yàn)值早了約0.11 s，這是因?yàn)橛?jì)算沒有考慮推進(jìn)劑的分解過程。計(jì)算的拐點(diǎn)壓力值與試驗(yàn)值非常接近，且穩(wěn)態(tài)壓力的計(jì)算值與試驗(yàn)值也比較吻合，這說明采用化學(xué)平衡計(jì)算得到的燃?xì)鈪?shù)與實(shí)際燃?xì)獗容^符合。

3.2 穩(wěn)態(tài)工作過程

發(fā)動機(jī)穩(wěn)定工作時(shí)，流場中的參數(shù)基本不再發(fā)生變化。圖8和圖9給出了推力室對稱面上的壓力分布云圖和沿著對稱軸的壓力分布。由于催化床為多孔介質(zhì)，其孔隙率影響床內(nèi)壓降，因前床的孔隙率較小所以壓降較大，后床的孔隙率較大所以壓降較小，這在壓力分布圖中得到了較好的體現(xiàn)。另外，從計(jì)算可得催化床下游的室壓計(jì)算值約為0.77 MPa，高于試驗(yàn)值(0.72 MPa)約7%。計(jì)算值高于試驗(yàn)值的原因，主要由于計(jì)算中沒有考慮反應(yīng)效率，認(rèn)為催化反應(yīng)處于平衡狀態(tài)；另一個(gè)原因是催化床的壓降小于實(shí)際壓降，這就需要調(diào)整流阻系數(shù)。

圖10給出了溫度分布云圖，從圖10可以看出推力室的最高溫度約為1 580 K；收斂段部位的燃?xì)鉁囟葹? 300～1 400 K；從推力室下游到噴管出口溫度一直下降。由于發(fā)動機(jī)殼體向外熱輻射，所以殼體溫度低于燃?xì)鉁囟取４送?，殼體在喉部的溫度高于其他部位，這一方面是因?yàn)楹聿繗んw與燃?xì)獾膶α鲹Q熱系數(shù)大，另一方面也是喉部表面積小，向環(huán)境輻射的熱量少。

圖11給出了軸向速度的幅值分布，通過對噴管出口的動量和壓差進(jìn)行積分得到計(jì)算的比沖為224 s，與理論計(jì)算值(228.1 s)接近。

4 結(jié)論

1) 采用本文建立的計(jì)算模型和模型參數(shù)計(jì)算方法，可以對單組元液體發(fā)動機(jī)進(jìn)行工作過程仿真，仿真得到的升壓曲線和外壁溫度與試驗(yàn)結(jié)果符合較好。

2) 催化床承載了較大的壓降，且壓降受孔隙率影響很大；前催化床孔隙率略小，所以前催化床壓降明顯大于后催化床壓降。

3) 單組元發(fā)動機(jī)起動過程中，室壓上升曲線分為兩個(gè)階段，第二階段時(shí)間長于第一階段。

4) 當(dāng)計(jì)算模型中不考慮液體推進(jìn)劑的反應(yīng)過程時(shí)，計(jì)算得到的初期升壓速率快于實(shí)際情況。

[1] 周漢申.單組元液體火箭發(fā)動機(jī)設(shè)計(jì)與研究[M].中國宇航出版社,2009.

[2] 白云峰，林慶國，金盛宇，等.過氧化氫單元催化分解火箭發(fā)動機(jī)研究[J].火箭推進(jìn),2006,32(4):15-20.

[3] 宋長青,徐萬武,張家奇,陳健.氧化亞氮推進(jìn)技術(shù)研究進(jìn)展[J].火箭推進(jìn),2014,40(2).

[4] 禹天福，賈月.一氧化二氮催化分解的研究與應(yīng)用[J].化學(xué)推進(jìn)劑與高分子材料,2006,4(2):6-10.

[5] 周悅,公緒濱,方濤.硝酸羥銨基無毒單組元推進(jìn)劑應(yīng)用探討[J].導(dǎo)彈與航天運(yùn)載技術(shù),2015(4).

[6] SHUTARO NISHIKIZAWA,TAKEHIRO OHIRA,HIRONORI SAHARA.Development of Mono /Bi-Propellant Propulsion System for Microsatellites and Planning of Space Demonstration[C].49th AIAA/ASME/SAE/ASEE Joint Propulsion Conference,San Jose,CA,2013.

[7] AMROUSSE R,HORI K,FETIMI W,et al.HAN and ADN as liquid ionic monopropellants:Thermal and catalytic decomposition processes[J].Applied Catalysis B:Environmental,2012,127:121-128.

[8] AMROUSSE R,KATSUMI T,AZUMA N,et al.Hydroxylammonium nitrate (HAN)-based green propellant as alternative energy resource for potential hydrazine substitution:From lab scale to pilot plant scale-up[J].Combustion and Flame 2017,176:334-348.

[9] KATSUMI T,INOUE T,NAKATSUKA J,et al.HAN-Based Green Propellant,Application,and Its Combustion Mechanism[J].Combustion,Explosion,and Shock Waves 2012,48(5):536-543.

[10] CHANG HWAN HWANG，SEUNG WOOK BAEK，SUNG JUNE CHO.Experimental investigation of decomposition and evaporation characteristics of HAN-based monopropellants[J].Combustion and Flame,2014,161:1109-1116.

[11] 陳永康,陳也弘,安振濤,等.HAN-基凝膠推進(jìn)劑的熱分解反應(yīng)動力學(xué)[J].火炸藥學(xué)報(bào),2016,39(4):77-81.

[12] 潘玉竹.HAN基液體推進(jìn)劑高壓燃燒特性的實(shí)驗(yàn)研究和數(shù)值模擬[D].南京:南京理工大學(xué),2013,4.

[13] 任曉光.硝酸羥胺基推進(jìn)劑的催化分解研究[D].大連:大連輕工業(yè)學(xué)院,2006.

[14] CHEN Jun,LI Guoxiu,ZHANG Tao,et al.Experimental investigation of the catalytic decomposition and combustion characteristics of a non-toxic ammonium dinitramide (ADN)-based monopropellant thruster.Acta Astronautica,2016，129:367-373.

[15] ZHANG Tao,LI Guoxiu,YU Yusong,et al.Numerical simulation of ammonium dinitramide (ADN)-based non-toxic aerospace propellant decomposition and combustion in a monopropellant thruster[J].Energy Conversion and Management,2014,(87):965-974.

[16] 李華樂.ADN/甲醇推進(jìn)劑燃燒反應(yīng)動力學(xué)模型及推力器燃燒過程仿真的研究[D].北京:北京交通大學(xué),2014.

[17] GUARDO A.,COUSSIRAT M.,RECASENS F,et al.CFD studies on particle-to-fluid mass and heat transfer in packed beds:Free convection effects in supercritical fluids[J].Chemical Engineering Science，2007,62：5503-5511.

[18] GUARDO A,COUSSIRAT M,RECASENS F,et al.CFD study on particle-to-fluid heat transfer in fixed bed reactors:Convective heat transfer at lowand high pressure[J].Chemical Engineering Science，2006,61:4341- 4353.

[19] 郭雪巖.CFD方法在固定床反應(yīng)器傳熱研究中的應(yīng)用[J].化工學(xué)報(bào),2008,59(8):1914-1922.

[20] KAMIUTO K,SAITOH S.Simultaneous Heat and Mass Transfer in Packed Bed Catalytic Reactors[J].Journal of Thermophysics and Heat Transfer,1995,9(3):524-530.

[21] ZHOU X,HITT D L.Numerical Modeling of Monopropellant Decomposition in a Micro-Catalyst Bed[C].35th AIAA Fluid Dynamics Conference and Exhibit,2005(6):1-10.

[22] ERGUN S.Fluid flow through packed columns [J].Chemical Engineering Progress,1952,48:89-94.

[23] 鄭坤燦,侯戎彬,王景甫,等.基于孔喉模型的多孔介質(zhì)對流換熱系數(shù)的分形研究[J].科學(xué)技術(shù)與工程,2014,14(28):44-49.

[24] 王宏偉,王建偉.AF-315液體單元推進(jìn)劑研究進(jìn)展[J].化學(xué)推進(jìn)劑與高分子材料,2010,8(5):6-9.

[25] SCHMITZ B W,SMETH W W.Development of Design and Scaling Criteria for Monopropellant Hydrazine Reactors Employing Shell 405 Spontaneous Catalyst[R].RRC Report 66-R-76,1967,11(1).