基于多元樣條函數(shù)的變體飛機非對稱氣動力建模

徐孝武, 張煒, 詹浩

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

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

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

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

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

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

1.1 質(zhì)心坐標系

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

t=〈v0,v1,…,vn〉

(1)

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

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

(2)

n元單形t中任何一點的質(zhì)心坐標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)描述,按照預設(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ì)心坐標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é)果表明,在變形速率不大時,變體飛機的氣動力可忽略非定常效應,按準定常計算[8,20]。變體飛機的氣動力僅受飛機當前的氣動外形和飛行狀態(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ù)對迎角和升降舵偏角的導數(shù);CD0和Cdi分別為零升阻力系數(shù)和升致阻力系數(shù);CYβ和CYδr分別為側(cè)力系數(shù)對側(cè)滑角和方向舵偏角的導數(shù);Clβ,Clδa,Clδr,Cl和Cl分別為滾轉(zhuǎn)力矩系數(shù)對側(cè)滑角、副翼偏角、方向舵偏角、無量綱滾轉(zhuǎn)角速度和無量綱偏航角速度的導數(shù);Cm0為零升俯仰力矩系數(shù);Cmα,Cmδe,Cm和Cm分別為俯仰力矩系數(shù)對迎角、升降舵偏角、無量綱俯仰角速度和無量綱迎角變化率的導數(shù);Cnβ,Cnδa,Cnδr,Cn和Cn分別為偏航力矩系數(shù)對側(cè)滑角、副翼偏角、方向舵偏角、無量綱滾轉(zhuǎn)角速度和無量綱偏航角速度的導數(shù)。

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

當變體飛機在變形過程中保持小迎角和小側(cè)滑角狀態(tài)時,仍可用形如(27)式的線性氣動模型表達當前狀態(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)選取直接影響模型的性能,最為復雜和重要,一般需要多次選取不同的維度對結(jié)果進行比較。三角測量的構(gòu)造受限于數(shù)據(jù)分布情況、數(shù)據(jù)量和計算能力,需要權(quán)衡選取。

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

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

1) 維度選擇

一般來講,提高維度會增加模型的精度,但同時也會使模型變得更加復雜,計算更加耗時,因此需要選擇合適的維度。文獻[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可知各頂點坐標: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ì)心坐標b(x)=(b0,b1,b2)。

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

ytj(1)

ytj(2)

?

?

(36)

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

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)對應的殘差值

4.2 確定最終模型

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

1) 0.01

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

3) 僅當提高維度有顯著的品質(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 總體坐標系中的表達式

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

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

(48)

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

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

5 結(jié) 論

本文采用多元樣條函數(shù)理論對變體飛機氣動力進行建模，以某型折疊機翼變體飛機為例，建立了基于多元樣條函數(shù)模型結(jié)構(gòu)的氣動力方程，給出了詳細的樣條函數(shù)模型結(jié)構(gòu)和估計流程，通過計算得到了該變體飛機的氣動力數(shù)學模型。通過誤差分析發(fā)現(xiàn)，當維度增加時，模型品質(zhì)未必一直增加，維度選取過高可能會造成擬合失真，因此不能單純的通過提高維度來提高模型品質(zhì)，當提高維度得到的模型無法滿足要求時需要考慮重新構(gòu)造三角測量。結(jié)果表明，該方法無需預知變形過程的氣動變化規(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