基于VB的ANSYS二次开发之超弹性材料模型算法
ANSYS现有版本虽包括了金属、橡胶、Drucker-Prager、混凝土等众多材料模型，但仍无法满足工程计算需要，为了弥补这一不足，ANSYS为用户提供了开发材料模型的接口，即UPFs。通过修改、编译连接相关用户子程序，可以得到各种符合用户需要的材料模型。ANSYS的TB，HYPER命令给用户提供了各种不可压缩和可压缩的超弹性材料模型，比如：Polynomid Form模型、Mooney-Rivlin模型、Neo-Hookean模型、Yeoh模型、Arruda-Boyce模型、Gent模型、Ogden模型、Hyperfoam模型以及Blatz-Ko模型等。但是对于需要使用另外模型的用户，则需要UserHyper用户子程序来编写自己的超弹性材料模型。
UserHyper用户子程序介绍
用户可以使用如下命令调用用户定义的超弹性材料模型：
TB,HYPER,,,,USER
可以使用所有支持超弹性材料的单元。
UserHyper用户子程序输入输出变量说明如下：
*deck.UserHyper USERDISTRIB parallel gal
subroutine UscrHyper(
＆ prophy,irtcomp,nprophy,invar,
＆ potential,plnvDer)
c*********************************************************
c 输入变量
c =================
c prophy (dp,ar(*),i) 材料特性数组
c npoophy (int,sc,i) 材料常数个数
c ivar dp,ar(3) 不变量
c
c
c 输出变量
c =================
c incomp (log,sc,i) 可压缩或不可压缩
c potential dp,sc 应变能
c pInvDer dp,ar(10) 应变能wtr对il,i2,j的导数
c l-应变能wtr对il的导数
c 2-应变能wtr对i2的导数
c 3-应变能wtr对il的二阶导数
c 4-应变能wtr对il和i2的导数
c 5-应变能wtr对i2的导数
c 6-应变能wtr对il和j的导数
c 7-应变能wtr对i2和j的导数
c 8-应变能wtr对j的导数
c 9-应变能wtr对j的二阶导数
模型算法
超弹性理论中一个很重要的量是应变能函数，UserHyper用户子程序的应变能函数适用于基于三个应变不变量的率无关且各向同性超弹性模型。如果用户想要定义率相关、非各向同性或非弹性材料模型，则需要利用前面介绍的usermat用户子程序。ANSYS安装文件夹中的UserHyper.F子程序基于Arruda-Boyce模型。这里将要介绍的是一种ANSYS没有的超弹性材料模型。
在超弹性理论中，应变能函数被划分为偏部分(deviatoric part)和体积部分(volumetric part)，因此一般情况下应变能函数的输入变量是3个修正的应变不变量(分别对应UserHyper用户子程序的输入变量invar(l)、invar(2)、invar(3））：
数值实施
在超弹性理论中，可以将应变能函数分成两部分处理，即偏应变能部分和体积应变能部分。本节使用的模型偏应变能函数采用Horgan-Saccomandi建议的形式：
UserHyper用户子程序的详细代码如下：
*deck.UserHyper USERDISTRIB parallel gal
subroutine UscrHyper(
＆ prophy,irtcomp,nprophy,invar,
＆ potential,plnvDer)
c*********************************************************
c 输入变量
c =================
c prophy (dp,ar(*),i) 材料特性数组
c npoophy (int,sc,i) 材料常数个数
c ivar dp,ar(3) 不变量
c
c
c 输出变量
c =================
c incomp (log,sc,i) 可压缩或不可压缩
c potential dp,sc 应变能
c pInvDer dp,ar(10) 应变能wtr对il,i2,j的导数
c l-应变能wtr对il的导数
c 2-应变能wtr对i2的导数
c 3-应变能wtr对il的二阶导数
c 4-应变能wtr对il和i2的导数
c 5-应变能wtr对i2的导数
c 6-应变能wtr对il和j的导数
c 7-应变能wtr对i2和j的导数
c 8-应变能wtr对j的导数
c 9-应变能wtr对j的二阶导数
c
c**************************************************************
c
c---参数
c
#include “impcom.inc”
DOUBLE PRECISION mu,jm,alpha,bulk
DOUBLE PRECISION il,i2,i3
DOUBLE PRECISION jdcnom, vpotcn
c
c ----变量列
c
INTEGER nprophy
DOUBLE PRECISION prophy(*),invar(*),
& potential,pInvDer(*)
LOGICAL incomp
c
c----局部变量
c
mu =prophy(l)/2.0d0
jm =prophy(2)
bulk =prophy(3)
alpha=prophy(4)
c
i1 = invar<l)
i2 = invar(2)
i3 = invar(3)
c
c-----计算偏应变能
c
jdenom=(jm*jm**2.0d0-jm**2.0d0*il+jm*i2-1.OdO)
potential =-mu*jm*log(jdenom /((jm-1.OdO)* !偏应变能Wd
&(jm - 1.0d0)**2.0d0))
plnvDer(1)= mu*jm*jm**2.0d0/jdenom !计算导数值
pInvDer(2)=-mu*jm**2.0d0/jdenom
pInvDer(3)= mu*jm*jm**2.0d0*jm**2.0d0/jdenom**2.0d0
pInvDer(4)=-mu*jm**2.0d0*jm**2.0d0/jdenom**2.0d0
pInvDcr(5)= mu*jm*jm**2.0d0/jdenom**2.0d0
c
c----计算体积应变能
c
pInvDer(6)= O.OdO
pInvDer(7)= O.OdO
pInvDer(8)= O.OdO
pInvDer(9)= O.OdO
IF(alpha.gt.0.0001d0)THEN !通过参数a判断材料是否可压缩
incomp=.FALSE. !可压缩
vpoten =bulk/alpha**2.0d0*(cosh(alpha*(i3-1.OdO))
!体积皮变能 Wv
& -l.OdO)
potential=potential+vpoten !总应变能
plnvDcr(8)= bulk/alpha**1.OdO*sinh(alpha*(i3-l.OdO))
!计算导数值
pInvDer(9)= bulk*cosh(alpha*(i3-l.OdO))
END IF
c
RETURN
END
生成并调用宏文件
在ANSYS中，宏是包含一系列ansys命令并且后缀为.MAC或.mac的命令文件。宏文件往往记录一系列频繁使用的ansys命令流，实现某种有限元分析或其他算法功能。宏文件在ansys中可以当作定义的ansys命令进行使用，可以带有宏输入参数，也可以有内部变量，同时在宏内部可以直接引用总体变量。除了执行一系列的ansys命令之外，宏还可以调用GUI函数或把值传递给参数。
利用*USE命令调用宏文件，并向宏文件传递参数：
*USE,Name,ARG1,ARG2,ARG3,ARG4,ARG5,ARG6,ARG7,ARG8,ARG9,AR10,AR11,AR12,AR13,AR14,AR15,AR16,AR17,AR18
其中，Name是宏文件名，ARGI到AR18是宏文件用到的参数值。
APDL实现过程
下面为两个简单的橡胶类材料受力分析的实例，目的是与ANSYS自带的Gent模型比较，以便验证前面建立的用户超弹性模型的正确性。通过模拟单轴拉伸试验考察Horgan-Saccomandi偏应变能函数，通过模拟静水压缩考察Bischoff体积应变能函数。
1．单轴拉伸
建立两个SOLID185单元，边界条件完全相同，只是使用的材料不同，如图8-7所示。命令流( \chp8\userhyper\userhyper_uniaxial.inp)如下：
finisb
/clear
！设置参数
STRE1TCH=4.9 !拉伸量参数
！进入前处理嚣
/prep7
et,l，185 ！设置单元类型
r,l
tb,hyper,1,1,4,user ！l号材料为用户材料
tbdata,,100,25,0,0 !a=0，所以是不可压缩材料
tb,hyper,2,1,3,gent !2号材料为ANSYS标准gent
tbdata,,lOO,25.4-3,0 !gent超弹材料参数
！几何建模
block，，1，，1，，1
block,,2,3,,1,,1
esize,,1 !仅划分一个单元
mat,1$vmesh,1 !划分网格
mat,2$vmesh,2
！设置边界条件
nsel,s,loc,x,0
nsel,s,loc,x,2
d,all,ux
nsel,s,loc,y,0
d,all,uy
nsel,s,loc,z,0
d,all,uz
nsel,s,loc,y,1
nsel,r,loc,x,0,1
cp,next,uy,all
MYPARM1-ndnext(0)
d*ndnext(0),uy,STRETCH-l !施加拉伸约束
nsel,s,loc,y,1
nsel,r,loc,x,2,3
cp,next,uy,all
d*ndnext(0),uy,STRETCH-l ！施加拉伸约束
allsel,all
finish
!进入求解器
/solu
antype,static ！静态分析类型
nlgeom,on ！打开大变形开关
outres,all,all ！输出所有计算结果
nsubst,100,1000,50 ！求解子步数设置
solve
save
finish
！1.定义工作文件名及文件标题
!(1)定义工作文件名：
/FILNAME,shiyan_thermal,0
!(2)定义文件标题
/TITLE,shiyan of experiment
!2.定义单元类型及材料属性
！（1）定义单元类型：
/PREP7
!*
ET,1,185
！（2）设置材料属性
!*
MPTEMP,,,,,,,,
MPTEMP,1,0
MPDATA,EX,1,,2e11
MPDATA,PRXY,1,,0.3
MPTEMP,,,,,,,,
MPTEMP,1,0
MPDATA,DENS,1,,7800
!3.建立几何模型
！（1）显示工作平面
WPSYYLE,,,,,,,,1
!(2)建立模型
BLOCK,0,10,0,20,0,2,
/REPLO
wpoff,3,5,0
/REPLO
CYLIND,1, ,0,10,0,360,
/REPLO
wpoff,0,10,0
CYLIND,1, ,0,10,0,360,
wpoff,0,-5,0
CYLIND,1, ,0,10,0,360,
wpoff,4,0,0
CYLIND,1, ,0,10,0,360,
wpoff,0,5,0
CYLIND,1, ,0,10,0,360,
wpoff,0,-10,0
CYLIND,1, ,0,10,0,360,
！（3）体布尔操作
FLST,3,6,6,ORDE,2
FITEM,3,2
FITEM,3,-7
VSBV, 1,P51X
！4.生成有限元模型
/REPLO
FINISH
/AUX12
FINISH
/PREP7
MSHAPE,1,3D
MSHKEY,0
!*
CM,_Y,VOLU
VSEL, , , , 8
CM,_Y1,VOLU
CHKMSH,'VOLU'
CMSEL,S,_Y
!*
VMESH,_Y1
!*
CMDELE,_Y
CMDELE,_Y1
CMDELE,_Y2
!*
！5.加载以及求解
！（1）加载
//REPLO
VPLOT
FINISH
/SOL
FLST,2,1,5,ORDE,1
FITEM,2,3
!*
/GO
DA,P51X,ALL,0
/REPLO
FLST,2,1,4,ORDE,1
FITEM,2,7
/GO
FLST,2,1,5,ORDE,1
FITEM,2,4
/GO
!*
SFA,P51X,1,PRES,1000
!*
/STATUS,SOLU
SOLVE
FINISH
/POST1
!*
PLDISP,2
!*
/EFACET,1
PLNSOL, S,EQV, 0,1.0
参考文献
[1] 师访编.ANSYS二次开发及应用实例讲解[M].中国水利水电出版社.2012.1
[2] 隋丽娜,迟剑,郭立峰编. Visual Basic范例开发大全[M].清华大学出版社.
[3] 张艳.Visual Basic程序设计教程[M].徐州:中国矿业大学出版社,2001.
[4] 龚曙光,谢桂兰. ANSYS操作命令与参数化编程[M].北京:机械工业出版社,2004.
当前暂无评论，小编等你评论哦！
![ANSYS/ABAQUS使用(带孔平板拉伸实例)[初识有限元CAE分析]](https://img.jishulink.com/cimage/0c683776c596fd2f7dde46fd773a666b_cdn.jpg?image_process=resize,fw_576,fh_320,)
