基于RANS/NLAS的火箭跨聲速脈動壓力環境預示①

任淑杰，張收運，閆桂榮

(西安交通大學航天航空學院，西安 710049)

0 引言

脈動壓力又稱為氣動噪聲，在飛行器設計階段具有重要地位。目前，飛行器表面脈動壓力環境的預測主要采用風洞實驗和數值計算方法。文獻[1－3]均采用風洞實驗方法，研究了飛行器表面的脈動壓力特性;文獻[4－7]分別采用混合 RANS/LES、LES/FWH、LES/APE、LES/Lighthill計算方法，進行了脈動壓力環境預示。但風洞實驗耗資大，且重復性較差，而現有數值計算方法均采用LES求解流場，LES對計算網格要求很嚴，造成流場計算耗時。

為實現對飛行器脈動壓力環境的快速求解，Morris等人在1996年首次提出了一種采用雷諾平均N-S方程(RANS)求解流場，再采用非線性噪聲求解方程(NLAS)求解聲場相結合的混合方法，并用于超聲速來流三維軸對稱噴流噪聲的預示，但忽略了粘性擾動項的影響[8]。Batten等人在2002年進行了改進，考慮了粘性擾動項的影響[4，9]，并用于汽車后視鏡和密封腔氣動噪聲的預測。此方法采用RANS求解流場，使得流場的計算網格較LES更稀疏，節約了計算資源，縮短了計算時間，但尚未用于大型結構的脈動壓力環境預示研究。

對于大型火箭結構，跨聲速階段的脈動壓力環境相當嚴重，且有明顯的低頻峰值信號。低頻脈動壓力可能與火箭箭體結構低階模態接近，進而導致火箭箭體結構的疲勞甚至破壞。因此，準確預測火箭表面的脈動壓力，對火箭結構安全至關重要。

然而，直接將改進的RANS/NLAS方法用于大型火箭結構的脈動壓力環境預示時，由于火箭結構尺寸較大，會造成流場和聲場求解區域較大。在實際計算時，難免遇到計算量大、計算時間長的問題。此外，要精確描述NLAS方程中的非線性項，進行RANS方程求解時，需采用一種考慮非線性的湍流模型。

針對這些關鍵問題，本文采取如下技術方案:削減了部分遠場計算區域，將遠場設置為吸收邊界，以消除遠場邊界反射;減小了近壁網格尺度，并采用壁面函數法求解壁面區域，以提高壁面計算精度;選用兩方程非線性k-ε湍流模型，以精確求解NLAS方程的非線性項;最后，進行了大型火箭跨聲速脈動壓力環境的預示研究，計算結果較好地反映了火箭所處的脈動壓力環境，為火箭跨聲速脈動壓力環境預示提供了一定的參考依據。

1 控制方程和數值計算方法

NLAS的控制方程是從Navier-Stokes方程導出的，它將每一項分成統計平均項和擾動項2部分，即φ=ˉφ+φ＇。代入Navier-Stokes方程，并重新安排擾動項和平均項，得到非線性擾動方程組[8－9]:

式中Q＇是瞬態擾動項是瞬態平均項;F＇i是線性無粘擾動項;F＇ni是非線性無粘擾動項;是無粘平均項;是粘性擾動項;是粘性平均項。

具體表達式如下:

式中 ρ是密度;ui(i=1，2，3)是x、y、z3 個坐標軸方向的速度;p是壓強;e是單位體積能。

剪切應力項

熱傳導項

RANS控制方程采用有限體積法求解，無粘項采用二階精度TVD格式離散，粘性項采用中心差分格式離散，時間推進選用隱式方法。湍流模型選用兩方程非線性k-ε湍流模型，以精確求解NLAS方程的非線性項。遠場邊界給定來流的壓強p∞、溫度T∞、速度u∞。NLAS控制方程采用有限差分方法求解，時間推進也選用隱式方法。遠場邊界為RANS求得的平均流場信息。

流場物面第1層網格尺度設置為2×10－5，聲場物面第1層網格尺度設置為10－4，物面邊界均采用粘性無滑移絕熱壁，并應用壁面函數法求解，以保證物面區域求解精度。遠場計算區域均為火箭長度的5～10倍，并設置為吸收邊界，以消除遠場邊界反射。

2 計算結果及其分析

2.1 計算模型

計算了以錐柱-船柱-裙柱為基本外形的火箭表面的脈動壓力，自由來流馬赫數分別為 0.8、0.825、0.85、0.875、0.9、0.925、0.95、0.975、1.025、1.075、1.1、1.125、1.15、1.175、1.2。圖1 給出了火箭幾何外形簡圖，并記火箭表面6 個轉折點為A、B、C、D、E、F。

圖1 火箭幾何外形簡圖Fig.1 Schematic diagram of rocket

圖2給出了錐柱肩部脈動壓力系數計算曲線，并與文獻[10]中的曲線進行了對比。由圖2可知，計算值和文獻值基本吻合。

圖2 錐柱肩部脈動壓力系數Fig.2 Fluctuating-pressure coefficient of cone-cylinder

2.2 脈動壓力產生機理分析

脈動壓力環境與飛行器表面的基本流動狀態密切相關。經分析速度矢量圖(圖3)和壓強云圖(圖4)可知，跨聲速階段，錐柱肩部和裙柱部區的氣流始終未發生分離，這是因為氣流是否發生分離與半錐角、裙壓縮角、裙前柱長度等密切相關;來流Ma≥0.875時，船尾倒錐區的氣流發生分離，Ma=0.975時，氣流分離最嚴重。

錐柱肩部和裙柱部區的脈動壓力是由轉折點附近的激波振蕩引起的，船尾倒錐區的脈動壓力主要是由轉折點附近的激波振蕩與繞流分離再附體相互作用引起的，且此時形成更加強烈的脈動壓力。

2.3 錐柱肩部脈動壓力

跨聲速階段激波處于不穩定狀態。隨著來流馬赫數增加，激波后移，強度減弱，激波前后壓強差減小，由激波振蕩引起的脈動壓力也就相應地減小，這從圖5錐柱肩部脈動壓力系數隨馬赫數增加峰值后移、幅值減小可看出。從脈動壓力系數量級上，可發現錐柱肩部產生抖振的馬赫數范圍是0.8～0.9，具有強局部特性，馬赫數大于0.9之后，脈動壓力系數曲線雖有峰值，但量級很小。

圖3 船尾倒錐區速度矢量圖Fig.3 Velocity vector of inverted cone

圖4 壓強云圖Fig.4 Pressure contours of rocket

2.4 船尾倒錐脈動壓力

對于各種實際類型的分離流場，其再附體點比分離點產生更大一些的脈動壓力[10]。然而，當分離流動和激波振蕩相互作用時，這個結論就不再適用，而是與來流馬赫數密切相關。

圖5 錐柱肩部脈動壓力系數曲線Fig.5 Fluctuating-pressure coefficient of cone-cylinder

來流馬赫數為0.8～0.85時，船尾倒錐區的脈動壓力主要由轉折點B、C附近的激波振蕩產生，且轉折點B處的脈動壓力系數比轉折點C處的要高。

來流馬赫數為0.875～1.075時，船尾倒錐區脈動壓力主要由激波振蕩和流動的分離再附體相互作用引起。來流馬赫數為0.875時，由于氣流剛開始分離，分離對脈動壓力的影響小于激波振蕩的影響，分離點的脈動壓力系數仍高于再附體點;來流馬赫數為0.9～0.975時，分離對脈動壓力的影響和激波振蕩的影響比重相同，分離點和再附體點的脈動壓力系數大小基本一致，此時脈動壓力環境最為嚴重;來流馬赫數為1.025～1.075時，分離對脈動壓力的影響大于激波振蕩的影響，分離點脈動壓力系數明顯低于再附體點。來流馬赫數為1.1～1.2時，分離點后移速度加快，分離點未和轉折點B附近的激波振蕩相互作用，以致形成3個脈動壓力峰值:第1個是由激波振蕩引起的;第2個是由分離點引起的;第3個是由再附體點和激波振蕩相互作用引起的。

2.5 裙柱部區脈動壓力

裙柱部區的脈動壓力主要是由轉折點D、E、F附近的弱激波振蕩引起的。來流馬赫數小于1時，E點附近激波強度最大，相應的脈動壓力系數也就最大。來流馬赫數大于1之后，E點附近的激波迅速衰減，D點附近的脈動壓力系數最大。因此，馬赫數為0.8～1.0時，裙柱部區脈動壓力環境最為嚴重。

觀察圖5～圖7脈動壓力系數量級可知，火箭表面脈動壓力環境最嚴重的區域就是船尾倒錐區。

2.6 典型點的聲壓級頻譜

圖8給出了Ma=0.975時轉折點C、E處聲壓級的1/3倍頻程頻譜。由圖8可看出，跨聲速階段脈動壓力能量主要集中在低頻(100 Hz附近)，且在90～100 Hz范圍內的聲壓級最大。C轉折點總聲壓級為164.7 dB，是由繞流分離和激波振蕩相互作用產生的，E轉折點總聲壓級為144.3 dB，是由激波振蕩產生的。

圖6 船尾倒錐區脈動壓力系數曲線Fig.6 Fluctuating-pressure coefficient of inverted cone

圖7 裙柱部區脈動壓力系數曲線Fig.7 Fluctuating-pressure coefficient of skirt-cylinder

圖8 聲壓級的1/3倍頻程頻譜Fig.8 1/3 octave band am p litudes of SPL

3 結論

(1)基于RANS/NLAS，并應用兩方程非線性k-ε湍流模型、遠場吸收邊界及壁面函數法，可成功地進行大型火箭跨聲速脈動壓力環境預示。

(2)脈動壓力環境與飛行器表面的基本流動狀態密切相關。錐柱肩部和裙柱部區的脈動壓力主要是由激波振蕩引起的，船尾倒錐區的脈動壓力主要是由激波振蕩與繞流分離再附體相互作用引起的。

(3)船尾倒錐區的脈動壓力環境比錐柱肩部、裙柱部區更為嚴重。當馬赫數為0.9～0.975時，船尾倒錐區的脈動壓力環境最為嚴重。

(4)傳統觀點認為，對各種實際類型的分離流場，其再附體點比分離點產生更大的脈動壓力。然而，當分離流動和激波振蕩相互作用時，這結論就不再適用，而是與來流馬赫數密切相關。

(5)在跨聲速階段，脈動壓力能量主要集中在低頻(100 Hz附近)。

[1] 王娜，高超.彈體脈動壓力特征的實驗研究[J].實驗流體力學，2010，24(1).

[2] 操小龍，羅金玲，周丹杰，等.錐-柱體外形脈動壓力及抖振載荷響應研究[J].戰術導彈技術，2010(1).

[3] 楊青，李勁杰，楊永年，等.邊條翼布局主要參數對其雙垂尾抖振響應影響的風洞實驗研究[J].西北工業大學學報，2006，24(3).

[4] Paul Batten，UrielGoldberg，Sukumar Chakravarthy.Reconstructed sub-grid methods for acoustics predictions at all reynolds numbers[R].AIAA 2002-2511.

[5] Ali Uzun，Yousuff M Hussaini.Noise generation in the nearnozzle region of a chevron nozzle jet flow[R].AIAA 2007-3596.

[6] Elmar R Groschel，Matthias Meinke，Wolfgang Schroder.Noise prediction for a turbulent jetusing an LES/CAAmethod[R].AIAA 2005-3039.

[7] Rembold B，Kleiser L.Noise prediction of a rectangular jet using large-eddy simulation[J].AIAA 2894-700.

[8] Philip JMorris，Lyle N Long，Ashok Bangalore，et al.A parallel three-dimensional computational aeroacoustics method using non-linear disturbance equations[J].Journal of Computational Physics，1997，133:56-74.

[9] Paul Batten，Enrico Ribaldone，Mauro Casella，et al.Towards a generalized non-linear acoustics solver[R].AIAA 2004-3001.

[10] Robertson JE.The prediction of fluctuating-pressure environments(induced by three-dimension protuberances included)[R].Wyle Laboratories Research Staff ReportWR 71-10，March，1971.