📄 elmt200.f90
字号:
IF( lnode( jnode ) .NE. 0 ) THEN
jloca = jloca + 1
stife( iloca, jloca ) = stife( iloca, jloca ) + &
shape( 3, inode ) * shape( 3, jnode )
END IF
END DO
DO istre = 1, 6
snode( istre, iloca ) = snode( istre, iloca ) + &
shape( 3, inode ) * strgp( istre )
END DO
END IF
END DO
darea = darea + wxgas( ixgas ) * &
xjaco * wygas( iygas )
END IF
END DO
END DO
IF( iswth .EQ. mswth + 7 ) THEN
knode = 0
DO inode = 1, nnode
IF( lnode( inode ) .NE. 0 ) knode = knode + 1
END DO
CALL Invers( stife, mdofe, knode, nerrc )
IF( nerrc .NE. 0 ) RETURN
DO istre = 1, 6
iloca = 0
DO inode = 1, nnode
ipoin = lnode( inode )
IF( ipoin .NE. 0 ) THEN
force( istre ) = 0.0D0
iloca = iloca + 1
jloca = 0
DO jnode = 1, nnode
IF( lnode( jnode ) .NE. 0 ) THEN
jloca = jloca + 1
force( istre ) = force( istre ) + &
stife( iloca, jloca ) * snode( istre, jloca )
END IF
END DO
strss( istre, ipoin ) = strss( istre, ipoin ) + &
force( istre ) * darea
IF( istre .EQ. 1 ) &
strss( 7, ipoin ) = strss( 7, ipoin ) + darea
END IF
END DO
END DO
END IF
END IF
2000 FORMAT( 20I5 )
2200 FORMAT( /' 致命错误：单元', I5, '为不可压缩材料' / &
' 单元类型200只适用可压缩材料!' )
3800 FORMAT( /20x, ' 单 元 ', i5 &
/ 7x, 'x', 10x, 'y', 10x, &
'S-xx', 10x, 'S-yy', 10x, 'S-zz' &
/ 5x, 'z(h)', 20x, 'S-yz', 10x, 'S-zx', 10x, 'S-xy ' )
4000 FORMAT( 2x, 2f10.3, 3e14.6 / 2x, f10.3, 10x, 3e14.6 )
RETURN
END
SUBROUTINE Bmat200( shape, bmatx, rcord, ietyp )
!.......
! 模块功能
! 计算平面单元的B矩阵.
!.......
IMPLICIT DOUBLE PRECISION( a-h, o-z )
DIMENSION shape( 3, 8 ), bmatx( 4, 16 )
IF( ietyp .LE. 3 ) THEN
!.........计算平面应力(变)单元和平面梁单元的B矩阵.
DO inode = 1, 8
locat = ( inode - 1 ) * 2
bmatx( 1, locat + 1 ) = shape( 1, inode )
bmatx( 2, locat + 1 ) = 0.0D0
bmatx( 3, locat + 1 ) = shape( 2, inode )
bmatx( 1, locat + 2 ) = 0.0D0
bmatx( 2, locat + 2 ) = shape( 2, inode )
bmatx( 3, locat + 2 ) = shape( 1, inode )
END DO
ELSE
!.........计算平面轴对称单元和轴对称壳单元的B矩阵.
DO inode = 1, 8
locat = ( inode - 1 ) * 2
bmatx( 1, locat + 1 ) = shape( 1, inode )
bmatx( 2, locat + 1 ) = shape( 3, inode ) / rcord
bmatx( 3, locat + 1 ) = 0.0D0
bmatx( 4, locat + 1 ) = shape( 2, inode )
bmatx( 1, locat + 2 ) = 0.0D0
bmatx( 2, locat + 2 ) = 0.0D0
bmatx( 3, locat + 2 ) = shape( 2, inode )
bmatx( 4, locat + 2 ) = shape( 1, inode )
END DO
END IF
RETURN
END
SUBROUTINE Dmat200( direc, dmatx, elast, poiss, iswth, ietyp )
!.......
! 模块功能
! 计算平面单元的D矩阵.
!.......
IMPLICIT DOUBLE PRECISION( a-h, o-z )
DIMENSION dmatx( 4, 4 ), trans( 4, 4 )
DIMENSION dmatw( 4, 4 ), direc( 3 )
!.......计算剪切模量, 初时化D矩阵.
shear = 0.5D0 * elast / ( 1.0D0 + poiss )
DO istre = 1, 4
DO jstre = 1, 4
trans( istre, jstre ) = 0.0D0
dmatw( istre, jstre ) = 0.0D0
dmatx( istre, jstre ) = 0.0D0
END DO
END DO
IF( ietyp .LE. 1 ) THEN
IF( ietyp .EQ. 0 ) THEN
!...........计算平面应力梁单元的D矩阵.
dmatx( 1, 1 ) = elast
dmatx( 2, 2 ) = elast * iswth
dmatx( 3, 3 ) = shear
ELSE
!...........计算平面应变梁单元的D矩阵.
dmatx( 1, 1 ) = elast / ( 1.0D0 - poiss * poiss )
dmatx( 2, 2 ) = iswth * dmatx( 1, 1 )
dmatx( 3, 3 ) = shear
END IF
trans( 1, 1 ) = direc( 1 ) * direc( 1 )
trans( 1, 2 ) = direc( 2 ) * direc( 2 )
trans( 1, 3 ) = direc( 1 ) * direc( 2 )
trans( 2, 1 ) = direc( 2 ) * direc( 2 )
trans( 2, 2 ) = direc( 1 ) * direc( 1 )
trans( 2, 3 ) =-direc( 1 ) * direc( 2 )
trans( 3, 1 ) =-trans( 1, 3 ) * 2.0D0
trans( 3, 2 ) = trans( 1, 3 ) * 2.0D0
trans( 3, 3 ) = trans( 1, 1 ) - trans( 1, 2 )
DO istre = 1, 3
DO jstre = 1, 3
dmatw( istre, jstre ) = 0.0D0
DO kstre = 1, 3
dmatw( istre, jstre ) = dmatw( istre, jstre ) + &
dmatx( istre, kstre ) * trans( kstre, jstre )
END DO
END DO
END DO
DO istre = 1, 3
DO jstre = 1, 3
dmatx( istre, jstre ) = 0.0D0
DO kstre = 1, 3
dmatx( istre, jstre ) = dmatx( istre, jstre ) + &
trans( kstre, istre ) * dmatw( kstre, jstre )
END DO
END DO
END DO
ELSE IF( ietyp .EQ. 2 ) THEN
!.........计算平面应力单元的D矩阵.
cnst1 = elast / ( 1.0D0 - poiss * poiss )
dmatx( 1, 1 ) = cnst1
dmatx( 2, 2 ) = cnst1
dmatx( 1, 2 ) = cnst1 * poiss
dmatx( 2, 1 ) = cnst1 * poiss
dmatx( 3, 3 ) = shear
ELSE IF( ietyp .EQ. 3 ) THEN
!.........计算平面应变单元的D矩阵.
cnst1 = elast * ( 1.0D0 - poiss ) / ( ( 1.0D0 + poiss ) * &
( 1.0D0 - 2.0D0 * poiss ) )
dmatx( 1, 1 ) = cnst1
dmatx( 2, 2 ) = cnst1
dmatx( 1, 2 ) = cnst1 * poiss / ( 1.0 - poiss )
dmatx( 2, 1 ) = cnst1 * poiss / ( 1.0 - poiss )
dmatx( 3, 3 ) = shear
ELSE IF( ietyp .EQ. 4 ) THEN
!.........计算平面轴对称单元的D矩阵.
cnst1 = elast * ( 1.0D0 - poiss ) / ( ( 1.0D0 + poiss ) * ( &
1.0D0 - 2.0D0 * poiss ) )
cnst2 = poiss / ( 1.0D0 - poiss )
cnst3 = 0.5D0 * ( 1.0D0 - 2.0D0 * poiss ) / ( 1.0D0 - poiss )
dmatx( 1, 1 ) = cnst1
dmatx( 2, 2 ) = cnst1
dmatx( 3, 3 ) = cnst1
dmatx( 4, 4 ) = cnst1 * cnst3
dmatx( 1, 2 ) = cnst1 * cnst2
dmatx( 2, 1 ) = dmatx( 1, 2 )
dmatx( 1, 3 ) = dmatx( 1, 2 )
dmatx( 3, 1 ) = dmatx( 1, 2 )
dmatx( 3, 2 ) = dmatx( 1, 2 )
dmatx( 2, 3 ) = dmatx( 1, 2 )
ELSE IF( ietyp .EQ. 5 ) THEN
!.........计算轴对称旋转壳单元的D矩阵.
!.........局部坐标下的材料系数
dmatx( 1, 1 ) = elast / ( 1.0D0 - poiss * poiss )
dmatx( 1, 2 ) = dmatx( 1, 1 ) * poiss
dmatx( 2, 1 ) = dmatx( 1, 2 )
dmatx( 2, 2 ) = dmatx( 1, 1 )
dmatx( 3, 3 ) = elast * iswth
dmatx( 4, 4 ) = elast * 0.5D0 / ( 1.0D0 + poiss )
!.........把D矩阵由局部坐标转换到整体坐标的变换矩阵
trans( 1, 1 ) = direc( 1 ) * direc( 1 )
trans( 1, 3 ) = direc( 2 ) * direc( 2 )
trans( 1, 4 ) = direc( 1 ) * direc( 2 )
trans( 2, 2 ) = 1.0D0
trans( 3, 1 ) = direc( 2 ) * direc( 2 )
trans( 3, 3 ) = direc( 1 ) * direc( 1 )
trans( 3, 4 ) =-direc( 1 ) * direc( 2 )
trans( 4, 1 ) =-trans( 1, 4 ) * 2.0D0
trans( 4, 3 ) = trans( 1, 4 ) * 2.0D0
trans( 4, 4 ) = trans( 1, 1 ) - trans( 1, 3 )
!.........D矩阵坐标变换
DO istre = 1, 4
DO jstre = 1, 4
dmatw( istre, jstre ) = 0.0D0
DO kstre = 1, 4
dmatw( istre, jstre ) = dmatw( istre, jstre ) + &
dmatx( istre, kstre ) * trans( kstre, jstre )
END DO
END DO
END DO
DO istre = 1, 4
DO jstre = 1, 4
dmatx( istre, jstre ) = 0.0D0
DO kstre = 1, 4
dmatx( istre, jstre ) = dmatx( istre, jstre ) + &
trans( kstre, istre ) * dmatw( kstre, jstre )
END DO
END DO
END DO
END IF
RETURN
END
SUBROUTINE Tan200( direc, xloca )
!........
! 模块功能
! 计算单元高斯点的切线方向
!........
USE CtrlData
USE ElmtData
IMPLICIT DOUBLE PRECISION( a-h, o-z )
DIMENSION direc( 3 ), lnodp( 3 )
DIMENSION corep( 9 ), shape( 2, 3 )
!.......计算中心线坐标
nnodp = 2
lnodp( 1 ) = 1
lnodp( 2 ) = 2
DO idime = 1, ndime
jdime = idime + ndime
corep( idime ) = ( coren( idime, 1 ) + &
coren( idime, 4 ) ) / 2.0D0
corep( jdime ) = ( coren( idime, 2 ) + &
coren( idime, 3 ) ) / 2.0D0
IF( nnode .GE. 7 .AND. lnode( 7 ) .NE. 0 ) THEN
nnodp = 3
lnodp( 3 ) = 3
jdime = idime + ndime + ndime
corep( jdime ) = ( coren( idime, 7 ) + &
coren( idime, 5 ) ) / 2.0D0
END IF
END DO
!.......计算中心线切线方向
iflag = 1
const = 0.0D0
CALL Shap1D( shape, corep, xloca, ndime, nnodp, lnodp, &
xjaco, iflag )
DO idime = 1, ndime
direc( idime ) = 0.0D0
DO inodp = 1, nnodp
locat = ( inodp - 1 ) * ndime + idime
direc( idime ) = direc( idime ) + &
corep( locat ) * shape( 1, inodp )
END DO
const = const + direc( idime ) * direc( idime )
END DO
const = DSQRT( const )
DO idime = 1, ndime
direc( idime ) = direc( idime ) / const
END DO
RETURN
END
SUBROUTINE RdMat200( imats )
USE CtrlData
USE MeshData
IMPLICIT DOUBLE PRECISION( a-h, o-z )
IF( nprop .LT. 13 ) nerrc = 10434
IF( nerrc .GT. 0 ) RETURN
READ( 11, 2200 ) elast, poiss, thick, denst, pdamp, &
sdamp, xbody, ybody, kbody, nxgas, &
nygas, ietyp
IF( nprnc .NE. 0 ) THEN
IF( ietyp .EQ. 0 ) WRITE( 12, 2400 )
IF( ietyp .EQ. 1 ) WRITE( 12, 2600 )
IF( ietyp .EQ. 2 ) WRITE( 12, 2800 )
IF( ietyp .EQ. 3 ) WRITE( 12, 3000 )
IF( ietyp .EQ. 4 ) WRITE( 12, 3200 )
IF( ietyp .EQ. 5 ) WRITE( 12, 3400 )
WRITE( 12, 3600 ) elast, poiss, thick, denst, pdamp, &
sdamp, xbody, ybody, kbody, nxgas, nygas
END IF
props( 2, imats ) = elast
props( 3, imats ) = poiss
props( 4, imats ) = thick
props( 5, imats ) = denst
props( 6, imats ) = xbody
props( 7, imats ) = ybody
props( 8, imats ) = nxgas
props( 9, imats ) = nygas
props( 10, imats ) = kbody
props( 11, imats ) = pdamp
props( 12, imats ) = sdamp
props( 13, imats ) = ietyp
2000 FORMAT( 2x, '致命错误：200号单元至少需要13个材料信息单元！' )
2200 FORMAT( 4F15.5 / 4F15.5 / 4I5 )
2400 FORMAT( 2x, '单元类型: 平面应力梁单元' )
2600 FORMAT( 2x, '单元类型: 平面应变梁单元' )
2800 FORMAT( 2x, '单元类型: 平面应力元' )
3000 FORMAT( 2x, '单元类型: 平面应变元' )
3200 FORMAT( 2x, '单元类型: 轴对称单元' )
3400 FORMAT( 2x, '单元类型: 轴对称旋转壳单元' )
3600 FORMAT( 2x, '杨氏模量:', E15.5/2x, '泊 松 比:', E15.5 / &
2x, '厚 度:', E15.5/2x, '质量密度:', E15.5 / &
2x, 'p-damp :', E15.5/2x, 's-damp :', E15.5 / &
2x, 'x-体积力:', E15.5/2x, 'y-体积力:', E15.5 / &
2x, '动载类型:', I15 /2x, 'x-高斯点:', I15 / &
2x, 'y-高斯点:', I15 )
RETURN
END
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -