《有限元基础教程》_【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,0.6); （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