基于Matlab语言的按平面三角形单元划分的结构有限元程序设计模板

有限元大作业程序设计学校：天津大学院系：建筑工程与力学学院专业：01级工程力学姓名：刘秀学号：\\\\\\\\\\\指导老师：连续体平面问题的有限元程序分析[题目]：如图所示的正方形薄板四周受均匀载荷的作用，该结构在边界上受正向分布压力，m kNp 1=，同时在沿对角线y 轴上受一对集中压力，载荷为2KN ，若取板厚1=t ，泊松比0=v 。

[分析过程]：由于连续平板的对称性，只需要取其在第一象限的四分之一部分参加分析，然后人为作出一些辅助线将平板“分割”成若干部分，再为每个部分选择分析单元。

采用将此模型化分为4个全等的直角三角型单元。

利用其对称性，四分之一部分的边界约束，载荷可等效如图所示。

[程序原理及实现]：用FORTRAN程序的实现。

由节点信息文件NODE.IN和单元信息文件ELEMENT.IN，经过计算分析后输出一个一般性的文件DATA.OUT。

模型基本信息由文件为BASIC.IN生成。

该程序的特点如下：问题类型：可用于计算弹性力学平面问题和平面应变问题单元类型：采用常应变三角形单元位移模式：用用线性位移模式载荷类型：节点载荷，非节点载荷应先换算为等效节点载荷材料性质：弹性体由单一的均匀材料组成约束方式：为“0”位移固定约束，为保证无刚体位移，弹性体至少应有对三个自由度的独立约束方程求解：针对半带宽刚度方程的Gauss消元法输入文件：由手工生成节点信息文件NODE.IN，和单元信息文件ELEMENT.IN结果文件：输出一般的结果文件DATA.OUT程序的原理如框图：（1）主要变量：ID：问题类型码，ID＝1时为平面应力问题，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)表示第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_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的格式如下表需读入的单元信息文件ELEMENT.IN的格式如下表输出结果文件DATA.OUT格式如下表[算例原始数据和程序分析]：（1）模型基本信息文件BASIC.IN的数据为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.02 0.0 1.03 1.0 1.04 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 04 （4）源程序文件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)D IMENSION 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_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,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 fileWRITE(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 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),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)E NDcc 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,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=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 IFR ETURNENDcC to form the stiffness matrix of elementSUBROUTINE 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,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*PTRETURNENDcc to form banded global stiffness matrixSUBROUTINE BAND_K(N_DOF,N_BAND,N_ELE,IE,N_NODE,IJK_ELE,X,Y,PE, & PR,PT,AK)DIMENSIONIJK_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)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 CONTINUERETURNENDcc to calculate the area of elementSUBROUTINE 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.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.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,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(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)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),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 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_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_BAND10 AK(II,JJ)=0.0DO 20 JJ=2,II20 AK(II-JJ+1,JJ)=0.0END IF30 CONTINUERETURNENDcc 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_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(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,'* * * * * 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,RESULT_N(2*I-1),RESULT_N(2*I)70 FORMAT(8X,I5,7X,2E15.6)RETURNENDcc calculate the stress components of element and nodeSUBROUTINECAL_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)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,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(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[算例结果]：chengxu.for所输出的数据文件DATA.OUT数据内容如下：=========PLANE STRESS PROBLEM========* * * * * RESULTS BY FEM2D * * * * *--DISPLACEMENT OF NODE--NODE NO X-DISP Y-DISP1 .000000E+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+00ELEMENT 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=-.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[结论与体会]:通过本次的课程设计，我对有限元的概念有了更加深刻的理解，同时也弥补了平时学习是疏忽的地方，充实了有限元知识。

基于Matlab说话的按平面三角形单元划分的构造有限元程序设计专业：建筑与土木匠程班级：建工研12-2姓名：韩志强学号：471220580基于Matlab说话的按平面三角形单元划分构造有限元程序设计一、有限单元发及Matlab说话概述1. 有限单元法跟着现代工业.临盆技巧的成长,不竭请求设计高质量.高程度的大型.庞杂和周详的机械及工程构造.为此目标,人们必须预先经由过程有用的盘算手腕,确实的猜测即将诞生的机械和工程构造,在将来工作时所产生的应力.应变和位移是以,须要追求一种简略而又准确的数值剖析办法.有限单元法恰是顺应这种请求而产生和成长起来的一种十分有用的数值盘算办法.有限元法把一个庞杂的构造分化成相对简略的“单元”,各单元之间经由过程结点互相衔接.单元内的物理量由单元结点上的物理量按必定的假设内插得到,如许就把一个庞杂构造从无穷多个自由度简化为有限个单元构成的构造.我们只要剖析每个单元的力学特征,然后按照有限元法的规矩把这些单元“拼装”成整体,就可以或许得到整体构造的力学特征.有限单元法根本步调如下：(1)构造离散：构造离散就是树立构造的有限元模子,又称为网格划分或单元划分,即将构造离散为由有限个单元构成的有限元模子.在该步调中,须要依据构造的几何特征.载荷情形等肯定单元体内随意率性一点的位移插值函数.(2)单元剖析：依据弹性力学的几何方程以及物理方程肯定单元的刚度矩阵.(3)整体剖析:把各个单元按本来的构造从新衔接起来,并在单元刚度矩阵的基本上肯定构造的总刚度矩阵,形成如下式所示的整体有限元线性方程：式中阵.(4)载荷移置:依据静力等效道理,将载荷移置到响应的节点上,形成节点载荷矩阵.(5)鸿沟前提处理:对式①所示的有限元线性方程进行鸿沟前提处理.(6)求解线性方程:求解式①所示的有限元线性方程,得到节点的位移.在该步调中,如有限元模子的节点越多,则线性方程的数目就越多,随之有限元剖析的盘算量也将越大.(7)求解单元应力及应变依据求出的节点位移求解单元的应力和应变.(8)成果处理与显示进入有限元剖析的后处理部分,对盘算出来的成果进行加工处理,并以各类情势将盘算成果显示出.在用有限元法进行构造剖析时,将会碰到大量的数值盘算,因而在实用上是必定要借助于盘算机和有限元程序,才干完成这些庞杂而沉重的数值盘算工作.而Matlab是当今国际科学界最具影响力和活气的软件.它来源于矩阵运算,并已经成长成一种高度集成的盘算机说话.它供给了壮大的科学盘算,灵巧的程序设计流程,高质量的图形可视化与界面设计,便捷的与其他程序和说话接口的功效.Matlab在列国高校与研讨单位起侧重大的感化.“工欲善其事,必先利其器”.假如有一种十分有用的对象能解决在教授教养与研讨中碰到的问题,那么Matlab说话恰是如许的一种对象.它可以将运用者从繁琐.无谓的底层编程中解放出来,把有限的珍贵时光更多地花在解决问题中,如许无疑会进步工作效力. 今朝,Matlab已经成为国际上最风行的科学与工程盘算的软件对象,如今的Matlab已经不但仅是一个“矩阵试验室”了,它已经成为了一种具有普遍运用远景的全新的盘算机高等编程说话了,有人称它为“第四代”盘算机说话,它在国表里高校和研讨部分正扮演侧重要的脚色.Matlab说话的功效也越来越壮大,不竭顺应新的请求提出新的解决办法.可以预感,在科学运算.主动掌握与科学画图范畴Matlab说话将长期保持其举世无双的地位.为此,本例采取Matlab说话编程,以运用其快捷壮大的矩阵数值盘算功效.二、问题描写一矩形薄板,一边固定,推却150kN分散荷载,构造简图及按平面三角形单元划分的有限元模子图如下所示.薄板厚在本例中,所取构造模子及数据重要用于程序设计理论剖析,与工程现实无关.三、参数输入：单元个数NELEM=4节点个数NNODE=6受束缚鸿沟点数NVFIX=2节点荷载个数NFORCE=1弹性模量YOUNG=2e8泊松比POISS=0.2厚度THICK=0.002单元节点编码数组节点坐标数组节点力数组FORCE =[4 0 -150]束缚信息数组以上数值数据为程序运行前输入的初始数据,存为“”文本格局,同时必须放在Matlab工作目次下,路径不合错误程序不克不及主动读取指定初始文件,运行出错.初始数据文本文件输入格局如下图：四、Matlab说话程序源代码：1.程序中变量解释NNODE 单元节点数NPION 总结点数NELEM 单元数NVFIX 受束缚鸿沟点数FIXED 束缚信息数组NFORCE 节点力数FORCE 节点力数组COORD 构造节点坐标数组LNODS 单元界说数组YOUNG 弹性模量POISS 泊松比THICK 厚度B 单元应变矩阵(3*6) D 单元弹性矩阵(3*3) S 单元应力矩阵(3*6) A单元面积ESTIF 单元刚度矩阵ASTIF 总体刚度矩阵ASLOD 总体荷载向量ASDISP 节点位移向量ELEDISP 单元节点位移向量STRESS 单元应力FP1 数据文件指针2.Matlab说话程序代码%****************************************************************************** %初始化及数据挪用format short e %设定输出类型clear %消除内存变量FP1=fopen('','rt'); %打开输入数据文件,读入掌握数据NELEM=fscanf(FP1,'%d',1), %单元个数（单元编码总数）NPION=fscanf(FP1,'%d',1), %结点个数（结点编码总数）NVFIX=fscanf(FP1,'%d',1)%受束缚鸿沟点数NFORCE=fscanf(FP1,'%d',1), %结点荷载个数YOUNG=fscanf(FP1,'%e',1), %弹性模量POISS=fscanf(FP1,'%f',1), %泊松比THICK=fscanf(FP1,'%f',1) %厚度LNODS=fscanf(FP1,'%d',[3,NELEM])' %单元界说数组（单元结点号）%响应为单元结点号（编码）.按逆时针次序输入COORD=fscanf(FP1,'%f',[2,NPION])' %结点坐标数组%坐标：x,y坐标（共NPOIN组）FORCE=fscanf(FP1,'%f',[3,NFORCE])' %结点力数组（受力结点编号,x偏向,y偏向）FIXED=fscanf(FP1,'%d',[3,NVFIX])' %束缚信息（束缚点,x束缚,y束缚）%有束缚为1,无束缚为0%***************************************************************************** %生成单元刚度矩阵并构成总体刚度矩阵ASTIF=zeros(2*NPION,2*NPION); %生成特定大小总体刚度矩阵并置0%***************************************************************************** for i=1:NELEM%生成弹性矩阵DD= [1 POISS 0;POISS 1 0;0 0 (1-POISS)/2]*YOUNG/(1-POISS^2)%***************************************************************************** %盘算当前单元的面积AA=det([1 COORD(LNODS(i,1),1) COORD(LNODS(i,1),2);1 COORD(LNODS(i,2),1) COORD(LNODS(i,2),2);1 COORD(LNODS(i,3),1) COORD(LNODS(i,3),2)])/2%***************************************************************************** %生成应变矩阵Bfor j=0:2b(j+1)= COORD(LNODS(i,(rem((j+1),3))+1),2)-COORD(LNODS(i,(rem((j+2),3))+1),2);c(j+1)=-COORD(LNODS(i,(rem((j+1),3))+1),1)+COORD(LNODS(i,(rem((j+2),3))+1),1);endB=[b(1) 0 b(2) 0 b(3) 0;...0 c(1) 0 c(2) 0 c(3);...c(1) b(1) c(2) b(2) c(3) b(3)]/(2*A);B1( :,:,i)=B;%***************************************************************************** %求应力矩阵S=D*BS=D*B;ESTIF=B'*S*THICK*A; %求解单元刚度矩阵a=LNODS(i,:); %暂时向量,用来记载当前单元的节点编号for j=1:3for k=1:3ASTIF((a(j)*2-1):a(j)*2,(a(k)*2-1):a(k)*2)…=ASTIF((a(j)*2-1):a(j)*2,(a(k)*2-1):a(k)*2)+ESTIF(j*2-1:j*2,k*2-1:k*2); %依据节点编号对应关系将单元刚度分块叠加到总刚%度矩阵中endendend%***************************************************************************** %将束缚信息参加总体刚度矩阵（对角元素改一法）for i=1:NVFIXif FIXED(i,2)==1ASTIF(:,(FIXED(i,1)*2-1))=0; %一列为零ASTIF((FIXED(i,1)*2-1),:)=0; %一行动零ASTIF((FIXED(i,1)*2-1),(FIXED(i,1)*2-1))=1;%对角元素为1endif FIXED(i,3)==1ASTIF( :,FIXED(i,1)*2)=0; %一列为零ASTIF(FIXED(i,1)*2,:)=0; %一行动零ASTIF(FIXED(i,1)*2 ,FIXED(i,1)*2)=1; %对角元素为1endend%***************************************************************************** %生成荷载向量ASLOD(1:2*NPION)=0; %总体荷载向量置零for i=1:NFORCEASLOD((FORCE(i,1)*2-1):FORCE(i,1)*2)=FORCE(i,2:3);end%***************************************************************************** %求解内力ASDISP=ASTIF\ASLOD' %盘算节点位移向量ELEDISP(1:6)=0;%当前单元节点位移向量for i=1:NELEMfor j=1:3ELEDISP(j*2-1:j*2)=ASDISP(LNODS(i,j)*2-1:LNODS(i,j)*2);%掏出当前单元的节点位移向量endiSTRESS=D*B1(:, :, i)*ELEDISP' %求内力endfclose(FP1); %封闭数据文件五、程序运行成果NELEM =4NPION =6NVFIX =2NFORCE =1YOUNG =200000000POISS =THICK =LNODS =1 2 62 3 42 4 52 5 6COORD =0 01 02 02 11 10 1FORCE =4 0 -150FIXED =1 1 16 1 1D =2.0833e+008 4.1667e+007 0 4.1667e+007 2.0833e+008 00 0 8.3333e+007A =D =2.0833e+008 4.1667e+007 0 4.1667e+007 2.0833e+008 00 0 8.3333e+007A =D =2.0833e+008 4.1667e+007 0 4.1667e+007 2.0833e+008 00 0 8.3333e+007A =D =2.0833e+008 4.1667e+007 0 4.1667e+007 2.0833e+008 0 0 0 8.3333e+007A =ASDISP =（解释：以上12个ASDISP输出成果数据暗示节点位移向量,即各节点分离在X,Y偏向上产生的位移.）i =1STRESS =-1.6793e+005-3.3586e+004-1.3207e+005i =2STRESS =-5.5503e+004-5.5503e+004-5.5503e+004i =3STRESS =5.5503e+004-4.9526e+004-9.4497e+004i =4STRESS =1.6793e+005-2.7040e+004-1.7932e+004（解释：以上四组STRESS输出成果数据暗示单元应力SX,SY,SXY,i为单元号.）六、ANSYS建模比较剖析运用ANSYS有限元剖析软件,完整按照前面Matlab程序设计的各项前提参数以及单元划分方法树立ANSYS模子,其求解成果与以上程序盘算成果比较,验证程序是否准确.ANSYS模子节点单元如下图所示：ANSYS模子求解变形图如下所示：ANSYS求解节点位移成果列表显示如下：（单位：m）PRINT U NODAL SOLUTION PER NODE***** POST1 NODAL DEGREE OF FREEDOM LISTING *****THE FOLLOWING DEGREE OF FREEDOM RESULTS ARE IN THE GLOBAL COORDINATESYSTEMNODE UX UY UZ USUM1 0.0000 0.0000 0.0000 0.00006 0.0000 0.0000 0.0000 0.0000MAXIMUM ABSOLUTE VALUESNODE 4 4 0 4ANSYS求解单元应力成果列表显示如下：（单位：KN/m2）PRINT S ELEMENT SOLUTION PER ELEMENT***** POST1 ELEMENT NODAL STRESS LISTING *****THE FOLLOWING X,Y,Z VALUES ARE IN GLOBAL COORDINATESELEMENT= 1 PLANE182NODE SX SY SZ SXY SYZ1 -0.16793E+06 -33586. 0.0000 -0.13207E+06 0.00002 -0.16793E+06 -33586. 0.0000 -0.13207E+06 0.0000 6 -0.16793E+06 -33586. 0.0000 -0.13207E+06 0.0000 ELEMENT= 2 PLANE182NODE SX SY SZ SXY SYZ2 -55503. -55503. 0.0000 -55503. 0.00003 -55503. -55503. 0.0000 -55503. 0.00004 -55503. -55503. 0.0000 -55503. 0.0000 ELEMENT= 3 PLANE182NODE SX SY SZ SXY SYZ2 55503. -49526. 0.0000 -94497. 0.00004 55503. -49526. 0.0000 -94497. 0.00005 55503. -49526. 0.0000 -94497. 0.0000 ELEMENT= 4 PLANE182NODE SX SY SZ SXY SYZ2 0.16793E+06 -27040. 0.0000 -17932. 0.00005 0.16793E+06 -27040. 0.0000 -17932. 0.00006 0.16793E+06 -27040. 0.0000 -17932. 0.0000结论经由过程比较Matlab说话设计程序运行成果和ANSYS建模剖析成果可知,两种方法算出的成果完整一致,说程序设计准确.所以本程序实用于按三角形单元划分的平面构造有限元剖析.format short eclearFP1=fopen('LinearTriangleElement of Thin plate.txt','rt');NELEM=fscanf(FP1,'%d',1)NPION=fscanf(FP1,'%d',1)NVFIX=fscanf(FP1,'%d',1)NFORCE=fscanf(FP1,'%d',1)YOUNG=fscanf(FP1,'%e',1)POISS=fscanf(FP1,'%f',1)THICK=fscanf(FP1,'%f',1)LNODS=fscanf(FP1,'%d',[3,NELEM])COORD=fscanf(FP1,'%f',[2,NPION])FORCE=fscanf(FP1,'%f',[3,NFORCE])FIXED=fscanf(FP1,'%d',[3,NVFIX])ASTIF=zeros(2*NPION,2*NPION); %生成特定大小总体刚度矩阵并置0for i=1:NELEM%生成弹性矩阵DD= YOUNG/(1-POISS^2)*[1 POISS 0;POISS 1 0;0 0 (1-POISS)/2]%盘算当前单元的面积AA=det([1 COORD(LNODS(i,1),1) COORD(LNODS(i,1),2);1 COORD(LNODS(i,2),1) COORD(LNODS(i,2),2);1 COORD(LNODS(i,3),1) COORD(LNODS(i,3),2)])/2%生成应变矩阵Bfor j=0:2b(j+1)= COORD(LNODS(i,(rem((j+1),3))+1),2)-COORD(LNODS(i,(rem((j+2),3))+1),2);c(j+1)=-COORD(LNODS(i,(rem((j+1),3))+1),1)+COORD(LNODS(i,(rem((j+2),3))+1),1); endB=[b(1) 0 b(2) 0 b(3) 0;...0 c(1) 0 c(2) 0 c(3);...c(1) b(1) c(2) b(2) c(3) b(3)]/(2*A);B1( :,:,i)=B;%求应力矩阵S=D*BS=D*B;ESTIF=B'*S*THICK*A;a=LNODS(i,:);for j=1:3for k=1:3ASTIF((a(j)*2-1):a(j)*2,(a(k)*2-1):a(k)*2)…=ASTIF((a(j)*2-1):a(j)*2,(a(k)*2-1):a(k)*2)+ESTIF(j*2-1:j*2,k*2-1:k*2);endendend%将束缚信息参加总体刚度矩阵for i=1:NVFIXif FIXED(i,2)==1ASTIF(:,(FIXED(i,1)*2-1))=0;ASTIF((FIXED(i,1)*2-1),:)=0;ASTIF((FIXED(i,1)*2-1),(FIXED(i,1)*2-1))=1;endif FIXED(i,3)==1ASTIF( :,FIXED(i,1)*2)=0;ASTIF(FIXED(i,1)*2,:)=0;ASTIF(FIXED(i,1)*2 ,FIXED(i,1)*2)=1;endend%生成荷载向量ASLOD(1:2*NPION)=0;for i=1:NFORCEASLOD((FORCE(i,1)*2-1):FORCE(i,1)*2)=FORCE(i,2:3)end%求解内力ASDISP=ASTIF\ASLOD'ELEDISP(1:6)=0;for i=1:NELEMfor j=1:3ELEDISP(j*2-1:j*2)=ASDISP(LNODS(i,j)*2-1:LNODS(i,j)*2);endiSTRESS=D*B1(:, :, i)*ELEDISP'endfclose(FP1);format short eclearFP1=fopen('LinearTriangleElement of Thin plate.txt','rt'); NELEM=fscanf(FP1,'%d',1)NPION=fscanf(FP1,'%d',1)NVFIX=fscanf(FP1,'%d',1)NFORCE=fscanf(FP1,'%d',1)YOUNG=fscanf(FP1,'%e',1)POISS=fscanf(FP1,'%f',1)THICK=fscanf(FP1,'%f',1)LNODS=fscanf(FP1,'%d',[NELEM,3])COORD=fscanf(FP1,'%f',[NPION,2])FORCE=fscanf(FP1,'%f',[NFORCE,3])FIXED=fscanf(FP1,'%d',[NVFIX,3])ASTIF=zeros(2*NPION,2*NPION); %生成特定大小总体刚度矩阵并置0 for i=1:NELEM%生成弹性矩阵DD= YOUNG/(1-POISS^2)*[1 POISS 0;POISS 1 0;0 0 (1-POISS)/2]%盘算当前单元的面积AA=det([1 COORD(LNODS(i,1),1) COORD(LNODS(i,1),2);1 COORD(LNODS(i,2),1) COORD(LNODS(i,2),2);1 COORD(LNODS(i,3),1) COORD(LNODS(i,3),2)])/2%生成应变矩阵Bfor j=0:2b(j+1)= COORD(LNODS(i,(rem((j+1),3))+1),2)-COORD(LNODS(i,(rem((j+2),3))+1),2);c(j+1)=-COORD(LNODS(i,(rem((j+1),3))+1),1)+COORD(LNODS(i,(rem((j+2),3))+1),1); endB=[b(1) 0 b(2) 0 b(3) 0;...0 c(1) 0 c(2) 0 c(3);...c(1) b(1) c(2) b(2) c(3) b(3)]/(2*A);B1( :,:,i)=B;%求应力矩阵S=D*BS=D*B;ESTIF=B'*S*THICK*A;a=LNODS(i,:);for j=1:3for k=1:3ASTIF((a(j)*2-1):a(j)*2,(a(k)*2-1):a(k)*2)…=ASTIF((a(j)*2-1):a(j)*2,(a(k)*2-1):a(k)*2)+ESTIF(j*2-1:j*2,k*2-1:k*2);endendend%将束缚信息参加总体刚度矩阵for i=1:NVFIXif FIXED(i,2)==1ASTIF(:,(FIXED(i,1)*2-1))=0;ASTIF((FIXED(i,1)*2-1),:)=0;ASTIF((FIXED(i,1)*2-1),(FIXED(i,1)*2-1))=1;endif FIXED(i,3)==1ASTIF( :,FIXED(i,1)*2)=0;ASTIF(FIXED(i,1)*2,:)=0;ASTIF(FIXED(i,1)*2 ,FIXED(i,1)*2)=1;endend%生成荷载向量ASLOD(1:2*NPION)=0;for i=1:NFORCEASLOD((FORCE(i,1)*2-1):FORCE(i,1)*2)=FORCE(i,2:3)end%求解内力ASDISP=ASTIF\ASLOD'ELEDISP(1:6)=0;for i=1:NELEMfor j=1:3ELEDISP(j*2-1:j*2)=ASDISP(LNODS(i,j)*2-1:LNODS(i,j)*2); endiSTRESS=D*B1(:, :, i)*ELEDISP'endfclose(FP1);。

平面三角形单元有限元程序设计P9 m 9 m一、题目如图1所示,一个厚度均匀的三角形薄板,在顶点作用沿板厚方向均匀分布的竖向载荷。

已知:P=150N/m,E=200GPa,=0、25,t=0、1m,忽略自重。

试计算薄板的位移及应力分布。

要求:1.编写有限元计算机程序,计算节点位移及单元应力。

(划分三角形单元,单元数不得少于30个);2.采用有限元软件分析该问题(有限元软件网格与程序设计网格必须一致),详细给出有限元软件每一步的操作过程,并将结果与程序计算结果进行对比(任选取三个点,对比位移值);3.提交程序编写过程的详细报告及计算机程序;4.所有同学参加答辩,并演示有限元计算程序。

有限元法中三节点三角形分析结构的步骤如下:1)整理原始数据,如材料性质、荷载条件、约束条件等,离散结构并进行单元编码、结点编码、结点位移编码、选取坐标系。

2)单元分析,建立单元刚度矩阵。

3)整体分析,建立总刚矩阵。

4)建立整体结构的等效节点荷载与总荷载矩阵5)边界条件处理。

6)解方程,求出节点位移。

7)求出各单元的单元应力。

8)计算结果整理。

一、程序设计网格划分如图,将薄板如图划分为6行,并建立坐标系,则刚度矩阵的集成建立与总刚度矩阵等维数的空矩阵,已变单元刚度矩阵的集成。

由单元分析已知节点、单元的排布规律,继而通过循环计算求得每个单元对应的节点序号。

通过循环逐个计算:(1)每个单元对应2种单元刚度矩阵中的哪一种; (2)该单元对应总刚度矩阵的那几行哪几列 (3)将该单元的单元刚度矩阵加入总刚度矩阵的对应行列循环又分为3层循环:(1)最外层:逐行计算(2)中间层:该行逐个计算(3)最里层:区分为第 奇/偶 数个计算XYPXYP节点编号单元编号单元刚度的集成:[][][][][][]''''''215656665656266256561661eZeeeZeZeeeekkkKkkkkkk+⋯++=⇓=⇒==⇒==⇒=⨯⨯⨯⨯⨯⨯M边界约束的处理:划0置1法适用:这种方法适用于边界节点位移分量为已知(含为0)的各种约束。

计算力学(有限元方法部分) 程序设计程序说明书程序1：平面问题的有限元分析文件名：h01.m算例文本：h01.txt输出文本：result1.txt使用前请先将h01.txt放入默认的文本读取路径（我的要求与m文件在同一文件夹内）！文本输入顺序：材料信息（编号、弹性模量、泊松比）注意：材料编号须按1、2、3、4……的顺序排列节点信息（编号、x坐标、y坐标）注意：节点编号须按1、2、3、4……的顺序排列约束信息（约束节点号、x方向有无约束、y方向有无约束、x方向位移、y 方向位移）有约束处填一正数，无约束处填0。

无约束处请勿输入位移。

单元信息（厚度、材料编号、节点编号，若为3节点单元，则第四个编号填0）注意：单元编号须按1、2、3、4……的顺序排列集中力（作用节点号、x方向力、y方向力）线荷载（作用边上的两个节点号、x方向分布力、y方向分布力）面荷载（作用单元号、x方向分布力、y方向分布力）程序可调部分：4-6行中可以指定输出哪些图像（按顺序依次为节点、单元图像，x、y方向位移图像，xx、yy、xy方向应力图像），第7行中可以指定输入的.txt文本名称。

文本输出内容：结点位移信息（节点号、x方向位移、y方向位移）单元形心处的应变信息（单元号、x方向正应变、y方向正应变、xy方向工程切应变）单元形心处的应力信息（单元号、x方向正应力、y方向正应力、xy方向切应力）本程序附有三角形单元自动加密前处理部分h01auto.m，其算例文本：h01coarse.txt，输出文本：h01new.txt。

它可以适用于本题的要求，在已给定本题所需全部信息的情况下将已有的单元加密为三角形单元。

其输出文本可直接作为上面程序的输入文本。

h01.m h01.txt h01auto.m h01coarse.txt欢迎交流与提问！附上邮箱：***************。

祝力学学习顺利！。

matlab有限元三单元编程MATLAB是一种功能强大的编程语言和环境，广泛用于工程和科学领域的数值计算、数据分析和可视化。

在有限元分析中，MATLAB的强大功能和直观的语法使其成为一个理想的选择。

在本文中，我们将讨论MATLAB中有限元三单元的编程方法和实践。

有限元分析是一种数值方法，用于解决连续介质的力学问题。

它将一个复杂的结构分解成更简单的有限元单元，然后通过求解线性代数方程组来得到结构的应力和位移解。

在有限元分析中，三角形和四边形单元是最常用的有限元单元之一。

本文将重点讨论三角形单元的编程实现。

首先，我们需要定义一个三角形单元的几何信息。

在三角形单元中，我们有三个顶点，每个顶点有两个坐标。

我们可以使用一个3x2的矩阵来表示这些坐标。

例如，定义三角形ABC的顶点坐标矩阵为P：P = [x_A, y_A;x_B, y_B;x_C, y_C];接下来，我们需要定义三角形单元的连接性信息。

在MATLAB中，我们可以使用一个3x1的向量来表示三个顶点的连接性。

例如，定义三角形ABC的连接性向量为E：E = [1;2;3];然后，我们可以定义三角形单元材料属性和载荷信息。

这些信息包括杨氏模量、泊松比和外力向量。

我们可以将这些信息存储在一个结构体中，例如：properties.E = 210e9; % 杨氏模量properties.v = 0.3; % 泊松比properties.f = [0; -1000; 0]; % 外力向量接下来，我们可以使用这些信息来计算三角形单元的刚度矩阵和力向量。

刚度矩阵是一个3x3的矩阵，力向量是一个3x1的向量。

我们可以使用以下公式来计算它们：K = (E/(1-v^2)) * [1 v 0;v 1 0;0 0 (1-v)/2] * A;f = (A/3) * [1; 1; 1] * properties.f;其中，A是三角形的面积。

我们可以使用以下公式来计算它：A = abs(det([1, 1, 1;x_A, x_B, x_C;y_A, y_B, y_C])) / 2;最后，我们可以使用这些信息来求解三角形单元的位移解。

matlab有限元三角形单元程序以下是一个简单的 MATLAB 有限元三角形单元程序的示例：```matlab% 定义模型参数E = 1000; % 弹性模量nu = 0.3; % 泊松比thickness = 1; % 板的厚度% 定义节点坐标nodes = [0, 0; 1, 0; 1, 1; 0, 1; 0.5, 0.5];% 定义单元连接关系connectivity = [1, 2, 5; 2, 3, 5; 3, 4, 5; 4, 1, 5];% 计算总节点数和总单元数numNodes = size(nodes, 1);numElements = size(connectivity, 1);% 初始化全局刚度矩阵和载荷向量K = zeros(numNodes);F = zeros(numNodes, 1);% 循环遍历每个单元for i = 1:numElements% 查找当前单元的节点编号nodesIndex = connectivity(i, :);% 根据节点编号从全局坐标矩阵中取出节点坐标coordinates = nodes(nodesIndex, :);% 计算当前单元的局部刚度矩阵localK = calculateLocalStiffness(E, nu, thickness, coordinates);% 组装局部刚度矩阵到全局刚度矩阵中K(nodesIndex, nodesIndex) = K(nodesIndex, nodesIndex) + localK;% 计算当前单元的局部载荷向量localF = calculateLocalLoad(thickness, coordinates);% 组装局部载荷向量到全局载荷向量中F(nodesIndex) = F(nodesIndex) + localF;end% 边界条件：节点1固定K(1, :) = 0;K(1, 1) = 1;F(1) = 0;% 解线性方程组U = K \ F;% 输出位移结果disp('节点位移：');disp(U);% 计算应力结果stress = calculateStress(E, nu, thickness, nodes, connectivity, U);% 输出应力结果disp('节点应力：');disp(stress);% 计算局部刚度矩阵的函数function localK = calculateLocalStiffness(E, nu, thickness, coordinates)% 计算单元的雅可比矩阵J = (1/2) * [coordinates(2,1)-coordinates(1,1), coordinates(3,1)-coordinates(1,1);coordinates(2,2)-coordinates(1,2), coordinates(3,2)-coordinates(1,2)];% 计算雅可比矩阵的逆矩阵invJ = inv(J);% 计算单元刚度矩阵B = invJ * [-1, 1, 0; -1, 0, 1];D = (E/(1-nu^2)) * [1, nu, 0; nu, 1, 0; 0, 0, (1-nu)/2]; localK = thickness * abs(det(J)) * (B' * D * B);end% 计算局部载荷向量的函数function localF = calculateLocalLoad(thickness, coordinates) localF = zeros(3, 1);midPoint = [sum(coordinates(:,1))/3,sum(coordinates(:,2))/3];localF(3) = thickness * 1 * det([coordinates(1,:); coordinates(2,:); midPoint]);end% 计算各节点应力的函数function stress = calculateStress(E, nu, thickness, nodes, connectivity, U)stress = zeros(size(nodes, 1), 3);for i = 1:size(connectivity, 1)nodesIndex = connectivity(i, :);coordinates = nodes(nodesIndex, :);Ke = calculateLocalStiffness(E, nu, thickness, coordinates);Ue = U(nodesIndex);stress(nodesIndex, :) = stress(nodesIndex, :) + (Ke * Ue)';endstress = stress / thickness;end```这个程序实现了一个简单的平面三角形单元有限元分析，包括定义节点坐标和单元连接关系、计算全局刚度矩阵和载荷向量、施加边界条件、解线性方程组、计算节点位移和应力等。

三角形常应变单元程序的编制与使用有限元法是求解微分方程边值问题的一种通用数值方法，该方法是一种基于变分法（或变分里兹法）而发展起来的求解微分方程的数值计算方法，以计算机为手段，采用分片近似，进而逼近整体的研究思想求解物理问题。

有限元分析的基本步骤可归纳为三大步：结构离散、单元分析和整体分析。

对于平面问题，结构离散常用的网格形状有三角形、矩形、任意四边形，以三个顶点为节点的三角形单元是最简单的平面单元，它较矩形或四边形对曲边边界有更好的适应性，而矩形或四边形单元较三节点三角形有更高的计算精度。

Matlab语言是进行矩阵运算的强大工具，因此，用Matlab语言编写有限元中平面问题的程序有优越性。

本章将详细介绍如何利用Matlab语言编制三角形常应变单元的计算程序，程序流程图见图1。

有限元法中三节点三角形分析结构的步骤如下：1）整理原始数据，如材料性质、荷载条件、约束条件等，离散结构并进行单元编码、结点编码、结点位移编码、选取坐标系。

2）单元分析，建立单元刚度矩阵。

3）整体分析，建立总刚矩阵。

4）建立整体结构的等效节点荷载和总荷载矩阵5）边界条件处理。

6）解方程，求出节点位移。

7）求出各单元的单元应力。

8）计算结果整理。

计算结果整理包括位移和应力两个方面；位移计算结果一般不需要特别的处理，利用计算出的节点位移分量，就可画出结构任意方向的位移云图；而应力解的误差表现在单元内部不满足平衡方程，单元与单元边界处应力一般不连续，在边界上应力解一般与力的边界条件不相符合。

图1 程序流程图1.1 程序说明%******************************************************************* % 三角形常应变单元求解结构主程序%******************************************************************* ●功能：运用有限元法中三角形常应变单元解平面问题的计算主程序。

一、平面三角形３结点有限元程序1、程序名：，FEM3.EXE2、程序功能本程序能计算弹性力学的平面应力问题和平面应变问题；能考虑自重和结点集中力两种荷载的作用，在计算自重时y轴取垂直向上为正；能处理非零已知位移，如支座沉降的作用。

主要输出的内容包括：结点位移、单元应力、主应力、第一主应力与x轴的夹角以及约束结点的支座反力。

程序采用Visual Fortran 5.0编制而成，输入数据全部采用自由格式。

3、程序流程及框图图１-1 程序流程图图１-2 程序框图其中，各子程序的功能如下：INPUT——输入结点坐标、单元信息和材料参数；MR——形成结点自由度序号矩阵；FORMMA——形成指标矩阵MA（N）并调用其他功能子程序，相当于主控程序；DIV——取出单元的3个结点号码和该单元的材料号并计算单元的b i，c i等；MGK——形成整体劲度矩阵并按一维存放在SK（NH）中；LOAD——形成整体结点荷载列阵F；OUTPUT——输出结点位移或结点荷载；TREAT——由于有非零已知位移，对K和F进行处理；DECOMP——整体劲度矩阵的分解运算；FOBA——前代、回代求出未知结点位移 ；ERFAC——计算约束结点的支座反力；KRS——计算单元劲度矩阵中的子块K rs。

4、输入数据及变量说明当程序开始运行时，按屏幕提示，键入数据文件的名字。

在运行程序之前，必须根据程序中输入要求建立一个存放原始数据的文件，这个文件的名字由少于8个字符或数字组成。

数据文件包括如下内容：⑴总控信息，共一条，9个数据NP，NE，NM，NR，NI，NL，NG，ND，NCNP——结点总数；NE——单元总数；NM——材料类型总数；NR ——约束结点总数；NI ——问题类型标识，0为平面应力问题，1为平面应变问题；NL ——受荷载作用的结点的数目；NG ——考虑自重作用为1，不计自重为0；ND ——非零已知位移结点的数目；NC ——要计算支座约束力的结点数目。

%变量说明%NPOIN NELEM NVFIX NFORCE NNODE%总结点数，单元数，受约束边界点数，结点力数, 单元结点数%COORD LNODS YOUNG POISS THICK %结构结点坐标数组，单元定义数组，弹性模量，泊松比，厚度%FORCE FIXED BMATX DMATX SMATX %节点力数组，约束信息数组，单元应变矩阵，单元弹性矩阵，单元应力矩阵%AREA HK ASLOD ASDISP FP1%单元面积，总体刚度矩阵，总体荷载向量，结点位移向量，数据文件指针format short eclearFP1=fopen('C:\Users\Administrator\Desktop\input.txt','rt');NPION=fscanf(FP1,'%d',1); %结点个数（结点编码总数）NELEM=fscanf(FP1,'%d',1); %单元个数（单元编码总数）NFORCE=fscanf(FP1,'%d',1); %结点荷载个数NVFIX=fscanf(FP1,'%d',1); %受约束边界点数YOUNG=fscanf(FP1,'%e',1); %弹性模量POISS=fscanf(FP1,'%f',1); %泊松比THICK=fscanf(FP1,'%d',1); %厚度LNODS=fscanf(FP1,'%d',[3,NELEM])'; %单元定义数组（单元结点号）COORD=fscanf(FP1,'%f',[2,NPION])'; %结点坐标数组FORCE=fscanf(FP1,'%f',[3,NFORCE])'; %结点力数组FIXED=fscanf(FP1,'%d',[3,inf])'; %约束信息数组%引用所需的全局变量%global NPION NELEM COORD LNODS YOUNG POISS %global BMATX DMATX SMATX AREA%生成弹性矩阵Da=YOUNG/(1-POISS^2);DMATX(1,1)=1*a;DMATX(1,2)=POISS*a;DMATX(2,1)=POISS*a;DMATX(2,2)=1*a;DMATX(3,3)=(1-POISS)*a/2;for i=1:NELEM; %i为当前所计算的单元号%计算当前单元的面积AREA=det([1 COORD(LNODS(i,1),1) COORD(LNODS(i,1),2);...1 COORD(LNODS(i,2),1) COORD(LNODS(i,2),2);...1 COORD(LNODS(i,3),1) COORD(LNODS(i,3),2);])/2; end%生成应变矩阵Bfor j=0:2b(j+1)=COORD(LNODS(i,(rem((j+1),3))+1),2)-COORD(LNODS(i,(rem((j+2),3))+1),2) ;c(j+1)=-COORD(LNODS(i,(rem((j+1),3))+1),1)+COORD(LNODS(i,(rem((j+2),3))+1), 1);endBMATX=[b(1) 0 b(2) 0 b(3) 0;...0 c(1) 0 c(2) 0 c(3);...c(1) b(1) c(2) b(2) c(3) b(3)]/(2*AREA);SMATX=DMATX*BMATX; %求应力矩阵S=D*B% 所引用的全局变量：%global NPION NELEM NVFIX LNODS ASTIF THICK%global BMATX SMATX AREA FIXEDHK=seros(2*NPION,2*NPION); %张成特定大小总刚矩阵并置0for i=1:NELEMEK=BMATX'*SMATX*THICK*AREA; %求解单元刚度矩阵a=LNODS(i,:); %临时向量,用来记录当前单元的节点编号for j=1:3for k=1:3HK((a(j)*2-1):a(j)*2,(a(k)*2-1):a(k)*2)=HK((a(j)*2-1):a(j)*2,(a(k)*2-1):a(k)*2)+EK(j*2-1:j*2,k*2-1:k*2);%跟据节点编号对应关系将单元刚度分块叠加到总刚矩阵中endendend%将约束信息加入总刚（置0置1法）NUM=1; %计数器(当前已分析的节点数)i=0; %计数器(当前已处理的约束数)tmp(NVFIX)=0; %临时存被约束的序号while i<NVFIXfor j=-1:0if FIXED(NUM,j+3)==1 %若发现约束i=i+1; %计数器自增tmp(i)=FIXED(NUM)*2+j; %求约束序号endendNUM=NUM+1;endfor i=1:NVFIXHK(1:2*NPION,tmp(i))=0; %将一约束序号处的总刚列向量清0 HK(tmp(i),1:2*NPION)=0; %.将一约束序号处的总刚行向量清0 HK(tmp(i),tmp(i))=1; %将行列交叉处的元素置为1end%所需引用的全局变量%global ASLOD NPION NFORCE FORCEASLOD(1:2*NPION)=0; %张成特定大小的0向量for i=1:NFORCEASLOD((FORCE(i,1)*2-1):FORCE(i,1)*2)=FORCE(i,2:3);end%将有约束处的荷载置为0NUM=1;i=0;tmp(NVFIX)=0;while i<NVFIXfor j=-1:0if FIXED(NUM,j+3)==1i=i+1;tmp(i)=FIXED(NUM)*2+j;endendNUM=NUM+1;endfor i=1:NVFIXASLOD(tmp(i))=0;endASDISP=HK\ASLOD'; %计算结点位移向量%所引用的全局变量%global NELEM NPION SMATX ASDISP LNODSELEDISP(1:6)=0; %当前单元的结点位移向量for i=1:NELEMfor j=1:3ELEDISP(j*2-1:j*2)=ASDISP(LNODS(i,j)*2-1:LNODS(i,j)*2);endSTRESS=SMATX*ELEDISP'; %求内力endTHANKS !!!致力为企业和个人提供合同协议，策划案计划书，学习课件等等打造全网一站式需求欢迎您的下载，资料仅供参考。

平面三角形３结点有限元程序1、程序名：，2、程序功能本程序能计算弹性力学的平面应力问题和平面应变问题；能考虑自重和结点集中力两种荷载的作用，在计算自重时y轴取垂直向上为正；能处理非零已知位移，如支座沉降的作用。

主要输出的内容包括：结点位移、单元应力、主应力、第一主应力与x轴的夹角以及约束结点的支座反力。

程序采用Visual Fortran 编制而成，输入数据全部采用自由格式。

3、程序流程及框图图１-1 程序流程图图１-2 程序框图其中，各子程序的功能如下：INPUT——输入结点坐标、单元信息和材料参数；MR——形成结点自由度序号矩阵；FORMMA——形成指标矩阵MA（N）并调用其他功能子程序，相当于主控程序；DIV——取出单元的3个结点号码和该单元的材料号并计算单元的b i，c i等；MGK——形成整体劲度矩阵并按一维存放在SK（NH）中；LOAD——形成整体结点荷载列阵F；OUTPUT——输出结点位移或结点荷载；TREAT——由于有非零已知位移，对K和F进行处理；DECOMP——整体劲度矩阵的分解运算；FOBA——前代、回代求出未知结点位移δ；ERFAC——计算约束结点的支座反力；KRS——计算单元劲度矩阵中的子块K rs。

4、输入数据及变量说明当程序开始运行时，按屏幕提示，键入数据文件的名字。

在运行程序之前，必须根据程序中输入要求建立一个存放原始数据的文件，这个文件的名字由少于8个字符或数字组成。

数据文件包括如下内容：⑴总控信息，共一条，9个数据NP，NE，NM，NR，NI，NL，NG，ND，NCNP——结点总数；NE——单元总数；NM——材料类型总数；NR——约束结点总数；NI ——问题类型标识，0为平面应力问题，1为平面应变问题；NL ——受荷载作用的结点的数目；NG ——考虑自重作用为1，不计自重为0；ND ——非零已知位移结点的数目；NC ——要计算支座约束反力的结点数目。

⑵ 材料信息，共NM 条，每条依次输入EO ，VO ，W ，tEO ——弹性模量（kN/m 2）；VO ——泊松比；W ——材料容重（kN/m 3）；t ——单元厚度（m ）。

matlab编译平面有限元计算编译平面有限元计算是一种常用的数值计算方法，可以用于求解各种复杂的工程问题。

在本文中，我们将介绍如何使用MATLAB编写平面有限元计算程序，并通过一个实例来说明其应用。

让我们了解一下有限元方法的基本原理。

有限元方法是一种将连续体划分为有限个单元，根据物理方程和边界条件，在每个单元上建立离散方程组，最终求解得到整个连续体的近似解的方法。

在平面问题中，连续体被划分为三角形或四边形单元。

MATLAB是一种功能强大的数值计算软件，它提供了丰富的工具箱和函数，可以方便地进行有限元计算。

下面我们将以求解平面弹性力学问题为例，介绍如何使用MATLAB编写平面有限元计算程序。

我们需要定义问题的几何信息、边界条件和材料参数。

在MATLAB中，可以通过定义节点坐标、单元连接关系、边界条件和材料参数来描述问题。

节点坐标可以用一个矩阵表示，其中每一行代表一个节点的坐标。

单元连接关系可以用一个矩阵表示，其中每一行代表一个单元的节点编号。

边界条件可以用一个向量表示，其中每个元素代表一个节点的边界条件。

材料参数可以用一个矩阵表示，其中每一行代表一个单元的材料参数。

接下来，我们需要建立有限元离散方程组。

在平面弹性力学问题中，离散方程组可以通过组装单元刚度矩阵和外力向量得到。

单元刚度矩阵可以通过单元的几何信息和材料参数计算得到。

外力向量可以根据边界条件和载荷信息计算得到。

最终，我们可以通过求解离散方程组得到节点的位移解。

我们可以根据节点的位移解计算出单元的应力和应变。

在平面弹性力学问题中，应力可以通过单元的几何信息和位移解计算得到。

应变可以通过应力和材料参数计算得到。

通过计算应力和应变，我们可以评估结构的稳定性和性能。

MATLAB编译平面有限元计算程序是一种强大的工具，可以用于求解各种复杂的工程问题。

通过定义几何信息、边界条件和材料参数，建立有限元离散方程组，求解得到节点的位移解，最终计算出单元的应力和应变。

有限元课程设计matlab一、课程目标知识目标：1. 学生能理解有限元分析的基本原理，掌握运用MATLAB进行有限元建模和求解的基本步骤。

2. 学生能够运用MATLAB软件进行简单物理场的有限元模拟，并解释模拟结果。

3. 学生掌握如何将实际问题抽象为有限元模型，并能够运用MATLAB进行模型参数的设定和调整。

技能目标：1. 学生能够独立操作MATLAB软件，进行有限元模型的构建和求解。

2. 学生能够通过MATLAB编程实现有限元模型的自动化处理，包括前处理、求解和后处理。

3. 学生通过解决实际问题，提高数值分析能力和计算机应用能力。

情感态度价值观目标：1. 学生培养对科学研究的兴趣，特别是在工程计算和仿真领域。

2. 学生通过解决实际问题，体会数学和工程结合的美，增强对工程问题的探究欲望。

3. 学生通过团队合作解决问题，培养协作精神和解决问题的能力。

本课程针对高年级本科生或研究生，他们具备一定的数学基础和编程能力。

课程性质偏重实践，旨在通过MATLAB这一工具将有限元理论应用于具体问题的求解。

课程目标旨在使学生不仅掌握理论知识，而且能够实际操作，将理论知识转化为解决实际问题的技能。

通过课程学习，学生应能够将所学知识应用于未来的学术研究或工程实践中。

二、教学内容1. 有限元方法基本原理回顾：包括有限元离散化、单元划分、形函数、刚度矩阵和载荷向量等概念。

- 教材章节：第二章 有限元方法基础2. MATLAB编程基础：介绍MATLAB的基本操作、数据结构、流程控制、函数编写等。

- 教材章节：第三章 MATLAB编程基础3. MATLAB中的有限元工具箱使用：学习如何使用MATLAB内置的有限元工具箱进行建模和求解。

- 教材章节：第四章 MATLAB有限元工具箱介绍4. 有限元模型构建与求解：结合实际问题，学习如何构建有限元模型，并进行求解。

- 教材章节：第五章 有限元模型构建与求解5. 实例分析与上机操作：通过案例分析，让学生实际操作MATLAB软件，解决具体的有限元问题。