箔片轉(zhuǎn)動數(shù)學(xué)建模及仿真分析

付俊博,金忠慶

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

1 引 言

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

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

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

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

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

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

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

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

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

圖2 箔片受力分析

箔片的動力學(xué)方程可以簡化為[8]:

(1)

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

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

(2)

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

(3)

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

(4)

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

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

(5)

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

(6)

(7)

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

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

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

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

(8)

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

圖3 氣動力和力矩分析

(9)

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

(10)

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

Ya=|Y|nYay/|nYay|

(11)

3.2 氣動阻尼的求解

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

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

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

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

(12)

(13)

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

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

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

(14)

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

(15)

4 箔片轉(zhuǎn)動實(shí)驗與模型驗證

4.1 箔片氣動力系數(shù)計算

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

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

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

表1 箔片氣動力系數(shù)計算結(jié)果

4.2 箔片轉(zhuǎn)動實(shí)驗

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

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

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

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

圖7 實(shí)驗結(jié)果與仿真結(jié)果對比

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

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

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

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

設(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ù)文中建立的模型對其運(yùn)動和轉(zhuǎn)動過程進(jìn)行求解,箔片2、箔片3、箔片4轉(zhuǎn)動角速度不變,分別為50 rad/s、100 rad/s以及200 rad/s。以箔片拋撒時刻其位置在水平面上的投影為坐標(biāo)原點(diǎn),載機(jī)飛行方向為正X軸建立坐標(biāo)系。比較四個箔片運(yùn)動和轉(zhuǎn)動過程,如圖8、圖9所示。

圖8 箔片運(yùn)動軌跡對比圖

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

圖9 箔片特性對比圖

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

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

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

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

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

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

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

6 結(jié) 論

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

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

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

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

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