问题描述:假设一圆柱半径为5mm,拉伸50mm,弹性模量E=200e3Mpa,泊松比V=0.3,本构曲线取线性强化,圆柱一端约束Z向变形和周向变形,另一端拉伸0.9mm,试对比ABAQUS自带本构与UMAT子程序的计算结果。
对比结果:abaqus结果见图左,umat子程序结果见图右
通过对比发现,Umat子程序与Abaqus自带本构得出结果值完全一致。
umat源码如下:
SUBROUTINE UMAT(STRESS,STATEV,DDSDDE,SSE,SPD,SCD,
1 RPL,DDSDDT,DRPLDE,DRPLDT,
2 STRAN,DSTRAN,TIME,DTIME,TEMP,DTEMP,PREDEF,DPRED,CMNAME,
3 NDI,NSHR,NTENS,NSTATV,PROPS,NPROPS,COORDS,DROT,PNEWDT,
4 CELENT,DFGRD0,DFGRD1,NOEL,NPT,LAYER,KSPT,JSTEP,KINC)
C
INCLUDE 'ABA_PARAM.INC'
C
CHARACTER*80 CMNAME
DIMENSION STRESS(NTENS),STATEV(NSTATV),
1 DDSDDE(NTENS,NTENS),DDSDDT(NTENS),DRPLDE(NTENS),
2 STRAN(NTENS),DSTRAN(NTENS),TIME(2),PREDEF(1),DPRED(1),
3 PROPS(NPROPS),COORDS(3),DROT(3,3),DFGRD0(3,3),DFGRD1(3,3),
4 JSTEP(4)
double precision E,V,yield0,H,lamdba,G,STATEV
double precision hydroStaticPressure,misesEqualStress
double precision dEqualPlasticStrain
double precision, DIMENSION(NTENS)::deviatoricStress
INTEGER i,j
C Elastic properties
E = PROPS(1)
V = PROPS(2)
yield0 = PROPS(3)
H = (PROPS(5)-PROPS(3))/(PROPS(6)-PROPS(4))
lamdba = E*V/((1+V)*(1-2*V))
G = E/(2*(1+V))
C Calculate the material jacobian matrix
do i=1, NTENS
do j=1, NTENS
DDSDDE(i,j)=0
end do
end do
do i=1,NDI
do j=1,NDI
DDSDDE(i,j)=lamdba
end do
DDSDDE(i,i)=lamdba+2*G
DDSDDE(i+NDI,i+NDI)=G
end do
C Calculate stress-try
do i=1,NTENS
do j=1,NTENS
stress(i)=stress(i)+DDSDDE(i,j)*DSTRAN(j)
end do
end do
C Calculate hydroStaticPressure
hydroStaticPressure=0
do i=1,NDI
hydroStaticPressure=hydroStaticPressure+stress(i)
end do
hydroStaticPressure=hydroStaticPressure/3
C Calculate the deviatoricStress-try
do i=1,NDI
deviatoricStress(i)=STRESS(i)-hydroStaticPressure
end do
do i=1+NDI,NTENS
deviatoricStress(i)=stress(i)
end do
C Calculate misesEqualStress
misesEqualStress=0
do i=1,NDI
misesEqualStress=misesEqualStress+deviatoricStress(i)**2
end do
misesEqualStress=sqrt(1.5*misesEqualStress)
C get the nth equalPlasticStrain from STATEV
equalPlasticStrain=STATEV(1)
C Calculate the nth yield
yieldn=yield0+H*equalPlasticStrain
if (misesEqualStress>yieldn) then
dEqualPlasticStrain=(misesEqualStress-yieldn)/(H+3*G)
do i=1,NTENS
deviatoricStress(i)=deviatoricStress(i)-3*G*dEqualPlasticStrain*deviatoricStress(i)/misesEqualStress
end do
do i=1,NDI
stress(i)=deviatoricStress(i)+hydroStaticPressure
end do
do i=1+NDI,NTENS
stress(i)=deviatoricStress(i)
end do
else
dEqualPlasticStrain=0
end if
STATEV(1)=equalPlasticStrain+dEqualPlasticStrain
RETURN
END
59~64:计算偏试应力: 67~71:计算Mises等效应力: 74:从STATEV状态变量中获取等效塑性应变 77:弹塑性线性强化本构关系, 79:判断是否屈服 80:计算等效塑性应变增量, 82~84:如果屈服,则更新偏试应力, 86~91:更新试应力 93:如果没有屈服,则不会发生塑性流动 96:将等效塑性应变储存在状态变量中 以上是我结合技 术 邻论坛上SuperID大神的模型,进一步讨论对弹塑性线性强化常刚度Umat的理解,下节内容将更新弹塑性线性强化切线刚度Umat的相关内容
来源:易木木响叮当