黎芷毓，蔣興良，韓興波，任曉東，王洋洋

(1. 輸配電裝備及系統(tǒng)安全與新技術(shù)國(guó)家重點(diǎn)實(shí)驗(yàn)室(重慶大學(xué))， 重慶 400044；2. 重慶交通大學(xué) 機(jī)電與車輛工程學(xué)院， 重慶 400074)

風(fēng)能作為可再生能源，越來(lái)越受到人們的重視，近年來(lái)，中國(guó)風(fēng)電行業(yè)發(fā)展迅速，風(fēng)力發(fā)電機(jī)裝機(jī)容量逐年攀升，目前已處于世界第一[1-3]。然而，中國(guó)近年來(lái)極端氣象災(zāi)害頻發(fā)，在重慶、云南、湖南等華中和西南地區(qū)內(nèi)陸風(fēng)電場(chǎng)，冬季低溫且潮濕，風(fēng)力發(fā)電機(jī)易遭到覆冰災(zāi)害的影響[4-5]。結(jié)冰的風(fēng)力機(jī)葉片不僅會(huì)使得輸出功率降低，嚴(yán)重的會(huì)導(dǎo)致停機(jī)[6-7]。文獻(xiàn)[8]顯示1996年至2002年期間因?yàn)槿~片覆冰導(dǎo)致芬蘭某風(fēng)電場(chǎng)風(fēng)力發(fā)電機(jī)出現(xiàn)問(wèn)題的時(shí)長(zhǎng)達(dá)到8 000小時(shí)以上，約占低溫故障總時(shí)長(zhǎng)的70%。研究發(fā)現(xiàn)，風(fēng)力發(fā)電機(jī)結(jié)冰導(dǎo)致年發(fā)電量(AEP)降低17%，空氣動(dòng)力學(xué)特性降低20%～50%[9]。

由于對(duì)風(fēng)力發(fā)電機(jī)葉片進(jìn)行人工覆冰試驗(yàn)耗費(fèi)大量物力及財(cái)力，近年來(lái)，利用計(jì)算機(jī)技術(shù)對(duì)風(fēng)力發(fā)電機(jī)覆冰進(jìn)行預(yù)測(cè)研究成為了研究熱點(diǎn)。利用計(jì)算機(jī)信息技術(shù)對(duì)葉片覆冰進(jìn)行模擬計(jì)算，不僅可以為風(fēng)力發(fā)電機(jī)翼型的防冰設(shè)計(jì)提供技術(shù)參考，對(duì)后續(xù)防冰除冰工作及設(shè)備維護(hù)也有重要的指導(dǎo)意義。早期芬蘭技術(shù)研究中心Makkonen等[10-11]借鑒飛機(jī)覆冰數(shù)值模型開(kāi)發(fā)了一套風(fēng)力機(jī)覆冰代碼TURBICE，該模型可計(jì)算得到不同溫度下風(fēng)力機(jī)的覆冰形狀及質(zhì)量。文獻(xiàn)[12]對(duì)大型水平軸“NREL 5MW”風(fēng)力發(fā)電機(jī)葉片的覆冰增長(zhǎng)進(jìn)行了數(shù)值計(jì)算發(fā)現(xiàn)，隨著葉片幾何參數(shù)的變化，覆冰的增長(zhǎng)速率及冰形都有一定的改變。文獻(xiàn)[13]對(duì)大型風(fēng)力發(fā)電機(jī)葉片覆冰進(jìn)行預(yù)測(cè)，并分析了覆冰后轉(zhuǎn)矩及功率的損失，但其預(yù)測(cè)仍為二維模型，無(wú)法獲得葉片整體覆冰增長(zhǎng)趨勢(shì)。文獻(xiàn)[14]運(yùn)用商用軟件Fluent對(duì)水平軸風(fēng)力發(fā)電機(jī)三維雨凇覆冰進(jìn)行了數(shù)值計(jì)算及人工覆冰試驗(yàn)，然而該項(xiàng)研究?jī)H針對(duì)小型風(fēng)力發(fā)電機(jī)。文獻(xiàn)[15]對(duì)常用于大型風(fēng)力發(fā)電機(jī)的翼型覆冰進(jìn)行了試驗(yàn)及仿真計(jì)算，但計(jì)算仍然停留在單翼型，無(wú)法了解整體葉片覆冰趨勢(shì)。

風(fēng)力發(fā)電機(jī)葉片的覆冰預(yù)測(cè)模型不僅要考慮其準(zhǔn)確性，同時(shí)也應(yīng)考慮時(shí)間和開(kāi)銷，目前研究其氣流場(chǎng)的方法有：有限元法、有限體積法、有限差分法及邊界元法。本文采用的邊界元法計(jì)算風(fēng)力發(fā)電機(jī)葉片覆冰流場(chǎng)，通過(guò)降維的方法，可將三/二維問(wèn)題化為二/一維問(wèn)題進(jìn)行分析，具有計(jì)算量小，計(jì)算時(shí)間短的優(yōu)點(diǎn)[16]。由于風(fēng)力發(fā)電機(jī)葉片外部氣流場(chǎng)可以近似看作二維流動(dòng)，而大型風(fēng)力發(fā)電機(jī)葉片的尺寸較大，進(jìn)行三維葉片空氣流場(chǎng)計(jì)算不僅耗費(fèi)時(shí)間長(zhǎng)，且需要較大的計(jì)算機(jī)存儲(chǔ)空間，導(dǎo)致其三維覆冰模擬較難實(shí)現(xiàn)[17]。本文提出了一種切片-重構(gòu)的方法，對(duì)風(fēng)力發(fā)電機(jī)葉片進(jìn)行切片得到數(shù)個(gè)二維截面，基于邊界元法對(duì)葉片外部氣流場(chǎng)進(jìn)行計(jì)算，使用拉格朗日法跟蹤水滴軌跡，通過(guò)優(yōu)化選擇時(shí)間步長(zhǎng)，仿真模擬各二維風(fēng)力發(fā)電機(jī)葉片切片的覆冰增長(zhǎng)，最終通過(guò)截面覆冰重構(gòu)，得到三維葉片霧凇覆冰圖像。

1 風(fēng)力發(fā)電機(jī)葉片氣液兩相流數(shù)值模型

當(dāng)風(fēng)力發(fā)電機(jī)葉片在運(yùn)轉(zhuǎn)時(shí)，其外部氣流場(chǎng)合速度為風(fēng)速與該截面所擁有的線速度的相反值之和，稱之為相對(duì)風(fēng)速，也即來(lái)流合速度。同時(shí)，風(fēng)力發(fā)電機(jī)在運(yùn)行的過(guò)程中，其外部空間存在氣液(空氣-水滴)兩相流，隨著溫度的降低，空氣中的過(guò)冷卻水滴撞擊到風(fēng)力機(jī)葉片表面凍結(jié)形成覆冰。研究葉片外部空氣流場(chǎng)是研究過(guò)冷卻水滴運(yùn)動(dòng)軌跡的基礎(chǔ)，因此，首先要對(duì)外部流場(chǎng)的速度分布進(jìn)行計(jì)算，再將水滴視為離散相，根據(jù)空氣外流場(chǎng)速度分布進(jìn)一步計(jì)算每個(gè)水滴的運(yùn)動(dòng)軌跡，得到水滴在葉片表面的局部碰撞率。

邊界層理論表明，結(jié)構(gòu)物大氣繞流場(chǎng)可分為邊界層內(nèi)的層流或湍流區(qū)域，以及邊界層以外的非黏性勢(shì)流區(qū)域[18-19]。

1.1 外部空氣流場(chǎng)求解

如圖1所示，假設(shè)P為外部流場(chǎng)中的任一點(diǎn)(a、b、c、d、e、f)，Q為邊界上的點(diǎn)，在邊界層以外的勢(shì)流區(qū)，流體可以看作是理想流體，根據(jù)拉普拉斯方程：

(1)

式中：φ為速度勢(shì)；vx，vy分別為沿x軸方向，y軸方向的來(lái)流速度。

圖1 風(fēng)力發(fā)電機(jī)葉片截面邊界元?jiǎng)澐?/p>

邊界條件：風(fēng)力發(fā)電機(jī)葉片邊界為壁面，流體不可穿透，因此流體在邊界的法向速度為0，即

(2)

(3)

根據(jù)格林公式得到如下邊界積分方程，即流場(chǎng)勢(shì)函數(shù)求解控制方程：

(4)

(5)

式中：φ*(P,Q)為拉普拉斯算子基本解，r(P,Q)為P點(diǎn)與Q點(diǎn)之間的距離。

需要確定邊界上的勢(shì)函數(shù)值，對(duì)于邊界點(diǎn)P，為避免出現(xiàn)P與Q重合而產(chǎn)生奇異性，式(5)化為

(6)

式中C(P)是與P點(diǎn)處的邊界幾何形狀有關(guān)的常數(shù)，對(duì)于光滑邊界，C(P)=1/2。

本文采用的邊界元法，將風(fēng)力發(fā)電機(jī)葉片的邊界分割成n個(gè)邊界單元，整個(gè)邊界上的積分以n個(gè)邊界單元上的積分的和來(lái)表示，且將節(jié)點(diǎn)取為微元中點(diǎn)。通過(guò)常量單元離散，假定每個(gè)邊界單元上的勢(shì)函數(shù)值與其法向?qū)?shù)值在邊界單元上為常數(shù)，且等于節(jié)點(diǎn)的值。式(6)可離散為

(7)

其中：

(8)

(9)

(10)

式中：i為P點(diǎn)所在微元序號(hào)，j為Q點(diǎn)所在微元序號(hào)，δij為克羅內(nèi)克δ。

對(duì)于n個(gè)節(jié)點(diǎn)，得到聯(lián)立的一次方程組，用矩陣形式可表示為

H·φ=G·Q

(11)

式中：H和G是n×n階的系數(shù)矩陣，φ和Q分別是邊界單元節(jié)點(diǎn)的擾動(dòng)勢(shì)函數(shù)值和擾動(dòng)勢(shì)函數(shù)的法向?qū)?shù)值的列向量。

對(duì)于非邊界單元，式(6)可離散為

(12)

通過(guò)式(2)～(12)，采用高斯積分公式求解邊界積分方程，即可求出流域內(nèi)任意位置的氣流速度。

1.2 邊界層內(nèi)流場(chǎng)求解

臨界邊界層動(dòng)量厚度雷諾數(shù)可作為判斷流體從層流發(fā)展為湍流的準(zhǔn)則數(shù)[17]：

(13)

式中：θ為風(fēng)力發(fā)電機(jī)葉片表面邊界層動(dòng)量厚度；θtr為層流發(fā)展為湍流的分界點(diǎn)的邊界層動(dòng)量厚度；μ為空氣黏度系數(shù)；ρa(bǔ)為空氣的密度；Ue為邊界層速度，此處定義為99%來(lái)流速度。

卡門邊界層動(dòng)量積分關(guān)系式[20]為

(14)

式中：p為壓強(qiáng)，δ為邊界層厚度。

邊界層厚度δ可由邊界層動(dòng)量厚度θl估算[21]：

δ=8.5θl

(15)

邊界層動(dòng)量厚度θl[22]的表達(dá)式為

(16)

式中：str表示層流與湍流分界點(diǎn)位置，θi(str)為分界點(diǎn)處的邊界層動(dòng)量厚度。

摩擦切應(yīng)力τ0表達(dá)式[18-19]為

(17)

邊界層以內(nèi)的氣流速度[23-24]為

(18)

式中：h為距離風(fēng)力發(fā)電機(jī)葉片邊界的垂直距離，ds為邊界微元單位長(zhǎng)度。

求解式(13)～(18)，即可求出邊界層內(nèi)氣流速度。

1.3 過(guò)冷卻水滴軌跡求解

本文計(jì)算水滴軌跡時(shí)利用拉格朗日法，為了得到過(guò)冷卻水滴在風(fēng)力發(fā)電機(jī)葉片表面的局部碰撞率，需要對(duì)過(guò)冷卻水滴運(yùn)動(dòng)軌跡進(jìn)行計(jì)算。假設(shè)水滴在距離障礙物較遠(yuǎn)時(shí)(設(shè)水滴從距離葉片前緣10倍弦長(zhǎng)處發(fā)射)與空氣擁有相同的速度，靠近障礙物時(shí)，氣流沿障礙物表面繞流，過(guò)冷卻水滴軌跡由于慣性的作用，與氣流軌跡產(chǎn)生偏差，此時(shí)對(duì)過(guò)冷卻水滴進(jìn)行受力分析后可知，過(guò)冷卻水滴在運(yùn)動(dòng)過(guò)程中主要受到重力FG、氣流浮力Fb、水滴前后的壓差阻力Fg、表觀質(zhì)量力FM、空氣黏性阻力Fd的作用，其中占主要作用的是空氣黏性阻力Fd，此處僅考慮空氣黏性阻力的作用[25]，根據(jù)牛頓第二定律，得

(19)

相對(duì)雷諾數(shù)表示為

(20)

將式(20)代入式(19)，得

(21)

(22)

k為運(yùn)動(dòng)水滴慣性的斯托克斯數(shù)，計(jì)算公式為

(23)

使用差分法對(duì)式(19)求解，根據(jù)速度和加速度的一階和二階差分格式可得到水滴在tn時(shí)刻沿x軸方向的速度及加速度的表達(dá)式為

(24)

式中：Δt為時(shí)間步長(zhǎng)，xn+1、xn和xn-1分別為水滴在tn+1時(shí)刻、tn時(shí)刻和tn-1時(shí)刻的位置的橫坐標(biāo)。同理可得水滴沿y軸方向的速度及加速度。

水滴軌跡追蹤方程：

(25)

式中：va_x,n和va_y,n分別為水滴在tn時(shí)刻所在位置的氣流沿x、y軸方向的速度，vd_x,n和vd_y,n分別為水滴在tn時(shí)刻沿x、y軸方向的速度，yn+1、yn和yn-1分別為水滴在tn+1時(shí)刻、tn時(shí)刻和tn-1時(shí)刻的位置的縱坐標(biāo)。

通過(guò)式(25)可根據(jù)水滴初始位置和初始速度求出水滴軌跡。

1.4 過(guò)冷卻水滴碰撞率求解

風(fēng)力發(fā)電機(jī)葉片水滴局部碰撞率α1計(jì)算公式：

(26)

式中：dL為兩相鄰水滴在葉片上碰撞位置在葉片表面的距離，dY為相鄰水滴縱坐標(biāo)之差。

基于以上算法計(jì)算得出的葉片表面氣流場(chǎng)和水滴運(yùn)動(dòng)軌跡如圖2所示。

2 風(fēng)力發(fā)電機(jī)葉片霧凇覆冰模型計(jì)算

2.1 葉片表面水滴凍結(jié)特性

在過(guò)冷卻水滴碰撞到風(fēng)力機(jī)葉片表面的過(guò)程中，忽略水滴飛濺的影響，因此捕獲率α2≈1。覆冰是一種熱交換的過(guò)程，在葉片表面水滴凍結(jié)過(guò)程中，可以通過(guò)傳質(zhì)傳熱平衡方程計(jì)算葉片表面的凍結(jié)率，傳質(zhì)傳熱方程表示為

(27)

其中：mim為撞擊到葉片表面水滴的質(zhì)量；min和mout為流入和流出的水質(zhì)量；meva為蒸發(fā)水質(zhì)量，mice為凍結(jié)水質(zhì)量；Qf為空氣摩擦葉片表面產(chǎn)生的熱量；Qim為碰撞到葉片的過(guò)冷卻水滴帶來(lái)的動(dòng)能轉(zhuǎn)化的熱量；Qi為水滴凍結(jié)成冰后釋放的熱量；Qin和Qout為流入水帶來(lái)的能量和流出水帶走的能量，設(shè)風(fēng)力發(fā)電機(jī)葉片表面水膜溫度Tw=0 ℃，則水膜的流動(dòng)在微元之間不會(huì)產(chǎn)生能量交換，Qin、Qout均為0；Qeva為蒸發(fā)升華熱損失；Qc為空氣對(duì)流換熱損失的熱量；Qw為將葉片表面捕獲的過(guò)冷卻水滴加熱到水膜溫度所需要的熱量；Qk為冰面?zhèn)鲗?dǎo)熱損失；Qr為長(zhǎng)波輻射熱損失。

忽略式(27)中的較小的熱量項(xiàng)摩擦熱Qf和傳導(dǎo)熱損失Qk，將各熱量控制項(xiàng)表達(dá)式代入式(27)后可得到凍結(jié)率的計(jì)算方法：

(28)

2.2 葉片冰形計(jì)算

風(fēng)力發(fā)電機(jī)葉片表面微元控制體在覆冰速率dm/dt計(jì)算方法如下：

dm/dt=V·α1·α2·α3·w·ds

(29)

式中：V為水滴碰撞在葉片上的速度；w為空氣液態(tài)水含量，kg/m3；ds為控制體微元的長(zhǎng)度。

在覆冰不脫落的情況下，覆冰沿葉片法向生長(zhǎng)，每個(gè)控制體微元的覆冰厚度增長(zhǎng)速率為

(30)

式中：dt為覆冰時(shí)間步長(zhǎng)；ρice為冰密度，其經(jīng)驗(yàn)公式[26]為

(31)

其中R為Macklin冰密度參數(shù)，R=-URd/TS，TS為葉片表面溫度。

本文基于某型300 kW風(fēng)力發(fā)電機(jī)葉片(葉片全長(zhǎng)14.5 m，額定風(fēng)速13 m/s)進(jìn)行建模，將葉片等分為數(shù)個(gè)翼型截面進(jìn)行分析，分別計(jì)算獲得氣-液兩相流軌跡，并根據(jù)不同的截面選取合適的覆冰時(shí)間步長(zhǎng)dt，計(jì)算覆冰厚度d，進(jìn)而得出覆冰冰形，迭代直至達(dá)到總覆冰時(shí)長(zhǎng)Tice=30 min，計(jì)算流程見(jiàn)圖3。

圖3 三維翼型覆冰仿真流程圖

3 模型計(jì)算結(jié)果及分析

3.1 算法驗(yàn)證

水滴局部碰撞系數(shù)是覆冰計(jì)算中較為關(guān)鍵的參數(shù)，為了驗(yàn)證本文提出的計(jì)算空氣流場(chǎng)和水滴軌跡最終獲得碰撞率的方法的有效性，采用文獻(xiàn)[27]中的相關(guān)參數(shù)對(duì)二維NACA0012翼型進(jìn)行仿真。其中，來(lái)流合速度為44.39 m/s，液態(tài)水含量為0.78 g/m3，溫度為-7.65 ℃，水滴中值直徑為20 μm。并將試驗(yàn)結(jié)果及本文計(jì)算出的碰撞率計(jì)算結(jié)果進(jìn)行比較。

試驗(yàn)結(jié)果及計(jì)算結(jié)果如圖4所示，s表示碰撞點(diǎn)距離駐點(diǎn)的位置，翼型的水滴局部碰撞率呈現(xiàn)中間大兩頭小的現(xiàn)象。試驗(yàn)中s>0即葉片下表面碰撞范圍略大于上表面而碰撞率略小于上表面，這是因?yàn)樗卧谶\(yùn)動(dòng)的過(guò)程中受重力的作用會(huì)稍向下偏移，但計(jì)算時(shí)由于忽略了重力的影響導(dǎo)致碰撞范圍和碰撞率是上下對(duì)稱的。通過(guò)對(duì)比可知，本文使用邊界元法數(shù)值計(jì)算出的碰撞率與實(shí)驗(yàn)結(jié)果吻合良好。計(jì)算值的最大局部碰撞率略大于試驗(yàn)值，其他部分基本吻合。但計(jì)算的碰撞范圍小于試驗(yàn)值，產(chǎn)生誤差的可能原因有二：一是因?yàn)楸疚牟捎玫碾x散方法求解流場(chǎng)積分的方程和求解水滴軌跡所取的時(shí)間步長(zhǎng)導(dǎo)致的計(jì)算誤差；二是因?yàn)槲墨I(xiàn)中的覆冰試驗(yàn)不易控制，測(cè)量精度有一定的分散性，此外，在實(shí)際的覆冰試驗(yàn)中，水滴直徑的取值分布在一定范圍內(nèi)，其中大尺寸水滴會(huì)導(dǎo)致試驗(yàn)中水滴碰撞范圍增大，而本文采用單一水滴直徑進(jìn)行計(jì)算，水滴碰撞范圍偏小，但從整體上看，本文數(shù)值計(jì)算的水滴局部碰撞率與試驗(yàn)值非常接近，表明使用該方法計(jì)算風(fēng)力發(fā)電機(jī)葉片外部空氣流場(chǎng)、水滴運(yùn)動(dòng)軌跡和水滴碰撞是合理的。

圖4 NACA0012翼型水滴局部碰撞率對(duì)比

3.2 風(fēng)力發(fā)電機(jī)葉片的覆冰增長(zhǎng)特性

為了分析沿葉展方向的整體葉片覆冰增長(zhǎng)特性和風(fēng)速變化對(duì)覆冰形態(tài)及覆冰量的影響，本文選取了3個(gè)風(fēng)速進(jìn)行了冰形計(jì)算，其中，風(fēng)速U分別為3、6、9 m/s，對(duì)應(yīng)的風(fēng)輪轉(zhuǎn)速分別為27、32、38 r/min，液態(tài)水含量為0.3 g/m3，溫度為-15 ℃，水滴中值直徑為20 μm。圖5～7為計(jì)算得到的三維覆冰葉片。圖5顯示，覆冰區(qū)域基本主要集中在迎風(fēng)面，P—P截面(葉尖部)處背風(fēng)面尚存在一定的覆冰，而J—J截面(葉中部)和E—E截面(葉根部)處，冰基本只存在于葉片的迎風(fēng)面，在圖5(b)中幾乎看不到冰的存在，這是因?yàn)檠刂~展方向翼型的不斷變化和攻角的逐漸減小。同時(shí)，仿真結(jié)果顯示，覆冰厚度沿著葉展的方向逐漸增大，距離葉根4 m處(E—E截面)的覆冰冰厚僅為0.15 mm，而距離葉根14.44 m處(P—P截面)的葉尖部位，冰厚達(dá)到22 mm。過(guò)冷卻水滴在空氣中受到慣性和阻力的影響向前運(yùn)動(dòng)，如果氣流阻力占主要作用，則過(guò)冷卻水滴將沿著氣流的方向運(yùn)動(dòng)而繞過(guò)葉片，沿著葉展方向的翼型弦長(zhǎng)和厚度逐漸減小，使阻力的作用逐漸減小而慣性的影響逐漸增大，同時(shí)由于沿著葉展方相對(duì)速度的增加，更加大了慣性對(duì)水滴運(yùn)動(dòng)的影響，最終導(dǎo)致了越靠近于葉尖部位截面的覆冰越厚。

圖5 覆冰后三維葉片成像(U=9 m/s)

(a)葉尖冰形 (b)葉中冰形 (c)葉根冰形

圖7 不同風(fēng)速下葉尖成像

3.3 風(fēng)速對(duì)風(fēng)力發(fā)電機(jī)覆冰的影響分析

圖8顯示了在不同風(fēng)速下沿葉展方向不同部位的覆冰情況。圖8表明，風(fēng)速的增大會(huì)使得水滴局部碰撞率增加，且單位時(shí)間內(nèi)將有更多的水滴碰撞在截面上，導(dǎo)致覆冰厚度的增加，同時(shí)也會(huì)改變覆冰的面積和區(qū)域，在風(fēng)速為3 m/s時(shí)，葉尖部位的覆冰區(qū)域明顯整體小于風(fēng)速為6 m/s和9 m/s時(shí)，這是由于隨著風(fēng)速的增大，葉片轉(zhuǎn)速也將逐漸增大，氣流合速度隨之增大，過(guò)冷卻水滴的慣性對(duì)水滴運(yùn)動(dòng)的影響逐漸增大，在風(fēng)速較小時(shí)本會(huì)繞過(guò)葉片運(yùn)動(dòng)的水滴在風(fēng)速增大后可能會(huì)碰撞在葉片上，導(dǎo)致碰撞面積的增大。但隨著風(fēng)速的增加，背風(fēng)面覆冰區(qū)域逐漸減小，而迎風(fēng)面覆冰區(qū)域逐漸增加，這是由風(fēng)力發(fā)電機(jī)葉片的空氣動(dòng)力學(xué)特性決定的，隨著風(fēng)速和轉(zhuǎn)速的增加，氣流的攻角也會(huì)逐漸增加，使水滴碰撞區(qū)域從背風(fēng)面向迎風(fēng)面移動(dòng)。同時(shí)，在葉中部和葉根部也出現(xiàn)了這種趨勢(shì)，隨著覆冰面積增大的同時(shí)，覆冰區(qū)域從背風(fēng)面向迎風(fēng)面移動(dòng)，且隨著截面位置向輪轂靠近，覆冰區(qū)域?qū)⒅淮嬖谟谌~片的迎風(fēng)面，但與葉尖相比，葉根附近受風(fēng)速的影響明顯降低，在距離輪轂8 m(I—I截面)處，尚能看到明顯的厚冰，而更靠近葉根的兩個(gè)截面(F—F、C—C截面)幾乎已經(jīng)看不到冰。另外，風(fēng)速的變化對(duì)冰形幾乎沒(méi)有影響，同一截面下，3種風(fēng)速形成的覆冰冰形基本相同。為了分析風(fēng)速對(duì)冰厚及覆冰增長(zhǎng)速率的影響，設(shè)置沿截面法向的最大冰厚度為該處截面的冰厚。

(a) O—O截面(葉尖) (b)L—L截面(葉尖-葉中) (c)I—I截面(葉中)

(d)F—F截面(葉中-葉根) (e)C—C截面(葉根)

圖9顯示了沿葉展方向的不同風(fēng)速的冰厚變化，在3種風(fēng)速下，覆冰厚度沿著葉展方向逐漸增大，在葉尖處達(dá)到最大厚度，從葉根部到葉片中部(距離輪轂0～6 m)，3種風(fēng)速下的冰厚度相差不大，即冰厚受風(fēng)速的影響較小，且冰厚增長(zhǎng)較緩，此時(shí)的過(guò)冷卻水滴運(yùn)動(dòng)主要受阻力作用的影響，同時(shí)相對(duì)氣流合速度差別不大，且該部分并沒(méi)有明顯的覆冰現(xiàn)象，從葉片中部到葉尖部分(距離輪轂6～14.6 m)，水滴慣性逐漸占據(jù)主要作用，冰厚逐漸增大。隨著風(fēng)速的增大，冰厚沿葉展的增長(zhǎng)也隨之加快，冰厚差增大，一直到葉尖位置，當(dāng)風(fēng)速?gòu)? m/s增加到9 m/s時(shí)，冰厚從16 mm增長(zhǎng)到22 mm，冰厚增長(zhǎng)速率從0.53 mm/min增長(zhǎng)到0.73 mm/min。

圖9 不同風(fēng)速下冰厚隨葉展的變化

綜合圖8、9，受覆冰影響最大的區(qū)域是從葉片中部到葉片尖部，此趨勢(shì)與文獻(xiàn)[12]對(duì)兆瓦級(jí)大型風(fēng)機(jī)的覆冰仿真中計(jì)算出的冰厚增長(zhǎng)趨勢(shì)基本一致。

4 結(jié) 論

基于切片-重構(gòu)的思想，將風(fēng)力發(fā)電機(jī)葉片沿葉展切為數(shù)個(gè)截面，對(duì)于每個(gè)截面，采用邊界元法對(duì)風(fēng)力發(fā)電機(jī)葉片周圍大氣流場(chǎng)進(jìn)行計(jì)算，并用拉格朗日法對(duì)過(guò)冷卻水滴的軌跡進(jìn)行追蹤。仿真得到每個(gè)截面的冰形，再將所有截面及冰形進(jìn)行復(fù)原重構(gòu)，最終建立了風(fēng)力發(fā)電機(jī)葉片三維霧凇計(jì)算模型，同時(shí)對(duì)邊界元算法用于計(jì)算大氣流場(chǎng)的有效性進(jìn)行了檢驗(yàn)，分析了風(fēng)速對(duì)葉片霧凇覆冰的影響規(guī)律。通過(guò)對(duì)比和分析可得出如下結(jié)論：

1)在風(fēng)力發(fā)電機(jī)葉片霧凇覆冰模擬中，碰撞率的計(jì)算至關(guān)重要，本文方法所得計(jì)算值和試驗(yàn)值十分接近，證明所提出的計(jì)算方法是合理且有效的。

2)對(duì)某型300 kW風(fēng)電設(shè)備葉片霧凇覆冰進(jìn)行了仿真模擬，結(jié)果表明，霧凇覆冰冰形呈現(xiàn)流線型。覆冰主要聚積在葉片的前緣和迎風(fēng)面靠近前緣位置，沿著葉展的方向，覆冰逐漸增厚，覆冰區(qū)域占比增大。

3)風(fēng)速對(duì)霧凇覆冰的影響主要在冰厚和覆冰區(qū)域，對(duì)冰形的變化無(wú)明顯影響，且風(fēng)速的變化對(duì)覆冰的影響在葉尖部位更加顯著，在靠近葉根的位置，不同風(fēng)速產(chǎn)生的冰厚都較薄，從葉中部到葉根，3種風(fēng)速下冰厚的差別越來(lái)越大，直到葉尖冰厚度相差6 mm，且隨著風(fēng)速的增加，覆冰區(qū)域變大的同時(shí)，覆冰區(qū)域從背風(fēng)面向迎風(fēng)面移動(dòng)。

4)葉片的幾何參數(shù)(弦長(zhǎng)，厚度和扭轉(zhuǎn)角)的變化，都會(huì)對(duì)覆冰的增長(zhǎng)有一定的影響，這意味著通過(guò)優(yōu)化葉片的幾何參數(shù)，可以一定程度上降低其覆冰速率，從而減少寒冷潮濕地區(qū)安裝在風(fēng)電設(shè)備上的防冰和除冰裝置的數(shù)量。