基于有限元方法的航空发动机叶片应力强度因子计算_李春旺

基金项目：国家自然科学基金(10804130) 收稿日期：2012-05-31 修回日期：2013-03-08 第一作者简介：李春旺，男，1977 年生，博士，西北工业大学航空学院，博士后（空军工程大学理学院，讲师）；研究方向——航空结构的振动、疲 劳与断裂。 E-mail：Li_chw@163.com 应 用 力 学 学 报 CHINESE JOURNAL OF APPLIED MECHANICS 第 30 卷 第 3 期 2013 年 6 月 Vol.30 No.3 Jun. 2013 文章编号：1000- 4939(2013) 03-0373-05 基于有限元方法的航空发动机叶片 应力强度因子计算 李春旺 1,2 罗秀芹 2 杨百愚 2 张贵斌 3 李亚智 1 张忠平 2 (1 西北工业大学航空学院 710072 西安；2 空军工程大学理学院 710051 西安； 3 西安航空动力股份公司军事代室 710021 西安) 摘要：为研究叶片裂纹尖端的应力奇异性，以某型航空发动机压气机叶片为例，利用有限元方法 研究了叶片裂纹尖端应力强度因子的计算方法，并研究了旋转叶片振动状态下裂尖应力强度因子 随裂纹长度的变化规律。建立计算模型时，在裂纹尖端划分了三维奇异单元，在裂尖外围划分了 过渡单元。计算结果表明：研究旋转叶片振动状态下的裂尖应力奇异性，仅利用 I 型应力强度因 子就具有足够的精度；对于同一裂纹，绝大多数情况下叶盆面应力强度因子大于叶背面应力强度 因子，故研究叶片应力强度因子时只需研究叶盆应力强度因子即可；随着裂纹扩展，叶盆面 I 型 应力强度因子不断增大。本文的研究方法及结论为进一步研究叶片的裂纹扩展规律及损伤容限奠 定了基础。 关键词：含裂纹叶片；应力强度因子；有限元方法；三维奇异单元；过渡单元 中图分类号：V235.1 文献标识码：A DOI：10.11776/cjam.30.03.B069 1 引 言 关于无裂纹叶片的振动特性，国内外专家学者 已进行了广泛深入的研究[1-6]，但对含裂纹叶片振动 状态下的应力-应变分布的研究较少，而且仅有的研 究在计算裂尖应力强度因子时，都没有充分考虑裂 纹尖端的奇异性[7-9]。 研究叶片裂尖应力强度因子的关键在于裂纹尖 端奇异单元的划分，对于振动导致的叶片裂纹，由于 叶片形状十分复杂，无论是采用网格加密还是采用二 维平面单元都无法反映裂尖的应力奇异性，此时，需 要采用三维奇异单元。基于这一思路，本文以某型航 空发动机压气机叶片为例，以叶片试验为背景，研究 叶片一阶弯曲振动过程中裂纹尖端应力强度因子的 计算方法以及应力强度因子随裂纹长度的变化规律。 2 平板的应力强度因子计算 2.1 平板的计算模型 为在裂纹尖端划分三维奇异单元并验证其正 确性，首先计算平板的应力强度因子。图 1 显示了 平板的形状及裂纹位置，平板长 300mm、宽 60mm、 厚 6mm。平板材料与叶片材料相同，都是锻造铝合 金 LD7-1，室温下弹性模量、质量密度、泊松比分 374 应 用 力 学 学 报 第 30 卷 别为 E= 107.10 10× N·m-2、 ρ = 32.73 10× kg·m-3； μ =0.33。 采用 20 节点等参单元 SOLID95 对图 1 进行网 格划分，在裂纹前缘采用四分之一节点退化奇异单 元，以适应裂尖的应力奇异性。同时，为了实现从 奇异单元到常规单元的平稳过渡、增加裂纹前缘的 网格密度，在奇异单元外围划分过渡单元。由于 ANSYS 不能自动生成三维奇异单元，所以采用手 动方式在裂尖附近划分网格，步骤如下。 1) 以裂纹正前方为 X ′轴、裂纹面法线方向为Y ′ 轴、裂纹前缘线为 Z ′轴建立局部坐标系，并在此坐 标系下以裂纹前缘线为中轴，在裂尖附近分割出用 于划分奇异单元的实心圆柱体和用于划分过渡单元 的空心圆柱体(见图 2)，并将空心圆柱体与整体叶片 粘贴在一起。 2) 分割实心圆柱体及其边界线，得到创建奇异单 元所需的节点后删除实心圆柱体。每个奇异单元的 节点分布如图 3 所示[10-11]，按照图 3 所示的顺序拾 取节点，即可直接用 20 个节点生成如图 4 所示的三 维奇异单元。模型中沿裂纹尖端线排列有 4 组奇异 单元，每组的 16 个奇异单元围绕裂纹前缘线分布， 它们有一个共用边，即裂纹前缘线，这样裂纹尖端 就完全被三维 SOLID95 奇异单元包围(见图 5)，并 且这些奇异单元正好紧贴外边的空心圆柱体。 3) 使用扫略法在空心圆柱体中划分过渡单元，并 进行控制，使过渡单元呈中心对称分布并将奇 异单元层层包围，每层过渡单元的数目及其分布与 奇异单元完全一致(见图 6)。这样，奇异单元和过渡 单元相同位置的节点就具有相同的 Z ′坐标，便于在 计算应力强度因子时定义积分路径。 4) 虽然奇异单元外层的节点和过渡单元内层的对 应节点在同一空间位置，但却彼此独立，为保证它 们具有相同的自由度，焊接两种单元对应节点，得 到叶片裂纹附近的网格分布(见图 7)。 图 1 平板形状及裂纹 图 2 局部坐标系及裂尖附近区域划分 Fig.1 The flat blade and crack Fig.2 Local coordinate and element parts 图 3 奇异单元节点及顺序 图 4 奇异单元 Fig.3 Node order of singular element Fig.4 Singular element 图 5 裂纹尖端奇异单元分布 图 6 过渡单元 Fig.5 Singular elements around crack tip Fig.6 Transition elements 图 7 平板裂纹附近网格分布 图 8 拉伸载荷下平板裂尖附近的应力 Fig.7 Elements distribution Fig.8 The stress distribution of the around the flat blade’s crack tip crack region under tensile loads 2.2 拉伸载荷下平板的应力强度因子计算 为了验证计算模型的正确性，这里需要计算拉 伸载荷作用下平板的应力强度因子。计算时，取裂 纹长度为 10mm，固支平板的下端面，并在上端面 施加 9730N 的拉伸载荷，计算结果见表 1。其中：d 为考察点与平板表面的距离，量纲为 mm；KΙ 、KⅡ、 KⅢ分别为 I、II、III 型应力强度因子，量纲为 N·m-3/2。 表 1 拉伸载荷下平板裂纹尖端的应力强度因子 Tab.1 Stress intensity factors of the flat blade’s crack tip d/mm 0 1 2 3 4 5 6 KΙ /×10 6 N·m-3/2 5.35 5.11 5.14 5.09 5.14 5.11 5.35 KⅡ/ N·m-3/2 70.3 29.6 12.4 2.2 7.6 18.0 87.9 KⅢ/N·m-3/2 209.7 266.6 202.7 260.1 203.1 268.2 190.5 图 8 显示了裂尖的应力分布，图中的“蝴蝶翅 膀”定性说明了计算模型的正确性，为了定量证明计 算模型的正确性，考虑如下理论公式[12] 1 .12 πK aσΙ = (1) 第 3 期 李春旺，等：基于有限元方法的航空发动机叶片应力强度因子计算 375 式中：a 是裂纹长度；σ 是作用于板端面的名义应 力。利用平板的横截面尺寸及作用载荷得到名义应 力σ 之后，再根据式(1)可得 KΙ 为 5.36 610× N·m-3/2， 这一数值与表 1 中的 KΙ 近似相等。另外，表 1 中的 数据显示：一方面，沿板的厚度方向 KΙ 几乎不变， 另一方面，与 KΙ 相比 KⅡ和 KⅢ可以忽略。这都充分 说明了计算模型的正确性，因此，可以利用所建模 型计算裂纹尖端的应力强度因子。 3 叶片应力强度因子的计算 3.1 计算模型 试验发现，当某型航空发动机转速达 5378r/min 时，发动机压气机叶片在外加激振力作用下发生一 阶弯曲共振，并在叶片排气边产生裂纹，最大应力 位于距离缘板 70mm 处，应力值为 137.8MPa。因此， 本文考虑最恶劣情况，假设叶片裂纹起源于排气边， 裂纹恰与叶片最大应力部位同高。显然，这样的裂 纹位置是最危险的，在该位置裂纹尖端应力强度因 子范围最大，裂纹扩展速率也最快。图 9 显示了叶 片的实体模型及裂纹位置，图 10 显示了裂尖附近的 网格划分。 图 9 叶片及裂纹位置 图 10 裂尖附近的网格 Fig.9 The blade and crack Fig.10 Elements around the crack tip 3.2 计算结果 计算过程中，依据发动机轮盘大小，令叶片绕发 动机轴以半径 160mm、速度 5378r/min 旋转，同时， 固定叶片榫头，调整叶尖排气边位移，保证叶身最 大应力部位距离缘板 70mm 且交变应力达到 137.8MPa，计算不同裂纹长度对应的应力强度因 子。裂纹长度为 35mm 时，叶片振动达到最大位移 时裂尖应力的分布见图 11，各个裂纹长度下 z 轴方 向的应力 zσ 及应力强度因子的计算结果见表 2，应 力强度因子随裂纹长度的变化见图 12、图 13。 表 2 裂尖应力及应力强度因子 Tab.2 Stresses of crack tip and stress intensity factors a/min 振动方向 (vibrating direction) 叶背面(the convex face) 叶盆面(the concave face) KI/ ×106N·m-3/2 KⅡ/ ×106N·m-3/2 KⅢ/ ×106N·m-3/2 zσ /MPa KⅠ/ ×106N·m-3/2 KⅡ/ ×106N·m-3/2 KⅢ/ ×106N·m-3/2 zσ /MPa Ke/ ×106N·m-3/2 δ /(%) 8 +y 6.409 0.728 0.450 124.20 10.382 1.153 0 242.47 10.382 0 -y 0 0.351 0.344 -40.11 0 0.547 0 -143.91 0 0 10 +y 6.143 0.693 0.506 118.93 11.557 1.326 0.140 271.66 11.558 0.009 -y 0 0.278 0.359 -16.76 0 0.531 0 -149.64 0 0 15 +y 4.854 0.585 0.611 92.27 14.012 1.707 0.373 326.41 14.022 0.07 -y 1.963 0.103 0.316 52.67 0 0.418 0 -146.11 0 0 20 +y 3.198 0.493 0.790 54.72 16.637 2.208 0.676 600.35 16.664 0.16 -y 5.858 0 0.197 140.69 0 0.245 0 -199.36 0 0 25 +y 1.512 0.431 1.140 13.98 19.878 2.869 1.003 470.17 19.928 0.25 -y 10.358 0.140 0 241.98 0 0 0.105 -96.56 0 0 30 +y 0 0.399 1.773 -21.45 24.285 3.762 1.046 572.88 24.330 0.19 -y 15.662 0.170 0.266 365.34 0 0.408 0.226 -38.39 0.318 — 35 +y 0 0.358 2.809 -39.70 30.649 4.961 2.004 732.19 30.778 0.42 -y 22.203 0.126 0.680 522.96 3.656 0.990 0.448 62.66 3.710 1.48 40 +y 0 0.253 4.445 -25.36 40.878 6.616 2.078 960.07 40.983 0.26 -y 31.287 0 1.335 718.43 11.555 1.883 0.832 237.00 11.614 0.51 45 +y 5.462 0 7.051 66.47 58.999 8.976 3.086 1403.9 59.159 0.27 -y 45.555 0.382 2.440 1043.0 26.837 3.249 1.050 594.52 26.878 0.15 50 +y 19.745 0.885 11.535 386.74 98.194 12.832 5.080 2345.3 98.454 0.26 -y 72.728 1.229 4.589 1715.0 62.061 5.616 2.088 1426.2 62.131 0.11 376 应 用 力 学 学 报 第 30 卷 计算应力 zσ 的目的在于确定 I 型应力强度因子。 当 0zσ > 时，裂纹张开， 0KΙ > ；当 0zσ < 时，裂纹 闭合， 0KΙ = 。通过对表 2 的数据及图 11~图 13 的对 比，可得到如下结论。 1) 叶片工作状态下的一阶弯曲振动，其叶盆与叶 背应力是离心应力与振动应力的合成。由表 2 可见， 叶盆、叶背应力相差很大，造成差异的原因是叶盆 与叶背表面形状的差异。 2) 对于同一裂纹，就 I 型应力强度因子而言，叶 背面与叶盆面的应力强度因子差距较大，而且绝大 多数情况下叶背面应力强度因子小于叶盆面应力强 度因子，这与表 1 中 KΙ 关于中面层对称的事实形成 鲜明对比。之所以这样，原因归结为平板形状的规 则性与叶片表面形状的复杂性，应力强度因子的差 异进一步表明了裂尖划分三维奇异单元的必要性。 另外，研究叶片应力强度因子时只需研究叶盆应力 强度因子即可。 3) 沿-y 方向振动时叶背应力强度因子与沿+y 方向 振动时叶盆应力强度因子的变化规律相近，都随着裂 纹扩展而不断增大。与此类似，沿+y 方向振动时叶背 应力强度因子与沿-y 方向振动时叶盆应力强度因子 的变化规律也相近，随着裂纹扩展，应力强度因子先 减小(不变)而后增大。这是由于沿-y、+y 方向振动时 叶背、叶盆的离心应力和弯曲振动应力都是拉应力； 而沿+y、-y 方向振动时叶背、叶盆的离心应力是拉应 力，而弯曲振动应力却是压应力。 图 11 裂纹部位应力分布 Fig.11 The stress distribution of the crack region 图 12 叶背应力强度因子随裂纹长度的变化 Fig.12 Curve of stress intensity factors on the suction side 图 13 叶盆应力强度因子随裂纹长度的变化 Fig.13 Curve of stress intensity factors on the pressure side 3.3 Ⅲ型应力强度因子的影响 叶片的振动疲劳试验结果表明，裂纹扩展方向与 叶身高度方向基本垂直。这说明决定裂纹扩展的是 I、 Ⅲ型应力强度因子，因为 I、Ⅲ型应力强度因子的作 用效果是使裂纹沿原(预置)方向扩展，而 II 型应力强 度因子的作用效果是使裂纹改变扩展方向[12]。 既然存在复合型应力强度因子，暂且引入等效 应力强度因子的概念。根据应变能密度因子理论， 与 I、Ⅲ型应力强度因子相关的等效应力强度因子 表达式为[12] 2 2 e 1 1 K K KμμΙ ΙΙΙ += + − (2) 式中 μ 是泊松比。对于所研究的叶片材料，μ 取值 为 0.33。 将表 2 中叶盆面的 KΙ 、KⅢ代入式(2)可得等效 应力强度因子 eK ，结果也列于表 2 中。由表可见， KΙ eK≈ ，其相对偏差δ 最大仅为 1.48%。这说明Ⅲ 型应力强度因子对裂纹扩展的影响可以不计，即计 算叶片工作状态下的应力强度因子时不需要引入等 效应力强度因子的概念，仅利用 I型应力强度因子 就已经具有足够的精度。 4 结 论 本文以某型航空发动机压气机叶片为例，通过 在裂纹尖端划分三维奇异单元、在裂尖外围划分过 渡单元，得到了叶片裂尖应力强度因子的计算方法， 研究了裂尖应力奇异性及裂尖应力强度因子随裂纹 长度的变化规律，得到如下结论。 1) 研究叶片工作状态下的裂尖应力奇异性，仅利 用 I 型应力强度因子就已经具有足够的精度。对于 所研究的叶片，在向叶盆面弯曲时，KI与 Ke的相对 最大偏差仅为 1.48%。 第 3 期 李春旺，等：基于有限元方法的航空发动机叶片应力强度因子计算 377 2) 对于同一裂纹，绝大多数情况下，叶背面应力强 度因子小于叶盆面应力强度因子，因此研究叶片应力 强度因子时只需研究叶盆面应力强度因子即可。 3) 随着裂纹扩展，叶盆面 I 型应力强度因子不断 增大。 参 考 文 献 [1] 孙强，汪波，柴桥，等．固支状态对压气机叶片自振频率计算结 果的影响[J].空军工程大学学报，2003，4(4)：8-10．(Sun Qiang， Wang Bo，Chai Qiao，et al. The effects of the compressor blade fixed states to the nature frequency calculated results[J].Journal of Air Force Engineering University，2003，4(4)：8-10．(in Chinese)) [2] 孙强，张忠平，柴桥，等．航空发动机压气机叶片振动频率与温 度的关系[J].应用力学学报，2004，21(4)：137-139．(Sun Qiang， Zhang Zhongping，Cai Qiao，et al．The relationship between the aero-engine compressor blade nature frequency and the tempreture[J]． Journal of Applied Mechanics， 2004， 21(4)： 137-139．(in Chinese)) [3] 杨文庆，孙强，马龙，等．某型航空发动机压气机叶片振动静频 与动频的关系[J].空军工程大学学报，2005，6(5)：5-7．(Yang Wenqing，Sun Qiang，Ma Long，et al. The relationship between the dynamic frequency and the static ones for an aero-engine compressor blade[J]．Journal of Air Force Engineering University，2005，6(5)： 5-7．(in Chinese)) [4] 孙强，张忠平，李春旺，等．航空发动机销钉式叶片固有振动频 率计算[J].应用力学学报，2006，23(3)：470-472．(Sun Qiang， Zhang Zhongping，Li Chunwang，et al.The study of the nature frequency for aero-engine pin fixed blade[J]．Chinese Journal of Applied Mechanics，2006，23(3)：470-472．(in Chinese)) [5] 张忠平，孙强，李春旺，等．航空发动机压气机叶片振动疲劳寿 命与 af 值的关系[J]．应用力学学报，2006，23(3)：459-461．(Zhang Zhongping，Sun Qiang，Li Chunwang，et al．The relationship between the fatigue life and the af value for aero-engine compressor blades[J]．Chinese Journal of Applied Mechanics，2006，23(3)： 459-461． [6] 谢永慧，王乐天，安宁，等．汽轮机叶片振动特性的三维有限元 分析[J]．机械强度，1997，19(4)： 1-5．(Xie Yonghui，Wang Letian， An Ning，et al. The three dimensional finite element anylasis for steam turbine blade vibration[J]. Journal of Mechnical Strength， 1997，19(4)：1-5．(in Chinese)) [7] 葛永庆，安连锁.裂纹参数对叶片固有频率影响的研究[J].动力工 程，2008，28(4)：519-522．(Ge Yongqing，An Liansuo.The effects of the crack parameters to the blade nature frequency[J].Journal of Dynamic Engineering，2008，28(4)：519-522．(in Chinese)) [8] 胡殿印，王荣桥，邓俊．基于有限元方法的裂纹扩展寿命预测[J]．机 械强度，2009，31(2)：264-268．(Hu Dianyin，Wang Rongqiao， Deng Jun. The prediction of the crack propagation life from the finite element method[J].Journal of Mechnical Strength，2009，31(2)： 264-268．(in Chinese)) [9] 祁顺彬．基于 ANSYS 的Ⅰ-Ⅱ复合型裂缝应力强度因子计算[J]. 水利与建筑工程学报，2008，6(2)：86-90．(Qi Shunbin.The stress intensity factor for Ⅰ-Ⅱ mixed crack from ANSYS[J]．Journal of Water Conservancy and Architectural Engineering，2008，6(2)： 86-90．(in Chinese)) [10] ANSYS Inc.Release 10.0 documentation for ANSYS[R]．Pittsburgh， US：ANSYS Inc，2005． [11] 樊鸿，张盛，王启智.复合型三维裂纹应力强度因子计算方法的研 究[J].四川大学学报：工程科学版，2009，41(4)：48-52．(Fan Hong， Zhang Sheng，Wang Qizhi.Study to the method of the stress intensity factors for three dimensional crack[J].Journal of Sichuan University， 2009，41(4)：48-52．(in Chinese)) [12] 洪起超.工程断裂力学基础[M].上海：上海交通大学出版社， 2000．(Hong Qichao.Engnerring fracture mechanics[M].Shanghai： The Publishing House of Shanghai University，2000．(in Chinese)) CHINESE JOURNAL OF APPLIED MECHANICS No.3 vii Keywords: transverse shear deformation, rectangular moderately thick shallow shells, small-deflection buckling, fundamental differential equations, double trigonometric series, characteristic equation, critical load. Dynamic modification and application of wing double-beam model He Erming1 Chen Yi1 Li Yulong1 Hu Xizhi2 Han Feng2 (1 School of Aeronautics, Northwestern Polytechnical University, 710072, Xian, China; 2 Shanghai Research Institute of Aircraft Design, 200232, Shanghai, China) Abstract: In order to investigate the characteristics of vibration load at each frame along airframe transmitted from the engine vibration though the wing structure, a new double-beam dynamical model with 28 ribs of airliner’s wing is presented in this paper. According to the results of the sensitivity analysis for unidentified parameters of the original double-beam model, the effective structural design variables, such as the bending stiffness, torsion stiffness of wing beams and ribs are chosen. Through modifying structural parameters of key wing ribs, the double-beam model’s dynamical characteristic has been improved, in which the maximum error in first six modal frequencies between the double-beam model and a all-size wing FEM model in Ref.[8] is lower than 11.11%. Then a half-full dynamical model with the Pylon-Wing-Airframe of airliner is built, and the vibration load at each frame along airframe transmitted from engine vibration through Mounting-Pylon-Wing is calculated based on it. The acceleration distribution along airframe at main engine feature frequencies under cruise condition is obtained. These results prove that the double-beam dynamical model is better to meet the research needs for the airframe vibration calculation and transmission path analysis of engine vibration, and provides the technical preparation for vibration reduction & isolated mounting system design of wing and pylon structures. Keywords: wing double-beam model, sensitivity analysis, dynamic modification, vibration transmission. Stress intensity factor calculation for aero-engine compressor blade with finite element method Li Chunwang1,2 Luo Xiuqin2 Yang Baiyu2 Zhang Guibin3 Li Yazhi1 Zhang Zhongping2 (1 Northwestern Polytechnic University, 710072, Xi’an, China; 2 Air Force Engineering University, 710051, Xi’an, China; 3 AVIC Xi’an Aero-Engine Limited, 710021, Xi’an, China) Abstract: In order to study the stress singularity at the blade crack tip, an aero-engine compressor blade with crack is taken as example, stress intensity factor from the finite element method is calculated. Then, the change of the stress intensity factor with the crack length is studied. To construct the calculation model, the three dimensional singular elements at crack tip and the transition elements around crack tip are meshed. The conclusions are that only mode I stress intensity factor of concave is adequate to describe the stress singularity for the blade, and that the mode I stress intensity factor of concave increases with the crack propagation. Both the stress intensity factor calculation method and the conclusions of the present paper provide the base for further investigation of the crack growth rate and the damage tolerance of the blade. Keywords: blade with crack, stress intensity factor, finite element method, three dimensional singular element, transition element.