结构有限元分析程序研究生阶段.docx

1、结构有限元分析程序研究生阶段中 南 大 学CENTRAL SOUTH UNIVERSITY题 目 有限单元法程序设计 学生姓名 陈 煜 学 号 114811114 指导老师 周文伟 学 院 土木工程学院 专 业 岩土工程 完成时间 2011年12月 一、PF程序1.源程序C 主程序 读入并输出控制数据，分配动态数组地址，并通过调用各C 子程序来调用整个计算C PF.FOR (A program for analysis of plane frame)C Version 4.3 1994C Main program reads the control data & organizes the w

2、holeC calculation by calling subroutines. DIMENSION W(20000) CHARACTER IDFN*20,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 IO

3、MJS (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=

4、,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 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

5、) 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,

6、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,NW100 CONTINUE WRITE (*,(/) END SUBROUTINE IOMJS (TITLE,E,NM,NJ,NS,NLC,IST,IEN, & AR,RI,X,Y,IS,VS)C Read data of members, joints, support

7、s & print them 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

8、) 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) WRITE (*,5) 5 FORMAT(/10X,The Joint Coordinates & /

9、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 SUBROUTINE LCVCT (NM,IST,IEN,LV,NJ,N)C Determine location

10、 vector of element 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)=J100 CONTINUE N=NJ*3 RETURN END SUBROUTINE LCDIA (NM,N,LV,MIN,IBDW,LD,MAXBDW,NA)C Determine location of diagonal elements of global stiffnessC m

11、atrix A DIMENSION LV(6,NM),MIN(N),IBDW(N),LD(N) DO 100 I=1,N MIN(I)=I100 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)=MINLV300 CONTINUE400 CONTINUE MAXBDW=0 DO 500 I=1,N IBDW(I)=I-MIN(I)+1 IF (IBDW

12、(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 SUBROUTINE RLCS (M,NM,NJ,IST,IEN,X,Y,RL,C,S)C Calculate length, cosine & sine of member 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

13、(X1*X1+Y1*Y1) C=X1/RL S=Y1/RL RETURN END SUBROUTINE KEBAR (M,E,NM,NJ,IST,IEN,AR,RI, & X,Y,C,S,E1,E2,E3,E4)C Calculate element stiffness matrix along local axes 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*E

14、2*RL E4=0.6666667*E3*RL RETURN END SUBROUTINE KE (M,E,NM,NJ,IST,IEN,AR,RI,X,Y,AE)C Calculate element stiffness matrix along global axes 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

15、*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 SUBROUTINE FORMA (E

16、,NM,NJ,N,NA,IST,IEN,AR,RI, & X,Y,LV,AE,LD,A)C Form global stiffness matrix A 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)-

17、LV(I,M)+LV(J,M) & =A(LD(LV(I,M)-LV(I,M)+LV(J,M)+AE(I,J) 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) END IF100 CONTINUE200 CONTINUE300 CONTINUE RETURN END SUBROUTINE AS (NS,N,NA,IS,LD,A)C Introduce support conditions into global stiffness matrix A DIMENSION IS(NS),LD(N)

18、DOUBLE PRECISION A(NA) DO 100 M=1,NS I=3*(IS(M)/10)-3+MOD(IS(M),10) A(LD(I)=1D22100 CONTINUE RETURN END SUBROUTINE LDLT (N,NA,LD,A,T)C Solve equations (1) - decomposition of matrix A DIMENSION LD(N) DOUBLE PRECISION A(NA),T(N),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) J

19、1=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 CONTINUE300 CONTINUE RETURN END SUBROUTINE SLVEQ (N,NA,MAXBDW,LD,A,B)C Solve equations (2) - forward & back substitution DIME

20、NSION LD(N) DOUBLE PRECISION 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 CONTINUE200 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.J

21、1) B(I)=B(I)-A(LDJ-J+I)*B(J)400 CONTINUE500 CONTINUE RETURN END SUBROUTINE B0 (LC,N,NLJ,B)C Initialize joint load vector B 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.0D0100 CONTINUE RETURN END SUBROUTINE IOL

22、JB (N,NLJ,ILJ,PX,PY,PM,B)C Read data of loading at joint & form joint load vector B DIMENSION ILJ(NLJ),PX(NLJ),PY(NLJ),PM(NLJ) DOUBLE PRECISION 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(

23、1X,I4,2F16.4,F18.5) DO 100 M=1,NLJ I=ILJ(M)*3 B(I-2)=PX(M) B(I-1)=PY(M) B(I)=PM(M)100 CONTINUE RETURN END SUBROUTINE F0 (NLM,NM,F)C Initialize terminal forces of members 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.0100 C

24、ONTINUE200 CONTINUE RETURN END SUBROUTINE IOLMFB (NM,NJ,N,NLM,ILM,ITL,PV,DST, & IST,IEN,X,Y,LV,F,B)C Read data of loading at member & calculate fixed-end forces,C add equivalent joint loads to vector B DIMENSION ILM(NLM),ITL(NLM),PV(NLM),DST(NLM),IST(NM), & IEN(NM),X(NJ),Y(NJ),LV(6,NM),F(6,NM) DOUBL

25、E PRECISION 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(

26、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 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

27、) 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

28、*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)-F6100 CONTINUE RETURN END SUBROUTINE BS (NS,N,IS,VS,B)C Introduce support conditions into joint load vector B DIMENSION IS(NS),VS(NS) DOUBLE PRECISION B(N) DO 100 M

29、=1,NS I=3*(IS(M)/10)-3+MOD(IS(M),10) B(I)=VS(M)*1D22100 CONTINUE RETURN END SUBROUTINE OJD (NJ,N,B)C Print joint displacements DOUBLE PRECISION B(N) WRITE (*,1) 1 FORMAT(/13X,The Results of Calculation & /10X,The Joint Displacements & /1X,joint,8X,u,15X,v,14X,phi) WRITE (*,2)(M,B(3*M-2),B(3*M-1),B(3*M),M=1,NJ) 2 FORMAT(1X,I4,1P3E16.6) RETURN END SUBROUTINE COTF (E,NM,NJ,N,IST,IEN,AR,RI,X,Y,LV,B,F)C Calculate & print terminal forces of members DIMENSION IST(NM),IEN(NM),AR(NM),RI(NM),X(NJ),Y(NJ), & LV(6,NM),F(6,NM) DOUBLE PRECISION B(N),U1,U2,U3,U4,U5,U6 WRITE (*,1) 1 FORMAT(/10X