cohesive model of Needleman was input into abaqus, but wrong???

I have written vumatcodesof a cohesive model of Needleman (whose paper entilted "Numerical simulations of fast crack growth in brittle solids"), and I have implemented it into ABAQUS/Explicit, but the predicted result is wrong. I don't know why. The corresponding fortran file is as follows:

subroutine vumat(

C Read only (unmodifiable)variables -

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 (modifiable) variables -

7 stressNew, stateNew, enerInternNew, enerInelasNew )

C

include 'vaba_param.inc'

C

dimension props(nprops), density(nblock), coordMp(nblock,*),

1 charLength(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),

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

9 fieldNew(nblock,nfieldv),

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

2 enerInternNew(nblock), enerInelasNew(nblock)

C

character*80 cmname

C

C six parameters

Parameter(EX=2.71828)

snmax=props(1) ! the maxinum normal stress

deltn=props(2)! the maxinum normal displacement

stmax=props(3)! the maxinum shear stress

delts=props(4)! the maxinumsheardisplacement

t=props(5) !depth of cohesive elements

r=props(6)

C

C)

C the maxinum energy

egyn=EX*snmax*deltn

egyt=SQRT(EX/2)*stmax*delts

q=egyt/egyn

C

C)

C

do 100 i=1, nblock

C

C))))calculatestrain

sta1=stateOld(i,1)+strainInc(i,1)

sta2=stateOld(i,2)+strainInc(i,2)

sta3=stateOld(i,3)+strainInc(i,3)

sta4=stateOld(i,4)+strainInc(i,4)

Ccalculatedisplacement

detn=sta1*t

dets=sta2*t

tn=detn/deltn

ts=dets/delts

C)calculate stress and energy

sig1=(egyn/deltn)*(EX**(-tn))*(tn*(EX**(-ts*ts))+((1-q)/(r-1))

1 *(1-(EX**(-ts*ts)))*(r-tn))

sig2=(egyn/deltn)*2*(deltn/delts)*ts*(q+((r-q)/(r-1))*tn)

2 *(EX**(-tn))*(EX**(-ts*ts))

sig3=0

sig4=0

ENGY=egyn+egyn*(EX**(-tn))*((1-r+tn)*((1-q)/(r-1))

1 -(q+((r-q)/(r-1))*tn)*(EX**(-ts*ts)))

C calculate the pure normal and shear energy

ENt=egyn+egyn*((1-r)*((1-q)/(r-1))

1 -q*(EX**(-ts*ts)))

ENn=egyn+egyn*(EX**(-tn))*((1-r+tn)*((1-q)/(r-1))

1 - (q + ((r q) / (r1)) * tn))

C

C................upate the stress ##################################,,

C

stressNew(i,1)=sig1

stressNew(i,2)=sig2

stressNew(i,3)=sig3

stressNew(i,4)=sig4

C) calculate the energy difference between the predicted energy and the maxinum energy

detENt=egyt-ENt

detENn=egyn-ENn

C.....,update the state variables

stateNew(i,1)=stateOld(i,1)+strainInc(i,1)

stateNew(我,2)= stateOld(我,2)+ strainInc(我,2)

stateNew(i,3)=stateOld(i,3)+strainInc(i,3)

stateNew(i,4)=stateOld(i,4)+strainInc(i,4)

Cupdate the state variables

stateNew(i,5)=2*sta1

stateNew(i,6)=dets

stateNew(i,7)=detn

stateNew(i,8)=ENGY

stateNew(i,9)=ENt

stateNew(i,10)=ENn

stateNew(i,11)=detENt

stateNew(i,12)=detENn

C..........predict the fracture........................,,

if(stateNew(i,11)>0.01) then

stateNew(i,13)=1

else

stateNew(i,13)=0

endif

C...........predict the fracture....................,,

if(stateNew(i,12)>0.01) then

stateNew(i,14)=1

else

stateNew(i,14)=0

endif

C..........predict the fracture......,,

if( (stateNew(i,13)==0) .or. (stateNew(i,14)==0) ) then

stateNew(i,15)=1

else

stateNew(i,15)=0

endif

C

100 continue

C

返回

end

the abaqus input file is at the attachment