有限元大作業(yè)matlab---課程設(shè)計報告例子

1、唸裁、丄丈才Y a n s h a n U ni versify有限元大作業(yè)程序設(shè)計學(xué)校：天津大學(xué)院系：建筑工程與力學(xué)學(xué)院專業(yè)：01級工程力學(xué)姓名：劉秀學(xué)號：指導(dǎo)老師：連續(xù)體平面問題的有限元程序分析題目:如圖所示的正方形薄板四周受均勻載荷的作用，該結(jié)構(gòu)在邊界上受正向分布壓力，Pi"%，同時在沿對角線y軸上受一對集中壓 力，載荷為2KN若取板厚t=1，泊松比v = 0。分析過程:由于連續(xù)平板的對稱性，只需要取其在第一象限的四分之一部分 參加分析，然后人為作出一些輔助線將平板“分割”成若干部分，再 為每個部分選擇分析單元。采用將此模型化分為4個全等的直角三角 型單元。利用其對稱性，四分之

2、一部分的邊界約束，載荷可等效如圖 所示。程序原理及實現(xiàn):用FORTRAIN序的實現(xiàn)。由節(jié)點信息文件NODE.IN和單元信息文件ELEMENT.IN經(jīng)過計算分析后輸出一個一般性的文件 DATA.OU。模型基本信息由文件為BASIC.IN生成。該程序的特點如下：問題類型：可用于計算彈性力學(xué)平面問題和平面應(yīng)變問題單元類型：采用常應(yīng)變?nèi)切螁卧灰颇Ｊ剑河糜镁€性位移模式載荷類型：節(jié)點載荷，非節(jié)點載荷應(yīng)先換算為等效節(jié)點載荷材料性質(zhì)：彈性體由單一的均勻材料組成約束方式：為“ 0”位移固定約束，為保證無剛體位移，彈性體至少應(yīng)有對三個自由度的獨立約束方程求解：針對半帶寬剛度方程的 Gauss消元法輸入文件：由

3、手工生成節(jié)點信息文件NODE.IN,和單元信息文件ELEMENT.IN結(jié)果文件：輸出一般的結(jié)果文件 DATA.OUT程序的原理如框圖：1) 主要變量：ID:問題類型碼，ID = 1時為平面應(yīng)力問題，ID=2時為平面應(yīng)變問題N_NOD：E 節(jié)點個數(shù)N_LOA:D 節(jié)點載荷個數(shù)N_DOF 自由度，N_DOF=N_NODE*平面問題)N_ELE:單元個數(shù)N_BAN:D 矩陣半帶寬N_BC:有約束的節(jié)點個數(shù)PE:彈性模量PR:泊松比PT:厚度LJK_ELE(I,3) :單元節(jié)點編號數(shù)組， LJK_ELE(I,1) ， LJK_ELE(I,2) ，LJK_ELE(I,3) 分別放單元 I 的三個節(jié)點的整

4、體編號X(N_NODE), 丫(N_NODE)節(jié)點坐標(biāo)數(shù)組，X(I),Y(I) 分別存放節(jié)點I的x, y 坐標(biāo)值P_LJK(N_BC,3): 節(jié)點載荷數(shù)組， P_LJK(I,1) 表示第 I 個作用有節(jié)點載荷的節(jié)點的編號， P_LJK(I,2) ， P_LJK(I,3) 分別為該節(jié)點沿 x,y 方向的節(jié)點載荷數(shù)值A(chǔ)K(N_DOF,N_BAND)AKE(6,6) ：整體剛度矩陣單元剛度矩陣BB(3,6) ：位移 應(yīng)變轉(zhuǎn)換矩陣(二節(jié)點單兀的幾何矩陣)DD(3,3) ：彈性矩陣SS(3,6) ；RESULT_N(N_NOF)存放節(jié)點位移DISP_E(6): ：應(yīng)力矩陣節(jié)點載荷數(shù)組，存放節(jié)點載荷向量，

5、解方程后該矩陣單元的節(jié)點位移向量STS_ELE(N_ELE,3)STS_ND(N_NODE,3)單元的應(yīng)力分量節(jié)點的應(yīng)力分量2) 子程序說明:CLA_BB：計算單兀位移應(yīng)變關(guān)系矩陣READ_IN : 讀入數(shù)據(jù) BAND_KFORM_KE: 計算單元剛度矩陣CAL_ARE A:計算單元面積CLA_DD:計算單元彈性矩陣CAL_STS計算單元和節(jié)點應(yīng)力: 形成半帶寬的整體剛度矩陣FORM_:P 計算節(jié)點載荷DO_B:C 處理邊界條件SOLVE : 計算節(jié)點位移3) 文件管理：源程序文件：chengxu.for程序需讀入的數(shù)據(jù)文件：BASIC.IN , NODE.IN, ELEMENTN（需要手工生

6、成）程序輸出的數(shù)據(jù)文件:DATA.OUT（4）數(shù)據(jù)文件格式:需讀入的模型 基本信息文件BASIC.IN的格式如下表欄目格式說明實際需輸入的數(shù)據(jù)基本模型數(shù)據(jù)第1行，每兩個數(shù)之間用“，”號 隔開問題類型，單元個數(shù)，節(jié) 點個數(shù)，有約束的節(jié)點 數(shù)，有載何的節(jié)點數(shù)材料性質(zhì)第2行，每兩個數(shù)之間用“，”號 隔開彈性模量，泊松比，單元 厚度節(jié)點約束信息在材料性質(zhì)輸入行之后另起行，每 兩個數(shù)之間用“，”號隔開LJK_U（N_BC,3）位移約束的節(jié)點編號，該 節(jié)點x方向約束代碼，該 節(jié)點y方向代碼，節(jié)點何載信息在節(jié)點約束信息輸入行之后另起 行，每兩個數(shù)之間用“，”號隔開P IJK（N LOAD,3）載何作用的節(jié)點

7、編號，該 節(jié)點x主向載荷，該節(jié)點 y方向載荷，需讀入的節(jié)點信息文件NODE.IN的格式如下表欄目格式說明實際需輸入的數(shù)據(jù)節(jié)點信息母仃為一個節(jié)點的信息 （每行三個數(shù)，每兩個數(shù) 之間用空格或“，”分開）ND_ANSYS（N_NIDE）節(jié)點號，該節(jié)點的x坐標(biāo)， 該節(jié)點y方向坐標(biāo)需讀入的單元信息文件ELEMENT.IN勺格式如下表欄目格式說明實際需輸入的數(shù)據(jù)單元信息母仃為一個單兀的信息（每行有14個整型數(shù)， 前4個為單元節(jié)點編號， 對于3節(jié)點編號，第4個 節(jié)點編號與第3個節(jié)點編 號相同，后10個數(shù)無用， 可輸入“ 0 ”，每兩個整 型數(shù)之間用至少一個空 格分開）NE ANSYS（N ELE,14） 單

8、元的節(jié)點號1 （空格） 單元的節(jié)點號2 （空格） 單元的節(jié)點號3 （空格） 單元的節(jié)點號4 （空格） 0 （空格）0 （空格）0 （空 格）0 （空格）0 （空格） 0 （空格）0 （空格）0 （空 格）0 （空格）0輸出結(jié)果文件DATA.OU格式如下表欄目實際輸出的數(shù)據(jù)節(jié)點位移I RESULT_N(2*I_ 1) RESULT_N(2*I) 節(jié)點號 x 方向位移y方向位移單元應(yīng)力的三個分量IESTE_ELE(IE,1)STE_ELE(IE,2)STE_ELE(IE,3)單元號x方向應(yīng)力 y 方向應(yīng)力剪切應(yīng)力節(jié)點應(yīng)力的三個分量ISTS-ND(I,1)STS-ND(I,2)STS-ND(I,3)

9、節(jié)點號x方向應(yīng)力y方向應(yīng)力剪切應(yīng)力算例原始數(shù)據(jù)和程序分析:(1)模型基本信息文件BASIC.IN的數(shù)據(jù)為1,4,6,5,31.,0.,1.1,1,021,0,4,1,1,5,0,1,6,0,11,-0.5,-1.5,3.,-1.,-1,6,-0.5,-0.5(2)手工準(zhǔn)備的節(jié)點信息文件NODE.IN的數(shù)據(jù)為10.02.020.01.031.01.040.0.51.00.62.00.(3)手工準(zhǔn)備的單元信息文件ELEMENT.IN勺數(shù)據(jù)為12330000 11 11 0124550000 11 11 0253220000 11 11 0335660000 11 11 04(4)源原程序文件che

10、n gxu.for為：PROGRAM FEM2DDIMENSION IJK_ELE(500,3),X(500),Y(500),IJK_U(50,3),P_IJK(50,3),& RESULT_N(500),AK(500,100)DIMENSION STS_ELE(500,3),STS_ND(500,3)OPEN(4,FILE='BASICN')OPEN(5,FILE='NODEN')OPEN(6,FILE='ELEMENTN')0PEN(8,FILE='DATA.0UT')0PEN(9,FILE='F0R_P0ST

11、.DAT')READ(4,*)ID,N_ELE,N_N0DE,N_BC,N_L0ADIF(ID.EQ.1)WRITE(8,20)IF(ID.EQ.2)WRITE(8,25)20F0RMAT(/5X,'=PLANE STRESS PROBLEM=')25F0RMAT(/5X,'=PLANE STRAIN PROBLEM=')CALL READ_IN(ID,N_ELE,N_NODE,N_BC,N_BAND,N_LOAD,PE,PR,PT,& IJK_ELE,X,Y,IJK_U,P_IJK)CALL BAND_K(N_DOF,N_BAND,N_ELE,

12、IE,N_NODE,&IJK_ELE,X,Y,PE,PR,PT,AK)CALL FORM_P(N_ELE,N_NODE,N_LOAD,N_DOF,IJK_ELE,X,Y,P_IJK,&RESULT_N)CALL DO_BC(N_BC,N_BAND,N_DOF,IJK_U,AK,RESULT_N)CALL SOLVE(N_NODE,N_DOF,N_BAND,AK,RESULT_N)CALL CAL_STS(N_ELE,N_NODE,N_DOF,PE,PR,IJK_ELE,X,Y,RESULT_N,& STS_ELE,STS_ND)c to putout a data fi

13、leWRITE(9,70)REAL(N_NODE),REAL(N_ELE)70 FORMAT(2f9.4)WRITE(9,71)(X(I),Y(I),RESULT_N(2*I-1),RESULT_N(2*I),& STS_ND(I,1),STS_ND(I,2),STS_ND(I,3),I=1,N_NODE)71 FORMAT(7F9.4)WRITE(9,72)(REAL(IJK_ELE(I,1),REAL(IJK_ELE(I,2),& REAL(IJK_ELE(I,3),REAL(IJK_ELE(I,3),& STS_ELE(I,1),STS_ELE(I,2),STS_

14、ELE(I,3),I=1, N_ELE)72 FORMAT(7f9.4)cCL0SE(4)CL0SE(5)CL0SE(6)CL0SE(8)CL0SE(9)ENDcc to get the orig inal data in order to model the problemSUBROUTINE READ_IN(ID,N_ELE,N_NODE,N_BC,N_BAND,N_LOAD,PE,PR,& PT,IJK_ELE,X,Y,IJK_U,P_IJK)DIMENSION IJK_ELE(500,3),X(N_NODE),Y(N_NODE),IJK_U(N_BC,3),& P_IJ

15、K(N_LOAD,3),NE_ANSYS(N_ELE,14)REAL ND_ANSYS(N_N0DE,3)READ(4,*)PE,PR,PTREAD(4,*)(IJK_U(I,J),J=1,3),I=1,N_BC)READ(4,*)(P_IJK(I,J),J=1,3),I=1,N_LOAD)READ(5,*)(ND_ANSYS(I,J),J=1,3),I=1,N_NODE)READ(6,*)(NE_ANSYS(I,J),J=1,14),I=1,N_ELE)DO 10 I=1,N_NODEX(I)=ND_ANSYS(I,2)Y(I)=ND_ANSYS(I,3)10 CONTINUEDO 11 I

16、=1,N_ELEDO 11 J=1,3IJK_ELE(I,J)=NE_ANSYS(I,J)11 CONTINUEN_BAND=0DO 20 IE=1,N_ELEDO 20 I=1,3DO 20 J=1,3IW=IABS(IJK_ELE(IE,I)-IJK_ELE(IE,J)IF(N_BAND. LT.IW)N_BAND=IW20 CONTINUEN_BAND=(N_BAND+1)*2IF(ID.EQ.1) THENELSEPE=PE/(1.0-PR*PR)PR=PR/(1.0-PR)END IFRETURNENDcC to form the stiffness matrix of eleme

17、ntSUBROUTINE FORM_KE(IE,N_NODE,N_ELE,IJK_ELE,X,Y,PE,PR,PT,AKE)DIMENSION IJK_ELE(5OO,3),X(N_NODE),Y(N_NODE),BB(3,6),DD(3,3),&AKE(6,6), SS(6,6)CALL CAL_DD(PE,PR,DD)CALL CAL_BB(IE,N_NODE,N_ELE,IJK_ELE,X,Y,AE,BB)DO 10 I=1,3DO 10 J=1,6SS(I,J)=0.0DO 10 K=1,310 SS(I,J)=SS(I,J)+DD(I,K)*BB(K,J)DO 20 I=1,

18、6DO 20 J=1,6AKE(I,J)=0.0DO 20 K=1,320 AKE(I,J)=AKE(I,J)+SS(K,I)*BB(K,J)*AE*PTRETURNENDcc to form ban ded global stiffness matrixSUBROUTINE BAND_K(N_DOF,N_BAND,N_ELE,IE,N_NODE,IJK_ELE,X,Y,PE,& PR,PT,AK)DIMENSION IJK_ELE(500,3),X(N_NODE),Y(N_NODE),AKE(6,6),AK(500,100) N_DOF=2*N_NODEDO 40 I=1,N_DOF

19、DO 40 J=1,N_BAND40AK(I,J)=0DO 50 IE=1,N_ELECALL FORM_KE(IE,N_NODE,N_ELE,IJK_ELE,X,Y,PE,PR,PT,AKE)DO 50 I=1,3DO 50 II=1,2IH=2*(I-1)+IIIDH=2*(IJK_ELE(IE,I)-1)+IIDO 50 J=1,3DO 50 JJ=1,2IL=2*(J-1)+JJIZL=2*(IJK_ELE(IE,J)-1)+JJIDL=IZL-IDH+1IF(IDL. LE.0) THENELSEAK(IDH,IDL)=AK(IDH,IDL)+AKE(IH,IL)END IF50 C

20、ONTINUERETURNENDcc to calculate the area of eleme ntSUBROUTINE CAL_AREA(IE,N_NODE,IJK_ELE,X,Y,AE)DIMENSION IJK_ELE(500,3),X(N_NODE),Y(N_NODE)I=IJK_ELE(IE,1)J=IJK_ELE(IE,2)K=IJK_ELE(IE,3)XIJ=X(J)-X(I)YIJ=Y(J)-Y(I)XIK=X(K)-X(I)YIK=Y(K)-Y(I)AE=(XIJ*YIK-XIK*YIJ)/2.0RETURNEND cc to calculate the elastic

21、matrix of elementSUBROUTINE CAL_DD(PE,PR,DD)DIMENSION DD(3,3)DO 10 I=1,3DO 10 J=1,310 DD(I,J)=0.0DD(1,1)=PE/(1.0-PR*PR)DD(1,2)=PE*PR/(1.0-PR*PR)DD(2,1)=DD(1,2)DD(2,2)=DD(1,1)DD(3,3)=PE/(1.0+PR)*2.0)RETURNENDcc to calculate the strain-displacement matrix of elementSUBROUTINE CAL_BB(IE,N_NODE,N_ELE,IJ

22、K_ELE,X,Y,AE,BB)DIMENSION IJK_ELE(500,3),X(N_NODE),Y(N_NODE),BB(3,6)I=IJK_ELE(IE,1)J=IJK_ELE(IE,2)K=IJK_ELE(IE,3)DO 10 II=1,3DO 10 JJ=1,310 BB(II,JJ)=0.0BB(1,1)=Y(J)-Y(K)BB(1,3)=Y(K)-Y(I)BB(1,5)=Y(I)-Y(J)BB(2,2)=X(K)-X(J)BB(2,4)=X(I)-X(K)BB(2,6)=X(J)-X(I)BB(3,1)=BB(2,2)BB(3,2)=BB(1,1)BB(3,3)=BB(2,4)

23、BB(3,4)=BB(1,3)BB(3,5)=BB(2,6)BB(3,6)=BB(1,5)CALL CAL_AREA(IE,N_NODE,IJK_ELE,X,Y,AE)DO 20 I1=1,3DO 20 J1=1,620 BB(I1,J1)=BB(I1,J1)/(2.0*AE)RETURNENDcc to form the global load matrixSUBROUTINE FORM_P(N_ELE,N_NODE,N_LOAD,N_DOF,IJK_ELE,X,Y,P_IJK,&RESULT_N)DIMENSION IJK_ELE(500,3),X(N_NODE),Y(N_NODE

24、),P_IJK(N_LOAD,3),&RESULT_N(N_DOF)DO 10 I=1,N_DOF10 RESULT_N(l)=0.0DO 20 I=1,N_LOADII=P_IJK(I,1)RESULT_N(2*II-1)=P_IJK(I,2)20 RESULT_N(2*II)=P_IJK(I,3)RETURNENDcmethodc to deal with BC(u) (here only for fixed displacement) using "1-0" SUBROUTINE DO_BC_BC,N_BAND,N_DOF,IJK_U,AK,RESULT_N)

25、DIMENSION RESULT_N(N_DOF),IJK_U(N_BC,3),AK(5OO,1OO)DO 30 I=1,N_BCIR=IJK_U(I,1)DO 30 J=2,3IF(IJK_U(I,J).EQ.O)THENELSEII=2*IR+J-3AK(II,1)=1.0RESULT_N(II)=0.0DO 10 JJ=2,N_BAND10AK(II,JJ)=0.0DO 20 JJ=2,II20AK(II-JJ+1,JJ)=0.0END IF30 CONTINUERETURNENDcc to solve the ban ded FEM equati on by GAUSS elim in

26、 ati onSUBROUTINE SOLVE(N_NODE,N_DOF,N_BAND,AK,RESULT_N) DIMENSION RESULT_N(N_DOF),AK(500,100)DO 20 K=1,N_DOF-1IF(N_DOF.GT.K+N_BAND-1)IM=K+N_BAND-1IF(N_DOF.LE.K+N_BAND-1)IM=N_DOFDO 20 I=K+1,IML=I-K+1C=AK(K,L)/AK(K,1)IW=N_BAND-L+1DO 10 J=1,IWM=J+I-K10 AK(I,J)=AK(I,J)-C*AK(K,M)20 RESULT_N(I)=RESULT_N(

27、I)-C*RESULT_N(K)RESULT_N(N_DOF)=RESULT_N(N_DOF)/AK(N_DOF,1)DO 40 I1=1,N_DOF-1I=N_DOF-I1IF(N_BAND.GT.N_DOF-I-1)JQ=N_DOF-I+1IF(N_BAND.LE.N_DOF-I-1)JQ=N_BANDDO 30 J=2,JQK=J+I-130 RESULT_N(I)=RESULT_N(I)-AK(I,J)*RESULT_N(K)40 RESULT_N(I)=RESULT_N(I)/AK(I,1)WRITE(8,50)50 FORMAT(/12X,'* * * * * RESULT

28、S BY FEM2D * * * * *',/8X,&'-DISPLACEMENT OF NODE-'/5X,'NODE NO',8X,'X-DISP',8X,'Y-DISP')DO 60 I=1,N_NODE60 WRITE(8,70) I,RESULT_N(2*I-1),RESULT_N(2*I)70 FORMAT(8X,I5,7X,2E15.6)RETURNENDcc calculate the stress components of element and nodeSUBROUTINE CAL_STS(N

29、_ELE,N_NODE,N_DOF,PE,PR,IJK_ELE,X,Y,RESULT_N, &STS_ELE,STS_ND)DIMENSION IJK_ELE(500,3),X(N_NODE),Y(N_NODE),DD(3,3),BB(3,6), &SS(3,6),RESULT_N(N_DOF),DISP_E(6)DIMENSION STS_ELE(500,3),STS_ND(500,3)WRITE(8,10)10 FORMAT(/8X,'-STRESSES OF ELEMENT-')CALL CAL_DD(PE,PR,DD)DO 50 IE=1,N_ELECA

30、LL CAL_BB(IE,N_NODE,N_ELE,IJK_ELE,X,Y,AE,BB)DO 20 I=1,3DO 20 J=1,6SS(I,J)=0.0DO 20 K=1,320 SS(I,J)=SS(I,J)+DD(I,K)*BB(K,J)DO 30 I=1,3DO 30 J=1,2IH=2*(I-1)+JIW=2*(IJK_ELE(IE,I)-1)+J30 DISP_E(IH)=RESULT_N(IW)STX=0STY=0TXY=0DO 40 J=1,6STX=STX+SS(1,J)*DISP_E(J)STY=STY+SS(2,J)*DISP_E(J)40 TXY=TXY+SS(3,J)

31、*DISP_E(J)STS_ELE(IE,1)=STXSTS_ELE(IE,2)=STYSTS_ELE(IE,3)=TXY50 WRITE(8,60)IE,STX,STY,TXY60 FORMAT(1X,'ELEMENT NO.=',I5/18X,'STX=',E12.6,5X,'STY=', &E12.6,2X,'TXY=',E12.6)c the following part is to calculate stress components of nodeWRITE(8,55)55 FORMAT(/8X,'-

32、STRESSES OF NODE-')DO 90 I=1,N_NODEA=0.B=0.C=0.II=0DO 70 K=1,N_ELEDO 70 J=1,3IF(IJK_ELE(K,J).EQ.I) THENII=II+1A=A+STS_ELE(K,1)B=B+STS_ELE(K,2)C=C+STS_ELE(K,3)END IF70 CONTINUESTS_ND(I,1)=A/IISTS_ND(I,2)=B/IISTS_ND(I,3)=C/IIWRITE(8,75)I,STS_ND(I,1),STS_ND(I,2),STS_ND(I,3)75 FORMAT(1X,'NODE NO.=',I5/18X,'STX=',E12.6,5X,'STY=', &E12.6,2X,'TXY=',E12.6)90 CONTINUERETURNENDc FEM2D programm end 算例結(jié)果 ：chengxu.for所輸出的數(shù)據(jù)文件DATA.OU數(shù)據(jù)內(nèi)容如下:PLANE STRESS PROBLEM* * * * * RESULTS BY FEM2D * * * * *-DISPLACEMENT