计算力学课程设计模板改副本概要Word下载.docx
《计算力学课程设计模板改副本概要Word下载.docx》由会员分享，可在线阅读，更多相关《计算力学课程设计模板改副本概要Word下载.docx（13页珍藏版）》请在冰豆网上搜索。
2.基本力学模型
如图所示的一个块体，在右端面上端点受集中力F作用。
基于MATLAB平台，计算各个节点位移，支座反力以及单元的应力。
取相关的参数E=1.5e10Pa,u=0.25,F=1.5e5N。
对该问题进行有限元分析的过程如下：
图一
（1）结构的离散化与编号
将结构离散为5个4节点四面体单元，单元编号和坐标如上图所示，连接关系见表一，节点的坐标见表二。
表一单元连接关系
单元号
节点号
1
1426
2
1437
3
6751
4
6784
5
1467
表二节点的坐标
节点
节点坐标/m
x
y
z
6
7
8
0.4
2.0
1.0
节点位移列阵
节点外载列阵
其中
约束的支反力列阵
总的节点载荷列阵
（2）计算各单元的刚度矩阵（以国际标准单位）
首先在MATLAB环境下，输入弹性模量E、泊松比NU，然后针对单元1和单元2，分别5次调用函数Tetrahedron3D4Node_Stiffness，就可以得到单元的刚度矩阵k1（6×
6）～k5（6×
6）。
>
E=2.5e10;
NU=0.25;
k1=Tetrahedron3D4Node_Stiffness（E,NU,0,0,0,0.4,2.0,0,0.4,0,0,0.4,0,1.0）;
k2=Tetrahedron3D4Node_Stiffness（E,NU,0,0,0,0.4,2.0,0,0,2.0,0,0,2.0,1.0）;
k3=Tetrahedron3D4Node_Stiffness（E,NU,0.4,0,1.0,0,2.0,1.0,0,0,1.0,0,0,0）;
k4=Tetrahedron3D4Node_Stiffness（E,NU,0.4,0,1.0,0,2.0,1.0,0.4,2.0,1.0,0.4,2.0,0）;
k5=Tetrahedron3D4Node_Stiffness（E,NU,0,0,0,0.4,2.0,0,0.4,0,1.0,0,2.0,1.0）;
（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个方向的位移将为零，即
。
因此，将针对节点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:
24]）;
p=[0,0,0,0,0,0,0,0,-2e5,0,0,-2e5]'
u=k\p
u=1.0e-003*
0.0958-0.0295-0.27690.0994-0.0430-0.2761[将列排成行排量[]___________________________________________________________________________________________________________________________]
0.09870.0500-0.28890.10020.0380-0.3000[将列排成行排量[]___________________________________________________________________________________________________________________________]
所求得的位移结果见表三。
表三空间块体的节点位移计算结果
＝0.0958×
10-3
＝0.0987×
＝-0.0295×
＝0.0500×
＝-0.2769×
＝-0.2889×
＝0.0994×
＝0.1002×
＝-0.0430×
＝0.0380×
＝-0.2761×
＝-0.3000×
（5）支反力的计算
在得到整个结构的节点位移后，由原整体刚度方程就可以计算出对应的支反力；
先将上面得到的位移结果与位移边界条件的节点位移进行组合（注意位置关系），可以得到整体的位移列阵U（24×
1），再代回原整体刚度方程，计算出所有的节点力P（24×
1），按式（4-192）的对应关系就可以找到对应的支反力。
U=[zeros（6,1）;
u（[1:
6]）;
zeros（6,1）;
u（7:
12）];
P=KK*U
P=1.0e+005*
0.92304.07170.3908-1.04813.92831.2071[将列排成行排量[]___________________________________________________________________________________________________________________________]
-0.0000-0.00000.00000.00000-0.0000[将列排成行排量[]___________________________________________________________________________________________________________________________]
-1.1625-4.07171.29651.2876-3.92831.1056[将列排成行排量[]___________________________________________________________________________________________________________________________]
-0.00000.0000-2.00000.0000-0.0000-2.0000[将列排成行排量[]___________________________________________________________________________________________________________________________]
由式（4-193）的对应关系，可以得到对应的支反力见表四。
表四空间块体的支反力计算结果
（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.4,2.0,0,0.4,0,0,0.4,0,1.0,u1）
stress1=1.0e+006*
-0.2150-0.6449-0.21500.4972-1.38070[将列排成行排量[]___________________________________________________________________________________________________________________________]
u2=[U（1:
U（7:
9）;
U（19:
21）];
stress2=Tetrahedron3D4Node_Stress（E,NU,0,0,0,0.4,2.0,0,0,2.0,0,0,2.0,1.0,u2
stress2=1.0e+005*
0.0885-4.7076-4.16401.4161-5.89410.4868[将列排成行排量[]___________________________________________________________________________________________________________________________]
u3=[U（16:
21）;
U（13:
15）;
U（1:
3）];
stress3=Tetrahedron3D4Node_Stress（E,NU,0.4,0,1.0,0,2.0,1.0,0,0,1.0,0,0,0,u3）
stress3=1.0e+006*
0.25000.75010.25000.4936-1.44470[将列排成行排量[]___________________________________________________________________________________________________________________________]
u4=[U（16:
U（22:
24）;
stress4=Tetrahedron3D4Node_Stress（E,NU,0.4,0,1.0,0,2.0,1.0,0.4,2.0,1.0,0.4,2.0,0,u4）
stress4=1.0e+005*
0.66673.7041-4.86282.0178-6.8964-2.6756[将列排成行排量[]___________________________________________________________________________________________________________________________]
u5=[U（1:
stress5=Tetrahedron3D4Node_Stress（E,NU,0,0,0,0.4,2.0,0,0.4,0,1.0,0,2.0,1.0,u5）
stress5=1.0e+005*
-0.1936-0.0239-1.6540-6.6710-9.47741.5635[将列排成行排量[]___________________________________________________________________________________________________________________________]
各个单元应力分量的计算结果列在表五中。
表五空间块体的各个单元应力分量的计算结果
1号单元
2号单元
3号单元
4号单元
5号单元
3计算结果分析
在如图一的悬臂式梁结构模型中，当受到如图的集中荷载作用时，1,2,5,6节点的y方向上的支反力较大，这符合竖向荷载作用下的支承力实际情况;
同时，通过表五对空间块体的各个单元应力分量的计算结果可以看出5号单元的剪切力最大，容易发生剪切破坏，因此需要对其进行加固处理。
4．参考文献
[1]曾攀.有限元基础教程[M].北京:
高等教育出版社,2009,7
[2]王勖成.有限单元法[M].北京:
清华大学出版社,2003,
[3]张国瑞有限元法北京机械工业出版社1991
[4]杨帆宋小军悬臂梁的有限元分析江西理工大学2008
[5]江见鲸陆新征混凝土结构有限元分析北京：
清华大学出版社2005
附录（程序）
1.Tetrahedron3D4Node_Stiffness函数
functionk=Tetrahedron3D4Node_Stiffness（E,NU,x1,y1,z1,x2,y2,z2,x3,y3,z3,x4,y4,z4）
%该函数计算单元的刚度矩阵
%输入弹性模量E，泊松比NU
%输入4个节点i、j、m、n的坐标xi,yi,zi,xj,yj,zj,xm,ym,zm,xn,yn,zn
%输出单元刚度矩阵k（12X12）
%---------------------------------------------------------------
%以下数组与书中的对应关系
xyz=[1x1y1z1;
1x2y2z2;
1x3y3z3;
1x4y4z4];
V=det（xyz）/6;
mbeta1=[1y2z2;
1y3z3;
1y4z4];
mbeta2=[1y1z1;
mbeta3=[1y1z1;
1y2z2;
mbeta4=[1y1z1;
1y3z3];
mgamma1=[1x2z2;
1x3z3;
1x4z4];
mgamma2=[1x1z1;
mgamma3=[1x1z1;
1x2z2;
mgamma4=[1x1z1;
1x3z3];
mdelta1=[1x2y2;
1x3y3;
1x4y4];
mdelta2=[1x1y1;
mdelta3=[1x1y1;
1x2y2;
mdelta4=[1x1y1;
1x3y3];
beta1=-1*det（mbeta1）;
beta2=det（mbeta2）;
beta3=-1*det（mbeta3）;
beta4=det（mbeta4）;
gamma1=det（mgamma1）;
gamma2=-1*det（mgamma2）;
gamma3=det（mgamma3）;
gamma4=-1*det（mgamma4）;
delta1=-1*det（mdelta1）;
delta2=det（mdelta2）;
delta3=-1*det（mdelta3）;
delta4=det（mdelta4）;
B1=[beta100;
0gamma10;
00delta1;
gamma1beta10;
0delta1gamma1;
delta10beta1];
B2=[beta200;
0gamma20;
00delta2;
gamma2beta20;
0delta2gamma2;
delta20beta2];
B3=[beta300;
0gamma30;
00delta3;
gamma3beta30;
0delta3gamma3;
delta30beta3];
B4=[beta400;
0gamma40;
00delta4;
gamma4beta40;
0delta4gamma4;
delta40beta4];
B=[B1B2B3B4]/（6*V）;
D=（E/（（1+NU）*（1-2*NU）））*[1-NUNUNU000;
NU1-NUNU000;
NUNU1-NU000;
000（1-2*NU）/200;
0000（1-2*NU）/20;
00000（1-2*NU）/2];
k=abs（V）*B'
*D*B;
2.Tetrahedron3D4Node_Assembly函数
functiony=Tetrahedron3D4Node_Assembly（KK,k,i,j,m,n）
%该函数进行单元刚度矩阵的组装
%输入单元刚度矩阵k
%输入单元的节点编号i、j、m、n
%输出整体刚度矩阵KK
DOF=[3*i-2:
3*i,3*j-2:
3*j,3*m-2:
3*m,3*n-2:
3*n];
forn1=1:
12
forn2=1:
KK（DOF（n1）,DOF（n2））=KK（DOF（n1）,DOF（n2））+k（n1,n2）;
end
end
y=KK;
3.Tetrahedron3D4Node_Stress函数
functiony=Tetrahedron3D4Node_Stress（E,NU,x1,y1,z1,x2,y2,z2,x3,y3,z3,x4,y4,z4,u）
%该函数计算单元的应力
%输入单元的位移列阵u（12X1）
%输出单元的应力stress（6X1）
%由于它为常应力单元，则单元的应力分量为Sx,Sy,Sz,Sxy,Syz,Szx
delta4=det（mdelta4
升级会员 