elmt003.f90

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

        SUBROUTINE Elmt003( ielem, imats, iswth, iuses )
!........
!       模块功能
!           二维梁单元.
!........
        USE CtrlData
        USE ElmtData
        USE MeshData
        USE FactorData
        IMPLICIT DOUBLE PRECISION( a-h, o-z )
        DIMENSION rotat( 6, 6 ), fwork( 6, 6 ), direc( 2 )
        IF( iswth .EQ. 0 ) THEN
          IF( nprnc .NE. 0 ) WRITE( 12, 2000 )
          READ( 11, 2200 ) elast, vbend, varea, qforc, denst
          IF( nprnc .NE. 0 ) WRITE( 12, 2400 )  elast, vbend,                  &
                                         varea, qforc, denst
          props( 2, imats ) = elast
          props( 3, imats ) = vbend
          props( 4, imats ) = varea
          props( 5, imats ) = qforc
          props( 6, imats ) = denst
          RETURN
        ELSE IF( iswth .EQ. 8 ) THEN
          nnode = 2
          ndofn = 3
          RETURN
        ELSE IF( iswth .EQ. 9 ) THEN
          IF( iuses .EQ. 0 ) RETURN
          WRITE( 15, 2600 ) 2, ielem, imats, 2, lnode( 1 ), lnode( 2 )
          RETURN
        ELSE IF( iswth .EQ. 21 ) THEN
          WRITE( 25 ) 0
          RETURN
        END IF

!.......取单元材料性质.
        elast = props( 2, imats )   !弹性模量
        vbend = props( 3, imats )   !弯曲刚度
        varea = props( 4, imats )   !截面面积
        qforc = props( 5, imats )   !分布载荷
        denst = props( 6, imats )   !质量密度
!.......计算轴线方向和单元长度.
        fleng = 0.0D0
        DO idime = 1, 2
          direc( idime ) = coren( idime, 2 ) - coren( idime, 1 )
          fleng = fleng + direc( idime ) * direc( idime )
        END DO
        fleng = DSQRT( fleng )
        direc( 1 ) = direc( 1 ) / fleng
        direc( 2 ) = direc( 2 ) / fleng
!.......矩阵初始化.
        DO idofe = 1, 6
          DO jdofe = 1, 6
            fwork( idofe, jdofe ) = 0.0D0
            rotat( idofe, jdofe ) = 0.0D0
            stife( idofe, jdofe ) = 0.0D0
          END DO
        END DO
        IF( iuses .EQ. 0 ) RETURN
!.......形成坐标旋转矩阵.
        rotat( 1, 1 ) = direc( 1 )
        rotat( 2, 2 ) = direc( 1 )
        rotat( 1, 2 ) = direc( 2 )
        rotat( 2, 1 ) =-direc( 2 )
        rotat( 4, 4 ) = direc( 1 )
        rotat( 5, 5 ) = direc( 1 )
        rotat( 4, 5 ) = direc( 2 )
        rotat( 5, 4 ) =-direc( 2 )
        rotat( 3, 3 ) = 1.0D0
        rotat( 6, 6 ) = 1.0D0
        IF( iswth .EQ. 1  .OR. iswth .EQ.  5 .OR.                              &
            iswth .EQ. 10 .OR. MOD( iswth, mswth ) .EQ. 7 ) THEN
!.........形成单元刚度矩阵.
          stife( 1, 1 ) = elast * varea / fleng
          stife( 3, 3 ) = 4.0D0 * elast * vbend / fleng
          stife( 6, 3 ) = 2.0D0 * elast * vbend / fleng
          stife( 2, 2 ) = 12.D0 * elast * vbend / fleng ** 3
          stife( 3, 2 ) = 6.0D0 * elast * vbend / fleng ** 2
          stife( 4, 4 ) = stife( 1, 1 )
          stife( 4, 1 ) =-stife( 1, 1 )
          stife( 5, 5 ) = stife( 2, 2 )
          stife( 5, 2 ) =-stife( 2, 2 )
          stife( 6, 2 ) = stife( 3, 2 )
          stife( 5, 3 ) =-stife( 3, 2 )
          stife( 6, 5 ) = stife( 5, 3 )
          stife( 6, 6 ) = stife( 3, 3 )
          IF( nincc .EQ. 2 .OR. iswth .EQ. 10 ) THEN
            DO idofe = 1, 6
              DO jdofe = 1, 6
                force( idofe ) = force( idofe ) -                              &
                dispe( jdofe ) * stife( idofe, jdofe )
              END DO
            END DO
          END IF
        ELSE IF( iswth .EQ. 2 ) THEN
!.........形成单元质量矩阵.
          const = fleng * varea * denst
          stife( 1, 1 ) = const / 3.0d0
          stife( 4, 4 ) = const / 3.0d0
          stife( 4, 1 ) = const / 6.0d0
          stife( 2, 2 ) = const * 13.d0 / 35.d0
          stife( 5, 5 ) = const * 13.d0 / 35.d0
          stife( 5, 2 ) = const * 9.0d0 / 70.d0
          stife( 3, 2 ) =-const * 11.d0 * fleng / 210.d0
          stife( 6, 2 ) = const * 13.d0 * fleng / 420.d0
          stife( 5, 3 ) =-const * 13.d0 * fleng / 420.d0
          stife( 6, 5 ) = const * 11.d0 * fleng / 210.d0
          stife( 3, 3 ) = const * fleng * fleng / 105.d0
          stife( 6, 6 ) = const * fleng * fleng / 105.d0
          stife( 6, 3 ) =-const * fleng * fleng / 140.d0
        ELSE IF( iswth .EQ. 3 ) THEN
        ELSE IF( iswth .EQ. 4 ) THEN
!.........形成单元几何矩阵.
          displ = ( dispe( 4 ) - dispe( 1 ) ) * direc( 1 ) +                   &
                  ( dispe( 5 ) - dispe( 2 ) ) * direc( 2 )
          faxis = displ * elast * varea / ( fleng * fleng )
          stife( 2, 2 ) = faxis * 6.0d0 / 5.0d0
          stife( 3, 2 ) = faxis * fleng / 10.d0
          stife( 6, 3 ) =-faxis * fleng * fleng / 30.d0
          stife( 3, 3 ) = faxis * fleng * fleng * 2.0d0 / 15.d0
          stife( 5, 5 ) = stife( 2, 2 )
          stife( 5, 2 ) =-stife( 2, 2 )
          stife( 6, 2 ) = stife( 3, 2 )
          stife( 5, 3 ) =-stife( 3, 2 )
          stife( 6, 5 ) =-stife( 3, 2 )
          stife( 6, 6 ) = stife( 3, 3 )
        ELSE IF( iswth .EQ. 6 ) THEN
!.........形成节点力矢量.
          fwork( 1, 1 ) = 0.0D0
          fwork( 4, 1 ) = 0.0D0
          fwork( 2, 1 ) = 0.5D0 * fleng * qforc
          fwork( 5, 1 ) = 0.5D0 * fleng * qforc
          fwork( 3, 1 ) = qforc * fleng * fleng / 12.D0
          fwork( 6, 1 ) =-qforc * fleng * fleng / 12.D0
!.........节点力矢量坐标变换.
          DO idofe = 1, 6
            force( idofe ) = 0.0D0
            DO jdofe = 1, 6
              force( idofe ) = force( idofe ) +                                &
              rotat( jdofe, idofe ) * fwork( jdofe, 1 )
            END DO
          END DO
          RETURN
        END IF
!.......形成对称部分.
        DO idime = 1, 6
          DO jdime = 1, idime
            stife( jdime, idime ) = stife( idime, jdime )
          END DO
        END DO
!.......单元矩阵坐标变换.
        DO idofe = 1, 6
          DO jdofe = 1, 6
            fwork( idofe, jdofe ) = 0.0D0
            DO kdofe = 1, 6
              fwork( idofe, jdofe ) = fwork( idofe, jdofe ) +                  &
              stife( idofe, kdofe ) * rotat( kdofe, jdofe )
            END DO
          END DO
        END DO
        IF( MOD( iswth, mswth ) .NE. 7 ) THEN
          DO idofe = 1, 6
            DO jdofe = 1, 6
              stife( idofe, jdofe ) = 0.0d0
              DO kdofe = 1, 6
                stife( idofe, jdofe ) = stife( idofe, jdofe ) +                &
                rotat( kdofe, idofe ) * fwork( kdofe, jdofe )
              END DO
            END DO
          END DO
        END IF
        DO idofe = 1, 6
          force( idofe ) = 0.0D0
          DO jdofe = 1, 6
            force( idofe ) = force( idofe ) +                                  &
            dispe( jdofe ) * fwork( idofe, jdofe )
          END DO
        END DO
        IF( MOD( iswth, mswth ) .EQ. 7 ) THEN
          WRITE( 12, 2800 ) ielem, force
        ELSE
          IF( nincc .EQ. 2 ) THEN
            DO idofe = 1, 6
              force( idofe ) =-force( idofe )
            END DO
          ELSE
            DO idofe = 1, 6
              force( idofe ) = 0.0D0
            END DO
          END IF
        END IF
2000    FORMAT( 2X, '单元类型:            二维梁单元' )
2200    FORMAT( 4E15.5 / 4E15.5 )
2400    FORMAT( 2x, '杨氏模量: ', E15.5/2x, '惯性矩:   ', E15.5 /              &
                2x, '截面面积: ', E15.5/2x, '均布载荷: ', E15.5 /              &
                2x, '质量密度: ', E15.5 )
2600    FORMAT( 10I5 )
2800    FORMAT( 2X, I5, 6E12.6 )
        RETURN
        END