利用matlab软件对板进行结构分析Word格式.docx

1、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,

2、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-NU2）*1 NU 0;NU

3、 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;endy=t*A*B*D*B; %根据虚功原理形成单元刚度矩阵该函数将连接节点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-

4、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

5、*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*

6、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

7、,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）;该函数计算在弹性模量为E、泊松比为NU、厚度为t、第一个节点坐标为（xi,yi）、第二个节点坐标为（xj,yj）、第三个节点坐标为（xm,ym）以及单元位移矢量为u时的单元应力。function y=LinearTriangleElementStre

8、sses（E,NU,xi,yi,xj,yj,xm,ym,p,u）if p=1 D=（E/（1-NU*NU）*1 NU 0;elseif p=2y=D*B*u; %返回单元应力该函数根据单元应力矢量sigma计算单元主应力。该函数返回3*1的矢量，其形式为sigma1,sigma2,thetaT。其中sigma1和sigma2为单元的主应力，theta为主应力方向角。function y = LinearTriangleElementPstresses（sigma）R=（sigma（1）+sigma（2）/2;Q=（sigma（1）-sigma（2）/2）2+sigma（3）*sigma（3）;

9、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;该函数用于显示应力云图，根据iStress分别取1、2、3、4、5，分别返回“x方向正应力”、“y方向正应力”、“剪应力”、“第一主应力”、“第二主应力”的应力云图。function PlotStress（iStress）switch iStresscase 1title =x方向正应力; %图形标题显示“x方向正应力”case 2y方向正应力 %图形标题显示“y方向正应力”case 3剪应力 %图形标

10、题显示“剪应力”case 4最大主应力 %图形标题显示“最大主应力”case 5最小主应力 %图形标题显示“最小主应力” endfigure; %创立图形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=

11、gElementcoordinate（ie,2）;gElementcoordinate（ie,4）;gElementcoordinate（ie,6）;c=gElementStress（ie,iStress）;gElementStress（ie,iStress）;gElementStress（ie,iStress）; patch（x,y,c）五MATLAB的command窗口中输入的命令流E=210e6 %弹性模量NU=0.1 %泊松比t=0.1 %板厚k1=LinearTriangleElementStiffness（E,NU,t,0,4,0,2,1,3,1）;k2=LinearTriangl

12、eElementStiffness（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=