📄 elmt200.f90
字号:
SUBROUTINE Elmt200( ielem, imats, iswth, iuses )
!.......
! 模块功能
! 4-8节点平面单元的元素矩阵.
! ietyp == 0 平面应力梁单元;
! ietyp == 1 平面应变梁单元;
! ietyp == 2 平面应力单元;
! ietyp == 3 平面应变单元;
! ietyp == 4 平面轴对称单元.
! ietyp == 5 轴对称旋转壳单元.
!.......
USE CtrlData
USE MeshData
USE ElmtData
USE GlobData
USE FactorData
IMPLICIT DOUBLE PRECISION( a-h, o-z )
DIMENSION strgp( 6 ), direc( 3 )
DIMENSION wygas( 6 ), snode( 6, 8 )
DIMENSION pygas( 6 ), dmatx( 4, 4 )
DIMENSION pxgas( 6 ), shape( 3, 8 )
DIMENSION wxgas( 6 ), bmatx( 4, 16 )
DIMENSION fbody( 2 ), smatx( 4, 16 )
IF( iswth .EQ. 0 ) THEN
CALL RdMat200( imats )
RETURN
ELSE IF( iswth .EQ. 8 ) THEN
nnode = 8
ndofn = 2
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
nstre = 3
elast = props( 2, imats )
poiss = props( 3, imats )
thick = props( 4, imats )
denst = props( 5, imats )
nxgas = props( 8, imats )
nygas = props( 9, imats )
nbody = props( 10, imats )
pdamp = props( 11, imats )
sdamp = props( 12, imats )
ietyp = props( 13, imats ) + 0.5D0
IF( ietyp .GE. 4 ) nstre = 4
IF( DABS( poiss - 0.5D0 ) .LT. 1.0D-5 ) THEN
WRITE( 12, 2200 ) ielem
nerrc = 28432
RETURN
END IF
IF( nxgas .LE. 0 ) nxgas = 3
IF( nygas .LE. 0 ) nygas = nxgas
IF( thick .LE. 1.0D-8 ) thick = 1.0D0
CALL Gauss( pxgas, wxgas, nxgas, nerrc )
IF( nerrc .GT. 0 ) RETURN
CALL Gauss( pygas, wygas, nygas, nerrc )
IF( nerrc .GT. 0 ) RETURN
CALL Factor( nbody )
IF( nerrc .GT. 0 ) RETURN
fbody( 1 ) = props( 6, imats ) * fctor
fbody( 2 ) = props( 7, imats ) * fctor
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
!.........计算单元刚度矩阵.
kgaus = 0
DO ixgas = 1, nxgas
!...........计算高斯点切线方向
xloca = pxgas( ixgas )
IF( ietyp .LE. 1 .OR. ietyp .EQ. 5 ) &
CALL Tan200( direc, xloca )
DO iygas = 1, nygas
kgaus = kgaus + 1
yloca = pygas( iygas )
CALL Shap2D( shape, coren, xloca, yloca, lnode, xjaco, &
ndime, mdime, nnode, ielem, iflag, nerrc )
IF( nerrc .GT. 0 ) RETURN
!.............计算高斯点坐标
rcord = 0.0D0
zcord = 0.0D0
DO inode = 1, nnode
rcord = rcord + shape( 3, inode ) * coren( 1, inode )
zcord = zcord + shape( 3, inode ) * coren( 2, inode )
END DO
IF( ietyp .GE. 4 ) thick = 8.0D0 * DATAN( 1.0D0 ) * rcord
CALL Bmat200( shape, bmatx, rcord, ietyp )
CALL Dmat200( direc, dmatx, elast, poiss, 1, ietyp )
dvolu = wxgas( ixgas ) * wygas( iygas ) * xjaco * thick
DO istre = 1, nstre
DO idofe = 1, ndofe
smatx( istre, idofe ) = 0.0d0
DO jstre = 1, nstre
smatx( istre, idofe ) = smatx( istre, idofe ) + &
dmatx( istre, jstre ) * bmatx( jstre, idofe )
END DO
END DO
END DO
DO idofe = 1, ndofe
DO jdofe = 1, ndofe
DO istre = 1, nstre
stife( idofe, jdofe ) = stife( idofe, jdofe ) + &
bmatx( istre, idofe ) * smatx( istre, jdofe ) * &
dvolu
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
!.........计算单元一致质量矩阵.
DO ixgas = 1, nxgas
DO iygas = 1, nygas
xloca = pxgas( ixgas )
yloca = pygas( iygas )
CALL shap2d( shape, coren, xloca, yloca, lnode, xjaco, &
ndime, mdime, nnode, ielem, iflag, nerrc )
IF( nerrc .GT. 0 ) RETURN
IF( ietyp .GE. 4 ) THEN
thick = 0.0D0
DO inode = 1, nnode
thick = thick + shape( 3, inode ) * coren( 1, inode )
END DO
thick = thick * 8.0D0 * DATAN( 1.0D0 )
END IF
const = xjaco * denst * thick
dvolu = wxgas( ixgas ) * wygas( iygas ) * const
DO inode = 1, nnode
iloca = ( inode - 1 ) * 2
DO jnode = 1, nnode
jloca = ( jnode - 1 ) * 2
const = shape( 3, inode ) * shape( 3, jnode ) * dvolu
stife( iloca + 1, jloca + 1 ) = &
stife( iloca + 1, jloca + 1 ) + const
stife( iloca + 2, jloca + 2 ) = &
stife( iloca + 2, jloca + 2 ) + const
END DO
END DO
END DO
END DO
ELSE IF( iswth .EQ. 3 ) THEN
!.........计算单元阻尼矩阵.
kgaus = 0
IF( DABS( pdamp ) .LT. 1.0D-8 .AND. &
DABS( sdamp ) .LT. 1.0D-8 ) RETURN
elast = elast * 2.0D0 * sdamp / ( 1.0D0 + poiss ) * &
( ( 3.0D0 * pdamp - 2.0D0 * sdamp ) + &
( 4.0D0 * sdamp - 3.0D0 * pdamp ) * poiss ) / &
( ( 2.0D0 * pdamp - sdamp ) + 2.0D0 * &
( sdamp - pdamp ) * poiss )
poiss = ( ( pdamp - sdamp ) + &
( 2.0D0 * sdamp - pdamp ) * poiss ) / &
( ( 2.0D0 * pdamp - sdamp ) + 2.0D0 * &
( sdamp - pdamp ) * poiss )
DO ixgas = 1, nxgas
!...........计算高斯点切线方向
xloca = pxgas( ixgas )
IF( ietyp .LE. 1 .OR. ietyp .EQ. 5 ) &
CALL Tan200( direc, xloca )
DO iygas = 1, nygas
kgaus = kgaus + 1
yloca = pygas( iygas )
CALL Shap2D( shape, coren, xloca, yloca, lnode, xjaco, &
ndime, mdime, nnode, ielem, iflag, nerrc )
IF( nerrc .GT. 0 ) RETURN
!.............计算高斯点坐标
rcord = 0.0D0
zcord = 0.0D0
DO inode = 1, nnode
rcord = rcord + shape( 3, inode ) * coren( 1, inode )
zcord = zcord + shape( 3, inode ) * coren( 2, inode )
END DO
IF( ietyp .GE. 4 ) thick = 8.0D0 * DATAN( 1.0D0 ) * rcord
CALL Bmat200( shape, bmatx, rcord, ietyp )
CALL Dmat200( direc, dmatx, elast, poiss, 0, ietyp )
dvolu = wxgas( ixgas ) * wygas( iygas ) * xjaco * thick
DO istre = 1, nstre
DO idofe = 1, ndofe
smatx( istre, idofe ) = 0.0D0
DO jstre = 1, nstre
smatx( istre, idofe ) = smatx( istre, idofe ) + &
dmatx( istre, jstre ) * bmatx( jstre, idofe )
END DO
END DO
END DO
DO idofe = 1, ndofe
DO jdofe = 1, ndofe
DO istre = 1, nstre
stife( idofe, jdofe ) = stife( idofe, jdofe ) + &
bmatx( istre, idofe ) * smatx( istre, jdofe ) * &
dvolu
END DO
END DO
END DO
END DO
END DO
ELSE IF( iswth .EQ. 4 ) THEN
ELSE IF( iswth .EQ. 5 ) THEN
ELSE IF( iswth .EQ. 6 ) THEN
!.........set nodal force from body force.
DO ixgas = 1, nxgas
DO iygas = 1, nygas
xloca = pxgas( ixgas )
yloca = pygas( iygas )
CALL Shap2d( shape, coren, xloca, yloca, lnode, xjaco, &
ndime, mdime, nnode, ielem, iflag, nerrc )
IF( nerrc .GT. 0 ) RETURN
rcord = 0.0D0
zcord = 0.0D0
DO inode = 1, nnode
rcord = rcord + shape( 3, inode ) * coren( 1, inode )
zcord = zcord + shape( 3, inode ) * coren( 2, inode )
END DO
IF( ietyp .GE. 4 ) &
thick = 8.0D0 * DATAN( 1.0D0 ) * rcord
dvolu = wxgas( ixgas ) * wygas( iygas ) * xjaco * thick
DO inode = 1, nnode
DO idofn = 1, ndofn
idofe = ( inode - 1 ) * ndofn + idofn
force( idofe ) = force( idofe ) + dvolu * &
fbody( idofn ) * shape( 3, inode )
END DO
END DO
END DO
END DO
ELSE IF( MOD( iswth, mswth ) .EQ. 7 ) THEN
!.........计算高斯点应力.
IF( iswth .EQ. 7 ) THEN
WRITE( 12, 3800 ) ielem
ELSE
darea = 0.0D0
DO inode = 1, nnode
DO jnode = 1, nnode
stife( inode, jnode ) = 0.0D0
END DO
DO istre = 1, 6
snode( istre, inode ) = 0.0D0
END DO
END DO
END IF
DO ixgas = 1, nxgas
xloca = pxgas( ixgas )
IF( ietyp .LE. 1 .OR. ietyp .EQ. 5 ) &
CALL Tan200( direc, xloca )
DO iygas = 1, nygas
kgaus = kgaus + 1
yloca = pygas( iygas )
CALL Shap2d( shape, coren, xloca, yloca, lnode, xjaco, &
ndime, mdime, nnode, ielem, iflag, nerrc )
IF( nerrc .GT. 0 ) RETURN
rcord = 0.0D0
zcord = 0.0D0
DO inode = 1, nnode
rcord = rcord + shape( 3, inode ) * coren( 1, inode )
zcord = zcord + shape( 3, inode ) * coren( 2, inode )
END DO
IF( ietyp .GE. 4 ) &
thick = 8.0D0 * DATAN( 1.0D0 ) * rcord
CALL Bmat200( shape, bmatx, rcord, ietyp )
CALL Dmat200( direc, dmatx, elast, poiss, 0, ietyp )
DO istre = 1, nstre
smatx( istre, 1 ) = 0.0D0
DO idofe = 1, ndofe
smatx( istre, 1 ) = smatx( istre, 1 ) + &
bmatx( istre, idofe ) * dispe( idofe )
END DO
END DO
DO istre = 1, nstre
strgp( istre ) = 0.0D0
DO jstre = 1, nstre
strgp( istre ) = strgp( istre ) + &
dmatx( istre, jstre ) * smatx( jstre, 1 )
END DO
END DO
IF( ietyp .LE. 3 ) THEN
strgp( 6 ) = strgp( 3 )
strgp( 3 ) = 0.0D0
strgp( 4 ) = 0.0D0
strgp( 5 ) = 0.0D0
ELSE
strgp( 5 ) = strgp( 4 )
strgp( 4 ) = 0.0D0
strgp( 6 ) = 0.0D0
END IF
IF( iswth .EQ. 7 ) THEN
WRITE( 12, 4000 ) rcord, zcord , strgp(1), strgp(2), &
strgp(3), 0.0D0, strgp(4), strgp(5), strgp(6)
ELSE
iloca = 0
DO inode = 1, nnode
IF( lnode( inode ) .NE. 0 ) THEN
iloca = iloca + 1
jloca = 0
DO jnode = 1, nnode
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -