isotropic-Mises-plasticity-model-modern-form.for

与ABAQUS用户帮助文档中的各向同性mises塑性模型相同，但是编写格式略作修改

      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 'ABA_PARAM.INC'
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)
C
      DIMENSION EELAS(6),EPLAS(6),FLOW(6)
      PARAMETER (ONE=1.0D0,TWO=2.0D0,THREE=3.0D0,SIX=6.0D0)
      DATA NEWTON,TOLER/10,1.D-6/
C
C -----------------------------------------------------------
C     UMAT FOR ISOTROPIC ELASTICITY AND ISOTROPIC PLASTICITY
C     J2 FLOW THEORY
C     CAN NOT BE USED FOR PLANE STRESS
C -----------------------------------------------------------
C     PROPS(1) - E
C     PROPS(2) - NU
C     PROPS(3) - SYIELD
C     CALLS AHARD FOR CURVE OF SYIELD VS. PEEQ
C -----------------------------------------------------------
C
      IF (NDI.NE.3) THEN
         WRITE(6,1)
 1       FORMAT(//,30X,'***ERROR - THIS UMAT MAY ONLY BE USED FOR ',
     1          'ELEMENTS WITH THREE DIRECT STRESS COMPONENTS')
      ENDIF
C
C     ELASTIC PROPERTIES
C
      EMOD=PROPS(1)
      ENU=PROPS(2)
      IF(ENU.GT.0.4999.AND.ENU.LT.0.5001) ENU=0.499
      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 K1=1,NTENS
        DO K2=1,NTENS
           DDSDDE(K2,K1)=0.0
        END DO
      END DO
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    CALCULATE STRESS FROM ELASTIC STRAINS
C
      DO K1=1,NTENS
        DO K2=1,NTENS
           STRESS(K2)=STRESS(K2)+DDSDDE(K2,K1)*DSTRAN(K1)
        END DO
      END DO
C
C    RECOVER ELASTIC AND PLASTIC STRAINS
C    状态向量STATEV：应变(1-6)、塑性应变(7-12)、等效塑性应变(13)
      DO K1=1,NTENS
         EELAS(K1) = STATEV(K1) + DSTRAN(K1)
         EPLAS(K1) = STATEV(K1 + NTENS)
      END DO
      EQPLAS = STATEV(1 + 2*NTENS)
C
C    IF NO YIELD STRESS IS GIVEN, MATERIAL IS TAKEN TO BE ELASTIC
C    如果没有屈服应力，材料选择弹性
      IF(NPROPS.GT.2.AND.PROPS(3).GT.0.0) THEN
C
C       MISES STRESS
C
        SMISES = SQRT(ONE/TWO*(
     1           (STRESS(1)-STRESS(2))**2 +
     1           (STRESS(2)-STRESS(3))**2 +
     1           (STRESS(3)-STRESS(1))**2 + 
     1       SIX*(STRESS(4)**2 +
     1            STRESS(5)**2 +
     1            STRESS(6)**2))) 

C
C       HARDENING CURVE, GET YIELD STRESS
C
        NVALUE=NPROPS/2-1
        CALL AHARD(SYIEL0,HARD,EQPLAS,PROPS(3),NVALUE)
C
C       DETERMINE IF ACTIVELY YIELDING
C
        IF (SMISES.GT.(1.0+TOLER)*SYIEL0) THEN
C
C        FLOW DIRECTION 流动方向
C        
         SHYDRO=(STRESS(1)+STRESS(2)+STRESS(3))/THREE
         FLOW(1) = (STRESS(1) - SHYDRO)/SMISES
         FLOW(2) = (STRESS(2) - SHYDRO)/SMISES
         FLOW(3) = (STRESS(3) - SHYDRO)/SMISES
         FLOW(4) = (STRESS(4))/SMISES
         FLOW(5) = (STRESS(5))/SMISES
         FLOW(6) = (STRESS(6))/SMISES
C
C       SOLVE FOR EQUIV STRESS, NEWTON ITERATION 牛顿迭代法
C
          SYIELD=SYIEL0
          DEQPL=0.0
          DO 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
                EFFHRD=EG3*HARD/(EG3+HARD)
             ELSE
                WRITE(6,2) NEWTON
 2              FORMAT(//,30X,'***WARNING - PLASTICITY ALGORITHM DID NOT ',
     1                 'CONVERGE AFTER ',I3,' ITERATIONS')
             END IF
          END DO

C
C       CALC STRESS AND UPDATE STRAINS 计算应力和更新应变
C
          DO 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
          END DO
C
          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
          SPD=DEQPL*(SYIEL0+SYIELD)/TWO
C
C       JACOBIAN
C
          EFFG=EG*SYIELD/SMISES
          EFFG2=TWO*EFFG
          EFFG3=THREE*EFFG2/TWO
          EFFLAM=(EBULK3-EFFG2)/THREE
          DO K1=1,NDI
             DO K2=1,NDI
                DDSDDE(K2,K1)=EFFLAM
             END DO
             DDSDDE(K1,K1)=EFFG2+EFFLAM
          END DO 
C
          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)+FLOW(K2)*FLOW(K1)
     1                                       *(EFFHRD-EFFG3)
             END DO
          END DO
        ENDIF
      ENDIF
C
C STORE STRAINS 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 AHARD(SYIELD,HARD,EQPLAS,TABLE,NVALUE)
C
      INCLUDE 'ABA_PARAM.INC'
      DIMENSION TABLE(2,NVALUE)
C
C    SET YIELD STRESS TO LAST VALUE OF TABLE, HARDENING TO ZERO
      SYIELD=TABLE(1,NVALUE) ! 最大屈服应力
      HARD=0.0
C
C   IF MORE THAN ONE ENTRY, SEARCH TABLE 如果塑性数据点数大于1
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(6,1)
 1           FORMAT(//,30X,'***ERROR - PLASTIC STRAIN MUST BE ',
     1                        'ENTERED IN ASCENDING 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
            EXIT
           ENDIF
        END DO
      ENDIF !------------------------------------------------------------

      RETURN
      END

      SUBROUTINE AHARD(SYIELD,HARD,EQPLAS,TABLE,NVALUE)
C
      INCLUDE 'ABA_PARAM.INC'
      DIMENSION TABLE(2,NVALUE)
C
C    SET YIELD STRESS TO LAST VALUE OF TABLE, HARDENING TO ZERO
      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 ',
     1                 'ENTERED IN ASCENDING 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