吳翰 王正平 周 洲 王睿

摘 要：

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

關鍵詞：???? ??垂起固定翼無人機; 遲滯環(huán); ONERA方程; 多體動力學; 凱恩方程

中圖分類號：??? ??TJ011;? V212 ??文獻標識碼：??? ?A ??文章編號：??? ??1673-5048（2019）02-0021-08

0 引? 言

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

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

非線性氣動力方程建立無人機垂起改平飛過渡過程中的非定常氣動力模型， 完善垂起改平飛過渡階段的動力學模型。

1 垂起固定翼無人機多剛體劃分

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

2 機翼遲滯效應分析

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

當無人機在做俯仰運動時， 機翼所產(chǎn)生的氣動力為非定常氣動力， 其控制方程為描述氣體非定常運動的NS方程[11]：

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

2.2 無人機機翼動導數(shù)求解

以遲滯環(huán)為基準， 推導無人機機翼俯仰運動動導數(shù)， 動導數(shù)是氣動載荷對狀態(tài)變化率（速度大?。┑膶?shù)， 是描述無人機動態(tài)變化過程的重要參數(shù)， 應用積分法求取無人機動導數(shù)。 假設機翼在俯仰自由度下進行小振幅強迫振蕩運動， 其迎角為

式中： k=（ωb/v）， 為振蕩運動的減縮頻率;? t為無量綱時間; θ=Δα=α-α0為機翼俯仰角。

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

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

3 ONERA非線性氣動力方程分析

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

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

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

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

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

通過上式可以得到 ΔCL α ， 將等式（12）和 ΔCL α 帶入式（10）中便可以得到相應情況下的升力系數(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ù)確定。 由于無人機在過渡階段對其側向力系數(shù)Cyβ影響較小， 因此側向力系數(shù)采用準定常假設進行求解即可。

4 機翼非定常動力學模型建立

4.1 機翼非定常氣動力模型

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

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

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

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

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

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

4.2 機翼非定常力矩模型

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

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為橫航向動導數(shù)， 表示該垂起固定翼無人機滾轉和偏航時左機翼所產(chǎn)生的非定常氣動力矩;? Cmθ， Cmθ? · 為無人機做俯仰運動時所產(chǎn)生非定常氣動力矩系數(shù);? S為機翼面積;? l， c分別為左機翼展長和弦長;? p， q， r為機翼的角速度分量。 需要特別注意的是主要通過無人機縱向動導數(shù)Cmθ， Cmθ? · 引入無人機俯仰運動中機翼所產(chǎn)生的遲滯效應， 而且遲滯效應僅在無人機迎角小于15°時才具有。

5 垂起固定翼無人機動力學模型建立

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

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

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

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

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

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

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

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

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

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

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

將式（18）和式（19）疊加便可得到該無人機的動力學模型[18]為

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

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

6 數(shù)值仿真

6.1 俯仰運動仿真

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

6.2 垂起改平飛過渡階段仿真

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

6.3 盤旋運動仿真

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

7 結? 論

通過對垂起固定翼無人機動力學建模和非定常氣動力的研究， 可以得到以下結論：

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

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

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

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

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

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

和設計過渡階段時， 應以引入非定常氣動力的動力學模型為基準， 這樣才能保證無人機的安全性。 圖10展示了引入非定常氣動力和未引入非定常氣動力模型的總能量曲線對比， 可以發(fā)現(xiàn)引入非定 常氣動力模型的總能量更小， ?其主要原因在于無

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

人機在過渡階段需要克服非定常氣動力做功， 進而導致其總能量減小。 垂起改平飛過渡過程的總能量曲線整體變化光滑， 該無人機能夠順利完成垂起改平飛過程。

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

況下采用ONERA非線性方程建立了升力系數(shù)與阻力系數(shù)模型， ?由于這種情況下的升力系數(shù)與阻力系數(shù)與未引入非定常氣動力的情況相比， 其值更小， 因此其波動的振幅會更小。 通過圖11與圖12的對比可以發(fā)現(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] ?李家樂，? 王正平. 基于Lagrange方法的單旋翼飛行器動力學建模[J]. 飛行力學，? 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] ?李家樂. 弧翼飛行器飛行控制與仿真[D]. 西安： 西北工業(yè)大學，? 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] ?曹蕓蕓. 傾轉旋翼飛行器飛行動力學數(shù)學建模方法研究[D]. 南京： 南京航空航天大學，? 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] ?喬宇航，? 馬東立，? 李陟. 螺旋槳/機翼相互干擾的非定常數(shù)值模擬[J]. 航空動力學報，? 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方法的俯仰靜、 動導數(shù)數(shù)值計算[J]. 空氣動力學學報，? 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] ?劉羅成，? 徐敏， 姚偉剛， 等. 翼型縱向俯仰非定常氣動力研究[J]. 飛行力學，? 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] ?林學海，? 任勇生. 基于ONERA氣動力模型的風力機葉片顫振時域分析[J]. 山東科技大學學報： 自然科學版，? 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]. 南寧： 廣西大學，? 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]. 機械工程學報，? 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.