硝酸羥胺基單組元發動機起動過程數值模擬

孫得川，金東洙，于澤游

(大連理工大學航空航天學院，工業裝備結構分析國家重點實驗室， 遼寧 大連 116024)

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

1 計算模型

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

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

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

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

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

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

6) 發動機外壁面只考慮輻射換熱。

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

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

1) 流體控制方程

在流體計算中，連續方程、能量方程和動量方程表示成統一的形式：

(1)

式中：φ為通用變量，Γ為廣義擴散系數，S為廣義源項。

2) 多孔介質的動量方程

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

(2)

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

(3)

(4)

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

3) 多孔介質能量方程

多孔介質模型對能量方程的影響體現在對時間導數和對流項的修正上，時間導數中，引入了固體區域對多孔介質的熱慣性效應，能量方程如下：

(5)

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

多孔介質中，氣相與固體的交界面處應該滿足能量守恒，即

(6)

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

2 計算條件

2.1 發動機參數

計算對象為某HAN基推進劑單組元試驗發動機。試驗測得該發動機的穩態工作壓強(絕對壓力)Pc=0.72 MPa。

發動機推力室結構如圖1所示。其中催化室由前床和后床組成；殼體材料為高溫合金GH1131。圖中T1、T2、T3點為試車時發動機外壁面的測溫點。

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

2.2 計算網格及邊界條件

圖2所示為對發動機簡化后的二維軸對稱計算網格，整個計算域網格總數為18 758個單元，其中流場計算域的網格為16 790個單元，固體計算域的網格數1 968個單元。

表1 前、后床顆粒的物性參數

推進劑入口邊界條件采用流量入口。根據試驗數據，擬合得到流量與時間之間的關系式：

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

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

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

在時間推進的雙時間法中，時間步長和時間步長內迭代的次數影響計算穩定性和準確性；經過試算驗證，本文中時間步長取為1×10-6s，內迭代次數取為30次，保證每個時間步長內殘差基本收斂到穩定值。

2.3 熱力計算

因為HAN推進劑的催化反應十分的復雜，故本文計算中假設推進劑分解后的氣體為化學平衡組分，采用熱力計算的方法對燃氣成分和發動機理論性能進行計算。

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

表2列出了根據化學平衡計算結果得到的燃氣組分數據(略去摩爾含量小于1%的組分)：

表2 燃氣成分

通過熱力計算得到了燃氣平均分子量為19.382，氣體常數為428.979 J/(mol·K)，動力黏性系數為5.397×10-5Pa·s，導熱系數為0.209 1 W/(m·K)；燃氣溫度為1 587 K，發動機的理論比沖為228.1 s。

3 計算結果討論

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

3.1 起動過程

在HAN基推進劑發動機試車時，需要對催化床進行預熱，以使推進劑能夠順利分解，催化床的有效預熱溫度為120 ℃。本文計算中，將計算區域的初始溫度統一設置為120 ℃。

圖3反映了起動過程中推力室對稱面上的溫度變化，圖4為t=2 s時其內部燃氣和催化床的軸向溫度曲線。

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

需要指出，因為本文的計算模型沒有考慮液體推進劑分解成為高溫燃氣的過程，所以計算得到的溫度在噴注器處呈現最大值，這與實際情況并不相符。實際狀態應該是噴注器處的溫度等于或略高于液體推進劑溫度，而后在下游液相分解反應時迅速升溫。

圖5給出了某次試驗的發動機外壁溫度變化，所示實驗值為圖1中T1處的溫度值。從試驗曲線看到，因為熱量傳遞到壁面處需要一定時間，所以0～1 s的溫度變化很小；當壁面溫度開始上升時，其升溫速率基本呈線性。本文計算了外壁面輻射系數為不同值時的外壁面升溫曲線，對比試驗數據可知當熱輻射系數設置為0.4時，計算得到的升溫速率和試驗值接近。

圖6給出了起動過程中推力室的溫度變化過程(測點位于催化床下游的燃燒室內，距內壁面0.5 mm)。從圖6可知：(1)起動時燃氣流過催化床的時間大約為0.094 s；(2)高溫燃氣在催化床內流動時向催化床顆粒傳熱，流出催化床時溫度較低，起動初期小于200 ℃，而后隨著時間推移逐漸升高，說明催化床的溫度也在升高，燃氣與催化床之間的換熱量逐漸減小。

圖7給出了起動過程中推力室的壓力Pc變化過程(測點位于催化床下游，距內壁面0.5 mm，與圖5的溫度測點相同)。可見，計算得到的升壓過程曲線與試驗曲線較為吻合，但計算的穩態值略高于試驗值。另外，升壓過程存在比較明顯的拐點(圖7中示出了拐點坐標)，在拐點前壓力上升非常快，而后則呈現緩慢爬升。升壓過程的這種表現具有單組元發動機的典型特點，與參考文獻[25]描述的肼發動機基本一致：即以壓力曲線的拐點為界可將升壓過程分為兩個階段；前面升壓快的部分對應著燃氣充填推力室的過程，燃氣在這個階段中迅速充滿推力室并使噴管達到壅塞流動狀態，但是因為燃氣向催化床傳熱，使燃氣熱能損失較大，導致拐點壓力并不高；而后面升壓慢的部分則對應著催化床的升溫過程，在這個階段中隨著燃氣對催化床的傳熱，使催化床溫度逐漸升高，進而燃氣的熱損失逐漸減小，壓力進一步升高。因為催化顆粒的熱容較大，所以催化床溫度上升較慢，升壓過程的第二階段時間較長。此外，在升壓曲線中沒有起動峰值出現，說明催化床的催化效能較好，沒有產生推進劑積存現象。

從計算結果來看，升壓曲線拐點的出現比試驗值早了約0.11 s，這是因為計算沒有考慮推進劑的分解過程。計算的拐點壓力值與試驗值非常接近，且穩態壓力的計算值與試驗值也比較吻合，這說明采用化學平衡計算得到的燃氣參數與實際燃氣比較符合。

3.2 穩態工作過程

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

圖10給出了溫度分布云圖，從圖10可以看出推力室的最高溫度約為1 580 K；收斂段部位的燃氣溫度為1 300～1 400 K；從推力室下游到噴管出口溫度一直下降。由于發動機殼體向外熱輻射，所以殼體溫度低于燃氣溫度。此外，殼體在喉部的溫度高于其他部位，這一方面是因為喉部殼體與燃氣的對流換熱系數大，另一方面也是喉部表面積小，向環境輻射的熱量少。

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

4 結論

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

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

3) 單組元發動機起動過程中，室壓上升曲線分為兩個階段，第二階段時間長于第一階段。

4) 當計算模型中不考慮液體推進劑的反應過程時，計算得到的初期升壓速率快于實際情況。

[1] 周漢申.單組元液體火箭發動機設計與研究[M].中國宇航出版社,2009.

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

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

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

[5] 周悅,公緒濱,方濤.硝酸羥銨基無毒單組元推進劑應用探討[J].導彈與航天運載技術,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].火炸藥學報,2016,39(4):77-81.

[12] 潘玉竹.HAN基液體推進劑高壓燃燒特性的實驗研究和數值模擬[D].南京:南京理工大學,2013,4.

[13] 任曉光.硝酸羥胺基推進劑的催化分解研究[D].大連:大連輕工業學院,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/甲醇推進劑燃燒反應動力學模型及推力器燃燒過程仿真的研究[D].北京:北京交通大學,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方法在固定床反應器傳熱研究中的應用[J].化工學報,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] 鄭坤燦,侯戎彬,王景甫,等.基于孔喉模型的多孔介質對流換熱系數的分形研究[J].科學技術與工程,2014,14(28):44-49.

[24] 王宏偉,王建偉.AF-315液體單元推進劑研究進展[J].化學推進劑與高分子材料,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).