羅嘯宇， 李秋明， 劉震卿， 熊世樹

(1. 廣東電網(wǎng)有限責(zé)任公司電力科學(xué)研究院， 廣東 廣州 510062；2. 華中科技大學(xué) 土木工程與力學(xué)學(xué)院， 湖北 武漢 430074)

近幾年風(fēng)力發(fā)電發(fā)展迅速，已成為新能源領(lǐng)域研究熱點。而上游機組產(chǎn)生的尾流減速與湍流增大效應(yīng)將降低下游機組發(fā)電效率、增加疲勞載荷，嚴(yán)重影響下游機組運行性能。風(fēng)機尾流的研究方法主要有風(fēng)洞試驗[1～6]和計算流體力學(xué)(Computational Fluid Dynamics，CFD)[7～13]兩種方法。Alfredsson等[2]利用熱線儀進行了風(fēng)機尾流風(fēng)速分布規(guī)律研究，揭示了尾流效應(yīng)對風(fēng)機出力的影響。Ebert等[3]也利用此技術(shù)，進行了大量風(fēng)機尾流實驗，分析了葉尖渦與流速分布規(guī)律。但是，風(fēng)洞實驗無法獲得全場信息，因而，部分學(xué)者轉(zhuǎn)而采用CFD方法開展研究。閆海津等[8]采用幾何貼體網(wǎng)格，得到了全流域風(fēng)速與壓強分布、流動分離等結(jié)果，但未考慮葉片旋轉(zhuǎn)效應(yīng)、地面效應(yīng)和風(fēng)速梯度影響。李少華等[9]同樣采用幾何貼體網(wǎng)格，研究了風(fēng)力機組在串并列情況下尾流分布特點。然而，幾何貼體網(wǎng)格的數(shù)量大、劃分復(fù)雜、網(wǎng)格系統(tǒng)存在較多交接面。楊瑞等[10]利用致動盤模型 (Actuator Disk Model，ADM) 開展了風(fēng)機尾流研究，該方法極大簡化了網(wǎng)格系統(tǒng)，但忽略了切向誘導(dǎo)力。任會來等[11,12]分別利用考慮切向誘導(dǎo)力的致動盤模型(ADM-R)和更精細的致動線模型(Actuator Line Model，ALM)，研究了不同葉尖速度比條件下的流場特點，但未考慮湍流來流影響。

本文基于三種致動方法 (ADM，ADM-R，ALM)進行風(fēng)機尾流研究，分析來流為湍流時不同葉尖風(fēng)速比(λ=5.52，9.69)情況下各致動模型的仿真性能，采用最為準(zhǔn)確的致動模型獲得全場信息，研究流場動量平衡中各平衡要素的分布特點，為后續(xù)風(fēng)機尾流解析模型的建立奠定基礎(chǔ)。

1 控制方程

大渦模擬(Large Eddy Simulation，LES)利用亞網(wǎng)格尺度模型模擬小尺度渦旋對大尺度渦旋的影響，以獲得濾波后的大尺度渦旋運動，濾波后的連續(xù)方程和動量方程如下：

(1)

(2)

(3)

(4)

(5)

Ls=min(κd,CsV1/3)

(6)

式中：Ls為亞格子混合尺度；κ為von Kármán常數(shù)(κ=0.42)；d為距地面最近距離；V為控制體體積；Cs為Smagorinsky常數(shù)，采用0.1[15]。

2 致動方法

2.1 標(biāo)準(zhǔn)致動盤方法(ADM)

(7)

式中：AD為葉輪掃過面積；Δx為網(wǎng)格順風(fēng)向長度(見圖1a，圖中Ω為葉輪轉(zhuǎn)速，r為葉根距，α為槳矩角)；U∞為遠場風(fēng)速。

圖1 葉素動量理論示意

2.2 帶旋轉(zhuǎn)的致動盤方法(ADM-R)

ADM-R方法基于葉素動量(Blade Element Momentum，BEM)理論，將葉片沿徑向分割為若干葉素，如圖1所示，每個葉素受到的升力(dL)與阻力(dD)由式(8)計算：

式中：Cl，Cd分別為升力系數(shù)與阻力系數(shù)，由升阻力系數(shù)曲線插值得到；c和r分別為截面弦長與截面半徑(見圖1b，另外圖中L，D分別為升力和阻力，γ，φ分別為截面扭角和入流角，F(xiàn)為合力，F(xiàn)x和Fθ分別為合力在軸向和切向的分力)；Vrel為葉片相對空氣流速，由圖1b所示的局部速度矢量圖求解：

(9)

式中：a為軸向誘導(dǎo)因子；a′為切向誘導(dǎo)因子。根據(jù)BEM理論，作用于葉素的推力(dT)和力矩(dQ)可分別由式(10)與式(11)計算：

dT=dLcosφ+dDsinφ

(10)

dQ=(dLcosφ+dDsinφ)r

(11)

式中：B為葉片數(shù)量。

而根據(jù)動量守恒定理，dT與dQ可分別由式(12)，(13)計算：

dT=2πrdrρU(1-a)2aU

(12)

dQ=2πrdrρU(1-a)2a′r2U

(13)

(14)

2.3 致動線方法(ALM)

ADM-R方法假定誘導(dǎo)力在風(fēng)輪環(huán)向均勻分布，ALM方法舍棄該假定，對于給定的葉尖速度比λ可求出Ω，進而確定任意時刻各葉素掃過的網(wǎng)格，然后在對應(yīng)網(wǎng)格上施加單位體積力源項。誘導(dǎo)力的求解思路與ADM-R方法一致，fi按式(15)計算：

(15)

式中：V為葉片掃過網(wǎng)格的體積。

3 風(fēng)洞試驗

本文以三菱重工MWt-1000的1∶100縮尺模型為研究對象[17]開展風(fēng)洞試驗，其風(fēng)輪半徑為 0.285 m，風(fēng)輪直徑為0.57 m，風(fēng)輪中心高度Hhub為0.7 m。試驗中設(shè)置尖劈與柵欄，以模擬不同湍流度的大氣邊界層，尖劈位于風(fēng)機上游6 m處，出風(fēng)口位于風(fēng)機下游5 m處，如圖2中所示。來流風(fēng)速設(shè)定10 m/s，調(diào)節(jié)模型風(fēng)機旋轉(zhuǎn)角速度以改變?nèi)~尖風(fēng)速比，試驗考慮了兩種不同葉尖風(fēng)速比工況，即風(fēng)機額定功率時低葉尖速度比(λ=5.52)和風(fēng)機最大發(fā)電功率時高葉尖速度比(λ=9.69)。在風(fēng)機位置(x=y=0)處，共布置17個測點(h=0.3，0.35，0.4，… m)，采樣時長為40 s。

圖2 實驗實物

4 數(shù)值模型

為驗證計算結(jié)果的準(zhǔn)確性，本文數(shù)值模型與實驗研究保持一致。計算域入口為速度入口條件 (u=10 m/s，?p/?n=0)；出口為壓力出口條件(壓強p=0，?u/?n=?v/?n=?w/?n=0)；頂面為對稱邊界條件 (?u/?n=?v/?n=0，w=0)；側(cè)面為對稱邊界條件 (?u/?n=?w/?n=0，v=0)，其余均為壁面邊界條件 (u=v=w=0 m/s，?p/?n=0)。采用三種致動模型分析兩種情況：低葉尖速度比(λ=5.52)、高葉尖速度比(λ=9.69)。來流風(fēng)速均為U0=10 m/s。工況如表1所示，其中工況0為不考慮風(fēng)機情況下湍流邊界層計算，用以獲得網(wǎng)格無關(guān)條件下的網(wǎng)格尺寸，網(wǎng)格無關(guān)性驗證在5.1節(jié)論述。最終采用的網(wǎng)格水平尺寸為20 mm，底面貼地網(wǎng)格的豎向尺寸為1 mm，z方向增長率為1.1，最大豎向尺寸為20 mm，網(wǎng)格總數(shù)約為250萬，圖3c為尖劈附近網(wǎng)格分布?？臻g離散采用有限體積法，時間離散方法為二階中心差分法，壓力速度解耦采用SIMPLE算法。計算時間步長設(shè)置為0.0001 s，每個時間步最大迭代次數(shù)為20 次，從第2 s開始統(tǒng)計，統(tǒng)計時長為7 s。

表1 工況匯總

圖3 計算模型示意

5 計算結(jié)果

5.1 網(wǎng)格無關(guān)解

5.2 平均風(fēng)速

圖5 不同葉尖速度比情況下三種致動模型y = 0平面歸一化平均風(fēng)速分布

圖6 不同葉尖速度比情況下三種致動模型y-z平面歸一化平均風(fēng)速分布

5.3 脈動風(fēng)速

在y=0平面，繪制三種致動模型計算得到的u′/U0分布圖(見圖7)。并在u′/U0分布圖上繪制位于x=2D，4D，6D，8D處u′/U0剖面。可以發(fā)現(xiàn)，經(jīng)過風(fēng)輪平面后，脈動風(fēng)速增加，在葉尖位置處脈動值最大，葉尖速度比越高風(fēng)機對流場擾動效果越明顯，計算結(jié)果與實驗結(jié)果趨勢一致。對比三種致動模型可見，無論葉尖速度比高低、尾流距離遠近，ADM方法計算結(jié)果與實驗值相差均較大。ADM方法最大相對誤差出現(xiàn)在高葉尖速度比x= 2D位置，相對誤差達到70%，在遠尾流區(qū)相對誤差有所減小，但也到達40%。因此，為了較為準(zhǔn)確地預(yù)測風(fēng)機尾流脈動風(fēng)場分布，不可采用ADM方法。ADM-R方法與ALM方法在遠尾流域表現(xiàn)基本一致，但在近尾流區(qū)對于高葉尖速度比情況，ALM方法過高估計脈動風(fēng)速，最大相對誤差達到19%。ADM-R方法預(yù)測的風(fēng)機尾流脈動風(fēng)場與實驗值最為吻合。根據(jù)三種致動模型平均風(fēng)場與脈動風(fēng)場預(yù)測的綜合表現(xiàn)，ALM方法結(jié)果最為可靠。因此，采用ALM方法得到的尾流風(fēng)場做進一步動量平衡分析。

圖7 不同葉尖速度比情況下三種致動模型y = 0平面歸一化脈動風(fēng)速分布

5.4 動量平衡

為了研究流體微元各動量成分平衡關(guān)系，將N-S方程做時均化處理，獲得雷諾平均后順風(fēng)向動量方程：

(16)

在風(fēng)機中心處，各動量平衡項沿順風(fēng)向的分布如圖8所示。對于低葉尖速度比λ=5.52，在0

圖8 不同葉尖速度比情況下動量平衡方程各項分布

圖9 不同葉尖速度比情況下對流項各子項分布

湍流應(yīng)力各子項包括x應(yīng)力項(?u′u′/?x)、y應(yīng)力項(?u′v′/?y)和z應(yīng)力項(?u′w′/?z)，其順風(fēng)向分布見圖10。對于低葉尖速度比，λ=5.52，x應(yīng)力項基本為0，y應(yīng)力項與z應(yīng)力項為湍流應(yīng)力主要成分，見圖10a。高葉尖速度比時，λ=9.69，在近尾流區(qū)x<1.5D，x湍流應(yīng)力出現(xiàn)明顯峰值，且大于y應(yīng)力項與z應(yīng)力項，不可忽略。當(dāng)x≥1.5D時，x應(yīng)力項貢獻基本消失，y應(yīng)力項與z應(yīng)力項為湍流應(yīng)力主要成分，見圖10b。綜合以上分析，在遠尾流區(qū)(x>2.5D)動量平衡方程可簡化為：

圖10 不同葉尖速度比情況下湍流應(yīng)力項各子項分布

(17)

而在近尾流區(qū)(x≤2.5D)，對于低葉尖速度比，動量平衡方程可簡化為：

(18)

在遠尾流區(qū)(x>2.5D)，對于低葉尖速度比，動量平衡方程可簡化為：

(19)

采用致動模型結(jié)合大渦模擬的數(shù)值模擬方法雖可得到滿意的尾流風(fēng)場預(yù)測，但當(dāng)采用人工智能方法快速優(yōu)化風(fēng)機布置時，需要反復(fù)多次進行尾流風(fēng)場計算，數(shù)值模擬方法顯然難以勝任。而由于N-S方程的強非線性與多物理場的強耦合性，采用未簡化的N-S方程難以獲得適用于風(fēng)機尾流評估的快速解析模型。本文所獲得的簡化動量平衡方程略去了若干非線性項，為解析模型的建立提供了可能，而解析模型可快速高效地獲得允許誤差范圍內(nèi)的風(fēng)機尾流流場分布，可為智能化風(fēng)機優(yōu)化布置提供更為可靠的輸入。

5.5 計算資源耗費

利用LES進行求解需要耗費較多計算資源。在18 核Intel core i9-7980XE 服務(wù)器上并行計算完成無風(fēng)力機作用下大氣邊界層數(shù)值模擬耗時120 h左右。當(dāng)考慮風(fēng)機對流場的作用時，計算時長將會增加。圖11列出了在相同網(wǎng)格系統(tǒng)、相同求解設(shè)置、相同配置條件下100個時間步的計算耗時?？梢钥闯?，ADM方法計算耗時最短，比未放置風(fēng)力機條件下增加4.34%；ADM-R方法增加10.45%；ALM方法耗時增加25.70%。這是由于ADM方法未引入切向誘導(dǎo)力，對流場的擾動較?。籄DM-R方法雖然引入了切向誘導(dǎo)力，但誘導(dǎo)力不隨時間變化；而ALM方法不僅引入了切向誘導(dǎo)力，且真實再現(xiàn)葉片各葉素誘導(dǎo)力的時空分布，計算得到的流場最為紊亂，計算耗時最長。

圖11 高葉尖速度比情況三種致動模型計算耗時比較

6 結(jié) 論

本文利用三種致動方法，開展了湍流來流條件下不同葉尖速度比風(fēng)機尾流數(shù)值模擬，獲得了平均風(fēng)場與脈動風(fēng)場信息，分析了動量平衡關(guān)系，得出以下結(jié)論：

(1)對于平均風(fēng)速，ALM方法數(shù)值模擬結(jié)果與實驗值最為接近。在遠尾流區(qū)ADM與ADM-R方法模擬結(jié)果基本可靠，但在近尾流區(qū)，模擬值與實驗值有較大偏離。

(2)對于脈動風(fēng)速，ADM-R方法數(shù)值模擬結(jié)果與實驗值最為接近。ADM方法過低估計尾流脈動風(fēng)速，而在近尾流區(qū)ALM方法計算結(jié)果較實驗結(jié)果偏大。

(3)在近尾流區(qū)，動量平衡主要發(fā)生在壓強梯度項、對流項與湍流應(yīng)力項之間；在遠尾流區(qū)，壓強梯度項基本為0，動量平衡主要表現(xiàn)在對流項與湍流應(yīng)力項之間。

(4)在動量平衡中，對流項無論葉尖速度比高低，x對流項始終占據(jù)主導(dǎo)，其余各子項貢獻均可忽略。遠尾流區(qū)x湍流應(yīng)力項可以忽略。分區(qū)域建立了簡化動量平衡方程，為后續(xù)風(fēng)機尾流解析模型的建立奠定基礎(chǔ)。

(5)在計算耗時方面，ALM方法、ADM-R方法，與ADM方法耗時依次減少。