基于數值模擬的復雜地形風場風資源評估方法

梁思超，張曉東，康 順，2，康雅蘭，趙永鋒

（1.華北電力大學電站設備狀態監測與控制教育部重點實驗室，北京 102206；2.西安現代控制技術研究所，陜西 西安 710065；3.中國電力工程顧問集團華北電力設計院工程有限公司，北京 100120；4.中國福霖風能工程有限公司，北京 100034）

0 引 言

風電場風資源評估和微觀選址的好壞將直接關系到風電場發電量的多少。經過大量工程實踐，基于線性模型的WAsP軟件被認為不適用于復雜地形風場。計算流體力學（CFD）作為一種很有前途的手段，在這一領域正不斷發展。

湍流數值模擬方法主要有直接模擬（DNS）、大渦模擬（LES）、分離渦模擬（DES）以及雷諾時均模擬（RANS）。從精度上看，在實際山地Bolund山測風實驗的盲評中，排名靠前的均為兩方程的雷諾時均數值模擬方法，優于大渦模擬和風洞實驗結果。綜合計算精度，運算時間，計算機資源，雷諾時均法還沒有受到其他任何模擬方法的挑戰［1］。

帶有壁面函數的k－ε湍流模型，由于同時考慮了風速，湍動能以及地表粗糙度，適于風工程的研究。為把這一用于工業流動問題的模型用于風工程研究，前人對湍流模型參數進行了修正［2－3］。2009 年 在NUMECA軟件平臺下，提出其他相關參數對模擬的影響，作出進一步修正［4］。使用該修正模型的Askervein山數值模擬結果很好地與實驗值吻合。

為了將這一模擬方法用于實際風場的風資源評估和微觀選址中去，以實際南澳島風場為研究對象。結合測繪數據和衛星數據等間隔生成不同風向下的計算域，進行全周數值模擬。

1 南澳島簡介

為了簡化問題，降低邊界條件設置難度，選取島嶼地形南澳島風場作為研究對象。南澳島瀕臨臺灣海峽，處于南海北部，距離汕頭市萊蕪島約10km，面積106.45km2，島上山地面積占93.6%。島上平均風速達8.44m／s，年有效發電小時數達7215h，年平均有效風能密度為678W／m2，有東北和西南兩個主風向，風力資源居世界最佳之列。

圖1 南澳島Fig.1 Nanao island

圖2 固定測風塔位置Fig.2 Fixed masts positions

圖3 6號塔風向玫瑰圖Fig.3 Wind rose of 6＃

南澳風電場資源勘測期間，設立了若干固定測風塔和移動測風塔。其中，1、2、3、4、6、7號測風塔為固定測風塔，測有全年風資源數據，位置如圖2所示。從6號測風塔70m高度處風向玫瑰圖（圖3）中可以看到，南澳島的主風向為40°左右。

2 數值方法

采用NUMECA的FINE／TURBO軟件包，該軟件采用時間相關法求解湍流Navier－Stokes方程、中心節點的有限體積離散、顯式Runge－Kutta法，全多重網格初場處理以及多重網格迭代加速收斂。選用的湍流模型為帶有壁面函數的k－ε湍流模型。

2.1 k－ε兩方程湍流模型

標準k－ε湍流模型三維標準輸送方程為：

k和ε為湍動能及耗散率，湍動能生成項Pk為：

式（3）中，uj和分別為j方向的平均風速和脈動速度。湍流粘度為：

湍流模型參數參數Cμ、Cε1、Cε2、σk、σε在模擬大氣邊界層時設置為Cμ＝0.03，Cε1＝1.21，Cε2＝1.92，σk＝1.0，σε＝1.3［5］。為排除軟件對流場中μtx／μ的影響，設置 MUCLIP＝1×109。

數值迭代時，由于湍動能是從入口不斷向出口發展的，迭代速度很慢，不利于計算的收斂。將參數KEGRID的值設置為大于全多重網格加速重數的最大值。此時，湍動能k會在最粗網格上進行迭代，為較細網格的迭代提供好的初場，加速收斂，防止發散。計算中KEGRID取4。

2.2 壁面函數

在FINE／TURBO軟件中，不考慮粘性底層的存在，粗糙度z0為粗糙地表上風速為零點到粗糙地表頂端的距離。壁面函數的表達式如下：其中，u＊為摩擦速度，與壁面剪切應力τwall有關，u＊卡門常數κ一般取0.41。z0為表面粗糙度，d為零偏移量，定義u（z0＋d）＝0。當d≠0時，相當于把壁面移到距真實壁面距離d 的地方。對于粗糙壁時，常數B取零。其中，ksy＋／y，y＋＝yu＊／v。計算中取d＝0，B＝0。

3 計算域與邊界條件設置

實際風工程中，由于經費有限，只測繪風場內有限的地形區域。而在CFD中，為了考慮周圍地形對風場空氣流動的影響，需要更大區域的地形數據來生成計算域和網格。針對這一問題，提出結合網上可免費獲取的SRTM衛星地形數據，補全測繪地形數據的方案。衛星數據絕對水平和高程精度分別為20m和16m，僅作為地形補充，精度足夠。南澳島風場二期的測繪地形數據與衛星數據如圖4、圖5所示。地形處理時考慮了WGS84與北京54坐標的轉換、二維插值、兩地形間的平滑［5］。

圖4 測繪地形圖Fig.4 Mapping terrain data

圖5 SRTM衛星地形圖Fig.5 SRTM terrain data

以30°為間隔，分別進行12個風向的數值模擬。圖6為0°、30°、60°風向下的計算域和網格。由于風向正交，如0°、90°、180°、270°，這四個風向采用相同的計算域和網格。于是，模擬中僅需要三個計算域，根據風向設置相應的邊界條件（即風速的大小和方向）即可。另做40°主風向計算域，這樣就有4個不同的計算域。

圖6 風向角為0°、30°、60°的網格和計算域Fig.6 Grids and domains of 0°，30°and 60°degree

計算域尺寸為18000m×18000m×4000m。水平方向網格單元尺寸為90m×90m。計算域x、y、z三個方向的網格數為200、200、76。垂直方向網格分布分為三層結構：底層、中層和頂層（圖7）。底層網格共高10m，共有3個網格單元，按等延展比分布，第一層網格高度2.3m。中層內每個網格高度相等，均為5m，共12個。頂層網格也按等延展比分布，隨著高度的增加網格單元尺度增大。由于中層和底層內每個網格單元具有固定的高度，在后處理時，可以在一定高度上生成一個與地面起伏趨勢相同的網格面，便于分析模擬結果。

以40°風向計算為例（圖8），設置邊界條件時，來流方向邊界為入口邊界，對面一側為出口邊界，兩側為鏡像邊界，計算域頂端為入口邊界，地面為粗糙固壁。將計算域劃分為若干塊，便于在不同固壁區域設置不同的地表粗糙度。設置海面粗糙度為0.002m，陸地粗糙度為0.3m。

圖7 垂直方向網格分布Fig.7 Vertical grid distribution

圖8 邊界條件設置 Fig.8 Boundary conditions

入口邊界上速度、湍動能、以及耗散率由以下模型確定［6］。

其中，κ取0.4，Cu＝0.03。hg為大氣邊界層高度，hg＝u＊／6f，f為科里奧利參數，f＝2Ωsin a（｜λ｜），其中，Ω為地球自轉的角速度，λ為風場所處緯度位置。認為大氣邊界層頂端為湍動能消失的地方，當高度超過大氣邊界層高度時，速度不再變化，湍動能和耗散率為零。計算中，Ω＝7.292×10－5rad，λ取23.5°，z0取0.002m。以50m為特征高度，以風速2、4、6、8、10、12、14、16、18、20m／s作為特征值，根據式（6）分別計算u＊，生成來流條件。

4 模擬結果分析

4.1 不同來流條件對計算結果的影響

為了研究不同來流風速對于模擬結果的影響，引入無量綱量風加速ΔS，定義式為：

Uloc（h）為距地表高h處的當地風速，Uref（h）為參考風速。取50m高處的風速進行比較。各向同性時，湍流強度IT與湍動能k和風速U的關系為：

以0°和90°風向計算為例，分別選取6個位置（如圖9）。其中，A點和B點位于背風坡；C點和D點位于山頂；E點和F點位于迎風坡。圖10和圖11分別為各個位置50m高度上的風加速和湍流強度與入口特征風速的關系。雖然各點位置不同，流動狀況也不同，但所處位置上的風加速和湍流強度基本上不隨來流發生變化，其余風向模擬結果均存在這一規律。這對于后處理帶來了極大的方便，每個方向僅作一個風況的特征計算即可。

圖9 0°和90°風向下特征位置Fig.9 Characteristic positions of 0°and 90°

4.2 計算結果比較

由于以上規律的存在，可以認為：同風向下、不同風速時，風場不同位置兩處風速比例關系不變。因此，在一個風向下，當知道某一位置的風速U1時，可采用公式（11）計算另一位置風速U2。其中，Us1和Us2為該風向下，對應位置模擬值。

為驗證數值模擬的精度，進行各測風塔間的風速互相推算比較。除6號塔外，其余各塔1、2、3、4、7均有50m高度處的數據。以7號塔風向為風向基準，在風向和風速穩定的時間段內選取3個小時平均數據，見表1。由式（11），以各測風塔風速計算其它測風塔位置處風速，并比較差異，得到表2～表4。

表1 樣本風速Table 1 Sample wind speed data

表2 30°風向計算結果比較Table 2 30°results comparisons

表3 60°風向計算結果比較Table 3 60°results comparisons

表4 90°風向計算結果比較Table 4 90°results comparisons

從三個風向的互推差異來看（表2、表3、表4），數值模擬結果存在較大誤差。產生誤差的原因可能來自兩方面。一是數值模擬網格的太稀，不能很好地表現地面的起伏狀況，簡單的壁面函數不足以描述復雜的粗糙地表。二是由于安裝，損壞，附近障礙物影響等原因造成的測量數據不準。但無論怎樣，這一方法在風資源評估當中是值得借鑒的。

4.3 主風向下的流場分析及微觀選址

由于40°為該風場的主風向，研究該風向下的流場分布狀況是十分必要的。圖12中標出的7個位置，綜合空氣流動狀況，風能密度和湍流強度，適于安裝風力機。最下面的區域由于面積較大，地形平坦，便于施工，應優先考慮。

風場區域內有4個位置（圖13）處于背風處，流動狀況復雜，風能密度小，流動不平穩。尤其是最下那個區域，處于高山的背風坡，十分陡峭，又由于左側山地的阻擋作用，一股氣流斜向流入，加劇了該處的流動分離。這些位置的湍流強度明顯大于周圍，因此，要避免在這些位置安裝風力機。

5 全周模擬結果的后處理

實際工程當中，為了簡化計算，多采用威布爾分布來擬合測風數據，描述風速的概率分布。但很多情況下，風速概率分布不滿足威布爾分布。這時，可采用處理時間序列的方式來評估風資源。即，每一個測風數據均進行累計計算。

處理時，選擇一位于高處、受地形影響較小的測風塔作為基準測風塔，用來判定該時刻的風況屬于哪個風向的數值計算。在風場區域內布置二維網格面，利用式（11）計算節點i采用第j個測風塔的數據點得到的風速Uij。由于測風塔往往有n個，可采用反距離加權法計算該結點的風速Ui，即：

其中dij為節點i到測風塔j的距離。用每一時刻測風數據計算得到風能密度，累計相加后平均，即得到每個網格點上的年平均風能密度。得到風場區域年平均風能密度云圖，如圖14。

圖1 4 年平均風能密度云圖Fig.1 4 Annual average wind energy density contour

6 結 論

在相同計算域下，地表附近風加速和湍流強度模擬值基本上不隨入口參考風速的變化而發生變化。因此，當知道風場某一位置的風速和湍流強度時，可根據模擬值計算出另一位置的風速和湍流強度。基于這一思想，建立起的數值模擬和后處理方法對于實際風場的風資源評估具有積極的指導作用。

［1］HANJALIC K，KENJERES S.Some developments in turbulence modeling for wind and environmental engineering［J］.Journal of Wind Engineering and Industrial Aerodynamics，2008，96（10－11）：1537－1570.

［2］BECHMANN A，SORENSEN N N，JOHNSEN S，et al.Hybrid RANS／LES method for high Reynolds numbers，applied to atmospheric flow over complex terrain［A］.The Science of Making Torque from Wind［C］.Lyngby，DENMARK，2007.No.012054.

［3］YANG Y，GU M，CHEN S Q，et al.New inflow boundary conditions for modelling the neutral equilibrium atmospheric boundary layer in computational wind engineering［J］.J.of Wind Engineering and Industrial Aerodynamics，2009，97：88－95.

［4］梁思超，張曉東，康順.復雜地形風場擾流數值模擬方法［J］.工程熱物理學報，2011，32（6）：945－948.

［5］梁思超.復雜地形風場數值處理及擾流模擬方法研［D］.［碩士論文］.北京：華北電力大學，2010.

［6］ZHANG X D.CFD simulation of neutral ABL flows［R］.Ris?－R－1688（EN），2009.