付俊博,金忠慶

(1.吉林交通職業(yè)技術(shù)學(xué)院,吉林 長(zhǎng)春 130000；2.空軍航空大學(xué),吉林 長(zhǎng)春 130000)

1 引 言

箔片型面源紅外誘餌由上千個(gè)箔片壓縮而成,誘餌發(fā)射后箔片迅速擴(kuò)散形成箔片云。箔片在運(yùn)動(dòng)過(guò)程中,只受重力和氣動(dòng)力的影響。由于箔片的氣動(dòng)中心與重心不重合,因此存在氣動(dòng)力矩,使得箔片在運(yùn)動(dòng)過(guò)程中還伴隨著轉(zhuǎn)動(dòng)。箔片的轉(zhuǎn)動(dòng)使得其所受氣動(dòng)力呈周期性變化,影響箔片的運(yùn)動(dòng)過(guò)程,并決定著箔片云的運(yùn)動(dòng)擴(kuò)散性能。但由于箔片的轉(zhuǎn)動(dòng)過(guò)程中存在著障礙轉(zhuǎn)動(dòng)運(yùn)動(dòng)進(jìn)行的氣動(dòng)阻尼,使得箔片的轉(zhuǎn)動(dòng)過(guò)程較為復(fù)雜。

目前,國(guó)內(nèi)外的學(xué)者對(duì)箔片的運(yùn)動(dòng)及轉(zhuǎn)動(dòng)進(jìn)行了相應(yīng)的建模研究。鄒濤將箔片運(yùn)動(dòng)過(guò)程中的轉(zhuǎn)動(dòng)角速度認(rèn)為是定值,并建立了箔片的運(yùn)動(dòng)模型[1]。黃蓓應(yīng)用非結(jié)構(gòu)動(dòng)網(wǎng)格方法研究了箔片非定常分離過(guò)程,并設(shè)計(jì)了實(shí)驗(yàn)對(duì)數(shù)值計(jì)算結(jié)果進(jìn)行驗(yàn)證,其仿真計(jì)算結(jié)果表明箔片在降落過(guò)程中存在滯空和轉(zhuǎn)動(dòng)運(yùn)動(dòng)[2]。黃蓓還應(yīng)用高速攝像機(jī)設(shè)計(jì)了一個(gè)箔片運(yùn)動(dòng)試驗(yàn),依靠火藥燃燒產(chǎn)生的推力將箔片推出發(fā)射筒,試驗(yàn)結(jié)果表明箔片在高速運(yùn)動(dòng)中還存在著高速轉(zhuǎn)動(dòng)[3]。

由于箔片的轉(zhuǎn)動(dòng)過(guò)程較為復(fù)雜,涉及多學(xué)科、交叉理論。其內(nèi)部既包含了數(shù)值積分算法、數(shù)學(xué)建模方法又涵蓋了空氣動(dòng)力學(xué)、飛行力學(xué)、流體力學(xué)、理論力學(xué)、非定常流動(dòng)等相關(guān)領(lǐng)域的知識(shí)[4-5]。同時(shí)箔片轉(zhuǎn)動(dòng)過(guò)程還需要考慮氣動(dòng)阻尼對(duì)箔片轉(zhuǎn)動(dòng)力矩的影響。因此,國(guó)內(nèi)外的學(xué)者對(duì)箔片的運(yùn)動(dòng)模型研究較為充分,而對(duì)箔片的轉(zhuǎn)動(dòng)過(guò)程及轉(zhuǎn)動(dòng)模型的研究還不夠深入,相關(guān)模型建立的較為簡(jiǎn)單,不能夠真實(shí)反應(yīng)箔片轉(zhuǎn)動(dòng)過(guò)程中其角速度的變化規(guī)律,進(jìn)而不能夠得到箔片云的真實(shí)擴(kuò)散運(yùn)動(dòng)過(guò)程。

基于此,本文對(duì)箔片的轉(zhuǎn)動(dòng)現(xiàn)象進(jìn)行深入研究,建立了相對(duì)精確的箔片轉(zhuǎn)動(dòng)模型,并將仿真結(jié)果與實(shí)驗(yàn)對(duì)比,以驗(yàn)證模型的準(zhǔn)確性。然后,應(yīng)用建立的模型對(duì)箔片的轉(zhuǎn)動(dòng)特性進(jìn)行分析,并仿真分析了箔片云的擴(kuò)散運(yùn)動(dòng)過(guò)程。

2 箔片運(yùn)動(dòng)建模

由于箔片的運(yùn)動(dòng)需要在地軸系內(nèi)計(jì)算,而箔片的動(dòng)力學(xué)方程需要在航跡坐標(biāo)系內(nèi)計(jì)算,因此首先需要建立箔片的地軸系Oxgygzg、航跡坐標(biāo)系Oxhyhzh以及速度坐標(biāo)系Oxayaza,各坐標(biāo)系的定義見(jiàn)文獻(xiàn)[6]～[7],坐標(biāo)系之間的關(guān)系如圖1所示。

圖1 箔片各坐標(biāo)系之間的關(guān)系

定義Oxh軸與水平面之間的夾角為航跡俯仰角θ,運(yùn)動(dòng)方向向上為正；Oya軸與Oyh軸之間的夾角為速度滾轉(zhuǎn)角γs,速度坐標(biāo)系繞Oxh軸向右旋轉(zhuǎn)時(shí),γs為正；Oxh軸水平面的投影與Oxgyg之間的夾角為航向角ψs,Oxh軸左偏為正。定義垂直于箔片平面且過(guò)圓心的單位向量為箔片軸心線(xiàn)向量nd,Oya軸與Oxa軸所組成的平面為氣流對(duì)稱(chēng)平面。則向量nd及箔片運(yùn)動(dòng)方向都在氣流對(duì)稱(chēng)平面內(nèi)。

由于箔片在整個(gè)運(yùn)動(dòng)過(guò)程中都是無(wú)動(dòng)力的,因此在箔片運(yùn)動(dòng)過(guò)程中只受氣動(dòng)力和重力的影響。由于箔片存在特殊的對(duì)稱(chēng)性,且其表面較為平滑,因此氣動(dòng)側(cè)力予以忽略。則其運(yùn)動(dòng)過(guò)程中只受到氣動(dòng)阻力X、氣動(dòng)升力Y的作用,且其受到的氣動(dòng)力都在氣流對(duì)稱(chēng)平面內(nèi),如圖2所示。

圖2 箔片受力分析

箔片的動(dòng)力學(xué)方程可以簡(jiǎn)化為[8]:

(1)

式中,X、Y為氣動(dòng)力;cx、cy為相應(yīng)的氣動(dòng)力系數(shù);V為箔片的運(yùn)動(dòng)速度;γs為箔片的速度滾轉(zhuǎn)角;ψs為箔片的航向角;θ為箔片航跡俯仰角;ρ為當(dāng)前環(huán)境的大氣密度。由于箔片存在特殊的對(duì)稱(chēng)性,因此其向各個(gè)方向運(yùn)動(dòng)時(shí),cx、cy都只與箔片的迎角和運(yùn)動(dòng)速度有關(guān)。

箔片的運(yùn)動(dòng)學(xué)方程為[9-10]:

(2)

式中,(xg,yg,zg)為箔片在地軸系中的坐標(biāo)。則:

(3)

式中,ndx、ndy、ndz為nd在地軸系中的三軸分量;nV為運(yùn)動(dòng)速度方向的單位向量;nVx、nVy、nVz為其在地軸系中的三軸分量。則:

(4)

式中,ψ為偏航角;?為俯仰角。

定義nY為氣動(dòng)升力的單位向量,nYx、nYy、nYz為其在地軸系上的三軸分量。則:

(5)

因此可以得到速度滾轉(zhuǎn)角γs、迎角α的表達(dá)式為:

(6)

(7)

3 箔片轉(zhuǎn)動(dòng)建模

箔片的運(yùn)動(dòng)過(guò)程中,由于其壓力中心與箔片的幾何中心不重合,因此在其運(yùn)動(dòng)過(guò)程中受到的合外力矩不為零。作用在其上的合外力矩使得箔片高速轉(zhuǎn)動(dòng),箔片的軸心線(xiàn)向量nd、迎角α等參數(shù)周期性地進(jìn)行變換,最終導(dǎo)致箔片所受的氣動(dòng)力不斷變化。

3.1 箔片轉(zhuǎn)動(dòng)模型

在箔片的運(yùn)動(dòng)過(guò)程中,其所受的合外力矩主要由升力產(chǎn)生,力矩方程在速度坐標(biāo)系內(nèi)可表示為[11]:

(8)

式中,J為箔片轉(zhuǎn)動(dòng)慣量；R為箔片半徑；xF為箔片壓力中心與箔片幾何中心之間的距離[12]；Md為氣動(dòng)阻尼力矩；ω為轉(zhuǎn)動(dòng)角速度；Ya為氣動(dòng)升力在速度坐標(biāo)系中的值。箔片轉(zhuǎn)動(dòng)過(guò)程中所受的氣動(dòng)力和力矩示意圖如圖3所示。

圖3 氣動(dòng)力和力矩分析

(9)

定義nYa為升力向量nY在速度坐標(biāo)系中的值,則其值可以表示為:

(10)

nYay為nYa在Oya軸上的值,則Ya可以表示為:

Ya=|Y|nYay/|nYay|

(11)

3.2 氣動(dòng)阻尼的求解

箔片的轉(zhuǎn)動(dòng)過(guò)程中,同時(shí)還伴隨著障礙轉(zhuǎn)動(dòng)運(yùn)動(dòng)進(jìn)行的氣動(dòng)阻尼力矩Md。如果將箔片的運(yùn)動(dòng)過(guò)程分解為質(zhì)心平移運(yùn)動(dòng)和繞質(zhì)心轉(zhuǎn)動(dòng)運(yùn)動(dòng),則箔片的質(zhì)心平移運(yùn)動(dòng)可以由箔片的運(yùn)動(dòng)方程求得,而箔片繞質(zhì)心的轉(zhuǎn)動(dòng)運(yùn)動(dòng)主要由氣動(dòng)升力產(chǎn)生的力矩M以及氣動(dòng)阻尼力矩Md求得。

設(shè)箔片的轉(zhuǎn)動(dòng)角速度為ω,dS為箔片表面到箔片中心距離為r的面源,如圖4所示。

圖4 箔片轉(zhuǎn)動(dòng)力矩分析

則轉(zhuǎn)動(dòng)過(guò)程中,dS所受的氣動(dòng)阻力為:

(12)

(13)

聯(lián)立以上公式,則可以得到箔片的運(yùn)動(dòng)模型和轉(zhuǎn)動(dòng)模型。

由于箔片的角速度垂直于氣流對(duì)稱(chēng)平面,氣流對(duì)稱(chēng)平面又重合于速度坐標(biāo)系的Oxaya平面,則Δt時(shí)刻后,軸心線(xiàn)向量nd在速度坐標(biāo)系內(nèi)可以表示為:

nda=[-sin(α+ωΔt),cos(α+ωΔt),0]

(14)

式中,nda為nd在速度坐標(biāo)系中的值,因此nd可以表示為:

(15)

4 箔片轉(zhuǎn)動(dòng)實(shí)驗(yàn)與模型驗(yàn)證

4.1 箔片氣動(dòng)力系數(shù)計(jì)算

本節(jié)采用CFD對(duì)箔片的氣動(dòng)力系數(shù)進(jìn)行數(shù)值計(jì)算。計(jì)算模型的選取如下:控制方程為Navier-Stokes方程,并采用PIM-PLE算法對(duì)其進(jìn)行求解；空間離散方法為空間二階精度線(xiàn)性插值算法和基于有限體積的空間離散算法；時(shí)間離散方法為二階精度的后向差分算法。初始條件設(shè)置如下:箔片表面設(shè)置成無(wú)滑移壁面,計(jì)算域的出口、入口以及其他外壁全部設(shè)置成壓力遠(yuǎn)場(chǎng),湍流計(jì)算模型選取SST湍流模型[13]。

采用多面體網(wǎng)格對(duì)箔片的氣動(dòng)力系數(shù)進(jìn)行計(jì)算,計(jì)算域?yàn)檫呴L(zhǎng)5000 mm的正方體,箔片位于計(jì)算域的正中心。箔片直徑為50 mm,厚度1 mm,邊界層一共設(shè)置12層,第一層厚度為10-6m,箔片表面及邊界層內(nèi)的網(wǎng)格如圖5所示。計(jì)算域內(nèi)的網(wǎng)格數(shù)總數(shù)為213萬(wàn),計(jì)算精度能夠滿(mǎn)足要求。則高度H=0 km、環(huán)境溫度T=288.15 K、大氣壓力P=101325 Pa、入口來(lái)流速度為100 m/s時(shí),不同迎角下箔片氣動(dòng)力系數(shù)CFD的計(jì)算結(jié)果如表1所示。

圖5 多面體網(wǎng)格及箔片邊界層網(wǎng)格示意圖

表1 箔片氣動(dòng)力系數(shù)計(jì)算結(jié)果

4.2 箔片轉(zhuǎn)動(dòng)實(shí)驗(yàn)

由于箔片運(yùn)動(dòng)過(guò)程中其姿態(tài)、位置變化較為明顯,且與其初始狀態(tài)有關(guān),實(shí)驗(yàn)中難以測(cè)得其運(yùn)動(dòng)過(guò)程中的角速度變化規(guī)律。因此本節(jié)設(shè)計(jì)了箔片在質(zhì)心靜止?fàn)顟B(tài)下的轉(zhuǎn)動(dòng)實(shí)驗(yàn),箔片在實(shí)驗(yàn)中只存在轉(zhuǎn)動(dòng)過(guò)程質(zhì)心不存在運(yùn)動(dòng),并將實(shí)驗(yàn)結(jié)果與建立的轉(zhuǎn)動(dòng)模型計(jì)算結(jié)果進(jìn)行對(duì)比。

箔片的直徑為50 mm,厚度1 mm,尺寸與CFD數(shù)值計(jì)算的尺寸保持一致。在箔片表面固定有一個(gè)穿過(guò)其圓形緊貼箔片表面的轉(zhuǎn)軸,轉(zhuǎn)軸固定在一個(gè)低速風(fēng)洞內(nèi),如圖6所致。則當(dāng)風(fēng)洞工作時(shí),箔片只能夠繞轉(zhuǎn)軸轉(zhuǎn)動(dòng)而不能夠自由運(yùn)動(dòng)。箔片的右側(cè)放置一個(gè)高速相機(jī),距其0.5 m。箔片轉(zhuǎn)動(dòng)時(shí),通過(guò)高速相機(jī)拍攝到的每一幀畫(huà)面,計(jì)算出箔片的轉(zhuǎn)動(dòng)角度,進(jìn)而得到其轉(zhuǎn)動(dòng)角速度。

圖6 箔片轉(zhuǎn)速測(cè)試試驗(yàn)示意圖

將實(shí)驗(yàn)風(fēng)速設(shè)置為20 m/s,為了與實(shí)驗(yàn)一致仿真中將箔片質(zhì)心位置設(shè)置為靜止不變,只考慮其轉(zhuǎn)動(dòng)問(wèn)題不考慮運(yùn)動(dòng)問(wèn)題,即認(rèn)為公式(1)中的V=20、θ=0、ψs=0不變,公式(2)中的xg、yg、zg始終等于零。則初始迎角α=90°和α=10°時(shí),實(shí)驗(yàn)結(jié)果和轉(zhuǎn)動(dòng)模型計(jì)算的結(jié)果對(duì)比圖如圖7所示。

圖7 實(shí)驗(yàn)結(jié)果與仿真結(jié)果對(duì)比

由圖7中的結(jié)果可知,箔片初始迎角α=90°時(shí),箔片基本處于靜止?fàn)顟B(tài),幾乎沒(méi)有轉(zhuǎn)動(dòng)運(yùn)動(dòng),其轉(zhuǎn)動(dòng)角度速度遠(yuǎn)低于α=10°。由此可知,箔片的轉(zhuǎn)動(dòng)角速度與其初始的角度關(guān)系很大。由于箔片在轉(zhuǎn)動(dòng)過(guò)程中其氣動(dòng)力矩周期性變化,因此箔片在轉(zhuǎn)動(dòng)過(guò)程中其角速度時(shí)刻變化,變化曲線(xiàn)類(lèi)似于正弦曲線(xiàn)。而箔片在轉(zhuǎn)動(dòng)過(guò)程中氣動(dòng)阻尼力矩與轉(zhuǎn)動(dòng)方向相反,使得箔片的轉(zhuǎn)動(dòng)角速度逐漸較小。

從圖7中的對(duì)比結(jié)果可知,轉(zhuǎn)動(dòng)模型仿真的角速度略大于實(shí)驗(yàn)結(jié)果,這主要與實(shí)驗(yàn)中箔片的轉(zhuǎn)軸與軸承之間還存在摩擦力有關(guān),仿真結(jié)果與實(shí)驗(yàn)結(jié)果總體一致,因此建立的轉(zhuǎn)動(dòng)模型較為可信,能夠滿(mǎn)足精度的要求。

5 仿真結(jié)果及現(xiàn)象分析

5.1 箔片轉(zhuǎn)動(dòng)仿真

設(shè)載機(jī)水平飛行,飛行高度為8 km,飛行馬赫數(shù)Ma=0.8。箔片1-4從載機(jī)上垂直向上拋撒,拋撒速度為20 m/s,箔片直徑均為50 mm。箔片1-4初始姿態(tài)角完全一致,初始俯仰角?=45°,初始偏航角ψ=-30°。箔片1發(fā)射后依據(jù)文中建立的模型對(duì)其運(yùn)動(dòng)和轉(zhuǎn)動(dòng)過(guò)程進(jìn)行求解,箔片2、箔片3、箔片4轉(zhuǎn)動(dòng)角速度不變,分別為50 rad/s、100 rad/s以及200 rad/s。以箔片拋撒時(shí)刻其位置在水平面上的投影為坐標(biāo)原點(diǎn),載機(jī)飛行方向?yàn)檎齒軸建立坐標(biāo)系。比較四個(gè)箔片運(yùn)動(dòng)和轉(zhuǎn)動(dòng)過(guò)程,如圖8、圖9所示。

圖8 箔片運(yùn)動(dòng)軌跡對(duì)比圖

由圖8的仿真結(jié)果可知,箔片的轉(zhuǎn)動(dòng)角速度對(duì)箔片的運(yùn)動(dòng)過(guò)程影響較大,箔片2、箔片3和箔片4相對(duì)于箔片1向Z軸運(yùn)動(dòng)的距離都較大,且箔片的轉(zhuǎn)動(dòng)角速度越小其向Z軸運(yùn)動(dòng)的距離越長(zhǎng)。由于箔片的初始偏航角為負(fù),因此箔片氣動(dòng)力在Z軸正方向具有一定的分量,導(dǎo)致箔片具有沿Z軸正方向運(yùn)動(dòng)的趨勢(shì)。箔片的初始轉(zhuǎn)速越小,其所受氣動(dòng)力周期性變化時(shí)間越長(zhǎng),最終導(dǎo)致箔片在Z軸上運(yùn)動(dòng)距離越長(zhǎng)。

圖9 箔片特性對(duì)比圖

圖9(a)的仿真結(jié)果可知,箔片拋撒后在氣動(dòng)阻力的作用下速度迅速降低,2.5 s后箔片幾乎停止運(yùn)動(dòng)。圖9(b)的仿真可知,箔片1的初始轉(zhuǎn)動(dòng)角速度較大,但隨著速度的下降轉(zhuǎn)動(dòng)角速度迅速減小。由于拋撒時(shí)刻箔片的運(yùn)動(dòng)速度較大,因此其所受的初始轉(zhuǎn)動(dòng)力矩很大,導(dǎo)致箔片的初始轉(zhuǎn)動(dòng)速度較大。隨著箔片運(yùn)動(dòng)速度的下降,箔片受到的轉(zhuǎn)動(dòng)力矩減小,最終導(dǎo)致箔片轉(zhuǎn)動(dòng)角速度越來(lái)越小。當(dāng)箔片轉(zhuǎn)動(dòng)到α<0°時(shí),氣動(dòng)力矩與轉(zhuǎn)動(dòng)角速度相反,在氣動(dòng)阻尼和氣動(dòng)力矩的作用下,箔片轉(zhuǎn)動(dòng)角速度迅速降低。因此從圖9(b)還可以看出,箔片1的轉(zhuǎn)動(dòng)角速度呈高頻大幅震蕩狀態(tài),且振幅隨著轉(zhuǎn)動(dòng)角速度的減小而下降。

5.2 箔片云運(yùn)動(dòng)擴(kuò)散仿真

面源紅外誘餌發(fā)射后由于大氣擾動(dòng)等隨機(jī)因素的影響,箔片的初始姿態(tài)有一定的差別。經(jīng)過(guò)與實(shí)測(cè)數(shù)據(jù)的對(duì)比分析,近似認(rèn)為箔片初始俯仰角?0服從U(0,π/2)的均勻分布,初始偏航角ψ0服從U(0,2π)的均勻分布。由于箔片的初始姿態(tài)已知,則根據(jù)文中建立的箔片運(yùn)動(dòng)模型和轉(zhuǎn)動(dòng)模型可以得到箔片云的整體運(yùn)動(dòng)擴(kuò)散過(guò)程。

設(shè)載機(jī)水平飛行,飛行高度為8 km,飛行馬赫數(shù)Ma=0.8。面源紅外誘餌由一千個(gè)箔片壓縮而成,箔片直徑均為50 mm,誘餌發(fā)射初始速度為20 m/s,速度相對(duì)于載機(jī)垂直向上。則誘餌發(fā)射后,箔片云的運(yùn)動(dòng)擴(kuò)散過(guò)程如圖10、圖11所示。

從圖10、圖11中的仿真結(jié)果可以看出,面源誘餌發(fā)射后,箔片云大致呈與X軸存在一定夾角的錐形分布。箔片云擴(kuò)散速度較快,0.5 s時(shí)刻,箔片云在X軸擴(kuò)散長(zhǎng)度就達(dá)到了59 m,Y軸長(zhǎng)度達(dá)37 m,Z軸擴(kuò)散長(zhǎng)度達(dá)40 m。1 s時(shí)刻,箔片云繼續(xù)擴(kuò)散,但擴(kuò)散速度明顯減慢。

圖10 0.5 s時(shí)刻箔片云擴(kuò)散圖

圖11 1 s時(shí)刻箔片云擴(kuò)散圖

6 結(jié) 論

本文研究了箔片在運(yùn)動(dòng)過(guò)程中的轉(zhuǎn)動(dòng)問(wèn)題。通過(guò)定義箔片的軸心線(xiàn)向量,確定了箔片迎角和速度滾轉(zhuǎn)角的求解方法,并建立了箔片的運(yùn)動(dòng)模型。然后,建立了箔片的轉(zhuǎn)動(dòng)模型,并研究了在轉(zhuǎn)動(dòng)過(guò)程中氣動(dòng)阻尼的求解方法。通過(guò)CFD流場(chǎng)計(jì)算方法,計(jì)算箔片的氣動(dòng)力系數(shù),并對(duì)比箔片轉(zhuǎn)速的實(shí)驗(yàn)結(jié)果和模型仿真結(jié)果,驗(yàn)證了模型的可信性。最終,仿真分析了箔片的轉(zhuǎn)動(dòng)特性及箔片云的運(yùn)動(dòng)擴(kuò)散特性。所得到的主要結(jié)論如下:

(1)箔片的轉(zhuǎn)動(dòng)角速度與其初始迎角有關(guān),初始迎角接近90°時(shí),其轉(zhuǎn)動(dòng)角速度最??；

(2)箔片的轉(zhuǎn)動(dòng)角速度與運(yùn)動(dòng)速度有關(guān),且隨著運(yùn)動(dòng)速度的下降箔片轉(zhuǎn)動(dòng)角速度越來(lái)越??；

(3)箔片在運(yùn)動(dòng)過(guò)程中其轉(zhuǎn)動(dòng)角速度呈高頻大幅震蕩狀態(tài),且當(dāng)來(lái)流風(fēng)速較小時(shí),其角速度變化曲線(xiàn)類(lèi)似于正弦函數(shù)；

(4)箔片云擴(kuò)散速度較快,誘餌發(fā)射1 s后箔片云已基本擴(kuò)散成型,箔片云大致呈與載機(jī)飛行方向存在一定夾角的錐形分布,且箔片云在載機(jī)飛行方向上擴(kuò)散能力最強(qiáng)。