快速检索
  气象   2012, Vol. 38 Issue (11): 1378-1389.  

论文

引用本文 [复制中英文]

方艳莹, 徐海明, 朱蓉, 等, 2012. 基于WRF和CFD软件结合的风能资源数值模拟试验研究[J]. 气象, 38(11): 1378-1389. DOI: .
[复制中文]
FANG Yanying, XU Haiming, ZHU Rong, et al, 2012. Study on Numerical Simulation of Wind Energy Resources Based on WRF and CFD Models[J]. Meteorological Monthly, 38(11): 1378-1389. DOI: .
[复制英文]

资助项目

科技部国际科技交流与合作专项 (2010DFA62830) 和科技部863课题 (2007AA05Z424) 共同资助

第一作者

方艳莹,主要从事风能资源数值模拟研究。

通信作者

朱蓉,主要从事大气边界层和风能资源数值模拟研究.Email:rongzhu@cma.gov.cn

文章历史

2011年11月01日收稿
2012年4月18日收修定稿
基于WRF和CFD软件结合的风能资源数值模拟试验研究
方艳莹 1,2, 徐海明 2, 朱蓉 1, 王鹏 1,3, 何晓凤 1, Didier Delaunay 4, 付斌 4, 王黎 4    
1. 中国气象局公共气象服务中心,北京 100081
2. 南京信息工程大学,南京 210044
3. 中国人民解放军理工大学,南京 210007
4. 法国美迪有限公司, 北京 100027
摘要:运用中尺度数值模式WRF与法国CFD软件Meteodyn WT相结合的方法 (WRF/WT),进行了广东省海陵岛地区的水平分辨率100 m×100 m的风能资源数值模拟试验,采用海陵岛上7座测风塔观测资料对WRF/WT模式的模拟风场进行误差检验,并与WRF/WAsP模式系统对单点风能参数模拟误差进行对比,研究WRF/WT模式系统在风电场微观选址和分散式风电开发利用中应用的可行性。结果表明:中尺度模式与CFD软件结合的数值模拟方法对区域风能资源分布趋势的模拟比单纯应用CFD软件更准确;WRF/WT模式系统应用于复杂地形风能资源数值模拟评估是可行的,其对区域风能资源参数分布模拟的准确率与WRF/WAsP对2 km范围内风能资源参数模拟的准确率相当;WRF/WT模式系统在风速频率分布不满足Weibull分布的情况下和陡峭地形条件下有较好的模拟效果,相对WRF/WAsP有明显优势。今后需进一步研究中尺度模式与CFD软件的衔接方法,以及对中尺度模式模拟结果的误差订正。
关键词中尺度数值模式    CFD模式    复杂地形    风电场微观选址    分散式风电开发    
Study on Numerical Simulation of Wind Energy Resources Based on WRF and CFD Models
FANG Yanying1,2, XU Haiming2, ZHU Rong1, WANG Peng1,3, HE Xiaofeng1, Didier Delaunay4, FU Bin4, WANG Li4    
1. Public Meteorological Service Centre, China Meteorological Administration, Beijing 100081;
2. Nanjing University of Information Science and Technology, Nanjing 210044;
3. PLA University of Science and Technology, Nanjing 210007;
4. Meteodyn China, Beijing 100027
Abstract: A combined model system (WRF/WT) of the mesoscale model WRF and the Meteodyn WT, a CFD model from France, was carried out the numerical simulation experiments of regional wind resources over Hailing Island of Guangdong Province with a horizontal resolution of 100 m × 100 m. The observational data from 7 wind towers in Hailing Island were used to test the results modeled by WRF/WT, and compared with the simulated errors to the single-point wind parameters, thus studying the feasibility of WRF/WT model system in micro-siting for wind farm and the application of distributed development of wind power and utilization. The results showed that, the combined model system of mesoscale model and CFD model in simulating the trends of the regional wind energy resource distribution is more accurate than CFD model used only; the WRF/WT model system used in complex terrain for wind resource numerical simulated evaluation is feasible, and its accuracy in simulating the parameters of the regional wind energy resource distribution is close to the result of WRF/WAsP simulating within 2 km. The WRF/WT model system compared to WRF/WAsP has obvious advantages when the frequency distribution does not meet the Weibull distribution or in a steep terrain. In the future, we need to further research in the interface methods of mesoscale model and CFD models, and the error correction of mesoscale model simulations.
Key words: mesoscale numerical model    CFD model    complex terrain    micro-siting for wind farm    distributed development of wind power    
引言

风能资源评估技术自20世纪70年代以来,经历了由基于观测资料统计分析的评估方法到应用数值模拟技术的发展历程,目前在国内外风能资源数值模拟已成为风能资源评估的主要技术手段[1-2]

20世纪90年代,丹麦Risoe国家实验室在Jackson he Hunt理论基础上,发展了一个用于风电场微观选址的资源分析工具软件——WAsP (Wind Atlas Analysis and Application Program)[3],在随后的10年中,又经过不断地改进[4-5]。进入21世纪,美国TrueWind Solutions公司在应用数值模式评估风能资源方面处于国际领先地位[6-7],其产品MesoMap和SiteWind风能资源评估系统在20多个国家和地区应用于风能资源评估[8-9]。澳大利亚科学院 (CSIRO) 也发展了类似的线性和非线性的小尺度风场模型 (WindScape)[10]。日本使用美国大气边界层模式RAMS (Regional Atmospheric Modeling System) 也开展了本国的高分辨率的风能资源数值模拟[11]。近年来,计算流体力学模式 (Computational Fluid Dynamics,CFD) 开始被用来模拟复杂地形风场,为风电场微观选址服务,如挪威的WindSim和法国的Meteodyn WT等[12-14]。我国也开展了一些CFD用于复杂地形风能资源评估的研究[15-16]

中尺度数值天气预报模式与小尺度数值模式相结合的风能资源数值模拟系统已被广泛地用于高分辨率的风能资源评估。在联合国环境规划署 (NUEP) 的SWERA项目中,美国国家可再生能源实验室用中尺度模式MASS与小尺度模式MesMap结合的模式系统完成了美国、巴西、中国东部等国家和地区水平分辨率1 km×1 km的风能资源模拟[17-18];法国能源部和环境部采用CFD模式Meteodyn WT结合美国中尺度模式RAMS对法国Provence-Alpes-Cote-d’Azur地区进行风能资源评估[19];王鹏等[20]采用WRF和小尺度动力诊断模式CALMET模式相结合的模式系统,对广东省海陵岛地区进行了风能资源数值模拟;杨罡等[21]在鄱阳湖地区对大气边界层特征进行过数值模拟研究,何晓凤等[22]采用MM5与WindSim结合进行了鄱阳湖区域风能资源的数值模拟研究。

当前,鉴于我国电网消纳能力对大规模风电开发的制约,国家能源局发出了因地制宜、开展分散式风电开发的通知。零散的可利用风能资源多位于复杂地形地区,分布点众多且覆盖面积小,不可能用设立测风塔的方式进行风能资源评估,中尺度模式的模拟精度也不够。因此,本文将采用中尺度模式WRF和法国CFD软件Meteodyn WT相结合的方法 (以下简写WRF/WT) 对复杂地形的风能资源进行数值模拟实验,并与WRF与丹麦微尺度风能资源模拟软件WAsP结合的模式系统 (以下简称WRF/WAsP) 的模拟结果进行对比,探索服务于分散式风电开发的风能资源数值模拟方法。

1 数值模式系统

Weather Research and Forecast (WRF) 模式系统是由美国国家大气研究中心 (NCAR),国家海洋与大气管理局[国家环境预报中心 (NCEP) 及预报系统实验室 (FSL)、空军气象局 (AFWA),海军研究实验室,奥克拉荷马大学及联邦航空管理局 (FAA)]共同参与开发的新一代中尺度预报模式和同化系统。该模式采用高度模块化、并行化和分层设计技术,集成了迄今为止在中尺度方面的研究成果。模拟和实时预报试验表明,WRF模式系统在预报各种天气中都具有较好的性能,具有广阔的应用前景,是目前比较成熟的数值预报模式之一。

Meteodyn WT是法国美迪公司开发的一套计算流体力学 (CFD) 模式[23-24]。它以大气湍流动量守恒方程和质量守恒方程 (N-S方程) 为动力框架,在定常不可压的情况下,方程变为

$\frac{{\partial U}}{{\partial x}} + \frac{{\partial V}}{{\partial y}} + \frac{{\partial W}}{{\partial \mathcal{z}}} = 0$ (1)
$\begin{array}{l} U\frac{{\partial U}}{{\partial x}} + V\frac{{\partial V}}{{\partial y}} + W\frac{{\partial W}}{{\partial \mathcal{z}}} = \\ - \frac{1}{\rho }\frac{{\partial P}}{{\partial x}} + \frac{{\partial \overline {{u'^{2}}} }}{{\partial x}} + \frac{{\partial \overline {{u'}{v'}} }}{{\partial y}} + \frac{{\partial \overline {{u'}{w'}} }}{{\partial \mathcal{z}}} \end{array}$ (2)
$\begin{array}{l} U\frac{{\partial V}}{{\partial x}} + V\frac{{\partial V}}{{\partial y}} + W\frac{{\partial V}}{{\partial \mathcal{z}}} = \\ - \frac{1}{\rho }\frac{{\partial P}}{{\partial y}} + \frac{{\partial \overline {{u'}{v'}} }}{{\partial x}} + \frac{{\partial \overline {{u'^{2}}} }}{{\partial y}} + \frac{{\partial \overline {{u'}{w'}} }}{{\partial \mathcal{z}}} \end{array}$ (3)
$\begin{array}{l} U\frac{{\partial W}}{{\partial x}} + V\frac{{\partial W}}{{\partial y}} + W\frac{{\partial W}}{{\partial \mathcal{z}}} = \\ - \frac{1}{\rho }\frac{{\partial P}}{{\partial \mathcal{z}}} + \frac{{\partial \overline {{u'}{w'}} }}{{\partial x}} + \frac{{\partial \overline {{u'}{w'}} }}{{\partial y}} + \frac{{\partial \overline {{w'^{2}}} }}{{\partial \mathcal{z}}} \end{array}$ (4)

利用K理论对 $\overline {{u'}{v'}} ,\overline {{u'}{w'}} ,\overline {{u'}{w'},} \overline {{u'^{2}}} ,\overline {{v'^{2}}} ,\overline {{m'^{2}}} $ 等湍流通量用转换成梯度输送的方式进行参数化,即

$\begin{array}{l} - \overline {{u'}{v'}} = - {v_T}\frac{{\partial U}}{{\partial y}} - \overline {{u'}{w'}} = - {v_T}\frac{{\partial U}}{{\partial z}}\\ - \overline {{u'}{w'}} = 0 \end{array}$ (5)

引入以下方程:

$\begin{array}{l} {U_j}\frac{{\partial k}}{{\partial {x_j}}} = {P_\lambda } - \varepsilon + \frac{\partial }{{\partial {x_j}}}\left[ {(\frac{{{v_T}}}{{{\sigma _k}}})\frac{{\partial k}}{{\partial {x_j}}}} \right]\\ 其中\left\{ \begin{array}{l} \varepsilon = {C_\mu }\frac{{{v_T}}}{{{L_T}}}k\\ {P_\lambda } = {v_T}\left( {\frac{{\partial {U_i}}}{{\partial {x_j}}} + \frac{{\partial {U_j}}}{{\partial {x_i}}}} \right)\frac{{\partial {U_j}}}{{\partial {x_j}}}\\ {v_T} = {k^{1/2}}{L_T} \end{array} \right. \end{array}$ (6)

根据Yamada[25]和Arritt[26]参数化方法,考虑了热力层结中的稳定或不稳定情况下,计算湍流尺度:

${L_T} = \sqrt 2 S_m^{3/2}l\left\{ \begin{array}{l} \frac{1}{l} = \left( {\frac{1}{{{l_0}}} + \frac{1}{{k\mathcal{z}}}} \right),\mathcal{z} = 高度\\ {C_\mu } = \frac{{4{S_m}}}{{{B_1}}}\\ {S_m} = \left\{ {\begin{array}{*{20}{l}} {1.96\frac{{(0.1912 - {R_{if}})(0.2341 - {R_{if}})}}{{(1 - {R_{if}})(0.2231 - {R_{if}})}},}&{{R_{if}} < 0.16}\\ {0.085,}&{{R_{if}} \ge 0.16} \end{array}} \right.\\ {B_1} = 16.6\\ {l_0} = 100 m\\ k = 0.41 \end{array} \right.$ (7)

其中Rif为通量Richarson数;k为卡门常数;LT是湍流混合长度;l0是一个临界值,当空间高度很大时,使l的值近似为一个常数;Cμ是湍流动能方程里湍流动能转化为热内能的部分;B1是从风动实验得到的经验性参数。法国Meteodyn WT采用分块结构对整个空间进行网格化。它在用户所定义的结果点 (即用户关心的地点) 处进行笛卡尔网格加密,网格单元维度的扩展以及长宽比受到控制,以避免收敛的不稳定性。网格不会自动根据地形来选择加密或者稀疏,如果想在一些变化复杂的地区提高准确度,需要在此区域附近定义结果点,网格则会在此处加密。如此的网格分布方式将显著地降低计算规模,图 1显示了本文研究区域的网格加密情况。

图 1 法国WT网格生成方法,(a) 图中圆点表示本文定义的结果点,(b) 在结果点加密以后的计算区域网格分布 Fig. 1 The method of mesh generation of Meteodyn WT (a) the dots represent the defined points, and (b) the mesh distribution after the intensifying defined points in computational domain

计算区域网格生成以后,采用与目前CFD软件常用的算法。按照16个等分的风向角,假设入口风速10 m·s-1,根据中性风速廓线分别进行定向模拟计算,由此得到标准入口风速、不同风向条件下的风场分布,即定向计算结果。这时如果有测风塔观测资料的话,将根据测风塔观测结果与计算区域内相同位置定向计算结果之间的统计关系,推算用户关心的结果点处的风速值和风能参数。Meteodyn WT模式可以允许输入多个测风塔的观测资料,模拟出整个区域的风能资源分布。

WAsP是由丹麦国家实验室 (Risoe) 开发出来的微尺度风能资源数值模拟软件。该软件通过测风塔观测资料计算地转风,假定一定范围内地转风速不会变化,再根据周围地区的地表特征,由地转风推算任意点的近地层风速和风能参数。由于是采用一个观测点外推,误差会随着与测风点距离的增加而加大。因此,不建议将WAsP直接用于较大区域的风能资源评估。另外,WAsP本身采用线性模型计算方法,有其一定的局限性,它会随着被计算流体经地形的复杂而带来计算结果的不确定性。所以WAsP对地形相对简单、地势较平坦的地区较为适用,但对较复杂地形,由于受许多边界条件等的限制,不很适合采用该软件[27]

WAsP是利用对数公式,假设在平坦地形下,风速较大时的风廓线为:

$u(\mathcal{z}) = \frac{{{u_*}}}{k}\ln \frac{\mathcal{z}}{{{\mathcal{z}_0}}}$ (8)

式中, $\mathcal{z}$ $\mathcal{z}$ 高度上的风速, $\mathcal{z}_0$ 是地面粗糙度,κ是卡曼常数,一般取值0.40,u*是摩擦速度,与地面压力有关:

$\left| \tau \right| = \rho u_*^2$ (9)

式中,ρ是空气密度。但即使是在正常的风速情况下,当 $\mathcal{z}$ 超出边界值几十米,风速也会偏离风廓线。发生偏离是由于浮力的动力扰动造成的,地面粗糙度不再是地面特征的唯一影响因素,地表热通量也会成为另一个影响地面特征的因子。在晚上,地面冷却时,扰动减少,导致风廓线随高度快速增长;相反,白天的加热作用使扰动增强,风廓线随高度基本不变。廓线表达式可概括为:

$u(\mathcal{z}) = \frac{{{u_*}}}{k}[\ln (\mathcal{z}/{\mathcal{z}_0}) -ψ(\mathcal{z}/L)]$ (10)

式中,ψ是经验函数,L是新引进的Monin-Obukhov长度:

$L = \frac{{{T_0}}}{{kg}}\frac{{{c_p}u_*^3}}{{{H_0}}}$ (11)

式中,T0H分别是地面的绝对温度和热通量,cP是常压下空气的热容量,g是重力加速度。

大气边界层的风主要由天气尺度系统造成气压不同而形成,如高压和低压系统的过境,当受到气压改变产生的强迫作用时,边界层结构也会发生变化,气压梯度力和摩擦力会达到一个适当的平衡状态,这个平衡是建立在静力稳定、均匀层结和正压的条件下,这个结果通常称为地转曳力定律,摩擦速度u*和地转风G之间的关系为:

$G = \frac{{{u_*}}}{k}\sqrt {{{\left[ {\ln (\frac{{{u_*}}}{{f{\mathcal{z}_0}}}) - A} \right]}^2} + {B^2}} $ (12)

WAsP的基本思路是,根据测风塔的测风数据推算出地转风,假设地转风在水平方向是不变的,因此测风塔周围任一点的风能参数就可由地转风和地表粗糙度推算出来。

2 数值模拟试验方案设计

本文选择广东省阳江市海陵岛为数值模拟试验区。海陵岛四面环海,总面积约108 km2,海岸线长约142 km。岛上多丘陵山地 (图 2),西部和中部山地高度超过300 m,东部山地不足200 m。岛上地势较低的地方以农田和居住区为主,植被较少,地势较高的山地上以大片的树林为主,另外全岛还分布有零星的湖泊和池塘。2003年9月1日至2004年8月31日广东省气候中心在海陵岛上开展了风能观测实验,设置了7座测风塔 (图 3)。其中1号和2号塔位于海陵岛西南侧接近山顶的地方,海拔高度201 m左右;3号塔位于海陵岛西北,高度239 m左右的山体顶部;4号塔位于海陵岛的北部海滩上;5号塔位于海陵岛对岸的大陆沿岸,海拔高度5 m;6号塔则位于海陵岛西部山体北侧的山腰上;7号塔处于海陵岛西部山体南侧的山坳中。各测风塔的观测高度和海拔高度见表 1

图 2 海陵岛地形图 Fig. 2 The topographic map of Hailing Island

图 3 Meteodyn WT模拟范围 (圆圈内) 和WRF模式网格点分布 (黑点) 以及测风塔分布 (黄色标记) Fig. 3 Simulation domain of Meteodyn WT (within circle), the mesh distribution of WRF (black points), and the distribution of observation towers (yellow marked)

表 1 测风塔海拔高度和观测高度 Table 1 The elevation above sea level and the observed height of towers

WRF模式采用三重嵌套网格,网格距分别为27、9和3 km;地表类型资料采用美国地质测量局USGS发布的30″全球资料,27和9 km网格地形资料采用USGS的10′数据,3 km网格采用美国SRTM3数据,分辨率为3″。同化全球环流模式再分析资料 (NCEP) 和常规地面和探空气象资料,生成WRF模式初值场和侧边界,边界层物理过程参数化采用Mellor-Yamada-Janjic方案。模拟结果逐小时输出,作为小尺度模式的输入数据。

法国Meteodyn WT模式水平分辨率设为100 m×100 m,垂直分辨率为8 m,水平扩展系数为1.1,垂直扩展系数为1.2,计算区域半径为18 km。图 3中的圆圈表示法国Meteodyn WT模式的计算范围,圆圈内的黑点表示WRF模式第三重、格距3 km的网格点,地形数据采用美国SRTM的30 m分辨率的数据,土地利用分类借助google earth工具手工设定。

在Meteodyn WT模式的计算范围内共有121个WRF模式网格点,将这121个网格点的数值模拟风速作为观测风速输入Meteodyn WT模式,Meteodyn WT模式具有综合分析计算范围内多个测风塔数据的功能。由于WAsP接受一个测风塔数据的输入和分析,所以将与实际测风塔距离最近的WRF网格点模拟风速输入WAsP,计算分析实际测风塔位置处的风能参数,再与实际测风数据进行对比。

3 数值模拟结果与分析 3.1 海陵岛风能资源数值模拟

Meteodyn WT模式首先在大气中性层结条件下进行16个方向的定向计算,得到不考虑天气背景和环境风场、不同来向的气流在地形和地表作用下产生的风速加速因子和湍流强度的空间分布。考虑到3号塔位于海拔209 m的山顶且周围地形平坦,认为其相对于其他6个测风塔来说,对海陵岛所处的区域环境风场的代表性更好一些。因此,采用3号塔60 m高度观测的2003年9月1日至2004年8月31日的风向频率分布,对定向计算得到的16个方位来向所产生的风速加速因子和湍流强度进行加权平均, 得到了风速加速因子和湍流强度分布 (图 4)。从图 4可以看出,在山顶区域风加速因子数值较大,湍流强度较弱。由于海陵岛全年的盛行风向以东北风和偏东风为主,因此在较高地形的西南侧和西侧,湍流强度较其他地方强,这是背风坡湍涡的作用。整个海陵岛基本呈东—西走向,但略偏一点西南—东北走向,北偏西的风向使海陵岛的迎风面最大,由于岛上山地较多,造成了岛的南偏东海域上湍流强度较大。海陵岛的最西端地势平坦,因此湍流强度不大。

图 4 Meteodyn WT模式的定向计算结果 (高度60 m,水平分辨率100 m×100 m)(a) 风加速因子分布,(b) 湍流强度分布 (图中右下角为3号塔年平均风玫瑰图) Fig. 4 The oriented computing results of Meteodyn WT, with height 60 m, and hori\mathcal{z}ontal resolution 100 m×100 m (a) the distribution of wind accelerating factor, and (b) the distribution of turbulence intensity

由于中尺度模式网格点上的计算值表示的是网格平均值,对于3 km的格距,相当于代表 9 km2范围的平均值。海陵岛地形复杂,离地面越近,风速分布就越不均匀。采用Meteodyn WT模式计算范围内121个WRF模式网格点的2003年9月1日至2004年8月31日逐小时输出的100 m高度风向和风速输入Meteodyn WT模式,对定向计算结果进行统计分析后,即可得到水平分辨率100 m×100 m的平均风速和风功率密度分布 (图 5c5d)。为了分析中尺度模拟结果的作用,本文将海陵岛7个测风塔2003年9月1日至2004年8月31日的观测数据输入Meteodyn WT模式,同样得到了水平分辨率100 m×100 m的平均风速和风功率密度分布 (图 5a5b)。

图 5 海陵岛100 m高度上2003年9月1日至2004年8月31日平均风速和风功率密度分布 (a) 和 (b) 分别为Meteodyn WT模式根据测风塔观测资料模拟的年平均风速和风功率密度分布,(c) 和 (d) 分别为Meteodyn WT模式根据WRF网格点上输出结果模拟的年平均风速和风功率密度分布 Fig. 5 The mean wind speed and power density of Hailing Island at 100 m height from 1 September 2003 to 31 August 2004. The mean wind speed (a) and power density (b) are modeled by Meteodyn WT according to the observating data of towers, and the mean wind speed (c) and power density (d) are modeled by Meteodyn WT according to the results of WRF

图 5可以看出,WRF/WT模式系统得到的年平均风速和风功率密度值总体大于Meteodyn WT模式根据测风塔资料计算的结果,说明WRF的近地层风速模拟结果值偏大;两种方法模拟出的岛内风能资源最大值出现的位置是相同的,都是位于山区地形较高处。

图 5图 6[28]对照来看,图 5c图 5d整体风能资源的分布形势与广东沿海及近海风能资源的分布形势一致,都是从海上到陆上风速有明显减小的变化规律。这是由于中尺度模式WRF能够模拟出海上到陆上风速的变化规律,输入进WT模式中的121个WRF网格点资料中包含的中尺度风场背景信息,因此图 5c图 5d中海陵岛上面向外海一侧的风能资源明显大于面向内陆一侧。

图 6 广东沿海及近海100 m高度年平均风速分布 (单位:m·s-1,引自文献[28]) Fig. 6 The distribution of average wind speed off China and the Guangdong coast at 100 m height (unit: m·s-1, from Ref.[28])

很显然,图 5a图 5b中的风能资源分布与环境背景场没有关系。这是由于海陵岛东西向约20 km、南北向约12 km,仅仅7个测风塔不能把握全岛范围的风速分布,况且这7个测风塔均分布在面向内陆一侧和岛的东西两端,因此对于类似海陵岛这样尺度10~20 km的复杂地形精细化风能资源评估,更适合采用中尺度模式与CFD模式相结合的风能资源数值模拟方法。由于中尺度数值模拟误差会影响CFD数值模拟准确率,因此,如果能根据实际测风资料对中尺度模拟结果进行订正后再输入CFD模式,定会明显减小模拟误差。

3.2 模拟结果误差分析

分析风能资源数值模拟的准确率时,除了检验平均风速以外,风速频率分布也是一个重要因素。因为风力机是在风速达到一定大小 (切入风速) 时才开始发电,随着风速的增大,发电量也增大,当风速大到一定值 (额定风速) 时,发电量达到满发,之后风速如何增加,发电量也不会变化。图 7给出了一个500 kW风力机发电量与风速的关系[29]。因此,发电量与风速在各个风速值区间中出现的时间长短有关,即一段时间内的发电量等于该时段内所有出现的各档风速所对应输出功率乘以该风速出现时间的总和。本文将根据图 7给出的风力机输出功率曲线计算发电量,以此检验模式系统对风速分布频率模拟的准确率。

图 7 发电量与风速的关系 Fig. 7 The relationship of electric energy production and wind speed

为了区分中尺度模式WRF与Meteodyn WT模式各自对WRF/WT模拟误差的贡献,本文采用代表性较好的3号测风塔资料,输入Meteodyn WT模式进行模拟计算,然后对其余6个测风塔位置上2003年10月和2004年1、4、7月的各月平均风速和发电量模拟值进行检验 (图 8),得到月平均风速和发电量的模拟值与实测值之间的相关系数RR的绝对值越大,表示两个量之间的相关性程度越高。从图 8可以看出模拟风速值和发电量的模拟值与实测值的相关系数分别为0.80和0.77,平均相对误差分别为6.6%和19.1%。同时,在Meteodyn WT模式根据121个WRF网格点数值模拟结果计算得到的风能资源分布参数的基础上,计算7个测风塔位置和各自观测高度上2003年10月和2004年1、4、7月的各月平均风速和发电量,然后与实测月平均风速和发电量进行对比检验 (图 9a9b)。结果表明WRF/WT得到的月平均风速和发电量模拟值与实测值的相关系数分别为0.86和0.79,平均相对误差分别为9.79%和31.68%。为了研究WRF/WT应用于复杂地形区域风能资源评估的可行性,本文还检验WRF/WAsP系统对测风塔位置上的风能参数模拟准确率,将与测风塔距离最近的WRF网格点上的逐小时模拟风速输入到WAsP中,计算测风塔位置上的月平均风速和发电量 (图 9c9d)。结果表明WRF/WAsP得到的月平均风速和发电量模拟值与实测值的相关系数分别为0.90和0.87,平均相对误差分别为8.08%和27.86%。

图 8 采用3号测风塔资料的Meteodyn WT对其余6个测风塔位置风能参数模拟值的检验结果 Fig. 8 The test results of rest 6 towers'wind parameters modeled by Meteodyn WT according to the data of No.3 tower

图 9 WRF/WT和WRF/WAsP模式系统的风能参数数值模拟结果的检验 Fig. 9 The test of the wind parameters-modeled results by WRF/WT (a, b) and WRF/WAsP (c, d)

从以上WRF/WT和WRF/WAsP的检验结果对比可以看出,模拟值与实测值的相关都非常好,WRF/WAsP的相关系数较大。这是由于WRF/WAsP并不进行大面积的计算,只是用距离测风塔最近的WRF网格点上的模拟结果推算测风塔处的风能参数。由于WRF模式最内重网格距是3 km,7个测风塔与最近WRF网格点的距离从0.2~1.7 km不等,因此WRF/WAsP模拟与实测的相关性检验结果完全依赖WRF的模拟能力。另一方面从图 9可见,总体来讲WRF模拟的平均风速要明显高于实际观测值,由此导致了WRF/WT和WRF/WAsP模拟得到的月平均风速和发电量的相对误差要高于图 8中Meteodyn WT的结果,其中发电量相对误差更大一些。这是因为发电量的模拟误差依赖对每个风速段出现时间频率的模拟准确率,这对模式的模拟能力要求更高。

以下通过个例分析进一步研究WRF/WT的模拟能力。2003年10月,6号塔WRF/WT模拟的平均风速相对误差是7.2%,WRF/WAsP则为0.2%,明显低于WRF/WT;但是对发电量的模拟,WRF/WT模拟误差为0.2%,反而远远小于WRF/WAsP的12.8%。图 10是2003年10月6号塔观测的和各个模式模拟的风速频率分布频率,可以看出6号塔的实测风速频率分布并不完全符合Weibull分布,风速频率分布呈一大一小双峰形态。因为WAsP并不是直接模拟风速频率分布,而是假定风速频率符合Weibull分布,通过求取Weibull参数的A和K值,间接得到风速频率分布,所以无法模拟出风速频率不符合Weibull分布的情况。在此个例中,WRF/WT成功地模拟出了风速频率一大一小双峰分布的形态,因此发电量的模拟误差就很小。这也说明,平均风速误差小,不一定代表发电量的预报也一定准确。

图 10 2003年10月6号塔不同模式系统模拟的风速频率分布 Fig. 10 The No.6 tower's distributions of wind speed modeled by different models in October 2003(a) observed, (b) by WRF/WT, (c) by WRF/WAsP and (d) by WRF

2003年10月,7号塔处WRF/WT和WRF/WAsP模拟的平均风速相对误差分别为0.3%和9.6%;对于发电量的模拟,WRF/WT和WRF/WAsP的模拟相对误差分别为12.9%和31.5%,可见WRF/WT的模拟效果均好于WRF/WAsP。这可能是由于7号塔位于海陵岛西部高度172 m山体南侧的山坳里,塔基座海拔高度49 m,只有东南方向面对大海、地势开阔,其他方向均被山体阻挡 (图 11)。2003年10月的主导风向为NNE,此时7号塔就刚好处于山的背风坡,气流过山后,产生了较大的湍流强度。WAsP是线性模式,更适用于平坦地形;而Meteodyn WT基于N-S方程采用有限元方法求解,充分考虑了近地层大气湍流的作用,并适用于陡峭地形,所以此时WRF/WT的模拟效果就明显优于WRF/WAsP。

图 11 7号塔周边地形及NNE方向定向计算的湍流强度分布结果 (图中黑色三角表示7号塔的位置)(a) 湍流强度分布,(b) 地形高度分布 Fig. 11 The topography around No.7 tower and the distribution of turbulence intensity oriented computed at NNE (a) the distribution of turbulence intensity, and (b) the topography (Dark triangle represents the location of No.7 tower)

2003年10月,1号塔处WRF/WT和WRF/WAsP模拟的平均风速相对误差分别为31.9%和23.4%,发电量模拟结果相对误差分别为127.0%和107.3%,两者的模拟误差都非常大。图 12为2003年10月1号塔处不同模式系统模拟的风向频率分布,可以看到1号塔实际主导风向以偏东风为主,但WRF模拟的主导风向则为偏北风。说明对于复杂地形的地区,来流方向不同,所产生的风能资源参数分布的差异是非常大的,因此WRF/WT模式系统用于复杂地形风能资源数值模拟时,要保证WRF模式对主导风向预报的准确率。

图 12 2003年10月1号塔处不同模式系统模拟的风向频率分布 Fig. 12 The No.1 tower's distributions of wind direction modeled by different models in October 2003(a) observed, (b) by WRF/WT, (c) by WRF/WAsP and (d) by WRF
4 结论与讨论

本文运用中尺度数值模式与CFD软件相结合的方法,在广东海陵岛地区开展复杂地形风能资源数值模拟研究,探讨针对风电场微观选址和分散式风电开发的风能资源数值模拟方法。通过采用海陵岛上7座测风塔观测资料对WRF/WT模式模拟风场的误差检验,以及与WRF/WAsP模式系统对单点风能参数模拟误差的对比,得到结论如下:

(1) 中尺度模式与CFD软件结合的模式系统应用于复杂地形风能资源数值模拟评估是可行的,可以在没有测风塔观测的条件下,得到高分辨率的区域风能资源分布,为风电场微观选址和分散式风电开发的前期工作提供科学依据。

(2) 中尺度模式与CFD软件结合的模式系统对区域风能资源分布趋势的模拟比单纯应用CFD模式更准确,因为中尺度模式为CFD模拟提供了环境背景风场。

(3) WRF/WT模式系统对较大区域风能资源参数分布模拟的准确率与WRF/WAsP对2 km范围内风能资源参数模拟的准确率相当。

(4) 相对WRF/WAsP,WRF/WT模式系统在风速频率分布不满足Weibull分布的情况下和陡峭地形条件下有明显优势。

在中尺度数值模拟的基础上采用CFD软件进行微尺度风场数值模拟,首先需要对中尺度数值模拟误差有一定的评估或进行必要的订正,怎样对中尺度数值模拟结果进行订正,尤其是在没有测风塔的情况下根据历史气象观测资料进行订正,是需要进一步研究的问题。此外,测风塔观测风速与中尺度数值模拟风速没有直接的可比性,中尺度模式网格点上的风速代表网格平均值,测风塔观测值只代表测点周围区域,依赖于地形的复杂程度,因此简单地将中尺度网格点作为测风塔放入CFD中进行数值模拟也是不恰当的。本文考虑到离地面越远,受下垫面影响越小,风速所代表的范围就越大,认为垂直层数越高的风速模拟值与实际风速值的代表范围越接近,所以采用WRF模式100m高度的风速模拟结果代入Meteodyn WT模式中。但是,采用多少高度的中尺度模拟结果合适,或者这种中尺度模式与CFD软件的方式是否合适,也是需要进一步研究的问题。

致谢:感谢宋丽莉研究员为本文提供广东海陵岛测风资料。

参考文献
李泽椿, 朱蓉, 何晓凤, 等, 2007. 风能资源评估技术方法研究[J]. 气象学报, 65(5): 708-717. DOI:10.11676/qxxb2007.066
Zhu Rong, Zhang De, Wang Yuedong, et al, 2009. Assessment of Wind Energy Potential in China[J]. Engineering Sciences, 7(2): 18-26.
Ib Troen, Erik L P. European Wind Atlas[R]. Risoe National Laboratory, 1989. https://rd.springer.com/content/pdf/10.1007%2FBF00122090.pdf
Mortensen N G, Bowen A J, Antoniou I. Improving WAsP predictions in (too) complex terrain.//Proceedings (online). 2006 European Wind Energy Conference and Exhibition, Athens (GR), 27 Feb-2 Mar 2006. (European Wind Energy Association, Brussels, 2006: 9.
Badger J, Mortensen N G, Hansen J C. The KAMM/WAsP numerical wind atlas-a powerful ingredient for wind energy planning.//Proceedings (CD-ROM). The Great Wall world renewable energy forum and exhibition 2006, Beijing (CN), 23-27 Oct 2007.
Brower M C, et al.2001, Applications and validations of the MesoMap wind mapping system in different climatic regimes[C]. Proc Amer Wind Energy Assoc, Windpower 2001, 6.
Ayotte K W, Davy R J, Coppin P A, 2001. A simple temporal and spatial analysis of flow in complex terrain in the context of wind energy modelling[J]. Boundary-Layer Meteorology, 98: 275-295. DOI:10.1023/A:1026583021740
Manwell J F, Rogers A L, McGowan J G, et al, 2011. An offshore wind resource assessment study for New England[J]. Renewable Energy, 27(2): 175-187.
Elliott D, Schwartz M. Development and validation of high-resolution state wind resource maps for the United States[R]. National Renewable Energy Laboratory. Technical Report NREL/TP-500-38127, July 2005.
Coppin P A, Ayotte K A, Steggel N. Wind resource assessment in Australia-a planners guide, Wind Energy Research Unit, CSIRO Land and Water, 2003.
Atsushi Yamaguchi, Takeshi Ishihara, Yozo Fujino. An assessment of offshore wind energy potential using mesoscale model and GIS[R]. European Wind Energy Conference & Exhibition, 2004, 22-25 November, London, UK.
Serozhah Milashuk, William A, 2011. Crane. Wind speed prediction accuracy and expected errors of RANS equations in low relief inland terrain for wind resource assessment purposes[J]. Environmental Modelling & Software, 26(4): 429-433.
Palma J M L M, Castro F A, Rodrigues A H, et al, 2008. Linear and nonlinear models in wind resource assessment and wind turbine micro-siting in complex terrain[J]. Journal of Wind Engineering and Industrial Aerodynamics, 96(12): 2308-2326. DOI:10.1016/j.jweia.2008.03.012
Hwang Y, Paek I, Yoon K, et al, 2010. Application of wind data from automated weather stations to wind resources estimation in Korea[J]. Journal of Mechanical Science and Technology, 24(10): 2017-2023. DOI:10.1007/s12206-010-0613-z
李磊, 张立杰, 张宁, 等, 2010. FLUENT在复杂地形风场精细模拟中的应用研究[J]. 高原气象, 29(3): 621-628.
肖仪清, 李朝, 欧进萍, 等, 2009. 复杂地形风能评估的CFD方法[J]. 华南理工大学学报 (自然科学版), 37(9): 30-35.
U.S. National Renewable Energy Laboratory, Technical Report-Central America Wind Energy Resource Assessment, http://swera.unep.net, 2006-Aug-21.
U.S. National Renewable Energy Laboratory, Technical Report-China Wind Energy Resource Assessment, http://swera.unep.net, 2006-Aug-21.
付斌. 针对PACA复杂区域的风资源绘图中尺度方法与微观尺度方法相结合. 2009第六届亚洲风能大会暨国际风能设备展览会, 2009年7月8-10日, 北京.
王鹏, 朱蓉, 方艳莹, 2011. 广东海陵岛风能资源高分辨率数值模拟研究[J]. 风能, No.5: 53-61.
杨罡, 刘树华, 朱蓉, 等, 2011. 鄱阳湖地区大气边界层特征的数值模拟[J]. 地球物理学报, 54(4): 896-908.
何晓凤, 周荣卫, 朱蓉, 2010. MM5与CFD软件相结合对复杂地形风资源模拟初探[J]. 资源科学, 32(4): 650-655.
法国美迪公司Meteodyn WT用户使用手册.
Delaunay D, Chantelot A, Guyader T, et al. Meteodyn WT: A software for wind resource assessment in complex terrain. European Wind Energy Conf., London (UK), November, 2004.
Yamada T, 1983. Simulations of nocturnal drainage flows by a q2l turbulence closure model[J]. J Atmos Sci, 40: 91-106. DOI:10.1175/1520-0469(1983)040<0091:SONDFB>2.0.CO;2
Arritt R W, 1987. The effect of water surface boundary layers[J]. Boundary-Layer Meteorology, 40: 101-125. DOI:10.1007/BF00140071
冯长青, 杜燕军, 包紫光, 等, 2010. 风能资源评估软件WAsP和WT的适用性[J]. 中国电力, 43(1): 61-65.
肖子牛, 朱蓉, 宋丽莉, 等, 2010. 中国风能资源评估2009[M]. 北京: 气象出版社.
Joshua Earnest, Tore Wizelius. Wind power plants and project development. PHI Learning Private Limited, 2011, New Kelhi:58. ISBN-978-81-203-3986-6.