冲击载荷下镁铝合金裂纹动态扩展过程的数值模拟*

郭历伦,钟卫洲,陈忠富,罗景润

(中国工程物理研究院总体工程研究所,四川 绵阳 621999)

冲击载荷下镁铝合金裂纹动态扩展过程的数值模拟*

郭历伦,钟卫洲,陈忠富,罗景润

(中国工程物理研究院总体工程研究所,四川 绵阳 621999)

采用基于黏聚裂纹模型的扩展有限元方法，开展了镁铝合金结构冲击破坏过程的数值模拟研究。通过镁铝合金三点弯曲试样冲击实验，获得了不同子弹撞击速度下试样的冲击破坏模式。在此基础上，建立了实验结构的扩展有限元模型，并采用最大主应力准则，以及含损伤型的本构关系模拟材料的冲击断裂行为。对于裂纹尖端附近区域，采用黏聚裂纹模型模拟裂纹的断裂过程。对子弹速度分别为12.2、15.1、26.3 m/s的3种工况下镁铝合金试样的动态破坏过程进行了数值模拟研究，获得了与实验相一致的断裂模式。计算结果表明，试样以Ⅰ型断裂模式为主，裂纹沿初始预制裂纹方向扩展。当裂纹扩展到一定程度后，在试样韧带区域被撞击端附近，由于应力波及边界效应导致该区域处于复杂应力状态，试样出现复合型断裂模式，裂纹偏离原扩展路径，与本文实验结果相吻合。

固体力学；动态破坏；扩展有限元；镁铝合金；黏聚裂纹模型；三点弯曲实验

结构中裂纹动态失稳扩展的数值模拟是一项具有挑战性的工作。采用传统有限元方法进行模拟时，伴随着裂纹的扩展，需不断地进行网格重构，对三维问题而言，其计算规模是难以接受的[1-2]。在商用有限元软件中，目前已存在部分可模拟裂纹动态扩展的数值方法，例如黏聚单元法(cohesive element)或基于表面的黏聚模型(surface-based cohesive behavior)等。这类方法中，需将黏聚表面(cohesive surfaces)布置在单元界面处，且预先指定裂纹扩展的路径，这对裂纹扩展模拟而言是一项极大的限制，因而传统有限元方法将无法给出裂纹动态扩展过程准确的模拟结果。

针对传统有限元方法进行裂纹扩展模拟时存在的弊端，近年来，提出了用于间断问题模拟的修正有限元方法----扩展有限元法(XFEM)[2-5]。其核心思想是在位移近似函数中增加具有间断特征的附加函数项(富集函数，enrichment functions)来模拟间断面的特性，采用level set 函数描述间断界面位置与几何特征，使间断面的描述能够独立于有限元网格。因此，应用XFEM方法模拟结构冲击破坏时，无须进行网格重构，并可模拟裂纹沿任意路径扩展。

本文中采用扩展有限元方法，开展镁铝合金结构冲击破坏过程的数值模拟研究。基于商用有限元软件Abaqus的基本求解技术，探讨结构冲击破坏的XFEM建模方法；开展镁铝合金标准三点弯曲实验，获得试样在不同撞击速度下的冲击破坏图像；在此基础上，建立镁铝合金结构的扩展有限元模型，并采用基于黏聚裂纹模型的扩展有限元方法进行结构破坏过程的数值模拟分析。通过与实验结果的比较，研究镁铝合金结构的裂纹动态扩展过程及结构动态破坏规律。

1 实验结构及扩展有限元建模

1.1 实验结构及加载条件

镁铝合金实验结构为标准三点弯曲试样，其几何尺寸如图1所示，厚度B=20 mm、高度W=40 mm、跨距S=160 mm、长度L=170 mm、初始裂纹长度a=20 mm。

图1 结构示意图Fig.1 Sketch of specimen

初始裂纹采用线切割和疲劳方法进行预制。利用线切割机对长方体试样沿高度W方向进行切割，当距离预定裂尖位置约3 mm时停止切割。线切割完成后，严格依据国家标准GB2038-91，利用MTS810试验机对试样进行高频加载，实现剩余3 mm的疲劳裂纹预制。采用Hopkinson压杆作为实验的加载系统，通过炮筒内发射出的短杆子弹以一定速度撞击入射杆，再由入射杆对试样进行加载，子弹的速度由激光分束测速装置测得。

1.2 结构的扩展有限元建模

为更好地与实验结果进行比较，网格建模时，使几何模型尽量与实验模型相一致。采用前处理软件Patran 2010建立了实验模型的有限元网格，整个模型共包含40 238个六面体单元和46 736个节点，如图2所示。

图2 实验结构的网格、边条及载荷Fig.2 Finite element model, boundary condition and loads in the experimental system

在Abaqus软件中，动态扩展有限元方法仅采用裂纹面富集函数，未采用裂尖富集函数，裂纹不能终止于单元内部，裂纹一旦扩展，一个时间增量步将扩展一个单元[4]。为避免因单个时间增量步裂纹扩展距离过大而导致计算误差，对裂纹可能的扩展区域进行了局部网格细化，如图3所示。

图3 试样的网格模型Fig.3 Finite element model of specimen

由于试样在制作时，预制裂纹由两部分构成，一部分为线切割裂纹，一部分为疲劳预制裂纹。线切割预制部分的裂纹，其裂纹面间的间隙达0.6 mm，而由疲劳方法所预制的裂纹，裂纹面间的间隙很小。因此，为使网格模型尽量与实际的实验结构相一致，避免模型误差所导致的计算误差，将线切割裂纹在有限元网格中考虑，如图3所示。而疲劳预制裂纹则不在几何模型中考虑，在扩展有限元建模时由level set函数定义，作为扩展有限元分析真正意义上的预制裂纹。

1.2.2 材料模型、载荷与边界条件

子弹和入射杆均为合金钢，由于撞击速度不高，其变形应均在弹性范围内，故取线弹性材料常数，杨氏模量为205 GPa，泊松比0.3，材料密度7 800 kg/m3。镁铝合金的弹性模量为43 GPa，泊松比为0.35，密度为1 780 kg/m3。

根据试样实际的支撑和约束情况，在试样被撞击端正对面的部分节点上施加平行于轴向(z向)的固定位移边界条件，在试样底面施加垂直于底面(y向)的固定位移边界条件。直接在子弹上施加实验时的实际初始速度作为其载荷条件，为模拟子弹的运动过程，在子弹与入射杆两撞击面间预留10 mm的初始距离，边界条件和载荷条件如图2所示。将子弹与入射杆、入射杆与试样之间的接触定义为面面接触，正碰可忽略摩擦效应。故接触行为采用软件默认值，法向硬接触、切向无摩擦。

1.2.3 扩展有限元建模

(3)属性rid表示指代中的照应要素(或照应事件)的顺序编号,属性anaphor表示指代中的照应要素(事件指代标注没有这个属性).

在结构网格及材料模型的基础上，在Abaqus软件中进行扩展有限元建模，定义XFEM裂纹、破坏准则和断裂参数、输出和求解控制等。

在模型的inp文件中通过关键字进行XFEM裂纹定义，将图3中的网格细化区域定义为裂纹扩展区域，XFEM裂纹面间的接触行为采用切向无摩擦、法向“硬接触”。 通过关键字*initial conditions指定裂纹面附近节点上的level set函数值，从而定义XFEM裂纹的初始位置(即试样的疲劳预制裂纹)，裂纹在inp文件中的定义方法如图4(a)所示，在模型中的显示如图4(b)所示。

图4 初始裂纹位置的定义及显示Fig.4 Definition and display of the initial crack

断裂准则在材料模块中定义，主要包括损伤起始准则和损伤演化律，需要给定的参数包括：最大主应力、临界应变能释放率、断裂准则因子的误差系数，以及保证计算收敛的黏聚裂纹黏性系数。由于三点弯曲实验中，Ⅰ型断裂模式占主导地位，因此损伤起始准则采用常用的最大主应力准则，裂纹扩展方向为与最大主应力相垂直的方向。材料手册中给出的镁铝合金的静态破坏强度在260～300 MPa之间，由文献[6]中关于铝合金动态强度与静态强度关系的计算公式，其动态强度为静态强度的1.6～1.8倍，故模型中最大主应力的临界值取为400 MPa，断裂准则因子的误差系数取默认值0.05。

损伤演化律采用基于能量法的BK模型，该模型在定义损伤量时需要给定材料的临界应变能释放率。镁铝合金的静态断裂韧性为30～40 MPa·m1/2，因而模型中动态断裂韧性取29～35 MPa·m1/2。由断裂力学公式[7]，计算获得其临界应变能释放率约为1.71×10-2～2.50×10-2MPa·m，黏聚裂纹黏性系数取为1×10-5。

由于裂纹扩展问题属强间断、强非线性问题，较其他问题而言，迭代计算的收敛难度更大。因此，需要根据间断问题的特征，增大迭代计算中允许的尝试步数和每步允许的迭代次数，以提高求解过程的收敛性。此外，为使软件后处理中能显示裂纹位置及扩展情况，在模型输出定义中选定裂纹面水平集函数PHILSM、裂纹前沿水平集函数PSILSM以及裂纹扩展状态参量STATUSXFEM。

2 结构冲击破坏过程的数值模拟

2.1 实验结果

采用Hopkinson实验系统开展了5次实验研究，获得了试样的破坏图像，实验概况如表1所示。5次实验中，试样均发生了裂纹扩展，除实验1因撞击速度较低试样未完全断开外，其余4次实验的试样均完全断裂，实验后试样的破坏情况如图5所示。

从图5可以看出，试样总体呈I型断裂模式，裂纹沿预制裂纹方向扩展，这与通常的认识一致。但在接近试样被撞击区域3～4 mm处，裂纹扩展过程中出现明显的拐弯现象，如图5(b)～(e)所示，说明在试样被撞击区域附近，由于应力波效应，存在复杂应力状态从而导致该区域呈现复合型断裂模式。从图5(f)可知，镁铝合金试样断面光滑，呈典型的脆性断裂模式。

表1 镁铝合金三点弯曲实验概况Table 1 Overview of the magnesium alloy three-point bending specimen

图5 试样破坏情况Fig.5 Fractured specimen

2.2 计算结果及分析

采用基于黏聚裂纹模型的扩展有限元方法[4-5]，对子弹打击速度分别为12.2、15.1和26.3 m/s的3种工况(实验1、2、4)进行了计算分析，获得了试样在不同子弹打击速度下的动态破坏过程。

图6为子弹速度为12.2 m/s时试样的破坏过程，计算获得了裂纹起裂到止裂的全过程。由图6可以看出，裂尖应力分布呈典型的蝴蝶状，结构止裂时共扩展了9个单元(即9 mm)。实验测得实际最终的裂纹扩展距离约为11 mm，计算结果与实验结果比较接近，如图6所示。

图7为子弹速度为15.1 m/s时试样的破坏过程。试样完全断开，与实验结果一致。裂纹扩展16个单元(即16 mm)后，在离被撞击端仅4 mm处裂纹开始偏离原扩展方向，发生偏离的位置与实验实际的位置非常接近。因此，虽然理论上讲，三点弯曲实验中试样的破坏模式应为Ⅰ型断裂模式。但是，在试样被撞击区域附近，由于应力波及试样的边界效应，导致该区域处于复杂应力状态，从而在该区域产生复合型断裂模式。文献[8]在采用扩展有限元方法研究双悬臂梁裂纹扩展过程时，也发现了类似的对称几何构型在对称载荷作用下裂纹扩展方向发生偏转的现象。

图8为子弹速度为26.3 m/s时试样的破坏过程。可以看出试样完全断开，并且也在接近试样被撞击端时裂纹偏离原扩展方向，计算所得的破坏形态也与实验结果较接近。

需要指出的是，实验中，试件韧带区域可能因应力波效应、局部塑形大变形等因素，导致其应力状态和断裂模式十分复杂。因此，对于韧带区裂纹偏转角度等，计算结果与实验结果间尚存在一定的偏差。

图6 子弹打击速度为12.2 m/s时试样裂纹扩展过程Fig.6 Fracture process at a bullet velocity of 12.2 m/s

图7 子弹打击速度为15.1 m/s时试样裂纹扩展过程Fig.7 Fracture process at a bullet velocity of 15.1 m/s

图8 子弹打击速度为26.3 m/s时试样裂纹扩展过程Fig.8 Fracture process at a bullet velocity of 26.3 m/s

3 结 论

采用基于黏聚裂纹模型的扩展有限元方法，对镁铝合金冲击破坏过程进行了数值模拟研究，获得的主要结论如下：

(1)镁铝合金三点弯曲结构的主要断裂模式为Ⅰ型断裂模式，裂纹沿初始预制裂纹方向扩展。当裂纹扩展到一定程度后，在试样韧带区域被撞击端附近，裂纹偏离原扩展路径，试样出现复合型断裂模式；

(2)基于本文数值模拟及实验结果，对试样被撞击端附近裂纹偏转现象的原因进行了分析，结果表明，由于应力波及边界效应导致试样韧带局部区域处于复杂应力状态，是试样在该区域发生复合型断裂模式的可能原因；

(3)计算获得的不同工况下镁铝合金三点弯曲结构的断裂模式均与实验结果较接近，表明基于黏聚裂纹模型的扩展有限元方法可实现冲击载荷下结构复合型断裂过程的模拟。

[1] 庄茁,柳占立,成斌斌,等.扩展有限单元法[M].北京:清华大学出版社,2012:1-7.

[2] 郭历伦,陈忠富,罗景润,等.扩展有限元方法及应用综述[J].力学季刊,2011,32(4):612-625. Guo Lilun, Chen Zhongfu, Luo Jingrun, et al. A review of the extended finite element method and its applications[J]. Chinese Quarterly of Mechanics, 2011,32(4):612-625.

[3] Belytschko T, Black T. Elastic crack growth in finite elements with minimal remeshing[J]. International Journal for Numerical Methods in Engineering, 1999，45(9):601-620.

[4] Remmers J J C, Borst R D, Needleman A. The simulation of dynamic crack propagation using the cohesive segments method[J]. Journal of the Mechanics and Physics of Solid, 2008,56(1):70-92.

[5] Song J H, Areias P M A, Belytschko T. A Method for dynamic crack and shear band propagation with phantom nodes[J]. International Journal for Numerical Methods in Engineering, 2006,67(6):868-893.

[6] 张守中.爆炸与冲击动力学[M].北京:兵器工业出版社,1993:248-248.

[7] 程靳,赵树山.断裂动力学[M].北京:科学出版社,2006:135-137.

[8] Mohammadi S. Extended finite element method for fracture analysis of structures[M]. Tehran, Iran: Blackwell Pulishing Ltd, 2008:187-188.

(责任编辑 曾月蓉)

本刊关于学术不端的声明

本刊注重原创性的学术成果，拒绝学术不端。学术不端包括剽窃(含观点或文字剽窃)、过度引用、伪造、一稿多投和重复发表等。

近年来，学术不端现象时有发生。学术不端损害了期刊形象和读者利益，极大地扰乱了编辑部的正常工作。为此，本刊郑重声明：稿件一经核实存在学术不端，立即终止稿件流程，在我刊网站公布相关信息(如作者姓名和单位等)，通报相关学术期刊和作者单位，并保留向作者追究各种损失(含经济赔偿等)的权利。

爆炸与冲击

2016年9月25日

Numerical research on dynamic fracture process of magnalium alloy under impact load

Guo Lilun, Zhong Weizhou, Chen Zhongfu, Luo Jingrun

(InstituteofSystemsEngineering,ChinaAcademyofEngineeringPhysics,Mianyang621999,Sichuan,China)

The impact fracture process of the magnalium alloy structure was investigated using the XFEM-based cohesive model. First, by the numerical modeling carried out in abaqus software based on XFEM, the fracture mode of magnalium alloy specimens at different bullet impact velocities were obtained from doing a three-point bending experiment. After this, the impact fracture process of experimental model under three different loads at respectively three bullet impact velocities of 12.2, 15.1 and 26.3 m/s was simulated using the XFEM, and the alloy’s failure pattern was obtained by performing numerical calculation, the results from which are consistent with those obtained from the experimental. The simulation results show that Mode Ⅰ is the major fracture mode of the specimen, and the crack propagates mostly along the initial crack direction. The crack makes a turn at a point 3～4 mm from the impacted part of the specimen, where the fixed fracture mode is dominant. This agrees with both the experimental results presented in this paper and with the calculated results found in the related literature. Finally, the reason for the fixed fracture mode in the specimen was also analyzed in the paper.

solid mechanics; dynamic fracture; XFEM; magnalium alloy; cohesive crack; three-point bending specimen

10.11883/1001-1455(2016)05-0648-07

2015-01-13； < class="emphasis_bold">修回日期：2015-06-01

2015-06-01

中国工程物理研究院科学技术发展基金项目(11302211)

郭历伦(1979— )，男，博士,副研究员,gllcr@163.com。

O347.3 <国标学科代码：13015 class="emphasis_bold"> 国标学科代码：13015 文献标志码：A国标学科代码：13015

A