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

维应力(nten), 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
Plain text iconUmat_mine2.txt 6.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

Recent comments

More comments

Syndicate

Subscribe to Syndicate