弹性力学上机题目教学文案

弹性力学上机题目精品文档弹性力学上机报告的要求：从以下六个题目中选取三个，完成每个题目后面的作业即可，上机报告的基本格式另见文件。弹性力学上机题目上机题目一 主应力的数值计算实验目的：进一步了解主应力、主应变、主方向。掌握用MATLAB计算主应力、主应变、主方向的方法理论介绍：主应力在某点的应力状态中，如果在某个斜面上的切应力等于零，则该斜面上的正应力称为该点的一个主应力。主应变在某点的应变状态中，如果在某个斜面上的切应变等于零，则该斜面上的线应变称为该点的一个主应变。主方向对应于主应力或者主应变的斜面的法线方向称为应力主方向或者应变主方向。对于同一点，应力主方向与应变主方向是一致的。二向状态下的主应力：21xyxy222xy2主方向：tan11x,tan2xyxy1x收集于网络，如有侵权请联系管理员删除精品文档三向状态下的主应力：主应力满足方向：320123式中：1xyz2222xyyzzxxyyzzx2222xyyzzx3xyzxyzyzxzxy其主方向需进一步讨论得到。通过研究发现，如果把应力分量写成应力矩阵：s

x yx zxxy y zyxz yz z则主应力及其相对应的主方向就是应力矩阵的特征值及其对应的特征向量。MATLAB中提供了计算矩阵特征值及特征向量的专门函数： eig()，该函数使用格式如下：[V,D]=EIG(X)其中D和V分别为X矩阵的特征值及特征向量。学生需要再编程将特征向量换算成主方向的方向余弦。平面应力状态计算：算例：x yxsxy y

10MPaMPa10MPaMPa2010这里的应力正负号采用弹性力学的规定方法。10 50具体命令如下：收集于网络，如有侵权请联系管理员删除精品文档>>A=[-20,10;10,-50] //定义应力矩阵A=-20 1010 -50>>[V,D]=eig(A) //调用EIG（）函数，V返回值是特征向量（按列，对应于特征值），也就是方向角与 x,y正轴间夹角（方向角）的余弦（方向余弦）； D返回特征值（即主应力大小，按从小到大排序，这与材料力学的规定正好相反）V=-0.2898 -0.95710.9571 -0.2898D=-53.0278 0-16.9722jiao=rad2deg(atan(V(2,:)./V(1,:)))//计算出两个主方向与X轴正向的交角（逆时针为正）jiao=-73.1550 16.8450>>jiao2=rad2deg(acos(V)) //计算出两个主方向的方向角jiao2=106.8450163.155016.8450106.8450收集于网络，如有侵权请联系管理员删除精品文档max以上计算结果告诉我们：该单元体的主应力为：min对应的主方向分别为：016.845o73.155o1若按方向角描述：两个主方向的方向角分别是：1 163.155o, 1 106.845o106.845o,216.845o有几何知识不难得到主方向的物理位置。三向应力状态计算：y30MPa40MPa70MPa50MPaxz（a）xyxzx70400sxyyzy40300xzyzz0050具体命令如下：A=[70,40,0;40,30,0;0,0,50]A=70 40 040 30 00 0 50

16.977253.0278收集于网络，如有侵权请联系管理员删除精品文档[V,D]=eig(A)V=0.52570-0.8507-0.85070-0.525701.00000D=5.278600050.00000094.7214jiao2=rad2deg(acos(V))jiao2=58.282590.0000148.2825148.282590.0000121.717590.0000090.0000以上计算结果表明：该单元体的主应力为： 1 94.7214, 2 50, 3 5.2786三个主应力对应的三个主方向的方向角分别是：1148.2825o,1121.7175o,190o290o,290o,30o58.2825o,3148.2825o,390o作业题：任选两个进行计算。收集于网络，如有侵权请联系管理员删除精品文档yy30MPa30MPa50MPa40MPa80MPa60MPa70MPa40MPa30MPax30MPax50MPa160MPazz收集于网络，如有侵权请联系管理员删除精品文档上机题目二 弹性力学公式的辅助推算实验目的：1初步了解MATLAB强大的符号推演功能。2用MATLAB的符号计算功能，导出弹性力学中的方程。例1：导出用应变表示应力的物理方程：Ex12xyyEyx12Exy xy21命令如下>>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.sysy=-(EE*ey+EE*ex*psu)/(psu^2-1)sxy=ls1.sxysxy=(EE*exy)/(2*psu+2)例2极坐标解答中，若设定应力函数为 f( )cos2 ，试用MATLAB 导出各应力分量。112222已知应力分量表示式：2，1收集于网络，如有侵权请联系管理员删除精品文档相容方程：4d4f()2d3f()9d2f()9df()。cos(2)4d32d230dd命令如下：>>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 //得到f 的表达式。>>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)srfsrf=-(2*sin(2*FI)*(C2*t^2-C4*t^4-3*C3*t^6+3*C5))/t^4可以将些结果与课本 P75页的结果进行对比，并继续算下去。作业题，请导出一至两个推导过程，例如：边界条件中待定系数的确定；三维问题中有关问题等。也可亲自验算以上两例。收集于网络，如有侵权请联系管理员删除精品文档上机题目三 应力变分法在平面问题中应用实验目的：掌握平面问题的应力变分方法的原理及计算过程。用MATLAB实现平面问题的应力变分方法的计算过程，并讨论形函数的收敛性。平面问题的应力变分方法的基本原理：平面问题的应变余能可以表示为： Vc

1

2222dxdyxyxy2E1222222引入应力函数后：Vcfxxfyy2dxdy2Ey2x2xy应力变分方程是：VcufxvfywfzdS由于面力是给定的，所以上式右端等于零。即 Vc 0。如果应力函数的构成是通过形函数的系数 Am构成，则余能的变分方程就写成如下的偏微分形式：Vc

0Am例题（P284，11-2）：正方形薄板，边长为 2a，如图所示在左右两边受有按抛物线分布的拉力，即求解：

2qy。试用应力变分法按如下的应力函数aaqy4x22y222y22x12a2qa1a21a2A1A2a2A3a2L。收集于网络，如有侵权请联系管理员删除精品文档q a a qaO xqaqy%clearall %清理已有变量cleara1a2a3a4a5symsaxyFIsxsysxyqa1a2a3a4 %定义符号变量，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^2VC=int(s4,'x','-a','a') %以下两式是进行余能的积分运算VC=int(VC,'y','-a','a')s1=diff(VC,a1) %对余能取变分A1=solve(s1,a1) %令余能变分等于零，计算待定系数 A1sx=diff(FI,2,y) %计算应力分量sx=subs(sx,'a1',A1) %将计算出的待定系数回代入应力分量表达式sx=subs(sx,'x',0) %令x=0sx=vpa(sx,2) %将应力表达式数值化并取 2个有效数字sx1=simple(sx) %简化并重新命名应力变量，为下一步计算做准备以下计算3阶收集于网络，如有侵权请联系管理员删除精品文档cleara1a2a3a4a5%clearallsymsaxyFIsxsysxyqa1a2a3a4FI=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^2VC=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=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%clearallsymsaxyFIsxsysxyqa1a2a3a4a5FI=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^2VC=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.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*y^2)/a^2sx3=(q*(0.137*a^4+0.804*a^2*y^2-0.357*y^4))/a^4sx5=(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,'-')收集于网络，如有侵权请联系管理员删除精品文档0.80.70.60.50.40.30.20.1-0.8-0.6-0.4-0.200.20.40.60.81-1结果发现，5阶的结果与3阶几乎没有多少改进，说明本题在 3阶时基本上已经收敛到精确解附近。作业题，11-1，11-3，或者亲自验算上例。收集于网络，如有侵权请联系管理员删除精品文档上机题目四 位移变分法在梁模型中的应用实验目的：掌握位移变分方法的原理及计算过程。用MATLAB实现梁的位移变分方法的计算过程，并讨论形函数的收敛性。位移变分方法的基本原理：将位移函数可以设为： u u0 Amum，v v0 Bmvm，w w0 Cmwm式中：u0,v0,w0满足给定的位移边界条件； um,vm,wm只需要满足边界为零的条件，并且称其为形函数。目前的情况下，位移函数称为 RIZE法形函数。若形函数在满足上述条件的同时，还能满足给定的应力边界条件，则又可称为伽辽金位移形函数。将形变势能V用位移函数表示后，位移变分公式就可以写成：VfxumdSfxumdxdydzAmVfyvmdSfyvmdxdydzBmVfzwmdSfzwmdxdydzCm解以上的方程组即可确定出待定系数。如果研究对象是一个梁，则方程又可以得到简化：首先位移可设vv0Bmvm，其变形能：2lM2EId2vEId2v2Vldx2ldxdxdx2dx202EI02EI0收集于网络，如有侵权请联系管理员删除精品文档RIZE（李滋）方程可写成为：VfyvmdlfyvmdxdyBm若忽略体力则：V（a）fyvmdlBm此方程组便可确定待定系数 Bm。其实，也可以采用最小热能原理的方法：VVfyvdl。式中V是系统总势能，V是系统应变能，fyvdl是外力做功的负值，就是外力的势能。根据最小势能原理，有：V(b)0Bm(a),(b)方程其实是等价的，读者可以证明。例1选用合适的形函数形式计算图示简支梁的挠曲线，并与精确解（yqxl32lx2x3）绘图对比。24EIqxyyxy本例中，边界条件是：x0,xl处，v0，所以其形函数可以设为：v0B1x(lx)B2x2(lx)2B3x3(lx)3...收集于网络，如有侵权请联系管理员删除精品文档不难看出v0 0。命令如下：clearallsymsb1b2b3vVVqLxyEI %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*EIpx=(0:0.05:1)pv1=subs(v1,'x',(0:0.05:1)*L)/q/L^4*EIplot(px,py,'.',px,pv1)legend('realsolve','v1')以下是二阶形函数计算：clearb1b2b3symsb1b2b3vVVqLxyEI %VV是总变形能，v是形函数v=b1*x*(L-x)+b2*x^2*(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.05:1)*L)/q/L^4*EIplot(px,py,'.',px,pv1,px,pv2)legend('realsolve','v1','v2')收集于网络，如有侵权请联系管理员删除精品文档0.0140.0120.010.0080.0060.0040.002

realsolvev1v200 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1从图中可以看出，二阶形函数已经与精确解重合了。例2仍是上题，如果选用vB1sinxB2sin2xB3sin3x...形函数，LLL试采用位移变分法，或者最小势能原理（两者公式上完全一致）计算待定系数，并与精确的挠曲线进行对比。只需要把例1中形函数的表达式修改一下，就可以了。三角级数收敛非常快，本题中只需要一项形函数，就可以达到理想精度。并且，计算后：B1=(4*L^4*q)/(EI*pi^5)B2=0例3：悬臂梁受集中力作用下，材料力学中，其挠曲线是2y Fx (3L x)。试采用幂级数形式的形函数计算其挠曲线，并与精确解对6EI比。收集于网络，如有侵权请联系管理员删除精品文档Fxy边界条件是，在x=0处，挠度和转角都等于零，所以，可设其形函数为：v B1x2 B2x3 B3x4 ...编程如下：clearallsymsb1b2b3vVVFLxyEI%VVê?×ü±?D??ü￡?vê?D?oˉêyv=b1*x^2VV=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*EIpx=(0:0.05:1)pv1=subs(v1,'x',(0:0.05:1)*L)/F/L^3*EIplot(px,py,'.',px,pv1)legend('realsolve','v1')clearb1b2b3symsb1b2b3vVVFLxyEI%VVê?×ü±?D??ü￡?vê?D?oˉêyv=b1*x^2+b2*x^3VV=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,'b2',B2)px=(0:0.05:1)收集于网络，如有侵权请联系管理员删除精品文档pv2=subs(v2,'x',(0:0.05:1)*L)/F/L^3*EIplot(px,py,'.',px,pv1,px,pv2)legend('realsolve','v1','v2')0.350.30.250.20.150.10.05

realsolvev1v200 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1从图中可以看出，一阶近似时有误差，二阶近似时，已经收敛到精度解。作业：亲自验算上例或者自行定义形函数再计算。或者计算其他形式的梁模型。收集于网络，如有侵权请联系管理员删除精品文档上机题目五 压杆临界压力计算实验目的：掌握最小势能的原理及计算过程。用MATLAB实现最小势能原理计算临界压力计算。一端固定一端自由：PLx基本原理（见变分方法讲义中的方法：）23ya1xa2xLL边界条件为：y(0)0,y(0)0,y(l)0。压杆弯曲变形的变形能：U1L2dx2EI23a1a221EIy3a13a220L杆上端由于弯曲变形而下降的高度为：dy2y21L2dxdx(1cos)dx[1cos(y)]dx[1(1)]2dx2y20所以外力的势能：收集于网络，如有侵权请联系管理员删除精品文档U2P1PL2dxP20a1245a1a227a22y2030L总势能：UU1U2EIa23aa3a2P20a245aa227a22L3112230L112按最小势能原理：U0:4EI4Pa16EI3Pa20a1L33LL32LU0:6EI3Pa112EI9Pa20a2L32LL35L由于两系数a1,a2不能同时为零，所以要注上述方程组的系数矩阵的行列式等于零：4EI4P6EI3PL33LL32L0，整理后：,3P12EI9P6EIL32LL35L3P2L4104PL2EI240E2E20所以：P2.49EI2EIEI，两者误差为：2。精确解为：Pcr20.24762L4LL0.92%。命令：clearallsymsb1b2vVVFLxyEI %最高算到二阶，因为通过计算发现三阶出现虚根b=[b1,b2]i=2 %定义的阶数， i=1时是一阶， i=2是二阶v=0forii=1:iv=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)forii=1:is(ii)=diff(VV,b(ii)) % 得到变分方程end收集于网络，如有侵权请联系管理员删除精品文档forii=1:iforij=1:iAAA(ii,ij)=diff(s(ii),b(ij)) % 得到变分的系数矩阵endenddA=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。作业：尝试采用三角函数做为形函数计算此例，或者另外其他类型的压杆如：两端铰支。上机题目六 对径对压薄圆板应力分析实验目的：收集于网络，如有侵权请联系管理员删除精品文档利用弹性力学解，分析对径受压薄圆板内应力分布情况。用MATLAB实现两条直径上应力