許東, 李欣, 徐力, 王亞輝

(1.北京航空航天大學 儀器科學與光電工程學院, 北京 100191; 2.北京航天自動控制研究所 宇航與智能控制技術(shù)國家級重點實驗室, 北京 100854)

球錐形頭罩高超聲速繞流流場數(shù)值模擬

許東1, 李欣1, 徐力2, 王亞輝2

(1.北京航空航天大學 儀器科學與光電工程學院, 北京 100191; 2.北京航天自動控制研究所 宇航與智能控制技術(shù)國家級重點實驗室, 北京 100854)

研究了不同飛行參數(shù)下球錐形光學頭罩的高超聲速繞流流場的數(shù)值模擬問題。首先,分別利用k-ε湍流模型和SA湍流模型求解N-S方程,獲得了球錐形光學頭罩繞流流場的模擬結(jié)果。然后,將模擬結(jié)果與風洞測試得到的實驗數(shù)據(jù)進行比較分析。結(jié)果表明，在高馬赫數(shù)、大迎角情況下采用k-ε湍流模型可以得到相對更準確的計算結(jié)果。最后,對流場的溫度、壓強、密度分布特性進行了分析,給出了球錐形頭罩繞流流場的壓強、溫度和密度隨飛行參數(shù)的變化規(guī)律。

球錐形頭罩; 流場; 數(shù)值模擬; 湍流模型

0 引言

隨著現(xiàn)代化戰(zhàn)爭對導彈飛行速度和突防能力要求的不斷提高,高超聲速導彈已成為當前的研究熱點。導彈頭罩與周圍氣體劇烈作用會形成復雜的流場,進而產(chǎn)生氣動光學效應[1],影響光學制導系統(tǒng)對目標的探測、識別和跟蹤的能力。球錐形光學頭罩具有較好的空氣動力學特性,可以減小導彈在高速飛行時受到的空氣阻力和摩擦熱,現(xiàn)已廣泛應用于空空導彈(如“霹靂”7、“楊樹”R-27T、“怪蛇”3等)。進行球錐形頭罩的繞流流場數(shù)值模擬,不僅對研究導彈的氣動特性具有重要意義,而且還是實現(xiàn)氣動光學效應校正、提高導彈制導精度的關鍵技術(shù)之一。

對流場參數(shù)計算的準確性在很大程度上取決于所選擇的湍流模型。目前對高速流場的數(shù)值模擬主要采用的是k-ε湍流模型和SA(Spalart-Allmaras)湍流模型,并且研究對象主要集中于飛行器整體、翼型、平板形側(cè)窗等[2-7]。如劉景源[4]、于子文[5]利用k-ε湍流模型分別對平板繞流、前飛旋翼流場進行了數(shù)值模擬,周大高[6]、段卓毅[7]利用SA湍流模型分別對翼型S825繞流、某機翼副翼舵面流場進行了數(shù)值模擬。這些研究都得到了與實驗數(shù)據(jù)吻合較好的計算結(jié)果，但是這兩種湍流模型對于球錐形頭罩高超聲速繞流流場是否適用并不確定。因此,本文將首先對這兩種湍流模型在計算球錐形頭罩高超聲速繞流流場數(shù)值模擬中的有效性進行討論;然后通過對模擬結(jié)果的分析,研究不同馬赫數(shù)和不同迎角條件下繞流流場的溫度、壓強、密度等參數(shù)的分布特性,為分析球錐形頭罩的氣動特性和氣動光學效應提供數(shù)據(jù)依據(jù)。

1 計算方法

1.1 控制方程

流體運動的基本方程是建立在質(zhì)量守恒定律、牛頓運動定律和能量守恒定律的基礎上的。習慣上將連續(xù)方程、動量方程和能量方程組成的流體運動方程稱為N-S方程組。在三維笛卡兒直角坐標系中,可壓縮粘性流體N-S方程組的守恒積分形式為:

(1)

式中,W為流場守恒矢量;V為任意控制體體積;?V為控制體邊界;H為矢量通量,可分解為無粘對流矢量通量Hc和粘性矢量通量Hv之差的形式;n為?V表面的單位外法矢量；S為面積。W,Hc,Hv的具體表達式可參見文獻[1]。

為封閉方程組,還需引入4個熱力學關系式，分別為完全氣體狀態(tài)方程:

p=ρRT

(2)

單位質(zhì)量氣體的總能量:

e=p/[(γ-1)ρ]+ujuj/2

(3)

由Sutherland公式計算得到的動力粘性系數(shù)為:

μ=C1T3/2/(T+C2)

(4)

式中,C1=1.458×10-6kg/m·s·K1/2;C2=110.4 K。

對于各向同性流體,其導熱系數(shù)為:

k=μCP/Pr=μγR/[(γ-1)Pr]

(5)

式中,γ為氣體比熱系數(shù),對于完全氣體通常γ=1.4;R為氣體常數(shù);Pr為普朗特數(shù),對于完全氣體通常Pr=0.72。

1.2 湍流模型

k-ε二方程湍流模型是由Launder等人在1972年提出的,因其穩(wěn)定性、經(jīng)濟性和計算精度較高，現(xiàn)已成為湍流模型中應用最廣泛、最為人熟知的。SA湍流模型是由Spalart和Allmaras在1992年提出的,是相對簡單的一方程模型,對流動分離區(qū)附近的計算效果較好,常用于模擬翼形、壁面邊界等流動,現(xiàn)已廣泛應用于航空領域。

在k-ε二方程湍流模型中,通過求解湍流動能k方程和湍動耗散率ε方程,得到k和ε的解,然后再用k和ε的值計算湍流粘度μt,最終通過Boussinesq假設得到雷諾應力的解。其相應的方程和函數(shù)表達式為[8-9]:

(6)

(7)

(8)

式中,模型系數(shù)Cμ=0.09,Cε1=1.44,Cε2=1.92,σk=1.0,σε=1.3。

(9)

式中,右邊四項分別為湍流粘性的生成項、耗散項、自定義源項和解體項。

(10)

(11)

式(9)～式(11)中的各項具體表達式和模型系數(shù)可參見文獻[10-11]。

1.3 網(wǎng)格劃分

圖1 計算網(wǎng)格示意圖Fig.1 Scheme of computational grids

1.4 邊界條件

上游邊界:采用無窮遠處自由可壓來流作為壓力遠場邊界條件,給定來流的靜壓、靜溫和馬赫數(shù);出口邊界:采用壓力出口邊界條件,對于達到高超聲速的流動,參數(shù)都由流場內(nèi)部通過插值外推得到;物面邊界:采用無滑移壁面邊界條件。

2 計算結(jié)果與分析

2.1 湍流模型有效性分析

以文獻[12]中風洞測試的實驗條件為輸入?yún)?shù),設來流的初始壓強為1.5 MPa,初始溫度為370 K,導彈飛行馬赫數(shù)為5和8,迎角α為0°,5°,10°,15°,對圖1中所示光學頭罩的繞流流場進行數(shù)值模擬。計算過程中分別采用k-ε湍流模型和SA湍流模型對N-S方程進行求解,并將計算得到的頭罩表面壓力系數(shù)與風洞實驗測得的頭罩表面壓力系數(shù)進行了比較。

圖2和圖3中給出了z=0對稱截面上頭罩表面壓力系數(shù)沿軸向距離的分布情況,其中Cp0表示風洞實驗測得的壓力系數(shù),Cp1表示采用k-ε湍流模型時計算得到的壓力系數(shù),Cp2表示采用SA湍流模型時計算得到的壓力系數(shù)。從圖中可以看出,當飛行馬赫數(shù)為5、迎角為0°～15°時,兩種模型的計算結(jié)果與實驗數(shù)據(jù)吻合程度都比較好;當飛行馬赫數(shù)為8、迎角為0°和5°時,兩種模型的吻合程度也都比較好;但當馬赫數(shù)為8、迎角為10°和15°時,k-ε湍流模型的計算結(jié)果比SA湍流模型的計算結(jié)果更為準確。這是由于SA模型運輸方程中的常系數(shù)是基于簡單流動分析得到的,在這些簡單流動中湍流運輸特性一般處于平衡狀態(tài)。但在高馬赫數(shù)、大迎角情況下湍流的非平衡運輸現(xiàn)象明顯,湍流粘性的生成和耗散之間的平衡關系發(fā)生了改變,導致計算結(jié)果出現(xiàn)了偏差?？梢娫陲w行條件較為苛刻的情況下,k-ε湍流模型比SA湍流模型更為適合。

圖2 Ma=5時壓力系數(shù)分布Fig.2 Distribution of the pressure coefficient when Ma=5

圖3 Ma=8時壓力系數(shù)分布Fig.3 Distribution of the pressure coefficient when Ma=8

2.2 流場參數(shù)的分布情況

以飛行高度為5 km、馬赫數(shù)分別為5和8、迎角分別為0°和15°的飛行參數(shù)為例,利用k-ε湍流模型對光學頭罩的繞流流場進行了數(shù)值模擬。圖4～圖6給出了z=0對稱截面上頭罩周圍流場的壓強、溫度、密度的等值分布云圖。從圖中可以看出,頭罩頂部區(qū)域的壓強、溫度、密度都最大,在Ma=5時壓強高達1.6 MPa,溫度高達1 500 K,密度高達4 kg/m3;在Ma=8時則分別高達4 MPa,3 500 K,4.5 kg/m3。這是由于來流與光學頭罩相遇時受到強烈壓縮,壓強、密度大為增加;空氣本身又具有粘性,與頭罩表面相接觸時受到阻滯,使得氣流的速度大為降低,在頭罩表面附近形成邊界層,在邊界層內(nèi)來流的動能幾乎全部轉(zhuǎn)化為內(nèi)能,使光學頭罩周圍的氣流溫度迅速升高。而頭罩后部的流場充分發(fā)展為湍流,形成弓形脫體激波,分子運動耗散掉湍動能,使得壓強、溫度、密度都有所降低。

同時從圖中還可以看出,隨著迎角的增加,駐點位置逐漸移至迎風面,頭罩弓形激波在迎風面更加貼近物面,脫體距離減少,同一位置對應的壓強、溫度、密度都有所增大;在背風面激波則遠離物面,脫體距離增大,同一位置對應的壓強、溫度、密度有所減小。另外,隨著馬赫數(shù)的增大,同一位置對應的壓強、溫度、密度都有所增大,激波脫體距離有所減小。

圖4 壓強分布云圖Fig.4 Contour map of pressure

圖5 溫度分布云圖Fig.5 Contour map of temperature

圖6 密度分布云圖Fig.6 Contour map of density

3 結(jié)束語

本文的分析表明,利用數(shù)值計算的方法模擬光學頭罩高超聲速繞流流場是有效的,且在高馬赫數(shù)、大迎角情況下采用k-ε湍流模型計算得到的結(jié)果比采用SA湍流模型計算得到的結(jié)果更為準確。通過分析不同飛行速度、不同迎角情況下流場壓強、溫度、密度等參數(shù)的分布情況可以看出,頭罩頂部駐點位置的壓強、溫度、密度最大。由此可知，此處的氣動加熱效應和熱輻射效應最為嚴重,對光傳輸及目標成像的影響最為明顯,需重點研究氣動光學效應的校正方法。此外,在大迎角飛行時,迎風面的壓強、溫度、密度較大,也需考慮氣動光學效應的影響。

總的來說,數(shù)值模擬方法與風洞實驗相比更具有靈活性和經(jīng)濟性,針對導彈不同的飛行參數(shù),只需改變計算的初始條件、邊界條件即可得到相應的數(shù)值解算結(jié)果,從而為定量研究球錐形頭罩的氣動特性和氣動光學效應提供數(shù)據(jù)支持。

[1] 殷興良.氣動光學原理[M].北京:中國宇航出版社,2003:55-62.

[2] Tilmann C P,Buter T A,Bowersox R D W. Characterization of the flowfield near a wrap-around fin at mach number 2.8 [R].AIAA-97-0522,1997.

[3] 陳澄,費錦東.側(cè)窗頭罩高速層流流場光學傳輸效應數(shù)值模擬[J].紅外與激光工程,2005,34(5):548-552.

[4] 劉景源,李椿萱.基于Reynolds平均的高超聲速二方程湍流模型[J].航空學報,2008,29(1):28-33.

[5] 于子文,曹義華.前飛旋翼三維湍流場的數(shù)值模擬[J].北京航空航天大學學報,2006,32(7):751-755.

[6] 周大高,柳陽威.改進SA模型對翼型分離流動的數(shù)值模擬[J].北京航空航天大學學報,2012,38(10):1384-1388.

[7] 段卓毅,王小震,陳迎春.三維舵面繞流的N-S方程計算研究[J].飛行力學,2007,25(3):67-70.

[8] Jones W P,Launder B E.The prediction of laminarization with a two-equation model of turbulence[J].Journal of Heat and Mass Transfer,1972,15(2):301-314.

[9] Jones W P,Launder B E.The calculation of low-reynolds number phenomena with a two-equation model of turbulence[J].Journal of Heat and Mass Transfer,1973,16(6):1119-1130.

[10] Spalart P R,Allmaras S R.A one-equation turbulence transport model for aerodynamic flows[R].AIAA-92-0439,1992.

[11] Paciorri R,Dieudonne W,Degrez G,et al.Validation of the spalart-allmaras turbulence model for application in hypersonic flows[R].AIAA-97-2323,1997.

[12] 李素循.典型外形高超聲速流動特性[M].北京:國防工業(yè)出版社,2007:80-109.

Numericalsimulationofhypersonicflowfieldsaroundsphere-cone-shapeddome

XU Dong1, LI Xin1, XU Li2, WANG Ya-hui2

(1.School of Instrument Science and Opto-electronics Engineering, BUAA, Beijing 100191, China; 2.National Key Laboratory of Science and Technology on Aerospace Intelligence Control, Beijing Aerospace Automatic Control Institute, Beijing 100854, China)

The numerical simulation of hypersonic flow fields around sphere-cone-shaped dome under different flight conditions was studied. Firstly, by separately applyingk-εturbulence model and SA (Spalart-Allmaras) turbulence model to solve the N-S equations, the flow fields around sphere-cone-shaped optical dome were simulated. Secondly, the simulation results were compared with related experimental data, and the comparison shows that, under the condition of high Mach number and attack-angle more accurate calculation results can be obtained when usingk-εturbulence model. Lastly, the distribution of pressure, temperature and density of the flow fields around sphere-cone-shaped dome was analyzed, and the variation of the three parameters with flight condition was also presented.

sphere-cone-shaped dome; flow fields; numerical simulation; turbulence model

V211.3

A

1002-0853(2013)06-0491-05

2013-03-29;

2013-09-01; < class="emphasis_bold">網(wǎng)絡出版時間

時間:2013-10-22 14:15

國家自然科學基金資助(61378077)

許東(1973-)，男，江蘇徐州人，副教授，博士，主要研究方向為氣動光學效應及目標探測識別；

李欣(1988-)，女，山東濟南人，碩士研究生，主要研究方向為氣動光學及流場仿真。

(編輯：姚妙慧)