自然風速下低雷諾數風力機翼型氣動特性隨機數值分析

唐新姿，袁可人，王效禹，彭銳濤

（湘潭大學機械工程學院，湖南湘潭 411105）

翼型氣動特性是風力機設計決定性因素，其特性參數通常在確定性條件下獲得.低風型風力機工作雷諾數較低，常年工作在地表層自然來流環境中，氣流極不穩定，風速大小和風向均存在隨機不確定特點，其翼型氣動特性受到不確定風速影響，給風力機實際運行性能與載荷預測帶來困難.研究隨機自然來流對低雷諾數風力機翼型氣動性能的影響，對于風電穩定可靠發展具有重要意義.

對于風力機翼型氣動特性研究，目前國內外主要基于確定性工況條件進行風洞試驗或數值模擬.黃宸武等[1]開展了S809 翼型低雷諾數下的氣動特性風洞試驗研究.Wang W C 等[2]通過風洞測試研究了有無湍流條件下風力機靜動態性能.李仁年等[3]采用數值計算方法分析了翼型后緣厚度對翼型升阻力系數的影響.姜鑫等[4]采用數值模擬計算了多種工況不同曲率葉片翼型升力系數.為考慮非定常流動，朱志斌等[5]采用大渦模擬方法計算低雷諾數下翼型升阻力系數.如果直接在數值模擬中加入不確定邊界條件計算耗時耗力，傳統確定性氣動數值分析已經不能滿足實際工程需求.

為考慮不確定因素，國內外學者將傳統確定性氣動數值分析和不確定理論相結合，發展了不確定數值分析方法解決氣動參數隨機問題.如Abdallah I等[6]將數值模擬與蒙特卡洛法相結合量化翼型升阻力系數不確定性.該方法依賴于抽樣數量，對于復雜流動分析計算效率不理想，計算時間較長，難以實際應用.為解決計算效率問題，董世充等[7]、琚亞平等[8]和劉智益等[9]基于譜展開法、嵌入式和非嵌入式混沌多項式方法，分別研究隨機表面粗糙度對風力機翼型氣動特性影響和隨機失諧葉片安裝角對風力機流場及氣動特性影響.趙軻等[10]基于混沌多項式展開并結合數值模擬方法分析了馬赫數對翼型氣動性能影響.鄔曉敬等[11]采用非嵌入式混沌多項式方法分析了飛行狀態不確定性對翼型氣動力影響.Zhong X P 等[12]與丁繼鋒等[13]采用多項式響應面模型進行不確定馬赫數條件下翼型穩健設計.

風力機隨機氣動問題正受到越來越多研究者的重視，研究發現自然來流風速不確定會導致風力機理論設計與實際運行產生偏差[14]，但由于自然隨機風速很難通過實驗方式產生，風速大小和方向不確定雙因素耦合作用對風力機翼型氣動特性影響及其不確定性在流場中的傳播機制亟待明確.

本文基于多項式混沌法和隨機配點法以及Transition SST 湍流方程，構建隨機翼型氣動數值分析模型，以S809 風力機翼型為研究對象，從不確定風速大小和方向均值、方差與氣動力等方面討論隨機自然來流對翼型升力和阻力不確定性影響，量化隨機風速大小和方向工況下風力機翼型氣動特性的不確定程度，為極限載荷和疲勞載荷分布譜分析和可靠性評估提供依據；獲得隨機風速條件下風力機翼型氣動力確定性和不確定性成分比例，通過壓力系數摩阻系數和湍動能分布揭示不確定性在流場中的傳播機制，為風電葉片氣動穩健性研究提供參考.

1 研究方法

1.1 低雷諾數數值分析方法及驗證

在低雷諾數流動下，準確判斷層流分離與轉捩位置至關重要.經過湍流模型對比分析，本文采用四方程轉捩模型Transition SST[15]預測層流分離.計算雷諾數為3×105，計算模型S809 翼型弦長c 為1 m，對應自然風速為4.38 m/s，采用二維“C”型結構網格，翼型上下及來流上游方向取15 倍翼型弦長，翼型下游方向取20 倍弦長作為計算域邊界.采用ICEM 生成結構化網格，近壁面進行局部加密，物面法向第一層網格高度為10-5倍弦長，保證壁面處y+＜1，網格增長率為1.1，近翼面網格如圖1 所示.流域左邊界設為速度入口，右邊界設為壓力出口，壓力設置為絕對壓力101 325 Pa，翼面為無滑移條件.采用二階迎風差分格式進行離散，壓力速度耦合采用SIMPLE 算法，收斂精度標準均為殘差小于1×10-5.

圖1 S809 翼型近翼面網格Fig.1 Local grid of S809

為驗證網格無關性，計算了風速4.38 m/s、風向6.11°時4 種不同網格數量的升阻力系數，計算結果對比如表1 所示.由表1 可知，網格數量達到89 700時，增加網格數量計算結果不再明顯變化，綜合考慮計算效率，最終確定網格總數為89 700.

表1 網格無關性驗證Tab.1 Grid independence verification

為進一步驗證計算方法的正確性，分別計算不同攻角下S809 翼型升阻力系數并與實驗值[16]進行對比，計算結果如圖2 所示.結果表明采用上述方法和網格計算所得的升阻力系數計算值與實驗值在低攻角下基本一致，驗證了該方法的可行性.

圖2 S809 翼型升阻力系數實驗值與計算值比較Fig.2 Comparison of calculated and experimental lift and drag coefficients of S809 airfoil

1.2 隨機分析方法與驗證

本文將多項式混沌法和隨機配點法結合，基于Transition SST 湍流方程，形成隨機翼型氣動分析方法.多項式混沌法和隨機配點法聯合又稱為非嵌入式概率配置點法（Non-intrusive Probabilistic Collocation Method，NIPRCM），參照文獻[17]，得出任意樣本空間上的變量關系如式（1）所示，其中yk（x，t）是由空間自變量x，時間t 和參數θ 組成的y 在第k 個配置點的值，Np為配置點的個數，hk（ξ）為對應的拉格朗日插值多項式混沌.

根據隨機變量分布函數，將選取的高斯積分點及對應的求積系數映射到隨機空間P 上，即得到配置點及其權重.得到權重系數后，隨機變量y 的均值y 和方差σ2可以通過式（2）和（3）計算，其中yk為物理參數的確定性計算結果，wk為權重系數.

根據文獻[18]可知，用公式（2）和（3）分別采用二階與三階方法計算翼型氣動參數所得計算結果基本一致，說明二階方法計算精度已滿足要求，因此后續采用二階NIPRCM 方法進行隨機分析.

為進一步驗證所采用的隨機分析方法，分別計算在風速服從高斯分布的不確定條件下S809 翼型升阻力系數并與實驗值[16]進行對比，計算雷諾數為3×105，攻角分別為1.99°、4.08°、6.11°、8.14°，標準差為10%均值，計算結果如圖3 所示.采用上述隨機方法計算所得的升阻力系數均值與確定性計算值基本一致，且兩者與實驗值吻合，驗證了該方法的可行性.

圖3 S809 翼型隨機分析方法驗證Fig.3 Validation of uncertain calculation of S809 airfoil

2 隨機風速對翼型氣動特性的影響

為量化風速大小隨機性對風力機翼型氣動性能的影響，計算雷諾數為3 × 105，典型攻角1.99°、6.11°、8.14°時風速大小標準差分別為均值的5%、10%、15%，即湍流強度為5%、10%、15%三種工況下S809 翼型升阻氣動特性.

2.1 升阻特性

圖4 為三種隨機風速條件下攻角為1.99°、6.11°和8.14°時S809 翼型的升阻力正態分布，經無關性分析確定隨機計算總頻次為5 000.標準差越大，即風速隨機性越大，升阻比不確定程度增加，但升阻比均值基本保持不變，且均值對應頻次減少.攻角8.14°時，隨機性對升阻比影響最大.由圖4（c）可知，風速大小標準差σv為15%時，翼型升阻比標準差分別是標準差為5%和10%時升阻比標準差的4.15 倍和1.75 倍，升阻比3σ 不確定度為39.576 7±11.947 5.

圖4 隨機風速下S809 翼型升阻比Fig.4 Lift drag ratio of S809 airfoil with stochastic wind speed

2.2 壓力分布

圖5 給出了攻角1.99°時隨機風速大小條件下S809 翼型上翼面壓力系數Cp的標準差分布，圖中橫坐標x/c 表示翼型表面點的相對弦長位置，即翼型表面點在弦長方向的投影坐標x 與弦長c 的比值.風速標準差越大，壓力系數標準差峰值越大；風速大小標準差由5%增加到15%，壓力系數標準差最大值由0.010 2 增加到0.047，且各標準差峰值位置始終出現在翼型0.6 倍弦長附近位置，這是由于在此位置翼型邊界層出現層流分離泡[15]，流動分離轉捩，對風速大小波動更為敏感.

圖5 隨機風速下攻角1.99°時S809 翼型上翼面壓力系數標準差Fig.5 Standard deviation of pressure coefficient of upper surface of S809 airfoil at attack of angle 1.99°under stochastic wind speed

2.3 摩阻分布

圖6 是攻角為1.99°時隨機風速條件下S809 翼型上翼面摩阻系數Cf的標準差分布.風速大小標準差由5%增加到15%，摩阻系數標準差最大值由5.58×10-4增加到1.66×10-3.翼型前緣出現極大值，且在0.6 倍弦長附近位置再次出現峰值，規律與前述壓力分布一致.

2.4 湍動能分布

圖7 是攻角為1.99°時隨機風速條件下S809 翼型上翼面湍動能k 的標準差分布.由圖7 可知，風速大小標準差由5%增加到15%，湍動能最大標準差由1.21×10-4增加到2.91×10-4.湍動能標準差峰值始終位于翼型0.6 倍弦長附近位置，與壓力和摩阻系數分布規律一致.

圖6 隨機風速下攻角1.99°時S809 翼型上翼面摩阻系數標準差Fig.6 Standard deviation of friction coefficient of upper surface of S809 airfoil at attack of angle 1.99°under stochastic wind speed

圖7 隨機風速下攻角1.99°時S809 翼型上翼面湍動能標準差Fig.7 Standard deviation of turbulent kinetic energy of upper surface of S809 airfoil at attack of angle 1.99°under stochastic wind speed

3 隨機風向對翼型氣動特性的影響

為分析風速方向隨機性對風力機翼型氣動性能的影響，在雷諾數為3×105，攻角1.99°、6.11°、8.14°工況，分別取標準差為均值的5%、10%、15%這三種風向分布，計算風速方向不確定S809 翼型氣動性能.

3.1 升阻特性

圖8 為三種隨機風向條件下攻角為1.99°、6.11°和8.14°時S809 翼型升阻比正態分布，計算總頻數為5 000.風向標準差越大，即風向隨機性越大，升阻比不確定程度增加.攻角8.14°時，風向隨機性對升阻比影響最大，風向標準差為15%時，翼型升阻比標準差σ 分別是5%和10%時的1.95 倍和1.18 倍，其升阻比3σ 不確定度為39.891 9±24.742 5.

圖8 隨機風向下S809 翼型升阻比Fig.8 Lift drag ratio of S809 airfoil with uncertain wind direction

3.2 壓力分布

圖9 給出了攻角1.99°時隨機風向條件下S809翼型上翼面壓力系數標準差分布.由圖9 可知，風向標準差由5%增加到15%，壓力系數最大標準差由0.032 3 增加到0.100 1.壓力系數最大標準差位于前緣，翼型0.6 倍弦長附近位置出現小峰值.

圖9 隨機風向下攻角1.99°時S809 翼型上翼面壓力系數標準差Fig.9 Standard deviation of pressure coefficient of upper surface of S809 airfoil at attack of angle 1.99°under stochastic wind direction

3.3 摩阻分布

圖10 是攻角為1.99°時隨機風向條件下S809翼型上翼面摩阻系數標準差分布.風向標準差由5%增加到15%，摩阻系數最大標準差由9.29×10-4增加到2.92×10-3，最大標準差位置位于前緣，且在翼型0.6 倍弦長附近位置再次出現峰值.

圖10 隨機風向下攻角1.99°時S809 翼型摩阻系數上翼面標準差Fig.10 Standard deviation of friction coefficient of upper surface of S809 airfoil at attack of angle 1.99°under stochastic wind direction

3.4 湍動能分布

圖11 是攻角為1.99°時隨機風向條件下S809翼型上翼面湍動能標準差分布.由圖11 可知，風向標準差由5%增加到15%，湍動能最大標準差由7.8 × 10-6增加到5.33 × 10-5，其峰值位置位于0.6倍弦長附近位置，敏感位置與壓力摩阻系數分析結論一致.

圖11 隨機風向下攻角1.99°時S809 翼型上翼面湍動能標準差Fig.11 Standard deviation of turbulent kinetic energy of upper surface of S809 airfoil at attack of angle 1.99°under stochastic wind direction

4 隨機風速風向耦合對翼型氣動特性的影響

為分析風速和風向同時存在隨機性對風力機翼型氣動性能的影響，在雷諾數為3×105，攻角1.99°、6.11°、8.14°工況，各計算標準差為5%均值，計算S809 翼型氣動性能.

4.1 升阻特性

圖12 給出了風速、風向單隨機因素以及風速和風向同時存在隨機性，攻角分別為1.99°、6.11°和8.14°時S809 翼型升阻比正態分布，計算總頻數為5 000.相同攻角下，隨機耦合作用下的升阻比標準差大于單隨機因素下的升阻比標準差.攻角8.14°時，風向不確定性對升阻比影響最大.由圖12（c）可知，風速風向隨機耦合作用下，翼型升阻比3σ 不確定度為39.022 3±13.708 5，即3σ 置信區間相對不確定度為±35.13%，且其標準差σv，α分別是風速（σv）和風向（σα）單隨機因素時的4.76 倍和1.08 倍.

圖12 隨機風速風向耦合下S809 翼型升阻比Fig.12 Lift drag ratio of S809 airfoil with uncertain wind speed and uncertain direction

4.2 壓力分布

圖13 為隨機風速、風向以及隨機風速和風向耦合作用下，攻角分別為1.99°、6.11°和8.14°時S809翼型上翼面壓力系數標準差分布.隨機風速風向耦合作用時，壓力系數不確定程度大于單隨機因素影響程度；且攻角增大，壓力系數不確定度增加.在隨機風速和風向耦合作用下，攻角1.99°、6.11°和8.14°時壓力系數最大標準差分別為0.035 2、0.214 1 和0.955 6，各壓力系數標準差最大值均出現在前緣.

圖13 隨機風速風向耦合下S809 翼型壓力系數標準差Fig.13 Standard deviation pressure coefficient of S809 airfoil with uncertain wind speed and uncertain direction

5 結論

基于NIPRCM 方法和Transition SST 轉捩模型，建立低雷諾數翼型氣動特性隨機數值分析模型，量化隨機風速大小和方向對風力機翼型升阻力和流場特性的影響，獲得自然來流條件下風力機翼型氣動特性不確定程度及流場非穩傳播規律.主要結論如下：

1）隨機風速風向對翼型升阻氣動因子不確定度影響顯著，且來流隨機性越大升阻比不確定程度越大，在計算攻角范圍內S809 翼型升阻比3σ 置信區間相對不確定度最大為±35.13%，該值遠大于常規設計標準靜態氣動因子評估不確定度，表明風力機設計應根據實際安裝位置氣流狀況，加強動態時變極限載荷和疲勞載荷風險評估.

2）隨機風速風向耦合作用時翼型升阻力系數不確定度比隨機風向或隨機風速時的不確定度明顯增大，且風向不確定影響更大.攻角8.14°隨機風速風向耦合作用且標準差為5%時，S809 翼型升阻比不確定度分別是風速和風向單隨機因素下的4.76 倍和1.08 倍.

3）在隨機自然風速條件下，翼型對來流不確定性敏感區域為前緣，且中部0.6 倍弦長位置轉捩處出現非穩流動加劇，可以考慮在翼型前緣部分進行氣動穩健性優化設計.