Fortran平面鋼架有限元分析

1、1有限元分析軟件的開發(fā)1.1 程序功能該程序為平面剛架靜力分析程序，能針對平面剛架問問題進行有限元計算，計算桿端位移及桿端力大小。程序從磁盤文件中讀取單元編號、節(jié)點編號及坐標、 材料屬性、荷載、邊界條件等信息；將桿端位移，桿端力等計算結(jié)果以磁盤文件 的形式輸出，采用等帶寬二維數(shù)組存儲整體剛度矩陣并使用高斯消去法進行求 解。1.2 程序結(jié)構(gòu)及流程1.3 程序的輸入與輸出詳細介紹輸入輸出數(shù)據(jù)的格式。如：數(shù)據(jù)文件分幾個部分，各有幾行，分別包含哪些容及其類型、先后次序，等等。輸入，共有九彳T。第一行：7,13,5,1,2,2。分另J為，7個結(jié)點，13個自由度，5個單元，1個類型，2個結(jié)點荷載，2個非結(jié)

2、點荷載。第二行：1,2,3,0.0,0.0,0,0,6.0,0.0分別為：一號結(jié)點的位移序號，x方向為1, y方向為2,轉(zhuǎn)角為3,坐標為（0.0,0.0）,因為二號結(jié)點固結(jié)在地面，所以二號結(jié)點的位移序號，x方向為0, y方向為0,轉(zhuǎn)角為0,坐標為（6.0,0.0）。第三行：4,5,6,0.0,6.0,4,5,7,0.0,6.0分別為：三號結(jié)點的位移序號，x方向為4, y方向為5,轉(zhuǎn)角為6,坐標為（0.0,6.0）,四號結(jié)點位移序號x方向和y相同， 轉(zhuǎn)角為7,坐標為（,0.0,6.0>第四行：8,9,10,6.0,6.0,0,0,11,0.0,12.0£號結(jié)點位移序號，x方向為8

3、, y方向為9,轉(zhuǎn)角為10,坐標為（6.0,6.0）。因為六號結(jié)點較接在地面，所以六號結(jié)點的位移序號，x方向和y方向為0,轉(zhuǎn)角為11,坐標為（0.0,12.0）。第五行：12,0,13,6.0,12.0.因為七號結(jié)點與地面用滑動支座固定，所以七號結(jié)點的位移序號，x方向為12, y方向為0,轉(zhuǎn)角為13,坐標（6.0,12.0）.第六行：1,2,1,1,3,1,4,5,1,3,6,1,5,7,1分別為，1號和2號結(jié)點組成的單元為1號類型。1號和3號結(jié)點組成的單元為1號類型，4號和5號結(jié)點組成的為1號 類型，3號和6號結(jié)點組成的單元為1號類型，5號和7號結(jié)點組成的單元為1號類型。第七行：分別為，彈性

4、模量為 E= 2X108 kN/m2 ,截面面積A=0.16m2,慣性矩 I=0.002m4 。第八行：1號結(jié)點轉(zhuǎn)角方向的集中力偶為-20.0kN ,3號結(jié)點集中力為10.0KN第九行：1號單元，受集中力（集中力型號為 3）,大小為15.0kN,到始端 的距離為3.0o 5號單元，受均布力（均布力型號為1）,大小為5.0kN,到端點 的距離為5.0。第十行：0為計算終止符。輸出：第一部分為輸入的數(shù)據(jù)。 RESULTSOF CALCULATION以下為輸出結(jié) 果，第二部分的第一段為4個結(jié)點的x,y方向的位移和轉(zhuǎn)角。第二段為1,2,3號 單元的軸力，剪力和彎矩。1.4程序求解中遇到的問題1對實例進

5、行計算時，坐標原點選用不同的點，會導(dǎo)致整個題目的坐標值發(fā) 生改變，輸入的容會有所不同，最后的結(jié)果也不相同2對結(jié)點荷載和非結(jié)點荷載的正負判斷不同，結(jié)點荷載的方向和整體坐標有 關(guān)，非結(jié)點荷載方向判斷和局部坐標有關(guān)。3在非結(jié)點荷載中，均布荷載和集中力到始端的距離判斷不同。2有限元分析算例2.1 算例說明已知圖示剛架，各桿的材料及截面均相同，彈性模量E= 2X108 kN/m2,A=0.16m2,慣性矩I=0.002m4 , q=5kN/m.，一號單元集中力為 15KN, 一號結(jié)點集中力偶為 20KN*M,三號結(jié)點集中力為10KN.試求剛架的力。節(jié)點編號如圖2.2 理論分析對所選取的力學(xué)問題進行理論分

6、析，要有詳細的推導(dǎo)過程和計算結(jié)果。1力計算1500因此對結(jié)構(gòu)進行分析，可以看出1,2,4單元組成的是二次超靜定結(jié)構(gòu)，3,5單元是靜定結(jié)構(gòu)。先對3,5單元組成的結(jié)構(gòu)進行分析。如上圖所示，可以根據(jù) x, y方向力平衡，對結(jié)點七力矩平衡算得支座反力。再畫出其彎矩, 剪力軸力圖。然后對 1,2,4單元組成的結(jié)構(gòu)分析。ft： x1.00用力法解超靜定，將結(jié)點六的約束解除，加上支座反力x1=1 , x2=1.畫出 M1 ,M2,MP 圖。MP圖M1圖M2圖肉 IXI+crl2X2 + AIP = 0 b21Xl + b22X2 + A2P = 0EI 2 EI1 (12 、1仃22 =*6*6*6 =72

7、 E八 23 JEI11?J440內(nèi)2 = o2l = (12*6 字 I2 + *12汩2*二"2)=EI23EIIP ,(3*100 町/+1*3*45 二 *3 + 3* 145*4.5)EI2 232542.5 EIA2P-(-*6*120*10)- -(3*100M2 + -*3*45* 12+3*145*12)El 2El213230EIXI，,0&犬2胃-14.06然后畫出其彎矩，剪力，軸力圖彎矩圖剪力圖軸力圖2位移計算計算結(jié)點1位移，x方向加單位力1根據(jù)公式:A丘乜區(qū)小十EA軸力圖圖乘A1x= (12.92*6*1) /EA= (12.92*6*1) /3*1

8、0A7*0.16=1.615e-5y方向加單位力1,忽略剪力的影響,彎矩圖圖乘Hy=-1/EI(3*15.05*3/2+1/2*3*2.82*2/3*3+1/2*1.27*(3+1.27/3)*17.87)+1/EI(1/2*1.73*24.32*(4.27+2*1.73/3)=1.611*E-5 (略小于程序結(jié)果)加單位力偶，剪力與軸力圖為零,彎矩圖為81=1/EI(3*15.05*1+1/2*3*2.82*1+1/2*1.27T17.87)- 1/EI(1/2*1.73*24.32*1)=0.6617e-3計算結(jié)點3的位移，x方向加單位力1.-&0C軸力彎矩1X=-1/EI(3*1

9、5.05*6+1/2*3*2.82*6+1/2*1.27*6*17.87)+1/EI(1/2*1.73*24.32*6)-1/EI(1/2*2.17*35.05*(3.29+2*2.71/3) )+1/EI(1/2*3.29*42.47*2.71/3)+(12.92*6*1) /EA =-0.675E-2.其它位移同理可得。2.3輸入輸出數(shù)據(jù)輸入:1平面剛架分析例題 著應(yīng) S, 1,2,2a6.J- J o o 6.o.J Jo 70012. 0LIL 2, 3, 0. 0, 0. 0, 0,0, 4,5, 6, 0, 0j 6. 0, 4, 5, a 9,10f 6. 0,6, 0, 0,

10、0 12, a 13, b 0.12. 0L £ L L 及 L 5 5, L 3. 0E7, 0.16, 0. 002 3, -20. 0. 4,10, 0 1, 3. L5. 3,。,5,1, 5 0輸出:WO= 1)平面剛架分析例題NJ二；7N=；13NE=;Mk; 1NPJ 二：2NPF 二;NO.E (1) (Z) m. KAT143611557113 067101113Y ,000 6. 000o o o o o o o o o o o o66. DOO,000.0006. 000$ 0006. 00012. 00012. 000闌.MAT ELASTIC MODULUS

11、 AREA1. 3000E+08.1600MOOT。INERTIA.2000E-02MODEL LOADSMO.DISF.VALUE3-20.000410.000HON-NOEEL LOADSMO.E NO. LOAD. MODEL AC1315.0003.000515.0005.000* * *料* *期 SULT S OF CALCJLATI0 M* * 址 * * 址* 科林* *MODEL DISPLACEHEITTNO.NX-DISPLACEMEITT¥-DISPLACEMENTANG. ROT. (RAD)116152E-04-1640&E-04.661 TIE-

12、032.00000.00000.000003-+ 67500E-02-17578E-04.29077E-03467500E-02-17E78E-04-14939E-02557875E-02t1875QB-Q4.30Q61E-Q26. 00000.0000018329E-02738324E-01.QOOOQ,60061E-02NO.E M(l)1 -12. 921222.9375970330. 000144 -14. 06241515.00000ELEMAMT FORCESQ(l) M。)M2)-537596815 0543012.92122-12. 92122-35.州始9- 53759701

13、5, 00000.1989305E-05-30.000147, 07883942.4730114,06241-30. 00001-90.00009-15. 00000QM-14. 0524024. 3201212. 92122-42. 47301-15. 0000090. 00001078839.2563000E-05.5722046E-05-. 7629395E-O52.4分析結(jié)果理論分析中，因為力計算應(yīng)用了力法，所以程序所得結(jié)果和理論結(jié)果一致 而對位移進行理論分析時忽略了剪力的影響，所以理論位移略小于程序所計算的結(jié)果?？梢钥闯鲕浖牡恼_性很高，但是此軟件只適用計算平面桿系結(jié)構(gòu)， 不 能解

14、決彈性力學(xué)問題.結(jié)點1位移uv9理論值1.615e-51.611*E-50.6617e-3程序值1.615e-51.640*E-50.6617e-33程序源代碼附上完整的程序源代碼。PROGRAM PFAPC ANALYSIS PROGRAM FOR PLANE FRAMEREAL K(200,200),KE(6,6),AKE(6,6),X(100),Y(100),AL(100),#EAI(3,100),PJ(100),PF(2,100),R(6,6),P(100),FF(6),#FE(6),D(100),ADE(6),DE(6),RT(6,6),AFE(6),F(3)INTEGER JE(2

15、,100),JN(3,100),JPJ(100),JPF(2,100),M(6),JEAI(100),NOOPEN(6,FILE='ht2.TXT')OPEN(8,FILE='ht.txt',STATUS='NEW)1 READ(6,*)NOIF(NO.EQ.0)STOPWRITE(8,'(/9X,A5,I3,A1)')'(NO=',NO,')'CALL READ(NJ,N,NEL,NM,NPJ,NPF,JN,X,Y,JE,JEAI,EAI,JPJ,PJ,JPF,PF)DO 5 I=1,NP(I) =0.0

16、DO 5 J=1,N5 K(I,J)=0.0DO 10 IE=1,NELCALL MKE(KE,IE,JE,JEAI,EAI,X,Y,AL)CALL MR(R,IE,JE,X,Y)CALL MAKE(KE,R,AKE)CALL CALM(M,IE,JN,JE)CALL MK(K,AKE,M)CONTINUEDO 20 IP=1,NPFCALL MR(R,JPF(1,IP),JE,X,Y)CALL TRAN(R,RT)CALL PE(FE,IP,JPF,PF,AL)CALL MULV6(RT,FE,AFE)CALL CALM(M,JPF(1,IP),JN,JE)CALL MF(P,AFE,M)C

17、ONTINUEDO 30 I=1,NPJP(JPJ(I)=P(JPJ(I)+PJ(I)CALL SOLV(K,P,D,N)WRITE(8,'(/2(26(1H*),A)')'RESULTS OF CALCULATION'WRITE(8,'(/28X,A)')'NODEL DISPLACEMENT'WRITE(8,40)FORMAT(9X,'NO.N',4X,'X-DISPLACEMENT',2X,#'Y-DISPLACEMENT',3X,'ANG.ROT.(RAD)'

18、)DO 60 KK=1,NJ10203040DO 50 II=1,3F(II)=0.0I1=JN(II,KK)50 IF(I1.GT.0)F(II)=D(I1)60 WRITE(8,70)KK,F(1),F(2),F(3)70 FORMAT(4X,I8,2X,3G16.5)WRITE(8,'(/30X,A)')'ELEMANT FORCES'WRITE(8,80)80 FORMAT(2X,'NO.E',4X,'N(1)',9X,'Q(1)',9X,'M(1)',#9X,'N(2)',

19、9X,'Q(2)',9X,'M(2)')DO 120 IE=1,NELCALL MADE(IE,JN,JE,D,ADE)CALL MKE(KE,IE,JE,JEAI,EAI,X,Y,AL)CALL MR(R,IE,JE,X,Y)CALL MULV6(R,ADE,DE)CALL MULV6(KE,DE,FF)DO 100 IP=1,NPFIF(JPF(1,IP).EQ.IE)THENCALL PE(FE,IP,JPF,PF,AL)DO 90 I=1,690 FF(I)=FF(I)-FE(I)ENDIF100 CONTINUEWRITE(8,110)IE,(FF(I

20、),I=1,6)110FORMAT(I5,2X,6G13.7)120 CONTINUEGOTO 1ENDSUBROUTINE READ(NJ,N,NEL,NM,NPJ,NPF,JN,X,Y,JE,JEAI,EAI,#JPJ,PJ,JPF,PF)REAL X(100),Y(100),EAI(3,100),PJ(100),PF(2,100)INTEGER JE(2,100),JN(3,100),JPJ(100),JPF(2,100),JEAI(100),#TITLE(20)READ(6,'(20A4)')(TITLE(I),I=1,20)WRITE(8,'(/9X,20A4

21、)')TITLEREAD(6,*)NJ,N,NEL,NM,NPJ,NPFWRITE(8,'(/3(8X,A4,1H:,I2)')'NJ=',NJ,#'N=',N,'NE=',NEL,'NM=',NM,'NPJ=',NPJ,'NPF=',NPFWRITE(8,10)10 FORMAT(/8X,'NO.N (1) (2) (3)',10X,'X',9X,'Y')READ(6,*)(JN(J,I),J=1,3),X(I),Y(I),I

22、=1,NJ)DO 20 I=1,NJ20 WRITE(8,'(8X,1H(,I2,1H),3I6,4X,2F10.3)')I,JN(1,I),JN(2,I),#JN(3,I),X(I),Y(I)READ(6,*)(JE(1,I),JE(2,I),JEAI(I),I=1,NEL)253040506070WRITE(8,25)FORMAT(/7X,'NO.E (1) (2) NO.MAT NO.E (1) (2) NO.MAT')N2=(NEL+1)/2DO 30 I=1,N2-1WRITE(8,40)I,(JE(J,I),J=1,2),JEAI(I),#I+N2,

23、(JE(J,I+N2),J=1,2),JEAI(I+N2)IF(N2*2.NE.NEL)WRITE(8,40)N2,JE(1,N2),JE(2,N2),JEAI(N2)IF(N2*2.EQ.NEL)WRITE(8,40)N2,JE(1,N2),JE(2,N2),JEAI(N2),#NEL,JE(1,NEL),JE(2,NEL),JEAI(NEL)FORMAT(4X,4I6,2X,4I6)READ(6,*)(EAI(I,J),I=1,3),J=1,NM)WRITE(8,50)(J,(EAI(I,J),I=1,3),J=1,NM)FORMAT(/3X,'NO.MAT',6X,

24、9;ELASTIC MODULUS',5X,#'AREA',7X,'MOMENT OF INERTIA'/(I6,5X,3G16.4)IF(NPJ.EQ.0)GOTO 70WRITE(8,'(/20X,12H NODEL LOADS)')WRITE(8,'(16XA)')'NO.DISP. VALUE'READ(6,*)(JPJ(I),PJ(I),I=1,NPJ)DO 60 I=1,NPJWRITE(8,'(14X,I7,F16.3)')JPJ(I),PJ(I)CONTINUEIF(NPF.

25、EQ.0)GOTO 100WRITE(8,'(/20X,16H NON-NODEL LOADS)')WRITE(8,'(7X,A,8X,A,9X,A)')'NO.E NO.LOAD.MODEL','A','C'READ(6,*)(JPF(1,I),JPF(2,I),PF(1,I),PF(2,I),I=1,NPF)DO 80 I=1,NPF80 WRITE(8,90)(JPF(J,I),J=1,2),PF(1,I),PF(2,I)90FORMAT(6X,I3,8X,I4,8X,2F10.3)100 CONTINUER

26、ETURNENDSUBROUTINE MKE(KE,IE,JE,JEAI,EAI,X,Y,AL)REAL KE(6,6),X(100),Y(100),EAI(3,100),AL(100),LINTEGER JE(2,100),JEAI(100)II=JE(1,IE)JJ=JE(2,IE)MT=JEAI(IE)L=SQRT(X(JJ)-X(II)*2+(Y(JJ)-Y(II)*2)AL(IE)=LA1=EAI(1,MT)*EAI(2,MT)/LA2=EAI(1,MT)*EAI(3,MT)/L*3A3=EAI(1,MT)*EAI(3,MT)/L*2A4=EAI(1,MT)*EAI(3,MT)/LK

27、E(1,1)=A1KE(1,4)=-A1KE(2,2)=12*A2KE(2,3)=6*A3KE(2,5)=-12*A2KE(2,6)=6*A3KE(3,3)=4*A4KE(3,5)=-6*A3KE(3,6)=2*A4KE(4,4)=A1KE(5,5)=12*A2KE(5,6)=-6*A3KE(6,6)=4*A4DO 10 I=1,6DO 10 K=1,610 KE(K,I)=KE(I,K)RETURNENDSUBROUTINE MR(R,IE,JE,X,Y)REAL R(6,6),X(100),Y(100),L,CX,CYINTEGER JE(2,100)I=JE(1,IE)J=JE(2,IE

28、)L=SQRT(X(J)-X(I)*2+(Y(J)-Y(I)*2)CX=(X(J)-X(I)/LCY=(Y(J)-Y(I)/LDO 10 J=1,6DO 10 I=1,610 R(I,J)=0.0DO 20 I=1,4,3R(I,I)=CXR(I,I+1)=CYR(I+1,I)=-CYR(I+1,I+1)=CX20 R(I+2,I+2)=1RETURNENDSUBROUTINE MAKE(KE,R,AKE)REAL KE(6,6),R(6,6),RT(6,6),TMP(6,6),AKE(6,6)CALL TRAN(R,RT)CALL MULV(RT,KE,TMP)CALL MULV(TMP,R

29、,AKE)RETURNENDSUBROUTINE CALM(M,IE,JN,JE)INTEGER M(6),JN(3,100),JE(2,100),IEDO 10 I=1,3M(I)=JN(I,JE(1,IE)10 M(I+3)=JN(I,JE(2,IE)RETURNENDSUBROUTINE MK(K,AKE,M)REAL K(200,200),AKE(6,6)INTEGER M(6)DO 10 I=1,6DO 10 J=1,6IF(M(I).NE.0.AND.M(J).NE.0)#K(M(I),M(J)尸K(M(I),M(J)+AKE(I,J)10 CONTINUERETURNENDSUB

30、ROUTINE PE(FE,IP,JPF,PF,AL)REAL FE(6),PF(2,100),AL(100),LINTEGER JPF(2,100)A=PF(1,IP)C=PF(2,IP)L=AL(JPF(1,IP)510203040IND=JPF(2,IP)DO 5 I=1,6FE(I)=0.0GOTO(10,20,30,40,50,60),INDFE(2)=(7.*A/20.+3.*C/20.)*LFE(3)=(A/20.+C/30.)*L*2FE(5)=(3.*A/20.+7*C/20.)*LFE(6)=-(A/30.+C/20.)*L*2RETURNFE(5)=A*C*3*(2.*L-C)/2./L*3FE(2)=A*C-FE(5)FE(3)=A*C*2*(6.*L*L-8.*C*L+3.*C*C)/12./L/LFE(6)=-A*C*3*(4.*L-3.*C)/12./L/LRETURNFE(2)=A*(L-C)*2*(L+2.*C)/L*3FE(3)=A*C*(L-C)*2/L*2FE(5)=A-FE(2)FE(6)=-A*C*2*(L-C)/L*2RETURNFE(2)=-6.*A*C*(L-C)/L*3FE(3)=A*(L-C)*(L-3.*C)/L*2FE(5)=-FE(2)FE(6)=A