Generation237
Student
- Jun 12, 2023
- 2
Hi,
i am having a problem with my vumat subroutine. It is based on a constitutive model for full range elastoplastic behavior. You can see the code below. My problem is that at the end of a simulation, the abaqus displays zero stress (zero for the smises,s11,s22,etc) everything thing after if the "if-loop" becomes zero. But once i delete the plastic part, abaqus does display the new stresses for the elastic part of the code. Do you have an idea for why that is? I really need your help.
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
C The state variables are stored as:
C state(*,1) = sback stress component 11 ! short-range backstress
C state(*,2) = sback stress component 22
C state(*,3) = sback stress component 33
C state(*,4) = sback stress component 12
C state(*,5) = sback stress component 13
C state(*,6) = sback stress component 23
C state(*,7) = equivalent plastic strain
C state(*,8) = center of the memory surface 11
C state(*,9) = center of the memory surface 22
C state(*,10) = center of the memory surface 33
C state(*,11) = center of the memory surface 12
C state(*,12) = center of the memory surface 13
C state(*,13) = center of the memory surface 23
C state(*,14) = radius of the memory surface
C state(*,15) = plastic strain 11
C state(*,16) = plastic strain 22
C state(*,17) = plastic strain 33
C state(*,18) = plastic strain 12
C state(*,19) = plastic strain 13
C state(*,20) = plastic strain 23
C state(*,21) = elastic strain 11
C state(*,22) = elastic strain 22
C state(*,23) = elastic strain 33
C state(*,24) = elastic strain 12
C state(*,25) = elastic strain 13
C state(*,26) = elastic strain 23
C state(*,27) = lback stress component 11 ! long-range backstress
C state(*,28) = lback stress component 22
C state(*,29) = lback stress component 33
C state(*,30) = lback stress component 12
C state(*,31) = lback stress component 13
C state(*,32) = lback stress component 23
C state(*,33) = yield index ! predefined variable
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),
4 defgradOld(nblock,ndir+nshr+nshr),
5 fieldOld(nblock,nfieldv), stressOld(nblock,ndir+nshr),
6 stateOld(nblock,nstatev), enerInternOld(nblock),
7 enerInelasOld(nblock), tempNew(nblock),
8 stretchNew(nblock,ndir+nshr),
9 defgradNew(nblock,ndir+nshr+nshr),
1 fieldNew(nblock,nfieldv),
2 stressNew(nblock,ndir+nshr), stateNew(nblock,nstatev),
3 enerInternNew(nblock), enerInelasNew(nblock)
C
character*80 cmname
C
parameter ( zero = 0.d0, one = 1.d0, two = 2.d0, three = 3.d0,
1 third = one/three, half = .5d0, twoThirds = two/three,
2 threeHalfs = three/two, six = 6.d0)
C
C Elastic constants
e = props(1) ! E-Mod
xnu = props(2) ! Poisson's ratio
yield0 = props(3) ! initial yield stress
epst = props(4) ! plastic strain at the end of the yield plateau in monotonic curve
epst_amp = props(5) ! threshold amplitude of plastic strain
cskin1 = props(6) ! short-range kinematic hardening parameter, c_s
cskin2 = props(7)
clkin1 = props(8) ! long-range kinematic hardening parameter, c_l
clkin2 = props(9)
sgamma1 = props(10) ! short-range kinematic hardening parameter, gamma_s
sgamma2 = props(11)
lgamma1 = props(12) ! long-range kinematic hardening parameter, gamma_l
lgamma2 = props(13)
qishs = props(14) ! short-range saturated value of the isotropic hardening component, Q
qishl = props(15) ! long-range
bsats = props(16) ! short-range rate in which the saturation is achieved, b_s
bsatl = props(17) ! long-range rate, b_l
cs_para = props(18) ! short-range scalar parameter, c_s
cl_para = props(19) ! long-range, c_l
C
mu = e / ( two * ( one + xnu ) ) ! mu (lamé constant) = G (shear modulus)
twomu = e / ( one + xnu ) ! 2G = 2mu
thremu = three * e / ( two * ( one + xnu ) ) ! 3G = 3mu
C alamda = twomu * ( e - twomu ) / ( sixmu - two * e ) ! lamé constant
xkap = e / ( three * ( one - two*xnu ) )
con = xnu / e
coef1 = sqrt(twoThirds)
coef2 = sqrt(threeHalfs)
C Deviatoric strain increment
do i = 1, nblock
trace = strainInc(i,1) + strainInc(i,2) + strainInc(i,3)
de1 = strainInc(i,1) - ( third * trace )
de2 = strainInc(i,2) - ( third * trace )
de3 = strainInc(i,3) - ( third * trace )
de4 = strainInc(i,4)
de5 = strainInc(i,5)
de6 = strainInc(i,6)
es1 = de5
es5 = trace
C Trial elastic Deviatoric stress incremment
s1_inc = twomu * de1
s2_inc = twomu * de2
s3_inc = twomu * de3
s4_inc = twomu * de4
s5_inc = twomu * de5
s6_inc = twomu * de6
es2 = twomu
es3 = s5_inc
C Predictor/Trial stress
sig1 = stressOld(i,1) + ( xkap*trace ) + s1_inc
sig2 = stressOld(i,2) + ( xkap*trace ) + s2_inc
sig3 = stressOld(i,3) + ( xkap*trace ) + s3_inc
sig4 = stressOld(i,4) + s4_inc
sig5 = stressOld(i,5) + s5_inc
sig6 = stressOld(i,6) + s6_inc
str1 = sig5
es4 = xkap
C Trial elastic and plastic strain
elas1 = stateOld(i,21) ! Elastic strain
elas2 = stateOld(i,22)
elas3 = stateOld(i,23)
elas4 = stateOld(i,24)
elas5 = stateOld(i,25)
elas6 = stateOld(i,26)
eltr1 = elas1 + strainInc(i,1) ! Trial elastic strain
eltr2 = elas2 + strainInc(i,2)
eltr3 = elas3 + strainInc(i,3)
eltr4 = elas4 + strainInc(i,4)
eltr5 = elas5 + strainInc(i,5)
eltr6 = elas6 + strainInc(i,6)
pls1 = stateOld(i,15) ! plastic strain
pls2 = stateOld(i,16)
pls3 = stateOld(i,17)
pls4 = stateOld(i,18)
pls5 = stateOld(i,19)
pls6 = stateOld(i,20)
if( stepTime .eq. zero ) then
stressNew(i,1) = sig1
stressNew(i,2) = sig2
stressNew(i,3) = sig3
stressNew(i,4) = sig4
stressNew(i,5) = sig5
stressNew(i,6) = sig6
else
C Trace of stress tensor
smean = third * ( sig1 + sig2 + sig3 )
str2 = smean
C Deviatoric stress from the backstress
s1 = sig1 - smean
s2 = sig2 - smean
s3 = sig3 - smean
s4 = sig4
s5 = sig5
s6 = sig6
C Equivalent plastic strain
eqps0 = stateOld(i,7) ! initial value should be zero
rad_ms0 = stateOld(i,14)
C Von Mises equivalent stress
smises = s1**2 + s2**2 + s3**2 + two*(s4**2 + s5**2 + s6**2)
smises = sqrt(threeHalfs*smises)
C Radius of yield
radpr0 = yield0 + qishs*(one - exp(-bsats*eqps0))
dradpr = bsats*(yield0 + qishs - radpr0)*eqps0 ! differentiation of R plateau region
yindex = zero
stateOld(i,33) = yindex
if (smises - yield0 .ge. zero) then
yindex = one
smises = smises + (one - yindex)
C On the bounding surface
delta = ( yindex*(smises - yield0) )/thremu ! equivalent plastic strain increment
nflow1 = s1 / (yield0 + (thremu*delta)) ! n_i+1
nflow2 = s2 / (yield0 + (thremu*delta))
nflow3 = s3 / (yield0 + (thremu*delta))
nflow4 = s4 / (yield0 + (thremu*delta))
nflow5 = s5 / (yield0 + (thremu*delta))
nflow6 = s6 / (yield0 + (thremu*delta))
ps_inc1 = threeHalfs * nflow1 * delta ! updated plastic strain increment
ps_inc2 = threeHalfs * nflow2 * delta
ps_inc3 = threeHalfs * nflow3 * delta
ps_inc4 = threeHalfs * nflow4 * delta
ps_inc5 = threeHalfs * nflow5 * delta
ps_inc6 = threeHalfs * nflow6 * delta
pls1 = stateOld(i,15) ! plastic strain
pls2 = stateOld(i,16)
pls3 = stateOld(i,17)
pls4 = stateOld(i,18)
pls5 = stateOld(i,19)
pls6 = stateOld(i,20)
stateNew(i,15) = pls1 + ps_inc1
stateNew(i,16) = pls2 + ps_inc2
stateNew(i,17) = pls3 + ps_inc3
stateNew(i,18) = pls4 + ps_inc4
stateNew(i,19) = pls5 + ps_inc5
stateNew(i,20) = pls6 + ps_inc6
ps1 = stateNew(i,15)
ps2 = stateNew(i,16)
ps3 = stateNew(i,17)
ps4 = stateNew(i,18)
ps5 = stateNew(i,19)
ps6 = stateNew(i,20)
ds1 = s1 - (twomu * ps_inc1) ! update deviatoric stress
ds2 = s2 - (twomu * ps_inc2)
ds3 = s3 - (twomu * ps_inc3)
ds4 = s4 - (twomu * ps_inc4)
ds5 = s5 - (twomu * ps_inc5)
ds6 = s6 - (twomu * ps_inc6)
epxi1 = ps1 - stateOld(i,8)
epxi2 = ps2 - stateOld(i,9)
epxi3 = ps3 - stateOld(i,10)
epxi4 = ps4 - stateOld(i,11)
epxi5 = ps5 - stateOld(i,12)
epxi6 = ps6 - stateOld(i,13)
term1 = epxi1**2 + epxi2**2 + epxi3**2 + two*(epxi4**2 + epxi5**2
1 + epxi6**2)
glamda = sqrt(twoThirds*term1) - rad_ms0
if (glamda .gt. zero) then
gdelta = glamda
else
gdelta = zero
end if
radpr1 = ( radpr0 + ( bsats*gdelta*(yield0 + qishs) ) ) ! updating R_i to R_i+1
1 / ( one + (bsats*gdelta) )
factor2 = one - ( radpr1 / yield0 )
trsig = ps_inc1 + ps_inc2 + ps_inc3
snew1 = sig1 - (xkap*trsig) - twomu*(ps_inc1 - (third*trsig)) ! new stress
snew2 = sig2 - (xkap*trsig) - twomu*(ps_inc2 - (third*trsig))
snew3 = sig3 - (xkap*trsig) - twomu*(ps_inc3 - (third*trsig))
snew4 = sig4 - twomu*(ps_inc4)
snew5 = sig5 - twomu*(ps_inc5)
snew6 = sig6 - twomu*(ps_inc6)
stressNew(i,1) = snew1
stressNew(i,2) = snew2
stressNew(i,3) = snew3
stressNew(i,4) = snew4
stressNew(i,5) = snew5
stressNew(i,6) = snew6
end if
C else
end if
end do
return
end
i am having a problem with my vumat subroutine. It is based on a constitutive model for full range elastoplastic behavior. You can see the code below. My problem is that at the end of a simulation, the abaqus displays zero stress (zero for the smises,s11,s22,etc) everything thing after if the "if-loop" becomes zero. But once i delete the plastic part, abaqus does display the new stresses for the elastic part of the code. Do you have an idea for why that is? I really need your help.
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
C The state variables are stored as:
C state(*,1) = sback stress component 11 ! short-range backstress
C state(*,2) = sback stress component 22
C state(*,3) = sback stress component 33
C state(*,4) = sback stress component 12
C state(*,5) = sback stress component 13
C state(*,6) = sback stress component 23
C state(*,7) = equivalent plastic strain
C state(*,8) = center of the memory surface 11
C state(*,9) = center of the memory surface 22
C state(*,10) = center of the memory surface 33
C state(*,11) = center of the memory surface 12
C state(*,12) = center of the memory surface 13
C state(*,13) = center of the memory surface 23
C state(*,14) = radius of the memory surface
C state(*,15) = plastic strain 11
C state(*,16) = plastic strain 22
C state(*,17) = plastic strain 33
C state(*,18) = plastic strain 12
C state(*,19) = plastic strain 13
C state(*,20) = plastic strain 23
C state(*,21) = elastic strain 11
C state(*,22) = elastic strain 22
C state(*,23) = elastic strain 33
C state(*,24) = elastic strain 12
C state(*,25) = elastic strain 13
C state(*,26) = elastic strain 23
C state(*,27) = lback stress component 11 ! long-range backstress
C state(*,28) = lback stress component 22
C state(*,29) = lback stress component 33
C state(*,30) = lback stress component 12
C state(*,31) = lback stress component 13
C state(*,32) = lback stress component 23
C state(*,33) = yield index ! predefined variable
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),
4 defgradOld(nblock,ndir+nshr+nshr),
5 fieldOld(nblock,nfieldv), stressOld(nblock,ndir+nshr),
6 stateOld(nblock,nstatev), enerInternOld(nblock),
7 enerInelasOld(nblock), tempNew(nblock),
8 stretchNew(nblock,ndir+nshr),
9 defgradNew(nblock,ndir+nshr+nshr),
1 fieldNew(nblock,nfieldv),
2 stressNew(nblock,ndir+nshr), stateNew(nblock,nstatev),
3 enerInternNew(nblock), enerInelasNew(nblock)
C
character*80 cmname
C
parameter ( zero = 0.d0, one = 1.d0, two = 2.d0, three = 3.d0,
1 third = one/three, half = .5d0, twoThirds = two/three,
2 threeHalfs = three/two, six = 6.d0)
C
C Elastic constants
e = props(1) ! E-Mod
xnu = props(2) ! Poisson's ratio
yield0 = props(3) ! initial yield stress
epst = props(4) ! plastic strain at the end of the yield plateau in monotonic curve
epst_amp = props(5) ! threshold amplitude of plastic strain
cskin1 = props(6) ! short-range kinematic hardening parameter, c_s
cskin2 = props(7)
clkin1 = props(8) ! long-range kinematic hardening parameter, c_l
clkin2 = props(9)
sgamma1 = props(10) ! short-range kinematic hardening parameter, gamma_s
sgamma2 = props(11)
lgamma1 = props(12) ! long-range kinematic hardening parameter, gamma_l
lgamma2 = props(13)
qishs = props(14) ! short-range saturated value of the isotropic hardening component, Q
qishl = props(15) ! long-range
bsats = props(16) ! short-range rate in which the saturation is achieved, b_s
bsatl = props(17) ! long-range rate, b_l
cs_para = props(18) ! short-range scalar parameter, c_s
cl_para = props(19) ! long-range, c_l
C
mu = e / ( two * ( one + xnu ) ) ! mu (lamé constant) = G (shear modulus)
twomu = e / ( one + xnu ) ! 2G = 2mu
thremu = three * e / ( two * ( one + xnu ) ) ! 3G = 3mu
C alamda = twomu * ( e - twomu ) / ( sixmu - two * e ) ! lamé constant
xkap = e / ( three * ( one - two*xnu ) )
con = xnu / e
coef1 = sqrt(twoThirds)
coef2 = sqrt(threeHalfs)
C Deviatoric strain increment
do i = 1, nblock
trace = strainInc(i,1) + strainInc(i,2) + strainInc(i,3)
de1 = strainInc(i,1) - ( third * trace )
de2 = strainInc(i,2) - ( third * trace )
de3 = strainInc(i,3) - ( third * trace )
de4 = strainInc(i,4)
de5 = strainInc(i,5)
de6 = strainInc(i,6)
es1 = de5
es5 = trace
C Trial elastic Deviatoric stress incremment
s1_inc = twomu * de1
s2_inc = twomu * de2
s3_inc = twomu * de3
s4_inc = twomu * de4
s5_inc = twomu * de5
s6_inc = twomu * de6
es2 = twomu
es3 = s5_inc
C Predictor/Trial stress
sig1 = stressOld(i,1) + ( xkap*trace ) + s1_inc
sig2 = stressOld(i,2) + ( xkap*trace ) + s2_inc
sig3 = stressOld(i,3) + ( xkap*trace ) + s3_inc
sig4 = stressOld(i,4) + s4_inc
sig5 = stressOld(i,5) + s5_inc
sig6 = stressOld(i,6) + s6_inc
str1 = sig5
es4 = xkap
C Trial elastic and plastic strain
elas1 = stateOld(i,21) ! Elastic strain
elas2 = stateOld(i,22)
elas3 = stateOld(i,23)
elas4 = stateOld(i,24)
elas5 = stateOld(i,25)
elas6 = stateOld(i,26)
eltr1 = elas1 + strainInc(i,1) ! Trial elastic strain
eltr2 = elas2 + strainInc(i,2)
eltr3 = elas3 + strainInc(i,3)
eltr4 = elas4 + strainInc(i,4)
eltr5 = elas5 + strainInc(i,5)
eltr6 = elas6 + strainInc(i,6)
pls1 = stateOld(i,15) ! plastic strain
pls2 = stateOld(i,16)
pls3 = stateOld(i,17)
pls4 = stateOld(i,18)
pls5 = stateOld(i,19)
pls6 = stateOld(i,20)
if( stepTime .eq. zero ) then
stressNew(i,1) = sig1
stressNew(i,2) = sig2
stressNew(i,3) = sig3
stressNew(i,4) = sig4
stressNew(i,5) = sig5
stressNew(i,6) = sig6
else
C Trace of stress tensor
smean = third * ( sig1 + sig2 + sig3 )
str2 = smean
C Deviatoric stress from the backstress
s1 = sig1 - smean
s2 = sig2 - smean
s3 = sig3 - smean
s4 = sig4
s5 = sig5
s6 = sig6
C Equivalent plastic strain
eqps0 = stateOld(i,7) ! initial value should be zero
rad_ms0 = stateOld(i,14)
C Von Mises equivalent stress
smises = s1**2 + s2**2 + s3**2 + two*(s4**2 + s5**2 + s6**2)
smises = sqrt(threeHalfs*smises)
C Radius of yield
radpr0 = yield0 + qishs*(one - exp(-bsats*eqps0))
dradpr = bsats*(yield0 + qishs - radpr0)*eqps0 ! differentiation of R plateau region
yindex = zero
stateOld(i,33) = yindex
if (smises - yield0 .ge. zero) then
yindex = one
smises = smises + (one - yindex)
C On the bounding surface
delta = ( yindex*(smises - yield0) )/thremu ! equivalent plastic strain increment
nflow1 = s1 / (yield0 + (thremu*delta)) ! n_i+1
nflow2 = s2 / (yield0 + (thremu*delta))
nflow3 = s3 / (yield0 + (thremu*delta))
nflow4 = s4 / (yield0 + (thremu*delta))
nflow5 = s5 / (yield0 + (thremu*delta))
nflow6 = s6 / (yield0 + (thremu*delta))
ps_inc1 = threeHalfs * nflow1 * delta ! updated plastic strain increment
ps_inc2 = threeHalfs * nflow2 * delta
ps_inc3 = threeHalfs * nflow3 * delta
ps_inc4 = threeHalfs * nflow4 * delta
ps_inc5 = threeHalfs * nflow5 * delta
ps_inc6 = threeHalfs * nflow6 * delta
pls1 = stateOld(i,15) ! plastic strain
pls2 = stateOld(i,16)
pls3 = stateOld(i,17)
pls4 = stateOld(i,18)
pls5 = stateOld(i,19)
pls6 = stateOld(i,20)
stateNew(i,15) = pls1 + ps_inc1
stateNew(i,16) = pls2 + ps_inc2
stateNew(i,17) = pls3 + ps_inc3
stateNew(i,18) = pls4 + ps_inc4
stateNew(i,19) = pls5 + ps_inc5
stateNew(i,20) = pls6 + ps_inc6
ps1 = stateNew(i,15)
ps2 = stateNew(i,16)
ps3 = stateNew(i,17)
ps4 = stateNew(i,18)
ps5 = stateNew(i,19)
ps6 = stateNew(i,20)
ds1 = s1 - (twomu * ps_inc1) ! update deviatoric stress
ds2 = s2 - (twomu * ps_inc2)
ds3 = s3 - (twomu * ps_inc3)
ds4 = s4 - (twomu * ps_inc4)
ds5 = s5 - (twomu * ps_inc5)
ds6 = s6 - (twomu * ps_inc6)
epxi1 = ps1 - stateOld(i,8)
epxi2 = ps2 - stateOld(i,9)
epxi3 = ps3 - stateOld(i,10)
epxi4 = ps4 - stateOld(i,11)
epxi5 = ps5 - stateOld(i,12)
epxi6 = ps6 - stateOld(i,13)
term1 = epxi1**2 + epxi2**2 + epxi3**2 + two*(epxi4**2 + epxi5**2
1 + epxi6**2)
glamda = sqrt(twoThirds*term1) - rad_ms0
if (glamda .gt. zero) then
gdelta = glamda
else
gdelta = zero
end if
radpr1 = ( radpr0 + ( bsats*gdelta*(yield0 + qishs) ) ) ! updating R_i to R_i+1
1 / ( one + (bsats*gdelta) )
factor2 = one - ( radpr1 / yield0 )
trsig = ps_inc1 + ps_inc2 + ps_inc3
snew1 = sig1 - (xkap*trsig) - twomu*(ps_inc1 - (third*trsig)) ! new stress
snew2 = sig2 - (xkap*trsig) - twomu*(ps_inc2 - (third*trsig))
snew3 = sig3 - (xkap*trsig) - twomu*(ps_inc3 - (third*trsig))
snew4 = sig4 - twomu*(ps_inc4)
snew5 = sig5 - twomu*(ps_inc5)
snew6 = sig6 - twomu*(ps_inc6)
stressNew(i,1) = snew1
stressNew(i,2) = snew2
stressNew(i,3) = snew3
stressNew(i,4) = snew4
stressNew(i,5) = snew5
stressNew(i,6) = snew6
end if
C else
end if
end do
return
end