Notifications
Clear all

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


genecodegen
Posts: 6
Topic starter
(@genecodegen)
Active Member
Joined: 11 years ago

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

Topic Tags
3 Replies
kiranckst
Posts: 93
(@kiranckst)
Trusted Member
Joined: 12 years ago

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.

Topic Tags
3 Replies
genecodegen
Posts: 6
Topic starter
(@genecodegen)
Active Member
Joined: 11 years ago

[QUOTE=lumpwood,8716]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.[/QUOTE]

Dear lumpwood

Thanks for your information, Ive modified the log part to ensure the value is greater than 0, and it seems right and its running now, I really appreciate your help.

Timmy

Reply
genecodegen
Posts: 6
Topic starter
(@genecodegen)
Active Member
Joined: 11 years ago

Dear all

I tried to use this subroutine in abaqus to replace the J-C model in a 2D compression model , but the simulation stops at the very beginning of due to the excessive distortion, and the stress distribution jumps between the two connecting elements. Can anybody come up with some suggestion? Thanks so much~

Tim

Reply
Share: