📄 elmt004.f90
字号:
SUBROUTINE Elmt004( ielem, imats, iswth, iuses )
!........
! 模块功能
! 空间梁单元
! iswth == 1 单元刚度矩阵.
! == 2 单元一致质量阵.
! == 3 单元阻尼矩阵.
! == 4 单元几何矩阵.
! == 5 单元残差应力节点力.
! == 6 单元分布载荷节点力.
! == 7 单元高斯点应力.
!........
USE CtrlData
USE ElmtData
USE MeshData
USE GlobData
USE FactorData
IMPLICIT DOUBLE PRECISION( a-h, o-z )
DIMENSION trans( 12, 12 ), fwork( 3 )
IF( iswth .EQ. 0 ) THEN
IF( nprnc .NE. 0 ) WRITE( 12, 2000 )
READ( 11, 2200 ) ( props( iprop, imats ), iprop = 2, 17 ), &
nfree, nstrp, ktype, &
( props( iprop, imats ), iprop = 21, &
20 + nstrp * 2 )
props( 18, imats ) = nfree
props( 19, imats ) = nstrp
props( 20, imats ) = ktype
IF( nprnc .NE. 0 ) &
WRITE( 12, 2400 ) ( props( iprop, imats ), iprop = 2, 17 ), &
( props( iprop, imats ), iprop =21, &
nstrp * 2 + 20 )
nstrh = MAX( nstrh, 12 )
RETURN
ELSE IF( iswth .EQ. 8 ) THEN
nnode = 2
ndofn = 6
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
!.......初时化数组
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
!.......检查单元的有效性
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 )
IF( fleng .LT. 1.0D-8 ) THEN
WRITE( 12, * ) ' 错误: 单元', ielem, '长度为零.'
RETURN
END IF
!.......提取单元的材料信息
elast = props( 2, imats ) ! 杨氏模量
shear = props( 3, imats ) ! 剪切模量
varea = props( 4, imats ) ! 截面面积
bendz = props( 5, imats ) ! XOZ平面弯曲惯性矩
bendy = props( 6, imats )
bendx = props( 7, imats )
facty = props( 8, imats )
factz = props( 9, imats )
qaxis = props( 10, imats )
qpxoy = props( 11, imats )
qpxoz = props( 12, imats )
denst = props( 13, imats )
kfree = props( 18, imats ) + 0.5D0
nstrp = props( 19, imats ) + 0.5D0
ktype = props( 20, imats ) + 0.5D0
CALL SetTheta( angle, imats )
CALL TransMatrix( trans, angle )
IF( iswth .EQ. 1 .OR. iswth .EQ. 5 .OR. &
iswth .EQ. 10 .OR. MOD( iswth, mswth ) .EQ. 7 ) THEN
!.........计算单元刚度
cnsty = 1.0D0 + facty
cnstz = 1.0D0 + factz
tempy = 4.0D0 + facty
tempz = 4.0D0 + factz
worky = 2.0D0 - facty
workz = 2.0D0 - factz
flen2 = fleng * fleng
flen3 = flen2 * fleng
stife( 1, 1 ) = elast * varea / fleng
stife( 4, 4 ) = shear * bendx / fleng
stife( 7, 7 ) = elast * varea / fleng
stife( 10, 10 ) = shear * bendx / fleng
stife( 10, 4 ) =-shear * bendx / fleng
stife( 7, 1 ) =-elast * varea / fleng
stife( 6, 2 ) = 6.0D0 * elast * bendz / flen2 / cnsty
stife( 12, 2 ) = 6.0D0 * elast * bendz / flen2 / cnsty
stife( 9, 5 ) = 6.0D0 * elast * bendy / flen2 / cnstz
stife( 11, 9 ) = 6.0D0 * elast * bendy / flen2 / cnstz
stife( 12, 8 ) =-6.0D0 * elast * bendz / flen2 / cnsty
stife( 5, 3 ) =-6.0D0 * elast * bendy / flen2 / cnstz
stife( 11, 3 ) =-6.0D0 * elast * bendy / flen2 / cnstz
stife( 8, 6 ) =-6.0D0 * elast * bendz / flen2 / cnsty
stife( 8, 8 ) = 12.D0 * elast * bendz / flen3 / cnsty
stife( 9, 9 ) = 12.D0 * elast * bendy / flen3 / cnstz
stife( 2, 2 ) = 12.D0 * elast * bendz / flen3 / cnsty
stife( 3, 3 ) = 12.D0 * elast * bendy / flen3 / cnstz
stife( 8, 2 ) =-12.D0 * elast * bendz / flen3 / cnsty
stife( 9, 3 ) =-12.D0 * elast * bendy / flen3 / cnstz
stife( 5, 5 ) = tempz * elast * bendy / fleng / cnstz
stife( 6, 6 ) = tempy * elast * bendz / fleng / cnsty
stife( 11, 11 ) = tempz * elast * bendy / fleng / cnstz
stife( 12, 12 ) = tempy * elast * bendz / fleng / cnsty
stife( 11, 5 ) = workz * elast * bendy / fleng / cnstz
stife( 12, 6 ) = worky * elast * bendz / fleng / cnsty
CALL ToGlobal( trans, iswth )
DO idofe = 1, 12
force( idofe ) = 0.0D0
DO jdofe = 1, 12
force( idofe ) = force( idofe ) + &
dispe( jdofe ) * stife( idofe, jdofe )
END DO
END DO
IF( MOD( iswth, mswth ) .EQ. 7 ) THEN
WRITE( 12, 2800 ) ielem
WRITE( 12, 3000 ) lnode(1), ( force(idofn), idofn = 1, 6 )
WRITE( 12, 3000 ) lnode(2), ( force(idofn), idofn = 7,12 )
IF( nstrp .GT. 0 .AND. nstrp .LE. 4 ) THEN
WRITE( 12, 3200 )
DO inode = 1, 2
locat = inode * 6 - 6
ipoin = lnode( inode )
DO istrp = 1, nstrp
yloca = props( 19 + istrp * 2, imats )
zloca = props( 20 + istrp * 2, imats )
strss = force( locat + 1 ) / varea + &
force( locat + 5 ) * zloca / bendy - &
force( locat + 6 ) * yloca / bendz
IF( inode .EQ. 1 ) strss =-strss
WRITE( 12, 3400 ) ipoin, yloca, zloca, strss
END DO
END DO
END IF
ELSE
DO idofe = 1, 12
IF( nincc .EQ. 2 .OR. iswth .EQ. 10 ) THEN
force( idofe ) =-force( idofe )
ELSE
force( idofe ) = 0.0D0
END IF
END DO
CALL LinkMatrix004
CALL LinkForce004
END IF
ELSE IF( iswth .EQ. 2 ) THEN
!.........质量矩阵
flen2 = fleng * fleng
weigh = fleng * varea * denst
stife( 1, 1 ) = weigh / 3.0D0
stife( 7, 7 ) = weigh / 3.0D0
stife( 2, 2 ) = weigh * 13.D0 / 35.D0 + &
weigh * 1.2D0 * bendz / varea / flen2
stife( 8, 8 ) = weigh * 13.D0 / 35.D0 + &
weigh * 1.2D0 * bendz / varea / flen2
stife( 3, 3 ) = weigh * 13.D0 / 35.D0 + &
weigh * 1.2D0 * bendy / varea / flen2
stife( 9, 9 ) = weigh * 13.D0 / 35.D0 + &
weigh * 1.2D0 * bendy / varea / flen2
stife( 4, 4 ) = weigh * bendx / 3.0D0 / varea
stife(10,10 ) = weigh * bendx / 3.0D0 / varea
stife( 5, 5 ) = weigh * flen2 / 105.0 + &
weigh * 2.0D0 * bendy / 15.D0 / varea
stife(11,11 ) = weigh * flen2 / 105.0 + &
weigh * 2.0D0 * bendy / 15.D0 / varea
stife( 6, 6 ) = weigh * flen2 / 105.0 + &
weigh * 2.0D0 * bendz / 15.D0 / varea
stife(12,12 ) = weigh * flen2 / 105.0 + &
weigh * 2.0D0 * bendz / 15.D0 / varea
stife( 6, 2 ) = weigh * 11.D0 * fleng / 210.0 + &
weigh * bendz / 10.D0 / varea / fleng
stife(12, 8 ) =-weigh * 11.D0 * fleng / 210.0 - &
weigh * bendz / 10.D0 / varea / fleng
stife( 5, 2 ) =-weigh * 11.D0 * fleng / 210.0 - &
weigh * bendy / 10.D0 / varea / fleng
stife(11, 8 ) = weigh * 11.D0 * fleng / 210.0 + &
weigh * bendy / 10.D0 / varea / fleng
stife( 8, 6 ) = weigh * 13.D0 * fleng / 420.0 - &
weigh * bendz / 10.D0 / varea / fleng
stife(12, 2 ) =-weigh * 13.D0 * fleng / 420.0 + &
weigh * bendz / 10.D0 / varea / fleng
stife( 9, 5 ) =-weigh * 13.D0 * fleng / 420.0 + &
weigh * bendy / 10.D0 / varea / fleng
stife(11, 3 ) = weigh * 13.D0 * fleng / 420.0 - &
weigh * bendy / 10.D0 / varea / fleng
stife( 7, 1 ) = weigh / 6.0D0
stife( 8, 2 ) = weigh * 9.0D0 / 70.0D0 - &
weigh * 1.2D0 * bendz / varea / flen2
stife( 9, 3 ) = weigh * 9.0D0 / 70.0D0 - &
weigh * 1.2D0 * bendy / varea / flen2
stife(10, 4 ) = weigh * bendx / 6.0D0 / varea
stife(11, 5 ) =-weigh * flen2 / 140.0 - &
weigh * bendy / 30.D0 / varea
stife(12, 6 ) =-weigh * flen2 / 140.0 - &
weigh * bendz / 30.D0 / varea
CALL ToGlobal( trans, iswth )
CALL LinkMatrix004
ELSE IF( iswth .EQ. 4 ) THEN
!.........几何矩阵
daxis = xleng * ( dispe( 7 ) - dispe( 1 ) ) / fleng + &
yleng * ( dispe( 8 ) - dispe( 2 ) ) / fleng + &
zleng * ( dispe( 9 ) - dispe( 3 ) ) / fleng
faxis = elast * varea * daxis / fleng / fleng
flen2 = fleng * fleng
const = bendy + bendz
stife( 2, 2 ) = faxis * 1.2D0
stife( 8, 8 ) = faxis * 1.2D0
stife( 3, 3 ) = faxis * 1.2D0
stife( 9, 9 ) = faxis * 1.2D0
stife( 8, 2 ) =-faxis * 1.2D0
stife( 9, 3 ) =-faxis * 1.2D0
stife( 4, 4 ) = faxis * const / varea
stife(10,10 ) = faxis * const / varea
stife(10, 4 ) =-faxis * const / varea
stife( 6, 2 ) = faxis * fleng / 10.D0
stife(12, 2 ) = faxis * fleng / 10.D0
stife( 8, 6 ) =-faxis * fleng / 10.D0
stife(12, 8 ) =-faxis * fleng / 10.D0
stife( 9, 5 ) =-faxis * fleng / 10.D0
stife(11, 9 ) =-faxis * fleng / 10.D0
stife(11, 5 ) =-faxis * flen2 / 30.D0
stife(12, 6 ) =-faxis * flen2 / 30.D0
stife( 5, 5 ) = faxis * 2.0D0 * flen2 / 15.D0
stife(11,11 ) = faxis * 2.0D0 * flen2 / 15.D0
stife( 6, 6 ) = faxis * 2.0D0 * flen2 / 15.D0
stife(12,12 ) = faxis * 2.0D0 * flen2 / 15.D0
CALL ToGlobal( trans, iswth )
CALL LinkMatrix004
ELSE IF( iswth .EQ. 6 ) THEN
!.........计算单元载荷
IF( ktype .NE. 0 ) THEN
fwork( 1 ) = qaxis
fwork( 2 ) = qpxoy
fwork( 3 ) = qpxoz
qaxis = trans( 1, 1 ) * fwork( 1 ) + &
trans( 1, 2 ) * fwork( 2 ) + &
trans( 1, 3 ) * fwork( 3 )
qpxoy = trans( 2, 1 ) * fwork( 1 ) + &
trans( 2, 2 ) * fwork( 2 ) + &
trans( 2, 3 ) * fwork( 3 )
qpxoz = trans( 3, 1 ) * fwork( 1 ) + &
trans( 3, 2 ) * fwork( 2 ) + &
trans( 3, 3 ) * fwork( 3 )
END IF
force( 1 ) = 0.5D0 * qaxis * fleng
force( 7 ) = 0.5D0 * qaxis * fleng
force( 2 ) = 0.5D0 * qpxoy * fleng
force( 3 ) = 0.5D0 * qpxoz * fleng
force( 8 ) = 0.5D0 * qpxoy * fleng
force( 9 ) = 0.5D0 * qpxoz * fleng
force( 5 ) = qpxoz * flen2 / 12.D0
force( 6 ) = qpxoy * flen2 / 12.D0
force( 11 ) =-qpxoz * flen2 / 12.D0
force( 12 ) =-qpxoy * flen2 / 12.D0
DO idofe = 1, 12
force( idofe ) = force( idofe ) - strsh( idofe, ielem )
END DO
CALL ToGlobal( trans, iswth )
CALL LinkForce004
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -