弹性力学上机题目教学文案.docx
《弹性力学上机题目教学文案.docx》由会员分享，可在线阅读，更多相关《弹性力学上机题目教学文案.docx（27页珍藏版）》请在冰豆网上搜索。
弹性力学上机题目教学文案
弹性力学上机题目
弹性力学上机报告的要求：
从以下六个题目中选取三个，完成每个题目后面的作业即可，上机报告的基本格式另见文件。
弹性力学上机题目
上机题目一主应力的数值计算
实验目的：
1进一步了解主应力、主应变、主方向。
2掌握用MATLAB计算主应力、主应变、主方向的方法
理论介绍：
主应力在某点的应力状态中，如果在某个斜面上的切应力等于零，则该斜面上的正应力称为该点的一个主应力。
主应变在某点的应变状态中，如果在某个斜面上的切应变等于零，则该斜面上的线应变称为该点的一个主应变。
主方向对应于主应力或者主应变的斜面的法线方向称为应力主方向或者应变主方向。
对于同一点，应力主方向与应变主方向是一致的。
二向状态下的主应力：
主方向：
三向状态下的主应力：
主应力满足方向：
式中：
其主方向需进一步讨论得到。
通过研究发现，如果把应力分量写成应力矩阵：
则主应力及其相对应的主方向就是应力矩阵的特征值及其对应的特征向量。
MATLAB中提供了计算矩阵特征值及特征向量的专门函数：
eig（），该函数使用格式如下：
[V,D]=EIG（X）
其中D和V分别为X矩阵的特征值及特征向量。
学生需要再编程将特征向量换算成主方向的方向余弦。
1平面应力状态计算：
算例：
这里的应力正负号采用弹性力学的规定方法。
具体命令如下：
>>A=[-20,10;10,-50]//定义应力矩阵
A=
-2010
10-50
>>[V,D]=eig（A）//调用EIG（）函数，V返回值是特征向量（按列，对应于特征值），也就是方向角与x,y正轴间夹角（方向角）的余弦（方向余弦）；D返回特征值（即主应力大小，按从小到大排序，这与材料力学的规定正好相反）
V=
-0.2898-0.9571
0.9571-0.2898
D=
-53.02780
0-16.9722
>>jiao=rad2deg（atan（V（2,:
）./V（1,:
）））//计算出两个主方向与X轴正向的交角（逆时针为正）
jiao=
-73.155016.8450
>>jiao2=rad2deg（acos（V））//计算出两个主方向的方向角
jiao2=
106.8450163.1550
16.8450106.8450
以上计算结果告诉我们：
该单元体的主应力为：
对应的主方向分别为：
若按方向角描述：
两个主方向的方向角分别是：
有几何知识不难得到主方向的物理位置。
2三向应力状态计算：
具体命令如下：
>>A=[70,40,0;40,30,0;0,0,50]
A=
70400
40300
0050
>>[V,D]=eig（A）
V=
0.52570-0.8507
-0.85070-0.5257
01.00000
D=
5.278600
050.00000
0094.7214
>>jiao2=rad2deg（acos（V））
jiao2=
58.282590.0000148.2825
148.282590.0000121.7175
90.0000090.0000
以上计算结果表明：
该单元体的主应力为：
三个主应力对应的三个主方向的方向角分别是：
作业题：
任选两个进行计算。
上机题目二弹性力学公式的辅助推算
实验目的：
1初步了解MATLAB强大的符号推演功能。
2用MATLAB的符号计算功能，导出弹性力学中的方程。
例1：
导出用应变表示应力的物理方程：
命令如下
>>symsexeyexysxsysxyuvxyEEpsus1s2s3ls1ls2ls3//定义符号变量，ex,ey,exy表示三个应变分量，sx,sy,sxy表示三个应力分量，u,v是位移分量,EE,psu是弹性常数，其余变量是计算所可能用到的中间变量。
>>s1='ex=1/EE*（sx-psu*sy）'//以下三条命令是定义物理方程
s1=
ex=1/EE*（sx-psu*sy）
>>s2='ey=1/EE*（sy-psu*sx）'
s2=
ey=1/EE*（sy-psu*sx）
>>s3='exy=2*（1+psu）/EE*sxy'
s3=
exy=2*（1+psu）/EE*sxy
>>ls1=solve（s1,s2,s3,'sx','sy','sxy'）//利用SOLVE函数从物理方程中反解出应力分量
ls1=
sx:
[1x1sym]
sxy:
[1x1sym]
sy:
[1x1sym]
>>sx=ls1.sx//查看应力分量
sx=
-（EE*ex+EE*ey*psu）/（psu^2-1）
>>sy=ls1.sy
sy=
-（EE*ey+EE*ex*psu）/（psu^2-1）
>>sxy=ls1.sxy
sxy=
（EE*exy）/（2*psu+2）
例2极坐标解答中，若设定应力函数为
，试用MATLAB导出各应力分量。
已知应力分量表示式：
，
相容方程：
。
命令如下：
>>symsFIStsrsfsrf//定义各符号变量：
FI代表
，S是中间变量，t代表半径，srsfsrf是三个应力分量
>>S=dsolve（'t^3*D4f+2*t^2*D3f-9*t*D2f+9*Df=0'）//利用微分方程求解函数求解相容方程。
S=
C2+C3*t^4+C4*t^2+C5/t^2//得到
的表达式。
>>SF=S*cos（2*FI）//应力函数
SF=
cos（2*FI）*（C2+C3*t^4+C4*t^2+C5/t^2）
>>sr=1/t*diff（SF,t）+1/t^2*diff（SF,2,FI）;//以下三式计算三个应力分量，分号表示不显示计算结果（因此结果较乱）。
>>sf=diff（SF,2,t）;
>>srf=-diff（1/t*diff（SF,FI）,t）;
>>sr=simple（sr）;//将较乱的结果化简
>>sf=simple（sf）;
>>srf=simple（srf）;
>>sr//查看应力分量。
sr=
-（2*cos（2*FI）*（C4*t^4+2*C2*t^2+3*C5））/t^4
>>sf
sf=
cos（2*FI）*（2*C4+12*C3*t^2+（6*C5）/t^4）
>>srf
srf=
-（2*sin（2*FI）*（C2*t^2-C4*t^4-3*C3*t^6+3*C5））/t^4
可以将些结果与课本P75页的结果进行对比，并继续算下去。
作业题，请导出一至两个推导过程，例如：
边界条件中待定系数的确定；三维问题中有关问题等。
也可亲自验算以上两例。
上机题目三应力变分法在平面问题中应用
实验目的：
1掌握平面问题的应力变分方法的原理及计算过程。
2用MATLAB实现平面问题的应力变分方法的计算过程，并讨论形函数的收敛性。
平面问题的应力变分方法的基本原理：
平面问题的应变余能可以表示为：
引入应力函数后：
应力变分方程是：
由于面力是给定的，所以上式右端等于零。
即
。
如果应力函数的构成是通过形函数的系数
构成，则余能的变分方程就写成如下的偏微分形式：
例题（P284，11-2）：
正方形薄板，边长为
，如图所示在左右两边受有按抛物线分布的拉力，即
。
试用应力变分法按如下的应力函数求解：
。
%clearall%清理已有变量
cleara1a2a3a4a5
symsaxyFIsxsysxyqa1a2a3a4%定义符号变量，FI为应力函数
FI=q*y^4/12/a^2+q*a^2*（1-x^2/a^2）^2*（1-y^2/a^2）^2*（a1）%将一阶形函数赋值给应力函数
s1=diff（FI,2,y）%以下四行是计算余能积分公式中的被积函数。
s2=diff（FI,2,x）
s3=diff（diff（FI,x）,y）
s4=s1^2+s2^2+2*s3^2
VC=int（s4,'x','-a','a'）%以下两式是进行余能的积分运算
VC=int（VC,'y','-a','a'）
s1=diff（VC,a1）%对余能取变分
A1=solve（s1,a1）%令余能变分等于零，计算待定系数A1
sx=diff（FI,2,y）%计算应力分量
sx=subs（sx,'a1',A1）%将计算出的待定系数回代入应力分量表达式
sx=subs（sx,'x',0）%令x=0
sx=vpa（sx,2）%将应力表达式数值化并取2个有效数字
sx1=simple（sx）%简化并重新命名应力变量，为下一步计算做准备
以下计算3阶
cleara1a2a3a4a5
%clearall
symsaxyFIsxsysxyqa1a2a3a4
FI=q*y^4/12/a^2+q*a^2*（1-x^2/a^2）^2*（1-y^2/a^2）^2*（a1+a2*x^2/a^2+a3*y^2/a^2）%三阶
s1=diff（FI,2,y）
s2=diff（FI,2,x）
s3=diff（diff（FI,x）,y）
s4=s1^2+s2^2+2*s3^2
VC=int（s4,'x','-a','a'）
VC=int（VC,'y','-a','a'）
ss1=diff（VC,a1）
ss2=diff（VC,a2）
ss3=diff（VC,a3）
S=solve（ss1,ss2,ss3,'a1','a2','a3'）
A1=S.a1
A2=S.a2
A3=S.a3
sx=diff（FI,2,y）
sx=subs（sx,'a1',A1）
sx=subs（sx,'a2',A2）
sx=subs（sx,'a3',A3）
sx=subs（sx,'x',0）
sx=vpa（sx,2）
sx3=simple（sx）
以下计算5阶
cleara1a2a3a4a5
%clearall
symsaxyFIsxsysxyqa1a2a3a4a5
FI=q*y^4/12/a^2+q*a^2*（1-x^2/a^2）^2*（1-y^2/a^2）^2*（a1+a2*x^2/a^2+a3*y^2/a^2+a4*x^4/a^4+a5*y^4/a^4）
s1=diff（FI,2,y）
s2=diff（FI,2,x）
s3=diff（diff（FI,x）,y）
s4=s1^2+s2^2+2*s3^2
VC=int（s4,'x','-a','a'）
VC=int（VC,'y','-a','a'）
ss1=diff（VC,a1）
ss2=diff（VC,a2）
ss3=diff（VC,a3）
ss4=diff（VC,a4）
ss5=diff（VC,a5）
S=solve（ss1,ss2,ss3,ss4,ss5,'a1','a2','a3','a4','a5'）
A1=S.a1
A2=S.a2
A3=S.a3
A4=S.a4
A5=S.a5
sx=diff（FI,2,y）
sx=subs（sx,'a1',A1）
sx=subs（sx,'a2',A2）
sx=subs（sx,'a3',A3）
sx=subs（sx,'a4',A4）
sx=subs（sx,'a5',A5）
sx=subs（sx,'x',0）
sx=vpa（sx,2）
sx5=simple（sx）
计算结果：
sx1=0.17*q+（0.49*q*y^2）/a^2
sx3=（q*（0.137*a^4+0.804*a^2*y^2-0.357*y^4））/a^4
sx5=（q*（0.137*a^6+0.7864*a^4*y^2-0.30124*a^2*y^4-0.03528*y^6））/a^6
绘图：
px=（-1:
0.02:
1）
psx1=subs（sx1,'y',px*a）
psx3=subs（sx3,'y',px*a）
psx5=subs（sx5,'y',px*a）
plot（px,psx1/q,'*',px,psx3/q,'.',px,psx5/q,'-'）
结果发现，5阶的结果与3阶几乎没有多少改进，说明本题在3阶时基本上已经收敛到精确解附近。
作业题，11-1，11-3，或者亲自验算上例。
上机题目四位移变分法在梁模型中的应用
实验目的：
1掌握位移变分方法的原理及计算过程。
2用MATLAB实现梁的位移变分方法的计算过程，并讨论形函数的收敛性。
位移变分方法的基本原理：
将位移函数可以设为：
，
，
式中：
满足给定的位移边界条件；
只需要满足边界为零的条件，并且称其为形函数。
目前的情况下，位移函数称为RIZE法形函数。
若形函数在满足上述条件的同时，还能满足给定的应力边界条件，则又可称为伽辽金位移形函数。
将形变势能
用位移函数表示后，位移变分公式就可以写成：
解以上的方程组即可确定出待定系数。
如果研究对象是一个梁，则方程又可以得到简化：
首先位移可设
，其变形能：
RIZE（李滋）方程可写成为：
若忽略体力则：
（a）
此方程组便可确定待定系数
。
其实，也可以采用最小热能原理的方法：
。
式中
是系统总势能，
是系统应变能，
是外力做功的负值，就是外力的势能。
根据最小势能原理，有：
（b）
（a）,（b）方程其实是等价的，读者可以证明。
例1选用合适的形函数形式计算图示简支梁的挠曲线，并与精确解（
）绘图对比。
本例中，边界条件是：
处，
，所以其形函数可以设为：
不难看出
。
命令如下：
clearall
symsb1b2b3vVVqLxyEI%VV是总变形能，v是形函数
v=b1*x*（L-x）
VV=EI/2*int（（diff（v,2,x））^2,x,0,L）-int（q*v,x,0,L）
s1=diff（VV,b1）
B1=solve（s1,b1）
v1=subs（v,'b1',B1）
y=q*x/24/EI*（L^3-2*L*x^2+x^3）
py=subs（y,'x',（0:
0.05:
1）*L）/q/L^4*EI
px=（0:
0.05:
1）
pv1=subs（v1,'x',（0:
0.05:
1）*L）/q/L^4*EI
plot（px,py,'.',px,pv1）
legend（'realsolve','v1'）
以下是二阶形函数计算：
clearb1b2b3
symsb1b2b3vVVqLxyEI%VV是总变形能，v是形函数
v=b1*x*（L-x）+b2*x^2*（L-x）^2
VV=EI/2*int（（diff（v,2,x））^2,x,0,L）-int（q*v,x,0,L）
s1=diff（VV,b1）
s2=diff（VV,b2）
B=solve（s1,s2,b1,b2）
B1=B.b1
B2=B.b2
v2=subs（v,'b1',B1）
v2=subs（v2,'b2',B2）
px=（0:
0.05:
1）
pv2=subs（v2,'x',（0:
0.05:
1）*L）/q/L^4*EI
plot（px,py,'.',px,pv1,px,pv2）
legend（'realsolve','v1','v2'）
从图中可以看出，二阶形函数已经与精确解重合了。
例2仍是上题，如果选用
形函数，试采用位移变分法，或者最小势能原理（两者公式上完全一致）计算待定系数，并与精确的挠曲线进行对比。
只需要把例1中形函数的表达式修改一下，就可以了。
三角级数收敛非常快，本题中只需要一项形函数，就可以达到理想精度。
并且，计算后：
B1=（4*L^4*q）/（EI*pi^5）
B2=0
例3：
悬臂梁受集中力作用下，材料力学中，其挠曲线是
。
试采用幂级数形式的形函数计算其挠曲线，并与精确解对比。
边界条件是，在x=0处，挠度和转角都等于零，所以，可设其形函数为：
编程如下：
clearall
symsb1b2b3vVVFLxyEI%VVÊÇ×ܱäÐÎÄÜ£¬vÊÇÐκ¯Êý
v=b1*x^2
VV=EI/2*int（（diff（v,2,x））^2,x,0,L）-F*subs（v,'x',L）%注意势能的计算与均布载荷不同
s1=diff（VV,b1）
B1=solve（s1,b1）
v1=subs（v,'b1',B1）
y=F*x^2/6/EI*（3*L-x）
py=subs（y,'x',（0:
0.05:
1）*L）/F/L^3*EI
px=（0:
0.05:
1）
pv1=subs（v1,'x',（0:
0.05:
1）*L）/F/L^3*EI
plot（px,py,'.',px,pv1）
legend（'realsolve','v1'）
clearb1b2b3
symsb1b2b3vVVFLxyEI%VVÊÇ×ܱäÐÎÄÜ£¬vÊÇÐκ¯Êý
v=b1*x^2+b2*x^3
VV=EI/2*int（（diff（v,2,x））^2,x,0,L）-F*subs（v,'x',L）
s1=diff（VV,b1）
s2=diff（VV,b2）
B=solve（s1,s2,b1,b2）
B1=B.b1
B2=B.b2
v2=subs（v,'b1',B1）
v2=subs（v2,'b2',B2）
px=（0:
0.05:
1）
pv2=subs（v2,'x',（0:
0.05:
1）*L）/F/L^3*EI
plot（px,py,'.',px,pv1,px,pv2）
legend（'realsolve','v1','v2'）
从图中可以看出，一阶近似时有误差，二阶近似时，已经收敛到精度解。
作业：
亲自验算上例或者自行定义形函数再计算。
或者计算其他形式的梁模型。
上机题目五压杆临界压力计算
实验目的：
1掌握最小势能的原理及计算过程。
2用MATLAB实现最小势能原理计算临界压力计算。
一端固定一端自由：
基本原理（见变分方法讲义中的方法：
）
边界条件为：
。
压杆弯曲变形的变形能：
杆上端由于弯曲变形而下降的高度为：
所以外力的势能：
总势能：
按最小势能原理：
由于两系数
不能同时为零，所以要注上述方程组的系数矩阵的行列式等于零：
，整理后：
所以：
。
精确解为：
，两者误差为：
。
命令：
clearall
symsb1b2vVVFLxyEI%最高算到二阶，因为通过计算发现三阶出现虚根
b=[b1,b2]
i=2%定义的阶数，i=1时是一阶，i=2是二阶
v=0
forii=1:
i
v=v+b（ii）*x^（ii+1）%定义形函数
end
VV=EI/2*int（（diff（v,2,x））^2,x,0,L）-F*int（1/2*（diff（v,x））^2,x,0,L）
forii=1:
i
s（ii）=diff（VV,b（ii））%得到变分方程
end
forii=1:
i
forij=1:
i
AAA（ii,ij）=diff（s（ii）,b（ij））%得到变分的系数矩阵
end
end
dA=det（AAA）
Fcr=solve（dA,F）
Fcrp=vpa（Fcr）
以上程序，当i=1及i=2时，可以分别计算出
Fcr=（3.0*EI）/L^2,及Fcr=（2.4859616991199415410147432028839*EI）/L^2。
作业：
尝试采用三角函数做为形函数计算此例，或者另外其他类型的压杆如：
两端铰支。
上机题目六对径对压薄圆板应力分析
实验目的：
1利用弹性力学解，分析对径受压薄圆板内应力分布情况。
2用MATLAB实现两条直径上应力变化规律。
3用MATLAB验证巴西劈裂实验的拉应力公式。
具体方法：
圆盘的直径为d，在一直径AB的两端点受到一对大小相同，方向相反的集中力F的作用，试求其应力。
弹性力学极坐标解答中，已经给出对径受压圆盘的应力解答为（见PPT）以下三式的叠加：
在直径CD上的应力解答为：
在直径AB上的应力解答为：
利用MATLAB的绘图功能给出应力分布变化规律的直观图。
命令如下：
d=1;
x=d/2*（-1:
0.01:
1）;
y=1-16*d^2*x.^2./（d^2+4*x.^2）.^2;
plot（x,y）
title（'stressonxdirc.onCD'）
ylabel（'stressx（2F/（paid））'）
xlabel（'x（d）'）
d=1;
x=d/2*（-1:
0.01:
1）;
y=1-4*d^4./（d^2+4*x.^2）.^2;
plot（x,y）
title（'stressonydirc.onCD'）
ylabel（'stressy（2F/（paid））'）
xlabel（'x（d）'）
AB直径：
d=1;
x=d/2*（-1:
0.01:
1）;
>>y=1-4*d^2./（d^2-4*x.^2）;
>>plot（x,y）
>>title（'stressonydirc.onAB'）
>>xlabel（'y（d）'）
>>ylabel（'stressy（2F/（paid））'）
升级会员 