基于代理模型的風(fēng)力機(jī)翼型動(dòng)態(tài)失速優(yōu)化設(shè)計(jì)

收稿日期：2021-07-19

基金項(xiàng)目：國家自然科學(xué)基金（51976131；52006148）；上海市“科技創(chuàng)新行動(dòng)計(jì)劃”地方院校能力建設(shè)項(xiàng)目（19060502200）

通信作者：李 春（1963—），男，博士、教授、博士生導(dǎo)師，主要從事動(dòng)力機(jī)械及流體機(jī)械優(yōu)化設(shè)計(jì)、風(fēng)能利用等方面的研究。

lichunusst@163.com

DOI：10.19912/j.0254-0096.tynxb.2021-0847 文章編號：0254-0096（2023）06-0343-08

摘 要：為改善風(fēng)力機(jī)翼型動(dòng)態(tài)失速性能，利用代理模型方法替代計(jì)算流體力學(xué)（CFD）方法開展翼型動(dòng)態(tài)失速特性優(yōu)化設(shè)計(jì)。通過CST參數(shù)化方法構(gòu)建翼型幾何外形，采用優(yōu)化的拉丁超立方抽樣進(jìn)行試驗(yàn)設(shè)計(jì)，獲得樣本點(diǎn)處的氣動(dòng)力參數(shù)，建立高斯過程回歸模型，依據(jù)改善期望最大準(zhǔn)則增加樣本點(diǎn)，不斷提高模型精度。以降低風(fēng)力機(jī)翼型的平均力矩與阻力系數(shù)為優(yōu)化目標(biāo)，以平均升力系數(shù)不降為限制條件，采用受自然啟發(fā)的全局進(jìn)化類遺傳算法進(jìn)行尋優(yōu)。結(jié)果表明：與原始翼型相比，優(yōu)化翼型綜合氣動(dòng)性能更優(yōu)，尤其是平均阻力與平均力矩系數(shù)，分別減小9.57%與16.6%；此外，優(yōu)化翼型可抑制后緣渦向前緣發(fā)展，在一定程度上改善動(dòng)態(tài)失速。

關(guān)鍵詞：風(fēng)力機(jī)翼型；動(dòng)態(tài)失速；CST參數(shù)化；代理模型；翼型優(yōu)化

中圖分類號：TK83 """"""""""" """" 文獻(xiàn)標(biāo)志碼：A

0 引 言

葉片作為風(fēng)力機(jī)獲得旋轉(zhuǎn)力矩的關(guān)鍵部件直接決定風(fēng)力機(jī)的氣動(dòng)效率。對于葉片整體的氣動(dòng)性能，其剖面翼型輪廓起到?jīng)Q定性作用，尤其是葉片中部和根部的幾何外形設(shè)計(jì)［1-2］。翼型設(shè)計(jì)主要分為基于壓力分布的反設(shè)計(jì)方法和直接氣動(dòng)優(yōu)化設(shè)計(jì)方法［3］。壓力反設(shè)計(jì)方法的優(yōu)點(diǎn)是可依據(jù)目標(biāo)壓力分布，通過迭代求解，直接有效地獲得滿足氣動(dòng)目標(biāo)的幾何外形，缺點(diǎn)是獲得滿足目標(biāo)工況的壓力分布需設(shè)計(jì)人員豐富的經(jīng)驗(yàn)［4］。對于直接氣動(dòng)優(yōu)化設(shè)計(jì)方法，一般是將翼型優(yōu)化目標(biāo)轉(zhuǎn)化為數(shù)學(xué)語言，在滿足一定限制條件的前提下，結(jié)合最優(yōu)化方法設(shè)計(jì)出滿足目標(biāo)要求的幾何外形［5］。

近些年，得益于計(jì)算機(jī)技術(shù)的發(fā)展，一種基于機(jī)器學(xué)習(xí)方法的代理模型方法被應(yīng)用于優(yōu)化流程，其目的是近似流場計(jì)算，這極大提高了優(yōu)化效率。例如，張玄武等［6］使用級聯(lián)前向神經(jīng)網(wǎng)絡(luò)構(gòu)造代理模型替代流場計(jì)算并使用遺傳算法優(yōu)化，結(jié)果表明與僅使用計(jì)算流體力學(xué)（computational fluid dynamic，CFD）方法進(jìn)行流場計(jì)算相比，優(yōu)化時(shí)間減少90%。Eric等［7］對比基于不同神經(jīng)網(wǎng)絡(luò)的代理模型，發(fā)現(xiàn)相較于CFD計(jì)算，預(yù)測結(jié)果均方根誤差約為1%，表明代理模型方法有助于實(shí)現(xiàn)風(fēng)力機(jī)的一些快速優(yōu)化設(shè)計(jì)。韓中華等［8］提出一種高保真度多級分層Kriging代理模型并通過加點(diǎn)準(zhǔn)則提高模型精度，將其應(yīng)用于NACA0012與ONERA M6翼型氣動(dòng)優(yōu)化設(shè)計(jì)，發(fā)現(xiàn)其預(yù)測精度更高并顯著提高優(yōu)化效率。

為保證風(fēng)力機(jī)在額定發(fā)電功率附近持續(xù)運(yùn)行，且為應(yīng)對風(fēng)速與風(fēng)向的變化，葉片會(huì)進(jìn)行變槳運(yùn)動(dòng)以穩(wěn)定輸出功率與抑制較大的流動(dòng)分離，葉片攻角的變化會(huì)導(dǎo)致葉片中部及根部極易出現(xiàn)動(dòng)態(tài)失速現(xiàn)象［9］，此時(shí)為減小翼型升力驟降與力矩波動(dòng)，需設(shè)計(jì)一種能改善動(dòng)態(tài)失速的翼型。Vishal等［10］研究發(fā)現(xiàn)動(dòng)態(tài)失速現(xiàn)象會(huì)加劇翼型前緣渦的形成，產(chǎn)生對葉片不利的氣動(dòng)力與力矩，但其之后通過基于代理模型的氣動(dòng)外形優(yōu)化技術(shù)可極大地降低俯仰力矩波動(dòng)。動(dòng)態(tài)失速現(xiàn)象嚴(yán)重影響葉片風(fēng)能利用率與結(jié)構(gòu)強(qiáng)度，故針對風(fēng)力機(jī)翼型動(dòng)態(tài)失速的翼型優(yōu)化設(shè)計(jì)具有重要的理論意義與現(xiàn)實(shí)緊迫性。

現(xiàn)階段將代理模型方法應(yīng)用于靜態(tài)風(fēng)力機(jī)翼型優(yōu)化設(shè)計(jì)已取得較多成果［11-13］，目前利用代理模型技術(shù)針對風(fēng)力機(jī)翼型動(dòng)態(tài)失速優(yōu)化設(shè)計(jì)文獻(xiàn)較少，亟需進(jìn)一步展開研究。本文采用可進(jìn)行全局優(yōu)化的遺傳算法展開對風(fēng)力機(jī)翼型動(dòng)態(tài)失速特性的優(yōu)化設(shè)計(jì)，并使用高斯過程回歸模型構(gòu)建近似模型，之后利用改善期望（expected improvement， EI）加點(diǎn)準(zhǔn)則提高模型精度以加快優(yōu)化效率，以減小翼型阻力與力矩系數(shù)平均值為目標(biāo)函數(shù)，確保翼型在發(fā)生動(dòng)態(tài)失速時(shí)仍具有良好的氣動(dòng)與結(jié)構(gòu)穩(wěn)定性。

1 數(shù)值模擬方法

使用商業(yè)軟件ICEM生成計(jì)算網(wǎng)格，整體流域網(wǎng)格及局部細(xì)節(jié)如圖1所示。為便于流域網(wǎng)格劃分，整體流域采用正方形拓?fù)浣Y(jié)構(gòu)，計(jì)算域尺寸為[50c×75c]（[c]為翼型弦長），翼型旋轉(zhuǎn)域半徑為2c，通過滑移網(wǎng)格實(shí)現(xiàn)流域中心圓形旋轉(zhuǎn)域的運(yùn)動(dòng)，邊界條件設(shè)置為速度入口和壓力出口，翼型表面設(shè)置為無滑移壁面。網(wǎng)格單元總數(shù)為126480，壁面邊界層第一層網(wǎng)格高度為0.01 mm，確保y+值約等于1，壁面附近網(wǎng)格增長率設(shè)置為1.02。

采用商業(yè)軟件Fluent求解不可壓縮URANS方程，壓力與速度耦合求解使用SIMPLEC算法。SST [k-ω]湍流模型對自由剪切流和適度分離流均有較好的計(jì)算精度，文獻(xiàn)［14］對比了實(shí)驗(yàn)值與各種湍流模型，計(jì)算結(jié)果表明SST [k-ω]湍流模型可更好地契合中低雷諾數(shù)流動(dòng)。利用Fluent求解獲得100個(gè)初始樣本點(diǎn)氣動(dòng)力數(shù)據(jù)用于建立代理模型。

為驗(yàn)證數(shù)值計(jì)算方法的有效性，選取OSU風(fēng)洞實(shí)驗(yàn)結(jié)果，對具有實(shí)驗(yàn)數(shù)據(jù)S809翼型進(jìn)行動(dòng)態(tài)失速數(shù)值模擬［15］。對翼型輕度失速與重度失速下分別進(jìn)行驗(yàn)證，所選實(shí)驗(yàn)工況為：[Re=1×106]，振蕩頻率為0.6、1.8 Hz，振蕩頻率越高失速程度越嚴(yán)重，平均來流攻角[α0]為14°，攻角振幅[α1]為10°，來流攻角為[α=α0+α1sinωt]，其中[ω]為振蕩角速度。數(shù)值模擬與實(shí)驗(yàn)值在不同振蕩頻率下升力系數(shù)遲滯回環(huán)對比如圖2所示，遲滯回環(huán)為翼型在一振蕩周期內(nèi)不同攻角（angle of attack， AOA）的氣動(dòng)力數(shù)據(jù)。

由圖2a可知，翼型輕度失速下升力系數(shù)遲滯回環(huán)模擬值與實(shí)驗(yàn)值較為吻合。圖2b中翼型重度失速時(shí)上仰階段與實(shí)

驗(yàn)數(shù)據(jù)吻合較為良好，而在下俯階段有較大偏差，下附階段是動(dòng)態(tài)失速再附著的過程［16］，僅僅依靠二維非定常CFD計(jì)算難以滿足計(jì)算精度，這是由于二維CFD自身缺陷導(dǎo)致的。動(dòng)態(tài)失速再附著過程是三維強(qiáng)非線性過程，二維建模雖然不足以精確表達(dá)該過程所有特征，但也能捕獲大部分流動(dòng)特點(diǎn)［17］?？紤]優(yōu)化設(shè)計(jì)的時(shí)間成本，本文采用二維數(shù)值模擬計(jì)算對翼型進(jìn)行動(dòng)態(tài)失速優(yōu)化設(shè)計(jì)。

2 優(yōu)化設(shè)計(jì)方法

2.1 翼型參數(shù)化方法

CST（class function/shape function transformation）參數(shù)化方法由Kulfan［18］提出，因其嚴(yán)格的數(shù)學(xué)表達(dá)與明確的幾何意義，使其具有設(shè)計(jì)變量少、魯棒性好與精度高的特點(diǎn)，廣泛應(yīng)用在翼型參數(shù)化過程中［19］。CST參數(shù)化公式具體數(shù)學(xué)表達(dá)為：

[y=Cx?Sx+x?Δyte]""" （1）

式中：[x、][y]——翼型表面坐標(biāo)點(diǎn)；[Cx、][Sx]——類函數(shù)與形函數(shù)；[Δyte]——尾緣厚度。

上下翼面坐標(biāo)具體展開為：

[yu=CN1N2xSux+xyteuyl=CN1N2xSlx+xytel]" （2）

式中：[u、][l]——翼型上下表面。

類型函數(shù)定義為：

[CN1N2x=xN11-xN2]"" （3）

式中：[N1、][N2]——翼型幾何外形的類別。

對于本文前圓后尖類翼型，[N1]為0.5，[N2]為1，上下翼面形函數(shù)定義為：

[Sux=0nuAuiBinuxSlx=0nlAliBinlx]""" （4）

式中：[i]——多項(xiàng)式指數(shù)；[n u、][n l]——上下翼面多項(xiàng)式階數(shù)；[A u]、[Al]——上下翼面形函數(shù)的系數(shù)，其構(gòu)成幾何外形的參數(shù)向量，也是本文的設(shè)計(jì)變量；[Bin]——Bernstein多項(xiàng)式。

[Bin]數(shù)學(xué)表達(dá)為：

[Bin=Kinxi1-xn-i]"" （5）

式中：[Kin]——組合系數(shù)，[Kin=n！i！n-i！]。

利用Bernstein多項(xiàng)式逼近形函數(shù)，通過最小二乘法使得擬合誤差平方和達(dá)到最小。根據(jù)參考文獻(xiàn)［20］的對比結(jié)果，采用5階Bernstein多項(xiàng)式對幾何外形的擬合精度已完全滿足氣動(dòng)力計(jì)算，故對翼型上下表面均采用5階Bernstein多項(xiàng)式進(jìn)行擬合，共計(jì)12個(gè)設(shè)計(jì)變量。擬合曲線如圖3所示，上下表面擬合誤差如圖4所示。由圖4可見，經(jīng)CST參數(shù)化方法擬合的翼型最大誤差在翼型中部，僅為[1.947×10-3]，說明選擇5階CST參數(shù)化方法對翼型進(jìn)行擬合已滿足本文優(yōu)化設(shè)計(jì)的精度。

2.2 高斯過程回歸模型

2.2.1 試驗(yàn)設(shè)計(jì)方法

高斯過程回歸模型的建立，首先需一定數(shù)量的初始樣本集用于模型的訓(xùn)練，初始樣本集的選取又稱為試驗(yàn)設(shè)計(jì)。常見試驗(yàn)設(shè)計(jì)包括正交數(shù)組、均勻設(shè)計(jì)與拉丁超立方抽樣等。傳統(tǒng)的拉丁超立方抽樣會(huì)出現(xiàn)采樣點(diǎn)相對集中的缺陷，會(huì)對代理模型的構(gòu)建產(chǎn)生不利影響［21］。本文選取改進(jìn)的拉丁超立方采樣方法：在拉丁超立方采樣的基礎(chǔ)上，利用基于最大最小距離準(zhǔn)則連續(xù)枚舉方法生成樣本［22］，以二維空間為例，改進(jìn)后的拉丁超立方和拉丁超立方采樣方法對比如圖5所示。由圖5可知，改進(jìn)后樣本點(diǎn)的空間分布更為均勻，有利提高代理模型精度。

Latin square samples

2.2.2 高斯過程回歸模型

高斯過程回歸（Gaussian process regression，GPR）是一種以貝葉斯定理為框架的機(jī)器學(xué)習(xí)方法，對于處理小樣本與復(fù)雜隨機(jī)過程有較好適應(yīng)性［23］。GPR模型是通過統(tǒng)計(jì)學(xué)原理構(gòu)建的輸入與輸出樣本點(diǎn)之間的一種映射關(guān)系。設(shè)有一組樣本集[D=xi，yii=1，2，…，n]，其中[xi]為樣本點(diǎn)輸入向量，[yi]是樣本點(diǎn)響應(yīng)值，其為標(biāo)量［24］。模型方程可表示為：

[yi=fxi+εi]""" （6）

式中：[εi～N0，σ2n]，是均值為0，方差為[σ2n]，服從獨(dú)立同分布的高斯噪聲。

[fxi]定義了高斯過程，即：

[fxi～GPμxi，covxi，x′i] （7）

高斯過程由均值[μxi]和協(xié)方差函數(shù)[covxi，xi′]確定，具體定義為：

[μxi=Efxi] （8）

[covxi，xi′=Efxi-μxifxi′-μxi′]" （9）

在無噪聲項(xiàng)時(shí)，輸出值[fxi]的先驗(yàn)分布是[fxi～GP0，cov]。對于測試集[xi*]，根據(jù)[yi]的先驗(yàn)分布建立訓(xùn)練集輸出與測試集輸出[yi*]之間的聯(lián)合高斯先驗(yàn)分布：

[yiyi*～GP0，covxi，xicovxi，xi*covxi，xi*covxi*，xi*]""" （10）

式中：[covxi，xi]——[n×n]階對稱正定的協(xié)方差矩陣；[covxi，xi*]——測試與訓(xùn)練輸入樣本間的協(xié)方差矩陣；[covxi*，xi*]——測試集樣本間的協(xié)方差矩陣。

在[xi]與訓(xùn)練樣本集[D]已知時(shí)，由后驗(yàn)概率公式計(jì)算得到[xi*]對應(yīng)的[yi*]，具體為：

[yi*xi*，"" D～GPμyi*，Σyi*] （11）

式中：

[μyi*=covxi*， xicovxi*， xi*+σ2nI-1yi] （12）

[Σyi*=covxi*，xi-covxi*，xicovxi*，xi*+σ2nI-1covxi，xi*]""""""""" （13）

若要計(jì)算后驗(yàn)分布的均值[μyi*]和方差[Σyi*]，需求解4個(gè)協(xié)方差。協(xié)方差矩陣又稱為核函數(shù)，有多種不同形式的核函數(shù)可應(yīng)用于回歸模型，本文使用matern5/2核函數(shù)。

邊緣概率[pyi*xi*]是GPR計(jì)算中的重要參數(shù)，其定義為：

[pyi*xi*=N0，covxi，xi+σ2nI] （14）

用極大似然法求解超參數(shù)：

[Lθ=-lnpyixi"""""""" =12yiTC-1yi+12lnC+n2ln2π]""" （15）

式中：[C=covxi，xi+σ2nI;][C]——[C]的行列式。

求解[Lθ]通過調(diào)用軟件中的一些模塊實(shí)現(xiàn)，如Matlab中fminunc與Python中scipy.optimize等。

2.2.3 加點(diǎn)準(zhǔn)則

經(jīng)試驗(yàn)設(shè)計(jì)的有限樣本數(shù)據(jù)構(gòu)建的高斯回歸模型精度較低，在某些點(diǎn)的預(yù)測結(jié)果與真實(shí)值相比誤差較大，針對該問題，發(fā)展出多種加點(diǎn)準(zhǔn)則［25］。本文介紹并使用Jones［26］發(fā)展的EI加點(diǎn)優(yōu)化準(zhǔn)則。

設(shè)優(yōu)化目標(biāo)函數(shù)最小值為[ymin]，GPR預(yù)測結(jié)果[yx]服從均值為[yx]，方差為[σ2]的高斯分布，期望改善量[Ix]表示為：

[Ix=ymin-yx]""" （16）

目標(biāo)函數(shù)期望改善準(zhǔn)則為：

[EIx=ymin-yxΦymin-yxs+s?ymin-yxs，sgt;00，s=0]"""""""" （17）

式中：[Φ]——標(biāo)準(zhǔn)正太分布累積函數(shù)；[?]——標(biāo)準(zhǔn)正太分布概率密度函數(shù)。

本文加點(diǎn)優(yōu)化策略為：在超參數(shù)優(yōu)化過程中，在訓(xùn)練好的近似模型上尋找EI值最大的樣本點(diǎn)加入訓(xùn)練集，重新構(gòu)建代理模型直至生成一個(gè)精度更高的近似模型用于翼型優(yōu)化流程。

2.3 優(yōu)化流程

優(yōu)化流程涉及到大量重復(fù)性工作，如網(wǎng)格劃分、氣動(dòng)力數(shù)據(jù)處理與加點(diǎn)尋優(yōu)，使用軟件自帶腳本語言進(jìn)行批處理操作，在Matlab軟件中實(shí)現(xiàn)自動(dòng)優(yōu)化，極大加快優(yōu)化流程。圖6為動(dòng)態(tài)失速翼型優(yōu)化流程，其具體步驟為：

1）使用CST參數(shù)化方法建模獲取表示翼型特征的輸入向量[xi]，由形函數(shù)系數(shù)[Au]與[Al]組成行向量，優(yōu)化空間為[Au]與[Al]變化20%，圖7為翼型優(yōu)化空間。

2）試驗(yàn)設(shè)計(jì)的目的是使用少量樣本刻畫設(shè)計(jì)空間的主要特征。采用基于最大距離準(zhǔn)則優(yōu)化的LHS方法在空間內(nèi)選取100個(gè)樣本，通過CFD方法求解的氣動(dòng)力即為訓(xùn)練樣本的輸出值[yi]。

3）利用初始樣本集[xi，yi]訓(xùn)練GPR模型，將優(yōu)化目標(biāo)轉(zhuǎn)化為數(shù)學(xué)語言，并給定設(shè)計(jì)變量的約束條件與優(yōu)化范圍，為避免陷入局部最優(yōu)，選擇具有全局優(yōu)化能力的遺傳算法尋優(yōu)。

4）通過EI加點(diǎn)優(yōu)化準(zhǔn)則在設(shè)計(jì)空間內(nèi)優(yōu)化潛力最大處增加樣本點(diǎn)，隨著樣本點(diǎn)的增加，最優(yōu)解范圍內(nèi)的模型精度逐漸提高，與真實(shí)最小值誤差越來越小，在小于給定誤差或滿足最大增加樣本數(shù)量時(shí)跳出循環(huán)。

3 優(yōu)化結(jié)果與分析

3.1 算例設(shè)置

以風(fēng)力機(jī)翼型S809為研究對象，雷諾數(shù)為1×106，使翼型繞氣動(dòng)中心（[x/c=0.25]）周期性運(yùn)動(dòng)。振蕩規(guī)律為[αt=10+14sin3.6πt]。翼型動(dòng)態(tài)失速情況下，阻力與力矩往往伴隨著劇烈變化，對葉片的穩(wěn)定產(chǎn)生嚴(yán)重影響，本文的優(yōu)化目標(biāo)為降低阻力系數(shù)與力矩系數(shù)平均值，將多目標(biāo)優(yōu)化分配不同權(quán)重轉(zhuǎn)化為單目標(biāo)優(yōu)化，具體目標(biāo)函數(shù)與約束條件為：

[min0.6i=1NCmi+0.4i=1NCdi]" （18）

[s.t.0.8T0max≤Tmax≤1.2T0maxi=1NCdi≤i=1NC0dii=1NCmi≤i=1NC0mii=1NCli≥i=1NC0li]"" （19）

式中：N——1個(gè)振蕩周期內(nèi)不同攻角氣動(dòng)力參數(shù)的數(shù)量，本文選取4°～24°攻角1112個(gè)數(shù)據(jù)計(jì)算平均升力、阻力及力矩系數(shù)；[Tmax]——翼型最大厚度；上標(biāo)0為原始翼型的氣動(dòng)力參數(shù)。

采用遺傳算法進(jìn)行尋優(yōu)，設(shè)置初始種群數(shù)為1000，進(jìn)化代數(shù)為300，目標(biāo)函數(shù)一般在120～180代內(nèi)收斂，利用高斯過程回歸模型建立代理模型，之后通過EI加點(diǎn)準(zhǔn)則，加入60個(gè)樣本。平均升力、阻力及力矩系數(shù)誤差約為3.21%、3.45%及2.52%。

3.2 結(jié)果分析

動(dòng)態(tài)優(yōu)化翼型與原始翼型對比如圖8所示。由圖8可見，優(yōu)化翼型有更大的前緣半徑，最大厚度位置后移，翼型下表面厚度增加，對減小翼型力矩十分有利。

圖9分別為升阻力系數(shù)以及力矩系數(shù)遲滯回線對比，其中力矩系數(shù)正負(fù)表示翼型扭轉(zhuǎn)方向的不同，本文為繞氣動(dòng)中心逆時(shí)針為正。由圖9a可見，在上仰運(yùn)動(dòng)時(shí)，優(yōu)化翼型升力有所提高，尤其在攻角較大時(shí)有顯著提升；在下附階段，升力的振蕩幅度有所減弱。由圖9b可知，在上仰階段，阻力系數(shù)在攻角較小時(shí)并未變化，在攻角大于18°后，阻力系數(shù)明顯降低；在下附階段，阻力系數(shù)不僅減小且振蕩幅度也減弱。由圖9c可知，優(yōu)化翼型的力矩系數(shù)在上仰與下附階段均有明顯減小。對比升力、阻力與力矩系數(shù)，優(yōu)化翼型的俯仰力矩減小最為明顯，升力與阻力系數(shù)在大攻角時(shí)改善較為明顯。為對比優(yōu)化翼型與基準(zhǔn)翼型具體氣動(dòng)力變化幅度，表1給出一振蕩周期內(nèi)平均氣動(dòng)力與厚度變化率。由表1可知，作為優(yōu)化目標(biāo)的阻力與力矩系數(shù)變化最為顯著，尤其是力矩系數(shù)，減小16.6%，阻力與扭轉(zhuǎn)力矩的減小可保證翼型的穩(wěn)定與結(jié)構(gòu)強(qiáng)度。

圖10分別為當(dāng)上仰來流攻角為20.65°時(shí)，原始翼型與優(yōu)化翼型的速度流線對比圖。由圖10流線可知，優(yōu)化翼型局部輪廓的改變雖未能改變后緣流動(dòng)分離現(xiàn)象，但其抑制了后緣渦向翼型前緣發(fā)展。前緣渦的減小從側(cè)面解釋了優(yōu)化翼型的升力、阻力及力矩系數(shù)遲滯回環(huán)比較平穩(wěn)的變化狀態(tài)。

b.優(yōu)化翼型

airfoil and optimized airfoil

本算例的優(yōu)化結(jié)果顯示，采用本文方法設(shè)計(jì)得到的新翼型氣動(dòng)性能有明顯改善，證明了該方法的可行性。

4 結(jié) 論

本文提出一種基于代理模型和啟發(fā)式全局優(yōu)化算法相結(jié)合的翼型動(dòng)態(tài)失速優(yōu)化設(shè)計(jì)方法，得到綜合性能更優(yōu)的新翼型，并結(jié)合流場進(jìn)行對比分析，得到主要結(jié)論如下：

1）優(yōu)化翼型的平均升力系數(shù)提升幅度較小，僅為3.51%，但其在大攻角時(shí)提升效果顯著；平均阻力與力矩系數(shù)改善較為明顯，分別減小9.57%與16.6%，尤其是力矩系數(shù)，在整個(gè)振蕩周期內(nèi)皆有較為明顯的減弱，可大幅提高翼型穩(wěn)定性與機(jī)構(gòu)強(qiáng)度。

2）對翼型局部輪廓的改變，可獲得在特定工況下氣動(dòng)性能更優(yōu)的氣動(dòng)外形，相較于原始翼型，優(yōu)化翼型在大攻角下對于后緣渦向前緣的發(fā)展有一定的抑制作用，在一定程度上減弱了翼型動(dòng)態(tài)失速現(xiàn)象。

［參考文獻(xiàn)］

［1］" 李春， 葉舟， 高偉， 等. 現(xiàn)代陸海風(fēng)力機(jī)計(jì)算與仿真［M］. 上海： 上?？茖W(xué)技術(shù)出版社， 2012.

LI C， YE Z， GAO W， et al. Computation and simulation of modern land-sea wind turbine［M］. Shanghai： Shanghai Scientific amp; technical publishers， 2012.

［2］" 陳進(jìn)， 陳剛， 謝翌. 考慮彎扭變形的風(fēng)力機(jī)葉片結(jié)構(gòu)優(yōu)化［J］. 太陽能學(xué)報(bào)， 2018， 39（4）： 1119-1126.

CHEN J， CHEN G， XIE Y. Structure optimization of wind turbine blade considering bend and twist deformation［J］. Acta energiae solaris sinica， 2018， 39（4）： 1119-1126.

［3］" 吳琪. 基于粘性伴隨方法的旋翼先進(jìn)氣動(dòng)外形優(yōu)化設(shè)計(jì)分析［D］. 南京： 南京航空航天大學(xué)， 2014.

WU Q. Optimal design and analysis on advanced aerodynamic shape of rotor based on a viscous adjoint method［D］. Nanjing： Nanjing University of Aerodynamics and Astronautics， 2014.

［4］" 劉俊， 宋文萍， 韓忠華， 等. Kriging模型在翼型反設(shè)計(jì)中的應(yīng)用研究［J］. 空氣動(dòng)力學(xué)學(xué)報(bào)， 2014， 32（4）： 518-526.

LIU J， SONG W P， HAN Z H， et al. Kriging-based airfoil inverse design［J］. Acta aerodynamica sinica， 2014， 32（4）： 518-526.

［5］" 張鑫帥， 劉俊， 羅世彬. 基于改進(jìn)多目標(biāo)布谷鳥搜索算法的翼型多目標(biāo)氣動(dòng)優(yōu)化設(shè)計(jì)［J］. 航空學(xué)報(bào)， 2019， 40（5）： 122550.

ZHANG X S， LIU J， LUO S B. An improved multi-objective cuckoo search algorithm for airfoil aerodynamic shape""" optimization""" design［J］."" Acta"" aeronautica""" et astronautica sinica， 2019， 40（5）： 122550.

［6］" 張玄武， 鄭耀， 楊波威， 等. 基于級聯(lián)前向網(wǎng)絡(luò)的翼型優(yōu)化設(shè)計(jì)［J］. 浙江大學(xué)學(xué)報(bào)（工學(xué)版）， 2017， 51（7）： 1405-1411.

ZHANG X W， ZHENG Y， YANG B W， et al. Aerodynamic optimization design of airfoil configurations based on cascade feedforward neural network［J］. Journal of Zhejiang University（engineering science edition）， 2017， 51（7）： 1405-1411.

［7］" ERIC R L， BENJAMIN V， GIRMA B， et al. Comparison of neural network types and architectures for generating a surrogate"" aerodynamic"" wind"" turbine" blade"" model［J］. Journal of wind engineering amp; industrial aerodynamics， 2021， 216： 104696.

［8］" HAN Z H， XU C Z， ZHANG L， et al. Efficient aerodynamic shape optimization using variable-fidelity surrogate models and multilevel computational grids［J］. China journal of aeronautics， 2020， 33（1）： 31-47.

［9］" DE T， FERRIRA C， VIRE A， et al. Controlling dynamic stall using vortex generators on a wind turbine airfoil［J］. Renewable energy， 2021， 172： 1194-1211.

［10］" VISUAL" R，" LEIFUR" L." Surrogate-based" aerodynamic shape optimization for delaying airfoil dynamic stall using kriging regression and infill criteria［J］. Aerospace science and technology， 2021， 111： 106555.

［11］" CHEN J， WANG Q， ZHANG S Q， et al. A new direct design method of wind turbine airfoils and wind tunnel experiment［J］. Applied mathematical modelling， 2016， 40（3）： 2002-2014.

［12］" PHOLDEE N， BUREERAT S， NUANTONG W. Kriging surrogate-based genetic algorithm optimization for blade design of a horizontal axis wind amp; sciences， 2021， 126（1）： 261-273.

［13］" OH S. Comparison of a response surface method and artificial neural network in predicting the aerodynamic performance of a wind turbine airfoil and its optimization［J］. Applied sciences-basel， 2020， 10（18）： 6277.

［14］" BELAMADI R， DJEMILI A， ILINCA A， et al. Aerodynamic performance analysis of slotted airfoils for application to wind turbine blades［J］. Journal of wind engineering amp; industrial aerodynamics， 2016， 151： 79-99.

［15］" RAMSAY R， HOFFMAN M J， GREGOREK G M. Effects of grit roughness and pitch oscillations on the S809 airfoil［R］. Golden， Colorado： National Renewable Energy Laboratory， 1999： 1-165.

［16］" CORKE T C， THMOAS F O. Dynamic stall in pitching airfoils： aerodynamic damping and compressibility effects［J］. Annual review of fluid mechanics， 2015， 47（1）： 479-505.

［17］" 喻伯平， 李高華， 謝亮， 等. 基于代理模型的旋翼翼型動(dòng)態(tài)失速優(yōu)化設(shè)計(jì)［J］. 浙江大學(xué)學(xué)報(bào)（工學(xué)版）， 2020， 54（4）： 833-842.

YU B P， LI G H， XIE L， et al. Dynamic stall optimization design of rotor aifoil based on surrogate model［J］. Journal of Zhejiang University（engineering science）， 2020， 54（4）： 833-842.

［18］" KULFAN B M. Universal parametric geometry representation method［J］. Journal of aircraft， 2008， 45（1）： 142-158.

［19］" 尹國慶， 王軍， 王威， 等. 基于CST參數(shù)化方法的軸流風(fēng)機(jī)多目標(biāo)優(yōu)化設(shè)計(jì)［J］. 風(fēng)機(jī)技術(shù)， 2020， 62（6）： 45-51.

YIN G Q， WANG J， WANG W， et al. Multi-objective optimization of axial flow fan based on CST parameterization""""" method［J］."""" Chinese"""" journal"""" of turbomachinery， 2020， 62（6）： 45-51.

［20］" 廖炎平， 劉莉， 龍騰. 幾種翼型參數(shù)化方法研究［J］. 彈箭與制導(dǎo)學(xué)報(bào)， 2011， 31（3）： 160-164.

LIAO Y P， LIU L， LONG T. The research on some parameterized methods for airfoil［J］. Journal of projectiles rockets missiles and guidance， 2011， 31（3）： 160-164.

［21］" 韓忠華. Kriging模型及代理優(yōu)化算法研究進(jìn)展［J］. 航空學(xué)報(bào)， 2016， 37（11）： 3197-3225.

HAN Z H. Kriging surrogate model and its application to design optimization： a review of recent progress［J］. Acta aeronautica et astronautica sinica， 2016， 37（11）： 3197-3225.

［22］" 葉鵬程， 潘光， 高山. 一種快速優(yōu)化拉丁超立方試驗(yàn)設(shè)計(jì)方法［J］. 西北工業(yè)大學(xué)學(xué)報(bào)， 2019， 37（4）： 714-723.

YE P C， PAN G， GAO S. Sampling design method of fast optimal"" Latin"" hypercube［J］." Journal"" of""" Northwestern Polytechnical University， 2019， 37（4）： 714-723.

［23］" 陳文英. 概率論與數(shù)理統(tǒng)計(jì)［M］. 北京： 科學(xué)出版社， 2012.

CHEN W Y. Probability and mathematical statistics［M］. Beijing： Science Press， 2012.

［24］" 吳寬展. 基于多輸出高斯過程回歸的超臨界翼型優(yōu)化［D］. 南京： 南京航空航天大學(xué)， 2015.

WU K Z. A supercritical airfoil design based on multi-output surrogate model［D］. Nanjing： Nanjing University of Aeronautics and Astronautics， 2015.

［25］" LIU J， HAN Z H， SONG W P. Comparison of infill sampling criteria in kriging-based aerodynamic optimization［C］//28th Congress of the International Council of the Aeronautical Sciences， Brisbane， Australia， 2012： 23-28.

［26］" JONES D R. Efficient global optimization of expensive black-box" functions［J］." Journal" of" global" optimization， 1998， 13（4）： 455-492.

OPTIMAL DESIGN OF DYNAMIC STALL OF WIND TURBINE

AIRFOIL BASED ON SURROGATE MODEL

Zhang Qiang1，Miao Weipao1，Chang Linsen1，Liu Qingsong1，Li Chun1，2，Zhang Wanfu1，2

（1. School of Energy and Power Engineering， University of Shanghai for Science and Technology， Shanghai 200093， China；

2. Shanghai Key Laboratory of Multiphase Flow and Heat Transfer in Power Engineering， Shanghai 200093， China）

Abstract：In order to improve the dynamic stall performance of wind turbine airfoil， the surrogate model method was used to replace CFD calculation to optimize the dynamic stall characteristics of wind turbine airfoil. The airfoil geometry profile was constructed by CST parameterization method， and the aerodynamic parameters at the sample points were obtained by using the optimized Latin hypercube sampling for experimental design. The Gaussian process regression model was established， and the sample points were added according to the maximum improvement expectation criterion to continuously improve the model accuracy. With the reduction of the average torque and drag coefficient of the wind turbine airfoil as the optimization objective and the non-reduction of the average lift coefficient as the limiting condition， Global evolutionary genetic algorithm inspired by nature is used to search for optimization. The results show that compared with the original airfoil， optimization of airfoil aerodynamic performance is better， especially the drag and torque coefficient are reduced 9.57% and 16.6%， respectively； In addition， the development of trailing edge vortex is inhibited by the optimized airfoil， and the dynamic stall is improved to some extent.

Keywords：wind turbine airfoil； dynamic stall； CST parameterization； surrogate model； airfoil optimization