ABAQUS用户子程序当用到某个用户子程序时,用户所关心的主要有两方面:一是ABAQUS提供的用户子程序的接口参数。
有些参数是ABAQUS传到用户子程序中的,例如SUBROUTINE DLOAD中的KSTEP,KINC,COORDS;有些是需要用户自己定义的,例如F。
二是ABAQUS何时调用该用户子程序,对于不同的用户子程序ABAQUS调用的时间是不同的。
有些是在每个STEP的开始,有的是STEP结尾,有的是在每个INCREMENT的开始等等。
当ABAQUS调用用户子程序是,都会把当前的STEP和INCREMENT利用用户子程序的两个实参KSTEP和KINC传给用户子程序,用户可编个小程序把它们输出到外部文件中,这样对ABAQUS何时调用该用户子程序就会有更深的了解。
(子程序中很重要的就是要知道由abaqus提供的那些参量的意义,如下)首先介绍几个子程序:一.SUBROUTINE DLOAD(F,KSTEP,KINC,TIME,NOEL,NPT,LAYER,KSPT,COORDS, JLTYP,SNAME)参数:1.F为用户定义的是每个积分点所作用的荷载的大小;2.KSTEP,KINC为ABAQUS传到用户子程序当前的STEP和INCREMENT值;3.TIME(1),TIME(2)为当前STEP TIME和INCREMENT TIME的值;4.NOEL,NPT为积分点所在单元的编号和积分点的编号;5.COORDS为当前积分点的坐标;6.除F外,所有参数的值都是ABAQUS传到用户子程序中的。
功能:1.荷载可以被定义为积分点坐标、时间、单元编号和单元节点编号的函数。
2.用户可以从其他程序的结果文件中进行相关操作来定义积分点F的大小。
例1:这个例子在每个积分点施加的荷载不仅是坐标的函数,而且是随STEP变化而变化的。
SUBROUTINE DLOAD(P,KSTEP,KINC,TIME,NOEL,NPT,LAYER,KSPT,COORDS, 1 JLTYP,SNAME)INCLUDE 'ABA_PARAM.INC' CDIMENSION TIME(2),COORDS(3)CHARACTER*80 SNAMEPARAMETER (PLOAD=100.E4)IF (KSTEP.EQ.1) THEN !当STEP=1时的荷载大小P=PLOADELSE IF (KSTEP.EQ.2) THEN !当STEP=2时的荷载大小P=COORDS(1)*PLOAD !施加在积分点的荷载P是坐标的函数ELSE IF (KSTEP.EQ.3) THEN !当STEP=3时的荷载大小P=COORDS(1)**2*PLOADELSE IF (KSTEP.EQ.4) THEN !当STEP=4时的荷载大小P=COORDS(1)**3*PLOADELSE IF (KSTEP.EQ.5) THEN !当STEP=5时的荷载大小P=COORDS(1)**4*PLOADEND IFRETURNENDUMAT 子程序具有强大的功能,使用UMAT 子程序:(1) 可以定义材料的本构关系,使用ABAQUS 材料库中没有包含的材料进行计算,扩充程序功能。
(2) 几乎可以用于力学行为分析的任何分析过程,几乎可以把用户材料属性赋予ABAQUS 中的任何单元;(3) 必须在UMAT 中提供材料本构模型的雅可比(Jacobian)矩阵,即应力增量对应变增量的变化率。
(4) 可以和用户子程序“USDFLD”联合使用,通过“USDFLD”重新定义单元每一物质点上传递到UMAT 中场变量的数值。
由于主程序与UMAT 之间存在数据传递,甚至共用一些变量,因此必须遵守有关UMAT 的书写格式,UMAT 中常用的变量在文件开头予以定义,通常格式为:SUBROUTINE UMA T(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,KSTEP,KINC)INCLUDE 'ABA_PARAM.INC'CHARACTER*80 CMNAMEDIMENSION STRESS(NTENS),STATEV(NSTA TV),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)user coding to define DDSDDE, STRESS, STA TEV, SSE, SPD, SCDand, if necessary, RPL, DDSDDT, DRPLDE, DRPLDT, PNEWDTRETURNENDUMAT 中的应力矩阵、应变矩阵以及矩阵DDSDDE ,DDSDDT ,DRPLDE 等,都是直接分量存储在前,剪切分量存储在后。
直接分量有NDI 个,剪切分量有NSHR 个。
各分量之间的顺序根据单元自由度的不同有一些差异,所以编写UMAT 时要考虑到所使用单元的类别。
下面对UMAT 中用到的一些变量进行说明:DDSDDE(NTENS,NTENS)是一个NTENS 维的方阵,称作雅可比矩阵,,是应力的增量,是应变的增量,DDSDDE(I,J)表示增量步结束时第J 个应变分量的改变引起的第I 个应力分量的变化。
通常雅可比是一个对称矩阵,除非在“*USER MATERIAL”语句中加入了“UNSYMM”参数。
STRESS (NTENS)应力张量矩阵,对应NDI 个直接分量和NSHR 个剪切分量。
在增量步的开始,应力张量矩阵中的数值通过UMAT 和主程序之间的接口传递到UMAT 中,在增量步的结束UMAT 将对应力张量矩阵更新。
对于包含刚体转动的有限应变问题,一个增量步调用UMAT 之前就已经对应力张量的进行了刚体转动,因此在UMA T 中只需处理应力张量的共旋部分。
UMAT 中应力张量的度量为柯西(真实)应力。
STATEV (NSTA TEV)用于存储状态变量的矩阵,在增量步开始时将数值传递到UMAT 中。
也可在子程序USDFLD或UEXPAN 中先更新数据,然后增量步开始时将更新后的数据传递到UMA T 中。
在增量步的结束必须更新状态变量矩阵中的数据。
和应力张量矩阵不同的是:对于有限应变问题,除了材料本构行为引起的数据更新以外,状态变量矩阵中的任何矢量或者张量都必须通过旋转来考虑材料的刚体运动。
状态变量矩阵的维数,等于关键字“*DEPV AR”定义的数值。
状态变量矩阵的维数通过ABAQUS 输入文件中的关键字“*DEPV AR”定义,关键字下面数据行的数值即为状态变量矩阵的维数。
材料常数的个数,等于关键字“*USER MATERIAL”中“CONSTANTS”常数设定的值。
PROPS (NPROPS)材料常数矩阵,矩阵中元素的数值对应于关键字“*USER MA TERIAL”下面的数据行。
SSE ,SPD ,SCD分别定义每一增量步的弹性应变能,塑性耗散和蠕变耗散。
它们对计算结果没有影响,仅仅作为能量输出。
其他变量:STRAN( NTENS) :应变矩阵;DSTRAN( NTENS) :应变增量矩阵;DTIME :增量步的时间增量;NDI :直接应力分量的个数;NSHR :剪切应力分量的个数;NTENS :总应力分量的个数,NTENS NDI NSHR = + 。
使用UMAT 时需要注意单元的沙漏控制刚度和横向剪切刚度。
通常减缩积分单元的沙漏控制刚度和板、壳、梁单元的横向剪切刚度是通过材料属性中的弹性性质定义的。
这些刚度基于材料初始剪切模量的值,通常在材料定义中通过“*ELASTIC”选项定义。
但是使用UMAT 的时候,ABAQUS 对程序输入文件进行预处理的时候得不到剪切模量的数值。
所以这时候用户必须使用“*HOURGLASS STIFFNESS”选项来定义具有沙漏模式的单元的沙漏控制刚度,使用“*TRANSVERSE SHEAR STIFFNESS”选项来定义板、壳、梁单元的横向剪切刚度。
几个关于子程序的问题及相应解答Q: 本人在用umat作本构模型时,*static,1,500,0.000001,0.1 此时要求的增量步很多,即每次增量要很小,*static1,500 时,在弹性向塑性过度时,出现错误,增量过大,出现尖点.?A: YOU CAN TRY AS FOLLOWS:*STEP,EXTRAPOLATION=NO,INC=2000000*STATIC0.001,500.0,0.00001,0.1。
Q: 在abaqus中,如果采用umat,利用自己的本构,如何让abaqus明白这种材料的弹塑性应变,也就是说,如何让程序返回弹性应变与塑性应变,好在output中输出,我曾想用最笨地方法,在uvarm中定义输出,利用getvrm获取材料点的值,但无法获取增量应力,材料常数等,研究了帮助中的例子,umatmst3.inp,umatmst3.for,他采用mises J2 流动理论,我在output history 显示他已进入塑性状态,但他的PE仍然为0!!?A: 用uvar( )勉强成功。
Q: 偶在umat中调用求主应力函数CALL SPRINC(STRESS,PS,LSTR,NDI,NSHR)后,存储主应力得数组PS中各个主应力排列顺序是什么?PS1>PS2>PS3 ?PS1<PS2<PS3 ?PS1>PS3>PS2 ?A: 第二个。
个人觉得:umat实现自己的本构没有固定的方法,对于不同的本构有可能必须采用不同的方法。
这要靠自己不断地摸索。
有可能一种方法对于简单加载问题还行,但有可能对于复杂问题并不收敛。
最重要一点,就是umat中采用的算法必须consistent.再就是ddsdde必须正确,(如果采用back_Euler 方法等一些算法,ddsdde错误有时不影响结果(对于简单加载问题没有影响,能收敛,),但对于复杂问题不收敛。
uptonow,你这个算法对于Mises,hill,J2,J2d等一类的屈服函数是正确的,但具体的本构还要灵活运用,这我也正学习,正在摸索。