哈工大材料力学大作业matlab编程.docx
《哈工大材料力学大作业matlab编程.docx》由会员分享，可在线阅读，更多相关《哈工大材料力学大作业matlab编程.docx（14页珍藏版）》请在冰点文库上搜索。
哈工大材料力学大作业matlab编程
HarbinInstituteofTechnology
材料力学上机作业
课程名称：
材料力学
设计题目：
应力状态分析
院系：
机电学院
班级：
分析者：
学号：
*******
设计时间：
2013年6月18日
哈尔滨工业大学
材料力学上机课设计说明书
一，设计题目
题目7应力状态分析
输入：
1.平面应力状态输入：
；某截面方位角
2.空间应力状态输入：
输出：
1.输出主应力
2.最大切应力（
）
3.如为平面应力状态则需要输出方位角
斜截面上的应力
及主方向角
4.画出应力圆示意图
二，程序计算设计过程
1.平面应力状态分析
对于任意平面应力状态，有
=
主应力为：
并且由
可求得主应力方向角
。
对于任意一个方位角
，有：
从而，输入任意角
，即可求得该截面的应力状态
并且
都是关于
的函数，上式即为应力圆的参数方程，参数为
。
将
从0到pi取一系列的值，则可以求出一系列的
，在坐标系中找到对应点，连接即可作出应力圆。
2.三向应力状态分析
解特征方程
即可求出主应力
其中：
再由
可求得最大切应力。
求解三向应力圆：
三个圆
的圆心分别为：
半径非别为：
由此可以求出三个应力圆的方程，从而作出三向应力圆。
三，程序代码
reg=input（'选择应力状态方式（1或2）：
'）;%1表示平面应力状态，2表示空间应力状态
ifreg==1%选择平面应力状态分析
%输入已知量，应力单位为MPa，转角单位为rad
cgmx=input（'输入x轴方向正应力cgmx='）;
cgmy=input（'输入y轴方向正应力cgmy='）;
txy=input（'输入切应力txy='）;
%求解主应力、主方向及最大剪应力并输出
cgm1=（cgmx+cgmy）/2+（（（cgmx-cgmy）/2）^2+txy^2）^（1/2）;
cgm2=0;
cgm3=（cgmx+cgmy）/2-（（（cgmx-cgmy）/2）^2+txy^2）^（1/2）;
tm=（cgm1-cgm3）/2;
aerfc=（1/2）*atan（2*txy/（cgmx-cgmy））;
cgmt=（cgmx+cgmy）/2+（cgmx-cgmy）*cos（2*aerfc）/2+txy*sin（2*aerfc）;
ifcgmt==cgm1;
aerfc1=aerfc;
aerfc3=aerfc+pi/2;
else
aerfc3=aerfc;
aerfc1=aerfc+pi/2;
end
display（'主应力为：
'）;
display（cgm1）;
display（cgm2）;
display（cgm3）;
display（'主方向为：
'）;
display（aerfc1）;
display（aerfc3）;
display（'最大切应力为：
'）;
display（tm）;
%求解任意截面上的应力
aerfa=input（'输入截面方位（以弧度表示）aerfa='）;
cgmr=（cgmx+cgmy）/2+（cgmx-cgmy）*cos（2*aerfa）/2+txy*sin（2*aerfa）;
tr=-（cgmx-cgmy）*sin（2*aerfa）/2+txy*cos（2*aerfa）;
display（'截面处应力状况：
'）;
display（'正应力：
'）;
display（cgmr）;
display（'切应力：
'）;
display（tr）;
%求解应力圆并作图
i=0;
fortheta=0:
pi/200:
pi
cgmt=（cgmx+cgmy）/2+（cgmx-cgmy）*cos（2*theta）/2+txy*sin（2*theta）;
tt=-（cgmx-cgmy）*sin（2*theta）/2+txy*cos（2*theta）;
i=i+1;
CG（i）=cgmt;TT（i）=tt;
plot（CG,TT）,axisequal;
title（'应力圆'）;xlabel（'正应力cgm/MPa'）;ylabel（'切应力t/MPa'）;gridon;
end
elseifreg==2%选择三向应力状态分析
%输入已知量，应力单位为MPa，转角单位为rad
cgmx=input（'输入x轴方向正应力cgmx='）;
cgmy=input（'输入y轴方向正应力cgmy='）;
cgmz=input（'输入y轴方向正应力cgmz='）;
txy=input（'输入切应力txy='）;
tyz=input（'输入切应力tyz='）;
tzx=input（'输入切应力tzx='）;
%求解主应力及最大剪应力并输出
I1=cgmx+cgmy+cgmz;
I2=det（[cgmx,txy;txy,cgmy]）+det（[cgmy,tyz;tyz,cgmz]）+det（[cgmz,tzx;tzx,cgmx]）;
I3=det（[cgmx,txy,tzx;txy,cgmy,tyz;tzx,tyz,cgmz]）;
symsx;
ffp=x^3-I1*x^2+I2*x-I3;
cgm=solve（ffp）;
cgm=eval（cgm）;
cgm1=max（cgm
（1）,cgm
（2））;
cgm1=max（cgm1,cgm（3））;
cgm3=min（cgm
（1）,cgm
（2））;
cgm3=min（cgm3,cgm（3））;
cgm2=cgm
（1）+cgm
（2）+cgm（3）-cgm1-cgm3;
tm=（cgm1-cgm3）/2;
display（'主应力为：
'）;
display（cgm1）;
display（cgm2）;
display（cgm3）;
display（'最大切应力为：
'）;
display（tm）;
%求解应力圆并作图
i=0;
r12=（cgm1-cgm2）/2;
r23=（cgm2-cgm3）/2;
r13=（cgm1-cgm3）/2;
x12=（cgm1+cgm2）/2;
x23=（cgm2+cgm3）/2;
x13=（cgm1+cgm3）/2;
fortheta=0:
pi/200:
2*pi
X12=x12+r12*cos（theta）;
Y12=r12*sin（theta）;
X23=x23+r23*cos（theta）;
Y23=r23*sin（theta）;
X13=x13+r13*cos（theta）;
Y13=r13*sin（theta）;
i=i+1;
XX12（i）=X12;YY12（i）=Y12;XX23（i）=X23;YY23（i）=Y23;XX13（i）=X13;YY13（i）=Y13;
plot（XX12,YY12,XX23,YY23,XX13,YY13）,axisequal;
title（'三向应力圆'）;xlabel（'正应力cgm/MPa'）;ylabel（'切应力t/MPa'）;gridon;
text（x12,0,'C12'）;text（x23,0,'C23'）;text（x13,0,'C13'）;
end
else
display（'选择方式错误！
'）;
end
四，程序说明
程序运行后，首先给reg变量赋值，选择应力状态方式，其中reg=1位平面应力状态，reg=2为三向应力状态。
若输入reg为其他值，则会显示“选择方式错误！
”。
1.平面应力状态
若选择平面应力状态，则需要输入：
正应力cgmx、cgmy以及切应力txy。
然后程序就会自动输出三个主应力：
cgm1、cgm2、cgm3以及主应力方向角：
aerfc1、aerfc3，和最大切应力：
tm，进一步输入任意截面方向角：
aerfa，即可求出该截面的正应力：
cgmr，切应力：
tr。
同时作出应力圆的图像。
2.三向应力状态
若选择三向应力状态，则需要输入：
正应力cgmx、cgmy、cgmz，以及切应力txy、tyz、tzx。
然后程序会自动输出三个主应力：
cgm1、cgm2、cgm3，以及最大切应力tm。
同时作出三向应力圆的图像。
五，举例验证
例1.选择平面应力状态，已知
，求主应力、主方向、最大切应力以及
斜截面上的应力，并作出应力圆。
例2.选择空间应力状态，已知：
求主应力及最大切应力，并作出应力圆。
（材料力学第二单元课后第13题）
收获感悟：
做这个材料大作业，虽然花费了很长时间，但是我感到收获很多。
在此过程中，我对MATLAB从一无所知到熟练编辑，甚至接触到对界面编辑，这将成为我一生的财富，MATLAB也必将成为我日后工作的得力助手。
我为能有一次这样锻炼的机会而感到幸运，希望以后还会更多这样的机会充实自己。
对于应力分析部分，我想也会因为这次程序的编写而理解的更加深刻。
不足之处在于因时间关系，没有完全编出界面。
所呈现的界面还不能进行运算。
界面相应程序见电子版。
升级会员 