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

（長(zhǎng)沙理工大學(xué) 電氣與信息工程學(xué)院，湖南 長(zhǎng)沙410114）

0 引言

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

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

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

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

1 數(shù)學(xué)模型

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

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

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

1.2 過(guò)冷水滴的運(yùn)動(dòng)軌跡

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

1.3 覆冰質(zhì)量計(jì)算

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

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

2 風(fēng)機(jī)覆冰模擬

2.1 風(fēng)機(jī)覆冰物理模型

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

圖1 風(fēng)機(jī)覆冰物理模型Fig.1 Physicalmodel of the VAWT icing

2.2 試驗(yàn)組參數(shù)

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

圖2 葉片翼型Fig.2 Blade airfoils

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

續(xù)表1

3 覆冰模擬結(jié)果及分析

3.1 冰型模擬結(jié)果及分析

將葉片的上下弧線(xiàn)作圖3所示的刻度標(biāo)記，以便于描述葉片表面的覆冰形貌。圖4，5分別為NACA0018和NACA6510翼型試驗(yàn)組覆冰模擬時(shí)長(zhǎng)20min后葉片表面的覆冰厚度曲線(xiàn)圖。

圖3 葉片弧線(xiàn)的標(biāo)記Fig.3 Blademarking

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

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

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

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

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

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

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

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

3.2 覆冰質(zhì)量模擬結(jié)果及分析

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

圖8 風(fēng)機(jī)整體所捕獲過(guò)冷卻水滴情況Fig.8 Overall icingmass of the VAWT

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

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

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

3.3 覆冰對(duì)風(fēng)機(jī)氣動(dòng)性能的影響

通過(guò)葉素理論分析風(fēng)機(jī)葉片的受力情況，流體對(duì)葉片的氣動(dòng)作用力可分解為升力FL和阻力FD。

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

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

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

則該葉片的圓周扭矩計(jì)算式為

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

由圖11可知，覆冰對(duì)三葉片型、五葉片型NACA0018翼型風(fēng)機(jī)的T均產(chǎn)生了明顯影響。在圖11（a）中，覆冰后各風(fēng)機(jī)的最大T和最小T均明顯下降，同時(shí)曲線(xiàn)的波動(dòng)振幅減小。在圖11（b）中，覆冰后各風(fēng)機(jī)的最大T和最小T均明顯下降，同時(shí)曲線(xiàn)的波動(dòng)振幅加大。

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

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

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

對(duì)比圖11，12可知，NACA6510五葉片翼型風(fēng)機(jī)在覆冰后其T為負(fù)的區(qū)間明顯增大，這會(huì)造成旋轉(zhuǎn)風(fēng)機(jī)的振動(dòng)加劇，縮短風(fēng)機(jī)的壽命。

4 結(jié)論

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

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

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

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

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