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

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

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

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

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

1 風電葉片數(shù)值計算方法

1.1 湍流模型修正

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

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

輸運方程中的σk和σω對風力發(fā)電機葉片翼型氣動力有影響[5]。通過對σk和σω進行修正，確定風力發(fā)電葉片周圍空氣湍流粘度。修正通過混合函數(shù)F1，F(xiàn)1將系數(shù)的近壁面貢獻系數(shù)與遠壁面結(jié)合在一起，表達式如下：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模型參數(shù)進行修正，確保風力發(fā)電機葉片翼型氣動評估模型滿足氣動設(shè)計要求。

1.2 風力發(fā)電機組入流風剪切模型

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

圖1 入流剪切風速定義

1.3 風力發(fā)電機組氣動分析模型

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

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

圖2 邊界條件定義

1.4 風力發(fā)電機組氣動分析模型網(wǎng)格劃分

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

圖3 葉片表面網(wǎng)格及周圍網(wǎng)格

2 計算結(jié)果分析

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

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

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

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

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

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

表1 軸功率計算結(jié)果對比

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

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

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