基于遺傳算法的潮流能水輪機(jī)葉片翼型優(yōu)化設(shè)計(jì)研究?

袁 鵬, 王旭超， 王樹杰， 譚俊哲， 司先才， 邊冰冰

(1.中國(guó)海洋大學(xué)工程學(xué)院，山東 青島 266100； 2. 山東省海洋工程重點(diǎn)實(shí)驗(yàn)室，山東 青島 266100)

葉片是潮流能水輪機(jī)直接承受水動(dòng)力并將其轉(zhuǎn)化為機(jī)械能的部件，其性能對(duì)潮流能水輪機(jī)的獲能效率有重要的影響。而葉片的水動(dòng)力學(xué)性能與葉片各剖面的翼型密切相關(guān)[1]。因此，開發(fā)適合潮流能水輪機(jī)性能要求的專用翼型對(duì)提高潮流能水輪機(jī)的獲能效率有重要意義。

目前，對(duì)翼型的研究多集中在飛行器和風(fēng)力機(jī)領(lǐng)域。飛行器上應(yīng)用的主要是美國(guó)國(guó)家航空航天局(NASA)開發(fā)的NACA系列翼型，如NACA44系列、NACA63-2系列、NACA64-4系列等。為了適應(yīng)風(fēng)力機(jī)的運(yùn)行特點(diǎn)，從1980年代開始，一些風(fēng)電技術(shù)發(fā)達(dá)的國(guó)家相繼開始了風(fēng)力機(jī)專用翼型的研究，并取得了顯著的成果?，F(xiàn)在應(yīng)用比較成熟的風(fēng)力機(jī)翼型主要包括美國(guó)可再生能源實(shí)驗(yàn)室開發(fā)的NREL-S系列翼型[2]、荷蘭代爾夫特工業(yè)大學(xué)(Delft university of technology)發(fā)展的DU系列翼型[3]、丹麥RIS國(guó)家重點(diǎn)實(shí)驗(yàn)室開發(fā)的RIS系列翼型[4]以及瑞典航空研究院(FFA)開發(fā)的FFA-W系列翼型[5]等。近幾年，國(guó)內(nèi)一些科研單位和科研院校也陸續(xù)開展了對(duì)風(fēng)力機(jī)專用翼型技術(shù)的研究，并且開發(fā)出了具有自主知識(shí)產(chǎn)權(quán)的風(fēng)力機(jī)專用翼型，主要包括中國(guó)科學(xué)院工程熱物理研究所開發(fā)的CAS-W系列翼型[6]、西北工業(yè)大學(xué)開發(fā)的NPU-WA系列翼型[7]、重慶大學(xué)設(shè)計(jì)的CQU-A系列翼型[8]等。

近年來，國(guó)內(nèi)外相關(guān)學(xué)者針對(duì)潮流能水輪機(jī)翼型開展了相關(guān)研究，Goundar等采用數(shù)值模擬和試驗(yàn)方法對(duì)HF-Sx、Fx63137、以及NACA63815翼型進(jìn)行了升阻特性對(duì)比[9]；英國(guó)南安普頓大學(xué)的Bahaj等對(duì)水平軸海流能發(fā)電機(jī)葉輪及葉輪翼型進(jìn)行了大量的試驗(yàn)及測(cè)試工作[10]；Grasso利用RFOIL軟件和序列二次規(guī)劃法相結(jié)合開展了翼型優(yōu)化設(shè)計(jì)[11]；Molland等利用XFOIL軟件對(duì)二維水翼的空泡問題進(jìn)行了研究[12]。任毅如、張?zhí)锾锏热嘶谶z傳算法建立了水輪機(jī)翼型優(yōu)化設(shè)計(jì)模型，并用XFOIL軟件對(duì)翼型的水動(dòng)力學(xué)性能進(jìn)行了評(píng)估[13]；朱國(guó)俊、馮建軍等人采用RBF神經(jīng)網(wǎng)絡(luò)與NSGA-Ⅱ相結(jié)合的現(xiàn)代優(yōu)化技術(shù)對(duì)潮流能水輪機(jī)的多工況優(yōu)化問題進(jìn)行了研究[14]；王儉超等人用CFD軟件Fluent對(duì)多種NACA系列翼型進(jìn)行了數(shù)值模擬研究，獲得了翼型的水動(dòng)力特性曲線[15]。這些研究極大地促進(jìn)了潮流能水輪機(jī)翼型技術(shù)的發(fā)展，本文在上述研究的基礎(chǔ)上，開展了潮流能水輪機(jī)翼型的優(yōu)化設(shè)計(jì)研究工作。

NACA4418翼型具有阻力系數(shù)小、最大升力系數(shù)高、扭距低等特點(diǎn)，在一些風(fēng)力機(jī)和潮流能水輪機(jī)中得到了應(yīng)用[16]，因此本文以NACA4418翼型為初始翼型，對(duì)初始翼型進(jìn)行參數(shù)化擬合，確定目標(biāo)函數(shù)、設(shè)計(jì)變量和約束條件，通過遺傳算法找到最優(yōu)解后，采用數(shù)值模擬的方法對(duì)初始翼型和優(yōu)化翼型進(jìn)行水動(dòng)力學(xué)分析，并討論、分析數(shù)值模擬結(jié)果。

1 潮流能水輪機(jī)翼型設(shè)計(jì)要求

目前，在潮流能水輪機(jī)葉片上主要應(yīng)用的是航空翼型和風(fēng)力機(jī)翼型，但是潮流能水輪機(jī)的工作介質(zhì)、工作條件與飛行器和風(fēng)力機(jī)有很大區(qū)別，本文將從潮流能水輪機(jī)翼型與航空翼型及潮流能水輪機(jī)翼型與風(fēng)力機(jī)翼型的區(qū)別方面來討論潮流能水輪機(jī)翼型的設(shè)計(jì)要求。

1.1 潮流能水輪機(jī)翼型與航空翼型的區(qū)別

潮流能水輪機(jī)翼型與航空翼型的主要區(qū)別在于以下幾個(gè)方面：

(1)飛行器主要是在高空條件下運(yùn)行，所以航空翼型是按表面光滑設(shè)計(jì)的；由于海水中浮游生物的附著會(huì)導(dǎo)致潮流能水輪機(jī)葉片表面變粗糙，所以在潮流能水輪機(jī)翼型設(shè)計(jì)過程中，必須考慮粗糙度敏感性對(duì)潮流能水輪機(jī)翼型水動(dòng)力學(xué)性能的影響。

(2)飛行器處于失速狀態(tài)時(shí)，會(huì)導(dǎo)致墜機(jī)的危險(xiǎn)，所以航空翼型在設(shè)計(jì)過程中應(yīng)盡量避免其處于失速狀態(tài)；而大攻角會(huì)引起流體與潮流能水輪機(jī)葉片翼型的分離，升力劇烈下降，導(dǎo)致水輪機(jī)的振動(dòng)[17]，因此潮流能水輪機(jī)翼型要具有良好的失速特性。

(3)飛行器一般在高空條件下運(yùn)行，高空空氣比較稀薄對(duì)葉片的沖擊較小，所以航空翼型的最大相對(duì)厚度一般不超過20%；而潮流能水輪機(jī)處于海水運(yùn)行條件下，海水對(duì)葉片的沖擊要比高空稀薄的空氣大的多，為了保證潮流能水輪機(jī)葉片的結(jié)構(gòu)和剛度要求，葉片根部翼型的厚度一般都比較大，有時(shí)最大相對(duì)厚度甚至達(dá)到40%以上[18]。

1.2 潮流能水輪機(jī)翼型與風(fēng)力機(jī)翼型的區(qū)別

潮流能水輪機(jī)翼型與風(fēng)力機(jī)翼型的主要區(qū)別在于以下幾個(gè)方面：

(1)風(fēng)力機(jī)葉片一般都比較細(xì)長(zhǎng)，葉片旋轉(zhuǎn)過程中會(huì)產(chǎn)生較大的扭轉(zhuǎn)力矩，所以力矩系數(shù)是風(fēng)力機(jī)翼型設(shè)計(jì)過程中必須需要考慮的一個(gè)重要參數(shù)；而潮流能水輪機(jī)葉片一般都比較剛硬，力矩系數(shù)并不是潮流能水輪機(jī)翼型設(shè)計(jì)過程中考慮的關(guān)鍵因素。

(2)由于風(fēng)速比較大，并且風(fēng)向和風(fēng)速變化都比較頻繁，在風(fēng)力機(jī)翼型設(shè)計(jì)過程中，一般選擇較高的設(shè)計(jì)升力系數(shù)來降低疲勞載荷和陣風(fēng)載荷；而潮流能水輪機(jī)所處的流場(chǎng)環(huán)境湍流度較低、流速和流向都具有較強(qiáng)的規(guī)律性，所以疲勞并不是潮流能水輪機(jī)翼型設(shè)計(jì)過程中考慮的主要因素。

(3)風(fēng)力機(jī)所處的運(yùn)行環(huán)境，一般風(fēng)速都比較大，并且變化頻繁，風(fēng)力機(jī)經(jīng)常在失速區(qū)運(yùn)轉(zhuǎn)，所以在風(fēng)力機(jī)翼型設(shè)計(jì)過程中，失速點(diǎn)的設(shè)計(jì)非常重要。而潮流能水輪機(jī)所處的運(yùn)行環(huán)境，雖然也會(huì)出現(xiàn)失速現(xiàn)象，但流速和流向具有較強(qiáng)的規(guī)律性，翼型失速點(diǎn)的設(shè)計(jì)并不是潮流能水輪機(jī)翼型設(shè)計(jì)過程中考慮的關(guān)鍵因素。

2 潮流能水輪機(jī)翼型設(shè)計(jì)方法

目前，常用的翼型優(yōu)化設(shè)計(jì)方法主要有正問題設(shè)計(jì)法和反問題設(shè)計(jì)法[19]。反問題設(shè)計(jì)法一般是以速度分布、壓力分布為目標(biāo)參數(shù)，通過迭代求解幾何和流動(dòng)控制方程來達(dá)到目標(biāo)的速度和壓力分布。這種方法難以處理多學(xué)科問題，并且翼型的控制點(diǎn)很難確定。正問題設(shè)計(jì)法以數(shù)值優(yōu)化方法為基礎(chǔ)，對(duì)初始翼型不斷進(jìn)行幾何修形，直至達(dá)到目標(biāo)升力或升阻比，這種設(shè)計(jì)方法可以改善多個(gè)設(shè)計(jì)參數(shù)并施加多種約束，隨著計(jì)算機(jī)技術(shù)的發(fā)展，正問題設(shè)計(jì)法得到了廣泛發(fā)展。

本文以在潮流能水輪機(jī)葉片上常用的NACA4418翼型為初始翼型，采用正問題設(shè)計(jì)法進(jìn)行翼型的優(yōu)化設(shè)計(jì)。用儒可夫斯基保角變換的方法對(duì)初始翼型進(jìn)行參數(shù)化擬合，確定目標(biāo)函數(shù)、設(shè)計(jì)變量和約束條件，通過尋優(yōu)算法尋找最優(yōu)解，找到最優(yōu)解后，對(duì)初始翼型和優(yōu)化翼型進(jìn)行水動(dòng)力學(xué)分析，并比較分析結(jié)果。具體的設(shè)計(jì)流程見圖1。

圖1 翼型設(shè)計(jì)流程圖 Fig.1 Flow chart of the hydrofoil design

3 潮流能水輪機(jī)翼型優(yōu)化設(shè)計(jì)模型的建立

3.1 翼型的參數(shù)化建模

翼型參數(shù)化方法的優(yōu)劣對(duì)翼型的設(shè)計(jì)效果有重要影響，它是決定翼型設(shè)計(jì)效率和翼型水動(dòng)力性能的關(guān)鍵因素。本文選用基于儒可夫斯基保角變換的方法對(duì)潮流能水輪機(jī)翼型進(jìn)行參數(shù)化表達(dá)。

利用儒可夫斯基保角變換式[20]：

ζ=f(x)=z+e2/z。

(1)

就可以將z平面上的一個(gè)圓變換為ζ平面上的一個(gè)翼型。其中e表示1/4翼型弦長(zhǎng)。

ζ平面上的翼型在笛卡爾坐標(biāo)系下就可以表示為：

(2)

其中:l表示翼型的矢徑長(zhǎng)度，將其表示為θ的函數(shù)：

l=e·exp(φ(θ))。

(3)

根據(jù)Taylor級(jí)數(shù)對(duì)的思想，將φ(θ)展開為：

φ(θ)=m1(1-cosθ)+n1sinθ+m2(1-cosθ)2+

n2sin2θ+L+mk(1-cosθ)k+nksinkθ+L,

k=1,2,3,L

(4)

從公式(4)可以看出當(dāng)θ=0時(shí)，φ(θ)=0，從而保證了翼型的尖尾緣特性。通過mk、nk的不同取值就可以得到不同的翼型。

3.2 目標(biāo)函數(shù)

潮流能水輪機(jī)葉片設(shè)計(jì)的主要目標(biāo)是使潮流能的功率系數(shù)CP達(dá)到最大，而CP通常是通過葉素理論來進(jìn)行計(jì)算。

選取距葉片旋轉(zhuǎn)中心距離為r處的葉素進(jìn)行研究，葉素上速度與受力分布見圖2。虛線為葉素所在的旋轉(zhuǎn)平面，V為來流速度，U為水流相對(duì)于葉素的合成來流速度，α為來流攻角，Ф是入流角，β為葉素扭角，c葉素弦長(zhǎng)。

圖2 葉素速度與受力圖Fig.2 Blade element velocity and load diagram

合成來流速度U為平行于旋轉(zhuǎn)平面的切向速度Vx與垂直于旋轉(zhuǎn)平面的法向速度Vy的合成速度，其中：

Vx=Ucosφ=(1+a′)Ωr。

(5)

式中：Ω為水輪轉(zhuǎn)動(dòng)角速度。

Vy=Usinφ=(1-a)V。

(6)

而水流流過水輪機(jī)產(chǎn)生的升力效應(yīng)會(huì)使水輪機(jī)旋轉(zhuǎn)做功，水流的軸向速度減小，并產(chǎn)生切向旋轉(zhuǎn)速度。定義a為軸向誘導(dǎo)因子，a′為切向誘導(dǎo)因子，分別表示水流軸向與切向速度的變化量，在穩(wěn)態(tài)運(yùn)行情況下，a、a′存在穩(wěn)定值[21]：

(7)

式中：λ為局部尖速比，為距旋轉(zhuǎn)中心r處的葉素旋轉(zhuǎn)速度與水流速度之比，即

(8)

根據(jù)庫(kù)塔儒可夫斯基升力方程可以推導(dǎo)出葉素微元的升阻力為：

(9)

(10)

為了使單個(gè)葉素的升、阻力的比值最大，即dL/dD最大，有公式(9)、(10)可知：

(11)

式中：CL為升力系數(shù)；CD為阻力系數(shù)。

由公式(11)可知要想使葉素的升、阻力的比值達(dá)到最大，就需要保證升阻比(CL/CD)最大，因此確定目標(biāo)函數(shù)為：

f(X)=max(CL/CD)。

(12)

3.3 設(shè)計(jì)變量

為了完整的表達(dá)翼型，同時(shí)減少計(jì)算量，本文選擇公式(4)中φ(θ)的前6項(xiàng)作為優(yōu)化設(shè)計(jì)變量，記作：

X=(x1,x2,x3,x4,x5,x6)。

(13)

3.4 約束條件

本文是針對(duì)NACA4418翼型進(jìn)行優(yōu)化，其最大相對(duì)厚度為18%，即施加約束：

t/c=18% 。

(14)

式中：t為翼型的最大厚度；c為翼型的弦長(zhǎng)。

NACA4418翼型的最大厚度所在的弦向位置在距前緣30%處，為了保證優(yōu)化翼型的幾何兼容性，施加約束：

Lmax=0.30。

(15)

考慮到翼型的強(qiáng)度要求，需要對(duì)它的截面積進(jìn)行約束：

(16)

式中：S為優(yōu)化翼型的截面面積；S0為初始翼型的截面面積。

減小翼型的前緣半徑可以使轉(zhuǎn)捩點(diǎn)位置前移[22]，從而保證光滑和粗糙條件下翼型的轉(zhuǎn)捩點(diǎn)位置基本相同，減小翼型的前緣粗糙度敏感性。因此，需要對(duì)翼型的前緣半徑進(jìn)行約束：

Le≤3% 。

(17)

潮流能水輪機(jī)作為一種大型旋轉(zhuǎn)類機(jī)械，需要對(duì)其進(jìn)行降噪處理，而翼型尾緣是主要的噪聲來源，根據(jù)翼型噪聲預(yù)測(cè)理論，其尾緣噪聲隨著尾緣厚度的增加而增大，因此，需要對(duì)翼型尾緣厚度進(jìn)行約束[23]：

yu,1-yl,1≤0.01。

(18)

式中：yu，1—x=1時(shí)翼型上表面的y坐標(biāo)值；yl，1—x=1時(shí)翼型下表面的y坐標(biāo)值。

增大翼型最大相對(duì)彎度可以提高翼型的升力系數(shù)和升阻比，因此本文對(duì)最大相對(duì)彎度進(jìn)行約束：

4%≤w/c≤6% 。

(19)

其中：w為翼型最大彎度；c為翼型弦長(zhǎng)。

翼型的最大相對(duì)彎度所處的弦向位置愈靠近翼型后緣，其升力系數(shù)愈高、阻力系數(shù)愈低、升阻比愈大，因此，需要對(duì)最大相對(duì)彎度所處的弦向位置進(jìn)行約束：

45%≤Lw≤50% 。

(20)

3.5 優(yōu)化設(shè)計(jì)程序

遺傳算法作為一種現(xiàn)代優(yōu)化方法，借鑒了達(dá)爾文的進(jìn)化論和孟德爾的遺傳學(xué)說，其本質(zhì)上是一種高效、全局搜索的方法[24]。本文采用遺傳算法進(jìn)行翼型的尋優(yōu)求解，以期獲得全局最優(yōu)解。該算法從一組初始種群開始求解，每一組解都采用二進(jìn)制形式進(jìn)行編碼，選擇一個(gè)適應(yīng)度函數(shù)，對(duì)每一個(gè)解的適應(yīng)度進(jìn)行評(píng)估，選出適應(yīng)度好的個(gè)體進(jìn)行選擇、交叉、變異，最后選出最優(yōu)個(gè)體。其具體的優(yōu)化求解過程見圖3。

圖3 翼型優(yōu)化設(shè)計(jì)流程圖Fig.3 Flow chart of the hydrofoil optimization and design

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

4.1 潮流能水輪機(jī)翼型優(yōu)化結(jié)果

本文以NACA4418翼型為初始翼型，通過在MATLAB中對(duì)遺傳算法進(jìn)行編程求解得到了一種新的翼型，圖4為初始翼型與優(yōu)化翼型的形狀對(duì)比。優(yōu)化后的翼型上翼面變厚，下翼面變薄翼型的彎度更大。

圖4 初始翼型與優(yōu)化翼型的形狀對(duì)比圖 Fig.4 Comparison diagram of the initial hydrofoil and optimized hydrofoil

4.2 潮流能水輪機(jī)翼型的水動(dòng)力學(xué)分析

優(yōu)化后的翼型需要進(jìn)一步通過水動(dòng)力學(xué)來確定其優(yōu)化效果，目前常用的翼型動(dòng)力學(xué)分析方法主要有兩種，一種是把粘性和無粘性相結(jié)合考慮的渦面元法，另外一種是CFD求解N-S方程的方法。

基于渦面元法和邊界層理論開發(fā)的XFOIL軟件，被廣泛應(yīng)用于翼型的設(shè)計(jì)和分析，并且用en法來預(yù)測(cè)轉(zhuǎn)捩點(diǎn)的位置，對(duì)薄翼型、小攻角的翼型水動(dòng)力學(xué)分析準(zhǔn)確性較高，但是在大攻角和翼型厚度比較厚的情況下準(zhǔn)確性降低。

CFD軟件Fluent被廣泛應(yīng)用于流體的動(dòng)力學(xué)分析，它不僅能準(zhǔn)確預(yù)測(cè)翼型的水動(dòng)力學(xué)性能，而且能得到一些相關(guān)的流場(chǎng)信息。

本文采用CFD軟件Fluent進(jìn)行翼型的水動(dòng)力學(xué)分析，其分析的關(guān)鍵是對(duì)翼型進(jìn)行網(wǎng)格劃分。本文用Gambit軟件進(jìn)行網(wǎng)格劃分，采用了C型結(jié)構(gòu)化網(wǎng)格，對(duì)翼型周圍進(jìn)行局部加密處理，翼型網(wǎng)格結(jié)構(gòu)見圖5。

在Fluent求解計(jì)算過程中假定來流速度1.5 m/s，通過改變x、y方向上的速度大小來計(jì)算-5°～25°攻角下的升、阻力系數(shù)。Standard k-ε湍流模型對(duì)二維翼型繞流數(shù)值計(jì)算具有較高的準(zhǔn)確性，尤其是對(duì)于大攻角分離流動(dòng)的計(jì)算精度高[25]，因此湍流模型選用Standard k-ε兩方程湍流模型。

圖5 NACA4418翼型整體網(wǎng)格圖和網(wǎng)格局部放大圖Fig.5 The whole grid diagram and the grid partial enlargement diagram of NACA4418 hydrofoil

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

通過CFD軟件Fluent的迭代求解計(jì)算，得到NACA4418翼型和優(yōu)化翼型在不同攻角下的升力系數(shù)曲線及升阻比曲線，見圖6、7。

圖6 NACA4418翼型和優(yōu)化翼型升力系數(shù)對(duì)比圖Fig.6 Lift coefficient comparison diagram of NACA4418 and the optimized hydrofoil

從圖6可以看出優(yōu)化后的翼型，在整個(gè)攻角范圍內(nèi)，都比初始翼型擁有更大升力系數(shù)。優(yōu)化后的翼型最大升力系數(shù)達(dá)到了1.017 8，相比于初始翼型的0.863 5提高了17.87%；并且，優(yōu)化翼型具有更緩和的失速特性。

從圖7可以看出優(yōu)化后的翼型相比于初始翼型在整個(gè)攻角范圍內(nèi)都擁有更高的升阻比。優(yōu)化后的翼型最大升阻比達(dá)到了15.351 3，相比于初始翼型的10.489 2提高了46.35%。

圖7 NACA4418翼型和優(yōu)化翼型升阻比對(duì)比圖Fig.7 Lift-to-drag ratio comparison diagram of NACA4418 and the optimized hydrofoil

由此可見，優(yōu)化后的翼型擁有更佳的水動(dòng)力學(xué)性能，有利于潮流能水輪機(jī)獲能效率的提高。

5 結(jié)論

通過分析潮流能水輪機(jī)翼型與航空翼型的區(qū)別以及潮流能水輪機(jī)翼型與風(fēng)力機(jī)翼型的區(qū)別，提出了潮流能水輪機(jī)翼型的設(shè)計(jì)要求，為潮流能水輪機(jī)翼型的設(shè)計(jì)提供了借鑒。

采用儒可夫斯基保角變換的方法對(duì)翼型進(jìn)行參數(shù)化表達(dá)，實(shí)現(xiàn)了翼型表達(dá)完備性和可控性的要求，并且保證了尖尾緣特性，為后續(xù)的翼型設(shè)計(jì)、研究提供了理論基礎(chǔ)。

采用遺傳算法進(jìn)行全局尋優(yōu)求解，通過設(shè)定多個(gè)約束條件來實(shí)現(xiàn)翼型的多種設(shè)計(jì)要求。求解出最優(yōu)解以后，通過CFD軟件Fluent進(jìn)行潮流能水輪機(jī)翼型繞流數(shù)值計(jì)算，驗(yàn)證了優(yōu)化翼型相比于初始翼型擁有更高的升阻比，證明了該翼型設(shè)計(jì)方法的可行性。