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

(長沙理工大學(xué)能源與動(dòng)力工程學(xué)院,能源高效清潔利用湖南省普通高等學(xué)校重點(diǎn)實(shí)驗(yàn)室,湖南長沙410114)

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

1 求解過程

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

1.1 流場(chǎng)計(jì)算

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

連續(xù)方程可表示為:

x方向的動(dòng)量方程可表示為:

y方向也有類似的表達(dá)式。

各項(xiàng)的具體表達(dá)式和方程的詳細(xì)求解方法參見文獻(xiàn)[7]。

1.2 水滴軌跡運(yùn)動(dòng)方程的求解

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

式中:Md為水滴質(zhì)量,ρd是水滴密度,ρa(bǔ)是空氣密度;V d是水滴體積;g是重力加速度;A d是水滴迎風(fēng)面積;C d是阻力系數(shù);u代表當(dāng)?shù)貧饬魉俣?u d表示水滴速度。

整理得:

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

基于流場(chǎng)的速度分布所給定的初始條件,很明顯,軌跡運(yùn)動(dòng)方程的求解可以看成是一個(gè)一階常微分方程的初值問題:

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

當(dāng)?shù)玫絫時(shí)刻水滴的速度分ud(t),νd(t)后,則t+Δt時(shí)刻水滴的位置可以表示為:

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

1.3 霜冰厚度的計(jì)算模型

對(duì)于霜冰情形,由于所有的水滴在碰撞時(shí)刻就完全凝結(jié)并且不再融化,所以可以不考慮相變熱傳導(dǎo)等能量導(dǎo)致的變化,只考慮質(zhì)量平衡就已經(jīng)足夠[9-10]。

圖1給出了進(jìn)入微元控制體中的質(zhì)量流;微元控制體收集的液態(tài)水的質(zhì)量m1;微元控制體內(nèi)凍結(jié)的水的質(zhì)量m2。微元控制體的質(zhì)量平衡方程為:

圖1 有限容積內(nèi)質(zhì)量平衡示意圖

將翼型表面分割成很小的一塊塊區(qū)域,這些小區(qū)域的面積就是在翼型表面上所對(duì)應(yīng)的弧長。前面由水滴運(yùn)動(dòng)方程的求解已經(jīng)計(jì)算到了某個(gè)水滴與翼型相碰撞的位置,當(dāng)水滴碰撞到某個(gè)編號(hào)為“i”的區(qū)域之內(nèi)時(shí),在計(jì)算程序處理中,令該區(qū)域所接收到的水滴數(shù)量上加1,即:

跟蹤計(jì)算完每個(gè)水滴后,最后撞擊到翼型上弧長為d si的一小區(qū)域之內(nèi)的水滴總數(shù)即為ni[11-13]。根據(jù)質(zhì)量守恒,冰的密度為300 kg/m3,水的密度為1000 kg/m3,則在已經(jīng)編號(hào)的區(qū)域之內(nèi),可得:

式中:V d是單個(gè)水滴的體積;hi是弧長為 Δsi的區(qū)域上的結(jié)冰厚度。

單個(gè)水滴的體積為:

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

2 結(jié)冰計(jì)算結(jié)果分析

當(dāng)翼型形狀有一定程度的改變后,又對(duì)新的結(jié)冰翼型的繞流流場(chǎng)重新進(jìn)行計(jì)算,以某廠家自主研發(fā)提供的翼型作為結(jié)冰表面,然后求解水滴軌跡運(yùn)動(dòng)方程和積冰厚度,如此反復(fù)迭代,直到需要的時(shí)間為止。

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

圖2 翼型一定時(shí)間后的積冰形狀

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

3 結(jié)冰翼型的氣動(dòng)模擬

3.1 控制方程與邊界條件

流場(chǎng)的連續(xù)方程:

動(dòng)量方程:

能量方程:

上面三個(gè)公式中ρ為空氣密度,v為速度矢量,t為時(shí)間,p為壓強(qiáng)為應(yīng)力張量,E為內(nèi)能,k′為熱傳導(dǎo)系數(shù)。

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

3.2 湍流模型

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

其中Gv是湍流粘性產(chǎn)生項(xiàng);Yv是由于壁面阻擋與粘性阻尼引起的湍流粘性的減少是用戶自定義源項(xiàng)和是常數(shù);v是分子運(yùn)動(dòng)粘性系數(shù)。湍流粘性系數(shù)用如下公式計(jì)算:

3.3 網(wǎng)格生成

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

圖3 結(jié)冰翼型網(wǎng)絡(luò)劃分

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

圖4 結(jié)冰翼型表面速度分布

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

4 結(jié)冰翼型氣動(dòng)模擬結(jié)果分析

4.1 結(jié)冰翼型的表面壓力

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

圖5 翼型的壓力分布曲線

4.2 結(jié)冰翼型的升阻力系數(shù)

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

圖6 翼型升力系數(shù)對(duì)比曲線

圖7 翼型阻力系數(shù)對(duì)比曲線

5 結(jié)論

本文采用四階龍格-庫塔法求解水滴運(yùn)動(dòng)軌跡,并且采用了全N-S方程對(duì)風(fēng)力機(jī)葉片翼型的積冰進(jìn)行了預(yù)測(cè),之后利用FLUENT軟件模擬了風(fēng)力機(jī)葉片翼型結(jié)冰后周圍流場(chǎng)的變化,并且與結(jié)冰前葉片翼型的氣動(dòng)性能進(jìn)行了對(duì)比。

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

(2)對(duì)于風(fēng)力機(jī)葉片的翼型來說,由于前沿發(fā)生的繞流,結(jié)冰現(xiàn)象大部分是發(fā)生在翼型前沿的下部,這和飛機(jī)翼型的結(jié)冰結(jié)論相吻合。

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

[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] 易賢等.翼型結(jié)冰的數(shù)值模擬[J].空氣動(dòng)力學(xué)學(xué)報(bào), 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ī)機(jī)翼表面水滴撞擊特性計(jì)算[J].北京航空航天大學(xué)學(xué)報(bào),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] 李巖,田川公太郎.葉片附著物對(duì)直線翼垂直軸風(fēng)力機(jī)性能的影響[J].動(dòng)力工程,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.