基于CFD風力發電機旋轉風輪氣動性能數值模擬

株洲時代新材料科技股份有限公司 李東旭 梁鵬程 鄧 航 譚 龍 易禮毅

作為可再生能源，風能在能源中所占比重不斷增加，預計2020年至2025年期間全球海上風能投資總額預計將達到2110億美元。隨著風電機組單機容量的不斷擴大，風力發電葉片的長度也發展到70m 以上。大型風力發電機旋轉風輪氣動性能評估成為行業主要研究課題。

風力發電機旋轉風輪氣動性能評估方法主要包括葉素動量理論、CFD 數值模擬、風洞試驗、風場現場測試。近年來，隨著計算流體力學和計算機硬件的不斷發展，數值模擬技術已逐漸成為認識風力發電機葉片翼型流動規律的重要研究手段。國內主要研究機構借助數值模擬和風洞測試方法對風力機翼型氣動性能做了大量研究，相關研究成果已廣泛應用于風電葉片氣動設計。

由于二維翼型氣動性能評估方法和風電葉片Bladed 氣動載荷計算軟件均無法準確預測風輪三維旋轉效應對葉片氣動性能的影響，因此國內外開始采用CFD 方法研究風力發電機風輪旋轉流場分析，Thanhtoan Tran 等運用計算流體力學對海上浮動風力發電機葉片氣動性能進行了研究，考慮海洋平臺周期性運動對氣動性能的影響[1]。李少華等借助Fluent 對風力發電機旋轉風輪氣動性能進行了評估，研究了旋轉風輪氣動性能影響因素[2-4]。隨著風電葉片大型化發展趨勢，迫切需要建立更加準確的風力發電機旋轉風輪氣動性能評估模型，考慮風剪切、塔筒、機艙、湍流模型、葉片附近區域湍流信息等對風力發電機葉片氣動性能影響。

為準確對風力發電機旋轉風輪氣動性能評估，本文借助流體動力學軟件STAR-CCM+軟件對公司新開發葉片開展三維旋轉流體動力學數值模擬研究，建立風力發電機風輪幾何模型，包括塔筒、機艙、輪轂、葉片等。建立風力發電機來流風剪切模型，對葉片表面及附近區域網格進行細化，并對SST k-ω 湍流模型進行修正，通過計算不同來流風速下風力發電機旋轉風輪氣動力，與Bladed 軟件氣動分析結果進行分析對比，本文建立的風力發電機風輪氣動分析模型計算精度控制在5%以內，為后續大型風力發電機葉片氣動性能評估提供技術支持。

1 風電葉片數值計算方法

1.1 湍流模型修正

風力發電機組氣動性能分析模擬運用剪切應力輸運(Shear-Stress Transport)k-ω 模型。SST kω 模型為了使標準k-ω 模型在近壁面區有更好的精度和算法穩定性，由Menter 對標準k-ω 模型進行改進而來。SST k-ω 湍流模型的方法是：在近壁面使用標準k-ω 模型，而在邊界層外部使用k-ε，在邊界層內則混合使用這兩種模型，并根據一個混合加權函數F_1的大小進行加權平均。

SST k-ω 湍流模型運用兩方程求解輸運方程，通過湍流動能k 和耗散率ω 來確定湍流粘度。湍動能k和耗散率ω 輸運方程如下：其中，為平均速度、m/s，μ為動力粘度、Pa·s，σk、σw為模型系數，Pk、Pw為產生項，fβ＊為自由剪切修改因子，fβ為漩渦拉伸修改因子，Sk、Sw為用戶源項，k0、ω0為周圍湍流值，用于抵消湍流衰減。

輸運方程中的σk和σω對風力發電機葉片翼型氣動力有影響[5]。通過對σk和σω進行修正，確定風力發電葉片周圍空氣湍流粘度。修正通過混合函數F1，F1將系數的近壁面貢獻系數與遠壁面結合在一起，表達式如下：500v/(d2ω))，2k/(d2CDkω)))4，式中，d 為到最近壁面的距離、m，v 為運動粘度、m2/s，CDkω=max(1/ω ▽k ▽ω，10-20)。

σk和σω修正公式如下：σk=F1σk1+(1-F1)σk2，σω=F1σω1+(1-F1)σω2，運用流體動力學軟件STAR－CCM+軟件對SST k-ω 湍流模型中σk1、σk2、σω1、σω2模型參數進行修正，確保風力發電機葉片翼型氣動評估模型滿足氣動設計要求。

1.2 風力發電機組入流風剪切模型

風力發電機風輪來流風速隨高度成指數變化[6]，高度為Z 處來流平均風速為：V（z）=Vhub（z/zhub）α，其中，V（z）為高度z 處的風速、m/s，z 為高度、m，zhub為輪轂高度、m，α 為冪指數（冪指數取0.2）。運用STAR－CCM+軟件Field Fuctions 功能定義入流風速，剪切風速分布如圖1。

圖1 入流剪切風速定義

1.3 風力發電機組氣動分析模型

建立3MW 風力發電機組幾何模型，風力發電機風輪直徑為121m、塔筒高度80米。運用流體動力學軟件STAR－CCM+建立3MW 風力發電機組流體動力學分析模型，模型包括風輪旋轉域和空氣靜止域兩部分，旋轉域包括葉片和輪轂，靜止域包括塔筒、機艙等。

邊界條件定義如圖2所示。入口采用速度入口邊界條件，考慮風速切變的影響。出口采用0壓力出口邊界條件，即設定出口壓力為大氣壓。為保證計算域入口和出口位置對風力發電機氣動性能沒有影響，入口距離風力發電機距離為5倍風輪直徑，出口距離風力發電機距離為10倍風輪直徑。計算域左右兩側距離風力發電機為5倍風輪直徑，并且定義對稱邊界條件。風輪葉片表面設置為無滑移壁面邊界條件，定義風輪繞旋轉軸以恒定轉速轉動。

圖2 邊界條件定義

1.4 風力發電機組氣動分析模型網格劃分

為有效控制風力發電機氣動分析模型網格數量，采用STAR－CCM+軟件切割體單元對計算域進行網格劃分，由于風電葉片最大弦長到葉尖區域對風電葉片氣動升力貢獻量在70%以上，為準確捕捉葉片該區域流場信息，采用STAR－CCM+軟件網格加密技術對旋轉域葉片最大弦長到葉尖區域網格進行局部加密，葉片表面網格尺寸由葉尖到最大弦長依次增大，旋轉域體網格尺寸由葉片往外逐漸增大。葉尖表面網格目標尺寸為0.016m，葉片其他區域表面尺寸為0.032m。葉片葉尖附近體網格尺寸為0.032m，葉片其他區域附近體網格目標尺寸為0.064m。旋轉域網格數為17262750，空氣靜止域網格數為6995449。葉片表面及周圍網格如圖3。

圖3 葉片表面網格及周圍網格

2 計算結果分析

分別運用旋轉坐標系方法和滑移網格方法，完成不同來流風速下3MW 風力發電機旋轉風輪氣動力計算，計算輸出風輪在恒定轉速運動下繞旋轉軸的扭矩，并轉換成風力發電機功率（繞旋轉軸機械功率），轉換公式如下：P=T×n×2π/60，式中：P 為3MW風力發電機功率、W，T 為風輪在恒定轉速運動下繞旋轉軸的扭矩、N·m，n 為葉輪轉速、r/min。

圖4為旋轉坐標系方法計算輸出的3MW 風力發電機風輪軸功率，由此可知：隨著來流風速提高，風輪軸功率逐漸增加。旋轉坐標系方法模擬3MW風力發電機風輪在特定位置所受到的氣動載荷，屬于穩態計算方法，經過5000迭代計算，軸功率計算結果基本穩定。

圖4 不同風速下風輪機械軸功率（旋轉坐標系方法）

圖5為滑移網格方法計算輸出的3MW 風力發電機風輪軸功率，由此可知：由于滑移網格方法為非穩態計算方法，風輪以恒定轉速旋轉過程中，由于存在風切變、動態旋轉和塔筒效應，葉片受到的氣動力會周期性變化，葉片旋轉到最高點附近區域所受到的氣動力最大，葉片旋轉到最下端附近區域時受到的氣動力最小。運用滑移網格方法計算時，計算迭代到第三周時，風輪旋轉過程中葉片受到的氣動力基本保持周期性變化。

圖5 不同風速下風輪機械軸功率（滑移網格方法）

為對比不同來流風速下氣動力計算數據，旋轉坐標系方法最后1000步氣動力結果取平均，滑移網格方法最后一周（第四周）氣動力取平均。由表1可知，采用滑移網格方法計算的風力機軸功率與Bladed 軟件計算的軸功率誤差控制在5%以內。采用旋轉坐標系方法計算的風力機軸功率與Bladed軟件計算的軸功率誤差控制在8%以內。

表1 軸功率計算結果對比

運用本文建立的風力發電機風輪氣動性能模型，評估風力發電機風輪旋轉過程中葉片葉根到最大弦長區域氣動性能，運用STAR-CCM+軟件制作該區域速度流線圖可知：5m/s 風速時由于葉根到最大弦長區域翼型厚度較大，同時風輪以恒定轉速旋轉時會產生離心力，導致葉根到最大弦長區域氣流會向葉尖方向偏離，因此吸力面后緣會產生失速現象，導致該區域氣動性能下降；9m/s 風速時由于風輪轉速提高，離心力會相應增大，導致葉根到最大弦長區域氣流會向葉尖方向偏離區域面積增加。風力發電機風輪氣動性能模型后續可對葉根到最大弦長區域進行氣動性能優化設計，通過添加渦流發生器來抑制氣動分離，進一步提高風力發電機風輪氣動效率。此外，該模型可輸出葉片表面均布載荷，用于葉片整體結構校核和局部結構優化設計。

綜上，本文運用CFD 方法對風力發電機旋轉風輪氣動性能進行仿真計算。研究發現：旋轉坐標系方法和滑移網格方法均能滿足風力發電機旋轉風輪氣動性能計算要求。通過對SST k-ω 湍流模型中的模型系數σk和σω進行修正，同時通過定義風剪切模型、葉片表面及周圍區域網格細化，可有效提高風輪旋轉過程中氣動力計算精度，滑移網格方法計算結果與Bladed 軟件精度控制在5%以內。

近年來，葉片氣動性能提升技術得到越來越多風力發電主機廠與業主的關注，通過建立準確的風力發電機組風輪三維旋轉氣動力計算模型，可準確評估渦流發生器、葉尖小翼等氣動性能提升組件的氣動性能提升效果，同時，風力發電機組風輪三維旋轉氣動力計算模型可為大型風電葉片氣動彈性設計提供準確載荷數據，為風電葉片大型化發展提供技術支撐。