基于增量學習的非定常氣動力參數化降階模型

陳志強，劉戰合，苗楠，馮偉

1.鄭州航空工業管理學院 航空工程學院，鄭州 450046 2.河南工業大學 機電工程學院，鄭州 450046

在飛機氣動彈性設計和分析中，需研究結構力、氣動力和慣性力的耦合，以確定飛行包線。由于流固相互作用對安全與效率的要求有著至關重要的影響，因此可靠、準確地預測氣動力顯得尤為重要。考慮到工業應用，非定常氣動力的建模主要采用高效的勢流方法。然而，在跨聲速飛行情況下，氣動力非線性現象表現明顯，經典方法在保真度上受到限制。使用現代計算流體力學(Computational Fluid Dynamics，CFD)方法可以提高精度和保真度，然而由于影響參數(結構特征模態、馬赫數、迎角和頻率)的廣泛變化，數值計算耗時大大增加。目前漫長的計算周期仍然限制了CFD方法在多學科計算方面的廣泛應用。

自20世紀90年代以來，研究者們提出了各種基于CFD仿真數據的氣動力降階模型[1-4](Reduced Order Modeling，ROM)。ROM是指采用線性或非線性模型來表達復雜氣動力系統中輸入與輸出間的映射關系。非定常氣動力ROM和結構模型可方便地進行耦合，進而得到整個氣動彈性系統ROM，可大幅提高跨聲速氣動彈性問題的數值仿真效率。根據數學方法、模型構造和數據源的不同，降階方法大致可分為兩類[3]：一類是基于投影方法提取流場的主要特征，例如本征正交分解(Proper Orthogonal Decomposition，POD)方法[5-6]、動態模態分解[7]、稀疏識別理論[8]、Koopman理論[9]等，另一類是基于系統識別方法反映流場非定常氣動力的輸入輸出特性，例如特征實現算法[10]、諧波平衡法[11]、神經網絡[12-14]、支持向量機[15-16]以及代理模型[17]等。

目前，已有多種ROM應用于跨聲速氣動彈性數值模擬[18-21]，能夠大幅提高計算效率，然而大多數ROM技術對飛行參數和結構參數缺乏魯棒性。在每次物理或建模參數變化時需要重新構建相應的ROM，但是在許多情況下，重建ROM需要從一組新的高保真度數值仿真獲得新系統的響應，對于追求高效率的目標來說，這種做法是不可行的。對于跨聲速流動，其穩定性對迎角、馬赫數甚至雷諾數等都很敏感，故跨聲速氣動彈性特性易受這些參數的影響，導致有時會存在較大的分散度和不確定性[22]。因此構建精度高、效率高的跨聲速非定常氣動力參數化降階模型對跨聲速氣動彈性分析和氣動彈性主動控制研究尤為重要。針對兩類降階方法，近些年國內外學者對ROM存在魯棒性差問題進行了研究。

特征提取方法是利用復雜動力學系統的密集數據提取特征模態，通過特征模態上的低階動力學行為來近似復雜系統動力學行為。對于流場特征結構模型, 需要發展具有自適應性和魯棒性的 ROM構造方法, 并通過并行化來提高計算效率。為了解決傳統氣動力降階方法缺乏自適應性和魯棒性的問題，Amsallem等[23-24]提出了一種基于微分幾何概念和結論的ROM，其能較好適應于物理或建模參數的變化，并且該方法適用于任何需要重構基向量的ROM。該方法在全機的氣動彈性分析中，考察了其對自由流馬赫數變化的適應性。與直接重建的降階模型、高保真非線性和線性仿真結果比較具有良好的一致性，證明了該方法具有較低的復雜度和較好的魯棒性。Chen等[25]提出了一種基于自適應POD和最小二乘支持向量回歸(Least Squares Support Vector Regression, LS-SVR)相結合的非定常氣動參數ROM方法。在給定的參數空間中，利用對Grassmann流形在切空間內的插值方法，生成了任意樣本點的模態矩陣。再通過LS-SVR得到激勵信號與相應POD系數之間的動態關系。與直接CFD模擬相比，該參數化ROM能夠有效地預測氣動力響應、顫振邊界和極限周期振蕩。

系統識別方法是利用流體系統的輸入和輸出關系來建立系統的一個低階傳遞函數或狀態空間模型代替原始的滿階模型，為高效、高精度模擬飛行器跨聲速氣動彈性力學與控制問題提供了解決途徑。基于系統辨識的ROM只需利用給定狀態下非定常CFD程序求解流場的輸入輸出特性進行構造，而無需對程序進行大的改動，方便研究者使用。傳統的系統識別方法在處理較強非線性和狀態參數變化等問題時通常顯得無能為力，因此近些年來針對上述問題，一些學者開始了相關研究[26-30]。Winter和Breitsamber[31]提出了一種基于模糊神經網絡模型的非線性系統辨識方法，從而構建非定常氣動參數化ROM。其用局部線性模型樹算法訓練了一個三層神經網絡使其適用于較寬馬赫數范圍，并分析了AGARD445.6機翼從亞聲速到高超聲速的顫振邊界，但該邊界不涉及極限環振蕩(Limit Cycle Oscillation, LCO)行為。Li等[32]提出了一種基于深度學習理論中長短期記憶(Long Short-Term Memory, LSTM)網絡的大訓練數據集和空間樣本的非定常空氣動力學模型。LSTM網絡在建模時間序列數據方面具有很大的潛力，非常適合捕捉非定常空氣動力學的時滯效應。與傳統的降階模型不同，基于LSTM網絡的模型不需要選擇延遲階數。計算結果表明，該模型能夠準確地描述氣動和氣動彈性系統在變流和變結構參數下的動態特性。

根據以上討論，本文提出了基于最小二乘支持向量回歸算法的跨聲速非定常氣動力參數化降階模型，建立結構激勵、飛行參數和相應的非定常氣動力之間的動態關系。為了提高模型的準確度，定義系統輸入包含當前激勵信號、適當的延遲激勵信號和響應反饋。為了解決模型構建時精度不足或需拓寬參數范圍的問題，采用了一種增量學習算法，算法執行中新增樣本后的矩陣求逆采用了迭代方法，不用對整個訓練數據集重新訓練，進而很大地提高了建模效率。為了驗證該方法的有效性和精確度，將該方法應用于預測NACA64A010翼型的跨聲速非定常氣動力響應和氣動彈性響應來評估其性能。

1 參數化降階建模方法

非定常氣動力模型可以看作一個多輸入多輸出的動力學系統，表示為

(1)

式中：k表示當前離散時間步；狀態函數f(k+1)表示當前狀態x(k)和輸入向量u(k)與新狀態x(k+1)的關系；輸出函數g(k+1)表示新狀態x(k+1)與系統輸出y(k+1)的關系；輸入向量u(t)可包含系統的參數、運動狀態等信息；輸出向量y(t)表示時域氣動力；狀態向量x(t)通常為流場信息，比如密度、速度、壓強等。基于CFD求解器的計算，x(k)的維數與流場網格點數成正比，通常可達到數十萬以上。式(1)表明全階CFD系統的計算量巨大，在此可由降階方法代替以達到高效且精確獲得氣動力的目的。

從線性系統識別理論可知，動態系統的輸出可以近似為當前激勵輸入、一定數量的延遲輸入和延遲系統響應的函數。對于降階模型，式(1)中的輸入和輸出關系可等價為非線性系統，表示為

yi(k)=Φi(u(k),u(k-1),…,u(k-m),

(2)

1.1 最小二乘法支持向量回歸

最小二乘支持向量回歸是在標準SVR基礎上進行了簡化求解過程的改進，把標準SVR中的不等式約束條件變為等式約束條件，將優化問題的求解變為線性方程組的求解，大大簡化了計算復雜性，適用于函數擬合和系統建模。

數據集通常具有非線性，利用非線性影射φ將輸入空間映射到高維特征空間, 使非線性函數擬合問題轉化為高維特征空間的線性函數擬合問題。回歸函數具有如下形式：

(3)

式中：w為特征空間的向量；b為偏置項。

回歸問題可表示為有約束條件的優化問題：

(4)

式中：優化目標函數J的兩項分別為模型的復雜度和誤差的范圍;ei是訓練數據誤差；C是懲罰因子，表示誤差項的懲罰程度大小。相應的Lagrange函數表示為

(5)

根據Karush-Kuhn-Tucker’s(KKT’s)最優條件：

(6)

得到如下線性方程：

(7)

(8)

LS-SVR的回歸模型表達式為

(9)

根據泛函的有關理論，任意滿足 Mercer條件的對稱函數可作為核函數，本文選擇高斯徑向基函數作為核函數。

1.2 增量學習算法

對于增量學習(Incremental Learning，IL)，隨著樣本集的增加，最小二乘法支持向量回歸中的核函數矩陣H、Lagrange乘子α和常值偏置項b都會發生變化。由式(8)可知，增量學習中關鍵問題是得到核函數矩陣H的逆H-1。由于非常氣動力數據量巨大，矩陣H的階數較高，所以求逆矩陣的計算量將很大。在此將利用分塊矩陣求逆的方法，利用原核函數矩陣的逆，快速得到需要的逆矩陣，進而快速完成增量學習。

設原訓練集中有M1個樣本，對應核函數矩陣H1，其階數為Nsp1，新增M2個樣本，對應核函數矩陣H2，其階數為Nsp2，可表示為

(10)

(11)

根據文獻[33]定理可知H1和H2可逆，則

(12)

實際上增量學習算法可采用離線學習的方式，在每次迭代過程中增加一個樣本數據, 直至包括所有的樣本數據，這樣可以充分利用前一次迭代的運算結果，減少離線學習的計算復雜性。

1.3 基于增量學習的降階建模

為了得到準確的降階模型，需要獲得充分的輸入輸出數據用于模型訓練。需要強調的是，只有由訓練數據激發的結構振幅和頻率范圍可以被ROM充分地復制。本文選擇濾波高斯白噪聲(Filtered White Gaussian Noise, FWGN)作為對氣動力系統的激勵信號，可通過選擇合適的時間步長、信號長度和幅值范圍，將氣動力系統的動特性激發出來。進行每個飛行工況下的非定常氣動力CFD計算，FWGN信號可同時激勵各個結構自由度，快速生成訓練數據。

基于上述最小二乘法支持向量回歸的增量學習算法，構建降階模型以高效獲取不同飛行參數和振動狀態下的非定常氣動力。建模步驟可分為：① 在定義的飛行參數空間內選擇一定數量的樣本點，在設定的FWGN信號激勵下獲得CFD計算的氣動力數據；② 確定參數m和n的值，收集所有樣本點處的輸入和輸出數據，根據要求整理得到訓練數據集；③ 基于訓練數據集，通過最小二乘法支持向量回歸算法訓練得到降階模型SVR-ROM；④ 拓展參數范圍，增加樣本點，根據第 ①② 步方法得到新訓練數據集；⑤基于新增訓練數據集，通過增量學習算法訓練得到新的降階模型IL-ROM；⑥ 測試得到的降階模型在定義的參數空間內的可靠性和精度。

為了說明該降階方法的精度，下面定義一個相對誤差指標：

(13)

式中：yROM是降階模型的預測值；yCFD是由CFD仿真得到的結果。

2 算例分析

2.1 模型描述

圖1 兩自由度二元翼型NACA64A010模型

α(t)=α0+αAsin(ωα·t)

(14)

式中：α0為平均迎角；αA為俯仰運動幅值。俯仰運動的頻率ωα由減縮頻率kr確定，其中減縮頻率kr與頻率ωα之間的關系為

(15)

式中：U∞為無窮遠來流的速度。

本文所用到訓練數據和測試數據均由CFD計算得到，使用的求解器為NASA Langley研究中心開發的開源計算工具CFL3D。CFL3D代碼使用單元中心格式的隱式有限體積算法來求解依賴時間的雷諾數平均Navier-Stokes方程。該算法基于單元中心格式，采用基于通量矢量分裂或通量差分裂的迎風格式差分方法，能夠快速捕獲激波。NACA64A010翼型流場采用2×257×129 C型結構網格，如圖2所示。

圖2 NACA64A010翼型流場網格

對于二元翼段，其氣動彈性方程為

(16)

(17)

其中：ρ∞為無窮遠處來流密度；V∞為無窮遠處來流速度；cR為平均氣動力弦長。

(18)

2.2 構建參數化降階模型

為了檢驗參數化降階模型隨參數變化時的精度和可靠性，選擇在馬赫數Ma=0.8～0.9、平均迎角α0=0°～3°范圍內預測氣動力響應。為了從氣動力系統中獲得所關注的幅值和頻率范圍的訓練數據，取無量綱流場時間步長為0.2，FWGN信號的數據長度為2 000，減縮頻率范圍為0.01～0.3，俯仰振幅限制在4°以內，沉浮振幅限制在無量綱長度0.1以內。圖3為選定的兩個自由度上的FWGN信號。在飛行參數范圍內選擇24個樣本點，如圖4所示，并在各樣點處進行CFD計算獲得激勵信號下的氣動力響應。基于樣本點處的氣動力數據，取m=3,n=2,整理得到訓練數據集，利用上述參數化降階建模方法分別構建升力系數CL和力矩系數Cm的降階模型。

圖3 NACA64A010翼型的FWGN激勵信號

圖4 參數空間及采樣點

2.3 SVR-ROM降階模型有效性檢驗

基于上述24個樣本點處的氣動力數據，構建SVR-ROM降階模型，其訓練時間為7.8 h。為了驗證該降階方法在參數空間內的效果，選取3個測試點:Ma=0.81,α0=0.5°;Ma=0.85,α0=2.5°;Ma=0.89,α0=1.5°;分別采用降階模型SVR-ROM和CFD方法預測在FWGN激勵下的氣動力響應。圖5給出了3個測試點處由SVR-ROM和CFD預測的升力系數和力矩系數時域比較。表1給出了3個測試點處SVR-ROM與CFD預測結果之間的相對誤差。計算結果表明，在FWGN激勵下由SVR-ROM預測的氣動力響應與CFD的計算結果吻合很好。

表1 FWGN激勵下SVR-ROM與CFD結果的相對誤差

圖5 FWGN激勵下SVR-ROM和CFD預測的氣動力響應

為了考察SVR-ROM降階模型在其他輸入下的可靠性，分別在3個測試點處使沉浮和俯仰自由度分別做不同頻率的間諧運動，采用降階模型SVR-ROM和CFD方法預測氣動力響應。如圖6所示，在Ma=0.81、α0=0.5°時，降階模型SVR-ROM預測的結果與CFD計算結果吻合較好。表2給出了3個測試點處SVR-ROM與CFD預測結果之間的相對誤差。

圖6 Ma=0.81、α0=0.5°時間諧激勵下氣動力響應

表2 間諧激勵下SVR-ROM與CFD結果的相對誤差

2.4 IL-ROM降階模型有效性檢驗

為了檢驗IL-ROM降階模型隨參數變化時的精度和可靠性，在2.2節參數范圍基礎上進行拓寬，選擇飛行參數范圍為馬赫數0.8～0.96，平均迎角0°～3°，在參數范圍內選擇36個樣本點，如圖7所示。對比圖4發現在馬赫數0.92～0.96中增加12個樣本點，所以在此分別采用LS-SVM和增量學習算法構建降階模型。本節采用與2.2節中相同的激勵信號獲得新增樣本點處的氣動力數據。整理新參數空間內36個樣本點處的氣動力數據，基于LS-SVM構建新的SVR-ROM降階模型。在2.3節SVR-ROM基礎上，補充增加的12個樣本點氣動力數據，基于增量學習算法構建IL-ROM降階模型。

圖7 新參數空間及采樣點

為了檢驗上述兩種降階模型的有效性，同樣選取3個測試點:Ma=0.83,α0=0.5°;Ma=0.89,α0=2.5°;Ma=0.95,α0=1.5°;分別采用SV-ROM、IL-ROM和CFD方法預測在FWGN激勵下的氣動力響應。圖8給出了3個測試點處由以上3種方法預測的升力系數和力矩系數時域比較。表3給出了3個測試點處SVR-ROM和IL-ROM 與CFD預測結果的相對誤差。對比結果表明，在FWGN激勵下由SVR-ROM和IL-ROM預測的氣動力響應誤差相差不大，且都與CFD的計算結果吻合很好。

表3 FWGN激勵下SVR-ROM和IL-ROM與CFD結果的相對誤差

圖8 FWGN激勵下SVR-ROM、IL-ROM和CFD預測的氣動力響應

采用與2.3節相同的方法，在3個測試點處，使沉浮和俯仰自由度分別做不同頻率的間諧運動，采用兩種降階模型和CFD方法預測氣動力響應。如圖9所示，在Ma=0.83、α0=0.5°時，降階模型SVR-ROM和IL-ROM預測的結果與CFD計算結果吻合較好。表4給出了3個測試點處SVR-ROM和IL-ROM與CFD預測結果之間的相對誤差。對比結果表明，在間諧激勵下由SVR-ROM和IL-ROM預測的氣動力響應誤差相差不大，且都與CFD的計算結果吻合的很好。

圖9 Ma=0.83、α0=0.5°時間諧激勵下氣動力響應

表4 間諧激勵下SVR-ROM和IL-ROM與CFD結果的相對誤差

以上算例均在具有8核CPU、3.2 GHz主頻和16 G內存的工作機上運行，基于36個樣本點處的氣動力數據，SVR-ROM和IL-ROM兩種降階模型的訓練時間分別為27.3 h和2.6 h，需要說明的是IL-ROM模型是在2.2節降階模型基礎上訓練得到的。對比結果表明，兩種降階模型的精度差不多，但增量學習算法訓練時間要快得多，建模效率得到很大提高。

2.5 氣動彈性仿真

圖10 時沉浮自由度的氣動彈性響應及 FFT 分析

圖11 時俯仰自由度的氣動彈性響應及 FFT 分析

圖12是基于CFD方法和IL-ROM方法計算機翼顫振邊界的對比，二者結果基本一致。從圖中可看出，顫振邊界表現為“雙臺階”型，分別在Ma∞=0.85,0.9處顫振速度出現兩次較大提升。產生這種現象的原因認為是激波運動導致氣動力對機翼做功降低，從而使顫振速度大幅提高。

圖12 CFD和IL-ROM方法預測NACA64A010翼型的顫振邊界

與CFD/CSD耦合方法相比，基于非定常氣動力降階方法的計算效率具有巨大優勢。本節氣動彈性仿真在同一臺計算機上完成，CFD/CSD耦合方法計算一個流速下氣動彈性響應大約需要37 min，而基于降階方法完成同樣氣動彈性響應預測只需約1.3 s。若是對于三維機翼或全機模型，氣動力降階方法的效率優勢將更加明顯，同時表明基于降階方法的氣動彈性模型適用于氣動伺服彈性分析。

3 結 論

本文提出了一種基于最小二乘支持向量回歸算法的跨聲速非定常氣動力參數化降階模型，并發展了增量學習算法，在增加樣本數據時可對降階模型進行快速二次訓練。該降階方法應用于兩自由度NACA64A010翼型，預測不同來流參數下的氣動力和氣動彈性響應，論證了其有效性，得出以下結論：

1)基于最小二乘支持向量回歸算法構建的參數化降階模型SVR-ROM可預測定義參數空間內隨機激勵下的跨聲速非定常氣動力，與CFD計算結果比較，整體相對誤差不超過4%。

2)拓寬參數空間中馬赫數的范圍，增加訓練數據時，基于增量學習算法構建IL-ROM模型，與SVR-ROM進行重新訓練相比，建模效率得到大幅提高，并且兩種模型誤差接近。

3)與CFD/CSD耦合方法相比，IL-ROM模型能夠準確地預測不同馬赫數的顫振邊界，但計算效率得到大幅提高。未來，該降階模型可以應用到更復雜的多自由模型進行氣動力預測和氣動彈性分析，也可與控制器耦合進行氣動伺服彈性分析。