Problem in damage model using VUMAT

[ami143]

Hi all,

I have implemented an isotropic damage model in VUMAT subroutine, and I would like to enquire about the problem in the code. The model is implemented with VUMAT subroutine and Mazars damage model was used. I always get strange results of the force-displacement curve. I am not sure about the definition of the utility subroutine VSPRINC. Can anyone look at the code and tell me what's wrong ?

Thanks in advance.

subroutine vumat(
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,
3 stressOld, stateOld, enerInternOld, enerInelasOld,
6 tempNew, stretchNew, defgradNew, fieldNew,
5 stressNew, stateNew, enerInternNew, enerInelasNew )

include 'vaba_param.inc'
C integer ndir, nshr, nstatev, nfieldv, nprops, lanneal
C real stepTime, totalTime, dt
C

dimension props(nprops), density(nblock),
1 coordMp(nblock,*),
2 charLength(*), strainInc(nblock,ndir+nshr),
3 relSpinInc(*), tempOld(*),C0(ndir+nshr,ndir+nshr),
4 stretchOld(*), defgradOld(*),eps0(nblock,ndir+nshr),
5 fieldOld(*), stressOld(nblock,ndir+nshr),eigVal(nblock,3),
6 stateOld(nblock,nstatev), enerInternOld(nblock),
7 enerInelasOld(nblock), tempNew(*), epsp(nblock,ndir+nshr),
8 stretchNew(*), defgradNew(*), fieldNew(*),
9 stressNew(nblock,ndir+nshr), stateNew(nblock,nstatev),
1 enerInternNew(nblock), enerInelasNew(nblock),
2 epsN(nblock,ndir+nshr), stressinc(nblock,ndir+nshr)
character*80 cmname
C
C dimension
C 1 epsN(4,1),stressinc(4,1)
*USER SUBROUTINE
C real(8) :: E, ENU
PARAMETER (ZERO=0.D0,ONE=1.D0,TWO=2.D0,THREE=3.D0,
+ SIX=6.D0, NINE=9.D0,S=4.18,n_svd_Required=2)
C
C
E=props(1)
ENU=props(2)
A=PROPS(3)
B=PROPS(4)
k0=PROPS(5)
C
C
A=ONE
C ----- plane stress
C0(1,1)=E/(ONE-ENU*ENU)
C0(1,2)=(E*ENU)/(ONE-ENU*ENU)
C0(1,3)=ZERO
C0(2,1)=(E*ENU)/(ONE-ENU*ENU)
C0(2,2)=E/(ONE-ENU*ENU)
C0(2,3)=ZERO
C0(3,1)=ZERO
C0(3,2)=ZERO
C0(3,3)=(E/(ONE-ENU*ENU))*(ONE-ENU)

epsN=strainInc
do 100 i = 1, nblock
if (stepTime.eq.zero) then
do j=1,3
stressNew(i,j)=stressOld(i,j)+C0(j,j)*epsN(i,j)
enddo
else
C
epsN=strainInc
epsstar1=epsstar
D1=D
C
CALL VSPRINC(nblock,epsN,eigVal,NDIR,NSHR)
do j=1,3
epsp(i,j)=eigVal(i,j)
enddo
C
epsstar=ZERO
do j=1,3
if (epsp(i,j).GT.ZERO) then
epsstar=epsstar+epsp(i,j)*epsp(i,j)
end if
enddo
C
epsstar=sqrt(epsstar)
epsstar_temp=epsstar
depsstar=epsstar-epsstar1
C print*, 'epsstar='
C write (*,*) epsstar
C
if ((epsstar-k0)>ZERO) then
f=epsstar-k0
else
f=ZERO
endif
C ---------------------- calculate damage increment and damage ----------------------
if ((epsstar-k0)>ZERO) then
D=ONE-A*exp(-B*f)
else
D=ZERO
endif
C ---------------------- limit the maximum of damage ----------------------
if (D>ONE) then
D=ONE
endif
if (D<D1) then
D=D1
endif
C Update the stress
do j=1,3
stressNew(i,j)=stressOld(i,j)+(ONE-D)*C0(j,j)*epsN(i,j)
enddo
C State variables
stateNew(i,1)=epsstar
stateNew(i,2)=D
endif
C endif
return
end