新型雙風輪風力機氣動特性的三維流場數值模擬

周云龍， 楊承志， 李律萬

（1.東北電力大學 能源與動力工程學院，吉林132012；2.中科院廣州能源研究所，廣州510640）

風力發電在潔凈可再生能源的開發利用中占有 重要的地位，風輪是風力發電機捕捉風能的關鍵部件，它的設計直接決定了風力機的出力和風能轉換效率，因此風輪的研究是風力機整機設計的重點.傳統的單風輪風力機葉片少則迎風面積小，葉片多則轉速低，都不能獲得較高的輸出功率，而新型的雙風輪風力機（見圖1）則彌補了單風輪風力機的缺陷，在發電機的兩側裝設前后2個風輪，后風輪捕捉前風輪的漏風，增大了風力機的捕風面積，并且各自保持了較高的轉速，同時應用了新型、高效的對轉式異步發電機［1］，前、后2個風輪分別帶動發電機的內、外轉子進行相對旋轉運動，增大了線圈切割磁感線的速度，從而獲得了較大的發電功率.

圖1 新型雙風輪風力機結構示意圖Fig.1 Structural diagram of a new double－rotor wind turbine

然而，目前雙風輪風力機作為一種新型的風力機，各國對其的研究還處于起步階段.日本的Kubo［2］以及國內的安源等［3］對2個風輪同側布置的雙風輪風力機已經有了初步的研究成果，但是這種布置方式在實際運行中存在缺陷，例如轉動部件摩擦嚴重、整個風力機維持受力平衡困難以及2個風輪由于間距小、可調性差而造成碰撞風險等，影響了風力機的運行安全和使用壽命.2個風輪異側布置的新型雙風輪風力機可以有效地解決同側布置存在的問題，并且可以同樣實現增大發電功率的目的.不過，這方面的研究剛開始，其流動特性尚不被人們熟知.近些年，隨著空氣動力學和CFD數值模擬理論的發展，Fluent軟件被廣泛應用于風力機流場的流動特性分析等領域，并得到了充分的驗證［4－7］.筆者利用Fluent軟件對這種2個風輪異側布置的新型雙風輪風力機的流場進行探索性的研究，并結合實驗測試結果加以分析.

1 風力機模型的建立

所研究的風力機為小型風力機，葉徑R為1.5 m，以吉林某地10年的風速數據為依據，確定風力機的額定風速為11.26m／s，選用NACA4412翼型數據，根據各種風輪的葉片數確定相應葉片的葉尖速比，利用Wilson葉片優化設計方法，基于Matlab自主編程完成風力機葉片的設計.風力機的塔架高度為3.3m，雙風輪風力機的前、后風輪間距為0.4 m，利用 Pro／E建立各類型風力機模型［8－9］，模型如圖2所示.

圖2 各類型風力機模型Fig.2 Models of different wind turbines

根據模擬的需要，構建了5種類型的風力機：3葉片單風輪風力機，記為N3；前輪3葉片、后輪3葉片雙風輪風力機，記為N前3，后3；前輪3葉片、后輪4葉片雙風輪風力機，記為N前3，后4；前輪3葉片、后輪5葉片雙風輪風力機，記為N前3，后5；前輪3葉片、后輪6葉片雙風輪風力機，記為N前3，后6.

2 數值計算方法

2.1 計算域和計算網格劃分

將風力機模型導入前處理軟件Gambit，并對流場進行網格劃分，為了能完整地查看風力機流場的氣動性能，計算域的橫向上風向長為10 R，下風向長為20 R，縱向關于風力機中心面對稱，寬為10 R.由于研究對象是扭曲葉片數目較多的風力機，物理模型比較復雜，盡量減少計算網格數有利于數值計算.在保證網格質量的前提下，在對計算域進行網格劃分時采取如下處理：一方面在不同的區域采用不同的網格生成方法，首先利用Size Function合理分布風力機各邊和面的網格，然后利用TGrid方法在風輪的旋轉區域以及風力機附近的小區域生成加密網格，并利用Cooper方法在風力機前后兩邊的大區域生成簡單網格；另一方面設置較大的網格節距，網格節距的增大在一定程度上增加了數值模擬結果曲線的拐點數目，大范圍的使用可能會降低曲線的平滑性，但保持合理的限度并不會影響數值計算的收斂性和結果的正確性.網格劃分結果如圖3所示，各類型風力機流場的總網格數均為350萬左右，Equiangle skew小于0.89，滿足了網格數量和質量的要求.

圖3 各類型風力機模型的網格劃分Fig.3 Grid division for models of different wind turbines

2.2 邊界條件

對各類型風力機在額定工況（風速v取額定風速11.26m／s）和5個工況（v分別取2m／s、5m／s、8m／s、10m／s和12m／s）下的運行情況進行數值模擬.模型入口為VELOCITY＿INLET，出口為OUTFLOW.風輪旋轉流體選用MRF模型，單風輪風力機的風輪附近為逆時針方向的旋轉域，雙風輪風力機的前、后風輪附近分別為逆時針和順時針方向的旋轉域.假設后風輪旋轉受前風輪的影響較小，根據式（1），可求得不同類型的風輪在不同風速下的轉速［10］，在Fluent軟件中進行旋轉域旋轉速度設定

式中：n為風輪轉速，r／min；λ為葉尖速比，3、4、5、6葉片風輪的λ一般分別取6、4.5、4、3；v為風輪上風向大氣流速，m／s；R 為風輪半徑，取1.5m.

2.3 求解方法

假定無傳熱現象，模擬采用Segregated（分離式）求解器隱式求解三維穩態不可壓縮流動，紊流模型使用SST k－ω 模型［11］，壓力－速度耦合采用Simplic算法，對流項差分格式采用二階迎風格式，參考壓力為1.01×105Pa.求解器的時間步長設為“自動調整時間步長”，最大迭代次數為2 000步，殘差類型為均方根，前后計算的殘差余量設為1×10－5.

2.4 模擬結果處理方法

風力機性能分析的主要參數包括風力機風輪輸出的機械功率P、風力機發電功率P0和風能利用系數CP，其計算公式［10－12］（對于單風輪風力機，可以認為是只有前風輪的雙風輪風力機進行計算）分別為

式中：M前、M后分別為前、后風輪的轉矩，N·m；n前、n后分別為前、后風輪的轉速，r／min；η1為發電機效率，%；η2為傳動效率，%，小型實驗性風力機一般取η1η2＝72%；ρ為通常狀況下的大氣密度，取值為1.205kg／m3；A 為風輪的掃掠面積，m2.

由式（4）可推導出雙風輪風力機較單風輪風力機的發電功率增加百分比αPo與風能利用系數增加百分比αCp的表達式為

式中：Po，單、Po，雙分別為單風輪、雙風輪風力機的發電功率，W.

3 模擬結果及分析

3.1 靜壓與速度云圖分析

在額定工況下，風力機風輪最有代表性的位置0.85 R處，5種類型風力機環形截面上的靜壓（相對壓力）云圖與速度（絕對速度）云圖見圖4.

對比圖4左側各類型風力機的靜壓云圖可知，大氣繞流各個葉片翼型時，在翼型的上表面附近大氣流動受阻，動能轉化為靜壓能，形成表壓為700 Pa左右的高壓區；在翼型的下表面邊界層發生分離，附近形成了表壓為－50Pa左右的低壓區，翼型這2個表面間的壓力差提供了葉片持續快速旋轉的升力.這樣，前風輪葉片后的低壓區與后風輪葉片前的高壓區綜合作用形成了2個風輪間相互影響的流場.當雙風輪風力機后風輪的葉片數目少時，這種影響作用較小，隨著后風輪葉片數目的增加，這種影響作用越明顯.

比較圖4右側各類型風力機的速度云圖可知，各風力機風輪旋轉區域入口處的大氣流速接近額定風速11.26m／s，而出口出現了風速為3m／s左右的低風速區域.與單風輪風力機相比，隨著雙風輪風力機后風輪葉片數目的增加，低風速區域增大，風速降低，并且無傳熱現象，這表明大氣通過對各個葉片的做功推動了風輪的旋轉，動能發生了形式上的轉化，雙風輪風力機轉化的能量較為明顯，從而使輸出較多的電能成為可能.

3.2 等速線圖分析

在額定工況下，5種類型風力機在沿大氣來流方向（x軸方向）與前風輪某葉片展向（葉片的葉根軸線方向）所形成的平面截面、沿大氣來流方向與前風輪該葉片展向垂直截面上的等速線圖見圖5.

圖4 各類型風力機的靜壓云圖（左）與速度云圖（右）Fig.4 Static pressure（left）and speed（right）contour of different wind turbines

圖5 各類型風力機的等速線圖Fig.5 Equal speed lines of different wind turbines

受風輪旋轉的影響，各類型風力機都在風輪前部出現了空氣壓縮流場（阻塞效應），但相差不大，可見雙風輪風力機雖然增加了后風輪，但只要葉片數目不是很多，其對前風輪流場的影響就不會很大，這對保持前風輪較大的輸出功率是非常重要的.各類型風力機的旋轉風輪都在葉片徑向外側的后部形成了漩渦，雙風輪風力機受到后風輪反向旋轉的影響，這部分渦流的大小比單風輪風力機的渦流有所減小，同時位置向下游延遲，這一區域為高速等值線區，風能較大，由于風力機葉片的出力主要是由這部分葉片貢獻，渦流的減小和延遲非常有利于后風輪捕捉前風輪的漏風.最后，大氣經過各類型雙風輪風力機的前風輪葉片時都會產生渦旋回流，并對后風輪造成一定程度的影響，在某些區域這種影響對后風輪的旋轉是有利的，在另外一些區域則不利.只有通過進一步的性能參數分析，確定后風輪葉片的合理數目，才能盡可能地增加對有利區域流場的利用，從而合理布置后風輪，使得后風輪可以較好地捕捉到前風輪的漏風，因而在增大迎風面積的同時，保持2個風輪較高的轉速，高效地實現風能的兩級利用.

3.3 湍動能云圖分析

在額定工況下，5種類型風力機流場的湍動能云圖如圖6所示.圖6中單風輪風力機的下風向流場的湍動能很快趨近于來流風速流場的湍動能E＝2.0m2／s2，而各類型雙風輪風力機的下游流場出現了湍動能為5.0m2／s2左右的高湍動能區域，且隨后風輪葉片數目的增加，該區域不斷擴大，N前3，后6類型雙風輪風力機下游流場部分區域的湍動能甚至達到了8.0m2／s2左右.由此可見，雙風輪風力機較單風輪風力機的湍動能有明顯地增大，且后風輪葉片數越多，流場的湍動能越大，相應的湍流強度也越大，N前3，后6類型雙風輪風力機的湍動能已經非常大了.所以，為了保證風力機的安全穩定運行，雙風輪風力機在設計加工時，后風輪葉片數目不能過多，而且在設計加工過程中應提高設備的機械強度.

圖6 各類型風力機的湍動能云圖Fig.6 Turbulent kinetic energy contours of different wind turbines

3.4 性能參數分析

通過數值計算，得到了各類型風力機風輪在不同工況下的輸出轉矩M，如表1所示.由表1可知，各類型的雙風輪風力機由于后風輪的影響，前風輪的轉矩比N3類型單風輪風力機風輪的轉矩小.但是，當后風輪葉片布置合理時，在大氣經過前風輪后的渦旋回流作用下，后風輪對前風輪漏風的利用能力有了較大的提高，其輸出轉矩甚至超過前風輪.前、后風輪的共同作用使雙風輪風力機總的輸出轉矩有所增加，進而在很大程度上提高了發電功率和增大了風能利用系數.

表1 各類型風力機風輪在不同工況下的輸出轉矩Tab.1 Output torque Mof wind rotors for different wind turbines under varying working conditions

將相應工況下各類型風力機風輪的轉速n分別代入式（2）～式（5），可得到各工況下評價各類型風力機的性能分析參數，即發電功率Po、風能利用系數CP和雙風輪較單風輪風力機的發電功率增加百分比及風能利用系數增加百分比，然后在Matlab中進行數據擬合得到各個參數的曲線（見圖7）.

由圖7可知，數值模擬得出的 N前3，后3、N前3，后4、N前3，后5與 N前3，后6類型雙風輪風力機的性能分析參數較N3類型單風輪風力機的性能分析參數均有一定的提高，其風能利用系數分別為29.4%、32.0%、34.3%和27.7%左右，相對于單風輪風力機分別增大了43.7%、56.3%、67.4%和35.3%，這與文獻［13］中的實驗測試結果一致，說明了模擬結果的正確性.

由數值模擬結果可知，由于后風輪葉片數目較多，N前3，后6類型雙風輪風力機流場的湍動能很大，使得運行穩定性較低，另外，受到前風輪的不利影響增加，其發電功率與風能利用系數相對于其他類型的雙風輪風力機反而有減小的趨勢，不適合實際應用，所以后風輪葉片數目應以5葉片以下為宜.其中，N前3，后5類型雙風輪風力機的流場相對平穩，發電功率及風能利用系數較大，可以作為主要研究方向.

圖7 模擬得出的各類型風力機的性能分析參數Fig.7 Performance parameters of various simulated wind turbines

4 結 論

（1）對比單風輪風力機的流場，當新型雙風輪風力機的后風輪葉片數目合理時，后風輪對前風輪的影響較小，而且能夠充分地利用大氣經過前風輪所產生的渦旋回流的有利影響，有效地捕捉到前風輪的漏風，這樣可以在增大迎風面積的同時，保持2個風輪較高的轉速，進而達到提高發電功率和增大風能利用系數的目的.當雙風輪風力機的后風輪在5葉片以內時，葉片數目越多，發電功率和風能利用系數的增加越明顯.

（2）雙風輪風力機比單風輪風力機的流場更加復雜，隨著后風輪葉片數目的增加，湍流強度變大，風力機運行穩定性變差.基于運行穩定性和高效性，N前3，后5類型雙風輪風力機將成為主要研究方向.

（3）數值模擬的性能參數變化趨勢與實驗測試的結果一致，說明與單風輪風力機相比，新型的雙風輪風力機在一定程度上提高了發電功率和增大了風能利用系數.

［1］BOOKER J D，MELLOR P H，WROBEL R，et al.A compact，high efficiency contra－rotating generator suitable for wind turbines in the urban environment［J］.Renewable Energy，2010，35（9）：2027－2033.

［2］KUBO K，KANEMOTO T.Development of intelligent wind turbine unite with tandem wind rotors and double rotational armatures［J］.Journal of Fluid Science and Technology，2008，3（3）：370－378.

［3］安源，韓鳳琴，久保田喬，等.反轉葉輪風力機非定常尾流對下游環境的影響［J］.水電能源科學，2010，28（10）：158－160.AN Yuan，HAN Fengqin，TAKASHI Kubota，et al.Effect of unsteady wake from counter－rotating wind rotor on downstream environment［J］.Water Resources and Power，2010，28（10）：158－160.

［4］BARTHELEMIE R J，RATHMANN O，FRANDSEN S T，et al.Modelling and measurements of wakes in large wind farms［J］.Journal of Physics：Conference Series，2007，75（3）：12－19.

［5］DUQUETTE M M，VISSER K D.Numerical implications of solidity and blade number on rotor performance of horizontal－axis wind turbines［J］.Journal of Solar Energy Engineering，2003，125（4）：425－432.

［6］李宇紅，張慶麟.風力機風輪三維流動特性與氣動性能的數值分析［J］.太陽能學報，2008，29（9）：1172－1174.LI Yuhong，ZHANG Qinglin.Numerical simulation of flow field and aerodynamic performance of a wind turbine blade［J］.Acta Energiae Solaris Sinica，2008，29（9）：1172－1174.

［7］祝賀，徐建源，騰云，等.風力機風輪氣動性能三維流場數值模擬［J］.中國電機工程學報，2010，30（17）：85－90.ZHU He，XU Jianyuan，TENG Yun，et al.3Dflow field numerical aerodynamic performance test of wind turbine rotor［J］.Proceedings of the CSEE，2010，30（17）：85－90.

［8］張果宇，蔣勁，劉長陸.風力發電機整機氣動性能數值模擬計算與仿真研究［J］.華東電力，2009，37（3）：449－451.ZHANG Guoyu，JIANG Jin，LIU Changlu.Numerical simulation of aerodynamic performance for wind turines［J］.East China Electric Power，2009，37（3）：449－451.

［9］李國寧，楊福贈，杜白石，等.基于 MATLAB與Pro／E的風力機風輪設計及造型［J］.機械設計，2009，26（6）：3－7.LI Guoning，YANG Fuzeng，DU Baishi，et al.Design and modeling of wind wheel of wind mill based on MATLAB and Pro／E［J］.Journal of Machine Design，2009，26（6）：3－7.

［10］何顯富，盧霞，楊躍進，等.風力機設計、制造與運行［M］.北京：化學工業出版社，2009：73－81.

［11］郭婷婷，吳殿文，王成蔭，等.風力發電機葉片預彎設計及其數值研究［J］.動力工程學報，2010，30（6）：450－455.GUO Tingting，WU Dianwen，WANG Chengyin，et al.Pre－bend design and numerical simulation of wind turbo－generator［J］.Journal of Chinese Society of Power Engineering，2010，30（6）：450－455.

［12］李少華，岳巍澎，匡青峰，等.雙機組風力機尾流互擾及陣列的數值模擬［J］.中國電機工程學報，2011，31（5）：101－107.LI Shaohua，YUE Weipeng，KANG Qingfeng，et al.Numerical simulation of wake interaction and array of double wind turbine［J］.Proceedings of the CSEE，2011，31（5）：101－107.

［13］周云龍，楊承志，岳巍澎.一種新型雙風輪風力發電裝置的特性分析［J］.中國電力，2011，44（4）：75－78.ZHOU Yunlong，YANG Chengzhi，YUE Weipeng.Characteristic analysis for a new type of double－rotor wind turbine［J］.Electric Power，2011，44（4）：75－78.