PDA

View Full Version : VUMAT for maxwell zener model-validation



sirig
2009-10-16, 15:55
Hello all ,

Here is the VUMAT I wrote for a standard linear solid (maxwell-zener model obtained from Mr.Frank Richter's thesis) .When I tried to test it on a single element ,though the code runs fine giving accurate values of stresses etc. , the stress and strain curves are not smooth.Attached is my code ,the input file I used and the plots for stresses and strains .

Please let me know if it is a valid response from a viscoelastic material.


Best Regards,
Sireesha

VUMAT


subroutine vumat(nblock,ndir,nshr,nstatev,nfieldv,nprops,lann eal,
1 stepTime,totalTime,dt,cmname,coordMp,charLength,pr ops,
2 density,strainInc,relSpinInc,tempOld,stretchOld,de fgradOld,
3 fieldOld,stressOld,stateOld,enerInternOld,enerInel asOld,tempNew,
4 stretchNew,defgradNew,fieldNew,stressNew,stateNew, enerInternNew,
5 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),fieldOld(nblock,nfiel dv),
4 defgradOld(nblock,ndir+nshr+nshr),stressOld(nblock ,ndir+nshr),
5 stateOld(nblock,nstatev),enerInternOld(nblock),
6 enerInelasOld(nblock),tempNew(nblock),enerInelasNe w(nblock),
7 stretchNew(nblock,ndir+nshr),defgradNew(nblock,ndi r+nshr+nshr),
8 fieldNew(nblock,nfieldv),stressNew(nblock,ndir+nsh r),
9 stateNew(nblock,nstatev),enerInternNew(nblock)
c
implicit none
integer i,j,j1,k,k1,statusMP
dimension StrainCurrent(nblock,ndir+nshr)
real A,B,D,F,G,H,lambda,checkf
real kk1,gk1,kk2,gk2,gv,E,nu,f1t
real omega1,omega2,psi1,psi2,mu1,mu2
real lambda1,lambda2,phi1,phi2,zeta1,zeta2,BM,mu
real STRtrace,SIGtrace,DSTRtrace,DSIGtrace,S1
C
C
character*80 cmname
c
C parameter(zero=0.,one=1.,two=2.,three=3.,four=4.,s ix=6.,nine=9.,half=.5,
C 1 third=one/three,sixth=one/six,ninth=one/nine,fourth=one/four)
c
c
c Model coefficients
c
kk1=props(1)
gk1=props(2)
kk2=props(3)
gk2=props(4)
gv=props(5)
E=props(6)
nu=props(7)
c kv (infinite-volume conserved flow in dashpot)
f1t=props(8)
c
omega1=0.0
omega2=1.0+kk2/kk1
psi1=1.0/(4.0*gv)
psi2=1.0/(4.0*gk1)
mu1=gk2/(2.0*gv)
mu2=0.5*(1.0+gk2/gk1)
lambda1=(-gk2/gv)/3.0
lambda2=(1.0/3.0)*(kk2/kk1-gk2/gk1)
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)
BM=E/(3.0*(1.0-2.0*nu))
mu=E/(2.0*(1.0+nu))
c
A=dt*lambda1+lambda2
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
lambda=BM-(2.0*mu/3.0)
C
c
if(totalTime.EQ.0) then
do k = 1,nblock
trace=strainInc(k,1)+strainInc(k,2)+strainInc(k,3)
C update normal stresses
stressNew(k,1)=stressOld(k,1)+(lambda*trace)
1 + (2.0*mu*strainInc(k,1))
stressNew(k,2)=stressOld(k,2)+(lambda*trace)
1 + (2.0*mu*strainInc(k,2))
stressNew(k,3)=stressOld(k,3)+(lambda*trace)
1 + (2.0*mu*strainInc(k,3))
C update shear stresses
stressNew(k,4)=stressOld(k,4)
1 + (2.0*mu*strainInc(k,4))
stressNew(k,5)=stressOld(k,5)
1 + (2.0*mu*strainInc(k,5))
stressNew(k,6)=stressOld(k,6)
1 + (2.0*mu*strainInc(k,6))
end do
else

do 100 i = 1,nblock
c
c current normal strain tensor components from deformation gradient
c
StrainCurrent(i,1)=0.5*(defgradNew(i,1)*defgradNew (i,1))
1 + (0.5*(defgradNew(i,7)*defgradNew(i,7)))
1 + (0.5*((defgradNew(i,6)*defgradNew(i,6))-1))
StrainCurrent(i,2)=0.5*(defgradNew(i,4)*defgradNew (i,4))
1 + (0.5*(defgradNew(i,2)*defgradNew(i,2)))
1 + (0.5*((defgradNew(i,8)*defgradNew(i,8))-1))
StrainCurrent(i,3)=0.5*(defgradNew(i,9)*defgradNew (i,9))
1 + (0.5*(defgradNew(i,5)*defgradNew(i,5)))
1 + (0.5*((defgradNew(i,3)*defgradNew(i,3))-1))
c
c current shear strain tensor components from deformation gradient
c
StrainCurrent(i,4)=0.5*((defgradNew(i,1)*defgradNe w(i,4))
1 + (defgradNew(i,7)*defgradNew(i,2))
1 + (defgradNew(i,6)*defgradNew(i,8)))
StrainCurrent(i,5)=0.5*((defgradNew(i,4)*defgradNe w(i,9))
1 + (defgradNew(i,2)*defgradNew(i,5))
1 + (defgradNew(i,8)*defgradNew(i,3)))
StrainCurrent(i,6)=0.5*((defgradNew(i,1)*defgradNe w(i,9))
1 + (defgradNew(i,7)*defgradNew(i,5))
1 + (defgradNew(i,6)*defgradNew(i,3)))
C
statusMp=stateNew(i,1)
c
STRtrace=StrainCurrent(i,1)+StrainCurrent(i,2)
1 + StrainCurrent(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)
1 + (omega1*STRtrace))
c
C Stress terms
c
c Update normal stresses
c
if(statusMP.EQ.1.0) then
do j=1,3
stressNew(i,j)= stressOld(i,j)+(A/F)*DSTRtrace
1 + (B/F)*strainInc(i,j)-(D/F)*DSIGtrace
2 + (dt/F)*((lamda1*STRtrace)+(2.0*mu1*StrainCurrent(i,j))
3 - (phi1*SIGtrace)-(2.0*psi1*stressOld(i,j)))
c Update shear stress
stressNew(i,j+3)=stressOld(i,j+3)+(B/F)*strainInc(i,j+3)
1 + ((2.0*dt)/F)*((mu1*StrainCurrent(i,j+3))
2 - (psi1*stressOld(i,j+3)))
end do
else
do j1=1,6
stressNew(i,j1)=0.0
end do
endif

c
c Tensile failure criteria
c
S1=stressNew(i,3)
checkf=S1/f1t
if(checkf.GE.1) then
stateNew(i,1)=0.0
else
stateNew(i,1)=1.0
endif
c
100 continue
endif
c
return
end


INPUT FILE

*Heading
** Job name: oct15 Model name: Model-1
*Preprint, echo=NO, model=NO, history=NO, contact=NO
**
** PARTS
**
*Part, name=Part-1
*Node
1, 20., 20., 40.
2, 20., -20., 40.
3, 20., 20., 0.
4, 20., -20., 0.
5, -20., 20., 40.
6, -20., -20., 40.
7, -20., 20., 0.
8, -20., -20., 0.
*Element, type=C3D8R
1, 5, 6, 8, 7, 1, 2, 4, 3
*Nset, nset=_PickedSet5, internal, generate
1, 8, 1
*Elset, elset=_PickedSet5, internal
1,
** Section: Section-1
*Solid Section, elset=_PickedSet5, controls=EC-1, material=Material-1
1.,
*End Part
**
**
** ASSEMBLY
**
*Assembly, name=Assembly
**
*Instance, name=Part-1-1, part=Part-1
*End Instance
**
*Nset, nset=_PickedSet13, internal, instance=Part-1-1
3, 4, 7, 8
*Elset, elset=__PickedSurf12_S3, internal, instance=Part-1-1
1,
*Surface, type=ELEMENT, name=_PickedSurf12, internal
__PickedSurf12_S3, S3
*End Assembly
**
** ELEMENT CONTROLS
**
*Section Controls, name=EC-1, ELEMENT DELETION=YES
1., 1., 1.
**
** MATERIALS
**
*Material, name=Material-1
*Density
1.44e-09,
*Depvar, delete=1
1,
*User Material, constants=8
933., 1400.,22700.,34000., 1.,90000., 0.3, 2610.
** ----------------------------------------------------------------
**
** STEP: applyload
**
*Step, name=applyload
*Dynamic, Explicit
, 0.0005
*Bulk Viscosity
0.06, 1.2
**
** BOUNDARY CONDITIONS
**
** Name: BC-2 Type: Displacement/Rotation
*Boundary
_PickedSet13, 1, 1
_PickedSet13, 2, 2
_PickedSet13, 3, 3
_PickedSet13, 4, 4
_PickedSet13, 5, 5
_PickedSet13, 6, 6
**
** LOADS
**
** Name: Load-1 Type: Surface traction
*Dsload
_PickedSurf12, TRVEC, 1300., 0., 0., 1.
**
** OUTPUT REQUESTS
**
*Restart, write, number interval=1, time marks=NO
**
** FIELD OUTPUT: F-Output-1
**
*Output, field
*Node Output
A, RF, U, V
*Element Output, directions=YES
LE, PE, PEEQ, S, STATUS
*Contact Output
CSTRESS,
**
** HISTORY OUTPUT: H-Output-1
**
*Output, history, variable=PRESELECT
*End Step

Jorgen
2009-10-17, 10:57
Since you are using an explicit simulation the resulting stress & strain values are often noisy. Are you sure you are running the Explicit anlysis correctly?

-Jorgen