曹九發(fā)，宋佺珉，王超群，朱衛(wèi)軍，柯世堂

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

0 引言

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

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

1 數(shù)值模型與方法

1.1 N-S/LES模型

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

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

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

τ由下式給出：

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

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

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

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

1.2 致動線模型AL

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

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

葉片翼型截面的當?shù)厮俣萚15]為：

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

翼型截面當?shù)厮俣扰c風輪平面夾角為入流角，

翼型截面的當?shù)毓ソ菫槿肓鹘菧p去槳距角γ，

確定了葉片翼型的當?shù)厮俣扰c攻角后，即可求出葉片單位展長的升力與阻力：

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

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

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

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

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

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

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

1.3 陣風

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

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

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

其中：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點為陣風周期內(nèi)風速最低的時間點位置，cg點為陣風周期內(nèi)風速最高的時間點位置，eg點表示陣風周期結(jié)束后風速回歸到正常風速值的位置。

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

2 結(jié)果與分析

2.1 算例驗證

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

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

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

圖4 計算域網(wǎng)格示意圖Fig. 4 Schematic of the computational domain and grid

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

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

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

2.2 NREL 5 MW陣風工況數(shù)值計算

2.2.1 邊界條件和計算域

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

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

圖6為流場計算區(qū)域示意圖。在長、寬、高分別為17D、3D、3D的長方體計算區(qū)域中串列布置3臺水平軸NREL 5 MW風力機，沿流向分別記為WT-1（第一臺風力機）、WT-2（第二臺風力機）、WT-3（第三臺風力機）。網(wǎng)格在風力機附近區(qū)域進行加密，加密區(qū)的長、寬、高分別為0.5D、1.5D、1.5D，加密區(qū)網(wǎng)格單元尺寸為2.04 m。為了捕捉到渦系的細微結(jié)構(gòu)，在風力機尾流區(qū)采用全加密方式，計算總網(wǎng)格量為1100萬。風力機葉片、機艙和塔筒均用體積力源項表示，其中葉片的分段數(shù)為30，塔筒的分段數(shù)為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可以觀察到，隨著軸向位置增加，風力機尾流區(qū)也隨之膨脹，并且風速呈現(xiàn)出時間和空間的瞬時變化的特性，尾流膨脹邊界在向下游發(fā)展的過程中變得越來越不清晰。由于受到風剪切的影響，整個尾流區(qū)中在縱向上出現(xiàn)了靠近地面的流速小于輪轂以上的流速，尤其是在第一臺到第二臺風力機之間表現(xiàn)得更為明顯。受下游風力機的影響，第一臺風力機的尾流來不及恢復就再次發(fā)生流速的降低，從而整個尾流區(qū)只有最后的區(qū)域出現(xiàn)了流速的恢復現(xiàn)象。此外，對比不同時刻的陣風速度云圖，可以發(fā)現(xiàn)，當陣風處于高點cg時，尾流膨脹較正常來流時更加明顯。

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

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

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

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

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

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

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

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