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

(1.北京航空航天大學(xué) 儀器科學(xué)與光電工程學(xué)院, 北京 100191; 2.北京航天自動(dòng)控制研究所 宇航與智能控制技術(shù)國(guó)家級(jí)重點(diǎn)實(shí)驗(yàn)室, 北京 100854)

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

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

(1.北京航空航天大學(xué) 儀器科學(xué)與光電工程學(xué)院, 北京 100191; 2.北京航天自動(dòng)控制研究所 宇航與智能控制技術(shù)國(guó)家級(jí)重點(diǎn)實(shí)驗(yàn)室, 北京 100854)

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

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

0 引言

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

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

1 計(jì)算方法

1.1 控制方程

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

(1)

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

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

p=ρRT

(2)

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

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

(3)

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

μ=C1T3/2/(T+C2)

(4)

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

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

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

(5)

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

1.2 湍流模型

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

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

(6)

(7)

(8)

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

(9)

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

(10)

(11)

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

1.3 網(wǎng)格劃分

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

1.4 邊界條件

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

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

2.1 湍流模型有效性分析

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

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

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

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

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

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

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

圖4 壓強(qiáng)分布云圖Fig.4 Contour map of pressure

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

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

3 結(jié)束語(yǔ)

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

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

[1] 殷興良.氣動(dòng)光學(xué)原理[M].北京:中國(guó)宇航出版社,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] 陳澄,費(fèi)錦東.側(cè)窗頭罩高速層流流場(chǎng)光學(xué)傳輸效應(yīng)數(shù)值模擬[J].紅外與激光工程,2005,34(5):548-552.

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

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

[6] 周大高,柳陽(yáng)威.改進(jìn)SA模型對(duì)翼型分離流動(dòng)的數(shù)值模擬[J].北京航空航天大學(xué)學(xué)報(bào),2012,38(10):1384-1388.

[7] 段卓毅,王小震,陳迎春.三維舵面繞流的N-S方程計(jì)算研究[J].飛行力學(xué),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] 李素循.典型外形高超聲速流動(dòng)特性[M].北京:國(guó)防工業(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)絡(luò)出版時(shí)間

時(shí)間:2013-10-22 14:15

國(guó)家自然科學(xué)基金資助(61378077)

許東(1973-)，男，江蘇徐州人，副教授，博士，主要研究方向?yàn)闅鈩?dòng)光學(xué)效應(yīng)及目標(biāo)探測(cè)識(shí)別；

李欣(1988-)，女，山東濟(jì)南人，碩士研究生，主要研究方向?yàn)闅鈩?dòng)光學(xué)及流場(chǎng)仿真。

(編輯：姚妙慧)