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,
3 stressOld, stateOld, enerInternOld, enerInelasOld,
6 tempNew, stretchNew, defgradNew, fieldNew,
C Write only -
5 stressNew, stateNew, enerInternNew, enerInelasNew )
C
include 'vaba_param.inc'
character*80 cmname
C
dimension props(nprops), density(nblock),
1coordMp(nblock,*),
2charLength(*), strainInc(nblock,ndir+nshr),
3relSpinInc(*), tempOld(*),
4stretchOld(*), defgradOld(*),
5fieldOld(*), stressOld(nblock,ndir+nshr),
6stateOld(nblock,nstatev), enerInternOld(nblock),
7enerInelasOld(nblock), tempNew(*),
8stretchNew(*), defgradNew(*), fieldNew(*),
9stressNew(nblock,ndir+nshr), stateNew(nblock,nstatev),
1enerInternNew(nblock), enerInelasNew(nblock),
2DDSDDE(ndir+nshr,ndir+nshr),stressLAO(nblock,ndir+nshr),
3stressXIN(nblock,ndir+nshr),EELAS(nblock,ndir+nshr),
4EPLAS(nblock,ndir+nshr), ALPHA(nblock,ndir+nshr),
4FLOW(nblock,ndir+nshr),OLDS(nblock,ndir+nshr),
7Yacobi1(5,5),Yacobi2(5,5), OLDPL(nblock,ndir+nshr),
8DSTRESSNEW(nblock,ndir+nshr)
C EPLAS - PLASTIC STRAINS
C ALPHA - SHIFT TENSOR
C FLOW - PLASTIC FLOW DIRECTIONS
C OLDS - STRESS AT START OF INCREMENT
C OLDPL - PLASTIC STRAIN AT START OF INCREMENT
C EQPLAS - EQUVILENT PLASTIC
C
C ----------------------------------------------------------------
C VUMAT FOR ISOTROPIC ELASTICITY AND MISES PLASTICITY
C WITH ISOTROPIC HARDENING - CANNOT BE USED FOR PLANE STRESS
C ----------------------------------------------------------------
C PROPS(1) - E
C PROPS(2) - NU
C PROPS(3..) - SYIELD AN HARDENING DATA
C ----------------------------------------------------------------
C
C
C RECOVER ELASTIC STRAIN, PLASTIC STRAIN AND SHIFT TENSOR AND ROTATE
C NOTE: USE CODE 1 FOR (TENSOR) STRESS, CODE 2 FOR (ENGINEERING) STRAIN
C
C DIMENSION EELAS(nblock,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
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
aa=PROPS(3)
rr=PROPS(4)
panduan=1
C
C ELASTIC STIFFNESS
C
DO i=1,nblock
C
DO k=1,ndir+nshr
EELAS(i,k)=stateOLD(i,k)
EPLAS(i,k)=stateOLD(i,k+6)
ENDDO
EQPLAS=stateOLD(i,13)
VGI=stateOLD(i,15)
VGIindex=stateOLD(i,16)
smcs=stateOLD(i,17)
C
IF(totaltime.EQ.dt) THEN
VGI=0.0
VGIindex=0.0
ENDIF
C
DO K1=1, ndir+nshr
DO K2=1, ndir+nshr
DDSDDE(K1, K2)=0.0
END DO
DstressNew(i, K1)=0.0
END DO
CCC define elastic matix
DO K1=1, ndir
DO K2=1, ndir
DDSDDE(K2, K1)=ELAM
END DO
DDSDDE(K1, K1)=EG2+ELAM
END DO
DO K1=ndir+1, ndir+nshr
DDSDDE(K1, K1)=EG2
END DO
CCC Calculate strss increment
DO K1=1, ndir+nshr
DO K2=1, ndir+nshr
DSTRESSNEW(i,K2)=DSTRESSNEW(i,K2)+DDSDDE(K2, K1)*STRAININC(i,K1)
END DO
END DO
CCC Calculate trial strss
DO K2=1, ndir+nshr
STRESSNEW(i,K2)=STRESSOLD(i,K2)+DSTRESSNEW(i,K2)
END DO
CCC Update trail elastic strain tensor
DO k=1,ndir
EELAS(i,k)=EELAS(i,k)+STRAININC(i,k)
ENDDO
DO k=ndir+1,ndir+nshr
EELAS(i,k)=EELAS(i,k)+STRAININC(i,k)*2
ENDDO
C
SMISES=(stressNew(i,1)-stressNew(i,2))**2
1 +(stressNew(i,2)-stressNew(i,3))**2
2 +(stressNew(i,3)-stressNew(i,1))**2
DO K1=ndir+1,ndir+nshr
SMISES=SMISES+SIX*stressNew(i,K1)**2
END DO
SMISES=SQRT(SMISES/TWO)
CCCCCCCCCCC
NVALUE=NPROPS/2-2
CALL UHARD(SYIEL0, HARD, EQPLAS,NVALUE,PROPS(5))
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(i,1)+STRESSNEW(i,2)+STRESSNEW(i,3))/THREE
DO K1=1,ndir
FLOW(i,K1)=(STRESSNEW(i,K1)-SHYDRO)/SMISES
END DO
DO K1=ndir+1, ndir+nshr
FLOW(i,K1)=STRESSNEW(i,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)
CALL UHARD(SYIELD,HARD,EQPLAS+DEQPL,NVALUE,PROPS(5))
IF(ABS(RHS).LT.TOLER*SYIEL0) GOTO 10
END DO
C
C WRITE WARNING MESSAGE TO THE .MSG FILE
C
C WRITE(7,2) NEWTON
C2 FORMAT(//,30X,'***WARNING - PLASTICITY ALGORITHM DID NOT
',
C 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(i,K1)=FLOW(i,K1)*SYIELD+SHYDRO
EPLAS(i,K1)=EPLAS(i,K1)+THREE/TWO*FLOW(i,K1)*DEQPL
EELAS(i,K1)=EELAS(i,K1)-THREE/TWO*FLOW(i,K1)*DEQPL
END DO
DO K1=NDIR+1,ndir+nshr
STRESSNEW(i,K1)=FLOW(i,K1)*SYIELD
EPLAS(i,K1)=EPLAS(i,K1)+THREE*FLOW(i,K1)*DEQPL
EELAS(i,K1)=EELAS(i,K1)-THREE*FLOW(i,K1)*DEQPL
END DO
EQPLAS=EQPLAS+DEQPL
C
ENDIF
C
SMM=(stressNew(i,1)+stressNew(i,2)+stressNew(i,3))/3.0
SMISES=(stressNew(i,1)-stressNew(i,2))**2
1 +(stressNew(i,2)-stressNew(i,3))**2
2 +(stressNew(i,3)-stressNew(i,1))**2
DO K1=ndir+1,ndir+nshr
SMISES=SMISES+SIX*stressNew(i,K1)**2
END DO
SMISES=SQRT(SMISES/TWO)
C CALCULATE PLASTIC DISSIPATION
CCC calculate strain energy
STRESSPOWER = 0.5 * (
*( STRESSOLD(i,1) + STRESSNEW(i,1) ) * STRAININC(i,1) +
*( STRESSOLD(i,2) + STRESSNEW(i,2) ) * STRAININC(i,2) +
*( STRESSOLD(i,3) + STRESSNEW(i,3) ) * STRAININC(i,3) )
DO K1=ndir+1,ndir+nshr
STRESSPOWER=STRESSPOWER+(STRESSOLD(i,k1)+STRESSNEW(i,k1))
1 *STRAININC(i,k1)
END DO
C
ENERINTERNNEW(i) = ENERINTERNOLD(i) + STRESSPOWER / DENSITY(i)
CCC Calculate dissipation energy
PLASTICWORKINC =SMISES*EQPLAS
ENERINELASNEW(i) = ENERINELASOLD(i) + PLASTICWORKINC /
DENSITY(i)
C
DO k=1,ndir+nshr
stateNew(i,k)=EELAS(i,k)
stateNew(i,k+6)=EPLAS(i,k)
END DO
stateNew(i,13)=EQPLAS
C
IF(SMISES.GT.0) THEN
RRSS=SMM/SMISES
ENDIF
SMCS=EQPLAS-aa*2.718**(-1.5*RRSS)
VGI=VGI+2.718**(1.5*RRSS)*DEQPL
VGIindex=VGI-rr
stateNew(i,14)=RRSS
stateNew(i,15)=VGI
stateNew(i,16)=VGIindex
stateNew(i,17)=SMCS
C
IF(panduan.EQ.1) THEN
C
IF(VGIindex.GT.0) THEN
stateNew(i,20)=0
c GOTO 200
ELSE
stateNew(i,20)=1
ENDIF
C
ELSEIF(panduan.EQ.2) THEN
IF(SMCS.GT.0) THEN
stateNew(i,20)=0
c GOTO 200
ELSE
stateNew(i,20)=1
ENDIF
c
ELSE
stateNew(i,20)=0
ENDIF
C
ENDDO
C
200 RETURN
END
C
SUBROUTINE UHARD(SYIELD,HARD,EQPLAS,NVALUE,TABLE)
c include 'vaba_param.inc'
c character*80 cmname
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=0.0
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
DEQPL=EQPL1-EQPL0
SYIEL0=TABLE(1, K1)
SYIEL1=TABLE(1, K1+1)
DSYIEL=SYIEL1-SYIEL0
HARD=DSYIEL/DEQPL
SYIELD=SYIEL0+(EQPLAS-EQPL0)*HARD
GOTO 10
ENDIF
END DO
10 CONTINUE
ENDIF
C
RETURN
END