鐘 陽,王良明,安亮亮
(南京理工大學能源與動力工程學院,江蘇 南京 210094)
美國陸軍裝備研發(fā)中心于上世紀90年代提出了低成本強力彈藥計劃,實現(xiàn)對現(xiàn)有庫存無控旋成彈丸低成本引信改造。該計劃經(jīng)歷了GPS 校射引信、一維修正引信以及二維修正引信的發(fā)展歷程。二維修正引信,具有縱向橫向修正能力、精度高、成本低和對常規(guī)庫存彈藥兼容性好等優(yōu)點,受到各國青睞,成為當前的研究熱點。目前在二維彈道修正組件研制方面,美國的精確制導組件(Precision Guidance Kit,PGK)、英國與法國合資共同研制的歐洲修正引信(European Correction Fuze,ECF)和法德圣路易斯研究所的方向修正引信(Course Correction Fuze,CCF)居于前列。二維彈道修正彈是加裝了二維修正引信的旋轉(zhuǎn)穩(wěn)定彈丸,具有射程和方向二維修正能力,大大提高了射擊精度,圓概率誤差(Circular Error Probable,CEP)可達20至50米。由于成本低、精度高和易改裝等優(yōu)點,受到各國青睞。目前我國對二維彈道修正彈還處于探索階段,開展了氣動、控制、彈道、穩(wěn)定性、結(jié)構(gòu)等大量研究課題。
固定舵二維修正彈的修正組件包含一對反向差動舵和一對同向控制舵,4片舵都與修正組件固連,控制舵提供控制力和力矩,差動舵使控制組件轉(zhuǎn)動,調(diào)整控制力和力矩的方向,改變彈丸姿態(tài),實現(xiàn)彈道修正。旋成彈丸在加裝了控制組件后,其氣動特性發(fā)生了很大的變化,原來的氣動模型已經(jīng)不再適用[1]。然而,氣動特性是彈道研究、穩(wěn)定性分析及控制系統(tǒng)設計的關(guān)鍵技術(shù)[2]。因此,針對二維修正彈的氣動模型研究逐漸成為熱點。為了更好地建立氣動模型,有必要對一些新的氣動現(xiàn)象加以研究[3],探究其機理。側(cè)向效應對旋轉(zhuǎn)穩(wěn)定彈丸的穩(wěn)定性有很大影響,基于此,本文以二維彈道修正彈為對象,采用流場數(shù)值模擬來探究其靜態(tài)側(cè)向氣動特性。為建立完善的二維彈道修正彈氣動模型提供參考。
通常獲取氣彈箭氣動力和力矩的方法有風洞實驗、飛行試驗和流場數(shù)值計算。其中,前兩者的成本較高,前期準備時間周期長。相較之下,數(shù)值仿真成本低,耗時相對較少,獲得數(shù)據(jù)信息全面,在工程中得到廣泛應用。
工程應用中通常采用Reynolds-Averaged Navier-Stokes(RANS )方程對流場進行求解。其積分形式為
1)連續(xù)性方程
(1)
式中,Ω為控制體體積,S為控制體表面面積,V=(u,v,w)T為流動速度矢量。
2)動量守恒方程
(2)
式中,?為張量積,τij為粘性應力張量,形式如下
(3)
3)能量守恒方程
(4)
式中,氣體單位質(zhì)量總能
(5)
熱流通量q=(qx,qy,qz)T各分量表達式為
(6)
式中,k=1,2,3,x1=x,x2=y,x3=z,對于空氣γ=1.4,層流普朗特數(shù)Pr 取0.72,湍流普朗特數(shù)PrT取0.9,溫度T通過氣體狀態(tài)方程計算。層流粘性系數(shù)μ由Sutherland公式給出,湍流粘性系數(shù)μT由湍流模型給出。
二維修正彈在飛行過程中一般會經(jīng)歷亞音速、跨音速和超音速三種狀態(tài)。簡單低耗散迎風矢通量分裂格式(Simple Low-Dissipation Advection Upstream Splitting Method,SLAU)[4]具有耗散低,格式簡單,在低馬赫數(shù)范圍內(nèi)無需調(diào)參數(shù)等優(yōu)點,其改進型SLAU2在壓力通量中加入了與來流馬赫數(shù)相關(guān)的耗散項,繼承了原格式優(yōu)點,同時實現(xiàn)了比其它全速格式更簡潔的形式?;诖?,本文采用SLAU2為空間離散格式。
時間離散采用經(jīng)典的三階顯示Runge-Kutta對流場進行推進。
本文仿真所用到的邊界條件有遠場邊界條件和壁面邊界條件。超音速進口遠場直接采用來流條件;超音速出口條件采用內(nèi)場插值獲得;采用黎曼不變量來構(gòu)造亞音速進口和出口遠場邊界條件。彈丸表面為無滑移壁面條件。具體形式參見文獻[5]。
Spalart-Allmaras (SA)湍流模型具有計算量小、穩(wěn)定性好、計算精度滿足工程要求等優(yōu)點,在葉輪機計算、航空航天等領(lǐng)域得到了廣泛應用。SA模型版本眾多,本文研究對象的運動速度能夠達到超音速,可以適當選擇一種可壓SA模型。本文采用了標準可壓SA-noft2[6]湍流模型,其積分形式為
(7)
湍流粘性系數(shù)計算如下
(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 。
采用M910子母彈來驗證本文數(shù)值計算方法獲取靜態(tài)氣動特性的準確性。計算模型和計算域網(wǎng)格見圖1。網(wǎng)格劃分技術(shù)參數(shù)見文獻[7]。
圖1 M910計算域軸向截面網(wǎng)格
對攻角為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等工況進行流場數(shù)值計算。
仿真結(jié)果如圖2-4所示。阻力系數(shù)在跨音速附近誤差較大,最大誤差出現(xiàn)在馬赫數(shù)為1.0處,約13%,其余馬赫數(shù)范圍內(nèi)符合較好;升力系數(shù)導數(shù)基本都在實驗值散點內(nèi)部,超音速結(jié)果比平均值稍偏大,誤差在10%以內(nèi);俯仰力矩系數(shù)導數(shù)與實驗值符合度較高,超音速結(jié)果較實驗值略小,整體誤差不超過6%。從總體仿真結(jié)果來看,最大誤差處是跨音速附近的阻力系數(shù)。造成的原因可能是跨音速流場具有雙曲和橢圓兩種特性,對數(shù)值計算造成較大困難,文獻[7]采用數(shù)值方法計算出的阻力系數(shù)和本文存在相同問題。本文方法總體仿真結(jié)果誤差在15%以內(nèi),滿足工程要求。
圖2 阻力系數(shù)隨馬赫數(shù)變化曲線
圖3 升力系數(shù)導數(shù)隨馬赫數(shù)變化曲線
圖4 靜力矩系數(shù)導數(shù)隨馬赫數(shù)變化曲線
以某二維修正彈為例,通過數(shù)值方法來模擬其外部空氣繞流特性,獲取并分析彈丸表面壓力系數(shù)分布和全彈靜態(tài)側(cè)向空氣動力系數(shù)。來流條件為:馬赫數(shù)2.0;壓力101325Pa;溫度288.15K;密度1.225kg/m3;攻角α為0°、±2°、±4°、5°、±6°、8°、10°。
幾何模型見圖5。仿真時保持控制組件為“+”形布局。差動舵位于鉛垂面內(nèi),鴨舵1繞Oy軸正向偏轉(zhuǎn)4°,鴨舵3繞Oy軸負向偏轉(zhuǎn)4°。差動舵產(chǎn)生繞Ox軸正向的導轉(zhuǎn)力矩;控制舵處于水平位置,鴨舵2和鴨舵4均繞Oz軸負向偏轉(zhuǎn)8°,控制舵舵面產(chǎn)生向上的力。
圖5 彈丸外形和計算坐標系Oxyz
采用六面體網(wǎng)格對計算域網(wǎng)格進行劃分;遠場邊界位置、近壁網(wǎng)格第1層厚度和延展比參考文獻[7];網(wǎng)格細節(jié)見圖6,網(wǎng)格數(shù)約為360萬。彈長978mm,質(zhì)心距離彈頂632.8mm處。
圖6 二維修正彈計算域三軸截面網(wǎng)格
數(shù)值仿真的流動情況和彈丸表面壓力分布如圖7所示。
圖7 攻角8°時彈丸表面壓力及差動舵附近流動情況
圖7左圖為鴨舵1,流線位置在y=0.034m附近;右圖為鴨舵3,流線位置在y=-0.034m附近。從圖7中可以看出,在有攻角的情況下,彈丸上下表面壓力不對稱,迎風區(qū)壓力明顯高于背風區(qū)壓力,特別是在彈丸頭部附近。
由于差動舵存在舵偏角,導致了氣流流向發(fā)生了變化,鴨舵1附近流動向Oz軸負向偏轉(zhuǎn),而鴨舵3附近流動向Oz軸正向偏轉(zhuǎn),都對其后的彈丸表面壓力分布產(chǎn)生了影響,造成了左右不對稱。為了便于定量分析,對不同彈軸位置的周向壓力系數(shù)分布進行觀察。定義彈丸表面網(wǎng)格點坐標矢量在Oyz上的投影與Oy軸之間的夾角為φ,繞Ox軸負向轉(zhuǎn)過的角度為正,如鴨舵1位置φ=0°,鴨舵2位置φ=90°,鴨舵3位置φ=180°,鴨舵4位置φ=270°。修正組件表面網(wǎng)格點橫坐標x∈[0,175)mm,圓錐段表面網(wǎng)格點橫坐標x∈[175,584)mm,圓柱段表面網(wǎng)格點橫坐標x∈[584,864)mm。不同x處彈丸表面周向壓力系數(shù)Cp隨φ變化關(guān)系曲線如圖8所示。
圖8 不同軸向位置周向壓力系數(shù)分布
從圖8可以看出,總體波形大致中間高兩邊低,即迎風區(qū)的壓力高,背風區(qū)的壓力低。由x=0.20m、x=0.25m和x=0.30m處壓力系數(shù)分布曲線看出,彈丸圓錐段迎風區(qū)受到鴨舵3的干擾,壓力分布形成了明顯不對稱,壓力峰值在φ=172°附近,即從彈丸正下方向Oz軸負向偏移。Cp最高可達0.33,且壓力衰減很快。
由x∈[0.35,0.50]m范圍內(nèi)的壓力系數(shù)分布曲線可以看出,隨著流動逐漸靠近圓柱段,壓力分布的不對稱性逐漸消失,說明差動舵對流場的不對稱干擾逐漸減弱,且壓力衰減與之前相比明顯減弱。
由x=0.6m和x=0.7m處壓力系數(shù)分布曲線可以看出,壓力呈現(xiàn)階梯式下降又有所回復,這是由于處于圓柱段膨脹區(qū)造成的。迎風區(qū)壓力分布的不對稱性幾乎消失。背風區(qū)在圓柱段開始呈現(xiàn)壓力不對稱,但隨著氣體流動到x=0.7m處時,壓力不對稱消失。
一般認為二維彈道修正彈的差動舵處于攻角平面內(nèi)時,不產(chǎn)生靜態(tài)側(cè)向氣動力和力矩。然而,由上述仿真分析知,差動舵會使得氣流發(fā)生偏轉(zhuǎn),相當于在后體的迎風區(qū)和背風區(qū)形成方向相反的附加側(cè)滑角。由于迎風區(qū)壓力比背風區(qū)大,后體附加側(cè)滑角導致的氣動力效應更加顯著,側(cè)向力和力矩由迎風區(qū)的附加側(cè)滑角主導。全彈的側(cè)向力系數(shù)和側(cè)向力矩系數(shù)數(shù)值計算結(jié)果如圖9所示。
圖9 +形布局不同攻角下側(cè)向力和力矩系數(shù)
從圖9側(cè)向力系數(shù)cz曲線可以看出,當攻角為0°時,幾乎不產(chǎn)生側(cè)向力;當攻角不為0°時,產(chǎn)生了側(cè)向力,且正攻角產(chǎn)生Oz軸正向側(cè)向力,負攻角產(chǎn)生沿Oz軸負向側(cè)向力。從對圖8的分析中可知,當攻角為8°時,壓力分布的峰值向Oz軸負向偏移(如圖8所示),從彈頭看向彈尾,左側(cè)的壓力高于右側(cè)的壓力,最終形成向Oz軸正向的側(cè)向力。彈丸表面壓力分布的規(guī)律和側(cè)向力系數(shù)計算結(jié)果一致。從曲線可知,在仿真條件內(nèi)側(cè)向力系數(shù)與攻角有近似線性關(guān)系。
從圖9中靜態(tài)側(cè)向力矩系my曲線可以看出,在本文仿真條件內(nèi),其隨攻角增大而增大,有近似線性關(guān)系。從仿真結(jié)果知,當α>0°時,my>0,又因本文研究的對象是靜不穩(wěn)定彈丸,因此當α>0°時流動綜合效果必然是向Oz軸正向拉偏,而此時背風區(qū)的鴨舵1將流動向Oz軸負向拉偏,迎風區(qū)處鴨舵3將流動向Oz軸正向拉偏,再次說明了靜態(tài)側(cè)向效應由迎風區(qū)鴨舵造成的偏流主導。
在旋轉(zhuǎn)穩(wěn)定彈的彈道計算中,動態(tài)側(cè)向效應,即馬格努斯效應是不可忽略的。根據(jù)文獻[7],當Ma=2.0且攻角α=3°時,馬格努斯力矩系數(shù)導數(shù)Cnpα≈0.2,來流速度大小v∞=680.4m/s,轉(zhuǎn)速ωx=7153rad/s,彈徑dM910=0.0162m。通過文獻中馬格努斯力矩系數(shù)導數(shù)公式Cnpα=my/(ωxdM910/2v∞)sinα反算出動態(tài)側(cè)向氣動力矩系數(shù)my=8.9186×10-4,與圖9Ma=2.0,α=3°時,靜態(tài)側(cè)向氣動力矩系數(shù)值2×10-4具有同等量級。由此可見,在二維彈道修正彈氣動建模時應合理加入靜態(tài)側(cè)向氣動特性。
本文對差動舵處于攻角面時的二維修正彈進行了數(shù)值仿真研究,重點研究分析了修正組件對彈丸側(cè)向氣動特性的影響,并得出以下結(jié)論:
1)差動舵能夠使得周圍氣流向自身方向進行偏轉(zhuǎn),在有攻角的情況下,加之迎風區(qū)和背風區(qū)的影響,彈丸上下差動舵附近相當于附加了方向相反大小不同的側(cè)滑角,形成側(cè)向壓力分布不對稱。
2)控制組件后局部區(qū)域壓力分布不對稱非常明顯,周向壓力系數(shù)分布峰值的偏向與迎風區(qū)的鴨舵舵偏角有關(guān)。從彈頭看向彈尾,如果鴨舵將氣流向右拉偏,則峰值出現(xiàn)在左側(cè),類似于產(chǎn)生一個附加側(cè)滑角。
3)在本文條件下,彈丸側(cè)向力和力矩系數(shù)隨攻角近似呈線性關(guān)系。側(cè)向力方向與上述第2條結(jié)論有關(guān),如果壓力峰值出現(xiàn)在Oz軸負向,則側(cè)向力指向Oz軸正向,cz>0,否則反之。側(cè)向力矩系數(shù)量級與馬格努斯力矩相當,在氣動建模時應予考慮。