平面三角形单元常应变单元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[结论与体会]:通过本次的课程设计，我对有限元的概念有了更加深刻的理解，同时也弥补了平时学习是疏忽的地方，充实了有限元知识。

基于MTLAB的3节点三角形单元求解平面弹性问题一．引言本程序相比其他参考资料的程序，在操作方面得到大大简化，只需要一次输入相关力学参量、单元总数、节点总数、节点坐标、单元节点编号，利用此程序就能够很快得到所有单元的单元刚度矩阵，以及系统的整体刚度矩阵。

为了便于后面计算单元的应力，在计算单元刚度矩阵时，MATLAB程序输出了弹性系数矩阵D和几何矩阵B。

此外，本程序在引入位移边界条件求解刚度方程时，采用了化零置1法，避免了手工化简整体刚度矩阵的麻烦。

当然，这在一定程度上加重了计算机的负担，增加了计算机的计算时间，但是这个延长的时间相对于手工化简刚度矩阵所需的时间来说，是微不足道的。

程序需要提前输入的参数有：（1）节点坐标node_gcoord（2）单元节点编号element_node_number （3）节点总数node_number （4）单元总数element_number （5）弹性模量E （6）泊松比NU （7）厚度t（8）平面问题性质指示参数ID，其中1为平面应力问题，2为平面应变问题该程序主要包括以下部分：1. 单元刚度矩阵程序：用于计算各单元的单元刚度矩阵；2. 整体刚度矩阵程序：用于组装整体刚度矩阵；3. 位移及支反力求解程序：用于求解未知位移和未知支反力4. 单元应力计算程序：用于计算各单元的单元应力二．问题描述为了验证程序的正确性，我们选取了参考文献[1]中的一个简单实例进行验证。

问题如下：如图所示为一矩形薄平板，在右端部受集中力F=100kN的作用，材料常数为：弹性模量E=1e7Pa，泊松比NU=1/3，板的厚度t=0.1m。

试按平面应力问题计算各节点位移及支座反力。

三．问题求解（1）结构的离散化与编号第1页共9页1）节点描述 2) 单元描述节点编号节点坐标x节点坐标y121节点i220单元编号3011240023（2）计算单元刚度矩阵KE(:,:,1) =1.0e+006 *0.2813 0 0 0.1875 -0.2813 -0.1875 00.0938 0.1875 0 -0.1875 -0.0938 0 0.1875 0.37500 -0.3750 -0.1875 0.1875 0 0 1.1250 -0.1875 -1.1250 -0.2813 -0.1875 -0.3750 -0.1875 0.6563 0.3750 -0.1875-0.0938 -0.1875 -1.1250 0.3750 1.2188KE(:,:,2) =1.0e+006 *0.2813 0 0 0.1875 -0.2813 -0.1875 00.0938 0.1875 0 -0.1875 -0.0938 0 0.1875 0.37500 -0.3750 -0.1875 0.1875 0 0 1.1250 -0.1875 -1.1250 -0.2813 -0.1875 -0.3750 -0.1875 0.6563 0.3750 -0.1875-0.0938 -0.1875 -1.1250 0.3750 1.2188 （3）组装整体刚度矩阵KK =1.0e+006 *0.6563 0.3750 -0.3750 -0.1875 -0.2813 -0.1875 00 0.3750 1.2188 -0.1875 -1.1250 -0.1875 -0.0938 00 -0.3750 -0.1875 0.6563 0 0 0.3750 -0.2813 -0.1875 -0.1875 -1.1250 0 1.2188 0.3750 0 -0.1875-0.0938 -0.2813 -0.1875 0 0.3750 0.6563 0 -0.3750-0.1875 -0.1875 -0.0938 0.3750 0 0 1.2188 -0.1875-1.1250 0 0 -0.2813 -0.1875 -0.3750 -0.1875 0.65630.3750 0 0 -0.1875 -0.0938 -0.1875 -1.1250 0.37501.2188第2页共9页节点j节点m3421 （4）求解位移及支反力 delta =0.1877 -0.8992 -0.1497 -0.8422 0 0 0 0 F =1.0e+004 *0 -0.5000 0 -0.5000 -2.0000 -0.0702 2.00001.0702（5）计算单元应力 stress(:,:,1) =1.0e+005 *-0.8419 -0.2806 -1.5791stress(:,:,2) =1.0e+004 *8.4187 -2.8953 -4.2094第3页共9页四．结果分析可以看出，计算得到的单元1的应力分量为σx =-84190Pa,σy=-28060Pa,τxy=-157910Pa,单元2的应力分量为σx =84187Pa,σy=-28953Pa,τxy=-42094Pa.此结果与参考文献[1]上给出的结果完全吻合。

计算力学(有限元方法部分) 程序设计程序说明书程序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;最后，我们可以使用这些信息来求解三角形单元的位移解。

%变量说明%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 !!!致力为企业和个人提供合同协议，策划案计划书，学习课件等等打造全网一站式需求欢迎您的下载，资料仅供参考。

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

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

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

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

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

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

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

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

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

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

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

8）计算结果整理。

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

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

●基本思想：单元结点按右手法则顺序编号。

●荷载类型：可计算结点荷载。

●说明：主程序的作用是通过赋值语句、读取和写入文件、函数调用等完成算法的全过程，即实现程序流程图的程序表达。

1 程序准备format short e %设定输出类型clear all %清除所有已定义变量clc %清屏●说明：format short e －设定计算过程中显示在屏幕上的数字类型为短格式、科学计数法；clear all －清除所有已定义变量，目的是在本程序的运行过程中，不会发生变量名相同等可能使计算出错的情况；clc －清屏，使屏幕在本程序运行开始时2 全局变量定义global NNODE NPION NELEM NVFIX NFORCE COORD LNODS YOUNG POISS THICKglobal FORCE FIXEDglobal BMATX DMATX SMATX AREAglobal ASTIF ASLOD ASDISPglobal FP1说明：NNODE—单元结点数，NPION—总结点数， NELEM—单元数，NVFIX—受约束边界点数，NFORCE—结点力数，COORD—结构结点坐标数组，LNODS —单元定义数组，YOUNG—弹性模量，POISS—泊松比，THICK—厚度FORCE —节点力数组(n,3) n:受力节点数目,(n,1):作用点,(n,2):x方向,(n,3):y方向； FIXED—约束信息数组(n,3) n:受约束节点数目, (n,1):约束点 (n,2)与(n,3)分别为约束点x方向和y方向的约束情况,受约束为1否则为0BMATX—单元应变矩阵(3*6)， DMATX—单元弹性矩阵(3*3)，SMATX—单元应力矩阵(3*6)，AREA—单元面积ASTIF—总体刚度矩阵，ASLOD—总体荷载向量，ASDISP—结点位移向量FP1—数据文件指针3 打开文件FP1=fopen('input.txt','rt'); %打开输入数据文件存放初始数据●说明：FP1=fopen('input.txt','rt'); －打开已存在的输入数据文件input.txt，且设置其为只读格式，使程序在执行过程中不能改变输入文件中的数值，并用文件句柄FP1来执行FP2=fopen('output.txt','wt'); －打开输出数据文件，该文件不存在时，通过此命令创建新文件，该文件存在时则将原有内容全部删除。

用常应变三角形单元解弹性力学平面问题的程序******************************************************************** ANALYSIS PROGTAM OF FINITE ELEMENT METHOD ** FOR PLANE STRESS/STRAIN OF TRIANGULAR ELEMENT ** ----- FEMT3.FOR ----- **------------------------------------------------------------- ** Subroutines: 1-SDATA, 2-STE, 3-ATE, 4-DTE, 5-BTE, 6-STIFF ** 7-EQUPE, 8-INSCD, 9-BGSMT, 10-SIGME ********************************************************************DIMENSION LND(50,3),X(100),Y(100),JR(20,3),PJ(20,3),P(200)REAL KS(200,100)OPEN(5,FILE='FEMT3.DAT')OPEN(6,FILE='FEMT3.OUT',STATUS='NEW')READ(5,*) NJ,NE,NS,NPJ,IPS(结点、单元、支承、荷载、类型)WRITE(6,*)' FINITE ELEMENT ANALYSIS IN PLANE PROBLEM'WRITE(6,*)' SOURCE DATA OUTPUT'WRITE(6,20) NJ,NE,NS,NPJ,IPS20 FORMAT(4X,'NJ',3X,'NE',3X,'NS',3X,'NPJ',2X,'IPS'/1X,5I5)IF(IPS.EQ.0) WRITE(6,*)' PLANE STRESS PROBLEM'IF(IPS.EQ.1) WRITE(6,*)' PLANE STRAIN PROBLEM'CALL SDATA(NJ,NE,NS,NW,NPJ,IPS,E,PR,T,V,LND,X,Y,JR,PJ)NJ2=2*NJWRITE(6,50) NJ250 FORMAT(/1X,'DEGREES OF FREEDOM=',I5)WRITE(6,60) NW60 FORMAT(1X,'BAND WIDTH=',I5)CALL STIFF(NJ,NE,NJ2,NW,LND,X,Y,E,PR,T,KS)(总刚6)CALL EQUPE(NJ,NE,NPJ,NJ2,T,V,LND,X,Y,PJ,P)({P}7)CALL INSCD(NS,NW,NJ2,JR,KS,P)(引入支承条件8)CALL BGSMT(NJ,NJ2,NW,KS,P)（解方程9）CALL SIGME(NE,NJ,NJ2,E,PR,LND,X,Y,P)（求应力10）CLOSE(5)CLOSE(6)END*--------------------------------------------------------C SUBPROGRAM-1C INPUT STRUCTURAL DATASUBROUTINE SDATA(NJ,NE,NS,NW,NPJ,IPS,E,PR,* T,V,LND,X,Y,JR,PJ)DIMENSION LND(NE,3),X(NJ),Y(NJ),JR(NS,3),PJ(NPJ,3)READ(5,*) E,PR,T,V(弹性模量、泊松比、厚度、容重)WRITE(6,10) E,PR,T,V10 FORMAT(/6X,'E',10X,'PR',9X,'T',9X,'V'/,4F10.2)READ(5,*)((LND(I,J),J=1,3),I=1,NE)（结点编码）WRITE(6,20)20 FORMAT(/1X,'ELEMENT INFORMATION'/3X,'ELEM',3X,* 'I J K'/)WRITE(6,30)(I,(LND(I,J),J=1,3),I=1,NE)30 FORMAT(1X,4I5)READ(5,*)(X(I),Y(I),I=1,NJ)（结点坐标）WRITE(6,40)40 FORMAT(/1X,'COORDINATES OF NODES'/3X,'NODES',* 8X,'X',13X,'Y')WRITE(6,50)(I,X(I),Y(I),I=1,NJ)50 FORMAT(1X,I5,2E15.6)READ(5,*)((JR(I,J),J=1,3),I=1,NS)（约束信息）WRITE(6,60)60 FORMAT(/1X,'CONSTRAINED NODES'/3X,'NODE',3X,'X',4X,'Y') WRITE(6,70)((JR(I,J),J=1,3),I=1,NS)70 FORMAT(1X,3I5)READ(5,*)((PJ(I,J),J=1,3),I=1,NPJ)（荷载信息）WRITE(6,80)80 FORMAT(/1X,'LOAD CASES'/3X,'NODE',8X,'X',13X,'Y')WRITE(6,90)((PJ(I,J),J=1,3),I=1,NPJ)90 FORMAT(1X,F5.0,2E15.6)100 NW=0（半带宽）DO 110 IE=1,NEDO 110 I=1,3DO 110 J=1,3IW=IABS(LND(IE,I)-LND(IE,J))IF(NW.LT.IW) THENNW=IWENDIF110 CONTINUENW=(NW+1)*2IF(IPS.NE.0) THENE=E/(1.0-PR*PR)PR=PR/(1.0-PR)ENDIFEND*---------------------------------------------------------C SUBPROGRAM-2C CALCULATE ELEMENT STIFFNESS MATRIXSUBROUTINE STE(IE,NJ,NE,LND,X,Y,E,PR,T,KE)DIMENSION LND(NE,3),X(NJ),Y(NJ),B(3,6),D(3,3)REAL KE(6,6)CALL ATE(IE,NJ,NE,LND,X,Y,AE)CALL DTE(E,PR,D)CALL BTE(IE,NJ,NE,LND,X,Y,AE,B)DO 10 I=1,6DO 10 J=1,6KE(I,J)=0.DO 10 K=1,3DO 10 K1=1,310 KE(I,J)=KE(I,J)+B(K,I)*D(K,K1)*B(K1,J)C=AE*TDO 30 I=1,6DO 30 J=1,630 KE(I,J)=KE(I,J)*CEND*------------------------------------------------ C SUBPROGRAM-3C CALCULATE ELEMENT AREASUBROUTINE ATE(IE,NJ,NE,LND,X,Y,AE)DIMENSION LND(NE,3),X(NJ),Y(NJ)I=LND(IE,1)J=LND(IE,2)K=LND(IE,3)XIJ=X(J)-X(I)YIJ=Y(J)-Y(I)XIK=X(K)-X(I)YIK=Y(K)-Y(I)AE=.5*(XIJ*YIK-XIK*YIJ)END*---------------------------------------------- C SUBPROGRAM-4C CALCULATE ELASTICITY MATRIXSUBROUTINE DTE(E,PR,D)DIMENSION D(3,3)DO 10 I=1,3DO 10 J=1,310 D(I,J)=0.D(1,1)=E/(1.-PR*PR)D(1,2)=E*PR/(1.-PR*PR)D(2,1)=D(1,2)D(2,2)=D(1,1)D(3,3)=.5*E/(1.+PR)END*------------------------------------------------ C SUBPROGRAM-5C CALCULATE MATRIX [B]SUBROUTINE BTE(IE,NJ,NE,LND,X,Y,AE,B)DIMENSION LND(NE,3),X(NJ),Y(NJ),B(3,6)I=LND(IE,1)J=LND(IE,2)K=LND(IE,3)DO 10 II=1,3DO 10 JJ=1,610 B(II,JJ)=0.B(1,1)=Y(J)-Y(K)B(1,3)=Y(K)-Y(I)B(1,5)=Y(I)-Y(J)B(2,2)=X(K)-X(J)B(2,4)=X(I)-X(K)B(2,6)=X(J)-X(I)B(3,1)=B(2,2)B(3,2)=B(1,1)B(3,3)=B(2,4)B(3,4)=B(1,3)B(3,5)=B(2,6)B(3,6)=B(1,5)DO 20 I1=1,3DO 20 J1=1,620 B(I1,J1)=.5/AE*B(I1,J1)END*------------------------------------------------------- C SUBPROGRAM-6C CALCULATE GLOBAL STIFFNESS MATRIXSUBROUTINE STIFF(NJ,NE,NJ2,NW,LND,X,Y,E,PR,T,KS) DIMENSION LND(NE,3),X(NJ),Y(NJ)REAL KS(NJ2,NW),KE(6,6)DO 5 I=1,NJ2DO 5 J=1,NW5 KS(I,J)=0.DO 10 IE=1,NECALL STE(IE,NJ,NE,LND,X,Y,E,PR,T,KE)DO 10 I=1,3IZ=LND(IE,I)DO 10 II=1,2IH =2*(I -1)+IIIDH=2*(IZ-1)+IIDO 10 J=1,3JZ=LND(IE,J)DO 10 JJ=1,2L =2*(J -1)+JJIL=2*(JZ-1)+JJIF(IL.GE.IDH) THENIDL=IL-IDH+1KS(IDH,IDL)=KS(IDH,IDL)+KE(IH,L)ENDIF10 CONTINUEEND*-------------------------------------------------------- C SUBPROGRAM-7C CALCULATE NODAL LOAD VECTORSUBROUTINE EQUPE(NJ,NE,NPJ,NJ2,T,V,LND,X,Y,PJ,P) DIMENSION LND(NE,3),X(NJ),Y(NJ),PJ(NPJ,3),P(NJ2) DO 10 I=1,NJ210 P(I)=0.DO 20 I=1,NPJII=PJ(I,1)P(2*II-1)=PJ(I,2)20 P(2*II)=PJ(I,3)30 IF(V.EQ.0.) GOTO 50DO 40 IE=1,NECALL ATE(IE,NJ,NE,LND,X,Y,AE)PE=-V*AE*T/3.DO 40 I=1,3II=LND(IE,I)40 P(2*II)=P(2*II)+PE50 RETURNEND*---------------------------------------------C SUBPROGRAM-8C INTRODUCE BOUNDARY CONDITIONSUBROUTINE INSCD(NS,NW,NJ2,JR,KS,P)DIMENSION P(NJ2),JR(NS,3)REAL KS(NJ2,NW)DO 30 I=1,NSIR=JR(I,1)DO 30 J=2,3IF(JR(I,J).EQ.0) GOTO 30II=2*IR+J-3KS(II,1)=1.DO 10 JJ=2,NW10 KS(II,JJ)=0.IF(II.GT.NW) JO=NWIF(II.LE.NW) JO=IIDO 20 JJ=2,JO20 KS(II-JJ+1,JJ)=0.P(II)=0.30 CONTINUEEND*-------------------------------------------C SUBPROGRAM-9C SOLVE EQUATIONS BY GS METHODSUBROUTINE BGSMT(NJ,NJ2,NW,KS,P)DIMENSION P(NJ2)REAL KS(NJ2,NW)NJ1=NJ2-1DO 20 K=1,NJ1IF(NJ2.GT.K+NW-1) IM=K+NW-1IF(NJ2.LE.K+NW-1) IM=NJ2K1=K+1DO 20 I=K1,IML=I-K+1C=KS(K,L)/KS(K,1)IW=NW-L+1DO 10 J=1,IWM=J+I-K10 KS(I,J)=KS(I,J)-C*KS(K,M)20 P(I)=P(I)-C*P(K)P(NJ2)=P(NJ2)/KS(NJ2,1)DO 40 I1=1,NJ1I=NJ2-I1IF(NW.GT.NJ2-I+1) JO=NJ2-I+1IF(NW.LE.NJ2-I+1) JO=NWDO 30 J=2,JOK=J+I-130 P(I)=P(I)-KS(I,J)*P(K)40 P(I)=P(I)/KS(I,1)WRITE(6,50)50 FORMAT(/1X,'NODAL DISPLACEMENTS'/3X,* 'NODE',5X,'X-DISP.',8X,'Y-DISP.')DO 60 I=1,NJ60 WRITE(6,70) I,P(2*I-1),P(2*I)70 FORMAT(1X,I5,2E15.6)END*--------------------------------------------------- C SUBPROGRAM-10C CALCULATE ELEMENT STRESS MATRIXSUBROUTINE SIGME(NE,NJ,NJ2,E,PR,LND,X,Y,P)DIMENSION LND(NE,3),X(NJ),Y(NJ),D(3,3),B(3,6), * S(3,6),ST(3),P(NJ2),DE(6)WRITE(6,*)WRITE(6,*)' ELEMENT STRESSES'CALL DTE(E,PR,D)DO 50 IE=1,NECALL ATE(IE,NJ,NE,LND,X,Y,AE)CALL BTE(IE,NJ,NE,LND,X,Y,AE,B)DO 10 I=1,3DO 10 J=1,6S(I,J)=0.DO 10 K=1,310 S(I,J)=S(I,J)+D(I,K)*B(K,J)DO 20 I=1,3DO 20 J=1,2IH=2*(I-1)+JIW=2*(LND(IE,I)-1)+J20 DE(IH)=P(IW)DO 30 I=1,3ST(I)=0.DO 30 J=1,630 ST(I)=ST(I)+S(I,J)*DE(J)SGX=ST(1)SGY=ST(2)TXY=ST(3)ASG=(SGX+SGY)*.5RSG=SQRT(.25*(SGX-SGY)**2+TXY*TXY)SGMA=ASG+RSGSGMI=ASG-RSGIF(SGY.EQ.SGMI) CETA=0.IF(SGY.NE.SGMI) CETA=90.-57.29578*ATAN* (TXY/(SGY-SGMI))50 WRITE(6,60) IE,SGX,SGY,TXY,SGMA,SGMI,CETA60 FORMAT(1X,'ELEMENT NO.=',I4/2X,'SIGX=',E10.4, * 2X,'SIGY=',E10.4,2X,'TXY =',E10.4/2X,'SGMA=', * E10.4,2X,'SGMI=',E10.4,2X,'CETA=',E10.4)END。

基于Matlab语言的按平面三角形单位划分的结构有限元法式设计之阿布丰王创作专业：建筑与土木工程班级：建工研12-2姓名：韩志强学号：471220580基于Matlab语言的按平面三角形单位划分结构有限元法式设计一、有限单位发及Matlab语言概述1. 有限单位法随着现代工业、生产技术的发展,不竭要求设计高质量、高水平的年夜型、复杂和精密的机械及工程结构.为此目的,人们必需预先通过有效的计算手段,确切的预测即将出生的机械和工程结构,在未来工作时所发生的应力、应变和位移因此,需要寻求一种简单而又精确的数值分析方法.有限单位法正是适应这种要求而发生和发展起来的一种十分有效的数值计算方法.有限元法把一个复杂的结构分解成相对简单的“单位”,各单位之间通过结点相互连接.单位内的物理量由单位结点上的物理量按一定的假设内插获得,这样就把一个复杂结构从无限多个自由度简化为有限个单位组成的结构.我们只要分析每个单位的力学特性,然后依照有限元法的规则把这些单位“拼装”成整体,就能够获得整体结构的力学特性.有限单位法基本步伐如下：(1)结构离散：结构离散就是建立结构的有限元模型,又称为网格划分或单位划分,即将结构离散为由有限个单位组成的有限元模型.在该步伐中,需要根据结构的几何特性、载荷情况等确订单位体内任意一点的位移插值函数.(2)单位分析：根据弹性力学的几何方程以及物理方程确订单位的刚度矩阵.(3)整体分析:把各个单位按原来的结构重新连接起来,并在单位刚度矩阵的基础上确定结构的总刚度矩阵,形成如下式所示的整体有限元线性方程：式中阵.(4)载荷移置:根据静力等效原理,将载荷移置到相应的节点上,形成节点载荷矩阵.(5)鸿沟条件处置:对式①所示的有限元线性方程进行鸿沟条件处置.(6)求解线性方程:求解式①所示的有限元线性方程,获得节点的位移.在该步伐中,若有限元模型的节点越多,则线性方程的数量就越多,随之有限元分析的计算量也将越年夜.(7)求解单位应力及应变根据求出的节点位移求解单位的应力和应变.(8)结果处置与显示进入有限元分析的后处置部份,对计算出来的结果进行加工处置,并以各种形式将计算结果显示出.2.Matlab简介在用有限元法进行结构分析时,将会遇到年夜量的数值计算,因而在实用上是一定要借助于计算机和有限元法式,才华完成这些复杂而繁重的数值计算工作.而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 Array单位节点编码数组节点坐标数组节点力数组FORCE =[4 0 -150]约束信息数组以上数值数据为法式运行前输入的初始数据,存为“471220580.txt”文本格式,同时必需放在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('471220580.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])' %结点坐标数组%坐标：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 =2.0000e-001THICK =2.0000e-003LNODS =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 =5.0000e-001D =2.0833e+008 4.1667e+007 0 4.1667e+007 2.0833e+008 00 0 8.3333e+007A =5.0000e-001D =2.0833e+008 4.1667e+007 0 4.1667e+007 2.0833e+008 0 0 0 8.3333e+007A =5.0000e-001D =2.0833e+008 4.1667e+007 04.1667e+007 2.0833e+008 00 0 8.3333e+007A =5.0000e-001ASDISP =-8.0607e-004-1.5848e-003-1.0281e-003-4.4727e-0031.1937e-003-4.6947e-0038.6670e-004-1.8880e-003（说明：以上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.00002 -0.80607E-03-0.15848E-02 0.0000 0.17780E-023 -0.10281E-02-0.44727E-02 0.0000 0.45893E-024 0.11937E-02-0.46947E-02 0.0000 0.48441E-025 0.86670E-03-0.18880E-02 0.0000 0.20774E-026 0.0000 0.0000 0.0000 0.0000MAXIMUM ABSOLUTE VALUESNODE 4 4 0 4VALUE 0.11937E-02-0.46947E-02 0.0000 0.48441E-02ANSYS求解单位应力结果列表显示如下：（单位：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.0000ELEMENT= 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'end。

第三章平面问题的有限元法本章通过三角形常应变单元，介绍有限元法应用于弹性体应力分析的基本原理和方法：包括弹性体的离散化，单元特性的分析，刚度矩阵的建立，等效节点力的计算，解答的收敛性以及实施步骤和注意事项，同时对形函数的性质也作简要的叙述。

第一节三角形常应变单元一、结构离散化用有限元法分析弹性力学平面问题，第一步就是把原来的连续的弹性体离散化。

(a) (b)将整个结构（平板）划分成有限个三角形。

这样的三角形称为单元（三角形单元）。

三角形单元的顶点取为节点。

3节点三角形单元用边界节点间的直线段来近似板的曲线边界。

这些三角形在其节点处相互连接，组成一个单元集合体，以代替原来的弹性体。

注：1. 全部节点和全部单元一般由1开始按自然顺序编号。

2. 节点编码：总码-----------用于整体分析，如1,2，…，按自然顺序编号局部码--------用于单元分析，i，j，m 要求按逆时针方向的顺序进行编码每个单元的节点局部码i，j，m和节点总码有一一对应的关系3. 单元间不能有重叠4. 一个单元的任一顶点不许为另一单元任一边的内点5. 所有作用在单元上的载荷，包括集中载荷、表面载荷和体积力，都按虚功等效的原则移置到节点上，成为等效节点载荷。

二、位移模式1. 单元节点位移列阵iu设单元e 的节点号码为i ，j ，m 。

由弹性力学平面问题可知，单元内任意一点有两个位移分量u ，v ，记为{}Tf uv ⎡⎤⎣⎦=故每个节点也有两个位移分量，因此称节点自由度为2。

节点i 位移分量记为{}Ti ii u v δ⎡⎤⎣⎦= （i ，j ，m 轮换）则3个节点的位移分量用列阵表示为{}i ei i e j j j m m m u v u v u v δδδδ⎧⎫⎪⎪⎪⎪⎧⎫⎪⎪⎪⎪⎪⎪⎪⎪⎨⎬⎨⎬⎪⎪⎪⎪⎪⎪⎪⎪⎩⎭⎪⎪⎪⎪⎩⎭==（3-1）称为单元节点位移列阵（向量）。

单元自由度是6。

2. 位移模式结构承受载荷以后发生位移，但位移分布事先并不知道。

matlab三角形单元扩展编程在MATLAB 中，扩展三角形单元的编程通常涉及到对现有的三角形单元进行修改或添加新的功能。

下面是一个简单的例子，演示如何创建一个用于扩展三角形单元的MATLAB 类。

假设你想要创建一个扩展的三角形单元类，其中包含了额外的属性和方法。

以下是一个可能的实现示例：```matlabclassdef ExtendedTriangle < handlepropertiesNodes % 顶点坐标Area % 三角形面积endmethods% 构造函数function obj = ExtendedTriangle(nodes)if nargin > 0obj.Nodes = nodes;obj.calculateArea();endend% 计算三角形面积function calculateArea(obj)% 在这里实现计算三角形面积的逻辑% 这可能涉及到对顶点坐标的操作% 这里仅作为示例，假设nodes 是一个3x2 的矩阵obj.Area = polyarea(obj.Nodes(:,1), obj.Nodes(:,2));end% 自定义的其他方法可以在这里添加endend```在这个示例中，`ExtendedTriangle` 类具有`Nodes` 属性表示三角形的顶点坐标，并且有一个`Area` 属性表示三角形的面积。

构造函数负责初始化这些属性，而`calculateArea` 方法用于计算三角形的面积。

你可以根据实际需求添加更多的属性和方法。

使用这个类的示例：```matlab% 创建一个三角形对象triangleNodes = [0, 0; 1, 0; 0, 1];triangle = ExtendedTriangle(triangleNodes);% 访问属性disp('Triangle Nodes:');disp(triangle.Nodes);disp(['Triangle Area: ', num2str(triangle.Area)]);```这只是一个简单的例子，你可以根据实际需要扩展和修改这个类。

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

有限元剖析的基本步骤可概括为三大步：构造失散、单元剖析和整体剖析。

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

开始Matlab 语言是进行矩阵运算的强盛工具，因此，用 Matlab 语言编写有限元中平面问题的程序输入初始数据有优胜性。

本章将详尽介绍怎样利用Matlab 语言编制三角形常应变单元的计算程序，程序流程图见生成单刚集成总刚图 1。

有限元法中三节点三角形剖析构造的步骤如施加拘束信息下：1）整理原始数据，如资料性质、荷载条件、约束条件等，失散构造并进行单元编码、结点生成荷载向量编码、结点位移编码、选用坐标系。

2）单元剖析，成立单元刚度矩阵。

界限条件办理3）整体剖析，成立总刚矩阵。

4）成立整体构造的等效节点荷载和总荷载矩计算结点位移阵5）界限条件办理。

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

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

8）计算结果整理。

计算结果整理包含位移和应计算结果整理力两个方面；位移计算结果一般不需要特别的办理，利用计算出的节点位移重量，便可结束画出构造随意方向的位移云图；而应力解的偏差表此刻单元内部不知足均衡方程，单元图 1 程序流程图与单元界限处应力一般不连续，在界限上应1.1 程序说明%*******************************************************************%三角形常应变单元求解构造主程序%*******************************************************************功能：运用有限元法中三角形常应变单元解平面问题的计算主程序。

《有限元法》实验报告专业班级力学（实验）1601 姓名田诗豪学号 10提交日期实验一（30分）一、实验内容编写一个计算平面3结点三角形单元的应变矩阵、应力矩阵和单元刚度矩阵的MATLAB 函数文件[B3,S3,K3] = ele_mat_tri3(xy3,mat)，其中：输入变量xy3为结点坐标数组，mat 为材料参数矩阵；输出变量B3为应变矩阵，S3为应力矩阵，K3为单元刚度矩阵。

（要求给出3个不同算例进行验证，并绘制出单元形状和结点号）二、程序代码通用函数function [B3,S3,K3] = ele_mat_tri3(xy3,mat)%生成平面3结点三角形单元的应变矩阵、应力矩阵和单元刚度矩阵的功能函数%*********变量说明****************%xy3------------------结点坐标数组%mat------------------材料参数矩阵（弹性模量，泊松比，壁厚）%B3-------------------应变矩阵%S3-------------------应力矩阵%K3-------------------单元刚度矩阵%*********************************xyh=[1,xy3(1,1),xy3(1,2);1,xy3(2,1),xy3(2,2);1,xy3(3,1),xy3(3,2)];A=*det(xyh);A=abs(A);D=mat(1)/(1-mat(2)^2)*[1,mat(2),0;mat(2),1,0;0,0,(1-mat(2))/2];b=zeros(1,3);c=zeros(1,3);%*********************************for i=1:3if i==1j=2;m=3;elseif i==2j=3;m=1;elsej=1;m=2;endb(i)=xy3(j,2)-xy3(m,2);c(i)=xy3(m,1)-xy3(j,1);end%*********************************B31=1/(2*A)*[b(1),0;0,c(1);c(1),b(1)];B32=1/(2*A)*[b(2),0;0,c(2);c(2),b(2)];B33=1/(2*A)*[b(3),0;0,c(3);c(3),b(3)];B3=[B31,B32,B33];%*********************************S3=D*B3;%*********************************K3=A*mat(3)*B3'*D*B3;主程序clear;clc;%*********输入结点坐标数组********xy3=[0,0;5,1;1,4];mat=[3e6,,];%****输入材料参数矩阵(弹性模量，泊松比，壁厚)****[B3,S3,K3]=ele_mat_tri3(xy3,mat)三、算例分析算例1：如图1所示三角形单元，结点坐标为1(0,0)，2(5,2)，3(1,4)，弹性模量为200GPa，泊松比为、厚度为。

姓名：刘刚学号：15平面应力应变分析有限元法Abstruct:本文通过对平面应力/应变问题的简要理论阐述，使读者对要分析的问题有大致的印象，然后结合两个实例，通过MATLAB软件的计算，将有限元分析平面应力/应变问题的过程形象的展示给读者，让人一目了然，快速了解有限元解决这类问题的方法和步骤！一.基本理论有限元法的基本思路和基本原则以结构力学中的位移法为基础，把复杂的结构或连续体看成有限个单元的组合，各单元彼此在节点出连接而组成整体。

把连续体分成有限个单元和节点，称为离散化。

先对单元进行特性分析，然后根据节点处的平衡和协调条件建立方程，综合后做整体分析。

这样一分一合，先离散再综合的过程，就是把复杂结构或连续体的计算问题转化简单单元分析与综合问题。

因此，一般的有限揭发包括三个主要步骤：离散化单元分析整体分析。

二.用到的函数1. LinearTriangleElementStiffness(E,NU,t,xi,yi,xj,yj,xm,ym,p)2.LinearBarAssemble(K k I f)3.LinearBarElementForces(k u)4.LinearBarElementStresses(k u A)5.LinearTriangleElementArea(E NU t)三.实例例1.考虑如图所示的受均布载荷作用的薄平板结构。

将平板离散化成两个线性三角元，假定E=200GPa，v=0.3，t=0.025m,w=3000kN/m.1.离散化2.写出单元刚度矩阵通过matlab 的LinearTriangleElementStiffness 函数，得到两个单元刚度矩阵1k 和2k ，每个矩阵都是6 6的。

>> E=210e6E =210000000>>k1=LinearTriangleElementStiffness(E,NU,t,0,0,0.5,0.25,0,0.25,1)k1 =1.0e+006 *Columns 1 through 52.0192 0 0 -1.0096 -2.01920 5.7692 -0.8654 0 0.86540 -0.8654 1.4423 0 -1.4423-1.0096 0 0 0.5048 1.0096-2.0192 0.8654 -1.4423 1.0096 3.46151.0096 -5.7692 0.8654 -0.5048 -1.8750Column 61.0096-5.76920.8654-0.5048-1.87506.2740>> NU=0.3NU =0.3000>> t=0.025t =0.0250>> k2=LinearTriangleElementStiffness(E,NU,t,0,0,0.5,0,0.5,0.25,1) k2 =1.0e+006 *Columns 1 through 51.4423 0 -1.4423 0.8654 00 0.5048 1.0096 -0.5048 -1.0096-1.4423 1.0096 3.4615 -1.8750 -2.01920.8654 -0.5048 -1.8750 6.2740 1.00960 -1.0096 -2.0192 1.0096 2.0192-0.8654 0 0.8654 -5.7692 0 Column 6-0.86540.8654-5.76925.76923.集成整体刚度矩阵8*8零矩阵K =0 0 0 0 0 0 0 00 0 0 0 0 0 0 00 0 0 0 0 0 0 00 0 0 0 0 0 0 00 0 0 0 0 0 0 00 0 0 0 0 0 0 00 0 0 0 0 0 0 00 0 0 0 0 0 0 0>> K=LinearTriangleAssemble(K,k1,1,3,4)K =1.0e+006 *Columns 1 through 52.0192 0 0 0 00 5.7692 0 0 -0.86540 0 0 0 00 0 0 0 00 -0.8654 0 0 1.4423-1.0096 0 0 0 0-2.0192 0.8654 0 0 -1.44231.0096 -5.7692 0 0 0.8654Columns 6 through 8-1.0096 -2.0192 1.00960 0.8654 -5.76920 0 00 0 00 -1.4423 0.86540.5048 1.0096 -0.50481.0096 3.4615 -1.8750-0.5048 -1.8750 6.2740>> K=LinearTriangleAssemble(K,k1,1,2,3)K =1.0e+007 *0.4038 0 0 -0.1010 -0.2019 0 -0.2019 0.10100 1.1538 -0.0865 0 0 -0.5769 0.0865 -0.57690 -0.0865 0.1442 0 -0.1442 0.0865 0 0-0.1010 0 0 0.0505 0.1010 -0.0505 0 0-0.2019 0 -0.1442 0.1010 0.4904 -0.1875 -0.1442 0.08650 -0.5769 0.0865 -0.0505 -0.1875 0.6779 0.1010 -0.0505-0.2019 0.0865 0 0 -0.1442 0.1010 0.3462 -0.18750.1010 -0.5769 0 0 0.0865 -0.0505 -0.1875 0.62744.引入边界条件.用上一步得到的整体刚度矩阵，可以得到该结构的方程组如下形式⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡4y 4X 3y 3X 2y 2X 1y 1X 4y 4X 3y 3X 2y 2X 1y 1X 6F F F F F F F F U U U U U U U U 6.2740 1.8750- 0.5048- 0.8654 0 0 5.7692- 1.0096 1.8750- 3.4615 1.0096 1.4423- 0 0 0.8654 2.0192- 0.5048- 1.0096 6.2740 0 5.7692- 0.8654 0 1.8750- 0.8654 1.4423- 0 3.4615 1.0096 2.0192- 1.8750- 0 0 0 5.7692- 1.0096 6.2740 1.8750- 0.5048- 0.8654 0 0 0.8654 2.0192- 1.8750- 3.4615 1.0096 1.4423- 5.7692- 0.8654 0 1.8750- 0.5048- 1.0096 6.2740 0 1.0096 2.0192- 1.8750- 0 0.8654 1.4423- 0 3.461510 本题的边界条件：04411====y x y x U U U U0,375.9,0,375.93322====y x y x F F F F将边界条件带入，得到：⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡4y 4X 1y 1X 3y 3X 2y 2X 6F F 09.37509.375F F 0 0 U U U U 0 0 6.2740 1.8750- 0.5048- 0.8654 0 0 5.7692- 1.0096 1.8750- 3.4615 1.0096 1.4423- 0 0 0.8654 2.0192- 0.5048- 1.0096 6.2740 0 5.7692- 0.8654 0 1.8750- 0.8654 1.4423- 0 3.4615 1.0096 2.0192- 1.8750- 0 0 0 5.7692- 1.0096 6.2740 1.8750- 0.5048- 0.8654 0 0 0.8654 2.0192- 1.8750- 3.4615 1.0096 1.4423- 5.7692- 0.8654 0 1.8750- 0.5048- 1.0096 6.2740 0 1.0096 2.0192- 1.8750- 0 0.8654 1.4423- 0 3.4615105.解方程分解上述方程组，提取总体刚度矩阵K 的第3-6行的第3-6列作为子矩阵⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎣⎡⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡ 09.37509.375 U U U U 6.2740 0 5.7692- 0.8654 0 3.4615 1.0096 2.0192- 5.7692- 1.0096 6.2740 1.8750- 0.8654 2.0192- 1.8750- 3.4615103y 3X 2y 2X 6 Matlab 命令>> k=K(3:6,3:6)k =1.0e+006 *3.4615 -1.8750 -2.0192 0.8654-1.8750 6.2740 1.0096 -5.7692-2.0192 1.0096 3.4615 00.8654 -5.7692 0 6.2740>> f=[9.375;0;9.375;0]f =9.37509.3750>> u=k\fu =1.0e-005 *0.71110.11150.65310.0045现在可以清楚的看出，节点2的水平位移和垂直位移分别是0.7111m和0.1115m。