📄 elmt002.f90
字号:
SUBROUTINE Elmt002( ielem, imats, iswth, iuses )
!........
! 模块功能
! 一维线性变截面轴单元
! iswth == 1 单元刚度矩阵
! == 2 单元一致质量矩阵
! == 3 单元阻尼矩阵
! == 4 单元几何矩阵(临界转速)
! == 5 单元内力等效节点力
! == 6 单元分布载荷等效节点力
! == 7 单元内力
!........
USE CtrlData
USE ElmtData
USE MeshData
USE FactorData
IMPLICIT DOUBLE PRECISION( a-h, o-z )
DOUBLE PRECISION lmatx( 4, 4 ), rmatx( 4, 4 )
DOUBLE PRECISION wmatx( 4, 4 ), hmatx( 4, 4 )
IF( iswth .EQ. 0 ) THEN
READ( 11, 2000 ) elast, denst, omiga, douti, dinni, &
doutj, dinnj, width, itype, ibeam, &
nlami, thick, const, alpha
IF( ibeam .NE. 0 ) ibeam = 1
IF( nprnc .NE. 0 ) THEN
SELECT CASE( itype )
CASE( 0 )
WRITE( 12, 2200 )
IF( ibeam .EQ. 0 ) WRITE( 12, 2400 )
IF( ibeam .NE. 0 ) WRITE( 12, 2401 )
WRITE( 12, 2600 ) elast, denst, omiga, douti, &
dinni, doutj, dinnj
CASE( 1 )
WRITE( 12, 2201 )
IF( ibeam .EQ. 0 ) WRITE( 12, 2400 )
IF( ibeam .NE. 0 ) WRITE( 12, 2401 )
WRITE( 12, 2601 ) elast, denst, omiga, douti, &
dinni, doutj, dinnj
CASE( 2 )
WRITE( 12, 2202 )
IF( ibeam .EQ. 0 ) WRITE( 12, 2400 )
IF( ibeam .NE. 0 ) WRITE( 12, 2401 )
WRITE( 12, 2602 ) denst, omiga, douti, dinni, &
doutj, dinnj, nlami, thick
CASE( 3 )
WRITE( 12, 2203 )
IF( ibeam .EQ. 0 ) WRITE( 12, 2400 )
IF( ibeam .NE. 0 ) WRITE( 12, 2401 )
WRITE( 12, 2603 ) denst, omiga, douti, dinni, &
doutj, dinnj, nlami, thick, &
const, alpha
CASE( 4 )
WRITE( 12, 2204 )
IF( ibeam .EQ. 0 ) WRITE( 12, 2400 )
IF( ibeam .NE. 0 ) WRITE( 12, 2401 )
WRITE( 12, 2604 ) denst, omiga, douti, &
dinni, doutj, dinnj, width
CASE DEFAULT
END SELECT
END IF
props( 2, imats ) = elast
props( 3, imats ) = denst
props( 4, imats ) = omiga
props( 5, imats ) = douti
props( 6, imats ) = dinni
props( 7, imats ) = doutj
props( 8, imats ) = dinnj
props( 9, imats ) = itype + 0.5D0
props( 10, imats ) = ibeam + 0.5D0
props( 11, imats ) = nlami + 0.5D0
props( 12, imats ) = thick
props( 13, imats ) = const
props( 14, imats ) = alpha
props( 15, imats ) = width
RETURN
ELSE IF( iswth .EQ. 8 ) THEN
nnode = 2
ndofn = 2
RETURN
ELSE IF( iswth .EQ. 9 ) THEN
IF( iuses .EQ. 0 ) RETURN
WRITE( 15, 2800 ) 2, ielem, imats, 2, lnode( 1 ), lnode( 2 )
RETURN
ELSE IF( iswth .EQ. 21 ) THEN
WRITE( 25 ) 0
RETURN
END IF
DO idofe = 1, ndofe
DO jdofe = 1, ndofe
stife( idofe, jdofe ) = 0.0D0
END DO
force( idofe ) = 0.0D0
END DO
IF( iuses .EQ. 0 ) RETURN
omiga = props( 4, imats )
itype = props( 9, imats )
IF( coren( 1, 2 ) - coren( 1, 1 ) .LT. 0.0D0 ) THEN
WRITE( 12, 3000 ) ielem
nerrc = 0020001
RETURN
END IF
IF( iswth .EQ. 1 .OR. iswth .EQ. 5 .OR. &
iswth .EQ. 7 .OR. iswth .EQ. 10 ) THEN
IF( itype .GT. 1 ) RETURN
CALL Lmat002( lmatx, imats )
CALL Hmat002( hmatx, imats )
IF( omiga .GE. 1.0D-8 ) THEN
CALL Rmat002( rmatx, imats )
DO idime = 1, 4
DO jdime = 1, 4
hmatx( idime, jdime ) = hmatx( idime, jdime ) - &
rmatx( idime, jdime ) * omiga * omiga
END DO
END DO
END IF
CALL MultMtr( hmatx, lmatx, wmatx, 0, 4 )
IF( nerrc .GT. 0 ) RETURN
CALL MultMtr( lmatx, wmatx, stife, 1, ndofe )
IF( nerrc .GT. 0 ) RETURN
IF( nincc .EQ. 2 .OR. iswth .EQ. 10 .OR. iswth .EQ. 7 ) THEN
DO idofe = 1, 4
force( idofe ) = 0.0D0
DO jdofe = 1, 4
force( idofe ) = force( idofe ) - &
dispe( jdofe ) * stife( idofe, jdofe )
END DO
END DO
IF( iswth .EQ. 7 ) THEN
WRITE( 12, 3200 ) ielem
WRITE( 12, 3400 ) (-force( idofe ), idofe = 1, 4 )
END IF
END IF
ELSE IF( iswth .EQ. 2 ) THEN
ELSE IF( iswth .EQ. 3 ) THEN
ELSE IF( iswth .EQ. 4 ) THEN
CALL Lmat002( lmatx, imats )
CALL Rmat002( rmatx, imats )
CALL MultMtr( rmatx, lmatx, wmatx, 0, 4 )
IF( nerrc .GT. 0 ) RETURN
CALL MultMtr( lmatx, wmatx, stife, 1, mdofe )
IF( nerrc .GT. 0 ) RETURN
IF( itype .EQ. 4 ) THEN
DO idofe = 1, 4
DO jdofe = 1, 4
rmatx( idofe, jdofe ) = stife( idofe, jdofe )
END DO
END DO
CALL Trans002( lmatx, props( 15, imats ) )
CALL MultMtr( rmatx, lmatx, wmatx, 0, 4 )
IF( nerrc .GT. 0 ) RETURN
CALL MultMtr( lmatx, wmatx, stife, 1, mdofe )
IF( nerrc .GT. 0 ) RETURN
END IF
ELSE IF( iswth .EQ. 6 ) THEN
END IF
2000 FORMAT( 4F15.5 / 4F15.5 / 3I5, 3F15.5 )
2200 FORMAT( 2x, '单元类型: 转轴单元(轴)' )
2201 FORMAT( 2x, '单元类型: 转轴单元(轮)' )
2202 FORMAT( 2x, '单元类型: 转轴单元(直叶片)' )
2203 FORMAT( 2x, '单元类型: 转轴单元(曲叶片)' )
2204 FORMAT( 2x, '单元类型: 转轴单元(外伸构件)' )
2400 FORMAT( 2x, ' (按轴分析)' )
2401 FORMAT( 2x, ' (按梁分析)' )
2600 FORMAT( 2x, '杨氏模量: ', E15.5/2x, '质量密度：', E15.5 / &
2x, '角 速 度: ', E15.5/2x, '节点外径：', F15.5 / &
2x, '节点内径: ', F15.5/2x, '节点外径：', F15.5 / &
2x, '节点内径: ', F15.5 )
2601 FORMAT( 2x, '杨氏模量: ', E15.5/2x, '质量密度：', E15.5 / &
2x, '角 速 度: ', E15.5/2x, '节点外径：', F15.5 / &
2x, '节点内径: ', F15.5/2x, '节点外径：', F15.5 / &
2x, '节点内径: ', F15.5 )
2602 FORMAT( 2x, '质量密度：', E15.5 /2x, '角 速 度: ', E15.5/ &
2x, '节点外径：', F15.5 /2x, '节点内径: ', F15.5/ &
2x, '节点外径：', F15.5 /2x, '节点内径: ', F15.5/ &
2x, '叶片数目：', I15 /2x, '叶片厚度：', F15.5 )
2603 FORMAT( 2x, '质量密度：', E15.5 /2x, '角 速 度: ', E15.5/ &
2x, '节点外径：', F15.5 /2x, '节点内径: ', F15.5/ &
2x, '节点外径：', F15.5 /2x, '节点内径: ', F15.5/ &
2x, '叶片数目：', I15 /2x, '叶片厚度：', F15.5/ &
2x, '常 数 c：', I15 /2x, '常 数 α：', F15.5 )
2604 FORMAT( 2x, '质量密度：', E15.5 /2x, '角 速 度: ', E15.5/ &
2x, '节点外径：', F15.5 /2x, '节点内径: ', F15.5/ &
2x, '节点外径：', F15.5 /2x, '节点内径: ', F15.5/ &
2x, '宽 度：', F15.5 )
2800 FORMAT( 10I5 )
3000 FORMAT( ///2x, '单元', I5, '没有按着自左而右的定义原则！'// )
3200 FORMAT( //10x, '单元', I5, ' 内力' )
3400 FORMAT( 2x, 4E15.5 )
RETURN
END
SUBROUTINE Lmat002( lmatx, imats )
USE ElmtData
USE MeshData
IMPLICIT DOUBLE PRECISION( a-h, o-z )
DOUBLE PRECISION lmatx( 4, 4 )
itype = props( 9, imats )
IF( itype .NE. 4 ) const = DABS( coren( 1, 2 ) - coren( 1, 1 ) )
IF( itype .EQ. 4 ) const = DABS( props( 15, imats ) )
lmatx( 1, 1 ) =-2.0D0
lmatx( 1, 2 ) =-const
lmatx( 1, 3 ) = 2.0D0
lmatx( 1, 4 ) =-const
lmatx( 2, 1 ) = 3.0D0
lmatx( 2, 2 ) = const
lmatx( 2, 3 ) =-3.0D0
lmatx( 2, 4 ) = const + const
lmatx( 3, 1 ) = 0.0D0
lmatx( 3, 2 ) = 0.0D0
lmatx( 3, 3 ) = 0.0D0
lmatx( 3, 4 ) =-const
lmatx( 4, 1 ) = 0.0D0
lmatx( 4, 2 ) = 0.0D0
lmatx( 4, 3 ) = 1.0D0
lmatx( 4, 4 ) = 0.0D0
RETURN
END
SUBROUTINE Rmat002( rmatx, imats )
USE CtrlData
USE ElmtData
USE MeshData
IMPLICIT DOUBLE PRECISION( a-h, o-z )
DIMENSION twork( 4, 4, 2 ), rmatx( 4, 4 )
DIMENSION swork( 4, 4, 2 ), cmatx( 2, 4 )
DIMENSION cwork( 4, 4, 2 ), rwork( 4, 4, 2 )
denst = props( 3, imats )
douti = props( 5, imats )
dinni = props( 6, imats )
doutj = props( 7, imats )
dinnj = props( 8, imats )
itype = props( 9, imats )
ibeam = props( 10, imats )
IF( ibeam .EQ. 0 ) ibeam = -1
IF( itype .EQ. 4 ) eleng = DABS( props( 15, imats ) )
IF( itype .NE. 4 ) eleng = DABS( coren( 1, 2 ) - coren( 1, 1 ) )
IF( itype .EQ. 2 ) THEN
nlami = props( 11, imats )
thick = props( 12, imats )
const = nlami * thick * denst * eleng / 2.0D0
DO idime = 1, 4
DO jdime = 1, 4
index = 8 - idime - jdime
rmatx( idime, jdime ) = ( F1( douti, doutj, index ) - &
F1( dinni, dinnj, index ) ) * &
const
END DO
END DO
const = const / ( 24.0D0 * eleng * eleng ) * ibeam
DO idime = 1, 3
DO jdime = 1, 3
index = 6 - idime - jdime
ipara = ( 4 - idime ) * ( 4 - jdime )
rmatx( idime, jdime ) = rmatx( idime, jdime ) + ( &
F3( douti, doutj, index ) - &
F3( dinni, dinnj, index ) ) * &
const * DFLOAT( ipara )
END DO
END DO
ELSE IF( itype .EQ. 3 ) THEN
DO idime = 1, 4
DO jdime = 1, 4
DO kdime = 1, 2
twork( idime, jdime, kdime ) = 0.0D0
swork( idime, jdime, kdime ) = 0.0D0
cwork( idime, jdime, kdime ) = 0.0D0
rwork( idime, jdime, kdime ) = 0.0D0
END DO
rmatx( idime, jdime ) = 0.0D0
END DO
END DO
iflag = 1
kcoun = 1
nitem = 1
halfs = 0.5D0
CALL C2002( cmatx, 0.0D0, eleng )
CALL CurveLamina002( 0.0D0, vbend, vtors, imats )
CALL CurveR002( rmatx, cmatx, vbend, vtors, halfs )
CALL C2002( cmatx, 1.0D0, eleng )
CALL CurveLamina002( 1.0D0, vbend, vtors, imats )
CALL CurveR002( rmatx, cmatx, vbend, vtors, halfs )
DO idime = 1, 4
DO jdime = 1, 4
twork( idime, jdime, 2 ) = rmatx( idime, jdime )
rmatx( idime, jdime ) = rmatx( idime, jdime ) * 0.5D0
END DO
END DO
DO WHILE( iflag .GT. 0 )
DO idime = 1, 4
DO jdime = 1, 4
twork( idime, jdime, 1 ) = twork( idime, jdime, 2 )
swork( idime, jdime, 1 ) = swork( idime, jdime, 2 )
cwork( idime, jdime, 1 ) = cwork( idime, jdime, 2 )
rwork( idime, jdime, 1 ) = rwork( idime, jdime, 2 )
END DO
END DO
DO iitem = 1, nitem
xloca = ( 2.0D0 * iitem - 1.0D0 ) * halfs
CALL C2002( cmatx, xloca, eleng )
CALL CurveLamina002( xloca, vbend, vtors, imats )
CALL CurveR002( rmatx, cmatx, vbend, vtors, halfs )
END DO
DO idime = 1, 4
DO jdime = 1, 4
twork( idime, jdime, 2 ) = rmatx( idime, jdime )
rmatx( idime, jdime ) = rmatx( idime, jdime ) * 0.5D0
END DO
END DO
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -