兆瓦级风机停机状态塔筒风致响应均值半经验计算方法

2021-09-27 07:05:28黄国庆

振动与冲击 2021年16期

黄国庆，袁酊，刘敏，许楠

(1.重庆大学土木工程学院，重庆 400045；2.中国广核新能源控股有限公司，北京 100000)

风能是一种可再生、无污染的绿色能源，储量丰富。我国陆地风力发电储量合计10 亿kW[1]。风电行业发展迅速，截至2019年底，我国陆上风电累计装机容量位列全球排名首位，累计陆上装机容量占全球陆上风电装机容量的37%。为提高风机的经济效益，单机容量变得越来越大，目前正向5 MW以上发展[2]。

当前风机的选型主要按照国际电工委员会发布的标准IEC 61400-1[3]进行。根据风机安装场地的风速和湍流度将场地分为Ⅰ级、Ⅱ级、Ⅲ级和S级以便选择风机叶片型号。但这种场地分级方法只能大致将实际场地的风速和湍流参数与上述标准分类匹配，而实际风场与规范规定风场存在的差异无法考虑。因此在精细化设计，特别是塔筒设计中还需要大量分析。在设计塔筒时，可以用数值仿真软件进行估计来获得比较准确的响应，但是数值仿真需要进行耗时的时程分析才能确定动力响应的统计特性，计算效率较低，难以适应风机快速、精细选型的需求。因此急需更加高效的风机塔筒响应计算方法。

目前设计规范中常采用等效静风荷载将复杂的动力分析问题转换为静力问题来评估风致结构响应(如构件的弯矩、轴力和变形等)。借鉴这一思想，Tempel[4]将风力机系统等效为一个单自由度体系，通过理论方法计算获得简化的风振动力响应。贺广零[5]提出了两自由度体系阵风荷载因子，其安全性较单自由度体系阵风荷载因子更高。柯世堂等[6]通过建立一体化有限元模型，揭示了风力发电塔轮系统风振响应背景与共振部分的分布特性和模态间耦合效应的作用机理。Xu等[7]推导了风机停机状态下顺风向和横风向基底弯矩及其极值的解析表达式。该解析表达式由均值和脉动两部分组成。在Xu等的研究中，基底弯矩极值的解析表达式是在准定常假设的基础上推导的，需要将叶片上各个截面的风荷载在叶片长度方向进行积分，形式复杂，不便于直接应用于设计。为了便于工程应用，Xu等[8]对多个不同型号的风机进行数值仿真后，通过参数拟合获得了简化的基底弯矩极值解析表达式，并将简化公式计算的响应与数值仿真得到的响应进行了对比[9]，结果表明在400 kW的小型风机上该简化公式的适用性较好。该简化公式可快速简便地获得响应预估结果，被日本《风力发电设备塔架结构设计指南及解说》[10]所采纳。

基底弯矩的均值在基底弯矩极值中占比可达50%，所以准确地估计其均值是计算基底弯矩极值的基础。需要指出的是在Xu等的研究中，计算基底弯矩均值时不考虑结构变形对基底弯矩的影响。该假设对于小型风机是合理的，但对于大型兆瓦级风机，由于叶片长度和塔筒高度都显著增加，结构整体刚度明显降低，叶片变形和塔筒变形相比小型风机明显增加，导致更大的基底弯矩等响应。因此基于柔性塔筒和刚性叶片假设的均值解析表达式对5 MW及以上的大型风机的适用性还有待进一步验证。

限于篇幅，本文以塔筒基底弯矩的均值为研究对象。针对风机的不利工况之一的停机状态，本文研究了塔筒响应均值计算解析表达式。首先基于以往学者的工作重新推导了风机基底弯矩均值解析表达式，然后以5 MW的变桨距控制风机为例，将均值解析表达式的计算结果与OpenFAST软件的计算结果进行比较，探讨了误差产生的原因。最后提出了考虑叶片和塔筒变形的基底弯矩均值半解析方法，以期减小解析公式计算误差，为风机塔筒响应均值的快速计算提供更精确的方法。

1 基底弯矩均值解析表达式

根据日本《风力发电设备塔架结构设计指南及解说》，对于变桨距控制风力机，当停机状态下作用于塔筒的最大风荷载发生在偏航角为90°，桨距角为90°，叶片1方位角为90°时，如图1所示。

图1 顺桨状态下变桨距控制风力机的姿势Fig.1 Parking state of pitch control wind turbine

来流风速包括水平方向分量平均风速U0和脉动风速u0以及竖向、横向的脉动风速w0和v0，如图2所示。由于各翼型的升阻力系数是基于垂直于叶片长度方向的截面定义的，因此计算叶片上作用的风荷载时需要将来流风速分解到沿叶片长度方向和垂直叶片长度方向。横向脉动风速v0始终垂直于叶片长度方向，因此v0不需要分解。将顺风向风速U0+u0与竖向脉动风速w0的矢量组合速度分解为沿叶片长度的分量W1和垂直叶片长度方向的分量V1。在图2中，U0+u0,w0,W1和V1处于同一平面内，其中分量W1引起沿叶片长度方向的荷载几乎为0，V1由均值部分U和脉动部分u组成，即V1=U+u。V1的均值部分U由平均风速U0的分解得到，如图3所示。叶片2相对于风向是倾斜的，叶片上平行于横截面的平均风速为U=U0cosβ，其中，β为叶片的方位角。V1的脉动部分u由脉动风速u0和w0分解得到，u=u0cosβ-w0sinβ。叶片3上风速分解的方法与叶片2的分解方法类似，不予赘述。

图2 来流风速分解示意图Fig.2 Wind speed decomposition

图3 平均风速分解示意图Fig.3 Mean wind speed decomposition

图4 叶片2截面所受升阻力示意图Fig.4 Lift force and drag force of a section of blade 2

(1)

(2)

式中：ρ为空气密度；CD(θ+θ′,r)和CL(θ+θ′,r)分别为在距离叶根r处的某截面在攻角为θ+θ′时的阻力系数和升力系数；c(r)为在距离叶根r处叶片的弦长。

(3)

将D和L分解到x,y方向，并在θ处泰勒展开，可得阻力fD为

fD(r)=D(θ+θ′,r)cosθ′-L(θ+θ′,r)sinθ′≈

(4)

相似的，升力fL可以表示为

(5)

(6)

(7)

(8)

(9)

每个叶片阻力和升力合力的均值分别为

(10)

(11)

式中，i=1,2,3。

计算塔筒基底弯矩时，先将分布在叶片上的风荷载简化为作用在轮毂高度处的集中荷载与附加弯矩，轮毂处集中荷载可以表示为

(12)

(13)

(14)

(15)

在实际情况中，塔筒底部的弯矩还包括由于风机质量的不对称分布导致的偏心弯矩，但本文中所述基底弯矩只考虑由等效到轮毂处的叶片总风荷载和塔筒上风荷载导致的弯矩。塔筒底部由顺风向风荷载和横风向风荷载产生的弯矩均值分别为

(16)

(17)

式中：H为轮毂高度；d(z)为塔筒直径；CD(z)为塔筒截面阻力系数；CL(z)为塔筒截面升力系数；U(z)为沿高度变化的平均风速。

2 案例分析

以美国国家可再生能源实验室(NREL)5 MW陆上风机为例，将均值解析表达式的计算结果与OpenFAST的计算结果进行比较，分析均值解析表达式产生误差的原因。NREL 5 MW陆上风机主要参数，如表1所示。

表1 5 MW风机计算中所用主要参数Tab.1 Some parameters of 5 MW wind turbine

5 MW风机的每个叶片有19个控制截面，对应8个不同的翼型，各控制截面与翼型的对应关系，如表2所示。其中部分翼型的升力系数和阻力系数随攻角的变化，如图5所示。

表2 5 MW风机叶片截面和翼型信息表Tab.2 Blade sections and airfoil information for 5 MW wind turbine

图5 风机叶片部分截面升力系数和阻力系数随攻角变化Fig.5 Lift coefficient and drag coefficient for different airfoils

2.1 风场

平均风速剖面采用指数率，U(z)=UH(z/H)α1，H为轮毂高度，UH为轮毂高度处平均风速，本案例中设置为30 m/s；α1为地面粗糙度指数。湍流度设置为0.18，不考虑湍流度沿高度的变化。脉动风速u0,v0,w03个方向的自功率谱均采用Kaimal谱，由式(18)给出

(18)

式中:k=u0，v0,w0；f为频率；Lk为各脉动风速成分的积分尺度，其取值由IEC规范给出

(19)

式中，ΛU为湍流积分尺度，其定义为

ΛU=0.7·min(60 m,H)

(20)

式中，min(60 m,H)为取轮毂高度H与60 m中较小者。标准差之间的关系定义为

σv0=0.8σu0,σw0=0.5σu0

(21)

相干函数cohi,j采用风力发电机设计要求(IEC规范)中u0方向风速分量的相干函数

(22)

式中：d为两点间的向量在与平均风向垂直的平面上的投影线的长度；a为衰减系数；Lc为相干尺度参数。在IEC规范中，参数a和Lc的取值为

a=12,

Lc=5.67×min(60 m,H)

(23)

IEC规范中没有指定横风向风速分量和竖向风速分量的相干函数。本文采用OpenFAST进行风场模拟时采用的横风向风速分量和竖向风速分量的相干函数[12]

(24)

2.2 攻角检验

从图5中可发现升力系数和阻力系数均对攻角敏感，对攻角的校核是风荷载计算的关键一步，故有必要对比OpenFAST和解析表达式计算的攻角均值。图6给出了叶片2由两种方法计算的叶片19个控制截面的平均攻角对比情况。可以看到两种方法得到的各截面攻角吻合良好，说明了采用解析表达式计算的攻角是正确的。

图6 OpenFAST和解析表达式计算的叶片截面的攻角θ均值对比Fig.6 Comparison of mean value of blade attack angle calculated by OpenFAST and analytic method

2.3 计算结果及误差分析

该停机工况在OpenFAST中运行一条600 s时程需要约10 min，为得到准确的统计量需要计算大量时程。使用解析表达式编程后大幅提升效率。

为了探究风机结构刚度对基底弯矩均值计算的影响，通过调整叶片和塔筒的自由度，分别设置了如表3所示的4个计算工况。在使用OpenFAST计算时，通过TurbSim模块模拟20个风场时程，每个风速时程长度为630 s，然后计算风机响应时程。计算基底弯矩均值时，为了避免瞬态效应，前30 s未使用，将这20个时程的均值再取平均后作为该工况的响应统计量，计算结果如表3所示。

表3 风机在不同刚柔性设置基底弯矩误差结果Tab.3 Error when different stiffness and flexibility

从表3中可得出以下结论。①当塔筒设置为刚性时，比较工况1和工况2，叶片的刚柔性对风荷载引起的基底弯矩的影响不大。②与工况1刚性塔筒刚性叶片相比，工况3和工况4中柔性塔筒工况下，由于塔筒顶部在风作用下产生位移，导致实际的偏心距与刚性塔筒工况中的初始偏心距略有差别，从而产生新的偏心弯矩值。而实际情况中的基底总弯矩中包含了风力弯矩、偏心弯矩和附加弯矩三部分(OpenFAST获得的弯矩包含各种荷载作用的响应)。在从中得到风力弯矩时，由于偏心弯矩发生变化，可能导致风力弯矩计算有所偏差。③在保持叶片为刚性的前提下，塔筒的刚度变化对结果有一定的影响，见工况1和工况3；计算中发现在刚性塔筒和柔性塔筒工况中，叶片上升力和阻力的合力相差不大(见表4所示工况1和工况3的结果对比)，说明塔筒变形导致的力臂变化可能是导致基底弯矩差别的主要原因。④当塔筒设置为柔性时，如工况3和工况4所示，柔性叶片工况产生了比刚性叶片工况更大的弯矩，可能的原因是：风机整体结构越柔，变形越大，力臂的变化也越大(图7显示叶片2各截面沿U方向的位移均值)。

表4 3个叶片上风荷载合力均值误差表Tab.4 Error of mean value of wind load on 3 blades

图7 叶片2各截面沿U方向的位移均值Fig.7 Mean displacement of blade 2 along the U direction

3 基底弯矩均值修正系数

从表3中的计算结果发现，现有的解析表达式计算得到的基底弯矩均值均小于OpenFAST计算得到的均值，所以按照现有规范设计偏危险，需要进行基底弯矩的修正。

根据IEC规范，轮毂高度处湍流度的变化范围约为0.11～0.18。柔性塔筒柔性叶片工况下，轮毂高度处平均风速UH=30 m/s，风剖面指数α1=0.2，湍流度分别为0.10和0.18时顺风向和横风向基底弯矩的均值，如表5所示。

表5 风机在不同湍流度下基底弯矩均值的相对差值Tab.5 Error of response mean under different turbulence

基底弯矩的误差与轮毂高度处平均风速、风剖面指数和湍流度相关，但从表5中发现，在不同湍流度下基底弯矩的差别不大。因此，基底弯矩均值的修正系数为轮毂处平均风速UH和风剖面指数α1的函数。

IEC规范将地面粗糙度分成了4个等级，分别对应的风剖面指数为0.10，0.15，0.20，0.27。由于在平均风速大于等于25 m/s时风机进入停机状态，同时在IEC规范中，以轮毂高度处的平均风速将风机的建设等级分为了3类，其中Ⅰ类风机的参考平均风速为50 m/s，因此将风速的计算范围确定为25～50 m/s。本文计算了在α1分别为0.10，0.15，0.20和0.27时的不同地面粗糙度场地上，轮毂高度处湍流度为0.18，轮毂处平均风速从25 m/s变化至50 m/s时，解析表达式与OpenFAST的顺风向和横风向基底风力弯矩均值误差。

定义顺风向和横风向基底弯矩均值修正系数分别为ΔD和ΔL，表示为

(25)

在不同的平均风速下顺风向和横风向基底弯矩修正系数，如图8和图9所示。