利用matlab软件对板进行结构分析Word格式.docx
《利用matlab软件对板进行结构分析Word格式.docx》由会员分享，可在线阅读，更多相关《利用matlab软件对板进行结构分析Word格式.docx（16页珍藏版）》请在冰豆网上搜索。
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;
%根据虚功原理形成单元刚度矩阵
该函数将连接节点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）;
该函数计算在弹性模量为E、泊松比为NU、厚度为t、第一个节点坐标为（xi,yi）、第二个节点坐标为（xj,yj）、第三个节点坐标为（xm,ym）以及单元位移矢量为u时的单元应力。
functiony=LinearTriangleElementStresses（E,NU,xi,yi,xj,yj,xm,ym,p,u）
ifp==1
D=（E/（1-NU*NU））*[1NU0;
elseifp==2
y=D*B*u;
%返回单元应力
该函数根据单元应力矢量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];
该函数用于显示应力云图，根据iStress分别取1、2、3、4、5，分别返回“x方向正应力”、“y方向正应力”、“剪应力”、“第一主应力”、“第二主应力”的应力云图。
functionPlotStress（iStress）
switchiStress
case1
title='
x方向正应力'
;
%图形标题显示“x方向正应力”
case2
y方向正应力'
%图形标题显示“y方向正应力”
case3
剪应力'
%图形标题显示“剪应力”
case4
最大主应力'
%图形标题显示“最大主应力”
case5
最小主应力'
%图形标题显示“最小主应力”
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）
五．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=
升级会员 