損傷Vumat子程序求助
下面是小弟寫的一個損傷vumat子程序（lemaitre本構）。現在在單胞測試中發現，D_c和depl都有合理的數值，但是在一次更新后，損傷量D和塑性應變一直是0。還請各位大神幫忙指點一下！
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),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
C
real*8 n,eqplas,E_0,D_C,xnu,E_D,E_R,a_1,a_2
CHARACTER*80 cmname
PARAMETER( zero=0.,one=1.,two=2.,three=3.,
1 third=one/three,half=.5,twoThirds=two/three,
2 threeHalfs=1.5,zt=.2 )
C
C
do 100 i = 1,nblock
e = props(1)
xnu = props(2)
D_C = props(3)
E_D = props(4)
E_R = props(5)
K = props(6) !強度系數
n = props(7) !硬化指數
E_0 = props(8) !預應變
c eqplas = stateOld(i,7)
c D = stateOld(i,8)
c e = (one-D)*eC 剪切和體積模量
G = e/(two*(one+xnu)) ! G=E/2(1+μ)
G2 = e/(one+xnu)
bulk = e/(one-two*xnu) * third ! K=E/3(1-2μ)
C
C
C 應變增量
trace = strainInc(i,1)+strainInc(i,2)+strainInc(i,3)
C 偏應變增量（deviatoric）
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)
C 預估彈性偏應力增量 柯西應力增量
s1_inc=two*G*de1
s2_inc=two*G*de2
s3_inc=two*G*de3
s4_inc=two*G*de4
s5_inc=two*G*de5
s6_inc=two*G*de6
C 靜壓增量
p_inc= -bulk*trace
C 考慮損傷下的試探應力
sig1 = stressOld(i,1) + (one-D)*(s1_inc - p_inc)
sig2 = stressOld(i,2) + (one-D)*(s2_inc - p_inc)
sig3 = stressOld(i,3) + (one-D)*(s3_inc - p_inc)
sig4 = stressOld(i,4) + (one-D)*s4_inc
sig5 = stressOld(i,5) + (one-D)*s5_inc
sig6 = stressOld(i,6) + (one-D)*s6_inc
print*,"損傷因子D=",D
C Set up old equivalent plastic strain
if(stepTime+totalTime .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 total trial hydrostatic stress
p = -third * (sig1+sig2+sig3)C 考慮平均應力的試探應力
s1 = sig1 + p
s2 = sig2 + p
s3 = sig3 + p
s4 = sig4
s5 = sig5
s6 = sig6
C mises等效應力
mises = s1**2 + s2**2 + s3**2
1 + two*s4**2 + two*s5**2 + two*s6**2
smises = threeHalfs*mises
smises = sqrt(smises)
print*,"米賽斯應力=",smises
C 流動應力和屈服應力
flow = K*(E_0+eqplas)**n
yield = K*(one-D)*(E_0+eqplas)**n
print*,"屈服應力=",yield
print*,"等效塑性應變=",eqplasC 判別屈服條件
if(smises .le. yield) then
stressNew(i,1) = s1
stressNew(i,2) = s2
stressNew(i,3) = s3
stressNew(i,4) = s4
stressNew(i,5) = s5
stressNew(i,6) = s6
else
c 求解塑性因子
U= D_C/(E_R-E_D)*(twoThirds*(one + xnu)
1 + three*(one - two*xnu)*((p/smises)**2))
2 *((flow/K)**2)
print*,"U=",U
H= n*K*(E_0 + eqplas)**(n-one)
print*,"H=",H
print*,"n=",n
A= twoThirds*sqrt(twoThirds)*U*H
print*,"A=",A
B= third*U*yield - G*(one-D)*(one + H/(three*G))
print*,"B=",B
C= smises - sqrt(twoThirds)*(one - D)*flow
print*,"C=",C
a_1 = (-B + sqrt(B**two-A*C))/A
a_2 = (-B - sqrt(B**two-A*C))/A
print*,"a_1=",a_1
print*,"a_2=",a_2
if(a_1>a_2)then
X=a_2
else
X=a_1
end if
print*,"X=",X
C 損傷增量
D_inc = sqrt(twoThirds)*U*X
print*,"D_inc=",D_inc
C 塑性應變增量
depl = sqrt(twoThirds)*X
C 尋找單位張量n
xnormal1 = threeHalfs*s1/smises
xnormal2 = threeHalfs*s2/smises
xnormal3 = threeHalfs*s3/smises
xnormal4 = threeHalfs*s4/smises
xnormal5 = threeHalfs*s5/smises
xnormal6 = threeHalfs*s6/smises
C 更新試探應力
s_trial1 = s_trial1 - two*G*X*xnormal1
s_trial2 = s_trial2 - two*G*X*xnormal2
s_trial3 = s_trial3 - two*G*X*xnormal3
s_trial4 = s_trial4 - two*G*X*xnormal4
s_trial5 = s_trial5 - two*G*X*xnormal5
s_trial6 = s_trial6 - two*G*X*xnormal6
C 更新應力，寫入stressNew中
stressNew(i,1) = s_trial1
stressNew(i,2) = s_trial2
stressNew(i,3) = s_trial3
stressNew(i,4) = s_trial4
stressNew(i,5) = s_trial5
stressNew(i,6) = s_trial6
end if
C 更新狀態變量
eqplas = stateOld(i,7) + depl
D = stateOld(i,8) + D_inc
C 更新狀態變量，寫入stateNew中
stateNew(i,7) = stateOld(i,7) + depl
stateNew(i,8) = stateOld(i,8) + D_inc
C 判斷材料是否開裂，狀態變量14控制單元刪除
if(D .ge. D_C)then
stateNew(i,14) = zero
else
stateNew(i,14) = one
end if
print*,"D_C=",D_CC 更新內能、耗散能
stressPower = half * (
1 ( stressOld(i,1)+stressNew(i,1) )*strainInc(i,1)
2 + ( stressOld(i,2)+stressNew(i,2) )*strainInc(i,2)
3 + ( stressOld(i,3)+stressNew(i,3) )*strainInc(i,3)
4 + two*( stressOld(i,4)+stressNew(i,4) )*strainInc(i,4)
5 + two*( stressOld(i,5)+stressNew(i,5) )*strainInc(i,5)
6 + two*( stressOld(i,6)+stressNew(i,6) )*strainInc(i,6))
C
enerInternNew(i) = enerInternOld(i)
1 + stressPower / density(i)
C Update the dissipated inelastic specific energy
equivStress = threeHalfs*((stressNew(i,1)-smean)**2
1 +(stressNew(i,2)-smean)**2
2 +(stressNew(i,3)-smean)**2
3 +two*stressNew(i,4)**2
4 +two*stressNew(i,5)**2
5 +two*stressNew(i,6)**2)
if(equivStress .ne. zero) then
equivStress = sqrt(equivStress)
end if
plasticWorkInc = equivStress*depl
end if
100 continue
return
end
小程序用戶_qne**3ld_ZZYK 
