基于非線(xiàn)性狀態(tài)空間辨識(shí)的氣動(dòng)彈性模型降階1)

張家銘 楊執(zhí)鈞 黃銳

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

引言

氣動(dòng)彈性力學(xué)主要研究空氣動(dòng)力、彈性力和慣性力相互耦合作用下的動(dòng)力學(xué)問(wèn)題,研究重點(diǎn)是彈性結(jié)構(gòu)在氣流中的顫振穩(wěn)定性[1-3].實(shí)現(xiàn)對(duì)飛行器結(jié)構(gòu)顫振邊界尤其是跨音速飛行條件下的顫振邊界的精確預(yù)計(jì),需要建立精確的氣動(dòng)彈性數(shù)學(xué)模型.傳統(tǒng)的線(xiàn)性化氣動(dòng)力模型難以準(zhǔn)確計(jì)算跨音速飛行時(shí)的非線(xiàn)性空氣動(dòng)力,進(jìn)而導(dǎo)致對(duì)飛行器結(jié)構(gòu)跨音速氣動(dòng)彈性力學(xué)行為的預(yù)測(cè)產(chǎn)生誤差.

計(jì)算流體力學(xué) (computational fluid dynamics,CFD)方法為精確模擬跨音速氣動(dòng)彈性行為提供了高保真分析手段.然而,要獲得精確的仿真結(jié)果,基于CFD 技術(shù)的跨音速氣動(dòng)彈性分析需要消耗大量的計(jì)算資源且難以用于控制律設(shè)計(jì),使得該方法的應(yīng)用受到了限制.為了應(yīng)對(duì)這項(xiàng)挑戰(zhàn),近年來(lái)國(guó)內(nèi)外學(xué)者圍繞如何建立高效且高精度的跨聲速氣動(dòng)彈性降階模型(reduced-order model,ROM)開(kāi)展研究,目的是發(fā)展低維、精確的氣動(dòng)力降階模型來(lái)進(jìn)行跨音速氣動(dòng)彈性分析,從而在數(shù)量級(jí)上減少計(jì)算成本.當(dāng)前的氣動(dòng)力降階模型研究可分為如下兩類(lèi):一類(lèi)是基于流場(chǎng)特征分析方法(如本征正交分解法[4-7]、均衡實(shí)現(xiàn)算法[8]、諧波平衡法[9-10]);另一類(lèi)是基于系統(tǒng)辨識(shí)方法(如遞歸局部線(xiàn)性神經(jīng)模糊模型[11]、Volterra級(jí)數(shù)法[12-13]、Wiener 模型法[14-15]、Kriging 代理法[16-18]、特征系統(tǒng)實(shí)現(xiàn)算法[8](eigensystem realization algorithm,ERA)).相對(duì)于基于流場(chǎng)特征分析方法,基于系統(tǒng)辨識(shí)方法受到了更為廣泛的關(guān)注.例如,西北工業(yè)大學(xué)張偉偉等[19]提出了一種混合降階模型,該模型基于完全耦合的子系統(tǒng),非線(xiàn)性輸出反饋到線(xiàn)性和非線(xiàn)性系統(tǒng)中,從而能夠更精確地?cái)M合非定常流中的物理現(xiàn)象.該團(tuán)隊(duì)還提出了一種基于長(zhǎng)短時(shí)記憶網(wǎng)絡(luò)的非定?？諝鈩?dòng)力ROM,該模型能夠處理大量的訓(xùn)練樣本,且能夠預(yù)測(cè)馬赫數(shù)在較大范圍內(nèi)的氣動(dòng)載荷[20].NASA 蘭利研究中心Silva[21]首次將Volterra 級(jí)數(shù)應(yīng)用到非線(xiàn)性空氣動(dòng)力系統(tǒng)的ROM研究中,實(shí)現(xiàn)了一種基于系統(tǒng)對(duì)脈沖輸入響應(yīng)的直接核識(shí)別方法.為了把Volterra 級(jí)數(shù)拓展應(yīng)用到多輸入/多輸出(multiple-input/multiple-output,MIMO)的情況,Balajewicz 等[22]利用多輸入Volterra 級(jí)數(shù),提出了一種適用于兩自由度氣動(dòng)彈性系統(tǒng)的非線(xiàn)性ROM.盡管如此,MIMO 系統(tǒng)的Volterra 核的辨識(shí)過(guò)程仍非常復(fù)雜且辨識(shí)效率較低.為了解決上述問(wèn)題,本文提出了一種基于非線(xiàn)性狀態(tài)空間辨識(shí)的非線(xiàn)性ROM方法,用來(lái)計(jì)算任意模態(tài)運(yùn)動(dòng)下的非定?？諝鈩?dòng)力響應(yīng).同基于Volterra 級(jí)數(shù)的ROM 辨識(shí)方法相比,該辨識(shí)算法更為簡(jiǎn)單,也能夠應(yīng)用到MIMO 空氣動(dòng)力系統(tǒng)的辨識(shí).同時(shí),由于該非線(xiàn)性ROM 方法是基于二級(jí)訓(xùn)練,可較為精確地預(yù)測(cè)跨音速線(xiàn)性和非線(xiàn)性氣動(dòng)彈性現(xiàn)象,進(jìn)而達(dá)到簡(jiǎn)化物理模型、節(jié)省計(jì)算時(shí)間和計(jì)算負(fù)荷的目的[23].

本工作采用非線(xiàn)性狀態(tài)空間辨識(shí)方法,推導(dǎo)出跨音速氣動(dòng)彈性系統(tǒng)的非線(xiàn)性降階模型,為了驗(yàn)證該非線(xiàn)性ROM 在預(yù)測(cè)跨聲速氣動(dòng)彈性問(wèn)題的精確性,給出了ROM 的辨識(shí)過(guò)程并分別進(jìn)行了氣動(dòng)力預(yù)測(cè)和氣動(dòng)彈性分析,以期為建立適用于跨聲速氣動(dòng)彈性系統(tǒng)的非線(xiàn)性降階模型提供參考依據(jù).

1 非線(xiàn)性狀態(tài)空間模型降階方法

不同于單自由度系統(tǒng),基于CFD 的跨聲速空氣動(dòng)力學(xué)系統(tǒng)是高維、非線(xiàn)性系統(tǒng),而且結(jié)構(gòu)多模態(tài)運(yùn)動(dòng)和廣義氣動(dòng)力輸出是典型的非線(xiàn)性多輸入/多輸出動(dòng)力學(xué)系統(tǒng).為了高效、高精度地預(yù)測(cè)跨音速非定常氣動(dòng)力、跨音速顫振邊界以及極限環(huán)振蕩(limit-cycle oscillations,LCO)現(xiàn)象,本節(jié)基于非線(xiàn)性狀態(tài)空間法建立多輸入/多輸出的跨聲速氣動(dòng)力降階模型,模型降階過(guò)程分為線(xiàn)性狀態(tài)空間模型辨識(shí)和非線(xiàn)性參數(shù)優(yōu)化兩部分.

1.1 線(xiàn)性狀態(tài)空間模型的辨識(shí)

線(xiàn)性系統(tǒng)的辨識(shí)方法有很多種,其中常見(jiàn)的有子空間辨識(shí)方法、ARMA 方法以及基于特征系統(tǒng)實(shí)現(xiàn)算法的辨識(shí)方法等.子空間辨識(shí)方法和ARMA 模型對(duì)訓(xùn)練信號(hào)的設(shè)計(jì)要求較高(需要設(shè)定帶寬和幅值),且?guī)捄头档倪x擇對(duì)結(jié)果影響較大.特征系統(tǒng)實(shí)現(xiàn)算法只需要非定常氣動(dòng)力的單位脈沖響應(yīng)數(shù)據(jù),過(guò)程簡(jiǎn)單且精度較高,因此本文采用基于特征系統(tǒng)實(shí)現(xiàn)算法的辨識(shí)方法來(lái)辨識(shí)線(xiàn)性狀態(tài)空間模型.

線(xiàn)性狀態(tài)空間模型形式如下[24]

其中,xa∈Rn是狀態(tài)向量,u∈Rm是輸入向量,ya∈Rl是輸出向量,n,m,l分別是狀態(tài)變量、輸入變量、輸出變量的維數(shù),狀態(tài)空間參數(shù)(Aa,Ba,Ca,Da)是待辨識(shí)的變量.基于ERA 方法的狀態(tài)空間系統(tǒng)辨識(shí)的過(guò)程如下:

考慮由微分方程表征的動(dòng)力學(xué)系統(tǒng)

該系統(tǒng)中的數(shù)據(jù)可以組成如下Hankel 矩陣

由于矩陣的每一行都與第一行線(xiàn)性相關(guān)(如第二行是第一行的eλ?t倍),該矩陣的秩r=1.類(lèi)似地,矩陣的每一列都與第一列線(xiàn)性相關(guān).在這種情況下,在延遲坐標(biāo)中增廣狀態(tài)空間不會(huì)增大Hankel 矩陣的秩,ROM 將始終只含有一階模態(tài).

對(duì)于形式如公式(1)所示的線(xiàn)性狀態(tài)空間模型,其結(jié)構(gòu)運(yùn)動(dòng)的離散時(shí)間脈沖響應(yīng)如下式給出

此時(shí),空氣動(dòng)力輸出可表示為

由空氣動(dòng)力輸出可分別建立如下形式的Hankel 矩陣H和H′

式中,s和m分別是Hankel 矩陣的行數(shù)和Hankel 矩陣中測(cè)量值的數(shù)目.通過(guò)對(duì)H進(jìn)行奇異值分解(singular value decomposition,SVD),可以得到脈沖響應(yīng)時(shí)間序列的主導(dǎo)時(shí)間模式

其中,r表示對(duì)矩陣做秩為r的截?cái)?U和V的列是特征時(shí)間序列,能夠描述跨越不同時(shí)間尺度的脈沖響應(yīng).通過(guò)對(duì)Hankel 矩陣H和H′進(jìn)行SVD,就可以得到A,B,C,D的低維近似

1.2 非線(xiàn)性部分的優(yōu)化

通過(guò)引入狀態(tài)變量和輸入變量的非線(xiàn)性函數(shù),可得到如下形式的非線(xiàn)性狀態(tài)空間模型

式中,Φ(xa(k))和Φ(u(k))分別是狀態(tài)變量和輸入變量的非線(xiàn)性函數(shù),非線(xiàn)性部分的系數(shù)(Ea,Fa,Ga,Ha)是待辨識(shí)的變量,可通過(guò)廣義優(yōu)化[25](Levenberg-Marquardt,LM)算法來(lái)確定非線(xiàn)性部分的系數(shù).本文選用神經(jīng)網(wǎng)絡(luò)辨識(shí)過(guò)程中最常見(jiàn)的雙曲正切函數(shù)tanh(·)來(lái)表征上述非線(xiàn)性函數(shù),該函數(shù)被認(rèn)為是一種能在線(xiàn)性和非線(xiàn)性行為間建立平衡的嚴(yán)格單調(diào)遞增函數(shù).需要說(shuō)明的是,本文提出的非線(xiàn)性模型與Mannarino 和Dowell 提出的非線(xiàn)性降階模型法相似[26].但在他們的研究中,只考慮了與狀態(tài)矢量有關(guān)的非線(xiàn)性貢獻(xiàn),當(dāng)空氣動(dòng)力非線(xiàn)性較強(qiáng)時(shí),該方法必須通過(guò)增加狀態(tài)空間的階數(shù)來(lái)提高辨識(shí)精度,進(jìn)而導(dǎo)致降階模型的維數(shù)過(guò)高且可能出現(xiàn)過(guò)擬合問(wèn)題.在本文的研究中,通過(guò)增加輸入變量的非線(xiàn)性函數(shù),可在不增加狀態(tài)空間維數(shù)的前提下顯著提高空氣動(dòng)力的辨識(shí)精度.

圖1 所示是非線(xiàn)性狀態(tài)空間模型的辨識(shí)過(guò)程.首先,以單位脈沖響應(yīng)數(shù)據(jù)為訓(xùn)練樣本,通過(guò)ERA 方法辨識(shí)線(xiàn)性狀態(tài)空間參數(shù).其次,以白噪聲激勵(lì)下的氣動(dòng)力輸出數(shù)據(jù)為樣本,通過(guò)LM 算法對(duì)非線(xiàn)性狀態(tài)空間中的非線(xiàn)性參數(shù)進(jìn)行優(yōu)化,進(jìn)而獲得最終的非線(xiàn)性降階模型.

圖1 空氣動(dòng)力系統(tǒng)非線(xiàn)性ROM 的辨識(shí)過(guò)程Fig.1 Identification process of nonlinear ROM for aerodynamic systems

1.3 基于非線(xiàn)性氣動(dòng)力降階模型的氣動(dòng)彈性建模

本節(jié)以BACT (benchmark active control technology)機(jī)翼為研究對(duì)象,給出基于非線(xiàn)性氣動(dòng)力降階模型的氣動(dòng)彈性建模過(guò)程.該機(jī)翼為矩形機(jī)翼,采用NACA 0012 翼型[27],其風(fēng)洞模型裝配了后緣控制面以及上下表面擾流板.本文研究中僅考慮后緣控制面作為氣動(dòng)伺服彈性輸入.

圖2 所示是BACT 機(jī)翼的表面氣動(dòng)網(wǎng)格,其結(jié)構(gòu)動(dòng)力學(xué)參數(shù)如表1 所示.表中M,Sα,Sβ,Iα,Sαβ,Kh,Kα的定義見(jiàn)參考文獻(xiàn)[28].

圖2 BACT 機(jī)翼模型Fig.2 BACT wing model

表1 BACT 機(jī)翼結(jié)構(gòu)參數(shù)Table 1 Structural parameters of BACT wing

采用Lagrange 方程方法和虛功原理,建立如下氣動(dòng)彈性運(yùn)動(dòng)方程

其中,h和α 分別是BACT 機(jī)翼的沉浮位移和俯仰角,β 是控制面偏轉(zhuǎn)角,q∞是飛行動(dòng)壓.由于本文提出的非線(xiàn)性氣動(dòng)力降階模型可用于預(yù)測(cè)考慮控制面輸入的非線(xiàn)性氣動(dòng)力,因此β 是作為輸入量.非定常空氣動(dòng)力Qh,Qα可通過(guò)下式計(jì)算

式中,xae是彈性軸的位置.基于模態(tài)坐標(biāo)變換,可將式(11)所示的氣動(dòng)彈性運(yùn)動(dòng)方程轉(zhuǎn)化為以廣義坐標(biāo)表示的氣動(dòng)彈性運(yùn)動(dòng)方程

式中,φ是振型矩陣,Mξξ,Cξξ,Kξξ分別是廣義質(zhì)量矩陣、廣義阻尼矩陣和廣義剛度矩陣.由于控制面的剛度不是無(wú)窮大,不能認(rèn)為控制面在控制輸入下會(huì)立即做出響應(yīng).在指定的控制輸入下,控制面的運(yùn)動(dòng)可認(rèn)為是一種有阻尼單自由度系統(tǒng)的受迫振動(dòng),其運(yùn)動(dòng)微分方程的數(shù)學(xué)表達(dá)式如下

其中,kβ,ζβ,ωβ,βc分別是增益、阻尼比、固有頻率和控制指令,參數(shù)值見(jiàn)文獻(xiàn)[28].式(13)表達(dá)的氣動(dòng)彈性模型在包含了后緣控制面執(zhí)行機(jī)構(gòu)后可改寫(xiě)為

如式(15)所示,廣義空氣動(dòng)力Qξξ(ξ1,ξ2,β)是以沉浮模態(tài)ξ1、俯仰模態(tài)ξ2和偏轉(zhuǎn)角β 為自變量的非線(xiàn)性函數(shù),因此本文的氣動(dòng)力降階模型中狀態(tài)變量包括沉浮模態(tài)、俯仰模態(tài)和控制面偏轉(zhuǎn)模態(tài).而在氣動(dòng)彈性仿真中,將β 角度置零.基于本文所提出的非線(xiàn)性ROM 可高效計(jì)算廣義氣動(dòng)力.

2 算例驗(yàn)證

本節(jié)通過(guò)識(shí)別基于CFD 的空氣動(dòng)力系統(tǒng)來(lái)研究非線(xiàn)性ROM 在給定馬赫數(shù)下對(duì)線(xiàn)性和非線(xiàn)性氣動(dòng)彈性行為預(yù)測(cè)的精確性.非定?？諝鈩?dòng)力訓(xùn)練樣本由基于雷諾平均Navier-Stokes 方程的開(kāi)源CFD求解器求出.CFD 求解器使用源代碼求解,當(dāng)馬赫數(shù)M=0.71 時(shí),正弦激勵(lì)信號(hào)幅值為0.05,頻率為21.11 Hz,流體的無(wú)量綱時(shí)間步長(zhǎng)?τ=5.0,轉(zhuǎn)化為物理時(shí)間步長(zhǎng)為9.769 4×10?4s.

2.1 非線(xiàn)性降階模型的辨識(shí)

首先,基于CFD 產(chǎn)生的非定?？諝鈩?dòng)力的單位脈沖響應(yīng)數(shù)據(jù),采用特征系統(tǒng)實(shí)現(xiàn)算法對(duì)非線(xiàn)性狀態(tài)空間模型的線(xiàn)性動(dòng)力學(xué)參數(shù)進(jìn)行系統(tǒng)辨識(shí).圖3(a)～圖3(d)分別給出了由狀態(tài)空間辨識(shí)方法和CFD 得到的空氣動(dòng)力單位脈沖響應(yīng),圖中Gij代表對(duì)應(yīng)第j個(gè)模態(tài)輸入對(duì)第i個(gè)廣義氣動(dòng)力輸出的影響.

如前文所述,建立該非線(xiàn)性ROM 的關(guān)鍵在于對(duì)線(xiàn)性動(dòng)力學(xué)部分的辨識(shí).在辨識(shí)過(guò)程中線(xiàn)性部分的系統(tǒng)階數(shù)的選擇對(duì)模型的精度具有較大影響.為了保證仿真中降階模型的收斂性并獲得線(xiàn)性部分的最優(yōu)狀態(tài)空間方程,需要通過(guò)ERA 算法在較大的階數(shù)范圍(最大可達(dá)20 階)內(nèi)進(jìn)行相應(yīng)的參數(shù)估計(jì),識(shí)別每個(gè)ERA 模型的動(dòng)力學(xué)線(xiàn)性部分,并計(jì)算這些已辨識(shí)的線(xiàn)性部分的相對(duì)誤差.真實(shí)信號(hào)和預(yù)估信號(hào)的百分比誤差如式(17)定義,選擇百分比誤差的最小值所對(duì)應(yīng)的系統(tǒng)階數(shù)來(lái)確定系統(tǒng)的最優(yōu)階數(shù),最后確定狀態(tài)空間方程.通過(guò)對(duì)動(dòng)力學(xué)系統(tǒng)線(xiàn)性部分的辨識(shí),得到模型的最優(yōu)階數(shù)r=10.

實(shí)現(xiàn)該方法的另一個(gè)關(guān)鍵問(wèn)題是激勵(lì)函數(shù)的選擇.已有研究表明,濾波高斯白噪聲(filtered Gaussian white noise,FWGN)函數(shù)是這類(lèi)應(yīng)用中的最優(yōu)函數(shù)[14].本文的非線(xiàn)性ROM 辨識(shí)是基于CFD 空氣動(dòng)力系統(tǒng)在FWGN 激勵(lì)下所產(chǎn)生的非定常訓(xùn)練數(shù)據(jù).當(dāng)馬赫數(shù)為0.71 時(shí),BACT 機(jī)翼受非定常模態(tài)激勵(lì)(如圖4 所示的白噪聲激勵(lì))時(shí)會(huì)產(chǎn)生激波振蕩.因此,非線(xiàn)性來(lái)源于激波運(yùn)動(dòng)產(chǎn)生的非線(xiàn)性空氣動(dòng)力.圖4 給出了線(xiàn)性/非線(xiàn)性ROM 預(yù)測(cè)的沉浮和俯仰模態(tài)所對(duì)應(yīng)的廣義空氣動(dòng)力,并與直接CFD 方法進(jìn)行了對(duì)比.如圖所示,線(xiàn)性狀態(tài)空間辨識(shí)方法雖然能夠較好的表達(dá)沉浮和俯仰模態(tài)受FWGN 激勵(lì)時(shí)的空氣動(dòng)力相位特征,但在廣義氣動(dòng)力較大峰值處的辨識(shí)精度要低于本文提出的非線(xiàn)性降階模型.為了驗(yàn)證非線(xiàn)性ROM 預(yù)測(cè)結(jié)構(gòu)模態(tài)受其他激勵(lì)下的空氣動(dòng)力學(xué)特征,本文研究了在M∞=0.71 時(shí)沉浮和俯仰模態(tài)受正弦激勵(lì)的廣義空氣動(dòng)力輸出響應(yīng).如圖5 所示,廣義空氣動(dòng)力的對(duì)比表明了該非線(xiàn)性狀態(tài)空間辨識(shí)方法能夠與直接CFD 計(jì)算結(jié)果相吻合.

圖3 脈沖激勵(lì)下的空氣動(dòng)力比較Fig.3 Comparison of aerodynamic forces under the impulse excitation

圖4 FWGN 激勵(lì)下直接CFD 模型輸出和非線(xiàn)性ROM 輸出對(duì)比Fig.4 Direct CFD model outputs vs nonlinear ROM outputs under the FWGN excitation

圖4 FWGN 激勵(lì)下直接CFD 模型輸出和非線(xiàn)性ROM 輸出對(duì)比(續(xù))Fig.4 Direct CFD model outputs vs nonlinear ROM outputs under the FWGN excitation(continued)

為了將線(xiàn)性與非線(xiàn)性方法預(yù)測(cè)能力進(jìn)行對(duì)比,引入CFD 結(jié)果與ROM 結(jié)果的方差(variance accounted for,VAF),其定義式如下

若兩種信號(hào)相同,則VAF 為100%;若兩種信號(hào)不同,則VAF 會(huì)變低.線(xiàn)性與非線(xiàn)性方法的VAF 值如表2所示.由表可知,非線(xiàn)性方法能夠提高辨識(shí)精度.

圖5 正弦激勵(lì)下直接CFD 模型輸出和非線(xiàn)性ROM 輸出對(duì)比Fig.5 Direct CFD model outputs vs nonlinear ROM outputs under the sinusoidal excitation

表2 線(xiàn)性與非線(xiàn)性方法的VAF 值Table 2 VAF values for linear and nonlinear methods

為了說(shuō)明非線(xiàn)性降階模型模型的計(jì)算效率,本文給出了白噪聲作為訓(xùn)練信號(hào)時(shí),對(duì)非線(xiàn)性ROM 和CFD 方法計(jì)算氣動(dòng)彈性響應(yīng)所需的時(shí)間進(jìn)行了比較,馬赫數(shù)為0.71 時(shí),采用CFD、線(xiàn)性ROM 和非線(xiàn)性ROM 方法的氣動(dòng)彈性響應(yīng)時(shí)間分別為187 h,39 s 和373 s.說(shuō)明氣動(dòng)彈性降階模型可在較為精確預(yù)測(cè)氣動(dòng)彈性響應(yīng)的前提下,大幅度縮短計(jì)算所需時(shí)間.由于Mannarino 和Dowell 提出的方法中,同樣分為線(xiàn)性狀態(tài)空間模型的辨識(shí)和非線(xiàn)性部分的優(yōu)化兩步,計(jì)算時(shí)間是兩者的疊加,因此計(jì)算效率與文中的ROM 方法是相當(dāng)?shù)?

圖6 krf=0.08 時(shí)諧波激勵(lì)下NACA 0012 翼型的ROM 驗(yàn)證Fig.6 The ROM validation of an NACA 0012 airfoil under harmonic excitation with krf=0.08

2.2 NACA0012 翼型非定常氣動(dòng)力驗(yàn)證

為了驗(yàn)證ROM 預(yù)測(cè)非定常氣動(dòng)力的精度,本文以NACA0012 翼型的非定常氣動(dòng)力計(jì)算為例,分別計(jì)算了折合頻率krf=0.08 和krf=0.2 的兩種諧波激勵(lì)下的升力系數(shù)和力矩系數(shù),并與直接CFD 方法和Mannarino 等提出的方法進(jìn)行對(duì)比.折合頻率krf的定義式為

圖7 krf=0.2 時(shí)諧波激勵(lì)下NACA 0012 翼型的ROM 驗(yàn)證Fig.7 The ROM validation of an NACA 0012 airfoil under harmonic excitation with krf=0.2

式中ω,c,V∞分別是俯仰角頻率、翼型弦長(zhǎng)和來(lái)流速度.如圖6 和圖7 所示,由非線(xiàn)性ROM 計(jì)算得到的升力系數(shù)Cl和力矩系數(shù)Cm與Mannarino 等[19]提出的方法相比,本文提出的方法與CFD 結(jié)果的吻合度更高,尤其是在力矩系數(shù)Cm的預(yù)測(cè)方面.

2.3 BACT 機(jī)翼跨音速顫振分析

為了研究非線(xiàn)性ROM 對(duì)BACT 機(jī)翼顫振邊界的預(yù)測(cè)精度,在不同馬赫數(shù)下(M∞=0.63,0.71,0.75,0.77,0.80,0.82,0.825)分別進(jìn)行了基于非線(xiàn)性狀態(tài)空間方法的系統(tǒng)辨識(shí),建立了不同的降階模型,并計(jì)算了顫振邊界.顫振計(jì)算方法為時(shí)域計(jì)算,求解過(guò)程框圖如圖8 所示,首先給定模態(tài)速度的初始擾動(dòng)代入到非線(xiàn)性ROM 方程(10)中求出氣動(dòng)力輸出,然后將來(lái)流動(dòng)壓q∞和氣動(dòng)力輸出代入到氣動(dòng)彈性方程(15),用預(yù)估?校正法進(jìn)行流固耦合計(jì)算,得到結(jié)構(gòu)位移ξ1,ξ2,再將其代入非線(xiàn)性ROM 方程中進(jìn)行迭代計(jì)算,最終得到收斂且穩(wěn)定的結(jié)構(gòu)位移矩陣.其中預(yù)估?校正法分為預(yù)估步和校正步,首先由預(yù)估步計(jì)算出模態(tài)解將預(yù)估步的解作用于結(jié)構(gòu)表面,在時(shí)間步為n+1 時(shí),流場(chǎng)收斂于該解,能夠計(jì)算出新的廣義力然后由校正步計(jì)算出xn+1.為了驗(yàn)證提出的ROM 法的有效性,將非線(xiàn)性ROM 的顫振結(jié)果與CFD 結(jié)果以及線(xiàn)性ROM結(jié)果進(jìn)行了比較,其中線(xiàn)性ROM 由相同的空氣動(dòng)力學(xué)系統(tǒng)訓(xùn)練數(shù)據(jù)的輸入和輸出來(lái)辨識(shí).圖9 表示通過(guò)直接CFD 模型、線(xiàn)性狀態(tài)空間辨識(shí)和非線(xiàn)性狀態(tài)空間辨識(shí)計(jì)算出的BACT 機(jī)翼顫振動(dòng)壓和顫振頻率的對(duì)比.如圖所示,通過(guò)非線(xiàn)性狀態(tài)空間辨識(shí)計(jì)算得到的顫振邊界與直接CFD 法計(jì)算得到的結(jié)果吻合較好.而與非線(xiàn)性狀態(tài)空間辨識(shí)相比,由線(xiàn)性狀態(tài)空間辨識(shí)確定的顫振邊界精度略低.從圖中還可以得到顫振動(dòng)壓和顫振頻率的相對(duì)誤差,相對(duì)誤差e定義式為[15,29]

圖8 顫振計(jì)算方法框圖Fig.8 Flutter calculation method block diagram

式中qCFD,qROM分別是通過(guò)直接CFD 法和線(xiàn)性/非線(xiàn)性狀態(tài)空間辨識(shí)計(jì)算出的顫振動(dòng)壓.圖9 表明了通過(guò)非線(xiàn)性狀態(tài)空間辨識(shí)計(jì)算出的顫振動(dòng)壓和顫振頻率的相對(duì)誤差始終低于2%,線(xiàn)性狀態(tài)空間辨識(shí)計(jì)算結(jié)果的相對(duì)誤差始終低于3%.該結(jié)果說(shuō)明非線(xiàn)性狀態(tài)空間辨識(shí)比線(xiàn)性狀態(tài)空間辨識(shí)更能準(zhǔn)確地表征BACT 模型的氣動(dòng)彈性行為.

圖9 BACT 機(jī)翼顫振邊界對(duì)比Fig.9 Comparison of the flutter boundaries of the BACT wing

2.4 極限環(huán)振蕩分析

極限環(huán)振蕩是另一種非常重要的非線(xiàn)性氣動(dòng)彈性現(xiàn)象[30-33],其非線(xiàn)性來(lái)源于氣動(dòng)力的非線(xiàn)性.由圖9 可知,當(dāng)馬赫數(shù)為0.71 時(shí),BACT 機(jī)翼的顫振動(dòng)壓為160 psf.因此,當(dāng)飛行動(dòng)壓高于顫振臨界動(dòng)壓(即分叉點(diǎn))時(shí),系統(tǒng)會(huì)發(fā)生極限環(huán)振蕩.當(dāng)前的非線(xiàn)性ROM 具有預(yù)測(cè)這一非線(xiàn)性氣動(dòng)彈性現(xiàn)象的能力,可將其應(yīng)用到BACT 機(jī)翼極限環(huán)振蕩現(xiàn)象的研究中.圖10 給出了在M∞=0.825,q∞=190 psf 下結(jié)構(gòu)兩階模態(tài)的LCO 響應(yīng).使用的模態(tài)速度的初始擾動(dòng)=0.001.由圖可見(jiàn),受空氣動(dòng)力非線(xiàn)性的影響,當(dāng)來(lái)流動(dòng)壓為190 psf 時(shí),BACT 機(jī)翼發(fā)生極限環(huán)振蕩.

圖10 非線(xiàn)性ROM 法的LCO 響應(yīng)Fig.10 LCO response via the nonlinear ROM method

圖11 極限環(huán)振蕩的幅值和頻率比較Fig.11 Comparison of the amplitude and frequency of the limit-cycle oscillation

為進(jìn)一步驗(yàn)證該非線(xiàn)性ROM 在預(yù)測(cè)BACT 機(jī)翼LCO 現(xiàn)象時(shí)的精度,將非線(xiàn)性ROM 計(jì)算出的LCO相圖和頻譜同直接CFD 方法計(jì)算出的結(jié)果進(jìn)行了對(duì)比.由圖11 可知,由非線(xiàn)性ROM 與CFD 方法計(jì)算出的兩階結(jié)構(gòu)模態(tài)LCO 頻率的相對(duì)誤差幾乎為零.而由兩階結(jié)構(gòu)模態(tài)相圖可以看出,與CFD 結(jié)果相比,兩階結(jié)構(gòu)模態(tài)的LCO 幅值相對(duì)誤差分別為19.0%和25.3%.極限環(huán)振蕩的幅值相對(duì)誤差較大的原因在于,在非線(xiàn)性狀態(tài)空間辨識(shí)的過(guò)程中,并未刻意設(shè)計(jì)訓(xùn)練信號(hào)的幅值和帶寬以提高其在預(yù)測(cè)LCO 現(xiàn)象時(shí)的精度.

3 結(jié)論

圖11 極限環(huán)振蕩的幅值和頻率比較(續(xù))Fig.11 Comparison of the amplitude and frequency of the limit-cycle oscillation(continued)

本文提出了一種基于非線(xiàn)性狀態(tài)空間辨識(shí)的氣動(dòng)彈性模型降階方法,實(shí)現(xiàn)跨音速氣動(dòng)彈性行為的高效預(yù)測(cè).為了驗(yàn)證該降階模型在預(yù)測(cè)跨音速氣動(dòng)彈性力學(xué)行為的精確性,本文以三維機(jī)翼為研究對(duì)象,分別從基于非線(xiàn)性降階模型的氣動(dòng)力辨識(shí)、跨聲速顫振邊界計(jì)算和極限環(huán)振蕩預(yù)測(cè)三方面進(jìn)行了算例驗(yàn)證.對(duì)比ROM 和CFD 在M∞=0.63,0.71,0.75,0.77,0.80,0.82,0.825 時(shí)的顫振邊界,并比較了非線(xiàn)性ROM 和直接CFD 計(jì)算出的LCO 響應(yīng).結(jié)果表明,該ROM 能夠與CFD 方法得到的空氣動(dòng)力學(xué)響應(yīng)相吻合,通過(guò)非線(xiàn)性ROM 計(jì)算出顫振動(dòng)壓和顫振頻率的相對(duì)誤差始終低于2%.在誤差允許范圍內(nèi),該ROM能夠得到與傳統(tǒng)CFD 方法計(jì)算氣動(dòng)彈性問(wèn)題相一致的結(jié)果,從而降低計(jì)算成本.在后續(xù)的工作中,需要進(jìn)一步提高極限環(huán)振蕩幅值的精度,計(jì)算在不同來(lái)流動(dòng)壓條件下極限環(huán)的變化趨勢(shì),并考慮變馬赫數(shù)的影響,研究具有馬赫數(shù)自適應(yīng)的非線(xiàn)性降階模型.