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

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

鋁冰發(fā)動機(jī)內(nèi)流場的數(shù)值計(jì)算①

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

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

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

鋁冰發(fā)動機(jī)；內(nèi)流場；數(shù)值計(jì)算

0 引言

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

與常規(guī)復(fù)合推進(jìn)劑發(fā)動機(jī)相比，鋁冰發(fā)動機(jī)主要不同在于推進(jìn)劑中有大量的金屬鋁顆粒。鋁顆粒在發(fā)動機(jī)中燃燒和運(yùn)動時(shí)要經(jīng)歷相變、團(tuán)聚、脫離燃面、點(diǎn)火、氣相反應(yīng)、表面反應(yīng)、氧化鋁凝結(jié)、顆粒成長、顆粒碰撞團(tuán)聚和顆粒破碎等過程[6]。研究鋁顆粒燃燒過程的理論模型主要有半經(jīng)驗(yàn)?zāi)Ｐ?、化學(xué)機(jī)理模型和CFD模型[7]。其中CFD模型可同時(shí)計(jì)算兩相流動和燃燒過程，得到整個(gè)內(nèi)流場的流動情況，用于發(fā)動機(jī)工作過程模擬。

為了計(jì)算鋁冰發(fā)動機(jī)的理論性能，驗(yàn)證計(jì)算模型并得到鋁冰發(fā)動機(jī)內(nèi)流場參數(shù)，本文用CFD模型對鋁冰發(fā)動機(jī)的燃燒和流動過程進(jìn)行計(jì)算，并把得到理論燃燒室壓強(qiáng)與實(shí)驗(yàn)比較，以驗(yàn)證數(shù)值模型在計(jì)算發(fā)動機(jī)性能(如燃燒室壓強(qiáng))時(shí)的可行性。

1 數(shù)值計(jì)算方法

1.1 物理模型

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

1.2 計(jì)算(CFD)模型

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

(1)

(2)

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

鋁顆粒在燃燒時(shí)同時(shí)發(fā)生鋁蒸氣的氣相反應(yīng)及顆粒表面反應(yīng)：

Alg+1.5H2O→0.5 Al2O3g+3H2↑

Als+1.5H2O→0.5 Al2O3s+3H2↑

從文獻(xiàn)[8]可知，參加2種反應(yīng)的鋁質(zhì)量比為4∶1?；瘜W(xué)反應(yīng)造成氣相中鋁、水蒸氣等組分含量發(fā)生變化，而各組分含量通過氣相組分輸運(yùn)方程求得：

(3)

在顆粒反應(yīng)的過程中，首先發(fā)生氣相反應(yīng)，然后再發(fā)生顆粒表面反應(yīng)。發(fā)生氣相反應(yīng)時(shí)，鋁顆粒吸熱蒸發(fā)為鋁氣，然后與水蒸氣反應(yīng)。反應(yīng)的速率由顆粒蒸發(fā)速率和化學(xué)動力學(xué)速率控制，其中顆粒蒸發(fā)速率的計(jì)算式為

-dmp/dt=A0fv，0mp，0

(4)

式中mp為顆粒質(zhì)量；mp，0為顆粒初始質(zhì)量；A0為蒸發(fā)速率常數(shù)，計(jì)算假定A0=5×105s-1；fv，0為氣相反應(yīng)鋁占鋁顆粒初始質(zhì)量的百分量，fv，0=0.8。

發(fā)生顆粒表面反應(yīng)時(shí)，鋁直接以凝相的形式與水蒸氣反應(yīng)，反應(yīng)速率由氣相組分?jǐn)U散速率和化學(xué)動力學(xué)速率控制。氣相反應(yīng)鋁已經(jīng)蒸發(fā)，因此在計(jì)算鋁顆粒表面反應(yīng)的過程中，假定不存在被氧化鋁包裹的未完全反應(yīng)的鋁，顆粒中鋁最終全部變?yōu)檠趸X。這一過程中假設(shè)顆粒粒徑不變，顆粒中鋁/氧化鋁的消耗/產(chǎn)生速率通過式(5)計(jì)算[9]：

(5)

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

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

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

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

1.3 發(fā)動機(jī)燃燒室壓強(qiáng)估算方法

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

發(fā)動機(jī)零維內(nèi)彈道理論用質(zhì)量守恒原則求解零維穩(wěn)態(tài)燃燒室壓強(qiáng)，由于計(jì)算簡單而被廣泛使用。零維平衡壓強(qiáng)公式最終形式如下[13]：

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

(6)

式中ρp、C*、a、n分別為推進(jìn)劑密度、特征速度、燃速系數(shù)和壓強(qiáng)指數(shù)；Ab和At分別為燃面和發(fā)動機(jī)噴管喉部面積。

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

2 計(jì)算結(jié)果分析

計(jì)算選用Fluent軟件中密度基穩(wěn)態(tài)求解器，通過反復(fù)迭代得到發(fā)動機(jī)內(nèi)流場穩(wěn)態(tài)解(流量誤差波動小于5%)，最終得到燃燒室壓強(qiáng)、燃燒室溫度、氣相反應(yīng)和顆粒表面反應(yīng)區(qū)厚度、計(jì)算監(jiān)測的顆粒數(shù)量及變化特征等參數(shù)如表1所示。

從表1可看出，發(fā)動機(jī)穩(wěn)態(tài)工作時(shí)燃燒室壓強(qiáng)在8.04～10.78 MPa變化，燃燒室壓強(qiáng)平均值為9.38 MPa，燃燒室溫度平均值為2950.6 K。

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

表1 不同發(fā)動機(jī)工作狀態(tài)下數(shù)值計(jì)算結(jié)果

注：1)計(jì)算按燃層厚分為3個(gè)狀態(tài)，已在圖1中標(biāo)出。

在鋁冰發(fā)動機(jī)中，貢獻(xiàn)壓強(qiáng)的工質(zhì)主要是氫氣，隨著發(fā)動機(jī)工作的進(jìn)行，燃燒室壓強(qiáng)呈下降趨勢，壓強(qiáng)值從10.78 MPa變化到8.04 MPa。燃?xì)馔ㄟ^噴管時(shí)膨脹做功，壓強(qiáng)迅速下降，在噴管出口下降為環(huán)境壓力。鋁和水在反應(yīng)過程中會放出大量熱量，建立燃燒室的高溫環(huán)境，計(jì)算得到燃燒產(chǎn)物總溫在發(fā)動機(jī)不同工作時(shí)刻的變化情況如圖4所示。

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

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

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

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

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

由圖8可見，由于顆粒粒徑較小反應(yīng)速率大，鋁顆粒在燃面附近很快就反應(yīng)完成變成氧化鋁。大部分顆粒隨燃?xì)鈴膰姽芰鞒?，但仍有少部分顆粒最終會沉積在燃燒室中，形成殘?jiān)?。?jì)算監(jiān)測的顆粒變化的最終結(jié)果已在表1中列出，可以看出，不同燃層厚下計(jì)算監(jiān)測的顆?？倲?shù)并不相同，而所有顆粒中包含氣相反應(yīng)顆粒、表面反應(yīng)顆粒和沉積的顆粒，顆粒的這些變化過程可用1.2節(jié)所給列出的模型近似計(jì)算。在燃燒室中沉積的顆粒會造成沉積區(qū)域附近流場有較大溫度變化(見圖4)。

3 計(jì)算結(jié)果的驗(yàn)證與討論

3.1 內(nèi)流場數(shù)值結(jié)果的驗(yàn)證

為了驗(yàn)證數(shù)值計(jì)算結(jié)果的有效性，下面給出與圖1相同尺寸的鋁冰發(fā)動機(jī)在相同的工作條件(水燃比1∶1)下實(shí)驗(yàn)測量得到的燃燒室壓強(qiáng)曲線，如圖9所示。

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

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

表2 不同方法得到的鋁冰發(fā)動機(jī)燃燒室壓強(qiáng)

上述結(jié)果證明，目前常用的零維內(nèi)彈道理論在鋁冰發(fā)動機(jī)中使用時(shí)燃燒室壓強(qiáng)的計(jì)算結(jié)果相比實(shí)驗(yàn)誤差較大，而用CFD模型計(jì)算的鋁冰發(fā)動機(jī)燃燒室壓強(qiáng)結(jié)果與實(shí)驗(yàn)相比誤差很低。這說明用CFD方法得到的鋁冰發(fā)動機(jī)內(nèi)流場參數(shù)是可信的，這一方法可推廣到其他高金屬含量發(fā)動機(jī)。

3.2 鋁冰發(fā)動機(jī)燃燒室壓強(qiáng)預(yù)測公式

為了得到簡單的公式以在鋁冰發(fā)動機(jī)實(shí)驗(yàn)前對發(fā)動機(jī)燃燒室壓強(qiáng)進(jìn)行準(zhǔn)確估算，需要對傳統(tǒng)壓強(qiáng)計(jì)算方法進(jìn)行修正。根據(jù)本文數(shù)值計(jì)算(CFD模型)得到的發(fā)動機(jī)不同工作狀態(tài)下的燃燒室壓強(qiáng)結(jié)果，與零維內(nèi)彈道理論得到的燃燒室壓強(qiáng)值進(jìn)行比較，并對后者進(jìn)行修正，得到鋁冰發(fā)動機(jī)零維燃燒室壓強(qiáng)pc，理的修正式如下：

pc，修正=0.7055pc，理

(7)

實(shí)驗(yàn)前可用該公式修正鋁冰發(fā)動機(jī)的零維燃燒室壓強(qiáng)計(jì)算結(jié)果，這樣就在壓強(qiáng)計(jì)算中近似考慮了推進(jìn)劑燃燒和兩相燃燒產(chǎn)物的流動過程對壓強(qiáng)的影響。對于同類鋁冰發(fā)動機(jī)，該公式有望提高壓強(qiáng)計(jì)算的準(zhǔn)確度并降低實(shí)驗(yàn)失敗的機(jī)率。

4 結(jié)論

(1)從鋁冰發(fā)動機(jī)的實(shí)驗(yàn)曲線和數(shù)值模擬結(jié)果可以看出，納米鋁的反應(yīng)活性很強(qiáng)，鋁和水蒸氣在推進(jìn)劑燃面附近就能迅速反應(yīng)且放出大量熱量，鋁冰發(fā)動機(jī)能夠正常工作。

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

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

[1] 龐愛民,黎小平.固體推進(jìn)劑技術(shù)的創(chuàng)新與發(fā)展規(guī)律[J].含能材料,2015,23(1):3-6.

[2] 胡潤芝,張小平.高能推進(jìn)劑燃燒效率研究和實(shí)測比沖預(yù)估[J].推進(jìn)技術(shù),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ìn)劑中鋁燃燒理論研究展望[J].飛航導(dǎo)彈,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ìn)劑中的應(yīng)用進(jìn)展[J].兵工學(xué)報(bào),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] 韓秀杰.鋁-冰固體推進(jìn)劑燃燒性能研究[D].哈爾濱:哈爾濱工程大學(xué),2013.

[13] 李宜敏,趙元修.固體火箭發(fā)動機(jī)原理[M]. 北京:國防工業(yè)出版社,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] 陳軍.固體火箭發(fā)動機(jī)零維兩相內(nèi)彈道研究[J].彈道學(xué)報(bào),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] 方丁酉.兩相流動力學(xué)[M].長沙:國防科技大學(xué)出版社,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—)，男，副教授，研究方向?yàn)榻饘偃剂习l(fā)動機(jī)。E-mail：liupingan631@126.com

V435

A

1006-2793(2017)04-0425-07

10.7673/j.issn.1006-2793.2017.04.005