基于POD和BPNN的流場快速計算方法

賈續毅, 龔春林, 李春娜

(西北工業大學 航天學院, 陜西 西安 710072)

在流場分析和飛行器設計優化等領域中,獲取氣動數據主要有3種途徑:經驗公式或半工程化的經驗方法;計算流體力學(computational fluid dyna-mics,CFD)方法;風洞等物理實驗和試飛的方法。其中,由于CFD相比于前者求解精度更高,相比于后者成本大大降低并且安全,成為了當前獲取高精度氣動數據的主要手段。但是針對復雜外形或復雜流場,CFD計算十分耗時,并且收斂較慢。當需要對大量外形和狀態進行CFD求解時,計算量相當可觀。如何提高CFD的分析效率,加快流場分析的收斂過程,是研究人員一直以來關注的重點之一[1-2]。

提高CFD分析效率的途徑可以歸為兩類:①從數值求解技術上對網格、格式等方面進行改進[3];②通過構建流場的近似模型,對未知狀態流場進行預測。流場數據一方面由于網格單元/節點數目多,具有數據量大的特點，另一方面具有特征冗余的特點。神經網絡(neural network,NN)模型適合于大數據量的建模,但是如果樣本數據存在大量冗余特征,會導致模型學習困難,泛化能力弱等問題[4-5]。

降階模型是一種將高維空間映射到低維空間的方法。采用降階模型,可以有效地去除數據的冗余特征,提高樣本量/特征數比。本征正交分解作為降階模型的一種,在流體力學中得到廣泛應用,例如,構建非結構網格有限元流體流動模型[6],Navier-Stokes方程降階[7]等。在流場分析中,POD可以把高維流場數據降階,映射到低維的正交基模態空間,從而分析出流場的主要特征及其對應的基模態系數。但是,單獨基于POD無法實現流場預測,而代理模型可以建立從設計變量到氣動數據的近似模型,并在氣動優化等領域得到廣泛應用。已有研究表明:將POD與代理模型結合使用,可以對流場進行預測[8-9]。邱亞松等[8]基于POD與RBF(radial basis function)、Kriging代理模型實現了變幾何翼型的定常流場預測;王晨等[9]使用POD-Kriging方法實現了雙圓柱非定常流場的預測。但是,傳統的Kriging、RBF代理模型對于流場預測這種多輸入多輸出問題的建模效率較低。

NN模型可以通過隱藏層、神經元個數的設置來建立多輸入參數到多輸出參數的模型,并且可以保證一定的預測精度和泛化能力,因此在流場分析領域得到廣泛應用[10-11]。Saleem等[12]采用了反向傳播神經網絡和RBF神經網絡,建立了翼型表面壓力系數分布到升力系數、阻力系數和力矩系數的神經網絡,并應用于翼型反設計優化。Zhu等[13]使用RBF神經網絡,利用NACA0012翼型的3個亞聲速流場樣本數據,實現了數據驅動的湍流模型建立。利用NN直接對高維流場數據(數據量與網格單元/節點數量有關)建模,會出現計算內存消耗大,訓練耗時長,建模復雜等問題。

本文結合降階模型與NN模型的特點,提出了一種基于POD和BPNN的流場快速計算方法。其創新點在于,通過引入BPNN,實現對POD基模態系數多輸入多輸出模型的構建。采用流場分區策略,降低POD建模難度,提高建模效率;采用聚類方法篩選建模樣本,降低模型訓練耗時,提高模型泛化能力。該方法被應用于亞/跨聲速下變幾何翼型的定常流場預測,使流場計算效率有了明顯提升。

1 建模方法

1.1 POD模型

POD是一種使降階后各維度上樣本方差最大化的線性降維方法,可以通過求解樣本協方差矩陣的特征值尋求方差最大軸[14]。設原始數據為X(i)(i=1,2,…,r),其中X(i)為n維向量,r為樣本數目。樣本的均值為

(1)

(2)

該協方差矩陣的大小為n×n,其前m階特征值從大到小依次為λ(1),λ(2),…,λ(m),對應的特征向量為q(1),q(2),…,q(m),即POD基。這m個特征向量構成m×n階的樣本線性映射矩陣U=[q(1),q(2),…,q(m)]T,則得到降階后的樣本向量矩陣為

(3)

Y為m×r階矩陣,實現了原始數據維度從n維到m維的降階。需要注意的是,為了使降階后的樣本數據在空間中是線性無關的,應保證降維后的數據特征數小于樣本數,即m

1.2 BPNN模型

BPNN由于具有較強的非線性映射能力和高度自學習和自適應能力,在分類、回歸等問題中得到廣泛應用。一個典型的BPNN包括輸入層、隱藏層和輸出層[15]。本文針對少量流場樣本建模,當訓練樣本相對較少時,淺層神經網絡在訓練以及預測中往往具有更強的穩定性[16],因此本文選取單隱層神經網絡,如圖1所示。

圖1 單隱層BPNN

該神經網絡建立了從m維輸入到n維輸出的映射,其中輸入為x=[x1,x2,…,xm]T, 隱藏層神經元個數為q,輸出為y=[y1,y2,…,yn]T。ω1,ω2為權值矩陣,其大小分別為q×m、n×q。

BPNN通過反向傳播誤差來更新連接權值,從而完成模型的訓練,具體訓練步驟如下:

1) 初始化權值,一般通過設定權值范圍隨機給出。

2) 計算隱藏層和輸出層的輸出值。

3) 計算各層反向誤差。

4) 判斷是否滿足訓練要求,若滿足,則表明模型訓練完畢,輸出模型;若不滿足,更新權值,重復第2)～3)步。該網絡的訓練流程如圖2所示。

圖2 BPNN訓練流程

2 流場快速計算方法

2.1 POD & BPNN

結合POD與BPNN模型,本文提出了一種流場快速計算方法,其流程如圖3所示。具體流程分為4步。

圖3 基于POD & BPNN模型的流場計算流程圖

1) 樣本流場生成

2) 對樣本流場數據降階

以壓強場數據Xp為例,其壓強場均值為

(4)

標準化后的樣本數據為

(5)

(6)

式中,每個特征向量表示壓強場的一個基模態,進而得到降階后的m×r階矩陣Yp

(7)

當取合適值時,可以保證

(8)

由于BPNN對于多輸出問題仍然適用,因此無需像傳統的Kriging和RBF代理模型那樣對m個基模態系數輸出建立m個映射模型。這m階基模態系數之間的量級差距較大,如果m取值較大,直接建立的BPNN難以對所有的基模態系數做出精準預測。因此,可以根據m的具體數值,例如每5～6個輸出建立一個BPNN模型,即可實現較好的模型擬合效果。

4) 預測流場

對于一個新的任意輸入參數Cs,通過BPNN模型可以預測出基模態系數Ys,然后對Ys進行POD重構,得到預測壓強場數據

(9)

2.2 建模問題

在飛行器設計和優化中,需要對大量計算狀態或者幾何外形進行CFD分析。針對翼型優化過程中幾何外形不斷變化,需要在優化迭代過程中進行大量CFD分析的問題,研究基于POD & BPNN的流場快速預測方法的效率和精度。

在對流場建模過程中,輸入參數取樣對應的是在變幾何翼型的翼型庫中選擇翼型樣本。本文的翼型庫通過Hicks-Henne參數化方法獲得。以典型的NACA、RAE系列翼型為基準翼型,在它的上下翼面各施加5個擾動,擾動的位置分別為0.10c,0.25c,0.45c,0.65c,0.80c,c為弦長,本文取1。10個擾動的大小取值范圍為[-0.005,0.005],是無量綱化的系數。因此,這10個無量綱化系數可以確定出一個新翼型的幾何形狀。本文通過拉丁超立方采樣方式在亞、跨聲速下分別獲取200和500組外形的系數,從而得到新翼型樣本庫。

以跨聲速為例,對網格收斂性和流場求解精度進行驗證。對RAE2822進行O型網格劃分,壁面第一層網格高度為1×10-5c,計算工況為:Ma=0.729,α=2.31°,Re=6.5×106,網格數目及相應的氣動計算結果如表1所示。δL,δD分別表示與網格Ⅳ的升阻力系數的相對誤差。可見網格Ⅲ(見圖4)可以滿足收斂性的要求。因此,本文算例采用網格Ⅲ設置生成其他樣本翼型網格。

圖4 RAE2822網格局部圖

表1 網格收斂性分析

2.3 分區策略

流場數據屬于高維數據。要減少建模數據的維度:①通過網格收斂性驗證可以在保證計算收斂的前提下有效減少網格數目;②可以通過基于輸入參數和局部流場信息對流場數據進行分區,實現區域內的降維[13]。在本文工作中,考慮到流場中不同區域的特性沿著附面層法向連續變化,將流場數據分為近壁區數據(包含全部附面層)、中間場數據和遠場區數據,每個區對應網格層數均為40(跨聲速下的壓強場數據分布見圖5)。這樣,流場的某一個物理量,需要利用POD&BPNN建立3個模型,來對流場信息進行預測。數據維度降低,但是模型數目增加,因此需要對模型的訓練時間和模型的泛化能力進行評估。

圖5 壓強場分區后的數據分布

首先對生成的樣本翼型的流場數據進行POD分析,可以得到基模態流場及對應的特征值。定義前q階基模態對應的特征值之和與所有特征值之和的比值為前q階基模態的能量占比,能量占比表征了前q階基模態所包含的流場信息量。圖6為亞聲速下分區后近壁區、中間場和遠場區以及未分區情況下壓強、馬赫數、密度3個物理量對應場數據的基模態的能量占比。

圖6 亞聲速分區與未分區情況下的基模態能量占比

從數值大小上看,亞聲速下前10階基模態的能量占比超過99.6%,同樣計算得到跨聲速下前30階基模態的能量占比超過99.7%,說明流場特征主要集中在前數階POD基模態中。

分別對亞聲速和跨聲速的樣本隨機抽取90%作為訓練樣本進行POD降階,將剩余10%的樣本作為測試樣本,進行POD重構。為了定量分析重構效果,定義測試樣本的重構(或預測)誤差為

(10)

對于T個測試樣本的整體誤差定義如下

(11)

圖7為亞/跨聲速下,選擇不同基模態個數對測試樣本壓強場進行重構的誤差。可以發現,亞聲速狀態的流場特征主要集中在前10階,與圖6的能量占比信息相吻合;跨聲速狀態的流場特征主要集中在前30階。此后,基模態數目的繼續增加對流場重構精度的提升不再有顯著影響。從是否分區上看,分區對亞聲速下模型精度的影響不大;對跨聲速情況影響會大一點,但是并不顯著,分區后最大誤差會變小,而平均誤差變化不大。

圖7 亞/跨聲速的壓強場重構誤差

POD建模的主要計算耗時來源于協方差矩陣特征值的求解。矩陣的特征值計算僅從浮點運算次數而言,其復雜度為O(n3),其中n為協方差矩陣的階數,即流場數據的維度。因此,矩陣階數的減少可以大幅降低計算耗時[17]。

本文所涉及的計算均在Inteli5-8500CPU、16G運行內存的PC上進行。以跨聲速算例為例,圖8為未分區和三分區下對單個物理量場數據的POD降階計算時間。其中,三分區的計算時間為建立3個分區POD模型的總時間。可見,對于POD模型,3個分區建模的效率約為未分區建模的3倍。此外,由于協方差矩陣的儲存空間與其階數n的平方成正比,單個分區建模的內存消耗僅為未分區建模的1/9。

因此,針對本文算例使用分區策略可以基本保持未分區模型的精度,有效縮短POD建模耗時,降低內存占用量。通過POD重構的誤差分析,可以確定亞、跨聲速狀態下建模的基模態個數分別為10,30以上。本文亞/跨聲速下選取的基模態個數為13和35。

2.4 聚類取樣策略

訓練樣本的數目以及選取方式對模型的訓練有較大影響。首先,訓練樣本應保證在設計空間中均勻分布,這樣模型才能較好地捕捉到空間中的全部特征。其次,訓練樣本量應適中,樣本量過少會導致訓練信息不足,容易發生欠擬合情況;樣本量過多會增加訓練難度,使模型訓練緩慢。因此,對于不確定合適訓練樣本數目的問題,可以首先生成一定數目的均勻分布的訓練樣本,再通過樣本量壓縮來減少訓練樣本數目。

樣本量壓縮的方式有多種,主要采用以下3種方式:①簡單刪減,即直接刪掉一部分樣本量;②逐次刪減,即逐次增加刪減個數或逐次減少刪減個數;③聚類分析,即將相似的樣本歸為一類,每類只保留一個樣本。考慮壓縮樣本量的同時不能降低模型精度,因此本文采用了聚類分析。

在POD&BPNN模型中,輸入參數是通過拉丁超立方采樣獲得的,其在空間中已經處于一種隨機均勻分布的狀態,因此不適用于聚類。而樣本流場通過POD分析得到的基模態系數,針對不同的流場特征會存在差異,相似的流場其基模態系數也相似,流場差異大基模態系數差異也大。因此,對基模態系數進行聚類分析,可以將流場信息相似的樣本歸為相同類。這樣在每個類中只取一個或少量樣本作為訓練樣本,可以在保持模型精度的前提下,實現訓練樣本量的壓縮。

下面對聚類取樣策略的效果進行分析。亞/跨聲速下聚類個數設置的范圍分別為[20,190]和[20,490],亞/跨聲速下訓練所得模型對10個測試樣本的預測誤差隨聚類個數的變化曲線如圖9所示。

圖9 聚類數目對模型預測誤差的影響分析

在亞聲速算例中,當聚類數目超過90時,預測誤差不再有明顯變化,說明在亞聲速狀態下通過聚類取樣獲得的90個訓練樣本可以基本表征出全部樣本的特征,從而有效地避免重復樣本對訓練時間的消耗。同樣,在跨聲速算例中,當聚類數目超過200時,預測誤差不再有明顯變化,說明在跨聲速狀態下通過聚類取樣獲得的200個訓練樣本可以基本表征出全部樣本的特征。因此在本文的亞聲速算例中,使用聚類取樣得到的90個樣本作為訓練集,跨聲速算例中,使用聚類取樣得到的200個樣本作為訓練集。

3 案例分析

3.1 亞聲速變幾何翼型

本文亞聲速算例計算工況為:Ma=0.35,α=3°,Re=6.5×106,基準翼型為NACA0012翼型。根據2.3與2.4節的誤差分析,選擇分區建模的POD&BPNN模型的基模態個數為13,聚類個數為90。因此,訓練樣本個數為90,隨機抽取10個測試樣本。測試樣本的升阻力系數的預測值與真實值的對比如圖10所示。10個測試樣本的升力系數的平均相對預測誤差為0.14%,阻力系數的平均相對預測誤差為0.37%。

圖10 亞聲速下測試樣本的氣動性能預測效果

取圖10中第4個樣本,對比預測和真實的翼面壓力系數、翼面摩擦力系數分布、等壓線分布，如圖11～13所示。可見,亞聲速下翼面壓力系數和摩擦力系數的預測值與真實值無明顯差異,流場等壓線重合的很好,說明模型能夠準確預測流場信息。

圖11 亞聲速下翼面壓力系數分布對比 圖12 亞聲速下翼面摩擦力系數分布對比 圖13 亞聲速下等壓線分布對比

表2中給出了亞聲速狀態下在相同的訓練和測試樣本條件下,是否使用分區策略的模型訓練效率與預測精度指標。其中,POD單次建模耗時是指對一個物理量(如壓強場)的單個區的數據完成POD降階的耗時;BPNN單次建模耗是指對一個物理場的單個區完成所需基模態個數的訓練耗時。可以看出,相比未分區,使用分區建模策略的模型總耗時有所減少,模型的精度略有提升。這種分區建模效率上的優勢對于3D構型的大規模網格更為突出。使用分區建模策略在亞聲速算例中是可行的。

表2 亞聲速下的模型訓練效率與預測精度指標對比

3.2 跨聲速變幾何翼型

本文中，跨聲速算例的計算工況為:Ma=0.75,α=2.79°,Re=5.7×106。根據2.3與2.4節的誤差分析,POD&BPNN建模的基模態個數為35,聚類個數為200。因此,訓練樣本個數為200,隨機抽取10個測試樣本。測試樣本的升阻力系數的預測值與真實值的對比如圖14所示。10個測試樣本的升力系數的平均相對預測誤差為0.40%,阻力系數的平均相對預測誤差為1.31%。以升阻力系數預測效果較差的第4個測試樣本為例,對比預測和真實的翼面壓力系數、摩擦力系數分布、等壓線分布，如圖15～17所示。可見,跨聲速工況中,下翼面壓力系數的預測值與真實值重合很好,上翼面的前緣與激波位置的預測與真實值稍有差異。翼面摩擦力系數在前緣和激波處的預測略有偏差。對于等壓線,激波位置預測很好,但是激波前后的流場預測精度略有下降。

圖14 跨聲速下測試樣本的氣動性能預測效果

圖15 跨聲速下翼面壓力系數分布對比 圖16 跨聲速下翼面摩擦力系數分布對比 圖17 跨聲速下等壓線分布對比

本文模型對激波誘導分離泡現象可以實現較好的預測,圖18為第4個測試樣本的真實和預測的激波誘導分離泡,通過對比發現,預測的激波誘導分離泡的位置和形態與真實情況相吻合。

圖18 激波誘導分離泡的真實與預測對比

表3中給出了跨聲速狀態下在相同的訓練和測試樣本條件下,是否使用分區策略的模型訓練效率與預測精度指標。可以看出,相比未分區,使用分區建模策略的模型總耗時有所減少,模型的最大相對誤差和平均相對誤差均值有所降低。同樣,這種分區建模效率上的優勢會在3D構型中體現出來,使用分區建模策略在跨聲速算例中是可行的。

表3 跨聲速下的模型訓練效率與預測精度指標對比

4 結 論

1) 本文提出了基于POD&BPNN模型的流場快速預測方法,相比于傳統的代理模型,BPNN可以更好地實現對多個基模態系數的同時建模。但當幾何外形的網格及其拓撲結構發生變化時,方法會失效。

2) 在POD建模中引入分區策略,可以在保證模型精度的前提下有效地縮短POD建模的時間。當流場數據量很大時(例如,百萬級的網格),分區策略的優勢會更加明顯。

3) 針對亞/跨聲速算例,通過對基模態系數進行聚類分析來確定訓練樣本集,可以有效地降低訓練樣本量,提高模型訓練效率。

4) 亞聲速和跨聲速算例表明,基于POD&BPNN模型預測流場的升阻力系數,其預測相對誤差不超過1.4%。對于跨聲速流場,所建模型能很好地預測激波位置和激波誘導分離泡等復雜流動特征。