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

(1.中國海洋大學(xué)工程學(xué)院，山東 青島 266100；2.青島市海洋可再生能源重點實驗室，山東 青島 266100)

潮流能是海洋可再生能源的一個重要門類，開發(fā)潮流能對于緩解中國的能源短缺和環(huán)境污染問題具有重要意義[1]。水平軸水輪機是目前應(yīng)用最為廣泛的一種潮流能捕獲裝置，其葉片直接承受水動力并將其轉(zhuǎn)化為轉(zhuǎn)子的機械能，葉片的性能與各剖面的翼型密切相關(guān)[2]。

國內(nèi)外研究機構(gòu)對風(fēng)力機葉片翼型進行了相關(guān)研究。荷蘭Delft大學(xué)設(shè)計出DU風(fēng)力機翼型系列[3]；丹麥RIS?國家實驗室研發(fā)出適用于風(fēng)力機的專用翼型，命名為RIS?翼型系列[4]；美國可再生能源實驗室(NREL)相繼開發(fā)出9個不同的NREL-S翼型系列[5]；瑞典航空研究院對翼型的研究也進行多年，研究開發(fā)出FFA-W1-XXX、FFA-W2-XXX和FFA-W3-XXX三個翼型族系列[6]。國內(nèi)也相繼開始對翼型進行了研究，西北工業(yè)大學(xué)開發(fā)出NPU-WA風(fēng)力機翼型族系列[7]；中國科學(xué)院工程熱物理研究開發(fā)出CAS-W1-XXX風(fēng)力機翼型族系列[8-9]。另外，重慶大學(xué)、北京航空航天大學(xué)等高校也在風(fēng)力機翼型研究方面取得一定的成果[10]。

在水平軸潮流能葉片設(shè)計方面，目前普遍采用的多是航空翼型或風(fēng)力機專用翼型[2]。為提高潮流能水輪機工作性能，開發(fā)適合潮流能水輪機工作環(huán)境的專用翼型是很有必要的。為提高潮流能水輪機的設(shè)計效率，中國海洋大學(xué)的王樹杰等[11]利用Visual C++開發(fā)了水平軸潮流能水輪機設(shè)計軟件(HAT design)。本文采用多目標優(yōu)化遺傳算法對潮流能葉片進行優(yōu)化設(shè)計，將優(yōu)化后的專用翼型數(shù)據(jù)應(yīng)用到水平軸潮流能水輪機設(shè)計軟件(HAT design)中，實現(xiàn)對專用翼型數(shù)據(jù)的調(diào)用，提高潮流能水輪機專用葉片的設(shè)計效率。

1 多目標遺傳算法

遺傳算法是模擬生物在自然環(huán)境中的遺傳和進化過程而形成的一種自適應(yīng)全局優(yōu)化概率搜索算法[12]。多目標問題的基本思想就是求解在滿足所有約束條件和各個目標函數(shù)的條件下的一組最優(yōu)解集[13]。其通用的數(shù)學(xué)表達式為：

minF(x)=[f1(x),f2(x),L,fm(x)]

(1)

式中：m,n-目標函數(shù)，設(shè)計變量個數(shù)；k,p-不等式約束，等式約束個數(shù)；ai、bi第i個設(shè)計變量的上限值和下限值[14]。遺傳算法步驟如下：

(1)初始種群選取。本文初始翼型的翼型坐標為初始化種群。設(shè)種群的大小為N，每一個個體按r+1遞增。選用三階B樣條曲線，因此控制點個數(shù)至少為4個，本文選取控制點個數(shù)為13個。

u=(1:n-k)=sort(G(1:n-k))。
n=length(G)-length(Q)。

其中：Q-初始數(shù)據(jù)序列；k-B樣條次數(shù)；length(G)-計算序列G中元素個數(shù)；sort(G)-對序列G排序；X(n:m)-序列X中第n個到第m個組成的序列[15]。

(2)適應(yīng)度函數(shù)。因翼型曲線上的控制點是隨著不同的參數(shù)選擇而變化的。擬合的翼型曲線不一定為最優(yōu)結(jié)果，所以需要另加約束，采用近弧長約束[15]，B(t)的弧長為：

(2)

化簡計算得：

(3)

該約束使所得到的擬合曲線為最優(yōu)擬合曲線。

(4)剪切和拼接算子。種群中的染色體可以隨著剪切拼接操作而變化，并產(chǎn)生新染色體，染色體也是可以變化的，這樣在擬合過程中產(chǎn)生更多的控制點，可以使樣條曲線擬合曲線匹配到更多的合適點。過程為：首先要選擇一個剪切點r1=Rand(1,L1);r2=Rand(1,L2)，其中Rand(n,m)表示隨機產(chǎn)生的[n,m]中的隨機的一個數(shù)，并且此數(shù)為整數(shù)[15]。在拼接和剪切過程中，容易產(chǎn)生非法染色體，為了避免此問題，需要對剪切點位置進行限制，因此需編寫程序?qū)崿F(xiàn)限制，具體程序如下：

L1=Length(G1);L2=Length(G2) ;

r1=Rand(1,L1);r2=Rand(1,L2);

while(

(r1+L2-r2)

r1=Rand(1,L1);

r2=Rand(1,L2);

)end

從而實現(xiàn)對翼型曲線弧長進行約束，使翼型曲線更加光順，避免擬合曲線擬合過程中出現(xiàn)劇烈震蕩。

(5)變異。通過改變每個個體的一些基因來恢復(fù)種群群體的多樣性，使搜索達到整個空間，這就是變異算子的操作[15]。

算法流程：

步驟1 輸入B樣條曲線次數(shù)K、平面有序坐標點{Pj(xj,yj),j=1,2,…,r}，種群大小為N，選擇概率p0，最大迭代次數(shù)tmax，初始種群迭代次數(shù)記為t=1[12]。取控制點個數(shù)為13，種群大小N=40，tmax=300。

步驟3 按照適應(yīng)度函數(shù)大則被選擇機率大的原則，對當(dāng)代種群進行隨機抽樣，得到N個新選擇的個體，計算每個個體的適應(yīng)度值，種群的平均適應(yīng)度值和最大、最小適應(yīng)度值[12]。

步驟4 將種群中個體隨機配對，計算交叉概率，運用線性交叉，得到新個體[12]，既為得到新的翼型坐標點，交叉概率取0.8。

步驟5 對交叉后的種群，先計算種群中每個個體的適應(yīng)度，種群的平均和最大適應(yīng)度值，從而計算出每個個體的變異概率pm，再對每個個體的每個基因執(zhí)行單點變異操作，變異概率取0.05。

步驟6 對變異后的種群，采用最差去除策略，將當(dāng)代種群中適應(yīng)度較小的個體去除，并把先前的最優(yōu)個體集添加到當(dāng)代種群中，這樣可保證算法最終收斂于全局最優(yōu)[12]。

步驟7 令t=t+1，重復(fù)步驟，直到滿足條件為止?？捎胻大于最大進化代數(shù)tmax或者誤差平方和SSE小于給定誤差作為迭代的終止條件，數(shù)據(jù)點誤差平方和為0.002。

2 翼型參數(shù)化建模

2.1 B樣條曲線擬合翼型曲線

翼型的曲線形狀為前圓后尖、表面光滑的流線圖形，要對翼型的曲線進行優(yōu)化，就要先了解分析翼型的幾何特性。采用B樣條曲線對翼型曲線進行擬合，首先要確定翼型曲線上的控制點，則此處設(shè)翼型上的坐標控制點為P={pj(xj,yj),j=1,K,r}，其次選擇參數(shù)向量，此處的參數(shù)向量為T={tj,j=1,K,r}，在非遞減節(jié)點向量上，定義一條k次B樣條曲線可表示為[16]:

U={u1=u2=…≤uk+1≤……≤un+k+1}。

(4)

節(jié)點向量定義為：

(5)

其中：D={di,i=1,K,n}表示控制點，Ni,k(t)為基函數(shù)，定義為：

(6)

選用三次B樣條曲線對翼型進行參數(shù)化表達，通過選擇翼型上的控制點就可以控制翼型曲線，用B樣條曲線擬合翼型曲線得到點能使得最小二乘擬合誤差平方和達到最小。則滿足條件的B樣條曲線寫成矩陣形式為[16]：

Q=ND。

(7)

其中：N為B樣條基函數(shù)陣，D為有序數(shù)據(jù)點陣。若事先確定參數(shù)向量T和節(jié)點向量U，可求得控制頂點陣D的近似解[12]：

D=(NTN)-1NTP。

(8)

則得到翼型擬合曲線：

PC=N(NTN)-1NTP。

(9)

(10)

如果求出控制頂點集D，使得：

p1=B(t1)=D1;
pr=B(tr)=Dn;
Pc=minSSE。

(11)

則稱該翼型擬合曲線Pc為插值于端點的最小二乘擬合曲線，此次選取13個控制點來進行翼型曲線的擬合[17]。

2.2 翼型曲線優(yōu)化模型

翼型設(shè)計與多個學(xué)科領(lǐng)域相關(guān)，比如可靠性設(shè)計、結(jié)構(gòu)設(shè)計、流體力學(xué)等。在設(shè)計過程中也會遇到翼型的結(jié)構(gòu)需求與水動力性能相矛盾的情況。例如：追求高的升力系數(shù)會得到更加敏感的前緣，翼型的升阻比會受到翼型相對厚度增加的限制[8]。多目標優(yōu)化問題的計算過程中，各子目標之間是相互競爭的，對翼型的水動力性能進行優(yōu)化設(shè)計必須考慮到將會出現(xiàn)的矛盾沖突，協(xié)調(diào)好將會出現(xiàn)的沖突和矛盾，從而得到最優(yōu)解[14]。

(1)設(shè)計變量：用B樣條曲線進行參數(shù)化建模，以翼型擬合曲線上的控制點為設(shè)計變量，設(shè)計變量為13個控制點，將設(shè)計變量表達為：

P=(P1，P2，P3，P4，P5，P6，P7，
P8，P9，P10，P11，P12，P13)。

(2)目標函數(shù)：設(shè)計目標主要是初始翼型取得最大升阻比和最大升力系數(shù)的攻角處，取得最優(yōu)解。翼型在設(shè)計攻角處的升阻比最大作為目標函數(shù)之一：

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

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

翼型在設(shè)計攻角處的升力系數(shù)最大也作為目標函數(shù)：

f1(X)=max(CL)。

(3)約束條件：多目標優(yōu)化時，目標之間有時會相互沖突，因此，在滿足優(yōu)化目標的同時還需約束其他因素。在設(shè)計攻角處，升阻比尋得最優(yōu)解時，升力系數(shù)也必須加以限制，不得小于初始翼型的升力系數(shù)，即CL≥X0；同理，當(dāng)設(shè)計攻角處升力系數(shù)最大時，升阻比不得小于初始翼型的升阻比值，即CL/CD≥Y0。

最大相對厚度為x，為了保證優(yōu)化后的翼型的相對厚度仍是x，施加約束：

t/c=x。

(12)

式中：t—翼型的最大厚度；c—翼型弦長。

初始翼型的最大厚度所在位置是距前緣35%～50%弦長位置處，因此施加約束[10]：

35%≤xt/c≤50%。

(13)

式中，xt—最大厚度所在位置距前緣的距離。

翼型的彎度與翼型的升阻比與強度有關(guān)，彎度增大可增大翼型升阻比，但又要對彎度進行限制[10]，因為彎度過大又影響翼型的強度，綜合以上因素施加約束[9]：

3%≤f/c≤5%。

(14)

式中，f為翼型彎度。

彎度位置適當(dāng)向后緣移動，有利于提高潮流能水輪機翼型的最大升阻比，施加約束[10]：

50%≤xj/c≤55%。

(15)

式中，xf為彎度位置距翼型前緣的距離。

考慮到翼型的強度要求，需要對它的截面積進行約束：

(16)

式中：S-優(yōu)化后翼型的截面積；S0-優(yōu)化前翼型的截面積。

翼型的尾緣是噪聲的主要來源，其尾緣噪聲隨著尾緣厚度的增加而增大，因此，需要對尾緣進行約束[18]：

yu,0.99-y1,0.99≤0.01。

(17)

式中yu,0.99、y1,0.99為默認弦長為1時，弦長位置處上翼面及下翼面的y坐標值[18]。

3 翼型設(shè)計結(jié)果

分別以FX60-126、NACA65019、NACA63-215、EPPLER 862 STRUT和EPPLER 864 STRUT為初始翼型，用B樣條曲線參數(shù)化建模，通過多目標遺傳算法求的五種新翼型，最大相對厚度分別為12.6%、19%、21.5%、32.4%和38.8%，分別命名為OUC-TT-1260、OUC-TT-1900、OUC-TT-2150、OUC-TT-3240和OUC-TT-3880。五種新翼型的幾何形狀如圖1所示。

圖1 OUC-TT-XXX0翼型族幾何形狀圖Fig.1 Geometrical shape of OUC-TT-XXX0 hydrofoil family

(1) OUC-TT-1260翼型優(yōu)化結(jié)果

在來流速度V=1.5 m/s工作環(huán)境下，因FX60-126翼型在α=15°時的升力系數(shù)最大，在α=2°時升阻比最大，因此選擇α=15°時的升力系數(shù)和α=2°時的最大升阻比作為優(yōu)化目標。以FX60-126翼型為初始翼型，通過用遺傳算法進行編程求解得到新的潮流能水輪機專用翼型OUC-TT-1260曲線坐標。通過對FX60-126翼型和翼型OUC-TT-1260進行網(wǎng)格劃分和CFD分析得到如圖2、3所示的仿真分析結(jié)果，水動力特性參數(shù)對比如圖4所示。

圖2 專用翼型多目標優(yōu)化前后翼型附近壓力云圖(V=1.5 m/s，α=2°)Fig.2 Pressure clouds around the original hydrofoil and optimized special hydrofoil(V=1.5 m/s，α=2°)

如圖4所示為OUC-TT-1260翼型與FX60-126翼型的水動力參數(shù)對比，如圖2、3所示為CFD仿真分析結(jié)果，由CFD仿真分析結(jié)果和水動力參數(shù)結(jié)果對比可知，OUC-TT-1260翼型的升阻比和最大升力系數(shù)比FX60-126翼型整體較高，攻角范圍在10°～20°時，OUC-TT-1260翼型的升力系數(shù)和升阻比明顯高于FX60-126翼型，則優(yōu)化結(jié)果符合要求。

圖3 專用翼型多目標優(yōu)化前后翼型附近壓力云圖(V=1.5 m/s，α=15°)Fig.3 Pressure clouds around the original hydrofoil and optimized special hydrofoil(V=1.5 m/s，α=15°)

圖4 OUC-TT-1260翼型與FX60-126翼型水動力參數(shù)對比Fig.4 Comparison of hydrodynamic parameters between OUC-TT-1260 and FX60-126 hydrofoil

(2) OUC-TT-1900翼型優(yōu)化結(jié)果

以NACA65019翼型為初始翼型，因?qū)ACA65019翼型進行分析得到的水動力性能參數(shù)，升力系數(shù)在α=20° 時取得最大值，在α=10° 時，升阻比取得最大值，因此選擇這兩種攻角環(huán)境下翼型的水動力學(xué)性能為優(yōu)化目標。

通過用遺傳算法進行編程求解得到新的潮流能水輪機專用翼型OUC-TT-1900曲線坐標，得到潮流能水輪機專用翼型OUC-TT-1900后進行網(wǎng)格劃分和CFD分析得到如圖5、圖6所示的仿真分析結(jié)果，水動力特性參數(shù)對比如圖7所示。

圖5 專用翼型多目標優(yōu)化前后翼型附近壓力云圖(V=1.5 m/s，α=10°)Fig.5 Pressure clouds around the original hydrofoil and optimized special hydrofoil (V=1.5 m/s，α=10°)

圖6 專用翼型多目標優(yōu)化前后翼型附近壓力云圖(V=1.5 m/s，α=20°)Fig.6 Pressure clouds around the original hydrofoil and optimized special hydrofoil (V=1.5 m/s，α=20°)

通過對如圖5、圖6所示的CFD仿真分析結(jié)果進行分析，對如圖7所示為OUC-TT-1900翼型與NACA65019翼型的水動力參數(shù)進行對比可知，OUC-TT-1900翼型的升阻比和最大升力系數(shù)比NACA 65019翼型整體較高，攻角范圍在15°～25°時，OUC-TT-1900翼型升力系數(shù)明顯高于NACA 65019翼型，優(yōu)化結(jié)果符合要求。

(3)OUC-TT-XXX0翼型族主要水動力參數(shù)結(jié)果

五種翼型的升力系數(shù)和升阻比都有相應(yīng)的提高，其中OUC-TT-1260翼型的升力系數(shù)提高最多，OUC-TT-1260翼型的升阻比提高最多。

由表1可知，優(yōu)化后的OUC-TT-XXX0翼型族具有較好的水動力特性。

圖7 OUC-TT-1900翼型與NACA65019翼型水動力參數(shù)對比Fig.7 Comparison of hydrodynamic parameters between OUC-TT-1900 and NACA65019 hydrofoil

表1 OUC-TT-XXX0翼型族主要水動力特性參數(shù)Table 1 Major parameters of hydrodynamic performance of OUC-TT-XXX0 hydrofoil family

4 專用翼型族數(shù)據(jù)庫設(shè)計及應(yīng)用

專用翼型族是針對潮流能水輪機工作環(huán)境而獲得的翼型系列。由以上研究可知，專用翼型族具有較好的水動力性能。水輪機的整體性能取決于翼型的水動力性能參數(shù)，翼型的水動力性能參數(shù)復(fù)雜繁多，搭建數(shù)據(jù)庫進行翼型數(shù)據(jù)管理，可以更加高效地獲取翼型的性能參數(shù)，從而快速地設(shè)計水輪機葉片[11]。

4.1 數(shù)據(jù)庫設(shè)計方案

為了實現(xiàn)在水平軸潮流能水輪機設(shè)計軟件中調(diào)用專用翼型族數(shù)據(jù)，將優(yōu)化后的潮流能水輪機專用翼型族數(shù)據(jù)進行整理，導(dǎo)入到SQL Server中，搭建專用翼型族數(shù)據(jù)庫。在水平軸潮流能水輪機設(shè)計軟件中對數(shù)據(jù)庫中的專用翼型族數(shù)據(jù)進行調(diào)用，完成后續(xù)的專用潮流能水輪機葉片的設(shè)計。水平軸潮流能水輪機設(shè)計軟件及其專用葉片翼型數(shù)據(jù)庫的整體設(shè)計方案如圖8所示。

圖8 專用翼型族數(shù)據(jù)庫設(shè)計方案Fig.8 Design scheme of special hydrofoil Database

4.2 運行結(jié)果

通過編程完成了水平軸潮流能水輪機設(shè)計軟件與專用翼型族數(shù)據(jù)庫之間的連接和數(shù)據(jù)調(diào)用。在水平軸潮流能水輪機設(shè)計軟件的翼型數(shù)據(jù)界面中，通過下拉框選擇專用翼型族數(shù)據(jù)庫中已有的翼型，從數(shù)據(jù)庫中調(diào)用專用翼型的幾何數(shù)據(jù)和水動力學(xué)參數(shù)等信息，通過二維圖形顯示翼型的外形，并顯示翼型的水動力特性曲線。潮流能水輪機設(shè)計軟件翼型數(shù)據(jù)界面如圖9所示。水平軸潮流能水輪機設(shè)計軟件也實現(xiàn)了葉片設(shè)計、葉片載荷分析、水輪機獲能分析等功能。葉片的載荷分析如圖10所示。水輪機葉片載荷主要是流體動力作用在水輪機葉片上的載荷，包括葉片受到的彎矩、扭矩及其它作用于葉片上的載荷。

圖9 專用翼型數(shù)據(jù)選擇與數(shù)據(jù)顯示界面Fig.9 Implementation interface of hydrofoil data selection from Database

圖10 載荷分析界面Fig.10 Load analysis interface

5 結(jié)語

本文介紹了潮流能水輪機專用翼型族設(shè)計采用的多目標遺傳算法思路、翼型的參數(shù)化建模方法及翼型的設(shè)計實例。得到五種潮流能水輪機專用葉片翼型，最大相對厚度范圍為12.6%～38.8%，分別命名為OUC-TT-1260、OUC-TT-1900、OUC-TT-2150、OUC-TT-3240和OUC-TT-3880。優(yōu)化后的潮流能葉片翼型族具有較好的水動力特性。整理優(yōu)化后的潮流能水輪機專用翼型族數(shù)據(jù)，搭建了潮流能水輪機專用翼型族數(shù)據(jù)庫，實現(xiàn)了潮流能水輪機設(shè)計軟件與專用翼型族數(shù)據(jù)庫的連接及翼型數(shù)據(jù)選取功能，提高了水輪機葉片的設(shè)計效率。