最近，我再將以前寫的一個隱式計算的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