后掠型风机叶片作为一种自适应被动降载的气动布局，近年来受到了研究者们的广泛关注^[1]。由于积叠线弯曲，叶片在运转状态下的三维流动显著，在结构方面弹性轴与扭转轴相分离，致使叶片在复杂气动载荷作用下的气弹耦合特性愈加显著。风力机在运行工况下，突遇极端风况，如风速或风向的突然变化，而控制系统(如变桨距系统)来不及响应的情况下，会使得叶片的气动载荷在极端风况作用期间和之后较正常风况有较大的变化。在IEC 61400-1^[2]标准和GL 2010^[3]准则提出的6种极端风况模型中，极端运行阵风是基于轮毂高度建立的阵风模型，在其作用下风力机载荷产生较大变化，可能会危及风力机正常运行并对和对电网造成冲击。因此极端风况下的后掠风力机气动特性研究对机组的安全稳定运行和优化设计具有重要意义。

对于风力机在极端阵风工况下叶片的响应特性，国内外进行了相关研究。李仁年等^[4-5]基于CFD (Computational Fluid Dynamics)方法研究了均匀来流和极端运行阵风对风力机非定常气动性能的影响，并利用Fluent技术对极端风况下的风力机流场进行模拟，分析叶片非稳态、稳态下失速差异形成。莫文威等^[6]基于多体模型进行了水平轴风力机气弹耦合分析；汪显能等^[7]在此基础上进行了非定常气动特性分析。张林伟等^[8]结合叶片气弹响应，在极端阵风条件下研究风力机动态失速、阵风上升时间、变桨速率等对叶片极限载荷的影响。张文广等^[9]研究了风力机智能叶片的非定常气动特性，发现尾缘襟翼可以减小动态失速效应，有利于降低疲劳载荷、抑制功率波动。丁建刚^[10]研究了在极端载荷作用下的风力机耦合振动稳定性。Bierbooms等^[11]对变桨调节风电机组在阵风极端上升时间下的极限载荷进行研究。后掠叶片是近年来新出现的叶片形式，关于后掠叶片在极端运行阵风下气动特性研究还鲜见报道。

本文通过“超级单元”^[12]方法将后掠叶片离散为若干个由运动副和力元联接的刚体，应用Roberson-Wittenburg递推建模方法，结合螺旋尾涡升力线模型(Helicoidal Wake Vortex Lift Line Model)，建立了后掠叶片非线性气弹耦合方程。以NREL 5 MW风力机后掠叶片为例，研究了极端运行阵风及其持续时间对后掠叶片气动性能的影响。

1 极端运行阵风模型(Extreme Operating Gust，EOG)

依据IEC 61400-1标准建立极端运行阵风模型，模拟指定高度处N年重复周期内的极端运行阵风风速变化。

(1) 标准风机级N年一遇轮毂高度阵风值 ${\upsilon _{{\rm{gust}},N}}$ 由式(1)给出。

${\upsilon _{{\rm{gust}},N}} = \beta \Bigg(\dfrac{{{\sigma _1}}}{{1 + 0.1(D/{{\varLambda} _1})}}\Bigg) $

(1)

式中： ${{\upsilon} _{{\rm{gust}},N}}$ 为极端运行阵风的风速，预计N年出现一次； ${\sigma _1}$ 为轮毂高度处的纵向风速标准偏差；D为风轮直径； ${\varLambda _1}$ 为湍流尺度参数， $\;\beta $ =4.8，N=1； $\;\beta $ =6.4，N=50。

(2) 轮毂高度处的纵向风速标准偏差由式(2)给出。

${\sigma _1} = {I_{15}}(15 + a{\upsilon _{{\rm{hub}}}})/(a + 1)$

(2)

式中： ${I_{15}}$ =0.18， $a$ =2，较高湍流强度的类别； ${I_{15}}$ =0.16， $a$ =3，较低湍流强度的类别； ${\upsilon _{{\rm{hub}}}}$ 为轮毂高度处的风速。

(3) ${\varLambda _1}$ 湍流尺度参数由式(3)给出。

${\varLambda _1} = \left\{ \begin{array}{l} 0.7{h_{{\rm{hub}}}}\;\;\;\;\;\;{h_{{\rm{hub}}}} < 30\;{\rm{m}} \\ 21\;{\rm{m}}\;\;\;\;\;\;\;\;\;\;{h_{{\rm{hub}}}} \geqslant 30\;{\rm{m}} \end{array} \right.$

(3)

式中 ${h_{{\rm{hub}}}}$ 为轮毂中心垂直高度。

(4) 极端运行阵风模型如式(4)所示。

$\upsilon (h,t)\! = \!\left\{ { \begin{aligned} & {\upsilon (h)\! - \!0.37{\upsilon _{{\rm{gust}},N}}\sin (3\pi (t - {t_0})/({t_1} - {t_0}))\text{×}}\\& {[1 \!-\! \cos (2\pi (t - {t_0})/({t_1} - {t_0})\;)],\;\;\;\;{t_0} \leqslant t \leqslant {t_1}}\\& {\upsilon (h),\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\qquad\quad{\rm{else}}} \end{aligned}} \right.$

(4)

式中： $\upsilon (h) = {\upsilon _{{\rm{hub}}}}(h/{h_{{\rm{hub}}}}){\;^\alpha }$ ，h为叶片位置处的垂直高度， $\alpha $ 为幂指数取0.12； ${t_0} \leqslant t \leqslant {t_1}$ 为阵风作用时间(让气弹耦合模型先计算t₀进入稳定状态再引入极端运行阵风模型，本文t₀取35 s)。

本文采用文献[13]的5 MW水平轴风力机后掠叶片为研究对象。考虑风速纵向分布规律(即风剪切)，轮毂中心风速取风力机额定风速11.4 m/s，忽略重力影响，选取极端运行阵风模型为1 a一遇、较高湍流强度类别，即 ${I_{15}}\!=\!0.18$ ， $a\!=\!2$ ， $\;\beta\!=\!4.8$ 。阵风作用时间设为8，10，12 s，极端运行阵风风速模型如图1所示。

2 风力机后掠叶片气弹耦合模型

将后掠叶片离散为由运动副和力元联接的多体系统^[14]，能用较少的自由度准确地展现后掠叶片的非线性振动特性，可实现程序化动力学模型建立并数值求解。为了准确描述后掠叶片的变形，引入“超级单元”^[12]对后掠叶片建立多体系统动力学模型。每个单元由4个刚体组成,刚体之间通过运动副、弹簧和阻尼器连接。

应用计算多体动力学中的R-W方法，取各铰的转动为广义坐标，则从叶片的拓扑构形可得其广义坐标阵为

${{q}}(t) = {({q_1}{\rm{(t),}} \cdot \cdot \cdot {\rm{,}}{q_n}(t))^{\rm{T}}}$

(5)

${{{q}}}(t)$ 对时间的一阶和二阶导数 $\dot {{q}}(t)$ 、 $\ddot {{q}}(t)$ 分别称为广义速度和广义加速度，n为系统自由度。叶片的动力学方程具体形式为

$ {{{Z}}\ddot { q}}+{ {\varPhi }}_{{q}}^{\rm{T}}{{\lambda} }={{z}} $

(6)

其中， ${{Z} }$ 称为广义质量阵， ${{z} }$ 称为广义力阵，包含气动力等相关项， ${{ \varPhi }_{{q}}}$ 为约束方程的雅可比矩阵， ${\lambda }$ 为拉格朗日乘子阵，拉格朗日乘子项 ${ \varPhi }_{{q}}^{\rm{T}}{{ {\lambda } }}$ 表示刚体间所有理想约束力的贡献，其具体表达形式可参阅文献[14]。求解时，需要运动约束联立求解。

式(6)中的广义力阵 ${{z}}$ 表达了叶片所受的空气动力等载荷，考虑到后掠叶片展向流动显著和叶片积叠线弯曲的特点，选用升力线理论^[15]结合刚性尾涡模型建立叶片气动模型，通过引入B-L^[16]动态失速模型进行非定常气动载荷计算。

升力线模型的核心思想是将风力机叶片三维流场中的涡量分布简化为一根涡量沿展向变化的线涡，表征叶片与流场之间的相互作用。该线涡代表了三维流场中真实的叶片环量分布，通过计算附着涡上控制点的诱导速度，获得该位置处的有效攻角，从而实现翼型气动力计算。

计算出控制点诱导速度后，通过二维翼型数据表插值得到升阻力系数C_L和C_D后可以通过沿着叶片径向积分得到轴向推力F、转动力矩Q和转动功率P。

$ F = \displaystyle\mathop \int \limits_{{\xi _{{\rm{hub}}}}}^1 \left( {\dfrac{1}{2}\rho NcW{'^2}} \right)\left( {{C_L}{\rm{cos}}\; \phi ' + {C_D}{\rm{sin}}\; \phi '} \right)R{\rm{cos}}\; \psi {\rm{d}}\xi $

$ Q \!= \!\displaystyle\mathop \int \limits_{{\xi _{{\rm{hub}}}}}^1 \!\!\left(\! {\dfrac{1}{2}\rho NcW{'^2}} \!\right)\!\left(\! {{C_L}{\rm{sin}}\;\phi ' \!-\! {C_D}{\rm{cos}}\;\phi '} \!\right){R^2}{\rm{cos}} \;\psi {\rm{cos}}\; \varLambda \xi {\rm{d}}\xi $

$ P \!= \!\mathop \int \limits_{{\xi _{{\rm{hub}}}}}^1 \!\!\!\left(\!\! {\dfrac{1}{2}\rho NcW{'^2}} \!\!\right)\!\!\left(\! {{C_L}{\rm{sin}}\;\phi ' \!\!-\! {C_D}{\rm{cos}}\; \phi '} \!\right)\!{R^2}\Omega {\rm{cos}} \; \psi {\rm{cos}} \; \varLambda \xi {\rm{d}}\xi $

$ \phi ' = {\rm{ta}}{{\rm{n}}^{ - 1}}\dfrac{{{V_0} + {V_h}}}{{\varOmega r'{\rm{cos}}\; \varLambda + {V_b}}} - {\alpha _i} $

$ W' = \frac{{\sqrt {{{({V_0} + {V_h})}^2} + {{\left( {\varOmega r + {V_b}} \right)}^2}} }}{{{\rm{cos}}\;{\alpha _i}}}\;\;\;{\alpha _i} = \frac{{{W_n}}}{W} $

(7)

其中，ρ为空气密度，N为叶片数，c为截面弦长， $ W $ 为未受扰动相对风速， $W' $ 为截面位置受扰动的相对风速， ${W_n}$ 为诱导速度， $\phi '$ 为实际入流角，R为风轮半径， $\psi $ 为风轮锥角， ${{{\varLambda}}}$ 为当地后掠角， $\xi $ 为叶片展向无量纲单位化参数，V_h和V_b为挥舞、摆振方向振动速度， ${\alpha _i}$ 为诱导攻角，Ω为风轮转速。

3 算例与分析

算例模拟文献[13]的5 MW水平轴风力机后掠叶片在额定工况运行下，在t₀时间突遇极端运行阵风，由式(4)分析阵风风速，通过式(6)和式(7)耦合联立求解，可得出阵风作用下的叶片动力响应。分析流程如图2所示。

3.1 叶根挥舞力矩

叶根力矩的时域响应能体现叶片整体的气动特性，图3是叶根处的挥舞力矩时域响应曲线。从图3可以知道50~80 s为风力机正常工作时间。风轮正常工作时叶根挥舞力矩的均值约为9×10⁶ N·m，叶根挥舞力矩具有周期性，同时存在明显的高频振动。

由于风轮以固定转速运动，在风剪切作用下气动载荷的周期性与风轮转动周期一致；在叶根附近的翼型处于失速状态，在计算气动载荷时耦合上叶片的振动和变形，造成了气动载荷的高频变化。在极端运行阵风作用期间，8，10，12 s极端运行阵风作用时间的叶根挥舞力矩的峰值为1.155×10⁷ ，1.130×10⁷，1.116×10⁷ N·m，与均值相比分别增长了近23.9%，23.6%，23.5%。从阵风用作时间上看，极端运行阵风作用时间越短，叶根挥舞力矩的峰值越大。

3.2 功率

图4为极端运行阵风下风轮功率时域响应曲线。在35 s开始计入极端运行阵风模型，在50 s后都基本恢复到正常工作状态。风轮在正常工作时，由于风剪切作用，功率具有明显的周期性，高频振动体现出叶片气动载荷计算与叶片结构振动的耦合和部分翼型的动态失速。风轮正常工作时的功率均值为5.15 MW，与设计值5 MW相近，验证了建立的风力机气弹耦合模型准确可信。在极端运行阵风作用期间，8，10，12 s极端运行阵风作用时间的风轮功率峰值分别为9.918，9.565，8.73 MW，与功率均值相比分别增长了92.6%，85.7%，69.5%，可以看出极端运行阵风作用时间长短对功率的影响比较明显，阵风作用时间越短，功率的峰值越大。

3.3 叶片展向65%位置处的攻角与升力系数

图5为后掠叶片65%展向位置处翼型截面的攻角时域响应曲线。50 s以后可以认为恢复到稳定状态，正常工作状态下该位置处的攻角在4.61°~6.12°之间变化，均值约为5.3°。曲线有明显周期性和高频振动，是由于气弹耦合和部分翼型失速造成的。在极端运行阵风作用期间，在8，10，12 s极端运行阵风作用时间下，攻角的峰值分别为9.91°，9.77°和9.91°，微小的差异可能是由于风剪切造成的。

图6为后掠叶片65%展向位置处翼型截面的升力系数时域响应曲线。在正常工作状态下，该翼型升力系数变化范围在1.04~1.22，均值约为1.14，曲线有明显周期性和高频振动。在极端运行阵风作用期间，在8，10，12 s极端运行阵风作用时间下，升力系数的峰值分别为1.47，1.46，1.46，峰值基本一致。该翼型静态失速攻角为9°，对应的升力系数为1.403，由于发生了动态失速延时，攻角在超过静态失速攻角9°后，升力系数没有随攻角增大而下降，还继续增大。

3.4 轴向推力

图7为风轮所受到的轴向推力时域响应曲线。在50 s以后恢复到正常工作状态，曲线有明显周期性和剧烈高频振动，变化范围在6.5×10⁵~7.5×10⁵ N之间，风轮轴向推力的均值为7.05×10⁵ N。在极端运行阵风作用期间，在8，10，12 s极端运行阵风作用时间下，轴向推力的峰值分别为9.47×10⁵，9.28×10⁵，9.15×10⁵ N，与均值相比，分别增长了34.3%，31.6%，29.8%，极端载荷较正常工况下的载荷变化较大，在风力机安全设计上需要充分考虑。

4 结论

(1) 通过引入“超级单元”建立风力机后掠叶片多刚体系统模型，结合升力线气动模型组成的风力机气弹耦合模型，能较好地展现出风剪切与风机叶片气动弹性耦合问题。研究工作对风力机的结构优化设计和安全稳定运行具有重要作用。

(2) 在极端阵风工况，风力机没来得及变桨距情况下，会极大影响风力机气动性能，部分翼型进入失速状态；极端载荷可能会对风力机叶片造成一定的损坏，在风力机安全设计上需充分考虑，同时加强天气预报工作。

(3) 极端运行阵风作用时间越短，极端载荷越大。后掠叶片变桨距控制系统需要结合当地风电场阵风特点，选择合理的变桨策略来降低阵风对风力机叶片载荷与电网的影响。