張 俊,高天宇,高璞清,王 剛

(1 上海海洋大學(xué)工程學(xué)院，上海 201306； 2 中國(guó)航天科工集團(tuán)第六研究院第41研究所，呼和浩特 010010)

0 引言

固體火箭發(fā)動(dòng)機(jī)的工作過程按照其燃燒室壓強(qiáng)變化可大致分為點(diǎn)火建壓期、穩(wěn)定工作期和工作末期。在發(fā)動(dòng)機(jī)點(diǎn)火使其燃燒室達(dá)到一定工作壓強(qiáng)后,進(jìn)入持續(xù)時(shí)間相對(duì)較長(zhǎng)的穩(wěn)定工作期。隨著推進(jìn)劑燃面的推移、肉厚不斷減少,襯層和絕熱層在高速高溫燃?xì)馍淞鞯臒g與沖刷作用下可能出現(xiàn)剝離、脫落、燒穿等現(xiàn)象。尤其對(duì)于含金屬的復(fù)合推進(jìn)劑而言,產(chǎn)生的高速高溫顆粒受到發(fā)動(dòng)機(jī)的結(jié)構(gòu)限制和飛行過載等條件的作用,將會(huì)使發(fā)動(dòng)機(jī)后封頭及噴管收斂段等部位的粒子沉積濃度升高明顯,進(jìn)一步加劇了這些部位的燒蝕。因此,研究SRM的工作末期的多相內(nèi)流場(chǎng)變化,確定流場(chǎng)中固相顆粒的分布規(guī)律、運(yùn)動(dòng)軌跡及沉積濃度等是一個(gè)必不可少的環(huán)節(jié)[1-3]。

目前,科研人員主要通過數(shù)值模擬研究這一問題,取得了一系列重要的研究成果[4-6]。文獻(xiàn)[7]對(duì)發(fā)動(dòng)機(jī)的燃燒室進(jìn)行冷流實(shí)驗(yàn),對(duì)徑向速度、脈動(dòng)速度等物理量進(jìn)行了實(shí)驗(yàn)測(cè)量;文獻(xiàn)[8]將SSTk-ω湍流模型及其改進(jìn)形式用于燃燒室湍流流場(chǎng)的數(shù)值模擬,并將計(jì)算結(jié)果與Wilcox和S-A湍流模型進(jìn)行了對(duì)比,得出SSTk-ω計(jì)算的燃燒室徑向速度分布與實(shí)驗(yàn)值一致性最好的結(jié)論;文獻(xiàn)[9]對(duì)凝相粒子加速度進(jìn)行了理論分析、數(shù)值計(jì)算及縮比發(fā)動(dòng)機(jī)旋轉(zhuǎn)模擬試驗(yàn),分析了地面模擬過載試驗(yàn)中燃燒室粒子聚集區(qū)及旋轉(zhuǎn)試驗(yàn)后的燒蝕部位;文獻(xiàn)[10-11]對(duì)高含鋁固體推進(jìn)劑低壓SRM尾流場(chǎng)復(fù)燃進(jìn)行了數(shù)值模擬和實(shí)驗(yàn)研究,得到了鋁粉引入、粒徑和燃燒室壓強(qiáng)對(duì)羽流溫度的影響規(guī)律,并與地面試車結(jié)果進(jìn)行了對(duì)比。

綜上所述,針對(duì)SRM工作末期、飛行過載條件下的三維兩相內(nèi)流場(chǎng)研究較少,尤其是在復(fù)合過載條件下的粒子沉積部位和運(yùn)動(dòng)規(guī)律有待進(jìn)一步深入研究。文中對(duì)比研究了SRM在兩個(gè)工作時(shí)刻下的內(nèi)流場(chǎng)壓強(qiáng)、溫度、速度、粒子濃度及粒子運(yùn)動(dòng)軌跡的分布,重點(diǎn)研究了飛行過載對(duì)粒子運(yùn)動(dòng)規(guī)律及沉積濃度的影響規(guī)律。

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

1.1 控制方程

在含金屬?gòu)?fù)合推進(jìn)劑的燃燒產(chǎn)物中,粒徑分布呈對(duì)數(shù)正態(tài)分布或雙峰分布。由于不同尺寸的顆粒有不同的控制方程組,如果按顆粒尺寸連續(xù)分布計(jì)算兩相流流場(chǎng),則顆粒相的控制方程過多而無法求解。一般處理方法是用多個(gè)離散的顆粒群代替真實(shí)的連續(xù)性分布,計(jì)算中不考慮顆粒相的燃燒、蒸發(fā)、破碎,不考慮兩相之間的組分變化、熱化學(xué)反應(yīng)。粒子從藥柱燃面拋出后,采用顆粒軌道模型計(jì)算粒子的慣性力、拖曳力和重力。

三維非定常可壓縮流動(dòng)的強(qiáng)守恒型N-S方程在直角坐標(biāo)系下為:

(1)

式中:U=(ρ,ρu,ρv,ρw,ρe)-1;E、F、G為矢通量;Ev、Fv、Gv為粘通量;引入加速度載荷及氣固兩相流動(dòng)的動(dòng)量交換后,源項(xiàng)Q為:

(2)

式中:ρ為燃?xì)庀嗝芏?V為燃?xì)庀嗨俣?Vp為粒子速度;u、v、w為3個(gè)方向上的速度分量;e為單位體積的總能量;gx、gy為加速度分量;?為網(wǎng)格體積;Np表示一個(gè)計(jì)算粒子所代表的物理粒子數(shù)目;下標(biāo)“p”表示顆粒相。

1.2 顆粒軌道模型

在顆粒軌道模型中,粒子接觸到絕熱層后不考慮其能量吸收和損耗,滿足動(dòng)量守恒定理:

(3)

式中:Fp、Fg、Fo分別為粒子受到的拖曳力、重力和過載力。Fpx、Fpy為拖曳力在x,y方向上的分量。

過載力和拖曳力為:

(4)

拖曳系數(shù)為:

(5)

積分得到粒子運(yùn)動(dòng)軌跡為:

(6)

1.3 湍流模型

SSTk-ω湍流模型屬于一種積分到壁面的兩方程渦粘性模型,能夠適用于逆壓梯度變化大的多種物理現(xiàn)象,可應(yīng)用于粘性內(nèi)層。對(duì)于雷諾數(shù)變化范圍較大,壓力梯度大的情況可提高計(jì)算精度和收斂速度。文獻(xiàn)[8]將多種湍流模型的計(jì)算結(jié)果和實(shí)驗(yàn)數(shù)據(jù)進(jìn)行了對(duì)比分析,結(jié)果表明,SSTk-ω湍流模型在SRM多相流場(chǎng)中的數(shù)值計(jì)算結(jié)果與實(shí)驗(yàn)值吻合最好。SSTk-ω湍流模型表示為:

φ3=F1φ1+(1-F1)φ2

(7)

式中:φ1代表標(biāo)準(zhǔn)k-ω兩方程湍流模型;φ2代表變形后的k-ε兩方程湍流模型;模式函數(shù)F1的計(jì)算參見文獻(xiàn)[8]。

2 數(shù)值計(jì)算模型

2.1 網(wǎng)格劃分及邊界條件

為充分考慮固相粒子對(duì)絕熱結(jié)構(gòu)的影響,選擇SRM工作末期的2個(gè)工況作為研究對(duì)象,發(fā)動(dòng)機(jī)工作時(shí)間為42 s,工況1為第38 s,工況2為第40 s。內(nèi)流場(chǎng)計(jì)算區(qū)域?yàn)榘l(fā)動(dòng)機(jī)外輪廓封閉空間去除藥柱結(jié)構(gòu)部分的區(qū)域。圖1為SRM工作末期兩個(gè)工況下的1/2三維幾何模型及邊界條件。

圖1 幾何模型及邊界條件

流場(chǎng)的計(jì)算域采用六面體結(jié)構(gòu)化網(wǎng)格劃分方法,對(duì)流場(chǎng)變化劇烈的區(qū)域進(jìn)行局部網(wǎng)格細(xì)化處理,在網(wǎng)格交接處采用網(wǎng)格交接面,在貼近壁面區(qū)域添加邊界層網(wǎng)格。圖2(a)、圖2(b)為內(nèi)流場(chǎng)的局部網(wǎng)格劃分結(jié)果。藥柱初始燃面為質(zhì)量流率入口邊界,包括氣相質(zhì)量流率入口和顆粒相質(zhì)量流率入口,噴管出口為壓力出口邊界,流場(chǎng)對(duì)稱面為軸對(duì)稱邊界,其余外邊界面為標(biāo)準(zhǔn)絕熱壁面邊界。

圖2 局部網(wǎng)格

2.2 計(jì)算方法及數(shù)據(jù)

離散相介質(zhì)的流入方向?yàn)榧淤|(zhì)面法線方向,并用完全氣體狀態(tài)方程封閉氣相計(jì)算模型。為獲得顆粒相和燃?xì)庀嗟鸟詈献饔?先對(duì)連續(xù)燃?xì)饬鲌?chǎng)計(jì)算得到收斂解,然后再加入顆粒相,并在一定時(shí)間步長(zhǎng)內(nèi)計(jì)算顆粒的速度及其運(yùn)動(dòng)規(guī)律,得到顆粒運(yùn)動(dòng)引起的質(zhì)量、動(dòng)量及能量通量。利用PSIC方法對(duì)兩相流場(chǎng)反復(fù)迭代計(jì)算直至獲得收斂解。

計(jì)算數(shù)據(jù):推進(jìn)劑密度1 795 kg/m3;燃?xì)舛▔罕葻? 900 J/(kg·K);燃燒室壓強(qiáng)7.9 MPa;燃?xì)饽栙|(zhì)量0.02 kg/mol;燃?xì)鉄釋?dǎo)率0.02 W/(m·K);燃?xì)馄骄ば韵禂?shù)5×10-5kg/(m·s);顆粒密度3 013 kg/m3;顆粒直徑40 μm;顆粒比熱1 437 J/(kg·K);顆粒質(zhì)量流率48.5 kg/s;氣相平均摩爾質(zhì)量0.02 kg/mol;顆粒平均摩爾質(zhì)量0.308 kg/mol;環(huán)境壓強(qiáng)35 651.5 Pa。飛行過載條件:軸向5g、橫向-10g。

3 數(shù)值計(jì)算結(jié)果分析

3.1 壓強(qiáng)、溫度及速度分布

圖3～圖5分別為SRM內(nèi)部對(duì)稱面上的壓強(qiáng)、溫度及速度分布云圖。圖中的(a)、(b)分別為工況1模型和工況2模型的計(jì)算結(jié)果。

圖3 兩種工況下的壓強(qiáng)分布

圖4 兩種工況下的溫度分布

圖5 兩種工況下的速度分布

結(jié)果表明:在噴管絕熱壁面上,由于高速燃?xì)馍淞鞯恼承宰铚饔?存在較大的速度變化梯度,噴管喉部壓強(qiáng)、溫度和速度變化梯度最大;兩種工況下的燃燒室壓強(qiáng)、溫度、速度分布規(guī)律相近;在燃?xì)庀嘀屑尤牍滔嗔Ｗ雍?由于顆粒相的粘性阻滯作用,兩相射流的速度增長(zhǎng)幅度小于單相射流,從而使兩相射流的最高速度小于單相射流的最高速度;由于顆粒相的存在,顆粒溫度變化比燃?xì)饩徛?這是由于顆粒相對(duì)燃?xì)庀鄿囟茸兓休^大阻滯作用,使得兩相流的整體溫度比單相流高;發(fā)動(dòng)機(jī)燃燒室最高壓強(qiáng)為7.9 MPa,峰值高溫度接近3 600 K;距發(fā)動(dòng)機(jī)頭部0.3 m范圍內(nèi),燃?xì)馑俣刃∮?0 m/s,距發(fā)動(dòng)機(jī)頭部2 m范圍內(nèi),燃?xì)馑俣刃∮?00 m/s,噴管出口最高速度達(dá)到3 300 m/s。

3.2 復(fù)合過載對(duì)粒子沉積率和運(yùn)動(dòng)規(guī)律的影響

圖6和圖7分別為發(fā)動(dòng)機(jī)后封頭和對(duì)稱面的固相粒子沉積濃度分布云圖。圖中的(a)、(b)分別為無過載和有過載的模擬結(jié)果。

圖6 后封頭的固相粒子沉積濃度分布

圖7 對(duì)稱面的固相粒子沉積濃度分布

表1為工況1在有過載、無過載兩種條件下的燃?xì)庀嗔骶€和固相粒子運(yùn)動(dòng)軌跡。

對(duì)比分析可知:在無過載時(shí),流線和粒子軌跡基本呈軸對(duì)稱分布,但粒子有撞擊壁面后反彈的現(xiàn)象,導(dǎo)致后封頭的粒子分布呈不規(guī)則狀態(tài);在復(fù)合過載條件下,燃?xì)庀嗪皖w粒相的速度呈不對(duì)稱分布規(guī)律,粒子滯后效應(yīng)明顯;在過載反方向聚集大量粒子,造成粒子聚集流及其兩側(cè)尺度較大的低速旋渦流動(dòng)現(xiàn)象;隨著過載反方向一側(cè)的固相粒子濃度增加,最高粒子沉積位置在發(fā)動(dòng)機(jī)后封頭的絕熱層內(nèi)壁;在噴管收斂段和擴(kuò)張段的氣相速度和溫度變化率要遠(yuǎn)大于顆粒相;由于靠近噴管的流場(chǎng)速度較高,該部位燃燒面進(jìn)入流場(chǎng)的顆粒在短時(shí)間內(nèi)就獲得了較高軸向速度,而沿徑向的運(yùn)動(dòng)距離則很短;因此,在慣性力和飛行過載作用下容易撞擊噴管的收斂段,在發(fā)動(dòng)機(jī)絕熱結(jié)構(gòu)設(shè)計(jì)中應(yīng)重點(diǎn)關(guān)注。為定量分析復(fù)合過載條件對(duì)顆粒沉積分布的影響規(guī)律,提取SRM對(duì)稱面的不同偏移距離的面平均粒子沉積濃度,如圖8所示。

表1 不同條件下的流線及粒子運(yùn)動(dòng)軌跡對(duì)比

圖8 不同偏移距離的面上平均粒子沉積濃度

對(duì)比分析可知:在飛行過載條件下,橫向過載反方向一側(cè)的粒子沉積率較高;隨著燃面不斷推進(jìn),在貼近絕熱層壁面位置處的面平均粒子沉積濃度呈升高趨勢(shì);工況1模型最高平均粒子沉積濃度為27.02 kg/m3,工況2模型最高平均粒子沉積濃度為30.2 kg/m3。

4 結(jié)論

通過仿真分析SRM工作末期在飛行過載條件下的氣固兩相內(nèi)流場(chǎng),得到以下結(jié)論:

1)在噴管絕熱壁面處,由于高速燃?xì)饬鞯恼承宰铚饔?存在較大速度梯度;距發(fā)動(dòng)機(jī)頭部0.3 m范圍內(nèi),氣流速度小于10 m/s,距發(fā)動(dòng)機(jī)頭部2 m范圍內(nèi),氣流速度小于100 m/s。噴管出口的氣流速度達(dá)到3 300 m/s。

2)在飛行過載條件下,裝藥表面承載方向集聚大量粒子,造成粒子聚集流及其兩側(cè)尺度較大的低速旋渦流動(dòng)現(xiàn)象,這種旋渦流動(dòng)結(jié)構(gòu)對(duì)顆粒運(yùn)動(dòng)軌跡的影響較大,隨著燃面不斷向絕熱層壁面推移,在貼近壁面位置處的面平均粒子沉積濃度呈升高趨勢(shì);粒子向著橫向過載方向的反方向運(yùn)動(dòng),滯后效應(yīng)明顯,容易撞擊噴管收斂段。