賈續(xù)毅，李春娜，常琦，季穩(wěn)

（1.西北工業(yè)大學(xué)陜西省空天飛行器設(shè)計(jì)重點(diǎn)實(shí)驗(yàn)室，西安 710072）（2.火箭軍裝備部 裝備項(xiàng)目管理中心，北京 100085）

0 引言

隨著計(jì)算機(jī)技術(shù)的發(fā)展，計(jì)算流體力學(xué)（Computational Fluid Dynamics，簡稱CFD）成為獲取高可信度氣動數(shù)據(jù)的主要手段并被大量應(yīng)用于氣動優(yōu)化設(shè)計(jì)。梯度優(yōu)化算法和啟發(fā)式優(yōu)化算法是當(dāng)前常用的氣動優(yōu)化設(shè)計(jì)算法。相比于梯度優(yōu)化算法，差分進(jìn)化（Differential Evolution，簡稱DE）算法、粒子群算法等啟發(fā)式優(yōu)化算法具有更好的全局搜索能力，從而在氣動優(yōu)化中得到廣泛使用。但是，該類算法在氣動優(yōu)化過程中需要進(jìn)行成千上萬次CFD仿真，導(dǎo)致優(yōu)化效率低下。為了兼顧全局搜索能力和優(yōu)化效率，代理優(yōu)化（Surrogatebased Optimization，簡稱SBO）應(yīng)運(yùn)而生。

SBO具有設(shè)計(jì)效率高、全局優(yōu)化能力好、魯棒性高等特點(diǎn)，在具有高維設(shè)計(jì)變量、多目標(biāo)等工程優(yōu)化問題中具有良好的應(yīng)用前景。目前已發(fā)展了Kriging模型、徑向基函數(shù)（Radical Basis Function，簡 稱RBF）、多 項(xiàng) 式 混 沌 展 開、神 經(jīng) 網(wǎng)絡(luò)等多種代理模型建模方法。在氣動SBO中，代理模型建立了從氣動外形參數(shù)到氣動力的映射關(guān)系。SBO本質(zhì)上降低了CFD計(jì)算次數(shù)但沒有改善單次CFD計(jì)算時間，且沒有有效使用CFD的流場信息。

機(jī)器學(xué)習(xí)技術(shù)在氣動力建模、流場降階、流動特征提取等流體力學(xué)領(lǐng)域得到大量應(yīng)用和發(fā)展，S.A.Renganathan等使用深度神經(jīng)網(wǎng)絡(luò)和高斯過程實(shí)現(xiàn)了氣動外形的優(yōu)化設(shè)計(jì)，但該方法仍為梯度優(yōu)化；Sun Z W等、Huang D M等使用具有特征提取能力的降階模型如本征正交分解（Proper Orthogonal Decomposition，簡稱POD）和具有高效建模能力的神經(jīng)網(wǎng)絡(luò)模型建立了從氣動外形參數(shù)到流場信息的映射關(guān)系，實(shí)現(xiàn)了流場預(yù)測。常用的神經(jīng)網(wǎng)絡(luò)包括反向傳播神經(jīng)網(wǎng)絡(luò)（Back Propagation based Neural Network，簡 稱BPNN）、徑 向 基 函 數(shù) 神 經(jīng) 網(wǎng) 絡(luò)等，Zhu L等使用徑向基函數(shù)神經(jīng)網(wǎng)絡(luò)構(gòu)建的渦黏預(yù)測模型代替了CFD求解中的湍流模型，通過與N-S方程耦合求解降低了單次CFD計(jì)算時間，但是該方法需要對CFD求解器進(jìn)行相應(yīng)修改，對于優(yōu)化設(shè)計(jì)人員，鮮見成熟的、已嵌入湍流建模的求解器可供使用。

本文以氣動優(yōu)化為研究對象，首先，結(jié)合SBO和機(jī)器學(xué)習(xí)流場預(yù)測模型的特點(diǎn)，提出一種基于本征正交分解—反向傳播神經(jīng)網(wǎng)絡(luò)（POD-BPNN）模型的熱啟動策略，并將其嵌入SBO；然后將該策略應(yīng)用于跨聲速翼型減阻優(yōu)化，并對優(yōu)化結(jié)果進(jìn)行分析；最后通過案例對優(yōu)化效率、POD-BPNN的流場預(yù)測性能、熱啟動策略的效果進(jìn)行分析。

1 優(yōu)化策略

1.1 優(yōu)化總流程

氣動優(yōu)化中所涉及的單目標(biāo)優(yōu)化問題可以表示為

針對以上氣動優(yōu)化問題，本文結(jié)合SBO和POD-BPNN流場預(yù)測模型，提出一種基于PODBPNN的熱啟動策略并將其應(yīng)用于氣動代理優(yōu)化。以翼型氣動減阻優(yōu)化為例，優(yōu)化策略的總流程如圖1所示，具體分為以下5個步驟。

圖1 優(yōu)化流程Fig.1 Procedure of the optimization

（1）生成初始樣本庫

使用拉丁超立方抽樣（Latin Hypercube Sampling，簡稱LHS）在幾何設(shè)計(jì)空間中抽取r-1組幾何參數(shù)，并與基準(zhǔn)外形共同構(gòu)成初始外形庫，共計(jì)r個外形。然后生成基準(zhǔn)外形的氣動網(wǎng)格，并基于該網(wǎng)格利用RBF網(wǎng)格變形技術(shù)生成其他r-1個外形的氣動網(wǎng)格。RBF網(wǎng)格變形方法具有良好的魯棒性和較高的計(jì)算效率，可以解決復(fù)雜的網(wǎng)格變形問題，并在后續(xù)氣動優(yōu)化過程中也使用該方法獲得氣動網(wǎng)格。最后通過CFD計(jì)算模塊得到r個外形的流場數(shù)據(jù)和氣動力。

（2）建立初始?xì)鈩宇A(yù)測模型

使用Kriging模型直接建立從幾何參數(shù)到氣動力的氣動力預(yù)測模型Model 1，使用POD和BPNN建立從幾何參數(shù)到流場數(shù)據(jù)的流場預(yù)測模型Model 2，具體建模流程如圖2所示。

圖2 氣動預(yù)測模型Fig.2 Procedure of the aerodynamic prediction model

（3）流場預(yù)測及CFD熱啟動計(jì)算

對于優(yōu)化模塊得到的兩組新的設(shè)計(jì)變量，通過Model 2預(yù)測相應(yīng)的流場數(shù)據(jù)Y′和Y′，預(yù)測過程詳見1.5節(jié)；并將Y′和Y′作為CFD熱啟動計(jì)算的初場進(jìn)行熱啟動計(jì)算，計(jì)算得到真實(shí)氣動力p、p，真實(shí)目標(biāo)函數(shù)y、y和真實(shí)流場數(shù)據(jù)Y、Y。

（4）更新氣動預(yù)測模型

每完成一輪CFD熱啟動計(jì)算，判斷優(yōu)化是否終止，若判斷出優(yōu)化仍將繼續(xù)，則將新樣本添加到樣本庫。由于Kriging建模效率高，每完成一輪CFD熱啟動計(jì)算，就將(x，p，y)和(x，p，y)添加到Kriging模型并重新建立Model 1；而PODBPNN模型的建模效率較低，每完成k輪CFD熱啟動 計(jì) 算，把2k個 樣 本 流 場 數(shù) 據(jù)(x，Y)，j=1，2，…，2k，添加到POD-BPNN模型并重新訓(xùn)練得到新的Model 2。

（5）重復(fù)進(jìn)行步驟（2）～步驟（4），直到優(yōu)化終止。

1.2 Kriging模型

Kriging模型具有較好的非線性擬合能力，能夠提供插值點(diǎn)的預(yù)測函數(shù)值和統(tǒng)計(jì)學(xué)誤差估計(jì)信息，該誤差估計(jì)信息可以用于動態(tài)加點(diǎn)，因此，本文使用Kriging模型用于氣動代理優(yōu)化過程。

Kriging模型可以表示為

Kriging模型的具體推導(dǎo)詳見文獻(xiàn)［5］。

1.3 POD模型

具有高維、復(fù)雜特征的流場數(shù)據(jù)可以通過降階模型實(shí)現(xiàn)其在低維空間的映射和流場特征的提取。在氣動優(yōu)化問題中會生成一定規(guī)模的高維流場數(shù)據(jù)，可以通過建立POD降階模型，將高維流場數(shù)據(jù)近似表示為若干階基模態(tài)的線性疊加。

當(dāng)前s階基模態(tài)的能量占比達(dá)到一定值時，X()

能夠以低誤差率近似表示為

式中：U′為U的偽逆矩陣。

通過POD降階模型實(shí)現(xiàn)了高維數(shù)據(jù)X到低維基模態(tài)系數(shù)ξ的映射。除了上述一般的POD求解方法，還有Snapshot-POD、基于SVD分解的POD等，它們的區(qū)別只是求解方法不同，但都會獲得相同的POD基模態(tài)等信息。通常，當(dāng)m?r時，使用一般的POD可以加快速度節(jié)省計(jì)算內(nèi)存；當(dāng)m?r時，使用Snapshot-POD可以加快速度節(jié)省計(jì)算內(nèi)存。

1.4 BPNN模型

BPNN作為神經(jīng)網(wǎng)絡(luò)的一種，對較低維度的多輸入多輸出數(shù)據(jù)具有較強(qiáng)的非線性擬合能力。典型的單隱含層BPNN架構(gòu)如圖3所示，圖中輸入層、隱含層和輸出層的神經(jīng)元分別用x(i=1，2，…，n)、z(j=1，2，…，b)和y(k=1，2，…，s)表示，該模型實(shí)現(xiàn)了n維輸入到s維輸出的建模。通過調(diào)整隱含層的層數(shù)、各層神經(jīng)元個數(shù)、各層間的激活函數(shù)類型、學(xué)習(xí)率等超參數(shù)，可以使模型具有較好的擬合能力和泛化能力，實(shí)現(xiàn)對新樣本輸出值的穩(wěn)健預(yù)測。BPNN的具體建模過程詳見文獻(xiàn)［27］。

圖3 BPNN架構(gòu)Fig.3 Architecture of BPNN

1.5 POD-BPNN流場預(yù)測模型

結(jié)合POD模型和BPNN模型，可以構(gòu)建出POD-BPNN流場預(yù)測模型。

針對在設(shè)計(jì)空間中抽樣獲取的r組二維外形，其幾何設(shè)計(jì)參數(shù)為x(i=1，2，…，r)，x為n維列向量，通過CFD計(jì)算得到對應(yīng)的流場數(shù)據(jù)：壓強(qiáng)P、速度U、速度V、溫度T、密度ρ，這五個物理場數(shù)據(jù)均為m×r階矩陣，m為氣動網(wǎng)格的網(wǎng)格數(shù)量。以 壓 強(qiáng) 場P=[P，P，…，P]（P為m維 列 向量）為例，對其建立POD模型，得到壓強(qiáng)場均值μ（m維列向量）、基模態(tài)系數(shù)ξ（s維列向量）、映射矩陣U（s×m階矩陣），且有：

對幾何設(shè)計(jì)參數(shù)x和基模態(tài)系數(shù)ξ建立BPNN模型，模型的輸入維度為n，輸出維度為s。把訓(xùn)練好的BPNN模型所表征的輸入輸出映射關(guān)系記為ψ，即ξ=ψ(x)，則通過POD和BPNN可以建立幾何設(shè)計(jì)變量x和壓強(qiáng)場P的映射關(guān)系：

式（7）為POD-BPNN流場預(yù)測模型。

對于設(shè)計(jì)空間的任一設(shè)計(jì)變量x，可以通過POD-BPNN模型中的BPNN預(yù)測其基模態(tài)系數(shù)ξ：

然后通過POD模型的μ和U重構(gòu)出壓強(qiáng)場的預(yù)測值P：

同樣可以通過對物理場U、V、T、ρ建立PODBPNN模型，從而快速獲得這些物理場的預(yù)測值，具體預(yù)測流程如圖4所示，模型的輸入為幾何設(shè)計(jì)變量，由優(yōu)化模塊給出；模型的輸出為各個物理量的流場預(yù)測數(shù)據(jù)，并作為計(jì)算初場傳遞給CFD熱啟動計(jì)算模塊。

圖4 使用POD-BPNN進(jìn)行流場預(yù)測Fig.4 Prediction of the flowfield by POD-BPNN

1.6 效率分析

對于氣動代理優(yōu)化中是否應(yīng)用基于PODBPNN的熱啟動策略進(jìn)行計(jì)算時間分析對比。優(yōu)化的初始樣本個數(shù)為r，優(yōu)化過程共調(diào)用λk+b次優(yōu)化模塊，其中λ，k∈N，b=0，1，…，k-1。當(dāng)不使用基于POD-BPNN的熱啟動策略時，優(yōu)化總時間T為

當(dāng)使用基于POD-BPNN的熱啟動策略時，優(yōu)化總時間T為

可以得到如下關(guān)系：

使用基于POD-BPNN的熱啟動策略，氣動代理優(yōu)化的時間節(jié)省量Δ為

2 案例分析

2.1 優(yōu)化問題描述

本文以RAE2822翼型為優(yōu)化的基準(zhǔn)翼型，優(yōu)化 設(shè) 計(jì) 工 況：Ma=0.734，α=2.79°，Re=6.5×10，優(yōu)化問題的數(shù)學(xué)模型為

式中：C、C和C分別為當(dāng)前翼型的升力、阻力和俯仰力矩系數(shù)；C和C分別為基準(zhǔn)翼型的升力和俯仰力矩系數(shù)；A和A分別為當(dāng)前翼型和基準(zhǔn)翼型的面積，為無量綱化數(shù)值。

2.2 參數(shù)化及網(wǎng)格收斂性分析

本文采用CST（Class-function Shape-function Transformation）參數(shù)化方法，采用6階伯恩斯坦多項(xiàng)式，共14個設(shè)計(jì)變量?；鶞?zhǔn)翼型網(wǎng)格收斂性分析如表1所示，表中網(wǎng)格均采用結(jié)構(gòu)/非結(jié)構(gòu)混合網(wǎng)格形式，這四套網(wǎng)格的物面網(wǎng)格疏密程度不同，根據(jù)C和C的數(shù)值對比確定出氣動優(yōu)化所使用的氣動網(wǎng)格為Ⅲ，網(wǎng)格數(shù)量約為4.7萬，遠(yuǎn)場邊界為70c（c為翼型弦長，這里c取無量綱數(shù)1）；物面第一層網(wǎng)格高度為5.0×10c。采用網(wǎng)格Ⅲ得到的CFD仿真結(jié)果與實(shí)驗(yàn)值以及文獻(xiàn)［7］中仿真結(jié)果對比如表2所示，翼面壓力系數(shù)分布與實(shí)驗(yàn)值的對比如圖5所示，可以看出：本文的升力系數(shù)計(jì)算結(jié)果與實(shí)驗(yàn)值相差不到3%；阻力系數(shù)與實(shí)驗(yàn)值相差不到5%；力矩系數(shù)與實(shí)驗(yàn)值相差不到11%，與文獻(xiàn)［7］相差不到3%；壓力系數(shù)分布曲線與實(shí)驗(yàn)值基本吻合，說明本文的CFD計(jì)算具有一定準(zhǔn)確性。

表1 基準(zhǔn)翼型的氣動網(wǎng)格收斂性分析Table 1 Grid convergence analysis for the baseline

表2 基準(zhǔn)翼型的力系數(shù)計(jì)算結(jié)果對比Table 2 Validation of CFD simulation for the baseline

圖5 RAE2822翼型壓力系數(shù)分布與實(shí)驗(yàn)值對比Fig.5 Comparison of pressure coefficient distribution of RAE2822

基準(zhǔn)網(wǎng)格為網(wǎng)格Ⅲ。使用RBF網(wǎng)格變形技術(shù)獲得其他外形的氣動網(wǎng)格，網(wǎng)格變形前后的氣動網(wǎng)格對比如圖6所示，可以看出：變形后網(wǎng)格的網(wǎng)格質(zhì)量與基準(zhǔn)網(wǎng)格基本保持一致。

圖6 網(wǎng)格對比Fig.6 Comparison of grids

2.3 優(yōu)化參數(shù)設(shè)置

初始樣本通過LSH獲取，抽取個數(shù)為n( )n-1 2，n為設(shè)計(jì)變量個數(shù)，剔除計(jì)算不收斂樣本，得到有效樣本89個。優(yōu)化中參數(shù)k的取值為5，即每完成5輪優(yōu)化模塊的調(diào)用，進(jìn)行一次PODBPNN預(yù)測模型的更新，既可以保證流場的預(yù)測精度，又能夠避免由于重復(fù)建立POD-BPNN而降低模型效率。

2.4 優(yōu)化結(jié)果及分析

目標(biāo)函數(shù)隨著加點(diǎn)個數(shù)的收斂歷程如圖7所示。優(yōu)化前后的氣動性能對比如表3所示，可以看出：對比基準(zhǔn)翼型，優(yōu)化后阻力系數(shù)C降低了約0.006 3，優(yōu)化得到的翼型相比于基準(zhǔn)翼型阻力系數(shù)減小了35.8%，升力系數(shù)、面積和力矩系數(shù)均滿足約束條件。

圖7 目標(biāo)函數(shù)的收斂曲線Fig.7 Convergence curve of goal function

表3 優(yōu)化前后的氣動性能對比Table 3 Compare of aerodynamic performance of baseline and optimum shape

優(yōu)化前后的翼型幾何形狀和翼面壓力系數(shù)分布對比如圖8～圖9所示。

圖8 優(yōu)化前后翼型對比Fig.8 Compare of baseline and optimum shape

圖9 優(yōu)化前后翼面壓力系數(shù)分布對比Fig.9 Comparison of pressure coefficient distribution of baseline and optimum shape

優(yōu)化前后翼型的壓強(qiáng)場云圖對比如圖10所示，可以看出：優(yōu)化后上表面的強(qiáng)激波得到有效削弱，減阻優(yōu)化效果明顯。

圖10 優(yōu)化前后壓強(qiáng)云圖對比Fig.10 Comparison of pressure field of baseline and optimum shape

由于POD-BPNN僅在SBO中加入樣本預(yù)測流場及CFD熱啟動計(jì)算，不會對SBO的相關(guān)參數(shù)、優(yōu)化過程和優(yōu)化結(jié)果產(chǎn)生影響。因此僅對本文策略和SBO進(jìn)行效率對比分析，結(jié)果如表4所示，可以看出：本文策略相比于SBO，單次CFD計(jì)算時間降低約68%，整體效率提升約37%。

表4 計(jì)算時間分析Table 4 Time-consuming analysis

在SBO中使用POD-BPNN對某一新加樣本點(diǎn)的流場預(yù)測效果如圖11所示，可以看出：PODBPNN在上翼面距離前緣60%c位置附近的壓強(qiáng)預(yù)測效果較差，而在非激波位置處的等壓線擬合均較好。

圖11 POD-BPNN對壓強(qiáng)場的預(yù)測效果Fig.11 Forecast result of pressure field by POD-BPNN

給出最優(yōu)外形下兩個方法的收斂曲線如圖12所示，可以看出：CFD需要約900步完成C和C的計(jì)算收斂，而使用基于熱啟動策略的CFD計(jì)算僅需200步左右就可以完成C和C的計(jì)算收斂。

圖12 兩種方法的收斂曲線Fig.12 The convergence curves of the two methods

在整個代理優(yōu)化過程中，初始樣本和新加樣本點(diǎn)的CFD計(jì)算的收斂步數(shù)如圖13所示，其中使用CFD計(jì)算的初始樣本收斂步數(shù)，數(shù)值集中在800步左右；采用基于POD-BPNN熱啟動策略的新加樣本點(diǎn)的計(jì)算收斂步數(shù)，數(shù)值集中在200步左右，新加樣本點(diǎn)在使用基于POD-BPNN模型的CFD熱啟動計(jì)算后，其收斂步數(shù)整體上低于使用CFD計(jì)算的初始樣本收斂步數(shù)，CFD熱啟動計(jì)算失敗的情況主要集中在加點(diǎn)的前期。熱啟動計(jì)算失敗的樣本量為9個，占整個熱啟動樣本量的5.8%，當(dāng)在計(jì)算中出現(xiàn)該情況，只需切換為CFD計(jì)算即可。

圖13 氣動計(jì)算的收斂步數(shù)分析Fig.13 Analysis of convergence steps for aerodynamic calculation

POD-BPNN模型對每個新加樣本點(diǎn)壓強(qiáng)場的預(yù)測均方根誤差如圖14所示，可以看出：預(yù)測誤差呈現(xiàn)初期大、后期小的趨勢，這與圖13中CFD熱啟動失敗集中在加點(diǎn)前期的現(xiàn)象是一致的。

圖14 POD-BPNN模型對壓強(qiáng)場的預(yù)測誤差Fig.14 The RMSE of pressure flowfiled by POD-BPNN

3 結(jié)論

（1）針對翼型跨聲速減阻優(yōu)化算例，使用基于POD-BPNN模型的熱啟動策略進(jìn)行氣動代理優(yōu)化，獲得了減阻明顯的優(yōu)化翼型，在優(yōu)化中單次CFD計(jì)算時間降低了68%，整體優(yōu)化效率提升了37%。

（2）由于POD-BPNN模型的建立相比于Kriging需要耗費(fèi)一定時間，但是Kriging模型隨著優(yōu)化模塊的調(diào)用實(shí)時更新；而POD-BPNN模型在調(diào)用多次優(yōu)化模塊后再進(jìn)行更新，可以有效提高建模效率。

（3）POD-BPNN流場預(yù)測模型中使用了POD，而POD的建模效率很高，因此POD-BPNN的建模效率與網(wǎng)格量的相關(guān)性較小。在氣動網(wǎng)格數(shù)目較多的三維外形氣動優(yōu)化中，該策略的引入使得優(yōu)化效率的提升更加顯著。

（4）使用POD-BPNN模型的預(yù)測流場作為初場進(jìn)行CFD熱啟動計(jì)算可能會產(chǎn)生計(jì)算不收斂的情況導(dǎo)致計(jì)算失敗，該比例不超過6%。可以通過提高模型預(yù)測精度來有效降低熱啟動計(jì)算失敗的概率。