Notifications
Clear all

Problem with UMAT subroutine for isotropic hardening


SharkOon
Posts: 2
Topic starter
(@SharkOon)
New Member
Joined: 3 years ago

Hi all,

I am a beginner when it comes to writing subroutines in abaqus, that is to say I am still working my way through already available codes. So, I am trying to make a subroutine for isotropic hardening work. While I have resolved other more obvious errors, one error that sustains is how the state variables would reset to zero after every iteration (This happens when it enters the plastic region). On furthur enquiry, I found it was because the solution was not converging. However, I am not able to find the cause for the non convergence and how to make it work. Kindly help me on this. Here is the UMAT code for your perusal.

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)
C
INCLUDE ABA_PARAM.INC
C
CHARACTER*8 CMNAME
C
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), TABLE(2,NPROPS/2-1)
C
C LOCAL ARRAYS
C ----------------------------------------------------------------
C EELAS - ELASTIC STRAINS
C EPLAS - PLASTIC STRAINS
C FLOW - DIRECTION OF PLASTIC FLOW
C ----------------------------------------------------------------
C
DIMENSION EELAS(6),EPLAS(6),FLOW(6), HARD(3)
C
PARAMETER(ZERO=0.D0, ONE=1.D0, TWO=2.D0, THREE=3.D0, SIX=6.D0,
1 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
C ELASTIC PROPERTIES
C
C ASSIGN VALUES TO TABLE
NVALUE=NPROPS/2-1
J=3
DO K1=1, 2
DO K2=1, NVALUE
TABLE(K1, K2)=PROPS(J)
J=J+1
END DO
END DO
EMOD=PROPS(1)
ENU=MIN(PROPS(2), ENUMAX)
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, 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
WRITE(7,*) statevariables2
WRITE(7,*) STATEV
C RECOVER ELASTIC AND PLASTIC STRAINS AND ROTATE FORWARD
C ALSO RECOVER EQUIVALENT PLASTIC STRAIN
C
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 CALCULATE PREDICTOR STRESS AND ELASTIC STRAIN
C
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
WRITE (7,*) stress
WRITE(7,*) STRESS
C
C CALCULATE EQUIVALENT VON MISES STRESS
C
SMISES=(STRESS(1)-STRESS(2))**2+(STRESS(2)-STRESS(3))**2
1 +(STRESS(3)-STRESS(1))**2
DO K1=NDI+1,NTENS
SMISES=SMISES+SIX*STRESS(K1)**2
END DO
SMISES=SQRT(SMISES/TWO)
WRITE(7,*) smises
WRITE(7,*) SMISES
C
C GET YIELD STRESS FROM THE SPECIFIED HARDENING CURVE
C

CALL UHARD(SYIEL0, HARD, EQPLAS, EQPLASRT,TIME,DTIME,TEMP,
1 DTEMP,NOEL,NPT,LAYER,KSPT,KSTEP,KINC,CMNAME,NSTATV,
2 STATEV,NUMFIELDV,PREDEF,DPRED,NVALUE,PROPS(3))
C
C DETERMINE IF ACTIVELY YIELDING
C

WRITE(7,*) syiel0
WRITE(7,*) SYIEL0
IF (SMISES.GT.(ONE+TOLER)*SYIEL0) THEN
C
C ACTIVELY YIELDING
C SEPARATE THE HYDROSTATIC FROM THE DEVIATORIC STRESS
C CALCULATE THE FLOW DIRECTION
C
SHYDRO=(STRESS(1)+STRESS(2)+STRESS(3))/THREE
WRITE (7,*) shydro
WRITE (7,*) SHYDRO
DO K1=1,NDI
FLOW(K1)=(STRESS(K1)-SHYDRO)/SMISES
END DO
DO K1=NDI+1, NTENS
FLOW(K1)=STRESS(K1)/SMISES
END DO
WRITE(7,*) flow
WRITE(7,*) FLOW
C
C SOLVE FOR EQUIVALENT VON MISES STRESS
C AND EQUIVALENT PLASTIC STRAIN INCREMENT USING NEWTON ITERATION
C
SYIELD=SYIEL0
DEQPL=ZERO
DO KEWTON=1, NEWTON
RHS=SMISES-EG3*DEQPL-SYIELD
WRITE(7,*) rhs
WRITE (7,*) RHS
WRITE(7,*) deqpl
WRITE (7,*) DEQPL
DEQPL=DEQPL+RHS/(EG3+HARD(1))
CALL UHARD(SYIEL0, HARD, EQPLAS, 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
C WRITE WARNING MESSAGE TO THE .MSG FILE
C
WRITE(7,2) NEWTON
2 FORMAT ( //,30X,***WARNING - PLASTICITY ALGORITHM DID NOT ,
1 CONVERGE AFTER ,I3, ITERATIONS )
10 CONTINUE
C
C UPDATE STRESS, ELASTIC AND PLASTIC STRAINS AND
C EQUIVALENT PLASTIC STRAIN
C
EQPLAS=EQPLAS+DEQPL
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
WRITE(7,*) eplas
WRITE(7,*) EPLAS
WRITE(7,*) eelas
WRITE(7,*) EELAS
WRITE(7,*) stress
WRITE(7,*) STRESS

WRITE(7,*) eqplas
WRITE(7,*) EQPLAS
C
C CALCULATE PLASTIC DISSIPATION
C
SPD=DEQPL*(SYIEL0+SYIELD)/TWO
C
C FORMULATE THE JACOBIAN (MATERIAL TANGENT)
C FIRST CALCULATE EFFECTIVE MODULI
C
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
C
C STORE ELASTIC AND (EQUIVALENT) PLASTIC STRAINS
C 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
STATEV(14)=15+STATEV(14)
WRITE(7,*) statevariables
WRITE(7,*) STATEV
C
RETURN
END
C
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(*),
1 PREDEF(NUMFIELDV),DPRED(*)
C
DIMENSION TABLE(2, NVALUE)
C
PARAMETER(ZERO=0.D0)
C
C SET YIELD STRESS TO LAST VALUE OF TABLE, HARDENING TO ZERO
C
SYIELD=TABLE(1, NVALUE)
HARD(1)=ZERO
C IF MORE THAN ONE ENTRY, SEARCH TABLE
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(7, 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(1)=DSYIEL/DEQPL
SYIELD=SYIEL0+(EQPLAS-EQPL0)*HARD(1)
WRITE(7,*) first
WRITE(7,*) SYIELD
GOTO 10
ENDIF
END DO
10 CONTINUE
ENDIF
RETURN
END

Thanks

Topic Tags
1 Reply
Jorgen
Posts: 3924
Moderator
(@jorgen)
Member
Joined: 2 years ago

I dont have time to review your code .
My guess is that the problem is related to the Jacobian.

-Jorgen

Topic Tags
1 Reply
Share: