基于LS-DYNA的理想弹塑性材料本构关系模型开发基本方法

李嘉文，李美求，李旭东，罗竞波

(长江大学机械工程学院, 湖北 荆州 434023)

基于LS-DYNA的理想弹塑性材料本构关系模型开发基本方法

李嘉文，李美求，李旭东，罗竞波

(长江大学机械工程学院, 湖北 荆州 434023)

LS-DYNA是一款功能强大的有限元分析计算软件，软件自身含有强大的材料库，更允许用户添加自定义的材料本构关系模型，用于得到更准确的仿真结果。基于该平台，介绍了该软件与用户自定义材料子程序之间的数据交互关系，并利用弹塑性力学相关知识，分别建立了理想弹塑性材料的弹性本构关系与塑性本构关系，最后通过算例验证了所构建材料本构关系的正确性，并编写了基于LS-DYNA平台的理想弹塑性材料本构关系，介绍了开发的基本方法，为应用自编写复杂本构关系材料的工程仿真提供了基础。

LS-DYNA；材料本构关系；二次开发

LS-DYNA是一款功能齐全的几何非线性、材料非线性以及摩擦和接触分离等界面状态非线性有限元数值计算软件。凡是涉及接触-碰撞、爆炸、穿甲与侵彻、应力波传播、金属加工、流固耦合等问题，LS-DYNA都可以进行求解。目前，LS-DYNA在国防军工、岩土工程、土木工程、建筑、汽车等领域中均获得广泛应用[1]。虽然LS-DYNA软件自身拥有较完善的材料库，涵盖了弹性、弹塑性、复合材料、蜂窝材料、丝织物、胶类物、生物肌体、炸药、推进剂、混凝土、土壤、地质、超弹性、橡胶、泡沫、玻璃、粘性流体、刚体等各种材料模型[2]，已经可以满足大部分的工程需要，但在实际工程中，针对特殊情况，用户更希望使用自定义的材料本构关系模型，如特殊的材料屈服条件，用以得到更准确的仿真结果。

2006年，金乾坤在TCK模型的基础上开发了一个用于混凝土或钢筋混凝土碰撞或侵彻的损伤模型，并用有关混凝土侵彻试验结果对模型参数进行了标定[3]。2011年，张安康等指出了采用数值算法求解本构方程组时，半隐式的图形返回算法具有良好的应用性[4]。2015年，庄蔚敏等建立了高强钢热成形损伤-相变本构关系模型,通过编写该力学模型的用户自定义材料子程序，将其应用于LS-DYNA软件的高强钢热成形数值模拟中，该试验验证了本构模型的有效性[5]。2016年，胡志强等二次开发功能，将冰材料模型嵌入LS-DYNA程序，并验证了该模型的准确性和适用性[6]。下面,笔者主要从编写材料本构关系的基础出发，介绍自定义本构关系模型开发的基本方法。

1 用户子程序的使用

图1 UMAT文件工作简化流程

LSTC公司在软件中提供了用于用户自定义二次开发的动态链接库文件(LSDYNA.LIB)并预设了程序接口数组实现数据传递。

1.1UMAT文件工作流程

UMAT文件工作流程分为3步：①用户需要利用以FORTRAN语言编写的UMAT文件，结合LS-DYNA.LIB数据库，在Fortran Compiler中编译得到新的可执行文件LSDYNA.EXE；②使用关键字*MAT_USER_DEFINED_MATERIAL_MODELS定义用户材料参数后生成K文件；③求解前在LS-DYNA MANAGER中选择Solver为生成的LSDYNA.EXE进行求解，实现UMAT文件的使用。UMAT文件简化流程图如图1所示。

1.2软件与UMAT文件的接口参数

用户需要在用F语言[7]编写的UMAT文件开头定义如下语言，作为数据接口，正确实现参数的相互传递，以LS-DYNA970版本为例，软件预设的接口为[1]：

Subroutine umat41 (cm,eps,sig,hisv,dt,capa,etype,time,temp,i,ixs,x,k,j)

不仅限于Subroutine UMAT41，Subroutine umat41～umat50都是用户可以自定义的材料，在LS-DYNA中，允许最多10个UMAT同时参与计算。现介绍其中主要的数组，cm是用户指定的材料参数数组；eps是当前时步的6个应变增量数组；sig是本子程序计算得到的6个应力；hisv为用户自定的历史变量；dt是当前时步的时间步长；capa为壳单元的横向剪应力系数；etype为单元类型；time是当前的计算时间。

在LS-DYNA中，一个正确的UMAT文件的主要任务是定义材料的应力-应变关系[8]，其具体的工作方式是由软件计算得到的应变增量Δε，在UMAT文件中计算得到应力增量Δσ，同时更新当前应力σ，并返回至软件中参与下一个时间步的应变增量计算。用户在编写UMAT文件时需要合理使用以上参数，实现数据在软件与UMAT文件间的正确输入与输出。

1.3K文件关键字定义

在LS-DYNA中，用户需要通过关键字*MAT_USER_DEFINED_MATERIAL_MODELS来定义用户材料参数[1]。其中的主要控制参数如表1所示。

表1主要控制参数及说明

用户需要在K文件中对应参数位置输入合理的数值，完成材料参数的定义。其中部分参数需要和UMAT文件结合使用，用户定义这些量时需要注意，避免出现计算错误。

2 理想弹塑性材料本构关系的编写

理想弹塑性材料在受逐渐增大的力作用时，材料的应力(σ)-应变(ε)行为分为弹性阶段和塑性阶段(屈服阶段)，在材料未进入塑性阶段时，材料可视为线弹性材料，其应力应变关系为：

σ=Eε

(1)

式中，E为弹性模量,MPa。

理想弹塑性材料在进入塑性阶段后，其屈服应力σs不再增加，某种程度上，理想弹塑性材料的屈服就意味着破坏。其应力应变关系为：

σ=σs

(2)

基于以上理论，UMAT需要对材料的状态阶段进行判断，从而选择不同的真实应力σ′计算方式，涉及到用户编程时对材料屈服条件的编写。用户可以根据具体需要编写合适的屈服条件。以金属材料为例，最常用的屈服准则有Tresca屈服条件与Mises屈服条件2种。

Tresca屈服条件可以表示为：

(3)

式中，J2为试应力偏量的第二不变量,MPa；θσ为Lode角。

Mises屈服条件则为：

(4)

用户在UMAT中需要按弹性阶段计算试应力σ*，并以此计算屈服条件中的相关参数，与设定的屈服强度σs进行比较，一旦试应力满足以上屈服条件，则认为材料进入了屈服阶段。下面均以Mises屈服条件为判断条件。

为了计算弹性试应力σ*，需要用到弹性力学中的广义胡克定律[9]：

(5)

为方便程序编写，可利用其增量形式进行表达，如式(6)所示：

(6)

其中，σx、σy、σz、τxy、τyz、τxz为6个应力分量，依次对应UMAT中的sig数组中的6个值,MPa； Δεx、Δεy、Δεz、Δγxy、Δγyz、Δγxz为6个应变增量，依次对应UMAT中的eps数组的6个值，1； Δεm为平均应变增量，Δεm=(Δεx+Δεy+Δεz)/6；σm为平均应力，σm=3KΔεm，K=E/3(1-2ν)为体积弹性模量；ν为材料泊松比；G为材料的剪切模量，MPa。

最终的试应力数组σ*为：

图2 应力更新算法

σ*=σ0+ Δσ

(7)

式中，σ0为当前时间步LS-DYNA输入UMAT的应力值，MPa。

σ′=σ*

(8)

若为假，则表明此时材料已经进入塑性阶段，此时应退回到初始应力状态σ0，即：

σ′=σ0

(9)

随后，用户需要将得到的应力值依次返回至sig数组中，更新应力状态，为下一时间的计算做准备。其实现应力更新过程的算法如图2所示。

3 算例验证

试件为5mm×5mm×5mm立方体，加载方式为下端面全约束，上端面施加速率为1mm/s、方向沿Y轴负方向的速度载荷。材料弹性模量E=200GPa，泊松比ν=0.3，密度ρ=6000kg/m3，屈服应力设置为σs=300MPa。使用ANSYS/LS-DYNA模块进行前处理[10]，使用Solid164单元划分单元模型。如图3所示。编译UMAT文件后导入求解，提取该单元的应力-时间、应变-时间数据，可绘制单元的应力-应变曲线，如图4所示。

图3 试件有限元模型 图4 单元应力-应变曲线

由图4可以看出，当应变ε≤1.5×10-3时，应力从零开始随着应变的增加，应力-应变关系为一定斜率直线，说明在该应变下材料处于弹性阶段，符合式(1)的表达形式。当应变ε>1.5×10-3时，此时的等效应力σ=300MPa为一定值，不随应变的增加而变化，且该定值等于屈服应力，符合式(2)的表达形式。材料的2个阶段的表现均符合上述理论分析结果，由此推断该理想弹塑性材料本构关系的UMAT完全正确。

4 结论

1)介绍了LS-DYNA的用户子程序的文件使用流程，并说明了UMAT文件与LS-DYNA软件间是通过应变增量-应力关系来传递更新数据。

2)从理论上分析并建立了理想弹塑性材料本构关系，材料在不同阶段采用不同的应力算法，给出了构建本构关系的基本的流程，用算例验证了本构关系的正确性，为复杂材料本构关系的编写提供了参考。

[1]白金泽.LS-DYNA3D理论基础与实例分析(附光盘)[M].北京：科学出版社, 2005.

[2]刘成清, 何斌, 陈驰.ANSYS/LS-DYNA工程结构抗震、抗撞击与抗连续倒塌分析[M].北京：中国建筑工业出版社, 2014.

[3]金乾坤.混凝土动态损伤与失效模型[J].兵工学报, 2006, 27(1):10～14.

[4]张安康, 陈士海.LS-DYNA用户自定义材料模型开发与验证[J].计算机应用与软件, 2011, 28(4):71～73.

[5]庄蔚敏, 解东旋, 余天明,等.基于损伤-相变本构模型的高强钢热成形数值模拟分析[J].吉林大学学报(工学版), 2015, 45(4):1206～1212.

[6]胡志强, 高岩, 姚琪.一种理想弹塑性模拟的冰材料本构模型[J].船舶与海洋工程, 2016(1):65～73.

[7]黄晓梅, 张伟林.FORTRAN90程序设计[M].合肥：安徽大学出版社, 2016.

[8]师访.ANSYS二次开发及应用实例详解[M].北京： 中国水利水电出版社, 2012.

[9]戴宏亮.弹塑性力学[M].长沙：湖南大学出版社, 2016.

[10]熊令芳, 胡凡金.ANSYS/LS-DYNA非线性动力分析方法与工程应用[M].北京：中国铁道出版社, 2015.

[编辑] 洪云飞

TG113.25

A

1673-1409(2017)17-0049-04

2017-06-18

中国石油科技创新基金项目计划(2015D-5006-0310)。

李嘉文(1992-)，男，硕士生，现主要从事机械结构强度方面的研究工作，519702602@qq.com。

[引著格式]李嘉文，李美求，李旭东，等.基于LS-DYNA的理想弹塑性材料本构关系模型开发基本方法[J].长江大学学报(自科版)，2017,14(17)：49～52.