如何将三维弹塑性本构应用于平面应力问题中

dearjj 6月7日浏览：875 收藏：4

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