于佳鑫， 王曉東， 陳江濤， 吳曉軍， 康 順

(1.電站能量傳遞轉(zhuǎn)化與系統(tǒng)教育部重點(diǎn)實驗室 (華北電力大學(xué))， 北京 102206；2.中國空氣動力研究與發(fā)展中心， 四川 綿陽 621000)

0 引 言

當(dāng)前風(fēng)力機(jī)葉片的尺寸越來越大，2030年風(fēng)輪直徑有望達(dá)275 m[1]，這必然會增加葉片的柔性。柔性葉片更易發(fā)生不穩(wěn)定振蕩。從經(jīng)濟(jì)性方面來看，顫振會導(dǎo)致葉片壽命降低，影響整個發(fā)電系統(tǒng)；從安全性角度來說，可能會引發(fā)葉片扭轉(zhuǎn)受損甚至安全事故。因此，在設(shè)計葉片時則不得不考慮振蕩的效應(yīng)[2]。

為了增加風(fēng)力機(jī)葉片氣彈性能，從氣動性能角度出發(fā)，學(xué)者們對提高氣動性能參數(shù)的魯棒性做了很多研究[8, 9]；從降低結(jié)構(gòu)疲勞破壞的角度來看，學(xué)者們對提升結(jié)構(gòu)材料性能[10]和結(jié)構(gòu)布局[11,12]也做了一些分析，然而，實際環(huán)境中的不確定因素可能恰恰是引發(fā)葉片不穩(wěn)定顫振的原因[13]。文獻(xiàn)[14]表明氣動力和結(jié)構(gòu)特性的隨機(jī)性影響著風(fēng)力機(jī)葉片不穩(wěn)定振蕩的發(fā)生，且結(jié)構(gòu)特性的隨機(jī)性會降低顫振的臨界速度。文獻(xiàn)[15]分析葉片剛度、材料等結(jié)構(gòu)不確定性對于葉片固有頻率的影響。因此，在翼型及葉片的設(shè)計階段就需要考慮不確定因素對其穩(wěn)定狀態(tài)的影響。

文獻(xiàn)[16]用蒙特卡洛方法(MC)研究風(fēng)力機(jī)葉片的隨機(jī)振蕩問題，分析了顫振發(fā)生的概率及臨界失效點(diǎn)，然而，MC方法效率很低，計算量巨大?；谧V分析的多項式混沌法效率高，通過與CFD進(jìn)行非嵌入式耦合，已經(jīng)被廣泛應(yīng)用于隨機(jī)氣動問題的求解[17]。文獻(xiàn)[18]采用基于稀疏網(wǎng)格的多項式混沌法(SGPC)對多維隨機(jī)幾何進(jìn)行氣動不確定性分析，相較于概率配置點(diǎn)法[19]在變量維數(shù)較多時計算量指數(shù)增加的情況，基于稀疏網(wǎng)格的多項式混沌法在保持結(jié)果精度的同時，計算效率也取得了顯著的提高。

因此本文針對翼型的俯仰沉浮耦合振蕩形式，采用能量法對顫振邊界進(jìn)行研究與分析。探究翼型運(yùn)行過程中誘導(dǎo)不穩(wěn)定顫振的主要不確定因素，并基于SGPC的不確定性分析方法量化它們的影響程度。

1 數(shù)值方法

1.1 幾何模型

DTU 10 MW風(fēng)力機(jī)葉片的氣動外形使用FFA-W3翼型族設(shè)計，該翼型族由瑞典航空研究院設(shè)計，包含7種翼型，厚度范圍為19.5%～36%，是非常有代表性的風(fēng)力機(jī)專用翼型。DTU 10 MW風(fēng)力機(jī)風(fēng)輪直徑達(dá)到178.3 m，單只葉片長86.466 m。在復(fù)雜的風(fēng)況下，葉片柔性大，導(dǎo)致葉尖位移較大，有很大可能面臨顫振風(fēng)險的可能。FFA-W3-241翼型的相對厚度小，升阻比高，被用于DTU 10 MW風(fēng)力機(jī)葉尖位置的氣動外形設(shè)計。因此本文選用FFA-W3-241翼型進(jìn)行研究，其幾何外形如圖1所示。

圖1 FFA-W3-241翼型幾何外形圖

1.2 計算域與網(wǎng)格

本文中CFD模擬采用嵌套網(wǎng)格方法。嵌套網(wǎng)格由兩套網(wǎng)格構(gòu)成，分別是背景網(wǎng)格和貼體網(wǎng)格。背景網(wǎng)格保證遠(yuǎn)場的計算，貼體網(wǎng)格隨翼型共同運(yùn)動，確保翼型周圍的網(wǎng)格質(zhì)量。背景網(wǎng)格的計算域為正方形，各邊長度為40倍弦長(40c);翼型旋轉(zhuǎn)軸位于正方形中心，各邊距離旋轉(zhuǎn)軸20c；貼體網(wǎng)格所在的區(qū)域為距離翼型旋轉(zhuǎn)軸3c的圓形。

為保證翼型壁面的計算精度，貼體網(wǎng)格包含了15層邊界層網(wǎng)格，膨脹率為1.2；邊界層外側(cè)膨脹率為1.5，并保證壁面網(wǎng)格y+<1；翼型上下表面各200個節(jié)點(diǎn)，F(xiàn)FA-W3-241翼型為鈍尾緣，尾緣設(shè)置50個節(jié)點(diǎn)。背景網(wǎng)格的各邊長的網(wǎng)格節(jié)點(diǎn)數(shù)為150，等距分布。貼體網(wǎng)格總數(shù)為50 958，背景網(wǎng)格總數(shù)為22 201，背景網(wǎng)格和貼體網(wǎng)格如圖2所示。

圖2 計算域網(wǎng)格

1.3 邊界條件

根據(jù)風(fēng)洞試驗條件[20]，CFD數(shù)值模擬采用的雷諾數(shù)為160萬，湍流強(qiáng)度為1.0%，左側(cè)邊界為速度入口，翼型表面設(shè)置為光滑的無滑移壁面，上下邊界為了實現(xiàn)平行流動，速度和壓力的邊界條件設(shè)置為零壓力梯度。計算域和邊界條件如圖3所示。求解器采用開源的OpenFOAM求解器。湍流模型采用k-ω SST模型。

圖3 計算域及邊界條件

2 數(shù)值方法確認(rèn)

2.1 動態(tài)時間步長、計算周期驗證

對于非定常計算，時間步長對計算結(jié)果和計算成本影響較大。較小的時間步長能精確的描述物理現(xiàn)象，但是增加了計算成本。因此，為了得到精確的物理現(xiàn)象和較低的計算成本，進(jìn)行了俯仰振蕩的時間步長驗證。俯仰振蕩的攻角變化為α=αmea+αamp·sin(ωt)。以平均攻角αmea為15.6°，振蕩幅值αamp為1.8°為例，在一個振蕩周期內(nèi)設(shè)計了4種時間步長，分別是一個周期內(nèi)計算360、720、1 080和1 440步。俯仰振蕩的折合頻率為0.093，折合頻率定義為

(1)

式中:f為振蕩頻率;c為弦長;U∞為自由來流速度;ω為角速度。

圖4展示了第1～5個周期內(nèi)各時間步長的升力系數(shù)時間歷程，其中第1個周期內(nèi)的結(jié)果未收斂，未在圖中展示。一個周期內(nèi)，采用360和720步的結(jié)果與采用1 080和1 440步的結(jié)果相差較大1 080和1 440步的結(jié)果幾乎重合，達(dá)到收斂。同時可以看出，采用720步長時，3個周期后的結(jié)果達(dá)到收斂。

圖4 升力系數(shù)的收斂歷史

2.2 數(shù)值驗證

圖5展示了未失速、輕度失速和深度失速三種情況下的升力系數(shù)曲線，并與文獻(xiàn)[20]中的實驗結(jié)果進(jìn)行對比。可以看出，穩(wěn)態(tài)計算時，附著流區(qū)域計算值與實驗吻合較好，然而在失速區(qū)內(nèi)，模擬結(jié)果對分離點(diǎn)的預(yù)估相比于實驗稍有推遲，升力系數(shù)計算結(jié)果偏高于實驗值。這是由于RANS對于分離渦的捕捉不夠精準(zhǔn)。在瞬態(tài)計算中，附著流區(qū)穩(wěn)態(tài)增長斜率較大，遲滯回環(huán)的斜率也較大。失速區(qū)內(nèi)，靜態(tài)升力系數(shù)曲線穿過動態(tài)遲滯回環(huán)，攻角15°時穩(wěn)態(tài)計算值與實驗值偏差相對較大，20°后差距逐漸減小，因此瞬態(tài)計算中，15°附近的遲滯回環(huán)計算偏大，20°附近差距減小。動態(tài)計算的誤差在可以接受的范圍內(nèi)。

圖5 三種工況下的升力系數(shù)曲線

3 兩自由度耦合振蕩結(jié)果分析

翼型的俯仰振蕩和沉浮振蕩分別對應(yīng)了葉片的揮舞與扭振，運(yùn)動方程如下式。

α=αmea+αamp·sin(ωt)

(2)

h=H·cos(ωt)

(3)

由于翼型運(yùn)動導(dǎo)致的弦線方向時刻改變，與來流方向夾角也時刻改變，因此存在實際等效攻角。根據(jù)文獻(xiàn)[21]，俯仰沉浮耦合振蕩時流場的有效攻角為式(4):

(4)

3.1 輕度失速顫振邊界

圖6 不同沉浮幅值下凈能量傳遞系數(shù)隨俯仰、沉浮折合頻率的變化(αmea=15°，αamp=4°)

圖7 凈能量傳遞系數(shù)等值線云圖(αmea=15°，αamp=4°)

3.2 深度失速顫振邊界

圖8 不同沉浮幅值下凈能量傳遞系數(shù)隨俯仰、沉浮折合頻率的變化(αmea=18°，αamp=4°)

圖9 凈能量傳遞系數(shù)等值線云圖(αmea=18°，αamp=4°)

3.3 輕度失速的不確定性分析

采用能量法分析顫振的穩(wěn)定性，由于能量法沒有考慮流固互相作用的過程，預(yù)先給定了振動的振幅和振型。而在真實工作環(huán)境下，存在來流風(fēng)速、幾何誤差等不確定性因素，使得失速顫振時振幅和振型不斷變化。由于3.1節(jié)中輕度失速的顫振邊界不明確，綜合不確定性因素影響，本節(jié)將振幅和振型作為不確定性量，對輕度失速的顫振發(fā)生進(jìn)行不確定性CFD數(shù)值模擬與分析。深度失速由于顫振邊界明確，將不再作不確定性討論。

考慮了俯仰平均攻角、俯仰攻角振幅、俯仰折合頻率、沉浮距離振幅、沉浮折合頻率五個不確定因素。假設(shè)這些不確定變量服從高斯分布，均值及標(biāo)準(zhǔn)差見表1。采用稀疏網(wǎng)格法構(gòu)造多項式混沌，對應(yīng)61個不確定工況?；谙∈杈W(wǎng)格的多項式混沌方法參見文獻(xiàn)[18]。

表1 不確定變量均值及標(biāo)準(zhǔn)差

敏感度分析采用的是文獻(xiàn)[22]的整體敏感性指數(shù)。圖11計算了5種不確定因素對凈能量傳遞系數(shù)的靈敏度指數(shù)，發(fā)現(xiàn)影響最大的因素是俯仰頻率，而不是本征值較大的攻角，其次影響較大的是沉浮頻率，沉浮距離、攻角均值及振幅對顫振的敏感度不高。

圖11 敏感性分析

4 結(jié) 論

本文采用CFD模擬對FFA-W3-214翼型進(jìn)行了俯仰與沉浮耦合運(yùn)動的數(shù)值計算，并用能量法找尋了輕度失速和深度失速工況下的顫振邊界，得出如下結(jié)論：

(1)對于輕度失速，當(dāng)沉浮距離幅值為一倍弦長時，在低折合頻率下有顫振發(fā)生，但邊界不明顯。當(dāng)沉浮距離增大至兩倍弦長，隨折合頻率的增加，顫振現(xiàn)象更加明顯；

(2)在深度失速工況下，顫振現(xiàn)象較輕度失速更明確，可找出明確的顫振邊界。且沉浮振幅增大使顫振邊界中俯仰折合頻率范圍減小，沉浮折合頻率范圍無明顯影響；

(3)通過對顫振邊界不明確的輕度失速工況運(yùn)動參數(shù)的不確定量化分析，可以預(yù)測該工況下顫振發(fā)生的概率達(dá)49.48%。此外，運(yùn)動參數(shù)中對顫振影響最大的因素是俯仰折合頻率，其次是沉浮折合頻率，沉浮距離、攻角均值及振幅對顫振的敏感度不高。