常 林，刘廷瑞

（山东科技大学 机械电子工程学院，山东 青岛 266590）

风能是一种具有较好发展前景的新型清洁能源[1]。风力机是捕获、转化风能的主要部分，传统风力机形式粗大，叶片刚度和强度较高，但是近年来随着材料科学和机械设计分析、制造水平的快速提升，风力机叶片的形式变得越来越细长，虽然风能捕获量得到了很大提高，随之而来的是刚度和强度较传统叶片降低了很多。叶片在高风速或达到失速条件时发生的动态不稳定现象称为失速颤振。叶片一旦发生颤振，带来的破坏是巨大的[2]。近年来发生在世界各地风场的风力机在大风速条件下的断裂失效也印证了此问题的广泛存在。2014年发生在湛江勇士风场的台风事故中[3]，风场33台风力机在台风“威马逊”的影响下几乎被全部吹倒，有5台风机完全损坏，其中叶片破坏是一种重要的破坏形式，且风力机叶片成本占风机总成本的15%～20%，所以叶片的保护至关重要。

风力机叶片所受载荷在气动力、弹性力、重力影响下发生耦合，极易在大风速条件下产生气动弹性不稳定现象。失速颤振作为气弹不稳定的一种重要形式，是叶片破坏的重要原因之一。叶片一旦发生失速颤振，便会表现出振幅不衰减的高频振动，而这种振动会在极短时间内使叶片失效。水平轴风力机叶片颤振分析涉及叶片气动弹性稳定性，即叶片在气动力作用下自身发生弹性变形，而这种变形反过来又影响周围流体，因此叶片上气动载荷在建模分析时应考虑其具有非线性[4]。此外，颤振分析还通常涉及弹性叶片的结构特征以及叶片附近非定常气流场的描述以及相互耦合的机理分析。所以叶片的颤振研究十分复杂。

叶片常用分析模型有弹性铰链模型、有限元模型、典型截面分析模型等。典型截面模型是一种在叶片气弹分析中经常采用的简化型叶片结构模型，典型截面模型可以简单、经济和快捷地获得失速颤振特性的主要结构和气动影响参数，是研究叶片挥舞/摆振气弹稳定性常用的一种有效结构模型，易于建模和模拟仿真。文献[5]基于定值转角下非线性气动力模型分析了挥舞-扭转以及挥舞-摆振的气弹不稳定问题，但过程未作线性化处理。文献[6]采用2D截面分析模型研究了挥舞-摆振耦合运动，考虑了结构阻尼以及耦合结果的稳定性。文献[7]将叶片简化为悬臂梁，考虑弯-弯耦合，在弹性梁建模的基础上引入材料阻尼从而得到运动方程。但在非线性气动载荷、大攻角下的失速颤振方面研究相对较少，且均未考虑风力机在实际工作过程中的塔影效应和风切变的影响。本文在课题组前期研究的基础上[8]，考虑风速的空间分布，主要基于风切变和塔影效应的影响，模拟非线性气动力的大攻角情形，进行叶片运动建模，提出MPC控制策略，抑制叶片非线性颤振。MPC控制策略与气弹模型结合的方法具有较高的可操作性且该方式易于实现，可以得到很好的控制效果；基于风切变和塔影效应的影响的研究更具现实意义，分析过程佐证了该方法不失一般性。

1 风速建模

风力机在实际工作过程中，所受风载并非是各处相同的，在方位角不断变化的情况下，叶片各点的作用风速也在变化，此处主要考虑两点影响：风切变效应和塔影效应。风切变是指风在运动过程中，因地面摩擦和热力因素的存在，使得作用在风力机上的风的速度随着高度的升高而下降的现象。由于风切变的存在，叶轮在运转过程中各处受载不均，甚至会出现叶片振动和功率波动，影响风力机的正常工作。图1为三叶片水平轴风力机示意图。

图1 风力机风速作用示意图

风切变在垂直高度方向的风速分布[9]，可以轮毂处中心风速为参考风速，得出某高度风速表达式

式中：V(H)是高度H处的风速，Vh为高度h处的参考风速，本文h取轮毂处高度。α为粗糙度指数，与地形相关。具体数值见表1。

表1 不同地形下粗糙度指数

在式（1）中以叶尖风速作为基准，考虑叶片的方位角θ以及叶轮半径R，可以将含H的表达式写为如下形式

风力机塔架对气流起到阻塞作用，改变风速的大小和方向，称为塔影效应。塔影效应主要存在于叶轮下方以塔架中心线为轴线的两侧120o的范围内，在塔架正前方叶片位置影响最大。建模时增加一个速度影响矢量因子以修正风速模型

其中：q=1+α(α-1)R2(8h2)-1，c0为风机塔筒支架半径，u0为分析截面微元距塔架中心轴线距离，称悬垂距离。

因此在综合考虑风切变和塔影效应之后，对叶片整体建立风速模型，以方位角和叶轮半径作为自变量可得风速模型如下

以小型风力机参数为输入参数，选取轮毂输入风速为15 m/s，根据《风力等级》国家标准，此风速对应疾风，对小型风力机破坏较大。塔架半径取为0.5 m，塔架高度取为6 m，风轮半径取为2 m，粗糙度指数因子选取为0.2，悬垂距离取为0.5 m进行仿真，由风速模型经仿真计算可得风速空间分布函数。将风速函数离散化为N位，取平均值作为叶片气弹模型的输入风速。

2 叶片结构模型及气弹方程

近年来该领域学者多以大型风力机作为研究对象，而对小型风机研究相对较少。本文选取小型风机结构模型数值进行仿真分析。考虑分析截面为叶片翼型典型对称面，叶片气动弹性中心与结构重心相距较近，因此弹性扭转变形可忽略。分析翼型截面如图3所示，其中U为输入风速，此处选为经风速模型处理后的平均风速；c为翼型截面弦长；v0为相对风速；α0为攻角；ϕ为变桨角；ψ为相对风速角；翼型上作用的非线性气动力包括非线性循环和非循环气动升力TNc与TC1；DN1为非线性中的气动阻力项。

图2 翼型截面参数示意图

如图2所示的运动中，y方向表示叶片垂直于叶轮旋转平面的挥舞运动；z0方向表示展长方向；z方向表示叶片在旋转平面上垂直于挥舞的摆振运动；x方向表示叶片转动。建模计算时综合考虑系统动能和势能，基于前期研究成果[10]，利用拉氏方程展开，忽略高阶项，可以得到挥舞和摆振运动方向运动方程

根据前期研究成果，采用大攻角下修正后的ONERA气动模型及非线性求解理论，利用傅氏变换提取气动力中的谐波项，可以得到气动力描述模型升力方程为

气动阻力方程为

其中：KF1、KF2、KDF为升力和阻尼分别对应的循环结构变量，比例系数m=0.5，j=0.014，循环变量方程表达为

其中：v0是相对风速，与输入风速和转速相关。为翼型截面回转半径，b为翼型截面半弦长，ρ为空气密度，SL、qL、ρL、LK分别为气动项系数。ΔLKC为静态气动升力曲线延长线值与非线性升力曲线的差值。与之相对地，ΔLKD为阻力曲线差值。其余气动参量，根据项目组前期研究成果，参考文献[11]，取值如下

式中其余参数取值如下

3 气弹方程求解仿真

上述气弹方程的仿真需通过降阶后化为系数矩阵状态空间形式，将2阶方程化为1阶，使方程易于仿真。同时在输出方程中取相应矩阵，通过观察输出值对应振动位移的振幅和频率来判断叶片相应方向的振动情形，从而判断系统响应是否发生发散，颤振是否发生。本文取常值变桨角进行仿真运算，观察风载作用下叶片挥舞-摆振方向的振动情形。

对气弹方程组进行降阶处理。定义Y=[XTX′T]T，将方程联立 ，可以得到 ：可右乘输入、控制矩阵，在时域内表示为u(t)。其中

系统稳定性的判定方法有许多种，诸如特征值判定法、伯德图、内奎斯特图、李雅普诺夫判定法等。在实际应用中，应综合分析系统响应状况。本文所述叶片系统，应从叶片响应振幅大小、稳定时间、静差等多方面分析，且往往出现某一个方向稳定而另外方向不稳定的现象，所以应观察各方向振动情形，且该种稳定情形需符合物理和工作实际，否则稳定性判别无现实意义。

使用4阶5级龙格-库塔方法进行计算仿真，模拟大攻角大风速条件下的系统响应。取恒叶尖速比系数为1.2，λ=LΩ/U，L为叶片长度，取L=0.9R，变桨角ϕ=π/6，由以上给定条件可计算出攻角数值。由于风速和转速的共同影响，攻角值或在大风速条件下表现出某一较大值，从而影响叶片气弹振动。

回转半径取在叶根处，模拟危险截面振动。一旦叶片发生颤振，就会表现出振幅不自衰减的振动，不断吸收能量维持发散振动。叶片发生颤振时极易破坏失效，且颤振的发生很难察觉，除已有报道的颤振发生前的“蜂鸣声”外几乎无征兆可循。所以颤振的抑制是风电技术发展的重要方面。模型结构参数选取匹配于小型风力机的结构数值。

在确定模型输入风速的前提下，代入小型风力机结构数值模拟大攻角状态下系统响应，通过4阶5级龙格-库塔法计算出风力机叶片振动位移。通过图3可以明显看出，在大攻角大风速条件下，挥舞和摆振两个方向都出现了振幅不自衰减的高频振动，也由此可判断出叶片此条件下发生颤振，需加以抑制。

表2 结构参数取值

4 振动控制策略

模型预测控制（MPC）是上世纪80年代发展起来的一种计算机控制算法。该种控制策略一经出现便运用在实际工作过程中，并且在工作实践中不断完善和发展。该种算法采用步长可给定的滚动优化和反馈校正等控制策略。MPC的控制思路是：根据模型未来的预测行为，基于误差校正和滚动优化选择对受控目标最合适的控制动作进行控制[12]。选取叶片气弹状态空间模型为受控模型，将该模型离散化处理可得

其中两式式尾新增两项为叶片气弹模型误差修正项。令ΔY(k)=Y(k)-Y(k-1)，

所以可以将式（14）化为增量形式

输出可写为增量形ΔY0(k+1)=Y0(k+1)-Y0(k)，所以式（15）可以写为

图3 未控制时挥舞、摆振方向振动情形

式中W1=010×10，W2=I10×10。定义新状态变量y(k)=[ΔYT(k)YT(K)]T，所以式（16）可以写为

式中各系数矩阵均可由式（16）推算出，此处不赘述。基于上式所列系统状态空间模型，可以得到系统在离散化后的某时刻系统轨迹。按照引入的系统模型误差以及预设目标值进行离散化的预测时域内的滚动优化和误差控制，按照MPC控制策略的一般规则，选取预测时域长度为11，选取控制时域长度为10，采样步长和采样次数分别为0.5、100，选取无约束条件下输入、输出，基于系统状态空间的模型预测控制策略运用MATLAB软件进行仿真，可以得出叶片控制后振动位移。

通过图4可以明显看出，经MPC策略控制后，挥舞和摆振方向振动最终稳定，很好抑制了叶片颤振问题。且挥舞和摆振方向的稳定位移减少为控制前的3%、5.83%，稳定时间少，静差可接受。

在实际工程中，有多种方法来实现该控制策略：运用贝加莱PLC可以实现从MATLAB/SIMULINK代码模型直接转化为PLC控制程序；运用三菱PLC可以实现从MATLAB代码-C语言-PLC语言的转换；运用西门子TIA技术可以直接实现从MATLAB语言到PLC程序的转换[13]。或通过对控制器编程实现MPC控制算法与传感系统和执行装置的结合，组成闭环振动控制系统。执行装置可选择变桨激励器、蒙皮内置作动器（如形状记忆合金电加热、压电材料电流控制等）进行物理实现，限于篇幅不再赘述。

5 结语

本文基于实际应用意义提出了可以实践于叶片气弹稳定性控制的控制策略和系统模型。经过仿真对比，MPC控制策略可以很好抑制叶片的气弹颤振：叶片两个方向上的位移发散得到根本抑制；振动的振幅、频率大大减小，且具有较短的稳定时间和极小的静差。并且该方法对受控模型精确度要求不高，具有控制效果好、鲁棒性强等特点，在某些行业已取得了一定的发展，但在风力机叶片振动抑制方面应用较少，本文在叶片气弹振动抑制技术和实践控制理论相结合方面做出了一点探索。

图4 经MPC控制后挥舞、摆振方向振动情形