Abaqus材料用户子程序UMAT基础知识及手册例子完整解释

1、为何需要使用用户材料子程序（User-Defined Material, UMAT ）？

很简单，当ABAQUS 没有提供我们需要的材料模型时。所以，在决定自己定义一种新的材料模型之前，最好对ABAQUS 已经提供的模型心中有数，并且尽量使用现有的模型，因为这些模型已经经过详细的验证，并被广泛接受。

UMAT 子程序具有强大的功能，使用UMAT 子程序：

(1)可以定义材料的本构关系，使用ABAQUS 材料库中没有包含的材料进行计算，扩充程序功能。

(2) 几乎可以用于力学行为分析的任何分析过程，几乎可以把用户材料属性赋予ABAQU S 中的任何单元。

(3) 必须在UMAT 中提供材料本构模型的雅可比（Jacobian ）矩阵，即应力增量对应变增量的变化率。

(4) 可以和用户子程序“USDFLD ”联合使用，通过“USDFLD ”重新定义单元每一物质点上传递到UMAT 中场变量的数值。

2、需要哪些基础知识？

先看一下ABAQUS 手册（ABAQUS Analysis User's Manual ）里的一段话：

Warning: The use of this option generally requires considerable expertise(一定的专业知识). The user is cautioned that the implementation （实现） of any realistic constitutive （基本） model requires extensive （广泛的） development and testing. Initial testing on a single eleme nt model with prescribed traction loading （指定拉伸载荷） is strongly recommended. 但这并不意味着非力学专业，或者力学基础知识不很丰富者就只能望洋兴叹，因为我们的任务不是开发一套完整的有限元软件，而只是提供一个描述材料力学性能的本构方程（Constitutive equation ）而已。当然，最基本的一些概念和知识还是要具备的，比如：

应力(stress),应变（strain ）及其分量； volumetric part 和deviatoric part ；模量（modul us ）、泊松比(Poisson’s ratio)、拉梅常数(Lame constant)；矩阵的加减乘除甚至求逆；还有一些高等数学知识如积分、微分等。

3、UMAT 的基本任务？

我们知道，有限元计算（增量方法）的基本问题是： 已知第n 步的结果（应力，应变等）n σ，n ε，然后给出一个应变增量1+n d ε，计算新的应力1+n σ。UMAT 要完成这一计算，并要计算Jacobian 矩阵DDSDDE(I,J) =εσΔ?Δ?/。σΔ是应力增量矩阵（张量或许更合适），εΔ是应变增量矩阵。DDSDDE(I,J) 定义了第J 个应变分量的微小变化对

第I 个应力分量带来的变化。该矩阵只影响收敛速度，不影响计算结果的准确性（当然，不收敛自然得不到结果）。

4、怎样建立自己的材料模型？

本构方程就是描述材料应力应变（增量）关系的数学公式，不是凭空想象出来的，而是根据实验结果作出的合理归纳。比如对弹性材料，实验发现应力和应变同步线性增长，所以用一个简单的数学公式描述。为了解释弹塑性材料的实验现象，又提出了一些弹塑性模型，并用数学公式表示出来。

对各向同性材料（Isotropic material）,经常采用的办法是先研究材料单向应力-应变规律（如单向拉伸、压缩试验），并用一数学公式加以描述，然后把该规律推广到各应力分量。这叫做“泛化“(generalization)。

5、一个完整的例子及解释

由于主程序与UMAT之间存在数据传递，甚至一些公共变量，因此必须遵循有关UMAT的书写格式，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,KSTEP,KINC)

INCLUDE 'ABA_PARAM.INC'

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)

user coding to define DDSDDE, STRESS, STATEV, SSE, SPD, SCD

and, if necessary, RPL, DDSDDT, DRPLDE, DRPLDT, PNEWDT

RETURN

END

COORDS 当前积分点的坐标

DDSDDE ( NTENS NTENS)大小为NTENS×NTENS的Jacobian矩阵(ε

σΔ?

?/)，D

Δ

DSDDE(I,J) 定义了第J个应变分量的微小变化对第I 个应

力分量带来的变化。通常Jacobian矩阵是一个对称矩阵，

除非在“*USER MATERIAL”语句中加入了“UNSYM

M”参数；需要更新

DROT对Finite strain问题，应变应该排除旋转部分，该矩阵提

供了旋转矩阵，详见下面的解释；已知

DSTRAN (NTENS)应变增量

dε，已知

n

1+

DTIME增量步的时间增量dt；已知

KSTEP,KINC 传到用户子程序当前的STEP和INCREMENT值

NDI直接应力、应变个数，对三维问题、轴对称问题自然是3

（11,22,33），平面问题是2(11,22)；已知

NOEL,NPT 积分点所在单元的编号和积分点的编号

NSHR剪切应力、应变个数，三维问题时3(12,13,23)，轴对称问

题是1(12)；已知

NTENS=NDI+ NSHR，总应力分量的个数；已知

PNEWDT可用来控制时间步的变化。如果设置为小于1的数，则程

序放弃当前计算，并用新的时间增量DTIME X PNEWDT

作为新的时间增量计算；这对时间相关的材料如聚合物等

有用；如果设为大余1的数，则下一个增量步加大DTIM

E为DTIME X PNEWDT。可以更新。

PROPS (NPROPS)材料常数数组，如模量啊，粘度系数等等；材料参数的个

数，等于关键词“*USER MATERIAL”中“CONSTANT

S”常数设定的值；矩阵中元素的数值对应于关键词“US

ER MATERIAL”下面的数据行。作为已知量传入；已知SSE,SPD,SCD 分别定义每一增量步的弹性应变能，塑性耗散和蠕变耗

散。它们对计算结果没有影响，仅仅作为能量输出STATEV (NSTATEV)状态变量矩阵，用来保存用户自己定义的一些变量，如累

计塑性应变，粘弹性应变等等。增量步开始时作为已知量

传入，增量步结束应该更新

STRAN (NTENS)当前应变数组

ε，已知

n

STRESS (NTENS)应力张量数组，对应NDI个直接分量和NSHR个剪切分

σ中的数值通过U

量。在增量步的开始，应力张量矩阵n

MAT和主程序之间的接口传递到UMAT中，在增量步的

σ。对于包含刚

结束UMAT将对应力张量矩阵更新为1+n

体转动的有限应变问题，一个增量步调用UMAT之前就

已经对应力张量进行了刚体转动，因此UMAT中只需处

理应力张量的共旋部分。UMAT中应力张量的度量为柯

西（真实）应力。

下面这个UMAT取自ABAQUS手册，是一个用于大变形下的弹塑性材料模型，注意的是这里需要了解J2理论。

SUBROUTINE UMAT(STRESS,STATEV,DDSDDE,SSE,SPD,SCD,RPL,DDSDDT,

1 DRPLDE,DRPLDT,STRAN,DSTRAN,TIME,DTIME,TEMP,DTEMP,PREDEF,DPRED,

2 CMNAME,NDI,NSHR,NTENS,NSTATV,PROPS,NPROPS,COORDS,DROT,

3 PNEWDT,CELENT,DFGRD0,DFGRD1,NOEL,NPT,LAYER,KSPT,KSTEP,KINC)

C

INCLUDE 'ABA_PARAM.INC'

定义了一些相关参数与变量什么，从ABAQUS安装目录下的子文件夹“…\site”中可找到

C

CHARACTER*8 CMNAME

C

DIMENSION STRESS(NTENS),STATEV(NSTATV),DDSDDE(NTENS,NTENS),

1 DDSDDT(NTENS)(应变矩阵),DRPLDE(NTENS),STRAN(NTENS),DSTRAN(NTENS)（应变增量矩阵）,

2 PREDEF(1),DPRED(1),PROPS(NPROPS)(材料常数矩),COORDS(3),DROT(3,3)（旋

转矩阵）,

3 DFGRD0(3,3),DFGRD1(3,3)

声明矩阵的尺寸

C

C LOCAL ARRAYS

C ----------------------------------------------------------------

C EELAS - ELASTIC STRAINS

C EPLAS - PLASTIC STRAINS

C FLOW - DIRECTION OF PLASTIC FLOW

C ----------------------------------------------------------------

C

局部变量，用来暂时保存弹性应变、塑性应变分量以及流动方向

DIMENSION EELAS(6),EPLAS(6),FLOW(6)

C

PARAMETER(ZERO=0.D0,ONE=1.D0,TWO=2.D0,THREE=3.D0,SIX=6.D0,

1 ENUMAX=.4999D0,NEWTON=10,TOLER=1.0D-6)

C

C ----------------------------------------------------------------

C UMAT FOR ISOTROPIC ELASTICITY AN

D ISOTROPIC MISES PLASTICITY

C CANNOT BE USE

D FOR PLAN

E STRESS

C ----------------------------------------------------------------

C PROPS(1) - E

C PROPS(2) - NU

C PROPS(3..) - SYIEL

D AN HARDENING DATA

C CALLS HARDSUB FOR CURVE OF YIEL

D STRESS VS. PLASTIC STRAIN

C ----------------------------------------------------------------

C

C ELASTIC PROPERTIES

C

获取杨氏模量，泊松比，作为已知量由PROPS 向量传入

EMOD=PROPS(1) E

ENU=PROPS(2) ν

EBULK3=EMOD/(ONE-TWO*ENU) 3K )21(3ν?=E

k

EG2=EMOD/(ONE+ENU) 2G )1(2υ+=E

G EG=EG2/TWO G )1(2υ+=E

G EG3=THREE*EG 3G

ELAM=(EBULK3-EG2)/THREE λ 3

)23(G k ?=λ DO K1=1,NTENS

DO K2=1,NTENS

DDSDDE(K1,K2)=ZERO

END DO

END DO 弹性部分，Jacobian 矩阵很容易计算

?????????

???????????+++=G G G J 2G 2G 2G λλλλλλλλλ 注意，在ABAQUS 中，剪切应变采用工程剪切应变的定义j,i i,j ij u u +=γ，所以剪切部分模量是G 而不是2G

C Jacobian 矩阵的程序

C ELASTIC STIFFNESS

C

DO K1=1,NDI

DO K2=1,NDI

DDSDDE(K2,K1)=ELAM

END DO

DDSDDE(K1,K1)=EG2+ELAM

END DO

DO K1=NDI+1,NTENS

DDSDDE(K1,K1)=EG

END DO

C

C RECOVER ELASTIC AN

D PLASTIC STRAINS AND ROTAT

E FORWARD

C ALSO RECOVER EQUIVALENT PLASTIC STRAIN

C

读取弹性应变分量，塑性应变分量，并旋转（调用了ROTSIG ），分别保存在EELAS 和EPLAS 中；

CALL ROTSIG(STATEV( 1),DROT,EELAS,2,NDI,NSHR)

CALL ROTSIG(STATEV(NTENS+1),DROT,EPLAS,2,NDI,NSHR)

读取等效塑性应变

EQPLAS=STATEV(1+2*NTENS)

先假设没有发生塑性流动，按完全弹性变形计算试算应力 .1σσσεJ σΔ+=Δ=Δ+n n

C

C CALCULATE PREDICTOR STRESS AN

D ELASTIC STRAIN

C

DO K1=1,NTENS

DO K2=1,NTENS

STRESS(K2)=STRESS(K2)+DDSDDE(K2,K1)*DSTRAN(K1)

END DO

EELAS(K1)=EELAS(K1)+DSTRAN(K1)弹性应变分量

END DO

C计算Mises应力

C CALCULATE EQUIVALENT VON MISES STRESS

C

SMISES=(STRESS(1)-STRESS(2))**2+(STRESS(2)-STRESS(3))**2

1 +(STRESS(3)-STRESS(1))**2

DO K1=NDI+1,NTENS

SMISES=SMISES+SIX*STRESS(K1)**2

END DO

SMISES=SQRT(SMISES/TWO)

C 根据当前等效塑性应变，调用HARDSUB得到对应的屈服应力

C GET YIEL

D STRESS FROM TH

E SPECIFIED HARDENING CURVE

C

NVALUE=NPROPS/2-1

CALL HARDSUB(SYIEL0,HARD,EQPLAS,PROPS(3),NVALUE)

C

C DETERMINE IF ACTIVELY YIELDING

C 如果Mises应力大余屈服应力，屈服发生，计算流动方向

IF (SMISES.GT.(ONE+TOLER)*SYIEL0) THEN

C

C ACTIVELY YIELDING

C SEPARATE THE HYDROSTATIC FROM THE DEVIATORIC STRESS C CALCULATE THE FLOW DIRECTION

C

SHYDRO=(STRESS(1)+STRESS(2)+STRESS(3))/THREE

DO K1=1,NDI

FLOW(K1)=(STRESS(K1)-SHYDRO)/SMISES

END DO

DO K1=NDI+1,NTENS

FLOW(K1)=STRESS(K1)/SMISES

END DO

C根据J2理论并应用Newton-Rampson方法求得等效塑性应变增量

C SOLVE FOR EQUIVALENT VON MISES STRESS

C AN

D EQUIVALENT PLASTIC STRAIN INCREMENT USING NEWTON ITERATIO N

C

SYIELD=SYIEL0

DEQPL=ZERO

DO KEWTON=1,NEWTON

RHS=SMISES-EG3*DEQPL-SYIELD

DEQPL=DEQPL+RHS/(EG3+HARD)

CALL HARDSUB(SYIELD,HARD,EQPLAS+DEQPL,PROPS(3),NVALUE)

IF(ABS(RHS).LT.TOLER*SYIEL0) GOTO 10

END DO

C

C WRITE WARNING MESSAGE TO THE .MSG FILE

C

WRITE(7,2) NEWTON

2 FORMAT(//,30X,'***WARNING - PLASTICITY ALGORITHM DID NOT ',

1 'CONVERGE AFTER ',I3,' ITERATIONS')

10 CONTINUE

σ，应变分量

C更新应力

n

1+

C UPDATE STRESS, ELASTIC AN

D PLASTIC STRAINS AND

C EQUIVALENT PLASTIC STRAIN

C

DO K1=1,NDI

STRESS(K1)=FLOW(K1)*SYIELD+SHYDRO

EPLAS(K1)=EPLAS(K1)+THREE/TWO*FLOW(K1)*DEQPL

EELAS(K1)=EELAS(K1)-THREE/TWO*FLOW(K1)*DEQPL

END DO

DO K1=NDI+1,NTENS

STRESS(K1)=FLOW(K1)*SYIELD

EPLAS(K1)=EPLAS(K1)+THREE*FLOW(K1)*DEQPL

EELAS(K1)=EELAS(K1)-THREE*FLOW(K1)*DEQPL

END DO

EQPLAS=EQPLAS+DEQPL

C

C CALCULATE PLASTIC DISSIPATION

C

SPD=DEQPL*(SYIEL0+SYIELD)/TWO

C

C 计算塑性变形下的Jacobian矩阵

FORMULATE THE JACOBIAN (MATERIAL TANGENT)

C FIRST CALCULATE EFFECTIVE MODULI

C

EFFG=EG*SYIELD/SMISES

EFFG2=TWO*EFFG

EFFG3=THREE/TWO*EFFG2

EFFLAM=(EBULK3-EFFG2)/THREE

EFFHRD=EG3*HARD/(EG3+HARD)-EFFG3

c...

if (props(7).lt..001) go to 99

c...

DO K1=1,NDI

DO K2=1,NDI

DDSDDE(K2,K1)=EFFLAM

END DO

DDSDDE(K1,K1)=EFFG2+EFFLAM

END DO

DO K1=NDI+1,NTENS

DDSDDE(K1,K1)=EFFG

END DO

DO K1=1,NTENS

DO K2=1,NTENS

DDSDDE(K2,K1)=DDSDDE(K2,K1)+EFFHRD*FLOW(K2)*FLOW(K1) END DO

END DO

c...

99 continue

c...

ENDIF

C将弹性应变，塑性应变分量保存到状态变量中，并传到下一个增量步

C STORE ELASTIC AN

D (EQUIVALENT) PLASTIC STRAINS

C IN STATE VARIABLE ARRAY

C

DO K1=1,NTENS

STATEV(K1)=EELAS(K1)

STATEV(K1+NTENS)=EPLAS(K1)

END DO

STATEV(1+2*NTENS)=EQPLAS

C

RETURN

END

c...

c...子程序，根据等效塑性应变，利用插值的方法得到对应的屈服应力

SUBROUTINE HARDSUB(SYIELD,HARD,EQPLAS,TABLE,NVALUE)

C

INCLUDE 'ABA_PARAM.INC'

C

DIMENSION TABLE(2,NVALUE)

C

PARAMETER(ZERO=0.D0)

C

C SET YIEL

D STRESS TO LAST VALU

E O

F TABLE, HARDENIN

G TO ZERO C

SYIELD=TABLE(1,NVALUE)

HARD=ZERO

C IF MORE THAN ONE ENTRY, SEARCH TABLE

C

IF(NVALUE.GT.1) THEN

DO K1=1,NVALUE-1

EQPL1=TABLE(2,K1+1)

IF(EQPLAS.LT.EQPL1) THEN

EQPL0=TABLE(2,K1)

IF(EQPL1.LE.EQPL0) THEN

WRITE(7,1)

1 FORMAT(//,30X,'***ERROR - PLASTIC STRAIN MUST BE `,

1 `ENTERED IN ASCENDING ORDER')

CALL XIT

ENDIF

C

C CURRENT YIEL

D STRESS AND HARDENING C

DEQPL=EQPL1-EQPL0

SYIEL0=TABLE(1,K1)

SYIEL1=TABLE(1,K1+1)

DSYIEL=SYIEL1-SYIEL0

HARD=DSYIEL/DEQPL

SYIELD=SYIEL0+(EQPLAS-EQPL0)*HARD

GOTO 10

ENDIF

END DO

10 CONTINUE

ENDIF

RETURN

END