材料属性千变万化,任何一个数值模拟软件的材料模型库也只能包含常见的材料本构模型,但是随着科技的发展以及极端服役环境对材料性能的要求越来越高,各种新的材料本构被不断开发出来,如何将新的本构模型转化为程序并植入到数值模拟软件中成为越来越重要的研究方向。
为满足用户自定义本构的需求,ABAQUS软件提供了非常方便的用户自定义材料本构的FORTRAN子程序接口UMAT(User-defined MATerial)和VUMAT,分别用于ABAQUS/Standard和Abaqus/Explicit 。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)
C user coding to define DDSDDE, STRESS, STATEV, SSE, SPD, SCD
C and, if necessary, RPL, DDSDDT, DRPLDE, DRPLDT, PNEWDT
RETURN
END
材料本构是一门大学问,阿信才疏学浅不敢妄谈,本文主要是以一个例子说明在UMAT中如何设置随时间变化的杨氏模量,对于已经很熟悉UMAT的朋友这自然是小事一桩,所以只针对刚入门的朋友。
材料杨氏模量服从如下公式:
上式中
为时间,这里以秒为单位。
话不多说,上代码!
具体代码核心部分如下,其余部分请自行补全。
E0 = 210000.0
XNUE = 0.3
E = E0 * (ONE - EXP(-0.009*TIME(1)))
STATEV(1) = E
C
C Set up elasticity matrix
C
EBULK3 = E/(ONE-TWO*XNUE)
EG2 = E/(ONE+XNUE)
EG = EG2/TWO
ELAM = (EBULK3-EG2)/THREE
C
CFull = ZERO
DO K1 = 1, 3
DO K2 = 1, 3
CFull(K2,K1) = ELAM
END DO
CFull(K1,K1) = EG2 + ELAM
END DO
CFull(4,4) = EG
CFull(5,5) = EG
CFull(6,6) = EG
C
C Compute the stress tensor
C
DO I=1, NTENS
DO J=1, NTENS
STRESS(I)=STRESS(I)+CFull(I,J)*DSTRAN(J)
END DO
END DO
C
C Determine jacobian
C
DO I=1, NTENS
DO J=1, NTENS
DDSDDE(I,J)=CFull(I,J)
END DO
END DO
在上述代码中,阿信采用了状态变量SDV1监控随时间变化的杨氏模量。采用单个单元进行子程序验证,单元类型为C3D8R,模型如下。
图1 UMAT验证model
下图为代表杨氏模量的SDV1的云图。
图2 自定义场变量SDV1(E)
杨氏模量与时间关系曲线如下图所示。
图3 杨氏模量与时间关系曲线
总的来说,采用UMAT进行自定义的材料本构灵活性很好,但是这需要更多的理论基础,尤其是涉及大量理论推导的本构模型,不但对数学和力学知识要求很高,另外也对编程能力要求较高。
另外,材料属性随空间变化和随机分布,以后有时间给大家分享吧。
水平有限,如有纰漏请指出。
祝各位平安!
如需转载开白,请留言。