固體火箭發(fā)動(dòng)機(jī)斜切噴管兩相流場(chǎng)特性數(shù)值模擬研究

溫雄飛，趙 瑜，馬新建，趙 志，孟 浩

(上海航天動(dòng)力技術(shù)研究所，上海 201109)

0 引言

為改善固體火箭發(fā)動(dòng)機(jī)工作性能,提高固體推進(jìn)劑能量特性,通常在固體推進(jìn)劑中加入金屬鋁粉[1],從而使推進(jìn)劑燃燒產(chǎn)物中含有氣相與凝相顆粒兩種成分[2],在發(fā)動(dòng)機(jī)燃燒室及噴管內(nèi)形成典型的兩相流動(dòng)。凝相顆粒的存在會(huì)對(duì)發(fā)動(dòng)機(jī)噴管內(nèi)流場(chǎng)特性產(chǎn)生較大影響[3],同時(shí)凝相顆粒流會(huì)沖擊噴管壁面,加速燒蝕而改變噴管型面和喉部面積[4],因此研究固體發(fā)動(dòng)機(jī)噴管中的兩相流動(dòng)情況,能夠?yàn)閲姽茉O(shè)計(jì)與熱防護(hù)提供指導(dǎo)。國(guó)外Hwang等[5]應(yīng)用拉格朗日方法與麥科馬克二步二階差分格式相結(jié)合的辦法來求解噴管中兩相流動(dòng)問題;Thakre等[6]應(yīng)用歐拉-拉格朗日方法結(jié)合多組分模型研究了HTPB復(fù)合推進(jìn)劑的燃燒產(chǎn)物在石墨噴管中的流動(dòng)與燒蝕問題;Greenfield[7]等針對(duì)超聲速微噴管兩相流問題,研究了兩相流對(duì)推力和噴管效率的影響。國(guó)內(nèi)淡林鵬等[8]應(yīng)用顆粒軌道模型研究了長(zhǎng)尾噴管中顆粒與壁面在不同碰撞情況下的粒子運(yùn)動(dòng)軌跡。于勇等[9]應(yīng)用顆粒軌道模型求解顆粒相,結(jié)合高精度空間離散格式,對(duì)噴管中不同尺寸顆粒運(yùn)動(dòng)規(guī)律及顆粒相對(duì)流場(chǎng)的影響進(jìn)行了研究;楊育文等[10]對(duì)高含鋁固體推進(jìn)劑低壓下固體火箭發(fā)動(dòng)機(jī)噴管尾流場(chǎng)進(jìn)行了數(shù)值模擬和實(shí)驗(yàn)研究;童悅等[11]對(duì)擺動(dòng)噴管進(jìn)行了三維兩相數(shù)值模擬研究,重點(diǎn)分析了擺角對(duì)粒子濃度分布與噴管推力系數(shù)的影響。斜切噴管主要應(yīng)用于推力矢量控制及導(dǎo)彈級(jí)間分離[12],非常規(guī)的幾何外形導(dǎo)致其流場(chǎng)特性要比一般直噴管更加復(fù)雜。國(guó)內(nèi)學(xué)者近年來逐步開展了對(duì)斜切噴管的研究工作:劉君等[13]應(yīng)用PHOENICS軟件,對(duì)有無(wú)斜切的雙噴管固體火箭發(fā)動(dòng)機(jī)內(nèi)流場(chǎng)進(jìn)行了二維數(shù)值模擬研究;徐瑋等[14]應(yīng)用有限體積法對(duì)某非軸對(duì)稱斜切噴管三維流場(chǎng)特性進(jìn)行了數(shù)值模擬研究,并分析了斜切角度對(duì)發(fā)動(dòng)機(jī)推力的影響;金賀龍等[15-16]采用非定?？蓧嚎sN-S方程與Realizablek-ε湍流模型相結(jié)合的方法,對(duì)不同角度斜切噴管的流場(chǎng)特性與推力特性進(jìn)行了數(shù)值模擬研究。

綜上,國(guó)內(nèi)外諸多學(xué)者對(duì)于固體火箭發(fā)動(dòng)機(jī)噴管內(nèi)兩相流動(dòng)的研究主要集中于常規(guī)直噴管,而針對(duì)斜切噴管兩相流場(chǎng)特性的研究工作相對(duì)較少[17]。文中基于固體火箭發(fā)動(dòng)機(jī)斜切噴管,運(yùn)用Realizablek-ε湍流模型和顆粒軌道模型,采用單元內(nèi)顆粒源法進(jìn)行兩相耦合計(jì)算方法,開展了斜切噴管內(nèi)三維兩相流動(dòng)數(shù)值模擬計(jì)算,分析了兩相流場(chǎng)特性及對(duì)推力的影響。

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

1.1 物理模型

文中研究的斜切噴管結(jié)構(gòu)如圖1所示,其主要由基本部分和斜切部分組成,基本結(jié)構(gòu)參數(shù)有:基本部分?jǐn)U張半角α,斜切部分?jǐn)U張半角β,斜切角δ,噴管入口處直徑L1,喉部直徑L2,噴管擴(kuò)張段軸線長(zhǎng)度L3,噴管出口截面在中心軸線方向的投影距離L4。由于斜切噴管的非對(duì)稱結(jié)構(gòu),其推力由兩部分組成:一部分為沿噴管軸線方向的軸向力Fx,另一部分為垂直于噴管軸線方向的側(cè)向力Fy,其合力與噴管軸線的夾角為推力偏轉(zhuǎn)角φ。選取噴管喉部直徑d=0.06 m作為特征長(zhǎng)度,該斜切噴管結(jié)構(gòu)參數(shù)見表1。

表1 噴管結(jié)構(gòu)參數(shù)表Table 1 Nozzle structure parameters

1.2 基本假設(shè)

由于固體火箭發(fā)動(dòng)機(jī)工作過程涉及復(fù)雜的物理和化學(xué)反應(yīng),為簡(jiǎn)化計(jì)算,又能夠真實(shí)模擬斜切噴管內(nèi)的兩相燃?xì)饬鲃?dòng)情況,引入如下假設(shè)條件:

1)忽略化學(xué)反應(yīng)過程,推進(jìn)劑燃燒生成的混合氣體滿足理想氣體狀態(tài)方程。

2)不考慮氣相、顆粒相以及固體壁面之間的輻射換熱現(xiàn)象。

3)不考慮兩相流對(duì)噴管型面的燒蝕和顆粒沉積。

1.3 數(shù)學(xué)模型

1.3.1 控制方程

氣相控制方程[18]可表示為:

(1)

式中:A=[ρ,ρu,ρv,ρw,e]T,B,C,D為守恒型通量,BV,CV,DV為粘性通量,Sp為固體顆粒對(duì)氣相產(chǎn)生的源項(xiàng)。

顆粒軌道求解控制方程[18]為(單位質(zhì)量下):

(2)

(3)

式中:vp表示顆粒速度;xp表示顆粒位置;u表示氣相的平均速度;u′表示氣相的瞬時(shí)湍流脈動(dòng)速度;g為重力加速度。FD為顆粒的阻力函數(shù):

(4)

式中:CD、ρp和dp分別表示曳力系數(shù)、顆粒密度和顆粒直徑。Re表示相對(duì)雷諾數(shù),其計(jì)算式為:

(5)

式中:ρ表示氣相密度;μ表示動(dòng)力粘度。

1.3.2 湍流模型

數(shù)值模擬選擇Realizablek-ε湍流模型,k方程和ε方程[19]分別為:

(6)

(7)

式中:Gk是由層流速度梯度引起的湍動(dòng)能k的產(chǎn)生項(xiàng);Gb是由浮力引起的湍動(dòng)能k的產(chǎn)生項(xiàng);YM是由于在可壓縮湍流中,過渡擴(kuò)散而產(chǎn)生的波動(dòng);其他參數(shù)為相關(guān)特征常量和自定義源項(xiàng)。

1.3.3 發(fā)動(dòng)機(jī)推力計(jì)算

固體火箭發(fā)動(dòng)機(jī)推力為發(fā)動(dòng)機(jī)內(nèi)、外表面全部作用力的合力,其計(jì)算公式[3]為:

(8)

式中:Pin是作用于發(fā)動(dòng)機(jī)內(nèi)表面的壓強(qiáng);Pout是作用于發(fā)動(dòng)機(jī)外表面的壓強(qiáng);dAs是表面微元面積。對(duì)于斜切噴管發(fā)動(dòng)機(jī),推力由式(8)對(duì)發(fā)動(dòng)機(jī)內(nèi)外表面的壓強(qiáng)積分得出。

1.4 邊界條件與計(jì)算方法

1)氣相邊界條件:選取鋁粉含量18%的HTPB復(fù)合推進(jìn)劑,采用質(zhì)量流量入口條件,給定入口總溫總壓,流動(dòng)方向垂直于入口邊界,出口采用壓力出口,壁面采用絕熱無(wú)滑移邊界條件。

2)顆粒相邊界條件:在入口處,取每個(gè)網(wǎng)格邊的中點(diǎn)作為顆粒的加入點(diǎn)。給定入口處顆粒的初始速度、溫度和顆粒質(zhì)量流量。經(jīng)大量實(shí)驗(yàn)證實(shí),固體火箭發(fā)動(dòng)機(jī)中顆粒相尺寸分布一般服從對(duì)數(shù)正態(tài)分布或者Rosin-Rammler分布,文中取平均粒徑為70 μm的Rosin-Rammler分布,顆粒與壁面碰撞后反彈,切向和法向恢復(fù)系數(shù)均為0.8,顆粒出口邊界不加任何限制條件,達(dá)到出口即逃逸。

3)采用標(biāo)準(zhǔn)壁面函數(shù)法作為近壁處理方式;采用第一層網(wǎng)格y+值為35,增長(zhǎng)率為1.2,共15層的邊界層處理方式。

4)顆粒相與氣相的耦合計(jì)算采用計(jì)算單元內(nèi)顆粒源項(xiàng)算法(PSIC),即先計(jì)算氣相壓強(qiáng)、速度、溫度等參數(shù),再計(jì)算顆粒軌道以及顆粒各項(xiàng)參數(shù)信息,最終將顆粒相代入氣相方程再次求解,將此過程循環(huán)迭代,直至收斂。計(jì)算過程采用二階迎風(fēng)格式,收斂精度為10-6。

5)邊界條件數(shù)值:氣相質(zhì)量流率為8.685 kg/s;氣相入口壓強(qiáng)為11.68 MPa;氣相入口溫度為3 556.3 K;氣相出口壓強(qiáng)為101 325 Pa;氣相出口溫度為293 K;顆粒相質(zhì)量流率為3.722 kg/s;顆粒相入口速度為0.1 m/s;顆粒溫度為3 253 K。

采用上述Realizablek-ε湍流模型、顆粒軌道模型與單元內(nèi)顆粒源法對(duì)固體火箭發(fā)動(dòng)機(jī)內(nèi)兩相流場(chǎng)特性開展研究的方法已被很多學(xué)者應(yīng)用并取得了良好效果[20-21]。

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

2.1 斜切噴管流場(chǎng)特性分析

首先選取顆粒相平均粒徑dp為70 μm的工況進(jìn)行計(jì)算,得到如圖2～圖4所示噴管對(duì)稱截面的氣相及兩相流場(chǎng)的特征參數(shù)分布,其中X、Y為模型位置坐標(biāo),單位為m。分析圖2～圖4可知,顆粒相的存在會(huì)對(duì)斜切噴管的溫度場(chǎng)和速度場(chǎng)產(chǎn)生重要影響,但對(duì)壓力場(chǎng)影響較小,壓力場(chǎng)分布基本不變。

圖2 斜切噴管溫度場(chǎng)分布云圖Fig.2 Temperature field distribution contour of angle-cut nozzle

圖3 斜切噴管壓力場(chǎng)分布云圖Fig.3 Pressure field distribution contour of angle-cut nozzle

圖4 斜切噴管速度場(chǎng)分布云圖Fig.4 Velocity field distribution contour of angle-cut nozzle

對(duì)比兩相流場(chǎng)與純氣相流場(chǎng)情況下,斜切噴管中心軸線處與壁面附近的流場(chǎng)分布,通過圖5、圖6可以發(fā)現(xiàn),顆粒相對(duì)噴管壁面附近的流場(chǎng)結(jié)構(gòu)影響較小,兩相流場(chǎng)速度略小于純氣相流場(chǎng)速度,而溫度基本相同,但顆粒相對(duì)噴管軸線附近區(qū)域的流場(chǎng)結(jié)構(gòu)影響較大。

圖5 斜切噴管距壁面0.02 m附近流場(chǎng)特征參數(shù)對(duì)比Fig.5 Comparison of flow field characteristic parameter curves along 0.02 m away from the angle-cut nozzle wall

圖6 斜切噴管軸線附近流場(chǎng)特征參數(shù)對(duì)比Fig.6 Comparison of flow field characteristic parameter curves along the angle-cut nozzle axis

通過圖7顆粒濃度分布可知顆粒相在軸線附近分布濃度較大,則噴管軸線附近區(qū)域受到顆粒阻力作用,燃?xì)饬鲃?dòng)速度明顯隨顆粒濃度增大而減小,存在一個(gè)低速區(qū);同時(shí)由于顆粒相對(duì)氣相作用的熱增量,導(dǎo)致噴管軸線附近也存在一個(gè)高溫區(qū),且此流場(chǎng)低速區(qū)和高溫區(qū)都起始于噴管喉部并沿噴管出口方向逐漸擴(kuò)大。

圖7 dp=70 μm粒徑下顆粒濃度分布云圖Fig.7 Particle concentration contour for 70 μm particle diameter

結(jié)合流場(chǎng)顆粒濃度分布,分析該現(xiàn)象原因?yàn)?噴管軸線附近顆粒濃度較大,顆粒隨氣體一起流入噴管過程中,顆粒相速度滯后于氣相,阻礙燃?xì)饬鲃?dòng),從而降低燃?xì)饬魉佟４送?由于顆粒溫度高于氣體溫度,對(duì)氣相有部分熱量傳遞,且顆粒之間互相碰撞會(huì)將部分動(dòng)能轉(zhuǎn)化為熱能,這兩部分熱增量導(dǎo)致該區(qū)域流場(chǎng)溫度升高。

2.2 顆粒直徑對(duì)斜切噴管兩相流場(chǎng)特性的影響

圖8為顆粒相最大粒徑為120 μm,最小粒徑為10 μm,平均粒徑分別為70 μm,60 μm,50 μm,40 μm的Rosin-Rammler分布時(shí),流場(chǎng)內(nèi)顆粒運(yùn)動(dòng)軌跡及粒徑分布。

圖8 顆粒軌跡及粒徑分布對(duì)比Fig.8 Comparison of particle tracks and distribution for different particle diameters

由圖可知在噴管基本部分?jǐn)U張段及斜切部分壁面附近存在一個(gè)無(wú)顆粒區(qū),且隨著平均粒徑的增大,斜切部分無(wú)顆粒區(qū)域范圍變大。分析圖9可知,顆粒運(yùn)動(dòng)速度在經(jīng)過喉部以后隨粒徑增大而減小。主要原因是顆粒直徑越小,其慣性越小,經(jīng)過喉部后氣流對(duì)顆粒加速效果越顯著,體現(xiàn)出小顆粒較好的隨流性,在斜切部分的運(yùn)動(dòng)軌跡分布更加分散,無(wú)顆粒區(qū)范圍較小;當(dāng)顆粒直徑變大時(shí),慣性增大導(dǎo)致顆粒運(yùn)動(dòng)受到氣流擾動(dòng)變小,運(yùn)動(dòng)軌跡更集中于軸線附近且比較平直,從而使無(wú)顆粒區(qū)域的范圍變大。

圖9 不同粒徑下噴管軸線顆粒速度分布對(duì)比Fig.9 Comparison of particle velocity along the angle-cut nozzle axis for different particle diameters

圖10給出純氣相以及不同粒徑時(shí),沿出口方向斜切噴管軸線處氣相溫度及速度分布的對(duì)比曲線。

圖10 不同粒徑下噴管軸線附近流場(chǎng)特征參數(shù)對(duì)比Fig.10 Comparison of flow field characteristic parameter curves along the angle-cut nozzle axis for different particle diameters

分析曲線圖可知,顆粒相的存在導(dǎo)致噴管中氣相溫度明顯升高,噴管軸線處氣相溫度先減小后增大,即在斜切噴管基本部分?jǐn)U張段起始處存在一個(gè)低溫區(qū),粒徑越大溫度變化越小,而在斜切部分增大到與收斂段相同的溫度。顆粒相的存在使噴管軸線處兩相流場(chǎng)的氣相速度小于純氣相流場(chǎng)速度,且噴管軸線處氣相速度先增大后減小,存在一個(gè)高速區(qū),粒徑越大速度變化越小。

分析原因認(rèn)為,粒徑越大,慣性越大,對(duì)氣相阻力作用越大,使氣相速度增大較少,速度變化較小;并且粒徑越大,其速度越小,滯留時(shí)間越長(zhǎng),從而與氣相的作用時(shí)間越長(zhǎng),兩相之間的動(dòng)量與熱量交換率增大,使氣相溫度降低較少,溫度變化較小;同時(shí)由于顆粒相存在熱慣性,將保持原來溫度,當(dāng)顆粒運(yùn)動(dòng)到噴管斜切部分時(shí),速度減小,與周圍燃?xì)鉄崃拷粨Q率增加,使氣相溫度升高。根據(jù)圖7可知由于通過噴管上下壁面的兩束濃度較大的顆粒流在斜切噴管基本部分?jǐn)U張段起始區(qū)域附近交匯,導(dǎo)致在此區(qū)域顆粒相與氣相相互作用增強(qiáng),從而動(dòng)量與熱量相互轉(zhuǎn)化率提升,形成如圖11所示的低溫區(qū)與高速區(qū)。

圖11 dp=40 μm粒徑下兩相流場(chǎng)特征Fig.11 Characteristics of two-phase flow field for 40 μm particle diameter

2.3 兩相流對(duì)斜切噴管推力特性的影響

定義Fx為噴管軸向推力;Fy為垂直于燃燒室軸線方向作用力,即噴管側(cè)向推力;φ為推力偏轉(zhuǎn)角;Pmax為流場(chǎng)最大壓強(qiáng);Pmin為流場(chǎng)最小壓強(qiáng)。根據(jù)固體火箭發(fā)動(dòng)機(jī)原理,發(fā)動(dòng)機(jī)的推力等于發(fā)動(dòng)機(jī)所有內(nèi)、外表面作用力的合力,通過對(duì)發(fā)動(dòng)機(jī)內(nèi)外表面作用力進(jìn)行積分計(jì)算可得推力F數(shù)值,計(jì)算結(jié)果如表2所示。

表2 不同粒徑下斜切噴管推力Table 2 Trusts of the angle-cut nozzle for different particle diameters

在顆粒相質(zhì)量分?jǐn)?shù)為30%工況下,通過表2中推力數(shù)據(jù)可知,顆粒相的存在使得噴管推力增大約1%～2%,且隨著顆粒直徑增大,推力減小,推力偏轉(zhuǎn)角減小。由圖12可知,顆粒相的加入使燃?xì)鉁囟葴?出口處燃?xì)鉁囟壬?出口壓強(qiáng)增大,且隨著粒徑增大溫度和壓強(qiáng)提升幅度減小,計(jì)算所得的推力變化趨勢(shì)與壓強(qiáng)變化趨勢(shì)相同。

圖12 噴管出口中心軸線附近燃?xì)鈪?shù)分布對(duì)比Fig.12 Comparison of gas parameter curves along the angle-cut nozzle outlet axis for different particle diameters

3 結(jié)論

文中進(jìn)行斜切噴管內(nèi)三維兩相流動(dòng)數(shù)值模擬,主要分析噴管內(nèi)兩相流場(chǎng)特性以及粒徑大小變化對(duì)流場(chǎng)與發(fā)動(dòng)機(jī)推力的影響,得到以下結(jié)論:

1)顆粒相的存在會(huì)對(duì)發(fā)動(dòng)機(jī)斜切噴管的流場(chǎng)結(jié)構(gòu)產(chǎn)生較大的影響,尤其是對(duì)噴管軸線附近區(qū)域流場(chǎng)結(jié)構(gòu)影響較大,但對(duì)斜切噴管邊界處流場(chǎng)結(jié)構(gòu)影響較小。噴管軸線附近存在一個(gè)燃?xì)饬鲃?dòng)速度較低,溫度較高的區(qū)域,且此區(qū)域起始于噴管喉部并沿噴管出口方向逐漸擴(kuò)大范圍。

2)斜切噴管兩相流場(chǎng)中,斜切部分的壁面附近存在無(wú)顆粒區(qū)域。顆粒直徑越大,運(yùn)動(dòng)越集中于軸線附近且軌跡比較平直,無(wú)顆粒區(qū)域越大。

3)顆粒相直徑的變化會(huì)對(duì)發(fā)動(dòng)機(jī)斜切噴管流場(chǎng)特征參數(shù)的分布產(chǎn)生影響,顆粒的加入使得沿噴管軸線方向氣相速度降低,溫度升高。隨著顆粒直徑增大,發(fā)動(dòng)機(jī)噴管軸線處氣相速度先增大后減小,燃?xì)鉁囟认葴p小后增大;且顆粒直徑越大,溫度和速度變化越小。

4)兩相流狀態(tài)下的噴管推力大于純氣相推力,且隨著顆粒直徑的增加,斜切噴管推力減小,推力偏轉(zhuǎn)角減小。