有限元计算结构力学大作业.docx
- 文档编号:10725052
- 上传时间:2023-05-27
- 格式:DOCX
- 页数:26
- 大小:1.29MB
有限元计算结构力学大作业.docx
《有限元计算结构力学大作业.docx》由会员分享,可在线阅读,更多相关《有限元计算结构力学大作业.docx(26页珍藏版)》请在冰点文库上搜索。
有限元计算结构力学大作业
本页仅作为文档页封面,使用时可以删除Thisdocumentisforreferenceonly-rar21year.March
有限元-计算结构力学-大作业
SHANGHAIJIAOTONGUNIVERSITY
平面应力问题解的Matlab实现
姓名:
heiya168
学号:
帆哥
班级:
指导老师:
1绪论
有限元方法(finiteelementmethod),是求取复杂微分方程近似解的一种非常有效的工具,是现代数字化科技的一种重要基础性原理。
将它用于在科学研究中,可成为探究物质客观规律的先进手段。
将它应用于工程技术中,可成为工程设计和分析的可靠工具。
弹性体在载荷作用下,其基本方程可写成以下的三类方程和两种边界条件。
平衡方程——应力与外载荷的关系;几何方程——应变位移关系;物理方程——应力应变关系;力的边界条件;几何边界条件。
应用最小位能原理,并利用上述关系,最终建立由刚度方程,节点位移和等效节点载荷所构成的求解方程。
带入边界条件求解方程,就可以得出弹性力学问题的一般性解答。
本次大作业基于有限元方法的基本原理,使用Matlab这一平台,针对平面应力问题,采用四节点四边形单元编写了求解单元节点位移的程序。
主要内容包括:
1)介绍有限元的基本原理;2)编程基本思路及流程介绍;3)程序原理及说明;4)具体算例这四个部分。
2平面问题的四节点四边形单元
2.1单元的构造
(1)单元的几何和节点描述
平面4节点矩形单元如图4-6所示,单元的节点位移共有8个自由度。
节点的编号为1、2、3、4,各自的位置坐标为(xi,yi),i=1,2,3,4,各个节点的位移(分别沿x方向和y方向)为(ui,vi),i=1,2,3,4。
图2.1平面4节点矩形单元
若采用无量纲坐标
(2.1)
则单元4个节点的几何位置为
(2.2)
将所有节点上的位移组成一个列阵,记作
;同样,将所有节点上的各个力也组成一个列阵,记作
,那么
(2.3)
若该单元承受分布外载,可以将其等效到节点上,也可以表示为如式(2.3)所示的节点力。
利用函数插值、几何方程、物理方程以及势能计算公式,可以将单元的所有力学参量用节点位移列阵及相关的插值函数
来表示;下面进行具体的推导。
(2)单元位移场的表达
从图2-1可以看出,节点条件共有8个,即x方向4个
,y方向4个
,因此,x和y方向的位移场可以各有4个待定系数,即取以下多项式作为单元的位移场模式
(2.4)
它们是具有完全一次项的非完全二次项,以上两式中右端的第四项是考虑到x方向和y方向的对称性而取的,除此外xy项还有个重要特点,就是“双线性”,当x或y不变时,沿y或x方向位移函数呈线性变化,这与前面的线性项最为相容,而
或
项是二次曲线变化的。
因此,未选
或
项。
由节点条件,在x=xi,y=yi处,有
(2.5)
将式(2.5)代入式(2.4)中,可以求解出待定系数a0,…,a3和b0,…,b3,然后代回式(4-52)中,经整理后有
(2.6)
其中
(2.7)
如以无量纲坐标系(2.1)来表达,式(2.7)可以写成
(2.8)
将式(2.6)写成矩阵形式,有
(2.9)
其中,
为该单元的形状函数矩阵。
(3)单元应变场的表达
由弹性力学平面问题的几何方程(矩阵形式),有单元应变的表达
(2.10)
其中几何矩阵B(x,y)为
=
(2.11)
式(4-59)中的子矩阵
为
(2.12)
(4)单元应力场的表达
由弹性力学中平面问题的物理方程,可得到单元的应力表达式
(2.13)
(5)单元应力场的表达
以上已将单元的三大基本变量
用基于节点位移列阵来进行表达,见式(2.9)、式(2.10)及式(2.13);将其代入单元的势能表达式中,有
(2.14)
其中
是4节点矩形单元的刚度矩阵。
将单元的势能对节点位移
取一阶极值,可得到单元的刚度方程
(2.15)
2.2等参变换
由前面的单元构造过程可以看出,一个单元的关键就是计算它的刚度矩阵,而由刚度矩阵的构成可知要实现两个坐标系中单元刚度矩阵的变换,必须计算两个坐标系之间的三种映射关系:
坐标映射
(2.16)
偏导数映射
(2.17)
面积映射
(2.18)
图2.2矩形单元映射为任意四边形单元
(1)两个坐标系之间的函数映射
设如图4-17所示的两个坐标系的坐标映射关系为
(2.19)
x和y方向上可以分别写出各包含有4个待定系数的多项式,即
(2.20)
其中待定系数a0,…,a3和b0,…,b3可由节点映射条件(2.19)来唯一确定。
对照前面4节点矩形单元的单元位移函数式(2.5),映射函数式(2.20)具有完全相同的形式,同样,将求出的待定系数再代回式(2.20)中,重写该式为
(2.21)
其中
(2.22)
(2.23)
这就可以实现两个坐标系间的映射。
(2)两个坐标系之间的偏导数映射
对物理坐标系(x,y)中的任意一个函数Φ(x,y),求它的偏导数,有
(2.24)
则偏导数的变换关系为
(2.25)
写成矩阵形式,有
(2.26)
其中
(2.27)
两个坐标系的偏导数映射关系
(2.28)
(3)两个坐标系之间的面积元映射
如图2-2所示,在物理坐标系(x,y)中,由dξ和dη所围成的微小平行四边形,其面积为
(2.29)
由于dξ和dη在物理坐标系(x,y)中的分量为
(2.30)
其中i和j分别为物理坐标系(x,y)中的x方向和y方向的单位向量。
由(2.29)式,则有面积微元的变换计算
(2.31)
1.1边界条件的处理——置“1”法
设边界条件为第r个自由度的位移为零,即
;可置整体刚度矩阵中所对应对角元素位置的
,而该行和该列的其它元素为零,即
,同时也置对应的载荷元素则
,则
(2.32)
进行以上设置后,这时方程(2.32)应等价于原方程加上了边界条件
,下面考察这种等价性,就(2.32)中的第r行,有
(2.33)
由于置
,则有
即为所要求的位移边界条件。
而式(5-53)中除第r行外,其它各行在对应于r列的位置上都置了“0”,这相当于考虑了
的影响,除此之外其余各项的影响不变;这恰好就是原方程加上了边界条件
的影响。
3有限元分析流程
3.1程序原理和流程
该程序的特点如下:
问题类型:
可用于弹性力学平面应力问题和平面应变问题的分析
单元类型:
四节点四边形单元
材料性质:
单一的均匀的弹性材料
节点信息:
节点信息文件node.txt
单元信息:
单元信息文件element.txt
载荷信息:
等效节点力文件force.txt
约束信息:
节点约束信息文件constrain.txt
结果文件:
输出节点位移文件node_displace.txt
程序流程图如下:
图3.1程序流程图
1.1使用的函数
Stiffness(E,NU,h,xi,yi,xj,yj,xm,ym,xp,yp):
计算单元的刚度矩阵(具体说明见附录)
Assembly(KK,k,i,j,m):
进行单元刚度矩阵的组装(具体说明见附录)
1.2文件管理
主程序源文件:
FEM.m
输入参数文件:
node.txt(节点信息文件)
element.txt(单元信息文件)
constrain.txt(位移边界条件文件)
force.txt(载荷状况文件)
输出位移文件:
node_displace.txt(节点位移文件)
3.2数据文件格式
4算例——开方孔的矩形板拉伸分析
4.1问题的具体参数与载荷
(a)对边拉伸的带孔矩形薄板(b)1/4模型的单元划分
图4-1受均匀荷载的作用的带孔矩形薄板薄板及有限元分析模型
如图4-1(a)所示的矩形带孔薄板对边受均匀荷载的作用,该结构在边界上受正向分布拉力0.3MPa。
若取板厚0.02m,弹性模量2.1E11,泊松比μ=0.3,采用FEM.m程序,和ANSYS分别分析该问题,最后给出各节点位移值。
4.2Matlab程序计算
(1)结构的离散化与编号
该薄板的荷载和几何形状关于x轴和轴对称,故可只取结构的1/4作为计算模型。
将此模型化分为四个全等的直角三角形单元,单元编号和节点编号如图4-1(b)所示。
(2)输入节点信息node.txt
(3)输入单元信息element.txt
(4)输入载荷信息force.txt
(5)输入载荷信息constrain.txt
(6)运行主函数FEM.m,输入E=2.1E11,NU-0.3,h=0.02。
(7)输出各节点位移
4.3ANSYS建模计算
选用plane42板单元,建立如下图所示模型,在2、3节点加载水平方向的载荷10000N.
图4-2受均匀荷载的作用的带孔矩形薄板薄板有限元分析模型
分析计算所得的最后结果为:
4.4误差分析
节点位移
2X
2Y
3X
3Y
4X
4Y
2SUM
3SUM
4SUM
Matlab
7.65E-06
2.53E-07
7.33E-06
-1.43E-06
1.87E-06
2.025E-07
7.66E-06
7.469E-06
1.88E-06
ANSYS
7.54E-06
-8.00E-07
7.81E-06
-2.58E-06
1.68E-06
2.45E-07
7.58E-06
8.22E-06
1.70E-06
误差
1.02%
9.18%
9.94%
表4-1Matlab程序与Ansys结果对比表
误差范围在10%以内,还是比较大的,导致误差的原因可能有以下几个方面:
1)刚度矩阵计算时,使用的积分方法不同,本程序中使用的是两点高斯积分。
2)刚度矩阵求解是,使用的算法不同,本程序使用的是高斯消去法。
3)约束处理时使用的方法不同,本程序使用的是置一法。
5总结
用了两周的时间,完成了这次大作业。
在做这次大作业之前,我对于有限元算法的了解还很肤浅,只是模模糊糊的能够记得刚度阵建立的过程,
矩阵表示的公式,对于其中具体的物理意义,各个量之间相互关系的表示不太清楚,在做完这次大作业后,对于四边形单元刚度阵的建立过程有了清晰深入的认识和掌握。
关于等参变化的认识,是这次大作业中的难点之一。
对于和位移差值函数相同的坐标差值函数,如何理解其中“坐标差值”物理意义,在差值完成,刚度矩阵建立后,刚度矩阵的位移模式是局部坐标系中的,还是整体坐标系中的位移,这个问题思考了很久,也和很多同学进行了争论和探讨。
最终得到了一致的结论,收获很多。
在完成大作业之前,我对于Matlab编程还是一窍不通,在编程的过程中,很多好的想法无法用计算机语言来实现是一件很痛苦的事情,但是在同学的帮助下和参考资料的借鉴下,完成了这次编程。
其中,刚度矩阵和组装矩阵是自己独立完成的,主程序借鉴了清华大学曾攀教授编写的《有限元分析基础教程》中平面三节点单元的主程序。
写这个总结的时候,也就是大作业完成的时候,其中还有很多不足的地方有待改进,比如两点高斯积分精确度较差的问题,比如无法直接输入面载荷和体载荷的问题。
参考文献
[1]王勖成.有限单元法[M].清华大学出版社.
[2]曾攀.有限元分析基础教程[M].清华大学出版社.
[3]陈杰.MATLAB宝典[M].电子工业出版社
附录
functionK=Stiffness(E,NU,h,xi,yi,xj,yj,xm,ym,xp,yp)
%该函数计算单元的刚度矩阵
%输入弹性模量E,泊松比NU,厚度h
%输入4个节点i、j、m、p的坐标xi,yi,xj,yj,xm,ym,xp,yp
%输出单元刚度矩阵k(8X8)
%---------------------------------------------------------------
symsst;
N1=1/4*(1-s)*(1-t);
N2=1/4*(1+s)*(1-t);
N3=1/4*(1+s)*(1+t);
N4=1/4*(1-s)*(1+t);%坐标变换函数,位移插值函数
X=N1*xi+N2*xj+N3*xm+N4*xp;
Y=N1*yi+N2*yj+N3*ym+N4*yp;%坐标变换
J11=diff(X,s);
J12=diff(Y,s);
J21=diff(X,t);
J22=diff(Y,t);
J=[J11J12;J21J22];%J矩阵的建立
H=1/det(J)*[J22-J1200;00-J21J11;-J21J11J22-J21];%B=HQ,H矩阵建立
q11=diff(N1,s);
q21=diff(N1,t);
q32=diff(N1,s);
q42=diff(N1,t);
Q1=[q110;q210;0q32;0q42];
q13=diff(N2,s);
q23=diff(N2,t);
q34=diff(N2,s);
q44=diff(N2,t);
Q2=[q130;q230;0q34;0q44];
q15=diff(N3,s);
q25=diff(N3,t);
q36=diff(N3,s);
q46=diff(N3,t);
Q3=[q150;q250;0q36;0q46];
q17=diff(N4,s);
q27=diff(N4,t);
q38=diff(N4,s);
q48=diff(N4,t);
Q4=[q170;q270;0q38;0q48];
Q=[Q1Q2Q3Q4];%Q矩阵的建立
B=H*Q;
D=(E/(1-NU^2))*[1NU0;NU10;00(1-NU)/2];
A=transpose(B)*D*B;
k=A*det(J)*h;
s=1/3*3^0.5;t=1/3*3^0.5;
k1=eval(k);
s=-1/3*3^0.5;t=1/3*3^0.5;
k2=eval(k);
s=1/3*3^0.5;t=-1/3*3^0.5;
k3=eval(k);
s=-1/3*3^0.5;t=-1/3*3^0.5;
k4=eval(k);
K=k1+k2+k3+k4;
functionz=Assembly(KK,k,i,j,m,p)
%该函数进行单元刚度矩阵的组装
%输入单元刚度矩阵k,单元的节点编号i、j、m、p
%输出整体刚度矩阵KK
%---------------------------------------------------------------
DOF
(1)=2*i-1;
DOF
(2)=2*i;
DOF(3)=2*j-1;
DOF(4)=2*j;
DOF(5)=2*m-1;
DOF(6)=2*m;
DOF(7)=2*p-1;
DOF(8)=2*p;
forn1=1:
8
forn2=1:
8
KK(DOF(n1),DOF(n2))=KK(DOF(n1),DOF(n2))+k(n1,n2);
end
end
z=KK;
%%%%%%%%%%%%FEM.m%%mainprogrambegin%%%%%%%%%%
%%该程序为计算平面应力问题的主程序%%
%%输入节点,单元,载荷,约束信息,调用Stiffness和Assembly%%
%%求得总刚,并解得节点位移%%
clear;
%读入节点,单元,载荷,约束信息
load('F:
\FEM\node.txt')
load('F:
\FEM\element.txt')
load('F:
\FEM\constrain.txt')
load('F:
\FEM\force.txt')
%确定节点、单元个数
[nnode,ntmp]=size(node);
[nelem,etmp]=size(element);
[nforce,ftmp]=size(force);
[nconstrain,ctmp]=size(constrain);
%预先设定总体刚度矩阵、节点力向量、节点位移向量
KKG=zeros(2*nnode);
FFG=zeros(2*nnode,1);
UUG=zeros(2*nnode,1);
k=zeros(8,8);%单元刚度矩阵
%给出相应材料及计算参数
E=input('弹性模量E=');
NU=input('泊松比NU=');
t=input('板厚=');
%单元循环形成总体刚度矩阵
fori=1:
nelem
K=Stiffness(E,NU,t,node(element(i,1),2),node(element(i,1),3),node(element(i,2),2),node(element(i,2),3),node(element(i,3),2),node(element(i,3),3),node(element(i,4),2),node(element(i,4),3));
KKG=Assembly(KKG,K,element(i,1),element(i,2),element(i,3),element(i,4));
end
%载荷的处理
fori=1:
nforce
m=force(i,1);
n=force(i,2);
FFG(2*(m-1)+n)=force(i,3);
end
%采用置”1”法处理边界条件
fori=1:
nconstrain
m=constrain(i,1);
n=constrain(i,2);
UUG(2*(m-1)+n)=constrain(i,3);
KKG(2*(m-1)+n,:
)=0;
KKG(:
2*(m-1)+n)=0;
KKG(2*(m-1)+n,2*(m-1)+n)=1;
FFG(2*(m-1)+n)=0;
end
%求解节点位移
UUG=KKG\FFG;
%输出节点位移值
fid=fopen('F:
\FEM\node_displace.txt','w');
fprintf(fid,'\n%s\n','--------------------NODEDISPLACEMENT---------------------------');
fprintf(fid,'\n%s\n','NODEX-dispY-disp');
fori=1:
nnode
fprintf(fid,'%d%18.10f%18.10f\n',node(i,1),UUG(2*(i-1)+1),UUG(2*(i-1)+2));
end
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 有限元 计算 结构 力学 作业