考虑气动阻尼的浮式风机频域响应分析

摘 要：分析了气动阻尼对浮式风机频域响应的影响.选取美国可再生能源实验室（NREL）提出的5兆瓦（MW）浮式风机模型作为算例，利用气动阻尼计算方法建立气动阻尼矩阵，再基于三维势流理论计算浮式平台的水动力系数，并将系泊系统视为线性弹簧以考虑其刚度，最后在频域内分别建立并求解考虑与不考虑气动阻尼两种情况下的浮式风机刚体运动方程.利用求解频域方程得到的幅频响应算子（response amplitude operators，RAOs）及结合JONSWAP海浪谱得到的响应谱，在频域内分析了气动阻尼对浮式风机刚体运动的影响.结果表明：作业工况下气动阻尼能有效地降低纵荡和纵摇运动RAOs的峰值，且能在一定范围内减小对应自由度上响应谱的幅值和零阶矩.

关键词：海上风机；频域分析；气动阻尼；幅频响应算子；响应谱

中图分类号：P752 文献标志码：A

海上风能有着风速大、湍流度低、储量丰富等特点，采用浮式风机是开发深水风能的首选技术手段.利用频域分析得到的RAOs及响应谱来评估浮式平台的动态性能，是其设计和优化中的重要内容：Bulder[1]等针对适用于5MW风机的三浮体式浮式平台，通过求解频域内浮式风机运动方程，得到了平台各运动自由度的RAOs及特定海况下的响应谱，并由此分析了该浮式风机在对应海洋环境下的适用性；唐友刚[2]等综合驳船式和单柱式（Spar式）平台的特点，提出了一种新型浮式平台方案，在频域内利用RAOs分析了波浪入射角和水深等因素对该平台运动的影响.而气动阻尼是多数学者在频域分析中忽略的因素.

Kühn [3]和Salzmann [4]指出，气动阻尼是海上风机支撑结构的振动与疲劳分析中的关键影响因素.邓露[5]等指出风机振动过大会对运行造成严重影响.Karimirad[6]等通过谱分析指出气动阻尼能有效抑制风机机舱在纵荡方向上的低频响应，且根据叶素动量理论指出浮式平台的刚体运动会受到气动阻尼的影响.但目前关于气动阻尼对浮式平台运动频域响应影响的研究较少：Roddier[7]等基于WindFloat浮式风机，通过模型实验得到浮式风机的RAOs，结果表明在研究中考虑与不考虑风荷载得到的RAOs相差不大，但该研究所涵盖的波浪周期和风速范围有限，且实验模型采用圆板代替风轮，不能准确地反映风机的气动特性；Wayman[8]等和Ramachandran[9]等联合利用分析软件WAMIT和FAST计算了多种浮式风机的RAOs并结合ISSC海浪谱计算了不同海况下各运动自由度的标准差，发现浮式风机不同于船舶和海洋平台，其RAOs的取值与海况有关，但该研究涵盖的工况较少，也没有深入探究影响气动阻尼的主要因素.

针对上述问题，本文在频域内分别建立并求解考虑与不考虑气动阻尼两种情况下的浮式风机运动方程.对比两种情况下RAOs和响应谱，在频域内按照不同的工况定量分析了气动阻尼对浮式风机响应的影响，并明确了影响气动阻尼的关键因素.

1 气动阻尼

1.1 气动阻尼的产生机理

气动阻尼源于风轮与空气的相互作用，如图1 [3]：叶素截面处于转速为Ω的风轮的径向r处，塔顶顺风向运动时，叶素处实际的轴向风速等于入流轴向风速V（1-α）减去塔顶速度top，导致叶片攻角α变小，在附着流的前提下并考虑升力系数CL，这将使得叶素受到的升、阻力变小，对应受到的推力dFx减小ΔdFx，整体上表现为风轮推力变小，故阻碍塔顶顺风向运动；同理，当塔顶逆风向运动时，风轮推力增大，进而阻碍塔顶的逆风向运动，则推力的变化始终阻碍着塔顶运动，这便是气动阻尼的机理.

1.2 气动阻尼的计算

利用NREL风机时域分析软件FAST建立风机模型，能考虑失速及小叶尖速比的情况，并能模拟具有复杂翼型的叶片[3-4].在FAST弹性动力模块中将结构阻尼设为零，并通过在软件中关闭浮式平台刚体运动自由度来限制平台运动以排除附加阻尼的影响，得到仅在气动阻尼作用下塔顶位移衰减曲线.

在FAST的气动模块中采用叶素动量理论计算气动荷载并考虑叶尖和叶根损失和采用Beddoes[10]等提出的动态失速模型.在FAST的控制模块中采用固定的风轮转速和浆距角，其具体值应根据不同的稳态风速确定[11]，并用指数模型来考虑竖直方向风剪切[12].通过时域模拟，得到了仅在气动阻尼作用下的塔顶振动衰减曲线，如图2所示（图中塔顶位移从对应风速下塔顶振动的平衡位置算起）.进而利用衰减曲线并结合式（1）和（2）计算对数衰减率和阻尼比 [3]：

式中：δ为塔顶位移的对数衰减率；ξaero为气动阻尼比；An，An+1 为相隔一个周期的两个位移峰值.

将风轮和机舱简化为塔架顶部的刚性质点，把塔架机舱转子体系视为有端部集中质量的悬臂梁，并建立对应的有限元模型，采用梁单元对塔架进行离散化建模[13]，通过有限元分析得到一阶模态质量和一阶自振频率并采用式（3）计算阻尼系数[3]：

caero=2M1ω1ξaero （3）

式中：caero为气动阻尼系数；M1，ω1分别为该有限元模型的一阶模态质量和一阶自振频率.

1.3 气动阻尼矩阵

在风向与风轮平面垂直的前提下，气动阻尼主要阻碍塔顶在风轮平面法向的振动，对其他自由度的影响不明显 [3，13].则浮式风机对应自由度上受到的气动阻尼力、力矩如式（4）和（5）所示：

Faero=-caerotop （4）

Maero=-caerotopL（5）

式中：Faero，Maero分别为气动阻尼力、力矩；L为Faero作用点到运动参考点的垂直距离，在风机俯仰角度不大的情况下，近似地认为L在运动过程中保持不变.

将整个浮式风机视作一个刚体.由于气动阻尼主要取决于叶片翼型的固有气动特性和风轮转速[14-15]，故这种情况下风机气动阻尼系数caero保持不变.当浮式平台发生运动并引起塔顶运动时，结合塔顶运动和浮式平台的运动关系，气动阻尼将阻碍平台的纵荡（surge）和纵摇（pitch）运动[6].利用式（6）所示的塔顶运动与浮式平台刚体运动的关系，将气动阻尼力、力矩写为式（7）和（8）：

2 水动力系数和系泊刚度

浮式平台为大尺度构件，宜采用三维线性势流理论计算水动力系数.根据拉普拉斯方程、海底边界条件、线性化自由液面的动力和运动边界条件可得线性入射势ΦI，在上述方程和边界条件的基础上补充物面边界条件，求解可得绕射势ΦD和辐射势ΦR.

再利用线性化的伯努利方程计算物体湿表面的压力分布，最终得到波浪作用下浮体上的波浪力、力矩.利用分析软件HydroD建立浮式风机质量模型和水动力模型，得到质量矩阵和包括附加质量、附加阻尼、波激力、静水回复力在内的水动力系数矩阵.

系泊系统提供的回复力与平台位移不成正比，即其回复刚度是非线性的，但可把系泊系统线性化并将其视为线性弹簧在运动方程中考虑其刚度[16].

3 频域内运动方程

得到气动阻尼、质量、水动力系数和系泊刚度矩阵后，在频域冉立式（11）所示运动方程：

-ω2M+Ma（ω）+iωCaero+Ca（ω）+

Khs+KmX（ω，β）=Fex（ω，β） （11）

式中：ω，β为入射波的频率和入射角；M为质量矩阵；Ma（ω） 为附加质量矩阵；Caero为气动阻尼矩阵，当不考虑气动阻尼时令Caero为零；Ca（ω） 为附加阻尼矩阵；Khs为静水回复刚度矩阵；Km为系泊刚度矩阵；X（ω，β）为浮式风机位移幅值矩阵；Fex（ω，β）为单位波幅规则波对应的波激力矩阵.

求解上述方程可得式（12）所示的浮式平台位移幅值关于波幅的传递函数H（ω，β），即RAOs（ω，β），其意义为特定频率单位波幅规则波作用下浮式风机的位移幅值，用于评价浮式风机的动态性能，如式（13）所示：

4 算 例

4.1 计算模型

选取NREL提出的5 MW风机和IEA（国际能源署）提出的OC3 HywindSpar浮式平台及对应的系泊系统作为算例.表1和表2分别列出了风机和浮式平台及系泊系统的主要参数[12，17].

对于OC3 HywindSpar平台，通过势流理论计算得到的附加阻尼不能反映平台受到的全部阻尼，还需在运动方程中添加与波浪频率无关的附加线性阻尼矩阵Ce[17].

Jonkman等利用FAST线性化功能，得到了该系泊系统的刚度矩阵Km，因FAST不能模拟该系泊中的三角连接，故另需在Km中添加附加艏摇刚度来考虑三角连接的贡献[17].

4.2 计算工况

代表性地选取作业和自存两种典型工况，并在各典型工况下细分了数组海况，如表3所示.选用JONSWAP海浪谱，取谱峰升高因子γ=3.3.考虑平台的对称性和为便于计算，取一个浪向角β=0°.

4.3 计算结果

根据前述章节，在时域模拟中排除了结构阻尼和附加阻尼的影响，得到了仅在气动阻尼作用下的塔顶位移衰减曲线，如图3所示（图中塔顶位移均从对应风速下塔顶振动的平衡位置算起）.代表性地给出了风速为7 m/s和30 m/s时的塔顶位移衰减曲线，作业工况其他风速下的衰减曲线与图3（a）相似，自存工况其他风速下的衰减曲线与（b）相似.基于该衰减曲线，并利用式（3）计算阻尼比，表4给出了不同风速下的气动阻尼比.

作业工况下的气动阻尼远大于自存工况下的气动阻尼，且与Tempel[18]给出的4%的估算值相当.作业工况下，风机在风速为7 m/s与11.4 m/s时，浆距角相同但风轮转速不同，风机在风速为11.4 m/s，16 m/s，20 m/s和24 m/s时，风轮转速相同但浆距角不同；而自存工况各风速下，风轮停止转动且叶片顺桨[3].结合上述结果与其他相关研究的结论[14-15]，易知风轮停止转动和叶片顺桨是导致自存工况下气动阻尼急剧减小的主要原因.

基于气动阻尼比建立不同工况下的气动阻尼矩阵，并结合水动力系数矩阵和系泊刚度矩阵，在频域内分别建立了考虑气动阻尼与否的运动方程，求解方程得到各工况下考虑与不考虑气动阻尼的RAOs，如图4所示.由于气动阻尼主要阻碍纵荡和纵摇运动[6]，故仅给出上述两个自由度的RAOs.主要频响范围在0.5 rad/s内，为便于识图仅绘出0.03～0.5 rad/s内的结果.

为分析气动阻尼影响的相对值，定义如式（15）所示的RAOs峰值相对差值.表5给出了不同风速下RAOs峰值的相对差值.

ζ=p0-paerop0（15）

式中：paero和p0分别为考虑与不考虑气动阻尼时RAOs的峰值；ζ为RAOs峰值的相对差值.

由图4可知，作业工况下气动阻尼能有效降低纵荡和纵摇运动RAOs的峰值，但对非峰值的影响很小.而自存工况下气动阻尼对RAOs的影响不大.再由表5可知，纵摇受气动阻尼的影响比纵荡更为明显，Roddier [7]的研究也给出了同样的结论.此外，分析得到的纵荡和纵摇运动固有频率分别约为0.05和0.213 rad/s，与

由图5可知，气动阻尼能在一定范围内降低该海况下纵荡和纵摇响应谱的幅值，且对响应谱中接近海浪谱峰频率的成分抑制效果更加明显.

表6和表7通过响应谱的零阶矩（即响应谱图形的面积）定量分析了气动阻尼对响应谱的影响.由表可知：作业工况下Tp大于5.5 s的海况中，响应谱零阶矩相对差值超过5%，气动阻尼对这种情况下的响应谱有一定影响，而对作业工况下Tp较小的海况和自存工况下的响应谱影响很小.

作业工况下，气动阻尼能大幅降低RAOs的峰值，但由于RAOs峰值频率和海浪谱峰频率相差较远，所以气动阻尼对响应谱的影响不如对RAOs的影响明显.作业工况下低海况中的Tp较小，海浪的能量在RAOs峰值频率范围内分布很少，故这种情况下气动阻尼对响应谱影响很小；而作业工况下高海况中的Tp较大，海浪的能量在低频范围内分布变多，故这种情况下气动阻尼能对响应谱产生一定影响.

[10]LEISHMAN J G， BEDDOES T S. A semiempirical model for dynamic stall[J]. Journal of the American Helicopter Society， 1989， 34（3）：3-17.

[11]JONKMAN J， BUTTERFIELD S， MUSIAL W， et al. Definition of a 5MW reference wind turbine for offshore system development[R]. Golden， Colorado： National Renewable Energy Laboratory， 2009： 5-16.

[12]IEC 61400-3 Wind turbines-Part 3： design requirements for offshore wind turbines[S]. Geneva： International Electrotechnical Commission， 2009： 21-32.

[13]刘雄， 李钢强， 陈严， 等. 水平轴风力机筒型塔架动态响应分析[J]. 太阳能学报， 2010， 31（4）： 412-417.

LIU Xiong， LI Gangqiang， CHEN Yan， et al. Dynamic response analysis of the tubular tower of horizontal axis wind turbines[J]. Acta Energiae Solaris Sinica， 2010， 31（4）：412-417. （In Chinese）

[14]FRERIS LEON L， FRERIS L L. Wind energy conversion systems[M]. London： Prentice Hall， 1990： 85-93.

[15]郭洪澈， 李强， 刘雄， 等. 气动阻尼对海上风力机筒形塔架的影响[J]. 太阳能学报， 2013， 34（8）： 1450-1457.

GUO Hongche， LI Gangqiang， LIU Xiong， et al. Influence of aerodynamic damping on tubular tower of offshore horizontal axis wind turbines[J]. Acta Energiae Solaris Sinica， 2013， 34（8）： 1450-1457. （In Chinese）

[16]李彬彬. 新型深吃水多立柱平台的水动力与运动响应研究[D]. 哈尔滨： 哈尔滨工业大学土木工程学院， 2011： 78-80.

LI Binbin. Investigation on hydrodynamics and motion performance of an innovative deep draft multispar platform[D]. Harbin： College of Civil Engineering， Harbin Institute of Technology， 2011： 78-80. （In Chinese）

[17]JONKMAN J M. Definition of the floating system for phase IV of OC3[R]. Golden， Colorado： National Renewable Energy Laboratory， 2010： 4-21.

[18]VAN DER TEMPEL J. Design of support structures for offshore wind turbines[D]. Delft： Delft Wind Energy Institute， Delft University of Technology， 2006： 107-132.