基于ABAQUS的某无人机机翼肋板裂纹扩展仿真

王海鹏，吴永东

(南昌航空大学，南昌 330063)

0 引言

无人机(Unmanned Aerial Vehicle，UAV)是一种由动力驱动、机上无人驾驶、可重复使用的航空器的简称。具有隐蔽性好，生命力强;造价低廉，不惧伤亡;起降简单，操作灵活等优点［1］。机翼各部件中存在的微小裂纹，在飞机的使用过程中会不断扩展，最终导致机翼的断裂破坏，严重威胁飞行安全。

裂纹扩展属于典型的强不连续问题，一直是工程力学分析的热点、重点和难点问题［2］。在模拟裂纹扩展时，Shephard M等［3］在每次裂纹扩后把裂纹作为边界条件的一部分，然后重新划分网格。Klein P等［4］引用可虚拟联接单元的概念，用虚拟联接模型来解决裂纹扩展问题。王慰军［5］编写了基于ABAQUS平台的裂纹扩展模拟软件。Ortiz M等［6］发展了三维有限变形的内聚元和一系列不可逆内聚力关系，该方法可以近似地模拟出裂尖的轨迹。娄路亮等［7］利用无网格法模拟了裂纹的扩展，得到的裂纹扩展路径与实验值吻合的很好。Nicolas Moёs等［8］利用扩展有限元法对裂纹扩展进行模拟，该方法允许裂纹与网格相对独立，因此在模拟过程中无需重新划分网格。目前，国内外对机翼的裂纹扩展研究较少，对机翼裂纹的扩展过程与形式认识不足。

研究采用Python语言对ABAQUS进行二次开发，实现模拟裂纹自动扩展过程中有限元网格的自动重划分功能，采用最大周向应力准则，计算裂纹扩展方向，应力强度因子等断裂参量。将机翼肋板简化为二维平板结构，在肋板与梁的连接处放置裂纹，模拟不同初始角度的裂纹的扩展过程，得到相关的力学参量，总结裂纹扩展路径的总体趋势。

1 裂纹扩展判断准则

在实际问题中，裂纹往往处于复合变形状态，必须采用复合型断裂准则，试验采用最大周向应力准则。基本假设是:(1)开裂方向是最大周向应力σθmax所在的方向;(2)该方向的周向应力达到临界值时，裂纹初始扩展。

裂纹尖端的周向应力为［9］:

式中:r为径向坐标，θ为角坐标，KⅠ和KⅡ分别为Ⅰ和Ⅱ型的应力强度因子。

裂纹扩展角度θ0可由式(1)对θ的一阶偏导得到，注意到当裂纹为纯Ⅰ型时，K=0，θ=Ⅱ00，可得

当沿θ0方向的周向应力达到σθC时，裂纹失稳扩展，即

临界值σθC一般由Ⅰ型开裂条件给出，即σθC= σθ(KⅠC，θ，θ0)，KⅠC为材料的断裂韧度。将 θ0代入式(1)和式(2)，得到裂纹失稳扩展条件为:

2 裂纹扩展模拟方法

2.1 裂纹自动扩展程序的设计思路及流程

断裂力学分析是ABAQUS的一个重要的功能模块，可以准确地计算出结构中裂纹的相关断裂力学参数，但却不具备模拟在未知路径下裂纹扩展的能力，裂纹每扩展一步，都需要将草图创建、裂纹定义、网格划分、计算提交、计算结果提取等这些步骤重复一遍，这就导致了人工工作量十分巨大。通过ABAQUS提供的二次开发接口，利用PYTHON语言，编写裂纹自动扩展的程序模块，极大的提高了分析效率。

本研究的裂纹自动扩展程序主要有以下部分组成:应力分析及裂纹扩展检测、计算分析判断裂纹扩展方向、网格重划分、计算结果写入指定数据文件。程序的核心思路是使用PYTHON控制seam(裂纹)尺寸，然后移动partition(裂尖区域划分)和网格。图1为程序的主流程图。

图1 裂纹自动扩展程序流程图Fig.1 Flow chart of crack automatic propagation program

带裂纹的有限元模型的前处理和有限元计算分别由ABAQUS/CAE和ABAQUS/Standard完成。裂纹判断模块，将根据计算结果，利用最大周向应力准则判断裂纹是否扩展并得出扩展角度。网格重划分模块将根据前一步模型中裂纹尖端的位置，通过给定的扩展增量和扩展角度，计算新裂纹尖端的位置，增加扩展方向上的裂纹，划分裂尖区域，并重新划分整个分析模型的网格，最后提交计算。这一过程会不断重复，直到分析结果满足一定的条件，从而实现对裂纹扩展的仿真。

在裂纹尖端，应力出现奇异点，为了确保得到的结果真实可靠，必须在裂尖放置奇异性单元。在ABAQUS中进行断裂计算时，要在裂纹的尖端布置一圈1/4奇异性单元，用来模拟应力应变的奇异性。引入的奇异性单元，可以看作是把四边形八结点等参单元的一条边上的3个结点重合，使该边退化为一个点，将与裂纹尖端相邻的两条边上的中点移动到距离裂纹尖端1/4处而形成［10］。理论和实际都表明，这种奇异单元具有的奇异性，可以满足裂纹尖端应力和应变奇异性的要求。因此，可在裂尖布置相应的1/4奇异性单元，如图2a所示，裂尖网格采用图2b所示的方式划分。

图2 裂尖网格划分示意图Fig.2 Diagram of crack tip mesh

2.2 程序验证

计算模型如图3所示，矩形板的2个圆孔的圆周上分别有一条斜裂纹，对这2条裂纹在远场拉伸应力作用下的扩展情况进行模拟。文献［8］利用扩展有限元法对上述问题进行了分析，本节将把分析结果与文献［8］中的结果进行对比，以此来验证程序的可靠性与准确性。

该裂纹板是边长为127mm的正方形板，两圆孔位于板的中心位置，两圆心相距25.4mm，圆孔半径为2mm，两条裂纹的初始长度均为2.54mm，裂纹与水平线的初始角度均为θ=45°，裂纹板的上下两表面均施加大小为13790 Pa的拉应力。材料的弹性模量和泊松比分别为E=689.5 MPa和 ν=0.3。文献［8］中裂纹的扩展增量为Δa=1.27mm，本研究分别模拟了 Δa=1.02、1.27、1.52mm情况下裂纹的扩展过程以及对应的断裂力学参数。为了便于与文献给出的结果进行对比，均采用左侧裂纹的值。

图3 双孔边裂纹板Fig.3 Plate with cracks emanating from two hole

图4是利用裂纹自动扩展程序模拟得到的裂纹扩展路径及与文献［8］的对比。可以看到，利用本研究的方法得到的裂纹扩展路径与文献给出的结果基本吻合，表明了程序的正确性。

图5是应力强度因子KⅠ和KⅡ随裂纹长度变化的规律。对于KⅠ，除了中间3步的结果与参考值有较大变化外，其它扩展步得到的结果都与文献结果吻合得很好，KⅡ在裂纹的整个扩展过程中与文献给出的值都基本符合，误差在可接受范围内。

图6是裂纹在不同步长下的扩展轨迹，由裂纹在每一步扩展过程中裂尖的坐标得出。可以看出，只要裂纹扩展的步长值在合理的范围之内，裂纹的扩展轨迹就不会发生较大偏差。因此，在分析实际问题时，当结构的尺寸相对于裂纹较大时，可以适当增加裂纹扩展的步长，这样能在不影响最终结果准确性的前提下减少运算量，节省运算时间。

图4 双孔边裂纹板裂纹扩展路径模拟结果Fig.4 Simulation result of a plate with cracks emanating from two holes

图5 应力强度因子变化规律Fig.5 Variation of stress intensity factor

图6 不同步长下的裂纹扩展轨迹Fig.6 Path of crack propagation under different step

3 无人机机翼肋板裂纹扩展模拟

3.1 机翼肋板结构分析建模与载荷确定

翼肋的形式多种多样，有构架式、腹板式、整体式等。不论使用何种形式，都是用来保持机翼外形、承受局部气动载荷。翼肋的厚度远小于其长度和宽度，而且与机翼上安置翼肋处的翼型在外形上基本一致。因此，可将翼肋简化为二维翼型平板结构进行分析［11］。采用的翼型为NACA－2412，肋板中间的孔为翼梁所在位置。

根据文献［12］的气动分析结果，可将翼肋所受的载荷分段施加，这样更符合实际情况(图7)。具体的载荷为:在翼肋上表面从左至右分5段施加均布载荷，大小依次为 0.12、0.096、0.06、0.096、0.12 MPa;翼肋下表面载荷同样处理，载荷大小依次为0.18、0.156、0.12、0.156、0.18MPa。计算时翼梁所在的位置采用固定约束。肋板材料为铝合金，参数为:E=72 GPa，ρ=2780 kg/m3，μ =0.33。

裂纹取在肋板与梁的连接处，如图8所示，其中θ为裂纹与x轴 (与肋板弦线平行)的夹角。本文模拟了初始夹角θ不同时裂纹的扩展形态，计算了每扩展一步对应的应力强度因子和裂纹扩展角。计算时设定裂纹的初始长度为3mm，裂纹扩展增量为0.2mm。

图7 简化翼肋的有限元模型和载荷模型Fig.7 Simplified finite element model of the wing rib and load model

图8 裂纹位置Fig.8 Location of crack

3.2 结果分析

分别对θ=－70°～70°的裂纹扩展情况进行模拟，图9是不同初始角度下裂纹的扩展形态。可以看到，裂纹的扩展轨迹可以分为两种情况。当裂纹初始角度θ=－70°～40°时，裂纹最终向着翼肋的下表面扩展;当裂纹初始角度θ=－50°～70°时，裂纹沿着初始方向向前扩展，不发生方向的改变。

图10与图11是裂纹初始角度θ=－70°～40°时裂纹断裂力学参量的变化规律。对于应力强度因子KⅠ，除θ=30°和θ=40°外，KⅠ值最终都以大致相同的速率增长，而对于 θ=30°和 θ=40°，KⅠ在经过较长的扩展步之后，其值也呈现与其他情况相同的变化趋势，且θ越大，对应的KⅠ值也越大。对于应力强度因子KⅡ，除了KⅡ值从第16步开始基本不再起伏变化，而其他的裂纹，从第34扩展步开始，其KⅡ值也在0上下小幅度变化。对于裂纹扩展角，θ=40°的裂纹扩展角从第34步开始稳定在0°，裂纹扩展方向基本不再变化，对于其他的裂纹，从第16步开始裂纹扩展的方向就不再变化。

图9 不同初始角度裂纹的扩展形态Fig.9 Crack propagation morphology of different initial crack angle

图10 应力强度因子变化规律 (θ=－70°～40°)Fig.10 Variation of stress intensity factor(θ=－70°～40°)

图11 裂纹扩展角变化规律 (θ=－70°～40°)Fig.11 Variation of crack propagation angle(θ=－70°～40°)

图12是裂纹初始角度时θ=50°～70°裂纹应力强度因子的变化规律。随着裂纹的增长，KⅠ的代数值呈下降趋势，KⅡ的值不断增加。对于θ=50°和θ=60°的裂纹，其应力强度因子的变化趋势都较平缓，值的增长或减小速率大致相同;而对于θ=70°的裂纹，应力强度因子的值则呈指数变化，迅速地增长或减小。

4 结论

1)通过ABAQUS提供的二次开发接口，利用PYTHON语言，编写裂纹自动扩展的程序，实现了对裂纹扩展的模拟，并与相关文献结果进行比较，验证了该模拟方法的有效性和准确性。

2)裂纹的扩展形态、扩展过程中的应力强度因子和裂纹扩展角的变化规律，与裂纹的初始角度有密切的关系，总体可分为两种情况:一种是裂纹沿着初始方向扩展，不发生偏转;另一种是裂纹在扩展的初始阶段发生较大的偏转，随后扩展方向基本不再变化。

图12 应力强度因子变化规律(θ=50°～70°)Fig.12 Variation of stress intensity factor(θ=50°～70°)

3)研究结果对研究机翼肋板裂纹的止裂以及裂纹扩展速率等有一定的参考价值。

［1］秦博，王蕾.无人机发展综述［J］.飞航导弹，2002(8):4－10.

［2］吴圣川，吴玉程.断裂分析及CAE软件的现状与发展［J］.计算机辅助工程，2011，20(1):1－2.

［3］Shephard M，Yehia N，Burd G，et al.Automatic crack propagation tracking［J］.Computers and Structures，1985，20(1):211－223.

［4］Klein P，Gao H.Crack nucleation and growth as strain localization in a virtual－bond continuum ［J］. Engineering Fracture Mechanics，1998，61(1):21－48.

［5］王慰军.基于ABAQUS的裂纹扩展仿真软件及应用［D］.杭州:浙江大学，2006:17－19.

［6］Ortiz M，Pandolfi A.用于三维裂纹扩展分析的有限变形不可逆内聚元［J］.力学进展，2008，38(5):630－640.

［7］娄路亮，曾攀，聂蕾.裂纹扩展的无网格数值模拟方法［J］.航空材料学报，2001，21(3):51－56.

［8］Moёs N，Dolbow J，Belytschko T.A finite element method for crack growth without remeshing［J］.International Journal for Numerical Methods in Engineering，1999，46:131－150.

［9］程靳，赵树山.断裂力学［M］.北京:科学出版社，1996:30－31.

［10］姜薇.基于ABAQUS的船舶典型结构裂纹扩展模拟及分析研究［D］.武汉:华中科技大学，2011:14－15.

［11］季武强.飞机机翼结构拓扑优化方法研究及其实现［D］.沈阳:沈阳航空航天大学，2013:39－40.

［12］郭远帅.某无人机机翼断裂力学分析［D］.南昌:南昌航空大学，2013:37－41.