一个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)
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)
C
      DIMENSION EELAS(6),EPLAS(6),FLOW(6),DVECT(4)
     REAL DAMAGE,DAMAGC,DAMAG0,BITA,EPDC,EPD0,EQPLAF,DDAMAG,FFDET
      PARAMETER (ONE=1.0D0,TWO=2.0D0,THREE=3.0D0,SIX=6.0D0)
      DATA NEWTON,TOLER/10,1.D-6/
C      
      IF(NDI.NE.3) THEN
      WRITE(6,1)
    1 FORMAT(//,30X,'***ERROR-THIS UMAT MAY ONLY BE USED FOR','ELEMENTS
     & WITH THREE DIRECT STRESS COMPONENTS')
      ENDIF
      WRITE(6,*) NTENS
   
C
C
C ELASTIC PROPERTIES
C 根据弹性模量，泊松比计算体积模量，剪切模量，拉梅常数
      EMOD=PROPS(1)
      ENU=PROPS(2)
c//避免体积模量无穷大
      IF(ENU.GT.0.4999.AND.ENU.LT.0.5001) ENU=0.499
c//3*体积模量
      EBULK3=EMOD/(ONE-TWO*ENU)                  
      EG2=EMOD/(ONE+ENU)
      EG=EG2/TWO                                 
      EG3=THREE*EG                           
      ELAM=(EBULK3-EG2)/THREE                     
C
C ELASTIC STIFFNESS
C 弹性刚度，先置零，后赋值
      DO 20 K1=1,NTENS
      DO 10 K2=1,NTENS
      DDSDDE(K2,K1)=0.0
   10 CONTINUE
   20 CONTINUE
C
      DO 40 K1=1,NDI
      DO 30 K2=1,NDI
      DDSDDE(K2,K1)=ELAM
   30 CONTINUE
      DDSDDE(K1,K1)=EG2+ELAM
   40 CONTINUE
      DO 50 K1=NDI+1,NTENS
      DDSDDE(K1,K1)=EG
   50 CONTINUE
     
C
C    CALCULATE STRESS FROM ELASTIC STRAINS
C   弹性应变引起的应力
      DO 70 K1=1,NTENS
      DO 60 K2=1,NTENS
      STRESS(K2)=STRESS(K2)+DDSDDE(K2,K1)*DSTRAN(K1)
   60 CONTINUE
   70 CONTINUE
C
C   RECVER ELASTIC AND PLASTIC STRAINS
C  更新弹性和塑性应变，说明STATEV(1~6)是6个方向弹性应变，7~12是塑性应变
C  等效应变EQPLAS是STATEV(13)
      DO 80 K1=1,NTENS
      EELAS(K1)=STATEV(K1)+DSTRAN(K1)
      EPLAS(K1)=STATEV(K1+NTENS)
   80 CONTINUE
      EQPLAS=STATEV(1+2*NTENS)
      WRITE(6,1000) KINC, NOEL, NPT, EQPLAS
1000 FORMAT(//,30X,"INCREMENT",2X,I2,2X,'ELEMENT NUMBER',
     & 2X,I3,2X,'NPT',2X,I1,2X,'CURRENT EQ PLASTIC STRAIN',2X,F12.6)
C
            
      IF(NPROPS.GT.2.AND.PROPS(3).GT.0.0) THEN
    SMISES=(STRESS(1)-STRESS(2))*(STRESS(1)-STRESS(2))+
     &(STRESS(2)-STRESS(3))*(STRESS(2)-STRESS(3))+
     &(STRESS(3)-STRESS(1))*(STRESS(3)-STRESS(1))
   DO 90 K1=NDI+1,NTENS
    SMISES=SMISES+SIX*STRESS(K1)*STRESS(K1)
   90 CONTINUE
       SMISES=SQRT(SMISES/TWO)
C    调用子程序得到屈服应力
      NVALUE=NPROPS/2-1
  CALL AHARD(SYIEL0,HARD,EQPLAS,PROPS(3),NVALUE)
C 判断是否屈服
      IF(SMISES.GT.(1.0+TOLER)*SYIEL0)THEN
C  屈服后的流动方向
    SHYDRO=(STRESS(1)+STRESS(2)+STRESS(3))/THREE
    ONESY=ONE/SMISES
   DO 110 K1=1,NDI
    FLOW(K1)=ONESY*(STRESS(K1)-SHYDRO)
  110 CONTINUE
      DO 120 K1=NDI+1,NTENS
    FLOW(K1)=STRESS(K1)*ONESY
  120 CONTINUE
C   牛顿迭代求等效塑性应变
C
      SYIELD=SYIEL0
    DEQPL=0.0
    DO 130 KEWTON=1,NEWTON
    RHS=SMISES-EG3*DEQPL-SYIELD
    DEQPL=DEQPL+RHS/(EG3+HARD)
   CALL AHARD(SYIELD,HARD,EQPLAS+DEQPL,PROPS(3),NVALUE)
   IF (ABS(RHS).LT.TOLER*SYIEL0) THEN
C  Output incremental equivalent plastic strain
       WRITE(6,2000) KINC,NOEL,NPT,DEQPL
2000  FORMAT(//,30X,"INCREMENT",2X,I2,2X,"ELEMENT NUMBER",2X,I3,2X,
     & 'NPT',2X,I1,2X,"EQ PLASTIC STRAIN INC",2X,F12.6)
       GOTO 140
      ENDIF
  130 CONTINUE
      WRITE(6,2) NEWTON
2    FORMAT(//,30X,'***WARNING-PLASTICITY ALGORITHM DID NOT',
     &'CONVERGE AFTER',I3,'ITERATIONS')
  140 CONTINUE
      EFFHRD=EG3*HARD/(EG3+HARD)
C         
C计算应力更新应变
      DO 150 K1=1,NDI
    STRESS(K1)=FLOW(K1)*SYIELD+SHYDRO
    EPLAS(K1)=EPLAS(K1)+THREE*FLOW(K1)*DEQPL/TWO
    EELAS(K1)=EELAS(K1)-THREE*FLOW(K1)*DEQPL/TWO  
  150 CONTINUE
      EQPLAS=EQPLAS+DEQPL
   SPD=DEQPL*(SYIEL0+SYIELD)/TWO
C更新J矩阵
      EFFG=EG*SYIELD/SMISES
   EFFG2=TWO*EFFG
   EFFG3=THREE*EFFG2/TWO
   EFFLAM=(EBULK3-EFFG2)/THREE
   DO 220 K1=1,NDI
     DO 210 K2=1,NDI
      DDSDDE(K2,K1)=EFFLAM
  210   CONTINUE
        DDSDDE(K1,K1)=EFFG2+EFFLAM
  220 CONTINUE
      DO 230 K1=NDI+1,NTENS
     DDSDDE(K1,K1)=EFFG
  230 CONTINUE
      DO 250 K1=1,NTENS
     DO 240 K2=1,NTENS
     DDSDDE(K2,K1)=DDSDDE(K2,K1)+FLOW(K2)*FLOW(K1)
     &  *(EFFHRD-EFFG3)
  240 CONTINUE
  250 CONTINUE
      ENDIF
   ENDIF
C储存状态变量
      DO 310 K1=1,NTENS
      STATEV(K1)=EELAS(K1)
      STATEV(K1+NTENS)=EPLAS(K1)
  310 CONTINUE
      STATEV(1+2*NTENS)=EQPLAS
C print equivalent plastic strain to .dat file
      WRITE(6,3000) KINC, NOEL,NPT, STATEV(1+2*NTENS)
3000 FORMAT(//,30X,"INCREMENT",2X,I2,2X,'ELEMENT NUMBER',
     & 2X,I3,2X,'NPT',2X,I1,2X,'TOTAL EQ PLASTIC STRAIN',2X,F12.6)
C
      RETURN
      END   
C
C***子程序****ahard***************
      SUBROUTINE AHARD(SYIELD,HARD,EQPLAS,TABLE,NVALUE)
      INCLUDE 'ABA_PARAM.INC'
      DIMENSION TABLE(2,NVALUE)
C 求EQPLAS在哪段斜率内，然后线性叠加求应力，返回SYIELD和斜率HARD
      SYIELD=TABLE(1,NVALUE)  
      HARD=0.0
C
C IF MORE THAN ONE ENTRY,SEARCH TABLE
C
      IF(NVALUE.GT.1) THEN
      DO 10 K1=1,NVALUE-1
      EQPL1=TABLE(2,K1+1)
      IF(EQPLAS.LT.EQPL1) THEN
      EQPL0=TABLE(2,K1)
      IF(EQPL1.LE.EQPL0) THEN
      WRITE(6,1)
    1 FORMAT(//,30X,'***ERROR-PLASTIC STRAIN MUST BE','ENTERED IN ASCEND
     &ING ORDER')
   CALL XIT
      ENDIF
C
C   CURRENT YIELD 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 20
      ENDIF
   10  CONTINUE
   20  CONTINUE
      ENDIF
      RETURN
      END   

同时请教论坛里各位高手，有个问题困扰我很久了。在程序中
1. 在“C    调用子程序得到屈服应力” 处
      NVALUE=NPROPS/2-1
  CALL AHARD(SYIEL0,HARD,EQPLAS,PROPS(3),NVALUE)

中的NVALUE=NPROPS/2-1是为了求什么呢？
2.还有在“C***子程序****ahard***************” 处
这个子程序中TABLE函数中的值在ABAQUS中如何赋值的，我找了很久都不会

还有这个子程序是什么意思呢，不知道应该看什么理论知识才能把这段程序看明白，请各位帮我解释一下啊，谢谢啦！！！

秦淮客

楼主现在的讲解正是我需要的:)

52cae

1. NVALUE=NPROPS/2-1是为了求什么呢？
    回答：UMAT中材料参数PROPS是一维数组形式给出，而在ahard子程序中则把这个一维数组的材料参数转为两行多列的二维TABLE给出的，第一行表示应力，第二行表示对应的塑性应变。所以调用ahard时得通过维度换算式找到这个TABLE的入口地址。
2. ahard子程序的原理
    回答：其实这里采用的就是最简单的分段线性插值的本构模型，先搜索到当前的等效塑性应变值EQPLAS位于TABLE中的两点之间，再线性插值得到对应的SYIELD值。

华夏岩土

52cae 发表于 2012-4-12 17:58
1. NVALUE=NPROPS/2-1是为了求什么呢？
    回答：UMAT中材料参数PROPS是一维数组形式给出，而在ahard子程 ...

多谢指点啊，那就意味着在ABAQUS中的mechanical constant中还要把应力和对就的等效塑性应变也填入表中了吗？那它们填入的顺序是填一个应力，接着填对应的应变吗？

52cae

华夏岩土 发表于 2012-4-13 10:20
多谢指点啊，那就意味着在ABAQUS中的mechanical constant中还要把应力和对就的等效塑性应变也填入表中了 ...

对于这个UMAT是这样的。

nuaalizhen

秦淮客 发表于 2012-4-12 17:53
楼主现在的讲解正是我需要的

文东，你可以看下这个：
http://forum.simwe.com/thread-1088788-1-1.html

hou2012

nuaalizhen 发表于 2013-11-1 21:16
文东，你可以看下这个：
http://forum.simwe.com/thread-1088788-1-1.html

nuaalizhen兄，请教一个问题，此贴楼主提供的这段代码对应的理论应该看哪些东西？我看了陈惠发的弹性和塑性力学，但是看到这段代码的  屈服后的流动方向  这个地方时就看不懂了，不知道那个flow向量对应塑性理论中的那个？求大神帮助

nuaalizhen

hou2012 发表于 2015-1-9 17:17
nuaalizhen兄，请教一个问题，此贴楼主提供的这段代码对应的理论应该看哪些东西？我看了陈惠发的弹性和塑 ...

你好，你可以看下我贴的网址，文档里有参考文献

qsw1990615

:(:(:(:(:(:(