鋁冰發動機內流場的數值計算①

劉平安，王 良，王 璐，王 革

(哈爾濱工程大學 航天與建筑工程學院，哈爾濱 150001)

鋁冰發動機內流場的數值計算①

劉平安，王 良，王 璐，王 革

(哈爾濱工程大學 航天與建筑工程學院，哈爾濱 150001)

為了使用數值模擬的方法計算鋁冰發動機的性能，用顆粒表面反應模型和氣相反應模型模擬鋁顆粒在鋁冰發動機燃燒室中與水蒸氣的燃燒過程，用歐拉-拉格朗日方法計算顆粒沿軌跡的參數，分析了數值模擬的結果，并進行了相同尺寸的鋁冰發動機實驗，把數值模擬結果與實驗結果進行了比較。數值計算得到的燃燒室穩態工作壓強約為9.38 MPa，與實驗結果接近，燃燒室平均溫度為2950.65 K，相比熱力計算得到的推進劑燃燒溫度略低。通過對鋁冰發動機的內流場數值計算，得到了與實驗相符合的結果，驗證了數值計算模型的有效性。

鋁冰發動機；內流場；數值計算

0 引言

近年來低溫固體推進劑成為新的研究方向，而金屬燃料特別是納米金屬燃料在推進劑中的應用也成為重要研究方向[1]。納米金屬燃料點火溫度低、燃燒效率高，在推進劑中作為添加劑時能大大提高其燃速[2]，其超大的比表面積又使其能充當粘結劑以提高推進劑的力學性能，最終可作為推進劑基體和主要成分而不僅是添加劑。包覆后的納米金屬燃料用于低溫固化的固體推進劑中能長期保持活性，燃燒產物無污染，有望被用于深空探測，目前已被實驗驗證[3]。鋁冰發動機是低溫固體推進劑和包覆納米金屬燃料結合的產物，配方中只有納米鋁顆粒和水，鋁在推進劑中充當基體、燃料和粘合劑。實驗表明，納米鋁水混合物的燃燒效率很高[4]，其理論性能參數相比某些復合推進劑要高很多[5]，雖然國外鋁冰發動機(鋁水比為0.71∶1)實驗的性能效率很低[3]，但國內的鋁冰發動機(鋁水比為1∶1)實驗結果很好，這證明了鋁冰發動機的實際可行性。在鋁冰發動機實驗中發現按照傳統理論(熱力計算和零維內彈道理論)設計的燃燒室壓強遠高于實驗壓強，這使發動機實驗屢屢失敗，浪費人力物力成本。燃燒室壓強理論與實驗的誤差主要是由推進劑燃燒和兩相不平衡效應造成的，因此對鋁冰發動機性能的預估，如燃燒室壓強的理論估算需要采用考慮燃燒和兩相不平衡效應的新方法。

與常規復合推進劑發動機相比，鋁冰發動機主要不同在于推進劑中有大量的金屬鋁顆粒。鋁顆粒在發動機中燃燒和運動時要經歷相變、團聚、脫離燃面、點火、氣相反應、表面反應、氧化鋁凝結、顆粒成長、顆粒碰撞團聚和顆粒破碎等過程[6]。研究鋁顆粒燃燒過程的理論模型主要有半經驗模型、化學機理模型和CFD模型[7]。其中CFD模型可同時計算兩相流動和燃燒過程，得到整個內流場的流動情況，用于發動機工作過程模擬。

為了計算鋁冰發動機的理論性能，驗證計算模型并得到鋁冰發動機內流場參數，本文用CFD模型對鋁冰發動機的燃燒和流動過程進行計算，并把得到理論燃燒室壓強與實驗比較，以驗證數值模型在計算發動機性能(如燃燒室壓強)時的可行性。

1 數值計算方法

1.1 物理模型

鋁冰發動機裝藥為錐柱形，內部結構如圖1所示。發動機的外徑160 mm，內徑100 mm，噴管喉徑18 mm，裝藥錐角為19°，在工作初始時刻(燃層厚w=0 mm)，裝藥最大長度為300 mm。計算假定初始時刻燃面全部被點燃，工作過程中各處燃速相同，因此在不同燃層厚時燃面位置如圖1所示燃層厚w=0、10、20 mm。

1.2 計算(CFD)模型

本文用Fluent軟件中的歐拉拉格朗日模型和可實現k-ε模型來模擬鋁冰發動機內流場中的湍流兩相流，用組分輸運和渦耗散概念(EDC)模型模擬氣相混合物流動以及鋁和水蒸氣的氣相反應過程，用顆粒表面反應模型模擬顆粒的表面燃燒過程。其中歐拉-拉格朗日模型假設顆粒不占體積，沿運動軌跡計算其參數，該模型的控制方程如下[8]：

(1)

(2)

其中各參數含義可參見文獻[8]，該方程組通過有限體積法數值求解，計算格式為一階迎風格式。湍流的參數也用相同方法求解k和ε的控制方程得到。

鋁顆粒在燃燒時同時發生鋁蒸氣的氣相反應及顆粒表面反應：

Alg+1.5H2O→0.5 Al2O3g+3H2↑

Als+1.5H2O→0.5 Al2O3s+3H2↑

從文獻[8]可知，參加2種反應的鋁質量比為4∶1。化學反應造成氣相中鋁、水蒸氣等組分含量發生變化，而各組分含量通過氣相組分輸運方程求得：

(3)

在顆粒反應的過程中，首先發生氣相反應，然后再發生顆粒表面反應。發生氣相反應時，鋁顆粒吸熱蒸發為鋁氣，然后與水蒸氣反應。反應的速率由顆粒蒸發速率和化學動力學速率控制，其中顆粒蒸發速率的計算式為

-dmp/dt=A0fv，0mp，0

(4)

式中mp為顆粒質量；mp，0為顆粒初始質量；A0為蒸發速率常數，計算假定A0=5×105s-1；fv，0為氣相反應鋁占鋁顆粒初始質量的百分量，fv，0=0.8。

發生顆粒表面反應時，鋁直接以凝相的形式與水蒸氣反應，反應速率由氣相組分擴散速率和化學動力學速率控制。氣相反應鋁已經蒸發，因此在計算鋁顆粒表面反應的過程中，假定不存在被氧化鋁包裹的未完全反應的鋁，顆粒中鋁最終全部變為氧化鋁。這一過程中假設顆粒粒徑不變，顆粒中鋁/氧化鋁的消耗/產生速率通過式(5)計算[9]：

(5)

式中Ap為顆粒表面積；Yj為顆粒組分j的質量含量；ηr為影響因子，ηr=1；Rj，r為由顆粒表面反應化學動力學速率、顆粒粒徑、顆粒溫度、氣相溫度、氣相壓強等共同確定的單元反應速率。

上述2種反應的化學動力學速率均由阿累尼烏斯公式計算，但是目前有關鋁反應的動力學參數的研究很少且并不準確，本文計算暫先假定氣相反應的活化能Ea，r1=100 J/mol，指前因子A1=1.14×107，不受溫度影響；顆粒表面反應活化能Ea，r2=7.9 J/mol，指前因子A2=2×10-3。氣相化學反應生成熱全部被氣相吸收，顆粒表面反應的生成熱被氣相吸收75%，剩余的全部被顆粒吸收。

在納米鋁冰發動機中，由于參與反應的鋁顆粒粒徑很小，燃燒時以單顆粒燃燒為主[10]。現有一發動機裝藥所用納米鋁粒徑為以80 nm(0.08 μm)為中心的分布，下面計算假定該粒徑分布如圖2所示，粒徑分布函數為Rosin-Rammler。

計算假設所有顆粒全部為球形顆粒，各顆粒內部溫度均勻，顆粒表面的包覆層在顆粒進入燃燒室時已經被去除(包覆層通常是在高溫下易分解的物質)。由于本文模擬的是發動機點火后裝藥穩態燃燒和內流場流動的情況，所以計算選定的初始條件為：pc=8 MPa，Tc=3000 K。

1.3 發動機燃燒室壓強估算方法

在發動機實驗前，由壓強預估結果設計的裝藥和噴管決定了發動機能否正常工作。通常這一估算用熱力計算和零維內彈道理論完成。熱力計算實質上是用最小Gibbs自由能方法或其他方法計算推進劑的完全燃燒產物參數，并用“兩相平衡流”假設計算噴管中的兩相流[11]。熱力計算得到鋁冰推進劑的配方為Al∶H2O= 1∶1時，推進劑的特征速度C*=1373.4 m/s[12]。該特征速度值通常會直接用于零維內彈道(壓強)計算。

發動機零維內彈道理論用質量守恒原則求解零維穩態燃燒室壓強，由于計算簡單而被廣泛使用。零維平衡壓強公式最終形式如下[13]：

pc=(ρpC*a·Ab/At)1/(1-n)

(6)

式中ρp、C*、a、n分別為推進劑密度、特征速度、燃速系數和壓強指數；Ab和At分別為燃面和發動機噴管喉部面積。

除零維內彈道理論外，燃燒室壓強的估算還可采用其他方法。這些方法通過流場計算求解燃燒室壓強，包括考慮燃燒的CFD方法[14]、考慮兩相流的數值方法[15]等。其中，CFD方法通過合適的簡化數學模型計算推進劑反應物的燃燒過程和兩相燃燒產物的流動過程，兩相流數值計算則僅計算兩相燃燒產物的流動過程，這2種方法相比零維內彈道理論考慮了兩相不平衡效應(以及推進劑的燃燒)對燃燒室壓強的影響，理論上結果更準確。

2 計算結果分析

計算選用Fluent軟件中密度基穩態求解器，通過反復迭代得到發動機內流場穩態解(流量誤差波動小于5%)，最終得到燃燒室壓強、燃燒室溫度、氣相反應和顆粒表面反應區厚度、計算監測的顆粒數量及變化特征等參數如表1所示。

從表1可看出，發動機穩態工作時燃燒室壓強在8.04～10.78 MPa變化，燃燒室壓強平均值為9.38 MPa，燃燒室溫度平均值為2950.6 K。

和所有固體火箭發動機一樣，鋁冰發動機工作過程中推進劑組元的分解/相變產物以較低溫度進入燃燒室，在燃面附近迅速發生反應放熱，由此建立起燃燒室的高溫高壓環境。在所選的發動機的3個工作狀態下，計算得到發動機內流場中氣相總壓沿發動機軸線的變化情況如圖3所示。

表1 不同發動機工作狀態下數值計算結果

注：1)計算按燃層厚分為3個狀態，已在圖1中標出。

在鋁冰發動機中，貢獻壓強的工質主要是氫氣，隨著發動機工作的進行，燃燒室壓強呈下降趨勢，壓強值從10.78 MPa變化到8.04 MPa。燃氣通過噴管時膨脹做功，壓強迅速下降，在噴管出口下降為環境壓力。鋁和水在反應過程中會放出大量熱量，建立燃燒室的高溫環境，計算得到燃燒產物總溫在發動機不同工作時刻的變化情況如圖4所示。

由圖4可見，推進劑反應物從燃面進入燃燒室時溫度很低，但通過反應區后溫度迅速升高，放出反應熱并在燃燒室內形成了一個較均勻的高溫環境。該高溫環境維持了推進劑的持續燃燒，而隨著燃面推移，燃燒室溫度在2900～3000 K之間稍有變化，平均值為2950.65 K。用熱力計算軟件CEA[11]對相同配方的鋁冰推進劑進行計算，在燃燒室壓強為9 MPa的情況下得到其絕熱燃燒溫度為3054.76 K。可見數值模擬得到的發動機燃燒室溫度相比推進劑絕熱燃燒溫度偏低，分析發現這是由于數值模擬忽略了氧化鋁的凝結放熱而造成的。鋁的氣相反應產物是氧化鋁氣，目前并沒有研究證明這一物質實際存在，在簡化分析和計算中，可用它來代替鋁的氣相燃燒產物[16]。實際上氧化鋁在流動過程中會因溫度下降而逐漸凝結為氧化鋁顆粒，并最終隨燃氣流出。研究表明，這些新凝結的顆粒相比由鋁顆粒直接變成的氧化鋁顆粒尺寸更小，更容易與燃氣達到兩相平衡，所以這些顆粒往往被研究者們用歐拉方法處理成氣相的一部分[8，17]。因此，本文為了簡化計算而忽略了這一部分氧化鋁的凝結過程。雖然計算得到的燃燒室溫度與實際有偏差，但是由于氧化鋁密度大，對壓強影響小，計算表明即使氧化鋁氣全部變為氧化鋁顆粒，壓強的計算結果變化也并不明顯。另外，從圖4還可看出，在發動機頭部某些位置存在溫度局部變化較大的情況，這是由于一些顆粒在發動機頭部運動不能向外排出，在這些區域聚集并放出大量熱量導致的。在計算選定的三個狀態下，內流場中氣相馬赫數沿發動機軸線的變化如圖5所示。

從圖5可見，燃燒室中氣流速度很低，通過噴管時流速迅速增大，并在噴管出口達到超音速(Ma>2.2)，在發動機工作的不同時刻，出口氣相馬赫數變化不大。鋁顆粒燃燒時，氣相反應和顆粒表面反應同時發生，在圖1中所選定的3條監測線上(這3條監測線分別從垂直于燃面方向遠離燃面，長度均為10 mm)，可給出氣相反應速率、顆粒表面反應速率和化學反應熱沿監測線遠離燃面方向的變化分別如圖6、圖7所示。

從圖6可看出，2個反應在監測線長度范圍內就已經完成，在監測線上，反應速率大于0的區域(反應區)很薄，其中氣相反應區的平均厚度約為7.7 mm，顆粒表面反應的反應區平均厚度約為1.4 mm，具體的厚度值已在表1中列出。鋁顆粒進入燃燒室后很快發生反應，反應速率最大值在貼近燃面的位置，最大氣相反應速率為4×105～1.2×106kmol/(m3·s)，最大表面反應速率為50～110 kmol/(m3·s)。推進劑反應物經過反應區快速完成反應后，其產物把反應所產生的熱量帶入燃燒室，因此相比反應速率而言，反應的放熱區向遠離燃面的方向有一定程度的推移和放大。

從圖7可看出，即使在反應速率為0的區域，反應熱仍不為0，反應速率最大的區域反應熱也并非最大，這正是由于反應物邊流動邊反應，把熱量帶出反應區，并向周圍擴散造成的。假定反應放熱量大于0的區域為反應放熱區，在不同燃層厚下反應放熱區厚度的近似值如表1所示，其平均厚度約為12 mm，可見放熱區厚度比反應區厚度大，這是由于燃燒產物的流動造成的。通過放熱區后混合物的溫度迅速升高，在反應放熱區附近形成很大的溫度梯度，如圖4所示。

顆粒的燃燒和運動過程對內流場影響較大，圖8給出了不同燃層厚下數值計算所監測的顆粒中鋁質量占顆粒總質量(鋁+氧化鋁)的百分數沿顆粒軌跡的變化過程。

由圖8可見，由于顆粒粒徑較小反應速率大，鋁顆粒在燃面附近很快就反應完成變成氧化鋁。大部分顆粒隨燃氣從噴管流出，但仍有少部分顆粒最終會沉積在燃燒室中，形成殘渣。計算監測的顆粒變化的最終結果已在表1中列出，可以看出，不同燃層厚下計算監測的顆粒總數并不相同，而所有顆粒中包含氣相反應顆粒、表面反應顆粒和沉積的顆粒，顆粒的這些變化過程可用1.2節所給列出的模型近似計算。在燃燒室中沉積的顆粒會造成沉積區域附近流場有較大溫度變化(見圖4)。

3 計算結果的驗證與討論

3.1 內流場數值結果的驗證

為了驗證數值計算結果的有效性，下面給出與圖1相同尺寸的鋁冰發動機在相同的工作條件(水燃比1∶1)下實驗測量得到的燃燒室壓強曲線，如圖9所示。

從圖9可看出，發動機穩態工作時間約為1 s，穩態工作時平均壓強約為9.33 MPa，這一壓強值與數值模擬得到的平均燃燒室壓強(9.38 MPa，如表1)很接近。數值計算主要是對發動機的3個穩定工作狀態進行的，假定在這3個狀態下燃燒室實驗壓強值如壓強曲線中3個“╋”點所示。為了比較流場數值計算得到的燃燒室壓強、零維內彈道理論得到的燃燒室壓強、鋁冰完全燃燒產物兩相流平衡流流場計算得到的燃燒室壓強、完全燃燒產物兩相非平衡流計算得到的燃燒室壓強與實驗測量得到的壓強的差別，把各個燃層厚下壓強結果列如表2所示(其中兩相流數學模型參見文獻[18])。

從表2可看出，數值模擬計算的燃燒室壓強相比零維內彈道理論的計算結果更接近實驗值，且考慮燃燒的CFD計算結果相比實驗的平均誤差僅為4.62%，小于工程誤差(5%)。使用兩相非平衡流方法計算鋁冰推進劑兩相燃燒產物的流動時，燃燒室壓強的計算結果與實驗誤差為8.52%，相比零維內彈道的計算結果，計算精度已有很大的提高。在推進劑熱力計算和零維內彈道理論中都使用了兩相平衡流假設，從表2壓強計算結果可看出，在兩相平衡流假設下計算的燃燒室壓強與實驗值的誤差為18.34%，相比零維內彈道計算結果的誤差(50.74%)已經有很大的提高，但其計算精度仍比兩相非平衡流模型和CFD模型低，這說明在計算中應該考慮兩相非平衡、推進劑的燃燒過程等因素的影響。

表2 不同方法得到的鋁冰發動機燃燒室壓強

上述結果證明，目前常用的零維內彈道理論在鋁冰發動機中使用時燃燒室壓強的計算結果相比實驗誤差較大，而用CFD模型計算的鋁冰發動機燃燒室壓強結果與實驗相比誤差很低。這說明用CFD方法得到的鋁冰發動機內流場參數是可信的，這一方法可推廣到其他高金屬含量發動機。

3.2 鋁冰發動機燃燒室壓強預測公式

為了得到簡單的公式以在鋁冰發動機實驗前對發動機燃燒室壓強進行準確估算，需要對傳統壓強計算方法進行修正。根據本文數值計算(CFD模型)得到的發動機不同工作狀態下的燃燒室壓強結果，與零維內彈道理論得到的燃燒室壓強值進行比較，并對后者進行修正，得到鋁冰發動機零維燃燒室壓強pc，理的修正式如下：

pc，修正=0.7055pc，理

(7)

實驗前可用該公式修正鋁冰發動機的零維燃燒室壓強計算結果，這樣就在壓強計算中近似考慮了推進劑燃燒和兩相燃燒產物的流動過程對壓強的影響。對于同類鋁冰發動機，該公式有望提高壓強計算的準確度并降低實驗失敗的機率。

4 結論

(1)從鋁冰發動機的實驗曲線和數值模擬結果可以看出，納米鋁的反應活性很強，鋁和水蒸氣在推進劑燃面附近就能迅速反應且放出大量熱量，鋁冰發動機能夠正常工作。

(2)在含金屬發動機中推進劑燃燒產物是兩相不平衡流，使用傳統零維內彈道理論計算燃燒室壓強時，誤差較大，而使用兩相流數值計算的方法，求解兩相流控制方程，能考慮兩相不平衡效應。使用數值計算方法(CFD模型)可有效預估鋁冰發動機的燃燒室壓強，該方法可用于工程計算。

(3)用壓強修正公式能夠簡單、快捷地在鋁冰發動機燃燒室壓強計算中考慮推進劑燃燒和兩相不平衡效應對壓強的影響，提高零維燃燒室壓強計算結果的準確度，從而提高實驗成功率，降低實驗成本。

[1] 龐愛民,黎小平.固體推進劑技術的創新與發展規律[J].含能材料,2015,23(1):3-6.

[2] 胡潤芝,張小平.高能推進劑燃燒效率研究和實測比沖預估[J].推進技術,2001,22(5):415-417+436.

[3] Risha G A,Connell Jr T L.Novel energetic materials for space propulsion[R].Nasa:DTIC Document:A546818,2011.

[4] Risha G A,Sabourin J L.Combustion and conversion efficiency of nanoaluminum-water mixtures[J].Combustion Science and Technology,2008,180(12):2127-2142.

[5] Ingenito A,Bruno C.Using aluminum for space propulsion[J].Journal of Propulsion and Power,2004,20(6):1056-1063.

[6] Geisler R L.A global view of the use of aluminum fuel in solid rocket motors[C]//38th AIAA/ASME/SAE/ASEE Joint Propulsion Conference & Exhibit,2002.

[7] 姜治深,梁晶晶.固體推進劑中鋁燃燒理論研究展望[J].飛航導彈,2013(10):88-91,96.

[8] Sabnis J S.Numerical simulation of distributed combustion in solid rocket motors with metalized propellant[J].Journal of Propulsion and Power,2003,19(1):48-55.

[9] Ansys_Inc.Fluent 15.0 Help[R].2013.

[10] 李穎,宋武林.納米鋁粉在固體推進劑中的應用進展[J].兵工學報,2005,26(1):121-125.

[11] Gordon S,Mcbride B J.Computer program for calculation of complex chemical equilibrium compositions and applications[R].Nasa,Cleveland,OH,United States:NASA Lewis Research Center,1994:58.

[12] 韓秀杰.鋁-冰固體推進劑燃燒性能研究[D].哈爾濱:哈爾濱工程大學,2013.

[13] 李宜敏,趙元修.固體火箭發動機原理[M]. 北京:國防工業出版社,1985.

[14] Puskulcu G.Analysis of 3-D grain burnback of solid propellant rocket motors and verification with rocket motor tests[D].MS.Thesis,Dept.of Mechanical Engineering,METU,2004.

[15] 陳軍.固體火箭發動機零維兩相內彈道研究[J].彈道學報,2013,25(2):39-43.

[16] Daniel E.Eulerian approach for unsteady two-phase solid rocket flows with aluminum particles[J].Journal of Propulsion and Power,2000,16(2):309-317.

[17] Najjar F M,Massa L.Effects of aluminum propellant loading and size distribution in BATES motors:A multiphysics computational analysis[C]//41st AIAA/ASME /SAE/ASEE Joint Propulsion Conference & Exhibit,2005.

[18] 方丁酉.兩相流動力學[M].長沙:國防科技大學出版社,1988:187,545.

(編輯：呂耀輝)

Numerical simulation of aluminum-ice solid propellant motor

LIU Ping-an，WANG Liang，WANG Lu，WANG Ge

(College of Aerospace and Civil Engineering，Harbin Engineering University，Harbin 150001，China)

To calculate the theoretical performance of the Aluminum-Ice Solid Rocket Motor(ALICE SRM),the particle surface reaction model and the gas phase reaction model were used to simulate the burning process of aluminum particle and the steam in the combustion chamber of the ALICE SRM.The Euler-Lagrange model was used to calculate the particle parameters along their trajectories.The numerical results were analyzed and compared with the experimental results of the same sized ALICE SRM.The numerical chamber pressure of the ALICE SRM is about 9.38 MPa,which is very close to the experimental value.The calculated average chamber temperature is 2950.65 K,which is lower than that of the thermo-chemistry calculation result.By the numerical simulation of the internal flow-field of the ALICE SRM,a good agreement between the numerical results and the experimental data was obtained,which clearly verified the numerical models in predicting ALICE SRM performance.

aluminum-ice solid rocket motor；internal flow-field；numerical calculation

2016-07-12；

2016-12-15。

中央高校基本科研基金(HEUCFD1404；HEUCFD1502)。

劉平安(1982—)，男，副教授，研究方向為金屬燃料發動機。E-mail：liupingan631@126.com

V435

A

1006-2793(2017)04-0425-07

10.7673/j.issn.1006-2793.2017.04.005