水平軸風力機葉片翼型結冰的數值模擬

鄧曉湖,盧緒祥,李錄平,李海波

(長沙理工大學能源與動力工程學院,能源高效清潔利用湖南省普通高等學校重點實驗室,湖南長沙410114)

風力機組一般安裝于高寒、沿海地帶,風力發電機的漿葉經常會出現覆冰現象,導致槳葉的性能以及風力機功率的輸出達不到設計要求。風力機葉片的結冰是由過冷水滴撞擊葉片表面后凍結成冰覆蓋在葉片表面形成的。在結冰過程中若撞擊到葉片表面的水滴立即凍結,生成的冰型就比較規則,稱為霜狀冰;若表面大部分區域的水滴不能立即完全凍結,未凍結的水就會在氣流驅動下沿表面流動,并在流動過程中逐漸凍結,形成的冰稱為瘤狀冰[1]。葉片的結冰影響風力機的氣動性能和運行安全,給風力機帶來極大危害。例如,覆冰會改變葉片的外部形狀和氣動性能,增大葉片阻力、減少升力,影響全機操縱性、穩定性;最終導致風能的轉化效率降低,嚴重的可導致葉片的損毀,發生運行事故,風力機葉片結冰已經成為保證風力機安全運行需要解決的迫切問題。早期對風力機結冰進行的研究都采用實驗的方法,例如Bose等人在1992年對處于寒冷氣候中的雙葉片水平軸風力機結冰進行了測量,但是并未指出結冰形狀和外部結冰條件的關系[2-3];Jasinski等人在1998年對NREL-S809翼型進行了冰風洞實驗,得出了該翼型結冰前后的性能參數[4],但是實驗的方法非常麻煩。目前計算流體力學(CFD)的模擬已經應用到了風力機葉片結冰模擬中,采用模擬方法對風力機葉片翼型結冰的研究通常分為三類[5]:對特定條件下的翼型結冰過程進行預測、對已知結冰形狀的翼型進行性能預測以及防冰除冰技術研究。不過國內這方面工作尚處于摸索階段,而且大多借鑒于對飛機機翼和輸電電纜結冰的相關研究經驗。本文將進行的是第一類研究,通過數值計算方法預測相關廠家提供的風力機葉片翼型前沿霜狀冰的形狀和增長過程,分析結冰對翼型繞流及氣動特性的影響。

1 求解過程

采用數值方法模擬葉片結冰一般分三步[6]:首先,求解流體力學的基本方程組從而得到物體繞流的流場解;其次,把流場解的結果作為定解條件進一步求解水滴的運動方程,從而確定水滴與物體的碰撞點;最后,按照一定的增長模型來確定與物體相碰撞的水滴的結冰以及冰的增長。

1.1 流場計算

采用SIMPLE算法來求解二維定常不可壓粘流的時均N-S方程:

連續方程可表示為:

x方向的動量方程可表示為:

y方向也有類似的表達式。

各項的具體表達式和方程的詳細求解方法參見文獻[7]。

1.2 水滴軌跡運動方程的求解

用CFD的方法進行翼型結冰數值模擬時,必須對每個水滴進行跟蹤,以確定水滴是否與翼型發生碰撞。在計算水滴軌跡時候做如下假設[8]:①水滴的體積保持不變 ,但是水滴的形狀可以改變 ,引入與水滴體積相等的當量球的概念,當量球的直徑為d eg;②水滴的密度在整個過程中保持不變;③水滴的初始速度與自由流的速度相等,水滴體積很小以至于它們的繞流不會影響流場的性質。水滴軌跡運動方程可以寫成[6]:

式中:Md為水滴質量,ρd是水滴密度,ρa是空氣密度;V d是水滴體積;g是重力加速度;A d是水滴迎風面積;C d是阻力系數;u代表當地氣流速度;u d表示水滴速度。

整理得:

上式為二階常微分方程,可以把它寫成:

基于流場的速度分布所給定的初始條件,很明顯,軌跡運動方程的求解可以看成是一個一階常微分方程的初值問題:

可以采用四階龍格-庫塔法求解,最常用的公式是[8]:

當得到t時刻水滴的速度分ud(t),νd(t)后,則t+Δt時刻水滴的位置可以表示為:

開始計算時,首先給定水滴的初始位置,然后計算Δt時間后水滴的新位置。對于每一個水滴要分別跟蹤,如此推進計算,直到水滴與翼面相碰撞或者水滴移出界定區域。

1.3 霜冰厚度的計算模型

對于霜冰情形,由于所有的水滴在碰撞時刻就完全凝結并且不再融化,所以可以不考慮相變熱傳導等能量導致的變化,只考慮質量平衡就已經足夠[9-10]。

圖1給出了進入微元控制體中的質量流;微元控制體收集的液態水的質量m1;微元控制體內凍結的水的質量m2。微元控制體的質量平衡方程為:

圖1 有限容積內質量平衡示意圖

將翼型表面分割成很小的一塊塊區域,這些小區域的面積就是在翼型表面上所對應的弧長。前面由水滴運動方程的求解已經計算到了某個水滴與翼型相碰撞的位置,當水滴碰撞到某個編號為“i”的區域之內時,在計算程序處理中,令該區域所接收到的水滴數量上加1,即:

跟蹤計算完每個水滴后,最后撞擊到翼型上弧長為d si的一小區域之內的水滴總數即為ni[11-13]。根據質量守恒,冰的密度為300 kg/m3,水的密度為1000 kg/m3,則在已經編號的區域之內,可得:

式中:V d是單個水滴的體積;hi是弧長為 Δsi的區域上的結冰厚度。

單個水滴的體積為:

將式(19)帶入式(18),則可以得到:

2 結冰計算結果分析

當翼型形狀有一定程度的改變后,又對新的結冰翼型的繞流流場重新進行計算,以某廠家自主研發提供的翼型作為結冰表面,然后求解水滴軌跡運動方程和積冰厚度,如此反復迭代,直到需要的時間為止。

本文中所模擬的風力機葉片參數和結冰氣象條件為:自由來流馬赫數Ma為0.029,風速取固定的10 m/s,時間 T1和 T2分別為6 m in和10 m in,雷諾數Re為2.6×106,環境壓力Ph為101300 Pa。環境溫度T h=-10℃,冰密度取300 kg/m3;葉片長3.75 m,攻角AOA為0°,水滴有效平均直徑M vd為20μm,液態水含量LWC為1 g/m3,假設冰層與翼面間絕熱。

圖2 翼型一定時間后的積冰形狀

圖2的(a)和(b)分別給出了6min和10min時風力機槳葉翼型的積冰分布的計算結果。由于翼型前沿的繞流,水滴大部分凝結在翼型的下部,從圖中可以看出,在上下結冰極限附近由于撞擊到結冰表面的水立即凍結,產生的冰形在滯止區域具有最大厚度,同時冰形比較規則,和文獻[2-3]的結果大致相符合。但是從圖中還可以看出,由于溫度原因,撞擊到翼型表面的過冷水滴還會有少量未完全凍結,這些未凍結的水會沿著冰面產生回流,然后再逐漸凍結,以至在滯止區下游形成棱角。總之,在霜狀冰結冰預測方面,本文所預測結果與文獻結果總體趨勢一致,說明本文的預測方法是有效的,預測結果對風力機防冰除冰以及除冰裝置的布置具有指導作用,同時也可以在此基礎上對風力機葉片在惡劣條件下氣動性能進行研究。

3 結冰翼型的氣動模擬

3.1 控制方程與邊界條件

流場的連續方程:

動量方程:

能量方程:

上面三個公式中ρ為空氣密度,v為速度矢量,t為時間,p為壓強為應力張量,E為內能,k′為熱傳導系數。

在對控制方程做雷諾平均后得到雷諾平均的NS方程,并采用最佳湍流模型,用于結冰計算。

3.2 湍流模型

本文選擇最為適合的Spalart-A llmaras湍流模型。該模型最早被用于有壁面限制情況的流動計算中,特別在存在逆壓梯度的流動區域內,對邊界層的計算效果較好,因此經常被用于流動分離區附近的計算,通過模型與壁面函數配合可以適用于網格比較粗糙的算例,并且可以在網格精度不高時得到比較精確的解。SpalartA llmaras模型的求解變量是,表征出了近壁區域以外的湍流運動粘性系數。的輸運方程為:

其中Gv是湍流粘性產生項;Yv是由于壁面阻擋與粘性阻尼引起的湍流粘性的減少是用戶自定義源項和是常數;v是分子運動粘性系數。湍流粘性系數用如下公式計算:

3.3 網格生成

在得出翼型結冰邊界后,對新邊界進行網絡劃分進而對其氣動性能進行分析和對比。計算中采用Hilgenstoke法生成技術生成計算網格。根據數值計算的需要,采用分區網格的技術和思路,利用各自的優勢和特點,生成貼體、與邊界正交的、適合N-S方程計算的高質量網格,并對網格疏密及光滑性進行適當控制。在對結冰后翼型進行網格劃分時可能會出現尖角或波動扭曲,這會對網格生成及數值計算帶來困難,采用保證控制體結冰量不變的原則,對新生成的邊界進行光順處理,圖3為結冰后網格形狀。

圖3 結冰翼型網絡劃分

數值CFD模擬計算中,采用壓力基耦合求解器,對來流采用遠場邊界條件;并且在初始化時使用Full M ultigrid(FMG)初始化,以獲得更好的初始流場數據。計算基于RANS方法,采用二維穩態分離解法的隱式解法,空間離散格式采用二階迎風格式,壓力—速度耦合采用SIMPLE解法。所采用的邊界條件為:固壁表面采用無滑移條件、進口、出口和Interior面由特征相容條件所確定、計算邊界由插值確定。

圖4 結冰翼型表面速度分布

圖4為結冰翼型翼面附近流場的速度分布的模擬結果,可以明顯的看出,由于翼型前緣積冰的不規則,翼型上的速度分布非常不均勻。這種不均勻的速度分布對翼型后部氣流形成持續干擾,最終增加了整個翼型邊界層的不穩定性,促使葉片上的流動提前分離,從而導致翼型壓差阻力增大,阻力系數也隨之增大。此外,翼型前緣不規則的冰型增加了翼型的表面粗糙度,從而增加了翼型的摩擦阻力。

4 結冰翼型氣動模擬結果分析

4.1 結冰翼型的表面壓力

圖5的(a)和(b)分別是6 m in和10 min時光潔翼型和結冰翼型的壓力分布對比曲線,圖中只對結冰翼型進行了標明。由光潔翼型的壓力分布曲線可以看出,在翼型上表面沿著曲線由前緣點至后緣點方向靜壓值變小,但是變化趨勢較平緩;在翼型下表面沿著翼型曲線由前緣點至后緣點方向靜壓值增大,變化曲率也較大,即變化較為迅速。將兩種積冰翼型的壓力分布曲線與光滑翼型進行比較,可以看出積冰后在下翼面壓力系數峰值附近,曲線的下降忽然變陡并出現振蕩,上翼面壓力系數前沿有所減小,因此,翼型前沿氣動特性已經被破壞,隨著結冰厚度的加劇,破壞也越明顯。

圖5 翼型的壓力分布曲線

4.2 結冰翼型的升阻力系數

圖6和圖7分別給出了翼型在潔凈和結冰狀態下的升阻力系數對比曲線,可以看出結冰翼型的升力系數普遍低于干凈翼型,在模擬氣象條件下,槳葉最大升力系數由原來的1.1876下降到了0.8642,下降了27.23%。同時,結冰后翼型的失速發生在攻角為11°的時候,比干凈翼型的15°提前了4°,也就是說結冰翼型比干凈翼型提前失速。由阻力系數可以看出,翼型前緣帶冰會增加風力機葉片的阻力。阻力系數平均增加了38.65%,特別是失速攻角大幅度下降,對阻力特性的影響最大,提前進入失速區也是導致阻力系數增大的主要原因。

圖6 翼型升力系數對比曲線

圖7 翼型阻力系數對比曲線

5 結論

本文采用四階龍格-庫塔法求解水滴運動軌跡,并且采用了全N-S方程對風力機葉片翼型的積冰進行了預測,之后利用FLUENT軟件模擬了風力機葉片翼型結冰后周圍流場的變化,并且與結冰前葉片翼型的氣動性能進行了對比。

(1)采用求解水滴軌跡的方程和結冰厚度計算模型的方法可以有效預測翼型的積冰發展過程和形狀,表明這種方法正確有效。

(2)對于風力機葉片的翼型來說,由于前沿發生的繞流,結冰現象大部分是發生在翼型前沿的下部,這和飛機翼型的結冰結論相吻合。

(3)通過對比升力阻力性能,發現模擬氣象條件下的結冰翼型與干凈翼型相比,結冰翼型的最大升力系數大約減少了27%,阻力系數增加了約38%,失速攻角降低了4°。結冰后翼型提前進入失速區是造成氣動性能惡化的主要原因。

[1] MAKKONEN L,LAAKSO T,MARJANIEM IM,et al.Modelling and p revention of ice accretion on wind turbine on w ind turbines[J].W ind Engineering,2001, 25(1):3-21.

[2] BOSE N.Icing on a small horizontal-axis w ind turbine-part I:glaze ice profiles[J].Journalofw ind Engineering and Industrial Aerodynam ics,1992,45:75-85.

[3] BOSE N.Icing on a small horizontal-axis w ind turbine-part II:three dimensional ice and wet snow formations[J].Journalof Wind Engineering and Industrial Aerodynam ics,1992,45:87-96.

[4] JASINSK IW J,NOE S C,SELIG M S,et al.Wind turbine performance under icing conditions[J].Journal of So lar Energy Engineering,1998,120:60-65.

[5] ADDY H E,POTAPCZUKJR M G,SHELDON D W.Modern Airfoil Ice Accretions[R].AIAA-97-0174,1997.

[6] 易賢等.翼型結冰的數值模擬[J].空氣動力學學報, 2002,20(4):428-433.

[7] ZHU Guolin,K ronstam.The Calculation of ground effect on a car flow field using tw o dimensionalNavier-Stokes equations[J].A cta Aerody Nam ica Sinica, 1993,11(1):.

[8] 韓鳳華.飛機機翼表面水滴撞擊特性計算[J].北京航空航天大學學報,1995,21(3):18-21.

[9] FU P,FARZANEH M,BOUCHARD G.Tw o-dimensional modeling of the ice accretion p rocess on transm ission line wires and cables[J].Cold Regions Science and Techno logy,2006,46(3):132-146.

[10] 李巖,田川公太郎.葉片附著物對直線翼垂直軸風力機性能的影響[J].動力工程,2009,29(3):292-296.

[11] MAKKONEN L,LAAKSO T,MARJANIEM I M, et al.Modelling and prevention of ice accretion on wind turbines[J].Wind Engineering,M ulti-Science Pub lishing Co Ltd.,2001,25(1):3-21.

[12] HOCHART C,FORTIN G,PERRON J.W ind turbine performance under icing conditions[J].Wind Energy,2008,11:319-33.

[13] WANG Kevin.Convective heat transfer and experimental icing aerodynam ics of wind turbine blades[D]. Winnipeg,Canada:University of Manitoba,2008.