《連續(xù)體平面問題的有限元程序分析-有限元課程設(shè)計(jì)例子》由會(huì)員分享，可在線閱讀，更多相關(guān)《連續(xù)體平面問題的有限元程序分析-有限元課程設(shè)計(jì)例子（16頁珍藏版）》請?jiān)谘b配圖網(wǎng)上搜索。

1、 有 限 元 分 析 程 序 設(shè) 計(jì) 學(xué)校： 院系：建筑工程與力學(xué)學(xué)院 專業(yè)： 姓名： 學(xué)號(hào)： 指導(dǎo)老師： 連續(xù)體平面問題的有限元程序分析 [題目]： 如圖所示的正方形薄板四周受均勻載荷的作用，該結(jié)構(gòu)在邊界上受正向分布壓力，，同時(shí)在沿對角線y軸上受一對集中壓力，載荷為2KN，若取板厚，泊松比。 2kN 2kN 1kN/m [分析過程]： 由于連續(xù)平板的對稱性，只需要取其

2、在第一象限的四分之一部分參加分析，然后人為作出一些輔助線將平板“分割”成若干部分，再為每個(gè)部分選擇分析單元。采用將此模型化分為4個(gè)全等的直角三角型單元。利用其對稱性，四分之一部分的邊界約束，載荷可等效如圖所示。 1kN/m [程序原理及實(shí)現(xiàn)]： 用FORTRAN程序的實(shí)現(xiàn)。由節(jié)點(diǎn)信息文件NODE.IN和單元信息文件ELEMENT.IN，經(jīng)過計(jì)算分析后輸出一個(gè)一般性的文件DATA.OUT。模型基本信息由文件為BASIC.IN生成。 該程序的特點(diǎn)如下： 問題類型：可用于計(jì)算彈性力學(xué)平面問題和平面應(yīng)變問題 單元類型：采用常應(yīng)變?nèi)切螁卧? 位移模

3、式：用用線性位移模式 載荷類型：節(jié)點(diǎn)載荷，非節(jié)點(diǎn)載荷應(yīng)先換算為等效節(jié)點(diǎn)載荷 材料性質(zhì)：彈性體由單一的均勻材料組成 約束方式：為“0”位移固定約束，為保證無剛體位移，彈性體至少應(yīng)有對三個(gè)自由度的獨(dú)立約束 方程求解：針對半帶寬剛度方程的Gauss消元法 輸入文件：由手工生成節(jié)點(diǎn)信息文件NODE.IN，和單元信息文件ELEMENT.IN 結(jié)果文件：輸出一般的結(jié)果文件DATA.OUT 程序的原理如框圖： 開始 輸入數(shù)據(jù)（子程序READ_IN） BASIC.IN（基本信息文件） NODE.IN（節(jié)點(diǎn)信息文件） ELEMENT.IN(單元信息文件) 形成單元?jiǎng)偠染仃嚕ㄗ映绦騀O

4、RM_KE） 以半帶存儲(chǔ)方式形成整體剛度矩陣（BAND_K） 形成節(jié)點(diǎn)載荷向量（子程序FORM_P） 處理邊界條件（子程序DO_BC） 求解方程獲得節(jié)點(diǎn)位移（子程序SOLVE） 計(jì)算單元及節(jié)點(diǎn)應(yīng)力（子程序） 結(jié)束 輸出方件 DATA.OUT （1）主要變量： ID： 問題類型碼，ID＝1時(shí)為平面應(yīng)力問題，ID=2時(shí)為平面應(yīng)變問題 N_NODE： 節(jié)點(diǎn)個(gè)數(shù) N_LOAD： 節(jié)點(diǎn)載荷個(gè)數(shù) N_DOF： 自由度，N_DOF=N_NODE*2（平面問題） N_ELE： 單元個(gè)數(shù) N

5、_BAND： 矩陣半帶寬 N_BC： 有約束的節(jié)點(diǎn)個(gè)數(shù) PE： 彈性模量 PR： 泊松比 PT： 厚度 LJK_ELE(I,3)： 單元節(jié)點(diǎn)編號(hào)數(shù)組，LJK_ELE(I,1)，LJK_ELE(I,2)，LJK_ELE(I,3)分別放單元I的三個(gè)節(jié)點(diǎn)的整體編號(hào) X(N_NODE), Y(N_NODE)：節(jié)點(diǎn)坐標(biāo)數(shù)組，X(I),Y(I)分別存放節(jié)點(diǎn)I的x，y坐標(biāo)值 P_LJK(N_BC,3)： 節(jié)點(diǎn)載荷數(shù)組，P_LJK(I,1)表示第I個(gè)作用有節(jié)點(diǎn)載荷的節(jié)點(diǎn)的編號(hào)，P_LJK(I,2)，P_LJK(I,3)分別為該節(jié)點(diǎn)沿x,y方向的節(jié)點(diǎn)載

6、荷數(shù)值 AK(N_DOF,N_BAND)： 整體剛度矩陣 AKE(6,6)： 單元?jiǎng)偠染仃? BB(3,6)： 位移……應(yīng)變轉(zhuǎn)換矩陣（三節(jié)點(diǎn)單元的幾何矩陣） DD(3,3)： 彈性矩陣 SS(3,6)； 應(yīng)力矩陣 RESULT_N(N_NOF)： 節(jié)點(diǎn)載荷數(shù)組，存放節(jié)點(diǎn)載荷向量，解方程后該矩陣存放節(jié)點(diǎn)位移 DISP_E(6):： 單元的節(jié)點(diǎn)位移向量 STS_ELE(N_ELE,3)： 單元的應(yīng)力分量 STS_ND(N_N

7、ODE,3)： 節(jié)點(diǎn)的應(yīng)力分量 （2）子程序說明： READ_IN： 讀入數(shù)據(jù) BAND_K： 形成半帶寬的整體剛度矩陣 FORM_KE： 計(jì)算單元?jiǎng)偠染仃? FORM_P： 計(jì)算節(jié)點(diǎn)載荷 CAL_AREA：計(jì)算單元面積 DO_BC： 處理邊界條件 CLA_DD： 計(jì)算單元彈性矩陣 SOLVE： 計(jì)算節(jié)點(diǎn)位移 CLA_BB： 計(jì)算單元位移……應(yīng)變關(guān)系矩陣 CAL_STS：計(jì)算單元和節(jié)點(diǎn)應(yīng)力 （3）文件管理： 源程序文件：

8、 chengxu.for 程序需讀入的數(shù)據(jù)文件： BASIC.IN，NODE.IN，ELEMENT.IN（需要手工生成） 程序輸出的數(shù)據(jù)文件： DATA.OUT （4）數(shù)據(jù)文件格式： 需讀入的模型 基本信息文件BASIC.IN的格式如下表 欄目 格式說明 實(shí)際需輸入的數(shù)據(jù) 基本模型數(shù)據(jù) 第1行，每兩個(gè)數(shù)之間用“，”號(hào)隔開 問題類型，單元個(gè)數(shù)，節(jié)點(diǎn)個(gè)數(shù)，有約束的節(jié)點(diǎn)數(shù)，有載何的節(jié)點(diǎn)數(shù) 材料性質(zhì) 第2行，每兩個(gè)數(shù)之間用“，”號(hào)隔開 彈性模量，泊松比，單元厚度 節(jié)點(diǎn)約束信息 在材料性質(zhì)輸入行之后另起行，每兩個(gè)數(shù)之間用“，”號(hào)隔開 LJK_U(N_BC,3) 位

9、移約束的節(jié)點(diǎn)編號(hào)，該節(jié)點(diǎn)x方向約束代碼，該節(jié)點(diǎn)y方向代碼， 節(jié)點(diǎn)荷載信息 在節(jié)點(diǎn)約束信息輸入行之后另起行，每兩個(gè)數(shù)之間用“，”號(hào)隔開 P_IJK(N_LOAD,3) 載荷作用的節(jié)點(diǎn)編號(hào)，該節(jié)點(diǎn)x主向載荷，該節(jié)點(diǎn)y方向載荷，…… 需讀入的節(jié)點(diǎn)信息文件NODE.IN的格式如下表 欄目 格式說明 實(shí)際需輸入的數(shù)據(jù) 節(jié)點(diǎn)信息 每行為一個(gè)節(jié)點(diǎn)的信息（每行三個(gè)數(shù)，每兩個(gè)數(shù)之間用空格或“，”分開） ND_ANSYS(N_NIDE) 節(jié)點(diǎn)號(hào)，該節(jié)點(diǎn)的x坐標(biāo)，該節(jié)點(diǎn)y方向坐標(biāo) 需讀入的單元信息文件ELEMENT.IN的格式如下表 欄目 格式說明 實(shí)際需輸入的數(shù)據(jù) 單元信

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

11、x方向位移 y方向位移 單元應(yīng)力的三個(gè)分量 IE STE_ELE(IE,1) STE_ELE(IE,2) STE_ELE(IE,3) 單元號(hào) x方向應(yīng)力 y方向應(yīng)力 剪切應(yīng)力 節(jié)點(diǎn)應(yīng)力的三個(gè)分量 I STS-ND(I,1) STS-ND(I,2) STS-ND(I,3) 節(jié)點(diǎn)號(hào) x方向應(yīng)力 y方向應(yīng)力 剪切應(yīng)力 [算例原始數(shù)據(jù)和程序分析]： （1）模型基本信息文件BASIC.IN的數(shù)據(jù)為 1,4,6,5,3 1.,0.,1. 1,1,0,2,1,0,4,1,1,5,0,1,6,0,1 1,-

12、0.5,-1.5,3.,-1.,-1,6,-0.5,-0.5 （2）手工準(zhǔn)備的節(jié)點(diǎn)信息文件NODE.IN的數(shù)據(jù)為 1 0.0 2.0 2 0.0 1.0 3 1.0 1.0 4 0. 0. 5 1.0 0. 6 2.0 0. （3）手工準(zhǔn)備的單元信息文件ELEMENT.IN的數(shù)據(jù)為 1 2 3 3 0 0 0 0 1 1 1 1 0

13、1 2 4 5 5 0 0 0 0 1 1 1 1 0 2 5 3 2 2 0 0 0 0 1 1 1 1 0 3 3 5 6 6 0 0 0 0 1 1 1 1 0 4 （4）源程序文件chengxu.for為： PROGRAM FEM2D DIMENSION IJK_ELE(500,3),X(

14、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='BASIC.IN') OPEN(5,FILE='NODE.IN') OPEN(6,FILE='ELEMENT.IN') OPEN(8,FILE='DATA.OUT') OPEN(9,FILE='FOR_POST.DAT') READ(4,*)ID,N_ELE,N_NODE,N_BC,N_LOAD IF(ID.E

15、Q.1)WRITE(8,20) IF(ID.EQ.2)WRITE(8,25) 20 FORMAT(/5X,'=========PLANE STRESS PROBLEM========') 25 FORMAT(/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,IE,N_NOD

16、E, & 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, &

17、 STS_ELE,STS_ND) c to putout a data file WRITE(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_E

18、LE(I,2)), &REAL(IJK_ELE(I,3)),REAL(IJK_ELE(I,3)), &STS_ELE(I,1),STS_ELE(I,2),STS_ELE(I,3),I=1, N_ELE) 72 FORMAT(7f9.4) c CLOSE(4) CLOSE(5) CLOSE(6) CLOSE(8) CLOSE(9) END c c to get the original data in order to model the problem SUBROUTINE READ_IN(ID

19、,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_IJK(N_LOAD,3),NE_ANSYS(N_ELE,14) REAL ND_ANSYS(N_NODE,3) READ(4,*)PE,PR,PT READ(4,*)((IJK_U(I,J),J=1,3),I=1,N_BC) READ(4,*)((P_IJK(I,J),J=1

20、,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_NODE X(I)=ND_ANSYS(I,2) Y(I)=ND_ANSYS(I,3) 10 CONTINUE DO 11 I=1,N_ELE DO 11 J=1,3 IJK_ELE(I,J)=NE_ANSYS(I,J) 11 CONTINUE N_BAND=0 DO 20 IE=1,N_E

21、LE DO 20 I=1,3 DO 20 J=1,3 IW=IABS(IJK_ELE(IE,I)-IJK_ELE(IE,J)) IF(N_BAND.LT.IW)N_BAND=IW 20 CONTINUE N_BAND=(N_BAND+1)*2 IF(ID.EQ.1) THEN ELSE PE=PE/(1.0-PR*PR) PR=PR/(1.0-PR) END IF RETURN END c C to form the stiffness matrix of element SUBRO

22、UTINE FORM_KE(IE,N_NODE,N_ELE,IJK_ELE,X,Y,PE,PR,PT,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,PR,DD) CALL CAL_BB(IE,N_NODE,N_ELE,IJK_ELE,X,Y,AE,BB) DO 10 I=1,3 DO 10 J=1,6 SS(I,J)=0.0 DO 10 K=1,3 10 SS(I,J)=SS(I

23、,J)+DD(I,K)*BB(K,J) DO 20 I=1,6 DO 20 J=1,6 AKE(I,J)=0.0 DO 20 K=1,3 20 AKE(I,J)=AKE(I,J)+SS(K,I)*BB(K,J)*AE*PT RETURN END c c to form banded global stiffness matrix SUBROUTINE BAND_K(N_DOF,N_BAND,N_ELE,IE,N_NODE,IJK_ELE,X,Y,PE, & PR,PT,AK) DIMENSI

24、ON IJK_ELE(500,3),X(N_NODE),Y(N_NODE),AKE(6,6),AK(500,100) N_DOF=2*N_NODE DO 40 I=1,N_DOF DO 40 J=1,N_BAND 40 AK(I,J)=0 DO 50 IE=1,N_ELE CALL FORM_KE(IE,N_NODE,N_ELE,IJK_ELE,X,Y,PE,PR,PT,AKE) DO 50 I=1,3 DO 50 II=1,2 IH=2*(I-1)+II IDH=2*(IJK_ELE(IE,I)

25、-1)+II DO 50 J=1,3 DO 50 JJ=1,2 IL=2*(J-1)+JJ IZL=2*(IJK_ELE(IE,J)-1)+JJ IDL=IZL-IDH+1 IF(IDL.LE.0) THEN ELSE AK(IDH,IDL)=AK(IDH,IDL)+AKE(IH,IL) END IF 50 CONTINUE RETURN END c c to calculate the area of element SU

26、BROUTINE 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.0 RETURN END c c to calculate the elastic matrix of element

27、 SUBROUTINE CAL_DD(PE,PR,DD) DIMENSION DD(3,3) DO 10 I=1,3 DO 10 J=1,3 10 DD(I,J)=0.0 DD(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) RETURN END c c to calculate the strain-displacement matrix of e

28、lement SUBROUTINE CAL_BB(IE,N_NODE,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,3 DO 10 JJ=1,3 10 BB(II,JJ)=0.0 BB(1,1)=Y(J)-Y(K) BB(1,3)=Y(K)-Y(I) BB(1,5)=Y(I)-Y(J)

29、 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) 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,3 DO 20 J1=1,6 20 BB(I1,J1)=BB(I1,J1)/(2.0*AE) RETURN END c

30、 c to form the global load matrix SUBROUTINE 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),P_IJK(N_LOAD,3), & RESULT_N(N_DOF) DO 10 I=1,N_DOF 10 RESULT_N(I)=0.0 DO 20 I=1,N_LOAD II=P_IJK(I,1) RESULT_N(2*II

31、-1)=P_IJK(I,2) 20 RESULT_N(2*II)=P_IJK(I,3) RETURN END c c to deal with BC(u) (here only for fixed displacement) using "1-0" method SUBROUTINE DO_BC(N_BC,N_BAND,N_DOF,IJK_U,AK,RESULT_N) DIMENSION RESULT_N(N_DOF),IJK_U(N_BC,3),AK(500,100) DO 30 I=1,N_BC IR=IJK_U(I,1)

32、DO 30 J=2,3 IF(IJK_U(I,J).EQ.0)THEN ELSE II=2*IR+J-3 AK(II,1)=1.0 RESULT_N(II)=0.0 DO 10 JJ=2,N_BAND 10 AK(II,JJ)=0.0 DO 20 JJ=2,II 20 AK(II-JJ+1,JJ)=0.0 END IF 30 CONTINUE RETURN

33、 END c c to solve the banded FEM equation by GAUSS elimination SUBROUTINE SOLVE(N_NODE,N_DOF,N_BAND,AK,RESULT_N) DIMENSION RESULT_N(N_DOF),AK(500,100) DO 20 K=1,N_DOF-1 IF(N_DOF.GT.K+N_BAND-1)IM=K+N_BAND-1 IF(N_DOF.LE.K+N_BAND-1)IM=N_DOF DO 20 I=K+1,IM L=I-K+1 C=AK(K,L)/AK(

34、K,1) IW=N_BAND-L+1 DO 10 J=1,IW M=J+I-K 10 AK(I,J)=AK(I,J)-C*AK(K,M) 20 RESULT_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-1 I=N_DOF-I1 IF(N_BAND.GT.N_DOF-I-1)JQ=N_DOF-I+1 IF(N_BAND.LE.N_DOF-I-1)JQ=N_BAND

35、 DO 30 J=2,JQ K=J+I-1 30 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,'* * * * * RESULTS BY FEM2D * * * * *',//8X, &'--DISPLACEMENT OF NODE--'//5X,'NODE NO',8X,'X-DISP',8X,'Y-DISP') DO 60 I=1,N_NODE 60

36、 WRITE(8,70) I,RESULT_N(2*I-1),RESULT_N(2*I) 70 FORMAT(8X,I5,7X,2E15.6) RETURN END c c calculate the stress components of element and node SUBROUTINE CAL_STS(N_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)

37、,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_ELE CALL CAL_BB(IE,N_NODE,N_ELE,IJK_ELE,X,Y,AE,BB) DO 20 I=1,3 DO 20 J

38、=1,6 SS(I,J)=0.0 DO 20 K=1,3 20 SS(I,J)=SS(I,J)+DD(I,K)*BB(K,J) DO 30 I=1,3 DO 30 J=1,2 IH=2*(I-1)+J IW=2*(IJK_ELE(IE,I)-1)+J 30 DISP_E(IH)=RESULT_N(IW) STX=0 STY=0 TXY=0 DO 40 J=1,6 STX=STX+SS(1,J)*DISP_E(J) STY=STY

39、+SS(2,J)*DISP_E(J) 40 TXY=TXY+SS(3,J)*DISP_E(J) STS_ELE(IE,1)=STX STS_ELE(IE,2)=STY STS_ELE(IE,3)=TXY 50 WRITE(8,60)IE,STX,STY,TXY 60 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 compon

40、ents of node WRITE(8,55) 55 FORMAT(//8X,'--STRESSES OF NODE--') DO 90 I=1,N_NODE A=0. B=0. C=0. II=0 DO 70 K=1,N_ELE DO 70 J=1,3 IF(IJK_ELE(K,J).EQ.I) THEN II=II+1 A=A+STS_ELE(K,1) B=B+STS_ELE(K,2) C=C+STS_ELE(K,3) END IF 70 CONTINUE S

41、TS_ND(I,1)=A/II STS_ND(I,2)=B/II STS_ND(I,3)=C/II WRITE(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 CONTINUE RETURN END c FEM2D programm end [算例結(jié)果]： chengxu.for所輸出的數(shù)據(jù)文件DATA.OUT數(shù)據(jù)內(nèi)容如下：

42、 =========PLANE STRESS PROBLEM======== * * * * * RESULTS BY FEM2D * * * * * --DISPLACEMENT OF NODE-- NODE NO X-DISP Y-DISP 1 .000000E+00 -.525275E+01 2 .000000E+00 -.225275E+01 3

43、 -.108791E+01 -.137363E+01 4 .000000E+00 .000000E+00 5 -.824176E+00 .000000E+00 6 -.182418E+01 .000000E+00 --STRESSES OF ELEMENT-- ELEMENT NO.= 1 STX=-.108791E+01 STY=-.300000

44、E+01 TXY= .439560E+00 ELEMENT NO.= 2 STX=-.824176E+00 STY=-.225275E+01 TXY= .000000E+00 ELEMENT NO.= 3 STX=-.108791E+01 STY=-.137363E+01 TXY= .307692E+00 ELEMENT NO.= 4 STX=-.100000E+01 STY=-.137363E+01 TXY=-.131

45、868E+00 --STRESSES OF NODE-- NODE NO.= 1 STX=-.108791E+01 STY=-.300000E+01 TXY= .439560E+00 NODE NO.= 2 STX=-.100000E+01 STY=-.220879E+01 TXY= .249084E+00 NODE NO.= 3 STX=-.105861E+01 STY=-.191575E+01

46、 TXY= .205128E+00 NODE NO.= 4 STX=-.824176E+00 STY=-.225275E+01 TXY= .000000E+00 NODE NO.= 5 STX=-.970696E+00 STY=-.166667E+01 TXY= .586081E-01 NODE NO.= 6 STX=-.100000E+01 STY=-.137363E+01 TXY=-.131868E+00 [結(jié)論與體會(huì)]: 通過本次的課程設(shè)計(jì)，我對有限元的概念有了更加深刻的理解，同時(shí)也彌補(bǔ)了平時(shí)學(xué)習(xí)是疏忽的地方，充實(shí)了有限元知識(shí)。通過編程及對chengxu.for的調(diào)試進(jìn)一步鞏固了FORTRAN語言知識(shí)，也鍛煉了自我，使自己在學(xué)習(xí)工作中更具耐心、細(xì)心。 [參考數(shù)目]： 《FORTRAN語言》清華大學(xué)出版社 譚浩強(qiáng)，田淑清 《彈性力學(xué)》 高等教育出版社 徐芝綸 《有限元分析及其應(yīng)用》清華大學(xué) 曾攀