参数耦合下风力发电机叶片机械颤振检测研究

张 纯，孙 山，韩煜航

（1.西安热工研究院有限公司，西安 710054；2.西安交通大学 能源与动力工程学院，西安 710049）

风力发电机是转化气流动能为机械能，最终转化为电能的主要动力机械设备[1]，目前已成为不可或缺的可再生能源发电设备，在全球风电产业不断发展下[2]，风力发电机的应用场景越发广泛，但同时暴露出许多技术问题。叶片机械颤振不仅会影响风力发电机组的运行效率和可靠性，还可能导致严重的结构损坏和安全事故。因此，对于风力发电机叶片机械颤振的检测和控制显得尤为重要。

文献[3]采用非线性谐波法对风力发电机叶片实行三维非定常流动分析，得到叶面脉动压力，通过流固转换程序获取叶片所受气动激振力。采用模态叠加法实现机械颤振检测。文献[4]通过构建风力发电机叶片三维壳体模型，采用有限元理论建立风力发电机叶片运动学方程，实现机械颤振检测。以上方法均存在一定的局限性，导致气动阻尼系数计算结果不理想、检测正确率低、计算时间长的问题。

为了解决上述方法中存在的问题，提出参数耦合下风力发电机叶片机械颤振检测方法。

1 风力发电机叶片结构动力学方程推导

基于悬臂旋转欧拉－伯努利梁运动微分方程和哈密顿原理建立风力发电机旋转叶片结构模型[5-6]，如下所示：

式中：U 表示应变能；T 表示动能；W 表示虚功；δ 表示变分。

考虑复合材料正交各向异性，改进风力发电机叶片挥舞弯曲和扭转控制微分方程为如下形式：

式中：x 表示自变量；R 表示定义域；L 表示气动载荷下风力发电机叶片气动力[7]；m 表示单位长度叶片质量；e 表示质心偏移；w 表示挥舞位移；GJ 表示扭转刚度；Rm表示极回转半径；ξ 表示扭转位移；ρ 表示叶片密度；A 表示叶片横截面积；ψ 表示风轮转速；Rm1和Rm2表示绕主中心轴和垂直于弹性轴的轴质量回转半径；EI 表示挥舞弯曲刚度；P 表示气动载荷下叶片变桨距[8]；′、·和″、··表示对应参数一阶和二阶导。

基于西奥道森理论，可得到非定态气动载荷下叶片气动力：

式中：VL表示升力系数；ρ∞表示空气密度；l 表示半弦长；d 表示弹性轴与弦中点距离；v 表示来流速度；C（k）表示西奥道森函数。

离散化叶片结构模型，并采用非保守系统Lagrange方程获取风力发电机叶片结构动力学方程[9]，如下所示：

式中：q 表示广义坐标矢量；M、C 和K 表示结构质量、阻尼矩阵和刚度矩阵；QF表示结构外激励。

通过振型叠加法求解结构动力学方程[10]，基于特征方程Kφi=Mλiφi求解特征值λ1，λ2，…，λM和特征向量φ1，φ2，…，φM，忽略系统高阶振型，采用保留的m 个低阶振型构建振型矩阵Ω＝［o1，o2，…，om］，由此将式（4）风力发电机叶片结构动力学方程近似为下式形式：

2 风力发电机叶片机械颤振检测方法

2.1 CFD/CSD 参数耦合算法

采用非保守系统Lagrange 方程，推导出叶片的结构动力学方程，能更准确地描述风力发电机叶片的动力学行为。将非稳态气动载荷与叶片的结构参数进行耦合求解，得到叶片的结构响应和振动情况。

在式（4）风力发电机叶片结构动力学方程中，右端项QF即为非稳态气动载荷。将欧拉方程表示为积分方式[11]，如下所示：

为了推进非稳态时间时，引入伪时间步τ，并利用双时间步长方法求解欧拉方程，获得叶片的结构响应，具体表达式为

在CFD 模拟中，使用CFD 技术对风流场数值模拟，计算出时间变化的气流场分布，并获得每个时间点上的气动载荷分布[12]。在CSD 模拟中，根据当前时间的气动载荷和先前时间步的结构响应计算叶片的位移、速度和加速度等参数。利用交错时间推进方式，将气动载荷传递给结构动力学模型，反馈得到叶片的结构响应。具体流程如图1 所示。

图1 CFD/CSD 参数耦合算法流程Fig.1 CFD/CSD parameter coupling algorithm flow chart

为了减小标准CFD/CSD 参数耦合算法误差，并提高计算精度，需要对该算法加以优化：在风力发电机叶片气动弹性计算中，非稳态气动力和结构位移通常表现出连续变化的趋势。由此对上图步骤（1）加以改进，基于tm-2、tm-1和tm3 个时刻气动力实行三点抛物插值，生成tm～tm+1内气动力关于时间的二次函数，并将其代入式（4）中，可实现对tm+1时刻的结构参数求解。通过这种优化方法，能够有效提升CFD/CSD 参数耦合算法的计算精度。这样的优化可以使非稳态气动力和结构位移在时间上更加连续，从而更准确地反映风力发电机叶片在风流中的响应特性。

2.2 叶片机械颤振边界计算

为了进一步提升CFD/CSD 参数耦合算法计算效率，提出基于沃尔泰拉级数的非稳态气动降阶模型[13]。沃尔泰拉级数为无穷级数，在系统受到小扰动时，可将非稳态气动力转换为二阶沃尔泰拉级数形式，其在时域内离散形式如下所示：

式中：u（·）表示输入值；n 表示离散时间变量；k、k1和k2表示时移；y（n）表示响应值；hi（·）表示系统的i阶沃尔泰拉核；h0表示稳态项。

基于阶跃响应对非稳态气动力加以辨识，构建CFD/CSD 参数耦合系统阶跃响应如下所示：

式中：δ（n）表示离散域内的单位脉冲信号；ζ0表示单位阶跃激励。

辨识后可得部分非对角分量和全部二阶核对角分量的近似一阶核（n），如下所示：

经由以上方法可生成非稳态气动降阶模型，将该模型作为广义位移输入，输出广义力，进而计算QF，并利用模态叠加法计算叶片结构动力学方程，通过调整叶片振型缩放系数，使每个振型的响应达到最大值，最终获取到叶片机械颤振边界。

2.3 叶片机械颤振检测

采用能量法在能量交换的角度分析风力发电机叶片是否存在机械颤振[14-15]，由于风力发电机叶片机械结构气动阻尼远远高于机械阻尼，因此，气动力做功可作为叶片振动系统自外界获取到的主要能量，通过检测非稳态气动力在一个振动周期内对风力发电机叶片的做功情况可确定是否存在叶片机械颤振。用Tb表示振动周期，（x，t）表示风力发电机叶片表面任意点，P（x，t）表示静压，n（x，t）表示单位外法向量，V（x，t）表示速度矢量，则周期Tb内单位面积上非稳态气动力做功W 如下所示：

阻尼的主要作用是消耗能量，为便于计算，简化振动系统中非线性阻尼为线形黏性阻尼。在流场绕流情况较好且系统中不存在较强非线性时，可采用黏性阻尼ζ 等效近似替代气动阻尼，如下所示：

式中：∑W 表示非稳态气动力做功之和；qcfd表示正则坐标下模态振幅。

基于沃尔泰拉级数的非稳态气动降阶模型和CFD/CSD 参数耦合算法实现非稳态计算后，进一步计算一个振动周期内非稳态气动力做功和等效模态气动阻尼比，若非稳态气动力做功和等效模态气动阻尼比为正，则风力发电机叶片存在机械颤振情况，需要对其进一步关注分析，以免对风力发电机安全运行产生威胁。

3 实验与结果

为了验证参数耦合下风力发电机叶片机械颤振检测研究整体有效性，需要测试参数耦合下风力发电机叶片机械颤振检测研究。

为模拟风力发电机叶片机械颤振效果，实验所用设备及实验流程如下所示：

（1）实验设备：鼓风机、调速器、叶片、压电薄 膜、LabVIEW 数据采集器、NI PXI-1042Q、上位机等。

（2）实验流程：利用鼓风机生成气流并通过风洞加速，采用调速器调节风速使其符合实验要求，在风洞出口0.1 m 处固定叶片，使用压电薄膜采集叶片信号，经由LabVIEW 数据采集器上传至上位机，计算叶片机械颤振情况作为实验实际值。

通过不断调节各项条件，采集不同实验环境下的实验样本，利用所提方法验证不同实验中叶片机械颤振情况，并与实际值加以对比，以评价方法在风力发电机叶片机械颤振检测中的性能。

通过所提方法获取不同叶间相角差下气动阻尼系数，并与实际值加以比较，结果如图2 所示。

图2 气动阻尼系数检测结果Fig.2 Aerodynamic damping coefficient test results

由图2 可以看出，所提方法得到的气动阻尼系数整体更接近于实际值。气动阻尼系数与风力发电机叶片机械颤振检测息息相关，准确的气动阻尼系数是有效识别机械颤振的基础。

判定不同条件下机械颤振情况是风力发电机叶片机械颤振检测的主要目的，采用检测正确率作为指标评估所提方法在颤振孕育状态和颤振爆发状态中的检测表现，结果如表1 所示。

表1 机械颤振检测结果Tab.1 Mechanical flutter test results

由表1 可以看出，所提方法检测正确率高达90%以上，适合用于实际检测之中。

验证所提方法在风力发电机叶片机械颤振检测中的响应效率，检测结果如表2 所示。

表2 计算时间检测结果Tab.2 Calculate time detection results

由表2 可以看出，所提方法在风力发电机叶片机械颤振检测中的响应效率明显较高。所提方法基于沃尔泰拉级数构建非稳态气动降阶模型，进而有效提升了CFD/CSD 参数耦合算法计算效率。

4 结语

风力发电机旋转叶片结构动力学方程，利用优化后CFD/CSD 参数耦合算法求解流体和结构动力学方程，建立非稳态气动降阶模型求解颤振边界，计算风力发电机旋转叶片动态载荷，利用能量法实现风力发电机机械颤振检测。该方法能够获取到更为准确的气动阻尼系数、检测正确率更高、计算时间更短，为风力发电机的进一步发展和应用奠定基础。