有限元程序設計.doc

1、.平面問題有限元程序設計理學院 學號 兔兔設計人 完成日期 一、 程序功能說明本程序適用于節(jié)點荷載作用下的桁架分析問題，當有節(jié)間荷載存在時可按照靜力等效原理將其轉化為節(jié)點荷載?？汕蠼馄矫骅旒茉陟o力荷載作用下的內力和位移。二、 框圖的設計開 始 輸 入 數(shù) 據(jù) 數(shù) 組 定 義計算各桿截面面積和半帶寬 調用形成單剛矩陣UNIT調用形成半帶寬存貯的結構原始剛度矩陣TOTAL有節(jié)點荷載否輸入節(jié)點荷載值，并將其送入相應的荷載列陣P（N）中考慮結構是否自重將桿自重引起的等效荷載疊加到P（N）中支座處理、解方程，并輸出U（N）、V（N）調用UNIT，求各單元桿端內力 結 束單元循環(huán)沒有有否是 三、 程序的標

2、識符及數(shù)組說明NPOIN 最大節(jié)點數(shù)NELEM 最大單元數(shù) NLOAD 節(jié)點的荷載總數(shù)NZERO 節(jié)點的約束位移總數(shù)WT 結構的自重EE 材料的彈性模量LL 一維數(shù)組，用于存放單元桿件的長度AA 一維數(shù)組，用于存放單元桿件的面積COORD 節(jié)點坐標數(shù)組LNODE 單元節(jié)點數(shù)組BH 二維數(shù)組，用于存放單元截面尺寸NRES 二維數(shù)組，用于存放約束的位移值JP 二維數(shù)組，用于存放節(jié)點的荷載值ESTIF 四維數(shù)組，用于存放整體坐標系下的單元剛度矩陣ASTIF 二十維數(shù)組，用于存放半帶寬結構原始剛度矩陣P 用于存放節(jié)點的荷載列陣U 用于存放節(jié)點x方向的位移值V 用于存放節(jié)點y方向的位移值四、 源程序IN

3、TEGER E,NELEM,Z,H REAL LL,ESTIF,ASTIF,JP DIMENSION COORD(3,2),LNODE(3,2),AA(200),BH(3,2),RES(3,2), &LL(200),ESTIF(4,4),ASTIF(400,20),JP(1,2),P(400),U(200), &V(200) OPEN(2,FILE='D:NMXJIA.DAT',STATUS='NEW') C 輸入已知數(shù)據(jù)DATA NPOIN,NELEM,NJP,NRES,EE,WT/3,3,1,3,21000,0/DATA COORD/0,6,

4、0,0,0,6/DATA LNODE/1,2,1,2,3,3/DATA BH/3*2,3*10/DATA RES/3*0,1,2,4/C 計算各單元面積DO 200 E=1,NELEMAA(E)=BH(E,1)*BH(E,2) CALL UNIT(E,EE,COORD,LNODE,AA,ESTIF,LL,CX,CY)200 CONTINUE C計算半帶寬 L2=2*NPOIN NHBW=0 DO 210 E=1,NELEM M=ABS(LNODE(E,1)-LNODE(E,2) IF(NHBW.LT.M) NHBW=M210 CONTINUEWRITE(2,*) '半帶寬' N

5、HBW=2*(NHBW+1) WRITE(2,220) NHBW220 FORMAT(1X,'NHBW=', I2)C單元循環(huán) DO 300 I1=1,L2 DO 300 J1=1,NHBW300 ASTIF(I1,J1)=0.0 DO 400 E=1,NELEM CALL UNIT(E,EE,COORD,LNODE,AA,ESTIF,LL,CX,CY) CALL TOTAL(E,LNODE,ESTIF,ASTIF) 400 CONTINUE DO 560 N=1,L2560 P(N)=0.0 IF(NJP.EQ.0) GOTO 650DATA JP/10,5/ DO 630

6、K1=1,NJP NN=JP(K1,2)+0.1 630 P(NN)=JP(K1,1) 650 IF(WT.LE.0.0) GOTO 750 DO 700 E=1,NELEM N1=LNODE(E,1) N2=LNODE(E,2) P(2*N1)=P(2*N1)-WT*AA(E)*LL(E)/2.0 P(2*N2)=P(2*N2)-WT*AA(E)*LL(E)/2.0700 CONTINUE WRITE(2,710)710 FORMAT(/4X,'荷載總數(shù)',8X,'水平荷載',8X,'鉛垂荷載') DO 730 K=1,NO730 WRITE(

7、2,740) K,P(2*K-1),P(2*K)740 FORMAT(4X,I2,8X,F8.3,8X,F8.3) 750 DO 800 I1=1,NRES Z=RES(I1,2)+1E-5 ASTIF(Z,1)=ASTIF(Z,1)*1E8 P(Z)=ASTIF(Z,1)*RES(I1,1)800 CONTINUE DO 850 K1=1,L2-1 IF(L2.GT.(K1+NHBW-1) THEN IM=K1+NHBW-1 ELSE IM=L2 ENDIF DO 850 I1=K1+1,IM L1=I1-K1+1 C1=ASTIF(K1,L1)/ASTIF(K1,1) DO 830 J1=

8、1,NHBW-L1+1 MM=J1+I1-K1 ASTIF(I1,J1)=ASTIF(I1,J1)-C1*ASTIF(K1,MM)830 CONTINUE P(I1)=P(I1)-C1*P(K1)850 CONTINUE P(L2)=P(L2)/ASTIF(L2,1) DO 900 I1=L2-1,1,-1 IF(NHBW.GT.(L2-I1+1) THEN JM=L2-I1+1 ELSE JM=NHBW ENDIF DO 880 J1=2,JM H=J1+I1-1 P(I1)=P(I1)-ASTIF(I1,J1)*P(H)880 CONTINUE P(I1)=P(I1)/ASTIF(I1,1

9、)900 CONTINUE WRITE(2,910)910 FORMAT(/10X,'節(jié)點位移',10X,'水平位移',10X,'鉛垂位移'/) DO 930 N=1,NO U(N)=P(2*N-1) V(N)=P(2*N)930 WRITE(2,950) N,U(N),V(N)950 FORMAT(15X,I2,6X,F12.7,6X,F12.7) WRITE(2,970)970 FORMAT(/4X,'單元號',8X,'節(jié)點號',8X,'N(KN)',8X,'Q(KN)'/) D

10、O 980 E=1,NELEM CALL UNIT(E,EE,COORD,LNODE,AA,ESTIF,LL,CX,CY) N1=LNODE(E,1) N2=LNODE(E,2) ULNODE=U(N1)-U(N2) VLNODE=V(N1)-V(N2) D1=ESTIF(1,1)*ULNODE+ESTIF(1,2)*VLNODE D2=ESTIF(1,2)*ULNODE+ESTIF(2,2)*VLNODE FI=CX*D1+CY*D2 FJ=-FI TI=-CY*D1+CX*D2 TJ=-TI WRITE(2,990) E,N1,FI,TI,N2,FJ,TJ990 FORMAT(4X,I2,

11、12X,I2,8X,F8.4,8X,F8.4/18X,I2,8X,F8.4,8X,F8.4)980 CONTINUE WRITE(2,1000)1000FORMAT(/28X,'結束',/15X,35('*')/) STOP END 子程序TOTAL形成總剛度矩陣SUBROUTINE TOTAL(E,LNODE,ESTIF,ASTIF) INTEGER E,DH,ZL,DL REAL ESTIF,ASTIF DIMENSION LNODE(3,2),ESTIF(4,4),ASTIF(400,20) DO 40 I1=1,2 DO 40 II=1,2 KH=2*(

12、I1-1)+II DH=2*(LNODE(E,I1)-1)+II DO 40 J1=1,2 DO 40 JJ=1,2 KL=2*(J1-1)+JJ ZL=2*(LNODE(E,J1)-1)+JJ DL=ZL-DH+1 IF(DL.GT.0) ASTIF(DH,DL)=ASTIF(DH,DL)+ESTIF(KH,KL)40 CONTINUE RETURN END 子程序UNIT形成單剛SUBROUTINE UNIT(E,EE,COORD,LNODE,AA,ESTIF,LL,CX,CY) INTEGER E REAL LL,ESTIF DIMENSION COORD(3,2),LNODE(3,2)

13、,AA(200),LL(200),ESTIF(4,4) N1=LNODE(E,1) N2=LNODE(E,2) CX=COORD(N2,1)-COORD(N1,1) CY=COORD(N2,2)-COORD(N1,2) LL(E)=SQRT(CX*CX+CY*CY) CX=CX/LL(E) CY=CY/LL(E) EAL=EE*AA(E)/LL(E) ESTIF(1,1)=EAL*CX*CX ESTIF(1,2)=EAL*CX*CY ESTIF(2,2)=EAL*CY*CY ESTIF(2,1)=ESTIF(1,2) DO 10 I=1,2 DO 10 J=1,2 ESTIF(I,J+2)=-

14、ESTIF(I,J) ESTIF(I+2,J)=-ESTIF(I,J)10 ESTIF(I+2,J+2)=ESTIF(I,J) RETURN END五、 算例如圖1所示桁架，已知桿件材料的彈性模量,桿件截面高度H=10cm，截面的寬度為b=2cm，不計各桿的自重，求在荷載作用下，各桿的軸力。圖1桁架例圖計算模型的輸入數(shù)據(jù)有以下各量：NPOINNRESNELEMNLOADYOGWT3331210000.0節(jié)點單元12LNODE(1,I)12LNODE(2,I)23LNODE(3,I)13編號數(shù)組123RES(I,1)0.00.00.0RES(I,2)124節(jié)點 I123COORD(I,1)0.06.00.0COORD(I，2)0.00.06.0JP(I,1)10.0JP(I,2)5單元號I123BH(I,1)101010BH(I,