参考资料：ABAQUS帮助文档
什么是UHYPER？
在ABAQUS中设置超弹性材料的本构模型时，有两种方法，如果模型适用性比较广泛，ABAQUS有内置的相关模型可以使用，只需要定义一些参数。但如果材料的本构ABAQSU中没有的话，ABAQSU提供了一些接口供使用者自己编写本构模型。
UHYPER就是其中的一种，它适用于各向同性超弹性材料，当然也有适用于各向异性超弹性材料的接口(UANISOHYPER)，后续文章也会涉及到。
UHYPER的功能与特点：
1.可以用来定义各向同性超弹性材料的应变势能；
2.可以定义依赖于场变量或状态变量的材料力学行为；
3.要求超弹性材料的应变能密度函数应根据偏应变不变量定义。
UHYPER如何定义本构模型？
根据连续介质力学的理论，超弹性材料的本构模型一般由能量密度函数给出，即用拉伸(stretch)或者不变量去表示能量密度函数，然后通过一些列转化得到应力应变关系。UHYPER是通过偏应变不变量(deviatoric strain invariant)去表示能量密度函数的。偏应变不变量与不变量以及拉伸(stretch)的转化关系如下图所示。
UHYPER的编程结构
UHYPER编写使用Fortran语言，采用Fortran90的固定格式文件，即后缀名为.for，需要遵从以下结构：
SUBROUTINE UHYPER(BI1,BI2,AJ,U,UI1,UI2,UI3,TEMP,NOEL,
1 CMNAME,INCMPFLAG,NUMSTATEV,STATEV,NUMFIELDV,FIELDV,
2 FIELDVINC,NUMPROPS,PROPS)
C
INCLUDE 'ABA_PARAM.INC'
C
CHARACTER*80 CMNAME
DIMENSION U(2),UI1(3),UI2(6),UI3(6),STATEV(*),FIELDV(*),
2 FIELDVINC(*),PROPS(*)
user coding to define U,UI1,UI2,UI3,STATEV
RETURN
END
下图是对程序架构的解释：
我们要做的是：
1.将目标自由能函数转化为用三个偏应变不变量表示，即用BI1，BI2，AJ表示U；
2.将自由能函数U分别对三个偏应变不变量求一阶、二阶、三阶偏导数，其中一阶有三项，二阶和三阶都是六项，具体见帮助文档;
3.将写好的子程序导入ABAQSU计算。
ABAQUS使用子程序时注意事项
1.分析步打开几何非线性，因为超弹性材料一般都涉及大变形；
2.网格类型必须调成全积分和混合计算，子程序不支持缩减积分计算；
3.如果使用状态变量必须在材料属性中定义Depvar；
4.设置超弹性材料的变量数量与输入的参数必须一一对应。
简单示例
我发现ABAQUS内置超弹性模型中有一种非常简单的模型是Mooney Rivlin模型，它能量密度的表达式如下：
$U=C_{10}(\bar{I}_1-3)+C_{01}(\bar{I}_2-3)+\frac{1}{D_1}(J-1)^2$
其中$U$是能量密度函数，$C_{10}、C_{01}和D_1$是三个材料参数。这个表达式三个偏应变不变量都是分开的，求导非常简单。
$\frac{\partial{U}}{\partial{\bar{I}_1}}=UI1(1)=C_{10}$, $\frac{\partial{U}}{\partial{\bar{I}_2}}=UI1(2)=C_{01}$, $\frac{\partial{U}}{\partial{J}}=UI1(3)=\frac{2}{D_1}(J-1)$
$\frac{\partial^2{U}}{\partial{J^2}}=UI2(3)=\frac{2}{D_1}$
即子程序如下：
subroutine uhyper(bi1,bi2,aj,u,ui1,ui2,ui3,temp,noel,
. cmname,incmpflag,numstatev,statev,
. numfieldv,fieldv,fieldvinc,numprops,props)
include 'aba_param.inc'
character*8 cmname
dimension ui1(3),ui2(6),ui3(6),statev(*),fieldv(*),
. fieldvinc(*),props(*)
c
c10 = props(1)
c01 = props(2)
d1 = props(3)
c
c
u = c10*(bi1-3.)+c01*(bi2-3.)+((aj-1.)**2)/d1
ui1(1) = c10
ui1(2) = c01
ui1(3) = 2./d1*(aj-1.)
ui2(1) = 0.
ui2(2) = 0.
ui2(3) = 2./d1
ui2(4) = 0.
ui2(5) = 0.
ui2(6) = 0.
ui3(1) = 0.
ui3(2) = 0.
ui3(3) = 0.
ui3(4) = 0.
ui3(5) = 0.
ui3(6) = 0.
c 以上为0的项可不写
return
end
程序编好了，接下来是导入ABAQUS计算。
不过在计算之前需要确保子程序没有语法错误，在进行ABAQSU关联操作的时候大家都装了Visual Studio，使用这个编译平台去编译以上的程序如下图：
子程序并不是完整的程序，无法编译成功，如图编译时有两个错误即表示子程序无语法错误，此时就可以进行后续操作了。
首先在ABAQUS中随便建立一个简单的拉伸模型，材料属性设置如下图：
由于子程序中PROPS有三个，所以材料参数数量设置为３。
我这个ABAQUS版本太老了，所以在输入这三个材料参数的时候需要在菜单栏选择Model——Editkeywords中修改参数如图：
在图中所示位置输入10,5,1e-5，代表$C_{10}=10,C_{01}=5,D_1=10^{-5}$
根据连续介质力学的理论，$D_1=10^{-5}$表示材料是不可压缩的。
在进行模拟计算时，由于直接使用不可压缩 的模拟方法计算量比较大，而改成可压缩之后计算量将大大减少。所以在进行复杂模拟计算时经常用可压缩的模式计算并引入压缩系数$D_1$，设置其值非常小使材料为近似不可压缩，与直接使用不可压缩的算法结果基本保持一致，但计算量大大减少。本文的这个模型太简单了，用不可压缩计算量也很小。
后续操作就是在ABAQUS中完成压缩模型的各项定义，记得在分析步中打开几何非线性，并把网格属性中的缩减积分关掉。另外，分析步中的最大步长建议调整为0.1，不然最后得到的应力应变曲线数据不够。
最后就是建立作业，需要在如下图所示的地方选择编好的子程序：
接下来就可以把作业提交了。
我建立的拉伸模型是将正方体块拉伸一倍，计算结果如图所示：
作出应力应变曲线如下图所示：
大家也可以将子程序得到的应力应变曲线与同样参数下的ABAQSU内置模型得到的曲线进行对比，材料参数设置如下图所示：
二者得到的应力应变曲线是一模一样的，这也验证了子程序编写的正确性。
至此，最简单的子程序编写与导入ABAQSU计算已经全部完成。当然我们需要编写的子程序肯定不会这么简单，但程序编写思路以及调试过程都是基本一致的，计算偏导数反而会成为其中最关键以及最难的一环，不过这考验的是大家的数学功底。
计算偏导数的一个小技巧：根据连续介质力学理论，能量密度函数的一阶偏导其实是对应的应力，所以求出$UI1$其实就可以得到应力应变曲线，二阶三阶偏导的作用可能只是辅助收敛（这里我也不是很确定）。
大家有什么问题可以通过qq或者发邮件联系我：191905466@qq.com
