elmt004.f90

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

        END IF
2000    FORMAT( 2X, '单元类型:            三维梁单元' )
2200    FORMAT( 4E15.5 / 4E15.5 / 4E15.5 / 4E15.5 / 3I5 /                      &
                4E15.5 / 4E15.5 / 4E15.5 / 4E15.5 )
2400    FORMAT( 2x, '杨氏模量: ', E15.5/2x, '剪切模量: ', E15.5/               &
                2x, '截面面积: ', E15.5/2x, 'z-惯性矩: ', E15.5/               &
                2x, 'y-惯性矩: ', E15.5/2x, 'x-极性矩: ', E15.5/               &
                2x, 'y-剪切比: ', E15.5/2x, 'z-剪切比: ', E15.5/               &
                2x, '均载 Qx : ', E15.5/2x, '均载 Qy : ', E15.5/               &
                2x, '均载 Qz : ', E15.5/2x, '质量密度: ', E15.5/               &
                2x, '主转动角: ', E15.5/2x, '主参考点: ',3F15.5/               &
                2x, '应 力 点: ',2F15.5/12x, 2F15.5/12x, 2F15.5/               &
               12x, 2F15.5 )
2600    FORMAT( 10I5 )
2800    FORMAT(///18X, ' 单元', I5, ' 内力和应力(局部坐标系)'///6X,            &
                '节点      轴向力-Nx      剪切力-Qy      剪切力-Qz'/6X         &
                '          扭  矩-Mx      弯  矩-My      弯  矩-Mz' )
3000    FORMAT( 5X, I5, 3E15.5 / 10X, 3E15.5 )
3200    FORMAT(/6x, '节点', 9x, 'y-坐标', 9x, 'z-坐标', 11x, '应力' )
3400    FORMAT( 5x, I5, 2F15.5, E15.5 )
        END

        SUBROUTINE SetTheta( angle, imats )
        USE CtrlData
        USE ElmtData
        USE MeshData
        IMPLICIT DOUBLE PRECISION( a-h, o-z )
!.......取或设定单元参数
        angle = props( 14, imats )
        xrefp = props( 15, imats )
        yrefp = props( 16, imats )
        zrefp = props( 17, imats )
        IF( DABS( angle ) .LT. 360.0 ) THEN
!.........如果输入的转动角有效则返回
          angle = angle * DATAN( 1.0D0 ) / 45.D0
          RETURN
        END IF
        xleng = coren( 1, 2 ) - coren( 1, 1 )
        yleng = coren( 2, 2 ) - coren( 2, 1 )
        zleng = coren( 3, 2 ) - coren( 3, 1 )
        xylen = DSQRT( xleng * xleng + yleng * yleng )
        vleng = DSQRT( xleng * xleng + yleng * yleng + zleng * zleng )
        IF( xylen .LT. 1.0D-8 ) THEN
!.........如果局部坐标系的x-坐标与整体坐标的z-坐标一致
          sinth = 1.0D0
          costh = 0.0D0
          sinna = 1.0D0
          cosna = 0.0D0
        ELSE
!.........如果局部坐标系的x-坐标与整体坐标的z-坐标不一致
          sinth = yleng / xylen
          costh = xleng / xylen
          sinna = zleng / vleng
          cosna = xylen / vleng
        END IF
!.......对参考点进行坐标变换
        xrotp = costh * cosna * xrefp + sinth *                                &
                cosna * yrefp + sinna * zrefp
        yrotp =-sinth * xrefp + costh * yrefp
        zrotp =-costh * sinna * xrefp - sinna *                                &
                sinth * yrefp + cosna * zrefp
!.......根据参考点计算主转动角
        vleng = DSQRT( yrotp * yrotp + zrotp * zrotp )
        sinth = zrotp / vleng
        costh = yrotp / vleng
        angle = DASIN( sinth )
        IF( costh .LT.-1.0D-8 ) angle = 4.0D0 * DATAN( 1.0D0 ) - angle
        RETURN
        END

        SUBROUTINE TransMatrix( trans, angle )
        USE CtrlData
        USE ElmtData
        IMPLICIT DOUBLE PRECISION( a-h, o-z )
        DIMENSION trans( 12, 12 )
!.......
        DO idofe = 1, 12
          DO jdofe = 1, 12
            trans( idofe, jdofe ) = 0.0D0
          END DO
        END DO
!.......
        xleng = coren( 1, 2 ) - coren( 1, 1 )
        yleng = coren( 2, 2 ) - coren( 2, 1 )
        zleng = coren( 3, 2 ) - coren( 3, 1 )
        ftemp = DSQRT( xleng * xleng + yleng * yleng )
        fleng = DSQRT( ftemp * ftemp + zleng * zleng )
!.......
        trans( 1, 1 ) = xleng / fleng
        trans( 1, 2 ) = yleng / fleng
        trans( 1, 3 ) = zleng / fleng
        IF( ftemp .LE. 1.0D-8 ) THEN
          trans( 2, 1 ) =-1.0D0
          trans( 2, 2 ) = 0.0D0
          trans( 2, 3 ) = 0.0D0
        ELSE
          trans( 2, 1 ) =-yleng / ftemp
          trans( 2, 2 ) = xleng / ftemp
          trans( 2, 3 ) = 0.0D0
        END IF
        trans( 3, 1 ) = trans( 1, 2 ) * trans( 2, 3 ) -                        &
                        trans( 1, 3 ) * trans( 2, 2 )
        trans( 3, 2 ) = trans( 1, 3 ) * trans( 2, 1 ) -                        &
                        trans( 1, 1 ) * trans( 2, 3 )
        trans( 3, 3 ) = trans( 1, 1 ) * trans( 2, 2 ) -                        &
                        trans( 1, 2 ) * trans( 2, 1 )
!.......
        IF( DABS( angle ) .GT. 1.0D-8 ) THEN
          trans( 2, 1 ) = trans( 2, 1 ) * DCOS( angle ) +                      &
                          trans( 3, 1 ) * DSIN( angle )
          trans( 2, 2 ) = trans( 2, 2 ) * DCOS( angle ) +                      &
                          trans( 3, 2 ) * DSIN( angle )
          trans( 2, 3 ) = trans( 2, 3 ) * DCOS( angle ) +                      &
                          trans( 3, 3 ) * DSIN( angle )
          trans( 3, 1 ) = trans( 1, 2 ) * trans( 2, 3 ) -                      &
                          trans( 1, 3 ) * trans( 2, 2 )
          trans( 3, 2 ) = trans( 1, 3 ) * trans( 2, 1 ) -                      &
                          trans( 1, 1 ) * trans( 2, 3 )
          trans( 3, 3 ) = trans( 1, 1 ) * trans( 2, 2 ) -                      &
                          trans( 1, 2 ) * trans( 2, 1 )
        END IF
!.......
        DO imatr = 2, 4
          DO idofn = 1, 3
            DO jdofn = 1, 3
              iloca = ( imatr - 1 ) * 3 + idofn
              jloca = ( imatr - 1 ) * 3 + jdofn
              trans( iloca, jloca ) = trans( idofn, jdofn )
            END DO
          END DO
        END DO
        RETURN
        END

        SUBROUTINE ToGlobal( trans, iswth )
!........
!       模块功能
!           元素矩阵, 载荷矢量转换为整体坐标
!........
        USE CtrlData
        USE ElmtData
        IMPLICIT DOUBLE PRECISION( a-h, o-z )
        DIMENSION trans( 12, 12 ), vwork( 12, 12 )
        IF( iswth .NE. 6 ) THEN
!.........形成元素矩阵的对称部分
          DO idofn = 1, 12
            DO jdofn = 1, idofn
              stife( jdofn, idofn ) = stife( idofn, jdofn )
            END DO
          END DO
!.........元素矩阵转换为整体坐标矩阵
          DO idofn = 1, 12
            DO jdofn = 1, 12
              vwork( idofn, jdofn ) = 0.0D0
              DO kdofn = 1, 12
                vwork( idofn, jdofn ) = vwork( idofn, jdofn ) +                &
                stife( idofn, kdofn ) * trans( kdofn, jdofn )
              END DO
            END DO
          END DO
          IF( MOD( iswth, mswth ) .EQ. 7 ) THEN
            DO idofe = 1, 12
              DO jdofe = 1, 12
                stife( idofe, jdofe ) = vwork( idofe, jdofe )
              END DO
            END DO
          ELSE
            DO idofn = 1, 12
              DO jdofn = 1, 12
                stife( idofn, jdofn ) = 0.0D0
                DO kdofn = 1, 12
                  stife( idofn, jdofn ) = stife( idofn, jdofn ) +              &
                  trans( kdofn, idofn ) * vwork( kdofn, jdofn )
                END DO
              END DO
            END DO
          END IF
        ELSE IF( iswth .EQ. 6 ) THEN
!.........载荷矢量转换为整体坐标载荷矢量
          DO idofe = 1, 12
            vwork( idofe, 1 ) = 0.0D0
            DO jdofe = 1, 12
              vwork( idofe,     1 ) = vwork( idofe, 1 ) +                      &
              trans( jdofe, idofe ) * force( jdofe    )
            END DO
          END DO
          DO idofe = 1, 12
            force( idofe ) = vwork( idofe, 1 )
          END DO
        END IF

        RETURN
        END

        SUBROUTINE LinkMatrix004
        USE CtrlData
        USE ElmtData
        IMPLICIT DOUBLE PRECISION( a-h, o-z )
        DIMENSION swork( 24, 24 )

        DO idofe = 1, 24
          DO jdofe = 1, 24
            swork( idofe, jdofe ) = 0.0D0
          END DO
        END DO

        DO inode = 1, 2
          DO jnode = 1, 2
            DO idofn = 1, 6
              DO jdofn = 1, 6
                idofe = ( inode - 1 ) * 6 + idofn
                jdofe = ( jnode - 1 ) * 6 + jdofn
                irefn = inode + 2
                jrefn = jnode + 2
                iloca = idofe
                jloca = jdofe
                IF( mnode .GE. irefn .AND. lnode( irefn ) .GT. 0 .AND.         &
                    idofn .LE. 3 ) iloca = ( irefn - 1 ) * 6 + idofn
                IF( mnode .GE. jrefn .AND. lnode( jrefn ) .GT. 0 .AND.         &
                    jdofn .LE. 3 ) jloca = ( jrefn - 1 ) * 6 + jdofn
                swork( iloca, jloca ) = stife( idofe, jdofe )
              END DO
            END DO
          END DO
        END DO

        DO idofe = 1, MIN( ndofe, 24 )
          DO jdofe = 1, MIN( ndofe, 24 )
            stife( idofe, jdofe ) = swork( idofe, jdofe )
          END DO
        END DO

        IF( mnode .GT. 3 .AND. lnode( 3 ) .GT. 0 ) THEN
          DO idofn = 1, 3
            jdofn = ndofn * 2 + idofn
            stife( idofn, idofn ) = stife( jdofn, jdofn )
            stife( idofn, jdofn ) =-stife( jdofn, jdofn )
            stife( jdofn, idofn ) =-stife( jdofn, jdofn )
            stife( jdofn, jdofn ) = stife( jdofn, jdofn ) * 2.0D0
          END DO
        END IF

        IF( mnode .GT. 4 .AND. lnode( 4 ) .GT. 0 ) THEN
          DO idofn = 7, 9
            jdofn = ndofn * 2 + idofn
            stife( idofn, idofn ) = stife( jdofn, jdofn )
            stife( idofn, jdofn ) =-stife( jdofn, jdofn )
            stife( jdofn, idofn ) =-stife( jdofn, jdofn )
            stife( jdofn, jdofn ) = stife( jdofn, jdofn ) * 2.0D0
          END DO
        END IF
        END

        SUBROUTINE LinkForce004
        USE CtrlData
        USE ElmtData
        IMPLICIT DOUBLE PRECISION( a-h, o-z )

        IF( mnode .GE. 3 .AND. lnode( 3 ) .GT. 0 ) THEN
          DO idofn = 1, 3
            idofe = idofn + ndofn + ndofn
            force( idofe ) = force( idofn )
            force( idofn ) = 0.0D0
          END DO
        END IF

        IF( mnode .GE. 4 .AND. lnode( 4 ) .GT. 0 ) THEN
          DO idofn = ndofn + 1, ndofn + 3
            idofe = idofn + ndofn + ndofn
            force( idofe ) = force( idofn )
            force( idofn ) = 0.0D0
          END DO
        END IF
        END