一个计算八节点实体元刚度矩阵的源程序

ilxy

!*********************************************************************                                                               
!*           ---文件-----                                                                                                   *
!*        材料参数文件：mat.dat                                                                    *
!*        实例的坐标文件：ELCODF.dat                                                                   *
!*        刚度矩阵（下三角存储）：ESTIF.txt                                                    *
!*    ---程序说明---                                                                                                   *
!*             主程序：CACUL_8_3D                                                                                   *
!*              弹性矩阵的子程序:MODPS(Y,P)                                                                   *
!*                形函数及其导数值子程序：SFR2(R,S,T)                                                   *
!*                雅可比矩阵子程序：JACOB2(DJACB,ELCOD)                                   *
!*        应变矩阵子程序：BMATPS                                                                       *
!*                                                                                                                                !*ESTIF(300) 下三角单刚矩阵
!ELCOD(3,8) element coordination vector
!*********************************************************************
        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,*) "-----------------------------------------------"
        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

ilxy

此程序是一个大程序的一部分，经过验证为正确的，欢迎使用。参数的说明在程序中都已标出。

cstarfeng

太好了!请email 一份:

xchen1@hotmail.com

谢谢

[ 本帖最后由 cstarfeng 于 2007-6-21 08:03 编辑 ]

lja999

楼主 还有后面的程序吗？ 只介绍了 弹性矩阵 应变矩阵 可是 还是没有刚度矩阵部分的程序来着

liuichini

楼主 还有后面的程序吗？ 只介绍了 弹性矩阵 应变矩阵 可是 还是没有刚度矩阵部分的程序来着
lja999 发表于 2011-4-6 21:48

认真看程序！

googlchina

请问楼主是使用什么语言写的？

jackieyaolu

应该是fortran吧

麦勒小一

谢谢楼主啊