長(zhǎng)時(shí)間小過(guò)載情況下固體火箭發(fā)動(dòng)機(jī)兩相流流場(chǎng)數(shù)值模擬①

那旭東，夏智勛，馬立坤，顏小婷

(國(guó)防科技大學(xué) 空天科學(xué)學(xué)院，長(zhǎng)沙 410073)

0 引言

近年來(lái)，隨著反導(dǎo)攔截能力需求的不斷提升，對(duì)高機(jī)動(dòng)性能導(dǎo)彈發(fā)動(dòng)機(jī)在過(guò)載情況下的可靠性和安全性均提出了更高要求。該類(lèi)導(dǎo)彈的動(dòng)力系統(tǒng)基本上采用高機(jī)動(dòng)性能的固體火箭發(fā)動(dòng)機(jī)，如美國(guó)的PAC-3、俄羅斯S-300V等，其主要技術(shù)特點(diǎn)是導(dǎo)彈在高機(jī)動(dòng)運(yùn)動(dòng)過(guò)程中，要承受長(zhǎng)時(shí)間小過(guò)載、中時(shí)間中過(guò)載、短時(shí)間大過(guò)載等過(guò)載情況[1-2]。高機(jī)動(dòng)性能的固體火箭發(fā)動(dòng)機(jī)大多采用含鋁復(fù)合固體推進(jìn)劑[3-4]，其燃燒產(chǎn)物中含有大量凝相顆粒Al2O3，發(fā)動(dòng)機(jī)內(nèi)流場(chǎng)是典型的氣固兩相流場(chǎng)。在過(guò)載條件下，高速凝相顆粒運(yùn)動(dòng)會(huì)引起燃燒室和噴管的局部聚集和沖刷，使局部燒蝕加劇，對(duì)發(fā)動(dòng)機(jī)熱防護(hù)結(jié)構(gòu)產(chǎn)生很大影響。嚴(yán)重時(shí)會(huì)造成內(nèi)絕熱防護(hù)失效、殼體燒穿，從而導(dǎo)致飛行失敗[5]。因此，研究和掌握過(guò)載條件下發(fā)動(dòng)機(jī)內(nèi)流場(chǎng)顆粒流的運(yùn)動(dòng)規(guī)律，對(duì)提升高機(jī)動(dòng)性能固體火箭發(fā)動(dòng)機(jī)的安全性及可靠性都具有重要的實(shí)現(xiàn)意義。

關(guān)于過(guò)載情況下固體火箭發(fā)動(dòng)機(jī)內(nèi)流場(chǎng)的兩相流流動(dòng)，國(guó)內(nèi)外都進(jìn)行了相關(guān)研究[6-18]，但對(duì)長(zhǎng)時(shí)間小過(guò)載情況下發(fā)動(dòng)機(jī)燃燒室及噴管的三維二相流流動(dòng)的研究報(bào)道則較少[19-20]。本文以某飛行試驗(yàn)失敗發(fā)動(dòng)機(jī)為研究對(duì)象，針對(duì)其長(zhǎng)時(shí)間小過(guò)載的工作特點(diǎn)，對(duì)燃燒室及潛入式噴管進(jìn)行三維兩相流流場(chǎng)的數(shù)值計(jì)算，在計(jì)算過(guò)程中，考慮了發(fā)動(dòng)機(jī)旋轉(zhuǎn)所帶來(lái)的科氏加速度的影響。通過(guò)仿真計(jì)算分析了發(fā)動(dòng)機(jī)內(nèi)部?jī)上嗔髁鲌?chǎng)，并給出顆粒沉積區(qū)域及對(duì)應(yīng)的顆粒沖刷角度、沖刷速度、顆粒粒徑及濃度分布，旨在研究在長(zhǎng)時(shí)間小過(guò)載條件下發(fā)動(dòng)機(jī)內(nèi)部顆粒的流動(dòng)及沉積情況，為該類(lèi)發(fā)動(dòng)機(jī)后續(xù)的熱防護(hù)設(shè)計(jì)提供指導(dǎo)依據(jù)。

1 計(jì)算模型和方法

1.1 發(fā)動(dòng)機(jī)構(gòu)型及網(wǎng)格劃分

以某飛行試驗(yàn)固體火箭發(fā)動(dòng)機(jī)為研究對(duì)象，該發(fā)動(dòng)機(jī)采用翼柱形裝藥。對(duì)該發(fā)動(dòng)機(jī)工作過(guò)程中3個(gè)典型時(shí)刻t1、t2和t3進(jìn)行定常計(jì)算。其中，t1時(shí)刻為法向過(guò)載最大時(shí)刻；t2時(shí)刻為法向過(guò)載變小的拐點(diǎn)；t3時(shí)刻為裝藥完全燒完時(shí)刻。t1～t3典型時(shí)刻發(fā)動(dòng)機(jī)幾何構(gòu)形如圖1所示，主要包括前后封頭、藥柱、發(fā)動(dòng)機(jī)殼體以及潛入式噴管。

采用六面體結(jié)構(gòu)網(wǎng)格對(duì)發(fā)動(dòng)機(jī)進(jìn)行網(wǎng)格劃分，并對(duì)壁面處及預(yù)估顆粒沉積量大的地方進(jìn)行網(wǎng)格加密。圖2給出了t1時(shí)刻其網(wǎng)格劃分結(jié)果，其中上面為采用結(jié)構(gòu)網(wǎng)格劃分的塊拓?fù)浣Y(jié)構(gòu)，下面為對(duì)應(yīng)的發(fā)動(dòng)機(jī)三維結(jié)構(gòu)網(wǎng)格。

圖1 發(fā)動(dòng)機(jī)幾何模型

1.2 數(shù)值方法及相關(guān)參數(shù)

采用Fluent軟件開(kāi)展仿真計(jì)算，兩相流流動(dòng)使用離散相模型DPM(Discrete Phase Model)求解。燃?xì)庾鳛檫B續(xù)相在歐拉坐標(biāo)系下求解，顆粒相作為離散相進(jìn)行處理，對(duì)每個(gè)顆粒在拉格朗日坐標(biāo)下進(jìn)行跟蹤。燃?xì)獾目刂品匠滩捎肗-S方程，顆粒相的運(yùn)動(dòng)軌跡則由拉格朗日方程確定[17]。

圖2 t1時(shí)刻發(fā)動(dòng)機(jī)塊拓?fù)浣Y(jié)構(gòu)及網(wǎng)格劃分

鑒于本文主要研究發(fā)動(dòng)機(jī)中氣相及顆粒相在長(zhǎng)時(shí)間小過(guò)載條件下的流動(dòng)特性，針對(duì)過(guò)載條件下發(fā)動(dòng)機(jī)流動(dòng)的特點(diǎn)，將復(fù)雜的兩相流流動(dòng)模型進(jìn)行適當(dāng)簡(jiǎn)化[21]：

(1)不考慮燃面退移過(guò)程，對(duì)選定的時(shí)刻進(jìn)行定常計(jì)算；

(2)不考慮化學(xué)反應(yīng)，將發(fā)動(dòng)機(jī)內(nèi)復(fù)雜的多組分氣體簡(jiǎn)化為單一燃?xì)?，典型時(shí)刻的相關(guān)燃?xì)鈪?shù)見(jiàn)表1；

表1 發(fā)動(dòng)機(jī)燃?xì)鈪?shù)

(3)燃?xì)饪刂品匠滩捎萌S可壓粘性N-S方程，采用標(biāo)準(zhǔn)k-ε湍流模型；

(4)固相顆粒為球形Al2O3，不考慮顆粒相的燃燒、蒸發(fā)、破碎過(guò)程以及顆粒間撞作用，忽略顆粒相變及密度變化的影響；

(5)忽略過(guò)載對(duì)氣相的影響，對(duì)固相的作用通過(guò)UDF施加體積力給出。

顆粒直徑分布采用Rosin-Rammler分布[18]，即假定顆粒直徑d與大于直徑的顆粒質(zhì)量分?jǐn)?shù)Yd之間有如下指數(shù)關(guān)系：

(1)

結(jié)合發(fā)動(dòng)機(jī)所采用推進(jìn)劑配方和文獻(xiàn)[20, 22-23]，選取粒度分布范圍為10～100 μm，顆粒直徑分布指數(shù)為3.5，粒徑平均直徑為55 μm，顆粒的直徑分布曲線(xiàn)如圖3所示。

由于燃燒溫度(見(jiàn)表1)明顯高于Al2O3熔點(diǎn)溫度2318 K ，所以固相顆粒的參數(shù)按液相Al2O3物性參數(shù)給定[24-27]，其中密度為2700 kg/m3、比定壓熱容1883 J/(kg·K)、熱導(dǎo)率7.4 W/(m·K)。

顆粒壁面碰撞模型對(duì)模擬固體火箭發(fā)動(dòng)機(jī)殼體內(nèi)顆粒的沉積分布至關(guān)重要。Mundo等[28-29]基于大量實(shí)驗(yàn)認(rèn)為，液滴撞擊到壁面后的狀態(tài)主要由液滴的動(dòng)能決定，并提出了以Sommerfeld數(shù)K=We0.5·Re0.25的大小為判定準(zhǔn)則，當(dāng)K≤57.7時(shí)，液滴粘附于絕熱材料壁面；反之，液滴在絕熱材料表面發(fā)生反彈飛濺，其中We與Re分別為入射液滴的韋伯?dāng)?shù)與雷諾數(shù)[30-31]。

本文采用以Sommerfeld數(shù)為判定準(zhǔn)則的顆粒壁面碰撞模型，對(duì)固體火箭發(fā)動(dòng)機(jī)內(nèi)顆粒碰撞壁面后狀態(tài)進(jìn)行判定。

1.3 邊界條件處理

計(jì)算過(guò)程中主要用到三類(lèi)邊界條件：加質(zhì)面邊界條件、固壁面邊界條件及出口邊界條件。

(1)加質(zhì)面邊界：根據(jù)發(fā)動(dòng)機(jī)實(shí)際設(shè)計(jì)情況給出當(dāng)前時(shí)刻推進(jìn)劑表面的質(zhì)量流量，其中凝相占總質(zhì)量流量的30%。顆粒從加質(zhì)壁面均勻加入，在推進(jìn)劑表面上應(yīng)用反彈模型，在前封頭、后封頭及噴管固壁面上通過(guò)Sommerfeld數(shù)判定顆粒在壁面處的行為，在噴管出口應(yīng)用消匿模型。

(2)固壁面邊界條件：采用無(wú)滑移條件，即認(rèn)為壁面處流體速度與該處壁面速度相同，法向壓力、溫度梯度為0。

(3)出口邊界條件：發(fā)動(dòng)機(jī)噴管出口氣流處于超聲速狀態(tài)，給定出口處背壓。

1.4 計(jì)算策略及驗(yàn)證方法

計(jì)算過(guò)程中，為節(jié)約計(jì)算時(shí)間，首先進(jìn)行燃?xì)饬鲌?chǎng)穩(wěn)態(tài)計(jì)算，燃?xì)饬鲌?chǎng)穩(wěn)定后再添加顆粒項(xiàng)進(jìn)行兩相流場(chǎng)計(jì)算。連續(xù)相每計(jì)算200次進(jìn)行一次顆粒的計(jì)算，顆粒計(jì)算的迭代步數(shù)為2000步。計(jì)算結(jié)果與表1中噴管喉部及噴管出口的壓強(qiáng)、溫度、馬赫數(shù)等參量進(jìn)行對(duì)比，驗(yàn)證數(shù)值結(jié)果的可靠性，并對(duì)網(wǎng)格無(wú)關(guān)性進(jìn)行驗(yàn)證。

2 彈體坐標(biāo)系及顆粒受力分析

2.1 彈體坐標(biāo)系

彈體坐標(biāo)系見(jiàn)圖4。

圖4 彈體坐標(biāo)系

圖4 中，坐標(biāo)系遵循右手定則，X軸為彈體軸向方向，指向彈頭方向?yàn)檎?，Y、Z軸垂直于X軸，Y軸位于Ⅲ軸上，Z軸位于Ⅳ軸上。軸向過(guò)載方向?yàn)閄軸方向，法向過(guò)載方向?yàn)閅軸方向，橫向過(guò)載方向?yàn)閆軸方向。

2.2 顆粒受力分析

在導(dǎo)彈飛行過(guò)程中，存在繞Z軸的角速度ω。當(dāng)牽連運(yùn)動(dòng)為轉(zhuǎn)動(dòng)時(shí)，動(dòng)點(diǎn)在某瞬時(shí)的絕對(duì)加速度等于該瞬時(shí)它的相對(duì)加速度、牽連加速度與科氏加速度的矢量和[32]。因此，發(fā)動(dòng)機(jī)中的固體顆粒除了存在橫向、軸向及法向的過(guò)載加速度外，還應(yīng)考慮顆粒所受的科氏加速度。

顆粒的絕對(duì)加速度為

a絕=a相+a牽+a科

式中a絕為質(zhì)點(diǎn)在慣性系中的加速度。

牽連加速度為

a牽=α×r+ω×(ω×r)

式中α為角加速度；r為顆粒位置矢量。

科氏加速度為

a科=2ω×v′

式中v′為顆粒相對(duì)于動(dòng)坐標(biāo)系的相對(duì)速度。

假定發(fā)動(dòng)機(jī)3個(gè)方向的過(guò)載加速度分別為a軸、a法和a橫，發(fā)動(dòng)機(jī)繞Z軸的角速度為ω，忽略角加速度α，則有：

相對(duì)加速度為

a相=-(a軸+a法+a橫)

牽相加速度為

a牽=ω×(ω×r)

科氏加速度為

a科=2ω×v′

將慣性坐標(biāo)系設(shè)定在發(fā)動(dòng)機(jī)上，在慣性坐標(biāo)系中研究發(fā)動(dòng)機(jī)中固相顆粒加速度，則固相顆粒所受的加速度為

a=a相+a牽+a科

=-a軸-a法-a橫+

ω×(ω×r)+2ω×v′

可得絕對(duì)加速度在慣性直角坐標(biāo)系下的3個(gè)分量分別為

ax= -a軸+[ωj(ωirj-ωjri)-ωk(ωirk-ωkri)]i+

ay= -a法+[ωi(ωirj-ωjri)-ωk(ωjrk-ωkrj)]j+

az= -a橫+[ωi(ωirk-ωkri)-ωj(ωjrk-ωkrj)]k+

由于導(dǎo)彈只存在Z軸方向上的角速度，因此ωi=0，ωj=0，所以絕對(duì)加速度在發(fā)動(dòng)機(jī)慣性直角坐標(biāo)系下的3個(gè)分量最終表達(dá)式為

az=-a橫

通過(guò)UDF的宏DEFINE_DPM_BODY_FORCE在計(jì)算中對(duì)顆粒施加過(guò)載加速度。所有變量單位采用國(guó)際單位制，角速度的單位為rad/s，加速度的單位為m/s2，速度的單位為m/s，顆粒位置矢量單位為m。t1～t3典型時(shí)刻的過(guò)載數(shù)值見(jiàn)表2。

表2 過(guò)載數(shù)值

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

對(duì)該發(fā)動(dòng)機(jī)在3個(gè)典型時(shí)刻下分別進(jìn)行定常計(jì)算。圖5分別為t1～t3時(shí)刻，Z=0平面處發(fā)動(dòng)機(jī)內(nèi)流場(chǎng)的馬赫數(shù)、速度、壓強(qiáng)及溫度場(chǎng)計(jì)算結(jié)果。將圖5中的數(shù)據(jù)與表1中噴管喉部及噴管出口的馬赫數(shù)、壓強(qiáng)和溫度的數(shù)值對(duì)比，其結(jié)果基本一致。圖5中，噴管二相流流場(chǎng)的速度和溫度分布受顆粒的影響較大，而壓力場(chǎng)受顆粒的影響較小。從速度場(chǎng)可看出，噴管軸心處的速度要小于接近壁面處的速度，而軸心處的溫度要高于接近壁面處的溫度，這是由于顆粒的滯后性造成的，是典型的氣固顆粒流流場(chǎng)特征[33]。

(a)t1時(shí)刻

(b)t2時(shí)刻

(c)t3時(shí)刻

圖6給出了3個(gè)典型時(shí)刻的顆粒軌跡線(xiàn)分布情況。從圖6可看出，顆粒從推進(jìn)劑表面進(jìn)入燃燒室后，在燃?xì)鈿饬鞯膸?dòng)下迅速向發(fā)動(dòng)機(jī)軸線(xiàn)方向靠近，與發(fā)動(dòng)機(jī)壁面幾乎不發(fā)生碰撞，而在離后封頭及噴管較近的某些位置進(jìn)入流場(chǎng)的顆粒比較容易撞擊壁面，從而形成沉積。此外，從噴管處顆粒的軌跡可看出，大尺寸的顆?；狙貒姽茌S線(xiàn)附近區(qū)域流動(dòng)，而小尺寸的顆粒遠(yuǎn)離噴管軸線(xiàn)區(qū)域流動(dòng)。這主要是小尺寸顆粒的隨流性較好，在氣流帶動(dòng)作用下軌跡發(fā)生了偏轉(zhuǎn)。

圖7給出了t1時(shí)刻發(fā)動(dòng)機(jī)顆粒沉積濃度分布情況?？煽闯?，濃度分布較大的位置主要集中在后封頭圓弧過(guò)渡段及潛入式噴管潛入段壁面處。

圖6 顆粒軌跡線(xiàn)

圖7 t1時(shí)刻顆粒沉積濃度分布

圖8為t1、t2和t3時(shí)刻發(fā)動(dòng)機(jī)后封頭及噴管處顆粒沉積濃度分布放大圖。其中，左圖為從發(fā)動(dòng)機(jī)外部視角觀(guān)測(cè)顆粒沉積濃度分布；右圖為從發(fā)動(dòng)機(jī)內(nèi)部視角觀(guān)測(cè)后封頭及潛入式噴管濃度分布。

從圖8(a)可看出，后封頭處的顆粒沉積濃度要大于噴管處，主要集中于后封頭圓弧過(guò)渡段處，且在承載方向上沉積濃度發(fā)生偏轉(zhuǎn)。從圖8(b)可看出，t2時(shí)刻其后封頭圓弧過(guò)渡段處顆粒沉積濃度分布沒(méi)有發(fā)生明顯的偏轉(zhuǎn)，但其仍為顆粒沉積濃度分布最大處，在噴管潛入段及與其相鄰的后封頭平面處顆粒沉積濃度在承載方向上發(fā)生了偏轉(zhuǎn)。從圖8(c)可看出，t3時(shí)刻其顆粒沉積濃度在后封頭圓弧過(guò)渡段處和噴管潛入段處在承載方向上發(fā)生了偏轉(zhuǎn)，且顆粒沉積濃度最大處仍處于圓弧過(guò)渡段處。

針對(duì)t1時(shí)刻的工況，除了考慮所有過(guò)載情況下的計(jì)算，還進(jìn)行了無(wú)過(guò)載計(jì)算及無(wú)橫向過(guò)載計(jì)算。圖9與圖10分別為無(wú)過(guò)載及無(wú)橫向過(guò)載條件下后封頭及噴管處顆粒沉積濃度分布。分析對(duì)比圖8(a)、圖9和圖10可知，對(duì)于后封頭圓弧過(guò)渡段，當(dāng)無(wú)過(guò)載時(shí)顆粒的沉積在周向方向上分布相對(duì)均勻，無(wú)明顯聚集偏轉(zhuǎn)，最大沉積量約為23 kg/m3；當(dāng)有過(guò)載時(shí)沉積的分布發(fā)生了明顯的偏轉(zhuǎn)，在承載方向上沉積量有所增加，考慮所有過(guò)載情況下時(shí)的最大沉積量約為42.5 kg/m3，約為無(wú)過(guò)載情況下的1.89倍；無(wú)橫向過(guò)載(0.5g)情況下的最大沉積量約為42.8 kg/m3，兩者之間差距并不明顯。

(c) t3時(shí)刻

由于顆粒的沖刷速度與沖刷角度對(duì)研究顆粒對(duì)絕熱層的沖擊燒蝕具有重要意義，因此對(duì)圖8(b)、(c)及圖9中后封頭圓弧過(guò)渡段處顆粒的沖刷速度及角度進(jìn)行統(tǒng)計(jì)，結(jié)果見(jiàn)圖11和圖12。從圖12可看出，顆粒在后封頭圓弧過(guò)渡段處的入射角度大多在20°～40°之間，屬于切向的剪切沖蝕。

此外，還對(duì)后封頭圓弧過(guò)渡段處顆粒尺寸進(jìn)行了統(tǒng)計(jì)，如圖13所示?？煽闯觯蠓忸^圓弧過(guò)渡段處沉積的顆粒主要為相對(duì)較大尺寸的顆粒。

圖9 t1時(shí)刻無(wú)過(guò)載條件下后封頭及噴管顆粒沉積濃度分布

圖10 t1時(shí)刻無(wú)橫向過(guò)載條件下后封頭及噴管顆粒沉積濃度分布

圖11 后封頭圓弧過(guò)渡段處顆粒的沖刷速度分布

圖12 后封頭圓弧過(guò)渡段處顆粒的沖刷角度分布

圖13 后封頭圓弧過(guò)渡段處顆粒的粒徑分布

從以上仿真結(jié)果可知，在進(jìn)行此類(lèi)發(fā)動(dòng)機(jī)絕熱層設(shè)計(jì)時(shí)，要充分考慮發(fā)動(dòng)機(jī)在長(zhǎng)時(shí)間小過(guò)載情況下所帶來(lái)的影響，對(duì)顆粒沉積量大的地方，要加大絕熱層的厚度，以防將發(fā)動(dòng)機(jī)殼體燒穿，造成發(fā)動(dòng)機(jī)失效。

4 結(jié)論

(1)本文以某飛行失效固體火箭發(fā)動(dòng)機(jī)為研究對(duì)象，針對(duì)其長(zhǎng)時(shí)間小過(guò)載的工作特點(diǎn)對(duì)其三維兩相流內(nèi)流場(chǎng)進(jìn)行數(shù)值模擬，在計(jì)算過(guò)程中考慮了發(fā)動(dòng)機(jī)旋轉(zhuǎn)所帶來(lái)的科氏加速度的影響。分析了兩相流流場(chǎng)特性，并給出了凝相顆粒沉積位置、沉積量及顆粒的沖刷速度、角度和粒徑的統(tǒng)計(jì)分布。

(2)分析3個(gè)典型時(shí)刻下的顆粒沉積分布，指出沉積量大的位置主要集中在后封頭圓弧過(guò)渡段處，且沉積的分布發(fā)生了的偏轉(zhuǎn)，在承載方向上沉積量有所增加。針對(duì)t1時(shí)刻，分析對(duì)比無(wú)過(guò)載、無(wú)橫向過(guò)載以及考慮所有過(guò)載的情況，相比無(wú)過(guò)載情況，考慮所有過(guò)載時(shí)，其最大沉積量為無(wú)過(guò)載情況下的1.89倍。

(3)在此類(lèi)發(fā)動(dòng)機(jī)后續(xù)設(shè)計(jì)過(guò)程中，要充分考慮長(zhǎng)時(shí)間小過(guò)載情況下兩相流流場(chǎng)對(duì)發(fā)動(dòng)機(jī)帶來(lái)的影響，在顆粒沉積量大的地方要進(jìn)行充分的熱防護(hù)。