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,
Here is the code
I would also request Dr. Jurgen to give me some insight if possible for this specific code problem.
**
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
