徐孝武, 張煒, 詹浩

(1.西北工業(yè)大學(xué) 航空學(xué)院, 陜西 西安 710072; 2.陜西省試驗飛機設(shè)計與試驗技術(shù)工程實驗室, 陜西 西安 710072)

變體飛機能夠主動改變氣動外形以適應(yīng)不同的飛行環(huán)境或飛行任務(wù)。經(jīng)過研究，變后掠飛機的非對稱變形有利于提高飛機的抗擾動能力(比如抗側(cè)風(fēng)能力)，能夠更有效地完成任務(wù)規(guī)劃。相比之下，折疊機翼變體飛機的變形量更大，更加有利于提高飛機的機動性和敏捷性，提高飛機的綜合性能。因此有必要專門針對非對稱變形引起的氣動力和動力學(xué)響應(yīng)特性進行研究。詳細的氣動模型在飛行仿真和控制系統(tǒng)設(shè)計中占據(jù)著非常重要的作用，變體飛機變形過程中氣動特性更加復(fù)雜，因此，如何建立精準(zhǔn)的氣動力數(shù)學(xué)模型至關(guān)重要。變形過程中氣動參數(shù)呈現(xiàn)強烈的非線性[1-5]，尤其是在非對稱變形時，又會至少增加一個變形參數(shù)，氣動參數(shù)的變化規(guī)律顯得更為復(fù)雜和難以描述，氣動數(shù)據(jù)模型復(fù)雜、高階且含有強非線性。目前針對變體飛機的研究文獻中在進行氣動力建模時多以數(shù)據(jù)插值的方法求解連續(xù)變形狀態(tài)的氣動力[6-8]。該方法不能得知變形參數(shù)與氣動系數(shù)之間的具體函數(shù)關(guān)系，計算精度也難以保證，無法得到變形過程的高精度氣動模型。

常用的氣動力數(shù)學(xué)模型有線性氣動模型和多項式模型。其中多項式模型是線性模型的直接推廣，在一定范圍內(nèi)能夠較好地描述非線性效應(yīng)，但其近似能力很大程度受限于階次選取。于是學(xué)者們又發(fā)明了多項式樣條函數(shù)，可以用低階項很好地逼近各種非線性，通過增加節(jié)點數(shù)增加多項式段數(shù)，從而提高近似能力[9]。文獻[10]采用多項式樣條函數(shù)研究輕型飛機在失速、過失速區(qū)的氣動力數(shù)學(xué)模型。多項式樣條函數(shù)模型與多項式模型一樣需要預(yù)設(shè)模型階次，選擇不合適仍然會有較大誤差，且無法擬合散亂數(shù)據(jù)[9]。近年來，學(xué)者們研究了一種由伯恩斯坦多項式組成的多元樣條函數(shù)[11-13]，并將該方法用于參數(shù)辨識[14-15]。研究表明，該方法用于飛行器氣動辨識能夠得到高精度的氣動力模型[16-18]。

本文提出了一種采用多元樣條函數(shù)理論進行變體飛機氣動力建模的方法。該方法無需預(yù)判具體模型結(jié)構(gòu)，通過數(shù)據(jù)擬合直接得到了變形參數(shù)和氣動系數(shù)之間的具體函數(shù)關(guān)系，能夠得到詳細描述變體飛機任意對稱和非對稱變形狀態(tài)下的氣動力模型。

1 多元樣條函數(shù)基礎(chǔ)理論

本節(jié)簡要介紹多元樣條函數(shù)基礎(chǔ)理論[13-14,19]。

1.1 質(zhì)心坐標(biāo)系

單形是組成多元樣條函數(shù)的基本單元。例如,二元單形是一個三角形，三元單形一個是四面體。n元單形t由(n+1)個頂點連接組成,描述如下:

t=〈v0,v1,…,vn〉

(1)

與單形固連的局部坐標(biāo)系稱為質(zhì)心坐標(biāo)系。質(zhì)心坐標(biāo)系上的每一個點x=(x1,x2,…,xn)都可以由單形t的頂點按權(quán)重線性相加得到,描述如下:

x=∑ni=0bivi, ∑ni=0bi=1

(2)

n元單形t中任何一點的質(zhì)心坐標(biāo)b(x)=(b0,b1,…,bn)可由下式計算得到:

b1

b2

?

bn=[(v1-v0) (v2-v0) …(vn-v0)]-1·

(x-v0)=Λ(x-v0)

(3)

b0=1-∑ni=1bi

(4)

1.2 樣條函數(shù)和樣條空間

三角測量Γ是將一個區(qū)域劃分為J個非交叉相連的單形,描述如下:

Γ:=∪Ji=1ti,ti∩tj∈{?,}, ?ti,tj∈Γ

(5)

給定一個區(qū)域通過三角測量得到J個單形,每一個單形tj均可由維度為d的多項式ptj(x)描述,按照預(yù)設(shè)連續(xù)條件組合成多元樣條函數(shù)s(x):

s(x)=∑Jj=1δj(x)ptj(x),δj(x)=1,x∈tj0,x?tj

(6)

樣條空間是指在三角測量Γ區(qū)域內(nèi),給定維度d和連續(xù)階數(shù)Cr的所有樣條函數(shù)s的集合。

(Γ):={s∈Cr(Γ):st∈Ρd,?t∈Γ}

(7)

式中,Ρd指所有維度為d的多項式集合。

1.3 B-form

質(zhì)心坐標(biāo)b(x)=(b0,b1,…,bn)寫成維度為d的伯恩斯坦多項式形式如下:

多元系數(shù)κ定義如下:

κ=(κ0,κ1,…,κn)∈Nn+1

(9)

κ的1-范數(shù)為:

|κ|=κ0+κ1+…+κn=d,d≥0

(10)

κ的階乘定義如下:

κ!=κ0!κ1!…κn!

(11)

κ的所有組合按照詞典排序法排序,排列順序如下:

κ∈{(d,0,0,…,0),(d-1,1,0,…,0),

(d-1,0,1,…,0),…,(0,…,0,1,d-1),

(0,…,0,0,d)}

(12)

κ的所有組合數(shù)量為:

=(d+n)!n!d!

(13)

給定三角測量Γ,任何多項式均可寫成如下形式:

(14)

B-系數(shù)可以描述為以下矢量形式:

∈R

(15)

基本多項式按照κ的組合排列可描述為矢量形式:

∈R

(16)

根據(jù)(15)和(16)式,(14)式的矢量形式可以描述如下:

B-系數(shù)的全局矢量c描述如下:

c=[ct1Tct2T…ctJT]T∈RJ·×1(18)

基本多項式的全局矢量描述如下：

…

根據(jù)(6)和(17)式，1組觀測值(x(i),y(i))可以用B-form多項式來描述如下：

根據(jù)(18)和(19)式，(20)式的矢量形式可以描述如下：

y(i)=Bd(b(x(i)))c+ε(i)=X(i)c+ε(i)(21)

綜上,對于所有觀測值,樣條函數(shù)的線性回歸模型描述如下:

Y=Xc+ε∈RN×1

(22)

1.4 連續(xù)條件

通過三角測量在給定區(qū)域得到J個單形,各單形之間相鄰邊緣之間的連續(xù)性通過連續(xù)條件來約束。給定連續(xù)階數(shù)為Cr時,每個相鄰邊緣的連續(xù)條件個數(shù)R按下式計算得到:

R=∑rm=0(d-m-n-1)!(n-1)!(d-m)!

(23)

定義H為平滑矩陣,可以得到全局區(qū)域內(nèi)的連續(xù)條件方程:

Hc=0,H∈RR·E×J·

(24)

式中,J為單形個數(shù),E為單形的相鄰邊緣個數(shù),和R分別由(13)和(23)式得到。

2 基于多元樣條函數(shù)的氣動力模型

2.1 幾何模型

以某型折疊機翼變體飛機為例,折疊段可以單獨向上折疊120°,同時外段機翼始終保持水平,外形如圖1所示,各翼段的主要幾何參數(shù)如表1所示。

圖1 折疊機翼變體飛機模型平面圖

幾何參數(shù)機身段折疊段外段翼段展長/m0．300．300．55參考面積/m20．2180．1310．165根弦長/m0．9000．7250．468梢弦長/m0．7250．4680．217前緣后掠角/(°)353535

2.2 氣動力模型

研究結(jié)果表明,在變形速率不大時,變體飛機的氣動力可忽略非定常效應(yīng),按準(zhǔn)定常計算[8,20]。變體飛機的氣動力僅受飛機當(dāng)前的氣動外形和飛行狀態(tài)影響。穩(wěn)定軸系繞自身橫軸轉(zhuǎn)動α(迎角)角度與體軸系重合,穩(wěn)定軸系上氣動力和力矩系數(shù)分別是:升力系數(shù)CL,阻力系數(shù)CD,側(cè)力系數(shù)CY,滾轉(zhuǎn)力矩系數(shù)Cl,俯仰力矩系數(shù)Cm,偏航力矩系數(shù)Cn。體軸系上的氣動力和力矩計算公式如下:

FA,x=QSw(CLsinα-CDcosα)

FA,y=QSwCY

FA,z=QSw(-CDsinα-CLcosα)

LA=QSwbw(Clcosα-Cnsinα)

MA=QSwcACm

NA=QSwbw(Cncosα+Clsinα)

(25)

式中,Q=12ρV2為動壓;ρ為空氣密度;V為空速;Sw為全機機翼參考面積;cA為全機機翼參考平均氣動弦長;bw為全機機翼參考翼展。

在非對稱變形過程中或變形完成后非對稱飛行時,還會產(chǎn)生附加側(cè)力Fu,y、滾轉(zhuǎn)力矩Lu和偏航力矩Nu,它們是變形參數(shù)μ=[μ1μ2]T的函數(shù)[21]。本文引入3個非對稱變形引起的氣動參數(shù):附加側(cè)力系數(shù)CYur、滾轉(zhuǎn)力矩系數(shù)Clur和偏航力矩系數(shù)Cnur,定義如下:

CYur=Fu,y(QSw)

Clur=Lu(QSwbw)

Cnur=Nu(QSwbw)

(26)

變體飛機變形過程的氣動系數(shù)線性化模型如下:

CL=CL0+CLα·α+CLδe·δe

CY=CYβ·β+CYδr·δr+CYur·CL

Cl=Clβ·β+Clδa·δa+Clδr·δr+

Cl·+Cl·+Clur·CL

Cm=Cm0+Cmα·α+Cmδe·δe+

Cm·+Cm·

Cn=Cnβ·β+Cnδa·δa+Cnδr·δr+

Cn·+Cn·+Cnur·CL

(27)

式中,β為側(cè)滑角;為無量綱迎角變化率;,和分別為無量綱滾轉(zhuǎn)、俯仰和偏航角速度;δa,δe和δr分別為副翼、升降舵和方向舵偏角;CL0為基本升力系數(shù);CLα和CLδe分別為升力系數(shù)對迎角和升降舵偏角的導(dǎo)數(shù);CD0和Cdi分別為零升阻力系數(shù)和升致阻力系數(shù);CYβ和CYδr分別為側(cè)力系數(shù)對側(cè)滑角和方向舵偏角的導(dǎo)數(shù);Clβ,Clδa,Clδr,Cl和Cl分別為滾轉(zhuǎn)力矩系數(shù)對側(cè)滑角、副翼偏角、方向舵偏角、無量綱滾轉(zhuǎn)角速度和無量綱偏航角速度的導(dǎo)數(shù);Cm0為零升俯仰力矩系數(shù);Cmα,Cmδe,Cm和Cm分別為俯仰力矩系數(shù)對迎角、升降舵偏角、無量綱俯仰角速度和無量綱迎角變化率的導(dǎo)數(shù);Cnβ,Cnδa,Cnδr,Cn和Cn分別為偏航力矩系數(shù)對側(cè)滑角、副翼偏角、方向舵偏角、無量綱滾轉(zhuǎn)角速度和無量綱偏航角速度的導(dǎo)數(shù)。

2.3 多元樣條函數(shù)表達形式

當(dāng)變體飛機在變形過程中保持小迎角和小側(cè)滑角狀態(tài)時,仍可用形如(27)式的線性氣動模型表達當(dāng)前狀態(tài)的氣動模型。其中的每一個氣動參數(shù)都是變形參數(shù)的非線性函數(shù),參數(shù)模型結(jié)構(gòu)未知。本文所示變體飛機的變形參數(shù)μ=[μ1μ2]T有2個變量,因此為二元樣條函數(shù)模型。在建模過程中,將每一個氣動參數(shù)都視為一個樣條函數(shù)模型進行建模,于是可得變體飛機氣動系數(shù)樣條函數(shù)模型最終表達形式如下:

s1=s11(μ1,μ2)+s12(μ1,μ2)·α+

s13(μ1,μ2)·δe

s3=s31(μ1,μ2)·β+s32(μ1,μ2)·δr+

s33(μ1,μ2)·CL

s4=s41(μ1,μ2)·β+s42(μ1,μ2)·δa+

s43(μ1,μ2)·δr+s44(μ1,μ2)·+

s45(μ1,μ2)·+s46(μ1,μ2)·CL

s5=s51(μ1,μ2)+s52(μ1,μ2)·α+s53(μ1,μ2)·δe+

s54(μ1,μ2)·+s55(μ1,μ2)·

s6=s61(μ1,μ2)·β+s62(μ1,μ2)·δa+

s63(μ1,μ2)·δr+s64(μ1,μ2)·+

s65(μ1,μ2)·+s66(μ1,μ2)·CL

(28)

根據(jù)(21)式,變體飛機氣動系數(shù)模型可以描述為線性回歸的形式。例如,CL描述如下:

CL(i)=

c12

c13

(29)

式中,cij為樣條函數(shù)sij的B-系數(shù)。由此,(28)式能夠表達出任意變形狀態(tài)的氣動力系數(shù)模型。

3 樣條函數(shù)模型估計

3.1 樣條函數(shù)模型估計流程

樣條函數(shù)模型估計流程如圖2所示,主要分為4大步驟:模型結(jié)構(gòu)選取、模型估計、模型驗證和模型確定。其中模型結(jié)構(gòu)選取直接影響模型的性能,最為復(fù)雜和重要,一般需要多次選取不同的維度對結(jié)果進行比較。三角測量的構(gòu)造受限于數(shù)據(jù)分布情況、數(shù)據(jù)量和計算能力,需要權(quán)衡選取。

圖2 樣條函數(shù)模型估計流程圖

3.2 模型結(jié)構(gòu)選取

1) 維度選擇

一般來講,提高維度會增加模型的精度,但同時也會使模型變得更加復(fù)雜,計算更加耗時,因此需要選擇合適的維度。文獻[11]研究表明,d≥3r+2時,樣條函數(shù)總是能夠得到較好的精度。本文考慮數(shù)據(jù)0階連續(xù),即r=0,因此初次選取維度為d=2。

2) 構(gòu)造三角測量;

圖3 三角測量(Γ4)

圖4 B-form排列(s∈Γ4))

3) 確定模型空間結(jié)構(gòu)

由(6)式和(17)式可得圖4所示B-form排列的樣條函數(shù)模型結(jié)構(gòu)為:

s(x)=δ1(x)pt1(x)+δ2(x)pt2(x)+

δ3(x)pt3(x)+δ4(x)pt4(x)=

由(12)式可得κ的排列為:

κ∈{(2,0,0),(1,1,0),(1,0,1),(0,2,0),

(0,1,1),(0,0,2)}

(31)

由(15)式可得B-系數(shù)ctj為:

由(16)式可得每一個三角的基本多項式矢量描述如下:

由(18)式可得B-系數(shù)的全局矢量c為:

(34)

由(19)式可得基本多項式的全局矢量為:

B2(b(x))=

3.3 模型估計

1) 計算B-form矩陣

由圖3可知各頂點坐標(biāo):v0=(0,0),v1=2π3,0,v2=2π3,2π3,v3=0,2π3,v4=π3,π3。三角形區(qū)域分別為:t1=〈v0v1v4〉,t2=〈v1v2v4〉,t3=〈v0v3v4〉,t4=〈v2v3v4〉。由(3)和(4)式容易得到每個三角形區(qū)域內(nèi)任一點x=(μ1,μ2)的質(zhì)心坐標(biāo)b(x)=(b0,b1,b2)。

由(21)式可得每個單形的線性回歸模型為:

ytj(1)

ytj(2)

?

?

(36)

式中,Nj為每個單形的數(shù)據(jù)組數(shù)。上式寫成標(biāo)準(zhǔn)形式如下:

Ytj=Xtjctj+εNj×1∈RNj×1

(37)

于是得到樣條函數(shù)空間內(nèi)所有數(shù)據(jù)的線性回歸模型如下:

Yt1

Yt2

Yt3

Yt4=Xt1000

0Xt200

00Xt30

000Xt4ct1

ct2

ct3

ct4+εN×1

(38)

至此,得到了形如(22)式的線性回歸模型,于是得到B-form矩陣X。

2) 計算平滑矩陣

于是可得平滑矩陣:

H=000-100100000000000000000

0000-10001000000000000000

00000-1000001000000000000

-100000000000100000000000

00-1000000000001000000000

00000-1000000000001000000

000000000-100000000100000

0000000000-10000000001000

00000000000-1000000000001

000000000000000-100000100

0000000000000000-10000010

3) B-系數(shù)估計

根據(jù)(22)式的線性回歸模型,可以用最小二乘估計方法估計B-系數(shù):

minc12(Y-Xc)T(Y-Xc)

(40)

給定約束條件Hc=0,可以得到最小二乘估計表達式:

H0+XTY

0=C1C2

C3C4XTY

0

(41)

可由下式計算出B-系數(shù):

=C1XTY

(42)

=-

0.7880.8000.7990.6040.7050.708 …

0.6040.6180.7050.4200.6110.708 …

0.7880.8000.7990.6040.7050.708 …

0.4200.6180.6110.6040.7050.708T

4 最終模型結(jié)構(gòu)確定

4.1 模型驗證

1) 誤差分析理論

χval為需要驗證的數(shù)據(jù)集合,表示為

χval=∪Ni=1xi

(43)

實際輸出與估計值之差,即殘差ε(χval)為

(44)

殘差均方根RMSε為

RMSε=1N∑Ni=1(ε(xi))2

(45)

相對殘差均方根RMSrelε為

RMSrelε=RMSε|maxC(χval)-minC(χval)|

(46)

最大相對殘差εrelmax為

εrelmax=max|ε(χval)||maxC(χval)-minC(χval)|

(47)

2) 計算結(jié)果分析

以Cmδe(μ1,μ2)為例,計算得出殘差各相關(guān)值見表2。

表2 Cmδe(μ1,μ2)對應(yīng)的殘差值

4.2 確定最終模型

模型最終確定依據(jù)以下3條基本準(zhǔn)則:

1) 0.01

2) RMSrelε<0.01時,模型品質(zhì)視為優(yōu)異,不再提高維度;

3) 僅當(dāng)提高維度有顯著的品質(zhì)提升時才選擇更高的維度。

所有氣動參數(shù)的最終估計結(jié)果如下:

Cl(μ1,μ2)∈Γ4);Clur(μ1,μ2)∈Γ4);

Cm(μ1,μ2)∈Γ4);Cnβ(μ1,μ2)∈Γ4);

Cn(μ1,μ2)∈Γ4);Cnur(μ1,μ2)∈Γ4)

確定的模型誤差分析結(jié)果如表3所示。

表3 模型誤差分析結(jié)果

Cmδe擬合數(shù)據(jù)圖示如下:

圖5 擬合數(shù)據(jù)(Cmδe(μ1,μ2)∈Γ4))與原始數(shù)據(jù)對比

4.3 總體坐標(biāo)系中的表達式

前面得到的模型是在質(zhì)心坐標(biāo)系中描述,質(zhì)心坐標(biāo)系為局部坐標(biāo)系,因此需要轉(zhuǎn)換到描述飛機受力情況和運動特性的總體坐標(biāo)系中。將質(zhì)心坐標(biāo)b(x)=(b0,b1,b2)和B-系數(shù)的估計值代入(17)式可得總體坐標(biāo)系中的多項式描述形式:

ptj(μ1,μ2)=aj0+aj1μ1+aj2μ2+aj3μ1μ2+

(48)

對應(yīng)各參數(shù)值見表4。

表4 Cmδe(μ1,μ2)∈Γ4)對應(yīng)的多項式系數(shù)

5 結(jié) 論

本文采用多元樣條函數(shù)理論對變體飛機氣動力進行建模，以某型折疊機翼變體飛機為例，建立了基于多元樣條函數(shù)模型結(jié)構(gòu)的氣動力方程，給出了詳細的樣條函數(shù)模型結(jié)構(gòu)和估計流程，通過計算得到了該變體飛機的氣動力數(shù)學(xué)模型。通過誤差分析發(fā)現(xiàn)，當(dāng)維度增加時，模型品質(zhì)未必一直增加，維度選取過高可能會造成擬合失真，因此不能單純的通過提高維度來提高模型品質(zhì)，當(dāng)提高維度得到的模型無法滿足要求時需要考慮重新構(gòu)造三角測量。結(jié)果表明，該方法無需預(yù)知變形過程的氣動變化規(guī)律和氣動模型具體結(jié)構(gòu)，通過誤差分析選擇合理維度，能夠得到詳細描述變體飛機任意對稱和非對稱變形狀態(tài)的氣動力模型，解決了變體飛機氣動模型結(jié)構(gòu)未知情況下的氣動力建模問題以及在非對稱狀態(tài)下的氣動力描述問題。

參考文獻：

[1] Abdulrahim M, Lind R. Control and Simulation of a Multi-Role Morphing Micro Air Vehicle[R]. AIAA-2005-6481

[2] Bourdin P, Gatto A, Friswell M I. The Application of Variable Cant Angle Winglets for Morphing Aircraft Control[R]. AIAA-2006-3660

[3] Grant D T, Abdulrahimy M, Lind R. Flight Dynamics of a Morphing Aircraft Utilizing Independent Multiple-Joint Wing Sweep[R]. AIAA-2006-6505

[4] Grant D T, Chakravarthy A, Lind R. Modal Interpretation of Time-Varying Eigenvectors of Morphing Aircraft[R]. AIAA-2009-5848

[5] Grant D T. Modeling and Dynamic Analysis of a Multi-Joint Morphing Aircraft[D]. Florida, University of Florida, 2009

[6] Seigler T M. Dynamics and Control of Morphing Aircraft[D]. Virginia, Virginia Polytechnic Institute and State University, 2005

[7] Seigler T M, Neal D A. Analysis of Transition Stability for Morphing Aircraft[J]. Journal of Guidance Control and Dynamics. 2009, 32(6): 1947-1953

[8] Yue T, Wang L, Ai J. Longitudinal Linear Parameter Varying Modeling and Simulation of Morphing Aircraft[J]. Journal of Aircraft, 2013, 50(6): 1673-1681

[9] De Visser C C, Van Kampen E, Chu Q P, et al. A New Framework for Aerodynamic Model Identification with Multivariate Splines[R]. AIAA-2013-4748

[10] Batterson J G. Analysis of Oscillatory Motion of a Light Airplane at High Value of Lift Coefficient[R]. NASA TM-84563, 1983

[11] Awanou G, Lai M, Wenston P. The Multivariate Spline Method for Scattered Data Fitting and Numerical Solutions of Partial Differential Equations[A]. Wavelets and Splines, Brentwoxl Nashboro, 2005: 24-75

[12] Lai M J. Multivariate Splines for Data Fitting and Approximation[C]∥12th Approximation Theory Conference, 2007

[13] Lai M J, Schumaker L L. Spline Function on Triangulations[M]. New York, Cambridge University Press, 2007: 18-23

[14] De Visser C C, Chu Q P, Mulder J A. A New Approach to Linear Regression with Multivariate Splines[J]. Automatica, 2009, 45(12): 2903-2909

[15] De Visser C C, Chu Q P, Mulder J A. Differential Constraints for Bounded Recursive Identification with Multivariate Splines[J]. Automatica, 2011, 47(9): 2059-2066

[16] De Visser C C, Mulder J A, Chu Q P. Validating the Multidimensional Spline Based Global Aerodynamic Model for the Cessna CitationⅡ[R]. AIAA-2011-6356

[17] Sun L G, De Visser C C, Chu Q P, et al. Online Aerodynamic Model Identification Using a Recursive Sequential Method for Multivariate Splines[J]. Journal of Guidance Control and Dynamics, 2013, 36(5): 1278-1288

[18] Tol H J, De Visser C C, Van Kampen E, et al. Nonlinear Multivariate Spline-Based Control Allocation for High-Performance Aircraft[J]. Journal of Guidance Control and Dynamics, 2014, 37(6): 1840-1862

[19] De Visser C C. Global Nonlinear Model Identification with Multivariate Splines[D]. Zuthpen, Delft University of Technology, 2011

[20] An J, Yan M, Zhou W, et al. Aircraft Dynamic Response to Variable Wing Sweep Geometry[J]. Journal of Aircraft, 1988, 25(3): 216-221

[21] Xu X W, Zhang W. Research on Methods of Dynamic Modeling and Simulation for Morphing Wing Aircraft[R]. ICAS 2012-698