微型撲翼飛行器動力學模型參數(shù)的靈敏度分析

畢富國，何廣平

(北方工業(yè)大學 機械與材料工程學院，北京 100144)

在1992年，美國國防高級研究計劃局(Defense Advanced Research Projects Agency,DARPA)首次提出了微型飛行器(MAV)的概念，并認為其對未來戰(zhàn)爭格局有深遠影響。在DAPRA計劃的影響下，各類微型飛行器在軍事中得到應用。而微型撲翼飛行器(FWMAV)相較于其他無人機[1]，因為其具有的高機動性、可懸停和隱蔽性好等優(yōu)點而備受關注。隨著微制造技術[2]的發(fā)展和昆蟲飛行動力學[3-4]的深入研究，目前開展微型撲翼飛行器的研究已經成為機器人學領域的一個熱點[5-7]。與傳統(tǒng)固定翼和旋翼飛行器相比。微型撲翼飛行器所具有的高氣動效率，抗干擾能力以及可懸停等特點依賴于非定常空氣動力學現(xiàn)象，如延遲失速、旋轉升力和尾跡捕捉[8]等。但由于尺度微小，微型撲翼飛行器的負載能力十分有限，系統(tǒng)設計中所需的傳感器、電源、通信設備、控制設備等機載器件對飛行器的負載能力提出了較高要求。在研制微型撲翼飛行器時，其涉及的低雷諾數(shù)非定常空氣動力學特點和性質是必須考慮的。

在小型撲翼飛行器研究方面，代表性地，如荷蘭代爾夫特理工大學的DelFly系列[9]，美國的Nano Hummingbird[10]等。目前，有關微型撲翼飛行器的研究主要以哈佛大學的RoboBee系列[11]為主，該微型撲翼飛行器由初期的80毫克升力逐步提升到253mg升力。在以往的研究中，RoboBee研究團隊主要通過實驗來定量研究機翼形態(tài)和慣性參數(shù)對撲翼飛行性能的影響。通過在不同機翼形狀和運動參數(shù)[12]下進行大量無源動態(tài)運動實驗，比較氣動力特性，逐步修正和優(yōu)化飛行器的氣動性能。本研究基于葉素理論建立氣動力模型，引入Sobol全局靈敏度分析方法，對建立的微型撲翼飛行器氣動模型中的設計參數(shù)進行定量分析，得到參數(shù)對升力的影響系數(shù)，利用數(shù)值方法研究參數(shù)與升力的關系，可以更有針對性地對主要的設計參數(shù)進行優(yōu)化，為微型撲翼飛行器的設計和控制提供參考。

1 氣動力建模

Sane[13]研究表明，基于準穩(wěn)態(tài)空氣動力學模型，根據(jù)實驗結果對氣動力系數(shù)進行修改，可以合理地預測由于延遲失速引起的瞬時氣動力，而無需使用調節(jié)撲翼飛行空氣動力學偏微分方程的復雜有限元解。因此，下面給出的空氣動力學模型是基于準穩(wěn)態(tài)方程的。準穩(wěn)態(tài)方程不能模擬尾跡捕獲這種復雜的非穩(wěn)態(tài)空氣動力學現(xiàn)象，但是已有研究表明，在類正弦運動中，尾跡捕捉所產生的氣動力對整體氣動力貢獻較小，準穩(wěn)態(tài)理論對昆蟲尺度的微型撲翼飛行器的氣動力分析是適用的。

在氣動力計算中采用葉素法來對氣動力進行計算，可以給出每個微元的氣動力計算公式[12]：

(1)

其中,U為機翼相對速度，如圖1所示；α為機翼與拍動平面的夾角；Caero(α)為關于攻角的氣動力系數(shù)。下面確定公式中的參數(shù)，從而能計算每個微元上的氣動力，如圖2所示。

圖2 參數(shù)化機翼示意圖

圖2中，R為翼根到機翼最遠點的長度；r為相對于翼根的距離；c(r)為距離翼根r處的展長分布；dr為沿展向的微元。

1.1 氣動力系數(shù)的確定

Dickinson[14]通過實驗研究總結出了關于昆蟲氣動力系數(shù)的經驗公式，并給出其應用范圍100

哈佛大學在研究過程中所采用的氣動力升力系數(shù)和阻力系數(shù)計算公式如下：

(2)

其中，CLmax=1.8，CD0=0.4，CDmax=3.4；將各參數(shù)代入式(2)可得：

CL(α)=1.8sin(2α)

CD(α)=1.9-1.5cos(2α)

并通過幾何關系得：

CN(α)=cos(α)CL+sin(α)CD

(3)

再將各個系數(shù)代入微元方程(1)中，積分后獲得單個機翅的升力、阻力及法向力。

伯克利的Deng等[15]基于實驗數(shù)據(jù)給出了如下的法向力和切向力系數(shù)，以表達考慮了延遲失速現(xiàn)象的氣動力系數(shù)：

(4)

并通過力平衡分析得到升力系數(shù)和阻力系數(shù)與法向力系數(shù)個切向力系數(shù)之間的關系：

CL(α)=cos(α)CN(α)-sin(α)CT(α)

CD(α)=sin(α)CN(α)+cos(α)CT(α)

(5)

進一步將各個系數(shù)代入微元方程(1)中，積分后獲得單個翅膀的法向力、切向力以及對應的升力和阻力。

本文通過給定攻角函數(shù)并代入兩種方法來對其準確性進行驗證，獲得結果如圖3所示。

圖3 攻角與氣動力的關系

可以看出，以上通過經驗公式得到的各組氣動力系數(shù)，當攻角在約45°時，氣動升力系數(shù)達到最大值。兩者之間的差值則主要由于實驗方法的不同導致的。

由于哈佛大學的Robobee系列已經展示了較好的飛行特性，因此，在下面均采用哈佛大學所給出的氣動力計算公式，且攻角公式由文獻[16]給出：

α(t)=αmaxsin(2πft+δ)

(6)

假設起始位置的平均攻角面偏差δ為0，αmax=45°。

由Dickinson[3]試驗獲得的平均升力系數(shù)表明，在發(fā)生超前和對稱運動時昆蟲獲得的升力較大，但在飛行器設計過程中，要實現(xiàn)超前運動存在很大困難，因此在本文中不考慮超前滯后現(xiàn)象，采用名義方程：

α(t)=αmaxsin(2πft+δ)

(7)

將式(7)代入氣動力系數(shù)計算公式(2)，取機翼向前撲動為正方向，假設機翼在行程變換時的旋轉時瞬間完成，得到圖4所示氣動力系數(shù)計算結果。

從計算結果可以看出，升、阻力系數(shù)在運動前半程和后半程先增大后減小并中間時刻出現(xiàn)最大值。由于在前半程和后半程運動變換時，機翅的方向發(fā)生改變。升力為機翅的切線方向，升力系數(shù)為0，阻力為機翅的法線方向，阻力系數(shù)不為0。

1.2 延展向方向弦長分布的確定

在使用葉素法計算升力時，需要求解沿機翼翼展向方向的弦長分布。

圖4 氣動力系數(shù)計算結果

(8)

一階半徑矩是翅膀的中心區(qū)域，取值在0.4～0.6[12]，二階半徑矩描述翅膀的面積分布。

并發(fā)現(xiàn)翅膀形狀可以用β分布來表示：

(9)

其中B(p,q)為β方程：

(10)

其中：

(11)

將無量綱方程進行推導，得到微型撲翼飛行器的翅翼輪廓方程：

(12)

1.3 機翼相對速度的確定

懸停狀態(tài)下，來流速度為0。因此，相對速度取機翅壓力中心速度[15]：

(13)

其中，R為機翼長度，ω為機翼拍動角速度。

將氣動力系數(shù)、弦長分布以及相對速度參數(shù)方程代入微元氣動力計算式(1)，得到

(14)

(15)

其中空氣密度ρ=1.205 kg/m3。

沿展向積分得到氣動力關于攻角的計算方程：

(16)

(17)

則單個翅膀的氣動力為升力和阻力的矢量和：

F=FL+FD

(18)

2 Sobol全局靈敏度分析

上一節(jié)所述氣動力模型呈高度非線性特征，且各參數(shù)之間存在耦合。傳統(tǒng)的局部靈敏度分析方法，如直接求導法、有限差分法和攝動法等是基于線性或非線性模型，對于對非線性較強的模型，局部靈敏度分析法無法提供合理的分析結果。全局靈敏度分析方法則適用于高非線性模型、非單調模型、可分析全范圍參數(shù)，并可以考察參數(shù)間相互作用，避免了局部分析方法的局限性。因此，在氣動力參數(shù)靈敏度分析過程中考慮其特點采用全局靈敏度分析方法，具有重要的實際應用價值。其中，Sobol方法是研究全局靈敏度常用的一種方法[18～20]。

2.1 Sobol分析方法理論

Sobol方法是一種基于方差分解的全局靈敏度定量分析方法。氣動力模型Y=f(X)，X=x1，x2，x3,…,xk,其中x1，x2，x3,…,xk為模型中待分析的輸入?yún)?shù)。Y為模型輸出值。可以將模型分解為如下形式：

f1,2,…,k(x1,x2,…,xk))

(19)

其中：

f0=E(Y)

fi(xi)=E(Y|xi)-f0

fij(xi，xj)=E(Y|xi，xj)-f0-fi-fj

從其中可以看出，fi是單獨改變xi的效果(稱為xi的主要效應)，并且fij表示除了它們各自變化的效果之外，同時改變xi和xj對輸出的效果，這被稱為二階交互。高階項具有類似的定義。

如果輸入?yún)?shù)之間是相互獨立的，則方程式(19)右側各項相互正交，協(xié)方差為0。對兩側求方差得：

(20)

其中：

Vi=Varxi(Ex～i(Y|xi))

Vij=Varx～ij(Ex～ij(Y|xi,xj))-Vi-Vj

x～i表示除了xi之外的所有變量的集合。

通過計算比率：

(21)

式中，Si是xi的一階靈敏度系數(shù)，表示參數(shù)xi對模型輸出的影響。

同理定義比率：

(22)

式中，STi是xi的全局靈敏度系數(shù)，表示考慮參數(shù)耦合作用后，xi對模型輸出的影響。

2.2 Sobol法計算過程

首先，進行因子篩選。根據(jù)定義設定展長為15 mm，對氣動力模型進行數(shù)值代入計算得到運動參數(shù)對氣動力的影響，如圖5所示。

圖5 運動參數(shù)對氣動力的影響

由圖5可以看出，在其他參數(shù)確定的情況下，氣動力隨最大撲動角的增大而增大；隨攻角的增大先增大后減小；隨撲動頻率增大而增大。在上述分析中只是分別對每個參數(shù)對氣動力的影響進行了控制變量分析，并未考慮在各個設計參數(shù)共同影響下氣動力的變化。在確定分析參數(shù)之后，對氣動力模型分析參數(shù)進行兩次采樣，得到AN·D和BN·D矩陣，其中N為抽樣參數(shù)，D為分析參數(shù)。另構造ABi(i=1,2，…，D)即用B的第i列替換A的第i列。

(23)

V(Y)=V(A)+V(B)

(24)

其中，根據(jù)文獻[21]圖6及表2給出的數(shù)據(jù)表明，在參數(shù)量為5的模型全局靈敏度分析中，取樣數(shù)在500左右即可達到收斂。在本次全局靈敏度分析過程中取樣數(shù)為1 000，滿足取樣要求。參數(shù)取值范圍如表1所示。

表1 各參數(shù)取值范圍

3 計算結果及分析

考慮到飛行器設計過程中機翼幾何參數(shù)的加工難度和加工條件要求較高，且在飛行器實際飛行過程中主要依賴于運動參數(shù)進行控制，因此在以往研究過程中關于參數(shù)靈敏度分析的研究受到更多關注。本文同時進行了幾何參數(shù)的靈敏度分析，從而為機翼設計提供參考。

通過Matlab計算得到氣動力模型中機翼幾何參數(shù)一階半徑矩、展弦比以及運動參數(shù)最大撲動角、最大攻角、撲動頻率的全參數(shù)一階靈敏度和全局靈敏度，計算結果分別是：總方差為4.6 e-2；一階半徑距方差為2.42 e-4；展弦比方差為1.19e-2；最大撲動角方差為1.579e-3； 最大攻角方差為1.576e-3；撲動頻率方差為1.7e-2。

計算可得：一階半徑矩的一階靈敏度為5.26e-3，全局靈敏度為2.43e-2；展弦比的一階靈敏度為2.59e-1，全局靈敏度為5.28e-1；最大撲動角的一階靈敏度為3.42e-2，全局靈敏度為1.41e-1；最大攻角的一階靈敏度為3.43e-2，全局靈敏度為1.08e-1；撲動頻率的一階靈敏度為3.68e-1，全局靈敏度為5.88e-1。分析結果如圖6所示。

圖6 全參數(shù)靈敏度分析結果

進一步，按照上述分析方法分別對機翼設計中設計的機翼幾何參數(shù)和運動控制參數(shù)進行靈敏度分析。

在獨立的幾何參數(shù)分析中設定最大撲動角為120°；最大攻角為45°；撲動頻率為120。得到分析結果為一階半徑矩的一階靈敏度為3.65e-2，全局靈敏度為4.48e-2；展弦比的一階靈敏度為9.95e-1，全局靈敏度為1.01。結果如圖7所示。

圖7 幾何參數(shù)靈敏度分析結果

在獨立的運動參數(shù)分析中設定一階半徑矩為0.56；展弦比為7。得到分析結果為最大撲動角的一階靈敏度為1.26e-1，全局靈敏度為2.32e-1；最大攻角的一階靈敏度為6.95e-2，全局靈敏度為1.09e-1；撲動頻率的一階靈敏度為5.82e-1，全局靈敏度為7.46e-1。結果如圖8所示。

圖8 運動參數(shù)靈敏度分析結果

通過對比可以看出，全參數(shù)分析結果與對幾何參數(shù)和運動參數(shù)分別分析得出的結果保持一致，從而證明了分析結果的有效性。全參數(shù)分析顯示，對于氣動力影響的因素大小依次為：撲動頻率>展弦比>最大撲動角>最大攻角>一階半徑矩。運動參數(shù)全局靈敏度為撲動頻率>最大撲動角>最大攻角。

相比于全局靈敏度，一階靈敏度的排序與全局靈敏度基本一致，同時也可以看出參數(shù)之間存在影響。在進行參數(shù)優(yōu)化設計過程中，全局靈敏度較小的參數(shù)可以選取固定值，這樣既可以減小在優(yōu)化時的計算成本，也為控制律設計提供了方便。

4 結論

建立了基于葉素法的微型撲翼飛行器的氣動力學模型，并采用Sobol全局靈敏度分析方法對氣動力模型進行靈敏度分析，獲得了機翼幾何參數(shù)和運動參數(shù)對氣動力的影響靈敏度。分析表明，幾何參數(shù)中展弦比的全局靈敏度最大，為5.28e-1，一階半徑距全局靈敏度最小，為2.44e-2；運動參數(shù)中撲動頻率的全局靈敏度最大，為5.88e-2，撲動角的全局靈敏度為1.41e-2，攻角的全局靈敏度為1.08e-1。可以看出，要提高微型撲翼飛行器的空氣動力性能，在結構參數(shù)方面，翅翼的展弦比需要進行精心設計，在運動參數(shù)方面，撲翼飛行器的攻角對撲翼飛行器的升力影響最大，在控制器設計中應優(yōu)先進行穩(wěn)定控制。