📄 dynamic.f90
字号:
SUBROUTINE Initial
!........
! 模块功能
! 读入初时条件(初时位移和初时速度).
!........
USE CtrlData
USE GlobData
USE ElmtData
IMPLICIT DOUBLE PRECISION( a-h, o-z )
IF( ntypc .EQ. 2 ) THEN
!.........初时化初时位移和初时速度
DO iload = 1, nload
DO idofs = 1, ndofs
disps( idofs, iload ) = 0.0D0
veloc( idofs, iload ) = 0.0D0
END DO
END DO
DO iload = 1, nload
!...........按各工况读入初时条件, 读入初时位移
label = 0
DO WHILE( label .EQ. 0 )
READ( 14, 2000 ) ipoin, idofn, value
IF( ipoin .GT. 0 .AND. ipoin .LE. npoin ) THEN
IF( idofn .GE. 1 .AND. idofn .LE. mdofn ) THEN
idofs = ( ipoin - 1 ) * mdofn + idofn
disps( idofs, iload ) = value
ELSE
nerrc = 1236
RETURN
END IF
ELSE
label = 1
END IF
END DO
!...........按各工况读入初时条件, 读入初时速度.
label = 0
DO WHILE( label .EQ. 0 )
READ( 14, 2000 ) ipoin, idofn, value
IF( ipoin .GT. 0 .AND. ipoin .LE. npoin ) THEN
IF( idofn .GE. 1 .AND. idofn .LE. mdofn ) THEN
idofs = ( ipoin - 1 ) * mdofn + idofn
veloc( idofs, iload ) = value
ELSE
nerrc = 1236
RETURN
END IF
ELSE
label = 1
END IF
END DO
END DO
!.........输出初时条件.
IF( nprnc .NE. 0 ) THEN
WRITE( 12, 2200 )
DO iload = 1, nload
WRITE( 12, 2400 ) iload
DO ipoin = 1, npoin
idofs = ( ipoin - 1 ) * mdofn + 1
jdofs = ipoin * mdofn
WRITE( 12, 2600 ) ipoin, ( disps( kdofs, iload ), &
kdofs = idofs, jdofs )
END DO
WRITE( 12, 2800 )
DO ipoin = 1, npoin
idofs = ( ipoin - 1 ) * mdofn + 1
jdofs = ipoin * mdofn
WRITE( 12, 2600 ) ipoin, ( veloc( kdofs, iload ), &
kdofs = idofs, jdofs )
END DO
END DO
END IF
END IF
IF( nhstr .GT. 0 ) THEN
label = 1
DO WHILE( label .EQ. 1 )
READ( 14, 3000 ) ielem, ( histe( ihstr ), ihstr = 1, nhstr )
IF( ielem .LE. 0 .OR. ielem .GT. nelem ) THEN
label = 0
ELSE
DO ihstr = 1, nhstr
histr( ihstr, ielem ) = histe( ihstr )
END DO
END IF
END DO
END IF
2000 FORMAT( 2I5, F15.5 )
2200 FORMAT( //20x, '= = = = 动力响应初时条件 = = = = ' )
2400 FORMAT( //20x, ' 工况(', I2, ')' &
/20X, ' 初时位移' )
2600 FORMAT( I5, 4E15.5 / ( 5X, 4E15.5 ) )
2800 FORMAT( /20X, ' 初时速度' )
3000 FORMAT( I5, 4E15.5, 10( /5x, 4E15.5 ) )
RETURN
END
SUBROUTINE DynFrc
USE CtrlData
USE MeshData
USE SolvData
USE GlobData
USE FrontData
IMPLICIT DOUBLE PRECISION( a-h, o-z )
IF( nerrc .GT. 0 ) RETURN
IF( ndync .GE. 1 .AND. ndync .LE. 3 ) THEN
IF( lsavc( 2 ) .EQ. 0 .OR. ldecc( 2 ) .NE. 1 ) THEN
nerrc = 2455
RETURN
END IF
DO idofs = 1, ndofs
vslvs( idofs ) = forcs( idofs ) + shift( 10 ) * ( &
forcs( idofs ) - frctm( idofs ) )
END DO
IF( nincc .LE. 1 ) THEN
DO idofs = 1, ndofs
fwork( idofs ) = shift( 4 ) * distm( idofs, kload ) + &
shift( 5 ) * veltm( idofs, kload ) + &
shift( 6 ) * acctm( idofs, kload )
END DO
ELSE IF( nincc .EQ. 2 ) THEN
DO idofs = 1, ndofs
fwork( idofs ) = shift( 4 ) * distm( idofs, kload ) + &
shift( 5 ) * veltm( idofs, kload ) + &
shift( 6 ) * acctm( idofs, kload ) - &
shift( 2 ) * vslvi( idofs )
END DO
ELSE IF( nincc .EQ. 3 ) THEN
DO idofs = 1, ndofs
fwork( idofs ) = shift( 4 ) * distm( idofs, kload ) + &
shift( 5 ) * veltm( idofs, kload ) + &
shift( 6 ) * acctm( idofs, kload ) - &
shift( 2 ) * distm( idofs, kload ) + &
vslvt( idofs )
END DO
END IF
CALL DotMat( fwork, vslvs, 2, 1, 1 )
IF( lsavc( 3 ) .NE. 0 ) THEN
IF( nincc .LE. 1 ) THEN
DO idofs = 1, ndofs
fwork( idofs ) = shift( 7 ) * distm( idofs, kload ) + &
shift( 8 ) * veltm( idofs, kload ) + &
shift( 9 ) * acctm( idofs, kload )
END DO
ELSE IF( nincc .EQ. 2 ) THEN
DO idofs = 1, ndofs
fwork( idofs ) = shift( 7 ) * distm( idofs, kload ) + &
shift( 8 ) * veltm( idofs, kload ) + &
shift( 9 ) * acctm( idofs, kload ) - &
shift( 3 ) * disps( idofs, kload )
END DO
ELSE IF( nincc .EQ. 3 ) THEN
DO idofs = 1, ndofs
fwork( idofs ) = shift( 7 ) * distm( idofs, kload ) + &
shift( 8 ) * veltm( idofs, kload ) + &
shift( 9 ) * acctm( idofs, kload ) - &
shift( 3 ) * distm( idofs, kload ) + &
veltm( idofs, kload )
END DO
END IF
CALL DotMat( fwork, vslvs, 3, 1, 1 )
END IF
ELSE
DO idofs = 1, ndofs
vslvs( idofs ) = forcs( idofs )
END DO
END IF
IF( nincc .EQ. 3 ) THEN
DO idofs = 1, ndofs
vslvs( idofs ) = vslvs( idofs ) - frctm( idofs )
END DO
END IF
RETURN
END
SUBROUTINE DynSlv
!........
!
! ndync == 1 Wilson
! == 2 Newmark
! == 3 Crank
!........
USE CtrlData
USE GlobData
IMPLICIT DOUBLE PRECISION( a-h, o-z )
IF( ndync .EQ. 1 ) THEN
cnst1 = 6.0D0 / ( theta * theta * theta * tstep * tstep )
cnst2 = 6.0D0 / ( theta * theta * tstep )
cnst3 = 1.0D0 - 3.0D0 / theta
cnst4 = tstep / 2.0D0
cnst5 = tstep * tstep / 6.0D0
DO idofs = 1, ndofs
accel( idofs, kload ) = cnst1 * disps( idofs, kload ) - &
cnst1 * distm( idofs, kload ) - &
cnst2 * veltm( idofs, kload ) + &
cnst3 * acctm( idofs, kload )
veloc( idofs, kload ) = veltm( idofs, kload ) + cnst4 * ( &
accel( idofs, kload ) + acctm( idofs, kload ) )
disps( idofs, kload ) = distm( idofs, kload ) + tstep * &
veltm( idofs, kload ) + cnst5 * ( &
accel( idofs, kload ) + 2.0D0 * acctm( idofs, kload ) )
END DO
ELSE IF( ndync .EQ. 2 ) THEN
cnst1 = 1.0D0 / ( belta * tstep * tstep )
cnst2 = 1.0D0 / ( belta * tstep )
cnst3 = 0.5D0 / belta - 1.0D0
cnst4 = alpha / ( belta * tstep )
cnst5 = alpha / belta - 1.0D0
cnst6 = tstep * ( 0.5D0 * alpha / belta - 1.0D0 )
DO idofs = 1, ndofs
accel( idofs, kload ) = cnst1 * disps( idofs, kload ) - &
cnst1 * distm( idofs, kload ) - &
cnst2 * veltm( idofs, kload ) - &
cnst3 * acctm( idofs, kload )
veloc( idofs, kload ) = cnst4 * disps( idofs, kload ) - &
cnst4 * distm( idofs, kload ) - &
cnst5 * veltm( idofs, kload ) - &
cnst6 * acctm( idofs, kload )
END DO
ELSE IF( ndync .EQ. 3 ) THEN
DO idofs = 1, ndofs
disps( idofs, kload ) = disps( idofs, kload ) * 2.0D0 - &
distm( idofs, kload )
END DO
ELSE
nerrc = 23234
END IF
RETURN
END
SUBROUTINE Update
USE CtrlData
USE MeshData
USE GlobData
USE SolvData
IMPLICIT DOUBLE PRECISION( a-h, o-z )
IF( nupdc .EQ. 1 ) nerrc = 23435
IF( nerrc .GT. 0 ) RETURN
nupdc = 1
timec = timec + tstep
DO idofs = 1, ndofs
fstrt( idofs ) = fstrs( idofs )
frctm( idofs ) = forcs( idofs )
distm( idofs, kload ) = disps( idofs, kload )
IF( ntypc .EQ. 2 ) THEN
veltm( idofs, kload ) = veloc( idofs, kload )
acctm( idofs, kload ) = accel( idofs, kload )
END IF
END DO
RETURN
END
SUBROUTINE Newmark
USE CtrlData
USE SolvData
USE FrontData
IMPLICIT DOUBLE PRECISION( a-h, o-z )
IF( ntypc .NE. 2 ) RETURN
shift( 1 ) = 1.0D0
shift( 2 ) = 1.0D0 / ( belta * tstep * tstep )
shift( 3 ) = alpha / ( belta * tstep )
shift( 4 ) = shift( 2 )
shift( 5 ) = 1.0D0 / ( belta * tstep )
shift( 6 ) = 0.5D0 / belta - 1.0D0
shift( 7 ) = shift( 3 )
shift( 8 ) = alpha / belta - 1.0D0
shift( 9 ) = 0.5D0 * alpha / belta - 1.0D0
shift( 10 ) = 0.0D0
DO imatx = 1, 3
IF( lsavc( imatx ) .NE. 0 .AND. ldecc( imatx ) .NE. 1 ) THEN
nerrc = 2334
RETURN
END IF
END DO
nshfc = 1
CALL ShftMat
RETURN
END
SUBROUTINE Wilson
USE CtrlData
USE SolvData
USE FrontData
IMPLICIT DOUBLE PRECISION( a-h, o-z )
IF( ntypc .NE. 2 ) RETURN
shift( 1 ) = 1.0D0
shift( 2 ) = 6.0D0 / ( theta * tstep ) ** 2
shift( 3 ) = 3.0D0 / ( theta * tstep )
shift( 4 ) = shift( 2 )
shift( 5 ) = 2.0D0 * shift( 3 )
shift( 6 ) = 2.0D0
shift( 7 ) = shift( 3 )
shift( 8 ) = 2.0D0
shift( 9 ) = theta * tstep / 2.0D0
shift( 10 ) = theta
DO imatx = 1, 3
IF( lsavc( imatx ) .NE. 0 .AND. &
ldecc( imatx ) .NE. 1 ) THEN
nerrc = 2334
RETURN
END IF
END DO
nshfc = 1
CALL ShftMat
RETURN
END
SUBROUTINE Crank
USE CtrlData
USE SolvData
USE FrontData
IMPLICIT DOUBLE PRECISION( a-h, o-z )
IF( ntypc .NE. 2 ) RETURN
shift( 1 ) = 1.0D0
shift( 2 ) = 2.0D0 / tstep
shift( 3 ) = 0.0D0
shift( 4 ) = 2.0D0 / tstep
shift( 5 ) = 0.0D0
shift( 6 ) = 0.0D0
shift( 7 ) = 0.0D0
shift( 8 ) = 0.0D0
shift( 9 ) = 0.0D0
shift( 10 ) =-0.5D0
DO imatx = 1, 3
IF( lsavc( imatx ) .NE. 0 .AND. ldecc( imatx ) .NE. 1 ) THEN
nerrc = 2334
RETURN
END IF
END DO
nshfc = 1
CALL ShftMat
RETURN
END
SUBROUTINE OptStep( kdofn, nstep )
!........
! 模块功能
! 根据上一步的计算结果选择下一步时间步长. 应当注意的是, 宏命
! 令 OPTSTEP必须在某一时间步计算完成但还没有执行UPDATE时执行, 否
! 则, 可能导致致命性错误.
!........
USE CtrlData
USE MeshData
USE GlobData
USE SolvData
USE FrontData
IMPLICIT DOUBLE PRECISION( a-h, o-z )
IF( ldecc( 2 ) .NE. 1 ) nerrc = 1224
IF( nupdc .EQ. 1 ) nerrc = 3243
IF( nerrc .GT. 0 ) RETURN
DO idofs = 1, ndofs
forcs( idofs ) = 0.0D0
fwork( idofs ) = disps( idofs, kload ) - &
distm( idofs, kload )
END DO
CALL DotMat( fwork, forcs, 2, 1, 1 )
upper = 0.0D0
downr = 0.0D0
DO ipoin = 1, npoin
DO idofn = 1, mdofn
IF( idofn .EQ. kdofn .OR. kdofn .LE. 0 ) THEN
idofs = ( ipoin - 1 ) * mdofn + idofn
fstep = fstrt( idofs ) - fstrs( idofs )
downr = downr + forcs( idofs ) * fwork( idofs )
upper = upper + fwork( idofs ) * fstep
END IF
END DO
END DO
omiga = DSQRT( upper / downr )
IF( nstep .LE. 0 ) nstep = 50
tstep = 8.0D0 * DATAN( 1.0D0 ) / ( omiga * nstep )
RETURN
END
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -