《有限元基础教程》-【MATLAB算例】4.8.1(1)基于4节点四面体单元的空间块体分析(Tetrahedron3D4Node).doc

MATLAB算例】4.8.1(1) 基于4节点四面体单元的空间块体分析(Tetrahedron3D4Node)
如图4-22所示的一个块体，在右端面上端点受集中力F作用。基于MATLAB平台，计算各个节点位移、支反力以及单元的应力。取相关参数为：，。

图4-22 一个空间块体的分析

解答：对该问题进行有限元分析的过程如下。
（1）结构的离散化与编号
将结构离散为5个4节点四面体单元，单元编号及节点编号和坐标如图4-22所示，连接关系见表4-8，节点的坐标见表4-9。
表4-8 单元连接关系
单元号
节点号
1
2
3
4
5
1 4 2 6
1 4 3 7
6 7 5 1
6 7 8 4
1 4 6 7

表4-9 节点的坐标
节点
节点坐标/m
x
y
z

1
2
3
4
5
6
7
8

0
0
0
0.2
0
0
0
0.8
0
0.2
0.8
0
0
0
0.6
0.2
0
0.6
0
0.8
0.6
0.2
0.8
0.6

节点位移列阵
(4-190)
节点外载列阵
(4-191)
其中

约束的支反力列阵
（4-192
其中

总的节点载荷列阵
(4-193)

（2）计算各单元的刚度矩阵(以国际标准单位)
首先在MATLAB环境下，输入弹性模量E、泊松比NU，然后针对单元1和单元2，分别5次调用函数Tetrahedron3D4Node_Stiffness，就可以得到单元的刚度矩阵k1(6×6) ～ k5(6×6)。
>> E=1e10;
>> NU=0.25;
>> k1 = Tetrahedron3D4Node_Stiffness(E,NU,0,0,0,0.2,0.8,0,0.2,0,0,0.2,0,0.6);
>> k2 = Tetrahedron3D4Node_Stiffness(E,NU,0,0,0,0.2,0.8,0,0,0.8,0,0,0.8,0.6);
>> k3 = Tetrahedron3D4Node_Stiffness(E,NU,0.2,0,0.6,0,0.8,0.6,0,0,0.6,0,0,0);
>> k4=Tetrahedron3D4Node_Stiffness(E,NU,0.2,0,0.6,0,0.8,0.6,0.2,0.8,0.6,0.2,0.8,0);
>> k5 = Tetrahedron3D4Node_Stiffness(E,NU,0,0,0,0.2,0.8,0,0.2,0,0.6,0,0.8,06);
（3） 建立整体刚度方程
由于该结构共有8个节点，则总共的自由度数为24，因此，结构总的刚度矩阵为KK(24×24)，先对KK清零，然后5次调用函数Tetrahedron3D4Node_Assembly进行刚度矩阵的组装。
>>KK = zeros(24);
>> KK = Tetrahedron3D4Node_Assembly(KK,k1,1,4,2,6);
>> KK = Tetrahedron3D4Node_Assembly(KK,k2,1,4,3,7);
>> KK = Tetrahedron3D4Node_Assembly(KK,k3,6,7,5,1);
>> KK = Tetrahedron3D4Node_Assembly(KK,k4,6,7,8,4);
>> KK = Tetrahedron3D4Node_Assembly(KK,k5,1,4,6,7);
（4） 边界条件的处理及刚度方程求解
由图4-22可以看出，节点1，2，5和6上3个方向的位移将为零，即。因此，将针对节点3，4，7和8的位移进行求解，节点1，2，5和6的位移将对应KK矩阵中的第1~6行，第13~18行和第1~6列，第13~18列，需从KK(24×24)中提出，置给k，然后生成对应的载荷列阵p，再采用高斯消去法进行求解。注意：MATLAB中的反斜线符号“\”就是采用高斯消去法。
>>k=KK([7:12,19:24],[7:12,19:24]);
>>p=[0,0,0,0,0,0,0,0,-1e5,0,0,-1e5]'
>>u=k\p
u = 1.0e-003 *
0.1249 -0.0485 -0.4024 0.1343 -0.0715 -0.4031 [将列排成行排量 []]
0.1314 0.0858 -0.4460 0.1353 0.0681 -0.4742 [将列排成行排量 []]
所求得的位移结果见表4-10。

表4-10 空间块体的节点位移计算结果
＝0.124 9×10-3
＝0.131 4×10-3
＝-0.048 5×10-3

＝0.085 8×10-3
＝-0.402 4×10-3
＝-0.446 0×10-3
＝0.134 3×10-3
＝0.135 3×10-3
＝-0.071 5×10-3
＝0.068 1×10-3
＝-0.403 1×10-3
＝-0.474 2×10-3
（5）支反力的计算
在得到整个结构的节点位移后，由原整体刚度方程就可以计算出对应的支反力；先将上面得到的位移结果与位移边界条件的节点位移进行组合(注意位置关系)，可以得到整体的位移列阵U(24×1)，再代回原整体刚度方程，计算出所有的节点力P(24×1)，按式(4-192)的对应关系就可以找到对应的支反力。
>>U=[zeros(6,1);u([1:6]);zeros(6,1);u(7:12)];
>>P=KK*U
P = 1.0e+005 *
0.3372 1.3774 0.1904 -0.4202 1.2892 0.4984 [将列排成行排量 []]
-0.0000 0.0000 0.0000 -0.0000 -0.0000 -0.0000 [将列排成行排量 []]
-0.4745 -1.3774 0.5604 0.5575 -1.2892 0.7509 [将列排成行排量 []]
-0.0000 -0.0000 -1.0000 -0.0000 0.0000 -1.0000 [将列排成行排量 []]
由式(4-193)的对应关系，可以得到对应的支反力见表4-11。

表4-11 空间块体的支反力计算结果

（6）各单元的应力计算
先从整体位移列阵U(24×1)中提取出单元的位移列阵，然后，调用计算单元应力的函数Tetrahedron3D4Node_Stress，就可以得到各个单元的应力分量。
>>u1=[U(1:3);U(10:12);U(4:6);U(16:18)];
>>stress1 = Tetrahedron3D4Node_Stress(E,NU,0,0,0,0.2,0.8,0,0.2,0,0,0.2,0,0.6,u1)
stress1 = 10e+006 *
-0.3574 -1.0721 -0.3574 0.6717 -2.0155 0 [将列排成行排量 []]
>>u2=[U(1:3);U(10:12);U(7:9);U(19:21)];
>> stress2= Tetrahedron3D4Node_Stress(E,NU,0,0,0,0.2,0.8,0,0,0.8,0,0,0.8,0.6,u2)
stress2 = 1.0e+006 *
0.0314 -0.8298 -0.9260 0.1649 -1.1170 0.0294 [将列排成行排量 []]
>> u3=[U(16:21);U(13:15);U(1:3)];
>> stress3=Tetrahedron3D4Node_Stress(E,NU,0.2,0,0.6,0,0.8,0.6,0,0,0.6,0,0,0,u3)
stress3 = 1.0e+006 *
0.4289 1.2867 0.4289 0.6568 -2.2301 0 [将列排成行排量 []]

>> u4=[U(16:21);U(22:24);U(10:12)];
>> stress4=Tetrahedron3D4Node_Stress(E,NU,0.2,0,0.6,0,0.8,0.6,0.2,0.8,0.6,0.2,0.8,0,u4)
stress4 = 1.0e+006 *
0.1046 0.6272 -1.0012 0.3233 -1.4402 -0.5562 [将列排成行排量 []]
>> u5=[U(1:3);U(10:12);U(16:21)];
>> stress5=Tetrahedron3D4Node_Stress(E,NU,0,0,0,0.2,0.8,0,0.2,0,0.6,0,0.8,0.6,u5)
stress5 = 1.0e+006 *
-0.0179 -0.0060 -0.3636 -0.9083 -1.5986 0.4192 [将列排成行排量 []]

各个单元应力分量的计算结果列在表4-12中。

表4-12 空间块体的各个单元应力分量的计算结果
1号单元

2号单元

3号单元

4号单元

5号单元