如何將三維彈塑性本構應用于平面應力問題中

dearjj

2024年6月7日 15:03

瀏覽：2563 評論：3 收藏：13

1 本構理論

本文講解如何將三維的率無關彈塑性理論應用到平面應力問題中。對于平面應變和軸對稱問題，由于是相應的應變分量為0，因為可以直接使用三維的本構，只需將相應的應變分量設為0作為本構的輸入即可。然后，對于平面應力問題，是相應的應力分量為0，由于本構是由應變驅動求得對應的應力，相應應力分量為0相當于對系統施加了相應的約束，因此三維的本構理論不可直接應用于平面應力問題中，需要將相應的約束考慮其中進行求解。

1.1 平面應力理論

對于線彈性情況，由三維本構方程推導平面應力方程如下：

如何將三維彈塑性本構應用于平面應力問題中的圖1

1.2 應力更新算法

采用一種嵌套迭代的方法進行應力更新。我們將平面外應變仍然作為本構的輸入，此時可調用三維的本構方程，得到對應的應力。如果得到的平面外應力不為0，則使用牛頓迭代法對平面外應變進行更新，持續此過程，直至滿足平面應力假設。算法流程如下：

如何將三維彈塑性本構應用于平面應力問題中的圖2

1.3 一致性切線剛度矩陣

如何將三維彈塑性本構應用于平面應力問題中的圖3

2 UMAT代碼

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

    ! 狀態變量
    real(8) :: statev_dummy(8)
    ! 調用三維塑性理論
    real(8) :: strain_n0(6)
    real(8) :: strain_inc(6)
    real(8) :: stress_n1(6)
    real(8) :: ddsdde_n1(6,6)

    ! 平面應力循環
    real(8) :: strain_33
    real(8) :: stress_33
    real(8) :: jac
    integer :: i,j
    real(8),parameter :: ITER_TOL = 0.1**12
    integer,parameter :: ITER_MAX = 100

    ! 計算一致性切線剛度模量所需變量
    real(8) :: D11(3,3)
    real(8) :: D12(3)
    real(8) :: D21(3)
    real(8) :: D1221(3,3)
    real(8) :: D22

    
    ! 從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

    ! 調用三維塑性本構
    strain_33 = statev(9)
    strain_n0(1:6) = 0.0d0
    strain_n0(1:2) = stran(1:2) ! 11,22
    strain_n0(3) = strain_33
    strain_n0(4) = stran(3) ! 12


    strain_inc(1:6) = 0.0d0 
    strain_inc(1:2) = dstran(1:2) ! 11,22
    strain_inc(4) = dstran(3) ! 12

    ! 平面應力循環
    do i=1,ITER_MAX
        statev_dummy(1:8) = statev(1:8)
        ! resume stress
        stress_n1(1:2) = stress(1:2)
        stress_n1(4) = stress(3)
        call plastic_iso_umat(mu,kappa, n_pt,hard_data, strain_n0,strain_inc, stress_n1,statev_dummy,ddsdde_n1)
        stress_33 = stress_n1(3)
        if (abs(stress_33) < ITER_TOL) then
            exit
        endif
        jac = ddsdde_n1(3,3)
        strain_33 = strain_33 - stress_33 / jac
        strain_n0(3) = strain_33
        stress_n1(3) = stress_33
    enddo

    ! 更新應力
    stress(1:2) = stress_n1(1:2)
    stress(3) = stress_n1(4)

    ! 計算一致性切線剛度模量
    D11(1:2,1:2) = ddsdde_n1(1:2,1:2)
    D11(1,3) = ddsdde_n1(1,4)
    D11(3,1) = D11(1,3)
    D11(2,3) = ddsdde_n1(2,4)
    D11(3,2) = D11(2,3)
    D11(3,3) = ddsdde_n1(4,4)

    D12(1) = ddsdde_n1(1,3)
    D12(2) = ddsdde_n1(2,3)
    D12(3) = ddsdde_n1(4,3)

    D21(1) = ddsdde_n1(3,1)
    D21(2) = ddsdde_n1(3,2)
    D21(3) = ddsdde_n1(3,4)

    do i=1,3
        do j=1,3
            D1221(i,j) = D12(i) * D21(j)
        enddo
    enddo
    D22 = ddsdde_n1(3,3)
    ddsdde = D11 - D1221 / D22

    ! 更新狀態變量
    statev(1:8) = statev_dummy(1:8)
    statev(9) = strain_33

    return
end subroutine UMAT