復(fù)雜地形風(fēng)電場(chǎng)尾流影響規(guī)律研究

王 欣 李 超 劉偉江

（1.浙江運(yùn)達(dá)風(fēng)電股份有限公司技術(shù)中心，浙江 杭州 310012；2.浙江省風(fēng)力發(fā)電技術(shù)重點(diǎn)實(shí)驗(yàn)室，浙江 杭州 311100；3.水利部產(chǎn)品質(zhì)量標(biāo)準(zhǔn)研究所 規(guī)劃咨詢處，浙江 杭州 310024）

0 引言

復(fù)雜地形風(fēng)電場(chǎng)尾流效應(yīng)仿真的準(zhǔn)確性直接影響風(fēng)電場(chǎng)微觀選址的可靠性。針對(duì)風(fēng)電場(chǎng)尾流影響分析，丹麥學(xué)者M(jìn)IKKELSEN基于BEM（blade element momentum）理論提出制動(dòng)模型，與線性疊加理論相比，該模型提高了尾流計(jì)算的精度，國(guó)內(nèi)外學(xué)者基于該思路延伸出制動(dòng)盤、制動(dòng)線以及制動(dòng)面模型。歐洲盲測(cè)報(bào)告曾指出，與只考慮湍動(dòng)能的一方程模型相比，考慮湍動(dòng)能和湍流耗散率的兩方程模型可以更準(zhǔn)確地捕捉大氣邊界層流動(dòng)。WAKES等人以中尺度或二維模擬結(jié)果作為三維流體域的邊界條件，該方法能使入口邊界來流風(fēng)況得到充分發(fā)展，增大了計(jì)算域。大氣邊界層流動(dòng)仿真是復(fù)雜地形流動(dòng)仿真的難點(diǎn)，國(guó)內(nèi)外學(xué)者對(duì)該問題進(jìn)行了大量研究，并基于CFD提出了精細(xì)化仿真方法。

該文結(jié)合三維Navier-Stokes（N-S）方程和制動(dòng)盤模型，并對(duì)大氣熱穩(wěn)定度的影響進(jìn)行計(jì)算，從而建立了可精確模擬包括風(fēng)力機(jī)尾流和大氣邊界層流動(dòng)的復(fù)雜地形風(fēng)電場(chǎng)尾流效應(yīng)計(jì)算的方法。同時(shí)，以某實(shí)際復(fù)雜地形風(fēng)電場(chǎng)作為分析對(duì)象，對(duì)受尾流影響的風(fēng)力機(jī)前后流場(chǎng)特性進(jìn)行數(shù)值計(jì)算，通過與實(shí)際運(yùn)行數(shù)據(jù)進(jìn)行對(duì)比，驗(yàn)證了該文所建立的計(jì)算方法的準(zhǔn)確性。

1 計(jì)算方法

1.1 制動(dòng)盤模型

假設(shè)大氣為單相干空氣，不考慮大氣中的塵粒和水汽。因此，該仿真求解定常、單相和三維不可壓縮雷諾時(shí)均N-S方程（簡(jiǎn)稱為RANS）。氣流經(jīng)過葉輪時(shí)受軸向的阻力和切向的誘導(dǎo)力作用，在該計(jì)算所采用的制動(dòng)盤模型中忽略切向力和尾流旋轉(zhuǎn)效應(yīng)，采用軸向阻力形成的壓差來描述葉輪在流場(chǎng)中的作用。

葉輪前后壓差如公式（1）所示。

式中：Δ為葉輪前后壓差；為葉輪后壓力；為葉輪前壓力；為空氣密度；為來流速度；為葉輪后速度。

經(jīng)過轉(zhuǎn)換和推導(dǎo)，公式（1）可表述為公式（2）。

式中：C為推力系數(shù)。

可以根據(jù)來流風(fēng)速和推力曲線獲得C，在計(jì)算時(shí)不考慮風(fēng)剪切和葉片徑向氣動(dòng)變化對(duì)C的影響。該計(jì)算根據(jù)風(fēng)輪前端的風(fēng)速變化并結(jié)合推力曲線動(dòng)態(tài)計(jì)算體積力源項(xiàng)。在數(shù)值計(jì)算時(shí)，通過迭代求解添加了體積力源項(xiàng)的動(dòng)量方程，以實(shí)現(xiàn)葉輪對(duì)氣流的作用。

1.2 求解器

該計(jì)算以O(shè)penFOAM中內(nèi)置兩方程的simpleFoam求解器為基礎(chǔ)進(jìn)行二次開發(fā)，添加了基于Monin-Obukhov長(zhǎng)度的浮力相關(guān)計(jì)算。simpleFoam采用壓力速度耦合的Simple算法，求解思路如下：1） 根據(jù)初始?jí)毫η蠼怆x散動(dòng)量方程，獲得初始速度。2） 因?yàn)槌跏妓俣葻o法滿足連續(xù)方程，所以構(gòu)建壓力修正量。3） 修正后的壓力通過動(dòng)量方程獲得修正后的速度。重復(fù)上述過程，直到速度滿足連續(xù)方程就可以完成迭代。仿真采用有限體積法在空間上對(duì)控制方程進(jìn)行離散。在離散格式中，gradSchemes采用Gauss linear，divSchemes采用Gauss upwind，laplacianSchemes采用Guass linear corrected。

1.3 計(jì)算域及網(wǎng)格處理

計(jì)算域如圖1所示，以機(jī)組WT#8為圓心截取5km×5km的地形范圍，計(jì)算域高度為1.2 km。采用二維模擬結(jié)果作為三維流體域的邊界條件，使來流到達(dá)對(duì)象區(qū)域時(shí)已形成穩(wěn)定風(fēng)況。

圖1 計(jì)算域及機(jī)組排布

計(jì)算網(wǎng)格采用OpenFOAM軟件內(nèi)的blockMesh模塊創(chuàng)建。高度方向grading參數(shù)給定90，網(wǎng)格高度0.3 m，對(duì)制動(dòng)盤所在區(qū)域進(jìn)行3次加密處理，總網(wǎng)格數(shù)為1.5×10。

1.4 計(jì)算設(shè)置

來流入口邊界條件采用考慮地表粗糙度高度（）和Monin-Obukhov 長(zhǎng)度（）的修正對(duì)數(shù)風(fēng)廓線（），如公式（3）所示。

式中：為地表粗糙度長(zhǎng)度；為垂直高度；為地表摩擦速度；為von Karman數(shù)，在該計(jì)算中=0.41；為經(jīng)驗(yàn)公式，無實(shí)際物理意義。

入口邊界湍動(dòng)能及耗散率如公式（4）、公式（5）所示。

式中：C為經(jīng)驗(yàn)常數(shù)，C=0.09。

計(jì)算域頂部邊界條件選擇壓力0和速度-；計(jì)算域底部邊界選擇壁面函數(shù)；利用decomposePar模塊對(duì)計(jì)算域進(jìn)行分塊，分配到計(jì)算機(jī)的多核心同步求解。

2 實(shí)例驗(yàn)證

2.1 結(jié)果驗(yàn)證

對(duì)機(jī)艙風(fēng)速風(fēng)向儀記錄的數(shù)據(jù)進(jìn)行風(fēng)向切片統(tǒng)計(jì)，建立尾流影響扇區(qū)（343 °～353 °）WT#7和WT#8機(jī)組處的風(fēng)速散點(diǎn)統(tǒng)計(jì)關(guān)系，如圖2所示。采用WT#7機(jī)組機(jī)艙風(fēng)速作為WT#7機(jī)組來流風(fēng)速，WT#8機(jī)組機(jī)艙風(fēng)速作為WT#7機(jī)組尾流區(qū)風(fēng)速。

圖2 WT#7機(jī)組和WT#8機(jī)組處風(fēng)速比

散點(diǎn)擬合給出WT#7機(jī)組實(shí)測(cè)風(fēng)輪前后風(fēng)速關(guān)系式，2點(diǎn)的風(fēng)速比為0.723。基于仿真計(jì)算結(jié)果擬合2點(diǎn)的風(fēng)速比為0.718，絕對(duì)誤差為0.005。

機(jī)艙風(fēng)速測(cè)量點(diǎn)位于風(fēng)輪后方，受葉輪影響，該位置的風(fēng)速小于葉輪前端的風(fēng)速，同時(shí)存在數(shù)值誤差和散點(diǎn)擬合誤差。因此導(dǎo)致實(shí)測(cè)結(jié)果曲線和仿真結(jié)果曲線存在一定偏移。

2.2 大氣熱穩(wěn)定度影響分析

在風(fēng)電場(chǎng)流動(dòng)仿真分析過程中的大氣熱穩(wěn)定度往往被假設(shè)為中性，其代表的熱效應(yīng)被忽略。大氣熱穩(wěn)定度可分為穩(wěn)態(tài)（stable）、中性（neutral）以及非穩(wěn)態(tài)（unstable）3種狀態(tài)類型。在不同狀態(tài)下的地表與大氣溫差會(huì)導(dǎo)致大氣在垂直方向上存在不同的對(duì)流形態(tài)，該垂直方向上的作用力導(dǎo)致氣流貼地流動(dòng)發(fā)生變化，尤其在遇到障礙后。在低風(fēng)速?gòu)?fù)雜地形上，空氣流動(dòng)緩慢且受起伏地形影響易形成持續(xù)漩渦，垂直方向上的浮力作用影響將變得顯著。

不同熱穩(wěn)定度下大氣邊界層的流動(dòng)狀態(tài)如圖3所示。區(qū)域1為來流平坦區(qū)域，熱穩(wěn)定度未產(chǎn)生較明顯影響，在穩(wěn)態(tài)、中性以及非穩(wěn)態(tài)3種狀態(tài)下氣流呈相似的貼地流動(dòng)。區(qū)域2為氣流在起伏山地上流動(dòng)，此時(shí)氣流貼地流動(dòng)性變差，并伴有減速現(xiàn)象，低速區(qū)域影響范圍在長(zhǎng)度上呈現(xiàn)非穩(wěn)態(tài)＞中性＞穩(wěn)態(tài)的規(guī)律。區(qū)域3為背風(fēng)面平坦區(qū)域，氣流流經(jīng)山地后在該區(qū)域呈現(xiàn)明顯的分層現(xiàn)象，且向斜上方飛逸，低速區(qū)域影響范圍在高度上同樣呈現(xiàn)非穩(wěn)態(tài)＞中性＞穩(wěn)態(tài)的規(guī)律。

圖3 不同熱穩(wěn)定度下風(fēng)速分布圖

2.3 尾流影響分析

在中性情況下，尾流區(qū)風(fēng)速分布如圖4所示，在WT#7機(jī)位點(diǎn)的山坡后面出現(xiàn)了1個(gè)低風(fēng)速區(qū)。WT#8機(jī)組不僅受WT#7的尾流影響，而且還受該低風(fēng)速區(qū)的影響。

圖5對(duì)圖4中制動(dòng)盤前后0.8（為機(jī)組風(fēng)輪直徑）距離的4個(gè)位置（P、P、P以及P）進(jìn)行風(fēng)廓線取樣分析。由圖5可知，當(dāng)高度高于150 m時(shí)，速度出現(xiàn)階梯形不光滑分布，這是由該高度網(wǎng)格分辨率開始降低引起的。該高度已高于制動(dòng)盤高度（制動(dòng)盤相對(duì)高度為38 m～113 m），并不影響結(jié)果分析。

圖4 尾流影響下風(fēng)速分布（中性）

圖5（a）給出了來流風(fēng)況特性。隨著相對(duì)高度的增加，風(fēng)速快速增大，在高度約為10 m～11 m附近達(dá)到7.0 m/s～7.5 m/s。當(dāng)相對(duì)高度為15 m～200 m時(shí)，前方山峰影響帶的作用顯著，風(fēng)廓線風(fēng)速呈降低的狀態(tài)。隨著高度的進(jìn)一步增加，脫離了影響帶的作用，風(fēng)廓線風(fēng)速回升并穩(wěn)定在11.3 m/s。

由圖5（a）和圖5（c）可知，在添加制動(dòng)盤后，風(fēng)廓線在制動(dòng)盤前端呈現(xiàn)微小的風(fēng)速降低的現(xiàn)象，這與理論上風(fēng)輪前端誘導(dǎo)區(qū)出現(xiàn)風(fēng)速減弱的特性相符，說明了該制動(dòng)盤模型的準(zhǔn)確性。在復(fù)雜地形上，風(fēng)廓線并不按標(biāo)準(zhǔn)的對(duì)數(shù)正態(tài)分布。由于山坡對(duì)來流的阻擋，在坡后方貼近地表處形成低速區(qū)。由圖5（c）可知，在沒有制動(dòng)盤的情況下，P位置的風(fēng)廓線呈“S”形。

由圖5（b）～圖5（d）可知，在制動(dòng)盤對(duì)應(yīng)高度下，風(fēng)廓線出現(xiàn)明顯的風(fēng)速衰減段。當(dāng)復(fù)雜地形影響和尾流影響疊加時(shí)，風(fēng)廓線的“S”形分布更加明顯，且風(fēng)速衰減段在高度上的相對(duì)位置發(fā)生相應(yīng)變化。

圖5 制動(dòng)盤周圍風(fēng)廓線（中性）

3 結(jié)語

目前，國(guó)內(nèi)風(fēng)電行業(yè)存在復(fù)雜地形風(fēng)資源計(jì)算不準(zhǔn)確的問題，該文基于三維N-S方程和制動(dòng)盤模型建立了精細(xì)化尾流數(shù)值仿真方法，分析了在地形、尾流和大氣狀態(tài)共同作用下的流動(dòng)特性，并采用實(shí)際運(yùn)行數(shù)據(jù)驗(yàn)證了該方法的準(zhǔn)確性。

該文對(duì)大氣熱穩(wěn)定度的影響進(jìn)行了分析，說明了在復(fù)雜地形上近地面流動(dòng)存在分層現(xiàn)象，且與大氣熱穩(wěn)定度具有相關(guān)性。同時(shí)，由于受地形的影響，因此到達(dá)風(fēng)力機(jī)的風(fēng)廓線呈現(xiàn)出明顯的“S”形分布，與平坦地形上的標(biāo)準(zhǔn)對(duì)數(shù)正態(tài)分布有較大區(qū)別。