User login

Navigation

You are here

UMAT subroutine for 3d solid elements

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

AttachmentSize
Plain text icon Umat_mine2.txt6.68 KB

Comments

Dear

Dear Afallah, 

I have seen the log file. it has no error as such with the same line of compilation error. One may try running it with one element model.

Subscribe to Comments for "UMAT subroutine for 3d solid elements"

Recent comments

More comments

Syndicate

Subscribe to Syndicate