| Newsgroups |
gmane.comp.mathematics.abaqus.user |
| Message-ID |
<[email protected]> |
Hello everyone,
I am writing a VUMAT of an elasto-plastic model for modified Mohr-Coulomb criterion with the isotropic hardening of cohesion based on the smooth hyperbolic approximation of Abbo-Sloan paper : "A smooth hyperbolic approximation to the Mohr-Coulomb yield criterion"
I am using semi-implicit stress itergration scheme. Actually, it could not get the convergence after entering the plastic step. Could anyone help me check my code? I have the code in UMAT but I want to use in explicit way so I have to convert into VUMAT.
Thank you in advance.
Here is my code
-------------------------------------------------------------------------------------------
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),defgradOld(nblock,ndir+nshr+nshr),
4 fieldOld(nblock, nfieldv), stressOld(nblock, ndir+nshr),
5 stateOld(nblock, nstatev), enerInternOld(nblock),
6 enerInelasOld(nblock), tempNew(nblock),
7 stretchNew(nblock, ndir+nshr),defgradNew(nblock,ndir+nshr+nshr),
8 fieldNew(nblock, nfieldv), stressNew(nblock,ndir+nshr),
9 stateNew(nblock, nstatev), enerInternNew(nblock),
1 enerInelasNew(nblock)
c
character*80 cmname
C
C
dimension devias(6), pstres(6), ce(6,6), plas(6),
+ dstran(6), dfds(6), dgds(6), ce_n(6)
real devias,pstres,sbar,invar2,invar3, ce, plas,
+ dfds, dgds, ce_n
c
parameter ( zero = 0.d0, one = 1.d0, two = 2.d0, three = 3.0d0,
* six = 6.0d0, enumax = 0.499e0, pi = 3.1415927,
* third = 1.d0 / 3.d0, half = 0.5d0, op5 = 1.5d0)
C
C Only for 3D cases using
C Mohr-Coulomb Plasticity with isotropic hardening of cohesion.
C
C The state variable is stored as:
C
C STATE(*,1) = equivalent plastic strain
C
C User needs to input
C props(1) Young's modulus, E
C props(2) Poisson's ratio, nu
C props(3) Friction angle, phi
C props(4) Dilatancy angle, pusi
C calls ahard for curve of cohesion vs. plastic strain
C
e = props(1)
xnu = props(2)
phi = props(3)
pusi = props(4)
c
twomu = e / ( one + xnu ) !2G
alamda = xnu * twomu / ( one - two * xnu )
thremu = op5 * twomu !3K
c
nvalue = (nprops-4)/2
ntens = ndir + nshr
c begin the loop
do 900 k = 1, nblock
eqplas = stateOld(k,1)
do k1 = 1, ntens
plas(k1) = stateOld(k,k1+1)
end do
c calculate elastic stiffness
do k1 = 1, ndir
do k2 = 1, ndir
ce(k2,k1) = alamda
end do
ce(k1,k1) = twomu + alamda
end do
do k1 = ndir+1, ntens
ce(k1,k1) = twomu
end do
c trial stress
pstres(:) = stressOld(k,:)
dstran(:) = strainInc(k,:)
do k1 = 1, ntens
do k2 = 1, ntens
pstres(k2) = pstres(k2) + ce(k2,k1)*dstran(k1)
end do
end do
do k1 = 1, ntens
dfds(k1) = zero
dgds(k1) = zero
end do
c calculate Invariants and parameters for yield and potential function
call ahard(yieldOld, hard, eqplas, props(5), nvalue)
yieldd = yieldOld
sigmam = third * ( pstres(1) + pstres(2) + pstres(3) )
devias(1) = pstres(1) - sigmam
devias(2) = pstres(2) - sigmam
devias(3) = pstres(3) - sigmam
invar2 = ( devias(1)**2 + devias(2)**2 + devias(3)**2 )/two
* + pstres(4)**2 + pstres(5)**2 + pstres(6)**2
sbar = sqrt(invar2)
invar3 = devias(1)*devias(2)*devias(3)+
* two*pstres(4)*pstres(5)*pstres(6)-devias(1)*pstres(5)**2-
* devias(2)*pstres(6)**2-devias(3)*pstres(4)**2
s3ta=-op5*sqrt(three)*invar3/(sbar**3)
if (s3ta .gt. 1.0) then
s3ta = 1.0
else if (s3ta .lt. -1.0) then
s3ta = -1.0
end if
theta=asin(s3ta)/three
thetat=(25.0/180.0)*pi
a1=1.432052062044e0
a2=0.406941858374e0
b1=0.544290524902e0
b2=0.673903324498e0
c calculate yield function
parama = 0.05 * yieldd
if (abs(theta).le. thetat) then
kf = cos(theta) - sin(phi)*sin(theta)/sqrt(three)
else
sgn = sign(one,theta)
const1 = a1 + a2*sgn*sin(phi)
const2 = b1*sgn + b2*sin(phi)
kf = const1 - const2*s3ta
end if
fb = sigmam*sin(phi)+
* sqrt( sbar*sbar*kf*kf + parama*parama*cos(phi)*cos(phi) ) -
* yieldd*cos(phi)
c calculate potential function
paramb = 0.2 * yieldd
if (abs(theta).lt.thetat) then
kg=cos(theta)-sin(pusi)*sin(theta)/sqrt(three)
else
sgn=sign(one,theta)
consta=a1+a2*sgn*sin(pusi)
constb=b1*sgn+b2*sin(pusi)
kg=consta-constb*sin(s3ta)
end if
gb = sigmam*sin(pusi)+
* sqrt( sbar*sbar*kg*kg + paramb*paramb*cos(pusi)*cos(pusi) )
iflag=0
c If fb>0 : Plasticity
c If fb<0 : Elasticity
if (fb .ge. 0.0d0) then
call ahard(yieldOld, hard, eqplas, props(5), nvalue)
yieldd = yieldOld
c Calculate the 1st derivative of yield function
alfa=sbar*kf
alfa=alfa/sqrt(alfa**2+parama**2*cos(phi)**2)
if (abs(theta).le.thetat) then
dkf=-sin(theta)-sin(phi)*cos(theta)/sqrt(three)
else
dkf=-3*const2*cos(three*theta)
end if
c1=sin(phi)
c2=kf-tan(three*theta)*dkf
c2=alfa*c2
if (abs(theta).le.thetat) then
c2=alfa*c2
c3=-sqrt(three)*dkf/(two*cos(three*theta)*sbar**2)
c3=alfa*c3
else
c2=kf
c2=alfa*c2
c3=zero
end if
dfds(1)=c1/three+c2*devias(1)/(two*sbar)+
1 c3*(devias(2)*devias(3)-pstres(5)**2+sbar**2/three)
dfds(2)=c1/three+c2*devias(2)/(two*sbar)+
1 c3*(devias(1)*devias(3)-pstres(6)**2+sbar**2/three)
dfds(3)=c1/three+c2*devias(3)/(two*sbar)+
1 c3*(devias(1)*devias(2)-pstres(4)**2+sbar**2/three)
dfds(4)= +c2*pstres(4)/sbar+
1 c3*two*(pstres(5)*pstres(6)-devias(3)*pstres(4))
dfds(5)= +c2*pstres(6)/sbar+
1 c3*two*(pstres(4)*pstres(5)-devias(2)*pstres(6))
dfds(6)= +c2*pstres(5)/sbar+
1 c3*two*(pstres(6)*pstres(4)-devias(1)*pstres(5))
c Calculate the 1st derivative of potential function
alfag=sbar*kg
alfag=alfag/sqrt(alfag**2+paramb**2*cos(pusi)**2)
if (abs(theta.lt.thetat)) then
dkg=-sin(theta)-sin(pusi)*cos(theta)/sqrt(three)
else
dkg=-3*constb*cos(three*theta)
end if
gc1=sin(pusi)
if (abs(theta).le.thetat) then
gc2=kg-tan(three*theta)*dkg
gc2=alfag*gc2
gc3=-sqrt(three)*dkg/(two*cos(three*theta)*sbar**2)
gc3=alfag*gc3
else
gc2=kg
gc2=alfag*gc2
gc3=zero
end if
dgds(1)=gc1/three+gc2*devias(1)/(two*sbar)+
1 gc3*(devias(2)*devias(3)-pstres(5)**2+sbar**2/three)
dgds(2)=gc1/three+gc2*devias(2)/(two*sbar)+
1 gc3*(devias(1)*devias(3)-pstres(6)**2+sbar**2/three)
dgds(3)=gc1/three+gc2*devias(3)/(two*sbar)+
1 gc3*(devias(1)*devias(2)-pstres(4)**2+sbar**2/three)
dgds(4)= +gc2*pstres(4)/sbar+
1 gc3*two*(pstres(5)*pstres(6)-devias(3)*pstres(4))
dgds(5)= +gc2*pstres(6)/sbar+
1 gc3*two*(pstres(4)*pstres(5)-devias(2)*pstres(6))
dgds(6)= +gc2*pstres(5)/sbar+
1 gc3*two*(pstres(6)*pstres(4)-devias(1)*pstres(5))
201 continue
c Check Newton-Ralphson convergence
iflag=iflag+1
if (abs(fb) .lt. yieldd*1e-6) then
go to 301
end if
if (iflag. ge. 500) then
write(*,*) 'divergence', iflag, fb
stop
go to 301
end if
c
do k1=1,ntens
do k2=1,ntens
ce_n(k2)=dgds(k2)+ce(k2,k1)*dgds(k1)
end do
end do
c
b=(dgds(1)-dgds(2))**2+
1 (dgds(2)-dgds(3))**2+(dgds(3)-dgds(1))**2
do k1=ndir+1,ntens
b=b+six*dgds(k1)**2
end do
b=sqrt(two*b)/three
nblfm2=hard*b
c
do k1=1,ntens
nblfm1=nblfm1+dfds(k1)*ce_n(k1)
end do
c
dnbl=fb/(nblfm1+nblfm2)
eqplas=eqplas+dnbl
do k1=1,ntens
pstres(k1)=pstres(k1)-dnbl*ce_n(k1)
plas(k1)=plas(k1)+dnbl*ce_n(k1)
end do
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
call ahard(yieldOld, hard, eqplas, props(5), nvalue)
yieldd = yieldOld
sigmam = third * ( pstres(1) + pstres(2) + pstres(3) )
devias(1) = pstres(1) - sigmam
devias(2) = pstres(2) - sigmam
devias(3) = pstres(3) - sigmam
invar2 = ( devias(1)**2 + devias(2)**2 + devias(3)**2 )/two
* + pstres(4)**2 + pstres(5)**2 + pstres(6)**2
sbar = sqrt(invar2)
invar3 = devias(1)*devias(2)*devias(3)+
* two*pstres(4)*pstres(5)*pstres(6)-devias(1)*pstres(5)**2-
* devias(2)*pstres(6)**2-devias(3)*pstres(4)**2
s3ta=-op5*sqrt(three)*invar3/(sbar**3)
if (s3ta .gt. 1.0) then
s3ta = 1.0
else if (s3ta .lt. -1.0) then
s3ta = -1.0
end if
theta=asin(s3ta)/three
thetat=(25.0/180.0)*pi
a1=1.432052062044e0
a2=0.406941858374e0
b1=0.544290524902e0
b2=0.673903324498e0
c calculate yield function
parama = 0.05 * yieldd
if (abs(theta).le. thetat) then
kf = cos(theta) - sin(phi)*sin(theta)/sqrt(three)
else
sgn = sign(one,theta)
const1 = a1 + a2*sgn*sin(phi)
const2 = b1*sgn + b2*sin(phi)
kf = const1 - const2*s3ta
end if
fb = sigmam*sin(phi)+
* sqrt( sbar*sbar*kf*kf + parama*parama*cos(phi)*cos(phi) ) -
* yieldd*cos(phi)
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
go to 201
301 continue
end if !End check plasticity
stressNew (k,:) = pstres(:)
stateNew (k,1) = eqplas
do k1=1, ntens
stateNew (k,2:7) = plas(k1)
end do
c
900 continue
return
end
c
c Subroutine ahard
c
subroutine ahard(syield, hard, eqplas, table, nvalue)
include 'vaba_param.inc'
c
dimension table(2, nvalue)
c
parameter(zero=0.d0)
c
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
c
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
c yield stress and hardening
c
deqpl=eqpl1-eqpl0
syiel0=table(1, k1)
syiel1=table(1, k1+1)
dsyiel=syiel1-syiel0
hard=dsyiel/deqpl
syield=syiel0+(eqplas-eqpl0)*hard
goto 10
endif
end do
10 continue
endif
return
end subroutine ahard