大型风力机叶片气弹颤振特性分析

曲珍壮，夏鸿建，李德源，林伟豪

（广东工业大学机电工程学院，广东 广州 510006）

0 引言

为了增强风力机的风能捕获能力，降低风力发电成本，风力机叶片的尺寸正逐渐变大，致使风力机叶片柔性增强，在惯性力、气动力的综合影响下，大型风力机叶片的弹性大变形和气动载荷耦合特性增强，气弹失稳导致的颤振问题突出。

近年来，对于大型风力机叶片的颤振特性分析，国内外学者进行了大量的研究,主要方法可分为时域分析方法和频域分析方法。时域分析方法通过判断振动位移的发散趋势，确定叶片颤振边界的风速和转速。如王晓东[1]采用FAST 软件，在时域条件下对NREL 5 MW 风力机叶片在运转和静止时的颤振性能进行了分析预测，并研究了扭转刚度、质心位置、结构阻尼和桨距角对颤振临界风速的影响。时域分析方法虽然可以快速判断颤振边界，但难以获得引起叶片颤振的结构振动模态和气动阻尼。频域分析方法，通过对气弹耦合方程线性化，并采用复特征值分析方法，判别颤振边界，能获得引起颤振的结构模态和气动阻尼变化。如HANSEN[2]采用（Multi-Blade Coordinates，MBC）坐标和模态法，结合叶素动量理论和B-L 模型，分析了NREL 5 MW 风力机颤振临界转速；李德源[3]采用多体动力学方法和叶素动量理论的方法建立叶片气弹耦合方程，研究了NREL 5 MW 风力机气弹稳定性问题；ZHOU[4]采用混合变分几何精确梁模型和叶素动量理论，分析了5 MW 风力机叶片的弯扭耦合颤振特性。

然而，为获得相对精确的气弹颤振边界，分析引起颤振的结构模态耦合特点和气弹耦合特性，需能高精度描述叶片大变形特点的非线性结构模型，以及结构振动与气动耦合的非定常气动模型。传统采用欧拉梁结构模型[5-6]或模态法[7]构建叶片结构模型，难以表达柔性叶片的几何非线性变形，以及转动产生的动力刚化现象。对转动叶片非线性结构建模，主要有几何精确梁模型[4,8]、绝对节点坐标模型和超级单元模型[3,9]等。其中，几何精确梁方法由于具有梁模型的自由度少、精度高，并且能准确表示非线性大变形的优点，在直升机叶片和风力机叶片获得广泛应用。然而，传统几何精确梁模型由于存在截面姿态角度插值和动力学线性化相对繁琐等问题，多数应用于大型叶片的时域响应分析。另一方面，大型柔性叶片振动会影响叶片气动载荷，气动载荷又反馈作用于结构，从而形成气弹耦合特性。传统的准静态的叶素动量理论采用静态气动系数，计算气动载荷，忽略了叶片振动引起的叶片气动非定常特性。为有效描述叶片的非定常特性，学者们分别提出了Boeing-Vertol 模型，ONERA 模型，Øye模型以及半经验B-L（Beddoes-Leishman）[10]等动态失速模型。由于B-L 模型具有求解方便，并且容易线性化的特点，获得广泛应用。

为此，本文针对大型风力机叶片结构非线性以及气动载荷的非定常特点，采用几何精确梁模型与集成B-L（Beddoes-Leishman）动态失速模型的叶素动量方法，建立风力机叶片的气弹耦合动力学方程。为进一步提高气弹模型的求解精度，便于对叶片结构方程线性化，由弱形式的求积元方法，高效构建叶片的气弹耦合动力学离散平衡方程；再基于变分原理，综合线性化气弹耦合动力学方程，由复特征值分析法，计算叶片的气弹频率和气弹阻尼比，辨识叶片气弹颤振区域；并以DTU-10MW 风力机叶片，研究大型风力机叶片的气弹颤振特性。

1 几何精确梁变形的描述

几何精确梁采用平截面假设，梁的位形由梁轴线的位移和相应截面的转动决定。如图1 有全局坐标系（不考虑叶片预弯预扭，以及截面偏心等因素），不变形位置坐标系（考虑叶片真实几何形状）和变形位置坐标系（考虑气弹耦合后叶片的变形）。

图1 梁的转动坐标系

在变形位置，梁上任意一点的位置向量可以写为：

式中：

u——剪切中心在全局坐标系下的坐标；

z——截面在梁轴线上的坐标，m；

x——截面在x轴的坐标，m；

y——截面在y轴的坐标，m；

Λ——截面转动矩阵；

——截面坐标系下的位置向量。

三维旋转矩阵Λ有9 个元素，包含6 个约束关系，难以直接作为姿态变量，欧拉四元数是描述截面转动非奇异的最小参数表达，本文采用欧拉四元数描述截面的转动[11]。

式中：

q0——四元数的标量部分；

q——四元数的向量部分。

四元数和转动矩阵Λ的关系为：

式中：

I——单位矩阵。

2 结构动力学方程的建立

2.1 几何精确梁本构关系

几何精确梁理论基于Reissner 应变，根据参考文献[12-13]，定义Reissner 应变Γ为平动应变向量，应变K为转动应变向量，截面等效应力N，M和应变Γ，K之间的关系为：

式中：

S——截面刚度矩阵。

1.3 样本采集及病原菌鉴定 观察病害症状，用数码相机拍照，填写样方调查表，采集典型样本压制成标本，带回实验室采用科赫法则对病原菌进行验证和鉴定。对白粉病、锈病等病症明显的病害，直接标本制片或以徒手切片镜检，参照相关资料[10-19]进行鉴定；对病症不明显的病害，进行病原菌常规组织分离[20]，获得纯培养的病原菌，再进一步鉴定并做致病性测定。

风力机叶片由于采用复合材料铺层和复杂的截面形状，截面耦合特性明显。本文采用专业的截面分析软件VABS 计算截面刚度矩阵S。

应变与变形的关系：

式中：

Γ0——梁的初始位置；

——梁的初始曲率参数；

（）'——在全局坐标中对截面在Z 轴上坐标的导数；

（）*——四元数的共轭；

‘°’——四元数乘法。

2.2 动力学平衡方程

首先，建立风力机叶片弱形式动力学平衡方程，然后通过数值积分和微分方法离散求解。

对内力沿单元轴线积分，得到内力虚功为：

式中：

δ——变分符号；

e——叶片单元；

le——叶片单元长度。

对气动力沿轴线积分，重力进行体积积分，得到单元的外力虚功：

式中：

δua——轴线气动节点的虚位移；

δθa——轴线节点气动点的虚转动；

fa——轴线节点的气动力向量；

ma——轴线节点的气动弯矩向量；

g——重力加速度向量；

ρ——叶片体密度，kg/m3。

各点的速度进行体积积分，得到单元的惯性力虚功：

式中：

——位置向量对时间求二阶导。

根据虚功原理，单元内力，外力和惯性力虚功满足：

为求解平衡方程，需要通过数值积分与数值微分方法，转换为节点求和形式，本文采用Gauss-Lobatto 积分[15]进行数值离散，Gauss-Lobatto 积分一般形式为：

式中：

积分区域需转换到标准区域[-1，1]；

ξ——标准化坐标；

f（ξ）——被积函数；

n——积分点个数；

ξi——积分点；

ωi——对应的积分权系数。平衡方程中还存在对叶片轴线坐标z 的导数项，需要采用微分求积法[15]，转换为单元节点加权求和形式，微分求积的一般形式为：

式中：

f（ξj）——函数在离散点上的值；

——微分求积权系数。

通过Gauss-Lobatto 积分和微分求积，得到基于节点虚位移表达的内力，外力和惯性力虚功[8]：

节点k虚位移表示为：

叶片的一个单元位移δde=［δd1δd2···δdK］。本文风力机叶片每个单元包含21 个节点。

把公式（13）代入公式（10），并在全局坐标系下进行组装，得到等效节点内力Rint、外力Rext和惯性力Rinert的离散平衡方程：

3 动力学方程线性化

为了分析叶片的振动频率和模态特性，需要对动力学平衡方程进行线性化处理，分为2 个部分：结构方程线性化和气动力线性化。

3.1 结构动力学线性化

分别对节点内力，外力和惯性力矩阵线性化得到[8]：

将公式（16）-（18）代入到平衡方程（15）可以得到：

式中：

K=Kint-Kext+Kinert，M=Minert，C=Cinert；

μ——弱形式求积长度系数；

B——微分求积定位矩阵。

3.2 气动力部分线性化

根据叶素动量定理，作用在叶片第false个截面的气动力和气动力矩可以写为：

式中：

c——截面弦长，m；

ρa——空气密度，kg/m3；

W——相对来流风速，m/s。

C dyn=（CL，CD，CM）T；其 中CL、CD、CM分 别是升力系数，阻力系数和扭矩系数，叶片翼型截面如图2 所示。

图2 叶片翼型截面

考虑叶片振动和入流角的影响，相对来流风速W和攻角α可以写为：

式中：

Vn——垂直旋转平面的来流风速，m/s；

——摆振方向振动速度，m/s；

——挥舞方向振动速度，m/s；

R——翼型截面剪心到轮毂中心的距离，m；

Ω——叶片旋转角速度，rad/s。

β——桨距角，（°）；

θ——叶片扭角，（°）。

对气动力求变分

式中：

系数矩阵Tx，Tc，Tα，Tα.的具体表达可以参考文献[5,16]；

——攻角对时间的导数；

α3/4——叶片翼型截面3/4 弦线处的攻角，求变分得：

式中：

L0——翼型截面3/4 弦点到弹性中心的距离。

为了考虑涡旋脱落和尾缘流动分离对气动特性的影响，引入基于状态变量xa=｛x1x2x3x4｝T的B-L动态失速模型[16]，修正叶片各翼型截面处的升力系数及其斜率。

对状态变量线性化得到

式中：

Ms，Ks，Cs——状态变量线性化引入的质量阵，刚度阵和阻尼阵；

Ad——系数矩阵可以参考文献[16]。

将式（24）-（26）代入（23）中得到气动载荷的线性化表示：

Ka，Ma，Ca——气动部分线性化得到的气动质量阵，刚度阵和阻尼阵。

联立式（19），（28）和（29），气弹线性平衡方程可以写为：

4 特征值复模态分析

根据式（30）建立风力机颤振分析的特征值平衡方程：

利用Matlab 软件提供的eig函数计算特征值和特征向量，得到叶片的气弹频率、气弹阻尼比和振型。其中特征值的虚部为气弹频率，特征值的实部为气弹阻尼比。

5 计算结果

5.1 结构频率验证

不考虑风力机叶片结构阻尼和气动力作用的特征方程为：

式中：

K——公式（19）中的结构质量矩阵；

M——公式（19）中的刚度矩阵。

通过求解特征方程可以获得叶片的固有频率ωn和模态。

为了验证所建立的叶片结构的动态精度，以DTU-10 MW 风力机叶片为研究对象，计算了叶片的固有频率，并与文献[17] 中发表的结果进行了比较。对比结果如表1 所示，最大误差为3.90%，说明本文建立的叶片模型是正确的。图3 还绘制了第一阶挥舞振型、第一阶摆振振型和第一阶扭转振型模态图。

表1 风力机叶片固有频率对比

图3 叶片的模态图

5.2 运行工况下气弹频率和气弹阻尼比

根据参考文献[17]提供的叶片数据和运行工况参数，计算了DTU-10MW 风力机叶片在运行工况下的低阶气弹频率和气弹阻尼比，并与HANSEN[18]的结果进行对比。从图4 中可以看出，在运行工况下，一阶摆振频率基本保持不变，和HANSEN 的结果也比较吻合。

图4 运行工况下一阶摆振气弹频率和阻尼比

一阶摆振气弹阻尼比和HANSEN 相比，在风速4～11 m/s 区间趋势基本相同，都处于缓慢上升趋势，风速11 m/s 之后本文的结果在略微下降后有比较明显的上升趋势，由于和HANSEN 采用的模型不同，以及截面耦合矩阵S的计算软件的不同，所以在结果上有一定差异。从图5可以看出，在运行工况下，一阶挥舞气弹阻尼，在风速4～7 m/s的保持稳定，在风速7～12 m/s的处于上升趋势，之后趋势保持平稳，此时叶片达到额定转速9.6 r/min，结果和HANSEN的较为接近，挥舞方向气弹频率在到达额定转速后，也保持稳定。

图5 运行工况下一阶挥舞气弹频率和阻尼比

5.3 定叶尖速比工况下颤振特性分析

为了研究风力机叶片的颤振极限，对风力机叶片在大变形和高转速颤振特性进行分析。将风力机叶片运行工况设为定叶尖速比10，即风力机叶尖速度为来流风速的10 倍，并将叶片的桨距角设为0。分析结果如图6，在叶片转速0.6～2.2 rad/s 之间，挥舞方向气弹频率和和摆振方向气弹频率都在增大，但挥舞方向变化率明显比摆振方向大，分别为106%和9%，这种现象也和之前对5 MW 风力机颤振特性的研究结果相似[3-4]。挥舞方向气弹阻尼比为较大的正值，摆振方向气弹阻尼比下降到较小的正值后保持稳定，扭转方向气弹频率先缓慢上升，然后呈下降趋势，扭转方向气弹阻尼比在转速1.77 rad/s 时出现负值，对应的气弹频率为8 Hz，且随转速增加持续下降，说明此时风力机叶片发生颤振失稳现象。

图6 定叶尖速比工况下气弹频率和气弹阻尼比

6 结论

本文采用非线性几何精确梁模型和叶素动量定理建立风力机叶片的气弹模型，在B-L（Beddoes-Leishman）动态失速模型的基础上对气弹模型线性化，构建特征平衡方程，通过特征值分析，计算风力机叶片的气弹频率和气弹阻尼比，主要结论如下。

在运行工况下，一阶摆振气弹频率保持稳定,受叶片模型和耦合等因素影响较小。一阶摆振气弹阻尼比变化趋势为：在额定转速位置上升到最高值，经历短暂的下降后继续上升。本文的结果和HANSEN相比，上升趋势明显，表明一阶摆振气弹阻尼比受模型和截面耦合影响较大。一阶挥舞气弹频率和气弹阻尼比都是在额定转速位置上升到最高值，然后保持稳定趋势。

在定叶尖速比工况下，随着转速的增加，风力机叶片挥舞方向气弹频率快速增加，摆振方向气弹频率变化较小，摆振方向气弹频率缓慢下降到较小的正值后保持平稳，没有出现负值。一阶扭转方向气弹频率在转速1.77 rad/s 时，出现负值，并且继续下降，风力机叶片出现颤振失稳现象。

本文建立的气弹模型，还可以用来分析风力机叶片其它的气弹问题。