完顏振海，馮順山，董永香，李順平

（北京理工大學(xué)爆炸科學(xué)與技術(shù)國家重點實驗室，北京 100081）

0 引言

當(dāng)處于運(yùn)動狀態(tài)的子母戰(zhàn)斗部拋撒子彈藥時，受拋撒力、氣動力、彈道諸元等諸多因素的影響，子彈藥會產(chǎn)生翻轉(zhuǎn)運(yùn)動［1］。翻轉(zhuǎn)運(yùn)動會嚴(yán)重影響子彈藥的彈道性能，進(jìn)而影響子彈藥的可靠性和散布合理性，從而影響終點效應(yīng)。

子彈藥在翻轉(zhuǎn)時，其攻角范圍為－180°～180°，因此需要對彈體的大攻角時氣動特性進(jìn)行研究。大攻角時彈身繞流具有顯著非線性特性，理論分析其氣動特性比較困難，因此一般采用風(fēng)洞試驗和數(shù)值計算方法。由于攻角范圍大，風(fēng)洞試驗需要變換支撐方式，試驗周期長、成本高。目前，采用求解N -S方程模擬細(xì)長體大攻角時氣動特性已經(jīng)得到了廣泛的應(yīng)用［2－4］，計算結(jié)果也與風(fēng)洞實驗比較吻合，但現(xiàn)有數(shù)值計算的攻角范圍一般小于90°，且一般只給出了特定攻角時的氣動特性。文中采用激波捕捉型和激波裝配型兩種網(wǎng)格對尖拱圓柱形彈進(jìn)行數(shù)值分析，由于彈體對稱性取攻角范圍為0°～180°。

1 數(shù)值計算方法

1.1 控制方程

積分形式的控制方程，其質(zhì)量、動量和能量方程分別如式（1）、式（2）、式（3）所示［5］：

其中：Ω、S 分別是控制體的體積和表面積；ρ、v、p、fe分別為氣體的密度、速度矢量、壓強(qiáng)和體積力矢量；E、H、T分別是單位體積氣體的總能、總焓和總溫；qH表示外部熱源為粘性摩擦應(yīng)力張量。

1.2 湍流模型

采用雷諾平均運(yùn)輸方程進(jìn)行數(shù)值計算，由于引入了雷諾應(yīng)力項，需要額外的湍流模型封閉方程。文中采用Shih等［6］提出的可實現(xiàn)k-ε兩方程模型，其在附著和分離流動方面有較好的結(jié)果。

式中：Gk表示由于平均速度梯度產(chǎn)生的湍流動能；k、ε、μt分別表示湍流動能、湍流耗散率和渦粘性系數(shù)；σk、σε分別為k和ε的湍流普朗特數(shù)；Cε2為常數(shù)，大小為1.92。

2 數(shù)值計算模型

2.1 模型幾何尺寸

文中研究子彈藥的尖拱段為正切型，即圓弧段與圓柱母線相切，其外形尺寸如圖1所示。

圖1 模型外形尺寸圖

2.2 邊界條件和計算網(wǎng)格

采用兩種不同的計算網(wǎng)格，一種是激波捕捉型網(wǎng)格，另一種是激波裝配型網(wǎng)格。

激波捕捉型網(wǎng)格采用的是遠(yuǎn)場邊界條件，用來獲取子彈藥在各攻角范圍內(nèi)的氣動力系數(shù)，其網(wǎng)格如圖2（a）所示，采用O O型混合網(wǎng)格，除了極軸附近為三棱柱網(wǎng)格外，其余為六面體網(wǎng)格，網(wǎng)格數(shù)為434700，其中環(huán)向69，徑向70，軸向90。

激波裝配型網(wǎng)格是在捕捉型網(wǎng)格計算結(jié)果的基礎(chǔ)上建立起來的，可以更準(zhǔn)確地反映激波形態(tài)，其網(wǎng)格如圖2（b）所示，采用超音速速度入口和出口邊界條件。由于子彈藥較短和計算資源限制，僅建立半模型，未考慮大攻角時非對稱性渦脫落問題。網(wǎng)格總數(shù)為151290，其中環(huán)向41，徑向41，軸向90。

物面均采用非滑移壁面條件，對于子彈藥來說，

圖2 計算網(wǎng)格圖

其受到的氣動力和力矩源于其表面的壓力分布和剪切應(yīng)力分布［5］，根據(jù)粘性流動的非滑移壁面條件，可以通過積分得到：

其中，R、L、D分別表示子彈受到的合力、升力和阻力；Sb表示子彈藥表面積。

3 結(jié)果與分析

數(shù)值計算結(jié)果采用 Missile Datco m（以下簡稱Datcom）進(jìn)行對比分析，Datcom是美國空軍研究實驗室在理論方法研究和幾十年積累的試驗數(shù)據(jù)的基礎(chǔ)上開發(fā)的，其有效攻角范圍0°～90°，并在α≤45°時與風(fēng)洞試驗結(jié)果符合很好［7－9］。

在超音速范圍內(nèi)，Datco m對單獨彈身的法向力N和壓心位置的計算采用的是Van Dyke混合理論或二階激波膨脹波理論，對零攻角時的軸向力A的計算采用的是經(jīng)驗公式和Van Dyke混合理論或二階激波膨脹波理論，對于非零攻角其采用的是Allen＆Per kins的改進(jìn)橫向流動理論［10］。

3.1 阻力特性

子彈藥阻力系數(shù)隨馬赫數(shù)的變化曲線如圖3所示，計算結(jié)果符合尖拱圓柱體的阻力系數(shù)分布規(guī)律。從圖中可以看出，超音速范圍內(nèi)，阻力系數(shù)先隨著馬赫數(shù)的增大而增大，在Ma＝1.3時達(dá)到峰值，隨后逐漸減小。

子彈藥在Ma＝1.5時阻力系數(shù)隨攻角變化曲線如圖4所示，數(shù)值計算和Datco m在α＝0°～45°范圍內(nèi)相差不到5%；α＞45°時，數(shù)值計算值略高，最大相差約為12%。從圖中可以看出，子彈藥的阻力系數(shù)隨攻角變化呈現(xiàn)出近似余弦曲線規(guī)律，阻力系數(shù)在攻角為85°時達(dá)到最大值，其大小約為零升阻力系數(shù)的13倍；α＝180°時（相當(dāng)于鈍頭尖尾），阻力系數(shù)大約是零攻角（即尖頭鈍尾）時的2.7倍。

圖3 阻力系數(shù)隨馬赫數(shù)變化曲線

圖4 Ma＝1.5時阻力系數(shù)隨攻角變化曲線

3.2 升力特性

子彈藥在Ma＝1.5時升力系數(shù)隨攻角的變化曲線如圖5所示。升力系數(shù)隨攻角變化呈現(xiàn)出近似正弦曲線規(guī)律。升力系數(shù)在α＝50°時達(dá)到最大值，α＝130°達(dá)到最小值，零值則發(fā)生在兩端以及α＝85°附近。

圖5 Ma＝1.5時升力系數(shù)隨攻角變化曲線

3.3 壓心位置

數(shù)值計算和Datco m計算的壓心位置對比曲線如圖6所示，隨著攻角的增大，壓心位置逐漸后移，并在α＝180°達(dá)到最大值。

圖6 Ma＝1.5時壓心位置隨攻角變化曲線

3.4 流場特性

流場采用數(shù)值紋影法顯示，其方法是給出密度梯度大小分布圖，即‖?ρ‖［11］。數(shù)值紋影法不僅能顯示間斷變化的密度場（如激波），還可以給出連續(xù)變化的密度場（如尾跡、渦等）。

圖7給出了子彈藥縱向截面的流場形態(tài)隨攻角變化圖，零攻角時彈體頭部形成斜激波，尖拱段末端和圓柱段后臺階處產(chǎn)生膨脹波，彈體的對稱性還導(dǎo)致尾部激波。隨著攻角的增大，子彈藥迎流面不斷擴(kuò)大，頭部激波的形態(tài)由圓錐形的激波變?yōu)榉菍ΨQ錐形激波，激波形式由斜激波逐漸變?yōu)楣渭げǎ瑫r彈體背流面發(fā)生渦脫落，如圖8所示。

圖7 縱向面流場數(shù)值紋影圖

圖8 攻角為30°時來流方向流場數(shù)值紋影圖

4 結(jié)論

從以上分析可以得出以下結(jié)論：1）子彈藥超音速氣動特性隨攻角的改變變化很大，攻角從0°增加到180°的過程中：阻力系數(shù)變化呈現(xiàn)出近似余弦曲線規(guī)律；升力系數(shù)的變化呈現(xiàn)正弦曲線規(guī)律；壓心位置則逐漸后移。2）隨著攻角的增大，尖拱形彈的頭部激波的形態(tài)由圓錐形的斜激波逐漸變?yōu)榉菍ΨQ錐形弓形激波，同時彈體背流面發(fā)生渦脫落。

［1］ 楊啟仁.子母彈飛行動力學(xué)［M］.北京：國防工業(yè)出版社，1999.

［2］ Hsieh T，Priolo F J，War dlaw A B.Calculations and comparisons of the flowfield about an ogive cylinder at M＝3.5［J］.Jour nal of Spacecraft and Rockets，1993，30（6）：665－673.

［3］ Josyula E.Co mputational Simulation i mprovements of supersonic high angle－of－attack missile flo ws［J］.Jour nal of Spacecraft and Rockets，1999，36（1）：59－66.

［4］ 鄧有奇，周乃春，胡漢東.尖拱彈身大攻角N S方程數(shù)值計算［J］.空氣動力學(xué)學(xué)報，2000（z1）：109－114.

［5］ Charles Hirsch.Numerical computation of internal and exter nal flows［M］.Oxfor d：Elsevier’s，2007.

［6］ Shih T H，Liou W W，Yang Zhigang，et al.A new k－ε eddy viscosity model f or high Reynolds number tur bulent flows［J］.Co mputers Fluids，1995，24（3）：227－238.

［7］ Blake W B.Missile datcom：1997 status and future plans，AIAA－1997－2280［R］.1997.

［8］ Simon J M，Blake W B.Missile datcom：high angle of attack capabilities，AIAA－1999－4258［R］.1998.

［9］ Abney E J，Mc Daniel M A.High angle of attack aerodynamics predictions using missile datcom，AIAA－2005－5086［R］.2005.

［10］ Blake W B.Missile datco m：apllications to projetiles，AIAA－1989－3370［R］.1989.

［11］ James DeSpiritio，Peter Plostins.CFD prediction of MP910 pr ojectile aerodynamics：unsteady wake effect on Magnus moment，AIAA－2007－6580［R］.2007.