有限元作业三角形单元求解.docx
《有限元作业三角形单元求解.docx》由会员分享，可在线阅读，更多相关《有限元作业三角形单元求解.docx（17页珍藏版）》请在冰豆网上搜索。
有限元作业三角形单元求解
《有限元作业》
年级
2015级
学院
机电工程学院
专业名称
班级学号
学生姓名
2016年05月
如下图所示为一受集中力P作用的结构，弹性模量E为常量，泊松比V=1/6,厚度为I=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.14290.8571-5.14292.1429
-3.00007.28572.1429-5.14290.8571-2.1429
-2.14292.14292.142900-2.1429
0.8571-5.142905.1429-0.85710
-5.14290.85710-0.85715.14290
2.1429-2.1429-2.1429002.1429
k2=
1.0e+06*
5.14290-5.14290.85710-0.8571
02.14292.1429-2.1429-2.14290
-5.14292.14297.2857-3.0000-2.14290.8571
0.8571-2.1429-3.00007.28572.1429-5.1429
0-2.1429-2.14292.14292.14290
-0.857100.8571-5.142905.1429
k3=
1.0e+06*
2.14290-2.1429-2.142902.1429
05.1429-0.8571-5.14290.85710
-2.1429-0.85717.28573.0000-5.1429-2.1429
-2.1429-5.14293.00007.2857-0.8571-2.1429
00.8571-5.1429-0.85715.14290
2.14290-2.1429-2.142902.1429
k4=
1.0e+06*
2.14290-2.1429-2.142902.1429
05.1429-0.8571-5.14290.85710
-2.1429-0.85717.28573.0000-5.1429-2.1429
-2.1429-5.14293.00007.2857-0.8571-2.1429
00.8571-5.1429-0.85715.14290
2.14290-2.1429-2.142902.1429
调用Triangle2D3Node_Assembly函数，求出总体刚度矩阵
求出的节点位移
U=
0
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
S3=
1.0e+03*
1.8907
-1.0601
2.1093
S4=
1.0e+03*
-1.8907
2.1093
1.8907
二、
（1）如图划分四边形单元，工分成四个分别为①②
（2）如图分别进行编号1、2、3、4、5、6，并建立坐标系
（3）编程进行求解，得出结果，其中假设力P=2000N
调用Quad2D4Node_Stiffness函数，求出单元刚度矩阵
调用Quad2D4Node_Assembly函数，求出求出总体刚度矩阵
求出节点位移
U=
0
0
0
0
0.0012
0.0017
-0.0012
0.0017
0.0016
0.0049
-0.0017
0.0052
调用Quad2D4Node_Stress函数，求出单元应力中的的S1、S2、S3分别为Sx,Sy,Sxy应力分量
S1=
1.0e+03*
0.0000
-0.2478
2.0000
S2=
1.0e+07*
0.6856
4.1135
-1.7137
程序附录
一、
1、三角形单元总程序：
E=1e7;
NU=1/6;
t=1;
ID=1;
%调用Triangle2D3Node_Stiffness函数，求出单元刚度矩阵
k1=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（4）;U（5）;U（6）];
u2=[U（3）;U（4）;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
u1=[U
（1）;U
（2）;U（3）;U（4）;U（5）;U（6）];
u2=[U（3）;U（4）;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（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;
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];
end
stress=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（4）;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*（1-s））/4;
B1=[a*（t-1）/4-b*（s-1）/40;0c*（s-1）/4-d*（t-1）/4;
c*（s-1）/4-d*（t-1）/4a*（t-1）/4-b*（s-1）/4];
B2=[a*（1-t）/4-b*（-1-s）/40;0c*（-1-s）/4-d*（1-t）/4;
c*（-1-s）/4-d*（1-t）/4a*（1-t）/4-b*（-1-s）/4];
B3=[a*（t+1）/4-b*（s+1）/40;0c*（s+1）/4-d*（t+1）/4;
c*（s+1）/4-d*（t+1）/4a*（t+1）/4-b*（s+1）/4];
B4=[a*（-1-t）/4-b*（1-s）/40;0c*（1-s）/4-d*（-1-t）/4;
c*（1-s）/4-d*（-1-t）/4a*（-1-t）/4-b*（1-s）/4];
Bfirst=[B1B2B3B4];
Jfirst=[01-tt-ss-1;t-10s+1-s-t;
s-t-s-10t+1;1-ss+t-t-10];
J=[xixjxmxp]*Jfirst*[yi;yj;ym;yp]/8;
B=Bfirst/J;
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
BD=J*transpose（B）*D*B;
r=int（int（BD,t,-1,1）,s,-1,1）;
z=h*r;
k=double（z）;
3、求总体刚度矩阵程序
functionz=Quad2D4Node_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;
4、求应力程序
functionstress=Quad2D4Node_Stress（E,NU,xi,yi,xj,yj,xm,ym,xp,yp,u,ID）
%该函数计算单元的应力
%输入弹性模量E，泊松比NU，厚度h，
%输入4个节点i、j、m、p的坐标xi,yi,xj,yj,xm,ym,xp,yp,
%输入平面问题性质指示参数ID（1为平面应力，2为平面应变）
%输入单元的位移列阵u（8X1）
%输出单元的应力stress（3X1）
%由于它为常应力单元，则单元的应力分量为Sx,Sy,Sxy
%---------------------------------------------------------------
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*（1-s））/4;
B1=[a*（t-1）/4-b*（s-1）/40;0c*（s-1）/4-d*（t-1）/4;
c*（s-1）/4-d*（t-1）/4a*（t-1）/4-b*（s-1）/4];
B2=[a*（1-t）/4-b*（-1-s）/40;0c*（-1-s）/4-d*（1-t）/4;
c*（-1-s）/4-d*（1-t）/4a*（1-t）/4-b*（-1-s）/4];
B3=[a*（t+1）/4-b*（s+1）/40;0c*（s+1）/4-d*（t+1）/4;
c*（s+1）/4-d*（t+1）/4a*（t+1）/4-b*（s+1）/4];
B4=[a*（-1-t）/4-b*（1-s）/40;0c*（1-s）/4-d*（-1-t）/4;
c*（1-s）/4-d*（-1-t）/4a*（-1-t）/4-b*（1-s）/4];
Bfirst=[B1B2B3B4];
Jfirst=[01-tt-ss-1;t-10s+1-s-t;
s-t-s-10t+1;1-ss+t-t-10];
J=[xixjxmxp]*Jfirst*[yi;yj;ym;yp]/8;
B=Bfirst/J;
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
str1=D*B*u;
str2=subs（str1,{s,t},{0,0}）;
stress=double（str2）;
升级会员 