徐 越 王跃方 李 聪 魏学敏 李盛文

（1.中国航空工业空气动力研究院；2.低速高雷诺数气动力航空科技重点实验室；3.大连理工大学）

0 引言

当前透平机械朝着高端化、节能、高可靠性的方向发展，远程运行监测、智能化控制及故障自愈化成为未来发展的趋势。叶片是工业轴流式压缩机的核心部件，承担功能转换的重要任务，运转时既承受离心力引起的高平均应力，又承受各种激振载荷引起的交变应力。当转-静子剐蹭或高阶模态被激振时，叶片边缘材料的初始缺陷易萌生微小裂纹，在转动过程中不断发展，形成宏观裂纹直至使叶片断裂。为保证旋转叶片长周期使役安全，研究初始裂纹萌生后的裂纹尖端应力场，判断应力强度因子是否超出材料断裂韧性而导致叶片断裂，对轴流叶片设计具有重要参考价值。

化工及空分装置中使用的轴流式压缩机叶片一般为薄壁金属件，离心力引起的拉应力在叶片横截面上的分布较为均匀，叶片受力接近于薄板平面应力状态。因此，可参照均质平板断裂力学分析方法计算叶片的应力强度因子[1]。对于含中心裂纹、处于平面应力状态的线弹性矩形薄板，中外学者大多采用权函数法推导裂纹尖端的I型应力强度因子理论解[2-5]。近年来，樊俊铃等[6-8]基于权函数法研究了平板含单边裂纹、双边裂纹和中心裂纹的应力强度因子，给出其与裂纹长度、位置、转速等参数之间的关系。用权函数法计算应力强度因子需要工程人员具备一定的断裂力学基础和解析推导能力。在数值解方面，尹奇志等[9]用有限元法计算了含裂纹平板的应力强度因子，表明采用仿真技术求解此类问题具有较好的精度。姜伟等[10]将动态有限元分析和相互作用积分相结合，计算了含裂纹平板受冲击载荷时的动态应力强度因子。大量研究表明，扩展有限元法(XFEM)在分析裂纹等不连续性问题给出了理想的计算结果。李录贤等[11]详细介绍了该方法并提供了应用实例，对以往数值结果作了修正。董玉文、余天堂等[12-13]研究表明扩展有限元法可以准确地求解应力强度因子。尽管扩展有限元法在求解裂纹扩展问题具有较大优势，但对于建模复杂、工况繁多等问题分析周期较长，仿真分析很难满足工程上对应力强度因子预测的效率要求。因此，探索如何有效地把理论分析与数值方法结合起来，通过少量仿真分析修正理论解，从而实现针对不同工况下的裂纹参数，快速、准确地计算应力强度因子。

本文开展压缩机轴流叶片单边贯穿型裂纹应力强度因子的计算方法研究。将离心力作用下的薄壁金属叶片简化为矩形平板，在应力强度因子权函数解法的基础上，通过扩展有限元数值分析修正理论解中的形状因子，并研究裂纹位置和板厚对理论解的影响，通过数据拟合提出计算叶片裂纹应力强度因子的半解析公式。进一步建立实际压缩机叶片模型，基于扩展有限元仿真分析修正基于薄板的半解析公式，建立实际叶片I 型应力强度因子的半解析解公式。本文将应力强度因子的经典理论解适用范围从中心型裂纹推广到非中心型裂纹，综合考虑叶片厚度和形状的影响，为快速评估实际轴流压缩机旋转叶片的应力强度因子提供新方法，也为后续叶片剩余强度及寿命分析提供分析基础。

1 矩形板I型裂纹应力强度因子

轴流压缩机叶片示意图如图1所示，由于叶片厚度远远小于长度和宽度，故可将其简化为xy面内、以角速度w旋转的矩形平板，简化模型如图2所示。设板的长度为h，宽度为b，厚度为W，板的底边坐标y=-h/2。考虑在平面应力状态下产生的裂纹以I型（张开型）为主，在左侧板边处(x=-b)设置了贯穿板厚度的长度为a的裂纹。与其他研究不同的是：本文设置的裂纹不限于板的中心位置（y=0），即起裂位置可以是任意值y=l-h/2-r，且旋转轴不是板底边，而是平行于x轴的任意AB(y=-h/2-r)轴。

图1 轴流压缩机叶片示意图Fig.1 Schematic diagram of axial flow compressor blade

图2 含单边贯穿型裂纹矩形板Fig.2 Rectangular plate of single-edged penetrating crack

依据断裂力学理论，裂纹尖端的I型应力强度因子的一般形式为[14]：

其中，σy为裂纹尖端名义应力，FI为形状因子。

当叶片旋转时，裂纹所在横截面处的离心力为：

式中，ρ为材料质量密度。于是，裂纹面上任意点x 处的拉应力为：

由此得到旋转叶片上裂纹应力强度因子：

对于含中心裂纹、平面应力状态的线弹性矩形薄板，大多采用权函数法推导I型裂纹尖端的应力强度因子理论解。基于此方法，在板中面(z=0)上由名义应力σy(x)产生的I型裂纹应力强度因子表达式为[14]：

其中，m(x，a)为权函数，其与载荷条件无关，由裂纹体结构自身的性质决定，且具有唯一性。E为材料弹性模量；vy(x，a)为外载荷引起的裂纹面张开位移。根据文献[3]，权函数取为x∈(0,a)内一致收敛的级数

式中，Mn为待定系数。根据文献[3]，级数展开式一般取到n=3即可达到精度要求。

在理论上，若已知互相独立的3个应力场对应的应力强度因子，则可根据式(5)、式(7)确定M1；M2；M3。对于中心型裂纹，由于位移场的对称性，裂纹面的位移场曲率为零，即

可知M2为常数，不影响应力强度因子解，M1；M3表达式为：

其中，F1和F2是矩形板在远场均布拉力和纯弯曲两个独立载荷作用下的形状因子，通过解析推导或数值计算得到。

对于受离心力作用的细长轴流叶片，可使用Brown等[15]提出的适于有限宽、无限长板条的公式计算形状因子F1。即

其中，a/b＜0.6 的中心型裂纹，上式计算误差不超过0.5%。也可采用Tada等[16]提出的理论解计算形状因子F1，即

对于a/b＜0.2 的中心型短裂纹，上式的计算误差不超过1%；对更长的裂纹，其误差不超过0.5%。

若纯弯加载情形，可使用文献[15]提出的公式计算形状因子F2，表达式为

对于a/b＜0.6 的中心型裂纹，上式计算误差不超过0.5%。将式(10)、(12)依次代入式(9)、(7)，再积分式，即可求出应力强度因子KI的理论解。

需要指出的是：上述理论解的前提条件是中心型裂纹假设，而实际裂纹可能萌生在叶片的任意位置，此时的权函数应与式(7)有所不同。此外，实际叶片具有一定厚度，裂纹尖端是一条曲线而不是平面情形下的单个点。因此，上述应力强度因子理论解在计算实际轴流叶片时存在一定误差，需对其进行验证并做必要的修正。

2 应力强度因子理论解修正

本节基于ABAQUS软件的扩展有限元法(XFEM)模块，选取不同的裂纹参数分析单边贯穿型裂纹的I型应力强度因子，验证形状因子并对其修正，讨论裂纹位置对理论解精度的影响。

扩展有限元法是当前数值上解决不连续力学问题的最佳选择之一，已广泛应用于含裂纹结构的应力强度因子、J积分等断裂力学性能分析。它在常规位移场插值的基础上引入了Heaviside函数描述裂纹两侧的间断位移，以及裂尖函数反映裂纹尖端的应力奇异[11]。扩展有限元方法优于传统有限元方法的一个重要方面是裂纹面无需创建在单元的边界，即裂纹独立于网格，并且无须随着裂纹的扩展而重新划分网格。在用扩展有限元计算应力强度因子时，积分路径原则上由裂纹下表面开始，沿封闭围道环绕裂纹尖端，终止于裂纹上表面。本文在实际计算时尝试了不同的围道、富集半径、裂尖网格密度，最终选出精度满意的参数组合，保证了应力强度因子解的收敛性。

2.1 形状因子验证

采用ABAQUS 软件，建立带裂纹的矩形平板仿真模型如图3(a)所示，模型包括50898 个节点，16765 个单元。如图3(b)所示，在裂纹尖端附近布置长度为1mm的奇应变三角形单元，采用环形区域积分，设置积分路径为4条，选择最大主应力准则为裂纹初始化准则。

图3 带裂纹的矩形板仿真模型Fig.3 Simulation model of a rectangular plate with crack

令板宽b=100mm，弹性模量E=200GPa，板顶边和底边受100MPa 均布拉力作用。给定50 组不同的板长度h和裂纹长度a，根据ABAQUS 仿真分析得到应力强度因子KI，再通过式反推出裂纹的形状因子。依据数值分析得到的形状因子与理论解对比，结果如图4所示。

图4 形状因子随裂纹长度的变化Fig.4 Variation of crack shape factor with crack length

从图4可以看出，形状因子数值解与式预测趋势一致，均随裂纹延伸而单调增大。其中，长宽比h/b越小，形状因子仿真值与理论值的偏差越大，对于长度大于0.05b 的裂纹，当h/b＜0.8 时，仿真值与理论值的偏差已不可忽略。这种情形多见于轴流压缩机的高压级短叶片，在以往文献中对此问题研究不够充分。

基于上述分析结果，对形状因子的理论解进行修正。本文采用LM法和通用全局优化算法（UGO），拟合形状因子的多项式解如下：其中，α=a/b，β=h/b，c0-c9为系数。对于本例，它们依次为：1.9472、9.1361、-4.4591、-1.9440、-14.3208、6.3277、9.7379、1.3786、6.5427、-2.6999。拟合结果与仿真结果的相关系数为0.998，说明式(13)具有较高精度。形状因子的函数曲面如图5所示。

图5 形状因子关于裂纹长度及板长宽比的曲面图Fig.5 Surface diagram of the crack shape factor with crack length and plate’s length to width ratio

2.2 裂纹位置对应力强度因子解的影响

设板长度h=180mm，宽度为b=100mm，裂纹位置l分别取10mm、36mm、60mm 和90mm，以5mm 为间隔计算不同裂纹长度对应的应力强度因子理论解和有限元解，结果如图6所示。

图6 不同裂纹位置的应力强度因子数值解和理论解Fig.6 Numerical and theoretical solutions of stress intensity factor at different crack locations

由图6(b)、(c)和(d)中可见，裂纹位置分别为l=36mm、60mm 和90mm 时，应力强度因子理论解与扩展有限元解吻合得较好。图6(a)中l=10mm 时，即裂纹足够接近板底边时，理论解与有限元解偏差较大，原因在于理论解建立在假设裂纹处于平板中心线上，计算精度随裂纹位置距离中心线越远而逐渐下降。

根据以上分析，采用拟合算法对式的应力强度因子理论解进行修正。以裂纹长度和位置为参数，用LM法和通用全局优化算法（UGO）拟合出应力强度因子的半解析解如下：

对本算例，系数d0～d9依次为189.2635、18.8499、6.3932、-0.1711、0.4155、-0.2272、0.01228、-0.00349、-0.00403、0.00166。上述拟合的相关系数为0.999，说明精度满足要求。

2.3 板厚对应力强度因子解的影响

分别取板厚W=6mm、10mm、20mm 和40mm，并分别设置长度a=10mm、30mm 的裂纹。由于裂纹尖端应力场关于板中面对称分布，故图7只展示从上表面到中面的I型应力强度因子XFEM解。其中横坐标为沿厚度方向的无量纲尺寸2z/W，图中水平线代表应力强度因子的理论解。

图7 应力强度因子沿板厚度的变化Fig.7 Variation of the stress intensity factor along plate thickness

由图7可见，应力强度因子随板厚的变化量大致不超过1%。在厚度方向上，中面(z=0)的应力强度因子值最大，上、下表面的因子值最小，与实际裂纹尖端接近于圆弧的规律相符。

考虑板厚度影响后，可对式的应力强度因子进行修正，得到受离心力作用的矩形平板I型应力强度因子半解析解：

依据仿真分析结果，修正系数dw在1.06至1.09之间波动，为保守起见，对实际轴流叶片建议取dw=1.1。

3 轴流叶片应力强度因子的半解析解

实际压缩机轴流叶片虽是薄壁结构，但其安装角沿叶片高度方向不断变化，从叶根至叶顶呈现出一定的空间扭曲，且横截面的厚度中间大、两边小。因此，基于平板公式计算的应力强度因子必然存在误差。以下针对含单边贯穿裂纹的实际压缩机轴流叶片，基于扩展有限元法对半解析公式进行进一步修正以提高其预测精度。

仿真分析对象为图1 所示的某轴流压缩机第6 级叶片，宽度b=100mm，高宽比为1.8。有限元模型网格尺寸为5mm，在叶片侧边及叶根处，网格尺寸加密为1mm，模型共划分127867 个单元。计算时分别在叶片上部和底边两个不同位置（l=36mm和l=90mm）设置不同长度的贯穿裂纹。图8给出了I型应力强度因子随裂纹扩展的变化趋势，同时给出了按式(15)算出的半解析解结果。

图8 实际叶片应力强度因子随裂纹长度的变化Fig.8 Variation of actual blade stress intensity factor with crack length

从图中看到，半解析解与扩展有限元解变化趋势基本一致，但数值上存在一定偏差，原因是实际叶片横截面积沿叶高方向逐渐减小，裂纹面实际承受的拉应力小于式计算值，对应的应力强度因子值必然小于等截面薄板情况。此外，由于叶片形状扭曲，受力后实际处于拉伸、剪切和扭转共存的组合应力状态，裂纹前缘变形复杂，同时存在I型、II型和III型应力强度因子，与推导理论解时假设的平面应力状态差别较大，故必然产生一定偏差。

考虑到完好叶片的宽、高、厚度均已确定，因此影响断裂力学性能的主要参数仅为裂纹长度及起裂位置。引入修正系数ζ(α，β)，其中β=l/h为无量纲裂纹位置，将式（16）的半解析解改写为

根据图7 的扩展有限元仿真结果，采用LM 法和UGO 法拟合出修正系数ζ(α，β)关于裂纹长度和裂纹位置的响应面。本例设置了75个(α，β)组合，得到修正系数ζ(α，β)的响应面如图9 所示，其中黑点为原始数据点。对应的函数表达式为:

图9 修正系数随裂纹长度和裂纹位置的变化Fig.9 The variation of modified parameter with crack lengths and crack locations

在本文算例中，表达式(17)的系数g0～g9依次是：0.19988,-0.05989,1.52846,2.26096,-0.22049,-4.93579,-2.30420,-0.34714,0.34261,5.30313。由此可得实际压缩机轴流叶片应力强度因子的半解析解，式（16）可用于计算不同位置、不同长度的贯穿型裂纹尖端应力强度因子。

分析图9 中曲面的变化特征可知，当裂纹长度a不变时，若改变裂纹位置l，则修正系数ζ(α，β)随裂纹位置l的增大而增大；当裂纹位置l不变时，改变裂纹长度a，则修正系数ζ(α，β)随裂纹长度a的增加而迅速增大，表明了基于数值仿真来修正半解析解的必要性。

尽管本例叶片的裂纹长度和裂纹位置的数据离散度均较大，但应力强度因子半解析解与扩展有限元法解的相对误差仍在5%至10%之间，表明精度满足工程要求。

4 结论

本文针对轴流压缩机行业快速评估含裂纹叶片应力强度因子的实际需求，根据结构和载荷特点建立了含单边贯穿性裂纹的矩形平板模型，结合断裂力学理论解公式，采用扩展有限元方法研究了形状因子、板厚及裂纹位置对I 型裂纹应力强度因子理论解精度的影响，通过拟合提出了改进后的半解析解公式。进一步开展实际轴流叶片的仿真分析，修正基于平板模型的半解析解公式，构建了用于计算包含不同位置、不同长度的贯穿型裂纹的实际轴流叶片应力强度因子的半解析解。主要结论如下：

1）依据断裂力学平板理论预测轴流叶片裂纹尖端的I 型应力强度因子存在一定误差，对于高压级短叶片，形状因子影响不可忽略；裂纹位置距离平板中心线越远，计算精度越低；应力强度因子随板厚的变化量大致不超过1%。

2）实际轴流叶片裂纹尖端的I 型应力强度因子与基于平板模型的计算结果存在差异，裂纹越长、裂纹距离板中心线越远，计算误差越大。

3）本文提出的实际轴流叶片I 型应力强度因子的半解析解，考虑了实际叶片厚度及叶片形状变化、不同裂纹位置和裂纹长度的影响，能够满足工程对计算精度的要求，与理论解和仿真分析相比，可以满足一般的工程使用需求,有效提升轴流叶片应力强度因子的评估精度及效率，节省断裂力学分析所需的工作量。

4）采用代理模型技术及优化设计算法有望进一步改善解的预测精度，此外，未来还可探索代理模型、深度学习等技术在本研究方向的工程应用。

致谢：本文得到国家自然科学基金及国家科技重大专项项目资助，关萍同志参加了有关分析计算以及审稿人对论文写作提出了宝贵意见，在此一并致谢。