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

1 本构理论

1.1 率形式

本构方程为：

当考虑混合硬化时，其应力应变曲线如下：

随动硬化导致屈服面的中心在移动，各向同性硬化导致屈服面在扩张。根据单轴试验得到两个硬化部分的曲线：

屈服条件为：

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

流动法则为：

屈服应力是等效塑性应变的函数：

背应力的演化法则为：

式中：

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

1.2 Return mapping算法应力更新

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

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

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

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

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

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

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

式中：

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

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

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

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

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

式中雅克比为：

式中：

是各向同性硬化曲线的梯度。

1.3 一致性切线刚度模量

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

求导可得：

式中：

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

式中指的是各向同性硬化曲线的梯度。

2 算法框图

Return mapping的算法如下：

3 umat源代码

include "plastic_combined_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_combined_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_i                ! 各向硬化曲线的数据点个数
    integer :: n_pt_k                ! 随动硬化曲线数据点的个数
    real(8),allocatable :: hard_data_i(:,:)     ! 各向硬化曲线的数据点表格
    real(8),allocatable :: hard_data_k(:,:)
    real(8) :: E,nu
    real(8) :: mu,kappa
    integer :: i,j
    
    ! 从props数组中读取材料参数
    E = props(1)
    nu = props(2)
    n_pt_i = int(props(3))
    n_pt_k = int(props(4))
    
    allocate(hard_data_i(n_pt_i,2))
    allocate(hard_data_k(n_pt_k,2))
    j=5
    do i = 1, n_pt_i
        hard_data_i(i,1) = props(j)
        hard_data_i(i,2) = props(j+1)
        j = j + 2
    enddo
    do i = 1, n_pt_k
        hard_data_k(i,1) = props(j)
        hard_data_k(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_combined_umat(mu,kappa, n_pt_i,hard_data_i, n_pt_k,hard_data_k, stran,dstran, stress,statev,ddsdde)


    return
end subroutine UMAT

设置的材料参数分别为杨氏模量、泊松比、各向同性硬化数据点个数、随动硬化数据点个数、各向同性硬化数据、随动硬化数据。umat采用Fortran90进行编写，其所需要的相关函数编写于一个单独的f90的module中，这样可以保持主程序的整洁。

4 测试

Abaqus的混合硬化中，背应力的演化准则采用的是Chaboche准则，我们这里的演化准则采用的是简单的线性Prager准则，因此在这里可能不与Abaqus的混合硬化结果进行对比。

文中的混合硬化两个极端分别为各向同性硬化和随动强化，构造数据满足上述两个极端并与Abaqus自带的塑性本构进行对比。

4.1 单个单元测试

首先使用该混合硬化的umat构造一个各向同性硬化的计算案例，设置材料参数为：

最后四个参数为两组随动屈服应力的数据点，这样设置能够使得插值得到的随动屈服应力为0（当然，最后一个数据点可以按需设置一个较大的等效塑性应变值），从而不含随动硬化现象。使用umat计算与Abaqus内置的本构计算对比如下：

应力应变曲线对比如下：

等效塑性应变的演化对比如下：

塑性耗散的对比如下：

然后使用该混合硬化的umat构造一个随动硬化的计算案例，设置材料参数为：

可见第5-8个参数能够使得屈服面大小的插值结果保持不变，即不各向同性硬化现象。使用umat计算与Abaqus内置的本构计算对比如下：

应力应变曲线对比如下：

等效塑性应变的演化对比如下：

塑性耗散的对比如下：

最后使用该混合硬化的umat构造一个混合硬化的计算案例，设置材料参数为：

应力应变曲线如下：

等效塑性应变的演化如下：

塑性耗散的如下：