典型風力機翼型結冰特性分析

程杰 來永斌 王龍 武浩 李亮

摘要：針對風力機葉片的結冰會嚴重影響風電機組的功率輸出、危害風電機組的安全運行的情況，采用計算流體力學方法研究了典型的風力機翼型的結冰特性。首先，對攻角分別為0°，4°，8°的NACA0012翼型水滴撞擊特性進行了計算，計算結果與文獻實驗結果吻合良好，表明了計算方法的正確性。在此基礎上，對不同風速和不同溫度下的典型風力機葉片表面的覆冰進行了模擬研究，得到了水收集系數和結冰厚度隨溫度和速度的變化趨勢。結果表明：風速增加導致葉片表面結冰厚度增加，結冰附著面積增加;溫度降低導致葉片表面結冰厚度增加;S801翼型的抗冰效果明顯大于S802翼型的抗冰效果。研究結果有利于在低溫環境中合理選擇翼型，為降低葉片結冰危害程度提供參考。

關鍵詞：發電工程;風力機;翼型結冰;水滴收集系數;計算流體力學

中圖分類號：V211.3文獻標志碼：A

CHENG Jie，LAI Yongbin，WANG Long，et al.Study on the icing characteristics of typical wind turbine[J].Journal of Hebei University of Science and Technology，2019，40（1）：9-14.Study on the icing characteristics of typical wind turbine

CHENG Jie1， LAI Yongbin1， WANG Long1， ?WU Hao1， LI Liang2

（1 School of Mechanical Engineering， Anhui University of Science and Technology， Huainan， Anhui 232001， China;2 School of Mechanics and Optoelectronic Physics， Anhui University of Science and Technology， ?Huainan， ?Anhui 232001， China）

Abstract：The icing of wind turbine blade will seriously affect the power output of wind turbine and endanger the safe operation of wind turbine. Computational fluid dynamics （CFD） is used to study the icing characteristics of typical wind wing. Firstly， the impact characteristics of water droplets on NACA0012 type turbine under different attack angle of 0°，4°，8°are calculated. The calculated results are in good agreement with the experimental results， which shows that the calculation method is correct， and， the ice coating on the blade surface of typical wind turbine with different wind speed and temperature was simulated， and the variation trend of water collection coefficient and ice thickness with temperature and velocity was obtained. The results show that： 1） with the increase of wind speed， the ice thickness of blade surface increases and the icing adhesion area increases; 2） with the decrease of temperature， the ice thickness of blade surface increases; and 3） the anti-ice effect of S801 airfoil is obviously greater than that of S802 airfoil. The research result helps logically select turbine under low temperature conditions， providing reference for lowering the blade icing hazard.

Keywords：project for power generation; wind turbine; turbine blade icing; water drop collection coefficient; computational fluid dynamics

風力發電是一種清潔無污染、可持續的新型能源，可用來緩解對常規能源的依賴。中國風電行業起步較晚，卻取得了驚人的成績，風電場遍布全國各地。風力機通常工作在野外，工作環境惡劣，尤其在中國北方地區，風電機組經常會在低溫環境下工作，并極有可能會遭遇結冰現象，增加了部件的不平衡性和疲勞程度。因此研究風力機翼型結冰特性，可為風力機防除冰提供數據支持，降低風力機在結冰條件下發生故障的可能性。

數值模擬方法是結冰研究最常用的方法，國外首先采用了數值模擬方法對結冰翼型進行了預測研究，國內數值模擬研究也在不斷發展。SHIN等[1]提供了一種翼型前緣結冰表面粗糙度的測量方法，證實了光滑和粗糙冰區的存在，并指出結冰表面粗糙度受溫度和空氣的液態含水量的影響較大，受風速的影響不大;SON等[2]對翼型結冰增長進行了定量分析，并開發出相關的翼型結冰數值模擬軟件;BOSE等[3-4]研究了小型水平軸風力機覆冰的過程， 對葉片表面的覆冰做了分析，給出了覆冰過程中風力機輸出功率的變化;FORTIN 等[5]利用CFD技術對風力機葉片旋轉過程中的覆冰過程進行了模擬，對覆冰后風力機葉片的性能進行了相應的分析。FU等[6]通過Fluent計算空氣和水的兩相流動，進行了霜冰的外形和質量積累的仿真計算。HOMOLA等[7]利用結冰軟件模擬翼型結冰，并通過Fluent軟件對比了結冰前后翼型氣動性能的變化。KRAJ等[8]通過實驗對渦輪機葉片結冰進行研究，指出葉片上冰的黏附力與附著面積與環境的溫度成比例關系。馬強等[9]采用Fluent中的歐拉多相流模型，通過添加質量傳遞以及對流換熱系數計算模型等，實現了對風力機葉片表面覆冰的二維數值模擬，得到了葉片表面覆冰狀況。孔滿昭等[10]應用fensapice結冰計算軟件，對NACA0012翼型進行了結冰過程的模擬研究，同時對翼型結冰前后的氣動特性作了比較和分析，發現翼型表面在結冰后，會導致翼型前緣氣流提前分離，并嚴重破壞其氣動性能。桑為民等[11]對于霜狀結冰， 在給定的結冰氣候條件下， 采用數值方法， 預測整個結冰過程，并分析不同形狀的冰型對翼型繞流及氣動特性的影響。邵曉海等[12]采用歐拉兩相流法研究了冷水滴撞擊翼型表面的收集特性，結果表明，平均水滴直徑越大，局部水滴收集系數和撞擊區域越大。李巖等[13]利用風洞試驗和數值模擬相結合的手段，研究了小型垂直軸風力機葉片在旋轉狀態下的結冰特性以及結冰后翼型氣動特性變化，指出葉片結冰導致不同旋轉角下葉片翼型周圍的壓力場和速度場發生了不同程度的變化，從而氣動特性發生變化，影響了風力機性能。

河北科技大學學報2019年第1期程杰，等：典型風力機翼型結冰特性分析筆者采用計算流體力學方法對風力機翼型表面進行結冰情況的模擬研究，分析溫度和來流速度對翼型表面冰層增長的影響，對比不同翼型在相同條件下的結冰情況。

1結冰計算的原理及過程

1.1空氣流場的計算

通過求解空氣流體的控制方程得到翼型的外繞流流場，采用無量綱化積分形式的N-S方程[14-15]：t∫∫ΩQdV+∫αΩF（Q）·ndS=∫αΩG（Q）·ndS，（1）式中：Ω是控制體;Ω是對應的邊界;Q是守恒變量;F是無黏通量;G是黏性通量;n是控制體外法線向量，其中：Q=ρ

1.2水滴運動方程與冰形計算

分別跟蹤每一個水滴的運動并做出如下假設：1）水滴在運動過程中，體積保持不變;2）水滴在撞擊過程中不破裂，不濺射;3）水滴初始速度與自由來流的速度相等，水滴在運動過程中不影響流場的性質。在計算水滴運動軌跡時，考慮作用在水滴上的重力、浮力和阻力的影響，則根據牛頓第二定律，水滴運動方程可寫作式（5）。Mdd2xddt2=12ρaAdCd|u-ud|（u-ud）+（ρd-ρa）Vdg。（5）上述各項含義可參見文獻＼[12＼]。

結冰過程采用Messinger 控制容積思想，結冰厚度計算可參見文獻＼[16＼]。

2計算結果和分析

2.1流場網格劃分

使用Ansys軟件劃分計算區域網格，對弦長為1 m的S801和S802翼型分別進行了網格的劃分，計算區域尺寸為30 m×18 m×18 m，其中入口邊界距離葉片10 m，出口邊界距離葉片20 m，上下邊界距離葉片9 m。網格總數約為50萬。網格的劃分結果如圖1所示。

2.2算例驗證

首先對NACA012翼型進行了水滴特性的計算，翼型弦長c=0533 4 m，來流馬赫數Ma=0302，攻角AOA分別選取0°，4°，8°，水滴的當量直徑MVD=20 μm，空氣中液態水含量為1 g/m3。圖2給出了不同攻角下的水滴表觀密度ρd的分布云圖。從圖2中可見，隨著攻角的增加，上翼型表面的陰影區區域逐漸增加，下翼型表面的陰影區域范圍逐漸減小。圖3給出了不同攻角下翼型的水收集系數的分布，橫坐標中的S為該單元的中心點到駐點的距離，本文將計算結果與文獻＼[17＼]結果進行比對。圖3中可見，水滴局部收集系數的計算結果與文獻結果吻合良好，驗證了本文計算方法的正確性。通過圖中不同攻角下的水收集系數的結果對比，可以發現最大水收集系數不受攻角變化的影響，均保持在0.7附近。攻角增大，最大水收集系數點將向右偏離駐點，實際情況則是沿著葉片上表面向上移動。

2.3S801和S802翼型計算結果及分析

在0°攻角時，對溫度為263.15 K，來流速度為10，15和20 m/s，以及來流速度為20 m/s、不同溫度兩種情況分別對S801和S802翼型進行模擬計算。模擬計算過程中其他計算參數如表1所示。

結冰時間/s弦長/m來流壓強/Pa空氣含水量/（g·m-3）水滴當量直徑/μm4201101 325120

2.3.1不同風速下的結冰對比

局部水收集系數對結冰計算至關重要，分別對不同風速下葉片局部的水收集系數進行了對比，圖4給出了S801和S802翼型在溫度為263.15 K，來流速度分別為10，15和20 m/s時，其他計算參數和表1相同條件下葉片的水收集系數變化趨勢。為方便比較，橫坐標采用翼型表面單元的中心點到駐點的距離與翼型弦長的比值，用“S/c”表示。從圖4中可以看出，隨著速度的增加水收集系數也在不斷增加，這與水滴計算的假設條件相關。在假設條件中，水滴不破裂，水滴不會濺射，隨著速度增大水滴的慣性也會增大，水滴撞擊翼型表面時的趨勢也會增加，最終導致收集系數增大[18]。因為只考慮0°攻角的情況，故水收集系數最大位置點均集中在駐點附近。另外，對比S801和S802翼型的水收集系數，可以看出S802翼型的水收集系數最大值在相同的條件下，均大于S801翼型的水收集系數的最大值，兩者水收集系數的最大值相差1.5倍左右。

圖4水收集系數

Fig.4Coefficient of local water collection

風速是影響葉片表面結冰的重要因素，圖5給出了在溫度為263.15 K時，來流速度分別為10，15和20 m/s，其他參數與表1相同的條件下的S801和S802翼型前緣表面的結冰厚度變化圖。從圖中可以看出，隨著速度不斷增加翼型的最大結冰厚度也在不斷增加，結冰最大位置點都在集中駐點附近。隨著風速的增大，單位時間內葉片表面收集到的水滴數就越大，從而導致結冰量增加。隨著風速增加，水滴在空氣剪切應力和壓力作用下沿翼型上下表面運動距離更遠，此外低溫導致水滴凝結成冰，從而冰形覆蓋范圍更廣[19]。對比不同翼型的最大結冰厚度可以明顯看出，在相同條件下S802翼型最大結冰厚度均比S801翼型的最大結冰厚度大， S801翼型前緣結冰附著面較大，S802翼型前緣結冰較集中，兩種翼型的最大結冰厚度相差1.5倍左右。

2.3.2不同溫度下的結冰對比

溫度是影響結冰厚度的重要因素，溫度的不同會對結冰速度產生影響。在來流風速為20 m/s，其他條件和表1相同，分別對S801和S802翼型在溫度不同的情況下進行了模擬研究，溫度條件設置為243.15，248.15，253.15，258.15，263.15和268.15 K。

圖6給出了S801和S802翼型葉片表面在不同溫度下結冰厚度的變化趨勢，從圖中可以看出，隨著溫度的降低翼型的表面最大結冰厚度在不斷增加，在263.15～273.15 K之間冰厚會變化比較明顯，計算條件溫度的變化使得結冰的類型發生了改變[20]，溫度的不斷降低使得結冰類型從明冰向霜冰發生轉變，冰面的粗糙程度不斷加大，最終導致了溫度在263.15和273.15 K時結冰厚度的劇烈變化。

筆者首先在算例部分給出了軟件計算與試驗對比結果，驗證了計算過程的正確性，再經過數值模擬方法對S801翼型葉片和S802翼型葉片進行模擬，計算了不同風速下兩種翼型水收集系數的變化趨勢，分析了在不同溫度和不同風速下的葉片前緣結冰厚度的變化，得到如下結論。

1）風速對風力機前緣葉片表面的水收集系數影響較大，隨著速度的增加，葉片表面水收集系數在不斷增大。對比兩種翼型的計算結果可以看出，S802翼型的水收集系數要大于S801翼型。

2）風速對風力機葉片前緣表面結冰厚度有較大的影響，風速的增加會導致葉片前緣結冰厚度增加，結冰附著面積加大，且前緣表面結冰厚度變化的趨勢和表面水收集系數變化趨勢一致，反映出水收集系數對葉片表面結冰的影響。

3）溫度是影響葉片表面結冰的另一個重要因素，溫度的降低會導致葉片結冰厚度的增加。

4）不同種翼型對結冰的抗性不同，在相同條件下S801翼型的抗冰性能要優于S802翼型的抗冰性能。

本文為簡化計算過程，采用二維截斷模型得到了溫度和來流風速對兩種翼型的影響結果，為解決工程實際問題提供了參考依據。但在計算過程中，水滴假設為球形，與實際情況可能存在一定偏差，對結果有一定的影響，后續工作尚需改進。

參考文獻/References：

[1]SHIN J. Characteristics of surface roughness associated with leading-edge ice accretion[J]. Journal of Aircraft， 2015， 33（2）：316-321.

[2]SON C， OH S， YEE K. Quantitative analysis of a two-dimensional ice accretion on airfoils[J]. Journal of Mechanical Science and Technology， 2012， 26（4）：1059-1071.

[3]BOSE N. Icing on a small horizontal-axis wind turbine part Ⅰ： Glaze-ice profiles[J]. Journal of Wind Engineering and Industrial Aerodynamics， 1992， 45（1）： 75-85.

[4]BOSE N. Icing on a small horizontal-axis wind tubine part Ⅱ： Three dimensional ice and wet snow formation[J].Journal of Wind Engineering and Industrial Aerodynamics，1992， 45（1）：87-96.

[5]FORTIN G， PERRON J. Wind Turbine Icing and De-Icing[R]. Florida： AIAA 2009-274，2009.

[6]FU P， FARZANEH M. A CFD approach for modeling the rime-ice accretion process on a horizontal-axis wind turbine[J]. Journal of Wind Engineering & Industrial Aerodynamics， 2010， 98（4）：181-188.

[7]HOMOLA M C， VIRK M S， WALLENIUS T， et al. Effect of atmospheric temperature and droplet size variation on ice accretion of wind turbine blades[J]. Journal of Wind Engineering & Industrial Aerodynamics， 2010， 98（12）：724-729.

[8]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.

[9]馬強， 吳曉敏， 陳彥. 風力機葉片覆冰的數值模擬[J]. 工程熱物理學報， 2014， 35（4）：770-773.

MA Qiang， WU Xiaomin， CHEN Yan. Numerical simulation of ice accretion on wind turbine blades[J]. Journal of Engineering Thermophysics， 2014， 35（4）：770-773.

[10]孔滿昭， 蒲宏斌. 翼型結冰過程數值模擬驗算與分析[J]. 航空工程進展， 2011， 2（3）：267-272.

KONG Manzhao， PU Hongbin. Numerical simulation check and analysis of ice accretion on airfoil[J]. Advances in Aeronautical Science and Engineering， 2011， 2（3）：267-272.

[11]桑為民， 蔣勝矩， 李鳳蔚. 翼型冰增長和結冰影響的數值模擬研究[J]. 應用力學學報， 2008， 25（3）：371-374.

SANG Weimin， JIANG Shengju， LI Fengwei. Numerical simulation for ice accretion and icing effects on airfoils[J]. Chinese Journal of Applied Mechanics， 2008， 25（3）：371-374.

[12]邵曉海， 劉星. 數值模擬過冷水滴撞擊翼型表面的收集特性[J]. 應用力學學報， 2015，32（2）：266-270.

SHAO Xiaohai， LIU Xing. Numerical method to simulate collection efficiency of supercooled droplets impinging on an airfoil[J]. Chinese Journal of Applied Mechanics， 2015，32（2）：266-270.

[13]李巖， 劉欽東， 王紹龍，等. 小型垂直軸風力機葉片結冰風洞試驗與數值計算[J]. 空氣動力學學報， 2016， 34（5）：568-572.

LI Yan， LIU Qiandong， WANG Shaolong， 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.

[14]李鑫， 白俊強， 王昆，等. 歐拉兩相流法在翼型積冰中的應用[J]. 中國科學：物理學 力學 天文學， 2014， 44（3）：258-266.

LI Xin， BAI Junqiang， WANG Kun， et al. The application of eulerian two-phase flow method in airfoil ice accretion[J]. Scientia Sinica Physica， Mechanica & Astronomica， 2014， 44（3）：258-266.

[15]朱程香， 王瓏， 孫志國，等. 風力機葉片翼型的結冰數值模擬研究[J]. 空氣動力學學報， 2011， 29（4）：522-528.

ZHU Chengxiang， WANG Long， SUN Zhiguo， et al. Numerical study of wind turbine blade airfoil ice accretion[J]. Acta Aerodynamica Sinica， 2011， 29（4）：522-528.

[16]GENT R W， DART N P， CANSDALE J T. Aircraft icing[J]. Philosophical Transactions Mathematical Physical & Engineering Sciences， 2000， 358：2873-2911.

[17]KIM J W， DENNIS P G， SANKAR L N， et al. Ice Accretion Modeling using an Eulerian Approach for Droplet Impingement[R]. Reston： AIAA 2013-0246，2013.

[18]魯天， 桑為民， 劉曉宇，等. 三維翼面結冰過程及其影響數值研究[J]. 航空動力學報， 2014， 29（6）：1339-1345.

LU Tian， SANG Weiming， LIU Xiaoyu， et al. Numerical simulation of icing accretion and influence on three-dimensional wing configuration[J]. Journal of Aerospace Power， 2014， 29（6）：1339-1345.

[19]陳曉圓. 飛機結冰多參數影響研究[D]. 南京：南京航空航天大學， 2015.

CHEN Xiaoyuan. Study of Multi Parameter Influence of Aircraft Icing[D]. Nanjing： Nanjing University of Aeronautics and Astronautics， 2015.

[20]任鵬飛. 結冰風力機葉片的空氣動力學特性數值研究[D]. 北京：中國科學院研究生院（工程熱物理研究所）， 2014.

REN Pengfei. Numerical Study on Aerodynamics of Wind Turbine Blades with Icing[D]. Beijing： Institute of Engineering Thermophysics Chinese Academy of Sciences， 2014.第40卷第1期河北科技大學學報Vol.40，No.1

2019年2月Journal of Hebei University of Science and TechnologyFeb. 2019