《有限元基础教程》-【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 结构的单元连接关系
单元号
节点号
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.5000 2.0000 0.5000 [将列排成行排量 []]
由式(4-189)的对应关系，可以得到对应的支反力为。

（6）各单元的应力计算
先从整体位移列阵U(12×1)中提取出单元的位移列阵，然后，调用计算单元应力的函数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的应力分量为 。