有限元大作业matlab---课程设计例子.doc

1、 有 限 元 大 作 业 程 序 设 计 学校：天津大学院系：建筑工程与力学学院专业：01 级工程力学姓名：刘秀学号：指导老师：连续体平面问题的有限元程序分析题目 ： 如图所示的正方形薄板四周受均匀载荷的作用，该结构在边界上受正向分布压力， ，同时在沿对角线 y 轴上受一对集中mkNp1压力，载荷为 2KN，若取板厚 ，泊松比 。t0v2kN2kN1kN/m分析过程 ：由于连续平板的对称性，只需要取其在第一象限的四分之一部分参加分析，然后人为作出一些辅助线将平板“分割”成若干部分，再为每个部分选择分析单元。采用将此模型化分为 4 个全等的直角三角型单元。利用其对称性，四分之一部分的边界约束，载

2、荷可等效如图所示。 程序原理及实现 ：用 FORTRAN 程序的实现。由节点信息文件 NODE.IN 和单元信息文件 ELEMENT.IN，经过计算分析后输出一个一般性的文件DATA.OUT。模型基本信息由文件为 BASIC.IN 生成。该程序的特点如下：问题类型：可用于计算弹性力学平面问题和平面应变问题单元类型：采用常应变三角形单元位移模式：用用线性位移模式载荷类型：节点载荷，非节点载荷应先换算为等效节点载荷 材料性质：弹性体由单一的均匀材料组成约束方式：为“0”位移固定约束，为保证无刚体位移，弹性体至少应有对三个自由度的独立约束方程求解：针对半带宽刚度方程的 Gauss 消元法1kN/m输

3、入文件：由手工生成节点信息文件 NODE.IN，和单元信息文件ELEMENT.IN结果文件：输出一般的结果文件 DATA.OUT程序的原理如框图：开始输入数据（子程序READ_IN）BASIC.IN（基本信息文件）NODE.IN（节点信息文件）ELEMENT.IN(单元信息文件)形成单元刚度矩阵（子程序FORM_KE）以半带存储方式形成整体刚度矩阵（ BAND_K）形成节点载荷向量（子程序 FORM_P）处理边界条件（子程序 DO_BC）求解方程获得节点位移（子程序SOLVE）计算单元及节点应力（子程序）结束输出方件DATA.OUT（1）主要变量：ID： 问题类型码， ID1 时为平面应力问题

4、， ID=2 时为平面应变问题N_NODE： 节点个数N_LOAD： 节点载荷个数N_DOF： 自由度，N_DOF=N_NODE*2 （平面问题）N_ELE： 单元个数N_BAND： 矩阵半带宽N_BC： 有约束的节点个数PE： 弹性模量PR： 泊松比PT： 厚度LJK_ELE(I,3)： 单元节点编号数组，LJK_ELE(I,1)，LJK_ELE(I,2)，LJK_ELE(I,3)分别放单元 I 的三个节点的整体编号X(N_NODE), Y(N_NODE)：节点坐标数组，X(I),Y(I)分别存放节点 I 的 x，y坐标值P_LJK(N_BC,3)： 节点载荷数组，P_LJK(I,1)表示第

5、 I 个作用有节点载荷的节点的编号，P_LJK(I,2)，P_LJK(I,3)分别为该节点沿 x,y 方向的节点载荷数值AK(N_DOF,N_BAND)： 整体刚度矩阵AKE(6,6)： 单元刚度矩阵BB(3,6)： 位移应变转换矩阵（三节点单元的几何矩阵）DD(3,3)： 弹性矩阵SS(3,6)； 应力矩阵RESULT_N(N_NOF)： 节点载荷数组，存放节点载荷向量，解方程后该矩阵存放节点位移DISP_E(6):： 单元的节点位移向量STS_ELE(N_ELE,3)： 单元的应力分量STS_ND(N_NODE,3)： 节点的应力分量（2）子程序说明：READ_IN： 读入数据 BAND_

6、K： 形成半带宽的整体刚度矩阵FORM_KE： 计算单元刚度矩阵 FORM_P： 计算节点载荷CAL_AREA：计算单元面积 DO_BC： 处理边界条件CLA_DD： 计算单元弹性矩阵 SOLVE： 计算节点位移CLA_BB： 计算单元位移应变关系矩阵 CAL_STS：计算单元和节点应力（3）文件管理：源程序文件：chengxu.for程序需读入的数据文件：BASIC.IN，NODE.IN，ELEMENT.IN（需要手工生成）程序输出的数据文件：DATA.OUT（4）数据文件格式：需读入的模型 基本信息文件 BASIC.IN 的格式如下表需读入的节点信息文件 NODE.IN 的格式如下表栏目

7、格式说明 实际需输入的数据节点信息每行为一个节点的信息（每行三个数，每两个数之间用空格或“， ”分开）ND_ANSYS(N_NIDE)节点号，该节点的 x 坐标，该节点 y 方向坐标需读入的单元信息文件 ELEMENT.IN 的格式如下表栏目 格式说明 实际需输入的数据栏目 格式说明 实际需输入的数据基本模型数据 第 1 行，每两个数之间用“， ”号隔开问题类型，单元个数，节点个数，有约束的节点数，有载何的节点数材料性质 第 2 行，每两个数之间用“， ”号隔开弹性模量，泊松比，单元厚度节点约束信息 在材料性质输入行之后另起行， 每两个数之间用“， ”号隔开LJK_U(N_BC,3)位移约束的

8、节点编号，该节点 x 方向约束代码，该节点 y 方向代码，节点荷载信息 在节点约束信息输入行之后另起行，每两个数之间用“， ”号隔开P_IJK(N_LOAD,3)载荷作用的节点编号，该节点 x 主向载荷，该节点 y 方向载荷，单元信息每行为一个单元的信息（每行有 14 个整型数，前 4 个为单元节点编号，对于 3 节点编号，第 4个节点编号与第 3 个节点编号相同，后 10 个数无用，可输入“0” ，每两 个整型数之间用至少一个空格分开）NE_ANSYS(N_ELE,14)单元的节点号 1（空格）单元的节点号 2（空格）单元的节点号 3（空格）单元的节点号 4（空格）0（空格）0（空格）0（空

9、格）0（空格）0（空格）0（空格）0（空格）0（空格）0（空格）0输出结果文件 DATA.OUT 格式如下表栏目 实际输出的数据节点位移 I RESULT_N(2*I_ 1) RESULT_N(2*I)节点号 x 方向位移 y 方向位移单元应力的三个分量IE STE_ELE(IE,1) STE_ELE(IE,2) STE_ELE(IE,3)单元号 x 方向应力 y 方向应力 剪切应力节点应力的三个分量 I STS-ND(I,1) STS-ND(I,2) STS-ND(I,3)节点号 x 方向应力 y 方向应力 剪切应力算例原始数据和程序分析 ：（1）模型基本信息文件 BASIC.IN 的数据为

10、1,4,6,5,31.,0.,1.1,1,0,2,1,0,4,1,1,5,0,1,6,0,11,-0.5,-1.5,3.,-1.,-1,6,-0.5,-0.5（2）手工准备的节点信息文件 NODE.IN 的数据为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）手工准备的单元信息文件 ELEMENT.IN 的数据为1 2 3 3 0 0 0 0 1 1 1 1 0 12 4 5 5 0 0 0 0 1 1 1 1 0 25 3 2 2 0 0 0 0 1 1 1 1 0 33 5 6 6 0 0 0 0 1 1 1 1 0 4（4）

11、源程序文件 chengxu.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=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

12、,N_LOADIF(ID.EQ.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_NODE,& IJK_ELE,X,Y,PE,PR,PT,AK) CALL FORM_P(N_ELE,N_NODE

13、,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 file WRITE(9,70)REAL(N_NODE),REAL(N_ELE)70 FORMAT(2f9.4)WRITE(9,71)(X(I)

14、,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_ELE(I,3),I=1, N_ELE)72 FORMAT(7f9.4)cCLOSE(4)CLOSE(5)CLOSE(6)CLOSE(8)CLOSE(9)ENDc c

15、 to get the original 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_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),

16、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=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,3

17、DO 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) THEN ELSEPE=PE/(1.0-PR*PR)PR=PR/(1.0-PR)END IFRETURNENDcC to form the stiffness matrix of element SUBROUTINE FORM_KE(IE,N_NODE,N_ELE,IJK_ELE,X,Y,PE,PR,PT,AKE)DIMENSION IJK_ELE(500,3),X(N_

18、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,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*PTRETURN ENDcc to f

19、orm banded 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_DOFDO 40 J=1,N_BAND40 AK(I,J)=0DO 50 IE=1,N_ELECALL FORM_KE(IE,N_NODE,N_ELE,IJK_ELE,X,Y,PE,PR,PT,AKE)D

20、O 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) THEN ELSEAK(IDH,IDL)=AK(IDH,IDL)+AKE(IH,IL) END IF50 CONTINUERETURNENDcc to calculate the area of element SUBROUTINE CAL_AREA(IE,N_NODE,IJK_ELE,X,Y,AE)

21、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.0RETURNENDcc to calculate the elastic 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.0D

22、D(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 ENDcc to calculate the strain-displacement matrix of element 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,

23、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)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,

24、3DO 20 J1=1,620 BB(I1,J1)=BB(I1,J1)/(2.0*AE)RETURN ENDcc 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_DOF10 RESULT_N(I)=0.0 DO 20 I=1,N_LOADII=P_IJK(I,

25、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 displacement) using “1-0“ methodSUBROUTINE 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_BCIR=IJK_U(I,1)DO 30 J=2,3IF(IJK_U(I,J).EQ.

26、0)THENELSEII=2*IR+J-3AK(II,1)=1.0RESULT_N(II)=0.0DO 10 JJ=2,N_BAND10 AK(II,JJ)=0.0DO 20 JJ=2,II20 AK(II-JJ+1,JJ)=0.0END IF 30 CONTINUERETURNENDc c to solve the banded FEM equation by GAUSS eliminationSUBROUTINE SOLVE(N_NODE,N_DOF,N_BAND,AK,RESULT_N)DIMENSION RESULT_N(N_DOF),AK(500,100)DO 20 K=1,N_DO

27、F-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(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-

28、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,* * * * * RESULTS 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,

29、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_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)DIMENS

30、ION 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,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)+J3

31、0 DISP_E(IH)=RESULT_N(IW)STX=0 STY=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(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 pa

32、rt 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,

33、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 CONTINUERETURN ENDc FEM2D programm end算例结果 ：chengxu.for 所输出的数据文件 DATA.OUT 数据内容如下：=PLANE STRESS PROBLEM=* * * * * RESULTS BY FEM2D * * * * *-DISPLACEMENT OF NODE-NODE NO X-DISP Y-DISP1 .00

34、0000E+00 -.525275E+012 .000000E+00 -.225275E+013 -.108791E+01 -.137363E+014 .000000E+00 .000000E+005 -.824176E+00 .000000E+006 -.182418E+01 .000000E+00-STRESSES OF ELEMENT-ELEMENT NO.= 1STX=-.108791E+01 STY=-.300000E+01 TXY= .439560E+00ELEMENT NO.= 2STX=-.824176E+00 STY=-.225275E+01 TXY= .000000E+00

35、ELEMENT NO.= 3STX=-.108791E+01 STY=-.137363E+01 TXY= .307692E+00ELEMENT NO.= 4STX=-.100000E+01 STY=-.137363E+01 TXY=-.131868E+00-STRESSES OF NODE-NODE NO.= 1STX=-.108791E+01 STY=-.300000E+01 TXY= .439560E+00NODE NO.= 2STX=-.100000E+01 STY=-.220879E+01 TXY= .249084E+00NODE NO.= 3STX=-.105861E+01 STY=

36、-.191575E+01 TXY= .205128E+00NODE NO.= 4STX=-.824176E+00 STY=-.225275E+01 TXY= .000000E+00NODE NO.= 5STX=-.970696E+00 STY=-.166667E+01 TXY= .586081E-01NODE NO.= 6STX=-.100000E+01 STY=-.137363E+01 TXY=-.131868E+00结论与体会:通过本次的课程设计，我对有限元的概念有了更加深刻的理解，同时也弥补了平时学习是疏忽的地方，充实了有限元知识。通过编程及对 chengxu.for 的调试进一步巩固了 FORTRAN 语言知识，也锻炼了自我，使自己在学习工作中更具耐心、细心。参考数目：FORTRAN 语言清华大学出版社 谭浩强，田淑清弹性力学 高等教育出版社 徐芝纶有限元分析及其应用清华大学 曾攀