alamda = twomu * ( e - twomu ) / ( sixmu - two * e )
term = one / ( twomu * ( one + hard/thremu ) )
con1 = sqrt( twoThirds )
C
do 100 i = 1,nblock
C
C Trial stress
SUBROUTINE VUMAT(
C Read only -
1 nblock, ndir, nshr, nstatev, nfieldv, nprops, lanneal,
2 stepTime, totalTime, dt, cmname, coordMp, charLength,
3 props, density, strainInc, relSpinInc,
1 enerInternNew(nblock), enerInelasNew(nblock)
C
character*80 cmname
C
parameter( zero = 0., one = 1., two = 2., three = 3.,
1 third = one/three, half = .5, twoThirds = two/three,
6
stateNew(i,4) = stateOld(i,4) + factor * s4
C
C Update the stress
factor = twomu * dgamma
stressNew(i,1) = sig1 - factor * ds1
stressNew(i,2) = sig2 - factor * ds2
C
STATE(*,4) = back stress component 12
4
C
STATE(*,5) = equivalent plastic strain
C
C
C All arrays dimensioned by (*) are not used in this algorithm
dimension props(nprops), density(nblock),
Author: xueweek@
1 材料本构模型：
在 Property 中定义材料时，在 General 下选中 User Material，输入 210000、0.3、200、 10000，以上两个数值代表 E、ν，、fy、E’。在用户子程序中代表着 PROPS(1)、PROPS(2)、 PROPS(3)、PROPS(4)。然后在 General 下选中 Depvar，由于该例子中使用了五个状态变量 （背应力张量和累积塑性应变变量），因此在第一项中输入大于 5 的数值即可。另外还需要 输入密度。
C
C The state variables are stored as:
C
STATE(*,1) = back stress component 11
C
STATE(*,2) = back stress component 22
C
STATE(*,3) = back stress component 33
1
建模大家都会，故省略
2 ABAQUS 中 STEP 的设置
由于 VUMAT 需要用到 Explicit 求解，因此需要在 step 步骤中设置 explicit 选项，如下图， 其设置可以用默认设置。
2
3 ABAQUS 调用 VUMAT 用户子程序
同 UMAT 用户子程序的调用方法。 在 Job Manager 中点击 Edit 选项，在 General 选项的最后一项中选择自己建立好的用户 子程序文件。（注：用户子程序文件可以使用文本编辑器进行编辑，当然也可以用 Fortran 编译器，如果对自己的用户子程序文件的语法不放心，可以先用 Fortan 编译器进行编译， 不过编译前要先建立 Project，关于 Fortran 编译，这里不再介绍）。 完成后，点击 submit 即可进行分析。
上一次发过《ABAQUS 初学者用户子程序小例子》，給学习 UMAT 的初学者带来了一定的帮助。现在用到 VUMAT，发现网上 这种小例子很少，关于 VUMAT 的资料也不多。摸索了一天，做个 VUMAT 的小例子供大家分享。
实例：简单的平面平板拉伸，材料本构模型采用随动强化模型， E＝210000MPa, ν=0.3，fy=200MPa, E’=10000MPa。左端约束， 右端施加位移载荷 V＝0.2mm。
4 结果 以下两张图分布是用户子程序和 ABAQUS 自带的材料模型（Standard 求
解器）得到的应力云图，可以看出两种图形基本相同。
3
5 VUMAT 子程序
对于初学者来说，需要注意的是，FORTRAN 对于程序语言格式上的要求。例如，对于 FORTRAN 语言，前六个字符必须空出来，等等。检查语法最好的方法就是在 FORTRAN 编 译器上进行编译。
4 tempOld, stretchOld, defgradOld, fieldOld,
3 stressOld, stateOld, enerInternOld, enerInelasOld,
6 tempNew, stretchNew, defgradNew, fieldNew,
C Write only -
5 stressNew, stateNew, enerInternNew, enerInelasNew )
C
include 'vaba_param.inc'
C
C J2 Mises Plasticity with kinematic hardening for plane
C strቤተ መጻሕፍቲ ባይዱin case.
C Elastic predictor, radial corrector algorithm.
1 coordMp(nblock,*),
2 charLength(*), strainInc(nblock,ndir+nshr),
3 relSpinInc(*), tempOld(*),
4 stretchOld(*), defgradOld(*),
5 fieldOld(*), stressOld(nblock,ndir+nshr),
1 + ( stressOld(i,2)+stressNew(i,2) )*strainInc(i,2)
1 + ( stressOld(i,3)+stressNew(i,3) )*strainInc(i,3)
1 + two*( stressOld(i,4)+stressNew(i,4) )*strainInc(i,4) )
trace = strainInc(i,1) + strainInc(i,2) + strainInc(i,3)
sig1 = stressOld(i,1) + alamda*trace + twomu*strainInc(i,1)
sig2 = stressOld(i,2) + alamda*trace + twomu*strainInc(i,2)
stressNew(i,3) = sig3 - factor * ds3
stressNew(i,4) = sig4 - factor * s4
C
C Update the specific internal energy -
stressPower = half * (
1 ( stressOld(i,1)+stressNew(i,1) )*strainInc(i,1)
2 threeHalfs = 1.5 )
C
e = props(1)
xnu = props(2)
yield = props(3)
hard = props(4)
C
twomu = e / ( one + xnu )
thremu = threeHalfs * twomu
sixmu = three * twomu
sig3 = stressOld(i,3) + alamda*trace + twomu*strainInc(i,3)
sig4 = stressOld(i,4)
+ twomu*strainInc(i,4)
C
C Trial stress measured from the back stress
5
s1 = sig1 - stateOld(i,1) s2 = sig2 - stateOld(i,2) s3 = sig3 - stateOld(i,3) s4 = sig4 - stateOld(i,4) C C Deviatoric part of trial stress measured from the back stress smean = third * ( s1 + s2 + s3 ) ds1 = s1 - smean ds2 = s2 - smean ds3 = s3 - smean C C Magnitude of the deviatoric trial stress difference dsmag = sqrt( ds1**2 + ds2**2 + ds3**2 + 2.*s4**2 ) C C Check for yield by determining the factor for plasticity, C zero for elastic, one for yield radius = con1 * yield facyld = zero if( dsmag - radius .ge. zero ) facyld = one C C Add a protective addition factor to prevent a divide by zero C when dsmag is zero. If dsmag is zero, we will not have exceeded C the yield stress and facyld will be zero. dsmag = dsmag + ( one - facyld ) C C Calculated increment in gamma (this explicitly includes the C time step) diff = dsmag - radius dgamma = facyld * term * diff C C Update equivalent plastic strain deqps = con1 * dgamma stateNew(i,5) = stateOld(i,5) + deqps C C Divide dgamma by dsmag so that the deviatoric stresses are C explicitly converted to tensors of unit magnitude in the C following calculations dgamma = dgamma / dsmag C C Update back stress factor = hard * dgamma * twoThirds stateNew(i,1) = stateOld(i,1) + factor * ds1 stateNew(i,2) = stateOld(i,2) + factor * ds2 stateNew(i,3) = stateOld(i,3) + factor * ds3