许可优化
软件
产品
abaqus是世界上功能最强大的有限元仿真软件之一，尤其是在固体非线性力学领域大有一骑绝尘之势，是固体力学CAE软件之翘楚。然而，再强大的仿真软件也不可能内置好能满足所有用户的解决方案，不可能面面俱到。为了解决部分用户无内置解决方案可用的困境，abaqus官方提供了子程序二次开发功能，便于用户开发自定义的解决方案。以下是常用的子程序功能简介，详细资料可参看abaqus官方帮助文档。
abaqus子程序会在恰当的时候被abaqus主程序调用，用户可以通过fortran 、C、C++三类语言实现子程序的编写。目前国内外CAEer通常使用fortran编写子程序，这是一个历史遗留问题，因为fortran统治着20世纪数值计算领域，abaqus最开始是由fortran编写，而后abaqus采用混合编程开发技术，才有了今天的巨无霸仿真软件。
虽然abaqus子程序支持三种语言，但是官方文档只提供了fortran语言的接口，网络上二次开发资料亦是以fortran为开发语言。由于大学入门语言通常是C语言、Python、java，同学们大多对fortran不了解，所以直接上手子程序有点困难，还需要恶补fortran入门子程序。此外，由于fortran的开发环境非常古老，非常不便于开发大型程序，即使是宇宙最强IDE——vs2022在面对fortran语言时也是显得心有力而气不足，其代码提示纠错功能非常差劲 。再次，fortran是面向过程的语言，其对大型可重复使用代码的包装性非常差。
本文旨在提供abaqus C++子程序的入门级教程，以umat程序为介绍对象，笔者学识有限，文中不免有错误，还望各位学友不吝赐教。
见名知意：umat=user+material即为用户自定义材料。官方文档提供的umat接口如下：
SUBROUTINE UMAT(STRESS,STATEV,DDSDDE,SSE,SPD,SCD,
1 RPL,DDSDDT,DRPLDE,DRPLDT,
2 STRAN,DSTRAN,TIME,DTIME,TEMP,DTEMP,PREDEF,DPRED,CMNAME,
3 NDI,NSHR,NTENS,NSTATV,PROPS,NPROPS,COORDS,DROT,PNEWDT,
4 CELENT,DFGRD0,DFGRD1,NOEL,NPT,LAYER,KSPT,JSTEP,KINC)
C
INCLUDE 'ABA_PARAM.INC'
C
CHARACTER*80 CMNAME
DIMENSION STRESS(NTENS),STATEV(NSTATV),
1 DDSDDE(NTENS,NTENS),DDSDDT(NTENS),DRPLDE(NTENS),
2 STRAN(NTENS),DSTRAN(NTENS),TIME(2),PREDEF(1),DPRED(1),
3 PROPS(NPROPS),COORDS(3),DROT(3,3),DFGRD0(3,3),DFGRD1(3,3),
4 JSTEP(4)
user coding to define DDSDDE, STRESS, STATEV, SSE, SPD, SCD
and, if necessary, RPL, DDSDDT, DRPLDE, DRPLDT, PNEWDT
RETURN
END子程序接口本质上就是一个fortran函数，其传入参数类型和数量是确定的，用户可以通过更新形参内容，完成自定义材料开发的功能。自定义材料的范围很广，本文只涉及各向同性线弹性材料，需要更新的形参如下：
线弹性材料在单轴应力下，其应力-应变曲线一定是一条过原点的直线，满足胡克定律。在多轴应力状态下，应力和应变满足广义胡克定律 ，即应力张量和应变张量符合线性张量方程，可以写成式(1)：
$$ \sigma_{i j}=D_{i j k l} \varepsilon_{k l}\tag{1}$$ 式中 \sigma_{ij} 为应力张量， \varepsilon_{ij} 为应变张量， D_{ijkl} 为刚度张量，为正定张量。
采用Viogit法则 可将式(1)写成如下的矩阵方程：
$$\left\{\begin{array}{c} \sigma_{11} \\ \sigma_{22} \\ \sigma_{33} \\ \sigma_{12} \\ \sigma_{13} \\ \sigma_{23} \end{array}\right\}=\left[\begin{array}{rrrrrr} D_{1111} & D_{1122} & D_{1133} & D_{1112} & D_{1113} & D_{1123} \\ & D_{2222} & D_{2233} & D_{2212} & D_{2213} & D_{2223} \\ & & D_{3333} & D_{3312} & D_{3313} & D_{3323} \\ & & & D_{1212} & D_{1213} & D_{1223} \\ & s y m & & & D_{1313} & D_{1323} \\ & & & & & D_{2323} \end{array}\right\}\left\{\begin{array}{c} \varepsilon_{11} \\ \varepsilon_{22} \\ \varepsilon_{33} \\ \gamma_{12} \\ \gamma_{13} \\ \gamma_{23} \end{array}\right\}\tag{2}$$
当材料满足各向同性特性时，材料在任意方位都表现出相同的本构行为，其应力应变关系可写成：
$$\left\{\begin{array}{c} \varepsilon_{11} \\ \varepsilon_{22} \\ \varepsilon_{33} \\ \gamma_{12} \\ \gamma_{13} \\ \gamma_{23} \end{array}\right\}=\left[\begin{array}{cccccc} 1 / E & -\nu / E & -\nu / E & 0 & 0 & 0 \\ -\nu / E & 1 / E & -\nu / E & 0 & 0 & 0 \\ -\nu / E & -\nu / E & 1 / E & 0 & 0 & 0 \\ 0 & 0 & 0 & 1 / G & 0 & 0 \\ 0 & 0 & 0 & 0 & 1 / G & 0 \\ 0 & 0 & 0 & 0 & 0 & 1 / G \end{array}\right]\left\{\begin{array}{c} \sigma_{11} \\ \sigma_{22} \\ \sigma_{33} \\ \sigma_{12} \\ \sigma_{13} \\ \sigma_{23} \end{array}\right\}\tag{3}$$式中 E 为杨氏模量 ， \nu 为泊松比 ， G 为剪切模量 ，三者满足如下约束条件： $$\left\{ \begin{array}{l} G = E/2(1 + \nu )\\ E > 0,G > 0, - 1 < \nu < 0.5 \end{array} \right.\tag{4}$$
使用C++编写umat子程序时，需要包含omi_for_c.h头文件，该文件存在于安装目录下的PublicInterfaces文件夹中，开发者们不用刻意找该文件的路径，只要电脑的abaqus二次开发环境配置好了，并且你的vs安装C++环境，编译时电脑会自动定位该头文件。
学友们，我已经给你们准备好了基于C++的abaqus umat接口，该接口是根据官网提供的fortran接口翻译过来的，接口所有的形参都是指针变量，共有double、int、char三种指针类型，笔者已利用该接口写过很多程序，大家可以放心使用。
#include<omi_for_c.h>
extern "C" void FOR_NAME(umat, UMAT)(double* stress, double* statev, double* ddsdde, double* sse, double* spd,
double* scd, double* rpl, double* ddsddt, double* drplde, double* drpldt,
double* strain, double* dstrain, double* time, double* dtime, double* temp,
double* dtemp, double* predef, double* dpred, char* cmname, int* ndi,
int* nshr, int* ntens, int* nstatv, double* props, int* nprops,
double* coords, double* drot, double* pnewdt, double* celent, double* dfgrd0,
double* dfgrd1, int* noel, int* npt, int* layer, int* kspt,
int* jstep, int* kinc)
{
//write your umat code
return;
}
在umat代码中，核心任务是要利用应变增量矩阵(dstrain)更新应力矩阵(stress)和雅克比矩阵(ddsdde)。说一千，道一万，不如拿眼看一看，请大家鉴赏如下代码，代码已经通过编译，并且与abaqus内置算法矫正过，值得信赖。
使用方法：将代码copy到任意的空cpp文件中，比如test.cpp，然后使用如下cmd命令编译，系统会自动将你的代码编译成test-std.obj二进制文件，其中-std后缀表示standard隐式求解器环境。
abaqus make library=test.cpp#include <iostream>
#include<omi_for_c.h>
extern "C" void FOR_NAME(umat, UMAT)(double* stress, double* statev, double* ddsdde, double* sse, double* spd,
double* scd, double* rpl, double* ddsddt, double* drplde, double* drpldt,
double* strain, double* dstrain, double* time, double* dtime, double* temp,
double* dtemp, double* predef, double* dpred, char* cmname, int* ndi,
int* nshr, int* ntens, int* nstatv, double* props, int* nprops,
double* coords, double* drot, double* pnewdt, double* celent, double* dfgrd0,
double* dfgrd1, int* noel, int* npt, int* layer, int* kspt,
int* jstep, int* kinc)
{
//get mat constants
if (nprops[0] < 2) return;
double E = props[0],v = props[1];
double G= E / (1 + v) / 2;
//Update stiffness matrix
double cfull[6][6]; memset(cfull, 0, sizeof(cfull));
for (int i = 0; i < 3; i++)
{
for (int j = 0; j < 3; j++)
{
cfull[i][j] = E * v / (1 + v) / (1 - 2 * v);
}
cfull[i][i] = E * (1 - v) / (1 + v) / (1 - 2 * v);
cfull[i+3][i+3] = G;
}
//Update stress
for (int i = 0; i < 6; i++)
{
for (int j = 0; j < 6; j++)
{
stress[j] +=cfull[j][i] * dstrain[i];
}
}
//Update Jacobian matrix
for (int i = 0; i < 6; i++)
{
for (int j = 0; j < 6; j++)
{
ddsdde[6*i+j] = cfull[i][j];
}
}
return;
}
为了验证umat程序的可靠性，笔者建立一个对比模型，采用内置线弹性算法，材料设置如图6所示。两个模型的结果如图7所示，由图可知，两者的Mises云图基本一致，这说明笔者提供的子程序是非常可靠的。
本文介绍了abaqus子程序的原理和开发方法，提供了基于C++语言的umat子程序代码接口，及其编译方法。在此基础上，利用C++ umat 接口实现了线弹性各向同材料本构方程的数值算法，并提供了可供读者使用的源代码。源代码在abaqus平台运行后效果甚佳，与内置算法相比相差无几，这说明笔者提供的算法具有较好的鲁棒性。
各位尊敬的读者，如果你们擅长C++，并fortran了解较少，不妨利用笔者提供的C++子程序接口开展abaqus二次开发工作，如果本文对你们有帮助，请你们多多点赞、收藏，并且分享给你们想分享的人，谢谢各位大佬支持。
技术文档
热门文章
155-2731-8020
