基于BP神經(jīng)網(wǎng)絡(luò)代理模型的翼型優(yōu)化及領(lǐng)域自適應(yīng)研究

陳晨銘， 郭雪巖， 常林森

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

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

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

領(lǐng)域自適應(yīng)可以將從一個(gè)樣本集學(xué)得的模型，遷移到另一個(gè)樣本集，從而完成泛化[9]。筆者針對青藏高原氣候特征的風(fēng)場(最高風(fēng)速為7 m/s左右)[10]，通過CST(Class function/Shape function Transformation)參數(shù)化方法實(shí)現(xiàn)翼型曲線與參數(shù)的轉(zhuǎn)換，使用深度前饋網(wǎng)絡(luò)構(gòu)建代理模型，結(jié)合遺傳算法對NACA64(3)-618風(fēng)力機(jī)翼型進(jìn)行優(yōu)化，并將代理模型用于NACA64(3)-218和NACA64(3)-418翼型，實(shí)現(xiàn)領(lǐng)域自適應(yīng)(或遷移學(xué)習(xí))。

1 機(jī)器學(xué)習(xí)方法

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

圖1 神經(jīng)網(wǎng)絡(luò)示意圖Fig.1 Schematic diagram of neural network

y=f(ωTx+θ)

(1)

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

型的參數(shù)。

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

(2)

(3)

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

(4)

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

2 翼型優(yōu)化設(shè)計(jì)

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

圖2 優(yōu)化設(shè)計(jì)流程圖Fig.2 Flow chart of optimal design

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

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

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

圖4 某一個(gè)樣本翼型的網(wǎng)格Fig.4 Mesh of a sample airfoil

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

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

表1 預(yù)測結(jié)果的最大相對誤差Tab.1 Maximum relative error of prediction results %

使用標(biāo)準(zhǔn)遺傳算法，優(yōu)化的目標(biāo)函數(shù)如下：

(5)

(6)

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

遺傳算法設(shè)置精度為10-10，種群最大迭代次數(shù)為50，共1 000個(gè)個(gè)體。式(5)中代理模型使用神經(jīng)網(wǎng)絡(luò)1和神經(jīng)網(wǎng)絡(luò)2分別進(jìn)行計(jì)算。式(6)中代理模型使用神經(jīng)網(wǎng)絡(luò)2進(jìn)行計(jì)算?？偣驳玫?個(gè)優(yōu)化翼型。

3 優(yōu)化翼型的氣動特性

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

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

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

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

圖8 優(yōu)化翼型3與擬合翼型壓力系數(shù)比較Fig.8 Comparison of pressure coefficient between optimized airfoil 3 and fitting airfoil

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

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

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

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

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

4 代理模型領(lǐng)域自適應(yīng)檢驗(yàn)

4.1 直接預(yù)測相似翼型

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

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

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

因?yàn)楣r相同、模擬方式相同，NACA64(3)-218、NACA64(3)-418和NACA64(3)-618翼型較為相似，所以阻力系數(shù)代理模型具有較好的精度。

4.2 修正系數(shù)調(diào)整

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

(7)

(8)

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

表3 領(lǐng)域自適應(yīng)方法的最大相對誤差Tab.3 Maximum relative error by domain adaption method %

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

4.3 結(jié)果分析

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

而翼型的升力系數(shù)是由上下表面的壓力差產(chǎn)生的，壓力差與翼型曲線關(guān)系較為復(fù)雜，故升力系數(shù)代理模型的精度較低。

5 結(jié) 論

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

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

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