最近，我再将以前写的一个隐式计算的USDFLS的子程序改写为VUSDFLD，结果遇到的问题

我使用vusdfld在实现单元的删除，在里面计算等效塑性应变的增量的时候，我令一个变量等于stateOld，然后用stateNew减去这个值来实现我的目的，为什么不行？也就是stateOld里面并不是存的上一步的等效塑性应变的结果？

这个问题困扰我很久了，能否有大神助一臂之力的？

附上子程序和inp文件，以及odb的结果

c
c User subroutine VUSDFLD for user-defined fields
c
      subroutine vusdfld(
c Read only -
     *   nblock, nstatev, nfieldv, nprops, ndir, nshr, 
     *   jElemUid, kIntPt, kLayer, kSecPt, 
     *   stepTime, totalTime, dt, cmname, 
     *   coordMp, direct, T, charLength, props, 
     *   stateOld, 
c Write only -
     *   stateNew, field )
c
      include 'vaba_param.inc'
c
      dimension props(nprops),
     *          jElemUid(nblock), coordMp(nblock, *), 
     *          direct(nblock, 3, 3), T(nblock,3,3), 
     *          charLength(nblock),
     *          stateOld(nblock, nstatev), 
     *          stateNew(nblock, nstatev),
     *          field(nblock, nfieldv)
      character*80 cmname
      parameter( nrData=6 )

      character*3 cData(maxblk*6)
      dimension jData(maxblk*nrData)
      dimension eqps(nblock,maxblk*6),stress(nblock,maxblk*6)
    	real dc,rc,plim,alpha,beta,eqpsOld,damageOld
c
      parameter ( zero = 0.d0 )
	  
	  
C 定义损伤的极限值
      parameter ( one = 1.d0 )
C 通过vgetvrm函数来获取等效塑性应变和应力分量
      jStatus = 1
      call vgetvrm( 'PEEQ', eqps, jData, cData, jStatus )      
	  call vgetvrm( 'S', stress, jData, cData, jStatus )      
C 循环进行进行求解响应的变量	  
      do k = 1, nblock
		 damageOld = stateOld(k,10)
c-----------------------------------------------------------------------c
C 主要是计算应力三轴度	  
C 首相将提取的应力分量赋给状态变量/变量，用来计算mises应力和静水压力	  
C		 stateNew(k,51) = stress(k,1)
		 stress1 = stress(k,1)
C		 stateNew(k,52) = stress(k,2)
		 stress2 = stress(k,2)
C		 stateNew(k,53) = stress(k,3)
		 stress3 = stress(k,3)
C		 stateNew(k,54) = stress(k,4)
		 stress4 = stress(k,4)
C		 stateNew(k,55) = stress(k,5)
		 stress5 = stress(k,5)
C		 stateNew(k,56) = stress(k,6)
		 stress6 = stress(k,6)
C 计算mises应力和静水压力		 
		 smises1 = (stress1-stress2)**2+(stress2-stress3)**2+(stress1-stress3)**2
		 smises = (0.5*smises1)**0.5		 
		 hpress = (stress1+stress2+stress3)/3.0
C 将计算的mises应力和静水压力赋给状态变量		 
		 stateNew(k,1) = smises
		 stateNew(k,2) = hpress
C 计算应力三轴度，并将其赋给状态变量
		 if(smises.le.stoler) then
			triax=0.d0
		 else
			triax=(1.0*hpress)/smises
		 endif
		 stateNew(k,3) = triax
c-----------------------------------------------------------------------c
C 主要是计算断裂应变
c	根据确定的系数和计算的应力三轴度，计算临界断裂应变，并将其赋给状态变量
		 frac_strain = 1.556*exp(-1.5*triax)+0.093
		 stateNew(k,4) = frac_strain
c-----------------------------------------------------------------------c
C 主要是计算损伤增量
C 上一步的等效塑性应变stateNew(k,5)，该步的等效塑性应变stateNew(k,6)，增量stateNew(k,7)		 
C		 if(stepTime.eq.zero) then
C			stateNew(k,5) = 0.0
C			stateNew(k,6) = 0.0
C		 else
			stateNew(k,5) = stateOld(k,6)
C		 endif
		 stateNew(k,6) = eqps(k,1)
		 stateNew(k,7) = stateNew(k,6)-stateNew(k,5)
C 该步的损伤增量		 
		 stateNew(k,8) = stateNew(k,7)/stateNew(k,4)
c-----------------------------------------------------------------------c
C 初始的损伤量		 
C		 if(stepTime.eq.zero) then
C			stateNew(k,9) = 0.0
C			stateNew(k,10) = 0.0
C		 else
			stateNew(k,9) = damageOld
C		 endif
C 该步的损伤增量		 
		 stateNew(k,10) = stateNew(k,9) + stateNew(k,8)
		 stateOld(k,10) = stateNew(k,10)
c-----------------------------------------------------------------------c
C 根据损伤量的大小来对网格进行删除		 
         if(stateNew(k,10).gt.one) then 
c Element Deletion
           stateNew(k,11) = zero
         endif
      end do      
c
      return
      end