超高强钢材料碰撞失效行为仿真预测技术研究
关键词：超高强钢;失效;GISSMO;
1 前言
节能与安全是汽车行业一直以来关注的两大话题。近年来，实现汽车燃油经济性目标，整车整备质量持续下降，车身轻量化重要性进一步凸显；另一方面，行业安全法规也在持续加严，对车身结构强度提出了更高的要求[1]。超高强钢材料兼具轻量化、性能与综合应用成本优势，近年来在汽车行业应用范围不断增加。随着汽车行业安全法规不断升级，超高强钢结构件在汽车碰撞过程中需要吸收更大的能量；同时，随着材料强度的上升，其韧性通常会有所下降[2]，因而导致车辆关键结构件在碰撞过程中开裂失效风险显著增加，严重影响车辆被动安全性能。为此，汽车行业普遍采用有限元仿真分析方法来预测超高强钢材料在碰撞工况下的失效行为[3,4,5,6,7]，为车辆结构与用材设计提供优化方向。
本研究介绍了目前超高强钢材料碰撞失效行为预测领域的最新研究进展，并针对目前常用的网格尺寸缩放问题进行了研究，可以为提高超高强钢材料在整车碰撞过程中的失效行为预测精度提供一定的参考。
2 失效模型选择
对于超高强钢等金属材料而言，韧性断裂是其最主要的失效形式，采用基于应变的失效模型可以更好地预测其失效行为。作为汽车碰撞仿真分析领域最主流的求解器，LS-DYNA提供了多种失效模型来应对不同的仿真预测精度需求[8]。
对于各向同性失效行为，*MAT＿24材料本构模型提供了基于恒定塑性应变的简单失效判据。这种简单的失效模型通常仅需要材料单轴拉伸试验数据即可以标定，试验量小，但由于没有考虑应力状态、网格尺寸等诸多关键因素，预测精度有限。*MAT＿123材料本构模型在此基础上新增了最大主应变和减薄率2个失效判据，一定程度上提高了失效行为预测精度。
对于金属材料而言，应变路径是影响失效行为的重要因素之一。*MAT＿15材料本构模型引入了损伤累积值概念，其定义如下所示。当损伤累积值达到上限时，相应的单元发生失效并进行删除。增量式损伤累积方法克服了应变失效判据应变路径依赖的问题，提高了失效行为预测的合理性。
式中，D为损伤累积值；Δεp为当前时间步内塑性应变增量；为当前时间步内对应的断裂极限应变；σ*为静水压力与等效应力比值；为应变率；T为温度；d1～d5为断裂极限准则待定参数。
*MAT＿15材料本构模型中材料断裂极限采用JOHNSON COOK准则，该准则是基于现象的经验型公式，对于金属材料而言，该准则在剪切、双拉等应力区间存在明显的不足。为克服*MAT＿15材料本构模型的局限性，*MAT＿224材料本构模型引入了自定义的断裂极限准则，并考虑了仿真网格尺寸对断裂极限应变的影响，进一步提高了失效行为预测精度，其失效准则定义如下：
式中，f (η,θ)为与应力三轴度η及罗德角θ有关的自定义准则；为与应变率有关的自定义修正项；h(T)为与温度相关的自定义修正项；i(Le,η)为与单元特征长度Le及应力三轴度η有关的自定义修正项。
另一方面，损伤与应力耦合效应也是金属材料失效行为不可忽视的现象之一。钢铁材料在某些加载条件下，随着应变增加，损伤不断累积，材料内部孔洞张大，材料有效承载面积下降，进而导致结构承载能力下降。如不考虑损伤累积过程对材料承载能力的影响，则会低估结构失效风险。为此，*MAT＿81材料本构模型引入了应力衰退现象，当单元应塑性变εp达到临界应变εc时，该单元硬化特性将进行缩放，缩放系数与损伤累积过程呈线性关系，其失效准则定义如下。*MAT＿81材料本构模型考虑了损伤与应力耦合效应，但模型失效准则简单，应用场景有限。
式中，εc为临界塑性应变；为单元修正后的应力值。
为克服不同的材料本构模型在预测失效行为的局限性，Neukamm等[9]针对LS-DYNA求解器开发了广义增量应力状态相关模型（Generalized In⁃cremental Stress State dependent Model,GISSMO）。GISSMO失效模型在MAT＿81基础上进一步引入了自定义失效准则、非线性损伤累积、非线性损伤-应力耦合以及单元尺寸、应变率效应修正等特性，可以更加准确地描述金属材料在变形过程中损伤行为。在GISSMO损伤模型中，单元损伤累积值D及失稳累积值F可表示为：
式中，n为损伤及失稳累积值累积指数；β(Le,η)为与单元特征长度Le及应力三轴度η有关的修正项；εc为当前时间增量步对应的材料失稳极限值。当单元的失稳累积值F达到设定上限时，通过以下方式对单元的进行修正：
式中，m为应力退化指数；Dc为单元达到失稳上限时刻对应的损伤累积值。
由上述的各个失效模型原理可知，GISSMO模型是目前LS-DYNA求解器功能最为完整的失效模型，在对超高强钢材料碰撞失效行为仿真预测精度有较高要求时，推荐使用该模型来表征材料失效行为（表1）。
表1 LS-DYNA各向同性失效模型特性 下载原图
3 断裂极限标定
对于完整的GISSMO失效型，需要标定的主要参数包括断裂极限准则、失稳极限准则、损伤累积指数、应力退化指数、网格修正系数等。其中断裂极限可通过实验方法进行标定，而其余参数不具备明确的物理意义，可以采用仿真方法进行反向标定。
金属材料的断裂极限应变通常与应力三轴度及罗德角这两个表征其受力状态的参数相关。修正Mohr-Coulomb(Modified Mohr-Coulomb,MMC）模型[10]是目前行业内应用最广泛的断裂极限准则之一，其基本公式如下所示：
式中，η为应力三轴度；θ为归一化罗德角；A、n是与材料硬化准则相关的参数；c1、c2是Mohr-Coulomb模型基本参数；c3为断裂极限与罗德角相关性调整参数。在整车碰撞工况有限元仿真建模时，结构件通常采用壳单元进行建模，在厚向应力为0假设下，罗德角θ与应力三轴度η存在对应关系，因而断裂极限曲面εf退化为断裂极限曲线。
为标定材料MMC断裂极限准则中相关的5个待定参数，需要设计了一系列不同加载路径的试验来获取材料在不同应力状态下的极限断裂应变。从应力三轴度稳定性、区分度以及试验重复性角度考虑，推荐采用剪切、单轴拉伸、中心孔拉伸、缺口拉伸、等双拉等5种不同应力状态的加载实验，其试样尺寸如下图1所示。材料断裂极限应变目前多采用数字图像相关法（Digital Image Corre⁃lation,DIC）进行测量，其测量精度与试样设计、散斑质量、加载及DIC拍摄速率均相关。在试验过程中，应合理匹配加载速度与DIC图像拍摄帧率，并在DIC结果后处理时按不同的DIC虚拟网格尺寸输出试样关键区域应变测量结果，为后续网格尺寸修正提供参考。
图1 断裂极限应变试验方案
试验及仿真结果表面，断裂极限应变测试过程中关键区域应变路径会随着变形而发生变化。因此，需要采用加载历程的平均应力三轴度ηavg来等效描述相应的加载过程，其计算方法如公式（8）所示。图2为QP1180超高强钢材料不同加载工况应变路径以及MMC断裂极限准则标定结果。
4 网格尺寸修正
有限元仿真模型网格尺寸是影响材料断裂失效行为的关键因素之一。网格尺寸过大时，无法准确描述零部件变形过程中的局部应力状态与应变梯度；网格尺寸过小时，则会显著提高仿真模型计算时间。为了提高大网格尺寸模型仿真预测精度，GISSMO中引入网格尺寸修正系数β(η,Le)来修正不同尺寸单元对应的断裂极限应变。β(η,Le)为与单元应力状态及特征长度Le相关的修正项，在早期版本的求解器中，通过如下的条件函数进行定义：
式中，β0、β1分别为应力三轴度η=0（纯剪）与η=2/3（等双拉）处的修正系数；k(Le)为单元特征长度对应的修正系数。仿真分析结果表明，金属材料在剪切及双拉应力状态下，失效行为通常对网格尺寸不敏感，因此β0与β1一般建议设置为1。在此基础上，采用单轴拉伸工况来最终确定网格尺寸修正系数k(Le)。
图2 QP1180应变路径及断裂极限准则拟合结果
为验证上述的设定方法在实际应用过程中的准确性，除单轴拉伸工况外，设计了剪切、缺口拉伸、等双拉等五种虚拟试验工况，通过对试样尺寸进行优化设计，确保各工况仿真模型最大网格尺寸覆盖3～4 mm，且应力三轴度覆盖[0,0.66]区间。采用GISSMO失效模型网格修正曲线，并将β0与β1设置为1，各个工况不同网格尺寸仿真力-位移特性如下图所示。仿真结果显示，采用上述的网格尺寸修正设置方法，在剪切(平均应力状态～0.1)、单轴（平均应力状态～0.44）以及等双拉（平均应力状态～0.66）工况下，1～4 mm网格模型失效行为仿真预测结果较为一致，但对于缺口拉伸三种工况（平均应力状态0.45～0.57），不同网格尺寸模型仿真结果差异性较大，大网格尺寸模型失效行为预测显著滞后。与1 mm网格模型相比，各工况4 mm网格模型仿真失效位移最大偏差达到85%（图3）。
图3 默认网格尺寸修正设置下不同网格尺寸仿真结果
为提高不同网格尺寸模型失效行为预测结果一致性，采用如图4所示的自定义应力三轴度修正曲线替代现有的分段线性修正曲线。该设置方法可以在R11.0版本之后的求解器中，通过定义网格尺寸修正曲面来实现。对于R11.0版本之前的求解器，则可以采用多个*MAT＿ADD＿EROSION关键字来等效实现。采用新网格尺寸修正设置方法后，6种加载工况不同网格尺寸仿真力-位移曲线如下图5所示。与默认网格修正设置方法相比，自定义应力三轴度修正曲线方式不同网格尺寸预测结果趋于一致，与1 mm网格模型相比，各工况4 mm网格模型仿真失效位移最大偏差由85%下降至9.8%。
图4 网格尺寸修正系数曲线
图5 自定义应力三轴度网格尺寸修正设置下不同网格尺寸仿真结果
5 结论
a.对比研究了LS-DYNA求解器典型的失效模型原理。GISSMO失效模型综合考虑了非线性加载历程、非线性损伤累积、损伤-应力耦合以及单元尺寸修正等问题，可以较好地描述材料损伤累积过程，在对超高强钢材料碰撞失效行为仿真预测精度有较高要求时，推荐使用该模型来进行表征。
b.材料断裂极限应变是失效模型中最为关键的参数，需要设计一系列不同加载路径的试验来进行测试，并采用加载历程的平均应力三轴度ηavg来描述相应的加载过程。
c.模型网格尺寸是影响失效行为仿真预测结果的关键因素。在采用GISSMO失效模型中网格尺寸修正曲线与剪切/双拉修正系数设置方法时，不同网格尺寸仿真模型对于0.45～0.57应力三轴度区间的失效预测结果存在较大差异。采用新的网格尺寸修正设置方法时，可以有效提高不同网格尺寸模型仿真预测结果一致性。
参考文献
[1] 中保研汽车技术研究院有限公司.中国保险汽车安全指数规程:C-IASI-SM.PI.SOT-A0[S]. 2017.
[2] 蒋浩民,陈新平,蔡宁,等.汽车车身用钢的发展趋势[J].锻压技术, 2018(7):56-61.
[3] 王栋,刘淼,王光耀,等.基于LS-DYNA的热成型钢断裂失效预测研究[J].固体力学学报, 2018, 39(2):197-202.
[4] 许伟,方刚,张钧萍,等.面向汽车碰撞安全的热成形钢断裂失效表征与验证[J].塑性工程学报, 2020, 27(6):121-128.
[5] 赖兴华,王磊,李洁,等.铝型材防撞梁的碰撞断裂失效表征[J].清华大学学报:自然科学版, 2017, 57(5):504-510.
[6] 梁宾,赵岩,赵清江,等.基于Gissmo失效模型的6016铝合金板材断裂行为研究及应用[J].机械工程学报,2019, 55(18):53:62.
[7] HUANG G, ZHU H, SRIRAM S, et al. Fracture Prediction and Correlation of AlSi Hot Stamped Steels with Different Models in LS-DYNA[C]//Proceedings of 13th International LS-DYNA Users Conference, Detroit, USA:LSTC,2014, 1-1:17.
[8] LSTC. LS-DYNA Keyword User’s Manual–Volume II:Material Models[M]. Livermore, 2016.
[9] NEUKAMM F, FEUCHT M, HAUFE A ROLL K. On closing the constitutive gap between forming and crash simulation[C]//Proceedings of 10th International LS-DYNA Users Conference, Detroit, USA:LSTC, 2008:21-31.
[10] BAI Y, WIERZBICKI T. A New Model of Metal Plasticity and Fracture with Pressure and Lode Dependence[J].International Journal of Plasticity, 2008, 24(6):1071-96.
文章来源：汽车工艺与材料
