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

【MATLAB算例】4.8.2(1)基于8节点六面体单元的空间块体分析(Hexahedral3D8Node)如图4-23所示的一个空间块体，在右端部受两个集中力F作用，其中的参数为：105110Pa,=0.25,=0.2m,=110NEtF。基于MATLAB平台，用一个空间8节点六面体单元计算各个节点位移、支座反力以及单元的应力。(a)问题描述(b)有限元分析模型图4-23右端部受集中力作用的空间块体解答：对该问题进行有限元分析的过程如下。（1）结构的离散化与编号将结构离散为一个8节点六面体单元，节点编号如图4-23(b)所示，节点的几何坐标见表4-13。表4-13节点的坐标节点节点坐标/mxyz10.20020.20.80300.80400050.200.660.20.80.6700.80.68000.6节点位移列阵T111222888(241)[]euvwuvwuvwq（4-194）总的节点载荷列阵T111222888(241)[]exyzxyzxyzPPPPPPPPPP（4-195）其中，节点外载567110zzPPFN；支反力为11xxPR，11yyPR，11zzPR，44xxPR，44yyPR，44zzPR，55xxPR，55yyPR，55zzPR，88xxPR，88yyPR，88zzPR；其余节点载荷分量为零。（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，因此，结构总的刚度矩阵为KK(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个方向的位移将为零，即1114445558880,0,0,0uvwuvwuvwuvw。因此，将针对节点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\pu=1.0e-003*0.0223-0.2769-0.6728-0.0223-0.2769-0.6728[将列排成行]-0.01290.3108-0.77740.01290.3108-0.7774[将列排成行]由此可以看出，所求得的结果为（单位为m）：2220.0223,-0.2769,-0.6728,uvw3330.0223,0.2769,0.6728,uvw6660.0129,0.3108,0.7774,uvw7770.0129,0.3108,0.7774,uvw（5）支反力的计算由方程KqP可知，在得到整个结构的节点位移后，由原整体刚度方程就可以计算出对应的支反力。先将上面得到的位移结果与位移边界条件的节点位移进行组合(注意位置关系)，可以得到整体的位移列阵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*UP=1.0e+005*-0.25091.33330.6938-0.00000.00000.0000[将列排成行]0.00000.0000-0.00000.25091.33330.6938[将列排成行]0.3455-1.33330.3062-0.0000-0.0000-1.0000[将列排成行]-0.00000.0000-1.0000-0.3455-1.33330.3062[将列排成行]由式（4-195）的对应关系，可以得到对应的支反力为（单位为N）：1110.2509,1.3333,0.6938,xyzRRR4440.2509,1.3333,0.6938,xyzRRR5550.3455,1.3333,0.3062,xyzRRR8880.3455,1.3333,0.3062,xyzRRR（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.01970.0000-0.8673-0.0000-1.6667-0.0000[将列排成行]可以看出：计算得到的单元1的中心点的应力分量为19700Pa,x0Pa,y867300Pa,z0Pa,xy1666700Pa,yz0Pazx