Abaqus-Umat子程序理论基础与代码编写

本期内容将详细讲解Abaqus二次开发之非线性强化Umat代码，通过这次的学习，只要是应用Mises屈服法则，率无关的任何非线性强化均可适用。

模型介绍：

弹性模量20GPa,泊松比0.3，延续上一节Umat的模型，硬化曲线方程：

$\sigma_{s}=200+400*(\varepsilon_{p})^{0.5}$

材料设置：

编写流程：

计算弹性雅克比矩阵、更新弹性试应力、读取状态变量（弹性应变和塑性应变以及等效塑性应变）、编写计算强化系数和等效Mises应力的子程序（方便随时调用）、屈服后计算等效塑性增量、计算当前应力、弹性应变、塑性应变、调用Mises子程序、计算弹塑性雅克比矩阵、更新状态变量。

本次代码中用到了13个状态变量，1~6个为弹性应变，7~12个为塑性应变，13为等效塑性应变。

代码讲解： 部分代码如下（篇幅有限）：

SUBROUTINE UMAT(STRESS,STATEV,DDSDDE,SSE,SPD,SCD,
     1  RPL,DDSDDT,DRPLDE,DRPLDT,STRAN,DSTRAN,
     2  TIME,DTIME,TEMP,DTEMP,PREDEF,DPRED,MATERL,NDI,NSHR,NTENS,
     3  NSTATV,PROPS,NPROPS,COORDS,DROT,PNEWDT,CELENT,
     4  DFGRD0,DFGRD1,NOEL,NPT,KSLAY,KSPT,KSTEP,KINC)
C
        INCLUDE &#39;ABA_PARAM.INC&#39;
C
        CHARACTER*80 MATERL
        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),
     4  DFGRD0(3,3),DFGRD1(3,3)
        double precision, DIMENSION(NTENS)::deviatoricStress
        double precision misesEqualStress,hydroStaticPressure,W
     1   ,ONESY,DEQPL,SYIELD
C
        DIMENSION EELAS(6),EPLAS(6),FLOW(6)
        PARAMETER (ONE=1.0D0,TWO=2.0D0,THREE=3.0D0,SIX=6.0D0)
        DATA NEWTON,TOLER/1000,1.D-6/
        INTEGER K1,K2,cn

C     PROPS(1) - E
C     PROPS(2) - NU
C     PROPS(3) - SYIELD
C     CALLS AHARD FOR CURVE OF SYIELD VS. PEEQ

C     ELASTIC PROPERTIES
C
        EMOD=PROPS(1)
        EMU=PROPS(2)
        EG2=EMOD/(ONE+EMU)
        EG=EG2/TWO
        EG3=THREE*EG
        ELAM=EMOD*EMU/((ONE+EMU)*(ONE-TWO*EMU))
C
C     ELASTIC STIFFNESS
C
        DO  K1=1,NTENS
        DO  K2=1,NTENS
        DDSDDE(K2,K1)=0.0
        enddo
        enddo
C
        DO K1=1,NDI
        DO K2=1,NDI
        DDSDDE(K2,K1)=ELAM
        enddo
        DDSDDE(K1,K1)=EG2+ELAM
        enddo
        DO K1=NDI+1,NTENS
        DDSDDE(K1,K1)=EG
        enddo
C
C    CALCULATE STRESS FROM ELASTIC STRAINS
C
        DO K1=1,NTENS
        DO K2=1,NTENS
        STRESS(K2)=STRESS(K2)+DDSDDE(K2,K1)*DSTRAN(K1)
        enddo
        enddo
C
C    RECOVER ELASTIC AND PLASTIC STRAINS
C
        DO K1=1,NTENS
        EELAS(K1)=STATEV(K1)+DSTRAN(K1)
        EPLAS(K1)=STATEV(K1+NTENS)
        enddo
        EQPLAS=STATEV(1+2*NTENS)
C
C       HARDENING CURVE, GET YIELD STRESS
C
        CALL AHARD(SYIEL0,HARD,EQPLAS,PROPS,NPROPS)
C
C       MISES STRESS
C
        call GetMisesEqualStress(STRESS,NTENS,NDI,misesEqualStress
     1  ,deviatoricStress,hydroStaticPressure)
C
C       DETERMINE IF ACTIVELY YIELDING
C
        IF (misesEqualStress.GT.(1.0+TOLER)*SYIEL0) THEN
C
        ONESY=ONE/misesEqualStress
        DO K1=1,NDI
            FLOW(K1)=ONESY*(STRESS(K1)-hydroStaticPressure)
......