I have written vumat codes of 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 maxinum shear displacement
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)))) calculate strain
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)
C calculate displacement
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)/(r-1))*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(i,2)=stateOld(i,2)+strainInc(i,2)
stateNew(i,3)=stateOld(i,3)+strainInc(i,3)
stateNew(i,4)=stateOld(i,4)+strainInc(i,4)
C update 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
return
end
the abaqus input file is at the attachment
| Attachment | Size |
|---|---|
| abaqus input file | 17.72 KB |