遠(yuǎn)程尾控制導(dǎo)炮彈氣動(dòng)特性及仿真分析

王桂奇,劉福朝,劉 寧,李存健

(1.北京信息科技大學(xué) 自動(dòng)化學(xué)院, 北京 100192; 2.北京信息科技大學(xué) 高動(dòng)態(tài)導(dǎo)航技術(shù)北京市重點(diǎn)實(shí)驗(yàn)室, 北京 100192)

0 引言

制導(dǎo)炮彈具有遠(yuǎn)程穿透、高速飛行、精準(zhǔn)打擊的優(yōu)點(diǎn),在現(xiàn)代化戰(zhàn)爭(zhēng)中扮演著至關(guān)重要的角色[1-3]。遠(yuǎn)程制導(dǎo)炮彈飛行過(guò)程中轉(zhuǎn)速高,控制難度大,如何降低轉(zhuǎn)速,提高測(cè)量精度、高效控制是當(dāng)前研究的熱門(mén)問(wèn)題[4]。

北京理工大學(xué)的陳慶森等[5]利用多剛體動(dòng)力學(xué)的方法研究了封控彈彈傘系統(tǒng)末彈道的動(dòng)態(tài)特性,使封控彈末彈道滿足設(shè)計(jì)要求,減旋葉片能使封控彈減旋到理想值。南京理工大學(xué)的逯麒等[6]提出了一種適用于電磁發(fā)射超高速制導(dǎo)炮彈的滾轉(zhuǎn)運(yùn)動(dòng)穩(wěn)定區(qū)域和減旋策略,主動(dòng)起旋與減旋止擺控制策略為電磁發(fā)射超高速制導(dǎo)炮彈飛行控制提供新穎可行的方法。王中原等[7]針對(duì)阻力環(huán)疊加減旋翼片的二維彈道修正技術(shù)方案,進(jìn)行橫向彈道修正能力及飛行穩(wěn)定性分析,并通過(guò)炮射試驗(yàn)驗(yàn)證該方法的準(zhǔn)確性和有效性。邱荷等[8]根據(jù)典型末敏彈的結(jié)構(gòu)特點(diǎn)及減速減旋段的運(yùn)動(dòng)特性,提出了一種利用Matlab和VR技術(shù)對(duì)其進(jìn)行可視化仿真研究的方法。中北大學(xué)的劉一鳴[9]、祝敬德[10]等設(shè)計(jì)了一種適用于彈載環(huán)境的減旋控制系統(tǒng),采用半捷聯(lián)模擬控制減旋電路,有效提高了測(cè)量精度。

本文中描述一種利用船尾裝置進(jìn)行彈體減旋的制導(dǎo)炮彈外形,這種減旋船尾采用軸承與船尾和前端非硬連接的裝配方式,船尾前端的軸承連接裝置,能夠保證彈體出膛后達(dá)到減旋效果。從而降低彈體的轉(zhuǎn)速,達(dá)到可控的目的。

1 物理模型

首先對(duì)具有減旋船尾的制導(dǎo)炮彈氣動(dòng)外形進(jìn)行設(shè)計(jì),采用三維建模軟件進(jìn)行結(jié)構(gòu)件設(shè)計(jì)及裝配,然后將建好的炮彈模型導(dǎo)入GAMBIT中。本文中研究的制導(dǎo)炮彈由彈頭、戰(zhàn)斗部、發(fā)動(dòng)機(jī)、密封圈、船尾、閉氣蓋等部分組成,制導(dǎo)炮彈整體氣動(dòng)外形結(jié)構(gòu)如圖1所示。

2 數(shù)值方法

2.1 工程計(jì)算方法

在空氣動(dòng)力學(xué)理論的基礎(chǔ)上進(jìn)行工程計(jì)算,并對(duì)炮彈結(jié)構(gòu)和流場(chǎng)進(jìn)行適當(dāng)?shù)募僭O(shè)和簡(jiǎn)化,最后進(jìn)行求解運(yùn)算。外彈道研究方法必須滿足計(jì)算速度快,計(jì)算簡(jiǎn)捷,計(jì)算精度高3個(gè)特點(diǎn)。所以,工程計(jì)算方法被廣泛應(yīng)用在外彈道研究中。相對(duì)于外形而言,工程計(jì)算方法可以使阻力、升力、翻轉(zhuǎn)力矩等系數(shù)的計(jì)算結(jié)果和正常值的誤差控制在10%以內(nèi),甚至更小。所以,工程計(jì)算方法比較適用于彈體結(jié)構(gòu)變化對(duì)氣動(dòng)特性的影響,并且通過(guò)計(jì)算可以獲得炮彈的氣動(dòng)參數(shù)等信息。

圖1 炮彈三維模型

根據(jù)文獻(xiàn)[11]中所推導(dǎo),炮彈組件的阻力系數(shù)表達(dá)式為:

CxBW=1.1(CxOBW+CxiBW)

(1)

升力系數(shù)為:

(2)

2.2 控制方程

連續(xù)性方程和動(dòng)量方程[11]為:

(3)

(4)

通常使用Boussinesq假設(shè)來(lái)求解能量方程,一般認(rèn)為雷諾應(yīng)力和平均速度梯度成正相關(guān),其表達(dá)式為[12]:

(5)

(6)

2.3 湍流模型

湍流作為流體的一種流動(dòng)形式,通常是無(wú)序的、多角度、不規(guī)則的流動(dòng)。一般而言,湍流具有良好的擴(kuò)散特性。在空氣動(dòng)力學(xué)中,湍流一般是指某特定區(qū)域、較小時(shí)間內(nèi)的風(fēng)速波動(dòng)。湍流的物理模型由不同大小的漩渦組成,并且漩渦大小、流向等特征都是無(wú)序的,無(wú)法用數(shù)學(xué)公式表示出來(lái)。

在空氣動(dòng)力學(xué)的研究中Spalart-Allmaras模型、k-ε模型、k-ω模型和雷諾應(yīng)力模型經(jīng)常使用,本文中采用Spalart-Allmaras模型。Spalart-Allmaras模型是用于航空領(lǐng)域的較為簡(jiǎn)便的單方程模型,用于解決墻壁束縛(wall-bounded)流動(dòng)方面的問(wèn)題,在解決逆壓梯度的邊界層問(wèn)題展現(xiàn)出很好的效果。在空氣動(dòng)力學(xué)中,一般使用該模型解決飛行器等繞流流場(chǎng)計(jì)算問(wèn)題。

(7)

2.4 網(wǎng)格劃分及邊界條件

為了精確獲得阻力系數(shù)、升力系數(shù)和力矩系數(shù)等的計(jì)算值,在GAMBIT劃分網(wǎng)格時(shí),需要在彈體附近布置較密的網(wǎng)格。針對(duì)彈體計(jì)算域進(jìn)行建模時(shí),也需要格外注意,一般選取計(jì)算域?yàn)榇髨A柱體和小圓柱體2個(gè)。對(duì)大圓柱體而言,長(zhǎng)度為彈體長(zhǎng)度的7.5倍,直徑為彈體直徑的40倍;小圓柱體即為內(nèi)部加密區(qū),長(zhǎng)度為彈體的2.5倍,直徑為彈體直徑的5倍,如圖2所示。在網(wǎng)格劃分時(shí)越靠近彈體的網(wǎng)格越密,遠(yuǎn)離彈體的網(wǎng)格可以相對(duì)稀疏一些,其中彈頭部分的網(wǎng)格越密,計(jì)算結(jié)果越準(zhǔn)確。為了加快計(jì)算機(jī)運(yùn)算速度,一般特地將計(jì)算域分割成兩個(gè)部分,只對(duì)一半進(jìn)行仿真計(jì)算。

圖2 計(jì)算域和網(wǎng)格劃分示意圖

模型網(wǎng)格劃分模型如圖3所示,整個(gè)計(jì)算域采用了網(wǎng)格加密的辦法,使得計(jì)算更加精確。整個(gè)計(jì)算域網(wǎng)格數(shù)在100萬(wàn)以上。將劃分好網(wǎng)格的炮彈流場(chǎng)域?qū)隖luent中進(jìn)行計(jì)算。

圖3 彈體周圍網(wǎng)格劃分示意圖

設(shè)置邊界條件的時(shí)候選用密度基耦合顯式求解器進(jìn)行求解,對(duì)來(lái)流和物面分別選用遠(yuǎn)場(chǎng)邊界條件和黏性邊界條件:

(8)

其中:ui=(i=1,2,3)為坐標(biāo)方向的速度分量;F為壓力;ρ為密度;n為法線方向。

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

使用Fluent軟件對(duì)減旋彈體進(jìn)行氣動(dòng)仿真分析時(shí)需要計(jì)算不同攻角和不同馬赫下的工況值,分別選取馬赫數(shù)為1Ma,1.5Ma,2Ma,2.5Ma,3Ma,3.5Ma,4Ma;攻角α取1°,2°,3°,4°;共28種狀態(tài),對(duì)應(yīng)升力系數(shù)、阻力系數(shù)、力矩系數(shù)、X方向速度分布云圖、全彈表面壓力分布云圖等計(jì)算仿真分析如下。

3.1 動(dòng)態(tài)監(jiān)控窗口分析

以4°攻角,2.5Ma為例,經(jīng)過(guò)Fluent 計(jì)算,得到以下阻力系數(shù),升力系數(shù)和力矩系數(shù)等。通過(guò)觀察可知,經(jīng)過(guò)3 000次迭代之后發(fā)現(xiàn),各系數(shù)均變化不大,說(shuō)明此時(shí)處于收斂狀態(tài)。

500次迭代完成之后, 4°攻角、1.5Ma時(shí)阻力系數(shù)的收斂值是0.174 0;4°攻角、2.0Ma時(shí)阻力系數(shù)收斂值是0.255 1;4°攻角、2.5Ma時(shí)阻力系數(shù)收斂值是0.324 6; 4°攻角、3.0Ma時(shí)阻力系數(shù)收斂值是0.466 5。通過(guò)3 000次迭代發(fā)現(xiàn),阻力系數(shù)的變化曲線圖在200次迭代,基本呈現(xiàn)直線狀態(tài),即此時(shí)在經(jīng)過(guò)迭代阻力系數(shù)的值基本不發(fā)生變化,說(shuō)明阻力系數(shù)的計(jì)算值達(dá)到收斂狀態(tài)。將4°攻角、1.5Ma,4°攻角、2Ma,4°攻角、2.5Ma,4°攻角、3Ma這4種不同馬赫的阻力系數(shù)值在Excel表格中對(duì)比發(fā)現(xiàn),在4°攻角一定的情況下,隨著馬赫數(shù)的升高,阻力系數(shù)值也在不斷升高。阻力系數(shù)隨迭代過(guò)程變化曲線如圖4所示。

500次迭代完成之后, 4°攻角、1.5Ma時(shí)升力系數(shù)的收斂值是0.066 99;4°攻角、2.0Ma時(shí)升力系數(shù)收斂值是0.136 2;4°攻角、2.5Ma時(shí)升力系數(shù)收斂值是0.230 3;3.0Ma時(shí)升力系數(shù)收斂值是0.408 8。通過(guò)3 000次迭代發(fā)現(xiàn),升力系數(shù)的變化曲線圖在200次迭代,基本呈現(xiàn)直線狀態(tài),即此時(shí)在經(jīng)過(guò)迭代升力系數(shù)的值基本不發(fā)生變化,說(shuō)明升力系數(shù)的計(jì)算值達(dá)到收斂狀態(tài)。將4°攻角、1.5Ma,4°攻角、2Ma,4°攻角、2.5Ma,4°攻角、3Ma這4種不同馬赫的升力系數(shù)值在Excel表格中對(duì)比發(fā)現(xiàn),在4°攻角一定的情況下,隨著馬赫數(shù)的升高,升力系數(shù)值也在不斷升高。升力系數(shù)隨迭代過(guò)程變化曲線如圖5所示。

圖4 阻力系數(shù)隨迭代過(guò)程變化曲線

圖5 升力系數(shù)隨迭代過(guò)程變化曲線

500次迭代完成之后, 4°攻角、1.5Ma時(shí)力矩系數(shù)的收斂值是0.086 68;4°攻角、2.0Ma時(shí)力矩系數(shù)收斂值是0.200 3;4°攻角、2.5Ma時(shí)力矩系數(shù)收斂值是0.355 0;4°攻角、3.0Ma時(shí)力矩系數(shù)收斂值是0.646 0。通過(guò)3 000次迭代發(fā)現(xiàn),力矩系數(shù)的變化曲線圖在200次迭代,基本呈現(xiàn)直線狀態(tài),即此時(shí)在經(jīng)過(guò)迭代力矩系數(shù)的值基本不發(fā)生變化,說(shuō)明力矩系數(shù)的計(jì)算值達(dá)到收斂狀態(tài)。將4°攻角、1.5Ma,4°攻角、2Ma,4°攻角、2.5Ma,4°攻角、3Ma這4種不同馬赫的力矩系數(shù)值在Excel表格中對(duì)比發(fā)現(xiàn),在4°攻角一定的情況下,隨著馬赫數(shù)的升高,力矩系數(shù)值也在不斷升高。力矩系數(shù)隨迭代過(guò)程變化曲線如圖6所示。

圖6 力矩系數(shù)隨迭代過(guò)程變化曲線

3.2 空氣繞炮彈流動(dòng)情況窗口分析

圖7(a)—圖7(d)四幅圖分別為4°攻角下1.5Ma、2Ma、2.5Ma、3Ma時(shí)X方向速度分布云圖。

當(dāng)馬赫數(shù)Ma>1時(shí),即速度V高于音速時(shí),高動(dòng)態(tài)下的彈體同時(shí)受到摩擦作用、渦阻作用以及激波作用。此時(shí),隨著周圍環(huán)境壓強(qiáng)Pa,密度ρ,溫度T的升高,彈體表面的流速會(huì)瞬間下降。在Fluent求解運(yùn)算中一般選擇滿足氣體狀態(tài)方程的理想氣體(ideal-gas),在數(shù)學(xué)理論中滿足氣體狀態(tài)方程的理想氣體的激波作用沒(méi)有厚度,是不連續(xù)的,但在高動(dòng)態(tài)彈體實(shí)際飛行中,還應(yīng)考慮密度和黏性的作用,實(shí)際的熱傳導(dǎo)速率也不容忽視。

雖然激波作用干擾高動(dòng)態(tài)飛行的彈體,但實(shí)際的影響非常小,所以在數(shù)值計(jì)算中,得到的激波值較小,并且隨著馬赫數(shù)的升高,激波值逐漸減小。高馬赫飛行的炮彈會(huì)受到氣流壓縮等影響,計(jì)算機(jī)計(jì)算的參數(shù)也會(huì)發(fā)生突變現(xiàn)象。仔細(xì)觀察圖像發(fā)現(xiàn),炮彈的頭部與尾部均可以看到頭部激波和尾部激波,并且馬赫數(shù)越高,產(chǎn)生的激波效果越明顯。

圖7 X方向速度分布云圖

圖8(a)—圖8(d)四幅圖分別為4°攻角下1.5Ma、2Ma、2.5Ma、3.0Ma時(shí)全彈表面壓力分布云圖。炮彈表面的壓力分布是炮彈設(shè)計(jì)的重要參考指標(biāo)。由圖8可以看出,炮彈在飛行中彈頭所受的壓力最大,1.5Ma時(shí)彈頭所受壓力值為2.27e+0.5 Pa,2Ma時(shí)彈頭所受壓力值為5.55e+0.5 Pa,2.5Ma時(shí)彈頭所受壓力值為7.47e+0.5 Pa,3.0Ma時(shí)彈頭所受壓力值為9.82e+0.5 Pa。船尾所受壓力較小,分別為1.03e+0.5、1.71e+0.5、6.02e+0.4、2.42e+0.3 Pa。

圖8 全彈表面的壓力分布云圖

由此可知,隨著馬赫數(shù)的增大,彈體表面壓力值逐漸增大,其中彈頭所承受的壓力值最大,在材料處理方面,也要對(duì)彈頭的材料進(jìn)行特殊處理,保證彈體飛行過(guò)程中的穩(wěn)定性。

3.3 炮彈模型受力分析

通過(guò)Fluent軟件仿真,在cd-history、cl-history和cm-history 中讀取可以得到最終收斂的阻力系數(shù)、升力系數(shù)和力矩系數(shù)值。為了節(jié)省計(jì)算機(jī)運(yùn)行時(shí)間和降低網(wǎng)格劃分難度,只對(duì)半個(gè)計(jì)算域進(jìn)行計(jì)算,這樣避免了流場(chǎng)中壓力波等值邊界的反射引起數(shù)值計(jì)算的不穩(wěn)定。所以,真實(shí)的阻力系數(shù)(Cd)、升力系數(shù)(Cl)和力矩系數(shù)(Cm)應(yīng)為計(jì)算機(jī)仿真得到計(jì)算值的2倍。

從圖9中可知,攻角一定時(shí),隨著馬赫數(shù)的增大,阻力值(Cd)也在不斷上升,且上升程度由迅速逐漸變?yōu)榫徛?但最終呈上升趨勢(shì);馬赫數(shù)一定時(shí),隨著攻角的增大,阻力值先下降后上升到峰值,然后下降,總體呈上升趨勢(shì)。

圖9 阻力隨馬赫數(shù)和攻角的變化曲線

從圖10中可知,攻角一定時(shí),隨著馬赫數(shù)的增大,升力值(Cl)也在不斷上升,雖然在0°攻角時(shí),上升率很小,但是升力系數(shù)在其余攻角時(shí)隨著馬赫數(shù)的增加而不斷增大,且攻角越大,升力系數(shù)的增長(zhǎng)趨勢(shì)越快。馬赫數(shù)一定時(shí),隨著攻角的增大,升力值也在不斷上升,且馬赫越大,升力系數(shù)的增長(zhǎng)趨勢(shì)越快。

圖10 升力隨馬赫數(shù)和攻角的變化曲線

從圖11中可知,攻角一定時(shí),隨著馬赫數(shù)的增大,力矩系數(shù)值(Cm)也在不斷上升,且隨著攻角的增大,力矩系數(shù)的增量差值也逐漸變大,攻角越大,升力系數(shù)的增長(zhǎng)趨勢(shì)越快。馬赫數(shù)一定時(shí),隨著攻角的增大,力矩系數(shù)逐漸呈上升趨勢(shì),且上升程度逐漸緩慢。

圖11 力矩隨馬赫數(shù)和攻角的變化曲線

4 結(jié)論

為解決高轉(zhuǎn)速導(dǎo)致遠(yuǎn)程制導(dǎo)炮彈控制難度大的問(wèn)題,設(shè)計(jì)了一種利用船尾裝置進(jìn)行彈體減旋的氣動(dòng)外形結(jié)構(gòu),對(duì)氣動(dòng)參數(shù)進(jìn)行了仿真和分析,根據(jù)仿真結(jié)果得出以下結(jié)論:

1) 具有減旋船尾結(jié)構(gòu)的制導(dǎo)炮彈氣動(dòng)外形設(shè)計(jì)合理,氣動(dòng)參數(shù)滿足穩(wěn)定飛行條件,可為遠(yuǎn)程尾控制導(dǎo)炮彈設(shè)計(jì)提供參考依據(jù)。

2) 遠(yuǎn)程尾控制導(dǎo)炮彈彈頭所承受的壓力值較大,在材料處理方面,需要對(duì)彈頭的材料進(jìn)行特殊處理,以便于更好地適應(yīng)飛行環(huán)境。