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

(內蒙古科技大學 機械工程學院，內蒙古 包頭 014010)

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

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

1 研究方法

1.1 翼型簡介

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

1.2 數值計算模型

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

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

1.3 網格劃分模型以及邊界條件

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

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

(1)

(2)

(3)

式中，T為溫度，℃；V為來流風速，m/s.

1.4 結冰粗糙度

結冰粗糙度主要依據美國國家宇航局結冰經驗公式得到：

(4)

(5)

0.245 7LWC+1.257 1(LWC)2.

(6)

(7)

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

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

表1 工況

2 結果分析與討論

2.1 流場分析

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

圖1 流場壓強圖

2.2 水滴場分析

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

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

圖2 水收集系數分布圖

2.3 覆冰分析

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

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

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

圖3 覆冰坐標轉換圖

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

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

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

(10)

經過計算得出MAPE=0.457 7，因此證明以上擬合結果可行.

2.4 氣動特性分析

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

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

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

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

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

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

3 結論

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

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

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

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

y=0.058-0.107x-0.001x2

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