不同覆冰形態下風力機葉片翼型氣動特性分析

徐琛苑 湯斯琦

DOI：10.20145/j.32.1894.20240205

作者簡介：徐琛苑（1996—），女，碩士研究生；研究方向：風力機葉片氣動特性，新能源電力系統規劃。xuchenyuan_grid@163.com

摘要：低溫環境下風力機葉片常面臨覆冰的危險，研究覆冰對翼型氣動特性的影響對覆冰葉片氣動特性分析具有重要意義。文章采用基于有限體積法的數值模擬算法對S809二維翼型氣動特性及覆冰影響進行了模擬分析。通過對比分析明冰與霜冰兩種覆冰形態對翼型氣動特性的影響程度，發現霜冰對翼型氣動特性影響不大，而明冰則會嚴重惡化翼型的氣動特性，甚至可能導致負阻力。

關鍵詞：風力機葉片；氣動特性；覆冰影響；數值模擬

中圖分類號： TK89? 文獻標志碼：? A

0? 引言

隨著全球變暖，極端氣候發生越來越頻繁，風力機冰凍災害問題也越來越突出。葉片表面覆冰是風力機冰凍災害最為突出的問題之一。葉片表面覆冰不僅給葉片增加了一個附加重力，還會改變葉片的幾何外形，從而影響其氣動特性［1-3］。對于MW級大型風力機超長柔性葉片，該問題尤為突出。葉片覆冰機理十分復雜，不僅與溫度、濕度、風速及風向變化等周圍環境有關，還與風力機葉片本身氣動特性及其工作狀態有關，目前尚無原型風力機三維葉片覆冰形成機理的理論解釋，二維翼型覆冰形成機制也在研究中。

根據葉片表面覆冰的幾何形狀，可以將葉片覆冰形態分為兩類：明冰［4-5］和霜冰［6-7］。當葉片周圍溫度低且濕度不是很高時，大氣中的小尺寸液滴在與葉片表面碰撞的過程中，瞬間凍結，形成不透明霜狀，此時的覆冰形態稱為霜冰［6-7］。霜冰形成過程中，液滴與葉片碰撞凍結時間極短，液滴碰撞凍結位置主要受葉片周圍流場及葉片表面邊界層流動影響，使得覆冰形態與葉片幾何外形類似，對葉片氣動特性影響并不是很大。當環境溫度相對較高但低于冰點，且濕度很高時，大氣中液滴的尺寸較大，大尺寸液滴與結構碰撞后凍結形成的透明冰型稱為明冰［4-5］。因為溫度不是很低且液滴尺寸較大，液滴與結構發生碰撞后并不能馬上完全凍結，只有一部分在碰撞處發生凍結，而另外一部分則沿著結構表面發生流動逐漸凍結，其在葉片表面的流動不僅與葉片表面邊界層流動有關，還與葉片運動狀態，如轉速、槳距角、方位角等有關，這導致明冰的幾何外形非常復雜且很難預測，通常會在液滴撞擊點附近形成一些尖角。明冰通常具有非流線型氣動外形，它將嚴重惡化結構的氣動特性，導致大幅度的流動分離從而使阻力大幅增加，升力大幅減小。

為了準確模擬覆冰對葉片氣動特性的影響，需要準確獲得葉片表面的壓力分布。一種方法是采用大渦模擬方法或者直接數值模擬方法，對葉片表面邊界層流動進行準確解析，但解析邊界層流動需要大量網格，計算成本過高。另一種方法則是采用湍流模型結合壁面函數的方法對邊界層流動進行模化。不同的湍流模型對應于不同的壁面函數，例如k-ε 湍流模型通常和標準壁面函數結合使用，增強型壁面函數則通常和k-ω 湍流模型結合使用。采用何種湍流模型與壁面函數則需要一定的經驗。例如，文獻［8］通過S809二維翼型的數值模擬發現，精確的預測轉捩分離點的位置是準確模擬翼型壓力分布曲線的關鍵。文獻［9］認為，湍流邊界層對過渡區非常敏感，針對渦粘模型對邊界層流動分離預測通常有所提前的弊端，提出了一種人工減小過渡區渦粘性的方法來延緩分離的發生，但該方法需要人為設定部分參數，且不同工況下參數并不相同，參數設定需要一定的經驗。

風力機葉片一般都是通過一系列翼型扭轉堆疊而成，研究二維翼型的氣動特性對三維葉片的氣動研究具有重要意義。此外，葉素理論［10］、制動線［11］、制動面［12］等葉片氣動特性計算方法都是基于二維翼型氣動特性附加一些修正所建立的，因此，研究二維翼型氣動特性覆冰影響分析具有重要意義。本文選取S809翼型為研究對象，首先通過對比不同的湍流模型與邊界層處理方法的模擬效果，驗證數值模擬算法并確定最優湍流模型與邊界層處理方法，進一步分析翼型表面不同覆冰形態對翼型氣動特性的影響規律。

1? 數值模擬參數

新能源科技

新能源科技

1.1? 幾何模型

本研究選用經典翼型——S809翼型作為研究對象，其幾何外形如圖1所示。模擬時取弦長為

C=600mm，計算域設置為矩形計算域，翼型前緣距離入口設為

10C，距離出口為20C，翼型中心線距離計算域上下邊界都為5C。邊界條件設置為速度入口、壓力出口，上下邊界為周期性邊界條件。來流選為均勻層流來流，來流速度設為U=51.7m/s，對應的雷諾數為Re=UC/ν=2×106，其中運動粘度

ν=1.55×10-6 m2/s。分別計算翼型在來流風攻角（Angle of attack，AOA）為0°、1.02°、5.13°、9.22°、14.24°和20.15°時的氣動特性。

1.2? 網格劃分

采用C形拓撲結構的結構化網格劃分方案進行網格劃分。為了模化邊界層流動，邊界層網格劃分需要滿足一定要求，且對于不同的壁面模型，邊界層網格劃分要求并不相同，即第一層網格的高度與邊界層內網格節點的數量需要滿足一定要求。邊界層網格劃分時需要保證邊界層內有足夠數量的網格節點以模擬邊界層的發展，一般需要在邊界層內至少布置15個網格節點。根據壁面定律，第一層網格的位置可以通過其與邊界層流動的關系來估算：

y1=y+νuτ（1）

式（1）中y1為第一層網格高度。uτ為摩擦速度，其與壁面剪切應力有關。在計算前并沒有壁面剪切應力信息，因此，需要采用經驗公式對其進行估算，經驗公式為：

uτ=0.5CfU2（2）

式（2）中Cf為摩擦系數，根據F.M.White算法，Cf=0.026/Re1/7。

對于不同的湍流模型及壁面模型，邊界層網格第一層網格高度不同，即y+取值不同。采用高雷諾數湍流模型時，結合標準壁面函數，要求第一層網格節點布置在湍流充分發展區，即要求y+≥30；采用低雷諾數湍流模型時，結合增強壁面函數，要求第一層網格節點布置在層流底層，即要求

y+≤5。針對本文模擬工況，當取y+=30時，由式（2）計算可得第一層網格高度

y1=0.2mm；當取

y+=5，由式（2）計算可得第一層網格高度

y1=0.03mm。需要注意的是，這里計算所得的第一層網格高度是根據經驗公式估算所得，可能并不滿足假設條件，一般需要根據模擬所得

y+的真實值進行調整，可以說，第一層網格高度的確定是一個重復調整的過程。來流風攻角為5.13°時的整體網格劃分結果及局部網格示意如圖2所示。

2? 數值驗證

通過對非覆冰狀態下的S809翼型進行模擬分析，驗證數值模擬方案的可靠性，并確定最優的湍流模型及壁面函數組合。

2.1? 湍流模型對比分析

根據前人實驗結果可知，S809翼型在攻角為14.24°和20.15°時處于失速區，邊界層流動涉及轉捩及流動分離，此時邊界層流動狀態的數值模擬對湍流模型、壁面函數的選擇及網格的劃分最為敏感。因此，通過模擬翼型在來流風攻角為14.24°和20.15°時的氣動特性與流場特征，可以最為可靠的對比不同湍流模型與壁面模型組合及網格劃分對模擬結果的影響。本研究共對比了Realize k-ε 湍流模型（以下簡稱“RKE湍流模型”）、SST

k-ω湍流模型和Spalart-Allmaras湍流模型（以下簡稱“S-A湍流模型”）3種湍流模型。其中RKE湍流模型為高雷諾數模型，可采用標準壁面函數，y+取值為

y+=30；SST k-ω 湍流模型為低雷諾數模型，

y+可取值為

y+=5；S-A湍流模型可采用增強壁面函數，

y+取值為

y+=5。

通過模擬所得y1值對網格進行驗證。經過反復模擬與調整發現，當第一層網格高度取值為

y1=0.6mm時，對應的第一層網格無量綱高度

y+=30；而當第一層網格的無量綱高度為

y+=5時，第一層網格厚度取值應為

y1=0.09mm。

y1=0.09mm時翼型上下表面不同位置處模擬計算所得

y+的結果如圖3所示，可以看出此時的

y1取值滿足要求。

通常采用無量綱化的翼型表面壓力系數分布曲線來描述翼型的氣動特性。翼型表面壓力系數計算公式如下：

Cp=P-P∞0.5ρU2（3）

式（3）中P為翼型表面壓力，P∞為環境壓力，U為來流風速，ρ為大氣密度。

.

不同湍流模型模擬所得翼型壓力系數曲線與實驗結果的對比圖如圖4所示。從圖4中可以看出，在來流風攻角為14.24°時，3種湍流模型模擬所得的壓力系數分布曲線與實驗所得壓力系數分布曲線都符合的很好；而在來流風攻角為20.15°時，RKE湍流模型模擬結果與實驗結果符合最好，SST k-ω 湍流模型模擬所得壓力曲線在翼型前緣上側有所偏低，即對邊界層分離預測有所提前，而S-A湍流模型模型模擬所得壓力曲線在翼型前緣上側有所偏大，在翼型前緣下側則有所偏低，即對邊界層分離的預測有所延遲。綜上所述，RKE湍流模型在3種湍流模型中對分離點的預測效果最好。不同湍流模型對邊界層分離的預測也可以從速度流場云圖中看出，如圖5所示。

注：左側圖攻角為14.24°，右側圖攻角為20.15°

綜合上述分析可得：在較小攻角情況下，3種湍流模型對翼型氣動特性模擬結果都很準確；但在攻角較大時，RKE湍流模型相對于SST k-ω 湍流模型和S-A湍流模型對翼型的氣動特性模擬相對更為

準確。

2.2? 數值模擬結果驗證

通過對比實驗與模擬所得的不同攻角下S809翼型的升阻力系數，對本文數值模擬方案進行驗證。基于上述對比分析，其他攻角下翼型氣動特性的模擬都采用RKE湍流模型進行模擬。升力系數及阻力系數計算公式如下：

CL=L0.5ρU2C（4）

CD=D0.5ρU2C（5）

式中L，D分別為翼型升力和阻力。

模擬結果與實驗結果符合的很好，驗證了本文所用數值方案的可靠性，如圖6所示。此外，從圖6中還可以看出，在來流風攻角為9.22°～20.15°時，由于翼型上側邊界層流動分離的發生，隨著來流攻角的增大，翼型升力系數增長變緩，在較大攻角時，翼型升力系數減小。對于阻力系數，隨著攻角的增大，阻力系數大幅增加。這與翼型進入失速區升阻力隨來流攻角的變化規律一致。

3? 覆冰影響分析

3.1? 覆冰狀態

本節對比分析明冰和霜冰兩種覆冰形態對翼型氣動特性的影響。基于大量葉片覆冰現場觀測數據對葉片覆冰形態進行測繪并進行一定簡化，確定翼型覆冰形態，其中霜冰與葉片幾何外形相似，其厚度分布為類拋物線分布，如圖7所示，在迎風點處覆冰厚度最大，本研究中霜冰最大厚度選取為0.05C，覆冰區域長度大致為0～0.3C；明冰覆冰形態較為復雜，如圖8所示，其尖角突出長度為0.1C，覆冰區域大致為0～0.2C。

3.2? 覆冰氣動影響分析

覆冰狀態下S809二維翼型氣動特性的模擬方案與非覆冰狀態下S809二維翼型氣動特性的模擬方案一致，即相同的計算域、相似的網格劃、同樣采用RKE湍流及標準壁面函數，具體模擬參數設置此處不再贅述。

霜冰覆冰形態下，模擬所得覆冰翼型壓力系數分布曲線與非覆冰狀態下實驗所得壓力系數分布曲線對比，如圖9所示。從圖9可以看出，覆冰狀態下翼型上下側壓差略有減小。在來流攻角較小時（0°、1.02°和5.13°），覆冰翼型非覆冰區域（0.3C～C）的壓力系數分布與非覆冰翼型相同。而在來流攻角為9.22°及14.24°時，由于覆冰的影響，在翼型上側0.3C～0.5C區域，壓力系數有所降低，在翼型上側0.5C～C區域壓力系數與非覆冰翼形一致，翼型下側非覆冰區域則不受覆冰的影響。值得一提的是，在這兩個攻角下，邊界層流動分離點發生的位置也沒有改變。在來流風攻角為20.15°時，非覆冰區域的壓力分布系數也不受覆冰的影響，但此時流動分離點由于覆冰的存在有所提前。

不同攻角下霜冰覆冰翼型繞流場速度云圖如圖10所示，可以看出，覆冰形態為霜冰時，速度流場與非覆冰翼型繞流場類似，這與霜冰覆冰翼型壓力系數分布曲線結果一致。此外，霜冰覆冰翼型表面邊界層流動分離點的位置也可以從速度流場云圖看出。從速度流場云圖所得分離點位置與通過壓力系數分布曲線所得分離點位置一致，且從速度流場云圖可以看出，在來流風攻角小于9.22°時，沒有流動分離發生，

在來流風攻角為9.22°、14.24°及20.15°時，翼型上側發生了穩定的流動分離，但沒有漩渦脫落產生，與非覆冰翼型類似。

明冰覆冰形態下模擬所得翼型壓力系數分布曲線與實驗所得非覆冰翼型壓力系數分布曲線對比，如圖11所示。從圖11中可以看出，由于明冰復雜幾何外形的影響，翼型壓力系數發生了劇烈變化，翼型上下側壓力差急劇減小，即明冰使得翼型的氣動特性發生了嚴重的惡化。在來流風攻角

≥5.13°時，翼型上側流動發生了完全分離。值得一提的是，在來流風攻角為20.15°時，由于非覆冰翼型分離點位于距離前緣0.2C處，非常接近迎風點，所以在翼型非覆冰區域（0.2C～C），明冰覆冰翼型壓力系數分布曲線與非覆冰翼型壓力系數分布曲線基本重合，而在覆冰區域（0～0.2C），明冰復雜的幾何外形導致流場十分復雜，壓力系數變化十分劇烈。

明冰覆冰狀態下模擬所得翼型周圍繞流場速度云圖如圖12所示。從圖12中可以看出，在來流風攻角為0°和1.02°時，明冰的復雜幾何外形在明冰尖角尾流區產生了兩個穩定的漩渦。當來流風攻角≥5.13°時，由于來流攻角的增大，使得翼型下側明冰尖角正對流線方向，翼型下側明冰尖角后部漩渦消失，而翼型上側則發生完全流動分離。

明覆冰翼型繞流場渦量云圖如圖13所示，從圖13中可以看出，在來流攻角為14.24°和20.15°時，在尾流區產生周期性的漩渦脫落。

模擬所得明冰覆冰形態及霜冰覆冰形態下S809翼型升阻力系數隨來流攻角的變化曲線如圖14所示。從圖14中可以看出，霜冰覆冰形態下，雖然覆冰使得翼型前緣上下側壓差有所降低，但覆冰增大了翼型的弦長，從而部分抵消了覆冰帶來的影響，使得霜冰翼型升阻力系數與非覆冰翼型升阻力系數差別不大，且在來流風攻角較小時，霜冰的存在使得升力系數有所提高，在來流攻角為0°時，升力系數提高了36.24%。不同覆冰狀態對翼型的升力系數與阻力系數影響程度如表1所示。同樣，由于霜冰在某種意義上增加了葉片的幾何長度，使得翼型表面摩擦阻力增大，從而使得霜冰覆冰狀態下翼型的阻力系數有所增大。相反，對于明冰覆冰翼型，由于明冰的非流線幾何外形，翼型表面邊界層流動分離大大提前，流動分離也大大加劇，使得翼型的升力系數大幅減小，最大降幅達83.38%，如表1所示；且由于明冰的非流線幾何外形，在大攻角下，覆冰翼型表面壓差阻力大大增加，導致明冰翼型阻力系數大幅增加；值得一提的是，在來流風攻角為0°、1.02°和5.12°時，由于明冰尖角回流區的影響，覆冰的存在導致了負阻力的產生。

4? 結語

本文通過S809二維翼型的數值模擬，對比分析了不同湍流模型及邊界層處理方法對模擬結果的影響。通過和實驗數據對比，發現

RKE湍流模型結合標準壁面函數可以很好地預測翼型表面流動的轉捩與分離，得到的壓力系數分布曲線與實驗結果符合很好；而當流動分離較大時，SST k-ω 湍流模型對流動分離點預測有所延遲，S-A湍流模型對流動分離點預測則大幅提前。

在此基礎上，本文對比分析了霜冰與明冰兩種覆冰形態對翼型的氣動特性的影響。主要結論如下：

（1）覆冰對翼型的氣動特性的影響程度不僅與覆冰形態有關，還與來流攻角有關。

（2）霜冰對翼型的氣動外形影響并不是很大，在小攻角時，霜冰對分離點幾乎無影響，此時翼型升力系數有所增大；大攻角時，分離點有所提前，此時升力系數有所減小。所有攻角下，阻力系數都有所增大。

（3）明冰復雜的非流線幾何外形導致流動分離的提前與加劇，嚴重惡化翼型的氣動特性，使得翼型的升力系數大幅減小，阻力系數大幅增加。

（4）在來流攻角為0°、1.02°和5.12°時，受明冰尖角后回流區的影響，明冰的存在將導致負阻力的產生。

［參考文獻］

［1］陳彥. 冷面結霜研究及風力機葉片覆冰的數值模擬［D］. 北京： 清華大學， 2012.

［2］馬強，吳曉敏，陳彥. 風力機葉片覆冰的數值模擬［J］. 工程熱物理學報， 2014（4）：770-773.

［3］郝艷捧，劉國特，陽林，等.風力機組葉片覆冰數值模擬及其氣動載荷特性研究［J］.電工技術學報， 2015（10）：292-300.

［4］ADDY H E. Ice accretions and icing effects for modern airfoils［M］. Washington D.C.：National Aeronautics Administration， Glenn Research Center， 2000.

［5］張強，曹義華，胡利. 翼型表面明冰的數值模擬［J］. 航空動力學報， 2009（1）： 91-97.

［6］MYERS T G， CHARPIN J PF， THOMPSON C P. Slowly accreting ice due to supercooled water impacting on a cold surface［J］. Physics of Fluids， 2002（1）：240-256.

［7］楊勝華，林貴平. 霜冰生長過程的數值模擬［J］. 計算機工程與設計， 2010 （1）： 191-194.

［8］SANEI M， RAZAGHI R. Numerical investigation of three turbulence simulation models for S809 wind turbine air-foil［J］. Proceedings of the Institution of Mechanical Engineers， Part A： Journal of Power and Energy， 2018（8）：1037-1048.

［9］ROCHA P C， ROCHA H B， CARNEIRO F M， et al. A case study on the calibration of the k–ω SST （shear stress transport） turbulence model for small scale wind turbines designed with cambered and symmetrical airfoils［J］. Energy， 2016（97）：144-150.

［10］? BURTON T， JENKINS N， SHARPE D， et al. Wind Energy Handbook［M］. Chichester： Wiley， 2011.

［11］? SORENSEN J N， SHEN W Z. Numerical modeling of wind turbine wakes［J］. Journal of Fluids Engineering， 2002（2）：393-399.

［12］? YANG X， SOTIROPOULOS F. A new class of actuator surface models for wind turbines［J］. Wind Energy， 2018（5）：285-302.

（編輯? 何琳）

Investigation on the aerodynamic characteristic of wind turbine blade airfoils

under different ice cover patterns

Xu? Chenyuan1，? Tang? Siqi2

（1.State Grid Jurong County Electric Power Supply Company， Zhenjiang 212400， China;

2.State Grid Zhenjiang Power Supply Company， Zhenjiang 212000， China）

Abstract：? The blades of wind turbine often face the danger of icing under low temperature environment， making it crucial to study the impact of icing on the aerodynamic characteristics of the wind turbine airfoil. This paper employs a numerical simulation based on the finite volume method to investigate the aerodynamic characteristics of the S809 airfoil and the effect of icing. Specifically， this study analyzes and compares the influence of two types of ice coverings， glaze ice and rime ice， on the aerodynamic characteristics of the airfoil. The finding reveal that the rime icing has negligible influence on the aerodynamic characteristics of the airfoil， whereas the glaze icing will significantly deteriorate the aerodynamic characteristics of the airfoil and may resulting negative drag under certain conditions.

Key words： wind turbine blade; aerodynamic characteristics; icing effects; numerical simulation