小型垂直軸風力機葉片結冰風洞試驗與數值計算

李巖，劉欽東，王紹龍，馮放，田川公太朗

(1.東北農業大學工程學院，黑龍江哈爾濱150030;2.東北農業大學理學院，黑龍江哈爾濱150030; 3.日本鳥取大學地域學部，鳥取6808551)

小型垂直軸風力機葉片結冰風洞試驗與數值計算

李巖1，*，劉欽東1，王紹龍1，馮放2，田川公太朗3

(1.東北農業大學工程學院，黑龍江哈爾濱150030;2.東北農業大學理學院，黑龍江哈爾濱150030; 3.日本鳥取大學地域學部，鳥取6808551)

利用風洞試驗和數值模擬相結合的手段，研究了小型垂直軸風力機葉片在旋轉狀態下的結冰特性以及結冰后翼型與風力機的氣動特性變化，以期為建立較復雜的大型水平軸風力機葉片旋轉試驗系統和研究其結冰機理、防除冰技術提供參考。試驗在東北農業大學自行設計的利用自然低溫的結冰風洞中進行，獲得了采用NACA0018翼型的小型2葉片垂直軸風力機風輪在5種尖速比下的結冰分布:風力機葉片結冰遍布葉片整個表面，隨著尖速比的增大，結冰形狀出現不對稱性。同時，數值模擬結果表明:葉片結冰后，隨著尖速比的增加和結冰量的增多，升力系數降低阻力系數增大的趨勢明顯，風力機的功率系數也隨之下降。分析發現，葉片結冰導致不同旋轉角下葉片翼型周圍的壓力場和速度場發生了不同程度的變化，從而氣動特性發生變化，影響了風力機性能。

結冰;垂直軸風力機;旋轉葉片;風洞試驗;數值模擬

0 引言

近年來，寒冷氣候條件對風力機的影響受到了國際風能界的廣泛關注［1-6］，研究風力機葉片結冰特性與防除冰技術成為了熱點領域［7］。當前國內外風力機葉片結冰的主要研究方法是結冰風洞試驗與數值模擬計算。大部分研究主要是針對靜止葉片翼型的二維結冰和小型水平軸風力機靜態結冰［8-12］。對于大型風力機的旋轉葉片的結冰試驗研究還未見報道，主要原因:一是尚未建立合理有效的結冰實驗相似準則;二是冰風洞試驗成本較高使得這方面的研究開展受到一定的限制。

本研究對風力機結冰試驗研究進行了新的探索。首先嘗試了利用北方寒冷地區的自然低溫實現在常規風洞中進行結冰試驗。通過前期研究確認了該方法的可行性［13］，研究小型直線翼垂直軸風力機葉片在靜態［14-15］與旋轉狀態下的結冰特性。選取該種風力機為研究對象的主要原因有兩個，一是這種風力機是目前國內外新型小型風力機的代表，研究其結冰特性有一定的實用價值;二是該種風力機葉片簡單，沒有水平軸風力機葉片變截面等變化，旋轉狀態下的葉片結冰特性較易通過實驗獲得。因此，可先通過研究這種風力機葉片結冰的試驗方法與特性，為探索建立較復雜的大型水平軸風力機葉片旋轉試驗系統提供參考，為研究水平軸風力機旋轉葉片結冰機理及風力機防除冰技術提供借鑒。通過實驗測試了不同尖速比下葉片結冰分布，其后將結冰后的翼型進行數值模擬計算，分析葉片及風輪的氣動特性變化，并從流場改變的角度分析氣動特性變化原因，為葉片防除冰技術提供依據。

1 實驗裝置與方法

1.1 風力機模型

根據實驗室風洞條件，設計了一臺小型直線翼垂直軸風力機，如圖1所示。風輪直徑為440 mm。葉片采用了直線翼垂直軸風力機常用的NACA0018翼型，葉片個數為2，材質為玻璃鋼，長度為250 mm，弦長為125mm。

圖1 試驗風力機模型Fig.1 W ind turbine model for test

1.2 結冰試驗系統與測試方法

本研究試驗在東北農業大學自行設計的利用自然低溫的結冰風洞［13］中進行，結冰試驗段截面為0.6m×0.6m，風速范圍為1～15m/s。為模擬室外結冰環境，在吹出口處安裝了水霧噴射系統提供結冰條件。圖2為風洞試驗系統示意圖。

圖2 測試系統Fig.2 Test system

試驗條件如下:LWC(Liquid Water Content，ISO 12494中規定大氣中物體結冰時的水滴含量參數)為2.32 g/m3;MVD(Median Volume Droplet，ISO 12494中規定大氣中物體結冰時的水滴粒子大小的參數)約為40μm;溫度在－12℃至－7℃之間;風速為4m/ s;風力機旋轉尖速比為五種:0.2、0.4、0.6、0.8、1.0。在固定轉速下，利用安裝在風輪上方的高速攝影記錄旋轉狀態下的葉片表面結冰。

2 數值模擬方法

2.1 氣動參數

計算涉及的葉片及風力機氣動參數定義如下:

式中，ρ為空氣密度，1.225 kg/m3;A為物體相對來流的投影面積，m2;U為來流速度，m/s。

即葉片尖端線速度與風速之比。式中，V是旋轉的風力機風輪外徑切線速度，m/s;U是風進葉輪前的速度，m/s;ω為風力機旋轉角速度，rad;R為風輪半徑，m。

式中，M為力矩，N·m;A為物體相對來流的投影面積，m2。

其中p為壓力，Pa。

2.2 求解方程

因為測試的空氣的雷諾數和馬赫數都很小，所以可以被視為不可壓縮流體。數值模擬是基于二維不可壓縮流N-S方程，連續性方程和動量方程如下所示:

式中，xi、xj分別為方向矢量;ui、uj分別為速度矢量。

將實驗得到的二維葉片結冰形狀進行數字化處理后導入流場中，網格劃分如圖3所示。

圖3 模型網格Fig.3 M esh around test model

3 結果與分析

3.1 結冰試驗結果

圖4所示為尖速比在0.2、0.6、1.0時葉片在30 min過程中的結冰情況。由圖可知，直線翼垂直軸風力機葉片結冰是遍布葉片整個表面的，隨著時間的增加，葉片結冰逐漸增厚，且結冰在葉片表面分布較為均勻。隨著尖速比的增大，結冰形狀出現了一定的不對稱性，在葉片尾緣外側以及前緣內側結冰情況更加明顯。

圖4 不同尖速比下30m in內葉片表面結冰分布Fig.4 Ice accretion on blade surface during 30 m inutes at different tip speed ratios

3.2 數值模擬結果

將上述結冰后的翼型進行數字化處理，對結冰前后的葉片升阻力系數、風力機功率系數、力矩系數以及葉片周圍流場等進行了數值模擬計算。圖5和圖6給出了結冰葉片的升力和阻力系數的變化，圖7給出了風力機功率系數的變化情況。

由圖5和圖6可知，隨著時間的增加，結冰量越大的情況下，葉片的升力系數呈明顯下降趨勢，而阻力系數則明顯增大。隨著結冰量的增加，升力的改變有一定的波動性，但整體趨勢一致。而阻力系數的降低基本上是隨時間呈線性增加。由此可知，結冰后的葉片氣動特性整體下降。

圖5 結冰后葉片升力系數Fig.5 Lift coefficient of iced airfoil

從圖7可知，葉片結冰后風力機功率系數在各個尖速比下都不同程度下降。由于尖速比在0.2時的風力機功率系數較小，所以結冰后的功率系數變化不大。當尖速比在0.4至0.8時，在前5分鐘內隨著結冰的出現，使得風力機功率系數迅速下降，然后隨著結冰的增加，風力機的功率變化維持在一定的水平。而當在尖速比為1.0時，可以明顯的看出風力機的功率系數隨結冰量的增加而下降。結冰30 min的實驗模型與未結冰的實驗模型相比較，功率系數下降了37%。由此可知，葉片結冰后風力機將不能按照原來的設計輸出功率，甚至要停機。

圖6 結冰后葉片阻力系數Fig.6 Drag coefficient of iced airfoil

圖7 風力機功率系數Fig.7 Power coefficient of w ind turbine w ith iced blade

為了進一步分析葉片結冰后對風力機氣動特性影響的機理，計算了葉片結冰后風輪在旋轉一周內的力矩系數變化，如圖8所示。

由圖8可知在未結冰時，風輪的力矩系數基本上都為正值，只有在攻角為110°到155°之間為負值。然而，葉片結冰后，風輪旋轉一周內的負力矩范圍明顯增多，且隨著時間的增加和結冰量的增大，負力矩的范圍和程度均呈增大的趨勢。結冰30 min分鐘的風輪力矩系數的影響最大。

為了進一步分析影響機理，針對結冰30min葉片和無結冰葉片風力機，選取三個典型的方位角0°、70°、120°，計算了結冰前后葉片周圍的流場和壓力場，如圖9～圖11所示。

圖8 有無結冰風輪旋轉一周內力矩變化(λ=1.0)Fig.8 Cmof rotors w ith and w ithout icing on blade(λ=1.0)

圖9 結冰前后葉片周圍流場(α=0°)Fig.9 Flow field around blade w ith or w ithout icing(α=0°)

由圖9可知，由于上面的葉片前緣結冰的存在導致前緣附近形狀變化，使得背面出現了較大的流動分離，從而影響到整個葉片周圍的流場。這一變化直接導致該葉片周圍壓力場的變化，較大的流動分離使得背面壓力增大，腹面壓力降低，從而導致升力下降，阻力升高。而對于圖中下面葉片來說，由于尾緣在前，結冰對其產生的影響相對較小;同時，在180°時葉片的升阻力特性原本也相對較差，因此，結冰對該葉片的影響程度相對較小。因此，在這一狀態下，風力機的整體氣動特性較結冰前大幅下降。

由圖10可以看出，無論是速度場還是壓力場，葉片結冰前后的流場變化并不是十分顯著，這與圖8中給出的風輪力矩變化情況基本一致。主要的變化還是來自左上方葉片結冰所產生葉片周圍流場變化。可以看到在該葉片腹面，由于結冰導致分離區增大，上下面的壓力變化增大，影響了升力的發揮。而右下側葉片周圍的壓力場也發生了一定的變化，因此從總體上影響了風力機的氣動特性。

由圖11可知，結冰導致左下方葉片背面的壓力顯著增大，以及右上方葉片前緣附近的壓力增大。在結冰時，該狀態下的風力機力矩為負值，因此，結冰后導致負力矩進一步加大，使得風力機的氣動特性進一步惡化。

圖10 結冰前后葉片周圍流場(α=70°)Fig.10 Flow field around blade w ith or w ithout icing(α=70°)

圖11 結冰前后葉片周圍流場(α=120°)Fig.11 Flow field around blade w ith or w ithout icing(α=120°)

4 結論

本文通過試驗測試了不同尖速比下葉片結冰分布，并對結冰后的翼型進行了數值模擬計算，分析了葉片及風輪的氣動特性變化，得到了如下結論:

1)試驗結果表明，采用對稱翼型的直線翼垂直軸風力機旋轉葉片的結冰在尖速比較小時均勻分布在整個葉片表面，隨著尖速比增大，冰層逐漸增厚，由于旋轉效應使結冰出現一定的不對稱性，在葉片尾緣外側及前緣內側的結冰不對性較明顯。

2)計算結果表明，旋轉葉片結冰后的升力系數降低，阻力系數增大，風輪功率系數下降，這種趨勢隨尖速比的增加和結冰量增多而更加顯著。

3)結冰后葉片翼型周圍流場發生改變，其變化程度與風輪旋轉角相關，這是導致結冰后的葉片氣動特性變化與風力機性能降低的主要原因。

［1］Kimura S，Sato T，Tsuboi K，et al.Numerical simulation of ice accretion on a wind turbine blade［J］.Journal of the Japanese Society of Snow＆Ice，2006，68(5):393-407.

［2］Mróz A，Holnicki-Szulc J，K?rn?T.Mitigation of ice loading on offshore wind turbines:feasibility study of a semi-active solution［J］.Computers＆Structures，2008，86(s 3-5):217-226.

［3］Virk M S，Homola M C，Nicklasson P J.Relation between angle of attack and atmospheric ice accretion on large wind turbine's blade［J］.Wind Engineering，2010，(6):607-614.

［4］Kraj A G，Bibeau E L.Impact of mitigation strategies on icing accumulation rate for wind turbines in cold climates［R］.Wind Energy-Technological Advances，2007.

［5］Kraj A G，Bibeau E L.Measurement method and results of ice adhesion force on the curved surface of a wind turbine blade［J］.Renewable Energy，2010，35(4):741-746.

［6］Homola M C，Virk M S，Nicklasson P J，et al.Performance lossesdue to ice accretion for a 5MW wind turbine［J］.Renewable Energy，2012，15，379-389.

［7］Bose Neil.Icing on a small horizontal axis wind turbine-Part 1: Glaze ice profiles［J］.Journal of Wind Engineering and Industrial Aerodynamics，1992，45(1):75-85.

［8］Kraj A G，Bibeau E L.Phases of icing on wind turbine blades characterized by ice accumulation［J］.Renewable Energy，2010，35(5):966-972.

［9］Han Yiqiang，Jose Palacios，Sven Schmitz.Scaled ice accretion experiments on a rotating wind turbine blade［J］.Wind Engineering and Industrial Aerodynamics，2012，109(11):55-67.

［10］Yi Xian，Zhu Guoling.Computation of glaze ice accretion on airfoil［J］.Acta Aerodynamica Sinica，2004，22(4):490-493.(in Chinese)易賢，朱國林.考慮傳質傳熱效應的翼型積冰計算［J］.空氣動力學學報，2004，22(4):490-493.

［11］Zhu Chengxiang，Wang Long，Fu Bin.Numerical study of wind turbine blade airfoil ice accretion［J］.Acta Aerodynamica Sinica，2011，29(4):522-528.(in Chinese)朱程香，王瓏，付斌.風力機葉片翼型的結冰數值模擬研究［J］.空氣動力學學報，2011，29(4):522-528.

［12］Deng Xiaohu，Lu Xuxiang，Li Luping.Numerical simulation of airfoil ice accretion process on horizontal-axis wind turbine blade［J］.Energy Technology，2010，31(5):266-271.(in Chinese)鄧曉湖，盧緒祥，李錄平.水平軸風力機葉片翼型結冰的數值模擬［J］.能源技術，2010，31(5):266-271.

［13］Li Y，Wang S L，Zheng Y F，et al.Design of wind tunnel experiment system for wind turbine icing by using natural low temperature［J］.Journal of Experiments in Fluid Mechanics，2016，30(2):54-58，66.(in Chinese)李巖，王紹龍，鄭玉芳，等.利用自然低溫的風力機結冰風洞實驗系統設計［J］.實驗流體力學，2016，30(2):54-58，66.

［14］Li Yan，Chi Yuan，Feng Fang.Wind tunnel test on icing on blade used for vertical axis wind turbine［J］.Journal of Engineering Thermophysics，2012，33(11):1872-1875.(in Chinese)李巖，遲媛，馮放.垂直軸風力機葉片表面結冰的風洞試驗［J］.工程熱物理學報，2012，33(11):1872-1875.

［15］Li Yan，Tagawa Kotaro，Feng Fang，et al.A wind tunnel experimental study of icing on wind turbine blade airfoil［J］.Energy Conversion and Management，2014，85(9):591-595.

W ind tunnel test and numerical simulation on blade icing of small-scaled vertical axis w ind turbine

Li Yan1，*，Liu Qindong1，Wang Shaolong1，Feng Fang2，Tagawa Kotaro3
(1.Engineering College，Northeast Agricultural University，Harbin 150030，China;2.College of Science，Northeast Agricultural University，Harbin 150030，China;3.Regional Faulty，Tottori University，Tottori 6808551，Japan)

The icing characteristics of rotating blade of a small-scaled Vertical Axis Wind Turbine (VAWT)and the aerodynamics performance changes of the iced blade airfoil and rotor were researched by wind tunnel tests and numerical simulation in order to provide technology reference for establishing the more complex test system for large scale Horizontal Axis Wind Turbine with rotating blade and study the icing mechanism and anti-icing and de-icing technology.The experiments were performed in an icing wind tunnel designed by Northeast Agricultural University using natural low temperature in cold winter.The test model was a small-scaled VAWT with 2 blades with NACA0018 airfoil.The icing accretions on rotating blade surface under 5 tip speed ratios were tested and recorded.Wind turbine blades were frozen over the entire surface of the blade.Ice shape asymmetry can be found on the blade with the increasing of tip speed ratio.Furthermore，the lift coefficient of iced blade decreased and drag coefficient increased with the increasing of tip speed ratio and the amount of ice based on the results of numerical simulation.The power coefficients of the rotor were also decreased.It can be also found that the change degree of pressure field and velocity field around the iced blade depended on rotating angles，which was the main reason affecting the aerodynamic performance of blade and the turbine.

icing;vertical axis wind turbine;rotating blade;wind tunnel test;numerical simulation

V211.753

A

10.7638/kqdlxxb-2015.0214

0258-1825(2016)05-0568-05

2015-12-18;

2016-01-14

國家自然科學基金(51576037)

李巖*(1972-)，男，黑龍江賓縣人，教授，研究方向:可再生能源綜合利用.E-mail:liyanneau@163.com

李巖，劉欽東，王紹龍，等.小型垂直軸風力機葉片結冰風洞試驗與數值計算［J］.空氣動力學學報，2016，34(5):568-572，586.

10.7638/kqdlxxb-2015.0214 Li Y，Liu Q D，Wang S L，et al.Wind tunnel test and numerical simulation on blade icing of small-scaled vertical axis wind turbine［J］.Acta Aerodynamica Sinica，2016，34(5):568-572，586.