随动硬化von Mises率无关弹塑性本构理论以及umat源代码

1 本构理论

1.1 率形式

本构方程为：

单轴拉伸的应力应变的硬化曲线如下：

根据单轴试验得到硬化部分的曲线：

当仅考虑随动硬化时，屈服面的中心在移动，而屈服面的大小不发生改变，即为常数。

屈服条件为：

增加了背应力来表示屈服面中心移动即随动硬化的效果。式中相对应力的表达式为：

流动法则为：

背应力的演化法则为：

式中：

上标k代表随动部分(kinematic)，表示以下随动硬化曲线的梯度：

1.2 Return mapping算法应力更新

在给定增量应变以及上一步的状态变量的情况下，首先计算试验状态下的量：

计算试验状态下的屈服函数值，判断是否发生屈服。当试验屈服函数值大于0时，说明需要进行塑性更正，反正则说明试验状态即为真实状态。以下对塑性更正环节进行详细说明。

当发生塑性流动时，需要求解以下非线性方程组：

可以将上述非线性方程组简化值一个非线性方程。首先将第一个关于应变的方程代入本构方程中，可得偏应力的表达式为：

减去方程组中的第二式可以得到相对应力的表达式为：

可以看出试验相对应力与真实相对应力之间满足关系式：

那么可以将相对应力表示为：

式中：

最后将该相对应力代入到方程组中最后一个一致性方程中可得：

该非线性方程的未知量为塑性乘子增量，可以用牛顿迭代法进行求解。另外可以利用公式：

代入替换可得背应力和相对应力的公式：

以及最后的非线性方程为：

使用牛顿迭代法求解的公式为：

式中雅克比为：

式中：

是各向同性硬化曲线的梯度。当仅考虑随动硬化时，该项为0。

1.3 一致性切线刚度模量

当没有发现塑性流动时，一致性切线刚度矩阵即为弹性矩阵；当发生塑性流动时，应力更新公式为：

求导可得：

式中：

最终可得一致性切线刚度矩阵表达式为：

计算该一致性切线刚度矩阵的代码如下：

!*******************************************************************************
! 计算一致性切线刚度矩阵
function plastic_kinematic_compute_consistent_modl(mu,kappa, n_pt,hard_data,alpha_n0,gamma_inc,q_n1_trial,n_trial) result(tant_modl)


  real(8),intent(in) :: mu,kappa
  integer,intent(in) :: n_pt
  real(8),intent(in) :: hard_data(n_pt,2)
  real(8),intent(in) :: alpha_n0
  real(8),intent(in) :: gamma_inc
  real(8),intent(in) :: q_n1_trial,n_trial(6)
  real(8)            :: tant_modl(6,6)


  ! 局部常数
  ! 四阶单位偏张量
  real(8),parameter :: I_dev(6,6) = reshape( [2.0/3.0, -1.0/3.0, -1.0/3.0, 0.0, 0.0, 0.0, &
                                              -1.0/3.0, 2.0/3.0, -1.0/3.0, 0.0, 0.0, 0.0, &
                                              -1.0/3.0, -1.0/3.0, 2.0/3.0, 0.0, 0.0, 0.0, &
                                              0.0, 0.0, 0.0, 0.5, 0.0, 0.0, &
                                              0.0, 0.0, 0.0, 0.0, 0.5, 0.0, &
                                              0.0, 0.0, 0.0, 0.0, 0.0, 0.5], [6,6])
  ! 四阶单位球张量
  real(8),parameter :: I_vol(6,6) = reshape( [1.0, 1.0, 1.0, 0.0, 0.0, 0.0, &
                                              1.0, 1.0, 1.0, 0.0, 0.0, 0.0, &
                                              1.0, 1.0, 1.0, 0.0, 0.0, 0.0, &
                                              0.0, 0.0, 0.0, 0.0, 0.0, 0.0, &
                                              0.0, 0.0, 0.0, 0.0, 0.0, 0.0, &
                                              0.0, 0.0, 0.0, 0.0, 0.0, 0.0], [6,6])
  ! 局部变量
  real(8) :: alpha_n1
  real(8) :: Hk
  real(8) :: nn(6,6)
  real(8) :: term1,term2
  integer :: i,j


  alpha_n1 = alpha_n0 + gamma_inc
  Hk = plastic_kinematic_compute_H(n_pt,hard_data,alpha_n1)


  nn = 0.0
  do i =1,6
    do j =1,6
      nn(i,j) = n_trial(i) * n_trial(j)
    enddo
  enddo


  term1 = 2.0 * mu * (1.0 - gamma_inc * 3.0 * mu / q_n1_trial)
  term2 = 6.0 * mu * mu * (gamma_inc / q_n1_trial - 1.0 / (3.0 * mu + Hk))
  
  tant_modl = 0.0
  tant_modl = term1 * I_dev + term2 * nn + kappa * I_vol


  return
end function plastic_kinematic_compute_consistent_modl

2 算法框图

Return mapping的算法如下：

3 umat源代码

umat采用Fortran90进行编写，其主程序的代码为：

include "plastic_kinematic_pack.f90"


subroutine UMAT(stress,statev,ddsdde,sse,spd,scd,                       &
                rpl,ddsddt,drplde,drpldt,                               &
                stran,dstran,time,dtime,temp,dtemp,predef,dpred,cmname, &
                ndi,nshr,ntens,nstatv,props,nprops,coords,drot,pnewdt,  &
                celent,dfgrd0,dfgrd1,noel,npt,layer,kspt,kstep,kinc)


    use plastic_kinematic_pack


    include 'aba_param.inc'


    character*80 cmname
    dimension    stress(ntens),statev(nstatv),ddsdde(ntens,ntens),   &
                ddsddt(ntens),drplde(ntens),                        &
                stran(ntens),dstran(ntens),time(2),                 &
                predef(1),dpred(1),props(nprops),coords(3),         &
                drot(3,3),dfgrd0(3,3),dfgrd1(3,3)


    !*******************************************************************************
    ! 材料参数
    integer :: n_pt                ! 硬化曲线的数据点个数
    real(8) :: hard_data(int((nprops-2)/2),2)   ! 硬化曲线的数据点表格
    real(8) :: E,nu
    real(8) :: mu,kappa
    real(8) :: sig_y0


    ! 从props数组中读取材料参数
    n_pt = int((nprops-3)/2)
    E = props(1)
    nu = props(2)
    sig_y0 = props(3)
    j=4
    do i = 1, n_pt
        hard_data(i,1) = props(j)
        hard_data(i,2) = props(j+1)
        j = j + 2
    enddo


    ! 更新应力，状态变量以及一致性切线刚度模量
    mu = E / (1.0 + nu) / 2.0
    kappa = E / (1.0 - 2.0*nu) / 3.0
    call plastic_kinematic_umat(mu,kappa,sig_y0, n_pt,hard_data, stran,dstran, stress,statev,ddsdde)


    return
end subroutine UMAT

相应的函数放在单独的一个f90文件的module中，用于调用，以实现主程序的整洁。

4 测试

4.1 一个单元加卸载测试

设置Abaqus自带线性随动硬化的本构为：

使用umat设置的材料参数为：

分别代表杨氏模量、泊松比，初始屈服应力，以及等效塑性应变与随动屈服应力的数据点。对于线性随动硬化模型，可以选取三个数据点，保证三点处于同一直线上，对最后一组数据点进行一个特殊处理，可以选取一个很大的塑性应变值，以保证计算过程中的等效塑性应变都落在这三个数据点点，由此插值得到便满足线性关系。

二者应力应变滞回曲线对比如图：

二者等效塑性应变演化曲线对比如图：

二者塑性耗散演化曲线对比如图：

4.2 带孔板拉伸测试

设置的材料参数为：

使用Abaqus计算的von Mises应力以及等效塑性应变的结果如下：

使用umat计算的von Mises应力以及等效塑性应变的结果如下：

可以看出，二者的结果是完全相同的。