基于Openfoam的懷來復雜地形風電場流動特性研究

葉昭良，閆 姝，史紹平，曾崇濟，陳新明，張 波，馬俊祥

（1.煤基清潔能源國家重點實驗室（中國華能集團清潔能源技術研究院有限公司），北京 102209； 2.華能河南清潔能源分公司，河南 鄭州 450000）

0 引言

復雜地形風電場的速度場流動特性的研究，對認識風電場場群的氣動變化規律、優化風力機控制和機組排布策略、提高風電場風能利用效率具有重要的意義。 在風工程計算模型中，大氣的邊界層流動是充分發展的，流場上下游的參數具有較好的一致性。 然而使用該模型對復雜地形風電場進行數值模擬時，進出、口廓線分布未達到良好的一致性，給復雜地形的湍流強度評估帶來了困難，影響到數值計算的結果評估[1]。

國內外許多學者對于復雜地形風電場的速度場流動特性進行了大量研究。 Richards P[2]提出了4 點大氣邊界層流動假設和入口邊界條件，采用k-湍流模型較好地模擬了中性大氣邊界層。 在此基礎上，Hargreaves D M[3]通過修改默認的粗糙度，實現了一定工況下的零流向梯度大氣邊界層數值模擬。 Yang Y[4]采用實驗導出的隨高度衰減的湍動能邊界條件，速度場流動計算結果與實驗數據有較好的一致性。 Parente A[5]采用上述類似的方法，通過添加湍動能和耗散率源項，結合基于粗糙長度的壁面函數，對均勻大氣邊界層和方體繞流的流動分離現象進行數值模擬，改進了大氣邊界層的流動模擬方法。 Balogh M[6]再次改進上述源項，在典型山丘地型和真實復雜地形中得到同樣結論。幾位學者從湍流模型參數、增加源項來提高大氣邊界層的自持性平衡屬性，但仍存在入口邊界條件設置問題。 Richards P J[7]除了考慮湍動能的平衡外，還考慮了壓力梯度和剪切應力，結合實驗結果進一步優化了大氣邊界層流動的計算。

部分文獻采用了計算流體動力學（CFD）方法進行復雜地形下風電場尾流場特性的研究。 盧玉華[8]采用雙方程k-模型方法，研究障礙物尺寸對下游山體流場的影響，得到障礙物尺寸距離山體8 倍或10 倍山高時，達到適宜風電場建設需求的湍流強度。 梁思超[9]簡化山體流動問題為方體擾流流動，通過DES模擬方法對擾流機理進行了細致的研究，得到高雷諾數下馬蹄渦系統的漩渦組成和發展機理。在此基礎上，閆姝[10]采用Fluent 軟件模擬了陡峭地形和風電機組前后的擾流現象，發現較大高寬比例誘導更大的山后流動分離，對場區選址有一定的影響。

綜上所述，學者采用不同的研究方法對平坦地形和簡化地形風電場下的的尾跡特性做了大量數值研究，但對于真實的復雜地形研究較少。本文提出一種基于Openfoam 的復雜地形模擬方法，并通過Askervein 山模型進行方法驗證，確認了方法的有效性。 以懷來某復雜地形風電場作為研究對象，快速結構網格建模，得到了不同機位點的速度分布，分析了一定機位處的尾流分布。計算方法為風電場布機設計和控制策略的優化提供了數值參考，對機組穩定運行具有一定的參考意義。

1 幾何模型與數值方法

1.1 幾何模型

Askervein 山位于蘇格蘭一座島嶼上，高度僅為116 m。80 年代初，國際上多家單位合作完成了Askervein 山測風工作，分別進行了239°和210°速度進口的測量，測風速度充分，在國際上被廣泛用于復雜地形模擬研究驗證工作[11]。 本文以Askervein 山為研究對象[圖1（a）]，進行 CFD 方法驗證。 懷來風電場場內布置了 2 座雷達[圖1（b）]，并選取T11，T12，T21 風機點位作為研究點（上述機位點距離雷達點位較近）。 圖1（c）為一個月內懷來風場雷達387 在40～230 m 高程風廓線分布，懷來地形的主風向為西北西和東北東。

圖1 驗證地形和懷來復雜地形實地圖及實際風場雷達風向圖Fig.1 Verification terrain and actual complex terrain and wind direction rose figure at HuaiLai

1.2 數值方法

Askervein 山計算域的設置見文獻[11]，地表粗糙度統一為0.03。 懷來地形的計算域如圖2 所示。 懷來測風的主風向為西北西292.5°和東南東110 °，在 X 方向和 Y 反向分別設置為 21.5 km 和17.6 km，沿地表高度方向設置為10 km，地表粗糙度按照懷來實際地形粗糙度設置。

圖2 懷來地形計算域Fig.2 HuaiLai terrian computational zone

網格劃分采用自主研發的網格生成器繪制三維地形網格，先繪制背景網格，后進行局部加密。 壁面 y+＜30，滿足湍流模型計算要求。 繪制三套懷來地形計算區域的網格，網格數分別為352萬、814 萬和1 452 萬，風機點位附近分辨率分別為 30 m×30 m，20 m×20 m，10 m×10 m（圖3）。從圖3（b），（c）可知，地表分辨率逐漸加大，輪廓更為清晰，貼體網格捕捉更為明顯。

圖3 懷來地形網格Fig.3 Huailai terrain mesh

本文在進行網格無關性計算后，選用814 萬網格。

邊界條件：Askervein 山主風向為210°，測量時氣流穩定，符合中性大氣速度廓線分布，邊界層厚度為500 m，對應位置風速為10 m/s，采用的摩擦速度為0.4 m/s。 在懷來復雜地形的計算中，按照一定扇區入口，以270°和292.5°風速方向為進口方向，入口設置速度分量，上表面壁邊界為自由滑移壁面，下表面為無滑移壁面，出口壓力設置為大氣壓力。入口邊界條件廓線分布采用式（1）進行設置。 在Askervein 山驗證算例中，采用WindSim（WS）類似的邊界條件進行比較，見式（2）。

式中：u*為摩擦速度；κ 為馮卡曼參數，取為0.4；z為近地面高度；z0為壁面粗糙度，在導入的地形文件中包含著粗糙度信息，幾何文件保留粗糙度信息，得到復雜地形下的不均勻特性；LBL，Lobukou分別設置為500 m 和1 000 m。

計算設置：數值計算采用基于Openfoam 開源CFD 工具，通過壓力-速度耦合的SIMPLE 算法求解RANS 方程； 空間離散采用二階精度的中心差分格式；湍流模型采用標準k-ε 模型，修正了湍流模型常數，適應復雜地形的計算要求。

2 結果與分析

本文首先以Askervein 山為研究對象，對研究方法進行驗證； 然后以懷來某復雜風電場地形作為研究對象，經過網格無關性的驗證后，研究了4組風向下的風場速度場和尾流特性分布。

2.1 結果驗證

圖4 為 Askervein 山在 AA-AA 和 B-B 山脊線的加速比分布圖。 加速比為當地風速和入口風速的差值除以入口風速。 圖中：B1，B2 分別為式（1），（2） 的入口邊界條件；M1 和 M2 的網格地表分辨率分別為 133×133×35，200×200×35。 湍流模型參數設置包含 C0，C1 和 C2 3 種（表1）。

圖4 沿著兩個不同的山脊線AA 和B 的加速比分布Fig.4 The profile of speedup along the line AA and B

表1 湍流模型參數設置方案Table 1 Turbulent model parameter setup

WindSim 模擬結果表明，在不同的湍流模型參數和邊界入流條件下，細網格的計算均較好地吻合試驗參數。網格數驗證后，計算了3 種入口配置方案，其中C0 為默認的k-ε 湍流模型參數設置，可知入口邊界條件B1 和湍流模型參數C2 的配置參數模擬效果較好，因此，后續計算均采用B1 和C2 的配置方案。

驗證Askervein 山的計算方案后，對懷來復雜地形的3 種網格配置方案進行驗證。 圖5 為3 套網格對比結果。 在上述工況下，不同網格差別較小，在輪轂高度處速度相差1%左右，均可以滿足網格計算要求。出于計算效率和精度考慮，下面的模擬結果分析均基于841 萬網格。

圖5 三種網格數在兩種扇區下T11，T12 機組的速度廓線分布Fig.5 The velocity profile at T11，T12 wind turbine under two velocity inflow conditions with three grid numbers

2.2 懷來地形不同扇區速度場分析

選取 841 萬地形網格，采用 Openfoam 和Fluent 數值模擬方法，計算來流進口方向分別為正北、正東、正南、正西（順時針方向）的數值模擬，對比云圖和風機點位的風速輪廓線。 圖6 為正北方向來流風電機組附近的流向截面速度云圖。 通過速度邊界層的對比分析可知，Openfoam 近地面的低速回流區域較多。

圖6 正北方向來流在機組附近流向截面速度云圖Fig.6 Velocity sectional contour nearby the wind turbines under the due north inflow

圖7 為兩臺機組在2 個風機點位 T12，T21的速度廓線，其中T12 機組的海拔相對T21 較低，兩臺機組此時處于山體的下游尾跡區。 由圖7可知，Openfoam 的速度廓線受到上游山體的分離流動影響，T12，T21 機位在輪轂高度處 （90 m）的速度相對于Fluent 分別相差了15%和12%，呈現一定的回流現象。

在正東主風向來流風況下，近地層風場受地形影響較為明顯，與山體北側來流匯集為向東運行（圖8）。 由圖8 可知：受地形影響，在小范圍內形成背風測的渦旋型或輻散型流場，流動分離較為明顯；該方向流動的大氣邊界層的自持性較差，流向速度梯度較大，流場需進一步的修正，以便得到更好的模擬結果；Fluent 得到的計算結果在山頂加速較明顯，比Openfoam 計算結果偏高。

圖7 正北方向來流在兩臺機組附近高度方向速度廓線分布Fig.7 Velocity profile nearby the wind turbines under the due north inflow

圖8 正東方向來流在機組附近流向截面速度云圖Fig.8 Velocity sectional contour nearby the wind turbines under due east inflow

圖9 為兩臺機組的速度廓線分布，速度廓線均受到山體背風面流動分離的影響，在輪轂高度處速度與Fluent 分別相差20%和24%，出現一定的回流現象，風能品質下降。 正西和正南方向風電機組都處于山體前方，同平坦地形流場類似，自持性較好。

通過上述結果發現：采用Openfoam 和Fluent模擬方法能夠比較合理地模擬出復雜地形上的近地層風場流動變化；受到局部不均勻粗糙度影響，Openfoam 邊界層發展較厚。

圖9 正東方向來流在兩臺機組附近高度方向速度廓線分布Fig.9 Velocity profile nearby the wind turbines under due east inflow

3 結論

本文基于開源的Openfoam 數值軟件進行了二次開發，采用自主研發的網格建模工具生成適應地形的帖體網格，通過Askervein 山的數值模擬驗證了計算方法的有效性。 以懷來復雜地形風電場為研究對象，在考慮地形粗糙度的前提下，研究了4 組來流工況下的復雜地形數值模擬，分析了對應風向下的風速和尾流特性的關系，得到以下結論。

①復雜地形風電場的流場受到地形起伏影響較大，當風向經過山體后再經過風電場時，會使得風場整體風能品質下降，表現為流場的速度較低，尾流損失較大。

②不同扇區的入流計算表明： 在復雜地形影響下，流場繞流現象較為復雜，特別是低速流動時； 在輪轂高度處，各個機位點的速度與Fluent計算結果對比，相差15%左右，主要是受到山后回流的影響。入口邊界條件修正、湍流模型參數修正等需要進一步完善。