張立啟 岳承宇 趙永輝

(南京航空航天大學機械結(jié)構(gòu)力學及控制國家重點實驗室，南京 210016)

引言

固定翼飛機一直面臨著高速和低速飛行性能要求的矛盾:后掠翼可以有效降低激波阻力，但在亞音速狀態(tài)下飛行效率低;平直機翼在低速飛行時具有較好的飛行性能，但在跨音速區(qū)域，受到的激波阻力卻急劇增長.變后掠機翼可同時兼顧高低速兩個方面的要求.從20 世紀30 年代以來，變后掠翼飛機經(jīng)歷了長達近百年的發(fā)展.早期的Pterodactyl IV 驗證機、F-14 戰(zhàn)斗機、Tornado 戰(zhàn)斗機和Tu-160 轟炸機均實現(xiàn)了主動改變機翼后掠角而帶來的飛行性能提升.然而，增加飛機重量和結(jié)構(gòu)復雜度的缺點使變后掠翼在20 世紀70 年代后的發(fā)展與應(yīng)用陷入停滯.2003 年，美國國防部先進研究計劃局(DAPRA)啟動“變體飛行器結(jié)構(gòu)(MAS) ”研究計劃，重新掀起了變體機翼的研究熱潮[1-3].受此計劃的鼓舞，一些關(guān)于變后掠翼的新概念性設(shè)計浮現(xiàn)出來，比如美國佛羅里達大學研制的獨立多關(guān)節(jié)變后掠翼，新一代航空技術(shù)公司研制的驗證無人機MFX-1 和MFX-2.與早期研究不同，當前變后掠翼的研究重點傾向于高性能材料(比如柔性蒙皮)和驅(qū)動裝置(比如智能結(jié)構(gòu)、柔順機構(gòu))的設(shè)計和使用.

目前，固定翼下的氣動彈性系統(tǒng)建模已經(jīng)成熟，ZAERO 軟件和NASTRAN 中的氣動彈性分析模塊均可對固定翼飛行器進行顫振速度的計算和動態(tài)響應(yīng)分析.而在可變后掠機翼的氣動彈性建模和分析方面的研究較少.飛行器在變體過程中，系統(tǒng)的動力學特性隨著變體控制參數(shù)的變化而變化，相應(yīng)的氣動彈性方程本質(zhì)上是隨參數(shù)時變的.對于參數(shù)慢變情形，可采用時間凍結(jié)法，獲得某一變體參數(shù)下的氣動彈性模型.

在變體飛機方面，Zhao 和Hu[4]和Huang 等[5]提出了一種利用子結(jié)構(gòu)綜合和偶極子網(wǎng)格法構(gòu)建折疊翼的參數(shù)化氣動彈性建模方法，可對不同構(gòu)型下的折疊翼實現(xiàn)快速的顫振分析，分析結(jié)果發(fā)現(xiàn)對于不同折疊角下的機翼其氣動彈性特征將發(fā)生顯著改變.Hu 等[6]在此基礎(chǔ)上采用Krigin 代理模型技術(shù)對不同后掠角下的空氣動力影響矩陣進行插值，獲得了參數(shù)化的氣動彈性模型并研究了折疊翼變形過程中的時變氣動彈性響應(yīng).Xu 等[7]利用ADAMS 對折疊翼進行柔性多體動力學建模，將不同后掠角下的非定常氣動力進行插值后與結(jié)構(gòu)進行耦合，計算分析了折疊翼快速展開過程中的瞬態(tài)響應(yīng).詹玖榆等[8]利用流形插值技術(shù)[9]對模態(tài)坐標下不同構(gòu)型折疊翼的剛度矩陣、質(zhì)量矩陣和振型矩陣進行插值，建立了參數(shù)化的結(jié)構(gòu)動力學模型，并將其與偶極子網(wǎng)格法相結(jié)合計算了折疊翼的氣動彈性行為.

對可變后掠翼，大量的研究集中在最優(yōu)變形規(guī)律以及變后掠翼氣動布局等方面，對于可變后掠翼的參變氣動彈性系統(tǒng)建模方法和慢時變氣動彈性方面的研究依然匱乏，變后掠過程中準確而高效地計算氣動彈性響應(yīng)仍是一項具有挑戰(zhàn)性的工作.由于變后掠翼在結(jié)構(gòu)變體過程中參數(shù)是變化的，使得參變氣動彈性系統(tǒng)的建模存在一定困難.線性參變系統(tǒng)(linear parameter varying，LPV)[10]模型可以用來描述受時變參數(shù)影響的線性系統(tǒng)的動力特征，也可以用一系列當?shù)鼐€性模型來表征非線性系統(tǒng).一般地，LPV 建模方法有兩種:全局建模方法[11-13]和當?shù)亟７椒╗14-15].利用全局模型能夠進行全局參變的研究，即任意激勵下結(jié)構(gòu)持續(xù)參變過程中的動態(tài)響應(yīng)研究.顯然，這一要求比較苛刻，一般很難實現(xiàn)[16];后者則是利用當?shù)亟＜夹g(shù)，采用時間凍結(jié)法，在不同固定參數(shù)下生成當?shù)氐木€性時不變(linear time invariant，LTI)模型，然后對這些LTI 模型進行合理的插值來建立線性參變模型.利用當?shù)亟＜夹g(shù)建立LPV 模型的方法要求所建立的當?shù)豅TI 模型必須具有一致的狀態(tài)空間形式.若直接對系統(tǒng)矩陣進行插值可能帶來很大風險甚至得到錯誤的狀態(tài)空間模型，這是由于同一系統(tǒng)的狀態(tài)空間實現(xiàn)并不唯一，導致不同參數(shù)點處的狀態(tài)空間表達并不具有一致性.上述研究均未能討論和考慮模型的一致性問題.

Paige[17]和Wassink 等[18]著重強調(diào)了不一致的LTI 狀態(tài)空間表達式對于當?shù)豅PV 模型實際應(yīng)用的影響和限制.因此，如何有效的解決當?shù)啬Ｐ椭g的不一致問題，便成為了當?shù)豅PV 模型能否成功運用的關(guān)鍵.Paijmans 等[19]提出了一種依賴于仿射參數(shù)空間實現(xiàn)的LPV 建模技術(shù)，該方法需要對當?shù)啬Ｐ偷臉O點和零點的軌跡進行排序和擬合.Caigy 等[20]提出了一種單輸入單輸出(single input single output，SISO)的LPV 建模方法，該方法需要劃分多個子模型.Caigy 等[21]將當?shù)豅TI 模型劃分為一系列當?shù)刈幽Ｐ蚚19-20]的方法推廣到多輸入多輸出(multiple input multiple output，MIMO)情況，對當?shù)啬Ｐ偷臉O點和零點進行排序，將原始的當?shù)啬Ｐ头纸鉃橐粋€增益與一階和二階子模型狀態(tài)空間串聯(lián)的形式，得到了一致的狀態(tài)空間表達式.Krolick 等[22]提出了一種基于數(shù)據(jù)驅(qū)動的模型降階方法，并探討了不同參數(shù)下模型的一致性問題，發(fā)現(xiàn)利用自回歸(autoregressive exogenous，ARX)方法建立的參數(shù)化模型一致性較好.Goizueta 等[23]利用Krylov 子空間來構(gòu)造線性轉(zhuǎn)換矩陣，解決了LTI 模型不一致問題.

以上方法均是從已知的、參數(shù)空間中離散點處的、不一致的狀態(tài)空間方程出發(fā)，通過對狀態(tài)空間方程的一致性處理來獲得參變系統(tǒng)狀態(tài)空間模型的.本文則從另一種角度出發(fā):即在構(gòu)造參變空間內(nèi)離散點處的狀態(tài)空間模型之前就解決可能存在的不一致性問題.因而后續(xù)所建立的、離散點處的狀態(tài)空間模型就是一致的，參變空間內(nèi)任意點處的狀態(tài)空間模型便可通過對一致的當?shù)啬Ｐ瓦M行插值得到.可通過數(shù)值結(jié)果驗證了本文方法的正確性和有效性.本文采用當?shù)貧鈩訌椥越＜夹g(shù)，深入研究了如何建立一致的當?shù)貧鈩訌椥誀顟B(tài)方程的可靠方法，以期為MIMO 的參變氣動彈性系統(tǒng)建模提供參考數(shù)據(jù).

1 變后掠翼的參數(shù)化結(jié)構(gòu)動力學建模

變后掠翼如圖1 所示.機翼采用旋轉(zhuǎn)-剪切變后掠設(shè)計，翼肋在變后掠過程中保持順氣流翼型不變.機翼的結(jié)構(gòu)設(shè)計參數(shù)為:翼根弦長758.3 mm，翼尖弦長283.3 mm，平均弦長346.4 mm，半展長1718.0 mm，機翼面積0.595 7 m2，重量20.23 kg.機翼的前緣后掠角χ可在 0°(完全平直狀態(tài))到 60 °(大后掠角狀態(tài))范圍內(nèi)變化.

圖1 變后掠翼的幾何模型Fig.1 Geometric model of the variable-sweep wing

后掠翼的運動方程依賴于后掠角χ這一控制參數(shù).常規(guī)有限元建模方法只適用于固定構(gòu)型情形.如需獲得不同后掠角下的有限元模型則需要進行大量的重復性建模工作，過程十分繁瑣，效率低下.本文采用如下參數(shù)化自動建模方案:僅在χ=8 °下建立原始結(jié)構(gòu)的有限元模型(如圖2 所示)，使用自編程序讀取該結(jié)構(gòu)結(jié)點和單元數(shù)據(jù);給定任意后掠角后，該角度下的有限元結(jié)點位置可由原始模型中的結(jié)點位置經(jīng)過坐標旋轉(zhuǎn)獲得，進而直接生成BDF 數(shù)據(jù)文件并提交至NASTRAN，自動進行結(jié)構(gòu)固有振動分析，為參變氣動彈性系統(tǒng)的建模提供結(jié)構(gòu)固有振動特性數(shù)據(jù).如圖2 所示，機翼共劃分為 3 369 個單元，結(jié)構(gòu)采用板、桿、梁和集中質(zhì)量單元構(gòu)建，翼梁和翼肋采用板單元和桿單元相結(jié)合的方式建模，蒙皮采用板單元.

圖2 變后掠翼的有限元模型(χ =8°)Fig.2 Finite element model of the variable-sweep wing (χ =8°)

機翼結(jié)構(gòu)的邊界條件為翼根固支.表1 給出了χ=8°時結(jié)構(gòu)的前6 階固有頻率，圖3 為相應(yīng)的固有振型.由于機翼展弦比較大，前6 階模態(tài)中包含面內(nèi)運動模態(tài)(第2 階).

表1 變后掠翼的前6 階固有模態(tài)(χ =8°)Table 1 The first six natural modes (χ =8°)

圖3 變后掠翼的前6 階固有振型(χ=8°)Fig.3 The first six mode shapes of the variable-sweep wing (χ=8°)

根據(jù)前述的結(jié)構(gòu)參數(shù)化建模方案，任意后掠角下的有限元模型都可以方便地利用NASTRAN 生成.建模的結(jié)果包含結(jié)構(gòu)整體質(zhì)量矩陣、剛度矩陣、固有模態(tài)矩陣以及結(jié)點坐標等模型信息.這些數(shù)據(jù)將在自編程序中進一步處理生成機翼的氣動彈性模型.在后掠角0°～ 60°的變化范圍內(nèi)，每隔 5 °計算出機翼的固有振動特性，機翼的前6 階固有頻率如表2 所示.

表2 不同后掠角下機翼的前6 階固有頻率Table 2 The first six order natural frequencies of the wing at different swept angles

2 固定后掠角下的氣動彈性模型

采用基于偶極子網(wǎng)格法(doublet lattice method，DLM)的非定常氣動力模型[24].DLM 將機翼視為無厚度升力面，并將升力面劃分為若干個與來流方向平行的空氣動力網(wǎng)格，如圖4 所示.與結(jié)構(gòu)有限元模型的參數(shù)化建模方法類似，在建立非定常氣動模型時，僅在χ=8 °下劃分氣動網(wǎng)格，共計799 個網(wǎng)格單元，其他任意后掠角下的氣動網(wǎng)格均可根據(jù)原始網(wǎng)格通過坐標旋轉(zhuǎn)獲得.

圖4 變后掠翼在不同后掠角下的空氣動力網(wǎng)格Fig.4 Aerodynamic boxes on the lifting surface of the variable-sweep wing

根據(jù)DLM 模型計算得到升力面的空氣動力影響系數(shù)(aerodynamic influence coefficient，AIC)矩陣，根據(jù)非穿透邊界條件建立簡諧運動情況下全部氣動網(wǎng)格上的壓力系數(shù)與網(wǎng)格控制點處無量綱下洗速度之間的關(guān)系.然后通過樣條插值將非定常氣動力等效到結(jié)構(gòu)結(jié)點上.變后掠翼的氣動彈性方程可寫為如下形式

式中，xs為結(jié)點位移向量，Ms，Cs和Ks分別為質(zhì)量矩陣、阻尼矩陣和剛度矩陣.fs，fg分別為結(jié)構(gòu)運動與突風擾動引起的非定常氣動力，具有如下形式[25]

式中，U∞為來流流速度，q∞為來流動壓，Gks為樣條矩陣.Skj為氣動網(wǎng)格積分矩陣，Djk為物質(zhì)導數(shù)矩陣. Φjg為突風模態(tài)向量，包含簡諧突風激勵對各氣動網(wǎng)格的作用比值.Aaic為AIC 矩陣，與減縮頻率k與馬赫數(shù)M∞相關(guān).wg為參考點處的突風速度擾動.

引入模態(tài)坐標變換xs=Φhξ，將式(1)轉(zhuǎn)換至模態(tài)坐標中，表示如下

式中，Φh為結(jié)構(gòu)模態(tài)矩陣，ξ 為相應(yīng)的結(jié)構(gòu)模態(tài)位移.Mhh，Chh和Khh分別為模態(tài)質(zhì)量、模態(tài)阻尼和模態(tài)剛度矩陣.Qhh，Qhg分別與結(jié)構(gòu)運動和突風擾動有關(guān)的廣義空氣動力(generalize aerodynamics force，GAF)矩陣.以上各矩陣的具體表達式為

式(4)中的GAF 矩陣為頻域形式，如果進行時域氣動彈性分析，可以通過有理函數(shù)擬合(rational functional approximation，RFA)技術(shù)將GAF 矩陣轉(zhuǎn)換至拉氏域，然后再轉(zhuǎn)化到時域.然而，當前固定后掠角下的氣動彈性模型并不適合在時變角度下做仿真計算，因為對每個時間步，在線計算結(jié)構(gòu)模態(tài)矩陣和非定常氣動力矩陣是非常低效的.下一節(jié)將提出一種基于當?shù)貧鈩訌椥阅Ｐ偷牟逯到＜夹g(shù)，其優(yōu)點是能夠快速獲得參數(shù)空間內(nèi)任意一點的氣動彈性模型.

3 基于當?shù)啬Ｐ筒逯档膮⒆兿到y(tǒng)建模方法

變后掠翼的氣動彈性分析和響應(yīng)計算均可在LPV 系統(tǒng)的框架內(nèi)進行，當?shù)亟＜夹g(shù)是構(gòu)造LPV系統(tǒng)的一種實用且高效的建模手段[26].首先對參數(shù)變化空間進行離散處理，然后在各離散參數(shù)點處分別建立一系列線性時不變(LTI)系統(tǒng)模型，最后對這些LTI 模型進行插值，從而快速獲得參數(shù)空間內(nèi)任意一點處的LTI 模型.在參數(shù)慢變的假設(shè)下，LTI 模型中缺少的動態(tài)參數(shù)依賴部分可以忽略不計[27]，本文采用此假設(shè).將插值步驟融入仿真過程，即可對可變后掠翼進行慢時變響應(yīng)計算.

然而，使用當?shù)亟＜夹g(shù)時，各當?shù)啬Ｐ捅仨毷孪忍幚頌橐恢碌男问?，否則將得到錯誤的插值后系統(tǒng).考察式(4)給出的氣動彈性方程，各矩陣均與結(jié)構(gòu)模態(tài)矩陣 Φh有關(guān)，而 Φh是隨后掠角變化的.一方面，對于參變系統(tǒng)來說，隨著參數(shù)的變化，可能存在所謂的“模態(tài)交叉”現(xiàn)象，此時交叉點前后同一模態(tài)分支(相似的振動形態(tài))的階次將發(fā)生變化，具體表現(xiàn)為矩陣 Φh的數(shù)值上的跳躍，此時不經(jīng)處理直接按照默認模態(tài)階次進行插值是錯誤的;另一方面，RFA 的近似過程中需迭代求解非線性最小二乘問題，算法具有數(shù)值不穩(wěn)定的特點，非常容易造成連續(xù)變化的后掠角對應(yīng)不連續(xù)的擬合系數(shù)矩陣的情況，在這種情形下生成的當?shù)啬Ｐ屯瑯硬荒苓M行直接插值.本文分別針對結(jié)構(gòu)和氣動這兩方面可能存在的不一致性問題給出了相應(yīng)的一致性處理方法.該方法能夠保證在參數(shù)空間內(nèi)生成一致的當?shù)貧鈩訌椥阅Ｐ?

3.1 模態(tài)匹配和對齊

為獲得后掠翼模態(tài)數(shù)據(jù)，在nt個離散的后掠角χ1，χ2，···，χnt處建立后掠翼的結(jié)構(gòu)動力學模型，求解以下廣義特征值問題

式中，Ωh，l為結(jié)構(gòu)固有頻率構(gòu)成的對角陣，矩陣Φh，l的列向量為按模態(tài)質(zhì)量歸一化的結(jié)構(gòu)固有振型.

一般情況下，特征值和特征向量通常簡單地按照大小排列順序.當存在模態(tài)交叉時，同一階次的模態(tài)振型由于模態(tài)性質(zhì)不同而不能直接用來插值.插值應(yīng)該在同一條模態(tài)分支上進行，而同一條模態(tài)分支上的模態(tài)具有相似的振動形態(tài)，因此對相鄰參數(shù)點下的固有頻率和振型進行比較，有望將 Φh，l中的振型調(diào)整為匹配的順序.考察兩個相鄰的后掠角χl-1和χl，定義如下距離度量來衡量相鄰的兩個角度下模態(tài)之間的接近程度

式中，ωi，l-1和 ωj，l分別為l-1和l點處的第i個和第j個固有頻率，nh為保留的結(jié)構(gòu)模態(tài)數(shù)量，MAC為模態(tài)置信準則，按照下式計算

式中，φi，l-1和 φj，l分別為l-1和l點處的第i個和第j個固有振型向量.

距離度量式(9)使用了兩個固有頻率之間的直線距離，并用MAC 加權(quán).MAC 在區(qū)間[0，1]之間取值，其值越大則表明固有振型之間的線性相關(guān)性越高.假設(shè)l-1 點處的模態(tài)數(shù)據(jù)具有正確的順序，模態(tài)匹配的任務(wù)是將l點處的模態(tài)數(shù)據(jù)與l-1 點進行配對，正確的配對應(yīng)使總的距離度量最小.基于圖論，將待考察的模態(tài)劃分成以下兩組頂點

定義邊集El=Vl-1×Vl，即任意兩個Vl-1和Vl中的頂點都有唯一的邊相連，從而構(gòu)成一個完全二分圖Gl=(Vl-1，Vl;El).距離度量式(9) 給出了每條邊對應(yīng)的費用.模態(tài)匹配問題可轉(zhuǎn)換為尋找圖Gl的一個完美匹配使費用總和最小，這個問題又稱為線性求和分配問題(linear sum assignment problem，LSAP)[28].匈牙利算法是求解LSAP 的標準方法，該算法由Kuhn 首先提出.目前廣泛采用Jonker 和Volgennant 開發(fā)的基于最短增廣路徑技術(shù)的改進算法[29]，使線性求和分配問題的求解達到O() 的最壞情況時間復雜度.

為了保證數(shù)值連續(xù)變化，固有振型向量的正負符號問題也應(yīng)順帶處理.將重新排列后的固有振型向量調(diào)整為ci，jφj，l，其中符號系數(shù)ci，j按照下式計算

通過模態(tài)匹配與對齊，能夠保證式(4)中建立的氣動彈性模型在全部后掠角下都有一致的形式.

3.2 一致的有理函數(shù)表達式

另一可能引起模型不一致的因素是非定常氣動力的有理函數(shù)近似.為獲得時域的氣動彈性方程，需要通過RFA 技術(shù)將頻域氣動力擴展至Laplace 域.RFA 技術(shù)利用一系列減縮頻率下的GAF 數(shù)據(jù)擬合特定形式的有理函數(shù)，其中最常用是下式給出的最小狀態(tài)近似(minimum state approximation，MSA)

式中，p為復減縮頻率，Q為GAF 矩陣Qhh和Qhg的組合，R為由空氣動力滯后根組成的對角陣.其余矩陣通過迭代求解非線性最小二乘問題獲得.雖然求解過程得到了確定的有理函數(shù)表達式，但是滯后項中D和E矩陣是不唯一的.一般地，對于任意可逆矩陣T，下式始終成立

在建立可變后掠角氣動彈性系統(tǒng)模型時，擬合后的系數(shù)矩陣的不唯一性會造成當?shù)貧鈩恿δＰ椭g的不一致，從而導致插值得到的系統(tǒng)是錯誤的.由于R矩陣事先給定且為對角陣，為保證R矩陣不變，易知T矩陣也必須限制為對角陣.因此，有理函數(shù)擬合的不唯一性表現(xiàn)在對D和E矩陣的列向量和行向量的任意縮放.假設(shè)R矩陣不隨參數(shù)變化或者隨參數(shù)連續(xù)變化，那么可以根據(jù)相鄰參數(shù)下的有理函數(shù)擬合結(jié)果構(gòu)造出唯一的縮放矩陣T，從而得到唯一且連續(xù)的有理函數(shù)表達式.

考慮相鄰兩個后掠角χl-1和χl，以l-1點處的Dl矩陣為標準參照，l點處的對角陣Tl根據(jù)下式進行構(gòu)造

式中，na為空氣動力滯后根的數(shù)量，對角元ti，l的計算方法如下

式中，di，l-1和di，l分別為矩陣Dl-1和Dl的第i列.得到Ti后根據(jù)式(14)調(diào)整矩陣Dl和El，調(diào)整后的矩陣Dl與參照矩陣Dl-1的列向量有相同的2 范數(shù)和一致的方向.

如后掠角和馬赫數(shù)均為可變參數(shù)，此時RFA 原始數(shù)據(jù)在二維參數(shù)網(wǎng)格上生成.這種情況下一致性處理方法應(yīng)該是:首先固定第一個馬赫數(shù)，按照后掠角順序依次調(diào)整系數(shù)矩陣;隨后固定各個后掠角，再按照馬赫數(shù)的順序使用相同的方法依次調(diào)整系數(shù)矩陣.

經(jīng)過模態(tài)匹配的氣動彈性方程式(4)與使用一致RFA 表達的式(13)相結(jié)合，再按照傳統(tǒng)氣動伺服彈性系統(tǒng)的構(gòu)造方式，最終生成如下時域狀態(tài)空間系統(tǒng)方程[25]

式中，參數(shù)向量 ρ 包括后掠角χ和飛行參數(shù).狀態(tài)變量xae包含結(jié)構(gòu)模態(tài)位移 ξ 和模態(tài)速度，以及有理函數(shù)近似的引入的空氣動力狀態(tài)變量.突風擾動為突風速度wg及其時間導數(shù)構(gòu)成的向量，多種突風模型可用于生成標準化的突風擾動信號[30].系統(tǒng)輸出yae可以是結(jié)構(gòu)某監(jiān)控結(jié)點的響應(yīng)或內(nèi)載荷響應(yīng)，輸出矩陣以及狀態(tài)空間矩陣具體的構(gòu)造方法參考文獻[31].

利用本文給出的模態(tài)匹配和RFA 系數(shù)處理方法，可解決變后掠翼當?shù)貧鈩訌椥阅Ｐ椭锌赡艹霈F(xiàn)的不一致性問題.不論是后掠角還是飛行參數(shù)的改變，各當?shù)叵到y(tǒng)矩陣都能保證隨參數(shù)連續(xù)變化.圖5 總結(jié)了本文提出了基于當?shù)亟＜夹g(shù)的變后掠翼參變氣動彈性系統(tǒng)的建模流程.對一致的當?shù)啬Ｐ瓦M行插值，能夠快速獲取任意后掠角下的狀態(tài)空間矩陣，大大提高時變氣動彈性分析與響應(yīng)計算的效率.

圖5 變后掠翼的參數(shù)化氣動彈性系統(tǒng)建模流程圖Fig.5 The flow chat of the parameter-varying aeroelastic modeling for a variable-sweep wing

4 參變氣動彈性建模方法的驗證

4.1 一致性驗證

為了利用當?shù)亟＜夹g(shù)建立變后掠翼的參變氣動彈性模型，首先在0°～ 60°的后掠角范圍內(nèi)每間隔5°選取13 個離散的后掠角參數(shù)點，在這些參數(shù)點上分別建立當?shù)貧鈩訌椥阅Ｐ?初始的氣動彈性模型包含式(4)中的結(jié)構(gòu)矩陣和GAF 矩陣，由于不同后掠角之間的結(jié)構(gòu)模態(tài)可能不匹配，因此原始當?shù)叵到y(tǒng)的一致性無法保證.根據(jù)本文一致性處理方法，選取前6 階結(jié)構(gòu)模態(tài)，隨后按照3.1 節(jié)給出的算法進行結(jié)構(gòu)模態(tài)匹配.

圖6 給出了經(jīng)過模態(tài)匹配與排序后，變后掠翼的前6 階固有頻率隨后掠角的變化規(guī)律.可以看出，在 55°～ 60 °后掠角之間，第5 階和第6 階結(jié)構(gòu)模態(tài)發(fā)生了交叉現(xiàn)象.模態(tài)匹配過程中，這兩個角度之間的具體MAC 數(shù)值如圖7 所示，圖中清晰地顯示出在 60°下的第5 階和第6 階模態(tài)應(yīng)與 55 °下的第6 階和第5 階模態(tài)相匹配.模態(tài)匹配算法給出了正確的匹配結(jié)果，確保每一個模態(tài)分支都隨參數(shù)連續(xù)變化.

圖6 前6 階固有頻率隨后掠角的變化規(guī)律Fig.6 The first six natural frequencies vs.the swept angles

圖7 MAC 值(55°和60°后掠角)Fig.7 The MAC values for 55° and 60° swept angles

為了建立時域狀態(tài)空間模型，接下來對式(4)中的GAF 矩陣進行有理函數(shù)擬合.假設(shè)來流為不可壓流動，GAF 矩陣在0.0～ 1.5 范圍內(nèi)的16 個減縮頻率處計算.使用MS 方法擬合得到式(13)所示的RFA表達式，氣動力滯后根由以下經(jīng)驗公式[32]得到

式中，kmax為減縮頻率的最大值.

生成GAF 矩陣時使用了相匹配的結(jié)構(gòu)模態(tài)，然而MS 算法并不保證各參數(shù)點下的RFA 系數(shù)矩陣具有一致性.如圖8 所示，初始的RFA 系數(shù)矩陣D和E在參數(shù)變化時數(shù)值的跳躍現(xiàn)象非常顯著.計算結(jié)果顯示離散參數(shù)點處的RFA 擬合結(jié)果是具有足夠精度的，然而跳躍的數(shù)值變化使參數(shù)點之間的插值無法進行.按照3.2 節(jié)給出的方法將系數(shù)矩陣處理為一致的形式后得到的結(jié)果如圖9 所示.可以看出，處理后的系數(shù)矩陣D和E隨著后掠角平滑、連續(xù)變化，具有一致性.需要強調(diào)的是，本文給出的處理方法不改變RFA 擬合結(jié)果，然而能夠?qū)FA 表達式轉(zhuǎn)化成為適合插值的形式.

圖8 不一致的RFA 系數(shù)矩陣隨后掠角的變化Fig.8 The incoherent RFA coefficients vs.the swept angles

圖9 一致性處理后RFA 系數(shù)矩陣隨后掠角的變化Fig.9 The coherent RFA coefficients vs.the swept angles

匹配的結(jié)構(gòu)模態(tài)和一致的RFA 表達式確保在構(gòu)建狀態(tài)空間模型時，各當?shù)叵到y(tǒng)矩陣能夠隨參數(shù)連續(xù)變化，因此整個參數(shù)空間內(nèi)任意點處的時域方程都能通過對當?shù)叵到y(tǒng)模型的插值得到.為了驗證插值模型的準確性，對后掠翼在突風擾動下的氣動彈性響應(yīng)進行計算.1-cos 離散突風是工程中常用的標準模型，當飛行器以速度U∞穿過突風場，突風參考點處受到的垂直擾動速度為

式中，為設(shè)計突風速度幅值，Lg為突風尺度，fg=U∞/Lg為突風激勵頻率.

在海平面大氣環(huán)境下，來流速度為 120 m/s，1-cos突風頻率為 3Hz，突風速度幅值為5 m/s.圖10 和圖11分別給出了后掠角為 28°和 48 °時后掠翼的時域動態(tài)響應(yīng)，圖中展示了翼梢加速度和翼根彎矩的計算結(jié)果. 28°后掠角下的系統(tǒng)模型通過對 25°和 30 °的當?shù)啬Ｐ途€性插值得到，48°的模型通過對 45°和 50 °的當?shù)啬Ｐ筒逯档玫?，圖中一并給出了各后掠角下精確的建模結(jié)果作為對比.結(jié)果顯示，由一致系統(tǒng)矩陣得到的插值模型與精確模型計算出的結(jié)果相同，驗證了對一致當?shù)啬Ｐ瓦M行插值的方法是準確可靠的，而根據(jù)不一致系統(tǒng)矩陣建立的插值模型誤差很大，動態(tài)響應(yīng)與精確結(jié)果相比相去甚遠.在 48 °時，不一致模型的響應(yīng)是發(fā)散的，插值得到了不穩(wěn)定的氣動彈性模型.

圖10 后掠翼的1-cos 離散突風響應(yīng)(χ=28°)Fig.10 1-cos gust response of wing at the swept angleχ=28°

圖11 后掠翼的1-cos 離散突風響應(yīng)(χ=48°)Fig.11 1-cos gust response of wing at the swept angleχ=48°

4.2 顫振分析

圖12 給出了固定馬赫數(shù)、不同后掠角下的顫振速度分布，隨著后掠角的增加，機翼的顫振速度不斷增加.當機翼后掠角大于 18 °時，機翼顫振速度迅速增加，而當后掠角超過21°時顫振速度增加減緩.顫振速度在20°和21°之間發(fā)生了跳躍現(xiàn)象.利用p-k法，計算了機翼在20°和21°后掠角下的顫振特性，如圖13 所示.發(fā)現(xiàn)當后掠角為20°時，隨著來流速度的增加，第4 階模態(tài)分支率先達到零點，系統(tǒng)發(fā)生顫振;而當后掠角增加到21°時，第1 階模態(tài)分支首先穿越零點，失穩(wěn)模態(tài)轉(zhuǎn)為第1 階模態(tài).失穩(wěn)模態(tài)的轉(zhuǎn)換導致了顫振速度的跳躍現(xiàn)象.

圖12 可變后掠翼顫振速度隨后掠角的變化Fig.12 Flutter speed vs.swept angle

圖13 不同后掠角下機翼的顫振分析Fig.13 Flutter analyses of the wing at different swept angle

為研究不同后掠角下系統(tǒng)特征值的變化規(guī)律，進行了固定馬赫數(shù)(M∞=0.6)和來流速度(U∞=200 m/s)下不同后掠角處氣動彈性系統(tǒng)的根軌跡分析，如圖14 所示，圖中箭頭的方向是后掠角增加的方向.隨著后掠角的增加，系統(tǒng)由失穩(wěn)狀態(tài)逐漸演變?yōu)榉€(wěn)定狀態(tài)，這是結(jié)構(gòu)和氣動相互耦合作用的結(jié)果.由圖6可以看出，當機翼的后掠角小于20°時，結(jié)構(gòu)的第3 階模態(tài)(第二階彎曲模態(tài))與第4 階模態(tài)(第一階扭轉(zhuǎn)模態(tài))的振動頻率比較接近，容易誘發(fā)顫振.事實上，從圖14 不同后掠角下氣動彈性系統(tǒng)的根軌跡中也能看出，當后掠角小于18°時，系統(tǒng)失穩(wěn)均發(fā)生在第4 階模態(tài)處，而當后掠角較大時，結(jié)構(gòu)的第3 階與第4 階固有頻率相距較遠，機翼不容易發(fā)生顫振.

圖14 可變后掠翼在不同后掠角下的根軌跡Fig.14 Root locus of the wing at various swept angles

4.3 時域仿真

為進一步分析可變后掠翼在變后掠過程中的時變動態(tài)特性，假定機翼后掠角的變化遵循如下規(guī)律

式中，χm=60°為最大后掠角，Tm=12 s 為從初始狀態(tài)到最大后掠角所需要的總時間.

以Dryden 連續(xù)突風作為擾動激勵，圖15 給出了飛行馬赫數(shù)為0.35，巡航速度為120 m/s 時，機翼連續(xù)改變后掠角過程中翼尖加速度和翼根彎矩的時間歷程，其中連續(xù)突風尺度Lw=267 m，均方根值σw=2 m/s.

圖15 機翼變后掠過程中的突風響應(yīng)Fig.15 Gust responses of the variable-sweep wing during the morphing process

由圖12 可以看出可變后掠翼的最小顫振速度發(fā)生在后掠角為 0 °時，其值為123.95 m/s.圖15 可以看出，機翼以120 m/s 巡航速度飛行時，其速度已接近顫振速度，整個氣動彈性系統(tǒng)的阻尼水平很低，所以開始遇到突風時，結(jié)構(gòu)響應(yīng)較大，但隨著后掠角的不斷增加，系統(tǒng)的阻尼水平增加，翼尖加速度與翼根彎矩響應(yīng)變減弱.

選取了飛行馬赫數(shù)為0.35 以及130 m/s 的巡航速度來研究變后掠翼變形過程中的時變氣動彈性響應(yīng).如圖16 所示.當后掠角小于8.27°時機翼處于氣動彈性失穩(wěn)狀態(tài)，在一個較小的擾動下結(jié)構(gòu)便發(fā)生了諧振蕩，且不斷從來流中獲取能量導致振幅不斷變大;而當后掠角大于8.27°時機翼結(jié)構(gòu)的顫振速度增加，穩(wěn)定性得到提升，結(jié)構(gòu)進入氣動彈性穩(wěn)定狀態(tài)，振幅不斷下降.

圖16 后掠翼變形過程中的氣動彈性響應(yīng)Fig.16 Aeroelastic responses of the variable-sweep wing during the morphing process

5 結(jié)論

本文針對可變后掠翼，提出了一種參變氣動彈性建模的實用方法.通過研究分析，可以得到如下結(jié)論.

(1)變后掠翼的氣動彈性建模問題涉及參變系統(tǒng)的建模技術(shù).對已獲得的參變空間內(nèi)的若干當?shù)貧鈩訌椥阅Ｐ瓦M行插值，是參變氣動彈性系統(tǒng)建模的高效而實用的方法.

(2)同一系統(tǒng)的狀態(tài)空間實現(xiàn)是不唯一的，這導致各當?shù)貧鈩訌椥阅Ｐ涂赡艽嬖诓灰恢滦裕苯訉@些不一致的當?shù)貧鈩訌椥阅Ｐ筒逯?，將得到新參?shù)點下錯誤的氣動彈性模型.因此，對當?shù)啬Ｐ瓦M行一致性處理至關(guān)重要.

(3)當?shù)貧鈩訌椥阅Ｐ偷牟灰恢滦詠碓从诮Y(jié)構(gòu)和氣動兩方面.本文首先利用模態(tài)跟蹤技術(shù)對結(jié)構(gòu)模態(tài)進行匹配和對齊，解決了結(jié)構(gòu)動力學方面的不一致性問題;在氣動方面，在一致的結(jié)構(gòu)模態(tài)基礎(chǔ)之上建立非定常空氣動力矩陣，然后對RFA 表達式中的系數(shù)矩陣進行合理的縮放處理，解決了空氣動力系數(shù)矩陣的不一致問題.這樣，最終的當?shù)貭顟B(tài)空間氣動彈模型便具有良好的一致性，任意參數(shù)下的氣動彈性模型可通過當?shù)啬Ｐ偷牟逯悼焖佾@得，大大提高了系統(tǒng)穩(wěn)定性分析和慢時變響應(yīng)計算的效率.

(4)本文研究為采用頻域氣動模型的參變系統(tǒng)的氣動彈性建模，提供了一種有效的解決方案.本文方法不僅使適用于變后掠情形，也適用于其他變體方式，具有一定的普適性.

(5)慢時仿真結(jié)果表明，機翼在變后掠過程中，其氣動彈性穩(wěn)定性和響應(yīng)特性均發(fā)生很大變化，這些變化規(guī)律對變后掠翼飛行器的設(shè)計具有重要意義.