首页/文章/ 详情

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

1年前浏览929

问题描述:假设一圆柱半径为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/3C     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年前
    易木木响叮当
    硕士 有限元爱好者
    获赞 187粉丝 173文章 282课程 2
    点赞
    收藏
    未登录
    1条评论
    仿了个真真
    秀不起来
    10月前
    老师好,请问子程序里面变量456都没定义,没问题吗?
    回复
    课程
    培训
    服务
    行家
    VIP会员 学习 福利任务 兑换礼品
    下载APP
    联系我们
    帮助与反馈