论坛 › 基础力学理论 › 平面弹性力学有限元源程序（FORTRAN）

hao1982 发表于 2006-10-9 12:48

平面弹性力学有限元源程序（FORTRAN）

$DEBUG
                PROGRAM PLANE
                IMPLICIT REAL*8(A-H,O-Z),INTEGER(I-N)
                ALLOCATABLE::IJK(:,:),XY(:,:),BCA(:,:),SK(:,:),STR(:,:),MB(:,:),ZB(:),B(:)
                ALLOCATABLE::DELD(:,:,:),TOD(:,:),DELST(:,:,:),TOST(:,:),DELSUP(:,:),TOTSUP(:)
                DIMENSION EK(6,6)
      CHARACTER PN*40,FN*12

      WRITE(*,'(A)') ' 本程序为计算平面问题的有限元程序'
      WRITE(*,'(A)') ' 特点：(1)采用三结点三角形单元；'
      WRITE(*,'(A)') '       (2)采用等带宽存贮技术；'
      WRITE(*,'(A)') '       (3)采用高斯消元法解线性方程组。'
      WRITE(*,'(/A)') ' 输入计算问题名(PN):'
      READ(*,'(A)') PN
      CALL FNAME(PN,'.DAT',FN)
      WRITE(*,'(2A)') '输入数据文件名为：',FN
      OPEN(5,FILE=FN,STATUS='OLD')
      CALL FNAME(PN,'.OUT',FN)
      WRITE(*,'(/2A)') ' 结果输出数据文件名为： ',FN
      OPEN(6,FILE=FN,STATUS='UNKNOWN')
      CALL FNAME(PN,'.OU1',FN)
      WRITE(*,'(/2A)') ' 参数输出数据文件名为： ',FN
      OPEN(7,FILE=FN,STATUS='UNKNOWN')

                READ(5,*) NG,NE,MC,NX,NB,EO,VO,DENSITY,T
                WRITE(6,120) NG,NE,MC,NX,NB
                WRITE(6,130) EO,VO,DENSITY,T
                READ(5,*) NWA,NWE,NWK,NWP
                NT=2*NG
                ALLOCATE (IJK(3,NE),XY(2,NG),BCA(7,NE),STR(3,NE),MB(2,NB),ZB(NB),B(NT))
                ALLOCATE (DELD(2,NG,NX),TOD(2,NG),DELST(3,NE,NX),TOST(3,NE),DELSUP(NB,NX),TOTSUP(NB))
                CALL CLEAR(2,NG,TOD)
                CALL CLEAR(3,NE,TOST)
                CALL CLEAR1(NB,TOTSUP)
                IF (NG.EQ.0) THEN
                  STOP 000
                ENDIF
                CALL INPUT(NG,NE,NB,IJK,XY,MB,ZB)
                CALL CALND(NG,NE,IJK,ND)
                ALLOCATE (SK(NT,ND))
                IF(MC.EQ.0) THEN
                E=EO
                V=VO
                ELSE
                E=EO/(1-VO*VO)
                V=VO/(1-VO)
                ENDIF
                IP=0
                NX1=NX
                A1=E/(1-V*V)/4.0
                A2=0.5*(1-V)
                CALL ABC(NG,NE,IJK,XY,BCA,NWA)
                CALL CLEAR(NT,ND,SK)
                DO 100 K=1,NE
                CALL KE(K,NE,T,V,A1,A2,IJK,BCA,EK,NWE)
                CALL SUMK(K,NE,NT,ND,IJK,EK,SK)
100                CONTINUE
                CALL CHECK(NT,ND,SK)
110                CONTINUE
                IP=IP+1
                WRITE(6,'(///1X,A,I4)') "荷载工况=",IP
                READ(5,*) NF,NP,NM
                DO I=1,NT
                B(I)=0.0
                ENDDO
                IF(NF.GT.0) THEN
                CALL PF(NF,NT,B)
                ENDIF
                IF(NP.GT.0) THEN
                CALL PP(NP,NG,NT,T,XY,B)
                ENDIF
                IF(NM.GT.0) THEN
                CALL PG(NM,NE,NT,DENSITY,T,IJK,BCA,B)
                ENDIF
                DO I=1,NB
                I1=2*(MB(1,I)-1)+2-MB(2,I)
                DELSUP(I,IP)=-B(I1)
                ENDDO
                IF(NF.NE.0.OR.NP.NE.0.OR.NM.NE.0) THEN
                CALL DBC(NB,ND,NT,NX,NX1,MB,ZB,SK,B)
                CALL GAUSS(NT,ND,NX1,NX,SK,B,NWK,NWP)
                CALL STRESS(NE,NT,V,A1,A2,IJK,BCA,B,STR)
                CALL TRES(IP,NG,NE,NX,NT,B,STR,DELD,TOD,DELST,TOST)
                CALL SUPT(IP,NG,NE,NB,NX,T,V,A1,A2,IJK,MB,BCA,DELD,DELSUP,TOTSUP)
                CALL OUTPUT(IP,NG,NE,NB,NX,MB,DELD,TOD,DELST,TOST,DELSUP,TOTSUP)
                ELSE
                WRITE(*,'(2x,A,I4,A)') "荷载工况=",IP,"没有荷载!"
                WRITE(6,'(2x,A,I4,A)') "荷载工况=",IP,"没有荷载!"
                ENDIF
                NX1=NX1-1
                IF (NX1.GT.0) GOTO 110
120                FORMAT(1X, '结点总数=', I3, 1X, '单元总数=', I3, 1X, '问题类型=', &
      & I1, 1X, '荷载工况数=', I2, 1X, '支承位移数=', I2)
130                FORMAT(/1X, '弹性模量=', E10.3, 1X, '泊松比=', F5.3, 1X, '密度=', E10.3, &
      & 1X, '厚度=', F6.3)
                END

                               
                SUBROUTINE GAUSS(NT,ND,NX1,NX,SK,B,NWK,NWP)
                IMPLICIT REAL*8(A-H,O-Z),INTEGER(I-N)
                DIMENSION SK(NT,ND),A(NT,NT),B(NT)
                CALL CLEAR(NT,NT,A)
                DO I=1,NT
               DO J=1,ND
                  IF((I+J-1).LE.NT) THEN
                  A(I,I+J-1)=SK(I,J)
                  ENDIF
               ENDDO
                ENDDO
                IF(NWK.EQ.1.AND.NX1.EQ.NX) THEN
            WRITE(7,'(A)') "结构总刚(列出右上三角元素)"
                DO I=1,NT
                WRITE(7,100) (A(I,J),J=1,NT)
                ENDDO
                ENDIF
                IF (NWP.EQ.1) THEN
                WRITE(7,'(1X,A,I3,A)') "第",NX-NX1+1,"荷载工况的荷载列阵"
                DO I=1,NT
                WRITE(7,'(1x,E11.4)') B(I)
                ENDDO
                ENDIF
                DO 50 K=1,NT-1
               DO 60 I=K+1,NT
                  C1=A(K,I)/A(K,K)
                  DO 70 J=I,NT
                  A(I,J)=A(I,J)-C1*A(K,J)
70                  CONTINUE
                  B(I)=B(I)-C1*B(K)
60               CONTINUE
50                CONTINUE
                B(NT)=B(NT)/A(NT,NT)
                DO 80 I=NT-1,1,-1
               DO 90 J=I+1,NT
               B(I)=B(I)-A(I,J)*B(J)
90               CONTINUE
               B(I)=B(I)/A(I,I)
80                CONTINUE
100                FORMAT(1x, 4000E10.3)
                END


                SUBROUTINE CHECK(NT,ND,SK)
                IMPLICIT REAL*8(A-H,O-Z),INTEGER(I-N)
                DIMENSION SK(NT,ND)
                M=0
                DO I=1,NT
               IF(SK(I,1).LE.0.0) THEN
               M=M+1
               ENDIF
                ENDDO
                IF (M.GT.0) THEN
                WRITE(*,*) "总刚矩阵的对角元素为负值!!"
                STOP 2222
                ENDIF
                END

                SUBROUTINE STRESS(NE,NT,V,A1,A2,IJK,BCA,B,STR)
                IMPLICIT REAL*8(A-H,O-Z),INTEGER(I-N)
                DIMENSION D1(6),D2(3),B(NT),S(3,6),IJK(3,NE),BCA(7,NE),STR(3,NE)
                CALL CLEAR(3,NE,STR)
                DO 10 K=1,NE
               DO I=1,3
               II=IJK(I,K)
               D1(2*(I-1)+1)=B(2*(II-1)+1)
               D1(2*(I-1)+2)=B(2*(II-1)+2)
               ENDDO
               CALL CLEAR(3,6,S)
               C1=2*A1/BCA(7,K)
               DO 20 I=1,3
                  S(1,2*(I-1)+1)=C1*BCA(I,K)
                  S(1,2*(I-1)+2)=C1*V*BCA(I+3,K)
                  S(2,2*(I-1)+1)=C1*V*BCA(I,K)
                  S(2,2*(I-1)+2)=C1*BCA(I+3,K)
                  S(3,2*(I-1)+1)=C1*A2*BCA(I+3,K)
                  S(3,2*(I-1)+2)=C1*A2*BCA(I,K)
20               CONTINUE
               CALL MATMUL(3,6,1,S,D1,D2)
               DO I=1,3
               STR(I,K)=D2(I)
               ENDDO
10                CONTINUE
                END

                SUBROUTINE OUTPUT(IP,NG,NE,NB,NX,MB,DELD,TOD,DELST,TOST,DELSUP,TOTSUP)
                IMPLICIT REAL*8(A-H,O-Z),INTEGER(I-N)
                DIMENSION MB(2,NB),DELD(2,NG,NX),TOD(2,NG)
                DIMENSION DELST(3,NE,NX),TOST(3,NE),DELSUP(NB,NX),TOTSUP(NB)
                CHARACTER VE*6
                WRITE(6,20) IP
                WRITE(6,30)
                DO I=1,NG
                WRITE(6,40) I,DELD(1,I,IP),TOD(1,I),DELD(2,I,IP),TOD(2,I)
                ENDDO
                WRITE(6,50)
                DO J=1,NE
                WRITE(6,60) J,(DELST(L,J,IP),TOST(L,J),L=1,3)
                ENDDO
                WRITE(6,70)
                DO J=1,NB
                IF (MB(2,J).EQ.1) THEN
                VE='x方向'
                ELSE
                VE='Y方向'
                ENDIF
                WRITE(6,80) MB(1,J),VE,DELSUP(J,IP),TOTSUP(J)
                ENDDO
20                FORMAT(/1X, '荷载工况=', I3, '的计算结果')
30                FORMAT(/, 1X, '结点号', 5X, 'X位移增量', 5X, 'X累计位移', 5x, &
      & 'Y位移增量', 5X, 'Y累计位移')
40                FORMAT(1X, I3, 6X, F10.4, 4x, F10.4, 4X, F10.4, 4x, F10.4)
50                FORMAT(/, 1X, '单元号', 6x, 'X应力增量', 6x, 'X累计应力', 6x, &
      & 'Y应力增量', 6x, 'Y累计应力', 6x, '剪应力增量', 6x, &
      & '累计剪应力')
60                FORMAT(2X, I3, 6X, E10.3, 5X, E10.3, 5X, E10.3, 5X, E10.3, 6X, E10.3, &
      & 6X, E10.3)
70                FORMAT(/, 1X, '支座结点号', 6x, '支反力方向', 6x, '支反力增量', &
      & 6x, '支反力累计量')
80                FORMAT(5X, I3, 12X, A, 6x, E10.3, 8X, E10.3)
          END                                                         


                SUBROUTINE MATMUL(M,N,L,A,B,C)
                IMPLICIT REAL*8(A-H,O-Z),INTEGER(I-N)
      DIMENSION A(M,N),B(N,L),C(M,L)
      DO 100 I=1,M
         DO 100 J=1,L
          C(I,J)=0.0
          DO 100 K=1,N
100   C(I,J)=C(I,J)+A(I,K)*B(K,J)
      END

                SUBROUTINE CLEAR(M,N,A)
                IMPLICIT REAL*8(A-H,O-Z),INTEGER(I-N)
                DIMENSION A(M,N)
                DO I=1,M
                DO J=1,N
                A(I,J)=0.0
                ENDDO
                ENDDO
                END

                SUBROUTINE ABC(NG,NE,IJK,XY,BCA,NWA)
                IMPLICIT REAL*8(A-H,O-Z),INTEGER(I-N)
                DIMENSION X(2,5),XY(2,NG),IJK(3,NE),B(7),BCA(7,NE)
                DO 10 I=1,NE
               DO 20 K=1,3
               DO 20 J=1,2
               X(J,K)=XY(J,IJK(K,I))
20               CONTINUE
               DO 30 J=1,2
               X(J,4)=X(J,1)
               X(J,5)=X(J,2)
30               CONTINUE
               DO 40 K=1,3
               B(K)=X(2,K+1)-X(2,K+2)
               B(K+3)=X(1,K+2)-X(1,K+1)
40               CONTINUE
               B(7)=(B(1)*B(5)-B(4)*B(2))*0.5
               IF (B(7).LE.0.0) THEN
               WRITE(*,*) "单元号=",I
               WRITE(*,*) "单元的结点I,J AND K编号：",(IJK(K,I),K=1,3)
               STOP 1111
               ELSE
               DO J=1,7
               BCA(J,I)=B(J)
               ENDDO
               ENDIF
10                CONTINUE
                IF(NWA.EQ.1) THEN
                WRITE(7,100)
                DO I=1,NE
                WRITE(7,110) I,(BCA(J,I),J=1,7)
                ENDDO
                ENDIF
100               FORMAT(1X, '单元号', 3x, 'bi', 7x, 'bj', 7x, 'bk', 7x, 'ci', 7x, &
      & 'cj', 7x, 'ck', 7x, '面积A')
110               FORMAT(2X, I3, 3X, F6.3, 3X, F6.3, 3X, F6.3, 3X, F6.3, 3X, F6.3, &
      & 3X, F6.3, 3X, F6.3)
             END

               SUBROUTINE SUMK(K,NE,NT,ND,IJK,EK,SK)
             IMPLICIT REAL*8(A-H,O-Z),INTEGER(I-N)
               DIMENSION        IJ(3),EK(6,6),IJK(3,NE),SK(NT,ND)
               DO I=1,3
               IJ(I)=IJK(I,K)
               ENDDO
               DO 10 I=1,3
                  DO 20 J=1,3
                  IF (IJ(I).GT.IJ(J)) GOTO 20
                  M=2*(IJ(I)-1)+1
                  N=2*(IJ(J)-1)+1-(2*(IJ(I)-1)+1)+1
                  MO=2*(I-1)+1
                  NO=2*(J-1)+1
                  SK(M,N)=SK(M,N)+EK(MO,NO)
                  SK(M,N+1)=SK(M,N+1)+EK(MO,NO+1)
                  SK(M+1,N)=SK(M+1,N)+EK(MO+1,NO+1)
                  IF (IJ(I).EQ.IJ(J)) GOTO 20
                  SK(M+1,N-1)=SK(M+1,N-1)+EK(MO+1,NO)
20                  CONTINUE
10             CONTINUE
               END

                SUBROUTINE KE(K,NE,T,V,A1,A2,IJK,BCA,EK,NWE)
                IMPLICIT REAL*8(A-H,O-Z),INTEGER(I-N)
                DIMENSIONIJ(3),IJK(3,NE),BCA(7,NE),EK(6,6)
                CALL CLEAR(6,6,EK)
                DO I=1,3
                IJ(I)=IJK(I,K)
                ENDDO
                DO 10 I=1,3
               DO 20 J=I,3
               IO=2*(I-1)
               JO=2*(J-1)
               EK(IO+1,JO+1)=BCA(I,K)*BCA(J,K)+A2*BCA(I+3,K)*BCA(J+3,K)
               EK(IO+1,JO+2)=V*BCA(I,K)*BCA(J+3,K)+A2*BCA(I+3,K)*BCA(J,K)
               EK(IO+2,JO+2)=BCA(I+3,K)*BCA(J+3,K)+A2*BCA(I,K)*BCA(J,K)
               IF (I.EQ.J) GOTO 20
               EK(IO+2,JO+1)=V*BCA(I+3,K)*BCA(J,K)+A2*BCA(I,K)*BCA(J+3,K)
20               CONTINUE
10                CONTINUE
                DO I=2,6
               DO J=1,I-1
               EK(I,J)=EK(J,I)
               ENDDO
                ENDDO
                DO I=1,6
               DO J=1,6
               EK(I,J)=EK(I,J)*A1*T/BCA(7,K)
               ENDDO
                ENDDO
                IF(NWE.EQ.1) THEN
           WRITE(7,40)        K, "号单元的单刚"
                DO I=1,6
               WRITE(7,50) (EK(I,J),J=1,6)
                ENDDO
                ENDIF
40                FORMAT(1X, I3, A)
50                FORMAT(1X, 6(E10.4, 4X))
                END

                SUBROUTINE PF(NF,NT,B)
                IMPLICIT REAL*8(A-H,O-Z),INTEGER(I-N)
                DIMENSION MF(2,NF),ZF(NF),B(NT)
                DO I=1,NF
                READ(5,*) K,MF(1,I),MF(2,I),ZF(I)
                ENDDO
                WRITE(6,'(/A,I3)') "集中荷载个数NF=",NF
                WRITE(6,20)
                DO I=1,NF
                WRITE(6,30) MF(1,I),MF(2,I),ZF(I)
                ENDDO
                DO 10 I=1,NF
               II=2*(MF(1,I)-1)-MF(2,I)+2
               B(II)=B(II)+ZF(I)
10                CONTINUE
20                FORMAT(1X, '结点号', 4x, '方向', 4x, '荷载值')
30                FORMAT(2X, I3, 5X, I3, 4X, E10.4)
                END               

      SUBROUTINE FNAME(PN,FN2,FN)
      CHARACTER PN*40,FN2*4,FN*12
!       去掉PN中前面的空格
      DO 10 I=1,40
         IF(PN(I:I).EQ.' ') GOTO 10
         IP=I
         GOTO 20
10      CONTINUE
20      CONTINUE
      FN(1:8)=PN(IP:IP+7)
!       去掉FN中后面的空格
      DO 30 I=8,1,-1
         IF(FN(I:I).EQ.' ') GOTO 30
         IL=I
         GOTO 40
30      CONTINUE
40      CONTINUE
!       生成文件名FN=PN+FN2
      FN(IL+1:IL+4)=FN2(1:4)
      END

                SUBROUTINE PP(NP,NG,NT,T,XY,B)
                IMPLICIT REAL*8(A-H,O-Z),INTEGER(I-N)
                DIMENSION MP(2,NP),ZP(NP),XY(2,NG),B(NT)
                DO I=1,NP
                READ(5,*) K,MP(1,I),MP(2,I),ZP(I)
                ENDDO
                WRITE(6,'(/A,I3)') "分布荷载数NP=",NP
                WRITE(6,20)
                DO I=1,NP
                WRITE(6,30) I,MP(1,I),MP(2,I),ZP(I)
                ENDDO
                DO 10 I=1,NP
               II=2*(MP(1,I)-1)
               JJ=2*(MP(2,I)-1)
               PX=0.5*ZP(I)*T*(XY(2,MP(1,I))-XY(2,MP(2,I)))
               PY=0.5*ZP(I)*T*(XY(1,MP(2,I))-XY(1,MP(1,I)))
               B(II+1)=B(II+1)+PX
               B(II+2)=B(II+2)+PY
               B(JJ+1)=B(JJ+1)+PX
               B(JJ+2)=B(JJ+2)+PY
10                CONTINUE
20                FORMAT(1X, '序号', 4x, '首结点号', 4x, '末结点号', 7x, '荷载值')
30                FORMAT(1X, I3, 8x, I3, 9X, I3, 7X, F10.3)
                END               

                SUBROUTINE PG(NM,NE,NT,DENSITY,T,IJK,BCA,B)
                IMPLICIT REAL*8(A-H,O-Z),INTEGER(I-N)
                DIMENSION IJK(3,NE),BCA(7,NE),MG(NM),B(NT)
                IF (NM.GE.NE) THEN
                NM=NE
                DO I=1,NM
                MG(I)=I
                ENDDO
                ELSE
                READ(5,*) (MG(I),I=1,NM)
                ENDIF
                WRITE(6,'(/A,I3)') "计自重单元数NM=",NM
                WRITE(6,'(A)') "计自重的单元号为："
                K1=0
                K2=0
                DO WHILE((NM-K2).GE.10)
                WRITE(6,100) (MG(K2+J),J=1,10)
                K1=K1+1
                K2=10*K1
                ENDDO
                WRITE(6,100) (MG(K2+J),J=1,NM-K2)
                DO 10 I=1,NM
               II=2*(IJK(1,MG(I))-1)
               JJ=2*(IJK(2,MG(I))-1)
               KK=2*(IJK(3,MG(I))-1)
               PE1=-DENSITY*T*BCA(7,MG(I))/3.0
               B(II+2)=B(II+2)+PE1
               B(JJ+2)=B(JJ+2)+PE1
               B(KK+2)=B(KK+2)+PE1
10                CONTINUE
100                FORMAT(10(1X,I3))
                END               

                SUBROUTINE INPUT(NG,NE,NB,IJK,XY,MB,ZB)
                IMPLICIT REAL*8(A-H,O-Z),INTEGER(I-N)
                DIMENSION IJK(3,NE),XY(2,NG),MB(2,NB),ZB(NB)
                DO I=1,NE
                READ(5,*) K,(IJK(J,I),J=1,3)
                ENDDO
                DO I=1,NG
                READ(5,*) K,XY(1,I),XY(2,I)
                ENDDO
                DO I=1,NB
               READ(5,*) K,MB(1,I),MB(2,I),ZB(I)
                ENDDO
                WRITE(6,15)
                DO K=1,NE
                WRITE(6,20) K,(IJK(J,K),J=1,3)
                ENDDO
                WRITE(6,30)
                DO K=1,NG
                WRITE(6,35) K,XY(1,K),XY(2,K)
                ENDDO
                WRITE(6,40)
                DO K=1,NB
               WRITE(6,45) K,MB(1,K),MB(2,K),ZB(K)
                ENDDO
15                FORMAT(/1X, '单元号', 4X, 'I结点号', 4x, 'J结点号', 4x, 'K结点号')
20                FORMAT(2X, I3, 7X, I3, 8X, I3, 8X, I3)
30                FORMAT(/1X, '结点号', 8X, 'X坐标', 8X, 'Y坐标')
35                FORMAT(2X, I3, 6X, F8.3, 5X, F8.3)
40                FORMAT(/1X, '序号', 4x, '支承结点号', 6x, '支承方向', 6x, &
      & '支承位移')
45                FORMAT(2X, I2, 8X, I4, 11X, I2, 7X, F10.3)
                END

                SUBROUTINE DBC(NB,ND,NT,NX,NX1,MB,ZB,SK,B)
                IMPLICIT REAL*8(A-H,O-Z),INTEGER(I-N)
                DIMENSIONMB(2,NB),ZB(NB),SK(NT,ND),B(NT)
                DO 5 I=1,NB
               N=2*MB(1,I)-MB(2,I)
               Z=ZB(I)
               IF (ABS(Z).LT.1.0E-10) THEN
                  IF (NX.NE.NX1) GOTO 30
                  SK(N,1)=1.0
                  DO J=2,ND
                  SK(N,J)=0.0
                  ENDDO
                  DO 10 K=2,ND
                  IF (N.LT.K) GOTO 20
                  M=N-K+1
                  SK(M,K)=0.0
10                  CONTINUE
20                  CONTINUE
30                  B(N)=0.0
               ELSE
                  IF (NX.NE.NX1) GOTO 40
                  SK(N,1)=SK(N,1)*1.0E+15
40                  B(N)=SK(N,1)*Z
               ENDIF
5                CONTINUE
                END

                SUBROUTINE CALND(NG,NE,IJK,ND)
                IMPLICIT REAL*8(A-H,O-Z),INTEGER(I-N)
                DIMENSIONIJK(3,NE)
                ND=0
                DO 10 I=1,NG
               M1=0
               DO 20 J=1,NE
                  DO 30 K=1,3
                  IF(IJK(K,J).EQ.I) THEN
                  M2=0
                   DO 40 L=1,3
                  M3=IJK(L,J)-I
                  IF(M3.GT.M2) THEN
                  M2=M3
                  ENDIF
40                   CONTINUE
                  GOTO 25
                  ENDIF
30                  CONTINUE
25               IF(M2.GT.M1) M1=M2
20                CONTINUE
                IF(M1.GT.ND)ND=M1
10                CONTINUE
                ND=2*(ND+1)
                WRITE(7,100) ND
100                FORMAT(1X,'总刚矩阵的半带宽ND为：',I4)
                END

                SUBROUTINE TRES(IP,NG,NE,NX,NT,B,STR,DELD,TOD,DELST,TOST)
                IMPLICIT REAL*8(A-H,O-Z),INTEGER(I-N)
                DIMENSION B(NT),STR(3,NE),DELD(2,NG,NX),TOD(2,NG),DELST(3,NE,NX),TOST(3,NE)
                DO I=1,NG
                IO=2*(I-1)
                DELD(1,I,IP)=B(IO+1)
                DELD(2,I,IP)=B(IO+2)
                TOD(1,I)=TOD(1,I)+DELD(1,I,IP)
                TOD(2,I)=TOD(2,I)+DELD(2,I,IP)
                ENDDO
                DO I=1,NE
                DO J=1,3
                DELST(J,I,IP)=STR(J,I)
                TOST(J,I)=TOST(J,I)+DELST(J,I,IP)
                ENDDO
                ENDDO
                END


                SUBROUTINE SUPT(IP,NG,NE,NB,NX,T,V,A1,A2,IJK,MB,BCA,DELD,DELSUP,TOTSUP)
                IMPLICIT REAL*8(A-H,O-Z),INTEGER(I-N)
                DIMENSION IJK(3,NE),MB(2,NB),BCA(7,NE),DELD(2,NG,NX),EK(6,6)
                DIMENSION D1(6),DELSUP(NB,NX),TOTSUP(NB)
                DO 10 K=1,NB
                I1=MB(1,K)
                I2=MB(2,K)
                RF=0.0
               DO 20 I=1,NE
                  DO 30 J=1,3
                  I3=IJK(J,I)
                  IF (I1.NE.I3) THEN
                   GOTO 30
                  ELSE
                   I4=2*(J-1)+2-I2
                   CALL KE(I,NE,T,V,A1,A2,IJK,BCA,EK,0)
                   DO M=1,3
                   IO=IJK(M,I)
                   D1(2*(M-1)+1)=DELD(1,IO,IP)
                   D1(2*(M-1)+2)=DELD(2,IO,IP)
                   ENDDO
                   DO L=1,6
                   RF=RF+EK(I4,L)*D1(L)
                   ENDDO
                   GOTO 20
                  ENDIF
30                  CONTINUE
20               CONTINUE
                DELSUP(K,IP)=DELSUP(K,IP)+RF
                TOTSUP(K)= TOTSUP(K)+DELSUP(K,IP)
10                CONTINUE
                END               
                                 
                SUBROUTINE CLEAR1(N,B)
                IMPLICIT REAL*8(A-H,O-Z),INTEGER(I-N)
                DIMENSION B(N)
                DO I=1,N
                B(I)=0.0
                ENDDO
                END


PFEM,EXE平面问题有限元程序的数据文件输入格式：

基本输入数据文件名：*.DAT
计算结果输出数据文件名：*.OUT
中间参数输出数据文件名：*.OU1
下面按在文件中的前后顺序进行说明：
1 基本控制参数信息：NG,NE,MC,NX,NB,EO,VO,DENSITY ,T(共计5个整形数，4个实型数)
NG：结构的结点总数；
NE：结构的单元总数；
MC：平面问题的类型，MC=0，为平面应力，MC=1，为平面应变；
NX：荷载工况数；
NB：支承位移数；
EO：材料弹性模量(Pa)；
VO：材料泊松比；
DENSITY ：容重(N/m3)
T ：材料厚度(m)；
2 打印输出控制参数：NWA,NEW,NWK,NWP(4个整形数)
等于1时，输出，否则不输出。
3 单元结点信息：(K,(IJK(I,K),I=1,3),K=1,NE) (每行4个整形数，共计NE行)
K：单元号；
IJK(1,K)：K单元I结点编号；
IJK(2,K)：K单元J结点编号；
IJK(3,K)：K单元K结点编号；
4结点坐标信息：((K,XY(1,K),XY(2,K)),K=1,NG)(每行3个整形数，共计NG行)
K：结点号
XY(1,K)：K结点X坐标；
XY(2,K)：K结点Y坐标；
5 支承信息：((K,MB(1,K),MB(2,K),ZB(K)),K=1,NB)(每行3个整形数，1个实型数，共计NB行)
K：支承位移序号；
MB(1,K)：第K个支承位移所在的结点号；
MB(2,K)：第K个支承位移的坐标方向；
ZB(K)：第K个支承位移的数值；
6 按NX荷载工况数输入荷载信息：每一荷载工况如下
(1) NF,NP,NM(3个整型数)
   NF：集中荷载个数；
   NP：分布荷载个数；
   NM：计自重单元数；
(2) 若NF≠0，则输入下面数据
      K,MF(1,K),MF(2,K),ZF(K)(每行3个整形数，1个实型数，共计NF行)
       K：集中荷载序号；
       MF(1,K)：第K个集中荷载作用的结点号；
       MF(2,K)：第K个集中荷载的坐标方向；
       ZF(K)：第K个集中荷载的数值；
(3) 若NP≠0，则输入下面数据
      K,MP(1,K),MP(2,K),ZP(K)(每行3个整形数，1个实型数，共计NP行)
       K：分布荷载序号；
       MP(1,K)：第K个分布荷载作用的结点号；
       MP(2,K)：第K个分布荷载的坐标方向；
       ZP(K)：第K个分布荷载的数值；
(4) 若NM≠0，则输入下面数据
      若NM≥NE，则表示计所有单元的自重，不需输入计自重的单元号；
      若NM<NE，则需要输入计自重的单元号；

例题验证

1 (选自徐芝伦编《弹性力学简明教程》(第二版))设有对角受压的正方形薄板，荷载沿厚度均匀分布，为2N/m，由于xz和yz面均为该薄板的对称面，所以只需取1/4部分作为计算对象，图1(a,b)，将该对象划分为4个单元，共有6个结点。对称面上的结点没有垂直于对称面的位移分量，因此，在1、2、4三个结点上设置了水平链杆支座，在4、5、6三个结点设置了铅直链杆支座，如图1c所示。

查看完整版本: 平面弹性力学有限元源程序（FORTRAN）