垂直軸風力機在不同運動狀態下的電加熱融冰策略研究

張漢文，羅日成，劉志勇

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

0引言

近年來我國中南、西南等地新建了多個風電基地，但這些地區春冬季節的濕冷氣候使得風力機面臨著嚴峻的覆冰問題[1]，[2]。覆冰后，葉片氣動性能下降，功率輸出降低，荷載增大造成主軸磨損加劇，以及其它安全事故[3]。目前，風力機防除冰技術主要有涂料、熱空氣、電加熱、電磁脈沖和超聲波振動等技術，其中電加熱被認為是最有效和最經濟的技術[4]。國內外相關學者在風力機電加熱防除冰技術的研究上已取得了一定成果。文獻[5]通過ANSYS軟件對風力機葉片表面電熱元件的布局形式、幾何形狀進行了研究，結果表明，圓形加熱片除冰速度更快，交錯布局可以降低葉片表面的熱應力。文獻[6]對風力機葉片表面的電阻絲布置方式進行了試驗研究。文獻[7]對風力機電加熱系統的控制理論進行了研究，通過仿真分析了周期性加熱對防除冰的影響。文獻[8]提出了一種風力機葉片循環控制加熱的除冰方法，并利用數值仿真進行了驗證。

垂直軸風力機具有低成本、占地空間小、無需偏航裝置、易于維護和安全性高等優點[9]。隨著垂直軸風力發電技術的不斷進步，已有一些技術較為完善的兆瓦級垂直軸風力機推向風電市場。現有的風力機電加熱防除冰技術主要是基于水平軸風力機的結構特點進行研發的，而垂直軸風力機的結構存在較大不同，在氣流中進行電加熱融冰時，兩者的對流散熱存在較大不同。因此，本文通過分析垂直軸風力機的空氣流場和葉片表面冰層的融冰過程，依據葉片在不同攻角下的融冰速度差異和葉片兩側的融冰速度差異，提出了垂直軸風力機在不同運動狀態下的電加熱融冰策略。該策略可以加快風力機的融冰速度，減少電能損失，具有一定的工程應用價值。

1 數理模型及控制方程

在風力機電加熱融冰過程中，電熱單元產生的焦耳熱逐漸向四周傳遞，附近冰層的溫度緩慢升高。一段時間后，冰層表面與周圍氣流間產生了溫度差，冰層表面出現對流散熱和輻射散熱。此外，當靠近葉片表面的冰層的溫度上升到0℃時，接觸面的冰層將逐漸融化，并產生空氣間隙[10]。當一定量的內部冰層被融化后，剩余冰層與葉片表面間的附著力極大降低，剩余冰層將極易在重力、離心力等作用下發生脫落。

基于上述物理過程，本文建立了風力機的電加熱融冰數值模型，分別對空氣流場、熱量傳遞過程以及相變過程進行計算。

1.1 空氣流場分布

不考慮外界氣流的波動以及散熱項對葉片周圍氣流流動的影響，本文假定風力機周圍的氣流場為穩態分布。空氣流場的控制方程選用低速粘流的雷諾平均N-S方程，湍流模型選用對邊界層湍流和自由剪切湍流均有較好模擬效果的SST k-ω湍流模型，該模型能夠較好地還原風力機周圍的空氣流場分布[11]。輸運變量為湍流動能k和比耗散率ω，其輸送方程為

式中：u為氣流速度場；ρ為空氣密度；μ為空氣動力粘度；t為時間；Pk為k的有效生成率；Pω為ω的有效生成率；β，β*，σk，σω和σω2均為經驗系數；F1為混合函數。

1.2 傳熱計算

在融冰過程中，電熱單元產生的焦耳熱主要有3個去向：①葉片、冰層吸收熱量溫度升高；②冰層表面與周圍環境間的對流散熱和輻射散熱；③部分冰層吸收熱量相變融化。

1.2.1溫度場分布

葉片、冰層的傳熱形式屬于固體域傳熱，周圍氣流的傳熱形式屬于流體域傳熱，這兩種傳熱形式存在一定差異。

周圍氣流的溫度場可由流體域內的時變傳熱微分方程描述。

式中：Tf為流體域的溫度場；Vf為流體域的速度場；Cp，λ分別為空氣流的比熱容和熱導系數。

在不考慮部分冰層相變融化的情形下，葉片、冰層的溫度場可由固體域內的時變傳熱微分方程描述。

式中：TS為固體域的溫度場；qin為固體域內的凈熱源。

當接觸面的冰層融化成液態水時，該區域的溫度場表達式為

式中：▽mice為融化冰層的質量；Lm為冰的融化潛熱，Lm=330 kJ/kg；qnet為該融化區域內的凈熱源。

1.2.2對流散熱項與輻射散熱項

當冰層表面的溫度為TW，環境的溫度為Ta時，冰層表面的對流散熱項Δqf和輻射散熱項Δqr分別為

式中：TW，Ta分別為冰層表面和環境的溫度；hc為冰層表面的對流換熱系數，其值與局部空氣流場有關；εi為冰層表面輻射率，取0.95；σ為斯特藩常量，取5.67×10-8W/（m2·K4）；S1為冰層表面與周圍環境的接觸面積。

1.2.3冰水相變過程

實際的冰水相變過程中，冰水混合物的溫度會一直保持在0℃，而數值計算模型會假設冰水相變過程發生在0℃附近一個較小的溫度區間，否則無法求解。本文使用顯熱容法來描述冰水相變過程（圖1）。

圖1 冰水相變數值模型Fig.1 Ice/water phase change numerical model

圖1中：體積分數θ1，θ2分別為冰和液態水的含量；冰相和水相之間的相變溫度Tpc，1→2設為0℃；轉變間隔ΔT1→2設為0.001℃。相變過程中潛熱L1→2設為330 kJ/kg，相變過程中涉及的一些等效物理參數計算如下：

式中：ρeq，Cp_eq，λeq分別為等效密度、等效熱容和等效導熱系數。

2 風力機電加熱融冰模型

2.1 幾何模型

文獻[12]對NACA0018翼型垂直軸風力機進行了覆冰研究，得到了多種覆冰工況下風力機葉片表面冰層的形貌和增長規律。這些不同覆冰工況下風力機葉片表面冰層的形貌均較為相似，覆冰均主要集中在葉片前緣區域。其中一種覆冰工況下的葉片表面冰層如圖2所示。

圖2 葉片表面冰層的形貌和增長規律Fig.2 Morphology and growth law of ice on blade surface

依據上述垂直軸風力機葉片表面冰層的形貌特點，通過多物理場耦合仿真軟件COMSOL建模，本文的仿真模型如圖3所示。

圖3 風力機電加熱融冰幾何模型Fig.3 VAWT's electric heating ice-melting geometric model

覆冰葉片由葉片自身、電熱元件和冰層（分為外冰層和內冰層，內冰層又分為內冰層_outside和內冰層_inside）構成。整個計算域尺寸為10 m×8 m，風力機直徑為3.6 m，葉片弦長為0.5 m，葉片翼型為NACA0018。模型存在如下假設條件：傳熱過程不影響葉片附近的空氣流場；忽略冰水相變過程中兩者密度的變化。

2.2 參數設置

融冰模型為二維模型，面外厚度為0.2 m，所有葉片的熱功率均設為60 W，對應的熱功率密度為580 W/m2。環境參數：來流風速為5 m/s，氣流溫度為-5℃，各物體的初始溫度均為-5℃，液態水含量（LWC）為0 g/m3，壓力出口為0 MPa，內冰層設為相變材料。模型中各材料的物理參數見表1。其中外冰層最厚處為14 mm，最薄處為2 mm。模型網格劃分如圖4所示。

表1 各材料的物理參數Table 1 Physical parameters of each material

圖4 網格劃分Fig.4 Mesh division of model

3 空氣流場分布

通過“凍結轉子”求解器得到風力機在λ為0和0.2兩種葉尖速比下的空氣流場（圖5）。由圖5可知，當風力機處在不同方位時，葉片附近區域的空氣流場存在差異，風力機處于靜止狀態（λ=0）和處于低轉速狀態（λ=0.2）的空氣流場存在一定程度的相似。

圖5 不同方位下風力機周圍的氣流流場分布Fig.5 In different directions，distribution of air flow field around the VAWT

以V外，V內分別計算葉片表面兩側的平均風速，在λ為0和0.2兩種葉尖速比下，葉片兩側的V外，V內變化曲線如圖6所示。

圖6 葉片兩側風速的變化規律Fig.6 Variation law of wind speed on both sides of blade

從圖6可知：當葉片位于迎風側時，葉片外表面的空氣流場流速高于葉片內表面；當葉片位于背風側時，葉片外表面的空氣流場流速多數情況下低于葉片內表面。

4 風力機處于靜止狀態下的融冰策略研究

4.1 溫升融化過程分析

加熱640 s后，不同方位和不同α下風力機周圍的溫度場分別如圖7，8所示。為了方便比較，所有溫度場云圖均采用相同標尺。

圖7 不同方位下風力機周圍的溫度場分布Fig.7 In different directions，distribution of temperature field around the VAWT

圖8 不同α下葉片周圍的溫度場分布Fig.8 Temperature distribution around blades at different angles of attack

從圖7，8可知：不同α下葉片周圍的溫度場存在較大差異，葉片前緣和葉片尾緣存在一條明顯的溫升帶；當α分別為240°和300°時，葉片的背風側均存在明顯的升溫區域，這是由于在加熱過程中，通過對流散熱和輻射散熱從葉片表面所釋放出來的熱量在葉片附近這些流速較低的區域積聚較多。

依據文獻[10]，本文假設當內冰層融化90%時，認為葉片表面剩余冰層可在重力、離心力等作用下脫落。從開始加熱至內冰層即將開始融化所經歷時長為內冰層的升溫期，從內冰層開始融化至融化90%所經歷時長為內冰層的90%融化期。不同方位下風力機葉片內冰層的受熱過程如圖9所示。不同α下風力機葉片內冰層的受熱過程如圖10所示。

圖9 不同方位下風力機的融冰情況Fig.9 Ice-melting of the VAWT in different directions

圖10 在不同α下葉片的融冰情況Fig.10 Ice-melting of blade at different angles of attack

從圖9可知，由于風力機的融冰時間取決于該風力機中所需加熱時間最長的那個葉片，則C方位下風力機的融冰速度最快，各葉片間的融冰速度差異也小，其次為B方位和D方位，A方位下風力機的融冰速度最慢。從圖10可知：α靠近0°的葉片的融冰速度較慢，當葉片位于迎風側時，葉片內表面的融冰速度快于葉片外表面；而當葉片位于背風側時，葉片外表面的融冰速度多數情況下快于葉片內表面，這是由于葉片兩側的對流散熱情況不同所導致的。

4.2 融冰策略分析

當對在氣流中處于靜止狀態下的垂直軸風力機進行電加熱融冰時，可采用下述幾種措施：①安裝特定機械裝置，利用該裝置將風力機調整至最佳融冰方位，C方位下風力機融冰速度最快，A方位下風力機融冰速度最慢；②α靠近0°的葉片的融冰速度較慢，通過控制系統來適當增大這些α下的葉片的電熱功率；③當葉片處于迎風側或背風側時，葉片表面兩側的融冰速度存在差異，并且這種差異會隨著風速的增大、環境溫度的下降而變大，通過控制系統來適當調節葉片表面兩側的電熱功率來縮短葉片兩側的融冰差異。

5 風力機處于低速旋轉狀態下的融冰策略研究

5.1 融冰策略分析

當風力機處于旋轉狀態時，葉片兩側會周期性交替出現高速、低速氣流場，將持續影響冰層表面的散熱項。本文提出一種電熱功率動態調節方式，通過跟蹤風力機中各葉片α的變化來動態調節葉片表面兩側的電熱功率：當葉片處于迎風側時，減少葉片外表面的對流散熱損耗；當葉片處于背風側時，減少葉片內表面的對流散熱損耗。

5.2 建立等效模型

考慮到在旋轉狀態下葉片兩側的空氣流場變化情況，圖11為等效簡化模型。利用兩個風速函數分別控制入口A和入口B的氣流場，以5 m/s的風速下風力機葉尖速比為0.2為例（圖6），高速氣流設為5 m/s，低速氣流設為0.2 m/s，周期為T，前T/2周期內葉片表面A側為高速氣流，B側為低速氣流，后T/2周期內A側為低速氣流，B側為高速氣流。

圖11 等效模型Fig.11 The equivalent model

5.3 仿真分析

葉片兩側的電熱功率分別為PA和PB，存在兩種功率調節方式：①當風速高時，電熱功率減弱，當風速低時，電熱功率加強；②當風速高時，電熱功率加強，當風速低時，電熱功率減弱。以第一種功率調節方式融冰時，在前T/2周期內，A側為低電熱功率，B側為高電熱功率，在后T/2周期內，A側為高電熱功率，B側為低電熱功率；以第二種功率調節方式融冰時，A，B兩側電熱功率與上述情況相反。加入恒功率組（PA，PB均為75 W）作為對照組，第一種功率調節方式和第二種功率調節方式的PA，PB分別在75 W的基礎上進行±25%的調整，且任意時刻PA+PB=150 W。

研究風力機在λ分別為0.2，0.1，0.05，0.025，0.012 5和0.006 25下的融冰過程，環境參數：高速氣流為5 m/s，低速氣流為0.2 m/s，氣流溫度為-8℃，各物體的初始溫度均設為-8℃，LWC為0 g/m3，壓力出口為0 MPa。

圖12為3種電熱功率調節方式下葉片的融冰情況。從圖12可知：當λ分別為0.2和0.1時，第二種功率調節方式的融冰速度均最快，第一種功率調節方式的融冰速度均明顯慢于恒功率方式；當λ為0.05時，第二種功率調節方式的融冰速度依舊最快，第一種功率調節方式快于恒功率方式；當λ分別為0.025，0.012 5和0.006 25時，3種功率調節方式的融冰速度均較為相近，恒功率方式具有微弱優勢。

圖12 3種電熱功率調節方式下葉片的融冰情況Fig.12 Ice-melting of blade under three electrothermal power regulation modes

6 結論

①風力機處于靜止狀態（λ=0）和處于低轉速狀態（λ=0.2）的空氣流場存在一定程度的相似。

②對在氣流中處于靜止狀態的風力機進行電加熱融冰時，葉片在不同α下的融冰速度存在明顯差異，其中α靠近0°時的葉片融冰速度明顯較慢；當葉片位于迎風側或背風側時，葉片兩側的融冰速度存在差異，當葉片位于迎風側時，葉片內表面的融冰速度快于葉片外表面，當葉片位于背風側時，葉片外表面快于葉片內表面。因此，安裝機械裝置對風力機進行方位調整，通過控制系統來適當調節特定攻角的葉片的電熱功率和葉片兩側的電熱功率分配可以提升融冰速度。

③基于垂直軸風力機在旋轉狀態下的葉片表面兩側的空氣流場差異提出了一種電熱功率動態調節方式，在該策略中，當葉片位于迎風側時，減少葉片外表面的對流散熱損耗；當葉片位于背風側時，減少葉片內表面的對流散熱損耗。利用所建立的等效模型仿真發現，當λ為0.05～0.2時，所提出的電熱功率動態調節方式的融冰速度快于恒功率方式。