有限元程序分析与算例

有限元分析程序与算例

平面四结点等参数单元分析程序

一、 源程序：

! 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 -

有限元分析程序与算例

RF(I)=0.0 DO 30 J=1,8 SKE(I,J)=0.0

30 CONTINUE 40 CONTINUE WRITE(*,*)'NE=',NEE DO 50 I=1,4 JB=NOD(I) DO 50 M=1,2 JJ=2*(I-1)+M NN(JJ)=JR(M,JB) 50 CONTINUE L=N

DO 80 I=1,8

IF(NN(I).EQ.0) GOTO 80 IF(NN(I).LT.L) L=NN(I) 80 CONTINUE

DO 85 M=1,8 JP=NN(M)

IF(JP.EQ.0) GOTO 85

IF(JP-L+1.GT.MA(JP)) MA(JP)=JP-L+1 85 CONTINUE

CALL STIF(AE)

IF(NK.EQ.0) GOTO 100 DO 95 IW=1,NK

READ(15,*) NNN,NFACE WRITE(16,650) NNN,NFACE Y0=WG(1,NNN) GAMA=WG(2,NNN) DO 90 JW=1,4 ND=NFACE(JW)

IF(ND.EQ.0) GOTO 90

CALL SURFOR (JW,AE,PR,Y0,GAMA,1) 90 CONTINUE 95 CONTINUE 100 CONTINUE

IF(NSF.EQ.0) GOTO 200 DO 110 IW=1,NSF READ(15,*) ND,PR WRITE(16,750) ND,PR

CALL SURFOR(ND,AE,PR,Y0,GAMA,0)

- 6 -

有限元分析程序与算例

110 CONTINUE 200 CONTINUE

IF(NST.EQ.0) GOTO 300 DO 210 IW=1,NST READ(15,*) ND,PR WRITE(16,750) ND,PR

CALL SURNST (ND,AE,PR) 210 CONTINUE 300 CONTINUE

!C WRITE (16,800) (RF(I),I=1,8) WRITE (12) NN,SKE

WRITE (7) NEE,NME,NOD,NN,XY

600 FORMAT(3X,' NEE= NME= NET= NK= NSF= NST= NOD='/3X,6I5,4I5) 650 FORMAT(3X,' SURFACE NEWS NNN=',I3,2X,4I1)

750 FORMAT(3X,' SURFACE NEWS ND=',I3,5X, 'PR=',2F12.2) 800 FORMAT(3X,'RF=',/(1X,8F10.4)) RETURN END

!C****************************************************************** SUBROUTINE STIF(AE)

!C IMPLICIT DOUBLE PRECISION(A-H,O-Z) DIMENSION AE(4,*),RJX(2,2),BV(8),D(4) 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 /C5/ FUN(4),P(2,4),XJR(2,2) COMMON /GAUSS/ RSTG(2)

E=AE(1,NME) U=AE(2,NME)

GAMMA=AE(3,NME) T=AE(4,NME)

D1=E*(1.0-U)/((1.0+U)*(1.0-2.0*U)) D2=E*U/((1.0+U)*(1.0-2.0*U)) D3=E*0.5/(1.0+U)

DO 99 IS=1,2 S=RSTG(IS) DO 98 IR=1,2

- 7 -

有限元分析程序与算例

R=RSTG(IR)

CALL RMSD (DET,R,S,RJX) IF(NET.EQ.0) GOTO 40 DO 30 I=1,4 J=2*I

RF(J)=RF(J)-FUN(I)*DET*GAMMA*T 30 CONTINUE 40 K2=0 DO 50 I=1,4 K2=K2+2 K1=K2-1

BV(K1)=B(1,K1) BV(K2)=B(2,K2) 50 CONTINUE DO 70 I=1,8 DO 60 J=I,8

SKE(I,J)=SKE(I,J)+BV(I)*BV(J)*DET*T 60 CONTINUE 70 CONTINUE 98 CONTINUE 99 CONTINUE

DO 104 I=2,8 M=I-1

DO 102 J=1,M SKE(I,J)=SKE(J,I) 102 CONTINUE 104 CONTINUE KK=-2

DO 120 I=1,4 KK=KK+2 K1=KK+1 K2=K1+1 DO 115 J=1,4 LL=(J-1)*2 L1=LL+1 L2=L1+1 IC=0

DO 110 K=1,2 M=K+KK DO 105 L=1,2 NL=L+LL IC=IC+1

- 8 -

有限元分析程序与算例

D(IC)=SKE(M,NL) 105 CONTINUE 110 CONTINUE

SKE(K1,L1)=D(1)*D1+D(4)*D3 SKE(K2,L2)=D(4)*D1+D(1)*D3 SKE(K1,L2)=D(2)*D2+D(3)*D3 SKE(K2,L1)=D(3)*D2+D(2)*D3 115 CONTINUE 120 CONTINUE

!C WRITE (16,500) ((SKE(I,J),I=1,4),J=1,4) 500 FORMAT (10X,'SKE='/(1X,8F10.2)) RETURN END

!C****************************************************************** SUBROUTINE RMSD(DET,R,S,RJX)

!C IMPLICIT DOUBLE PRECISION(A-H,O-Z) DIMENSION RJX(2,2)

COMMON /C4/ NEE,NME,NET,NK,NSF,NST COMMON /C5/ FUN(4),P(2,4),XJR(2,2)

COMMON /C3/ SKE(8,8),NN(8),RF(8),B(3,8),XY(2,4)

CALL FPJD(R,S,DET)

REC=1.0/DET

RJX(1,1)=REC*XJR(2,2) RJX(2,1)=-REC*XJR(2,1) RJX(1,2)=-REC*XJR(1,2) RJX(2,2)=REC*XJR(1,1)

K2=0

DO 30 K=1,4 K2=K2+2 K1=K2-1 DO 20 L=1,2 B(L,K1)=0.0 B(L,K2)=0.0

20 CONTINUE DO 25 I=1,2

B(1,K1)=B(1,K1)+RJX(1,I)*P(I,K) B(2,K2)=B(2,K2)+RJX(2,I)*P(I,K) 25 CONTINUE B(3,K1)=B(2,K2) B(3,K2)=B(1,K1)

- 9 -

有限元分析程序与算例

30 CONTINUE

RETURN END

!C****************************************************************** SUBROUTINE FPJD(R,S,DET)

!C IMPLICIT DOUBLE PRECISION(A-H,O-Z) DIMENSION XI(4),ETA(4)

COMMON /C4/ NEE,NME,NET,NK,NSF,NST COMMON /C5/ FUN(4),P(2,4),XJR(2,2)

COMMON /C3/ SKE(8,8),NN(8),RF(8),B(3,8),XY(2,4) DATA XI/-1.0,1.0,1.0,-1./ DATA ETA/-1.0,-1.0,1.0,1.0/ DO 10 I=1,4

FUN(I)=0.25*(1.0+XI(I)*R)*(1.0+ETA(I)*S) 10 CONTINUE

DO 50 I=1,4

P(1,I)=0.25*XI(I)*(1.0+ETA(I)*S) P(2,I)=0.25*ETA(I)*(1.0+XI(I)*R) 50 CONTINUE

DO 80 I=1,2 DO 75 J=1,2 DET=0.

DO 70 K=1,4

DET=DET+P(I,K)*XY(J,K) 70 CONTINUE XJR(I,J)=DET

75 CONTINUE 80 CONTINUE

DET=XJR(1,1)*XJR(2,2)-XJR(1,2)*XJR(2,1) !C WRITE (0,500) DET IF (DET.LT.1.0D-5) GOTO 100 RETURN

500 FORMAT (10X,'DET=',F15.3) 100 WRITE(16,600) NEE,R,S,DET

600 FORMAT(1X,'ERROR OF NEGTIVE OR ZERO', & 'JACOBIAN DETERMINANT CONPUTED FOR',& 'ELEMENT (',I5,')'/5X,'R=S=DET=',3F10.5) RETURN END

!C******************************************************************

- 10 -

有限元分析程序与算例

SUBROUTINE SURFOR(ND,AE,PR,Y0,GAMA,NSI) !C IMPLICIT DOUBLE PRECISION(A-H,O-Z) DIMENSION PR(2),KCRD(4),RST(2),AE(4,*), &

KFACE(2,4),IPRM(2),FVAL(4),NODES(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 /C5/ FUN(4),P(2,4),XJR(2,2) COMMON /GAUSS/ RSTG(2) DATA KCRD/1,1,2,2/

DATA KFACE/2,3,1,4,3,4,1,2/ DATA IPRM/2,1/

DATA FVAL/1.0,-1.0,1.0,-1./

T=AE(4,NME) FACT=-FVAL(ND)

DO 20 I=1,2 J=KFACE(I,ND) NODES(I)=J

IF (NSI.EQ.0) GOTO 20 Y=Y0-XY(2,J) PR(I)=0.0

IF (Y.GT.0.0) PR(I)=Y*GAMA 20 CONTINUE

ML=KCRD(ND) MM=IPRM(ML) RST(ML)=FVAL(ND) DO 70 LX=1,2

RST(MM)=RSTG(LX)

CALL FPJD(RST(1),RST(2),DET) PXY=0.0 DO 25 I=1,2 J=NODES(I)

PXY=PXY+FUN(J)*PR(I) 25 CONTINUE

A1=XJR(MM,2)*(-1)**MM A2=-XJR(MM,1)*(-1)**MM

DO 60 I=1,2 J=NODES(I)

- 11 -

有限元分析程序与算例

K2=2*J K1=K2-1

Q=PXY*FUN(J)*FACT*T RF(K1)=RF(K1)+Q*A1 RF(K2)=RF(K2)+Q*A2 60 CONTINUE 70 CONTINUE RETURN END

!C****************************************************************** SUBROUTINE SURNST(ND,AE,PR)

DIMENSION PR(2),KCRD(4),RST(2),AE(4,*),KFACE(2,4), & IPRM(2),FVAL(4),NODES(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 /C5/ FUN(4),P(2,4),XJR(2,2) COMMON /GAUSS/ RSTG(2) DATA KCRD/1,1,2,2/

DATA KFACE/2,3,1,4,3,4,1,2/ DATA IPRM/2,1/

DATA FVAL/1.0,-1.0,1.0,-1.0/ T=AE(4,NME) DO 20 I=1,2

20 NODES(I)=J ML=KCRD(ND) MM=IPRM(ML)

RST(ML)=FVAL(LX) DO 70 LX=1,2

RST(MM)=RSTG(LX)

CALL FPJD(RST(1),RST(2),DET) RXY=0.0

DO 25 I=1,2 J=NODES(I)

PXY=PXY+FUN(J)*PR(I) 25 CONTINUE

A2=XJR(MM,2)*(-1)**MM A1=-XJR(MM,1)*(-1)**MM DO 60 I=1,2 J=NODES(I) K2=2*J K1=K2-1

- 12 -

有限元分析程序与算例

Q=PXY*FUN(J)*T RF(K1)=RF(K1)+Q*A1 RF(K2)=RF(K2)+Q*A2 60 CONTINUE 70 CONTINUE RETURN END

!C****************************************************************** SUBROUTINE ASLOAD (R,LP)

!C IMPLICIT DOUBLEPRECISION (A-H,O-Z) DIMENSION R(*)

COMMON /C1/ NP,NE,NM,NR,NW,NF,NI,NDP

COMMON /C3/ SKE(8,8),NN(8),RF(8),B(3,8),XY(2,4) DO 20 I=1,LP L=NN(I)

IF(L.EQ.0) GOTO 20 R(L)=R(L)+RF(I) 20 CONTINUE RETURN END

!C****************************************************************** SUBROUTINE OUTPUT (JR,F,NS)

!C IMPLICIT DOUBLEPRECISION (A-H,O-Z) DIMENSION JR(2,*),F(*),B(2)

COMMON /C1/ NP,NE,NM,NR,NW,NF,NI,NDP If(ns.eq.0)write(16,300) If(ns.eq.1)write(16,400) WRITE (16,500) DO 40 I=1,NP DO 20 J=1,2 20 B(J)=0. DO 30 J=1,2 L=JR(J,I)

IF(L.EQ.0) GOTO 30 B(J)=F(L)

30 CONTINUE WRITE (16,550) I,B

IF(NS.EQ.1) WRITE(10,650) I,B 40 CONTINUE

300 FORMAT (/28X,'LOAD OF NODES'/)

400 FORMAT (/25X,'DISPLACEMENT OF NODES'/)

500 FORMAT (10X,'NODE NO.',8X,'X-DIRECTION',8X,'Y-DIRECTION'/) 550 FORMAT (13X,I5,2(8X,E12.4)) 650 FORMAT (3X,I4,2E12.4)

- 13 -

有限元分析程序与算例

RETURN END

!C****************************************************************** SUBROUTINE ASESK(SK,MA)

!C IMPLICIT DOUBLEPRECISION (A-H,O-Z) DIMENSION SK(*),MA(*)

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) DO 10 I=1,NH SK(I)=0.0

10 CONTINUE REWIND 12

DO 95 IE=1,NE READ (12) NN,SKE DO 50 JJ=1,8 I2=NN(JJ)

IF (I2.GT.0.AND.I2.LE.N) GOTO 60 50 CONTINUE GOTO 95 60 DO 70 L=1,8 DO 70 M=1,8 JJ=NN(L) JK=NN(M)

IF(JJ.GE.JK.AND.JK.GT.0.AND.JJ.GT.0.AND.JJ.LE.N) GOTO 65 GOTO 70

65 J3=MA(JJ)-JJ+JK SK(J3)=SK(J3)+SKE(L,M) 70 CONTINUE 95 CONTINUE RETURN END

!C****************************************************************** SUBROUTINE TREAT(SK,MA,R,JR,NDI,DV)

!C IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION SK(*),MA(*),R(*),JR(2,*),NDI(*),DV(2,*) COMMON /C1/ NP,NE,NM,NR,NW,NF,NI,NDP WRITE(16,500) DO 5 I=1,NDP

READ(15,*)NDI(I),(DV(J,I),J=1,2) WRITE(16,600)NDI(I),(DV(J,I),J=1,2) 5 CONTINUE

DO 30 I=1,NDP

- 14 -

有限元分析程序与算例

JJ=NDI(I) DO 20 J=1,2 L=JR(J,JJ) JN=MA(L)

IF(DV(J,I).EQ.0)GOTO 20 SK(JN)=1.E30

R(L)=DV(J,I)*1.E30 20 CONTINUE 30 CONTINUE

500 FORMAT(/1X,'KNOWN DISPLACEMENT NODES AND (X,Y) VALUES') 600 FORMAT(10X,I8,10X,F10.6,10X,F10.6) RETURN END

!C****************************************************************** SUBROUTINE DECOMP(SK,MA)

!C IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION SK(*),MA(*)

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) DO 65 I=2,N

L=I-MA(I)+MA(I-1)+1 K=I-1 L1=L+1

IF (L1.GT.K) GOTO 55 DO 50 J1=L1,K IJ=MA(I)-I+J1

M=MA(J1-1)-MA(J1)+J1+1 IF(L.GT.M) M=L MP=J1-1

IF(M.GT.MP) GOTO 50 DO 40 LP=M,MP IP=MA(I)-I+LP JP=MA(J1)-J1+LP

SK(IJ)=SK(IJ)-SK(IP)*SK(JP) 40 CONTINUE 50 CONTINUE

55 IF(L.GT.K) GOTO 65 DO 60 LP=L,K IP=MA(I)-I+LP LPP=MA(LP)

SK(IP)=SK(IP)/SK(LPP)

- 15 -

有限元分析程序与算例

II=MA(I)

SK(II)=SK(II)-SK(IP)*SK(IP)*SK(LPP) 60 CONTINUE 65 CONTINUE

RETURN END

!C****************************************************************** SUBROUTINE FOBA(SK,MA,R)

!C IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION SK(*),MA(*),R(*)

COMMON /C1/ NP,NE,NM,NR,NW,NF,NI,NDP COMMON /C2/N,MX,NH

DO 10 I=2,N

L=I-MA(I)+MA(I-1)+1 K=I-1

IF(L.GT.K) GOTO 10 DO 5 LP=L,K IP=MA(I)-I+LP

R(I)=R(I)-SK(IP)*R(LP) 5 CONTINUE 10 CONTINUE DO 20 I=1,N II=MA(I)

R(I)=R(I)/SK(II) 20 CONTINUE DO 30 J1=2,N I=2+N-J1

L=I-MA(I)+MA(I-1)+1 K=I-1

IF(L.GT.K) GOTO 30 DO 40 J=L,K IJ=MA(I)-I+J

R(J)=R(J)-SK(IJ)*R(I) 40 CONTINUE 30 CONTINUE RETURN END

!C****************************************************************** SUBROUTINE STRESS (AE,R,ND,SNOD) !C IMPLICIT REAL*8(A-H,O-Z)

DIMENSION ME(4),XYG(2),AE(4,*),RJX(2,2),SG(6),S(3,8),ND(*),R(*),SNOD(6,*)

- 16 -

有限元分析程序与算例

COMMON /C1/ NP,NE,NM,NR,NW,NF,NI,NDP

COMMON /C3/SKE(8,8),NN(8),RF(8),B(3,8),XY(2,4) COMMON /C5/FUN(4),P(2,4),XJR(2,2) DO 1 I=1,NP 1 ND(I)=0 DO 5 I=1,NP DO 5 J=1,6

5 SNOD(J,I)=0.0 WRITE(16,300) REWIND 7

DO 105 IE=1,NE

READ(7) NEE,NME,ME,NN,XY E=AE(1,NME) U=AE(2,NME)

D1=E*(1.-U)/((1.+U)*(1.-2.*U)) D2=E*U/((1.+U)*(1.-2.*U)) D3=.5*E/(1.+U) DO 10 ID=1,4 LD=ME(ID)

10 ND(LD)=ND(LD)+1 SS=0. TT=0.

CALL RMSD (DET,SS,TT,RJX) J2=0

DO 30 I=1,4 J1=J2+1 J2=J1+1 BI=B(1,J1) CI=B(2,J2)

S(1,J1)=D1*BI S(2,J1)=D2*BI S(3,J1)=D3*CI S(1,J2)=D2*CI S(2,J2)=D1*CI S(3,J2)=D3*BI 30 CONTINUE DO 20 J=1,2 XYG(J)=0. DO 20 K=1,4

XYG(J)=XYG(J)+FUN(K)*XY(J,K) 20 CONTINUE DO 55 II=1,6 55 SG(II)=0. DO 70 K=1,8

- 17 -

有限元分析程序与算例

NA=NN(K)

IF (NA.EQ.0) GOTO 70 DO 60 L=1,3

SG(L)=SG(L)+S(L,K)*R(NA) 60 CONTINUE 70 CONTINUE CALL PRI(SG)

WRITE(16,500) NEE,(SG(I),I=1,6),(XYG(J),J=1,2) WRITE(13,550) IE,SG,XYG DO 80 JJ=1,4 LD=ME(JJ) DO 80 II=1,3

SNOD (II,LD)=SNOD(II,LD)+SG(II) 80 CONTINUE 105 CONTINUE

DO 120 I=1,NP LL=ND(I)

IF(LL.EQ.0) GOTO 120 DO 110 K=1,3

SNOD(K,I)=SNOD(K,I)/LL 110 CONTINUE 120 CONTINUE

WRITE(16,700) DO 240 J=1,NP DO 215 K=1,3 SG(K)=SNOD(K,J) 215 CONTINUE CALL PRI (SG)

WRITE (16,660) J,(SG(K),K=1,6) WRITE (11,550) J,SG 240 CONTINUE

300 FORMAT (/25X,'ELEMENT MID-POINT STRESS'/'ELE NO.', & 'SIG--XX',2X,'SIG--YY',2X,'SIG--XY',3X,'SIG-- & 1',3X, 'SIG--2 ANGLE',3X,'X--CO.',2X,'Y--CO.'/) 500 FORMAT (2X,I4,6F9.2,2X,2F7.2) 550 FORMAT (1X,I4,6F9.2,2X,2F7.2)

700 FORMAT (/25X,'NODES & PRINCIPAL STRESS'/2X,'NODE NO.',& 4X,'SIG-XX',5X,'SIG-YY',5X,'SIG-XY',6X,'SIG-1',6X,& 'SIG-2',6X,'ANGLE'/) 660 FORMAT (3X,I5,1X,6F11.3) RETURN

- 18 -

有限元分析程序与算例

END

!C****************************************************************** SUBROUTINE PRI(A)

!C IMPLICIT REAL*8(A-H,O-Z) DIMENSION A(6)

H1=A(1)+A(2) H2=A(1)-A(2)

H2=SQRT(H2*H2+4. *A(3)*A(3)) A(4)=(H1+H2)/2. A(5)=(H1-H2)/2.

IF (ABS(A(3)).GT.1.0E-4) GOTO 230 IF (A(1).GT.A(2)) GOTO 220 A(6)=90. GOTO 250

220 A(6)=0. GOTO 250

230 A(6)=ATAN((A(4)-A(1))/A(3))*57.2958 250 RETURN END

!C----------------NODE4.FOR PROGRAM HAS BEEN END ED

二、数据输入：

15 8 1 4 1 0 1 0 1 0 1 6 0 1 11 0 1 15 1 0 1 0.0 1.0 2 2.0 1.0 3 4.0 1.0 4 6.0 1.0 5 8.0 1.0 6 0.0 0.5 7 2.0 0.5 8 4.0 0.5 9 6.0 0.5 10 8.0 0.5 11 0.0 0.0 12 2.0 0.0 13 4.0 0.0 14 6.0 0.0 15 8.0 0.0

- 19 -

有限元分析程序与算例

20000000.0 0.2 0.0 1.0 11.0 1.0

1 6 7 2 1 1 0 1 0 0 1 0 0 1 0

2 7 8 3 2 1 0 1 0 0 1 0 0 1 0

3 8 9 4 3 1 0 1 0 0 1 0 0 1 0

4 9 10 5 4 1 0 1 0 0 1 0 0 1 0

5 11 12 7 6 1 0 0 0 0 6 12 13 8 7 1 0 0 0 0 7 13 14 9 8 1 0 0 0 0 8 14 15 10 9 1 0 0 0 0

三、结果输出：

NUMBER OF NODE………………………………………………NP= 15 NUMBER OF ELEMENT …………………………………………NE= 8 NUMBER OF MATERIAL…………………………………………NM= 1 NUMBER OF CONSTAINT ………………………………………NR= 4 NUMBER OF WATER PRESS KIND………………………………NW= 1 NUMBER OF CONCENTRATE LOAD………………………………NF= 0 PLANE STRESS OR PLANE STAIN ……………………………NI= 1

NUMBER OF KNOWN-DISPLACEMENT……………………………NDP= 0 NSPACE= 114813

NODAL INFORMATION CONSTRINED MESSAGE NODE NO.STATE

1 01 6 01 11 01 15 10

NODE NO. X-FREEDOM Y-FREEDOM X-COORDINAT Y-COORDINAT 1 0 1 .0000 1.0000 2 2 3 2.0000 1.0000 3 4 5 4.0000 1.0000 4 6 7 6.0000 1.0000 5 8 9 8.0000 1.0000 6 0 10 .0000 .5000 7 11 12 2.0000 .5000 8 13 14 4.0000 .5000 9 15 16 6.0000 .5000

- 20 -

有限元分析程序与算例

10 17 18 8.0000 .5000 11 0 19 .0000 .0000 12 20 21 2.0000 .0000 13 22 23 4.0000 .0000 14 24 25 6.0000 .0000 15 26 0 8.0000 .0000

MATERIAL PROPERTIES

N.M. YOUNGS MODULUS POISION RATIO UNIT WEIGHT WIDTH 1 .2000E+08 .2000E+00 .0000E+00 .1000E+01

PARAMETERES OF WATER ANDSILT PRESSURE

N.P. ZERO-PRESSURESURFACE UNIT WEIGHT 1 11.00000 1.00000

NEE= NME= NET= NK= NSF= NST= NOD=

1 1 0 1 0 0 6 7 2 1 SURFACE NEWS NNN= 1 0010

NEE= NME= NET= NK= NSF= NST= NOD=

2 1 0 1 0 0 7 8 3 2 SURFACE NEWS NNN= 1 0010

NEE= NME= NET= NK= NSF= NST= NOD=

3 1 0 1 0 0 8 9 4 3 SURFACE NEWS NNN= 1 0010

NEE= NME= NET= NK= NSF= NST= NOD=

4 1 0 1 0 0 9 10 5 4 SURFACE NEWS NNN= 1 0010

NEE= NME= NET= NK= NSF= NST= NOD=

5 1 0 0 0 0 11 12 7 6 NEE= NME= NET= NK= NSF= NST= NOD=

6 1 0 0 0 0 12 13 8 7 NEE= NME= NET= NK= NSF= NST= NOD=

7 1 0 0 0 0 13 14 9 8 NEE= NME= NET= NK= NSF= NST= NOD=

8 1 0 0 0 0 14 15 10 9 TOTAL DEGREES OF FREEDOM.…………………N= 26 MAX-SEMI-BANDWIDTH ………………………… MX= 13 TOTAL-STORAGE …………………………………… NH= 230

LOAD OF NODES

NODE NO. X-DIRECTION Y-DIRECTION

1 .0000E+00 -.1000E+02 2 .0000E+00 -.2000E+02

- 21 -

有限元分析程序与算例

3 .0000E+00 -.2000E+02 4 .0000E+00 -.2000E+02 5 .0000E+00 -.1000E+02 6 .0000E+00 .0000E+00 7 .0000E+00 .0000E+00 8 .0000E+00 .0000E+00 9 .0000E+00 .0000E+00 10 .0000E+00 .0000E+00 11 .0000E+00 .0000E+00 12 .0000E+00 .0000E+00 13 .0000E+00 .0000E+00 14 .0000E+00 .0000E+00 15 .0000E+00 .0000E+00

DISPLACEMENT OF NODES

NODE NO. X-DIRECTION Y-DIRECTION

1 .0000E+00 -.1905E-02 2 -.6940E-04 -.1764E-02 3 -.1299E-03 -.1358E-02 4 -.1723E-03 -.7407E-03 5 -.1875E-03 -.3574E-05 6 .0000E+00 -.1907E-02 7 .8459E-07 -.1766E-02 8 .2524E-06 -.1359E-02 9 .1587E-06 -.7414E-03 10 .6309E-06 -.2797E-05 11 .0000E+00 -.1905E-02 12 .6961E-04 -.1764E-02 13 .1302E-03 -.1358E-02 14 .1728E-03 -.7412E-03 15 .1894E-03 .0000E+00

ELEMENT MID-POINT STRESS

ELE NO.SIG--XX SIG--YY SIG--XY SIG-- 1 SIG--2 ANGLE X--CO. Y--CO.

1 -348.24 -8.26 10.23 -7.95 -348.55 88.28 1.00 .75 2 -303.17 -8.66 29.51 -5.73 -306.10 84.33 3.00 .75 3 -214.19 -6.33 51.28 5.63 -226.15 76.87 5.00 .75 4 -76.58 -16.76 69.22 28.73 -122.08 56.68 7.00 .75 5 348.24 -1.23 9.78 348.52 -1.50 1.60 1.00 .25 6 303.17 -4.26 30.59 306.18 -7.27 5.63 3.00 .25 7 214.19 8.88 48.82 225.21 -2.14 12.72 5.00 .25

- 22 -

有限元分析程序与算例

8 76.58 -43.42 70.89 109.45 -76.29 24.88 7.00 .25

NODES & PRINCIPAL STRESS

NODE NO. SIG-XX SIG-YY SIG-XY SIG-1 SIG-2 ANGLE

1 -348.242 -8.255 10.233 -7.948 -348.550 88.278 2 -325.706 -8.459 19.869 -7.219 -326.946 86.430 3 -258.680 -7.497 40.392 -1.161 -265.016 81.086 4 -145.387 5 -76.584 6 .000 7 .000 8 -.001 9 -.001 10 .000 11 348.243 12 325.706 13 258.678 14 145.386 15 76.585

-11.548 -16.764 -4.742 -5.601 -2.593 -14.408 -30.090 -1.228 -2.744 2.311 -17.268 -43.417 60.249 69.221 10.004 20.026 40.049 60.051 70.053 9.775 20.183 39.707 59.854 70.886 - 23 -

11.578 28.732 7.910 17.421 38.773 53.278 56.606 348.516 326.942 264.687 165.037 109.454 -168.513 -122.080 -12.652 -23.022 -41.367 -67.686 -86.696 -1.501 -3.979 -3.698 -36.919 -76.286 69.001 56.685 38.334 41.020 44.073 41.580 38.939 1.601 3.503 8.606 18.176 24.877

有限元分析程序与算例

平面刚架静力分析程序

一、 源程序：

! PF.F90(A program form analysis of plane frmae) ! Version 4.3 2007 ! Main program reads the control date & organizes the ! whole calculation by calling subroutines

IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION W(20000) CHARACTER IDFN*12,TITLE(5)*72 READ (*,'(A12)')IDFN OPEN (3,FILE=IDFN,STATUS='OLD') READ (3,'(A72)')(TITLE(M),M=1,5) READ (3,*)E,NM,NJ,NS,NLC L1=1 L2=L1+NM L3=L2+NM L4=L3+NM L11=L4+NM L12=L11+NJ L21=L12+NJ L22=L21+NS L31=L22+NS L41=L31+6*NM CALL IOMJS (TITLE,E,NM,NJ,NS,NLC,W(L1),W(L2),W(L3),W(L4),& W(L11),W(L12),W(L21),W(L22))

CALL LCVCT (NM,W(L1),W(L2),W(L31),NJ,N) CALL LCDIA (NM,N,W(L31),W(L41),W(L41),W(L41),MAXBDW,NA) L51=L41+N L52=L51+36 L53=L52+NA*2 L54=L53 L61=L54+N*2 NW=L61+6*NM-1 WRITE (*,1)NA,NW

1 FORMAT(/40X,'( NA=',I6,' )' /40X,'( NW=',I6,' )') CALL FORMA (E,NM,NJ,N,NA,W(L1),W(L2),W(L3),W(L4),W(L11),W(L12),W(L31),W(L51),W(L41),W(L52)) CALL AS (NS,N,NA,W(L21),W(L41),W(L52)) CALL LDLT (N,NA,W(L41),W(L52),W(L53)) DO 100 LC=1,NLC

- 1 -

有限元分析程序与算例

READ (3,*)NLJ L62=L61+NLJ L63=L62+NLJ L64=L63+NLJ L71=L61 L81=L71+6*NM CALL B0 (LC,N,NLJ,W(L54)) IF (NLJ.NE.0) CALL IOLJB (N,NLJ,W(L61),W(L62),W(L63),W(L64),W(L54)) READ (3,*)NLM L82=L81+NLM L83=L82+NLM L84=L83+NLM CALL F0 (NLM,NM,W(L71)) IF(NLM.NE.0) CALL IOLMFB (NM,NJ,N,NLM,W(L81),W(L82),W(L83),W(L84),& W(L1),W(L2),W(L11),W(L12),W(L31),W(L71),W(L54)) CALL BS (NS,N,W(L21),W(L22),W(L54)) CALL SLVEQ (N,NA,MAXBDW,W(L41),W(L52),W(L54)) CALL OJD (NJ,N,W(L54)) CALL COTF (E,NM,NJ,N,W(L1),W(L2),W(L3),W(L4),W(L11),W(L12),W(L31),W(L54),W(L71)) NW=L84+NLM-1 WRITE (*,1)NA,NW 100 CONTINUE WRITE (*,'(/)') END

! Read data of members,joints,supports & print them SUBROUTINE IOMJS (TITLE,E,NM,NJ,NS,NLC,IST,IEN,AR,RI,X,Y,IS,VS) IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION IST(NM),IEN(NM),AR(NM),RI(NM),X(NJ),Y(NJ),IS(NS),VS(NS) CHARACTER TITLE(5)*72 WRITE (*,'(/)') WRITE (*,1)(TITLE(M),M=1,5) 1 FORMAT(1X,A72) WRITE (*,2)E,NM,NJ,NS,NLC

2 FORMAT(/13X,'The Input Data'//10X,'The General Information'& //6X,'E',9X,'NM',5X,'NJ',5X,'NS',5X,'NLC' /1X,1PE10.3,4I7) READ (3,*)(IST(M),IEN(M),AR(M),RI(M),M=1,NM) WRITE (*,3)

3 FORMAT(/10X,'The Information of Members'//1X,'member',2X,'start',& 2X,'end',9X,'A',15X,'I') WRITE (*,4)(M,IST(M),IEN(M),AR(M),RI(M),M=1,NM) 4 FORMAT(1X,I4,I8,I6,1P2E16.6) READ (3,*)(X(M),Y(M),M=1,NJ)

- 2 -

有限元分析程序与算例

WRITE (*,5)

5 FORMAT(/10X,'The Joint Coordinates'//1X,'joint',11X,'X',17X,'Y') WRITE (*,6)(M,X(M),Y(M),M=1,NJ) 6 FORMAT(1X,I4,2F18.6) READ (3,*)(IS(M),VS(M),M=1,NS) WRITE (*,7)

7 FORMAT(/10X,'The Information of Supports'//4X,'IS',9X,'VS') WRITE (*,8)(IS(M),VS(M),M=1,NS) 8 FORMAT(1X,I5,F16.6) RETURN END

! Determin location vector of element SUBROUTINE LCVCT (NM,IST,IEN,LV,NJ,N) IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION IST(NM),IEN(NM),LV(6,NM) DO 100 M=1,NM I=IST(M)*3 J=IEN(M)*3 LV(1,M)=I-2 LV(2,M)=I-1 LV(3,M)=I LV(4,M)=J-2 LV(5,M)=J-1 LV(6,M)=J 100 CONTINUE N=NJ*3 RETURN END

! Determine location of diagonal elements of global stiffness SUBROUTINE LCDIA (NM,N,LV,MIN,IBDW,LD,MAXBDW,NA) IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION LV(6,NM),MIN(N),IBDW(N),LD(N) DO 100 I=1,N MIN(I)=I 100 CONTINUE DO 400 M=1,NM MINLV=LV(1,M) DO 200 I=2,6 IF (LV(I,M).LT.MINLV) MINLV=LV(I,M) 200 CONTINUE DO 300 I=1,6 IF (MINLV.LT.MIN(LV(I,M))) MIN(LV(I,M))=MINLV

- 3 -

有限元分析程序与算例

300 CONTINUE 400 CONTINUE MAXBDW=0 DO 500 I=1,N IBDW(I)=I-MIN(I)+1 IF (IBDW(I).GT.MAXBDW) MAXBDW=IBDW(I) 500 CONTINUE LD(1)=IBDW(1) DO 600 I=2,N LD(I)=LD(I-1)+IBDW(I) 600 CONTINUE NA=LD(N) RETURN END

! Calculate length, cosine & sine of member SUBROUTINE RLCS (M,NM,NJ,IST,IEN,X,Y,RL,C,S) IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION IST(NM),IEN(NM),X(NJ),Y(NJ) I=IST(M) J=IEN(M) X1=X(J)-X(I) Y1=Y(J)-Y(I) RL=SQRT(X1*X1+Y1*Y1) C=X1/RL S=Y1/RL RETURN END

! Calculate element stiffness matrix along local axes SUBROUTINE KEBAR (M,E,NM,NJ,IST,IEN,AR,RI,X,Y,C,S,E1,E2,E3,E4) IMPLICIT DOUBLE PRECISION (A-H,O-Z)

DIMENSION IST(NM),IEN(NM),X(NJ),Y(NJ),AR(NM),RI(NM) CALL RLCS (M,NM,NJ,IST,IEN,X,Y,RL,C,S) E1=E*AR(M)/RL E2=12.0*E*RI(M)/(RL*RL*RL) E3=0.5*E2*RL E4=0.6666667*E3*RL RETURN END

! Calculate element stiffness matrix along global axes SUBROUTINE KE (M,E,NM,NJ,IST,IEN,AR,RI,X,Y,AE) IMPLICIT DOUBLE PRECISION (A-H,O-Z)

- 4 -

有限元分析程序与算例

DIMENSION IST(NM),IEN(NM),AR(NM),RI(NM),X(NJ),Y(NJ),AE(6,6) CALL KEBAR (M,E,NM,NJ,IST,IEN,AR,RI,X,Y,C,S,E1,E2,E3,E4) A1=E1*C*C+E2*S*S A2=(E1-E2)*C*S A3=E1*S*S+E2*C*C A4=E3*S A5=E3*C A6=E4 AE(1,1)=A1 AE(2,1)=A2 AE(2,2)=A3 AE(3,1)=-A4 AE(3,2)=A5 AE(3,3)=A6 AE(4,1)=-A1 AE(4,2)=-A2 AE(4,3)=A4 AE(4,4)=A1 AE(5,1)=-A2 AE(5,2)=-A3 AE(5,3)=-A5 AE(5,4)=A2 AE(5,5)=A3 AE(6,1)=-A4 AE(6,2)=A5 AE(6,3)=0.5*A6 AE(6,4)=A4 AE(6,5)=-A5 AE(6,6)=A6 RETURN END

! Form global stiffness matrix A SUBROUTINE FORMA (E,NM,NJ,N,NA,IST,IEN,AR,RI,X,Y,LV,AE,LD,A) IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION

IST(NM),IEN(NM),AR(NM),RI(NM),X(NJ),Y(NJ),LV(6,NM),AE(6,6),LD(N) DOUBLE PRECISION A(NA) DO 300 M=1,NM CALL KE (M,E,NM,NJ,IST,IEN,AR,RI,X,Y,AE) DO 200 I=1,6 DO 100 J=1,I IF (LV(I,M).GE.LV(J,M)) THEN A(LD(LV(I,M))-LV(I,M)+LV(J,M))=A(LD(LV(I,M))-LV(I,M)+LV(J,M))+AE(I,J)

- 5 -

有限元分析程序与算例

ELSE A(LD(LV(J,M))-LV(J,M)+LV(I,M))=A(LD(LV(J,M))-LV(J,M)+LV(I,M))+AE(I,J) ENDIF

100 CONTINUE 200 CONTINUE 300 CONTINUE RETURN END

! Introduce support conditions into global stiffness matrix A SUBROUTINE AS(NS,N,NA,IS,LD,A)

IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION IS(NS),LD(N),A(NA) DO 100 M=1,NS I=3*(IS(M)/10)-3+MOD(IS(M),10) A(LD(I))=1D22 100 CONTINUE RETURN END

! Solve equations(1)-decomposition of matrix A SUBROUTINE LDLT(N,NA,LD,A,T)

IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION LD(N),A(NA),T(N) DOUBLE PRECISION SUM DO 300 I=2,N LDI=LD(I) I1=I-LDI+LD(I-1)+1 DO 200 J=I1,I-1 LDJ=LD(J) J1=J-LDJ+LD(J-1)+1 IF(I1.GT.J1) J1=I1 SUM=0.0D0 DO 100 K=J1,J-1 SUM=SUM+T(K)*A(LDJ-J+K) 100 CONTINUE T(J)=A(LDI-I+J)-SUM A(LDI-I+J)=T(J)/A(LDJ) A(LDI)=A(LDI)-T(J)*A(LDI-I+J) 200 CONTINUE 300 CONTINUE RETURN END

! Solve equations(2)-forward & back substitution SUBROUTINE SLVEQ (N,NA,MAXBDW,LD,A,B)

- 6 -

有限元分析程序与算例

IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION LD(N),A(NA),B(N) DO 200 I=2,N LDI=LD(I) I1=I-LDI+LD(I-1)+1 DO 100 J=I1,I-1 B(I)=B(I)-A(LDI-I+J)*B(J) 100 CONTINUE 200 CONTINUE DO 300 I=1,N B(I)=B(I)/A(LD(I)) 300 CONTINUE DO 500 I=N-1,1,-1 IMIN=I+MAXBDW IF(IMIN.GT.N) IMIN=N DO 400 J=I+1,IMIN LDJ=LD(J) J1=J-LDJ+LD(J-1)+1 IF(I.GE.J1) B(I)=B(I)-A(LDJ-J+I)*B(J) 400 CONTINUE 500 CONTINUE RETURN END

! Initialize joint load vector B SUBROUTINE B0 (LC,N,NLJ,B)

IMPLICIT DOUBLE PRECISION (A-H,O-Z) DOUBLE PRECISION B(N) WRITE (*,1)LC,NLJ

1 FORMAT(/15X,'Loading Case',I3 //10X,'The Loadings at Joints'//17X,'NLJ=',I4) DO 100 I=1,N B(I)=0.0D0 100 CONTINUE RETURN END

! Read data of loading at joint & form joint load vector B SUBROUTINE IOLJB (N,NLJ,ILJ,PX,PY,PM,B) IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION ILJ(NLJ),PX(NLJ),PY(NLJ),PM(NLJ),B(N) READ (3,*)(ILJ(M),PX(M),PY(M),PM(M),M=1,NLJ) WRITE (*,1)

1 FORMAT(/2X,'ILJ',11X,'PX',14X,'PY',15X,'PM') WRITE (*,2)(ILJ(M),PX(M),PY(M),PM(M),M=1,NLJ) 2 FORMAT(1X,I4,2F16.4,F18.5) DO 100 M=1,NLJ

- 7 -

有限元分析程序与算例

I=ILJ(M)*3 B(I-2)=PX(M) B(I-1)=PY(M) B(I)=PM(M) 100 CONTINUE RETURN END

! Initialize terminal forces of members SUBROUTINE F0 (NLM,NM,F)

IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION F(6,NM) WRITE (*,1) NLM

1 FORMAT(/10X,'The Loadings at Members'//17X,'NLM=',I4) DO 200 J=1,NM DO 100 I=1,6 F(I,J)=0.0 100 CONTINUE 200 CONTINUE RETURN END

! Read data of loading at member & calculate fixed-end forces SUBROUTINE IOLMFB (NM,NJ,N,NLM,ILM,ITL,PV,DST,IST,IEN,X,Y,LV,F,B) IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION

ILM(NLM),ITL(NLM),PV(NLM),DST(NLM),IST(NM),IEN(NM),X(NJ),& Y(NJ),LV(6,NM),F(6,NM),B(N) READ (3,*)(ILM(M),ITL(M),PV(M),DST(M),M=1,NLM) WRITE (*,1)

1 FORMAT(/2X,'ILM',2X,'ITL',11X,'PV',12X,'DST') WRITE (*,2)(ILM(M),ITL(M),PV(M),DST(M),M=1,NLM) 2 FORMAT(1X,I4,I5,F16.4,F16.6) DO 100 M=1,NLM L=ILM(M) CALL RLCS (L,NM,NJ,IST,IEN,X,Y,RL,C,S) D1=DST(M) D2=RL-D1 IF (ITL(M).EQ.1.OR.ITL(M).EQ.3) THEN P1=PV(M)*C P2=-PV(M)*S END IF IF (ITL(M).EQ.2.OR.ITL(M).EQ.4) THEN P1=PV(M)*S P2=PV(M)*C END IF

- 8 -

有限元分析程序与算例

IF (ITL(M).EQ.1.OR.ITL(M).EQ.2) THEN F1=-P1*D2/RL F4=-P1-F1 F2=-P2*D2*D2*(RL+2.0*D1)/(RL*RL*RL) F5=-P2-F2 F3=-P2*D1*D2*D2/(RL*RL) F6=P2*D1*D1*D2/(RL*RL) END IF IF (ITL(M).EQ.3.OR.ITL(M).EQ.4) THEN G=P2*D1*D1/(12.0*RL*RL) F3=-G*((6.0*RL-8.0*D1)*RL+3.0*D1*D1) F6=G*D1*(4.0*RL-3.0*D1) F5=-6.0*G*D1*(2.0-D1/RL) F2=-P2*D1-F5 F4=-P1*D1*D1/(2.0*RL) F1=-P1*D1-F4 END IF F(1,L)=F(1,L)+F1 F(2,L)=F(2,L)+F2 F(3,L)=F(3,L)+F3 F(4,L)=F(4,L)+F4 F(5,L)=F(5,L)+F5 F(6,L)=F(6,L)+F6 B(LV(1,L))=B(LV(1,L))-F1*C+F2*S B(LV(2,L))=B(LV(2,L))-F1*S-F2*C B(LV(3,L))=B(LV(3,L))-F3 B(LV(4,L))=B(LV(4,L))-F4*C+F5*S B(LV(5,L))=B(LV(5,L))-F4*S-F5*C B(LV(6,L))=B(LV(6,L))-F6 100 CONTINUE RETURN END

! Introduce support conditions into load vector B SUBROUTINE BS (NS,N,IS,VS,B)

IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION IS(NS),VS(NS),B(N) DO 100 M=1,NS I=3*(IS(M)/10)-3+MOD(IS(M),10) B(I)=VS(M)*1D22 100 CONTINUE RETURN END

! print joint displacements SUBROUTINE OJD (NJ,N,B)

- 9 -

有限元分析程序与算例

IMPLICIT DOUBLE PRECISION (A-H,O-Z) DOUBLE PRECISION B(N) WRITE (*,1)

1 FORMAT(/13X,'The Results of Calculation' //10X,'The Joint Displacements'& //1X,'joint',8X,'u',18X,'v',14X,'phi') WRITE (*,2)(M,B(3*M-2),B(3*M-1),B(3*M),M=1,NJ) 2 FORMAT(1X,I4,1P3E18.8) RETURN END

! Calculate & print terminal forces of members SUBROUTINE COTF (E,NM,NJ,N,IST,IEN,AR,RI,X,Y,LV,B,F) IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION

IST(NM),IEN(NM),AR(NM),RI(NM),X(NJ),Y(NJ),LV(6,NM),F(6,NM),B(N) DOUBLE PRECISION U1,U2,U3,U4,U5,U6 WRITE (*,1)

1 FORMAT(/10X,'The Terminal Forces'//1X,'member',4X,'N(st)',6X,'Q(st)',7X,& 'M(st)', 6X,'N(en)',6X,'Q(en)',7X,'M(en)') DO 100 M=1,NM CALL KEBAR (M,E,NM,NJ,IST,IEN,AR,RI,X,Y,C,S,E1,E2,E3,E4) U1=B(LV(1,M))*C+B(LV(2,M))*S U2=-B(LV(1,M))*S+B(LV(2,M))*C U3=B(LV(3,M)) U4=B(LV(4,M))*C+B(LV(5,M))*S U5=-B(LV(4,M))*S+B(LV(5,M))*C U6=B(LV(6,M)) F(1,M)=F(1,M)+E1*(U1-U4) F(2,M)=F(2,M)+E2*(U2-U5)+E3*(U3+U6) F(3,M)=F(3,M)+E3*(U2-U5)+E4*(U3+0.5*U6) F(4,M)=F(4,M)+E1*(U4-U1) F(5,M)=F(5,M)+E2*(U5-U2)-E3*(U3+U6) F(6,M)=F(6,M)+E3*(U2-U5)+E4*(0.5*U3+U6) WRITE (*,2)M,F(1,M),F(2,M),F(3,M),F(4,M),F(5,M),F(6,M) 2 FORMAT(1X,I4,2(2F12.3,F12.3)) 100 CONTINUE RETURN END

二、数据输入：

********************************** * * * PLANE FRAME EX.1 2007.6.14 * * *

- 10 -

有限元分析程序与算例

********************************** 2.1E8 2 3 3 1

1 2 102E-4 32240E-8 2 3 102E-4 32240E-8 0 0 4 3 8 3 11 0 12 -0.1 32 0 1

2 0 -20 0 2

2 2 -40 2 2 4 -10 4

三、结果输出：

********************************** * * * PLANE FRAME EX.1 2007.6.14 * * * **********************************

The Input Data

The General Information

E NM NJ NS NLC 2.100E+08 2 3 3 1

The Information of Members

member start end A I

1 1 2 1.020000E-02 3.224000E-04 2 2 3 1.020000E-02 3.224000E-04

The Joint Coordinates

joint X Y

1 .000000 .000000 2 4.000000 3.000000 3 8.000000 3.000000

- 11 -

有限元分析程序与算例

The Information of Supports

IS VS

11 .000000 12 -.100000 32 .000000

( NA= 36 ) ( NW= 179 ) Loading Case 1

The Loadings at Joints

NLJ= 1

ILJ PX PY PM 2 .0000 -20.0000 .00000

The Loadings at Members

NLM= 2

ILM ITL PV DST 2 2 -40.0000 2.000000 2 4 -10.0000 4.000000

The Results of Calculation

The Joint Displacements

joint u v phi

1 1.14352090E-34 -1.00000000E-01 7.87516325E-03 2 -2.80901559E-02 -6.26164867E-02 1.23062152E-02 3 -2.80901559E-02 -3.00000175E-21 1.78204134E-02

The Terminal Forces

member N(st) Q(st) M(st) N(en) Q(en) 1 18.000 24.000 .000 -18.000 -24.000 2 .000 10.000 -120.000 .000 70.000

( NA= 36 ) ( NW= 187 )

- 12 -

M(en) 120.000 .000

本文来源：https://www.bwwdw.com/article/y5o5.html