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

(上海航天動(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é)合的辦法來(lái)求解噴管中兩相流動(dòng)問(wèn)題;Thakre等[6]應(yīng)用歐拉-拉格朗日方法結(jié)合多組分模型研究了HTPB復(fù)合推進(jìn)劑的燃燒產(chǎn)物在石墨噴管中的流動(dòng)與燒蝕問(wèn)題;Greenfield[7]等針對(duì)超聲速微噴管兩相流問(wèn)題,研究了兩相流對(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é)者近年來(lái)逐步開(kāi)展了對(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ì)算方法,開(kāi)展了斜切噴管內(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ù)見(jiàn)表1。

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

1.2 基本假設(shè)

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

1)忽略化學(xué)反應(yīng)過(guò)程,推進(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是由于在可壓縮湍流中,過(guò)渡擴(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ù)信息,最終將顆粒相代入氣相方程再次求解,將此過(guò)程循環(huán)迭代,直至收斂。計(jì)算過(guò)程采用二階迎風(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)特性開(kāi)展研究的方法已被很多學(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)分布,通過(guò)圖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

通過(guò)圖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)?噴管軸線附近顆粒濃度較大,顆粒隨氣體一起流入噴管過(guò)程中,顆粒相速度滯后于氣相,阻礙燃?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)過(guò)喉部以后隨粒徑增大而減小。主要原因是顆粒直徑越小,其慣性越小,經(jīng)過(guò)喉部后氣流對(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í)由于顆粒相存在熱慣性,將保持原來(lái)溫度,當(dāng)顆粒運(yùn)動(dòng)到噴管斜切部分時(shí),速度減小,與周圍燃?xì)鉄崃拷粨Q率增加,使氣相溫度升高。根據(jù)圖7可知由于通過(guò)噴管上下壁面的兩束濃度較大的顆粒流在斜切噴管基本部分?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)、外表面作用力的合力,通過(guò)對(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%工況下,通過(guò)表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)角減小。