李辰奇， 許 昌， 楊建川， 甄文忠， 劉德有， 鄭 源

（1.河海大學(xué)能源與電氣學(xué)院，南京 211100；2.丹麥科技大學(xué)風(fēng)能系，林格比 2800）

基于CFD的復(fù)雜地形風(fēng)能分布研究

李辰奇1， 許 昌1， 楊建川1， 甄文忠2， 劉德有1， 鄭 源1

（1.河海大學(xué)能源與電氣學(xué)院，南京 211100；2.丹麥科技大學(xué)風(fēng)能系，林格比 2800）

提出一種通過計算流體力學(xué)（CFD）計算風(fēng)電場風(fēng)速和風(fēng)能分布的方法，該方法通過Argis軟件對風(fēng)電場CAD等高線圖進(jìn)行離散處理，按照風(fēng)速玫瑰圖分布，利用Fluent軟件對復(fù)雜地形條件下風(fēng)場按照12分度進(jìn)行數(shù)值計算，并按照風(fēng)速風(fēng)向分布概率由Tecplot軟件進(jìn)行后處理，得到整個風(fēng)場的風(fēng)速分布特點和風(fēng)能分布.將計算結(jié)果與WAsP軟件作比較，結(jié)合地形特點進(jìn)行分析，認(rèn)為本文計算方法得到的結(jié)論較為可靠.

風(fēng)電場；復(fù)雜地形；風(fēng)能資源評估；CFD數(shù)值模擬

近年來，風(fēng)電能源成為發(fā)展速度最快的可再生能源［1-3］，風(fēng)能資源評估是建設(shè)風(fēng)電場的基礎(chǔ)性工作，發(fā)電量計算結(jié)果直接關(guān)系到風(fēng)電場的經(jīng)濟(jì)指標(biāo)和整體效益.目前風(fēng)能資源評估和發(fā)電量計算基本上采用丹麥Riso實驗室開發(fā)的WAsP軟件，該軟件采用線性模型對風(fēng)能資源和風(fēng)能產(chǎn)量進(jìn)行評估，對于平坦地形條件，預(yù)測精度較高，而對于復(fù)雜地形條件，預(yù)測精度有待于進(jìn)一步提高［4-8］.隨著風(fēng)能資源開發(fā)的不斷深人，風(fēng)電場規(guī)劃選址開始從風(fēng)能穩(wěn)定、施工條件好的平坦地形向高湍流、施工難度大的復(fù)雜地形發(fā)展［9-11］.這樣，以WAsP為計算引擎的風(fēng)能資源評估軟件不適合處理復(fù)雜地形的風(fēng)資源情況.

目前我國已有學(xué)者開展針對復(fù)雜地形的風(fēng)資源預(yù)測的研究工作，但其研究水平與國外有一定差距，主要體現(xiàn)在地形數(shù)據(jù)的處理、風(fēng)電場空氣動力場的數(shù)值計算方法和空氣動力場數(shù)值計算結(jié)果的處理上［12-16］.本文嘗試開展適用于復(fù)雜地形風(fēng)資源評估的小尺度非線性數(shù)值模式研究，采用近年來快速發(fā)展的計算流體力學(xué)（CFD）技術(shù)，計算分析風(fēng)電場風(fēng)速和風(fēng)能分布，并與WAsP軟件的計算結(jié)果進(jìn)行比較，并定性定量分析本文方法的可靠性.

1 地形數(shù)字化

目前風(fēng)電場復(fù)雜地形通常是通過AutoCAD等高線方式表示，而通過CFD技術(shù)計算風(fēng)電場空氣動力場，需要建立包括地形在內(nèi)的物理模型，物理模型的建立首先需要對復(fù)雜地形進(jìn)行數(shù)字化.Argis是一款集地理信息系統(tǒng)開發(fā)、分析、地圖數(shù)字化、地理信息的采集等功能于一體的軟件，目前在地理信息系統(tǒng)行業(yè)該軟件的功能最為強大，應(yīng)用最為廣泛.本文實例中，使用Argis將AutoCAD中的等高線分離出來生成tin文件，從而得到轉(zhuǎn)化成坐標(biāo)點的dem文件，對dem文件進(jìn)行離散處理，風(fēng)電場等高線如圖1所示，而在Argis離散處理后地形如圖2所示.通過圖形的對比分析，認(rèn)為Argis軟件離散后的地形圖基本達(dá)到了原等高線圖的精度要求.

圖1 風(fēng)電場等高線圖Fig.1 Wind farm contour map

圖2 Argis地形顯示圖Fig.2 Argis mapping

根據(jù)復(fù)雜地形離散后的數(shù)字文件，在gambit中建立地形物理模型.根據(jù)地形資料，選取的計算區(qū)域為7 000 m×7 000 m×500 m，如圖3所示.其中，高度為500 m，選擇非結(jié)構(gòu)網(wǎng)格，地面面網(wǎng)格為30 m× 30 m，垂直方向上網(wǎng)格將地表到500 m高度分為3層，分別為0～50 m，50～200 m，200～500 m，而網(wǎng)格間隔分別為5，10，30 m.

圖3 風(fēng)場地形計算區(qū)域Fig.3 Wind farm calculation area

2 數(shù)學(xué)模型和求解

流體動力學(xué)標(biāo)準(zhǔn)的Navier-Stokes方程組（直角坐標(biāo)系）為［17］

本文選用穩(wěn)態(tài)、常物性控制方程.選用的標(biāo)準(zhǔn)κ-ε湍流模型［18］為

式中，u為絕對速度；ρ為空氣密度；u，v為水平風(fēng)速；w為垂直風(fēng)速；p為空氣壓強；μ為動力粘度；κ為湍流動能；ε為湍動耗散率；Γκ為κ的擴(kuò)散系數(shù)；G為紊流能量的生成率，也是動能的耗散率；G-ρε為方程的凈源項；C1ε和C2ε為模型常數(shù).

方程求解采用二階上風(fēng)方法離散，底面邊界條件是壁面函數(shù)，入口為給定的速度入口，出口為給定的出流口，而其它的面采用對稱條件.風(fēng)速按照12分度的風(fēng)向玫瑰圖，計算出每個分區(qū)的平均風(fēng)速，在垂直地面方向采用的模型為

式中，u*為表面摩擦速度系數(shù)；K為Karman常數(shù)，一般取0.4；z0為地表面粗糙度長度.

該模型通過Fluent的UDF編輯實現(xiàn).其中，0～30°風(fēng)向的入口速度輸入如圖4所示，可以看出，速度進(jìn)口沿著高度方向的變化趨勢，地面風(fēng)速接近為零，隨著高度的增加，根據(jù)網(wǎng)格劃分的情況，出現(xiàn)層狀的顏色分布.

圖4 0～30°方向速度進(jìn)口局部圖Fig.4 0～30°speed import local map

3 結(jié)果分析

3.1 風(fēng)能計算結(jié)果

風(fēng)能密度是氣流在單位時間內(nèi)垂直通過單位截面積的風(fēng)能［19］.從而風(fēng)功率密度公式，也稱風(fēng)能密度公式為

式中，W為風(fēng)功率；V為給定位置的風(fēng)速大小.

通過計算，在Fluent軟件中得到12個方向上的風(fēng)速分布，將這些風(fēng)速在Tecplot中按照風(fēng)速概率分布擬合，計算出風(fēng)電場中的風(fēng)能分布.

圖5 風(fēng)能分布Fig.5 Distribution of wind energy

在圖5（a）中，風(fēng)能隨著離地高度的增加也逐漸增加，風(fēng)能已經(jīng)在300～600 W／m2之間，也可以看出風(fēng)速主風(fēng)向的趨勢，大概是在45°左右.在圖5（b）和圖5（c）中，風(fēng)能分布趨勢已大致穩(wěn)定，反映出海拔高度變化平緩的區(qū)域風(fēng)速梯度變化小，海拔高度變化差異大的區(qū)域風(fēng)速梯度變化大.圖6和圖7是在地形圖的位置y=-3 000 m，x=2 500～5 000 m范圍內(nèi)的某一山坡下的風(fēng)能分布圖和速度分布圖.其中，圖6是將12個風(fēng)向上的風(fēng)速綜合而成，圖7是在0～30°風(fēng)向上從右往左吹過山坡.從圖7中可以看出，迎風(fēng)坡的風(fēng)速比背風(fēng)坡的風(fēng)速大，這是因為受到山坡的阻礙，風(fēng)受到削減，能量損失，造成風(fēng)速減小.

圖6 z＝1 700 m風(fēng)能分布Fig.6 Distribution of wind energy at z＝1 700 m

3.2 數(shù)值模擬與WAsP軟件進(jìn)行比較

圖7 z＝1 800 m風(fēng)速分布Fig.7 Distribution of wind energy at z＝1 600 m

表1 CFD與WAsP在20個計算點的風(fēng)能分布Tab.1 Distribution of wind energy at 20 calculation poins by using CFD and WAsP

在地圖上，沿著主風(fēng)向，取地圖對角線上均勻分布的20個點為參考點，通過CFD數(shù)值模擬方法和WAsP軟件方法，分別計算這些點的風(fēng)能大小.可以看到，CFD計算得到的點的風(fēng)能大小變化很大，而WAsP得到的風(fēng)能大小變化范圍較小，這主要是因為WAsP只考慮到高度對風(fēng)速的影響，而忽略了復(fù)雜地形下地形對風(fēng)的削減改變也有很強的作用，而CFD數(shù)值模擬方法將地形對風(fēng)速風(fēng)能的影響也考慮在內(nèi)，從而得到了更符合實際的風(fēng)能分布結(jié)果.從表1中可以看出，第6點開始到第20個點，CFD比WAsP的結(jié)果要小100～200 W／m2，而前面5個點可能靠近地形的邊緣，結(jié)果差的變化范圍較大.從圖8可以直觀地看到，只有第2個計算點的CFD計算風(fēng)能分布大于WAsP的計算結(jié)果，其它計算點都是WAsP計算結(jié)果比CFD結(jié)果大，這也與目前復(fù)雜地形利用WAsP計算結(jié)果偏大的結(jié)果吻合.表1是CFD與WAsP在20個計算點的風(fēng)能分布.

圖8 WAsP與CFD計算20個點的風(fēng)能圖Fig.8 Wind energy diagram at 20 calculation points by using CFD and WAsP

4 結(jié) 論

a.在復(fù)雜地形下，將CFD計算方法與WAsP軟件計算得到的風(fēng)能分布相比較，CFD計算方法更客觀、更接近實際地反映出風(fēng)經(jīng)過復(fù)雜地形下的風(fēng)能分布，同時考慮到了山坡等障礙物對風(fēng)的削弱，而不僅僅局限于高度要素的影響.

b.由于傳統(tǒng)風(fēng)資源評估軟件WAsP對復(fù)雜地形計算誤差較大，通常偏大，本文采用CFD數(shù)值模擬計算方法，考慮了復(fù)雜地形和高度對風(fēng)速變化的影響，計算結(jié)果比WAsP結(jié)果小，這與目前復(fù)雜地形利用WAsP計算結(jié)果偏大的結(jié)果吻合.

［1］ 韓春福.風(fēng)能資源評估方法的分析及應(yīng)用［J］.節(jié)能，2009，322（5）：22-24.

［2］ 胡毅，張健.風(fēng)資源評估中風(fēng)速方法研究［J］.內(nèi)蒙古科技與經(jīng)濟(jì)，2010，233（21）：76-78.

［3］ 梁思超，張曉東，康順，等.基于數(shù)值模擬的復(fù)雜地形風(fēng)場風(fēng)資源評估方法［J］.空氣動力學(xué)學(xué)報，2012，30（3）：415-421.

［4］ Okulov V I，Sorensen J N.Instability of a vortex wake behind wind turbines［J］.Doklady Physics，2004，49（12）：772-777.

［5］ Johansen J，Sorensen N N.Aerofoil characteristics from 3D CFD rotor computations［J］.Wind Energy，2004，7（4）：283-294，.

［6］ 賴永倫，巫卿.WAsP軟件在貴州四格風(fēng)電場風(fēng)資源評估中的應(yīng)用分析［J］.紅水河，2009，28（4）：106 -108.

［7］ 楊振斌，薛桁，桑建國.復(fù)雜地形風(fēng)能資源評估研究初探［J］.太陽能學(xué)報，2004，25（6）：744-749.

［8］ Yang X Y，Xiao Y，Chen S Y.Wind speed and generated power forecasting in wind farm［J］. Proceedings of the Chinese Socity for Electrical Engineering，2005，25（11）：1-2.

［9］ Berge E，Nyhammer F，Tallhaug L，et al.An evaluation of the WAsP model at a coastal mountainous site in Norway［J］.Wind Energy，2006，9（1）：131-140.

［10］ 王遠(yuǎn)成，吳文權(quán).不同形狀建筑物周圍風(fēng)環(huán)境的研究［J］.上海理工大學(xué)學(xué)報，2004，26（1）：19-23.

［11］ Kusiak A，Zheng H Y.Optimization of wind turbine energy and power factor with an evolutionary computation algorithm［J］.Renewable Energy，2010，35（3）：1324-1332

［12］ 馮長青，杜燕軍，包紫光，等.風(fēng)能資源評估軟件WAsP和WT的適用性［J］.中國電力，2010，43（1）：61 -65.

［13］ Madsen H A，Larsen G C，Larsen T J.Validation of the dynamic wake meandering model for implementation in an aeroelastic code［J］.Journal of Solar Energy Engineering，2010，132（4）：41-71.

［14］ Berg J，Mann J，Bechmann A.The bolund experiment，Part I：Flow over a steep，three-dimensional hill［J］. Boundary-Layer Meteorology，2011，141（2）：219.

［15］ Bingol F，Mann J，F(xiàn)oussekis D.Conically scanning lidar error in complex terrain［J］.Meteorologische Zeitschrift，2009，18（2）：189-196.

［16］ Kindler D，Oldroyd A，MacAskill A，et al.An eight month test campaign of the Qinetiq ZephIR system：preliminary results［J］.Meteorologische Zeitschrift，2007，16（5）：479-480

［17］ 王美琳，羅勇，周榮衛(wèi).WindSim軟件在復(fù)雜地形風(fēng)電場風(fēng)能資源評估中的應(yīng)用［J］.氣象，2010，36（2）：113 -119.

［18］ 王福軍.計算流體動力學(xué)分析［M］.北京：清華大學(xué)出版社，2004.

［19］ 田子嬋，楊勇平，劉永前.復(fù)雜地形的風(fēng)資源評估研究［D］.北京：華北電力大學(xué)，2009.

（編輯：石 瑛）

Prediction of Wind Energy Distribution in Complex Ter rain Using CFD

LIChen-qi1， XUChang1， YANGJian-chuan1， ZHENWen-zhong2， LIUDe-you1， ZHENGYuan1
（1.College of Energy and Electrical Engineering，Hohai University，Nanjing 211100，China；2.Department of Wind Energy，Technical University of Denmark，Lyngby 2800，Denmark）

Numerical simulations were carried out by using the Fluent software with a mesh generated by using the Argis software to predict the wind speed distribution in complex terrain. The post-processing by Tecplot was used to get the wind speed distribution characteristics and the wind energy distribution in the whole wind field.The obtained results were compared to the calculated results with WAsP software.

wind farm；complex terrain；wind energy resources evaluation；computational fluid dynamics numerical simulation

TM 344.1

A

1007-6735（2013）03-0270-05

2012-10-21

教育部留學(xué)回國人員科研啟動基金資助項目（2012-940）；國家重點基礎(chǔ)研究發(fā)展計劃（973計劃）資助項目（2010CB227102）；中央高?；究蒲袠I(yè)務(wù)費項目（2009B19414）

李辰奇（1988-），男，碩士研究生.研究方向：流體機械及工程.E-mail：lichenqi19881123＠163.com

許 昌（1973-），男，副教授.研究方向：從事風(fēng)力發(fā)電和太陽能發(fā)電研究.E-mail：zhuifengxu＠163.com