水環(huán)境下固體火箭發(fā)動(dòng)機(jī)噴流流場(chǎng)數(shù)值研究①

燕國(guó)軍，權(quán)曉波，王占瑩，魏海鵬

(北京宇航系統(tǒng)工程研究所，北京　100076)

?

水環(huán)境下固體火箭發(fā)動(dòng)機(jī)噴流流場(chǎng)數(shù)值研究①

燕國(guó)軍，權(quán)曉波，王占瑩，魏海鵬

(北京宇航系統(tǒng)工程研究所，北京100076)

為了研究不同環(huán)境介質(zhì)下發(fā)動(dòng)機(jī)噴流流場(chǎng)特征，基于VOF多相流模型和SSTk-ω湍流模型，結(jié)合動(dòng)網(wǎng)格方法，建立了水環(huán)境下固體火箭發(fā)動(dòng)機(jī)運(yùn)動(dòng)過程噴流流場(chǎng)數(shù)值仿真模型。通過數(shù)值仿真分析獲得了航行體尾部附著氣泡環(huán)境下發(fā)動(dòng)機(jī)工作過程推力以及尾部氣泡內(nèi)壓力變化特征，燃?xì)馍淞鲗⑹沟梦膊繗馀莩霈F(xiàn)振蕩加強(qiáng)、再次斷裂等復(fù)雜的非定常流動(dòng)特征。同時(shí)獲得了環(huán)境背壓和環(huán)境介質(zhì)對(duì)發(fā)動(dòng)機(jī)推力的影響規(guī)律，發(fā)現(xiàn)環(huán)境介質(zhì)的不同會(huì)對(duì)發(fā)動(dòng)機(jī)工作初始過程的推力產(chǎn)生一定影響。

固體發(fā)動(dòng)機(jī)；水環(huán)境；氣泡；數(shù)值仿真

0　引言

固體火箭發(fā)動(dòng)機(jī)在水環(huán)境中工作可滿足水下工程設(shè)計(jì)時(shí)動(dòng)力系統(tǒng)需求。對(duì)于采用垂直彈射的水下航行體來說，安裝在航行體尾部的發(fā)動(dòng)機(jī)工作前的外部環(huán)境不僅是水介質(zhì)，還有彈射燃?xì)庑纬傻母街诤叫畜w尾部的附體氣泡，深入了解尾部附體氣泡對(duì)固體火箭發(fā)動(dòng)機(jī)水下工作過程噴流流場(chǎng)的影響具有重要的研究?jī)r(jià)值。

采用數(shù)值仿真或試驗(yàn)手段，仍是當(dāng)前研究發(fā)動(dòng)機(jī)在水環(huán)境中工作的主要方法。數(shù)值仿真方面，有學(xué)者[1-2]通過對(duì)燃?xì)馍淞鬟M(jìn)行等壓球形氣泡假設(shè)開展數(shù)值計(jì)算，但更多學(xué)者[3-8]采用基于雷諾時(shí)均化的均質(zhì)平衡多相流模型對(duì)水下燃?xì)馍淞鬟M(jìn)行數(shù)值仿真，研究了發(fā)動(dòng)機(jī)在水環(huán)境中工作推力特性和流場(chǎng)特征；試驗(yàn)方面，學(xué)者們[9-10]通過搭建地面試驗(yàn)系統(tǒng)研究獲得了發(fā)動(dòng)機(jī)在水環(huán)境中的推力矢量特性和燃?xì)馍淞鳉馀莸纳L(zhǎng)及壓力傳播過程。對(duì)于航行體水下運(yùn)動(dòng)過程燃?xì)馀莸难芯匡@示，基于多相流模型的數(shù)值仿真[11]，或采用空泡獨(dú)立膨脹原理結(jié)合Rayleigh-Plesset方程的方法[12]都能較好地模擬水中燃?xì)馀莸男螒B(tài)和壓力變化過程。但上述工作均未將發(fā)動(dòng)機(jī)噴流和水中航行體尾部附體氣泡環(huán)境聯(lián)合起來研究，而本文研究重點(diǎn)關(guān)注的是水環(huán)境中固體火箭發(fā)動(dòng)機(jī)噴流流場(chǎng)在尾部附體氣泡內(nèi)的工作過程和流場(chǎng)演化特征。

本文采用VOF多相流模型和SSTk-ω湍流模型建立控制方程，運(yùn)用層變動(dòng)網(wǎng)格模型進(jìn)行運(yùn)動(dòng)網(wǎng)格更新，對(duì)水環(huán)境中的固體火箭發(fā)動(dòng)機(jī)噴流流場(chǎng)進(jìn)行了非定常數(shù)值仿真(模擬航行體水中運(yùn)動(dòng)過程)，獲得了噴流流場(chǎng)演化特征，分析了環(huán)境背壓和環(huán)境介質(zhì)對(duì)發(fā)動(dòng)機(jī)推力特性的影響。

1　數(shù)值仿真模型

1.1多相流控制方程

首先，對(duì)流體介質(zhì)進(jìn)行適當(dāng)簡(jiǎn)化，將發(fā)動(dòng)機(jī)尾部附體氣泡內(nèi)燃?xì)狻l(fā)動(dòng)機(jī)燃?xì)夂涂諝庖暈橐环N介質(zhì)，為理想氣體，其屬性按燃?xì)鈪?shù)定義，并定義為主相，另一種介質(zhì)則為不可壓的水。采用雷諾時(shí)均化平衡均質(zhì)多相流VOF模型作為流動(dòng)控制方程，包括連續(xù)性方程、動(dòng)量方程、能量方程、體積分?jǐn)?shù)輸運(yùn)方程和狀態(tài)方程。

連續(xù)性方程：

(1)

動(dòng)量方程：

(2)

能量方程：

(3)

水相體積分?jǐn)?shù)方程：

(4)

狀態(tài)方程：

(5)

ρw=const

(6)

通過SSTk-ω湍流模型[13]對(duì)上述方程進(jìn)行封閉。采用有限體積法和SIMPLEC算法對(duì)控制方程進(jìn)行離散求解。

1.2計(jì)算網(wǎng)格、動(dòng)網(wǎng)格方法和邊界條件

本文計(jì)算模型為航行體彈射出筒后運(yùn)動(dòng)一段距離時(shí)發(fā)動(dòng)機(jī)點(diǎn)火工作，工作初始發(fā)動(dòng)機(jī)所處的航行體尾部為彈射燃?xì)庑纬傻母襟w氣泡，計(jì)算過程僅考慮航行體一維的垂直向上運(yùn)動(dòng)過程。由于本文重點(diǎn)分析發(fā)動(dòng)機(jī)水下工作過程，因此簡(jiǎn)化了對(duì)航行體出筒過程的描述。圖1給出了初始計(jì)算流場(chǎng)的計(jì)算域，包括發(fā)動(dòng)機(jī)、模擬航行體結(jié)構(gòu)和外界計(jì)算域。為適應(yīng)水環(huán)境中發(fā)動(dòng)機(jī)工作需求，位于航行體底部位置的發(fā)動(dòng)機(jī)采用小擴(kuò)張比設(shè)計(jì)，其面積擴(kuò)張比(噴管出口面積Ae/喉部面積At)為5.3。全域均劃分為結(jié)構(gòu)化網(wǎng)格，如圖2所示。并采用滑移網(wǎng)格方法和層變動(dòng)網(wǎng)格模型更新運(yùn)動(dòng)網(wǎng)格，分別在發(fā)動(dòng)機(jī)底部和水域中后部設(shè)置了網(wǎng)格參考面，隨著航行體運(yùn)動(dòng)，參考面位置處將不斷分裂生成新的網(wǎng)格。

圖1　初始計(jì)算域

(a)初始全域網(wǎng)格　　　　(b)發(fā)動(dòng)機(jī)網(wǎng)格

外界計(jì)算域邊界除底部為無滑移壁面外，其他均為壓力出口邊界(靜壓滿足重力場(chǎng)分布)。發(fā)動(dòng)機(jī)壁面及模擬航行體壁面為無滑移壁面邊界。當(dāng)發(fā)動(dòng)機(jī)處于H0(對(duì)應(yīng)時(shí)刻t0，發(fā)動(dòng)機(jī)底部環(huán)境壓力為p0)，發(fā)動(dòng)機(jī)開始工作時(shí)，燃燒室入口由無滑移壁面邊界改為燃?xì)饪倝喝肟谶吔鐥l件。發(fā)動(dòng)機(jī)工作總壓pc采用試驗(yàn)測(cè)量結(jié)果作為輸入，如圖3所示(圖中時(shí)間零點(diǎn)為喉部打開時(shí)刻，pc,max為燃燒室最大壓力，tall為發(fā)動(dòng)機(jī)工作時(shí)長(zhǎng)，工作總溫為Tc)。

圖3　無量綱化燃燒室總壓曲線

1.3推力計(jì)算方法

發(fā)動(dòng)機(jī)推力計(jì)算采用內(nèi)外壁面壓力積分的方法。發(fā)動(dòng)機(jī)推力Fthrust計(jì)算公式如下：

Fthrust=pcAc+∫prdsr+∫ppdsp-paAe

(7)

式中pcAc+∫prdsr為燃燒室內(nèi)壁面推力；pc為燃燒室入口總壓；Ac為燃燒室橫截面積；pr為如圖2中除壓力入口以外燃燒室內(nèi)壁面網(wǎng)格單元上的壓力；dsr為對(duì)應(yīng)網(wǎng)格單元沿推力方向的投影面積；∫ppdsp為噴管內(nèi)壁面推力；pp為噴管內(nèi)壁面網(wǎng)格單元上的壓力；dsp為對(duì)應(yīng)網(wǎng)格沿推力方向的投影面積；-paAe則為發(fā)動(dòng)機(jī)外壁面推力;pa為環(huán)境背壓，Ae為噴管出口面積。

1.4計(jì)算模型驗(yàn)證

以發(fā)動(dòng)機(jī)喉部打開時(shí)刻為時(shí)間零點(diǎn)，仿真獲得的推力計(jì)算結(jié)果和試驗(yàn)結(jié)果的比較如圖4所示(圖中Fthrust,max為非定常計(jì)算獲得的最大推力)，計(jì)算值與試驗(yàn)值吻合良好。以航行體彈射時(shí)刻為零點(diǎn)，航行體底部壓力測(cè)點(diǎn)(測(cè)點(diǎn)位置見圖1)計(jì)算結(jié)果與試驗(yàn)對(duì)比情況如圖5所示，整體來看，呈現(xiàn)出底部壓力隨所處水深位置的減小而降低的變化趨勢(shì)，計(jì)算結(jié)果和試驗(yàn)結(jié)果具有較好的可比擬性，反映出了試驗(yàn)獲得的3個(gè)階段壓力特征(本文重點(diǎn)分析發(fā)動(dòng)機(jī)點(diǎn)火工作后對(duì)尾部附著氣泡流場(chǎng)影響，并忽略了航行體彈射段過程流體動(dòng)力現(xiàn)象，對(duì)尾部附著氣泡生成演化過程不重點(diǎn)描述，詳細(xì)可見文獻(xiàn)[11])：T1階段為航行體出筒后尾部附著氣泡生成演化，并出現(xiàn)壓力周期性振蕩，發(fā)動(dòng)機(jī)工作后T2段的壓力振蕩加強(qiáng)(圖5中T2段相對(duì)于T1段壓力振蕩幅值增大)，以及T3段噴流、氣泡和水環(huán)境相互作用，形成一種較穩(wěn)定的小幅壓力振蕩。

從上述兩方面來看，本文所建立的仿真模型能較準(zhǔn)確地預(yù)示航行體垂直發(fā)射尾空泡條件下發(fā)動(dòng)機(jī)水下工作過程的尾部流動(dòng)特征和推力變化過程。

圖4　無量綱化推力仿真值與試驗(yàn)值比較

圖5　無量綱化航行體底部壓力測(cè)點(diǎn)仿真值與試驗(yàn)值比較

2　尾部流場(chǎng)分析

圖6和圖7利用軸對(duì)稱特征，分別給出了發(fā)動(dòng)機(jī)工作過程典型時(shí)刻航行體尾部壓力-體積分?jǐn)?shù)云圖分布和溫度-馬赫數(shù)云圖分布(圖中t0時(shí)刻為發(fā)動(dòng)機(jī)工作初始時(shí)刻)。T1階段為發(fā)動(dòng)機(jī)工作前航行體(發(fā)動(dòng)機(jī))尾部生成附體氣泡的過程，其流場(chǎng)演化特征與文獻(xiàn)[11]中闡述的過程相似，彈射出筒后航行體尾部附著氣泡與環(huán)境水介質(zhì)相互作用經(jīng)膨脹、收縮和斷裂后形成航行體尾部附體氣泡(見(a)圖)。發(fā)動(dòng)機(jī)工作后的T2段(見(b)～(e)圖)，燃?xì)馍淞髯⑷敫襟w氣泡，不僅抑制了航行體尾部附體氣泡底部回射水流的發(fā)展，且促進(jìn)了附體氣泡體積的增大，進(jìn)而使得附體氣泡的壓力振蕩幅度加大。但由于環(huán)境水的作用，附體氣泡中部再次發(fā)生了收縮和多處斷裂，并卷入水沖擊航行體尾部，形成航行體尾部端面緊鄰噴管出口小范圍的壓力升高，但隨著氣水的相互摻混，沖擊逐漸衰減。發(fā)生斷裂后的附體氣泡體積進(jìn)一步減小，航行體尾部主要以燃?xì)馍淞鹘Y(jié)構(gòu)為主，進(jìn)而進(jìn)入壓力小幅振蕩的T3階段(典型時(shí)刻見(f)圖)，燃?xì)馍淞髋c水相互作用而不斷發(fā)生斷裂脫落。

3　推力特性分析

為了對(duì)航行體尾部附著氣泡條件下固體火箭發(fā)動(dòng)機(jī)水環(huán)境中工作過程推力特性開展進(jìn)一步分析，采用上述發(fā)動(dòng)機(jī)和相同的總壓、總溫輸入條件。固定外界環(huán)境背壓為發(fā)動(dòng)機(jī)工作水深H0的靜壓，構(gòu)造了另外2個(gè)仿真對(duì)比工況，一個(gè)環(huán)境介質(zhì)為空氣，另一個(gè)則為水，并基于相同的數(shù)值計(jì)算方法獲得其仿真結(jié)果，通過與尾部附著氣泡環(huán)境仿真結(jié)果三者之間的對(duì)比，分析環(huán)境背壓和環(huán)境介質(zhì)對(duì)發(fā)動(dòng)機(jī)水下工作過程推力特性的影響。

3.1環(huán)境背壓影響

比較3個(gè)工況發(fā)動(dòng)機(jī)水下工作全過程的推力曲線，如圖8所示。從圖8中可看出，3個(gè)工況的推力變化趨勢(shì)基本一致，與燃燒室的總壓變化趨勢(shì)相同。從量值上來看，固定環(huán)境背壓工況下，環(huán)境水介質(zhì)和氣介質(zhì)的推力量值在發(fā)動(dòng)機(jī)額定工作時(shí)相當(dāng)；對(duì)于尾部附著氣泡條件下發(fā)動(dòng)機(jī)水下運(yùn)動(dòng)工況，由于發(fā)動(dòng)機(jī)環(huán)境背壓隨著垂直向上運(yùn)動(dòng)而逐漸減小，導(dǎo)致發(fā)動(dòng)機(jī)平穩(wěn)段工作時(shí)的推力逐漸高于固定背壓工況。

圖6　尾部壓力-體積分?jǐn)?shù)云圖

圖7　尾部溫度-馬赫數(shù)云圖

圖8　3個(gè)工況工作全過程推力無量綱化時(shí)變曲線

3.2環(huán)境介質(zhì)影響

通過對(duì)比3個(gè)工況在發(fā)動(dòng)機(jī)工作之初(0～8 ms)的推力特征，如圖9所示?？砂l(fā)現(xiàn)，尾部附著氣泡條件下，發(fā)動(dòng)機(jī)工作之初的推力特征與在空氣介質(zhì)中工作時(shí)結(jié)果非常接近，而在0.5～6 ms階段尾部附著氣泡條件下，發(fā)動(dòng)機(jī)推力小于發(fā)動(dòng)機(jī)在水介質(zhì)中工作的工況?？梢?，外界介質(zhì)對(duì)噴管起動(dòng)過程的推力量值存在影響。且3個(gè)工況均在起始時(shí)刻，由于燃燒室壓力較低，出現(xiàn)了短暫的負(fù)推力現(xiàn)象。

針對(duì)固定背壓的工況，進(jìn)一步將發(fā)動(dòng)機(jī)的推力構(gòu)成進(jìn)行分解，比較該過程發(fā)動(dòng)機(jī)燃燒室內(nèi)壁面和噴管內(nèi)壁面的壓力積分結(jié)果分別如圖10和圖11所示?？煽闯?，環(huán)境介質(zhì)的不同主要影響噴管內(nèi)壁面的壓力，對(duì)燃燒室內(nèi)壁面幾乎沒有影響；相比于環(huán)境介質(zhì)為氣體介質(zhì)的工況，環(huán)境介質(zhì)為重密度的水時(shí)，燃?xì)馍淞饕驘o法瞬時(shí)沖出噴管，其短暫滯留效應(yīng)會(huì)使噴管內(nèi)壁壓力升高，進(jìn)而形成噴管內(nèi)壁面推力的升高，這也是燃?xì)馍淞髦苯由淙胨?，形成推力峰的原因?/p>

圖9　3個(gè)工況工作之初的推力無量綱化時(shí)變曲線

圖10　不同環(huán)境介質(zhì)的燃燒室內(nèi)壁面壓力積分無量綱化時(shí)變曲線

圖11　不同環(huán)境介質(zhì)的噴管內(nèi)壁面壓力積分無量綱化時(shí)變曲線

4　結(jié)論

(1)計(jì)算獲得的發(fā)動(dòng)機(jī)推力和壓力分布結(jié)果與試驗(yàn)值吻合較好，驗(yàn)證了本文所建方法的正確性。

(2)固體火箭發(fā)動(dòng)機(jī)在航行體尾部附體氣泡中工作，使得附體氣泡出現(xiàn)振蕩加強(qiáng)、再次斷裂等復(fù)雜的非定常流動(dòng)特征，噴流流場(chǎng)由附體氣泡逐漸演化為燃?xì)馍淞鹘Y(jié)構(gòu)。

(3)相比于固定環(huán)境背壓水環(huán)境中工作，航行體尾部附體氣泡條件下的運(yùn)動(dòng)發(fā)動(dòng)機(jī)在水環(huán)境工作過程中，發(fā)動(dòng)機(jī)額定工作推力將隨著環(huán)境背壓的降低而增大。

(4)環(huán)境介質(zhì)則主要影響發(fā)動(dòng)機(jī)工作初始過程的推力，當(dāng)環(huán)境介質(zhì)為重密度的水時(shí)，燃?xì)馍淞鞫虝簻粼趪姽軆?nèi)將使噴管內(nèi)壁面壓力升高，進(jìn)而導(dǎo)致總推力短時(shí)增大。

[1]魯傳敬,陳方,樊泓,等.導(dǎo)彈水下點(diǎn)火的流體動(dòng)力學(xué)研究[J].航空學(xué)報(bào),1992,13(4):B124-B130.

[2]張有為,王曉宏,楊舉賢.利用球形氣泡模型研究導(dǎo)彈水下點(diǎn)火瞬間的推力狀況[J].水動(dòng)力學(xué)研究與進(jìn)展,2005,20(5):636-640.

[3]甘曉松,賈有軍,魯傳敬,等.水下燃?xì)馍淞髁鲌?chǎng)數(shù)值研究[J].固體火箭技術(shù),2009,32(1):23-26.

[4]魏海鵬,郭風(fēng)美,權(quán)曉波.水下氣體射流數(shù)值研究[J].導(dǎo)彈與航天運(yùn)載技術(shù),2009(2):37-39,47.

[5]郝宗睿,王樂勤,吳大轉(zhuǎn).水下噴氣推進(jìn)高速射流場(chǎng)數(shù)值研究[J].浙江大學(xué)學(xué)報(bào)(工學(xué)版),2010,44(2):408-412.

[6]胡勇,陳鑫,魯傳敬,等.水下航行體尾噴燃?xì)馀c通氣超空泡相互作用的研究[J].水動(dòng)力學(xué)研究與進(jìn)展,2008,23(4):438-445.

[7]曹嘉怡,魯傳敬,陳鑫,等.導(dǎo)彈水下熱發(fā)射出筒過程流動(dòng)特性[J].固體火箭技術(shù),2011,34(3):281-284,294.

[8]唐嘉寧,李世鵬,王寧飛.水下固體火箭發(fā)動(dòng)機(jī)的負(fù)推力現(xiàn)象研究[J].固體火箭技術(shù),2012,35(3):325-329,343.

[9]王寶壽,許晟,易淑群,等.水下推力矢量特性試驗(yàn)研究[J].船舶力學(xué),2000,4(5):9-15.

[10]湯龍生,劉宇,吳智鋒,等.水下超聲速燃?xì)馍淞鳉馀莸纳L(zhǎng)及壓力傳播特性試驗(yàn)研究[J].推進(jìn)技術(shù),2011,32(3):417-420.

[11]權(quán)曉波,燕國(guó)軍 李巖,等.水下航行體垂直發(fā)射尾空泡生成演化過程三維數(shù)值研究[J].船舶力學(xué),2014,18(7):739-745.

[12]程少華,權(quán)曉波,王占瑩,等.水下航行體垂直發(fā)射尾部空泡形態(tài)與壓力預(yù)示方法研究[J].水動(dòng)力學(xué)研究與進(jìn)展(A輯),2015(3):299-305.

[13]Menter F R,Kuntz M,Lantry R.Ten years of experience with the SST turbulence model[M].Turbulence,Heart and Mass Transfer 4,Begell House Inc.,2003.

(編輯：薛永利)

Numerical simulation research on the jet flow of solid rocket engine in water environment

YAN Guo-jun, QUAN Xiao-bo, WANG Zhan-ying, WEI Hai-peng

(Beijing Institute of Astronautical Systems Engineering,Beijing100076,China)

Based on the VOF multiphase model,SSTk-ωturbulence model and dynamic-mesh method,a simulation model for the solid rocket engine ignition was established to investigate the unsteady characteristics of the jet flow of engine in different environment.The calculated results of thrust curve and measuring pressure distributing were obtained,and the unsteady characteristics of the jet flow in trailing cavity were also analyzed.The gas jet can strengthen pressure oscillation of the cavity and make it fracture once more.The influence of ambient pressure and medium on the engine thrust were studied simultaneously.The results indicate the working medium influences the initial thrust of solid rocket engine.

solid rocket engine;water environment;cavity;numerical simulation

2015-05-20；

2015-08-20。

燕國(guó)軍(1987—)，男，工程師，主要從事流體力學(xué)數(shù)值仿真和理論研究工作。E-mail：yanguojun_calt@126.com

V435

A

1006-2793(2016)03-0312-05

10.7673/j.issn.1006-2793.2016.03.003