大型風力機結冰過程水滴收集率三維計算

易 賢，王開春，馬洪林，朱國林

（中國空氣動力研究與發展中心 空氣動力學國家重點實驗室，四川 綿陽 621000）

0 引 言

風力機在寒冷的氣候（溫度低于0℃）環境下運行時，如果遇到含有過冷水滴的氣流，或者遭遇雨雪天氣，在風力機表面（包括葉片、輪轂、機艙、塔架、風速風向儀等風力機所有部件表面）會發生結冰現象。結冰會改變葉片的氣動外形，破壞風力機的流場特性，使風能利用系數降低，減小出力；同時，結冰還會引起額外的過載及振動，冰載荷的不平衡會增加機組部件的疲勞，輕者，使風力機停止運行，帶來經濟損失，重者，造成葉片損壞甚至導致整個風機倒塌，危及人們的生命財產安全。鑒于風力機結冰的危害，近年來人們越來越重視這方面的研究［1］。

風力機結冰包括沉積結冰（Precipitation Icing，主要由凍雨和濕雪引起）和過冷水滴結冰（In－Cloud Icing，主要由過冷水滴撞擊在風力機表面所引起）兩種形式［2］，其中過冷水滴結冰是低溫氣象條件下常見的一種結冰形式。水滴收集率（也稱水滴局部收集系數）定義為風力機某局部區域實際所收集的水量與該區域可能收集的水量最大值之比，它表征了部件表面的水滴撞擊范圍以及撞擊區域內水量的分布。水滴收集率是結冰研究中最重要的參數之一，只有獲得水滴收集率，才能進行結冰預測、探測以及防／除冰系統設計。對于過冷水滴結冰，常用的研究手段有工程估算、結冰風洞試驗和數值計算。其中試驗研究是早期主要采用的方法［3－4］，近年來，開始陸續出現用數值計算的手段開展研究的報道［1，5－7］，由于旋轉條件下風力機流場的復雜性，這些計算多采用拉格朗日法計算水滴收集率，或者采用葉片的截面翼型進行二維簡化計算。本文在前期發展的飛機結冰三維計算方法基礎上［8－10］，建立了適合于大型風力機結冰過程水滴收集率計算的三維數值方法和計算程序，并對某1.5MW級水平軸風力機的水滴收集率進行了計算，獲得了不同直徑水滴在風力機表面的撞擊特征以及水滴收集率在葉片上的分布規律，為開展風力機結冰防護研究打下了較好基礎。

1 水平軸風力機空氣流場計算方法

風力機的空氣繞流場是典型的包含旋轉機械的流場，對于這類流場計算，有單旋轉坐標系方法、多參考坐標系（Multiple Reference Frame，簡稱 MRF）方法［11］、滑移網格方法、動態重疊網格方法等，其中，當流場中同時存在旋轉部件和非旋轉部件（如轉子／定子）時，MRF方法由于其魯棒、簡便、高效的特點，在工程中得到了廣泛的使用，一些CFD商業軟件如FLUENT、STAR－CD、CFX等，都集成了該方法。

MRF方法的基本思路為：將計算區域分為旋轉區和非旋轉區，在非旋轉區，流動的控制方程為非旋轉坐標系（如慣性坐標系）下的方程，在旋轉區域，流動的控制方程為旋轉坐標系下的控制方程。對不同區域的流動分別進行求解，在區域的交界面上，要求流場的參數完全一致。

本文在原有的非旋轉流場計算研究的基礎上，建立了計算風力機流場的MRF方法，現概述如下。

1.1 空氣流場計算的控制方程

非旋轉區域空氣流場計算的控制方程為低速粘流的時均N－S方程，其通用形式為：

方程（1）中，φ為輸運變量，ρa為空氣密度，ua為空氣速度，Γφ為擴散系數，qφ為源項，φ、Γφ和qφ取不同的值，可代表空氣的連續性方程、動量方程和其他標量（如湍動能等）的輸運方程，方程中各項的物理意義及具體表達式可參見文獻［8］，此處不再詳述。

在旋轉區域內，控制方程也可以寫成輸運方程的形式：

方程（2）中，uar為空氣在旋轉坐標系中的速度，qφr為旋轉坐標系中的源項。當φ＝1、qφr＝0時，方程（2）為連續性方程，當φ＝uar時，方程（2）為動量方程，此時源項的表達式為：

其中，ω為旋轉角速度，r為流體質點在旋轉坐標系中的矢徑。

在旋轉區域和非旋轉區域的交界面上，要求流場參數一致，兩個坐標系內的速度關系為：

1.2 控制方程的有限體積離散

采用有限體積法離散求解方程（1）、（2）。以點P為中心點的六面體為例，方程中各項的離散方法為：

（1）對流項的離散。

式（5）中，nb＝e，w，n，s，t，b，分別代表以 P點為中心的控制體的六個面，Snb為各面對應的面積，對于空氣相的控制方程有：

控制體界面上的變量值fnb用迎風插值和線性插值相結合的方法計算，表達式為：

其中為迎風插值的量，為線性插值的量，ε為混合因子，0≤ε≤1。

（2）源項的離散。輸運方程的源項對不同的方程都有不同的表達式，為了使其離散化表達式盡可能逼近源項本身，加強代數方程的主對角優勢，提高代數方程組求解的穩定性，對源項采用線性化方式做進行處理，即令

則源項的離散形式為：

式中δV為控制體的體積。

（3）時間項離散。方程（1）、（2）的時間項采用如下二階精度的隱式離散：

其中，上標n＋1、n和n－1分別代表t＋Δt時刻、t時刻和t－Δt時刻的值。

2 水滴項的控制方程及求解方法

水滴收集率計算在空氣流場計算完成之后進行，主要有兩種方法［12－13］，一種是拉格朗日法，即用拉格朗日的觀點，計算跟蹤流場中單個水滴的運動，并判斷其是否與物面碰撞，從而得到物面各處的水滴收集量及其分布；另一種是歐拉法，該方法將水滴看作與空氣單向作用的、在空間連續分布的液態相，通過流體力學的觀點求解水滴相的控制方程，得到流場中水滴容積分數分布，進而獲得水滴收集率。鑒于拉格朗日法在三維計算中的不足，本文發展了三維水滴收集率計算的歐拉方法。

2.1 水滴相流場計算的控制方程

水滴相方程的求解統一在慣性坐標系中進行，引入水滴容積分數α，其定義為空間微團中水滴相所占的體積比例，則可以建立水滴相的控制方程，包括連續方程和動量方程，分別為［8］：

方程（12）、（13）中，ud為水滴速度，ρd為水滴密度，g為重力加速度，K為慣性因子，其表達式為：

式（14）中，μa為空氣動力粘性系數，d為水滴直徑，CD為水滴阻力系數，Re為相對雷諾數，其表達式為：

水滴運動過程中，阻力隨相對雷諾數變化而變化，本文采用如下公式確定阻力［8］：

與空氣控制方程類似，水滴相控制方程可以統一寫成輸運方程的形式：

式中，qφ為源項，φ取1、ud、vd或wd分別代表連續方程和x、y、z方向的動量方程。

2.2 水滴相控制方程的離散方法

對于方程（17），對流項和源項的離散方法與方程（1）、（2）一致，時間項采用一階顯式離散，即：

2.3 壁面邊界條件

與空氣流場計算在壁面采用無滑移壁面邊界條件不同，對于水滴相計算，采用壁面吸入邊界條件，即如果水滴與物面碰撞，則認為水滴從碰撞點流出。

3 水滴收集率的計算

由于水滴容積分數α較?。?0－6量級），可以認為空氣和水滴是單向作用，即只考慮空氣對水滴的作用，忽略水滴對空氣的作用。因此，水滴收集率的計算步驟可概括為：首先，計算空氣場，本文采用SIMPLE方法計算空氣場，湍流模型為標準k－ε湍流模型［14］；其次，在得到空氣流場分布的基礎上，求解水滴相控制方程；最后，水滴收集率β可在獲得當地水滴容積分數α和水滴與物面的相對速度udr之后，由以下表達式得到：

其中，α∞為遠場水滴容積分數，u∞為遠場來流速度，n為物面碰撞點處的單位法線向量。

由于udr＝ud－ω×r，式（19）可以寫成：

從式（20）可以看出，旋轉狀態下的水滴收集率β，由兩部分組成，一部分是風力機旋轉速度為0時的水滴收集率β0，另一部分是由于旋轉而產生的水滴收集率βr。在慣性坐標系下，水滴收集率的物理意義定義為：物面某局部區域所收集的水量與該區域所收集的水量最大值之比，因此β0值應該在［0，1］區間，而βr的值則與物面某區域的旋轉半徑有關，半徑越大，βr的值也越大。因此，在旋轉條件下，由式（19）計算得到的β值將會出現大于1的情況。

4 大型風力機水滴收集率計算

根據以上數值方法，開發了相應的計算程序，并針對某1.5MW級的水平軸風力機進行了應用計算。

4.1 計算外形及網格

計算外形選用某1.5MW級的水平軸風力機，半徑R＝41m，忽略塔架的影響。如圖1，坐標原點取為風力機輪轂中心，坐標軸方向為：x軸與遠場來流方向一致，y軸垂直向上，z軸按右手系確定。計算采用多塊對接網格，網格規模為1800萬。圖2顯示了風力機表面輪轂附近的網格。

圖1 計算外形及坐標方向Fig.1 Wind turbine configuration and frame

圖2 輪轂附近表面網格 Fig.2 Surface grid near the hub

4.2 風力機空氣流場計算結果

本文計算了兩種來流條件下風力機的氣動特性，并與設計值進行了對比，計算條件為：

計算條件1：來流速度8m／s，風輪轉速16.7r／min，漿距角0°，葉尖速比為9。

計算條件2：來流速度11m／s，風輪轉速17.4r／min，漿距角3.2°，葉尖速比為6.8。

圖3給出了速度為8m／s時葉片上的壓力計算結果。圖3（a）給出了葉片上壓力分布云圖，圖3（b）～圖3（f）顯示的是葉片上R＝4m、13m、22m、31m和40m五個不同位置截面上的壓力系數曲線。計算顯示，在靠近輪轂附近（R＝4m），由于葉片截面還不是翼型形狀，迎風面和背風面的壓力呈現相反的變化趨勢，在30%弦長處，吸力面在迎風和背風面之間發生轉換；葉片中部（R＝22m）一直到葉尖附近（R＝40m），壓力分布呈相似的曲線輪廓。

圖3 葉片表面壓力分布（v＝8m／s）Fig.3 Pressure distribution on the blade（v＝8m／s）

表1給出了兩種速度條件下，本文計算的軸功率與設計值的比較，可見本文計算的功率與設計值接近，計算值比設計值只分別低0.6%和1.3%。

表1 計算的軸功率與設計值的對比Table 1 Comparison of computed power and designed power

采用固定8m／s的來流速度、改變風輪轉速的方式，計算了不同葉尖速比條件下風輪的氣動特性。圖4（a）～圖4（c）分別給出了λ＝9、4和2三種葉尖速比條件下，吸力面上的極限流線分布情況。圖中顯示，葉尖速比λ＝9時，吸力面的流動為附著流，未出現分離，整個葉片表面近似為二維弦向流動，葉尖速比λ＝4時，在葉片的中部首先出現分離，當葉尖速比降至2時，整個吸力面上的流動已完全分離。

圖4 吸力面極限流線分布Fig.4 Limiting streamlines distribution on suction surface

4.3 風力機水滴收集率計算結果

在風力機空氣流場計算的基礎上，開展了水滴收集率計算，對比了不同直徑水滴的收集率在風力機表面的分布情況。計算條件為：來流速度8m／s，密度1.225kg／m3，環境壓力101325Pa，風輪轉速16.7r／min，漿距角0°，葉尖速比為9，水滴直徑選取了三種，從大到小分別為100μm、40μm和20μm。

圖5～圖8給出了水滴直徑d＝100μm時的計算結果。圖5為風力機表面的水滴收集率分布云圖，由圖可見，水滴主要撞擊在輪轂表面和葉片前緣；圖6為收集率云圖在輪轂附近的放大，圖7為單個葉片表面的分布云圖，圖8給出的是葉片表面不同展向位置的水滴收集率曲線。計算顯示：（1）從葉片根部到葉尖，水滴收集率逐漸增大，葉尖附近（R＝40m）的最大收集率達到7.5以上，而葉片根部附近（R＝4m）的收集率僅為0.2左右；（2）輪轂上的水滴收集率極小（不到0.1），小于葉片根部的收集率。

根據以上水滴收集率計算結果，可以對輪轂表面的結冰厚度進行估算，物體表面某處結冰的厚度h可以表達為［8］：

圖5 風力機表面水滴收集率分布云圖（d＝100μm）Fig.5 Contours of collection efficiency（d＝100μm）

圖6 輪轂附近的水滴收集率分布云圖（d＝100μm）Fig.6 Collection efficiency near hub（d＝100μm）

圖7 葉片表面的水滴收集率分布云圖（d＝100μm）Fig.7 Collection efficiency on blade（d＝100μm）

圖8 葉片上不同位置水滴收集率曲線（d＝100μm）Fig.8 Collection efficiency distribution on blade（d＝100μm）

其中f為凍結系數（0≤f≤1），β為當地的水滴收集率，LWC為空氣中的液態水含量［15］，V為來流速度，dt為結冰時間，ρi為冰的密度。

根據式（20），輪轂表面β的最大值取為0.1，取f為1，V取為8m／s，LWC取0.5g／m3，ρi取700kg／m3，則1h的結冰厚度約為0.2cm，連續50h的結冰厚度也僅為10cm。由于結冰主要影響風力機的氣動特性和葉片的載荷強度及疲勞特性，因此在實際分析時，可以忽略輪轂的結冰，主要考慮葉片結冰的影響。

圖9、圖10給出了水滴直徑d＝40μm時，風力機表面的水滴收集率分布；圖11、圖12給出的是水滴直徑為20μm的結果?？梢钥闯?，水滴直徑變小之后，水滴撞擊特性的變化規律為：（1）水滴收集率隨水滴直徑變小而減??；（2）不同直徑水滴在風力機表面的撞擊區域以及水滴收集率在葉片上的分布特征保持類似。

圖9 輪轂附近的水滴收集率分布云圖（d＝40μm）Fig.9 Collection efficiency near hub（d＝40μm）

圖10 葉片上不同位置水滴收集率曲線（d＝40μm）Fig.10 Collection efficiency distribution on blade（d＝40μm）

圖11 輪轂附近的水滴收集率分布云圖（d＝20μm）Fig.11 Collection efficiency near hub（d＝20μm）

圖12 葉片上不同位置水滴收集率曲線（d＝20μm）Fig.12 Collection efficiency distribution on blade（d＝20μm）

以上結果說明對于過冷水滴在風力機葉片表面的結冰，其基本規律有兩個，一是冰主要結在葉片前緣，二是葉尖結冰遠比葉根處嚴重。因此，在進行結冰分析、設計結冰探測方法和防除冰方案時，葉片前緣，尤其是葉尖附近的葉片前緣，需要重點關注。

5 結 論

建立了適合于大型風力機結冰過程水滴收集率計算的三維數值方法和計算程序，并對某1.5MW級水平軸風力機的水滴收集率進行了計算，得到如下結論：

（1）本文的工作基于前期發展的飛機結冰三維計算方法，旋轉和非旋轉區的空氣流場、水滴相的控制方程均采用統一的有限體積離散方法，方便了計算程序的編制，對風力機流場和水滴收集率的計算結果初步表明本文發展的算法以及相應的計算程序是有效的。

（2）對某1.5MW級水平軸風力機結冰過程水滴收集率的計算表明，水滴主要撞擊在輪轂表面和葉片前緣，輪轂上的水滴收集率極小，從葉片根部到葉尖，水滴收集率逐漸增加，且葉尖處收集率遠比葉根處高。

（3）在進行結冰分析、設計結冰探測方法和防除冰方案時，可以忽略輪轂和葉片根部附近結冰的影響，葉片前緣，尤其是葉尖附近的葉片前緣，需要重點關注。

（4）風力機結冰過程涉及流動、傳質、傳熱和相變等多種物理現象的耦合，對該問題進行三維數值模擬是一項復雜的工作，本文發展的風力機流場和水滴收集率求解方法，為進行三維結冰計算打下了較好的基礎。

［1］FORTIN G ，PERRON J.Wind turbine icing and deicing［R］.AIAA－2009－0274，2009

［2］FROHBOESE P，ANDERS A.Effects of icing on wind turbine fatigue loads［J］.JournalofPhysics：ConferenceSeries，2007，（75）：012061.

［3］JASINSKI W J，NOE S C，SELIG M S，et al.Wind turbine performance under icing conditions［R］.AIAA－1997－0977，1997.

［4］HOCHART C，FORTIN G，PERRON J.Icing simulation of wind－turbine blades［R］.AIAA－2007－1373，2007.

［5］MAKKONEN L，LAAKSO T，MARJANIEMI M，et al.Modeling and prevention of ice accretion on wind turbines［J］.WindEngineering，2001，25（1）：3－21.

［6］BATTISTI L，FEDRIZZI R，RIALTI M，et al.A model for the design of hot－air based wind turbine ice prevention system［A］.Proceedings of the WREC［C］.Aberdeen，2005.

［7］朱程香，王瓏，孫志國，等.風力機葉片翼型的結冰數值模擬研究［J］.空氣動力學學報，2011，29（4）：522－528.

［8］易賢，王開春，桂業偉，等.結冰面水滴收集率歐拉計算方法研究及應用［J］.空氣動力學學報，2010，28（5）：596－601.

［9］易賢，桂業偉，朱國林.飛機三維結冰模型及其數值求解方法［J］.航空學報，2010，31（11）：2152－2158.

［10］易賢，桂業偉，朱國林，杜雁霞.運輸機翼型結冰的計算和實驗［J］.航空動力學報，2011，26（4）：808－813.

［11］LUO J Y，ISSA R I，GOSMAN A D.Prediction of impeller－induced flows in mixing vessels using multiple frames of reference［J］.InIChemESymposiumSeries，1994，136：549－556.

［12］BRAGG M B，GREGOREK G M，SHAW R J.An analytical approach to airfoil icing［R］.AIAA－81－0403，1981.

［13］BOURGAULT Y，HABASHI W G，DOMPIERRE J，et al.An eulerian approach to supercooled droplets impingement calculations［R］.AIAA－97－0176，1997.

［14］ZHU G L，KRONAST M.The calculation of ground effect on a car flow field using two dimensional Navier－Stokes equations［J］.ACTAAerodynamicSinica，1993，11（1）：35－40.

［15］易賢.飛機積冰的數值計算與積冰試驗相似準則研究［D］.綿陽：中國空氣動力研究與發展中心，2007.