風(fēng)機(jī)尾流區(qū)湍流場(chǎng)特性數(shù)值仿真

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

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

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

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

1 控制方程

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

(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 致動(dòng)方法

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

(7)

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

圖1 葉素動(dòng)量理論示意

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

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

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

(9)

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

dT=dLcosφ+dDsinφ

(10)

dQ=(dLcosφ+dDsinφ)r

(11)

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

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

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

(12)

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

(13)

(14)

2.3 致動(dòng)線方法(ALM)

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

(15)

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

3 風(fēng)洞試驗(yàn)

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

圖2 實(shí)驗(yàn)實(shí)物

4 數(shù)值模型

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

表1 工況匯總

圖3 計(jì)算模型示意

5 計(jì)算結(jié)果

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

5.2 平均風(fēng)速

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

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

5.3 脈動(dòng)風(fēng)速

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

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

5.4 動(dòng)量平衡

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

(16)

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

圖8 不同葉尖速度比情況下動(dòng)量平衡方程各項(xiàng)分布

圖9 不同葉尖速度比情況下對(duì)流項(xiàng)各子項(xiàng)分布

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

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

(17)

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

(18)

在遠(yuǎn)尾流區(qū)(x>2.5D)，對(duì)于低葉尖速度比，動(dòng)量平衡方程可簡(jiǎn)化為：

(19)

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

5.5 計(jì)算資源耗費(fèi)

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

圖11 高葉尖速度比情況三種致動(dòng)模型計(jì)算耗時(shí)比較

6 結(jié) 論

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

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

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

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

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

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