SUBROUTINE VUMAT(

C Pass-in arguments from Abaqus/Explicit

1 NBLOCK, NDIR, NSHR, NSTATV, NPROPS, PROPS,

2 COORDS, JCRYS, JMAC,

3 KSTEP, KINC, TIME, DTIME, TEMPOLD, DTEMP,

4 CMNAME,

5 DEFORMATIONGRADIENTOLD, DEFORMATIONGRADIENTNEW,

6 STRAININC, RELSPININC,

7 STRESSOLD, STATEOLD, ENERINTERNOLD, ENERINELASOLD,

C Pass-out arguments to Abaqus/Explicit

8 STRESSNEW, STATENEW, ENERINTERNNEW, ENERINELASNEW,

C User-defined arguments (optional, if using USDFLD, etc.)

9 FIELDOLD, FIELDNEW, LFIELD, NFFIELD)

IMPLICIT NONE

C ABA_PARAM.INC defines constants related to precision.

INCLUDE 'ABA_PARAM.INC'

C Declare pass-in arguments

INTEGER NBLOCK, NDIR, NSHR, NSTATV, NPROPS

INTEGER KSTEP, KINC

INTEGER JCRYS(NBLOCK), JMAC(NBLOCK)

INTEGER LFIELD, NFFIELD

CHARACTER*80 CMNAME

DOUBLE PRECISION PROPS(NPROPS)

DOUBLE PRECISION COORDS(NBLOCK,NDIR)

DOUBLE PRECISION TIME(2), DTIME

DOUBLE PRECISION TEMPOLD(NBLOCK), DTEMP(NBLOCK)

DOUBLE PRECISION DEFORMATIONGRADIENTOLD(NBLOCK,NDIR,NDIR)

DOUBLE PRECISION DEFORMATIONGRADIENTNEW(NBLOCK,NDIR,NDIR)

DOUBLE PRECISION STRAININC(NBLOCK,NDIR+NSHR)

DOUBLE PRECISION RELSPININC(NBLOCK,NSHR)

DOUBLE PRECISION STRESSOLD(NBLOCK,NDIR+NSHR)

DOUBLE PRECISION STATEOLD(NBLOCK,NSTATV)

DOUBLE PRECISION ENERINTERNOLD(NBLOCK), ENERINELASOLD(NBLOCK)

DOUBLE PRECISION FIELDOLD(NBLOCK,LFIELD), FIELDNEW(NBLOCK,LFIELD)

C Declare pass-out arguments

DOUBLE PRECISION STRESSNEW(NBLOCK,NDIR+NSHR)

DOUBLE PRECISION STATENEW(NBLOCK,NSTATV)

DOUBLE PRECISION ENERINTERNNEW(NBLOCK), ENERINELASNEW(NBLOCK)

C Local variables

INTEGER K, I, J

DOUBLE PRECISION EMOD, ENU, SYIELD

DOUBLE PRECISION LAMBDA, MU

DOUBLE PRECISION DEPS_PLAS(NDIR+NSHR) ! Plastic strain increment

DOUBLE PRECISION EP_EQ_OLD, EP_EQ_NEW

DOUBLE PRECISION PLAS_STRAIN_OLD(NDIR+NSHR)

DOUBLE PRECISION PLAS_STRAIN_NEW(NDIR+NSHR)

DOUBLE PRECISION STRESS_TRIAL(NDIR+NSHR)

DOUBLE PRECISION DEVIATORIC_STRESS(NDIR+NSHR)

DOUBLE PRECISION HYDRO_STRESS

DOUBLE PRECISION VON_MISES_STRESS

DOUBLE PRECISION FLOW_DIRECTION(NDIR+NSHR)

DOUBLE PRECISION D_GAMMA ! Plastic multiplier

C Material properties (example: E, nu, YieldStress)

C PROPS(1) = Young's Modulus

C PROPS(2) = Poisson's Ratio

C PROPS(3) = Initial Yield Stress

EMOD = PROPS(1)

ENU = PROPS(2)

SYIELD = PROPS(3)

LAMBDA = (EMOD * ENU) / ((1.0D0 + ENU) * (1.0D0 - 2.0D0 * ENU))

MU = EMOD / (2.0D0 * (1.0D0 + ENU))

C --- START OF LOOP OVER MATERIAL POINTS IN THE BLOCK ---

DO K = 1, NBLOCK

C Initialize energies for this point

ENERINTERNNEW(K) = ENERINTERNOLD(K)

ENERINELASNEW(K) = ENERINELASOLD(K)

C Retrieve plastic strain components from state variables

C STATEOLD(K,1) = plastic strain_11

C STATEOLD(K,2) = plastic strain_22

C STATEOLD(K,3) = plastic strain_33

C STATEOLD(K,4) = plastic shear strain_12 (engineering form gamma_12)

C STATEOLD(K,5) = plastic shear strain_23 (engineering form gamma_23)

C STATEOLD(K,6) = plastic shear strain_13 (engineering form gamma_31)

C STATEOLD(K,7) = equivalent plastic strain (scalar)

DO I = 1, NDIR+NSHR

PLAS_STRAIN_OLD(I) = STATEOLD(K,I)

END DO

EP_EQ_OLD = STATEOLD(K,NDIR+NSHR+1)

C ******************************************************************

C * ELASTIC PREDICTOR STEP *

C * Stress and strain components are in the corotational system. *

C * No explicit rotation is needed for these quantities here. *

C ******************************************************************

DO I = 1, NDIR ! Direct components

STRESS_TRIAL(I) = STRESSOLD(K,I) + LAMBDA * (

& STRAININC(K,1) + STRAININC(K,2) + STRAININC(K,3)) +

& 2.0D0 * MU * STRAININC(K,I)

END DO

DO I = NDIR+1, NDIR+NSHR ! Shear components

STRESS_TRIAL(I) = STRESSOLD(K,I) + 2.0D0 * MU * STRAININC(K,I)

END DO

C ******************************************************************

C * PLASTICITY CALCULATION (CONCEPTUAL) *

C * This is a placeholder for a real plasticity model (e.g., Von Mises)

C ******************************************************************

C Calculate Von Mises equivalent stress for the trial stress

HYDRO_STRESS = (STRESS_TRIAL(1)+STRESS_TRIAL(2)+STRESS_TRIAL(3))

& / 3.0D0

DO I = 1, NDIR

DEVIATORIC_STRESS(I) = STRESS_TRIAL(I) - HYDRO_STRESS

END DO

DO I = NDIR+1, NDIR+NSHR

DEVIATORIC_STRESS(I) = STRESS_TRIAL(I)

END DO

VON_MISES_STRESS = SQRT( (3.0D0/2.0D0) * (

& DEVIATORIC_STRESS(1)**2 + DEVIATORIC_STRESS(2)**2 +

& DEVIATORIC_STRESS(3)**2 +

& 2.0D0 * (DEVIATORIC_STRESS(4)**2 +

& DEVIATORIC_STRESS(5)**2 +

& DEVIATORIC_STRESS(6)**2) ) )

C Yield condition (very basic, no hardening for simplicity here)

IF (VON_MISES_STRESS .GT. SYIELD) THEN

C Material has yielded - Perform plastic correction (conceptual)

D_GAMMA = (VON_MISES_STRESS - SYIELD) / (3.0D0 * MU) ! Simplified

C Calculate flow direction (normalized deviatoric stress)

DO I = 1, NDIR+NSHR

FLOW_DIRECTION(I) = (3.0D0/2.0D0) * & DEVIATORIC_STRESS(I) / VON_MISES_STRESS

END DO

C Calculate plastic strain increment components

DO I = 1, NDIR+NSHR

DEPS_PLAS(I) = D_GAMMA * FLOW_DIRECTION(I)

END DO

C Update stress tensor (return to yield surface)

DO I = 1, NDIR

STRESSNEW(K,I) = STRESS_TRIAL(I) - 2.0D0 * MU * D_GAMMA * & FLOW_DIRECTION(I)

END DO

DO I = NDIR+1, NDIR+NSHR

STRESSNEW(K,I) = STRESS_TRIAL(I) - 2.0D0 * MU * D_GAMMA *

& FLOW_DIRECTION(I)

END DO

C Update plastic strains

C **************************************************************

C * IMPORTANT: Updating Plastic Strain State Variables *

C * These plastic strain components are in the corotational *

C * system. The DEPS_PLAS components were calculated in this *

C * same system. Therefore, the update is a direct addition. *

C * NO EXPLICIT ROTATION of PLAS_STRAIN_OLD or DEPS_PLAS *

C * (like using ROTSIG in UMAT) is needed here because the *

C * entire calculation framework, including the state *

C * variables, is already in the coordinate system that *

C * rotates with the material. *

C **************************************************************

DO I = 1, NDIR+NSHR

PLAS_STRAIN_NEW(I) = PLAS_STRAIN_OLD(I) + DEPS_PLAS(I)

END DO

C Update equivalent plastic strain (conceptual)

EP_EQ_NEW = EP_EQ_OLD + D_GAMMA ! Simplified, assumes D_GAMMA is related to dep_eq

C Update plastic work (simplified)

DO I = 1, NDIR+NSHR

ENERINELASNEW(K) = ENERINELASNEW(K) +

& STRESSNEW(K,I)*DEPS_PLAS(I)

IF(I.LE.NDIR) THEN ! direct stress/strain

ENERINELASNEW(K) = ENERINELASNEW(K) +

& STRESSNEW(K,I)*DEPS_PLAS(I)

ELSE ! shear stress/strain (assuming DEPS_PLAS is tensor shear)

ENERINELASNEW(K) = ENERINELASNEW(K) +

& STRESSNEW(K,I)*DEPS_PLAS(I)

ENDIF

ENDDO

ELSE

C Material behaves elastically

DO I = 1, NDIR+NSHR

STRESSNEW(K,I) = STRESS_TRIAL(I)

PLAS_STRAIN_NEW(I) = PLAS_STRAIN_OLD(I) ! No change

END DO

EP_EQ_NEW = EP_EQ_OLD ! No change

END IF

C Store updated state variables

DO I = 1, NDIR+NSHR

STATENEW(K,I) = PLAS_STRAIN_NEW(I)

END DO

STATENEW(K,NDIR+NSHR+1) = EP_EQ_NEW

C Update elastic strain energy (simplified, total strain energy - plastic work)

C A more rigorous update would use elastic strain components.

ENERINTERNNEW(K) = 0.0D0

DO I = 1, NDIR+NSHR

ENERINTERNNEW(K) = ENERINTERNNEW(K) +

& 0.5D0 * (STRESSNEW(K,I) + STRESSOLD(K,I)) *

& STRAININC(K,I)

ENDDO

ENERINTERNNEW(K) = ENERINTERNOLD(K) + ENERINTERNNEW(K) -

& (ENERINELASNEW(K) - ENERINELASOLD(K))

END DO ! K - Loop over material points

C --- END OF LOOP OVER MATERIAL POINTS ---

RETURN

END

Makenzie Diaz I hope this is helpful to you.

have you checked the good abaqus documentation? state variables are calculated at material points, hence they rotate with material coordinate system (which abaqus takes care of)

For a VUMAT, you must do it manually by setting up the rotation matrix and conducting the rotation as a two-fold matrix multiplication. I typically do this by writing a subroutine dedicated for rotation. I derived an elementary version from my one and enclosed it to this answer. Most likely you need to adapt it to fit the input/output structure you want.