利用matlab软件对板进行结构分析.docx
- 文档编号:29607302
- 上传时间:2023-07-25
- 格式:DOCX
- 页数:19
- 大小:259.98KB
利用matlab软件对板进行结构分析.docx
《利用matlab软件对板进行结构分析.docx》由会员分享,可在线阅读,更多相关《利用matlab软件对板进行结构分析.docx(19页珍藏版)》请在冰豆网上搜索。
利用matlab软件对板进行结构分析
一.算例
选取一块矩形形薄板,荷载沿厚度均匀分布,受100KN/m均布荷载,板的各部分尺寸如下列图。
板厚t=0.1,弹性模量E=100,泊松比
=0.1,容重
=0
二.划分单元
划分单元,标出单元号码及节点码,选取坐标,如上图示。
三.程序选择
选择使用MATLAB作为开发程序,下面根据程序的构造特点进展表达
四.子程序〔m文件〕
子程序〔m文件〕在MATLAB程序中具有重要的地位,通过它可以方便的设置自己的函数,这些函数可以在MATLAB命令窗口被方便的调用。
本例共使用了五个子程序〔m文件〕,分别是:
1.LinearTriangleElementStiffness(E,NU,t,xi,yi,xj,yj,xm,ym,p)
2.LinearTriangleAssemble(K,k,i,j,m)
3.LinearTriangleElementStresses(E,NU,xi,yi,xj,yj,xm,ym,p,u〕
4.LinearTriangleElementPstresses(sigma)
5.PlotStress(iStress)
详细程序和解释说明如下:
1.LinearTriangleElementStiffness(E,NU,t,xi,yi,xj,yj,xm,ym,p)
该函数用于计算弹性模量为E、泊松比为NU、厚度为t、第一个节点坐标为(xi,yi)、第二个节点坐标为(xj,yj)、第三个节点坐标为(xm,ym)时的线性三角形元的单元刚度矩阵。
P=1说明函数用于平面应力情况。
P=2说明函数用于平面应变情况。
该函数返回6*6的单元刚度矩阵k。
程序如下:
functiony=LinearTriangleElementStiffness(E,NU,t,xi,yi,xj,yj,xm,ym,p)
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;%形成单元应变矩阵B
0gammai0gammaj0gammam;
gammaibetaigammajbetajgammambetam]/(2*A);
ifp==1%对于平面应力问题形成弹性矩阵D
D=(E/(1-NU^2))*[1NU0;NU10;00(1-NU)/2];
elseifp==2%对于平面应变问题形成弹性矩阵D
D=(E/(1+NU)/(1-2*NU))*[1-NUNU0;NU1-NU0;00(1-2*NU)/2];
end
y=t*A*B'*D*B;%根据虚功原理形成单元刚度矩阵
2.LinearTriangleAssemble(K,k,i,j,m)
该函数将连接节点i和节点j的线性三角形元的单元刚度矩阵k集成到整体刚度矩阵K。
每集成一个单元,该函数都返回2n*2n的整体刚度矩阵K。
functiony=LinearTriangleAssemble(K,k,i,j,m)
K(2*i-1,2*i-1)=K(2*i-1,2*i-1)+k(1,1);
K(2*i-1,2*i)=K(2*i-1,2*i)+k(1,2);
K(2*i-1,2*j-1)=K(2*i-1,2*j-1)+k(1,3);
K(2*i-1,2*j)=K(2*i-1,2*j)+k(1,4);
K(2*i-1,2*m-1)=K(2*i-1,2*m-1)+k(1,5);
K(2*i-1,2*m)=K(2*i-1,2*i-1)+k(2,1);
K(2*i,2*i-1)=K(2*i,2*i-1)+k(2,1);
K(2*i,2*i)=K(2*i,2*i)+k(2,2);
K(2*i,2*j-1)=K(2*j,2*j-1)+k(2,3);
K(2*i,2*j)=K(2*i,2*j)+k(2,4);
K(2*i,2*m-1)=K(2*i,2*m-1)+k(2,5);
K(2*i,2*m)=K(2*i,2*m)+k(2,6);
K(2*j-1,2*i-1)=K(2*j-1,2*i-1)+k(3,1);
K(2*j-1,2*j)=K(2*j-1,2*i)+k(3,2);
K(2*j-1,2*j-1)=K(2*j-1,2*j-1)+k(3,3);
K(2*j-1,2*j)=K(2*j-1,2*j)+k(3,4);
K(2*j-1,2*m-1)=K(2*j-1,2*m-1)+k(3,5);
K(2*j-1,2*m)=K(2*j-1,2*m)+k(3,6);
K(2*j,2*i-1)=K(2*j,2*i-1)+k(4,1);
K(2*j,2*i)=K(2*j,2*i)+k(4,2);
K(2*j,2*j-1)=K(2*j,2*j-1)+k(4,3);
K(2*j,2*j)=K(2*j,2*j)+k(4,4);
K(2*j,2*m-1)=K(2*j,2*m-1)+k(4,5);
K(2*j,2*m)=K(2*j,2*m)+k(4,6);
K(2*m-1,2*i-1)=K(2*m-1,2*i-1)+k(5,1);
K(2*m-1,2*i)=K(2*m-1,2*i)+k(5,2);
K(2*m-1,2*j-1)=K(2*m-1,2*j-1)+k(5,3);
K(2*m-1,2*j)=K(2*m-1,2*j)+k(5,4);
K(2*m-1,2*m-1)=K(2*m-1,2*m-1)+k(5,5);
K(2*m-1,2*m)=K(2*m-1,2*m)+k(5,6);
K(2*m,2*i-1)=K(2*m,2*i-1)+k(6,1);
K(2*m,2*i)=K(2*m,2*i)+k(6,2);
K(2*m,2*j-1)=K(2*m,2*j-1)+k(6,3);
K(2*m,2*)=K(2*i-1,2*i-1)+k(1,1);
end
3.LinearTriangleElementStresses(E,NU,xi,yi,xj,yj,xm,ym,p,u〕
该函数计算在弹性模量为E、泊松比为NU、厚度为t、第一个节点坐标为(xi,yi)、第二个节点坐标为(xj,yj)、第三个节点坐标为(xm,ym)以及单元位移矢量为u时的单元应力。
P=1说明函数用于平面应力情况。
P=2说明函数用于平面应变情况。
functiony=LinearTriangleElementStresses(E,NU,xi,yi,xj,yj,xm,ym,p,u)
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);
ifp==1
D=(E/(1-NU*NU))*[1NU0;NU10;00(1-NU)/2];
elseifp==2
D=(E/(1+NU)/(1-2*NU))*[1-NUNU0;NU1-NU0;00(1-2*NU)/2];
end
y=D*B*u;%返回单元应力
4.LinearTriangleElementPstresses(sigma)
该函数根据单元应力矢量sigma计算单元主应力。
该函数返回3*1的矢量,其形式为[sigma1,sigma2,theta]T。
其中sigma1和sigma2为单元的主应力,theta为主应力方向角。
functiony=LinearTriangleElementPstresses(sigma)
R=(sigma
(1)+sigma
(2))/2;
Q=((sigma
(1)-sigma
(2))/2)^2+sigma(3)*sigma(3);
M=2*sigma(3)/(sigma
(1)-sigma
(2));
s1=R+sqrt(Q);
s2=R-sqrt(Q);
theta=(atan(M)/2)*180/pi;
y=[s1;s2;theta];
5.PlotStress(iStress)
该函数用于显示应力云图,根据iStress分别取1、2、3、4、5,分别返回“x方向正应力〞、“y方向正应力〞、“剪应力〞、“第一主应力〞、“第二主应力〞的应力云图。
functionPlotStress(iStress)
switchiStress
case1
title='x方向正应力';%图形标题显示“x方向正应力〞
case2
title='y方向正应力';%图形标题显示“y方向正应力〞
case3
title='剪应力';%图形标题显示“剪应力〞
case4
title='最大主应力';%图形标题显示“最大主应力〞
case5
title='最小主应力';%图形标题显示“最小主应力〞
end
figure;%创立图形
axisequal;%均分坐标轴
axisoff;%关闭坐标轴
set(gcf,'NumberTitle','off');%关闭NumberTitle
set(gcf,'Name',title);%图形标题显示title的返回值
forie=1:
1:
16%绘制16个单元的相应应力云图
x=[gElementcoordinate(ie,1);gElementcoordinate(ie,3);gElementcoordinate(ie,5)];
y=[gElementcoordinate(ie,2);gElementcoordinate(ie,4);gElementcoordinate(ie,6)];
c=[gElementStress(ie,iStress);gElementStress(ie,iStress);gElementStress(ie,iStress)];
patch(x,y,c)
end
五.MATLAB的command窗口中输入的命令流
E=210e6%弹性模量
NU=0.1%泊松比
t=0.1%板厚
k1=LinearTriangleElementStiffness(E,NU,t,0,4,0,2,1,3,1);
k2=LinearTriangleElementStiffness(E,NU,t,0,4,1,3,2,4,1);
k3=LinearTriangleElementStiffness(E,NU,t,2,4,1,3,2,2,1);
k4=LinearTriangleElementStiffness(E,NU,t,1,3,0,2,2,2,1);
k5=LinearTriangleElementStiffness(E,NU,t,2,4,2,2,3,3,1);k6=LinearTriangleElementStiffness(E,NU,t,2,4,3,3,4,4,1);
利用LinearTriangleElementStiffness
子函数计算单元刚度矩阵
k7=LinearTriangleElementStiffness(E,NU,t,3,3,4,2,4,4,1);
k8=LinearTriangleElementStiffness(E,NU,t,3,3,2,2,4,2,1);
k9=LinearTriangleElementStiffness(E,NU,t,2,2,2,0,3,1,1);
k10=LinearTriangleElementStiffness(E,NU,t,2,2,3,1,4,2,1);
k11=LinearTriangleElementStiffness(E,NU,t,3,1,4,0,4,2,1);
k12=LinearTriangleElementStiffness(E,NU,t,3,1,2,0,4,0,1);
k13=LinearTriangleElementStiffness(E,NU,t,4,4,4,2,5,3,1);
k14=LinearTriangleElementStiffness(E,NU,t,4,4,5,3,6,4,1);
k15=LinearTriangleElementStiffness(E,NU,t,5,3,6,2,6,4,1);
k16=LinearTriangleElementStiffness(E,NU,t,4,2,6,2,5,3,1);
K=zeros(28,28);建立一个28*28的零矩阵
K=LinearTriangleAssemble(K,k1,1,4,3);
K=LinearTriangleAssemble(K,k2,1,3,2);
K=LinearTriangleAssemble(K,k3,2,3,5);
K=LinearTriangleAssemble(K,k4,3,4,5);
K=LinearTriangleAssemble(K,k5,2,5,6);
K=LinearTriangleAssemble(K,k6,2,6,8);
利用LinearTriangleAssemble子函数形成总刚,K=LinearTriangleAssemble(K,k16,2,3,15)
即为总刚
K=LinearTriangleAssemble(K,k7,6,7,8);
K=LinearTriangleAssemble(K,k8,6,5,7);
K=LinearTriangleAssemble(K,k9,5,10,9);
K=LinearTriangleAssemble(K,k10,5,9,7);
K=LinearTriangleAssemble(K,k11,9,11,7);
K=LinearTriangleAssemble(K,k12,9,10,11);
K=LinearTriangleAssemble(K,k13,8,7,12);
K=LinearTriangleAssemble(K,k14,8,12,14);
K=LinearTriangleAssemble(K,k15,12,13,14);
K=LinearTriangleAssemble(K,k16,7,13,12)
引入边界条件,在相应对角元素上乘大数,消除总刚的奇异性
K(1,1)=K(1,1)*1e15;
K(2,2)=K(2,2)*1e15;
K(19,19)=K(19,19)*1e15;
K(20,20)=K(20,20)*1e15;
K(21,21)=K(21,21)*1e15;
K(22,22)=K(22,22)*1e15;
F=[0;-100;0;-200;0;0;0;0;0;0;0;0;0;0;0;-200;0;0;0;0;0;0;0;0;0;0;0;-100];
U=K\F;计算得出各个节点位移,形成节点位向量
u1=[U
(1);U
(2);U(7);U(8);U(5);U(6)];
u2=[U
(1);U
(2);U(5);U(6);U(3);U(4)];
u3=[U(3);U(4);U(5);U(6);U(9);U(10)];
u4=[U(5);U(6);U(7);U(8);U(9);U(10)];
u5=[U(3);U(4);U(9);U(10);U(11);U(12)];
将各个三角形单元的6个节点位移从节点位移向量U中取出,形成各个三角形单元的节点位移向量,作为备用
u6=[U(3);U(4);U(11);U(12);U(15);U(16)];
u7=[U(11);U(12);U(13);U(14);U(15);U(16)];
u8=[U(11);U(12);U(9);U(10);U(13);U(14)];
u9=[U(9);U(10);U(19);U(20);U(17);U(18)];
u10=[U(9);U(10);U(17);U(18);U(13);U(14)];
u11=[U(17);U(18);U(21);U(22);U(13);U(14)];
u12=[U(17);U(18);U(19);U(20);U(21);U(22)];
u13=[U(15);U(16);U(13);U(14);U(23);U(24)];
u14=[U(15);U(16);U(23);U(24);U(27);U(28)];
u15=[U(23);U(24);U(25);U(26);U(27);U(28)];
u16=[U(13);U(14);U(25);U(26);U(23);U(24)]
sigma1=LinearTriangleElementStresses(E,NU,0,4,0,2,1,3,1,u1);
sigma2=LinearTriangleElementStresses(E,NU,0,4,1,3,2,4,1,u2);
sigma3=LinearTriangleElementStresses(E,NU,2,4,1,3,2,2,1,u3);
sigma4=LinearTriangleElementStresses(E,NU,1,3,0,2,2,2,1,u4);
sigma5=LinearTriangleElementStresses(E,NU,2,4,2,2,3,3,1,u5);
sigma6=LinearTriangleElementStresses(E,NU,2,4,3,3,4,4,1,u6);
sigma7=LinearTriangleElementStresses(E,NU,3,3,4,2,4,4,1,u7);
sigma8=LinearTriangleElementStresses(E,NU,3,3,2,2,4,2,1,u8);
sigma9=LinearTriangleElementStresses(E,NU,2,2,2,0,3,1,1,u9);
sigma10=LinearTriangleElementStresses(E,NU,2,2,3,1,4,2,1,u10);
sigma11=LinearTriangleElementStresses(E,NU,3,1,4,0,4,2,1,u11);
sigma12=LinearTriangleElementStresses(E,NU,3,1,2,0,4,0,1,u12);
sigma13=LinearTriangleElementStresses(E,NU,4,4,4,2,5,3,1,u13);
sigma14=LinearTriangleElementStresses(E,NU,4,4,5,3,6,4,1,u14);
sigma15=LinearTriangleElementStresses(E,NU,5,3,6,2,6,4,1,u15);
sigma16=LinearTriangleElementStresses(E,NU,4,2,6,2,5,3,1,u16);
s1=LinearTriangleElementPStresses(sigma1);
s2=LinearTriangleElementPStresses(sigma2);
s3=LinearTriangleElementPStresses(sigma3);
s4=LinearTriangleElementPStresses(sigma4);
s5=LinearTriangleElementPStresses(sigma5);
利用LinearTriangleElementPStresses
子函数计算每个单元的主应力,返回值的意义是,如:
S5=[第一主应力;第二主应力;主应力方向角]
s6=LinearTriangleElementPStresses(sigma6);
s7=LinearTriangleElementPStresses(sigma7);
s8=LinearTriangleElementPStresses(sigma8);
s9=LinearTriangleElementPStresses(sigma9);
s10=LinearTriangleElementPStresses(sigma10);
s11=LinearTriangleElementPStresses(sigma11);
s12=LinearTriangleElementPStresses(sigma12);
s13=LinearTriangleElementPStresses(sigma13);
s14=LinearTriangleElementPStresses(sigma14);
s15=LinearTriangleElementPStresses(sigma15);
s16=LinearTriangleElementPStresses(sigma16);
yipslong1=LinearTriangleElementStresses(E,NU,0,4,0,2,1,3,2,u1);
yipslong2=LinearTriangleElementStresses(E,NU,0,4,1,3,2,4,2,u2);
yipslong3=LinearTriangleElementStresses(E,NU,2,4,1,3,2,2,2,u3);
yipslong4=LinearTriangleElementStresses(E,NU,1,3,0,2,2,2,2,u4);
yipslong5=LinearTriangleElementStresses(E,NU,2,4,2,2,3,3,2,u5);
利用LinearTriangleElementStresses子函数计算每个单元的应变
yipslong6=LinearTriangleElementStresses(E,NU,2,4,3,3,4,4,2,u6);
yipslong7=LinearTriangleElementStresses(E,NU,3,3,4,2,4,4,2,u7);
yipslong8=LinearTriangleElementStresses(E,NU,3,3,2,2,4,2,2,u8);
yipslong9=LinearTriangleElementStresses(E,NU,2,2,2,0,3,1,2,u9);
yipslong10=LinearTriangleElementStresses(E,NU,2,2,3,1,4,2,2,u10);
yipslong11=LinearTriangleElementStresses(E,NU,3,1,4,0,4,2,2,u11);
yipslong12=LinearTriangleElementStresses(E,NU,3,1,2,0,4,0,2,u12);
yipslong13=LinearTriangle
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 利用 matlab 软件 进行 结构 分析