考慮大氣穩(wěn)定度的風電場等效粗糙度模型研究

摘要：目前模型未充分考慮大氣穩(wěn)定度的影響，針對該問題，在不同大氣穩(wěn)定度下，通過建立風電場邊界層的動量平衡關(guān)系和揭示風電場的流動不均勻性，提出一種適用于不同大氣穩(wěn)定度的風電場等效粗糙度模型。采用大渦模擬方法對所提模型進行驗證，結(jié)果顯示該方法能有效評估大氣穩(wěn)定度對風電場流動不均勻性以及等效粗糙度的影響，所得等效粗糙度平均誤差約10%。

關(guān)鍵詞：風電場；大氣邊界層；尾流；大渦數(shù)值模擬；等效粗糙度模型

中圖分類號：TM614文獻標志碼：A

0引言

隨著風電場建設(shè)規(guī)模的擴大，大型風電場尾流可顯著改變大氣邊界層的速度剖面[1、湍流度2等物理統(tǒng)計量。在風電場中尺度模擬中，通常將風電場簡化為具有附加粗糙度的下墊面。粗糙度是反映地面遮擋物平均粗糙程度的指標，其大小等于風速對數(shù)廓線中速度為零的高度。在風電場簡化處理中，可將風電場視為一種特殊的“粗糙元”陣列，其等效粗糙度的大小取決于機組運行狀態(tài)、機組幾何尺寸、風電場疏密程度、地面粗糙度以及大氣穩(wěn)定度等參量。

Lettau3首次提出僅考慮遮擋物幾何尺寸的等效粗糙度模型，開啟了粗糙度模型的量化研究，并被Barrie等4用于評估大型風電場對全球天氣影響的研究中；Frandsen？以輪轂高度為界，將風電場邊界層分為上下兩層，假定兩層均符合風速對數(shù)廓線，并通過推力系數(shù)將機組運行狀態(tài)引入到等效粗糙度模型中；Calaf等1基于大渦數(shù)值模擬（large eddysimulation，LES）揭示了充分發(fā)展風電場的三應力層結(jié)構(gòu)，進一步提高了該類模型的預測精度。Stevens等6將尾流模型與粗糙度模型進行耦合，所提耦合尾流邊界層模型與大渦模擬結(jié)果吻合較好。張歡等7參考楊曉雷等I？對風電場不均勻度的定義，將不均勻度引入到Calaf模型1中，提高了粗糙度模型對不同排布間距風電場的普適性。

但上述模型均未考慮大氣穩(wěn)定度的影響，無法適用于不同大氣穩(wěn)定狀態(tài)下風電場等效粗糙度的預測。通常情況下，大氣邊界層在強迫因素（如地面熱通量）的作用下會呈現(xiàn)出穩(wěn)定、中性和不穩(wěn)定3種大氣熱穩(wěn)定度交替變化的日循環(huán)9。并且由于海、陸下墊面性質(zhì)的差異，長期大氣穩(wěn)定度也會在不同趨向上輕微偏離中性大氣[10-11]。因此，部分學者嘗試在粗糙度模型中引入穩(wěn)定度修正，以提高模型的適用范圍，如Pena等12在Emeis等[13所提兩應力層的基礎(chǔ)上單向引入穩(wěn)定度修正，但忽略了風電場應力層變化對大氣穩(wěn)定度的反作用；Sescu等14]在三應力層的基礎(chǔ)上解決了風電場與大氣穩(wěn)定度的耦合問題，但未考慮風電場排布間距的影響，普適性較低。

準確預測不同大氣穩(wěn)定度下的風電場等效粗糙度對大型風電場中尺度數(shù)值模擬、功率預測以及效能評估均具有較高的實用價值，但當前風電場等效粗糙度計算方法均存在一定缺陷。因此，本文綜合考慮大氣穩(wěn)定度和風電場流動不均勻性，提出一種普適性的大型風電場等效粗糙度計算模型，并使用風電場大渦數(shù)值模擬驗證所提模型的適用性。

1模型推導

大氣熱力分層狀態(tài)是影響大型風電場等效粗糙度的重要因素，兩者相互耦合。一方面，大型風電場對氣流的阻滯作用增強了大氣邊界層向下的動量輸運，打破了原本平衡的動量關(guān)系，進而改變了邊界層熱穩(wěn)定性；另一方面，大氣熱力分層狀態(tài)可增強或抑制流體微團的垂向運動能力，進而影響風電機組的尾流演化和風電場等效粗糙度。本節(jié)將基于兩者的耦合關(guān)系，建立普適性的風電場等效粗糙度求解模型。

1.1大氣穩(wěn)定度

對于大型風電場，當其流向距離超過邊界層高度一個數(shù)量級時（10 km），風電場流動將進入充分發(fā)展狀態(tài)，此時可近似為理想的無限大風電場。該理想狀態(tài)為模型的建立提供了有利條件。根據(jù)Calaf等1結(jié)果，如圖1所示，中性大氣條件下，充分發(fā)展風電場自下而上可分為內(nèi)層、尾流層和外層3個應力層，即三應力層結(jié)構(gòu)。由應力平衡可知，外層的等效摩擦應力應等于內(nèi)層的摩擦應力與尾流層的阻滯應力之和，即：

式中：p——空氣密度，kg/m3;u——風電場摩擦速度，m/s;u——地面摩擦速度，m/s;Ca——風電機組摩擦系數(shù)，Cn=πGr/（4S.S，），其中C為風電機組的推力系數(shù)，S。和S，——風電場中相鄰風電機組的無量納流向間距和展向間距（見圖2）;α——輪轂高度風速的不均勻度；（u）——風電場輪轂高度處的時空平均風速，m/s。

對于一般性的大氣邊界層，依照莫寧-奧布霍夫相似理論，考慮大氣穩(wěn)定度修正的風速梯度為：

式中：K——馮·卡門常數(shù)，k=0.4;z——距離地面的高度，m;u.——摩擦速度，m/s;Kzu.——湍流粘度，可合寫為v，m2/s;（u〉——某高度處平面的時空平均風速；m/s;φ（z/L）——大氣穩(wěn)定度的動量修正系數(shù)；L——莫寧-奧布霍夫長度，m。

風電場邊界層的內(nèi)層和外層仍符合一般風電場的物理特性，滿足式（2）。將式（2）從z。到z進行定積分，并忽略小項可分別得到內(nèi)層和外層的穩(wěn)定度修正的風速對數(shù)廓線，即：

式中：201?！孛娲植诙?，m;20，i——風電場等效粗糙度，m;L。和L——輪轂高度以下和以上的莫寧-奧布霍夫長度，m，L？=-u30。/（kgq），L=-u30。/（kgq），其中g(shù)為重力加速度（取9.81 m/s2），θ。為地面參考位溫（K）;q，為地面熱流通量（K·m/s）。

本文將尾流層范圍取做（z-3D/4，zn+D/R），其中z為輪轂高度，m，穩(wěn)定大氣和中性大氣下R=2，不穩(wěn)定大氣下R=4。尾流層以下即為內(nèi)層，尾流層以上即為外層。大氣穩(wěn)定度動量修正項ψa（z/L）=J[[1-φn（z/L）]dIn（z/L），是根據(jù)實際測量數(shù)據(jù)歸納的經(jīng)驗函數(shù)，根據(jù)文獻[15]，大氣穩(wěn)定度動量修正項ψm（z/L）的具體表達式為：

式中：am=6.1;b=2.5;a。=0.33;b？=0.41;c=（-2（a.L）;ψo=-Inaa+32a1bπ/6。

多臺風電機組的尾流相互摻混使尾流層的平均風速減小、湍流動能增大，因此反映動量交換劇烈程度的湍流黏度vi需增加額外的附加渦黏系數(shù)v，此時vr=（kzu。+v）。因此在尾流層中，式（2）可改寫為：

上述尾流層的湍流動能增大是由風電機組吸收的動量導致的，而風電機組吸收的來流動量正比于，由此估算湍流的速度尺度為α（u）JCa/2。長度尺度可為風輪直徑D，因此附加渦粘系數(shù)v。可估算為α（u）DJC？/2。

對式（6）進行不定積分得：

式中：vw=v./（kznu.）;C——待確定的積分常數(shù)。

在尾流層與內(nèi)層的交界處風速應相等，據(jù)此聯(lián)立式（7）與式（3）可得：

同理，聯(lián)立式（7）與式（4）可得：

輪轂高度風速應保持連續(xù)，因此在z=z處聯(lián)立式（8）和式（9）可得：

式中：

將式（10）代入式（1）中，即可得到適用不同大氣穩(wěn)定度的風電場等效粗糙度解析模型：

在α已知的前提下，式（11）中仍包含3個獨立變量u、u和z？，但該模型推導過程中僅包含式（1）和式（10）兩個獨立的方程。因此，需引入額外的獨立方程使模型封閉。在風電場上游，來流的風速廓線可描述為：

式中：u——來流的摩擦速度，m/s;L。——對應的大氣穩(wěn)定度，m。

該模型在未來實際工程應用中，上游大氣邊界層穩(wěn)定度參數(shù)可通過測量至少兩個高度處的風速和溫度數(shù)據(jù)161，利用最小二乘法擬合式（12）得到。假定在風電場邊界層頂部δ處風速保持不變，則聯(lián)立式（4）和式（12）可得：

至此，可通過聯(lián)立式（1）、式（10）和式（13）進行求解。

1.2入流不均勻度

文獻[7-8]認為，在充分發(fā)展區(qū)域，機組的入流風速ur并不等于風電場輪轂高度處的時空平均風速（u）。為提高模型的普適性，適應不同排布間距的風電場，張歡等7引入了輪轂高度風速的不均勻度α，定義為：

類似地，本文認為當風電場處于不同大氣穩(wěn)定度時，其尾流膨脹率的變化也會導致不均勻度的變化。因此，本文將不均勻度引入到所提模型中。與張歡等7不同的是，本文采用Frandsen等17基于動量守恒推導的一維尾流模型（Frandsen尾流模型示意圖如圖3所示）來評估流場的不均勻性。該尾流模型所描述的機組下游x位置處的尾流速度u。為：

式中：u?！獊砹黠L速，m/s;β——尾流擴張后的初始截面積與風輪掃掠面積之比，β=（1-a）/（1-2a），其中a為軸向誘導因子；k——尾流膨脹率，用以表征尾流直徑d.隨下游距離x變化的快慢，是一個待定參數(shù)。

采用平方和疊加方法計算多臺機組尾流損失，假定流場中某位置處，其上游所有的機組集合為A，則該位置處的風速u可計算為：

在輪轂高度平面，將機組的投影線段L（圖2中表示風輪的線段）均勻地劃分為N，個線段，分辨率dy=6m，計算所有子線段中點的平均風速作為機組的入流風速ur，計算式為：

式中：u;——第i個線段中點的風速，m/s。

在輪轂高度平面，將每臺風電機組所占區(qū)域S（圖2中陰影區(qū)域）均勻地劃分為N，個網(wǎng)格，分辨率dx=dy=6m，計算所有子網(wǎng)格中點的平均風速作為充分發(fā)展區(qū)域輪轂高度的平均風速（ū），計算式為：

式中：u;——第j個網(wǎng)格中點位置處的風速，m/s。

至此，利用式（17）和式（18）即可計算風電場入流不均勻度α。

1.3模型耦合求解

從1.2節(jié)和1.3節(jié)的推導過程可知，風電場等效粗糙度模型中含有待定參數(shù)不均勻度α，尾流模型中含有待定參數(shù)尾流膨脹率k?？衫幂嗇灨叨忍幤骄L速相等的條件將兩者耦合，迭代求解得到不均勻度α和尾流膨脹率k。，即可獲得適應不同大氣穩(wěn)定度的風電場等效粗糙度。具體計算流程如下：

1）利用式（17）、式（18）計算得到特定排布下，一系列尾流膨脹系數(shù)k所對應的輪轂高度入流風速ur/u。和平面平均風速（u）lu。。計算時，為保證風電場充分發(fā)展，按照原有排布間距將風電機組上游擴充100行，兩側(cè)各擴充4列。利用式（14）即可得到（u）/u=f（k.）和α=h（lk.）兩條對應關(guān)系曲線。

2）給定尾流膨脹系數(shù)的初值k。

3）根據(jù)α=h（k.）的關(guān)系曲線，得到初始的尾流膨脹系數(shù)相應的初始不均勻度α。

4）將得到的不均勻度α代入所述等效粗糙度模型，求解此時的風電場等效粗糙度。

5）利用式（12）和式（9）分別計算輪轂高度的來流風速和充分發(fā)展區(qū)域的平面平均風速，兩者比值為（u）lu。

6）從式（15）可知，（π）lu。=f（2.）為單調(diào)遞增關(guān)系，因此可將步驟5）所得（z）u.代入反函數(shù)k，=9（u）tu），得到更新的尾流膨脹系數(shù)k。

7）根據(jù)步驟6）所得k，更新相應的不均勻度α。

8）重復步驟4）～7）直至不均勻度α變化小于1%，最終得到風電場等效粗糙度。

2模型驗證

2.1大渦數(shù)值模擬方法及算例設(shè)置

本文采用開源的大渦模擬求解器LESGO進行風電場數(shù)值模擬并用以驗證模型準確性。

2.1.1大渦數(shù)值模擬方法

風電場大氣邊界層的大渦數(shù)值模擬需求解無量綱的空間濾波不可壓縮N-S方程、連續(xù)性方程和位溫的對流擴散方程。高雷諾數(shù)下忽略粘性項的影響，因此控制方程組為：

式中：？;——可解尺度的速度；上標～——網(wǎng)格尺度△的空間濾波；p——修正后的可解尺度壓強，p=p+pT從/3;r——亞格子偏應力，采用拉格朗日動態(tài)亞格子模型進行封閉[183，r=T，-δ，r/3;，下標i或j可取1、2、3，分別代表流向x、展向y和垂向z;0——可解尺度的位溫；〈·——水平方向的空間平均；δ？——Kronecker算子；f——風電機組引發(fā)的體積力，采用修正致動盤模型，詳見文獻[20];F，——流向的強制壓力梯度，F(xiàn)，=u2/8，其中本文取u=0.45，δ為邊界層高度，m;g;——亞格子熱流通量，依據(jù)文獻[19]進行建模。

在方程求解中，在計算域的水平方向上采用均勻分布的網(wǎng)格，采用偽譜法進行方程的離散21];在垂直方向上采用交錯網(wǎng)格；采用二階中心差分方法進行離散；時間推進采用二階精度的Adams-Bashforth方法。在邊界條件設(shè)置中，頂部邊界條件采用無應力邊界條件；水平方向采用周期性邊界條件，用以模擬充分發(fā)展風電場；底部采用基于莫寧-奧布霍夫相似理論的剛性壁面模型[22]。

2.1.2算例設(shè)置

共設(shè)置3類大氣穩(wěn)定度的6個算例，通過調(diào)節(jié)地面熱通量模擬具有不同熱穩(wěn)定性的大氣邊界層。地面粗糙度取0.1 m;風電機組的輪轂高度和風輪直徑均設(shè)置為100m;流向間距SD和展向間距S，D分別為πD/4和πD/6;風輪的推力系數(shù)取0.75;穩(wěn)定大氣（SBL）、中性大氣（NBL）和不穩(wěn)定大氣（CBL）3種條件下，風電場邊界層高度δ分別取600、850和1100 m。不同大氣穩(wěn)定度大渦模擬算例設(shè)置見表1。風電場排布示意圖如圖2所示。

2.2結(jié)果分析和對比

圖4為輪轂高度不同大氣穩(wěn)定度下的時均風速等值線圖。由圖4可見，由于浮力的差異造成了尾流摻混的程度不同，不同大氣穩(wěn)定度下速度分布的均勻程度不同。對于穩(wěn)定大氣條件，兩列風電機組之間形成了氣流的高速通道，導致風電機組的入流風速ur低于平面平均風速（u），因此αlt;1;不穩(wěn)定大氣則相反，風電機組的入流風速ur略高于平面平均風速（u）;中性大氣則介于兩者之間。圖5為所提模型和LES所得不均勻度對比。由圖5可見，本文模型能較好地捕獲不同大氣穩(wěn)定度下不均勻度α的變化趨勢，使模型具有更高的精度。

圖6為粗糙度模型和LES所得等效粗糙度的對比。LES所得等效粗糙度是根據(jù)2z，位置處的平均風速，由式（4）求解得到的。各算例的不均勻度α以及等效粗糙度詳見表2?？梢园l(fā)現(xiàn)，未考慮大氣穩(wěn)定度的等效粗糙度模型在非中性條件下具有較大誤差，所提模型能夠較好地捕獲風電場等效粗糙度隨大氣穩(wěn)定度的變化趨勢，模型計算結(jié)果與大渦模擬結(jié)果符合較好，平均誤差約為10%。

3結(jié)論

隨著風電場建設(shè)規(guī)模的擴大，大型風電場尾流可顯著改變邊界層的物理統(tǒng)計量，影響其自身以及下游風電場的發(fā)電能力。在風電場中尺度模擬中，風電場等效粗糙度常被用來表征大型風電基地，但當前模型未能充分考慮大氣穩(wěn)定度的影響，無法適用于不同大氣穩(wěn)定狀態(tài)下風電場等效粗糙度的預測。針對該問題，本文綜合考慮風電場邊界層的動量平衡關(guān)系和風電場的流動不均勻性，提出一種適用于不同大氣穩(wěn)定度的風電場等效粗糙度模型。采用大渦模擬方法對所提模型進行驗證。結(jié)果顯示，所提等效粗糙度模型可有效捕獲不均勻度隨大氣穩(wěn)定度的變化趨勢，與大渦模擬所得不均勻度具有較高的一致性。同時，較好地預測了不同大氣穩(wěn)定度下的風電場等效粗糙度，平均誤差約為10%。該模型可快速、準確地計算出不同大氣穩(wěn)定度下的風電場等效粗糙度，有望對風電場中尺度模擬、大型風電場尾流效應評估等方面提供重要依據(jù)。

[參考文獻]

[1]CALAF M，MENEVEAU C，MEYERS J.Largeeddysimulation study of fully developed wind-turbine arrayboundary layers [J].Physics of fluids，2010，22（1）：015110.

[2]CALAF M，PARLANGE M B，MENEVEAU C.Largeeddy simulation study of scalar transport in fully developedwind-turbine array boundary layers [J].Physics of fluids，2011，23（12）：126603.

[3]LETTAUH.Note on aerodynamic roughness-parameterestimation on the basis of roughness-element description[J].Journal of applied meteorology，1969，8（5）：828-832.

[4]BARRIE D B，KIRK-DAVIDOFF D B.Weatherresponseto management of a large wind turbine array [J].Atmospheric chemistry and physics discussions，2009，9（1）：2917-2931.

[5]FRANDSENS.On the wind speed reduction in the centerof large clusters of wind turbines [J].Journal of windengineering and industrial aerodynamics，1992，39（1/2/3）：251-265.

[6]STEVENS R J A M，GAYME D F，MENEVEAU C.Coupled wake boundary layer model of wind-farms [J].Journal of renewable and sustainable energy，2015，7（2）：023115.

[7]ZHANG H，GE M W，LIU Y Q，etal.A new coupledmodel for the equivalent roughness heights of wind farms[J].Renewable energy，2021，171：34-46.

[8]YANGXL，KANGS，SOTIROPOULOS F.Computationalstudy and modeling of turbine spacing effects infinitealigned wind farms[J].Physics of fluids，2012，24（11）：115107.

[9]ABKAR M，SHARIFI A，PORTé-AGEL F.Wake flow ina wind farm during a diurnal cycle[J].Journal ofturbulence，2016，17（4）：420-441.

[10]PENA A，HAHMANN A N.Atmospheric stability andturbulence fluxes at Horns Rev—an intercomparison ofsonic，bulk and WRF model data[J].Wind energy，2012，15（5）：717-731.

[11]TROEN I，LUNDTANG PETERSEN E.Europeanwindatlas[M].Roskilde：Riso National Laboratory，1989.

[12]PENA A，RATHMANN 0.Atmospheric stability-dependent infinite wind-farm models and the wake-decaycoefficient[J].Wind energy，2014，17（8）：1269-1285.

[13]EMEIS S，F(xiàn)RANDSEN S.Reduction of horizontal windspeed in a boundary layer with obstacles [J].Boundary-layer meteorology，1993，64（3）：297-305.

[14]SESCU A，MENEVEAU C.Large-eddy simulation andsingle-column modeling of thermally stratified wind turbinearrays for fully developed，stationaryatmosphericconditions [J].Journal of atmospheric and oceanictechnology，2015，32（6）：1144-1162.

[15]BRUTSAERTW.Hydrology：an introduction[M].Cambridge：CambridgeUniveristy Press，2005.

[16]劉永前，馬曉梅，高小力，等.基于中性等效風速的海上風廓線建模方法研究[J].華北電力大學學報（自然科學版），2021，48（6）：97-105.

LIU YQ，MAXM，GAOXL，etal.Study on wind profilemodeling of offshore based on neutral equivalent windspeed [J].Journal of North China Electric PowerUniversity（natural science edition），2021，48（6）：97-105.

[17]FRANDSEN S，BARTHELMIE R，PRYOR S，etal.Analytical modelling of wind speed deficit in large offshorewind farms[J].Wind energy，2006，9（1-2）：39-53.

[18]BOU-ZEID E，MENEVEAU C，PARLANGE M.A scale-dependent Lagrangian dynamic model for large eddysimulation of complex turbulent flows [J].Physics offluids，2005，17（2）：025105.

[19]DU B W，GE M W，ZENG C J，etal.Influenceofatmospheric stability on wind-turbine wakes with a certainhub-height turbulence intensity[J].Physics of fluids，2021，33（5）：055111.

[20]SHAPIRO C R，GAYME DF，MENEVEAU C.Filteredactuatordisks：theory and application to wind turbinemodels in large eddy simulation[J].Wind energy，2019，22（10）：1414-1420.

[21]MOENG C H.A large-eddy-simulation model for the studyof planetary boundary-layer turbulence[J].Journal of theatmospheric sciences，1984，41（13）：2052-2062.

[22]MONIN A S，OBUKHOV A M.Basic laws of turbulentmixing in the surface layer of the atmosphere [EB/OL].https：//www.mcnaughty.com/keith/papers/Monin_and_Obukhov_1954.pdf

STUDY ON EFFECTIVE ROUGHNESS MODEL FOR WIND FARMSCONSIDERING ATMOSPHERIC STABILITY

Li Baoliang1，Zhang Ziliang2，Ge Mingwei1，Wang Luo2，Liu Yongqian1

（1.State Key Laboratory for Alternate Electrical Power System with Renewable Energy Sources，North China Electric Power University，Beijing 102206，China;2.Science and Technology Research Institute，China Three Gorges Corporation，Beijing 101199，China）

Abstract：Current wind farm effective roughness models do not fully consider the influence of atmospheric stabilities.To address thisproblem，an effective roughness model for wind farms with different atmospheric stabilities is proposed by establishing the momentumbalance relationship in the wind farm boundary layer and revealing the flow inhomogeneity of wind farms under different atmosphericstabilities.The proposed model is validated by using the large eddy simulation method，and the results show that the proposed model canffectively evaluate the effect of atmospheric stabilities on the flow inhomogeneity and the wind farm effective roughness，andtheaverage error of the obtained effective roughness is about 10%.

Keywords：windfarm;atmospheric boundary layer;wakes;large eddy simulation;effective roughness model