Skip to main content

UMAT subroutine for 3d solid elements

Submitted by hermatician on

Dear All,
I am okay in ABAQUS although trying subroutines for the first time.
I am trying to modify some codes available online, but the primary UMAT which I am posting here, is not working.
It shows no error in .log file but the job monitor says "Problem during compilation - D:\Abaqus Work Directory\Trial_user_subroutine\try\Umat_mine.for"

Please help me this, its very urgent.
Thanking you all,
P.S. The code is attached in .txt format
Here is the code


**
        SUBROUTINE UMAT(STRESS, STATEV, DDSDDE, SSE, SPD, SCD, RPL,
     1 DDSDDT, DRPLDE, DRPLDT, STRAN, DSTRAN, TIME, DTIME, TEMP, DTEMP,
     2 PREDEF, DPRED, CMNAME, NDI, NSHR, NTENS, NSTATV, PROPS, NPROPS,
     3 COORDS, DROT, PNEWDT, CELENT, DFGRD0, DFGRD1, NOEL, NPT, LAYER,
     4 KSPT, KSTEP, KINC)
   
        INCLUDE 'ABA_PARAM.INC'
   
        CHARACTER*80 CMNAME
   
        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    LOCAL ARRAYS
C ----------------------------------------------------------------
C  EELAS - ELASTIC STRAINS
C  EPLAS - PLASTIC STRAINS
C  FLOW - DIRECTION OF PLASTIC FLOW
C -------------------------------------------

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

      PARAMETER(ZERO=0.D0, ONE=1.D0, TWO=2.D0, THREE=3.D0, SIX=6.D0,
     & ENUMAX=.4999D0, NEWTON=10, TOLER=1.0D-6)
C
C ----------------------------------------------------------------
C     UMAT FOR ISOTROPIC ELASTICITY AND ISOTROPIC MISES PLASTICITY
C
CANNOT BE USED FOR PLANE STRESS
C ----------------------------------------------------------------
C PROPS(1) - E
C PROPS(2) - NU
C PROPS(3..) - SYIELD AN HARDENING DATA
C
CALLS UHARD FOR CURVE OF YIELD STRESS VS. PLASTIC STRAIN
C ----------------------------------------------------
C ELASTIC PROPERTIES
      EMOD=PROPS(1)
      ENU= PROPS(2)
      EBULK3=EMOD/(ONE-TWO*ENU)
      EG2=EMOD/(ONE+ENU)
      EG=EG2/TWO
      EG3=THREE*EG
      ELAM=(EBULK3-EG2)/THREE


C ELASTIC STIFFNESS
      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 Example 5: Isotropic Hardening Plasticity RECOVER ELASTIC AND PLASTIC STRAINS AND ROTATE FORWARD
C  ALSO RECOVER EQUIVALENT PLASTIC STRAIN
      CALL ROTSIG(STATEV(   1), DROT, EELAS, 2, NDI, NSHR)
      CALL ROTSIG(STATEV(NTENS+1), DROT, EPLAS, 2, NDI, NSHR)
      EQPLAS=STATEV(1+2*NTENS)
C
C
C
CALCULATE PREDICTOR STRESS AND ELASTIC STRAIN
        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 CALCULATE EQUIVALENT VON MISES STRESS
        SMISES=(STRESS(1)-STRESS(2))**2+(STRESS(2)-STRESS(3))**2
     &  +(STRESS(3)-STRESS(1))**2
        DO K1=NDI+1,NTENS
        SMISES=SMISES+SIX*STRESS(K1)**2
        END DO
        SMISES=SQRT(SMISES/TWO)

C
C Example 5: Isotropic Hardening Plasticity GET YIELD STRESS FROM THE SPECIFIED HARDENING CURVE
        NVALUE=NPROPS/2-1
        CALL UHARD(SYIEL0, HARD, EQPLAS, EQPLASRT,TIME,DTIME,TEMP,
     &             DTEMP,NOEL,NPT,LAYER,KSPT,KSTEP,KINC,CMNAME,NSTATV,
     &             STATEV,NUMFIELDV,PREDEF,DPRED,NVALUE,PROPS(3))
     
     
        IF (SMISES.GT.(ONE+TOLER)*SYIEL0) THEN
       
C  ACTIVELY YIELDING SEPARATE THE HYDROSTATIC FROM THE DEVIATORIC STRESS  CALCULATE THE FLOW DIRECTION
        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  SOLVE FOR EQUIVALENT VON MISES STRESS AND EQUIVALENT PLASTIC STRAIN INCREMENT USING NEWTON ITERATION





        SYIELD=SYIEL0
        DEQPL=ZERO
        DO KEWTON=1, NEWTON
        RHS=SMISES-EG3*DEQPL-SYIELD
        DEQPL=DEQPL+RHS/(EG3+HARD(1))
        CALL UHARD(SYIELD,HARD,EQPLAS+DEQPL,EQPLASRT,TIME,DTIME,TEMP,
     1  DTEMP,NOEL,NPT,LAYER,KSPT,KSTEP,KINC,CMNAME,NSTATV,
     2   STATEV,NUMFIELDV,PREDEF,DPRED,NVALUE,PROPS(3))
        IF(ABS(RHS).LT.TOLER*SYIEL0) GOTO 10
        END DO
C        WRITE WARNING MESSAGE TO THE .MSG FILE
       
C   Example 5: Isotropic Hardening Plasticity
        WRITE(7,2) NEWTON
     2  FORMAT(//,30X,'***WARNING - PLASTICITY ALGORITHM DID NOT ',
     1   'CONVERGE AFTER ',I3,' ITERATIONS')
  10      CONTINUE
       
        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
       
        SPD=DEQPL*(SYIEL0+SYIELD)/TWO
        EFFG=EG*SYIELD/SMISES
        EFFG2=TWO*EFFG
        EFFG3=THREE/TWO*EFFG2
        EFFLAM=(EBULK3-EFFG2)/THREE
        EFFHRD=EG3*HARD(1)/(EG3+HARD(1))-EFFG3
        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
        ENDIF
        DO K1=1, NTENS
        STATEV(K1)=EELAS(K1)
        STATEV(K1+NTENS)=EPLAS(K1)
        END DO
        STATEV(1+2*NTENS)=EQPLAS
C
        RETURN
        END


        SUBROUTINE UHARD(SYIELD,HARD,EQPLAS,EQPLASRT,TIME,DTIME,TEMP,
     1  DTEMP,NOEL,NPT,LAYER,KSPT,KSTEP,KINC,
     2  CMNAME,NSTATV,STATEV,NUMFIELDV,
     3  PREDEF,DPRED,NVALUE,TABLE)
        INCLUDE 'ABA_PARAM.INC'
        CHARACTER*80 CMNAME
       
       
        DIMENSION HARD(3),STATEV(NSTATV),TIME(*),
     $          PREDEF(NUMFIELDV),DPRED(*),PROPS(*)
C
C     USES LUDWIG EQUATION: SIGMA_YIELD = SIGMA_0 + K*EQPLAS^N
C                           d(SIGMA)/d(EQPLAS) = K*N*EQPLAS^(N-1)
C     USER INPUT:
     
C    THE MATERIAL UNKNOWNS PROPS 1 2 AND 3 ARE DEFINED IN ABAQUS NOT IN CODE, SO NO NEED TO DEFINE PROPS 1 2 3
C
        PARAMETER(ZERO=0.D0, ONE=1.D0)
C
C     CALCULATE
C
          SYIELD=PROPS(4) + PROPS(5)*EQPLAS**PROPS(6)
          HARD(1)=PROPS(5)*PROPS(6)*EQPLAS**(PROPS(6)-ONE)
          HARD(2)=ZERO
          HARD(3)=ZERO
        GOTO 10
        RETURN
      END

Attachment Size
Umat_mine2.txt 6.68 KB