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

【MATLAB算例】4.8.2(1) 基于8节点面体单元。基于MATLAB平台，用一个空间8节点六面体单元计算各个节点位移、支座反力以及单元的应力。 (a) 问题描述 (b) 有限元分析模型 图4-23 右端部受集中力作用的空间块体 解答：对该问题进行有限元分析的过程如下。 （1）结构的离散化与编号 将结构离散为一个8节点六面体单元，节点编号如图4-23(b)所示，节点的几何坐标见表4-13。 表4-13 节点的坐标 节点 节点坐标/m x y z 1 0.2 0 0 2 0.2 0.8 0 3 0 0.8 0 4 0 0 0 5 0.2 0 0.6 6 0.2 0.8 0.6 7 0 0.8 0.6 8 0 0 0.6 节点位移列阵 （4-194） 总的节点载荷列阵 （4-195） 其中，节点外载；支反力为，，，，，，，，，，，；其余节点载荷分量为零。 （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个方向的位移将为零，即。因此，将针对节点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）： （5）支反力的计算 由方程可知，在得到整个结构的节点位移后，由原整体刚度方程就可以计算出对应的支反力。先将上面得到的位移结果与位移边界条件的节点位移进行组合(注意位置关系)，可以得到整体的位移列阵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 [将列