弹性力学上机题目

弹性力学上机报告的要求：从以下六个题目中选取三个，完成每个题目后面的作业即可，上机报告的基本格式另见文件。弹性力学上机题目上机题目一主应力的数值计算实验目的：1进一步了解主应力、主应变、主方向。2掌握用MATLAB计算主应力、主应变、主方向的方法理论介绍：主应力在某点的应力状态中，如果在某个斜面上的切应力等于零，则该斜面上的正应力称为该点的一个主应力。主应变在某点的应变状态中，如果在某个斜面上的切应变等于零，则该斜面上的线应变称为该点的一个主应变。主方向对应于主应力或者主应变的斜面的法线方向称为应力主方向或者应变主方向。对于同一点，应力主方向与应变主方向是一致的。二向状态下的主应力：212222xyxyxy主方向：1121tan,tanxyxxyx三向状态下的主应力：主应力满足方向：321230式中：1222222232xyzxyyzzxxyyzzxxyzxyzyzxzxyxyyzzx其主方向需进一步讨论得到。通过研究发现，如果把应力分量写成应力矩阵：xyxzxxyyzyxzyzzs则主应力及其相对应的主方向就是应力矩阵的特征值及其对应的特征向量。MATLAB中提供了计算矩阵特征值及特征向量的专门函数：eig()，该函数使用格式如下：[V,D]=EIG(X)其中D和V分别为X矩阵的特征值及特征向量。学生需要再编程将特征向量换算成主方向的方向余弦。1平面应力状态计算：算例：50MPa20MPa10MPa10MPa20101050xyxxyys这里的应力正负号采用弹性力学的规定方法。具体命令如下：A=[-20,10;10,-50]//定义应力矩阵A=-201010-50[V,D]=eig(A)//调用EIG（）函数，V返回值是特征向量（按列，对应于特征值），也就是方向角与x,y正轴间夹角（方向角）的余弦（方向余弦）；D返回特征值（即主应力大小，按从小到大排序，这与材料力学的规定正好相反）V=-0.2898-0.95710.9571-0.2898D=-53.027800-16.9722jiao=rad2deg(atan(V(2,:)./V(1,:)))//计算出两个主方向与X轴正向的交角（逆时针为正）jiao=-73.155016.8450jiao2=rad2deg(acos(V))//计算出两个主方向的方向角jiao2=106.8450163.155016.8450106.8450以上计算结果告诉我们：该单元体的主应力为：maxmin16.977253.0278对应的主方向分别为：0116.84573.155若按方向角描述：两个主方向的方向角分别是：11163.155,106.84522106.845,16.845有几何知识不难得到主方向的物理位置。2三向应力状态计算：30MPaxyza（）70MPa40MPa50MPa70400403000050xyxzxxyyzyxzyzzs具体命令如下：A=[70,40,0;40,30,0;0,0,50]A=70400403000050[V,D]=eig(A)V=0.52570-0.8507-0.85070-0.525701.00000D=5.278600050.000000094.7214jiao2=rad2deg(acos(V))jiao2=58.282590.0000148.2825148.282590.0000121.717590.0000090.0000以上计算结果表明：该单元体的主应力为：12394.7214,50,5.2786三个主应力对应的三个主方向的方向角分别是：111148.2825,121.7175,9022390,90,033358.2825,148.2825,90作业题：任选两个进行计算。160MPa80MPa30MPa30MPaxyz40MPa50MPa60MPa30MPaxyz40MPa50MPa70MPa30MPa上机题目二弹性力学公式的辅助推算实验目的：1初步了解MATLAB强大的符号推演功能。2用MATLAB的符号计算功能，导出弹性力学中的方程。例1：导出用应变表示应力的物理方程：221121xxyyyxxyxyEEE命令如下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*sxyls1=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极坐标解答中，若设定应力函数为()cos2f，试用MATLAB导出各应力分量。已知应力分量表示式：22222111，相容方程：432443223()2()9()9()cos(2)0dfdfdfdfdddd。命令如下：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^4sfsf=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页的结果进行对比，并继续算下去。作业题，请导出一至两个推导过程，例如：边界条件中待定系数的确定；三维问题中有关问题等。也可亲自验算以上两例。上机题目三应力变分法在平面问题中应用实验目的：1掌握平面问题的应力变分方法的原理及计算过程。2用MATLAB实现平面问题的应力变分方法的计算过程，并讨论形函数的收敛性。平面问题的应力变分方法的基本原理：平面问题的应变余能可以表示为：222122cxyxyVdxdyE引入应力函数后：22222222122cxyVfxfydxdyEyxxy应力变分方程是：cxyzVufvfwfdS由于面力是给定的，所以上式右端等于零。即0cV。如果应力函数的构成是通过形函数的系数mA构成，则余能的变分方程就写成如下的偏微分形式：0cmVA例题（P284，11-2）：正方形薄板，边长为2a，如图所示在左右两边受有按抛物线分布的拉力，即2xxayqa。试用应力变分法按如下的应力函数求解：22422222123222221112qyxyxyqaAAAaaaaa。aaaaOxyqqqq%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