唐嘉寧，李世鵬，王寧飛

(北京理工大學(xué)宇航學(xué)院，北京 100081)

水下固體火箭發(fā)動機的負(fù)推力現(xiàn)象研究①

唐嘉寧，李世鵬，王寧飛

(北京理工大學(xué)宇航學(xué)院，北京 100081)

針對水下固體火箭發(fā)動機工作環(huán)境壓強高的特點，結(jié)合固體推進劑的燃燒特性，采用UDF方法定義噴管入口邊界條件，建立了固體推進劑燃氣質(zhì)量生成與水下超音速氣體射流的耦合計算模型。將該模型的計算結(jié)果與水下固體火箭發(fā)動機的實驗測量結(jié)果進行對比，驗證了該模型的合理性。研究發(fā)現(xiàn)，水下固體火箭發(fā)動機在點火初期會出現(xiàn)負(fù)推力現(xiàn)象，負(fù)推力產(chǎn)生的原因是發(fā)動機點火初期，噴管內(nèi)被過度壓縮的燃氣沖出噴管后，在噴管尾部形成一個超音速燃氣泡，超音速流動使泡內(nèi)壓強降低;同時受到流動慣性作用的影響，氣泡持續(xù)膨脹使泡內(nèi)壓強進一步大幅降低，發(fā)動機前后端面上的壓差最終導(dǎo)致負(fù)推力現(xiàn)象產(chǎn)生。

水下固體火箭發(fā)動機;負(fù)推力;水下氣體射流;多相流

0 引言

固體火箭發(fā)動機以其結(jié)構(gòu)簡單、操作使用方便、快速反應(yīng)能力強等突出優(yōu)點，成為各類主動攻擊型深海武器的首選動力。發(fā)動機在水下工作時，由于外部水介質(zhì)的密度大、環(huán)境壓強高等特點，使其工作過程變得極為復(fù)雜［1］。例如，發(fā)動機點火后，水環(huán)境的慣性作用嚴(yán)重制約噴管內(nèi)氣體的噴出，可能使發(fā)動機工作故障或失敗;高速燃氣射流與周圍水環(huán)境之間的相互作用，導(dǎo)致推力難以預(yù)估;深水情況下推力和阻力存在耦合關(guān)系，使得深水工作推進系統(tǒng)的實際推力與理論計算值相差較大;深水發(fā)動機推力損失異常，難以用簡單的背壓變化機理解釋等。因此，深入了解水下發(fā)動機推力特性，具有十分重要的理論意義和工程應(yīng)用價值。

在早期研究中，Rogers［2］提出了球形氣泡模型，并運用該模型對水下發(fā)動機點火瞬間產(chǎn)生的最大推力進行了預(yù)測。Labotz［3］在球形氣泡模型的基礎(chǔ)上，對發(fā)動機水下點火后達到平衡壓強的最短時間進行了估計。Tang［4］等將球形氣泡模型的計算結(jié)果與CFD計算結(jié)果進行對比，發(fā)現(xiàn)球形氣泡模型結(jié)構(gòu)簡單，計算速度快，并能有效反映發(fā)動機點火瞬間的推力峰值特性。但由于該模型假設(shè)氣體在整個射流過程中以單一氣泡結(jié)構(gòu)存在，并沿徑向等速無限膨脹，因此不能用于說明較長時間內(nèi)的壓強及推力特性。Lindau［5］等對考慮空泡作用下的高速超空化水下推進過程進行了數(shù)值研究，提出了包括燃氣、水和水蒸氣的三相流計算模型，但沒有對水下推進過程中的流動特性和發(fā)動機工作性能做進一步分析。

Tang等［6］建立了水下超音速氣體射流數(shù)學(xué)模型，并將數(shù)值計算結(jié)果與實驗結(jié)果進行對比，驗證了該模型的合理性和準(zhǔn)確性，為水下固體火箭發(fā)動機工作過程的研究提供了基礎(chǔ)。該水下超音速氣體射流數(shù)學(xué)模型中采用固定壓力入口邊界條件，即噴管入口處的總壓從工作開始到結(jié)束都保持恒定，不會受到下游流場的影響。然而，在實際的水下發(fā)動機工作過程中，一方面存在燃燒室建壓延遲過程;另一方面，外部水環(huán)境高壓可能對燃燒室內(nèi)工作情況產(chǎn)生影響。因此，噴管入口處的總壓并不能保持恒定。例如，外部水介質(zhì)在點火瞬間產(chǎn)生的慣性作用會對噴管內(nèi)流動產(chǎn)生阻滯，使燃燒室內(nèi)部壓強升高;高速燃氣與水體發(fā)生強烈撞擊后產(chǎn)生的強激波可能反傳入噴管內(nèi)，使上游燃燒室壓強發(fā)生振蕩，從而導(dǎo)致固體推進劑燃速和燃氣生成率發(fā)生變化。當(dāng)發(fā)動機所處水深較大、外部水環(huán)境壓強較高時，此影響更為嚴(yán)重。因此，Tang等［6］提出的水下超音速氣體射流數(shù)學(xué)模型，適用于噴管喉部達到音速流或不考慮發(fā)動機點火初期工作特性的情況。

為準(zhǔn)確模擬水下固體火箭發(fā)動機工作過程，并深入分析點火初期的工作特性，本文將在Tang等［6］提出的水下超音速氣體射流數(shù)學(xué)模型的基礎(chǔ)上，結(jié)合固體推進劑的燃燒特性，對燃燒室內(nèi)固體推進劑燃氣質(zhì)量生成過程和水下超音速氣體射流過程進行耦合求解。

1 建模與計算方法

1.1 控制方程

在水下固體火箭發(fā)動機的工作過程中，將氣體視為目標(biāo)流體，基于Volume of Fluid(VOF)多相流模型的質(zhì)量、動量、能量守恒方程及體積分?jǐn)?shù)輸運方程分別為

其中，φm為混合相的物理性質(zhì)，由 φm=φgαg+φw(1-αg)計算得到，φ可用來表示密度、粘性系數(shù)等物理量;下標(biāo)w表示水，g表示氣體。氣液兩相的體積分?jǐn)?shù)關(guān)系為αw+αg=1。模型中，將氣體假設(shè)為理想氣體，水為不可壓縮相。由于水下氣體射流過程中的Froude數(shù)約為0(103)，因此在動量守恒方程(2)中，可忽略體積力Fi的作用。在能量方程(3)中，Sh為源項;keff為有效熱傳導(dǎo)率(k+kt)，其中kt是湍流熱傳導(dǎo)率，根據(jù)所采用的湍流模型定義;能量Ek可由下式?jīng)Q定:

在本文的研究中，假設(shè)火箭發(fā)動機固定不動，整個射流過程屬于高壓范疇，不考慮由于超低壓作用產(chǎn)生的空泡現(xiàn)象，且不考慮氣-液相之間的相變作用。因此，體積分?jǐn)?shù)輸運方程(4)右端項為0。

1.2 湍流模型

選用標(biāo)準(zhǔn)k-ε湍流模型，對雷諾平均N-S方程進行封閉。該湍流模型適合完全湍流流動，是針對高雷諾數(shù)的湍流計算模型。湍流動能k及湍流能量耗散率ε由以式(6)及式(7)求得:

式中Gk和Gb分別表示由時均速度梯度和浮力產(chǎn)生的湍流動能;YM體現(xiàn)了流體可壓縮性，定義為YM=，其中Mt=(k/ γRT)1/2為湍流馬赫數(shù);Sk和Sε是用戶定義的源項。

由于不考慮體積力作用，且流場中無源項存在，所以Gb、Sk和Sε均為0。湍流粘性系數(shù)定義為

其中，C1ε、C2ε、Cμ、σk和 σε均為常數(shù)，取值如表 1 所示。

表1 湍流模型中的常值參數(shù)Table 1 Constant parameters of the turbulence model

1.3 計算域及邊界條件

采用二維軸對稱模型，對水下固體火箭發(fā)動機工作過程進行計算。計算域包括噴管內(nèi)流動區(qū)域和噴管外射流區(qū)域。計算域如圖1所示，總長度為噴管出口直徑的600倍，寬度為噴管出口直徑的35倍，以克服邊界條件對主射流區(qū)流動的影響。計算中，將整個計算區(qū)域分解為多個區(qū)域，并對各個區(qū)域分別生成計算網(wǎng)格。同時，根據(jù)流動物理特征分析，在壁面、噴管喉部及射流核心區(qū)等物理量變化梯度較大區(qū)域進行網(wǎng)格加密。

圖1 軸對稱計算域及水下發(fā)動機結(jié)構(gòu)Fig.1 The axis-symmetric computational domain together with the geometry of the underwater engine

由于在水下發(fā)動機的工作過程中，固體推進劑燃速會影響發(fā)動機燃燒室壓強的變化過程;同時，燃燒室壓強又將反過來影響推進劑燃速，引起噴管質(zhì)量流量和下游氣-液兩相流場的變化。因此，為了結(jié)合固體推進劑的燃燒特性，研究中借助用戶自定義函數(shù)(UDF)，使用C語言編寫擴展FLUENT的程序代碼，然后動態(tài)加載到FLUENT環(huán)境中，對噴管入口邊界條件進行定義。噴管入口處采用質(zhì)量入口邊界條件，入口質(zhì)量流量為

式中表示噴管入口處氣體的質(zhì)量流量為固體推進劑燃氣質(zhì)量生成率，由推進劑燃燒特性和燃燒室壓強共同決定;為固體推進劑燃速;表示入口處水的質(zhì)量流量，因為入口處介質(zhì)為燃氣，所以=0;燃燒室壓強pc會受到發(fā)動機外部氣-液兩相流場的影響，并不保持恒定。

在自由出口邊界上，總壓等于環(huán)境壓強，由噴管所處的水深決定;總溫等于水環(huán)境的溫度，假設(shè)為300 K。壁面采用固壁面邊界條件，無滑移且絕熱;流場中心采用軸對稱邊界條件。

初始時刻，噴管內(nèi)為氣相區(qū)，噴管外為水環(huán)境區(qū)，兩個區(qū)域均處于靜止?fàn)顟B(tài)。噴管入口處的初始壓強等于實驗中給定的發(fā)動機燃燒室初始壓強pinitial，水環(huán)境區(qū)域的初始壓強等于實驗中發(fā)動機所處水深處的環(huán)境壓強ρwgh+p1。其中，h為水深;p1為大氣壓強，取值等于0.101 MPa。計算開始時，噴管內(nèi)外2個區(qū)域瞬間聯(lián)通，流動建立。

1.4 耦合求解過程

耦合求解過程如圖2所示。首先，用給定的發(fā)動機燃燒室初始壓強進行初始化，開始第一步迭代計算，燃燒室入口邊界處的燃氣質(zhì)量流量根據(jù)初始壓強下的推進劑燃速計算得到;然后，對分離求解后得到的計算結(jié)果進行判斷，若計算不收斂，則繼續(xù)迭代，若計算收斂，則本次循環(huán)終止;同時，運用UDF讀取本次循環(huán)結(jié)束后得到的燃燒室入口邊界處的壓強值，并根據(jù)新讀取的壓強值計算得到新的燃氣生成率，將其作為下一步迭代計算的入口質(zhì)量流量值，開始下一次迭代，計算時間步長Δt=10-6s。本文研究中，采用有限體積法對控制方程進行離散，并運用SIMPLE算法對流場進行求解?？紤]到氣-液兩相流數(shù)值計算量大的特點，選用一階迎風(fēng)格式進行插值。

圖2 燃氣質(zhì)量生成與水下超音速氣體射流的耦合求解流程Fig.2 The coupling simulation method between the process of gas generation and the process of gaseous jets injected into water

1.5 推力計算方法

火箭發(fā)動機在水下工作時，由于高速燃氣與水之間的相互作用以及射流過程中的膨脹、脹鼓、頸縮/斷裂和回?fù)衄F(xiàn)象的出現(xiàn)，整個工作過程為非穩(wěn)態(tài)的湍流過程［6］。另一方面，發(fā)動機外壁所受到的外部介質(zhì)的靜壓作用非均勻分布，且發(fā)動機尾部所受到的背壓不等于環(huán)境壓強，pa≠pB。因此，水下火箭發(fā)動機的推力不可按空氣中的簡化推力公式［7］來計算。根據(jù)圖3所示，對發(fā)動機內(nèi)外表面所受到的壓力進行積分，可得水下發(fā)動機推力計算公式為

式中re、rs分別表示噴管出口和發(fā)動機外壁所對應(yīng)的半徑。

由式(9)可知，水下火箭發(fā)動機推力由3部分組成:第一部分為動量推力，由高速燃氣產(chǎn)生的動量變化率提供;第二部分為壓力推力，由噴管出口截面上的壓強pe與外界環(huán)境壓強pa之間的壓差產(chǎn)生;第三部分為壓力推力，由作用于發(fā)動機前后端面上環(huán)境壓強pa與背壓pB之間壓力差產(chǎn)生。

圖3 發(fā)動機推力計算方法示意圖Fig.3 Computation schematic of the rocket engine thrust and pressure distribution

2 計算結(jié)果及分析

在本文的計算中，選取的裝藥為端燃藥柱，燃面Ab在燃燒過程中保持不變。推進劑燃速為r==cm/s(pc單位為 kg/cm2)，推進劑密度ρP=1.5×103kg/m3。不同的燃燒室平衡壓強可通過調(diào)整裝藥燃燒面積來實現(xiàn)。發(fā)動機尺寸和詳細工作參數(shù)如表2所示。其中，dt、de、ds分別表示噴管喉部、出口及發(fā)動機外壁所對應(yīng)的直徑;Ae/At為噴管擴張比。

表2 發(fā)動機結(jié)構(gòu)及工作條件Table 2 The geometries and operation conditions of the simulated engine

2.1 計算結(jié)果與實驗對比

運用本文建立的燃氣質(zhì)量生成與水下超音速氣體射流的耦合計算模型，對表2中的H-5型發(fā)動機工作過程進行模擬，得到推力-時間曲線如圖4(a)所示。實驗測量得到的推力-時間曲線如圖4(b)所示。圖4(a)和(b)中，縱坐標(biāo)為無量綱化后的推力F/Fref，無量綱參數(shù)Fref為發(fā)動機燃燒室壓強達到平衡后的平均推力。比較兩圖可知，水下發(fā)動機在點火初期會出現(xiàn)推力峰值，峰值后均出現(xiàn)了負(fù)推力現(xiàn)象，之后推力緩慢爬升達到平衡，在整個過程中伴隨有強烈的推力振蕩。推力爬升是建壓過程中燃燒室壓強不斷升高引起的，推力振蕩是水下超音速氣體射流過程中膨脹、脹鼓、頸縮/斷裂和回?fù)衄F(xiàn)象［6］循環(huán)出現(xiàn)的結(jié)果。通過耦合計算方法得到的發(fā)動機推力在開始瞬間為負(fù)值，如圖4(a)所示，而通過實驗方法測量得到的推力在開始瞬間為0，如圖4(b)所示。這是由于實驗中的發(fā)動機噴管出口處裝有尾蓋，只有當(dāng)尾蓋左側(cè)(噴管內(nèi))的壓強比尾蓋右側(cè)(噴管外)的壓強高出一定值時，尾蓋才會被吹開，以防止外部流體倒流，因此實驗中的發(fā)動機推力在開始瞬間為0。耦合計算過程中，沒有考慮噴管尾蓋的影響，即計算開始時，噴管內(nèi)外的兩個區(qū)域瞬間聯(lián)通。因此，在外部環(huán)境壓強較高的情況下，噴管出口處的流體會在初始瞬間出現(xiàn)倒流現(xiàn)象，使推力在初始瞬間為負(fù)值。由圖4(a)和(b)可知，通過耦合計算方法和實驗方法得到的推力-時間曲線變化趨勢吻合較好，說明了本文提出的耦合計算模型的合理性。

圖4 水下固體火箭發(fā)動機推力隨時間的變化曲線Fig.4 Time vs thrust profile of the underwater solid rocket engine

圖5(a)為通過耦合計算方法得到的H-5型發(fā)動機噴管內(nèi)的壓強-時間變化曲線，測壓點位于噴管擴張段內(nèi)距離噴管出口0.1 m處。圖5(b)為通過實驗方法測量得到的壓強-時間曲線。兩圖中，縱坐標(biāo)為無量綱化后的壓強p/pref，參考壓強pref為發(fā)動機燃燒室壓強達到平衡后測壓點處的平均壓強值。比較兩圖可知，在發(fā)動機點火初期，噴管擴張段內(nèi)會出現(xiàn)很高的壓強峰值，壓強從峰值回落后將出現(xiàn)短期振蕩，然后緩慢爬升至平衡值。壓強峰出現(xiàn)的原因為發(fā)動機點火初期，大量燃氣沖出噴管與水發(fā)生撞擊，在水環(huán)境的靜壓和慣性作用的影響下，噴管擴張段內(nèi)產(chǎn)生了反傳壓力波，當(dāng)壓力波到達測壓點時，便會引起壓力峰值的出現(xiàn)。發(fā)動機水下點火過程中出現(xiàn)在噴管內(nèi)的壓強峰是水下工作的特別現(xiàn)象，相同發(fā)動機在地面實驗中，并沒有觀察到該現(xiàn)象。實驗測量得到的壓強在爬升之后，還會出現(xiàn)輕微振蕩(如圖5(b)所示)，主要原因是當(dāng)噴管擴張段內(nèi)壓強較低、環(huán)境壓強較高，或測壓點離噴管口較近時，測壓點工作情況容易受到外部氣-液兩相流動的影響，使噴管內(nèi)壓力出現(xiàn)波動。由圖5(a)和(b)可知，通過耦合計算方法和實驗方法得到噴管擴張段內(nèi)的壓強-時間曲線變化趨勢吻合較好，說明了本文提出的耦合計算模型合理性。

圖5 水下固體火箭發(fā)動機噴管內(nèi)壓強隨時間的變化曲線Fig.5 Time vs pressure profile tested in the nozzle of the underwater solid rocket engine

2.2 負(fù)推力現(xiàn)象及其產(chǎn)生的原因

基于本文的數(shù)值模擬結(jié)果，對發(fā)動機點火初期的負(fù)推力現(xiàn)象進行深入研究。由圖4(a)可知，發(fā)動機點火初期首先出現(xiàn)較高的推力峰值，對應(yīng)工作時間為0.015 s。峰值之后，推力迅速降低到負(fù)值，出現(xiàn)負(fù)推力現(xiàn)象，對應(yīng)時間為0.022 s。為了分析負(fù)推力現(xiàn)象出現(xiàn)的原因，圖6和圖7中分別給出了工作時間為0.015 s(推力峰值出現(xiàn)時刻)和0.022 s(負(fù)推力出現(xiàn)時刻)時的流場液體體積分?jǐn)?shù)云圖和壓強分布云圖。首先，從圖6(b)中可知，當(dāng)工作時間t=0.015 s時，燃燒室產(chǎn)生的高速燃氣與慣性作用很強的水發(fā)生強烈撞擊，在噴管尾部產(chǎn)生一個高壓區(qū)。此時，噴管出口壓強和背壓都急劇升高，從而導(dǎo)致推力峰值的出現(xiàn)，如圖4(a)所示。由于此時燃燒室壓強尚未達到平衡壓強值，噴管內(nèi)部壓強較低，外部環(huán)境高壓會把周圍的水介質(zhì)反推入噴管內(nèi)，使噴管內(nèi)流動受阻并達不到超音速流動，如圖6(a)所示。

圖6 H-5型水下發(fā)動機在t=0.015 s時的流場結(jié)構(gòu)云圖Fig.6 Flow structures of engine H-5 in the water at t=0.015 s

圖7 H-5型水下發(fā)動機在t=0.022 s時的流場結(jié)構(gòu)云圖Fig.7 Flow structures of engine H-5 in the water at t=0.022 s

隨著燃燒室壓強的升高，以及外部水流場流動的建立，噴管內(nèi)流動受到的阻礙逐漸減弱，燃氣在被過度壓縮之后，進行快速膨脹，噴管中的超音速氣流終于沖出噴管口，并在噴管尾部形成一個超音速燃氣泡，如圖7(a)所示。燃氣泡內(nèi)的超音速流動，使泡內(nèi)壓強降得很低，在流動慣性作用的影響下，氣泡將持續(xù)膨脹，使泡內(nèi)壓強進一步降低，并低于環(huán)境壓強，如圖7(b)所示。

3 結(jié)論

(1)建立了固體推進劑燃氣質(zhì)量生成與水下超音速氣體射流流動的耦合計算模型，通過與水下固體火箭發(fā)動機的實驗數(shù)據(jù)對比，驗證了模型的合理性。模型可反映出水下發(fā)動機外部氣-液兩相流場對燃燒室壓強和推進劑燃速的影響，為深水發(fā)動機的研究提供了基礎(chǔ)。

(2)水下發(fā)動機在點火初期會出現(xiàn)推力爬升和負(fù)推力現(xiàn)象，推力爬升是建壓過程中燃燒室壓強不斷升高引起的。

(3)水下發(fā)動機點火初期產(chǎn)生負(fù)推力現(xiàn)象的原因是發(fā)動機點火初期，噴管內(nèi)被過度壓縮的燃氣沖出噴管后，在噴管尾部形成一個超音速燃氣泡，超音速流動使泡內(nèi)壓強降低;同時，由于受到流動慣性作用的影響，氣泡持續(xù)膨脹使泡內(nèi)壓強進一步大幅降低，從而導(dǎo)致發(fā)動機產(chǎn)生的推力變?yōu)槲?，產(chǎn)生負(fù)推力現(xiàn)象。

［1］張有為，王曉宏.導(dǎo)彈水下點火推力峰值問題的數(shù)值研究［J］.應(yīng)用力學(xué)學(xué)報，2007，24(1):298-301.

［2］Rogers K W.A theoretical and experimental investigation of the transient phase of underwater rocket motor firing［R］.U-niversity of Southern California Engineering Center Report，1962.

［3］Labotz R J.Hydrodynamic consideration and limitations in submerged rocket firings［J］.Journal of Spacecraft and Rocket，1965，2(3):320-324.

［4］Tang J N，Li S P，Wang N F，et al.Flow structures of gaseous jet injected into liquid for underwater propulsion［C］//46th AIAA/ASME/SAE/ASEE Joint Propulsion Conference and Exhibit，Nashville，TN.2010，AIAA 2010-6911.

［5］Lindau J W，Venkateswaran S，Kunz R F，et al.Multiphase computations for underwater propulsive flows［C］//16th AIAA Computational Fluid Dynamics Conference，Orlando，F(xiàn)lorida，AIAA 2003-4105.

［6］Tang Jia-ning，Wang Ning-fei，Wei Shyy.Flow structures of gaseous jets injected into water for underwater propulsion［J］.Acta Mechanica Sinica，2011，27(4):461-472.

［7］董師顏，張兆良.固體火箭發(fā)動機原理［M］.北京:北京理工大學(xué)出版社，1996.

Study on the negative thrust of the underwater solid rocket engines

TANG Jia-ning，LI Shi-peng，WANG Ning-fei
(School of Aerospace Engineering，Beijing Institute of Technology，Beijing 100081，China)

To simulate the working process of the underwater solid rocket engines，a coupling model between the process of gas generation and the process of gaseous jets injected into water was proposed.The burning rate of the solid propellant was considered in this coupling model and the User Defined Functions(UDF)were used to define the boundary condition of the nozzle inlet.The numerical method is validated by the experimental data.The results show that the negative thrust appears when the engines work in the water.The reason is as follows:at the initial stage of the propulsion process，the high compressed gas jets into water to form a supersonic“gas bag”behind the nozzle exit，causing the pressure in this area to drop suddenly.Furthermore，the inertia effect of the water causes the pressure in the“gas bag”to keep dropping，creating the negative thrust due to the pressure difference between the rear and front parts of the engine.

underwater solid rocket engine;thrust;submerged gaseous jets;multiphase flow

V435

A

1006-2793(2012)03-0325-05

2011-11-09;

2011-12-01。

唐嘉寧(1984—)，女，博士生，研究方向為水下固體火箭發(fā)動機。E-mail:tjn1216@163.com

(編輯:崔賢彬)