Dear Jorgen . I met a vumat problem,can you help me?

Dear prof

I wrote a subroutine vumat about Johnson-cook constitutive model, but 693 error occurs, some people telling me Im passing variables

which abaqus couldnt read. But I can not find any of those. Here is my vumat code about JC model. Can you help me out of it?

Thanks so much!

c User subroutine VUMAT

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 coordMp(nblock,*), charLength(nblock), props(nprops),

1 density(nblock), strainInc(nblock,ndir+nshr),

2 relSpinInc(nblock,nshr), tempOld(nblock),

3 stretchOld(nblock,ndir+nshr),

4 defgradOld(nblock,ndir+nshr+nshr),

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

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

7 enerInelasOld(nblock), tempNew(nblock),

8 stretchNew(nblock,ndir+nshr),

9 defgradNew(nblock,ndir+nshr+nshr),

1 fieldNew(nblock,nfieldv),

2 stressNew(nblock,ndir+nshr), stateNew(nblock,nstatev),

3 enerInternNew(nblock), enerInelasNew(nblock)

c

character*80 cmname

c

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

1 third = 1.d0 / 3.d0, half = 0.5d0, op5 = 1.5d0 )

c

c For plane strain, axisymmetric, and 3D cases using

c the J2 Mises Plasticity with linear hardening.

c The state variable is stored as:

c STATE(*,1) = equivalent plastic strain

c STATE(*,2) = plastic strain rate

c

c

c User needs to input

c props(1) Youngs modulus

c props(2) Poissons ratio

c props(3) A

c prpps(4) B

c prpps(5) n

c prpps(6) C

c prpps(7) m

c

c

e = props(1)

xnu = props(2)

A = props(3)

B = props(4)

EN = props(5)

C = props(6)

EM = props(7)

twomu = e / ( one + xnu )

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

thremu = op5 * twomu

1

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)

1 + twomu * strainInc(k,1) + alamda * trace

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

1 + twomu * strainInc(k,2) + alamda * trace

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

1 + twomu * strainInc(k,3) + alamda * trace

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

if ( nshr .gt. 1 ) then

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

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

end if

end do

else

do k = 1, nblock

c

c calculate hard and yieldOld

c

peeqOld=stateOld(k,1)

c

if(peeqOld.eq.zero) then

yieldOld = A

else

hard = EN * B * peeqOld ** (EN-1)

yieldOld = A + B * peeqOld ** EN

end if

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

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

1 trace

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

1 trace

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

1 trace

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

if ( nshr .gt. 1 ) then

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

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

end if

c

smean = third * ( s11 + s22 + s33 )

s11 = s11 - smean

s22 = s22 - smean

s33 = s33 - smean

if ( nshr .eq. 1 ) then

vmises = sqrt( op5*

1 (s11*s11+s22*s22+s33*s33+two*s12*s12) )

else

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

1 two * s12 * s12 + two * s13 * s13 + two * s23 *

2 s23 ) )

end if

c

tvp = C * log(stateOld(k,2))

tvp1 = tvp + one

hard1 = hard * tvp1 + yieldOld * C / stateOld(k,2)

yieldOld = yieldOld * tvp1

sigdif = vmises - yieldOld

facyld = zero

if ( sigdif .gt. zero ) facyld = one

deqps = facyld * sigdif / (thremu + hard1)

c

c Update the state variable

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

stateNew(k,2) = deqps/dt

c

c

c

c

c Update the stress

yieldNew = yieldOld + hard1 * deqps

factor = yieldNew / ( yieldNew + thremu * deqps )

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

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

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

stressNew(k,4) = s12 * factor

if ( nshr .gt. 1 ) then

stressNew(k,5) = s13 * factor

stressNew(k,6) = s23 * factor

end if

c

c Update the specific internal energy -

c

if ( nshr .eq. 1 ) then

stressPower = half * (

1 ( stressOld(k,1) + stressNew(k,1) ) * strainInc(k,1) +

2 ( stressOld(k,2) + stressNew(k,2) ) * strainInc(k,2) +

3 ( stressOld(k,3) + stressNew(k,3) ) * strainInc(k,3) ) +

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

else

stressPower = half * (

1 ( stressOld(k,1) + stressNew(k,1) ) * strainInc(k,1) +

2 ( stressOld(k,2) + stressNew(k,2) ) * strainInc(k,2) +

3 ( stressOld(k,3) + stressNew(k,3) ) * strainInc(k,3) ) +

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

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

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

end if

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

c

c Update the dissipated inelastic specific energy -

plasticWorkInc = yieldNew * deqps

enerInelasNew(k) = enerInelasOld(k)

1 + plasticWorkInc / density(k)

end do

end if

return

end

c

c

These errors usually come from numerical problems in the coding of the VUMAT itself, i.e. a divide by zero. Looking through your code, if stateold(k,2) ever takes a value of 0, which seems to be possible in your first increment, youll get a divide by zero and a log of zero.