基于運(yùn)動(dòng)捕捉系統(tǒng)的小型有動(dòng)力無人機(jī)氣動(dòng)特性分析

趙嶷飛, 張悅宸, 趙欣宇

(1.中國民航大學(xué)空中交通管理學(xué)院，天津 300300；2.中國民航大學(xué)空中交通管理學(xué)院，天津 300300；3.民航機(jī)場成都電子工程設(shè)計(jì)有限責(zé)任公司北京分公司，北京 100015)

隨著無人機(jī)逐漸成為社會(huì)熱點(diǎn)，科研工作者們對(duì)無人機(jī)的研究日漸成熟，然而小型固定翼飛機(jī)的空氣動(dòng)力學(xué)性能尚未有系統(tǒng)標(biāo)準(zhǔn)的獲取方法。小型無人機(jī)飛行特性需要更準(zhǔn)確地研究數(shù)據(jù)來進(jìn)一步了解其空氣動(dòng)力學(xué)模型，然而因?yàn)轶w積小、質(zhì)量輕，阻礙了常規(guī)檢測系統(tǒng)的使用，使得氣動(dòng)數(shù)據(jù)難以測量?，F(xiàn)有的方式多以風(fēng)洞實(shí)驗(yàn)為主，但操作復(fù)雜、成本較高、對(duì)所需硬件要求嚴(yán)格等原因成為很多研究者難以克服的困難。因此，如何方便、經(jīng)濟(jì)地獲取小型有動(dòng)力無人機(jī)的空氣動(dòng)力學(xué)數(shù)據(jù)變得至關(guān)重要，現(xiàn)將在此背景下提出一種普遍易行的獲取方法。

運(yùn)動(dòng)捕捉系統(tǒng)用于測定非定常效應(yīng)，而不需要風(fēng)洞實(shí)驗(yàn)中所需要的復(fù)雜儀器。文獻(xiàn)[1]描述了一種直接從收集的軌跡中提取飛機(jī)氣動(dòng)特性的技術(shù)。文獻(xiàn)[2]介紹了一種從航跡數(shù)據(jù)中提取飛行器氣動(dòng)特性的方法。利用現(xiàn)代運(yùn)動(dòng)跟蹤技術(shù)，可以從飛行數(shù)據(jù)確定作用在飛機(jī)上的力。文獻(xiàn)[3]提出了一套物理啟發(fā)的基礎(chǔ)函數(shù)，使系統(tǒng)識(shí)別一個(gè)非線性空氣動(dòng)力學(xué)模型沿棲息軌跡。數(shù)據(jù)的收集使用動(dòng)作捕捉系統(tǒng)，關(guān)鍵是允許自由飛行數(shù)據(jù)從真正的系統(tǒng)軌跡收集。文獻(xiàn)[4]采用在運(yùn)動(dòng)捕捉系統(tǒng)中通過真實(shí)的運(yùn)動(dòng)學(xué)飛行數(shù)據(jù)獲得準(zhǔn)確的飛機(jī)模型。文獻(xiàn)[5]中一個(gè)舷外運(yùn)動(dòng)跟蹤系統(tǒng)捕獲了一架小型無線電控制飛機(jī)飛行時(shí)的運(yùn)動(dòng)軌跡，對(duì)記錄的運(yùn)動(dòng)軌跡時(shí)間歷程進(jìn)行了分析，確定了飛機(jī)在無動(dòng)力飛行時(shí)的氣動(dòng)特性。文獻(xiàn)[6]為了收集微型飛行器的非定常飛行數(shù)據(jù)，研究人員使用了一個(gè)離體運(yùn)動(dòng)跟蹤系統(tǒng)來捕獲飛行器的自由飛行軌跡，該參數(shù)模擬了攻角快速變化時(shí)的氣動(dòng)滯后，從而捕捉到了升力、阻力和力矩系數(shù)數(shù)據(jù)中動(dòng)態(tài)失速的影響。文獻(xiàn)[7]實(shí)驗(yàn)期間飛行數(shù)據(jù)由傳感器數(shù)據(jù)采集系統(tǒng)記錄，給出了用于分析飛行試驗(yàn)數(shù)據(jù)的數(shù)據(jù)分析方法。

對(duì)于上述已有的研究，多數(shù)研究對(duì)象為小型無動(dòng)力無人機(jī)，對(duì)于有動(dòng)力無人機(jī)的氣動(dòng)特性在不適用風(fēng)洞模型的情況下尚未有明確可供參考且可供通用的測量方法。因此，在小型無動(dòng)力無人機(jī)氣動(dòng)參數(shù)分析的基礎(chǔ)上，對(duì)小型有動(dòng)力無人機(jī)在不同功率輸出下的推力值做出了測量，得到無人機(jī)的推力曲線。建立小型有動(dòng)力無人機(jī)運(yùn)動(dòng)模型，使用MATLAB仿真得到飛機(jī)的空氣動(dòng)力學(xué)參數(shù)，為小型固定翼飛機(jī)的空氣動(dòng)力學(xué)性能提供方便、經(jīng)濟(jì)的新思路。

1 無人機(jī)物理參數(shù)

采用偉力F959后退式固定翼滑翔機(jī)(圖1)。F959翼展75.0 cm，質(zhì)量118 g，三視圖如圖2所示。機(jī)身由輕質(zhì)泡沫材質(zhì)構(gòu)成，手工組裝而成。飛機(jī)裝配前，分別測量其各部件的質(zhì)量和位置，并且根據(jù)質(zhì)量慣性矩公式[式(1)～式(4)]，估算其慣性力矩，結(jié)果如表1所示。

表1 飛機(jī)物理特性

圖1 實(shí)驗(yàn)飛機(jī)實(shí)體

圖2 實(shí)驗(yàn)飛機(jī)三視圖

(1)

(2)

(3)

(4)

式中：Ixx、Iyy、Izz分別為航空器沿x、y、z軸的慣性矩；Ixz為航空器在xz平面的慣性積;ρ為剛體密度；xi、yi為質(zhì)量微元mi在x、y軸上的位置信息。

2 小型有動(dòng)力無人機(jī)空氣動(dòng)力學(xué)參數(shù)分析

2.1 推力測量

F959使用的電池容量為7.4 V，330 mA。由于飛行過程中發(fā)動(dòng)機(jī)的輸出功率不同，導(dǎo)致飛機(jī)在不同時(shí)刻受到的發(fā)動(dòng)機(jī)推力都不盡相同，從而飛機(jī)每一時(shí)刻所受到的升阻力也會(huì)隨之變化。因此，如何獲得飛機(jī)行進(jìn)時(shí)的不同推力值，成為分析小型有動(dòng)力無人機(jī)空氣動(dòng)力學(xué)參數(shù)的關(guān)鍵。

采用直接測量的方式對(duì)飛機(jī)所受推力進(jìn)行測量。F959小型無人機(jī)遙控器如圖3所示，其中液晶屏所示CH3通道可顯示油門桿輸出效率，并在實(shí)驗(yàn)結(jié)束后選擇性輸出其飛行進(jìn)程中所產(chǎn)生的數(shù)據(jù)，即在飛行實(shí)驗(yàn)后可導(dǎo)出油門桿實(shí)時(shí)輸出功率，如在已知油門桿在每一節(jié)點(diǎn)輸出功率對(duì)應(yīng)的飛機(jī)推力的情況下，即可獲得小飛機(jī)在飛行過程中所受發(fā)動(dòng)機(jī)推力。獲得每一節(jié)點(diǎn)的飛機(jī)推力方法如下。

圖3 遙控器

首先將輕質(zhì)硬紙片固定在飛機(jī)機(jī)頭處，使得小飛機(jī)機(jī)頭垂直朝下仍可立與水平地面。選擇測量精度為0.1 g(即10-3N)的電子秤，去除輕質(zhì)紙片質(zhì)量及飛機(jī)自重后將小型飛機(jī)立與水平電子秤測量板上，保持飛機(jī)豎直穩(wěn)定并且不受外力影響，緩慢推動(dòng)遙控器油門桿，電子秤顯示數(shù)值即為飛機(jī)所受實(shí)時(shí)推力。以1%為單位記錄小型飛機(jī)所受推力數(shù)值，并繪制飛機(jī)推力曲線(圖4)。

圖4 F959推力曲線

從推力曲線可以看出，推力隨輸出效率的增長而增長，符合發(fā)動(dòng)機(jī)運(yùn)行規(guī)律。在輸出效率為0～30%的推力增長幅度大于30%～90%，90%后的推力變化不明顯。

2.2 建立氣動(dòng)參數(shù)模型

對(duì)于試驗(yàn)后期數(shù)據(jù)的處理，首先是通過基于軌跡時(shí)程中相鄰點(diǎn)的線性插值估計(jì)姿態(tài)和位置來填補(bǔ)Vicon系統(tǒng)由于相機(jī)擺放位置、角度等問題導(dǎo)致的觀測數(shù)據(jù)缺失。

在填充這些丟失的數(shù)據(jù)點(diǎn)之后，就可以對(duì)原始數(shù)據(jù)進(jìn)行平滑和微分。由于Vicon系統(tǒng)不是100%穩(wěn)定的，系統(tǒng)的原始測量數(shù)據(jù)包括了飛行實(shí)驗(yàn)中不可避免的散射、不確定性和噪聲。原始測量結(jié)果后期用于計(jì)算速度和加速度，并進(jìn)一步確定氣動(dòng)系數(shù)。在這些計(jì)算中，噪聲將被放大和傳播。因此，在對(duì)原始數(shù)據(jù)進(jìn)行微分之前，需要對(duì)其進(jìn)行平滑處理，以減少噪聲的影響。

使用三階多項(xiàng)式回歸方法對(duì)原始測量數(shù)據(jù)進(jìn)行平滑處理。S-G(Savitzky Golay)濾波法是一種多項(xiàng)式回歸方法，利用給定次數(shù)的多項(xiàng)式進(jìn)行非加權(quán)線性最小二乘擬合。

使用多項(xiàng)式擬合一組數(shù)量為2M+1，且以n=0為中心的數(shù)列x(n)，即

(5)

式(5)中：f(n)為目標(biāo)函數(shù)，n為0～+∞的隨機(jī)變量；ak為待定系數(shù)，k=1,2,…,N。

最小二乘擬合的殘差(ε)為

(6)

利用卷積運(yùn)算，

(7)

式(7)中：h[m]、x[m]均為卷積變量。

若要ε最小，ε對(duì)各個(gè)參數(shù)的偏導(dǎo)數(shù)都應(yīng)為0，即

(8)

化簡后得

(9)

對(duì)用于平滑的三階多項(xiàng)式擬合進(jìn)行兩次微分運(yùn)算，計(jì)算出地球坐標(biāo)系下的一階、二階導(dǎo)數(shù)。一階導(dǎo)數(shù)得到速度、角速度，二階導(dǎo)數(shù)得到加速度、角加速度。根據(jù)以往數(shù)據(jù)結(jié)論，采用四階有限差分法的微分方法對(duì)實(shí)驗(yàn)獲取的位置、姿態(tài)信息進(jìn)行微分運(yùn)算。一階導(dǎo)數(shù)的四階離散近似表達(dá)式為

x′(t)=[-x(t+2Δt)+8x(t+Δt)-8x(t-Δt)+x(t-2Δt)](12Δt)-1

(10)

式(10)中：x為變量；t為時(shí)間。

獲得結(jié)果后利用多項(xiàng)式回歸結(jié)果的導(dǎo)數(shù)來確定作用在飛機(jī)上的力和力矩，完成空氣動(dòng)力學(xué)分析機(jī)身坐標(biāo)系，如圖5所示。

圖5 機(jī)身坐標(biāo)系

作用在飛機(jī)上的總外力是空氣動(dòng)力、推力和重力的組合。

Fe=Fa+FG+FT

(11)

通過在總外力中減去重力(FG)和推力(FT)，可以得到氣動(dòng)力為

Fa=Fe-FG-FT

(12)

而重力為一常數(shù)，即

(13)

式(13)中：φ為滾轉(zhuǎn)角;θ為俯仰角。

作用在飛機(jī)上的總外力可通過將飛機(jī)的質(zhì)量乘以由數(shù)據(jù)微分給出的固定軸加速度(ax,ay,az)得

(14)

將空氣動(dòng)力(Fe)在機(jī)身坐標(biāo)系下的元素轉(zhuǎn)換為風(fēng)坐標(biāo)系，以產(chǎn)生升力(L)和阻力(D)的表達(dá)式為

L=-Fzcosα+Fxsinα

(15)

D=-Fzsinαcosβ-Fxcosβcosα-Fysinβ

(16)

根據(jù)圖5示意的幾何關(guān)系，得到迎角(α)和側(cè)滑角(β)與速度沿xb、yb、zb三軸的分量Vx、Vy、Vz的關(guān)系為

(17)

升力系數(shù)(CL)和阻力系數(shù)(CD)可表達(dá)為

(18)

式(18)中：ρ為空氣密度；V為空速；Sref為無人機(jī)基準(zhǔn)面面積。

利用轉(zhuǎn)動(dòng)運(yùn)動(dòng)方程可計(jì)算出用于小飛機(jī)上的滾轉(zhuǎn)力矩、俯仰力矩和偏航力矩分別為

(19)

滾轉(zhuǎn)力矩系數(shù)(Cl)、俯仰力矩系數(shù)(Cm)和偏航力矩系數(shù)(Cn)分別為

(20)

2.3 飛行試驗(yàn)

試驗(yàn)中，Vicon系統(tǒng)提供的軟件對(duì)貼有反射標(biāo)記的運(yùn)動(dòng)物體位置進(jìn)行三角測量，8臺(tái)Vicon紅外攝像機(jī)跟蹤圓形標(biāo)記的反射。捕捉范圍為3.61 m的矩形場地，并延伸至3 m的高度。攝像機(jī)分布在捕捉區(qū)域的上邊緣，跟蹤飛機(jī)上的6個(gè)反光標(biāo)記。反光標(biāo)記為直徑約5 mm的圓形球體。

飛行試驗(yàn)結(jié)束后，導(dǎo)出遙控器中推力輸出數(shù)值，生成試驗(yàn)推力曲線，導(dǎo)入MATLAB對(duì)空氣動(dòng)力學(xué)參數(shù)進(jìn)行分析，結(jié)果如圖6所示。

圖6(a)所示為飛機(jī)飛行的3D運(yùn)動(dòng)軌跡。圖6(b)、圖6(c)為飛機(jī)自身的速度和加速度曲線，可以看出飛機(jī)有輕微的側(cè)向運(yùn)動(dòng)，這是由于人工發(fā)射不能準(zhǔn)確按照預(yù)定軌跡飛行所致。還可得出在水平方向減速，在垂直方向上呈下降趨勢的結(jié)論，這完全符合無動(dòng)力飛機(jī)的運(yùn)動(dòng)規(guī)律。

由圖6(d)可得，飛機(jī)在迎角為0°時(shí)升力系數(shù)最小，約為0.21，在迎角為30°時(shí)取得最大升力系數(shù)1.5，升力隨著迎角的增大而增大。圖6(e)所示為飛機(jī)的阻力曲線，當(dāng)迎角為2°時(shí)阻力系數(shù)最小，約為0.12，飛機(jī)在迎角為30°時(shí)同樣取得最大阻力系數(shù)0.25，阻力隨著迎角的增大而增大。

圖6(f)所示為飛機(jī)滾轉(zhuǎn)力矩系數(shù)隨側(cè)滑角的變化曲線，在側(cè)滑角變化范圍為-8°～10°，滾轉(zhuǎn)力矩系數(shù)隨側(cè)滑角的增大而減小，飛機(jī)具有滾轉(zhuǎn)穩(wěn)定性。圖6(g)所示為飛機(jī)俯仰力矩系數(shù)隨攻角的變化曲線，在攻角變化范圍為-10°～23°俯仰力矩隨迎角的增大而減小，飛機(jī)具有縱向靜穩(wěn)定性。圖6(h)所示為飛機(jī)偏航力矩系數(shù)隨側(cè)滑角的變化曲線，在側(cè)滑角變化范圍為-10°～10°偏航力矩隨側(cè)滑角的增大而增大，飛機(jī)具有航向靜穩(wěn)定性。因此小飛機(jī)在此次實(shí)驗(yàn)過程中是靜態(tài)穩(wěn)定的。

圖6 飛行實(shí)驗(yàn)結(jié)果

3 結(jié)論

描述了一架后退式固定翼滑翔機(jī)空氣動(dòng)力學(xué)系數(shù)的獲取過程。首先對(duì)飛機(jī)在不同輸出功率下的推力做了測試及記錄，得到飛機(jī)推力曲線。飛機(jī)在動(dòng)作捕捉系統(tǒng)的記錄下進(jìn)行飛行試驗(yàn)，將遙控器記錄的輸出效率導(dǎo)出獲取飛機(jī)實(shí)時(shí)推力，將位置、姿態(tài)的基礎(chǔ)上進(jìn)行多項(xiàng)式回歸結(jié)果的導(dǎo)數(shù)來確定作用在飛機(jī)上的力和力矩，利用MATLAB建模得到飛機(jī)的氣動(dòng)曲線。提出了一種方便易行的獲取小型有動(dòng)力無人機(jī)空氣動(dòng)力學(xué)數(shù)據(jù)的方法。然而，小型無人機(jī)產(chǎn)生的推力值較小，使用電子秤測量推力可能產(chǎn)生一定的誤差，因此在未來的工作中，可嘗試對(duì)螺旋槳系數(shù)進(jìn)行研究，使用轉(zhuǎn)速儀測試螺旋槳實(shí)時(shí)轉(zhuǎn)速，根據(jù)螺旋槳性能曲線獲取發(fā)動(dòng)機(jī)輸出推力值，對(duì)比其與過往方法的誤差值，分析空氣動(dòng)力學(xué)性能。