📄 cacul_8.for
字号:
!*********************************************************************
!* 这是一个形成三维八节点单元刚度矩阵的程序。 *
!* *
!* ---文件----- *
!* 材料参数文件：mat.dat *
!* 实例的坐标文件：ELCODF.dat *
!* 刚度矩阵（下三角存储）：ESTIF.txt *
!* ---程序说明--- *
!* 主程序：CACUL_8_3D *
!* 弹性矩阵的子程序:MODPS(Y,P) *
!* 形函数及其导数值子程序：SFR2(R,S,T) *
!* 雅可比矩阵子程序：JACOB2(DJACB,ELCOD) *
!* 应变矩阵子程序：BMATPS *
!* *
!* 北京交大桥结所 *
!* 2003.9.1 *
!*********************************************************************
program CACUL_8_3D
IMPLICIT REAL*8(A-H,O-Z)
DIMENSION PROPS(2)
DIMENSION POSGP(4),WEIGP(4),ESTIF(300),ELCOD(3,8),DBMAT(6)
COMMON/BMDMP/BMATX(6,24),DMATX(6,6)
OPEN(1,FILE='mat.dat',status='old')
OPEN(2,FILE='ELCODF.dat',status='old')
OPEN(3,FILE='GAOSP.txt',status='replace')
OPEN(4,FILE='ESTIF.txt',status='replace')
! 读入材料参数
read(1,*) PROPS
read(2,*) ELCOD
! write (3,*) ELCOD(2,7)
! 以下进行的是高斯积分点位置赋值。
POSGP(1)=-0.861136311594053
POSGP(2)=-0.339981043584856
POSGP(3)=-POSGP(2)
POSGP(4)=-POSGP(1)
! 以下进行的是加权系数赋值。
WEIGP(1)=0.347854845137454
WEIGP(2)=0.652145154862546
WEIGP(3)=WEIGP(2)
WEIGP(4)=WEIGP(1)
! 以下进行的是应力(应变)分量个数赋值。
NSTRE=6
! 以下进行的是形成弹性矩阵[D]。
YOUNG=PROPS(1) !弹性模量
POISS=PROPS(2) !泊松比
CALL MODPS(YOUNG,POISS)
! 以下进行的是存放[K]的一维数组赋初值零。
KEVAB=0
DO 20 I=1,24
DO 20 J=1,I
KEVAB=KEVAB+1
ESTIF(KEVAB)=0.0
20 CONTINUE
! 以下进行的是计算积分点处形函数
! 及其导数之值以及对整体坐标的偏导数之值。
KGASP=0
DO 80 IGAUS=1,4 !高斯积分点数取4
DO 80 JGAUS=1,4
DO 80 KGAUS=1,4
KGASP=KGASP+1
EXISP=POSGP(IGAUS)
ETASP=POSGP(JGAUS)
EYBSP=POSGP(KGAUS)
WRITE(3,705) KGASP
705 FORMAT(6X,'高斯积分点编号为：',I5)
CALL SFR2(EXISP,ETASP,EYBSP)
CALL JACOB2(DJACB,ELCOD)
DVOLU=DJACB*WEIGP(IGAUS)*WEIGP(JGAUS)*WEIGP(KGAUS)
CALL BMATPS
! 以下进行的是弹性矩阵[D]赋初值零。
KEVAB=0
DO 70 IEVAB=1,24
DO 31 ISTRE=1,NSTRE
DBMAT(ISTRE)=0.0
DO 30 JSTRE=1,NSTRE
DBMAT(ISTRE)=DBMAT(ISTRE)+BMATX(JSTRE,IEVAB)*DMATX(JSTRE,ISTRE)
30 CONTINUE
31 CONTINUE
DO 60 JEVAB=1,IEVAB
KEVAB=KEVAB+1
BTDBM=0.0
DO 50 ISTRE=1,NSTRE
BTDBM=BTDBM+DBMAT(ISTRE)*BMATX(ISTRE,JEVAB)
50 CONTINUE
ESTIF(KEVAB)=ESTIF(KEVAB)+BTDBM*DVOLU
60 CONTINUE
70 CONTINUE
80 CONTINUE
! 以下为输出及其格式
write(4,*) "采用下三角存储，以下是每行对应元素的刚度"
write(4,*) "-----------------------------------------------"
write(4,*) "有问题请发邮件到bj-jd@sohu.com或电话：010-51681171"
DO 102,J=1,24
write(4,*) " "
write(4,5001) J
5001 FORMAT(8X,'---------------第',I2,'行-----------------')
write(4,*) " "
DO 103,I=1,J
WRITE (4,5002) I,ESTIF(I*(I-1)*0.5+J)
5002 FORMAT(14X,'第',I2,'列:',e15.3)
103 CONTINUE
102 CONTINUE
write(4,*) " ----------------文件结束!-------"
END program CACUL_8_3D
!**************************************
! 这是一个形成弹性矩阵的子程序。*
!**************************************
SUBROUTINE MODPS(Y,P)
IMPLICIT REAL*8(A-H,O-Z)
COMMON/BMDMP/BMATX(6,24),DMATX(6,6)
DO 100 I=1,6
DO 100 J=1,6
DMATX(I,J)=0.0
100 CONTINUE
CONST1=Y*(1.0-P)/((1.0+P)*(1.0-2.0*p))
CONST2=Y*P/((1.0+P)*(1.0-2.0*p))
CONST3=Y/(2.0*(1.0+P))
DMATX(1,1)=CONST1
DMATX(2,2)=CONST1
DMATX(3,3)=CONST1
DMATX(1,2)=CONST2
DMATX(2,1)=CONST2
DMATX(3,1)=CONST2
DMATX(1,3)=CONST2
DMATX(2,3)=CONST2
DMATX(3,2)=CONST2
DMATX(4,4)=CONST3
DMATX(5,5)=CONST3
DMATX(6,6)=CONST3
RETURN
END
!*******************************************************************
! 计算形函数当前积分点及形函数对局部坐标的导数值子程序。 *
!*******************************************************************
SUBROUTINE SFR2(R,S,T)
IMPLICIT REAL*8(A-H,O-Z)
COMMON/SD/SHAPEN(8),DERIV(3,8)
! 以下进行的是给形函数赋值。
SHAPEN(1)=(1-R)*(1-S)*(1-T)/8.0
SHAPEN(2)=(1+R)*(1-S)*(1-T)/8.0
SHAPEN(3)=(1+R)*(1+S)*(1-T)/8.0
SHAPEN(4)=(1-R)*(1+S)*(1-T)/8.0
SHAPEN(5)=(1-R)*(1-S)*(1+T)/8.0
SHAPEN(6)=(1+R)*(1-S)*(1+T)/8.0
SHAPEN(7)=(1+R)*(1+S)*(1+T)/8.0
SHAPEN(8)=(1-R)*(1+S)*(1+T)/8.0
! 以下进行的是给形函数对局部坐标的导数赋值。
DERIV(1,1)=(-1+S)*(1-T)/8.0
DERIV(1,2)=(1-S)*(1-T)/8.0
DERIV(1,3)=(1+S)*(1-T)/8.0
DERIV(1,4)=(-1-S)*(1-T)/8.0
DERIV(1,5)=(-1+S)*(1+T)/8.0
DERIV(1,6)=(1-S)*(1+T)/8.0
DERIV(1,7)=(1+S)*(1+T)/8.0
DERIV(1,8)=(-1-S)*(1+T)/8.0
DERIV(2,1)=(-1+R)*(1-T)/8.0
DERIV(2,2)=(-1-R)*(1-T)/8.0
DERIV(2,3)=(1+R)*(1-T)/8.0
DERIV(2,4)=(1-R)*(1-T)/8.0
DERIV(2,5)=(-1+R)*(1+T)/8.0
DERIV(2,6)=(-1-R)*(1+T)/8.0
DERIV(2,7)=(1+R)*(1+T)/8.0
DERIV(2,8)=(1-R)*(1+T)/8.0
DERIV(3,1)=(-1+R)*(1-S)/8.0
DERIV(3,2)=(-1-R)*(1-S)/8.0
DERIV(3,3)=(-1-R)*(1+S)/8.0
DERIV(3,4)=(-1+R)*(1+S)/8.0
DERIV(3,5)=(1-R)*(1-S)/8.0
DERIV(3,6)=(1+R)*(1-S)/8.0
DERIV(3,7)=(1+R)*(1+S)/8.0
DERIV(3,8)=(1-R)*(1+S)/8.0
RETURN
END
!*******************************************
! 这是一个形成雅可比矩阵的子程序。 *
!*******************************************
SUBROUTINE JACOB2(DJACB,ELCOD)
IMPLICIT REAL*8(A-H,O-Z)
DIMENSION XJACM(3,3),XJACI(3,3),ELCOD(3,8)
COMMON/SD/SHAPEN(8),DERIV(3,8)
COMMON/CAR/CARTD(3,8)
! 以下进行的是形成雅可比矩阵[J]的各元素。
DO 2000 I=1,3
DO 2000 J=1,3
XJACM(I,J)=0.0
DO 2000 K=1,8
XJACM(I,J)=XJACM(I,J)+DERIV(I,K)*ELCOD(J,K)
2000 CONTINUE
! 以下进行的是计算雅可比行列式┃Ｊ┃的值。
DJACB=XJACM(1,1)*XJACM(2,2)*XJACM(3,3)+XJACM(1,2)*XJACM(2,3)
1 *XJACM(3,1)+XJACM(1,3)*XJACM(3,2)*XJACM(2,1)-XJACM(3,1)
2 *XJACM(2,2)*XJACM(1,3)-XJACM(1,2)*XJACM(2,1)*XJACM(3,3)
3 -XJACM(1,1)*XJACM(3,2)*XJACM(2,3)
WRITE(3,6000) DJACB
6000 FORMAT(14X,'雅可比行列式┃Ｊ┃的值为：',F12.5)
IF(DJACB.LT.1.E-6)THEN
WRITE(3,6100)
6100 FORMAT('雅可比行列式的值小于或等于零！')
STOP '****** 程序运行被中断于子程序JACOB2 ******'
END IF
XJACI(1,1)=(XJACM(2,2)*XJACM(3,3)-XJACM(2,3)*XJACM(3,2))/DJACB
XJACI(1,2)=(XJACM(3,1)*XJACM(2,3)-XJACM(2,1)*XJACM(3,3))/DJACB
XJACI(1,3)=(XJACM(2,1)*XJACM(3,2)-XJACM(2,2)*XJACM(3,1))/DJACB
XJACI(2,1)=(XJACM(1,3)*XJACM(3,2)-XJACM(1,2)*XJACM(3,3))/DJACB
XJACI(2,2)=(XJACM(1,1)*XJACM(3,3)-XJACM(1,3)*XJACM(3,1))/DJACB
XJACI(2,3)=(XJACM(1,2)*XJACM(3,1)-XJACM(1,1)*XJACM(3,2))/DJACB
XJACI(3,1)=(XJACM(1,3)*XJACM(3,2)-XJACM(1,2)*XJACM(3,3))/DJACB
XJACI(3,2)=(XJACM(3,1)*XJACM(2,2)-XJACM(1,2)*XJACM(2,3))/DJACB
XJACI(3,3)=(XJACM(1,1)*XJACM(2,2)-XJACM(1,2)*XJACM(2,1))/DJACB
! 以下进行的是计算形函数对整体坐标的导数。
DO 300 I=1,3
DO 300 K=1,8
CARTD(I,K)=0.0
DO 300 J=1,3
CARTD(I,K)=CARTD(I,K)+XJACI(I,J)*DERIV(J,K)
300 CONTINUE
RETURN
END
!******************************************
! 这是一个形成应变矩阵的子程序。 *
!******************************************
SUBROUTINE BMATPS
IMPLICIT REAL*8(A-H,O-Z)
COMMON/CAR/CARTD(3,8)
COMMON/SD/SHAPEN(8),DERIV(3,8)
COMMON/BMDMP/BMATX(6,24),DMATX(6,6)
! 以下进行的是应变矩阵[B]赋初值零。
DO 108 I=1,24
DO 108 J=1,6
108 BMATX(J,I)=0.0
! 以下进行的是形成应变矩阵[B]。
KGASH=0
DO 200 I=1,8
MGASH=KGASH+1
NGASH=MGASH+1
KGASH=NGASH+1
BMATX(1,MGASH)=CARTD(1,I)
BMATX(2,NGASH)=CARTD(2,I)
BMATX(3,KGASH)=CARTD(3,I)
BMATX(4,MGASH)=CARTD(2,I)
BMATX(4,NGASH)=CARTD(1,I)
BMATX(5,NGASH)=CARTD(3,I)
BMATX(5,KGASH)=CARTD(2,I)
BMATX(6,MGASH)=CARTD(3,I)
BMATX(6,KGASH)=CARTD(1,I)
200 CONTINUE
RETURN
END⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -