（華北電力大學(xué) 電站能量傳遞轉(zhuǎn)化與系統(tǒng)教育部重點(diǎn)實(shí)驗(yàn)室）

0 引言

風(fēng)力機(jī)的工作環(huán)境復(fù)雜。例如，高原、寒冷地區(qū)以及山脊、山頂處的風(fēng)資源通常較為豐富，具有很好的開發(fā)價(jià)值。然而這些地方海拔高、溫度低、濕度大。工作在這些地區(qū)的風(fēng)力機(jī)很容易遇到含有過(guò)冷液態(tài)水滴的空氣或是雨雪天氣，最終致使葉片出現(xiàn)覆冰現(xiàn)象。覆冰會(huì)改變?nèi)~片的幾何形狀，可能會(huì)降低葉片的氣動(dòng)性能，導(dǎo)致風(fēng)力機(jī)組出力下降。其次，覆冰會(huì)增加葉片的質(zhì)量，改變?nèi)~片的載荷分布，對(duì)機(jī)組的安全運(yùn)行造成嚴(yán)重危害甚至出現(xiàn)運(yùn)營(yíng)事故[1]。因此，研究風(fēng)力機(jī)結(jié)冰問(wèn)題對(duì)寒冷、潮濕地區(qū)的風(fēng)力機(jī)設(shè)計(jì)具有重要的意義。

結(jié)冰過(guò)程與葉片的迎風(fēng)面積，以及空氣溫度、風(fēng)速、空氣中水含量和液滴直徑有關(guān)，涉及到多相流與傳熱，過(guò)程復(fù)雜。實(shí)驗(yàn)是研究風(fēng)力機(jī)葉片結(jié)冰問(wèn)題的直接手段。結(jié)冰的實(shí)驗(yàn)研究通常采用自然環(huán)境下的外場(chǎng)實(shí)驗(yàn)或者特殊的冰風(fēng)洞實(shí)驗(yàn)。Bose對(duì)某小型風(fēng)力機(jī)在自然環(huán)境下的結(jié)冰過(guò)程進(jìn)行了研究，獲得了冰型輪廓和質(zhì)量，并分析了其對(duì)風(fēng)力機(jī)性能的影響[2]。Jasinski等采用低溫風(fēng)洞研究了結(jié)冰對(duì)S908翼型段氣動(dòng)性能的影響，發(fā)現(xiàn)不同的冰型對(duì)風(fēng)力機(jī)的功率的影響有明顯不同[3]。Hochart等人采用風(fēng)洞實(shí)驗(yàn)研究了氣象條件對(duì)風(fēng)力機(jī)葉片結(jié)冰的影響，同時(shí)也研究了結(jié)冰對(duì)風(fēng)力機(jī)的氣動(dòng)性能的影響[4]。國(guó)內(nèi)的結(jié)冰風(fēng)洞較少，僅有空氣動(dòng)力研究中心、東北農(nóng)業(yè)大學(xué)、南航等幾家單位建設(shè)有結(jié)冰風(fēng)洞[5-6]，因此，相關(guān)的實(shí)驗(yàn)研究很少。東北農(nóng)業(yè)大學(xué)的李巖教授團(tuán)隊(duì)利用冬季的自然條件，采用開口式風(fēng)洞，開展了垂直軸風(fēng)力機(jī)[7]與旋轉(zhuǎn)葉片結(jié)冰的風(fēng)洞實(shí)驗(yàn)[8]。

相比受制于風(fēng)洞條件的試驗(yàn)研究，翼型和葉片結(jié)冰的數(shù)值模擬方法應(yīng)用更為廣泛。從上世紀(jì)50年代開始，國(guó)外研究者采用數(shù)值模擬方法對(duì)風(fēng)力機(jī)翼型結(jié)冰問(wèn)題進(jìn)行了研究。Messinger[9]將結(jié)冰表面的傳熱過(guò)程假定為準(zhǔn)定常過(guò)程，基于能量守恒，建立了計(jì)算結(jié)冰成長(zhǎng)率的Messinger模型。該模型在飛機(jī)機(jī)翼的結(jié)冰研究中得到了廣泛使用。MacArthur[10]建立了計(jì)算二維翼型的數(shù)值方法。他將模擬過(guò)程分為流場(chǎng)求解、水滴運(yùn)動(dòng)軌跡求解和結(jié)冰熱力學(xué)模型求解三個(gè)步驟，在三個(gè)步驟之間進(jìn)行迭代求解。這一思路將流場(chǎng)計(jì)算、水滴運(yùn)動(dòng)和傳熱進(jìn)行了解耦，簡(jiǎn)化了求解過(guò)程，被當(dāng)前的結(jié)冰模擬軟件和程序所廣泛使用，例如NASA的LEWICE程序。Kwon和Sankar將上述方法擴(kuò)展到了三維，并在計(jì)算中考慮了三維真實(shí)冰型[11]。Homola和Nicklasson[12]等利用CFD方法模擬了NREL 5MW風(fēng)力機(jī)葉片表面的結(jié)冰過(guò)程。Alberts的研究發(fā)現(xiàn)濕度大、溫度低于零攝氏度時(shí)，風(fēng)力機(jī)葉片結(jié)冰的概率會(huì)增大[13]。我國(guó)研究人員在翼型結(jié)冰的數(shù)值模擬研究也取得了很多成果。朱程香等人和胡良權(quán)等人采用RANS模擬分別研究了NACA63618翼型與S809翼型表面的結(jié)冰情況。易賢和王開春等針對(duì)過(guò)冷水滴在風(fēng)力機(jī)葉片表面撞擊并結(jié)冰的現(xiàn)象，建立了適用于風(fēng)力機(jī)葉片結(jié)冰過(guò)程中水滴收集率計(jì)算的三維數(shù)值計(jì)算方法，并開發(fā)了計(jì)算程序[14]。蔣維等人利用商用CFD軟件對(duì)風(fēng)力機(jī)葉片和風(fēng)輪的結(jié)冰過(guò)程進(jìn)行了數(shù)值模擬[15-16]。

以上的風(fēng)力機(jī)葉片翼型結(jié)冰的相關(guān)研究中，并沒有考慮風(fēng)力機(jī)工作的來(lái)流條件具有一定的隨機(jī)性。其來(lái)流風(fēng)速和風(fēng)向等參數(shù)均是服從一定概率分布的隨機(jī)變量，如高斯分布和Weibull分布。由于結(jié)冰的冰型與來(lái)流條件直接相關(guān)，結(jié)冰后葉片或翼型的幾何表面的變化具有一定的不確定性，進(jìn)而對(duì)其氣動(dòng)特性的影響也具有不確定性。在風(fēng)力機(jī)氣動(dòng)性能處于敏感區(qū)間時(shí)，例如在翼型失速點(diǎn)附近，結(jié)冰對(duì)風(fēng)力機(jī)氣動(dòng)性能影響會(huì)有顯著的變化。

因此，本文對(duì)考慮來(lái)流風(fēng)況隨機(jī)變化情況下風(fēng)力機(jī)翼型的結(jié)冰過(guò)程進(jìn)行模擬，并研究這種不確定性在流場(chǎng)中的傳播，以及翼型氣動(dòng)特性的響應(yīng)變化。由于傳統(tǒng)的采樣法，如蒙特卡洛法，需要大量的樣本，計(jì)算量大，效率低。因此本文采用基于譜展開的多項(xiàng)式混沌法。根據(jù)與CFD模擬耦合方式的不同，多項(xiàng)式混沌方法分為嵌入式和非嵌入式方法[17]。王曉東等采用概率配點(diǎn)法(NIPRCM，Non-intrusiveProbabilistic Collocation method)構(gòu)造多項(xiàng)式混沌，研究了在工作條件含有不確定性參數(shù)時(shí)翼型和轉(zhuǎn)子氣動(dòng)特性的變化，以及不確定性擾動(dòng)在流場(chǎng)中的傳播[18-19]。研究結(jié)果顯示該方法可以在保證計(jì)算精度一定的情況下大大減少確定性計(jì)算點(diǎn)的數(shù)量，從而提高了計(jì)算效率。

本文首先對(duì)確定性來(lái)流參數(shù)下翼型結(jié)冰進(jìn)行了模擬，獲得確定性參數(shù)下的冰型和翼型氣動(dòng)參數(shù)的變化。在此基礎(chǔ)上，采用概率配點(diǎn)法結(jié)合CFD方法模擬了來(lái)流風(fēng)速和攻角發(fā)生隨機(jī)波動(dòng)時(shí)，翼型結(jié)冰冰型的變化和統(tǒng)計(jì)特性，及翼型氣動(dòng)特性對(duì)此不確定性參數(shù)的響應(yīng)特性。

1 翼型結(jié)冰過(guò)程的確定性CFD模擬

本文研究的模型為NREL S825翼型，其型線如圖1所示，圖中x，y為型線坐標(biāo)，c為弦長(zhǎng)，在本文研究中，弦長(zhǎng)設(shè)置為1m。

圖1 S825翼型幾何型線Fig.1 Geometry profile of S825 airfoil

1.1 結(jié)冰模擬計(jì)算設(shè)置

采用FENSAP求解器求解定常雷諾平均N-S方程，離散方法采用中心差分格式，湍流模型為Spalart-Allmaras方程模型。給定來(lái)流風(fēng)速為28.6m/s，大氣壓力，氣溫為-10.3℃，給定空氣中液態(tài)水含量（LWC）為7.5×10-4kg/m3，水滴直徑20μm，密度為1 000kg/m3。翼型表面所結(jié)冰的溫度為-5℃，結(jié)冰時(shí)間設(shè)為420s。給定翼型壁面為絕熱無(wú)滑移條件。

1.2 不同攻角下冰型

圖2給出了420s后，0°～12°攻角下翼型表面結(jié)冰后的廓線，圖中AoA表示攻角，clean airfoil表示沒有結(jié)冰的翼型，iced表示結(jié)冰的翼型，0deg表示0度攻角，其他數(shù)字以此類推。從圖2可看出在其他各個(gè)參數(shù)相同的條件下，結(jié)冰現(xiàn)象主要出現(xiàn)在翼型前緣處，并且結(jié)冰形狀多帶有棱角狀，因此可以斷定所結(jié)冰為明冰。

圖2 不同攻角下翼型結(jié)冰冰型Fig.2 Icing profiles of airfoil under different angles of attack

隨攻角的不同，引起翼型表面結(jié)冰的形狀、位置也有所不同。當(dāng)攻角為0°時(shí)，水滴撞擊到翼型后，沒有被立刻凍結(jié)，由于翼型前緣的“阻擋”，水滴會(huì)分別沿著翼型上、下表面繼續(xù)向后流動(dòng)，因此翼型吸力面和壓力面都有冰層的分布，形狀接近對(duì)稱。當(dāng)攻角為2°時(shí)，沒有被凍結(jié)的水滴沿著翼型吸力面繼續(xù)向后流動(dòng)，結(jié)冰分布主要集中在翼型上表面。但是隨著攻角的繼續(xù)增大，翼型前緣與壓力面成為水滴主要撞擊的位置，冰層的位置向前緣和壓力面移動(dòng)。由于壓力面前緣附近存在逆壓梯度，沒有凍結(jié)的水滴無(wú)法沿著壓力面繼續(xù)流動(dòng)，因此翼型壓力面前緣附近出現(xiàn)了較厚的冰層。

2 結(jié)冰對(duì)翼型氣動(dòng)力的影響

2.1 計(jì)算設(shè)置

結(jié)冰模擬提供了結(jié)冰后翼型的幾何。為了研究結(jié)冰對(duì)翼型氣動(dòng)力的影響，采用CFD對(duì)結(jié)冰后翼型的繞流場(chǎng)進(jìn)行模擬。NUMECA的AutoGrid5模塊生成結(jié)構(gòu)化網(wǎng)格，每個(gè)攻角的網(wǎng)格均采用相同的拓?fù)浣Y(jié)構(gòu)，網(wǎng)格總數(shù)均為14萬(wàn)左右。計(jì)算采用FINETM/TURBO求解器求解定常雷諾平均N-S方程，采用二階中心差分格式，湍流模型為Spalart-Allmaras方程模型。計(jì)算域外邊界設(shè)置為遠(yuǎn)場(chǎng)條件，給定來(lái)流速度、靜壓和靜溫等，壁面為絕熱無(wú)滑移條件。計(jì)算全局殘差下降6個(gè)量級(jí)以上為收斂。圖3分別給出了2°，5°，8°和12°攻角下前緣的冰型翼型與網(wǎng)格局部放大圖。

圖3 翼型前緣冰型與網(wǎng)格局部放大圖Fig.3 Ice profile and closed view of the mesh close to leading edge

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

圖4 翼型結(jié)冰對(duì)氣動(dòng)力影響Fig.4 Influence of Ice accretion of airfoil on the aerodynamic coefficients

圖4給出了翼型表面結(jié)冰后升阻力系數(shù)的變化比較，圖中橫坐標(biāo)。對(duì)比看出，結(jié)冰問(wèn)題顯著降低了該翼型的氣動(dòng)性能。結(jié)冰后，各個(gè)攻角下的升力系數(shù)明顯下降，在11°攻角附近最大下降40%左右。阻力系數(shù)則有明顯增加，最大差別也出現(xiàn)在11°攻角左右。但翼型結(jié)冰后，失速點(diǎn)略微后移，在12°攻角時(shí)，升力系數(shù)達(dá)到最大。這是因?yàn)橛捎诓还忭槺停沟靡硇臀γ媲熬壐浇桶l(fā)生了流動(dòng)分離。圖5給出了5°攻角與12°攻角下，壓力系數(shù)的分布云圖，cp表示壓力系數(shù)。可以看到，隨著攻角的增大，分離區(qū)的范圍擴(kuò)大。攻角達(dá)到12°時(shí)，分離區(qū)擴(kuò)展到吸力面三分之一弦長(zhǎng)處。而在潔凈翼型的分離流動(dòng)從吸力面尾緣處開始。在12°攻角時(shí)分離流覆蓋了吸力面后半部分。

圖5 前緣附近壓力系數(shù)云圖Fig.5 Contours of pressure coefficient close to leading edge

3 考慮攻角和風(fēng)速不確定性的CFD模擬

本節(jié)采用概率配置點(diǎn)法來(lái)研究當(dāng)攻角和來(lái)流風(fēng)速發(fā)生不確定性波動(dòng)時(shí)，對(duì)結(jié)冰形成的影響以及結(jié)冰對(duì)氣動(dòng)性能影響，有關(guān)非嵌入式概率配點(diǎn)法的詳細(xì)介紹可以參考文獻(xiàn)[19-21]。

3.1 計(jì)算設(shè)置

仍選取5°和12°攻角為研究工況，分別假設(shè)來(lái)流風(fēng)速、攻角分別為服從高斯分布的單個(gè)不確定性變量，且風(fēng)速和攻角同時(shí)為服從高斯分布的變量。高斯分布的統(tǒng)計(jì)均值與確定性算例的值一致，標(biāo)準(zhǔn)差為設(shè)置為均值的10%，計(jì)算采用二階概率配置點(diǎn)法，配置點(diǎn)見表1。

表1 不確定性計(jì)算參數(shù)設(shè)置Tab.1 Settings of stochastic calculations

3.2 計(jì)算結(jié)果及分析

1）速度不確定性計(jì)算結(jié)果

圖6 速度不確定性對(duì)翼型結(jié)冰冰型的影響Fig.6 Influence of uncertainty of wind speed on the ice profile of airfoil

圖6給出了Case 1和Case 2兩個(gè)工況下，即存在速度不確定性時(shí)翼型前緣結(jié)冰冰型的變化。圖中給出了三個(gè)配置點(diǎn)的計(jì)算結(jié)果與統(tǒng)計(jì)均值的結(jié)果，mean value表示統(tǒng)計(jì)均值?？梢园l(fā)現(xiàn)，其中從圖中可以看出在小攻角下，當(dāng)來(lái)流風(fēng)速發(fā)生波動(dòng)時(shí)，結(jié)冰的位置和形狀變化不大，但冰層的厚度隨著來(lái)流風(fēng)速的增加而明顯增大。統(tǒng)計(jì)均值的結(jié)果與確定性計(jì)算的結(jié)果幾乎完全重合，介于大風(fēng)速和小風(fēng)速情況之間。在12°攻角下，結(jié)冰的位置變化不大，但是冰型有了一定的變化。大風(fēng)速下冰層更向壓力面?zhèn)榷逊e，前緣處的冰層與其他風(fēng)速下的冰層厚度基本相同。因此，在大攻角下，冰型對(duì)來(lái)流速度的隨機(jī)變化更敏感。

圖7給出了5°和12°兩個(gè)特征工況下升阻力系數(shù)的統(tǒng)計(jì)均值及置信度為99%的置信區(qū)間。圖中deterministic表示確定性計(jì)算，non-deterministic表示不確定性計(jì)算。兩個(gè)特征工況下升力系數(shù)的統(tǒng)計(jì)均值與確定性結(jié)果基本重合，置信區(qū)間關(guān)于均值對(duì)稱。5°攻角下，升力系數(shù)的置信區(qū)間上下變化約在0.08。隨著風(fēng)速的攻角的增加，上下變化達(dá)到0.14。而阻力系數(shù)的統(tǒng)計(jì)均值略高于確定性結(jié)果，置信區(qū)間關(guān)于均值基本對(duì)稱。5°攻角下，置信區(qū)間上下變化約在0.008。12°攻角時(shí)，置信區(qū)間上下變化達(dá)到了0.012。

圖7 速度不確定性對(duì)升阻力特性的影響（99%置信區(qū)間）Fig.7 Influence of uncertainty of wind speed on the lift and drag coefficients(99%confidence interval)

2）攻角不確定性計(jì)算結(jié)果

圖8給出了5°和12°工況下攻角不確定性對(duì)翼型結(jié)冰狀況的影響。從圖中可以看出當(dāng)攻角發(fā)生波動(dòng)時(shí)，即使在小攻角下，翼型上結(jié)冰的位置和冰型明顯發(fā)生了變化。結(jié)冰幾何外形的不確定性計(jì)算統(tǒng)計(jì)均值與確定性計(jì)算也不再重合，產(chǎn)生了一定的偏差。在下攻角下，冰型的偏差主要體現(xiàn)在最大厚度位置逐漸向翼型壓力面移動(dòng)。在大攻角下，冰型的偏差主要體現(xiàn)在冰層覆蓋面積。隨著攻角增大，冰層厚度變化趨于平緩。

圖8 攻角不確定性對(duì)翼型結(jié)冰冰型的影響Fig.8 Influence of uncertainty of angle of attack on ice profile of airfoil

圖9給出了來(lái)流攻角具有不確定性時(shí)兩個(gè)特征工況下升阻力系數(shù)的統(tǒng)計(jì)均值及99%置信區(qū)間。從圖中可以看出，在5°攻角下，升力系數(shù)的統(tǒng)計(jì)均值與確定性計(jì)算結(jié)果基本重合，其置信區(qū)間關(guān)于均值對(duì)稱，上下變化范圍約在0.04，明顯小于Case 1的情況。而在12°攻角下，統(tǒng)計(jì)均值比確定性計(jì)算結(jié)果低約0.06，置信區(qū)間不再關(guān)于均值對(duì)稱，置信區(qū)間的上限略大于確定性計(jì)算結(jié)果，而下限比確定性計(jì)算降低了約0.23，有明顯的降低。因此，在接近失速的大攻角下，確定性計(jì)算低估了結(jié)冰對(duì)翼型升力的影響。確定性計(jì)算結(jié)果已表明在給定來(lái)流風(fēng)速的條件下，在12°攻角時(shí)翼型升力系數(shù)已接近最大（結(jié)合圖4），來(lái)流攻角發(fā)生微小波動(dòng)都會(huì)造成升力系數(shù)的下降。這與圖7的置信區(qū)間變化有明顯不同，這是因?yàn)樵贑ase 1中，隨著來(lái)流風(fēng)速不同，升力曲線不再是28.6m/s來(lái)流對(duì)應(yīng)的升力曲線。而來(lái)流一定的情況下（圖8），升力曲線只有一條。從圖8(b)中可以看到，攻角的不確定性波動(dòng)對(duì)阻力系數(shù)影響較小。大攻角下阻力系數(shù)的置信區(qū)間的變化大于小攻角下的置信區(qū)間的變化。

圖9 攻角不確定性對(duì)升阻力特性的影響（99%置信區(qū)間）Fig.9 Influence of uncertainty of angle of attack on the lift and drag coefficients(99%confidence interval)

3）速度和攻角同為不確定性變量的計(jì)算結(jié)果

Case 3主要比較同時(shí)存在多個(gè)獨(dú)立的隨機(jī)變量時(shí)，各變量的作用是否存在相互影響。圖10比較了單個(gè)不確定量作用結(jié)果線性疊加后與兩個(gè)不確定性量同時(shí)作用的結(jié)果下升力系數(shù)和阻力系數(shù)的標(biāo)準(zhǔn)差。由圖可見，兩個(gè)不確定性量的同時(shí)作用的結(jié)果略小于直接線性疊加的結(jié)果。這說(shuō)明速度和攻角雖然是獨(dú)立的隨機(jī)變量，但二者之間存在物理的耦合作用，削弱了不確定性對(duì)翼型氣動(dòng)特性結(jié)果的影響。

圖10 不確定性量獨(dú)立作用與同時(shí)作用影響比較Fig.10 Comparison of effects of individual and multi-uncertainties

4 結(jié)論

本文采用CFD模擬研究了S825翼型結(jié)冰對(duì)翼型升阻力特性的影響，采用研究了來(lái)流風(fēng)速和攻角的隨機(jī)變化對(duì)翼型表面結(jié)冰形狀以及結(jié)冰后翼型的氣動(dòng)性能的影響。結(jié)果表明：

1）與未結(jié)冰的潔凈翼型相比，結(jié)冰后的翼型氣動(dòng)性能嚴(yán)重下降。大攻角下，結(jié)冰更容易在葉片前緣與壓力面前緣堆積。

2）在來(lái)流風(fēng)速存在隨機(jī)變化時(shí)，冰型的位置和形狀變化不大。來(lái)流攻角發(fā)生變化時(shí)，冰型的位置和形狀會(huì)發(fā)生明顯變化。

3）大攻角下，來(lái)流風(fēng)速的不確定波動(dòng)對(duì)結(jié)冰后翼型的氣動(dòng)性能影響比來(lái)流攻角的不確定性波動(dòng)影響大。在大攻角下，確定性計(jì)算會(huì)低估攻角對(duì)結(jié)冰的影響，進(jìn)而低估對(duì)升力系數(shù)的影響。風(fēng)速和攻角的耦合作用削弱了不確定性對(duì)翼型氣動(dòng)特性的影響。