有限元課程設(shè)計例子分析

1、有限兀分析程序設(shè)計學校：燕山大學 院系：建筑工程與力學學院 專業(yè)：01級工程力學 姓名：張任良 學號：010107030025 指導老師：杜國君完成時間2004年12月31日連續(xù)體平面問題的有限元程序分析題目:如圖所示的正方形薄板四周受均勻載荷的作用，該結(jié)構(gòu)在邊界 上受正向分布壓力，Pik%，同時在沿對角線y軸上受一對集中壓 力，載荷為2KN若取板厚f，泊松比v=0。分析過程:由于連續(xù)平板的對稱性，只需要取其在第一象限的四分之一部分參加分析，然后人為作出一些輔助線將平板“分割”成若干部分，再 為每個部分選擇分析單元。采用將此模型化分為4個全等的直角三角 型單元。利用其對稱性，四分之一部分的邊界

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

3、點信息文件NODE.IN和單元信息文件ELEMENTN結(jié)果文件：輸出一般的結(jié)果文件 DATA.OUT程序的原理如框圖:開始輸入數(shù)據(jù)（子程序 READ_IN）BASIC.IN （基本信息文件）NODE.IN （節(jié)點信息文件）ELEMENT.IN（單元信息文件）1形成單元剛度矩陣（子程序 FORM_KE）ID：N_NODE： N_LOAD ： N_DOF： N_ELE： N_BAND ： N_BC： PE： PR： PT：整體剛度矩陣單元剛度矩陣位移應變轉(zhuǎn)換矩陣(三節(jié)點單元的幾何矩陣) 彈性矩陣應力矩陣節(jié)點載荷數(shù)組， 存放節(jié)點載荷向量， 解方程后該矩單元的節(jié)點位移向量 單元的應力分量 節(jié)點的應力分

4、量讀入數(shù)據(jù) 計算單元剛度矩陣 ：計算單元面積 計算單元彈性矩陣 計算單元位移 應變關(guān)系矩陣READ_IN：FORM_KE ：CAL_AREACLA_DD ：CLA_BB ：CAL_STS :計算單元和節(jié)點應力BAND_K ：形成半帶寬的整體剛度矩陣FORM_P：DO_BC：SOLVE：計算節(jié)點載荷 處理邊界條件 計算節(jié)點位移1)主要變量：問題類型碼，ID = 1時為平面應力問題，ID=2時為平面應變問題 節(jié)點個數(shù)節(jié)點載荷個數(shù)自由度， N_DOF=N_NODE*2 (平面問題) 單元個數(shù)矩陣半帶寬 有約束的節(jié)點個數(shù) 彈性模量 泊松比 厚度LJK_ELE(I,3) ： 單元節(jié)點編號數(shù)組， LJK_

5、ELE(I,1) ， LJK_ELE(I,2) ， LJK_ELE(I,3) 分別放單元 I 的三個節(jié)點的整體編號X(N_NODE), Y(N_NODE) ：節(jié)點坐標數(shù)組， X(I),Y(I) 分別存放節(jié)點 I 的 x， y 坐標值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)：DD(3,3)：SS(3,6)；RESULT_N(N_NOF) ： 陣存放節(jié)點位移DISP_E(6

6、):： STS_ELE(N_ELE,3)：STS_ND(N_NODE,3)：2)子程序說明：3)文件管理：源程序文件：chengxu.for程序需讀入的數(shù)據(jù)文件：BASIC.IN ，NODE.IN ， ELEMENT.IN （需要手工生成）程序輸出的數(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ì)輸入行之后另起行

7、，每 兩個數(shù)之間用“，”號隔開LJK U（N BC,3）位移約束的節(jié)點編號，該 節(jié)點x方向約束代碼，該 節(jié)點y方向代碼，節(jié)點荷載信息在節(jié)點約束信息輸入行之后另起 行，每兩個數(shù)之間用“，”號隔開P IJK（N LOAD,3）載荷作用的節(jié)點編號，該 節(jié)點x主向載何，該節(jié)點 y方向載何，需讀入的節(jié)點信息文件NODE.IN的格式如下表欄目格式說明實際需輸入的數(shù)據(jù)節(jié)點信息每行為一個節(jié)點的信息 （每行三個數(shù)，每兩個數(shù) 之間用空格或“，”分開）ND_ANSYS（N_NIDE）節(jié)點號，該節(jié)點的x坐標， 該節(jié)點y方向坐標需讀入的單元信息文件ELEMENT.IN的格式如下表欄目格式說明實際需輸入的數(shù)據(jù)單元信息每行

8、為一個單元的信息（每行有14個整型數(shù)， 前4個為單兀節(jié)點編號， 對于3節(jié)點編號，第4個 節(jié)點編號與第3個節(jié)點編 號相同，后10個數(shù)無用， 可輸入“ 0 ”，每兩個整 型數(shù)之間用至少一個空 格分開）NE_ANSYS（N_ELE,14） 單元的節(jié)點號1 （空格） 單元的節(jié)點號2 （空格） 單元的節(jié)點號3 （空格） 單元的節(jié)點號4 （空格） 0 （空格）0 （空格）0 （空 格）0 （空格）0 （空格） 0 （空格）0 （空格）0 （空 格）0 （空格）0輸出結(jié)果文件DATA.OUT格式如下表欄目實際輸出的數(shù)據(jù)節(jié)點位移IRESULT_N(2T_ 1)RESULT_N(2*I)節(jié)點號x方向位移y方向位

9、移單元應力的三個分量IESTE_ELE(IE,1)STE_ELE(IE,2)STE_ELE(IE,3)單兀號 x方向應力y方向應力剪切應力節(jié)點應力的三個分量I STS-ND(I,1)STS-ND(I,2)STS-ND(I,3)節(jié)點號 x方向應力y方向應力剪切應力算例原始數(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)手工準備的節(jié)點信息文件NODE.IN的數(shù)據(jù)為10.02.020.01.031.01.040.0.51.00.

10、62.00.(3)手工準備的單元信息文件ELEMENTN的數(shù)據(jù)為12532435352635260000000000000000111111111111111100001234(4)源程序文件che ngxu.for 為:P ROGRAM 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)OP EN(4,FILE=BASIC.IN)OP EN(5,FILE=NODE.IN)OP E

11、N(6,FILE=ELEMENTN)OP EN(8,FILE=DATA.OUT)OPEN(9,FILE=FOR_POST.DA T)READ(4,*)ID,N_ELE,N_NODE,N_BC,N_LOADIF(ID.EQ.1)WRITE(8,20)IF(ID.EQ.2)WRITE(8,25)20FORMAT(/5X,=PLANE STRESS P ROBLEM=)25FORMAT(/5X,=PLANE STRAIN P ROBLEM=)CALL READ_IN(ID,N_ELE,N_NODE,N_BC,N_BAND,N_LOAD ,PE,PR,PT,&IJK_ELE,X,Y ,IJK_U,P_

12、IJK)CALL BAND_K(N_DOF,N_BAND,N_ELE,IE,N_NODE,&IJK_ELE,X,Y, PE,PR,P T,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_DO F,PE,P R,IJK_ELE,X,Y,RESULT_N,&STS_ELE,

13、STS_ND)c to pu tout a data fileWRITE(9,70)REAL(N_NODE),REAL(N_ELE)707172FORMAT(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)FORMA T(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)

14、,STS_ELE(I,2),STS_ELE(I,3),I=1, N_ELE)FORMAT(7f9.4)CLOSE (4)CLOSE(5)CLOSE(6)CLOSE(8)CLOSE(9)ENDcc to get the orig inal data in order to model the p roblemSUBROUTINE 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(5OO,3),X(N_NODE),Y(N_NODE),IJK_U(N_BC

15、,3), & P_IJK(N_LOAD,3),NE_ANSYS(N_ELE,14)REAL ND_ANSYS(N_NODE,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 CON

16、TINUEDO 11 I=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.I) THENELSEPE=P E/(1.0-PR*PR)PR=P R/(1.0-PR)END IFRETURNENDcC to form the stiffness mat

17、rix of eleme ntSUBROUTINE FORM_KE(IE,N_NODE,N_ELE,IJK_ELE,X,Y,PE,PR,P T,AKE)DIMENSION IJK_ELE(500,3),X(N_NODE),Y(N_NODE),BB(3,6),DD(3,3),&AKE(6,6), SS(6,6)CALL CAL_DD (PE,P R,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

18、,J)DO 20 I=1,6DO 20 J=1,6AKE(I,J)=0.0DO 20 K=1,320AKE(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,P T,AK)DIMENSIONIJK_ELE(500,3),X(N_NODE),Y(N_NODE),AKE(6,6),AK(500,100)N DOF=2*N NODEDO 40 I

19、=1,N_DOFDO 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,P T,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.O) THENELSEAK(IDH,IDL)=AK(IDH,IDL)+AKE(IH,IL)E

20、ND IF50 CONTINUERETURNENDcc to calculate the area of eleme ntSUBROUTINE CAL_AREA(IE,N_NODE,IJK_ELE,X,Y,AE)DIMENSION IJK_ELE(5OO,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

21、 elastic matrix of element SUBROUTINE 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)RETURNENDc to calculate the strain-displacement matrix of elementSUBROUTINE CAL_BB(IE,N_NOD

22、E,N_ELE,IJK_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

23、(1,1) BB(3,3)=BB(2,4) 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)RETURNENDc 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(

24、500,3),X(N_NODE),Y(N_NODE),P_IJK(N_LOAD,3),& RESULT_N(N_DOF)DO 10 I=1,N_DOF10 RESULT_N(I)=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)RETURNENDcc to deal with BC(u) (here only for fixed dis placeme nt) using 1-0 method SUBROUTINE DO_BC (N_BC,N_BAND,N_DOF,IJ

25、K_U,AK,RESULT_N) 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.0)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 IF30CONTINUERETURNEND cc to solve the ban ded FEM equati o

26、n by GAUSS elim in ati onSUBROUTINE SOLVE(N_NODE,N_DOF,N_BAND,AK,RESULT_N)DIMENSION RESULT_N(N_DOF),AK(5OO,1OO)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 RESU

27、LT_N(I)=RESULT_N(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+1 IF(N_BAND.LE.N_DOF-I-1)JQ=N_BAND DO 30 J=2,JQK=J+I-1RESULT_N(I)=RESULT_N(I)-AK(I,J)*RESULT_N(K) RESULT_N(I)=RESULT_N(I)/AK(I,1)WRITE(8,50)FORMAT(/12X,* * * * *

28、 RESULTS BY FEM2D * * * * *,/8X,&-DISPLACEMENT OF NODE-/5X,NODE NO,8X,X-DISP,8X,Y-DISP) DO 60 I=1,N_NODE60 70304050WRITE(8,70) I,RESULT_N(2*I-1),RESULT_N(2*I)FORMAT(8X,I5,7X,2E15.6)RETURNENDcc calculate the stress components of element and node SUBROUTINECAL_STS(N_ELE,N_NODE,N_DOF,PE,PR,IJK_ELE,X,Y,

29、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_ELECALL CAL_BB(IE,N_NODE,N_ELE,IJK_ELE,X,Y,AE,BB)DO 20 I=1

30、,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)*DISP_E(J)STS_ELE(IE,1)=STXSTS_ELE(IE,2)=STYSTS_ELE(

31、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,-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