《有限元基础教程》_【MATLAB算例】4.8.2(1) 基于8节点六面体单元的空间块体分析(Hexahedral3D8Node)

如图4-23所示的一个空间块体，在右端部受两个集中力F作用，其中的参数为：

E11010Pa,=0.25,t=0.2m,F=1105N。基于MATLAB平台，用一个空间8节点六

面体单元计算各个节点位移、支座反力以及单元的应力。

(a) 问题描述 (b) 有限元分析模型

图4-23 右端部受集中力作用的空间块体

解答：对该问题进行有限元分析的过程如下。 （1）结构的离散化与编号

将结构离散为一个8节点六面体单元，节点编号如图4-23(b)所示，节点的几何坐标见表4-13。

表4-13 节点的坐标

节点 1 2 3 4 5 6 7 8

节点位移列阵

qe[u1(241)节点坐标/m x 0.2 0.2 0 0 0.2 0.2 0 0

y 0 0.8 0.8 0 0 0.8 0.8 0

z 0 0 0 0 0.6 0.6 0.6 0.6

v1w1u2v2w2u8v8w8]T （4-194）

总的节点载荷列阵

(241)Pe[Px1Py1Pz1Px2Py2Pz2Px8Py8Pz8]T （4-195）

5其中，节点外载Pz6Pz7F110N；支反力为Px1Rx1，Py1Ry1，Pz1Rz1，

Px4Rx4，Py4Ry4，Pz4Rz4，Px5Rx5，Py5Ry5，Pz5Rz5，Px8Rx8，Py8Ry8，Pz8Rz8；其余节点载荷分量为零。

（2）计算单元的刚度矩阵(以国际标准单位)

首先在MATLAB环境下，输入弹性模量E和泊松比NU，然后针对题中单元节点坐标，调用函数Hexahedral3D8Node_Stiffness，就可以得到单元的刚度矩阵k1 (24×24)。

>>E=1.0e10; >>NU=0.25; >>lx=0.2; >>ly=0.8; >>lz=0.6;

>>k1=Hexahedral3D8Node_Stiffness(E,NU,lx,0,0,lx,ly,0,0,ly,0,0,0,0,lx,0,lz,lx,ly,lz,0,ly,lz,0,0,lz);

（3） 建立整体刚度方程

由于该结构共有8个节点，则总共的自由度数为24，因此，结构总的刚度矩阵为K K(24×24)，先对KK清零，然后调用函数Hexahedral3D8Node_Assembly进行刚度矩阵的组装。由于本题中只用了一个单元，因此总体刚度矩阵KK与单元刚度k1相同，此处不再列出，调用函数的过程如下：

>>KK=zeros(24,24);

>>KK=Hexahedral3D8Node_Assembly(KK,k1,1,2,3,4,5,6,7,8);

（4） 边界条件的处理及刚度方程求解

由图4-23(b)可以看出，节点1，4，5和8的3个方向的位移将为零，即

u1v1w10,u4v4w40,u5v5w50,u8v8w80。因此，将针对节点2，

3，6和7的位移进行求解，这4个节点的位移将对应KK矩阵中的4～9行、4～9列，4～9行、16～21列，16～21行、4～9列，以及16～21行、16～21列，则需从KK(24×24)中提取相应行和列的数据，置给k，然后生成对应的载荷列阵p，再采用高斯消去法进行求解，即MATLAB中的反斜线符号“\\”求解。

>> k=[KK(4:9,4:9),KK(4:9,16:21);KK(16:21,4:9),KK(16:21,16:21)]; >> p=[0;0;0;0;0;0;0;0;-1e5;0;0;-1e5]; >> u=k\\p

u = 1.0e-003 *

0.0223 -0.2769 -0.6728 -0.0223 -0.2769 -0.6728 [将列排成行] -0.0129 0.3108 -0.7774 0.0129 0.3108 -0.7774 [将列排成行]

由此可以看出，所求得的结果为（单位为m）：

u20.022 3,v2-0.276 9,w2-0.672 8, u30.022 3,v30.276 9,w30.672 8,

u60.012 9,v60.310 8,w60.777 4, u70.012 9,v70.310 8,w70.777 4,

（5）支反力的计算

由方程KqP可知，在得到整个结构的节点位移后，由原整体刚度方程就可以计算出对应的支反力。先将上面得到的位移结果与位移边界条件的节点位移进行组合(注意位置关系)，可以得到整体的位移列阵U(24×1)，再代回原整体刚度方程，计算出所有的节点力P(24×1)，按式的对应关系就可以找到对应的支反力。

>> U=[0;0;0;u(1:6);0;0;0;0;0;0;u(7:12);0;0;0]; >> P=KK*U P = 1.0e+005 *

-0.2509 1.3333 0.6938 -0.0000 0.0000 0.0000 [将列排成行] 0.0000 0.0000 -0.0000 0.2509 1.3333 0.6938 [将列排成行] 0.3455 -1.3333 0.3062 -0.0000 -0.0000 -1.0000 [将列排成行] -0.0000 0.0000 -1.0000 -0.3455 -1.3333 0.3062 [将列排成行]

由式（4-195）的对应关系，可以得到对应的支反力为（单位为N）：

Rx10.250 9,Ry11.333 3,Rz10.693 8, Rx40.250 9,Ry41.333 3,Rz40.693 8,

Rx50.345 5,Ry51.333 3,Rz5 0.306 2,

Rx80.345 5,Ry81.333 3,Rz8 0.306 2,

（6）各单元的应力计算

先从整体位移列阵U(24×1)中提取出单元的位移列阵，然后，调用计算单元应力的函数Hexahedral3D8Node_Stress，就可以得到各个单元的应力分量。

>> u1=U(1:24);

>>stress1=Hexahedral3D8Node_Stress(E,NU,lx,0,0,lx,ly,0,0,ly,0,0,0,0,lx,0,lz,lx,ly,lz,0,ly,lz,0,0,lz,u1) stress1 = 1.0e+006 *

0.0197 0.0000 -0.8673 -0.0000 -1.6667 -0.0000 [将列排成行]

可以看出：计算得到的单元1的中心点的应力分量为

x19 700Pa,y0Pa, z867 300Pa,xy0Pa,yz1 666 700Pa,zx0Pa