*********************************************************************************** ** UMAT, FOR ABAQUS/STANDARD INCORPORATING ELASTIC-PLASTIC LINEAR ** ** ISOTROPIC HARDENING. LARGE DEFORMATION FORMULATION FOR PLANE STRAIN ** ** AND AXI-SYMMETRIC ELEMENTS. IMPLICIT INTEGRATION WITH CONSISTENT JACOBIAN ** *********************************************************************************** *********************************************************************************** ** ** ** 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