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

【MATLAB算例】4.8.1(1)基于4节点四面体单元的空间块体分析(Tetrahedron3D4Node)如图4-22所示的一个块体，在右端面上端点受集中力F作用。基于MATLAB平台，计算各个节点位移、支反力以及单元的应力。取相关参数为：10110Pa,=0.25E，5=110NF。图4-22一个空间块体的分析解答：对该问题进行有限元分析的过程如下。（1）结构的离散化与编号将结构离散为5个4节点四面体单元，单元编号及节点编号和坐标如图4-22所示，连接关系见表4-8，节点的坐标见表4-9。表4-8单元连接关系单元号节点号1234514261437675167841467表4-9节点的坐标节点节点坐标/mxyz123456780000.20000.800.20.80000.60.200.600.80.60.20.80.6节点位移列阵111222888Tuvwuvwuvwq(4-190)节点外载列阵34780000TTTTTFFFFF(4-191)其中3478500000110NFFFF约束的支反力列阵12560000TTTTTRRRRR（4-192其中1256112255661256xxxxyyyyzzzzRRRRRRRRRRRRRRRR总的节点载荷列阵12345678TTTTTTTTTPFRRRRFFRRFF(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个方向的位移将为零，即1112225556660uvwuvwuvwuvw。因此，将针对节点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.40240.1343-0.0715-0.4031[将列排成行]0.13140.0858-0.44600.13530.0681-0.4742[将列排成行]所求得的位移结果见表4-10。表4-10空间块体的节点位移计算结果3u＝0.1249×10-37u＝0.1314×10-33v＝-0.0485×10-37v＝0.0858×10-33w＝-0.4024×10-37w＝-0.4460×10-34u＝0.1343×10-38u＝0.1353×10-34v＝-0.0715×10-38v＝0.0681×10-34w＝-0.4031×10-38w＝-0.4742×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*UP=1.0e+005*0.33721.37740.1904-0.42021.28920.4984[将列排成行]-0.00000.00000.0000-0.0000-0.0000-0.0000[将列排成行]-0.4745-1.37740.56040.5575-1.28920.7509[将列排成行]-0.0000-0.0000-1.0000-0.00000.0000-1.0000[将列排成行]由式(4-193)的对应关系，可以得到对应的支反力见表4-11。表4-11空间块体的支反力计算结果510.337210xRN550.474510xRN511.377410yRN551.377410yRN510.190410zRN550.560410zRN520.420210xRN560.557510xRN521.289210yRN561.289210yRN520.498410zRN560.750910zRN（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=1.0e+006*-0.3574-1.0721-0.35740.6717-2.01550[将列排成行]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.92600.1649-1.11700.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.42891.28670.42890.6568-2.23010[将列排成行]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.10460.6272-1.00120.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.59860.4192[将列排成行]各个单元应力分量的计算结果列在表4-12中。表4-12空间块体的各个单元应力分量的计算结果1号单元0.3574xMPa1.0721yMPa0.3574zMPa0.6717xyMPa2.0155yzMPa0zxMPa2号单元0.0314xMPa0.8298yMPa0.9260zMPa0.1649xyMPa1.1170yzMPa0.0294zxMPa3号单元0.4289xMPa1.2867yMPa0.4289zMPa0.6568xyMPa2.2301yzMPa0zxMPa4号单元0.1046xMPa0.6272yMPa1.0012zMPa0.3233xyMPa1.4402yzMPa0.5562zxMPa5号单元0.0179xMPa0.0060yMPa0.3636zMPa0.9083xyMPa1.5986yzMPa0.4192zxMPa