強(qiáng)迫振動(dòng)和極限環(huán)振動(dòng)時(shí)域響應(yīng)的高效預(yù)測(cè)方法研究

劉　艷 ， 白俊強(qiáng) ， 華　俊 ， 劉　南

(1.西北工業(yè)大學(xué) 航空學(xué)院，西安　710072；2. 中國(guó)航空研究院，北京　100012)

強(qiáng)迫振動(dòng)和極限環(huán)振動(dòng)時(shí)域響應(yīng)的高效預(yù)測(cè)方法研究

劉艷1， 白俊強(qiáng)1， 華俊2， 劉南1

(1.西北工業(yè)大學(xué) 航空學(xué)院，西安710072；2. 中國(guó)航空研究院，北京100012)

建立一種基于改進(jìn)Kriging的KSBRF(Kriging-Surrogate-Based Recurrence Framework)降階模型(Reduced-Order Model, ROM)，用于高效地預(yù)測(cè)非線性非定常氣動(dòng)力及力矩、極限環(huán)振動(dòng)(Limit Cycle Oscillations, LCO)等。首先基于Kriging代理模型建立非線性系統(tǒng)輸入-輸出關(guān)系的循環(huán)預(yù)測(cè)框架。然后對(duì)翼型做沉浮/俯仰組合運(yùn)動(dòng)時(shí)的非線性非定常氣動(dòng)力進(jìn)行辨識(shí)。結(jié)果表明：在固定來(lái)流馬赫數(shù)情況下，KSBRF預(yù)測(cè)結(jié)果與CFD計(jì)算結(jié)果吻合良好，阻力系數(shù)及俯仰力矩系數(shù)等的平均預(yù)測(cè)誤差均在2.0%以內(nèi)；在變來(lái)流馬赫數(shù)情況下，升力系數(shù)和俯仰力矩系數(shù)的平均預(yù)測(cè)誤差在2.5%以內(nèi)，而阻力系數(shù)的預(yù)測(cè)誤差則稍大，但也不超過(guò)7.0%。通過(guò)研究不同m,n取值對(duì)模型精度的影響，得出考慮歷史效應(yīng)有利于提高當(dāng)前時(shí)間步的模型預(yù)測(cè)精度。除此之外，還對(duì)NACA64A010翼型的LCO振動(dòng)形態(tài)進(jìn)行預(yù)測(cè)，預(yù)測(cè)結(jié)果與CFD計(jì)算結(jié)果吻合良好，誤差均保證均小于5.17%，節(jié)省近三分之二的計(jì)算時(shí)長(zhǎng)。

代理模型；降階模型；Kriging插值；非定常氣動(dòng)力；極限環(huán)振動(dòng)

目前，一些非線性流動(dòng)現(xiàn)象分析往往需要高精度的計(jì)算流體力學(xué)(Computational Fluid Dynamics，CFD)方法。例如在跨音速、接近失速等狀態(tài)下，非定常氣動(dòng)力響應(yīng)分析、顫振分析、極限環(huán)振動(dòng)(Limit Cycle Oscillations，LCO)分析等分析過(guò)程復(fù)雜，應(yīng)用高精度的CFD分析方法的計(jì)算代價(jià)太過(guò)巨大。在實(shí)際應(yīng)用中通常會(huì)采用一些基于線性假設(shè)條件的CFD方法。這些CFD方法在計(jì)算時(shí)間上是可以接受的，但是對(duì)于非線性較強(qiáng)的流動(dòng)現(xiàn)象，計(jì)算精度往往無(wú)法達(dá)到應(yīng)用要求。

為了提高氣動(dòng)力計(jì)算效率，而不損失氣動(dòng)力的非線性特性，一些氣動(dòng)力降階模型(Reduced-Order Modelling, ROM)方法，得到了廣泛應(yīng)用，大大節(jié)省了計(jì)算時(shí)間。典型的降階模型有，Proper Orthogonal Decomposition (POD)[1-5]，Volterra級(jí)數(shù)[4,6]，高階諧波平衡[7]，基于代理模型法[8-12]等。Volterra級(jí)數(shù)[4,6]適用于弱非線性問(wèn)題，對(duì)于跨音速問(wèn)題則需要辨識(shí)出高階Volterra核，計(jì)算耗費(fèi)非常大。POD降階方法[1-5]最近發(fā)展較快，但是在可壓流和黏性問(wèn)題的發(fā)展較為緩慢，高階非線性項(xiàng)的Galerkin投影處理難度很高。高階諧波平衡[7]適用于周期性運(yùn)動(dòng)問(wèn)題分析，但是需要對(duì)計(jì)算程序進(jìn)行較大的改動(dòng)，復(fù)雜度很高。而基于代理模型的降階方法[8-12]則利用代理模型辨識(shí)出CFD的輸入輸出關(guān)系，從而實(shí)現(xiàn)非定常氣動(dòng)力的快速預(yù)測(cè)。其思路與Volterra級(jí)數(shù)類似，區(qū)別在于Volterra級(jí)數(shù)利用脈沖激勵(lì)-響應(yīng)之間的關(guān)系，而代理模型則利用大量CFD樣本中輸入-輸出之間的關(guān)系。相比而言，代理模型具有更好的適用性，及非線性預(yù)測(cè)能力。Liea等[1,13]應(yīng)用無(wú)黏POD降階模型對(duì)F-16飛機(jī)進(jìn)行了時(shí)域彈性響應(yīng)分析，預(yù)測(cè)變來(lái)流馬赫數(shù)及變迎角條件下的時(shí)域氣動(dòng)彈性響應(yīng)分析。王博斌等[14]采用輸入反饋ANN網(wǎng)絡(luò)建立了NACA64A010翼型的非定常氣動(dòng)力模型。Glaz等[15-18]建立了SBRF(Srrogate-Based Recurrence Framework)方法，并應(yīng)用該方法對(duì)翼型沉浮、俯仰、沉浮俯仰等條件下進(jìn)行了時(shí)域法非定常氣動(dòng)力的預(yù)測(cè)，并進(jìn)一步將這個(gè)方法應(yīng)用在旋翼機(jī)設(shè)計(jì)及其氣動(dòng)彈性分析中。

如今代理模型方法發(fā)展眾多分支，如應(yīng)用廣泛的Kriging[19]、BP神經(jīng)網(wǎng)絡(luò)[20]、RBF神經(jīng)網(wǎng)絡(luò)[21]、支持向量機(jī)[22]等，一些學(xué)者進(jìn)行了綜述和比較[23-24]。其中，Kriging方法不僅可以得到具有高階非線性和無(wú)偏特性的預(yù)測(cè)結(jié)果，同時(shí)可獲得預(yù)測(cè)值的方差量，因此被認(rèn)為是一種精準(zhǔn)可靠和具有吸引力的代理模型方法。

本文主要基于改進(jìn)的Kriging插值建立KSBRF(Kriging-Srrogate-Based Recurrence Framework)降階模型，研究KSBRF降階模型預(yù)測(cè)翼型在固定來(lái)流馬赫數(shù)和變來(lái)流馬赫數(shù)兩種情況下，做俯仰/沉浮組合運(yùn)動(dòng)時(shí)的氣動(dòng)力響應(yīng)歷程，證明KSBRF降階模型預(yù)測(cè)非線性非定常氣動(dòng)力響應(yīng)歷程的能力；除此之外，還應(yīng)用KSBRF降階模型預(yù)測(cè)了NACA64A010翼型的LCO振動(dòng)特性。

本文的具體工作內(nèi)容如下：

(1) 基于改進(jìn)的Kriging插值發(fā)展一種氣動(dòng)力KSBRF降階模型，對(duì)二維翼型做俯仰/沉浮組合運(yùn)動(dòng)的非定常氣動(dòng)力進(jìn)行預(yù)測(cè)。證明KSBRF降階模型預(yù)測(cè)非定常氣動(dòng)力(升力、阻力、力矩)的能力。分析不同參數(shù)對(duì)KSBRF降階模型精度的影響。

(2) 建立KSBRF降階模型，通過(guò)改變輸入變量中參數(shù)，研究對(duì)建立KSBRF降階模型的精度影響。

(3) 以NACA64A010翼型為例，應(yīng)用KSBRF降階模型預(yù)測(cè)其LCO的振動(dòng)特性，證明KSBRF降階模型的預(yù)測(cè)能力。

1　理論基礎(chǔ)

1.1輸入輸出關(guān)系的動(dòng)力學(xué)表示

代理模型可以將時(shí)域內(nèi)變量之間的輸入輸出關(guān)系轉(zhuǎn)化成一個(gè)動(dòng)力系統(tǒng)狀態(tài)空間方程的形式，例如非定常氣動(dòng)力、振動(dòng)形態(tài)等。這個(gè)動(dòng)力系統(tǒng)不僅是時(shí)間的函數(shù)，同時(shí)也是之前時(shí)刻氣動(dòng)力輸出結(jié)果的函數(shù)。通常非線性時(shí)間離散動(dòng)力系統(tǒng)[26-27]表示如下：

q(t)=Gq(q(t-Δt),u(t-Δt))

y(t)=Gy(q(t))

(1)

式中：向量q為系統(tǒng)狀態(tài)，向量u為系統(tǒng)的外部輸入，Gq為前一個(gè)時(shí)間步的瞬時(shí)狀態(tài)函數(shù)，y(t)為當(dāng)前時(shí)刻t的輸出，Gy為系統(tǒng)狀態(tài)到輸出的映射關(guān)系。

以非定常氣動(dòng)力模型為例，方程(1)表示的非線性系統(tǒng)方程代表離散的Navier-Stokes(NS)方程。向量q可以是包含NS方程計(jì)算域內(nèi)每個(gè)網(wǎng)格點(diǎn)上的流場(chǎng)參數(shù)(例如，密度、速度、能量等)其維度是網(wǎng)格量的倍數(shù)，與全階模型的階數(shù)一致。而系統(tǒng)的輸出為：

y(t)=Cl(t),Cm(t),orCd(t)

(2)

式中：Cl(t)，Cm(t)，Cd(t)分別為非定常升力系數(shù)、俯仰力矩系數(shù)和阻力系數(shù)。根據(jù)線性氣動(dòng)力理論，非定常升力與力矩是俯仰角、沉浮位移及其關(guān)于時(shí)間導(dǎo)數(shù)的函數(shù)，因此，相關(guān)的外部輸入為：

(3)

因?yàn)閼?yīng)用全階計(jì)算系統(tǒng)完成方程(1)的計(jì)算是非常耗時(shí)的，所以建立一種基于代理模型的降階模型方法可以高效地計(jì)算出系統(tǒng)輸出變量y(t)。在不進(jìn)行方程(1)所描述的非線性方程求解的情況下，給定一個(gè)輸入變量，通過(guò)代理模型建立的輸入輸出關(guān)系就可以快速地得到結(jié)果y(t)。在文獻(xiàn)[27-28]中提到，如果向量q用有限次狀態(tài)的疊加表示，那么方程(1)等價(jià)于方程(4)：

y(t)=φ(u(t),u(t-Δt),…,u(t-mΔt),

y(t-Δt),…,y(t-nΔt))

(4)

式中：函數(shù)φ代表輸入與輸出的映射函數(shù)，m,n分別代表需要考慮的時(shí)間歷史影響的數(shù)量，m,n的選擇在第3節(jié)算例驗(yàn)證中進(jìn)行詳細(xì)討論。

1.2非線性映射函數(shù)的近似表示

(1) 應(yīng)用試驗(yàn)設(shè)計(jì)方法選取有限個(gè)訓(xùn)練樣本進(jìn)行CFD計(jì)算。每個(gè)訓(xùn)練樣本包括做俯仰/沉浮振動(dòng)：

(5)

(6)

方程(5)、(6)給定的振動(dòng)頻率為

(7)

式中：V0來(lái)流速度，b為參考長(zhǎng)度。

采用優(yōu)化的拉丁超立方方法[27](Optimal Latin Hypercube, OLH)進(jìn)行樣本采樣，保證樣本點(diǎn)更加均勻地填充整個(gè)樣本空間。

(2) 通過(guò)CFD分析獲取每個(gè)樣本點(diǎn)非定常氣動(dòng)力的時(shí)間響應(yīng)歷程，從結(jié)果中提取方程式(4)所需要的時(shí)域響應(yīng)結(jié)果。

(8)

方程(4)中的對(duì)應(yīng)樣本輸出變量y的輸入變量x表示如下：

x(t)=[uj(t)uj(t-Δt)…uj(t-mΔt)

yj(t-Δt)…yj(t-nΔt)]j=1,2,…,NT

(9)

式中：yj，uj分別由方程(2)、(3)給出第j個(gè)樣本點(diǎn)的輸出、外部輸入向量，向量X的維度為k*(m+1)+n，其中k是由方程(3)給定的外部輸入函數(shù)確定的。

為了更有效的預(yù)測(cè)非線性氣動(dòng)力，要求在構(gòu)建代理映射函數(shù)時(shí)使用的插值方法具有預(yù)測(cè)非線性函數(shù)的能力。因此，本文應(yīng)用滿足以上要求的改進(jìn)Kriging插值方法構(gòu)建代理映射函數(shù)。

(10)

式中：f(x)=[f1(x),f2(x),…,fp(x]T為回歸多項(xiàng)式，β=[β1,β2,…,βp]T為回歸參數(shù)。

假設(shè)隨機(jī)項(xiàng)z(x)的均值為0，z(w)和z(x)之間的相關(guān)度如下：

E[z(w)z(x)]=σ2R(θ,w,x)

(11)

式中：σ2是隨機(jī)項(xiàng)的方差，R(θ,w,x)為點(diǎn)w和x之間的相關(guān)函數(shù)。根據(jù)已知樣本點(diǎn)，將式(10)寫(xiě)成矩陣形式如下：

(12)

式中：

(13)

F(xi)=[f(x1),f(x2),…f(xm)]T

(14)

B=[β]

(15)

根據(jù)無(wú)偏估計(jì)和最大似然估計(jì)可得：

β=(FTRF)-1FTR-1Y

(16)

可以建立代理模型：

(17)

式中：y是已知樣本點(diǎn)的函數(shù)值。R是已知樣本點(diǎn)處的相關(guān)矩陣，r是未知點(diǎn)和已知樣本點(diǎn)之間的相關(guān)向量，在工程應(yīng)用中最常采用的相關(guān)函數(shù)為Gauss函數(shù)，如下所示：

(18)

式中：θ=[θ1,θ2,…,θn]T是空間相關(guān)參數(shù)矢量，可以通過(guò)調(diào)整θ的大小從而改變預(yù)測(cè)精度，本文采用粒子群優(yōu)化算法，以Kriging模型的誤差為適應(yīng)值，對(duì)θ進(jìn)行優(yōu)化，從而提高代理模型的精度，但是優(yōu)化過(guò)程中涉及到大量的矩陣運(yùn)算，增加了構(gòu)建代理模型時(shí)長(zhǎng)。

1.3應(yīng)用KSBRF構(gòu)建時(shí)域降階模型的流程

(19)

(20)

圖(1)為KSBRF降階模型的框架圖，同時(shí)也描述出KSBRF降階模型過(guò)程需要預(yù)先進(jìn)行一些時(shí)間步長(zhǎng)為Δt的預(yù)測(cè)。值得注意的是，提前預(yù)測(cè)n個(gè)時(shí)間步的響應(yīng)是KSBRF初始化所必需的。

圖1　KSBRF降階模型框架Fig.1 Surrogate-based recurrence framework

2　氣動(dòng)彈性控制方程的時(shí)域解法

應(yīng)用Lagrange方程，飛行器的結(jié)構(gòu)運(yùn)動(dòng)方程可表示成矩陣形式[31]，

(21)

式中，M是廣義質(zhì)量矩陣，D是廣義阻尼矩陣，K是廣義剛度矩陣，F(xiàn)是廣義氣動(dòng)力。

為了便于利用Runge-Kutta法求解結(jié)構(gòu)運(yùn)動(dòng)方程，引入狀態(tài)變量X，

(22)

然后將結(jié)構(gòu)運(yùn)動(dòng)方程(20)變?yōu)椋?/p>

(23)

對(duì)于兩自由度二元機(jī)翼的氣彈模型，一個(gè)自由度是剛心的垂直位移h，向下為正；另一個(gè)是機(jī)翼繞剛心的轉(zhuǎn)角α，迎風(fēng)抬頭為正。Cl、Cm分別為升力系數(shù)和俯仰力矩系數(shù)。rα為繞彈性軸的旋轉(zhuǎn)半徑。xα為重心與鉸鏈軸之間的無(wú)量綱距離。ωh、ωα分別為沉浮模態(tài)、俯仰模態(tài)的自然頻率。b、μ、V、V*分別為半弦長(zhǎng)、質(zhì)量比、自由來(lái)流速度和減縮速度。對(duì)于兩自由度二元機(jī)翼的氣彈模型，無(wú)量綱時(shí)間為τ=ωαt，則可通過(guò)以下公式計(jì)算得到需要的矩陣及向量，

以NACA64A010翼型[32]為例，相關(guān)參數(shù)：Ma=0.8，α=0.0 deg，xα=0.25，rα2=0.75，ωh/ωα=0.5，μ=75，V*=0.67。研究了初始擾動(dòng)對(duì)LCO擾動(dòng)形態(tài)的影響。

由圖2可得，初始擾動(dòng)對(duì)LCO的收斂速度具有很重要的影響，但不同的初始擾動(dòng)均能達(dá)到相同的收斂狀態(tài)。給定大初始擾動(dòng)較給定小初始擾動(dòng)分析LCO的收斂速度要快。在對(duì)未知問(wèn)題進(jìn)行分析時(shí)，如果能夠估計(jì)極限環(huán)的振動(dòng)幅值，則可以使LCO的計(jì)算更快地收斂。

圖2　不同初始擾動(dòng)下LCOs響應(yīng)對(duì)比(V*=0.67)Fig.2 Comparison of LCO responses between different initial disturbances

3　算例驗(yàn)證

應(yīng)用式(24)[16]的誤差矩陣判斷方法評(píng)估KSBRF時(shí)域降階模型與CFD分析結(jié)果的誤差，式中代表預(yù)測(cè)響應(yīng)區(qū)間的平均誤差值。

(24)

3.1非定常氣動(dòng)力預(yù)測(cè)

以NACA0012翼型為例，預(yù)測(cè)其做俯仰/沉浮組合運(yùn)動(dòng)時(shí)非定常氣動(dòng)力隨時(shí)間的響應(yīng)歷程。考慮確定馬赫數(shù)和變馬赫數(shù)兩種情況下，證明KSBRF降階模型的預(yù)測(cè)能力。需要注意的是，應(yīng)在已經(jīng)收斂的定常流場(chǎng)解得基礎(chǔ)上進(jìn)行非定常氣動(dòng)力分析。

3.1.1固定來(lái)流馬赫數(shù)翼型俯仰/沉浮組合運(yùn)動(dòng)

本算例中采用NACA0012翼型進(jìn)行俯仰/沉浮組合運(yùn)動(dòng)時(shí)分析KSBRF降階模型對(duì)非定常氣動(dòng)力及力矩隨時(shí)間響應(yīng)的預(yù)測(cè)能力。應(yīng)用KSBRF降階模型預(yù)測(cè)500個(gè)時(shí)間步內(nèi)的非定常氣動(dòng)力，圖3中給出500個(gè)時(shí)間步的預(yù)測(cè)結(jié)果。初始計(jì)算狀態(tài)為：Ma=0.6，Re=4 800 000，初始攻角α0=0°。

(25)

0.01≤kθ≤0.3

(26)

(27)

0.01≤kh≤0.3

(28)

(29)

分別選取三組方程(9)中m，n構(gòu)建KSBRF降階模型，得到的對(duì)測(cè)試樣本點(diǎn)的預(yù)測(cè)誤差見(jiàn)表1。圖3對(duì)比了非定常氣動(dòng)力隨時(shí)間響應(yīng)的歷程的非線性ROM預(yù)測(cè)值和CFD分析值。給出三組不同的m、n取值，通過(guò)表1和圖3的結(jié)果，KSBRF模型對(duì)20個(gè)測(cè)試樣本的升力系數(shù)Cl、阻力系數(shù)Cd及俯仰力矩系數(shù)Cmy的預(yù)測(cè)結(jié)果，可以得出，在n值相同的情況下，不同的m值對(duì)結(jié)果的影響明顯，尤其是對(duì)較難準(zhǔn)確預(yù)測(cè)的阻力系數(shù)和俯仰力矩系數(shù)，m=2時(shí)KSBRF模型的預(yù)測(cè)精度明顯高于m=0。當(dāng)m=0、n=2時(shí)，對(duì)氣動(dòng)力及力矩幅值處的預(yù)測(cè)值較CFD計(jì)算值偏大，而且具有鋸齒形。m=2、n=2或n=3時(shí)，預(yù)測(cè)結(jié)果與CFD計(jì)算結(jié)果均吻合良好。但n=3時(shí)的預(yù)測(cè)結(jié)果稍優(yōu)于n=3。說(shuō)明考慮歷史結(jié)果對(duì)模型精度的提高有益。

除此之外，將KSBRF模型預(yù)測(cè)的結(jié)果與基于RBF代理模型的SBRF-RBF模型預(yù)測(cè)結(jié)果進(jìn)行對(duì)比，結(jié)果如圖4所示。由圖可得，兩種方法對(duì)升力系數(shù)的預(yù)測(cè)結(jié)果均與CFD結(jié)算吻合較好，而對(duì)阻力系數(shù)和俯仰力矩系數(shù)的預(yù)測(cè)KSBRF的結(jié)果較好，誤差分別為0.76%、0.06%，而SBRF-RBF的預(yù)測(cè)誤差分別為0.83%、0.08%。采用SBRF-RBF方法預(yù)測(cè)10個(gè)測(cè)試樣本耗時(shí)307 s，而采用KSBRF方法需要845 s。

值得注意的是，這些分析只是針對(duì)本算例開(kāi)展的定量分析，針對(duì)不同問(wèn)題還應(yīng)從實(shí)際情況出發(fā)。

圖3　非定常Cl、Cd和CmyFig.3 Unsteady lift, moment, and drag coefficients

圖4　KSBRF結(jié)果與SBRF-RBF結(jié)果對(duì)比Fig.4 Comparison of results between KSBRF and

模型Cl誤差Cd誤差Cmy誤差m=0n=20.025-4.28%(均值0.96%)0.07-7.40%(均值1.60%)0.002-12.94%(均值0.74%)m=2n=20.02-3.81%(均值0.49%)0.11-5.80%(均值1.79%)0.001-2.12%(均值0.14%)m=2n=30.008-2.28%(均值0.247%)0.096-5.66%(均值1.56%)0.0008-1.69%(均值0.103%)

3.1.2變馬赫數(shù)下翼型俯仰/沉浮組合運(yùn)動(dòng)

本節(jié)將測(cè)試KSBRF時(shí)域降階模型在考慮變來(lái)流馬赫數(shù)情況下，翼型做俯仰/沉浮組合運(yùn)動(dòng)時(shí)對(duì)氣動(dòng)力及俯仰力矩系數(shù)的預(yù)測(cè)能力。初始計(jì)算狀態(tài)及沉浮/俯仰組合振動(dòng)形式同3.1.A節(jié)，考慮馬赫數(shù)的變化區(qū)間為：

0.6≤M≤1.1

(30)

應(yīng)用OLH取樣方法生成200樣本，隨機(jī)選取180個(gè)樣本點(diǎn)為訓(xùn)練樣本，其余20個(gè)作為測(cè)試樣本。

應(yīng)用KSBRF時(shí)域降階方法(m=2，n=2)預(yù)測(cè)20個(gè)樣本的預(yù)測(cè)誤差見(jiàn)表2。本算例中，隨著馬赫數(shù)的變化，亞音速流場(chǎng)是近似線性的，而跨音速流動(dòng)是具有強(qiáng)烈的非線性的特征。由表2的誤差結(jié)果可以看出，升力系數(shù)、阻力系數(shù)及俯仰力矩系數(shù)的平均誤差分別為2.482%、6.964%和0.832%。證明建立的KSBRF降階模型對(duì)不同馬赫數(shù)條件下的氣動(dòng)力及力矩具有較高的預(yù)測(cè)精度。

圖5中給出其中2個(gè)測(cè)試樣本全階CFD分析結(jié)果與KSBRF預(yù)測(cè)結(jié)果的對(duì)比曲線。圖(a)樣本的升力系數(shù)、阻力系數(shù)和俯仰力矩系數(shù)的預(yù)測(cè)誤差分別為0.283%、0.995%、0.054%。圖(b)樣本的升力系數(shù)、阻力系數(shù)和俯仰力矩系數(shù)的預(yù)測(cè)誤差分別為0.616%、2.228%、0.165%。

(a) θ=5.15,kθ=0.0278,h=0.008,kh=0.0327,M=1.049(b) θ=6.15,kθ=0.0062,h=0.0158,kh=0.0164,M=0.878圖5 變來(lái)流馬赫數(shù)條件下非定常Cl、Cd和CmyFig.5Unsteadylift,drag,andpitchmomentcoefficientsundervariableMachnumber

表2　考慮馬赫數(shù)不確定性的KSBRF降階模型對(duì)10個(gè)測(cè)試樣本的預(yù)測(cè)誤差

3.2LCO振動(dòng)形態(tài)預(yù)測(cè)

以上算例證明了KSBRF降階模型對(duì)非定常氣動(dòng)力的預(yù)測(cè)能力。在本算例中將進(jìn)一步研究KSBRF降階模型對(duì)跨音速區(qū)的LCO振動(dòng)形態(tài)的預(yù)測(cè)能力。以NACA64A010翼型為例，相關(guān)參數(shù)見(jiàn)第2節(jié)，研究在不同無(wú)量綱速度下LCO狀態(tài)。八個(gè)訓(xùn)練樣本點(diǎn)：V*=0.670，0.675，0.680，0.690，0.720，0.740，0.770，0.800，七個(gè)測(cè)試樣本見(jiàn)表4。

圖6給出兩個(gè)不同無(wú)量綱速度下LCO沉浮振動(dòng)、俯仰振動(dòng)、沉浮振動(dòng)速率和俯仰振動(dòng)速率的收斂曲線對(duì)比圖，左側(cè)圖為全體結(jié)果對(duì)比，右側(cè)圖后1 000步結(jié)果放大圖?？傻?，非線性ROM預(yù)測(cè)結(jié)果與CFD計(jì)算結(jié)果吻合得很好，僅在一些區(qū)域的振動(dòng)幅值處有小的差別，總體上對(duì)LCO振動(dòng)周期、振動(dòng)幅值預(yù)測(cè)得較準(zhǔn)。

(a1)V*=0.6725時(shí)全體結(jié)果圖(a2)V*=0.6725時(shí)局部結(jié)果圖

(b1)V*=0.700時(shí)全體結(jié)果圖(b2)V*=0.700時(shí)局部結(jié)果圖圖6 LCO響應(yīng)值對(duì)比Fig.6ComparisonofLCOsresponses

圖7為隨無(wú)量綱速度的LCO響應(yīng)曲線，同樣顯示出非線性ROM預(yù)測(cè)結(jié)果與CFD計(jì)算結(jié)果吻合很好，證明非線性ROM具有很高的預(yù)測(cè)精度。

表3為非線性ROM計(jì)算與直接CFD計(jì)算的耗時(shí)對(duì)比，可得非線性ROM方法比直接CFD計(jì)算要省時(shí)得多。計(jì)算8個(gè)訓(xùn)練樣本，預(yù)測(cè)14個(gè)狀態(tài)下的LCO響應(yīng)，比直接進(jìn)行CFD計(jì)算節(jié)省近2倍的時(shí)間。值得注意的是，模型一經(jīng)建立，只需要45 s就可以預(yù)測(cè)給定V*時(shí)的LCO。

表4為各測(cè)試樣本的預(yù)測(cè)誤差。除第7個(gè)測(cè)試樣本外，預(yù)測(cè)誤差均在3.5%之內(nèi)。與之前得到相同的結(jié)論，說(shuō)明此非線性ROM方法的精度很高。

圖7　CFD分析與ROM預(yù)測(cè)的LCO響應(yīng)結(jié)果對(duì)比Fig.7 Comparison of LCOs behaviors between CFD analysis and ROM identification

方法時(shí)間步長(zhǎng)時(shí)間步長(zhǎng)CFD計(jì)算22個(gè)計(jì)算狀態(tài),每個(gè)計(jì)算狀態(tài)耗時(shí)67minTCFD=22×67min=1474min=24h34minROM1.計(jì)算8訓(xùn)練樣本2.預(yù)測(cè)14狀態(tài)的LCO,每個(gè)狀態(tài)建模到預(yù)測(cè)結(jié)果耗時(shí)3sTROM=8×67min+14×45s=546min30s=9h6min42s

表4　各測(cè)試樣本預(yù)測(cè)誤差

4　結(jié)　論

經(jīng)過(guò)對(duì)基于KSBRF的非線性ROM與直接CFD計(jì)算的對(duì)比研究，得到如下結(jié)論：

(1) 應(yīng)用非線性ROM模型KSBRF對(duì)固定來(lái)流馬赫數(shù)下，翼型做沉浮/俯仰運(yùn)動(dòng)的非定常氣動(dòng)力進(jìn)行預(yù)測(cè)。不同m,n取值對(duì)模型精度有顯著影響，考慮歷史效應(yīng)有利于提高當(dāng)前時(shí)間步的模型預(yù)測(cè)精度。

(2) 應(yīng)用非線性ROM模型預(yù)測(cè)考慮變來(lái)流馬赫數(shù)情況，翼型做俯仰/沉浮組合運(yùn)動(dòng)時(shí)的非定常氣動(dòng)力，結(jié)果與CFD分析結(jié)果相比，最大平均誤差為1.79%。

(3) 非線性ROM能夠預(yù)測(cè)出LCO的振動(dòng)周期及振動(dòng)幅值，與直接CFD分析結(jié)果相比，升力系數(shù)和俯仰力矩系數(shù)的平均預(yù)測(cè)誤差在2.5%以內(nèi)，阻力系數(shù)不超過(guò)7.0%。

(4) 非線性ROM方法的計(jì)算效率較直接CFD分析提高顯著，只占直接CFD分析時(shí)長(zhǎng)的近三分之一，并且只需要45 s就可預(yù)測(cè)未知V*時(shí)的LCO。

[1] Lieu T, Farhat C. Adaptation of aeroelastic reduced-order models and application to an F-16 configuration[J]. AIAA Journal,2007,45 (6):1244-1257.

[2] Hall K C, Thomas J P, Dowell E H. Proper orthogonal decomposition technique for transonic unsteady aerodynamic flows[J]. AIAA Journal,2000 38 (10):1853-1862.

[3] Thomas J P, Dowell E H, Hall K C. Using automatic differentiation to create a nonlinear reduced-order-model aerodynamic solver[J]. AIAA Journal,1020,48(1):19-24.

[4] Lucia D J, Beran P S, Silva W A. Reduced-order modeling: new approaches for computational physics[J]. Progress in Aerospace Sciences,2004,40 (1/2):51-117.

[5] 邱亞松，白俊強(qiáng)，華俊. 基于本征正交分解和代理模型的流場(chǎng)預(yù)測(cè)方法[J].航空學(xué)報(bào),2013,34 (6):1249-1260.

QIU Ya-song, BAI Jun-qiang, HUA Jun. Flow field estimation method based on proper orthogonal decomposition and surrogate model[J]. Acta Aeronautica et Astronautica Sinica,2013,34 (6):1249-1260.

[6] Silva W. Identification of nonlinear aeroelastic systems based on the Volterra theory: progress and opportunities[J]. Nonlinear Dynamics,2005,39:25-62.

[7] Thomas J P, Custer C H, Dowell E H, et al. F-16 fighter aeroelastic computations using a harmonic balance implementation of the OVERFLOW2 flow solver[C]//AIAA Paper 2010-2632. Presented at the 51stAIAA/ASME/ASCE/AHS/ASC Structures, Structural Dynamics, and Meterials Conference.

[8] Liu L, Friedmann P P, Padthe A K. An approximate unsteady aerodynamic model for flapped airfoils including improved drag predictions[C]//Proceedings of the 34th European Rotorcraft Forum,Liverpool,UK,2008:1037-1081.

[9] Trizila P C, Kang C, Visbal M R, et al. Unsteady fluid physics and surrogate modeling of low Reynolds number, flapping airfoils[C]//38th Fluid Dynamics Conference and Exhibit, Seattle, WA, AIAA Paper 2008-3821,2008.

[10] Trizila P C, Kang C, Visbal M R, et al. A surrogate model approach in 2D versus 3D flapping wing aerodynamic analysis[C]//12th AIAA/ISSMO Multidisciplinary Analysis and Optimization Conference, Victoria,British Columbia, AIAA Paper 2008-5914,2008.

[11] Suresh S, Omkar S N, Mani V, et al. Lift coefficient prediction at high angle of attack using recurrent neural network[J]. Aerospace Science and Technology,2003,7(8):595-602.

[12] Marques F D, Anderson J. Identification and prediction of unsteady transonic aerodynamic loads by Multi-Layer Functionals[J]. Journal of Fluids and Structures,2001,15:83-106.

[13] Lieu T, Farhat C. Adaptation of POD-based aeroelastic ROMs for varying mach number and angle of attack: application to a complete F-16 configuration[C]//2005 U. S. Air Force T&E Days.2005-7666.

[14] Glaz B, Liu L, Friedmann P P. Reduced order nonlinear unsteady aerodynamic modeling using a Surrogate Based Recurrence Framework[J]. AIAA Journal, 2010, 48 (10):2418-2429.

[15] Glaz B, Liu L, Friedmann P P, et al. A Surrogate-Based approach to reduced-order dynamic stall modeling[J]. Journal of the American Helicopter Society, 2012, 57 (2): 1-9.

[16] Glaz B, Liu L, Friedmann P P, et al. A Surrogate Based Approach to reduced-order dynamic stall modeling[C]//51st AIAA/ASME/ASCHE/AHS/ASC Structures, Structural Dynamics and Materials Conference, Orlando, FL, AIAA Paper 2010-3042,2010,1-24.

[17] Glaz B, Friedmann P P, Liu L. Reduced order dynamic stall modeling with swept flow effects using a Surrogate Based Recurrence Framework[J]. AIAA Journal, 2013, 51 (4):910-921.

[18] 王博斌，張偉偉，葉正寅.基于神經(jīng)網(wǎng)絡(luò)模型的動(dòng)態(tài)非線氣動(dòng)力辨識(shí)方法[J]. 航空學(xué)報(bào),2010, 31(7):1379-1388.

WANG Bo-bin, ZHANG Wei-wei, YE Zheng-yin. Unsteady nonlinear aerodynamics identification based on neural network model[J]. Acta Aeronautica et Astronautica Sinica,2010,31(7):1379-1388.

[19] Sacks J, Schiller S B, Welch W J. Designs for computer experiments[J]. Technometrics,1989,31 (1):41-47.

[20] 閻平凡，張長(zhǎng)水. 人工神經(jīng)網(wǎng)絡(luò)與模擬進(jìn)化計(jì)算[M].北京：清華大學(xué)出版社,2005.

[21] Moody J E, Darken C J. Fast learning in networks of locally-tuned processing units[J]. Neural Computation,1989,1(2):281-294.

[22] Vapnik V N, Somla A. Support vector method for function approximation[J], Regression Estimation and Signal Processing,1997.

[23] Sacks J, Welch W J, Mitchell T J, et al. Design and analysis of computer experiments[J]. Statistical Science, 1989,4(4):409-423.

[24] Simpson T W, Toropov V, Balabanov V, et al. Design and analysis of computer experiments in multidisciplinary design optimization: a review of how we have come or not[C]//12th AIAA/ISSMO Multidisciplinary Analysis and Optimization Conference, Victoria, British Colombia, 10-12 September 2008.

[25] Queipo N V, Haftka R T, Shyy W, et al. Surrogate-Based analysis and optimization[J]. Progress in Aerospace Sciences, 2005, 41:1-28.

[26] Leontaritis I J, Billings S A, Input-output parametric models for nonlinear systems[J]. International Journal of Control,1985,41 (2): 03-344.

[27] Levin A U, Narendra K S. Control of nonlinear dynamical systems using neural networks—Part II: observability, identification, and control[J]. IEEE Transactions on Neural Networks,1996,7 (1):30-42.

[28] Jin R, Chen W, Sudjianto A. An efficient algorithm for constructing optimal design of computer experiments[J]. Journal of Statistical Planning and Inference,2005,134(1):268-287.

[29] Martin J, Simpson T. Use of Kriging models to approximate deterministic computer models[J]. AIAA Journal, 2005,43(4):853-863.

[30] Jones D R. A taxonomy of global optimization methods based on response surfaces[J]. Journal of Global Optimization,2001, 21:345-383.

[31] 趙永輝. 氣動(dòng)彈性力學(xué)與控制[M]. 北京：科學(xué)出版社. 2007.

[32] Thomas J P, Dowell E H, Hall K C. Nonlinear inviscid aerodynamic effects on transonic divergence, flutter, and limit-cycle oscillations[J]. AIAA Journal,2002,40(4):638-646.

Efficient prediction for time domain responses of forced oscillations and limit cycle oscillations

LIU Yan1, BAI Jun-qiang1, HUA Jun2, LIU Nan1

(1. School of Aeronautics, Northwestern Polytechnical University, Xi’an 710072, China；2. Chinese Aeronautical Research Institute, Beijing 100012, China)

A nonlinear reduced-order model (ROM) using Kriging surrogate-based recurrence framework (KSBRF) was built. Because Kriging interpolation has the ability of estimating nonlinear input-output relationship, nonlinear aerodynamic forces and characteristics of LCO (limit cycle oscillations) could be estimated with KSBRF ROM. Firstly, the relationship between input and output of a nonlinear system was built using KSBRF ROM. Then, the identification of nonlinear unsteady aerodynamic forces was conducted when an airfoil had pitch/plunge oscillations. It was shown that under the same free-stream’s Mach number, the prediction results with KSBRF ROM agree well with those with CFD, the mean prediction errors of drag coefficient and pitch moment coefficient are within 2.0%; considering a variable free-stream Mach number, the mean prediction errors of lift coefficient and pitch moment coefficient are within 2.5%, the preceding error of drag coefficient is less than 7%; it is concluded that considering the history effect is helpful to improve the accuracy of ROM through studying the effects of differentm,nvalues on the accuracy of ROM; the LCO characteristics of NACA64A010 airfoil are predicted with the nonlinear ROM, the prediction results with KSBRF ROM agree well with those with CFD, the errors are less than 5.17%.

surrogate model; reduced-order model; kriging interpolation; unsteady aerodynamic force; limit cycle oscillations

10.13465/j.cnki.jvs.2016.13.023

國(guó)家“973”計(jì)劃(2014CB744804)

2014-12-30修改稿收到日期：2015-06-15

劉艷 女，博士生，1988年3月生

白俊強(qiáng) 男，教授，1971年1月生

E-mail:junqiang@nwpu.edu.cn

V215.3+2；V211.3

A