1 本構理論

1.1 率形式

本構方程為：

單軸拉伸的應力應變的硬化曲線如下：

根據單軸試驗得到硬化部分的曲線：

當僅考慮隨動硬化時，屈服面的中心在移動，而屈服面的大小不發(fā)生改變，即為常數。

屈服條件為：

增加了背應力來表示屈服面中心移動即隨動硬化的效果。式中相對應力的表達式為：

流動法則為：

背應力的演化法則為：

式中：

上標k代表隨動部分(kinematic)，表示以下隨動硬化曲線的梯度：

1.2 Return mapping算法應力更新

在給定增量應變以及上一步的狀態(tài)變量的情況下，首先計算試驗狀態(tài)下的量：

計算試驗狀態(tài)下的屈服函數值，判斷是否發(fā)生屈服。當試驗屈服函數值大于0時，說明需要進行塑性更正，反正則說明試驗狀態(tài)即為真實狀態(tài)。以下對塑性更正環(huán)節(jié)進行詳細說明。

當發(fā)生塑性流動時，需要求解以下非線性方程組：

可以將上述非線性方程組簡化值一個非線性方程。首先將第一個關于應變的方程代入本構方程中，可得偏應力的表達式為：

減去方程組中的第二式可以得到相對應力的表達式為：

可以看出試驗相對應力與真實相對應力之間滿足關系式：

那么可以將相對應力表示為：

式中：

最后將該相對應力代入到方程組中最后一個一致性方程中可得：

該非線性方程的未知量為塑性乘子增量，可以用牛頓迭代法進行求解。另外可以利用公式：

代入替換可得背應力和相對應力的公式：

以及最后的非線性方程為：

使用牛頓迭代法求解的公式為：

式中雅克比為：

式中：

是各向同性硬化曲線的梯度。當僅考慮隨動硬化時，該項為0。

1.3 一致性切線剛度模量

當沒有發(fā)現塑性流動時，一致性切線剛度矩陣即為彈性矩陣；當發(fā)生塑性流動時，應力更新公式為：

求導可得：

式中：

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

計算該一致性切線剛度矩陣的代碼如下：

!*******************************************************************************
! 計算一致性切線剛度矩陣
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


    ! 更新應力，狀態(tài)變量以及一致性切線剛度模量
    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應力以及等效塑性應變的結果如下：

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