有限元分析课程作业.docx
《有限元分析课程作业.docx》由会员分享,可在线阅读,更多相关《有限元分析课程作业.docx(33页珍藏版)》请在冰豆网上搜索。
有限元分析课程作业
有限元分析》课程作业
任课教师:
徐亚兰
学生姓名:
陈新杰
学号:
班级:
1304012时间:
2016-01-05
、问题描述及分析
问题:
如图1所示,有一矩形平板,在右侧受到P=10KN/m的分布力,材料常数为:
弹性模量E1107pa;泊松比1/3;
板的厚度为t=;试按平面应力问题利用三角形与矩形单元分别计算各个节点位移及支座反力。
图1平面矩形结构的有限元分析
分析:
使用两种方案:
一、基于3节点三角形单元的有
限元建模,将矩形划分为两个3节点三角形单元;二、基于4节点矩形单元的有限元建模,使用一个4节点矩形单元。
利用MATLAB^件计算出各要求量,再将两种方案的计算结果进行比较、分析、得出结论。
二、有限元建模及分析
1、基于3节点三角形单元的有限元建模及分析
(1)结构的离散化与编号
元。
单元①三个节点的编号为1,2,4,单元②三个节点的
编号为3,4,2,各个节点的位置坐标为xi,yi,i1,2,3,4,各个
节点的位移(分别沿x方向和y方向)为Ui,Vi,i1,2,3,4。
图2方案一:
使用两个3节点三角形单元
(2)各单元的刚度矩阵及刚度方程
a.单元的几何和节点描述
单元①有6个节点位移自由度(DOF。
将所有节点上的位移组成一个列阵,记作q⑴;同样,将所有节点上的各个力也组成一个列阵,记作F⑴,则有
q⑴(Ui,w,U2,V2,U4,V4)
同理,对于单元②,有
q
(2)(U3,V3,U4,V4,U2,V2)
b.单元的位移场描述
对于单元①,设位移函数
u(x,y)aa/a?
y
v(x,y)bobxb2y
由节点条件,在xx,yyi处,有
(1-1)
u(Xi,yi)uii1,2,4
v(Xj,yJVi
(1-2)
将式(1-1)代入节点条件式(1-2)中,可求出式(
1-1)
中待定系数,即
a。
a’
a2
bo
b1
b2
U1X1
1
—u2X2
2A
U4X4
1
2A
U1
U2
U4
1
1
—1
2A
1
1
2A
X1
X2
X4
GW
y1
y2
y4
y1
y2
y4
U1
U2
U4
a2v2
1
2A
(ay
a?
U2a3U4)
2A
(dm
b2U2b3U4)
a3V4)(1-6)
药冲1b2V2b3V4)(1-7)
1
2a(C1V1C2V2C3V4)(1-8)
2A
(1-3)
(1-4)
(1-5)
在式(1-3)~式(1-8)中
1X1y1
1X2y2
1X4y
ai
x4y4
上式中的符号(1,2,3)表示下标轮换,如12,23,31同时更
换。
将单元①各节点的位置坐标&0,y10必1,y20,&0,y4
代入得
将系数式(1-3)〜式(1-8)式代入(1-1)中,重写位移函
数,并以节点位移的形式进行表示,有
u(x,y)N^x,y)uN2(x,y)U2v(x,y)N1(x,y)vN2(x,yM
(1-12)
位移函数式(1-11)写成矩阵形式,有
u
(1)(x,y)u!
x,y)
(26)v(x,y)
V1
U2
v
U4
V4
(1-13)
对于单元②,
过程同上,
形状函数矩阵
⑵/、
H(x,y)
(1-14)
位移函数
(2)
u(x,y)
(26)
u(x,y)
v(x,y)
1+x
0
U3
V3
U4
V4
U2
V
(1-15)
c.单元的应变场描述
对于单元①,
应变函数
xx
⑴(x,y)
(31)
yy
xy
u
x
v
y
uv
yx
u(x,y)v(x,y)
Ui
Vi
x0y0u2
V2
U4
V4
-101000
=0-10001
-1-10110
V4
其中几何矩阵
(1)/、
B(x,y)
(3"6)
-101000
0-10001
-1-10110
(1-17)
对于单元②,应变方程
—0
x
1xy01x01y0
0
y01xy01x01y
1
0
-1
0
0
0
V3
0
1
0
0
0
-1
U4
1
1
0
-1
-1
0
V4
U2
V
U3
(1-18)
U2
其中几何矩阵
B(x,y)0
(36)
1
0-1000
1000-1
10-1-10
(1-19)
d.单元的应力场描述
xx
(x,y,Z)yy
(31)
xy
0
0
1
"~2
xx
yy
xy
=D
(33)
(1)
(3
(1)
1)
(1-20)
其中,弹性系数矩阵
(1)
D)
0
0
1
2
107
~~2
1-1
=3.75
106
(1-21
(1)
(31)
(1)
(1)
(1)
(1)
(1-22)
其中应力函数矩阵
s
(1)为
(1)
(1)
(1)
D)(B)=375
3
1
0
-1
0
1
0
0
0
-3
-1
3
0
0
1
1061
3
0
0
-1
0
0
0
1=3.75
106-1
-3
1
0
0
3
0
0
1
-1
-1
0
1
1
0
-1
-1
0
1
1
0
(1-23)
应力方程⑴为
-3
-1
3
0
0
1
-3
-1
3
0
0
1
(1)=3.75106-1
-3
1
0
0
3
(1)
q=3.75
106-1
-3
1
0
0
3
(31)
-1
-1
0
1
1
0
(61)
-1
-1
0
1
1
0
u2v2
v1
Pa
u4
v4
(1-24)
对于单元②,过程同上
弹性系数矩阵D
(2)为
(2)
D=3.7510
(33)
31
0
613
0
00
1
应力函数矩阵
s
(2)为
31
-30
0
-1
(2)
S
(2)=13
-10
0
-3
(36)11
0-1
-1
0
(1-25)
(1-26)
应力方程
(2)为
u3
3
1
-3
0
0
-1
3
1
-3
0
0
-1
(2)6
=3.751061
3
-1
0
0
-3
(2)
q=3.75
1061
3
-1
0
0
-3
(31)
1
1
0
-1
-1
0
(61)
1
1
0
-1
-1
0
v3
u4
4Pav4u2
(1-27)
e.单元的势能表达
K
(1)是单元刚度矩阵,即
K
(1)
(1)B
(1)TD
(1)B
(1)d
(1)B
(1)TD
(1)B
(1)dAt(1-28)
其中薄板厚度t0.1m。
将式(1-17)、式(1-21)代入式(1-29),
得到单元①的刚度阵
K
(1)B
(1)TD
(1)B
(1)tA3.75
106
0.1
-1
0
1
0
0
0
0
-1
0
0
0
1
-1
-1
0
1
1
0
-1
0
-1
0
-1
-1
计算得
同理,
u1
v1
u2
v2
u4
v4
7.5
3.75
5.625
1.875
1.875
1.875
u1
3.75
7.5
1.875
1.875
1.875
5.625
v1
(1)55.625
K
(1)105
-1.875
5.625
0
0
1.875
u2
-1.875
1.875
0
1.875
1.875
0
v2
1.875
1.875
0
1.875
1.875
0
u4
1.875
5.625
1.875
0
0
5.625
v4
得到单元②的刚度阵为
u3
v3
u4
v4
u2
v2
7.5
3.75
5.625
1.875
1.875
1.875
u1
3.75
7.5
1.875
1.875
1.875
5.625
v1
(2)55.625
(2)105
-1.875
5.625
0
0
1.875
u2
-1.875
1.875
0
1.875
1.875
0
v2
1.875
1.875
0
1.875
1.875
0
u4
1.875
5.625
1.875
0
0
5.625
v4
K
将两个单元按节点位移所对应的位置进行组装,得到
总体刚度矩阵为
KKK
(1)K
(2)
节点力
F
FPx3FPx3
FRx4
10.110103(FRx101010Frx40)
0.5103(Frx101010FRx40)N
系统的势能
UW=1qTKq-FTq(计算结果在下面呈现)
(4)边界条件的处理及方程求解
边界条件为UiViU4V40。
因此,将针对节点2和节点3的位移求解,节点2和节点3对应总体刚度阵KK中的第3行到第6行、第3列到第6列,贝U需从KK88中提出,置给k,然后生成对应的载荷列阵p,再采用高斯消去法进行求解。
>>k=KK(3:
6,3:
6);
>>p=[500;0;500;0];
>>u=k\p
u=[将列排成了行]
再计算支反力。
在得到整个结构的节点位移后,由原整体刚度方程就可以计算出对应的支反力;先将上面得到的位移结果与边界条件的节点位移进行组合,得到整体的位移列阵U(81),再代回原整体刚度方程,计算出所有的节点力,按照位置关系找出对应的支反力。
>>U=[0;0;u;0;0]
U=0000
[将列排成了行]
>>P=KK*U
P=-50050005000-500
[将列排成了行]所以,节点1的支反力为FRx1500N,FRy1-176.4706N,节点2的支反力为FRx2500N,FRy2176.4706N。
根据已求得的位移和支反力计算系统的势能。
>>A=*U'*KK*U-P'*U
A=
(5)结果分析上述支反力计算结果满足静力平衡,验证了以上求解过程及MATLAB算法的正确性。
2、基于四节点四边形单元的有限元建模及分析
(1)结构的离散化与编号
如图3所示一个4节点矩形单元,单元的节点位移共有8个自由度(DOF。
节点编号为123,4,各自的位置坐标为x*,i123,4,各个节点的位移(分别沿x方向和y方向)
为ui,vi,i1,2,3,4。
4o
图3方案二:
使用一个4节点矩形单元
(2)局部坐标系下单元的描述a.单元的几何和节点描述
采用无量纲坐标
=Xy
a'b
其中a0.5,b0.5。
则单元四个节点的几何位置为
i,i1
1,21
1,31
1,41
将所有节点上的位移组成一个列阵,记作q;同样,将
所有节点上的各个力也组成一个列阵,记作F,则有
q
(81)
(F)
(U1V1U2V2U3V3U4V4)T
(FX1Fy1Fx2Fy2Fx3Fy3Fx4Fy4)
b.单元的位移场描述
设位移函数为
u(x,y)aoaxa?
yasxy
v(x,y)bo
b|Xb2yb3xy
由节点条件,
在xx’yyi处,有
u(Xi,yJUi
i123,4
v(Xj,yJVi
将位移试函数代入节点条件中,求出待定系数
和bo,bl,b2,b3,再代入位移函数中,整理后得
u(x,y)
Ndx,y)u1
n
2(x,y)U2
N3(x,y)U3
N4(x,y)U4
v(x,y)
N,x,y)v1
n
2(x,y)v2
N3(x,y)v3
N4(x,y)v4
其中
M(x,y)
1
—12x
4
1
2y
N2(x,y)
112x
4
1
2y
Wx,y)
112x
4
1
2y
N4(x,y)
112x
4
1
2y
如以无量纲坐标来表达,可写成
1
Ni-1i1i,i123,4
4
将11,1121,21,31,3141,4
a0,ai,a2,a3
1带入上式
得到形状函数矩阵
1
N1
—
1
1
4
N2
1
1-
1
4
N3
1
1
1
4
N4
1
1-
1-
4
写成矩阵形式,有
Ui
Vi
U2
N10N20N30N40v2
0N10N20N30N4u3Nq
C.单元的应变场描述
单元应变为
其中几何矩阵B(x,y)为
应力表达式为
(3i)D)(3i)dm羸儿其中,应力函数矩阵SDB。
e.单元的势能表达
以上已将单元的三大基本变量U,,用基于节点的位移
其中K为4节点矩形单元的刚度矩阵,即
其中,t为薄板的厚度,t0.1m,上式的各个字块矩阵为
检0.1bTdBdr,s1,2,3,4
f.单元刚度阵及刚度方程
单元刚度阵在上面已经列出。
将单元的势能对节点位移
q取一阶极值,可得到单元的刚度方程
(4)边界条件的处理及方程求解
处理方法与3节点三角形单元一致,利用上述求解程序具有的可移植性,简化了求解过程。
>>k=K(3:
6,3:
6);
>>p=[500;0;500;0];
>>u=k\p
u=*
[将列排成了行]
再计算支反力。
同样注意按照位置关系找出对应的支反
力。
>>U=[0;0;u;0;0]
0000
[将列排成了行]
>>P=K*U
P=
-50050005000-500
[将列排成了行]所以,节点1的支反力为FRx1500N,FRy1-111.1111N,节点
2的支反力为FRx2500N,FRy2111.1111N。
根据已求得的位移和支反力计算系统的势能。
>>A=*U'*K*U-P'*U
A=
(5)结果分析
基于4节点矩形单元计算出的势能小于基于3节点三角
形单元计算出的结果,若将该系统分为更多的单元,计算精度也将提高。
3.两种方案的比较与分析从以上计算可以看出,用三角形单元计算时,由于形函数是完全一次式,因而其应变场在单元内均为常数;而对于四边形单元,其形函数带有二次式,计算得到的应变场和应力场都是坐标的一次函数,但不是完全的一次函数,对提高精度有一定作用;根据最小势能原理,势能越小,则整体计算精度越高,比较两种单元计算得到的系统势能,可看出,在相同的节点自由度情况下,矩形单元的计算精度要比三角形单元高。
三、基于MATLAB勺编程实现
1.基于3节点三角形单元的有限元编程实现
(1)程序编写说明
Triangle2D3Node_Stiffness(E,NU,t,xi,yi,xj,yj,xm,ym,ID)
该函数计算单元的刚度矩阵,输入模量E,泊松比NU,
厚度t,三个节点i,j,m平面问题性质指示参数(1为平面应力,2为平面应变),输出单元刚度矩阵k(66)。
Triangle2D3Node_Assemble(KK,k,i,j,m)
该函数进行单元刚度矩阵的组装,输入单元刚度矩阵k,单元的节点编号i,j,m输出整体刚度矩阵KKo
Triangle2D3Node_Stress(E,NU,xi,yi,xj,yj,xm,ym,u
ID)
该函数计算单元的应力,输入弹性模量E,泊松比NU,
厚度t,三个节点i,j,m平面问题性质指示参数(1为平面应力,2为平面应变),单元的位移列阵u(61),输出单元的应力,由于它为常应力单元,则单元的应力分量为Sx,Sy,Sxy。
2)程序清单
%%%%%Triangle2D3Node%%%begin%%%%%%%function
k=Triangle2D3Node_Stiffness(E,NU,t,xi,yi,xj,yj,xm,y
m,ID)
%该函数计算单元的刚度矩阵
%俞入弹性模量E、泊松比NU和厚度t
%俞入3个节点i,j,m的坐标xi,yi,xj,yj,xm,ym
%俞入平面问题性质指示参数ID(1位平面应力,2为平
面应变)
%俞入单元刚度矩阵k(6*6)
%
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;0
0(1-2*NU)/2];
end
k=t*A*B'*D*B;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%functionz=Triangle2D3Node_Assemble(KK,k,i,j,m)%该函数进行单元刚度矩阵的组装
%输入单元刚度矩阵k
%输入单元的节点编号i,j,m
%输入整体刚度矩阵KKrf
%
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;
forn1=1:
6
forn2=1:
6
KK(DOF(n1),DOF(n2))=KK(DOF(n1),DOF(n2))+k(n1,n2);
end
end
z=KK;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%function
stress=Triangle2D3Node_Stress(E,NU,xi,yi,xj,yj,xm,y
m,u,ID)
%该函数计算单元的应力
%俞入弹性模量E、泊松比NU和厚度t
%输入平面问题性质指示参数ID(1位平面应力,2为平
面应变),单元的位移列阵u(6*1)
%俞出单元的应力stress(3*1),由于它为常应力单元,则单元的应力分量Sx,Sy,Sxy
%
A=(xi*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=[batai0bataj0betam0;
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;0
0(1-2*NU)/2];
endstress=D*B*u;%%%%%%%%%%%%%%%%%%%%%%%%%%
>>E=1E7;
>>NU=1/3;
>>t=;
>>c=Triangle2D3Node_Stiffness(E,NU,t,0,0,1,0,0,1,2)
>>CC=zeros(8,8);
>>CC=Triangle2D3Node_Assemble(KK,k1,1,2,4);
>>CC=Triangle2D3Node_Assemble(KK,k1,3,4,2)
>>k1=Triangle2D3Node_Stiffness(E,NU,t,0,0,1,0,0,1,1)
>>KK=zeros(8,8);
>>KK=Triangle2D3Node_Assemble(KK,k1,1,2,4)
>>KK=Triangle2D3Node_Assemble(KK,k1,3,4,2)
>>k=KK(3:
6,3:
6);
>>p=[500;0;500;0];
>>u=k\p
>>U=[0;0;u;0;0]
>>P=KK*U
>>A=*U'*KK*U-P'*U
(3)计算结果
①应变
CC=
+05*
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
②位移
U=
0
0
3节点力
P=
0
0其中,节点1的支反力为
FRx1500N,FRy1-176.4706N,节点2的支反力为FRx2500N,FRy2176.4706N。
4势能A=
5单元刚度阵
KK=
+05*
00
00
00
00
00
00
00
00
2.基于四节点四边形单元的有限元建模及分析
(1)程序编写说明
Quad2D4Node_Stiffness(E,NU,h,xi,yi,xj,yj,xm,ym,xp,y
p,ID)
该函数计算单元的刚度矩阵,输入模量E,泊松