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

【MATLAB 算例】4.8.1(1) 基于 4 节点四面体单元的空间块体分析 (Tetrahedron3D4Node)如图 4-22 所示的一个块体，在右端面上端点受集中力 F 作用。基于 MATLAB 平台，计算各个节点位移、支反力以及单元的应力。取相关参数为： ，10Pa,=.25E。5=10NF图 4-22 一个空间块体的分析解答：对该问题进行有限元分析的过程如下。(1)结构的离散化与编号将结构离散为 5 个 4 节点四面体单元，单元编号及节点编号和坐标如图 4-22 所示，连接关系见表 4-8，节点的坐标见表 4-9。表 4-8 单元连接关系单元号 节点号123451 4 2 61 4 3 76 7 5 16 7 8 41 4 6 7表 4-9 节点的坐标节点 节点坐标/mx y z123456780 0 00.2 0 00 0.8 00.2 0.8 00 0 0.60.2 0 0.60 0.8 0.60.2 0.8 0.6节点位移列阵(4-190) 112288 Tuvwvuvwq节点外载列阵(4-191)34780 0 TTFF其中 34785 010N约束的支反力列阵(4-1921256TTT RR其中 1256 xxxxyyyyzzzzRR总的节点载荷列阵(4-193)12345678 TTTT PFFF(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 个方向的位移将为零，即。因此，将针对节点112 60uvwuvuvwuv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\pu = 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 -33u ＝0.131 4×10 -37u＝-0.048 5×10 -3v ＝0.085 8×10 -3v＝-0.402 4×10 -33w＝-0.446 0×10 -37w＝0.134 3×10 -34 ＝0.135 3×10 -38＝-0.071 5×10 -3v ＝0.068 1×10 -3v＝-0.403 1×10 -34 ＝-0.474 2×10 -38(5)支反力的计算在得到整个结构的节点位移后，由原整体刚度方程就可以计算出对应的支反力；先将上面得到的位移结果与位移边界条件的节点位移进行组合( 注意位置关系) ，可以得到整体的位移列阵 U(24×1)，再代回原整体刚度方程，计算出所有的节点力 P(24×1)，按式(4-192)的对应关系就可以找到对应的支反力。>>U=[zeros(6,1);u([1:6]);zeros(6,1);u(7:12)];>>P=KK*UP = 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 空间块体的支反力计算结果510.37 21xRN 550.47 1xRN4y 3y51.9 z 55.6 z2xx521.89 0yRN 561.289 0yRN4z 7z(6)各单元的应力计算先从整体位移列阵 U(24×1)中提取出单元的位移列阵，然后，调用计算单元应力的函数