eigenp.f90

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

        SUBROUTINE EigenP( toler, iflag )
!........
!        模块功能
!            逆幂法求解特征值问题.
!........
!       buffs 缓冲区数组(I)             shmod 特征向量(O)
!       diagk 刚度对角元(I)             diagm 质量对角元(I)
!       vslvs 载荷矢量(W)               statf 静力节点载荷(I)
!       vfixs 给定位移矢量(I)           lbcks 回代信息(I)
!       lfixs 位移约束信息(I)           lnods 单元定义(I)
!       lopts 优化单元次序(I)           lpnte 节点单元数(W)
!       lpntw 波阵节点信息(W)           lrcrd 记录信息(I)
!       lsavc 存储信息(I)               toler 迭代精度(I)
!       mdofn 最大节点自由度(I)         mdofw 最大波阵自由度(I)
!       mnode 最大单元节点数(I)         mpntw 最大波阵节点数(I)
!       nbufs 缓冲区大小(I)             ndime 问题维数(I)
!       ndofs 问题总自由度数(I)         nelem 单元总数(I)
!       nerrc 出错控制参数(I/O)         nexch 交换缓冲区大小(I)
!       iflag 特征向量输出控制          nmode 特征对数(I)
!         == 2 质量阵乘质量阵规一化     npoin 节点总数(I)
!         == 1 质量阵规一化             nrcrd 一致矩阵记录数(I)
!         == 0 规一化
!........
        USE CtrlData
        USE MeshData
        USE GlobData
        USE SolvData
        USE FrontData
        IMPLICIT DOUBLE PRECISION( a-h, o-z )
        INTEGER * 2 dummy, rseed, secnd
!.......求解特征值问题.
        OPEN( 26, FORM = 'UNFORMATTED', STATUS = 'SCRATCH' )
        DO imode = 1, nmode
          label = 0
          rseed = 0
          vxtkx = 0.0D0
          eigen = 1.0D8
          DO WHILE( rseed .EQ. 0 )
            CALL GETTIM( dummy, dummy, secnd, rseed )
          END DO
          vrand = rseed + secnd * 13.0D0
          DO WHILE( vrand .GT. 1.0D0 )
            vrand = vrand / ( 4.0D0 * DATAN( 1.0D0 ) )
          END DO
          DO idofs = 1, ndofs
            vtemp = vrand * vrand
            DO WHILE( vtemp .LT. 1.0D0 )
              vtemp = vtemp * 10.D0
            END DO
            rseed = vtemp + 0.5D0
            vrand = vtemp - rseed
            shmod( idofs, imode ) = vrand
          END DO
          DO WHILE( label .EQ. 0 )
            iswth = 2
            DO idofs = 1, ndofs
              vslvs( idofs ) = 0.0D0
            END DO
            CALL DotMat( shmod, vslvs, iswth, imode, nmode )
            IF( nerrc .GT. 0 ) RETURN
            CALL NormVec( shmod, vslvs, freqs, vxtkx, imode, ndofs,            &
                          nmode, eigen, toler, label )
            IF( label .EQ. 0 ) THEN
              iswth = 1
              DO idofs = 1, ndofs
                shmod( idofs, imode ) = vslvs( idofs )
              END DO
              CALL BckSub( 1 )
              vxtkx = 0.0D0
              DO idofs = 1, ndofs
                vxtkx = vxtkx + shmod( idofs, imode ) * vslvs( idofs )
                shmod( idofs, imode ) = vslvs( idofs )
              END DO
            END IF
          END DO

          imatx = 3
          vxtcx = 0.0D0
          IF( lsavc( 3 ) .NE. 0 ) THEN
            DO idofs = 1, ndofs
              fwork( idofs ) = 0.0D0
            END DO
            CALL DotMat( shmod, fwork, imatx, imode, nmode )
            DO idofs = 1, ndofs
              vxtcx = vxtcx + fwork( idofs ) * shmod( idofs, imode )
            END DO
          END IF

          imatx = 2
          vxtmx = 0.0D0
          DO idofs = 1, ndofs
            fwork( idofs ) = 0.0D0
          END DO
          CALL DotMat( shmod, fwork, imatx, imode, nmode )
          DO idofs = 1, ndofs
            vxtmx = vxtmx + fwork( idofs ) * shmod( idofs, imode )
          END DO

          vflag = 1.0D0
          vdamp = 0.0D0
          eigen = eigen - shift( 2 )
          IF( eigen .LT. 1.0D-8 ) vflag =-1.0D0
          omiga = vflag * DSQRT( DABS( eigen ) )
          vfreq = omiga / ( 8.0D0 * DATAN( 1.0D0 ) )
          IF( DABS( vfreq ) .GT. 1.0D-4 )                                      &
          vdamp = vxtcx / ( 2.0D0 * vxtmx * omiga )
          IF( ntypc .EQ. 3 ) WRITE( 12, 1800 ) imode, eigen, vdamp
          IF( ntypc .NE. 3 ) WRITE( 12, 2000 ) imode, vfreq, vdamp
          freqs( imode ) = omiga

          DO ipoin = 1, npoin
            jpoin = llnks( ipoin )
            idofs = ( jpoin - 1 ) * mdofn + 1
            jdofs =   jpoin       * mdofn
            WRITE( 12, 2200 ) ipoin, ( shmod( kdofs, imode ),                  &
                                      kdofs = idofs, jdofs )
          END DO
          vnorm = 0.0D0
          DO idofs = 1, ndofs
            vnorm = vnorm + shmod( idofs, imode ) * vslvs( idofs )
          END DO
          vnorm = 1.0D0 / DSQRT( DABS( vnorm ) )
          IF( iflag .EQ. 1 ) THEN
            DO idofs = 1, ndofs
              shmod( idofs, imode ) = shmod( idofs, imode ) * vnorm
            END DO
          END IF
          IF( iflag .NE. 2 ) THEN
            WRITE( 26 ) ( shmod( idofs, imode ), idofs = 1, ndofs )
          END IF
          DO idofs = 1, ndofs
            shmod( idofs, imode ) = vslvs( idofs ) * vnorm
          END DO
        END DO
        IF( iflag .NE. 2 ) THEN
          REWIND( 26 )
          DO imode = 1, nmode
            READ( 26 ) ( shmod( idofs, imode ), idofs = 1, ndofs )
          END DO
        END IF
        CLOSE( 26 )
        DO ipoin = 1, npoin
          jpoin = llnks( ipoin )
          IF( ipoin .NE. jpoin ) THEN
            DO idofn = 1, mdofn
              idofs = ( ipoin - 1 ) * mdofn + idofn
              jdofs = ( jpoin - 1 ) * mdofn + idofn
              DO imode = 1, nmode
                shmod( idofs, imode ) = shmod( jdofs, imode )
              END DO
            END DO
          END IF
        END DO
1800    FORMAT( //20x, '第', I3, '阶失稳模态',                                 &
                 /10x, '载荷因子', E15.6, 5x, '阻尼比', E15.6,                 &
                //20x, '失  稳  模  态' )
2000    FORMAT( //20x, '第', I3, '阶固有模态',                                 &
                 /10x, '频率', E15.6, '赫兹(Hz)', 5x, '阻尼比', E15.6,         &
                //20x, '振          型' )
2200    FORMAT( 5X, I5, 4E15.5 / ( 10X, 4E15.5 ) )
        RETURN
        END
        SUBROUTINE NormVec( shmod, vslvs, freqs, vxtkx, kmode, ndofs,          &
                            nmode, eigen, toler, label )
!........
!       模块功能
!           迭代特征矢量规范化, 正交化, 控制迭代精度.
!........
!       shmod 特征向量数组(I/O)   vslvs 节点载荷矢量(W)
!       kmode 当前特征对号(I)     ndofs 自由度总数目(I)
!       nmode 特征对总数目(I)     eigen 当前特征对号(O)
!       toler 迭代计算精度(I)     label 收敛控制参数(O)
!........
        IMPLICIT DOUBLE PRECISION( a-h, o-z )
        DIMENSION freqs( nmode )
        DIMENSION vslvs( ndofs ), shmod( ndofs, nmode )
        vxtmx = 0.0D0
        DO idofs = 1, ndofs
          vxtmx = vxtmx + shmod( idofs, kmode ) * vslvs( idofs )
        END DO
!.......迭代特征矢量正交化.
        DO imode = 1, kmode - 1
          vnorm = 0.0D0
          DO idofs = 1, ndofs
            vnorm = vnorm + shmod( idofs, imode ) *                            &
                            shmod( idofs, kmode )
          END DO
          IF( freqs( imode ) .LT.-1.0D-5 ) vnorm =-vnorm
          DO idofs = 1, ndofs
            vslvs( idofs ) = vslvs( idofs ) -                                  &
                     vnorm * shmod( idofs, imode )
          END DO
        END DO
!.......迭代特征矢量规范化.
        vmaxi = 0.0D0
        DO idofs = 1, ndofs
          IF( DABS( shmod( idofs, kmode ) ) .GT. DABS( vmaxi ) )               &
          vmaxi = shmod( idofs, kmode )
        END DO
        vmaxi = 1.0D0 / vmaxi
        DO idofs = 1, ndofs
          vslvs( idofs        ) = vslvs( idofs        ) * vmaxi
          shmod( idofs, kmode ) = shmod( idofs, kmode ) * vmaxi
        END DO
!.......计算固有频率.
        vmaxi = vxtkx / vxtmx
        IF( DABS( vmaxi ) .LT. 1.0D-8 ) THEN
          tolec = DABS( vmaxi - eigen )
        ELSE
          tolec = DSQRT( DABS( ( vmaxi - eigen ) / vmaxi ) )
        END IF
        IF( tolec .LT. toler ) label = 1
        eigen = vmaxi

        RETURN
        END