?? 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 + -