Gurson塑性本构理论和umat源代码

1 一般塑性理论

从一般的塑性理论入手，然后将Gurson塑性本构代入到该一般理论框架中。从应变分解开始：

从热力学角度推导本构方程和耗散不等式。自由能是全应变、塑性应变（作为内变量）和其他内变量的函数：

代入到Clausius–Duhem不等式中可得：

因此可以得到本构方程：

以及耗散不等式：

屈服函数表示为热力学共轭力的函数：

流动法则和硬化法则的表达式如下：

最后总结为：

Return-mapping算法框架如下：

一致性切线刚度模量的推导如下。对于上框中的非线性return-mapping方程，将系统输入应变也当做一个自变量，然后对方程求微分可得：

由表达式：

求逆可得：

式中各线性算子的表达式如下：

由于自由能的分解：

有：

式中：

因此可得到方程：

求逆的：

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

2 本构理论

此理论部分来自于Abaqus官方理论文档。

2.1 屈服函数

屈服函数为：

有关系式：

式中q为von Mises应力：

S为偏应力：

p为静水压：

式中：

为matrix相（solid相）的屈服应力，为matrix等效累积塑性应变的函数。

2.2 流动法则和硬化法则

塑性应变的流动法则为：

对比上节中的一般塑性理论，即可得：

为塑性流动方向。对于两个内变量：体积分数和基体等效累积塑性应变，由上节的一般的塑性理论可以写为如下的由硬化模量表示的表达式：

matrix的等效累积塑性应变的硬化法则如下：

可得：

孔隙体积分数的演化分为两部分，分别void growth和void nucleation：

因此可以得到表达式：

2.3 应力更新算法

首先计算试验状态下的值（elastic predictor），

如果试验状态下屈服函数的值小于0，则代表试验状态即为真实状态；如果试验状态下屈服函数的值大于0，则需要进行塑性更正（plastic corrector）。根据一般塑性理论框架可得需要求解的方程为：

未知量为：

使用牛顿迭代法求解上述非线性方程组。其雅克比矩阵为（略去下标n+1）：

式中涉及的导数依次如下：

和：

和：

和：

3 UMAT代码

采用Fortran90进行编写，其主程序代码框架如下：

include "Gurson_umat_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 Gurson_umat_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)

    !*******************************************************************************
    ! 材料参数
    real(8) :: E,nu
    real(8) :: q_1,q_2, f_N,s_N,e_N
    integer :: n_pt                ! 硬化曲线的数据点个数
    real(8) :: hard_data(int((nprops-7)/2),2)   ! 硬化曲线的数据点表格
    real(8) :: mu,kappa
    real(8) :: strain_n1(6)
    integer :: i,j
    

    ! 从props数组中读取材料参数
    E = props(1)
    nu = props(2)
    q_1 = props(3)
    q_2 = props(4)
    f_N = props(5)
    s_N = props(6)
    e_N = props(7)
    n_pt = int((nprops-7)/2)

    ! 读取硬化数据
    j=8
    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
    strain_n1(1:6) = stran(1:6) + dstran(1:6)
    call GT_umat(n_pt, hard_data, E,nu,q_1,q_2, f_N,s_N,e_N, strain_n1, stress, nstatev, statev, ddsdde)
    ! write(*,"(A,T20,6ES16.8)") "stress_n1 = ",stress(1:6)

    return
end subroutine UMAT

相应的函数单独放于文件"Gurson_umat_pack.f90"中，以保证umat主程序的整洁。

4 测试

4.1 一个单元单轴拉伸

对一个单元进行单轴拉伸测试，应力应变曲线如下：

孔洞体积分数的演化曲线如下：

基体相等效塑性应变的演化曲线如下：

4.2 一个单元加卸载

对一个单元进行加卸载试验，应力应变曲线如下：

孔洞体积分数的演化曲线如下：

基体相等效塑性应变的演化曲线如下：

4.3 带孔板拉伸

对一带孔板进行单轴拉伸测试，对比结果如下（左图为Abaqus内置Gurson本构，右图为umat结果）

Von Mise应力的结果如下：

孔洞体积分数的结果如下：

基体等效塑性应变的结果如下：

塑性应变11分量如下：

塑性应变22分量如下：

结果与Abaqus基本相同

5 参考书籍

Neto, E. A. de Souza, D. R. J. Owen, and D. Peric. , 'Computational Methods for Plasticity: Theory and Applications'