基于多島遺傳算法的垂直軸風機翼型優(yōu)化設(shè)計

趙杰， 汪志成， 邵華梅

(1.東華理工大學機械與電子工程學院， 南昌 330013； 2.江西省新能源工藝及裝備工程技術(shù)研究中心， 南昌 330013)

垂直軸風力機葉片翼型的優(yōu)化設(shè)計是改善風機氣動性能、提升風機風能利用率的有效手段，大致可以分為直接優(yōu)化法和反設(shè)計法兩種。直接設(shè)計法通常采用梯度優(yōu)化算法和全局搜索算法來正向求解空氣動力學問題；反設(shè)計法需要對翼型優(yōu)化后氣動參數(shù)建模求解和迭代計算，是一種逆向求解的方法[1]。賴怡等[2]基于模擬退火算法對翼型進行優(yōu)化設(shè)計。Wickramsinghe等[3]使用粒子群優(yōu)化算法對低速翼型進行研究，開發(fā)了一個新的流體評估方案。劉春等[4]對比幾種算法對NACA0012翼型氣動性進行分析，認為改進的Coupled算法對氣動性能預(yù)測有幫助。Chen等[5]在考慮了前緣粗糙度敏感性、失速性能等因素影響的情況下，基于粒子群算法構(gòu)建多目標優(yōu)化算法。Timnak等[6]基于改進的遺傳算法對翼型進行優(yōu)化設(shè)計，計算速度加快并且優(yōu)化效果更好。Le-Duc等[7]將計算流體動力學(computational fluid dynamics，CFD)方法與改進的遺傳算法聯(lián)合仿真，并在特定條件下對風力機翼型進行了優(yōu)化研究。余剛等[8]以遺傳算法和復(fù)合型法為基礎(chǔ)，開發(fā)并編寫流場計算程序，研發(fā)了更符合實際應(yīng)用的翼型優(yōu)化方法，無論是單點設(shè)計還是多點設(shè)計，該方法都能提高翼型的氣動性能。陳進等[9]構(gòu)建了多目標翼型優(yōu)化模型，并采用改進的粒子群算法對初始翼型進行優(yōu)化，得到的優(yōu)化翼型的升阻比和升力系數(shù)均有所提高。汪泉等[10]采用類函數(shù)與B樣條結(jié)合的翼型參數(shù)化表達方法，以切向力系數(shù)之和作為優(yōu)化的目標函數(shù)，并使用粒子群算法和氣動性能預(yù)測軟件RFOIL對翼型氣動外形進行優(yōu)化設(shè)計，分別從速度分布、功率系數(shù)、渦量分布和功率系數(shù)這4個方面對優(yōu)化翼型進行，結(jié)果顯示優(yōu)化后的翼型有效提高了力矩系數(shù)和功率系數(shù)。綜上所述，目前研究人員對垂直軸風力機翼型的氣動優(yōu)化進行了大量研究，但對垂直軸風力機翼型依靠CFD方法和全局搜索算法進行氣動性能優(yōu)化研究的分析還不夠深入。針對以上問題，現(xiàn)以MATLAB實現(xiàn)參數(shù)化程序設(shè)計，基于CFD方法中的網(wǎng)格自動生成和流程自動計算技術(shù)，結(jié)合多島遺傳算法構(gòu)建翼型目標函數(shù)自動優(yōu)化程序，以期實現(xiàn)基于CFD和遺傳算法的垂直軸風力機翼型優(yōu)化。

1 多島遺傳算法

1.1 多島遺傳算法基本原理

多島遺傳算法(multi-island genetic algorithm，MIGA)是基于遺傳算法改進的，比傳統(tǒng)的優(yōu)化算法具有更好的全局求解能力和計算效率。MIGA將整個進化種群再被劃分為若干個子種群，被劃分的子種群稱為“島嶼”，然后對每個島嶼上的子種群分別進行傳統(tǒng)遺傳算法的遺傳算子操作(選擇、交叉、變異)，如圖1所示。MIGA會每隔一段時間隨機選取若干個個體進行“遷移”操作，把這些選中的個體遷移到其他島嶼上，這種遷移操作有效地增強了種群的多樣性，從而防止“早熟”問題。此外，MIGA的計算速度高于傳統(tǒng)的遺傳算法，更加節(jié)省研究時間。

MIGA是從遺傳算法(genetic algorithm，GA)發(fā)展而來的，在參數(shù)以及算法的計算流程上基本一致，但MIGA在參數(shù)名稱和數(shù)目上略有區(qū)別。MIGA同樣采用是適應(yīng)度函數(shù)來評價個體的優(yōu)劣情況，適應(yīng)度值大的個體發(fā)生繁殖行為的概率較大，這樣保持了產(chǎn)生的新群體的平均適應(yīng)度值高于舊群體。在基因編碼方面，MIGA采用格雷碼編碼，表達式為

(1)

式(1)中：gm為m位的格雷碼;bm為m位的二進制碼。

種群、選擇、交叉和變異操作與傳統(tǒng)遺傳算法一致，但MIGA的子種群稱之為“島”，并且增加“遷移”操作。具體流程如圖2所示。

圖1 多島遺傳算法示意圖Fig.1 Schematic diagram of multi-island genetic algorithm

圖2 MIGA相鄰兩代進化過程Fig.2 The evolutionary process of the two adjacent generations of MIGA

1.2 多島遺傳算法參數(shù)

在多島遺傳算法中，有10個參數(shù)可以調(diào)節(jié)，其中3個為基本條件參數(shù)，另外7個為高級條件參數(shù)。對配置參數(shù)進行簡單說明如下。

(1)子群規(guī)模數(shù)。

(2)總體規(guī)模數(shù)，總體規(guī)模數(shù)等于子群規(guī)模數(shù)與島數(shù)的乘積，大的群體規(guī)模能夠改經(jīng)GA的搜索的質(zhì)量。

(3)島的個數(shù)。

(4)進化代數(shù)。

(5)交叉概率Pc，顧名思義即表示交叉算子的使用頻率。在每一代群體中，需要對Pc×n個個體的染色體結(jié)構(gòu)進行交叉操作。交叉概率與產(chǎn)生新一代的速度有關(guān)，通常取值較大，但Pc取值過大會破壞群體的優(yōu)良模式，過小會導(dǎo)致進化速度過慢。因此，一般取交叉概率Pc為0.6～1.0。

(6)變異概率Pm。變異概率是群體保持多樣性的保障，表示個體染色體上基因產(chǎn)生變異的概率，交叉算子結(jié)束后，處在繁衍階段的全部個體的每位等位基因以Pm隨機改變，故在每代當中大約會產(chǎn)生Pm×n×L次變異。變異概率同樣會影響新一代的產(chǎn)生速度和質(zhì)量，一般取Pm為0.005～0.01。

(7)島間遷移率。島與島之間種群交換所遵循的比率，一般取值在0～0.5。

(8)遷移間隔代數(shù)。即兩次遷移之間的間隔。

(9)競賽子群體比率。從父代的所有子種群中以一定比例隨機選擇個體進行比較，將種群當中最好的個體保留到下一點。

(10)精英個體數(shù)量。即通過競爭后，保留到下一代的最好個體的數(shù)量。

在優(yōu)化研究過程中，對多島遺傳算法的配置參數(shù)設(shè)計如表1所示，其他參數(shù)則按照默認設(shè)置。

表1 多島遺傳算法配置參數(shù)設(shè)置Table 1 MIGA configuration parameter settings

2 優(yōu)化算法設(shè)計

以NACA0012作為初始翼型，對風速10 m/s、攻角8°時的翼型進行優(yōu)化，找出滿足該條件下的最大升阻比的最佳翼型輪廓，并分析優(yōu)化翼型氣動性能。

2.1 目標函數(shù)

升阻比是風力機葉片性能的重要指標，并且高升阻比的翼型可以在低風速的地區(qū)，滿足風力機的低風速要求。在優(yōu)化研究中，以升阻比作為優(yōu)化設(shè)計的目標函數(shù)。翼型的優(yōu)化目標函數(shù)為

(2)

式(2)中：Cl、Cd分別為優(yōu)化翼型的升、阻力系數(shù)。

2.2 設(shè)計變量

采用Hicks-Henne參數(shù)法對翼型輪廓描述，并基于MATALB編寫翼型參數(shù)化計算程序。本優(yōu)化以控制翼型形狀的系數(shù)ci(i=0,1,…,13)為設(shè)計變量，通過多島遺傳算法改變設(shè)計變量的取值來控制翼型形狀的變化。

2.3 約束條件

在翼型的優(yōu)化約束中，一般分為兩種：幾何條件約束與氣動性能約束。前者主要是指翼型在參數(shù)化建模過程中可以修改翼型輪廓曲線的相關(guān)控制參數(shù)的取值范圍；后者則是指翼型的氣動性能相關(guān)參數(shù)，如升阻比、升力系數(shù)、阻力系數(shù)等，具體的約束需要更具翼型的設(shè)計要求來確定。優(yōu)化設(shè)計中主要約束包括：設(shè)計變量的取值范圍約束；初始翼型升、阻力系數(shù)的約束。優(yōu)化設(shè)計中的變量約束范圍如表2所示。

表2 設(shè)計變量ci取值范圍Table 2 Value range of design variable ci

此外，優(yōu)化翼型的升、阻力系數(shù)需要滿足以下條件。

(3)

式(3)中：Cl0和Cd0為原始升、阻力系數(shù)。經(jīng)數(shù)值模擬計算結(jié)果，Cl0=0.808 0，Cd0=0.021 8。

2.4 優(yōu)化設(shè)計程序

優(yōu)化設(shè)計程序是以多島遺傳算法作為核心算法來計算目標函數(shù)值的。如圖3所示，開發(fā)了用于優(yōu)化的自動目標函數(shù)評估流程，該流程將翼型參數(shù)化、網(wǎng)格生成和目標函數(shù)評估與CFD仿真相結(jié)合。采用宏命令，順序調(diào)用所有軟件。自動優(yōu)化過程開始于在多島遺傳算法中生成14個控制變量，即Hicks-Henne型函數(shù)的控制系數(shù)。然后基于該型函數(shù)的參數(shù)化，生成新的翼型輪廓數(shù)據(jù)點并保存到腳本文件中。利用ICEM的網(wǎng)格自動生成技術(shù)完成每一次修改后翼型形狀的網(wǎng)格生成，以便在Fluent中調(diào)用日志文件進行流場自動計算。接著對數(shù)值計算結(jié)果進行后處理以計算翼型的升、阻力系數(shù)及升阻比。最后，如果滿足全局最優(yōu)解，則評估過程終止；否則，遺傳算子會產(chǎn)生的一組新的變量重復(fù)這個過程，自動優(yōu)化重復(fù)執(zhí)行。

圖3 自動優(yōu)化評估過程Fig.3 Automatic optimization evaluation process

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

3.1 翼型幾何輪廓

以升阻比為目標函數(shù)的優(yōu)化研究，如圖4所示，經(jīng)過40代的計算得到了收斂的優(yōu)化結(jié)果，設(shè)計變量的最終取值如表3所示。

圖4 翼型優(yōu)化迭代過程Fig.4 Airfoil optimization iterative process

表3 優(yōu)化翼型設(shè)計變量取值Table 3 Values of optimized airfoil design variables

設(shè)計變量取值對應(yīng)的翼型輪廓即為優(yōu)化翼型，在本文中命名為“優(yōu)化翼型”。翼型優(yōu)化前后輪廓曲線對比如圖5所示，翼型輪廓的弦線位置用(x，y)表示，c為弦長，圖中橫坐標為翼型輪廓弦線橫坐標值與弦長的比值，縱坐標為翼型輪廓弦線縱坐標與弦長的比值。從圖5中可以看出，優(yōu)化翼型的中弧線輪廓曲線已經(jīng)是非對稱圖形，中弧線不再與弦長重合，優(yōu)化翼型有明顯的彎度變化。圖6中，將翼型與風輪旋轉(zhuǎn)軸直接連接的一側(cè)定為翼型下翼面，另一側(cè)則為翼型上翼面。相比初始翼型NACA0012，優(yōu)化翼型的上翼面厚度增加，下翼面厚度減少，且下翼面厚度減少幅度大于上翼面的增加幅度。根據(jù)中弧線可以看出，靠近優(yōu)化翼型前緣位置的下翼面厚度略大于上翼面厚度。

圖5 翼型優(yōu)化前后輪廓Fig.5 Airfoil optimized front and rear profile

如圖6所示，初始翼型的相對厚度大于優(yōu)化翼型，初始翼型的最大相對厚度為0.12c，其對應(yīng)弦線上的位置0.30c；優(yōu)化翼型的最大相對厚度為0.11c，其對應(yīng)弦線的位置為0.264c。與初始翼型相比，優(yōu)化翼型的最大相對厚度降低了0.01c，其對應(yīng)弦線的位置向翼型前緣方向移動了12%；優(yōu)化翼型的前緣部分與初始翼型基本一致，但后緣的夾角明顯減小了。比較而言，采用自適應(yīng)模擬退火算法對NACA4418 翼型優(yōu)化后，翼型厚度變化在5%以內(nèi)，最大厚度位于整個弦長的 0.24～0.35位置處[11]，基于CFD技術(shù)與遺傳算法對DU97-W-300翼型優(yōu)化后，相對厚度變化為1%，最大厚度的弦向位置由0.3c前移至0.294c。本文算法與這兩種算法比較，最大厚度對應(yīng)弦長位置均前移，其范圍基本一致，但本文算法對于翼型厚度的優(yōu)化更加明顯，針對NACA0012的翼型厚度變化達10%左右。

圖6 翼型優(yōu)化前后的相對厚度分布Fig.6 Relative thickness distribution before and after airfoil optimization

3.2 翼型氣動參數(shù)

優(yōu)化目標為翼型的最大升阻比，并對升、阻力系數(shù)設(shè)置了約束條件，根據(jù)優(yōu)化結(jié)果，優(yōu)化翼型與初始翼型的氣動數(shù)據(jù)如表4所示。

表4 翼型優(yōu)化前后的氣動數(shù)據(jù)Table 4 Aerodynamic data before and after airfoil optimization

與初始翼型相比，8°攻角下優(yōu)化翼型的升阻比和升力系數(shù)都得到了明顯提高，升阻比提高了17.78%，升力系數(shù)提高了16.7%。由于在優(yōu)化過程中對阻力系數(shù)進行了約束，因此阻力系數(shù)的變化可以忽略不計。比較而言，采用標準模擬退火算法對NACA4418 翼型優(yōu)化后，翼型阻力系數(shù)減小10.6%，升力系數(shù)提高1.45%，升阻比提高13.48%[11]，而本文算法通過抑制阻力系數(shù)，可獲得更高的升力系數(shù)和升阻比。

圖7、圖8分別給出了初始翼型和優(yōu)化翼型的壓力和速度云圖。從優(yōu)化前后的壓力云圖來看，相比初始翼型，優(yōu)化翼型的上翼面壓力增加，下翼面減少了，但上翼面的壓力分布更加均勻，并且尾緣位置的渦流現(xiàn)象降低了；從速度云圖來看，上翼面速度增加下翼面速度降低，尾緣附近速度有所增加。翼型的壓力與速度在優(yōu)化前后的變化可能是翼型的彎度增加導(dǎo)致的，翼型的彎度增加，氣體在前緣的流動變大，吸力增大；翼型最大厚度對應(yīng)在弦長位置上向前移動，前緣和尾緣處的氣流變細，從而導(dǎo)致氣體的流速加快，壓力變小，升力系數(shù)變大。在優(yōu)化中對阻力系數(shù)設(shè)置了約束條件，阻力系數(shù)幾乎不發(fā)生變化，因此升阻比增大。

圖7 翼型優(yōu)化前后壓力云圖Fig.7 Pressure cloud diagram before and after airfoil optimization

圖8 翼型優(yōu)化前后速度云圖Fig.8 Speed cloud diagram before and after airfoil optimization

另外分別計算了攻角為0°、6°、12°初始翼型與優(yōu)化翼型的壓力系數(shù)，計算結(jié)果如圖9所示。當攻角為0°[圖9(a)]時，已經(jīng)可以明顯看出兩者的差異了，因為初始翼型屬于對稱翼型，上下翼面的壓力分布對稱，此時無升力。但優(yōu)化翼型并非對稱形狀，因此上下翼面的壓力存在壓差，此時的優(yōu)化翼型具有升力。當攻角為6°[圖9(b)]時，初始翼型與優(yōu)化翼型的最大正壓力系數(shù)值均位于翼型的前緣處，值約為1，在弦長0.1c內(nèi)，初始翼型上下翼面壓力系數(shù)的絕對值大于優(yōu)化翼型，但在0.1c到1c的范圍內(nèi)情況則相反。當攻角為12°[圖9(c)]時，優(yōu)化翼型最大負壓力系數(shù)明顯高于初始翼型，從圖中可以看出初始翼型負壓力系數(shù)的最下端為-4，優(yōu)化翼型則為-2.6左右。

圖9 攻角為0°、6°、12°下的壓力系數(shù)變化情況Fig.9 Pressure coefficient changes at angles of attack of 0°, 6° and 12°

4 結(jié)論

首先對多島遺傳算法原理作了基本闡述，構(gòu)建了翼型優(yōu)化數(shù)學模型，建立基于多島遺傳算法與CFD耦合的垂直軸風力機翼型優(yōu)化程序。該程序主要依靠多島遺傳算法順序調(diào)用翼型參數(shù)化程序、ICEM網(wǎng)格劃分、Fluent流場計算，并對流場計算結(jié)果進行處理。以NACA0012翼型作為初始研究對象，使用該算法進行優(yōu)化算例分析，主要分析結(jié)論具體如下。

(1)相比初始翼型，優(yōu)化翼型已經(jīng)是非對稱翼型，其最大相對厚度降低了1%，為0.11c；其最大相對厚度對應(yīng)在弦上的位置向翼型的前緣移動了0.036c。并且優(yōu)化翼型的中弧線不再與弦線重合，翼型彎度明顯增加，中弧線偏移量最大值為0.017 4c。優(yōu)化翼型的前緣半徑變化不大，但翼型尾緣夾角明顯變小，翼型上翼面輪廓線相對平緩。

(2)優(yōu)化翼型的升阻比、升力系數(shù)均有提高，升阻比提高了17.78%，升力系數(shù)提高了16.7%。

(3)分析了攻角0°、6°、12°時翼型優(yōu)化前后的上下翼面壓力系數(shù)變化情況。發(fā)現(xiàn)當攻角為0°時，優(yōu)化翼型相比對稱翼型已有升力，到攻角為12°時，優(yōu)化翼型的負壓力系數(shù)明顯高于初始翼型。