液態(tài)水含量對寒區(qū)海上風力機葉片霜冰結冰的影響

曹慧清 白旭 尹群

(江蘇科技大學,船舶與海洋工程學院,江蘇 鎮(zhèn)江 212003)

0 引言

近年,風力發(fā)電行業(yè)發(fā)展迅猛,陸上風力發(fā)電已趨近飽和,許多企業(yè)紛紛把目光轉向海上風力機的建設。而處于寒冷地區(qū)的海上風力機在冬季運行時將存在結冰的風險,風力機結冰(如圖1)影響其氣動性能和運行安全,給風力機帶來了極大危害。例如,風力機表面結冰會降低風能利用系數,減小出力;影響風力機的結構強度與疲勞壽命,產生安全隱患。因此,研究海上風力機結冰預報與防除冰技術至關重要。

圖1 風力機葉片結冰Fig.1.Wind turbine blade icing

為了精確地對寒區(qū)海上風力機的結冰進行預報,首先需要研究影響風力機結冰的因素,包括環(huán)境溫度、液態(tài)水含量(Liquid Water Content,LWC)、平均水滴直徑(Medium Volume Droplet Diameter,MVD)、風速等。其中液態(tài)水含量(單位體積的空氣中所有液態(tài)水的質量)是影響結冰形狀的重要參數,結冰形狀改變將影響葉片的氣動性能,最終影響風能的轉化效率。空氣中液態(tài)水含量首先取決于地域特性,潮濕氣候條件下,液態(tài)水含量較高;其次取決于環(huán)境溫度,溫度低于–10℃,空氣中過冷水滴凝結為冰晶并降落到地面,使得空氣中液態(tài)水含量降低[1]。同時,已有研究表明,空氣中的氣態(tài)水含量幾乎不會對葉片結冰造成影響[2]。因此,本文以空氣中液態(tài)水含量為主要研究對象。

目前,已有國內外學者對風力機葉片結冰進行了相關試驗研究。孫策[3]研究了霜冰結冰條件下,結冰對風力機葉片的影響,初步提出了一套適用于小尺度、小縮比和霜冰條件下的葉片結冰形狀相似準則,利用自行搭建的結冰風洞試驗系統(tǒng)完成了對該準則的驗證性試驗。Shin 和Bond[4]研究了不同的液態(tài)水含量、風速等參數對NACA0012 翼型結冰的影響,并做了大量結冰試驗。Bose[5]對直徑為1.05 m 的水平軸風力機做了降水結冰的試驗,得到了葉片不同位置處明冰結冰的形狀,發(fā)現大多數的明冰都發(fā)生在葉片前緣和翼尖處。Hu 等[6]研究了葉片受不同參數影響的結冰分布。結果表明,冰的質量和厚度從葉根到葉尖呈線性增長。較高的風速、較小的攻角、較高的液態(tài)水含量和較低的溫度均對葉片結冰有較大影響。

相比于風力機風洞試驗研究,數值模擬研究更加廣泛。Homola 等[7]研究了NREL 5 MW 風力機模型因結冰造成的性能損失,并在霜冰條件下對沿葉片的五個截面結冰進行了模擬,并用計算流體力學(Computational Fluid Dynamics,CFD)方法模擬分析了潔凈和結冰翼型截面的氣動特性。顧聲龍等[8]利用Fluent UDF(User Defined Function)對水平軸風力機葉片翼型的明冰結冰進行了數值研究,得出了明冰對翼型氣動性能的影響。鄧曉湖等[9]通過CFD 方法對水平軸風力機葉片翼型的霜冰結冰進行了數值模擬,結果表明,結冰后翼型提前進入失速區(qū)是造成葉片氣動性能惡化的主要原因。結冰模型用于模擬在結冰條件下某結冰對象在流體與結構之間的相互作用關系。目前存在2 種常用冰積模型:一種是結冰經驗模型,最常見的是Makkonen 模型[10],該模型是基于準三維、較靈活的框架模型,是其他大多數結冰經驗模型的基礎;另一種是結冰理論模型,分為二維和三維葉片結冰模型,著名的結冰仿真商業(yè)軟件LEWICE、FENSAP-ICE 和TURBICE 都是此類模型,其中TURBICE 和LEWICE 為二維葉片模型,FENSAP-ICE 為三維葉片模型?，F階段,二維模型求解空氣流場使用的是面元法,主要用于求解勢流模型;而三維模型求解空氣流場使用的是空間網格法,用于求解N-S 方程或歐拉方程。

海上風力機有著區(qū)別于陸上的獨特氣象條件,陸上的結冰預報方法是否適用仍需驗證。由于葉片是風力機接受風能最重要的部件,本文針對寒區(qū)海上風力機葉片,采用Fluent 和FENSAP-ICE相結合的方法,選取液態(tài)水含量作為敏感參數,分析液態(tài)水含量的變化對結冰厚度、結冰量和結冰外形的影響。

1 結冰氣象條件與結冰理論

1.1 結冰氣象條件

空氣中的水滴只會在低于–40℃時自發(fā)結冰[11],而溫度在0～–40℃之間時,水滴處在亞穩(wěn)定狀態(tài)即過冷狀態(tài),經過很長一段時間后產生結晶核才會結冰。但是,若水滴與葉片發(fā)生碰撞,葉片便會開始結冰。

風力機葉片結冰按其結冰性質可分為霜冰、明冰、濕雪和混合冰。圖2 為霜冰、明冰與風速和溫度的關系。本文海上風力機的氣象條件以Drage 和Hauge[12]的研究為參考,他們通過實驗確定了挪威西海岸Brosviks?ta 山大氣結冰的氣象條件范圍(見表1)。結合表1 和圖2 可知,本文研究的結冰類型為霜冰。同時,世界氣象組織指出,海平面以上16 m 是海洋飛沫的上限[13],本文研究的海上風機葉片葉尖離海平面最近處約30 m,所以本文暫不考慮海洋飛沫對葉片結冰的影響。

圖2 冰型與氣溫和風速的關系[14]Fig.2.Relationship between ice type and temperature and wind speed[14]

表1 氣象條件Table 1.Meteorological condition

1.2 結冰機理及模擬

1.2.1 結冰機理

結冰是一種常見的自然現象,從熱力學角度來看,水的溫度降到0℃以下成為過冷水,過冷水并不穩(wěn)定,而是處于亞穩(wěn)態(tài),該狀態(tài)的解除需要大于臨界尺寸的冰核的形成。當過冷水中出現尺寸大于臨界尺寸的冰核時,結冰過程開始,冰核在過冷水中長大,最終成為宏觀意義上的冰[15]。風力機在運行過程中,當葉片表面周圍溫度低于0℃并與過冷水滴發(fā)生碰撞時,便會出現結冰現象。在結冰過程中,伴隨復雜的傳熱、傳質現象,撞擊到葉片表面的過冷水滴與空氣會產生對流換熱、自身伴隨蒸發(fā)和升華等,同時釋放潛熱,根據外界條件差異決定是否產生結冰表面的液態(tài)水回流等現象[16]。

霜冰主要出現在環(huán)境溫度較低、液態(tài)水含量較低和液滴大小較小情況下,這些條件下過冷水滴撞擊到葉片表面立即凍結。結冰過程中釋放的熱量不足以使葉片表面溫度上升到冰點以上,因此冰沒有發(fā)生融化,也無向后的溢流。所以,結冰發(fā)生在葉片前緣狹窄區(qū)域,外形較規(guī)則且不透明。

1.2.2 結冰數值模擬

采用數值模擬方法進行結冰計算,主要包括3 個步驟:葉片流場計算、水滴軌跡計算和冰層生長計算。

(1)葉片流場計算。計算葉片周圍的流場分布,以得到葉片周圍空氣流速等參數。

控制方程為低速黏流N-S 方程,通用形式為:

其連續(xù)方程為:

式中,ρ為空氣密度,t為時間,φ為運輸變量,?為空氣速度,φΓ 為擴散系數,qφ為源項。

風力機工作時葉片處于旋轉狀態(tài),此時空氣處于復雜繞流運動狀態(tài),因此本文使用葉素動量(Blade Element Momentum,BEM)理論來解決葉片旋轉和多葉片對空氣流場的影響[17]。

BEM 理論是分析風力機葉片氣動性能主要方法。空氣繞流葉片的運動是具有渦旋的復雜流動,BEM 理論是將風輪葉片沿展向分成多個微段,這些微段為葉素,通過引入軸向、切向誘導因子反映氣流沿軸向、切向速度的變化,使基于動量守恒原理和氣動理論的分析方法成為可能。

數值模擬包括兩個主要參數,相對風速(Vrel)和迎角(α0)。這些參數取決于來流風速(V∞)、軸向和切向誘導因子(a和γ)、截面半徑(r)、局部速比(λr)和該截面處的葉片扭轉角(?)。每個截面的相對風速和迎角通過下式計算:

其中:

其中,R是風輪半徑,CL為升力系數,B為葉片個數,c為葉片翼型弦長,λ為葉尖速比,λr為局部速比,β0是流動角,β0=90?α0?φ。

(2)水滴軌跡計算。水滴軌跡的計算使用的是FENSAP-ICE 軟件中DROP3D 模塊,在葉片流場計算結果的基礎上,求解過冷水滴的運動方程,對每個水滴進行跟蹤,以確定水滴是否與葉片發(fā)生碰撞。用歐拉法建立氣-液兩相控制方程,然后用有限體積法求解控制方程,從而得到水滴運動軌跡和葉片表面的撞擊特性。

水滴軌跡運動方程如下:

式中,ρa和ρd分別為氣流密度和水滴密度,Vrel為相對風速,和分別為氣流速度和水滴速度,d為水滴直徑,aμ為運動黏度,L∞為特征長度,g∞為來流重力加速度，→為重力加速度。

(3)冰層生長計算。冰層生長計算運用FENSAP-ICE 軟件中ICE3D 模塊,此模塊使用的是Bourgault 結冰模型,利用質量和能量守恒方程求解葉片表面的結冰熱力學模型,得到葉片的表面溫度和結冰量的分布。

質量守恒方程:

能量守恒方程:

式中,f表示水膜,為瞬時蒸發(fā)質量,為瞬時結冰質量,fρ、fc、sc、ε、σ、Levap和Lfusion分別代表流體的物理性能,Qanti-icing為防冰熱通量,為表面法線向量,LWC為局部液態(tài)水含量,β局部水滴收集系數。各項的具體表達式和方程的詳細求解方法參見文獻[18]。

當結冰類型為霜冰時,液滴撞擊葉片表面會直接結冰,則不需考慮傳熱及相變,只需考慮質量守恒。其質量守恒方程可以寫為:

式中,ρice為冰的密度,tΔ 為時間步長,為結冰增長速度,Δhice結冰表面位移,LWC∞為來流液態(tài)水含量。

結合公式(16)～(18)可以看出,在其他參數條件給定時,LWC的改變會直接影響到葉片的結冰質量和結冰生長速度。

2 數值預報方法及驗證

2.1 計算外形及網格

本文計算外形選用某5 MW 級海上風力機[19],其單葉片長63 m。風力機葉片(見圖3)上距離翼根85%～100%之間的區(qū)域對于風力發(fā)電而言最重要,本文采用位于風力機葉片上距離翼根90%的位置的翼型作為研究對象[20]。此處的翼型為NACA64-A17,翼型弦長為2.086 m,厚度為0.05 m。

光滑翼型與帶冰翼型均采用ICEM 劃分的C型結構化網格(見圖4),翼型固壁面采用無滑移和絕熱條件。

圖4 計算域網格劃分及局部放大圖。a) 潔凈翼型計算域網格劃分情況;b) 潔凈翼型局部網格劃分情況;c) 結冰翼型局部網格劃分情況Fig.4.Grids in computational zone and part grid of the airfoil.a) grids in computational zone of the clean airfoil;b) partial grid of a clean airfoil;c) partial grid of an iced airfoil

2.2 計算流程

本文結合Fluent 和FENSAP-ICE 軟件進行風力機葉片的結冰數值模擬研究。利用ICEM 軟件劃分潔凈翼型截面網格;將劃分好的網格導入Fluent 軟件中分析葉片周圍流場;將流場分析結果導入FENSAP-ICE 中,在FENSAP-ICE 軟件中DROP3D模塊計算水滴軌跡,ICE3D 模塊計算結冰生長并導出結冰后翼型;由于霜冰主要發(fā)生在翼型的前緣,翼型前緣需加密處理,重新在ICEM 中劃分結冰后的翼型截面網格并導入Fluent,重復上述過程直到給定的截止時間。圖5 為本文的計算流程圖。

圖5 計算流程Fig.5.Calculation process

2.3 數值驗證

為了驗證本文數值模擬結果的準確性,本文將數值模擬結果與試驗結果進行對比。其中,NACA0012 翼型結冰的數值模擬是最經典的算例,有著豐富的冰風洞試驗數據。Shin 和Bond[21]在美國國家航空咨詢委員會(NACA)的Lewis 研究中心對NACA0012 翼型結冰進行了多次冰風洞試驗,積累了大量的試驗數據,本文選取其中一組試驗結果進行對比。

文獻中的計算條件為:翼型弦長c=0.5334 m,攻角AOA=4°,環(huán)境溫度T=–28.3℃,來流速度U∞=67.05 m·s–1,平均水滴直徑MVD=20 μm,液態(tài)水含量LWC=1 g·m–3,結冰時間t=360 s。

由圖6 可知,本文預測冰形與試驗結果吻合較好,結冰外形為霜冰,霜冰主要集中在翼型前緣,結冰極限與結冰厚度基本一致,結冰冰形大部分重合,由此驗證了本文計算方法的正確性。

圖6 本文計算結果與試驗結果對比Fig.6.Comparison of the ice shape between this paper and the experiment

3 液態(tài)水含量對葉片結冰的影響

3.1 計算條件

本文采用了如表2 所示的環(huán)境條件和氣象條件,以此來計算翼型表面的結冰量和結冰外形,并分析不同液態(tài)水含量和在1 h、3 h 和10 h 三種不同結冰時間條件下對葉片結冰的影響。

3.2 時間步長選擇

在結冰數值模擬的計算中,時間步長的選擇會對數值模擬的計算結果造成干擾。隨著翼型結冰后表面局部收集系數β的改變,結冰數值模擬的計算結果也受到影響。對于CFD 分析軟件來說,時間步長選取的越小,計算結果越精確,但對計算機的性能要求越高,計算時間也隨之增加,因此有必要選取合適的時間步長。

本文選取表2 中一組結冰條件進行研究,其中液態(tài)水含量為0.15 g·m–3,結冰時間為3 h,其他結冰條件不變。時間步長分別選取5 min、15 min和30 min,表3 為三種不同時間步長時的數值計算結果。從表3 中可知,相較于時間步長為5 min 時的結冰量,15 min 和30 min 時間步長下結冰量的變化幅值百分比分別為 1.30%和4.84%。三種時間步長模擬結果相近,但在計算中發(fā)現,隨著結冰時間的增加,15 min 和30 min相較于5 min 時間步長下結冰量的變化增幅呈不斷增大的趨勢。因此,考慮到三種時間步長下的計算效率和計算精度,本文選取15 min 作為時間步長。

表2 不同液態(tài)水含量對結冰影響的計算參數Table 2.Calculation parameters of influence of different liquid water content on icing

表3 三種時間步長數值計算結果Table 3.Three time-step numerical calculation results

3.3 不同液態(tài)水含量下的計算結果分析

圖7 至圖9 分別為液態(tài)水含量為0.05 g·m–3、0.10 g·m–3、0.15 g·m–3、0.20 g·m–3和0.25 g·m–3時,在1 h、3 h 和10 h 三種不同結冰時間條件下,葉片所取截面的結冰量、結冰最大厚度和結冰形狀圖。

圖7 不同液態(tài)水含量下葉片表面的結冰量Fig.7.Amount of icing on blade surface under different liquid water content

從圖7 和圖8 可以看出,在相同的結冰條件下,隨著空氣中液態(tài)水含量的增加,翼型表面的結冰量和結冰厚度也逐漸增加。這是由于液態(tài)水含量代表單位體積空氣中所含有的液態(tài)水的質量,在其他條件相同的情況下,隨著液態(tài)水含量的增加,單位體積空氣中存在著更多的水滴,將有更多的水滴撞擊到葉片表面,從而使葉片表面微元中收集到的液態(tài)水更多,并在之后的傳質和傳熱過程中凍結,形成更多的冰。同時,從圖中可以看出,隨著結冰時長的增加,不同濃度液態(tài)水含量下結冰量的增幅呈不斷變大的趨勢。這是由于隨著結冰時間的增加,不同濃度液態(tài)水含量下結冰外形差異更加明顯,結冰表面的局部收集系數β也隨之變化,因此結冰量不再隨液態(tài)水含量變化呈線性增加。

圖8 不同液態(tài)水含量下葉片表面的結冰最大厚度Fig.8.Maximum thickness of icing on blade surface under different liquid water content

從圖9 中可以看出,不同濃度液態(tài)水含量下生成的冰形均為霜冰,未發(fā)現明顯溢流現象。液態(tài)水含量的改變并沒有改變葉片表面的結冰類型,也未影響結冰的生長趨勢。不同濃度液態(tài)水含量下,冰與葉片的吸附面沒有產生明顯變化,這是因為盡管空氣中液態(tài)水含量增加了,但風速、平均水滴直徑等參數并未發(fā)生改變,水滴的運動軌跡也沒有改變。

圖9 不同液態(tài)水含量下葉片表面的結冰外形。a) 結冰1 h;b) 結冰3 h;c) 結冰10 hFig.9.Shapes of icing on blade surface under different liquid water content.a) 1 h-icing;b) 3 h-icing;c) 10 h-icing

4 結論

本文以寒區(qū)海上風力機葉片結冰為背景,采用Fluent 和FENSAP-ICE 相結合的方法對葉片結冰開展數值模擬,選取液態(tài)水含量作為敏感參數,分析了液態(tài)水含量的變化對葉片結冰厚度、結冰量和結冰外形的影響,得出如下結論。

(1)液態(tài)水含量在0.05～0.25 g·m–3范圍內,隨著液態(tài)水含量的增加,葉片表面的結冰量和結冰厚度隨之增加。

(2)液態(tài)水含量在0.05～0.25 g·m–3范圍內,液態(tài)水含量并不影響葉片表面的結冰類型和結冰生長趨勢,同時冰與葉片的吸附面也不隨液態(tài)水含量變化而改變。

(3)采用不同的液態(tài)水含量對葉片結冰進行模擬,雖然模擬出了寒區(qū)海上液態(tài)水含量的變化范圍,但是液態(tài)水含量會隨海拔高度變化而變化,本文研究的海上風力機塔架高達90 m,葉尖高度可達海平面以上153 m,葉片旋轉時,短時間內會掠過不同液態(tài)水含量的空氣,對液態(tài)水含量時刻變化的空氣的模擬還有待改進。

(4)針對不同結冰時長、不同液態(tài)水含量的結冰進行了分析。若風速、平均水滴直徑等參數發(fā)生改變,對不同液態(tài)水含量下的結冰進行對比分析,是否能得到類似的結論還有待研究。