User login

Navigation

You are here

vumat for isotropic perfect plastic

Gayan Aravinda's picture

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,

Comments

aslan mohammadpour's picture

you didn't define "third" .... third should have value in this line ...

smean = third*(s11+s22+s33)

may be there are some other mistakes .....

Frank Richter's picture

Hello

 

see Abaqus Verification Manual

4.1.28 VUMAT: rotating cylinder

You need to modify the code if you simulate in 3D.

 

Good luck

 

Frank

 

------------------------------------------
Ruhr-University
Bochum
Germany

Gayan Aravinda's picture

Thank you very much all..i got correct it..

Hi everyone. I would like to know how to get started with a VUMAT subroutine for multiphases for example a subroutine for pore water and sand. It could be helpful when anyone could share VUMAT subroutines for learning such complex programmings. Could be appreciated when some one teaches with equations. thank you.

Subscribe to Comments for "vumat for isotropic perfect plastic"

More comments

Syndicate

Subscribe to Syndicate