基于MATLAB的有限元法分析平面应力应变问题--刘刚

姓名：刘刚学号：15平面应力应变分析有限元法Abstruct:本文通过对平面应力/应变问题的简要理论阐述，使读者对要分析的问题有大致的印象，然后结合两个实例，通过MATLAB软件的计算，将有限元分析平面应力/应变问题的过程形象的展示给读者，让人一目了然，快速了解有限元解决这类问题的方法和步骤！一.基本理论有限元法的基本思路和基本原则以结构力学中的位移法为基础，把复杂的结构或连续体看成有限个单元的组合，各单元彼此在节点出连接而组成整体。把连续体分成有限个单元和节点，称为离散化。先对单元进行特性分析，然后根据节点处的平衡和协调条件建立方程，综合后做整体分析。这样一分一合，先离散再综合的过程，就是把复杂结构或连续体的计算问题转化简单单元分析与综合问题。因此，一般的有限揭发包括三个主要步骤：离散化单元分析整体分析。二.用到的函数1.LinearTriangleElementStiffness(E,NU,t,xi,yi,xj,yj,xm,ym,p)2.LinearBarAssemble(KkIf)3.LinearBarElementForces(ku)4.LinearBarElementStresses(kuA)5.LinearTriangleElementArea(ENUt)三.实例例1.考虑如图所示的受均布载荷作用的薄平板结构。将平板离散化成两个线性三角元，假定E=200GPa，v=0.3，t=0.025m,w=3000kN/m.1.离散化2.写出单元刚度矩阵通过matlab的LinearTriangleElementStiffness函数，得到两个单元刚度矩阵1k和2k，每个矩阵都是66的。E=210e6E=210000000k1=LinearTriangleElementStiffness(E,NU,t,0,0,0.5,0.25,0,0.25,1)k1=1.0e+006*Columns1through52.019200-1.0096-2.019205.7692-0.865400.86540-0.86541.44230-1.4423-1.0096000.50481.0096-2.01920.8654-1.44231.00963.46151.0096-5.76920.8654-0.5048-1.8750Column61.0096-5.76920.8654-0.5048-1.87506.2740NU=0.3NU=0.3000t=0.025t=0.0250k2=LinearTriangleElementStiffness(E,NU,t,0,0,0.5,0,0.5,0.25,1)k2=1.0e+006*Columns1through51.44230-1.44230.8654000.50481.0096-0.5048-1.0096-1.44231.00963.4615-1.8750-2.01920.8654-0.5048-1.87506.27401.00960-1.0096-2.01921.00962.0192-0.865400.8654-5.76920Column6-0.865400.8654-5.769205.76923.集成整体刚度矩阵8*8零矩阵K=0000000000000000000000000000000000000000000000000000000000000000K=LinearTriangleAssemble(K,k1,1,3,4)K=1.0e+006*Columns1through52.0192000005.769200-0.865400000000000-0.8654001.4423-1.00960000-2.01920.865400-1.44231.0096-5.7692000.8654Columns6through8-1.0096-2.01921.009600.8654-5.76920000000-1.44230.86540.50481.0096-0.50481.00963.4615-1.8750-0.5048-1.87506.2740K=LinearTriangleAssemble(K,k1,1,2,3)K=1.0e+007*0.403800-0.1010-0.20190-0.20190.101001.1538-0.086500-0.57690.0865-0.57690-0.08650.14420-0.14420.086500-0.1010000.05050.1010-0.050500-0.20190-0.14420.10100.4904-0.1875-0.14420.08650-0.57690.0865-0.0505-0.18750.67790.1010-0.0505-0.20190.086500-0.14420.10100.3462-0.18750.1010-0.5769000.0865-0.0505-0.18750.62744.引入边界条件.用上一步得到的整体刚度矩阵，可以得到该结构的方程组如下形式4y4X3y3X2y2X1y1X4y4X3y3X2y2X1y1X6FFFFFFFFUUUUUUUU6.27401.8750-0.5048-0.8654005.7692-1.00961.8750-3.46151.00961.4423-000.86542.0192-0.5048-1.00966.274005.7692-0.865401.8750-0.86541.4423-03.46151.00962.0192-1.8750-0005.7692-1.00966.27401.8750-0.5048-0.8654000.86542.0192-1.8750-3.46151.00961.4423-5.7692-0.865401.8750-0.5048-1.00966.274001.00962.0192-1.8750-00.86541.4423-03.461510本题的边界条件：04411yxyxUUUU0,375.9,0,375.93322yxyxFFFF将边界条件带入，得到：4y4X1y1X3y3X2y2X6FF09.37509.375FF00UUUU006.27401.8750-0.5048-0.8654005.7692-1.00961.8750-3.46151.00961.4423-000.86542.0192-0.5048-1.00966.274005.7692-0.865401.8750-0.86541.4423-03.46151.00962.0192-1.8750-0005.7692-1.00966.27401.8750-0.5048-0.8654000.86542.0192-1.8750-3.46151.00961.4423-5.7692-0.865401.8750-0.5048-1.00966.274001.00962.0192-1.8750-00.86541.4423-03.4615105.解方程分解上述方程组，提取总体刚度矩阵K的第3-6行的第3-6列作为子矩阵09.37509.375UUUU6.274005.7692-0.865403.46151.00962.0192-5.7692-1.00966.27401.8750-0.86542.0192-1.8750-3.4615103y3X2y2X6Matlab命令k=K(3:6,3:6)k=1.0e+006*3.4615-1.8750-2.01920.8654-1.87506.27401.0096-5.7692-2.01921.00963.461500.8654-5.769206.2740f=[9.375;0;9.375;0]f=9.375009.37500u=k\fu=1.0e-005*0.71110.11150.65310.0045现在可以清楚的看出，节点2的水平位移和垂直位移分别是0.7111m和0.1115m。节点3的水平位移和垂直位移分别是0.6531m和0.0045m。6.后处理用matlab命令求出节点1和节点4的支反力以及每个单元的应力。首先建立总体节点位移矢量U，U=[0;0;u;0;0]U=1.0e-005*000.71110.11150.65310.004500F=K*UF=-9.3750-5.62959.37500.00009.37500.0000-9.37505.6295由以上知，节点1的水平反力和垂直反力分别是9.375kn（指向左边）和5.6295kn（作用力方向向下），节点4的水平反力和垂直反力分别是9.375kn（指向左边）和5.6295kn（作用力方向向下）.满足力平衡条件。接着，建立单元节点位移矢量21u和u，然后调用matlab命令LinearTriangleElementStresses计算单元应力sigma1和sigma2u1=[U(1);U(2);U(5);U(6);U(7);U(8)]u1=1.0e-005*000.65310.004500u2=[U(1);U(2);U(3);U(4);U(5);U(6)]u2=1.0e-005*000.71110.11150.65310.0045sigma1=LinearTriangleElementStresses(E,NU,0.025,0,0,0.5,0.25,0,0.25,1,u1)sigma1=1.0e+003*3.01440.90430.0072sigma2=LinearTriangleElementStresses(E,NU,0.025,0,0,0.5,0,0.5,0.25,1,u2)sigma2=1.0e+003*2.9856-0.0036-0.0072由以上可知，单元1的应力)(0144.3拉应力MPax，)(9043.0y拉应力MPa，)(0072.0y正值MPax。单元2的应力是)(9856.2拉应力MPax)(0036.0y压应力MPa)(0072.0y负值MPax。显然，在x方向的应力（拉应力）接近于正确的值3MPa（拉应力）。接着调用LinearTriangleElementStresses函数计算每个单元的主应力和主应力方向角。s1=LinearTriangleElementPStresses(sigma1)s1=1.0e+003*3.01440.90430.0002s2=LinearTriangleElementPStresses(sigma2)s2=1.0e+003*2.9856-0.0036-0.0001)(0144.31拉应力MPa，)(9043.02拉应力MPa主应力方向角2.0p)(M9856.21拉应力Pa，)(0036.02压应力MPa，1.0p例2.考虑如图3.1所示的由均匀分布载荷和集中载荷作用的薄平板结构。将平板离散化成12个线性三角单元，如图4所示。假定E=210GPa，v=0.3，t=0.025m，w=100kN/m和P=12.5kN。1.离散化2.写出单元刚度矩阵E=201e6;NU=0.3;t=0.025;k1=LinearTriangleElementStiffness(E,NU,t,0,0.5,0.125,0.375,0.25,0.5,1);k2=LinearTriangleElementStiffness(E,NU,t,0,0.5,0,0.25,0.