abaqus 线弹性各向同性umat开发（附代码）

概述

  帖子介绍了abaqus线弹性各向同性umat的开发过程，目的是科普umat二次开发最基础的工作流程，不涉及复杂高深的推导本构工作，仅仅包括umat子程序接口关键参数介绍、切线刚度矩阵和应力更新计算。并且设计一个立方体表面受动荷载算例，统计了位移、应力等数据，均与abaqus自带的材料模型保持一致。最后给出了子程序计算代码，拿过来就能直接用的那种。

如何调用umat子程序

  调用umat子程序的步骤比较简单，在abaqus的cae界面中就可以实现，以下是具体步骤。  打开箭头1的Edit Material界面，在箭头2的General中点击3的User Material，然后下方会出现4的Mechanical Constants，我们需要在这个地方输入自己用到的参数，在这个帖子的线弹性各向同性材料中，第一个参数弹性模量    为    ，第二个参数为泊松比    为    。
  上面的操作在inp文件中表现为

*Material, name=Material-1
**密度不会传进umat子程序，但是动力计算需要
*Density
2000.,
**User Material关键字的数据会传进umat中
*User Material, constants=2
 1e+10, 0.25

abaqus umat接口介绍

  umat子程序是abaqus平台提供的用于实现自己本构模型的接口，官方文档上可以找到该子程序的介绍，在介绍的最开始，官方文档先整了个活，给出了一个“Warning”

  Warning: The use of this subroutine generally requires considerable expertise. You are cautioned that the implementation of any realistic constitutive model requires extensive development and testing. Initial testing on a single-element model with prescribed traction loading is strongly recommended.

  我不得不说，这段话起到了开门见山、首尾呼应、吸引读者阅读兴趣、画龙点睛、blabla略略略等的作用，但事实上我第一次看这个子程序的时候并没有注意到这段话，后来才发现这是官方文档在给我们打预防针：小崽子们小心点儿！这个子程序很难！
  但其实大家完全不用担心，umat子程序本身并不难，总共就那么几个参数，但是数学难！难于上青天！推本构是一项十分有挑战性的工作，谁搞谁知道，劝退umat二次开发的不是umat接口本身，而是推公式。
  好了，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)

      user coding to define DDSDDE, STRESS, STATEV, SSE, SPD, SCD
      and, if necessary, RPL, DDSDDT, DRPLDE, DRPLDT, PNEWDT

      RETURN
      END

  为了让记忆更加的有逻辑，官方文档大致将上面的参数分为了两大类：需要用户更新的参数和abaqus传进子程序的参数。
  在线弹性各向同性材料中，需要更新的参数有应力项STRESS和切线刚度矩阵DDSDDE。
  三维问题的应力项STRESS可表示为

  注意这里的应力排序，这是abaqus规定好的，一定要按着这个顺序写程序，不然会报错。我在用umat给用户自定义单元可视化中遇到过这个问题，当时我的应变-位移算子顺序跟abaqus规定的不一样，老是算不准，后来搞了很长时间才校正回来。
  三维问题的切线刚度矩阵DDSDDE表达方式不唯一，这里给出两种，后面的代码中分别将两种表达方式封装成了两个函数，可以直接调用。第一种表达方式为

  其中

  第二种表达方式为

  其中

  abaqus传进子程序的参数比较多，这里重点讲解线弹性各项同性材料中涉及的参数。
  参数STRAN(NTENS)为当前增量步开始时刻的应变项，括号中NTENS为应变项的个数，三维问题中为6，与应力项对应。
  参数DSTRAN(NTENS)为当前增量步结束时刻的应变增量，是abaqus估算的一个量，我们可以用它来计算当前增量步的应力增量。
  参数NDI是当前材料的正应力或正应变的个数，对应的是切应力或者切应变的个数NSHR，这两个参数的和为NTENS，即

  参数PROPS(NPROPS)是材料属性数组，是我们在cae界面中定义的参数，NPROPS是定义的参数个数，在这个帖子中，PROPS（1）是杨氏模量，PROPS（2）是泊松比。

算例

  模型尺寸为    ，杨氏模量为    ，泊松比为    ，采用C3D8单元离散，单元总数为125，节点总数为216，网格图为  模型一端设置固定边界，另一端施加竖直向面力，边界条件和荷载示意图为  荷载采用正弦荷载，计算总时长为10s，增量步长设置为0.1s，采用dynamic implicit分析步，施加的荷载幅值为1e5，幅值曲线为  统计了3.1s竖直向位移云图对比为（左边为abaqus自带材料模型计算结果，右边为umat计算结果）  统计了3.1s的    应力云图对比为  统计了3.1s的    应力云图对比为  统计了3.1s的    应力云图对比为  统计了7.9s竖直向位移云图对比为  统计了7.9s的    应力云图对比为  统计了7.9s的    应力云图对比为  统计了7.9s的    应力云图对比为  统计了加载面节点加载向位移时程曲线对比为  统计了加载面节点    时程曲线对比为  统计了加载面节点    时程曲线对比为  统计了加载面节点    时程曲线对比为  可以看到，以上结果吻合得跟造数据一样好！其实就算我复 制一下ODB文件，然后再比较，你们也不知道，这可怎么是好！
  很简单，我把代码扔出来就好了！

代码

  inp文件太长了，我就不放在帖子里面了，只把umat代码放出来，如果需要inp文件的可以私聊我，我发给你，但是验证出来数据有问题别吭气儿，人艰不拆，小声儿私信我，我给你解决。

      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,DFGRDO,DFGRD1,NOEL,NPT,LAYER,KSPT,JSTEP,KINC)
     
      INCLUDE 'ABA_PARAM.INC'
      
      CHARACTER*80CMNAME
      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),DFGRD(3,3),DFGRD1(3,3),
     4 JSTEP(4)
     
     
      young=props(1)
      poiss=props(2)    

      DDSDDE=0.D0     
!这里调用的是切线刚度矩阵的第一种写法，第二种也可以
      call elas1(ddsdde,young,poiss)   
            
      STRESS=STRESS+MATMUL(DDSDDE,DSTRAN)
                       
      RETURN
      END
      
!切线刚度矩阵的第一种写法   

      subroutine elas1(ddsdde,young,poiss) 
      
      INCLUDE 'ABA_PARAM.INC'    
      
      double precision ddsdde(6,6)
      double precision young,poiss
      double precision ar   
      
      ar=young/((1+poiss)*(1-2*poiss))
      
      do i=1,3
       do ii=1,3
         DDSDDE(i,ii)=poiss*ar
       enddo
       DDSDDE(i,i)=(1-poiss)*ar       
      enddo
      
      do i=4,6
       DDSDDE(I,I)=((1-2.0D0*poiss)/2.0D0)*ar
      enddo   
      
      return
      end
      
!切线刚度矩阵的第二种写法            
      subroutine elas2(ddsdde,young,poiss) 
      
      INCLUDE 'ABA_PARAM.INC'
      
      DOUBLE PRECISION DDSDDE(6,6)
      DOUBLE PRECISION young,poiss
      DOUBLE PRECISION KK,GG  
        
      kk=young/(3.0D0*(1-2.0D0*poiss))
      gg=young/(2.0D0*(1+poiss))   
      
      do i=1,3
       do ii=1,3
         DDSDDE(i,ii)=kk-(2.0D0/3.0D0)*gg
       enddo
       DDSDDE(i,i)=kk+(4.0D0/3.0D0)*gg       
      enddo
      
      do i=4,6
       DDSDDE(i,i)=gg
      enddo         
      
      return
      end

来源：有限元先生