抑制风力发电机叶片高频振动的尾缘襟翼主动控制方法

2021-10-17刘廷瑞孙长乐李善耀张晓林刘桂芳

工程科学与技术 2021年5期

刘廷瑞，孙长乐，李善耀，张晓林，刘桂芳

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

长期以来，学者们认为基于失速颤振和经典颤振的挥舞/摆振/扭转耦合位移的发散不稳定运动是回转叶片（如风力机和直升机叶片）结构断裂失效的重要原因，其中挥舞/扭转（弯扭）运动造成的失效是其中表现之一[1–2]。而近年来，尾缘襟翼结构被广泛用于叶片气弹减载和颤振控制。Haselbach[3]提出了一种基于壳体单元模型的先进叶片建模方法，其中尾缘区域中的粘结线通过实心砖元件离散化，这可以提高风轮机叶片的结构预测的可靠性。Llorente等[4]提出了一种新的实验方法，重点研究了具有尾缘锯齿的风力机叶片的空气动力性能的变化，以及具体的设计过程和预测方法。Chen等[5]设计了一种尾缘流量控制装置并讨论了其性能，包括通过求解旋转框架内的3维雷诺平均方程，在风轮上安装微插片并在尾缘设计微射流单元，以达到减载目的。张文广等[6]以NREL–5MW风力机为研究对象，实现了智能叶片风力机建模及多目标尾缘襟翼控制，研究了尾缘襟翼在风力机主动降载和功率控制方面的效果。张广兴[7]将风力机叶片描述为展向分布尾缘襟翼的复合材料悬臂梁叠层板，通过瑞利–里茨法实现了结构建模，建立了气弹模型，实现了振动频率分析和主动振动控制。穆安乐等[8]将风力机叶片简化为展向分布尾缘襟翼的复合材料悬臂梁叠层板，通过瑞利–里茨法实现了结构建模，并结合Theodorsen气动力，建立了气弹模型，最后采用模型预测控制算法实现了主动振动控制。

然而很少有学者注意到造成叶片断裂失效的隐性故障：在常规工况下，一种由于长期处于“低幅–高频率”振动[9–10]条件下的叶片，会形成隐性的裂纹缺陷，该缺陷使得叶片体在极端颤振条件下，更容易造成断裂失效。该高频振动在湍流时的尾缘处容易产生[10]，或者在柔性叶片的气动力矩表现为低幅、高频振荡[11]时的失速状态下容易产生。Zuheir等[12]也基于有限元分析，证实了尾缘处存在的低幅、高频振动效应。

作者项目组成员在前期工作中发现：高频信号的可能频率范围确定为复合材料层合结构的2阶至5阶固有频率之间。其一，当外因激励频率小于复合材料叶片的1阶固有频率时，不会发生持续的等幅高频振动（1阶频率时的共振）现象，而当外因激励频率等于复合材料叶片的1阶固有频率时，可以发生共振，但不会出现微幅共振现象；其二，当外因激励频率大于叶片的6阶固有频率时，将不会有共振现象发生，而可能发生发散不稳定现象，同时在解耦算法Galerkin法[12]的应用中，对于风力机复合材料叶片而言，Galerkin法的保留振型在5阶范围内即可以充分满足精度要求。

由于振动的高频特性，诸多智能控制算法由于迭代次数太大或者寻优难度问题而失效；一般的线性控制方法，也仅仅能进一步实现幅度的衰减，在抑制高频方面效果较弱。作者基于一种沿叶片展长分布的尾缘襟翼结构，经由步进电机驱动，并通过基于LMI设计的H∞控制方法，实现了“低幅–高频”振动控制。其实质是：通过加入步进电机硬件，并结合鲁棒控制方法，改变了整个气弹系统的频率结构，在相同的外界激励情况下，使得新系统能够避开原来气弹系统的2阶至5阶固有频率范围内的共振。

1 气弹系统方程

为简化复合材料结构参数的计算，叶片采用固定弦长的实心复合材料叶片梁，翼型采用二次拟合的NACA64418翼型[13]，如图1所示。

图1 叶片梁结构及翼型运动参数Fig. 1 Structure of blade beam and airfoil motion parameters

图1中：L为叶片展长；r为翼型截面到旋转中心的位移；y、z分别为摆振和挥舞方向的运动（位移），并假设x为展长方向；θ为扭转位移；α为攻角；U为风速；c为叶片弦长；β为襟翼摆动角度，L0为襟翼2D截面（沿弦向）长度；V0为入流风速；λ0为速比系数；Ω为常值转速。襟翼采用轻质太空铝材料，襟翼两侧与叶片梁母体采用松配合铰接方式连接。襟翼通过平头键由步进电机带动旋转，空心襟翼的中段由叶片梁通过轻质凸轮摇摆结构通过光滑接触而驱动，起到平稳传动的作用（所有轻质结构质量不计，步进电机固定在实心叶片中部靠近尾缘的位置）。

复合材料铺层，以弦向为对称，采用基于周向反对称刚度（CAS）的铺层设计。在CAS铺层下，叶片展现单一的弯扭耦合位移，即挥舞/扭转耦合位移[14]。则单一截面的挥舞z/扭转θ 耦合运动可描述为：

式中，GJ、EI、Kc分别为相应的扭转刚度、弯曲刚度、弯扭耦合刚度，mc为单位长度翼型质量，KA、Km为转动惯量相关系数，M、F分别为相应的准稳态气动扭矩和气动升力。

文献[15]展示了一种悬臂振动2D翼型的基于尾缘襟翼的准稳态气动力。作者将文献[15]中襟翼激励气动方程中的绝对风速改为风力机的相对风速，将其中的变桨角度替换为风力机叶片的攻角α，并增加相对转动带来的襟翼激励气动项，在襟翼2D截面（沿弦向）特定长度为L0=c/6的基础上，得到适合于襟翼结构的旋转翼型的准稳态气动力表达式，可以描述为：

式中：ρa为空气密度；b=c/2；C1α、Cmα、C1β、Cmβ分别为相应的弦向襟翼匹配系数，并且，满足C1α=6.28，Cmα=(0.5+c/6)C1α，C1β=3.358，Cmβ=–0.635。

2 方程求解及有效性检验

Galerkin法常用来离散化并求解系统（1）的偏微分方程组[12]，弯扭耦合位移分别描述为：

式中，zi(x)、θi(x)为振型函数。本设计采用一种满足悬臂梁边界问题的振型函数以简化分析[16]：

式中，M0、C0、K0分别为 2N×2N质量、阻尼、刚度矩阵，Q0为 2N×1数值矩阵。

将式（8）转换为1阶系统，得到：

式中，状态变量Y=[X]T，且满足：

2.1 Galerkin法有效性检验

振型函数zi(x)、θi(x)的精确性对于Galerkin法求解的有效性至关重要，故本设计的Galerkin法有效性需要进一步检验。经过仿真可知，本设计中振型保留阶数N=4时，求解精度完全满足要求，由式（9）中矩阵A0的特征值λ可以求得挥舞/扭转振动的前4阶频率。令 β=0,U=0，则式（1）～（2）退化为悬臂叶片梁的自由振动方程，此时求解对应式（9）中矩阵A0的特征值λ可以获得挥舞、扭转自由振动各自的4阶频率。在CAS设计的铺层角度为60°时，复合材料结构参数数值计算如下：GJ=3.3935×107N·m2，Kc=−4.8869×105N·m2，EI=1.8822×107N·m6,KA=1.418×10−1，mc=1.703×102kg，Km=1.239×10−1，Km1=3.023×10−2，Km2=1.201×10−1。

叶片结构参数和运动参数如下：L=5 m，U=15 m/s，c=0.6 m， Ω=λ0U/L= 1 ×U/L。文献[17]展示了在CAS设计的铺层角度为60°时，系统呈现以挥舞弯曲为主模态的振动，本设计恰恰体现了这一特征。经Galerkin法（GM）求解，得到主模态前4阶频率f及其相应的特征值λ如表1所示。

另外，Stefan等[17]也提出了一种精确计算悬臂结构自由振动主模态频率f的计算方法（CM），描述为：

式中，参数g的两种形式为：g=jk1，k1=(2m+1)π/2/L，m=0,1,2，3；g=k2，k2为 c os(k2L)cosh(k2L)的解。经式（10）精确计算的挥舞运动的前4阶固有频率f，CM法计算的结果如表1所示。

表1 两种方法前4阶频率f 比较Tab. 1 Comparisons of the first four-order frequencies f based on two approaches

从表1的结果对比可以看出：在60°铺层角时，两种方法的前4阶频率结果相近，特别是高阶次的2～4阶频率结果非常接近，高阶频率的误差不超过0.68%。同样在铺层角度（0°～90°）范围内，每隔15°计算了复合材料结构参数，分析了前3阶自由振动的频率，并将GM与CM法进行了对比，如图2所示。由图2可以看出：1阶频率均比较接近，而高阶频率的误差不超过0.7%。GM法的振型函数选择的有效性及保留振型项数N=4的有效性，为高阶振动的颤振抑制提供了保证。

图2 前3阶自然频率对比Fig. 2 Comparisons of the first three-order natural frequencies

2.2 基于工程应用的方程改进

图1中的轻质结构质量可以忽略，但步进电机质量和动态驱动效应不能忽略。假设步进电机质量为ms=2kg，将其视为一种弹性旋转体的payload载荷[18]，则能对耦合运动的振动起到抑制作用，但对于控制输入而言，它亦相当于一个扰动信号，其基于仿真时间的扰动模型可以近似表达为ω=2ms(L/2)·sin(2πf0t)/(t+1)。其中：t为仿真时间[17]；f0=5/(2π)为步进电机驱动角度β变化的脉冲频率，该频率以一种正弦信号效应来影响控制系统的输入[19]。故考虑步进电机的实际驱动效应，系统方程（9）可以改进为：

式中：C1和D1为H∞控制设定的系数矩阵；B1=[01×2N11×2N]T,其相对于位移部分的分量为0，而相对于速度部分的分量为1。这是由于全反馈过程中，状态变量的测量是利用加速度传感器测量挥舞/扭转加速度而得到的[20]。由测得的加速度信号；经过滤波滤去低频扰动后，积分操作得到速度信号；然后第二次滤波后，再次积分操作得到位移信号。所以，来自输入的扰动信号首先影响加速度状态变量，其次影响速度状态变量。

改进后的系统方程（11），无法通过常规的智能控制方法来求解，本设计采用一种H∞控制理论可实现优化控制。

3 H∞控制及结果分析

对于系统方程（11），给定鲁棒性能系数 γ＞0，存在矩阵P1=＞0以及P2，根据Schur余引理，如果满足线性矩阵不等式（LMI）[18]：

则H∞控制的状态反馈控制器可以设计为：

式中，K=[k1k2···k4N]。

控制目标设计为：对于任意扰动 ω＞0闭环系统，鲁棒性可以表示为：式中：鲁棒控制参数ρ与控制输入信号联系，是影响控制输入的因素；γ是待优化的鲁棒性能参数。

在控制输出表达式Z=C1Y+D1β 中，C1=[I4N×4N04N×1]T，D1=[01×4Nρ]T，则有：

因此，

所以，使得闭环系统内部稳定的且满足式（16）的H∞控制器则正好是满足控制目标设计（14）的控制器K。

LMI不等式（12）的求解是一个复杂的问题，这是一个线性矩阵不等式约束的最优化问题，可以采用MATLAB的YALMIP工具箱简化求解过程，在选定鲁棒控制参数ρ的基础上，优化不确定的鲁棒性能参数γ，使得被控位移和控制输入β保持在一个合理的范围。以上为基于LMI的H∞控制的完整设计流程。

为降低全状态反馈时状态变量检测的误差影响，提出了一种利用状态重构和状态观察器来改善系统性能的思路。由于本挥舞/扭转耦合系统的能观性矩阵始终满秩，故可以设置系统观察器。基于状态观察器的状态重构方案如图3所示。

图3 基于状态观察器的状态重构方案Fig. 3 State reconstructing scheme based on state observer

图3中，观察器增益为G。选择G矩阵，使A0−GC1的所有特征值具有负实部, 采用一种固定的期望特征值的选择范式：先求取矩阵A0的特征值，实部为负值则保留，实部为正值则直接将其实部取反之后保留；将所得的全部特征值作为期望特征值，构成特征多项式；将该多项式作为状态观察器的特征多项式。由此，（16）可以修改为：

从而使得γ更具有鲁棒性，且在不同的鲁棒控制参数条件下，均可以得到合理的鲁棒性能参数γ。

3.1 理论控制与实际控制结果比较

图4为无控制情况下的挥舞/扭转位移。由图4可以看出两个位移均处“低幅–高频”振动状态。虽然对于扭转位移而言，振动频率更高，但其幅度数量级非常小（1 0−18），以致现有的任何工业化传感器都无法测出其数值，完全可以忽略，所以不会形成断裂失效的潜在故障；对于挥舞位移而言，其振动幅度在高频情况下可造成断裂失效的潜在故障，这也正是本文控制“高频振动”的目的所在，在后续主动控制中，将主要讨论挥舞振动的控制效果。另外，由于扭转位移非常小，所以不会影响沿展长分布的尾缘襟翼机械结构的扭转运动，为尾缘襟翼的摆动控制实现提供了条件。

图4 无控制的挥舞/扭转高频振动位移Fig. 4 Flap-wise/twist displacements of high frequency vibration without control

针对图4中的高频振动案例，取鲁棒控制参数ρ=1，当 ω =0（即无步进电机驱动的payload载荷效应和步进驱动效应）时，实现H∞理论控制，将其结果与 ω≠0时的H∞实际控制结果对比，如图5所示。图5（a）为H∞控制下的被控挥舞/扭转位移对比结果，图5（b）为控制性能对比，包括襟翼摆角β和鲁棒性能参数γ波动的对比。从图5（a）中的被控位移1对比结果可以看出：与图4中的未控制案例相比，理论控制和实际控制不仅在幅度抑制上取得了明显的效果，而且在高频振动控制上取得了更为显著的效果。以扭转位移为例，未控制时的扭转位移在10 s后仍然维持高频振动，而理论控制的扭转位移在1 s后便退出高频振动，更显著的是理论控制和实际控制的挥舞位移一开始便直接进入了低频振动状态。从理论控制效果和实际控制效果对比结果可以看出，扭转位移的实际控制效果要优越于理论控制的效果，而挥舞位移的实际控制效果要次于理论控制的效果（但已达到控制目标）。这可能是因为：步进电机驱动下，摆角β的控制能直接作用和影响扭转位移θ，而间接影响挥舞位移z的控制。需要说明的是：针对挥舞控制而言，与图4相比，其控制前后的频率数值均可直观展示，需专门分析频谱，这也是H∞控制的一个优点。

图5 H∞ 控制下的位移响应以及控制性能Fig. 5 Displacement responses and control performance under H∞ controller

图5（b）为控制性能对比，包括控制输入信号和鲁棒性能两方面的对比。从控制器输入（襟翼摆角β）的波动可以看出：理论控制的摆角范围超过了±5rad的范围即±287°；这显然是不可能实现的，其原因是前述H∞控制理论中没有涉及到“输入限制”问题，完全是一个理论的数值仿真问题。加入步进电机驱动后，实际上扰动模型ω的存在起到了限制输入的作用，使得襟翼摆角β的范围限定在 ±0.4rad的范围即±23°，这恰好在物理可实现的范围之内，这也表明另一个关于H∞控制的特点——不是所有的扰动信号都具有完全负面的作用。从图5（b）的鲁棒性能（即参数γ的波动）对比结果可以看出：理论控制下参数γ急速衰减而收敛于一个极小值范围，体现了良好的鲁棒性能，这也反映了H∞控制的理论优越性；实际上由于需要限制摆角β的范围，实际控制的鲁棒性能参数γ的波动虽然没有收敛，但其幅度维持在一个合理的波动范围内，且波动平稳，体现了良好的物理实现性能。

3.2 鲁棒控制参数讨论

鲁棒控制参数ρ对控制结果有着重要的影响。本设计以挥舞运动的高频振动控制为目的，图6为不同ρ参数下实际控制的襟翼摆角（图6（a）），及被控挥舞位移（图6（b））。从实际被控挥舞位移来看，不同参数ρ的影响并不大，均能保持良好的控制效果；从实际控制的襟翼摆角看， ρ=3 时，控制输入β维持在±0.2rad的范围即±12°，这是一个适中的范围。因为过大的输入摆角β不仅意味着过大的功耗，且容易带来物理实现中不同构件的物理干涉；而过小的输入摆角β则意味着步进电机脉冲驱动的累积误差相对变大，从而影响控制精度。所以鲁棒控制参数ρ的取值，能直接影响物理性的输入。

图6 不同ρ参数条件下实际控制的襟翼摆角以及被控挥舞位移Fig. 6 Trailing-edge flap angles and the controlled flapwise displacements of practical control under different parametersρ

另外，本设计中鲁棒控制参数ρ的取值，具有普适性。本文风速取值为U=15 m/s，已足够验证参数ρ的取值。这是由于考虑到常规叶片梁机构的实际运行工况，在湍流强度不变及特定地表粗糙度情况下，在80 m高度以内的湍流分形特性中，风速一般不超过10 m/s[21]。图7为在最优鲁棒控制参数 ρ=3及U=10 m/s条件下实际控制的被控挥舞位移z、性能参数γ波动以及控制输入的襟翼摆角β。由图7可以看出：被控挥舞位移幅度适中，频率大大降低；鲁棒性能参数波动平稳，且处于合理的幅度范围；襟翼摆角β的范围限定在 ±0.25rad的范围内，便于物理实现。综上，满足H∞鲁棒控制的性能要求。

图7 ρ=3、U=10 m∙s–1条件下实际控制的被控挥舞位移z、性能参数γ、及襟翼摆角βFig. 7 Controlled flap-wise displacement, performance parameter γ, and the trailing-edge flap angle β of practical control under conditions of ρ=3, U=10 m/s

需要说明的是：鉴于改进后的系统方程（11）的控制难度，难以找到适合的其他主动控制方法，故本设计缺少与其他主动控制方案的对比，然而与一些优秀的被动控制方案相比，可以间接论证本方案的有效性。作者所在项目组成员在文献[22]中，研究了基于形状记忆合金（SMA）的被动控制，也能获得较好的幅度控制结果；随着SMA在复合材料母体中的含量的增加，幅度控制效果甚至更加显著，但无论何种情况，基于SMA的被动控制的时域响应时间会延长，此为不利因素。故相对于本文的主动控制，其效果在伯仲之间。

4 基于工程应用的过程控制实验

由于常规工控设备的控制器（如PLC）CPU主频和内存等性能的局限、控制器独立程序编制的难度，大多数高性能的智能控制控制算法无法直接在控制器硬件中实现，因此存在着工程应用的局限。OPC技术是一种能把MATLAB/SIMULINK (MS) 仿真环境与PLC连接起来实时数据通信的技术[23–24]。为检验过程控制性能，即检验本文提出的H∞控制算法在实际的控制器硬件系统中的有效性、实时性和可靠性，本设计采用基于S7–300PLC、WinCC组态软件、MS环境的“实时OPC技术”来实现基于硬件的控制算法工程应用、数据通讯、虚拟仿真技术。

常规的OPC技术仅仅实现MS仿真环境和S7–300PLC的数据通讯功能，实时OPC技术能将智能控制算法在控制器硬件中完美实现，能将基于SIMULINK构架的算法程序转换为PLC控制器程序，从而实现与仿真环境的联机通讯和算法的执行。具体步骤如下：1）采用Simulink® PLC CoderTM工具，将嵌入SIMULINK的MATLAB代码（以S-function形式嵌入）作为一个整体，直接编译成与硬件无关的IEC 61131结构化文本（ST）；2）以ST形式生成的代码具有专业手写代码的清晰性和高效性。在此基础上作适当修改，修改为适合于S7–300PLC运行的“新ST”文本；3）将新ST文本编译并部署到S7–300PLC上直接应用，或者转化成梯形图程序后再应用；4）由于PLC扫描周期的限制，如果扫描周期过长，则可以启动冗余控制器系统来分部运行控制算法。5）H∞控制理论完全在S7–300控制器硬件中运行；叶片系统的虚拟仿真过程完全在MS环境实现；WinCC作为OPC服务器，一方面与客户端S7–300通讯，另一方面与客户端MS环境的叶片系统模型通讯。MS中的OPC toolbox模块库，具有相关的OPC读写模块，借助于OPC服务器，可实现与S7–300的读写操作。

以图5中鲁棒控制参数 ρ=1情况下的实际控制为例，在WinCC的组态界面中可以实时监控鲁棒性能参数γ和襟翼摆角β。图8为S7–300控制柜系统，包括电源模块、控制器及备用冗余控制系统；WinCC组态界面中的过程控制性能检验结果，展示了5 s内的鲁棒性能参数波动以及实时控制下的襟翼摆角的波动，对比图5中的相关项目，可以发现：襟翼摆角的实时波动与图5中的仿真结果相当一致，而鲁棒性能参数的实时波动与图5中的仿真结果基本保持一致：波动趋势与振动频率完全一致，但波动幅度位置的极大值有差别。这种差别来源于WinCC对S7–300硬件的采样间隔最小值Tmin的局限以及OPC toolbox中读写模块所定义的采样间隔时间T0的不同。其中，T0为选定常数，一般要介于叶片系统虚拟仿真的平均采样时间和Tmin之间，既能保证适中的算法运行速度，也能保证适中的控制精度。因为过高的精度，必然会消耗算法的运行时间，并可能使得S7–300控制器硬件出现内存溢出、循环扫描周期过长等缺陷，或由于累积误差过大而使控制过程失效。但从图8展示的襟翼摆角看，本实验过程的控制性能可得到有效的保证。