有限元程序設(shè)計001

1、有限元程序設(shè)計平面三角形結(jié)點有限元程序1、程序功能本程序能計算彈性力學的平面應(yīng)力問題和平面應(yīng)變問題；能考慮自重和結(jié)點集中力兩種荷載的作用，在計算自重時軸取垂直向上為正； 能處理非零已知位移，如支座沉降的作用。主要輸出的內(nèi)容包括：結(jié)點位移、單元應(yīng)力、主應(yīng)力、第一主應(yīng)力與軸的夾角以及約束結(jié)點的支座反力。程序采用Visual Fortran 5.0編制而成，輸入數(shù)據(jù)全部采用自由格式。子程序的功能如下：INPUT輸入結(jié)點坐標、單元信息和材料參數(shù)；MR形成結(jié)點自由度序號矩陣；FORMMA形成指標矩陣MA（N）并調(diào)用其他功能子程序，相當于主控程序；DIV取出單元的3個結(jié)點號碼和該單元的材料號并計算單元的，

2、等；MGK形成整體勁度矩陣并按一維存放在SK（NH）中；LOAD形成整體結(jié)點荷載列陣F；OUTPUT輸出結(jié)點位移或結(jié)點荷載；TREAT由于有非零已知位移，對K和F進行處理；DECOMP整體勁度矩陣的分解運算；FOBA前代、回代求出未知結(jié)點位移；ERFAC計算約束結(jié)點的支座反力；KRS計算單元勁度矩陣中的子塊Krs。2、程序流程及框圖啟 動輸入數(shù)據(jù)信息 形成K 形成R解K R，輸出 計算，輸出應(yīng)力 停 機圖-1 程序流程圖圖-2 程序框圖3、輸入數(shù)據(jù)及變量說明當程序開始運行時，按屏幕提示，鍵入數(shù)據(jù)文件的名字。在運行程序之前，必須根據(jù)程序中輸入要求建立一個存放原始數(shù)據(jù)的文件，這個文件的名字由少于8

3、個字符或數(shù)字組成。數(shù)據(jù)文件包括如下內(nèi)容：總控信息，共一條，9個數(shù)據(jù)NP，NE，NM，NR，NI，NL，NG，ND，NCNP結(jié)點總數(shù)；NE單元總數(shù)；NM材料類型總數(shù)；NR約束結(jié)點總數(shù)；NI問題類型標識，0為平面應(yīng)力問題，1為平面應(yīng)變問題；NL受荷載作用的結(jié)點的數(shù)目；NG考慮自重作用為1，不計自重為0；ND非零已知位移結(jié)點的數(shù)目；NC要計算支座約束反力的結(jié)點數(shù)目。材料信息，共NM條，每條依次輸入EO，VO，W，tEO彈性模量（）；VO泊松比；W材料容重（）；t單元厚度（m）。這些信息都存放在數(shù)組AE（4，NM）中。坐標信息，共NP條，每條依次輸入IP，X，YIP結(jié)點號；X，Y分別為結(jié)點的坐標和坐標

4、。坐標信息存放在數(shù)組X(2，NP)中。單元信息，共NE條，每條依次輸入JE，L ，Io，Jo，MoJE單元號；L為該單元的材料類型號。Io，Jo，Mo分別為該單元i，j，m的整體編碼。單元信息存放在數(shù)組MEO(4,NE)中。約束信息共NR條，每條依次輸入一個數(shù) IP IP結(jié)點號；，分別為該結(jié)點的約束情況，如果方向受約束時填0，如果自由則填1。荷載信息，共NL條，每條依次輸入IP，IP結(jié)點號；，分別為該結(jié)點的，方向的荷載分量（kN）。結(jié)點號存放在數(shù)組NF（NL）中，結(jié)點荷載分量存放在數(shù)組FV（2，NL）中。若ND0，輸入非零已知位移信息，共ND條，每條依次輸入IP，IP結(jié)點號；，分別為該結(jié)點，方

5、向已知位移分量（m），若其中某方向為自由，則其相應(yīng)分量為0。結(jié)點號存放在數(shù)NDI（ND）中，已知位移分量存放在數(shù)組DV(2，ND)中。支座反力信息，共NC條，每條依次輸入IP，M1，M2，M3，M4IP 支座結(jié)點號；M1，M2，M3，M4 為與該支座結(jié)點相關(guān)的單元號，若不足4個，則用0補充。支座結(jié)點號存放在數(shù)組NCI(NC)中，相關(guān)單元號存放在數(shù)組NCE(4，NC)中。以上數(shù)據(jù)須按如上順序存放在數(shù)據(jù)文件中。除此之外，程序中還用到其他一些主要變量和數(shù)組，說明如下：N結(jié)構(gòu)自由度總數(shù)；NH 按一維存貯的整體勁度矩陣的總?cè)萘?；NX 最大半帶寬；SK(10000)維存貯的勁度矩陣；R(1000)開始存放

6、等效結(jié)點荷載，求解方程以后，用來存放結(jié)點位移；B(6)存放單元應(yīng)力， ，；MA(1000)主元素序號指標矩陣；JR(2，500)結(jié)點自由度序號矩陣；ME(3)存放單元結(jié)點i，j，m的整體編碼；NN(6)單元結(jié)點自由度序號；BI(3)，CI(3)單元勁度矩陣計算公式中的，和，；S三角形單元的面積；H11，H12，H21，H22單元勁度矩陣中子塊Krs的4個元素。4、源程序C FINITE ELEMENT PROGRAM FOR TWO DIMENSIONALC TRIANGLE ELEMENTC DIMENSION K(800000),COOR(2,3000),AE(4,11),* MEL(5,

7、2000),MA(6000) CHARACTER*32 dat COMMON /CA/ NP,NE,NM,NR,NI,NL,NG,ND,NC WRITE(*,300) 300 FORMAT(/ * :*/+PLEASE INPUT FILE NAME OF DATA) READ(*,*)data OPEN (4,FILE=data,STATUS=OLD) OPEN (7,FILE=OUT,STATUS=UNKNOWN) READ (4,*) NP,NE,NM,NR,NI,NL,NG,ND,NC C WRITE (*,400) NP,NE,NM,NR,NI,NL,NG,ND,NC C WRITE

8、 (7,400) NP,NE,NM,NR,NI,NL,NG,ND,NC CALL INPUT (JR,COOR,MEL,AE) CALL CBAND (MA,JR,MEL) CALL SK0(SK,R,COOR,MEL,MA,JR,AE) CALL LOAD (COOR,MEL,R,JR,AE) IF (ND.GT.0) CALL TREAT (SK,MA,JR,R) CALL DECOMP (SK,MA) CALL FOBA (SK,MA,R) WRITE(*,650) WRITE(7,650) CALL OUTPUT(JR,R) WRITE(*,700) WRITE(7,700) CALL

9、 CES (COOR,MEL,JR,R,AE) IF(NC.GT.0) call ERFAC (COOR,MEL,JR,R,AE) 400 FORMAT (/2X,NP=,I3,2X,NE=,I3,2X,NM= *,I3,2X,NR=,I3,2X,NI=I3,2X,NL=,I3,2X, *NG=,I3,2X,ND=,I3,2X,NC=,I3)500 FORMAT(1X,TOTAL DEGREES OF FREEDOM N=, *I4,1X,TOTAL-STORAGE ,NH=,I5,1X, *MAX-SEMI-BANDWIDTH MX=,I3)550 FORMAT(/20X,TOTAL STO

10、RAGE IS *GREATER THAN 50000)600 FORMAT(30X,NODAL FORCES/8X,NODE, *11X,X-COMP.,14X,Y-COMP.)650 FORMAT(/30X,NODAL DISPLACEMENTS/8X, *NODE,13X,X-COMP.,12X,Y-COMP.)700 FORMAT(/30X,ELEMENT STRESSES/5X, *ELEMENT,5X,X-STRESS,3X,Y-STRESS, *2X,XY-STRESS,1X,MAX-STRESS,1X, *MIN-STRESS,6X,ANGLE/) STOP END C * S

11、UBROUTINE KRS (BR,BS,CR,CS) COMMON /CB/ EO,VO,W,T,A,H11,H12,H21,H22 *,ME(3),BI(3),CI(3) ET=EO*T/(1.0-VO*VO)/A/4.0 V=(1.0-VO)/2.0 H11=ET*(BR*BS+V*CR*CS) H12=ET*(VO*BR*CS+V*BS*CR) H21=ET*(VO*CR*BS+V*BR*CS) H22=ET*(CR*CS+V*BR*BS) RETURN END C * SUBROUTINE INPUT (JR,COOR,MEL,AE) DIMENSION JR(2,*),COOR(2

12、,*),AE(4,*),MEL(3,*) COMMON /CA/ NP,NE,NM,NR COMMON /CC/ N,MX,NH DO 70 I=1,NP READ(4,*) IP,X,Y COOR(1,IP)=X COOR(2,IP)=Y 70 CONTINUE DO 11 J=1,NE READ(4,*)NEE,NME,(MEL(I,NEE),I=1,3) MEL(3,NEE)=NME11 CONTINUE DO 10 I=1,NP DO 10 J=1,2 10 JR(J,I)=1 DO 20 I=1,NR READ(4,*) IP,IX,IY JR(1,IP)=IX JR(2,IP)=I

13、Y 20 CONTINUE 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 55 J=1,NM READ (4,*)jj,(AE(I,jj),I=1,4)C WRITE(*,910)jj,(AE(I,jj),I=1,4)IF(NI.eq.1)then AE(1,jj)=AE(1,jj)/(1.0-AE(2,jj)*AE(2,jj) AE(2,jj)=AE(2,jj)/(1.0-AE(2,jj) endif55 CONTINUE910 FORMAT (/20X,MATERIAL

14、 PROPERTIES/(3X,I5,4(1x,E8.3) RETURN ENDC * SUBROUTINE CBAND (MA,JR,MEL) DIMENSION MA(*),JR(2,*),MEL(3,*),NN(6) COMMON /CA/ NP,NE,NM,NR COMMON /CC/ N,MX,NH DO 65 I=1,N65 MA(I)=0 DO 90 IE=1,NE DO 75 K=1,3 IEK=MEL(K,IE) DO 95 M=1,2 JJ=2*(K-1)+M NN(JJ)=JR(M,IEK)95 CONTINUE75 CONTINUE L=N DO 80 I=1,6 NN

15、I=NN(I) IF(NNI.EQ.0) GO TO 80 IF(NNI.LT.L) L=NNI 80 CONTINUE DO 85 M=1,6 JP=NN(M) IF(JP.EQ.0) GO TO 85 JPL=JP-L+1 IF(JPL.GT.MA(JP) MA(JP)=JPL 85 CONTINUE 90 CONTINUE MX=0 MA(1)=1 DO 10 I=2,N IF(MA(I).GT.MX) MX=MA(I) MA(I)=MA(I)+MA(I-1) 10 CONTINUE NH=MA(N) WRITE (*,500) N,MX,NH WRITE (7,500) N,MX,NH

16、500 FORMAT (/5X,FREEDOM N= *,I5,3X,SEMI-BANDWI. MX=,I5,3X, * STORAGE NH=,I7) RETURN END C* SUBROUTINE SK0(SK,R,COOR,MEL,MA,JR,AE) DIMENSION AE(4,*),COOR(2,*),MEL(3,*),JR(2,*),R(*), *MA(*),SK(*),SKE(6,6),NN(6) COMMON /CA/ NP,NE,NM,NR,NI,NL,NG,ND,NC COMMON /CB/ EO,VO,W,T,A,H11,H12,H21,H22, *ME(3),BI(3

17、),CI(3) COMMON /CC/ N,NH DO 10 I=1,NH 10 SK(I)=0.0 DO 70 IE=1,NE CALL DIV (IE,COOR,MEL,AE) DO 30 I=1,3 DO 30 J=1,3 CALL KRS (BI(I),BI(J),CI(I),CI(J) SKE(2*I-1,2*J-1)=H11 SKE(2*I-1,2*J)=H12 SKE(2*I,2*J-1)=H21 SKE(2*I,2*J)=H22 30 CONTINUE DO 40 I=1,3 J2=ME(I) DO 40 J=1,2 J3=2*(I-1)+J NN(J3)=JR(J,J2) 4

18、0 CONTINUE DO 60 I=1,6 DO 60 J=1,6 IF(NN(J).EQ.0.OR.NN(I).LT.NN(J) GO TO 60 JJ=NN(I) JK=NN(J) JL=MA(JJ) JM=JJ-JK JN=JL-JM SK(JN)=SK(JN)+SKE(I,J) 60 CONTINUE 70 CONTINUE c WRITE(0,500) (SK(I),I=1,20) 500 FORMAT(/10X,SK=/(6F12.5) RETURN END C * SUBROUTINE LOAD (COOR,MEL,R,JR,AE) DIMENSION AE(4,*),COOR

19、(2,*),MEL(3,*),R(*),JR(2,*), * NF(50),FV(2,50) COMMON /CA/ NP,NE,NM,NR,NI,NL,NG,ND,NC COMMON /CB/ EO,VO,W,T,A,A1(4),ME(3),BB(6) COMMON /CC/ N,NH DO 10 I=1,N 10 R(I)=0.0 IF(NG) 70,70,30 30 DO 60 IE=1,NE CALL DIV (IE,COOR,MEL,AE) DO 50 I=1,3 J2=ME(I) J3=JR(2,J2) IF(J3) 50,50,40 40 R(J3)=R(J3)-T*W*A/3.

20、0 50 CONTINUE 60 CONTINUE 70 IF(NL) 110,110,80 80 READ(4,*) (NF(I),I=1,NL) READ(4,*) (FV(I,J),I=1,2),J=1,NL) C WRITE(*,500) (NF(I),I=1,NL) C WRITE(7,500) (NF(I),I=1,NL)C WRITE(*,600) (FV(I,J),I=1,2),J=1,NL) C WRITE(7,600) (FV(I,J),I=1,2),J=1,NL) DO 100 I=1,NL JJ=NF(I) J=JR(1,JJ) M=JR(2,JJ) IF (J.GT.

21、0) R(J)=R(J)+FV(1,I) IF (M.GT.0) R(M)=R(M)+FV(2,I) 100 CONTINUE 110 RETURN 500 FORMAT(/20X,NODES OF APPLIED LOAD*NF= */(1X,10I8)600 FORMAT(/30X,LUMPED-LOADS*FV= */(5X,5F15.3) END C * SUBROUTINE TREAT (SK,MA,JR,R) DIMENSION SK(*),MA(*),NDI(75),DV(2,75),JR(2,*),R(*) COMMON /CA/ NP,NE,NM,NR,NI,NL,NG,ND

22、,NC COMMON /CC/ N,NH READ(4,*) (NDI(J),J=1,ND) READ(4,*) (DV(I,J),I=1,2),J=1,ND) C WRITE(*,500) (NDI(J),J=1,ND) C WRITE(7,500) (NDI(J),J=1,ND)C WRITE(*,550) (DV(I,J),I=1,2),J=1,ND) C WRITE(7,550) (DV(I,J),I=1,2),J=1,ND) DO 20 I=1,ND DO 20 J=1,2 IF(DV(J,I) 10,20,10 10 JJ=NDI(I) L=JR(J,JJ) JN=MA(L) SK

23、(JN)=1.0E30 R(L)=DV(J,I)*1.0E30 20 CONTINUE 500 FORMAT(/25X,NODE NO.*NDI=/(1X,10I8)550 FORMAT(/25X,DISPLACEMENT-VALUES*DV=/ *(10X,6F10.6) RETURN END C * SUBROUTINE DECOMP (SK,MA) DIMENSION SK(*),MA(*) COMMON /CC/ N,NH DO 50 I=2,N L=I-MA(I)+MA(I-1)+1 K=I-1 L1=L+1 IF (L1.GT.K) GO TO 30 DO 20 J=L1,K IJ

24、=MA(I)-I+J M=J-MA(J)+MA(J-1)+1 IF (L.GT.M) M=L MP=J-1 IF (M.GT.MP) GO TO 20 DO 10 LP=M,MP IP=MA(I)-I+LP JP=MA(J)-J+LP SK(IJ)=SK(IJ)-SK(IP)*SK(JP) 10 CONTINUE 20 CONTINUE 30 IF (L.GT.K) GO TO 50 DO 40 LP=L,K IP=MA(I)-I+LP LPP=MA(LP) SK(IP)=SK(IP)/SK(LPP) II=MA(I) SK(II)=SK(II)-SK(IP)*SK(IP)*SK(LPP) 40 CONTINUE 50 CONTINUE C WRITE(0,500) (SK(I),I=1,20) 500 FORMAT(/10X,SK=/(1X,6F12.4)