基于多島遺傳算法的垂直軸風機翼型優化設計

趙杰， 汪志成， 邵華梅

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

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

1 多島遺傳算法

1.1 多島遺傳算法基本原理

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

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

(1)

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

種群、選擇、交叉和變異操作與傳統遺傳算法一致，但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 多島遺傳算法參數

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

(1)子群規模數。

(2)總體規模數，總體規模數等于子群規模數與島數的乘積，大的群體規模能夠改經GA的搜索的質量。

(3)島的個數。

(4)進化代數。

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

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

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

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

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

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

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

表1 多島遺傳算法配置參數設置Table 1 MIGA configuration parameter settings

2 優化算法設計

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

2.1 目標函數

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

(2)

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

2.2 設計變量

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

2.3 約束條件

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

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

此外，優化翼型的升、阻力系數需要滿足以下條件。

(3)

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

2.4 優化設計程序

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

圖3 自動優化評估過程Fig.3 Automatic optimization evaluation process

3 優化結果分析

3.1 翼型幾何輪廓

以升阻比為目標函數的優化研究，如圖4所示，經過40代的計算得到了收斂的優化結果，設計變量的最終取值如表3所示。

圖4 翼型優化迭代過程Fig.4 Airfoil optimization iterative process

表3 優化翼型設計變量取值Table 3 Values of optimized airfoil design variables

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

圖5 翼型優化前后輪廓Fig.5 Airfoil optimized front and rear profile

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

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

3.2 翼型氣動參數

優化目標為翼型的最大升阻比，并對升、阻力系數設置了約束條件，根據優化結果，優化翼型與初始翼型的氣動數據如表4所示。

表4 翼型優化前后的氣動數據Table 4 Aerodynamic data before and after airfoil optimization

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

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

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

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

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

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

4 結論

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

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

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

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