东海某风场4 MW风机塔筒结构疲劳特性

2022-11-01和庆冬杨逢杰余泉苏凯

科学技术与工程 2022年26期

和庆冬，杨逢杰，余泉，苏凯*

(1.国家电投集团江苏海上风力发电有限公司，盐城 224005； 2.武汉大学水资源与水电工程科学国家重点实验室，武汉 430072； 3.浙江华东工程咨询有限公司，杭州 310014)

针对传统化石燃料消耗带来的气候变化、能源短缺、环境破坏等问题，大力开发可再生能源的呼声不断高涨。风能以其清洁可再生、成本相对较低、运行可靠安全、技术较成熟等优点，逐渐被世界各国重视，拥有广阔的发展前景。风电机组塔筒是风力发电机中连接桩基础与机舱的管状承载结构，不仅支撑着机舱、轮毂、叶片等上部结构，还需要承受外界各种荷载的作用，并为相关人员提供安装、维修等工作的平台。受地表环境、施工条件等方面因素的限制，钢筋混凝土塔筒无法广泛地运用到海上风电场中，因此现阶段海上风电场绝大部分塔筒结构都采用锥台型钢制塔筒，横截面为薄圆环。

海上风力机所处的自然条件与陆地风力机相比更为恶劣，其在整个运行期内将受到大量循环荷载的作用，这些交变荷载会使风力机塔筒产生振动，应力响应也在不断发生变化，进而造成结构的疲劳损伤，缩短风力机组的使用寿命，因此对于塔筒的疲劳研究十分重要。

21世纪以来，随着大型有限元软件的应用和数据监测技术的提高，风机结构疲劳分析方法发展迅速，针对塔筒结构疲劳方面的研究也取得了一定的成果。赵艳[1]根据现场测量的风机振动特征数据和风力数据，提出了基于实测数据的风机塔架疲劳寿命计算方法。霍涛[2]基于气象站观测数据，定义若干个风工况，研究了风速风向对塔筒结构疲劳寿命的影响。张广隶[3]基于主S-N曲线法(主疲劳应力范围-疲劳寿命曲线法)，评估并对比了三个塔筒模型风致疲劳寿命。刘文等[4]在静力分析的基础上，对塔筒两个关键点进行疲劳分析。程良灿[5]对比分析了塔筒-叶片整体模型和只建塔筒的简化模型两者的疲劳性能，并根据中国、美国、欧洲钢结构设计规范要求计算疲劳损伤。陈旺[6]建立了塔架-地基-叶片模型，计算了塔架在各种荷载作用下的疲劳损伤。柯世堂等[7]建立了塔架-叶片耦合模型，对风力机全机关键部位的风振疲劳寿命进行估算。Stavridou等[8]针对两个高度相同、壁厚有明显差异的塔筒，采用详细时程加载和解析公式评估两种方法计算疲劳寿命，结果表明厚度是决定塔筒抗疲劳性能的重要因素。欧阳儒贤[9]分析了螺栓断裂前后两种状况下风机塔筒的疲劳寿命，发现螺栓断裂后疲劳损伤急剧增加，寿命严重缩短。

综合以上研究成果可以发现，风机疲劳分析较为普遍的做法是：在确定结构危险部位的基础上，划分不同的工况对结构的疲劳损伤进行计算。然而，研究者在对塔筒结构进行建模时过于简化，通常默认塔筒的直径、壁厚是沿高度均匀变化的，但是实际工程中并非如此，每个筒节都有各自的直径和壁厚，而且两种建模方式之间的差别鲜有学者开展深入对比。鉴于此，现依托江苏某海上风电项目，建立塔筒三维精细化有限元模型，采用中国常用的疲劳规范[10]，从整体上详细讲述风力机塔筒疲劳研究的具体方法和步骤，并与壁厚随高度均匀变化的简化塔筒模型进行对比，以期为实际工程中塔筒结构设计和抗疲劳防护提供指导，为相关研究人员提供一定的借鉴和参考。

1 研究方法及有限元计算模型建立

1.1 研究方法

以江苏东海某风场4 MW风电机组塔筒结构为例，建立三维精细化有限元模型，通过静力分析获取塔筒结构若干个控制点；其次，利用模拟得到的每个风速工况的荷载时程，对各风速工况下的塔筒结构展开动力分析，获取各个控制点的应力时程；最后，结合雨流计数法和Goodman直线方程，根据相关规范选取S-N曲线和疲劳累积损伤法则，计算各个控制点的疲劳总损伤和塔筒结构的疲劳寿命，评估塔筒是否满足设计寿命要求。

上述疲劳分析的方法流程可用图1来表示。

图1 疲劳分析技术路线Fig.1 Technical route of fatigue analysis

1.2 结构基本参数及建模方式

本风场风力发电机为三叶片水平轴式，轮毂高度为95 m，塔筒总长为75.99 m，由3个筒段组成，每个筒段由若干个筒节组成，各筒段之间由法兰螺栓连接，各筒节之间通过焊接方式连接。塔筒顶部直径3 120 mm，底部直径5 500 mm，不同筒节的壁厚不同，最大壁厚65 mm，最小壁厚17 mm。塔筒结构采用Q355ND高强度低碳钢，弹性模量、泊松比、密度、抗拉强度分别为210 GPa、0.3、7 850 kg/m3、490 MPa。机舱、轮毂、单个叶片的质量分别为142、56.2、19.6 t。

在建立风力机塔筒结构时，忽略筒段之间的连接法兰、螺栓、塔筒底部门洞以及塔筒内部的附属结构等构件，也不考虑焊缝可能存在缺陷造成的影响，将塔筒视作统一整体。塔筒采用四节点壳单元模拟；机舱、轮毂、叶片简化成一个集中质量单元进行模拟，该质量单元处于轮毂高度上的塔筒横截面中心位置，并与塔筒顶部节点建立刚性连接以传递荷载。塔筒有限元模型采用笛卡儿直角坐标系，规定Z轴为竖直方向，以向上为正；X轴为顺风向方向；Y轴为垂直风向方向；坐标原点位于塔筒底部的中心位置。在塔筒迎风面节点上(即X<0的节点)建立表面效应单元，便于施加单方向风荷载。塔筒底部施加全约束。

根据本项目塔筒各个筒节的实际尺寸，所建立的塔筒三维精细化有限元模型如图2所示。图2中左侧放大部位为塔筒顶部的刚性连接区域；右侧三种不同颜色分别代表塔筒的三个筒段范围，小矩形的高度和宽度分别表示塔筒各个筒节的长度L和壁厚t，具体数值均已给出。

图2 塔筒精细化有限元模型图Fig.2 Refined finite element model of tubular tower

1.3 塔筒荷载分析

海上风电机组塔筒结构所经受的荷载主要包括：自重，机舱、轮毂、叶片等上部结构的重量，沿塔筒高度方向的风压荷载和风机荷载。

沿塔筒高度方向的风压荷载与平均风速有关，平均风速与高度的关系可以根据式(1)得到，风压荷载与风速之间的转换式如式(2)所示。

(1)

式(1)中：v1、v2分别为h1、h2高度上的风速；α为风切变指数。

(2)

风机荷载是风轮叶片在正常运转时空气对叶片的作用力，此处作简化处理，仅考虑风压荷载的影响，即风机荷载为作用在叶片迎风面积上的均布风压产生的顺风向力，风机荷载计算公式为

F=W0S

(3)

式(3)中：F为风机荷载；S为叶片迎风总表面积。

由于风力机起动、关机、紧急关机等工况发生的时间相对较短，所造成的损伤仅占总疲劳损伤的很小一部分[2, 11]，不到1%，因此在本次塔筒结构疲劳分析中不考虑以上工况产生的瞬时荷载。

建立好有限元模型并施加约束和荷载之后，通过静力分析可以获得塔筒模型中等效应力较大的10个控制点，即等效应力较大的10个节点。上述控制点(节点)的位置和编号如图3所示。由图3可以看出，以上控制点基本都位于塔筒第2筒段背风侧，而且都在筒节之间的连接处或附近。

图3 精细化塔筒模型控制点位置示意图Fig.3 Schematic diagram of the location of the control points of the refined tubular tower model

2 风速工况划分

编制塔筒疲劳应力谱是结构优化设计和疲劳寿命评估的基础[12]，对疲劳分析结果的精度也有大影响[13]。因此，准确获取塔筒随机循环应力，编制合适的疲劳应力谱，是分析结构寿命的重要环节。

诸如参数法等传统的疲劳应力谱外推方法，需要采取实测或仿真方法来获取大量的样本数据，在此基础上进行统计外推[14]。然而，由于风力机所受到的荷载大小主要与风速有关，而每个风速在一年中出现的概率有所不同，在采集样本的过程中难免会受到外部环境、主观因素等条件的限制[15]，不能兼顾各种风速条件下产生的荷载，导致无法全面客观地反映风力机塔筒结构实际所承受的荷载情况。

因此，风力机疲劳应力谱的编制不应该按照传统方法进行，需要在整个风力机工作范围划分成若干个风速段，形成对应数目的风速工况，分别针对每个风速工况下产生的应力进行统计处理[16]，最后按照各个工况出现的概率叠加得到风力机一年内的疲劳应力谱。工程中一般采用双参数威布尔分布来描述一个地区各风速段出现的频率，其概率密度函数和分布函数如下：

(4)

fw(v)=1-e-(v/c)k

(5)

式中：pw(v)为威布尔分布的概率密度函数；v为风速；k为形状参数；c为尺度参数；fw(v)为威布尔分布的分布函数。

根据当地气象站实测资料可以推得威布尔分布两个参数分别为k=2.0，c=9.2 m/s。在本项目风力机3～25 m/s的工作范围内，以2 m/s的区间将整个风速范围划分成11个风速工况，并以每个风速工况的区间中值来表征该工况的平均风速，各风速工况一年内出现的概率、时长如表1所示。

表1 各风速工况一年内出现的概率、时长Table 1 Probability and duration of occurrence of each wind speed condition within one year

其中某个风速工况一年内出现的小时数可由式(6)计算而得。

Ti=24×365Pi

(6)

式(6)中：Ti、Pi分别为风速工况i一年内出现的小时数、发生概率。

为便于后续的动力分析计算，需要编写程序获得空间中某一点的瞬时风速时程曲线。本文瞬时风速模拟选择线性滤波法中的自回归模型，脉动风速功率谱采用Davenport谱，每一种工况的模拟时长均为600 s。借助MATLAB软件模拟得到风速工况10下的瞬时风速时程曲线如图4所示，限于篇幅限制，其他风速工况模拟时程曲线不再赘述。

图4 风速工况10下瞬时风速Fig.4 Instantaneous wind speed under wind speed condition 10

3 塔筒结构疲劳分析

3.1 应力循环计数

结合各个风速工况下的瞬时风速时程曲线和塔筒的受力模型，在有限元软件中展开塔筒结构的瞬态分析，得到各控制点在不同风速工况下的应力响应时程。图5显示了风速工况10时1737号节点的应力时程(1737号节点位于第2筒段第1筒节和第2筒节交界处附近，高程约31.46 m)。

由图5可以发现，该控制点的应力时程曲线中均包含多次应力循环，本文选择在疲劳分析领域中运用最广泛的雨流计数法来提取每个循环应力的均值和幅值[17]。

图5 风速工况10时1737号节点等效应力时程Fig.5 Time history of equivalent stress of No.1737 node under wind speed condition 10

雨流计数法的计数过程与材料真实应力-应变特性相符[18]。雨流计数法的示意图如图6所示，其主要步骤如下。

图6 雨流计数法Fig.6 Rain-flow counting method

(1)根据获取到的控制点应力时程曲线，以绝对值最大的点为起点，对其进行重新排序。

(2)雨流依次从每个峰谷值开始流动，即从1、2、3等点位开始。

(3)雨流流到峰(谷)值处竖直下滴，一直流到对面出现比起始值更高的峰值(或更低的谷值)为止。

(4)当雨流遇到从上方流下的雨时便停止流动，构成一个应力循环。

(5)根据雨滴流动的起点和终点，逐一取出所有的应力循环，记录其峰值和谷值。

(6)每个循环峰值与谷值之差的一半作为该循环的幅值，峰值与谷值的平均作为该循环的均值。

在MATLAB中导入塔筒瞬态分析得到的应力结果进行雨流计数，可计算得到在仿真时间600 s内不同工况下塔筒任一控制点的雨流矩阵图。图7显示了风速工况10时1737号节点的雨流矩阵图。

图7 风速工况10时1737号节点雨流矩阵图Fig.7 Rain-flow matrix diagram of node 1737 under wind speed condition 10

运用式(7)～式(9)将仿真时间600 s内每一个风速工况下的均幅值的统计结果扩展至一年当中，然后进行叠加，可以得到任一控制点一年的雨流矩阵图。图8显示了1737号节点一年内的雨流矩阵图。

(7)

T′i=3 600Ti

(8)

(9)

式中：Nj为第j个均幅值区域一年出现的频次；T′i为风速工况i一年内出现的秒数；fij为第i个风速工况下第j个均幅值区域出现的频率；nij为第i个风速工况下第j个均幅值区域仿真时间内出现的频次；ti为第i个风速工况的仿真时间。

由图8可以看出，1737号节点一年内循环应力的均值大多集中在10～30 MPa范围内，幅值绝大部分集中在0～30 MPa范围内，远小于材料的屈服强度355 MPa，而且总循环次数超过106次，因此风力机在风荷载作用下通常发生高周疲劳。

图8 1737号节点一年雨流矩阵图Fig.8 Rain-flow matrix diagram of node 1737 in a year

运用Goodman直线方程(10)将非零均值的循环应力转化成零均值的循环应力，并统计各个幅值的超值累计频次，绘制成控制点的疲劳应力谱如图9所示。

图9 1737号节点应力谱Fig.9 Stress spectrum of node 1737

(10)

式(10)中：σa表示均值不为0情况下的幅值；σN(R=-1)表示均值为0情况下的幅值，即疲劳强度；σm为均值；σu为材料的抗拉强度。

3.2 疲劳计算结果

塔筒相邻的筒节之间是通过焊接连接的。对于具有不同名义厚度、直径的圆管对接焊缝，其本身和周围母材金属可能会产生应力集中。在本项目中，塔筒制造商生产检测工艺完备，基本不会造成筒节间椭圆度不一致和偏心，且相邻筒节之间焊接采用中心对齐，厚度过渡处进行了良好的削坡工作，因此本文在进行疲劳计算时没有考虑塔筒的应力集中效应。

海上风机塔筒结构焊缝及其周边的母材金属属于非管节点，长期暴露在空气中。根据中国船级社《海洋工程结构物疲劳强度评估指南》的推荐做法[10]，本项目塔筒结构疲劳分析的S-N曲线采用非管节点在空气中的E曲线。该S-N曲线在双对数坐标下可分为两段直线，当材料失效时的应力循环次数N≤107时，N和应力幅σa的关系可表示为式(11)。当N≥107时，N和σa的关系可表示为式(12)。

lgN=lg(1.04×1012)-3.0lg(2σa)

(11)

lgN=lg(2.30×1015)-5.0lg(2σa)

(12)

式中：N为在应力幅σa作用下产生疲劳破坏的循环次数；σa为应力幅，即一个循环中最大应力与最小应力之差的一半。

在工程设计中通常规定：钢制试件常温下经历107次应力循环仍未发生疲劳破坏时的最大应力幅值规定为钢材的疲劳极限，即只要应力水平不超过疲劳极限，就可以认为钢材的疲劳寿命趋近于无限次数，不会产生疲劳损伤。由式(11)和式(12)可以计算得107次对应的应力幅为23.51 MPa，即塔筒材料的疲劳极限为23.51 MPa。

在考虑筒壁厚度修正的基础上，根据Miner法则，塔筒各控制点在设计寿命20 a内的疲劳总损伤可由式(13)计算。各控制点的疲劳寿命可由式(14)计算。

(13)

(14)

式中：D为设计寿命内的疲劳总损伤；L0为设计寿命；k为应力水平总个数；ni为一年内在第i个应力水平σi下的实际应力循环数；Ni为在第i个应力水平σi下材料失效时的循环周次；Sftg为疲劳强度安全系数；L为疲劳寿命。

根据参考文献[10]可取Sftg=2.0。塔筒各个控制点的疲劳总损伤及疲劳寿命如表2所示。

由表2可以看出，以上10个控制点的疲劳总损伤D均小于0.5，满足疲劳设计要求，所以塔筒在设计寿命内不会产生疲劳破坏。其中2088号节点的疲劳寿命最短，为40.90 a，即整个塔筒结构的疲劳寿命为40.90 a，大于其设计寿命20 a要求。

表2 精细化塔筒模型控制点疲劳计算结果Table 2 Fatigue calculation results of control points of refined tubular tower model

4 塔筒有限元计算模型对比

相比于塔筒精细化有限元模型，直径或壁厚随高度均匀改变的简化塔筒模型具有建模方法简单、速度快的优势，被广泛采用。本文在上述研究的基础上，探讨简化塔筒模型对疲劳分析结果的影响。

本文的简化塔筒模型中，筒壁厚度随高度均匀变化，底部壁厚65 mm，顶部壁厚17 mm，其他条件与精细化塔筒模型保持一致，所建立的简化塔筒模型如图10所示。

图10 塔筒简化有限元模型图Fig.10 Simplified finite element model of tubular tower

对简化塔筒模型开展静力分析，得到等效应力较大的10个控制点位置如图11所示。由图11可以看出，简化塔筒模型控制点的位置与精细化塔筒模型完全不同，主要位于塔筒底部附近。简化塔筒模型相当于一根竖直向悬臂梁结构，靠近固定端处的应力比较大，且最大等效应力只有精细化塔筒模型的一半，因此塔筒有限元计算模型的选择对静力计算结果影响较大。

图11 简化塔筒模型控制点位置示意图Fig.11 Schematic diagram of the location of the control points of the simplified tubular tower model

按照疲劳分析步骤可计算得到简化塔筒模型各个控制点20 a的疲劳总损伤如表3所示。

对比表2和表3的疲劳计算结果可知，精细化塔筒模型控制点20 a的最大疲劳总损伤值为0.244 5，简化塔筒模型为0.001 0，两者相差了接近245倍。

表3 简化塔筒模型控制点疲劳计算结果Table 3 Fatigue calculation results of control points of simplified tubular tower model

综上所述，两种塔筒有限元模型在开展疲劳分析时所选取的控制点位置有很大的区别，而且最大总损伤值相差240多倍。简化塔筒模型的壁厚均匀过渡，减小了筒节之间壁厚突变可能产生的应力集中效应，是造成计算结果巨大差异的主要原因。因此，采用壁厚均匀变化的简化塔筒模型会使疲劳计算结果偏危险。在对风机塔筒进行受力分析和疲劳分析时，建议采用实际工程的塔筒结构尺寸(厚度)进行精细化建模计算。