Hi all
I am a new ABAQUS VUMAT user. Recently, I try to substitute the explicit intergation method in VUMAT by back EULER intergation method. But the calculated equivalent stress is nonuniform for simple compression test. Would anyone pls help me find the reasons? Many thanks. Attached pls find the VUMAT subroutine and INP file.
J.G. Liu
Ph.D
HK Polyu
C C Coding for Isotropic
C
C Coding for Isotropic Hardening Plasticity VUMAT
C
SUBROUTINE VUMAT(
C Read only -
1 NBLOCK, NDIR, NSHR, NSTATEV, NFIELDV, NPROPS, LANNEAL,
2 STEPTIME, TOTALTIME, DT, CMNAME, COORDMP, CHARLENGTH,
3 PROPS, DENSITY, STRAININC, RELSPININC,
4 TEMPOLD, STRETCHOLD, DEFGRADOLD, FIELDOLD,
5 STRESSOLD, STATEOLD, ENERINTERNOLD, ENERINELASOLD,
6 TEMPNEW, STRETCHNEW, DEFGRADNEW, FIELDNEW,
C Write only -
7 STRESSNEW, STATENEW, ENERINTERNNEW, ENERINELASNEW)
C
INCLUDE 'VABA_PARAM.INC'
C
DIMENSION PROPS(NPROPS), DENSITY(NBLOCK), COORDMP(NBLOCK),
1 CHARLENGTH(NBLOCK), STRAININC(NBLOCK, NDIR+NSHR),
2 RELSPININC(NBLOCK, NSHR), TEMPOLD(NBLOCK),
3 STRETCHOLD(NBLOCK, NDIR+NSHR),DEFGRADOLD(NBLOCK,NDIR+NSHR+NSHR),
4 FIELDOLD(NBLOCK, NFIELDV), STRESSOLD(NBLOCK, NDIR+NSHR),
5 STATEOLD(NBLOCK, NSTATEV), ENERINTERNOLD(NBLOCK),
6 ENERINELASOLD(NBLOCK), TEMPNEW(NBLOCK),
7 STRETCHNEW(NBLOCK, NDIR+NSHR),DEFGRADNEW(NBLOCK,NDIR+NSHR+NSHR),
8 FIELDNEW(NBLOCK, NFIELDV), STRESSNEW(NBLOCK,NDIR+NSHR),
9 STATENEW(NBLOCK, NSTATEV), ENERINTERNNEW(NBLOCK),
1 ENERINELASNEW(NBLOCK)
CHARACTER*80 CMNAME
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,HALF=0.5D0)
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
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
NVALUE=NPROPS/2-1
C
C CALCULATE PREDICTOR STRESS AND ELASTIC STRAIN
C
IF ( STEPTIME .EQ. ZERO ) THEN
DO K = 1, NBLOCK
TRACE = STRAININC(K,1) + STRAININC(K,2) + STRAININC(K,3)
STRESSNEW(K,1) = STRESSOLD(K,1)
* + EG2 * STRAININC(K,1) + ELAM * TRACE
STRESSNEW(K,2) = STRESSOLD(K,2)
* + EG2 * STRAININC(K,2) + ELAM * TRACE
STRESSNEW(K,3) = STRESSOLD(K,3)
* + EG2 * STRAININC(K,3) + ELAM * TRACE
STRESSNEW(K,4)=STRESSOLD(K,4) + EG2 * STRAININC(K,4)
IF ( NSHR .GT. 1 ) THEN
STRESSNEW(K,5)=STRESSOLD(K,5) + EG2 * STRAININC(K,5)
STRESSNEW(K,6)=STRESSOLD(K,6) + EG2 * STRAININC(K,6)
END IF
END DO
ELSE
C
C111111111111
C
DO K = 1, NBLOCK
C
TRACE = STRAININC(K,1) + STRAININC(K,2) + STRAININC(K,3)
STRESSNEW(K,1) = STRESSOLD(K,1)
* + EG2 * STRAININC(K,1) + ELAM * TRACE
STRESSNEW(K,2) = STRESSOLD(K,2)
* + EG2 * STRAININC(K,2) + ELAM * TRACE
STRESSNEW(K,3) = STRESSOLD(K,3)
* + EG2 * STRAININC(K,3) + ELAM * TRACE
STRESSNEW(K,4)=STRESSOLD(K,4) + EG2 * STRAININC(K,4)
IF ( NSHR .GT. 1 ) THEN
STRESSNEW(K,5)=STRESSOLD(K,5) + EG2 * STRAININC(K,5)
STRESSNEW(K,6)=STRESSOLD(K,6) + EG2 * STRAININC(K,6)
END IF
C
SMISES=(STRESSNEW(K,1)-STRESSNEW(K,2))**2 +
1 (STRESSNEW(K,2)-STRESSNEW(K,3))**2 +
2 (STRESSNEW(K,3)-STRESSNEW(K,1))**2
DO 90 K1=NDIR+1,NDIR+NSHR
SMISES=SMISES+SIX*STRESSNEW(K,K1)**2
90 CONTINUE
SMISES=SQRT(SMISES/TWO)
C
C RECOVER ELASTIC AND PLASTIC STRAINS AND ROTATE FORWARD
C ALSO RECOVER EQUIVALENT PLASTIC STRAIN
C
C CALL ROTSIG(STATEOLD(K,1), DROT, EELAS, 2, NDIR, NSHR)
C CALL ROTSIG(STATEOLD(K,1+NDIR+NSHR), DROT, EPLAS, 2, NDIR, NSHR)
DO K1=1, NDIR+NSHR
EELAS(K1)=STATEOLD(K,K1)
EPLAS(K1)=STATEOLD(K,K1+NDIR+NSHR)
END DO
EQPLAS=STATEOLD(K,1+2*(NDIR+NSHR))!!!
C
C GET YIELD STRESS FROM THE SPECIFIED HARDENING CURVE
C
CALL VUHARD(SYIEL0, HARD, EQPLAS, EQPLASRT,TOTALTIME,DT,TEMP,
1 DTEMP,NOEL,NPT,LAYER,KSPT,KSTEP,KINC,CMNAME,
2 NUMFIELDV,NVALUE,PROPS(3))
C
C DETERMINE IF ACTIVELY YIELDING
C
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=(STRESSNEW(K,1)+STRESSNEW(K,2)+STRESSNEW(K,3))/THREE
DO K1=1,NDIR
FLOW(K1)=(STRESSNEW(K,K1)-SHYDRO)/SMISES
END DO
DO K1=NDIR+1,NDIR+NSHR
FLOW(K1)=STRESSNEW(K,K1)/SMISES
END DO
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
DEQPL=DEQPL+RHS/(EG3+HARD(1))
CALL VUHARD(SYIELD,HARD,EQPLAS+DEQPL,EQPLASRT,TOTALTIME,DT,TEMP,
1 DTEMP,NOEL,NPT,LAYER,KSPT,KSTEP,KINC,CMNAME,NUMFIELDV,
2 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
DO K1=1,NDIR
STRESSNEW(K,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=NDIR+1,NDIR+NSHR
STRESSNEW(K,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
C
C CALCULATE PLASTIC DISSIPATION
C
SPD=DEQPL*(SYIEL0+SYIELD)/TWO
ENDIF
C
C STORE ELASTIC AND (EQUIVALENT) PLASTIC STRAINS
C IN STATE VARIABLE ARRAY
C
DO K1=1, NDIR+NSHR
STATENEW(K,K1)=EELAS(K1)
STATENEW(K,K1+NDIR+NSHR)=EPLAS(K1)
END DO
STATENEW(K,1+2*(NDIR+NSHR)) = STATEOLD(K,1+2*(NDIR+NSHR)) + DEQPL
C
C UPDATE THE SPECIFIC INTERNAL ENERGY -
C
IF ( NSHR .EQ. 1 ) THEN
STRESSPOWER = HALF * (
* ( STRESSOLD(K,1) + STRESSNEW(K,1) ) * STRAININC(K,1) +
* ( STRESSOLD(K,2) + STRESSNEW(K,2) ) * STRAININC(K,2) +
* ( STRESSOLD(K,3) + STRESSNEW(K,3) ) * STRAININC(K,3) ) +
* ( STRESSOLD(K,4) + STRESSNEW(K,4) ) * STRAININC(K,4)
ELSE
STRESSPOWER = HALF * (
* ( STRESSOLD(K,1) + STRESSNEW(K,1) ) * STRAININC(K,1) +
* ( STRESSOLD(K,2) + STRESSNEW(K,2) ) * STRAININC(K,2) +
* ( STRESSOLD(K,3) + STRESSNEW(K,3) ) * STRAININC(K,3) ) +
* ( STRESSOLD(K,4) + STRESSNEW(K,4) ) * STRAININC(K,4) +
* ( STRESSOLD(K,5) + STRESSNEW(K,5) ) * STRAININC(K,5) +
* ( STRESSOLD(K,6) + STRESSNEW(K,6) ) * STRAININC(K,6)
END IF
ENERINTERNNEW(K) = ENERINTERNOLD(K) + STRESSPOWER / DENSITY(K)
C
C UPDATE THE DISSIPATED INELASTIC SPECIFIC ENERGY -
C
PLASTICWORKINC = HALF * ( SYIEL0 + SYIELD ) * DEQPL
ENERINELASNEW(K) = ENERINELASOLD(K) + PLASTICWORKINC / DENSITY(K)
ENDDO
END IF
C
RETURN
END
SUBROUTINE VUHARD(SYIELD,HARD,EQPLAS,EQPLASRT,TIME,DTIME,TEMP,
1 DTEMP,NOEL,NPT,LAYER,KSPT,KSTEP,KINC,CMNAME,
2 NUMFIELDV, NVALUE,TABLE)
INCLUDE 'VABA_PARAM.INC'
CHARACTER*80 CMNAME
DIMENSION HARD(3)
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
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 EXIT
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)
GOTO 10
ENDIF
END DO
10 CONTINUE
ENDIF
RETURN
END
I don't have time to check
I don't have time to check your code closely. One way to debug your subroutine is to simulate one element with your vumat and output all stresses, elastic strains, current hardening coeff., equivalent plastic strain, and current strength. Then you can check whether abaqus updates those internal variable correctly. Meanwhile, simulate one element using Abaqus's isotropic plasticity model (assign the same plastic behavior as you do in vumat, check keyword '*plastic' in abaqus) and see whether you get the same stress from your vumat.
thanks
Dear Dr.wei
Thanks for your suggestion. I have find the bug in my VUMAT subroutine. Your suggestion will also be very important to me to debugging the more complicated VUMAT.
Thanks again.
Hello Dr.Liu I am trying
Hello Dr.Liu
I am trying to implement isotropic hardenining in VUMAT. I have a modified Jhonson cook model which I have to use as constitutive equation for my material behavior.
Do you think I can call this equation through UHARD as you have done in the above case and use it in VUMAT ? Or can you suggest me a better way of implementing this modified JC equation in VUMAT ?