有限單元法FORTRAN程序設計作業(yè)

1、有限元分析是從二十世紀發(fā)展起來的一種新的結構分析方法，其分析問題的主要思想是把連續(xù)彈性體離散成為一群僅在節(jié)點處相互連接的有限單元的集合，通過把無限連續(xù)問題有限化進行進一步的分析研究。有限元最早是美國波音公司工程師特納采用三角形和矩形單元，把結構力學中的位移法應用于飛機結構的分析中；后來，人們逐漸認識到，有限元是一種基于虛功原理的廣義里茨法；到二十世紀七十年代后，隨著計算機技術的迅猛發(fā)展，人們采用各種成熟的單元編制了一系列結構分析程序，例如威爾遜教授等編制的SAP系列、拜塞等編制的ADINA系列、ANSYS公司研制的ANSYS系列等。通過本學期課程的學習，我們進一步學習了Fortran90編制空

2、間桁架的基本方法和理念，并對以下程序進行了調(diào)試，結合實例進行了計算。程序1：平面三節(jié)點有限元分析程序（）試用程序計算圖2所示懸臂梁結構，其受均布荷載q=1N/mm2 作用，×105 N/mm2 ，厚度h=10mm，有限元法分析其位移及應力。圖2 懸臂梁單元劃分及荷載簡圖梁視為平面應力狀態(tài)，按圖2尺寸劃分為均勻的三角形網(wǎng)格，共80個單元，55個節(jié)點，坐標軸及單元與節(jié)點的編號如圖。將均布荷載分配到相應的節(jié)點上，把有約束的節(jié)點51、52、53、54、55視作固定鉸鏈，建立如圖2所示的離散化計算模型。源程序：PROGRAM MAINDIMENSIONSK(300,30),EK(12,12),

3、Q(300),MC(55),XY(3,100),XYE(3,4),& QE(12),NX(4,100)OPEN(7,FILE='INPUT.TXT')REWIND 7READ (7,*) NF,NE,NN,MB,ND,LE,LSREAD (7,*) E,UM,T10 FORMAT(7I5)12 FORMAT(3F15.2)WRITE(*,600) NF,NE,NN,MB,ND,LE,LS,E,UM,TME=NE*NFMS=NN*NFCALL INPUT (XY,Q,NX,MC,LS,NN,MS,NE,LE,ND)WRITE (*,102) (XY(I,J),I=1,LS

4、),J=1,NN)102 FORMAT (10X,'XY'/,(2X,6F12.3)WRITE (*,101) (Q(I),I=1,MS)101 FORMAT (10X,'Q'/,(2X,6F12.3)WRITE (*,500) (NX(I,J),I=1,NE),J=1,LE)500 FORMAT (10X,'NX'/,(2X,12I6) WRITE (*,400) (MC(I),I=1,ND)600 FORMAT (10X,'NF NE NN MB ND LE LS E UM T'/7(2X,I4),3(2X,F8.4)400

5、FORMAT (10X,'MC'/,(2X,10I6)CALL STIFS (SK,EK,Q,NX,XY,XYE,MC,MS,MB,ME,ND,LE,NE,NF,& NN,LS,E,UM,T)CALL SOLVE (SK,Q,MS,MB)OPEN (9,FILE='OUT.DAT')REWIND 9WRITE (9,200)WRITE (9,250) (Q(I),I=1,MS)200 FORMAT (5X,'DISPLACEMENT')250 FORMAT (2X, 6E14.5)CALL STRES (Q,QE,NX,XY,XYE,MS

6、,ME,NE,LE,NF,NN,LS,E,UM,T)STOP 1000ENDSUBROUTINE INPUT (XY,Q,NX,MC,LS,NN,MS,NE,LE,ND)DIMENSION XY(LS,NN),Q(MS),NX(NE,LE),MC(ND)READ (7,*) XYREAD (7,*) QREAD (7,*) NXREAD (7,*) MCCLOSE(7)10 FORMAT(6F11.2)20FORMAT(12I5)RETURNENDSUBROUTINE STIFS(SK,EK,Q,NX,XY,XYE,MC,MS,MB,ME,ND,LE,NE,NF,& NN,LS,E,U

7、M,T)DIMENSION SK(MS,MB),EK(ME,ME),Q(MS),NX(NE,LE),MC(ND),& XY(LS,NN),XYE(LS,NE)DO 35 I=1,MSDO 35 J=1,MB35 SK (I,J)=0DO 200 L=1,LEDO 40 J=1,NELJ=NX(J,L)DO 40 I=1,LS40 XYE(I,J)=XY(I,LJ)DO 50 I=1,MEDO 50 J=1,ME50 CALL STIFE(EK,XYE,ME,NE,NF,LS,E,UM,T)IF(L.EQ.1) WRITE(*,70)EK70 FORMAT (10X,'EK

8、9;/,(6E14.5)DO 200 I=1,NEDO 200 II=1,NFM=NF*(I-1)+IIM1=NF*(NX(I,L)-1)+IIDO 200 J =1,NEDO 200 JJ=1,NF N=NF*(J-1)+JJN1=NF*(NX(J,L)-1)+JJMN=N1-M1+1IF(MN) 200,200,150150 SK(M1,MN)=SK(M1,MN)+EK(M,N)200 CONTINUEDO 220 I=1,NDM=MC(I)220 SK(M,1)=SK(M,1)*1E8RETURNENDSUBROUTINE SOLVE(SK,Q,MS,MB)DIMENSION SK(MS

9、,MB),Q(MS)K1=MS-1DO 125 K=1,K1IF(K+MB-1-MS) 105,106,106105 N=K+MB-1GOTO 110106 N=MS110 I1=K+1DO 125 I=I1,NL=I-K+1C=SK(K,L)/SK(K,1)J1=MB-L+1DO 122 J=1,J1M=J+I-K122 SK(I,J)=SK(I,J)-C*SK(K,M)125 Q(I)=Q(I)-C*Q(K)Q(MS)=Q(MS)/SK(MS,1)M=MS-1DO 145 I1=1,MI=MS-I1IF(MS-I+1-MB) 135,136,136135 N=MS-I+1GOTO 1401

10、36 N=MB140 DO 142 J=2,NL=J+I-1142 Q(I)=Q(I)-SK(I,J)*Q(L)145 Q(I)=Q(I)/SK(I,1)WRITE(*,147)147 FORMAT (5X, 'DISPLACEMENT' )WRITE( *, 150) (Q(I) ,I=1 ,MS)150 FORMAT(2X, 6E14.5)RETURNENDSUBROUTINE STRES (Q,QE,NX,XY,XYE,MS,ME,NE,LE,NF,NN,LS,E,UM,T)DIMENSION Q(MS),QE(ME),NX(NE,LE),XY(LS,NN),XYE(LS

11、,NE)DO 400 L=1,LEDO 160 I=1,NEDO 160 J=1,NFN=NF*(I-1)+JN1=NF*(NX(I,L)-1)+J160 QE(N)=Q(N1)WRITE (*,165) LWRITE (*,170) (QE(I),I=1,ME)165 FORMAT (4X,'L=',I4)170 FORMAT (6E14.5)DO 200 J=1,NELJ=NX(J,L)DO 200 I=1,LS200 XYE(I,J)=XY(I,LJ)CALL STE (XYE,QE,NE,LS,ME,E,UM,T)400 CONTINUERETURNENDSUBROUT

12、INE STIFE (EK,XYE,ME,NE,NF,LS,E,UM,T)DIMENSION EK(ME,ME),XYE(LS,NE),B(3),C(3)B(1)=XYE(2,2)-XYE(2,3)B(2)=XYE(2,3)-XYE(2,1)B(3)=XYE(2,1)-XYE(2,2)C(1)=XYE (1,3)-XYE (1,2)C(2)=XYE(1,1)-XYE(1,3)C(3)=XYE(1,2)-XYE(1,1)AE=(B(2)*C(3)-B(3)*C(2)/2D=E*T/4/(1-UM*2)/AEDO 30 I=1,3DO 30 J=1,3M=2*(I-1)+1N=2*(J-1)+1E

13、K(M,N)=D*(B(I)*B(J)+C(I)*C(J)*(1-UM)/2)EK(M,N+1)=D*(UM*B(I)*C(J)+C(I)*B(J)*(1-UM)/2)EK(M+1,N)=D*(UM*C(1)*B(J)+B(I)*C(J)*(1-UM)/2)30 EK(M+1,N+1)=D*(C(I)*C(J)+B(I)*B(J)*(1-UM)/2)RETURNENDSUBROUTINE STE (XYE,QE,NE,LS,ME,E,UM,T)DIMENSION XYE(LS,NE),QE(ME),SG(3),B(3),C(3),S(3,6)B(1)=XYE(2,2)-XYE(2,3)B(2)

14、=XYE(2,3)-XYE(2,1)B(3)=XYE(2,1)-XYE(2,2)C(1)=XYE(1,3)-XYE(1,2)C(2)=XYE(1,1)-XYE(1,3)C(3)=XYE(1,2)-XYE(1,1)AE=(B(2)*C(3)-B(3)*C(2)/2D=E/2/(1-UM*2)/AEDO 220 I=1,3S(1,2*I-1)=D*B(I)S(2,2*I-1)=D*UM*B(I)S(3,2*I-1)=D*(1-UM)*C(I)/2S(1,2*I)=D*UM*C(I)S(2,2*I)=D*C(I)220 S(3,2*I)=D*(1-UM)*B(I)/2DO 250 I=1,3SG(I

15、)=0DO 250 K=1,ME250 SG(I)=SG(I)+S(I,K)*QE(K)WRITE (9,300) SG300 FORMAT (5X,'SEGMA',3E20.4)RETURNEND輸入文本：輸出文本：程序2：平面剛架靜力分析程序（）試用程序計算圖1所示結構。梁和柱均用45a號工字鋼，截面面積A=102cm2，截面二次距I=32240cm4，彈性模量為×108kPa。圖1 平面鋼架荷載簡圖源程序：! PF.FOR (A program for the analysis of plane frame)! Version 6.3 2004! MAIN PR

16、OGRAM ! Reading and printing the control data,allocating the arrays ! and organizing the whole calculation by calling subroutinesDIMENSION W(80000)CHARACTER IDFN*12,OUTFN*12,TITLE(5)*72WRITE(*,1)1 FORMAT(1X,'Input Data File Name:')READ(*,'(A12)')IDFNOPEN(3,FILE=IDFN,STATUS='OLD&#

17、39;)WRITE(*,2)2 FORMAT(1X,'Output File Name:')READ(*,'(A12)')OUTFNOPEN(4,FILE=OUTFN,STATUS='NEW')WRITE(4,3)IDFN3 FORMAT(1X,'Input Data File Name:',1X,A12)WRITE(4,4)OUTFN4 FORMAT(5X,'Output File Name:',1X,A12)READ(3,'(A72)') (TITLE(M),M=1,5)WRITE(4,5) (

18、TITLE(M),M=1,5)5 FORMAT(1X,A72)WRITE(4,6)6 FORMAT(/13X,'The Input Data'/10X,'The General& Information'/6X,'E',9X,'NM',5X,'NJ',5X,'NS',5X,'NLC')READ(3,*)E,NM,NJ,NS,NLCWRITE(4,7)E,NM,NJ,NS,NLC7 FORMAT(1X,1PE10.3,4I7)WRITE(4,8)8 FORMAT(/10X,&#

19、39;The Information of Members' &9 /1X,'member',2X,'start',2X,'end',9X,'A',15X,'I')L1=1L2=L1+NML3=L2+NML4=L3+NML11=L4+NML12=L11+NJL21=L12+NJL22=L21+NSL31=L22+NSL32=L31+3*NJL41=L32+6*NMCALL IOMJS (NM,NJ,NS,NLC,W(L1),W(L2),W(L3),W(L4),&W(L11),W(L12),W

20、(L21),W(L22)CALL IDUM (NJ,N,W(L31),W(L11),W(L12)CALL LCVCT (NM,W(L1),W(L2),W(L32),NJ,W(L31)CALL LCDIA (NM,N,W(L32),W(L41),W(L41),W(L41),MAXBDW,NA)L51=L41+NL52=L51+72L53=L52+NA*2L54=L53L61=L54+N*2CALL FORMA (E,NM,NJ,N,NA,W(L1),W(L2),W(L3),W(L4),&W(L11),W(L12),W(L32),W(L51),W(L41),W(L52)CALL AS(NS

21、,NJ,N,NA,W(L21),W(L31),W(L41),W(L52)CALLLDLT(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*2 CALL B0(LC,N,NLJ,W(L54) IF(NLJ.NE.0) CALL IOLJB(N,NJ,NLJ,W(L31),W(L61),W(L62),&W(L63),W(L64),W(L54) READ(3,*)NLM L82=L81+NLM L83=L82+NLM

22、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(L32),W(L71),W(L54) CALLBS(NS,NJ,N,W(L21),W(L22),W(L31),W(L54) CALLSLVEQ(N,NA,MAXBDW,W(L41),W(L52),W(L54) CALL OJD(NJ,N,W(L31),W(L54) CALL COTF(E,NM,NJ,N,W(L1),W(L2),W(L

23、3),W(L4),W(L11),&W(L12),W(L32),W(L54),W(L71)100CONTINUEENDSUBROUTINE IOMJS(NM,NJ,NS,NLC,IST,IEN,AR,RI,X,Y,IS,VS)! Reading and printing the data of members, joints and supports DIMENSION IST(NM),IEN(NM),AR(NM),RI(NM),X(NJ),&Y(NJ),IS(NS),VS(NS)READ(3,*)(IST(M),IEN(M),AR(M),RI(M),M=1,NM)WRITE(4

24、,1)(M,IST(M),IEN(M),AR(M),RI(M),M=1,NM)1FORMAT(1X,I4,I8,I6,1P2E16.6)WRITE(4,2)2FORMAT (/10X,'The Joint Coordinates'/1X,'joint',8X,'X',15X,'Y')READ(3,*)(X(M),Y(M),M=1,NJ)WRITE(4,3)(M,X(M),Y(M),M=1,NJ)3FORMAT(1X,I4,2F16.6) WRITE(4,4)4FORMAT(/10X,'The Information of

25、Supports'/4X,'IS',8X,'VS')READ(3,*)(IS(M),VS(M),M=1,NS)WRITE(4,5)(IS(M),VS(M),M=1,NS)5FORMAT(1X,I5,F15.6)RETURNENDSUBROUTINE IDUN(NJ,N,IU,X,Y)! Determining the numbers ofindependent unknowns DIMENSION IU(3,NJ),X(NJ),Y(NJ)IU(1,1)=1IU(2,1)=2IU(3,1)=3N=3DO 200 J=2,NJ DO 100 I=J-1,1,

26、-1 DX=ABS(X(J)-X(I) DY=ABS(Y(J)-Y(I) IF(DX.LT.1E-4.AND.DY.LT.1E-4)THEN IU(1,J)=IU(1,I) IU(2,J)=IU(2,I) IU(3,J)=N+1 N=IU(3,J) GOTO 150 END IF100 CONTINUE IU(1,J)=N+1 IU(2,J)=N+2 IU(3,J)=N+3 N=IU(3,J)150 CONTINUE200 CONTINUERETURNENDSUBROUTINE LCVCT(NM,IST,IEN,LV,NJ,IU)! Determine the location vector

27、of elementDIMENSION IST(NM),IEN(NM),LV(6,NM),IU(3,NJ)DO 100 M=1,NM I=IST(M) J=IEN(M)LV(1,M)=IU(1,I)LV(2,M)=IU(2,I)LV(3,M)=IU(3,I)LV(4,M)=IU(1,J)LV(5,M)=IU(2,J) LV(6,M)=IU(3,J)100 CONTINUERETURNENDSUBROUTINE LCDIA(NM,N,LV,MIN,IBDW,LD,MAXBDW,NA)! Determine the location of diagonal elements of global s

28、tiffness! matrix ADIMENSIONLV(6,NM),MIN(N),IBDW(N),LD(N)DO 100 I=1,NMIN(I)=I100CONTINUEDO 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 CONTINUE400CONTINUEMAXBDW=0DO 500 I=1,N IBDW(I)=I-MIN(I)+1 IF(IBDW

29、(I).GT.MAXBDW) MAXBDW=IBDW(I)500CONTINUELD(1)=IBDW(1)DO 600 I=2,N LD(I)=LD(I-1)+IBDW(I)600CONTINUENA=LD(N)RETURNENDSUBROUTINE LCS(M,NM,NJ,IST,IEN,X,Y,RL,C,S)! Calculating length, cosine and sine of memberDIMENSION IST(NM),IEN(NM),X(NJ),Y(NJ)DOUBLE PRECISION RL,C,S,X1,Y1I=IST(M)J=IEN(M)X1=X(J)-X(I) Y

30、1=Y(J)-Y(I)RL=SQRT(X1*X1+Y1*Y1)C=X1/RL S=Y1/RLRETURNENDSUBROUTINEKEBAR (M,E,NM,NJ,IST,IEN,AR,RI,X,Y,C,S,E1,E2,E3,E4)! Calculating element stiffness matrix in local coordinatesDIMENSION IST(NM),IEN(NM),X(NJ),Y(NJ),AR(NM),RI(NM)DOUBLE PRECISION RL,C,S,E1,E2,E3,E4CALL LCS (M,NM,NJ,IST,IEN,X,Y,RL,C,S)E1

31、=E*AR(M)/RLE2=12.0*E*RI(M)/(RL*RL*RL)E3=0.5*E2*RLE4=4.0*E*RI(M)/RLRETURNENDSUBROUTINE KE (M,E,NM,NJ,IST,IEN,AR,RI,X,Y,AE)! Calculating element stiffness matrix in global coordinatesDIMENSION IST(NM),IEN(NM),AR(NM),RI(NM),X(NJ),Y(NJ)DOUBLE PRECISION AE(6,6),C,S,E1,E2,E3,E4,A1,A2,A3,A4,A5,A6CALL KEBAR

32、 (M,E,NM,NJ,IST,IEN,AR,RI,X,Y,C,S,E1,E2,E3,E4)A1=E1*C*C+E2*S*SA2=(E1-E2)*C*SA3=E1*S*S+E2*C*CA4=E3*SA5=E3*CA6=E4AE(1,1)=A1AE(2,1)=A2AE(2,2)=A3AE(3,1)=-A4AE(3,2)=A5AE(3,3)=A6AE(4,1)=-A1AE(4,2)=-A2AE(4,3)=A4AE(4,4)=A1AE(5,1)=-A2AE(5,2)=-A3AE(5,3)=-A5AE(5,4)=A2AE(5,5)=A3AE(6,1)=-A4AE(6,2)=A5AE(6,3)=0.5*

33、A6AE(6,4)=A4AE(6,5)=-A5AE(6,6)=A6RETURNENDSUBROUTINE FORMA (E,NM,NJ,N,NA,IST,IEN,AR,RI,X,Y,LV,AE,LD,A)! Forming theglobal stiffness matrix ADIMENSION IST(NM),IEN(NM),AR(NM),RI(NM),X(NJ),Y(NJ),&LV(6,NM),LD(N)DOUBLE PRECISION A(NA),AE(6,6)DO 300 M=1,NM CALL KE (M,E,NM,NJ,IST,IEN,AR,RI,X,Y,AE) DO 2

34、00 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)ELSEA(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 IF100CONTINUE200CONTINUE300CONTINUERETURNENDSUBROUTINE AS (NS,NJ,N,NA,IS,IU,LD,A)! Introducing support

35、 conditions into the global stiffness matrix ADIMENSION IS(NS),IU(3,NJ),LD(N)DOUBLE PRECISION A(NA)DO 100 M=1,NS J=IS(M)/10 I=MOD(IS(M),10) K=IU(I,J) A(LD(K)=1D22100CONTINUERETURNENDSUBROUTINE LDLT (N,NA,LD,A,T)! Solve equations (1) - decomposition of matrix ADIMENSION LD(N)DOUBLE PRECISION A(NA),T(

36、N),SUMDO 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 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 CONTINUE300CONTINUERETURNENDSUBROUTINE SLVEQ (N,NA,MAXBDW,LD

37、,A,B)! Solve equations (2) - forward and back substitutionDIMENSION 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 CONTINUE200CONTINUEDO 300 I=1,N B(I)=B(I)/A(LD(I)300CONTINUEDO 500 I=N-1,1,-1 IMIN=I+MAXBDW IF(IMIN.GT.N) IMIN=N

38、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 CONTINUE500CONTINUERETURNENDSUBROUTINE B0 (LC,N,NLJ,B)! Initializing joint load vector B to zeroDOUBLE PRECISION B(N) WRITE (4,1)LC,NLJ1FORMAT(/15X,'Loading Case',I3/10X,'The Loadings at Joints'/1

39、7X,'NLJ=',I4)DO 100 I=1,N100CONTINUE RETURNENDSUBROUTINE IOLJB (N,NJ,NLJ,IU,LJ,FX,FY,FM,B)! Reading and printing the data of loadings at ! joints and forming joint load vector BDIMENSIONIU(3,NJ),LJ(NLJ),FX(NLJ),FY(NLJ),FM(NLJ)DOUBLEPRECISION B(N)WRITE (4,1)1 FORMAT(/1X,'joint',8X,

40、9;FX',14X,'FY',16X,'FM') READ (3,*)(LJ(M),FX(M),FY(M),FM(M),M=1,NLJ) WRITE (4,2)(LJ(M),FX(M),FY(M),FM(M),M=1,NLJ)2 FORMAT(1X,I4,2F16.6,F18.6)DO 100 M=1,NLJ J=LJ(M) B(IU(1,J)=FX(M) B(IU(2,J)=FY(M) B(IU(3,J)=FM(M)100CONTINUERETURN ENDSUBROUTINE F0(NLM,NM,F)! Initializing the termin

41、al forces of members to zeroDOUBLE PRECISION F(6,NM)WRITE (4,1)NLM1FORMAT(/10X,'The Loadings at Members'/17X,'NLM=',I4)DO 200 J=1,NM DO 100 I=1,6100 CONTINUE200CONTINUE RETURNENDSUBROUTINE IOLMFB(NM,NJ,N,NLM,LM,LT,VF,DST,IST,IEN,X,Y,LV,F,B)! Reading and printing the data of loadings

42、at members,calculating fixed-end ! forces and equivalent joint loads,adding the latter to the joint load vector BDIMENSION LM(NLM),LT(NLM),VF(NLM),DST(NLM),IST(NM), &IEN(NM),X(NJ),Y(NJ),LV(6,NM)DOUBLE PRECISION RL,C,S,D1,D2,P1,P2,F1,F2,F3,F4,&F5,F6,G,B(N),F(6,NM)WRITE (4,1)1FORMAT(/1X,'m

43、ember',2X,'type',9X,'VF',14X,'DST')READ (3,*) (LM(M),LT(M),VF(M),DST(M),M=1,NLM)WRITE (4,2)(LM(M),LT(M),VF(M),DST(M),M=1,NLM)2FORMAT(1X,I4,I7,1X,2F16.6)DO 100 M=1,NLM L=LM(M) CALL LCS (L,NM,NJ,IST,IEN,X,Y,RL,C,S) D1=DST(M) D2=RL-D1 IF (LT(M).EQ.1.OR.LT(M).EQ.3)THEN P1=VF(

44、M)*C P2=-VF(M)*S ENDIF IF(LT(M).EQ.2.OR.LT(M).EQ.4)THEN P1=VF(M)*S P2=VF(M)*C ENDIF IF(LT(M).EQ.1.OR.LT(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) ENDIF IF(LT(M).EQ.3.OR.LT(M).EQ.4)THEN G=P2*D1*D1/(12.0*RL*RL) F3=-

45、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 ENDIF 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(L

46、V(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)-F6100CONTINUERETURNENDSUBROUTINE BS(NS,NJ,N,IS,VS,IU,B)! Introducing support conditions into the joint load vector BDIMENSION IS(NS),VS(NS),IU(3,NJ)DOUBLE PRECISION B(N)DO 100 M=1,NS J=IS(M)/10 I=MOD(

47、IS(M),10) K=IU(I,J) B(K)=VS(M)*1D22100 CONTINUERETURNENDSUBROUTINE OJD(NJ,N,IU,B)! Printing the joint displacementsDIMENSION IU(3,NJ)DOUBLE PRECISION B(N)WRITE (4,1)1 FORMAT(/13X,'The Results of Calculation'/10X,'The Joint Displacements'& /1X,'joint',8X,'u',15X,&#

48、39;v',12X,'rotation')WRITE (4,2)(J,B(IU(1,J),B(IU(2,J),B(IU(3,J),J=1,NJ)2 FORMAT(1X,I4,1P3E16.6)RETURNENDSUBROUTINE COTF(E,NM,NJ,N,IST,IEN,AR,RI,X,Y,LV,B,F)! Calculating and printing terminal forces of membersDIMENSION IST(NM),IEN(NM),AR(NM),RI(NM),X(NJ),Y(NJ),LV(6,NM)DOUBLE PRECISION C,

49、S,E1,E2,E3,E4,U1,U2,U3,U4,U5,U6,B(N),F(6,NM)WRITE (4,1)1 FORMAT(/10X,'The Terminal Forces'&2 /1X,'member',18X,'FN',14X,'FS',16X,'M')DO 100 M=1,NMCALLKEBAR(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)*SU2=-B(LV(1,M)*S+B(LV(2,M)*C U3=

50、B(LV(3,M)U4=B(LV(4,M)*C+B(LV(5,M)*SU5=-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(4,2) M,IST(M),F(1

51、,M),F(2,M),F(3,M),IEN(M),F(4,M),F(5,M),F(6,M)2 FORMAT(1X,I4,2X,'start',I4,2F16.6,F18.6,/9X,'end',I4,2F16.6,F18.6)100 CONTINUERETURNEND輸入文本：輸出文本：接下頁程序3：薄板彎曲分析程序（PLANE.FOR*）試用程序PLANE.FOR*計算圖3所示結構，取其1/4如圖4所示。圖3為兩對邊簡支、另外兩對邊固支的矩形薄板。其上受有均布壓力q的作用，已知：a=1000mm，b=500mm，板厚t=20mm，橫向荷載q=-1N/mm2，&

52、#215;108kPa，。求對稱截面上的撓度及彎矩Mx、My分布曲線。圖3矩形薄板四邊支撐示意圖圖41/4板單元編號圖源程序：PROGRAM MAINDIMENSION SK(90,21),EK(12,12),Q(90),MC(35),XY(2,30),XYE(2,4),QE(12),NX(4,30)OPEN(7,FILE='INP.DAT')REWIND 7READ(7,*) NF,NE,NN,MB,ND,LE,LSREAD(7,*) E,UM,T10 FORMAT(7I5)12 FORMAT(3F15.2)WRITE(*,600) NF,NE,NN,MB,ND,LE,LS,

53、E,UM,TME=NE*NFMS=NN*NFCALL INPUT(XY,Q,NX,MC,LS,NN,MS,NE,LE,ND)WRITE(*,102)(XY(I,J),I=1,LS),J=1,NN)102 FORMAT(10X,'XY'/,(2X,6F12.3)WRITE(*,101)(Q(I),I=1,MS)101 FORMAT(10X,'Q'/,(2X,6F12.3)WRITE(*,500)(NX(I,J),I=1,NE),J=1,LE)WRITE(*,400)(MC(I),I=1,ND)500 FORMAT(10X,'NX'/,(2X,12I

54、6)600 FORMAT(10X,'NF NE NN MB ND LE LS E UM T'/7(2X,I4),3(2X,F8.4)400 FORMAT(10X,'MC'/,(2X,10I6)CALL STIFS(SK,EK,Q,NX,XY,XYE,MC,MS,MB,ME,ND,LE,NE,NF,NN,& LS,E,UM,T)CALL SOLVE(SK,Q,MS,MB)OPEN(9,FILE='OUT.DAT')REWIND 9WRITE(9,200)WRITE(9,250)(Q(I),I=1,MS)200 FORMAT(5X,'

55、DISPLACEMENT')250 FORMAT(2X,6E14.5)CALL STRES(Q,QE,NX,XY,XYE,MS,ME,NE,LE,NF,NN,LS,E,UM,T)STOP 1000ENDSUBROUTINE INPUT(XY,Q,NX,MC,LS,NN,MS,NE,LE,ND)DIMENSION XY(LS,NN),Q(MS),NX(NE,LE),MC(ND)READ(7,*)XYREAD(7,*)QREAD(7,*)NXREAD(7,*)MCCLOSE(7)10 FORMAT(6F11.2)20 FORMAT(12I5)RETURNENDSUBROUTINE STIFS(SK,EK,Q,NX,XY,XYE,MC,MS,MB,ME