張漢文，羅日成，劉志勇，姜貽哲

（長沙理工大學 電氣與信息工程學院，湖南 長沙410114）

0 引言

近年來，我國風電場建設速度較快，在中南、西南等地區(qū)新建了多個風電場基地，但這些地區(qū)春冬季節(jié)的濕冷氣候也使得風力發(fā)電機面臨著嚴峻的覆冰問題[1]，[2]。風力發(fā)電機覆冰的主要危害體現(xiàn)在風機葉片氣動性能下降導致風機出力降低，風機荷載增大造成主軸磨損加劇以及覆冰可能會在風機旋轉過程中甩出發(fā)生安全事故[3]，[4]。

風機葉片覆冰是指葉片在旋轉過程中持續(xù)捕獲空氣中的過冷卻水滴而在其表面發(fā)生凍結的物理現(xiàn)象。這個過程中影響風機葉片覆冰的因素包括環(huán)境因素和自身因素。環(huán)境因素包括風速、氣溫、濕度等。自身因素包括翼型、葉片尺寸、葉片數(shù)量、葉片表面涂層、葉尖速比λ等。

由于試驗條件、經(jīng)費等因素限制，目前國內外學者主要采用數(shù)值模擬法來研究風機葉片的覆冰問題。數(shù)值模擬法是指通過計算機軟件建立風機的幾何模型，利用有限元法劃分網(wǎng)格，再分別對空氣流場、水滴碰撞以及質量和熱量的傳遞過程進行計算，從而獲得葉片表面的覆冰增長過程。文獻[5]通過LEWICE系統(tǒng)將三維旋轉葉片覆冰問題簡化為二維葉片覆冰問題，并通過試驗驗證了該模型模擬霧凇覆冰的準確性。文獻[6]基于二維覆冰模型完成了對風機的雨淞、霧凇覆冰的模擬，并計算分析了風機在覆冰工況下的轉矩損失。文獻[7]基于氣液兩相流和Messinger模型建立了水平軸風機的三維雨淞覆冰預測模型，并對小型水平軸風力機進行了人工覆冰試驗驗證。文獻[8]利用拉格朗日法研究了三維旋轉坐標系下的水平軸風力機葉片的水滴撞擊特性。

目前，我國在垂直軸風機葉片覆冰問題的研究成果較少。本文采用數(shù)值模擬法研究不同葉片翼型、葉片數(shù)量以及λ對垂直軸風機葉片覆冰的影響。該研究有助于幫助人們了解垂直軸風機葉片的覆冰特性，實現(xiàn)對特定氣候下的風機葉片表面覆冰形貌預測，為風電場的防冰、除冰系統(tǒng)研發(fā)提供幫助。

1 數(shù)學模型

1.1 空氣流場計算

旋轉風機的空氣流場屬于低速流場，可近似認為空氣是不可壓縮流體，空氣的粘性會導致葉片表面氣流邊界層分離。因此，本文空氣流場的控制方程選用低速粘流的雷諾平均N-S方程，湍流模型采用對邊界層湍流和自由剪切湍流均有較好模擬效果的SST k-ω湍流模型[9]。該模型在近壁面利用k-ω模型的魯棒性以捕捉到粘性邊界層的流動，而在主流區(qū)域采用k-ε模型，k-ε模型能夠較準確地模擬旋轉風機的空氣流場，已經(jīng)廣泛用于旋轉機械問題的求解[10]。本文湍流模型的輸運變量為湍流動能k和比耗散率ω，其輸送方程為

式中：u為氣流速度場；ρ為空氣密度；μ為空氣粘性系數(shù)；t為時間；Pk為k的有效生成率；Pω為ω的有效生成率；β，β*，σk，σω，σω2為經(jīng)驗系數(shù)；F1為混合函數(shù)。

1.2 過冷水滴的運動軌跡

為簡化模型，假設過冷卻水滴在空氣流場中是均勻分布，且是直徑、密度相等的球體，運動過程中僅受空氣阻力、重力和浮力的作用，且重力和浮力相等。本文基于離散顆粒模型，利用拉格朗日法對單個過冷卻水滴建立受力平衡運動方程；借助空氣流場條件和過冷卻水滴的初始釋放條件，將運動方程轉化為常微分方程；通過龍格庫塔法求解水滴的運動微分方程，得到水滴在空氣流場內的運動軌跡，求得風機葉片表面所捕獲的過冷卻水滴數(shù)量及分布情況。單個過冷卻水滴的受力分析如下：

1.3 覆冰質量計算

實際上，風機葉片的覆冰過程是一個非定常過程，即隨著葉片表面覆冰厚度不斷增長，空氣流場的解也會發(fā)生變化，從而又影響到過冷卻水滴的運動軌跡，這須要不斷對覆冰邊界進行重構。為了簡化模型，假定風機葉片的整個覆冰過程為定常過程。同時，假定過冷卻水滴與風機葉片碰撞瞬間發(fā)生結冰，且冰沿著葉片表面法向方向增長。參照Makkonen覆冰模型，在旋轉風機持續(xù)捕獲過冷水滴的過程中，其葉片表面覆冰的增長速率為

式中：α為過冷水滴碰撞系數(shù)；ω為空氣中液態(tài)水含量；v為過冷水滴的速度；A為垂直于風速矢量的截面積。

2 風機覆冰模擬

2.1 風機覆冰物理模型

本文主要利用COMSOL軟件中的湍流模塊和粒子追蹤模塊來對垂直軸風機葉片的覆冰問題進行模擬研究，垂直軸風機覆冰的物理模型如圖1所示。過冷卻水滴由入口處釋放并受到空氣流場的作用運動，旋轉域中的葉片不斷捕獲來流過冷卻水滴，被捕獲的水滴黏附在葉片壁面，未被捕獲的水滴由兩側的對稱邊界和出口處流出。模型中的主要參數(shù)：來流風速為10m/s，出口為壓力邊界，壓力為0MPa，環(huán)境溫度為-5℃，水滴直徑為20μm，入口邊界處每間隔0.01 s均勻釋放5000個水滴顆粒，旋轉域轉速由試驗組參數(shù)決定。整個空氣流場空間大小為10m×5m，風機直徑為2.2 m。

圖1 風機覆冰物理模型Fig.1 Physicalmodel of the VAWT icing

2.2 試驗組參數(shù)

以升力型垂直軸風機常用的NACA0018和NACA6510翼型葉片為仿真研究對象，葉片弦長均為500mm（圖2）。根據(jù)不同的葉片翼型、葉片數(shù)量以及λ，將仿真共分為12組（表1）。

圖2 葉片翼型Fig.2 Blade airfoils

表1 試驗組參數(shù)Table 1 Parameters of experimental groups

續(xù)表1

3 覆冰模擬結果及分析

3.1 冰型模擬結果及分析

將葉片的上下弧線作圖3所示的刻度標記，以便于描述葉片表面的覆冰形貌。圖4，5分別為NACA0018和NACA6510翼型試驗組覆冰模擬時長20min后葉片表面的覆冰厚度曲線圖。

圖3 葉片弧線的標記Fig.3 Blademarking

依據(jù)所得到的覆冰厚度曲線分別繪制了NACA0018和NACA6510翼型試驗組葉片的覆冰形貌，如圖6，7所示。

圖6 NACA0018翼型試驗組葉片覆冰形貌Fig.6 Icingmorphology of the NACA0018 airfoil experimental groups

圖7 NACA6510翼型試驗組葉片覆冰形貌Fig.7 Icingmorphology of the NACA6510 airfoil experimental groups

由圖4，6可知：覆冰主要集中在葉片前緣區(qū)域；λ為0.35 的A01和A11組的葉片后緣存在少量覆冰層；在A01，A02和A03三葉片型風機中，隨著λ的增大，葉片表面的最大覆冰厚度逐漸增大；在A11，A12和A13五葉片型風機中，隨著λ的增大，葉片表面的最大覆冰厚度先增大，后減?。籄01，A03和A12組的最大覆冰厚度位于葉片尖端偏下弧線區(qū)域；A02，A11和A13組的最大覆冰厚度位于葉片尖端頂點處。

由圖5，7可知：覆冰主要集中在葉片前緣區(qū)域；λ為0.35 的B01和B11組的葉片后緣存在少量覆冰層；在B01，B02和B03三葉片型風機中，隨著λ的增大，葉片表面的最大覆冰厚度逐漸增大；在B11，B12和B13五葉片型風機中，隨著λ的增大，葉片表面的最大覆冰厚度先增大，后無顯著變化；六組葉片的最大覆冰厚度均位于葉片尖端頂點處。

采用進口多相流泵和專有微氣泡釋放技術為核心的微氣泡生成系統(tǒng)。微氣泡采用大尺寸管道釋放，可調節(jié)，不易結垢和堵塞;注氣量小，僅為處理水量的2%～3%(體積比，標態(tài)氣體)，不需要配套龐大的輔助制氮系統(tǒng);微氣泡粒徑小(5～30μm)，均勻性好，可實現(xiàn)微氣泡純物理破乳，無需化學藥劑，不產生含油污泥。

圖5 NACA6510翼型試驗組葉片覆冰厚度Fig.5 Ice thickness curve of the NACA6510 airfoil experimental groups

由圖4～7可知，當風機的λ增大時，風機葉片表面的覆冰更多地向葉片前緣尖端區(qū)域集中，葉片后緣的覆冰量減少。

3.2 覆冰質量模擬結果及分析

覆冰模擬20min后，12組風機所捕獲過冷卻水滴的數(shù)量情況如圖8，9所示。

圖8 風機整體所捕獲過冷卻水滴情況Fig.8 Overall icingmass of the VAWT

由圖8可知：NACA0018五葉片翼型風機的覆冰量最多，NACA6510三葉片翼型風機的覆冰量最少；隨著λ的增大，三葉片型風機的覆冰量先減少后增加，五葉片型風機的覆冰量先增加后減少；在相同λ和葉片數(shù)量條件下，NACA0018翼型風機的覆冰量高于NACA6510翼型風機。

由圖9可知：NACA0018三葉片翼型風機葉片的覆冰量最多，NACA6510五葉片翼型風機葉片的覆冰量最少；隨著λ的增大，三葉片型風機葉片的覆冰量先減少后增加，五葉片型風機葉片的覆冰量先增加后減少；在相同λ和葉片數(shù)量條件下，NACA0018翼型風機葉片的覆冰量高于NACA6510翼型。

圖9 單個葉片所捕獲過冷卻水滴情況Fig.9 Icingmass of single blade of the VAWT

3.3 覆冰對風機氣動性能的影響

通過葉素理論分析風機葉片的受力情況，流體對葉片的氣動作用力可分解為升力FL和阻力FD。

式中：U∞為風速；A為翼型的面積；CL，CD分別為升力系數(shù)、阻力系數(shù)，其值取決于攻角大小和葉片外形。

葉片受力分析如圖10所示。

圖10 葉片受力分析Fig.10 Force analysis of blade

則該葉片的圓周扭矩計算式為

風機動態(tài)扭矩特性的研究必須結合真型風力機的轉動數(shù)據(jù)。真型風機的實際轉動過程較為復雜，存在間歇性加減速的特征，而風機靜態(tài)扭矩T的研究相對容易，并且T也能夠在一定程度上反映風機整體覆冰后的氣動性能變化情況。因此本文僅研究風機覆冰后的T特性。根據(jù)上述研究所得到的12組覆冰風機模型，將其導入COMSOL軟件中進行流場計算，則NACA0018和NACA6510翼型覆冰風機試驗組的T與方位角φ的關系曲線分別如圖11，12所示。例如，風機為三葉片型時，風機的T為風機的3個葉片（φ，φ+120°，φ+240°）T1之和。

由圖11可知，覆冰對三葉片型、五葉片型NACA0018翼型風機的T均產生了明顯影響。在圖11（a）中，覆冰后各風機的最大T和最小T均明顯下降，同時曲線的波動振幅減小。在圖11（b）中，覆冰后各風機的最大T和最小T均明顯下降，同時曲線的波動振幅加大。

圖11 NACA0018翼型風機靜態(tài)扭矩特性曲線Fig.11 Static torque characteristic curves of the NACA0018 airfoil experimental groups

由圖12可知：覆冰對三葉片型、五葉片型NACA6510翼型風機的T均產生了一定程度影響；在圖12（a）中，覆冰后各風機的最大T和最小T出現(xiàn)不同程度的增大，同時曲線的波動振幅略微減小；在圖12（b）中，覆冰后各風機的最大T和最小T出現(xiàn)不同程度的下降，同時曲線的波動振幅不變或略微加大。

圖12 NACA6510翼型風機靜態(tài)扭矩特性曲線Fig.12 Static torque characteristic curves of the NACA6510 airfoil experimental groups

對比圖11，12可知，NACA6510五葉片翼型風機在覆冰后其T為負的區(qū)間明顯增大，這會造成旋轉風機的振動加劇，縮短風機的壽命。

4 結論

①垂直軸風機葉片的覆冰主要集中在葉片前緣區(qū)域。當風機的λ增大時，風機葉片表面的覆冰更多地向葉片前緣尖端區(qū)域集中，葉片后緣的覆冰量減少。

②不同葉片翼型、葉片數(shù)量和λ對垂直軸風機的整體覆冰量、單個葉片的覆冰量均產生重要影響。在相同λ和葉片數(shù)量下，NACA0018翼型的覆冰量高于NACA6510翼型。隨著λ的增大，三葉片型風機的整體覆冰量、單個葉片覆冰量均是先減少后增加，而五葉片型風機則是先增加后減少。

③不同葉片翼型、葉片數(shù)量和λ對風機葉片表面的最大覆冰厚度產生重要影響。隨著λ的增大，三葉片型風機葉片表面的最大覆冰厚度逐漸增大，NACA0018五葉片翼型風機葉片表面的最大覆冰厚度先增加后減少，NACA6510五葉片翼型風機葉片表面的最大覆冰厚度先增加后無顯著變化。

④NACA0018翼型風機葉片表面最大覆冰厚度的位置受不同葉片數(shù)量、λ的影響可能發(fā)生移動，而NACA6510翼型風機葉片表面最大覆冰厚度的位置始終位于葉片尖端頂點處。

⑤不同葉片翼型、葉片數(shù)量和λ的風機的靜態(tài)扭矩特性受覆冰影響差異大。