吳翰 王正平 周 洲 王睿

摘 要：

由于垂起固定翼無(wú)人機(jī)結(jié)構(gòu)復(fù)雜， 無(wú)論是采用數(shù)值模擬得到其非定常氣動(dòng)力還是采用單剛體建模方法建立其動(dòng)力學(xué)模型都較為困難。 為了解決該問(wèn)題， 通過(guò)CFD數(shù)值模擬機(jī)翼翼型的遲滯環(huán)線， 引入無(wú)人機(jī)俯仰運(yùn)動(dòng)所產(chǎn)生的遲滯效應(yīng)， 通過(guò)ONERA方程建立無(wú)人機(jī)垂起改平飛過(guò)程的非定常氣動(dòng)力模型， 最終基于多體動(dòng)力學(xué)， 將無(wú)人機(jī)劃分成機(jī)翼、 機(jī)身、 旋翼、 舵面的多剛體系統(tǒng)， 通過(guò)凱恩方程推導(dǎo)并建立無(wú)人機(jī)動(dòng)力學(xué)模型， 通過(guò)數(shù)值仿真可以發(fā)現(xiàn)： 多體動(dòng)力學(xué)的方法適用于這類結(jié)構(gòu)復(fù)雜無(wú)人機(jī)的建模， 引入機(jī)翼的遲滯環(huán)對(duì)無(wú)人機(jī)的俯仰運(yùn)動(dòng)起阻尼作用; 引入非定常氣動(dòng)力與未引入非定常氣動(dòng)力相比垂起速度存在8%的差異， 垂起改平飛過(guò)程結(jié)束時(shí)兩者上升速度存在40%的差異。

關(guān)鍵詞：???? ??垂起固定翼無(wú)人機(jī); 遲滯環(huán); ONERA方程; 多體動(dòng)力學(xué); 凱恩方程

中圖分類號(hào)：??? ??TJ011;? V212 ??文獻(xiàn)標(biāo)識(shí)碼：??? ?A ??文章編號(hào)：??? ??1673-5048（2019）02-0021-08

0 引? 言

垂起固定翼無(wú)人機(jī)一般由旋翼和固定翼所組成， 兼具垂直起降和高巡航效率的性能， 具有很大的研究?jī)r(jià)值和市場(chǎng)前景， 然而對(duì)于這類無(wú)人機(jī)的動(dòng)力學(xué)建模仍然是個(gè)難點(diǎn)。 文獻(xiàn)[1-2]分別采用牛頓-歐拉方程建立旋翼無(wú)人機(jī)動(dòng)力學(xué)模型， 文獻(xiàn)[3-4]分別采用拉格朗日方程[5]建立固定翼無(wú)人機(jī)動(dòng)力學(xué)模型， 文獻(xiàn)[6]采用牛頓-歐拉方程建立傾轉(zhuǎn)旋翼無(wú)人機(jī)的動(dòng)力學(xué)模型， 以上文獻(xiàn)是將無(wú)人機(jī)的氣動(dòng)力看成準(zhǔn)定常模型， 未系統(tǒng)地考慮無(wú)人機(jī)的非定常氣動(dòng)力模型。 本文主要對(duì)垂起固定翼無(wú)人機(jī)的動(dòng)力學(xué)建模與非定常氣動(dòng)力進(jìn)行了研究， 針對(duì)垂起固定翼這類特殊布局的無(wú)人機(jī)， 提出三個(gè)值得研究的問(wèn)題： （1）現(xiàn)有單剛體建模方法， 在建立復(fù)雜結(jié)構(gòu)無(wú)人機(jī)的動(dòng)力學(xué)模型時(shí)較為繁瑣， 能否采用新的建模方法， 以便更合理、 更快速地建立這類結(jié)構(gòu)復(fù)雜的無(wú)人機(jī)的動(dòng)力學(xué)模型; （2）以前的動(dòng)力學(xué)模型都是將無(wú)人機(jī)看成準(zhǔn)定常模型進(jìn)行處理， 在實(shí)際情況中， 由于垂起固定翼無(wú)人機(jī)的機(jī)翼面積較大， 機(jī)翼的遲滯效應(yīng)將對(duì)無(wú)人機(jī)的飛行動(dòng)力學(xué)產(chǎn)生較大影響， 因此， 如何將遲滯效應(yīng)與無(wú)人機(jī)的動(dòng)力學(xué)模型耦合起來(lái)， 值得研究; （3）垂起固定翼無(wú)人機(jī)最具有研究?jī)r(jià)值的環(huán)節(jié)為垂起改平飛的過(guò)渡過(guò)程。 在該過(guò)程中無(wú)人機(jī)的迎角變化幅度較大， 原有氣動(dòng)模型不再適用， 而如果直接采用CFD進(jìn)行該階段的數(shù)值模擬， 則運(yùn)算量太大不適用于工程領(lǐng)域， 因此， 能否采用更為便利的方法建立這一過(guò)程無(wú)人機(jī)的非定常氣動(dòng)力模型。

為了解決上述的三個(gè)問(wèn)題， 對(duì)比現(xiàn)有的動(dòng)力學(xué)建模方法， 選取凱恩方程[7]基于多體動(dòng)力學(xué)[8] 建立該垂起固定翼無(wú)人機(jī)的動(dòng)力學(xué)模型; ?采用CFD數(shù)值模擬[9]該無(wú)人機(jī)機(jī)翼力矩系數(shù)的遲滯環(huán)線， 采用積分法得到無(wú)人機(jī)縱向動(dòng)導(dǎo)數(shù)， 通過(guò)縱向動(dòng)導(dǎo)數(shù)引入無(wú)人機(jī)機(jī)翼的遲滯效應(yīng); ?采用ONERA

非線性氣動(dòng)力方程建立無(wú)人機(jī)垂起改平飛過(guò)渡過(guò)程中的非定常氣動(dòng)力模型， 完善垂起改平飛過(guò)渡階段的動(dòng)力學(xué)模型。

1 垂起固定翼無(wú)人機(jī)多剛體劃分

該垂起固定翼無(wú)人機(jī)采用常規(guī)布局， “T”形尾， 其主要由四個(gè)垂起旋翼和兩個(gè)前拉螺旋槳組成。 按照多體動(dòng)力學(xué)的思路[10]將其劃分為左右機(jī)翼、 機(jī)身、 旋翼、 螺旋槳、 垂尾、 左右平尾、 電機(jī)等18個(gè)剛體， 其中機(jī)翼、 垂尾、 電機(jī)與機(jī)身固結(jié)， 平尾與垂尾固結(jié)， 各舵面與翼面鉸接。 忽略無(wú)人機(jī)在空中飛行時(shí)的柔性變形， 其廣義坐標(biāo)為 q =[ x y z φ θ ψ ]T， 表示無(wú)人機(jī)質(zhì)心處的位置和姿態(tài)， 以各剛體質(zhì)心為原點(diǎn)建立右手坐標(biāo)系， 具體坐標(biāo)系建立如圖1所示。

2 機(jī)翼遲滯效應(yīng)分析

2.1 無(wú)人機(jī)翼型遲滯環(huán)模擬

當(dāng)無(wú)人機(jī)在做俯仰運(yùn)動(dòng)時(shí)， 機(jī)翼所產(chǎn)生的氣動(dòng)力為非定常氣動(dòng)力， 其控制方程為描述氣體非定常運(yùn)動(dòng)的NS方程[11]：

由于該無(wú)人機(jī)展弦比大且是平直機(jī)翼， 因此， 可以用機(jī)翼二維翼型的數(shù)值模擬結(jié)果來(lái)代替整個(gè)三維機(jī)翼， 該機(jī)翼翼型為FX-60-126-1， 劃分網(wǎng)格如圖2所示。 強(qiáng)迫其繞剛心做0°到15°的俯仰運(yùn)動(dòng)， 得到不同飛行速度下力矩系數(shù)的遲滯環(huán)線， 如圖3所示。

2.2 無(wú)人機(jī)機(jī)翼動(dòng)導(dǎo)數(shù)求解

以遲滯環(huán)為基準(zhǔn)， 推導(dǎo)無(wú)人機(jī)機(jī)翼俯仰運(yùn)動(dòng)動(dòng)導(dǎo)數(shù)， 動(dòng)導(dǎo)數(shù)是氣動(dòng)載荷對(duì)狀態(tài)變化率（速度大小）的導(dǎo)數(shù)， 是描述無(wú)人機(jī)動(dòng)態(tài)變化過(guò)程的重要參數(shù)， 應(yīng)用積分法求取無(wú)人機(jī)動(dòng)導(dǎo)數(shù)。 假設(shè)機(jī)翼在俯仰自由度下進(jìn)行小振幅強(qiáng)迫振蕩運(yùn)動(dòng)， 其迎角為

式中： k=（ωb/v）， 為振蕩運(yùn)動(dòng)的減縮頻率;? t為無(wú)量綱時(shí)間; θ=Δα=α-α0為機(jī)翼俯仰角。

式中： Cmθ0為俯仰剛度;? Cmθ? · 0為俯仰阻尼導(dǎo)數(shù)， 利用上文得到的周期性簡(jiǎn)諧強(qiáng)迫運(yùn)動(dòng)的氣動(dòng)力矩， 通過(guò)數(shù)值辨識(shí)即可直接獲得俯仰阻尼導(dǎo)數(shù)。 當(dāng)強(qiáng)迫運(yùn)動(dòng)為周期性正弦運(yùn)動(dòng)時(shí)， 對(duì)式（6）兩邊同時(shí)乘以sin（kt）和cos（kt）， 再進(jìn)行積分可得

式中： ts為簡(jiǎn)諧強(qiáng)迫運(yùn)動(dòng)達(dá)到周期性解后的任意時(shí)刻; T為簡(jiǎn)諧強(qiáng)迫運(yùn)動(dòng)周期; Cm（t）的數(shù)據(jù)由CFD算得;? Cmθ， Cmθ? · 分別為無(wú)人機(jī)俯仰運(yùn)動(dòng)所產(chǎn)生的動(dòng)導(dǎo)數(shù)。

3 ONERA非線性氣動(dòng)力方程分析

ONERA方程[13]目前主要被用于解決機(jī)翼的氣動(dòng)彈性問(wèn)題。 該方程為非線性氣動(dòng)力方程， 適用于大迎角的情形。 ONERA方程的線性部分是對(duì)Theodorsen氣動(dòng)力的模擬， 其非線性部分則考慮了靜態(tài)失速的影響。 通過(guò)該方程的推導(dǎo)以及推導(dǎo)的前提[14]可以發(fā)現(xiàn)該方程適用于無(wú)人機(jī)垂起改平飛的過(guò)渡過(guò)程， 選取ONERA方程建立該過(guò)程的非定常氣動(dòng)力模型。

首先給出升力系數(shù)為

CL=CLa+CLbCLa=Sz1×α? · +Sz2×α? ¨ +Sz3×α? · +CLγC? · Lb+rz1×C? · Lb+rz2×CLb=-rz2×ΔCL-

rz3× ΔCL α ×α? ·? ? （10）

式中：α為瞬時(shí)迎角; CLa， CLb為升力系數(shù)的組成部分; ΔCL為ONERA方程線性部分與非線性部分的差值; CLγ為非線性部分的升力系數(shù)。 式中rz1， rz2， rz3均與雷諾數(shù)有關(guān)， 需要通過(guò)雷諾數(shù)進(jìn)行確定。 由于ONERA方程的線性部分是對(duì)經(jīng)典Theodorsen氣動(dòng)力模型的擬合， 因此可以得到如下參數(shù)值：

sz1=π， sz2=π/2， sz3=0 （11）

同樣可以得到ONERA方程線性部分與非線性部分的差值ΔCL為

ΔCL=

0 0≤α≤0.139 6 rad6.322 84×（α-0.139 6） ?0.139 6 rad<

α≤0.314 2 rad6.322 84×（α-0.139 6）-0.422 84×

（α-0.314 2）? 0.314 2 rad<α ? ΔCL= 0 -0.139 6 rad≤α≤06.322 84×（α-0.139 6） ?-0.314 2 rad≤

α<-0.139 6 rad6.322 84×（α-0.139 6）-0.422 84×

（α-0.314 2）? ?α<-0.314 2 rad

（12）

通過(guò)上式可以得到 ΔCL α ， 將等式（12）和 ΔCL α 帶入式（10）中便可以得到相應(yīng)情況下的升力系數(shù)。

ONERA方程的阻力系數(shù)為

CD=CDa+CDbC? ¨ Db+rD1×C? · Db=-r2D×ΔCD-rD3×α? · ΔCD=-0.042×α-0.147 3×α2-4.923×α3 ? （13）

式中： CDa， CDb為阻力系數(shù)的組成部分， ΔCD為線性部分阻力系數(shù)與非線性部分阻力系數(shù)的差值。 rD1， rD2， rD3同樣由雷諾數(shù)確定。 由于無(wú)人機(jī)在過(guò)渡階段對(duì)其側(cè)向力系數(shù)Cyβ影響較小， 因此側(cè)向力系數(shù)采用準(zhǔn)定常假設(shè)進(jìn)行求解即可。

4 機(jī)翼非定常動(dòng)力學(xué)模型建立

4.1 機(jī)翼非定常氣動(dòng)力模型

以左機(jī)翼為例基于上文得到的動(dòng)導(dǎo)數(shù)和ONERA方程建立其非定常氣動(dòng)力模型為

f z= C zl - 1 2 ρV2SCD（Re， α， θ）- 1 2 ρV2SCyβ（Re， α， θ）- 1 2 ρV2SCL（Re， α， θ）（14）

式中： ?C zl為氣流向左機(jī)翼的坐標(biāo)系轉(zhuǎn)換矩陣; V為左機(jī)翼質(zhì)心處的合速度， CD， Cyβ， CL分別為左機(jī)翼的氣動(dòng)力系數(shù)。

左機(jī)翼質(zhì)心處的力模型為

F z= f z+? 00mzg?（15）

式中： ?F z為左機(jī)翼的力模型; mz為左機(jī)翼的質(zhì)量; ?f z為左機(jī)翼的非定常氣動(dòng)力模型。

4.2 機(jī)翼非定常力矩模型

垂起固定翼無(wú)人機(jī)的力矩主要由兩部分產(chǎn)生， 一部分是各剛體質(zhì)心處的廣義主動(dòng)力相對(duì)于無(wú)人機(jī)質(zhì)心所產(chǎn)生的力矩， 該部分力矩將通過(guò)凱恩方程的偏速度矩陣進(jìn)行引入。 另一部分是由于無(wú)人機(jī)在飛行過(guò)程中， 滾轉(zhuǎn)、 偏航以及俯仰運(yùn)動(dòng)所引起的非定常氣動(dòng)力矩。 直接給出左機(jī)翼非定常力矩模型如下：

M z=?? Clpp +Clrr?? 1 2 ρV 2Sc

（Cmθθ+Cmθ? · q? - ） 1 2 ρV 2Sl Cnpp +Cnrr?? 1 2 ρV 2Sc?（16）

式中： p? - = pc 2V ， q? - = qc 2V ， r = rc 2V ， Clp， Clr， Cnp， Cnr為橫航向動(dòng)導(dǎo)數(shù)， 表示該垂起固定翼無(wú)人機(jī)滾轉(zhuǎn)和偏航時(shí)左機(jī)翼所產(chǎn)生的非定常氣動(dòng)力矩;? Cmθ， Cmθ? · 為無(wú)人機(jī)做俯仰運(yùn)動(dòng)時(shí)所產(chǎn)生非定常氣動(dòng)力矩系數(shù);? S為機(jī)翼面積;? l， c分別為左機(jī)翼展長(zhǎng)和弦長(zhǎng);? p， q， r為機(jī)翼的角速度分量。 需要特別注意的是主要通過(guò)無(wú)人機(jī)縱向動(dòng)導(dǎo)數(shù)Cmθ， Cmθ? · 引入無(wú)人機(jī)俯仰運(yùn)動(dòng)中機(jī)翼所產(chǎn)生的遲滯效應(yīng)， 而且遲滯效應(yīng)僅在無(wú)人機(jī)迎角小于15°時(shí)才具有。

5 垂起固定翼無(wú)人機(jī)動(dòng)力學(xué)模型建立

多體動(dòng)力學(xué)建模方法主要應(yīng)用于機(jī)器人[15]等結(jié)構(gòu)較為復(fù)雜的對(duì)象， 這種建模方法是將建模對(duì)象劃分為多剛體系統(tǒng)， 針對(duì)每個(gè)剛體建立其動(dòng)力學(xué)模型， 然后通過(guò)方程自身的約束條件， 將各剛體的動(dòng)力學(xué)模型引入整個(gè)對(duì)象的質(zhì)心， 進(jìn)而完成動(dòng)力學(xué)建模。 對(duì)于該垂起固定翼無(wú)人機(jī)， 這種建模方法較為適用。

基于多體動(dòng)力學(xué)中的凱恩方程建立該垂起固定翼無(wú)人機(jī)動(dòng)力學(xué)模型， 凱恩方程所建模型為一階偏微分方程組， 可減少仿真計(jì)算量， 其建模的基準(zhǔn)為系統(tǒng)廣義主動(dòng)力與系統(tǒng)廣義慣性力之和為零， 系統(tǒng)廣義主動(dòng)力矩與系統(tǒng)廣義慣性力矩之和為零。

首先建立無(wú)人機(jī)各剛體的廣義慣性力和廣義慣性力矩模型， 依據(jù)無(wú)人機(jī)各剛體質(zhì)心的線加速度和角加速度， 建立各剛體廣義慣性力和廣義慣性力矩模型如下：

R i=- m i a i， ??R *i=- J i ω?? · i （17）

式中： ?m i， ?J i分別為剛體i的質(zhì)量矩陣和轉(zhuǎn)動(dòng)慣量矩陣;? ?a i， ?ω?? · i分別為剛體i質(zhì)心處的線加速度和角加速度矩陣;? ?R i， ?R *i分別為剛體i質(zhì)量所產(chǎn)生的廣義慣性力和轉(zhuǎn)動(dòng)慣量所產(chǎn)生的廣義慣性力矩。

無(wú)人機(jī)整機(jī)的廣義慣性力和廣義慣性力矩模型即為各剛體廣義慣性力和廣義慣性力矩模型的疊加， 具體形式如下[6]：

K =∑ i??V iq?? ·?? R i+ω iq?? ·?? R *i（18）

式中： i代表剛體的編號(hào)， 由于無(wú)人機(jī)被劃分為18個(gè)剛體， 因此該方程為18個(gè)子方程的疊加;? V iq?? ·? 和ω iq?? ·? 分別為剛體i的偏線速度矩陣和偏角速度矩陣。

無(wú)人機(jī)的廣義主動(dòng)力為各剛體質(zhì)心力模型的疊加， 無(wú)人機(jī)系統(tǒng)的廣義主動(dòng)力矩由兩部分組成， 一部分是各剛體質(zhì)心力模型相對(duì)于全機(jī)質(zhì)心所產(chǎn)生的力矩; 另一部分是各剛體自身運(yùn)動(dòng)時(shí)繞其自身質(zhì)心所產(chǎn)生的力矩。 采用左機(jī)翼動(dòng)力學(xué)模型， 建立無(wú)人機(jī)右機(jī)翼、 機(jī)身、 左右平尾、 垂尾的力模型和非定常力矩模型， 分別為（ F y， ?F b， ?F p， ?F q， ?F r）和 （M y， ?M b， ?M p， ?M q， ?M r）; 對(duì)于無(wú)人機(jī)的前拉螺旋槳和垂起旋翼， 采用文獻(xiàn)[1]中的方法通過(guò)流入比和前進(jìn)比建立其力模型和力矩模型， 分別為（ F c， ?F f， ?F a， ?F n， ?F d， ?F g）和（ M c， ?M f， ?M a ?M n， ?M d， ?M g）; 對(duì)于無(wú)人機(jī)的電機(jī)采用文獻(xiàn)[16]中的方法建立其力和力矩模型， 分別為（ F 1， ?F 2， ?F 3， ?F 4， ?F 5， ?F 6）和（ M 1， ?M 2， ?M 3， ?M 4， ?M 5， ?M 6）。 將上述的力和力矩模型以及左機(jī)翼模型進(jìn)行疊加， 可得到該垂起固定翼無(wú)人機(jī)的廣義主動(dòng)力和廣義主動(dòng)力矩模型， 其具體形式如下[17]：

K *=∑? i???V iq?? ·?? （ F i）+?ω iq?? ·?? （ M i）（19）

式中： i表示各個(gè)剛體的編號(hào); ?F i中包含剛體i自身的重力;? ?M i為剛體i在機(jī)體坐標(biāo)系下的合力矩， 其包含剛體i質(zhì)心氣動(dòng)力所產(chǎn)生的力矩以及剛體i的自身運(yùn)動(dòng)所產(chǎn)生的非定常力矩。 各剛體質(zhì)心氣動(dòng)力和各剛體重力相對(duì)于無(wú)人機(jī)質(zhì)心所產(chǎn)生的力矩， 由偏線速度矩陣與力 F i的乘積引入。

將式（18）和式（19）疊加便可得到該無(wú)人機(jī)的動(dòng)力學(xué)模型[18]為

K + K *=∑? i???V iq?? ·?? （ F i+ R i）+?ω iq?? ·?? （ M i+ R *i） =0? （20）

式（20）即為該垂起無(wú)人機(jī)的六自由度動(dòng)力學(xué)模型， 其中?V iq?? ·?? ， ?? ω iq?? ·?? 均為6×3的矩陣，F(xiàn) i， ?M i， ?R i， ?R *i均為3×1矩陣， 因此通過(guò)式（20）可以得到6個(gè)子方程， 再加上位移與速度和姿態(tài)角與角速度關(guān)系式， 共計(jì)12個(gè)未知量12個(gè)方程， 可直接采用Matlab進(jìn)行數(shù)值仿真。

6 數(shù)值仿真

6.1 俯仰運(yùn)動(dòng)仿真

為了研究無(wú)人機(jī)機(jī)翼的遲滯環(huán)對(duì)無(wú)人機(jī)運(yùn)動(dòng)的影響， 對(duì)無(wú)人機(jī)小迎角俯仰運(yùn)動(dòng)進(jìn)行仿真。 該仿真初始條件為： 垂起固定翼無(wú)人機(jī)飛行高度120 m， ?前飛速度為10 m/s， 四個(gè)垂起旋翼轉(zhuǎn)速均為2 000 r/min， 前拉螺旋槳轉(zhuǎn)速為1 500 r/min， 初始迎角為4°， 其余狀態(tài)量初始時(shí)刻均為0。 有遲滯環(huán)修正的模型仿真曲線對(duì)應(yīng)有遲滯環(huán)， 無(wú)遲滯環(huán)修正的對(duì)應(yīng)無(wú)遲滯環(huán)， 在該飛行過(guò)程當(dāng)中無(wú)橫航向運(yùn)動(dòng)。

6.2 垂起改平飛過(guò)渡階段仿真

對(duì)于垂起固定翼無(wú)人機(jī)而言， 其垂起改平飛的過(guò)渡階段是數(shù)值仿真和研究的難點(diǎn)。 采用非定常氣動(dòng)力模型和未引入非定常氣動(dòng)力（即迎角小于15°時(shí)為定常氣動(dòng)模型， 迎角大于15°時(shí)， 機(jī)翼不產(chǎn)生氣動(dòng)力）的動(dòng)力學(xué)模型進(jìn)行對(duì)比。 仿真過(guò)程為， 在前40 s無(wú)人機(jī)進(jìn)行垂直起飛， 四個(gè)垂起旋翼轉(zhuǎn)速均為2 000 r/min， 并且保持不變， 40～45 s， 無(wú)人機(jī)四個(gè)垂起旋翼轉(zhuǎn)速逐漸減為零， 兩個(gè)前拉螺旋槳轉(zhuǎn)速均由零增加為2 000 r/min， 45 s后， 四個(gè)垂起旋翼轉(zhuǎn)速為零， 兩個(gè)前拉螺旋槳轉(zhuǎn)速保持不變， 直到60 s無(wú)人機(jī)前飛速度達(dá)到15 m/s， 無(wú)人機(jī)完成垂起改平飛的過(guò)渡過(guò)程。

6.3 盤旋運(yùn)動(dòng)仿真

無(wú)人機(jī)的盤旋運(yùn)動(dòng)是一個(gè)橫航向與縱向耦合的運(yùn)動(dòng)， 無(wú)人機(jī)在空中執(zhí)行任務(wù)時(shí)， 常會(huì)發(fā)生這樣的耦合運(yùn)動(dòng)。 對(duì)垂起固定翼無(wú)人機(jī)的盤旋運(yùn)動(dòng)進(jìn)行數(shù)值仿真以研究非定常氣動(dòng)力對(duì)其動(dòng)力學(xué)的影響。 仿真初始條件為： 無(wú)人機(jī)前飛速度8 m/s， 迎角為2°， 無(wú)人機(jī)完成垂起改平飛過(guò)程后進(jìn)行盤旋， 四個(gè)垂起旋翼轉(zhuǎn)速為零， 兩個(gè)前拉螺旋槳轉(zhuǎn)速分別為1 000 r/min和2 000 r/min， 通過(guò)兩個(gè)前拉螺旋槳的拉力差形成力矩進(jìn)而發(fā)生盤旋運(yùn)動(dòng)， 其余狀態(tài)量初始時(shí)刻均為零。

7 結(jié)? 論

通過(guò)對(duì)垂起固定翼無(wú)人機(jī)動(dòng)力學(xué)建模和非定常氣動(dòng)力的研究， 可以得到以下結(jié)論：

（1） 采用多體動(dòng)力學(xué)基于凱恩方程能夠系統(tǒng)和方便地建立垂起固定翼無(wú)人機(jī)的動(dòng)力學(xué)模型， 通過(guò)遲滯環(huán)和ONERA方程能夠建立垂起固定翼無(wú)人機(jī)的非定常氣動(dòng)力模型。

（2） 通過(guò)仿真曲線圖4～6可以發(fā)現(xiàn)， 引入遲滯環(huán)將導(dǎo)致無(wú)人機(jī)俯仰振蕩速率減小， 其主要原因在于： 由于無(wú)人機(jī)機(jī)翼的氣動(dòng)中心位于其重心之后， 當(dāng)無(wú)人機(jī)做俯仰運(yùn)動(dòng)， 俯仰角速度為正值時(shí)， ?機(jī)翼氣動(dòng)力起阻尼作用， ?俯仰角速度為負(fù)值

時(shí)， 機(jī)翼氣動(dòng)力起促進(jìn)作用， 其產(chǎn)生的阻尼作用比促進(jìn)作用大， 因此， 整體呈現(xiàn)阻尼效果， 對(duì)于該類無(wú)人機(jī)進(jìn)行俯仰運(yùn)動(dòng)時(shí)引入遲滯環(huán)進(jìn)行修正是合理的。 引入遲滯環(huán)的本質(zhì)其實(shí)是引入氣流在流動(dòng)過(guò)

程中所產(chǎn)生的非定常效應(yīng)， 這種非定常效應(yīng)不

僅在無(wú)人機(jī)做俯仰等縱向運(yùn)動(dòng)時(shí)產(chǎn)生阻尼效果， 而且當(dāng)無(wú)人機(jī)進(jìn)行盤旋等橫航向與縱向耦合的運(yùn)動(dòng)時(shí)同樣也將產(chǎn)生阻尼效果， 其主要原因在于引入遲滯環(huán)將導(dǎo)致無(wú)人機(jī)的力矩系數(shù)隨運(yùn)動(dòng)而發(fā)生變化， 進(jìn)而使得無(wú)人機(jī)在運(yùn)動(dòng)過(guò)程中產(chǎn)生較大的阻尼力矩。

（3） 通過(guò)仿真曲線圖7～10可以發(fā)現(xiàn)， 兩者仿真結(jié)果整體變化趨勢(shì)一致， 但是變化過(guò)程有差異。 從圖9可以明顯發(fā)現(xiàn)， 無(wú)人機(jī)垂起速度引入非定常氣動(dòng)力與未引入有8%的差異， 引入非定常氣動(dòng)力模型的垂起速度更小， 在圖9的動(dòng)態(tài)過(guò)渡過(guò)程中， 即40～60 s， 兩者也具有較大差異， 引入非定常氣動(dòng)力模型與未引入非定常氣動(dòng)力模型兩者的最終上升速度存在40%的差異， 而上升速度是垂起改平飛過(guò)程中需要重點(diǎn)關(guān)注的物理量， 因此在控制

和設(shè)計(jì)過(guò)渡階段時(shí)， 應(yīng)以引入非定常氣動(dòng)力的動(dòng)力學(xué)模型為基準(zhǔn)， 這樣才能保證無(wú)人機(jī)的安全性。 圖10展示了引入非定常氣動(dòng)力和未引入非定常氣動(dòng)力模型的總能量曲線對(duì)比， 可以發(fā)現(xiàn)引入非定 常氣動(dòng)力模型的總能量更小， ?其主要原因在于無(wú)

式中各變量的定義和具體的非定常動(dòng)態(tài)流場(chǎng)數(shù)值方法可參考文獻(xiàn)[11-12]。

人機(jī)在過(guò)渡階段需要克服非定常氣動(dòng)力做功， 進(jìn)而導(dǎo)致其總能量減小。 垂起改平飛過(guò)渡過(guò)程的總能量曲線整體變化光滑， 該無(wú)人機(jī)能夠順利完成垂起改平飛過(guò)程。

（4） 通過(guò)仿真曲線圖11～12的對(duì)比可以得到兩個(gè)結(jié)論： a. 垂起固定翼無(wú)人機(jī)在盤旋運(yùn)動(dòng)中引入非定常氣動(dòng)力與未引入非定常氣動(dòng)力仿真結(jié)果最終趨于一致， 證明了本文采用ONERA方程的準(zhǔn)確性， 主要在于ONERA方程在小迎角情況下為線性函數(shù)， 其與傳統(tǒng)氣動(dòng)力模型相差不大， 兩者仿真結(jié)果在小迎角情況下是一致的; b. 通過(guò)圖12可以發(fā)現(xiàn)， 垂起固定翼無(wú)人機(jī)在大側(cè)滑角的情況下引入非定常氣動(dòng)力后其波動(dòng)的振幅更小， 收斂速度更快。 盤旋運(yùn)動(dòng)是一個(gè)耦合運(yùn)動(dòng)， 在大側(cè)滑角情

況下采用ONERA非線性方程建立了升力系數(shù)與阻力系數(shù)模型， ?由于這種情況下的升力系數(shù)與阻力系數(shù)與未引入非定常氣動(dòng)力的情況相比， 其值更小， 因此其波動(dòng)的振幅會(huì)更小。 通過(guò)圖11與圖12的對(duì)比可以發(fā)現(xiàn)， 盤旋運(yùn)動(dòng)過(guò)程中大側(cè)滑角情況其實(shí)就對(duì)應(yīng)著大迎角情況。

參考文獻(xiàn)：

[1] ?Kim S K，? Tilbury D M. Mathematical Modeling and Experimental Identification of a Model Helicopter[C]∥AIAA Modeling and Simulation Technologies Conference ?and Enhibit，? AIAA-98-4357，? Boston， 1998，? 203-213.

[2] ?DeTore J，? Conway S. Technology Needs for High Speed Rotorcraft（3）， NASA-CR-177592[R].?? 1991.

[3] ?李家樂(lè)，? 王正平. 基于Lagrange方法的單旋翼飛行器動(dòng)力學(xué)建模[J]. 飛行力學(xué)，? 2016，? 34（4）： 15-18.

Li Jiale，? Wang Zhengping. Dynamics Modeling for Monowing Rotorcraft Using Lagrange Method[J]. Flight Dynamics，? 2016，? 34（4）： 15-18.（in Chinese）

[4] ?李家樂(lè). 弧翼飛行器飛行控制與仿真[D]. 西安： 西北工業(yè)大學(xué)，? 2016： 32-43.

Li Jiale. Research on Control and Simulation for a Cambered Wing Vehicle[D]. Xian： Northwestern Polytechnical University，? 2016： 32-43.（in Chinese）

[5] ?Hogan F R，? Forbes J R. Modeling of Spherical Robots Rolling on Generic Surfaces[J]. Multibody System Dynamics，? 2015，? 35（1）： 91-109.

[6] ?曹蕓蕓. 傾轉(zhuǎn)旋翼飛行器飛行動(dòng)力學(xué)數(shù)學(xué)建模方法研究[D]. 南京： 南京航空航天大學(xué)，? 2012： 29-63.

Cao Yunyun. Research on Mathematical Modeling Method for Tilt Rotor Aircraft Flight Dynamics[D]. Nanjing： Nanjing University of Aeronautics and Astronautics，? 2012： 29-63.（in Chinese）

[7] ?Kane T R，? Likins P W，? Levinson D A. Spacecraft Dynamics[M]. New York： McGrawHill Book Company，? 1983.

[8] ?Zhao Zhenjun，? Ren Gexue. Multibody Dynamic Approach of Flight Dynamics and Nonlinear Aeroelasticity of Flexible Aircraft [J]. AIAA Journal，? 2011，? 49 （1）： 41-54.

[9] ?喬宇航，? 馬東立，? 李陟. 螺旋槳/機(jī)翼相互干擾的非定常數(shù)值模擬[J]. 航空動(dòng)力學(xué)報(bào)，? 2015， 30（6）： 1366-1373.

Qiao Yuhang，? Ma Dongli，? Li Zhi. Unsteady Numerical Simulation of Propeller/Wing Interaction[J]. Journal of Aerospace Power，? 2015， 30（6）： 1366-1373.（in Chinese）

[10] ?袁先旭，? 張涵信，? 謝昱飛. 基于CFD方法的俯仰靜、 動(dòng)導(dǎo)數(shù)數(shù)值計(jì)算[J]. 空氣動(dòng)力學(xué)學(xué)報(bào)，? 2005，? 23（4）： 458-463.

Yuan Xianxu，? Zhang Hanxin，? Xie Yufei. The Pitching Static/Dynamic Derivatives Computation Based on CFD Methods[J]. Acta Aerodynamica Sinica，? 2005，? 23（4）： 458-463. （in Chinese）

[11] ?East R A，? Hutt G R. Comparison of Predictions and Experimental Data for Hypersonic Pitching Motion Stability [J]. Journal of Spacecraft and Rockets，? 1998，? 25（3）：225-233.

[12] ?劉羅成，? 徐敏， 姚偉剛， 等. 翼型縱向俯仰非定常氣動(dòng)力研究[J]. 飛行力學(xué)，? 2011，? 29（4）： 28-31.

Liu Luocheng，? Xu Min，? Yao Weigang，? et al. Study of Hypersonic Airfoil Pitch Oscillation Unsteady Aerodynamic Force[J]. Flight Dynamics，? 2011，? 29（4）： 28-31.（in Chinese）

[13] ?Kim T， Dugundji J. Nonlinear Large Amplitude Aeroelastic Behavior of Composite Rotor Blades[J]. AIAA Journal， ?1993，? 31（8）： 1489-1497.

[14] ?林學(xué)海，? 任勇生. 基于ONERA氣動(dòng)力模型的風(fēng)力機(jī)葉片顫振時(shí)域分析[J]. 山東科技大學(xué)學(xué)報(bào)： 自然科學(xué)版，? 2009，? 28（3）： 56-60.

Lin Xuehai，? Ren Yongsheng. Analysis of Flutter TimeDomain for Blades of Wind Turbine Based on the ONERA Aerodynamic Model[J]. Journal of Shandong University of Science and Technology： Natural Science，? 2009，? 28（3）： 56-60.（in Chinese）

[15] ?Forbes J R，? Barfoot T D，? Damaren C J. Dynamic Modeling and Stability Analysis of a PowerGenerating Tumbleweed Rover[J]. Multibody System Dynamics，? 2010，? 24（4）： 413-439.

[16] ?蔣達(dá). 新型無(wú)刷直流電動(dòng)機(jī)建模與仿真研究[D]. 南寧： 廣西大學(xué)，? 2012：? 17-19.

Jiang Da. Research on Modeling and Simulation of New Type Brushless DC Motor[D]. Nanning：? Guangxi University，? 2012：? 17-19.（in Chinese）

[17] ?夏丹， 陳維山， 劉軍考， 等. 基于Kane方法的仿魚機(jī)器人波狀游動(dòng)的動(dòng)力學(xué)建模[J]. 機(jī)械工程學(xué)報(bào)，? 2009， 45（6）： 41-49.

Xia Dan，? Chen Weishan，? Liu Junkao，? et al. Dynamic Modeling of a Fishlike Robot with Undulatory Motion Based on Kanes Method[J]. Journal of Mechanical Engineering，? 2009，? 45（6）： 41-49.（in Chinese）

[18] ?Tarn T J，? Shoults G A，? Yang S P. A Dynamic Model of Underwater Vehicle with a Robotics Manipulator Using Kanes Method[J]. Autonomous Robots，? 1996，? 3（2-3）： 269-283.