基于大渦模擬的風力機尾流特性研究

楊從新，張亞光，張旭耀，何 攀

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

1 引言

風力機尾流中存在葉尖渦、葉根渦等復雜的渦系，它們之間存在大尺度相干結構決定的分離流動，對風力機組的安全運行提出挑戰，研究風力機尾流特性，對合理開發與研究新型風力機和合理利用風資源有重要指導作用。近年來，大渦模擬在風力機尾流與風場的高精度數值模擬中的應用越來越廣泛。文獻［1］首次耦合大渦模擬與致動線方法研究風力機尾流。文獻［2］用致動線方法研究了均勻來流條件下的風力機近尾流特性，對渦結構、速度虧損、擬渦能進行分析，驗證了致動線模型的準確性。文獻［3］將致動線方法的數值模擬結果與BEM理論、CFD方法以及風洞實驗進行了系統的比較，驗證了致動線方法用于風力機氣動數值模擬的可行性。文獻［4］基于Lagrangian動力亞格子模型對風力機尾流場進行數值模擬，分析了風力機近尾流場與遠尾流場的流動特性。文獻［5］利用致動線方法分析了不同入流條件下風切變、湍流強度等因素對風力機尾流特性的影響。文獻［6］對充分發展的風電場邊界層進行數值模擬，發現高雷諾數流比低雷諾數時對SGS模型的選擇更敏感。文獻［7］引入周期性擾動，研究了小型風力機葉尖渦的穩定性及其演變，發現周期性擾動的存在對渦旋的發展與破碎過程有明顯的影響；在渦旋破裂前，還觀察到了明顯的渦旋配對現象，并通過改變強迫擾動的幅值和頻率進行了參數研究，評價其對葉尖渦的影響。文獻［8］對水平軸風力機在水槽中的近尾流進行實驗研究，探討了水平軸風力機尾跡的旋渦相互作用及螺旋形渦絲的穩定性。文獻［9］對單臺風力機進行數值模擬研究，受葉尖附近的低振幅激發源的擾動，沿螺旋運動的波的放大會觸發不穩定，導致尾流破碎，重點研究了從葉片后緣脫落渦的穩定性。

以上研究使用單一湍流模型對風力機進行模擬，不同亞格子模型的選取會對研究對象數值模擬的結果產生影響。這里以流場中的單一風輪為研究對象，排除塔架、機艙等結構對風力機尾流的影響，探討均勻入流條件下采用三種亞格子模型的風力機尾流特性。

2 研究對象與計算方案

2.1 研究對象

以一臺33kW兩葉片水平軸風力機為研究對象，由于尾流的旋轉效應，葉尖使用Glauert修正，風力機基本參數如表1所示。

表1 33 kW風力機基本參數Tab.1 The Parameters of 33 kW Wind Turbine

2.2 計算方案

計算域及加密區域網格尺寸，如圖1所示。圖中：D-風輪直徑。風輪位于進口后方2D處，為減少壁面效應，根據文獻［10］提出的準則，設置計算域尺寸為22D×6D×6D；同時，為減小流動方向網格加密界限處不同網格尺寸對計算結果的影響，網格加密區域為入口到出口，沿風輪徑向經過二次加密后的六面體網格總數為4.356×107，最小網格尺寸為0.25m。

圖1 計算域及網格加密Fig.1 Computing Domain and Mesh Refinement

計算時間步長為0.005 s；大渦模擬對亞網格尺度流動進行模擬，重點關注尾流渦結構，無需進行網格無關性驗證。這里數值計算中風力機葉尖速比l分別取值13.17、8.23、5.98，對應入流速度分別為5、8、11 m·s-1。

3 數值模型

耦合致動線與大渦模擬方法，采用開源CFD軟件Open-FOAM對處于均勻入流條件下的風力機做數值模擬研究。

3.1 致動線方法

風力機的模擬用致動線方法，其基本概念為每個風力機葉片用一條承受體積力的線表示，將葉片沿展向分割為若干葉素，分別求解葉素上的體積力。葉素上的升力、阻力由下式計算為：

式中：eL-升力方向向量；eD-阻力方向向量。

各致動元上離散的點力通過高斯正則化函數將致動力反作用于流體：

點（x，y，z）處的體積力為：

式中：di-第i個致動點（xi，yi，zi）的中心與投射點（x，y，z）的距離；ε-高斯光順參數。

致動點處網格尺寸為0.25m，高斯光順參數ε=0.5，風輪直徑上分布致動點個數為n=60，以使致動點間的距離小于風輪處網格尺寸。

3.2 大渦模擬

對Navier-Stokes方程施加卷積過濾運算，濾波后的LES控制方程如下：

式中：ρ-空氣密度；ν-分子粘度濾波后的速度分量；p-修正壓力；fε-由于風輪的存在作用于空氣的體積力；τij-SGS應力，定義為

3.3 湍流模型

3.3.1 Smagorinsky模型

計算亞格子應力τij有不同的方法，常用的方法是利用Boussinesq渦粘假設。將τij表示為正應力與偏向應力，正應力部分被認為各向同性，求解Smagorinsky模型的表達式為：

正應力和亞格子動能kSGS聯系起來，則：

所以，亞格子應力可表示為：

式中：δij-克羅內克符號，當i=j時δij=1，i≠j時δij=0；

νSGS-亞格子尺度粘度應變率張量。

Smagorinsky模型中，νSGS表示為：

將上式代入方程（8），則方程封閉參數為kSGS，

式中：Ck、Ce-Smagorinsky模型無量綱系 數，Sm模型取值Ck=0.0676、Ce=0.93；Δ-濾波器過濾尺度。

3.3.2 Lagrangian動力模型

亞格子動力模型對Navier-Stokes方程施加兩次過濾尺度不同的過濾運算，從解析流中動態計算模型參數，并允許模型參數在空間和時間上發生變化，亞格子動力模型基于Germano等式：

求SGS模型系數Cs的動態過程使Germano等式誤差最小化：

由文獻［11］等人提出的初始動態模型滿足?ijSij=0得到Cs。文獻［12］發現，當?ij在最小二乘法上最小化時，方程表現良好，從而得到結果：

式中：“＜＞”-系綜平均。

文獻［13］提出了Lagrangian動力亞格子模型，該模型沿流體質點運動軌跡做統計平均最小化Germano等式誤差。

本次數值實驗采用以下亞格子模型進行數值計算：Smagorinsky模型（Smagorinsky model）、Lagrangian動力模型（Lagrangian dynamic model）；另外當νSGS=0時沒有明確表示湍流粘性，此時引起動能耗散的唯一影響因素是數值耗散。以上三種模型在這里中分別以Sm（Smagorinsky model）、dyL（Lagrangian dynamic model）、Nom（No model，νSGS=0）表示。

3.3.3 數值離散與邊界條件

采用開源軟件OpenFOAM進行數值計算，N-S方程中的原始變量用有限體積法進行離散，求解方程時使用瞬態壓力全隱分離式（PISO）算法。時間項采用二階有界隱式Crank-Nicholoson格式，梯度項、散度項等空間項采用高斯線性離散。進口為速度入口，方向沿X軸，出口為壓力出口，壁面處采用滑移邊界條件。

4 計算結果分析

本次數值實驗采用甘肅省計算中心高性能計算集群的兩節點48核心計算，其中，計算工況λ=13.17時，設置總計算時間為110s，三種亞格子模型分別花費計算時間180.7h（dyL）、140.2h（Sm）、158.8h（Nom）。亞格子動力模式需要進行統計平均，理論上需要進行系綜平均，非常花費計算時間；但拉格朗日動力模式沿質點軌跡平均確定模型系數，增加的計算量不多。當尖速比l分別取值13.17、8.23、5.98時，對應的風輪輸出功率為8.016、16.623、33.267 kW，額定風速下的功率誤差為0.809%，滿足計算要求。當葉尖速比l=13.17時，三種SGS模型在計算過程中的最大庫郎數Co（Courant Number）分別為0.135（dyL）、0.136（Sm）、0.149（Nom）。

4.1 尾流發展

Lagrangian動力模式允許在不調整任何參數的情況下，根據流動在時間和空間上的變化，對模型系數和參數進行動態計算。由于渦粘模型對渦生成的直接影響，不考慮輪轂、塔架、機艙等對風力機尾流的影響。

由圖2（a）可知，尾流從葉尖渦與葉根渦以近似對稱的渦旋結構向下游傳播，在遠尾流區，由于流體具有粘性，在周圍氣體的作用下尾流的對稱結構被打破。葉片失速前，λ=8.23、l=5.98時隨葉尖速比減小，攻角增大，葉片吸力面與壓力面間的壓差增大，導致葉尖渦更強具有較高的穩定性。λ=13.17時，尾流區呈現完整的近尾流區、尾流發展區與遠尾流區，從尾流發展的角度看具有代表性，將主要對此工況下不同SGS模型的風力機尾流進行研究。

如圖2（b）所示，為λ=13.17時使用不同SGS模型獲得的垂直于風輪旋轉平面的渦旋云圖。近尾流區，尾渦從葉尖與葉根處向下游傳播；之后三種亞格子模型下的葉尖渦均出現了K-H（Kelvin Helmholz）不穩定性現象，與文獻［14］中實驗觀測到的葉尖渦橫截面尾跡相似；最后由于湍流的能量級串，大尺度尾渦在遠尾流區耗散變成小尺度渦。

圖2 使用dyL模型垂直于風輪旋轉平面的渦旋圖Fig.2 2D Snapshots with Vorticity Contour in the Vertical Plane in dyL Model

4.2 尾流速度分布及其湍流特性

實際流體中，由于流體具有粘性，尾流終將恢復；數值模擬中，由于數值耗散，尾流也終會恢復。由圖3中使用dyL模型在不同葉尖速比下的速度分布曲線可知，λ=5.98時，直到風力機尾流下游11D處速度有微弱減小；λ=13.17時尾流恢復最快。由于周圍空氣的作用，尾流逐漸膨脹發展為完全湍流。從三種SGS模型預測的尾流速度分布曲線（圖4）可以看出，0D-2D之間的近尾流區，不同SGS模型預測的尾流速度廓線幾乎相同，在尾流發展區與遠尾流區，三種SGS模型預測的尾流速度廓線交替上升，不能說明哪種亞格子模型的尾流預測能力更突出，其中尾流速度廓線在7D后均由原來的“倒鐘形”發展成為近高斯分布。

圖3 使用dyL模型不同葉尖速比下風力機尾流時均軸向速度分布曲線Fig.3 Profiles of the Time Averaged Axial Velocity for Different Tip-Speed Ratios in dyL Model

圖4 λ=13.17時不同亞格子模型下風力機尾流均軸向速度分布曲線Fig.4 Profiles of the Time Averaged Axial Velocity for Different SGS Models in Tip-Speed Ratio of λ=13.17

為分析尾流場中湍流的流動特性，引入無量綱的雷諾應力與湍流切應力分別定義為，引入湍動能：

式中：u′、v′、w′-軸向脈動速度、橫向脈動速度、垂向脈動速度。

對于Nom模型，νSGS=0，所預測的雷諾應力與湍流切應力為零，圖5中不顯示。λ=13.17工況下的風力機尾流中，雷諾應力與湍流切應力在葉尖渦與葉根渦區有明顯的最大值與極大值；隨著尾流向下游發展，尾流逐漸恢復，雷諾應力與湍流切應力逐漸變小；湍動能與雷諾應力的變化趨勢一致，說明尾流場中正應力對湍動能的貢獻最大，含有更多的能量；在近尾流區，由dyL模型產生的雷諾應力與湍流切應力大于Sm模型產生的雷諾應力與湍流切應力，隨尾流發展兩模型產生的應力趨于接近，表明尾跡隨下游位置的變化呈現出逐漸增大的各向同性。比較圖4～圖6，在13D處，雷諾應力與切應力接近于零，平均速度虧損比產生的湍流更持久。

圖5 葉尖速比λ=13.17時不同亞格子模型風力機尾流的雷諾應力（5a）切應力（5b）分布曲線Fig.5 Profiles of the Reynolds Stress（5a）and Shear Stress（5b）for Different SGS Models in Tip Velocity Ratio of λ=13.17

圖6 λ=13.17時不同亞格子模型風力機尾流湍動能分布曲線Fig.6 Profiles of the TKE for Different SGS Models in Tip Velocity Ratio of λ=13.17

4.3 軸向速度干擾因子與速度環量

如圖7所示，為軸向速度干擾因子ax=1-Vx/V∞沿葉片徑向的時均分布曲線。由圖可知，葉尖速比l=13.17時軸向干擾因子從葉尖到葉根增大過快，在0.85R處過大，其值接近0.5，此處的速度環量也最大，此時流動處于湍流狀態，處理此問題通過傳統的BEM理論無法解決，只能通過經驗進行修正；l=8.23時，在0.2R～0.9R的葉片重載區，軸向干擾因子增長緩慢；l=5.98時，相比于大葉尖速比工況，ax變化范圍最小。

圖7 軸向干擾因子沿葉片的徑向分布Fig.7 Radial Distribution of the Axial Interference Factor along the Blades

如圖8所示，為環量Γ=L/（ρ/Vrel）沿葉片徑向的時均分布。圖中，l=13.17時葉尖與葉根部分存在較大速度梯度，對應風輪處則存在較強的葉尖渦與葉根渦；且λ=13.17、λ=8.23時葉片重載區的環量變化較小；λ=5.98時環量沿葉片展向不斷變化，在葉片中部部分區域大于λ=8.23工況，反映到風輪上表示兩種工況下葉片推動風輪旋轉的主要升力貢獻區域發生變化，相比λ=13.17工況，葉尖渦與葉根渦強度較弱。使用三種亞格子模型所得到的軸向干擾因子與環量沿葉片徑向的分布幾乎相同，在葉片重載區有微小波動，且與λ=8.23、λ=5.98工況下得到的結論一致。

圖8 環量沿葉片的徑向分布Fig.8 Radial Distribution of Circulation along the Blades

4.4 尾渦分布

使用三種亞格子模型均能很好地獲得風力機尾流場結構。如圖9所示，為使用Q準則所得使用dyL模型在λ=13.17工況下風力機速度梯度第二不變量的不同Q值三維等值面，由圖可知隨著Q值增大，尾渦的圓柱形膨脹區直徑逐漸變小，位于尾渦內部的高渦量層仍然能夠被很好地捕捉，但風力機下游強度較弱的渦旋逐漸消失；相同位置處的尾渦由不同強度的渦旋組成，強度高的渦旋其對應的渦管直徑較小，尾渦膨脹之前能量較高，能量較小的渦旋是尾渦組成中的主要部分。

圖9 λ=13.17時使用dyL模型尾流不同Q值三維等值面Fig.9 3D Snapshots of Q Using dyL Models in λ=13.17

由曲線圍成的閉環不規則圖形（Q值二維等值面，圖10）代表三維Q值等值面與Y-Z面的交線，其物理意義為渦量的模，適合提取邊界層之外的渦結構；在遠尾流區，根據其統計特性可以獲得尾渦中各渦旋的大小與渦流間隙。

圖10 使用dyL模型在λ=13.17工況下風力機尾流不同位置處的速度梯度第二不變量Q值二維等值面（Q=0.005）Fig.10 2D Snapshots of Q at Different Downstream Positions Using dyL Models in l=13.17（Q=0.005）

尾流從葉尖渦與葉根渦向下游傳播，由圖10可知在3D后尾渦開始潰散，不同Y-Z截面處的二維Q值等值面在均勻流下可視為近似圓形，隨尾流位置的變化，圓的直徑不斷變小，說明尾流發展的過程中，強度較大的渦集中在尾渦中心，強度較小的渦在尾渦外圍分布；若風力機在無限長區域內運行，隨著尾流逐漸恢復，尾流中的湍流將逐步轉化為入流狀態的均勻流，不同強度的渦將逐漸消失。

5 結論

（1）雷諾應力與湍流切應力在葉尖渦與葉根渦區有明顯的最大值與極大值，且由Lagrangian動力亞格子模型比Smagorinsky模型預測的應力大，隨著尾流向下游發展，雷諾應力與湍流切應力逐漸變小，尾跡隨下游位置的變化呈現出逐漸增大的各向同性；湍動能中正應力的貢獻最大；在均勻入流條件下，平均速度虧損比產生的湍流更持久。

（2）相同位置處的尾渦由不同強度的渦旋組成，能量較小的渦旋是尾渦組成中的主要部分。尾流發展的過程中，強度較大的渦集中在尾渦中心，強度較小的渦在尾渦外圍分布。

（3）對于所使用的三種亞格子模型，能預測出相似的尾流效應，亞格子模型的選擇對尾流的模擬影響較小；νSGS=0時，求解過程只有數值耗散，仍能獲得很好的模擬結果。