各向同性硬化von Mises率無關彈塑性本構理論以及umat源代碼

1 本構理論

1.1 率形式

對于各向同性線彈性材料，其本構方程為：

式中假設了應變張量可以分解為彈性應變和塑性應變兩部分：

因此塑性本構的關鍵在于計算塑性應變的演化。對于率無關彈塑性的本構理論，需要確定以下三個部分：

(1)：屈服條件

(2)：流動法則

(3)：硬化法則

在此采用的是 von Mises 屈服條件：

式中后繼屈服應力是等效塑性應變的函數：

流動法則為：

式中流動方向的表達式為：

硬化法則為：

1.2 Return-mapping算法

上述的本構方程均為率形式。在增量步中，給定增量應變：

首先假設該增量應變全為彈性應變，計算試驗狀態下的一些物理量：

試驗狀態下的應力

試驗狀態下的屈服函數值：

利用該試驗屈服函數值來判斷在該增量步下是否發生了塑性屈服。如果：

則說明試驗狀態即為真實狀態，即可進行更新：

反之則需要進行塑性更正，即需要計算塑性乘子的增量，利用以下非線性方程組進行計算：

可以將該非線性方程組簡化至一個非線性方程，過程如下，將該方程組中的第一式分解為球量和偏量兩部分：

因此可以計算應力為：

將上式中的第二式整理得到：

可以得到兩個張量的方向相同：

因此偏應力可以用試驗狀態的信息表示出來：

代入到最后一個一致性方程中可得：

即可利用牛頓迭代法對上述非線性方程進行求解，得到塑性乘子增量。求解得到塑性乘子增量之后，即可更新：

也可以更新塑性應變：

1.3 一致性切線剛度矩陣

umat除了要求更新應力以及狀態變量之外，還需要更新算法的一致性切線剛度模量。當沒有發生塑性屈服時，一致性切線剛度矩陣即為彈性矩陣。當發生塑性屈服時，根據應力更新的公式：

可以計算：

式中：

式中：

為硬化曲線的梯度，各向同性硬化曲線使用數據點擬合為分段線性函數：

因此最終可得到一致性切線剛度矩陣的表達式為：

2 算法描述

整個 return-mapping 算法框圖如下：

牛頓迭代法算法框圖如下：

3 UMAT代碼

umat使用F90編寫，其主程序如下：

include "plastic_iso_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_iso_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
    


    ! 從props數組中讀取材料參數
    n_pt = int((nprops-2)/2)
    E = props(1)
    nu = props(2)
    j=3
    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_iso_umat(mu,kappa, n_pt,hard_data, stran,dstran, stress,statev,ddsdde)


    return
end subroutine UMAT

umat所需的支持函數編寫于單獨的module中（plastic_iso_pack.f90）用于調用，該方式可以使得主程序十分簡潔。

4 測試

4.1 一個單元單軸加卸載試驗

對一個單元進行單軸的加卸載試驗，繪制應力應變滯回曲線。Abaqus設置的材料參數為，彈性部分的性質：

塑性部分的性質：

umat設置相同的材料參數為：

以及需要設置狀態變量的個數為8。二者計算的結果對比如下：

二者等效塑性應變的演化對比圖如下：

二者塑性耗散的演化對比圖為：

4.2 帶孔板拉伸實驗

利用相同的材料性質計算一帶孔板的拉伸實驗。Abaqus計算的結果如下：

使用umat計算的結果如下：

可看到二者的結果是完全相同的。

5 參考書籍

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