Hi, I am new to UMAT. I want to implement UMAT for isotropic hardening plasticity. I came accross this UMAT and corresponding UHARD program.
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) 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=0.4999D0,NEWTON=10,TOLER=1.0D-6) C C ---------------------------------------------------------------- C UMAT FOR ISOTRPIC 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 F YIELD STRESS VS. PLASTIC STRAIN C ---------------------------------------------------------------- C C ELASTICS PROPORTIES 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 C C ELASTIC STIFFNESS C DO K1=1, NDT 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 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 pridictor 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 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) C C Get yield stress from the specified hardening curve C NVALUE=NPROPS/2-1 C C UHARD is found in (1.1.36 UHARD)- Abaqus user Subroutines C Reference Guide,but SYIELD NOT SYIEL0, and PROPS not PROPS(3) 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,NUMPROPS,PROPS(3)) C C Determine if actively yielding C IF (SMISES .GT. ONE+TOLER*SYIEL0) THEN C C C Actively yielding C Separate the Hydrostatic from Deviatoric stress C Calculate the flow direction C SHYDRO=(STRESS(1)+STRESS(2)+STRESS(3))/3 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 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 UHARD(SYIELD,HARD,EQPLAS+DEQPL,EQPLASRT,TIME,DTIME,TEMP, 1 DTEMP,NOEL,NPT,LAYER,KSPT,KSTEP,KINC,CMNAME,NSTATV, 2 STATEV,NUMFIELDV,PREDEF,DPRED,NUMPROPS,PROPS(3)) C IF (ABS(RHS) .LT. TOLER*SYIEL0) GOTO 10 C END DO C C Write warning message to the .msg file WRITE (7,2) NEWTON 2 FORMAT (//,30X,'***Warning - Plastcity algorithm did not', 1 'Converge after ',I3,'iterations') 10 Continue C C Update Stress, Elastic and Plastic Strain and C Equivalent Plastic Strain C 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 C C Calculate Plastic Dissipation C SPD=DEQPL*(SYIEL0+SYIELD)/TWO C C Formulate the Jacobian (Materail 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 C RETURN END SUBROUTINE UHARD(SYIELD,HARD,EQPLAS,EQPLASRT,TIME,DTIME,TEMP, 1 DTEMP,NOEL,NPT,LAYER,KSPT,KSTEP,KINC,CMNAME,NSTATV, 2 STATEV,NUMFIELDV,PREDEF,DPRED,NVALUE,TABLE) C INCLUDE 'ABA_PARAM.INC' C 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 SYIELD=SYIEL0+(EQPLAS-EQPL0)*HARD(1) GOTO 10 ENDIF ENDDO 10 CONTINUE ENDIF RETURN END
Here how can I specify the yield stress and plastic strain as table, as we normally specify in Abaqus CAE, when this UMAT being used?. And how to obtain the deformation gradient (DFGRAD) as output from the abaqus?
