《利用MATLAB软件对板进行结构分析》由会员分享，可在线阅读，更多相关《利用MATLAB软件对板进行结构分析（10页珍藏版）》请在金锄头文库上搜索。
1、一算例选取一块矩形形薄板，荷载沿厚度均匀分布，受 100KN/m 均布荷载，板的各部分尺寸如图所示。板厚 t=0.1，弹性模量 E=100，泊松比=0.1，容重=0KN/myx02m2m2m2m2m142351671413812101192345671314151681012119结构离散化二划分单元划分单元，标出单元号码及节点码，选取坐标，如上图示。三程序选择选择使用 MATLAB 作为开发程序，下面根据程序的结构特点进行叙述四子程序（m 文件）子程序（m 文件）在 MATLAB 程序中具有重要的地位，通过它可以方便的设置自己的函数，这些函数可以在 MATLAB 命令窗口被方便的调用。本例共
2、使用了五个子程序（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
3、,xm,ym,p) 该函数用于计算弹性模量为E、泊松比为NU、厚度为t、第一个节点坐标为(xi,yi)、第二个节点坐标为 (xj,yj)、第三个节点坐标为(xm,ym)时的线性三角形元的单元刚度矩阵。P=1表明函数用于平面应力情 况。P=2表明函数用于平面应变情况。该函数返回6*6的单元刚度矩阵k。程序如下：function y = 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; %计算单元应变矩阵中的元素b
4、etaj=ym-yi; %计算单元应变矩阵中的元素betam=yi-yj; %计算单元应变矩阵中的元素gammai=xm-xj; %计算单元应变矩阵中的元素gammaj=xi-xm; %计算单元应变矩阵中的元素gammam=xj-xi; %计算单元应变矩阵中的元素B=betai 0 betaj 0 betam 0; %形成单元应变矩阵B0 gammai 0 gammaj 0 gammam;gammai betai gammaj betaj gammam betam/(2*A); if p=1 %对于平面应力问题形成弹性矩阵DD=(E/(1-NU2)*1 NU 0;NU 1 0;0 0 (1-N
5、U)/2; elseif p=2 %对于平面应变问题形成弹性矩阵DD=(E/(1+NU)/(1-2*NU)*1-NU NU 0;NU 1-NU 0;0 0 (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。 function y=LinearTriangleAssemble(K,k,i,j,m) K(2*i-1,2*i-1)=K(2*i-1,2*i-1)
6、+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)+
7、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*
8、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(
9、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);end3. LinearTriangleElementStresses(E,NU,xi,yi,xj,yj,xm,ym,p,u）该函数计算在弹性模量为E、泊松比为NU、厚度为t、第一个节点坐标为(x
10、i,yi)、第二个节点坐标为 (xj,yj)、第三个节点坐标为(xm,ym)以及单元位移矢量为u时的单元应力。P=1表明函数用于平面应力 情况。P=2表明函数用于平面应变情况。function y=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=betai 0 betaj 0 betam 0;0 gam
11、mai 0 gammaj 0 gammam;gammai betai gammaj betaj gammam betam/(2*A);if p=1D=(E/(1-NU*NU)*1 NU 0;NU 1 0;0 0 (1-NU)/2;elseif p=2D=(E/(1+NU)/(1-2*NU)*1-NU NU 0;NU 1-NU 0;0 0 (1-2*NU)/2;end y=D*B*u; %返回单元应力4. LinearTriangleElementPstresses(sigma)该函数根据单元应力矢量 sigma 计算单元主应力。该函数返回 3*1 的矢量，其形式为sigma1,sigma2,t
12、hetaT。其中 sigma1 和 sigma2 为单元的主应力，theta 为主应力方向角。function y = 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)该函数用于显示应力云图，根据 iStr
13、ess 分别取 1、2、3、4、5，分别返回“x 方向正应力” 、 “y 方向正应力” 、 “剪应力” 、 “第一主应力” 、 “第二主应力”的应力云图。function PlotStress(iStress)switch iStresscase 1 title =x方向正应力; %图形标题显示“x方向正应力”case 2 title =y方向正应力; %图形标题显示“y方向正应力”case 3 title =剪应力; %图形标题显示“剪应力”case 4 title =最大主应力; %图形标题显示“最大主应力”case 5 title =最小主应力; %图形标题显示“最小主应力”end fi
14、gure; %创立图形axis equal ; %均分坐标轴axis off ; %关闭坐标轴set(gcf, NumberTitle,off);%关闭NumberTitleset(gcf,Name,title) ; %图形标题显示title的返回值for ie=1:1:16 %绘制16个单元的相应应力云图x=gElementcoordinate(ie,1);gElementcoordinate(ie,3);gElementcoordinate(ie,5);y=gElementcoordinate(ie,2);gElementcoordinate(ie,4);gElementcoordinat
15、e(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
