固體火箭發動機水下超音速射流數值研究

王利利， 劉影， 李達欽， 吳欽， 王國玉

(北京理工大學 機械與車輛學院, 北京 100081)

0 引言

潛射導彈具有隱蔽性好、系統簡便靈活、體積小、質量輕等優點，作為可靠的二次核打擊力量，越來越受到各國的重視[1]。固體火箭發動機由于結構簡單、可靠性高、快速反應能力強等優點，逐漸成為潛射導彈的首選動力源[2]。然而，由于發動機水下工作環境介質密度遠大于空氣，且處于高溫、高壓、超音速工作狀態，其射流過程十分復雜[3]：一方面，大量的高溫燃氣噴入水中，會在導彈尾部形成瞬態燃氣泡，高溫燃氣泡與水之間存在氣- 水摻混、換熱和相變等復雜物理現象并伴有激烈的流體動力干擾[4]；另一方面，燃氣泡中存在的復雜激波結構[5-6]以及激波結構的運動[7]均會使流場壓力劇烈變化，從而導致推力性能難以預測。近年來，為了提高發動機動力系統的推力性能，其采用的噴管擴張比越來越大，導致其穩定性大大降低。因此，深入研究大擴張比水下發動機的工作特性及其相關兩相流動機理，對水下發射技術發展具有重要意義。

國內外學者對水下超音速氣體射流過程開展了大量研究。王柏懿等[8]、Shi等[9]、施紅輝等[10-11]通過試驗對擴張比為1.2、馬赫數為1～2的水下氣體射流過程進行了觀察研究，并利用探針排測量了尾流場及噴管壁面上的壓力脈動情況，試驗結果表明，氣體中出現的大幅值壓力脈沖與在噴管出口兩倍直徑處的射流形貌存在關聯性，噴口端面處的回擊現象與流場氣相區中的壓力脈動相關，射流的脹鼓與回擊均可引起壓力脈動,回擊頻率為脹鼓頻率的2～3倍。Fronzeo等[12]針對擴張比為1.7的噴管在常溫下不同密度液相環境下的射流過程進行了數值研究，研究發現液體密度是制約射流貫穿距的主要因素，氣- 液界面的湍流水平隨著環境液體密度的增大而急劇增加，導致混合層變化劇烈、激波胞室結構較早破裂。在此基礎上， Li等[13]結合試驗與數值計算方法，對含有蒸汽的氣體射流提出了一種預測射流貫穿距動量平衡模型。李婷婷等[14]用高速攝像機觀測了不同流量下環形噴管的射流夾斷過程，結果表明，頸縮導致氣體內部壓力突增，從而引發射流夾斷現象，夾斷的頻率隨著氣體流量的增加而降低。Wilson等[15]對1 800 K氣體、馬赫數分別為1和2的水下超音速射流進行了試驗研究，分別考察了圓形、矩形、橢圓形和三角形噴嘴的射流情況，研究表明在使用橢圓形及矩形噴管時，氣體射流的傳播速度及氣- 水之間的混合要比使用圓形噴管時快，而與周圍低溫水環境之間強烈的混合作用會使射流總溫降得較低。由于試驗中所用的射流入口壓強較低(111～237 Pa)，且沒有對速度、壓力等參數進行測量，無法充分表征發動機水下工作特性。

水下超音速氣體射流過程的推力特性及影響因素一直是研究的重點與難點。王寶壽等[16]通過試驗測量了噴管擴張比為3、滯止壓強為9 MPa、氣體初始溫度為2 700 K、水深為10～40 m條件下固體火箭發動機水下工作時的推力和側向力,獲得了推力、側向力、推力矢量角隨發動機噴口堵塞面積比、水深的變化規律。張有為[3]利用球形氣泡模型對發動機水下推力特征的影響因素進行了研究，結果表明，推力峰值隨著燃氣泡初始半徑的增大而迅速降低，隨著燃面面積的減小而降低。

已有研究表明，水下超音速射流往往發生周期性的變化，表現為形態上的周期性變化，壓力場、速度場等出現循環振蕩，該現象對推力將產生較大影響。然而現有研究大多關注流場整體特征的變化，且對大擴張比噴管在超高溫、超高壓工作條件下的超音速水下射流研究較少。本文采用計算流體力學(CFD)方法，針對高溫高壓下兩種擴張比噴管的射流開展模擬研究，揭示射流流場變化與推力變化機理,以獲得擴張比對流場結構與推力的影響規律,所得研究成果對研究水下氣體射流與固體火箭發動機水下推進具有重要意義。

1 數值計算方法

1.1 流體體積方法及其控制方程

將氣體假設為理想氣體，水為不可壓縮相，采用流體體積方法對氣- 液運動界面進行追蹤。連續方程、動量方程、能量方程、體積輸運方程和狀態方程分別如下：

(1)

(2)

式中：ρm為汽- 液混合密度；μm為汽- 液混合黏性系數；ρm、μm可以統一用物質屬性φm來表示，φm=φgαg+φwαw，φg、φw分別表示氣相密度或黏性系數、液相密度或黏性系數，αg、αw分別表示氣相體積分數、液相體積分數，ui、uj為速度分量；p為壓強；Fi為體積力，忽略重浮力作用下Fi=0.

(3)

式中：氣- 液兩相體積分數關系為αg+αw=1；能量En和溫度Tt為各項的質量平均變量；γeff為有效熱傳導率；Sh為源項。

(4)

(5)

式中：R為理想氣體常數。

1.2 基于濾波器的修正湍流模型

本文采用基于濾波器的修正湍流模型[17]，在標準k-ε湍流模型中加入濾波函數，對尺度小于濾波器尺寸的湍流采用標準k-ε模型(見(6)式)；對尺度大于濾波器尺寸的湍流結構則采用直接計算方法求解(見(7)式)。

(6)

(7)

式中：ρ為流體密度；k為湍動能；ε為動能耗散率；Pk為湍動能生成項；模型常數σk=1.0，Cε1=1.44，Cε2=1.92，σε=1.30；μt為湍流黏性系數，

(8)

模型常數Cμ=0.09，F為濾波函數，由濾波器尺寸λ和湍流長度比值的大小決定，

(9)

式中：模型常數C3=1.0.

由(9)式可知，當湍流尺度較大時，湍流黏性系數表達為

(10)

1.3 計算模型及邊界條件

圖1 邊界條件及幾何示意Fig.1 Boundary conditions and geometry

圖2 入口無量綱壓力與入口溫度函數曲線Fig.2 Curves of inlet pressure and temperature

根據流動物理特征分析，在壁面、噴管喉部及射流核心區等物理量變化梯度較大的區域進行網格加密，其局部放大圖如圖3所示，計算時間步長Δt=10-6s.

圖3 網格劃分結果及其放大圖Fig.3 Mesh generation and its enlarged view

擴張比η為出口截面面積與喉部截面面積之比，本文以擴張比η=14與η=3.4的噴管為例分析流場和推力變化規律，二者受縮段、喉部段尺寸一致，控制出口直徑de以獲得不同的擴張比。工作參數如表1所示，其中Ma為設計工況出口馬赫數。

表1 噴管工作參數

1.4 數值結果驗證

本文根據文獻[9]試驗結果進行數值方法驗證，噴管結構尺寸與工作條件均與文獻[9]一致。

表2為射流發展過程形貌對比。從表2中可以看出，數值模擬結果與試驗結果呈現出相同的周期性變化特征，數值計算結果對頸縮位置、脹鼓泡和回擊流的捕捉均與文獻[9]較為一致。

表2 射流發展形貌的數值模擬結果與試驗結果[9]對比

Tab.2 Comparison of numerically simulated and experimental results of jet pattern

圖4 尾流場靜壓均方根數值模擬結果與試驗結果[9] 對比Fig.4 Comparison of numerically simulated and experimental results[9] about flow field static pressure RMS

1.5 推力計算方法

固體火箭發動機水下所受推力示意圖如圖5所示，其中As為發動機外殼橫截面積，pb為發動機尾部背壓，ve為出口排氣速度，me為出口排氣質量。由圖5可見，固體火箭發動機在水下工作時，由于整個工作過程為非穩態湍流過程，發動機外壁所受到的外部介質靜壓作用非均勻分布，且發動機尾部所受到的背壓不等于環境壓強，即pb≠pa. 因此，固體水下火箭發動機的推力不可按空氣中的簡化推力公式[18]來計算。

圖5 推力示意圖Fig.5 Schematic diagram of thrust

考慮背壓pb后的發動機推力計算公式為

(11)

由(9)式可知，推力由3部分組成：

2) 壓差推力Ae(pe-pa)，由噴管排氣面上的壓強pe與外界壓強pa之間的不平衡產生，取決于排氣截面的燃氣參數pe、Ae和發動機所處的環境壓強pa.

3) 壓差推力(As-Ae)(pb-pa)，由作用于發動機前后端面上的環境壓強pa與背壓pb之間的壓力差產生。

由于不同發動機的外壁尺寸As差異較大，As取不同值時對推力公式中的最后一項影響較大，故省略(11)式最后一項，化為積分形式：

(12)

式中：re為噴管出口半徑；ρe為出口截面上流體的密度；rs為出口截面半徑。

2 結果與討論

2.1 不同擴張比噴管射流瞬態演化過程

表3給出了不同擴張比噴管在初始膨脹階段和充分發展階段的典型水氣云圖與壓力云圖，p*=(p-pa)/0.5Aeρrefu2，以入口中心處為坐標零點，橫坐標x*=x/de, 縱坐標y*=y/de.

表3 射流發展過程云圖

Tab.3 Water-vapor contour of jet evolution process

由表3可見：

初始膨脹時期，在氣囊推開液相作用下，氣- 液界面壓力略有升高(如虛線框所示)；由于氣囊處于膨脹階段，一部分流線輻射向外，另一部分流線在靠近管口處形成回流。

在充分發展時期，壓力波動集中在噴管出口外(如虛線框所示)，隨著射流向下游發展，壓力逐漸與環境一致；射流下游形成了大尺寸旋渦，這是因為膨脹氣體隨主射流通道向下游噴射，造成水與氣體裹挾與摻混。對比兩種噴管射流形態可知，η=14.0的噴管射流軸向貫穿距遠大于3.4的噴管，這是因為擴張段的大擴張比增加了氣體膨脹慣性。

為了進一步研究射流發展的瞬態演變規律，根據充分發展時期射流特征周期性變化特點，將其分為頸縮、脹鼓和回擊3個階段，表4、表5給出了兩種噴管在各階段的壓力云圖與水氣輪廓線、流線圖與速度云圖、渦流結構判別標準Q云圖，Q以其最大值進行無量綱化處理。Q定義為Q=1/2(|Ω|2-|S|2)，其中Ω為渦度張量，S為應變張量變化率。在頸縮階段，氣流在垂直射流的方向上被壓縮，射流通道變窄，由于主射流通道阻塞，噴管口形成了較大范圍的高壓區，導致噴管壁壓力突增。噴管口流線呈輻射狀，未出現氣- 液倒灌現象，表明高壓高溫條件有效遏制了流動分離現象。在脹鼓階段，氣體阻塞在頸縮點上游，發生膨脹，且由于過膨脹氣體推開水的能力較弱，脹鼓位置緊貼噴管口外壁；隨著脹鼓泡內壓力能的釋放，脹鼓泡膨脹，出口速度增加，噴管內部高壓區消失；在脹鼓泡內部，圍繞主射流通道形成較高速度的回流，出現了Q的最高值，如表4、表5中脹鼓階段Q云圖所示。在回擊階段，脹鼓泡內累積能量與液相能量平衡最終沖破頸縮點，射流通道打開，形成了穩定的高速射流，此時出口速度達到設計馬赫數。由表4中回擊階段Q云圖可以看到，形成了大范圍的回流，回流氣體來自脹鼓泡內被主射流帶動形成負速度的氣體。同時，噴管口附近形成低壓區，如表4、表5中壓力云圖所示。由于回流對主射流通道有剪切作用、高壓水環境對射流通道的擠壓作用，射流通道在噴管口附近再次發生頸縮，上述過程將周期性地重復。

表4η=14.0時射流發展過程云圖各階段水氣分布圖、壓力云圖、Q云圖(局部放大)

Tab.4 Partially enlarged contours of water vapor distribution, pressure andQforη=14.0

表5 η=3.4時射流發展過程云圖各階段水氣分布圖、壓力云圖、Q云圖(局部放大)

對比擴張比η=14.0與η=3.4的噴管射流形貌可知：η=3.4的噴管在水氣輪廓圖上未發生較大變化，頸縮階段沒有出現明顯的主射流通道變窄現象，這是因為擴張比小，氣體離開噴管后延徑向擴張的慣性較小，故仍在軸向保持較大的動能與壓力能；η=3.4的噴管出口維持穩定的超音速射流，除了在頸縮點出現小范圍高壓區，出口一直維持高速低壓狀態，而η=14.0噴管在頸縮階段出口速度為亞音速；η=3.4噴管出口處沒有觀察到紊亂的流線，且該處氣體徑向膨脹的趨勢小，故流線未呈輻射狀。脹鼓階段在水氣輪廓圖上并未出現明顯脹鼓泡，出口速度增加，由高壓轉為低壓，回流產生在液相區。在回擊階段出口速度增加，出現大范圍低壓區；從流線圖上可以看到，回流并未直接回擊到管壁，表明該噴管回擊帶來的壓力振蕩與腐蝕較小。高速氣體射流氣- 液邊界的相對運動較穩定，氣- 液摻混較為緩和。

2.2 不同擴張比噴管的推力特性

為了進一步研究不同擴張比噴管推力特性，給出了噴管出口中心處監測點A(見圖6)的瞬態壓力變化，如圖7所示，其中橫坐標為時間無量綱化結果t*=t/T. 由圖7可知：η=14.0時，在初始膨脹時期呈現高頻隨機性震蕩，在充分發展時期，由于激波推入管內，發生明顯的脈沖振蕩，平均周期為0.010 12 s；對比發現，η=3.4的壓力曲線基本呈現一條直線，未出現脈沖，這是因為擴張比較小，出口壓力較大，形成了穩定射流，激波一直未進入噴管內。

圖6 監測點A位置Fig.6 Location of Point A

圖7 監測點A壓力變化曲線Fig.7 Pressure at Point A

圖8 單個周期內動量推力與壓差推力曲線Fig.8 Momentum thrust and differential pressure thrust curves in a single cycle

圖9 噴管口附近壓力云圖Fig.9 Pressure contours near nozzle

2.3 不同擴張比噴管背壓振蕩特性

作為推力的影響參數，背壓振蕩與推力振蕩有相同的變化規律。為了研究背壓振蕩規律，在出口側壁上布置4個監測點B1、B2、B3、B4如圖10所示，每個監測點相距10 mm. 圖11和圖12分別給出了B1點壓力振蕩曲線及相應快速傅里葉變換頻譜圖；表6給出了各點的主頻。從表6中可知，各點主頻相近，η=14.0的噴管均在100 Hz左右，η=3.4的噴管均在200 Hz左右。由于高頻噪聲在水下衰減更快，η=3.4的噴管在水下推進時隱蔽性更好。

圖10 監測點位置及放大圖Fig.10 Locations of monitoring points

圖11 監測點B1壓力曲線Fig.11 Pressure curve of B1

圖12 B1點頻譜圖Fig.12 Spectrogram of B1

表7為兩種擴張比噴管的推力主要參數。從表7中可以看到：兩個擴張比噴管均呈現周期性特征；η=14.0的噴管振蕩頻率低于η=3.4的噴管，后者周期短，壓力振蕩頻率高；η=14.0的噴管呈現明顯的脹鼓現象，出現大范圍的高壓區；η=3.4的噴管周期性特點較弱，脹鼓程度小，因此脈動幅值小，射流穩定；η=3.4的噴管負推力小，平均推力大。

表6 B1、B2、B3、B4點主頻

表7 η=3.4與η=14.0推力主要參數

3 結論

本文對比了兩種擴張比噴管的流場結構與推力變化，并分析了二者之間的關系。得到結論如下：

1) 固體火箭發動機水下射流分為初始膨脹時期和周期性發展時期，其中周期性發展時期包括頸縮、脹鼓、回擊3個階段：頸縮階段，激波由噴管外向管內移動，噴管口處壓力在頸縮末期達到峰值，此時動量推力最小；脹鼓階段，激波向外移動，噴口速度增加，動量推力上升，壓差推力下降； 出口處達到穩定的高速射流；回擊階段，出口壓差推力不變，動量推力小幅度震蕩。

2) 不同擴張比的噴管射流形態不同。η=3.4的噴管各階段射流氣體形貌變化較弱，η=14.0的噴管各階段形貌特性明顯；η=3.4的噴管出口維持著穩定的超音速射流，而η=14.0的噴管在頸縮階段出口速度為亞音速；η=3.4的噴管高速氣體射流氣液邊界的相對運動較穩定，氣- 液摻混較為緩和。

3) 不同擴張比的噴管推力特性不同。η=3.4的噴管脹鼓程度小，因此脈動幅值小，射流穩定；η=3.4的噴管負推力小，平均推力大；η=3.4的噴管蕩頻率高于η=14.0的噴管，因此工作噪聲衰減更快、隱蔽性更好。故η=3.4的噴管比η=14.0的噴管更有利于水下推進。