固體火箭發(fā)動機(jī)水下超音速射流數(shù)值研究

王利利， 劉影， 李達(dá)欽， 吳欽， 王國玉

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

0 引言

潛射導(dǎo)彈具有隱蔽性好、系統(tǒng)簡便靈活、體積小、質(zhì)量輕等優(yōu)點(diǎn)，作為可靠的二次核打擊力量，越來越受到各國的重視[1]。固體火箭發(fā)動機(jī)由于結(jié)構(gòu)簡單、可靠性高、快速反應(yīng)能力強(qiáng)等優(yōu)點(diǎn)，逐漸成為潛射導(dǎo)彈的首選動力源[2]。然而，由于發(fā)動機(jī)水下工作環(huán)境介質(zhì)密度遠(yuǎn)大于空氣，且處于高溫、高壓、超音速工作狀態(tài)，其射流過程十分復(fù)雜[3]：一方面，大量的高溫燃?xì)鈬娙胨校瑫趯?dǎo)彈尾部形成瞬態(tài)燃?xì)馀?，高溫燃?xì)馀菖c水之間存在氣- 水摻混、換熱和相變等復(fù)雜物理現(xiàn)象并伴有激烈的流體動力干擾[4]；另一方面，燃?xì)馀葜写嬖诘膹?fù)雜激波結(jié)構(gòu)[5-6]以及激波結(jié)構(gòu)的運(yùn)動[7]均會使流場壓力劇烈變化，從而導(dǎo)致推力性能難以預(yù)測。近年來，為了提高發(fā)動機(jī)動力系統(tǒng)的推力性能，其采用的噴管擴(kuò)張比越來越大，導(dǎo)致其穩(wěn)定性大大降低。因此，深入研究大擴(kuò)張比水下發(fā)動機(jī)的工作特性及其相關(guān)兩相流動機(jī)理，對水下發(fā)射技術(shù)發(fā)展具有重要意義。

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

水下超音速氣體射流過程的推力特性及影響因素一直是研究的重點(diǎn)與難點(diǎn)。王寶壽等[16]通過試驗測量了噴管擴(kuò)張比為3、滯止壓強(qiáng)為9 MPa、氣體初始溫度為2 700 K、水深為10～40 m條件下固體火箭發(fā)動機(jī)水下工作時的推力和側(cè)向力,獲得了推力、側(cè)向力、推力矢量角隨發(fā)動機(jī)噴口堵塞面積比、水深的變化規(guī)律。張有為[3]利用球形氣泡模型對發(fā)動機(jī)水下推力特征的影響因素進(jìn)行了研究，結(jié)果表明，推力峰值隨著燃?xì)馀莩跏及霃降脑龃蠖杆俳档?，隨著燃面面積的減小而降低。

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

1 數(shù)值計算方法

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

將氣體假設(shè)為理想氣體，水為不可壓縮相，采用流體體積方法對氣- 液運(yùn)動界面進(jìn)行追蹤。連續(xù)方程、動量方程、能量方程、體積輸運(yùn)方程和狀態(tài)方程分別如下：

(1)

(2)

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

(3)

式中：氣- 液兩相體積分?jǐn)?shù)關(guān)系為αg+αw=1；能量En和溫度Tt為各項的質(zhì)量平均變量；γeff為有效熱傳導(dǎo)率；Sh為源項。

(4)

(5)

式中：R為理想氣體常數(shù)。

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

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

(6)

(7)

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

(8)

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

(9)

式中：模型常數(shù)C3=1.0.

由(9)式可知，當(dāng)湍流尺度較大時，湍流黏性系數(shù)表達(dá)為

(10)

1.3 計算模型及邊界條件

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

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

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

圖3 網(wǎng)格劃分結(jié)果及其放大圖Fig.3 Mesh generation and its enlarged view

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

表1 噴管工作參數(shù)

1.4 數(shù)值結(jié)果驗證

本文根據(jù)文獻(xiàn)[9]試驗結(jié)果進(jìn)行數(shù)值方法驗證，噴管結(jié)構(gòu)尺寸與工作條件均與文獻(xiàn)[9]一致。

表2為射流發(fā)展過程形貌對比。從表2中可以看出，數(shù)值模擬結(jié)果與試驗結(jié)果呈現(xiàn)出相同的周期性變化特征，數(shù)值計算結(jié)果對頸縮位置、脹鼓泡和回?fù)袅鞯牟蹲骄c文獻(xiàn)[9]較為一致。

表2 射流發(fā)展形貌的數(shù)值模擬結(jié)果與試驗結(jié)果[9]對比

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

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

1.5 推力計算方法

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

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

考慮背壓pb后的發(fā)動機(jī)推力計算公式為

(11)

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

2) 壓差推力Ae(pe-pa)，由噴管排氣面上的壓強(qiáng)pe與外界壓強(qiáng)pa之間的不平衡產(chǎn)生，取決于排氣截面的燃?xì)鈪?shù)pe、Ae和發(fā)動機(jī)所處的環(huán)境壓強(qiáng)pa.

3) 壓差推力(As-Ae)(pb-pa)，由作用于發(fā)動機(jī)前后端面上的環(huán)境壓強(qiáng)pa與背壓pb之間的壓力差產(chǎn)生。

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

(12)

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

2 結(jié)果與討論

2.1 不同擴(kuò)張比噴管射流瞬態(tài)演化過程

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

表3 射流發(fā)展過程云圖

Tab.3 Water-vapor contour of jet evolution process

由表3可見：

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

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

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

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

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

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

對比擴(kuò)張比η=14.0與η=3.4的噴管射流形貌可知：η=3.4的噴管在水氣輪廓圖上未發(fā)生較大變化，頸縮階段沒有出現(xiàn)明顯的主射流通道變窄現(xiàn)象，這是因為擴(kuò)張比小，氣體離開噴管后延徑向擴(kuò)張的慣性較小，故仍在軸向保持較大的動能與壓力能；η=3.4的噴管出口維持穩(wěn)定的超音速射流，除了在頸縮點(diǎn)出現(xiàn)小范圍高壓區(qū)，出口一直維持高速低壓狀態(tài)，而η=14.0噴管在頸縮階段出口速度為亞音速；η=3.4噴管出口處沒有觀察到紊亂的流線，且該處氣體徑向膨脹的趨勢小，故流線未呈輻射狀。脹鼓階段在水氣輪廓圖上并未出現(xiàn)明顯脹鼓泡，出口速度增加，由高壓轉(zhuǎn)為低壓，回流產(chǎn)生在液相區(qū)。在回?fù)綦A段出口速度增加，出現(xiàn)大范圍低壓區(qū)；從流線圖上可以看到，回流并未直接回?fù)舻焦鼙?，表明該噴管回?fù)魩淼膲毫φ袷幣c腐蝕較小。高速氣體射流氣- 液邊界的相對運(yùn)動較穩(wěn)定，氣- 液摻混較為緩和。

2.2 不同擴(kuò)張比噴管的推力特性

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

圖6 監(jiān)測點(diǎn)A位置Fig.6 Location of Point A

圖7 監(jiān)測點(diǎn)A壓力變化曲線Fig.7 Pressure at Point A

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

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

2.3 不同擴(kuò)張比噴管背壓振蕩特性

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

圖10 監(jiān)測點(diǎn)位置及放大圖Fig.10 Locations of monitoring points

圖11 監(jiān)測點(diǎn)B1壓力曲線Fig.11 Pressure curve of B1

圖12 B1點(diǎn)頻譜圖Fig.12 Spectrogram of B1

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

表6 B1、B2、B3、B4點(diǎn)主頻

表7 η=3.4與η=14.0推力主要參數(shù)

3 結(jié)論

本文對比了兩種擴(kuò)張比噴管的流場結(jié)構(gòu)與推力變化，并分析了二者之間的關(guān)系。得到結(jié)論如下：

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

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

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