固體火箭發(fā)動機噴管擴張段型面參數(shù)對其性能影響仿真分析①

穆 旭，田維平，董新剛，王德鑫，褚佑彪

(1.中國航天科技集團有限公司四院四十一所，西安 710025；2.中國航天科技集團有限公司第四研究院，西安 710025)

0 引言

噴管是固體火箭發(fā)動機的能量轉換裝置，將燃燒室高溫高壓燃氣的熱能轉換為動能產生推力。噴管擴張段型面影響燃氣流動及壁面壓力，關系噴管推力效率，選取合適的型面參數(shù)是發(fā)動機工程設計的重要步驟[1-2]。固體火箭發(fā)動機推進劑燃燒后生成顆粒(質量分數(shù)達30%)導致噴管內實際流動為復雜的氣體-顆粒兩相湍流，須采用兩相流仿真[3-4]，且燃氣溫度、壓強等參數(shù)在噴管內變化較大，豐富的流場特性和機理也值得探究。噴管型面設計需考慮發(fā)動機推力要求、燃氣熱力學特性、外界環(huán)境等因素，是固體火箭發(fā)動機設計的前沿課題，國內外學者對此開展了大量研究。

噴管最初依靠試驗和經(jīng)驗公式設計，但是受到發(fā)動機尺寸和研究成本等方面的限制，數(shù)值計算方法開始逐漸應用到噴管型面設計。RAO最早采用變分法在擴張段長度和質量流量一定的條件下設計了最大推力噴管型面，但未考慮兩相流損失[5]。后來HOFFMAN基于各種假設通過氣-粒兩相流理論計算分析了噴管性能與型面的關系并提出了一種軸對稱噴管型面最優(yōu)設計方法，未進行仿真且設計復雜，不適于工程應用[6-7]。ALLMAN采用直接優(yōu)化法比較一系列初始擴張角不變的二階多項式型面噴管的效率損失，選取性能最佳的噴管并與RAO的最大推力型面對比，誤差較小，但也采用理論計算方法[8]。

隨著數(shù)值方法和計算機的發(fā)展，計算流體力學被廣泛應用到噴管型面設計，縮短了研究周期且節(jié)約經(jīng)費。王成軒通過一維兩相平衡流計算設計了擴張段長度最短的固體火箭發(fā)動機最佳內型面[9]，研究主要聚焦理論分析，忽略了徑向分量影響。陳林泉通過對固發(fā)噴管流場數(shù)值計算，研究了噴管收斂段與喉部對流量的影響，但未考慮兩相流影響[10]，后來研究了擴張段長徑比對發(fā)動機比沖的影響[11]，但未涉及擴張半角的影響探究。王一白等研究了初始膨脹圓弧半徑、入口角和出口角對拋物線型噴管流動分離的影響[12]，對流場特性做出細致分析，對推力等性能的分析較少。上述研究提出了各種型面設計優(yōu)化方法，探究了不同型面噴管流場參數(shù)的變化規(guī)律，為固體火箭發(fā)動機噴管設計提供了堅實的理論基礎和設計準則，對工程應用有重要參考價值，但未全面系統(tǒng)地結合流場激波等特性探究擴張段型面參數(shù)對噴管性能的影響，存在一定局限性。

本文以橢圓-三次曲線型固體火箭發(fā)動機噴管為研究對象，采用數(shù)值計算方法，考慮流場激波影響，通過兩相湍流的性能變化規(guī)律研究和流場特征分析，從初始擴張半角和出口半角改進基準型面噴管，進一步探究擴張段長徑比及擴張比的影響，對噴管型面多維優(yōu)化和工程設計具有一定參考意義。

1 計算模型和方法

1.1 物理模型

根據(jù)仿真模擬和實際工況，噴管內燃氣參數(shù)軸向分量遠大于徑向及周向分量，且三維對稱模型與二維軸對稱模型流場分布和出口參數(shù)差異很小。本文研究高空飛行工況的上面級發(fā)動機噴管，出口靜壓大于空氣背壓，燃氣始終處于欠膨脹狀態(tài)，仿真對比也驗證了有無尾流場對出口截面及噴管內流場參數(shù)的影響較小，見圖1。另外，為精確模擬噴管流動，在收斂段上游適當增加一段燃燒室區(qū)域[13]。綜上，為提高計算速度，最終采用圖2所示幾何模型，噴管型面為橢圓-三次曲線型。收斂段為橢圓(式(1))，參數(shù)有收斂比和收斂段長度；喉部為短圓柱段，參數(shù)有直徑和長度；擴張段分為初始擴張段和主擴張段，初始擴張段為圓弧(式(2))，與喉部柱段和主擴張段相切；主擴張段為三次曲線(式(3))，設計變量包括初始擴張半角、出口半角、長徑比和擴張比。為驗證計算可靠性，以某噴管試驗結果作對比，將其型面設為參考基準型面(見表1)，保持收斂段及喉部參數(shù)不變，分別對各擴張段型面參數(shù)擴展范圍后仿真分析。

(1)

(X-X0)2+(Y-Y0)2=R2

(2)

Y=AX3+BX2+CX+D

(3)

圖1 有/無尾流場馬赫數(shù)云圖

燃氣熱傳導率由動能理論計算，粘性系數(shù)采用三系數(shù)Sutherland定律(式(4))計算；顆粒為慣性顆粒，初始溫度和速度與入口流體相同，直徑服從Rosin-Rammler分布(式(5))，均勻性指數(shù)取3.5，最大、最小及平均直徑分別為100、1、50 μm，其概率密度分布見圖3。

(4)

(5)

圖2 噴管幾何模型圖

表1 參考基準型面參數(shù)

圖3 顆粒直徑概率密度分布

1.2 控制方程

采用歐拉-拉格朗日法求解兩相流場，氣相在歐拉坐標下基于雷諾平均N-S方程求解：

(6)

(7)

顆粒相在拉格朗日坐標系下采用DPM模型求解，X方向顆粒單位質量受力平衡方程為

(8)

(9)

(10)

式中FD(u-up)為曳力項；u、up分別為氣相及顆粒速度；μ為流體分子粘度；ρ、ρp分別為流體及顆粒密度；dp為顆粒直徑；Fx為附加力項。

1.3 計算方法與邊界條件

計算采用壓力基求解器，壓力-速度以Coupled算法耦合求解，壓強、密度及能量等變量離散格式為二階(迎風)格式[14-15]，利于提高計算精度。湍流求解采用標準k-ε二方程模型和標準壁面函數(shù)；DPM模型中以隨機軌道模型考慮顆粒湍流分散，忽略顆粒間作用及顆粒對流場湍流渦影響；量級分析證明顆粒所受附加質量力、壓力梯度力、重力等相較曳力可忽略。研究高空飛行工況下噴管性能參數(shù)，流場邊界條件：

(1)質量流入口(Mass-flow Inlet)?？倻?700 K，純氣相及兩相中燃氣的質量流率分別為95 kg/s和66.5 kg/s，湍流強度5.1%，湍流粘度比3.7，來流方向垂直入口邊界，離散相邊界類型為逃逸。

(2)壓力出口(Pressure outlet)。靜壓22 Pa，湍流強度6.0%，湍流粘度比270.3，回流總溫3700 K，方向垂直出口邊界，離散相邊界類型為逃逸。

(3)固體壁面(Wall)。無滑移絕熱固壁，離散相邊界類型為沉積。

(4)軸線邊界(Axis)。軸對稱邊界條件。

1.4 模型驗證

1.4.1 網(wǎng)格無關性驗證

為確定合理的計算網(wǎng)格分布，兼顧計算準確性和效率，對初始擴張半角為38°、出口半角為9°、擴張比為48、長徑比為1.2的噴管分別采用不同數(shù)量的網(wǎng)格作單相流仿真。圖4為噴管出口馬赫數(shù)及推力隨網(wǎng)格數(shù)變化曲線，當網(wǎng)格數(shù)超過80 000后，馬赫數(shù)和推力變化很小，不同網(wǎng)格的流場分布也相當接近(圖5為馬赫數(shù)云圖比較)，可證明仿真結果網(wǎng)格無關性。故后續(xù)計算采用網(wǎng)格數(shù)為80 000的網(wǎng)格，具有較小誤差和較好收斂性，網(wǎng)格為結構化網(wǎng)格，在壁面附近局部加密，見圖6。

(a)Mach number (b)Thrust

圖5 兩種網(wǎng)格數(shù)下流場馬赫數(shù)分布比較

圖6 網(wǎng)格劃分示意圖

1.4.2 兩相流實驗驗證

對文獻[16]中的噴管進行氣體-顆粒兩相數(shù)值模擬并與實驗數(shù)據(jù)對比來驗證模型可靠性。計算模型及求解算法不變，噴管構型、邊界條件及顆粒參數(shù)與文獻一致，入口壓強為392 kPa，溫度為293 K，顆粒速度為20 m/s；出口壓強為101 325 Pa；氣體為理想氣體，顆粒為聚苯乙烯，直徑分別為545、1095 μm。圖7為噴管擴張段顆粒速度沿軸線分布，圓點和五角星表示文獻中兩種不同直徑顆粒的實驗結果，實線和虛線為對應的本文仿真結果，二者吻合良好，驗證了模型的可靠性。

圖7 顆粒速度分布對比

2 計算結果及分析

為研究各型面參數(shù)對噴管流場和性能的影響，采用控制變量法，擴充基準型面參數(shù)選取若干點分別數(shù)值模擬，參數(shù)見表2。

表2 噴管型面參數(shù)

先依次改進出口半角(Study 1)及初始擴張半角(Study 2)；再基于最佳型面分別研究長徑比(Study 3)和擴張比(Study 4)對性能的影響。通過噴管推力表征其性能[17]：

(11)

對各參數(shù)積分取平均值后，代入公式得各型面對應推力。

2.1 基準噴管流場仿真分析

基準噴管型面為某短程試驗噴管初始型面，發(fā)動機工作過程可忽略喉部燒蝕，擴張段型面基本不變。對其作純氣相及兩相流仿真，并同試驗數(shù)據(jù)對比(表3)，誤差較小，可驗證計算模型準確性。

表3 基準噴管結果對比

圖8為噴管純氣相馬赫數(shù)等值線圖，擴張段內激波ED和DB交于D點，導致BCD區(qū)域燃氣被二次壓縮，影響噴管性能。表明二者壓強差異主要在燃燒室，相差約6.93%，因為總質量流率不變，而兩相流顆粒對壓強無貢獻。

圖8 基準噴管純氣相馬赫數(shù)等值線圖

圖9為純氣相與兩相流燃氣參數(shù)對比。由圖9可知，純氣相流場變化更均勻。圖9(a)和圖10為兩相流顆粒參數(shù)，直徑越大的顆粒對流體跟隨性越差，更易保持慣性集中在軸線附近，與流體的阻力和熱交換也更大，此區(qū)域速度滯后和溫度滯后更嚴重，導致顆粒溫度偏高，速度偏低，與圖9(b)、(c)氣相分布規(guī)律相似。因此，兩相流同一截面參數(shù)不均勻度增加，噴管推力下降。

(a)Static pressure (b)Static temperature (c)Mach number

2.2 出口半角對噴管性能的影響

對表2中Study 1系列不同出口半角噴管分別仿真，圖11所示其參數(shù)變化規(guī)律和激波結構相同，僅分布存在差異，對噴管性能影響較小，出口參數(shù)及推力變化見圖12。

對比圖12(a)、(d)，噴管出口軸向速度決定動量推力，均隨出口半角增加先增大后減小，但兩相流極值點低于純氣相；對比圖12(b)、(e)，噴管出口靜壓決定靜推力，均隨出口半角增加而逐漸降低，且變化速率也隨之減小。圖12(c)所示噴管總推力包含圖12(d)、(e)分別對應的動量推力和靜推力，兩相流中還有圖12(f)中顆粒動量貢獻的推力，其中動量推力占比最大，在單相流及兩相流推力分別達到94%和83%，是噴管總推力性能的主導因素。

圖11 不同出口半角噴管馬赫數(shù)分布

圖12(c)變化趨勢顯示，出口半角大于極值點13°時，噴管推力隨出口半角增大而減小，因為出口半角越大，氣流偏離軸線程度加劇，燃氣徑向膨脹所受抑制減弱，軸向速度分量減小，推力減??；當出口半角過小時，噴管型面曲率變化更劇烈，壁面附近的燃氣在靠近出口區(qū)域受到壓縮，噴管性能反而下降，而且凝聚相顆粒容易碰撞壁面造成燒蝕。圖12(f)表明，兩相流出口截面上顆粒動量貢獻的推力隨出口半角增大而增大，說明燃氣對顆粒的曳力作用增強。

由上述分析可知，在Study 1系列中出口半角為13°的噴管性能最佳，將其同基準噴管的兩相流計算結果進行比較，見表4。由于最佳型面與初始型面較為接近，推力提升較小，也反映了出口半角對噴管推力的影響較小。在極值點13°附近，出口半角減小，單位角度比增大，單位角度對噴管推力的影響更大。

表4 基準型面與Study 1系列最佳型面兩相流性能對比

(a)Axial velocity (b)Static pressure (c)Total thrust

(d)Momentum thrust (e)Static thrust (f)Particle thrust

2.3 初始擴張半角對噴管性能的影響

初始擴張半角與噴管擴張段燃氣膨脹和激波結構相關，影響壁面壓強分布。仿真表2中Study 2系列不同初始擴張半角的噴管，同樣出口軸向速度決定動量推力，出口壓強決定靜推力，且動量推力起主導作用，其推力變化見圖13。圖13(a)表明，噴管推力隨初始擴張半角增大先逐漸增加，當初始擴張半角大于38°后又迅速減小，即存在最佳初始擴張半角，對應燃氣膨脹均勻，氣動損失較小，圖13(b)動量推力變化規(guī)律與之類似。由圖13(c)可知，隨初始擴張半角增加，噴管靜壓力先逐漸減小，達到最小值后又不斷增加，但兩相流對應的初始擴張半角極值點大于純氣相。圖13(d)表明，兩相流顆粒力隨初始擴張半角增大先減小后基本不變，說明燃氣對顆粒的曳力作用逐漸減弱。

(a)Total thrust (b)Momentum thrust

(c)Static thrust (d)Particle thrust

當初始擴張半角從38°開始減小，噴管擴張段曲率增大，型面趨近錐形，非軸向損失增加，噴管性能降低；初始擴張半角過小時，喉部附近內凹，性能迅速下降；但當擴張半角過大時，擴張段長度和擴張比固定不變，導致型面曲率半徑過大，容易發(fā)生流動分離，性能反而降低。就初始擴張半角而言，兩相流較純氣相推力降低約6.8%，但變化規(guī)律基本一致，因而求解最佳初始擴張半角時，可通過單相流計算確定大致范圍，再進行兩相流優(yōu)化。

從流場特征分析，圖14、圖15分別為初始擴張半角為25°、38°和50°的流場馬赫數(shù)云圖和出口軸向速度分布。

圖14 不同初始擴張半角下馬赫數(shù)云圖對比

圖15 不同初始擴張半角下出口軸向速度分布

當初始擴張半角較小(25°)時，內激波在噴管內部相交，B1-C1-D1區(qū)域燃氣被二次壓縮，B1-C1段出口軸向速度大幅降低，造成氣動損失；當初始擴張半角過大(50°)時，內激波更靠近壁面，激波強度增強且與軸線夾角增加，導致A3-B3-E3區(qū)域燃氣受壓縮程度加劇，B3附近出口軸向速度驟降，造成氣動損失；當初始擴張半角取最佳值(38°)時，僅產生一條內激波，強度較弱且基本與軸線平行，燃氣受壓縮強度低，B2-E2兩側燃氣參數(shù)梯度較小，氣動損失小。因此，為提高噴管性能，設計初始擴張半角時，應避免內激波在噴管內相交造成燃氣二次壓縮，同時盡量使內激波平行軸線，降低強度以減小氣動損失。綜上分析可知，study 2系列中初始擴張半角為38°的噴管性能最佳，同基準噴管兩相流計算結果比較見表5。

表5 基準型面與Study 2系列最佳型面兩相流性能對比

在極值點38°附近，初始擴張半角減小單位角度比增大單位角度對噴管推力的影響更大。而且初始擴張半角對噴管推力的影響高于出口半角，且推力增加主要因為占主導的動量推力增加，故設計時，應重點提高出口軸向平均速度。通過參數(shù)研究發(fā)現(xiàn)，本文計算條件下，初始擴張半角為38°，出口半角為13°的型面推力最大，提高了2.012 kN，可作為長徑比和擴張比探究的初始型面，也為后續(xù)型面多維優(yōu)化提供參考。

2.4 長徑比對噴管性能的影響

研究表2中Study 3系列不同擴張段長徑比的噴管性能差異，推力變化見圖16。圖16(a)、(b)說明長徑比小于1.2時，隨長徑比增加，噴管動量推力及總推力快速增加，因為噴管擴張段長度增加，型面曲線更加平滑，流體膨脹更充分；當長徑比為1.2～1.3時，推力隨之增長速度變緩；當長徑比為1.3～1.6時，噴管推力開始逐漸降低，存在最佳長徑比使噴管性能最佳，而兩相流對應的最佳長徑比小于純氣相。圖16(c)說明顆粒相使噴管兩相流出口壓強變化規(guī)律與純氣相不同，引起靜推力間的差異；圖16(d)顯示，顆粒力與噴管長徑比成正比。另外，噴管長度增加會導致質量增加，推重比下降，實際設計中需加以考慮。

(a)Total thrust (b)Momentum thrust

從流場特征分析，圖17、圖18分別給出了長徑比為1.2、1.4和1.6的流場馬赫數(shù)云圖和出口軸向速度分布。

圖17 不同長徑比下馬赫數(shù)云圖對比

圖18 不同長徑比下出口軸向速度分布

以長徑比1.2的圖為例，B1-C1-D1-E1為高馬赫區(qū)，A1-B1-E1為低馬赫區(qū)，不同長徑比噴管流場內激波形式相同，強度和位置存在差異，導致燃氣壓縮前高馬赫區(qū)與壓縮后低馬赫區(qū)占比不同。長徑比越大，內激波與出口截面交點越靠近軸線，激波強度也越強，導致交點附近出口軸向速度下降地越快，而且次激波強度也隨長徑比增大而增強，導致氣動損失加劇。但長徑比越大的噴管激波前高馬赫區(qū)未受壓縮的燃氣膨脹的更充分，靠近軸線段出口軸向速度更大。而噴管動量推力取決于整個出口截面的平均軸向速度，故存在最佳長徑比在1.3～1.4之間，使圖18中出口軸向速度曲線積分最大，對應噴管性能最佳。

2.5 擴張比對噴管性能的影響

圖19為表2中Study 4系列不同擴張比的噴管計算結果。保持長徑比不變，當噴管擴張比增加時，擴張段長度及出口面積隨之增加，燃氣膨脹更充分，因此噴管推力增加；圖19(b)展示了推力增長率隨擴張比變化的情況，在擴張比整個變化范圍內，單相流和兩相流推力變化分別達到3.325%和3.501%，說明擴張比的變化對噴管性能影響很大，而單相流和兩相流的增長趨勢和幅度較為接近；圖19(c)說明，出口處顆粒相對燃氣速度滯后和溫度滯后嚴重，且隨擴張比增加而增加，導致兩相流損失增加，噴管效率降低，但擴張比增加，使得噴管推力系數(shù)增加。綜合來看，擴張比越大，噴管推力越大。而實際設計時，還需考慮噴管質量增加，以確定性能最佳的噴管擴張比。

(a)Total thrust (b)Growth rate of thrust (c)Lag

3 結論

采用歐拉-拉格朗日數(shù)值方法，模擬了橢圓-三次曲線型噴管在擴張段不同型面參數(shù)下(出口半角、初始擴張半角、長徑比和擴張比)的兩相流場，并與試驗數(shù)據(jù)對比，驗證了計算模型和方法的可靠性。通過積分出口截面參數(shù)得到噴管推力各分量：動量推力、靜推力和顆粒力。采用控制變量法，對擴張段型面主要參數(shù)進行了逐一對比研究。

(1)噴管性能與流場變化緊密聯(lián)系：隨擴張段型面變化，流場內激波結構相應改變，尤其內激波相交會對燃氣二次壓縮，顯著影響噴管出口軸向速度。

(2)擴張段出口半角對噴管推力影響較小，而初始擴張半角的影響相對明顯。通過流場提取的參數(shù)對比表明，本文計算條件下，出口半角為13°，初始擴張半角為38°時，噴管性能最優(yōu)，其推力較基準型面提升2 kN。

(3)隨長徑比增大，噴管推力先增大后逐漸減小，長徑比小于1.2時，隨長徑比增加，推力收益增速明顯。擴張比越大，噴管推力系數(shù)越大，但顆粒滯后和兩相流損失隨之增大，導致噴管效率降低。綜合來講，噴管推力仍呈上升趨勢。