徐旺丁， 張 兵， 王華畢

（合肥工業(yè)大學(xué) 機(jī)械工程學(xué)院，合肥230009）

1 引 言

非定常氣動(dòng)力建模從理論研究到基于CFD（Computational Fluid Dynamics）技術(shù)的非定常流場數(shù)值模擬，再到基于CFD技術(shù)的非定常氣動(dòng)力降階模型，是氣動(dòng)彈性力學(xué)研究的重要組成部分和關(guān)鍵技術(shù)之一。基于CFD技術(shù)的氣動(dòng)彈性分析對(duì)氣體流動(dòng)的刻畫越細(xì)致，則氣動(dòng)彈性系統(tǒng)的維數(shù)就越高，從而導(dǎo)致氣動(dòng)彈性模擬計(jì)算時(shí)間耗費(fèi)巨大，特別是涉及氣動(dòng)彈性優(yōu)化問題時(shí)，有時(shí)能達(dá)到難以接受的地步。因此，發(fā)展非定常氣動(dòng)力降階模型ROM（Reduced Order Model）已成為一個(gè)熱門的研究領(lǐng)域?；跉鈩?dòng)力降階模型的氣動(dòng)彈性分析，可以快速求解該氣動(dòng)彈性系統(tǒng)對(duì)任意輸入的響應(yīng)，僅需利用CFD求解器求解給定狀態(tài)下流場的輸入輸出特性進(jìn)行構(gòu)造，無需對(duì)程序進(jìn)行大的改動(dòng)。

目前，基于系統(tǒng)辨識(shí)的模型降階方法為非定常氣動(dòng)力建模提供了新思路。系統(tǒng)辨識(shí)根據(jù)系統(tǒng)的輸入輸出數(shù)據(jù)來確定描述系統(tǒng)行為的數(shù)學(xué)模型，其基本原理和依據(jù)是現(xiàn)代系統(tǒng)辨識(shí)理論，對(duì)系統(tǒng)進(jìn)行分析的主要問題是根據(jù)輸入信號(hào)和系統(tǒng)的特性來確定輸出信號(hào)。在氣動(dòng)彈性系統(tǒng)中，較為常見的基于系統(tǒng)辨識(shí)的降階方法有 Volterra級(jí)數(shù)法［1，2］、神經(jīng)網(wǎng)絡(luò)模 型［3］、ARMA（Auto-Regressive Moving Average）模 型［4，5］和 支 持 向 量 機(jī) SVM （Support Vector Machines）方法［6］等。非定常氣動(dòng)力 ROM和結(jié)構(gòu)模型可以耦合得到整個(gè)氣動(dòng)彈性系統(tǒng)ROM，因此可以方便地用于氣動(dòng)彈性領(lǐng)域的顫振分析［7，8］。

Volterra／ROM在模擬較強(qiáng)非線性的氣動(dòng)彈性現(xiàn)象時(shí)會(huì)遇到模型維數(shù)災(zāi)難的問題；而ARMA／ROM是線性系統(tǒng)辨識(shí)模型，存在難以對(duì)強(qiáng)非線性特性進(jìn)行建模的理論限制，且需要反復(fù)調(diào)制輸入信號(hào)來確定激勵(lì)特征頻率范圍；以神經(jīng)網(wǎng)絡(luò)等為代表的響應(yīng)面模型則難以建立傳遞函數(shù)或狀態(tài)空間模型，且易出現(xiàn)過學(xué)習(xí)或欠學(xué)習(xí)的困難，限制了其應(yīng)用范圍［9］；支持向量機(jī)在其訓(xùn)練階段運(yùn)算量非常大，特別是對(duì)于非線性氣動(dòng)彈性分析，需要采用的支持向量機(jī)數(shù)目非常多，造成計(jì)算量急劇增大。

隨機(jī)森林RF（Random Forest）算法是一種比較新的機(jī)器學(xué)習(xí)算法，其計(jì)算開銷小，具有較好的大數(shù)據(jù)集回歸預(yù)測(cè)效果。根據(jù)大數(shù)據(jù)集理論，隨機(jī)森林不存在過擬合的問題［10，11］，在很多現(xiàn)實(shí)任務(wù)中展現(xiàn)出強(qiáng)大的性能，是代表集成學(xué)習(xí)技術(shù)水平的方法［12］，該方法在航空航天領(lǐng)域已有較好的應(yīng)用［13，14］。本文嘗試將RF算法應(yīng)用于建立非定常氣動(dòng)力模型ROM，并將所得到的模型命名為RF／ROM。RF／ROM也是一種基于系統(tǒng)辨識(shí)的氣動(dòng)力建模方法，采用CFD方法計(jì)算訓(xùn)練信號(hào)的非定常氣動(dòng)力，僅以一組訓(xùn)練樣本對(duì)RF模型進(jìn)行訓(xùn)練，采用簡諧信號(hào)對(duì)該RF模型進(jìn)行測(cè)試，最后將測(cè)試結(jié)果與CFD仿真結(jié)果對(duì)比，檢驗(yàn)RF／ROM對(duì)非定常氣動(dòng)力的預(yù)測(cè)精度和效率，并將RF／ROM應(yīng)用于顫振邊界以及極限環(huán)振蕩（Limit Cycle Oscillation，LCO）特性的預(yù)測(cè)。

2 隨機(jī)森林的基本原理

隨機(jī)森林是由多顆決策樹組成的集成學(xué)習(xí)模型，在處理回歸問題時(shí)，隨機(jī)森林以回歸樹為基學(xué)習(xí)器。隨機(jī)森林是在Bagging的基礎(chǔ)上引入了隨機(jī)屬性選擇的思想，回歸效果優(yōu)于Bagging。根據(jù)大數(shù)據(jù)集理論，隨機(jī)森林不存在過擬合的問題［10］。

令fi（X）（i＝1，2，…，m）表示m 個(gè)子模型，δi表示fi（X）的權(quán)值，那么集成模型fE（X）可表示為

隨機(jī)森林算法如圖1所示，先采用Bootstrap方法從訓(xùn)練樣本集TTrain＝｛（Xk，Yk）中按照隨機(jī)的、有放回的和重新選擇的方式，產(chǎn)生m個(gè)相互獨(dú)立的樣本子集Ti＝ ｛（Xk，Yk）（i＝1，2，…，m）。對(duì)于任意的Ti，NB（NB＜N）個(gè)樣本服從統(tǒng)一分布且相互獨(dú)立。然后，在每個(gè)Ti上都建立一個(gè)回歸樹子模型，并通過隨機(jī)屬性選擇的方式分裂根結(jié)點(diǎn)或父結(jié)點(diǎn)。對(duì)于任意樣本X，m個(gè)子模型將產(chǎn)生m 個(gè)預(yù)測(cè)值P1，P2，…，Pm，最后將m個(gè)預(yù)測(cè)值去平均值后作為隨機(jī)森林的預(yù)測(cè)值PE。

從式（2）可以看出，所有子模型的權(quán)值相同，δi＝1／m（i＝1，2，…，m），但實(shí)際上子模型的樣本子集不同，準(zhǔn)確度也不同。如果給模型分配相同的權(quán)值，則會(huì)降低隨機(jī)森林的預(yù)測(cè)精度。因此需要對(duì)子模型進(jìn)行融合處理，本文采用GEM（Generalized Ensemble Method）算法［15］融合子模型。該算法是以fE（X）的均方誤差最小為目標(biāo)函數(shù)，作為約束條件，通過子模型的偏差構(gòu)造相關(guān)矩陣，從而直接求解權(quán)值向量。GEM算法具體過程如下。

圖1 隨機(jī)森林算法Fig.1 Diagram of random forest algorithm

令g（X）代表真實(shí)函數(shù)，di（X）＝fi（X）－g（X）代表子模型與真實(shí)函數(shù)的偏差函數(shù)，則集成模型還可以表示為

設(shè)A為m×m的相關(guān)矩陣，且

式中E為單位矩陣。對(duì)有限個(gè)訓(xùn)練樣本近似估計(jì)，則

則集成模型的誤差可表示為

最優(yōu)權(quán)值δ可通過求解最優(yōu)化問題式（7）得到，

由拉格朗日定理可知，

式中λ為拉格朗日乘子。獲得最優(yōu)條件：

式中 1v＝［1，1，…，1］T。解得最優(yōu)解為

3 基于RF／ROM的氣動(dòng)彈性建模

3.1 結(jié)構(gòu)模型

為驗(yàn)證RF／ROM，選擇 NACA0012和 NACA64A010翼型模型，模型結(jié)構(gòu)示意圖如圖2所示。結(jié)構(gòu)動(dòng)力學(xué)系統(tǒng)包含繞彈性軸的俯仰運(yùn)動(dòng)α（順時(shí)針為正）和y方向的沉浮運(yùn)動(dòng)h（向下為正）兩個(gè)自由度。結(jié)構(gòu)動(dòng)力學(xué)方程如下。

式中 Sα＝mxαb為對(duì)彈性軸的靜矩，fy為沉浮方向的外力合力，Mz為對(duì)z軸的合力矩。定義無量綱時(shí)間τ＝ωαt，二維翼型氣動(dòng)彈性方程可以寫成無量綱形式為

式中 無量綱速度v＊＝V／（ωαbμ1／2），V 為來流速度，質(zhì)量比μ＝m／（πρb2），h＝h／b為結(jié)構(gòu)響應(yīng)廣義位移，rα為翼型圍繞彈性軸的回轉(zhuǎn)半徑，ωh和ωα分別為彎曲和扭轉(zhuǎn)模態(tài)的固有頻率，CL和CM分別為升力系數(shù)和俯仰運(yùn)動(dòng)力矩系數(shù)。

3.2 基于CFD技術(shù)的ROM訓(xùn)練過程

RF／ROM的訓(xùn)練過程主要是基于CFD求解器的耦合計(jì)算和設(shè)計(jì)來訓(xùn)練輸入。對(duì)于CFD耦合計(jì)算，非定常氣動(dòng)力的計(jì)算需要花費(fèi)大量的計(jì)算成本，引入降階模型用于非定常流場的模擬可以提高氣動(dòng)彈性計(jì)算的效率。在氣動(dòng)彈性系統(tǒng)中，結(jié)構(gòu)狀態(tài)值u和v作為輸入值，廣義氣動(dòng)力系數(shù)向量f作為輸出，采樣時(shí)間序列數(shù)據(jù)（u，v，f）可以由CFD耦合求計(jì)算。如兩個(gè)自由度的氣動(dòng)彈性系統(tǒng)的輸入為（h，α），輸出為CL和CM。為了構(gòu)建圖1所示的基于RF的ROM，需要精心準(zhǔn)備采樣數(shù)據(jù)作為訓(xùn)練樣本集。一個(gè)高質(zhì)量的RF／ROM不僅需要精確的CFD訓(xùn)練計(jì)算，還要求訓(xùn)練輸入對(duì)系統(tǒng)的動(dòng)態(tài)特性進(jìn)行充分激勵(lì)。因此，需要精心設(shè)計(jì)訓(xùn)練輸入信號(hào)以保證激勵(lì)具有足夠的頻帶寬度。圖3為本文設(shè)計(jì)的以一種過濾的高斯白噪聲信號(hào)作為輸入訓(xùn)練信號(hào)，因?yàn)榘自肼曅盘?hào)可以很好地代表非線性系統(tǒng)的自然動(dòng)態(tài)特性。該信號(hào)是通過采用數(shù)字生成的高斯分布隨機(jī)時(shí)間序列來創(chuàng)建，并對(duì)其進(jìn)行濾波處理，可以看出，該訓(xùn)練信號(hào)既含有大振幅成分，也含有小振幅成分，頻譜也比較寬，適合作為RF／ROM的訓(xùn)練輸入信號(hào)。

由于非定常流場時(shí)間具有延遲效應(yīng)，所以可以選擇xi＝［（ui，vi），（ui－1，vi－1），…，（ui－r，vi－r），fi－1，…，fi－s）作為輸入信號(hào)，r和s是由用戶選擇的時(shí)間延遲。因此，可以從非定常CFD模擬的結(jié)果構(gòu)建RF／ROM 的新訓(xùn)練樣本集（xi，fi）。

圖2 二維翼型氣動(dòng)彈性模型Fig.2 Diagram of two-DOF aeroelastic model

3.3 RF／ROM的構(gòu)建

基于RF的氣動(dòng)彈性ROM的構(gòu)建流程如圖4所示，共有5個(gè)主要步驟。（1）設(shè)計(jì)多組輸入信號(hào)，并通過非定常CFD求解器計(jì)算相應(yīng)的非線性氣動(dòng)力響應(yīng)；（2）將輸入信號(hào)和非定常氣動(dòng)力響應(yīng)組合成訓(xùn)練樣本集（xi，fi）；（3）執(zhí)行 RF算法獲得辨識(shí)模型RF／ROM，并用測(cè)試集驗(yàn)證模型的精度；（4）將RF／ROM應(yīng)用于預(yù)測(cè)不同工況下的非定常氣動(dòng)力；（5）耦合結(jié)構(gòu)動(dòng)力學(xué)系統(tǒng)，根據(jù)圖4的虛線路徑，將RF／ROM用于氣動(dòng)彈性響應(yīng)的預(yù)測(cè)。

4 算例驗(yàn)證

為了驗(yàn)證氣動(dòng)力RF／ROM的精確性和高效性，應(yīng)用3.2節(jié)設(shè)計(jì)的過濾的高斯白噪聲信號(hào)位移輸入對(duì)CFD求解器進(jìn)行激勵(lì)，選擇r＝5和s＝4為RF／ROM的準(zhǔn)備訓(xùn)練樣本集，根據(jù)3.3節(jié)流程建立RF／ROM。本文主要做了兩方面工作，一是對(duì)NACA0012翼型的不同馬赫數(shù)工況下的顫振邊界進(jìn)行了計(jì)算，并將RF／ROM的計(jì)算結(jié)果和CFD仿真結(jié)果以及試驗(yàn)結(jié)果進(jìn)行對(duì)比；二是將RF／ROM用于NACA64A010翼型的跨音速極限環(huán)顫振特性預(yù)測(cè)，并與CFD仿真結(jié)果對(duì)比。

4.1 顫振邊界預(yù)測(cè)

圖3 訓(xùn)練輸入信號(hào)Fig.3 Training input signal

圖4 RF／ROM構(gòu)建流程Fig.4 Workflow of RF-based ROM

顫振邊界算例驗(yàn)證選用的NACA0012機(jī)翼結(jié)構(gòu)模型、來流Ma、Re以及顫振邊界實(shí)驗(yàn)數(shù)據(jù)來自文獻(xiàn)［16］，選 擇 0.45Ma，0.61Ma，0.71Ma，0.77Ma和0.8Ma共5個(gè)工況進(jìn)行顫振邊界分析計(jì)算，初始迎角為0°。模型結(jié)構(gòu)參數(shù)為m＝87.07kg，c＝0.4064m，Kh＝38821N／m，Kα＝3928N／m，a＝0，xα＝0，Iα＝3.68kg·m2。

在0.71Ma工況下，訓(xùn)練輸出的結(jié)果如圖5所示，可以看出，RF／ROM預(yù)測(cè)升力系數(shù)和力矩系數(shù)的精度很高。為進(jìn)一步驗(yàn)證RF／ROM的正確性，使用簡諧信號(hào)輸入下俯仰運(yùn)動(dòng)的氣動(dòng)力系數(shù)作為測(cè)試信號(hào)，圖6為在最后一個(gè)周期中，ROM預(yù)測(cè)結(jié)果和CFD計(jì)算結(jié)果的對(duì)比。

使用已建立的RF／ROM可以方便地開展高效顫振邊界預(yù)測(cè)，利用變密度的方式調(diào)整來流動(dòng)壓，可以快速預(yù)測(cè)氣動(dòng)彈性系統(tǒng)的顫振邊界。如在Ma＝0.71工況下，計(jì)算動(dòng)壓q分別設(shè)為0.95qe，0.85qe和0.9qe時(shí)（qe為顫振動(dòng)壓的實(shí)驗(yàn)值），翼型的結(jié)構(gòu)響應(yīng)如圖7所示。圖7（a）中廣義位移的幅值逐漸增大，說明此時(shí)結(jié)構(gòu)振動(dòng)是發(fā)散的，即已經(jīng)發(fā)生顫振；圖7（b）中廣義位移的幅值逐漸減小，說明此時(shí)結(jié)構(gòu)振動(dòng)是收斂的，尚未達(dá)到顫振邊界；圖7（c）中廣義位移表現(xiàn)為等幅振蕩，正好位于顫振臨界邊界，此時(shí)的動(dòng)壓為計(jì)算顫振動(dòng)壓，相應(yīng)的振動(dòng)頻率為計(jì)算顫振頻率。

圖5 NACA0012翼型CFD直接模擬和RF／ROM預(yù)測(cè)的輸出對(duì)比Fig.5 Results of NACA0012airfoil comparison between CFD and RF／ROM

對(duì)5個(gè)工況下的翼型顫振邊界進(jìn)行計(jì)算，計(jì)算結(jié)果和CFD直接模擬計(jì)算以及文獻(xiàn)［16］的試驗(yàn)值進(jìn)行對(duì)比，如圖8所示。可以看出，RF／ROM的預(yù)測(cè)結(jié)果和CFD直接模擬計(jì)算的結(jié)果吻合良好，表明RF／ROM預(yù)測(cè)顫振邊界特性具有較高的精度。然而，模擬值和實(shí)驗(yàn)值匹配不佳的主要原因是部分結(jié)構(gòu)參數(shù)，如重心位置等對(duì)顫振結(jié)果存在顯著的影響，且實(shí)測(cè)數(shù)據(jù)存在著一定的誤差，因此顫振實(shí)驗(yàn)值的離散度較大［17］。

在上述算例中，其中一種馬赫數(shù)工況下，CFD直接模擬顫振特性通常需要進(jìn)行4～5次的CFD／CSD耦合計(jì)算方可得到顫振臨界點(diǎn)，每次CFD／CSD耦合計(jì)算需要計(jì)算5～8個(gè)周期，約1250～2000個(gè)時(shí)間步，共需5000～10000個(gè)時(shí)間步。而對(duì)于RF／ROM方法，在某種馬赫數(shù)工況下，只需要進(jìn)行一次訓(xùn)練過程，約500個(gè)時(shí)間步，除此之外的RF／ROM計(jì)算時(shí)間可以忽略不計(jì)。顯然，在定馬赫數(shù)工況下，基于RF／ROM的顫振計(jì)算效率提高了約10～20倍。氣動(dòng)網(wǎng)格數(shù)目越多，則CFD計(jì)算量越大，RF／ROM方法的效率越明顯。

圖6 最后一個(gè)周期氣動(dòng)力系數(shù)對(duì)比Fig.6 Comparison of aerodynamic coefficients for the last cycle

4.2 LCO預(yù)測(cè)

由于非線性因素，彈性結(jié)構(gòu)在流場中的運(yùn)動(dòng)形式通常表現(xiàn)為極限環(huán)振蕩。針對(duì)因跨音速氣動(dòng)力非線性而產(chǎn)生的極限環(huán)型顫振分析研究［18］，選擇NACA64A010翼型作為結(jié)構(gòu)模型，結(jié)構(gòu)參數(shù)為xα＝0.25，r2α＝0.75，a＝－0.6，μ＝75，ωh／ωα＝0.5。先運(yùn)用CFD方法計(jì)算該模型的LCO特性，然后建立RF／ROM，并開展LCO特性計(jì)算，將兩種方法的計(jì)算結(jié)果進(jìn)行對(duì)比，驗(yàn)證RF／ROM預(yù)測(cè)LCO的正確性和計(jì)算精度。

在 Ma＝0.8，v＊＝0.7工況下，給定初始沉浮運(yùn)動(dòng)和俯仰運(yùn)動(dòng)速度值均為0.1，CFD耦合計(jì)算LCO特性和RF／ROM預(yù)測(cè)LCO特性的結(jié)果對(duì)比如圖9所示。

為了進(jìn)一步驗(yàn)證RF／ROM預(yù)測(cè)LCO特性的性能，計(jì)算了馬赫數(shù)為0.76，0.8，0.825，0.85和0.875共5種情況下的氣動(dòng)彈性系統(tǒng)的LCO特性，所有馬赫數(shù)情況下均取v＊＝0.75。運(yùn)用CFD方法和RF／ROM方法計(jì)算的LCO振幅如圖10所示，并引入?yún)⒖嘉墨I(xiàn)［19］的CFD和SVM／ROM計(jì)算結(jié)果作為參考和對(duì)比?？梢钥闯?，RF／ROM和CFD計(jì)算的結(jié)果很一致，平均誤差在1 0%以內(nèi)。CFD方法捕獲到LCO振幅需要大約1小時(shí)，而RF／ROM方法預(yù)測(cè)相同的LCO特性不超過5分鐘。故RF／ROM的計(jì)算效率較高，這對(duì)實(shí)時(shí)飛行模擬和飛行控制器的設(shè)計(jì)有重要意義。

圖7 0.71Ma，不同動(dòng)壓下的NACA0012翼型結(jié)構(gòu)響應(yīng)Fig.7 0.71Ma，NACA0012airfoil structure response under different dynamic pressure

圖8 NACA0012翼型的顫振邊界計(jì)算結(jié)果Fig.8 Flutter boundary of the NACA0012airfoil

圖9 NACA64A010翼型在 Ma＝0.8，v＊＝0.7工況下的LCO響應(yīng)對(duì)比Fig.9 Comparison of LCO response at Ma＝0.8，v＊＝0.7

圖10 NACA64A010翼型在不同馬赫數(shù)工況下LCO振幅對(duì)比Fig.10 Comparison of LCO amplitude with Mach number

5 結(jié) 語

基于CFD技術(shù)，將隨機(jī)森林算法引入非定常氣動(dòng)力降階模型建模研究領(lǐng)域，詳細(xì)說明了RF／ROM的構(gòu)建過程，采用所構(gòu)建的RF／ROM實(shí)現(xiàn)對(duì)二維翼型模型的顫振邊界和LCO特性的預(yù)測(cè)，并與CFD直接耦合計(jì)算的結(jié)果進(jìn)行對(duì)比。結(jié)果表明，RF／ROM可以用于非定常流場中的氣動(dòng)彈性特性預(yù)測(cè)，計(jì)算效率也提高了10～20倍。如何進(jìn)一步提高ROM的精度和效率，特別是對(duì)具有變參數(shù)的模型，將會(huì)成為未來氣動(dòng)力建模研究領(lǐng)域的重點(diǎn)。