首页/文章/ 详情

弹塑性线性强化常刚度Umat(1)

1年前浏览1117

问题描述:假设一圆柱半径为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,jC     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


接下来是代码的Step-by-Step解读:
1~14:子程序的框架可以到Aabqus帮助文档里面Copy
15~19:声明语句,其中第18行声明了6维数组
21~24:从abaqus读取材料属性,强化系数H=屈服应力差值/塑性应变
26~27:拉梅系数,  ,  


30~34:初始化雅克比矩阵(十分重要!!!)
36~49:计算初始雅克比矩阵与试应力,








此处说一下试应力的作用:若试应力在屈服面外面,则会产生塑性流动,将试应力拉回屈服面上。示意图如下:




51~56:计算静水压力,计算之前一定要初始化,  

59~64:计算偏试应力:      


67~71:计算Mises等效应力:  


74:从STATEV状态变量中获取等效塑性应变

77:弹塑性线性强化本构关系,  


79:判断是否屈服

80:计算等效塑性应变增量,  


82~84:如果屈服,则更新偏试应力,  


86~91:更新试应力

93:如果没有屈服,则不会发生塑性流动

96:将等效塑性应变储存在状态变量中


以上是我结合技 术 邻论坛上SuperID大神的模型,进一步讨论对弹塑性线性强化常刚度Umat的理解,下节内容将更新弹塑性线性强化切线刚度Umat的相关内容

来源:易木木响叮当

ACTAbaqusUM材料
著作权归作者所有,欢迎分享,未经许可,不得转载
首次发布时间:2023-05-31
最近编辑:1年前
易木木响叮当
硕士 有限元爱好者
获赞 218粉丝 255文章 348课程 2
点赞
收藏
未登录
1条评论
仿了个真真
秀不起来
1年前
老师好,请问子程序里面变量456都没定义,没问题吗?
回复
课程
培训
服务
行家
VIP会员 学习 福利任务 兑换礼品
下载APP
联系我们
帮助与反馈