箔片轉動數學建模及仿真分析

付俊博,金忠慶

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

1 引 言

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

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

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

基于此,本文對箔片的轉動現象進行深入研究,建立了相對精確的箔片轉動模型,并將仿真結果與實驗對比,以驗證模型的準確性。然后,應用建立的模型對箔片的轉動特性進行分析,并仿真分析了箔片云的擴散運動過程。

2 箔片運動建模

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

圖1 箔片各坐標系之間的關系

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

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

圖2 箔片受力分析

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

(1)

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

箔片的運動學方程為[9-10]:

(2)

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

(3)

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

(4)

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

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

(5)

因此可以得到速度滾轉角γs、迎角α的表達式為:

(6)

(7)

3 箔片轉動建模

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

3.1 箔片轉動模型

在箔片的運動過程中,其所受的合外力矩主要由升力產生,力矩方程在速度坐標系內可表示為[11]:

(8)

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

圖3 氣動力和力矩分析

(9)

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

(10)

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

Ya=|Y|nYay/|nYay|

(11)

3.2 氣動阻尼的求解

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

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

圖4 箔片轉動力矩分析

則轉動過程中,dS所受的氣動阻力為:

(12)

(13)

聯立以上公式,則可以得到箔片的運動模型和轉動模型。

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

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

(14)

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

(15)

4 箔片轉動實驗與模型驗證

4.1 箔片氣動力系數計算

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

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

圖5 多面體網格及箔片邊界層網格示意圖

表1 箔片氣動力系數計算結果

4.2 箔片轉動實驗

由于箔片運動過程中其姿態、位置變化較為明顯,且與其初始狀態有關,實驗中難以測得其運動過程中的角速度變化規律。因此本節設計了箔片在質心靜止狀態下的轉動實驗,箔片在實驗中只存在轉動過程質心不存在運動,并將實驗結果與建立的轉動模型計算結果進行對比。

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

圖6 箔片轉速測試試驗示意圖

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

圖7 實驗結果與仿真結果對比

由圖7中的結果可知,箔片初始迎角α=90°時,箔片基本處于靜止狀態,幾乎沒有轉動運動,其轉動角度速度遠低于α=10°。由此可知,箔片的轉動角速度與其初始的角度關系很大。由于箔片在轉動過程中其氣動力矩周期性變化,因此箔片在轉動過程中其角速度時刻變化,變化曲線類似于正弦曲線。而箔片在轉動過程中氣動阻尼力矩與轉動方向相反,使得箔片的轉動角速度逐漸較小。

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

5 仿真結果及現象分析

5.1 箔片轉動仿真

設載機水平飛行,飛行高度為8 km,飛行馬赫數Ma=0.8。箔片1-4從載機上垂直向上拋撒,拋撒速度為20 m/s,箔片直徑均為50 mm。箔片1-4初始姿態角完全一致,初始俯仰角?=45°,初始偏航角ψ=-30°。箔片1發射后依據文中建立的模型對其運動和轉動過程進行求解,箔片2、箔片3、箔片4轉動角速度不變,分別為50 rad/s、100 rad/s以及200 rad/s。以箔片拋撒時刻其位置在水平面上的投影為坐標原點,載機飛行方向為正X軸建立坐標系。比較四個箔片運動和轉動過程,如圖8、圖9所示。

圖8 箔片運動軌跡對比圖

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

圖9 箔片特性對比圖

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

5.2 箔片云運動擴散仿真

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

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

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

圖10 0.5 s時刻箔片云擴散圖

圖11 1 s時刻箔片云擴散圖

6 結 論

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

(1)箔片的轉動角速度與其初始迎角有關,初始迎角接近90°時,其轉動角速度最??；

(2)箔片的轉動角速度與運動速度有關,且隨著運動速度的下降箔片轉動角速度越來越小；

(3)箔片在運動過程中其轉動角速度呈高頻大幅震蕩狀態,且當來流風速較小時,其角速度變化曲線類似于正弦函數；

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