大型風(fēng)力機(jī)結(jié)冰過程水滴收集率三維計算

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

（中國空氣動力研究與發(fā)展中心 空氣動力學(xué)國家重點(diǎn)實(shí)驗室，四川 綿陽 621000）

0 引 言

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

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

1 水平軸風(fēng)力機(jī)空氣流場計算方法

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

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

本文在原有的非旋轉(zhuǎn)流場計算研究的基礎(chǔ)上，建立了計算風(fēng)力機(jī)流場的MRF方法，現(xiàn)概述如下。

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

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

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

在旋轉(zhuǎn)區(qū)域內(nèi)，控制方程也可以寫成輸運(yùn)方程的形式：

方程（2）中，uar為空氣在旋轉(zhuǎn)坐標(biāo)系中的速度，qφr為旋轉(zhuǎn)坐標(biāo)系中的源項。當(dāng)φ＝1、qφr＝0時，方程（2）為連續(xù)性方程，當(dāng)φ＝uar時，方程（2）為動量方程，此時源項的表達(dá)式為：

其中，ω為旋轉(zhuǎn)角速度，r為流體質(zhì)點(diǎn)在旋轉(zhuǎn)坐標(biāo)系中的矢徑。

在旋轉(zhuǎn)區(qū)域和非旋轉(zhuǎn)區(qū)域的交界面上，要求流場參數(shù)一致，兩個坐標(biāo)系內(nèi)的速度關(guān)系為：

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

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

（1）對流項的離散。

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

控制體界面上的變量值fnb用迎風(fēng)插值和線性插值相結(jié)合的方法計算，表達(dá)式為：

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

（2）源項的離散。輸運(yùn)方程的源項對不同的方程都有不同的表達(dá)式，為了使其離散化表達(dá)式盡可能逼近源項本身，加強(qiáng)代數(shù)方程的主對角優(yōu)勢，提高代數(shù)方程組求解的穩(wěn)定性，對源項采用線性化方式做進(jìn)行處理，即令

則源項的離散形式為：

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

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

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

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

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

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

水滴相方程的求解統(tǒng)一在慣性坐標(biāo)系中進(jìn)行，引入水滴容積分?jǐn)?shù)α，其定義為空間微團(tuán)中水滴相所占的體積比例，則可以建立水滴相的控制方程，包括連續(xù)方程和動量方程，分別為［8］：

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

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

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

與空氣控制方程類似，水滴相控制方程可以統(tǒng)一寫成輸運(yùn)方程的形式：

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

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

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

2.3 壁面邊界條件

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

3 水滴收集率的計算

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

其中，α∞為遠(yuǎn)場水滴容積分?jǐn)?shù)，u∞為遠(yuǎn)場來流速度，n為物面碰撞點(diǎn)處的單位法線向量。

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

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

4 大型風(fēng)力機(jī)水滴收集率計算

根據(jù)以上數(shù)值方法，開發(fā)了相應(yīng)的計算程序，并針對某1.5MW級的水平軸風(fēng)力機(jī)進(jìn)行了應(yīng)用計算。

4.1 計算外形及網(wǎng)格

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

圖1 計算外形及坐標(biāo)方向Fig.1 Wind turbine configuration and frame

圖2 輪轂附近表面網(wǎng)格 Fig.2 Surface grid near the hub

4.2 風(fēng)力機(jī)空氣流場計算結(jié)果

本文計算了兩種來流條件下風(fēng)力機(jī)的氣動特性，并與設(shè)計值進(jìn)行了對比，計算條件為：

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

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

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

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

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

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

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

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

4.3 風(fēng)力機(jī)水滴收集率計算結(jié)果

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

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

根據(jù)以上水滴收集率計算結(jié)果，可以對輪轂表面的結(jié)冰厚度進(jìn)行估算，物體表面某處結(jié)冰的厚度h可以表達(dá)為［8］：

圖5 風(fēng)力機(jī)表面水滴收集率分布云圖（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為凍結(jié)系數(shù)（0≤f≤1），β為當(dāng)?shù)氐乃问占剩琇WC為空氣中的液態(tài)水含量［15］，V為來流速度，dt為結(jié)冰時間，ρi為冰的密度。

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

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

圖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）

以上結(jié)果說明對于過冷水滴在風(fēng)力機(jī)葉片表面的結(jié)冰，其基本規(guī)律有兩個，一是冰主要結(jié)在葉片前緣，二是葉尖結(jié)冰遠(yuǎn)比葉根處嚴(yán)重。因此，在進(jìn)行結(jié)冰分析、設(shè)計結(jié)冰探測方法和防除冰方案時，葉片前緣，尤其是葉尖附近的葉片前緣，需要重點(diǎn)關(guān)注。

5 結(jié) 論

建立了適合于大型風(fēng)力機(jī)結(jié)冰過程水滴收集率計算的三維數(shù)值方法和計算程序，并對某1.5MW級水平軸風(fēng)力機(jī)的水滴收集率進(jìn)行了計算，得到如下結(jié)論：

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

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

（3）在進(jìn)行結(jié)冰分析、設(shè)計結(jié)冰探測方法和防除冰方案時，可以忽略輪轂和葉片根部附近結(jié)冰的影響，葉片前緣，尤其是葉尖附近的葉片前緣，需要重點(diǎn)關(guān)注。

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

［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，F(xiàn)ORTIN 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，F(xiàn)EDRIZZI 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］朱程香，王瓏，孫志國，等.風(fēng)力機(jī)葉片翼型的結(jié)冰數(shù)值模擬研究［J］.空氣動力學(xué)學(xué)報，2011，29（4）：522－528.

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

［9］易賢，桂業(yè)偉，朱國林.飛機(jī)三維結(jié)冰模型及其數(shù)值求解方法［J］.航空學(xué)報，2010，31（11）：2152－2158.

［10］易賢，桂業(yè)偉，朱國林，杜雁霞.運(yùn)輸機(jī)翼型結(jié)冰的計算和實(shí)驗［J］.航空動力學(xué)報，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］易賢.飛機(jī)積冰的數(shù)值計算與積冰試驗相似準(zhǔn)則研究［D］.綿陽：中國空氣動力研究與發(fā)展中心，2007.