S809風(fēng)力發(fā)電機(jī)葉片翼型霧凇覆冰程度及氣動性能模擬*

李震，萬濤，郝宇超，于躍

(內(nèi)蒙古科技大學(xué) 機(jī)械工程學(xué)院，內(nèi)蒙古 包頭 014010)

李翰濤等結(jié)果表明雨凇覆冰大幅降低二維翼型的氣動特性[1]；鄧杰等研究了覆冰以及除冰后葉片表面氣動特性，翼型前緣覆冰對氣動性能影響顯著，且相對厚度越小的翼型影響越大[2]；李巖等對NACA0018翼型進(jìn)行風(fēng)洞試驗，試驗分析結(jié)果表明攻角不同會導(dǎo)致冰形不同，來流方向的投影面積越大，結(jié)冰量與速率越大[3]；JAIWON S對NACA0012翼型實驗和冰形計算表明隨著冰面形狀轉(zhuǎn)變?yōu)橛员?，葉片的阻力系數(shù)在接近冰點時急劇增加，并且發(fā)現(xiàn)冰層表面的等效顆粒粗糙高度主要受液態(tài)水含量(LWC)、環(huán)境溫度(TS)、水滴中值直徑(MVD)等因素影響，受來流速度的影響較?。籎ANUSZ S顯示純霧氣結(jié)冰發(fā)生在比釉面結(jié)冰更低的溫度下，并導(dǎo)致沉積冰表面的粗糙度更高[4].

覆冰類型主要分為雨凇以及霧凇[5]，主要基于ANSYS FENSAP-ICE軟件研究S809翼型硬霧凇覆冰后的覆冰程度以及氣動特性，主要考慮環(huán)境溫度(TS)對覆冰厚度的影響和攻角對氣動特性的影響，并對覆冰模擬數(shù)據(jù)進(jìn)行一定程度擬合并且進(jìn)行數(shù)據(jù)驗證.

1 研究方法

1.1 翼型簡介

S809翼型是NERL S系列風(fēng)力發(fā)電機(jī)專用翼型之一，相比NACA系列翼型可使年能源利用率得到顯著提升，失速調(diào)節(jié)風(fēng)力機(jī)提高23%～35%，變槳距風(fēng)力機(jī)提高8%～20%，變轉(zhuǎn)速風(fēng)力機(jī)提高8%～10%[6].

1.2 數(shù)值計算模型

利用FENSAP-ICE軟件依次進(jìn)行空氣繞流、水收集系數(shù)以及后續(xù)的覆冰模擬計算，風(fēng)力機(jī)翼型流場控制方程采用時均Navier-Stokes方程.

選取Spalart-Allmaras湍流模型，S-A模型是1個一階模型，于二階模型而言，具有計算量不大、穩(wěn)定性良好以及計算精度較高的特點，在航空動力學(xué)計算領(lǐng)域中廣泛應(yīng)用[7]，并且在參考文獻(xiàn)[7]中，此模型的適用性已經(jīng)得到驗證.

1.3 網(wǎng)格劃分模型以及邊界條件

模擬地理環(huán)境參數(shù)是基于我國內(nèi)蒙某風(fēng)場進(jìn)行設(shè)置的，擬定風(fēng)力機(jī)功率為1.5 MW，根據(jù)風(fēng)力發(fā)電機(jī)設(shè)計手冊[8]計算公式，計算得葉片直徑為70 m，輪轂直徑設(shè)定為2 m，單葉片長度為34 m，海拔高度為2 170 m，風(fēng)速設(shè)定為15 m/s.根據(jù)參考文獻(xiàn)[9]，[10]，風(fēng)力發(fā)電故障高發(fā)時間段為每年的12 m、次年1 m，于中國氣象網(wǎng)歷年數(shù)據(jù)處，查出此處每年12 m、次年1 m風(fēng)場得相對濕度以及平均氣溫，根據(jù)飽和空氣濕度對照表查得出空氣液態(tài)水含量為0.8 g/m3.根據(jù)參考文獻(xiàn)[6]，[11]，[12]，雨凇(glaze)、硬霧凇(hard rime)、軟霧凇(soft rime)形成條件可用式(1)～(3)近似表達(dá)，取液滴中值直徑15 mm，計算域采用(c為弦長)，通過計算得出在攻角為6°時，弦長為2.60 m[13,8]，將左側(cè)圓弧以及上下邊界設(shè)置為速度入口，右側(cè)邊界設(shè)置為壓力出口，壁面設(shè)置為無滑移壁面.

網(wǎng)格劃分四邊形共計60 643個，節(jié)點數(shù)共計61 137個，經(jīng)檢測網(wǎng)格質(zhì)量符合模擬需求，滿足網(wǎng)格無關(guān)性要求，計算工況如表1所示.

(1)

(2)

(3)

式中，T為溫度，℃；V為來流風(fēng)速，m/s.

1.4 結(jié)冰粗糙度

結(jié)冰粗糙度主要依據(jù)美國國家宇航局結(jié)冰經(jīng)驗公式得到：

(4)

(5)

0.245 7LWC+1.257 1(LWC)2.

(6)

(7)

式中：基準(zhǔn)值(ks/c)base=0.001 177，ks：表面粗糙度，mm；MVD：液態(tài)水含量，g/m3；TS：環(huán)境溫度，℃；c：翼型弦長，m.

粗糙度的設(shè)定主要是為了計算表面結(jié)冰的流動解，同樣會增加熱通量以增加結(jié)冰速率和厚度.設(shè)置葉片表面粗糙度為0.5 mm，此數(shù)值適用范圍廣泛，采用此默認(rèn)數(shù)值并且設(shè)為初始表面粗糙度.利用式(4)經(jīng)過計算在環(huán)境溫度為-15 ℃工況為表1時，表面粗糙度為0.68 mm，此數(shù)值設(shè)定為-15 ℃霧凇結(jié)冰狀態(tài)下的表面粗糙度，以正交對比進(jìn)行后續(xù)模擬.

表1 工況

2 結(jié)果分析與討論

2.1 流場分析

首先對S809翼型進(jìn)行流場計算，計算溫度從-2 ℃開始由高到低，每2攝氏度計算一次，直至計算到-20 ℃，流場計算結(jié)果、壓強(qiáng)數(shù)據(jù)圖如圖1，從翼型表面壓力分布可以看出隨著溫度的降低，上翼型表面壓力隨溫度一同出現(xiàn)梯度較少狀況，下翼型表面的壓力同樣呈現(xiàn)此變化.由于來流空氣視為不可壓縮氣體，根據(jù)查理定律：一定質(zhì)量的空氣在體積不變的情況下，壓強(qiáng)與熱力學(xué)溫度成正比，圖示結(jié)果與此結(jié)論相符，由此判定模擬結(jié)果是無誤的.

圖1 流場壓強(qiáng)圖

2.2 水滴場分析

在空氣流場的基礎(chǔ)上進(jìn)行水滴場的計算，計算結(jié)果如圖2，為水收集系數(shù)數(shù)據(jù)分布圖.由水收集系數(shù)數(shù)據(jù)全局圖可以看出翼型下表面水滴主要撞擊在翼型前緣與弦長的2/5之間，上表面水滴主要撞擊在翼型前緣與弦長的3/20之間.

由圖2可以看出，上翼型下表面水滴收集量從-2 ℃開始至-20 ℃隨著溫度的降低呈現(xiàn)階梯狀增長.另外可以明顯地看出水滴撞擊主要集中在翼型的下表面，迎風(fēng)區(qū)撞擊范圍遠(yuǎn)大于背風(fēng)區(qū)[14]，上翼型上表面水收集系數(shù)在溫度-8 ℃～-2 ℃和-20 ℃～-10 ℃區(qū)間內(nèi)，均表現(xiàn)均勻上升，但是后者水滴收集范圍更廣，水滴收集量更多.

圖2 水收集系數(shù)分布圖

2.3 覆冰分析

最后在水滴場基礎(chǔ)上進(jìn)行覆冰模擬，由圖3可知，模擬覆冰結(jié)果中冰形主要集中在翼型的前緣以及下表面.隨著溫度的降低，覆冰厚度呈現(xiàn)階梯式均勻增長，造成此現(xiàn)象的原因是在結(jié)冰交界處，冰點溫度保持零攝氏度不變，而外界溫度的降低使得交界處的水分無法放熱，因此只能通過結(jié)冰來達(dá)到溫度一致，由此導(dǎo)致交界處的水滴結(jié)冰使冰的厚度增加.

李聲茂定義了無量綱常數(shù)冰厚比：最大覆冰厚度與弦長的比值以此來表征覆冰程度[15].由于覆冰模擬是在笛卡爾坐標(biāo)系下進(jìn)行的，而覆冰的主要區(qū)域也是翼型的前緣以及迎風(fēng)區(qū)，僅僅使用最大覆冰厚度與現(xiàn)場的比值雖然可以一定程度上表達(dá)覆冰程度，但是不夠準(zhǔn)確，翼型的覆冰還要考慮冰的密度以及體積對機(jī)翼造成的相關(guān)影響.

文章采用坐標(biāo)轉(zhuǎn)換法重新表達(dá)覆冰程度以更準(zhǔn)確表征覆冰程度并進(jìn)行量化表征.將在笛卡爾坐標(biāo)系中的覆冰翼型轉(zhuǎn)換到極坐標(biāo)系中，由此再將覆冰翼型角度坐標(biāo)與初始翼型角度坐標(biāo)相減得到覆冰區(qū)域的Δθ坐標(biāo)，計算覆冰的角度范圍.接著再將覆冰翼型長度與初始翼型長度相減，得到各個角度上的覆冰增量Δρ，將某一溫度下的覆冰增長量累加再除以此溫度下的覆冰角度范圍得到此溫度下的覆冰程度(覆冰角度范圍可從極坐標(biāo)軸讀出).按照此方法進(jìn)行分析得到如圖3結(jié)果，并經(jīng)過計算得到表2數(shù)據(jù).

圖3 覆冰坐標(biāo)轉(zhuǎn)換圖

表2 不同溫度覆冰總厚度及均值

在表1工況條件下，計算得到的覆冰增長量數(shù)據(jù)表，從數(shù)據(jù)表不難看出覆冰厚度-2 ℃～-20 ℃基本呈現(xiàn)規(guī)律增長，每2個溫度間隔冰厚增長量呈現(xiàn)由快之慢的現(xiàn)象.對數(shù)據(jù)表2用origin軟件進(jìn)行公式擬合，origin對多項式擬合是否準(zhǔn)確判定，依據(jù)是相關(guān)系數(shù)R平方值(回歸方程整體擬合度)是否接近1，越接近1則多項式擬合越準(zhǔn)確，一般認(rèn)定R平方大于0.99即為準(zhǔn)確，此處計算R2>0.999.考慮到隨著外界溫度的降低，結(jié)冰類型改變[16]，冰的表面粗糙度不斷增大，從明冰轉(zhuǎn)變?yōu)樗鵞17]，最終導(dǎo)致結(jié)冰平均厚度驟增，對比擬合的一次、二次、三次多項式，此處選擇二次多項式(8)并進(jìn)行簡化為最終結(jié)果式(9)：

y=(0.058±0.010 82)+(-0.107 55±0.002 26)x+(-0.001 14±1.001 16E-4)x2

y=0.058-0.107 55x-0.001 14x2.

(9)

式中：x為自變量溫度，℃；y為因變量冰厚平均值，mm.

為了驗證表達(dá)式的準(zhǔn)確與否，再次計算溫度-9，-13，-19和-22 ℃條件下的覆冰厚度，經(jīng)過模擬以及計算，得出如表4結(jié)果，由平均絕對百分比誤差MAPE(Mean Absolute Percentage Error)公式(10)，對以上數(shù)據(jù)進(jìn)行預(yù)測度誤差分析，計算數(shù)值越接近0模型越準(zhǔn)確：

(10)

經(jīng)過計算得出MAPE=0.457 7，因此證明以上擬合結(jié)果可行.

2.4 氣動特性分析

通過模擬計算得到覆冰外形，計算在初始翼型和-15 ℃覆冰外形進(jìn)行靜態(tài)以及此兩種狀態(tài)下葉尖速比為3的氣動特性，攻角范圍為-10°～20°，計算結(jié)果如圖4～7所示.由圖4可知，在初始翼型-10°～4°攻角范圍內(nèi)，葉尖速比3狀態(tài)下翼型氣動特性優(yōu)于靜止?fàn)顟B(tài)下氣動特性，在其余區(qū)間內(nèi)則是靜止?fàn)顟B(tài)氣動特性優(yōu)于轉(zhuǎn)動狀態(tài)；由圖5可知，-15 ℃狀態(tài)下結(jié)冰翼型3°～20°攻角范圍內(nèi)靜態(tài)氣動特性優(yōu)于旋轉(zhuǎn)狀態(tài)，并且在超過20°攻角時靜止?fàn)顟B(tài)氣動特性降幅超過旋轉(zhuǎn)狀態(tài)；由圖6可知，靜止?fàn)顟B(tài)下結(jié)冰翼型在1°～20°范圍內(nèi)氣動性能優(yōu)于初始翼型，升阻比趨勢表明初始翼型氣動特性在超過20°攻角時將優(yōu)于結(jié)冰翼型；由圖7可知，-10°～5°攻角區(qū)間內(nèi)霧凇狀態(tài)下氣動特性優(yōu)于初始狀態(tài)，另外5°～20°攻角范圍內(nèi)2種狀態(tài)下氣動特性基本一致，由升力系數(shù)以及阻力系數(shù)漲幅趨勢可得到：在大于20°的攻角范圍結(jié)冰初始翼型氣動特性將優(yōu)于霧凇狀態(tài).

圖4 初始翼型氣動性能曲線(a)阻力系數(shù)；(b)升力系數(shù)；(c)升阻比

圖5 結(jié)冰翼型氣動性能曲線(a)升力系數(shù)；(b)阻力系數(shù)；(c)升阻比

圖6 靜態(tài)翼型氣動性能對比曲線(a)升力系數(shù)；(b)阻力系數(shù)；(c)升阻比

圖7 動態(tài)翼型氣動性能對比曲線(a)升力系數(shù)；(b)阻力系數(shù)；(c)升阻比

(氣動特性圖中字母S表示恒定來流風(fēng)速的靜態(tài)，M表示葉尖速比為3的轉(zhuǎn)動狀態(tài))

3 結(jié)論

文章利用ANSYS FENSAP-ICE軟件對S809翼型進(jìn)行不同溫度下覆冰狀態(tài)以及不同攻角下潔凈與覆冰葉片進(jìn)行模擬，得到以下結(jié)論：

(1)對于潔凈S809翼型的空氣流場，翼型表面的壓強(qiáng)隨著溫度的降低而呈現(xiàn)均勻降低，符合查理定律.

(2)對于翼型的覆冰程度，可采用將笛卡爾坐標(biāo)系轉(zhuǎn)換為極坐標(biāo)系的方法，進(jìn)行相同θ度上ρ相減得到Δρ，接著求出重復(fù)覆冰的θ范圍角，累加Δρ并除以累加角得到覆冰均值，接著求出覆冰增長公式，對結(jié)冰相關(guān)計算有較大幫助.

(3)對潔凈S809翼型，在0 ℃～-20 ℃溫度范圍內(nèi)，覆冰是呈現(xiàn)規(guī)律增長的，并且滿足覆冰增長公式：

y=0.058-0.107x-0.001x2

(4)對比翼型初始狀態(tài)與霧凇狀態(tài)工作條件的氣動性能，在-10°～4°以及大于20°攻角范圍內(nèi)初始翼型氣動特性優(yōu)于霧凇狀態(tài).