基于BP神經網絡代理模型的翼型優化及領域自適應研究

陳晨銘， 郭雪巖， 常林森

(上海理工大學 能源與動力工程學院， 上海 200093)

風力機中的風輪葉片是獲得較大風力發電功率的關鍵，翼型設計是風力發電機組葉片研發的核心任務。航空翼型通常不能完全用作風力機翼型，一般要求風力機翼型在設計工況下具有較高的升力系數和升阻比，在其他工況下有良好的失速性能[1]。

翼型優化方法主要有2種：一種是正問題設計法，先設計翼型后再進行模擬，然后根據結果調整翼型；另一種是反問題設計法，根據需要的結果，推算出翼型[2]。韓忠華[3]總結了Kriging模型及代理優化算法的相關研究進展，如通過增加樣本點來提高代理模型精度的EI準則等。劉穎超等[4]通過Kriging代理模型優化汽輪機葉片頂端跨音速葉型。何磊等[5-6]通過深度學習的卷積神經網絡，以壓力圖像作為樣本輸入、翼型曲線作為樣本輸出，實現了反問題設計。Tao等[7]通過主成分分析-深度信念(PCA-DBN)網絡構造代理模型，采用拉丁超立方-適應度距離比-代理粒子群(LHS-FDR-SOPSO)算法對層流翼型和跨聲速機翼進行優化。李潤澤等[8]通過指定一個策略樹，使用強化學習方法，從策略樹中尋找最優策略，使得該策略下生成的翼型的氣動性能符合要求。

領域自適應可以將從一個樣本集學得的模型，遷移到另一個樣本集，從而完成泛化[9]。筆者針對青藏高原氣候特征的風場(最高風速為7 m/s左右)[10]，通過CST(Class function/Shape function Transformation)參數化方法實現翼型曲線與參數的轉換，使用深度前饋網絡構建代理模型，結合遺傳算法對NACA64(3)-618風力機翼型進行優化，并將代理模型用于NACA64(3)-218和NACA64(3)-418翼型，實現領域自適應(或遷移學習)。

1 機器學習方法

某個神經網絡如圖1所示，每一層都包含神經元，第一層為輸入層或可見層(v1～v3)，其中每一個神經元對應輸入的特征向量的一個值，通過連接輸送給下一層，中間層稱為隱藏層(h1和h2)，其神經元接受與其連接的上一層神經元傳遞的值，重新構成向量，通過M-P神經元模型進行計算，如式(1)所示，得到的值傳遞給下一層，并通過最后一層(即輸出層(o1))進行輸出。理論上，輸出可以是任意維度的向量，包含單個值。深度前饋網絡由含有多個隱藏層的神經網絡構成[9,11-12]。

圖1 神經網絡示意圖Fig.1 Schematic diagram of neural network

y=f(ωTx+θ)

(1)

式中：y為該神經元的輸出值；f(·)為某一函數，本文使用sigmoid函數；x為輸入的向量；ω和θ為模

型的參數。

選定f(·)后，ω和θ就是需要調整的參數，調整參數的過程又稱為訓練神經網絡。誤差逆傳播(error Back Propagation，BP)算法可以訓練一個神經網絡，但是由于標準BP算法只使用了梯度下降方法，而梯度下降方法可能收斂于局部最優而非接近全局最優，可能在梯度不明顯的地方反復振蕩，也可能在梯度突變的地方發散，所以BP算法也存在這一缺點[9,11-12]。LM(Levenberg-Marquardt)-BP算法是近似于牛頓法的正則化方法，牛頓法在一定條件下一直沿下降方向移動，有可能在其他情況下不收斂，而用了LM方法卻能達到收斂[13-14]。貝葉斯正則化(Bayesian regularization)-BP算法通過LM方法更新權值和偏值，最小化平方誤差和權重的線性組合，然后不斷修正該線性組合進行迭代，以獲得較好的神經網絡泛化性[15]。LM方法的平方誤差和權重的組合如式(2)～式(4)[16-17]所示。

(2)

(3)

M=αEW(w|A)+βED(D|w,A)

(4)

式中：D為含有N個樣本的樣本集，D={xm,tm}，其中xm和tm分別為第m個樣本的輸入特征向量和輸出向量；w為神經網絡參數；i表示維度；A為神經網絡結構；ED(D|w,A)為由以w和A構成的神經網絡計算得到的樣本集D預測值的0.5倍平方誤差和；y(xm,w,A)為以w和A構成的神經網絡計算得到的預測值；EW(w|A)為以A為結構的神經網絡中w的0.5倍L2范數平方；M為貝葉斯正則化函數；α和β為對應線性項的權重。

2 翼型優化設計

本文方法的優化設計流程圖如圖2所示，包括CST參數化方法、機器學習方法、代理模型和智能優化方法。首先通過CST參數化方法、B樣條曲線和拉丁超立方采樣生成1 000個翼型曲線樣本，然后使用計算流體力學(CFD)方法進行模擬計算，獲得升力系數Cl和阻力系數Cd，使用其中150個翼型樣本來訓練深度前饋網絡代理模型，通過代理模型和遺傳算法在給定區間內對升阻比Cl/Cd進行尋優，并通過CFD再次驗證獲得的最優解的精度。

圖2 優化設計流程圖Fig.2 Flow chart of optimal design

筆者選用NACA64(3)-618翼型，采用4階CST參數化方法進行參數擬合。在一定區間，根據拉丁超立方采樣，獲得1 000個樣本。之后通過CST參數化方法還原出301個坐標點數據，再將這些坐標點作為控制點，用B樣條曲線進行曲線的近似和光滑化，得到2 001個坐標點數據。樣本翼型區域與擬合翼型如圖3所示，其中X、Y分別為相對橫坐標和相對縱坐標。

圖3 樣本翼型區域與擬合翼型Fig.3 Sample airfoil region and fitting airfoil

本文使用的是C型網格，使用Pointwise軟件生成網格。由于參數化時生成的翼型是B樣條曲線，且有2 001個坐標點，不使用軟件自帶的基于Bézier原理生成的三次樣條插值方法[18]，其效果更好，通過設置參數，最后獲得翼型網格。某一個樣本翼型的網格如圖4所示。

圖4 某一個樣本翼型的網格Fig.4 Mesh of a sample airfoil

使用SSTk-ω湍流模型，以空氣作為流體，入口絕對速度為7.303 7 m/s，攻角為6°。使用SIMPLEC格式，動量方程使用三階MUSCL離散格式，湍動能和比耗散率均使用二階迎風離散格式，其余采用默認設置。迭代次數為2 000。

本文使用2種神經網絡，每種神經網絡再分為2個，分別輸出升力系數和阻力系數，總共生成2種、4個神經網絡代理模型。神經網絡1使用1層隱藏層，為10個節點，通過LM-BP算法進行訓練。神經網絡2使用2層隱藏層，各5個節點，通過貝葉斯正則化-BP算法進行訓練。150個樣本分為訓練集、驗證集和泛化集，每個集對應的樣本數分別為150個樣本數的70%、15%和15%。通過4階CST參數化方法將光滑化后的翼型曲線再次擬合成參數，作為神經網絡的輸入特征向量，采用之前CFD計算的升力系數和阻力系數分別作為輸出。最終得到神經網絡1升力系數代理模型、神經網絡1阻力系數代理模型以及對應的神經網絡2升力系數代理模型和神經網絡2阻力系數代理模型。4個神經網絡代理模型預測的最大相對誤差如表1所示。

表1 預測結果的最大相對誤差Tab.1 Maximum relative error of prediction results %

使用標準遺傳算法，優化的目標函數如下：

(5)

(6)

式中：b為CST參數；bmin和bmax為b的最小值和最大值；fl和fd分別為之前訓練的升力系數和阻力系數代理模型中的神經網絡函數，使用CST參數b作為輸入，升力系數Cl和阻力系數Cd分別作為輸出；f(b)為需要最優化的升阻比，符號為負；A1為翼型面積；A0為原翼型面積。

遺傳算法設置精度為10-10，種群最大迭代次數為50，共1 000個個體。式(5)中代理模型使用神經網絡1和神經網絡2分別進行計算。式(6)中代理模型使用神經網絡2進行計算。總共得到3個優化翼型。

3 優化翼型的氣動特性

通過神經網絡1和式(5)得出的優化翼型1與擬合翼型的對比曲線如圖5所示。通過神經網絡2和式(5)得出的優化翼型2與擬合翼型的對比曲線如圖6所示。通過神經網絡2和式(6)得出的優化翼型3與擬合翼型的對比曲線如圖7所示。采用CFD方法后，優化翼型3與擬合翼型壓力系數(Cp)對比如圖8所示，優化翼型3與擬合翼型上、下翼面摩阻系數(Cf)對比如圖9和圖10所示，優化翼型3的流線與壓力如圖11所示。上翼面平均壓力系數降低，下翼面平均壓力系數提高，導致升力系數提高。由于是低速翼型，僅在上翼面末端出現流動分離，上翼面末端和前駐點的摩阻系數有較大波動，阻力系數降低的原因較為復雜。

圖5 優化翼型1與擬合翼型比較Fig.5 Comparison between optimized airfoil 1 and fitting airfoil

圖6 優化翼型2與擬合翼型比較Fig.6 Comparison between optimized airfoil 2 and fitting airfoil

圖7 優化翼型3與擬合翼型比較Fig.7 Comparison between optimized airfoil 3 and fitting airfoil

圖8 優化翼型3與擬合翼型壓力系數比較Fig.8 Comparison of pressure coefficient between optimized airfoil 3 and fitting airfoil

圖9 優化翼型3與擬合翼型上翼面摩阻系數比較Fig.9 Comparison of skin friction coefficient on the upper surface between optimized airfoil 3 and fitting airfoil

圖10 優化翼型3與擬合翼型下翼面摩阻系數比較Fig.10 Comparison of skin friction coefficient on the lower surface between optimized airfoil 3 and fitting airfoil

圖11 優化翼型3的流線和壓力Fig.11 Streamline and pressure of optimized airfoil 3

由CFD模擬計算可得，擬合翼型的升力系數為1.022 01，阻力系數為0.021 37，升阻比為47.817 9。優化翼型1的升力系數為1.062 97，阻力系數為0.020 59，升阻比為51.638 0。與擬合翼型相比，優化翼型1的升力系數提高了4.01%，阻力系數降低了3.69%，升阻比提高了7.99%。優化翼型2的升力系數為1.065 04，阻力系數為0.020 60，升阻比為51.704 2。與擬合翼型相比，優化翼型2的升力系數提高了4.21%，阻力系數降低了3.62%，升阻比提高了8.13%。優化翼型3的升力系數為1.063 41，阻力系數為0.021 28，升阻比為49.980 4。與擬合翼型相比，優化翼型3的升力系數提高了4.05%，阻力系數降低了0.42%，升阻比提高了4.52%。

由代理模型計算得到：優化翼型1的升力系數為1.062 07，阻力系數為0.020 60，升阻比為51.561 4，升阻比相對誤差為0.148%；優化翼型2的升力系數為1.065 16，阻力系數為0.020 61，升阻比為51.685 1，升阻比相對誤差為0.037%；優化翼型3的升力系數為1.063 26，阻力系數為0.021 29，升阻比為49.951 9，升阻比相對誤差為0.057%。

4 代理模型領域自適應檢驗

4.1 直接預測相似翼型

直接用NACA64(3)-618翼型的神經網絡2的2個代理模型對NACA64(3)-218和NACA64(3)-418翼型的升力系數和阻力系數進行預測分析，在相同工況和CFD求解設置下對2種翼型進行模擬，并對結果進行分析。這種方法屬于領域自適應，通過一個樣本集生成的模型去預測另一個樣本集[9]。

通過LM-BP算法生成的代理模型預測得到NACA64(3)-218、NACA64(3)-418和NACA64(3)-618翼型在本文工況下的升力系數分別為0.901 9、0.921 4和1.016 3，相對誤差分別為18.45%、3.68%和5.852 3×10-5；阻力系數分別為0.018 2、0.019 5和0.021 3，相對誤差分別為0.75%、0.28%和0.02%。通過貝葉斯正則化-BP算法生成的代理模型預測得到的3種翼型的升力系數分別為0.898 0、0.915 6和1.016 4，相對誤差分別為17.94%、3.02%和1.710 0×10-5；阻力系數分別為0.018 2、0.019 5和0.021 3，相對誤差分別為0.40%、0.23%和0.03%。通過CFD模擬得到的3種翼型的升力系數分別為0.761 4、0.888 7和1.016 4，阻力系數分別為0.018 1、0.019 5和0.021 3。2種阻力系數代理模型的相對誤差如表2所示。

表2 阻力系數相對誤差Tab.2 Relative error of drag coefficient %

因為工況相同、模擬方式相同，NACA64(3)-218、NACA64(3)-418和NACA64(3)-618翼型較為相似，所以阻力系數代理模型具有較好的精度。

4.2 修正系數調整

根據前文優化流程，對NACA64(3)-218和NACA64(3)-418翼型各生成50個翼型的樣本集。對之前獲得的代理模型通過修正系數法進行調整，如式(7)和式(8)所示。

(7)

(8)

通過對NACA64(3)-218和NACA64(3)-418 2種翼型擬合后的翼型計算修正系數，分別對包含50個樣本的樣本集進行預測和CFD模擬檢驗，結果如表3所示。

表3 領域自適應方法的最大相對誤差Tab.3 Maximum relative error by domain adaption method %

通過修正系數法，阻力系數代理模型的精度更高，升力系數代理模型精度也有一定提升。但升阻比的精度不佳，主要原因是升力系數代理模型不精確，部分數據甚至不如直接用擬合后的翼型相關數據替換。

4.3 結果分析

由于翼型的阻力系數主要是由湍流邊界層的摩擦導致的。不可壓縮流動平板邊界層中表面的總摩阻系數與平板特征尺度有關，不可壓縮流動沿曲面邊界層速度分布與Pohlhausen型參數λ有關，其中λ為邊界層坐標點橫坐標x的函數，與物面的形狀密切關聯[19]，這些可以作為本阻力系數代理模型精度較高的一個參考依據。

而翼型的升力系數是由上下表面的壓力差產生的，壓力差與翼型曲線關系較為復雜，故升力系數代理模型的精度較低。

5 結 論

(1) NACA64(3)-618翼型的優化結果表明，優化后翼型的升力系數可以達到1.063 41，阻力系數為0.021 28，升阻比為49.980 4，升力系數提高了4.05%，阻力系數降低了0.42%，升阻比提高了4.52%。

(2) 本文優化方法將高維數據降低為低維數據，且保留了比較完整的翼型數據。

(3) 代理模型可以大幅提高預測的效率，并保持一定的精度。同時，代理模型也可以用于類似的NACA翼型的升力系數和阻力系數預測，其中阻力系數代理模型對其他翼型有較好的泛化性。代理模型能夠在較短時間內和較嚴苛條件下得到較好的翼型。