基于隨機森林算法的非定常氣動力建模研究

徐旺丁， 張 兵， 王華畢

（合肥工業大學 機械工程學院，合肥230009）

1 引 言

非定常氣動力建模從理論研究到基于CFD（Computational Fluid Dynamics）技術的非定常流場數值模擬，再到基于CFD技術的非定常氣動力降階模型，是氣動彈性力學研究的重要組成部分和關鍵技術之一?；贑FD技術的氣動彈性分析對氣體流動的刻畫越細致，則氣動彈性系統的維數就越高，從而導致氣動彈性模擬計算時間耗費巨大，特別是涉及氣動彈性優化問題時，有時能達到難以接受的地步。因此，發展非定常氣動力降階模型ROM（Reduced Order Model）已成為一個熱門的研究領域?；跉鈩恿惦A模型的氣動彈性分析，可以快速求解該氣動彈性系統對任意輸入的響應，僅需利用CFD求解器求解給定狀態下流場的輸入輸出特性進行構造，無需對程序進行大的改動。

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

Volterra／ROM在模擬較強非線性的氣動彈性現象時會遇到模型維數災難的問題；而ARMA／ROM是線性系統辨識模型，存在難以對強非線性特性進行建模的理論限制，且需要反復調制輸入信號來確定激勵特征頻率范圍；以神經網絡等為代表的響應面模型則難以建立傳遞函數或狀態空間模型，且易出現過學習或欠學習的困難，限制了其應用范圍［9］；支持向量機在其訓練階段運算量非常大，特別是對于非線性氣動彈性分析，需要采用的支持向量機數目非常多，造成計算量急劇增大。

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

2 隨機森林的基本原理

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

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

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

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

圖1 隨機森林算法Fig.1 Diagram of random forest algorithm

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

設A為m×m的相關矩陣，且

式中E為單位矩陣。對有限個訓練樣本近似估計，則

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

最優權值δ可通過求解最優化問題式（7）得到，

由拉格朗日定理可知，

式中λ為拉格朗日乘子。獲得最優條件：

式中 1v＝［1，1，…，1］T。解得最優解為

3 基于RF／ROM的氣動彈性建模

3.1 結構模型

為驗證RF／ROM，選擇 NACA0012和 NACA64A010翼型模型，模型結構示意圖如圖2所示。結構動力學系統包含繞彈性軸的俯仰運動α（順時針為正）和y方向的沉浮運動h（向下為正）兩個自由度。結構動力學方程如下。

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

式中 無量綱速度v＊＝V／（ωαbμ1／2），V 為來流速度，質量比μ＝m／（πρb2），h＝h／b為結構響應廣義位移，rα為翼型圍繞彈性軸的回轉半徑，ωh和ωα分別為彎曲和扭轉模態的固有頻率，CL和CM分別為升力系數和俯仰運動力矩系數。

3.2 基于CFD技術的ROM訓練過程

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

由于非定常流場時間具有延遲效應，所以可以選擇xi＝［（ui，vi），（ui－1，vi－1），…，（ui－r，vi－r），fi－1，…，fi－s）作為輸入信號，r和s是由用戶選擇的時間延遲。因此，可以從非定常CFD模擬的結果構建RF／ROM 的新訓練樣本集（xi，fi）。

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

3.3 RF／ROM的構建

基于RF的氣動彈性ROM的構建流程如圖4所示，共有5個主要步驟。（1）設計多組輸入信號，并通過非定常CFD求解器計算相應的非線性氣動力響應；（2）將輸入信號和非定常氣動力響應組合成訓練樣本集（xi，fi）；（3）執行 RF算法獲得辨識模型RF／ROM，并用測試集驗證模型的精度；（4）將RF／ROM應用于預測不同工況下的非定常氣動力；（5）耦合結構動力學系統，根據圖4的虛線路徑，將RF／ROM用于氣動彈性響應的預測。

4 算例驗證

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

4.1 顫振邊界預測

圖3 訓練輸入信號Fig.3 Training input signal

圖4 RF／ROM構建流程Fig.4 Workflow of RF-based ROM

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

在0.71Ma工況下，訓練輸出的結果如圖5所示，可以看出，RF／ROM預測升力系數和力矩系數的精度很高。為進一步驗證RF／ROM的正確性，使用簡諧信號輸入下俯仰運動的氣動力系數作為測試信號，圖6為在最后一個周期中，ROM預測結果和CFD計算結果的對比。

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

圖5 NACA0012翼型CFD直接模擬和RF／ROM預測的輸出對比Fig.5 Results of NACA0012airfoil comparison between CFD and RF／ROM

對5個工況下的翼型顫振邊界進行計算，計算結果和CFD直接模擬計算以及文獻［16］的試驗值進行對比，如圖8所示?？梢钥闯?，RF／ROM的預測結果和CFD直接模擬計算的結果吻合良好，表明RF／ROM預測顫振邊界特性具有較高的精度。然而，模擬值和實驗值匹配不佳的主要原因是部分結構參數，如重心位置等對顫振結果存在顯著的影響，且實測數據存在著一定的誤差，因此顫振實驗值的離散度較大［17］。

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

圖6 最后一個周期氣動力系數對比Fig.6 Comparison of aerodynamic coefficients for the last cycle

4.2 LCO預測

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

在 Ma＝0.8，v＊＝0.7工況下，給定初始沉浮運動和俯仰運動速度值均為0.1，CFD耦合計算LCO特性和RF／ROM預測LCO特性的結果對比如圖9所示。

為了進一步驗證RF／ROM預測LCO特性的性能，計算了馬赫數為0.76，0.8，0.825，0.85和0.875共5種情況下的氣動彈性系統的LCO特性，所有馬赫數情況下均取v＊＝0.75。運用CFD方法和RF／ROM方法計算的LCO振幅如圖10所示，并引入參考文獻［19］的CFD和SVM／ROM計算結果作為參考和對比?？梢钥闯觯琑F／ROM和CFD計算的結果很一致，平均誤差在1 0%以內。CFD方法捕獲到LCO振幅需要大約1小時，而RF／ROM方法預測相同的LCO特性不超過5分鐘。故RF／ROM的計算效率較高，這對實時飛行模擬和飛行控制器的設計有重要意義。

圖7 0.71Ma，不同動壓下的NACA0012翼型結構響應Fig.7 0.71Ma，NACA0012airfoil structure response under different dynamic pressure

圖8 NACA0012翼型的顫振邊界計算結果Fig.8 Flutter boundary of the NACA0012airfoil

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

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

5 結 語

基于CFD技術，將隨機森林算法引入非定常氣動力降階模型建模研究領域，詳細說明了RF／ROM的構建過程，采用所構建的RF／ROM實現對二維翼型模型的顫振邊界和LCO特性的預測，并與CFD直接耦合計算的結果進行對比。結果表明，RF／ROM可以用于非定常流場中的氣動彈性特性預測，計算效率也提高了10～20倍。如何進一步提高ROM的精度和效率，特別是對具有變參數的模型，將會成為未來氣動力建模研究領域的重點。