風(fēng)力機(jī)遠(yuǎn)尾流的計(jì)算研究

田琳琳，趙 寧，鐘 偉，胡 偶

(南京航空航天大學(xué)航空宇航學(xué)院，江蘇南京 210016)

0 引言

風(fēng)力機(jī)是將自然界的風(fēng)能轉(zhuǎn)化為機(jī)械能并獲得電能的裝置。由于能量的轉(zhuǎn)移，風(fēng)經(jīng)過旋轉(zhuǎn)的風(fēng)力機(jī)之后，流動(dòng)情況發(fā)生了很大的變化:風(fēng)速減小、湍流強(qiáng)度增加、出現(xiàn)了明顯的風(fēng)剪切層。風(fēng)速減小會(huì)使下游風(fēng)力機(jī)的輸出功率降低，尾跡附加的風(fēng)剪切和強(qiáng)湍流會(huì)影響下游風(fēng)力機(jī)的疲勞載荷、使用壽命和結(jié)構(gòu)性能。經(jīng)過一段距離之后，在周圍氣流的作用下，風(fēng)速逐漸得到恢復(fù)。這就是風(fēng)力機(jī)的尾流效應(yīng)［1］。為了使有限面積的風(fēng)場產(chǎn)生最大的經(jīng)濟(jì)效益，風(fēng)電場開發(fā)商盡可能多在風(fēng)場內(nèi)布置風(fēng)力機(jī)，這樣上游機(jī)組會(huì)對下游機(jī)組產(chǎn)生較大的尾流效應(yīng)，導(dǎo)致下游機(jī)組發(fā)電量減少，進(jìn)而降低整個(gè)風(fēng)場的總輸出功率。因此，開展尾流效應(yīng)對風(fēng)力機(jī)性能影響的研究對合理布置風(fēng)力機(jī)，減小尾流對下游的影響，提高風(fēng)能資源的利用效率，減小風(fēng)電場的占地面積，最終使風(fēng)電場的經(jīng)濟(jì)性達(dá)到最佳起著重要的作用。

風(fēng)力機(jī)的尾流區(qū)分為近尾流區(qū)和遠(yuǎn)尾流區(qū)。近尾流區(qū)的研究著眼于功率提取的物理過程和風(fēng)機(jī)性能，而遠(yuǎn)尾流區(qū)的研究重點(diǎn)是尾流模型、地形影響、湍流模型的恰當(dāng)選取，更著重于研究風(fēng)電場機(jī)組間的相互干擾［1］。

許多專家學(xué)者對風(fēng)力機(jī)的尾流效應(yīng)開展了廣泛的計(jì)算研究?？梢詫⑦@些研究大致分為兩類:一類是尾流模型研究，這些模型是由學(xué)者提出的簡化的尾流模型，然后利用實(shí)驗(yàn)數(shù)據(jù)檢驗(yàn)?zāi)Ｐ筒⑦M(jìn)行修正，主要用來解決工程問題。如風(fēng)資源分析軟件WAsP采用的Park模型［2］以及Ainslie提出的渦粘性尾流模型［3］等等。另一類是流動(dòng)機(jī)理研究，該部分主要是回歸基礎(chǔ)研究，旨在揭示流動(dòng)的內(nèi)在規(guī)律和本質(zhì)。如早期的動(dòng)量－葉素理論［4］和渦尾跡方法［5］。之后隨著計(jì)算機(jī)的發(fā)展，又出現(xiàn)了很多基于CFD的方法，例如Ivanell等學(xué)者［6］使用了制動(dòng)線(ACL)方法研究了單個(gè)風(fēng)力機(jī)的近尾跡區(qū)的流動(dòng)狀況，Sankar等學(xué)者構(gòu)造了N－S方程、全勢能方程與自由渦尾跡相混合的方法用于求解水平軸風(fēng)力機(jī)周圍的非定常粘性流場［7］。S?rense［8］等使用 RANS 方程進(jìn)行流場的全三維數(shù)值模擬。

本文提出了兩種方法用來研究單個(gè)風(fēng)力機(jī)的遠(yuǎn)尾跡特征。第一種方法為修正的Park模型，由初始的一維Park模型進(jìn)行三角函數(shù)修正得到。Park模型由于計(jì)算量小且相對準(zhǔn)確［9］，在工程上得到了廣泛的應(yīng)用，經(jīng)過修正的Park模型計(jì)算量也很小，且在徑向方向的速度型更接近于真實(shí)的流場分布情況。第二種方法是結(jié)合制動(dòng)盤理論和CFD［10］，將風(fēng)力機(jī)簡化成為一個(gè)無限薄的產(chǎn)生壓力躍升的制動(dòng)盤，采用FLUENT商用軟件提供的Fan邊界對流場進(jìn)行數(shù)值模擬。該方法較之建立真實(shí)的風(fēng)力機(jī)模型并進(jìn)行模擬來說，降低了計(jì)算量和計(jì)算復(fù)雜度，為風(fēng)力機(jī)的微觀選址提供很好的理論支持。最后將本文的兩種數(shù)值方法的結(jié)果與風(fēng)洞實(shí)驗(yàn)結(jié)果以及Park模型的預(yù)估結(jié)果進(jìn)行對比分析，發(fā)現(xiàn)這兩種方法都能夠很好地預(yù)測風(fēng)力機(jī)的尾流流場。

1 研究模型

1989年 Garrad Hassan在 Marchwood工程實(shí)驗(yàn)室(MEL)開展了一系列的風(fēng)力機(jī)尾流的實(shí)驗(yàn)研究［11］，所選用的實(shí)驗(yàn)對象為模型縮比尺度為1/160的水平軸風(fēng)力機(jī)。對于真實(shí)尺度的風(fēng)力機(jī)來說，能反映其性能的重要的參數(shù)就是推力系數(shù)，所以在風(fēng)洞試驗(yàn)中，觀察了風(fēng)力機(jī)在三個(gè)不同尖速比下的氣動(dòng)性能。三個(gè)不同尖速比 λ =2.9，λ =4.0 ，λ =5.1 對應(yīng)的軸向推力系數(shù)分別為 CT=0.62，CT=0.79 ，CT=0.85 。

本文選用的計(jì)算模型為Hassan實(shí)驗(yàn)的縮比模型所對應(yīng)的真實(shí)模型，即風(fēng)輪直徑為43.2m，輪轂高度為50m，地表粗糙度為0.075m，相當(dāng)于風(fēng)力機(jī)樹立于短的農(nóng)作物覆蓋的真實(shí)地形之上。

2 計(jì)算方法

2.1 Park模型的修正

Park模型是由學(xué)者Jensen于1986年提出的一維尾流模型［2］，該模型假設(shè)在風(fēng)力機(jī)下游，尾流是線性膨脹的，Park模型的原理圖如圖1所示。根據(jù)質(zhì)量守恒定理推導(dǎo)出在風(fēng)力機(jī)下游x位置處速度為:

其中u0為來流風(fēng)速，a為軸流誘導(dǎo)因子，與推力系數(shù)CT有關(guān)，r1為下游x位置處尾流區(qū)域的半徑，α為尾流擴(kuò)散系數(shù)，表示尾流的膨脹速率。a、r1和α的公式如下:

其中，rr為風(fēng)輪半徑，z為風(fēng)力機(jī)輪轂高度，z0為地表粗糙度。

圖1 尾流模型原理圖Fig.1 Schematic of wake model

大量的實(shí)驗(yàn)結(jié)果以及數(shù)值模擬表明，在真實(shí)的尾流流場中，沿風(fēng)力機(jī)徑向的速度分布應(yīng)該為拋物型的［12］，這與三角函數(shù)的圖形有著很大的相似性，故在此使用三角函數(shù)對一維的Park模型進(jìn)行修正，使其為二維模型，保證修正后的徑向速度型為拋物型。修正后的Park模型示意圖如圖1所示。

令y方向?yàn)閺较蚍较颍僭O(shè)修正后的速度分布符合標(biāo)準(zhǔn)的三角函數(shù)關(guān)系式:

下面將分別求解A，k以及b。根據(jù)余弦函數(shù)周期的定義可得:

假設(shè)余弦函數(shù)取最大值時(shí)所對應(yīng)的徑向位置為r1即膨脹尾流的半徑，則有:

基于修正前后動(dòng)量守恒的原則，也就是說修正后的總動(dòng)量與Park模型推導(dǎo)出的總動(dòng)量相等，即速度型下所包含的面積相等，有:

聯(lián)立(6)～(8)式分別得到:

可得u*的表達(dá)式為:

所以經(jīng)過修正的Park模型分為兩步:

2.2 制動(dòng)盤理論結(jié)合CFD方法

當(dāng)空氣經(jīng)過風(fēng)輪盤時(shí)，由于風(fēng)力機(jī)的阻塞作用，使得流管內(nèi)的風(fēng)速降低，根據(jù)質(zhì)量守恒原理，流管的橫截面積就會(huì)膨脹以適應(yīng)減速的空氣，我們將風(fēng)輪所在平面稱為制動(dòng)盤［12］。

由經(jīng)典的Rankine－Froude原理可知，風(fēng)經(jīng)過制動(dòng)盤時(shí)，會(huì)在制動(dòng)盤前后產(chǎn)生壓差，壓差的表達(dá)式如下:

其中為△P空氣經(jīng)過制動(dòng)盤時(shí)產(chǎn)生的壓力降，ρ為空氣密度，U0為來流速度5.3m/s，a為軸流誘導(dǎo)因子，見式(2)。

把風(fēng)力機(jī)簡化為產(chǎn)生壓力降的制動(dòng)盤，這與商用軟件FLUENT提供的Fan的物理模型具有相似性，所以可借助Fan邊界來模擬風(fēng)輪，計(jì)算Fan對整個(gè)流場的影響。在FLUNT中，把Fan假設(shè)為一個(gè)無限薄的面，流體經(jīng)過這個(gè)面時(shí)出現(xiàn)壓力躍降，而躍降的大小是速度的多項(xiàng)式函數(shù)，壓強(qiáng)躍降函數(shù)的表達(dá)式為:

在式(13)中△P為壓強(qiáng)躍降量，fn為多項(xiàng)式系數(shù)，v為速度。這樣做的目的是把Fan前后的速度變化轉(zhuǎn)化成壓差項(xiàng)反映到流動(dòng)方程中。

在本文的計(jì)算中，將風(fēng)力機(jī)簡稱為一個(gè)無限薄的制動(dòng)盤，并根據(jù)風(fēng)力機(jī)的推力系數(shù)與風(fēng)速關(guān)系的特征曲線設(shè)定制動(dòng)盤前后的壓差。然后在全流場內(nèi)求解帶湍流模型的N－S方程。

3 計(jì)算結(jié)果與分析

3.1 風(fēng)輪中心軸向速度分布

能否準(zhǔn)確地計(jì)算尾流軸向動(dòng)量虧損是判斷一個(gè)模型是否能夠很好地預(yù)測尾流發(fā)展情況的一個(gè)重要指標(biāo)。圖2給出了在三個(gè)不同尖速比下，制動(dòng)盤結(jié)合CFD方法、Park模型以及風(fēng)洞實(shí)驗(yàn)這三種方法分別得到的沿風(fēng)輪中心的速度輪廓圖。在圖中，x方向是來流速度方向，D為風(fēng)輪直徑，U是尾流的當(dāng)?shù)厮俣?，U0是無窮遠(yuǎn)處的來流速度。由于Park模型在風(fēng)輪下游的兩倍風(fēng)輪直徑之后才有效，所以本文選擇從風(fēng)力機(jī)下游的2D距離之后對三種方法進(jìn)行分析比較。圖中很明顯地反映了風(fēng)力機(jī)的下游速度會(huì)下降，經(jīng)過一段距離之后，逐步恢復(fù)到來流速度。且隨著尖速比的增加，推力系數(shù)就會(huì)增加，尾流效應(yīng)就越明顯。在下游6D之前，Park模型和本文的制動(dòng)盤結(jié)合CFD方法都偏大很多，是因?yàn)檫@兩種方法都沒有考慮到風(fēng)力機(jī)的旋轉(zhuǎn)效應(yīng)，而真實(shí)的風(fēng)洞試驗(yàn)中，風(fēng)力機(jī)是旋轉(zhuǎn)的，旋轉(zhuǎn)會(huì)引起風(fēng)力機(jī)的下游產(chǎn)生很強(qiáng)的湍流以及尾渦，尾渦通過發(fā)展、破碎等過程會(huì)耗散掉空氣中很多的能量，導(dǎo)致下游的速度進(jìn)一步降低。Helmis等學(xué)者［13］指出風(fēng)洞實(shí)驗(yàn)過高預(yù)估了近尾跡區(qū)的尾流效應(yīng)，使分析得出的風(fēng)輪中心的速度虧損較大。這也是本文方法的結(jié)果與風(fēng)洞試驗(yàn)結(jié)果在近尾流區(qū)相差較大的原因之一。

開展尾流效應(yīng)研究的目的是在風(fēng)電場內(nèi)合理布置風(fēng)力機(jī)，風(fēng)力機(jī)的安裝間距要滿足風(fēng)場總體效益最大化的目標(biāo)。通過對國內(nèi)外風(fēng)場多年的建設(shè)經(jīng)驗(yàn)分析，風(fēng)力發(fā)電機(jī)組安裝間距在盛行風(fēng)向上選擇為6～10倍風(fēng)輪直徑，在垂直盛行風(fēng)向上選擇為3～5倍的風(fēng)輪直徑較為合適。在下游距離為6～15倍風(fēng)輪直徑之間，本文的制動(dòng)盤結(jié)合CFD方法的計(jì)算結(jié)果與實(shí)驗(yàn)值吻合的很好，而Park模型的計(jì)算結(jié)果偏高。這說明本文的計(jì)算方法可以為風(fēng)電場的微觀選址提供很好的理論依據(jù)。

圖2 風(fēng)輪中心軸向速度分布Fig.2 Velocity deficit on the wake centerline

在風(fēng)力機(jī)下游15倍風(fēng)輪直徑之后，Park模型的計(jì)算值比實(shí)驗(yàn)值偏大，本文的計(jì)算結(jié)果偏小，說明尾流還沒有來得及完全恢復(fù)。

3.2 下游位置的徑向速度分布

圖3 尖速比2.9時(shí)，風(fēng)力機(jī)下游的不同位置處的徑向速度分布Fig.3 The crosswind profiles of mean velocity at hub height for the tip speed ratio of 2.9

圖3和圖4分別給出了不同尖速比下風(fēng)力機(jī)下游的不同橫截面在風(fēng)輪高度處的徑向速度分布。該圖明顯的顯示了尾流效應(yīng)，且尖速比越大，尾流效應(yīng)越明顯。從整體來看，Park模型以及本文的制動(dòng)盤結(jié)合CFD方法的計(jì)算結(jié)果都大于實(shí)驗(yàn)結(jié)果，這是因?yàn)樵陲L(fēng)力機(jī)下游的一段距離之后，壓差不再在整個(gè)流動(dòng)中起主導(dǎo)作用，而是湍流強(qiáng)度以及尾渦的變化主導(dǎo)著整個(gè)流場的能量變化，在真實(shí)情況中，風(fēng)力機(jī)下游的尾渦會(huì)逐漸耗散掉，導(dǎo)致能量的降低，進(jìn)而速度會(huì)降低。在越靠近下游的地方，制動(dòng)盤結(jié)合CFD方法以及Park模型結(jié)果會(huì)更接近風(fēng)洞試驗(yàn)值，尤其是在尖速比為2.9時(shí)，風(fēng)力機(jī)下游10倍風(fēng)輪直徑的地方，本文結(jié)果與風(fēng)洞試驗(yàn)值幾乎吻合。從圖中還可以看出經(jīng)過修正的Park模型與風(fēng)洞試驗(yàn)結(jié)果和制動(dòng)盤理論結(jié)合CDF方法計(jì)算結(jié)果吻合的較好，尤其是在風(fēng)輪中心兩側(cè)，與實(shí)驗(yàn)值更加接近。由于其計(jì)算量小，計(jì)算復(fù)雜度低，因此在今后可作為一項(xiàng)計(jì)算工具應(yīng)用于工程實(shí)際。圖3和圖4還反映了制動(dòng)盤結(jié)合CFD方法計(jì)算處的尾流寬度略大于Park模型的尾流寬度，更接近于實(shí)驗(yàn)值。一個(gè)比較有趣的現(xiàn)象是在流管邊緣的地方，制動(dòng)盤結(jié)合CFD方法計(jì)算的結(jié)果比實(shí)驗(yàn)值和Park模型的值偏大2%，風(fēng)出現(xiàn)了稍微的加速現(xiàn)象，這是因?yàn)樵谥苿?dòng)盤的周圍拖出來的尾渦從周圍流體中吸收了能量，導(dǎo)致速度的稍微上升。

圖4 尖速比5.1時(shí)，風(fēng)力機(jī)下游的不同位置處的徑向速度分布Fig 4 The crosswind profiles of mean velocity at hub height for the tip speed ratio of 5.1

4 總結(jié)

本文根據(jù)GH公司提供的風(fēng)力機(jī)數(shù)據(jù)，通過兩種方法對單個(gè)風(fēng)力機(jī)的尾流進(jìn)行了計(jì)算模擬。一種方法是構(gòu)造新的尾流模型計(jì)算流場速度分布，另一種方法是結(jié)合制動(dòng)盤理論和CFD，在全流場進(jìn)行數(shù)值模擬。當(dāng)風(fēng)吹過風(fēng)力機(jī)時(shí)，風(fēng)速會(huì)下降，經(jīng)過下游一段距離慢慢恢復(fù)到來流速度，且尖速比越大，尾流效應(yīng)越明顯。從風(fēng)輪中心軸向速度分布圖來看，在下游6倍風(fēng)輪直徑之后，本文結(jié)果與風(fēng)洞試驗(yàn)值的最大誤差約為4%。徑向速度分布顯示了尾流的膨脹系數(shù)大概為0.08，越靠近下游，本文計(jì)算方法的精度就越高。雖然Park模型沒有精確預(yù)測尾渦的形狀，但是它近似給出了尾渦的發(fā)展過程，而本文構(gòu)造的新尾流計(jì)算模型與風(fēng)洞實(shí)驗(yàn)值吻合較好。

將本文兩種計(jì)算方法的結(jié)果與風(fēng)洞實(shí)驗(yàn)以及Park尾流模型進(jìn)行分析比較，顯示了本文能夠定性反映風(fēng)力機(jī)周圍的物理流動(dòng)現(xiàn)象。以上結(jié)果表明，本文的兩種計(jì)算方法適用于風(fēng)力機(jī)的流場模擬，可作為工程應(yīng)用的工具，用于風(fēng)電場的流場模擬，為風(fēng)力機(jī)的微觀選址提供一定的參考依據(jù)。

［1］VERMEER L J，S?RENSEN J N，CRESPO A.Wind turbine wake aerodynamics［J］.Progress in Aerospace(Sciences)，2003，39(6 －7):467 －510.

［2］JENSEN N O.A note on wind generator interaction［R］//Ris?National Laboratory，1983.

［3］AINSLIE J F.Development of an eddy viscosity model for wind turbine wakes［A］.Proceeding of 7 th BWEA Wind Energy Conference［C］，Oxford，1985.

［4］GLAUERT H.Airplane propellers［A］.In Aerodynamics theory.Dover Publication Inc:New York［M］，1963，4(Division L):169－360.

［5］BAREISS R，GUIDATI G，WAGNER S.Wake simulation for wind turbines with a free，prescribed and hybrid wake method［A］.Proc.10th IEA Symp.On the Aerodynamics of Wind Turbines［C］，Edinburgh，1996:23－28.

［6］IVANELL S，MIKKELSEN R，S?RENSEN J，HENNINGSON D.A stability analysis of the tip vortices of a wind turbine［J］.Wind Energy，2009.

［7］XU G，SANKER L N.Computational study of horizontal axis wind turbines［R］.AIAA Paper 99－0042，1999.

［8］GLAUERT H.Airplane propellers［A］.In Aerodynamics theory［M］.Dover Publication Inc:New York，1963，4(Division L):169－360.

［9］BARTHELMIE R，F(xiàn)OLKERTS，LARSEN G C，et al.Comparison of wake model simulation with offshore wind turbine wake profiles measured by sodar［J］.Journal of Atmospheric and Oceanic Technology，23(7):888-901.

［10］DZUIBINSKI Z，WIENZYSLAW S.Vortex ring state simulation using actuator disc［EB/OL］.http://www.scs－europe.net/conf/ecms2007/ecms2007－cd/ecms2007/ecms2007%20pdf/ese_0121.pdf.

［11］SCHLEZ W，TINDAL A，QUARTON D.GH windfarmer validation report［R］.Bristol:Garrad Hassan and Partners Ltd，2003.

［12］EBERLE W R.Wind flow characteristics in the wake of large wind turbines，Vol I，analytical development［R］.DOE/NASA CR－165380.

［13］賀德馨，等.風(fēng)工程與工業(yè)空氣動(dòng)力學(xué)［M］.北京:國防工業(yè)出版社，2006.

［14］HELMIS C G，PAPADOPOULOS K H，ASIMAKOPOULOS D N，PAPAGEORGAS P G，SOILEMES A T.An experimental study of the near wake structure of a wind turbine operating over complex terrain［J］.Solar Energy，1995，54(6):413-28.