LS-DYNA 的材料模型二次开发过程

LS-DYNA 作为一个大型的通用有限元程序，对于多重非线性的大规模问题的解决具有独特的优势，在实际工程中也得到非常广泛的应用。目前材料库有 300 种材料模型，其中多数都提供二维平面应力和三维应力两个版本。LS-DYNA 提供完整的用户材料模型的开发模板，让用户可以开发自己的材料模型。与一般的隐式算法相比，显式有限元分析的时间步长很小，计算规模大，导致对用户子程序的调用非常频繁。LS-DYNA 为减少子程序的调用，内部采用批处理的方式调用用户子程序，要求一次调用能处理几百个单元，这也为用户子程序实现矢量化计算提供了方便。因此，考虑到大变形，LS-DYNA 对用户子程序的特殊要求也增加了用户开发的复杂度。另外，对于一个对初次接触 LS-DYNA 的用户来说，主程序的执行码不带调试信息，较难在源程序上跟踪调试，加大了二次开发中的程序查错的难度。

本文以一个简单大变形下的线弹性材料模型为例，演示在新的开发环境下的完整的开发，调试和验证过程，包括编译，连接，动态加载，源程序跟踪调试，以及模型验证等环节。

1. 线弹性材料物理模型 

1.1 应力应变关系 

在 LS-DYNA 中应力和应变都是 6 个分量。

应力排列顺序为：

\[\sigma = \begin{Bmatrix} \sigma _x \\ \sigma _y \\ \sigma _z \\ \tau _{xy}\\ \tau _{yz} \\ \tau _{zx} \end{Bmatrix}\]

应变排列顺序为：

\[\varepsilon =\begin{Bmatrix} \varepsilon _x\\ \varepsilon _y\\ \varepsilon _z\\ \gamma _{xy}\\ \gamma _{yz}\\ \gamma_{zx} \end{Bmatrix}\]

线弹性材料方程为 \(\sigma=D \varepsilon \)

\[D=\frac{E}{(1+\nu)(1-2\nu))}\begin{Bmatrix} 1-\nu & \nu & \nu& 0 & 0& 0\\ \nu& 1-\nu & \nu & 0 &0 &0 \\ \nu& \nu&1-\nu & 0 &0 &0 \\ 0&0 &0 &\frac{1-2\nu}{2} & 0 & 0\\ 0&0 &0 &0 & \frac{1-2\nu}{2} &0 \\ 0& 0& 0& 0 &0 & \frac{1-2\nu}{2} \end{Bmatrix}\]

其中 E 和 v 分别为杨氏模量和泊松比。

1.2 验证算例 

本文的验证算例为只有一个 8 节点体单元的模型，如图一所示，长为 \(L=2.0m\)，宽和高均为 \(b=1.0m\)。加载条件为在 X 方向单拉，而在 Y 及 Z 方向的位移为零。上述应力应变关系在小变形的情况下则简化为 

应变为：\[\varepsilon _x=\varepsilon =\frac{u}{L}\] \[\varepsilon_y=0\] \[\varepsilon_z=0\]

应力为：\[\sigma _x=\sigma =\frac{E(1-\nu)}{(1+\nu)(1-2\nu) }\varepsilon\] \[\sigma_y=\frac{\nu}{1-\nu}\sigma\] \[\sigma_z=\frac{\nu}{1-\nu}\sigma\]

2. LS-DYNA 用户子程序开发和调试 

2.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 开发包的目录后，进行如下几个操作步骤：

• 把 umat41 子程序从 dyn21.f 中删除，并复制到另一个新的源文件 umat41.f 

subroutine umat41 (cm,eps,sig,epsp,hsv,dt1,capa,etype,tt, temper,failel,crv,nnpcrv,cma,qmat,elsiz,idele,reject) 

• 编辑开发包中的 Makefile，把 umat41.f 的 obj 文件加到 MY_OBJS 变量中 

MY_OBJS = dyn21.o dyn21b.o init_dyn21.o couple2other_user.o dynrfn_user.o umat41.o 

• 选择编译器 

原来的编译器设置是和主程序一致的，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 的二次开发环境及应用 的“LS-DYNA 用户子程序的编译和连接”一节。