施加一阶挥舞振型后的风轮气动特性研究

施加一阶挥舞振型后的风轮气动特性研究
王志光;高志鹰;汪建文;陈天阁;孙晓颖
【摘要】文章采用Fluent软件,对半径为0.7 m的水平轴风力机进行了三维非定常气动特性数值模拟,应用动网格技术,为风轮模型施加了一阶挥舞振动.对比分析有无振动情况下的计算结果,得到了振动叶片表面动压、转矩及轴向推力的周期变化规律,并发现:当叶片振动后,动压变化频率约为振动频率的2倍;转矩和推力较振动出现相位差,且推力与振动相位差0.3π.通过为风轮叶片施加一阶振动,实现了对振动状态下风力机气动特性的初步探索,为更加准确地模拟风力机气动特性提供了方法参考.%The Fluent is used to simulate 3D unsteady aerodynamic characteristics of a horizontal axis wind turbine of R=0.7 m.A one order waving vibration is applied to the wind turbine model using the dynamic mesh method.After comparative analysis of the simulation results between vibrating and novibrating conditions,change law of blade surface's dynamic pressures,torque and axial thrust are obtained,and it is found out that while the blades are vibrating,the frequency of dynamic pressure change is approximately two times higher than that of vibration change,phase difference exists between vibration and torque or thrust,the phase difference between vibration and thrust is about 0.3π.A preliminary exploration of wind turbine aerodynamics characteristics under vibration condition is obtained after one order vibration is applied to wind turbine,and a method for more accurate simulation of wind turbine is provided also in this article.
【期刊名称】《可再生能源》
【年(卷),期】2014(032)006
【总页数】5页(P787-791)
【关键词】水平轴风力机;风轮;气动特性;一阶振型;动网格
【作者】王志光;高志鹰;汪建文;陈天阁;孙晓颖
【作者单位】内蒙古工业大学能源与动力工程学院,内蒙古呼和浩特010051;内蒙古工业大学能源与动力工程学院,内蒙古呼和浩特010051;风能太阳能利用技术省部共建教育部重点实验室,内蒙古呼和浩特010051;内蒙古工业大学能源与动力工程学院,内蒙古呼和浩特010051;风能太阳能利用技术省部共建教育部重点实验室,内蒙古呼和浩特010051;内蒙古工业大学能源与动力工程学院,内蒙古呼和浩特010051;辽宁省能源研究所,辽宁营口 115003
【正文语种】中文
【中图分类】TK83
0 引言
风轮在运转的过程中，由于受到外界风载荷的作用，叶片会产生振动。

振动叶片不仅会对流经的空气产生扰动，还将影响到风力机的气动特性。

因此，对振动情况下的风力机气动特性进行研究，不仅能够了解风轮在实际运行过程中风能吸收能力及尾迹流场特性，还能对整机性能进行评估。

风力机气动性能的实验研究耗资巨大，且对实验环境要求较为苛刻。

随着计算机的快速发展，计算精度的不断提高，对风机气动特性模拟则体现出前所未有的优势[1]。

目前，在风力机气动特性的模拟中，多使用多参考系模型和滑移网格模型来实现风轮的旋转，计算过程中未考虑叶片的振动[2]～[4]。

少数学者考虑了振动问题，但
仅对叶片某一截面翼型在振动情况下的气动性能进行了二维数值研究。

2012年，德国的Sayed Mohamed A[5]应用非定常的雷诺平均N-S方程（URANS），对NREL改进型风力机翼型，在非定常低速流动中的空气动力学特
性进行了数值模拟。

2011年，北京航空航天大学的钟国华[6]应用侵入边界法建立了振动翼型的数学模型，并就浅层次的流固耦合问题给出了初步的分析结果。

2009年，布拉格查里大学的L.Dubcová·M.Feistauer[7]应用动网格技术，对既有俯仰运动又有往复运动的翼型模拟研究了二维粘性流动与振动翼型之间的相互作用。

由于风轮转动过程中，振动出现在叶片的各个表面，因此仅对某截面进行考虑显然无法真实反映出风轮的整体气动性能。

本文采用商业CFD软件Fluent 6.3，应用
动网格技术，模拟了水平轴风力机在运行过程中，一阶挥舞振动下的叶片表面动压力、转矩及轴向推力，初步揭示了叶片振动对风力机气动特性的影响。

通过模拟一阶振动状态下的风轮气动特性，为进一步研究实际振动形态下的风力机气动特性提供方法，为风力机整体气动性能做了数据补充。

1 模型建立
1.1 GAMBIT建模
本次计算所采用的风轮模型与实际风轮比例为1:1，风轮直径为1.4 m，叶片数为
3个，定义X轴为风轮旋转轴。

计算域模型如图1所示。

主要包括加密区①和非加密区②，加密区包含着风轮旋
转域③。

在风轮旋转域③内应用动网格技术，以实现叶片的一阶挥舞振动。

计算区域的尺寸选择考虑了壁面边界层及入口、出口对风轮尾迹流场的影响，设定计算域直径D=3 m，计算域长L=11 m，风轮距入口段4 m，距出口段7 m。

为了减小
叶片变形对网格质量的影响，模型选用非结构四面体网格进行划分，总网格数量约350万。

叶片转动控制方法选用滑移网格技术，在额定风速为8 m/s的情况下，
选取叶尖速比λ=5，转速为545 r/min。

图1 计算域模型Fig.1 Computational domain model
1.2 数学模型
流动过程中的控制方程主要由连续方程和动量方程（N-S方程）来描述。

连续方程：
动量方程：
2 振动控制方法
叶片振动控制方法选用动网格技术（Dynamic Mesh Method），叶片在挥舞面
上做一阶对称模态振动。

叶片振动控制方程推导如下。

设叶片上各节点振动位移随时间的变化关系为
式中：X为叶片上某节点在挥舞面上的实际位移，m；A为该节点的振幅，m；ω为叶片振动圆频率，r/s，实验测得叶片一阶挥舞振动固有频率f=20 Hz，故此处
取ω=40π。

叶片上各节点的振幅随径向距离的增大而增大，振型方程由叶片静态振动实验数据确定。

假设风轮叶尖处振幅为±5 mm，叶片与轮毂连接处振动位移为0。

最终振动方程为
式中：R为叶片上该节点距离旋转轴的距离，R=m；y和z分别为叶片上节点的 Y、Z坐标；ai为待定系数，由一阶振动特征向量确定。

ai取值：a1=0.124 4，a2=-0.255 8，a3=0.201 8，a4=-0.05918，
a5=0.007811，a6=-0.0003201。

图2所示为叶片振动效果放大图。

图2 振动效果放大图Fig.2 Vibration effect enlargement
图 2（a）为ωt=36π时，叶片处于平衡位置时的状态；图 2（b）为时，叶片处
于正半轴最值点时的状态；图 2（c）为时，叶片处于负半轴最值点时的状态。

3 结果分析
本文在来流风速v=8 m/s，尖速比λ=5的工况下，对比分析了叶片在振动与不振动时的计算结果。

选取叶片振动到负半轴最值点后一个周期内的全部数据进行分析，以得到叶片振动后对风轮气动特性影响的全部信息。

3.1 叶片表面动压分析
当气流流经风机时，由于受到叶片的影响，导致在叶片的压力面和吸力面附面层内，气流的流速出现差异，此时作用于叶片表面的动压值也出现相应的变化。

由于动压值等于单位质量流体所具有的动能，因此对叶片表面动压值的把握能够了解到叶片在运行过程中，叶片表面气流流速的准确信息。

图3所示为叶片表面动压随时间变化图。

图中3条曲线分别代表了无振动时动压
线（DP-ST线）、振动时动压线（DP-SH线）及振动速度线（Velocity线）。

图3 叶片表面动压随时间变化图Fig.3 Blade surface dynamic pressure vs.time 对比曲线DP-ST和DP-SH后可以发现，当叶片不振动时，叶片表面动压最大值
较稳定，该值在附近变化；当叶片振动后，叶片表面动压值随时间出现周期变化，最小值为最大值为随叶片振动，叶片表面动压变化最大可达171.15 Pa。

对比曲线DP-SH和Velocity可得动压随振动速度的变化规律：当叶片由向运动时，振动速度与来流方向相同，随叶片振动速度先增大再减小，作用于叶片表面的动压值则先减小后增大，当振动速度达到最大值时，动压值最小；当叶片由向ωt=运
动时，振动速度与来流方向相反，动压随振动速度的变化出现先减小后增大的变化趋势。

总体而言，叶片表面动压值与振动速度幅值呈反相关关系：当速度幅值最大时动压最小，而速度幅值最小时动压最大。

以此导致动压的变化频率近似为速度变化频率的2倍。

3.2 风轮转矩分析
风轮的转矩表征着驱使风轮转动的力的大小。

通常情况下，转矩的大小直接决定了风力机的轴功率。

在风轮的实际运行过程中，叶片的振动将对施加在风轮上的转矩产生影响。

本文对比了振动与无振动情况下风机所受转矩的变化，以此揭示叶片振动对风机转矩的影响。

图4 风轮转矩随时间变化图Fig.4 Wind turbine torque vs.time
图4所示为风轮转矩随时间变化图。

图中3条曲线分别代表了无振动时转矩线（T-ST线）、振动时转矩线（T-SH线）及振动速度线（Velocity线）。

对比曲线T-ST和T-SH可知：当叶片不振动时，风轮转矩较为稳定，此时转矩值近似为TT-ST=4.05 N·m；叶片振动后，气流受到振动叶片的影响，作用于风轮上的转矩相应出现周期变化，其转矩最大值为转矩最小值为3.71 N·m。

对比曲线T-SH与Velocity可得风轮转矩随叶片振动速度的变化规律：当叶片振动后，随着叶尖振动速度按余弦规律变化，作用于风轮上的转矩出现近似的余弦变化趋势，且其变化趋势与振动速度变化趋势相同。

转矩变化周期与振动周期相同，相位较振动有少许增加。

造成该现象的原因：以风轮为参照物，当叶片振动后流经叶片的等效风速发生变化，由此导致相对风速和攻角发生变化。

3.3 风轮轴向推力分析
当气流流经风轮时，受风压的影响，风轮会受到与来流风速同向的推力，而该推力将由塔架承受。

在风机实际运行过程中，风机所受的脉动推力将导致在塔架上出现
相应的脉动应力。

因此对风轮轴向推力的研究将有助于了解塔架的实际脉动受力情况。

本文对比了振动和无振动情况下风轮的轴向推力，以此揭示叶片振动对风轮轴向推力的影响。

图5为风轮轴向推力随时间变化图。

图中3条曲线分别代表了无振动时推力线（F-ST线）、振动时推力线（F-SH线）及振动速度线（Velocity线）。

图5 风轮轴向推力随时间变化图Fig.5 Wind turbine axial thrust vs.time
对比曲线F-ST和F-SH可知：当叶片不振动时，推力大小在FF-ST=4.05 N附近小幅变动；叶片振动后，风轮轴向推力出现周期变化，推力最大值为最小值为不同时刻推力变化最大可达ΔF=7.59 N。

对比曲线F-SH和Velocity可得轴向推力随振动速度的变化规律。

就总体趋势而言，推力变化周期与振动周期相同，其变化趋势与振动速度呈反相关关系。

当叶片振动后，随着叶尖振动速度按余弦规律变化，作用于风轮上的轴向推力变化趋势近似与振动速度变化趋势相反。

观察推力曲线和速度曲线的极值可发现：当振动速度最大值时，所对应的推力尚未达到最小值，而是向后有所推移。

计算得到推力与速度的相位差为φ=0.3π。

这种现象是由于叶片的变形对近壁面气流导流情况改变所造成的。

因此体现在塔架上的脉动应力与叶片振动速度相位相差0.3π。

4 结论
本文应用结合动网格技术的大涡模拟方法，在风速v=8 m/s、叶尖速比λ=5的条件下，对比分析了有无振动情况下水平轴风力机叶表动压力、转矩及轴向推力得到以下结论。

①当为叶片施加一阶挥舞振动后，与无振动计算结果相比，表面动压、转矩及推力均出现周期变化。

②叶片振动后，风轮表面动压值与振动速度幅值呈反相关关系，动压变化频率为振
动频率的2倍。

③叶片振动后，风轮转矩及推力变化周期与振动周期相同；转矩与振动相比，相位出现少许差异；推力与振动相比，相位相差0.3π。

本文初步探索了一阶模态振动下运行风轮的气动特性。

而在实际运行中叶片的挥舞振动是多阶振动形式的耦合，因此对于挥舞振动下风轮的气动特性研究需要在后续的工作中逐步实现。

参考文献：
[1]Hansen M O L.风力机空气动力学[M].北京：中国电力出版社，2009.16-22.
[1]Hansen M O L.Aerodynamics of Wind Turbine[M].Beijing： China Electric Power Press，2009.16-22.
[2]Thomas Hahm，Jurgen Kroning.In the wake of a wind turbine[N].Fluent NEWS，2000-4-26.
[3]WuBow S， Sitzki L， Hahm T.3D-simulation of the turbulent wake behind a wind turbine[A].The Second Conference on the Science of Making Torque from Wind[C].Lyngby:IOP Publishing，2007.12-33.
[4]Sezer-Uzol N，Gupta A，and Long L.N.3-D Time Accurate Inviscid and Viscous CFD Simulations of Wind Turbine Rotor Flow Fields[A].Proceedings of PARALEL CFD 2007，International Conference on Parallel Computational Fluid Dynamics[C].Antalya:Springer Berlin
Heidelberg,2007.457-464.
[5]Sayed， Mohamed A， Kandil， et al.Dynamic Stall Analysis of Horizontal-Axis-Wind-Turbine Blades Using Computational Fluid Dynamics[Z].Melaka:MA.4th International Meeting of Advances in Thermofluids（IMAT）， 2011.
[6]Zhong Guohua， Du Lin， Sun Xiaofeng.Numerical investigation of an oscillating airfoil using immersed boundary method [J].Journal of Thermal Science,2011，20（5）：413-422.
[7]Lenka Dubcová， Miloslav Feistauer.Numerical simulation of interaction between turbulent flow and a vibrating airfoil[J].Comput Visual Sci,2009，12 （5）：207-225.
[8]Whale J，Anderson C G.An experimental and numerical study of the vortex structure in the wake of a wind turbine[J].Journal of Wind Engineering and Industrial Aerodynamics， 2000， 84（1）：1-21.
[9]张立茹，汪建文，曲立群.风力机叶尖有无小翼三维流场的数值研究[J].太阳能学报，2011，32（11）：1599-1603.
[9]Zhang Liru， Wang Jianwen， Qu Liqun.A three-dimensional numerical simulation on flow field of wind turbine with and without tip vane[J].Acta Energiae Solaris Sinica，2011，32（11）：1599-1603.。