二次开发
LS-DYNA的材料模型二次开发过程
Zhidong Han, Brian Wainscott
Livermore Software Technology Corp.
摘要
上期我们介绍了LS-DYNA 新一代二次开发环境的编译连接过程和新增功能，多个用户子程序可以通过动态连接库的方式同时动态加载。本文是其续篇，将以一个简单的材料模型来演示在新的环境下的一个完整的开发过程，包括编译，连接，动态加载，源程序跟踪调试，以及模型验证等环节。
引言
LS-DYNA 作为一个大型的通用有限元程序，对于多重非线性的大规模问题的解决具有独特的优势，在实际工程中也得到非常广泛的应用。目前材料库有300 种材料模型，其中多数都提供二维平面应力和三维应力两个版本。LS-DYNA 提供完整的用户材料模型的开发模板，让用户可以开发自己的材料模型。与一般的隐式算法相比，显式有限元分析的时间步长很小，计算规模大，导致对用户子程序的调用非常频繁。LS-DYNA 为减少子程序的调用，内部采用批处理的方式调用用户子程序，要求一次调用能处理几百个单元，这也为用户子程序实现矢量化计算提供了方便。因此，考虑到大变形，LS-DYNA 对用户子程序的特殊要求也增加了用户开发的复杂度。另外，对于一个对初次接触LS-DYNA 的用户来说，主程序的执行码不带调试信息，较难在源程序上跟踪调试，加大了二次开发中的程序查错的难度。
本文以一个简单大变形下的线弹性材料模型为例，演示在新的开发环境下的完整的开发，调试和验证过程。
线弹性材料物理模型
1）应力应变关系
在LS-DYNA 中应力和应变都是6 个分量，排列顺序为
其中E和v分别为杨氏模量和泊松比。
2）验证算例
本文的验证算例为只有一个8 节点体单元的模型，如图一所示，长为L=2.0m，宽和高均为b=1.0m。加载条件为在X 方向单拉，而在Y 及Z 方向的位移为零。上述应力应变关系在小变形的情况下则简化为
LS-DYNA 用户子程序开发和调试
1）UMAT 子程序的编译和连接
LS-DYNA 中的用户材料号是从41 号到50 号，对应的第一级用户入口子程序是dyn21.f 中的usrmat，受关键字*MAT_USER_DEFINED_MATERIAL_MODELS 控制。
subroutine usrmat (lft,llt,cm,bqs,capa,eltype,mt,ipt,
. npc,plc,crv,nnpcrv,rcoor,scoor,tcoor,nnm1,nip,ipt_thk)
进入这个子程序后，再根据不同的单元类型选择不同的第二级材料子程序
·urmathn: 体单元的三维材料模型
·urmats: 壳单元的二维平面应力材料模型
·urmatb, urmatd, urmatt: 三种不同的梁单元模型
这三个不同子程序根据各自的单元特点对应力应变进行相应的第二级处理之后，再进入的第三级的用户子程序umat41, umat42, … , umat50。这10 个子程序是标准的串行版本模板，演示不同类型的材料模型。一般的开发过程中，第一级和第二级入口程序都不需要改动,只需要从第三级这10 个子程序模板中选一个较为贴近的开始。本文选择umat41 这个子程序。
进入LS-DYNA 开发包的目录后，进行如下几个操作步骤：
1. 把umat41 子程序从dyn21.f 中删除，并复制到另一个新的源文件umat41.f
subroutine umat41 (cm,eps,sig,epsp,hsv,dt1,capa,etype,tt,
1 temper,failel,crv,nnpcrv,cma,qmat,elsiz,idele,reject)
2. 编辑开发包中的Makefile，把umat41.f的obj文件加到MY_OBJS变量中
MY_OBJS = dyn21.o dyn21b.o init_dyn21.o couple2other_user.o dynrfn_user.o umat41.o
3. 选择编译器
原来的编译器设置是和主程序一致的，LINUX系统一般是INTEL或者PGI的Fortran编译器，这些商业编译器的执行代码一般来说效率比较高。在LS-DYNA新的编译环境下，用户子程序的编译器不要求和主程序一致，本文采用开源的gfortran来演示编译过程。编译环境为：
Linux系统：OpenSUSE LEAP 42.1
编译器：gfortran 4.8.5
MPI: platformmpi Community Edition 9.1.2
( http://www-03.ibm.com/systems/platformcomputing/products/mpi/ )
将Makefile中的编译变量设置为
MY_FLAG = -g -fPIC -fcray-pointer -I/opt/platform_mpi/include
FC = /usr/bin/gfortran
LD = /usr/bin/gfortran -shared
export MPI_F77 := /usr/bin/gfortran
MY_TARGET = gnu.so
其中-g是让编译的用户模块带源程序的调试跟踪信息。这些变量的详细解释请参阅上期的“LS-DYNA用户子程序的编译和连接”一节。
4.用make命令编译，生成gnu.so，就完成了编辑和连接。
2）UMAT子程序的调用
上面编译好的gnu.so可以做为开发好的用户模块配合模型使用。这个模块和LS-DYNA主执行程序是分开的，即使将来LS-DYNA主程序的版本升级也不影响这个模块。调用的方法是在模型的.k文件里面加入三行
*MODULE_LOAD
myumat41
gnu.so
其中：第一行是关键字，第二行是这个模块在这个模型的id，第三行是这个模块的编译后文件。然后就可以按照原来的方法执行LSDYNA主程序就可以了。这个关键字有很多匹配规则，详见上期的“LS-DYNA用户子程序的动态连接库的调用”一节。本文演示的是MPP版本的主程序，单个单元模型只能用一个CPU来运行：
mppdyna i=demo.k
3）UMAT子程序的跟踪调试
当子程序运行遇到问题的时候，简单直接方法是用打印命令，与LS-DYNA的59号文件对应的是message文件，对于MPP程序，每个CPU都有一个自己的message文件，因此打印方法不容易混乱，也很方便。比如：
WRITE(59,*)’sig=’,sig(1),sig(2),sig(3)
WRITE(59,*)’hsv=’,hsv(1),hsv(2)
有些情况下，还是要进入到源程序里面，在源程序上进行跟踪调试。本文以gdb为例，启动调试程序，进行以下步骤：
1. 调入主程序
gdb mppdyna
2. 设置断点
b umat41
(此时会显示umat41 不存在，可能要用set breakpoint pending on 激活在调用时补设断点)
3. 运行程序
r i=demo.k
4. 程序在进入umat41 后，就会停下来。比如，打印变量
p cm(1)
这是* MAT_USER_DEFINED_MATERIAL_MODELS 卡片输入的第一个材料常数P1。
材料模型的验证
根据前面定义的物理模型，利用LS-Prepost 建立一个有限元模型，包含一个八节点的实体单元，长度L=2m，宽和高为b=1m。并用LS-Prepost 施加所有的边界条件和给定位移以及材料属性后，有限元模型如图二所示：
所有的节点都固定，只有2，4，6，8 节点在X 方向有指定速度为1.0m/s。分析时间为1s， 则节点位移和工程应变时间曲线为：
与图三的结果吻合得很好。
材料模型为线弹性，E=150x109 及v=0.25，则应力时间曲线为：
σ(t)=180 x109ε(t)=90x109t
而图四中计算结果给出的最大应力则只有σmax =73 x109Pa
原因是LS-DYNA 中的应变都是真实应变，而不是上面计算的工程应变。真实应变的计算方法为：
得到最大真实应变emax=ln(1.5)=0.405
代入应力计算公式可以得到最大真实应力的理论值是73 x109Pa，与有限元的分析结果完全吻合。同时计算结果也表明，当时间t=0.001, 工程应变ε(0.001)=0.005，仅为千分之五， 而此时的真实应力为σ(0.001)=89.7x109t, 与线性公式相差0.3%，而当时间t=0.2时，工程应变达到百分之十，ε(0.2)=0.1，真实应力为σ(0.2)=171 x109t, 与线性公式相差4.7%。
在开发LS-DYNA 的用户子程序的时候，一定要考虑大变形情况下的本构关系。一般来说， 把线性小变形的本构关系直接放进LS-DYNA，真实应力会有很大的差别。这点与一般的小变形下的本构模型开发有很大的不同，也加大了LS-DYNA 二次开发的难度。
作者简介
*韩志东/Zhidong Han博士1998年毕业于清华大学计算固体力学专业，于2011年加入LSTC。他目前从事材料损伤断裂分析及厚壳单元等方面研发。
