趙 瑜，汪浩平，何 快，張建欣

（1.上海航天動力技術(shù)研究所，上海 201109；2.上海航天技術(shù)研究院，上海 201109）

0 引言

固液發(fā)動機是介于液體火箭發(fā)動機與固體火箭發(fā)動機間的一種推進系統(tǒng)，由氧化劑供應(yīng)系統(tǒng)和固體燃料系統(tǒng)組成，兩者單獨存貯，相互隔離，因此不僅具有很高的安全性和可靠性，而且通過控制氧化劑的供應(yīng)，可實現(xiàn)發(fā)動機的多次起動與關(guān)機及推力的調(diào)節(jié)［1－2］。固液發(fā)動機被視為火箭推進系統(tǒng)的一個重要發(fā)展方向，成為液體和固體火箭發(fā)動機的有力競爭者，特別是其中的N2O／HTPB體系固液發(fā)動機已在工程中獲得應(yīng)用并顯示出優(yōu)異的推進性能。該類型發(fā)動機的氧化劑無毒、無污染，可長期自然溫度貯存；飽和蒸汽壓高，便于實現(xiàn)氧化劑供應(yīng)和簡化管路系統(tǒng)設(shè)計，其具有的自持續(xù)分解放熱效應(yīng)為實現(xiàn)固液發(fā)動機多次重復(fù)起動提供了便利。此外，N2O／HTPB體系固液發(fā)動機的氧燃比適用范圍寬，氧燃比變化對發(fā)動機總體性能影響較小，這為燃燒過程組織提供了更大的自由度。

與固體和液體發(fā)動機相比，固液發(fā)動機優(yōu)勢實現(xiàn)的關(guān)鍵是實現(xiàn)高效燃燒和高比沖，提高系統(tǒng)裝填系數(shù)，降低系統(tǒng)殘藥量。這些問題的核心是高效燃燒組織，即將氧化劑的霧化、氣化、分解與固體燃料的氣化、分解進行有效的組織，實現(xiàn)兩者的高效摻混。因此，目前國內(nèi)外在固液發(fā)動機研究中燃燒和流動仍是其中的核心。但在發(fā)動機試車過程中能測得的數(shù)據(jù)相當(dāng)有限，且對內(nèi)部的燃燒和流場結(jié)構(gòu)仍無有效手段可直接觀測。因此，數(shù)值仿真成為普遍采用以及更傾向的一種技術(shù)途徑。隨著近年來仿真技術(shù)的發(fā)展，它越來越成為預(yù)測發(fā)動機性能、優(yōu)化發(fā)動機總體設(shè)計的可靠工具，尤其是在新方案的應(yīng)用和新技術(shù)的發(fā)展中更是必不可少。本文對某N2O／HTPB體系固液發(fā)動機工作過程的仿真進行了研究。

1 發(fā)動機工作過程數(shù)值模型

固液發(fā)動機的工作機理如圖1所示。由圖可知：發(fā)動機中燃料與氧化劑在物理與相態(tài)上彼此分離，其燃燒過程由兩者的宏觀擴散過程控制，即經(jīng)過噴注器后形成的氣態(tài)氧化劑與固體燃料高溫?zé)峤夂笮纬傻臍鈶B(tài)燃料在宏觀擴散作用下，于邊界層內(nèi)相遇而形成有效燃燒區(qū)并不斷向氣相和固相傳遞熱量，從而維持燃燒的持續(xù)進行。因此，為模擬固液發(fā)動機的工作過程需對其中的氧化劑噴注和霧化、氣相燃燒、燃料注入，以及流動過程等進行數(shù)學(xué)模型的建立和選擇。但在此之前，計算所用控制方程需提前確定。此時，鑒于固液發(fā)動機內(nèi)流場具有黏性和可壓縮的特性，并考慮對燃燒過程的引入，控制方程選取帶組分的N－S方程。方程中各組分質(zhì)量生成率的求解是數(shù)值模擬要解決的關(guān)鍵之一，且其不同計算方法就代表了不同的燃燒模型。

圖1 固液發(fā)動機工作機理Fig.1 Operating principle of hybrid rocket motor

1.1 氣相化學(xué)反應(yīng)動力學(xué)模型

研究表明：固體燃料HTPB在受熱時不發(fā)生固相反應(yīng)，而是產(chǎn)生高溫分解［3－4］。因收集到的氣相產(chǎn)物大部分為C4H6，故通常假設(shè)C4H6為HTPB的高溫分解產(chǎn)物。此時，N2O／HTPB體系的氣相反應(yīng)物應(yīng)為C4H6，N2O。它們在7MPa下進行絕熱燃燒，分析后可知當(dāng)氧燃比（O／F）為7附近時其燃溫最高（約3 542K），此時的主要反應(yīng)產(chǎn)物（摩爾分數(shù)在10－4以上的組分）見表1。

根據(jù)表1，經(jīng)原子配平后可得在O／F＝7時C4H6與N2O的反應(yīng)方程為

由于不能確定反應(yīng)式（1）的化學(xué)動力學(xué)參數(shù)，如活化能和指前因子等，在應(yīng)用該式對N2O／HTPB體系固液發(fā)動機的氣相燃燒進行數(shù)值模擬時，化學(xué)反應(yīng)速率將采用渦耗散模型（Eddy－Dissipation Modal）計算，因這在確定化學(xué)反應(yīng)速率時無需化學(xué)動力學(xué)參數(shù)。具體來說，在渦耗散模型中當(dāng)湍流出現(xiàn)，即κ／ε＞0時燃燒就開始進行。此時，化學(xué)反應(yīng)速率取決于大渦摻混時間尺寸κ／ε，即反應(yīng)方程r內(nèi)組分i的生成率Ri，r取

中較小的一個。此處：v′i，r為反應(yīng)式r中反應(yīng)物i的化學(xué)當(dāng)量系數(shù)；v″j，r為反應(yīng)式r中生成物j的化學(xué)當(dāng)量系數(shù)；Mw，i為組分i的分子量；YP為任意燃燒產(chǎn)物P的質(zhì)量分數(shù)；YR為特定反應(yīng)物R的質(zhì)量分數(shù)；A為經(jīng)驗常數(shù)，取4.0；B為經(jīng)驗常數(shù)，取B＝0.5。

1.2 固體燃料注入方式

固液發(fā)動機工作時，氣態(tài)燃料（N2O／HTPB體系中為C4H6）由燃面噴出。但目前通用軟件中尚無對應(yīng)的、可直接應(yīng)用的此類邊界條件。因此，該邊界條件在多數(shù)數(shù)值模擬中被處理成壁面＋質(zhì)量源項的方式［5－6］。為與真實情況相符，質(zhì)量源項的添加位置應(yīng)為圖2所示緊貼壁面的一個薄層中。需說明的是，三維模型中固體燃料的型面較復(fù)雜，在處理加質(zhì)薄層時很難像二維模型那樣方便進行單獨的CAD建模。因此，面對三維模型本次數(shù)值模擬中用UDF方式先標記所有緊貼燃面的單元，然后以體積加權(quán)的方法將所需添加的質(zhì)量源項依次分配到各單元內(nèi)。編寫的UDF中主要采用DEFINE＿ADJUST函數(shù)，DEFINE＿SOURCE函數(shù)實現(xiàn)。同時，為節(jié)省計算時間，該UDF是以并行的方式編寫和運行。

表1 C4H6，N2O反應(yīng)的生成物Tab.1 Product of reaction between C4H6and N2O

圖2 燃料質(zhì)量源項添加Fig.2 Demonstration of fuel source’s addition

1.3 氧化劑（N2O）噴注系統(tǒng)模型

與液體火箭發(fā)動機類似，固液混合火箭發(fā)動機中的氧化劑噴注系統(tǒng)多樣，本文對圓孔噴注和自擊霧化方式進行了研究。在數(shù)值模擬中用實心錐噴注模型，并輔以碰撞、二級破碎和蒸發(fā)模型描述噴注系統(tǒng)的工作過程。模型中所需的N2O飽和蒸汽壓可表示為

式中：A＝9.670 86；B＝1 429.91K；C＝－44K；pscal＝0.1MPa。

或者，根據(jù)文獻［7］中的試驗數(shù)據(jù)對N2O的飽和蒸汽壓進行設(shè)置。此外，模型中所需其他參數(shù)為：分子量44kg／kmol；密度780kg／m3；定壓比熱容2.269 5×103J／（kg·K）；汽化潛熱1.86×105J／kg；黏性系數(shù)9×10－5Pa·s；導(dǎo)熱系數(shù)（參考水）0.6W／（m·K）；汽化溫度273K；沸點315.5K；擴散 系 數(shù) （參 考 水 ）3.05×10－5m2／s；表 面 張 力0.005n／m。其中部分參數(shù)無確鑿可資利用的數(shù)據(jù)，故參考了N2O其他狀態(tài)下的數(shù)值或采用了相近物質(zhì)的數(shù)據(jù)予以替代。

2 計算結(jié)果與分析

用本文建立的模型，對所設(shè)計固液發(fā)動機的工作過程進行數(shù)值模擬。因發(fā)動機的氧化劑噴注面板和燃燒室（包括藥柱）為周期旋轉(zhuǎn)對稱結(jié)構(gòu)，計算時只選取了1／4結(jié)構(gòu)，并對其進行了網(wǎng)格劃分（單元總數(shù)約250萬）和邊界條件設(shè)置（如圖3所示）。

圖3 三維發(fā)動機模型邊界條件設(shè)置Fig.3 Boundary conditions of 3－D model

計算時，設(shè)氧化劑總流量1kg／s，則模型中每個氧化劑噴注孔（共64個孔）的質(zhì)量流量0.015 6kg／s，所得氧化劑的霧化效果如圖4所示。由圖4可知：N2O液滴經(jīng)噴注孔噴出后確實形成了自擊霧化，與初始直徑1mm相比自擊后其直徑降低了至少2個量級，由此證實設(shè)計噴注器的霧化效果較好。

圖4 N2O液滴直徑分布Fig.4 Distribution of N2O droplet’s diameter

考察發(fā)動機工作時燃燒室的各主要參數(shù)。模型中N2O與C4H6間的氧燃比取最佳氧燃比7∶1，所得燃燒室壓強、溫度和平均分子量分別如圖5～7所示。由圖可知：燃燒室的壓強和溫度約為6.6MPa，3 600K，該數(shù)值與熱力學(xué)計算得到的理論值相近。將仿真所得的燃燒室溫度3 600K和噴管喉部處的平均分子量26.4g／mol代入零維內(nèi)彈道公式

在已給定的流量和噴管喉部直徑（18.5mm）下所得燃燒室壓強6.98MPa，與仿真得到的6.6MPa較接近。此處：為質(zhì)量流率；R為氣體常數(shù)；pc為燃燒室壓力；Tf為燃燒室溫度；At為發(fā)動機喉部面積。其中：k為比熱比。因此，由熱力學(xué)計算和內(nèi)彈道分析可知，所建模型以及使用的仿真方法較可靠和精確。此外，由氣態(tài)N2O的濃度分布（圖8）可知：N2O在發(fā)動機尾噴管出口處的濃度已基本接近零，表明氧化劑和燃料在燃燒室內(nèi)的燃燒較充分，由此證實設(shè)計中采用前預(yù)混室，擾流環(huán)和補燃室提高燃燒效率的方法可行且有效。

圖5 三維固液混合發(fā)動機工作過程模擬壓強分布Fig.5 Pressure simulation result of hybrid rocket motor’s operating process

圖6 三維固液混合發(fā)動機工作過程模擬溫度分布Fig.6 Temperature simulation result of hybrid rocket motor’s operating process

圖7 三維固液混合發(fā)動機工作過程模擬平均分子量分布Fig.7 Average molecular weight simulation result of hybrid rocket motor’s operating process

圖8 燃燒室內(nèi)N2O質(zhì)量分數(shù)分布Fig.8 Distribution of N2O’s mass fraction in combustion chamber

需說明的是，由圖6可知，與傳統(tǒng)固體發(fā)動機相比固液發(fā)動機工作時溫度場極不均勻。其原因一是固液發(fā)動機工作時其燃燒為擴散燃燒，氧化劑和燃料發(fā)生化學(xué)反應(yīng)前需要一個摻混過程，二是因固液發(fā)動機的藥型和擾流結(jié)構(gòu)將流場變得復(fù)雜，極大地影響了摻混過程。可預(yù)見該不均勻現(xiàn)象在三維效應(yīng)下將更突出，而這也體現(xiàn)了固液發(fā)動機工作過程仿真采用三維模型的必要性。

3 結(jié)束語

根據(jù)對N2O／HTPB體系固液發(fā)動機工作時涉及的各物理過程建模，本文對發(fā)動機燃燒室內(nèi)的流動、燃燒過程進行了三維數(shù)值模擬。由計算結(jié)果可知，所建模型和方法可靠且較精確，可用于預(yù)示固液發(fā)動機的內(nèi)彈道和性能。

［1］ 廖少英.固液火箭發(fā)動機性能特征［J］.上海航天，2004，21（5）：8－14.

［2］ 楊玉新，胡春波，何國強，等.固液混合火箭發(fā)動機中的關(guān)鍵技術(shù)及其發(fā)展［J］.宇航學(xué)報，2008，29（5）：1616－1621.

［3］ CHIAVERINI M J，KUO K K.Fundamentals of hybrid rocket combustion and propulsion［M］.Virginia：American Institute of Aeronautics and Astronautics Inc，2007.

［4］ CHIAVERINI M J，HARTING G C，LU Yeucherng，et al.Pyrolysis behavior of hybrid rocket solid fuels under rapid heating conditions［J］.Journal of Propulsion and Power，1999，15（6）：888－895.

［5］ LI Xin－tian，TIAN Hui，CAI Guo－biao.Numerical analysis of fuel regression rate distribution characteristics in hybrid rocket motors with different fuel types［J］.Science China，2013（7）：1807－1817.

［6］ BELLOMO N，LAZZARIN M，BARATO F.Numerical Investigation of the effect of a diaphragm on the performance of a hybrid rocket motor［C］／／46th AIAA／ASME／SAE／ASEE Joint Propulsion Conference ＆Exhibit.Nashville：AIAA，2010：1－16.

［7］ JAMIESON D T，BEATON C F，AMBROSE D，et al.Thermophysical properties of nitrous oxide［R］.Engineering Science Data Unit，1991，91022.