sirig
2009-08-18, 21:03
Hello ,
I have written a VUMAT for a viscoelastic model (I obtained it from one Mr.Frank ritcher's thesis) which best represents fabric behavior in a ballistic impact simulation.I have attached the constitutive relation in the jpg file.
VUMAT
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(*), coordMp(nblock,*),
1 charLength(*), strainInc(nblock,ndir+nshr),
2 relSpinInc(*), tempOld(*),
3 stretchOld(*),
4 defgradOld(nblock,ndir+nshr+nshr),
5 fieldOld(*), stressOld(nblock,ndir+nshr),
6 stateOld(nblock,1), enerInternOld(*),
7 enerInelasOld(*), tempNew(*),
8 stretchNew(*),
8 defgradNew(nblock,ndir+nshr+nshr),
9 fieldNew(*),
1 stressNew(nblock,ndir+nshr), stateNew(nblock,1),
2 enerInternNew(*), enerInelasNew(*),
C
character*80 cmname
c
real,kk1,gk1,kk2,gk2,gv,kv
real,A,B,D,F,G,H
real,phi1,phi2,zeta1,zeta2
real,psi1,psi2,omega1,omega2
real,mu1,mu2,lamda1,lamda2
real,STRtrace,SIGtrace,DSTRtrace,DSIGtrace
real,Term1,Term2,Term3,Term4,Term5
real,SIG1
c
c Model coefficients
c
kk1=9.33E+08
gk1=1.40E+09
kk2=2.27E+10
gk2=3.40E+10
gv=1.00E+06
c kv (infinite-volume conserved flow in dashpot)
f1t=2610
c
omega1=0.0
omega2=1.0+kk2/kk1
psi1=1/(4.0*gv)
psi2=1/(4.0*gk1)
mu1=gk2/(2.0*gv)
mu2=(1.0+gk2/gk1)/2.0
lamda1=(-gk2/gv)/3.0
lamda2=(kk2/kk1-gk2/gk1)/3.0
phi1=-1.0/(6.0*gv)
phi2=1.0/(9.0*kk1)-1.0/(6.0*gk1)
zeta1=0.0
zeta2=1.0/(3.0*kk1)
c
c
A=dt*lamda1+lamda2
B=2.0*mu1*dt+2.0*mu2
D=phi1*dt+phi2
F=2.0*psi1*dt+2.0*psi2
G=zeta1*dt+zeta2
H=omega1*dt+omega2
c
Real,DIMENSION(nblock,ndir+nshr)::strainCurrent
c
save
c
do 100 i= 1,nblock
c
c Trace of current strain tensor from deformation gradient tensor
strainCurrent(i,1)=0.5*((defgradOld(i,1)*defgradOl d(i,1))+(defgradOld(i,7)*defgradOld(i,7))+(defgrad Old(i,6)*defgradOld(i,6))-1)
strainCurrent(i,2)=0.5*((defgradOld(i,4)*defgradOl d(i,4))+(defgradOld(i,2)*defgradOld(i,2))+(defgrad Old(i,8)*defgradOld(i,8))-1)
strainCurrent(i,3)=0.5*((defgradOld(i,9)*defgradOl d(i,9))+(defgradOld(i,5)*defgradOld(i,5))+(defgrad Old(i,3)*defgradOld(i,3))-1)
STRtrace=strainCurrent(i,1)+strainCurrent(i,2)+str ainCurrent(i,3)
c
c Trace of current stress tensor
SIGtrace=stressOld(i,1)+stressOld(i,2)+stressOld(i ,3)
c
c Trace of current strain increment tensor
DSTRtrace=strainInc(i,1)+strainInc(i,2)+strainInc( i,3)
c
c Evaluate the trace of the stress increment
DSIGtrace=((H/G)*STRtrace)+(dt/G)*((-zeta1*SIGtrace)+(omega1*STRtrace))
c
C Sress terms
c
do j=1,ndir
c
Term1=D*DSIGtrace
Term2=((phi1*SIGtrace)+(2*psi1*stressOld(i,j))
c
c Strain terms
c
Term3=A*DSTRtrace
Term4=B*strainInc(i,j)
Term5=(lamda1*STRtrace)+(2*mu1*strainCurrent(i,j))
c
c update normal stresses
c
stressNew(i,j)=(1/F)*((stressOld(i,j)+Term3+Term4-Term1+(dt*(Term5-Term2)))
End do
c
c Shear strain components from deformation gradient
c
strainCurrent(i,4)=0.5*((defgradOld(i,1)*defgradOl d(i,4))+(defgradOld(i,7)*defgradOld(i,2))+(defgrad Old(i,6)*defgradOld(i,8)))
strainCurrent(i,5)=0.5*((defgradOld(i,4)*defgradOl d(i,9))+(defgradOld(i,2)*defgradOld(i,5))+(defgrad Old(i,8)*defgradOld(i,3)))
strainCurrent(i,6)=0.5*((defgradOld(i,1)*defgradOl d(i,9))+(defgradOld(i,7)*defgradOld(i,5))+(defgrad Old(i,6)*defgradOld(i,3)))
c
c update shear stress
c
stressNew(i,4)=stressOld(i,4)+((B/F)*strainInc(i,4))+((2*dt/F)*((mu1*strainCurrent(i,4))-(psi1*stressold(i,4))))
stressNew(i,5)=stressOld(i,5)+((B/F)*strainInc(i,5))+((2*dt/F)*((mu1*strainCurrent(i,5))-(psi1*stressold(i,5))))
stressNew(i,6)=stressOld(i,6)+((B/F)*strainInc(i,6))+((2*dt/F)*((mu1*strainCurrent(i,6))-(psi1*stressold(i,6))))
c
c
c Tensile failure criteria
c
SIG1=stressNew(i,1)
if(SIG1.gt.zero)
f=SIG1/f1t
if(f.gt.1)then
stateNew(i,1)=zero*stateOld(i,1)
else stateNew(i,1)=stateOld(i,1)
endif
endif
c
100 continue
return
end
Please let me know if my approach is correct.
I have obtained the current strain components from the deformation gradient
Also,
1.Do i need to update the deformation gradient tensor at the end of each time step?If yes ,how do I do it?
2.I have defined local material orientations in my model.Are the variable tensors like stresses etc provided by abaqus in global cs?If yes how do I change it to local coordinate systems?
3.The failure is assumed to occur only via tensile failure mode.
I have written a VUMAT for a viscoelastic model (I obtained it from one Mr.Frank ritcher's thesis) which best represents fabric behavior in a ballistic impact simulation.I have attached the constitutive relation in the jpg file.
VUMAT
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(*), coordMp(nblock,*),
1 charLength(*), strainInc(nblock,ndir+nshr),
2 relSpinInc(*), tempOld(*),
3 stretchOld(*),
4 defgradOld(nblock,ndir+nshr+nshr),
5 fieldOld(*), stressOld(nblock,ndir+nshr),
6 stateOld(nblock,1), enerInternOld(*),
7 enerInelasOld(*), tempNew(*),
8 stretchNew(*),
8 defgradNew(nblock,ndir+nshr+nshr),
9 fieldNew(*),
1 stressNew(nblock,ndir+nshr), stateNew(nblock,1),
2 enerInternNew(*), enerInelasNew(*),
C
character*80 cmname
c
real,kk1,gk1,kk2,gk2,gv,kv
real,A,B,D,F,G,H
real,phi1,phi2,zeta1,zeta2
real,psi1,psi2,omega1,omega2
real,mu1,mu2,lamda1,lamda2
real,STRtrace,SIGtrace,DSTRtrace,DSIGtrace
real,Term1,Term2,Term3,Term4,Term5
real,SIG1
c
c Model coefficients
c
kk1=9.33E+08
gk1=1.40E+09
kk2=2.27E+10
gk2=3.40E+10
gv=1.00E+06
c kv (infinite-volume conserved flow in dashpot)
f1t=2610
c
omega1=0.0
omega2=1.0+kk2/kk1
psi1=1/(4.0*gv)
psi2=1/(4.0*gk1)
mu1=gk2/(2.0*gv)
mu2=(1.0+gk2/gk1)/2.0
lamda1=(-gk2/gv)/3.0
lamda2=(kk2/kk1-gk2/gk1)/3.0
phi1=-1.0/(6.0*gv)
phi2=1.0/(9.0*kk1)-1.0/(6.0*gk1)
zeta1=0.0
zeta2=1.0/(3.0*kk1)
c
c
A=dt*lamda1+lamda2
B=2.0*mu1*dt+2.0*mu2
D=phi1*dt+phi2
F=2.0*psi1*dt+2.0*psi2
G=zeta1*dt+zeta2
H=omega1*dt+omega2
c
Real,DIMENSION(nblock,ndir+nshr)::strainCurrent
c
save
c
do 100 i= 1,nblock
c
c Trace of current strain tensor from deformation gradient tensor
strainCurrent(i,1)=0.5*((defgradOld(i,1)*defgradOl d(i,1))+(defgradOld(i,7)*defgradOld(i,7))+(defgrad Old(i,6)*defgradOld(i,6))-1)
strainCurrent(i,2)=0.5*((defgradOld(i,4)*defgradOl d(i,4))+(defgradOld(i,2)*defgradOld(i,2))+(defgrad Old(i,8)*defgradOld(i,8))-1)
strainCurrent(i,3)=0.5*((defgradOld(i,9)*defgradOl d(i,9))+(defgradOld(i,5)*defgradOld(i,5))+(defgrad Old(i,3)*defgradOld(i,3))-1)
STRtrace=strainCurrent(i,1)+strainCurrent(i,2)+str ainCurrent(i,3)
c
c Trace of current stress tensor
SIGtrace=stressOld(i,1)+stressOld(i,2)+stressOld(i ,3)
c
c Trace of current strain increment tensor
DSTRtrace=strainInc(i,1)+strainInc(i,2)+strainInc( i,3)
c
c Evaluate the trace of the stress increment
DSIGtrace=((H/G)*STRtrace)+(dt/G)*((-zeta1*SIGtrace)+(omega1*STRtrace))
c
C Sress terms
c
do j=1,ndir
c
Term1=D*DSIGtrace
Term2=((phi1*SIGtrace)+(2*psi1*stressOld(i,j))
c
c Strain terms
c
Term3=A*DSTRtrace
Term4=B*strainInc(i,j)
Term5=(lamda1*STRtrace)+(2*mu1*strainCurrent(i,j))
c
c update normal stresses
c
stressNew(i,j)=(1/F)*((stressOld(i,j)+Term3+Term4-Term1+(dt*(Term5-Term2)))
End do
c
c Shear strain components from deformation gradient
c
strainCurrent(i,4)=0.5*((defgradOld(i,1)*defgradOl d(i,4))+(defgradOld(i,7)*defgradOld(i,2))+(defgrad Old(i,6)*defgradOld(i,8)))
strainCurrent(i,5)=0.5*((defgradOld(i,4)*defgradOl d(i,9))+(defgradOld(i,2)*defgradOld(i,5))+(defgrad Old(i,8)*defgradOld(i,3)))
strainCurrent(i,6)=0.5*((defgradOld(i,1)*defgradOl d(i,9))+(defgradOld(i,7)*defgradOld(i,5))+(defgrad Old(i,6)*defgradOld(i,3)))
c
c update shear stress
c
stressNew(i,4)=stressOld(i,4)+((B/F)*strainInc(i,4))+((2*dt/F)*((mu1*strainCurrent(i,4))-(psi1*stressold(i,4))))
stressNew(i,5)=stressOld(i,5)+((B/F)*strainInc(i,5))+((2*dt/F)*((mu1*strainCurrent(i,5))-(psi1*stressold(i,5))))
stressNew(i,6)=stressOld(i,6)+((B/F)*strainInc(i,6))+((2*dt/F)*((mu1*strainCurrent(i,6))-(psi1*stressold(i,6))))
c
c
c Tensile failure criteria
c
SIG1=stressNew(i,1)
if(SIG1.gt.zero)
f=SIG1/f1t
if(f.gt.1)then
stateNew(i,1)=zero*stateOld(i,1)
else stateNew(i,1)=stateOld(i,1)
endif
endif
c
100 continue
return
end
Please let me know if my approach is correct.
I have obtained the current strain components from the deformation gradient
Also,
1.Do i need to update the deformation gradient tensor at the end of each time step?If yes ,how do I do it?
2.I have defined local material orientations in my model.Are the variable tensors like stresses etc provided by abaqus in global cs?If yes how do I change it to local coordinate systems?
3.The failure is assumed to occur only via tensile failure mode.