水平軸風力機葉片翼型光冰結冰的數值研究

顧聲龍，吳玉帥，袁曉偉，鄭仙佩，解宏偉

(青海大學 水利電力學院，青海 西寧810016)

青藏高原風能資源十分豐富，風能發展在青藏高原上如火如荼的進行著.青海選定的風力廠都處于高寒地區，風力機葉片的覆冰現象不可避免，這將導致葉片的氣動性能以及風力機功率的輸出達不到設計要求.風力機葉片的結冰是由過冷水滴撞擊葉片表面后凍結成冰覆蓋在葉片表面上.根據過冷水滴凍結的過程，將結冰類型分成霜冰(rime ice)和光冰(glaze ice).霜冰是指撞擊到葉片表面的過冷水滴立即凍結所形成的比較規則的冰;而光冰是指過冷水滴不能立即完全凍結，未凍結的水被風吹動沿表面流動，并在流動過程中逐漸凍結而形成的冰.

葉片結冰影響風力機的氣動性能和運行安全，給風力機帶來了極大危害.例如，覆冰會改變葉片的外部形狀和氣動性能，增大葉片阻力、減少升力，影響全機操縱性、穩定性，最終導致風能的轉化效率降低，嚴重的甚至會導致葉片的損毀，發生運行事故.風力機葉片結冰已經成為保證風力機安全運行需要解決的迫切問題.

對于葉片翼型結冰的研究，最理想的方法是風洞實驗.如Bose等人在1992年對處于寒冷氣候中的雙葉片水平軸風力機結冰進行了測量，但是并未指出結冰形狀和外部結冰條件的關系.Jasinski等人在1998年對NREL－S809翼型進行了冰風洞實驗，得出了該翼型結冰前后的性能參數［1］，但其實驗方法非常麻煩且相當昂貴.目前數值模擬已經應用到了風力機葉片結冰模擬中，但大多都是借鑒于對飛機機翼和輸電電纜結冰的相關研究經驗［1－6］.近幾年，國內對風力機葉片翼型的研究工作已逐步開展起來［7－9］.

文中將對在不同環境參數下，光冰結冰的形狀、結冰過程及結冰后葉片翼型的氣動性能進行研究.

1 求解方法

葉片翼型結冰的數值模擬一般采用基于拉格朗日方法的離散相模型DPM(discrete phase model).其過程為:1)計算連續相(空氣場)，即得到葉片翼型繞流的流場;2)基于連續相的流場進一步求解水滴的運動方程，從而確定水滴與物體的碰撞效率;3)按照一定的增長模型來確定與物體相碰撞的水滴的結冰以及冰的增長量.

對于連續相和離散相的求解，文中利用FLUENT軟件實現，利用FLUENT用戶自定義函數UDF自主開發出進行碰撞效率以及結冰增長的程序.

1.1 連續相的計算

對于風力機葉片翼型繞流來說，其流動雷諾數在百萬量級，屬于湍流.本項目采用DES來處理湍流粘性，求解雷諾平均方程，獲得空氣流場分布.而DES中采用Spalart－Allmaras(S－A)RANS模型，S－A模型是基于修正湍流運動粘滯率ν~的輸運方程，即

式中:ρ為密度;uj為j方向的速度分量;xj為j方向的曲線坐標;t為時間;Gv是湍流粘度生成項;Yν是在近壁區湍流粘度消散項;Sν~是用戶定義項，具體表達式可以參看 FLUENT User’s Guide.

1.2 離散相的計算

對葉片翼型結冰數值模擬時，必須對每個水滴進行跟蹤，以確定水滴是否與翼型發生碰撞.在計算水滴軌跡時需做如下假設:1)水滴均勻分布，以球形存在，運動過程中水滴尺寸保持不變;2)水滴的初始速度與自由流的速度相等，水滴體積很小以至于它們的繞流不會影響流場的性質;3)懸浮在運動空氣中水滴的附加質量力、壓差力、Saffman升力等與氣體阻力和重力相比可忽略.根據牛頓第二定律，水滴運動方程為:

式中:md為水滴質量;ρ為空氣密度;Cd為阻力系數;u為當地氣流速度;ud為水滴速度;Ad為水滴的迎風面積.對上式積分即可獲得水滴運動軌跡.

1.3 光冰結冰量的計算

光冰結冰量的計算方法與霜冰不同，霜冰結冰認為水滴立即結冰，不考慮熱力學變化，但光冰的計算需考慮水滴的熱力學變化(即部分結冰、部分液態、部分蒸發等現像)來計算控制體內的結冰厚度.

文中采用控制體內質量守恒和能量守恒原理［10］.

1)質量守恒

2)能量守恒

文中采用Messinger積冰傳熱模型，即

將式(1，2)聯立，即可求得結冰量，從而得到結冰厚度，即

文中利用FLUENT UDF來編寫水滴碰撞效率計算程序及結冰計算程序.其中結冰形狀的計算方法采用楊勝華等人關于機翼結冰過程模擬的方法——角平分線增長法［2］.

由于缺乏風力機葉片翼型結冰實驗數據，所以文中采用鄧曉湖［1］論文中所引用的Bose論文中的工況進行計算，并驗證本程序的正確性，如圖1所示.圖中C為弦長，單位為m(以下同)點劃線為文中計算結果，虛線為文獻數據，可以看到文中結果與文獻結果大致符合，由此說明了文中計算方法的正確性.

圖1 計算程序的驗證Fig.1 Verification of the procedure

2 結果分析

文中研究對象為NURLS－S809翼型.風力機葉片上距離翼根85%～100% 之間的區域對于風力發電而言最重要［8］，文中采用位于風力機葉片上距離翼根90%的位置的翼型作為研究對象，即葉片翼型弦長為1.52m.根據750kW風力機額定風速范圍以及結構的大氣結冰(ISO12494:2012)中關于大氣結冰的條件，設定光冰結冰的環境參數為風速10m/s，LWC(liquid water content)1g/m3，環境溫度－5℃，文中選取水滴直徑為30μm，來流攻角為0°、4°、6°和8°.

2.1 結冰過程及形狀

圖2給出了攻角為0°和6°，水滴直徑為30 μm，風速為10m/s，環境溫度為－5℃時葉片翼型結冰的形狀及過程.

圖2 光冰結冰形狀及過程Fig.2 Iced shape for glaze Ice

從圖2可以看到，葉片翼型光冰結冰的形狀在結冰初期是比較規則的，其增長方向基本上沿一個固定方向，即迎著空氣來流方向.隨著結冰過程的繼續，結冰的形狀逐漸出現棱角，變得不規則起來.這種冰型會破壞翼型自身的流線型，盡管只有前緣結冰，但會嚴重影響氣流的運動，使邊界層分離，流動不再穩定，影響氣動性能.

2.2 結冰翼型的氣動性能分析

2.2.1 流場分析

圖3為風速10m/s，LWC=1.0g/m3，水滴直徑30 μm，攻角分別為 0°、4°、6°和 8°來流條件下，潔凈翼型和結冰翼型的速度場以及翼型周邊的流線.從圖中可發現，隨著攻角增大，潔凈翼型表面會出現邊界層分離，而結冰后，這種分離區尺度加大，位置提前，尤其是在攻角0°下，由于冰型分叉，對翼型繞流非常不利，導致在翼型尾緣產生相當大的分離區，這將直接影響翼型的氣動性能.

圖3 攻角為0°，4°，6°，8°的潔凈翼型和結冰翼型流場比較Fig.3 Flow field of clean airfoil and iced airfoil at AOA 0°，4°，6°，8°

2.2.2 氣動性能分析

圖4給出了2.2.1中所述工況下，潔凈翼型和結冰翼型的升阻比的比較結果.圖中CL，CD分別表示升力系數和阻力系數，其表達式為:

式中:L，D分別為翼型所受的升力和阻力，其中升力方向垂直于來流方向，阻力方向沿著來流方向.從圖中可以看到，潔凈翼型的最大升阻比在攻角(θ)6°附近，約為30.而結冰翼型的最大升阻比降到18，下降約40%，且發生在攻角4°，較潔凈翼型前移 2°.在 0°、4°、6°和 8°攻角下，升阻比都有顯著下降，這說明結冰后升阻比下降，將使翼型的氣動性能下降，導致葉片翼型失速提前發生且使風能的轉化效率降低.

圖4 潔凈翼型與結冰翼型的氣動性能Fig.4 Aerodynamic performance of clean airfoil nd ced

2.3 結冰厚度分析

結冰厚度對葉片翼型上邊界層分離的影響至關重要，為此，文中又選取風速6～20 m/s，環境溫度－2℃ ～ －5℃，LWC=1.0 g/m3、水滴直徑為 30 μm、攻角0°條件下，風速、溫度與結冰厚度的關系，定義結冰厚度為翼型表面法線方向的冰層厚度.圖5給出了結冰中最大結冰厚度(Hmax)與環境溫度(TA)、風速(vA)的關系，如圖5所示.

從圖5可以看到，在同一風速下，光冰結冰厚度隨著環境溫度的升高而降低，這是因為環境溫度高，水滴撞擊翼型表面后保持液態的可能性就大，在翼型表面隨氣流向下游運動，最終飛離翼型，這樣結冰量就減少.但在同一環境溫度下，隨著風速增加而結冰厚度增加，基本呈線性關系.其主要原因是氣動結冰中，水滴的慣性占主要作用，速度越大，與翼型碰撞的機率就越大，結冰量就增加.

圖5 環境溫度、風速與結冰厚度的關系Fig.5 Relation between ice height，air velocity and air temperature

3 結論

對水平軸風力機葉片翼型光冰結冰的情況進行了數值研究，得到不同風速、不同環境溫度、不同攻角下葉片翼型光冰結冰的形狀、冰增長的過程.結果表明:光冰結冰厚度、結冰量與風速的關系較大，隨著風速增大，其結冰量也加大，基本呈線性關系，而在同一風速下，不同環境溫度與光冰結冰量呈反比關系.對結冰翼型氣動性能的分析結果表明，邊界層在結冰翼型尾緣提前分離，并產生較大的分離渦，這將減少翼型的升力，同時增大阻力，其最大升阻比下降約40%.葉片結冰對氣動性能影響甚大，使升阻比顯著下降，失速提前，這將嚴重影響風力機將風能轉化為電能的轉換效率的同時，對葉片本身安全造成影響.

References)

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

［2］ 楊勝華，林貴平.機翼結冰過程的數值模擬［J］.航空動力學報，2011，26(2):323 －330.Yang Shenghua，Lin Guiping.Numerical simulation of ice accretion on airfoils［J］.Journal of Aerospace Power，2011，26(2):323－330.(in Chinese)

［3］ 陳維建，張大林.瘤狀冰結冰過程的數值模擬［J］.航空動力學報，2005，20(3):472 －476.Chen Weijian，Zhang Dalin.Numerical simulation of glaze ice accretion［J］.Journal of Aerospace Power，2005，20(3):472－476.(in Chinese)

［4］ 陳維建，張大林.飛機機翼結冰過程的數值模擬［J］.航空動力學報，2005，20(6):1010 －1017.Chen Weijian，Zhang Dalin.Numerical simulation of ice accretion on airfoils［J］.Journal of Aerospace Power，2005，20(6):1010 －1017.(in Chinese)

［5］ 張大林，陳維建.飛機機翼表面霜狀冰結冰過程的數值模擬［J］.航空動力學報，2004，19(1):137－141.Zhang Dalin，Chen Weijian.Numerical simulation of rime ice accretion process on airfoils［J］.Journal of Aerospace Power，2004，19(1):137 －141.(in Chinese)

［6］ 周志宏，李鳳蔚，李廣寧.基于兩相流歐拉方法的翼型結冰數值模擬［J］.西北工業大學學報，2010，28(1):138－142.Zhou Zhihong，Li Fengwei，Li Guangning.Applying eulerian droplet impingement model to numerically simulating ice accretion but with some improvements［J］.Journal of Northwestern Polytechnical University，2010，28(1):138－142.(in Chinese)

［7］ 李聲茂，李巖.垂直軸風力機葉片翼型結冰的數值計算［J］.動力工程學報，2011，31(3):214－219.Li Shengmao，Li Yan.Numerical simulation on icing of a blade aero foil for vertical axis wind turbines［J］.Journal of Chinese Society of Power Engineering，2011，31(3):214－219.(in Chinese)

［8］ 朱程香，王瓏，孫志國，等.風力機葉片翼型的結冰數值模擬研究［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.(in Chinese)

［9］ 任鵬飛.結冰風力機葉片的空氣動力學特性數值研究［D］.北京:中國科學院研究生院(工程熱物理研究所)，2014.

［10］ 易賢.飛機積冰的數值計算與積冰試驗相似準則研究［D］.四川綿陽:中國空氣動力研究與發展中心，2007.