plastic.f90

边界元程序,供力学工作者参考,希望对大家有所帮助

!       ===============================================================
!       |  文件功能：
!       |      本文件的各个子程序为实现弹塑性分析需要的公共模块
!       |      创建者：凌道盛
!       |      时  间：2002/04/13
!       |      变量说明：
!       |                sflow — 流动矢量(屈服函数对应力矢量的偏导数)
!       |                pflow — 塑性流动矢量(屈服函数对塑性应力矢量的偏导数)
!       |                dstrn — 应变增量
!       |                strsp — 预应力
!       |                dmate — 弹性矩阵
!       |                dmatp — 弹塑性矩阵
!       |                hardk — 强化规律(屈服函数对内时变量κ的偏导数)
!       ===============================================================


        DOUBLE PRECISION FUNCTION Lamda( dmate, strsg, strsp, dstrn,           &
                                  sflow, pflow, hardk, ietyp, itype )
!       计算dλ
        USE CtrlData
        IMPLICIT DOUBLE PRECISION( a-h, o-z )
        DIMENSION sflow( 6 ), pflow( 6 )
        DIMENSION vwork( 6 ), dstrn( 6 )
        DIMENSION strsp( 6 ), strsg( 6 ), dmate( 6, 6 )
!       计算(塑性)流动矢量转置与弹性矩阵的乘积
        DO istrs = 1, 6
          vwork( istrs ) = 0.0D0
          DO jstrs = 1, 6
            vwork( istrs ) = vwork( istrs ) +                                  &
            sflow( jstrs ) * dmate( jstrs, istrs )
          END DO
        END DO
!       先计算dλ的分子部分
        partn = 0.0D0
        DO istrs = 1, 6
          partn = partn + vwork( istrs ) * dstrn( istrs )
        END DO
!       再计算dλ的分母部分
        partd = 0.0D0
        DO istrs = 1, 6
          partd = partd + vwork( istrs ) * sflow( istrs )
          partd = partd - vwork( istrs ) * pflow( istrs )
        END DO
        dkdla = DkDlamda( sflow, strsg, strsp, ietyp, itype )
        IF( nerrc .NE. 0 ) RETURN
        partd = partd - hardk * dkdla
        Lamda = partn / partd
        IF( Lamda .LT. 0.0D0 ) Lamda = 0.0D0
        END

        FUNCTION DkDlamda( sflow, strsg, strsp, ietyp, itype )
!       ================================================================
!       | 模块功能：计算dκ/dλ
!       ================================================================
        USE CtrlData
        IMPLICIT DOUBLE PRECISION( a-h, o-z )
        DIMENSION strsg( 6 ), vwork( 6 )
        DIMENSION sflow( 6 ), vtemp( 6 )
        DIMENSION strsp( 6 ), vtran( 6 )

        CALL ReorderStress( strsg, vwork, ietyp, 0 )
        CALL ReorderStress( strsp, vtran, ietyp, 0 )
        CALL ReorderStress( sflow, vtemp, ietyp, 0 )
        DO istrs = 1, 6
          vwork( istrs ) = vwork( istrs ) + vtran( istrs )
        END DO

        reslt = 0.0D0
        SELECT CASE( itype )
          CASE( 1 )
!...........塑性功
            DO istrs = 1, 6
              reslt = reslt + vwork( istrs ) * vtemp( istrs )
            END DO
          CASE( 2 )
!...........塑性体积应变
            reslt = vtemp(1) + vtemp(2) + vtemp(3)
          CASE( 3 )
!...........塑性等效应变
            DO istrs = 1, 6
              reslt = reslt + vtemp( istrs ) * vtemp( istrs )
            END DO
            reslt = DSQRT( reslt )
          CASE DEFAULT
            nerrc = 3403546
            WRITE( 12, 2000 )
        END SELECT
        DkDlamda = reslt
2000    FORMAT( //2x, '致命错误：选择了无效的内时变量！' )
        END

        SUBROUTINE DEPMatix( dmate, dmatp, strsg, strsp, sflow,                &
                             pflow, hardk, ietyp, itype )
!       计算弹塑性矩阵Dep
        USE CtrlData
        IMPLICIT DOUBLE PRECISION( a-h, o-z )
        DIMENSION dmate( 6, 6 ), strsg( 6 ), sflow( 6 )
        DIMENSION dmatp( 6, 6 ), pflow( 6 ), vwork( 6 ), strsp( 6 )
!       计算(塑性)流动矢量转置与弹性矩阵的乘积
        DO istrs = 1, 6
          vwork( istrs ) = 0.0D0
          DO jstrs = 1, 6
            vwork( istrs ) = vwork( istrs ) +                                  &
            sflow( jstrs ) * dmate( jstrs, istrs )
            dmatp( istrs, jstrs ) = dmate( istrs, jstrs )
          END DO
        END DO
!       计算A的分母部分
        apara = 0.0D0
        DO istrs = 1, 6
          apara = apara + vwork( istrs ) * sflow( istrs )
          apara = apara - vwork( istrs ) * pflow( istrs )
        END DO
        dkdla = DkDlamda( sflow, strsg, strsp, ietyp, itype )
        IF( nerrc .NE. 0 ) RETURN
        apara = apara - hardk * dkdla
        DO istrs = 1, 6
          DO jstrs = 1, 6
            dmatp( istrs, jstrs ) = dmatp( istrs, jstrs ) -                    &
            vwork( istrs ) * vwork( jstrs ) / apara
          END DO
        END DO
        END