基于POD-RBFN降階模型的串列葉柵流場預測

尚珣，劉漢儒，杜亦璨，胡之頡

（1.西北工業大學 太倉長三角研究院，太倉 215400）（2.西北工業大學 動力與能源學院，西安 710129）

0 引言

隨著現代航空發動機對壓縮系統載荷、性能、穩定裕度提出的要求越來越高，傳統的單葉片設計已經逐漸無法滿足現代壓氣機的需求。串列葉片作為一種被動的流動控制技術，因其結構簡單、便于工程應用等特點而得到廣泛關注。在研究串列葉型優化設計時，通常采用數值模擬的方法，但對于串列葉柵這種流場結構較復雜或存在多個設計變量的優化問題，傳統CFD所需要的計算成本往往是不可接受的。為了快速確定氣動參數的函數響應，提高優化設計效率，降階預測模型應運而生。

本征正交分解（Proper Orthogonal Decomposition，簡稱POD）方法是一種高效的數據降維分解方法，可以提取出數據的主導模態，因此可以用于流場分析及降階預測模型的構建。例如，Duan Y等對Gappy POD降階方法進行了研究，將其用 于 葉 型 的 優 化 設 計；Gong Helin等對Gappy POD降階模型進行了進一步發展，提高了方法的穩定性，并將其應用于功率分布的預測。然而，Gappy POD是通過線性響應POD基函數來實現數據填充預測的，這雖然使其具有簡單高效、可以同時響應多個輸出參數等優點，但其對于強非線性響應問題是無效的。

近年來，神經網絡方法逐漸得到了流體力學領域研究者們的關注。例如，Zhu Linyang等將徑向基函數神經網絡（Radial Basis Function Network，簡稱RBFN）用于湍流模型的構建，實現了渦黏的預測；Sun Gang等使用翼型數據庫訓練人工神經網絡，實現了翼型/機翼的快速反設計。在流體力學研究中，流場信息在流動機理分析時是不可或缺的?？紤]到流場信息的維度較高，將神經網絡方法直接用于復雜流場信息的獲取會產生維度災難問題，因此需要先降低流場的維度，POD方法就是一種不錯的選擇，且已經在葉輪機械內流領域得到了應用。例如，H.Kato等基于POD方法構建代理模型，并將其應用于渦輪幾何外形優化；李春娜等采用基于POD-BPNN的代理模型實現了翼型的快速反設計。綜上所述，現有研究大多是通過幾何與氣動參數的一一對應關系，實現變幾何的性能預測或變性能參數的幾何反設計，鮮有研究關注非定常流場的演化過程。

POD流場降階需要基于一定的樣本流場來進行，且原始樣本的選擇直接影響著流場的降階質量，進而影響后續降階模型的預測精度。目前各學科常用的靜態抽樣方法包括隨機抽樣、均勻抽樣、拉丁超立方抽樣等，但靜態抽樣方法只考慮了樣本在采樣區間中的空間分布，因此當采樣區間中的樣本質量及敏感性程度分布不均勻時，靜態抽樣方法就無法取得令人滿意的結果。近年來，基于函數響應偏差的自適應抽樣方法已得到迅速發展，該方法在進行采樣時考慮了樣本對于需求的敏感性，通過函數響應偏差確定敏感空間，因此能夠有效提高抽樣效率與樣本質量。

本文以串列葉柵為研究對象，結合POD方法的流場降階能力和徑向基函數神經網絡可以響應強非線性問題的特性，構造POD-RBFN混合降階模型，以實現串列葉柵非定常流場的快速預測；而后將基于函數響應偏差的自適應抽樣方法應用于降階模型初始樣本的獲取，以提高降階模型的采樣效率。

1 數值方法

1.1 網格及計算設置

本文選取J.Eckel等設計的高載荷壓氣機串列葉柵為研究對象，葉柵幾何及主要參數分別如表1和圖1所示。該串列葉柵前后兩個葉片弦長均為31.5 mm，整體弦長62.4 mm，設計工況氣流轉折角為50°，前后兩個葉片的設計載荷相等。

圖1 葉柵幾何示意圖Fig.1 Geometric description of tandem cascade

表1 串列葉柵幾何參數Table 1 Geometric parameters of tandem cascade

采用AutoGrid 5軟件生成計算網格，前、后葉片網格均采用H-O-H結構拓撲。葉柵進口段計算域延伸2倍弦長，出口段延伸2.5倍弦長。壁面附近網格進行加密處理，以保證近壁面第一層網格y值小于1。最終計算網格總數為39 896，此時網格正交性、延展比、長寬比均已滿足計算要求。葉柵通道及局部網格示意圖如圖2所示。

圖2 網格示意圖Fig.2 Schematic diagram of single-channel grid

為了對網格獨立性進行驗證，本文對三套網格（Mesh1，Mesh2，Mesh3）進行對比，網格數分別為22 480，39 896，79 906。0°攻角、來流馬赫數0.6工況下的葉表靜壓升系數（C）分布如圖3所示，橫軸z為垂直于柵距方向的距離。C的定義為C=（P-P）/（P-P），其 中P為 靜 壓，P為 進 口 靜壓，P

圖3 葉表靜壓升系數分布圖Fig.3 Tandem cascade CP distribution

*為進口總壓。

從圖3可以看出：隨著網格數量的增大，Mesh1的靜壓升系數整體略微偏大，而Mesh2與Mesh3的靜壓升系數分布基本一致。綜合考慮計算精度與計算效率，本文采用Mesh2進行后續研究。

采用CFX-Solver進行流場數值求解。進口邊界條件給定總溫為288.15 K，總壓為101 325 Pa。葉片表面采用無滑移邊界條件和絕熱壁面條件，上、下兩個邊界為周期性平移邊界。通過調節進口軸向和節距方向氣流速度分量來實現不同的氣流攻角，調節出口平均靜壓保證進口馬赫數為0.6。采用SST兩方程湍流模型，以實現對流場細節較為準確地捕捉。非定常計算采用URANS方法，物理時間步長為2×10s，每個時間步進行10次內迭代。將定常計算結果作為非定常計算的初場以節省計算時間。

1.2 POD及RBFN理論

早期POD的奇異值分解需要定義域中的所有解集，這使得其在特征值的求解及系統穩定性方面均表現欠佳，從而難以在工程中應用。直到快照POD方法的提出，改善了計算量及穩定性方面的問題，POD方法才逐漸被應用于各個領域?？煺誔OD的基本原理如下：

首先根據一組解集得到快照矩陣U=[u，u，…，u]，矩陣中每一個元素均為一個流場解。然后對矩陣進行去中心化處理：

POD方法的核心工作就是尋找一組最優正交基Φ=[φ，φ，…，φ]，以及與之對應的基系數矩陣A=[a，a，…，a]，使得其中的元素滿足：

快照脈動矩陣的最優正交基可通過求解如下最優值問題得到：

使用奇異值分解（SVD）求解該最優值問題。獲取最優正交基Φ后，將各個流場快照分別映射到Φ上，得到與各個初始流場解對應的基系數a。

按照基模態“能量”的定義，對基模態進行排序。前幾階模態包含了流場中的大部分“能量”，因此使用較少的模態階數即可實現原始流場的降階重構。某一階基模態的“能量”定義為

式中：λ為自相關矩陣R（R=VV）的第i個特征值。

求解得到最優正交基Φ之后，原始流場可以根據系數矩陣A及式（2）精確重構。因此，如果可以得到某非樣本流場解對應的系數向量，就能根據已知的正交基Φ對該流場進行預測。

非樣本流場POD系數的響應通常采用線性方法，但考慮到流場時空演化中可能存在的非線性規律，本文采用徑向基函數神經網絡來響應POD系數。徑向基函數神經網絡（RBFN）是由D.S.Broomhead等提出的一種收斂快、結構簡單，且對非線性系統具有一致逼近性的神經網絡。RBFN是經典的三層前向網絡，第一層為輸入層，第二層為隱藏層，第三層為輸出層。隱藏層使用徑向基函數作為激活函數的神經元。

至此，POD-RBFN混合模型構建完成，該降階模型可以根據一組同類型的流場解集（可以是同一幾何在不同工況下的流場解，也可以是同一工況下不同葉型幾何參數的流場解）獲取流場基模態，并對目標流場（非樣本流場）的POD系數進行非線性響應，從而實現目標流場的預測。本文的混合預測模型是通過編寫計算機程序實現的。

1.3 自適應抽樣方法

在降階模型重構過程中，原始樣本流場的質量直接影響POD基的生成，進而影響重構精度。為了提高獲取的POD基的穩定性，考慮采用自適應抽樣方法獲取原始樣本流場。近年來，基于函數響應偏差的自適應抽樣方法已得到迅速發展，故本文采用四叉樹法進行原始樣本抽取。四叉樹法自適應抽樣的基本原理為：首先將樣本空間中的邊界角點、空間中心點、子空間中心點作為初始樣本集；其次通過計算子空間中心點的函數響應偏差來確定目標子空間；然后在目標子空間中抽取其下一代子空間的中心點，并將其增加到初始樣本集中；最后在新樣本集中重新確定目標子空間。重復該步驟，直到最終所有采樣點的響應偏差均滿足要求或總樣本數目滿足要求。

這一過程中的函數響應偏差由式（4）計算得到：

式中：E（u）為第i個流場快照u的響應偏差；D（u）為第i個流場樣本u與其他樣本之間的相關性強弱。

式中：n為樣本總數，代表樣本集合中第i個樣本流場與第j個樣本流場之間的線性相關性系數。

式（5）的含義為，樣本流場u與其他樣本流場的不相關性（1-ρ）的均值，相關性越強則D（u）越小。

流場快照u的響應偏差E（u）是由“棄一法”（Leave-one-out，簡稱LOO）交叉驗證得到的。該方法的基本原理為：首先從初始快照集中去除掉樣本快照u，然后基于剩余樣本對缺失樣本u進行預測，預測得到流場的整體誤差即為該樣本流場快照的響應偏差。

綜上可知，如果某一個樣本與其余樣本的相關性越強，或者該樣本的響應偏差越小，則說明該樣本對于生成POD基及流場重構的貢獻不大，而計算得到的該樣本的函數響應偏差也就越小。因此可以根據函數響應偏差的大小來確定目標子空間。編寫計算機程序并與上文得到的降階模型進行適配，最終完成采用自適應抽樣的POD-RBFN降階預測模型的構建。

2 結果與討論

2.1 CFD計算結果

首先通過定常數值模擬對該串列葉柵的氣動特性進行研究。分別在進口一倍弦長和出口0.5倍弦長處設立計算面，計算得到的葉柵壓比、損失攻角特性曲線如圖4所示，左邊縱軸P為葉柵的靜壓比，右邊縱軸ω為葉柵損失系數，圖中灰色圓點為實驗所測得的葉柵損失。從圖4可以看出：隨著進口氣流攻角的增大，葉柵靜壓比呈現出先增大，后急劇減小的趨勢；而葉柵損失系數ω的變化趨勢則大致與P變化趨勢相反。這是因為隨著氣流攻角的增大，葉盆分離流動逐漸減弱，故此時靜壓比逐漸增大，損失逐漸減??；而當氣流攻角進一步增大，此時會在葉片背面發生氣流分離，使得流場惡化，靜壓比急劇減小，且損失急劇增大。

圖4 葉柵攻角特性Fig.4 Incidence characteristic of tandem cascade

綜合考慮不同攻角下的兩種參數值，可知該串列葉柵在攻角［-10°，4°］范圍內，當攻角為-2°時的氣動性能最佳，而在攻角為4°時性能最差。為了突出流場中葉背分離流動的不穩定流動特征，對大攻角分離工況下的流動機理進行研究，選擇2°攻角工況開展后續研究。

非定常計算得到的葉片附近流場馬赫數云圖及流線分布如圖5所示，可以看出：在較大的氣流攻角下，葉背的大范圍分離流動主要發生在前葉片。這是由于一方面，與常規單列葉柵相比，串列葉柵單個前葉片的弦長較短且氣流轉折角較大，使得前葉片更易發生分離；另一方面，由于前、后葉片間隙射流的存在，前葉片尾緣分離產生的低能流體不易與后葉片氣流發生摻混作用，從而使得前葉片的分離更為嚴重。與前葉片相反，由于串列葉柵縫隙射流的吹除作用，后葉片一直處于一個較為穩定的來流條件，故未發生明顯的分離流動。

圖5 流場馬赫數及流線分布Fig.5 Mach number and streamline

2.2 POD流場分解

非定常數值計算得到的通常是一組離散時刻的流場解，而兩個時刻之間的流場時空演化規律可能是非線性的，因此若想獲得某非結果時刻的流場解時，采用簡單的線性插值可能是不準確的。通過對串列葉柵非定常流場解的預測，來驗證構造出的POD-RBFN混合模型的準確性。

當非定常計算整體殘差穩定在0.000 01以下，且進、出口流量隨時間小幅周期性波動時，認為此時計算已經收斂，保存100個連續時間步的非定常流場解作為原始樣本集合，相鄰時間步的時間間隔為0.000 02 s。本文使用流場渦量分布構造流場快照矩陣，對流場快照矩陣進行POD分解，得到相應的POD基函數和與之對應的POD系數。

通過求解自相關矩陣R的特征值λ，將特征值按照大小排列。這樣就可以根據式（3）得到各階POD基模態的“能量”占比，進而可以確定前多少階模態可以包含原始流場中的絕大部分信息。計算得到的前10階POD基模態所占流場的能量份額，以及前n階模態總能量份額的大小如表2所示。前30階各階POD基模態所占能量份額的變化趨勢如圖6所示，橫坐標為基模態的階數，縱坐標為對數形式的各階模態所占能量的比例E。

圖6 各階模態能量分布Fig.6 Energy distribution of each POD mode

表2 各階模態所占能量比例Table 2 The energy proportion of each mode

從表2及圖6可以看出：從8階開始，后續各階模態能量占比均小于2%；前10階模態占總模態能量的比例已經高達93.96%，因此由前10階基模態重構得到的流場已經可以包含流場中的絕大多數信息。

限于篇幅，本文只給出第1、3、5、7階POD基模態分布云圖，如圖7所示，可以看出：POD基模態渦核結構沿流向正負交替出現，前葉片葉背分離渦和尾緣脫落渦結構的波動均明顯強于后葉片，占據流場中的絕大部分渦系能量，說明前葉片的流動不穩定性高于后葉片；此外，串列葉柵間隙射流在一定程度上阻礙了前葉片尾緣脫落渦與后葉片葉背分離渦之間的干涉作用，從而避免了后葉片葉背分離渦的增強，使得后葉片即使在較大的來流攻角下也能保持較好的流動狀態；隨著階數增大，高階模態云圖中相鄰渦核結構的距離縮小，表明高階模態可能捕捉到了串列葉柵流場中更小尺度的渦系結構，渦核之間的相互作用與影響也更為復雜。

圖7 各模態無量綱化能量分布Fig.7 Dimensionless energy distribution of each mode

t=1時間步原始流場渦量云圖及使用低階模態重構得到的渦量云圖如圖8所示，可以看出：使用前4階POD基模態重構的流場已經可以反映原始流場的絕大多數流動結構，前8階基模態重構流場與原始流場僅存在細微差別，隨著模態階數進一步增大，重構流場對原始流場細微結構的表達也更加清楚。這一結果與根據模態能量占比得出的結論一致。

圖8 原始流場及重構流場渦量分布圖Fig.8 Original& reconstructed vorticity distribution

2.3 流場預測結果分析

從上述原始樣本集（100個連續時間步流場）中隨機抽取60個樣本流場快照，用于獲取POD基函數及相應的系數。而后將60個樣本流場的時間步作為輸入，將POD模態系數作為輸出，構造RBF神經網絡降階模型，從而實現對非樣本流場POD模態系數的預測。

使用POD-RBFN混合模型對非樣本流場進行預測，定義流場整體重構誤差為

式中：下標m為重構的目標流場時間步數；N為流場馬赫數分布云圖中的節點數；Ma為該時刻重構流場中第i個流場節點的當地馬赫數值；Ma為該時刻原始流場第i個節點的當地馬赫數值。

分別采用不同的模態階數對某時刻流場（非樣本流場）進行預測，得到整體重構誤差隨模態階數的變化關系，如圖9所示，可以看出：當模態階數較少時，重構誤差隨模態階數的增加急劇減小，但當模態階數大于22后，重構誤差隨模態階數增加而變化的幅值很小。這是由于前20階基模態已經包含了流場98%以上的能量，更高階的基模態所占據的流場能量份額極少，對流場重構的影響也很小。

圖9 不同階模態重構流場整體誤差Fig.9 Overall reconstruction error with different modes

使用前10階模態對某時刻流場（非樣本流場）進行預測的結果如圖10所示，重構流場的誤差分布云圖如圖11所示，可以看出：使用前10階模態預測的馬赫數流場與原始流場僅存在細微差異；重構流場的誤差主要分布在前、后葉片葉背分離渦及尾緣脫落渦渦核位置處，這是由于渦系結構的相互干涉、演變，使得該區域流動較為復雜多變，進而影響了流場的預測精度。

圖10 原始流場及重構流場馬赫數云圖Fig.10 Mach number distribution of original and reconstructed flow field

圖11 重構流場誤差分布云圖Fig.11 Reconstruction error distribution

2.4 自適應抽樣結果分析

本文構造的POD-RBFN混合模型需要基于一定的樣本數才能進行流場預測，而函數響應偏差的計算需要進行流場預測，故需要先從原始樣本空間的100個流場快照中等距抽取11個樣本作為自適應抽樣的初始樣本，然后采用基于函數響應偏差的自適應抽樣方法進行抽樣。抽樣對象為一維時間步序列，故每次抽樣后會得到兩個新的樣本點，總共進行20次自適應抽樣。

分別使用自適應抽樣和隨機抽樣得到的樣本點對隨機25個時間步的流場進行預測，并求取重構誤差的均值，平均重構誤差隨樣本數增加而變化的曲線如圖12所示，可以看出：經過自適應抽樣之后，相同樣本數下的重構精度得以提高，達到同一精度所需樣本數量明顯降低，比如圖中采用自適應抽樣之后，只需要45個樣本就能達到隨機抽樣60個樣本的重構精度，樣本數量減少了25%。

圖12 重構誤差與樣本數的關系曲線Fig.12 Relation curve between reconstruction error and sample number

分別使用自適應抽樣和隨機抽樣得到45個樣本，對某時間步（非樣本時間步）的流場馬赫數分布進行預測，兩種抽樣方法進行重構的誤差分布云圖如圖13所示，可以看出：自適應抽樣方法明顯降低了渦系結構發生、發展區域的流場重構誤差，尤其是葉片下游尾跡處的流場重構誤差，從而提高了流場的整體重構精度。

圖13 兩種抽樣方法的重構誤差分布云圖Fig.13 Comparison of reconstruction error distribution

3 結論

（1）低階POD模態包含了流場中的絕大部分能量；串列葉柵間隙射流的存在一定程度上阻礙了前葉片尾緣脫落渦與后葉片葉背分離渦之間的干涉作用，改善了后葉片的流動狀態；高階POD模態捕捉到了串列葉柵流場中更小尺度的復雜渦系結構。

（2）POD-RBFN混合模型能夠有效捕捉流場的主要模態特征，并通過低階重構的方式對非樣本流場進行快速預測，預測誤差主要分布在前、后葉片葉背分離渦及尾緣脫落渦渦核位置處。

（3）采用基于函數響應偏差的自適應抽樣方法能夠有效提高降階模型的采樣效率，與靜態隨機采樣相比，達到同樣重構精度所需的樣本數降低了25%左右。