基于FLUENT的風(fēng)力機(jī)葉片翼型NREL S809覆冰的數(shù)值模擬

中南大學(xué)能源科學(xué)與工程學(xué)院 ■ 朱名巖 鄧勝祥

0 引言

隨著科學(xué)技術(shù)的發(fā)展，能源消耗與日俱增，傳統(tǒng)化石能源的不可再生性及其燃燒對環(huán)境的污染性，已經(jīng)成為世界各國共同關(guān)注的難題[1]。風(fēng)能作為可再生、無污染的新能源，分布廣泛、儲量豐富，已成為重要的化石能源替代品，世界各國越來越重視風(fēng)能的發(fā)展。

我國風(fēng)資源儲量高達(dá)2.5億kWh，風(fēng)資源雖然分布廣泛，但主要集中在我國的西北、東北、華北地區(qū)，以及東南沿海地區(qū)，由于“三北”等高緯度地區(qū)冬季溫度很容易達(dá)到0 ℃以下，所以風(fēng)電場中的風(fēng)力機(jī)常受到葉片結(jié)冰的困擾，防冰、除冰的研究工作迫在眉睫。易賢等[2]采用數(shù)值模擬計算的方法預(yù)測了翼型NACA0012前緣的霜冰積冰過程，并得出霜冰積冰的形狀。陳維建等[3]通過研究了粗糙度對翼型NACA0012的明冰積冰過程的影響，并得到翼型在攻角4°下的瘤狀冰形。鄧曉湖等[4]研究了霜冰的形成過程，以及覆冰前后流場的變化，指出覆冰導(dǎo)致速度沿翼型的分布紊亂，翼型的氣動性能急劇惡化。李聲茂等[5]對翼型NACA0015的結(jié)冰過程進(jìn)行了仿真模擬計算，并引入離散相模型DPM，最后指出在翼型的結(jié)冰過程中，過冷水滴量與風(fēng)速起著決定性的作用。朱程香等[6]的研究表明，葉片覆冰造成葉片的翼型幾何形狀被破壞，引起升力減小，阻力增大，升阻比降低，最嚴(yán)重時的降低幅度可達(dá)61%；覆冰導(dǎo)致翼型失速提前，嚴(yán)重破壞了翼型的氣動性能。Hochart等[7]通過數(shù)值模擬方法預(yù)測了水滴的運動軌跡、求解翼型覆冰形狀，研究了覆冰對翼型的氣動性能的影響，并對防冰、除冰系統(tǒng)的開發(fā)提出了建議。張聘亭等[8]利用ANSYS CFD軟件研究了風(fēng)力機(jī)專用翼型在4種覆冰形態(tài)下的動、靜態(tài)流場，并分析了4種覆冰形態(tài)對翼型靜態(tài)流場、動態(tài)流場的影響。

本文主要研究了覆冰對風(fēng)力機(jī)葉片翼型的氣動性能的影響，以經(jīng)典風(fēng)力機(jī)翼型NREL S809為例進(jìn)行仿真研究，建立了仿真模型，對計算域模型進(jìn)行網(wǎng)格劃分，通過網(wǎng)格無關(guān)性驗證，得到空氣流場的計算結(jié)果。利用FLUENT軟件自帶的離散相模型DPM及用戶自定義函數(shù)UDF，求解翼型表面的水滴運動軌跡方程、水滴的撞擊極限和局部收集系數(shù)；選用改進(jìn)的Messinger積冰傳熱模型計算積冰厚度，得到積冰形狀。分析在不同的攻角下，潔凈翼型與覆冰翼型的氣動特性，通過仿真計算得到與之對應(yīng)的氣動特性參數(shù)，并獲得相應(yīng)的速度云圖、流線圖、升阻力系數(shù)曲線圖等，進(jìn)而對結(jié)果進(jìn)行對比分析，得出覆冰對翼型的氣動性能的影響。

1 覆冰數(shù)值模擬方法

1.1 流場計算

本文選擇Spalart-Allmaras湍流模型，該模型是一種較為有效的低雷諾數(shù)下的湍流模型，適用于求解邊界層受粘性影響的區(qū)域，具有魯棒性好、計算速度快且準(zhǔn)確度高的優(yōu)點，適用于風(fēng)力機(jī)翼型流場的計算[9]。

本文以風(fēng)力機(jī)葉片翼型NREL S809為例，選用ICEM CFD軟件進(jìn)行建模，為了使翼型在流場中的計算不受到干擾且計算精確，計算域要比翼型大很多。如圖1所示，左側(cè)來流方向距離翼型前緣為15c(c為葉片翼弦長)，上下邊界也為15c，為了防止回流的干擾，尾緣到出口的邊界要長一點，為25c。本文為了準(zhǔn)確求解風(fēng)力機(jī)葉片翼型的仿真結(jié)果，近壁面邊界層網(wǎng)格劃分得很細(xì)。經(jīng)過網(wǎng)格的無關(guān)性驗證，最終確定網(wǎng)格的劃分，此時經(jīng)驗值y+的值為1。

圖1 全局和局部網(wǎng)格放大圖

邊界條件的設(shè)置[10]為：翼型表面和流域外邊界設(shè)為無滑移的壁面邊界；采用速度入口邊界，設(shè)置來流速度值并通過來流攻角確定來流的速度方向；設(shè)置壓力出口邊界條件，壓力值為1個大氣壓。

求解參數(shù)的設(shè)置為：來流按不可壓縮流處理，壓力-速度耦合雙精度選用SIMPLE算法求解。

1.2 水滴撞擊特性計算

1.2.1 水滴軌跡運動方程

為進(jìn)行過冷水滴的受力分析，進(jìn)行如下假設(shè)：1)來流中的水滴均勻分布，運動過程中其體積保持不變，以球形存在且直徑為deq；2)在整個過程中水滴的物性參數(shù)保持不變，且水滴不發(fā)生熱交換；3)以水滴初始速度值選擇來流速度值，水滴體積很小，不造成繞流的流場分布的變化；4)僅考慮作用在水滴上的黏性阻力、重力和浮力。

選用離散相模型DPM，根據(jù)拉格朗日法確定水滴運動方程，結(jié)合水滴受力的假設(shè)和分析，以及牛頓第二定律，確定水滴軌跡運動方程[11]為：

式中，ρa為空氣的密度；為重力加速度；Ad為水滴迎風(fēng)面積；Vd為水滴的體積；CD為阻力系數(shù)；為本地空氣流速；為水滴的速度；Md為水滴平均質(zhì)量；ρd為水滴的密度；為x方向的位移；t為水滴撞擊翼型所用時間。

采用近似求解法，經(jīng)過推導(dǎo)計算可得到方程的精確解[12]。

令：

求解方程的解，通過流場分布情況和定義水滴初始運動條件，求解水滴在每個時刻的位置，以及葉片和水滴的撞擊情況。

1.2.2 水滴撞擊極限

水滴撞擊極限[12]是指風(fēng)力機(jī)葉片在旋轉(zhuǎn)過程中，水滴所能撞擊到葉片翼型表面上下兩條相切軌跡包圍的表面的長度S與葉片翼型弦長c的比值，用Sm表示，其表達(dá)式為：

1.2.3 水滴收集系數(shù)

局部水收集系數(shù)β[12]是指微元表面的實際水收集量Wβ與該微元表面上最大能夠收集的水量Wβmax之比，因此它是表征微元表面的水收集能力的一個重要參數(shù)。β可表示為：

式中，ds為水滴上、下界的距離；dy為無窮遠(yuǎn)處的兩相鄰水滴的距離。

1.3 積冰熱力學(xué)模型

本文采用經(jīng)典Messinger積冰傳熱模型[13]的思想，依據(jù)質(zhì)量守恒定律，用進(jìn)入該控制體的質(zhì)量減去離開該控制體的質(zhì)量，即為該控制體內(nèi)冰的質(zhì)量；考慮覆冰表面上某一控制體積內(nèi)的能量傳遞情況，能量的傳遞項分成6項，假設(shè)不存在內(nèi)熱源，根據(jù)熱力學(xué)第一定律，控制體積表面的質(zhì)量守恒和能量平衡方程[12]可表示為：

式中，mim為水滴的質(zhì)量總和；mice為控制體內(nèi)水滴結(jié)成冰的總質(zhì)量；min、mva、mout分別為流入、蒸發(fā)損失、流出的液態(tài)水質(zhì)量；Qc為對流換熱；Qf為摩擦換熱；Qva為蒸發(fā)熱量；Qim為水滴自身熱能；QE為轉(zhuǎn)化熱量；Qso為相變熱量。

覆冰密度是計算覆冰形態(tài)的主要參數(shù)，根據(jù)每個控制體內(nèi)的覆冰質(zhì)量，可計算該控制體的覆冰厚度di，其可表示為[14]：

式中，Δt為覆冰時間步長；ρi為冰的密度；Ai為控制體長度。

在每個控制體覆冰厚度確定之后，沿表面法向增長，最終可確定冰形。

2 仿真結(jié)果及分析

2.1 覆冰計算結(jié)果

本文采用上節(jié)所述的方法進(jìn)行計算，得到覆冰計算結(jié)果，并進(jìn)行了適當(dāng)?shù)男拚?。?為覆冰計算工況參數(shù)表，圖2為覆冰計算結(jié)果圖。

表1 覆冰計算工況參數(shù)表

圖2 覆冰計算結(jié)果圖

2.2 升力、阻力系數(shù)和升阻比分析

圖3 是潔凈翼型與覆冰翼型的升力系數(shù)(圖3a)、阻力系數(shù)(圖3b)和升阻比(圖3c)曲線圖。由圖3可知，隨著攻角α的增大，覆冰翼型的升力系數(shù)和阻力系數(shù)與潔凈翼型的差距不斷增大。覆冰導(dǎo)致了最大升阻比下降，覆冰前翼型的最大升阻比為90.69，覆冰后下降為12.72，最大降幅為85.97%。由此可知，風(fēng)力機(jī)葉片覆冰會嚴(yán)重破壞翼型的升阻比，造成風(fēng)力機(jī)工作效率下降，因此，風(fēng)電場防冰、除冰工作勢在必行。

圖3 升、阻力系數(shù)和升阻比曲線圖

2.3 流線分析

圖4 分別為潔凈翼型(圖4a)和覆冰翼型(圖4b)在攻角為18e時的流線圖。從圖4可知，在攻角為18e時，潔凈翼型尾緣發(fā)生了分離，局部形成了較小的分離渦。覆冰翼型尾緣產(chǎn)生了更為明顯的分離，翼型的失速狀態(tài)更為嚴(yán)重，形成了更大的分離渦；同時在翼型的前緣開始發(fā)生失速現(xiàn)象，并產(chǎn)生了較小的分離渦。覆冰翼型整體的失速狀態(tài)更為嚴(yán)重，氣動性能產(chǎn)生了較為嚴(yán)重的惡化。

圖4 攻角為18e時翼型覆冰前后流線圖

2.4 壓力系數(shù)分析

圖5 是潔凈翼型與覆冰翼型分別在攻角0°、攻角5e和攻角14e時的壓力系數(shù)曲線圖。分析可知，覆冰翼型前緣幾何外形發(fā)生畸變，致使在翼型的前緣區(qū)域壓力系數(shù)Cp曲線出現(xiàn)震蕩現(xiàn)象；由于尾緣未受到覆冰的影響，在小攻角下，翼型尾緣的壓力系數(shù)Cp曲線與潔凈翼型吻合，但隨著攻角的增大，翼型尾緣壓力系數(shù)逐漸產(chǎn)生畸變；在攻角α=14e時，覆冰翼型在厚度與弦長的比值x/c=0.75時壓力系數(shù)開始出現(xiàn)偏差，這表明翼型尾緣邊界層開始發(fā)生分離，尾緣失速狀態(tài)提前，開始產(chǎn)生分離渦。

圖5 不同攻角下的壓力系數(shù)曲線圖

2.5 速度云圖分析

圖6 為潔凈翼型(見圖6a～圖6c)與覆冰翼型(見圖6d～圖6f)在攻角分別為0°、5°、14°時流場的速度分布云圖，可更直觀、更形象地了解覆冰對翼型周圍流場的影響。由圖6可知，對于潔凈翼型而言，隨著攻角的增大，駐點逐漸后移，且翼型尾緣部分出現(xiàn)的低速區(qū)逐漸擴(kuò)大；尤其在翼型進(jìn)入失速階段后，上表面的低速區(qū)擴(kuò)大，上表面壓力上升，翼型的升力下降；大攻角時會出現(xiàn)流動分離現(xiàn)象，翼型上表面流場紊亂，嚴(yán)重影響了翼型的升力。相較于潔凈翼型，覆冰翼型破壞了原翼型的幾何結(jié)構(gòu)，增加了翼型的弦長，導(dǎo)致流場被破壞，繼而導(dǎo)致翼型失速提前；翼型上表面由于覆冰的影響使低速區(qū)域增大且流場紊亂，造成翼型上表面壓力上升，嚴(yán)重破壞了翼型的氣動性能。

圖6 不同攻角下流場的速度分布云圖

3 結(jié)論

本文以風(fēng)力機(jī)葉片翼型NRELS809為研究對象，借助ANSYS FLUENT進(jìn)行外流場的數(shù)值模擬計算，并通過仿真計算得到覆冰冰形，對比覆冰前、后不同攻角下潔凈翼型與覆冰翼型的升力、阻力系數(shù)曲線圖，對翼型流線圖、速度云圖等進(jìn)行分析，得出覆冰對翼型的氣動性能的影響為：

1)在相同攻角下，覆冰造成翼型阻力系數(shù)增大而升力系數(shù)減小；隨著攻角的增大，覆冰翼型與潔凈翼型的升力系數(shù)和阻力系數(shù)的差距不斷增大。覆冰導(dǎo)致最大升阻比下降，覆冰前翼型的最大升阻比為90.69，覆冰后下降為12.72，最大降幅為85.97%。

2)覆冰翼型的失速狀態(tài)更為嚴(yán)重，尾緣產(chǎn)生更為明顯的分離渦，并且翼型前緣開始發(fā)生失速現(xiàn)象，產(chǎn)生了較小的分離渦。

3)相較于潔凈翼型，覆冰破壞了翼型的幾何外形，導(dǎo)致翼型失速狀態(tài)提前，翼型的氣動性能惡化，風(fēng)能利用率降低。覆冰還會導(dǎo)致風(fēng)力機(jī)荷載增大，增加其部件的不平衡性和疲勞程度，導(dǎo)致風(fēng)力機(jī)無法正常工作，嚴(yán)重影響了風(fēng)力機(jī)的發(fā)電效率。