《有限元基础教程》-【MATLAB算例】4.7.2(1)基于4节点四边形单元的矩形薄板分析(Quad2D4Node).doc

【MATLAB算例】4.7.2(1) 基于4节点四边形单元的矩形薄板分析(Quad2D4Node) 如图4-21所示的一个薄平板，在右端部受集中力F作用，其中的参数为：。基于MATLAB平台，按平面应力问题计算各个节点位移、支座反力以及单元的应力。 (a) 问题描述 (b) 有限元分析模型图4-21 右端部受集中力作用的薄平板解答：对该问题进行有限元分析的过程如下。（1）结构的离散化与编号 将结构离散为二个4节点矩形单元，单元编号及节点编号如图4-21(b)所示，连接关系见表4-6，节点的几何坐标见表4-7，载荷F按静力等效原则向节点1，2移置。表4-6 结构的单元连接关系单元号节点号123 5 6 41 3 4 2表4-7 节点的坐标节点节点坐标/mx y 1234562 12 0 1 11 00 10 0节点位移列阵 (4-186)节点外载列阵 (4-187)约束的支反力列阵 (4-188)总的节点载荷列阵 (4-189)其中，和分别为节点5和节点6的两个方向的支反力。（2）计算各单元的刚度矩阵(以国际标准单位) 首先在MATLAB环境下，输入弹性模量E、泊松比NU、薄板厚度h和平面应力问题性质指示参数ID，然后针对单元1和单元2，分别两次调用函数Quad2D4Node_Stiffness，就可以得到单元的刚度矩阵k1(88)和k2(88)。 E=1e7; NU=1/3; t=0.1; ID=1; k1= Quad2D4Node_Stiffness(E,NU,t, 1,1, 0,1,0,0,1,0, ID);k2=Quad2D4Node_Stiffness(E,NU,t, 2,1,1,1,1,0,2,0, ID);（3） 建立整体刚度方程由于该结构共有6个节点，则总共的自由度数为12，因此，结构总的刚度矩阵为KK(1212)，先对KK清零，然后两次调用函数Quad2D4Node_Assembly进行刚度矩阵的组装。KK = zeros(12,12);KK= Quad2D4Node_Assembly(KK,k1, 3,5,6,4);KK= Quad2D4Node_Assembly(KK,k2, 1,3,4,2);（4） 边界条件的处理及刚度方程求解 由图4-21(b)可以看出，节点5和节点6的两个方向的位移将为零，即，因此，将针对节点1，2，3和4的位移进行求解。节点1，2，3和4的位移将对应KK矩阵中的前8行和前8列，则需从KK(1212)中提出，置给k，然后生成对应的载荷列阵p，再采用高斯消去法进行求解。注意：MATLAB中的反斜线符号“”就是采用高斯消去法。k=KK(1:8,1:8); p=0;-50000;0;-50000;0;0;0;0;u=kpu = 0.8000 -2.5333 -0.8000 -2.5333 0.6000 -0.8667 -0.6000 -0.8667 将列排成行排量 由此可以看出，所求得的结果为： 。（5）支反力的计算在得到整个结构的节点位移后，由原整体刚度方程就可以计算出对应的支反力。先将上面得到的位移结果与位移边界条件的节点位移进行组合(注意位置关系)，可以得到整体的位移列阵U(121)，再代回原整体刚度方程，计算出所有的节点力P(121)，按式(4-189)的对应关系就可以找到对应的支反力。U=u;0;0;0;0;P=KK*UP = 1.0e+005 * -0.0000 -0.5000 0.0000 -0.5000 0.0000 0.0000 将列排成行排量 0.0000 -0.0000 -2.0000 0.5000 2.0000 0.5000 将列排成行排量 由式(4-189)的对应关系，可以得到对应的支反力为。（6）各单元的应力计算先从整体位移列阵U(121)中提取出单元的位移列阵，然后，调用计算单元应力的函数Quad2D4Node_Stress，就可以得到各个单元的应力分量。u1= U(5);U(6) ;U(9);U(10);U(11);U(12);U(7);U(8)u1 = 0.6000 -0.8667 0 0 0 0 -0.6000 -0.8667 将列排成行排量 stress1=Quad2D4Node_Stress(E,NU, 1,1,0,1,0,0,1,0, u1,ID)stress1 = 1.0e+005 *0.0000 0.0000 -10.0000 将列排成行排量 u2= U(1);U(2) ;U(5);U(6); U(7);U(8);U(3);U(4)u2 = 0.8000 -2.5333 0.6000 -0.8667 将列排成行排量 -0.6000 -0.8667 -0.8000 -2.5333 将列排成行排量 stress2=Quad2D4Node_Stress(E,NU, 2,1,1,1,1,0,2,0, u2,ID)stress2 = 1.0e+005 * 0.0000 0.0000 -10.0000 将列排成行排量 可以看出：计算得到的单元1的应力分量为 ；单元2的应力分量为 。