收稿日期：2022-09-14

基金項目：中國華能集團清潔能源技術(shù)研究院青年基金（TO-21-CERI01）；上海市科技創(chuàng)新行動計劃（20dz1205300）

通信作者：杜靜宇（1995—），女，博士、工程師，主要從事相變傳熱方面的研究。 jy_du@qny.chng.com.cn

DOI：10.19912/j.0254-0096.tynxb.2022-1392 文章編號：0254-0096（2023）12-0298-08

摘 要：該文首先建立歐拉法水滴收集系數(shù)模型和風(fēng)力機葉片翼型結(jié)冰模型，并以NACA0012翼型結(jié)冰算例驗證該模型的有效性，而后開展NREL 5 MW風(fēng)力機葉片外部翼型（NACA64618）結(jié)冰數(shù)值模擬，對風(fēng)速、液態(tài)水含量和水滴直徑等參數(shù)在結(jié)冰速率的影響方面進行了討論。數(shù)值模擬結(jié)果表明：隨攻角增大，NACA64618翼型結(jié)冰區(qū)域移向壓力面，結(jié)冰厚度增加（最大增加量約為167%）；增大風(fēng)速、液態(tài)水含量和水滴直徑均加劇NACA64618翼型結(jié)冰，結(jié)冰區(qū)域擴大、結(jié)冰量與結(jié)冰厚度增加。

關(guān)鍵詞：風(fēng)力機；葉片；翼型；結(jié)冰；數(shù)值模擬；水滴撞擊特性

中圖分類號：TK83"nbsp;"""""""" """"""" ""文獻標(biāo)志碼：A

0 引 言

高海拔、低溫區(qū)域風(fēng)能資源更為充裕，風(fēng)力機組常安裝在高地山區(qū)［1］。但這些地區(qū)氣候惡劣、濕度大、常出現(xiàn)凍雨凍雪等現(xiàn)象，導(dǎo)致風(fēng)力機葉片極易結(jié)冰［2-4］。這一現(xiàn)象不僅使葉片氣動外形發(fā)生改變，令氣動性能惡化，阻力系數(shù)升高、升力系數(shù)降低、失速攻角提前［5-10］，還會降低機組運行穩(wěn)定性［11］。

國內(nèi)外已有許多學(xué)者開展了風(fēng)力機葉片結(jié)冰的相關(guān)實驗與仿真研究。

在實驗研究方面，Barber等［12］研究不同葉片結(jié)冰位置對風(fēng)力機氣動性能的影響程度，發(fā)現(xiàn)葉尖區(qū)域結(jié)冰影響最大，因此對葉片外部結(jié)冰研究是相關(guān)研究的重點。李巖等［13］實驗研究垂直軸風(fēng)力機結(jié)冰性質(zhì)，考察了迎風(fēng)面積對葉片結(jié)冰量的影響。文獻［14］開展NREL Phase VI型風(fēng)力機結(jié)冰實驗，并使用LEWICE軟件對冰形進行預(yù)測，指出預(yù)測冰形與實驗冰形仍有差異。這些研究指明了風(fēng)力機葉片結(jié)冰關(guān)鍵部位，實驗研究了結(jié)冰的部分影響參數(shù)，是開展結(jié)冰數(shù)值模擬的重要參照。

在仿真研究方面，易賢等［15］針對水平軸風(fēng)力機開發(fā)三維結(jié)冰模擬程序，指出葉尖結(jié)冰量最大，說明葉片結(jié)冰應(yīng)更關(guān)注葉尖部位，這與Barber等人的實驗結(jié)論一致。Pedersen等［16-17］基于ANSYS Fluent的多相流模型和動網(wǎng)格重構(gòu)幾何技術(shù)開展二維翼型結(jié)冰研究，并模擬了結(jié)霜冰的情況。鄧曉湖等［18］建立水滴收集系數(shù)的拉格朗日法計算模型，并模擬了翼型結(jié)霜冰外形。Homola等［19］使用TURBICE軟件對結(jié)冰后的翼型外形特點及其氣動性能變化進行分析。任鵬飛等［20］針對NREL Phase VI型風(fēng)力機開展葉片結(jié)冰數(shù)值模擬，分析了霜冰對該風(fēng)力機性能的影響。這些研究分析了風(fēng)力機葉片結(jié)冰的基本原理，建立了葉片結(jié)冰數(shù)值模擬的基本流程方法，并完成了部分翼型的結(jié)冰模擬，探討了結(jié)冰對翼型氣動性能的影響；但目前的研究尚未形成較通用的結(jié)冰數(shù)值仿真模型，對于結(jié)冰的影響因素研究還不夠深入。

本文首先采用歐拉法建立水滴收集系數(shù)計算模型，在此基礎(chǔ)上建立翼型結(jié)冰模型，并以NACA0012翼型結(jié)霜冰驗證此方法的有效性，而后開展NREL 5 MW風(fēng)力機葉片外部翼型（NACA64618）結(jié)冰模擬，最后分析風(fēng)速、水滴直徑和液態(tài)水含量等參數(shù)對結(jié)冰的影響。

1 水滴收集系數(shù)模型建立與驗證

1.1 水滴撞擊特性

在寒冷潮濕的條件下，濕空氣中攜帶的過冷液滴易在風(fēng)力機運行的葉片前緣發(fā)成撞擊。通常，把水滴在葉片/翼型前緣的撞擊范圍、撞擊水量和水量的分布等統(tǒng)稱為水滴撞擊特性［21］。水滴撞擊特性的求解是后續(xù)結(jié)冰模擬的前提。

水滴收集系數(shù)[β]是表征水滴撞擊特性的關(guān)鍵參數(shù)，其定義為翼型表面局部區(qū)域?qū)嶋H水滴撞擊量與最大可能水滴撞擊量之比。水滴收集系數(shù)體現(xiàn)翼型表面的水滴收集能力，反映水滴在翼型表面的分布情況。

1.2 歐拉法水滴收集系數(shù)模型建立

1.2.1 水滴守恒方程

水滴守恒方程的建立基于以下假設(shè)［22-23］：

1）空氣相與水滴相的相互作用?？紤]到水滴體積分?jǐn)?shù)很小，忽略水滴對空氣的作用；假設(shè)水滴相與空氣相不存在質(zhì)量與熱量交換。

2）水滴相的基本狀態(tài)。假設(shè)水滴相是均勻的；假設(shè)水滴的初始速度與自由來流速度一致。

基于以上假設(shè)，空氣相和水滴相實現(xiàn)了解耦［24］，同時可將水滴相視為連續(xù)相。設(shè)控制體內(nèi)水滴的體積分?jǐn)?shù)為α，則水滴相的連續(xù)性方程和動量方程［25］為：

[?（αρw）?t+▽?（αρwuw）=0]""" （1）

[?（αρwuw）?t+▽?（αρwuw?uw）=αρwg+αρwK（ua-uw）] （2）

式中：[ρ]——密度，kg/m3；[u]——速度，m/s；[g]——重力加速度，m/s2；[K]——空氣-水滴動量交換系數(shù)［26］。下標(biāo)：w——水滴參數(shù)；a——空氣相參數(shù)，下同。

當(dāng)水滴與翼型表面相撞時，水滴收集系數(shù)[β]采用式（3）計算［25，27］：

[β=αnα∞?un?nu∞] （3）

式中：[n]——翼表單位法向量。下標(biāo)：n——翼型表面處參數(shù)；∞——來流參數(shù)。

1.2.2 守恒方程求解

通過編寫自定義函數(shù)（user defined function，UDF），ANSYS Fluent可求解形如式（4）所示的自定義標(biāo)量（user defined scales， UDS）方程［28］。

[?ρ?k?t非穩(wěn)態(tài)項+??xiρui?k對流項-Γk??k?xi擴散項=S?k源項， k=1，…，Nscalar]" （4）

水滴守恒方程共包含3個標(biāo)量：[α、uwx]和[uwy]，這些參數(shù)可看作UDS嵌入到ANSYS Fluent中，然后通過UDS方程進行求解，為此在UDF中共使用5種宏函數(shù)，其名稱與作用如表1所示。

1.3 水滴收集系數(shù)計算方法驗證

本文以NACA0012翼型［27］算例為驗證算例，以證明水滴收集系數(shù)歐拉法計算模型的準(zhǔn)確性。工況見表2。

采用C型拓撲結(jié)構(gòu)化網(wǎng)格，取計算域為10倍翼型弦長，加密近壁面網(wǎng)格。網(wǎng)格無關(guān)性驗證以升力系數(shù)為指標(biāo)，結(jié)果如圖1所示。由圖1可知，當(dāng)網(wǎng)格總數(shù)大于1×105時，升力系數(shù)的變化很小，達到網(wǎng)格無關(guān)性要求。最終采用如圖2所示的翼型網(wǎng)格，表面共布置360個網(wǎng)格點，網(wǎng)格數(shù)量為1.24×105。

在計算水滴相之前須計算空氣相?？諝庀嗟臄?shù)值方法如下：計算采用穩(wěn)態(tài)的壓力基求解器，求解算法為二階迎風(fēng)格的SIMPLE算法，湍流模型采用k-ω SST模型；邊界條件采用速度進口、壓力出口，翼表無滑移，上下邊界對稱。待空氣相計算收斂后，再開展水滴相的求解。水滴相控制方程離散使用二階迎風(fēng)格式，以10-8為收斂殘差。邊界條件為：初始進口[α∞=0.3、uwx=135.48 m/s]和[uwy=11.85 m/s]，出口和上下邊界參數(shù)梯度為零。

圖3為本方法計算的結(jié)果與文獻中NACA0012水滴收集系數(shù)的實驗數(shù)據(jù)對比，圖3中[x]表示水滴撞擊曲線長度，[c]表示弦長，以翼型吸力面為正、壓力面為負。由圖3可知，水滴收集系數(shù)計算值與實驗值和文獻值均能較好地吻合。經(jīng)定量比較，最大水滴收集系數(shù)的誤差為3.58%，說明本文采用的水滴收集系數(shù)方法具有良好的準(zhǔn)確性。

2 翼型結(jié)冰模型與驗證

2.1 翼型結(jié)冰模型

結(jié)冰過程中翼型表面任一控制體均遵守質(zhì)量和能量守恒，其示意圖如圖4所示。

2.1.1 結(jié)冰過程中的質(zhì)量守恒

由圖4a可知，控制體質(zhì)量守恒方程可表示為：

[Mimp+Min=Meva+Mout+Mice]" （5）

式中：[Mimp]——撞擊水量，kg；[Min]——控制體流入水量， kg；[Meva]——蒸發(fā)水量，kg；[Mout]——控制體流出水量， kg；[Mice]——結(jié)冰水量，kg。式（5）中各子項計算方法參見表3。

表中：[u∞]——來流速度，m/s；[ρLWC]——液態(tài)水含量，g/m3；[ds]——控制體面積，m2；[hc]——對流換熱系數(shù)，W/（m2·K）；[R]——水蒸氣氣體常數(shù)，J/（kg·K）；[ρa]——空氣密度，kg/m3；[cp，a]——空氣定壓比熱容，J/（kg·K）；[Lew]——Lewis數(shù)；[psat（Ts）]——平衡溫度Ts下的飽和蒸氣壓，Pa；[φ]——空氣相對濕度；[psat（T∞）]——來流溫度[T∞]下的飽和蒸氣壓，Pa；[f]——凍結(jié)系數(shù)，[f=Mice/（Mimp+Min）]，[f]的取值范圍介于0～1之間。

2.1.2 結(jié)冰過程中的能量守恒

由圖4b可知，控制體能量守恒方程可表示為：

[Qimp+Qin+Qair=Qhtc+Qeva+Qclh+Qout+Qice]"""""" （6）

式中：[Qimp]——撞擊液態(tài)水能量，J；[Qin]——控制體流入換熱量，J；[Qair]——氣動加熱量，J；[Qhtc]——表面對流傳熱量，J；[Qeva]——蒸發(fā)熱量，J；[Qclh]——結(jié)冰放熱量，J；[Qout]——控制體流出換熱量，J；[Qice]——結(jié)冰存儲能量，J。式（6）中各子項計算方法參見表4。

2.2 翼型結(jié)冰模型的驗證

本文以NACA0012翼型［29］結(jié)霜冰算例為翼型結(jié)冰模型驗證算例，結(jié)冰條件如表5所示。需注意的是，翼型表面的實際結(jié)冰過程是瞬態(tài)的，但為減少計算消耗，本文中空氣和水滴相均采用穩(wěn)態(tài)求解器［30］。

表中：[cp，w]——水定壓比熱容，J/（kg·K）；[cp，i]——冰定壓比熱容，J/（kg·K）；[Tm]——前一個控制體的表面平衡溫度，K；[Lf]——冰融解潛熱，J/kg；[r]——附面層恢復(fù)系數(shù)，層流時，[r]可取為0.85，湍流時，[r]可取為0.89；[Le]——水的蒸發(fā)潛熱，J/kg。

圖5顯示了本文模擬結(jié)霜冰形狀與實驗和LEWICE預(yù)測結(jié)果的對比。由圖5可知，總體而言，本文模擬結(jié)果與實驗值和LEWICE預(yù)測結(jié)果均能較好地吻合，且預(yù)測精度略超

過LEWICE。為定量地分析模擬與實驗之間的誤差，選取結(jié)冰面積和最大結(jié)冰厚度作為評價指標(biāo)，計算得霜冰的結(jié)冰面積誤差達7.59%，最大結(jié)冰厚度誤差達6.64%。綜上，本文建立的結(jié)冰數(shù)值模型能準(zhǔn)確預(yù)測翼型結(jié)霜冰。

3 NACA64618翼型結(jié)冰模擬

實地觀測、實驗與計算結(jié)果均顯示外部區(qū)域是風(fēng)力機葉片的最易結(jié)冰區(qū)，且該區(qū)域?qū)τ秒姸擞绊戄^大［31-32］，因此，本節(jié)以NREL 5 MW型風(fēng)力機葉片外部翼型（NACA64618）為研究對象展開結(jié)冰模擬。

3.1 NACA64618翼型水滴收集系數(shù)計算

考慮到在實際工作中，攻角會隨風(fēng)力機變槳而發(fā)生變化，在翼型的結(jié)冰模擬中應(yīng)考慮變攻角條件。因此本文共計算4個攻角下（0°、4°、8°和12°）NACA64618翼型的水滴收集系數(shù)，計算工況如表6所示。其中風(fēng)速為60 m/s，雷諾數(shù)為4.1×106，能滿足風(fēng)力機翼型低雷諾數(shù)的要求［33］。

NACA64618翼型的網(wǎng)格劃分、數(shù)值方法與邊界條件設(shè)置等與NACA0012翼型類似，在此不再贅述。圖6為NACA64618翼型網(wǎng)格無關(guān)性驗證結(jié)果，可知網(wǎng)格數(shù)量為0.85×105"" 時已滿足無關(guān)性要求。圖7為最終使用的NACA64618翼型計算網(wǎng)格，網(wǎng)格總數(shù)為1.26×105。

圖8和表7為攻角0°～12°時NACA64618翼型的水滴收集系數(shù)計算結(jié)果。隨攻角增大，水滴收集系數(shù)峰值降低，壓力面撞擊極限增大，吸力面撞擊極限減小。當(dāng)攻角由0°增至12°時，水滴收集系數(shù)峰值降低11.76%，壓力面撞擊極限增加48.18%，吸力面撞擊極限減小84.75%。

3.2 NACA64618翼型表面結(jié)霜冰模擬

在模擬NACA64618翼型表面結(jié)霜冰時，需考慮以下問題：1）計算工況應(yīng)當(dāng)是變攻角的，以模擬實際風(fēng)力機變槳情況；2）風(fēng)場中葉片表面結(jié)冰往往持續(xù)幾個小時或幾天，在數(shù)值模擬時需對實際結(jié)冰條件進行縮比處理［34］；3）實際結(jié)冰參數(shù)是時間和空間的函數(shù)，難以精確測量。本文模擬時采用的水滴直徑、液態(tài)水含量和結(jié)冰時間等條件參考Reid等［35-36］的研究結(jié)果，NACA64618翼型結(jié)冰的計算時間為20 min，液態(tài)含水量為0.5 g/m3，其他參數(shù)與表7相同。

圖9顯示了不同攻角下NACA64618翼型表面結(jié)霜冰形狀。由圖9可知，攻角由0°增至12°時，結(jié)冰區(qū)域移向壓力面；結(jié)冰厚度增加，最大結(jié)冰厚度約增加167%。這是因為翼型表面水滴的附著與收集是結(jié)冰的基礎(chǔ)條件，水滴撞擊極限直接影響結(jié)冰范圍。由前文可知，攻角增大時，壓力面撞擊極限增大、吸力面撞擊極限減小，因此結(jié)冰區(qū)域向壓力面遷移。

3.3 結(jié)冰影響參數(shù)分析

在風(fēng)力機實際工作中，風(fēng)速（wind speed， WS）、水滴直徑（medium volume droplet diameter， MVD）和液態(tài)水含量（liquid water content， LWC）等是較易發(fā)生變化的環(huán)境條件［37］，因此本節(jié)研究這3種因素對結(jié)冰的影響。參考結(jié)冰條件如表8所示，除此之外，翼型弦長均為1 m，攻角均為0°，溫度均為260 K。

圖10～圖12為上述算例中NACA64618翼型表面結(jié)冰冰形模擬結(jié)果。由圖10～圖12可知，風(fēng)速、水滴直徑和液態(tài)水含量的增大均可使得NACA64618翼型表面結(jié)冰強度增強，結(jié)冰區(qū)域擴大、結(jié)冰厚度與結(jié)冰量增加。風(fēng)速和水滴直徑的增大使水滴動量增加，使得水滴更易撞擊到翼型前緣，導(dǎo)致水滴收集系數(shù)增大；增加液態(tài)水含量則直接導(dǎo)致水滴收集系數(shù)增大。這就使得結(jié)冰的物質(zhì)基礎(chǔ)增強，NACA64618翼型表面結(jié)冰強度隨之增強。

4 結(jié) 論

本文首先建立了水滴收集系數(shù)的歐拉法計算方法，在此基礎(chǔ)上建立了風(fēng)力機葉片翼型結(jié)冰模型，并以NACA0012翼型結(jié)霜冰算例驗證了該模型的有效性，而后開展了NREL 5 MW風(fēng)力機葉片外部翼型（NACA64618）結(jié)冰模擬，最后分析了風(fēng)速、液態(tài)水含量和水滴直徑等因素對NACA64618翼型結(jié)冰的影響。得到的主要結(jié)論如下：

1）NACA0012的水滴收集系數(shù)與翼型結(jié)冰驗證算例計算結(jié)果與實驗值及文獻值的對比，證明本文建立的水滴收集系數(shù)計算方法與翼型結(jié)冰模擬方法具有足夠的有效性與準(zhǔn)確性。

2）攻角由0°增大至12°時，NACA64618翼型表面結(jié)霜冰區(qū)域移向壓力面；結(jié)冰厚度增加，最大結(jié)冰厚度約增加167%。

3）增大風(fēng)速、液態(tài)水含量和水滴直徑均加劇NACA64618翼型表面結(jié)冰，結(jié)冰區(qū)域擴大、結(jié)冰量與結(jié)冰厚度增加。

［參考文獻］

［1］"""" BATTISTI L. Wind turbines in cold climates： icing impacts and mitigation systems［M］. Berlin： Springer， 2015.

［2］"""" PRYOR S C， BARTHELMIE R J. Climate change impacts on wind energy： a review［J］. Renewable and sustainable energy reviews， 2010， 14（1）： 430-437.

［3］"""" TAMMELIN B， SEIFERT H. Large wind turbines go into cold climate regions［C］//EWEC. Copenhagen， Denmark， 2001.

［4］"""" CATTIN" R." Icing" of" wind" turbines［C］//Onshore" Oamp;M Forum， Hamburg， Germany， 2013.

［5］"""" DROUBI M G， ISLAM S Z. Investigation of the effect of ice accretion on wind turbine aerodynamics and performance［C］// UK-China Workshop on Wind Power in Cold Climate. Mianyang， China， 2016.

［6］"""" DAVIS N N， PINSON P， HAHMANN A N， et al. Identifying and characterizing the impact of turbine icing on wind farm power generation［J］. Wind energy， 2016， 19（8）： 1503-1518.

［7］"""" BURTCH D G. Using reanalysis data for the prediction of seasonal wind turbine power losses due to icing［D］. North Dakota： University of North Dakota， 2014.

［8］"""" WALLENIUS T， LEHTOM?KI V. Overview of cold climate wind energy： challenges， solutions， and future needs［J］. WIREs energy and environment， 2016， 5（2）： 128-135.

［9］"""" SUNDEN B， WU Z. On icing and icing mitigation of wind turbine" blades" in" cold" climate［J］." Journal" of" energy resources technology， 2015， 137（5）： 051203.

［10］""" MALMSTEN J K. Wind turbine production losses in cold climate［D］. Gotland： Gotland University， 2011.

［11］""" ILINCA A. Analysis and mitigation of icing effects on wind turbines［J］. Wind turbines， 2011， 4（4）： 183-214.

［12］""" BARBER S， WANG Y， JAFARI S， et al. The impact of ice formation on wind turbine performance and aerodynamics［J］. Journal of solar energy engineering， 2011， 133（1）： 1.

［13］""" 李巖， 田川公太朗. 葉片附著物對直線翼垂直軸風(fēng)力機性能的影響［J］. 動力工程， 2009， 29（3）： 292-296.

LI Y， KOTARO T. Influence of blade attachment on performance of the straight-bladed vertical axis wind turbine［J］. Journal of power engineering， 2009， 29（3）： 292-296.

［14］""" HAN Y Q， PALACIOS J， SCHMITZ S. Scaled ice accretion experiments on a rotating wind turbine blade［J］. Journal of wind engineering and industrial aerodynamics， 2012， 109： 55-67.

［15］""" 易賢， 王開春， 馬洪林， 等. 水平軸風(fēng)力機結(jié)冰及其影響計算分析［J］. 太陽能學(xué)報， 2014， 35（6）： 1052-1058.

YI X， WANG K C， MA H L， et al. Computation of icing and its effect of horizontal axis wind turbine［J］. Acta energiae solaris sinica， 2014， 35（6）： 1052-1058.

［16］""" PEDERSEN M C， S?RENSEN H. Towards a CFD model for prediction of wind turbine power losses due to icing in cold climate［C］// International Symposium on Transport Phenomena and Dynamics of Rotating Machinery. Honolulu， HI， USA， 2016.

［17］""" PEDERSEN M C， S?RENSEN H， SWYTINK-BINNEMA N， et al. Measurements from a cold climate site in Canada： boundary conditions and verification methods for CFD icing models for wind turbines［J］. Cold regions science and technology， 2018， 147： 11-21.

［18］""" 鄧曉湖， 盧緒祥， 李錄平， 等. 水平軸風(fēng)力機葉片翼型結(jié)冰的數(shù)值模擬［J］. 能源技術(shù)， 2010， 31（5）： 266-271.

DENG X H， LU X X， LI L P， et al. Numerical simulation of airfoil ice accretion process on horizontal-axis wind turbine blade［J］. Energy technology， 2010， 31（5）： 266-271.

［19］""" HOMOLA M C， WALLENIUS T， MAKKONEN L， et al. Turbine size and temperature dependence of icing on wind turbine blades［J］. Wind engineering， 2010， 34（6）： 615-627.

［20］""" 任鵬飛， 徐宇， 宋娟娟， 等. 霜冰條件對翼型氣動性能影響數(shù)值研究［J］. 工程熱物理學(xué)報， 2014， 35（4）： 663-668.

REN P F， XU Y， SONG J J， et al. Numerical study about the influence of rime ice conditions on airfoil［J］. Journal of engineering thermophysics， 2014， 35（4）： 663-668.

［21］""" 趙秋月. 航空發(fā)動機進口支板及整流帽罩水滴撞擊特性的計算分析［D］. 上海： 上海交通大學(xué)， 2011.

ZHAO Q Y. Computational analysis of water droplet impingerment property for the inlet struct and the cone［D］. Shanghai： Shanghai Jiao Tong University， 2011.

［22］""" 楊勝華， 林貴平， 申曉斌. 三維復(fù)雜表面水滴撞擊特性計算［J］. 航空動力學(xué)報， 2010， 25（2）： 284-290.

YANG S H， LIN G P， SHEN X B. Water droplet impingement prediction for three-dimensional complex surfaces［J］. Journal of aerospace power， 2010， 25（2）： 284-290.

［23］""" 易賢， 王開春， 桂業(yè)偉， 等. 結(jié)冰面水滴收集率歐拉計算方法研究及應(yīng)用［J］. 空氣動力學(xué)學(xué)報， 2010， 28（5）： 596-601， 608.

YI X， WANG K C， GUI Y W， et al. Study on Eulerian method for icing collection efficiency computation and its application［J］. Acta aerodynamica sinica， 2010， 28（5）： 596-601， 608.

［24］""" KIM J， DENNIS P G， SANKAR L， et al. Ice accretion modeling using an eulerian approach for droplet impingement［C］//51st AIAA Aerospace Sciences Meeting including the New Horizons Forum and Aerospace Exposition. Grapevine， TX， USA， 2013： 246.

［25］""" WIROGO S， SRIRAMBHATLA S. An eulerian method to calculate the collection efficiency on two and three dimensional bodies［C］//41st Aerospace Sciences Meeting and Exhibit. Reno， NV， USA， 2003： 1073.

［26］""" 胡良權(quán)， 陳進格， 沈昕， 等. 風(fēng)力機葉片結(jié)冰水滴收集系數(shù)計算［J］. 太陽能學(xué)報， 2020， 41（3）： 22-28.

HU L Q， CHEN J G， SHEN X， et al. Calculation of water droplets collection efficiency on wind turbine blades［J］. Acta energiae solaris sinica， 2020， 41（3）： 22-28.

［27］""" TONG X L， LUKE E. Eulerian simulations of icing collection efficiency using a singularity diffusion model［C］//43rd AIAA Aerospace Sciences Meeting and Exhibit， Reno， NV， USA， 2005： 1246.

［28］""" ANSYS Inc. ANSYS FLUENT UDF Manual. ANSYS FLUENT 14.0［M］ California：ANSYS Inc， 2011.

［29］""" SHIN J， BOND T. Results of an icing test on a NACA 0012 airfoil in the NASA Lewis Icing Research Tunnel［C］//30th Aerospace Sciences Meeting and Exhibit. Reno， NV， USA， 1992： 647.

［30］""" HASANZADEH K， LAURENDEAU E， PARASCHIVOIU I. Quasi-steady convergence of multistep Navier-Stokes icing simulations［J］. Journal of aircraft， 2013， 50（4）： 1261-1274.

［31］""" 戰(zhàn)培國. 國外寒冷地區(qū)風(fēng)力機結(jié)冰問題研究［J］. 航空科學(xué)技術(shù)， 2016， 27（2）： 1-6.

ZHAN P G. Review of the wind turbine icing in overseas cold regions［J］. Aeronautical science amp; technology， 2016， 27（2）： 1-6.

［32］""" BURTON T， SHARPE D， JENKINS N， et al. Wind energy" handbook［M］." Chichester：" John" Wiley" amp;" Sons Ltd， 2001.

［33］""" 劉磊. 風(fēng)力機葉片非定常氣動特性的研究［D］. 北京： 中國科學(xué)院研究生院（工程熱物理研究所）， 2012.

LIU L. Research on the unsteady aerodynamic characteristics"" of"" wind"" turbine"" blades［D］."" Beijing： Institute of Engineering Thermophysics， Chinese Academy of Sciences， 2012.

［34］""" VILLALPANDO F， REGGIO M， ILINCA A. Prediction of ice accretion and anti-icing heating power on wind turbine blades using standard commercial software［J］. Energy， 2016， 114： 1041-1052.

［35］""" REID T， BARUZZI G， OZCER I， et al. FENSAP-ICE simulation of icing on wind turbine blades， part 1： performance""" degradation［C］//51st"""" AIAA"""" Aerospace Sciences Meeting including the New Horizons Forum and Aerospace Exposition. Grapevine， TX， USA， 2013： 750.

［36］""" REID T， BARUZZI G， OZCER I， et al. FENSAP-ICE simulation of icing on wind turbine blades， part 2： ice protection"" system"" design［C］//51st"" AIAA"" Aerospace Sciences Meeting including the New Horizons Forum and Aerospace Exposition， Grapevine， TX， USA， 2013： 751.

［37］""" 胡良權(quán). 水平軸風(fēng)力機葉片結(jié)冰問題研究［D］. 上海： 上海交通大學(xué)， 2019.

HU L Q. Study of horizontal axis wind turbine blade icing［D］. Shanghai： Shanghai Jiao Tong University， 2019.

ICING SIMULATION OF AIRFOIL OF NREL 5 MW WIND TURBINE BLADE

Du Jingyu1，Hu Liangquan2，Ren Xin1，Shen Xin2

（1. Huaneng Clean Energy Research Institute， Beijing 102209， China；

2. School of Mechanical Engineering， Shanghai Jiao Tong University， Shanghai 200240， China）

Abstract：In this paper， the Euler method water droplet collection coefficient model and the airfoil icing model are established， and the NACA0012 airfoil icing calculation example is used to verify the effectiveness of the model. The numerical simulation of NREL 5 MW wind turbine blade outer airfoil icing （NACA64618） is carried out. The effects of wind speed， liquid water content and droplet diameter on the icing of the airfoil are studied. The numerical simulation results show that with the increase of the angle of attack， the icing area of NACA64618 airfoil moves towards the pressure surface， and the icing thickness increases （the maximum increase is about 167%）. The increase of wind speed， liquid water content and droplet diameter can all aggravate the icing of NACA64618 airfoil， and increase the icing area， icing amount and icing thickness.

Keywords：wind turbines； blades； airfoil； freezing； numerical simulation； water droplets impingement characteristics