IDDES 方法模擬風切變下大型水平軸風力機流場特性

楊從新，張宇婷，岳念西

（1.蘭州理工大學能源與動力工程學院，甘肅 蘭州 730050；2.甘肅省流體機械及系統重點實驗室,甘肅 蘭州 730050）

大氣邊界層內的氣流受地形、地表植物、建筑等障礙物的影響，流速隨高度而發生變化，越遠離地面流速越大且到達某一高度后保持相對穩定。隨著風電機組大型化，輪轂高度和風輪直徑越來越大，這種不均勻的剪切流對風力機的影響越來越明顯，導致風力發電機組實際輸出功率難以達到預期值，而且在風力機后方形成較大的速度損失和復雜的渦結構等；因此，風力機氣動性能分析和流場分析很有必要考慮風剪切的影響。目前，已有不少國內外學者用不同方法分析了風力機的氣動性能和流場狀況。文獻[1-5]以小型風力機為研究對象，利用風洞實驗測量數據，建立了較詳細的風力機氣動信息數據庫，為分析風力機的流場特征、驗證數值模擬方法提供了可靠的依據。文獻[6]用BEM 方法和RANS(雷諾平均模擬)方法對飛行器專用風機在高空環境中的氣動特性進行了瞬態模擬比較，發現在設計工況下兩種方法有較好的一致性。文獻[7-8]基于致動線模型結合RANS 方法和LES(大渦模擬)方法模擬小型風力機尾流中的湍流結構，結果發現均勻入流與剪切入流條件下風力機流場的主要區別是尾流場速度虧損的對稱性，另外二者在距地面較高處雷諾應力表現出相似性。文獻[9]基于致動面模型結合N-S 方程分析了單臺Nibe A 型風力機尾流區域內的風速變化、湍流強度、渦結構等。文獻[10-11]采用RANS 方法對風力機進行全尺度數值模擬，分析了風切變系數對風力機近尾流區域特性的影響，以及不同葉尖速比下風力機載荷分布和機艙、塔架附近流場結構。以上文獻調研可以看出風力機流場分析方法的多樣化，其中實驗方法是研究風力機空氣動力學問題的基礎，能比較全面地獲取流場信息，促進人類對水平軸風力機流場狀態的認知，但其耗時耗力，實現難度較大。致動線、致動面模型結合CFD 方法將風力機簡化為線或面，氣流經過風力機時無法精確描述葉片表面的流動狀態，對流場細節的描述還有待改進。數值模擬中的RANS 方法對定常工況描述具有優勢，但風力機運行大多處于非定常狀態，顯然該方法的計算精度略低。LES 方法計算精度足夠高，但它對網格密度要求高、計算量大，已有的研究大多是針對小型風力機。針對以上問題，有人構建RANS-LES 混合模型用于瞬態流動分析，從原始的DES(分離渦模擬)方法到DDES(延遲分離渦模擬)方法再到Spalart 等[12]提出的IDDES 方法，一步步優化改進模擬方法的適用性和精確性。目前，IDDES 方法已被廣泛應用于多個領域，比如高速列車氣動性能分析[13]、氫燃料燃燒模擬[14]、水泵內部流動分析[15]以及海上垂直軸風力機[16]的空氣動力學問題等，以上研究充分發揮了IDDES 方法在解決機械運動非定常問題的適用性和精確度。鑒于此，本文采用高精度的IDDES 方法，分析風切變來流條件下大型水平軸風力機流場特征。

1 計算模型

1.1 幾何模型

選擇NH 1 500 kW 風力機作為計算模型，風力機的具體參數見表1。對模型進行簡化，忽略風輪的錐角和仰角，該風力機的幾何模型如圖1 所示。

圖1 風力機幾何模型Fig.1 Geometric model of wind turbine

表1 NH 1 500 kW 風力機參數Tab.1 NH 1 500 kW Wind turbine parameters

1.2 數學模型

1.2.1 數值方法

IDDES 方法結合壁面邊界層區域附近的雷諾時均模擬和非穩定分離區域大渦模擬的特性。計算流動問題過程中用到的三維不可壓縮流體的N-S 方程結合IDDES 方法得到新的湍流輸運方程式為：

式中：Pk為由于黏性力引起的湍流生成項；μt為湍流黏度。μt定義為

IDDES 方法湍流長度尺度lIDDES定義為

式中：β*為經驗系數，取值0.09；F1為混合函數；CDES1與CDES2為通過混合函數計算的模型系數，分別取值0.78 和0.61；Δ為網格尺度，表達式為

式中：dw為計算點到壁面的距離；hmax為單元最大邊長；hwn為垂直于壁面方向的網格尺度。上述IDDES 方法引入lIDDES較好地解決了非定常來流的適用性問題，從而更精確地預測瞬態流動湍動能運輸及耗散情況。

1.2.2 風速沿高度方向的指數律分布

用指數律描述風速沿高度變化規律的分布函數表達式為

式中：V為距離地面高度H處的風速大小，m/s；V0為距離地面參考高度H0處 的風速大小，m/s；α為風剪切指數。本文中H0取輪轂中心高度65 m，V0取額定風速10.4 m/s，α取0.16。

1.3 計算域的劃分

風力機的計算域由圖2 所示的旋轉域和靜止域兩部分組成，原點位于風輪中心，來流速度垂直于XZ平面指向Y軸負方向，計算域入口距離風力機2D(D為風輪直徑)，出口距離風力機5D，左右邊界距離風輪中心各1D。對計算域進行結構化網格劃分，其中旋轉域內葉片表面和邊界層網格作加密處理，網格加密的局部放大如圖3 所示。

圖2 計算域劃分示意圖Fig.2 Computational domain partition diagram

圖3 旋轉域網格加密局部放大圖Fig.3 Local magnification of rotating domain dense mesh

1.4 邊界條件及求解器設置

進口邊界條件選擇速度入口(velocity inlet)，使用用戶自定義函數(user defined functions，簡稱UDF)改變風速在高度方向上的變化，實現剪切入流方式；出口邊界條件選擇自由出流(outflow)；壁面邊界條件為風力機塔筒和機艙設置為無滑移壁面(no slipping wall)，輪轂和葉片設置為旋轉壁面(moving wall)。求解器選擇壓力求解器，采用SIMPLEC 算法實現壓力和速度的耦合。二階迎風格式離散方程，先用RANS 方法計算穩態，待結果收斂后以該結果作為IDDES 方法瞬態計算的初始狀態。由風力機額定轉速17.2 r/min 計算風輪旋轉1°所需時間為0.009 689 9 s，即為瞬態計算的時間步長。監測殘差收斂曲線，當達到設置的收斂標準并且其余監測量穩定變化時，可認為計算結果收斂。

2 計算驗證

2.1 網格無關性驗證

數值模擬中網格密度對計算結果的準確性影響很大，所以做網格無關性驗證排除網格密度對結果的干擾。圖4 給出了額定風速10.4 m/s 下由不同網格數目計算得到的風力機輸出功率以及相對誤差大小。由圖可知，網格數目越多，風力機的輸出功率越大，同時越接近額定功率。當網格數目約2 500 萬時，相對誤差在5%以內，再加密網格相對誤差越來越小。綜合考慮計算時長、計算資源、計算誤差等因素，選定3 200 萬網格數目的數值模擬結果進行分析。

圖4 不同網格數下機組輸出功率Fig.4 The output power of wind turbine under different mesh number

2.2 準確性驗證

為了驗證湍流模型的準確性，將額定風速下模擬得到的推力系數CT、轉矩系數CM與文獻[18]中商用軟件GH Bladed 校核結果作對比，結果如表2所示。從表中可以看出IDDES 方法計算結果稍大于GH Bladed 校核結果，推力系數和轉矩系數的相對誤差分別為6.34%和3.2%，均在10%以內，由此證明IDDES 方法模擬風力機氣動性能可靠性較高。

表2 風力機推力、轉矩系數對比Tab.2 Comparison of wind turbine thrust and torque coefficients

2.3 數值方法對比驗證

從兩種模擬方法的瞬時速度云圖（圖5）對比可以看出，RANS 方法模擬得到的流場在風輪正對高度處速度云圖基本呈對稱狀，而IDDES 方法模擬的流場波動性更強，能捕捉到流場中葉尖渦、葉片尾緣渦和中心渦以及大面積不規則的速度虧損，這種現象更符合風力機真實流場狀況，由此可見IDDES 方法在模擬流場細節方面的優勢遠大于RANS 方法。

圖5 YZ 平面(X=0)瞬時速度云圖Fig.5 YZ plane(X=0) instantaneous velocity contour plots

3 計算結果與分析

為了描述方便，定義風輪旋轉方位角，圖6 為風輪旋轉方位角示意圖。坐標原點位于風輪中心，當圖中1#葉片與X軸垂直并指向Z軸正方向時的方位角記為0°，該葉片繞軸順時針(風力機正前方視角)旋轉一周記為風輪的一個旋轉周期，需要3.488 364 s。

圖6 風輪旋轉方位角示意圖Fig.6 Wind wheel rotation azimuth diagram

3.1 流場速度

速度是描述流場特性的一個重要物理量，流場中任意點的速度可以分解為直角坐標系3 個方向上的速度分量。圖7 給出了風輪旋轉9 周到達0°方位角時流場中速度分量隨高度的變化情況。沿軸向在風力機前后共選取了9 個位置，Z/R=0 對應風輪中心位置，Z/R=±1 對應風輪旋轉上下邊界。

圖7 0°方位角風力機YZ 平面(X=0)速度分量分布Fig.7 Velocity component distribution in YZ plane(X=0) of wind turbine at 0° azimuth

從圖7 可以看出：風力機前后方的速度變化差異很大，前方1D和0.5D處的橫向速度和縱向速度沿高度方向變化量均很小，軸向速度沿高度方向的分布與計算域入口設置的剪切入流曲線幾乎重合，說明前方1D和0.5D范圍內風力機對氣流流速的影響特別小，基本可以忽略；0.1D處3 種速度較另外兩處變化幅度變大，橫向速度沿高度方向的變化呈現“W”型，平均增長15.3%，最大增長39.4%，風輪中心下方的變化幅度小于上方，縱向速度沿高度方向的變化曲線大致呈對稱狀，平均增長4.5%，最大增長94.8%，軸向速度在正對風輪位置處小于來流設定值，而在風輪上方略大于來流設定值，平均虧損率為6%，說明風輪的旋轉效應對機組前方0.1D處的氣流造成擾動。

在風力機后方同一軸向位置處由于組成葉片翼型的氣動特性和翼型所處位置的旋轉線速度不同，所以速度分量沿高度方向發生不同程度的變化，后方0.1D處3 種速度分量變化最明顯，尤其在風輪中心下方(Z/R<0)，由于該位置位于塔筒正后方且距離塔筒很近，受到風輪旋轉和塔影效應的綜合影響，導致該位置處的速度波動極大，下游處風輪中心下方的速度波動程度大于上方。其主要是因為塔筒后方形成一系列渦結構向下游發展，導致速度波動明顯，軸向速度分量在經過風輪之后，后方的速度大小明顯低于來流，表現出較強的尾流速度虧損，越靠近風力機位置，虧損越明顯，軸向速度分量在風輪后發生虧損用于氣流靜壓的恢復，由此形成尾流效應。在風輪旋轉邊界上方，軸向速度分量出現加速現象，主要是由于旋轉的風輪對葉尖處氣流的擠壓作用使得該區域的湍流發展充分，氣流流動速度加快。

總的來說，由于IDDES 方法在模擬瞬態流場方面有很強的優勢，所以用這種方法模擬得到的流場速度波動性較強。越靠近風力機，速度變化越明顯，橫向速度分量與縱向速度分量都變大，而軸向速度分量相較于來流會減小。原因是氣流經過風輪旋轉平面時，由于旋轉平面的誘導作用使得氣流速度在軸向方向發生虧損，軸向速度在風輪前發生虧損用于提高氣流靜壓推動風輪轉動作功，軸向速度分量在風輪后發生虧損用于氣流靜壓的恢復，由此形成尾流效應。

3.2 脈動速度頻譜分析

為了描述風輪旋轉時流場的運動特性，選3 個典型高度(Z=41.5、0、-41.5 m)處風力機前后方不同位置(圖8(a))的瞬時風速在4 個旋轉周期(17.44182～31.395 276 s) 內的時域變化，監測結果如圖8(b)所示。

圖8 流場監測點分布(a)與監測點速度(b)Fig.8 Monitoring point distribution(a) and velocity(b) in flow field

風力機前方的速度變化較平穩，后方速度波動明顯，41.5 m 高度處的速度波動呈現周期性，而0 m高度處受風輪旋轉產生中心渦和-41.5 m 高度處受葉片葉尖渦與塔筒的塔影效應綜合影響使得速度變化不存在周期性。用脈動風速進行快速傅里葉變換來描述監測點的湍流功率譜特性，脈動風速的功率譜結果如圖9 所示。由風力機額定轉速n=17.2 r/min 計算可得風輪轉頻fn為0.286 7 Hz，風輪旋轉一周，葉片通過的頻率fp為3 倍的風輪轉頻，即0.86 Hz。從圖9 可以看出：風輪前的脈動風速功率譜在任意高度處都存在明顯波峰，對應頻率均在fp的整數倍處，風輪后的脈動風速功率譜只在41.5 m 高度處有顯著波峰，0 與-41.5 m 高度處的監測點處存在的大尺度旋渦影響了功率譜的譜峰分布；速度波動程度越大對應功率譜密度越大，湍動能越大，具體表現為越靠近風力機脈動速度的功率譜密度越大；風輪前的功率譜密度比風輪后方的低，對比風輪后方0.5D距離處功率譜沿高度方向的變化發現，Z=0 高度的脈動速度功率譜密度最大，說明中心渦對功率譜密度的影響最明顯。

圖9 脈動風速的功率譜Fig.9 Power spectrum of fluctuating wind velocity

3.3 輸出功率分析

來流條件的不穩定決定了風力機的輸出功率會隨時間不斷發生變化，但這種變化并不是毫無規律的。對比分析一個旋轉周期內剪切入流和均勻入流條件下的輸出功率隨方位角波動情況，結果如圖10 所示。由圖可以看出兩種入流條件下功率輸出都呈現余弦變化，在1#葉片的方位角分別處于60°、180°、300°時輸出功率出現極小值，此時均勻入流和剪切入流條件下的輸出功率分別只有額定功率的95.6%和94.1%。這是因為此時3 支葉片正分別經過塔筒正前方，繞葉片流動的氣流與繞塔筒流動的氣流相互摻雜影響了風力機氣動力的產生，總的來看均勻流下的輸出功率大于剪切流。該段時間內均勻流和剪切流的平均輸出功率分別為1 510 kW 和1 490.425 kW，風能利用系數分別為0.426 和0.421。為了比較兩種入流條件下風力機輸出功率的功率譜，選定風輪4 個旋轉周期內的轉矩數據作為樣本，采樣時間間隔為0.009 689 9 s，即風輪旋轉1°采樣一次，得到風輪轉矩隨時間波動如圖11 所示。由圖可以看出風輪轉矩呈現出脈動特性，在采集樣本期間出現了12 個極小值，正好對應每個葉片經過塔筒正前方，說明塔筒的繞流效應會影響風輪的受力狀況。

圖10 輸出功率隨方位角的變化Fig.10 Variation of output power with azimuth angle

圖11 4 個旋轉周期內風輪轉矩變化Fig.11 The torque of the wind wheel during four rotation cycles

為了反映風力機流場能量輸運規律，圖12 給出了兩種入流情況下的功率譜密度。由圖可以觀察到兩種情況具有相近的頻譜變化，功率譜密度隨著頻率的增大而減小。氣流受風輪旋轉擾動以及剪切風速的不均勻性影響，流場內的湍流充分發展，其能量譜分為含能區、慣性子區和耗散區3 個部分。其中，含能區生成的大尺度旋渦含有湍動能的絕大部分能量，在慣性子區大尺度旋渦輸出能量擴散成眾多靠慣性輸入能量的小尺度旋渦，并伴隨著湍動能逐級連續傳遞，呈現常見的-5/3 衰減斜率，最終旋渦尺度太小且受流體黏滯力的影響能量損失不斷增大，在耗散區湍流能量耗散消失。

圖12 輸出功率的頻譜特性Fig.12 Frequency spectrum characteristic of output power

4 結論

本文通過IDDES 方法模擬大型風力機運行的流場特性，得出以下結論。

1)在風力機前方，風輪旋轉對氣流造成的擾動沿高度方向分布比較規律，其中橫向速度分量和縱向速度分量都大于入口設定值且呈對稱分布，軸向速度分量在風輪旋轉平面高度內相較于入口設定變小而在旋轉邊界以外會增大。在風力機前方速度虧損提高氣流靜壓推動風輪旋轉做功，后方流場中出現速度虧損用于恢復氣流靜壓。

2)風力機前方流場內的速度隨時間變化較平穩，后方速度呈現明顯的脈動特性且離地面越近變化幅度越大。脈動風速的頻譜圖顯示低頻區譜峰的出現時刻與葉片轉頻的整數倍有關，速度波動越大對應的功率譜密度越大。

3)風力機的輸出功率呈余弦變化，風切變的存在使得葉片表面受力不均勻從而機組功率和風輪載荷高度波動，在三支葉片分別經過塔筒正前方時出現極小值。同樣條件下均勻流的輸出功率大于剪切流，兩種入流方式的能量變化趨勢一致。在含能區大尺度旋渦含有大量湍動能，在慣性子區大尺度旋渦輸出能量擴散成小尺度旋渦并滿足-5/3 衰減斜率，在耗散區小尺度旋渦受黏滯力影響湍動能耗散消失。

本文總結了剪切入流條件下大型水平軸風力機的流場和功率特性，充分驗證了高精度的IDDES方法模擬計算大型水平軸風力機的可行性，為未來該方法應用于更加復雜的流動條件提供參考。