【VUMAT第1课】线弹性VUMAT

前言

我们在上篇-子程序开发基础、中篇-各向同性线弹性本构和下篇-各向同性线弹性UMAT实践文章中已经探讨了线弹性本构的UMAT实现。UMAT适用于ABAQUS隐式分析（如静力通用分析步），若要使用ABAQUS显式求解器，则用户自定义的本构需要采用VUMAT子程序实现。

线弹性本构

首先，我们先回顾下线弹性本构及其离散形式，这部分内容其实我们在【UMAT1】中篇-各向同性线弹性本构其实已经做过详细介绍，这里不再赘述，只给出结论：

VUMAT编程实现

基于前述理论，编写的各向同性线弹性VUMAT代码如下：

      subroutine vumat(
C Read only (unmodifiable)variables -
     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 Writeonly (modifiable) variables -
     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),
     8  defgradNew(nblock,ndir+nshr+nshr),
     9  fieldNew(nblock,nfieldv),
     1  stressNew(nblock,ndir+nshr), stateNew(nblock,nstatev),
     2  enerInternNew(nblock), enerInelasNew(nblock)
C
      character*80 cmname
C
      parameter (half=0.5d0, one=1.d0, two=2.d0)
C     elastic properties
      amod=props(1)
      anu=props(2)
      g=amod/(two*(one+anu))
      alam=(anu*amod)/((one+anu)*(one-two*anu))
      twog=two*g
C
      do100 km = 1,nblock
C ----------start: update stress-------------------------------
C           trace of strain increment
            dstrainTrace=strainInc(km,1)+strainInc(km,2)+
     1                   strainInc(km,3)
            stressNew(km,1)=twog*strainInc(km,1)+alam*dstrainTrace+
     1                      stressOld(km,1)
            stressNew(km,2)=twog*strainInc(km,2)+alam*dstrainTrace+
     1                      stressOld(km,2)
            stressNew(km,3)=twog*strainInc(km,3)+alam*dstrainTrace+
     1                      stressOld(km,3)
C
            stressNew(km,4)=twog*strainInc(km,4)+stressOld(km,4)
            stressNew(km,5)=twog*strainInc(km,5)+stressOld(km,5)
            stressNew(km,6)=twog*strainInc(km,6)+stressOld(km,6)
C ----------end: update stress-------------------------------------
C
C ----------start: update elastic internal energy------------------
            stressPower=zero
            do k1=1, ndir
                  stressPower=stressPower+
     1            (stressOld(km,k1)+stressNew(km,k1))*half*
     2            strainInc(km,k1)
            enddo
C
            do k1=ndir+1, ndir+nshr
                  stressPower=stressPower+
     1            (stressOld(km,k1)+stressNew(km,k1))*
     2            strainInc(km,k1)
            enddo
            enerInternNew(km)=enerInternOld(km)+stressPower/density(km)
100   continue
      return
      end

仿真结果

下图左侧为使用Abaqus中的内置模块的仿真结果，右侧为使用二次开发的子程序仿真结果。两者结果一致。