📄 elmt001.f90
字号:
SUBROUTINE Elmt001( ielem, imats, iswth, iuses )
!........
! 模块功能
! 杆单元, 同时包括一维, 二维, 三维杆单元.
!........
USE CtrlData
USE ElmtData
USE GlobData
USE MeshData
USE FactorData
IMPLICIT DOUBLE PRECISION( a-h, o-z )
DIMENSION direc( 3, 3 ), bmatx( 6 )
IF( iswth .EQ. 0 ) THEN
nprnc = ielem
IF( nstrh .LT. 1 ) nstrh = 1
IF( nprnc .NE. 0 ) WRITE( 12, 2000 )
READ( 11, 2200 ) elast, varea, denst, itype, nlinc, &
xbody, ybody, zbody, nbody
IF( nprnc .NE. 0 ) &
WRITE( 12, 2400 ) elast, varea, denst, itype, nlinc, &
xbody, ybody, zbody, nbody
props( 2, imats ) = elast
props( 3, imats ) = varea
props( 4, imats ) = denst
props( 5, imats ) = itype
props( 6, imats ) = nlinc
props( 7, imats ) = xbody
props( 8, imats ) = ybody
props( 9, imats ) = zbody
props( 10, imats ) = nbody
IF( nlinc .NE. 0 ) THEN
IF( nincc .EQ. 3 ) nerrc = 45864
nincc = 2
END IF
RETURN
ELSE IF( iswth .EQ. 8 ) THEN
nnode = 2
ndofn = ndime
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
force( idofe ) = 0.0D0
DO jdofe = 1, ndofe
stife( idofe, jdofe ) = 0.0D0
END DO
END DO
IF( iuses .EQ. 0 ) RETURN
!.......
elast = props( 2, imats )
varea = props( 3, imats )
denst = props( 4, imats )
itype = props( 5, imats ) + 0.5D0
nlinc = props( 6, imats ) + 0.5D0
if( itype .GT. 0 ) itype = 1
if( itype .LT. 0 ) itype =-1
CALL Axis001( direc, fleng )
CALL UpdateElast( elast, nlinc )
CALL Strs001( direc, fleng, elast, vstrs, nlinc )
IF( iswth .EQ. 1 .OR. iswth .EQ. 5 .OR. iswth .EQ. 10 ) THEN
const = elast * varea / fleng
CALL BMat001( bmatx, direc, nlinc, fleng )
DO idofe = 1, ndime * 2
DO jdofe = 1, ndime * 2
stife( idofe, jdofe ) = bmatx( idofe ) * &
const * bmatx( jdofe )
END DO
END DO
IF( nlinc .NE. 0 ) THEN
const = vstrs * varea / fleng
DO idime = 1, ndime
idofe = idime + ndime
stife( idime, idime ) = stife( idime, idime ) + const
stife( idime, idofe ) = stife( idime, idofe ) - const
stife( idofe, idime ) = stife( idofe, idime ) - const
stife( idofe, idofe ) = stife( idofe, idofe ) + const
END DO
END IF
CALL Trans001( direc, 1, 0 )
IF( nincc .EQ. 2 .OR. iswth .EQ. 10 ) THEN
const = varea * vstrs
DO idofe = 1, 2 * ndime
force( idofe ) = force( idofe ) - bmatx( idofe ) * const
END DO
CALL Trans001( direc, 0, 1 )
IF( DABS( strsh( 1, ielem ) ) .GT. 1.0D-8 ) THEN
vleng = 0.0D0
DO idime = 1, ndime
bmatx( idime ) = coren( idime, 2 ) - coren( idime, 1 ) + &
dispe( ndofn + idime ) - dispe( idime )
vleng = vleng + bmatx( idime ) * bmatx( idime )
END DO
vleng = DSQRT( vleng )
DO idime = 1, ndime
bmatx( idime ) = bmatx( idime ) / vleng
END DO
DO idime = 1, ndime
idofe = ndofn + idime
const = bmatx( idime ) * strsh( 1, ielem ) * varea
force( idime ) = force( idime ) + const
force( idofe ) = force( idofe ) - const
END DO
END IF
END IF
ELSE IF( iswth .EQ. 2 ) THEN
vmass = fleng * varea * denst / 2.0D0
DO idofe = 1, ndofe
stife( idofe, idofe ) = vmass
END DO
ELSE IF( iswth .EQ. 3 ) THEN
IF( itype .EQ. 0 ) RETURN
IF( itype .EQ.-1 ) iloca = 0
IF( itype .EQ. 1 ) iloca = ndime
const = DSQRT( elast * denst ) * varea
DO idime = 1, ndime
DO jdime = 1, ndime
idofe = iloca + idime
jdofe = iloca + jdime
stife( idofe, jdofe ) = stife( idofe, jdofe ) + const * &
direc( idime, 1 ) * direc( jdime, 1 )
END DO
END DO
ELSE IF( iswth .EQ. 4 ) THEN
faxis = vstrs * varea / fleng
DO idime = 1, ndime
DO jdime = 1, ndime
const = direc( idime, 1 ) * direc( jdime, 1 )
stife( idime, jdime ) =-const * faxis
END DO
stife( idime, idime ) = stife( idime, idime ) + faxis
END DO
DO idime = 1, ndime
DO jdime = 1, ndime
iloca = idime + ndime
jloca = jdime + ndime
stife( iloca, jdime ) =-stife( idime, jdime )
stife( idime, jloca ) =-stife( idime, jdime )
stife( iloca, jloca ) = stife( idime, jdime )
END DO
END DO
ELSE IF( iswth .EQ. 6 ) THEN
nbody = props( 10, imats ) + 0.5D0
CALL Factor( nbody )
IF( nerrc .GT. 0 ) RETURN
dvolu = varea * fleng * fctor / 2.0D0
idofe = 0
DO inode = 1, nnode
DO idime = 1, ndime
idofe = idofe + 1
force( idofe ) = dvolu * props( idime + 6, imats )
END DO
END DO
ELSE IF( MOD( iswth, mswth ) .EQ. 7 ) THEN
faxis = vstrs * varea
WRITE( 12, 2800 ) ielem, vstrs, faxis
END IF
2000 FORMAT( 2x, '单元类型: 杆单元' )
2200 FORMAT( 3F15.5, 2I5 / 3F15.5, I5 )
2400 FORMAT( 2x, '杨氏模量: ', E15.5/2x, '截面面积: ', E15.5 / &
2x, '质量密度: ', E15.5/2x, '边界约束: ', I15 / &
2x, '几何线性: ', I15 /2x, 'x-体积力: ', F15.5 / &
2x, 'y-体积力: ', F15.5/2x, 'z-体积力: ', F15.5 / &
2x, '动载类型: ', I15 )
2600 FORMAT( I3, 1X, I5, 1X, I4, 1X, I3, 8( 1X, I5 ) )
2800 FORMAT( 2x, '单元', I5, ' 应力', E15.6, ' 轴力', E15.6 )
RETURN
END
SUBROUTINE Axis001( direc, eleng )
!........
! 模块功能
! 计算杆单元的轴线方向和单元长度
!........
USE CtrlData
USE ElmtData
IMPLICIT DOUBLE PRECISION( a-h, o-z )
DIMENSION direc( 3, 3 )
DO idime = 1, 3
DO jdime = 1, 3
direc( idime, jdime ) = 0.0D0
END DO
END DO
eleng = 0.0D0
DO idime = 1, ndime
direc( idime, 1 ) = coren( idime, 2 ) - coren( idime, 1 )
eleng = eleng + direc( idime, 1 ) * direc( idime, 1 )
END DO
eleng = DSQRT( eleng )
DO idime = 1, ndime
direc( idime, 1 ) = direc( idime, 1 ) / eleng
END DO
fleng = DSQRT( direc( 1, 1 ) * direc( 1, 1 ) + &
direc( 2, 1 ) * direc( 2, 1 ) )
IF( fleng .LT. 1.0D-8 ) THEN
direc( 1, 2 ) = direc( 3, 1 )
direc( 2, 3 ) = 1.0D0
ELSE
direc( 1, 2 ) =-direc( 2, 1 ) / fleng
direc( 2, 2 ) = direc( 1, 1 ) / fleng
direc( 1, 3 ) =-direc( 2, 2 ) * direc( 3, 1 )
direc( 2, 3 ) = direc( 1, 2 ) * direc( 3, 1 )
direc( 3, 3 ) = direc( 1, 1 ) * direc( 2, 2 ) - &
direc( 1, 2 ) * direc( 2, 1 )
fleng = DSQRT( direc( 1, 3 ) * direc( 1, 3 ) + &
direc( 2, 3 ) * direc( 2, 3 ) + &
direc( 3, 3 ) * direc( 3, 3 ) )
direc( 1, 3 ) = direc( 1, 3 ) / fleng
direc( 2, 3 ) = direc( 2, 3 ) / fleng
direc( 3, 3 ) = direc( 3, 3 ) / fleng
END IF
RETURN
END
SUBROUTINE BMat001( bmatx, direc, nlinc, fleng )
USE CtrlData
USE ElmtData
IMPLICIT DOUBLE PRECISION( a-h, o-z )
DIMENSION bmatx( 6 ), direc( 3, 3 ), displ( 3, 2 )
DO idofe = 1, 6
bmatx( idofe ) = 0.0D0
END DO
bmatx( 1 ) = -1.0D0
bmatx( ndime + 1 ) = 1.0D0
IF( nlinc .EQ. 0 ) RETURN
DO inode = 1, 2
DO idime = 1, ndime
displ( idime, inode ) = 0.0D0
DO jdime = 1, ndime
idofe = ( inode - 1 ) * ndime + jdime
displ( idime, inode ) = displ( idime, inode ) + &
direc( jdime, idime ) * dispe( idofe )
END DO
END DO
END DO
DO idime = 1, ndime
idofe = idime + ndime
const = displ( idime, 2 ) - displ( idime, 1 )
bmatx( idime ) = bmatx( idime ) - const / fleng
bmatx( idofe ) = bmatx( idofe ) + const / fleng
END DO
RETURN
END
SUBROUTINE Strs001( direc, fleng, elast, vstrs, nlinc )
USE CtrlData
USE ElmtData
IMPLICIT DOUBLE PRECISION( a-h, o-z )
DIMENSION direc( 3, 3 ), displ( 3, 2 )
DO inode = 1, 2
DO idime = 1, ndime
displ( idime, inode ) = 0.0D0
DO jdime = 1, ndime
idofe = ( inode - 1 ) * ndime + jdime
displ( idime, inode ) = displ( idime, inode ) + &
direc( jdime, idime ) * dispe( idofe )
END DO
END DO
END DO
DO idime = 1, ndime
displ( idime, 1 ) = displ( idime, 2 ) - displ( idime, 1 )
displ( idime, 1 ) = displ( idime, 1 ) / fleng
END DO
vstrn = displ( 1, 1 )
IF( nlinc .NE. 0 ) THEN
!.........如果考虑几何非线性，先计算非线性应变
DO idime = 1, ndime
vstrn = vstrn + displ( idime, 1 ) * &
0.5D0 * displ( idime, 1 )
END DO
END IF
vstrs = vstrn * elast
RETURN
END
SUBROUTINE Trans001( direc, nstif, nforc )
USE CtrlData
USE ElmtData
IMPLICIT DOUBLE PRECISION( a-h, o-z )
DIMENSION direc( 3, 3 ), vwork( 3, 3 ), vtemp( 3 )
IF( nstif .NE. 0 ) THEN
DO inode = 1, 2
DO jnode = 1, 2
iloca = ( inode - 1 ) * ndime
jloca = ( jnode - 1 ) * ndime
DO idime = 1, ndime
DO jdime = 1, ndime
vwork( idime, jdime ) = 0.0D0
DO kdime = 1, ndime
idofe = iloca + idime
jdofe = jloca + kdime
vwork( idime, jdime ) = vwork( idime, jdime ) + &
stife( idofe, jdofe ) * direc( jdime, kdime )
END DO
END DO
END DO
DO idime = 1, ndime
DO jdime = 1, ndime
idofe = iloca + idime
jdofe = jloca + jdime
stife( idofe, jdofe ) = 0.0D0
DO kdime = 1, ndime
stife( idofe, jdofe ) = stife( idofe, jdofe ) + &
direc( idime, kdime ) * vwork( kdime, jdime )
END DO
END DO
END DO
END DO
END DO
END IF
IF( nforc .NE. 0 ) THEN
DO inode = 1, 2
DO idime = 1, ndime
vtemp( idime ) = 0.0D0
DO jdime = 1, ndime
idofe = ( inode - 1 ) * ndime + jdime
vtemp( idime ) = vtemp( idime ) + &
force( idofe ) * direc( idime, jdime )
END DO
END DO
DO idime = 1, ndime
idofe = ( inode - 1 ) * ndime + idime
force( idofe ) = vtemp( idime )
END DO
END DO
END IF
RETURN
END
SUBROUTINE UpdateElast( elast, nlinc )
USE CtrlData
USE ElmtData
IMPLICIT DOUBLE PRECISION( a-h, o-z )
IF( nlinc .EQ. 0 ) RETURN
dprev = 0.0D0
dcurr = 0.0D0
DO idime = 1, ndime
dtemp = coren( idime, 2 ) - coren( idime, 1 )
dwork = dispe( ndofn + idime ) - dispe( idime ) + dtemp
dprev = dprev + dtemp * dtemp
dcurr = dcurr + dwork * dwork
END DO
epson = ( dcurr - dprev ) / dprev
!.........如果给出的是现实构型中模量或大变形小应变
elast = elast / ( 1.0D0 + epson )
RETURN
END
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -