基于大渦模擬的風(fēng)力機尾流特性研究

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

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

1 引言

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

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

2 研究對象與計算方案

2.1 研究對象

以一臺33kW兩葉片水平軸風(fēng)力機為研究對象，由于尾流的旋轉(zhuǎn)效應(yīng)，葉尖使用Glauert修正，風(fēng)力機基本參數(shù)如表1所示。

表1 33 kW風(fēng)力機基本參數(shù)Tab.1 The Parameters of 33 kW Wind Turbine

2.2 計算方案

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

圖1 計算域及網(wǎng)格加密Fig.1 Computing Domain and Mesh Refinement

計算時間步長為0.005 s；大渦模擬對亞網(wǎng)格尺度流動進行模擬，重點關(guān)注尾流渦結(jié)構(gòu)，無需進行網(wǎng)格無關(guān)性驗證。這里數(shù)值計算中風(fēng)力機葉尖速比l分別取值13.17、8.23、5.98，對應(yīng)入流速度分別為5、8、11 m·s-1。

3 數(shù)值模型

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

3.1 致動線方法

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

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

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

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

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

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

3.2 大渦模擬

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

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

3.3 湍流模型

3.3.1 Smagorinsky模型

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

正應(yīng)力和亞格子動能kSGS聯(lián)系起來，則：

所以，亞格子應(yīng)力可表示為：

式中：δij-克羅內(nèi)克符號，當(dāng)i=j時δij=1，i≠j時δij=0；

νSGS-亞格子尺度粘度應(yīng)變率張量。

Smagorinsky模型中，νSGS表示為：

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

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

3.3.2 Lagrangian動力模型

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

求SGS模型系數(shù)Cs的動態(tài)過程使Germano等式誤差最小化：

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

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

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

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

3.3.3 數(shù)值離散與邊界條件

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

4 計算結(jié)果分析

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

4.1 尾流發(fā)展

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

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

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

圖2 使用dyL模型垂直于風(fēng)輪旋轉(zhuǎn)平面的渦旋圖Fig.2 2D Snapshots with Vorticity Contour in the Vertical Plane in dyL Model

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

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

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

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

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

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

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

圖5 葉尖速比λ=13.17時不同亞格子模型風(fēng)力機尾流的雷諾應(yīng)力（5a）切應(yīng)力（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時不同亞格子模型風(fēng)力機尾流湍動能分布曲線Fig.6 Profiles of the TKE for Different SGS Models in Tip Velocity Ratio of λ=13.17

4.3 軸向速度干擾因子與速度環(huán)量

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

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

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

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

4.4 尾渦分布

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

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

由曲線圍成的閉環(huán)不規(guī)則圖形（Q值二維等值面，圖10）代表三維Q值等值面與Y-Z面的交線，其物理意義為渦量的模，適合提取邊界層之外的渦結(jié)構(gòu)；在遠尾流區(qū)，根據(jù)其統(tǒng)計特性可以獲得尾渦中各渦旋的大小與渦流間隙。

圖10 使用dyL模型在λ=13.17工況下風(fēng)力機尾流不同位置處的速度梯度第二不變量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值等值面在均勻流下可視為近似圓形，隨尾流位置的變化，圓的直徑不斷變小，說明尾流發(fā)展的過程中，強度較大的渦集中在尾渦中心，強度較小的渦在尾渦外圍分布；若風(fēng)力機在無限長區(qū)域內(nèi)運行，隨著尾流逐漸恢復(fù)，尾流中的湍流將逐步轉(zhuǎn)化為入流狀態(tài)的均勻流，不同強度的渦將逐漸消失。

5 結(jié)論

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

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

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