弹性力学上机题目

1、弹性力学上机报告的要求：从以下六个题目中选取三个，完成每个题目后面的作业即可，上机报告的基本格式另见文件。弹性力学上机题目上机题目一 主应力的数值计算实验目的：1 进一步了解主应力、主应变、主方向。2 掌握用MATLAB计算主应力、主应变、主方向的方法理论介绍：主应力 在某点的应力状态中，如果在某个斜面上的切应力等于零，则该斜面上的正应力称为该点的一个主应力。主应变 在某点的应变状态中，如果在某个斜面上的切应变等于零，则该斜面上的线应变称为该点的一个主应变。主方向 对应于主应力或者主应变的斜面的法线方向称为应力主方向或者应变主方向。对于同一点，应力主方向与应变主方向是一致的。二向状态下的主应力

2、：主方向：三向状态下的主应力：主应力满足方向：式中：其主方向需进一步讨论得到。通过研究发现，如果把应力分量写成应力矩阵：则主应力及其相对应的主方向就是应力矩阵的特征值及其对应的特征向量。MATLAB中提供了计算矩阵特征值及特征向量的专门函数：eig()，该函数使用格式如下：V,D = EIG(X) 其中D和V分别为X矩阵的特征值及特征向量。学生需要再编程将特征向量换算成主方向的方向余弦。1 平面应力状态计算：算例：这里的应力正负号采用弹性力学的规定方法。具体命令如下：>> A=-20,10;10, -50 /定义应力矩阵A = -20 10 10 -50>> V,D=e

3、ig(A) /调用EIG（）函数，V返回值是特征向量（按列，对应于特征值），也就是方向角与x,y正轴间夹角（方向角）的余弦（方向余弦）；D返回特征值（即主应力大小，按从小到大排序，这与材料力学的规定正好相反）V = -0.2898 -0.9571 0.9571 -0.2898D = -53.0278 0 0 -16.9722>> jiao=rad2deg(atan(V(2,:)./V(1,:) /计算出两个主方向与X轴正向的交角（逆时针为正）jiao = -73.1550 16.8450>> jiao2=rad2deg(acos(V) /计算出两个主方向的方向角jiao

4、2 = 106.8450 163.1550 16.8450 106.8450以上计算结果告诉我们：该单元体的主应力为：对应的主方向分别为：若按方向角描述：两个主方向的方向角分别是：有几何知识不难得到主方向的物理位置。2 三向应力状态计算：具体命令如下：>> A=70,40,0;40,30,0;0,0,50A = 70 40 0 40 30 0 0 0 50>> V,D=eig(A)V = 0.5257 0 -0.8507 -0.8507 0 -0.5257 0 1.0000 0D = 5.2786 0 0 0 50.0000 0 0 0 94.7214>>

5、jiao2=rad2deg(acos(V)jiao2 = 58.2825 90.0000 148.2825 148.2825 90.0000 121.7175 90.0000 0 90.0000以上计算结果表明：该单元体的主应力为：三个主应力对应的三个主方向的方向角分别是：作业题：任选两个进行计算。 上机题目二 弹性力学公式的辅助推算实验目的：1 初步了解MATLAB强大的符号推演功能。2 用MATLAB的符号计算功能，导出弹性力学中的方程。例1：导出用应变表示应力的物理方程：命令如下>> syms ex ey exy sx sy sxy u v x y EE psu s1 s2

6、s3 ls1 ls2 ls3 /定义符号变量，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)

7、/EE*sxy>> ls1=solve(s1,s2,s3,'sx','sy','sxy') /利用SOLVE函数从物理方程中反解出应力分量ls1 = sx: 1x1 sym sxy: 1x1 sym sy: 1x1 sym>> sx=ls1.sx /查看应力分量 sx = -(EE*ex + EE*ey*psu)/(psu2 - 1) >> sy=ls1.sy sy = -(EE*ey + EE*ex*psu)/(psu2 - 1) >> sxy=ls1.sxy sxy = (EE*exy)/(2*

8、psu + 2)例2极坐标解答中，若设定应力函数为，试用MATLAB导出各应力分量。已知应力分量表示式：，相容方程：。命令如下：>> syms FI S t sr sf srf /定义各符号变量：FI代表，S是中间变量，t代表半径，sr sf srf是三个应力分量>> S=dsolve('t3*D4f+2*t2*D3f-9*t*D2f+9*Df=0') /利用微分方程求解函数求解相容方程。 S = C2 + C3*t4 + C4*t2 + C5/t2 /得到的表达式。 >> SF=S*cos(2*FI) /应力函数 SF = cos(2*FI

9、)*(C2 + C3*t4 + C4*t2 + C5/t2) >> sr=1/t*diff(SF,t)+1/t2*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*

10、t4 + 2*C2*t2 + 3*C5)/t4 >> sf sf = cos(2*FI)*(2*C4 + 12*C3*t2 + (6*C5)/t4) >> srf srf = -(2*sin(2*FI)*(C2*t2 - C4*t4 - 3*C3*t6 + 3*C5)/t4可以将些结果与课本P75页的结果进行对比，并继续算下去。作业题，请导出一至两个推导过程，例如：边界条件中待定系数的确定；三维问题中有关问题等。也可亲自验算以上两例。上机题目三 应力变分法在平面问题中应用实验目的：1 掌握平面问题的应力变分方法的原理及计算过程。2 用MATLAB实现平面问题的应力变分方

11、法的计算过程，并讨论形函数的收敛性。平面问题的应力变分方法的基本原理：平面问题的应变余能可以表示为：引入应力函数后：应力变分方程是：由于面力是给定的，所以上式右端等于零。即。如果应力函数的构成是通过形函数的系数构成，则余能的变分方程就写成如下的偏微分形式：例题（P284，11-2）：正方形薄板，边长为，如图所示在左右两边受有按抛物线分布的拉力，即。试用应力变分法按如下的应力函数求解：。%clear all %清理已有变量clear a1 a2 a3 a4 a5syms a x y FI sx sy sxy q a1 a2 a3 a4 %定义符号变量，FI为应力函数FI=q*y4/12/a2+q

12、*a2*(1-x2/a2)2*(1-y2/a2)2*(a1) %将一阶形函数赋值给应力函数s1=diff(FI,2,y) %以下四行是计算余能积分公式中的被积函数。s2=diff(FI,2,x)s3=diff(diff(FI,x),y)s4=s12+s22+2*s32VC=int(s4,'x','-a','a') %以下两式是进行余能的积分运算VC=int(VC,'y','-a','a')s1=diff(VC,a1) %对余能取变分A1=solve(s1,a1) %令余能变分等于零，计算待定系数A1

13、sx=diff(FI,2,y) %计算应力分量sx=subs(sx,'a1',A1) %将计算出的待定系数回代入应力分量表达式sx=subs(sx,'x',0) %令x=0sx=vpa(sx,2) %将应力表达式数值化并取2个有效数字sx1=simple(sx) %简化并重新命名应力变量，为下一步计算做准备以下计算3阶clear a1 a2 a3 a4 a5%clear allsyms a x y FI sx sy sxy q a1 a2 a3 a4FI=q*y4/12/a2+q*a2*(1-x2/a2)2*(1-y2/a2)2*(a1+a2*x2/a2+a3*

14、y2/a2) %三阶s1=diff(FI,2,y)s2=diff(FI,2,x)s3=diff(diff(FI,x),y)s4=s12+s22+2*s32VC=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.a1A2=S.a2A3=S.a3sx=d

15、iff(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阶clear a1 a2 a3 a4 a5%clear allsyms a x y FI sx sy sxy q a1 a2 a3 a4 a5FI=q*y4/12/a2+q*a2*(1-x2/a2)2*(1-y2/a2)2*(a1+a2*x2/a2+a3*y2/a2+a4*x4/a4+a5*y4/a4)

16、s1=diff(FI,2,y)s2=diff(FI,2,x)s3=diff(diff(FI,x),y)s4=s12+s22+2*s32VC=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

17、9;,'a4','a5')A1=S.a1A2=S.a2A3=S.a3A4=S.a4A5=S.a5sx=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*y2)/a2sx

18、3 =(q*(0.137*a4 + 0.804*a2*y2 - 0.357*y4)/a4sx5 =(q*(0.137*a6 + 0.7864*a4*y2 - 0.30124*a2*y4 - 0.03528*y6)/a6绘图：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阶的结果

19、与3阶几乎没有多少改进，说明本题在3阶时基本上已经收敛到精确解附近。作业题，11-1，11-3，或者亲自验算上例。上机题目四 位移变分法在梁模型中的应用实验目的：1 掌握位移变分方法的原理及计算过程。2 用MATLAB实现梁的位移变分方法的计算过程，并讨论形函数的收敛性。位移变分方法的基本原理：将位移函数可以设为：，式中：满足给定的位移边界条件；只需要满足边界为零的条件，并且称其为形函数。目前的情况下，位移函数称为RIZE法形函数。若形函数在满足上述条件的同时，还能满足给定的应力边界条件，则又可称为伽辽金位移形函数。将形变势能用位移函数表示后，位移变分公式就可以写成：解以上的方程组即可确定出待

20、定系数。如果研究对象是一个梁，则方程又可以得到简化：首先位移可设，其变形能：RIZE（李滋）方程可写成为：若忽略体力则： （a）此方程组便可确定待定系数。其实，也可以采用最小热能原理的方法：。式中是系统总势能，是系统应变能，是外力做功的负值，就是外力的势能。根据最小势能原理，有： (b)(a),(b)方程其实是等价的，读者可以证明。例1选用合适的形函数形式计算图示简支梁的挠曲线，并与精确解（）绘图对比。本例中，边界条件是：处，所以其形函数可以设为：不难看出。命令如下：clear allsyms b1 b2 b3 v VV q L x y EI %VV是总变形能，v是形函数v=b1*x*(L-x

21、)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*(L3-2*L*x2+x3)py=subs(y,'x',(0:0.05:1)*L)/q/L4*EIpx=(0:0.05:1)pv1=subs(v1,'x',(0:0.05:1)*L)/q/L4*EIplot(px,py,'.',px,pv1)legend('real solve','v1')

22、以下是二阶形函数计算：clear b1 b2 b3syms b1 b2 b3 v VV q L x y EI % VV是总变形能，v是形函数v=b1*x*(L-x)+b2*x2*(L-x)2VV=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.b1B2=B.b2v2=subs(v,'b1',B1)v2=subs(v2,'b2',B2)px=(0:0.05:1)pv2=subs(v2,'x',(0:0.0

23、5:1)*L)/q/L4*EIplot(px,py,'.',px,pv1,px,pv2)legend('real solve','v1','v2')从图中可以看出，二阶形函数已经与精确解重合了。例2 仍是上题，如果选用形函数，试采用位移变分法，或者最小势能原理（两者公式上完全一致）计算待定系数，并与精确的挠曲线进行对比。只需要把例1中形函数的表达式修改一下，就可以了。三角级数收敛非常快，本题中只需要一项形函数，就可以达到理想精度。并且，计算后：B1 =(4*L4*q)/(EI*pi5)B2 =0例3：悬臂梁受集中力作用下，材料力学

24、中，其挠曲线是。试采用幂级数形式的形函数计算其挠曲线，并与精确解对比。边界条件是，在x=0处，挠度和转角都等于零，所以，可设其形函数为：编程如下：clear allsyms b1 b2 b3 v VV F L x y EI %VVÊÇ×ܱäÐÎÄÜ£¬vÊÇÐκ¯Êýv=b1*x2VV=EI/2*int(diff(v,2,x)2,x,0,L)-F*subs(v,'x',L) %

25、注意势能的计算与均布载荷不同s1=diff(VV,b1)B1=solve(s1,b1)v1=subs(v,'b1',B1)y=F*x2/6/EI*(3*L-x)py=subs(y,'x',(0:0.05:1)*L)/F/L3*EIpx=(0:0.05:1)pv1=subs(v1,'x',(0:0.05:1)*L)/F/L3*EIplot(px,py,'.',px,pv1)legend('real solve','v1') clear b1 b2 b3syms b1 b2 b3 v VV F L x

26、y EI %VVÊÇ×ܱäÐÎÄÜ£¬vÊÇÐκ¯Êýv=b1*x2+b2*x3VV=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.b1B2=B.b2v2=subs(v,'b1',B1)v2=subs(v2,

27、9;b2',B2)px=(0:0.05:1)pv2=subs(v2,'x',(0:0.05:1)*L)/F/L3*EIplot(px,py,'.',px,pv1,px,pv2)legend('real solve','v1','v2')从图中可以看出，一阶近似时有误差，二阶近似时，已经收敛到精度解。作业：亲自验算上例或者自行定义形函数再计算。或者计算其他形式的梁模型。上机题目五 压杆临界压力计算实验目的：1 掌握最小势能的原理及计算过程。2 用MATLAB实现最小势能原理计算临界压力计算。一端固定一端自由：

28、基本原理（见变分方法讲义中的方法：）边界条件为：。压杆弯曲变形的变形能：杆上端由于弯曲变形而下降的高度为：所以外力的势能：总势能：按最小势能原理：由于两系数不能同时为零，所以要注上述方程组的系数矩阵的行列式等于零：，整理后：所以：。精确解为：，两者误差为：。命令：clear allsyms b1 b2 v VV F L x y EI %最高算到二阶，因为通过计算发现三阶出现虚根b=b1,b2i=2 %定义的阶数，i=1时是一阶，i=2是二阶v=0for ii=1:i v=v+b(ii)*x(ii+1) %定义形函数endVV=EI/2*int(diff(v,2,x)2,x,0,L)-F*int(1/2*(diff(v,x)2,x,0,L)for ii=1:i s(ii)=diff(VV,b(ii) %得到变分方程endfor ii=1:i for ij=1:i AAA(ii,ij)=diff(s(ii),b(ij)