复合材料风力机叶片表面裂纹扩展的模拟研究

焦念鹏，宋 力，2，3，陈永艳，2，3，焦晓峰，冯 瑞，樊 亮

（1.内蒙古工业大学 能源与动力工程学院，内蒙古 呼和浩特 010051；2.风能太阳能利用技术教育部重点实验室，内蒙古 呼和浩特 010051；3.内蒙古自治区可再生能源重点实验室，内蒙古 呼和浩特 010051；4.内蒙古电力（集团）有限责任公司内蒙古电力科学研究院分公司，内蒙古 呼和浩特 010051；5.国水集团化德风电有限公司，内蒙古 乌兰察布 013350）

0 引言

裂纹是风力机叶片损伤的主要形式之一，风力机在运行过程中受到复杂载荷的影响，容易在应力集中部位产生裂纹。在应力的积累过程中，裂纹会产生进一步扩展，进而导致叶片的断裂。因此，研究叶片表面裂纹的扩展机理具有重要意义。

研究风力机叶片表面裂纹，首先要研究叶片的受力特性，通常基于流固耦合理论对叶片的受力进行研究。Himayat Ullah[1]计算了复杂载荷下风力机叶片的应力分布规律，找到了容易导致叶片失效的高应力区域。周勃[2]基于流固耦合理论，通过分析叶片表面裂纹的应力强度因子，研究其扩展规律。由于风力机叶片由复合材料制成，其表面裂纹扩展的机理极为复杂。近年来，随着扩展有限元（XFEM）方法的发展，其在复合材料裂纹扩展的研究领域得到了推广。杨宇宙[3]基于XFEM法分析了内冲击载荷循环作用下复合材料的疲劳裂纹扩展特性。彭英[4]针对某平板结构，研究了裂纹初始角度、初始长度和所受载荷大小等因素对裂纹扩展路径的影响。Matthias Holl[5]提出了一种研究三维物体表面裂纹扩展的多尺度方法，并将该方法应用于燃气轮机叶片表面裂纹的研究。胡舵[6]以具有初始裂纹损伤的复合材料压力容器为研究对象，通过试验和数值模拟相结合的方法，对复合材料裂纹扩展机理进行了分析。

综上所述，在风力机受力特性的研究领域，国内外专家学者已进行了深入的研究，但对于叶片表面裂纹扩展的研究尚不全面。因此，本文将XFEM理论引入叶片表面裂纹的研究之中。本文以1.5MW水平轴风力机为研究对象，基于流固耦合理论，首先分析叶片的受力特性，得到叶片的受力集中部位。进而通过XFEM法，分析叶片表面裂纹的扩展规律。

1 理论模型

1.1 流固耦合控制方程

流固耦合法是研究流场与流场中的固体之间相互作用的分析方法，在计算流体力学中，通常把流体和固体的控制方程直接耦合到同一个方程矩阵中，基于同一求解器同时求解流体的控制方程和固体的控制方程[7]。

式中：Aff为流场的系统矩阵；Asf为固体耦合矩阵；Afs为流体耦合矩阵；Ass为固体的系统矩阵；，分别为流体和固体的待求解量；k为迭代时间步Bf为流体外部作用力；Bs为固体外部作用力。

1.2 XFEM动态方程

扩展有限元的基本原理是单位分解法（PUM）。对于如图1所示的有限元网格中的任意一条裂纹，为了对其扩展过程进行更精确的求解，确保裂尖在扩展过程中仅在单元边界之间移动，需要引入阶跃增强函数[4]。

图1 含任意一条裂纹的网格Fig.1 Mesh containing any one crack

对 于 该 裂 纹，其 位 移 场uh（x）可 按 式（2）表述[8]。

式中：I为单元节点；NI（x）为节点的位移形函数；uI为节点的常规自由度；H（x）为广义Heaviside（强不连续）函数；K为单元内全部点的集合；Kr为被裂纹面切割的全部点的集合。

2 基于流固耦合的叶片受力分析

2.1 叶片模型的建立

本文以1.5MW大型水平轴风力机为研究对象，结合内蒙古西部某风电场提供的真实风力机模型设计叶片的气动外形。其中，单支叶片长度设为37.5m，轮毂直径为2m，轮毂中心距地面高度为65m，定义风力机的额定风速为12m/s。

基于Ansys软件的ACP模块，参考文献[9]，选择叶片材料为玻璃钢，对风力机叶片进行铺层设计。沿展向将叶片分成了5个区域，具体铺层结构如图2所示。

图2 铺层结构示意图Fig.2 Layer structure schematic

真实风力机叶片受到切变风的作用，其受力从叶根到叶尖依次递减，因此实际叶片的铺层厚度也是从叶根到叶尖依次递减，各层材料铺层角度如表1所示。

表1 叶片各区域铺层方式Table1Layout mode of each area of the blade

2.2 流场模型建立

建立长为300m，宽、高均为150m的计算域（图3）[9]，设置计算域入口为速度入口，出口为压力出口。在距入口65m处设置直径为80m，长为4m的圆柱形流场域。在流场域的中心设置风轮模型，轮毂的中心与旋转域中心重合。设轮毂中心距地面高度为65m，采用四面体网格的形式进行网格划分。

图3 计算域网格Fig.3 Computing domain mesh

2.3 边界条件设置与网格无关性验证

本文通过编译UDF函数模拟风切变来流变化，以此模拟受高度影响的实际风速的变化情况。切变来流的控制方程[10]为

式中：VZ为高度Z处的风速；Zref为轮毂中心所在高度；Vref为轮毂中心风速；α为风剪切系数，一般取0.2。

对流场模型的网格数量进行无关性验证。以风速为19m/s的工况为研究对象，分析叶片最大应力随网格数量的变化（表2）。当网格总数达到304万以上时，叶片表面应力值不再发生变化，故采用的网格划分单元总数为304万。本文基于滑移网格稳态计算法对旋转域进行设置，基于设计参数设定顺时针风轮旋转正向，风轮转数为19.8 r/min，计算时同时考虑了气动、重力和离心力载荷的影响。

表2 叶片最大应力随网格数量的变化Table2Stress at crack tip varies with the number of meshes

文章以风场所在的内蒙古地区12月份的气象参数为参考，设定来流密度为1.147kg/m3，湍流度为7%，空气温度为256.45K。同时模拟选用SSTk-ω湍流模型，算法设置为SIMPLE，差分方式为二阶迎风。设置残差为10-4，计算步数为2000步。

2.4 计算结果分析

风力机叶片在强风风速下更容易产生损伤，一般来说，强风指的是风速大于17m/s[11]。因此，本文选用19m/s的风速为计算工况（图4）。叶片主要在迎风面产生应力集中，最大值为23.86MPa。

图4 叶片应力分布云图Fig.4 Blade stress distribution nephogram

定义叶片长度为R，叶根至叶片任意翼型截面的距离为r，图5所示为叶片应力分布曲线。

图5 叶片应力分布曲线Fig.5 Blade stress distribution curve

由图4，5可知：在位置r/R=0.1和r/R=0.6附近存在应力集中，表明叶片表面的受力主要集中于叶根和叶中部位；叶片后缘的应力明显大于前缘，这是由于叶片后缘的相对厚度薄，更容易产生大的变形。因此，叶根和叶中（尤其是叶中）是裂纹的高发区域，且叶片后缘更容易产生裂纹。

3 叶片表面裂纹扩展分析

3.1 物理模型的建立

叶中部位的相对厚度较薄，位于该处的裂纹对叶片危害更高。截取叶片r/R=0.6截面位置处一长度为1.5m的片段。因叶片中部位于铺层区域4，参照表1对各层纤维方向进行设置。XFEM网格与结构内部物理界面无关，因此可以克服裂纹扩展时裂尖变形带来的网格划分困难问题。故本文基于XFEM进行研究，在叶片后缘建立长250mm，深10mm的 裂 纹（图6）。

图6 裂纹与约束设置示意图Fig.6 Crack and constraint setting schematic

由于叶片的主要变形方向为展向，故可以将叶片简化为梁结构，其变形可看成是由单一展向拉应力导致的[12]。在叶片段的一端添加固定约束，在另一端添加拉应力载荷约束。

3.2 角度对裂纹扩展规律的影响

选用各向异性均质材料和复合材料的两种叶片段进行分析研究，两种材料叶片段的大小以及表面裂纹的尺寸和设定位置保持一致。如图7所示，在翼型段表面分别设置裂纹角度为30°，60°和90°，研究裂纹角度变化对裂纹扩展路径的影响。

图7 裂纹角度设置Fig.7 Crack angle setting

因在2.4中计算了风速为19m/s下的无裂纹叶片的受力情况，故沿叶展方向提取距r/R=0.6截面0.75m处截面的平均应力（图8），作为施加载荷的依据。基于此，在叶片段右侧施加30.05 MPa的拉应力。

图8 应力提取位置Fig.8 Location of stress extraction

图9所示为叶片段的应力分布随裂纹角度的变化。

图9 叶片段的应力分布随裂纹角度的变化Fig.9 The stress distribution of blade segment varies with crack angle

在 拉 力 载 荷 的 作 用 下，30°，60°和90°的 裂纹均会沿叶片弦线方向进行扩展，也就是按典型的Ⅰ型裂纹方式进行扩展。从叶片段表面应力值来进行分析：对于均质材料结构，当裂纹角度从30°变化到90°，表面应力值从90.03MPa增加到127.50MPa；对于复合材料结构，表面应力值从55.66MPa增加到65.34MPa。裂纹角度越靠近90°，扩展趋势越强，扩展时裂尖应力集中现象越明显。

图10为裂纹生长量随时间步长的变化曲线。

图10 裂纹生长量随时间步长的变化Fig.10 Crack growth changes along the length direction

由图10可知：沿90°分布的裂纹，在时间步长为30时开始扩展；沿60°分布的裂纹，扩展时时间步长为20；沿30°分布的裂纹，在时间步长为10时就已产生扩展。角度为90°的裂纹，其尖端距离右侧加载端最远，起裂时间最长。但由于Ⅰ型裂纹的性质，结合图9可知，受到拉伸载荷作用时，角度越靠近90°的裂纹，其尖端越容易产生应力集中。因此，相同时间内，90°裂纹的扩展速率最快。以90°裂纹为例，均质和复合材料叶片裂纹在长度上的最大生长量分别为298.74mm和458.204 mm，但在深度上，均质叶片裂纹的生长量在时间步长为100时为26.42mm，复合材料叶片裂纹则几乎不存在扩展。结合云图分析，证明了裂纹在叶片内部的扩展更容易导致应力集中，复合材料的铺层结构可以抑制裂纹对叶片内部结构的破坏。

3.3 尺寸对裂纹扩展规律的影响

选取不同拉力载荷下，复合材料表面不同尺寸的裂纹为研究对象。

保持裂纹深度为10mm，改变其长度[图11（a）]，以30MPa拉 力 载 荷 工 况 为 例，时 间 步 长 为100时，随裂纹的初始长度从200mm增加到350 mm，其起裂应力从52.50MPa下降到32.42MPa，最大长度生长量由397.134mm增加到521.078 mm。随裂纹的初始长度增加，其起裂时所需应力减小，扩展速率提高。

图11 裂纹特征值随初始尺寸的变化Fig.11 The variation of crack characteristic values with the initial size

保持裂纹长度为250mm，改变其深度[图11（b）]，随裂纹的初始深度增大，其扩展需要破坏的纤维层数随之增加，起裂所需应力随之增大。在30MPa拉力载荷下，时间步长为100时，初始深度为15mm和20mm的裂纹，其起裂应力从33.34MPa增加到78.82MPa。同时，裂纹的扩展速率随其深度增加呈下降趋势。裂纹扩展时，裂尖应力值与裂纹尺寸成正比，裂纹尺寸越大，对叶片的危害程度越高。

4 结论

本文以风电场实测数据为基础，研究1.5MW风力机的应力分布，找到了叶片的应力集中区域。从该区域截取一片段，在片段上预设裂纹，对叶片表面裂纹的扩展进行分析，发现了裂纹扩展随叶片材料、裂纹角度、尺寸等因素变化的规律，得到以下结论。

①叶中和叶根部位，尤其是叶中，容易产生应力集中。风速为19m/s以下时，叶中部位应力最大，其值为23.86MPa。可判断裂纹主要集中于叶根和叶中部位，且叶片后缘受力更大，因此叶片后缘产生裂纹的概率远大于前缘。

②对于均质材料叶片表面裂纹，其在受力时会同时沿长度和深度方向进行扩展。对于复合材料叶片表面裂纹，其在受力时仅沿长度方向扩展。90°裂纹的扩展速率最快，时间步长为100时，均质和复合材料叶片表面90°裂纹扩展时裂尖的最大应力分别为127.50MPa和65.34MPa。复合材料的铺层结构能降低裂纹扩展时裂尖的应力集中，阻碍裂纹扩展对叶片的破坏。裂纹在叶片内部扩展时，裂纹更易集中应力，对叶片的损害更大。

③裂纹扩展时的路径基本沿叶片弦向，其裂尖处应力值大小与裂纹初始尺寸成正比。裂纹初始长度越大，起裂所需应力越小，扩展速率越快。但裂纹的深度越大，起裂所需应力越大，扩展速率越慢。