哈工大材料力学大作业--matlab编程

H a r b i n I n s t i t u t e o f T e c h n o l o g y

材料力学上机作业

课程名称：材料力学

设计题目：应力状态分析

院系：机电学院

班级：

分析者：

学号：

指导教师：***

设计时间：2013年6月18日

哈尔滨工业大学

材料力学上机课设计说明书

一， 设计题目

题目7 应力状态分析 输入：

1. 平面应力状态输入：

x y xy σστ（,,）；某截面方位角α

2. 空间应力状态输入：

,x y z xy yz zx σσστττ（,,,,）

输出： 1. 输出主应力123σσσ（,,）

2. 最大切应力（13

max 132

σσττ-==

）

3.

如为平面应力状态则需要输出方位角α斜截面上的应力αα

στ、及主方向角*σα

4. 画出应力圆示意图

二， 程序计算设计过程

1. 平面应力状态分析

对于任意平面应力状态，有

max min σσ

=

2x y σσ+±

主应力为：

1max 23min ,0,σσσσσ===

并且由 2tan 2xy

x y

στασσ=

-

可求得主应力方向角13σσαα、。 对于任意一个方位角α，有：

=

cos 2sin 222

sin 2cos 22

x y

x y

xy x y

xy αασσσσσατα

σστατα

+++

+-=-

+

从而，输入任意角α，即可求得该截面的应力状态

ααστ、

并且

ααστ、都是关于α的函数，上式即为应力圆的参数方程，参数为α。

将α从0到pi 取一系列的值，则可以求出一系列的ααστ、，在坐标系中找

到对应点，连接即可作出应力圆。 2. 三向应力状态分析

解特征方程 32

1230I I I σσσ-+-=

即可求出主应力123σσσ、、 其中：

123||||||||x y z x

yx y zy z xz xy y yz z zx x x yx zx xy

y zy xz

yz z I I I σσσστστσττστστσστττστττσ=-+⎛⎫⎛⎫⎛⎫

=++

⎪ ⎪ ⎪⎝⎭⎝⎭⎝⎭⎛⎫

⎪

= ⎪ ⎪⎝⎭

再由 13

max 132

σσττ-== 可求得最大切应力。

求解三向应力圆：

三个圆121323C 、C 、C 的圆心分别为：

231312122313,0,0,0222C C C σσσσσσ+++⎛⎫⎛⎫⎛⎫

⎪ ⎪ ⎪⎝⎭⎝⎭⎝⎭

、、

半径非别为：

23

13

12

122313r =

,r =

,r =

2

2

2

σσσσσσ---

由此可以求出三个应力圆的方程，从而作出三向应力圆。

三， 程序代码

reg=input('选择应力状态方式(1或2)：');%1表示平面应力状态，2表示空间应力状态

if reg==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);

if cgmt==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;

for theta=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;