📄 elmt300.f90
字号:
props( 15, imats ) = nbody
props( 16, imats ) = imtyp
IF( imtyp .EQ. 1 ) THEN
IF( nincc .EQ. 1 .OR. nincc .EQ. 3 ) nerrc = 38556
IF( ngaus .LT. 1 ) ngaus = 27
IF( nerrc .NE. 0 ) RETURN
khstr = ngaus * 8
nincc = 2
IF( nhstr .LT. khstr ) nhstr = khstr
END IF
IF( nprnc .NE. 0 ) THEN
WRITE( 12, 5400 )
IF( imtyp .EQ. 0 ) WRITE( 12, 5600 )
IF( imtyp .EQ. 1 ) WRITE( 12, 5601 )
IF( imtyp .EQ. 0 ) THEN
IF( qmt33 .LT. 1.0D-20 ) THEN
WRITE( 12, 5800 ) qmt11, qmt22
ELSE
WRITE( 12, 6000 ) qmt11, qmt22, qmt33, qmt12, qmt23, &
qmt13, qmt44, qmt55, qmt66
END IF
ELSE IF( imtyp .EQ. 1 ) THEN
WRITE( 12, 6200 ) qmt11, qmt22, qmt33, qmt12, qmt23, &
qmt13
END IF
WRITE( 12, 6400 ) denst, xbody, ybody, zbody, nbody
END IF
END IF
2000 FORMAT( 12I5 )
2200 FORMAT( 2x, '分块数目: ', I5 )
2400 FORMAT( 2x, '几何特性: 几何线性' )
2401 FORMAT( 2x, '几何特性: 几何非线性' )
2600 FORMAT( 2x, '致命错误: 材料性质参数数至少为', I15 )
2800 FORMAT( 2x, '各块特性: 序号 几何 材料' )
3000 FORMAT(10x, I4, 4x, I4, 4x, I4 )
3200 FORMAT( 2x, '错误：输入了无效的几何类型序号', I5 )
3400 FORMAT( 2x, '错误：输入了无效的材料类型序号', I5 )
3600 FORMAT( 2x, '单元类型: 三维实体元(板壳元)分域几何' )
3800 FORMAT( 2x, '单元类型: 三维实体元' )
3801 FORMAT( 2x, '单元类型: 一般壳单元' )
3802 FORMAT( 2x, '单元类型: 薄板壳单元' )
3803 FORMAT( 2x, '单元类型: 三维梁单元' )
3804 FORMAT( 2x, '单元类型: 三维细长梁单元' )
3805 FORMAT( 2x, '单元类型: 三维杆件单元' )
4000 FORMAT( 2x, 'x-高斯点:', I15 / 2x, 'y-高斯点:', I15 / &
2x, 'z-高斯点:', I15 )
4200 FORMAT( 3F15.5 )
4400 FORMAT( 2X, '节点坐标:', I15, 3F10.5 )
4600 FORMAT( 12I5 )
4800 FORMAT( 2X, '分块定义:', 10X, 10I5 / 22X, 10I5 )
5000 FORMAT( 2x, '错误：输入了无效的内部节点序号', I5 )
5200 FORMAT( 4F15.5 / 4F15.5 / 4F15.5 / F15.5, 4I5 )
5400 FORMAT( 2x, '单元类型: 三维实体元(板壳元)分层材料' )
5600 FORMAT( 2x, '单元材料: 线弹性材料' )
5601 FORMAT( 2x, '单元材料: 蠕变材料' )
5800 FORMAT( 2x, '弹性模量:', E15.5/2X, '泊松比 :', F15.5 )
6000 FORMAT( 2X, '模量 D11:', E15.5/2X, '模量 D22:', E15.5/ &
2X, '模量 D33:', E15.5/2X, '模量 D12:', E15.5/ &
2X, '模量 D23:', E15.5/2X, '模量 D13:', E15.5/ &
2X, '模量 D44:', E15.5/2X, '模量 D55:', E15.5/ &
2X, '模量 D66:', E15.5 )
6200 FORMAT( 2X, '短期体积模量:', E15.5/2X, '长期体积模量:', E15.5/ &
2X, '体积粘性系数:', E15.5/2X, '短期剪切模量:', E15.5/ &
2X, '长期剪切模量:', E15.5/2X, '剪切粘性系数:', E15.5 )
6400 FORMAT( 2X, '质量密度:', E15.5/2X, 'x-体积力:', E15.5/ &
2X, 'y-体积力:', E15.5/2X, 'z-体积力:', E15.5/ &
2X, '动载类型:', I15 )
RETURN
END
SUBROUTINE MidPln300( coreb, corep, lnodb, ielem, nnodb )
USE CtrlData
USE ElmtData
IMPLICIT DOUBLE PRECISION( a-h, o-z )
DIMENSION xloct( 8 ), corep( 3, 8 ), coreb( 3, 20 )
DIMENSION yloct( 8 ), shape( 4, 21 ), lnodb( 20 )
DATA xloct / -1.0D0, 1.0D0, 1.0D0,-1.0D0, &
0.0D0, 1.0D0, 0.0D0,-1.0D0 /
DATA yloct / -1.0D0,-1.0D0, 1.0D0, 1.0D0, &
-1.0D0, 0.0D0, 1.0D0, 0.0D0 /
iflag = 1
DO inode = 1, 8
xloca = xloct( inode )
yloca = yloct( inode )
zloca = 0.0D0
CALL Shap3D( xloca, yloca, zloca, coreb, shape, xjaco, &
3, nnodb, lnodb, ielem, iflag, nerrc )
IF( nerrc .NE. 0 ) RETURN
xloca = 0.0D0
yloca = 0.0D0
zloca = 0.0D0
DO inodb = 1, nnodb
xloca = xloca + shape( 4, inodb ) * coreb( 1, inodb )
yloca = yloca + shape( 4, inodb ) * coreb( 2, inodb )
zloca = zloca + shape( 4, inodb ) * coreb( 3, inodb )
END DO
CALL Shap3D( xloca, yloca, zloca, coren, shape, xjaco, &
ndime, nnode, lnode, ielem, iflag, nerrc )
IF( nerrc .NE. 0 ) RETURN
DO idime = 1, ndime
corep( idime, inode ) = 0.0D0
DO jnode = 1, nnode
corep( idime, inode ) = corep( idime, inode ) + &
shape( 4, jnode ) * coren( idime, jnode )
END DO
END DO
END DO
RETURN
END
SUBROUTINE Bloc300( coreb, lnodb, igeob, nnodb )
USE CtrlData
USE MeshData
IMPLICIT DOUBLE PRECISION( a-h, o-z )
DIMENSION coreb( 3, 20 ), lnodb( 20 )
DO inodb = 1, 20
lnodb( inodb ) = 0
coreb( 1, inodb ) = 0.0D0
coreb( 2, inodb ) = 0.0D0
coreb( 3, inodb ) = 0.0D0
END DO
npntb = props( 3, igeob ) + 0.5D0
nnodb = props( 4, igeob ) + 0.5D0
DO inodb = 1, nnodb
iloca = 3 * npntb + inodb + 7
ipntb = props( iloca, igeob ) + 0.5D0
lnodb( inodb ) = ipntb
IF( ipntb .GT. 0 ) THEN
coreb( 1, inodb ) = props( ipntb * 3 + 5, igeob )
coreb( 2, inodb ) = props( ipntb * 3 + 6, igeob )
coreb( 3, inodb ) = props( ipntb * 3 + 7, igeob )
END IF
END DO
RETURN
END
SUBROUTINE Shape300( shape, coreb, lnodb, xloca, yloca, &
zloca, xjaco, nnodb, ielem, iflag )
USE CtrlData
USE ElmtData
IMPLICIT DOUBLE PRECISION( a-h, o-z )
DIMENSION shape( 4, 21 )
DIMENSION coreb( 3, 20 ), lnodb( 20 )
CALL Shap3D( xloca, yloca, zloca, coreb, shape, trans, &
ndime, nnodb, lnodb, ielem, iflag, nerrc )
IF( nerrc .GT. 0 ) RETURN
xloca = 0.0D0
yloca = 0.0D0
zloca = 0.0D0
DO inodb = 1, nnodb
xloca = xloca + shape( 4, inodb ) * coreb( 1, inodb )
yloca = yloca + shape( 4, inodb ) * coreb( 2, inodb )
zloca = zloca + shape( 4, inodb ) * coreb( 3, inodb )
END DO
CALL Shap3D( xloca, yloca, zloca, coren, shape, xjaco, &
ndime, nnode, lnode, ielem, iflag, nerrc )
xjaco = xjaco * trans
RETURN
END
SUBROUTINE BMat300( shape, bmatx, nlinc, iswth )
!........
! 模块功能
! 计算B矩阵
! nlinc 不等于零时考虑几何非线性
! iswth == 1 计算应变的B矩阵；
! == 2 计算切线刚度阵的B矩阵；
!........
USE ElmtData
IMPLICIT DOUBLE PRECISION ( a-h , o-z )
DIMENSION shape( 4, 21 ), bmatx( 6, 63 )
DIMENSION amatx( 6, 9 ), gmatx( 9, 63 ), agmat( 6, 63 )
DO idofe = 1, 63
DO istre = 1, 6
bmatx( istre, idofe ) = 0.0D0
END DO
END DO
DO inode = 1, nnode
locat = ( inode - 1 ) * 3
bmatx( 1, locat + 1 ) = shape( 1, inode )
bmatx( 2, locat + 2 ) = shape( 2, inode )
bmatx( 3, locat + 3 ) = shape( 3, inode )
bmatx( 4, locat + 2 ) = shape( 3, inode )
bmatx( 4, locat + 3 ) = shape( 2, inode )
bmatx( 5, locat + 1 ) = shape( 3, inode )
bmatx( 5, locat + 3 ) = shape( 1, inode )
bmatx( 6, locat + 1 ) = shape( 2, inode )
bmatx( 6, locat + 2 ) = shape( 1, inode )
END DO
IF( nlinc .NE. 0 ) THEN
const = 0.5D0 * iswth
CALL AMat300( amatx, shape )
CALL GMat300( gmatx, shape, nnode )
CALL AGmat300( amatx, agmat, gmatx, ndofe )
DO istre = 1, 6
DO idofe = 1, ndofe
bmatx( istre, idofe ) = bmatx( istre, idofe ) + &
agmat( istre, idofe ) * const
END DO
END DO
END IF
RETURN
END
SUBROUTINE DMat300( dmatx, imats, corep, xloca, yloca, &
ielem, igeob, imatb, iswth )
USE CtrlData
USE MeshData
IMPLICIT DOUBLE PRECISION ( a-h , o-z )
DIMENSION dmatx( 6, 6 ), corep( 3, 8 )
IF( iswth .EQ. 0 ) pcons = 0.0D0
IF( iswth .NE. 0 ) pcons = 1.0D3
DO istre = 1, 6
DO jstre = 1, 6
dmatx( istre, jstre ) = 0.0D0
END DO
END DO
qmt11 = props( 2, imatb )
qmt22 = props( 3, imatb )
qmt33 = props( 4, imatb )
qmt12 = props( 5, imatb )
qmt23 = props( 6, imatb )
qmt13 = props( 7, imatb )
qmt44 = props( 8, imatb )
qmt55 = props( 9, imatb )
qmt66 = props( 10, imatb )
imtyp = props( 16, imatb ) + 0.5D0
ietyp = props( 2, igeob ) + 0.5D0
IF( imtyp .EQ. 1 .OR. qmt33 .LT. 1.0D-10 ) THEN
IF( imtyp .EQ. 1 ) THEN
CALL Vsel300( imatb, elast, shear, poiss, vnorm )
ELSE IF( qmt33 .LT. 1.0D-20 ) THEN
elast = qmt11
poiss = qmt22
vnorm = elast
shear = elast / ( 1.0D0 + poiss ) / 2.0D0
END IF
IF( ietyp .EQ. 0 ) THEN
const = 2.0D0 * shear / ( 1.0D0 - 2.0D0 * poiss )
qmt11 = ( 1.0D0 - poiss ) * const
qmt22 = ( 1.0D0 - poiss ) * const
qmt33 = ( 1.0D0 - poiss ) * const
qmt12 = poiss * const
qmt13 = poiss * const
qmt23 = poiss * const
qmt44 = shear
qmt55 = shear
qmt66 = shear
ELSE IF( ietyp .EQ. 1 .OR. ietyp .EQ. 2 ) THEN
qmt11 = elast / ( 1.0D0 - poiss * poiss )
qmt22 = elast / ( 1.0D0 - poiss * poiss )
qmt33 = vnorm * pcons
qmt12 = elast * poiss / ( 1.0D0 - poiss * poiss )
qmt13 = 0.0D0
qmt23 = 0.0D0
qmt44 = shear
qmt55 = shear
qmt66 = shear
ELSE IF( ietyp .EQ. 3 .OR. ietyp .EQ. 4 ) THEN
qmt11 = elast
qmt22 = vnorm * pcons
qmt33 = vnorm * pcons
qmt12 = 0.0D0
qmt13 = 0.0D0
qmt23 = 0.0D0
qmt44 = shear
qmt55 = shear
qmt66 = shear
ELSE IF( ietyp .EQ. 5 ) THEN
qmt11 = elast
qmt22 = elast * pcons
qmt33 = elast * pcons
qmt12 = 0.0D0
qmt13 = 0.0D0
qmt23 = 0.0D0
qmt44 = 0.0D0
qmt55 = 0.0D0
qmt66 = 0.0D0
END IF
END IF
IF( ietyp .EQ. 2 ) THEN
qmt44 = qmt44 * pcons
qmt55 = qmt55 * pcons
ELSE IF( ietyp .EQ. 4 ) THEN
qmt55 = qmt55 * pcons
qmt66 = qmt66 * pcons
END IF
dmatx( 1, 1 ) = qmt11
dmatx( 2, 2 ) = qmt22
dmatx( 3, 3 ) = qmt33
dmatx( 1, 2 ) = qmt12
dmatx( 2, 1 ) = qmt12
dmatx( 2, 3 ) = qmt23
dmatx( 3, 2 ) = qmt23
dmatx( 1, 3 ) = qmt13
dmatx( 3, 1 ) = qmt13
dmatx( 4, 4 ) = qmt44
dmatx( 5, 5 ) = qmt55
dmatx( 6, 6 ) = qmt66
IF( ietyp .EQ. 0 ) RETURN
CALL Trans300( dmatx, corep, xloca, yloca, ielem )
RETURN
END
SUBROUTINE Trans300( dmatx, corep, xloca, yloca, ielem )
IMPLICIT DOUBLE PRECISION( a-h, o-z )
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -