有限元作业三角形单元求解.docx
《有限元作业三角形单元求解.docx》由会员分享,可在线阅读,更多相关《有限元作业三角形单元求解.docx(15页珍藏版)》请在冰豆网上搜索。
![有限元作业三角形单元求解.docx](https://file1.bdocx.com/fileroot1/2023-1/31/2a65882f-68d3-4c6f-a41f-44f636cccf6f/2a65882f-68d3-4c6f-a41f-44f636cccf6f1.gif)
有限元作业三角形单元求解
《有限元作业》
年级
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=
+06*
00
00
00
00
k2=
+06*
00
00
00
00
k3=
+06*
00
00
00
00
k4=
+06*
00
00
00
00
调用Triangle2D3Node_Assembly函数,求出总体刚度矩阵
求出的节点位移
U=
0
0
0
0
调用Triangle2D3Node_Stress函数,求出应力,S1、S2、S3、中求出的分别为Sx,Sy,Sxy
S1=
+03*
S2=
+03*
S3=
+03*
S4=
+03*
二、
(1)如图划分四边形单元,工分成四个分别为
(2)如图分别进行编号1、2、3、4、5、6,并建立坐标系
(3)编程进行求解,得出结果,其中假设力P=2000N
调用Quad2D4Node_Stiffness函数,求出单元刚度矩阵
调用Quad2D4Node_Assembly函数,求出求出总体刚度矩阵
求出节点位移
U=
0
0
0
0
调用Quad2D4Node_Stress函数,求出单元应力中的的S1、S2、S3分别为Sx,Sy,Sxy应力分量
S1=
+03*
S2=
+07*
程序附录
一、
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);