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

畢富國(guó)，何廣平

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

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

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

1 氣動(dòng)力建模

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

在氣動(dòng)力計(jì)算中采用葉素法來(lái)對(duì)氣動(dòng)力進(jìn)行計(jì)算，可以給出每個(gè)微元的氣動(dòng)力計(jì)算公式[12]：

(1)

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

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

圖2中，R為翼根到機(jī)翼最遠(yuǎn)點(diǎn)的長(zhǎng)度；r為相對(duì)于翼根的距離；c(r)為距離翼根r處的展長(zhǎng)分布；dr為沿展向的微元。

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

Dickinson[14]通過(guò)實(shí)驗(yàn)研究總結(jié)出了關(guān)于昆蟲(chóng)氣動(dòng)力系數(shù)的經(jīng)驗(yàn)公式，并給出其應(yīng)用范圍100

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

(2)

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

CL(α)=1.8sin(2α)

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

并通過(guò)幾何關(guān)系得：

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

(3)

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

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

(4)

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

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

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

(5)

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

本文通過(guò)給定攻角函數(shù)并代入兩種方法來(lái)對(duì)其準(zhǔn)確性進(jìn)行驗(yàn)證，獲得結(jié)果如圖3所示。

圖3 攻角與氣動(dòng)力的關(guān)系

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

由于哈佛大學(xué)的Robobee系列已經(jīng)展示了較好的飛行特性，因此，在下面均采用哈佛大學(xué)所給出的氣動(dòng)力計(jì)算公式，且攻角公式由文獻(xiàn)[16]給出：

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

(6)

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

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

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

(7)

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

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

1.2 延展向方向弦長(zhǎng)分布的確定

在使用葉素法計(jì)算升力時(shí)，需要求解沿機(jī)翼翼展向方向的弦長(zhǎng)分布。

圖4 氣動(dòng)力系數(shù)計(jì)算結(jié)果

(8)

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

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

(9)

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

(10)

其中：

(11)

將無(wú)量綱方程進(jìn)行推導(dǎo)，得到微型撲翼飛行器的翅翼輪廓方程：

(12)

1.3 機(jī)翼相對(duì)速度的確定

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

(13)

其中，R為機(jī)翼長(zhǎng)度，ω為機(jī)翼拍動(dòng)角速度。

將氣動(dòng)力系數(shù)、弦長(zhǎng)分布以及相對(duì)速度參數(shù)方程代入微元?dú)鈩?dòng)力計(jì)算式(1)，得到

(14)

(15)

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

沿展向積分得到氣動(dòng)力關(guān)于攻角的計(jì)算方程：

(16)

(17)

則單個(gè)翅膀的氣動(dòng)力為升力和阻力的矢量和：

F=FL+FD

(18)

2 Sobol全局靈敏度分析

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

2.1 Sobol分析方法理論

Sobol方法是一種基于方差分解的全局靈敏度定量分析方法。氣動(dòng)力模型Y=f(X)，X=x1，x2，x3,…,xk,其中x1，x2，x3,…,xk為模型中待分析的輸入?yún)?shù)。Y為模型輸出值?？梢詫⒛Ｐ头纸鉃槿缦滦问剑?/p>

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是單獨(dú)改變xi的效果(稱(chēng)為xi的主要效應(yīng))，并且fij表示除了它們各自變化的效果之外，同時(shí)改變xi和xj對(duì)輸出的效果，這被稱(chēng)為二階交互。高階項(xiàng)具有類(lèi)似的定義。

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

(20)

其中：

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

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

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

通過(guò)計(jì)算比率：

(21)

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

同理定義比率：

(22)

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

2.2 Sobol法計(jì)算過(guò)程

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

圖5 運(yùn)動(dòng)參數(shù)對(duì)氣動(dòng)力的影響

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

(23)

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

(24)

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

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

3 計(jì)算結(jié)果及分析

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

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

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

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

進(jìn)一步，按照上述分析方法分別對(duì)機(jī)翼設(shè)計(jì)中設(shè)計(jì)的機(jī)翼幾何參數(shù)和運(yùn)動(dòng)控制參數(shù)進(jìn)行靈敏度分析。

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

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

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

圖8 運(yùn)動(dòng)參數(shù)靈敏度分析結(jié)果

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

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

4 結(jié)論

建立了基于葉素法的微型撲翼飛行器的氣動(dòng)力學(xué)模型，并采用Sobol全局靈敏度分析方法對(duì)氣動(dòng)力模型進(jìn)行靈敏度分析，獲得了機(jī)翼幾何參數(shù)和運(yùn)動(dòng)參數(shù)對(duì)氣動(dòng)力的影響靈敏度。分析表明，幾何參數(shù)中展弦比的全局靈敏度最大，為5.28e-1，一階半徑距全局靈敏度最小，為2.44e-2；運(yùn)動(dòng)參數(shù)中撲動(dòng)頻率的全局靈敏度最大，為5.88e-2，撲動(dòng)角的全局靈敏度為1.41e-2，攻角的全局靈敏度為1.08e-1?？梢钥闯觯岣呶⑿蛽湟盹w行器的空氣動(dòng)力性能，在結(jié)構(gòu)參數(shù)方面，翅翼的展弦比需要進(jìn)行精心設(shè)計(jì)，在運(yùn)動(dòng)參數(shù)方面，撲翼飛行器的攻角對(duì)撲翼飛行器的升力影響最大，在控制器設(shè)計(jì)中應(yīng)優(yōu)先進(jìn)行穩(wěn)定控制。