有限元课程设计.docx
- 文档编号:10043368
- 上传时间:2023-05-23
- 格式:DOCX
- 页数:26
- 大小:278.51KB
有限元课程设计.docx
《有限元课程设计.docx》由会员分享,可在线阅读,更多相关《有限元课程设计.docx(26页珍藏版)》请在冰点文库上搜索。
有限元课程设计
《有限元法》课程设计报告
班级:
10级工程力学2班
姓名:
刘海波
学号:
1002060221
指导教师:
谢献忠
湖南科技大学土木工程学院
2013年6月12日
(一)问题描述
题2:
图示薄板左边固定,右边受均布压力P=100Kn/m作用,板厚度为0.3cm;试采用如下方案,对其进行有限元分析,并对结果进行比较。
1)三节点常应变单元;(2个和200个单元)
2)四节点矩形单元;(1个和50个单元)
(二)有限元计算模型
因为该题荷载和边界都关于横向对称,故从中间取一半进行分析。
网格划分方案如下图所示,将作用在单元上的均布载荷按虚功等效的原则移置到节点上,成为等效节点载荷。
因取一半,故下方需在节点位置加上竖向链杆支座。
(三)基本理论
1.三角形单元的线性位移模式
单元应变矩阵
对于平面应力问题,弹性矩阵[D]为
总刚方程[K]{}={R}
2.矩形单元
位移模式
单元刚度矩阵
如果单元厚度t是常量,则
(四)计算结果及分析
两单元的计算结果及200单元的计算结果见Matlab中
每个单元的应变
2单元的应力δxδyΤxysigm=
1.0e+009*
0.0605-0.3760
0.24191.1627
-0.4013-0.2759
100单元的位移分析:
(放大1000倍后)
结果数据过多,见matlab
矩形单元也进行了分析
放大1000倍的变形位移ans=
-0.00000.0000
-0.05110.0112
-0.00000.0000
-0.0468-0.0000
-0.0000-0.0000
-0.0511-0.0112
(五)多方案计算比较
取一半分析时,是三角形单元,2单元的没有100个单元的精度高,网格划分越多,结果也越精确,矩形双线性单元比三角形单元具有更高的精度。
从定性分析来看,编程结果也基本符合实际物体受力。
(六)建议与体会
经过这些天的编程,我学到了很多关于matlab的技巧,对有限元也有了进一步的理解,当把所学的用程序表达出来时,虽然在所学过程中遇到了困难,花费了不少时间,当看到最后的编程成果时,内心还是充满了成就感的,这些天的编程使我对三角形单元有了熟练的掌握,但对矩形双线性单元还有些困惑,深感自己在有限元课程上还有很长的道路要走。
感谢老师的悉心教诲。
回顾这个学期所学这门课程,有限元是门很有趣,很实用的课程。
(七)源程序
(1)
%课程设计:
第一问:
2个单元的三角形单元,这里取一半来分析,然后对称
%生成原始数据
clc
clear
Gao=1.0;%板高1.0m
Kuan=1.5;%板宽1.5m
DG=1.0;%沿高单元网格尺寸
DK=1.5;%沿宽单元网格尺寸
NG=Gao/DG+1;%沿高节点数
(2)
NK=Kuan/DK+1;%沿寛节点数
(2)
t=0.003;%板厚3mm
P=-100;%均布荷载100KN/m(100N/mm)
E=1e9;
v=0.25;
fori=1:
NG%NG=2
forj=1:
NK%NK=2
x((i-1)*NK+j,:
)=[(j-1)*DKGao-(i-1)*DG];%节点坐标
end
end
%节点整体编码与局部编码的对应关系
id(1,:
)=[132];
id(2,:
)=[234];
nelm=(NG-1)*(NK-1)*2;%单元数
(2)
nd1=NG*NK;%节点数(4)
%计算单元的刚度矩阵并组集总刚
K=zeros(2*nd1,2*nd1);%整体刚度矩阵赋零初值
fori=1:
nelm%单元常数
(2)
mA=[1x(id(i,1),1)x(id(i,1),2)
1x(id(i,2),1)x(id(i,2),2)
1x(id(i,3),1)x(id(i,3),2)];%单元节点坐标矩阵
A(i)=0.5*det(mA);%单元面积
mb=[1x(id(i,2),2)%bi
1x(id(i,3),2)];
b(i,1)=-det(mb);
mb=[1x(id(i,3),2)%bj
1x(id(i,1),2)];
b(i,2)=-det(mb);
mb=[1x(id(i,1),2)%bm
1x(id(i,2),2)];
b(i,3)=-det(mb);
mc=[1x(id(i,2),1)%ci
1x(id(i,3),1)];
c(i,1)=det(mc);
mc=[1x(id(i,3),1)%cj
1x(id(i,1),1)];
c(i,2)=det(mc);
mc=[1x(id(i,1),1)%cm
1x(id(i,2),1)];
c(i,3)=det(mc);
%应变矩阵B
mB=[b(i,1),0
0,c(i,1)
c(i,1),b(i,1)];
mB=0.5/A(i)*mB;%Bi
B(:
1:
2,i)=mB;%行(3行),列(1,2),页
mB=[b(i,2),0
0,c(i,2)
c(i,2),b(i,2)];
mB=0.5/A(i)*mB;%Bj
B(:
3:
4,i)=mB;
mB=[b(i,3),0
0,c(i,3)
c(i,3),b(i,3)];
mB=0.5/A(i)*mB;%Bm
B(:
5:
6,i)=mB;
clearmB
%弹性矩阵
mD=[1,v,0
v,1,0
0,0,(1-v)/2];
D(:
:
i)=E/(1-v*v)*mD;%平面应力问题弹性矩阵D
clearmDmcmbGaoKuanDGDKNGNKEEvv
%单元刚度矩阵
Ke(:
:
i)=B(:
:
i)'*D(:
:
i)*B(:
:
i)*t*A(i);%单元刚度矩阵Ke
%单元定位矩阵
L=zeros(6,2*nd1);%6,2*4
L(1:
2,(id(i,1)-1)*2+1:
id(i,1)*2)=eye(2,2);
L(3:
4,(id(i,2)-1)*2+1:
id(i,2)*2)=eye(2,2);%单元定位矩阵
L(5:
6,(id(i,3)-1)*2+1:
id(i,3)*2)=eye(2,2);
%组集整体刚度矩阵
K=K+L'*Ke(:
:
i)*L;
end
%形成整体节点荷载列阵
R=zeros(2*nd1,1);
R(3,1)=P/2;
R(7,1)=P/2;
%加零位移约束,方程缩减法
fori=[8,6,5,2,1]
K(i,:
)=[];
K(:
i)=[];
R(i,:
)=[];
end
uvw=inv(K)*R;
wy=zeros(8,1);
zzz=[3,4,7];
wy(zzz)=uvw;%节点位移矩阵
clearzzz
figure
holdon
hwy=wy*5000;%位移放大1000倍
hx=[x(:
1)+hwy(1:
2:
2*nd1-1),x(:
2)+hwy(2:
2:
2*nd1)]%变形后的节点坐标
fori=1:
nd1
plot(x(i,1),x(i,2),'bx')
plot(x(i,1)+hwy(i*2-1),x(i,2)+hwy(i*2),'ro')
end
fori=1:
nelm
line=[x(id(i,1),:
);
x(id(i,2),:
);
x(id(i,3),:
);
x(id(i,1),:
)];
plot(line(:
1),line(:
2),'b-')
line=[hx(id(i,1),:
);
hx(id(i,2),:
);
hx(id(i,3),:
);
hx(id(i,1),:
)];
plot(line(:
1),line(:
2),'r:
')
end
u=zeros(nd1,1);
w=zeros(nd1,1);
u=hx(1:
2:
2*nd1-1);%水平位移
w=hx(2:
2:
2*nd1);%垂直位移
%建立每一个单元的节点位移
det=zeros(6,nelm);
fori=1:
nelm
det(1,i)=u(id(i,1));
det(2,i)=w(id(i,1));
det(3,i)=u(id(i,2));
det(4,i)=w(id(i,2));
det(5,i)=u(id(i,3));
det(6,i)=w(id(i,3));
end
%计算每一个单元的应变
fori=1:
nelm
epsl(:
i)=B(:
:
i)*det(:
i);
end
%计算每一个单元的应力
fori=1:
nelm
sigm(:
i)=D(:
:
i)*epsl(:
i);
end
(2)
%课程设计:
第一问:
200个单元的三角形单元,这里取一半来分析,然后对称
clc
clear
Gao=1.0;%板高1.0m
Kuan=1.5;%板宽1.5m
DG=0.1;%沿高单元网格尺寸
DK=0.3;%沿宽单元网格尺寸
NG=Gao/DG+1;%沿高节点数11个
NK=Kuan/DK+1;%沿寛节点数6个
t=0.003;%板厚3mm
P=-100;%均布荷载100KN/m(100N/mm)
E=1e9;
v=0.25;
fori=1:
NG%有11个
forj=1:
NK%有6个
x((i-1)*NK+j,:
)=[(j-1)*DKGao-(i-1)*DG];%节点坐标
end
end
%节点整体编码与局部编码的对应关系
n1=1;n2=2;n3=3;n4=4;n5=5;n6=6;n7=7;n8=8;n9=9;n10=10;n11=11;n12=12;
fori=1:
NG-1
id((i-1)*10+1,:
)=[n1n7n2];
id((i-1)*10+2,:
)=[n2n7n8];
id((i-1)*10+3,:
)=[n2n8n3];%节点整体编码与局部编码的对应关系
id((i-1)*10+4,:
)=[n3n8n9];
id((i-1)*10+5,:
)=[n3n9n4];
id((i-1)*10+6,:
)=[n4n9n10];
id((i-1)*10+7,:
)=[n4n10n5];
id((i-1)*10+8,:
)=[n5n10n11];
id((i-1)*10+9,:
)=[n5n11n6];
id((i-1)*10+10,:
)=[n6n11n12];
n1=n1+6;n2=n2+6;n3=n3+6;n4=n4+6;n5=n5+6;n6=n6+6;
n7=n7+6;n8=n8+6;n9=n9+6;n10=n10+6;n11=n11+6;n12=n12+6;
end
clearn1n2n3n4n5n6n7n8n9n10n11n12
nelm=(NG-1)*(NK-1)*2;%单元数100个
nd1=NG*NK;%节点数66个
%计算单元的刚度矩阵并组集总刚
K=zeros(2*nd1,2*nd1);%整体刚度矩阵赋零初值
fori=1:
nelm%单元常数
mA=[1x(id(i,1),1)x(id(i,1),2)
1x(id(i,2),1)x(id(i,2),2)
1x(id(i,3),1)x(id(i,3),2)];%单元节点坐标矩阵
A(i)=0.5*det(mA);%单元面积
mb=[1x(id(i,2),2)
1x(id(i,3),2)];
b(i,1)=-det(mb);
mb=[1x(id(i,3),2)
1x(id(i,1),2)];
b(i,2)=-det(mb);%参数b
mb=[1x(id(i,1),2)
1x(id(i,2),2)];
b(i,3)=-det(mb);
mc=[1x(id(i,2),1)
1x(id(i,3),1)];
c(i,1)=det(mc);
mc=[1x(id(i,3),1)
1x(id(i,1),1)];
c(i,2)=det(mc);%参数c
mc=[1x(id(i,1),1)
1x(id(i,2),1)];
c(i,3)=det(mc);
%应变矩阵B
mB=[b(i,1),0
0,c(i,1)
c(i,1),b(i,1)];
mB=0.5/A(i)*mB;
B(:
1:
2,i)=mB;
mB=[b(i,2),0
0,c(i,2)
c(i,2),b(i,2)];
mB=0.5/A(i)*mB;
B(:
3:
4,i)=mB;
mB=[b(i,3),0
0,c(i,3)
c(i,3),b(i,3)];
mB=0.5/A(i)*mB;
B(:
5:
6,i)=mB;
clearmB
%弹性矩阵
mD=[1,v,0
v,1,0
0,0,(1-v)/2];
D(:
:
i)=E/(1-v*v)*mD;%平面应力问题弹性矩阵D
clearmDmcmbKuanDGDKEEvv
%单元刚度矩阵
Ke(:
:
i)=B(:
:
i)'*D(:
:
i)*B(:
:
i)*t*A(i);%单元刚度矩阵Ke
%单元定位矩阵
L=zeros(6,2*nd1);
L(1:
2,(id(i,1)-1)*2+1:
id(i,1)*2)=eye(2,2);
L(3:
4,(id(i,2)-1)*2+1:
id(i,2)*2)=eye(2,2);%单元定位矩阵
L(5:
6,(id(i,3)-1)*2+1:
id(i,3)*2)=eye(2,2);
%组集整体刚度矩阵
K=K+L'*Ke(:
:
i)*L;
end
%形成整体节点荷载列阵
R=zeros(2*nd1,1);%nd1=66
q1=Gao/NG;
fori=1:
NG
ifi==1|i==NG
R(i*NK*2-1)=0.5*P*q1;
else
R(i*NK*2-1)=P*q1;
end
end
clearq1
%加零位移约束,乘大数法
xyues=1:
NK:
NK*(NG-1)+1;
yyues=[xyues,NK*(NG-1)+2:
nd1];%有约束的节点
%约束的自由度
yszyd=[xyues*2-1,yyues*2];%x方向约束与y方向约束
R(yszyd)=0;
fori=yszyd
K(i,i)=K(i,i)*1e10;
end
deta=inv(K)*R;%节点变形位移
figure
holdon
fori=1:
nelm%单元数100
line=[x(id(i,1),:
);
x(id(i,2),:
);
x(id(i,3),:
);
x(id(i,1),:
)];
plot(line(:
1),line(:
2),'b-')
end%变形前的单元网格
hx=zeros(nd1,2);%所有点位移后的坐标
deta1=deta*1000;%放大1000倍
hx(:
1)=x(:
1)+deta(1:
2:
2*nd1-1);
hx(:
2)=x(:
2)+deta(2:
2:
2*nd1);
figure
holdon
fori=1:
nd1
plot(x(i,1),x(i,2),'bx')
plot(x(i,1)+deta1(i*2-1),x(i,2)+deta1(i*2),'ro')
end
fori=1:
nelm
line=[x(id(i,1),:
);
x(id(i,2),:
);
x(id(i,3),:
);
x(id(i,1),:
)];
plot(line(:
1),line(:
2),'b-')
line=[hx(id(i,1),:
);
hx(id(i,2),:
);
hx(id(i,3),:
);
hx(id(i,1),:
)];
plot(line(:
1),line(:
2),'r:
')
end
u=zeros(nd1,1);
w=zeros(nd1,1);
u=hx(1:
2:
2*nd1-1);%水平位移
w=hx(2:
2:
2*nd1);%垂直位移
%建立每一个单元的节点位移
det=zeros(6,nelm);
fori=1:
nelm
det(1,i)=u(id(i,1));
det(2,i)=w(id(i,1));
det(3,i)=u(id(i,2));
det(4,i)=w(id(i,2));
det(5,i)=u(id(i,3));
det(6,i)=w(id(i,3));
end
%计算每一个单元的应变
fori=1:
nelm
epsl(:
i)=B(:
:
i)*det(:
i);
end
%计算每一个单元的应力
fori=1:
nelm
sigm(:
i)=D(:
:
i)*epsl(:
i);
end
(3)
%课程设计第二问矩形双线性单元,这时以整个图形为分析对象
clc
clear
Gao=2.0;%板高2.0m
Kuan=1.5;%板宽1.5m
he=2;%横向节点数
sh=3;%竖向节点数
nelm=(he-1)*(sh-1);%单元数
nd1=he*sh;%节点数
bb=Gao/(sh-1)/2;%局部坐标与正则坐标转换的系数b
aa=Kuan/(he-1)/2;%局部坐标与正则坐标转换的系数a
t=0.003;%板厚3mm
P=-100;%均布荷载100KN/m(100N/mm)
E=1e9;
v=0.25;
x=[0,2;
1.5,2;
0,1;
1.5,1;
0,0;
1.5,0];%变形前的节点直角坐标
id=[3,4,2,1;
5,6,4,3];
symsssnn%先定义2个符号变量,由课本p137的形函数公式来计算形函数
N=1/4*[(ss-1)*(nn-1);
-(ss+1)*(nn-1);
(ss+1)*(nn+1);
-(ss-1)*(nn+1)];
ndx=1/aa*diff(N,ss);
ndy=1/bb*diff(N,nn);
B(3,2,4)=sym('0');
fori=1:
4%每个单元有4个节点,3行2列4页
B(:
:
i)=[ndx(i),0;
0,ndy(i);
ndy(i),ndx(i)];
end
D=E/(1-v*v)*[1,v,0
v,1,0
0,0,(1-v)/2];
K=zeros(12,12);
forie=1:
nelm%这里单元数为2
fori=1:
4
forj=1:
4
K(id(ie,i)*2-1:
id(ie,i)*2,id(ie,j)*2-1:
id(ie,j)*2)=...
K(id(ie,i)*2-1:
id(ie,i)*2,id(ie,j)*2-1:
id(ie,j)*2)+...
aa*bb*t*double(int(int(B(:
:
i)'*D*B(:
:
j),ss,-1,1),nn,-1,1));
end
end
end
R=zeros(12,1);
R(3)=-50;
R(7)=-100;
R(11)=R(3);
yszyd=[1,2,5,6,9,10];
fori=yszyd
R(i)=0;
K(i,i)=K(i,i)*1e10;
end
deta=inv(K)*R;
deta1=deta*1000;%位移放大1000倍
hx=[x(:
1)+deta1(1:
2:
2*nd1-1),x(:
2)+deta1(2:
2:
2*nd1)]%变形后的节点坐标
figure
holdon
fori=1:
nelm
line=[x(id(i,1),:
);
x(id(i,2),:
);
x(id(i,3),:
);
x(id(i,4),:
);
x(id(i,1),:
)];
plot(line(:
1),line(:
2),'b-')
line=[hx(id(i,1),:
);
hx(id(i,2),:
);
hx(id(i,3),:
);
hx(id(i,4),:
);
hx(id(i,1),:
)];
plot(line(:
1),line(:
2),'r*')
end
u=zeros(nd1,1);
w=zeros(nd1,1);
u=hx(1:
2:
2*nd1-1);%水平位移
w=hx(2:
2:
2*nd1);%垂直位移
%建立每一个单元的节点位移
det=zeros(8,nelm);
fori=1:
nelm
det(1,i)=u(id(i,1));
det(2,i)=w(id(i,1));
det(3,i)=u(id(i,2));
det(4,i)=w(id(i,2));
det(5,i)=u(id(i,3));
det(6,i)=w(id(i,3));
det(7,i)=w(id(i,4));
det(8,i)=w(id(i,4));
end
%计算每一个单元的应变
fori=1:
nelm
epsl(:
i)=B(:
:
i)*det(:
i);
end
%计算每一个单元的应力
fori=1:
nelm
sigm(:
i)=D(:
:
i)*epsl(:
i);
end
2013.6.12
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 有限元 课程设计