You are here
VUMAT of Drucker Prager, help needed please..
Hi Everyone,
I have written Drucker Prager in VUMAT version as shown below..but still there are some deviation with the original in built version of ABAQUS. I have all the notes how i derived these equation.could you please show me where the code has gone wrong or have i missed a critical point in here..i tried my best to spot the error but could not..may be you all experts could find it. Thank you very much.
subroutine vumat(
c Read only -
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 -
7 stressNew, stateNew, enerInternNew, enerInelasNew)
c
include 'vaba_param.inc'
C
dimension props(nprops), density(nblock),
1 coordMp(nblock,*), nElement(nblock), nMatPoint(nblock),
2 charLength(*), strainInc(nblock,ndir+nshr),
3 relSpinInc(*), tempOld(*),nLayer(nblock),
4 stretchOld(*), defgradOld(*),nSecPoint(nblock),
5 fieldOld(*), stressOld(nblock,ndir+nshr),
6 stateOld(nblock,nstatev), enerInternOld(nblock),
7 enerInelasOld(nblock), tempNew(*),
8 stretchNew(*), defgradNew(*),
9 fieldNew(*), stressNew(nblock,ndir+nshr),
1 stateNew(nblock,nstatev),
2 enerInternNew(nblock), enerInelasNew(nblock)
C
character*80 cmname
C
parameter(zero = 0.d0, one = 1.d0, two = 2.d0, half = .5d0,
& Three = 3.d0, op5 = 1.5d0,third = one/Three)
C Assume purely elastic material at the beginning of the analysis
C
E = props(1)
xnu = props(2)
beta=props(3)
sigmaT=props(4)
K=props(5)
yields1= props(6)
eqpl1 = props(7)
yields2= props(8)
eqpl2 = props(9)
twomu = E /(one + xnu)
alambda = xnu*twomu/(one-two*xnu)
thremu = op5 * twomu
term1=one/K
nvalue = (nprops/three)-1
if (stepTime.eq.zero) then
do k = 1,nblock
trace=strainInc(k,1)+strainInc(k,2)+strainInc(k,3)
stressNew(k,1)=stressOld(k,1)
& + twomu*strainInc(k,1)+ alambda*trace
stressNew(k,2)=stressOld(k,2)
& + twomu*strainInc(k,2)+ alambda*trace
stressNew(k,3)=stressOld(k,3)
& + twomu*strainInc(k,3)+ alambda*trace
stressNew(k,4)=stressOld(k,4)+twomu*strainInc(k,4)
if (nshr.gt.1) then
stressNew(k,5)=stressOld(k,5)+twomu*strainInc(k,5)
stressNew(k,6)=stressOld(k,6)+twomu*strainInc(k,6)
end if
end do
else
do k=1,nblock
c
peeqOld = stateOld(k,1)
call vuhard(yieldOld, hard, peeqOld, props(6), nvalue)
c
trace=strainInc(k,1)+strainInc(k,2)+strainInc(k,3)
c
C Trial Stress
s11=stressOld(k,1)+twomu*strainInc(k,1)
& + alambda*trace
s22=stressOld(k,2)+twomu*strainInc(k,2)
& + alambda*trace
s33=stressOld(k,3)+twomu*strainInc(k,3)
& + alambda*trace
s12=stressOld(k,4)+twomu*strainInc(k,4)
if (nshr.gt.1) then
s13=stressOld(k,5)+twomu*strainInc(k,5)
s23=stressOld(k,6)+twomu*strainInc(k,6)
end if
C
smean = (one/Three)*(s11+s22+s33)
C
s11 = s11 - smean
s22 = s22 - smean
s33 = s33 - smean
C
if (nshr.eq.1) then
vmises = sqrt(op5*(s11*s11+s22*s22+s33*s33
& +two*s12*s12))
else
vmises = sqrt(op5*(s11*s11+s22*s22+s33*s33
& +two*s12*s12+two*s13*s13+two*s23*s23))
end if
t=vmises
d=sigmaT*(term1+(third*tan(beta)))
DPYield=t-(-smean*tan(beta))-d
a=term1+(third*tan(beta))
c
sigdif=DPYield-yieldOld
facyld=zero
if (sigdif.le.zero)then
factor =one
c deqps=facyld*sigdif/(thremu+hard)
deqps=(sqrt(two)/three)*sqrt(((dps11-dps22)**two)+
& ((dps22-dps33)**two)+((dps11-dps33)**two)+six*
& ((dps12**two)+(dps23**two)+(dps13**two)))
yieldNew = yieldOld + hard * deqps
factor = yieldNew/(yieldNew + thremu * deqps)
else
c facyld=one
c deqps=facyld*sigdif/(thremu+hard)
c equivalant plastic strain expression
deqps=(sqrt(two)/three)*sqrt(((dps11-dps22)**two)+
& ((dps22-dps33)**two)+((dps11-dps33)**two)+six*
& ((dps12**two)+(dps23**two)+(dps13**two)))
c plastic strain components for Drucker Prager criterion ( derived these equations from the flow rule shown in abaqus section 4.4.2)
dps11=((deqps/(a*two*vmises))*((two*s11)-s33-s22)
& - third*tan(beta))
dps22=((deqps/(a*two*vmises))*((two*s22)-s11-s33)
& - third*tan(beta))
dps33=((deqps/(a*two*vmises))*((two*s33)-s11-s22)
& - third*tan(beta))
dps12=(deqps/(a*two*vmises))*six*s12
dps13=(deqps/(a*two*vmises))*six*s13
dps23=(deqps/(a*two*vmises))*six*s23
c
yieldNew = yieldOld + hard * deqps
factor = yieldNew/(yieldNew + thremu * deqps)
end if
C update the stress
stressNew(k,1) = s11 * factor + smean
stressNew(k,2) = s22 * factor + smean
stressNew(k,3) = s33 * factor + smean
stressNew(k,4) = s12 * factor
if (nshr.gt.1)then
stressNew(k,5) = s13 * factor
stressNew(k,6) = s23 * factor
end if
C Update the state variables
stateNew(k,1) = stateOld(k,1) + deqps
C Update the specific internal energy
if (nshr.eq.1) then
stressPower = half * (
& ( stressOld(k,1) + stressNew(k,1) ) * strainInc(k,1) +
& ( stressOld(k,2) + stressNew(k,2) ) * strainInc(k,2) +
& ( stressOld(k,3) + stressNew(k,3) ) * strainInc(k,3))+
& ( stressOld(k,4) + stressNew(k,4) ) * strainInc(k,4)
else
stressPower = half * (
& ( stressOld(k,1) + stressNew(k,1) ) * strainInc(k,1) +
& ( stressOld(k,2) + stressNew(k,2) ) * strainInc(k,2) +
& ( stressOld(k,3) + stressNew(k,3) ) * strainInc(k,3))+
& ( stressOld(k,4) + stressNew(k,4) ) * strainInc(k,4) +
& ( stressOld(k,5) + stressNew(k,5) ) * strainInc(k,5) +
& ( stressOld(k,6) + stressNew(k,6) ) * strainInc(k,6)
end if
C
enerInternNew(k) = enerInternOld(k) + stressPower/density(k)
c
C Update the dissipated inelastic specific energy
C
plasticWorkInc = half*(yieldOld + yieldNew)*deqps
c
enerInelasNew(k) = enerInelasOld(k)
& + plasticWorkInc / density(k)
end do
end if
return
end
c
c
subroutine vuhard(syield, hard, eqplas, table, nvalue)
include 'vaba_param.inc'
dimension table(2,nvalue)
parameter(zero=0.d0)
C set yield stress to last value of table, hardening to zero
C
syield=table(1, nvalue)
hard=zero
C
C if more than one entry, search table
if (nvalue.gt.1)then
do k1=1,nvalue-1
eqpl1 = table(2,k1+1)
if(eqplas.lt.eqpl1)then
eqpl0 = table(2,k1)
C yield stress and hardening
deqpl=eqpl1-eqpl0
syield0=table(1,k1)
syield1=table(1,k1+1)
dsyiel=syield1-syield0
hard=dsyiel/deqpl
syield=syield0+(eqplas-eqpl0)*hard
goto 40
end if
end do
40 continue
end if
return
end
- Gayan Aravinda's blog
- Log in or register to post comments
- 9464 reads
Comments
Dear user
Dear user
you had a post before and I answered your post but you didn't reply that post, Plase let other know you have solved your problem...
Sometimes "small" deviation in your answere and other Fe codes are quiet normal, I will tell you my experience ...
I wrote a full FE code in Fortran two years ago(Drucker prager model), and I compared my code by ansys and I found some deviation...
My code was correct but the deviation was because of the number 1 in below. there may be many reason for deviation:
1- Ansys is using different shape function for elements(extra shape functions) ( I used classical shape function)
2- Ansys may use different number of integration point
3- Ansys may use B_bar method for calculation of volumetric constrains...
After I changed the setting of the ansys to have same option as my code, my result became accurate up to four floating point.
In your case may be the same happens in abaqus, may be when you use Vumat the abaqus disables some option for elements...
The first step I suggest is to set initial yield condition vey high in order that plastic flow wont happen... then check your elastic answer with the one in original built-in abaqus... this helps you to make sure the deviation is not because of element options ....
Hello
Hi, Aslan...thank you very much...i may forgotten soory about it...i solved the VUMAT for perfect plastic case.....thank you very much for the suggestions ..i will have a look at it...
VUMAT code
Hi Gayan Aravinda,
I am new with VUMAT.
I kindly ask you for the VUMAT in the corrected version if possible.
I really would appreciate it.
ciprian.zub@student.upt.ro