助推器頭部跨音速脈動壓力CFD數(shù)值仿真

沈 丹，蘇 虹，徐珊姝

(北京宇航系統(tǒng)工程研究所，北京100076)

1 引言

大型運載火箭往往采用捆綁式布局，火箭在主動段以跨聲速飛行時，在箭體外形突變處，會產(chǎn)生局部的非定常氣動力，即分離流和附體流交替的流場以及激波振蕩[1-3]。 基于這個原因，助推器頭部和整流罩都會在跨音速階段產(chǎn)生強烈的低頻脈動壓力。 其中，助推器的頭錐不但面臨自身“錐-柱”交界面形狀突變帶來的跨音速效應(yīng)，同時還面臨著與火箭芯級之間的氣動干擾，即助推器錐角處的激波在芯級發(fā)生反射。 大型捆綁火箭由于尺寸大，導(dǎo)致風洞模型的縮放比例較大，因此針對邊界層、激波-邊界層干擾等影響脈動壓力的局部物理現(xiàn)象難以模擬。

關(guān)于火箭脈動壓力的研究，針對整流罩開展的計算較多。 任淑杰等[4]通過求解流場、求解非線性聲場相結(jié)合的技術(shù)，數(shù)值求解了火箭跨聲速表面錐柱肩部、船尾倒錐、裙柱部區(qū)典型位置處脈動壓力。 趙瑞等[5]基于5 階WENO 格式構(gòu)造隱式大渦模擬方法，對跨聲速來流條件下火箭整流罩脈動壓力環(huán)境進行數(shù)值模擬，擬合出干擾區(qū)脈動壓力的經(jīng)驗公式，并提出一種脈動壓力環(huán)境的工程預(yù)測流程，在保證計算效率的同時，完善并提高工程經(jīng)驗公式的預(yù)測精度。 NASA 早在上世紀70 年就已摸索出一套關(guān)于火箭頭部脈動壓力的工程預(yù)示方法，并給出航天飛機脈動壓力環(huán)境[6]。 然而，這些研究并未關(guān)注助推器肩部的跨聲速脈動壓力以及助推器與芯級之間干擾對脈動壓力的影響。

本文使用計算流體力學(xué)(Computational Fluid Dynamics，CFD)方法數(shù)值求解大型捆綁火箭跨音速非定常流場，特別關(guān)注助推器頭部氣動干擾區(qū)域，針對直頭錐、斜頭錐2 種助推器頭部形狀，提取頭錐附近壁面壓力的時間歷程，獲取脈動壓力特性與馬赫數(shù)、頭錐形狀等參數(shù)之間的規(guī)律。

2 數(shù)值仿真

2.1 模型及網(wǎng)格

在全箭流場范圍內(nèi)構(gòu)造結(jié)構(gòu)化網(wǎng)格，拓撲結(jié)構(gòu)共分為20 個塊，每個助推器獨立成塊。 在全箭頭部用O 型網(wǎng)格拓撲伸展出6 個，尾部O 型網(wǎng)格拓撲伸展出6 個，助推器、芯級包裹6 個塊(4 個助推4 塊，1 個芯級1 塊，縫隙1 塊)，助推器外側(cè)2 塊。網(wǎng)格和拓撲結(jié)構(gòu)如圖1 所示。 計算網(wǎng)格數(shù)500 萬～700 萬，在壁面第一層網(wǎng)格尺寸為0.01 mm。

2.2 數(shù)值計算方法

為了提高模擬非定常流動隱式推進過程的時間精度，采用偽時間子迭代技術(shù)(Pseudo Time Sub-Iteration， τTS，也稱作雙時間推進技術(shù)，Dual Time Stepping)，將隱式高斯賽德爾(Lower-Upper Symmetric Gauss-Seidel，LU-SGS)格式的時間精度提高到2 階[7]。 根據(jù)經(jīng)驗并考慮計算精度和效率因素，取子迭代步數(shù)為3 步。

在湍流計算方面，采用基于剪切應(yīng)力輸運模型(Menter's Shear Stress Transport， SST)的雷諾平均/大渦模擬( Reynolds Averaged Navier-Stocks/Large Eddy Simulation，RANS/LES)混合方法，在高頻、小尺度運動占主導(dǎo)地位的近壁區(qū)域使用高效可靠的湍流模式，而在低頻大尺度運動占優(yōu)的非定常分離流動區(qū)域則使用LES 類方法[8-10]。 模擬跨聲速火箭助推器頭部的非定常激波/邊界層干擾流場，包括流動分離、激波振蕩等現(xiàn)象。

圖1 計算網(wǎng)格和拓撲結(jié)構(gòu)Fig.1 The mesh and topology structure

在并行計算方面，采用基于并行計算(Message Passing Interface，MPI)的數(shù)值算法，可以獲得更大規(guī)模的網(wǎng)格分辨率，以期得到更符合實際的數(shù)值預(yù)測結(jié)果。

2.3 工況及初邊值

按照運載火箭實際飛行彈道規(guī)律，挑選如表1 所示的工況作為主要的計算工況，隨馬赫數(shù)變化，攻角和雷諾數(shù)按照典型彈道得到。

表1 計算工況Table 1 Calculation cases

同時對直頭錐、斜頭錐2 種頭部構(gòu)型的助推器進行對比。 其中，直頭錐助推器廣泛應(yīng)用于我國現(xiàn)役火箭，而斜頭錐則應(yīng)用于我國新一代大型運載火箭。 本文中，直頭錐助推器的半錐角為15°，斜頭錐助推器的最大半錐角為30°。 2種幾何構(gòu)型的示意圖如圖2 所示。 工況1～6 的來流參數(shù)均是根據(jù)火箭飛行中的典型時刻確定的。

圖2 剖面示意Fig.2 Section plane position

圖2同時給出了A-J 共10 個壓力數(shù)據(jù)提取剖面的位置，A、B、C、D、G、H 為助推器上的6 個半圓，E、F 為助推表面的1/4 圓，I、J 為芯級表面的1/4 圓。 半圓上按順時針方向分別在0°、45°、90°、135°、180°5 個點布置測壓點，1/4 圓上按順時針方向分別在0°、45°、90° 3個點布置測壓點。

3 計算結(jié)果

3.1 流場特征

直頭、斜頭錐助推器的時均壓力云圖如圖3所示。 由圖可見，在跨音速時期助推器頭部的自由來流馬赫數(shù)小于1，氣流經(jīng)過頭錐壓縮增速，在拐角處接近音速，經(jīng)過拐角后膨脹加速，在肩前部產(chǎn)生1 道斜激波，在肩后部產(chǎn)生1 道正激波的結(jié)尾激波，2 道激波前后跳動振蕩。 所不同的是，直頭錐的拐角角度小，氣流經(jīng)過拐角時產(chǎn)生一系列波系與芯級壁面相互作用。 而對于斜頭錐，外側(cè)拐角角度大(是斜頭錐的2 倍)，氣流在拐角處膨脹嚴重，形成較大的分離，流動在超聲速附體流和亞聲速分離流之間反復(fù)交替[11]。

圖3 Ma=0.8 時均壓力云圖Fig.3 Pressure contour(Ma=0.8)

直頭、斜頭錐助推器的Q ＝-0.0001 等值面壓力云圖如圖4 所示。 其中Q 為速度梯度張量第2 不變量，其等值面可以用來表現(xiàn)三維流場的渦結(jié)構(gòu)。 由圖可見，直頭錐的肩部渦系對芯級表面影響很大，二者之間的干擾嚴重。 相對而言，斜頭錐除背風面以外，助推與芯級之間少有渦結(jié)構(gòu)出現(xiàn)，二者之間的干擾輕微。

3.2 均方根壓力系數(shù)

均方根脈動壓力數(shù)據(jù)也是火箭脈動壓力最重要的設(shè)計指標，將提供給載荷計算所用。 均方根脈動壓力系數(shù)按式(1)計算：

式中，PLU為局部非定常壓力，PLS為局部非定常壓力的時間平均值，T 為取樣時間，q 為來流動壓。 則針對上述仿真結(jié)果進行統(tǒng)計，得到各個工況下均方根脈動壓力系數(shù)的中位數(shù)值、峰值及其所在剖面，如表2 所示。

圖4 Ma=0.8 等Q 面壓力云圖Fig.4 Pressure contour of Q(Ma=0.8)

表2 各工況均方根脈動壓力系數(shù)統(tǒng)計Table 2 CP-rmsstatistics

由此可見，工況1(Ma ＝0.7、a ＝4.0°)時，直頭錐助推器首先迎來脈動特征。 此時雖然來流馬赫數(shù)僅有0.7，但助推器肩部局部已經(jīng)進入跨音速區(qū)域，跨音速激波不穩(wěn)，在助推肩部振蕩。 此時各剖面的脈動壓力分布規(guī)律如圖5 所示，圖中脈動壓力的高值都出現(xiàn)在助推與芯級間的狹縫處，是助推直頭錐肩部波系與芯級之間的干擾所致。

工況2(Ma ＝0.8、a ＝6.0°)時，直頭錐和斜頭錐均反映出一定的脈動壓力特性，且分布比較隨機。

圖5 工況1 直頭錐均方根脈動壓力分布規(guī)律Fig.5 Case 1，Prmsdistribution on positive nosecone

工況3(Ma ＝0.9、a ＝8.0°)開始，直頭助推器脈動壓力顯著減小，而斜頭助推器在此時迎來脈動壓力極值，最高值達0.19821，位于D 剖面的135°位置上，是助推器頭錐外側(cè)。 顯然是由助推器頭錐肩部折角處產(chǎn)生激波振蕩所主導(dǎo)的脈動，計算中各截面的激波振蕩頻率在15～50 Hz 之間。 此時各剖面的脈動壓力分布如圖6所示。 極值出現(xiàn)的截面D-135°測點上，出現(xiàn)了圖7 所示的無量綱功率譜密度，存在33 Hz 的主頻。

圖6 工況3 斜頭錐均方根脈動壓力分布規(guī)律Fig.6 Case 3，Prmsdistribution on oblique nose cone

圖7 工況3D-135°測點無量綱功率譜密度Fig.7 Case 3，Amplitude spectrum on D-135°

隨著馬赫數(shù)的進一步增加，2 種構(gòu)型的跨音速效應(yīng)均明顯減少，脈動壓力的均方根隨之減小。對于工況5 和工況6，無論直頭錐還是斜頭錐，所引起的脈動壓力均可忽略。

這一結(jié)果與文獻[6]中的結(jié)論是相符的，即：在一般的運載火箭上升飛行過程中，對于15°錐角的錐柱交界面，在Ma ＝0.74 左右發(fā)生由激波振蕩引發(fā)脈動壓力極值，而對于30°錐角的錐柱交界面，則在Ma ＝0.90 左右發(fā)生由激波振動引發(fā)脈動壓力極值。

從本文的計算可見，直頭錐助推器先于斜頭錐發(fā)生激波振蕩，按照一般的運載火箭軌跡，在Ma ＝(0.7，1.0)這一飛行時段，動壓的變化范圍并不大，均在20 000～25 000 N/m2量級左右，即來流的能量相差不大。 但15°錐角脈動壓力系數(shù)極值為0.119，30°錐角脈動壓力系數(shù)極值為0.198，二者相差顯著。 可見在飛行中，大的外折角出現(xiàn)脈動壓力極值的時間靠后，但強度更強。

除膨脹轉(zhuǎn)折角度的差異以外，直頭錐和斜頭錐還有1 個顯著不同，即與芯級的干擾關(guān)系。顯然，直頭助推器的折角靠近芯級，而斜頭助推器的平直段靠近芯級，前者與芯級表面的干擾更為顯著。 在來流速度不斷增加、全場完全進入超音速流動后，直頭助推器頭激波在芯級表面上發(fā)生反射，工況4～5 中，芯級表面(IJ 剖面)的脈動壓力顯著回升。 而對于斜頭錐，此時助推器肩部出現(xiàn)脈動極值，而芯級表面卻仍然十分平靜。

3.3 與試驗值的比對

符合工況2 和4 的斜頭錐外形風洞試驗在中國空氣動力研究與發(fā)展中心的FL-24 風洞進行，試驗使用壓阻式脈動壓力傳感器，傳感器布置在上述各截面處，采樣頻率為5 kHz。 將數(shù)值計算得到的均方根脈動壓力系數(shù)與試驗數(shù)據(jù)進行對比，如圖8 所示。

對于工況2(Ma ＝0.8)，計算與試驗值在量級和分布規(guī)律上比較接近。 而對于工況4，計算得到某些大峰值，這些峰值出現(xiàn)在背風助推器或助推器的背風面，而試驗值卻并沒有體現(xiàn)。 當Ma在0.9～1.0 時，背風面出現(xiàn)大的激波震蕩，這種震蕩由激波和邊界層干擾，對來流馬赫數(shù)和表面邊界層的細微變化十分敏感，因此在試驗中未必能夠捕捉到。

圖8 均方根脈動壓力系數(shù)計算值與試驗值比較Fig.8 Comparison of calculated and experimental CP-rms values

4 結(jié)論

1)針對直頭錐和斜頭錐2 種不同的助推器外形，由于外折的膨脹角度不同，迎來最大均方根脈動壓力系數(shù)的先后和量級不同。 直頭錐在Ma ＝0.7左右率先達到再附條件并實現(xiàn)激波振蕩，斜頭錐在Ma ＝0.9 左右才達到脈動壓力極值，但脈動壓力系數(shù)極值更大；

2)直頭助推器與芯級的干擾效應(yīng)更顯著，將在芯級表面由于激波反射引發(fā)脈動壓力；

3)隨著來流速度的提高，跨音速效應(yīng)明顯減少，脈動壓力均方根隨之減??；

4)計算和試驗的均方根脈動壓力數(shù)值符合良好，驗證了本文中數(shù)值算法在工程應(yīng)用中的可行性。

參考文獻(References)

[1]NASA.Buffeting during atmosphere ascent[R].NASA SP-8001，1964.

[2]鮑賢棟.對兩種頭部外形氣動噪聲的預(yù)測分析[J].強度與環(huán)境， 1980(2)：31-42.Bao X D.Prediction of aerodynamic noise on two different nose cone shape [J].Structure ＆ Environment Engineering，1980(2)：31-42.(in Chinese)

[3]沈丹，吳彥森，岑拯.芯級與助推器頭部氣動干擾流場數(shù)值模擬[J].導(dǎo)彈與航天運載技術(shù)， 2013(6)：42-46.Shen D，Wu Y S，Cen Z.Numerical simulation of aerodynamic interaction characteristics between the rocket and a booster's nosecone [J].Missile and Space Vehicles， 2013(6)：42-46.(in Chinese)

[4]任淑杰，張收運，閆桂榮.基于RANS/NLAS 的火箭跨聲速脈動壓力環(huán)境預(yù)示[J].固體火箭技術(shù)，2011，42(4)：418-422.Ren S J，Zhang S Y， Min G R.A prediction of fluctuation pressure conditions with transonic rocket by RANS /NLAS method [J].Journal of Solid Rocket Technology，2011，42(4)： 418-422.(in Chinese)

[5]趙瑞榮， 吉利， 王超，等.跨聲速條件下火箭整流罩表面脈動壓力環(huán)境分析與預(yù)測[C]/ /中國力學(xué)大會， 上海，2015：158.Zhao R R， Ji L， Wang C H，et al.Analysis and prediction of fluctuating pressure environment on the surface of rocket fairing under transonic condition [C]/ /CCTAM， Shanghai，2015：158.(in Chinese)

[6]Kenneth J P， Robertson J E.Prediction of space shuttle fluctuating pressure environments， including rocket plume effects[R].NASA-CR-124347.

[7]趙慧勇，樂嘉陵.雙時間步方法的應(yīng)用分析[J].計算物理， 2008， 25(3)：253-258.Zhao H Y， Le J L.Application analysis on dual-time stepping[J].Chinese Journal of Computational Physics， 2008， 25(3)：253-258.(in Chinese)

[8]肖志祥.湍流模型在復(fù)雜流場數(shù)值模擬中的應(yīng)用[J].計算物理， 2003，20(4)：335-340.Xiao Z X.Applications of turbulence models in simulation of complex flows [J].Chinese Journal of Computation Physics，2003， 20(4)：335-340.(in Chinese)

[9]肖志祥，陳海昕.采用RANS/LES 混合方法研究分離流動[J].空氣動力學(xué)學(xué)報， 2006，24(6)：218-222.Xiao Z X， Chen H X.Simulation of separation flows with RANS/LES hybrid methods [J].Acta Aerodynamica Sinica，2006， 24(6)：218-222.(in Chinese)

[10]肖志祥， 陳海昕.激波/邊界層干擾流動數(shù)值模擬格式的應(yīng)用研究[J].西北工業(yè)大學(xué)學(xué)報，2004， 22(6)：800-806.Xiao Z X， Chen H X.The simulation of shock/boundary interaction flows with several temporal and spatial schemes [J].Journal of Northwestern Poly Technical University， 2004， 22(6)：800-806.(in Chinese)

[11]操小龍.錐-柱體外形脈動壓力及抖振載荷響應(yīng)研究[J].戰(zhàn)術(shù)導(dǎo)彈技術(shù)，2010(1)：67-72.Cao X L.The Research of fluctuating pressure and buffeting load on the cone-column body[J].Tactical Missile Technology， 2010(1)：67-72.(in Chinese)