📄 elmt400.f90
字号:
SUBROUTINE Elmt400( ielem, imats, iswth, iuses )
!.......
! 模块功能
! 4-8节点弹性地基单元(两参数Winkler地基模型).
!.......
USE CtrlData
USE MeshData
USE ElmtData
USE GlobData
USE FactorData
IMPLICIT DOUBLE PRECISION( a-h, o-z )
DIMENSION vnorm( 3 )
DIMENSION pygas( 6 ), wygas( 6 ), crden( 2, 8 )
DIMENSION pxgas( 6 ), wxgas( 6 ), shape( 3, 8 )
IF( iswth .EQ. 0 ) THEN
CALL RdMat400( imats )
RETURN
ELSE IF( iswth .EQ. 8 ) THEN
nnode = 8
ndofn = ndime
RETURN
ELSE IF( iswth .EQ. 9 ) THEN
knode = MIN0( nnode, 8 )
IF( iuses .EQ. 0 ) RETURN
WRITE( 15, 2000 ) 4, ielem, imats, knode, &
( lnode( inode ), inode = 1, knode )
RETURN
ELSE IF( iswth .EQ. 21 ) THEN
WRITE( 25 ) 0
RETURN
END IF
iflag = 2
mdime = 2
epara = props( 2, imats )
spara = props( 3, imats )
nxgas = props( 7, imats ) + 0.5D0
nygas = props( 8, imats ) + 0.5D0
vnorm( 1 ) = props( 4, imats )
vnorm( 2 ) = props( 5, imats )
vnorm( 3 ) = props( 6, imats )
CALL Trans400( crden, vnorm )
CALL Gauss( pxgas, wxgas, nxgas, nerrc )
CALL Gauss( pygas, wygas, nygas, nerrc )
IF( nerrc .GT. 0 ) RETURN
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
IF( iswth .EQ. 1 .OR. iswth .EQ. 5 .OR. iswth .EQ. 10 ) THEN
!.........计算单元刚度矩阵.
DO ixgas = 1, nxgas
DO iygas = 1, nygas
xloca = pxgas( ixgas )
yloca = pygas( iygas )
CALL Shap2D( shape, crden, xloca, yloca, lnode, xjaco, &
mdime, mdime, nnode, ielem, iflag, nerrc )
IF( nerrc .GT. 0 ) RETURN
dvolu = wxgas( ixgas ) * wygas( iygas ) * xjaco
DO inode = 1, nnode
DO jnode = 1, nnode
const = epara * shape(3, inode) * shape(3, jnode) + &
spara * shape(1, inode) * shape(1, jnode) + &
spara * shape(2, inode) * shape(2, jnode)
DO idofn = 1, ndofn
DO jdofn = 1, ndofn
iloca = ( inode - 1 ) * ndofn + idofn
jloca = ( jnode - 1 ) * ndofn + jdofn
stife( iloca, jloca ) = stife( iloca, jloca ) + &
vnorm( idofn ) * vnorm( jdofn ) * const * dvolu
END DO
END DO
END DO
END DO
END DO
END DO
!.........计算内力等效节点力
IF( nincc .EQ. 2 .OR. iswth .EQ. 10 ) THEN
DO idofe = 1, ndofe
DO jdofe = 1, ndofe
force( idofe ) = force( idofe ) - &
stife( idofe, jdofe ) * dispe( jdofe )
END DO
END DO
END IF
ELSE IF( iswth .EQ. 2 ) THEN
ELSE IF( iswth .EQ. 3 ) THEN
ELSE IF( iswth .EQ. 4 ) THEN
ELSE IF( iswth .EQ. 5 ) THEN
ELSE IF( iswth .EQ. 6 ) THEN
ELSE IF( MOD( iswth, mswth ) .EQ. 7 ) THEN
END IF
2000 FORMAT( 20I5 )
2200 FORMAT( /' 致命错误：单元', I5, '为不可压缩材料' / &
' 单元类型20-24只适用可压缩材料!' )
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 ' )
RETURN
END
SUBROUTINE RdMat400( imats )
USE CtrlData
USE MeshData
IMPLICIT DOUBLE PRECISION( a-h, o-z )
IF( nprop .LT. 8 ) THEN
WRITE( 12, 2000 ) 8
nerrc = 56450
RETURN
END IF
READ( 11, 2200 ) epara, spara, xnorm, ynorm, znorm, &
nxgas, nygas
vnorm = xnorm * xnorm + ynorm * ynorm + znorm * znorm
xnorm = xnorm / DSQRT( vnorm )
ynorm = ynorm / DSQRT( vnorm )
znorm = znorm / DSQRT( vnorm )
IF( nxgas .LE. 0 ) nxgas = 3
IF( nygas .LE. 0 ) nxgas = nxgas
props( 2, imats ) = epara
props( 3, imats ) = spara
props( 4, imats ) = xnorm
props( 5, imats ) = ynorm
props( 6, imats ) = znorm
props( 7, imats ) = nxgas
props( 8, imats ) = nygas
IF( nprnc .NE. 0 ) THEN
WRITE( 12, 2400 )
WRITE( 12, 2600 ) epara, spara, nxgas, nygas, &
xnorm, ynorm, znorm
END IF
2000 FORMAT( 2x, "致命错误：材料参数数目至少为", I5 )
2200 FORMAT( 4E15.5 / E15.5, 2I5 )
2400 FORMAT( 2x, "单元类型： 文克尔地基单元" )
2600 FORMAT( 2x, '基床系数:', E15.5/2x, '剪切模量:', E15.5 / &
2x, 'x-高斯点:', I15 /2x, 'y-高斯点:', I15 / &
2x, '法向矢量:', 5x, 3F10.5 )
RETURN
END
SUBROUTINE Trans400( crden, vnorm )
USE CtrlData
USE ElmtData
IMPLICIT DOUBLE PRECISION( a-h, o-z )
DIMENSION crden( 2, 8 ), vnorm( 3 )
DIMENSION fwork( 3, 8 ), vaxis( 3, 2 )
IF( ndime .EQ. 2 ) THEN
DO inode = 1, nnode
DO idime = 1, ndime
crden( idime, inode ) = coren( idime, inode )
END DO
END DO
ELSE
DO inode = 1, nnode
DO idime = 1, ndime
fwork( idime, inode ) = coren( idime, inode ) - &
coren( idime, 1 )
END DO
END DO
xaxis = fwork( 1, 2 ) * vnorm( 1 ) + &
fwork( 2, 2 ) * vnorm( 2 ) + &
fwork( 3, 2 ) * vnorm( 3 )
vaxis( 1, 1 ) = fwork( 1, 2 ) - xaxis * vnorm( 1 )
vaxis( 2, 1 ) = fwork( 2, 2 ) - xaxis * vnorm( 2 )
vaxis( 3, 1 ) = fwork( 3, 2 ) - xaxis * vnorm( 3 )
xaxis = DSQRT( vaxis( 1, 1 ) * vaxis( 1, 1 ) + &
vaxis( 2, 1 ) * vaxis( 2, 1 ) + &
vaxis( 3, 1 ) * vaxis( 3, 1 ) )
vaxis( 1, 1 ) = vaxis( 1, 1 ) / xaxis
vaxis( 2, 1 ) = vaxis( 2, 1 ) / xaxis
vaxis( 3, 1 ) = vaxis( 3, 1 ) / xaxis
vaxis( 1, 2 ) = vnorm( 2 ) * vaxis( 3, 1 ) - &
vnorm( 3 ) * vaxis( 2, 1 )
vaxis( 2, 2 ) = vnorm( 3 ) * vaxis( 1, 1 ) - &
vnorm( 1 ) * vaxis( 3, 1 )
vaxis( 3, 2 ) = vnorm( 1 ) * vaxis( 2, 1 ) - &
vnorm( 2 ) * vaxis( 1, 1 )
yaxis = DSQRT( vaxis( 1, 2 ) * vaxis( 1, 2 ) + &
vaxis( 2, 2 ) * vaxis( 2, 2 ) + &
vaxis( 3, 2 ) * vaxis( 3, 2 ) )
vaxis( 1, 2 ) = vaxis( 1, 2 ) / yaxis
vaxis( 2, 2 ) = vaxis( 2, 2 ) / yaxis
vaxis( 3, 2 ) = vaxis( 3, 2 ) / yaxis
DO inode = 1, nnode
crden( 1, inode ) = 0.0D0
crden( 2, inode ) = 0.0D0
DO idime = 1, ndime
crden( 1, inode ) = crden( 1, inode ) + &
vaxis( idime, 1 ) * coren( idime, inode )
crden( 2, inode ) = crden( 2, inode ) + &
vaxis( idime, 2 ) * coren( idime, inode )
END DO
END DO
END IF
RETURN
END
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -