关于vusdfld中stateOld和stateNew的问题

Tao_Tony

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

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

这个问题困扰我很久了，能否有大神助一臂之力的？
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

Lockie

楼主问题解决了么？？

Nonwhy

楼主解决了吗可以一块商量一下，我也在搞这个方向