5 MW风力机塔影效应双向流固耦合分析

李 丽， 潘 扬， 胡丹梅

(1.国网湖南省电力有限公司防灾减灾中心(电网输变电设备防灾减灾国家重点实验室)，长沙 410129；2.上海电力学院 能源与机械工程学院，上海 200090)

近年来，由于风力机单机容量不断增大，风力机叶轮的扫风面积和工作高度也随之提高，导致作为支撑结构的塔筒规格不断增大，塔筒引起的塔影效应也更明显。塔影效应是指塔架对流场产生干扰，而流场将此干扰传递到风力机叶片上，导致风力机的气动特性发生变化。由于风力机塔架对空气有阻挡和排挤作用，使得其附近的流场发生变化，空气作用于叶片上的气动载荷发生剧烈波动，造成叶片的形变与其应力也会发生相应变化。

对于目前的兆瓦级大型风力机而言，塔影效应已经成为一个不可忽视的影响因素。研究表明，塔影效应使单个叶片载荷波动20%～40%,风轮载荷波动6%～12%[1]，说明塔架会对风力机组的运行产生影响，因此塔影效应具有一定的研究价值。

1968年，研究人员就提出了风力机塔影效应的概念[2]。进入21世纪后，随着风电产业的兴起，风电从业者开始认识到塔架对风力机的影响不可忽略。Das等[3]建立了关于塔影效应和分切效应的时域模型，用于计算塔架造成的风力机周期波动。Sescu等[4]利用商业CFD软件对塔影效应进行了仿真计算和分析。李少林等[5]在设计动态风力机模拟器时对塔影效应进行了模拟。

流固耦合是指流体与固体之间发生相互干涉，如果这种干涉既包括流体对固体的干涉，又包括固体对流体的干涉，则称为双向流固耦合。国内外学者针对风力机双向流固耦合情况下的工作状态进行了研究。Hsu等[6]和Bazilevs等[7-8]利用2种方法对NREL 5 MW海上风力机[9]进行了单根叶片和整机的双向流固耦合数值计算，研究内容包括风力机固体结构和流场。潘旭[10]利用Ansys软件对兆瓦级风力机叶片进行了双向流固耦合数值模拟，分析了蒙皮材料与叶厚等因素对叶片工作表现的影响。何福添[11]对3 MW风力机叶片进行了双向流固耦合数值模拟，并根据模拟结果进行了叶片设计优化。李媛[12]对二维翼型和三维叶片进行了流固耦合，并研究了不同风速条件下风力机叶片的流固耦合特性。目前，研究人员在研究风力机双向流固耦合时往往仅关注单根叶片或单独叶轮，对风力机整机塔影效应的研究相对较少。

笔者以NREL 5 MW风力机[9]为模型，通过数据交互模块System Coupling将流体模拟软件Fluent与固体分析软件Transient Structural进行联动，实现风力机固体结构与其周围流场的双向耦合，展开了关于塔影效应的流固耦合研究。

1 计算模型与网格划分

1.1 物理模型与网格

风力机模型叶轮直径为126 m，额定风速为11.4 m/s，额定转速为12.1 r/min，将转速取整为12 r/min，叶轮旋转周期为5 s。塔筒为圆台状，上底面直径为3.87 m，下底面直径为6 m，高为87.6 m。首先利用3D建模软件ProE建立风力机整机模型和立方体流场域模型，再根据布尔运算得到风力机的计算流场。流场分为旋转域和静止域，旋转域是半径为70 m、高为8 m的圆柱体，叶轮位于旋转域中央，被其完全包裹；静止域为长方体，截面高为400 m，宽为600 m，入口距叶轮旋转平面为125 m，出口距叶轮旋转平面为500 m，机舱与塔筒位于静止域内，流场域整体模型如图1所示。通过二次开发端口将模型直接导入Ansys Workbench，对风力机叶轮进行厚度赋值，使其成为厚度为8 cm的壳结构，将机舱与塔筒视为刚体。

1.1.2 划分计算网格

利用流固耦合的Fluent流体模块(FFF)自带的Mesh功能对流场进行网格划分，为能启用Fluent中动网格的弹簧光顺功能和网格重构功能，采用四面体网格，并对固体表面进行加密，网格单元个数为259万，最大网格偏斜小于0.85，满足计算要求。利用Workbench中的Mechanical Model对风力机叶轮进行壳单元网格划分，生成1.6万个六面体网格，如图2所示。叶片采用45°铺层各向异性玻璃纤维环氧树脂复合材料，其材料力学性能如表1所示，其中E1、E2分别为垂直纤维方向与沿纤维方向的弹性模量，G12为剪切模量，ρ为密度，σ12为泊松比[13]。

图2 叶轮壳网格

E1/GPaE2/GPaG12/GPaσ12ρ/(g·cm-3)398.63.80.282.1

1.1.3 可信度验证

为验证计算的准确性，分别利用627万、190万结构化网格和415万、259万非结构化网格在非流固耦合情况下进行计算，所得功率偏差分别为0.28%、3.6%、3.2%和5.3%。由于流固耦合计算非常占用计算机资源，为节省计算时间和提高计算机工作稳定性，最终采用259万非结构化四面体网格，计算误差约为5%，满足计算要求。

1.2 计算方法与耦合方程

由于NREL 5 MW风力机的几何尺寸较大，为了使网格数量保持在计算机可计算范围内，需对网格的最小尺寸进行控制，这导致固体表面的网格尺寸难以模拟出边界层。因此，湍流模型选用带有湍流漩涡修正和低雷诺数修正的RNGk-ε两方程模型，以达到利用标准壁面函数[14]处理边界层的目的。研究表明，k-ε模型对于流动的预测具有一定准确性[15-16]。采用滑移网格法进行流体计算，设定旋转域转速为12 r/min，时间步长为0.1 s，流场入口风速为11.4 m/s，出口为自由出口，离散格式为二阶迎风，利用Simple算法进行求解。对于叶轮固体部分，通过添加约束限制轮毂自由度，只保留旋转方向自由度，并在叶轮上加载角速度与标准重力，将整个叶轮外表面设为流固耦合面。

Ansys-Fluent的流固耦合方程如下[11]:

(1)

式中：A为系数矩阵，其中下标F代表流体域，S代表固体域，I代表流固交界面；U为流场中的速度；p为流场中的压力；δ为叶片位移；R为残差。

利用SPSS软件进行描述性统计分析；利用Aquachem水化学软件绘制Piper三线图进行水化学类型分析；选取超标较严重的硫酸盐(SO2-4)、亚硝酸盐(NO-2)、硝酸盐(NO-3)、总硬度和溶解性总固体(TDS)共5项评价指标，利用熵权密切值法对地下水水质进行评价。评价的标准采用《地下水质量标准》(GB/T14848—93)将地下水质量划分为5类(Ⅰ—Ⅴ)，见表1。

如图3所示，利用Workbench中的System Coupling连接Fluent模块和Transient Structural，并在System Coupling中设置Fluent的计算顺序为1，Transient Structural的计算顺序为2，使流体先作用于固体，固体再反作用于流体。将固体耦合面和流体耦合面配对，设置时间步长为0.1 s，进行耦合计算。

图3 Workbench模块结构图

2 计算结果分析

为研究塔影效应对叶轮的影响，笔者在相同的工况条件下对有塔筒和无塔筒这2种情况进行计算，并进行对比分析。

2.1 风力机叶轮气动载荷

图4和图5分别为风力机叶轮在25～30 s内转矩和所受轴向推力的变化曲线。在25～30 s内叶轮旋转了1周，转矩和轴向推力均出现3次峰值和3次谷值，验证了风力机的3P闪变[17]。由图4和图5可知，在有塔筒的情况下2种曲线的波动程度均比无塔筒时更剧烈。在风力机工况基本稳定的时间段内(25～30 s)，在有塔筒的情况下转矩峰值约为4 030 kN·m，谷值约为3 800 kN·m，峰谷值之差为230 kN·m；无塔筒时转矩输出的峰值与谷值分别为4 000 kN·m和3 930 kN·m，峰谷值之差为70 kN·m。对于叶轮所受轴向推力，在有塔筒的情况下轴向推力的峰值和谷值分别约为787 kN和743 kN，最大与最小值之差为44 kN；在无塔筒的情况下轴向推力的峰值和谷值分别约为779 kN和759 kN，最大与最小值之差为20 kN。综上，塔筒对风力机叶轮所受气动载荷有明显影响。以本计算为例，在塔影效应影响下叶轮所受转矩和轴向推力的波动幅度是未受塔影效应影响时的3倍以上，说明对于兆瓦级大型风力机，塔影效应的影响不可忽视。

图4 叶轮转矩曲线

图5 叶轮轴向推力曲线

2.2 叶轮形变与应力

图6给出了模拟时间t为10 s、20 s和30 s时风力机叶轮的形变放大图，放大倍数约为19倍。由于叶轮旋转周期T为5 s，在时间为5 s的整数倍时，z方向的叶片位于塔筒正前方，如图6(b)中叶片1所示。由图6可以看出，在10 s时风力机叶轮工况尚未稳定，叶片的形变较大，位于塔筒前方叶片的轴向挥舞形变和周向摆振形变均小于其他2根叶片；30 s时叶轮工况基本稳定，3根叶片的摆振形变程度较为接近，但挥舞形变仍有一定差异，此现象与文献[18]相符。

(a) 10 s时叶片轴向挥舞形变图

(b) 10 s时叶片周向摆振形变图

(c) 20 s时叶片轴向挥舞形变图

(d) 20 s时叶片周向摆振形变图

(e) 30 s时叶片轴向挥舞形变图

(f) 30 s时叶片周向摆振形变图

为直观反映塔影效应对叶片形变的影响，表2给出了模拟时间t为10 s、20 s和30 s时在有塔筒和无塔筒情况下3根叶片尖部的位移。无塔筒时3根叶片的叶尖位移基本保持一致，偏差不超过5%，在有塔筒干涉时塔筒正前方叶片叶尖的位移明显小于其余2根叶片，且另外2根叶片间也存在较大差异。

表2 叶片尖部位移

图7和图8分别是风力机叶轮的最大形变和最大等效应力曲线，在有塔筒时叶轮的形变和应力的波动程度和数值均高于无塔筒时。与无塔筒时相比，有塔筒时叶轮处于较不稳定的状态。当经过塔筒时，叶片都会受到一次塔筒干涉，这相当于在风力机运行过程中对叶轮施加了一个周期性的气动干扰，使叶轮运行的受力情况更为恶劣，这说明了塔影效应会对风力机组运行的稳定性产生负面影响。

图7 叶轮最大形变曲线

图8 叶轮最大等效应力曲线

2.3 流固耦合与非流固耦合结果对比

图9为流固耦合和非流固耦合下叶高为0.6R处叶片附近的流线图。由图9可以看出，流固耦合下流线的总体趋势与非流固耦合情况相似，但在叶片表面附近的流线变化较大。无塔筒时，在流固耦合下叶片前缘和后缘处的流线均出现明显弯折和回流现象；有塔筒时，叶片压力面甚至出现了涡旋，说明流固耦合下的塔影效应比非耦合情况更为明显。

(a) 非耦合，无塔筒

(b) 耦合，无塔筒

(c) 非耦合，有塔筒

(d) 耦合，有塔筒

无论有无塔筒，在流固耦合下叶片对流场的扰动均大于非耦合情况，造成此现象的原因可能是叶片的形变与位移使得叶片表面压力发生变化。为研究叶片表面的压力变化，图10给出了在耦合和非耦合情况下叶片表面的压力分布。由图10可知，无论是压力面还是吸力面，在耦合情况下叶片表面压力均略低于非耦合情况下，且在耦合情况下压力面的压力分布情况比非耦合情况下更复杂，这说明形变对叶片的气动受力具有明显影响。

3 结 论

(1)在靠近塔筒的过程中叶片的气动载荷会发生较大突变，以本计算为例，叶轮受到的周向转矩与轴向推力的波动幅度是无塔筒时的3倍以上。

(2)有塔筒时叶片的形变程度大于无塔筒时，且有塔筒时3根叶片的叶尖位移差异较大，在塔筒正前方叶片的叶尖位移明显小于其他2根叶片。

(3)在风力机运行过程中塔筒会不断对叶轮产生气动干扰，导致整个叶片的形变和应力要劣于无塔筒时，塔影效应对风力机运行的稳定性有不利影响。

(4)考虑流固耦合的计算更符合风力机的实际工作条件。与非流固耦合相比，在流固耦合下叶片对流场的扰动较大，叶片附近流场流线的涡旋现象更明显，且考虑流固耦合时叶片表面的压力分布更复杂。

(5)对于减小塔影效应的对策，现在较为主流的解决方法是独立变桨距控制，但学者们提出的控制策略大多为反馈控制，其变桨具有一定滞后性。为了探寻减弱塔影效应对风力机影响的方法，笔者曾尝试将塔筒替换为桁架并进行了分析，结果表明桁架式风力机的气动载荷波动程度小于塔筒式风力机。采用桁架可在一定程度上提高机组运行的稳定性与安全性，但由于桁架结构复杂，安装和维护成本较高，其实用性仍需要进一步探讨。

(a) 非耦合，压力面

(b) 耦合，压力面

(c) 非耦合，吸力面

(d) 耦合，吸力面

：

[1] 范忠瑶, 康顺, 赵萍. 上风向风力机塔影效应的数值模拟研究[J].工程热物理学报, 2012, 33(10): 1707-1710.

FAN Zhongyao, KANG Shun, ZHAO Ping. Numerical simulations of upwind wind turbine tower shadow effects[J].JournalofEngineeringThermophysics, 2012, 33(10): 1707-1710.

[2] CERMAK J E, HORN J D. Tower shadow effect[J].JournalofGeophysicalResearch, 1968, 73(6): 1869-1876.

[3] DAS S, KARNIK N, SANTOSO S. Time-domain modeling of tower shadow and wind shear in wind turbines[J].ISRNRenewableEnergy, 2011, 2011: 1-11.

[4] SESCU A, ANDERSEN B, AFJEH A A, et al. Computational investigation of tower shadow effects on wind turbines[C]//ASME2011InternationalMechanicalEngineeringCongressandExposition. Denver, Colorado, USA:American Society of Mechanical Engineers, 2011.

[5] 李少林, 张兴, 杨淑英, 等. 风力机动态模拟器的仿真研究[J].太阳能学报, 2010, 31(10): 1366-1372.

LI Shaolin, ZHANG Xing, YANG Shuying, et al. Simulation research of dynamic wind turbine simulator[J].ActaEnergiaeSolarisSinica, 2010, 31(10): 1366-1372.

[6] HSU M C, BAZILEVS Y. Fluid-structure interaction modeling of wind turbines: simulating the fullmachine[J].ComputationalMechanics,2012, 50(6): 821-833.

[7] BAZILEVS Y, HSU M C, TAKIZAWA K, et al. ALE-VMS and ST-VMS methods for computer modeling of wind-turbine rotor aerodynamics and fluid-structure interaction[J].MathematicalModelsandMethodsinAppliedSciences, 2012, 22(S2): 1-61.

[8] BAZILEVS Y, TAKIZAWA K, TEZDUYAR T E, et al. Aerodynamic and FSI analysis of wind turbines with the ALE-VMS and ST-VMS methods[J].ArchivesofComputationalMethodsinEngineering, 2014, 21(4): 359-398.

[9] JONKMAN J, BUTTERFIELD S, MUSIAL W, et al. Definition of a 5-MW reference wind turbine for offshore system development[R]. Colorado, USA: National Renewable Energy Laboratory, 2009.

[10] 潘旭. MW级风力发电机风轮叶片流固耦合场强度分析[D]. 郑州: 郑州大学, 2011.

[11] 何福添. 水平轴风力机叶片的设计及流固耦合计算[D]. 武汉: 武汉理工大学, 2013.

[12] 李媛. 风力机叶片流固耦合数值模拟[D]. 北京: 华北电力大学, 2013.

[13] 胡丹梅, 张志超, 孙凯, 等. 风力机叶片流固耦合计算分析[J].中国电机工程学报, 2013, 33(17): 98-104.

HU Danmei, ZHANG Zhichao, SUN Kai, et al. Computational analysis of wind turbine blades based on fluid-structure interaction[J].ProceedingsoftheCSEE, 2013, 33(17): 98-104.

[14] 方平治, 顾明, 谈建国, 等. 数值模拟大气边界层中解决壁面函数问题方法研究[J].振动与冲击, 2015, 34(2): 85-90.

FANG Pingzhi, GU Ming, TAN Jianguo, et al. Method to solve the wall function problem in simulation of atmospheric boundary layer[J].JournalofVibrationandShock, 2015, 34(2): 85-90.

[15] 金新阳, 杨伟, 金海, 等. 数值风工程应用中湍流模型的比较研究[J].建筑科学, 2006, 22(5): 1-5, 28.

JIN Xinyang, YANG Wei, JIN Hai, et al. Comparative study of turbulence model on prediction of wind pressure and wind velocity around rectangular buildings[J].BuildingScience, 2006, 22(5): 1-5, 28.

[16] 张群峰, 何鸿涛. 不同湍流模型数值模拟三维轴对称凸体分离流动的比较[J].科学技术与工程, 2009, 9(13): 3693-3697, 3703.

ZHANG Qunfeng, HE Hongtao. Numerical study on three-dimensional separated flow on an axisymmetric bump by different turbulent models[J].ScienceTechnologyandEngineering, 2009, 9(13): 3693-3697, 3703.

[17] 胡煜, 伍青安, 袁越, 等. 风电引起3p闪变的仿真分析[J].电力自动化设备, 2013, 33(3): 108-111.

HU Yu, WU Qing'an, YUAN Yue, et al. Simulative analysis of 3p voltage flicker caused by wind farm integration[J].ElectricPowerAutomationEquipment, 2013, 33(3): 108-111.

[18] 胡斌, 刘双. 风电场风速对风机叶片性能影响的模拟分析[J].动力工程学报, 2016, 36(12): 993-999.

HU Bin, LIU Shuang. Influence of wind speed on the performance of wind turbine blades[J].JournalofChineseSocietyofPowerEngineering, 2016, 36(12): 993-999.