Tresca本构模型VUMAT（2D/轴对称）（附源代码和详细注释）
一、ABAQUS自带Tresca本构与VUMAT对比
二、Tresca本构模型介绍
以下， 粗体符号表示向量或矩阵，上标“T”表示向量转置。
当屈服函数f（σ）的值为零时，材料屈服。这里σ是应力张量（为列矩阵）。如采用相关联的流动法则，则无穷小的塑性应变增量为
式中，a是屈服函数的梯度，dλ是一个非负塑性乘子（non-negative plastic multiplier）。
如果在分析增量步j 和 j+1之间施加一个应变增量Δε，则应力根据下式进行更新
式中，εj 是增量步 j 之后的总应变，Dep是无穷小弹塑性本构矩阵。
由于Dep高度依赖于σ，因此上式（4.2）通常取近似解。应变增量dε由弹性增量dεe和塑性增量dεp组成，因此：
根据胡克定律得到弹性应力增量如下：
式中D是弹性本构矩阵。
结合式（4.1）与式（4.4）可得：
综上可得下图中的径向映射法则（radial return mapping）：
上图中，在增量开始时，考虑点处的应力为σA，其在，弹性区域（f<0，f 是屈服函数），当然其同样可以位于屈服面上（f=0）。弹性预测增量为Δσe，由此预测应力为σB。通过塑性修正增量−Δσp将应力返回到屈服面上的最终应力点σC。其中Δσp由下式得出：
通常，Δσp由下式进行数值计算：
其中，下标P表示积分路径上要计算a的点。在后向欧拉算法中，由于P对应于未知的最终应力状态σC，因此需要迭代。对于主应力呈线性的Tresca屈服准则，a沿积分路径是常数。
以下介绍Tresca屈服准则在主应力空间中的表示。屈服面采用与线σ1=σ2=σ3对齐的六角棱镜的形式，可由六个屈服函数定义：
对于i=1,…,6，方程fi（σ，k）=0对应于主应力空间中的平面Si ,这六个屈服面如下图所示。
三、Tresca VUMAT 算法实现
输入：
输出：
1：
2: 转化至主应力空间
3:
4:
6：
四、Tresca VUMAT 源代码（带注释）
C ABAQUS 自带变量
subroutine vumat(
C 只读-
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 只写-
7 stressNew,stateNew,enerInternNew,enerInelasNew)
C
include 'vaba_param.inc'
dimension props(nprops),density(nblock),
1 coordMp(nblock,*),
2 charLength(*),strainInc(nblock,ndir+nshr),
3 relSpinInc(*),tempOld(*),
4 stretchOld(*),defgradOld(*),
5 fieldOld(*),stressOld(nblock,ndir+nshr),
6 stateOld(nblock,nstatev),enerInternOld(nblock),
7 enerInelasOld(nblock),tempNew(*),
8 stretchNew(*),defgradNew(*),fieldNew(*),
9 stressNew(nblock,ndir+nshr),stateNew(nblock,nstatev),
1 enerInternNew(nblock),enerInelasNew(nblock)
C
character*80 cmname
* !
* ! ************************** 用户自定义变量 ************* *************
* !
integer :: block_count
dimension DElas (4 ,4) , test(4,4),
1 stressNewInt (4 ,1) ,
2 strainIncInt (4 ,1) ,
3 strPr (3 ,1) ,
4 a1 (3 ,1) , a2 (3 ,2) ,
5 dStressPlas (3 ,1),
6 stressinc(4,1)
* ! ************************************************** ****************
C ! ** 检查ndir和nshr是否符合2D或axisymmetry（ndir =3 and nshr =1）
if (( ndir.ne.3) . or . ( nshr . ne . 1)) then
call xplb_exit
end if
* !
C ! ** 构建弹性刚度矩阵
* !
C PROPS(1) - E ,为弹性模量
C PROPS(2) - v ,为泊松比
C PROPS(3) - c, 极限剪应力
E = props (1)
v = props (2)
c=props(3)
zero=0.0D0
second=0.5D0
One=1.0D0
* !
tempOne = 1.0 D0 - v
tempTwo = 1.0 D0 - 2.0 D0 *v
* !
DElas (1 ,1:4) = (/ tempOne , v , v , zero /)
DElas (2 ,1:4) = (/ v , tempOne , v , zero /)
DElas (3 ,1:4) = (/ v , v , tempOne , zero /)
DElas (4 ,1:4) = (/ zero , zero , zero ,second * tempTwo/)
* !
DElas = ( E /((1.0 D0 + v )*( tempTwo )) ) * DElas
* !
C ! ************ 开始主循环 ************
* !
do 100 block_count = 1, nblock
* !
C ! 使用工程剪应变
strainIncInt ( 1:3 , 1 ) = strainInc( block_count,1:3 )
strainIncInt ( 4 , 1 ) = 2.0D0*strainInc( block_count,4 )
* ! Elastic predictor stress 弹性预测应力
stressinc=matmul(DElas,strainIncInt)
stressNewInt(1:4,1)=stressOld(block_count,1:4)+stressinc(1:4,1)
* !
* ! Principal values of elastic trial stress 弹性预测主应力
stressMean = 0.5 D0 * ( stressNewInt (1 ,1) +
1 stressNewInt (2 , 1) )
RMohr = sqrt ((0.5 D0 *( stressNewInt (1 ,1) -
1 stressNewInt (2 ,1)))**2
2 + ( stressNewInt (4 ,1))**2)
* !
strPr (1 ,1) = stressMean + RMohr
strPr (2 ,1) = stressMean - RMohr
strPr (3 ,1) = stressNewInt (3 ,1)
twoTheta = atan2 ( stressNewInt (4 ,1) , 0.5 D0
1 *( stressNewIn t (1 ,1) - stressNewInt (2 ,1)) )
* !
* ! Evaluate yield functions 计算屈服函数
tempOne = 2.0 D0 * c
f1 = strPr (1 ,1) - strPr (2 ,1) - tempOne
f2 = strPr (2 ,1) - strPr (1 ,1) - tempOne
f3 = strPr (2 ,1) - strPr (3 ,1) - tempOne
f4 = strPr (3 ,1) - strPr (2 ,1) - tempOne
f5 = strPr (3 ,1) - strPr (1 ,1) - tempOne
f6 = strPr (1 ,1) - strPr (3 ,1) - tempOne
* !
* ! 根据屈服函数值进行修正
if ( ( max ( f1 ,f2 , f3 , f4 ,f5 , f6 ) <= 0.0 D0 ) . OR .
1 ( totalTime == 0.0 D0 )) then
* ! ** ELASTIC 弹性
* ! 无需修正
* !
else if ( ( f1 >= 2.0 D0 * f4 ).AND.( f1 >= 2.0 D0 * f6 ) ) then
* ! ** Return to single surface (1) 修正至面（1）
a1 (1:3 ,1) = (/ -0.5 D0 , 0.5 D0 , 0.0 D0 /)
strPr = strPr + a1 * f1
* !
else if ( ( f4 >= 2.0 D0 * f1 ).AND.( f4 >= 2.0 D0 * f5 ) ) then
* ! ** Return to single surface (4) 修正至面（4）
a1 (1:3 ,1) = (/ 0.0 D0 , 0.5 D0 , -0.5 D0 /)
strPr = strPr + a1 * f4
* !
else if ( ( f6 >= 2.0 D0 * f1 ).AND.( f6 >= 2.0 D0 * f3 ) ) then
* ! ** Return to single surface (6) 修正至面（6）
a1 (1:3 ,1) = (/ -0.5 D0 , 0.0 D0 , 0.5 D0 /)
strPr = strPr + a1 * f6
* !
else if ( f6 < 2.0 D0 * f3 ) then
* ! ** Return to line (6/3) 修正至线（6/3）
a2 (1:3 ,1) =(/ 1.0 D0 /3.0 D0 ,-2.0 D0 /3.0D0,1.0 D0 /3.0 D0 /)
a2 (1:3 ,2) = (/-2.0 D0 /3.0 D0,1.0 D0 /3.0 D0,1.0 D0 /3.0 D0 /)
strPr = strPr +matmul ( a2 , reshape ( (/ f3 , f6 /) ,(/2,1 /)))
* !
else if (( f1 < 2.0 D0 * f6 ) . AND . ( f6 < 2.0 D0 * f1)) then
* ! ** Return to line (1/6) 修正至线（1/6）
a2 (1:3 ,1)=(/-1.0 D0 /3.0 D0,2.0 D0/3.0 D0 ,-1.0 D0 /3.0 D0/)
a2 (1:3 ,2)=(/ -1.0 D0/3.0 D0,-1.0 D0/3.0 D0 ,2.0 D0 /3.0 D0 /)
strPr = strPr + matmul ( a2 ,reshape ( (/ f1, f6/) ,(/ 2,1/) ) )
* !
else if ( ( f1 < 2.0 D0 * f4 ) . AND . ( f4 < 2.0 D0 * f1)) then
* ! ** Return to line (4/1) 修正至线（4/1）
a2 (1:3 ,1) = (/-2.0 D0 /3.0 D0 ,1.0 D0 /3.0 D0,1.0 D0 /3.0D0 /)
a2 (1:3 ,2) = (/ 1.0 D0 /3.0 D0 , 1.0 D0 /3.0 D0,-2.0D0/3.0D0 /)
strPr = strPr + matmul ( a2 , reshape ( (/ f1 ,f4/) ,(/2,1/)))
* !
else if ( f4 < 2.0 D0 * f5 ) then
* ! ** Return to line (5/4) 修正至线（5/4）
a2 (1:3 ,1) = (/ -1.0 D0 /3.0D0,2.0 D0/3.0 D0 ,1.0 D0 /3.0 D0 /)
a2 (1:3 ,2) = (/2.0 D0 /3.0D0 , -1.0D0 /3.0D0,-1.0 D0 /3.0 D0 /)
strPr = strPr + matmul(a2,reshape ((/ f4 , f5 /) ,(/ 2, 1 /) ) )
* !
else
* ! ** 错误
call xplb_exit
end if
* !
* ! ** Transform principal stresses back to Cartesian 将主应力转换回笛卡尔坐标
twoTheta = - twoTheta
stressMean = 0.5 D0 *( strPr (1 ,1) + strPr (2 ,1))
offset = 0.5 D0 *( strPr (1 ,1) - strPr (2 ,1)) * cos (twoTheta)
stressNewInt (1 ,1) = stressMean + offset
stressNewInt (2 ,1) = stressMean - offset
stressNewInt (3 ,1) = strPr (3 ,1)
stressNewInt (4 ,1) = 0.5 D0 *( strPr (2 ,1) - strPr (1 ,1)) *
1 sin( twoTheta )
* !
* ! ** Prepare output 输出，更新应力
stressNew ( block_count , 1:4 ) = stressNewInt ( 1:4 , 1 )
* !
100 continue
end
程序下载：
