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

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

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

葉片結(jié)冰影響風(fēng)力機的氣動性能和運行安全，給風(fēng)力機帶來了極大危害.例如，覆冰會改變?nèi)~片的外部形狀和氣動性能，增大葉片阻力、減少升力，影響全機操縱性、穩(wěn)定性，最終導(dǎo)致風(fēng)能的轉(zhuǎn)化效率降低，嚴重的甚至?xí)?dǎo)致葉片的損毀，發(fā)生運行事故.風(fēng)力機葉片結(jié)冰已經(jīng)成為保證風(fēng)力機安全運行需要解決的迫切問題.

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

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

1 求解方法

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

對于連續(xù)相和離散相的求解，文中利用FLUENT軟件實現(xiàn)，利用FLUENT用戶自定義函數(shù)UDF自主開發(fā)出進行碰撞效率以及結(jié)冰增長的程序.

1.1 連續(xù)相的計算

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

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

1.2 離散相的計算

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

式中:md為水滴質(zhì)量;ρ為空氣密度;Cd為阻力系數(shù);u為當(dāng)?shù)貧饬魉俣?ud為水滴速度;Ad為水滴的迎風(fēng)面積.對上式積分即可獲得水滴運動軌跡.

1.3 光冰結(jié)冰量的計算

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

文中采用控制體內(nèi)質(zhì)量守恒和能量守恒原理［10］.

1)質(zhì)量守恒

2)能量守恒

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

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

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

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

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

2 結(jié)果分析

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

2.1 結(jié)冰過程及形狀

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

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

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

2.2 結(jié)冰翼型的氣動性能分析

2.2.1 流場分析

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

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

2.2.2 氣動性能分析

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

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

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

2.3 結(jié)冰厚度分析

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

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

圖5 環(huán)境溫度、風(fēng)速與結(jié)冰厚度的關(guān)系Fig.5 Relation between ice height，air velocity and air temperature

3 結(jié)論

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

References)

［1］ 鄧曉湖，盧緒祥，李錄平，等.水平軸風(fēng)力機葉片翼型結(jié)冰的數(shù)值模擬［J］.能源技術(shù)，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］ 楊勝華，林貴平.機翼結(jié)冰過程的數(shù)值模擬［J］.航空動力學(xué)報，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］ 陳維建，張大林.瘤狀冰結(jié)冰過程的數(shù)值模擬［J］.航空動力學(xué)報，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］ 陳維建，張大林.飛機機翼結(jié)冰過程的數(shù)值模擬［J］.航空動力學(xué)報，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］ 張大林，陳維建.飛機機翼表面霜狀冰結(jié)冰過程的數(shù)值模擬［J］.航空動力學(xué)報，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］ 周志宏，李鳳蔚，李廣寧.基于兩相流歐拉方法的翼型結(jié)冰數(shù)值模擬［J］.西北工業(yè)大學(xué)學(xué)報，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］ 李聲茂，李巖.垂直軸風(fēng)力機葉片翼型結(jié)冰的數(shù)值計算［J］.動力工程學(xué)報，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］ 朱程香，王瓏，孫志國，等.風(fēng)力機葉片翼型的結(jié)冰數(shù)值模擬研究［J］.空氣動力學(xué)學(xué)報，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］ 任鵬飛.結(jié)冰風(fēng)力機葉片的空氣動力學(xué)特性數(shù)值研究［D］.北京:中國科學(xué)院研究生院(工程熱物理研究所)，2014.

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