《有限元基础教程》_【MATLAB算例】4.8.docx

【算8(基于8节点六面体单元的空间块体分析(Hexahedral3D8Node)如图4-23所示的一个空间块体，在右端部受两个集中力F作用，其中的参数为：

E=1×101?Pa,μ=0.25,t=0.2m,F=1×10?N。基于MATLAB平台，用一个空间8节点六面体单元计算各个节点位移、支座反力以及单元的应力。

(a)问题描述(b)有限元分析模型图4-23右端部受集中力作用的空间块体

解答：对该问题进行有限元分析的过程如下。

(1)结构的离散化与编号

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

表4-13。

节点位移列阵

总的节点载荷列阵

表4-13节点的坐标

节点

X

节点坐标/my

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)

P=[P?P?P?P?P?P2…P?Py?P?8](4-195)

其中，节点外载P??=P27=-F=-1×10?N;支反力为P?=R,P?=R?,P?=R?,

Px?=Rx?,Py?=Ry?,P??=R?4,P?=R?,P,?=R?,P?5=R?5,P?=R?,Py?=Ry?,

P?=R?;其余节点载荷分量为零。

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

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

E=1.0e10;NU=0.25;lx=0.2;

ly=0.8;lz=0.6;

kl=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个方向的位移将为零，即

u?=v?=w?=0,u=v?=w?=0,u?=v?=w?=0,ug=vg=w?=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=klp

u=..0e2-0203*-0.2769-0.6728-0.0223-0.2769-0.6728[将列排成行]

-0012003108-077740012003108-07774「将列排成行1

[将列排成]

由此可以看出，所求得的结果为(单位为m):

u?=0.0223,v?=-0.2769,w?=-0.6728,u?=-0.0223,v?=-0.2769,w?=-0.6728,u?=-0.0129,v?=0.3108,w?=-0.7774,u?=0.0129,v?=0.3108,w?=-0.7774,

(5)支反力的计算

由方程Kq=P可知，在得到整个结构的节点位移后，由原整体刚度方程就可以计算出

对应的支反力。先将上面得到的位移结果与位移边界条件的节点位移进行组合(注意位置关系),可以