基于matlab的有限元法分析平面应力应变问题刘刚
有限元法MATLAB作业
例题：如图所示的受力均匀分布载荷作用的薄平板结构，将平板离散成两个线性三角元，如图所示，假定2210,0.3,0.025,3000/m E GPa v t m w KN ====。
求 （1）该结构的整体刚度矩阵。
（2）节点2和节点3的水平位移和垂直位移。
（3）节点1和节点4的支反力。
（4）每个单元的应力。
（5）每个单元的主应力和主应力方向角。
图2. 用双线性三角元离散化的薄木板 解：（1）离散化域我们将平板分为两个单元，4个节点，如图2所示，分布载荷的总作用力平均分给节点2和节点3，由于结构是薄平板，所以假定其属于平面应力情况。
MATLAB 中采用的计算单位是KN 和m 。
表1给出了该题的单元连通性。
（2）单元刚度矩阵通过调用MATLAB 的LinearTriangleElementStiffness 函数，得到两个单元矩阵k 1和k 2，每个矩阵都是66⨯的。
k1=gangdujuzhen(@LinearTriangleElementStiffness, 210000000, 0.3, 0.025, 0, 0, 0.5, 0.25, 0, 0.25, 1)的源程序：function k1= gangdujuzhen( f,E,NU,t,xi,yi,xj,yj,xm,ym,p ) %UNTITLED4 此处显示有关此函数的摘要% 此处显示详细说明k1=gangdujuzhen(@LinearTriangleElementStiffness, 210000000, 0.3, 0.025, 0, 0, 0.5, 0.25, 0, 0.25, 1)图1. 薄平板结构k1= f(E,NU,t,xi,yi,xj,yj,xm,ym,p);>>k1=gangdujuzhen(@LinearTriangleElementStiffness, 210000000, 0.3, 0.025, 0, 0, 0.5, 0.25, 0, 0.25, 1)k1 =1.0e+06 *2.0192 0 0 -1.0096 -2.0192 1.00960 5.7692 -0.8654 0 0.8654 -5.76920 -0.8654 1.4423 0 -1.4423 0.8654-1.0096 0 0 0.5048 1.0096 -0.5048-2.0192 0.8654 -1.4423 1.0096 3.4615 -1.87501.0096 -5.7692 0.8654 -0.5048 -1.8750 6.2740k2=gangdujuzhen2(@LinearTriangleElementStiffness, 210000000, 0.3, 0.025, 0, 0, 0.5, 0, 0.5, 0.25, 1)的源程序：function k2= gangdujuzhen2( f,E,NU,t,xi,yi,xj,yj,xm,ym,p )%UNTITLED3 此处显示有关此函数的摘要% 此处显示详细说明k2=gangdujuzhen2(@LinearTriangleElementStiffness, 210000000, 0.3, 0.025, 0, 0, 0.5, 0, 0.5, 0.25, 1)k2= f(E,NU,t,xi,yi,xj,yj,xm,ym,p);>>k2=gangdujuzhen2(@LinearTriangleElementStiffness, 210000000, 0.3, 0.025, 0, 0, 0.5, 0, 0.5, 0.25, 1)k2 =1.0e+06 *1.4423 0 -1.4423 0.8654 0 -0.86540 0.5048 1.0096 -0.5048 -1.0096 0-1.4423 1.0096 3.4615 -1.8750 -2.0192 0.86540.8654 -0.5048 -1.8750 6.2740 1.0096 -5.76920 -1.0096 -2.0192 1.0096 2.0192 0-0.8654 0 0.8654 -5.7692 0 5.7692（3）集成整体刚度矩阵⨯的，因此为了得到整体刚度矩阵K，我们由于结构有4个节点，所以刚度矩阵是88⨯的零矩阵。
(完整版)有限元大作业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的平面应力问题研究
基于matlab的平面应力问题研究近年来，随着科学技术的发展，研究者不断对平面应力问题进行研究，并采用不同的方法来模拟和解决这一问题。
其中，借助 MATLAB 术可以高效地完成这项研究。
本文首先简要介绍了MATLAB的基本知识，然后根据现有研究成果，总结出了应用MATLAB解决平面应力问题的典型方法，并结合实际应用，进一步分析其优缺点。
此外，本文还提出了在计算机程序中避免平面应力问题的策略以及可行性研究。
MATLAB是计算机软件，适用于数值计算、科学计算、可视化和交互式编程等多种方面。
它的实用性应用和强大的数据处理能力，使得MATLAB成为当前应用最广泛的计算软件。
目前，MATLAB于平面应力问题的研究，主要有以下两种方法：第一种是使用 MATLAB非线性方程求解方法，即数值求解法。
要解决一个复杂的平面应力问题，需要先建立一个非线性方程组，然后使用 MATLAB行求解。
这种方法的优点是：计算快速，求解精度高，并且可以解决大规模、非线性问题。
但是，它的缺点是：由于非线性方程的求解非常复杂，需要花费较长的时间才能获得准确的结果。
第二种研究方法是使用MATLAB的有限元分析方法。
有限元分析是一种采用数学建模技术来对物体受力过程详细分析的方法。
它可以用来分析平面应力问题，从而更好地了解平面应力的变化趋势，为后续的控制策略提供理论依据。
优点是：可以在计算机上快速求解非线性方程，以及通过二维或三维数值模拟，更准确地分析问题；缺点是：求解精度较低，计算时间较长，技术要求较高。
在计算机程序开发过程中，为了避免出现平面应力问题，应坚持以下步骤：首先，要持续关注平面应力的变化趋势，并分析各种可能影响应力变化的因素，以便找到平面应力的增加的原因；其次，要确定加强应力控制的策略，以避免发生平面应力大小的变化；最后，结合实际情况，应控制计算机程序的开发进度，确保软件系统的性能达到预期目标。
本文分析了MATLAB技术解决平面应力问题的典型方法，以及计算机程序中避免平面应力问题的策略，对现有的研究工作提供了一定的参考。
基于matlab的平面应力问题研究
基于matlab的平面应力问题研究
把握平面应力问题的发展趋势，不仅可以更好地解决工程问题，而且可以更精准地预测结构安全性。
Matlab作为一种计算数学软件，可以模拟平面应力问题，并有效解决相关问题，因此基于Matlab技术的平面应力分析研究受到社会各界的高度重视。
一、Matlab技术在平面应力分析中的应用
（1）Matlab技术可以求解平面应力分析问题，实现非线性状态的力学分析，可以模拟结构的动态变化，以最大程度地测量结构的稳定性。
（2）Matlab技术可以模拟复杂的平面应力分析问题，预测应力的变化情况，确定各种不同环境平面应力的影响因素，以确保结构的安全性。
（3）Matlab技术可以在非线性应力状态下进行模拟，实现对材料强度及变形性能的准确评估，以分析材料对结构安全性的影响，从而更好地保证安全可靠。
二、基于Matlab的平面应力分析的研究方法
（1）首先，建立数学模型，包括应力场模型、变形模型等，以用于模拟应力场分布，研究基于Matlab的平面应力分析问题；
（2）其次，运用Matlab对建立的数学模型进行对应的数值求解，将平面应力分析问题转化为适当的线性或非线性方程；
（3）最后，运用Matlab进行数据分析，并提出有效的解决方案，实现平面应力分析问题的有效解决。
三、结论
Matlab技术在平面应力分析中的应用使其能够实现非线性状态的力学分析，借助Matlab对复杂的平面应力分析问题进行有效模拟，研究分析实现了结构安全性的预测。
Matlab技术为平面应力分析提供了有效的计算方法，使得分析更准确，有助于提高工程的安全性及可靠性。
在未来，基于Matlab 的平面应力分析研究将不断发展，以更好地解决工程应用中的问题。
基于matlab的平面应力问题研究
基于matlab的平面应力问题研究近年来，随着研究的深入及科技的迅速发展，应力学问题受到了越来越多的关注。
应力学问题是工程和物理科学研究的基本内容，它主要涉及金属、机械、材料、结构力学等多个领域。
应力学问题可以从结构力学分析、材料性能预测及有关地震、水力、气动力学等几方面考虑，其中研究应力分布最为重要。
由于应力学问题涉及到多个领域，要研究其解析解和数值模拟解是一项极具挑战性的工作，针对应力分布的研究，主要是利用应力分析的计算机软件，如ANSYS、Comsol等进行模拟研究。
本研究采用Matlab作为计算工具，利用其强大的编程功能以及各种数值算法，构建应力分布模拟系统，对平面应力问题进行研究分析。
首先通过详细的理论推导，导出二维平面应力的基本的计算公式。
然后，利用Matlab的编程，构建一个应力分布模拟系统，根据一维变形规律，分析平面应力的分布情况，采用数值的方法研究不同数据的变化规律。
在此基础上，结合经典的应力模型，探讨不同材料的应力分布情况，并对相关的细节参数进行深入研究，为进行模型模拟提供更加准确的参数。
既然研究计算机模拟，则需要考虑实际情况下的问题，本研究中采用Matlab作为编程语言，首先利用Matlab的视觉技术、绘图工具和图表化工具，构建出一个多维的动态模拟系统，实时展现出平面应力随着外界影响变化的情况，以便分析应力分布中出现的问题。
其次，利用Matlab的数值计算工具，对计算出来的数据进行处理，最后可以求出应力随着外界因素变化的参数分析结果，从而为实际工程提供设计依据。
本文结合Matlab软件与其他相关软件，通过数值模拟的方式，研究了平面应力问题的分布规律，从而确定应力的大小及受力是如何变化的，从而为实际工程设计提供参考。
本研究所得出的结果可以提供一种新的应力分析方法，为现有应力分析手段的改进提供一种有效性的检验方法，能够更加准确、快捷的对平面应力问题进行分析。
综上所述，本文基于Matlab软件，利用其强大的编程功能和数值算法，构建应力分布模拟系统，研究了平面应力问题的分布情况，并对不同材料的应力分布情况作了详细的分析，为工程设计提供了参考。
根据MATLAB的有限元法分析平面应力应变问答刘刚
姓名：刘刚学号：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=210e6 E =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.8750 Column 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.8654 Columns 6 through 8-1.0096 -2.0192 1.00960 0.8654 -5.76920 0 0* *0 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 01.8750- 0.8654 1.4423- 0 3.4615 1.00962.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- 03.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 01.8750- 0.8654 1.4423- 0 3.4615 1.00962.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- 03.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 0 0.8654 -5.7692 0 6.2740>> f=[9.375;0;9.375;0] f =9.3750 0 9.3750 0>> u=k\fu =1.0e-005 *0.71110.11150.65310.0045现在可以清楚的看出，节点2的水平位移和垂直位移分别是0.7111m和0.1115m。
基于有限元法的物体受力变形和应力研究
基于有限元法的物体受力变形和应力研究作者：王重来源：《科技视界》 2014年第17期王重（上海飞机设计研究院环控氧气系统设计研究部，中国上海 201210）【摘要】本文以工程中常见的弯曲悬臂梁作为研究对象，基于有限元分析方法对悬臂梁进行计算模型建立，并求解出悬臂梁在受力后的变形轮廓和应力分布，为其设计和优化提供了验证和指导。
结果表明，该悬臂梁在受力时中部位置应力较大，为该工况下的薄弱区域，需要进行局部加强，以防止在使用中发生失效。
【关键词】弯曲悬臂梁；有限元；刚度矩阵；变形；应力；应变0 引言物体受力变形和应力分析是工程中的常见问题，它在受力零件设计过程中是不可或缺的重要工作，例如在飞机承力梁、高压导管支架等的设计阶段对其将来在飞行中承受载荷后会出现的变形量、变形后轮廓进行估计以及完成强度校核。
这类问题的实质是经典弹性问题，它们的数学模型一般都是一组具有相应边界条件和初值的微分方程，这些微分方程组的解析解能够向我们展示出精确且完整的系统行为，也就是我们所需要的分析结果[1-2]。
但由于这些微分方程组的复杂性，我们又往往无法得到它们的解析解，这时我们就需要利用数值方法来求出近似解，这时在系统中各“节点”的数值解近似于解析解。
因此我们在使用数值方法进行求解前需要对“节点”和“单元”进行合理划分和定义，也就是“离散化”[2-3]。
在此过程之后我们再使用数值解法对问题进行求解。
有限元法是工程中常用的一种数值解法，它使用积分方法来建立系统的代数方程组，用一个连续的函数来近似描述每个单元的解，正因为每个单元的边界是连续的，因此整个系统的解可以由每一个单个的解“组装”起来[2-3]。
不难理解，当我们所划分的单元趋近于无穷多时，使用有限元法所得的解会趋近于精确解。
本文将基于有限元法的基本原理，对具体受力物体进行变形和应力的分析和研究，求解出物体的形变轮廓和相应的应力分布。
1 物体受力工况图1所示为工程中常见的弯曲悬臂梁。
基于matlab的平面应力问题研究
基于matlab的平面应力问题研究近年来，随着物理模拟技术和计算机技术的发展，应用Matlab 从事物理模拟及仿真的研究已经成为重要的学科。
Matlab综合多种数学运算功能及图像处理技术，能够快速、有效地解决各种复杂物理模拟及仿真问题，因此，在电力、能源及制造等多种领域得到了广泛的应用。
本文以Matlab为工具，从物理学的角度，结合一些物理模型，研究平面应力问题。
平面应力问题是几何力学理论中最基本的问题之一，它是研究物体内部力学系统状态的重要数学方法。
根据力学原理，平面应力问题可分为三种，分别是平面压强问题、拉力问题和推力问题，它们的计算主要着重于求解物体内部的力学系统的各种分量的模型。
首先，基于Matlab应用程序，利用拉格朗日乘子法和有限元法，构建平面应力问题的数学模型，并对模型进行参数估计，以求出最优解。
然后，利用Matlab中内置的物理仿真引擎对模型进行数值仿真。
最后，通过Matlab软件对模型进行可视化和调试，实现更为直观的应力分布及力学分析。
Matlab可以解决复杂的物理模拟问题，因此，基于Matlab的平面应力分析在工程实践中广泛使用。
近年来，许多工程师和数学家均使用Matlab研究平面应力问题，在结构力学、机械设计及结构抗震等方面取得重大进展，丰富和拓展了Matlab应用领域，提高了结构性能，为结构抗震和可靠性分析提供了有力的技术支持。
除此之外，本文还介绍了Matlab的其它应用，如图像处理、声音识别等，可用于平面应力问题的分析和计算。
基于Matlab的图像处理技术可以模拟和模拟任何形状的结构以及物体的动态运动，为平面应力问题提供了一种新的解决方法，大大提高了计算的精确性和可靠性。
另外，基于Matlab的声音处理技术也可用于研究平面应力问题，通过捕捉声音信号，可以使这些信号可视化，从而更好地理解平面应力问题产生的原因。
本文介绍了基于Matlab的平面应力分析的基本方法及其在工程实践中的重要性，同时介绍了Matlab的其它应用，如图像处理和声音处理等。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
姓名：刘刚学号：15平面应力应变分析有限元法Abstruct:本文通过对平面应力/应变问题的简要理论阐述，使读者对要分析的问题有大致的印象，然后结合两个实例，通过MATLAB软件的计算，将有限元分析平面应力/应变问题的过程形象的展示给读者，让人一目了然，快速了解有限元解决这类问题的方法和步骤！一.基本理论有限元法的基本思路和基本原则以结构力学中的位移法为基础，把复杂的结构或连续体看成有限个单元的组合，各单元彼此在节点出连接而组成整体。
把连续体分成有限个单元和节点，称为离散化。
先对单元进行特性分析，然后根据节点处的平衡和协调条件建立方程，综合后做整体分析。
这样一分一合，先离散再综合的过程，就是把复杂结构或连续体的计算问题转化简单单元分析与综合问题。
因此，一般的有限揭发包括三个主要步骤：离散化单元分析整体分析。
二.用到的函数1. LinearTriangleElementStiffness(E,NU,t,xi,yi,xj,yj,xm,ym,p)(K k I f)(k u)(k u A)(E NU t)三.实例例1.考虑如图所示的受均布载荷作用的薄平板结构。
将平板离散化成两个线性三角元，假定E=200GPa，v=，t=0.025m,w=3000kN/m.1.离散化2.写出单元刚度矩阵通过matlab 的LinearTriangleElementStiffness 函数，得到两个单元刚度矩阵1k 和2k ，每个矩阵都是6 6的。
>> E=210e6 E =>> k1=LinearTriangleElementStiffness(E,NU,t,0,0,,,0,,1) k1 = +006 *Columns 1 through 50 0 0 0 0 0 0 0Column 6 >> NU= NU = >> t= t =>> k2=LinearTriangleElementStiffness(E,NU,t,0,0,,0,,,1)k2 =+006 *Columns 1 through 50 00 0Column 63.集成整体刚度矩阵 8*8零矩阵K =0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 >> K=LinearTriangleAssemble(K,k1,1,3,4)K =+006 *Columns 1 through 50 0 0 00 0 00 0 0 0 00 0 0 0 0 0 0 00 0 0 0 0 0 0 0Columns 6 through 80 0 0 0 0 0 0>> K=LinearTriangleAssemble(K,k1,1,2,3) K =+007 *0 0 0 0 0 00 0 0 0 0 0 0 0 0 00 0 0 04.引入边界条件.用上一步得到的整体刚度矩阵，可以得到该结构的方程组如下形式⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡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 01.8750- 0.8654 1.4423- 0 3.4615 1.00962.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- 03.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 01.8750- 0.8654 1.4423- 0 3.4615 1.00962.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- 03.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 =+006 *0 0>> f=[;0;;0] f = 0 0 >> u=k\f u =*现在可以清楚的看出，节点2的水平位移和垂直位移分别是0.7111m 和0.1115m 。
节点3的水平位移和垂直位移分别是0.6531m 和0.0045m 。
6.后处理用matlab 命令求出节点1和节点4的支反力以及每个单元的应力。
首先建立总体节点位移矢量U ，U=[0;0;u;0;0] U = *0 0 0 0 >> F=K*U F =由以上知，节点1的水平反力和垂直反力分别是（指向左边）和（作用力方向向下），节点4的水平反力和垂直反力分别是（指向左边）和（作用力方向向下）.满足力平衡条件。
接着，建立单元节点位移矢量21u 和u ，然后调用matlab 命令LinearTriangleElementStresses 计算单元应力sigma1和sigma2>> u1=[U(1);U(2);U(5);U(6);U(7);U(8)]u1 =*>> u2=[U(1);U(2);U(3);U(4);U(5);U(6)]u2 =*>> sigma1=LinearTriangleElementStresses(E,NU,,0,0,,,0,,1,u1) sigma1 =+003 *>> sigma2=LinearTriangleElementStresses(E,NU,,0,0,,0,,,1,u2) sigma2 =+003 *由以上可知，单元1的应力)(0144.3拉应力MPa x =σ，)(9043.0y 拉应力MPa =σ，)(0072.0y 正值MPa x =τ 。
单元2的应力是)(9856.2拉应力MPa x =σ)(0036.0y 压应力MPa =σ)(0072.0y 负值MPa x =τ。
显然，在x方向的应力（拉应力）接近于正确的值3MPa （拉应力）。
接着调用LinearTriangleElementStresses 函数计算每个单元的主应力和主应力方向角。
>> s1= LinearTriangleElementPStresses(sigma1) s1 = +003 *>> s2= LinearTriangleElementPStresses(sigma2) s2 = +003 *)(0144.31拉应力MPa =σ，)(9043.02拉应力MPa =σ主应力方向角ο2.0=p θ )(M 9856.21拉应力Pa =σ，)(0036.02压应力MPa =σ，ο1.0-=p θ例 2.考虑如图所示的由均匀分布载荷和集中载荷作用的薄平板结构。
将平板离散化成12个线性三角单元，如图4所示。
假定E=210GPa，v=，t=0.025m，w=100kN/m和P=。
1.离散化2.写出单元刚度矩阵>> E=201e6;>> NU=;>> t=;>> k1= LinearTriangleElementStiffness(E,NU,t,0,,,,,,1); >> k2= LinearTriangleElementStiffness(E,NU,t,0,,0,,,,1); >> k3= LinearTriangleElementStiffness(E,NU,t,,,,,,,1); >> k4= LinearTriangleElementStiffness(E,NU,t,,,0,,,,1); >> k5= LinearTriangleElementStiffness(E,NU,t,0,,,,,,1); >> k6= LinearTriangleElementStiffness(E,NU,t,0,,0,0,,,1); >> k7= LinearTriangleElementStiffness(E,NU,t,,,,,,0,1); >> k8= LinearTriangleElementStiffness(E,NU,t,,,0,0,,0,1); >> k9= LinearTriangleElementStiffness(E,NU,t,025,,,0,,,1); >> k10= LinearTriangleElementStiffness(E,NU,t,,,,,,,1); >> k11= LinearTriangleElementStiffness(E,NU,t,,0,,0,,,1); >> k12= LinearTriangleElementStiffness(E,NU,t,,,,0,,,1)k1 =+006 *3.集成整体刚度矩阵：>>K=zero(22,22);>>K=LinearTriangleAssemble(K,k1,1,3,2);>>K=LinearTriangleAssemble(K,k2,1,4,3);>>K=LinearTriangleAssemble(K,k3,3,5,2);>>K=LinearTriangleAssemble(K,k4,3,4,5);>>K=LinearTriangleAssemble(K,k5,4,6,5);>>K=LinearTriangleAssemble(K,k6,4,7,6);>>K=LinearTriangleAssemble(K,k7,5,6,8);>>K=LinearTriangleAssemble(K,k8,6,7,8);>>K=LinearTriangleAssemble(K,k9,5,8,9);>>K=LinearTriangleAssemble(K,k10,5,9,10);>>K=LinearTriangleAssemble(K,k11,8,11,9);>>K=LinearTriangleAssemble(K,k12,9,11,10)运行得+008 *Columns 1 through 70 00 0 00 00 0 00 0 0 0 0 00 0 0 0 0 00 0 0 0 0 00 0 0 0 0 00 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0Columns 8 through 140 0 0 0 0 00 0 0 0 0 00 0 0 0 00 0 0 0 00 0 0 00 0 0 00 00 0 00 00 00 00 0 0 0 00 0 0 0 00 0 0 0 00 0 0 0 00 0 0 0 0 0 0 0 0 0 0 0 0 0Columns 15 through 210 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 00 0 0 0 00 0 0 0 00 0 0 0 00 0 0 0 00 00 00 00 0Column 224.引入边界条件：U1x= U1y= U4x= U4y=U7x=U7y=0F2x= F2y= F3x= F3y=F6x=F6y=F8x= F8y= F9x= F9y=F10x=F10y= F11x= F11y= 0F5x= 0,F5y=5.解方程：>>k=[K(3:6,3:6),K(3:6,9:12),K(3:6,15:22);K(9:12,3:6),K(9:12,9:12),K(9:12,15:22) ;K(15:22,3:6),K(15:22,9:12) ,K(15:22,15:22)];K =+008 *Columns 1 through 80 00 00 0 00 0 00 0 00 0 00 0 0 0 0 00 0 0 0 0 00 0 0 0 0 00 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 00 0 0 0 0 0 0 0Columns 9 through 160 0 0 0 0 0 0 00 0 0 0 0 0 0 00 0 0 0 0 00 0 0 0 0 00 0 0 0 0 00 0 0 0 0 00 0 00 0 00 00 00 00 00 0 0 00 0 0 00 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0Columns 17 through 220 0 0 0 0 00 0 0 0 0 00 0 0 0 0 00 0 0 0 0 00 0 0 0 0 00 0 0 0 0 00 0 0 0 0 00 0 0 0 0 00 00 00 0 0 0 0 00 0 0 0 0 00 0 0 0 0 00 0 0 0 0 00 00 0>>f=[0;0;0;0;0;;0;0;0;0;0;0;0;0;0;0];>>u=k\fu=*[ ]T6.后处理>>U=[0;0;u(1:4);0;0;u(5:8);0;0;u(9:16)];>>F=K*UF=[ 0 0 ]T三.小结通过这次练习，对有限元结合MATLAB软件解决平面应力/应变问题的方法和过程逐渐清晰，特别是在应用MATLAB软件的过程中学到了很多东西。