張立啟 岳承宇 趙永輝

(南京航空航天大學(xué)機(jī)械結(jié)構(gòu)力學(xué)及控制國(guó)家重點(diǎn)實(shí)驗(yàn)室，南京 210016)

引言

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

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

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

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

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

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

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

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

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

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

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

機(jī)翼結(jié)構(gòu)的邊界條件為翼根固支.表1 給出了χ=8°時(shí)結(jié)構(gòu)的前6 階固有頻率，圖3 為相應(yīng)的固有振型.由于機(jī)翼展弦比較大，前6 階模態(tài)中包含面內(nèi)運(yùn)動(dòng)模態(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é)點(diǎn)坐標(biāo)等模型信息.這些數(shù)據(jù)將在自編程序中進(jìn)一步處理生成機(jī)翼的氣動(dòng)彈性模型.在后掠角0°～ 60°的變化范圍內(nèi)，每隔 5 °計(jì)算出機(jī)翼的固有振動(dòng)特性，機(jī)翼的前6 階固有頻率如表2 所示.

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

2 固定后掠角下的氣動(dòng)彈性模型

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

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

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

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

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

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

式中，Φ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)運(yùn)動(dòng)和突風(fēng)擾動(dòng)有關(guān)的廣義空氣動(dòng)力(generalize aerodynamics force，GAF)矩陣.以上各矩陣的具體表達(dá)式為

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

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

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

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

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

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

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

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

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

式中，φi，l-1和 φj，l分別為l-1和l點(diǎn)處的第i個(gè)和第j個(gè)固有振型向量.

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

4 參變氣動(dòng)彈性建模方法的驗(yàn)證

4.1 一致性驗(yàn)證

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

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

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

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

生成GAF 矩陣時(shí)使用了相匹配的結(jié)構(gòu)模態(tài)，然而MS 算法并不保證各參數(shù)點(diǎn)下的RFA 系數(shù)矩陣具有一致性.如圖8 所示，初始的RFA 系數(shù)矩陣D和E在參數(shù)變化時(shí)數(shù)值的跳躍現(xiàn)象非常顯著.計(jì)算結(jié)果顯示離散參數(shù)點(diǎn)處的RFA 擬合結(jié)果是具有足夠精度的，然而跳躍的數(shù)值變化使參數(shù)點(diǎn)之間的插值無(wú)法進(jìn)行.按照3.2 節(jié)給出的方法將系數(shù)矩陣處理為一致的形式后得到的結(jié)果如圖9 所示.可以看出，處理后的系數(shù)矩陣D和E隨著后掠角平滑、連續(xù)變化，具有一致性.需要強(qiáng)調(diào)的是，本文給出的處理方法不改變RFA 擬合結(jié)果，然而能夠?qū)FA 表達(dá)式轉(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 表達(dá)式確保在構(gòu)建狀態(tài)空間模型時(shí)，各當(dāng)?shù)叵到y(tǒng)矩陣能夠隨參數(shù)連續(xù)變化，因此整個(gè)參數(shù)空間內(nèi)任意點(diǎn)處的時(shí)域方程都能通過(guò)對(duì)當(dāng)?shù)叵到y(tǒng)模型的插值得到.為了驗(yàn)證插值模型的準(zhǔn)確性，對(duì)后掠翼在突風(fēng)擾動(dòng)下的氣動(dòng)彈性響應(yīng)進(jìn)行計(jì)算.1-cos 離散突風(fēng)是工程中常用的標(biāo)準(zhǔn)模型，當(dāng)飛行器以速度U∞穿過(guò)突風(fēng)場(chǎng)，突風(fēng)參考點(diǎn)處受到的垂直擾動(dòng)速度為

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

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

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

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

4.2 顫振分析

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

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

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

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

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

4.3 時(shí)域仿真

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

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

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

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

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

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

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

5 結(jié)論

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

(1)變后掠翼的氣動(dòng)彈性建模問(wèn)題涉及參變系統(tǒng)的建模技術(shù).對(duì)已獲得的參變空間內(nèi)的若干當(dāng)?shù)貧鈩?dòng)彈性模型進(jìn)行插值，是參變氣動(dòng)彈性系統(tǒng)建模的高效而實(shí)用的方法.

(2)同一系統(tǒng)的狀態(tài)空間實(shí)現(xiàn)是不唯一的，這導(dǎo)致各當(dāng)?shù)貧鈩?dòng)彈性模型可能存在不一致性，直接對(duì)這些不一致的當(dāng)?shù)貧鈩?dòng)彈性模型插值，將得到新參數(shù)點(diǎn)下錯(cuò)誤的氣動(dòng)彈性模型.因此，對(duì)當(dāng)?shù)啬Ｐ瓦M(jìn)行一致性處理至關(guān)重要.

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

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

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