有限元分析程序与算例
平面四结点等参数单元分析程序
一、 源程序：
! PROGRAM NODE4.FOR
! IMPLICIT DOUBLE PRECISION(A-H,O-Z) DIMENSION A(115000) INTEGER IA(115000) EQUIVALENCE(IA,A) CHARACTER*12 AR,BR;
COMMON /C1/ NP,NE,NM,NR,NW,NF,NI,NDP COMMON /C2/ N,MX,NH
COMMON /C3/ SKE(8,8),NN(8),RF(8),B(3,8),XY(2,4) WRITE(0,250)
250 FORMAT(///'PLEASE INPUT FILE NAME OF DATA=') READ(*,400)AR 400 FORMAT(A12)
OPEN(15,FILE=AR,STATUS='OLD') WRITE(*,350)
350 FORMAT(///'PLEASE OUTPUT FILE NAME OF DATA=') READ(*,400)BR
OPEN(16,FILE=BR,STATUS='UNKNOWN')
CALL CONTROL NBY=1
N1=1
N2=N1+NP*2 N3=N2+NP*2 N4=N3+NP*3
N5=N4+NP*2*NBY N6=N5+NP*NBY*2 N7=N6+4*NM*NBY N8=N7+2*NW*NBY N9=N8+NP N10=N9+NDP
N11=N10+6*NDP*NBY NSPACE=115000-N11
!C --N1=MA N2=JR N3=NRR N4=COOR N5=R N6=AE !C --N7=WG N8=IV N9=NDP N10=DV N11=SK WRITE(16,800)NSPACE
- 1 -
有限元分析程序与算例
800 FORMAT(40X,'NSPACE=',I8)
CALL INPUT(IA(N1),IA(N2),A(N4),A(N5),A(N6),A(N7),IA(N3))
OPEN(12,FILE='ELEMENT',STATUS='UNKNOWN',FORM='UNFORMATTED') OPEN(7,FILE='LOAD',STATUS='UNKNOWN',FORM='UNFORMATTED') OPEN(10,FILE='GRAPH1',STATUS='UNKNOWN') OPEN(11,FILE='GRAPH2',STATUS='UNKNOWN') OPEN(13,FILE='GRAPH3',STATUS='UNKNOWN') DO 10 IE=1,NE
CALL SDTK4(IA(N1),IA(N2),A(N4),A(N6),A(N7)) CALL ASLOAD(A(N5),8)
10 CONTINUE CALL CBAND(IA(N1))
CALL OUTPUT(IA(N2),A(N5),0) CALL ASESK(A(N11),IA(N1))
IF(NDP.GT.0)CALL TREAT(A(N11),IA(N1),A(N5),IA(N2),IA(N9),A(N10)) CALL DECOMP(A(N11),IA(N1)) CALL FOBA(A(N11),IA(N1),A(N5)) CALL OUTPUT(IA(N2),A(N5),1)
CALL STRESS(A(N6),A(N5),IA(N1),A(N9)) CLOSE(12,STATUS='DELETE') WRITE(16,700)
700 FORMAT(2X,'PROGRAM HAS BEEN ENDED') STOP END
!C****************************************************************** SUBROUTINE CONTROL
!C IMPLICIT DOUBLE PRECISION(A-H,O-Z) COMMON/C1/NP,NE,NM,NR,NW,NI,NF,NDP COMMON/C2/N,MX,NH
READ(15,*)NP,NE,NM,NR,NW,NI,NF,NDP WRITE(16,600)NP,NE,NM,NR,NW,NI,NF,NDP
600 FORMAT('NUMBER OF NODE………………………………………………NP=',I5,/ & 1X,'NUMBER OF ELEMENT …………………………………………NE=',I5,/& 1X,'NUMBER OF MATERIAL…………………………………………NM=',I5,/ & 1X,'NUMBER OF CONSTAINT ………………………………………NR=',I5,/ & 1X,'NUMBER OF WATER PRESS KIND………………………………NW=',I5,/ & 1X,'NUMBER OF CONCENTRATE LOAD………………………………NF=',I5,/ & 1X,'PLANE STRESS OR PLANE STAIN ……………………………NI=',I5,/ &
1X,'NUMBER OF KNOWN-DISPLACEMENT……………………………NDP=',I5) RETURN END
- 2 -
有限元分析程序与算例
!C****************************************************************** SUBROUTINE INPUT(MA,JR,COOR,R,AE,WG,NRR) !C IMPLICIT DOUBLE PRECISION(A-H,O-Z)
DIMENSION JR(2,*),COOR(2,*),AE(4,*),NRR(3,*),XY(2),IR(2),WG(2,*),R(*),MA(*) COMMON /C1/ NP,NE,NM,NR,NW,NF,NI,NDP COMMON /C2/ N,MX,NH
DO 10 I=1,NP DO 10 J=1,2 10 JR(J,I)=1 WRITE(16,500)
DO 20 I=1,NR
READ(15,*) NN,IR NRR(1,I)=NN NRR(2,I)=IR(1) NRR(3,I)=IR(2) DO 15 J=1,2
JR(J,NN)=IR(J) 15 CONTINUE
20 WRITE(16,600) ((NRR(I,J),I=1,3),J=I,I) N=0
DO 30 I=1,NP DO 30 J=1,2
IF(JR(J,I))30,30,25 25 N=N+1 JR(J,I)=N
30 CONTINUE
DO 40 I=1,N MA(I)=0
40 CONTINUE
DO 70 I=1,NP
READ(15,*) NN,XY
IF(NN.NE.I) GOTO 60 DO 45 J=1,2
45 COOR(J,NN)=XY(J) GOTO 70
60 WRITE(16,750) NN,I !C STOP 333 70 CONTINUE
- 3 -
有限元分析程序与算例
WRITE(16,800) (NN,(JR(I,NN),I=1,2), (COOR(J,NN),J=1,2),NN=1,NP) READ(15,*) ((AE(I,J),I=1,4),J=1,NM)
WRITE(16,910) (J,(AE(I,J),I=1,4),J=1,NM) IF(NI.EQ.0) GOTO 75 DO 73 J=1,NM
AE(2,J)=AE(2,J)/(1.+AE(2,J))
AE(1,J)=AE(1,J)*(1.-AE(2,J)*AE(2,J)) 73 CONTINUE 75 CONTINUE IF(NW.EQ.0) GOTO 80
READ(15,*) ((WG (I,J),I=1,2),J=1,NW) WRITE(16,960) (J,(WG(I,J),I=1,2),J=1,NW) 80 DO 85 I=1,N R(I)=0.0
85 CONTINUE IF(NF.EQ.0) GOTO 200 WRITE(16,980)
DO 100 I=1,NF
READ(15,*) NN,XY WRITE(16,990) NN,XY DO 95 J=1,2 L=JR(J,NN)
IF(L.EQ.0) GOTO 95 R(L)=R(L)+XY(J) 95 CONTINUE 100 CONTINUE
200 CONTINUE
500 FORMAT(/25X,'NODAL INFORMATION'/10X,'CONSTRINED& MESSAGE NODE NO.STATE'/) 600 FORMAT(6X,8(I5,2X,2I1))
750 FORMAT(1X,'***FATAL ERROR ***',/,'CARDS',& 'INPUT',I5,'IS NOT EQUAL TO',I5)
800 FORMAT(1X,'NODE NO.',2X,'X-FREEDOM',2X,&
'Y-FREEDOM',2X,'X-COORDINAT',2X,'Y-COORDINAT', & /(1X,I5,2I10,4X,2F14.4))
910 FORMAT(/20X,'MATERIAL PROPERTIES',/,2X, &
'N.M.',5X,'YOUNGS MODULUS',5X,'POISION RATIO', & 4X,'UNIT WEIGHT WIDTH'/(1X,I5,4E16.4))
960 FORMAT(/5X,'PARAMETERES OF WATER AND', & 'SILT PRESSURE'/2X,'N.P.',2X,'ZERO-PRESSURE', & 'SURFACE',8X,'UNIT WEIGHT'/(1X,I5,2F15.5)) 980 FORMAT(/20X,'CONCENTRAED FORCES'/1X,&
- 4 -
有限元分析程序与算例
'NODE NO.',8X,'X-DIRECTION',8X,'Y-DIRECTION'/) 990 FORMAT(1X,I5,2F20.6) RETURN END
!C****************************************************************** SUBROUTINE CBAND(MA) !C $LARGE:MA
!C IMPLICIT DOUBLE PRECISION(A-H,O-Z) DIMENSION MA(*)
COMMON /C2/ N,MX,NH MX=0 MA(1)=1 DO 10 I=2,N
IF(MA(I).GT.MX) MX=MA(I) 10 MA(I)=MA(I)+MA(I-1) CONTINUE NH=MA(N)
WRITE(16,500) N,MX,NH
500 FORMAT(4X,'TOTAL DEGREES OF FREEDOM.…………………N=',I6/ & 4X,'MAX-SEMI-BANDWIDTH ………………………… MX=',I6/ & 4X,'TOTAL-STORAGE …………………………………… NH=',I6) RETURN END
!C****************************************************************** SUBROUTINE SDTK4(MA,JR,COOR,AE,WG)
!C IMPLICIT DOUBLE PRECISION(A-H,O-Z)
DIMENSION MA(*),JR(2,*),COOR(2,*),AE(4,*),WG(2,*),NOD(4),NFACE(4),PR(2) COMMON /C1/ NP,NE,NM,NR,NW,NF,NI,NDP COMMON /C2/N,MX,NH
COMMON /C3/ SKE(8,8),NN(8),RF(8),B(3,8),XY(2,4) COMMON /C4/ NEE,NME,NET,NK,NSF,NST COMMON /GAUSS/ RSTG(2)
RSTG(1)=-0.5773503 RSTG(2)=0.5773503
READ(15,*) NEE,NOD,NME,NET,NK,NSF,NST WRITE(16,600) NEE,NME,NET,NK,NSF,NST,NOD DO 10 I=1,4 J=NOD(I) DO 10 L=1,2
XY(L,I)=COOR(L,J) 10 CONTINUE
DO 40 I=1,8
- 5 -