有限元作业三角形单元求解.docx
- 文档编号:7244797
- 上传时间:2023-01-22
- 格式:DOCX
- 页数:17
- 大小:41.64KB
有限元作业三角形单元求解.docx
《有限元作业三角形单元求解.docx》由会员分享,可在线阅读,更多相关《有限元作业三角形单元求解.docx(17页珍藏版)》请在冰豆网上搜索。
有限元作业三角形单元求解
有限元作业»
年级2015级
学院机电工程学院专业名称
班级学号
学生姓名
2016年05月
如下图所示为一受集中力P作用的结构,弹性模量E为常量,泊松比V=1/6,厚度为1=1。
按平面应力问题计算,运用有限元方法,分别采用三角形及四边
形单元求解,求节点位移及单元应力(要求三角形单元数量不少于4个,四边形单元不少于2个)
图
(一)
图
(二)三角形单元求解
图(三)四边形单元求解
(1)如图划分三角形单元,工分成四个分别为④
(2)如图分别进行编号1、2、3、4、5、6,并建立坐标系
(3)编程进行求解,得出结果,其中假设力P=2000N
调用Triangle2D3Node_Stiffness函数,求出单元刚度矩阵k1=
1.0e+06
7.2857
-3.0000
-2.1429
0.8571
-5.1429
2.1429
-3.0000
7.2857
2.1429
-5.1429
0.8571
-2.1429
-2.1429
2.1429
2.1429
0
0
-2.1429
0.8571
-5.1429
0
5.1429
-0.8571
0
-5.1429
0.8571
0
-0.8571
5.1429
0
2.1429
-2.1429
-2.1429
0
0
2.1429
k2=
1.0e+06
5.1429
0
-5.1429
0.8571
0
-0.8571
0
2.1429
2.1429
-2.1429
-2.1429
0
5.1429
2.1429
7.2857
-3.0000
-2.1429
0.8571
0.8571
-2.1429
-3.0000
7.2857
2.1429
-5.1429
0
-2.1429
-2.1429
2.1429
2.1429
0
0.8571
0
0.8571
-5.1429
0
5.1429
k3=
1.0e+06
2.1429
0
-2.1429
-2.1429
0
2.1429
0
5.1429
-0.8571
-5.1429
0.8571
0
-2.1429
-0.8571
7.2857
3.0000
-5.1429
-2.1429
-2.1429
-5.1429
3.0000
7.2857
-0.8571
-2.1429
0
0.8571
-5.1429
-0.8571
5.1429
0
2.1429
0
-2.1429
-2.1429
0
2.1429
k4=
1.0e+06
2.1429
0
-2.1429
-2.1429
0
2.1429
0
5.1429
-0.8571
-5.1429
0.8571
0
-2.1429
-0.8571
7.2857
3.0000
-5.1429
-2.1429
-2.1429
-5.1429
3.0000
7.2857
-0.8571
-2.1429
0
0.8571
-5.1429
-0.8571
5.1429
0
2.1429
0
-2.1429
-2.1429
0
2.1429
函数,求出总体刚度矩阵
调用Triangle2D3Node_Assembly
0.7286
-0,1000
-0.2K3
9.0857
-0.5143
CL2143
0
0
0
0
0
0
-0.3000
0.72即
0.2M3
00S57
-0.214.3
0
Q
g
0
0
0
-0,2143
0-2U3
0-72S6
0
Q
-0,3G00
-0=3143
D-085?
Q
0
0
0
CL0057
-0-5143
D
-a.3ooo
0
0.2143
Th2hB
a
0
0
Q
-0/5143
0adB57
0
-0.3QOO
1.4571
0
-0.42册
0
0
a.3ooo
-0-5113
-0.0357
0,2143
-0.2143
-0.3WD
0
0
1,4571
0
-1.02W
(L300(1
0
-0-2143
-0,2143
0
0
-0,5343
-a伽
0
I.4E71
4
-«S5I43
-0.5143
0
0
0
0
山0857
-Ou2143
a
-|a02BS
0
-ft.OSS"
-0.2143
0
0
0
0
0
0
0
CL3000
-0-5143
-0..0857
0.728S
0
-0-2113
-0.21^3
0
0
0
0
仇3000
0
-0.2H3
0
0.7266
-0.0S57
-0r5H3
0
0
0
0
75143
-0,3143
Q
0
-Ou2143
-0,065:
Q.12B8
乳3000
0
0
D
0
■仇D857
-CL2U3
0
0
-0.2143
-0.SU3
0.3M0
D.72A6
求出的节点位移
0
0
0
-0.0004
0.0008
0.0005
0.0010
0.0007
0.0023
-0.0007
0.0026
中求出的分别为
调用Triangle2D3Node_Stress函数,求出应力,S1、S2、S3、
Sx,Sy,Sxy
S1=
1.0e+03*
-4.4086
-0.7348
3.5914
S2=
1.0e+03
4.4086
-0.6405
0.4086
53=
1.0e+03*
1.8907
-1.0601
2.1093
54=
1.0e+03*
-1.8907
2.1093
1.8907
二二、
(1)如图划分四边形单元,工分成四个分别为
(2)如图分别进行编号1、2、3、4、5、6,并建立坐标系
(3)编程进行求解,得出结果,其中假设力P=2000N
调用Quad2D4Node_Stiffness函数,求出单元刚度矩阵
4.8S71
-1.5H0
0.2857
-0.6429
-2.4286
1.5000
-2.7143
0.6429
-1.5000
18571
0.6429
-2.T143
L5QOO
-2.4286
-0.642d
山2867
0=2857
0.M29
4.8571
L50QQ
-2.7143
-CL6429
-Z4286
-L50W
-0.6429
-2.7143
1.5000
4.8571
0.6129
0.2857
-L5000
-2.4286
-2.4286
1.6000
-2.7143
0.石竝9
4.3571
-1-5000
CL2857
-0.6428
L5000
-2,42B6
-0.C429
0,2B57
^IrSOOO
4・B57i
0,9429
一Z7143
-2.7143
-0.6t29
4286
-L5000
0.28B7
0.6429
4.8671
1.5000
Q.6429
Du2857
-L5000
-N42鯛
-0.6i2S
-2-7U3
L5Q00
4.S57J
k2二
LOe+06*■
4・8S71
1.5Q0C
-2.7143
-0.6429
-2.4286
-1.5000
0,28S7
Ou€429
1B5000
札3571
0.6120
0.2B57
-LSOOO
-2.4236
-Du6429
-2/7143
-2.7143
0.0429
4.S571
-i.noo
0.2857
-0.6429
-2,4286
LSOOO
-0.6429
0.2867
-L.EOOO
4.8671
0.6429
*2.7143
1.6000
-2.4286
-2.4236
-L5000
0.235:
4.8571
1-5000
-2-7143
-0.6^25
-1.5000
-2.4286
-0.6429
-2.7143
1,5000
4・S57I
0.S429
0l2857
0.2S5T
-0.642B
-2.1286
L500S
-2.7143
0.6429
4・8571
-LSOOO
Q.£429
-2.7U3
hSQOQ
-N42S6
-Q+6429
0-2857
-扎5Q0Q
4.8571
1.Ofl+Ofl*
4.8571
-1=5000
CL2S57
■山@439
-2...42BS
K5000
-2.71113
0.0429
0
0
0
0
-1.5400
4.8571
CL6429
-2,7143
LGflOO
■242B6
-0.6129
0,2657
0
0
0
0
D.2S6T
0.6i2»
4・8B71
]■.5000
-2.7143
-0.6429
-2.42M
0
0
D
Q
-0.6129
-2.7143
1.5000
4.8S71
(L6429
0.2557
-1.5000
-2.428fi
0
0
0
0
-2.^283
£.5000
-瓷7143
0.6420
9-.7143
0
O・571也
-2.7H3
-CL開2B
-£42&6
-1.50DO
K500Q
-2.4280
-CLB42-9
CL2857
0
9.7U3
0
-5.428A
0.6429
CL2857
-L5000
-2u42BG
-2.7143
7*429
-2.42SC
-3,5000
OU5714
0
9.7143
0
■監428C
L50G0
-2.7143
56429
D■仙
0.2B67
-I.S004
-2.4S86
0
-6.42B6
0
9.7143
1.5000
-2a42&6
7、6429
也2857
0
0
0
-2.7143
0.6429
-2.428fl
h&DQO
4,9571
-1.EOOfl
0,2357
-a.6429
0
0
0
-0US42B
0.2S5"
1.5000
-2.
-I,5000
LB571
OL6429
-2-7143
a
0
fl
0
-i42B6
-K5G00
-2.7U3
-D.A42»
0.2B57
CLB42B
4.8571
K50D0
Q
0
0
0
-LEOOO
临428E
0.C42B
0.39E7
79429
临7143
L5000
4.8571
求出节点位移
0
0
0
0.0012
0.0017-0.0012
0.0017
0.0016
0.0049
-0.0017
0.0052
S3分别为
调用Quad2D4Node_Stress函数,求出单元应力中的的S1、S2、
Sx,Sy,Sxy应力分量
S1=
1.0e+03*
0.0000-0.2478
2.0000
S2=
1.0e+07
0.68564.1135
-1.7137
程序附录
、
1、三角形单元总程序:
E=1e7;
NU=1/6;
t=1;
ID=1;
%调用Triangle2D3Node_Stiffness函数,求出单元刚度矩阵
k仁Triangle2D3Node_Stiffness(E,NU,t,0,1,0,0,1,1,ID)
k2=Triangle2D3Node_Stiffness(E,NU,t,0,0,1,0,1,1,ID)
k3=Triangle2D3Node_Stiffness(E,NU,t,1,1,1,0,2,0,ID)
k4=Triangle2D3Node_Stiffness(E,NU,t,2,0,2,1,1,1,ID)
%调用Triangle2D3Node_Assembly函数,求出总体刚度矩阵
KK=zeros(12,12);
KK=Triangle2D3Node_Assembly(KK,k1,1,2,3);
KK=Triangle2D3Node_Assembly(KK,k2,2,4,3);
KK=Triangle2D3Node_Assembly(KK,k3,3,4,5);
KK=Triangle2D3Node_Assembly(KK,k4,5,6,3)
%边界条件的处理及刚度方程求解
k=KK(5:
12,5:
12)
p=[0;0;0;0;0;0;0;2000]
u=k\p
%支反力的计算
U=[0;0;0;0;u]%为节点位移
P=KK*U
%调用Triangle2D3Node_Strain函数,求出应变SN1、SN2、SN3中求出的分
别为SNx,SNy,SNxy
u1=[U
(1);U
(2);U(3);U⑷;U(5);U(6)];
u2=[U(3);U⑷;出7);U(8);U(5);U(6)];
u3=[U(5);U(6);U(7);U(8);U(9);U(10)];
u4=[U(9);U(10);U(11);U(12);U(5);U(6)];
SN1=Triangle2D3Node_Strain(0,1,0,0,1,1,u1)
SN2=Triangle2D3Node_Strain(0,0,1,0,1,1,u2)
SN3=Triangle2D3Node_Strain(1,1,1,0,2,0,u3)
SN4=Triangle2D3Node_Strain(2,0,2,1,1,1,u4)
%调用Triangle2D3Node_Stress函数,求出应力,S1、S2、S3、中求出的分别
为Sx,Sy,Sxy
u仁[U
(1);U
(2);U(3);U⑷;U(5);U(6)];
u2=[U(3);U⑷;出7);U(8);U(5);U(6)];
u3=[U(5);U(6);U(7);U(8);U(9);U(10)];
u4=[U(9);U(10);U(11);U(12);U(5);U(6)];
S1=Triangle2D3Node_Stress(E,NU,0,1,0,0,1,1,u1,ID)
S2=Triangle2D3Node_Stress(E,NU,0,0,1,0,1,1,u2,ID)
S3=Triangle2D3Node_Stress(E,NU,1,1,1,0,2,0,u3,ID)
S4=Triangle2D3Node_Stress(E,NU,2,0,2,1,1,1,u4,ID)
2、求刚度矩阵程序
functionk=Triangle2D3Node_Stiffness(E,NU,t,xi,yi,xj,yj,xm,ym,ID)
%该函数计算单元的刚度矩阵
%输入弹性模量E,泊松比NU,厚度t
%输入三个节点i、j、m的坐标xi,yi,xj,yj,xm,ym
%输入平面问题性质指示参数ID(1为平面应力,2为平面应变)
%输出单元刚度矩阵k(6X6)
%
A=(xi*(yj-ym)+xj*(ym-yi)+xm*(yi-yj))/2;
betai=yj-ym;
betaj=ym-yi;
betam=yi-yj;
gammai=xm-xj;
gammaj=xi-xm;
gammam=xj-xi;
B=[betai0betaj0betam0;
0gammai0gammaj0gammam;
gammaibetaigammajbetajgammambetam]/(2*A);
ifID==1
D=(E/(1-NU*NU))*[1NU0;NU10;00(1-NU)/2];
elseifID==2
D=(E/(1+NU)/(1-2*NU))*[1-NUNU0;NU1-NU0;00(1-2*NU)/2];end
k=t*A*B'*D*B;
3、求整体刚度矩阵
functionz=Triangle2D3Node_Assembly(KK,k,i,j,m)
%该函数进行单元刚度矩阵的组装
%输入单元刚度矩阵k
%输入单元的节点编号I、j、m
%输出整体刚度矩阵KK
%
DOF
(1)=2*i-1;
DOF
(2)=2*i;
DOF(3)=2*j-1;
DOF⑷=2*j;
DOF(5)=2*m-1;
DOF(6)=2*m;
forn仁1:
6
forn2=1:
6
KK(DOF(n1),DOF(n2))=KK(DOF(n1),DOF(n2))+k(n1,n2);
end
end
z=KK;
4、求应变程序
functionstrain=Triangle2D3Node_Strain(xi,yi,xj,yj,xm,ym,u)
%该函数计算单元的应变
%输入三个节点i、j、m的坐标xi,yi,xj,yj,xm,ym
%输入单元的位移列阵u(6X1)
%输出单元的应力strain(3X1),由于它为常应变单元,则单元的应变分量为
SNx,SNy,SNz
%
A=(xi*(yj-ym)+xj*(ym-yi)+xm*(yi-yj))/2;
betai=yj-ym;
betaj=ym-yi;
betam=yi-yj;
gammai=xm-xj;
gammaj=xi-xm;
gammam=xj-xi;
B=[betai0betaj0betam0;
0gammai0gammaj0gammam;
gammaibetaigammajbetajgammambetam]/(2*A);
strain=B*u;
5、求应力程序
functionstress=Triangle2D3Node_Stress(E,NU,xi,yi,xj,yj,xm,ym,u,ID)
%该函数计算单元的应力
%输入弹性模量E,泊松比NU,厚度t
%输入三个节点i、j、m的坐标xi,yi,xj,yj,xm,ym
%输入平面问题性质指示参数ID(1为平面应力,2为平面应变),单元的位移列
阵u(6X1)
%输出单元的应力stress(3X1),由于它为常应力单元,则单元的应力分量为
Sx,Sy,Sxy
%
A=(xi*(yj-ym)+xj*(ym-yi)+xm*(yi-yj))/2;
betai=yj-ym;
betaj=ym-yi;
betam=yi-yj;
gammai=xm-xj;
gammaj=xi-xm;
gammam=xj-xi;
B=[betai0betaj0betam0;
0gammai0gammaj0gammam;
gammaibetaigammajbetajgammambetam]/(2*A);
ifID==1
D=(E/(1-NU*NU))*[1NU0;NU10;00(1-NU)/2];
elseifID==2
D=(E/(1+NU)/(1-2*NU))*[1-NUNU0;NU1-NU0;00(1-2*NU)/2];endstress=D*B*u;
1、四边形单元总程序:
E=1e7;
NU=1/6;
h=1;
ID=1;
%调用Quad2D4Node_Stiffness函数,求出单元刚度矩阵
k1=Quad2D4Node_Stiffness(E,NU,h,0,1,0,0,1,0,1,1,ID)
k2=Quad2D4Node_Stiffness(E,NU,h,1,0,2,0,2,1,1,1,ID)
%调用Quad2D4Node_Assembly函数,求出求出总体刚度矩阵
KK=zeros(12,12);
KK=Quad2D4Node_Assembly(KK,k1,1,2,3,4);
KK=Quad2D4Node_Assembly(KK,k2,3,5,6,4)
%边界条件的处理及刚度方程求解
k=KK(5:
12,5:
12)
p=[0;0;0;0;0;0;0;2000]
u=k\p
%支反力的计算
U=[0;0;0;0;u]%为节点位移
P=KK*U
%调用Quad2D4Node_Stress函数,求出单元应力中的的S1、S2、S3分别为
Sx,Sy,Sxy应力分量
u1=[U
(1);U
(2);U(3);U⑷;U(5);U(6);U(7);U(8)];
u2=[U(5);U(6);U(9);U(10);U(11);U(12);U(7);(8)];
S1=Quad2D4Node_Stress(E,NU,0,1,0,0,1,0,1,1,u1,ID)
S2=Quad2D4Node_Stress(E,NU,1,0,2,0,2,1,1,1,u2,ID)
2、求刚度矩阵程序
functionk=Quad2D4Node_Stiffness(E,NU,h,xi,yi,xj,yj,xm,ym,xp,yp,ID)
%该函数计算单元的刚度矩阵
%输入弹性模量E,泊松比NU,厚度h
%输入4个节点i、j、m、p的坐标xi,yi,xj,yj,xm,ym,xp,yp
%输入平面问题性质指示参数ID(1为平面应力,2为平面应变)
%输出单元刚度矩阵k(8X8)
%
symsst;
a=(yi*(s-1)+yj*(-1-s)+ym*(1+s)+yp*(1-s))/4;
b=(yi*(t-1)+yj*(1-t)+ym*(1+t)+yp*(-1-t))/4;
c=(xi*(t-1)+xj*(1-t)+xm*(1+t)+xp*(-1-t))/4;
d=(xi*(s-1)+xj*(-1-s)+xm*(1+s)+xp*
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 有限元 作业 三角形 单元 求解