陣風工況下多臺風力機尾流效應的非定常特性

曹九發，宋佺珉，王超群，朱衛軍，柯世堂

(1. 揚州大學 電氣與能源動力工程學院，揚州 225127；2. 丹麥科技大學 風能系，丹麥 哥本哈根 靈比 2800；3. 南京航空航天大學 民航學院，南京 210016)

0 引言

風力機尾流影響是風電場的選址和布局設計中的一個重要考量因素。在風電場的整體布局中，一部分風力機不可避免地會處于上游風力機的尾跡中，由此帶來的直接影響是位于尾跡中的風力機發電量降低，同時尾流效應帶來的風切效應和湍流增強等因素增加了尾跡區風力機的疲勞載荷。丹麥某個海上風電場的測試研究中發現，由于尾流效應的影響，整個風電場的發電功率損失可達到總裝機功率的10%～20%[1]。風力機工作在非穩態的復雜來流環境中，其性能會受到較大影響，尤其是極端風況陣風（即風速在短時間內產生劇烈波動），會導致風力機氣動載荷發生破壞性波動，影響風力機的運行安全和使用壽命。因此，研究風力機尾流特性及上下游風力機的載荷特性，著眼于提高發電效率、減少尾流對風力機性能的影響，對風電場高效健康運行具有十分重要的意義。

目前風力機尾流研究主要有實驗研究和數值模擬研究。國內外很多學者[2-4]對風力機尾流進行了PIV實驗研究，得到了尾流流場數據，但是實驗研究大多為縮比模型的風洞實驗，而風力機尾流的場外實驗，特別是針對多臺風力機在復雜非穩態工況下的尾流干擾實驗研究，成本較高，難度大。隨著計算方法及設備性能的日益發展，CFD已經成為研究風力機尾流效應和尾流干擾的常用手段[5-9]。近年來，致動理論耦合CFD方法成為風力機尾流數值模擬方法的熱點[10-12]。該耦合方法將風力機葉片的氣動力作為體積力源項加載到流場中，然后通過求解N-S（Navier Stokes）方程來獲得流場信息，這種方法相對于傳統具有葉片物面網格的方法，可以大大提高數值計算效率。致動理論模型主要有致動盤（Actuator Disk，AD）、致動線（ Actuator Line，AL）和致動面（Actuator Surface，AS）模型。丹麥科技大學的S?rensen[13]在致動盤模型的基礎上發展出了風力機尾流模擬的致動線模型，并采用自主研發的EllipSys3D求解器，對比了風力機在大氣湍流和均勻來流下的尾流特性。Zhu等[14]采用EllipSys3D求解器，研究了在葉片根部、風輪前方增加導流盤后風力機功率的變化。Shen等[15]將致動線模型與N-S/LES（Large Eddy Simulation）耦合，然后基于MEXICO風力機實驗數據，開展了風力機尾流數值模擬方法的改進研究。朱翀等[16]以NH1500 葉片為研究對象，采用RANS（Reynolds Averaged Navier-Stokes）結合k-ωSST湍流模型的方法，從葉片載荷分布和功率系數兩個方面將致動線方法、葉素動量理論、葉片物面網格方法進行了對比，同時分析了風力機尾流速度變化特性與氣動載荷特性。Qian等[17]采用 ALM-LES 方法，模擬了風力機尾流的高精度湍流流動結構及氣動載荷。卞鳳嬌等[18]基于OpenFOAM 平臺，對比了PISO求解的致動線模型和多重參考穩態求解的CFD實體模型，進而采用致動線模型數值研究了NREL 5 MW風力機在均勻來流中的尾流場特性。Nathan等[19]比較了采用OpenFOAM和EllipSys3D求解器求解致動線模型的結果，并將其與NEW MEXICO實驗數據進行對比，結果表明兩者均能在近尾流處較好地吻合實驗數據。周洋等[20]比較了致動線方法中均勻分布和高斯分布兩種不同的體積力分布方式，并分別采用這兩種方式對Nibe A型風力機進行了數值模擬，得到了高斯分布方式更接近實驗數據的結論。目前風力機尾流效應研究和尾流干擾研究主要是針對穩態工況開展的，本文則計劃針對風電場非定常陣風入流工況，采用LES耦合致動線模型的方法，開展多臺風機的尾流特性、尾流效應及風力機動態氣動載荷的相關研究。

1 數值模型與方法

1.1 N-S/LES模型

數值模擬的控制方程采用三維不可壓N-S方程。連續性方程為：

對于三維不可壓縮牛頓流體，動量守恒方程為：

其中：v為 慣性坐標系下的速度矢量；t為時間；ρ為密度；p為靜壓；τ為應力張量；S為添加的體積力動量源項，本文主要是風力機葉片產生的氣動力。

τ由下式給出：

其中 μ為分子黏性系數，I為單位張量。

LES大渦模擬是求解瞬態流動的主要方法之一，它是在空間域上對瞬態Navier-Stokes方程進行濾波處理。濾波過程可將小于濾波器寬度或計算網格尺度的小尺度渦過濾出去，從而形成關于大尺度渦的控制方程。在LES方法中，通過使用濾波器函數，將每個變量都分成兩部分—大尺度分量和小尺度分量。大尺度部分是濾波后的變量，是在LES模擬時直接計算的部分；小尺度部分是需要通過模型來表示的。濾波處理后，瞬態狀態下的N-S方程包含了如下新的應力項：

式 中， τij被 定 義 為 亞 格 子 尺 度 應 力（Sub-Grid-scale Stress，SGS），體現了小尺度渦的運動對所求解的運動方程的影響。 τij是一個對稱張量，其包括6個獨立的未知變量。對 τij進行不同的建模，就可以得到不同的SGS模型。本文采用Smagorinsky-Lilly模型進行風力機尾流流場的LES數值模擬。整個風力機流場數值求解流程如圖1所示。

圖1 致動模型與N-S方程耦合計算流程圖Fig. 1 Flow chat of the actuator line model coupled with the N-S equations

1.2 致動線模型AL

風力機的致動線模型是基于風力機的葉素動量理論延展出的等效風力機葉片模型[15]。如圖2所示，將葉片分成一個個微元段（葉素），根據當地的流場信息與已知的翼型性能數據、葉片外形數據，可求得每個葉素上的氣動力。圖2最右邊為某個葉素的截面圖，x為來流方向。

圖2 致動線模型示意圖Fig. 2 Schematic of the actuator line model

葉片翼型截面的當地速度[15]為：

其中， ?為風輪轉動的角速度，r為截面翼型到葉根的位置，Vx和Vθ分別為翼型的軸向速度和切向速度。

翼型截面當地速度與風輪平面夾角為入流角，

翼型截面的當地攻角為入流角減去槳距角γ，

確定了葉片翼型的當地速度與攻角后，即可求出葉片單位展長的升力與阻力：

其中，nb為 葉片數，CL和CD分別為升力、阻力系數，c為翼型弦長，ρ為空氣密度。

為了考慮機艙和塔架對風力機尾流流場的影響，也可把機艙等效于體積力源項來模擬。機艙體積力計算公式為：

其中：CD,nac為機艙阻力系數（取值范圍一般為0.8～1.2），本文取值1.0；Anac為機艙的橫截面積。

參考風力機葉片葉素理論方法，可以把塔架同樣分成一段段單個長度為dh的微元。任意一段微元處的軸向阻力為：

式中：CD,tower為塔架阻力系數（取值范圍一般在0.8～1.2），本文取值1.0；dtower為微元段塔架直徑。

確定了所有體積力之后，需要將計算出來的力光順到葉素點周圍的網格中去。為了防止求解時發生數值震蕩，本文采用三維高斯分布將體積力光順過渡到周圍的網格上：

其中：dm為葉素點到網格型心的距離；ε為高斯分布因子，該參數取值與風輪旋轉加密區網格單元尺寸Dm有 關，即ε=nDm，n為倍數，本文取n= 2。

1.3 陣風

風況是結構完整性設計主要考慮的外部條件。從載荷和安全角度考慮，風況可分為風力發電機組運行期間通常的正常風況以及按1年或50年重復周期確定的極端風況。極端風況包括由于風向和風速急速劇烈變化產生的峰值風速，因此在設計時需要對其影響加以考慮。一般采用極端風況確定風力機的極端載荷。

對標準等級風力機組，根據IEC61400-1風力發電設計要求規范，輪轂高度處的重復周期為N年的陣風幅值關系式如下：

式中：I15= 0.18，a= 2，為較高湍流強度類別；I15= 0.16，a= 3，為較低湍流強度類別。本文選取較高湍流類別。 Λ1為湍流尺度參數，取值由下式給出：

其中：Drotor是風輪直徑；β = 4.8，對應N= 1； β = 6.4，對應N= 50。

當重復周期為N年時，由下列方程確定風速：

式中：Vz選用風切；T= 10.5 s，對應N= 1；T= 14.0 s，對應N= 50。本文選取T= 10.5 s。

圖3為極端運行陣風示意圖，假設正常風速為11.4 m/s（亦為本文作為邊界條件的輸入風況）。ag點表示正常風速值的位置，bg、dg點為陣風周期內風速最低的時間點位置，cg點為陣風周期內風速最高的時間點位置，eg點表示陣風周期結束后風速回歸到正常風速值的位置。

圖3 陣風示意圖Fig. 3 Schematic of the gust wind

2 結果與分析

2.1 算例驗證

本文以丹麥Nibe風力機作為致動線模型數值模擬的驗證算例。Nibe風力機的主要參數見表1，具體的葉片、輪轂和塔架參數等見文獻[21]。

表1 Nibe B風力機主要參數信息表Table 1 Key parameters of the Nibe B wind turbine

計算采用如圖4所示的長方體計算域，計算域的長、寬、高分別為24.5R、6R、6R。對轉子所在的區域進行局部加密，加密范圍為：x方向，風輪平面前后0.5R；z方向，輪轂中心左右1.5R；y方向，輪轂向上1.5R。加密區以外的單元，網格按照一定的比例向外延伸。風力機轉子中心在（0，45，0）位置。整個計算域總網格量為430萬，均采用六面體正交網格。

圖4 計算域網格示意圖Fig. 4 Schematic of the computational domain and grid

風力機的葉片被分為20段，每一葉素段長度為1 m，恰好與加密區的網格單元尺寸相同。在計算時，取三維高斯分布參數ε= 2 m。Nibe風力機位于丹麥北部沿海地區，西側為淺水區，東側為以矮草為主的平坦地形，可認為來流風速不受地形影響。因此，根據實驗[21]工況，入口速度邊界速度為8.0～9.1 m/s工況（本文選取8.55 m/s），湍流強度為10%～15%（本文選取10%），其中入口風速考慮風剪切效應，出口邊界設置為壓力出口，下邊界為物面，上邊界為滑移邊界，左右邊界設置為滑移邊界。對采用笛卡爾正交坐標系的致動線模型而言，非定常數值計算的時間步長由風輪轉速確定，即單個時間步內葉尖旋轉的距離不能大于所在區域的網格單元尺寸，Dm： ωrdt≤Dm。在本算例中，時間步長為dt≤ 0.006 s，總共模擬400 s，尾流風速由200～400 s數據取平均值得到。

圖5為Nibe風力機在軸向上不同位置處的速度剖面圖。通過與實驗測量值對比，可以看出：風力機尾流速度在軸向位置上呈現出拋物線特性的變化趨勢，輪轂中間位置存在速度最低點；隨著軸向位置的增加，風力機尾流速度不斷恢復，這是由于隨著軸向距離增加，周圍大氣向尾流中心不斷進行能源補充，使得速度逐漸向來流速度值恢復。從圖中也可以看出，本文的計算方法，不僅近尾流區的速度計算結果有較好的準確度，而且在比較重要的遠尾區，計算結果與實驗值吻合較好。因此，本文的非定常致動線模型耦合LES的風力機尾流數值計算方法，具有較好的準確性和可行性。

圖5 不同位置的風力機尾流的剖面速度對比圖Fig. 5 Velocity profile comparison at different locations of the wind turbine wake

2.2 NREL 5 MW陣風工況數值計算

2.2.1 邊界條件和計算域

NREL 5 MW風力機主要參數見表2，具體的葉片和翼型數據可參閱文獻[22]。

表2 NREL 5 MW風力機主要參數信息表Table 2 Key parameters of the NREL 5 MW wind turbine

圖6為流場計算區域示意圖。在長、寬、高分別為17D、3D、3D的長方體計算區域中串列布置3臺水平軸NREL 5 MW風力機，沿流向分別記為WT-1（第一臺風力機）、WT-2（第二臺風力機）、WT-3（第三臺風力機）。網格在風力機附近區域進行加密，加密區的長、寬、高分別為0.5D、1.5D、1.5D，加密區網格單元尺寸為2.04 m。為了捕捉到渦系的細微結構，在風力機尾流區采用全加密方式，計算總網格量為1100萬。風力機葉片、機艙和塔筒均用體積力源項表示，其中葉片的分段數為30，塔筒的分段數為45。三臺風力機中，WT-1距離入口2D，WT-3距離出口5D，每臺風力機之間的距離固定為5D。

圖6 計算域示意圖Fig. 6 Schematic of the computational domain

2.2.2 非定常陣風計算工況

計算采用如圖3所示的來流陣風工況，正常平均風速為11.4 m/s，陣風周期為10.5 s，計算時間段分為兩部分：第一部分仿真時間0～220 s，參考風速為輪轂風速11.4 m/s；第二部分仿真時間220 ～500 s，入口邊界條件為陣風條件。選取五個時間點所對應的縱向瞬時速度云圖（圖7），分別為：陣風來流前215 s、陣風遇到風輪的時刻282.5 s（bg點）、陣風遇到風輪的時刻285 s（cg點）、陣風遇到風輪的時刻287.5 s（dg點）。

圖7 不同陣風時刻的縱向中心處的速度截面云圖Fig. 7 Velocity contours in the slice along the longitudinal center at different time instances under the gust wind condition

從圖7可以觀察到，隨著軸向位置增加，風力機尾流區也隨之膨脹，并且風速呈現出時間和空間的瞬時變化的特性，尾流膨脹邊界在向下游發展的過程中變得越來越不清晰。由于受到風剪切的影響，整個尾流區中在縱向上出現了靠近地面的流速小于輪轂以上的流速，尤其是在第一臺到第二臺風力機之間表現得更為明顯。受下游風力機的影響，第一臺風力機的尾流來不及恢復就再次發生流速的降低，從而整個尾流區只有最后的區域出現了流速的恢復現象。此外，對比不同時刻的陣風速度云圖，可以發現，當陣風處于高點cg時，尾流膨脹較正常來流時更加明顯。

圖8為陣風經過時的尾流流場渦量云圖，采用的是Q準則，Q值取0.0005。圖8（a）為三臺風力機在沒有陣風入流工況下的尾流渦系結構圖，可以看出，葉尖渦卷起并逐漸向下游發展。近尾流的葉尖渦系結構明顯。隨著渦系的發展，出現了渦混合和渦破碎現象，特別是在第二臺到第三臺風力機位置，渦已經完全混合，并且隨著尾流中的湍流強度的增加，尾流中的渦系受到風剪切的影響，導致上下渦傾斜往下游發展。隨著時間的推移，陣風入流開始影響尾流渦系結構，圖8（b）可見葉尖渦間距開始發生改變，圖8（c）中葉尖渦系開始變大并向下游移動，并且導致尾流區出現更大的膨脹；圖8（d）陣風渦系也同樣發生了渦破碎和渦融合現象，同時產生尾流速度瞬時變化和湍流強度變化。對比圖8（a）可以發現，陣風渦系的葉尖渦渦核更大，并且葉尖渦之間的間距更寬。

圖8 不同陣風時刻的由速度染色的渦等值面（Q = 0.0005 ）Fig. 8 Vortex iso-surface (Q = 0.000 5) colored by the velocity at different time instances under the gust wind condition

圖9是遭遇陣風之后三臺風力機的輸出功率及風輪推力對比圖，數值模擬中假定風力機來不及變槳。由圖可見，受陣風影響，功率及軸向推力較平穩周期的響應值出現了較大幅值的周期性波動，尤其是第一臺風力機的功率及軸向推力曲線出現了明顯的一個峰值和兩個谷值。對于第二臺風力機，由于上下游風力機之間存在尾流效應和尾流干擾，因此陣風入流后，第二臺風力機氣動載荷（功率和推力曲線）出現了兩個比較明顯的峰值突變，此時陣風的響應時間在第335 s左右，陣風響應結束時間大概在第390 s，是入流陣風周期的5倍左右，可見，陣風速度和尾流效應的耦合作用導致了第二臺風力機的陣風氣動載荷作用周期延長，同時曲線會出現兩個峰值，這個結論對于評估風力機載荷性能具有重要參考意義。第三臺風力機的輸出功率相對于第二臺出現了回升，這是因為，隨著尾流的發展，尾流與外界的空氣混合及能量交換加快，從而導致風速恢復速度變快。對比不同風力機的功率和推力時間曲線圖，可看出，第三臺風力機的氣動載荷波動幅度較大，這是由于風力機尾流發展帶來的湍流強度的增加，使得第三臺風力機來流速度的湍流強度增加，同時當陣風傳遞到第三臺風力機時，風力機對陣風的響應周期繼續拉長。由表3可見，第三臺風力機發電功率值回升，在有/無陣風的情況下，第一臺風力機承受的推力載荷最大，第二臺和第三臺風力機的推力載荷值分別增加9.1%和21%。

表3 有無陣風功率和推力的RMS與STD標準方差Table 3 RMS and STD values of the power and thrust under the gust and non-gust wind conditions

圖9 三臺風力機輸出功率和軸向推力對比圖Fig. 9 Comparison of the power and thrust of the three wind turbines

3 結論

本論文建立了 LES 耦合致動線模型的方法，并加入非定常修正，考慮了風力機機艙、塔架和風剪切等因素影響，實現了風力機非定常工況的數值模擬仿真，開展了陣風入流工況下的多臺風力機尾流的渦系結構演變規律和氣動載荷變化規律研究。得到結論如下：

1）基于Nibe風力機實驗數據，驗證了 LES 耦合致動線模型方法的可行性及有效性，模擬結果也表明，該方法在近尾流區和遠尾流區都可以獲得較好的尾流速度計算結果。

2）隨著軸向位置的增加，風力機尾流區也隨之膨脹，尾流速度出現先降低再恢復的發展特性，尾流區膨脹邊界在向下游發展的過程中變得越來越不清晰。陣風使得尾流區會比正常來流時尾流膨脹得更加明顯，即尾流周向直徑變大。

3）成功捕捉到了陣風渦系結構。陣風會導致葉尖渦的渦間距變大，并且在逐漸向下游發展的過程中葉尖渦核也變得更大，渦破碎和渦混合更易發生。

4）陣風入流對第一臺風力機的影響比較明顯。而第二臺風力機在尾流效應和陣風耦合的影響下，氣動載荷曲線出現了兩個峰值，下游陣風響應周期增加為入流陣風周期的5倍。雖然第一臺風力機推力載荷最大，但是第二臺和第三臺風力機在陣風的影響下，推力載荷也分別增加了9.1%和21%。

5）本文研究內容及結論對風電場布局和風力機載荷強度設計具有指導意義。