## You are here

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

- Gayan Aravinda's blog
- Log in or register to post comments
- 5183 reads

## Comments

## you didn't define "third"

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

smean = third*(s11+s22+s33)

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

## code is in manual

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

## hello

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

## How to write a VUMAT for multiphases

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.