huang隐式程序修改为显式及计算案例

黄永刚原始晶体塑性具有良好的收敛性，以及高效的计算效率，在一般变形下无需修改，即可直接使用。然而一些特殊的工况，如切削，轧制，冲压等隐式存在收敛性问题。因此通常使用显示程序进行计算。但从头完成显式晶体塑性构造对于一般学者显然难度过高，一个简单的想法就是直接将现成的黄永刚隐式程序改成显式。abaqus里这是可以实现的。其基本的步骤是：

1，加入vumat接口程序（见附录abaqus官网有）

2，对nblock进行循环，计算应力和状态变量

3，更新应力与状态变量，重复计算直到增量结束。

值得注意的是，umat与vumat程序里面剪应力分量定义顺序与应力不同

umat：12，13，23（工程剪应变）

vumat：12，23，13（2*工程剪应变）

同时采用该方法计算时计算效率显著高于完全显式，并允许较大的时间增量。为评估模型计算效率，采用1000个晶粒80000个单元的二维模型进行20%的压缩模拟。耗时3小时，计算结果与隐式结果类似。模型结果如下：

附录

      subroutine vumat(
C Read only (unmodifiable)variables -
     1  nblock, ndir, nshr, nstatev, nfieldv, nprops, jInfoArray,
     2  stepTime, totalTime, dtArray, 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 (modifiable) variables -
     7  stressNew, stateNew, enerInternNew, enerInelasNew )
C
      include 'vaba_param.inc'
      parameter (i_info_AnnealFlag = 1, 
     *     i_info_Intpt    = 2, ! Integration station number
     *     i_info_layer  = 3, ! Layer number
     *     i_info_kspt   = 4, ! Section point number in current layer
     *     i_info_effModDefn = 5, ! =1 if Bulk/ShearMod need to be defined
     *     i_info_ElemNumStartLoc   = 6) ! Start loc of user element number
C
      dimension props(nprops), density(nblock), coordMp(nblock,*),
     1  charLength(nblock), dtArray(2*(nblock)+1), 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),
     8  defgradNew(nblock,ndir+nshr+nshr),
     9  fieldNew(nblock,nfieldv),
     1  stressNew(nblock,ndir+nshr), stateNew(nblock,nstatev),
     2  enerInternNew(nblock), enerInelasNew(nblock), jInfoArray(*)
C
      character*80 cmname
C
      pointer (ptrjElemNum, jElemNum)
      dimension jElemNum(nblock)
C
      lAnneal = jInfoArray(i_info_AnnealFlag) 
      iLayer = jInfoArray(i_info_layer)
      kspt   = jInfoArray(i_info_kspt)
      intPt  = jInfoArray(i_info_Intpt)
      iUpdateEffMod = jInfoArray(i_info_effModDefn)
      iElemNumStartLoc = jInfoArray(i_info_ElemNumStartLoc)
      ptrjElemNum = loc(jInfoArray(iElemNumStartLoc))

      do 100 km = 1,nblock
        
  100 continue

      return