ABAQUS材料子程序学习（20年12月2日）

@ABAQUS材料子程序学习（20年12月2日）

---

前言

继续记录自己学习过程，本文针对《非线性本构关系在ABAQUS中的实现》第三章“黏弹性“本构的学习，UMAT子程序，有VS code编写。与abaqus帮助文档里的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,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)

      DIMENSION DSTRESS(NTENS),KK(NTENS,NTENS)

      E1=PROPS(1) !储存模量（100~800MPa）
      E2=PROPS(2)  !损耗模量（0~100MPa）
      ETA=PROPS(3)  !黏度（0~1000MPa）
      MU=PROPS(4)   !泊松比

      MU_BAR=ETA/(E1+E2)
      LMD=(MU*E1*E2)/((1.0+MU)*(1-2.0*MU)*(E1+E2))
      LMD_BAR=LMD*ETA/E2
      EG=0.5*E1*E2/((1.0+MU)*(E1+E2))
      EG_BAR=EG*ETA/E2

C     应力更新
      D=0.5*DTIME+MU_BAR
      D_1=1.0/D 
      A=EG*DTIME+2.0*EG_BAR
      B=0.5*DTIME*LMD+LMD_BAR

      STRAN_V=0.0
      DSTRAN_V=0.0
      DO I1=1,NDI
        STRAN_V=STRAN_V+STRAN(I1)
        DSTRAN_V=DSTRAN_V+STRAN(I1)
      ENDDO

      DO I1=1,NDI
        X=LMD*STRAN_V+2.0*EG*STRAN(I1)-STRESS(I1)
        DSTRESS(I1)=(A*DSTRAN(I1)+B*DSTRAN_V+X*DTIME)*D_1
        STRESS(I1)=STRESS(I1)+DSTRESS(I1)
      ENDDO
      DO I1=NDI+1,NTENS
        Y=2.0*EG*STRAN(I1)-STRESS(I1)
        DSTRESS(I1)=(A*DSTRAN(I1)+Y*DTIME)*D_1
        STRESS(I1)=STRESS(I1)+DSTRESS(I1)
      ENDDO

C     一致性切线模量
      T1=(A+B)*D_1
      T2=B*D_1
      T3=0.5*A*D_1
      DDSDDE=0.0
      DO I1=1,NDI
        DDSDDE(I1,I1)=T1
      ENDDO
      DO I1=2,NDI
        DO I2=1,I1-1
          DDSDDE(I1,I2)=T2
          DDSDDE(I2,I1)=T2
        ENDDO
      ENDDO
      DO I1=1+NDI,NTENS
        DDSDDE(I1,I1)=T3
      ENDDO


C     状态量更新
      DU=0.0
      DO I1=1,NTENS
        DU=DU+(STRESS(I1)+0.5*DSTRESS(I1))*DSTRAN(I1)     
      ENDDO
      KK=0.0
      DO I1=1,NDI
        KK(I1,I1)=LMD+2.0*EG
      ENDDO
      DO I1=2,NDI
        DO I2=1,I1-1
          KK(I1,I2)=EG
          KK(I2,I1)=EG
        ENDDO
      ENDDO
      DO I1=1+NDI,NTENS
        KK(I1,I1)=EG
      ENDDO
      DUE=0.0
      DO I1=1,NTENS
        DO I2=1,NTENS
          DUE=DUE+(STRAN(I1)+0.5*DSTRAN(I1))*KK(I1,I2)*DSTRAN(I2)     
        ENDDO
      ENDDO
      SSE=SSE+DUE
      SCD=SCD+DU-DUE

      RETURN
      END

结果