《有限元基础教程》_【MATLAB算例】基于节点四边形单元的矩形薄板分析(QuadDNode).docVIP

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