二維彈道修正彈靜態(tài)側(cè)向氣動(dòng)特性仿真研究

鐘 陽(yáng)，王良明，安亮亮

(南京理工大學(xué)能源與動(dòng)力工程學(xué)院，江蘇 南京 210094)

1 引言

美國(guó)陸軍裝備研發(fā)中心于上世紀(jì)90年代提出了低成本強(qiáng)力彈藥計(jì)劃，實(shí)現(xiàn)對(duì)現(xiàn)有庫(kù)存無(wú)控旋成彈丸低成本引信改造。該計(jì)劃經(jīng)歷了GPS 校射引信、一維修正引信以及二維修正引信的發(fā)展歷程。二維修正引信，具有縱向橫向修正能力、精度高、成本低和對(duì)常規(guī)庫(kù)存彈藥兼容性好等優(yōu)點(diǎn)，受到各國(guó)青睞，成為當(dāng)前的研究熱點(diǎn)。目前在二維彈道修正組件研制方面，美國(guó)的精確制導(dǎo)組件(Precision Guidance Kit，PGK)、英國(guó)與法國(guó)合資共同研制的歐洲修正引信(European Correction Fuze，ECF)和法德圣路易斯研究所的方向修正引信(Course Correction Fuze，CCF)居于前列。二維彈道修正彈是加裝了二維修正引信的旋轉(zhuǎn)穩(wěn)定彈丸，具有射程和方向二維修正能力，大大提高了射擊精度，圓概率誤差(Circular Error Probable，CEP)可達(dá)20至50米。由于成本低、精度高和易改裝等優(yōu)點(diǎn)，受到各國(guó)青睞。目前我國(guó)對(duì)二維彈道修正彈還處于探索階段，開(kāi)展了氣動(dòng)、控制、彈道、穩(wěn)定性、結(jié)構(gòu)等大量研究課題。

固定舵二維修正彈的修正組件包含一對(duì)反向差動(dòng)舵和一對(duì)同向控制舵，4片舵都與修正組件固連，控制舵提供控制力和力矩，差動(dòng)舵使控制組件轉(zhuǎn)動(dòng)，調(diào)整控制力和力矩的方向，改變彈丸姿態(tài)，實(shí)現(xiàn)彈道修正。旋成彈丸在加裝了控制組件后，其氣動(dòng)特性發(fā)生了很大的變化，原來(lái)的氣動(dòng)模型已經(jīng)不再適用[1]。然而，氣動(dòng)特性是彈道研究、穩(wěn)定性分析及控制系統(tǒng)設(shè)計(jì)的關(guān)鍵技術(shù)[2]。因此，針對(duì)二維修正彈的氣動(dòng)模型研究逐漸成為熱點(diǎn)。為了更好地建立氣動(dòng)模型，有必要對(duì)一些新的氣動(dòng)現(xiàn)象加以研究[3]，探究其機(jī)理。側(cè)向效應(yīng)對(duì)旋轉(zhuǎn)穩(wěn)定彈丸的穩(wěn)定性有很大影響，基于此，本文以二維彈道修正彈為對(duì)象，采用流場(chǎng)數(shù)值模擬來(lái)探究其靜態(tài)側(cè)向氣動(dòng)特性。為建立完善的二維彈道修正彈氣動(dòng)模型提供參考。

2 數(shù)值計(jì)算方法

通常獲取氣彈箭氣動(dòng)力和力矩的方法有風(fēng)洞實(shí)驗(yàn)、飛行試驗(yàn)和流場(chǎng)數(shù)值計(jì)算。其中，前兩者的成本較高，前期準(zhǔn)備時(shí)間周期長(zhǎng)。相較之下，數(shù)值仿真成本低，耗時(shí)相對(duì)較少，獲得數(shù)據(jù)信息全面，在工程中得到廣泛應(yīng)用。

2.1 控制方程和湍流模型

工程應(yīng)用中通常采用Reynolds-Averaged Navier-Stokes(RANS )方程對(duì)流場(chǎng)進(jìn)行求解。其積分形式為

1)連續(xù)性方程

(1)

式中，Ω為控制體體積，S為控制體表面面積，V=(u，v，w)T為流動(dòng)速度矢量。

2)動(dòng)量守恒方程

(2)

式中，?為張量積，τij為粘性應(yīng)力張量，形式如下

(3)

3)能量守恒方程

(4)

式中，氣體單位質(zhì)量總能

(5)

熱流通量q=(qx，qy，qz)T各分量表達(dá)式為

(6)

式中，k=1，2，3，x1=x，x2=y，x3=z，對(duì)于空氣γ=1.4，層流普朗特?cái)?shù)Pr 取0.72，湍流普朗特?cái)?shù)PrT取0.9，溫度T通過(guò)氣體狀態(tài)方程計(jì)算。層流粘性系數(shù)μ由Sutherland公式給出，湍流粘性系數(shù)μT由湍流模型給出。

2.2 空間離散和時(shí)間離散

二維修正彈在飛行過(guò)程中一般會(huì)經(jīng)歷亞音速、跨音速和超音速三種狀態(tài)。簡(jiǎn)單低耗散迎風(fēng)矢通量分裂格式(Simple Low-Dissipation Advection Upstream Splitting Method，SLAU)[4]具有耗散低，格式簡(jiǎn)單，在低馬赫數(shù)范圍內(nèi)無(wú)需調(diào)參數(shù)等優(yōu)點(diǎn)，其改進(jìn)型SLAU2在壓力通量中加入了與來(lái)流馬赫數(shù)相關(guān)的耗散項(xiàng)，繼承了原格式優(yōu)點(diǎn)，同時(shí)實(shí)現(xiàn)了比其它全速格式更簡(jiǎn)潔的形式?；诖?，本文采用SLAU2為空間離散格式。

時(shí)間離散采用經(jīng)典的三階顯示Runge-Kutta對(duì)流場(chǎng)進(jìn)行推進(jìn)。

2.3 邊界條件

本文仿真所用到的邊界條件有遠(yuǎn)場(chǎng)邊界條件和壁面邊界條件。超音速進(jìn)口遠(yuǎn)場(chǎng)直接采用來(lái)流條件；超音速出口條件采用內(nèi)場(chǎng)插值獲得；采用黎曼不變量來(lái)構(gòu)造亞音速進(jìn)口和出口遠(yuǎn)場(chǎng)邊界條件。彈丸表面為無(wú)滑移壁面條件。具體形式參見(jiàn)文獻(xiàn)[5]。

2.4 湍流模型

Spalart-Allmaras (SA)湍流模型具有計(jì)算量小、穩(wěn)定性好、計(jì)算精度滿足工程要求等優(yōu)點(diǎn)，在葉輪機(jī)計(jì)算、航空航天等領(lǐng)域得到了廣泛應(yīng)用。SA模型版本眾多，本文研究對(duì)象的運(yùn)動(dòng)速度能夠達(dá)到超音速，可以適當(dāng)選擇一種可壓SA模型。本文采用了標(biāo)準(zhǔn)可壓SA-noft2[6]湍流模型，其積分形式為

(7)

湍流粘性系數(shù)計(jì)算如下

(8)

(9)

式中，d為壁面距離，fν2=1-χ/(1+χfν1)，S表示渦量大小，形式如下

(10)

模型(7)中，輔助關(guān)系式

cw1=cb1/κ2+(1+cb2)/σ

(11)

其余封閉系數(shù)：cb1=0.1355，cb2=0.622，cv1=7.1，σ=2/3，cw2=0.3，cw3=2，κ=0.41 。

3 算例驗(yàn)證

采用M910子母彈來(lái)驗(yàn)證本文數(shù)值計(jì)算方法獲取靜態(tài)氣動(dòng)特性的準(zhǔn)確性。計(jì)算模型和計(jì)算域網(wǎng)格見(jiàn)圖1。網(wǎng)格劃分技術(shù)參數(shù)見(jiàn)文獻(xiàn)[7]。

圖1 M910計(jì)算域軸向截面網(wǎng)格

對(duì)攻角為3°，馬赫數(shù)分別為0.4、0.5、0.7、0.9、1.0、1.2、1.5、2.0、2.5、3.0、3.5、4.0等工況進(jìn)行流場(chǎng)數(shù)值計(jì)算。

仿真結(jié)果如圖2-4所示。阻力系數(shù)在跨音速附近誤差較大，最大誤差出現(xiàn)在馬赫數(shù)為1.0處，約13%，其余馬赫數(shù)范圍內(nèi)符合較好；升力系數(shù)導(dǎo)數(shù)基本都在實(shí)驗(yàn)值散點(diǎn)內(nèi)部，超音速結(jié)果比平均值稍偏大，誤差在10%以內(nèi)；俯仰力矩系數(shù)導(dǎo)數(shù)與實(shí)驗(yàn)值符合度較高，超音速結(jié)果較實(shí)驗(yàn)值略小，整體誤差不超過(guò)6%。從總體仿真結(jié)果來(lái)看，最大誤差處是跨音速附近的阻力系數(shù)。造成的原因可能是跨音速流場(chǎng)具有雙曲和橢圓兩種特性，對(duì)數(shù)值計(jì)算造成較大困難，文獻(xiàn)[7]采用數(shù)值方法計(jì)算出的阻力系數(shù)和本文存在相同問(wèn)題。本文方法總體仿真結(jié)果誤差在15%以內(nèi)，滿足工程要求。

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

圖3 升力系數(shù)導(dǎo)數(shù)隨馬赫數(shù)變化曲線

圖4 靜力矩系數(shù)導(dǎo)數(shù)隨馬赫數(shù)變化曲線

4 二維修正彈側(cè)向氣動(dòng)特性仿真結(jié)果分析

以某二維修正彈為例，通過(guò)數(shù)值方法來(lái)模擬其外部空氣繞流特性，獲取并分析彈丸表面壓力系數(shù)分布和全彈靜態(tài)側(cè)向空氣動(dòng)力系數(shù)。來(lái)流條件為：馬赫數(shù)2.0；壓力101325Pa；溫度288.15K；密度1.225kg/m3；攻角α為0°、±2°、±4°、5°、±6°、8°、10°。

4.1 幾何模型與網(wǎng)格劃分

幾何模型見(jiàn)圖5。仿真時(shí)保持控制組件為“+”形布局。差動(dòng)舵位于鉛垂面內(nèi)，鴨舵1繞Oy軸正向偏轉(zhuǎn)4°，鴨舵3繞Oy軸負(fù)向偏轉(zhuǎn)4°。差動(dòng)舵產(chǎn)生繞Ox軸正向的導(dǎo)轉(zhuǎn)力矩；控制舵處于水平位置，鴨舵2和鴨舵4均繞Oz軸負(fù)向偏轉(zhuǎn)8°，控制舵舵面產(chǎn)生向上的力。

圖5 彈丸外形和計(jì)算坐標(biāo)系Oxyz

采用六面體網(wǎng)格對(duì)計(jì)算域網(wǎng)格進(jìn)行劃分；遠(yuǎn)場(chǎng)邊界位置、近壁網(wǎng)格第1層厚度和延展比參考文獻(xiàn)[7]；網(wǎng)格細(xì)節(jié)見(jiàn)圖6，網(wǎng)格數(shù)約為360萬(wàn)。彈長(zhǎng)978mm，質(zhì)心距離彈頂632.8mm處。

圖6 二維修正彈計(jì)算域三軸截面網(wǎng)格

4.2 仿真結(jié)果與分析

數(shù)值仿真的流動(dòng)情況和彈丸表面壓力分布如圖7所示。

圖7 攻角8°時(shí)彈丸表面壓力及差動(dòng)舵附近流動(dòng)情況

圖7左圖為鴨舵1，流線位置在y=0.034m附近；右圖為鴨舵3，流線位置在y=-0.034m附近。從圖7中可以看出，在有攻角的情況下，彈丸上下表面壓力不對(duì)稱，迎風(fēng)區(qū)壓力明顯高于背風(fēng)區(qū)壓力，特別是在彈丸頭部附近。

由于差動(dòng)舵存在舵偏角，導(dǎo)致了氣流流向發(fā)生了變化，鴨舵1附近流動(dòng)向Oz軸負(fù)向偏轉(zhuǎn)，而鴨舵3附近流動(dòng)向Oz軸正向偏轉(zhuǎn)，都對(duì)其后的彈丸表面壓力分布產(chǎn)生了影響，造成了左右不對(duì)稱。為了便于定量分析，對(duì)不同彈軸位置的周向壓力系數(shù)分布進(jìn)行觀察。定義彈丸表面網(wǎng)格點(diǎn)坐標(biāo)矢量在Oyz上的投影與Oy軸之間的夾角為φ，繞Ox軸負(fù)向轉(zhuǎn)過(guò)的角度為正，如鴨舵1位置φ=0°，鴨舵2位置φ=90°，鴨舵3位置φ=180°，鴨舵4位置φ=270°。修正組件表面網(wǎng)格點(diǎn)橫坐標(biāo)x∈[0，175)mm，圓錐段表面網(wǎng)格點(diǎn)橫坐標(biāo)x∈[175，584)mm，圓柱段表面網(wǎng)格點(diǎn)橫坐標(biāo)x∈[584，864)mm。不同x處彈丸表面周向壓力系數(shù)Cp隨φ變化關(guān)系曲線如圖8所示。

圖8 不同軸向位置周向壓力系數(shù)分布

從圖8可以看出，總體波形大致中間高兩邊低，即迎風(fēng)區(qū)的壓力高，背風(fēng)區(qū)的壓力低。由x=0.20m、x=0.25m和x=0.30m處壓力系數(shù)分布曲線看出，彈丸圓錐段迎風(fēng)區(qū)受到鴨舵3的干擾，壓力分布形成了明顯不對(duì)稱，壓力峰值在φ=172°附近，即從彈丸正下方向Oz軸負(fù)向偏移。Cp最高可達(dá)0.33，且壓力衰減很快。

由x∈[0.35，0.50]m范圍內(nèi)的壓力系數(shù)分布曲線可以看出，隨著流動(dòng)逐漸靠近圓柱段，壓力分布的不對(duì)稱性逐漸消失，說(shuō)明差動(dòng)舵對(duì)流場(chǎng)的不對(duì)稱干擾逐漸減弱，且壓力衰減與之前相比明顯減弱。

由x=0.6m和x=0.7m處壓力系數(shù)分布曲線可以看出，壓力呈現(xiàn)階梯式下降又有所回復(fù)，這是由于處于圓柱段膨脹區(qū)造成的。迎風(fēng)區(qū)壓力分布的不對(duì)稱性幾乎消失。背風(fēng)區(qū)在圓柱段開(kāi)始呈現(xiàn)壓力不對(duì)稱，但隨著氣體流動(dòng)到x=0.7m處時(shí)，壓力不對(duì)稱消失。

一般認(rèn)為二維彈道修正彈的差動(dòng)舵處于攻角平面內(nèi)時(shí)，不產(chǎn)生靜態(tài)側(cè)向氣動(dòng)力和力矩。然而，由上述仿真分析知，差動(dòng)舵會(huì)使得氣流發(fā)生偏轉(zhuǎn)，相當(dāng)于在后體的迎風(fēng)區(qū)和背風(fēng)區(qū)形成方向相反的附加側(cè)滑角。由于迎風(fēng)區(qū)壓力比背風(fēng)區(qū)大，后體附加側(cè)滑角導(dǎo)致的氣動(dòng)力效應(yīng)更加顯著，側(cè)向力和力矩由迎風(fēng)區(qū)的附加側(cè)滑角主導(dǎo)。全彈的側(cè)向力系數(shù)和側(cè)向力矩系數(shù)數(shù)值計(jì)算結(jié)果如圖9所示。

圖9 +形布局不同攻角下側(cè)向力和力矩系數(shù)

從圖9側(cè)向力系數(shù)cz曲線可以看出，當(dāng)攻角為0°時(shí)，幾乎不產(chǎn)生側(cè)向力；當(dāng)攻角不為0°時(shí)，產(chǎn)生了側(cè)向力，且正攻角產(chǎn)生Oz軸正向側(cè)向力，負(fù)攻角產(chǎn)生沿Oz軸負(fù)向側(cè)向力。從對(duì)圖8的分析中可知，當(dāng)攻角為8°時(shí)，壓力分布的峰值向Oz軸負(fù)向偏移(如圖8所示)，從彈頭看向彈尾，左側(cè)的壓力高于右側(cè)的壓力，最終形成向Oz軸正向的側(cè)向力。彈丸表面壓力分布的規(guī)律和側(cè)向力系數(shù)計(jì)算結(jié)果一致。從曲線可知，在仿真條件內(nèi)側(cè)向力系數(shù)與攻角有近似線性關(guān)系。

從圖9中靜態(tài)側(cè)向力矩系my曲線可以看出，在本文仿真條件內(nèi)，其隨攻角增大而增大，有近似線性關(guān)系。從仿真結(jié)果知，當(dāng)α>0°時(shí)，my>0，又因本文研究的對(duì)象是靜不穩(wěn)定彈丸，因此當(dāng)α>0°時(shí)流動(dòng)綜合效果必然是向Oz軸正向拉偏，而此時(shí)背風(fēng)區(qū)的鴨舵1將流動(dòng)向Oz軸負(fù)向拉偏，迎風(fēng)區(qū)處鴨舵3將流動(dòng)向Oz軸正向拉偏，再次說(shuō)明了靜態(tài)側(cè)向效應(yīng)由迎風(fēng)區(qū)鴨舵造成的偏流主導(dǎo)。

在旋轉(zhuǎn)穩(wěn)定彈的彈道計(jì)算中，動(dòng)態(tài)側(cè)向效應(yīng)，即馬格努斯效應(yīng)是不可忽略的。根據(jù)文獻(xiàn)[7]，當(dāng)Ma=2.0且攻角α=3°時(shí)，馬格努斯力矩系數(shù)導(dǎo)數(shù)Cnpα≈0.2，來(lái)流速度大小v∞=680.4m/s，轉(zhuǎn)速ωx=7153rad/s，彈徑dM910=0.0162m。通過(guò)文獻(xiàn)中馬格努斯力矩系數(shù)導(dǎo)數(shù)公式Cnpα=my/(ωxdM910/2v∞)sinα反算出動(dòng)態(tài)側(cè)向氣動(dòng)力矩系數(shù)my=8.9186×10-4，與圖9Ma=2.0，α=3°時(shí)，靜態(tài)側(cè)向氣動(dòng)力矩系數(shù)值2×10-4具有同等量級(jí)。由此可見(jiàn)，在二維彈道修正彈氣動(dòng)建模時(shí)應(yīng)合理加入靜態(tài)側(cè)向氣動(dòng)特性。

5 結(jié)論

本文對(duì)差動(dòng)舵處于攻角面時(shí)的二維修正彈進(jìn)行了數(shù)值仿真研究，重點(diǎn)研究分析了修正組件對(duì)彈丸側(cè)向氣動(dòng)特性的影響，并得出以下結(jié)論：

1)差動(dòng)舵能夠使得周圍氣流向自身方向進(jìn)行偏轉(zhuǎn)，在有攻角的情況下，加之迎風(fēng)區(qū)和背風(fēng)區(qū)的影響，彈丸上下差動(dòng)舵附近相當(dāng)于附加了方向相反大小不同的側(cè)滑角，形成側(cè)向壓力分布不對(duì)稱。

2)控制組件后局部區(qū)域壓力分布不對(duì)稱非常明顯，周向壓力系數(shù)分布峰值的偏向與迎風(fēng)區(qū)的鴨舵舵偏角有關(guān)。從彈頭看向彈尾，如果鴨舵將氣流向右拉偏，則峰值出現(xiàn)在左側(cè)，類似于產(chǎn)生一個(gè)附加側(cè)滑角。

3)在本文條件下，彈丸側(cè)向力和力矩系數(shù)隨攻角近似呈線性關(guān)系。側(cè)向力方向與上述第2條結(jié)論有關(guān)，如果壓力峰值出現(xiàn)在Oz軸負(fù)向，則側(cè)向力指向Oz軸正向，cz>0，否則反之。側(cè)向力矩系數(shù)量級(jí)與馬格努斯力矩相當(dāng)，在氣動(dòng)建模時(shí)應(yīng)予考慮。