基于OpenFOAM的层结性对风力机尾流大气参数的影响

吴正人，张智博，刘 梅，杜思源，张 雷

（1.华北电力大学 动力工程系，河北 保定 071003；2.华北电力大学经济管理系，河北 保定 071003）

0 引言

随着全球经济的快速发展，优化能源结构已成为可持续发展的重要组成部分。风能为优化全球整体能源结构和满足世界整体能源供应做出了巨大贡献。风力机作为开发风能的发电装置，其运行必定会破坏大气层的物质能量循环，因此研究风力机运行对大气参数的影响具有重要的意义。袁仁育[1]采用动态入流条件对风力机运行进行数值模拟，开发多尺度下的风力机模型，发现风电场会对风力机尾流10 km内的大气边界层产生影响。刚蕾[2]研究了大气稳定性对风力机尾流的影响，其提出的Park-Gauss模型对尾流区风速模拟取得了较好的效果。Xie S[3]通过大涡模拟分析不同温度层结下的风力机尾流，发现尾流的风速、湿度和温度变化受温度层结影响较大。Sharma V[4]利用大涡模拟研究发现，区域大气层在风电场影响下升高了200 m。徐荣会[5]通过研究苏尼特右旗某一风电场区域的气候变化情况，发现风速和陆地地表温度会随着风向变化逐渐下降，风力机尾流区域的近地层会产生降温增湿的现象。通过对苏格兰Black Law风电场的气象数据进行分析，Armstrong A[6]发现风力机周围环境温度升高了0.18℃、湿度升高了0.03 g/m3。Lu H[7]通过大涡模拟与风力机致动线技术结合的方法，发现科里奥利效应和叶片运动引起的旋转效应会导致边界层高度增加，其表面动量通量减小了30%，表面浮力通量减小了15%。Emanuel K[8]通过全球气候模型设想在全球各地部署足够的风力机来开发整个风能资源的极端情况，发现平均表面温度可能变化几摄氏度，两极可能降低多达10℃。张亚光[9]通过模拟3种不同的大气稳定度，发现随着稳定度的增加，尾流湍流强度和剪切应力下降。杨祥生[10]通过新修正的Park-Gauss和Park-Polynomial模型来探究大气边界层稳定性对风电场微观选址的影响，发现大气边界层越不稳定，越能促进风力机尾流的恢复。

本文采用开源CFD平台OpenFoam，实现对PimpleDyMFoam求解器的二次开发，并结合大气能量守恒方程编译求解温度场的分布，采用codedFixedValue边界条件建立大气边界层的风廓线模型和大气温度层结模型，从而建立不同温度层结下的风力机模型，并在此基础上研究风力机运行对大气边界层参数的影响。

1 模型建立及数值模拟方法

1.1 控制方程组

不可压缩流动大涡模拟的连续方程和N-S方程可表示为

式中：υθ为热扩散系数；Rj为xj方向的净辐射分量；cp为湿空气的定压比热；E为单位时间相变的水量。

方程左边第一项表示能量的贮存，第二项表示平流作用。方程右边第一项表示分子热传导作用，第二项表示辐射散度，描述xj方向净辐射通量密度的梯度引起的空气层温度的变化，第三项表示因相变水汽凝结或蒸发时所产生的能量变化。

本文将大气作为不可压缩均质流体处理，同时不考虑气压的变化对大气温度的影响，所以采用温度T代替位温θ，根据大气能量守恒方程简化的温度方程为

1.2 几何模型及网格划分

以1.2 MW水平轴风力机为研究对象，应用坐标变换理论得到叶片截面的坐标数据，从而建立风力机叶片实体模型，进而完成风力机流场的建模工作。风力机参数如表1所示[12]。

表1 风轮及风况参数Table 1 Wind wheel and wind condition parameters

风力机运行时承载负荷的主要区域在轮毂和桨叶连接处，所以选用NACA634翼型，而叶尖为了提升叶片的气动特性选用FX66S196翼型。风力机d为60 m，桨叶安装角和风轮锥角均为0°，轮毂直径为2.8 m，转轴倾角为4°。通过Gambit软件建立的叶片实体模型如图1所示。

图1 叶片实体模型图Fig.1 Solid model diagram of blade

流场的计算域尺寸设置如图2所示。

图2 整体流场区域尺寸Fig.2 Overall flow field area size

流向沿Z轴，风轮中心高度为60 m。为了使入口风得到充分发展，把风力机设置在流场入口后3d的位置处。叶片区域网格划分如图3所示。

图3 叶片区域网格划分Fig.3 Grid division of blade area

通过比较模型中心线上湍流动能变化曲线来验证网格无关性（图4）。对比网格数分别为595万、652万和724万的模型，发现数目较大的两个模型的计算结果接近，而数目较少的网格模型计算结果有所偏差，所以选用652万的网格模型进行计算。

图4 网格无关性验证Fig.4 Grid independence verification

1.3 模拟条件设置

为了更好地适应模拟的实际工况，本文采用codedFixedValue边界条件对入口风速的垂直分布以及入口的温度进行编译。入口压力设置为zeroGradient边界条件；对于流场区域的出口，压力采用toalPressure边界条件，速度采用pressureInletOutletVelocity边界条件；在风力机所处的旋流区域，因为与流场的其他区域进行物质交换，采用cyclicAMI边界条件；对于流场区域的顶部和两侧采用symmetry边界条件，而底部需采用无滑移边界条件。

风力机的来流风速由风廓线方程定义。

式中：u（zr）为参考风速，u（zr）=10.0；zr为参考高度，zr=10.0；α为指数幂，α=0.12。

入口风温为305 K，为了加快收敛速度防止发散，采用adjustableRunTime，即变时间步长的控制方法，时间步长设为10-4s，所有参数的残差设置为10-4。在LES下的S-A湍流模型基础上，对时间项采用欧拉离散，梯度项、散度项和拉普拉斯项均采用高斯方法，线性插值为二阶离散。风力机周围小区域内流体随着叶片及轮毂共同以19.27 r/min的转速匀速旋转，其余外部区域的流场静止。

1.4 模拟方法验证

将风力机尾流z=7d以及z=10d处的速度剖面和瑞典航空研究院（FFA）风洞实验数据结果[13]进行对比（图5），该风洞采用L1回流式风洞，长、宽、高分别为6，2.4 m和1.4 m，风洞自由来流的湍流水平为0.5%，风速为8 m/s；风力机模型直径D为0.25 m，轮毂高度为1D。

由图5可知：数值模拟的数据与风洞实验的数据在分布规律上具有一定的相似性，风力机的近尾流区速度亏损明显；风洞实验数据较数值模拟数据偏小。

图5 模拟结果与风洞实验结果尾流速度剖面比较Fig.5 Comparison of wake velocity profile between simulation results and wind tunnel test results

2 计算结果与分析

2.1 速度场的对比分析

对风力机尾流影响最大的因素是当地的天气条件和风电场入口自由风的状态，自由风被风力机吸入后动能减小，使得尾流区原本正常的湍流输运和交换循环发生改变。风力机轮毂高度处的速度分布如图6所示。

图6 轮毂高度处的速度分布Fig.6 Velocity distribution at hub height

由图6可知，两种温度层结状态下的速度分布规律基本一致，但稳定层结下尾流对区域大气边界层的影响范围明显比不稳定层结远。这是因为在不稳定层结下，高处空气的温度低、密度大，从而有空气向下流动及流速增大的现象。另外尾流会与不稳定边界层中的大尺度旋涡产生剧烈的相互作用，使得尾流剪切层的形状被快速破坏，且空间分布更加均匀，速度的补充也较快。

在不稳定层结下，尾流在3d后呈现蜿蜒特性，而在稳定层结下，尾流的形状相对固定，这反映了稳定层结对剪切速度的抑制作用。

图7为尾流区不同距离处风速的垂直分布曲线。图中向下凸起的范围即为尾流的高剪切层，也是尾涡结构的主要部分。

图7 尾流区不同距离风速的垂直分布Fig.7 Vertical distribution of wind speed at different distances in the wake region

由图7可知：在稳定层结下，大气的速度损失最大为7 m/s，而不稳定层结下为5 m/s；稳定层结和不稳定层结下风速的变化规律基本相似。不稳定层结下的速度损失较小是由于自由风的风速较大，动能交换剧烈。观察不稳定层结z=1d处可以发现，叶片尖端处的风速比轮毂处要小，这正好描述了尾流中高剪切层是呈圆柱形分布的，而稳定层结中该现象并不明显，体现了稳定层结的自由风分布较稳定，并对剪切速度有明显的抑制作用。其它曲线由于逐渐远离风力机，损失的动能逐渐从周围空气中得到补充，所以不断贴近入口风速的分布曲线。

2.2 温度场的对比分析

温度作为研究大气边界层的关键参数，其影响程度虽不及风速，但仍值得考虑。对比风廓线方程和温廓线方程可以看出，风廓线的变化率要远大于温度，使得局部区域的温度分布差异很小。同时伴随着尾涡结构的破裂，对流换热的作用在远离风力机的位置也逐渐减弱，但并不影响从中探讨温度参数的变化。图8为轮毂高度处的温度分布图。

图8 轮毂高度处的温度分布Fig.8 Temperature distribution at hub height

由图8可知：在远尾流区，不稳定层结的温度较来流温度的变化比稳定层结明显；在不稳定层结下，远尾流区域一侧温度较来流温度下降约为0.7 K，这是因为由于风轮的转动，高处的低温空气与低处的高温空气进行混合，出现了风力机一侧温度较低的现象；而稳定层结下的温度变化规律与不稳定层结相反，其温度升高约为0.15 K；在远尾流区域，风轮轮毂高度处的温度较来流温度的变化程度均会随着距离的增大而增大。

图9为尾流区不同距离温度的垂直分布曲线。

图9 尾流区不同距离温度的垂直分布Fig.9 Vertical distribution of temperature at different distances in the wake region

由图9可知：在高度为1.8d附近，温度有明显的界限点，高于此值的尾流区温度变大，而且作用的范围超过1d以上，并且随着尾流的发展，其温度升高的幅度减小；在1.8d高度以下的尾流区，温度较低，这说明风力机的运行会使其处于降温的环境中；在高度为0.5d以下时，尾流区各温度曲线的变化相似，说明在此区域的温度分布不受尾流发展的影响，在高度为3d以上时也有类似规律。

在稳定层结下，高度为0.5d以下区域的温度升高，而高度为0.5d以上区域的温度则下降，并且温度升高的幅度小于温度降低的幅度，这反映了风轮的旋转会使高温空气向下扩散，并与下垫面的低温空气混合，使得下垫面出现温度升高的现象。在尾流1d处降幅达到0.15 K，5d处降幅达到0.075 K，在高度为3.5d以上的区域温度已不发生明显变化。在两种层结状态下，随着尾流的发展，曲线弯曲程度变小，说明尾涡结构的不断分解使其作用的范围也在变小。

上述分析印证了文献[6]中得到的风电场下游区域的温度比上游区域的温度明显升高的现象，并由此得出夜晚风电场会显著影响尾流区下垫面空气温度的结论。

3 结论

本文采用OpenFoam对1.2 MW风力机进行了模拟分析，得到以下结论。

①通过开发编译特定环境下的OpenFoam求解器，以实现不同层结状态下的风力机模型，与固定求解器相比，其模拟方式有较好的灵活性和针对性。

②在稳定层结下，风力机尾流对大气边界层的影响范围比不稳定层结远，不稳定层结下的风速损失比稳定层结小。但两种层结下风力机尾流的风速变化趋势基本相似，随着尾流距离的增加其变化的趋势变小，且速度逐渐向入口风速恢复。

③在不稳定层结下，远尾流区温度较来流温度的变化比稳定层结明显。在不稳定层结下，轮毂高度处远尾流区的温度会下降，而稳定层结下的温度则会上升，并且这种变化趋势会随着距离的增加不断加大。

④随着高度的增加，不稳定层结下的风力机尾流温度呈先下降再上升后下降趋势，而稳定层结下的尾流温度变化趋势与其相反；在不稳定层结下，风轮高度以下区域的尾流温度分布不受尾流发展特性变化的影响，且高度为3d以上的区域有相似的规律。