elmt300.f90

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

                    DO jnode = 1, nnode
                      idofe = ( jnode - 1 ) * ndofn + idofn
                      dispx( idofx, kload ) = dispx( idofx, kload ) +          &
                      shape(     4, jnode ) * dispe( idofe        )
                    END DO
                  END IF
                END IF
              END DO
            END IF
          END DO
          kelem = kelem + 1
        END DO
        END

        SUBROUTINE GMat300( gmatx, shape, nnode )
!........
!         P R O G R A M
!             To calculate G matrix of 3-dimensional element.
!........
        IMPLICIT DOUBLE PRECISION ( a-h, o-z )
        DIMENSION shape ( 4, 21 ), gmatx ( 9, 63 )
        DO idofe = 1, 63
          DO istre = 1, 9
            gmatx( istre, idofe ) = 0.0D0
          END DO
        END DO
        DO inode = 1, nnode
          locat = ( inode - 1 ) * 3
          gmatx( 1, locat + 1 ) = shape( 1, inode )
          gmatx( 2, locat + 2 ) = shape( 1, inode )
          gmatx( 3, locat + 3 ) = shape( 1, inode )
          gmatx( 4, locat + 1 ) = shape( 2, inode )
          gmatx( 5, locat + 2 ) = shape( 2, inode )
          gmatx( 6, locat + 3 ) = shape( 2, inode )
          gmatx( 7, locat + 1 ) = shape( 3, inode )
          gmatx( 8, locat + 2 ) = shape( 3, inode )
          gmatx( 9, locat + 3 ) = shape( 3, inode )
        END DO
        RETURN
        END

        SUBROUTINE AMat300( amatx, shape )
!.........
!       P R O G R A M
!             To calculate A matrix of 3-dimensional element.
!.........
        USE ElmtData
        IMPLICIT DOUBLE PRECISION ( a-h, o-z )
        DIMENSION amatx( 6, 9 ), shape( 4, 21 ), awork( 3, 3 )

        DO istre = 1, 6
          DO jstre = 1, 9
            amatx( istre, jstre ) = 0.0D0
          END DO
        END DO

        DO idime = 1, 3
          DO jdime = 1, 3
            awork( idime, jdime ) = 0.0D0
            DO inode = 1, nnode
              locat = ( inode - 1 ) * 3
              awork( idime, jdime ) = awork( idime,  jdime ) +                 &
              shape( idime, inode ) * dispe( locat + jdime )
            END DO
          END DO
        END DO

        DO jdime = 1, 3
          amatx( 1, jdime     ) = awork( 1, jdime )
          amatx( 2, jdime + 3 ) = awork( 2, jdime )
          amatx( 3, jdime + 6 ) = awork( 3, jdime )
          amatx( 4, jdime + 3 ) = awork( 3, jdime )
          amatx( 4, jdime + 6 ) = awork( 2, jdime )
          amatx( 5, jdime     ) = awork( 3, jdime )
          amatx( 5, jdime + 6 ) = awork( 1, jdime )
          amatx( 6, jdime     ) = awork( 2, jdime )
          amatx( 6, jdime + 3 ) = awork( 1, jdime )
        END DO
        RETURN
        END

        SUBROUTINE SMat300( smatx, strgp )
!.......
!         P R O G R A M
!             To calculate S matrix of 3-dimensional element.
!.......
        IMPLICIT DOUBLE PRECISION ( a-h, o-z )
        DIMENSION smatx( 9, 9 ), strgp ( 6 )
        DO istre = 1, 9
          DO jstre = 1, 9
            smatx ( istre , jstre ) = 0.0D0
          END DO
        END DO

        DO istrs = 1, 3
          smatx( istrs    , istrs     ) = strgp( 1 )
          smatx( istrs + 3, istrs + 3 ) = strgp( 2 )
          smatx( istrs + 6, istrs + 6 ) = strgp( 3 )
          smatx( istrs    , istrs + 3 ) = strgp( 6 )
          smatx( istrs    , istrs + 6 ) = strgp( 5 )
          smatx( istrs + 3, istrs + 6 ) = strgp( 4 )
        END DO

        DO istre = 1, 9
          DO jstre = 1, istre - 1
            smatx ( istre , jstre ) = smatx (jstre , istre )
          END DO
        END DO

        RETURN
        END

        SUBROUTINE AGmat300( amatx, agmat, gmatx, ndofe )
!.........
!       P R O G R A M
!             To calculate AG matrix of 3-dimensional element.
!.........
        IMPLICIT DOUBLE PRECISION ( a-h, o-z )
        DIMENSION gmatx( 9, 63 ), amatx( 6, 9 ),agmat( 6, 63 )
        DO istre = 1, 6
          DO idofe = 1, ndofe
            agmat( istre, idofe ) = 0.0D0
            DO jstre = 1, 9
              agmat( istre, idofe ) = agmat( istre, idofe ) +                  &
                                      amatx( istre, jstre ) *                  &
                                      gmatx( jstre, idofe )
            END DO
          END DO
        END DO
        RETURN
        END

        SUBROUTINE Vsel300( imatb, elast, shear, poiss, vnorm )
        USE CtrlData
        USE MeshData
        IMPLICIT DOUBLE PRECISION( a-h, o-z )

        qmt11 = props(  2, imatb )
        qmt22 = props(  3, imatb )
        qmt33 = props(  4, imatb )
        qmt12 = props(  5, imatb )
        qmt23 = props(  6, imatb )
        qmt13 = props(  7, imatb )

        pcntd = qmt33
        pcnt1 = qmt11
        gcnt0 = qmt22
        gcnt1 = qmt11 - qmt22
        IF( timec .LT. 1.0D-12 ) THEN
          qcnt2 = 1.0D0
        ELSE IF( pcntd .GT. 1.0D-12 ) THEN
          pcnt2 = gcnt0 * pcnt1 / ( pcnt1 - gcnt0 )
          belta = ( pcnt1 + pcnt2 ) / pcntd * tstep
          qcnt2 = ( 1.0D0 - DEXP( -belta ) ) / belta
        ELSE
          qcnt2 = 0.0D0
        END IF
        pmodu =   gcnt0 + gcnt1 * qcnt2

        ecntd = qmt13
        ecnt1 = qmt12
        gcnt0 = qmt23
        gcnt1 = qmt12 - qmt23
        IF( timec .LT. 1.0D-12 ) THEN
          qcnt2 = 1.0D0
        ELSE IF( ecntd .GT. 1.0D-12 ) THEN
          ecnt2 = gcnt0 * ecnt1 / ( ecnt1 - gcnt0 )
          belta = ( ecnt1 + ecnt2 ) / ecntd * tstep
          qcnt2 = ( 1.0D0 - DEXP( -belta ) ) / belta
        ELSE
          qcnt2 = 0.0D0
        END IF
        shear = gcnt0 + gcnt1 * qcnt2
        elast = 9.0D0 * pmodu * shear / ( 3.0D0 * pmodu + shear )
        poiss = ( 0.5D0 * elast - shear ) / shear
        vnorm = elast
        RETURN
        END

        SUBROUTINE Strn300( cordg, strng, shape, nlinc, iswth )
!........
!       模块功能
!           计算Gauss点的应变，iswth == 1 当前位移应变(dispe)
!                                    == 2 结构变化前的应变(dishe)
!                                    == 3 上述应变和
!........
        USE CtrlData
        USE ElmtData
        IMPLICIT DOUBLE PRECISION( a-h,o-z )
        DIMENSION strng( 6 ), bmatx( 6, 63 )
        DIMENSION cordg( 3 ), shape( 4, 21 )

        CALL Bmat300( shape, bmatx, nlinc, 1 )
        IF( nerrc .GT. 0 ) RETURN

        DO idime = 1, ndime
          cordg( idime ) = 0.0D0
          DO inode = 1, nnode
            cordg(    idime ) = cordg(        idime ) +                        &
            shape( 4, inode ) * coren( idime, inode )
          END DO
        END DO

        DO istre = 1, 6
          strng( istre ) = 0.0D0
          IF( iswth .EQ. 1 .OR. iswth .EQ. 3 ) THEN
            DO idofe = 1, ndofe
              strng( istre        ) =   strng( istre ) +                       &
              bmatx( istre, idofe ) *   dispe( idofe )
            END DO
          END IF

          IF( iswth .EQ. 2 .OR. iswth .EQ. 3 ) THEN
            DO idofe = 1, ndofe
              strng( istre        ) = strng( istre ) +                         &
              bmatx( istre, idofe ) * dishe( idofe )
            END DO
          END IF
        END DO
        RETURN
        END

        SUBROUTINE Strs300( corep, strng, strsg, dmatx, xloca,                 &
                           yloca, imtyp, kgaus, kelem, imatb,                  &
                           ietyp, iswth )
        USE CtrlData
        USE MeshData
        USE ExtraMesh
        IMPLICIT DOUBLE PRECISION( a-h, o-z )
        DIMENSION strng( 6 ), corep( 3, 8 ), xaxis( 3 )
        DIMENSION strsg( 6 ), dmatx( 6, 6 ), zaxis( 3 )
        DIMENSION vwork( 6 ), tmatx( 6, 6 ), yaxis( 3 )

        IF( imtyp .EQ. 0 ) THEN
          DO istre = 1, 6
            strsg( istre ) = 0.0D0
            DO jstre = 1, 6
              strsg( istre        ) = strsg( istre ) +                         &
              dmatx( istre, jstre ) * strng( jstre )
            END DO
          END DO
        ELSE IF( imtyp .EQ. 1 ) THEN
          iloca = ( kgaus - 1 ) * 8
          IF( ietyp .GT. 0 ) THEN
            CALL Axis300( xaxis, yaxis, zaxis, xloca, yloca,                   &
                         corep, kelem )
            CALL TMatrix300( xaxis, yaxis, zaxis, tmatx )
            DO istre = 1, 6
              vwork( istre ) = 0.0D0
              DO jstre = 1, 6
                vwork( istre        ) = vwork( istre ) +                       &
                tmatx( istre, jstre ) * strng( jstre )
              END DO
            END DO
            DO istre = 1, 6
              strng( istre ) = vwork( istre )
            END DO
            strng( 3 ) = histx( iloca + 8, kelem )
            IF( ietyp .GE. 3 ) strng( 2 ) = strng( 3 )
          END IF

          qmt11 = props(  2, imatb )
          qmt22 = props(  3, imatb )
          qmt33 = props(  4, imatb )
          qmt12 = props(  5, imatb )
          qmt23 = props(  6, imatb )
          qmt13 = props(  7, imatb )

          pvold = qmt33
          pvol1 = qmt11
          gvol0 = qmt22
          gvol1 = qmt11 - qmt22
          IF( iswth .EQ. 0 ) THEN
            rvol2 = 1.0D0
          ELSE IF( pvold .GT. 1.0D-12 ) THEN
            pvol2 =   gvol0 * pvol1 / ( pvol1 - gvol0 )
            alpha = ( pvol1 + pvol2 ) / pvold * tstep
            rvol2 = DEXP( -alpha )
          ELSE
            rvol2 = 0.0D0
          END IF

          pshrd = qmt13
          pshr1 = qmt12
          gshr0 = qmt23
          gshr1 = qmt12 - qmt23
          IF( iswth .EQ. 0 ) THEN
            rshr2 = 1.0D0
          ELSE IF( pshrd .GT. 1.0D-12 ) THEN
            pshr2 = gshr0 * pshr1 / ( pshr1 - gshr0 )
            belta = ( pshr1 + pshr2 ) / pshrd * tstep
            rshr2 = DEXP( -belta )
          ELSE
            rshr2 = 0.0D0
          END IF

          cnst1 = gvol0 + 4.0D0 * gshr0 / 3.0D0
          cnst2 = gvol0 - 2.0D0 * gshr0 / 3.0D0
          cnst3 = 4.0D0 * gshr1 * rshr2 / 3.0D0
          cnst4 = 2.0D0 * gshr1 * rshr2 / 3.0D0
          cnst5 = gvol1 * rvol2
          cnst6 = gshr1 * rshr2

          DO istre = 1, 3
            jstre = istre + 1
            kstre = jstre + 1
            IF( jstre .GT. 3 ) jstre = jstre - 3
            IF( kstre .GT. 3 ) kstre = kstre - 3
            strsg( istre ) = cnst1 * strng( istre ) +                          &
                             cnst2 * strng( jstre ) +                          &
                             cnst2 * strng( kstre ) +                          &
                             cnst3 * histx( iloca + istre, kelem ) -           &
                             cnst4 *