运行状态风力机地震响应的解耦分析方法

2022-09-03席仁强许成顺杜修力

工程力学 2022年9期

席仁强，许成顺，杜修力，许坤

(1. 常州大学机械与轨道交通学院，江苏，常州 213164；2. 北京工业大学城市与工程安全减灾教育部重点实验室，北京 100124)

为应对环境挑战，风能在现代能源体系中所占比重不断扩大[1]。20 世纪90 年代，风力发电率先在欧洲兴起；进入21 世纪，随着兆瓦级风力机的应用，风力发电从北欧扩展至世界各地[2-3]。近年来，越来越多的风电场在高地震危险性区域建设，风力机抗震受到研究者的关注[4]。对于多兆瓦风力机，叶片柔度大，当叶尖位移峰值为1 m 时，其法向振动速度峰值可达6 m/s[5]。因此，叶轮气动荷载计算需考虑气弹性效应[6]。根据气弹性模型，风力机结构地震响应有耦合、解耦两类分析方法。

耦合方法直接考虑风力机动力响应的气动-伺服-弹性耦合行为。Witcher[7]利用BLADED 软件，开展了风-地震联合激励下风力机动力反应分析，结果表明：该软件能够模拟气弹性和控制器的影响。随后，Prowell 等[8]对开源软件FAST 进行二次开发，使其能够实施风-地震作用下风力机动力反应分析，他们还通过振动台试验和数值模拟验证该方法的可靠性[9]。改进的FAST 程序被广泛用于分析风力机地震反应[10-12]，相关工作[13-14]表明风荷载与地震作用效应具有显著的耦合行为。然而，现有风力机气弹性软件的建模能力有限，难以模拟复杂基础结构和岩土体非线性行为；同时，由于风速场、地震动具有随机性，耦合方法计算量较大。这些困难推动了风力机解耦动力响应分析方法的发展。

气动-弹性解耦地震响应分析方法的应用要早于耦合方法。针对千瓦级风力机，叶轮-机舱(RNA)被简化为塔顶处的集中质量，气动力则被当做静荷载[15]。根据叶素动量(BEM)理论，Zhao和Maißer[16]将叶轮气动荷载分解为平均推力和脉动推力，后者可等效为气动阻尼。鉴于多兆瓦风力机叶片结构模型对系统的地震响应至关重要[17]，Asareh 等[18]通过梁单元模拟叶片，将叶轮气动阻尼力等效为塔架模态气动阻尼比，提出了一种风力机地震响应解耦分析方法，并开展NREL 5 MW风力机地震易损性分析。Santangelo 等[19]通过耦合、解耦方法计算结果的比较，评价了这种解耦方法的可靠性。解耦方法计算效率高，被广泛用于风力机地震响应分析。然而，该解耦方法将不同模态和风速的气动阻尼比取为相同值。这一假定与水平轴风力机气动阻尼研究的最新研究结论并不一致[20]。同时，Xi 等[21]选取一组包括近场、远场和近场脉冲地震动的强震记录，发现该方法不具有普适性，可能产生较大误差。

考虑到现有解耦方法的气动阻尼模型存在错误、其仿真结果也不满足工程精度要求，本文基于风力机模态气动阻尼理论模型，通过最小二乘法，简化气动阻尼比与平均风速关系，形成一种改进的风力机地震响应解耦分析方法，通过耦合、解耦方法计算结果的比较，验证该方法的可靠性，并评估解耦方法的计算效率。

1 气动阻尼简化模型

图1(a)为典型的三叶片水平轴风力机，是工业界应用最广泛的机型。图示的惯性参考系xyz是固定于地面的笛卡尔坐标系，x轴沿风力机前后向，指向风速场下风向，z轴竖直向上，y轴沿风力机侧向。

1.1 模态气动阻尼理论

根据文献[22]建议，不考虑叶片变形，依据叶素动量理论，半径r处叶素的气动荷载为：

叶轮旋转使得风力机成为时变体系，为建立系统气动阻尼模型，参照文献[22]，将叶片近似为刚体，塔架变形引起的叶片振动速度可表示为：

这些气动荷载增量与塔顶平动和转动速度大小成正比，方向相反，对风力机具有阻尼效应，因此，被称为气动阻尼力。Xi 等[20]将式(9)～式(10)所示的叶片分布气动阻力向塔顶截面形心处简化，并根据风力机前后向、侧向模态，得到水平轴风力机模态气动阻尼比，具体的推导过程可参照该文献。

NREL 5 MW 风力机是应用最广泛的参考机型，主要技术参数如表1 所示，更多的结构、控制系统参数可参考其研究报告[23]。本文采用欧拉梁单元模拟叶片和塔架动力行为，如图1(b)所示，将塔架等分为50 个梁单元。叶片划分为17 个梁单元，如图1(c)所示，长度依次为：自叶根起，前3 个2.73 m，中间11 个4.1 m，其余3 个2.73 m。该风力机前后向、侧向前二阶振型的模态气动阻尼比如图2 所示。对于前后向模态，风速小于额定风速时，气动阻尼比单调递增；风速超过额定风速后，气动阻尼比变化较小，呈单调递减趋势。对于侧向模态，气动阻尼比单调递增。同时，前后向和侧向二阶模态的气动阻尼比显著小于一阶模态。因此，现有解耦方法将各模态的气动阻尼比取为常数是不恰当的。

表1 NREL 5 MW 风机的主要性质Table 1 Properties of NREL 5 MW wind turbines

图2 NREL 5 MW 风力机模态气动阻尼比Fig. 2 Modal aerodynamic damping ratios of NREL 5 MW wind turbine

1.2 双线性气动阻尼模型

水平轴风力机的模态气动阻尼比与平均风速有关，呈非线性变化，不利于其在工程中的应用。因此，本节通过最小二乘法，对不同风速条件下的风力机模态气动阻尼比进行简化。

根据图2，对于前后向模态，平均风速介于切入风速、额定风速之间时，假定气动阻尼比呈线性变化；平均风速介于额定风速、切出风速之间时，假定气动阻尼比为常数。因此，风力机前后向一阶、二阶模态气动阻尼比均可表示为：

表2 拟合曲线参数Table 2 Parameters of fitting curves

图3 和图4 分别为NREL 5 MW 风力机前后向、侧向模态气动阻尼比。拟合值与理论解具有良好的一致性，因此，可采用这一简化模型描述气动阻尼比随轮毂高度处10 min 平均风速的变化规律。

图3 NREL 5 MW 风力机前后向模态气动阻尼比Fig. 3 Modal aerodynamic damping ratios of fore-aft modes for NREL 5 MW wind turbine

图4 NREL 5 MW 风力机侧向模态气动阻尼比Fig. 4 Modal aerodynamic damping ratios of side-side modes for NREL 5 MW wind turbine

2 风力机地震响应的解耦分析方法

考虑气动-弹性效应，半径r处叶素的法向气动荷载为：

根据式(15)～式(19)，叶片振动速度影响叶轮气动荷载，风-地震激励下的风力机动力响应属于典型的气弹性问题。现代风力机广泛采用了变速-变桨距控制技术[2]。因此，需要采用气动-伺服-弹性耦合模型分析风-地震共同激励下的风力机动力响应。以FAST 为例，该软件基于Kane 方法，建立的风力机系统运动控制方程为：

3 算例与分析

研究表明：现有解耦方法的计算精度与地震动有关[21]，因此，参照FEMA 建议[25]，选取一组包括近场、远场和脉冲地震动的强震记录作为输入地震动。考虑到本文旨在探讨气弹性解耦地震响应分析方法的可靠性，将这些原始记录作为输入地震动。

一个完整的强震记录包含两个水平分量和一个竖向分量。风力机的气动-弹性行为主要在其前后向，包括其他两个分量会低估解耦模型的误差。因此，对于所选强震记录，采用附录A 所列的水平分量为输入地震动，并假定地震动和风速方向一致，不考虑偏航，这是气弹性效应最为突出的工况。总模拟时长取600 s，时间步长取0.002 s，为避免初始条件的影响，地震动400 s 时开始输入。验证性算例表明计算参数的取值能够保证数值模拟的精度与稳定性。

3.1 稳态风-地震激励

鉴于1.1 节所述气动阻尼模型针对均匀稳态风，本节探讨了稳态风-地震动激励下新解耦方法分析风力机地震响应的可靠性。附表A1 所选的2 号和4 号地震动分量分别来自1999 年土耳其Duzce 地震的Bolu 记录和台湾Chi-Chi 地震的Chy101 记录，各自的加速度时程和反应谱如图5所示。其中，无量纲加速度反应谱值q通过式(23)获得：

图5 Bolu 强震记录的一个水平分量Fig. 5 One horizontal component of Bolu seismic record

将风速取为11.4 m/s，首先，以Bolu 波为输入地震动。稳态风、地震动单独作用引起的塔底弯矩时程如图6 所示。稳态风单独激励下，由于气弹性效应，塔底弯矩均值为60 MN·m，塔顶加速度峰值为0.01 m/s2，因此，其对塔顶加速度的贡献可以忽略。图7 为耦合、新解耦方法得到的塔顶加速度和塔底弯矩时程。整体而言，两种方法所得结果具有较好的一致性。耦合、新解耦方法得到的塔顶加速度峰值分别为2.21 m/s2和2.23 m/s2。对于塔底弯矩，耦合、新解耦方法得到的峰值分别为156 MN·m 和158 MN·m。相对误差不超过5 %，满足工程精度要求。

图6 Bolu 强震记录单独激励下的风力发电机塔底弯矩Fig. 6 Tower-base bending moment of wind turbine excited by Bolu seismic record only

图7 稳态风-Bolu 波共同激励下的风力发电机动力响应Fig. 7 Dynamic response of wind turbine excited by combined steady wind and Bolu wave

随后，以Chy101 波为输入地震动，风速仍取为11.4 m/s。地震动单独作用引起的塔底弯矩时程如图8 所示。图9 为耦合、新解耦方法得到的塔顶加速度和塔底弯矩时程。耦合方法和新解耦方法得到的塔顶加速度峰值分别为3.14 m/s2和3.29 m/s2。对于塔底弯矩，耦合方法和解耦方法得到的峰值分别为173 MN·m 和185 MN·m。相对误差小于10 %，满足工程精度要求。

图8 Chy101 强震记录单独激励下的风力发电机塔底弯矩Fig. 8 Tower-base bending moment of wind turbine excited by Chy101 seismic record only

图9 稳态风-Chy101 波激励下的风力发电机动力响应Fig. 9 Dynamic response of wind turbine excited by steady wind and Chy101 wave

考虑篇幅，本节仅列出上述两个算例。实际上，对于附表A1 的所有输入地震动，稳态风-地震激励下，新解耦方法分析风力机动力响应具有足够的精度，满足工程要求。

3.2 湍流风-地震激励

风具有显著的湍流特性，脉动风功率谱密度函数采用Kaimal 谱：

式中：Δr为两节点的距离；b为衰减系数，取为12；LC为湍流积分尺度，取为340.2 m。将湍流风速场总时长取为650 s，利用Turbsim 软件，生成风速场样本。

根据已有研究，轮毂高度处平均风速等于额定风速时，风力机受到的气动载荷最大。首先，输入地震动为Bolu 波，将轮毂高度处平均风速取为11.4 m/s，采用上述模型，生成5 个随机风速场，分别采用耦合、新解耦和现有解耦方法分析风力机地震响应，塔架响应峰值如图10 所示。由于生成了5 个风速场，此处的峰值是关于风速场样本的平均。对于塔架位移、加速度、剪力和弯矩峰值，新解耦方法与耦合方法的误差不超过10%，满足工程精度要求；而现有解耦方法具有显著的误差。由于脉动风的影响，湍流风速随时间变化。现有解耦方法采用的模态阻尼和瑞利阻尼模型均为时不变的，也为了保持解耦方法的易实施性，因此，根据轮毂高度处10 min 平均风速值确定模态气动阻尼比。

图10 Bolu 波-湍流风激励下的塔身响应峰值Fig. 10 Response amplitudes of the tower when the wind turbine excited by Bolu wave and turbulent wind

为涵整个运行状态，将轮毂高度处10 min 平均风速分别取为5 m/s、7 m/s、9 m/s、11.4 m/s、13 m/s、15 m/s、18 m/s、21 m/s 和24 m/s，图11为塔顶加速度和塔底弯矩幅值。新解耦模型与耦合模型的相对误差不超过10%，满足工程精度要求；现有解耦方法则存在较大误差。

图11 Bolu 波-风激励下风力发电机响应峰值Fig. 11 Response amplitude of wind turbines excited by Bolu seismic wave and wind

随后，输入地震动为Chy101 波，轮毂高度处平均风速等于11.4 m/s 的塔架响应峰值如图12 所示。新解耦方法得到的塔架位移、加速度、剪力和弯矩峰值与基准方法的误差不超过10%。图13为不同平均风速条件下的塔顶加速度和塔底弯矩峰值。新解耦模型与耦合模型的相对误差不超过10%，满足工程精度要求。对于Bolu 波和Chy101波，湍流风-地震动激励下，新解耦模型分析风力机动力响应的精度满足工程需要，且显著高于现有解耦方法。

图12 Chy101 波-湍流风激励下的塔身响应峰值Fig. 12 Response amplitude of the tower induced by Chy101 wave and turbulent wind

图13 Chy101 波-风激励下的风力发电机响应峰值Fig. 13 Response amplitude of wind turbines excited by Chy101 seismic wave and wind

图14 为叶片各单元法向气动荷载，分布特征与轮毂高度处平均风速有关。对于叶根起的前3 个单元，气动荷载较小，因此，叶轮气动合力受它们的影响较小。轮毂高度处10 min 平均风速等于切入风速(3 m/s)时，相对风速最小，此时，地震引起的叶片振动速度平方与相对风速平方之比η 最大：

图14 叶片各单元法向气动荷载Fig. 14 Aerodynamic forces of blade element in the normal direction

输入地震动为Bolu 波、Chy101 波的叶片振动速度平方比如图15 所示。对于前4 个节点，相对风速较小，因此，速度平方比η 较大。此时，尽管式(22)会产生一定的截断误差，然而，这些节点所在单元对叶轮气动荷载贡献小，因此，对于支撑结构动力响应的影响较小。对于其他地震动，地震引起的叶片振动速度平方与相对风速平方之比η 也具有这一特征。实际上，新解耦方法和耦合方法仿真结果的比较也间接表明了这一假定的合理性。

图15 叶片振动速度平方比Fig. 15 Ratio for squared blade vibration velocity

3.3 解耦模型误差分析

为评估解耦模型分析风力机地震响应的可靠性，支撑结构响应误差δ 定义为：

轮毂高度处平均风速等于11.4 m/s 时，塔顶加速度和塔底弯矩误差如图16 所示。对于所有地震动，风力机支撑结构响应峰值误差均小于15%，与Xi 等[21]评估的解耦方法相比，本文改进的解耦模型具有更高的精度。

图16 平均风速等于11.4 m/s 的解耦模型误差Fig. 16 Errors of uncoupled model when mean wind speed is 11.4 m/s

轮毂高度处平均风速等于5 m/s 和18 m/s 时，塔顶加速度和塔底弯矩误差分别如图17 和图18所示。对于所有地震动，支撑结构响应误差均小于15 %。实际上，平均风速取其他值时，新解耦模型预测风力机支撑结构响应峰值的误差也不超过15 %，显著小于Xi 等[20]评估的解耦方法。

图17 平均风速等于5 m/s 的解耦模型误差Fig. 17 Errors of uncoupled model when mean wind speed is 5 m/s

图18 平均风速等于18 m/s 的解耦模型误差Fig. 18 Errors of uncoupled model when mean wind speed is 18 m/s

3.4 计算效率的比较

GL 标准[27]要求将地震作用与所有可能出现的环境荷载组合。考虑到这些作用的随机性，风力机地震响应分析涉及多种荷载工况。以3.3 节的湍流风-地震动激励下风力机动力响应分析为例，涉及风速场样本45 个、强震记录50 个。耦合方法需要开展2 250 次时域地震响应分析，解耦方法仅需进行95 次时域分析。当然，解耦方法还需将风作用效应与地震作用效应时程线性叠加。整体而言，解耦方法的计算量小于耦合方法。对于本算例，使用的计算机置为：CPU Intel I7 7700k，内存16 GB。耦合方法消耗的CPU 时长为11 250 min，解耦方法消耗的时长为670 min。因此，解耦方法计算效率高于耦合方法。

4 结论与展望

本文基于水平轴风力机气动阻尼理论，通过最小二乘法将气动阻尼模型简化，改进风力机地震响应的解耦分析方法；随后，以NREL 5 MW风力机为例，探讨稳态风、湍流风分别与地震联合激励下风力机地震响应。可得以下结论：

(1) 线性化模态气动阻尼模型与风力机气动阻尼理论解具有良好的一致性，因此，可采用这一简化模型描述水平轴风力机模态气动阻尼比随风力机轮毂高度处10 min 平均风速的变化规律。

(2) 对于本文选取的所有输入地震动，解耦方法分析风-地震联合激励下风力机动力响应时，塔顶加速度和塔底弯矩峰值误差均不超过15 %，基本满足工程精度需要。

(3) 对于本文3.3 节的验证性算例，风力机耦合地震响应分析方法需开展的时域模拟数量超过解耦方法的20 倍，由于解耦方法需要在时域将风、地震作用效应叠加，耦合方法消耗的CPU 时长约解耦方法的16 倍。

为保证解耦、耦合两类模型的差异仅在于气动-弹性处理方法，本文的数值模拟工作均在FAST 软件中实施。实际上，这种解耦方法可用于商用有限元软件。需要指出，本文将风力机气动阻尼效应以模态气动阻尼施加于塔架，还可以通过塔顶集中阻尼器或者叶片分布阻尼器的方式实现，这些方法待进一步探讨。同时，NREL 发布的FAST 软件尚不能考虑土-结相互作用，后续工作应基于FAST 程序开发桩-土非线性相互作用模型，进一步验证上述解耦方法在考虑桩-土非线性相互作用后的适用性。

附录：