冀 鵬，雷凡培，梁樹強，肖明杰

（1. 西安航天動力研究所液體火箭發(fā)動機重點實驗室，西安，710100；2. 中國船舶集團有限公司，北京，100044； 3. 航天推進技術研究院，西安，710100）

0 引 言

肼類單組元發(fā)動機是運載火箭和各類航天器姿態(tài)控制、軌道控制使用最廣泛的推進系統(tǒng)，其性能以及穩(wěn)定性直接影響著運載火箭和航天器的控制精度和可靠性[1,2]。近年來隨著航天器飛行任務的增多，對單組元發(fā)動機的環(huán)境適應性愈加嚴苛，單組元肼發(fā)動機在冷起動時，時常出現(xiàn)起動響應延遲長、起動壓力峰大的現(xiàn)象。一方面較長的起動響應延遲難以滿足總體對推進系統(tǒng)的響應要求，另一方面較高的起動壓力峰也降低了發(fā)動機的結構可靠性，可能導致發(fā)動機結構破壞。亟待對肼類單組元催化分解發(fā)動機冷起動工作過程開展深入的研究工作。

采用催化分解產生熱燃氣的肼類單組元發(fā)動機，推進劑進入推力室后涉及在催化床這種多孔介質中的多相多尺度流動、均相熱分解和異相催化分解等復雜的物理化學過程，目前對該過程的研究主要依靠試驗。Hearn[3]在研究脈沖工作對催化床壽命影響時發(fā)現(xiàn)，在脈沖工作模式下冷起動壓力峰不會發(fā)生在較短的脈沖過程中。Kagawa等[4,5]采用可視化技術研究催化床內液態(tài)肼流動狀態(tài)，發(fā)現(xiàn)相較于溫啟動和熱啟動，冷起動時液態(tài)肼的滲透面積更大、滲透距離更遠。此外，火星著陸飛船[6]、阿里安5號運載火箭[7]、NASA空間飛行器[8,9]等所采用單組元發(fā)動機研制過程中，也主要采用試驗手段來考核和鑒定發(fā)動機的冷起動能力并評判其對發(fā)動機性能的影響，發(fā)現(xiàn)過低的起動邊界溫度會導致肼推力器的使用壽命大大降低。單組元發(fā)動機的冷起動仿真研究方面目前尚未有公開文獻發(fā)表。

一般認為溫度較低時化學反應速率會減慢，但在單組元肼類發(fā)動機冷起動過程中，會出現(xiàn)數(shù)倍于溫啟動或熱啟動的響應延遲，并且延遲過程中室壓維持不變，然后突然產生明顯的壓力峰，該現(xiàn)象與低溫導致化學反應變慢機理并不契合。因此考慮該現(xiàn)象與微觀催化分解過程有關。肼分解用催化劑并不是簡單的球形或柱形顆粒結構，而是由復雜的浸漬、焙燒等制備工藝得到的具有眾多微孔通道的層狀骨架結構[10]。Sanglovanni等[11]通過試驗觀察到在低溫情況下，液態(tài)肼滴定到催化劑顆粒后迅速浸濕了疏松的催化劑顆粒，在試驗后發(fā)現(xiàn)催化劑顆粒發(fā)生了明顯的破碎。趙許群[12]在對低溫下催化劑的失活機理進行分析時提到，在低溫下肼會在催化劑微孔內短時劇烈分解產生較大內應力導致催化劑顆粒破碎。

綜上，現(xiàn)在仍然缺乏對肼類單組元發(fā)動機冷起動過程的完整過程認識和可靠的數(shù)值仿真方法，導致在發(fā)動機研制過程中主要依賴試驗來預判和評價發(fā)動機的低溫起動特性。因此結合催化劑微孔通道內流動-反應和經典燃燒時滯理論，開展肼類單組元發(fā)動機冷起動過程數(shù)值仿真分析，研究初始溫度對冷起動壓力與延遲時間的影響，揭示發(fā)動機冷起動延遲長、壓力峰高的主導因素。

1 發(fā)動機典型組件的動力學模型

對于恒壓擠壓式的單組元姿控發(fā)動機，工作過程中貯箱壓力基本保持不變，因此假設貯箱內氣體壓力恒定，主要建立起液路推進劑供應系統(tǒng)以及催化床的動力學模型。

1.1 管路模型

在建立液體管路模型時不考慮管路的變形，認為管道為剛性的；管路內液體是一維的，管道橫截面上的速度分布是均勻的。本文主要關注的發(fā)動機起動過程屬于低頻動力學問題，因此采用分段集中參數(shù)的管路模型即可滿足精度要求。根據(jù)經驗，一種比較保守的方法是將分段管長取波長的4%以內，即管路分段長度應滿足：Lmax≤0.04a/fmax。此外，毛細管、閥門流道等在一定程度上也可以等效為液體管路。

對于特定的一段長度、直徑、截面積和體積分別為L、d、A、V的管路，考慮慣性、粘性時管路瞬變流：

式中qm，p1，p2分別為管路流量，管路入口和出口處壓力；j為管內液體的慣性系數(shù)，ALj/=；ξ為總流阻系數(shù)，由管路沿程損失和局部損失組成，取決于流體管路的形狀和尺寸、管壁粗糙度及流體在管路的流動狀態(tài)等，其中流體流動狀態(tài)常采用雷諾數(shù)來表征。

考慮液體壓縮性的方程：

式中z為管路的流容系數(shù)，z=V/a2，表征管路流體的壓縮性。

1.2 催化床模型

本文采用反映平均效應的零維模型。在建立催化床模型時作如下假設：分解產物為理想氣體；不考慮氣體的不均勻性和波動過程；不考慮催化床內催化劑活性點處失活的影響。

1.2.1 冷起動延遲模型

在單組元發(fā)動機冷起動時，考慮到此時催化床和推進劑的溫度都比較低，液體推進劑噴入催化床內后難以快速蒸發(fā)而以液態(tài)形式流入催化劑顆粒的間隙中，同時催化劑上微孔直徑較小（直徑范圍從10A到104A），在毛細作用下滲透到催化劑顆粒微孔的內部。進入微孔后推進劑緩慢蒸發(fā)，在催化劑微孔內部形成氣液分界面，部分氣相推進劑通過質量擴散到微孔壁面的活性位置分解得到氨氣、氫氣和氮氣等。這些分解氣體被催化劑外的液體推進劑堵塞在微孔內，隨著分解產物逐步增多微孔內部壓力逐漸升高，推進劑分界面逐步向微孔外部移動，微孔內部壓力最終會大于毛細壓力而使得分界面被擠壓到顆粒外部。大量的分解產物從微孔中溢出，氣相推進劑得以實現(xiàn)催化劑間隙到微孔活性位置的快速擴散分解，在這兩者共同的作用下催化床內產生了較大的壓力峰。將冷起動延遲時間定義為從液體推進劑噴入催化床內到氣體產物溢出催化劑微孔這段時間。

基于上述冷起動過程催化床內部復雜過程的分析，下面將對催化劑顆粒內部進行建模分析，計算氣液分界面在顆粒微孔內部壓力的作用下逐漸推移到微孔外過程，得到單組元發(fā)動機冷起動延遲時間。

首先以催化劑顆粒內部氣體為研究對象，分析催化劑顆粒微孔內部壓升過程。根據(jù)微孔毛細理論[13]，微孔內部氣相產物沿顆粒徑向的平均速度通過泊肅葉方程描述如下：

式中vm為微孔內部氣相的平均速度；R為大孔的特征半徑；μg為氣相產物的粘性系數(shù)；δ為催化劑微孔彎曲因子；P為微孔內壓力；r為流體徑向半徑。

式中J為氣相產物徑向單位面積上的質量流量；γν為空隙分數(shù)，指單位體積顆粒內氣相填充的比例；M為氣體產物的平均摩爾質量；Ru為氣體常數(shù)；T為氣體溫度，假設在較短的冷起動延遲時間內氣相產物的溫度與初始溫度相差不大。

在微孔內考慮氣體的連續(xù)方程：

聯(lián)立式（4）和式（5）可得催化劑顆粒內部壓力：

在邊界條件方面，初始時刻微孔內壓力為環(huán)境壓力；而在顆粒中心的壓力梯度為零，還有一個邊界條件是在微孔內部移動的氣液分界面氣相質量守恒，即微孔內部氣相產物的質量變化等于在分界面上蒸發(fā)和顆粒內部擴散氣態(tài)肼的質量變化。則初邊值條件如下：

式中DA為氣相推進劑在其分解產物中的擴散系數(shù)，。

為了求解內部壓力方程，還需要知道分界面隨時間的變化ai(t)，以及在分界面上氣相肼的摩爾質量沿微孔徑向分布。顆粒內部氣液分界面的徑向位置可用描述液體在長管中流動規(guī)律的泊肅葉方程得到：

式中μl為液態(tài)推進劑粘性系數(shù)；Pm為毛細壓力，；σ為液態(tài)推進劑表面張力，隨溫度變化；θ為接觸角。

根據(jù)在分界面后局部氣態(tài)肼考慮分解反應的質量守恒方程得到：

式中y為液體推進劑滲入催化劑顆微孔內的長度；為液體推進劑的密度；ks(T)為推進劑的異相反應速率常數(shù)，ks=1010exp(-1389/T)。

這樣得到了完整的控制方程與初邊值條件，經過無量綱化和小擾動法計算得到分界面位置ai(tlign)=0時，即氣體分解產物從微孔中溢出，此時得到了冷起動延遲時間tlign為[14]

式中 (ρA)vp為推進劑的飽和蒸汽密度，與溫度相關[1]。

1.2.2 催化床建壓過程模型

考慮到催化床內復雜的物理化學過程，引入燃燒室時滯模型[15]，假設在該轉化時間內推進劑的轉化速率均勻，且分解后的氣體在任何瞬時分布均勻，即氨的解離過程也近似認為與DT-3的分解過程同時進行。 則液體推進劑及氣體分解產物在催化床中的質量變化分別為

式中mout是催化床出口的氣體質量流量，；At為噴管的喉部截面積；ml、mg分別為催化床內積存的液體推進劑質量和氣體質量；ql、分別為流入催化床內的液體推進劑質量流量和流出推力室的氣體質量流量；τ為推進劑的轉化時間，可根據(jù)推進劑的分解延遲期試驗[16]得到。

由于催化床內裝載著緊密的催化劑，因此推進劑通過催化床時需要克服阻力稱為床流阻，床流阻常采用經驗公式計算[1]。

1.3 其他組件模型

用于單組元姿控發(fā)動機冷起動過程仿真的其他組件，如限流圈流阻模型，電磁閥考慮了閥芯動作過程的Γ-m模型，電磁閥后集液腔考慮氣體多變過程的充填排氣模型等[17]。限于文章篇幅不再詳述。

2 仿真模型的建立及校驗

某次單臺發(fā)動機地面試車考核了發(fā)動機在低溫下的工作特性，試驗系統(tǒng)如圖1所示。在該單機地面試驗系統(tǒng)中所采用的推進劑為液體單推-3[18]（DT-3）。由于試驗系統(tǒng)中流量計的采樣頻率較低，不能反映發(fā)動機起動過程中的快速流量變化，故在此以推力室的室壓表征單組元發(fā)動機的冷起動特性。

圖1 單臺發(fā)動機地面試驗系統(tǒng)示意Fig.1 Engine Ground Test System Schematic

基于第1節(jié)所建立的典型組件數(shù)學模型，在MWorks軟件平臺上搭建該型單臺發(fā)動機的仿真模型，數(shù)值計算方法上采用4階Range-Kutta方法。計算結果和試驗測量結果如圖2所示。

圖2 單臺發(fā)動機冷起動過程仿真與試驗結果對比Fig.2 Simulation Results Compare with Ground Test Results of Engine Cold Starting Process

圖2中試驗值與計算值的比較結果表明，發(fā)動機室壓的計算值與試驗值變化趨勢吻合較好，但室壓計算值較試驗值偏高，引起這一誤差的主要原因是計算值是包含了催化床流阻部分的整個催化床平均壓力，而試驗測量值僅為推力室噴管入口處的壓力；并且發(fā)動機在冷起動時催化床分解效率是隨發(fā)動機工作時間逐漸升高的，而仿真中將分解效率設為恒定的穩(wěn)態(tài)分解效率，這2個因素導致了室壓計算值較高。其次，計算得到的室壓建壓曲線相比試驗建壓曲線更為陡峭，壓力峰的計算值稍快于試驗值，則是因為催化床模型采用了分解時滯模型反映推進劑的分解過程，認為積存推進劑全部在瞬間發(fā)生快速分解，而發(fā)動機實際過程中積存的推進劑也是根據(jù)進入催化床順序先后逐步發(fā)生快速分解的，同時測量系統(tǒng)可能無法完全捕捉到室壓快速上升過程，考慮到這兩方面因素從而導致計算建壓過程快于試驗測量的建壓過程。

國外TRW公司對自發(fā)型催化劑Shell-405及其單組元肼發(fā)動機進行了充分的真空、地面點火試驗研究[19]，得到了不同初始溫度下發(fā)動機的冷起動響應時間，但文獻中未能給出試驗所用發(fā)動機系統(tǒng)的具體參數(shù)，因此選取起動響應時間作為關注量，起動響應時間指推進劑進入催化床到推力室內壓力達到穩(wěn)態(tài)室壓的5%所需時間。圖3給出了TRW公司地面試驗用發(fā)動機和本文所研究單組元發(fā)動機的冷起動響應時間隨初始溫度變化的試驗值和仿真計算曲線。這里的溫度是指推進劑與催化床的初始溫度，考慮發(fā)動機地面試驗時第1次起動的催化劑與推進劑溫度相差不大，在此假設這兩者溫度一致。

圖3 響應時間隨初始溫度變化對比Fig.3 Comparison of Test Data and Calculated Values of Response Time with Initial Temperature

本文的計算結果與文獻19中試驗數(shù)據(jù)的吻合情況較好。在初始溫度降至-20 ℃以下時，計算值與文獻19試驗值的誤差有增大的趨勢，分析認為是該溫度范圍接近推進劑的冰點-30 ℃，數(shù)學模型中所采用的推進劑物性參數(shù)計算公式在該范圍內適用性變差所導致。

3 仿真結果分析

3.1 單組元發(fā)動機冷起動過程研究

以第2節(jié)驗證過的仿真模型為基礎，以發(fā)動機設計溫度下的狀態(tài)參數(shù)值為參考，討論和研究初始溫度對發(fā)動機冷起動過程中室壓和響應時間的影響。圖4給出了不同初始溫度下發(fā)動機冷起動過程中無量綱室壓cP隨時間的變化。由圖4可知，催化床與推進劑的初始溫度對發(fā)動機的冷起動過程影響較大。初始溫度越低，起動響應時間越長，冷起動壓力峰越高，發(fā)動機的冷起動性能越差，降低了發(fā)動機的結構可靠性。還可以看到，在出現(xiàn)冷起動壓力峰后壓力曲線還會伴隨“凹坑”。這種“凹坑”現(xiàn)象的產生是由于高壓力峰的存在，使得整個系統(tǒng)壓差變化，發(fā)動機供應系統(tǒng)出現(xiàn)短暫的推進劑供應不足導致推力室產生了較大的壓力波動。冷起動壓力峰越高，所產生的壓力“凹坑”越突出。在初始溫度為10 ℃時冷起動壓力峰已經高達設計值的4.5倍左右，并且伴隨了更為明顯的“凹坑”，壓力的最低值已經接近0值而后緩慢爬升至設計穩(wěn)定值，導致發(fā)動機達到穩(wěn)態(tài)工作壓力的時間明顯延長，嚴重影響了發(fā)動機的工作可靠性。

圖4 初始溫度對冷起動過程的影響Fig.4 Influence of Initial Temperature on Cold Starting Process

進一步模擬不同初始溫度下發(fā)動機的冷起動過程，圖5為初始溫度在0~30 ℃推力裝置的無量綱冷起動壓力峰值及響應時間隨初始溫度的變化?？梢钥吹?，冷起動壓力峰值和冷起動響應時間呈現(xiàn)正相關的關系。這意味著冷起動響應越快，冷起動壓力峰越低。

圖5 壓力峰值及峰值時刻隨初始溫度的變化Fig.5 Influence of Initial Temperature on Pressure Spike and Response Time

隨著初始溫度的升高，冷起動壓力峰和響應時間隨溫度的變化近似呈指數(shù)變化。對于所研究地面試驗系統(tǒng)中的單組元發(fā)動機，當催化床和推進劑的初始溫度從0 ℃升高到30 ℃時，壓力峰值降低了85.8%，響應時間減小了67.5%。催化床和推進劑的初始溫度顯著影響著單組元發(fā)動機的冷起動過程，尤其是在初始溫度較低時，即便是較小程度的提高初始溫度也可以明顯地改善發(fā)動機冷起動延遲長、壓力峰高的現(xiàn)象。

3.2 初始溫度對冷起動延遲的作用機制分析

由冷起動延遲過程分析可知，初始溫度主要影響著包括催化分解反應、液相推進劑蒸發(fā)以及氣相推進劑擴散這3種過程以及液相推進劑的表面張力，為進一步比較分析4種因素對冷起動延遲的影響程度，以20 ℃時各參數(shù)值為參照，對式（10）進行無量綱化后再取以為底的對數(shù)得到：

由于log函數(shù)不改變原函數(shù)的單調性，因此式（14）中右側分別代表著表面張力、初始溫度、飽和密度、氣相推進劑擴散系數(shù)以及分解反應速率對冷起動延遲時間的影響因子。由于推進劑的冰點為-30 ℃，因此計算初始溫度在-20~30 ℃情況下各影響因子的變化如圖6所示。

圖6 不同初始溫度下各項的影響因子Fig.6 Influencing Factors at Different Initial Temperature

由圖6可知，隨著初始溫度的升高，飽和密度項和分解反應項呈現(xiàn)一定的下降趨勢，而表面張力、擴散系數(shù)和溫度3項相應的影響程度逐漸提高。其中表面張力的作用隨溫度上升提高明顯，在大約10 ℃其作用已經超過了擴散系數(shù)和溫度項的影響，但是這3項的影響因子都小于5%。

在所研究的溫度范圍內，飽和密度項的影響因子一直處在較高值，并且超過其他4項之和的2倍以上。因此液體推進劑的蒸發(fā)過程是單組元肼發(fā)動機冷起動過程中導致延遲較長、壓力峰較高的主導因素；其次分解反應項的影響因子超過了10%，表明分解反應過程在冷起動過程中的影響也不容忽略。

3.3 不同推進劑的冷起動溫度敏感性

肼-70、無水肼和DT-3是3種常用肼類單組元推進劑[16,20]，本節(jié)討論相同條件下采用不同推進劑的單組元發(fā)動機冷起動特性。圖7為采用3種推進劑的發(fā)動機在-15~30 ℃初始溫度條件下的冷起動延遲時間。由圖7可知，在所計算溫度范圍內，無水肼發(fā)動機冷起動延遲最長，且隨初始溫度的降低，冷起動延遲的差距愈加明顯，-15 ℃時無水肼和DT-3的差值已達到1倍以上?？紤]3.1節(jié)得到的冷起動壓力峰與冷起動延遲時間正相關，故認為采用DT-3的單組元發(fā)動機在相同條件下較無水肼和肼-70具有更好的冷起動特性。

圖7 冷起動延遲時間隨初始溫度變化Fig.7 Influence of Initial Temperature on Response Time of Cold Starting

4 結 論

本文建立了反映冷起動延遲的單組元姿控發(fā)動機典型組件數(shù)學模型，分析了單組元姿控軌發(fā)動機在不同初始溫度下的冷起動特性，得到以下結論：a）單組元姿控發(fā)動機冷起動過程仿真模型的計算值與試驗值吻合較好，冷起動延遲時間與文獻中試驗值趨勢相同，驗證了所建立催化床冷起動延遲數(shù)學模型和相應仿真模型，可為工程設計提供理論支持。b）肼類單組元發(fā)動機冷起動過程中壓力峰與起動延遲正相關，并都與初始溫度近似指數(shù)變化的關系。c）相比于低溫下氣態(tài)推進劑的擴散過程和催化分解過程，液體推進劑進入催化床后的蒸發(fā)氣化過程是單組元發(fā)動機在較低溫度起動時延遲長、壓力峰高的主導因素。

相較于肼-70和無水肼，采用單推-3的單組元催化分解發(fā)動機具有更好的冷起動特性。