目前, 随着风力机容量的不断增大, 塔架在风载荷作用下的可靠性愈加受到国内外学者们的关注。刘雄等人^[1]分析了塔架的气动阻尼和结构阻尼, 并将塔架简化成悬臂结构, 得到了其在变载荷作用下的动态响应。闫海津等人^[2]考虑塔架效应, 利用FLUENT软件研究了风力机塔架附近的流场。张湘伟和文武^[3]研究了大型风力机塔架在脉动风下的动力响应特性。刘新喜等人^[4]基于ABAQUS有限元软件, 建立了4种不同塔架形式的风力机一体化模型, 探讨了不同塔架形式风力机的抗台风性能。王介龙等人^[5]研究了风机不同部件之间的耦合问题, 建立了桨叶、机舱、塔架的动力学方程, 并利用子空间迭代方法计算了结构的动力响应。ZHANG J P等人^[6-7]在考虑流固耦合作用的情况下, 讨论了不同来流风速及湍流效应对风力机叶片动力响应的影响。

本文针对海上风力机建立了塔架结构的三维实体模型, 基于ANSYS开展了数值模拟, 分析了不同平均风速作用下考虑几何非线性的塔架动力响应, 以期为优化风机结构提供参考。

从实际的工程计算来看, 一般通过风压来表示风对结构的作用。由伯努利方程可知, 风速在单位面积上的风压F_e(t)为

$ F_{\mathrm{e}}(t)=\frac{1}{2} \rho v^{2} $

(1)

式中:ρ——空气密度;

v——平均风速。

由式(1)可得到风载荷计算公式^[8]为

$ F(t)=\frac{1}{2} C_{\mathrm{P}} \rho A v^{2} $

(2)

式中:C_P——风压分布系数;

A——作用面积。

几何非线性是指放弃小位移假设, 从几何上严格分析单元体的尺寸和形状变化, 得到非线性的几何运动方程, 由此引起基本控制方程的非线性问题。对于结构的几何非线性分析, 可以归结为节点位移的非线性方程。在有限元分析中, 根据系统的几何非线性特性对系统的平衡方程进行不断的修正, 在新的一步增量求解前, 对坐标系进行修正, 然后去求解方程, 并计算几何非线性对刚度矩阵和载荷阵的修正。有限元方法主要在单元体离散化的基础上实现微分方程的求解。一般情况下, 塔架结构所满足的几何非线性离散运动微分方程为

$[M][\ddot{x}]_{t}+[C][\dot{x}]_{t}+[\bar{P}]_{t}=[F(t)]$

(3)

式中:[M], [C]——质量矩阵和阻尼矩阵;

几何非线性项, 与刚度矩阵相关;

结构的速度和加速度;

[F(t)]——塔架受到的风载荷列向量。

此外, 在计算时不考虑阻尼, 即阻尼矩阵[C]=[0]。

利用Newmark法, 推导结构在风载荷作用下的几何非线性动力增量平衡方程。在有限元分析中, 结构控制方程对应于某一时刻t+Δt的平衡方程, 可表示为

$[F(t)]_{t+\Delta t}-[M][\ddot{x}(t)]_{t+\Delta t}-[\bar{P}]_{t+\Delta t}=0$

(4)

式中:Δt——时间步长。

式(4)的几何非线性项可线性表示为

$[\bar{P}]_{t+\Delta t}=[\bar{P}]_{t}+\left[K_{\mathrm{T}}\right]_{t}\left([x]_{t+\Delta t}-[x]_{t}\right)$

(5)

式中:[K_T]_t——t时刻的节点位移算出的切向刚度矩阵;

[x]_t——结构位移, [x]=[u, vw]^T。

将式(5)代入式(4), 可得到

$[M][\ddot{x}]_{t+\Delta t}+[\bar{P}]_{t}+$ $\left[K_{\mathrm{T}}\right]_{t}\left([x]_{t+\Delta t}-[x]_{t}\right)=[F(t)]_{t+\Delta t}$

(6)

令

$[x]_{t+\Delta t}-[x]_{t}=\Delta[x]$

(7)

通过结构的初始位移、速度和加速度计算出位移增量Δ[x], 利用Newton-Raphson迭代法进行迭代计算来解决由线性化带来的误差, 最终计算出结构新的位移、速度和加速度, 即可得到结构整个动力响应的时程。

根据文献[9], 用于数值模拟的海上风力机塔架几何参数如表 1所示。

海上风力机塔架材料通常使用的是Q345钢材。本文中计算模拟的塔架材料Q345钢的基本材料参数如表 2所示。

根据表 1和表 2给出的塔架参数, 采用SolidWorks软件进行三维建模, 并将建立的塔架实体模型导入到ANSYS中, 所得模型如图 1所示。

塔架实体模型建立后, 划分网格成为一个重要问题。依次增大网格密度, 结果发现, 塔架顶端最大位移随着网格单元数的增加是收敛的, 如表 3所示。

为了保证计算精度, 同时节约计算时间, 这里采用23 265个网格单元数进行数值计算分析, 网格划分后的模型如图 2所示。

基于动力学分析方法, 结合ANSYS Workbench软件对风力机塔架进行振动特性分析, 选定8 m/s, 12 m/s, 16 m/s 3种平均风速(这里的平均风速指的是风的大小和方向都不随时间变化), 由此开展塔架动态特性分析, 得出塔架的位移和Mises应力的动态特性。

图 3为3种平均风速工况下, 塔架顶部在风载荷作用方向上的最大位移响应曲线。

由图 3可知, 其中, 8 m/s工况下塔架的最大位移为1.51×10^-3 m, 12 m/s, 16 m/s工况下的最大位移分别为3.39×10^-3 m和4.62×10^-3 m, 这表明随着平均风速的增加, 风力机塔架的位移也随之增大。

图 4给出了塔架最大位移随平均风速变化的曲线。考虑风速范围为8~45 m/s, 由图 4可以看出, 最大位移随着平均风速的增大而增大, 呈非线性递增关系。

图 5至图 7分别为3种平均风速工况下整个风力机塔架最大Mises应力时间响应曲线。

3种工况下, 塔架承受的最大Mises应力由8 m/s工况下的0.74 MPa增加到12 m/s工况下的1.65 MPa, 16 m/s工况下的最大Mises应力是8 m/s工况下的3.04倍, 表明随着风速的增加, 塔架最大Mises应力也随之而增加。

与图 4工况相同, 给出了近海风力机塔架平均风速作用下最大Mises应力与平均风速的关系, 如图 8所示。由图 8可以看出, 最大Mises应力随着平均风速的增大而增大, 基本呈非线性递增关系。

(1) 针对海上风力机塔架, 建立了考虑几何非线性的离散运动微分方程, 并基于Newmark和Newton-Raphson法, 实现了对平均风速作用下塔架位移和Mises应力响应的求解。

(2) 从3种不同平均风速下风力机塔架的位移、Mises应力响应曲线可以发现, 塔架最大位移随平均风速的增加而增加, 最大Mises应力也是如此。

(3) 考虑风速范围为8~45 m/s, 通过分析塔架最大位移、Mises应力随平均风速变化的曲线可以看出, 塔架最大位移和Mises应力均随平均风速的增大而增大, 且呈非线性递增关系。