于 雷,滕 鵬,魯 藝,周志剛
(空軍工程大學(xué)工程學(xué)院,西安 710038)
在現(xiàn)代生產(chǎn)設(shè)計(jì)和科學(xué)研究中,實(shí)驗(yàn)研究與數(shù)值模擬是兩個(gè)重要的手段,二者相輔相成。尤其是在航空航天及武器系統(tǒng)的設(shè)計(jì)研究中,存在著許多可以通過(guò)數(shù)值模擬來(lái)解決的問(wèn)題,諸如飛行器復(fù)雜流場(chǎng)及氣動(dòng)力、飛行器與外掛分離氣動(dòng)干擾等[1]。通過(guò)對(duì)這些問(wèn)題的數(shù)值模擬,可以更加深刻地理解問(wèn)題產(chǎn)生的機(jī)理,為實(shí)驗(yàn)研究和結(jié)果整理提供指導(dǎo),節(jié)省實(shí)驗(yàn)所需的大量的人力、物力和時(shí)間。Fluent是一款基于有限體積法,集網(wǎng)格生成、流場(chǎng)求解、數(shù)據(jù)可視化等于一體的國(guó)內(nèi)外使用最廣泛的商用計(jì)算流體動(dòng)力學(xué)(CFD)分析軟件系統(tǒng);Gambit作為Fluent的前處理軟件,可以讀入多種CAD軟件的三維幾何模型和多種CAE軟件的網(wǎng)格模型[2]。
某型航空子母彈是一種用于對(duì)付地面集群目標(biāo)、封鎖機(jī)場(chǎng)跑道、電站電線的對(duì)地攻擊武器,其母彈采用的是某型電視制導(dǎo)炸彈。本文以該電視制導(dǎo)炸彈(以下簡(jiǎn)稱炸彈)的外形為基礎(chǔ),應(yīng)用CFD的相關(guān)理論、方法以及軟件,對(duì)炸彈在亞、跨、超音速以及不同攻角條件下的三維流場(chǎng)進(jìn)行數(shù)值模擬,對(duì)不同工況下的三維流場(chǎng)結(jié)構(gòu)特性以及炸彈氣動(dòng)特性進(jìn)行分析研究。
利用Gambit 2.2軟件對(duì)該炸彈以及繞流流場(chǎng)進(jìn)行三維建模,設(shè)置邊界類型并進(jìn)行網(wǎng)格劃分[3]。炸彈表面網(wǎng)格如圖1所示,整個(gè)流場(chǎng)計(jì)算區(qū)域網(wǎng)格劃分如圖2所示。
圖1 母彈外形尺寸及表面網(wǎng)格Fig.1 The three-dimensional shape and surface grid of TV-guided bomb
圖2 三維繞流流場(chǎng)及網(wǎng)格劃分Fig.2 Three-dimensional grid models of the flow field
劃分網(wǎng)格時(shí)采用混合網(wǎng)格技術(shù)[2,4-5],即在彈頭部及流動(dòng)規(guī)律復(fù)雜的彈體周圍區(qū)域采用半結(jié)構(gòu)網(wǎng)格并且網(wǎng)格布置得較密集,而在遠(yuǎn)離彈體流動(dòng)較穩(wěn)定的區(qū)域采用非結(jié)構(gòu)網(wǎng)格并且網(wǎng)格布置得較稀疏,這樣生成速度較快,生成質(zhì)量也比較高。最終劃分網(wǎng)格總數(shù)為267997,沒(méi)有負(fù)網(wǎng)格,而且正交性、長(zhǎng)寬比、擴(kuò)展比都達(dá)到一定要求。根據(jù)流體力學(xué)的相似性理論[6],可以通過(guò)改變來(lái)流的方向,使之與炸彈成相應(yīng)的夾角,便可實(shí)現(xiàn)炸彈具有一定攻角飛行的繞流流場(chǎng)數(shù)值模擬,并不需要改變炸彈姿態(tài)與模型,以及整個(gè)流場(chǎng)的網(wǎng)格結(jié)構(gòu)。
積分形式的 Navier- Stokes方程為[2,4-5]
式中:矢量H為源項(xiàng);ρ為密度;u、υ、w為時(shí)均速度;E為單位質(zhì)量的總能;p為流體壓力;τ為粘性應(yīng)力張量;q為熱流通量;A為表面面積矢量。
單位質(zhì)量總能E與總焓H的關(guān)系為
湍流模型為標(biāo)準(zhǔn) k-ε 兩方程模型[2,4-5],湍流動(dòng)能k及其耗散率ε從下列輸運(yùn)方程得到:
方程中Gk代表由平均速度梯度引起的湍動(dòng)能生成項(xiàng);YM是可壓縮性影響項(xiàng);C1ε和 C2ε是常量;σk和 σε分別是k和ε的湍流普朗特?cái)?shù)。
方程差分求解的定解條件以及初始條件如下[7-9]:1)來(lái)流邊界條件,對(duì)于超音速入口,給定無(wú)窮遠(yuǎn)來(lái)流條件;對(duì)于亞音速入口,給定來(lái)流的總壓和總溫邊界條件;2)下游邊界條件,對(duì)于超音速出口,邊界條件采用外推;對(duì)于亞音速出口,給定出口壓力;3)邊界條件,物面邊界采用絕熱壁假設(shè)和無(wú)滑移條件,并且速度采用無(wú)穿透條件;外邊界取壓力遠(yuǎn)場(chǎng)邊界條件;4)初始條件,整個(gè)初始流場(chǎng)采用無(wú)限遠(yuǎn)來(lái)流。
采用耦合求解法,在求解連續(xù)、動(dòng)量和能量方程的基礎(chǔ)上分別求解其他的標(biāo)量方程。應(yīng)用Gauss-Seidel迭代法求解由有限體積法離散基本控制方程組得到的代數(shù)方程組[2,4],為了便于收斂和滿足計(jì)算精度對(duì)控制方程各量采用二階差分格式,對(duì)湍動(dòng)能及其耗散率采用一階差分格式。本文采用Fluent 6.2軟件求解控制方程,完成對(duì)炸彈三維繞流流場(chǎng)的數(shù)值模擬,其中物性參數(shù)及湍流模型參數(shù)采用默認(rèn)值[3]。
對(duì)炸彈三維繞流流場(chǎng)進(jìn)行數(shù)值模擬,得到來(lái)流馬赫數(shù)Ma=0.7~1.4,攻角α=0°~15°時(shí)的流場(chǎng)靜壓力分布圖、溫度分布云圖、馬赫數(shù)分布圖分別如圖3~圖11所示。其中計(jì)算來(lái)流符合標(biāo)準(zhǔn)大氣條件,氣溫T=288.15 K。
由圖3~圖5可知,當(dāng)來(lái)流馬赫數(shù)Ma=0.8,在炸彈頭部形成一個(gè)高壓區(qū),在炸彈頭部與圓柱部連接處氣流發(fā)生膨脹,該處的馬赫數(shù)亦相應(yīng)增大;在彈體表面由于空氣的粘性而形成摩擦阻力,其能量損失轉(zhuǎn)化為內(nèi)能從而使得彈體以及附面層的溫度升高;在炸彈底部形成一個(gè)低壓區(qū),由于周圍空氣向底部補(bǔ)充時(shí)流動(dòng)非常紊亂,產(chǎn)生許多小的漩渦形成底部渦流,氣流速度的大小與方向在此區(qū)域均發(fā)生較大變化,大部分能量轉(zhuǎn)化為渦流的動(dòng)能和內(nèi)能使得尾部溫度升高尤為顯著;流場(chǎng)內(nèi)無(wú)激波產(chǎn)生。進(jìn)一步的數(shù)值模擬表明,在亞音速范圍內(nèi)炸彈繞流流場(chǎng)有相似的結(jié)構(gòu)形式。
圖3 靜壓分布圖(Ma=0.8,α=0°)Fig.3 The pressure distribution(Ma=0.8,α =0°)
圖4 溫度分布圖(Ma=0.8,α=0°)Fig.4 The temperature distribution(Ma=0.8,α =0°)
圖5 馬赫數(shù)分布圖(Ma=0.8,α=0°)Fig.5 The mach number distribution(Ma=0.8,α =0°)
分析圖6~圖11可知,在跨音速條件下(來(lái)流馬赫數(shù)Ma=1.0、1.1),當(dāng)前方空氣流到激波面上時(shí),其壓強(qiáng)突然增大,必然使得被壓縮的空氣溫度突然升高、內(nèi)能增大,激波阻力所消耗的動(dòng)能轉(zhuǎn)化為空氣內(nèi)能;并且由于表面摩擦和底部渦流的原因使得彈體周圍和炸彈底部的溫度升高,尤其是后者溫度升高更為顯著。在馬赫數(shù)較小時(shí)(Ma=1.0、1.1),炸彈頭部產(chǎn)生脫體正激波,隨著馬赫數(shù)的增大(Ma=1.2),激波角逐漸減小;在炸彈底亦有尾部激波存在,并且隨著馬赫數(shù)的增大(Ma>1.2),激波傾斜角亦逐漸減小。同時(shí),當(dāng)來(lái)流馬赫數(shù)Ma>1.2,在炸彈頭部與圓柱部連接處、炸彈尾部有兩道膨脹波產(chǎn)生(如圖9~圖11所示),這是由于超音速氣流流過(guò)上述連接處后發(fā)生膨脹、加速所致,使得兩處氣流的密度減小,壓強(qiáng)下降,溫度下降,流速增大;并且隨著馬赫數(shù)的增大,膨脹波角度亦逐漸減小。如圖6~圖8所示,當(dāng)來(lái)流與炸彈之間具有一定攻角時(shí),所產(chǎn)生的激波與膨脹波發(fā)生相應(yīng)的傾斜,并且其強(qiáng)度也會(huì)在背風(fēng)面與迎風(fēng)面產(chǎn)生不同的差異。
圖6 靜壓分布圖(Ma=1.0,α=15°)Fig.6 The pressure distribution(Ma=1.0,α =15°)
圖7 溫度分布圖(Ma=1.0,α=15°)Fig.7 The temperature distribution(Ma=1.0,α =15°)
圖8 馬赫數(shù)分布圖(Ma=1.0,α=15°)Fig.8 The mach number distribution(Ma=1.0,α =15°)
圖9 靜壓分布圖(Ma=1.2,α=0°)Fig.9 The pressure distribution(Ma=1.2,α =0°)
圖10 溫度分布圖(Ma=1.2,α=0°)Fig.10 The temperature distribution(Ma=1.2,α =0°)
圖11 馬赫數(shù)分布圖(Ma=1.2,α=0°)Fig.11 The mach number distribution(Ma=1.2,α =0°)
通過(guò)大量三維繞流流場(chǎng)的數(shù)值模擬,得到炸彈在來(lái)流馬赫數(shù)為0.7~1.4范圍內(nèi)(炸彈的有效工作段)的阻力系數(shù)與升力系數(shù)。其中,不同攻角條件下的阻力系數(shù)Cx與馬赫數(shù)的關(guān)系如圖12所示??梢钥闯?,當(dāng)攻角 α 為5°、10°或者15°時(shí),升阻力系數(shù) Cx變化趨勢(shì)與零升阻力系數(shù)Cxo基本相同——來(lái)流馬赫數(shù)Ma在0.7~1.1范圍內(nèi),阻力系數(shù)隨來(lái)流馬赫數(shù)的增大而迅速增大;當(dāng)來(lái)流Ma在1.1~1.4范圍內(nèi),阻力系數(shù)基本保持不變。不同攻角條件下的升力系數(shù)Cy與馬赫數(shù)的關(guān)系如圖13所示??梢钥闯?,升力系數(shù)Cy有如下變化趨勢(shì):在亞、跨音速階段(Ma=0.7~1.1),隨著馬赫數(shù)的增加,升力系數(shù)Cy基本保持緩慢增加;Ma在1.1~1.4范圍內(nèi),升力系數(shù)Cy略有減小。
圖12 阻力系數(shù)曲線Fig.12 The drag coefficient curve
圖13 升力系數(shù)曲線Fig.13 The lift coefficient curve
本文應(yīng)用計(jì)算流體力學(xué)軟件Fluent 6.2進(jìn)行數(shù)值模擬所得的仿真結(jié)果很好地展現(xiàn)了某型電視制導(dǎo)炸彈三維繞流流場(chǎng)的流動(dòng)狀態(tài)效果,得出了炸彈繞流流場(chǎng)宏觀結(jié)構(gòu)的變化規(guī)律和炸彈空氣動(dòng)力特性參數(shù),本文給出的流場(chǎng)波譜為進(jìn)一步等離子體激勵(lì)改善該電視制導(dǎo)炸彈氣動(dòng)特性研究提供一定的參考與借鑒。
[1]朱自強(qiáng),吳子牛,李津.應(yīng)用計(jì)算流體力學(xué)[M].北京:北京航空航天大學(xué)出版社,2001.
[2]王福軍.計(jì)算流體動(dòng)力學(xué)分析[M].北京:清華大學(xué)出版社,2004.
[3]韓占忠,王敬,蘭小平.流體工程仿真計(jì)算實(shí)例與應(yīng)用[M].北京:北京理工大學(xué)出版社,2004.
[4]ANDERSON J D Jr.Computational fluid dynamics:the basics with application[M].Beijing:Tsinghua University Press,2002.
[5]張涵信,沈孟育.計(jì)算流體力學(xué)—差分方法的原理和應(yīng)用[M].北京:國(guó)防工業(yè)出版社,2003.
[6]惲起麟.風(fēng)洞實(shí)驗(yàn)[M].北京:國(guó)防工業(yè)出版社,2000.
[7]陳大偉,馬小亮,楊國(guó)偉.超聲速大攻角旋成體迎風(fēng)激波數(shù)值模擬[J].力學(xué)學(xué)報(bào),2006(6):721-732.
[8]姜毅,郝繼光,傅德彬.導(dǎo)彈發(fā)射過(guò)程三維非定常數(shù)值模擬[J]. 兵工學(xué)報(bào),2008(8):911-915.
[9]雷娟棉,吳甲生.機(jī)載布撒器繞流場(chǎng)數(shù)值模擬[J].彈道學(xué)報(bào),2004(4):33-37.