水下固體火箭發動機的負推力現象研究①

唐嘉寧，李世鵬，王寧飛

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

水下固體火箭發動機的負推力現象研究①

唐嘉寧，李世鵬，王寧飛

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

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

水下固體火箭發動機;負推力;水下氣體射流;多相流

0 引言

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

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

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

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

1 建模與計算方法

1.1 控制方程

在水下固體火箭發動機的工作過程中，將氣體視為目標流體，基于Volume of Fluid(VOF)多相流模型的質量、動量、能量守恒方程及體積分數輸運方程分別為

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

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

1.2 湍流模型

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

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

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

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

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

1.3 計算域及邊界條件

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

圖1 軸對稱計算域及水下發動機結構Fig.1 The axis-symmetric computational domain together with the geometry of the underwater engine

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

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

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

初始時刻，噴管內為氣相區，噴管外為水環境區，兩個區域均處于靜止狀態。噴管入口處的初始壓強等于實驗中給定的發動機燃燒室初始壓強pinitial，水環境區域的初始壓強等于實驗中發動機所處水深處的環境壓強ρwgh+p1。其中，h為水深;p1為大氣壓強，取值等于0.101 MPa。計算開始時，噴管內外2個區域瞬間聯通，流動建立。

1.4 耦合求解過程

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

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

1.5 推力計算方法

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

式中re、rs分別表示噴管出口和發動機外壁所對應的半徑。

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

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

2 計算結果及分析

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

表2 發動機結構及工作條件Table 2 The geometries and operation conditions of the simulated engine

2.1 計算結果與實驗對比

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

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

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

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

2.2 負推力現象及其產生的原因

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

圖6 H-5型水下發動機在t=0.015 s時的流場結構云圖Fig.6 Flow structures of engine H-5 in the water at t=0.015 s

圖7 H-5型水下發動機在t=0.022 s時的流場結構云圖Fig.7 Flow structures of engine H-5 in the water at t=0.022 s

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

3 結論

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

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

(3)水下發動機點火初期產生負推力現象的原因是發動機點火初期，噴管內被過度壓縮的燃氣沖出噴管后，在噴管尾部形成一個超音速燃氣泡，超音速流動使泡內壓強降低;同時，由于受到流動慣性作用的影響，氣泡持續膨脹使泡內壓強進一步大幅降低，從而導致發動機產生的推力變為吸力，產生負推力現象。

［1］張有為，王曉宏.導彈水下點火推力峰值問題的數值研究［J］.應用力學學報，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，Florida，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］董師顏，張兆良.固體火箭發動機原理［M］.北京:北京理工大學出版社，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—)，女，博士生，研究方向為水下固體火箭發動機。E-mail:tjn1216@163.com

(編輯:崔賢彬)