利用MATLAB软件对板进行结构剖析.docVIP

一．算例 选取一块矩形形薄板，荷载沿厚度均匀分布，受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。程序如下： 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; %计算单元应变矩阵中的元素 betaj=ym-yi; %计算单元应变矩阵中的元素 betam=yi-yj; %计算单元应变矩阵中的元素 gammai=xm-xj; %计算单元应变矩阵中的元素 gammaj=xi-xm; %计算单元应变矩阵中的元素 gammam=xj-xi; %计算单元应变矩阵中的元素 B=[betai 0 betaj 0 betam 0; %形成单元应变矩阵B 0 gammai 0 gammaj 0 gammam; gammai betai gammaj betaj gammam betam]/(2*A); if p==1 %对于平面应力问题形成弹性矩阵D D=(E/(1-NU^2))*[1 NU 0;NU 1 0;0 0 (1-NU)/2]; elseif p==2 %对于平面应变问题形成弹性矩阵D D=(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)+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-