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

徐孝武, 張煒, 詹浩

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

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

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

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

1 多元樣條函數基礎理論

本節簡要介紹多元樣條函數基礎理論[13-14,19]。

1.1 質心坐標系

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

t=〈v0,v1,…,vn〉

(1)

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

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

(2)

n元單形t中任何一點的質心坐標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 樣條函數和樣條空間

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

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

(5)

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

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

(6)

樣條空間是指在三角測量Γ區域內,給定維度d和連續階數Cr的所有樣條函數s的集合。

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

(7)

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

1.3 B-form

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

多元系數κ定義如下:

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

(9)

κ的1-范數為:

|κ|=κ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)

κ的所有組合數量為:

=(d+n)!n!d!

(13)

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

(14)

B-系數可以描述為以下矢量形式:

∈R

(15)

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

∈R

(16)

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

B-系數的全局矢量c描述如下:

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

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

…

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

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

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

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

Y=Xc+ε∈RN×1

(22)

1.4 連續條件

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

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

(23)

定義H為平滑矩陣,可以得到全局區域內的連續條件方程:

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

(24)

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

2 基于多元樣條函數的氣動力模型

2.1 幾何模型

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

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

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

2.2 氣動力模型

研究結果表明,在變形速率不大時,變體飛機的氣動力可忽略非定常效應,按準定常計算[8,20]。變體飛機的氣動力僅受飛機當前的氣動外形和飛行狀態影響。穩定軸系繞自身橫軸轉動α(迎角)角度與體軸系重合,穩定軸系上氣動力和力矩系數分別是:升力系數CL,阻力系數CD,側力系數CY,滾轉力矩系數Cl,俯仰力矩系數Cm,偏航力矩系數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為全機機翼參考翼展。

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

CYur=Fu,y(QSw)

Clur=Lu(QSwbw)

Cnur=Nu(QSwbw)

(26)

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

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)

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

2.3 多元樣條函數表達形式

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

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)

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

CL(i)=

c12

c13

(29)

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

3 樣條函數模型估計

3.1 樣條函數模型估計流程

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

圖2 樣條函數模型估計流程圖

3.2 模型結構選取

1) 維度選擇

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

2) 構造三角測量;

圖3 三角測量(Γ4)

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

3) 確定模型空間結構

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

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-系數ctj為:

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

由(18)式可得B-系數的全局矢量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。三角形區域分別為:t1=〈v0v1v4〉,t2=〈v1v2v4〉,t3=〈v0v3v4〉,t4=〈v2v3v4〉。由(3)和(4)式容易得到每個三角形區域內任一點x=(μ1,μ2)的質心坐標b(x)=(b0,b1,b2)。

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

ytj(1)

ytj(2)

?

?

(36)

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

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

(37)

于是得到樣條函數空間內所有數據的線性回歸模型如下:

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-系數估計

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

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

(40)

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

H0+XTY

0=C1C2

C3C4XTY

0

(41)

可由下式計算出B-系數:

=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 最終模型結構確定

4.1 模型驗證

1) 誤差分析理論

χval為需要驗證的數據集合,表示為

χ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) 計算結果分析

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

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

4.2 確定最終模型

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

1) 0.01

2) RMSrelε<0.01時,模型品質視為優異,不再提高維度;

3) 僅當提高維度有顯著的品質提升時才選擇更高的維度。

所有氣動參數的最終估計結果如下:

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

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

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

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

表3 模型誤差分析結果

Cmδe擬合數據圖示如下:

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

4.3 總體坐標系中的表達式

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

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

(48)

對應各參數值見表4。

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

5 結 論

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

參考文獻：

[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