基于优化CFD模型的风电场风速数值模拟
马致远;孙凯丰;侯金锁;于殿富;汪濙海;高乐
【摘 要】为了提高风资源普查的精度,更好地针对我国地形及风况,文章优化了现有风资源计算流体力学模型,并编写了相应计算模块.优化模型包括:①贴合复杂山地地形的网格化分器,可以对任意地形进行网格划分;②通过分析测风数据自动计算湍流模型系数;③增加温度运输方程,将大气边界层热稳定度耦合到动量方程和湍流模型中;④与实际大气边界层热稳定度分层效应一致的入口条件及壁面函数.为了验证优化后的风资源计算方法的精度,文章对一待开发风电场进行了风资源计算.计算结果显示,使用优化的模块可以更精确地计算风速,与优化前相比,可以将误差至少降低10%.%To improve the wind resource assessment accuracy,especially for China complex terrain,this paper optimized the current Computational fluid dynamic (CFD)wind resource assessment models,and programmed corresponding computing modules. The optimized functions include:①Mesh tool for complex terrain to make the mesh more attached to the real terrain surface,can mesh any terrain; ②Automatic calculation of turbulent model parameters based on the analysis of the wind mast tower data; ③ Adds temperature transport equation,coupled the atmospheric boundary layer thermal stratification effects into momentum equation as well as turbulent models; ④Consistent inlet boundary condition and wall function considering atmospheric boundary layer thermal stratification. To evaluate the new wind resource assessment method,this paper conducted an under-developing wind farm resource assessment. The results show the optimized modules can reduce the error by more than 10% compared with the one without these modules.honeys中国
【期刊名称】《可再生能源》
关典史【年(卷),期】2018(036)003
【总页数】7页(P446-452)
【关键词】风电场开发;风资源数值计算;大气边界层稳定度;湍流模型
【作 者】马致远;孙凯丰;侯金锁;于殿富;汪濙海;高乐
【作者单位】中电投电力工程有限公司 新能源分公司,上海220023;中电投电力工程有限公司 新能源分公司,上海220023;中电投电力工程有限公司 新能源分公司,上海220023;中电投电力工程有限公司 新能源分公司,上海220023;中电投电力工程有限公司 新能源分公司,上海220023;中电投电力工程有限公司 新能源分公司,上海220023
【正文语种】中 文
【中图分类】TK81
0 前言
随着计算机性能的提高，高精度的流体力学计算（CFD）已经被广泛的用于工业设备的设计制造[1]，[2]。近年来，风电场的前期风资源普查也越来越依赖于CFD计算 [3]～[7]。目前工业界主流商用CFD 软件如 Fluent，Star-CCM+，NUMECA 等均可用于风场相关研究[7]～[9]。然而，因为专业性和易操作性，目前大部分风资源计算软件市场被WindSim和 Meteodyn WT占领，文献[10]～[12]比较了这几款软件的优缺点以及在我国风电场微观选址中的应用。
开平方机器
孙志刚事件目前，风资源计算软件的核心流程是确定风场范围后，计算不同风向条件下的流场（一般按圆周平均划分为16份），得出风加速因子ai=（其中，Vfan，CFD，i为风向 i时风机点位的CFD 计算风速，Vmast，CFD，i为测风塔点位的 CFD 计算风速）；之后带入测风数据风向和风速时间序列（时间间隔一般为10 min，整个数据满足一完整 年 [13]）， 加 权 得 出 风 机 点 位 风 速 Vfan，j=（j为时间序列中数据点，i为不同测风塔，fi为不同测风塔对应的权重）；再结合风力发电机组功率曲线（Cp-V曲线）得出发电量。
大量的工程应用和研究表明，这些软件并不完全适用于我国目前待开发风电场。目前待开发地区的地形普遍较为复杂，大气边界层热稳定度对风切变有极大的影响[14]～[16]。此外，随着智能运维系统的日益更新，对风参数预报的需求也随之提高，商用软件底层均为黑箱，难以与现有实时天气预报（WRF）系统耦合[17]。各大风电机组制造商和高校均开始了自主软件的开发工作[18]～[20]。
本文开发了一款基于OpenFOAM 4.1版本的风资源计算软件，可并行于Ubuntu16.04操作系统的工作站上。该软件包括了测风数据处理、网格化分、后处理和不可压稳态雷诺平均方程求解器。其中，前三者基于python和C++，由作者独立编写，后者基于OpenFOAM开
x射线机源代码使用C++进行了改编。该计算软件通过优化网格、湍流模型参数和大气稳定度模型，提高了计算精度。
1 方程
为了考虑大气边界层稳定度的影响，本文在不可压稳态雷诺平均方程的基础上加入了使用Boussinesq近似计算的浮力项：
式中：Ui为速度；P为压力；T为温度；ν为空气粘度系数；β为热膨胀系数；g为重力加速度；θ为温度势场。
湍流粘度系数 νt=Cμκ2/ε，κ，ε 分别为湍流动能及湍流动能扩散率，本文使用κ-ε运输方程求解。
在风场计算中，Cμ的值取为 0.01～0.09[3]，Cμ也可以通过分析测风数据得到估算值[6]，[21]，[22]。 在大气边界层流动区域，假设湍流动能生成与耗散达到平衡（ε=Pk/ρ），并结合 Prandtl混合长度理论可得。因为表面层中剪应力变化较小，可认为。为了确定摩擦速度u*，选取测风塔处接近中性大气的数据进行拟合，选择条件为
吹填工程
①3个测风高度（40，60，70 m）平均风速大于3 m/s的数据；
②湍流强度比较接近中性大气条件下经验公式所得值的0.5～2.0倍的数据，经验公式为Iu=其中，z选为70 m，z0为测风塔位置地表粗糙度，约为0.051 7。
本文计算的风电场地区测风数据为逐10 min记录，数据较完整，使用多年测风数据进行插补得到完整年数据。其中，一座测风塔风瑰图如图1所示。图1中百分比为相应圆周线处所占比例，如该地区主风向为西北方向（315°，约占全年测风数据的18.41%），平均风速约为7.6 m/s，其余时间点上风向主要分布在西南至北。
图1 测风数据风瑰图Fig.1 Wind data rose plot
图2 湍流强度在一天中的变化Fig.2 Turbulence intensity in one day
图2为1 d湍流强度的变化图，从图2可知，在 17：00-翌日 7：00，因地面温度较低，大气边界层分层处于稳定状态，因此，湍流强度较小。而在7：00-17：00，因地面温度较高，大气边界层底部气流受热力影响与上部气流交互混合，湍流强度较大。因此，在处理测风数据时选取湍流强度较小的时段，并假设测风高度均在大气边界层的表面层以内，湍流动
能取3个测风高度上的平均值。通过以上公式可得Cμ约为0.083。同时根据边界层理论，κ 是 Von Karman 常数，代入Cμ即可求得C2ε约为1.97。本文也比较了其他文献中使用的κ-ε模型参数。
大气边界层稳定度对风速有极大的影响，例如当大气边界层处于不稳定态时，底部气流受热力影响与上部气流加剧混合，从而降低上部气流风速。在通常的风资源计算中，默认大气边界层为中性稳定度，但在实际情况中，每个季节或者每一天中，边界层稳定度都是不断改变的。本文通过在入口风廓线和壁面函数中使用Monin-Obukov长度来体现其对计算的影响 [21]～[24]：
式中：z为离地面高度；z0为地表粗糙度；L为Monin-Obukov长度。
在没有温度数据和中尺度计算结果的情况下，本文通过拟合测风数据近地曲线来估计大气边界层稳定度特征长度[19]，[20]。同时考虑到平均划分圆周做风电场的定向计算不够精确，因此，在计算前根据测风数据中的风向频率在主风向上进行加密划分。
2 计算域及求解器设定
本文开发了根据地形自动化分网格的网格编辑器。该模块可自动读入STL格式的地形文件，并根据提取点（例如风机和测风塔点位）坐标或用户自定义数值来定义计算域水平大小，计算域的高则根据 h=zmin+max[5（zmax-zmin），2 800]计算，其中zmin和zmax分别为最低和最高高度。为了使得网格更加的贴合实际地形，采用非结构六面体网格划分，并在计算核心区域进行了加密，第一层网格底面法向垂直于地面。整个计算域为长方体，根据划分的风向旋转计算域，选择相应的入口、出口和侧边界。其中，入口为根据式（6）定义的速度入口边界，出口为压力出口边界（p=0，∂u/∂n=0，n 为该面法线方向），两侧使用零梯度边界（∂u/∂n=0 ），地面设置为无滑移边界（u=0），天空设置为零梯度边界（∂u/∂n=0）。网格化分及边界如图3所示，其中风向角度以正北风（0°）为起点顺时针计算。在该网格上使用有限体积法离散求解公式 （1）～（5），离散的矩阵通过多层网格嵌套提高求解速度。
图3 入流角度为11.25°时网格划分及边界示意图Fig.3 Mesh and boundary at the wind direction of 11.25°
本文开发的软件底层为C++和Python编写，方程求解模块基于作者改编的OpenFOAM 4.1
开源程序，后处理主要为风资源计算相关模块，包括计算风电场风速图谱、风机点位风速、湍流强度、入流角、风切变和发电量，可多核并行于Ubuntu 16.04操作系统的工作站上。
本文的风电场数值模拟是为了验证开发的软件在风速计算上的准确性，因此，只进行了测风塔风速互推。选取了某待开发风电场作为计算域，该待开发地区有两座相距5.5 km的测风塔，可使用测风数据对本文开发的软件计算精度进行验证。计算中，使用1∶2 000地形图及300 m分辨率粗糙度地图。核心计算域为8 km×8 km，网格为均匀划分，分辨率为50 m。在核心计算域的基础上外扩10 km，网格扩展率为1.05，总计算域为28 km×28 km。在海拔方向，计算域高5 400 m，第一层网格高度为4 m，垂直扩展率为1.05，总网格数量约为911万。
