计算力学课程设计模板改副本概要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