复杂山地风力机流场及气动性能

邓宗伟^{1, 2}，彭文春³，高乾丰⁴，董辉⁴，朱志祥^{1, 2}

(1. 湖南城市学院 土木工程学院，湖南 益阳，413000；

2. 中南大学 土木工程学院，湖南　长沙，410075；

3. 中国水电顾问集团中南勘测设计研究院，湖南 长沙，410014；

4. 湘潭大学 土木工程与力学学院，湖南 湘潭，411105)

摘要：为研究复杂山地风力机流场的运动规律及相关气动性能，以具体的风速及地形图为依据，基于pro/E与Gambit软件建立XE93-2000型风力机整机模型和复杂山地流场模型，利用Fluent软件对复杂山区地形条件下的水平轴风力机外流场及气动性能进行分析。研究结果表明：在山顶上空60~200 m高度范围内风剪切较小，特别适合安装大型风力机，与主风向相切的山脊也是风力机的较佳安装场址；气流在绕过风力机后速度骤减且恢复过程较长，其后500 m之内仍受影响；塔筒受到的风荷载约占整个风电机组水平风荷载的15.6%；驱动力矩和阻力是风轮的主要气动荷载，均随风速的提高而显著增大。

关键词：山地风场；水平轴风力机；流场分析；气动性能；整机模型；CFD模拟

中图分类号：TK8          文献标志码：A         文章编号：1672-7207(2013)10-4294-07

Flow field and aerodynamic performance of wind turbine in complex terrain

DENG Zongwei^{1, 2}, PENG Wenchun³, GAO Qianfeng⁴, DONG Hui⁴, ZHU Zhixiang^{1, 2}

(1. School of Civil Engineering, Hunan City University, Yiyang 413000, China;

2. School of Civil Engineering, Central South University, Changsha 410075, China;

3. Mid-South Design and Research Institute,

China Hydropower Engineering Consulting Group Co., Changsha 410014, China;

4. College of Civil Engineering and Mechanics, Xiangtan University, Xiangtan 411105, China)

Abstract: In order to study the motion characteristics of flow field and related aerodynamic performance for complex mountain wind turbine, based on wind speed and topographic data, using pro/E and Gambit software, the models of XE93-2000 type wind turbine and complex mountain flow field were built. Using fluent software, the external flow field and aerodynamic performance of horizontal axis wind turbine were analyzed. The results show that due to small wind shear, the large-scale wind turbine is particularly suitable to be installed on the mountain top within the height range of 60-200 m, in addition, the ridge tangent to the main wind direction is also a good installation site for wind turbine. Airflow speed drops after the wind turbine and the recovery process is long, wind speed in 500 m scope after wind turbine is still affected. Wind load on the tower accounts for about 15.6% of the total wind load of wind turbine, driving moment and drag force are the main aerodynamic loads of wind wheel, and their values increase with the increase of wind speed significantly.

Key words: mountain wind field; horizontal axis wind turbine; flow field analysis; aerodynamic performance; whole model; CFD simulation

近年来，我国许多大中城市频繁出现空气质量问题，以煤炭为基础的火力发电弊端日益凸显，而水电开发利用程度接近饱和，因此，以风力发电为代表的新能源成为电力资源的有力补充。过去我国风资源的开发主要集中在三北地区及东南沿海一带，为了缩小能源分布的地区差距，促进中部及南部地区经济的可持续发展，目前，这些区域山区的风力发电也迎来了投资热潮。然而，山区地形复杂多样，风力机绕流流场的运动规律不易掌握，而风场的流动状况直接决定了风力机性能与经济效益，因此，弄清复杂山地风力机的绕流流场及气动性能对山区风电场的选址、风力机的研发和改进都具有重要意义。计算流体力学(computational fluid dynamics，CFD)是利用计算机和数值方法求解流体流动问题的一种计算技术。CFD软件模拟与实验研究相比可节约大量成本，因此，基于CFD理论的流场模拟及风力机气动性能分析越来越多地被采用，如许多学者采用各种CFD软件对物理地形风场^[1-3]、叶片气动性能^[4-6]及风力机各部分间相互作用^[7]等方面进行了大量研究，并取得了许多成果。一些研究者对真实山地风场、风力机整机绕流流场及气动性能进行了研究，如：邓院昌等^[8]以苏格兰Askervein小山为例，基于Fluent软件并结合NURBS地形建模方法，模拟了实际地形下的风场分布；李磊等^[9]尝试将Fluent软件用于复杂地形风场的精细模拟研究，确定Fluent完全可以用于小尺度复杂地形上风场的精细模拟；陈爱等^[10]用CFD软件分析了某油田井组简化三维地形的流场，得到了不同来流风向时的风场分布及山丘表面的风速分布图，并最终确定了风力机最佳安装位置及塔架高度；张果宇^[11]等采用Fluent对风力机整机的外流场进行了数值模拟，得到了风力机流场的速度分布、压力分布以及叶片受力情况。上述研究主要对原地形流场和风力机绕流流场进行单独模拟分析，但实际上，风力机的扰动绕流与原地形流场相互影响，两者不是独立的，且草原和海上风力机的研究成果亦不能对山区风电场的选址和风力机的设计提供正确的指导，因此，有针对性地对复杂山地和风力机进行整体建模分析就显得尤为重要。在此，本文作者以湖南郴州桥市风电场为背景，建立山地真实地面模型和大型水平轴风力机整机模型，运用Fluent计算流体软件对复杂山地风力机绕流流场及气动性能进行整体研究。模拟时，采用RNG k-ε紊流模型、移动参考坐标系模型MRF、考虑风剪切的速度输入以尽可能真实地模拟风力机绕流特性。经分析得到流场的速度分布和速度矢量以及风力机机舱、塔筒和风轮的受力、扭转力矩等气动载荷结果，以了解山地风力机绕流流场特有的运动规律，以便为山地水平轴风力机的微观选址和研发提供针对性指导，亦为进一步研究风力机的结构动力特性提供参考依据。

1  Fluent流场模拟

采用Fluent流体计算软件进行风场模拟，其过程包括前处理、模拟计算和后处理3个部分。首先由实际地形图和风机参数通过Gambit或其他前处理器建立几何模型，并用Gambit进行流场网格划分与边界设置，然后，将输出文件导入Fluent求解器进行求解设置和迭代计算，最后对计算结果进行后处理。图1所示为整个模拟过程的流程图。

图1  Fluent模拟流程图

Fig. 1  Flow chart of FLUENT simulation

2  实测风速特性

湖南郴州桥市风电场位于郴州北湖区鲁塘镇、桂阳县荷叶镇及临武县金江镇沿S214公路西侧一带，该区域地形突兀，山顶浑圆，属丘陵地带，海拔为400~800 m，属亚热带湿润季风气候，四季分明，气候温和，雨量充沛，光热充足，风源汇集，是湖南省风能资源最丰富的地区之一。该风电场一期工程于2010-11-23正式开工建设，是湖南省目前并网规模最大、单机扫风面积最大、建设时间最短的风电项目。一期工程总装机24台湘电风能生产的XE93-2000 型永磁直驱风力发电机组，采用叶片加长型风机，以降低风力机发电对风速的要求。目前，一期工程已并网发电，二期工程也正在施工建设中。本文的研究对象为该风电场二期工程的9号风力机。

通过对大量实地测得的风速风向进行处理，整理出9号风力机场址附近2009-10-12至2010-10-11的10 min平均风速风向。图2和图3所示分别为10 m高度处的风速、风向频率图，得到距地表10 m高度处的年平均风速为5.97 m/s，主风向为东南方向且比较稳定，说明该场地是风力机的合理场址。

采用指数模型，对不同高度处实测年平均风速加以整理并进行曲线拟合，得到如图4所示的线性函数，其表达式为

         (1)

式中：v_H为H高度处的风速；v₁₀为10 m高处风速。

图2  风速直方图及威布尔图

Fig. 2  Histogram and weibull diagram of wind speed

图3  风向玫瑰图

Fig. 3  The rose map of wind direction

α=0.200 8为风切变指数，对比规范指数取值表，该场地地表状态介于一般不平滑(短草地、庄稼地、乡村)与不平滑(树林、市郊)之间，这显然与实际情况相符。

图4  风切变指数拟合

Fig. 4  Fitting curve of wind shear exponent

3  复杂山地流场分析

3.1  地形建模与前处理

根据实测资料，选取山顶为原点以主风向为x轴正方向建立直角坐标系，将实际地形高程数据导入Gambit2.4.6生成风电场局部三维地形曲面，如图5所示，其水平投影范围为1 000 m×700 m。

图5  1 000 m×700 m三维地形

Fig. 5  1 000 m×700 m 3D terrain

为了尽量减少边界对地表附近流场的影响，流域高度取地表至y轴正方向以上500 m处。由于地表起伏不平，且考虑到计算机条件的限制，对上边界以下450 m范围采用较粗的结构化网格，地表附近则采用加密的非结构化网格，两部分粗细网格循序过渡，并通过多次尝试选取最佳的网格划分方式，最后得到737 180个网格体。

利用实测风速数据，设置入口边界为速度进口边界(velocity inlet)；求解前出口处流速和压力未知，设置出口为出流边界(outflow)；其他边界均设为壁面边界(wall)，等价于非光滑的无滑移边界条件；定义连续介质为FLUID。

3.2  模拟计算

利用Fluent6.3.26进行求解。选择压力基、定常、隐式求解器，操作环境考虑重力加速度g=9.8 m/s²。入口边界假设风速的垂直分布满足幂指数风廓线规律，风切变指数采用实测值α=0.200 8，通过用户自定义函数(user defined function，UDF)模块编程实现。湍流模型选用重整化群RNGk-ε模型，此模型由Yakhot和Orszag^[12]于1986年提出，它考虑了湍流中涡流因素的影响和低雷诺数效应，因而与标准k-ε相比具有更高的可信度和精度。其输运方程^[13]为

  (2)

式中：k和ε分别为湍动能和湍动耗散率；α_k和α_ε分别为k和ε的反向有效普朗特数；G_K为平均速度梯度引起的湍动能的产生项；和为经验常数；μ_eff为扩散系数；ρ为空气密度。

求解控制时，流动方程中压力-速度耦合采用SIMPLE算法。为了提高求解精度，离散格式均选择高精度的二阶迎风格式。在求解过程中，控制3个方向的速度分量、湍动能k、湍动耗散率ε的收敛。最后，初始化风场，进行迭代求解，当迭代到214次时残差收敛。

3.3  结果分析

山区地表起伏多变，近地面风场表现出很不均匀性。图6所示为距地表10 m 高度处的风速矢量图。从图6可以看出：(1) 进口距地表10 m高处风速为5.97 m/s，气流从进口到山脚一定高度因受到地表的摩擦作用，能量逐渐耗散，风速减小至进口速度的75%左右；(2) 气流从山脚攀升到山顶过程中速度逐渐增大，这是由于气流受到山坡的挤压而增速，即存在伯努利效应，山顶附近风速较大，最大为16.9 m/s；(3) 经过小山时，气流受到阻挡而向山峰两侧绕流，在与主风向相切的山峰两侧形成加速区；(4) 气流绕过山顶后速度回落到迎风面半山腰风速，为9~11 m/s，不过山背低洼处风速降低更明显。

图6  距地表10 m高处风速矢量图

Fig. 6  Wind vector at height of 10 m from ground

图7和图8所示分别为z=0剖面流场的风速云图和山顶上空风速沿高度分布的散点图。由图7和图8可知：(1) 速度入口边界以指数模型输入风速，迎风面山体高度以下范围内的气流速度等势线接近水平；(2) 由于受到山体绕流的影响，山顶上空60 m范围内气流速度急剧变化，风剪切偏离指数分布规律而表现出严重的“负切变”现象，而在60~200 m高度区域风剪切几乎为0，因此，山顶特别适合安装大型风力机；(3) 气流流经山丘的迎风面进入背风面时，气流受到扰动减速，但由于小山地表坡度较缓，在背风面没有形成明显涡流。

图7  z=0剖面流场风速云图

Fig. 7  Velocity contour at z=0 section

图8  山顶上空风速沿高度分布散点图

Fig. 8  Scatter plot of wind speed along height above the mountain top

可见：山顶和与主风向相切的山脊能使风速提高，是风力机的最佳安装场址；地表粗糙度较大的区域、低洼地带和背风坡会使风速降低，不宜安装风力机。

4  风机绕流流场分析

4.1  XE93-2000型风机建模

以湘电风能XE93-2000 型变桨距风力发电机进行研究，其总体参数见表1。由于风力发电机几何形状复杂，用Gambit直接建模很困难，故采用三维CAD软件Pro/E进行风力机整机建模，然后，通过Pro/E导出stp格式文件与Gambit联接进行网格划分等前处理工作。

表1  XE93-2000型风力机特性

Table1 Characteristics of XE93-2000 wind turbine

4.2  前处理设置

将风力机安装在计算域中心小山的山顶，且使风轮轴线与主风向对中。因风力机表面及地表形状复杂，对风力机各部分表面网格及附近的体网格、旋转域内的体网格及地表附近体网格采用非结构化网格，并进行加密；其他静止域采用相对较粗的结构化网格，粗细网格逐渐过渡，最终体网格总数达160万。图9所示为加密的风机表面非结构化网格。

风力机表面设定为无滑移壁面边界(wall)；由于计算域的尺度为风轮直径的5倍以上，因此，可忽略外边界对计算的影响，其边界条件亦按壁面边界来处理；旋转域与静止域的重合面设为interface边界；分别定义旋转域、静止域连续介质，其他设置与前面的相同。

4.3  模拟计算

求解器及湍流模型等的选择亦与3.2节保持一致。本模型中包含风轮旋转域和周围静止域，两者之间的藕合采用移动参考坐标系模型MRF，旋转速度取18 r/min。创建旋转域与静止域交界面，使Fluent自动对交界面两侧的压力和风速进行差值传递。仍按前面相同设置进行求解控制，然后，初始化风场，进行迭代求解。

图9  风机表面加密网格

Fig. 9  Encrypted meshes of wind turbine surface

4.4  结果分析

4.4.1  风机绕流流场

图10所示为风力机在停机状态和运行状态流场横剖面的速度云图。从图10可以看出：(1) 在停机状态下，塔影效应明显，但叶片之间没有相互影响，气流在叶片根部附近速度较大，最大达23.6 m/s；(2) 在运行状态下，风轮旋转引起扫掠面内的气流增速，高速区主要集中在靠近叶尖的1/2叶展附近，叶片尾缘气流产生分离；(3) 叶片周围的速度等值线密集，数值变化较大，说明叶片表面附近的压力梯度较大。

图10  x=0剖面流场速度云图

Fig. 10  Velocity contour at x=0 section

图11所示为风力机在停机状态和运行状态流场纵剖面的速度云图。由图11可知：(1) 风力机前后的流速变化较快，风轮在停机状态的阻风效应比运行时明显；(2) 风通过轮毂时在叶片间空隙绕流加速，形成了局部高速区；(3) 气流在绕过风力机后速度骤减且恢复过程漫长，风力机后500 m内仍受到影响，但是，高空处气流的速度比地表层气流的速度恢复快，此现象对确定塔筒高度和塔间间距具有参考价值。

4.4.2  风机气动载荷

(1) 轮毂及机舱表面压力。机舱尾外轮廓比较圆滑，其外表面受到的压力比较均匀；轮毂及发电机壳的前表面气流滞止处形成了相对高压区，而在侧表面形成了负压区，因此，在设计风力机各部件外形时应考虑平滑过渡以降低表面压力。图12所示为风力机轮毂与机舱表面的压力云图。

图11  z=0剖面流场速度云图

Fig. 11  Velocity contour at z=0 section

图12  轮毂与机舱表面压力云图

Fig. 12  Pressure in hub and nacelle for wind turbine

(2) 塔筒气动载荷。图13和图14所示分别为塔筒表面各点的压力随高度分布和顺风向分布的散点图。由图13和图14可知：1) 距塔底18 m高度范围内塔筒壁的压力明显低于此高度以上塔筒壁的压力，前者仅为后者的一半左右；2) 塔筒壁除了迎风面小范围内受到正压力外，绝大多数部位均承受负压力，且各截面最大负压力绝对值为最大正压力的1.2~2.0倍，塔筒壁最大正压力为214 Pa，最大负压力为-448 Pa；3) 计算得到塔筒水平风荷载合力为25.92 kN，约占整个风电机组顺风向水平风荷载的15.6%。

图13  塔筒壁压力沿高度分布

Fig. 13  Pressure distribution of tower wall with height

图14  塔筒壁压力沿顺风向分布

Fig. 14  Pressure distribution of tower wall along windward

(3) 风轮气动载荷。改变入口边界风速并考虑风剪切重复计算，得到不同来流风速条件下风轮的气动荷载情况，如图15所示。从图15可见：阻力和绕x轴的力矩(扫掠面内的扭转力矩)是风轮的主要气动荷载；当入口风速为5.97 m/s时，风轮阻力为140.65 kN，扫掠面扭转力矩为5 870.3 kN·m；随着风速的增大，风轮阻力和扫掠面扭转力矩亦逐渐增大，且风速越大其增加的幅度也越显著；而风轮升力和绕y轴和z轴的力矩均较小，随风速增大变化甚小，风轮整体升力为负值，说明风轮受到的升力向下。

图15  不同入口风速下风轮的受力情况

Fig. 15  Force of wind wheel under different inlet wind speeds

5  结论

(1) 气流经过山丘时受到挤压，山顶上空60 m内风速急剧变化，表现出严重的“负切变”现象，而在60~200 m高度区域风剪切几乎为0，因此，山顶特别适合安装大型风力机。与主风向相切的山脊能使风速提高，也是风力机的较佳场址；地表粗糙度较大的区域、低洼地带和背风坡使风速降低，不宜安装风力机。

(2) 风力机运行时叶片、塔筒之间相互影响，但此时风轮的阻风作用较停机时明显减小；气流在绕过风力机后速度骤减且恢复过程较长，风力机后500 m内仍受到影响，此现象对合理安排山地风力机间距具有参考价值。风轮旋转时引起扫掠面气流增速，叶尖旋转速度最大、表面形状突变，其尾流产生分离。

(3) 塔筒壁除了迎风面小范围内受到正压力外，绝大多数部位均承受负压力，且负压力的绝对值大于正压力。塔筒水平风荷载合力为25.92 kN，约占整个风电机组顺风向水平风荷载的15.6%。阻力和扫掠面内的扭转力矩是风轮受到的主要气动荷载，当入口风速为5.97 m/s时风轮阻力为140.65 kN，扫掠面扭转力矩为5 870.3 kN·m；随着风速的提高，风轮阻力和扫掠面扭转力矩显著增大。

参考文献：

[1] Yu F L, Akashi M. Numerical simulation of flow over topographic features by revised k-ε models[J]. Journal of Wind Engineering and Industrial Aerodynamics, 2003, 91: 231-245.

[2] 魏慧荣, 康顺. 风电场地形绕流的CFD结果确认研究[J]. 工程热物理学报, 2007, 28(4): 577-579.

WEI Huirong, KANG Shun. The validate and research of CFD result about the wind flow in wind farm terrain[J]. Journal of Engineering Thermophysics, 2007, 28(4): 577-579.

[3] 陈平. 地形对山地丘陵风场影响的数值研究[D]. 杭州: 浙江大学建筑工程学院, 2007: 21-24.

CHEN Ping. Numerical study of terrain influence on the airflow over hilly land[D]. Hangzhou: Zhejiang University. College of Civil Engineering and Architecture, 2007: 21-24.

[4] 任年鑫, 欧进萍. 大型风力机二维翼型气动性能数值模拟[J]. 太阳能学报, 2009, 30(8): 1087-1091.

REN Nianxin, OU Jinping. Numerical simulation for pneumatic characteristics for two-dimensional airfoils large wind turbine[J]. Acta Energiae Solaris Sinica, 2009, 30(8): 1087-1091.

[5] 范忠瑶, 康顺, 王建录. NREL Phase VI叶片气动性能数值分析[J]. 工程热物理学报, 2009, 30(10): 1665-1668.

FAN Zhongyao, KANG Shun, WANG Jianlu. Numerical investigation of aerodynamic performance of NREL phase VI wind turbine blades[J]. Journal of Engineering Thermophysics, 2009, 30(10): 1665-1668.

[6] 王琦峰, 陈伯君, 安亦然. 大型风力机三维空气动力学数值模拟[J]. 空气动力学学报, 2011, 29(6): 810-814.

WANG Qifeng, CHEN Bojun, AN Yiran. 3D numerical simulation for aerodynamic performance of large scale wind turbine[J]. Acta Aerodynamica Sinica, 2011, 29(6): 810-814.

[7] Khaled A, Christian M, Peter J E. 2D and 3D numerical simulation atmospheric boundary layer of the wind-rotor/nacelle interaction in an atmospheric boundary layer[J]. Journal of Wind Engineering and Industrial Aerodynamics, 2011, 99: 833-844.

[8] 邓院昌, 刘沙, 余志, 等. 实际地形风场CFD模拟中粗糙度的影响分析[J]. 太阳能学报, 2010, 31(12): 1644-1648.

DENG Yuanchang, LIU Sha, YU Zhi. Effect of roughness on CFD wind field simulation over natural terrain[J]. Acta Energiae Solaris Sinica, 2010, 31(12): 1644-1648.

[9] 李磊, 张立杰, 张宁, 等. FLUENT 在复杂地形风场精细模拟中的应用研究[J]. 高原气象, 2010, 29(3): 621-628.

LI Lei, ZHANG Lijie, ZHANG Ning. Application of FLUENT on the fine-scale simulation of the wind field over complex terrain[J]. Plateau Meteorology, 2010, 29(3): 621-628.

[10] 陈爱, 刘宏昭, 杨迎超, 等. 复杂地形条件下风力机微观选址[J]. 太阳能学报, 2012, 33(5): 782-788.

CHEN Ai, LIU Hongzhao, YANG Yingchao. Micro-siting technique for wind turbines under complex terrain[J]. Acta Energiae Solaris Sinica, 2012, 33(5): 782-788.

[11] 张果宇, 蒋劲, 刘长陆. 风力发电机整机气动性数值模拟计算与仿真研究[J]. 华东电力, 2009, 37(3): 449-452.

ZHANG Guoyu, JIANG Jin, LIU Changlu. Numerical simulation of aerodynamic performance for wind turbines[J]. East China Electric Power, 2009, 37(3): 449-452.

[12] Yakhot V, Orszag S A. Renormalization-group analysis of turbulence[J]. Journal of Scientific Computing, 1986, 1(1): 3-51.

[13] 王福军. 计算流体动力学分析: CFD软件原理与应用[M]. 北京: 清华大学出版社, 2004: 1-20.

WANG Fujun. Computational fluid dynamics analysis: The principle and application of CFD software[M]. Beijing: Tsinghua University Press, 2004: 1-20.

(编辑  陈灿华)

收稿日期：2013-01-12；修回日期：2013-03-05

基金项目：湖南省科技计划项目(2013GK3086)；中国水电顾问集团中南勘测设计研究院科技项目(YJ2012⒉6)

通信作者：邓宗伟(1972-)，男，湖南安化人，博士，副教授，从事岩土与基础工程研究；电话：13973760738；E-mail：teapotd@163.com