海上风力机随机风场模拟及风振响应分析

柯世堂，王同光，曹九发，王珑（1. 南京航空航天大学 土木工程系，江苏 南京，10016；

海上风力机随机风场模拟及风振响应分析

柯世堂1， 2，王同光2，曹九发2，王珑2
（1. 南京航空航天大学 土木工程系，江苏 南京，210016；

2. 南京航空航天大学 江苏省风力机设计高技术研究重点实验室，江苏 南京，210016）

为研究海上风力机的风振特性，进行风力发电塔-轮系统的随机风场模拟和风振动力响应计算。采用谐波叠加法模拟塔架和风轮的来流风速时程，进而基于改进的叶素-动量理论（MBEM）模拟考虑风轮和塔架相干效应、风轮旋转效应的风轮脉动风速时程。结合已提出的柔性结构风振精细化频域计算方法“一致耦合法”，对海上风力发电系统结构进行风振动力响应和风振系数计算。研究结果表明：海上风力发电塔-轮系统的风振动力响应以共振效应为主，但背景响应和耦合项不能忽略，风振呈现多模态耦合和多振型响应2个显著特征；系统风振系数的分布差异较大，其中风轮尖部最大（2.35），塔架中下部位最小（1.40）。

海上风力机；随机风场模拟；风振响应；风振系数；一致耦合法

海上风电场技术相对于内陆风电场具有独特优势，未来必然会成为重要的可持续能源［1-3］。典型海上水平轴风力发电机组的显著特点是风轮、传动链、发电机等部件都放置在高耸细长的塔架之上，塔架的高度通常都要超过百米，且系统具有质量轻、阻尼小、自振频率低等结构特性，属于典型的风敏感结构［4-6］。针对海上风力发电机组的抗风研究，陈小波等［7］采用谐波叠加法进行了考虑空间相关性的海上风力发电机组随机风场模拟；OSAMU 等［8］采用有限元方法求解了风力发电塔架在风荷载和地震作用下的动力响应；李德源等［9］对海上风机塔架在风波联合作用下的动力响应进行了数值分析，获得了塔架的动态响应特征；徐亚洲等［10］建立了近海风力发电高塔的风浪相互作用模型，并对塔架在风浪作用下的随机响应进行分析；贺广零等［11］对风浪联合作用下的海上单桩基础风机进行动力响应分析。已有研究成果为海上风力发电机组的抗风设计提供了参考，但还有2个方面的问题值得进一步探讨：1） 海上风力发电塔-轮系统的精细化随机风场模拟，需要全面考虑到风轮和塔架之间的相干性、风轮自身的旋转效应；2） 海上风力发电系统结构的风振机理探讨和风振系数取值，特别是共振现象和耦合机理的研究。鉴于此，本文作者结合谐波叠加法和改进的叶素-动量理论，模拟考虑风轮和塔架之间相干效应、风轮自身旋转效应影响的海上风力发电塔-轮系统随机风场。建立某5 MW海上风力发电塔-轮系统“风轮-机舱-塔架”一体化有限元模型，基于作者已提出的风振精细化频域计算方法“一致耦合法”进行系统的风振动力响应、风振系数和耦合效应分析。

1　工程简介及风场模拟

某 5 MW三桨叶海上风力发电机组系统塔高为124 m，海平面高度为40 m，底径为4.8 m，顶径为2.6 m，塔体通长为变厚度结构。机舱长为12 m，宽为4.6 m，高为4.2 m，总质量为140.2×103kg。各风轮之间成120°夹角，沿周向平均分布，风轮直径为120 m，宽度为2.4 m，厚度为0.38 m，长度为60 m，偏航角为0°，额定转速为17 r/min。

海上风机脉动风场由2部分构成：塔架脉动风场和风轮脉动风场。图1所示为风力发电机组系统的风场图。从图1可以看出：在海上风机运营时风轮和塔架的风场之间相互影响，尤其是处于上风向的风轮对塔架的影响更为显著；而相比塔架的相干性，风轮本身的风场受到自身旋转效应和气弹效应影响更加明显。顺风向脉动风场能量相比侧向和垂直向脉动风场来说是风轮和塔架主要承受的载荷，因此本文主要进行风轮和塔架的顺风向脉动风速时程的数值模拟和风振动力响应计算分析。

图1　海上风力发电机组风场特性图Fig. 1　Wind characteristics of offshore wind turbine system

1.1塔架的脉动风速模拟

塔架和风轮的来流风速时程均采用谐波叠加法模拟［12］，该方法是对平稳随机过程进行离散化数值模拟的一种研究方法。由于风轮和塔架在3个方向上均存在相干性，采用 Davenport相关系数考虑风轮和塔架之间的相关性。

式中：Cx，Cy和Cz分别为风轮和塔架上任意两点横向、顺风向和垂直向的衰减系数；ω为脉动风频率；v（H）为H处平均风速。若只考虑垂直方向的相关性影响，Cx=Cy=0，则取Cz=10。在固定的脉动风频率下各点的相关性随着距离的增大逐渐变小。本文采用Davenport水平脉动风速谱为

式中：n表示模拟的频率取值；f=1 200n/U10， U10表示高度z=10 m的平均风速；u*为折算风速。

定义风机上风速模拟节点，假定均为零均值的平稳高斯过程，其风谱密度函数矩阵为

式中：Sii（ω）为节点脉动风自功率谱，采用式（2）中的风谱模型计算；Sij（ω）为互功率谱，其表达式需要用到塔架和风轮、塔架上任意两点之间的相干性，计算公式为

其中：风轮只需考虑旋转平面内各点间的相干性，风轮和塔架之间需要还需要考虑顺风向的相干性影响。再将S（ω）进行Cholesky分解得到

式中：H*（ω）T为H（ω）的共轭转置。H（ω）三维矩阵表达式为

此时风力发电塔架上的任何1个节点脉动风速时程可以由其功率谱决定，根据Shinozuka理论，模拟的风速时程可以表达为

其中：风谱在频率范围内划分成 N个相同部分；Δω=ω/N为频率增量；|Hjm（ωl）|为基于Davenport来流风谱矩阵进行Cholesky分解获得的下三角矩阵的模；θml为介于 0和 2π之间均匀分布的随机数，可采用Matlab的随机数生成函数，建议每次生成随机数后应恢复初始状态；ωl=l·Δω是频域的递增变量；ψjm（ωl）为2个不同作用点之间的相位角，是由Hjm（ωl）的虚部和实部的比值确定。

1.2风轮的脉动风速模拟

与塔架的风场不同的是，风轮的风速时程模拟需要考虑风轮自身的旋转效应和气弹效应。目前解决风轮气动性能有 3种方法［13］：叶素-动量理论、涡尾迹方法、CFD方法。其中涡尾迹方法适合模拟风力机风轮的复杂风场，能准确地计算出叶片风荷载的分布，虽然计算量比CFD小很多，但不能满足风力机日常快速计算的要求；CFD方法是最能精确计算风轮气动特性的方法，但是计算量太大、耗时太长。因此本文采用叶素-动量理论进行风轮气动载荷的模拟，实现气动和风轮结构的双向耦合。

叶素动量理论（BEM）是进行风力机气动载荷计算最为经典的方法［14］。不仅简便快捷，而且在具备准确风轮数据的条件下能够提供满意的计算结果。本文采用修正的BEM 理论［15］，引入叶根损失和叶尖损失，在轴向诱导因子较大时使用损失因子的经验模型，并加入动态入流和动态失速模型。使用该方法，可以计算风力机在不同风速、转速、桨距角及偏航角情况下的动态载荷。进而获取作用在风轮上的脉动风速时程。

根据BEM理论，风轮上的相对风速vrel采用下式计算：

式中：vox和 voy分别为沿顺风向和横风向的来流脉动风速，采用式（7）谐波叠加法计算；vrot为叶片旋转导致的线速度；W为诱导速度；vbx和vby分别为叶片振动速度。各个速度的关系如图2所示。

图2　流经某一风轮的局部速度三角形Fig. 2　Local velocity triangle through a wind wheel

诱导速度W可由下式表示：

式中：B为叶片数；L是指升力；φ为入流角；ρ为空气密度；r为叶片截面的展向位置；n为推力方向的单位向量；F为普朗特叶尖损失因子；fg为Glauert修正。同时，本文还采用动态入流模型和动态失速模型，修正风轮运转的非定常效应。

根据下式计算风轮攻角α

式中：β为桨矩角；θtwist为风轮剖面几何扭角，计算公式为

通过风轮翼型插值方法，可以得到升力系数 Cl和阻力系数Cd，从而计算出升力L和阻力D

这样得到风轮的法向载荷Fn和切向载荷Ft

综上所述，基于式（7）采用Matlab程序模拟风轮的来流风速时程，再采用式（9）计算每个风速时程样本对应的诱导速度，如此循环计算最终获得风轮上的脉动风速时程。同时也可以计算得到风轮的升力系数和阻力系数，以及受到的法向载荷和切向载荷分布。

1.3风力发电塔-轮系统随机风场模拟结果

基于上述风场模拟方法，采用Matlab语言编制相应的数值模拟程序。在计算过程中脉动风上限频率取为 2π，脉动风频率分割点数取为 2 048，频率增量Δω=0.003 07 Hz，取当地10 m高平均风速为24 m/s。来流脉动风速谱和相干函数均取式（1）和式（2）的Davenport模型。图3和图4所示分别为风轮和塔架中部2点的顺风向脉动风速时程和功率谱曲线图。其中：坐标采用对数坐标，对比谱为Davenport风谱。

图3　风轮脉动风速时程模拟结果Fig. 3　Simulating result of fluctuating wind velocity of rotor

图4　塔架中部脉动风速时程模拟结果Fig. 4　Simulating result of fluctuating wind velocity of tower

从图3和图4可以看出：风轮的脉动风速功率谱曲线在高频处存在较大的能量和脉动特性，应该是由风轮的旋转和气弹效应引起的高频能量，在风力发电机组的风振动力分析中应引起重视；塔架中部的脉动风速模拟曲线和Davenport风谱吻合较好，但在高频处由于受到风轮相干性的影响数值有微弱的波动。因此根据本文模拟过程和对比分析可知：采用本文方法可以很好地模拟考虑风轮和塔架相干效应、风轮自身旋转的脉动风速时程，为进一步的风振动力响应提供输入参数。

2　海上风力发电塔-轮系统动力特性

基于ANSYS软件平台，建立了海上风力发电塔-轮系统的“风轮-机舱-塔体”一体化有限元模型。其中风轮和塔体采用 SHELL91单元，机舱及其内部结构可作为整体采用梁单元BEAM189模拟。通过多点约束单元耦合命令将各部分连接在一起，形成整体的海上风力发电机组系统。依据效率和精度均衡的原则，模型一共划分了2 436个单元。模态分析时把风轮旋转产生的离心力作为预应力均匀施加在风轮上，计算的频率和模态信息均考虑风轮转动带来的离心力效应。

图5和图6所示分别为塔-轮系统和风轮各自的模态振型，图7所示为考虑/不考虑风轮离心力作用2种工况下结构前100阶自振频率的分布曲线。

从图5～7可以看出：基于ANSYS的塔-轮系统一体化建模时可以同时考虑风轮和塔架的耦合模态。当考虑风轮转动引起的离心力作用时，塔-轮系统的基频要略大于不考虑离心力作用下的系统频率，并且随着模态数目的增加，离心力效应带来的频率影响越来越大。本文的后续计算均采用考虑风轮离心力作用的更加真实的模态参数。系统的基频（0.29 Hz）很低，第100阶模态频率为21.27 Hz，模态之间间隔较小。

图5　风轮的第1阶模态振型示意图Fig. 5　Sketch of first order vibration mode ofwind rotors systems

图6　风轮的第1阶模态振型示意图Fig. 6　Sketch of first order vibration mode of wind rotors systems

图7　塔-轮系统固有频率分布图Fig. 7　Scattergram of natural frequence for ower-rotor system

3　海上风力发电塔-轮系统风振特性

采用柯世堂等［16］提出的完全考虑结构背景、共振及背景和共振模态之间交叉项的一致耦合方法（CCM）来计算海上风力发电塔-轮系统风振响应，可以很好地考虑风轮和塔架之间的耦合模态、系统各共振模态之间的耦合效应。图8所示为CCM方法的计算流程图，具体推理及验证过程见文献［16］。

图9和图10所示分别为风轮尖部和塔架顶部位移响应功率谱密度函数曲线。图9和图10中目标响应功率谱密度函数曲线由背景和共振分量2部分构成，且共振响应明显占据主导地位；风轮和塔架作为一个系统承受考虑旋转和气弹效应的随机复杂风场时，会激发系统多个模态的共振效应；且与一般高耸/高层建筑不同的是，风振响应中第1阶模态激发的共振能量并不是最大，而分别是由第9阶和第11阶模态激发的共振响应最大，其对应振型是风轮的前后舞动和塔架的弯曲变形耦合形态，除此之外，系统第1，8和12阶模态也可能激发共振效应。海上风力发电塔-轮系统的风振响应分析时需要考虑多阶风轮和塔体的耦合振型，建模时应该采用风轮和塔体一体化模型。

图11所示为风力发电塔-轮系统5个不同目标风振响应中背景、共振及背景和共振耦合分量的分布图。从图11可见：横坐标为5个节点编号，纵坐标表示各分量对位移响应的贡献数值，值得注意的是耦合分量数值有正负之分，负值表示当忽略改耦合效应时就高估其数值，正值则相反。5个节点分别为风轮尖部位移/m、风轮根部位移/m、塔架顶部位移/m、塔架中部轴力/103N和塔架底部弯矩/104Nm。

图8　CCM方法计算流程示意图Fig. 8　Flow sketch map of CCM

图9　风轮尖部位移响应功率谱图Fig. 9　PSD of displacement responses for rotor point

从图11可以发现：对于海上风力发电机组系统的不同目标风振响应，共振分量均要大于背景分量和背景和共振之间耦合项；背景响应在塔架底部和中部所占比例较高，而在风轮的风振响应中要远远小于共振分量，这在风轮系统的抗风设计时要加以重视。

图10　塔架顶部位移响应功率谱图Fig. 10 PSD of displacement responses for tower top

风振系数是海上风力发电塔-轮系统抗风设计的关键参数，是工程设计人员最容易理解和应用的设计思路。目前，一般按照 GB50135—2006 “高耸结构设计规范”［17］的经验公式进行计算，获得一个单一常数来进行整个系统的抗风动力设计。然而在很多风敏感柔性结构（大跨空间结构、冷却塔、输电塔等）中已经证明采用这个单一参数来进行整体结构设计并不合理，特别对于海上风力发电机组来说，风轮和塔架本身就是受力形式不同的结构，相关风机设计规范也没有给出精确的分区风振系数取值，因此本文进行海上风力发电塔-轮系统的风振系数取值探讨。其中风振系数的计算公式为

式中：Ri，Rei和Rfi分别指节点i的总响应、平均响应和脉动响应；g为保证系数，或称为峰值因子。当脉动风响应的概率分布为正态分布时，g可表示为

式中：T为最大值相应的时距，我国荷载规范规定平均风的时距为10 min，因此T取600 s；γ为欧拉常数，通常取0.577 2；v为水平跨越数，可通过风振动力响应时程计算获得。

图11 典型节点目标响应各分量贡献图Fig. 11 Contributions of respective component for typical nodes

图12所示为海上风力发电塔-轮系统的风振系数数值三维分布图，并在3个风轮尖部、塔架顶部、海平面部位和底部标出了风振系数。从图12可知：风力发电塔-轮系统各个不同的部位计算得到的风振系数相差较大，若采用统一常数取值并不合理，而且由于风轮受到的风场更加复杂，自身结构也比较柔，其风振系数要比塔架的大；而塔架不同部位的风振系数浮动也比较明显，例如塔顶风振系数较小，仅为1.59，而塔顶部位风振系数达到1.96，建议在设计时应该分段取值。

图12　海上风力发电塔-轮系统风振系数分布图Fig. 12 Wind vibration coefficient of offshore wind turbine tower system

4　结论

1） 采用本文提出的基于谐波叠加法和改进的叶素动量理论方法可以有效地模拟风轮和塔架的脉动风速时程。

2） 海上风力发电机组的动力特性分析需要考虑风轮和塔架的耦合模型，且风轮旋转引起的离心力会增大风力机系统的结构频率，并且随着模态数目的增大影响愈加明显。

3） 海上风力机系统的风振动力响应以共振分量为主，背景分量和交叉耦合项不能忽略，主要呈现多模态耦合和多振型响应2个特性。系统风振系数分布数值差异较大，其中风轮的风振系数最大可达到2.35，塔架的风振系数在1.4～2.0之间，且塔顶最大。

［1］ NAESS A， GAIDAI O， HAVER S. Efficient estimation of extreme response of drag-dominated offshore structures by Monte Carlo simulation［J］. Ocean Engineering， 2007， 34（16）：2188-3197.

［2］ RONOLD K O， LARSEN G C. Optimization of a design code for wind turbine rotor blades in fatigue［J］. Engineering Structure，2001， 23（8）： 993-1002.

［3］ KE S T， WANG T G， GE Y G， et al. Wind-induced responses and equivalent static wind loads of tower-blade coupled large wind turbine system［J］. Structural Engineering and Mechanics，an International Journal， 2014， 52（3）： 485-505.

［4］ KE S T， GE Y J， WANG T G， et al. Wind field simulation and wind-induced responses of large wind turbine tower-blade coupled structure［J］. The Structural Design of Tall and Special Buildings， 2015， 24（8）： 571-590.

［5］ 柯世堂， 王同光， 曹九发， 等. 考虑土-结相互作用大型风力发电结构风致响应分析［J］. 土木工程学报， 2015， 48（2）： 35-44. KE Shitang， WANG Tongguang， CAO Jiufa， et al. Analysis of wind-induced responses on large wind power structuresconsidering soil-structure interaction［J］. Journal of Civil Engineering， 2015， 48（2）： 35-44.

［6］ 柯世堂， 王同光， 陈少林， 等. 大型风力机全机风振响应和等效静力风荷载［J］. 浙江大学学报（工学版）， 2014， 48（4）：686-692. KE Shitang， WANG Tongguang， CHEN Shaolin， et al. Wind-induced responses and equivalent static wind load of large wind turbine system［J］. Journal of Zhejiang University （Engineering Science）， 2014， 48（4）： 686-692.

［7］ 陈小波， 陈建云. 海上风力发电塔脉动风速时程数值模拟［J］.中国电机工程学报， 2008， 28（32）： 111-116. CHEN Xiaobo， CHEN Jianyun. Numerical simulation of fluctuating wind velocity time series of offshore wind turbine［J］. Proceedings of the Chinese Society of Electrical Engineering，2008， 28（32）： 111-116.

［8］ OSAMU K， TATSUOMI R. Dynamic response analysis of onshore wind energy power units during earthquakes and wind［C］//Proceedings of the twelfth international offshore and polar engineering conference. Kitakyushu， Japan： The International Society of Offshore and Polar Engineers， 2002：520-527.

［9］ 李德源， 刘胜祥. 风波联合作用下的风力机塔架疲劳特性分析［J］. 太阳能学报， 2009， 30（10）： 1250-1256. LI Deyuan， LIU Shengxiang. Fatigue characteristics analysis of wind turbine tower under wind-wave combined effect［J］. Acta Energiae Solaris Sinica， 2009， 30（10）： 1250-1256.

［10］ 徐亚洲， 李杰. 风浪相互作用 Stokes模型［J］. 水科学进展，2009， 20（2）： 281-286. XU Yazhou， LI Jie. Stokes model for wind-wave interaction［J］. Advances in Water Science， 2009， 20（2）： 281-286.

［11］ 贺广零， 仲政. 风浪联合作用下的海上单桩基础风力发电机组动力响应分析［J］. 电力建设， 2012， 33（5）： 1-7. HE Guangling， ZHONG Zheng. Dynamic responses analysis of offshore wind turbine systems with monopole foundation under combined wind & wave loads［J］. Electric Power Construction，2012， 33（5）： 1-7.

［12］ KAREEM A. Numerical simulation of wind effects：a probabilistic perspective［J］. Journal of Wind Engineering and Industrial Aerodynamics， 2008（10）： 1472-1497.

［13］ 赵峰， 段巍. 基于叶素-动量理论及有限元方法的风力机叶片载荷分析和强度计算［J］. 机械设计与制造， 2010（8）： 42-45. ZHAO Feng， DUAN Wei. Loading analysis and strength calculation of wind turbine blade based on blade element momentum theory and finite element method［J］. Machinery Design & Manufacture， 2010（8）： 42-45.

［14］ 伍艳， 谢华， 王同光. 风力机叶片的非定常气动特性计算方法的改进［J］. 工程力学， 2008， 25（10）： 54-60. WU Yan， XIE Hua， WANG Tongguang. Modification of calculating unsteady aerodynamic characteristics of wind turbine blades［J］. Engineering Mechanics， 2008， 25（10）： 54-60.

［15］ WANG T G， COTON F N. Prediction of the unsteady aerodynamic characteristics of horizontal-axis wind turbines including three-dimensional effects［J］. Proceedings of the Institution of Mechanical Engineers Part A： Journal of Power and Energy， 2000， 214（A5）： 385-400.

［16］ 柯世堂， 葛耀君， 赵林， 等. 一致耦合方法的提出及其在大跨空间结构风振分析中的应用［J］. 中南大学学报 （自然科学版），2012， 43（11）： 4457-4463. KE Shitang， GE Yaojun， ZHAO Lin， et al. Proposition and application of consistent coupling method in wind-induced response of long span structures［J］. Journal of Central South University （Science and Technology）， 2012， 43（11）： 4457-4463.

［17］ GB 50135—2006. 高耸结构设计规范［S］. GB 50135—2006. Code for design of high-rising structures［S］.

（编辑 罗金花）

Simulation of stochastic wind field and
wind-induced responses of offshore wind turbines

KE Shitang1， 2， WANG Tongguang2， CAO Jiufa2， WANG Long2

（1. Department of Civil Engineering， Nanjing University of Aeronautics and Astronautics， Nanjing 210016， China；2. Jiangsu Key Laboratory of Hi-Tech Research for Wind Turbine Design，Nanjing University of Aeronautics and Astronautics， Nanjing 210016， China）

In order to analyze the wind-induced characteristics of offshore wind turbine system， the simulation of stochastic wind field and wind-induced responses of wind turbine tower-rotor system were finished. The fluctuating wind velocity time series was simulated by the modified harmony superposition method. The incoming fluctuating wind velocity time series of wind rotors was simulated by the harmony superposition method， then the fluctuating wind velocity time series considering the rotor rotational effect and aero-elastic effect was simulated with modified blade element momentum. Finally， using the domain calculating method “consistent coupled method”， the wind-induced responses and coefficients for the offshore wind turbine system were calculated. The results show that the resonant effect of wind-induced dynamic responses is obvious for offshore wind turbine system， and is characterized with multimode coupling effect and multiple mode of vibration responses； the wind vibration coefficients are fluctuating in different components， and the value of wind vibration in rotor top is 2.35， the middle section of tower is 1.40.

offshore wind turbine system； stochastic wind field simulation； wind-induced response； wind vibration coefficient； consistent coupled method

TK83；TU279.7+4

A

1672-7207（2016）04-1245-08

10.11817/j.issn.1672-7207.2016.04.022

2015-04-09；

2015-06-09

（Foundation item）：国家重点基础研究发展（973计划）项目（2014CB046200）；国家博士后科学基金资助项目（2015T80551）（Project （2014CB046200） supported by the National Basic Research Development Program （973 Program） of China； Project （2015T80551） supported by the National Science Foundation for Post-doctoral Scientists of China）

柯世堂，博士，副教授，从事结构风工程研究；E-mail：keshitang@163.com