vumat for isotropic perfect plastic

Hi all,

        I wrote VUMAT subroutine using fortran 77 for perfect plasticity unisng von mises plasticity..when analysed it using unidirectional tensile simulation for a simple cube model...it seems to me that it only goes through elastic deformation and not plastic...i tried to find out what is the error..but could not find it..could anyone please find what is the error in this code..i have attached the VUMAT code for the reference...thank you very much...

 

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),

1 coordMp(nblock,*), nElement(nblock), nMatPoint(nblock),

2 charLength(*), strainInc(nblock,ndir+nshr),

3 relSpinInc(*), tempOld(*),nLayer(nblock),

4 stretchOld(*), defgradOld(*),nSecPoint(nblock),

5 fieldOld(*), stressOld(nblock,ndir+nshr),

6 stateOld(nblock,nstatev), enerInternOld(nblock),

7 enerInelasOld(nblock), tempNew(*),

8 stretchNew(*), defgradNew(*),

9 fieldNew(*), stressNew(nblock,ndir+nshr),

1 stateNew(nblock,nstatev),

2 enerInternNew(nblock), enerInelasNew(nblock)

 

 

C

 

 

character*80 cmname

 

 

C

 

 

parameter(zero = 0.d0, one = 1.d0, two = 2.d0, half = .5d0,

& Three = 3.d0, op5 = 1.5d0)

 

 

C Assume purely elastic material at the beginning of the analysis

 

C

 

 

E = props(1)

xnu = props(2)

yield1 = props(3)

pls1 = props(4)

twomu = E /( one + xnu )

alambda =xnu*twomu/(one-two*xnu)

thremu = op5 * twomu

 

 

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)

& + twomu*strainInc(k,1)+alambda*trace

stressNew(k,2)=stressOld(k,2)

& + twomu*strainInc(k,2)+alambda*trace

stressNew(k,3)=stressOld(k,3)

& + twomu*strainInc(k,3)+alambda*trace

stressNew(k,4)=stressOld(k,4)+twomu*strainInc(k,4)

stressNew(k,5)=stressOld(k,5)+twomu*strainInc(k,5)

stressNew(k,6)=stressOld(k,6)+twomu*strainInc(k,6)

end do

 

 

c else

 

 

do k = 1,nblock

trace=strainInc(k,1)+strainInc(k,2)+strainInc(k,3)

 

 

C Trial Stress

 

 

s11=stressOld(k,1)+twomu*strainInc(k,1)

& + alambda*trace

s22=stressOld(k,2)+twomu*strainInc(k,2)

& + alambda*trace

s33=stressOld(k,3)+twomu*strainInc(k,3)

& + alambda*trace

s12=stressOld(k,4)+twomu*strainInc(k,4)

s13=stressOld(k,5)+twomu*strainInc(k,5)

s23=stressOld(k,6)+twomu*strainInc(k,6)

 

 

C

 

 

smean = third*(s11+s22+s33)

 

 

C

 

 

s11 = s11 - smean

s22 = s22 - smean

s33 = s33 - smean

 

 

C

 

 

vmises = sqrt(op5*((s11*s11)+(s22*s22)+(s33*s33)

& +(two*s12*s12)+(two*s13*s13)+(two*s23*s23)))

 

 

c print*,'value is',vmises

 

C

 

 

sigdif=vmises-yieldNew

facyld=zero

if (sigdif.eq.vmises)then

factor=one

 

 

C update the stress

 

 

else if (sigdif.lt.vmises)then

yieldNew=yield1

facyld=one

deqps = facyld * sigdif /thremu

 

 

C update the stress

 

 

factor = yieldNew / ( yieldNew + thremu * deqps)

 

 

c

 

 

stressNew(k,1) = s11 * factor + smean

stressNew(k,2) = s22 * factor + smean

stressNew(k,3) = s33 * factor + smean

stressNew(k,4) = s12 * factor

stressNew(k,5) = s13 * factor

stressNew(k,6) = s23 * factor

end if

 

 

C

 

C Update the state variables

 

C

 

 

stateNew(k,1) = stateOld(k,1) + deqps

 

 

C

 

C Update the specific internal energy

 

C

 

C

 

 

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)

 

 

C

 

 

enerInternNew(k) = enerInternOld(k) + (stressPower/density(k))

 

 

c

 

C Update the dissipated inelastic specific energy

 

C

 

 

plasticWorkInc = yield1* deqps

enerInelasNew(k) = enerInelasOld(k)

& + plasticWorkInc / density(k)

end do

 

 

c end if

 

 

return

end

 

 

C

Regards,