损伤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