计算力学课程设计模板改副本概要.docx
《计算力学课程设计模板改副本概要.docx》由会员分享，可在线阅读，更多相关《计算力学课程设计模板改副本概要.docx（13页珍藏版）》请在冰豆网上搜索。
计算力学课程设计模板改副本概要
本科课程设计
计算力学课程设计说明书
专业班级
姓名
学号
指导教师
完成日期
目录
1.引言
有限元法是一种高效能、常用的计算方法。
有限元法在早期是以变分原理为基础发展起来的，所以它广泛地应用于以拉普拉斯方程和泊松方程所描述的各类物理场中（这类场与泛函的极值问题有着紧密的联系）。
它最早应用于结构力学，后来随着计算机的发展慢慢用于流体力学的数值模拟。
在有限元方法中，把计算域离散剖分为有限个互不重叠且相互连接的单元，在每个单元内选择基函数，用单元基函数的线形组合来逼近单元中的真解，整个计算域上总体的基函数可以看为由每个单元基函数组成的，则整个计算域内的解可以看作是由所有单元上的近似解构成。
在河道数值模拟中，常见的有限元计算方法是由变分法和加权余量法发展而来的里兹法和伽辽金法、最小二乘法等。
根据所采用的权函数和插值函数的不同，有限元方法也分为多种计算格式。
从权函数的选择来说，有配置法、矩量法、最小二乘法和伽辽金法，从计算单元网格的形状来划分，有三角形网格、四边形网格和多边形网格，从插值函数的精度来划分，又分为线性插值函数和高次插值函数等。
不同的组合同样构成不同的有限元计算格式。
对于权函数，伽辽金（Galerkin）法是将权函数取为逼近函数中的基函数；最小二乘法是令权函数等于余量本身，而内积的极小值则为对代求系数的平方误差最小；在配置法中，先在计算域内选取N个配置点。
令近似解在选定的N个配置点上严格满足微分方程，即在配置点上令方程余量为0。
插值函数一般由不同次幂的多项式组成，但也有采用三角函数或指数函数组成的乘积表示，但最常用的多项式插值函数。
有限元插值函数分为两大类，一类只要求插值多项式本身在插值点取已知值，称为拉格朗日（Lagrange）多项式插值；另一种不仅要求插值多项式本身，还要求它的导数值在插值点取已知值，称为哈密特（Hermite）多项式插值。
单元坐标有笛卡尔直角坐标系和无因次自然坐标，有对称和不对称等。
常采用的无因次坐标是一种局部坐标系，它的定义取决于单元的几何形状，一维看作长度比，二维看作面积比，三维看作体积比。
在二维有限元中，三角形单元应用的最早，近来四边形等参元的应用也越来越广。
对于二维三角形和四边形电源单元，常采用的插值函数为有Lagrange插值直角坐标系中的线性插值函数及二阶或更高阶插值函数、面积坐标系中的线性插值函数、二阶或更高阶插值函数等。
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
1
2
3
4
5
6
7
8
0
0
0
0.4
0
0
0
2.0
0
0.4
2.0
0
0
0
1.0
0.4
0
1.0
0
2.0
1.0
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:
12,19:
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×10-3
＝-0.0295×10-3
＝0.0500×10-3
＝-0.2769×10-3
＝-0.2889×10-3
＝0.0994×10-3
＝0.1002×10-3
＝-0.0430×10-3
＝0.0380×10-3
＝-0.2761×10-3
＝-0.3000×10-3
（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:
3）;U（10:
12）;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:
21）;U（22:
24）;U（10:
12）];
>>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:
3）;U（10:
12）;U（16:
21）];
>>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;1y3z3;1y4z4];
mbeta3=[1y1z1;1y2z2;1y4z4];
mbeta4=[1y1z1;1y2z2;1y3z3];
mgamma1=[1x2z2;1x3z3;1x4z4];
mgamma2=[1x1z1;1x3z3;1x4z4];
mgamma3=[1x1z1;1x2z2;1x4z4];
mgamma4=[1x1z1;1x2z2;1x3z3];
mdelta1=[1x2y2;1x3y3;1x4y4];
mdelta2=[1x1y1;1x3y3;1x4y4];
mdelta3=[1x1y1;1x2y2;1x4y4];
mdelta4=[1x1y1;1x2y2;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:
12
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）
%该函数计算单元的应力
%输入弹性模量E，泊松比NU
%输入4个节点i、j、m、n的坐标xi,yi,zi,xj,yj,zj,xm,ym,zm,xn,yn,zn
%输入单元的位移列阵u（12X1）
%输出单元的应力stress（6X1）
%由于它为常应力单元，则单元的应力分量为Sx,Sy,Sz,Sxy,Syz,Szx
%---------------------------------------------------------------
xyz=[1x1y1z1;1x2y2z2;1x3y3z3;1x4y4z4];
V=det（xyz）/6;
mbeta1=[1y2z2;1y3z3;1y4z4];
mbeta2=[1y1z1;1y3z3;1y4z4];
mbeta3=[1y1z1;1y2z2;1y4z4];
mbeta4=[1y1z1;1y2z2;1y3z3];
mgamma1=[1x2z2;1x3z3;1x4z4];
mgamma2=[1x1z1;1x3z3;1x4z4];
mgamma3=[1x1z1;1x2z2;1x4z4];
mgamma4=[1x1z1;1x2z2;1x3z3];
mdelta1=[1x2y2;1x3y3;1x4y4];
mdelta2=[1x1y1;1x3y3;1x4y4];
mdelta3=[1x1y1;1x2y2;1x4y4];
mdelta4=[1x1y1;1x2y2;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
升级会员 